湖南楠木次生林生長與收獲預(yù)估模型

龔召松，曾思齊，賀東北，石振威，龍時勝，姜興艷，陳亞文

（1.四川省林業(yè)勘察設(shè)計研究院，四川 成都 610081；2.中南林業(yè)科技大學(xué) 林學(xué)院，湖南 長沙 410004；3.國家林業(yè)和草原局 中南調(diào)查規(guī)劃設(shè)計院，湖南 長沙 410014；4.貴州省安順市西秀區(qū)林業(yè)局，貴州 安順 561000；5.四川省林業(yè)和草原調(diào)查規(guī)劃院，四川 成都 610081）

楠木Phoebe zhennan，為亞熱帶常綠闊葉喬木，是一種極其珍貴的樹種，數(shù)量稀少，生長緩慢，現(xiàn)為國家Ⅱ級珍稀漸危種，已被列入中國國家重點保護野生植物，主要分布在我國南方地區(qū)的湖北、貴州、湖南及四川等地。林分的生長和收獲模型能夠?qū)崿F(xiàn)對林分生長量和收獲量的預(yù)估[1]。林分生長量和收獲量是從兩個角度定量說明森林的變化狀況，林分收獲量是林分生長量積累的結(jié)果，而生長量又是森林的生產(chǎn)速度，它體現(xiàn)了特定期間的收獲量概念。1721年，德國科學(xué)家Reanmur 最早對林分生長與收獲模型進行了研究，近年來Adame 等[2]、Budhathoki 等[3]和Alegria 等[4]也對其進行了深入的研究。國內(nèi)對林分生長與收獲模型的研究最早始于20世紀(jì)80年代，邵國凡[5]對黑龍江紅松生長模型進行了研究，隨后，胥輝[6]對思茅松、華偉平等[7]對黃山松、盧軍等[8]和馬武等[9]對蒙古櫟的生長與收獲模型進行了研究。大量研究表明，傳統(tǒng)的建模方法存在許多的局限性，不能充分考慮樹種和立地之間的差異，而含有隨機參數(shù)的混合效應(yīng)模型比傳統(tǒng)應(yīng)用最小二乘法的建模方法取得了更高的預(yù)測精度[10-13]。

本研究以湖南省1989—2014年6 期一類清查數(shù)據(jù)中的楠木次生林樣地為研究對象，建立其相容性生長與收獲預(yù)估模型，在此基礎(chǔ)上，引入以樣地為隨機效應(yīng)的混合效應(yīng)聯(lián)立方程系統(tǒng)，對自然狀態(tài)下林分的生長動態(tài)進行預(yù)測和模擬，以期精準(zhǔn)預(yù)測其斷面積和蓄積量的生長，為湖南省楠木次生林的提質(zhì)增量提供參考依據(jù)以及制定科學(xué)合理的可持續(xù)經(jīng)營方案。

1 研究區(qū)概況

湖南省地處中國中部、長江中游，地理坐標(biāo)為108°47′～114°15′E，24°38′～30°08′N，海 拔在24～2 122.35 m 之間，地貌主要以海拔800 m以下的低山和丘陵為主，氣候為大陸性亞熱帶季風(fēng)濕潤氣候，年日照時數(shù)為1 300～1 800 h，年均氣溫為15～18oC，雨量充沛，年降水量為1 200～1 700 mm。土壤主要為紅壤（51.0%）和黃壤（12.62%），森林資源豐富，全省林業(yè)用地面積為1 299.8 萬hm2，森林覆蓋率達59.57%。主要樹種有：柏木Cupressus funebris、紅花木蓮Manglietia insignis、樂昌含笑Michelia chapensis、馬尾松Pinus massonana、楠木、青岡櫟Cyclobalanopsis glauca、濕地松Pinus elliottii、杉木Cunninghamialanceolata和樟樹Cinnamomum camphora等。

2 研究方法

2.1 數(shù)據(jù)來源與處理

本研究以湖南省1989—2014年6 期的一類清查數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，從中篩選出117 塊有楠木分布的樣地，并篩選出以楠木為主要樹種的樣地55塊，以前一期的數(shù)據(jù)作為林分的期初生長狀況，后一期的數(shù)據(jù)作為林分的期末生長狀況，間隔期為5 a，6 期數(shù)據(jù)共得到275 組數(shù)據(jù)，因林分生長量是關(guān)于年齡單調(diào)非減的函數(shù)，需對獲取的數(shù)據(jù)進行篩選，期末的林分因子必須大于或等于期初的，經(jīng)篩選得到170 組數(shù)據(jù)，將其分成兩組獨立樣本，包括建模樣本（115 組）和檢驗樣本（55 組）。對樣本中過大或過小的異常數(shù)據(jù)應(yīng)進行剔除，以提高基礎(chǔ)數(shù)據(jù)的可靠性，保證模型的預(yù)測精度。

建模所需數(shù)據(jù)獲取方法如下。

期初年齡（t1）：根據(jù)龍時勝等[14]基于林木多期直徑數(shù)據(jù)的異齡林年齡估計的方法得到各混交樹種的年齡，再計算得到樣地的平均年齡。期末年齡（t2）：t2=t1+5（間隔期，a）。地位指數(shù)（SI）：利用所編地位指數(shù)表獲取[15]。平均胸徑（D）、每公頃斷面積（G）和每公頃蓄積（M）：對一類清查數(shù)據(jù)進行簡單計算和統(tǒng)計可得。

對樣地各林分因子進行描述性統(tǒng)計（表1），林分的年齡為12～65 a，平均胸徑為6.4～16.8 cm，地位指數(shù)為9.7～22.7 cm；公頃株數(shù)為805～2 385 株，斷面積為6.140～66.990 m2/hm2，蓄積量為15.235～481.845 m3/hm2。

2.2 聯(lián)立方程組模型的構(gòu)建

建立林分生長和收獲模型的方法有很多，通常將其分為三類：全林分模型、徑階分布模型和單木生長模型（分為與距離有關(guān)和無關(guān)的兩種）。本研究擬建立湖南省楠木次生林的可變密度相容性林分生長與收獲模型[16]，根據(jù)Sullivan 和Clutter 提出的方法，假設(shè)3 個基本方程：

從而可以形成以現(xiàn)在林分的一些變量及預(yù)測年齡來預(yù)測未來收獲量的公式，可以表達為生長和收獲的聯(lián)立方程組：

表1 林分因子描述性統(tǒng)計Table 1 The descriptive statistics of stand factors

式（4）～（6）中：M1為期初林分收獲量；SI 為地位指數(shù)；t1為期初林分年齡；G1為期初林分?jǐn)嗝娣e；G2為期末林分?jǐn)嗝娣e；t2為期末林分年齡；M2為期末林分蓄積；b1、b2、b3、b4、b5、b6為待定參數(shù)。

在聯(lián)立方程組中，因變量lnM1在另一個方程中以自變量的形式出現(xiàn)，即有些因變量既是自變量又是因變量。因此，在聯(lián)立方程組中，自變量和因變量的劃分不能按常規(guī)方法進行。為了能夠明確劃分，采用內(nèi)生變量（含隨機誤差）和外生變量（不含隨機誤差）來代替通常使用的自變量和因變量。在聯(lián)立方程組中，lnM1、lnM2和lnG2為內(nèi)生變量，SI、t1、t2和lnG1為外生變量。由于聯(lián)立方程組中各方程之間的隨機誤差具有相關(guān)性，其參數(shù)估計不能采用普通的最小二乘法，而應(yīng)采用二步最小二乘法或三步最小二乘法，中國林業(yè)科學(xué)研究院開發(fā)的Forstat 2.0 軟件提供了這些估計方法，本研究采用EView 9.0 軟件對聯(lián)立方程組中的參數(shù)進行估計。

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

本研究考慮在固定效應(yīng)模型的基礎(chǔ)上，增加樣地水平的隨機效應(yīng)，綜合考慮海拔、坡度、土壤等立地因子和溫度、降水、光照等氣象因子對林分各生長量的影響，構(gòu)建基于樣地水平的湖南楠木次生林生長與收獲混合效應(yīng)聯(lián)立方程組模型。構(gòu)建混合效應(yīng)模型一般要考慮以下三個問題：

1）確定模型的固定效應(yīng)參數(shù)和隨機效應(yīng)參數(shù)

將基礎(chǔ)模型參數(shù)的所有組合作為混合參數(shù)進行模型擬合，對所有收斂的組合，通過其赤池信息準(zhǔn)則（AIC）、貝葉斯信息準(zhǔn)則（BIC）的大小，比較不同模擬過程的擬合效果，AIC、BIC 值越小，表明模型的擬合效果越好，為避免模型參數(shù)過多化，對不同參數(shù)個數(shù)的模擬進行似然比檢驗（LRT），當(dāng)P＜0.000 1 時，認(rèn)為模型之間差異達到顯著水平，當(dāng)差異不顯著時，優(yōu)先選取參數(shù)較少的模型。

式（7）～（9）中：n為模型擬合樣本數(shù)；p為參數(shù)個數(shù)；SSE 為模型的殘差平方和；LL1和LL2分別為兩個比較模型的對數(shù)似然值。

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

本研究的林分?jǐn)嗝娣e和蓄積量數(shù)據(jù)在調(diào)查時間上存在自相關(guān)性，而不同的樣地間也存在一定的差異，可通過混合效應(yīng)模型中的樣地內(nèi)方差-協(xié)方差結(jié)構(gòu)對其進行修正。

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

不同樣地內(nèi)各林分變量的生長差異可通過隨機效應(yīng)的方差-協(xié)方差結(jié)構(gòu)來反映。本研究采用廣義正定矩陣來描述隨機效應(yīng)的方差-協(xié)方差結(jié)構(gòu)，以包含3 個隨機參數(shù)的方差-協(xié)方差矩陣為例，其結(jié)構(gòu)為：

2.4 模型評價與檢驗指標(biāo)

用建模時未使用的55 個獨立樣本數(shù)據(jù)對模型的擬合效果進行綜合評價，包括統(tǒng)計檢驗（Stastitical tests）和視圖分析（Visual techniques）。其中視圖分析為繪制模型的預(yù)估值與殘差值之間的散點圖，分析其殘差分布的差異；統(tǒng)計檢驗指標(biāo)包括AIC、BIC、確定系數(shù)（R2）、平均誤差（Bias）、均方根誤差（Root mean square error，RMSE）等。

式中：yi為實際值；為預(yù)測值；為平均預(yù)估值；n為檢驗樣本數(shù)。

3 結(jié)果與分析

3.1 林分生長與收獲模型的建立

利用Eview 9.0 軟件對聯(lián)立方程組的參數(shù)進行估計，從結(jié)果（表2）可以看出聯(lián)立方程組中3 個方程的擬合相關(guān)系數(shù)（R2）均在0.77 以上，殘差平方和（SSE）均較小，擬合效果較好。

表2 參數(shù)估計結(jié)果Table 2 Results of parameter estimation

以參數(shù)估計的結(jié)果，建立湖南省楠木次生林相容性林分生長與收獲模型系統(tǒng)。

林分?jǐn)嗝娣e生長預(yù)測模型為：

相應(yīng)的林分蓄積量生長預(yù)測模型為：

當(dāng)t2=t1=t時（即預(yù)測間隔期為0 a），此時G2=G1=G，可得到與現(xiàn)在林分一致的收獲量方程：

方程（14）～（16）可預(yù)測未來林分的生長和收獲量。以研究區(qū)的第5497 號樣地為例，其現(xiàn)實林分年齡t1=20 a，地位指數(shù)SI=12，每公頃斷面積G1=16.668 m2，利用方程（16）可得到現(xiàn)實林分的蓄積M1，利用方程（14）和方程（15）可以實現(xiàn)對未來林分的斷面積（G2）和蓄積量（M2）的預(yù)測，并分別計算其平均生長量（Zt）、連年生長量（θt）及生長率（Pt），以5 a 為一個齡階，至60 a 的預(yù)測值見表3。

表3 第5497 號樣地生長與收獲預(yù)測值?Table 3 The predicted value of growth and yield of sample No.5497

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

以廣義正定矩陣來描述隨機效應(yīng)的方差-協(xié)方差結(jié)構(gòu)，在此基礎(chǔ)上擬合包含不同隨機參數(shù)個數(shù)的混合效應(yīng)模型，從中選出最優(yōu)模型。模擬結(jié)果（表4）顯示，式（16）共有10 種收斂的模擬，通過對比其AIC、BIC 值，發(fā)現(xiàn)所有含隨機參數(shù)的模擬，其擬合效果均優(yōu)于基礎(chǔ)模型。當(dāng)只有1 個隨機參數(shù)時，有4 種收斂的情況，模擬1 的效果優(yōu)于其他3 個模擬；當(dāng)包含2 個隨機參數(shù)時，也有4 種收斂的情況，模擬6 的擬合效果優(yōu)于其余模擬；當(dāng)考慮3 個隨機參數(shù)時，有2 種收斂的情況，模擬10 的擬合效果優(yōu)于模擬9。為避免參數(shù)過多化，對不同參數(shù)個數(shù)效果最好的模擬進行似然比檢驗，結(jié)果顯示，模擬6的擬合效果顯著優(yōu)于模擬1（P＜0.000 1），而模擬10 與模擬6 之間無顯著差異（P=0.020 8），將參數(shù)更少的模擬6 作為式（14）混合效應(yīng)擬合的最優(yōu)結(jié)果，其隨機參數(shù)為β1、β3；同理，式（14）的隨機參數(shù)為β6，式（15）的隨機參數(shù)為β1、β3、β6。

表4 混合效應(yīng)模型擬合結(jié)果?Table 4 Fitting result of mixed effect model

混合效應(yīng)模型中隨機參數(shù)的模擬過程在林業(yè)統(tǒng)計軟件Forstat 2.0 的非線性混合模型模塊中進行，同時還可計算得到樣地間的方差-協(xié)方差矩陣和隨機效應(yīng)的方差-協(xié)方差矩陣（表5），最終，基于樣地效應(yīng)和誤差結(jié)構(gòu)矩陣的混合效應(yīng)模型形式如下：

式中：D為樣地間隨機效應(yīng)方差-協(xié)方差矩陣；Ri為樣地內(nèi)誤差效應(yīng)方差-協(xié)方差矩陣；β1i、β3i、β6i為隨機效應(yīng)參數(shù)；Γi(ρ)為自相關(guān)矩陣；i為樣地數(shù)。

3.3 模型擬合效果檢驗

相比基礎(chǔ)模型，混合效應(yīng)模型不僅在參數(shù)構(gòu)造和方差組成上有所差異，在各擬合統(tǒng)計量上也有較大的差異，除AIC、BIC 值比基礎(chǔ)模型小外，其平均誤差和均方根誤差值在數(shù)量上均有較大幅度的下降（表5）。其中M1的Bias 值下降了40.53%，RMSE 值下降了52.08%；G2的Bias 值下降了45.00%，而RMSE 值下降幅度較?。?.03%）；M2的Bias 值下降了38.87%，RMSE 值下降了38.04%，表明混合效應(yīng)模型在預(yù)測過程中誤差減小，精度提升。對比兩種模型的確定系數(shù)，發(fā)現(xiàn)混合效應(yīng)模型的R2值均較基礎(chǔ)模型有所提升，這在一定程度上保證了模型的適用性，其中M1和M2的R2提升較大，分別從0.860 提升到0.968、從0.778 提升到0.915，而G2的確定系數(shù)提升較小，其原因可能是由于其基礎(chǔ)模型的R2值已較高（0.907）。

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

為更加直觀地對比基礎(chǔ)模型與混合效應(yīng)模型的擬合效果差異，分別作其殘差分布（圖1），可以看出相比基礎(chǔ)模型（圖1A、C、E），所構(gòu)建的混合效應(yīng)模型其殘差分布（圖1B、D、F）范圍更小，殘差值的分布也更加均勻，沒有出現(xiàn)明顯的不規(guī)則形狀，說明該模型不存在異方差。綜上，混合效應(yīng)模型的預(yù)測精度和適用性均優(yōu)于基礎(chǔ)模型，能夠更加準(zhǔn)確地預(yù)估湖南楠木次生林的生長與收獲情況。

4 結(jié)論與討論

楠木為我國傳統(tǒng)的珍貴樹種，具有極高的經(jīng)濟價值，其天然林的數(shù)量也是日漸減少，因此對楠木次生林提出科學(xué)合理的經(jīng)營規(guī)劃措施迫在眉睫。本研究以湖南省1989—2014年共6 期的一類清查數(shù)據(jù)中楠木為優(yōu)勢樹種或主要樹種的55 塊樣地數(shù)據(jù)為基礎(chǔ)，統(tǒng)計出每塊樣地的平均胸徑（D）、林分?jǐn)嗝娣e（G）和林分蓄積量（M），以龍時勝等[14]基于林木多期直徑數(shù)據(jù)的異齡林年齡估計的方法得到林分的年齡（t），建立湖南省楠木次生林相容性林分生長與收獲模型，對該林分未來一段時間內(nèi)斷面積和蓄積的生長進行預(yù)測。基礎(chǔ)模型雖有著較高的預(yù)測精度，其對M1、G2、M2擬合的確定系數(shù)（R2）分別為0.860、0.907、0.778，預(yù)測精度（P）分別為85.125%、89.785%、98.707%，但仍存在著許多的不足，未能充分考慮樣地和樹種之間的差異。大量研究表明，加入啞變量或混合效應(yīng)能夠顯著提升模型的預(yù)估精度和適用性，而本研究所構(gòu)建的基于樣地效應(yīng)的混合效應(yīng)聯(lián)立方程系統(tǒng)，在模型的確定系數(shù)和預(yù)測精度上均有著一定的提升，除G2（0.913，94.101%）的提升幅度較小外，M1（0.968，99.497%）、M2（0.915，98.707%）的提升均很明顯，在模型預(yù)測的誤差上，同樣較基礎(chǔ)模型有著大幅地降低，3 個林分變量預(yù)測的平均誤差（Bias）和均方根誤差（RMSE）均減小了40%左右，模型的適用性明顯提升。已有學(xué)者的相關(guān)研究表明非線性混合效應(yīng)模型在擬合效果上具有一定的優(yōu)越性[17-29]，近年來，非線性混合效應(yīng)模型已廣泛應(yīng)用于林業(yè)中，主要包括落葉松、高山松、杉木和栓皮櫟等樹種[20-23]，提高了其在樹高、蓄積、生物量和冠幅等模型的預(yù)測精度[24-26]，本研究建立的湖南楠木次生林林分?jǐn)嗝娣e和蓄積生長與收獲模型，優(yōu)化了模型效果，提高了預(yù)測精度，可為該林分的可持續(xù)經(jīng)營提供技術(shù)指導(dǎo)。

圖1 基礎(chǔ)模型與混合效應(yīng)模型殘差分布Fig.1 Residual distributions of basic model and mixed effect model

本研究受數(shù)據(jù)資料的限制，模型預(yù)測的間隔期為5 a，能夠?qū)﹂敬紊侄唐趦?nèi)的生長與收獲進行預(yù)估，但較長時期的預(yù)估尚沒有對其準(zhǔn)確性做過檢驗，因此造成對模型的應(yīng)用以分時間段進行為好，再將各段累加為累積生長量。這樣做對模型的總體性和準(zhǔn)確性可能會有一定影響，其改進的辦法需要增加長期觀測的樣地資料。此外，所構(gòu)建的是相容性模型組，該建模方法保證了生長量模型與收獲量模型之間的相容性，以及未來與現(xiàn)在收獲模型之間的統(tǒng)一性，但在聯(lián)立方程組中，一些因變量在另一個方程中以自變量的形式出現(xiàn)，即這些變量既是因變量又是自變量，這也導(dǎo)致了建模過程中為了滿足方程組整體的預(yù)測精度，而降低了單個模型的預(yù)估精度。本研究對其降低程度沒有進行分析，今后會在此基礎(chǔ)上通過對比其與單個模型的擬合效果，對模型進行優(yōu)化和完善。