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

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

(東北林業大學，哈爾濱，150040)

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

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

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

1 研究區概況

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

表1 實測單木結構參數的基本統計量

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

2 研究方法

2.1 研究路線

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

圖1 研究區彩色航空影像

2.2 三維單木分割方法

2.2.1 局部最大值點探測

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

(1)

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

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

G=G5×5×C。

(2)

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

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

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

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

(3)

化簡式(3)可得

(4)

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

(5)

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

(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)所在位置的灰度值大??；an為混合高斯曲面擬合后的樹冠范圍。

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

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

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

③重復步驟②，直到所有樣本點被劃分至相應簇時停止，實現單木三維分割。

2.3 精度評價方法

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

①采用召回率r、正確率p、調和值F3個指標衡量單木分割精度。

r=TP/(TP+FN)；

(13)

P=TP/(TP+FP)；

(14)

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

(15)

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

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

(16)

式中：n為正確分割的三維單木數量；Mi為分割單木結構參數；Ni為實測單木結構參數。

3 結果與分析

3.1 三維單木分割結果

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

圖3 樹冠頂點識別結果

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

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

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模型的局部最大值檢測進行單木分割，提高了11.3%準確率。

表2 高斯模型聚類法和分水嶺算法準確性評估

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

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

3.2 單木樹高和冠幅參數提取結果

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

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

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

4 結束語

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

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