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

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

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

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

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

1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4.1 樣地蓄積計算 首先按照表3進行樹種組歸并，然后基于一元材積公式計算樣地總蓄積。蓄積量計算通用公式(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 激光雷達點云數(shù)據(jù)預處理 激光雷達原始點云數(shù)據(jù)經(jīng)姿態(tài)校正、噪聲點剔除、坐標轉(zhuǎn)換、航帶拼接、系統(tǒng)差改正等預處理形成1∶10 000比例尺分幅成果數(shù)據(jù)。本研究在分幅成果數(shù)據(jù)上利用LiDAR360軟件進行點云數(shù)據(jù)拼接、點云濾波、點云分類等，最終獲得用于建模的激光雷達點云數(shù)據(jù)。

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

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

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

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

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

3 建模方法

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

3.1 隨機森林算法

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

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

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

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

3.2 變量篩選

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

3.3 模型評價

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

4 結(jié)果與分析

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

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

表4 輸入變量Tab.4 Input variables

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

由圖4可知，當mtry=17時，模型具有最小的RMSE(18.01 m3·hm-2)、最小的rRMSE(14.02%)、最小的MAE(13.07 m3·hm-2)和相對較小的MRE(16.28%)； 當mtry=1時，模型具有最大的RMSE(20.17 m3·hm-2)、最大的rRMSE(15.97%)、最大的MAE(15.14 m3·hm-2)和最大的MRE(17.93%)。整體趨勢是： 當mtry≤17時，RMSE、rRMSE和MAE不斷減??； 當mtry>17時，RMSE、rRMSE、MAE和MRE顯著升高。所有模型的R2均在0.96附近，變化不明顯。故選取mtry=17作為最優(yōu)模型參數(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ù)和郁閉度建模，參與計算的變量為24個，1≤mtry≤24。圖6為24個模型的評價指標。

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

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

圖6 聯(lián)合高度參數(shù)和郁閉度建模評價指標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)隨機森林模型在訓練樣本和驗證樣本中的精度。由圖8可知，僅用高度參數(shù)建模，訓練階段的估測精度為R2=0.96、RMSE=18.01 m3·hm-2、MAE=13.07 m3·hm-2、rRMSE=14.02%、MRE=16.28%； 驗證階段的估測精度為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ù)和郁閉度建模，訓練階段的估測精度為R2=0.97、RMSE=16.94 m3·hm-2、MAE=12.44 m3·hm-2、rRMSE=13.18%、MRE=15.32%； 驗證階段的估測精度為R2=0.79、RMSE=36.23 m3·hm-2、MAE=26.16 m3·hm-2、rRMSE=32.73%、MRE=38.35%。

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

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

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

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

運用R軟件VSURF包對上述高度參數(shù)和郁閉度變量數(shù)據(jù)集進行變量篩選，過程如圖10所示。通過計算各變量重要性均值、變量重要性標準差和模型OOB誤差，在所有提取24個變量中，最終篩選出7個變量用于建模，分別為最大高(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給出了隨機森林模型在訓練樣本和驗證樣本中的估測精度。訓練階段的估測精度為R2=0.97、RMSE=17.24 m3·hm-2、MAE=12.76 m3·hm-2、rRMSE=13.42%、MRE=14.76%； 驗證階段的估測精度為R2=0.79、RMSE=36.50 m3·hm-2、MAE=26.08 m3·hm-2、rRMSE=32.97%、MRE=38.05%。

對比變量篩選前后的建模結(jié)果(圖9、圖12)可知，在模型訓練階段，經(jīng)變量篩選后的模型R2未降低，但是RMSE、rRMSE、MAE上升； 在模型驗證階段，經(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)變量篩選后，模型精度雖有變化，但是差別不大，因此可直接用篩選后的變量進行建模。

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

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

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

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

5 討論

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

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

6 結(jié)論

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