基于增強(qiáng)回歸樹(shù)的杉木人工林林分?jǐn)嗝娣e模型研究

臧 顥，黃錦程，劉洪生，歐陽(yáng)勛志，寧金魁*

（1.江西農(nóng)業(yè)大學(xué) 林學(xué)院，江西 南昌 330045；2.江西農(nóng)業(yè)大學(xué) 鄱陽(yáng)湖流域森林生態(tài)系統(tǒng)保護(hù)與修復(fù)國(guó)家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330045；3.江西省崇義縣林業(yè)局，江西 贛州 341300）

【研究意義】杉木（Cunninghamia lanceolata（Lamb.）Hook.）是我國(guó)最重要的用材樹(shù)種之一。據(jù)第9 次全國(guó)森林資源清查結(jié)果[1]，杉木人工林總面積達(dá)990.20 萬(wàn)hm2、蓄積為7.55 億m3，占全國(guó)人工林總面積的17.33%和總蓄積的22.30%。林分?jǐn)嗝娣e是林業(yè)調(diào)查中的主要因子以及林分生長(zhǎng)和收獲預(yù)估體系中的重要變量[2-3]，因而，研究杉木林分?jǐn)嗝娣e生長(zhǎng)對(duì)預(yù)估杉木林的收獲有重要意義?！厩叭搜芯窟M(jìn)展】為了構(gòu)建林分?jǐn)嗝娣e模型，一些學(xué)者采用多種統(tǒng)計(jì)方法進(jìn)行了研究?；诮魇?1 塊樣地?cái)?shù)據(jù)，杜紀(jì)山等[4]采用Richards 和Schumacher 方程作為基礎(chǔ)模型，以年齡、林分密度指標(biāo)和立地質(zhì)量指標(biāo)為自變量構(gòu)建了杉木林分?jǐn)嗝娣e生長(zhǎng)模型，并采用最小二乘法估計(jì)了模型參數(shù)。倪成才等[5]基于美國(guó)佐治亞州148 個(gè)火炬松人工林固定樣地的調(diào)查數(shù)據(jù)，采用代數(shù)差分方程構(gòu)建了林分?jǐn)嗝娣e模型，并采用最大似然和限制性最大似然估計(jì)了模型的參數(shù)。Zhao 等[6]收集了江西省、廣東省和廣西省共583 個(gè)杉木人工林固定樣地?cái)?shù)據(jù)，采用最小二乘法和3 水平混合效應(yīng)模型對(duì)含優(yōu)勢(shì)高和每公頃株數(shù)的林分?jǐn)嗝娣e生長(zhǎng)模型進(jìn)行了參數(shù)估計(jì)。Zhang 等[7]以江西省分宜縣杉木人工林為研究對(duì)象，采用貝葉斯模型平均法將林分模型、單木模型和直徑分布模型聯(lián)合起來(lái)，并對(duì)林分?jǐn)嗝娣e進(jìn)行了預(yù)測(cè)。張雄清等[8]采用廣義代數(shù)差分法構(gòu)建杉木林分?jǐn)嗝娣e生長(zhǎng)模型，并采用最小二乘法和貝葉斯法分別估計(jì)了模型參數(shù)?！颈狙芯壳腥朦c(diǎn)】盡管最小二乘法、混合效應(yīng)模型、貝葉斯等方法被廣泛用于構(gòu)建林分?jǐn)嗝娣e生長(zhǎng)模型，但這些統(tǒng)計(jì)方法在應(yīng)用時(shí)都存在著一定的假設(shè)前提，諸如參數(shù)的特定分布形式、自變量的獨(dú)立性等，且需要預(yù)先給定因變量和自變量之間的關(guān)系式[9-11]。而森林的生長(zhǎng)是一個(gè)復(fù)雜的過(guò)程，它受到立地條件[12-14]、氣候變化[15-16]、林分結(jié)構(gòu)[17-18]及其交互作用的綜合影響，因而上述假設(shè)通常難以滿足[19]。由于能處理復(fù)雜的變量關(guān)系，且不需要預(yù)設(shè)模型的形式和輸入數(shù)據(jù)的分布等優(yōu)點(diǎn)[20-21]，近些年來(lái)，機(jī)器學(xué)習(xí)被逐漸用于模擬森林生長(zhǎng)模型中，而增強(qiáng)回歸樹(shù)作為機(jī)器學(xué)習(xí)的一種，其預(yù)估精度較高且估計(jì)結(jié)果比較穩(wěn)定[22-24]，因而在林業(yè)上有著較多的應(yīng)用?！緮M解決的關(guān)鍵問(wèn)題】本研究以江西省崇義縣杉木人工林為研究對(duì)象，采用增強(qiáng)回歸樹(shù)構(gòu)建杉木人工林林分?jǐn)嗝娣e生長(zhǎng)模型，探討了增強(qiáng)回歸樹(shù)算法中不同參數(shù)的變化對(duì)模型的影響，并分析了不同林分密度指標(biāo)和不同立地指標(biāo)及其交互作用對(duì)林分?jǐn)嗝娣e的模擬影響。

1 研究區(qū)概況

研究區(qū)位于江西省贛州市崇義縣，屬于南嶺的北端。崇義縣的地理坐標(biāo)為113°55′E～114°38′E，25°24′N(xiāo)～25°55′N(xiāo)。全縣屬亞熱帶季風(fēng)氣候，年平均氣溫17.9 ℃，極端最高39.9 ℃，最低-8 ℃，年均降水量1 603.1 mm，無(wú)霜期307 d。全縣總面積2 206.27 km2，主要地貌為低山、丘陵，土壤類(lèi)型有黃壤、紅壤等。崇義縣屬常綠闊葉林生物氣候帶，主要樹(shù)種有絲栗栲（Castanopsis fargesiiFranch.）、米櫧（Castanopsis carlesii（Hemsl.）Hayata）、杉木（Cunninghamia lanceolata（Lamb.）Hook.）、馬尾松（Pinus massonianaLamb.）、苦楝（Melia azedarachL.）、南酸棗（Choerospondias axillaris（Roxb.）Burtt.et Hill）等。

2 材料與方法

2.1 數(shù)據(jù)來(lái)源

數(shù)據(jù)來(lái)源于2016 年江西省贛州市崇義縣的森林資源二類(lèi)調(diào)查數(shù)據(jù)，本研究選擇了起源為人工且優(yōu)勢(shì)樹(shù)種為杉木的喬木林小班進(jìn)行分析，共計(jì)8 780個(gè)小班，主要因子統(tǒng)計(jì)信息見(jiàn)表1。

表1 主要因子統(tǒng)計(jì)Tab.1 Summary statistics of primary factors

2.2 研究方法

2.2.1 確定自變量 基于前人的研究結(jié)果[25-27]，立地、密度和年齡是林分?jǐn)嗝娣e的主要影響因子。因此，本研究的自變量也按這3方面進(jìn)行選擇，其中，立地指標(biāo)選取了地位指數(shù)和地位級(jí)指數(shù)，密度指標(biāo)選取了每公頃株數(shù)和林分密度指數(shù)，年齡為平均年齡。地位級(jí)指數(shù)、地位指數(shù)和林分密度指數(shù)的計(jì)算如下：

（1）地位指數(shù)。依據(jù)蔡學(xué)林等[28]對(duì)贛南山地杉木人工林構(gòu)建的數(shù)量化地位指數(shù)模型和小班的地形因子信息，計(jì)算各小班的地位指數(shù)，基準(zhǔn)年齡為20 a。

（2）地位級(jí)指數(shù)。地位級(jí)指數(shù)是基準(zhǔn)年齡下的林分平均高，可采用差分方程進(jìn)行計(jì)算[29]。計(jì)算步驟為：1）采用理論生長(zhǎng)方程構(gòu)建林分平均高模型，如式①；2）假定立地條件的高低僅影響林分平均高所能達(dá)到的最大值，即方程中反映理論最大值的參數(shù)為自由參數(shù)，此時(shí)，年齡取基準(zhǔn)年齡則方程計(jì)算得到的就是地位級(jí)指數(shù)，如式②所示，進(jìn)而可采用差分方程對(duì)該參數(shù)進(jìn)行差分，最終形式如式③所示。

式中，Hi和ti分別表示第i個(gè)小班的林分平均高和年齡，SCIi表示第i個(gè)小班的地位級(jí)指數(shù)，t0為基準(zhǔn)年齡，本研究中取20年，β0為反映理論最大值的參數(shù)，βj為其他待估參數(shù)，β0i表示隨小班立地條件變化的自由參數(shù)，f（）表示理論方程，g（）表示對(duì)應(yīng)的差分方程。

本研究中選取了常見(jiàn)的理論生長(zhǎng)方程[30]，各方程形式如下：

式中，β1、β2為待估參數(shù)，其他符號(hào)如上所述。

（3）林分密度指數(shù)。依據(jù)Reineke[31]給出的公式計(jì)算各小班的林分密度指數(shù)，計(jì)算公式如下：

式中：SDIi、Ni和Dgi為第i個(gè)小班的林分密度指數(shù)、每公頃株數(shù)和林分平均胸徑，D0表示標(biāo)準(zhǔn)平均直徑，本研究中取25 cm[32]。

Pretzsch 等[33]通過(guò)對(duì)歐洲中部的4 種森林類(lèi)型進(jìn)行分析，發(fā)現(xiàn)樹(shù)木的徑向生長(zhǎng)會(huì)受到立地和密度的交互作用影響，這一影響隨樹(shù)種的變化而變化，因此，本研究在構(gòu)建林分?jǐn)嗝娣e模型時(shí)考慮了不同立地指標(biāo)和不同密度指標(biāo)的交互作用，并于不考慮交互作用的模型進(jìn)行對(duì)比，以判斷立地和密度的交互作用對(duì)杉木人工林林分?jǐn)嗝娣e的影響。通過(guò)8種自變量組合（如表2所示），分別模擬林分?jǐn)嗝娣e，構(gòu)成8種模型，以模型的模擬效果對(duì)立地指標(biāo)、密度指標(biāo)及交互作用的影響。

表2 不同的自變量組合Tab.2 Combinations of different independent variables

2.2.2 增強(qiáng)回歸樹(shù)構(gòu)建林分?jǐn)嗝娣e模型 增強(qiáng)回歸樹(shù)是一種基于集成學(xué)習(xí)的決策樹(shù)模型，通過(guò)對(duì)訓(xùn)練樣本進(jìn)行多次不放回隨機(jī)抽樣，生成一系列訓(xùn)練子樣本，并用回歸樹(shù)對(duì)每一個(gè)子樣本進(jìn)行擬合，擬合過(guò)程按順序進(jìn)行，將當(dāng)前的模型添加到之前的模型中，該過(guò)程一直持續(xù)到指定的迭代次數(shù)為止，最終的預(yù)測(cè)結(jié)果則依據(jù)各回歸樹(shù)擬合效果進(jìn)行加權(quán)[11,24]。

增強(qiáng)回歸樹(shù)建模有5 個(gè)輸入?yún)?shù)需要設(shè)定，分別為損失函數(shù)（distribution）、交互深度（interaction.depth）、收縮參數(shù)（shrinkage）、回歸樹(shù)數(shù)量（n.trees）、子取樣比例（bag.fraction）。針對(duì)回歸問(wèn)題，損失函數(shù)一般設(shè)為平方誤差。結(jié)合前人的研究結(jié)果[24,34]，子取樣比例取0.5。而交互深度、收縮參數(shù)和回歸樹(shù)數(shù)量的確定，不同的研究中無(wú)統(tǒng)一的結(jié)果[11,21,24]。為探討交互深度、收縮參數(shù)和回歸樹(shù)的數(shù)量的設(shè)置對(duì)擬合結(jié)果的影響，本研究中對(duì)這3個(gè)參數(shù)設(shè)置不同的數(shù)值，通過(guò)對(duì)比不同數(shù)值組合的模型評(píng)價(jià)指標(biāo)，對(duì)模型進(jìn)行優(yōu)化：交互深度取1、2、3、4、5、6、7、8、9、10，收縮參數(shù)取0.1、0.05、0.01、0.005、0.001，回歸樹(shù)的數(shù)量取100、200、300、400、500、600、700、800、900、1000。

2.2.3 模型評(píng)價(jià) 本研究采用5折交叉驗(yàn)證進(jìn)行模型評(píng)價(jià)：1）將數(shù)據(jù)隨機(jī)分成5個(gè)子樣本；2）每次用其中4個(gè)子樣本作為建模樣本，最后1個(gè)子樣本作為檢驗(yàn)樣本，依據(jù)檢驗(yàn)樣本計(jì)算評(píng)價(jià)指標(biāo)；3）重復(fù)上一步，直到5個(gè)子樣本都進(jìn)行了模型檢驗(yàn)；4）計(jì)算5次模型檢驗(yàn)得到的評(píng)價(jià)指標(biāo)的平均值，并以此判斷模型優(yōu)劣。

采用3種指標(biāo)評(píng)價(jià)模型的預(yù)估精度：MAB，平均絕對(duì)誤差；RMA，平均相對(duì)誤差絕對(duì)值；RMSE，均方根誤差，計(jì)算公式如下[21,35-36]：

式中，Gi和分別表示第i個(gè)小班的林分?jǐn)嗝娣e實(shí)測(cè)值和預(yù)測(cè)值，n為樣本數(shù)。

本研究中的所有計(jì)算均在R 軟件中完成，其中，增強(qiáng)回歸樹(shù)調(diào)用“gbm”包中的“gbm”函數(shù)實(shí)現(xiàn)的，增強(qiáng)回歸樹(shù)的參數(shù)優(yōu)化是在“caret”包中的“train”函數(shù)實(shí)現(xiàn)的。

3 結(jié)果與分析

3.1 地位級(jí)指數(shù)的計(jì)算

由表3可知，7種理論方程擬合的效果差異不大，其中，式⑩擬合的林分平均高效果最好。相較其他6 種理論方程，式⑩的MAB減少了0.06%～7.22%，RMSE減少了0.01%～3.91%，因此選擇式⑩基礎(chǔ)模型構(gòu)建差分方程，并計(jì)算地位級(jí)指數(shù)。計(jì)算地位級(jí)指數(shù)的差分形式如下：

式中各符號(hào)如前所述。

3.2 增強(qiáng)回歸樹(shù)參數(shù)優(yōu)化

圖1 顯示了8 種林分?jǐn)嗝娣e模型的參數(shù)優(yōu)化情況?？偟膩?lái)看，隨著回歸樹(shù)數(shù)量的增加，8 種模型的MAB、RMA和RMSE均呈現(xiàn)減小的趨勢(shì)。這一趨勢(shì)在不同交互深度下沒(méi)有明顯差異，但隨著收縮參數(shù)增加，這一趨勢(shì)逐漸減弱，當(dāng)收縮參數(shù)取0.05 和0.1 時(shí)，隨著回歸樹(shù)數(shù)量的增加，MAB、RMA和RMSE的變化無(wú)明顯趨勢(shì)。

綜合考慮MAB、RMA和RMSE，8 種模型的最優(yōu)參數(shù)值及其交叉驗(yàn)證的評(píng)價(jià)指標(biāo)見(jiàn)表4。可知：含有林分密度指數(shù)的模型明顯優(yōu)于含有每公頃株數(shù)的模型，含有林分密度指數(shù)的模型的MAB減少了80.05%～84.15%，RMA減少了80.42%～84.66%，RMSE減少了77.94%～82.27%；含有地位級(jí)指數(shù)的模型優(yōu)于含有地位指數(shù)的模型，含有地位級(jí)指數(shù)的模型的MAB減少了11.32%～25.84%，RMA減少了11.68%～22.96%，RMSE 減少了6.95%～25.24%；含有立地和密度交互作用的模型略優(yōu)于不含交互作用的模型，含交互作用的模型的MAB減少了4.56%～12.48%，RMA減少了5.88%～15.48%，RMSE減少了1.07%～9.34%，僅M3（地位指數(shù)與每公頃株數(shù)的無(wú)交互模型）優(yōu)于M7（地位指數(shù)與每公頃株數(shù)的有交互模型）。M6（考慮了地位級(jí)指數(shù)、林分密度指數(shù)及其交互作用）的各項(xiàng)評(píng)價(jià)指標(biāo)均為最小，可用于預(yù)估崇義縣杉木人工林林分?jǐn)嗝娣e。

表3 基于交叉驗(yàn)證的林分平均高模型評(píng)價(jià)Tab.3 The evaluation of stand mean height model using cross-validation

圖1 8種林分?jǐn)嗝娣e模型的增強(qiáng)回歸樹(shù)參數(shù)優(yōu)化Fig.1 Tuning boosted regression tree for 8 stand basal area models

表4 8種林分?jǐn)嗝娣e模型的參數(shù)優(yōu)化結(jié)果Tab.4 The tuning results for 8 stand basal area models

4 結(jié)論與討論

本研究以江西省崇義縣的杉木人工林為對(duì)象，采用增強(qiáng)回歸樹(shù)構(gòu)建了林分?jǐn)嗝娣e模型，通過(guò)交叉驗(yàn)證分析了不同立地指標(biāo)、不同密度指標(biāo)及其交互作用對(duì)模擬的影響。研究發(fā)現(xiàn)，經(jīng)過(guò)參數(shù)優(yōu)化的增強(qiáng)回歸樹(shù)能較好的擬合林分?jǐn)嗝娣e，最優(yōu)模型的MAB為0.504 7 m2/hm2，RMA為2.29%，RMSE為0.900 8 m2/hm2。

地位指數(shù)是基準(zhǔn)年齡時(shí)的優(yōu)勢(shì)高[37]，被廣泛應(yīng)用于林分?jǐn)嗝娣e模型中以反映立地質(zhì)量對(duì)林分?jǐn)嗝娣e的影響[4,9,27]。因?yàn)檫B續(xù)調(diào)查固定樣地?cái)?shù)據(jù)通常缺乏林分優(yōu)勢(shì)高數(shù)據(jù)，只有林分平均高數(shù)據(jù)[38-39]，因此有學(xué)者使用基于林分平均高計(jì)算的地位級(jí)指數(shù)評(píng)價(jià)立地質(zhì)量[29]。同時(shí)，也有學(xué)者研究發(fā)現(xiàn)林分平均高和優(yōu)勢(shì)高之間存在著相當(dāng)穩(wěn)定的數(shù)量關(guān)系[40-41]，因而使用基于林分平均高的地位級(jí)指數(shù)反映立地質(zhì)量是一種可行的方法。本研究對(duì)比了地位指數(shù)和地位級(jí)指數(shù)對(duì)林分?jǐn)嗝娣e模擬的效果，發(fā)現(xiàn)地位級(jí)指數(shù)對(duì)林分?jǐn)嗝娣e的模擬效果良好，這說(shuō)明用基于林分平均高的地位級(jí)指數(shù)能有效的反映立地質(zhì)量對(duì)林分?jǐn)嗝娣e的影響。此外，盡管本研究中，地位級(jí)指數(shù)對(duì)林分?jǐn)嗝娣e的模擬優(yōu)于地位指數(shù)，這可能是因?yàn)楸狙芯渴遣捎们叭藰?gòu)建的模型推算得到的地位指數(shù)，即存在著一定的誤差，進(jìn)而影響了地位指數(shù)對(duì)林分?jǐn)嗝娣e的模擬效果。

每公頃株數(shù)和林分密度指數(shù)是被廣泛應(yīng)用于森林生長(zhǎng)收獲模型中[17,42-43]，本研究對(duì)比了每公頃株數(shù)和林分密度指數(shù)對(duì)杉木人工林林分?jǐn)嗝娣e的模擬效果，結(jié)果顯示林分密度指數(shù)對(duì)林分?jǐn)嗝娣e的模擬效果明顯優(yōu)于每公頃株數(shù)，這與一些學(xué)者的研究結(jié)果相同[9,27]。盡管每公頃株數(shù)描述了單位面積上的林木株數(shù)，但卻沒(méi)有包含樹(shù)木大小信息，隨著樹(shù)木的徑向生長(zhǎng)，可能存在每公頃株數(shù)較少但由于樹(shù)木胸徑較大，導(dǎo)致樹(shù)木個(gè)體之間競(jìng)爭(zhēng)更強(qiáng)的情況，這在一定程度上也會(huì)影響到林分?jǐn)嗝娣e的生長(zhǎng)。林分密度指數(shù)則是將當(dāng)前的株數(shù)轉(zhuǎn)換到標(biāo)準(zhǔn)胸徑下的株數(shù)，即考慮了林分平均胸徑的影響，這可能是林分密度指數(shù)的模擬效果優(yōu)于每公頃株數(shù)的原因。

為了探討立地和競(jìng)爭(zhēng)是否存在對(duì)林分?jǐn)嗝娣e的交互作用，本研究對(duì)比了含交互作用項(xiàng)與不含交互作用項(xiàng)的林分?jǐn)嗝娣e模型的模擬效果，結(jié)果顯示，除含地位指數(shù)和每公頃株數(shù)的交互作用項(xiàng)的模型外，其余含交互作用項(xiàng)的林分?jǐn)嗝娣e模型擬合精度均略高于不含交互作用項(xiàng)的模型，這說(shuō)明了立地和林分密度的交互作用對(duì)杉木人工林林分?jǐn)嗝娣e有一定影響。這與Pretzsch 等[33]的研究發(fā)現(xiàn)類(lèi)似。因此，如果需要精確模擬森林的徑向生長(zhǎng)時(shí)，有必要考慮立地和林分密度的交互作用。

由于機(jī)器學(xué)習(xí)算法對(duì)模型的形式和輸入數(shù)據(jù)的分布沒(méi)有假設(shè)前提，且能處理復(fù)雜的變量關(guān)系、擬合效果好[11,20]，所以逐漸應(yīng)用于森林生長(zhǎng)收獲預(yù)估中。但也有學(xué)者認(rèn)為機(jī)器學(xué)習(xí)存在“黑箱”問(wèn)題，無(wú)法給出因變量和自變量的關(guān)系式，但其實(shí)每種機(jī)器學(xué)習(xí)算法存在各自的工作原理，而因變量和自變量的關(guān)系也可以通過(guò)模擬的方式給出，偏依賴圖也可以給出各自變量對(duì)因變量的重要性[21]。相對(duì)傳統(tǒng)的統(tǒng)計(jì)方法更關(guān)注于模型的基本形式（線性、非線性等）和變量的關(guān)系（是否存在共線性、自變量參數(shù)是否顯著）等問(wèn)題，采用機(jī)器學(xué)習(xí)方法構(gòu)建森林生長(zhǎng)模型時(shí)需要更多的考慮參數(shù)值的確定和過(guò)擬合問(wèn)題，參數(shù)值的確定可以通過(guò)對(duì)比不同參數(shù)取值組合的擬合效果來(lái)確定，但過(guò)于追求擬合效果也可能會(huì)導(dǎo)致過(guò)擬合問(wèn)題，即對(duì)建模樣本較好，但在檢驗(yàn)樣本表現(xiàn)較差的情況[11]，因此建議使用機(jī)器學(xué)習(xí)構(gòu)建模型時(shí)需要進(jìn)行模型檢驗(yàn)。

綜上所述，增強(qiáng)回歸樹(shù)作為一種機(jī)器學(xué)習(xí)算法，可應(yīng)用于林分?jǐn)嗝娣e模型的構(gòu)建，不同的交互深度、收縮參數(shù)及回歸樹(shù)數(shù)量對(duì)擬合精度影響較大，有必要進(jìn)行參數(shù)優(yōu)化，同時(shí)不同的自變量組合對(duì)擬合精度也有較大的影響，基于林分平均高的地位級(jí)指數(shù)和林分密度指數(shù)及其交互作用可作為模型自變量以提高模型的預(yù)測(cè)效果。