孫忠秋 高金萍 吳發(fā)云 高顯連 胡 楊 高劍新

(1.國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院 北京 100714； 2.寧夏大學(xué)生態(tài)環(huán)境學(xué)院 西北土地退化與生態(tài)恢復(fù)國家重點(diǎn)實(shí)驗(yàn)室培育基地西北退化生態(tài)系統(tǒng)恢復(fù)與重建教育部重點(diǎn)實(shí)驗(yàn)室 銀川 750021)

森林是陸地生態(tài)系統(tǒng)的主體，具有調(diào)節(jié)氣候、涵養(yǎng)水源、防風(fēng)固沙、減少污染、改善生態(tài)環(huán)境等重要作用，在山水林田湖草生命共同體中處于不可或缺的基礎(chǔ)地位。森林蓄積量指一定森林面積上全部樹木材積的總和(孟憲宇， 1996； Nilsson， 1996)，是反映一個(gè)國家或地區(qū)森林資源總規(guī)模和水平的基本指標(biāo)之一，與木材安全、氣候變化、動(dòng)物棲息等密切相關(guān)，可為制定森林經(jīng)營管理方案提供科學(xué)依據(jù)(李崇貴等， 2006； 徐濟(jì)德， 2014)，準(zhǔn)確估測(cè)森林蓄積量對(duì)提高森林資源管理水平和生態(tài)環(huán)境保護(hù)建設(shè)具有重要意義(N?sset， 2002； 陳新云等, 2019； 閆飛， 2014)。傳統(tǒng)大面積森林蓄積量估測(cè)主要根據(jù)國家森林資源規(guī)劃設(shè)計(jì)調(diào)查技術(shù)規(guī)程，通過對(duì)標(biāo)準(zhǔn)樣地相關(guān)因子的實(shí)測(cè)或用角規(guī)測(cè)量的方法，基于不同樹種材積公式計(jì)算樣地或小班蓄積量(李崇貴等， 2006)，進(jìn)而推算區(qū)域森林蓄積量； 此項(xiàng)工作精度高，但耗時(shí)耗力，在地形復(fù)雜的林區(qū)往往無法開展。近年來，隨著遙感應(yīng)用技術(shù)不斷發(fā)展，利用遙感影像結(jié)合地面樣地信息估測(cè)森林蓄積量成為可能。目前，大部分研究基于不同衛(wèi)星遙感影像數(shù)據(jù)源，如Landsat、Sentinel-2、高分系列、MODIS等(Chenetal.， 2012； Gireeetal.， 2013； 王月婷等， 2015； 楊柳等， 2017； 王佳等， 2014)，結(jié)合少量地面樣地調(diào)查數(shù)據(jù)，應(yīng)用數(shù)學(xué)算法建立森林蓄積量估測(cè)模型； 然而，由于光學(xué)遙感影像只能獲取森林的水平結(jié)構(gòu)信息，不具備獲取森林垂直結(jié)構(gòu)信息的能力，因此導(dǎo)致森林蓄積量估測(cè)精度普遍偏低。如李世波等(2019)基于GF-1影像估測(cè)森林蓄積量，模型估測(cè)精度(R2)僅0.50左右。

激光雷達(dá)是一種新興的主動(dòng)遙感技術(shù)，其突破了傳統(tǒng)光學(xué)遙感的光譜局限性，能夠穿透森林冠層，獲取森林三維結(jié)構(gòu)信息，在森林資源監(jiān)測(cè)中正逐漸發(fā)揮作用(Nilsson， 1996； N?sset， 1997； 曹林等， 2013； 李增元等， 2016)。學(xué)者們利用機(jī)載激光雷達(dá)數(shù)據(jù)反演林分平均高、生物量、郁閉度等森林參數(shù)，均取得了比傳統(tǒng)光學(xué)遙感精度更高的結(jié)果(N?ssetetal.， 2001； 2002； 2005； 段祝庚等， 2016； 耿林等， 2018)。湯旭光(2013)基于ICESat-GLAS大光斑激光雷達(dá)數(shù)據(jù)先估測(cè)林分平均高，再進(jìn)行森林生物量建模，得出針葉林樣地的估測(cè)精度(R2)為0.82。劉美爽等(2014)采用ICESat-GLAS數(shù)據(jù)對(duì)吉林省汪清林業(yè)局林區(qū)的林分冠層高度進(jìn)行估測(cè)，得出該區(qū)域的估測(cè)精度(R2)為0.84。吳迪等(2014)基于ICESat-GLAS數(shù)據(jù)，結(jié)合黑龍江省塔河林場(chǎng)109塊標(biāo)準(zhǔn)地調(diào)查數(shù)據(jù)，采用隨機(jī)森林算法對(duì)該地區(qū)林分平均高進(jìn)行估測(cè)，R2為0.72，RMSE為1.83 m。Hollaus等(2007)利用小光斑激光雷達(dá)數(shù)據(jù)估測(cè)奧地利福拉爾貝格州高山地區(qū)森林蓄積量，R2= 0.89，RMSE = 90.90 m3·hm-2。劉琪璟等(2008)基于日本長崎縣小光斑激光雷達(dá)數(shù)據(jù)估測(cè)林分平均高，其誤差為0.4～0.5 m。高婷等(2017)使用小光斑激光雷達(dá)數(shù)據(jù)估測(cè)甘肅張掖大野口林區(qū)林分平均高，R2=0.81。Sheridan等(2014)基于一元和多元線性回歸模型探索小光斑激光雷達(dá)數(shù)據(jù)估測(cè)美國俄勒岡州東部邁哈爾國家森林公園森林蓄積量的能力，得出一元和多元線性回歸的蓄積量估測(cè)模型R2分別為0.83和0.88。

由于ICESat-GLAS大光斑激光雷達(dá)數(shù)據(jù)可以免費(fèi)獲取，目前基于大光斑激光雷達(dá)數(shù)據(jù)進(jìn)行森林參數(shù)提取的研究較多，使用小光斑激光雷達(dá)數(shù)據(jù)進(jìn)行森林參數(shù)提取的研究較少。本研究以全覆蓋的有人機(jī)機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)和每木檢尺的地面樣地?cái)?shù)據(jù)為數(shù)據(jù)源，提取樣地點(diǎn)云高度參數(shù)和郁閉度等特征，采用隨機(jī)森林算法構(gòu)建森林蓄積量估測(cè)模型，同時(shí)對(duì)激光雷達(dá)生成的森林參數(shù)進(jìn)行變量篩選，確定蓄積量估測(cè)中重要的森林參數(shù)，最終確立以樣地為基本單元的森林蓄積量估測(cè)模型，分析機(jī)載激光雷達(dá)數(shù)據(jù)在森林蓄積量反演方面的潛力，以期為森林蓄積量高效準(zhǔn)確估測(cè)提供方法依據(jù)。

1 研究區(qū)概況

以大興溝林業(yè)局為研究區(qū)，該區(qū)位于吉林省東部、延邊朝鮮族自治州東北部，屬圖們江流域嘎呀河支流的中游(129°05′—130°01′E，43°20′—43°40′N)，總面積128 097 hm2。地勢(shì)東西兩端山高坡陡，山脈起伏較大，中部逐漸降低，東西溝中下游地段稍為平緩，平均坡度15°左右。屬溫帶大陸性季風(fēng)氣候，年平均氣溫2 ℃左右，無霜期105～125天。境內(nèi)植被屬長白山植物區(qū)系，大多為闊葉林和針闊混交林。人工林以落葉松(Larixspp.)林居多，其次為紅松(Pinuskoraiensis)林。常見的天然喬木樹種有紅松、云杉(Piceaasperata)、冷杉(Abiesfabri)等針葉樹種，黃檗(Phellodendronamurense)、水曲柳(Fraxinusmandshurica)、胡桃楸(Juglansmandshurica)、椴樹(Tiliatuan)、蒙古櫟(Quercusmongolica)、白樺(Betulaplatyphylla)等闊葉樹種。研究區(qū)內(nèi)森林資源豐富，是東北林區(qū)林業(yè)研究的重點(diǎn)區(qū)域之一(圖1)。

圖1 研究區(qū)森林及樣地分布Fig. 1 Forest and sample plot distribution of research area

2 數(shù)據(jù)及預(yù)處理

2.1 樣地調(diào)查數(shù)據(jù)

采用2018年10月“陸地碳衛(wèi)星吉林重點(diǎn)林區(qū)綜合試驗(yàn)地面樣地調(diào)查”項(xiàng)目獲取的數(shù)據(jù)，其中大興溝林業(yè)局區(qū)域共調(diào)查232塊半徑15 m的圓形樣地。

2.1.1 樣地位置 樣地位置對(duì)后續(xù)建模精度影響較大，為保證調(diào)查樣地類型具有代表性，樣地選擇主要基于以下原則： 1) 依據(jù)2016—2018年東北內(nèi)蒙古重點(diǎn)國有林區(qū)森林資源規(guī)劃設(shè)計(jì)調(diào)查數(shù)據(jù)，對(duì)小班中布設(shè)角規(guī)樣地樹高數(shù)據(jù)進(jìn)行分析，得到蓄積量排名前10的森林類型(云杉、冷杉、落葉松、樺木、楊樹、椴樹、櫟類、針葉混、針闊混和闊葉混)，并按照5個(gè)樹高級(jí)梯度(最大值與最小值之間劃分出5個(gè)區(qū)間)和3個(gè)郁閉度級(jí)梯度(低0.20～0.39、中0.40～0.69、高0.70以上)進(jìn)行劃分； 2) 選擇的樣地不在同一坡面或坡向； 3) 2塊樣地之間距離超過500 m； 4) 同類型單元樣地布設(shè)在不同起源的森林、不同地形的森林。

2.1.2 樣地調(diào)查 為保證后續(xù)樣地與激光雷達(dá)點(diǎn)云在空間上精確配準(zhǔn)，采用三基站聯(lián)合差分定位技術(shù)對(duì)樣地中心和樣木定位，并記錄坐標(biāo)。調(diào)查內(nèi)容包括胸徑、樹高、枝下高、冠幅、郁閉度、起源、樹種組成等信息。利用R軟件從232塊樣地中隨機(jī)抽取70%數(shù)據(jù)作為訓(xùn)練樣本(164塊樣地)，剩余30%數(shù)據(jù)作為驗(yàn)證樣本(68塊樣地)，樣地統(tǒng)計(jì)信息和分組情況見表1。

表1 樣地蓄積量統(tǒng)計(jì)信息Tab.1 Sample stand stock volume statistics

2.2 有人機(jī)機(jī)載激光雷達(dá)數(shù)據(jù)

有人機(jī)機(jī)載激光雷達(dá)數(shù)據(jù)于2018年8月由搭載在塞斯納208B有人機(jī)平臺(tái)上的RIEGL-VQ-1560i激光雷達(dá)航攝儀獲取。該設(shè)備是超高性能、高度集成的雙通道機(jī)載激光掃描儀系統(tǒng)，能夠在不同高度飛行作業(yè)獲取大范圍高密度點(diǎn)云，適合對(duì)大面積區(qū)域和復(fù)雜環(huán)境進(jìn)行空中測(cè)繪。本研究中，飛機(jī)設(shè)計(jì)航飛相對(duì)航高1 800 m，飛行速度(相對(duì)地面)240 km·h-1，激光器發(fā)射頻率2 000 kHz(兩通道同時(shí)工作，單通道頻率1 000 kHz)，此參數(shù)下獲取的激光點(diǎn)云密度為每平方米10個(gè)點(diǎn)； 但由于地形起伏因素影響，最終成果約每平方米13個(gè)點(diǎn)。有人機(jī)機(jī)載激光雷達(dá)數(shù)據(jù)概況見表2。

表2 有人機(jī)機(jī)載激光雷達(dá)數(shù)據(jù)概況Tab.2 An overview of airborne LiDAR data

數(shù)據(jù)獲取時(shí)采用的測(cè)繪基準(zhǔn)如下： 1) 坐標(biāo)系統(tǒng)， 2000國家大地坐標(biāo)系； 2) 高程基準(zhǔn)， 1985國家高程基準(zhǔn)； 3) 投影方式， 高斯克呂格投影，3°分帶，東偏500 km，加帶號(hào)，中央子午線129°。激光雷達(dá)數(shù)據(jù)對(duì)研究區(qū)全覆蓋，總面積約1.2×105hm2。數(shù)據(jù)存儲(chǔ)采用LAS 1.2格式的點(diǎn)云，數(shù)據(jù)量約1.72 TB。

2.3 數(shù)據(jù)準(zhǔn)備和建模流程

在數(shù)據(jù)準(zhǔn)備階段，首先對(duì)研究區(qū)點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，去除異常點(diǎn)并分類，分離地面點(diǎn)； 然后進(jìn)行樣地和點(diǎn)云的幾何配準(zhǔn)，該過程需選取大量同名點(diǎn)位，因研究區(qū)森林覆蓋率高，同名點(diǎn)選取困難，耗時(shí)較長。完成樣地點(diǎn)云數(shù)據(jù)裁切后，即進(jìn)入數(shù)據(jù)建模階段，分別提取訓(xùn)練樣地和驗(yàn)證樣地點(diǎn)云的高度參數(shù)和郁閉度，采用隨機(jī)森林算法建模并進(jìn)行精度評(píng)價(jià)。數(shù)據(jù)準(zhǔn)備和建模流程如圖2所示。

圖2 蓄積量估測(cè)流程Fig. 2 Flow chart of forest stock volume estimation

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

2.4.1 樣地蓄積計(jì)算 首先按照表3進(jìn)行樹種組歸并，然后基于一元材積公式計(jì)算樣地總蓄積。蓄積量計(jì)算通用公式(E·N·楚里克等， 1989)如下：

表3 樹種組基本信息Tab.3 Tree species group basic information

式中：V為樣地總蓄積；N為樣木株數(shù)；D為樣木胸徑；a、b、c、d、e、k為材積公式中的系數(shù)，根據(jù)《中國立木材積表》吉林省立木材積表(劉琪璟， 2017)確定。

2.4.2 激光雷達(dá)點(diǎn)云數(shù)據(jù)預(yù)處理 激光雷達(dá)原始點(diǎn)云數(shù)據(jù)經(jīng)姿態(tài)校正、噪聲點(diǎn)剔除、坐標(biāo)轉(zhuǎn)換、航帶拼接、系統(tǒng)差改正等預(yù)處理形成1∶10 000比例尺分幅成果數(shù)據(jù)。本研究在分幅成果數(shù)據(jù)上利用LiDAR360軟件進(jìn)行點(diǎn)云數(shù)據(jù)拼接、點(diǎn)云濾波、點(diǎn)云分類等，最終獲得用于建模的激光雷達(dá)點(diǎn)云數(shù)據(jù)。

2.4.3 樣地點(diǎn)云數(shù)據(jù)提取高度參數(shù) 基于樣地點(diǎn)云數(shù)據(jù)，采用數(shù)學(xué)統(tǒng)計(jì)方法可以提取多個(gè)與高度相關(guān)的參數(shù)，即根據(jù)指定的高度間隔將其進(jìn)一步分割成不同的“層”，統(tǒng)計(jì)各層的點(diǎn)數(shù)。建模中常用的高度參數(shù)有最大高(Hmax)、最小高(Hmin)、平均高(Hmean)、高度中位數(shù)(Hmedian)、高度百分位數(shù)(H%)和高度標(biāo)準(zhǔn)差(Hstd)等。其中H%的計(jì)算方法如下： 某一統(tǒng)計(jì)單元內(nèi)，將其內(nèi)部所有歸一化的激光雷達(dá)點(diǎn)云按高度排序，計(jì)算每一統(tǒng)計(jì)單元內(nèi)X%的點(diǎn)所在高度即為該統(tǒng)計(jì)單元的高度百分位數(shù)(圖3)。建模中使用的高度百分位數(shù)包含15個(gè)，即1%、5%、10%、20%、25%、30%、40%、50%、60%、70%、75%、80%、90%、95%和99%。

圖3 激光雷達(dá)點(diǎn)云高度百分位數(shù)Fig. 3 Height percentile of LiDAR point cloud

2.4.4 樣地點(diǎn)云數(shù)據(jù)提取郁閉度 樣地點(diǎn)云的郁閉度定義為植被回波點(diǎn)數(shù)與總點(diǎn)數(shù)的比值。本研究將植被高度閾值設(shè)為2 m，在計(jì)算過程中大于高度閾值的點(diǎn)均被認(rèn)為是植被點(diǎn)，公式如下：

式中：Pc為郁閉度；Nveg為植被回波點(diǎn)數(shù)；Ntotal為總點(diǎn)數(shù)。

輸出值范圍為0(沒有林冠層覆蓋/完全裸露)～1(全植被覆蓋)。提取的郁閉度與樣地實(shí)測(cè)郁閉度作差，差值均值為0.08，標(biāo)準(zhǔn)差為0.14，二者具有很好的一致性。

3 建模方法

蓄積量模型的構(gòu)建一般使用兩大類方法(Shaoetal.， 2017； Takmaetal.， 2012)。一是參數(shù)化方法，構(gòu)建由有限數(shù)量的參數(shù)定義或參數(shù)化的模型，該方法需要作出某些先驗(yàn)假設(shè)，且采用測(cè)試數(shù)據(jù)以確保不違反這些假設(shè)，有時(shí)還需要對(duì)變量進(jìn)行適當(dāng)轉(zhuǎn)換。參數(shù)化方法能夠很好解釋待測(cè)參數(shù)與變量之間的相關(guān)關(guān)系，但缺點(diǎn)是獲取新數(shù)據(jù)必須重新構(gòu)建新模型。在眾多蓄積量估測(cè)研究中，多元線性回歸是最常用的方法，通過對(duì)實(shí)測(cè)蓄積量與激光雷達(dá)提取變量之間的關(guān)系進(jìn)行回歸分析，得到蓄積量估測(cè)模型，利用驗(yàn)證樣地?cái)?shù)據(jù)和一系列檢驗(yàn)指標(biāo)可以驗(yàn)證模型精度。二是非參數(shù)化方法，與參數(shù)化方法相比，該方法無需先驗(yàn)假設(shè)，模型構(gòu)建更便捷。機(jī)器學(xué)習(xí)算法是典型的非參數(shù)化模型構(gòu)建方法，盡管該類算法不能得出具體模型，但并不影響算法的回歸預(yù)測(cè)，且機(jī)器學(xué)習(xí)算法的預(yù)測(cè)結(jié)果往往高于傳統(tǒng)的線性回歸方法(García-Gutiérrezetal.， 2015； Yuetal.， 2008； 2011)。因此，本研究采用機(jī)器學(xué)習(xí)算法中的隨機(jī)森林算法構(gòu)建森林蓄積量估測(cè)模型。

3.1 隨機(jī)森林算法

隨機(jī)森林算法是由Breiman(2001)提出的，其具體建模步驟如下： 對(duì)輸入隨機(jī)森林模型的訓(xùn)練樣本進(jìn)行隨機(jī)采樣，包括行(單個(gè)樣本)和列(特征變量)2個(gè)維度。行采樣使用Bootstrap算法，列采樣為從M個(gè)特征變量中隨機(jī)選擇mtry(mtry≤M)個(gè)特征變量?；谏鲜霾襟E，構(gòu)建k株決策樹，通過求取平均值得到最終預(yù)測(cè)結(jié)果。隨機(jī)森林回歸公式可表示為Y=Eθh(X,θ)。利用袋外數(shù)據(jù)(out of bag, OOB)計(jì)算每株決策樹的預(yù)測(cè)誤差，每株決策樹的平均泛化誤差(generalization error,GE)計(jì)算公式為：

GE=EθEX,Y[Y-h(X,θ)]2。

式中：θ為隨機(jī)變量；Eθ為期望函數(shù)；X、Y為訓(xùn)練集抽取的隨機(jī)變量；h為決策樹預(yù)測(cè)函數(shù)；EX,Y為X、Y的聯(lián)合期望函數(shù)。

使用R軟件Random Forest 包，經(jīng)多次試驗(yàn)確定決策樹數(shù)目(ntree)和節(jié)點(diǎn)分裂時(shí)變量個(gè)數(shù)(mtry)，代入隨機(jī)森林回歸模型，并利用回歸模型對(duì)驗(yàn)證樣本進(jìn)行預(yù)估。

3.2 變量篩選

建模過程中變量個(gè)數(shù)越多，運(yùn)算量越大，耗時(shí)越長，通常需要進(jìn)行有效的變量篩選，以在不損失建模精度的前提下降低運(yùn)算量。本研究運(yùn)用R軟件VSURF(variable selection using random forests)包進(jìn)行變量篩選，主要包含3個(gè)步驟。1) 閾值處理： 首先，在ntree和mtry默認(rèn)設(shè)置下，采用隨機(jī)森林算法的重要值函數(shù)計(jì)算參數(shù)nfor.thres rf； 然后，按變量的平均變量重要性(variable importance, VI)降序排序； 接下來，計(jì)算閾值min.thres，修剪的CART樹的最小預(yù)測(cè)值擬合到VI的標(biāo)準(zhǔn)偏差曲線； 最后，計(jì)算閾值，僅保留平均VI大于nmin*min.thres的變量。2) 解釋： 考慮步驟1選擇的變量，采用nfor.interp函數(shù)，首先，只選擇最重要的變量，直到計(jì)算完第一步選擇的所有變量結(jié)束； 然后，采用err.min函數(shù)計(jì)算模型的最小平均袋外(OOB)誤差及其相關(guān)的標(biāo)準(zhǔn)偏差sd.min； 最后，選擇平均OOB誤差小于err.min+nsd*sd.min的最小模型(及其相應(yīng)的變量)。3) 預(yù)測(cè)： 起點(diǎn)與步驟2相同，但是現(xiàn)在變量逐步添加到模型中，使用步驟2遺漏的變量和mean.jump計(jì)算平均跳躍值，并設(shè)置為一個(gè)模型的平均OOB誤差與其第一個(gè)跟隨模型之間的平均絕對(duì)差值，如果平均OOB誤差減小大于nmj*mean.jump，則變量包含在模型中。通過VSURF包設(shè)置的3個(gè)步驟，對(duì)蓄積量估測(cè)相對(duì)不重要的變量會(huì)被移除，篩選后的變量用于構(gòu)建森林蓄積量估測(cè)模型。

3.3 模型評(píng)價(jià)

采用決定系數(shù)(coefficient of determination,R2)、均方根誤差(root mean square error, RMSE)評(píng)價(jià)模型擬合優(yōu)度(Huetal.， 2019)，采用相對(duì)均方根誤差(relative RMSE, rRMSE)、平均絕對(duì)誤差(mean absolute error, MAE)和平均相對(duì)誤差(mean relative error, MRE)評(píng)價(jià)模型估測(cè)精度。決定系數(shù)(R2)為自變量變異程度占總變異程度的比例，R2越大，表示模型擬合能力越強(qiáng)。均方根誤差(RMSE)為標(biāo)準(zhǔn)誤差的算術(shù)平方根，用于衡量預(yù)測(cè)值與真實(shí)值之間的偏差，RMSE越小，表示模型預(yù)測(cè)能力越強(qiáng)。相對(duì)均方根誤差(rRMSE)為RMSE與估測(cè)結(jié)果算術(shù)平均值的比值，與評(píng)價(jià)量本身數(shù)量級(jí)水平無關(guān)，能夠更好體現(xiàn)總體值域差別較大的模型預(yù)測(cè)精度(張瑞英等， 2016)，rRMSE越小，表示模型預(yù)測(cè)效果越好。平均絕對(duì)誤差(MAE)為絕對(duì)誤差的平均值，可以反映預(yù)測(cè)值誤差的實(shí)際情況。平均相對(duì)誤差(MRE)可以反映模型的預(yù)估精度。5個(gè)指標(biāo)的計(jì)算公式如下：

4 結(jié)果與分析

4.1 輸入?yún)?shù)與最優(yōu)模型參數(shù)選擇

本研究對(duì)比2種情況下的建模效果，具體輸入變量如表4所示。

表4 輸入變量Tab.4 Input variables

隨機(jī)森林模型好壞由mtry和ntree決定。對(duì)于回歸問題，mtry 默認(rèn)設(shè)置為全部自變量數(shù)量的1/3(取整)，ntree一般取500(本研究也取值500)。通常，mtry 取默認(rèn)值不一定能獲取最優(yōu)模型，選擇合適的mtry可以降低隨機(jī)森林模型的預(yù)測(cè)誤差(歐強(qiáng)新等， 2019)，因此本研究對(duì)mtry進(jìn)行調(diào)優(yōu)。利用高度參數(shù)建模，參與計(jì)算的變量為23個(gè)，1≤mtry≤23。圖4為23個(gè)模型的評(píng)價(jià)指標(biāo)。

由圖4可知，當(dāng)mtry=17時(shí)，模型具有最小的RMSE(18.01 m3·hm-2)、最小的rRMSE(14.02%)、最小的MAE(13.07 m3·hm-2)和相對(duì)較小的MRE(16.28%)； 當(dāng)mtry=1時(shí)，模型具有最大的RMSE(20.17 m3·hm-2)、最大的rRMSE(15.97%)、最大的MAE(15.14 m3·hm-2)和最大的MRE(17.93%)。整體趨勢(shì)是： 當(dāng)mtry≤17時(shí)，RMSE、rRMSE和MAE不斷減?。?當(dāng)mtry>17時(shí)，RMSE、rRMSE、MAE和MRE顯著升高。所有模型的R2均在0.96附近，變化不明顯。故選取mtry=17作為最優(yōu)模型參數(shù)，此時(shí)，ntree最優(yōu)參數(shù)值為64(圖5)。

圖5 高度參數(shù)建模ntree最優(yōu)參數(shù)值確定Fig. 5 Optimal ntree parameter of height parameter modeling

聯(lián)合高度參數(shù)和郁閉度建模，參與計(jì)算的變量為24個(gè)，1≤mtry≤24。圖6為24個(gè)模型的評(píng)價(jià)指標(biāo)。

據(jù)《財(cái)經(jīng)》報(bào)道，阿里內(nèi)部對(duì)于楊偉東的調(diào)查已經(jīng)持續(xù)一段時(shí)間，楊偉東此次涉嫌貪腐的項(xiàng)目主要集中在優(yōu)酷于2018年推出的“這就是”系列綜藝，主要是關(guān)于綜藝項(xiàng)目的收支問題。該系列綜藝以《這！就是街舞》打頭陣，后陸續(xù)推出了《這！就是灌籃》等。

由圖6可知，當(dāng)mtry=22時(shí)，模型具有最小的RMSE(16.94 m3·hm-2)、相對(duì)較小的rRMSE(13.18%)、最小的MAE(12.44 m3·hm-2)和最小的MRE(15.32%)； 當(dāng)mtry=1時(shí)，模型具有最大的RMSE(19.46 m3·hm-2)、最大的rRMSE(15.45%)、最大的MAE(14.47 m3·hm-2)和最大的MRE(17.59%)。整體趨勢(shì)是：R2變化不明顯，當(dāng)mtry≤22時(shí)，RMSE、rRMSE、MAE和MRE不斷減??； 當(dāng)mtry>22時(shí)，RMSE、rRMSE、MAE和MRE顯著升高。故選取mtry=22作為最優(yōu)模型參數(shù)，此時(shí)，ntree最優(yōu)參數(shù)值為406(圖7)。

圖6 聯(lián)合高度參數(shù)和郁閉度建模評(píng)價(jià)指標(biāo)Fig. 6 Evaluation indicators of height parameter and crown density modeling

圖7 聯(lián)合高度參數(shù)和郁閉度建模ntree最優(yōu)參數(shù)值確定Fig. 7 Optimal ntree parameter of height parameter and crown density modeling

4.2 建模結(jié)果

圖8、9分別給出了上述最優(yōu)隨機(jī)森林模型在訓(xùn)練樣本和驗(yàn)證樣本中的精度。由圖8可知，僅用高度參數(shù)建模，訓(xùn)練階段的估測(cè)精度為R2=0.96、RMSE=18.01 m3·hm-2、MAE=13.07 m3·hm-2、rRMSE=14.02%、MRE=16.28%； 驗(yàn)證階段的估測(cè)精度為R2=0.75、RMSE=40.07 m3·hm-2、MAE=29.21 m3·hm-2、rRMSE=36.20%、MRE=49.40%。由圖9可知，聯(lián)合高度參數(shù)和郁閉度建模，訓(xùn)練階段的估測(cè)精度為R2=0.97、RMSE=16.94 m3·hm-2、MAE=12.44 m3·hm-2、rRMSE=13.18%、MRE=15.32%； 驗(yàn)證階段的估測(cè)精度為R2=0.79、RMSE=36.23 m3·hm-2、MAE=26.16 m3·hm-2、rRMSE=32.73%、MRE=38.35%。

圖8 僅用高度參數(shù)(23個(gè)變量)估測(cè)森林蓄積量(左： 訓(xùn)練模型結(jié)果； 右： 驗(yàn)證模型結(jié)果)Fig. 8 Estimation forest stock volume only using height parameters(left： training result; right： validation result)

圖9 聯(lián)合高度參數(shù)和郁閉度(24個(gè)變量)估測(cè)森林蓄積量(左： 訓(xùn)練模型結(jié)果； 右： 驗(yàn)證模型結(jié)果)Fig. 9 Estimation forest stock volume using height parameters and canopy density(left： training result; right： validation result)

對(duì)比分析2種情況下的建模結(jié)果可知，增加郁閉度信息后，模型R2升高，RMSE、rRMSE、MAE和MRE均顯著下降?？梢娫谏中罘e量建模時(shí)，增加林分郁閉度信息能夠提升模型精度。

4.3 變量篩選與建模結(jié)果

運(yùn)用R軟件VSURF包對(duì)上述高度參數(shù)和郁閉度變量數(shù)據(jù)集進(jìn)行變量篩選，過程如圖10所示。通過計(jì)算各變量重要性均值、變量重要性標(biāo)準(zhǔn)差和模型OOB誤差，在所有提取24個(gè)變量中，最終篩選出7個(gè)變量用于建模，分別為最大高(Hmax)、平均高(Hmean)、郁閉度(Pc)、50%高度百分位數(shù)(H%8)、60%高度百分位數(shù)(H%9)、高度的二次冪平均(Hsqrt_mean_sq)和高度方差(Hvar)。篩選后變量的重要性排序如圖11所示。

圖10 基于VSURF包篩選變量Fig. 10 Variable selection based on VSURF package

圖11 篩選后變量的重要性排序Fig. 11 The importance of selecting variables

在最優(yōu)篩選變量情況下，圖12給出了隨機(jī)森林模型在訓(xùn)練樣本和驗(yàn)證樣本中的估測(cè)精度。訓(xùn)練階段的估測(cè)精度為R2=0.97、RMSE=17.24 m3·hm-2、MAE=12.76 m3·hm-2、rRMSE=13.42%、MRE=14.76%； 驗(yàn)證階段的估測(cè)精度為R2=0.79、RMSE=36.50 m3·hm-2、MAE=26.08 m3·hm-2、rRMSE=32.97%、MRE=38.05%。

對(duì)比變量篩選前后的建模結(jié)果(圖9、圖12)可知，在模型訓(xùn)練階段，經(jīng)變量篩選后的模型R2未降低，但是RMSE、rRMSE、MAE上升； 在模型驗(yàn)證階段，經(jīng)變量篩選后的模型R2未變化，RMSE從36.23 m3·hm-2升至36.50 m3·hm-2，rRMSE從31.92%升至32.97%，MAE從26.16 m3·hm-2降至26.08 m3·hm-2，MRE從38.35%降至38.05%?？梢?，經(jīng)變量篩選后，模型精度雖有變化，但是差別不大，因此可直接用篩選后的變量進(jìn)行建模。

為了進(jìn)一步驗(yàn)證隨機(jī)森林算法在訓(xùn)練和驗(yàn)證階段所得結(jié)果的穩(wěn)定性，本研究額外增加10次獨(dú)立重復(fù)試驗(yàn)與圖12結(jié)果進(jìn)行對(duì)比，10次獨(dú)立重復(fù)試驗(yàn)的模型訓(xùn)練階段結(jié)果如表5所示、驗(yàn)證階段結(jié)果如表6所示。由表5可知，R2的平均值為0.96，RMSE的平均值為17.36 m3·hm-2，MAE的平均值為12.56 m3·hm-2，rRMSE的平均值為13.84%，MRE的平均值為14.00%，與圖12訓(xùn)練模型結(jié)果基本一致； 由表6可知，R2的平均值為0.78，RMSE的平均值為40.30 m3·hm-2，MAE的平均值為28.74 m3·hm-2，rRMSE的平均值為33.57%，MRE的平均值為34.39%，與圖12驗(yàn)證模型結(jié)果也基本一致。這說明，數(shù)據(jù)隨機(jī)分組后在隨機(jī)森林算法下的建模預(yù)測(cè)結(jié)果具有很好的一致性。

圖12 基于篩選變量估測(cè)森林蓄積量(左： 訓(xùn)練模型結(jié)果； 右： 驗(yàn)證模型結(jié)果)Fig. 12 Estimation forest stock volume using selection variables(left： training result; right： validation result)

表5 10次建模訓(xùn)練結(jié)果Tab.5 Ten times results in the training phases

表6 10次建模驗(yàn)證結(jié)果Tab.6 Ten times results in the validation phases

5 討論

采用隨機(jī)森林算法對(duì)基于機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)提取的不同變量進(jìn)行建模，模型精度均很高，隨機(jī)森林算法的適應(yīng)性較強(qiáng)。僅用高度參數(shù)建模的估測(cè)精度為R2=0.75、RMSE=40.07 m3·hm-2、MAE=29.21 m3·hm-2、MRE=49.40%，聯(lián)合高度參數(shù)和郁閉度建模的估測(cè)精度為R2=0.79、RMSE=36.23 m3·hm-2、MAE=26.16 m3·hm-2、MRE=38.35%，說明基于機(jī)載激光雷達(dá)點(diǎn)云估測(cè)森林蓄積量時(shí)，增加林分郁閉度信息是提高建模精度的重要途徑。Hu等(2020)基于哨兵2號(hào)光譜數(shù)據(jù)聯(lián)合機(jī)器學(xué)習(xí)算法和多元線性回歸方法估測(cè)森林蓄積量，得到的最優(yōu)建模精度僅為R2=0.58、RMSE=65.03 m3·hm-2。Chrysafis等(2017)探討Sentinel-2和Landsat-8估測(cè)森林蓄積量的能力，結(jié)果分別為R2=0.63、RMSE=63.11 m3·hm-2和R2=0.62、RMSE=64.40 m3·hm-2。相較其他研究，本研究在基于激光雷達(dá)點(diǎn)云數(shù)據(jù)提取的變量下，不同形式的變量組合均取得了較好結(jié)果。與傳統(tǒng)光學(xué)遙感影像相比，取得較好結(jié)果可能基于以下3點(diǎn)原因： 一是激光雷達(dá)穿透能力較強(qiáng)，能穿透林分到達(dá)地面，直接獲取林分高度等垂直結(jié)構(gòu)信息，與光學(xué)遙感影像只能獲取林分水平結(jié)構(gòu)信息不同，垂直結(jié)構(gòu)信息估測(cè)森林蓄積量更有效(Limetal.， 2009)； 二是本研究采用參數(shù)優(yōu)化的隨機(jī)森林模型，以RMSE為標(biāo)準(zhǔn)分別計(jì)算模型的最優(yōu)變量，使得模型具有更好的估測(cè)能力(歐強(qiáng)新等， 2019)； 三是本研究區(qū)位于我國東北地區(qū)，地形和森林結(jié)構(gòu)相較于南方等研究區(qū)可能比較簡(jiǎn)單(李崇貴等, 2006)。目前也有一些基于激光雷達(dá)點(diǎn)云數(shù)據(jù)的森林參數(shù)研究，如龐勇等(2011)采用小光斑激光雷達(dá)數(shù)據(jù)、ICESat-GLAS大光斑激光雷達(dá)數(shù)據(jù)和光學(xué)遙感影像數(shù)據(jù)，以78塊地面樣地作為真值，對(duì)大湄公河次區(qū)域森林生物量進(jìn)行估測(cè)，建模估測(cè)相關(guān)系數(shù)為0.70，相比本研究結(jié)果略差，可能是南方地區(qū)復(fù)雜的地形環(huán)境和較少的樣本量造成的。曹林等(2014)以江蘇常熟虞山林場(chǎng)為研究區(qū)，采用逐步回歸方法得到森林蓄積量最優(yōu)估測(cè)結(jié)果的決定系數(shù)僅為0.55，其精度較低的原因主要有3點(diǎn)： 1) 模型選擇得不好，與隨機(jī)森林算法相比，逐步回歸方法的表現(xiàn)能力較差，這是造成其估測(cè)精度較低的最主要原因； 2) 樣本量較少，其研究?jī)H有73塊地面樣地?cái)?shù)據(jù)，數(shù)據(jù)在建模時(shí)代表性可能不足； 3) 相較于北方地區(qū)，江蘇常熟虞山林場(chǎng)地形可能更復(fù)雜。劉浩等(2018)基于江蘇東臺(tái)林場(chǎng)55塊地面樣地?cái)?shù)據(jù)，得出小光斑激光雷達(dá)數(shù)據(jù)估測(cè)該地區(qū)人工林森林蓄積量精度的調(diào)整R2=0.84、rRMSE=14.27%，結(jié)果要略好于本研究，主要原因是其研究對(duì)象為人工林，林型相對(duì)較簡(jiǎn)單，估測(cè)信息提取較好。Holmgren(2004)對(duì)比研究不同點(diǎn)云密度估測(cè)瑞典西南部地區(qū)森林蓄積量的能力，得到蓄積量的最優(yōu)估測(cè)均方根誤差為31.00 m3·hm-2，略優(yōu)于本研究(36.23 m3·hm-2)，說明點(diǎn)云密度是影響蓄積量估測(cè)結(jié)果的一個(gè)重要因素。

此外，本研究充分評(píng)估了郁閉度對(duì)森林蓄積量估測(cè)的影響，結(jié)果發(fā)現(xiàn)，增加林分郁閉度信息可提高模型估測(cè)森林蓄積量的能力，R2由0.75提高到0.79，RMSE從40.07 m3·hm-2降至36.23 m3·hm-2，說明采用機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)反演森林蓄積量時(shí)，輔助林分郁閉度增加數(shù)據(jù)水平結(jié)構(gòu)信息，能夠取得更好的估測(cè)結(jié)果。同時(shí)，本研究還評(píng)估了篩選變量對(duì)森林蓄積量估測(cè)的影響，結(jié)果表明，雖然通過變量篩選能夠降低參數(shù)數(shù)量，由原來的24個(gè)減少至7個(gè)，可極大提高運(yùn)算效率，但R2未變化，RMSE從36.23 m3·hm-2升至36.50 m3·hm-2，rRMSE從31.92%升至32.97%，MAE從26.16 m3·hm-2降至26.08 m3·hm-2，MRE從38.35%降至38.05%，說明經(jīng)變量篩選后，模型精度雖有變化，但是差別不大，可直接用篩選后的變量進(jìn)行建模。因此，處理區(qū)域大尺度問題時(shí)，在數(shù)據(jù)量大或運(yùn)算能力不足的情況下，建議使用篩選變量建模，這樣可在稍微犧牲精度的情況下盡可能縮短計(jì)算時(shí)間。

6 結(jié)論

本研究基于機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)提取的森林高度參數(shù)和郁閉度，結(jié)合分層地面樣地調(diào)查數(shù)據(jù)，采用隨機(jī)森林算法構(gòu)建森林蓄積量估測(cè)模型，結(jié)果表明，增加林分郁閉度信息可顯著提高森林蓄積量估測(cè)精度。通過變量篩選，雖然能夠降低參數(shù)數(shù)量，但對(duì)模型精度具有一定影響。在建模精度要求較高的情況下，建議使用全變量進(jìn)行蓄積量估測(cè)； 而在數(shù)據(jù)量較大的情況下，建議使用篩選變量進(jìn)行蓄積量估測(cè)?；跈C(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)估測(cè)森林蓄積量顯著優(yōu)于光學(xué)遙感數(shù)據(jù)，可為森林蓄積量高效準(zhǔn)確估測(cè)提供方法依據(jù)，能夠滿足大范圍森林蓄積量快速反演需求。