李洋 彭道黎 袁鈺娜

(森林資源和環(huán)境管理國家林草局重點實驗室(北京林業(yè)大學),北京,100083)(珠江水利委員會珠江流域水土保持監(jiān)測中心站)

森林生物量是評價森林生態(tài)系統(tǒng)固碳和碳平衡能力的重要變量,準確估算森林生物量對于研究大面積陸地生態(tài)系統(tǒng)的碳循環(huán)尤為重要[1]。森林地上生物量占陸地生態(tài)系統(tǒng)總量的90%[2],作為森林結(jié)構參數(shù)的重要組成部分,間接地反映森林的固碳能力,在森林碳儲量評估中充當著重要的評價因子。因此,快速、準確地估測森林地上生物量已經(jīng)成為全球氣候變化研究領域日益關注的問題。

國內(nèi)外對森林地上生物量反演進行了大量的研究,傳統(tǒng)的生物量估計通常是根據(jù)固定樣地連續(xù)觀測數(shù)據(jù),通過異速生長方程預測樣地內(nèi)單木生物量求和得到樣地生物量,從而進一步得到更大尺度的森林生物量[3]。但是這種測量方法存在較大的局限性,不僅耗時多、成本高、工作量大,而且對生態(tài)系統(tǒng)具有破壞性。隨著遙感與地理信息技術的快速發(fā)展,利用主被動遙感影像技術(波段信息、植被指數(shù)等)與森林地上生物量建立參數(shù)化及非參數(shù)化模型研究越來越多[4-7]。然而光學被動遙感影像技術易受森林植被物候、天氣等影響,并且存在植被信號飽和問題。隨著激光雷達技術的不斷發(fā)展及日趨成熟,森林參數(shù)的定量獲取取得了顯著突破。

激光雷達依據(jù)地面采樣點激光回波脈沖相對于發(fā)射激光主波之間的時間延遲獲得傳感器到地面采樣點的距離[8],獲取點云高度、密度、分布、強度及波形信息,從而得到樣地的三維結(jié)構參數(shù),并且具有極高的測距分辨能力和抗干擾能力。因此,相比于其他遙感技術,激光雷達在林木及林分的垂直、水平信息獲取上更加準確便捷。目前,激光雷達技術已被廣泛應用于測量植被高度[9]、葉面積指數(shù)[10]、郁閉度[11]、地上生物量[12]等林分特征參數(shù)[13]。

20世紀80年代中期,激光雷達技術開始應用于估測森林生物量的研究。激光雷達反演生物量的估測方法主要分為傳統(tǒng)的統(tǒng)計回歸模型(簡單線性回歸、多元逐步回歸等)和機器學習模型(隨機森林、神經(jīng)網(wǎng)絡、支持向量機等)。其中傳統(tǒng)的統(tǒng)計回歸方法常常需滿足一定的假設前提,但森林生長數(shù)據(jù)的連續(xù)觀測和層次性難以滿足以上假設。隨著機器學習技術的快速發(fā)展,為森林生長收獲預估提供了新的途徑。機器學習算法能夠深度挖掘數(shù)據(jù)中的有效信息,通過快速處理自變量與因變量的復雜關系建立預測模型,已經(jīng)廣泛應用于遙感和生態(tài)領域[14]。李旺等[15]基于激光雷達數(shù)據(jù)及實測單木結(jié)構信息,分別從樣地尺度和單木尺度對森林地上生物量進行估算,結(jié)果表明兩種尺度下模型的預測結(jié)果與地面實測值都具有明顯的相關性;龐勇等[16]采用高度變量和密度變量,基于多元線性回歸方法對3種不同森林類型的生物量進行估測,結(jié)果表明3種森林類型估測的決定系數(shù)(R2)均在0.8以上。雖然統(tǒng)計回歸方法模型公式較為簡單、便于計算,但機器學習算法能更好的發(fā)揮遙感大數(shù)據(jù)量的優(yōu)勢,捕捉生物量與各種特征變量之間的復雜非線性關系[17]。Gleason et al.[18]基于激光雷達數(shù)據(jù)的高度變量,比較了線性混合效應回歸、隨機森林、支持向量回歸和Cubist模型,結(jié)果表明在樣地水平上支持向量回歸的生物量模型精度最高。Cao et al.[19]采用兩種遙感數(shù)據(jù)源,結(jié)合5種機器學習算法建立生物量估測模型,結(jié)果發(fā)現(xiàn)RF算法的模型最優(yōu)。由此可見,機器學習方法在林業(yè)科學領域的研究有廣闊的發(fā)展前景。

近年來,作為一種新興的深度機器學習優(yōu)化算法,極限梯度提升(XGBoost)算法[20]能適應復雜的非線性關系,模型具有更強的并行處理能力,可以有效解決在機器學習回歸模型中出現(xiàn)的過擬合問題。黃宇玲等[21]基于逐步回歸的XGBoost方法建立森林蓄積量的估測模型;張亦然等[22]基于XGBoost算法進行了草甸地上生物量的反演研究,該算法都表現(xiàn)出更好的效果。然而,該方法在森林地上生物量的反演建模中的研究尚不充分,并且對于XGBoost模型的解釋性并未分析。因此,本文利用機載激光雷達數(shù)據(jù)和樣地實測數(shù)據(jù),通過多元逐步回歸、隨機森林、支持向量機和XGBoost算法分別對根河林業(yè)局開拉氣林場研究區(qū)森林地上生物量進行估測比較,并利用樣本中每個特征所分配到的沙普利加和解釋(SHAP)值對XGBoost模型的可解釋性進行分析,進一步驗證XGBoost算法在森林地上生物量模型構建中的可行性和適用性。

1 研究區(qū)概況

研究區(qū)位于地處內(nèi)蒙古呼倫貝爾盟東北部開拉氣林場,隸屬根河市林業(yè)局(120°12′～122°55′E,50°20′～52°30′N)(見圖1)。地處大興安嶺西坡中段,海拔800～1 000 m。該地區(qū)屬寒溫帶濕潤型森林氣候,并具有大陸季風性氣候的特征。年降水量400～550 mm,降水多集中在7—8月份,年平均氣溫-5.3 ℃,結(jié)凍期210 d以上,境內(nèi)遍布永凍層,個別地段30 cm以下即為永凍層。研究區(qū)主要的優(yōu)勢樹種為興安落葉松(Larixgmelinii),伴生樹種有白樺(Betulaplatyphylla)和山楊(Populusdavidiana)等。

2 研究方法

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

研究區(qū)共建立了125塊圓形樣地,每塊樣地半徑為13.82 m,面積約為600 m2,其中包括針闊混交林69塊、闊葉林30塊、針葉林26塊。在每個樣地的中心點設立明顯標志,并采用差分定位方法,準確定位樣地中心位置,保證樣地橫縱坐標的定位精度達到1 m以內(nèi),以確保地面樣地數(shù)據(jù)與激光雷達數(shù)據(jù)的匹配,準確提取遙感信息。地面樣地調(diào)查時間為2019年9—11月,調(diào)查內(nèi)容包括3大類:(1)林分因子(優(yōu)勢樹種組、起源、齡組、郁閉度等);(2)林木因子(樣木編號、立木類型、樹種、胸徑、樹高、方位角、水平距等);(3)其他信息(樣地號、位置坐標、海拔、坡度、坡位、坡向、調(diào)查員、調(diào)查日期等)。

根據(jù)樣地調(diào)查數(shù)據(jù)應用異速生長方程[23]對樣地內(nèi)每株單木的生物量計算,求和得到樣地尺度內(nèi)森林地上生物量,進而換算成單位面積的地上生物量。樣地數(shù)據(jù)基本信息及研究區(qū)主要樹種異速生長方程見表1、表2。

表1 調(diào)查樣地數(shù)據(jù)基本統(tǒng)計結(jié)果

表2 調(diào)查樣地樹種異速生長方程

2.2 機載LiDAR數(shù)據(jù)

機載激光雷達數(shù)據(jù)在2019年9月6日—10月16日航攝獲取,傳感器型號為RIEGL VUX-1UAV機載激光掃描儀,掃描儀最大脈沖發(fā)射頻率為550 kHz,光束發(fā)散角為0.5 mrad,光斑直徑為50 mm,平均點密度約4點·m-2,平均地面點距離約1 m,測量精度10 mm;搭載平臺為中型旋翼無人機,機長2.2 m,翼長3.3 m,飛行高度約100 m,平飛速度70～110 km·h-1。

圖1 根河研究區(qū)地面樣地點分布

本研究應用國產(chǎn)LiDAR360軟件對原始點云數(shù)據(jù)進行處理。首先,根據(jù)航跡線和激光發(fā)射器搭載平臺的參數(shù)信息對LiDAR點云進行航帶平差和去噪預處理;其次,對去噪處理后的點云數(shù)據(jù)進行分類,分為地面點(0 m)、低矮植被點(0～1 m)、中等植被點(1～20 m)、高植被點(>20 m)4類;再次,通過不規(guī)則三角網(wǎng)算法(TIN)結(jié)合地面點插值生成數(shù)值高程模型(DEM)[24],對點云進行高度歸一化處理,去除地形影響;最后,根據(jù)地面樣地中心點和地面樣地半徑對歸一化點云進行裁切,得到與125塊地面樣地對應的LiDAR點云集。

2.3 特征變量提取和篩選

根據(jù)激光雷達提取的數(shù)據(jù),特征空間可以分為點云和波形兩種,其中點云特征數(shù)據(jù)可以根據(jù)點云高度、密度、強度特征等直接進行分析,計算相關統(tǒng)計或指數(shù)特征量(不同的高度分位數(shù)、平均值、標準差、峰度、偏度等[4,25])。在林業(yè)研究中,高度變量是指點云高程值相關的統(tǒng)計參數(shù),直接反映了樹木的垂直結(jié)構和高度信息;密度變量反映了樹木點云的返回密度;強度變量則體現(xiàn)了樹木點云返回激光雷達傳感器的能量大小。這些指標已經(jīng)被廣泛應用于林業(yè)領域中單木及林分結(jié)構參數(shù)的估算[26-27]。

在本研究中,使用國產(chǎn)LiDAR360軟件對高度歸一化的點云數(shù)據(jù)集進行處理,提取得到樣地尺度的LiDAR點云特征變量。主要選取包括主要的16個高度變量、5個密度變量和12個強度變量等共33個特征變量[28]作為建模自變量(見表3)。

表3 激光雷達特征變量

續(xù)(表3)

許多變量都可以參與構建生物量模型,但這些變量之間往往存在高度相關性,并且與生物量之間的相關性不高[29]。因此,并非所有的變量都可用于生物量模型的構建。為了避免數(shù)據(jù)冗余,降低模型復雜度,提高模型精度,需要對原始特征變量數(shù)據(jù)集進行篩選。

在本研究中,對于多元線性回歸統(tǒng)計模型,采用皮爾森相關系數(shù)對選取的自變量進行篩選,提取與生物量相關極顯著的特征因子,再通過逐步回歸方法對每一個引入的自變量進行篩選剔除。對于機器學習模型,采用遞歸特征消除算法(RFE)進行特征選擇[30]。RFE算法是一個遞歸過程,它根據(jù)不同的重要性度量對特征變量進行排序。其核心思想是:在每次迭代中都會評估特征的重要性并刪除非重要特征,最終得到用于構建模型的最優(yōu)特征子集[31]。常用的RFE算法有支持向量機-遞歸特征消除算法(SVM-RFE)和隨機森林-遞歸特征消除算法(RF-RFE)[32-33]。最后,分別利用多元線性回歸(MLR)、隨機森林(RF)、支持向量機(SVM)和極限梯度提升(XGBoost)算法構建樣地生物量估測模型(算法均通過R4.1.1軟件實現(xiàn))。

2.4 生物量模型建立

多元逐步回歸模型:多元逐步回歸是一種是以向前引入為主,變量可進可出的采用多個自變量對因變量進行預測的統(tǒng)計方法。綜合了向前引入法和向后剔除法的優(yōu)點,可以修正多重共線性,廣泛應用于各領域的回歸模型[34]。公式如下:y=β0+β1x1+β2x2+…+βnxn+ε。式中:y為因變量(生物量);xi為自變量(激光雷達特征變量),βi為模型參數(shù),i=0、1、2、…、n;ε為誤差項。

激光雷達特征變量能間接反映出生物量的變化,當多個自變量和因變量形成線性關系時,即建立生物量與激光雷達變量的多元線性回歸模型。

隨機森林模型:隨機森林(RF)是一種統(tǒng)計學習理論,它是利用自舉(bootsrap)重抽樣方法從原始樣本中抽取多個樣本,對每個自舉(bootsrap)樣本進行決策樹建模,然后組合多棵決策樹的預測,通過投票得出最終預測結(jié)果[35]。在整個采樣過程中,有些樣本可能會被多次采樣,而有些訓練數(shù)據(jù)不會被采樣,這部分訓練數(shù)據(jù)稱為“袋外”(OOB)數(shù)據(jù);OOB數(shù)據(jù)不參與模型擬合過程,但用于檢查模型的泛化。由于隨機性可以有效地降低模型方差,RF算法可以實現(xiàn)良好的泛化性和低方差抗性,而無需對決策樹進行額外的“修剪”[36]。利用R軟件中的隨機森林包構建隨機森林模型。

支持向量機模型:支持向量機(SVM)算法的基本思想是將樣本空間通過非線性映射到一個高維特征空間,使在特征空間中可以應用線性學習機的方法解決樣本空間中的高度非線性分類和回歸等問題[37]。支持向量回歸利用結(jié)構風險最小化的同時優(yōu)化模型的性能和泛化能力,并且能夠找到非線性和唯一的解決方法[38]。本研究應用R軟件中的e1071包構建SVR模型,核函數(shù)選擇徑向基核函數(shù)。

極限梯度提升模型:極限梯度提升(XGBoost)算法是在傳統(tǒng)的推進(Boosting)算法基礎上引入正則化項,進一步控制模型的復雜度,也是梯度提升(Gradient Boosting)算法的實現(xiàn),具有運算速度快、效果好、易于調(diào)整參數(shù)和海量數(shù)據(jù)處理等優(yōu)勢[20]。并且相比于其他機器學習算法,具有更強的可解釋性。XGBoost的核心算法思想是在逐步添加樹到模型的過程中,使得模型整體的預測效果提升。其目標是要使得樹群的預測值盡量接近真實值,而且有盡量大的泛化能力。XGBoost目標函數(shù)分為損失函數(shù)和正則化項,損失函數(shù)揭示訓練誤差(即預測值和實測值的誤差),正則化定義復雜度,避免過擬合。應用R軟件中的極限梯度提升包建立模型。

2.5 模型驗證

本研究采用以下評價指標:均方根誤差(RMSE)、決定系數(shù)(R2)、平均絕對偏差(MAE)。R2值越接近于1,模型的擬合精度越好,均方根誤差和平均絕對偏差的值越小,模型的預測精度越高。3個評價指標計算公式如下:

2.6 機器學習模型超參數(shù)優(yōu)選

算法模型的性能由其超參數(shù)決定,超參數(shù)設定是構建模型的關鍵。相同的算法,不同的參數(shù),得到的結(jié)果具有差異性,參數(shù)調(diào)整對模型性能有很大影響,優(yōu)化算法參數(shù)能顯著提高模型預測精度[39]。采用網(wǎng)格搜索技術對3種機器學習算法進行超參數(shù)調(diào)優(yōu),對訓練集數(shù)據(jù)進行5次10倍交叉驗證法得到穩(wěn)定的模型結(jié)果。利用R軟件包對SVR、RF和XGBoost算法進行超參數(shù)調(diào)優(yōu)(見表4)。

3 結(jié)果與分析

3.1 變量篩選

3.1.1 皮爾森相關系數(shù)及逐步回歸的特征變量選擇

由表5可知,應用SPSS19.0軟件,對提取的33個激光雷達特征變量和地面樣地的森林地上生物量進行皮爾森相關性分析,得到與地上生物量極顯著相關(P<0.01)的自變量特征因子共計27個,其中包括16個高度特征變量、5個點云密度特征變量和6個點云強度特征變量。使用R4.1.1軟件包,在皮爾森相關性分析的基礎上,再通過逐步回歸方法,基于最小信息統(tǒng)計量(AIC)為準則進行自變量再優(yōu)選,最后得到多元線性回歸建模的最佳自變量均為高度變量共4個(包括點云高度平均值、高度偏度、點云高度的二次冪平均、點云高度的三次冪平均)。篩選結(jié)果剔除了點云密度變量和強度變量,僅保留點云高度變量,樣地的點云高度變量與地上生物量存在極顯著的相關性(P<0.01),說明點云高度變量是回歸分析最佳建模因子。

表4 3種機器學習算法的超參數(shù)調(diào)整范圍

表5 自變量因子與樣地地上生物量的Pearson相關系數(shù)

3.1.2 遞歸特征消除法的特征變量選擇

由圖2可知,在自變量達到4個時,RMSE的值最小。按其重要值排序依次為:點云高度的二次冪平均、點云高度平均值、點云高度的三次冪平均、高度垂直分布的25%分位數(shù)。遞歸特征消除特征變量的篩選結(jié)果與Pearson及逐步回歸有三個自變量重合。

圖2 遞歸特征消除方法篩選特征變量

3.2 模型精度評價

由表6可知,按照3∶1的比例隨機劃分為訓練集(93塊)和測試集(32塊),4種模型在訓練集和測試集進行擬合,不同模型的訓練集和測試集擬合精度不同。在訓練集上,RF模型擬合效果最好(RMSE=9.98 t·hm-2,R2=0.93,MAE=5.69 t·hm-2),其次為XGBoost模型(RMSE=10.80 t·hm-2,R2=0.89,MAE=7.24 t·hm-2),而MLR模型(RMSE=15.92 t·hm-2,R2=0.81,MAE=10.58 t·hm-2)和SVM模型(RMSE=16.49 t·hm-2,R2=0.81,MAE=10.15 t·hm-2)擬合精度相近,4種模型在訓練集上的表現(xiàn)都很好,R2都在0.8以上,但這并不能說明模型的泛化能力。在測試集上,XGBoost模型的擬合效果最好(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2),其次是SVR模型(RMSE=12.88 t·hm-2,R2=0.69,MAE=9.31 t·hm-2),而MLR模型的精度相對最低(RMSE=13.99 t·hm-2,R2=0.68,MAE=10.21 t·hm-2)。綜合上看,XGBoost模型的兩種數(shù)據(jù)集擬合精度都高且RMSE、R2、MAE差距最小,并且在測試集上明顯優(yōu)于MLR、RF和SVR模型,具有最佳的穩(wěn)定性和泛化能力。

表6 不同模型的精度評價結(jié)果

3.3 XGBoost重要值分析及模型解釋

同樣應用R軟件中的xgboost包功能,進一步分析經(jīng)遞歸特征消除算法(RFE)篩選后的4種特征變量與模型生物量的相關性和重要性,并利用SHAP值(樣本中每個特征所分配到的數(shù)值)解釋XGBoost模型。SHAP在2017年提出,用于解釋XGBoost等“黑箱”模型[40],其基本定義是對于每個預測樣本,模型都產(chǎn)生一個預測值,SHAP值是該樣本中每個特征所分配到的數(shù)值。相比于傳統(tǒng)的特征值重要性排序,SHAP值可以進一步反映出每個樣本中特征變量的影響力及正負性[41]。

圖3是SHAP特征密度散點圖,圖上的每個點都是一個特征和一個樣本的SHAP值,該值代表了這個特征對單個預測的貢獻,點的集合代表了特征整體上對預測結(jié)果影響的大小和趨勢。y軸上的位置從高到低由特征重要性決定,x軸上的位置由SHAP值決定,顏色從淺到深代表特征值從小到大,該圖結(jié)合特征重要性對影響地上生物量的因素進行分析。由圖3可知,參與建模的4種變量重要值排序依次為:點云高度的平均值、點云高度的三次冪平均、點云高度的二次冪平均、點云高度垂直分布的25%分位數(shù)的重要值從大到小依次為11.684、7.058、2.284、1.324。SHAP給出的影響地上生物量的最重要特征為點云高度的平均值。并且該特征與地上生物量呈正相關關系,即隨點云高度平均高的增加,估測的林分地上生物量越大?？傮w上看,4種特征變量與生物量基本都呈正相關關系。

Hm為點云高度平均值;Hsq為點云高度的二次冪平均;Htq為點云高度的三次冪平均;H25為點云高度垂直的25%分位數(shù)。

4 結(jié)論與討論

本研究以機載LIDAR數(shù)據(jù)和林分樣地調(diào)查數(shù)據(jù)作為數(shù)據(jù)源,單位地上生物量為研究對象,根據(jù)LIDAR數(shù)據(jù)源提取的特征變量,分別采用MLR、RF、SVR和XGBoost共兩類4種算法建立研究區(qū)的森林地上生物量反演模型,探討了XGBoost算法在森林地上生物量反演模型中的適用能力。

總體而言,XGBoost模型和MLR、RF、SVR模型相比,XGBoost在訓練集(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2)和測試集(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2)的精度評價指標接近,綜合表現(xiàn)最優(yōu),而其他三種模型在測試集的評價指標相較于訓練集都存在一定程度上的降低,說明XGBoost模型具有更強的泛化能力。張亦然等[22]在利用XGBoost算法建立草地地上生物量模型時,發(fā)現(xiàn)整體上利用XGBoost算法要優(yōu)于MLR和RF算法建立的模型。Li et al.[42]使用Landsat 8和Sentinel-A影像估算湖南省亞熱帶森林生物量的研究結(jié)果表明,在3種遙感數(shù)據(jù)集下XGBoost模型表現(xiàn)均優(yōu)于MLR和RF模型。而且與以往其他學者在該地區(qū)的地上生物量反演的研究成果相比[43],本研究的XGBoost模型精度要明顯優(yōu)于其基于k-NN算法的研究結(jié)果。因此,根據(jù)XGBoost算法構建的模型精度最高,具有更好的穩(wěn)定性和準確性。

在篩選建模因子方面,初始的LIDAR特征變量集包含高度、密度和強度等3種類型變量。而兩種不同篩選方法都同時剔除了密度變量和強度變量,只保留高度變量參與模型建立,說明剔除的兩種變量與生物量的相關性不強,不適于參與構建生物量模型。這是由于變量易受發(fā)射功率、范圍、入射角、環(huán)境參數(shù)和目標結(jié)構特征的影響[44],導致在不同情況下同一特征的獲取值差異較大,難以反映目標的真實特征。例如,機載LiDAR的回波信號不易穿透高郁閉度林分,從而獲取的點云數(shù)據(jù)多分布于冠層表面[45]。另外,篩選得到的高度變量都包含點云高度的二次冪平均(Hsq)、點云高度平均值(Hm)和點云高度的三次冪平均(Htq),證明三種高度變量都與生物量有顯著的相關性,點云高度能夠很好地反映林分的平均高度信息[27]。

在模型解釋性分析方面,本研究基于SHAP值對XGBoost算法構建生物量模型結(jié)果進行簡單分析。4種特征變量都對地上生物量有積極影響。但由于研究中樣本數(shù)較少,數(shù)據(jù)類型單一,尚未顯示出SHAP模型對于多源數(shù)據(jù)和復雜模型關系的解釋分析能力。隨著多元遙感數(shù)據(jù)融合與機器學習算法在林業(yè)領域中的應用越來越廣泛,可以采用SHAP模型對復雜變量及“黑箱”模型進行解釋,為林業(yè)數(shù)據(jù)挖掘與相關分析提供新途徑。