貴州馬尾松單株木二元材種出材率表的編制

金 明，丁貴杰

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

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

1 建模及編表資料

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

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

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

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

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

2.1 擬合方法

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

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

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

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

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

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

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

2.3 削度方程

2.3.1 試驗(yàn)方程 選用3個(gè)具有代表性的削度方程，分別是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)計(jì)上不顯著（P=0.1317），變動(dòng)系數(shù)等于66.33%，需對(duì)該參數(shù)進(jìn)行優(yōu)化。經(jīng)去除參數(shù)b2后，得到新模型（4*）：

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

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

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

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

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

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

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

2.4 累積材積方程

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

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

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

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

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

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

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

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

2.5 樹皮率方程

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

經(jīng)試驗(yàn)，采用表現(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 單株木二元材種出材率表的編制與精度檢驗(yàn)

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

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

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

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

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

4 小結(jié)與討論

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

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

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

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

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

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

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

[2]江希鈿，楊錦昌，溫素平.馬尾松可變參數(shù)削度方程及應(yīng)用[J].福建林學(xué)院學(xué)報(bào)，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è)計(jì)，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é)報(bào)，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é)報(bào)，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]中華人民共和國(guó)國(guó)家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局，中國(guó)國(guó)家標(biāo)準(zhǔn)化管理委員會(huì).GB/T 20381-2006材種出材率表編制技術(shù)規(guī)程[S].北京：中國(guó)標(biāo)準(zhǔn)出版社，2006.

[10]孟憲宇.削度方程和林分直徑結(jié)構(gòu)在編制材種表中的重要意義[J].北京林業(yè)大學(xué)學(xué)報(bào)，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.