陳思宇，劉憲釗，王懿祥，梁 丹

（1.浙江農(nóng)林大學(xué) 環(huán)境與資源學(xué)院，浙江 杭州 311300；2.浙江農(nóng)林大學(xué) 省部共建亞熱帶森林培育國家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300；3.浙江農(nóng)林大學(xué) 浙江省森林生態(tài)系統(tǒng)碳循環(huán)與固碳減排重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300；4.中國林業(yè)科學(xué)研究院 資源信息研究所 國家林業(yè)和草原局森林經(jīng)營與生長模擬重點(diǎn)實(shí)驗(yàn)室，北京 100091）

華北落葉松Larix principis-rupprechtii是中國華北地區(qū)第一大樹種，具有較高的經(jīng)濟(jì)和生態(tài)價(jià)值，探討華北落葉松單木參數(shù)的提取、林分蓄積量的估算具有重要意義[1]。高精度的單木識別是獲取單木和森林結(jié)構(gòu)參數(shù)的重要前提，單木識別包括單木位置識別和樹冠分割[2-4]。相比高分辨率衛(wèi)星影像和光學(xué)影像，機(jī)載LiDAR技術(shù)不受區(qū)域信號影響，分辨率較高，是一種主動遙感技術(shù)，能獲取林木的水平和垂直結(jié)構(gòu)信息[5-6]。利用機(jī)載LiDAR研究單木識別的方法主要有2類：基于冠層高度模型(CHM)識別方法和基于歸一化后的點(diǎn)云識別方法[3]。其中基于CHM的單木位置識別方法即在灰度圖像中探測單木樹頂，主要實(shí)現(xiàn)算法包括局部最大值法和多尺度分割法。局部最大值法假設(shè)樹冠頂點(diǎn)高于樹冠邊緣，通過設(shè)置固定尺寸[7-10]或動態(tài)尺寸[11-12]的搜索窗口從CHM中識別樹冠頂點(diǎn)特征，單木位置探測率可達(dá)85%。多尺度分割法以像元為單位，根據(jù)相似特征逐漸外擴(kuò)分割，關(guān)鍵在于最優(yōu)分割尺度的選取[13-17]。基于歸一化后的點(diǎn)云單木位置識別方法以點(diǎn)云的結(jié)構(gòu)信息為基礎(chǔ)，主流實(shí)現(xiàn)算法包括區(qū)域生長法和均值漂移法。區(qū)域生長法通過選擇合適的種子點(diǎn)，定義距離生長準(zhǔn)則，判別不同屬類[18-20]，選擇的種子點(diǎn)和制定的距離生長規(guī)則不同，取得的單木識別精度也不同。均值漂移法是一種基于密度梯度上升的非參數(shù)方法，通過迭代運(yùn)算求得同簇的數(shù)據(jù)點(diǎn)[21]，關(guān)鍵在于漂移均值向量的計(jì)算。單木識別時(shí)可利用林分樹冠的距離和多維特征空間為向量來完成[22-23]，其在中低郁閉度林分中的單木識別精度較高，但在高郁閉度林分條件下，由于樹冠重疊度增加，漂移均值向量計(jì)算會有偏差，識別結(jié)果不盡合理。為提高識別精度，結(jié)合樹冠冠型特點(diǎn)，本研究提出基于機(jī)載LiDAR點(diǎn)云空間特征的高斯核函數(shù)改進(jìn)的均值漂移單木位置識別方法(MSP)。為驗(yàn)證方法的有效性，對比了基于點(diǎn)云空間特征的區(qū)域生長圖像分割方法(RGP)、基于CHM的局部最大值單木探測方法(LMC)和基于CHM的多尺度分割方法(MSC)等3種目前應(yīng)用較為廣泛的單木位置識別方法。最后，在單木位置識別結(jié)果上，進(jìn)行單木樹冠的提取。

1 研究區(qū)概況

忻州市五寨縣羊圈溝林場 (38°44′～39°01′N，111°48′～112°03′E)位于山西省中北部，屬亞寒帶大陸性氣候，季節(jié)變化明顯，海拔為1 350～2 700 m，年均氣溫為5 ℃，年降水量為453.9 mm；華北落葉松-云杉Picea asperata天然次生林與華北落葉松人工林帶狀混交。本研究在華北落葉松人工林帶內(nèi)設(shè)置8個(gè)20 m×20 m 的正方形樣地，植株密度為 1 875～2 925 株·hm-2，郁閉度為 0.8～1.0，是典型高郁閉度森林。

2 數(shù)據(jù)采集與研究方法

2.1 數(shù)據(jù)采集

標(biāo)記并編號8個(gè)樣地內(nèi)胸徑大于5 cm 的華北落葉松，調(diào)查記錄每木冠幅、樹高和胸徑；利用實(shí)時(shí)動態(tài)定位和全站型電子測距儀獲取每木的位置坐標(biāo)。樣地林分基本情況如表1。

表1 樣地林分基本特征Table 1 Investigation results of basic characteristics of forest stand in sample plot

2018年6月2日，使用GV 1500六旋翼無人機(jī)搭載RIGEL VUX-1激光掃描儀獲取機(jī)載LiDAR點(diǎn)云數(shù)據(jù)；設(shè)置飛行高度為200 m，飛行速度為6 m·s-1，離散點(diǎn)云平均密度為70個(gè)·m-2，坐標(biāo)系為WGS_1984_UTM_Zone_49N。

在同一天，使用GV 120六旋翼無人機(jī)搭載傳感器SONY Alpha7 RⅡ采集可見光遙感影像。通過Limapper V 2.0軟件生成數(shù)字正射影像作為輔助數(shù)據(jù)，用于目視解譯與可視化數(shù)據(jù)展示。

2.2 數(shù)據(jù)預(yù)處理

首先，對原始點(diǎn)云數(shù)據(jù)進(jìn)行2 m緩沖帶裁剪，獲取8個(gè)樣地。然后，對樣地點(diǎn)云數(shù)據(jù)進(jìn)行以下處理：①采用距離閾值法移除噪聲點(diǎn)。②采用漸進(jìn)加密三角網(wǎng)濾波算法提取地面點(diǎn)云。為提高后期生成數(shù)字高程模型(DEM)的精度，設(shè)置迭代距離與迭代角度為1.4 m和8°，以便分離地面點(diǎn)云和非地面點(diǎn)云。③對去噪和分類后的點(diǎn)云數(shù)據(jù)進(jìn)行克里金插值[23]。通過定義該點(diǎn)鄰域中采樣點(diǎn)的距離加權(quán)平均計(jì)算像元值。對地面點(diǎn)云和非地面點(diǎn)云分別插值，得到DEM和數(shù)字地表模型(DSM)，兩者相減獲得冠層高度模型(CHM) (圖1)。圖1C中，高程值較低的條帶為華北落葉松人工林帶，較高的為華北落葉松-云杉天然次生林。④基于LiDAR360 V 3.2軟件，借助DEM對去噪后的點(diǎn)云進(jìn)行高度信息歸一化，避免地形起伏對點(diǎn)云高程的影響。

圖1 研究數(shù)據(jù)Figure 1 Study data

2.3 單木位置識別

2.3.1 基于點(diǎn)云空間特征的單木位置識別 均值漂移算法是一種基于非參數(shù)概率密度梯度估計(jì)的迭代聚類過程[24-25]。單木點(diǎn)云分布具有與高斯核函數(shù)相似的幾何形狀特征，可用高斯核函數(shù)改進(jìn)的均值偏移算法對點(diǎn)云數(shù)據(jù)進(jìn)行聚類，進(jìn)而獲取單木位置。設(shè)置迭代運(yùn)算的收斂閾值(ε)，基于點(diǎn)云空間特征的改進(jìn)均值漂移單木位置識別方法，算法如下：(1) 在點(diǎn)云數(shù)據(jù)點(diǎn)中隨機(jī)選擇1個(gè)點(diǎn)作為初始中心點(diǎn)，找出離中心點(diǎn)距離在帶寬h之內(nèi)的所有點(diǎn)，記做集合M，認(rèn)為這些點(diǎn)屬于同一類簇；(2)選用點(diǎn)云三維空間特征，計(jì)算中心點(diǎn)到集合M中每個(gè)數(shù)據(jù)點(diǎn)的向量之和，得到偏移向量Mh(x)；(3)令中心點(diǎn)沿偏移向量方向移動，移動距離是偏移向量Mh(x)的模；(4)重復(fù)步驟(2)～(3)，直到偏移向量小于ε，記住此時(shí)的中心點(diǎn)；(5)重復(fù)步驟(1)～(4)，直到所有數(shù)據(jù)點(diǎn)都被歸類；(6)提取單木位置，每一個(gè)類簇的中心點(diǎn)即為單木位置。

2.3.2 基于CHM的單木位置識別 局部最大值法基于單木樹冠頂部高程大于樹冠邊緣高程的假設(shè)，設(shè)置搜索窗口進(jìn)行局部最大值探測。本研究將搜索窗口的起始尺寸設(shè)置為林分的平均冠幅大小，對窗口尺寸按增加或減少1個(gè)像素進(jìn)行動態(tài)調(diào)整，搜索出樹冠的局部最大值。基于CHM的圖像分割法通過圖像分割實(shí)現(xiàn)單木分割，根據(jù)區(qū)域中心確定單木位置。標(biāo)記控制分水嶺分割法是最常用的分割方法，但在本研究中分割精度不理想，因此選用面向?qū)ο蟮亩喑叨确指睢＠胑Cognition Developer V 8.9軟件的Estimation of Scale Parameter 2 (ESP2)工具包得到合適的分割尺度。本研究中，分割尺度為 66，形狀參數(shù)為0.8，緊致度參數(shù)為0.9，樣地分割效果最好。

2.4 單木樹冠提取

在單木位置識別的基礎(chǔ)上，劃分出同一類點(diǎn)云的投影范圍即為該株華北落葉松的樹冠范圍?？紤]到冠層提取過程中，林分樹冠間隙部分容易產(chǎn)生過分割現(xiàn)象，需要對其結(jié)果進(jìn)行區(qū)域合并[26]，本研究將形狀異質(zhì)性作為類別之間的合并準(zhǔn)則。選擇eCognition軟件對點(diǎn)云初始分類結(jié)果進(jìn)行區(qū)域合并，分割尺度為58，形狀參數(shù)和緊致度參數(shù)分別為0.8和0.8。

2.5 精度評價(jià)

2.5.1 單木位置識別精度評價(jià) 根據(jù)衛(wèi)星定位測量獲取坐標(biāo)的定位精度，考慮樹頂與樹根的位置偏差，將華北落葉松位置誤差預(yù)計(jì)為1 m；以真值坐標(biāo)點(diǎn)為圓心，創(chuàng)建半徑為1 m的緩沖區(qū)。單木位置識別結(jié)果落入此緩沖區(qū)內(nèi)，則認(rèn)為正確。測算正確識別的單木總數(shù)(NC)、實(shí)地調(diào)查單木實(shí)際總數(shù)(NF)、識別出的單木總數(shù)(ND)、重復(fù)識別的單木總數(shù)(NR)等數(shù)據(jù)。計(jì)算精度評價(jià)指標(biāo)，并使用單因素方差分析檢驗(yàn)不同單木位置識別方法在各個(gè)指標(biāo)方面的效果。①總體精度(AO)即正確識別的單木總數(shù)占實(shí)地調(diào)查單木實(shí)際總數(shù)的比值，公式為：AO=NC/NF×100%。②識別精度(AD)即正確識別的單木總數(shù)占識別出的單木總數(shù)的比值，公式為：AD=NC/ND×100%。③漏分誤差(EO)即未識別到的單木總數(shù)占實(shí)地調(diào)查單木實(shí)際總數(shù)的比值，公式為：EO=(NF-NC)/NF×100%。④錯(cuò)分誤差(EC)即被錯(cuò)誤識別為單木的總數(shù)占識別出的單木總數(shù)的比值，公式為：EC=(ND-NC)/ND×100%。⑤重復(fù)性誤差(ER)即對同一株單木重復(fù)識別的總數(shù)占識別出的單木總數(shù)的比值，公式為：ER=NR/ND×100%。

2.5.2 單木樹冠提取精度評價(jià) 采用單株冠幅提取精度 (ω)來進(jìn)行評價(jià)，公式為：ω=100%-∣(C1-C0)/C0∣×100%。其中：C0為實(shí)測冠幅，C1為單木對應(yīng)的提取冠幅。

3 結(jié)果與分析

3.1 不同單木位置識別方法的識別效果

利用4種方法進(jìn)行單木位置識別，空間分布圖(圖2)顯示：4種方法均存在漏分、錯(cuò)分和重復(fù)性識別的現(xiàn)象，結(jié)合表2可知：4種方法的識別精度(AD)由高到低依次為 MSP (89.30%)、LMC (85.60%)、RGP (77.50%)、MSC (70.00%)；MSP 法的單木位置識別效果要好于其他3種方法。4種方法的識別精度均小于總體精度(AO)，說明存在單木錯(cuò)誤識別現(xiàn)象；相比之下，MSP法錯(cuò)分誤差最低(EC為8.0%),，MSC法錯(cuò)分誤差值最高(EC為30.1%)，說明依靠鄰近相似特征識別單木的方法更容易產(chǎn)生將單木識別錯(cuò)誤。相比RGP和MSC法，MSP法的漏分誤差(EO)、錯(cuò)分誤差(EC)和重復(fù)性誤差(ER)均較低，說明高郁閉度針葉林中，按照鄰域分屬識別單木位置會受到樹冠重疊度的影響，并不適用華北落葉松的單木位置識別。盡管LMC法和MSP法的總體精度都高于90%，但后者錯(cuò)分誤差(Ec)較低，說明MSP法更能避免將非樹對象錯(cuò)誤識別為單木而增加識別總數(shù)的現(xiàn)象。

表2 4 種單木位置識別方法的精度統(tǒng)計(jì)Table 2 Accuracy statistics of four methods on individual tree detection

圖2 4 種單木位置識別的空間分布Figure 2 Spatial distribution of the result of individual tree position detection based on four methods

3.2 改進(jìn)均值漂移法的樹冠提取效果

采用均值漂移方法提取華北落葉松單株樹冠(圖3)。以單株華北落葉松實(shí)測冠幅為真值，對利用MSP法提取的單株冠幅進(jìn)行精度分析，結(jié)合表3可知：MSP法總體平均提取精度為90.18%。單株提取精度小于80%的占總體結(jié)果的10.20%，單株提取精度大于80%的占總體結(jié)果的89.80%。單株樹冠分割偏小的結(jié)果占總體的62.07%，要遠(yuǎn)遠(yuǎn)高于提取冠幅偏大的結(jié)果。

表3 MSP 法提取單株冠幅精度Table 3 Precision distribution of extracting individual canopy diameter by MSP

圖3 改進(jìn)均值漂移法樹冠提取結(jié)果 (局部)Figure 3 Crown extraction results of improved mean shift method (local)

4 討論與結(jié)論

4.1 討論

本研究提出的改進(jìn)均值漂移法對高郁閉度華北落葉松林的單木提取總體精度達(dá)91.3%，這與張怡卓等[25]對中低密度人工林的研究結(jié)果(總體精度達(dá)89.0%)相似，遠(yuǎn)高于JING等[27]對加拿大高郁閉林分的總體精度(69.0%)。相比于其他3種提取方法，MSP法有2個(gè)優(yōu)勢。一是在林木樹冠重疊的情況下，MSP法直接利用了點(diǎn)云空間特征，避免了生成CHM過程中產(chǎn)生的柵格誤差；二是該方法利用統(tǒng)計(jì)迭代的方式確定聚類點(diǎn)位，根據(jù)訪問頻率分類，比直接比較鄰域特征而進(jìn)行分類穩(wěn)健性更好。

本研究發(fā)現(xiàn)：MSP法提取單木冠幅的平均提取精度為90.18%，比同等郁閉度林分的樹冠提取精度要高出很多[28]，可滿足森林調(diào)查需求。約10%的單株冠幅提取精度低于80%，這些點(diǎn)主要位于林分下層。由于上層遮擋，下層點(diǎn)云密度遠(yuǎn)低于上層，而點(diǎn)云聚類算法和點(diǎn)云密度息息相關(guān)，這也是本研究發(fā)現(xiàn)冠幅提取值偏小的原因。具體來說，由于上層冠層存在重疊，人工目視解譯可以考慮到冠層重疊的影響并準(zhǔn)確提取冠幅，而在三維點(diǎn)云圖上，樹冠界限分類時(shí)存在距離閾值類似的情況，影響提取結(jié)果；實(shí)測中，研究者往往通過東—西、南—北方向來測量冠幅，但是這只是一種理想狀態(tài)下的測量，在現(xiàn)實(shí)中，華北落葉松林木樹冠是不規(guī)則的，測量因此會存在誤差；同時(shí)研究區(qū)域的坡度、地形等的限制，也會對人工測量造成影響。

4.2 結(jié)論

本研究基于點(diǎn)云單一空間特征，引入高斯核函數(shù)簡化并改進(jìn)了均值漂移量的計(jì)算過程，同時(shí)所提方法不需要先驗(yàn)知識，算法的穩(wěn)健性高，解決了高郁閉度的華北落葉松林受復(fù)雜冠層影響的問題，提高了在該林分內(nèi)機(jī)載LiDAR單木識別的精度，為高郁閉度華北落葉松林單木識別提供了新的方法。