李曉偉， 吳保國*， 蘇曉慧， 陳玉玲， 彭意欽， 于永輝， 范小虎

(1.北京林業(yè)大學(xué)信息學(xué)院,林業(yè)信息化研究所,北京 100083; 2.廣西壯族自治區(qū)國有高峰林場，南寧 530000)

林分生產(chǎn)力的準(zhǔn)確預(yù)測對森林經(jīng)營及其方案編制意義重大，是實現(xiàn)人工林可持續(xù)經(jīng)營的基礎(chǔ)[1]。林分蓄積是衡量小班林分生產(chǎn)力的重要指標(biāo)，國內(nèi)外學(xué)者對蓄積預(yù)估模型已有大量研究，構(gòu)建林分蓄積模型的方法很多，主要為傳統(tǒng)統(tǒng)計建模方法和機(jī)器學(xué)習(xí)方法。傳統(tǒng)方法包括多元非線性回歸方程[2]、借助樹高和斷面積等中間變量的聯(lián)立方程組法[3]、混合效應(yīng)模型法[4-5]。傳統(tǒng)方法是在滿足數(shù)據(jù)獨(dú)立、正態(tài)分布和方差齊性等假設(shè)前提下進(jìn)行的，但是由于森林?jǐn)?shù)據(jù)的固有變異性，上述假設(shè)通常難以滿足[6]。此外，林分蓄積還受林分密度及經(jīng)營水平等因素影響，林分生長系統(tǒng)具有非線性、復(fù)雜性本質(zhì)特征，使得傳統(tǒng)的數(shù)學(xué)公式模型難以精準(zhǔn)表述[6]。而機(jī)器學(xué)習(xí)方法可對非線性關(guān)系進(jìn)行良好模擬，不僅在預(yù)測精度上有明顯優(yōu)勢，而且應(yīng)用方便[7]，因此，構(gòu)建機(jī)器學(xué)習(xí)模型進(jìn)行林分蓄積量預(yù)估的研究也越來越多。目前，機(jī)器學(xué)習(xí)方法中應(yīng)用較多的是人工神經(jīng)網(wǎng)絡(luò)[8-12]。對于非神經(jīng)網(wǎng)絡(luò)方法，林卓等[13]用支持向量機(jī)構(gòu)建福建西北地區(qū)杉木人工林蓄積模型，高若楠[14]使用小班數(shù)據(jù)采用隨機(jī)森林方法對東北天然闊葉林蓄積生長量進(jìn)行預(yù)估。決策樹機(jī)器學(xué)習(xí)方法在預(yù)測方面表現(xiàn)更好的解釋性，并且通過集成學(xué)習(xí)方式進(jìn)行整合，既可以提高預(yù)測精度又可以減少產(chǎn)生過擬合的幾率，在生態(tài)研究等領(lǐng)域中應(yīng)用較多[15-16]。但以小班數(shù)據(jù)為基礎(chǔ)，采用多種決策樹為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)方法對考慮立地條件和密度因子的桉樹(EucalyptusrobustaSmith)人工林進(jìn)行蓄積預(yù)測的研究比較欠缺。本研究以廣西國有高峰林場速生桉為研究對象，構(gòu)建多種以決策樹為基礎(chǔ)結(jié)構(gòu)器的非集成、集成學(xué)習(xí)模型，結(jié)合年齡、立地條件、密度因素研究蓄積預(yù)估。為采用蓄積預(yù)估模型作為生產(chǎn)力判斷的泛化研究進(jìn)行一定探索，并結(jié)合具體研究地區(qū)為林場速生桉的造林決策提供支持。

1 材料與方法

1.1 研究區(qū)概況

研究地點(diǎn)為廣西壯族自治區(qū)國營高峰林場，該林場為廣西最大國有林場，位于22°49′N～23°15′N，108°08′E～108°53′E，屬南亞熱帶氣候，夏長冬短，光熱充足，雨量充沛，年平均氣溫在21 ℃左右，極端最高溫40 ℃，最低溫-2 ℃，積溫為7 500 ℃左右，年降雨量為1 200～1 500 mm，多集中在每年6—9月[17-18]。地貌主要由山丘和丘陵構(gòu)成。海拔為150～400 m，坡度為20～30°。土壤以赤紅壤為主，質(zhì)地為中壤土或輕粘土。森林植被主要有桉樹、馬尾松(PinusmassonianaLamb)、杉木(CunninghamialanceolataHook)等樹種，其中桉樹類全部為速生桉樹種，主要包括尾葉桉(EucalyptusurophyllaS.T.Blake)、巨尾桉(Eucalyptusgrandis×urophylla)、尾巨桉(Eucalyptusurophylla×E.grandis)。

1.2 研究數(shù)據(jù)

本研究以廣西壯族自治區(qū)桉樹人工林為研究對象，選擇2008年高峰林場森林資源規(guī)劃設(shè)計調(diào)查(簡稱二類調(diào)查)中優(yōu)勢樹種為桉樹的3 500個小班以及2010年廣西國家森林資源連續(xù)清查(簡稱一類清查)的200個桉樹固定樣地兩個數(shù)據(jù)源作為研究數(shù)據(jù)。參照《森林資源規(guī)劃設(shè)計調(diào)查技術(shù)規(guī)程》[19]與《森林資源數(shù)據(jù)采集技術(shù)規(guī)范》[20]確定二類調(diào)查、一類清查各因子。提取小班和固定樣地數(shù)據(jù)中共有的立地因子(海拔高度、坡向、坡位、坡度、枯枝落葉厚度、腐殖質(zhì)層厚度、土壤種類)和林分因子(林齡、公頃蓄積、公頃株數(shù)、平均胸徑、平均樹高)進(jìn)行整理，其中二類調(diào)查數(shù)據(jù)用作模型訓(xùn)練與驗證，一類清查數(shù)據(jù)用作模型泛化測試。用于此建模數(shù)據(jù)的因子統(tǒng)計量見表1與表2。

表1 建模數(shù)據(jù)林分因子統(tǒng)計

表2 建模數(shù)據(jù)立地因子統(tǒng)計

1.3 模型構(gòu)建

集成學(xué)習(xí)主要分為并行與串行兩種，本文選擇決策樹集成學(xué)習(xí)構(gòu)建2種并行算法(Bagging和隨機(jī)森林)和2種串行算法(梯度提升決策樹和XGboost)來構(gòu)建桉樹小班蓄積模型。Bagging(bootstrap aggregating)算法[21]采用隨機(jī)有放回的選擇訓(xùn)練數(shù)據(jù)。隨機(jī)森林[22]在以決策樹為基學(xué)習(xí)器構(gòu)建Bagging集成的基礎(chǔ)上，進(jìn)一步在決策樹的訓(xùn)練過程中引入了隨機(jī)屬性(如小班數(shù)據(jù)中坡度、土壤厚度等字段)選擇。相比Bagging，隨機(jī)森林在當(dāng)前結(jié)點(diǎn)的屬性集合中選擇一個包含m個屬性的集合。隨機(jī)森林的基學(xué)習(xí)器，分別采用CART模型與Ctree模型。梯度提升決策樹GBDT(gradient boosting decision tree)算法[23]基于回歸樹模型，采用迭代的方法最小化損失函數(shù)，進(jìn)而得到最優(yōu)解。GBDT中每一棵回歸樹結(jié)點(diǎn)劃分都是基于之前多個回歸樹的預(yù)測結(jié)果、采用梯度迭代法訓(xùn)練新弱分類器，使得新回歸樹的預(yù)測結(jié)果與實際值之間的損失函數(shù)達(dá)到最小?；貧w問題損失函數(shù)設(shè)定為均方差函數(shù)，其梯度結(jié)果與殘差公式一致?；貧w樹劃分節(jié)點(diǎn)后樣本的取值采用平均值法。XGboost[24]方法作為GBDT的改進(jìn)，基于梯度提升部分目標(biāo)函數(shù)為損失函數(shù)加正則項，損失函數(shù)為訓(xùn)練誤差，正則項為訓(xùn)練樹復(fù)雜度，通過正則項控制模型復(fù)雜程度，模型更為精準(zhǔn)。

模型構(gòu)建分為數(shù)據(jù)清洗與預(yù)處理、特征選擇與數(shù)據(jù)集劃分和訓(xùn)練模型構(gòu)建與調(diào)優(yōu)三步，具體流程如圖1所示。采用統(tǒng)計軟件R 3.5.4調(diào)用相關(guān)包進(jìn)行模型訓(xùn)練、驗證與結(jié)果分析。

圖1 模型構(gòu)建流程

Step1：數(shù)據(jù)清洗與預(yù)處理。數(shù)據(jù)預(yù)處理按順序包括干擾數(shù)據(jù)清除、定性因子量化處理與數(shù)據(jù)標(biāo)準(zhǔn)化處理。廣西壯族自治區(qū)速生桉樹的輪伐期為5～6 a，因此，首先剔除數(shù)據(jù)中樹木平均年齡大于 6 a 的數(shù)據(jù)，然后將土壤種類、坡位、坡向定性因子啞變量化[25]轉(zhuǎn)化為定量因子，最后，由于min-max標(biāo)準(zhǔn)化會受極端值影響，對數(shù)據(jù)進(jìn)行Z-score標(biāo)準(zhǔn)化處理，以消除量綱帶來的影響。

Step2：特征選擇與數(shù)據(jù)集劃分。林分生長的主要影響因子包括立地質(zhì)量、林分平均年齡、林分密度，因此選擇2個林分因子(年齡和公頃株數(shù)，表1)和8個立地因子(表2)作為輸入變量，公頃蓄積量為輸出變量。數(shù)據(jù)集劃分方式為70%訓(xùn)練數(shù)據(jù)集，30%測試數(shù)據(jù)集，對于驗證數(shù)據(jù)集部分，不單獨(dú)劃分，而是統(tǒng)一采用十折交叉驗證的方法3次重復(fù)進(jìn)行模型驗證。

Step3：參數(shù)單值預(yù)訓(xùn)練、參數(shù)全范圍模型訓(xùn)練與調(diào)優(yōu)。采用任一訓(xùn)練模型，初步選擇參數(shù)取值范圍中間值作為參數(shù)進(jìn)行模型預(yù)訓(xùn)練，將其結(jié)果做理論性可用性驗證。在理論可行性驗證成功基礎(chǔ)上進(jìn)行參數(shù)全范圍訓(xùn)練模型構(gòu)建與調(diào)優(yōu)，選擇非集成與集成兩種方式共9種模型，調(diào)用caret軟件包的train函數(shù)進(jìn)行模型訓(xùn)練。模型構(gòu)建過程三個方面核心參數(shù)、參數(shù)訓(xùn)練范圍見表3。

表3 模型構(gòu)建與調(diào)優(yōu)參數(shù)統(tǒng)計表

Step4：模型泛化預(yù)估測試。選擇Step3對比結(jié)果最優(yōu)模型，以2015年廣西壯族自治區(qū)部分一類清查桉樹數(shù)據(jù)進(jìn)行泛化測試，判斷其對廣西所有地區(qū)桉樹蓄積預(yù)估是否具有泛化性。對于訓(xùn)練好的最優(yōu)模型，首先用R語言封裝模型成為函數(shù)，之后使利用Java語言構(gòu)建工程，輸入一類清查數(shù)據(jù)文件，選擇Step2中特征選擇的10個自變量，使用Rserve模式調(diào)用并執(zhí)行模型程序包文件，計算樣地蓄積預(yù)測值[7]，與觀測值進(jìn)行線性回歸分析計算R2與P值，其結(jié)果作為泛化測試結(jié)果。

1.4 理論性可用性驗證

蓄積預(yù)估屬于生長收獲模型，變化應(yīng)符合生物S型生長曲線規(guī)律。隨機(jī)選擇一個小班，采用單變量法固定1.3 Step2中8個立地因子以及公頃株數(shù)，以時間為自變量，蓄積預(yù)估量為因變量，以1.3 Step3預(yù)訓(xùn)練模型作生長曲線圖，評價模型是否符合樹木生長規(guī)律，進(jìn)而判斷模型的理論可用性。若理論可用性通過，進(jìn)入?yún)?shù)全范圍的模型訓(xùn)練與調(diào)優(yōu)、模型評價階段。若理論可用性不通過，回到數(shù)據(jù)預(yù)處理部分對針對業(yè)務(wù)問題對數(shù)據(jù)再處理。

1.5 模型評價

對于數(shù)字評價，交叉驗證模型評價方法采用3個指標(biāo)進(jìn)行評價和檢驗,分別是決定系數(shù)(R2)、均方根誤差(RMSE)以及平均絕對偏差(MAE)。

R2越接近1表示模型擬合效果越好，MAE與RMSE相對越低表示模型誤差越小。對于模型訓(xùn)練集效果，選擇R2與RMSE，對于模型測試集效果選擇RMSE與MAE。

2 結(jié)果與分析

2.1 理論可用性分析

預(yù)訓(xùn)練和修正訓(xùn)練結(jié)果得到圖2所示。發(fā)現(xiàn)第5 a蓄積量相較第4 a增長量較小，整體不符合S型生長曲線，因此返回數(shù)據(jù)預(yù)處理階段進(jìn)行數(shù)據(jù)再處理。分析原因，該地區(qū)桉樹主伐年齡為 5 a，在二類調(diào)查時對n<年齡

圖2 理論可用性預(yù)訓(xùn)練與修正訓(xùn)練結(jié)果

2.2 模型結(jié)果與對比分析

表4與表5分別顯示9個模型以最優(yōu)R2為標(biāo)準(zhǔn)的最優(yōu)參數(shù)組合以及模型在訓(xùn)練集與測試集中的相關(guān)精度指標(biāo)。

表4 模型訓(xùn)練最優(yōu)參數(shù)組合

由表5結(jié)果可知，相同方法、不同基學(xué)習(xí)器的模型訓(xùn)練結(jié)果不同。其中單棵樹模型Ctree基學(xué)習(xí)器相比CART基學(xué)習(xí)器在訓(xùn)練集上的R2與RMSE幾乎相同，但是在測試集上Ctree基學(xué)習(xí)器的R2與RMSE明顯降低。無論訓(xùn)練集還是測試集，并行集成學(xué)習(xí)模型R2與RMSE均明顯低于采用CART基學(xué)習(xí)器模型；而在隨機(jī)森林方法中實驗結(jié)果出現(xiàn)反轉(zhuǎn)，即CART基學(xué)習(xí)器模型各指標(biāo)在訓(xùn)練集與非測試集結(jié)果均好于Ctree基學(xué)習(xí)器模型，并且模型整體精度與誤差數(shù)字均為相對最低，模型整體效果為前六種中相對最優(yōu)。

表5 模型訓(xùn)練與測試結(jié)果

非集成學(xué)習(xí)與集成學(xué)習(xí)的模型評價量化指標(biāo)明顯不同，集成學(xué)習(xí)模型整體優(yōu)于非集成模型。無論訓(xùn)練集還是測試集，并行集成學(xué)習(xí)中采用相同基學(xué)習(xí)器的隨機(jī)森林模型相關(guān)指標(biāo)優(yōu)于一般并行Bagging方法。串行集成學(xué)習(xí)方法中模型評價指標(biāo)R2與RMSE最優(yōu)的依次為XGboost模型、增強(qiáng)回歸樹模型。對于串行與并行兩種集成學(xué)習(xí)方式，同時選擇CART作為基學(xué)習(xí)器，結(jié)合模型評價指標(biāo)，模型訓(xùn)練集結(jié)果由好到壞依次為XGboost、隨機(jī)森林、增強(qiáng)回歸樹與Cubist，模型測試集結(jié)果由好到次依次為XGboost、隨機(jī)森林、Cubist與增強(qiáng)回歸樹。

2.3 模型變量重要性評估與解釋

2.3.1變量重要性評估 選擇非集成決策樹類、并行集成類以及串行集成類中效果最優(yōu)的CART的單棵決策樹、隨機(jī)森林以及XGboost三種模型，計算得到各自變量相對于因變量每公頃蓄積的影響重要性。由圖3可知，對于立地、樹木以及密度三類自變量因子，三種模型的自變量重要性排名幾乎一致，特別是前5名重要性自變量三個模型相同，說明重要的幾個變量對于蓄積預(yù)估的影響程度不會因模型不同而改變，只是在貢獻(xiàn)度次序上有相對的改變。對于各自變量重要性排名，占比最高的為樹木因子中的年齡，其對于因變量重要性顯著高于其他因子，三種模型分別為達(dá)到86.5%、83.5%、78.0%。其次對于因變量存在影響的三個自變量在CART、隨機(jī)森林以及XGboost模型中的重要性數(shù)值排序分別為：土層厚度(3.0%)、腐殖質(zhì)層厚度(2.9%)、海拔(2.4%)；土層厚度(4.5%)、密度(3.0%)、海拔(2.4%)；海拔(4.9%)、土層厚度(3.8%)、密度(3.2%)。而坡向、坡位、坡度、土壤種類等因素重要性在三個模型中占比都低于1%，影響極小。

注：圖中編號1～10分別為海拔(m)、坡向、坡位、坡度(°)、枯枝落葉厚度(cm)、腐殖質(zhì)層厚度(cm)、土層厚度(cm)、土壤種類、年齡(a)、密度(hm-2)。

2.3.2模型解釋 相比于其他機(jī)器學(xué)習(xí)模型，決策樹具有一定的模型可解釋性，從而根據(jù)其樹杈分支的決策過程可以從中探究各自變量屬性在蓄積預(yù)估模型的劃分順序與劃分節(jié)點(diǎn)數(shù)字，從而對機(jī)器學(xué)習(xí)方法應(yīng)用于蓄積預(yù)估提供更好的解釋性。單棵樹模型相比集成學(xué)習(xí)方法，從實驗工程的角度可以還原決策的屬性劃分過程，具有更好的模型解釋性。屬性劃分為五級，第一、二級為年齡，第三級為土層厚度與枯枝落葉厚度，第四級為密度與海拔，層級越小劃分屬性影響越大。年齡值、土層與枯枝落葉厚度、密度與海拔值越大，蓄積越大。

2.4 模型泛化預(yù)估

模型調(diào)用結(jié)果如表6所示。由表6可知，桉樹在1～6 a蓄積平均值持續(xù)增長，且從第3 a開始蓄積增長值明顯擴(kuò)大，但是相對應(yīng)標(biāo)準(zhǔn)差也顯著地提升。殘差整體殘差分布較為均勻，不存在異方差情況(圖4)。模型經(jīng)檢驗后調(diào)整后R2為0.785，P值為2.2e-16，符合檢驗標(biāo)準(zhǔn)(圖5)。

表6 蓄積預(yù)估結(jié)果

圖4 XGboost模型在泛化測試數(shù)據(jù)集上殘差圖

圖5 XGboost模型在泛化測試數(shù)據(jù)集上值的散點(diǎn)圖

3 討論

本研究選擇二類調(diào)查數(shù)據(jù)用作模型擬合與驗證，選擇一類清查數(shù)據(jù)用作模型泛化。二類調(diào)查數(shù)據(jù)特點(diǎn)是數(shù)據(jù)量大，研究范圍相對集中，可以滿足該地區(qū)決策樹集成學(xué)習(xí)等機(jī)器學(xué)習(xí)模型對于數(shù)據(jù)訓(xùn)練量的基本要求，并且較多訓(xùn)練數(shù)據(jù)量可以在一定程度上增加模型擬合精度。泛化測試即模型對于未知數(shù)據(jù)預(yù)測的相對準(zhǔn)確性，兩個核心要求體現(xiàn)在待測試的未知數(shù)據(jù)不同于訓(xùn)練數(shù)據(jù)集所屬地區(qū)以防止過擬合以及未知數(shù)據(jù)精度有一定保證從而可以反饋泛化測試結(jié)果。一類清查數(shù)據(jù)精度相對較高、調(diào)查范圍大而分散的兩個特點(diǎn)符合模型泛化測試核心要求。

對于模型判定數(shù)字指標(biāo)，基于二調(diào)數(shù)據(jù)，XGboost在森林蓄積預(yù)估問題上有著最優(yōu)的模型效果,R2超過0.8，同比陳玉玲等[7]采用BP神經(jīng)網(wǎng)絡(luò)模型進(jìn)行華北落葉松蓄積預(yù)估精度有一定的提高。模型泛化應(yīng)用結(jié)果R2為0.785，P值小于0.000 1，符合檢驗，說明該模型精度達(dá)到應(yīng)用水平。本研究以海拔、坡向、坡位、坡度以及土壤相關(guān)等立地屬性來作為立地因素，以每公頃種植株數(shù)作為密度因素。林卓等[13]、王少杰等[5]分別使用斷面積、計算密度指數(shù)作為密度因素。僅考慮變量獨(dú)立性以及與其他自變量較小的交互性，株數(shù)相對更合適。受年齡影響的林分因子是否可以作為自變量引入模型仍有待討論，李宗俊[26]認(rèn)為可以先預(yù)估林分因子再代入模型進(jìn)行二次訓(xùn)練的方法也是參考選項但可能面臨自變量多重共線性等問題。

決策樹模型相比神經(jīng)網(wǎng)絡(luò)可以避免潛在的數(shù)據(jù)的過度擬合等特點(diǎn)[7]。集成學(xué)習(xí)方法精度顯著高于非集成方法，原因在于單一的決策樹模型具有一些不足，如模型的不穩(wěn)定性(數(shù)據(jù)中微小的變動可能會引起樹的巨大變化，從而影響解釋性)、次優(yōu)的預(yù)測能力等[27]。作為基學(xué)習(xí)器的模型，在集成學(xué)習(xí)方法上CART效果好，非集成學(xué)習(xí)上條件推斷樹更好。采用串行集成學(xué)習(xí)方法普遍好于一般并行集成學(xué)習(xí)，但是若并行集成考慮到隨機(jī)性而采用隨機(jī)森林的方式則與串行集成學(xué)習(xí)差異不明顯，該結(jié)果與Fernández-Delgado等[28]關(guān)于回歸問題采用的多種機(jī)器學(xué)習(xí)方法對比研究結(jié)果一致，原因是相比于條件推斷樹對于切分結(jié)果偏度的糾正，隨機(jī)森林在切分屬性的選擇上的隨機(jī)性改變更有效的減小了誤差，說明基于屬性劃分的選擇對于誤差的影響效果要大于切分后切分值的偏度修改。模型參數(shù)訓(xùn)練范圍與歐強(qiáng)新[19]的研究類似，但是最優(yōu)參數(shù)組合存在明顯不同，說明因變量、自變量屬性的不同會影響同類模型的訓(xùn)練情況。對于集成學(xué)習(xí)，本研究探討了并行與串行兩種方式，后續(xù)可以采用混合集成方式進(jìn)一步研究。