一種基于圖割的機載LiDAR單木識別方法*

王 濮，邢艷秋，王 成，習曉環(huán)

(1 東北林業(yè)大學 森林作業(yè)與環(huán)境研究中心，哈爾濱 150040； 2 中國科學院遙感與數(shù)字地球研究所 數(shù)字地球重點實驗室， 北京 100094)

機載激光雷達(light detection and ranging, LiDAR)可以直接、高效地獲取高精度、高密度的森林植被三維結(jié)構(gòu)信息，而且其光斑直徑較小、全數(shù)字化記錄回波信息(回波次數(shù)、回波強度)[1]，適于單木位置的精準探測與樹冠的精確分割，可為森林三維結(jié)構(gòu)定量、精細描述[2](如樹高[3]、冠幅、樹種等)提供數(shù)據(jù)源保障，因而也成為當前森林資源調(diào)查最有效的手段之一[4]。

國內(nèi)外學者利用機載LiDAR在單木結(jié)構(gòu)參數(shù)，特別是垂直結(jié)構(gòu)參數(shù)方面提出了多種單木自動探測方法，其中基于冠層高度模型(canopy height model, CHM)進行單木分割是主流方法之一，主要包括局部最大值法、區(qū)域增長算法、分水嶺算法等。通常情況下可以認為從樹頂反射的激光點是單木冠層的最高點，即樹頂通常表現(xiàn)為局部最高點。因此，局部最大值算法通常以固定或可變窗口尺寸來檢測單木樹頂。李響等[5]采用動態(tài)窗口局部最大值法對郁閉度較高的針葉林進行單木位置自動提取，精度可達85%。區(qū)域增長算法從種子像素開始，以迭代的方式從鄰域不斷生長，直至滿足生長停止條件。甄貞等[6]采用標記控制區(qū)域增長方法勾繪樹冠邊界，但針葉林易出現(xiàn)過度生長情況。分水嶺算法是區(qū)域增長算法的一種特例，將冠層模型倒置，并將其想象成一個盆地，以局部最小值為入口開始注水，水位上升直至盆地的邊界。Chen等[7]提出一種基于冠層極值模型(canopy maxima model, CMM)的樹頂探測方法，并基于此應用分水嶺算法分割冠層。這些方法存在3點不足之處：1)CHM更多地描述冠層上部，缺少對冠層下部及垂直剖面信息的描述；2)插值方法和格網(wǎng)間距的選擇易引入空間誤差和不確定性；3)CHM表面平滑程度難以自動控制，可能會引起樹高的高估或低估[8-9]。

近年來機載LiDAR系統(tǒng)飛速發(fā)展，獲取的點云數(shù)據(jù)密度和精度都極大提高，更多研究趨向于直接利用點云數(shù)據(jù)進行單木探測[8]，不僅提高了冠層上部的單木探測準確率，而且提升了次生冠層單木探測的可能性。Li等[10]提出基于空間距離的單木分割方法，避免CHM插值誤差和不確定性，并在針葉混交林得到較好的應用效果?？紤]到樹冠的各向異性，Zhang等[9]提出基于滑動膨脹環(huán)線圈方法(the donut expanding and sliding algorithm)分割單木，但效率較低且難以探測到下層單木。此外，計算機視覺和圖像處理領域的一些經(jīng)典算法(K-means聚類算法、體素化分割算法等)也被用于激光點云數(shù)據(jù)的單木分割。Kandare等[11]利用二維和三維K-means聚類算法從點云數(shù)據(jù)中檢測冠層上層、下層樹木，但閾值不易控制，易造成過分割。Wang等[12]用體素結(jié)構(gòu)和一種分層的形態(tài)學方法，先在同一高度層生成冠層區(qū)域，最后從不同高度層的樹冠區(qū)域中提取單木，但該方法易受激光點密度影響。

前人的方法雖然取得了一定的進展，但在應用于密集林區(qū)的重疊冠層分割、中下層單木探測等時效果甚微，而茂密植被區(qū)域的單木測量更是林業(yè)野外調(diào)查的難點。本研究試圖從森林冠層結(jié)構(gòu)分層的特點出發(fā)，結(jié)合歸一化割(normalized cut, Ncut)方法提取單木，以提高密集林區(qū)的單木探測精度。

1 方法

與冠層表面單木不同，密集林區(qū)的中下層單木受到抑制，特征不凸顯，無法直接對其進行探測?？紤]到森林冠層垂直結(jié)構(gòu)這種分層的特點，本研究首先利用局部最大值方法和Ncut方法對冠層激光點云數(shù)據(jù)探測并分割冠層上部的凸顯單木；并在此基礎之上，以樹冠激光點數(shù)閾值和樹冠形狀指數(shù)作為約束條件，全局最大值替代局部最大值，再次利用Ncut算法進行漏檢單木探測，識別森林中下層單木；由于缺乏單木位置的地面實測數(shù)據(jù)，單木探測的精度評價選用目視判讀所勾繪的單木位置作為參考，具體處理流程如圖1所示。

圖1 技術(shù)路線流程框圖Fig.1 Flow chart of single-tree segmentation

1.1 點云數(shù)據(jù)預處理

在LiDAR數(shù)據(jù)采集過程中，由于各種測量誤差和隨機誤差，原始點云數(shù)據(jù)中存在噪聲點，嚴重影響數(shù)據(jù)處理和單木分割效果，因而對點云數(shù)據(jù)處理前需要去除噪聲點。本文采用噪聲類型適應性較強的離群點去除方法去除點云噪聲[13]。然后，利用具有較強地形適應能力的布料濾波算法[14]進行地面點和非地面點分類，為高程歸一化等后續(xù)處理做準備工作。最后，為降低地形對單木樹高提取精度的影響，采用Lee等[15]方法對植被點云進行歸一化處理，每個點的歸一化高程值表示該點到地面的高度，因此位于樹頂?shù)募す恻c高程值可認為是此單木樹高。

1.2 樹頂探測

采用局部最大值方法探測冠層表面明顯的樹頂。首先分別對地面點和非地面點采用反距離加權(quán)法(公式(1)和(2))插值生成數(shù)字高程模型(digital elevation model, DEM)和數(shù)字表面模型(digital surface model, DSM)；然后，基于DSM和DEM生成CHM，并利用文獻[16]中方法填充CHM冠層中高程值突變的凹坑；最后，采用窗口為5×5的高斯濾波平滑CHM以減少樹頂誤判，并采用局部最大值濾波器進行樹頂識別。

(1)

(2)

式中：(X,Y,Z)為插值點坐標，(xi,yi,zi) 為搜索樣本點坐標，p為權(quán)重，q為冪次，i為樣本點坐標編號，n為搜索區(qū)域內(nèi)樣本點總數(shù)。

在冠層形狀規(guī)則和中低密集度的林區(qū)該方法具有良好的適用性，但應用于樹冠形狀復雜和高密集度的林區(qū)易發(fā)生誤判。本研究通過設定樹頂間的水平距離閾值來降低誤判率，即在排除灌木影響的情況下如果兩樹頂間的水平距離小于設定閾值，則對高程較低的樹頂進行剔除。

1.3 Ncut方法分割冠層

Ncut算法即基于圖論的圖像分割方法——歸一化割，通過由節(jié)點間相似性構(gòu)建的權(quán)值矩陣的特征向量分割圖像/點集[17]。Reiberger等[18]和Yao等[19]以樹頂、樹干的位置和歸一化割方法相結(jié)合用于激光點云數(shù)據(jù)單木分割，雖然冠層中下層單木探測率明顯提高，但樹干識別精度易受激光點云密度的影響。本文直接應用于激光點云，先計算局部最大值以確定冠層表面的明顯樹頂，以此作為單木位置的先驗知識并采用Ncut方法進行冠層初始分割，然后以全局最大值代替局部最大值作為Ncut方法的先驗，并作規(guī)則約束，完成對漏檢單木(主要冠層中下層單木)的進一步探測，進而實現(xiàn)森林單木的精確探測。算法描述如下：

將三維點云映射到無向圖G=(V,E)中，V為表示激光點的節(jié)點，E為連接每一對節(jié)點的邊。兩節(jié)點{i,j}∈V之間的相似性通過權(quán)重Wij來描述，即權(quán)重Wij表示為兩節(jié)點{i,j}∈V對應的激光點之間特征相似程度，有

(3)

Ncut分割通過切除圖G=(V,E)中的一些邊將圖G分割成為兩個非連接性點集A、B，且滿足A∩B=? 和A∪B=V，即類內(nèi)相似度最大、類間相似度最小。歸一化割準則如下

(4)

求解NCut(A,B)最優(yōu)解是一個NP完全問題(non-deterministic polynomial complete problem)，本研究通過對矩陣的特征值與特征向量的求解實現(xiàn)近似求解，這也是目前比較常用的方式。

(D-W)y=λDy,

(5)

式中:令n=|V|，權(quán)重矩陣W為n×n對稱矩陣，W(i,j)=Wij，{i,j}∈V；D為n×n對角矩陣，且D(i,i)=∑jW(i,j)。由Rayleigh 商的基本性質(zhì)可知，求解NCut(A,B)的最小值問題轉(zhuǎn)化為求解該特征系統(tǒng)的第二最小特征向量，即可利用該特征向量對圖進行劃分。但由于第二最小特征向量的元素一般呈現(xiàn)為連續(xù)實數(shù)值，因此需要引入一個分離點來二分，通常用0或特征向量元素的中值作為分離點。為使NCut(A,B)最小，即圖G得到最優(yōu)化分，采用試探法尋找最優(yōu)分離點[18]。

1.4 漏檢單木的探測

本文樹頂探測算法可能出現(xiàn)漏檢情況(如圖2)，即根據(jù)局部最大值算法只能檢測到冠層頂部明顯的樹頂。為降低失誤率，提高本文方法的普適性，將全局最大值代替局部最大值作為Ncut方法的先驗知識，并進行規(guī)則約束，對未分類點云中可能存在的單木或樹梢進行迭代檢測，方法描述如下：

輸入：Pts— 未分類的林區(qū)冠層點云數(shù)據(jù)

輸出：TreePts— 單木冠層激光點云

1) 定義Pts激光點集為U，做其高程值降序處理；

2) 選取最高點作為種子點，定義為pi(pi∈U)，統(tǒng)計半徑為dXY的領域點數(shù)N，若N小于閾值TD，則將其歸屬到其最近單木中；

3) 利用1.3中方法進行冠層分割得到種子點的單木冠層點云Ti，此時，U=U-Ti；

4) 利用文獻[20]的方法計算Ti的形狀指數(shù)(shape index, SI)，若判斷Ti為樹梢，則其屬于其周圍單木的一部分，否則，TreePts=TreePts+Ti；

5) 重復2)、3)、4)，直至Ui為空。

圖2 多株單木相鄰較近Fig.2 Example of multiple adjacent trees

2 結(jié)果與分析

利用美國華盛頓州西部的Capitol State森林藍嶺地區(qū)(The Blue Ridge area)的機載點云數(shù)據(jù)驗證本文方法的有效性。研究區(qū)以丘陵為主，海拔260～425 m，地面坡度0°～45°；主要樹種以針葉林為主，如道格拉斯冷杉、西部鐵杉和西部紅雪松等。機載點云數(shù)據(jù)由直升機搭載TopoEyeAB的TopEye系統(tǒng)獲取，時間為1999年春季，平均點密度約為4.86 pts/m2。本文選取3塊100 m×100 m的樣地作為研究樣地，標記為Plot 1、Plot 2、Plot 3(圖3)。

圖3 研究區(qū)高分辨率影像Fig.3 High-resolution image of the study area

圖4 原始點云與濾波結(jié)果Fig.4 Original point cloud and filter result

圖5 點云生成的DSM(a)和DEM(b)Fig.5 DSM (a) and DEM (b) generated using LiDAR data

然后，采用1.2所描述方法對研究區(qū)域濾波后的點云進行樹頂初步探測，提取出位于冠層上部的明顯樹頂。由于植被葉片間的空隙及遮擋性等其他因素影響，CHM影像上易出現(xiàn)凹坑。凹坑的存在，使得CHM不能正確表達冠層表面的結(jié)構(gòu)形態(tài)，同時可能導致樹高的低估，因此很有必要去除CHM局部區(qū)域凹坑。本文利用文獻[16]中方法，目視確定最佳閾值，進行凹坑去除。用窗口為5×5的局部最大值方法進行探測樹頂，為排除可能存在的干擾值，降低樹頂?shù)腻e誤率，引入樹頂間的水平歐式距離閾值；該閾值若設置過小，則樹頂錯誤率增加；若設置過大，則樹頂漏檢率增加。研究區(qū)域的實地調(diào)查，可作為該閾值調(diào)整的依據(jù)，本文該閾值設為4 m。圖6分別為3個樣方CHM影像和樹頂?shù)某醪教綔y結(jié)果。

將以上樹頂探測結(jié)果作為1.3中算法運行的先驗知識，進行迭代式分割處理，并將探測到的樹頂數(shù)作為迭代次數(shù)上限，每次分割出一棵單木。由于冠層的遮擋性和結(jié)構(gòu)分層的特點，利用1.2中方法只能探測到冠層表面主要單木的樹頂，以此為先驗知識的1.3中的冠層分割方法也只能分割出冠層表面主要單木的冠層，而位于冠層中下層的單木存在漏檢情況。對于漏檢單木，在未分類點集中采用1.4中方法以激光點高程全局最大值代替局部高程最大值，作為Ncut分割方法的先驗條件，從高到低依次進一步探測，最終冠層分割結(jié)果如圖7所示。

為評價本文算法精度，通過對激光點云和高

圖6 局部最大值法探測3個區(qū)域的樹頂結(jié)果Fig.6 Treetop detection results using the local maximum method

圖7 3個區(qū)域的單木冠層分割結(jié)果Fig.7 Segmentation results of three plots

空間分辨率航空影像目視判讀方式得到本研究單木位置精度的驗證數(shù)據(jù)，采用召回率r、正確率p、綜合考慮r和p的調(diào)和值F[10]來評價。召回率表示有效單木探測株數(shù)占真實參考數(shù)據(jù)的比例，正確率表示有效單木探測株數(shù)占整個提取結(jié)果的比率。單木分割算法將點云或影像分割成多個部分，每個部分可以看作一棵單木，分割的結(jié)果大致分為3種情況:正確分割(TP)、欠分割(FN)和過分割(FP)。因此，r，p，F(xiàn)可描述如下：

(6)

利用上述評價標準，本研究方法對漏檢單木處理前和處理后的分割結(jié)果分別進行評價，并與經(jīng)典的分水嶺方法比較，如表1所示。因本文測試數(shù)據(jù)Plot 3的株密度較小，不參與算法結(jié)果的定量比較，但仍然可以定性地說明本文方法在該株密度條件下具有較好的適用性。已知Plot 1的株密度大于Plot 2，對比可知，基于CHM的分水嶺算法在林區(qū)環(huán)境比較簡單的Plot 2研究區(qū)的結(jié)果與本文方法相當，但在株密度較高的Plot 1研究區(qū)域中的分割精度遠低于本文方法，有近27%的單木未探測到。繼續(xù)分析發(fā)現(xiàn)，本文研究方法在對漏檢單木處理前的分割結(jié)果與基于CHM的分水嶺算法相當，且在株密度較高的Plot 1研究區(qū)域中分割精度較差，原因是本文樹頂探測方法是基于CHM影像，缺少冠層垂直結(jié)構(gòu)信息，因此未能有效地探測到冠層中下層單木。根據(jù)表1，本文考慮到林區(qū)冠層結(jié)構(gòu)分層特點的單木探測方法，在茂密林區(qū)取得非常好的單木檢測率，并能有效地降低單木漏檢率。

3 結(jié)束語

本研究結(jié)合密集林區(qū)冠層結(jié)構(gòu)分層特點提出基于圖割的單木分割方法，實現(xiàn)冠層凸顯單木的分離，試驗結(jié)果表明本文方法在高株密度林區(qū)具有較高的單木探測精度。盡管如此，但當本文方法應用于高密度點云數(shù)據(jù)時，會增加Ncut的求解負擔，降低算法效率，自適應尺度的體素化高密度點云或許是一種解決辦法。另外，本研究只選擇了針葉林區(qū)域進行算法測試和分析，需要進一步嘗試算法對冠層結(jié)構(gòu)更加復雜的闊葉林或針闊混交林的適應性。這也是目前相關(guān)研究工作者所面臨的一個難題，利用分層分割的混合方法或新型手段(如：點云數(shù)據(jù)+高光譜影像/多光譜影像等)來解決此問題將是未來的研究方向。

表1 本文方法與分水嶺算法的對比Table 1 Comparison between the proposed method and the watershed algorithm