基于非線性混合效應(yīng)模型的3個針葉樹種削度方程研究

李春明 付 卓

（1.中國林業(yè)科學(xué)研究院資源信息研究所， 北京 100091；2.生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心， 北京 100094）

編制科學(xué)合理的森林經(jīng)營方案，需要編制通用性的材種出材率表[1-2]。削度是描述樹干直徑沿其樹干向上隨樹干直徑位置的升高而逐漸減小變化的程度，是用來說明干形變化的一種指標(biāo)[3]。削度方程可以為森林經(jīng)營者提供樹干任意部位的直徑、樹干任意部分直徑的高度、全樹干材積以及從地面到任意高度的商品材材積[4]。構(gòu)建樹干削度方程已成為編制材種出材率表的首選方法和基礎(chǔ)工作。在歐美等國家，樹干削度方程已經(jīng)逐漸取代材積表和材積方程[5]。因此，構(gòu)建一個精度高的削度方程，可用較少的樣木資料取得較為理想的結(jié)果，以滿足生產(chǎn)上不同材種規(guī)格的需要。削度方程其主要形式有簡單削度方程[6-9]、分段削度方程[10-14]和可變指數(shù)削度方程[15-19]。任何一個削度方程都不可能完滿地描述所有樹種樹干形狀的變化，同時也不會完全適應(yīng)某一樹種的所有林分。

在構(gòu)建削度方程時，通常選擇具有典型代表性樣地中的多株樹木，然后進行樹干解析獲取數(shù)據(jù)。這些樹木由于所處環(huán)境不同，樹木之間其形狀會存在大的差異，另外數(shù)據(jù)存在著重復(fù)測量等特點。在構(gòu)建削度方程時，存在著數(shù)據(jù)和統(tǒng)計上的挑戰(zhàn)?；旌闲?yīng)模型方法能滿足獨立同分布的假設(shè)，能得到最好線性無偏估計，即能反映總體平均水平又能夠體現(xiàn)林分或樹木之間的差異，可以很好的解決削度方程中存在的上述問題。目前，國內(nèi)外很多學(xué)者在構(gòu)建削度方程時，已經(jīng)考慮數(shù)據(jù)間存在的這些問題。例如基于混合效應(yīng)模型方法的簡單削度方程[20]，基于混合效應(yīng)模型方法的分段削度方程[3,14,20-21]，基于混合效應(yīng)模型方法的可變指數(shù)削度方程[22-25]。

本研究以吉林省汪清林業(yè)局96塊云冷杉針闊混交林局級固定樣地為例，構(gòu)建云杉（Picea asperata）、冷杉（Abies fabri）和紅松（Pinus koraiensis）等優(yōu)勢樹種的削度方程。把3個樹種的解析木數(shù)據(jù)分為兩部分，一部分為模擬數(shù)據(jù)，一部分為驗證數(shù)據(jù)（每個樹種選擇5株）。選擇2種常用的削度方程，利用傳統(tǒng)的擬合方法進行擬合，并在此基礎(chǔ)上考慮樣木水平的隨機效應(yīng)，構(gòu)建基于混合效應(yīng)模型方法的樹干削度方程。選擇均方根誤差、絕對平均殘差和確定系數(shù)等指標(biāo)，對基于混合效應(yīng)模型方法的削度方程與傳統(tǒng)的削度方程擬合方法進行精度比較，最后確定理想的削度方程。

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

本研究區(qū)位于吉林省汪清林業(yè)局金溝嶺林場境內(nèi)（43°17′～43°25′N，130°05′～130°20′E）。在1987年設(shè)立了96塊云冷杉針闊混交林樣地（20 m×30 m）。在設(shè)立的云冷杉針闊混交林中，以云杉、冷杉為主要優(yōu)勢樹種，其他主要樹種包括紅松、白樺（Betula platyphylla）、水曲柳（Fraxinus mandschurica）、椴樹（Tilia tuan）和色木槭（Acer mono）等。在每個樣地中按紅松、云杉和冷杉3種樹種分別選擇1株標(biāo)準(zhǔn)木作為解析木。解析木在伐倒后量測其胸徑、樹高和基部年齡。樹木采用2 m 區(qū)分段，其各段分別在伐根0、1.3、3.6、5.6、7.6 m等位置截取圓盤，在最后一段不足2 m時，在高于上一段1 m處截取最后一個圓盤。對每1株樹木分段測量各圓盤的帶皮胸徑和去皮胸徑。由于一些解析木存在數(shù)據(jù)質(zhì)量問題，最后共選擇了82株云杉、95株冷杉和71株紅松解析木數(shù)據(jù)。每個樹種選擇5株樹木作為驗證數(shù)據(jù)，其余作為模擬數(shù)據(jù)。樹干材積采用平均斷面區(qū)分求積式進行計算，3個樹種的胸徑、樹高及年齡等信息見表1。

表1 云杉、冷杉和紅松等3種樹種樹干解析統(tǒng)計Table 1 The statistical results of stem analysis for 3 species including P.asperata,A.fabri and P.koraiensis

2 研究方法

2.1 基礎(chǔ)模型

本研究采用了2個常用的削度方程形式。第1個是Kozak等[6]在1969年提出的簡單多邊形削度方程，具體見式（1），文中簡稱F1。Max等[10]于1976年在Kozak等簡單削度方程基礎(chǔ)上發(fā)展出分段削度方程，具體見式（2），文中簡稱F2。2個削度方程的具體形式如下：

式中：d為樹干h高處的帶皮直徑；D為帶皮胸徑；H為全樹高；h為從地面起算到某一直徑位置的高度。b1、b2、b3、b4為模型待估參數(shù)，a1、a2為樹干下部和上部拐點處的相對高度。e為模型誤差項。

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

非線性混合效應(yīng)模型包括固定效應(yīng)參數(shù)和隨機效應(yīng)參數(shù)兩部分。在利用混合效應(yīng)模型時，候選的削度方程一般表達形式如式（3）。

2.3 模型評價和檢驗指標(biāo)

通常選擇AIC信息準(zhǔn)則（AIC）、BIC信息準(zhǔn)則（BIC）和-2倍對數(shù)似然值（-2logL）等3個指標(biāo)來比較模型間的擬合效果，評估模型對特定數(shù)據(jù)集的擬合優(yōu)度，表示模型預(yù)測值與數(shù)據(jù)集中觀測值之間的距離。指標(biāo)值越小，說明模型的模擬效果越好[26]。選擇似然比卡方檢驗（LRT）來比較模型之間的差異程度。

在對模型進行精度驗證時，選擇確定系數(shù)（R2）、均方根誤差（RMSE）和絕對平均殘差（）進行評價，如式（4）～（6）。

式中：esti為直徑的預(yù)測值，obji為直徑的測量值，為實測直徑的平均值。N為所有解析木的總體圓盤數(shù)量，m為解析木數(shù)量，n為某一樣木的圓盤數(shù)。

2.4 模型驗證

在對模型進行驗證時，傳統(tǒng)方法需要根據(jù)估計出的參數(shù)來計算驗證數(shù)據(jù)的均方根誤差、絕對平均殘差和確定系數(shù)。而混合效應(yīng)模型需要計算驗證數(shù)據(jù)中具體解析木的隨機效應(yīng)參數(shù)值，然后再計算均方根誤差、絕對平均殘差和確定系數(shù)。隨機效應(yīng)參數(shù)的計算方法可參考Vonesh等[27]，具體的計算公式如式（7）。

3 結(jié)果與分析

3.1 模擬結(jié)果

利用SAS的NLMIXED模塊對Kozak 等[6]和Max等[10]的2個削度方程進行參數(shù)的擬合，具體結(jié)果見表2。利用確定系數(shù)、均方根誤差和絕對平均殘差計算誤差，并利用LRT指標(biāo)對傳統(tǒng)削度方程模擬結(jié)果和基于非線性混合效應(yīng)模型模擬結(jié)果進行比較，具體結(jié)果見表3。其中，Kozak等[6]的簡單削度方程簡稱M1，Max等[10]的分段削度方程簡稱M2，Kozak等[6]的簡單削度方程考慮隨機效應(yīng)的方程簡稱M1a，Max等[10]的分段削度方程考慮隨機效應(yīng)的方程簡稱M2a。

由表2～3可知，這2個削度方程均顯示出較好的模擬效果，確定系數(shù)最小的為0.961 0，最大的為0.991 2。無論是云杉、冷杉還是紅松，F(xiàn)2計算出來的確定系數(shù)都要高于F1，而均方根誤差和絕對平均殘差都要小于F1，這說明F2模擬效果要好于F1。這2個削度方程，無論擬合哪個樹種，當(dāng)考慮單木水平的隨機效應(yīng)后，均只有1個隨機參數(shù)收斂，考慮2個或2個以上參數(shù)時，模型均不能夠收斂。F1在參數(shù)b1上考慮隨機效應(yīng)的模擬效果要好于b2上考慮隨機效應(yīng)。F2在3個樹種上的表現(xiàn)各不相同，在模擬紅松時，b1的模擬效果最好。在模擬冷杉時，b4的模擬效果最好。在模擬云杉時，b2的模擬效果最好。所有3個樹種，無論F1還是F2，在考慮隨機效應(yīng)后模擬效果均有所提高（AIC、BIC和-2logL等3個值降低），LRT值和p值 顯示效果均達顯著程度（α=0.05），并且確定系數(shù)、均方根誤差和絕對平均殘差也支持AIC、BIC和-2logL所獲的結(jié)論。

表2 云杉、冷杉和紅松3種樹種不同削度方程的參數(shù)模擬結(jié)果Table 2 The parameters simulation results of different taper equations for 3 species including P.asperata,A.fabri and P.koraiensis

表3 云杉、冷杉和紅松3種樹種不同削度方程的模擬效果比較Table 3 The comparison of accuracy of different taper equations for 3 species including P.asperata,A.fabri and P.koraiensis

根據(jù)模擬結(jié)果，以相對高（樹干任意高度/樹高）為橫坐標(biāo)，相對直徑（樹干任意高度/樹高）為縱坐標(biāo)，選擇最優(yōu)的削度方程，繪制3個樹種的樹干曲線，具體見圖1。很顯然，對于3個樹種來說，分段削度方程擬合效果很好。

圖1 云杉、冷杉和紅松3種樹種最優(yōu)削度方程的樹干曲線圖Fig.1 Relative height plotted against relative diameter with a segment taper equation of 3 species including P.asperata,A.fabri and P.koraiensis

驗證結(jié)果表明，與不考慮隨機效應(yīng)相比，當(dāng)考慮隨機效應(yīng)后，除了紅松的F2的絕對平均殘差增大外，其余的確定系數(shù)均提高，而均方根誤差和絕對平均殘差都降低，符合上述的模擬結(jié)論。這3個樹種無論哪個樹種，F(xiàn)2的驗證精度都要高于F1的驗證精度，與模擬結(jié)論也一致。

3.2 模型驗證

為了驗證上述模擬結(jié)果，利用確定系數(shù)（R2）、均方根誤差（RMSE）和絕對平均殘差（）3個模型精度評價指標(biāo)對驗證數(shù)據(jù)進行計算。具體結(jié)果見表4。其中不考慮隨機效應(yīng)的模型（M1、M2）直接利用表2的參數(shù)代入公式結(jié)果。而考慮隨機效應(yīng)的模型（M1a、M2a），首先選擇5株驗證數(shù)據(jù)來估計每個樣木的隨機效應(yīng)值，然后再計算3個評價指標(biāo)。隨機效應(yīng)參數(shù)值的計算方法參考式（7）。

表4 云杉、冷杉和紅松3種樹種2種削度方程的驗證結(jié)果Table 4 The validation result of 2 taper equations for 3 species including P.asperata,A.fabri and P.koraiensis

4 結(jié)論與討論

本研究以云冷杉針闊混交林為例，選擇云杉、冷杉和紅松等3個主要樹種的樹干解析數(shù)據(jù)，采用了2種常用的削度方程形式，基于混合效應(yīng)模型方法來構(gòu)建削度方程。由于每個樣地只選擇了1株標(biāo)準(zhǔn)木，因此只選擇了樹木效應(yīng)，將不同隨機效應(yīng)參數(shù)進行組合然后模擬，利用AIC、BIC和-2logL 3個指標(biāo)來評價混合效應(yīng)模型與傳統(tǒng)方法的效果。結(jié)果表明：在2個以上參數(shù)進行隨機效應(yīng)組合時，模型均不能夠收斂。這3個樹種無論考慮哪種削度方程，在考慮混合效應(yīng)模型方法后，AIC、BIC和-2logL 3個指標(biāo)值均明顯下降，說明混合模型的模擬效果都比傳統(tǒng)模型的模擬精度高，利用確定系數(shù)、均方根誤差及絕對平均殘差 3個指標(biāo)進行比較，無論是模擬數(shù)據(jù)還是驗證數(shù)據(jù)均支持此結(jié)論。

1株樹木從基部到頂部，可能由凹曲線體、圓柱體、拋物線體、圓錐體等一部分或幾部分組成。因此，這3個樹種，無論考慮隨機效應(yīng)與否，考慮2個拐點的Max等[10]的分段削度方程在模擬和驗證中的精度都要高于Kozak等[6]的簡單削度方程的模擬和驗證精度。

本研究限于樣本量，沒有考慮樣地水平的隨機效應(yīng)，也沒有同時考慮兩層次隨機效應(yīng)。另外本研究考慮了樣木多次測量的異方差和自相關(guān)，但模型均不能夠收斂，可能受限于樣本量的大小，今后如果數(shù)據(jù)滿足要求，可增加內(nèi)容研究。另外樹木的形狀和林分密度、經(jīng)營措施、樹冠大小，甚至氣候變化都有一定的關(guān)系，隨著數(shù)據(jù)積累的增多，需要增加這方面的研究。