金 明，丁貴杰

（貴州大學(xué) 造林生態(tài)研究所，貴州 貴陽 550025）

馬尾松Pinus massoniana是中國主要用材樹種之一，其用材林和人工林面積分別占全國的16.6%和16.9%。在貴州馬尾松的面積和蓄積已分別占到喬木林的26.95%和36.10%[1]。研究和編制馬尾松材種出材率表在指導(dǎo)生產(chǎn)實踐和森林資產(chǎn)評估等方面具有重要意義。21世紀(jì)初，馬尾松材種出材率表的研究受到一些學(xué)者的重視。江希鈿等[2]使用改進的Demearschalk削度方程對馬尾松干形進行了研究，隨后又做了馬尾松二元材種出材率表的編制研究[3]，但并未應(yīng)用其之前進行的削度方程研究成果，而僅使用材積比方程。同年，林劍峰[4]和王鵬程等[5]也做了類似研究，前者仍采用材積比方程編表，后者則采用削度-材積系統(tǒng)編表。本研究采用類似削度-材積系統(tǒng)的編表方式，在建立削度方程、累積材積方程和樹皮材積方程基礎(chǔ)上，編制二元材種出材率表。

1 建模及編表資料

共使用406株解析木（下文簡稱 “樣木”）數(shù)據(jù)，其中26 cm徑階以下樣木來自貴州省龍里林場，28～66 cm徑階樣木來自貴州省第3次森林資源規(guī)劃設(shè)計調(diào)查（由貴州省林業(yè)調(diào)查規(guī)劃設(shè)計院提供）。全部樣木來自貴州29個縣區(qū)，分屬貴陽、安順、六盤水、畢節(jié)、遵義、銅仁、黔南和黔東南等8個地市。來源地以黔南和黔東南為主，分別占總株數(shù)的53.7%和31.8%。

來自龍里林場的解析木共163株。每株樣木數(shù)據(jù)包括樹干全高（H）和11個分別位于1.3 m，0～9/10 H高處的斷面的帶去皮徑，以及現(xiàn)場造材數(shù)據(jù)。另外，243株樣木數(shù)據(jù)包括樹干全高H和13個分別位于1.3 m，0～9/10 H，1/4 H和3/4 H高處的斷面的帶去皮徑，無造材數(shù)據(jù)。

全部樣木分布在6～66 cm共31個徑階內(nèi)。其中15株的徑階21個，占徑階總數(shù)的2/3略強，低于10株的徑階6個，11～14株的徑階4個。在盡量保證各徑階建模樣木滿足10株的前提下，將全部樣木分成2組，分別用于模型擬合和模型檢驗。數(shù)據(jù)詳情見表1。

表1 樣木數(shù)據(jù)基本統(tǒng)計特征Table1 Mean sample tree characteristics of Pinus massoniana

2 模型擬合與優(yōu)選

2.1 擬合方法

方程的擬合在R軟件[6]中實現(xiàn)。由功能包nlrwr[7]中的nls非線性回歸函數(shù)進行擬合，并結(jié)合擬合效果診斷，適當(dāng)采用boxcox.nls函數(shù)，對原方程進行Box-Cox函數(shù)轉(zhuǎn)換，以削弱異方差性，且不改變原模型變量之間的關(guān)系。Box-Cox轉(zhuǎn)換函數(shù)[7]表達式如下：

式（1）中： y=f（x，β），β 為模型參數(shù)， λ 為待估參數(shù)。

另外， Mak-Burkhart方程[式（2）]和 Kozak方程[式（3）]含有需在模型回歸之前預(yù)先給定的常數(shù)（a1，a2和p；下文稱 “特征參數(shù)”）。對此特征參數(shù)采用取值空間格網(wǎng)搜索的方法，以擬合殘差平方和最小為準(zhǔn)則確定。

2.2 模型優(yōu)選原則

2.2.1 多指標(biāo)綜合比較 各模型的優(yōu)選在考慮回歸模型常用評價指標(biāo)（剩余平方和、剩余標(biāo)準(zhǔn)差、復(fù)相關(guān)系數(shù)、修正復(fù)相關(guān)系數(shù)和信息量準(zhǔn)則）的同時，還考察了總相對誤差、平均系統(tǒng)誤差、平均相對誤差絕對值和預(yù)估精度等指標(biāo)[8]。

2.2.2 進行適合性檢驗 對削度方程進行適合性檢驗，即利用檢驗數(shù)據(jù)，分別對各方程在指定樹干5個相對位置（0.1H，0.2H，0.5H，0.8H和0.9H）的預(yù)估效果，以偏差（W），絕對偏差（|W|），相對偏差（MSD）和平均相對偏差（ME）為指標(biāo)， 進行綜合評分[9]比較。

2.2.3 模型優(yōu)化 對模型參數(shù)擬合結(jié)果不能全部滿足統(tǒng)計顯著性要求的模型進行參數(shù)優(yōu)化。

2.3 削度方程

2.3.1 試驗方程 選用3個具有代表性的削度方程，分別是Demearschalk削度方程[式（2）]，Mak-Burkhart的分段削度方程[式（3）]和Kozak的可變指數(shù)削度方程[式（4）]。其形式如下：

以上各式中，D為帶皮胸徑，H為樹干全高，h為斷面高度，d為h高度處斷面的去皮直徑，ai（i=0，1，2）和bi（i=0，1，…，5）為待估參數(shù)；徑和高的單位分別為cm和m。

2.3.2 模型擬合 經(jīng)預(yù)估，式（3）特征參數(shù)a1和a2分別取值0.08和0.73，式（4）特征參數(shù)p取值0.05。由于式（4）的b2參數(shù)擬合結(jié)果統(tǒng)計上不顯著（P=0.1317），變動系數(shù)等于66.33%，需對該參數(shù)進行優(yōu)化。經(jīng)去除參數(shù)b2后，得到新模型（4*）：

式（4*）中各符號含義同前文。用該模型取代式（4）做進一步的模型評價與比較。Demearschalk方程[式（2）]原為一致性削度方程[10]，即其積分材積與山本式材積一致。為保留其一致性，采用分步擬合的方法，即先擬合山本式材積方程的3個參數(shù)，再擬合第4個參數(shù)a1。

估計過程中， 分別對式（2）（3）（4）和（4*）進行了 Box-Cox函數(shù)轉(zhuǎn)換， λ 取值均為 0.4。 經(jīng) Box-Cox函數(shù)轉(zhuǎn)換后，異方差得到式一定程度的削弱，殘差的正態(tài)性得到改善。各方程參數(shù)擬合結(jié)果見表2。

表2 削度方程參數(shù)估計結(jié)果Table2 Parameters estimation taper function models

2.3.3 模型優(yōu)選 參數(shù)估計結(jié)果（表2）顯示，式（2）（3）和（4*）各參數(shù)估計值的變動系數(shù)均小于11%，各參數(shù)估計值穩(wěn)定性很高。參數(shù)估計的t檢驗也顯示，全部參數(shù)在0.001水平上顯著。用建模樣木數(shù)據(jù)，計算并分析各統(tǒng)計指標(biāo)（表3和表4）。由表3各統(tǒng)計指標(biāo)計算結(jié)果可知，式（4*）最優(yōu)，式（3）次之。進一步考察表4，式（4*）在整體及全部徑級的平均相對誤差絕對值（RMA）和預(yù)估精度（P）這2項指標(biāo)上，均表現(xiàn)最優(yōu)，在總相對誤差（RS）指標(biāo)上，式（4*）也表現(xiàn)最優(yōu)。經(jīng)綜合分析，可以確認式（4*）最優(yōu)。

另用檢驗數(shù)據(jù)（樣木數(shù)116）對3個方程進行了適合性檢驗，經(jīng)綜合分析比較，仍以式（4*）最優(yōu)。

表3 削度方程模型擬合效果評價指標(biāo)Table3 Estimation indexes for regression of taper function models

表4 削度方程模型分徑級評價指標(biāo)Table4 Estimation indexes calculated by diameter class for taper function models

2.4 累積材積方程

從理論上講，削度方程沿整個樹干區(qū)間的定積分即為樹干材積。但通過其積分而估算的方法為間接方法，其精度需要驗證。為此，以優(yōu)選出的削度式（4*）的積分式為基本型。另外，構(gòu)建了線性修正式（7）和冪函數(shù)修正式（8），加上山本材積式（5），組成以下4個方程，用于比較和檢驗：

以上式（5）～（8）中，V為帶皮材積，V′為去皮材積，Vk（h）為直接通過kozak方程積分而得的0～h段去皮材積，V′（h）為0～h段去皮材積，h為斷面高度，fk（D，H，x）指Kozak方程高處去皮斷面徑的預(yù)估函數(shù)，b0和b1為待估參數(shù)。

擬合過程中所用的材積觀測值用插值的方法，以Stineman函數(shù)[11]為插值函數(shù)，以1 cm為積分區(qū)間長度，對樹干全長或材段區(qū)間積分求得。Eerikainen[12]曾在研究思茅松Pinus kesiya時采用類似方法。

參數(shù)估計 （表5）顯示，各模型參數(shù)估計值的變動系數(shù)均處于小于4%的低水平。各參數(shù)估計的t檢驗也表明各方程參數(shù)均在0.001水平上顯著。

經(jīng)比較（表 6 和表 7）， 式（5）對全干材積的預(yù)測效果最好，式（6）和（7）表現(xiàn)相當(dāng)， 而式（8）最差。 即便如此， 式（6）和（7）由建模樣本和檢驗樣本得出的系統(tǒng)誤差E均小于式（5），前兩者的預(yù)估精度僅比后者少0.02個百分點。就式（6）和（7）而言，對建模樣本全干和材段去皮材積，前者的修正復(fù)相關(guān)系數(shù)略大于后者，而對建模樣本和檢驗樣本預(yù)估的系統(tǒng)誤差，則后者均略低于前者，其余指標(biāo)兩者也基本一致。綜合而言，兩者相關(guān)不大。由于式 （5）不具預(yù)估材段材積的能力。而在材種出材量估算中，材段材積的預(yù)估能力居很重要地位，因此，選用式（7）作為累積材積方程。

表5 去皮材積方程參數(shù)估計Table5 Parameter estimation of under-bark volume models

表6 全干及材段材積方程比較Table6 Comparation of volume models

表7 全干、材段去皮材積方程模型評價指標(biāo)Table7 Model test indexes for under-bark volume functions

2.5 樹皮率方程

樹皮率是指全干樹皮材積占全干帶皮材積的百分比。由于削度方程以及由其發(fā)展的累積材積方程均針對去皮材段，而材種出材率針對全干帶皮材積，因此，需要估算帶皮材積。由于樹皮材積模型的預(yù)測效果不如去皮材積模型，故以修正削度方程積分材積加樹皮材積作為全干帶皮材積，樹皮率及材種出材率以該帶皮材積計算。

經(jīng)試驗，采用表現(xiàn)優(yōu)秀的式（6）作為樹皮材積方程。該方程形式上與山本材積式相同。

式（9）中：Vb為樹皮材積，其他如前文。

樹皮材積方程參數(shù)擬合結(jié)果（Box-Cox函數(shù)變換參數(shù)λ=-0.2）為：b0=-4.645120（變異系數(shù)=2.55%），b1=1.981140（變異系數(shù) =2.81%），b2=0.478215（變異系數(shù) =29.09%）；剩余標(biāo)準(zhǔn)差 =0.075468，修正復(fù)相關(guān)系數(shù)=0.754505，平均系統(tǒng)誤差=-6.47%，預(yù)估精度為93.91%。

3 單株木二元材種出材率表的編制與精度檢驗

造材和編表程序在R軟件平臺上實現(xiàn)。按照先造大材，后造小材的原則進行造材。材種規(guī)格按照大、中、小徑材小頭直徑分別不小于26.0，20.0，6.0 cm，其材長均不小于2.0 m，短小材小頭直徑不小于4.0 cm且不大于12.0 cm，其材長不小于1.0 m且不大于4.8 m的要求進行。理論造材采用的伐樁高度為10.0 cm。馬尾松單株木二元材種出材率表樣表見表8。全干帶皮材積和樹皮率分別采用式（10）和式（11）：

以上式（10）（11）中，V為全干帶皮材積，V′（H）為全干去皮材積，Vb樹皮材積，Pb（%）為樹皮率。

表8 馬尾松單株木二元材種出材率表簡表Table8 Two-dimention merchantable volume yielding volume rate tables of Pinus massoniana

編表精度的檢驗采用式（12），使用精度的檢驗采用置信橢圓F檢驗[9]。各材種的觀測值采用模擬造材的方法獲取。該方法以樣木為單位，以各斷面帶去皮徑為基礎(chǔ)點，用Stine插值法獲取帶去皮樹干的模擬曲線，再按照造材規(guī)格要求進行造材。精度檢驗的結(jié)果見表9。從表中可以看出，大、中、小徑材3個材種均滿足國家標(biāo)準(zhǔn)GB/T 20381-2006“材種出材率表編制技術(shù)規(guī)程”（下文簡稱 “國家標(biāo)準(zhǔn)”）不超過±5%的要求，經(jīng)濟材（包括大、中、小徑材）平均系統(tǒng)誤差為0.58%，滿足國家標(biāo)準(zhǔn)不超過±3%的要求。短小材和薪材的平均系統(tǒng)誤差較大，均超出±5%。置信橢圓F檢驗表明，除小徑材和經(jīng)濟材未通過檢驗（0.01＜P＜0.05）外， 其他各材種均通過檢驗（P＜0.01）。

式（12）中：E為編表精度，yi為具體某一材種觀察值，為該材種估計值。

4 小結(jié)與討論

經(jīng)參數(shù)優(yōu)化的Kozak削度式（4*）， 累積材積式（7）和樹皮材積式（9），共同組成了馬尾松單株木材種出材率模型系統(tǒng)。該系統(tǒng)在R軟件平臺上形成計算機造材程序。對材種出材率表的精度檢驗表明，基于該模型系統(tǒng)編制的馬尾松二元材種出材率表對主要材種（大、中、小徑材和三者共同構(gòu)成的經(jīng)濟材）滿足國家標(biāo)準(zhǔn)要求，且大、中徑材、短小材和薪材通過置信橢圓F檢驗，僅小徑材和由大、中、小徑材組成的經(jīng)濟材未通過此檢驗。因此，可以得出以下2點結(jié)論：①該馬尾松單株木材種出材率模型系統(tǒng)具有一定的可靠性和預(yù)估精度；②該馬尾松二元材種出材率表滿足國家標(biāo)準(zhǔn)技術(shù)要求，具有應(yīng)用和推廣價值。

表9 馬尾松材種出材率表的精度檢驗Table9 Accuracy test for merchantable volume yielding rate tables of Pinus massoniana

就削度方程而言，經(jīng)參數(shù)優(yōu)化的Kozak削度方程式（4*）在所考察的3個方程中表現(xiàn)最優(yōu)，它能很好地體現(xiàn)馬尾松中部飽滿、基部膨大的干形特點，且具有很高的預(yù)估精度（＞99%）。Mak-Burkhart方程表現(xiàn)居中，Demearschalk方程表現(xiàn)表現(xiàn)最差。

累積材積方程的優(yōu)選表明，把非一致性Kozak削度方程引入該方程的做法是可行的。用此方法獲得的累積材積方程，可描述從樹干基部到樹干任意斷面處的材段去皮材積，且在樹干基部和頂部兩端點的取值均有意義。

在樹皮材積模型方面，形如山本材積式的樹皮材積模型在一定程度上提高了樹皮材積的預(yù)估效果。該方法預(yù)估樹皮材積，有別于直接預(yù)估樹皮率的做法，具有一定的借鑒意義。

總之，Kozak削度方程對馬尾松干形的描述表現(xiàn)優(yōu)異，以此方程為基礎(chǔ)構(gòu)建的累積材積方程具有很好的數(shù)學(xué)特性和預(yù)估效果，兩者再加上形如山本材積方程的樹皮材積方程，共同組成具有一定預(yù)估精度和應(yīng)用價值的馬尾松單株木二元材種出材率模型系統(tǒng)。

[1]朱松，夏忠勝.貴州省第3次森林資源規(guī)劃設(shè)計調(diào)查森林資源報告[R]//貴州省林業(yè)廳.貴州省第3次森林資源規(guī)劃設(shè)計調(diào)查報告集.貴陽：貴州大學(xué)，2009：1-80.

[2]江希鈿，楊錦昌，溫素平.馬尾松可變參數(shù)削度方程及應(yīng)用[J].福建林學(xué)院學(xué)報，2000，20（4）：294-297.JIANG Xidian， YANG Jinchang， WEN Suping.Study and application of taper equation of variable parameter of Pinus massoniana [J].J Fujian Coll For， 2000， 20 （4）： 294-297.

[3]江希鈿，林文龍，劉玉明，等.馬尾松人工林材種出材率表的編制[J].林業(yè)勘察設(shè)計，2001（2）：10-13.JIANG Xidian， LIN Wenlong， LIU Yuming， et al.The development of merchantable volume yielding rate tables of masson pine from artifitial forest[J].For Invest Des， 2001 （2）： 10-13.

[4]林劍峰.馬尾松人工林材種出材率表的研究[J].北京林業(yè)大學(xué)學(xué)報，2001，23（4）：35-38.LIN Jianfeng.Research on log rule of timbers from artificial masson pine forest [J].J Beijing For Univ， 2001， 23（4）： 35-38.

[5]王鵬程，莊爾奇，涂炳坤，等.湖北省馬尾松人工林削度方程及材種出材率表的研究[J].華中農(nóng)業(yè)大學(xué)學(xué)報，2001，20（1）： 67-72.WANG Pengcheng， ZHUANG Erqi， TU Bingkun， et al.A study on the taper function and merchantable volume yielding rate table of mason pine in Hubei Province [J].J Huazhong Agric Univ， 2001， 20 （1）： 67-72.

[6]R DEVELOPMENT CORE TEAM.R： A Language and Environment for Statistical Computing [CP/OL].[2010-12-16].http：//ftp.ctex.org/mirrors/CRAN/.

[7]RITZ C， STREIBIG J C.Nonlinear Regression with R [M].New York： Springer， 2008.

[8]曾偉生，駱期邦，賀東北.論加權(quán)回歸與建模[J].林業(yè)科學(xué)，1999，35（5）：5-11.ZENG Weisheng， LUO Qibang， HE Dongbei.Research on weighting regression and modelling [J].Sci Silv Sin， 1999， 35（5）： 5-11.

[9]中華人民共和國國家質(zhì)量監(jiān)督檢驗檢疫總局，中國國家標(biāo)準(zhǔn)化管理委員會.GB/T 20381-2006材種出材率表編制技術(shù)規(guī)程[S].北京：中國標(biāo)準(zhǔn)出版社，2006.

[10]孟憲宇.削度方程和林分直徑結(jié)構(gòu)在編制材種表中的重要意義[J].北京林業(yè)大學(xué)學(xué)報，1991，13（2）：14-19.MENG Xianyu.The significance of taper equations and stand diameter structures in developing merchantable volume tables [J].J Beijing For Univ， 1991， 13 （2）： 14-19.

[11]TOMAS J， HALLDOR B， ICELANDIC M O， et al.Stinepack： Stineman， A Consistently Well Behaved Method of Interpolation[CP/OL].[2009-03-20].http： //ftp.ctex.org/mirrors/CRAN/.

[12]EERIK?INEN K.Stem volume models with random coefficients for Pinus kesiya in Tanzania， Zambia， and Zimbabwe [J].Can J For Res， 2001， 31 （5）： 879-888.