龍 飛,岳彩榮,金 京,李春干,羅洪斌,徐婉婷

(1.西南林業(yè)大學(xué) 林學(xué)院,昆明 650224；2.廣西大學(xué)林學(xué)院,南寧 530004)

森林生態(tài)系統(tǒng)在水平和垂直維度上都體現(xiàn)出較強(qiáng)的異質(zhì)性,不僅在維持全球生物多樣性上有重要作用,同時(shí)還參與組成全球生物碳循環(huán)系統(tǒng)[1-2]。森林生物量則是作為評(píng)價(jià)森林生產(chǎn)力的重要評(píng)價(jià)指標(biāo),故實(shí)現(xiàn)森林生物量快速高效地監(jiān)測(cè),對(duì)全球生物碳循環(huán)有著重要的現(xiàn)實(shí)意義[3-4]。傳統(tǒng)估測(cè)森林生物量的主要手段是通過(guò)大規(guī)模的森林資源調(diào)查來(lái)實(shí)現(xiàn),但此類調(diào)查方式弊端明顯,耗財(cái)耗力的同時(shí),甚至需要破壞原有的森林生態(tài)系統(tǒng)來(lái)完成取樣,不僅人為產(chǎn)生的誤差相對(duì)較大,還會(huì)因?yàn)榄h(huán)境條件限制部分地方無(wú)法實(shí)地測(cè)量[5-6]。機(jī)載激光雷達(dá)(Airborne Light Laser Detection and Ranging,LiDAR)的出現(xiàn)很好地彌補(bǔ)和解決了這一問(wèn)題,給未來(lái)的森林資源調(diào)查提供一種新的選擇。

機(jī)載LiDAR相較于光學(xué)遙感具有較強(qiáng)的穿透性和主動(dòng)性,精度上也相對(duì)較高,可以獲取森林的三維結(jié)構(gòu)信息,此類信息在林業(yè)方面得到廣泛應(yīng)用[7-8]。胸高斷面積(Basal Area,BA)是指樹(shù)木距地面1.3m處樹(shù)干的橫切面面積,一般是通過(guò)測(cè)量樹(shù)干橫切面的直徑得到,是預(yù)測(cè)林分材積和森林生長(zhǎng)的良好指標(biāo)[9-11],也是林業(yè)工作者在進(jìn)行森林資源調(diào)查時(shí)必不可少的林木參數(shù)之一[12]。前期,有關(guān)機(jī)載LiDAR數(shù)據(jù)在森林特征變量提取方面的研究主要是圍繞估測(cè)樹(shù)高、冠幅、林分生物量、材積和碳儲(chǔ)量等方面[13-15],用機(jī)載LiDAR數(shù)據(jù)反演林分胸高斷面積的研究相對(duì)較少。對(duì)林木胸高斷面積的估測(cè),主要集中在運(yùn)用不同類型的激光雷達(dá)進(jìn)行多元線性回歸模型[16]、點(diǎn)云密度[17]、k近鄰算法[18]、均方平均冠層高度[19]等方面,但采用機(jī)載LiDAR數(shù)據(jù)考慮濾波以及隨機(jī)森林(Random Forest,RF)和迭代的決策樹(shù)算法(Gradient Boost Regression Tree,GBRT)估測(cè)林木胸高斷面積的研究較少。

本研究以廣西壯族自治區(qū)的高峰林場(chǎng)為實(shí)驗(yàn)區(qū),以VUX1LR激光雷達(dá)系統(tǒng)獲取的LiDAR數(shù)據(jù)和實(shí)測(cè)的樣地?cái)?shù)據(jù)為基礎(chǔ),比較漸進(jìn)不規(guī)則三角網(wǎng)(Progressive Triangulated Irregular Network,PTIN)、布料模擬濾波算法(Cloth Simulation filter,CSF)、漸進(jìn)形態(tài)學(xué)濾波算法(Progressive Morphological Filter,PMF)和基于插值的濾波方法(Interpolation-Based Filtering,IBF)對(duì)點(diǎn)云濾波效果的影響,同時(shí)采用隨機(jī)森林和迭代的決策樹(shù)進(jìn)行建模,選取模型精度高、適應(yīng)性好的模型完成森林胸高斷面積制圖。通過(guò)本文的研究,旨在探索一種快速高效獲取林木胸高斷面積的方法,以期為今后進(jìn)一步研究提供參考。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于廣西壯族自治區(qū)國(guó)營(yíng)高峰林場(chǎng),地理位置為22°49′～23°15′N,108°08′～108°53′E,是目前廣西最大的國(guó)有林場(chǎng),整體是一個(gè)近矩形區(qū)域(長(zhǎng)11.2km,寬4.2km,面積約為47.04km2)(圖1)。實(shí)驗(yàn)樣地主要以高峰林場(chǎng)界碑、東升、延河3個(gè)分場(chǎng)為主,此外,還包括武鳴區(qū)、興寧區(qū)、西鄉(xiāng)塘區(qū)的少量林地[20]。高峰林場(chǎng):屬低山地貌,高程范圍在88.76～462.38m,坡度多數(shù)在20～35°,海拔由西南向東北逐漸升高；屬南亞熱帶季風(fēng)氣候區(qū),該氣候條件下植被多為南亞熱帶常綠闊葉林；森林覆蓋率較好,保持在90%以上；天然林占比較低,約為5%,其余為人工林；按照樹(shù)木生長(zhǎng)周期可將其劃分為短周期經(jīng)營(yíng)的桉樹(shù)人工林和長(zhǎng)周期經(jīng)營(yíng)杉木人工林。由2016年森林資源調(diào)查數(shù)據(jù)可知[20],研究區(qū)主要樹(shù)種為尾葉桉(EucalyptusurohpyllaS.T.Blake)、巨尾桉(Eucalyptusgran-dis×E.urohpylla)、馬尾松(PinusmassonianaLamb)、濕地松(PinuselliottiiEngelmann)、杉木(Cunninghamialanceolata(Lamb.)Hook)、八角(llliciumverumHook.f.)等。樣地情況如圖1所示。

圖1 研究區(qū)樣地分布情況示意圖Fig.1 Schematic diagram of the sample plots distribution in the study area

1.2 樣地?cái)?shù)據(jù)獲取

樣地?cái)?shù)據(jù)獲取于2016年5月至12月?？紤]到森林類型、林齡和樹(shù)種組成等因素對(duì)實(shí)驗(yàn)準(zhǔn)確性的影響,采用手持雙頻差分GPS接收器記錄樣地四個(gè)角點(diǎn)的地理坐標(biāo),用激光測(cè)距儀和森林羅盤儀設(shè)置30m×30m的方形樣地共計(jì)105塊。對(duì)樣地內(nèi)胸徑(DBH)≥5cm的樹(shù)木進(jìn)行每木檢尺并記錄樹(shù)種、胸徑、樹(shù)高、枝下高、冠幅等實(shí)驗(yàn)數(shù)據(jù),優(yōu)勢(shì)樹(shù)種主要有杉木、馬尾松、巨尾桉和一般闊葉樹(shù)。

1.3 機(jī)載LiDAR數(shù)據(jù)的獲取

2016年9月,采用R4直升機(jī)搭載集激光測(cè)距、全球定位系統(tǒng)和慣性導(dǎo)航系統(tǒng)為一體的VUX1LR激光雷達(dá)系統(tǒng)來(lái)獲取樣地點(diǎn)云數(shù)據(jù),其中,激光器波長(zhǎng)為1 550nm,激光發(fā)射角0.5mrad,發(fā)射頻率820KHz,點(diǎn)云密度約為2.9point/m2。采用CCD相機(jī)來(lái)獲取研究區(qū)航空影像,其中,像素為3 000萬(wàn),空間分辨率為0.2m。R4直升機(jī)獲取樣地時(shí)的實(shí)際飛行高度約為1 000m,覆蓋范圍為23°0′28″～22°56′3″N,108°20′9″～108°27′33″E,飛行3次共計(jì)15h,覆蓋面積約為55km2。

2 研究方法

2.1 樣地?cái)?shù)據(jù)處理

按照樣地植物葉的形態(tài)特征將森林類型劃分為針葉林和闊葉林,并計(jì)算林分平均高、樣地蓄積和林分胸高斷面積。通過(guò)樣地胸高斷面積和平均樹(shù)高,查閱形高表計(jì)算樣地的蓄積量。樣地?cái)?shù)據(jù)統(tǒng)計(jì)情況如表1所示。

表1 樣地?cái)?shù)據(jù)統(tǒng)計(jì)表Tab.1 Plot data statistics table

2.2 機(jī)載LiDAR數(shù)據(jù)的處理

通過(guò)去噪、濾波和分類對(duì)機(jī)載LiDAR數(shù)據(jù)進(jìn)行預(yù)處理。為提高地面點(diǎn)分類效果,濾波時(shí)就漸進(jìn)不規(guī)則三角網(wǎng)(PTIN)、漸進(jìn)形態(tài)學(xué)濾波算法(PMF)、布料模擬濾波算法(CSF)和基于插值的濾波方法(IBF)在樣地中的濾波效果進(jìn)行對(duì)比,同時(shí)結(jié)合實(shí)驗(yàn)區(qū)環(huán)境選擇可適性較高的方法完成點(diǎn)云濾波。點(diǎn)云濾波方法如下:

1)漸進(jìn)不規(guī)則三角網(wǎng)算法(PTIN)在森林點(diǎn)云數(shù)據(jù)濾波研究中常被使用?；舅悸?首先,根據(jù)一定策略選擇初始地面點(diǎn)完成初始地形的建立；其次,采用相應(yīng)準(zhǔn)則不斷迭代增加地面點(diǎn)；最后,將相關(guān)數(shù)據(jù)疊加至所有的點(diǎn)完成歸類。經(jīng)典的三角網(wǎng)加密方法在坡度陡峭的區(qū)域?yàn)V波效果不佳,故此本研究使用郭慶華團(tuán)隊(duì)提出改進(jìn)的漸進(jìn)三角網(wǎng)方法進(jìn)行濾波[21]。

2)漸進(jìn)形態(tài)學(xué)濾波算法(PMF)的理論基礎(chǔ)是形態(tài)學(xué)。該方法主要是利用具有特定形狀的結(jié)構(gòu)元素提取圖像中的拓?fù)浣Y(jié)構(gòu)信息,從而用于圖像分類,其基本操作為膨脹和腐蝕[22]。

3)布料模擬濾波算法(CSF)是基于“l(fā)idR”包實(shí)現(xiàn),主要將激光LiDAR數(shù)據(jù)進(jìn)行反轉(zhuǎn),并用硬布覆蓋反轉(zhuǎn)表面?；玖鞒?首先,通過(guò)分析分布節(jié)點(diǎn)和相應(yīng)激光LiDAR數(shù)據(jù)點(diǎn)之間的相互作用,從而確定分布節(jié)點(diǎn)來(lái)生成近似地面點(diǎn)；其次,通過(guò)對(duì)比原始激光雷達(dá)點(diǎn)與生成的表面之間的關(guān)系來(lái)提取地面點(diǎn)；最后,布料模擬本身是由質(zhì)量和相互連接的粒子組成的網(wǎng)格,這些粒子共同確定布料三維位置和形狀[23]。

4)基于插值的濾波方法(IBF)是通過(guò)ENVI的BCAL-tool插件來(lái)實(shí)現(xiàn)。一般流程:首先,依靠一定的規(guī)則形成一個(gè)初始地表面；其次,采用特定規(guī)則對(duì)未分類的數(shù)據(jù)進(jìn)行判定；最后,不斷迭代加入新的地面點(diǎn),迭代到符合停止加密的規(guī)則為止。

點(diǎn)云數(shù)據(jù)濾波效果評(píng)價(jià)主要采用Ⅰ類誤差(地面點(diǎn)誤分為非地面點(diǎn))、Ⅱ類誤差(非地面點(diǎn)誤分為地面點(diǎn))和kp系數(shù)為主要指標(biāo)(表2)[2],利用LiDAR360軟件進(jìn)行人工檢查和修改。受到森林結(jié)構(gòu)特征影響,在復(fù)雜地形條件或高郁閉度條件情況時(shí),濾波算法很難準(zhǔn)確完成地面點(diǎn)分類,因此在濾波處理后,結(jié)合目視解譯才能得到研究區(qū)準(zhǔn)確地面點(diǎn)。將機(jī)載LiDAR數(shù)據(jù)分為地面點(diǎn)和非地面點(diǎn),針對(duì)提取的地面點(diǎn),采用插值的方法生成數(shù)字高程模型(Digital Elevation Model,DEM),其分辨率為0.5m。使用DEM對(duì)點(diǎn)云進(jìn)行歸一化處理,提取特征變量。

表2 Ⅰ類誤差 Ⅱ類誤差 總誤差和kp系數(shù)的計(jì)算Tab.2 Calculation of typeⅠ error,typeⅡ error,total error and kp coefficient

2.3 特征變量提取

利用樣地矢量數(shù)據(jù)對(duì)預(yù)處理后的點(diǎn)云進(jìn)行裁剪,從裁剪后的樣地點(diǎn)云數(shù)據(jù)中提取34個(gè)特征變量(23個(gè)高度變量,10個(gè)密度變量和1個(gè)郁閉度變量)。各點(diǎn)云特征變量描述如表3所示。

表3 機(jī)載激光雷達(dá)數(shù)據(jù)中提取的特征變量指標(biāo)Tab.3 Feature variable indexes extracted from ALS data

2.4 建模方法

通過(guò)Python語(yǔ)言將機(jī)載LiDAR數(shù)據(jù)獲取的34個(gè)變量,利用隨機(jī)森林算法(RF)和迭代的決策樹(shù)算法(GBRT)結(jié)合樣地實(shí)測(cè)數(shù)據(jù)分別建模[25-26]。隨機(jī)森林算法(RF)由Breiman最早提出,是一種集成學(xué)習(xí)方法。其主要思想是為了讓機(jī)器學(xué)習(xí)具有更強(qiáng)的泛適性,從而在將無(wú)數(shù)的弱機(jī)器學(xué)習(xí)組合在一起形成的合集[27]。迭代的決策樹(shù)算法(GBRT)是一種由梯度提升(Gradient Boosting)、縮減(Shrinkage)、和回歸決策樹(shù)(Regression Decision Tree)等3個(gè)概念組成的新型迭代決策樹(shù)算法。該算法的最終結(jié)果是多棵決策樹(shù)的結(jié)果之和,最早由Jerome Friedman在1999年提出[28]。

通過(guò)上述兩種算法完成建模后采用獨(dú)立樣本檢驗(yàn)的方式來(lái)進(jìn)行模型校驗(yàn),隨機(jī)將樣本數(shù)據(jù)的70%用來(lái)建模,剩下的30%則作為檢驗(yàn)樣本。通過(guò)十折交叉驗(yàn)證法和網(wǎng)格尋優(yōu)法對(duì)建模參數(shù)進(jìn)行選擇,采用決定系數(shù)(R2)、均方根誤差(RMSE)和相對(duì)均方根誤差(rRMSE)來(lái)對(duì)模型進(jìn)行檢驗(yàn),從而比較隨機(jī)森林算法和新型迭代決策樹(shù)算法在機(jī)載LiDAR數(shù)據(jù)反演林分胸高斷面積時(shí)的適應(yīng)能力,以便選擇最優(yōu)樣本來(lái)參與建模。R2,RMSE和rRMSE公式:

(1)

(2)

(3)

3 結(jié)果分析

3.1 濾波方法對(duì)比

對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行濾波處理是激光LiDAR數(shù)據(jù)應(yīng)用的基礎(chǔ),旨在正確有效地區(qū)分地面點(diǎn)和非地面點(diǎn)。濾波效果的好壞將對(duì)后期提取的點(diǎn)云特征變量和森林參數(shù)估測(cè)產(chǎn)生影響。本研究隨機(jī)選擇4塊樣地,樣地信息如表4所示。

表4 樣地信息Tab.4 Sample plot information

通過(guò)比較PTIN,PMF,CSF和IBF濾波方法在4塊樣地中的濾波效果,來(lái)選擇適合具有普適性的方法對(duì)研究區(qū)點(diǎn)云進(jìn)行濾波。濾波效果如表5所示,當(dāng)林地坡度較大、林分郁閉度較低時(shí),IBF濾波效果相較于其他算法較好；當(dāng)?shù)匦纹骄?郁閉度較低時(shí),4種濾波方法的精度有一定差異。森林的郁閉度、坡度和點(diǎn)云密度對(duì)濾波效果的好壞都會(huì)產(chǎn)生影響,選擇濾波算法時(shí),應(yīng)考慮研究區(qū)環(huán)境特征。研究區(qū)大部分坡度在25～33°之間,相對(duì)平緩,郁閉度在0.5～0.7之間,雷達(dá)獲取數(shù)據(jù)時(shí)的點(diǎn)云密度約為2.9 point/m2,樣地2與研究區(qū)的情況基本一致。結(jié)合表5可知,樣地2中,CSF濾波算法對(duì)該樣地濾波效果最好。綜上,選用CSF濾波算法對(duì)研究區(qū)點(diǎn)云數(shù)據(jù)進(jìn)行濾波處理可得到最好效果。

表5 4種濾波算法的精度評(píng)價(jià)Tab.5 Accuracy evaluation of 4 filtering algorithms

3.2 特征選擇

從LiDAR數(shù)據(jù)中提取34個(gè)特征變量用于構(gòu)建RF和GBRT模型。研究表明,采用RF進(jìn)行變量篩選方法被廣泛應(yīng)用,且效果較好[2,20,30],故采用RF對(duì)LiDAR特征變量的重要性進(jìn)行排序篩選,胸高斷面積反演時(shí),選擇累計(jì)貢獻(xiàn)值達(dá)到80%的特征變量參與胸高斷面積反演。從圖2看出,高度變量占所選變量的比例較大,其中Hmax和H99重要性都達(dá)到0.17,累計(jì)貢獻(xiàn)到34%,說(shuō)明Hmax和H99與胸高斷面積具有較強(qiáng)的相關(guān)性。

圖2 特征變量重要性Fig.2 Importance of characteristic variables

3.3 模型效果分析

構(gòu)建反演模型時(shí),采用網(wǎng)格尋優(yōu)法和十折交叉驗(yàn)證法對(duì)RF反演參數(shù)進(jìn)行尋優(yōu)。結(jié)果顯示:當(dāng)構(gòu)造決策樹(shù)的數(shù)量為100,樹(shù)節(jié)點(diǎn)分裂的最小特征數(shù)為2,最小葉子節(jié)點(diǎn)為1時(shí),模型效果最優(yōu)。胸高斷面積模型的建模精度和驗(yàn)證精度如表6所示。

表6 胸高斷面積反演模型建模和檢驗(yàn)精度Tab.6 Modeling and test accuracy of inversion model of BA

從表6得出,建模結(jié)果中RF和GBRT模型的R2分別為0.77和0.83,RMSE為3.99m2/hm2和3.42m2/hm2,rRMSE為17.76%和15.22%；但在RF和GBRT模型進(jìn)行獨(dú)立樣本檢驗(yàn)的結(jié)果中，R2分別為0.66和0.64,RMSE為3.27m2/hm2和3.37m2/hm2,rRMSE為14.73%和15.18%。綜合考慮模型建模和檢驗(yàn)精度的R2,RMSE和rRMSE之間值的關(guān)系,RF模型建模精度和檢驗(yàn)精度相對(duì)穩(wěn)定,表明該模型的泛化能力優(yōu)于GBRT模型,故選用RF模型完成胸高斷面積估測(cè)。

3.4 胸高斷面積反演

綜合考慮模型的建模精度、檢驗(yàn)精度和模型的泛化能力等條件,最終決定采用RF模型實(shí)現(xiàn)研究區(qū)林木胸高斷面積的反演。胸高斷面積反演結(jié)果如圖3所示。

從圖3中看出:胸高斷面積在9～12m2/hm2區(qū)間的林分呈離散分布；胸高斷面積在12～18m2/hm2和18～23m2/hm2區(qū)間的林分面積最大,且兩者交錯(cuò)分布；胸高斷面積在23～30m2/hm2的林分總體面積不大且呈零星塊狀分布；胸高斷面積在>30(m2/hm2)區(qū)間的林分最少。上述結(jié)果與林地實(shí)際情況基本吻合。

圖3 胸高斷面積反演制圖Fig.3 Inversion mapping of BA

4 結(jié)論與討論

基于機(jī)載LiDAR數(shù)據(jù)進(jìn)行廣西高峰林場(chǎng)林分胸高斷面積估測(cè)研究,在點(diǎn)云數(shù)據(jù)的濾波處理中,對(duì)比了漸進(jìn)不規(guī)則三角網(wǎng)(PTIN)、漸進(jìn)形態(tài)學(xué)濾波算法(PMF)、布料模擬濾波算法(CSF)和基于插值的濾波方法(IBF),結(jié)果表明,CSF濾波算法效果最好。在濾波基礎(chǔ)上提取了34個(gè)點(diǎn)云特征變量,利用隨機(jī)森林算法(RF)和迭代的決策樹(shù)算法(GBRT)進(jìn)行林分胸高斷面積估測(cè),研究表明,RF模型估測(cè)效果優(yōu)于GBRT模型。在RF模型反演中,Hmax和H99兩個(gè)特征變量貢獻(xiàn)值最大,表明點(diǎn)云最大及99%分位高度變量與林分胸高斷面積呈顯著相關(guān)。這一結(jié)論與Silva等[16]利用機(jī)載LiDAR數(shù)據(jù)對(duì)火炬松人工林胸高斷面積進(jìn)行估測(cè)研究結(jié)論有一致性。

相關(guān)學(xué)者針對(duì)激光雷達(dá)估測(cè)森林胸高斷面積的研究表明,胸高斷面積不僅受到高度變量影響,還與林齡結(jié)構(gòu)、林分密度、平均樹(shù)高和平均冠幅等變量有關(guān)[29-30]。同時(shí),從測(cè)樹(shù)學(xué)角度來(lái)看,胸高斷面積的大小同林分密度有顯著相關(guān)性,這與本研究結(jié)論有差異,究其原因可以歸納為:1)高峰林場(chǎng)相較于其他學(xué)者的研究區(qū)在地理環(huán)境上存在差異,因此林齡結(jié)構(gòu)、林分密度、平均樹(shù)高等森林結(jié)構(gòu)因子有一定差異,導(dǎo)致所提取的主要特征變量存在差異；2)由于研究區(qū)多數(shù)為人工速生豐產(chǎn)林(如桉樹(shù)林),林分密度無(wú)較大差異,故在采用隨機(jī)森林對(duì)胸高斷面積進(jìn)行反演時(shí)體現(xiàn)林分密度的變量貢獻(xiàn)值較低,而Hmax,H95和H99等3個(gè)高度變量貢獻(xiàn)率較大。

本研究利用機(jī)器學(xué)習(xí)的方法進(jìn)行胸高斷面積反演,反演方法在模型機(jī)理表達(dá)方面欠佳,且研究區(qū)95%的林木為人工林,天然林范圍小,致使本研究在天然林以及大面積混交林上是否適用有待研究,期望在未來(lái)的研究中探索更多森林類型以及建模方法(如線性回歸和聯(lián)立方程組等)在胸高斷面積反演研究中的應(yīng)用潛力。

致謝:廣西壯族自治區(qū)林業(yè)勘測(cè)設(shè)計(jì)院,為本文的實(shí)驗(yàn)提供了機(jī)載 LiDAR 點(diǎn)云數(shù)據(jù)及野外調(diào)查地面數(shù)據(jù)。在此,表示衷心的感謝！