徐志揚，劉浩棟，陳永富，陳巧＊，李華玉,3，王娟,3

(1.中國林業(yè)科學(xué)研究院資源信息研究所，國家林業(yè)和草原局林業(yè)遙感與信息技術(shù)重點實驗室，北京 100091；

2.國家林業(yè)和草原局華東調(diào)查規(guī)劃設(shè)計院，浙江 杭州 310019；

3.西南林業(yè)大學(xué)林學(xué)院，云南 昆明 650224)

杉木（Cunninghamia lanceolata(Lamb.) Hook.）是我國南方重要用材樹種，分布廣泛，杉木上冠在杉木全冠占主要地位，決定了杉木的光合作用、生長狀況以及冠形特征，因此對杉木樹冠上部外輪廓模擬具有重要意義，準確描繪樹冠結(jié)構(gòu)可以為樹種鑒別提供參考。樹冠結(jié)構(gòu)因子包括冠幅、冠長、冠基高、樹冠形狀等，其中樹冠形狀可以用樹冠輪廓模型表達[1]。國內(nèi)外進行了很多樹冠輪廓模型的研究，有些學(xué)者使用規(guī)則幾何體表達[2-3]，可是很難用簡單幾何體對自然界各種形態(tài)各異的樹冠進行特定描述，因此學(xué)者們也采用相對復(fù)雜的數(shù)學(xué)模型模擬樹冠輪廓。郭艷榮等[4]、吳丹子等[5]采用多項式和指數(shù)函數(shù)模擬杉木樹冠外輪廓。由于樹冠上部和下部往往具有不同的幾何形狀，一些學(xué)者認為采用分段形式能夠更為精確的描繪樹冠外輪廓。Dong等[6]分上冠、下冠、全冠對福建順昌杉木采用9個模型擬合并選擇最佳模型估計樹冠形狀。高慧淋等[7-9]利用東北林區(qū)人工紅松（Pinus koraiensisSieb.et Zucc.）、人工長白落葉松（Larix olgensisHenry）、人工樟子松（Pinus sylvestrisvar.mongolicaLitv.）枝解析數(shù)據(jù)，采用分段拋物線方程、分段冪函數(shù)方程、分段單分子式方程、修正Kozak 方程和修正Weibull 方程對樹冠輪廓進行模擬，構(gòu)建3個樹種的樹冠輪廓非線性混合效應(yīng)模型。MARSHALL等[2]、盧康寧等[10]、王小明等[11]均將樹冠從最大半徑位置處分為上部和下部兩部分，分別采用線性或非線性形式單獨建立方程模擬，取得不錯的效果。

樹冠輪廓建模變量的測量一般采取樹木伐倒再枝解析方式解決，精確度較高，但是效率不夠高。主動遙感作為傳統(tǒng)遙感的重要補充，為林業(yè)調(diào)查以及測量提供了新的途徑。激光雷達（LiDAR）是一種主動遙感技術(shù)，能夠穿透樹冠獲取其三維結(jié)構(gòu)信息，在林分及單木水平上都得到了廣泛應(yīng)用[12-13]。無人機（UAV）作為一種新興的低空遙感平臺，能夠靈活、高效地獲取遙感數(shù)據(jù)，無人機和激光雷達的系統(tǒng)集成已成為森林精細化調(diào)查的有力支撐[1,14]。利用無人機激光雷達能夠快速獲取杉木的樹冠垂直結(jié)構(gòu)參數(shù)，但是，無人機激光雷達在杉木樹冠輪廓模擬以及可視化方面的研究卻鮮有報道。

本研究采用半自動化的方式，利用無人機激光雷達數(shù)據(jù)提取杉木單木點云進行杉木樹冠上部外輪廓建模，以相對著枝深度為自變量、枝長為因變量，分別采用多項式、冪函數(shù)、指數(shù)函數(shù)建立杉木樹冠上部外輪廓模型，并選擇最優(yōu)模型進行可視化，旨在為研究樹種識別和單木預(yù)測提供參考。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)位于江西省分宜縣境內(nèi)的中國林業(yè)科學(xué)研究院亞熱帶林業(yè)實驗中心年珠實驗林場（27°30′～27°45′ N，114°30′～114°45′ E），屬低山丘陵，海拔220～1 092 m，土壤為黃紅壤，年均氣溫16.8℃，日最高氣溫39.9℃，最低氣溫?8.3℃，年降水量1 950.9 mm，集中在3—6月，無霜期252 d。主要森林類型為杉木林（Cunninghamia lanceolata(Lamb.) Hook.）、毛竹林（Phyllostachys edulis(Carr.) H.de Lehaie）、闊葉林等。闊葉林中主要樹種有樟樹（Cinnamomum camphora(L.)Presl.）、大葉錐栗（Castanea henryi(Skan) Rehd.et Wils.）、鉤錐（Castanopsis tibetanaHance）、甜櫧（Castanopsis eyrei(Champ.ex Benth.)Tutch.）、苦櫧（Castanopsis sclerophylla(Lindl.et Paxton) Schottky）、木荷（Schima superbaGardn.et Champ.）、刨花潤楠（Machilus pauhoiKanehira）等。

1.2 數(shù)據(jù)來源

無人機搭載RIEGL VUX-1LR 激光雷達傳感器，通過近紅外（1 550 nm）激光束和旋轉(zhuǎn)鏡330°視場角快速掃描實現(xiàn)數(shù)據(jù)的高速獲取，精度為15 mm，數(shù)據(jù)采集于2019年植被生長旺季，采用仿地飛行模式，以地形表面為基準設(shè)置相對高度，依次在40、70、100、130、160 m 共5個不同高度掃描獲取激光雷達數(shù)據(jù)。此外，無人機搭載Sony ILCE-6000 微單相機采集可見光數(shù)據(jù)以供后續(xù)生成高分辨率正射影像（DOM）、搭載RedEdge M 快照式多光譜相機采集5 波段多光譜數(shù)據(jù)以供后續(xù)生成多光譜影像。使用160 m 航高激光雷達數(shù)據(jù)，正射影像作為激光雷達數(shù)據(jù)的輔助參考材料。

地面樣地數(shù)據(jù)的采集也是在2019年植被生長旺季完成，依據(jù)研究區(qū)范圍內(nèi)森林資源規(guī)劃設(shè)計調(diào)查（“二類調(diào)查”）數(shù)據(jù)中的樹種信息，在無人機有效飛行范圍內(nèi)設(shè)置面積為0.04 hm2的圓形樣地，分別記錄樣地中心點經(jīng)緯度與3株定位樹經(jīng)緯度，進行每木檢尺，記載5 cm 以上樹木的每木相對位置、樹種、胸徑、樹高、枝下高、林木分級、枯死情況、東西與南北向冠幅等。選擇樹種組成為杉木純林的2個樣地（46 號、47 號樣地）進行樣地的單株杉木提取，樣地所在林分為杉木近成熟林。

2 研究方法

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

通過數(shù)字綠土公司LiDAR360 軟件進行激光雷達點云數(shù)據(jù)預(yù)處理。首先去除點云中的噪聲點，它一般為孤立點，包括明顯高于地表物和低于地表面的激光點，根據(jù)絕對高程或閾值去除較為明顯的異常點；其次進行點云分類，將點云分為地面點和非地面點，本研究非地面點為森林的激光雷達反射脈沖點；然后，利用Kriging 插值法對地面點進行插值，得到0.5 m 分辨率的數(shù)字高程模型（DEM），進而對點云進行歸一化，將去噪分類后的點云高程值減去對應(yīng)DEM 像元值得到點云相對地面點的絕對高度。為減少CHM 孔洞影響，參照相關(guān)研究的孔洞填充算法[15-16]，利用LASTools 開源工具對歸一化后的點云生成無孔洞的冠層高度模型（Pit-Free CHM）。

對于地面樣地數(shù)據(jù)，需要將屬性表示的每木數(shù)據(jù)轉(zhuǎn)換為矢量數(shù)據(jù)。首先在ArcGIS 中添加經(jīng)緯度記錄的中心點和定位樹位置數(shù)據(jù)，定義地理坐標(biāo)系為GCS_WGS_1984，生成經(jīng)緯度表示的SHP 矢量圖層，然后進行投影變換，定義投影坐標(biāo)系為WGS_1984_UTM_Zone_50N，轉(zhuǎn)換為橫縱坐標(biāo)表示的SHP 矢量數(shù)據(jù)，以中心點橫縱坐標(biāo)值為原點對樣地內(nèi)每木相對位置進行三角函數(shù)解算以得到每木的橫縱坐標(biāo)，最后以橫縱坐標(biāo)為XY值將每木數(shù)據(jù)添加到ArcGIS 中生成每木矢量位置數(shù)據(jù)，采用投影變換后的定位樹位置對每木位置進行檢查控制，從而得到樣地內(nèi)單木矢量數(shù)據(jù)。

2.2 杉木參數(shù)提取

杉木參數(shù)提取中，在LiDAR360 軟件下完成單木分割，在ArcGIS 下完成單木匹配與補充建模樣本，單木點云純化與參數(shù)（樹高、冠幅、上部冠長）提取等均采用Python3.6 編程自動化完成。

2.2.1 單木分割與純化 采用基于冠層高度模型（CHM）種子點的點云分割算法，先利用局部最大值濾波器對生成的無孔洞的冠層高度模型進行單木樹頂探測生成種子點，根據(jù)種子點完成單木點云分割。

采用選定的2個杉木純林樣地對分割后的單木進行匹配檢驗。借鑒Reitberger等[17]的匹配原則來匹配分割后的單株樹木與實測樣木，在ArcGIS 中疊加CHM、檢測到的樣木、實測樣木，經(jīng)多次對比試驗，探測到的樹頂與參考樹頂距離若小于樣地中所有樣木平均距離的60%、而且樹高差異小于樣地中最大的樹高20% 的單木為1:1 匹配的樣木，若一個參考的樹頂對應(yīng)了多個探測的樹頂，則將探測的最小距離的樹頂作為1∶1 匹配的樣木，然后從匹配的樣木中剔除枯死木、下層被壓木、非杉木，剩余43株杉木作為單株模擬樣本。

為彌補單株杉木樣本量不足的缺陷，以DOM和CHM 為底圖，結(jié)合“二類調(diào)查”矢量數(shù)據(jù)優(yōu)勢樹種屬性和先驗知識，在樣地周邊以及大片杉木純林中從單株探測結(jié)果矢量中目視選擇補充杉木單株樣本到150個以上。

對于本研究的密集林分下分割后的單株點云，仍然存在著“欠分割”樣本，這將導(dǎo)致提取的單株樹冠特征因子異常與無效，因此需要對單株杉木點云再進行純化。為有效提取樹冠的特征因子，引入Kim等[18]單株分割點云純化思路并改進優(yōu)化，首先，以單株杉木樹頂為中心，在水平方向上，按照90°間隔將單株點云分為4個分區(qū)，各分區(qū)的點云再按照0.5 m 間隔從樹木中心向外逐層分段緩沖，統(tǒng)計在垂直方向上各分區(qū)、各段點云的平均高度；其次，以實地調(diào)查的單株杉木樹高與冠幅數(shù)據(jù)建立樹高-冠幅模型（CW=HT0.3446?1.0645，R2為0.458 3，其中CW為冠幅，HT為樹高），利用該模型結(jié)合單木CHM 最大值得到一個冠幅預(yù)測值，再將該值乘以1.25 倍系數(shù)，以此對各分區(qū)的分段緩沖進行冠幅約束。正確分割的單株杉木點云各分區(qū)中逐層緩沖的高度值從樹木中心點向樹冠邊緣逐漸減小，而當(dāng)單株點云“欠分割”時，就會發(fā)生平均高度減小到一定水平（拐點）然后增大的現(xiàn)象，因此，“欠分割”單株點云的該拐點處應(yīng)當(dāng)是可識別的樹冠邊緣，故舍棄拐點處以外的點云，保留各段內(nèi)樹木中心到拐點分區(qū)的點云，形成純化后的目標(biāo)單株點云樣本。單木點云純化原理如圖1 所示（XY頂視圖）。

圖1 單木點云純化示意圖Fig.1 Schematic diagram of individual tree point clouds purification

2.2.2 樹高提取 單株杉木點云純化后完整的單木點云為樹冠點與樹干點的非地面點集合，單株點云的最高點即為樹冠頂點，樹冠頂點高度就是單木樹高。為減少噪聲點對樹冠頂點的干擾，使用單株杉木CHM 最大值對單株點云的最高點進行檢查，二者差異絕對值大于固定閾值則對單株點云輸出二維X-Z圖像進行目視檢查，確定是否存在噪聲點，本研究此閾值設(shè)置為0.5 m。

2.2.3 冠幅提取 參考Reitberger等[17]的分層處置方法，模擬冠幅的定義提取單木冠幅。在垂直方向上，對純化后的完整單株點云，自下而上以0.5 m為間隔分層，統(tǒng)計每層點云數(shù)量與占整株點云數(shù)的百分比；同時，各分層在水平方向上，以樹冠頂點垂直投影為中心，按照5°間隔對點云進行分區(qū)，共得到72個分區(qū)，求取分區(qū)中點云距離樹木中心的最大值作為分區(qū)的樹冠半徑，取各分區(qū)樹冠半徑均值的2 倍為該層樹冠平均直徑。經(jīng)多次實驗比較，取自下而上點云累計數(shù)量占整株點云數(shù)2%以上、高度5 m 以上的最大樹冠直徑為冠幅。如圖2 所示，(a)對單木點云進行分層，(b)為各層點云數(shù)量占單木總點云數(shù)量比例，(c)為對各分層點云進行分區(qū)求取樹冠平均直徑。

圖2 冠幅提取示意圖Fig.2 Schematic diagram of crown width extraction

2.2.4 上部冠長提取 從樹冠頂點到最大樹冠位置的部分即為樹冠上部，單木樹高減去最大樹冠半徑所在高度就可得到上部冠長（CLU）。

2.3 樹冠上部輪廓模擬

構(gòu)建樹冠上部外輪廓模型需要獲取單株樹冠的外部輪廓點。本研究參考相關(guān)研究采用寬度分位數(shù)表示樹冠外輪廓的方法[7]，對純化后的單株樹冠點云自最大樹冠位置處向上按照0.5 m 間隔分層提取，分別統(tǒng)計各層點云寬度百分位數(shù)85%、90%、95%、98%當(dāng)作相應(yīng)層外輪廓點，目視對比樹冠上部外輪廓模型精度，最終選擇95%分位數(shù)建立樹冠外輪廓模型。

將處理后的單株樣木按照4:1 隨機分成建模樣本與檢驗樣本，并且對建模和檢驗的樹冠外輪廓點樣本分別按照相對著枝深度分層，使用3 倍標(biāo)準差法剔除外部樹冠半徑異常的樹冠上部外輪廓點。以相對著枝深度（RDINC）和外部樹冠半徑（OR）建立模型，相對著枝深度為樹冠輪廓點到樹冠頂點在垂直方向上的相對距離，等于著枝深度與上部冠長之比，外部樹冠半徑等于輪廓點到樹冠中心的水平距離。建模變量模型參數(shù)見圖3 所示。

圖3 單木樹冠上部外輪廓模型參數(shù)示意圖Fig.3 Schematic diagram of individual tree upper crown profile model’s parameters

構(gòu)建模型時，選取樹冠輪廓模擬較為常用的二次拋物線、冪函數(shù)和對數(shù)函數(shù)作為基礎(chǔ)模型[4-6,19]，方程形式如下：

采用決定系數(shù)（R2）和均方根誤差（RMSE）評價擬合優(yōu)度，使用平均偏差（ME）和平均絕對偏差（MAE）檢驗?zāi)Ｐ汀?/p>

樹冠上部輪廓模擬時，外輪廓點參數(shù)提取、單木樣本隨機分配和異常輪廓點剔除、模型建立均采用Python3.6 編程完成。

3 結(jié)果與分析

3.1 單木探測及空間格局

兩個杉木林圓形樣地中，46 號樣地31株樹木中無人機激光雷達探測到23株全部為杉木，47 號樣地23株樹木中探測出20株全部為杉木，兩個樣地探測率分別為74.19% 和86.96%，綜合探測率79.63%，如圖4 所示，可以看出樣地中杉木水平分布并不均勻。圖4 中的單株樹冠邊界顯示的是提取后的冠幅。

圖4 無人機激光雷達單木探測結(jié)果Fig.4 Individual tree detection of UAV-LiDAR

3.2 單木樹冠特征因子提取結(jié)果比較

圖5 所示的是43株檢測到的杉木樣本提取樹高與實測樹高的散點圖與回歸關(guān)系，回歸相關(guān)系數(shù)R2為0.890 5，單木樹高提取精度均值達到95.25%，表明無人機激光雷達能夠獲得與地面實測相關(guān)性較好且精度較高的預(yù)測樹高值。

圖5 提取樹高與實測樹高的回歸關(guān)系Fig.5 Regression equation of inversed tree height and measured tree height

圖6 所示的是43株檢測到的杉木樣本提取冠幅與實測冠幅的散點圖與回歸關(guān)系，回歸相關(guān)系數(shù)R2為0.845 6，單木冠幅提取精度均值達到94.32%，表明本研究單木冠幅提取算法較好的反映了樹冠真實大小和分布。

圖6 提取冠幅與實測冠幅的回歸關(guān)系Fig.6 Regression equation of the inversed crown and measured crown

3.3 樹冠上部外輪廓擬合模型

采用95% 分位數(shù)方法提取出141株建模杉木樣本的2 538個樹冠上部外輪廓點，再按照3 倍標(biāo)準差法剔除異常外輪廓點得到2 152個樹冠上部外輪廓點，據(jù)此擬合出了3個杉木樹冠上部外輪廓模型，擬合結(jié)果如圖7 所示。

圖7 無人機激光雷達擬合杉木樹冠上部輪廓Fig.7 The fitted upper canopy profile of C.lanceolata with UAV-LiDAR

由表1 可看出，參數(shù)估計值標(biāo)準誤均較低，可見參數(shù)估計值較為穩(wěn)定。建模樣本的R2均達到0.8 以上，二次拋物線和對數(shù)函數(shù)擬合R2較近，冪函數(shù)擬合R2明顯高于二者，達到0.817 0，且冪函數(shù)的均方根誤差（RMSE）均低于二者，從驗證樣本的平均偏差（ME）和平均絕對偏差（MAE）來看，冪函數(shù)亦明顯優(yōu)于二次拋物線和對數(shù)函數(shù)，因此本研究中冪函數(shù)擬合效果最優(yōu)。表1 中模型參數(shù)估計值、建模評價與驗證指標(biāo)編程得到的結(jié)果與SPSS 19.0 所得結(jié)果一致，估計值標(biāo)準誤從SPSS19.0 得到。

表1 無人機激光雷達杉木樹冠上部輪廓擬合與模型檢驗Table 1 Upper canopy profile model fitting and validating of C.lanceolata with UAV-LiDAR

3.4 樹冠上部三維可視化模擬

在Python3.6 下編程讀入已提取的單株點云的樹高、最大樹冠半徑和上部冠長，將單木參數(shù)代入冪函數(shù)外輪廓模型后可以生成一組高度和樹冠半徑數(shù)據(jù)，再以杉木中心為原點將樹冠半徑轉(zhuǎn)入X-Y平面中進行三角函數(shù)解算，將外輪廓二維數(shù)據(jù)轉(zhuǎn)換為三維模擬數(shù)據(jù)，加載樣地內(nèi)所有樣木的三維模擬數(shù)據(jù)，調(diào)用matplotlib 庫對兩個建模杉木林樣地進行樣地尺度的三維重建，為減少地理橫、縱坐標(biāo)大量的數(shù)字冗余顯示，在三維坐標(biāo)系中保持Z軸為樹木高度而在水平面上對樣地中樣木地理位置的橫、縱坐標(biāo)進行整體平移，在樣地范圍內(nèi)將地理坐標(biāo)原點移動到樣地橫、縱坐標(biāo)最小值位置，并用X-Y軸分別表示相對地理橫、縱坐標(biāo)，結(jié)果如圖8 所示。

圖8 2個樣地樹冠上部三維模型重建Fig.8 Three dimensions model reconstruction of 2 plot’s tree upper canopy

4 討論

在單木提取方面，本研究2個圓形樣地54株樹木檢測出1∶1 匹配的樣本43株（全為杉木），綜合探測率為79.63%。與同樣采用無人機激光雷達技術(shù)進行單木探測結(jié)果比較，Wallace等[20]分別采用5種提取算法對4年生桉樹林進行單木提取，提取精度均達到90%以上，最高為98.00%，與之相比本研究探測率明顯偏低，這除了與點云密度有關(guān)外，還與提取對象的林分特征有關(guān)，桉樹林中樹木在水平方向上更為均勻且垂直方向上更為整齊，桉樹林的這種林分特征能夠非常明顯的提高單木提取精度，而本研究的杉木林單株樹木分布不均勻，甚至部分樹木群團聚集在一起，46 號樣地因此現(xiàn)象有5株樹木未被檢測出，此外，杉木近成熟林中部分杉木樹冠衰老、枝葉凋落造成點云質(zhì)量不佳也影響了單木檢測結(jié)果；在國內(nèi)，文獻[1]基于CHM和標(biāo)記控制分水嶺法對東北落葉松林進行單木提取，兩個樣地323株共計檢測出1∶1 匹配的樹木235株，綜合精度為72.76%，本研究的杉木林較落葉松林更為復(fù)雜但結(jié)果仍然略高于其提取精度，原因主要是本研究采用基于CHM種子點的單木點云分割方法，其次應(yīng)與樣地大小（樣本量）有關(guān)，樣本較少時得出的探測率意義并不強，這也是下一步工作中需要注意的地方。

本研究對象為密集林分的單株對象，單木點云分割后仍然存在“欠分割”現(xiàn)象，為準確提取單株樹冠的建模特征因子，研究中改進Kim[18]的單株點云純化思想，利用理想狀態(tài)的單木樹冠輪廓點高度值呈現(xiàn)出從樹木中心點向外逐漸減小的規(guī)律，并結(jié)合樹高-冠幅模型進行樹冠點云約束，對單木點云進行篩選純化。因杉木為針葉樹，冠形規(guī)律性較強，該算法效果較為明顯，但是，對于其他樹種則不一定適用。

本研究中3個模型的建模擬合R2均達到了0.80 以上，較文獻[1]基于無人機激光雷達的0.752要高，這與對樹冠上部外輪廓點樣本按照相對著枝深度分層進行3 倍標(biāo)準差法剔除外部樹冠半徑異常點所起作用有關(guān)；與郭艷榮等[4]基于實測值建模的杉木中齡林擬合R2結(jié)果（0.80～0.85）較為接近，略高于吳丹子等[5]基于杉木實測值建模擬合的R2（0.76～0.80），表明在特定樹種和林分條件下，采用無人機激光雷達模擬樹冠輪廓，理論上可以達到實測值的擬合效果。

本研究僅對樹冠上部進行外輪廓模擬，一是為基于無人機多源遙感數(shù)據(jù)樹種分類識別提供樹木冠形方面的參考，二是亞熱帶密集林分條件下獲取的樹冠上部點云尚能夠反映樹木的冠形狀態(tài)，而樹冠下部點云因林分條件影響較為稀疏或者對樹木冠形描繪無任何作用。

大多數(shù)樹冠輪廓模型與可視化以實測的樹木枝干解析數(shù)據(jù)為樣本材料，往往通過選取平均木伐倒進行枝條測量，精度較高但是效率有限，且人力物力耗費較大，而無人機激光雷達能夠提高工作效率，本研究也表明了利用無人機激光雷達表達杉木樹冠上部外輪廓冠型結(jié)構(gòu)的可行性。

本研究杉木人工同齡林為近成熟林，而同一樹種在不同齡組的樹冠形態(tài)并不完全一致，今后可以分別齡組開展研究。

5 結(jié)論

本研究中單株樹木點云外輪廓參數(shù)提取是關(guān)鍵，采用單木點云分割與純化相結(jié)合得到純化后的單木點云數(shù)據(jù)，進而提取外輪廓參數(shù)，程序?qū)崿F(xiàn)便捷，但適用對象有限，因此研究在高密度林分條件下單木點云有效提取純化方法對于客觀描繪樹冠三維形狀具有重要價值?；跓o人機激光雷達建立的3個杉木樹冠上部模型均取得較好的擬合效果，而冪函數(shù)模型表現(xiàn)更優(yōu)，模擬結(jié)果反映了杉木的樹冠上部形態(tài)，它為杉木的分類識別提供了參考依據(jù)。