包含密度因子的杉木可變指數(shù)削度方程研制

張森森，段愛國,2★，張建國,2，孫建軍

（1.中國林業(yè)科學(xué)研究院林業(yè)研究所·國家林業(yè)和草原局林木培育重點實驗室，北京100091；2.南京林業(yè)大學(xué)·南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心，江蘇 南京210037；3.中國林業(yè)科學(xué)研究院·亞熱帶林業(yè)實驗中心，江西 分宜336600）

樹干削度方程可以較好地描述樹干直徑隨著樹干高度變化的變異規(guī)律，精準(zhǔn)估算樹干商品材材積、總材積與林分蓄積，對經(jīng)營管理森林，評估森林質(zhì)量等具有重要實踐意義和指導(dǎo)作用?？勺冎笖?shù)削度方程[1-4]是通過連續(xù)函數(shù)中自變量指數(shù)的變化來更好地估計樹木干形，理論上可描述任何形狀的樹干?？勺冎笖?shù)削度方程已證明比簡單削度方程和分段式削度方程具有更優(yōu)的擬合效果及應(yīng)用前景[1,4-6]。Kozak于2004年建立了一個結(jié)構(gòu)簡單且能高精度預(yù)估樹干商品材的可變指數(shù)削度方程，并被普遍應(yīng)用[7]。

除樹種特性外，樹干削度受年齡、立地、間伐，密度等林分經(jīng)營因子以及冠長、樹冠高度、枝下高等樹干特征變量的明顯影響，其中林分密度對樹干削度的影響非常大[4,8-9]，在削度方程中加入林分密度變量可解釋密度對樹干干形形成的作用，而相對于其他林分變量因子，林分密度因子更容易獲取。將林分密度變量構(gòu)建入可變指數(shù)削度方程后，發(fā)現(xiàn)可提高黑云杉（Picea mariana）樹干干形模擬精度，且擬合結(jié)果能解釋高密度林分中林木干形削度小、低密度林分中林木干形削度大的現(xiàn)象[4,10]。在研究落葉松（Larix gmelinii）的削度方程時，將冠長率和密度因子構(gòu)建到Max和Burkhart分段削度方程中[11]，結(jié)果表明包含冠長率和密度因子的削度方程可以提高模型的預(yù)測精度，而且密度越大，單木冠長率越小，樹干的削度就越小[12]。

杉木（Cunninghamia lanceolata）是我國南方最重要的鄉(xiāng)土針葉用材樹種[13-14]，其削度方程的構(gòu)建一直是研究熱點問題，但仍缺乏可變指數(shù)削度方程，尤其是引入密度因子的可變指數(shù)削度方程的研制[15-16]。鑒此，以成熟齡杉木密度試驗林內(nèi)的樹干解析資料為基礎(chǔ)，構(gòu)建了杉木可變密度可變指數(shù)削度方程，以精準(zhǔn)預(yù)估初植密度變化條件下杉木商品材材積，為杉木林高效培育提供科學(xué)定性及定量依據(jù)。

1 材料與方法

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

試驗林布設(shè)于江西省分宜縣大崗山年珠林場，地處羅霄山脈北端的武功山支脈，位于27°34′N，114°33′E，低山，海拔250 m，母巖為砂頁巖，黃棕壤，年平均氣溫16.8℃，降水量1 656 mm，年蒸發(fā)量1 503 mm。采用1年生實生苗于1981年春造林，設(shè)置5種造林密度處理，每種造林密度重復(fù)3次，完全隨機(jī)區(qū)組設(shè)計，同一區(qū)組坡向和坡位基本一致，林分保持自然生長狀態(tài)，不進(jìn)行間伐。2008年冰凍雪壓造成其中4種相對高的造林密度林分中發(fā)生大量倒木（表1），對不同造林密度中的183株倒木進(jìn)行了樹干解析，并測量樹干在0.2 m、1 m、1.3 m、2 m和2 m以上每間隔1 m處的直徑，將183株倒木按密度劃分，每個密度中按徑階隨機(jī)分為建模數(shù)據(jù)和檢驗數(shù)據(jù)。不同造林密度林分的建模數(shù)據(jù)和檢驗數(shù)據(jù)的樹高及胸徑數(shù)據(jù)統(tǒng)計量如表1所示。

表1杉木密度試驗林建模數(shù)據(jù)和檢驗數(shù)據(jù)概況Tab.1 Summary statistics for fit data and validation data of density test stand of Chinese fir

1.2 基礎(chǔ)模型

選擇Kozak可變指數(shù)削度方程為基礎(chǔ)模型，模型表達(dá)式為：

1.3 密度因子構(gòu)建

在Kozak削度方程中，描述樹干直徑隨著樹高的變化而變化的部分是削度方程中的指數(shù)部分，因此密度因子可以添加在削度方程的指數(shù)部分，來反映密度因子對樹干干形的影響[4]。為了檢驗密度因子對樹干削度的影響，設(shè)置了8種密度因子形式并引入削度方程進(jìn)行擬合，結(jié)果發(fā)現(xiàn)其中5種收斂，3種不收斂，5種 收 斂 的 表 達(dá) 形 式 為

1.4 混合效應(yīng)模型建立

考慮樣木效應(yīng)建立混合效應(yīng)模型時，由于過多的隨機(jī)參數(shù)，可能會導(dǎo)致過度參數(shù)化或收斂的問題[17]，對選取的基礎(chǔ)模型（1）中1～2個隨機(jī)參數(shù)組合的可能性都進(jìn)行擬合，通過比較統(tǒng)計量AIC(Akaike information criteria)、BIC(Bayesian information criteria)和-2LL的最低值來決定最優(yōu)的固定參數(shù)和隨機(jī)參數(shù)的組合。

1.5 評價標(biāo)準(zhǔn)

式中：yi為觀測值，y^i為預(yù)測值，yˉi為觀測值的平均值，n為觀察值數(shù)量，p為模型參數(shù)個數(shù)。

2 結(jié)果與分析

2.1 包含密度因子的削度方程構(gòu)建

表2列出了5種密度因子表達(dá)形式的參數(shù)估計和擬合統(tǒng)計量。發(fā)現(xiàn)方程1中參數(shù)b10在b10/sd和b10/sd2形式中的標(biāo)準(zhǔn)差過大，分別為25.8784和63296.2。參數(shù)的標(biāo)準(zhǔn)差過大說明參數(shù)不穩(wěn)定，故而排除這兩種形式。在剩下的其他3種形式中(b10/ √s d ，b10/ √3s d，b10/log(sd))，參數(shù)b10的標(biāo)準(zhǔn)差均較小，分別為0.7066、0.2510和0.3430，說明參數(shù)的數(shù)據(jù)整齊穩(wěn)定，且由參數(shù)的標(biāo)準(zhǔn)誤差可知，各個參數(shù)在0.05水平上顯著；這3種形式的擬合統(tǒng)計量基本一致，R2adj、Bias、RMSE分別為0.9838、0.0040～0.0044、0.4450～0.4453，而由于b10/ √s d 相對于其他兩種形式更為簡潔，容易求算，故選擇b10/√s d 為最優(yōu)的包含密度的指數(shù)形式。引入密度因子后的Kozak可變指數(shù)削度方程表達(dá)式如下：

式中各參數(shù)同式1，sd為林分的初植密度(株·hm-2)，b10為估計參數(shù)。

2.2 包含密度因子的混合效應(yīng)模型建立

以含有初植密度變量的削度方程為基礎(chǔ)模型建立基于樣木效應(yīng)的混合效應(yīng)模型，首先確定隨機(jī)參數(shù)。在式6中的10個參數(shù)中選擇1～2個不同隨機(jī)參數(shù)的所有組合，利用S-Plus軟件的nlme模塊進(jìn)行擬合，結(jié)果見表3。表3顯示了隨機(jī)參數(shù)的所有收斂組合(共44對組合)的AIC、BIC和-2LL的值，其中b6、b8隨機(jī)參數(shù)組合的值最小，因此選擇b6、b8為基于樣木效應(yīng)的混合效應(yīng)模型的隨機(jī)參數(shù)。

利用S-Plus軟件分別對不包含密度變量和包含密度變量的Kozak削度方程用非線性回歸和混合效應(yīng)模型兩種方法分別進(jìn)行擬合，共4種不同形式，分別為不包含密度削度方程的非線性回歸擬合(NLS)，包含密度削度方程的非線性回歸擬合(NLSSD)，不包含密度變量削度方程的非線性混合效應(yīng)擬合(NLME)，包含密度變量削度方程的非線性混合效應(yīng)擬合(NLMESD)。表4顯示了4種不同形式的基于樣木效應(yīng)的Kozak削度方程的參數(shù)估計值和擬合統(tǒng)計量（括號內(nèi)為標(biāo)準(zhǔn)差）。其中，NLSSD模型的調(diào)整決定系數(shù)相對于NLS模型的調(diào)整決定系數(shù)提高了0.0001，NLSSD模型的平均偏差（Bias=0.0037）相對于NLS模型的平均偏差（Bias=0.0041）降低了0.0004，NLSSD模型的均方根誤差（RMSE=0.4450）相對于NLS模型的均方根誤差（RMSE=0.4470）降低了0.0020，NLSSD模型的平均絕對偏差（MAD=0.3359）相對于NLS模型的平均絕對偏差（MAD=0.3379）降低了0.0020。而NLMESD模型相對于NLME模型的調(diào)整決定系數(shù)（），平均偏差（Bias）和均方根誤差（RMSE）都無顯著變化，均分別為0.9945，0.0006，0.2603，但NLMESD模型的平均絕對偏差（MAD=0.3359）相對于NLME模型的平均絕對偏差（MAD=0.3379）降低了0.0020?？梢园l(fā)現(xiàn)，引入密度因子后，無論是非線性回歸還是混合效應(yīng)模型，削度方程的擬合精度均略有提高，比較而言，采用混合效應(yīng)模型建立的削度方程較非線性回歸模型具有更高的調(diào)整決定系數(shù)及更低的平均偏差、均方根誤差和平均絕對偏差。

圖1為建模數(shù)據(jù)中的不同密度林分平均木，其中包括密度為3 333株·hm-2（D=13 cm，H=13.8 m）、5 000株·hm-2（D=11 cm，H=12.7 m）、6 667株·hm-2（D=10.4 cm，H=12.2 m）和10 000株·hm-2（D=9.7cm，H=11.7 m），用包含密度因子的基于樣木效應(yīng)的Kozak混合效應(yīng)模型擬合的模型輪廓圖。發(fā)現(xiàn)不同密度林分平均木不同高度處直徑的模擬值與實測值十分接近，且當(dāng)林分密度為3 333株·hm-2時，樹干削度最大，而當(dāng)林分密度為10 000株·hm-2時，樹干削度最小，總體上，樹干削度隨著密度的增大而逐漸降低，降低的程度也隨密度的增大而減小。

圖1不同密度林分平均木模擬樹干輪廓圖Fig.1 Taper profiles for mean sample tree of different density with tree effects NLME model using fit data

2.3 包含密度因子的模型驗證

表5列出了60株檢驗樣本擬合不同形式下Kozak可變指數(shù)削度方程的統(tǒng)計量。在非線性回歸擬合中，包含密度因子的削度方程的調(diào)整決定系數(shù)（=0.9804）高于不包含密度因子的削度方程的調(diào)整決定系數(shù)（=0.9801），且包含密度因子的削度方程的平均偏差（Bias=-0.0374）、均方根誤差（RMSE=0.5063）、平均絕對偏差（MAD=0.3527）分別比不包含密度因子的削度方程的平均偏差（Bias=-0.0473）、均方根誤差（0.5104）、平均絕對偏差（0.3537）有所減小。在混合效應(yīng)模型中，包含密度因子的削度方程的調(diào)整決定系數(shù)（=0.9888）高于不包含密度因子的削度方程的調(diào)整決定系數(shù)（=0.9887），而且包含密度因子的削度方程的平均偏差（Bias=0.0184）、均方根誤差（RMSE=0.3828）、平均絕對偏差（MAD=0.2534）分別比不包含密度因子削度方程的平均偏差（Bias=-0.0271）、均方根誤差（0.4401）、平均絕對偏差（0.2568）有所減小。因此，在可變指數(shù)削度方程中構(gòu)建密度因子，是可以提高模型的預(yù)估精度，降低模型預(yù)估偏差，并且利用非線性混合效應(yīng)擬合的包含密度因子的Kozak可變指數(shù)削度方程的預(yù)估精度最高，調(diào)整決定系數(shù)達(dá)到0.9888。最終選擇的模型方程式如下式：

表5 4種不同形式的Kozak(2004)模型的檢驗統(tǒng)計量Tab.5 Goodness-of-fit statistics of Kozak(2004)for the four different forms using validation data

式中變量與式（6）相同。

為評價包含密度因子的Kozak模型對分段樹干直徑的預(yù)測能力，將檢驗數(shù)據(jù)中每株樹的相對高度分為10部分。對于每部分直徑的預(yù)測平均偏差和標(biāo)準(zhǔn)差進(jìn)行分段統(tǒng)計（表6）。從表6可以看出，包含密度因子的Kozak基于樣木效應(yīng)的混合效應(yīng)模型在0.1＜h/H≤0.5，0.7＜h/H≤1處的平均偏差為負(fù)值，說明模型預(yù)測值偏高，在樹干中部（0.0＜h/H≤0.1）和頂端（0.5＜h/H≤0.7）處的平均偏差為正值，說明預(yù)測值偏低，而且在樹干的0.9＜h/H≤1.0處為最高值，在0.7＜h/H≤0.8處為最低值。平均絕對偏差在樹干的基部處的值為0.4833，是整個分段最高值，其次是樹干的頂端處的值為0.2921，而在樹干中部的平均絕對偏差較低，說明包含密度因子的Kozak模型在樹干基部和樹干頂端的預(yù)測值偏離真實值程度最差，在樹干中部的預(yù)測值偏離真實值的程度較低。檢驗值的標(biāo)準(zhǔn)差在樹干基部和頂端處存在最高值0.7540和0.6234，而其他部分的標(biāo)準(zhǔn)差值相接近，說明包含密度因子的Kozak模型在樹干基部和頂端處預(yù)測值的離散程度較高，而其他部分預(yù)測值的離散程度較低?？傮w來說，包含密度因子的Kozak模型對整個樹干的預(yù)測能力較高，但是對樹干基部和頂端的預(yù)測精度稍低于樹干中部。

為了評價包含密度因子的Kozak模型對不同密度林分內(nèi)平均木的預(yù)估效果，選取檢驗數(shù)據(jù)中每個密度林分中的平均木，其中包括B密度（D=12.6 cm，H=13.5 m）、C密度（D=10.8 cm，H=13.0 m）、D密度（D=10.2 cm，H=11.6 m）和E密度（D=10.0 cm，H=12.0 m）。用包含密度因子的Kozak基于樣木效應(yīng)的混合效應(yīng)模型對4個密度的平均木進(jìn)行樹干輪廓預(yù)估（圖2）。由圖2可知，不同初植密度林分中平均木的樹干表現(xiàn)出不同的削度。初植密度為3 333株·hm-2的林分平均木的尖削度最大，并且隨著密度的增大，林分平均木的削度在組間減小，減小的程度也隨著密度的增大而逐漸減小，當(dāng)初植密度為10 000株·hm-2時的林分平均木的尖削度最小。

表6考慮樣木效應(yīng)包含密度因子的Kozak混合模型的分段檢驗統(tǒng)計量Tab.6 Examined statistics for Kozak with density of mixed model with tree effects in the different relative height classes

圖2不同密度下平均木的Kozak模型檢驗樹干輪廓Fig.2 Taper profiles for mean sample tree of different density in Kozak with tree effects NLME model using validation data

3 討論與結(jié)論

大多數(shù)的削度方程以胸徑（D）、全樹高（H）、任意樹干高度（h）和相對樹高（h/H）為自變量，然而一些學(xué)者在削度方程中加入其他自變量，如冠長率、冠幅等[2]，研究結(jié)果表明加入其他自變量之后，并沒有對削度方程的擬合精度產(chǎn)生顯著影響，并且模型模擬的精度取決于樹干直徑測量的準(zhǔn)確度。但亦有研究表明，在削度方程中加入林分經(jīng)營因子或其他自變量之后，削度方程的精度有顯著提高[16,18-20]，本研究結(jié)果亦表明林分密度因子引入后擬合精度有所提高。此外，林分密度因子引入后，通過削度方程預(yù)估的不同密度林分內(nèi)林木干形隨密度的增大，削度逐漸減小，且減小的程度也隨著密度的增大而逐漸降低，這與Sharma and Zhang的研究結(jié)果相同[4]，表明在初植密度較大的林分中，由于林木激烈競爭而使樹冠更多地朝向頂端收縮導(dǎo)致樹干的干形更為通直飽滿，而在初植密度較小的林分中，競爭相對較小，樹冠占全樹高比重更大，樹干尖削度則相對較大。

本研究以Kozak削度方程為原型，探討了杉木包含密度因子的可變指數(shù)削度方程的構(gòu)建過程及效果，比較得到了最優(yōu)的包含密度的指數(shù)表達(dá)形式，發(fā)現(xiàn)引入密度因子后，無論是非線性回歸還是混合效應(yīng)模型，削度方程的擬合精度均有所提高，而應(yīng)用混合效應(yīng)模型對包含密度因子的削度方程進(jìn)行預(yù)估時較非線性回歸模型更優(yōu)。構(gòu)建的杉木可變密度可變指數(shù)削度方程在樹干中部預(yù)估值的平均絕對偏差較低，而在樹干基部和頂端的預(yù)測效果較差，且樹干削度隨初植密度的增大而降低，降低的程度隨著密度的增大而減少。所建立的杉木可變密度可變指數(shù)削度方程拓展了杉木干形預(yù)估的靈活性與準(zhǔn)確性，可滿足杉木單株商品材及林分材積的高精度預(yù)測需求，亦能為杉木材積表研制提供可靠途徑。