基于機載激光雷達點云的桉樹林分蓄積量估算模型構(gòu)建

鄧 焯，李 斌，范光鵬，趙天忠，于永輝

（1. 北京林業(yè)大學 信息學院，北京 100083；2. 北京林業(yè)大學 國家林業(yè)草原林業(yè)智能信息處理工程技術(shù)研究中心，北京 100083；3. 北京林業(yè)大學 林業(yè)信息化研究所，北京 100083；4. 廣西壯族自治區(qū)國有高峰林場，廣西 南寧 530001）

森林是陸地生態(tài)系統(tǒng)的主體，在全球碳循環(huán)中起著十分重要的作用。森林蓄積量作為衡量森林數(shù)量的重要的指標之一，能夠直觀反映森林資源數(shù)量和質(zhì)量，獲取森林蓄積量是推算森林生物量和碳儲量的必要環(huán)節(jié)[1?3]?！丁笆奈濉绷謽I(yè)草原保護發(fā)展規(guī)劃綱要》明確指出：將增加中國森林蓄積量作為主要目標之一，并將森林蓄積量與森林覆蓋率作為兩大約束性指標。因此，實現(xiàn)森林資源動態(tài)檢測、準確獲取森林蓄積量等參數(shù)信息，已成為當前森林資源調(diào)查的迫切需求。激光雷達(light detection and ranging,LiDAR)作為一種新興的主動式遙感技術(shù)，能夠從不同空間尺度對生態(tài)系統(tǒng)進行高效精準的監(jiān)測[4]。機載LiDAR發(fā)射的激光脈沖能穿透森林冠層并獲取樹冠上部枝葉的空間信息，通過林分內(nèi)的空隙，測量冠層結(jié)構(gòu)信息和林下地形信息[5]。而傳統(tǒng)方法通過樣木、樣地等抽樣調(diào)查進行林分森林參數(shù)的推算，調(diào)查周期長、成本高且調(diào)查范圍有限[6?7]。激光雷達在獲取森林空間結(jié)構(gòu)因子和地形因子等信息方面具有精度高、范圍廣等優(yōu)勢。

近年來，基于機載激光雷達的林分蓄積量反演已有許多研究成果，集中在構(gòu)建參數(shù)方法和非參數(shù)方法的蓄積量反演模型。PAWE等[8]基于機載LiDAR數(shù)據(jù)對波蘭東南部林區(qū)建立多元線性回歸蓄積量反演模型，其中均方根誤差(root mean square error, RMSE)為15.2%；CHIRICI等[9]以Landsat 5TM、衛(wèi)星LiDAR數(shù)據(jù)等預測變量聯(lián)合氣溫、降水和地形等輔助變量進行意大利中部地區(qū)蓄積量的大尺度反演，其中隨機森林回歸模型最優(yōu)，決定系數(shù)(R2)為0.69、RMSE為37.2%；陳松等[10]基于Sentinel-2與機載LiDAR數(shù)據(jù)采用不同回歸方法對廣西高峰林場界牌、東升分場進行蓄積量反演，構(gòu)建MLR-Logistic聯(lián)立模型精度優(yōu)于隨機森林等機器學習方法，R2為0.60、相對均方根誤差(relative root mean-squared error,RRMSE)為29.29%；曾偉生等[11]基于機載LiDAR數(shù)據(jù)，采用線性和非線性參數(shù)回歸方法對東北林區(qū)進行蓄積量反演，其中非線性回歸模型R2為0.71 ～0.82，略優(yōu)于線性回歸。已有多位研究者基于機載LiDAR數(shù)據(jù)進行林分蓄積量反演模型研究，但對于建模方法中參數(shù)回歸和非參數(shù)回歸模型的比較國內(nèi)研究較少。本研究以廣西國有高峰林場桉Eucalyptus樹人工林為研究對象，基于機載激光雷達數(shù)據(jù)及地面調(diào)查數(shù)據(jù)，采用逐步回歸、偏最小二乘回歸等參數(shù)回歸和隨機森林、支持向量回歸等非參數(shù)回歸進行蓄積量反演模型研建，并通過模型評價指標對以上4種方法進行模型評估，進而選擇出擬合優(yōu)度、泛化能力最優(yōu)模型。

1 數(shù)據(jù)獲取與預處理

1.1 研究區(qū)概況

研究區(qū)位于廣西壯族自治區(qū)南寧市興寧區(qū)的國有高峰林場，22°51′ ～23°02′N，108°06′ ～108°31′E，該區(qū)地處亞熱帶地區(qū)，年平均氣溫約21 ℃，年平均降水量為1 200 ～1 500 mm，相對濕度為79%，屬丘陵地貌，海拔為100 ～460 m，坡度為6° ～35°，具有較厚的赤紅壤，適宜亞熱帶和熱帶樹種生長，森林覆蓋率達87%，主要樹種為杉木Cunninghamia lanceolata、巨尾桉Eucalyptus grandis×E. urophylla、馬尾松Pinus massoniana等。

1.2 樣地數(shù)據(jù)

調(diào)查時間為2018年1—2月，研究區(qū)內(nèi)共設置71塊桉樹樣地(圖1)，其中激光雷達覆蓋范圍內(nèi)共57塊樣地。樣地大小為20 m×20 m、25 m×25 m和25 m×50 m。采用實時動態(tài)差分技術(shù)(real-time Kinematic，RTK)進行樣地定位，記錄樣地中心點及樣地角點。采用每木檢尺的方法，使用胸徑尺、激光測高儀和皮尺逐一測量樣地內(nèi)樹木的胸徑、樹高等數(shù)據(jù)。統(tǒng)計樣地調(diào)查數(shù)據(jù)得到樣地算數(shù)平均樹高、算數(shù)平均胸徑、樣地面積(表1)。通過廣西地區(qū)桉樹二元材積表對單木材積量進行計算[12]，進而計算得到樣地尺度的公頃蓄積量值 (V樣地)：V樣地=V公頃S樣地/10 000。其中：V公頃為通過二元材積公式計算得的每公頃蓄積量，S樣地為桉樹樣地面積。

圖1 研究區(qū)樣地分布示意圖Figure 1 Location of the sample plots distribution

1.3 LiDAR數(shù)據(jù)及預處理

于2018年1月采集機載LiDAR數(shù)據(jù)，使用有人機并搭載LMS-Q680i激光雷達掃描儀，實際飛行高度為1 000 m，最大掃描角度為 30°，波長為1 550 nm，激光脈沖長度為3 ns，采樣間隔為1 ns，最大掃描頻率為400 KHz，垂直分辨率為0.15 m，點云密度為 10 pt·m?2。

使用LiDAR 360軟件對機載激光雷達點云進行點云拼接、去噪、地面點分類、基于地面點歸一化等預處理，結(jié)果如圖2所示。根據(jù)樣地角點以中心點坐標對點云進行裁剪，提取出基于樣地尺度特征變量共48個，包括37個點云高度參數(shù)、10個密度參數(shù)、郁閉度等點云特征變量，具體參數(shù)見表2。

圖2 樣地點云預處理結(jié)果示例Figure 2 Example of point cloud preprocessing results of sample plots

表2 提取點云特征變量Table 2 Extracting point cloud feature variables

1.4 建模數(shù)據(jù)準備

使用MATLAB激光雷達覆蓋范圍內(nèi)57個樣本進行隨機抽樣，按照3∶1的比例選取42個作為訓練樣本，15個作為驗證樣本。以樣地實測公頃蓄積量為因變量，篩選后點云特征為自變量，采用逐步回歸、偏最小二乘回歸、隨機森林回歸、支持向量回歸模型進行擬合 。

2 研究方法

2.1 變量篩選

2.1.1 基于逐步篩選法優(yōu)選變量 對于逐步回歸，采用逐步篩選法對所有特征變量進行變量篩選；再對篩選出特征變量進行多重共線性檢驗，計算各變量間容忍度或特征變量間方差膨脹因子(FVI) 2個統(tǒng)計量[13]，對于容忍度≤0.2或FVI≥5的變量進行進一步討論，確定最優(yōu)特征變量子集。FVI=1/(1?R2)。其中：R為特征變量間相關(guān)系數(shù)。

2.1.2 基于隨機森林優(yōu)選變量 對于偏最小二乘回歸、支持向量機回歸、隨機森林回歸等方法，使用隨機森林中重要性排序?qū)μ卣髯兞窟M行優(yōu)選。其主要原理為隨機森林算法在構(gòu)建各決策樹時，對某一特征變量進行取舍，若此時均方誤差(mean square error, MSE)有較大變化，則該特征變量重要性高，最終得到所有特征變量重要性排序。

2.2 蓄積量模型建立方法

2.2.1 逐步回歸逐步回歸 (stepwise regression,SR)可用于篩選并剔除引起多重共線性的變量，逐步回歸建立模型一般形式為：Y=β0+β1X1+β2X2+···+βiXi+ε。逐步篩選法結(jié)合了向前選擇變量法和向后選擇變量法的優(yōu)點，對i個自變量X分別與因變量Y建立一元回歸模型，計算各變量所對應F值，其中，β0為常數(shù)，βi為回歸系數(shù)，ε為誤差項。在建立逐步回歸模型時，選擇當前未加入模型的預測變量中F的最大值所對應的Xi加入模型，再對已選入預測變量逐個進行t檢驗，若存在已選入預測變量不再顯著，則將其剔除。重復以上步驟，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止，將保留下來的最優(yōu)子集作為建立模型的特征變量，并將候選變量對應系數(shù)保留。

2.2.2 偏最小二乘回歸 偏最小二乘回歸(partial least squares regression, PLSR)結(jié)合了主成分分析、典型相關(guān)分析的優(yōu)點，實現(xiàn)了數(shù)據(jù)降維、信息綜合[14]。首先對訓練樣本進行數(shù)據(jù)標準化處理，調(diào)用plsregress函數(shù)提取主成分與自變量和因變量組合。分別計算自變量和因變量提取出成分的貢獻率進而計算得到累計貢獻率，當度量因子Q2h≥0.097 5停止主成分的提取[15]。統(tǒng)計主成分對個數(shù)并建立各自變量與因變量的線性表達式，最后根據(jù)所建立的各主成分對模型，整理得到PLSR模型。

2.2.3 支持向量機回歸 支持向量機算法(support vector machine, SVM)利用內(nèi)積核函數(shù)代替高維空間的非線性映射, 不涉及大數(shù)定律和概率測度等有關(guān)問題, 且SVM的決策函數(shù)僅有少數(shù)支持向量決定，該算法在解決小樣本多維度回歸和分類問題時泛化能力強，具備較好的“魯棒性”。VAPNIK[16]在SVM分類的基礎上引入了不敏感損失函數(shù)，得到了支持向量機回歸算法(support vector regression, SVR)。

使用Libsvm工具箱實現(xiàn)SVR算法，采用網(wǎng)格搜索法(grid search)對常用的4種核函數(shù)進行參數(shù)尋優(yōu)[17]，即線性、多項式、RBF (徑向基核函數(shù))、sigmoid (多層感知機核函數(shù))。同時進行十折交叉驗證，保證回歸模型中懲罰系數(shù)(C)與gamma值(g)達到最優(yōu)。

2.2.4 隨機森林回歸 隨機森林回歸(RFR)算法采用自助采樣法(bootstrap sampling)，在以決策樹構(gòu)建Bagging集成的基礎上，對樣本和特征變量進行隨機選擇[18]。調(diào)用TreeBagger函數(shù)進行RFR算法建模，通過對決策樹數(shù)量(ntree)和最小葉子點數(shù)進行參數(shù)尋優(yōu)，直到袋外(out-of-bag，OOB)誤差的MSE達到最小，保證模型預測性能達到最優(yōu)，將尋優(yōu)結(jié)果作為RFR模型的建模參數(shù)，用于模型構(gòu)建。

2.3 模型評估方法

本研究使用決定系數(shù)(R2)、均方根誤差(RMSE)和平均絕對誤差(MAE)對林分蓄積量估測模型進行定量的精度驗證和評價。其計算公式為：

3 結(jié)果與分析

3.1 變量篩選結(jié)果

3.1.1 逐步篩選法 在采用逐步篩選法進行偏F檢驗時，取偏F檢驗拒絕域的臨界值為F進＞F出，F(xiàn)進為選入變量時的臨界值、F出為刪除變量時的臨界值。本研究設定F進為0.10，F(xiàn)出為0.11，篩選特征變量結(jié)果為D9、HP95、Hmax、Hkurtosis。

對篩選出的變量進行雙變量相關(guān)性分析，根據(jù)相關(guān)系數(shù)(R)計算各變量間方差膨脹因子(FVI)。計算結(jié)果如表3，Hmax與HP95之間FVI為 250.25，F(xiàn)VI大于5，認為兩者之間存在共線性。因此分別以D9、Hmax、Hkurtosis和D9、HP95、Hkurtos為建模因子。

表3 各特征變量間的方差膨脹因子Table 3 FVI calculation of each variable

3.1.2 基于隨機森林篩選變量 對所有的48個特征變量進行重要性排序，調(diào)用隨機森林中OOBPermuted VarDeltaError參數(shù)，得到所有特征變量的重要性。選擇重要性大于0.2的變量作為建模因子，篩選結(jié)果如圖3，包括8個高度變量、2個密度變量，其中Hcurt重要性最高，達0.52；D6重要性相對最低，為0.22。

圖3 點云特征變量重要性排序Figure 3 Importance ranking of point cloud characteristic variables

3.2 模型建立與評價

3.2.1 逐步回歸法 根據(jù)特征變量優(yōu)選結(jié)果，分別建立以D9、Hmax、Hkurtosis和D9、HP95、Hkurtos為建模因子的多元線性模型：

使用SPSS 26對Y1、Y2模型進行初步評估，計算2種模型的相關(guān)系數(shù)(R)、R2以及標準估算誤差(SE)，結(jié)果如表4。結(jié)果表明Y2中R2、SE均優(yōu)于Y1，因此選定Y2作為逐步回歸模型。

表4 逐步回歸模型初步評估Table 4 Preliminary evaluation of stepwise regression model

3.2.2 偏最小二乘回歸 利用PLSR提取的主成分F1、F2與因變量Y關(guān)系分別為F1=5.840Y+b和F2=0.271Y+b，b為常數(shù)，得到最小二乘回歸模型為：

3.2.3 支持向量機回歸 調(diào)用meshgrid函數(shù)對C、gamma進行參數(shù)尋優(yōu)并采用十折交叉驗證，調(diào)用svmtrain函數(shù)，分別構(gòu)建4種不同核函數(shù)的SVR模型并進行模型訓練，訓練結(jié)果如表5。得到擬合結(jié)果最優(yōu)模型為RBF-SVR，其C為8，gamma為0.125。RBF-SVR 模型R2為 0.85，RMSE 為 29.24 m3·hm?2，MAE為94.98 m3·hm?2，選定該模型作為本研究的SVR模型。

表5 SVR不同核函數(shù)擬合結(jié)果Table 5 Fitting results of different kernel functions of SVR

3.2.4 隨機森林回歸 本研究采用窮舉法對決策樹數(shù)量和最小葉子點數(shù)進行參數(shù)尋優(yōu)。分別設置決策樹數(shù)量和最小葉子點數(shù)最小值為50和2，每次遞增10和1，對尋優(yōu)過程循環(huán)，通過觀察OOB誤差的MSE變化，直到尋找到本模型最優(yōu)參數(shù)，尋優(yōu)結(jié)果如圖4。最終確定決策樹數(shù)量為90，葉子節(jié)點數(shù)量為2，將尋優(yōu)結(jié)果作為建模參數(shù)，輸入訓練樣本，完成RFR模型構(gòu)建。

圖4 隨機森林參數(shù)尋優(yōu)結(jié)果Figure 4 Evaluation of growing stock volume inversion model

3.2.5 模型評價 使用驗證樣本對各模型的預測性能進行評估，結(jié)果如表6、圖5。可以看出不同的回歸模型均表現(xiàn)出較好的擬合結(jié)果，其中擬合結(jié)果最優(yōu)模型為RFR模型，模型評價結(jié)果R2為0.95，RMSE 為 12.64 m3·hm?2，MAE 為 8.00 m3·hm?2，RBF-SVR 模型其次，R2為 0.94，RMSE 為 13.09 m3·hm?2，MAE為 11.65 m3·hm?2。將驗證樣本帶入模型，檢驗模型的預測能力，結(jié)果顯示各模型預測能力與模型的擬合效果一致。通過傳統(tǒng)方法的劃分訓練樣本與驗證樣本，存在一定的偶然性與不確定性，為了進一步確保各模型的穩(wěn)定性與泛化能力，采用留一法交叉驗證(leave-one-out cross validation，LOOCV)對本研究中各模型進行再次評估[19]，如圖6所示：RBF-SVR模型表現(xiàn)最優(yōu)，R2為0.88，RMSE為21.35 m3·hm?2，MAE 為 16.62 m3·hm?2，與其他模型相比，R2高出 0.03 ～0.07，RMSE 減少 2.58 ～5.17 m3·hm?2，MAE 減少 0.79 ～3.32 m3·hm?2。

表6 蓄積量反演模型評估Table 6 Evaluation of growing stock volume inversion model

圖5 各模型訓練結(jié)果與驗證結(jié)果散點圖Figure 5 Scatter diagram of training results and verification results of each model

圖6 各模型十折交叉驗證結(jié)果Figure 6 Ten-fold cross validation results of each model

4 討論

本研究采用SR、PLSR等參數(shù)回歸方法與RFR、SVR等非參數(shù)回歸方法進行廣西高峰林場桉樹人工林的蓄積量反演模型研建，通過對點云特征變量篩選、建模因子優(yōu)選、參數(shù)尋優(yōu)等方法保證各模型達到最優(yōu)性能，在此基礎上采用留一法對各模型進行交叉驗證，保證了模型的穩(wěn)定性與泛化能力，并將該結(jié)果作為本研究各模型性能評估的最終結(jié)果。研究結(jié)果如下：①核函數(shù)為RBF的SVR在4種模型中預測精度最高 (ΔR2為 0.03 ～0.07、ΔRMSE 為 2.58 ～5.17 m3·hm?2、ΔMAE 為 0.79 ～3.32 m3·hm?2)，采用 SVR 模型在解決林業(yè)激光雷達領域的回歸預測問題已有廣泛應用，肖越[20]結(jié)合Landsat 8、高分2號(GF-2)共2種光學遙感數(shù)據(jù)與極化SAR數(shù)據(jù)進行旺業(yè)甸林場尺度的蓄積量反演，采用SVR模型擬合精度最高(R2為0.48，RMSE 為 57.27 m3·hm?2)，與 MLR、RFR 等方法相比，ΔR2為 0.10 ～0.11、ΔRMSE 為 5.30 ～6.00 m3·hm?2。趙勛等[21]基于機載激光雷達進行廣西高峰林場林分平均樹高估測中，通過隨機森林特征變量篩選，采用SVR相比較RFR方法，R2增加0.01、RMSE減少0.06 m。以上研究結(jié)果與本研究一致，進一步表明解決基于LiDAR點云特征變量的林分蓄積量反演問題時，結(jié)合隨機森林篩選特征變量與支持向量機回歸可作為有效建模方法。

②本研究采用2種方法對特征變量進行篩選，逐步回歸法保證入選的特征變量顯著且各變量之間不存在共線關(guān)系，隨機森林則通過計算各特征變量對建模時的貢獻度排序進而篩選變量。肖越[20]利用逐步回歸法從遙感特征變量中篩選出建模因子用于蓄積量回歸模型構(gòu)建；周蓉等[22]在基于Landsat 8 遙感影像反演地上生物量的研究中，采用隨機森林重要性排序方法從遙感因子中篩選出特征變量用于構(gòu)建模型。這2種方法均廣泛應用于構(gòu)建模型前的變量篩選，合理地選擇模型篩選方法，更有利于回歸模型構(gòu)建。

③本研究旨在選取最優(yōu)桉樹蓄積量反演模型，參數(shù)方法SR中逐步篩選法已對變量進行篩選，故SR模型不再采用隨機森林進行特征變量篩選。在本研究模型評估中，SR模型評估結(jié)果均優(yōu)于PLSR模型，袁鈺娜等[23]對東北林區(qū)4種不同針葉林蓄積量反演所建立的回歸模型中，PLSR模型的擬合結(jié)果和預測精度均優(yōu)于SR (ΔR2為0.05 ～0.15，ΔRMSE為2.6% ～4.2%)，與本研究結(jié)果不一致，但由于2種回歸方法并未使用同一變量篩選方法，因此，不能直接認為在解決預測回歸問題時前者更加可靠，對于以上2種回歸方法的模型性能與預測能力需進一步討論。

④本研究在林場尺度通過較少樣本進行桉樹林分蓄積量反演，非參數(shù)方法總體預測精度較高且RBF-SVR模型略優(yōu)于RFR模型。SVR核心思想為將低維空間中的向量用非線性函數(shù)映射到一個高維特征空間，進而尋求線性回歸超平面并解決低維空間中的非線性問題[24]，支持向量機回歸解決多維度小樣本回歸預測問題相較隨機森林回歸更有優(yōu)勢。在反演尺度較大且訓練樣本較多的林分蓄積量反演問題時，采用支持向量機回歸方法訓練效率將明顯下降，此時選擇隨機森林回歸將更適用。

⑤本研究采用逐步篩選法篩選出的4個特征變量中，點云密度變量被選擇1次；在采用隨機森林篩選的10個變量中，點云密度變量被選擇2次，且密度變量D9在重要性排序中位于第3位。點云密度變量作為本研究中重要解釋變量之一，該類型變量能夠描述林分中樹木的水平結(jié)構(gòu)信息，與樹高(垂直結(jié)構(gòu)變量)等參數(shù)結(jié)合可解釋林分空間結(jié)構(gòu)信息。孫忠秋等[25]加入林分郁閉度作為水平結(jié)構(gòu)的解釋變量，與僅用點云高度參數(shù)相比較，RMSE下降0.27 m3·hm?2、MAE下降0.08 m3·hm?2，該變量對模型預測性能有一定提升。本研究將郁閉度作為候選變量，但在逐步篩選法中又將其剔除，因為在基于隨機森林法的重要性排序中，林分郁閉度參數(shù)在48個候選變量中排名第44位。由此可見，在林業(yè)激光雷達應用中，將林分密度用于衡量林分水平結(jié)構(gòu)信息更加可靠。

5 結(jié)論

本研究通過2種參數(shù)回歸方法(逐步回歸、偏最小二乘回歸)與2種非參數(shù)方法(隨機森林回歸、支持向量機回歸)對林分蓄積量進行估測。在采用傳統(tǒng)方法劃分訓練樣本與驗證樣本的結(jié)果中，非參數(shù)方法模型精度均優(yōu)于參數(shù)回歸方法；采用留一法對各模型預測精度進行交叉驗證時，表現(xiàn)最優(yōu)模型RBFSVR屬于非參數(shù)回歸方法，表明解決林業(yè)激光雷達領域中的回歸預測問題時，非參數(shù)方法相較參數(shù)方法更有優(yōu)勢。本研究中蓄積量反演模型已在廣西高峰林場內(nèi)取得較好模型精度，但對于大尺度反演森林蓄積量的適用性有待進一步論證。