高郁閉度人工林無人機激光雷達單木分割方法優化*

朱泊東 羅洪斌 金 京 岳彩榮

(西南林業大學林學院 昆明 650224)

森林覆蓋地球表面的1/3，是世界上一半以上陸地生物的家園，也是全球75%淡水的來源，在碳循環、水循環、生物多樣性保護和可持續發展中發揮著至關重要的作用(Siepielskietal.， 2017)。單木是森林的基本結構單元，其三維結構、生長狀況、空間分布等特性是森林資源調查、監管和生態建模研究等所需的關鍵因子(李增元等， 2016)，準確的單木分割是獲取精細化森林參數的前提條件(Meietal.， 2004； Korpelaetal.， 2007； Zhaoetal.， 2018； Nuijtenetal.， 2019)。

傳統森林資源調查多以人工為主，存在投入大、成本高、效率低、時空機動性差等明顯缺點。激光雷達(light detection and ranging， LiDAR)是一種快速獲取目標地物三維空間信息的主動遙感技術，其通過發射激光脈沖到樹木的枝葉等組分上，接收包含三維空間位置等信息的散射回波信號，并可生成點云數據，對精確估算冠層三維結構信息尤其在林分尺度和單木尺度的森林垂直結構探測中具有顯著優勢(Pedersenetal.， 2012； 劉魯霞等， 2014)。目前，基于機載LiDAR數據進行單木分割的方法大致可分為2類。一是冠層高度模型(canopy height model， CHM)單木法，利用2D柵格化CHM提取單木特征(Brandtbergetal.， 2003； Brandtberg， 2007； Kochetal.， 2006； Popescu，2007； Yuetal.， 2011； 劉清旺等， 2010)，其基本思路是將點云分為地面點和植被點，采用空間插值方法生成相同空間分辨率的數字表面模型(digital surface model， DSM)和數字高程模型(digital elevation model， DEM)，二者相減得到CHM，通過分析CHM局部最大值判斷樹冠頂點，再依據樹冠頂點提取樹冠邊界(Popescuetal.， 2003； Johanetal.， 2003； Maltamoetal.， 2004)，或先提取樹冠邊界，再將樹冠內的高度最大值作為樹冠頂點(Brandtbergetal.， 2003； Zimbleetal.， 2003)。常見的CHM分割方法包括分水嶺分割(watershed segmentation)(Tangetal.， 2007)、標記控制分水嶺(marker controlled watershed)(Chenetal.， 2006； 郭昱杉等， 2016)、區域增長法分割(regional growth segmentation)(Hyypp?， 1999)等，該類方法的缺點是僅利用冠層表面和地面信息，對于復雜林分的亞冠層、中下層喬木和灌木存在探測和提取困難問題，且對點云數據的三維信息挖掘深度不夠，點云柵格化處理過程中會造成點云空間信息丟失，分割精度也受空間插值方法和柵格分辨率大小影響(霍朗寧等， 2021)。二是歸一化點云(normalized point cloud， NPC)單木法，利用NPC數據，根據點云之間的空間位置關系設定判別準則對點云進行聚類，從而實現單木分割并提取單木特征(Alexander， 2009； Reitbergeretal.， 2009； Lietal.， 2012)。常見的NPC分割方法包括k均值聚類(Zimbleetal.， 2003)、均值漂移聚類(Ferrazetal.， 2012)、自適應距離聚類(Leeetal.， 2010)等，該類方法的缺點是在高郁閉度林分條件下對中下層林木提取能力較弱，分割精度受點云密度、郁閉度以及林分中每株樹水平分布和垂直結構的影響，表現為林分密度或郁閉度較大時分割精度較低(Lietal.， 2012)，郁閉度較小時分割精度較高(Falkowskietal.， 2008； Kaartinenetal.， 2012)。Forzieri等(2009)研究表明，在株數密度較高的林分內，現有單木分割方法對冠層結構復雜和郁閉度高的林分分割效果均不佳。Vauhkonen等(2011)采用6種單木分割方法對針葉林、闊葉林和混交林進行分割，結果發現不同方法在不同林分條件下的分割效果差異較大。Yang等(2019)探討植被特征(森林類型、葉面積指數、冠層覆蓋度、株數密度和高度變異系數)對單木分割(CHM分割、無坑CHM分割、基于點云分割和基于層堆疊種子點分割)精度的影響，得出植被特征差異會對分割效果產生較大影響，多數植被特征均與分割精度呈負相關，但其機制尚不明確。以往研究大多基于單一分割方法對一種或多種復雜森林場景進行單木分割，視森林參數變化采用不同分割方法的單木分割報道較少； 且多數研究雖有對分割方法進行參數優化和方法改進使之更適用于某一林分類型，但不同分割方法在不同林分條件下的分割效果各有優劣，某一林分條件下采用與之對應的分割方法能夠進一步提升分割精度； 同時，據目前報道，對高郁閉度林分條件下亞熱帶人工針葉林內同一樹種在不同齡組下單木分割精度差異的相關研究仍較少。

無人機LiDAR技術是獲取精細化森林參數的重要技術手段，利用LiDAR點云數據進行單木分割是提取精細化森林參數的基礎，森林參數獲取精確與否取決于單木分割是否準確。針對高郁閉度林分條件下基于LiDAR點云數據單木分割林木提取困難、總體精度較低等問題，本研究以云南思茅松(Pinuskesiyavar.langbianensis)人工林為研究對象，以無人機LiDAR點云數據為基礎，開展多種點云分割方法試驗，提出一種基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法，以期為開展無人機LiDAR技術森林資源調查提供技術參考，為提高單木分割總體精度提供新策略。

1 研究區概況與數據

1.1 研究區概況

萬掌山林場位于云南省普洱市思茅區北部，思茅區地處云南省南部、普洱市中南部、瀾滄江中下游，北鄰無量山，東連江城縣，22°27′—23°06′N、100°19′—101°27′E。思茅區屬低緯高原南亞熱帶季風氣候區，境內立體氣候明顯，有北熱帶、南亞熱帶、中亞熱帶和北亞熱帶4種氣候類型，具有低緯、高溫、多雨、濕潤、靜風的特點，年均氣溫17.8 ℃，年均降雨量1 524.4 mm，無霜期318 天，冬季嚴寒，夏季酷暑，四季溫和(歐光龍等， 2014)。萬掌山林場占地面積1.95萬 hm2，其中包含生態公益林4 777 hm2、商品林5 866 hm2(潘瑤等， 2016)，以思茅松人工純林為主。研究區位置及樣地分布如圖1所示。

1.2 遙感數據

以縱橫大鵬CW-30油電混合中型垂直起降固定翼無人機為遙感平臺，平臺搭載CW-30 LiDAR離散回波形激光雷達傳感器和JOUAV CA-102 4200萬像素全畫幅航攝相機，用于獲取LiDAR數據和DOM影像。無人機飛行高度900 m，航寬180 m，旁向重疊度≥60%。激光雷達脈沖發射頻率820 kHz，點云平均密度19.21 pts·m-2，DOM影像空間分辨率0.1 m，DEM空間分辨率1 m。機載LiDAR數據獲取時間為2020年1月5日。

1.3 目視解譯樣地數據

基于無人機高分辨率影像，結合激光點云的三維形態并實地踏查輔助視覺驗證確定目視解譯的單木樹冠中心位置，以此作為驗證數據。在研究區內設置35塊目視解譯樣地(圖1)，包含幼齡林、中齡林和近熟林3個齡組，其中幼齡林樣地13塊、中齡林樣地14塊、近熟林樣地9塊，共1 427株樹，樣地大小均為20 m×20 m。本研究采用的目視解譯結果僅包含上層林木。

2 研究方法

將高郁閉度思茅松林分按齡組劃分為幼齡林、中齡林和近熟林，采用分水嶺算法、基于點云的局部最大值聚類算法進行單木分割，對比2種方法的分割效果，找到不同林分條件下的最佳分割方法和適用條件； 在此基礎上提出了基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法，并分析分水嶺算法中CHM空間分辨率和DSM插值方法對分割效果的影響。ULS數據單木分割方法對比及精度分析技術路線如圖 2 所示。

圖2 ULS數據單木分割方法對比及精度分析技術路線Fig.2 Comparison of individual tree segmentation methods for ULS data and technical route of accuracy analysis

2.1 基于點云數據的CHM生成

在對LiDAR點云數據進行去噪、地面點分類(濾波)的基礎上，用地面點插值生成DEM，用非地面點插值生成DSM，試驗中CHM由不同空間分辨率下采用不同插值方法(克里金插值、反距離權重插值、不規則三角網插值)生成的DSM與同一DEM作差得到。克里金(Kriging)插值是一種基于統計學的插值方法，其計算優化的協方差，利用高斯過程進行插值，插值的整個過程相當于對未知區域化變量的期望值進行加權滑動平均，并將該值作為柵格單元值。反距離權重(inverse distance weighting， IDW)插值基于相近相似原理，利用附近點計算柵格單元值，并以點距柵格單元中心點距離冪次方的倒數為權重進行加權平均，距柵格單元中心點距離越近的點被賦予的權重越大。不規則三角網(triangulated irregular network， TIN)插值基于地面點數據的三角形不規則網絡推導每個三角形的雙變量函數，并將其用于估計未采樣位置(已知地面點之間)的值，每個生成的三角形平面用于插值。

2.2 單木分割方法

2.2.1 分水嶺算法 分水嶺算法(watershed algorithm)是一種基于拓撲理論的數學形態學分割方法(Meyeretal.， 1990； Vincentetal.， 1991)，其基本思想是把圖像看作測地學上的拓撲地貌，圖像中每一點像素的灰度值表示地形起伏高低，對灰度值取反； 把圖像中一些局部極小值點看作盆地底部，在每個盆地底部鉆一個洞，此時地下涌出的水不斷上升并逐漸淹沒整個地形； 當2個相鄰盆地中的水即將匯合到一起時，在水交匯的邊緣上修筑一道堤壩阻止盆地間水的交匯，該堤壩即對應的分割線，當眾多堤壩將每盆地都包圍時，堤壩的集合就是分水嶺，即圖像分割邊界(韓懸等， 2019)。分水嶺算法對微弱邊緣具有較強的響應，由于圖像質量不斷提高，圖像紋理特征逐漸增多，圖像中的噪聲、物體表面細微灰度變化均有可能產生過分割現象，因此柵格分辨率大小對圖像分割效果具有較大影響(孫釗等， 2021)。本研究采用分水嶺算法先確定每株樹的樹冠邊界范圍，再在該范圍中尋找局部最大值，將局部最大值看作樹冠頂部(Wulderetal.， 2000)。

2.2.2 基于點云的局部最大值聚類算法 基于點云的局部最大值聚類算法假定樹木之間存在一定間距，通過分析點的高程以及該點到其他點的距離，自上而下地對點云進行聚類； 將歸一化點云中的全局最高點作為最高木頂點，基于間距閾值和最小間距判別準則進行迭代判斷，設置一個閾值，將離目標單木頂點水平距離小于閾值的點云作為目標單木的點云集群，反之則移除該集群，最終實現單木分割(Lietal.， 2012)。

2.2.3 基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法 高郁閉度林分條件下，齡組差異對單木分割精度具有一定影響，考慮到實際應用中大范圍、精確的齡組信息獲取難度較大，本研究采用相關分析方法挖掘反映齡組差異的激光雷達點云特征。首先基于激光點云數據提取59個點云特征，包括高度變量(46個)、點云強度變量(10個)、郁閉度、葉面積指數和間隙率，然后分析點云特征與分割精度的相關性，結果發現，在59個點云特征中，冠層起伏率(canopy relief ratio， CRR)與分割精度顯著相關，故本研究以冠層起伏率差異作為分割方法選擇依據。冠層起伏率(Parkeretal.， 2004)計算公式如下：

(1)

式中： mean為每一統計單元內所有點的平均高度； min為每一統計單元內所有點的最小高度； max為每一統計單元內所有點的最大高度。

對機載LiDAR點云數據進行去噪、濾波、點云分類和歸一化，提取歸一化點云的冠層起伏率，采用四分位數法將冠層起伏率由低到高分為4級，每一級應用對應分割方法進行單木分割。

2.3 精度驗證

以探測率、準確率和F得分為指標對單木分割精度進行驗證和評價(Goutteetal.， 2005； Sokolovaetal.， 2006)，計算公式如下：

(2)

(3)

(4)

式中： TP表示被正確探測到的樹的數量(正確分割)； FN表示未被探測到的樹的數量(欠分割)； FP表示被算法識別到但實際不存在的樹的數量(過分割)；r為單木探測率，p為被探測到單木的正確率，r和p可以近似看作制圖精度和用戶精度；F得分為綜合考慮r與p計算得到的加權調和平均值。

3 結果與分析

3.1 基于齡組的單木分割效果

圖3所示為分水嶺算法和基于點云的局部最大值聚類算法在不同齡組樣地的分割效果及其精度驗證，表1所示為3個齡組樣地的基本概況及2種分割方法在各齡組對應的單木分割精度。

圖3 2種分割方法在不同齡組樣地的分割效果及其精度驗證Fig.3 Results and accuracy of two methods of individual tree segmentation for the forest with different age groups

由表1可知，分水嶺算法的F得分由高到低分別為幼齡林(0.78)、中齡林(0.68)和近熟林(0.67)，基于點云的局部最大值聚類算法的F得分由高到低分別為近熟林(0.77)、中齡林(0.66)和幼齡林(0.64)，分水嶺算法和基于點云的局部最大值聚類算法在各齡組對應的單木分割精度不同且存在較大差異，但2種方法在分割精度的變化趨勢上剛好相反，具有互補性。

表1 3個齡組樣地的基本概況及2種分割方法在各齡組對應的單木分割精度Tab.1 The basic profile of the sample plots of the 3 age groups and the segmentation accuracy of individual trees corresponding to the 2 segmentation algorithms in each age group

不同年齡階段，思茅松樹冠形狀不同，冠層表面起伏程度不同，林分內部結構也存在較大差異。在幼齡林中，分水嶺算法相較基于點云的局部最大值聚類算法有著更高的單木探測率和準確率，這是因為幼齡林大多為高度差異較小、株數密度較高的單層林，冠層點云幾乎處于同一水平面，通過探測局部最大值點云的基于點云的局部最大值聚類算法容易造成欠分割，導致分割精度大大降低，而基于拓撲理論的分水嶺2D形態學分割方法通過對CHM進行濾波增強圖像紋理有助于提升分割精度，但當林分郁閉度過高時2種分割方法均無法得到較好的分割效果，錯分割和欠分割現象非常嚴重。在中齡林中，分水嶺算法與基于點云的局部最大值聚類算法的F得分差異不明顯，分水嶺算法相較基于點云的局部最大值聚類算法單木探測率更高，但準確率較低，基于點云的局部最大值聚類算法與之相反。在近熟林中，基于點云的局部最大值聚類算法相較分水嶺算法對上層喬木有著更高的探測率和準確率，這是因為近熟林大多為高度差異較大的復層林，樹冠起伏明顯，但樹木之間冠層交織嚴重，且許多思茅松近熟林樹枝分叉較多，分水嶺2D形態學分割方法利用CHM進行局部最值探測時易將樹枝識別為樹冠導致過分割，而基于點云的局部最大值聚類算法(3D分割方法)利用點云空間分布對點云進行迭代聚類，能夠有效降低將樹枝錯分為樹冠的概率。

3.2 基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法

從圖4可以看出，高郁閉度思茅松幼齡林點云大多分布在冠層，導致樣地內點云平均高度mean上移，由式(1)可知幼齡林的冠層起伏率最大； 中齡林和近熟林因其垂直結構復雜，點云在垂直方向上分布較為均勻，不會導致mean過大，近熟林的冠層起伏率略高于中齡林。

圖4 不同齡組樣地點云剖面Fig.4 Point cloud profile of different age groups sample

為更好揭示不同齡組樣地間冠層起伏率的分布規律，構建小提琴圖(圖5)。小提琴圖基本上是核密度圖以鏡像方式在箱線圖上的疊加，白點表示中位數，白色盒型范圍為下四分位點到上四分位點，細白線表示須，外部形狀即核密度估計。圖5顯示，思茅松幼齡林和近熟林的冠層起伏率近似正態分布，中齡林的冠層起伏率近似雙峰偏態分布，對比三者分布情況，幼齡林冠層起伏率的中位數、上四分位數和下四分位數均高于中齡林和近熟林； 核密度估計結果亦顯示，幼齡林與中齡林和近熟林存在顯著差異； 近熟林冠層起伏率的下四分位數、中位數、上四分位數均高于中齡林，說明隨著年齡增加林分內點云的垂直結構分布發生改變。

圖5 3個齡組樣地冠層起伏率的小提琴圖Fig.5 Violin plot of canopy relief ratio of three age groups in sample plot

圖4、5分析表明，冠層起伏率與齡組之間存在較好的對應關系，以齡組作為分層分割參數可以有效提高分割精度，但在實際應用中，一些區域大范圍精細化的齡組信息獲取難度較大，且因立地條件差異常導致相同齡組的林分分割效果也存在較大差異，為此，本研究提出對冠層起伏率進行分級代替齡組作為分層分割的重要依據。根據試驗結果，采用四分位數法可將冠層起伏率由低到高分為4級： 1級對應中齡林，2級對應近熟林，3～4級對應幼齡林。冠層起伏率1～2級樣地采用基于點云的局部最大值聚類算法進行單木分割，冠層起伏率3～4級樣地采用分水嶺算法進行單木分割，幼齡林與中齡林、近熟林劃分的正確率分別為0.67和0.94，幼齡林劃分結果包含12塊幼齡林樣地、2塊中齡林樣地和4塊近熟林樣地； 中齡林、近熟林齡組劃分結果包括1塊幼齡林樣地、11塊中齡林樣地和5塊近熟林樣地。

為評價基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法的分割效果，本研究對比不分層單一分割方法(分水嶺算法和基于點云的局部最大值聚類算法)的分割精度(表2)，基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法的分割精度(0.75)相較分水嶺算法(0.71)或基于點云的局部最大值聚類算法(0.68)有明顯提升。

表2 單木分割精度驗證Tab.2 Accuracy assessments of individual tree segmentation

3.3 CHM空間分辨率和插值方法對分割精度的影響

分水嶺算法基于CHM進行圖像分割，為評價不同方法生成的CHM對分割精度的影響，采用反距離權重、克里金和不規則三角網3種常用空間插值方法得到的DSM生成CHM，選取0.2 m×0.2 m、0.5 m×0.5 m、0.8 m×0.8 m和1.0 m×1.0 m 4種CHM空間分辨率進行單木分割，并對比其分割精度(圖6)。

在3種常用空間插值方法中，當CHM空間分辨率為0.5 m×0.5 m時，反距離權重插值得到DSM生成的CHM分割精度最高(0.81)； 當CHM空間分辨率為0.2 m×0.2 m時，分割精度由高到低依次為反距離權重(0.44)、不規則三角網(0.42)和克里金(0.40)； 當CHM空間分辨率為0.5 m×0.5 m時，分割精度由高到低依次為反距離權重(0.81)、不規則三角網(0.80)和克里金(0.77)； 而當CHM空間分辨率為0.8 m×0.8 m和1.0 m×1.0 m時，不規則三角網插值得到的分割精度最高，由高到低依次是不規則三角網、反距離權重和克里金，0.8 m×0.8 m時F分別為0.59、0.47和0.47，1.0 m×1.0 m時F分別為0.46、0.37和0.37(表3)。

表3 分水嶺算法中CHM空間分辨率和插值方法對分割精度的影響Tab.3 Influence of CHM resolution and interpolation method on segmentation accuracy in watershed algorithm

圖6 3種空間插值方法在不同分辨率下的分割精度Fig.6 Segmentation accuracy of three spatial interpolation methods at different resolutions

4 討論

基于高郁閉度思茅松林分無人機激光雷達數據提取59個點云特征，分析點云特征與分割精度的相關性發現，在59個點云特征中，除冠層起伏率外，高度變異系數和累計高度百分位數等垂直結構特征與齡組也表現出一定相關性，但冠層起伏率對齡組的區分度優于其他特征參數。研究曾對樣地冠層起伏率嘗試多種分級方式，如相等間隔法、自然斷點分級法、分位數法和標準差法等，結果表明采用分位數法將冠層起伏率分為4級與樣地齡組的一致性最好。

分水嶺算法探測準確率較低，過分割現象較多； 基于點云的局部最大值聚類算法探測率較低，欠分割現象較嚴重，該結果與李平昊等(2018)的研究結論一致，表明分水嶺算法和基于點云的局部最大值聚類算法單木分割適用的林分條件不同，并表現出一定互補性，這也是本研究將2種算法相結合的主要依據。

研究發現，當林分郁閉度為1時，冠層起伏率變化出現2種情況： 一是當點云最小高度為0時，冠層起伏率偏大，統計單元內有少數激光點云透過樹冠間間隙打到地面上，大部分點云集中在樹冠上部，統計單元內所有點云平均高度偏大，導致點云平均高度與最小高度的差變大，即式(1)分子變大，而統計單元內最大高度和最小高度的差不變，即式(1)分母不變，致使冠層起伏率出現異常增大； 二是當統計單元內沒有激光點云打到地面上時，冠層起伏率變化規律與齡組無關，此時統計單元內點云平均高度和最小高度均變大，即式(1)分子和分母同時變大，冠層起伏率可能變大、變小或不變，其變化趨勢主要取決于統計單元內最低點高度。此外，當郁閉度為1時，上述2種分割方法無法準確進行單木探測和提取，故不討論其分層分割精度。

本研究主要對高郁閉度林分上層林木進行探測和分割，未涉及中、下層林木，且受點云數據限制，僅針對思茅松林分開展試驗，其分割優化方法是否適用于其他森林類型有待進一步研究。分割精度方面，如果能將機載CCD相機獲取的高分辨DOM影像與LiDAR數據相結合，在高密度點云基礎上加入高空間分辨率的光譜信息或紋理特征，則有望提升其單木分割精度。此外，應引入更先進的單木分割方法，如深度學習等。

因缺少足夠數量具有精確地理坐標的實測單木檢驗樣本，本研究單木分割精度檢驗主要通過無人機DOM影像目視判讀進行，這一定程度上影響了評價結果的可靠性。林內精確定位單木位置坐標一直是森林調查中的瓶頸，如果能解決高郁閉度林分林下單木定位精度低的問題，將對準確評價單木分割精度、改進分類方法具有較好促進作用。

5 結論

針對高郁閉度林分條件下基于LiDAR點云數據單木分割林木提取困難、總體精度較低等問題，本研究比較分水嶺算法和基于點云的局部最大值聚類算法的單木分割精度，提出基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法。通過提取樣地冠層起伏率，確定分水嶺算法和基于點云的局部最大值聚類算法單木分割適用的林分條件，可拓寬單一單木分割方法在不同林分條件下的優勢，提升單木分割精度。此外，還通過改變柵格化冠層高度模型(CHM)的空間分辨率和數字表面模型(DSM)的空間插值方法，對分水嶺算法單木分割精度敏感性進行分析。結果表明：

1) 在幼齡林中，冠層起伏率較大，分水嶺算法對單層林的分割效果優于基于點云的局部最大值聚類算法； 在中齡林和近熟林中，冠層起伏率較小，分水嶺算法易將思茅松樹枝識別為樹冠，基于點云的局部最大值聚類算法的分割效果優于分水嶺算法；

2) 基于冠層起伏率結合分水嶺算法和基于點云的局部最大值聚類算法的分層分割法充分考慮不同齡組的林分結構差異，精度最高(F= 0.75)，優于分水嶺算法(F= 0.71)和基于點云的局部最大值聚類算法(F= 0.68)；

3) 在分水嶺算法中，當分辨率為0.5 m×0.5 m時采用反距離權重法(IDW)插值得到的CHM單木分割精度最高(r= 0.70，p= 0.94，F= 0.81)，高于克里金插值(F=0.80)和不規則三角網插值(F=0.77)。