基于混合效應(yīng)的湖南馬尾松次生林單木生長模型

陳哲夫，肖化順，龍時勝

（1.湖南文理學(xué)院 資源環(huán)境與旅游學(xué)院，湖南 常德 415000；2.中南林業(yè)科技大學(xué) 林學(xué)院，湖南 長沙 410004）

馬尾松Pinus massoniana分布極廣，產(chǎn)于秦嶺、淮河流域以南，東起沿海低山丘陵，西至川西大相嶺東坡，南達(dá)華南南部，遍布于華中、華南，是我國南方集體林區(qū)經(jīng)營歷史悠久、最喜聞樂見的樹種之一，在我國國民經(jīng)濟(jì)發(fā)展、國土生態(tài)安全和應(yīng)對全球氣候變化等方面發(fā)揮著特別重要的作用。根據(jù)第七、八、九次全國森林資源連續(xù)清查結(jié)果，我國馬尾松林面積占比均在4.5%以上，蓄積占比均在3.5%以上，在所有優(yōu)勢樹種中面積與蓄積占比均處于前列[1]。

我國有關(guān)馬尾松生長模型的研究主要集中在人工林[2-3]，包括林分優(yōu)勢高[4]和自稀疏[5]模型等，關(guān)于馬尾松天然林的研究則相對較少，主要包括材積、胸徑[6]和樹高[7]等生長模型的研究。單木生長模型是林分生長模型的基礎(chǔ)，建立林木的單木生長模型，不僅對預(yù)測林木個體的生長過程具有重要意義，還可為林分整體的經(jīng)營管理提供經(jīng)營決策方案，從而提升森林的經(jīng)營水平[8-11]。傳統(tǒng)的建模方法雖然能夠取得一定的預(yù)測精度，但未能考慮立地之間的差異，導(dǎo)致其存在一定的局限性，而含有隨機(jī)參數(shù)的建模方法則可以很好地解決這一問題，其相比傳統(tǒng)應(yīng)用最小二乘法的建模方法取得了較高的預(yù)測精度[12-14]。近年來，混合效應(yīng)模型已廣泛應(yīng)用于林業(yè)中[15-16]，包括改進(jìn)Gompertz 模型對栓皮櫟樹高與胸徑的研究[17]，以及對華北落葉松單木冠幅[18]和樹高模型[19]的研究等，均取得了一定的成果，提升了模型的預(yù)測精度。本研究以湖南省馬尾松次生林為研究對象，建立其單木斷面積和材積的基礎(chǔ)生長模型，在此基礎(chǔ)上加入隨機(jī)參數(shù)，構(gòu)建基于樣地水平的湖南馬尾松次生林單木斷面積和材積生長模型，以期準(zhǔn)確預(yù)估該林分的生長過程，為其生長預(yù)測和可持續(xù)經(jīng)營提供技術(shù)指導(dǎo)。

1 研究區(qū)概況

湖南省位于我國中南部，長江中游，地處108°47′～114°15′E，24°38′～30°08′N，土地總面積2 118 萬hm2，其中林地面積為1 300 萬hm2，森林覆蓋率為59.57%，活立木蓄積5.05 億m3[1]。海拔24～2 122 m，地貌以800 m 以下的低山和丘陵為主，氣候?yàn)榈湫偷拇箨懶詠啛釒Ъ撅L(fēng)濕潤氣候，年均氣溫15～18oC，年日照數(shù)1 300～1 800 h，年降水量1 200～1 700 mm，土壤呈垂直性地帶性分布，主要為紅壤和黃壤。主要的喬木樹種有馬尾松、青岡櫟Cyclobalanopsis glauca、杉木Cunninghamia lanceolata、樟樹Cinnamomumbodinieri、楓香Liquidambar formosana和木荷Schima superba等，主要的灌木樹種包括厚皮香Ternstroemia gymnanthera、山茶Camellia japonica、鹿角杜鵑Rhododendron latoucheae和細(xì)枝柃Eurya loquaiana等，主要的草本植物為蘭花Cymbidium、麥冬Ophiopogon japonicus、芒萁Dicranopteris pedata和蕨Pteridium aquilinumvar.latiusculum等。

2 材料與方法

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

以湖南省第八次國家森林資源連續(xù)清查樣地數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，選取20 塊林分起源為天然林、樹種組成中馬尾松占比80%以上且林分郁閉在0.6以上的馬尾松次生林，樣地分布于湖南省9 個地（市、州）區(qū)，包括張家界、岳陽、常德、郴州、婁底、邵陽、衡陽、永州市和湘西土家族苗族自治州，能夠代表湖南省馬尾松次生林的整體生長情況。

馬尾松樣地為0.066 7 hm2的方形樣地，立地因子包括海拔、坡度、坡位、坡向、土壤類型、土層厚度、腐殖質(zhì)厚度和枯枝落葉厚度等，林分因子包括樹種、郁閉度、胸徑（D）、優(yōu)勢高（H）及其它衍生數(shù)據(jù)等，此外，還以龍時勝等[20]基于林木多期測定數(shù)據(jù)的異齡林年齡估計方法計算出樣地內(nèi)每株林木的年齡（A）。

對各樣地概況進(jìn)行統(tǒng)計（表1），其中海拔范圍為80～510 m，平均年齡為18～31 a，平均胸徑為8.4～15.3 cm，樣地內(nèi)的主要喬木樹種除馬尾松外，還包括櫟類、樟木、杉木、楓香和泡桐Paulownia等。

表1 樣地概況Table 1 The basic facts of sample plots

篩選出所有樣地中的馬尾松共993 株，以隨機(jī)抽樣的方法選取總樣本的2/3 作為建模數(shù)據(jù)，剩余1/3 的數(shù)據(jù)作為檢驗(yàn)樣本，分別對建模樣本和檢驗(yàn)樣本的各林分變量進(jìn)行特征統(tǒng)計（表2）。其中建模樣本的年齡變化為7～81 a，平均值為26 a；胸徑變化范圍為5.0～32.2 cm，平均值為12.6 cm；單木斷面積（Basal area，G）變化范圍為0.002 0～0.081 4 m2，平均值為0.014 5 m2；單株材積（Volume，V）變化范圍為0.004～0.548 m3，平均值為0.065 m3。檢驗(yàn)樣本的年齡變化為7～103a，平均值為26 a；胸徑變化范圍為5.1～30.4 cm，平均值為12.3 cm；單木斷面積變化范圍為0.002 0～0.072 6 m2，平均值為0.014 0 m2；單株材積變化范圍為0.005～0.472 m3，平均值為0.065 m3。

表2 建模數(shù)據(jù)與檢驗(yàn)數(shù)據(jù)特征統(tǒng)計Table 2 Characteristic statistics of modeling data and testing data

2.2 單木生長模型的選擇

選擇5 個具有生物學(xué)意義的理論方程，分別構(gòu)建湖南馬尾松次生林單木斷面積和材積生長的基礎(chǔ)模型，各模型及其表達(dá)式見表3，模型的擬合及參數(shù)計算在林業(yè)統(tǒng)計軟件Forstat 2.0 中的非線性回歸模塊中進(jìn)行。

2.3 混合效應(yīng)模型的構(gòu)建

在建立固定效應(yīng)模型的基礎(chǔ)上，考慮不同樣地間的差異對馬尾松林木斷面積、材積生長的影響，構(gòu)建基于樣地水平的混合效應(yīng)模型，其模型的一般表達(dá)形式如下：

表3 基礎(chǔ)模型與表達(dá)式?Table 3 The basic model and expression

式（1）中：Gij、Vij、Aij分別為第i塊樣地第j株林木的斷面積、材積和年齡；f(ijφ,Aij)為描述斷面積-年齡、材積-年齡關(guān)系的函數(shù)；ijφ為r×1 維的形參向量（r為形參個數(shù)）；β為p×1 維的固定效應(yīng)向量（p為固定參數(shù)個數(shù)）；εij為隨機(jī)誤差項(xiàng)；μi是服從期望為0、方差-協(xié)方差矩陣為Ψi(q×1)的獨(dú)立隨機(jī)效應(yīng)向量，且假定μi與εij間相互獨(dú)立（q為模型隨機(jī)參數(shù)個數(shù)）；Pij和Qij分別為對應(yīng)向量的設(shè)計矩陣；Rij為樣地內(nèi)誤差效應(yīng)的方差-協(xié)方差結(jié)構(gòu)矩陣。

構(gòu)建混合效應(yīng)模型，通常需要考慮以下三個方面的問題：

1）參數(shù)效應(yīng)的確定

混合效應(yīng)模型的參數(shù)包括固定效應(yīng)參數(shù)和隨機(jī)效應(yīng)參數(shù)，本研究將所有可能的隨機(jī)參數(shù)組合加入到模型中進(jìn)行混合效應(yīng)模擬，不同隨機(jī)參數(shù)的混合效應(yīng)模型擬合效果對比采用赤池信息準(zhǔn)則（AIC）和貝葉斯信息準(zhǔn)則（BIC），AIC 和BIC 值越小表示模擬效果越好，對于不同參數(shù)個數(shù)的模擬過程，還要進(jìn)行似然比檢驗(yàn)（LRT），當(dāng)P＜0.000 1 時，代表二者之間差異顯著，當(dāng)二者之間差異不顯著時，選擇參數(shù)更少的模擬以避免參數(shù)過多化的現(xiàn)象。

式（2）～（4）中：n為樣本數(shù)；p為參數(shù)個數(shù)；LL1和LL2分別為2 個對比模型的極大似然值。

2）確定樣地內(nèi)方差-協(xié)方差結(jié)構(gòu)（Rij）

樣地內(nèi)方差-協(xié)方差結(jié)構(gòu)能夠反映誤差異方差和自相關(guān)問題，本文中所獲取的研究數(shù)據(jù)為非連續(xù)性測量數(shù)據(jù)，因此在時間上不存在自相關(guān)，因而采用常用的誤差效應(yīng)方差-協(xié)方差結(jié)構(gòu)來描述：

式中：σ2為模型的誤差方差；In為描述樣地內(nèi)誤差的n×n維方差矩陣。

3）隨機(jī)效應(yīng)的方差-協(xié)方差結(jié)構(gòu)的確定（D）

隨機(jī)效應(yīng)的方差-協(xié)方差結(jié)構(gòu)矩陣可以反映不同樣地內(nèi)馬尾松斷面積和材積與年齡之間關(guān)系的差異性，其結(jié)構(gòu)會根據(jù)隨機(jī)參數(shù)個數(shù)的變化而變化，本研究以廣義正定矩陣來描述隨機(jī)效應(yīng)的方差-協(xié)方差結(jié)構(gòu)，以包含3 個隨機(jī)參數(shù)的方差-協(xié)方差結(jié)構(gòu)矩陣為例，其結(jié)構(gòu)如式（6）。

2.4 模型評價指標(biāo)的選取

基礎(chǔ)模型的選取指標(biāo)為確定系數(shù)（R2）、預(yù)估精度（P）和殘差平方和（SSE），其中R2和P越大，SSE 越小，模型的擬合效果越好；對于基礎(chǔ)模型與混合效應(yīng)模型的擬合效果對比，采用確定系數(shù)、平均誤差（Bias）和均方根誤差（RMSE）等指標(biāo)來評價，其中R2越大，Bias 和RMSE 越小，說明模型擬合效果越好。

式（7）～（11）中：yi為實(shí)際值；為預(yù)測值；為平均預(yù)估值；t0.05為置信水平a=0.05 時的t分布值；n為樣本數(shù)。

3 結(jié)果與分析

3.1 基礎(chǔ)生長模型的選擇

馬尾松單木斷面積基礎(chǔ)模型的擬合結(jié)果見表4。由表4可以看出，5 個模型除模型1 的擬合效果較差（R2=0.639，P=87.16%，SSE=0.036）外，其余4 個模型的擬合效果尚可，其確定系數(shù)均在0.73 以上，擬合精度也在96.5%以上，選擇其中擬合效果最佳的模型2（Logistic 模型）作為斷面積生長的基礎(chǔ)模型（R2=0.746，P=98.13%，SSE=0.025）；同理，將馬尾松單木材積生長模型確定為模型4（Richards 模型），其確定系數(shù)（R2=0.703）和預(yù)測精度（P=97.20）均為5 個模型中最大，其殘差平方和（SSE =1.034）為最?。ū?）。湖南馬尾松次生林單木斷面積和材積基礎(chǔ)模型如式（12）～（13）：

表4 斷面積生長方程擬合結(jié)果Table 4 Fitting results of basal area growth equation

表5 材積生長方程擬合結(jié)果Table 5 Fitting results of volume growth equation

3.2 斷面積和材積生長特征分析

根據(jù)所得到的斷面積生長方程，可繪制出湖南馬尾松次生林單木斷面積總生長量曲線，再根據(jù)總生長量計算出斷面積的平均生長量和連年生長量，并繪制二者之間的變化規(guī)律曲線（圖1）。由圖1可以看出，斷面積在0～20 a 生長較為緩慢，在21～50 a 生長速度較快，51 a 以后生長趨于平穩(wěn)，在90 a 時生長已基本停止，此時的斷面積達(dá)到最大值（0.046 4 m2）；斷面積連年生長量在35 a時達(dá)到最大值，為0.001 25 m2/a，隨后逐漸減小，在90 a 時其值趨近于0；平均生長量在50 a 時達(dá)到最大，為0.000 78 m2/a，此時的斷面積平均生長量與連年生長量相等，隨后斷面積平均生長量逐漸減小，其值也一直大于連年生長量。

圖1 斷面積生長特征曲線Fig.1 Growth characteristic curve of basal area

同理，對材積生長特性進(jìn)行分析，發(fā)現(xiàn)材積在0～60 a 均有較快的生長速度，61 a 以后生長速度減緩，生長趨于平穩(wěn)，在150 a 時生長已基本停止，此時的材積達(dá)到最大值0.356 5 m3；材積連年生長量在40 a 時達(dá)到最大值（0.00 669 m3/a），隨后逐漸減小，在150 a 時其值趨近于0；平均生長量在65 a 時達(dá)到最大，為0.004 25 m3/a，此時的材積平均生長量與連年生長量相等，隨后材積平均生長量逐漸減小。

圖2 材積生長特征曲線Fig.2 Growth characteristic curve of volume

3.3 基于混合效應(yīng)的單木生長模型

在確定基礎(chǔ)模型的前提下，加入含樣地效應(yīng)的隨機(jī)參數(shù)，構(gòu)建基于混合效應(yīng)的湖南馬尾松次生林單木斷面積和材積生長模型。首先，在式（12）中加入所有可能的隨機(jī)參數(shù)組合，發(fā)現(xiàn)共有6 種收斂的模擬過程，擬合結(jié)果（表6）顯示，所有含隨機(jī)參數(shù)的模擬其擬合效果均優(yōu)于基礎(chǔ)模型（AIC=-5 417.850，BIC=-5 401.364），含1 個隨機(jī)參數(shù)的模擬共有3 種收斂的情況，其中模擬4（AIC=-5 680.460，BIC=-5 658.007）的結(jié)果優(yōu)于模擬2 和模擬3，含2 個隨機(jī)參數(shù)的擬合有2 種收斂的結(jié)果，模擬5（AIC=-6 348.868，BIC=-6 317.433）的效果優(yōu)于模擬6，含3個隨機(jī)參數(shù)的模擬僅有一種結(jié)果模擬7（AIC=-6 458.433，BIC=-6 413.526）。對于含不同隨機(jī)參數(shù)個數(shù)的模擬過程進(jìn)行似然比檢驗(yàn)，發(fā)現(xiàn)模擬5 的結(jié)果顯著優(yōu)于模擬4（LRT=672.408，P＜0.0 001），而模擬7 的結(jié)果顯著優(yōu)于模擬5（LRT=115.566，P＜0.000 1），最終將含隨機(jī)參數(shù)μ1、μ2、μ3的混合效應(yīng)模型作為最優(yōu)模型。同理，材積的隨機(jī)參數(shù)模擬結(jié)果顯示模擬6（AIC=-4 211.681，BIC=-4 166.774）的擬合效果最佳（表7），將其隨機(jī)參數(shù)確定為μ1、μ2、μ3，因此，基于混合效應(yīng)的湖南馬尾松次生林單木斷面積和材積生長模型一般表達(dá)式為：

式（14）～（15）中：μ1、μ2、μ3為樣地水平的隨機(jī)效應(yīng)參數(shù)。

表6 斷面積混合效應(yīng)模型擬合結(jié)果比較Table 6 Comparison of fitting results of basal area mixed effect model

表7 材積混合效應(yīng)模型擬合結(jié)果比較Table 7 Comparison of fitting results of basal area mixed effect model

3.4 模型評價分析

通常對比基礎(chǔ)模型與混合效應(yīng)模型的殘差分布可以直觀地判斷模型擬合效果的差異及誤差的異方差性（圖3）。為確定所擬合的模型是否存在異方差，分別繪制了馬尾松斷面積、材積生長基礎(chǔ)模型與混合效應(yīng)模型的殘差分布（圖3～4）。結(jié)果（圖3～4）顯示，基礎(chǔ)模型的殘差分布范圍大（圖3A，圖4A），分布不均勻，隨著預(yù)測值的逐漸增大，其殘差值也隨之增大，存在一定的異方差；而混合效應(yīng)模型的殘差分布范圍明顯減?。▓D3B，圖4B），且分布較為均勻，沒有出現(xiàn)明顯的不規(guī)則形狀（如啞鈴型、喇叭形和拋物線形等），表明混合效應(yīng)模型的異方差已基本消除。

圖3 斷面積殘差分布Fig.3 Residual plots of basal area

圖4 材積殘差分布Fig.4 Residual plots of volume

對馬尾松斷面積和材積生長的基礎(chǔ)模型與混合效應(yīng)模型進(jìn)行參數(shù)估計與模型評價，各模型的參數(shù)估計值與方差組成見表8。各模型的擬合統(tǒng)計指標(biāo)結(jié)果顯示，混合效應(yīng)模型的AIC、BIC、平均誤差（Bias）和均方根誤差（RMSE）值均小于基礎(chǔ)模型，在數(shù)量上，模型14 的Bias 值由模型12 的0.000 26 降低到0.000 01，其RMSE 值由0.006 40降低到0.001 95，模型15 的該兩項(xiàng)指標(biāo)也較模型13 有顯著降低，分別從0.001 73 降低到0.000 13、由0.038 20 降低到0.000 20；混合效應(yīng)模型的確定系數(shù)（R2）相比基礎(chǔ)模型有很大的提升，其中模型14 的R2較模型12 提升了36.2%，模型15 的R2較模型13 提升了37.0%，相比基礎(chǔ)模型，混合效應(yīng)模型的預(yù)測精度（P）也有一定的提升，分別提升了1.84%和2.67%。綜上所述，含以樣地為隨機(jī)效應(yīng)的混合效應(yīng)模型擬合效果優(yōu)于基礎(chǔ)生長模型，能夠更好地預(yù)測馬尾松的生長過程。

表8 基礎(chǔ)模型與混合效應(yīng)模型擬合效果對比Table 8 Fitting result of basic model and mixed effect model

4 結(jié)論與討論

本研究分別以馬尾松單木斷面積和材積為因變量，構(gòu)建其關(guān)于年齡的生長模型，在5 個基礎(chǔ)模型中Logistic 方程能夠最好地反映斷面積的生長規(guī)律，該模型有最大的預(yù)測精度（P=98.13%）和確定系數(shù)（R2=0.746），同時其殘差平方和（SSE=0.025）最??；而材積生長的最優(yōu)基礎(chǔ)模型為Richards 方程，其R2和P最大（0.703，97.20%），SSE 最?。?.034）。在此基礎(chǔ)上，構(gòu)建了基于樣地水平的混合效應(yīng)模型，以期消除立地間的差異，對不同隨機(jī)參數(shù)的混合效應(yīng)模型進(jìn)行模擬，對比不同隨機(jī)效應(yīng)模型的赤池信息準(zhǔn)則（AIC）、貝葉斯信息準(zhǔn)則（BIC），對不同隨機(jī)參數(shù)個數(shù)的模擬進(jìn)行似然比檢驗(yàn)（LRT），以避免參數(shù)過多化的現(xiàn)象，發(fā)現(xiàn)斷面積和材積生長的混合效應(yīng)模型其隨機(jī)參數(shù)均為μ1、μ2、μ3，其擬合效果均顯著優(yōu)于其他模型過程（P＜0.000 1）。相比基礎(chǔ)模型，混合效應(yīng)模型的擬合效果明顯更優(yōu)，其平均誤差（0.000 01，0.000 13）和均方根誤差（0.001 95，0.000 20）顯著降低，確定系數(shù)（0.974，0.984）大幅提升，預(yù)測精度（99.936%，99.798%）也有所提高?；旌闲?yīng)模型有更好的擬合效果，能夠更精確地預(yù)測馬尾松次生林的生長規(guī)律，為其提供更有效的經(jīng)營措施。

本研究的研究數(shù)據(jù)來自湖南省9 個地區(qū)，具有很強(qiáng)的代表性，其所構(gòu)建的生長模型能夠很好地反映湖南馬尾松次生林的整體生長規(guī)律，所形成的生長預(yù)測方程也能夠廣泛應(yīng)用于湖南省各地區(qū)。研究結(jié)果與前人構(gòu)建的啞變量模型[21]和混合效應(yīng)模型[22-23]所得到的結(jié)果相似，均提示了模型的預(yù)測精度，具有更強(qiáng)的適用性。通常，單木生長與年齡、立地質(zhì)量和林分密度（或競爭因子）密切相關(guān)，本研究選取林分密度相似的馬尾松次生林樣地，以年齡為自變量、樣地為隨機(jī)效應(yīng)構(gòu)建單木生長模型，但未考慮到林木之間的競爭關(guān)系，為研究中存在的不足之處，今后將在建模過程中加入林木競爭因子，以期更加精確地反映林木的生長規(guī)律。