張中偉，周根苗，易 烜，齊戰(zhàn)濤，呂 勇

（1. 中南林業(yè)科技大學(xué) 林學(xué)院，湖南 長(zhǎng)沙 410004；2. 湖南省農(nóng)林工業(yè)勘察設(shè)計(jì)研究總院，湖南 長(zhǎng)沙 410007；3. 湖南省青羊湖國(guó)有林場(chǎng)，湖南 寧鄉(xiāng) 410600；4. 湘潭市自然資源和規(guī)劃局，湖南 湘潭 411100）

削度是描述樹干直徑沿其樹干向上隨干徑位置的升高而逐漸減小變化程度的指標(biāo)，削度方程是削度的數(shù)學(xué)表達(dá)式[1]，是以樹干上不同部位的樹干直徑為因變量，全樹高、樹干直徑及不同部位時(shí)距地面高等變量為自變量的數(shù)學(xué)函數(shù)，根據(jù)研究方法的不同，削度方程又可分為簡(jiǎn)單削度方程、分段削度方程和可變指數(shù)削度方程[2]。同時(shí)，削度方程能為森林經(jīng)營(yíng)者提供樹干任意部位直徑的高度、全樹干材積以及從地面到任意高度的商品材材積[3]，因此構(gòu)建樹干削度方程已成為編制材種出材率表的首選方法和基礎(chǔ)工作。

樹干的削度受樹種、年齡、林分密度、氣候因素等多種因素的影響[4-6]。目前，大量學(xué)者研究了林分密度因子[7]和林地因子[8]對(duì)樹干削度的影響，而對(duì)氣候因子影響樹干削度的研究較少。但樹木的生長(zhǎng)對(duì)溫度、降水量等氣候因子非常敏感[9]，因此，近年來有學(xué)者開始逐漸研究氣候因子對(duì)樹干削度的影響。Gordon 等[10]在研究氣候因子對(duì)加拿大樟子松Pinus sylvestnisvar.mongolica樹干削度的影響時(shí)，發(fā)現(xiàn)年平均降水量和無霜期均能顯著影響樹干削度的變化，且對(duì)樹干上部削度的影響更大；Liu 等[11]在研究興安落葉松Larix olgensis的樹干削度時(shí)，考慮氣候因子對(duì)樹干干形的影響，將氣候因子添加到削度方程中，建立了含氣候因子的樹干削度方程，使模型的精度有顯著的提升。然而，在以往的杉木Cunninghamia lanceolata樹干削度方程中，缺少代表氣候因子的變量，故無法反映樹干削度隨氣候因子變化而變化的過程。因此，構(gòu)建含氣候因子的杉木樹干削度方程具有非常重要的意義。

隨著數(shù)學(xué)建模方法的發(fā)展，在構(gòu)建削度方程模型時(shí)，混合效應(yīng)模型方法成為了更多學(xué)者的選擇，因?yàn)閭鹘y(tǒng)的回歸分析方法即最小二乘法是假定數(shù)據(jù)間相互獨(dú)立且非異質(zhì)性[12-13]，因此無法體現(xiàn)出不同因子之間的差異對(duì)樹干生長(zhǎng)的影響，從而使模型的精度不夠高，但混合效應(yīng)模型方法卻能很好地解決上述問題[14-16]。如張森森[8]、聶璐毅等[17]使用了混合效應(yīng)模型方法分別構(gòu)建了江西省杉木和長(zhǎng)白落葉松的樹干削度方程，研究結(jié)果表明，使用混合效應(yīng)模型方法所構(gòu)建的模型精度均優(yōu)于傳統(tǒng)方法所構(gòu)建的模型，且誤差均有明顯的降低。

杉木是我國(guó)主要的用材樹種，也是我國(guó)南方人工造林的主要樹種，其具有分布區(qū)域廣、經(jīng)濟(jì)效益高等特點(diǎn)，廣泛應(yīng)用于建筑、家具等領(lǐng)域。本研究以湖南省永州市杉木人工林為研究對(duì)象，使用混合效應(yīng)模型方法構(gòu)建了含氣候因子的杉木樹干削度方程，研究氣候因子對(duì)杉木樹干削度的影響，可為森林經(jīng)營(yíng)者在氣候影響的條件下采取的杉木經(jīng)營(yíng)措施提供理論依據(jù)。

1 材料與方法

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

本研究的數(shù)據(jù)來源于湖南省永州市東安縣、道縣和零陵區(qū)的杉木人工林固定樣地及臨時(shí)樣地（108°47′～114°15′E，24°38′～30°08′N）。采集地區(qū)的氣候?yàn)榇箨懶詠啛釒Ъ撅L(fēng)濕潤(rùn)氣候，雨量充沛，水熱充足，年均溫度為18 ℃，年均降水量為1 500 mm，土壤肥沃，主要樹種有：杉木、馬尾松Pinus massoniana、木荷Schima superba、青岡Cyclobalanopsis glauca等。

本次一共調(diào)查了75 塊樣地的解析木數(shù)據(jù)，樣地大小均為20 m×20 m，其中東安縣有20 塊樣地，道縣有13 塊樣地，零陵區(qū)有42 塊樣地。每塊樣地選擇一株平均木伐倒，按照樹干解析的標(biāo)準(zhǔn)對(duì)其進(jìn)行樹干解析，測(cè)量其樹高（H）、胸徑（D）、樹干上各部位的直徑（d）和該干徑位置距地面的高（h），同時(shí)記錄對(duì)應(yīng)樣木的經(jīng)緯度坐標(biāo)和海拔（HB），共計(jì)550 組數(shù)據(jù)，其因子基本情況詳見表1。

表1 各因子基本情況?Table 1 Basic information table of factors

氣候數(shù)據(jù)是從Wang 等[18]編寫的可獲取亞洲區(qū)域氣候數(shù)據(jù)的ClimateAP（2019）軟件中[19-21]獲得，其中包括年積溫（DD5）、夏季平均最高溫（Txmax）、最冷月平均溫度（TMmin）等共11 類氣候因子。

為保證模型構(gòu)建的合理性，本研究先使用三倍標(biāo)準(zhǔn)差法對(duì)異常數(shù)據(jù)進(jìn)行剔除，剔除后剩476組數(shù)據(jù)，再將這476 組數(shù)據(jù)按照2∶1 的比例分為建模數(shù)據(jù)317 組與檢驗(yàn)數(shù)據(jù)159 組，用以杉木人工林樹干削度方程的構(gòu)建與檢驗(yàn)。

1.2 研究方法

1.2.1 基礎(chǔ)模型的篩選

在三類削度方程中，可變指數(shù)方程有著更好的擬合效果[8]，能更好地用來描述樹干削度的變化[2]。參考前人的研究，本研究選擇了4 個(gè)描述杉木樹干削度較好的可變指數(shù)削度方程作為備選的基礎(chǔ)模型[22-25]，各方程的表達(dá)式如下：

式中：h為不同部位時(shí)距地面高；d為不同部位的樹干直徑；D為樹干胸徑；H為總樹高；b1、b2、b3、b4、b5、b6、b7、b8、b9為參數(shù)；T=h/H，p=1.3/H，

為了后續(xù)表達(dá)方便，將模型（1）稱為M1，模型（2）稱為M2，模型（3）稱為M3，模型（4）稱為M4。

1.2.2 自變量的篩選

本研究使用Forstat 中的數(shù)量化方法Ⅰ對(duì)氣候因子進(jìn)行顯著篩選，以“P＞F”為標(biāo)準(zhǔn)，剔除P＞0.05 的因子，從而保留影響顯著的因子。參考前人的研究，林分密度因子對(duì)樹干削度的影響非常顯著[7-8]，所以本文直接使用林分密度因子作為顯著影響因子。

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

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

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

1.2.4 模型的評(píng)價(jià)指標(biāo)

本研究選用調(diào)整決定系數(shù)Ra2、均方根誤差RMSE、平均絕對(duì)誤差MAE、赤池信息準(zhǔn)則AIC與貝葉斯信息準(zhǔn)則BIC 對(duì)構(gòu)建的模型進(jìn)行評(píng)價(jià)，表達(dá)式如下：

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

2 結(jié)果與分析

2.1 基礎(chǔ)模型的篩選

利用R 軟件的nls 板塊將4 個(gè)備選基礎(chǔ)模型進(jìn)行擬合，以Ra2、MAE 和RMSE 等評(píng)價(jià)指標(biāo)對(duì)各基礎(chǔ)模型進(jìn)行評(píng)價(jià)，結(jié)果詳見表2。

表2 候選基礎(chǔ)模型評(píng)價(jià)Table 2 Candidate basic models

通過擬合，對(duì)比Ra2、MAE、RMSE 這3 個(gè)指標(biāo)，可以得出，在備選的4 個(gè)模型中，模型M3 的各項(xiàng)指標(biāo)均為最優(yōu)，其Ra2最大為0.959 5，RMSE 最小為0.894 7，MAE 最小為0.658 9，很多學(xué)者的研究也得出了一樣的結(jié)論[28-29]，說明Kozak（2004）模型的普適性較好，在可變指數(shù)模型中很有代表性[30]。因此將模型M3 作為后續(xù)研究的基礎(chǔ)模型。

2.2 因子顯著性分析

從獲取的11 個(gè)影響樹干削度的氣候因子中，根據(jù)取值范圍的不同，參考齊戰(zhàn)濤等[31]的劃分標(biāo)準(zhǔn)，對(duì)氣候因子進(jìn)行等級(jí)劃分，結(jié)果詳見表3。

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

使用Forstat 中的數(shù)量化方法Ⅰ對(duì)各已分好級(jí)的氣候因子進(jìn)行顯著性分析，并以“P＞F”為標(biāo)準(zhǔn)，剔除P＞0.05 的因子。由于削度是描述樹干直徑沿其樹干向上隨干徑位置的升高而逐漸減小變化程度的指標(biāo)，為考慮距地面高（h）、樹干胸徑（D）不同部位的樹干直徑（d）的影響，在篩選時(shí)，選取不同部位時(shí)距地面高（h）、樹干胸徑（D）作為輔助變量進(jìn)行分析。各因子分析結(jié)果詳見表4。

由表4 可得，生長(zhǎng)積溫（DD5）、夏季平均溫度（Txave）、夏季平均最高溫（Txmax）為影響樹干上各部位的直徑（d）變化的顯著因子。由于Txave=（Txmax+Txmin）/2，因此夏季平均溫度（Txave）與夏季平均最高溫（Txmax）之間可能存在自相關(guān)、共線性問題，為避免這些問題對(duì)模型的構(gòu)建造成影響，因此本次研究暫時(shí)不考慮夏季平均溫度，僅選取生長(zhǎng)積溫、夏季平均最高溫為顯著影響樹干上各部位的直徑（d）變化的氣候因子。其中生長(zhǎng)積溫的第一個(gè)等級(jí)為3 200 ≤DD5＜3 400；夏季平均最高溫的第一個(gè)等級(jí)為25 ≤Txmax＜26。

2.3 含林分密度因子的非線性混合模型構(gòu)建

針對(duì)選中的基礎(chǔ)模型M3，使用非線性混合效應(yīng)模型方法構(gòu)建含林分密度因子的杉木樹干削度方程，利用R 語言的nlme 包，將林分密度（DM）作為隨機(jī)效應(yīng)分別添加至參數(shù)b1、b2、b3、b4、b5、b6、b7、b8、b9上進(jìn)行擬合，剔除不收斂的組合形式，根據(jù)AIC、BIC 對(duì)剩下的組合形式進(jìn)行評(píng)選，最后篩選出最好的組合形式，結(jié)果詳見表5。

表5 林分密度因子隨機(jī)效應(yīng)精度評(píng)價(jià)結(jié)果Table 5 Evaluation results of the random effects of stand density factors

由表5 可得，當(dāng)林分密度（DM）作為隨機(jī)效應(yīng)添加在參數(shù)b8上時(shí)模型的AIC、BIC 值最低，即效果最佳。為后續(xù)表達(dá)方便，將構(gòu)建的含林分密度因子的削度方程模型簡(jiǎn)稱為M11。

因此，構(gòu)建的含林分密度因子的混合效應(yīng)模型表達(dá)式如下，簡(jiǎn)稱M11：

式中：d、D、H、h、X、T分別為不同部位的樹干直徑、樹干胸徑、總樹高、不同部位時(shí)距地面高、X、T值；a8i為DM的隨機(jī)效應(yīng)參數(shù)；參數(shù)b1、b2、b3、b4、b5、b6、b7、b8、b9為固定參數(shù)。

2.4 含氣候因子和林分密度因子的非線性混合模型構(gòu)建

為比較氣候因子和林分密度因子對(duì)杉木樹干削度的影響，針對(duì)選中的基礎(chǔ)模型M3，使用非線性混合效應(yīng)模型方法構(gòu)建含雙因子的杉木樹干削度方程。

運(yùn)用R 語言的nlme 包，分別以篩選所得的對(duì)樹干削度有顯著影響的2 個(gè)氣候因子夏季平均最高溫、生長(zhǎng)積溫和林分密度及其組合形式作為隨機(jī)效應(yīng)，將其分別添加至參數(shù)b1、b2、b3、b4、b5、b6、b7、b8、b9及其組合形式上進(jìn)行擬合，剔除不收斂的組合形式，根據(jù)AIC、BIC 對(duì)剩下的組合形式進(jìn)行評(píng)選，最后篩選出最好的組合形式，結(jié)果詳見表6。

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

從表6 可得，將Txmax添加到參數(shù)b5、DM添加到參數(shù)b7、DD5添加到參數(shù)b8上時(shí)其AIC、BIC值最低，即效果最好。因此，構(gòu)建的含雙因子的混合效應(yīng)模型表達(dá)式如下，簡(jiǎn)稱M12：

式中：d、D、H、h、X、T分別為不同部位的樹干直徑、樹干胸徑、總樹高、不同部位時(shí)距地面高、X、T值；a5i、a7i、a8i分別為Txmax、DM與DD5的隨機(jī)效應(yīng)參數(shù)；參數(shù)b1、b2、b3、b4、b5、b6、b7、b8、b9為固定參數(shù)。

2.5 模型模擬分析

為分析模型M3、M11、M12 之間擬合表現(xiàn)的優(yōu)劣，利用R 軟件對(duì)模型M3、M11、M12 進(jìn)行模擬，各擬合結(jié)果詳見表7。

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

由表7 可得，模型M12 在AIC、BIC 上的擬合結(jié)果均低于模型M3 和M11，說明加入雙因子的模型M12 的擬合效果要優(yōu)于基礎(chǔ)模型M3 和只加入林分密度因子的M11，擬合效果提升明顯。從Ra2、RMSE、MAE 上的表現(xiàn)來看，模型M12無論是在建模樣本還是檢驗(yàn)樣本上的擬合效果均優(yōu)于模型M3 和M11。其建模樣本調(diào)整后決定系數(shù)Ra2的排序?yàn)镸12 ＞M11 ＞M3，其中M12 相較于M3 的增幅約為2.3%，相較于M11 的增幅約為1.1%；均方根誤差RMSE 的排序?yàn)镸12 ＜M11 ＜M3，其中M12 相較于M3 的降幅約為32.56%，相較于M11 的降幅約為20.31%；平均誤差MAE 的排序?yàn)镸12 ＜M11 ＜M3，其中M12相較于M3 的降幅約為34.56%，相較于M11 的降幅約為21.12%。

為了檢驗(yàn)基礎(chǔ)模型M3 和混合模型M12 是否存在異方差效應(yīng)，建立了模型M12 和M3 的殘差分布（圖1 ～2）。

圖1 模型M12 殘差分布Fig. 1 Residual distribution diagram of model M12

圖2 模型M3 殘差分布Fig. 2 Residual distribution diagram of model M3

由圖1 ～2 可得， 無論是模型M3 還是M12，它們的殘差都較集中且隨機(jī)分布在橫坐標(biāo)軸兩側(cè)，說明模型M3 和M12 都沒有明顯的異方差效應(yīng)，故不用考慮。

3 結(jié)論與討論

3.1 討 論

研究杉木人工林樹干削度方程對(duì)于杉木人工林的生長(zhǎng)與收獲的預(yù)測(cè)和森林經(jīng)營(yíng)方案的編制具有重要意義。本研究使用混合效應(yīng)模型方法構(gòu)建了含雙因子和單含林分密度因子的杉木樹干削度方程，以對(duì)比分析氣候因子對(duì)樹干削度的影響。姜立春等[32]使用了同樣的方法研究了樣地效應(yīng)和樹木效應(yīng)對(duì)落葉松樹干削度的影響，雖然研究的對(duì)象不同，但模型的精度都有顯著提升；李春明等[33]使用混合效應(yīng)模型方法構(gòu)建了紅松、冷杉和云杉的樹干削度方程，同時(shí)將其與傳統(tǒng)擬合方法所構(gòu)建的樹干削度方程進(jìn)行比較，最終得出無論是哪種樹種，混合模型都要優(yōu)于傳統(tǒng)擬合方法所構(gòu)建的模型，因此構(gòu)建削度方程時(shí)使用混合效應(yīng)模型方法非常合理。

最終添至模型中的氣候因子為夏季平均最高溫和生長(zhǎng)積溫，相較于基礎(chǔ)模型，加入氣候因子的混合模型的精度有明顯提升，說明氣候因子對(duì)杉木樹干削度有顯著影響。Liu 等[11]在研究氣候因子對(duì)興安落葉松樹干削度的影響時(shí)，得出影響樹干削度顯著的氣候因子為年平均溫度和年平均降水量，造成這種不同的原因可能是不同樹種對(duì)不同氣候因子的敏感程度不同；Gordon 等[10]在研究加拿大松樹干削度時(shí)，引入氣候因子，發(fā)現(xiàn)氣候?qū)涓缮喜康南鞫扔绊懽畲?，認(rèn)為削度和氣候之間的聯(lián)系是通過樹冠來實(shí)現(xiàn)的，而樹冠又是影響削度變化的一個(gè)非常重要的因子[5,34]，因此氣候因子對(duì)樹干削度的影響非常顯著。

為減少模型的參數(shù)數(shù)量，本研究只考慮了氣候因子和林分密度因子，沒有將樣地因子考慮在內(nèi)。在后續(xù)對(duì)樹干削度的研究中，可以將樣地因子也一并考慮，加入到削度方程中，或許可以使所構(gòu)建模型的擬合精度有進(jìn)一步的提升。

3.2 結(jié) 論

本研究以湖南省杉木人工林的樹干削度為研究對(duì)象，基于永州市75 塊杉木人工林樣地的實(shí)測(cè)數(shù)據(jù)，選取4 個(gè)備選基礎(chǔ)模型，通過擬合，確定Kozak（2004）模型為后續(xù)研究的基礎(chǔ)模型（Ra2=0.959 5，RMSE=0.894 7，MAE=0.658 9）。為研究氣候因子對(duì)樹干削度的影響，通過數(shù)量化方法Ⅰ對(duì)氣候因子進(jìn)行顯著性分析，結(jié)果顯示，夏季平均最高溫和生長(zhǎng)積溫為影響杉木樹干削度顯著的氣候因子。使用非線性混合效應(yīng)模型方法構(gòu)建了含林分密度因子和含雙因子的杉木人工林樹干削度方程，確定了最優(yōu)的隨機(jī)效應(yīng)參數(shù)構(gòu)造形式。從Ra2、RMSE、MAE 上的表現(xiàn)來說，含林分密度因子的模型M11 和含雙因子的模型M12 相較于基礎(chǔ)模型M3 均有明顯提升，其中模型M12的Ra2最大，為0.981 6，相較于M3 的增幅為2.3%，相較于M11 的增幅為1.1%；RMSE 最小為0.603 4，相較于M3 的降幅為32.56%，相較于M11 的降幅為20.31%；MAE 最小為0.431 2，相較于M3 的降幅為34.56%，相較于M11 的降幅為21.12%。上述研究結(jié)果表明，氣候因子對(duì)杉木樹干削度有顯著影響，為今后在研究樹干削度時(shí)考慮氣候因子對(duì)樹干削度的影響提供了理論依據(jù)，其既提高了模型的擬合精度，同時(shí)也有利于杉木人工林森林經(jīng)營(yíng)管理措施的制定和林業(yè)數(shù)表的編制。