應(yīng)用高斯聚類的單木分割及樹高和冠幅的提取1)

張怡卓 呂阿康 蔣大鵬 陳金浩 王克奇

(東北林業(yè)大學(xué)，哈爾濱，150040)

樹高和冠幅不僅是反映林分結(jié)構(gòu)特征的重要指標(biāo)，而且在生產(chǎn)潛力預(yù)估、生物量預(yù)測(cè)、林木生長(zhǎng)預(yù)測(cè)等研究中也是最常用的觀測(cè)變量。因此，單木分割及樹高冠幅提取在森林資源調(diào)查中具有重要的意義。機(jī)載激光雷達(dá)作為一種主動(dòng)遙感技術(shù)，具有較強(qiáng)的穿透能力，并不易受天氣及光照條件影響，可獲得高精度的地表及地物的垂直結(jié)構(gòu)信息。運(yùn)用機(jī)載激光雷達(dá)技術(shù)，可有效提取單木及樹高冠幅等參數(shù)信息[1-2]。

目前，單木分割與樹高冠幅提取方法可分為兩類：基于柵格化的冠層高度模型(CHM)和基于點(diǎn)云聚類的分割方法。其中，基于柵格化的CHM-分割方法是通過(guò)確定樹冠邊界或局部最大值識(shí)別樹冠頂點(diǎn)，采用區(qū)域生長(zhǎng)或圖像分割技術(shù)進(jìn)行單木分割與結(jié)構(gòu)參數(shù)提取。王軼夫等[3]在CHM上采用固定窗口識(shí)別樹頂點(diǎn)提取林木樹高，提取的樣地平均樹高與實(shí)測(cè)值線性回歸關(guān)系系數(shù)為0.694，固定窗口大小直接決定識(shí)別局部最大值的結(jié)果。李響等[4]、甄貞等[5]利用動(dòng)態(tài)窗口局部最大值法對(duì)單木位置進(jìn)行探測(cè)，采用標(biāo)記控制區(qū)域生長(zhǎng)法進(jìn)行冠幅勾繪的樹冠面積相對(duì)誤差平均值為8.74%，由于樹木生長(zhǎng)呈階段性變化，方法應(yīng)用具有一定局限性。耿林等[6]針對(duì)CHM表面存在的孔洞問題，采用線性平滑方法使樹冠邊緣連續(xù)，樹冠內(nèi)部得到填充，并利用局部動(dòng)態(tài)窗口探測(cè)樹冠頂點(diǎn)進(jìn)而提取樹高，但是冠層表面平滑程度難以控制，導(dǎo)致單木結(jié)構(gòu)參數(shù)提取精度不高?；邳c(diǎn)云聚類的分割方法，采用聚類算法直接對(duì)原始點(diǎn)云進(jìn)行單木分割與樹高冠幅提取，消除了柵格化CHM丟失數(shù)據(jù)的現(xiàn)象。Sandeep et al.[7]采用K-means聚類法實(shí)現(xiàn)對(duì)原始激光雷達(dá)點(diǎn)云數(shù)據(jù)的單木分割，但單木參數(shù)提取結(jié)果依賴于初始值的選取。Ferraz et al.[8]提出基于均值偏移算法的聚類方法實(shí)現(xiàn)單木分割，該算法不需要選擇聚類的數(shù)量，但因森林環(huán)境復(fù)雜，樹木大小不均勻，均值滑動(dòng)窗口大小難以確定。Li et al.[9]將全局最大值點(diǎn)作為最高樹的頂點(diǎn)，根據(jù)點(diǎn)云空間距離判別進(jìn)行聚類，在單木樹高提取方面達(dá)到較好效果，但未能準(zhǔn)確勾畫出樹木冠層部分。Lu et al.[10]根據(jù)點(diǎn)云的強(qiáng)度和回波特征分離樹干，將點(diǎn)識(shí)別為樹木點(diǎn)和噪聲點(diǎn)，但該算法需要更加密集的點(diǎn)云數(shù)據(jù)，數(shù)據(jù)采集成本高。趙晨陽(yáng)等[11]在確定聚類數(shù)量基礎(chǔ)上，以樹冠頂點(diǎn)為聚類中心，經(jīng)過(guò)迭代后估測(cè)單木冠幅的相關(guān)性系數(shù)為0.76，冠幅分割穩(wěn)定性較差。Ayrey et al.[12]提出基于原始點(diǎn)云的分層疊加聚類方法，該方法對(duì)整個(gè)森林點(diǎn)云進(jìn)行高度間隔切片，并對(duì)每一層進(jìn)行聚類實(shí)現(xiàn)樹冠形狀的三維重建，但方法依賴聚類條件。

基于柵格化的冠層高度模型(CHM)的單木分割及樹高冠幅提取方法易造成點(diǎn)云丟失，且樹冠部分不完整導(dǎo)致噪聲點(diǎn)出現(xiàn)；而利用原始三維點(diǎn)云進(jìn)行單木分割及樹高冠幅提取避免點(diǎn)云信息丟失，但分割結(jié)果需以單木數(shù)量為已知搜索條件，故該方法存在局限性。本研究以針葉材為研究對(duì)象，利用機(jī)載激光雷達(dá)獲取點(diǎn)云數(shù)據(jù)，提出一種基于高斯模型聚類的單木分割及樹高和冠幅的提取方法，根據(jù)單木形狀結(jié)構(gòu)建立高斯模型，通過(guò)確定單木位置及初步樹冠范圍，利用高斯輪廓約束聚類得到樹高和冠幅的精準(zhǔn)提取。

1 研究區(qū)概況

研究區(qū)位于美國(guó)華盛頓州國(guó)會(huì)森林藍(lán)嶺地區(qū)，該區(qū)域地貌主要以丘陵為主，海拔260～425m，地面坡度0°～45°。研究區(qū)域內(nèi)以針葉林為主，森林植被樹種主要為道格拉斯冷杉(Pseudotsugamenziesii)、西部鐵杉(Tsugaheterophylla)和西部紅雪松(Thujaplicata)。根據(jù)華盛頓州自然資源部和華盛頓大學(xué)對(duì)藍(lán)嶺地區(qū)347個(gè)野外調(diào)查點(diǎn)進(jìn)行地形調(diào)查實(shí)測(cè)單木結(jié)構(gòu)參數(shù)的基本統(tǒng)計(jì)量見表1。

表1 實(shí)測(cè)單木結(jié)構(gòu)參數(shù)的基本統(tǒng)計(jì)量

本研究數(shù)據(jù)采集時(shí)間為1999年春季，采用Saab TopEye LiDAR系統(tǒng)獲取小光斑機(jī)載雷達(dá)點(diǎn)云數(shù)據(jù)。遙感平臺(tái)飛行高度為200 m，飛行傾角為8°，采樣密度4次/m2脈沖回波，激光脈沖速度7 000點(diǎn)/s，最大回波記錄為4次，光斑直徑為40 cm，直升機(jī)飛行速度為25 m/s，掃描寬度70 m。實(shí)驗(yàn)選取6塊半徑30 m的圓形純林樣地(實(shí)測(cè)單木共893株)，樣地面積共1.69 hm2，平均點(diǎn)密度約為4.86點(diǎn)/m2，標(biāo)記為Plot1、Plot2、Plot3、Plot4、Plot5、Plot6(樣地位置如圖1所示)。該研究區(qū)彩色航空影像以藍(lán)嶺地區(qū)通過(guò)自相關(guān)技術(shù)衍生的偽樹冠表面模型利用軟考貝攝影測(cè)量系統(tǒng)完成正攝校正，最終圖像分辨率為0.3 m。

2 研究方法

2.1 研究路線

高斯模型聚類的單木分割及其參數(shù)提取方法技術(shù)路線如圖2所示。針對(duì)區(qū)域LiDAR點(diǎn)云數(shù)據(jù)，使用Kriging插值法生成冠層高度模型(CHM)，通過(guò)形態(tài)學(xué)開運(yùn)算和高斯濾波得到GCMM模型；利用局部最大值法和最速下降法的高斯曲面擬合得到單木位置及初步冠幅范圍；將歸一化點(diǎn)云數(shù)據(jù)采用最小二維歐式距離聚類得到三維單木及樹高；最后，利用實(shí)測(cè)數(shù)據(jù)分別對(duì)提取的冠幅和樹高進(jìn)行精度驗(yàn)證和分析。

圖1 研究區(qū)彩色航空影像

2.2 三維單木分割方法

2.2.1 局部最大值點(diǎn)探測(cè)

首先將原始點(diǎn)云分為非地面點(diǎn)和地面點(diǎn)[13]，分別進(jìn)行插值運(yùn)算[14-18]，得到數(shù)字地面模型(DSM)和數(shù)字高程模型(DEM)，兩者作差值運(yùn)算求得CHM[19]；由于原始的CHM中存在黑色或灰色的孔洞，進(jìn)而影響樹頂點(diǎn)的探測(cè)。對(duì)CHM模型表面進(jìn)行形態(tài)學(xué)開運(yùn)算，將CHM模型中的孔洞替換為鄰域內(nèi)點(diǎn)云高度最大值點(diǎn)，形成冠層最大模型(CMM)模型，即

(1)

圖2 高斯模型聚類的單木分割及其參數(shù)提取方法技術(shù)路線

CMM模型可增強(qiáng)樹冠邊緣信息使樹冠表面更平滑，運(yùn)用高斯濾波方法對(duì)CMM模型進(jìn)行線性平滑處理得到GCMM模型，即

G=G5×5×C。

(2)

式中：G為GCMM模型；C為CMM模型；G5×5為5×5(像元)的高斯矩陣。

固定窗口能準(zhǔn)確地識(shí)別樹頂點(diǎn)，通過(guò)實(shí)驗(yàn)分析選用大小5×5(像元)窗口可以找到GCMM模型局部最大值點(diǎn)，使用GCMM模型可消除偽局部最大值，減少樹梢識(shí)別中錯(cuò)誤檢測(cè)。

2.2.2 最速下降法的高斯曲面擬合

通過(guò)GCMM模型的柵格點(diǎn)進(jìn)行高斯曲面擬合，該模型灰度圖像矩陣為EM×N(1≤i≤M,1≤j≤N)，在三維空間中當(dāng)x=i，y=j時(shí)z=eij，z代表圖像的灰度值大小，且滿足條件eij>0，其高斯曲面模型表達(dá)式為

(3)

化簡(jiǎn)式(3)可得

(4)

令B=lnA，則式(4)可改寫為

(5)

則上述問題轉(zhuǎn)化為求參數(shù)B、x0、y0、a(參數(shù)B、x0、y0、a都是1×n的數(shù)組)，使得目標(biāo)函數(shù)為

(7)

實(shí)驗(yàn)采用最速下降法的算法計(jì)算X=(B，x0，y0，a)T的流程如下。

①給定初始值X0，ε>0，k=0計(jì)算

F0=F(X0),g0=g(X0)=F(X0)。

(8)

式中：g(x)=(g1、g2、g3、g4)T。

②計(jì)算

(9)

Xk+1=Xk-tkg(Xk)；

(10)

Fk+1=F(Xk+1),gk+1=g(Xk+1)。

(11)

式中：H(X)=2F(X)=(hij)(4×4)。

運(yùn)用最速下降法求解出擬合后的混合高斯模型為

(12)

式中：(xn，yn)為混合高斯曲面擬合后的單木位置；An為(xn，yn)所在位置的灰度值大小；an為混合高斯曲面擬合后的樹冠范圍。

2.2.3 基于最小二維歐式距離聚類

基于混合高斯模型的最小二維歐式距離聚類是將高程歸一化點(diǎn)云進(jìn)行距離迭代判斷，劃分臨近點(diǎn)云歸屬情況，實(shí)現(xiàn)三維單木準(zhǔn)確分割。其利用最小二維歐式距離聚類算法如下：

①定義(xn，yn)所在空間位置為簇的中心線，提取樣本與中心線距離小于或等于4an的點(diǎn)云集群；

③重復(fù)步驟②，直到所有樣本點(diǎn)被劃分至相應(yīng)簇時(shí)停止，實(shí)現(xiàn)單木三維分割。

2.3 精度評(píng)價(jià)方法

本研究從二個(gè)方面對(duì)機(jī)載雷達(dá)點(diǎn)云提取單木及樹高冠幅能力進(jìn)行評(píng)價(jià)，通過(guò)實(shí)測(cè)地面數(shù)據(jù)和高空間分辨率圖片進(jìn)行目視解譯得到單木分割精度的驗(yàn)證，具體方法如下。

①采用召回率r、正確率p、調(diào)和值F3個(gè)指標(biāo)衡量單木分割精度。

r=TP/(TP+FN)；

(13)

P=TP/(TP+FP)；

(14)

F=2×(r×p)/(r+p)。

(15)

式中：TP為正確分割樹分割精度；FN為欠分割樹分割精度；FP為過(guò)分割樹分割精度。

②對(duì)正確分割的三維單木進(jìn)行樹高冠幅精度評(píng)價(jià)和線性回歸分析[20]，其精度評(píng)價(jià)計(jì)算公式為

(16)

式中：n為正確分割的三維單木數(shù)量；Mi為分割單木結(jié)構(gòu)參數(shù)；Ni為實(shí)測(cè)單木結(jié)構(gòu)參數(shù)。

3 結(jié)果與分析

3.1 三維單木分割結(jié)果

在6塊樣地中選取半徑為10 m的圓形樣地，圓形樣地共有實(shí)測(cè)樹木12株，基于兩種模型樹冠頂點(diǎn)識(shí)別結(jié)果如圖3所示?；诰植孔畲笾捣z測(cè)樹冠頂點(diǎn)時(shí)，采用GCMM模型可消除12個(gè)偽局部最大值，且樹冠頂點(diǎn)識(shí)別未受到影響。

圖3 樹冠頂點(diǎn)識(shí)別結(jié)果

采用不同大小窗口分別對(duì)6塊樣地的GCMM模型進(jìn)行局部最大值探測(cè)找到初步的樹冠頂點(diǎn)。6塊樣地實(shí)際單木為893株，相比實(shí)際單木值3×3(像元)窗口共檢測(cè)到1 153個(gè)樹冠頂點(diǎn)，多檢測(cè)到260個(gè)樹冠頂點(diǎn)，而7×7(像元)窗口共檢測(cè)到537個(gè)樹冠頂點(diǎn)，存在未識(shí)別樹冠頂點(diǎn)356個(gè)，因此實(shí)驗(yàn)選用5×5(像元)窗口找到GCMM模型中局部最大值點(diǎn)。

采用高斯模型聚類法和分水嶺方法分別對(duì)6塊樣地點(diǎn)云進(jìn)行單木分割，其分割結(jié)果如圖4所示。Plot5和Plot6樣地點(diǎn)云密度高，樹木數(shù)量較多，運(yùn)用分水嶺方法進(jìn)行單木分割存在嚴(yán)重的過(guò)分割現(xiàn)象，并且提取到樹冠閉合邊緣存在殘缺現(xiàn)象。由于針葉冠層俯視形似圓形，單木點(diǎn)云在GCMM模型中符合高斯分布，高斯模型聚類法初步提取到的樹冠范圍更加接近真實(shí)樹木形態(tài)。高斯模型聚類方法利用局部最大值法在GCMM模型上探測(cè)初步樹木位置，減少了單木過(guò)分割數(shù)量。通過(guò)與實(shí)測(cè)樹木位置(圖4中綠色圓點(diǎn)位置)進(jìn)行比對(duì)可知，該方法分割準(zhǔn)確性更高。

a1、b1、c1、d1、e1、f1為分水嶺算法；a2、b2、c2、d2、e2、f2為高斯模型聚類法。

如表2所示，6塊樣地共有893株，高斯模型聚類方法分割了860株，其中正確分割790株，過(guò)分割103株，欠分割70株?；诜炙畮X方法分割了867株，該算法正確分割為709株，過(guò)分割184株，欠分割158株。與分水嶺方法相比，高斯模型聚類法是基于GCMM模型的局部最大值檢測(cè)進(jìn)行單木分割，提高了11.3%準(zhǔn)確率。

表2 高斯模型聚類法和分水嶺算法準(zhǔn)確性評(píng)估

如表3所示，分水嶺算法平均調(diào)和值僅為0.80，是基于CHM模型參考空間位置和灰度值相近的像素點(diǎn)互相連接起來(lái)構(gòu)成一個(gè)封閉輪廓，形成單木冠層易導(dǎo)致過(guò)分割現(xiàn)象。而高斯模型聚類法是基于GCMM模型采用局部最大值法探測(cè)得到樹冠頂點(diǎn)，在樹冠頂點(diǎn)處使用最速下降法的高斯曲面擬合，最大程度擬合樹木真實(shí)形態(tài)，保證單木分割準(zhǔn)確性；但由于plot5和plot6林分郁閉度較高，激光信號(hào)無(wú)法識(shí)別林下較小樹木，平均F值為0.89。相比于分水嶺算法，高斯模型聚類法在6塊樣地分別得到r平均值為0.87，p平均值為0.91，證明該方法的優(yōu)越性。

表3 兩種方法分割精度得分

3.2 單木樹高和冠幅參數(shù)提取結(jié)果

使用高斯模型聚類方法和分水嶺方法對(duì)6塊樣地進(jìn)行單木樹高和冠幅提取，樹高提取精度分別為95%和90%，冠幅提取精度分別為91%和86%。兩種方法正確分割三維單木提取的冠幅與實(shí)測(cè)冠幅的線性回歸關(guān)系如圖5c、圖5d所示，高斯模型聚類方法R2=0.84，但分水嶺方法冠幅提取結(jié)果較弱，R2=0.73；單木樹高的線性回歸關(guān)系如圖5a、圖5b所示，分水嶺方法R2=0.84，相比于該方法高斯模型聚類方法樹高提取結(jié)果較好,R2=0.92。

從實(shí)驗(yàn)結(jié)果得出，高斯模型聚類法和分水嶺法樹高提取的平均誤差分別為-0.83、-1.41 m，平均誤差為負(fù)值，表明樹高存在低估現(xiàn)象，可能是由于激光點(diǎn)采樣密度不夠大，導(dǎo)致系統(tǒng)無(wú)法獲取全部樹木信息。利用兩種算法提取冠幅的平均誤差分別為-0.42、1.05 m，分析發(fā)現(xiàn)：本研究樣地樹木為針葉類型，樹冠俯視形狀近似圓形，高斯曲面擬合法比分水嶺算法提取到樹冠更接近真實(shí)情況，冠幅提取平均誤差減小0.63 m；相比于分水嶺方法，高斯模型聚類法直接對(duì)原始點(diǎn)云進(jìn)行聚類，避免了CHM模型采用柵格內(nèi)最大值法導(dǎo)致樹木信息丟失，保證點(diǎn)云信息完整性，提取樹高平均誤差減少0.58 m。

a.高斯模型聚類法樹高；b.分水嶺算法樹高；c.高斯模型聚類法冠幅；d.分水嶺算法冠幅。

4 結(jié)束語(yǔ)

高斯模型聚類算法性能良好，相比于分水嶺算法單木分割召回率r為0.87，正確率p為0.91，綜合r和p的調(diào)和值F為0.89，該算法提取冠幅與樹高精度分別為95%，91%，平均誤差分別為-0.42、-0.83 m，表明該算法具有良好的應(yīng)用前景。

現(xiàn)階段的單木分割及樹高冠幅的提取研究中，針葉純林仍占據(jù)重要地位，由于針葉林形狀大多為中心高，四周低的傘狀，分布較為規(guī)則，易在遙感影像中辨認(rèn)。在今后研究中基于高斯模型聚類方法可嘗試加入更多有效的約束規(guī)則來(lái)提取復(fù)雜闊葉林或針闊混交林的森林結(jié)構(gòu)參數(shù)，如生物量[21]、葉面積指數(shù)和蓄積量等。