含氣候效應(yīng)的湖南杉木人工林?jǐn)嗝娣e生長(zhǎng)模型

（中南林業(yè)科技大學(xué) 林學(xué)院，湖南 長(zhǎng)沙 410004）

斷面積是評(píng)價(jià)森林質(zhì)量高低、計(jì)算收益率的重要指標(biāo)之一[1-2]，有著易于測(cè)量、數(shù)據(jù)準(zhǔn)確等優(yōu)勢(shì)，是制定森林經(jīng)營(yíng)措施、采伐計(jì)劃的依據(jù)。林分?jǐn)嗝娣e作為反映林分生長(zhǎng)過(guò)程的指標(biāo)還影響著林分蓄積[3]，其模型精度更是直接關(guān)系全林整體模型的預(yù)測(cè)精度，因此成為國(guó)內(nèi)外研究的熱點(diǎn)之一[3]。

杉木Cunninghamia lanceolata是我國(guó)南方的主要造林與用材樹(shù)種之一[5]，也是我國(guó)亞熱帶地區(qū)特有的優(yōu)良用材樹(shù)種[6]。第八次全國(guó)森林資源清查表明，杉木人工林現(xiàn)有種植面積約1.24×103萬(wàn)hm2，約占全國(guó)人工林面積的26.6%，在商品用材中的占有率約為1/4，其在滿足國(guó)民經(jīng)濟(jì)發(fā)展和人民對(duì)森林多種效益的需求上具有重要的地位和作用[7]。

林分年齡、立地質(zhì)量及林分密度是構(gòu)建林分?jǐn)嗝娣e模型的變量，目前對(duì)林分?jǐn)嗝娣e模型的模擬在立地質(zhì)量上的處理多數(shù)都局限于地形地貌等因子[8-9]。高東啟等[10]根據(jù)油松林分平均高推算出優(yōu)勢(shì)木平均高，并以此為立地指數(shù)，運(yùn)用啞變量的方法構(gòu)建了油松林分?jǐn)嗝娣e生長(zhǎng)模型；胡松[3]以湖南櫟類天然林為研究對(duì)象，運(yùn)用混合效應(yīng)模型方法，以標(biāo)準(zhǔn)年齡為20 a時(shí)林分優(yōu)勢(shì)高作為立地指數(shù)，構(gòu)建了湖南櫟類天然林林分?jǐn)嗝娣e生長(zhǎng)模型；顏偉等[11]考慮到林地立地復(fù)雜所產(chǎn)生的誤差，選擇不含立地指標(biāo)的基礎(chǔ)模型構(gòu)建了楊樹(shù)與櫟類林分?jǐn)嗝娣e模型。上述模型雖解決了相應(yīng)問(wèn)題，但大多假定氣候不變，無(wú)法考慮在氣候因子影響下的林分生長(zhǎng)。而《第二次氣候變化國(guó)家評(píng)估報(bào)告》顯示，至21世紀(jì)末，我國(guó)年平均溫度、年降水量將分別提升3.5℃、4.2%[12]。IPCC也在第五次評(píng)估報(bào)告中提出，過(guò)去50 a的升溫率幾乎是過(guò)去100 a的兩倍[13-14]，氣候變化速率已顯著加劇?？紤]到杉木對(duì)氣候的敏感性，本研究在使用林分優(yōu)勢(shì)木平均高作為立地質(zhì)量指標(biāo)的基礎(chǔ)上，在后續(xù)建模過(guò)程中添加溫度、降水量等氣候因子，分析了其對(duì)杉木林分?jǐn)嗝娣e生長(zhǎng)的影響，為森林經(jīng)營(yíng)者在氣候影響下采取的杉木經(jīng)營(yíng)措施提供了理論依據(jù)。

1 研究區(qū)概況

湖南地處云貴高原向江南丘陵、南嶺山脈向江漢平原過(guò)渡的地帶，三面環(huán)山，呈朝北開(kāi)口的馬蹄形地貌，位于24°38′～30°08N′，108°47′～114°15′E，東臨江西，西接重慶、貴州，南毗廣東、廣西，北連湖北，總面積21.18萬(wàn)km2。氣候類型為亞熱帶季風(fēng)濕潤(rùn)氣候，春秋氣溫多變，冬寒冷夏酷熱，春夏多雨，秋冬干旱。大部分地區(qū)日平均氣溫穩(wěn)定，年均溫16～19℃，無(wú)霜期253～311 d，全省年均降水量為1 200～1 700 mm，雨量充沛，水熱充足。

2 材料與方法

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

本研究構(gòu)建模型所用數(shù)據(jù)來(lái)源于湖南省各市縣1980年杉木人工林樣地調(diào)查數(shù)據(jù)，樣地大小667 m2左右，主要分布在湖南的東（長(zhǎng)沙、株洲、湘潭）、北（常德、益陽(yáng)）、南（永州、郴州、衡陽(yáng)）。經(jīng)緯度坐標(biāo)、海拔、土壤質(zhì)地、坡位、林分平均胸徑、坡形、坡度、土壤厚度、林分年齡、每公頃株數(shù)、優(yōu)勢(shì)木平均高等為主要的樣地調(diào)查因子，其中所需數(shù)據(jù)完整的樣地共638塊。氣候數(shù)據(jù)是根據(jù)各樣地經(jīng)緯度坐標(biāo)及海拔從Wang等[12]編寫(xiě)的可提取亞太地區(qū)氣候數(shù)據(jù)的ClimateAP (2019)中獲得，共12類氣候因子（表1）。

表1 氣候變量含義說(shuō)明Table 1 Definition of climate variables

為提高建模的精準(zhǔn)性與模型評(píng)價(jià)的合理性，本研究將樣地?cái)?shù)據(jù)按3:1劃分為建模數(shù)據(jù)478組與檢驗(yàn)數(shù)據(jù)160組，用來(lái)模擬與檢驗(yàn)杉木人工林林分?jǐn)嗝娣e生長(zhǎng)模型，相關(guān)統(tǒng)計(jì)量如表2所示。

表2 建模和檢驗(yàn)數(shù)據(jù)統(tǒng)計(jì)?Table 2 Modeling and verifying data statistics

2.2 研究方法

2.2.1 自變量的篩選

本研究使用多元逐步回歸分析對(duì)氣候因子進(jìn)行自變量篩選，根據(jù)方差膨脹因子（VIF）剔除共線性嚴(yán)重的因子，從而保留共線性弱且影響顯著的因子，具體步驟[3]如下：

1）給定顯著性α，其對(duì)應(yīng)臨界值記為F（1）。對(duì)回歸模型中的i個(gè)自變量Xi分別同因變量進(jìn)行一元線性回歸分析。計(jì)算各自變量回歸系數(shù)的F檢驗(yàn)統(tǒng)計(jì)量Fi（1），選其最大值Fi1（1），當(dāng)Fi1（1）≥F（1）時(shí)，則將Xi1引入回歸模型[3]。

2）給定顯著性α，其對(duì)應(yīng)臨界值記為F（2），建立因變量與自變量子集的二元回歸模型，并計(jì)算自變量回歸系數(shù)的F檢驗(yàn)統(tǒng)計(jì)量值Fj（2），選取其中的最大值Fj2（2），且當(dāng)Fj2（2）≥F（2）時(shí)，將Xj2引入回歸模型，否則終止[3]。

3）重復(fù)步驟2，將各自變量逐一代入模型進(jìn)行F檢驗(yàn)，當(dāng)原自變量因后引入的自變量變得不再顯著時(shí)，將其剔除，以確保模型中只含顯著性變量[3]。

2.2.2 林分?jǐn)嗝娣e模型的構(gòu)建

林分生長(zhǎng)發(fā)育時(shí)，深受其林分年齡、立地質(zhì)量以及對(duì)林分資源利用程度的影響[15-16]，因此構(gòu)建林分?jǐn)嗝娣e生長(zhǎng)模型中需包含立地質(zhì)量指標(biāo)、密度指標(biāo)與林分年齡3種變量[17]。本研究選擇常用的4種模型作為候選基礎(chǔ)模型，在基礎(chǔ)模型的構(gòu)建中，立地質(zhì)量指標(biāo)選用林分優(yōu)勢(shì)木平均高HD。密度指標(biāo)在對(duì)比各基礎(chǔ)模型精度差異后，從林分密度指數(shù)SDI與每公頃株數(shù)N中擇優(yōu)選擇[7]，具體模型形式如下：

式中：T、BA、SDI、HD分別為林分年齡、林分?jǐn)嗝娣e、林分密度指數(shù)、林分優(yōu)勢(shì)木平均高；a、b、c、d、f、g均為模型固定參數(shù)。

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

依據(jù)回歸函數(shù)同時(shí)依賴于固定效應(yīng)參數(shù)和隨機(jī)效應(yīng)參數(shù)的回歸關(guān)系而建立的模型稱為混合效應(yīng)模型（Mixed effects model），其一般形式[3,18]如下：

式中：yi與xi分別為第i個(gè)樣地的因變量向量和自變量向量；εi為誤差項(xiàng)；β與μi分別為固定效應(yīng)參數(shù)向量和隨機(jī)參數(shù)向量。

在構(gòu)建混合效應(yīng)模型的過(guò)程中，關(guān)鍵步驟是模型兩大效應(yīng)——隨機(jī)效應(yīng)與固定效應(yīng)的參數(shù)構(gòu)造，即將所有與研究對(duì)象有關(guān)系的自變量以排列組合的形式添加到模型各參數(shù)上作為隨機(jī)效應(yīng)擬合，根據(jù)AIC、BIC對(duì)模型進(jìn)行評(píng)價(jià)，AIC、BIC值越小，模型的擬合效果越好。然而過(guò)多的自變量以及參數(shù)均會(huì)造成模型的不收斂[19-20]。為避免該問(wèn)題，選取出參數(shù)較少且收斂形式的參數(shù)構(gòu)造，本研究先進(jìn)行研究對(duì)象的顯著性因子篩選，再利用最優(yōu)基礎(chǔ)模型進(jìn)行多參數(shù)效應(yīng)模擬檢驗(yàn)[3]。

2.3 模型精度評(píng)價(jià)

本研究使用調(diào)整決定系數(shù)Ra2、均方根誤差RMSE、平均相對(duì)誤差絕對(duì)值MARE、赤池信息量AIC、貝葉斯信息量BIC進(jìn)行模型評(píng)價(jià)，表達(dá)式如下：

式中：yi為第i樣本實(shí)測(cè)值；是第i樣本預(yù)估值；為平均實(shí)測(cè)值；p為模型中參數(shù)的個(gè)數(shù)；n為樣本數(shù)；l為模型極大似然函數(shù)值。

3 結(jié)果與分析

3.1 氣候因子篩選及分級(jí)

選取年均溫（TMA）、年積溫（DD5）、最熱月均溫（TMWM）、最冷月均溫（TMCM）、均溫差（TD）、年降水量（PMA）、干燥指數(shù)（AHM）、哈格里夫斯氣候水汽虧損（CMD）、夏季平均氣溫（Tave）、夏季平均最高溫（Tmax）、夏季平均最低溫（Tmin）、夏季平均降水量（PPT）等12個(gè)影響杉木生長(zhǎng)的氣候因子，為避免定量因子間嚴(yán)重多重共線性問(wèn)題，使用SPSS 22.0軟件中的多元逐步回歸對(duì)其進(jìn)行分析，利用方差擴(kuò)大因子剔除共線性嚴(yán)重（VIF＞5）的自變量[21]。由表3可知，共線性弱且貢獻(xiàn)度大的氣候因子為T(mén)MCM、PPT、Tmax、CMD，其對(duì)應(yīng)的標(biāo)準(zhǔn)化系數(shù)分別為0.268、0.135、-0.597和-0.322。由此可得TMCM、PPT與斷面積生長(zhǎng)呈正相關(guān)，而Tmax、CMD與斷面積生長(zhǎng)呈負(fù)相關(guān)。

表3 氣候因子多元逐步回歸分析結(jié)果Table 3 Results of multiple stepwise regression analysis of climate factors

選取篩選所得的顯著氣候因子，根據(jù)其各自的取值范圍，按比例劃分等級(jí)（表4）。最終最冷月均溫（TMCM）被劃分為8級(jí)，哈格里夫斯氣候水汽虧損（CMD）被劃分為10級(jí)，夏季均降水量（PPT）被劃分為9級(jí)，夏季均最高溫（Tmax）被劃分為12級(jí)。

表4 氣候因子等級(jí)劃分Table 4 The division of climatic factors grades

3.2 基礎(chǔ)模型

3.2.1 林分密度指數(shù)的計(jì)算

構(gòu)建林分?jǐn)嗝娣e或蓄積量生長(zhǎng)模型時(shí)，林分密度指數(shù)為常用密度指標(biāo)之一[22]，其表達(dá)式如下：

式中：SDI為林分密度指數(shù)；N為林分每公頃株數(shù)；D0為標(biāo)準(zhǔn)平均胸徑；D為林分平均胸徑；β為自然稀疏率。

為了確定SDI值，必須先估算β。本研究使用二次剔除不足立木度的樣地的方法來(lái)估算β[23]。先取全部樣本，建立回歸方程 lnN=a1-b1lnDg，剔除 lnN＜a1-b1lnDg的樣地。然后用剩余的樣地建立回歸方程 lnN=a2-b2lnDg，再剔除lnN＜a2-b2lnDg的樣地。最后用剩余的樣地建立回歸方程[23]：

取全部樣本數(shù)據(jù)，在Forstat 2.2軟件中按上述步驟對(duì)各回歸模型進(jìn)行非線性擬合，最終得上述回歸方程表達(dá)式如下，其調(diào)整決定系數(shù)Ra2=0.71。

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

3.2.2 模型的選取

用R軟件對(duì)上述基礎(chǔ)模型（1）～（4）進(jìn)行非線性擬合，最終得擬合與精度評(píng)價(jià)結(jié)果見(jiàn)表5。

表5 候選模型的參數(shù)擬合與精度評(píng)價(jià)Table 5 Parameter fitting and precision evaluation of candidate models

模型（1）、模型（3）各指標(biāo)均優(yōu)于模型（2）、模型（4），說(shuō)明相比用每公頃株數(shù)去表示密度指標(biāo)，林分密度指數(shù)會(huì)更合適。其中，模型（1）模擬效果最優(yōu)秀，其建模數(shù)據(jù)Ra2最高，各誤差指標(biāo)最低，故選擇模型（1）為模擬林分?jǐn)嗝娣e生長(zhǎng)的基礎(chǔ)模型。

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

為比較分析4種顯著的氣候因子各自差異及其綜合差異對(duì)林分?jǐn)嗝娣e生長(zhǎng)的影響，針對(duì)模型（1），構(gòu)建含隨機(jī)效應(yīng)的林分?jǐn)嗝娣e生長(zhǎng)模型。

運(yùn)用R軟件中的混合模型模塊，分別以篩選所得無(wú)嚴(yán)重共線性且對(duì)斷面積影響顯著的4個(gè)氣候因子TMCM、PPT、Tmax、CMD及其組合形式為隨機(jī)效應(yīng)，并將該效應(yīng)分別添至參數(shù)a、b、c、d、f及其組合形式（1 150種）上進(jìn)行隨機(jī)效應(yīng)模擬，剔除不收斂的類型，根據(jù)AIC、BIC對(duì)各模型形式進(jìn)行評(píng)價(jià)，從而選出最優(yōu)隨機(jī)效應(yīng)構(gòu)造形式，結(jié)果見(jiàn)表6。

表6 氣候因子隨機(jī)效應(yīng)精度評(píng)價(jià)結(jié)果?Table 6 Evaluation results of the random effects of climatic factors

由表6可知，將Tmax添加至參數(shù)a、CMD添加至參數(shù)b的隨機(jī)效應(yīng)構(gòu)造方式構(gòu)建的混合效應(yīng)模型效果最佳。因此，最終構(gòu)建的含氣候隨機(jī)效應(yīng)的林分?jǐn)嗝娣e生長(zhǎng)模型表達(dá)式如下：

式中：BAij、HDij、Tij、SDIij分別表示第i級(jí)夏季均最高溫（Tmax）、第j級(jí)哈格里夫斯氣候水汽虧損（CMD）的林分?jǐn)嗝娣e、優(yōu)勢(shì)木平均高、林分年齡與林分密度指數(shù)，ai0、bi0分別為T(mén)max與CMD的隨機(jī)效應(yīng)參數(shù)，a、b、c、d、f為模型固定參數(shù)。

3.4 模型模擬分析

使用R軟件對(duì)模型（14）進(jìn)行非線性混合效應(yīng)模擬，所得參數(shù)擬合結(jié)果與模型檢驗(yàn)結(jié)果見(jiàn)表7。

表7 模型參數(shù)擬合結(jié)果Table 7 Fitting results of model parameters

由表7可得，各模型參數(shù)擬合結(jié)果均為顯著。對(duì)比AIC、BIC可知，加入隨機(jī)效應(yīng)后的模型（14）擬合效果優(yōu)于基礎(chǔ)模型（1），隨機(jī)效應(yīng)效果明顯。從Ra2、RMSE、MARE來(lái)看，加入隨機(jī)效應(yīng)后的模型（14）在建模與檢驗(yàn)精度上均優(yōu)于基礎(chǔ)模型（1）。其建模樣本的調(diào)整決定系數(shù)Ra2從0.835 5提高至0.892 1，增幅為6.77%，均方根誤差RMSE降低了19.04%，平均相對(duì)誤差絕對(duì)值MARE降低了15.95%；檢驗(yàn)樣本均方根誤差RMSE降低了17.15%，平均相對(duì)誤差絕對(duì)值MARE降低了12.33%。

4 結(jié)論與討論

4.1 討 論

研究林分?jǐn)嗝娣e生長(zhǎng)預(yù)估模型對(duì)預(yù)測(cè)森林生長(zhǎng)與收獲、統(tǒng)籌森林全局經(jīng)營(yíng)具有重要意義，而林分年齡、立地質(zhì)量指標(biāo)及密度指標(biāo)的表達(dá)是構(gòu)建林分?jǐn)嗝娣e生長(zhǎng)模型的必要部分。本研究運(yùn)用混合模型將各氣候因子以隨機(jī)效應(yīng)的形式添加至基礎(chǔ)模型，以分析各水平氣候因子對(duì)林分?jǐn)嗝娣e生長(zhǎng)的隨機(jī)影響。胡松[3]用相同方法分析了不同林分類型與立地類型差異對(duì)櫟類林分?jǐn)嗝娣e生長(zhǎng)的影響，雖研究對(duì)象不一致，但模型精度均有顯著提升；李春明等[24]在對(duì)比傳統(tǒng)的回歸模型方法與混合模型方法構(gòu)建落葉松云冷杉林分?jǐn)嗝娣e模型之后，得出混合模型方法精度更高的結(jié)果，說(shuō)明運(yùn)用混合模型方法構(gòu)建林分?jǐn)嗝娣e模型是合理且有效的。

最終添加至模型隨機(jī)效應(yīng)的氣候因子為夏季平均最高溫Tmax與哈格里夫斯氣候水汽虧損CMD，其中Tmax影響斷面積生長(zhǎng)的最大值，CMD影響斷面積的生長(zhǎng)速率。Tmax與CMD的標(biāo)準(zhǔn)化系數(shù)均為負(fù)值，分別為-0.597、-0.322，說(shuō)明斷面積最大值、生長(zhǎng)速率分別與Tmax、CMD呈負(fù)相關(guān)，這和朱安明[25]在氣候因子對(duì)不同種源的杉木樹(shù)輪中的研究結(jié)果相吻合?？赡芤蚝蠟閬啛釒Ъ撅L(fēng)氣候，在生長(zhǎng)季時(shí)溫度普遍較高，此時(shí)溫度不再是樹(shù)木生長(zhǎng)的主要需求，反而過(guò)高溫會(huì)加劇蒸騰作用，且水汽的虧損同樣也加劇林木水分的缺失，導(dǎo)致杉木缺水生長(zhǎng)受阻礙。與此同時(shí)，本研究運(yùn)用多元逐步回歸分析所得的共線性弱且貢獻(xiàn)度大的氣候因子為T(mén)MCM、PPT、Tmax、CMD，這與臧顥[12]在落葉松立地指數(shù)模型上用相同方法篩選所得的顯著氣候因子也有相似之處；呂振剛等[26]在對(duì)蒙古櫟優(yōu)勢(shì)樹(shù)種適宜區(qū)分布的研究中，運(yùn)用刀切法得出最熱月平均氣溫、濕季降水量、年積溫均對(duì)蒙古櫟影響顯著，都說(shuō)明氣候因子對(duì)林木生長(zhǎng)的影響顯著。

為方便建模，構(gòu)建的斷面積模型中僅使用了優(yōu)勢(shì)木平均高表示立地質(zhì)量指標(biāo)，沒(méi)有直接考慮地形、地貌等立地因子的影響。在后續(xù)研究中，可考慮使用地形、地貌數(shù)據(jù)先構(gòu)建杉木人工林立地指數(shù)模型，求得各樣地位指數(shù)，進(jìn)而考慮氣候因子構(gòu)建斷面積模型，或許模型精度會(huì)進(jìn)一步提升。

4.2 結(jié) 論

本研究以湖南地區(qū)的638塊杉木人工林樣地為研究對(duì)象，估算了林分密度指數(shù)SDI（β=-0.960 53）。比較不同基礎(chǔ)斷面積生長(zhǎng)方程的擬合結(jié)果，確定Richards模型效果最優(yōu)（Ra2=0.835 5，RMSE=3.803 5，MARE=11.779 6）。為考慮不同氣候因子的影響，運(yùn)用多元逐步回歸分析篩選了共線性弱且貢獻(xiàn)度大的氣候因子為T(mén)MCM、PPT、Tmax、CMD，其中前兩者與斷面積生長(zhǎng)呈正相關(guān)，后兩者呈負(fù)相關(guān)。運(yùn)用混合效應(yīng)模型的方法構(gòu)建了含氣候隨機(jī)效應(yīng)的杉木人工林林分?jǐn)嗝娣e生長(zhǎng)模型，確定了最優(yōu)隨機(jī)效應(yīng)參數(shù)構(gòu)造形式。相比基礎(chǔ)模型，含氣候隨機(jī)效應(yīng)的模型精度有明顯提升，其建模精度Ra2=0.892 1，提升了6.77%；RMSE=3.079 2，降低了19.04%；MARE=9.901 1，降低了15.95%。研究結(jié)果說(shuō)明了氣候?qū)α址稚L(zhǎng)影響顯著，為在生長(zhǎng)模型中添加了氣候因子的合理性提供了支撐，其在提高精度的同時(shí)，也有利于杉木人工林區(qū)域性的森林經(jīng)營(yíng)。