杜 志，陳振雄，孟京輝，王建軍，王景弟，李志剛

（1.國(guó)家林業(yè)和草原局 中南調(diào)查規(guī)劃設(shè)計(jì)院，湖南 長(zhǎng)沙 410014；2.北京林業(yè)大學(xué)國(guó)家林業(yè)和草原局森林經(jīng)營(yíng)工程技術(shù)研究中心，北京 100083；3.慈利縣林業(yè)局，湖南 慈利 427200；4.慈利縣天心閣林業(yè)工作站，湖南 慈利 427200）

森林生長(zhǎng)和收獲模型是描述林分或林木生長(zhǎng)規(guī)律的數(shù)學(xué)表達(dá)式，它是森林經(jīng)營(yíng)管理的重要工具[1]。森林生長(zhǎng)和收獲模型可以分為三類，即全林分模型、徑階分布模型以及單木生長(zhǎng)模型[2]。單木生長(zhǎng)模型是以單木生長(zhǎng)量為因變量，林木大小、林木間的競(jìng)爭(zhēng)以及立地條件為自變量的數(shù)學(xué)表達(dá)式[3-4]，它不僅能夠準(zhǔn)確地預(yù)估人工純林的生長(zhǎng)與收獲，還廣泛地應(yīng)用于異齡混交林和天然林中。此外，單木模型還可以模擬不同經(jīng)營(yíng)措施下林木的生長(zhǎng)動(dòng)態(tài)，它的模擬結(jié)果比其他幾類生長(zhǎng)模型更準(zhǔn)確，從而為森林質(zhì)量的精準(zhǔn)提升提供理論支持[5-8]。國(guó)內(nèi)外學(xué)者對(duì)多種樹(shù)種構(gòu)建了生長(zhǎng)模型，如孟憲宇等[9]以華北落葉松Larix principisrupprechtii人工林為研究對(duì)象，采用生長(zhǎng)量修正分析方法構(gòu)建了與距離無(wú)關(guān)的單木生長(zhǎng)模型；杜紀(jì)山[10]以單木的直徑為因變量，林木大小、林木間的競(jìng)爭(zhēng)以及立地條件為自變量，構(gòu)建了與距離無(wú)關(guān)的單木直徑生長(zhǎng)量模型。大部分學(xué)者構(gòu)建生長(zhǎng)模型采用的方法是傳統(tǒng)回歸方法，即最小二乘法。林分生長(zhǎng)和收獲模擬數(shù)據(jù)往往具有層次結(jié)構(gòu)、重復(fù)測(cè)量等特點(diǎn)[11-12]，表現(xiàn)為數(shù)據(jù)之間存在異方差和自相關(guān)，因而不滿足普通回歸分析中的獨(dú)立性假設(shè)，即獨(dú)立、等方差、正態(tài)分布，會(huì)得到有偏的參數(shù)估計(jì)[13-15]。而混合模型作為一種處理分組數(shù)據(jù)的有力工具，可以靈活地解決這一問(wèn)題，從而提高預(yù)測(cè)精度并解釋隨機(jī)誤差的來(lái)源[16-18]。

馬尾松Pinus massoniana是陽(yáng)性樹(shù)種，不耐庇蔭，喜光、喜溫，適合生于年均溫13～22 ℃、年降水量800～1 800 mm的地區(qū)。馬尾松分布范圍廣，是松屬分布最廣的樹(shù)種，東從東部平原地區(qū)，西至云貴高原和四川盆地，北自山東及河南南部，南至廣西、廣東、臺(tái)灣、海南等地，遍布全國(guó)16個(gè)省（區(qū)、市）。它是我國(guó)主要的用材樹(shù)種之一，在全國(guó)1 400多個(gè)喬木樹(shù)種的重要值排序中位居第三，在我國(guó)的森林資源總量中占有很重要的地位。根據(jù)第九次全國(guó)森林資源清查，馬尾松林面積和蓄積分別為804.30萬(wàn)hm2和6.26億m3，占我國(guó)喬木林面積和蓄積的4.47%和3.67%。馬尾松具有重要的生態(tài)價(jià)值和經(jīng)濟(jì)效益。因此，構(gòu)建精度較高的生長(zhǎng)預(yù)估模型對(duì)于馬尾松的經(jīng)營(yíng)管理十分必要。

本研究基于湘西地區(qū)2004、2009、2014年三期全國(guó)森林資源一類清查固定樣地?cái)?shù)據(jù)，運(yùn)用逐步回歸法，構(gòu)建馬尾松林胸高斷面積預(yù)估模型。同時(shí)在構(gòu)建模型的基礎(chǔ)上引入樣地隨機(jī)效應(yīng)，運(yùn)用單水平線性混合模型的方法構(gòu)建不同參數(shù)個(gè)數(shù)的混合效應(yīng)模型，篩選最佳的變量組合，同時(shí)利用獨(dú)立樣本數(shù)據(jù)對(duì)胸高斷面積預(yù)估模型進(jìn)行檢驗(yàn)，以期構(gòu)建的模型能為湘西地區(qū)馬尾松的生長(zhǎng)和經(jīng)營(yíng)提供參考依據(jù)。

1 數(shù)據(jù)和方法

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

試驗(yàn)數(shù)據(jù)來(lái)源于湘西地區(qū)2004、2009、2014年三期全國(guó)森林資源一類清查固定樣地，本次篩選的樣地?cái)?shù)量為180塊，樣地為面積667 m2的正方形。樣地信息包括樣地測(cè)量因子和樣木測(cè)量因子，樣木測(cè)量因子主要包括樹(shù)種名稱、胸徑、材積等。樣地測(cè)量因子包括坡度、坡向、坡位、海拔、土層厚度等立地因子以及優(yōu)勢(shì)樹(shù)種、郁閉度、林種、齡級(jí)等林分因子。剔除不滿足條件及數(shù)據(jù)異常的樣地后，3期數(shù)據(jù)共得到樣地100塊，馬尾松的株數(shù)為3 737株。

1.2 方 法

1.2.1 變量選擇及基礎(chǔ)模型的構(gòu)建

單木斷面積生長(zhǎng)模型常常描述為林木大小、競(jìng)爭(zhēng)及立地條件的線性函數(shù)，本研究充分考慮以上因素對(duì)于馬尾松單木斷面積的生長(zhǎng)，選擇的變量如下：

1）林木大?。褐饕ㄆ诔跣貜剑―BH）及其變形形式，即期初胸徑（DBH）、期初胸徑的平方（DBH2）、期初胸徑的對(duì)數(shù)（logDBH）以及期初胸徑的倒數(shù)（1/DBH）。

2）林分競(jìng)爭(zhēng)：為了使模型具有更好的通用性或適用性，選擇了以下與距離無(wú)關(guān)的競(jìng)爭(zhēng)指標(biāo)，即每公頃株數(shù)（NT）、每公頃斷面積（BA）、大于對(duì)象木的胸高斷面積和（BAL）和對(duì)象木胸徑與林分平均平方胸徑之比（RD）。

3）立地條件：充分考慮了海拔（EL）、坡度（SL）及坡向（ASP）的影響，同時(shí)依據(jù)Stage[19]提出的公式對(duì)坡度和坡向進(jìn)行了組合，即SlCos=SL×cos (ASP)。

不同形式的因變量會(huì)導(dǎo)致不同的模型表現(xiàn)。在預(yù)分析中，對(duì)不同形式的因變量進(jìn)行了測(cè)試，根據(jù)殘差圖和評(píng)價(jià)指標(biāo)，發(fā)現(xiàn)log[(DBH22-DBH2)+1]效果最好，因此選擇log[(DBH22-DBH2)+1]作為模型的因變量。其中，DBH、DBH2分別指期初胸徑和期末胸徑。因變量加入常數(shù)1的作用是避免斷面積生長(zhǎng)量為0時(shí)無(wú)法參與模型的計(jì)算。

基礎(chǔ)模型的構(gòu)建依據(jù)相關(guān)性分析，選擇與因變量相關(guān)性強(qiáng)的變量進(jìn)入模型，對(duì)進(jìn)入模型的變量進(jìn)行回歸系數(shù)顯著性檢驗(yàn)和共線性檢驗(yàn)，最終回歸系數(shù)顯著（P＜0.05）且方差膨脹因子（VIF）小于5的變量進(jìn)入模型。

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

采用單水平線性混合效應(yīng)模型的方法構(gòu)建湘西地區(qū)馬尾松單木斷面積預(yù)估模型，模型的一般形式如下：

式中：yi為第i個(gè)樣地的觀測(cè)值；Xi為已知固定效應(yīng)設(shè)計(jì)矩陣；β為未知的固定效應(yīng)參數(shù)向量；Zi為已知隨機(jī)效應(yīng)設(shè)計(jì)矩陣；bi為未知的隨機(jī)效應(yīng)參數(shù)向量；ei為誤差向量；n為樣地個(gè)數(shù)；D為隨機(jī)效應(yīng)協(xié)方差矩陣；σ2為模型的誤差方差值；Ri為組內(nèi)方差協(xié)方差結(jié)構(gòu)。在混合效應(yīng)模型構(gòu)建過(guò)程中需明確以下問(wèn)題：

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

混合模型構(gòu)建過(guò)程中最重要的一步是確定參數(shù)效應(yīng)，即哪個(gè)參數(shù)為固定參數(shù)，哪個(gè)參數(shù)為混合參數(shù)。本文中參數(shù)效應(yīng)的確定采用的是把基礎(chǔ)模型中所有參數(shù)均作為混合參數(shù)，對(duì)其不同組合進(jìn)行模擬，最后根據(jù)擬合統(tǒng)計(jì)量進(jìn)行選擇。

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

測(cè)試林業(yè)上常用的3種隨機(jī)效應(yīng)方差-協(xié)方差結(jié)構(gòu)即對(duì)角矩陣（DM）、復(fù)合對(duì)稱矩陣(CS)和廣義正定矩陣（General positive-definite matrix）來(lái)確定最優(yōu)的隨機(jī)效應(yīng)方差-協(xié)方差結(jié)構(gòu)。

3）誤差方差-協(xié)方差結(jié)構(gòu)的確定

由于建模數(shù)據(jù)來(lái)源于固定樣地的復(fù)測(cè)數(shù)據(jù)，這類數(shù)據(jù)往往存在異方差和自相關(guān)性的問(wèn)題，為了確定誤差的方差-協(xié)方差結(jié)構(gòu)，必須考慮以上兩方面問(wèn)題。在林業(yè)上通常采用式（2）來(lái)描述該問(wèn)題：

式中：σ2為誤差方差；Gi為誤差方差異質(zhì)性矩陣；Γi為誤差相關(guān)性結(jié)構(gòu)。

文中方差異質(zhì)性的描述選擇指數(shù)函數(shù)、冪函數(shù)和常數(shù)加冪函數(shù)3種，誤差相關(guān)性的表達(dá)選擇林業(yè)上常用的復(fù)合對(duì)稱（CS）、一階自回歸[AR(1)]與一階自回歸移動(dòng)平均結(jié)構(gòu)[ARMA(1,1)]。

1.2.3 模型的選擇與檢驗(yàn)

模型的篩選依據(jù)赤池信息準(zhǔn)則（AIC）、貝葉斯信息準(zhǔn)則（BIC）越小越好、對(duì)數(shù)似然值（Loglikelihood）越大越好的原則，同時(shí)要對(duì)篩選出的混合模型進(jìn)行似然比檢驗(yàn)（LRT），避免過(guò)多參數(shù)化問(wèn)題的產(chǎn)生。

利用篩選出的模型進(jìn)行獨(dú)立樣本數(shù)據(jù)檢驗(yàn)之前，需要利用式（3）計(jì)算隨機(jī)效應(yīng)參數(shù)值。

式中：D為隨機(jī)效應(yīng)參數(shù)方差-協(xié)方差矩陣；Zk為已知設(shè)計(jì)矩陣；Rk為誤差方差-協(xié)方差矩陣；ek為誤差向量。

模型檢驗(yàn)的結(jié)果通過(guò)平均絕對(duì)誤差（Bias）、均方根誤差（RMSE）和決定系數(shù)（R2）進(jìn)行評(píng)價(jià)。

2 結(jié)果與分析

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

本研究采用逐步回歸，依據(jù)相關(guān)性分析，選擇與因變量相關(guān)性強(qiáng)的變量進(jìn)入模型，對(duì)進(jìn)入模型的變量進(jìn)行回歸系數(shù)顯著性檢驗(yàn)和共線性檢驗(yàn)，最終模型形式如下：

參數(shù)估計(jì)及其統(tǒng)計(jì)描述如表1所示。1/DBH、BA、RD和EL對(duì)馬尾松單木斷面積的生長(zhǎng)有顯著影響（P＜0.05）且VIF小于5。該基礎(chǔ)模型的R2為0.364，從模型的殘差圖（圖1）可以看出存在異方差性。

2.2 混合效應(yīng)模型

考慮不同隨機(jī)參數(shù)的組合，對(duì)等式（4）進(jìn)行擬合，基于樣地效應(yīng)的馬尾松單木斷面積生長(zhǎng)混合效應(yīng)模型收斂的情況共有31種（表2）。根據(jù)AIC、BIC和Log-likelihood可以看出，混合效應(yīng)模型的模擬效果都優(yōu)于傳統(tǒng)模型。此外，只含有一個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，模擬3效果最好；當(dāng)含有2個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，模擬9效果最好；當(dāng)包含3個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，模擬20效果最好；當(dāng)含有4個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，模擬27效果最好；當(dāng)含有5個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，得到模擬32。與此同時(shí)，可以看出模擬27的擬合效果優(yōu)于其他任意一個(gè)組合的擬合效果，模擬精度最高。

表1 傳統(tǒng)方法模擬馬尾松單木胸高斷面積預(yù)估模型的結(jié)果Table 1 The simulated result of Pinus massoniana individual tree basal area prediction model based on conventional method

圖1 基礎(chǔ)模型的殘差圖和QQ圖Fig.1 Residual plot and QQ plot of the basic model

表2 胸高斷面積生長(zhǎng)量線性混合模型收斂結(jié)果?Table 2 The simulated results of basal area growth based on linear mixed model

為了比較變量參數(shù)個(gè)數(shù)的差異，本研究對(duì)以上5種模擬進(jìn)行了似然比檢驗(yàn)，其結(jié)果見(jiàn)表3。根據(jù)結(jié)果模擬27即等式（5）作為最佳參數(shù)組合的混合效應(yīng)模型。

式中：β1～β5是固定效應(yīng)；b1～b4是隨機(jī)效應(yīng)。

3種隨機(jī)效應(yīng)方差-協(xié)方差結(jié)果如表4所示。根據(jù)AIC、BIC、Log-likelihood和LRT，選擇廣義正定矩陣為最終的隨機(jī)效應(yīng)方差-協(xié)方差結(jié)構(gòu)。

由于基礎(chǔ)模型的殘差圖表明存在異方差性，因此使用指數(shù)函數(shù)、冪函數(shù)和常數(shù)加冪函數(shù)對(duì)模型誤差的方差結(jié)構(gòu)進(jìn)行模擬。從表5可以看出，引入異方差函數(shù)后，混合效應(yīng)模型在AIC、BIC和Loglikelihood方面得到了改善。根據(jù)LRT結(jié)果，在這3個(gè)異方差函數(shù)中，冪函數(shù)表現(xiàn)最佳，作為該模型最終的異方差函數(shù)。此外，由于加入自相關(guān)結(jié)構(gòu)后模型未能收斂，因此在最終的模型中未考慮。

表3 不同參數(shù)變量的模型似然比Table 3 Model likelihood ratios for different parameter variables

表4 胸高斷面積線性混合效應(yīng)模型隨機(jī)效應(yīng)參數(shù)的方差-協(xié)方差矩陣Table 4 The linear mixed-effects individual-tree basal area increment model using the random effects variance-covariance structures

表5 考慮異方差結(jié)構(gòu)矩陣后混合效應(yīng)模型的比較結(jié)果Table 5 The comparison results of mixed models with heteroscedasticity structure in different variance functions

確定了參數(shù)效應(yīng)、隨機(jī)效應(yīng)方差-協(xié)方差結(jié)構(gòu)以及誤差方差-協(xié)方差結(jié)構(gòu)，采用極大似然法在R語(yǔ)言中獲得的馬尾松單木斷面積生長(zhǎng)混合效應(yīng)模型，最終結(jié)果如下：

其中：

此外，混合效應(yīng)模型即式（6）的殘差圖和QQ圖如圖2所示。由圖2可以看出，引入隨機(jī)效應(yīng)后，模型的殘差圖得到了改善。

2.3 模型的檢驗(yàn)

本研究采用80%的數(shù)據(jù)共有2 857株馬尾松構(gòu)建了馬尾松的胸高斷面積預(yù)估模型，采用剩余的20%馬尾松數(shù)據(jù)對(duì)構(gòu)建的模型進(jìn)行檢驗(yàn)，比較線性混合模型與傳統(tǒng)最小二乘方法的精度。表6顯示，傳統(tǒng)最小二乘方法的3個(gè)指標(biāo)值計(jì)算結(jié)果Bias為0.328 1，RMSE為0.426 4，確定系數(shù)R2為0.333 2?？紤]樣地效應(yīng)的混合效應(yīng)模型的Bias為0.223 4，RMSE為0.315 2，確定系數(shù)R2為0.604 6?？梢钥闯龌旌闲?yīng)的平均絕對(duì)殘差與均方根誤差低于一般的線性模型，確定系數(shù)R2由0.333 2提高到0.604 6。

圖2 混合效應(yīng)模型的殘差圖與QQ圖Fig.2 Residual plot and QQ plot of final mixed-effects model

表6 傳統(tǒng)模型與混合效應(yīng)模型驗(yàn)證擬合統(tǒng)計(jì)Table 6 Fitting statistics for validation of the basic model and mixed-effects model

3 結(jié)論與討論

從林木本身大?。―BH）、林分的胸高斷面積（BA）、林分競(jìng)爭(zhēng)指數(shù)（RD）以及立地條件海拔（EL）對(duì)湘西地區(qū)馬尾松林生長(zhǎng)影響效果明顯。其中林分的生長(zhǎng)與胸徑的倒數(shù)呈負(fù)相關(guān)，與林木的胸徑呈正相關(guān)，即林木個(gè)體的林木胸徑越大，林分的競(jìng)爭(zhēng)能力也越強(qiáng)，林分林木生長(zhǎng)量也越大。王建軍等[2]、Zhao等[20]、Cannell等[21]分別以杉木Cunninghamia lanceolata、硬闊樹(shù)種以及云杉Picea asperata和松樹(shù)得到了相同的結(jié)論。BA是一個(gè)間接反映樹(shù)木間競(jìng)爭(zhēng)的指標(biāo)，胸高斷面積生長(zhǎng)量與林分的胸高斷面積呈反比，BA越大樹(shù)木間的競(jìng)爭(zhēng)壓力越大，胸高斷面積生長(zhǎng)量越小。RD是指林木的胸徑與平均胸徑的比值，它是反映林木競(jìng)爭(zhēng)強(qiáng)弱的指標(biāo)，其值越大說(shuō)明林木的競(jìng)爭(zhēng)力越強(qiáng)，樹(shù)木獲取光照、水分以及營(yíng)養(yǎng)物質(zhì)的能力越強(qiáng)，因此RD與胸高斷面積生長(zhǎng)量呈正比。

由于林木生長(zhǎng)和收獲數(shù)據(jù)不滿足獨(dú)立同分布的假設(shè)，本研究運(yùn)用層次混合效應(yīng)模型方法，構(gòu)建了湘西地區(qū)馬尾松林的胸高斷面積預(yù)估模型。根據(jù)AIC、BIC和Log-likelihood可以看出，考慮不同隨機(jī)參數(shù)組合混合效應(yīng)模型的模擬效果都優(yōu)于傳統(tǒng)模型，通過(guò)比較不同隨機(jī)參數(shù)組合的混合效應(yīng)模型，發(fā)現(xiàn)包含int、1/DBH、BA、RD的4個(gè)變量的組合效果最佳。同時(shí)選擇廣義正定矩陣為最終的隨機(jī)效應(yīng)方差-協(xié)方差結(jié)構(gòu)，考慮基礎(chǔ)模型方差存在異方差性，在消除異方差時(shí)，指數(shù)函數(shù)、冪函數(shù)和常數(shù)加冪函數(shù)均能在一定程度上縮小殘差的變化范圍，但冪函數(shù)形式的擬合效果較好。在模型檢驗(yàn)中，與傳統(tǒng)的模型相比，混合效應(yīng)模型的Bias、RMSE其相對(duì)值均顯著減少，模型的決定系數(shù)從0.333 2提高到了0.604 6?；旌闲?yīng)模型大大提高了預(yù)估的精度，混合效應(yīng)模型優(yōu)于傳統(tǒng)模型。

馬尾松為我國(guó)重要用材樹(shù)種，分布范圍廣，生長(zhǎng)立地條件復(fù)雜，以湘西地區(qū)為研究區(qū)域時(shí)，海拔對(duì)馬尾松單木斷面積的生長(zhǎng)呈現(xiàn)顯著影響。然而，在省或更大區(qū)域上，溫度、降水、土壤類型等因子可能是影響馬尾松生長(zhǎng)的重要參數(shù)。將來(lái)要充分利用好全國(guó)第九次森林資源連續(xù)清查成果，分析我國(guó)馬尾松分布地區(qū)，篩選生長(zhǎng)典型區(qū)域，綜合考慮地區(qū)氣候環(huán)境、立地條件等多因素，客觀全面地構(gòu)建單木斷面積預(yù)估模型，為馬尾松生長(zhǎng)預(yù)估和經(jīng)營(yíng)管理提供科學(xué)依據(jù)。