基于機載LiDAR點云多層聚類的單木信息提取及其精度評價*

霍朗寧 張曉麗

(北京林業大學精準林業北京市重點實驗室 北京 100083)

激光雷達(light detection and ranging，LiDAR)是一種集激光、全球定位系統(GPS)和慣性導航系統(inertial navigation systems，INS)于一身的用于獲取地面及地面目標三維空間信息的主動式雷達探測技術系統。將LiDAR應用于林業領域，不僅能夠生成地面地形特征，還能夠得到林木樹冠表面和垂直結構信息，如樹高、冠幅、胸徑、生物量等，可為森林生物量、生長模型等研究提供數據支撐(N?sset，2002；Hyypp?etal.，2001)。近年來，隨著科學技術快速發展，LiDAR數據獲取成本不斷下降，應用于林業領域的研究也越來越多，特別是機載LiDAR，探測森林垂直結構的優勢日益突出，在林分和單木尺度森林信息提取方面的探索逐步深入(Pedersenetal.，2012；Hyypp?etal.，2012；Packalenetal.，2015)。在單木尺度方面，如果能夠通過機載LiDAR數據獲取林分中單木位置、樹高、冠幅等信息，將極大提高森林經營管理效率，為林木生長動態評估、林分擇伐作業、森林可視化經營等提供數據支持。

在LiDAR單木提取研究中，樹冠2D分割和3D空間信息分割是2種主要的分割方式(Hyypp?etal.，1999；Solbergetal.，2006；Tangetal.，2007；Ferrazetal.，2012)。其中，常用的2D分割方法有分水嶺分割(Tangetal.，2007；Eneetal.，2012)、區域增長法分割(Hyypp?etal.，1999；2001；Solbergetal.，2006)等，但這些方法存在一個共通的問題，即依賴CHM進行分割操作，CHM的濾波效果會極大影響最終提取效果(Eneetal.，2012)，且很可能導致信息丟失。從原理上分析，基于CHM構建并進行分割的方法仍為基于表面模型信息的分割，往往是在對最上層點云進行分割、劃分單木位置和冠幅范圍后將下層點云與之疊加，由于中下層點云利用率較低，數據信息沒有充分挖掘，復雜林分結構條件下幾乎無法提取到中下層林木，因此研究結果應用范圍受限。相比之下，3D空間信息分割則利用每一個點云的空間信息，單木提取更全面，信息也更完整(Lietal.，2012)。目前，常用的3D分割方法主要包括K均值聚類(Morsdorfetal.，2004；Guptaetal.，2010)、均值移動聚類(Ferrazetal.，2012)、自適應距離聚類(Leeetal.，2010；Guptaetal.，2010)等，然而這些方法并不一定都能獲得很好的提取效果，一些研究指出其中的單木提取方法誤差較大(Kaartinenetal.，2012；Eysnetal.，2015)，而對于提取單木正確率能夠達到多少的問題，很多研究避而不談(Guptaetal.，2010；Tangetal.，2013)。3D分割方法存在的困難主要有以下幾點：1)對高郁閉度林分下層林木提取能力較低，林分密度、郁閉度、林分結構異質性程度嚴重影響提取效果，提取效果較好的往往是林分密度較低的(Lietal.，2012)；2)提取林木有效性問題，由于3D分割方法對下層林木提取能力較低，導致能夠準確提取的林木數量占總體比例偏低(Kandareetal.，2016)，如果沒有地面實測數據進行匹配，將不能得知是哪部分林木被準確提取；3)點云密度問題，盡管有研究指出點云密度對單木樹高、冠幅精度影響不大(Reitbergeretal.，2009；Kaartinenetal.，2012；Yaoetal.，2014)，但在復雜林分結構條件下，點云密度過低將造成下層林木信息不完整，影響更新木的提取效果。

此外，以往研究在反演森林結構方面沒有過多關注，LiDAR反演森林結構的探索仍停留在林分尺度，且構造的指標僅適用于LiDAR數據，并沒有將其與實際生產實踐中常用的森林結構指標進行轉換，從而導致應用性較低。常用的森林結構指標如樹高基尼指數、樹高變異系數、聚集指數、競爭指數等，在林分空間結構優化、撫育改造等森林經營管理工作中具有重要指導意義，如果能夠通過LiDAR數據反演出高精度的指標，則可進一步將遙感調查結果向生產實踐應用轉化。

鑒于此，本研究聚焦上述問題，改進單木提取策略和算法，主要特色和創新點如下：1)挑選1塊林分密度高、空間結構異質性大的復層林作為研究對象，重點關注中下層林木的提取效果；2)在聚類前將點云劃分為多個水平層，增加下層林木探測能力，同時制定合理的點云簇融合方法，重點分析水平層設置情況對信息提取精度的影響；3)構建合理的精度評價指標，重點關注提取林木的正確性、有效性，對提取效果進行綜合評價，提升成果的可用性。

1 研究區概況與研究方法

1.1 研究區概況

研究區地處我國甘肅省祁連山國家級自然保護區肅南裕固族自治縣西水林場，為祁連山山區天然次生林、黑河流域水源涵養林，海拔2 700～3 000 m，屬溫帶高寒半干旱、半濕潤山地森林草原氣候，陽坡為山地草原，陰坡為森林景觀，林區主要森林類型為青海云杉(Piceacrassifolia)天然純林(成熟林)。研究區位置、樣地布置和原始點云示意如圖1所示。

圖1 研究區位置、樣地布置和原始點云示意

1.2 數據獲取與預處理

1.2.1 機載LiDAR數據與預處理 機載LiDAR數據獲取于2008年6月23日，無人機平均飛行海拔約3 560 m，距離地面飛行高度約760 m；搭載LiteMapper5600 LiDAR系統，激光掃描儀為RiegILMS-Q560，波長1 550 nm，激光脈沖長度3.5 ns，激光脈沖發散角小于等于0.5 m·rad，地面平均光斑直徑38 cm，可分辨目標最小間隔0.6 m，脈沖重復頻率50 kHz。

從原始點云中裁出與地面實測數據相匹配的100 m×100 m區域，采用Terrasolid軟件進行點云去噪、地面點與非地面點分離、植被標準化，并去掉高度小于2 m的點云，獲得樣地林層點云作為算法的試驗數據。平均點云密度為2.86 m-2。

1.2.2 地面調查數據 地面調查時間為2008年6月1—13日，沿山坡走勢布設一塊100 m×100 m超級樣地，方位角約122°，按照25 m×25 m尺寸劃分為16塊子樣地(圖1)。對樣地中林木進行每木檢尺，記錄林木胸徑、樹高、冠幅和枝下高，通過全站儀獲取單木位置，共得到1 435株單木信息。本研究重點討論LiDAR反演單木樹高、冠幅、位置以及林分結構參數的精度，相關統計量如表1所示。

1.3 單木分割方法

聚類是一種無監督學習，其將相似的對象歸到同一個簇中，簇內對象越相似，聚類效果越好。本研究應用K均值聚類算法，首先從數據集中隨機選取K個點作為初始聚類中心，然后計算各樣本到聚類中心的距離，將樣本歸到距離其最近的聚類中心所在的類。計算新形成的每一個聚類數據對象的平均值得到新的聚類中心，如果相鄰2次的聚類中心沒有任何變化，樣本調整結束，說明聚類準則函數已收斂。基于MATLAB中的聚類算法，完成單木分割、信息提取和精度評價。分層聚類方法提取單木工作流程如圖2所示。

圖2 分層聚類方法提取單木工作流程

1)點云水平切片 對于郁閉度高、株數密度大的復層林，直接對所有點云進行聚類分析存在2個問題：一是數據量大、運算速度慢，每加入一個新的聚類數據，則需要重新計算其聚類中心，將點云直接聚為2類的時間遠大于聚為4類后對相近2類進行融合的時間；二是初始聚類中心的選擇依賴于局部最大值位置，對于樹高異質性大的林分，只對整體提取局部最大值會不可避免錯過被壓木的樹冠頂點，從而無法分割優勢層下的低矮木和更新木。通過對點云進行水平切片(圖3)，可提取到各林層的局部最大值，尤其更利于下層林木頂點的偵測。切片片層高可根據林分分層情況人為設置，也可依照點云在垂直方向的分布規律自動分層，保證每層點云數量相同。

圖3 點云水平切片和局部最大值位置提取

2)點云局部最大值位置提取 決定聚類算法運行效率和聚類效果的主要因素之一是初始聚類中心選擇，點云數量龐大，初始聚類中心越接近最終聚類中心，運算效率越高。本研究區主要樹種為青海云杉，其形態近似圓錐體，樹頂點大致位于該樹垂直投影的形狀中心，因此通過提取點云局部最大值位置確定樹頂(圖3)，并將其作為初始聚類中心。利用MATLAB中的imregionalmax函數提取局部最大值點，其原理是通過一系列腐蝕、膨脹運算，從1個點的26個鄰近點判斷該點是否為局部最大值。由于點云數量龐大、異質性高，因此需對點云進行平滑濾波后再提取局部最大值。本研究采用高斯平滑濾波器對點云進行多次濾波，再提取點云局部最大值，所得結果為單木樹頂點坐標。

3)點云縱向壓縮K均值聚類算法適合于三維特征空間內大體為球體、簇密度均勻的數據集聚類分析，本研究區青海云杉林形態近似圓錐體，因此需對其進行縱向壓縮。采用3倍縮放比例，即點云x、y坐標不變，z坐標縮放為原來的1/3。

4)點云層聚類 計算每個縮放后的點云數據到每個聚類中心的歐氏距離，將點云歸到距離其最近的聚類中心所在的類，并重新計算聚類中心。重復該步驟至聚類中心不再變化，即完成點云聚類，可獲得標注不同類別屬性的點云簇。

5)隸屬不同點云層的點云簇融合 對接近分割線且聚類中心距離滿足一定條件的點云簇進行合并。合并條件具體設置和計算過程如下：經凸包處理求得點云簇面積，將點云簇視為圓求圓形半徑作為點云簇半徑；融合條件設置為上下2層點云簇中心的水平距離不超過下層點云簇半徑。合并順序如下：從最上層開始，先合并1、2兩個點云層的點云簇，得到新的點云簇和聚類中心，再與第3層點云簇合并，以此類推。合并完所有符合條件的點云簇后，重新計算點云簇中心，作為單木位置信息。

1.4 參數反演

1.4.1 單木參數反演 單木分類工作完成后，可獲得標注不同類別屬性的點云，具有同一屬性的點云為同一株單木的點云簇。對點云簇的x、y進行平均，作為單木位置信息；提取點云簇中z的最大值，作為單木樹高。在冠幅提取方面，由于點云密度較低，尤其在提取最下層被壓木時，存在因點云數量過少、點云共線而無法形成凸包的情況，因此不能應用凸包算法提取林木冠幅。本研究通過點云密度、數量計算冠幅面積的方式求得冠幅，即：

(1)

式中：c為冠幅；d為點云密度；n為點云簇中的點云個數。

1.4.2 林分結構參數反演 獲得所有單木的位置、樹高和冠幅后，進一步計算林分結構參數。

1)林分樹高分布 在林分中，不同樹高林木按樹高組的分配狀態，以2 m為一級，分別統計各樹高級范圍內的林木數量。為了量化LiDAR提取林分樹高分布與實際情況的差異性，計算各樹高級林木數量均方根誤差：

(2)

式中：n為樹高級數；N-test、N-ref分別為各樹高級中LiDAR提取數量和地面實測數量。

2)林分冠幅分布 在林分中，不同冠幅林木按樹高組的分配狀態，原理同樹高分布。

3)樹高基尼指數(Gc) 樹高基尼指數是目前成功應用于樹高不均等性測度的極少數指標之一(Gini，1912)，為樣地中所有個體間某一測定指標成對比較的差值絕對值的算術平均，范圍在0～1之間，用于表示樹高指標的不均勻性，即群落的樹高多樣性，指數值越接近1，則林分樹高越不均勻，指數值越接近0，則越均勻。計算公式為：

(3)

式中：Xi為第i株樹的樹高；i為研究區所有林木按高度從小到大排序后的順序；N為林木總數。

4)樹高變異系數(Cv) 樹高變異系數可反映林分垂直結構的變化程度，為樣地中所有樹高標準差和平均值的比較。計算公式為：

(4)

5)聚集指數(Clark-Evans’s aggregation indes)聚集指數是常用的林木空間分布格局指數(Clorketal., 1954)，為樣地種群內每個個體和其最近鄰體間的平均距離與其為隨機分布時的期望平均距離比較，所得比率為一隨機性偏離測度值，可以根據其統計顯著性判斷種群分布格局類型，取值范圍在0～2.149 6之間，若等于1為隨機分布，大于1為均勻分布，小于1為聚集分布。計算公式為：

(5)

式中：r0為樣地中每個個體和其最近鄰近體間的平均距離；ri為第i個個體和其最鄰近體間的距離；E(r)為樣地中所有個體隨機分布時r0的期望值；A為樣地面積(m2)。

1.5 精度評價

1.5.1 參數名稱設置 為便于結果呈現，本研究涉及的單木參數名稱及縮寫如表2所示。

表2 本研究涉及的單木參數名稱及縮寫

1.5.2 單木聚類效果評價指標 聚類經常會出現欠分類或過分類的情況，在聚類效果評價中，數量精度評價尤為重要。從匹配原理來說，由于匹配時較看重位置和樹高因素，而過分類對樹高影響不大，只增加了林木位置點，更容易將原本不能匹配的林木匹配在一起，當聚類數量增加甚至超過實際數量時，可以成功匹配實測數據的林木數量往往隨之增加，因此片面提升成功匹配林木比例或單獨關注聚類數量與實際數量的差異，都不能全面反映聚類效果。本研究設置3個精度指標，具體如下。

1)聚類數量精度(A-QT) 即單位林分中聚類出的單木數量與地面實測單木數量的差異性，代表聚類效果的數量精確程度，公式如下：

A-QT=1-|N-test-N-ref|/N-ref。

(6)

式中：N-test為聚類的單木數量；N-ref為地面實測單木數量。

2)聚類質量精度(A-QL) 即聚類出的單木中可與地面實測單木相匹配數量占總體的比例，代表聚類質量的好壞，公式如下：

A-QL=1-|N-test-N-match|/N-test。

(7)

式中：N-match為成功匹配的單木數量。

3)聚類有效精度(A-ED) 即單位林分中可匹配的聚類單木數量與地面實測單木數量的差異性，代表聚類效果相對實際情況而言有效部分所占比例，公式如下：

A-QL=1-|N-ref-N-match|/N-ref。

(8)

以數量精度、質量精度和有效精度三者的平均值作為總體聚類效果的評價指標，定義為聚類效果指數(A-k),代表聚類方法對單木的偵測和提取能力，公式如下：

A-k=(A-QT+ A-QL+ A-ED)/3。

(9)

2 結果與分析

2.1 不同分層方式對精度的影響

為了對比不同分層方式對精度的影響，暫將點云分為2層并設置分層界限高分別為3、4、5…12 m，統計聚類效果指數(A-k)、單木樹高精度(A-h)、單木冠幅精度(A-c)、林分平均高精度(A-H)、樹高基尼指數精度(A-Gc)如圖4a所示。可以看出，單木樹高精度、單木冠幅精度隨層高增大變動幅度較小，規律性不大；而樹高基尼指數精度、林分平均高精度、聚類效果指數隨層高增大呈單峰曲線變化，尤其是樹高基尼指數精度在50%～85%范圍內劇烈變化，三者的最佳點分別對應層高5、7和9 m，且存在此消彼長現象，無法選出最優情況，故將點云分為2層不能同時滿足各方面要求，需進行多層聚類。

圖4 不同分層方式對精度的影響

在多層聚類中，為了避免人為設置層高對精度的影響，統一采用均分百分位數方式進行分層，即保證每層點云數量相同。由4b可知，分層數對林分水平參數精度的影響很大，尤其是樹高基尼指數精度變動最大，由65%升至90%，林分平均高精度和聚類效果指數也有明顯提升；而單木樹高精度變化不明顯，單木冠幅精度顯著下降。

2.2 參數提取精度

由表3可知，無論分層數設置多少，單木樹高精度均在90%以上，可滿足林業調查對精度的要求，而單木冠幅精度較低，在60%左右，說明整體冠幅提取難度較大，需進一步改進；此外，聚類效果指數均在80%以上，最高可達88.28%，說明提取有效性較高，而林分平均高精度最高可達99.84%，樹高基尼指數精度最高可達89.65%，說明對林分整體信息和結構反演能力較強。

表3 分層數設置對聚類效果和精度的影響

為了解具體影響情況，按樹高組的分配狀態,以5 m為一級，分別對各樹高級的單木樹高精度和聚類效果指數進行分析。由圖5a可知，無論分層數設置多少，精度曲線大致趨勢相同，即10～20 m主林層單木樹高精度較高，達94%左右，而少數樹高達20 m以上的林木，由于點云密度過低，可能錯過樹頂點造成樹高精度較低。

對比不同分層數的單樹高精度和聚類效果(圖5b)可知，分層數對5 m以上主林層單木樹高精度影響不大，整體單木樹高精度差異主要是對5 m以下更新層林木的聚類。隨分層數不斷增加，對2～5 m更新層林木的提取能力顯著增強，由最低的34.86%升至74.43%，而其他林層聚類效果變動不大，直接導致整體LiDAR提取能力較顯著提升，但也由此帶來了單木樹高、冠幅精度的波動。更新層點云密度比主林層低，平均每株林木僅有4.9個點，聚類誤差引起1個點的聚類錯誤，即會對單木整體造成較大影響。對更新層林木的提取也是所有LiDAR單木提取研究中較為困難的一部分。

圖5 不同分層數對單木樹高精度和聚類效果的影響

2.3 提取有效性

以分2層聚類結果為例，將LiDAR提取出的林木與地面實測林木以5 m為一級，分別統計各樹高級范圍內的林木數量，計算二者比值代表該樹高級的提取率，超過100%的部分視為無效。不同層高對應的各樹高級提取率如圖6所示。

圖6 層高3～12 m時各樹高級提取率

以7、10 m為界，各樹高級提取率大致分2種情況：1)分層高小于7 m或大于等于10 m時，10 m以下林木提取率較低，10 m以上林木提取率虛高，均超過100%，最高甚至達194%，即有近一半林木提取錯誤，沒有實際意義；而決定最終聚類效果指數的往往是10 m以下林木提取的有效性；2)分層高在7～10 m之間，情況正相反?？傮w來說，提取率曲線在100%以下的面積可代表整體提取有效性高低，如以7 m為分層高，各樹高級提取率較為均衡，有效性最高。

2.4 林分結構信息反演

對聚類結果和地面實測林木按樹高統計數量，結果如圖7所示?？梢钥闯?，不分層時樹高近似呈正態分布，而實際樹高分布為雙峰曲線，聚類結果與實際情況嚴重不符；相比之下，分層聚類對樹高分布的預測能力更為突出，尤其對10 m以上主林層，無論是數量還是比例均與實際情況十分吻合，而樹高2～8 m時還存在一定誤差，可通過提升對下層林木的提取能力來攻克這一難題。此外由上述可知，林分平均高精度最高達99.84%，樹高基尼指數精度最高達89.65%，樹高變異系數最高達89.23%，同時計算出最鄰近體間平均距離為1.53 m，與實際的1.31 m較為接近，均體現了分層聚類算法在林分結構反演上的突出優勢。

圖7 聚類結果與地面實測林木樹高分布情況對比

3 討論

3.1 LiDAR提取林木的有效性及評價指標

以往LiDAR單木提取研究中，評判算法優劣時較看重提取率、匹配精度和單木樹高精度等指標。以分水嶺分割為代表的分割表面模型算法以及不分層、直接分割的相關算法由于對林下層點云利用率低，造成林木分割數量少于實際數量，而基于LiDAR提取林木只要不超過樹高、距離閾值均視為有效，因此LiDAR提取林木的有效性問題并不突出。但本研究采用的分層K均值聚類算法，很有可能因過分類出現LiDAR提取林木數量多于實際數量的情況，在實際應用過程中，如果沒有地面實測數據匹配驗證，有效林木和無效林木將無法區分，并可能進一步影響整體結果的有效性。

以往研究表明，林木提取率(提取數量與實際數量比)在20%～60%之間不等(Morsdorfetal.，2004；Kandareetal., 2016；尹艷豹，2010)，大多數研究更加關注提取出并匹配的林木在三維信息方面的精度，忽視了提取有效性問題，即沒有驗證提取出的林木究竟是主林層、優勢木還是包括多少比例被壓木，因此盡管提取出的部分林木信息很精確，但由于不清楚其與整體的關系，不能夠確定這部分信息是否可以代表總體情況，在實際生產應用中意義不大。在LiDAR單木提取研究中，追求能夠“更準確”地提取出“更多”的林木，其中“更準確”以單木信息精度來反映，而“更多”則以提取有效性來體現，即本研究中設置的“聚類效果指數”，用該指數綜合判斷“LiDAR提取了多少林木”“提取的林木中有多少比例是正確的”等問題，充分考慮了結果的可用性。聚類效果指數越高，提取有效性越高，只有達到一定標準，才能說明算法可以擺脫對地面實測數據的依賴，能夠直接應用于生產實踐中。

3.2 LiDAR點云密度對單木及林分結構信息的影響

點云密度對信息提取的影響是近年來比較熱門的研究內容之一。在單木信息精度方面，研究指出點云密度高于20 m-2帶來的精度提升微乎其微(Reitbergeretal.，2009；Kaartinenetal.，2012；Yaoetal.，2014)，點云密度2、4、8 m-2(Kaartinenetal.，2012)的單木信息提取精度差異不顯著，因此，很多后續研究采用先稀釋點云再進行信息提取的方法(Magnussenetal.，2010；Hansenetal.，2015)。本研究中，林冠層平均點云密度為2.86 m-2，所得單木樹高精度達90%以上，單木冠幅精度達60%以上，與大多數研究結果一致，驗證了該密度水平的點云數據可以很好完成單木信息提取工作。

然而在林分結構方面，有研究對60 m-2點云先進行不同程度稀釋后再反演林分結構信息，包括林分密度、水平結構和垂直結構，結果發現60 m-2與6 m-2點云的效果存在顯著差異，極高密度點云在林分結構反演方面更具有優勢(Kandareetal., 2016)。本研究中，通過點云多次分層再聚類的方法反演林分結構信息也得到很好效果，當分層數大于8時，樹高基尼指數、樹高變異系數精度在85%以上，不以增加10倍信息量為代價也能夠較好反演森林結構參數。

在提取有效性方面，目前還未有研究專門探討點云密度對提取有效性的影響。本研究中，通過點云分層可提升對下層林木的提取能力，但受困于下層林木的點云密度過低，平均每株林木僅有4.9個點云，在聚類時1個點分類錯誤就有可能引起整體信息的劇烈變化，因此猜測如能提高點云密度，將會提升對下層林木的提取能力，從而提升整體的提取有效性。5 m以上的單木樹高精度達90%以上，提取有效性也在90%以上，可反映整體林分除更新層外其他林木的生長情況；如果能夠通過提高點云密度從而攻克更新層林木的提取問題，將在森林精準量測、林木生長模型構建、可視化經營等方面起到一定的推進作用。

4 結論

本研究提出一種基于機載LiDAR點云多層聚類的單木信息提取方法，并將其應用于1塊林木密度高、空間結構異質性大的復層林進行驗證，同時對比不同匹配方法對同一聚類結果的匹配精度。結果表明：1)分多層次K均值聚類進行單木提取具有很好效果，尤其在下層林木的提取方面表現突出；2)可在低點云密度條件下完成工作，降低高密度點云數據信息的冗余；3)在單木信息精度方面，分層K均值聚類提取的單木樹高精度達90%以上，單木冠幅精度在60%左右；4)在林分結構方面，采用多次分層方式可提升對林分結構的反演能力，林分垂直結構、水平結構、樹高分布結構等方面均有較好反映；5)通過構建LiDAR林木提取有效性指標，著重關注算法的適用性，提取林木有效性最高達88.70%，進一步說明該方法可以擺脫對地面實測數據的依賴，能夠應用于生產實踐中。