廣西杉木人工林多形地位指數(shù)模型構(gòu)建1)

張藝超 趙天忠 蘇曉慧

(北京林業(yè)大學(xué)，北京，100083)

森林的立地質(zhì)量影響著森林生產(chǎn)經(jīng)營(yíng)的多個(gè)方面，而立地質(zhì)量評(píng)價(jià)是實(shí)現(xiàn)科學(xué)造林以及經(jīng)營(yíng)森林的關(guān)鍵[1]。立地質(zhì)量評(píng)價(jià)有多種方法，地位指數(shù)法是以基準(zhǔn)年齡下優(yōu)勢(shì)木平均高作為衡量立地質(zhì)量的標(biāo)準(zhǔn)，優(yōu)勢(shì)木平均高受林分密度影響較小，而且相對(duì)易于獲取，因此應(yīng)用較為廣泛。

地位指數(shù)模型有單形和多形兩種。前者由導(dǎo)向曲線通過(guò)一定方法展開得到，所有地位指數(shù)下樹木生長(zhǎng)曲線形狀相同；后者不同的地位指數(shù)下生長(zhǎng)曲線形狀不同，即不同地位指數(shù)下生長(zhǎng)過(guò)程不同。多形地位指數(shù)模型由Trousdell K B, et al提出[2]，被證明更符合樹木生長(zhǎng)的客觀規(guī)律。段愛(ài)國(guó)等[3]采用差分法以6種生長(zhǎng)方程為基礎(chǔ)，探討了擬合多形地位指數(shù)方程最佳表達(dá)式；Dario Martín Benito, et al[4]應(yīng)用自適應(yīng)差分進(jìn)化算法，以Richards模型為基礎(chǔ)，建立了黑松的非線性混合模型；高光芹[5]等以神經(jīng)網(wǎng)絡(luò)模型構(gòu)建馬尾松多形地位指數(shù)模型；施恭明等[6]采用Korf生長(zhǎng)方程，通過(guò)改進(jìn)的單純形法，建立了馬尾松多形地位指數(shù)模型。

綜上所述，雖然國(guó)內(nèi)對(duì)多形地位指數(shù)模型已有較多研究，但多以現(xiàn)有的生長(zhǎng)方程為基礎(chǔ)，構(gòu)建參數(shù)模型進(jìn)行擬合，所得模型都需經(jīng)過(guò)迭代才能求出地位指數(shù)，應(yīng)用上存在不便，且針對(duì)廣西杉木人工林的多形地位指數(shù)模型研究還較少。本文以廣西高峰林場(chǎng)為研究區(qū)域，分別以Richards方程和隨機(jī)森林回歸算法為基礎(chǔ)構(gòu)建多形地位指數(shù)模型，并對(duì)兩個(gè)模型進(jìn)行對(duì)比分析，力求在不同的場(chǎng)景下應(yīng)用最適宜的模型，為廣西杉木人工林立地質(zhì)量評(píng)價(jià)提供參考。

1 研究區(qū)域概況

廣西國(guó)有高峰林場(chǎng)是自治區(qū)林業(yè)局直屬國(guó)有林場(chǎng)，林場(chǎng)場(chǎng)部位于南寧市興寧區(qū)，屬于大明山脈，地理坐標(biāo)為北緯22°49′～23°15′，東經(jīng)108°8′～108°53′。林場(chǎng)經(jīng)營(yíng)面積8.7萬(wàn)余hm2，森林蓄積量570多萬(wàn)m3，森林覆蓋率達(dá)87%,是廣西規(guī)模最大的國(guó)有林場(chǎng)。林場(chǎng)地貌多為丘陵與山地，分別占總面積的55.5%與38.7%。林場(chǎng)氣候?qū)儆趤啛釒Ъ撅L(fēng)氣候，年平均氣溫21 ℃，降水量1 200～1 500 mm[7]。廣西高峰林場(chǎng)植被主要以次生人工林為主，主要樹種為杉木(Cunninghamialanceolata)，巨尾桉(Eucalyptusgrandis×urophylla)，馬尾松(Pinusmassoniana)等。

2 研究方法

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

本文研究數(shù)據(jù)來(lái)自廣西壯族自治區(qū)高峰林場(chǎng)提供的2010年二類調(diào)查數(shù)據(jù)，其中以杉木為優(yōu)勢(shì)樹種的小班共有342個(gè)。本文應(yīng)用廣西高峰林場(chǎng)杉木樹種因子與解析木數(shù)據(jù)建立多形地位指數(shù)模型。小班樹種因子詳細(xì)情況如表1所示。

表1 小班樹種因子統(tǒng)計(jì)

2.2 模型的評(píng)價(jià)標(biāo)準(zhǔn)

本文分別基于Richards方程與隨機(jī)森林回歸算法進(jìn)行多形地位指數(shù)模型構(gòu)建，為了便于模型比較，將數(shù)據(jù)分為模型構(gòu)建數(shù)據(jù)集與測(cè)試檢驗(yàn)數(shù)據(jù)集，其中測(cè)試檢驗(yàn)數(shù)據(jù)集占總數(shù)據(jù)集的25%。模型的評(píng)價(jià)選用ERMS，EMA，R23個(gè)衡量指標(biāo)，評(píng)價(jià)指標(biāo)公式如式(1)、(2)、(3)所示。

均方根誤差(ERMS)：

(1)

平均絕對(duì)誤差(EMA)：

(2)

(3)

2.3 應(yīng)用Richards方程的多形地位指數(shù)模型

Richards方程為“S”型曲線方程，其基本形式如下：

H=a(1-e-bT)c。

(4)

式中:H為樹高(m)，T為林齡(a)，a、b、c為待定系數(shù)。

Richards方程為“S”型曲線，其拐點(diǎn)是可變的，相較于Logistic等固定拐點(diǎn)的方程更符合樹木的生長(zhǎng)特性，是目前應(yīng)用最廣泛的樹木生長(zhǎng)擬合方程之一。許多研究運(yùn)用Richards方程模擬樹木生長(zhǎng)，并取得較好的效果[3，8-11]，因此本文將以該模型為基礎(chǔ)進(jìn)行模型構(gòu)建。

傳統(tǒng)的同形地位指數(shù)模型是由生長(zhǎng)方程擬合的導(dǎo)向曲線展開得來(lái)的，不同地位指數(shù)的生長(zhǎng)曲線形狀相同，人為忽略了不同立地質(zhì)量下樹木生長(zhǎng)曲線不同這一客觀規(guī)律，因此多形地位指數(shù)模型優(yōu)于同形地位指數(shù)模型[12]。樹木生長(zhǎng)曲線形狀隨立地質(zhì)量變化，兩者之間滿足一定函數(shù)關(guān)系。立地質(zhì)量由地位指數(shù)(IS)來(lái)衡量，不同立地質(zhì)量曲線形狀由RT表示，則有關(guān)系函數(shù)：

RT=f(IS)。

(5)

本文以Richards方程來(lái)模擬樹木生長(zhǎng)，由式(4)可知，方程有3個(gè)未知參數(shù)a、b、c，這3個(gè)參數(shù)決定了方程所代表的地位指數(shù)及其生長(zhǎng)趨勢(shì)，結(jié)合式(5)可知此3個(gè)參數(shù)均與IS存在函數(shù)關(guān)系，設(shè)a=f1(IS)、b=f2(IS)、c=f3(IS)，可得：

H=f1(IS)(1-e-f2(IS)T)f3(IS)。

(6)

式(6)給出了不同立地質(zhì)量下樹高生長(zhǎng)模型的基本公式，找出關(guān)系f1，f2，f3即可得到多形地位指數(shù)模型。

2.4 應(yīng)用隨機(jī)森林回歸算法的地位指數(shù)模型

由2.3可知，樹高(H)，年齡(T)，地位指數(shù)(IS)之間是滿足某一函數(shù)關(guān)系的，而此關(guān)系為非線性關(guān)系。本小節(jié)將應(yīng)用集成學(xué)習(xí)方法中的隨機(jī)森林回歸算法來(lái)構(gòu)建此3個(gè)變量間的關(guān)系模型。

隨機(jī)森林回歸模型：隨機(jī)森林算法由Breiman于1995年提出[15]，該算法對(duì)Bagging算法進(jìn)行了改進(jìn)，利用Bootstrap方法從原始樣本中隨機(jī)抽取若干數(shù)量相等樣本，從所有輸入特征中選擇若干子特征集，對(duì)每個(gè)樣本使用分類與回歸樹(CART)作為弱學(xué)習(xí)器進(jìn)行建模。對(duì)于回歸算法，會(huì)取所有弱學(xué)習(xí)器結(jié)果的均值作為最終結(jié)果。

本次試驗(yàn)采用Scikit-learn方法庫(kù)實(shí)現(xiàn)隨機(jī)森林回歸算法。Scikit-learn是一個(gè)Python第3方提供的非常強(qiáng)力的機(jī)器學(xué)習(xí)庫(kù)，它包含了從數(shù)據(jù)預(yù)處理到訓(xùn)練模型的各個(gè)方面的算法實(shí)現(xiàn)。

模型的構(gòu)建：2.3中，最終得到形如H=f(IS，T)函數(shù)關(guān)系式(6)，當(dāng)需要求在一定地位指數(shù)下，不同年齡對(duì)應(yīng)的樹高時(shí)，可直接求得。但是式(6)并不能寫成IS=f(H，T)的形式，因此想要直接求得IS必須通過(guò)迭代求解。為了便于應(yīng)用，本文分別以(IS，T)為輸入變量，以H為輸出變量、以(H，T)為輸入變量，以IS為輸出變量進(jìn)行隨進(jìn)森林模型訓(xùn)練。

在Scikit-learn開發(fā)包中，采用sklearn. ensemble模塊下的Random Forest Regressor類實(shí)現(xiàn)隨機(jī)森林回歸算法。Random Forest Regressor類中有兩個(gè)主要參數(shù)：決策樹數(shù)目(n_estimators)、尋找最佳分割時(shí)需要考慮的特征數(shù)目(max_features)。決策樹數(shù)目越大越好，但相應(yīng)占用的內(nèi)存與訓(xùn)練和預(yù)測(cè)的時(shí)長(zhǎng)將增加，達(dá)到一定值后對(duì)訓(xùn)練效果提升將持續(xù)減弱，本次試驗(yàn)以測(cè)試集R2為衡量標(biāo)準(zhǔn)，由1開始依次測(cè)試發(fā)現(xiàn)決策樹數(shù)目接近300時(shí)，R2不再變動(dòng)，因此選取決策樹數(shù)目為300；本次試驗(yàn)訓(xùn)練數(shù)據(jù)特征維度較低，參數(shù)特征數(shù)目對(duì)結(jié)果影響較小，因此采用默認(rèn)值，即取特征總數(shù)的開方。

最后進(jìn)行模型訓(xùn)練，求出關(guān)系H=f(IS，T)與IS=f(H，T)的模型。

3 結(jié)果與分析

為了確定各小班地位指數(shù)，將小班平均樹高與該小班平均年齡下各解析木樹高進(jìn)行差值運(yùn)算，最終找到差值最小的解析木。以基準(zhǔn)年齡為20 a[16]，將基準(zhǔn)年齡下該解析木樹高作為小班的地位指數(shù)。重復(fù)以上步驟可以得到所有小班地位指數(shù)，為了減小誤差，如果該小班樹高與所有解析木樹高差值都大于1 m時(shí)，應(yīng)舍棄該小班不參與建模。

3.1 應(yīng)用Richards方程的模型擬合結(jié)果

將數(shù)據(jù)中的地位指數(shù)按級(jí)距1 m進(jìn)行分組，各組分別利用式(4)進(jìn)行非線性回歸，所用軟件為SPSS 25.0，各組所得未知參數(shù)值及決定系數(shù)如表2所示。

表2 各地位指數(shù)組參數(shù)擬合結(jié)果

利用SPSS軟件的曲線估計(jì)功能，選取幾種常用的基本初等函數(shù)來(lái)擬合關(guān)系f1，f2，f3，進(jìn)行量化探究，曲線估計(jì)結(jié)果如表3所示。

表3 參數(shù)擬合函數(shù)擬合決定系數(shù)

根據(jù)表2繪制散點(diǎn)圖，并綜合表3選取擬合函數(shù)。在滿足精度的前提下，為了降低模型復(fù)雜度，優(yōu)先選用形式更簡(jiǎn)單的基本函數(shù)。最終選取線性函數(shù)擬合f1，指數(shù)函數(shù)擬合f2，二次函數(shù)擬合f3。

將選定的參數(shù)擬合函數(shù)代入式(6)最終得出模型的最終式為：

(7)

式中：K1、K2、K3、K4、K5、K6、K7為7個(gè)待定系數(shù)，IS為地位指數(shù)，H為樹高，T為年齡。

運(yùn)用SPSS軟件，進(jìn)行公式(7)的擬合，擬合結(jié)果為K1=0.498、K2=11.401、K3=0.035、K4=0.071、K5=-0.002、K6=0.007、K7=1.867。由檢驗(yàn)數(shù)據(jù)所求得各項(xiàng)檢驗(yàn)指標(biāo)為：ERMS=0.649，EMA=0.415，R2=0.906。

除上述評(píng)價(jià)標(biāo)準(zhǔn)外，根據(jù)地位指數(shù)的定義可知，當(dāng)T為基準(zhǔn)年齡時(shí)，H等于IS。由此可得基準(zhǔn)年齡下的樹高預(yù)測(cè)值，并可計(jì)算與地位指數(shù)之間的均方根誤差。以2m為級(jí)距，地位指數(shù)由6～20m的樹高預(yù)測(cè)值分別為：6.662、8.304、10.151、12.145、14.205、16.244、18.176、19.945m。計(jì)算得到與地位指數(shù)之間的均方根誤差為0.298。

令I(lǐng)S=6、8、10、12、14、16、18、20，以年齡T為橫軸，以H為縱軸可繪制各地位指數(shù)下樹高生長(zhǎng)曲線，結(jié)果如圖1所示。

3.2 隨機(jī)森林回歸模型結(jié)果

分別以IS或H為輸出變量，由驗(yàn)證數(shù)據(jù)所得評(píng)價(jià)結(jié)果如表4所示。

表4 隨機(jī)森林回歸模型檢驗(yàn)結(jié)果

同理，由關(guān)系H=f(IS,T)可計(jì)算基準(zhǔn)年齡下各地位指數(shù)等級(jí)樹高預(yù)測(cè)值，并計(jì)算與地位指數(shù)之間的均方根誤差。以2m為級(jí)距，地位指數(shù)由6～20m的樹高預(yù)測(cè)值分別為：6.061、7.927、10.327、12.016、14.020、15.746、17.931、19.866m。計(jì)算得到與地位指數(shù)之間的均方根誤差為0.159。

固定IS值，由關(guān)系H=f(IS，T)將各年齡下樹高求出，以年齡T為橫軸，以H為縱軸繪制散點(diǎn)圖依次連接，也可得到樹高生長(zhǎng)曲線。但由于其無(wú)確定的表達(dá)式，應(yīng)用上不如Richards方程模型方便。

3.3 兩種模型結(jié)果對(duì)比分析

Richards方程模型最終擬合結(jié)果R2=0.906，ERMS=0.649，EMA=0.415；隨機(jī)森林回歸算法以H和IS為輸出變量所構(gòu)建模型R2均大于0.930，ERMS為0.510與0.530，EMA為0.352與0.281。除了傳統(tǒng)檢驗(yàn)方法，本文根據(jù)地位指數(shù)的定義，計(jì)算了基準(zhǔn)年齡下兩模型樹高預(yù)測(cè)值的均方根誤差來(lái)驗(yàn)證結(jié)果，最終得到預(yù)測(cè)誤差值分別為0.298與0.159，均在較小范圍內(nèi)，進(jìn)一步驗(yàn)證模型可行性。

與隨機(jī)森林回歸的這樣的非參數(shù)模型相比，Richards方程模型具有明確的表達(dá)式，當(dāng)需要預(yù)測(cè)樹木生長(zhǎng)趨勢(shì)，并繪制樹木生長(zhǎng)曲線時(shí)，僅需要將相應(yīng)地位指數(shù)代入公式，便可以得到生長(zhǎng)曲線表達(dá)式，從而快速繪制平滑的曲線；隨機(jī)森林回歸模型比Richards方程模型有更低的誤差值，ERMS降低了21%，EMA降低了15%，且具有更高的R2。該模型的另一個(gè)優(yōu)勢(shì)在于可以直接求得地位指數(shù)而無(wú)需通過(guò)迭代計(jì)算，求解效率更高。

4 結(jié)論

本文以廣西高峰林場(chǎng)杉木人工林為研究對(duì)象，利用二類調(diào)查小班數(shù)據(jù)與解析木數(shù)據(jù)，分別采用Richards方程與隨機(jī)森林回歸算法構(gòu)建多形地位指數(shù)模型。

結(jié)果表明，兩種模型均具有可行性，且各有優(yōu)勢(shì)。Richards方程模型更加適合樹木生長(zhǎng)趨勢(shì)預(yù)測(cè)，可以方便繪制樹木生長(zhǎng)曲線；隨機(jī)森林回歸模型優(yōu)勢(shì)在于具有更高的準(zhǔn)確率與更低的誤差值，可直接求得地位指數(shù)而無(wú)需迭代運(yùn)算。后續(xù)將利用構(gòu)建的模型，開發(fā)立地質(zhì)量評(píng)價(jià)系統(tǒng)。根據(jù)兩個(gè)模型的特點(diǎn)，結(jié)合不同需求選取不同的模型，以發(fā)揮每個(gè)模型的優(yōu)勢(shì)。

本文所構(gòu)建的模型僅適用于有林地，后續(xù)將利用坡向、坡度、坡位、土壤因子等立地因子替代模型中的地位指數(shù)，構(gòu)建適用于無(wú)林地的評(píng)價(jià)模型。