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

李春明 付 卓

（1.中國林業科學研究院資源信息研究所， 北京 100091；2.生態環境部衛星環境應用中心， 北京 100094）

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

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

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

1 數據來源

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

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

2 研究方法

2.1 基礎模型

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

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

2.2 非線性混合效應模型

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

2.3 模型評價和檢驗指標

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

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

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

2.4 模型驗證

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

3 結果與分析

3.1 模擬結果

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

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

表2 云杉、冷杉和紅松3種樹種不同削度方程的參數模擬結果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

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

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

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

3.2 模型驗證

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

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

4 結論與討論

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

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

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