吳旭平,呂 勇,張雄清,易 烜,朱光玉

(1.中南林業(yè)科技大學(xué),長沙 410004；2.中國林業(yè)科學(xué)研究院林業(yè)研究所,北京 100091；3.青羊湖國有林場,長沙 410600)

0 引言

林分斷面積是指林分中所有林木胸高斷面積之和,作為常用的林分密度指標之一,具有較高的穩(wěn)定性和預(yù)估性[1]；因而,研究林分斷面積生長對森林的收獲預(yù)估具有重要意義[2]。

為了構(gòu)建斷面積生長模型,過去諸多學(xué)者對此進行了大量的研究,基礎(chǔ)模型常常選用Richards和Schumacher方程,模型自變量主要包括年齡、密度、立地指數(shù)等因子[3]。其中李永慈等[4]研究得到:林分斷面積生長模型是全林生長模型的核心,其精度直接影響到系統(tǒng)整體的預(yù)測精度。符利勇等[5]通過構(gòu)建混合效應(yīng)模型，有效解決了傳統(tǒng)的建模方法無法反映不同立地水平及林分水平對蒙古櫟(QuercusmongolicaFisch.ex Ledeb)生長的隨機影響,從而導(dǎo)致精度低的問題。但目前將該方法應(yīng)用于建立基于立地類型的杉木人工林林分斷面積生長模型的研究較少。

湖南杉木(Cunninghamialanceolata)人工林面積和蓄積分別占湖南全省人工林面積和蓄積總量的33%和41%[6]。湖南是杉木人工林主要分布的區(qū)域之一,以該區(qū)域杉木人工林為研究對象,參考前人的研究,利用方差分析、回歸分析法(最小二乘法)、非線性混合效應(yīng)模型、K-means聚類等方法構(gòu)建了基于立地效應(yīng)的林分斷面積生長模型,客觀解釋了立地因子對杉木人工林斷面積生長的影響規(guī)律,提高了林分斷面積生長模型的模擬精度,為湖南杉木人工林的科學(xué)經(jīng)營管理提供了理論支撐。

1 研究區(qū)概況

湖南省位于我國中南部,地理坐標24°38′～30°08′N,108°47′～114°15′E。土地總面積21.18萬km2,林地面積約1 300萬hm2,活立木蓄積5.05億m3,森林覆蓋率59.57%[7]。海拔為24～2 122m,大部分地區(qū)海拔高度100～800m；年日照時數(shù) 1 300～1 800h；年平均降水量1 200～1 700mm,雨量充沛,水熱充足。土壤類型以紅壤和黃壤為主。研究區(qū)主要的喬木樹種為杉木、馬尾松(PinusmassonianaLamb)、濕地松(Pinuselliottii)、柏木(Cupressusfunebris)、樟樹(Cinnamomumbodinieri)、木荷(Schimasuperba)、櫸樹(Zelkovaserrata)等。

2 材料與方法

2.1 數(shù)據(jù)來源

以湖南省150塊杉木人工林樣地為研究對象,在湖南省各市縣設(shè)置杉木林樣地進行調(diào)查。調(diào)查因子主要包括地理坐標、海拔、坡度、坡向、坡位、土壤類型、胸徑、樹高、林分年齡等。各因子統(tǒng)計如表1所示。

表1 杉木林分調(diào)查因子統(tǒng)計Tab.1 Statistics of stand factors for Cunninghamia lanceolata

2.2 立地類型劃分

1)立地因子等級劃分。

選用調(diào)查得到的海拔(HB)、坡向(PX)、坡位(PW)、坡度(PD)、土壤厚度(TH)、土壤類型(TL)等6個立地因子,并參照《測樹學(xué)》[8]中的各類因子分級標準進行分級。其中,為了分析不同海拔對杉木生長的影響,將海拔按照100m進行分級[1]。具體劃分標準如表2所示。

表2 立地因子等級劃分Tab.2 Site factor classification

2)顯著性分析及立地類型劃分。

采用方差分析,以林分胸高斷面積為因變量,以年齡和立地因子為自變量,進行顯著性因子篩選[9]；然后以顯著性因子為對象,通過將分級之后的立地因子進行組合來劃分立地類型(Site type),簡稱ST。

2.3 基礎(chǔ)模型篩選

基礎(chǔ)模型的選擇影響林分斷面積生長模型擬合的結(jié)果[10-11]。本文參考朱光玉等[12]構(gòu)建林分斷面積時所選用的基礎(chǔ)模型,來作為本次湖南杉木人工林林分斷面積生長預(yù)測的基礎(chǔ)模型。候選基礎(chǔ)模型如表3所示。

表3 候選基礎(chǔ)模型Tab.3 Basic models candidates

2.4 非線性混合效應(yīng)模型

非線性混合效應(yīng)模型是回歸函數(shù)依賴于固定效應(yīng)參數(shù)和隨機效應(yīng)參數(shù)的非線性關(guān)系而建立的模型[13],其一般形式為:

Yi=f(β,μi,Xi)+εi

(1)

式中:Yi與Xi分別代表第i個樣地的因變量向量和自變量向量；εi為誤差變量；β和μi分別為固定效應(yīng)參數(shù)向量和隨機效應(yīng)參數(shù)向量[14]。

2.5 模型精度評價

使用Forstat 2.2和R語言等軟件進行數(shù)據(jù)統(tǒng)計,得到模型參數(shù)。用確定系數(shù)(R2)、平均絕對誤差(MAE)、預(yù)估精度(P)、均方根誤差(RMSE)、貝葉斯信息準則(BIC)及赤池信息準則(AIC)對模型的各個預(yù)測效果進行精度評價,其中MAE,RMSE,AIC,BIC值越小,且R2和P值越大,表明預(yù)測結(jié)果越好[15]。

確定系數(shù):

(2)

平均絕對誤差:

(3)

均方根誤差:

(4)

預(yù)估精度:

(5)

赤池信息準則:

AIC=-2×lnl+2×p

(6)

貝葉斯信息準則:

BIC=-2×lnl+lnn×p

(7)

3 結(jié)果與分析

3.1 顯著性分析及立地類型劃分結(jié)果

采用方差分析,以林分胸高斷面積為因變量,以立地因子和年齡為自變量,對立地因子進行顯著性因子篩選[1]。當“Pr>F”的值大于0.05,即認為該因子對林分斷面積影響不顯著,否則影響顯著,結(jié)果如表4所示。

表4 立地因子顯著性檢驗Tab.4 Site factor significance test

由表4可得:HB,PD,PW,TH,TL對林分斷面積的生長具有顯著影響,其顯著性順序為TL>TH>PD>PW>HB。

以方差分析篩選所得到的HB,PD,PW,TH,TL為對象,通過對顯著性因子的等級值進行組合劃分,可將150個樣地劃分為107個初始立地類型(ST)。

3.2 基礎(chǔ)模型模擬結(jié)果

1)林分密度指數(shù)計算。

林分密度指數(shù)為常用密度指標之一,常用于構(gòu)建林分斷面積[16],其表達式為:

SDI=N×(D1/D)β

(8)

式中:β為自然稀疏率,N為林分每公頃株數(shù),D為林分平均胸徑,D1為標準平均胸徑(我國一般取D1=10cm),SDI為密度指數(shù)。

使用二次剔除不足立木度的樣地方法來估算β。首先,建立回歸方程lnN=a1-b1lnD,剔除lnN

具體公式為:

lnN=α-βlnD

(9)

在Forstat 2.2中對上述模型進行非線性擬合,最終得模型表達式如(10)式所示,其調(diào)整確定系數(shù)Ra2=0.71。

lnN=4.5273-0.96053lnD

(10)

將所得自然稀疏率β=-0.96053結(jié)合相關(guān)變量代入SDI表達式即可計算各樣地林分密度指數(shù)。

2)基礎(chǔ)模型的擬合及最優(yōu)模型確定。

以湖南省150塊杉木人工林斷面積-年齡數(shù)據(jù)為基礎(chǔ),利用Forstat 2.2對候選基礎(chǔ)模型M1—M8進行參數(shù)擬合與篩選,得到不同基礎(chǔ)模型參數(shù)估計值,擬合結(jié)果如表5所示。

由表5可知,模型M5:確定系數(shù)(R2=0.763 6)及P值(97%)最大,平均絕對誤差(MAE=3.08)及均方根誤差(RMSE=4.40)最小。可見,8個候選基礎(chǔ)模型擬合結(jié)果中,模型M5擬合效果最好,因此,確定M5為最優(yōu)斷面積生長模型。公式為:

表5 候選基礎(chǔ)模型擬合結(jié)果Tab.5 The fitting results of candidate basic models

G=HT(a+b/T)×(SDI/1000)(c+d/T)×e(q+w/T)

(11)

式中:G表示林分斷面積；T表示林分年齡；SDI表示林分密度指數(shù)；HT表示林分優(yōu)勢高；a,b,c,d,q,w表示模型參數(shù)。

3.3 基于立地類型的非線性混合效應(yīng)模型構(gòu)建

由于隨機效應(yīng)構(gòu)造受隨機因子個數(shù)及其水平數(shù)的影響,可以按任意方式相互組合衍生出多種類型的非線性混合效應(yīng)模型。本文以主導(dǎo)因子劃分的立地類型(ST)作為隨機效應(yīng),構(gòu)建非線性混合效應(yīng)模型進行分析[18-19]。結(jié)合基礎(chǔ)模型參數(shù)個數(shù)和立地類型,運用R語言的非線性混合效應(yīng)模型板塊,將立地類型分別加在各個參數(shù)上及同時加在多個參數(shù)上進行隨機效應(yīng)擬合,并利用評價指標AIC和BIC等確定最優(yōu)模型結(jié)構(gòu),擬合結(jié)果如表6所示。

表6 基于立地類型的非線性混合效應(yīng)模型擬合Tab.6 Nonlinear mixed-effect model fitting based on site type

由表6不同隨機效應(yīng)因子組合類型的擬合結(jié)果顯示,構(gòu)建單個參數(shù)及多個參數(shù)的混合效應(yīng)模型的AIC,BIC相對基礎(chǔ)模型M5有所降低,R2大幅提高。其中,立地類型(ST)作為隨機效應(yīng)加在參數(shù)c上面的混合效應(yīng)模型R2(0.895 1)最高,AIC(872.41)、BIC(896.50)最低。因為基礎(chǔ)模型是總體平均模型,而混合模型考慮了各種立地因子的影響。結(jié)果表明,基于立地類型的混合效應(yīng)模型明顯優(yōu)于基礎(chǔ)模型。

3.4 基于立地類型組的非線性混合效應(yīng)模型構(gòu)建

依據(jù)HB,PD,PW,TH,TL這5個主導(dǎo)因子,可將湖南杉木人工林劃分為107個立地類型(ST),相應(yīng)的混合模型隨機效應(yīng)參數(shù)有107個水平值。由于立地類型數(shù)過多,不利于混合模型的有效應(yīng)用[1]。為了簡化混合效應(yīng)模型和進一步提高模型模擬精度,本文將初始立地類型(ST)應(yīng)用到模型(11)式擬合的隨機效應(yīng)參數(shù)值進行K-means聚類得到立地類型組(Site Type Group),簡稱STG。進而構(gòu)建基于立地類型組的非線性混合效應(yīng)模型。

1)立地類型組的劃分。

本文以聚類精度≥90%為標準,將107個立地類型的隨機效應(yīng)參數(shù)值進行聚類,隨機效應(yīng)參數(shù)值接近的立地類型合并成為立地類型組(STG),結(jié)果如表7所示。

表7 立地類型聚類結(jié)果Tab.7 Site type clustering results

由表7可得:聚類數(shù)為5的時候,聚類精度達到了91.0%,達到聚類精度要求。

2)基于立地類型組的非線性混合效應(yīng)模型擬合結(jié)果。

運用R語言的非線性混合效應(yīng)模型板塊,將立地類型組分別加在各個參數(shù)上及同時加在多個參數(shù)上進行隨機效應(yīng)擬合,擬合結(jié)果如表8所示。

表8的擬合結(jié)果顯示,基于立地類型組的非線性混合效應(yīng)模型的AIC,BIC相對模型(11)式及基于立地類型的非線性混合效應(yīng)模型均有所降低,而R2大幅提高。其中,立地類型組(STG)作為隨機效應(yīng)加在參數(shù)c上面的混合效應(yīng)模型R2(0.920 2)最高,AIC(749.40)、BIC(770.48)最低。這說明,基于立地類型組的混合效應(yīng)模型明顯優(yōu)于基礎(chǔ)模型及基于立地類型的混合效應(yīng)模型。最終模型為:

表8 基于立地類型組的非線性混合效應(yīng)模型擬合Tab.8 Fitting of nonlinear mixed effects model based on site type group

G=HT(a+b/T)×(SDI/1000)(c+ci+d/T)×e(q+w/T)+ε

(12)

式中:G為林分斷面積；HT為林分優(yōu)勢高；T為林分年齡；SDI為林分密度指數(shù)；ci為立地類型組的隨機效應(yīng)參數(shù)；ε為誤差項。

3.5 模型精度評價

將聚類后的立地類型組(STG)作為隨機效應(yīng)進行非線性混合效應(yīng)模擬,利用MAE,RMSE,P,R2等4個評價指標進行模型評價,并與基礎(chǔ)模型(None)、初始立地類型(ST)模擬結(jié)果進行對比分析,其結(jié)果如表9所示。

從表9中可以看出,在將聚類后的立地類型組作為隨機效應(yīng)加入模型后,林分斷面積生長模型的R2值從0.763 7提高到0.920 2、提高了20.49%,MAE降低了17.21%,RMSE降低了14.77%。

表9 模型精度評價Tab.9 Model evaluation of three models

為了更加直觀地對比基礎(chǔ)模型與混合效應(yīng)模型的模擬效果,分別建立最優(yōu)基礎(chǔ)模型(11)式和最優(yōu)非線性混合效應(yīng)模型(12)式的殘差關(guān)系圖,其結(jié)果如圖1、圖2所示。

圖1 基礎(chǔ)模型(11)式殘差分布情況Fig.1 Residual distribution of the basic model(11)

圖2 混合效應(yīng)模型(12)式殘差分布情況Fig.2 Residential distribution mixed effects model(12) type

從圖1和圖2可以明顯看出:基礎(chǔ)模型與混合效應(yīng)模型的殘差分布都較隨機地分布于坐標軸兩側(cè),無明顯異方差現(xiàn)象；相比較基礎(chǔ)模型(圖1),混合效應(yīng)模型的殘差分布范圍更小,且更加集中于坐標軸兩側(cè)。可見:基于立地隨機效應(yīng)的林分斷面積模型可以極大提高模型精度,同時利用K-means聚類劃分立地類型組的方法,可以進一步提高模型模擬精度且解決了復(fù)雜立地類型的模型使用問題。

4 結(jié)論與討論

采用方差分析、回歸分析法(最小二乘法)、非線性混合效應(yīng)模型、K-means聚類等方法構(gòu)建含立地效應(yīng)的湖南杉木人工林林分斷面積生長模型,客觀解釋了立地因子對杉木人工林斷面積生長的影響,提高了林分斷面積生長模型的模擬精度。

1)采用方差分析,以林分胸高斷面積為因變量,以立地因子和年齡為自變量,對立地因子進行顯著性篩選,其結(jié)果表明HB,PD,PW,TH,TL對林分斷面積的生長具有顯著影響,其顯著性順序為TL>TH>PD>PW>HB。

2)通過使用回歸分析法對8個候選基礎(chǔ)模型進行模擬,其中Schumacher(M5)的確定系數(shù)最高(R2=0.763 6),因此,選擇此方程作為林分斷面積生長模擬的基礎(chǔ)模型。

3)本文在確定基礎(chǔ)模型之后,構(gòu)建非線性混合效應(yīng)模型,期望能得到湖南省杉木人工林林分斷面積生長規(guī)律。結(jié)果表明,構(gòu)建混合效應(yīng)模型可以顯著提高林分斷面積建模精度,其確定系數(shù)(R2)從0.763 6提高0.895 1。胡松[1]用相同方法分析了不同林分類型與立地類型差異對櫟類林分斷面積生長的影響；李春明等[3]在對比傳統(tǒng)的回歸模型方法與混合模型方法構(gòu)建落葉松云冷杉林分斷面積模型之后,得到混合模型方法精度更高的結(jié)果。這充分說明運用混合模型方法構(gòu)建林分斷面積模型是合理且有效的。為方便建模,構(gòu)建的斷面積模型中僅考慮地形、地貌等立地因子的影響,未能考慮土壤、氣候等因子的影響。在后續(xù)研究中,可考慮使用地形、地貌數(shù)據(jù)先構(gòu)建杉木人工林立地指數(shù)模型求得各樣地的地位指數(shù),進而再考慮土壤因子構(gòu)建斷面積模型,或許模型精度會進一步提升。

4)通過K-means對立地類型進行聚類得到立地類型組,然后進行混合效應(yīng)模擬。劉飛虎[20]采用此方法分析了不同林層對櫟類次生林斷面積生長的影響。由于聚類是指將不同的樣本總體劃分成不同的類型,及各類型的樣本個體之間的差異盡可能小。所以本文構(gòu)建的基于立地類型組的非線性混合效應(yīng)模型可以進一步提高建模精度,其確定系數(shù)(R2)從0.895 1提高到0.920 2。