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

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

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

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

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

1 研究區(qū)概況

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

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

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

2 研究方法

2.1 研究路線

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

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

2.2 三維單木分割方法

2.2.1 局部最大值點探測

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

(1)

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

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

G=G5×5×C。

(2)

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

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

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

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

(3)

化簡式(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)

實驗采用最速下降法的算法計算X=(B，x0，y0，a)T的流程如下。

①給定初始值X0，ε>0，k=0計算

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

(8)

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

②計算

(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)。

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

(12)

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

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

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

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

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

2.3 精度評價方法

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

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

r=TP/(TP+FN)；

(13)

P=TP/(TP+FP)；

(14)

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

(15)

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

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

(16)

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

3 結(jié)果與分析

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

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

圖3 樹冠頂點識別結(jié)果

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

高斯模型聚類算法性能良好，相比于分水嶺算法單木分割召回率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ī)則來提取復(fù)雜闊葉林或針闊混交林的森林結(jié)構(gòu)參數(shù)，如生物量[21]、葉面積指數(shù)和蓄積量等。