浙江省毛竹竹稈材積模型*

沈錢勇 湯孟平

(浙江農(nóng)林大學環(huán)境與資源學院 省部共建亞熱帶森林培育國家重點實驗室 杭州 311300)

毛竹(Phyllostachysedulis)是我國南方重要的森林資源之一，在增加經(jīng)濟效益、發(fā)揮生態(tài)服務(wù)功能等方面具有積極作用(周國模等， 2004； 楊麒麟等， 2017)。浙江省是我國毛竹主產(chǎn)區(qū)，根據(jù)2016年浙江省森林資源與生態(tài)狀況年度監(jiān)測，全省森林面積605.91萬hm2，森林蓄積31 529.17萬m3，其中毛竹林面積80.71萬hm2，占森林面積的13.32%，毛竹林每公頃立竹量3 413株。森林蓄積是當前評價森林生產(chǎn)力和林木生長狀況的重要指標(唐守正等， 2000； 孟憲宇， 2006)，對毛竹林蓄積僅采用每公頃立竹量表示是不精確的，在我國木質(zhì)資源短缺、倡導(dǎo)以竹代木的發(fā)展趨勢下，通過單株毛竹材積測定從而準確估算毛竹材積具有重要意義。

目前，國內(nèi)外對立木材積的估算方法主要有2種(劉鏡婷等， 2016)： 一是采用線性或非線性回歸擬合一元或二元材積模型(Snorrasonetal.， 2006； Caseetal.， 2008)； 二是利用二元材積方程和樹高曲線方程構(gòu)建含誤差變量的相容性聯(lián)立方程組模型進行參數(shù)估計(胥輝，1999； Tangetal.， 2001； 2002； 曾鳴等， 2013)。毛竹竹稈中空、尖削度大，其材積計算不能簡單地套用立木材積方程(陳雙林等， 2008)。20世紀50年代，虞岳世等(1958)研究毛竹材積與胸徑、竹高的關(guān)系，采用竹稈分段測量壁厚方法測定竹壁材積，建立了一元(胸徑)和二元(胸徑、竹高)材積模型； 此后，吳丁(1987)、張振瀛等(1990)也采用類似方法分別建立了一元和二元材積模型。然而，以往建立的材積模型都忽視了竹隔材積，為此，張剛?cè)A等(2007)改進測量方法，除分段測量直徑、壁厚外還測量每個竹隔厚度，并建立了包含竹隔在內(nèi)的一元竹稈材積模型。

應(yīng)當指出，現(xiàn)有毛竹竹稈材積模型均未充分考慮毛竹的竹節(jié)、竹壁和竹隔的形狀不規(guī)則性，不能準確測量竹稈材積，建立的竹稈材積模型不夠精確。而且，毛竹竹桿彎曲生長特性以及冠層枝葉遮擋對竹高的準確測量也有較大影響，造成基于竹高所建立的材積模型存在應(yīng)用的局限性。鑒于此，本研究以浙江省毛竹林為研究對象，采用可測量不規(guī)則物體體積的排水法準確測量竹稈材積，引入胸徑、竹高和胸高節(jié)長因子，通過初步篩選選擇5個材積模型進行比較研究，并通過模型誤差結(jié)構(gòu)及模型擬合優(yōu)度和預(yù)估精度的評價分析，確定適用于浙江省的毛竹竹稈材積模型，以期為毛竹林經(jīng)營管理提供依據(jù)。

1 材料與方法

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

2016年7月至2018年7月，在浙江臨安、慶元、武義、常山、寧海、安吉、諸暨、余姚、黃巖和泰順10個縣(區(qū)、市)，選擇近年未受人為經(jīng)營管理的毛竹林，隨機選取梢頭完整、竹稈通直、斷面近似圓形、無破損和病蟲害的樣竹216株，測量胸徑、竹高和胸高竹節(jié)長(圖1)。樣竹株數(shù)徑階分布如圖2所示，經(jīng)檢驗，樣竹株數(shù)按徑階呈正態(tài)分布(P<0.01)。

圖1 毛竹竹稈材積建模變量示意

伐倒樣竹，采用排水法測定竹稈材積。首先，將水注入定制水桶內(nèi)至水龍頭齊平處，排出桶內(nèi)多余的水； 然后，將竹稈分成竹條放入水桶內(nèi)，用電子提秤稱量通過水龍頭排出的水量； 最后，計算毛竹竹稈材積。公式如下：

(1)

式中：Vs為竹稈材積(dm3)；Mw為排水質(zhì)量(kg)；ρ為水的密度(1 g·cm-3)。

樣竹按胸徑(D)、竹高(H)、胸高節(jié)長(L)和竹稈材積(V)等主要指標統(tǒng)計情況如表1所示?？梢姡穸挷姆e的變異系數(shù)最大，胸高節(jié)長的變異系數(shù)最小。

圖2 樣竹徑階分布

表1 毛竹樣竹實測數(shù)據(jù)統(tǒng)計

1.2 毛竹竹稈材積模型研建

1.2.1 竹稈材積模型選擇 立木材積方程一般包括4類，即組合變量方程、對數(shù)材積方程、Horner材積方程和廣義材積方程。通過初步篩選和比較，本研究基于廣義材積方程，選擇冪函數(shù)形式的一元模型及山本式和寺崎渡方程的4個二元模型進行分析比較。一元模型基于以胸徑為自變量的式(2)，二元模型基于以胸徑、竹高和胸高節(jié)長為自變量的式(3)～(6)(孟憲宇， 2006)，5個毛竹竹稈材積方程為：

M1：V=a0Da1；

(2)

M2：V=a0Da1Ha2；

(3)

M3：V=a0Da1La2；

(4)

M4：V=a0Da1ea2/H；

(5)

M5：V=a0Da1ea2/L。

(6)

式中：V為竹稈材積(dm3)；D為胸徑(cm)；H為竹高(m)；L為胸高節(jié)長(cm)；ai為模型參數(shù)，i=0、1、2。

基于以上5個材積方程建立模型需考慮誤差項，有加性和乘性之分。對于含乘性誤差項的模型，利用對數(shù)回歸進行參數(shù)估計后需要轉(zhuǎn)化成原始結(jié)構(gòu)，且需對誤差項進行修正(Baskerville, 1974； 曾偉生等， 2011a)。

在對數(shù)回歸模型研究中,學者們提出過許多校正因子(Finney，1941； Baskerville，1974； Snowdon，1991)，其中以Baskerville(1974)提出的校正因子[CF=exp(s2/2)](s為回歸估計標準差)應(yīng)用最多。曾偉生等(2011b)則基于對數(shù)轉(zhuǎn)換特點提出了新的校正因子[CF=(1+s2/2)]。通過比較發(fā)現(xiàn)，采用不同校正因子修正后的模型各項評價指標并無差異(Dongetal.， 2014)，本研究選擇的校正因子為CF=(1+s2/2)。預(yù)估材積模型如下：

(7)

(8)

(9)

(10)

(11)

1.2.2 毛竹竹稈材積模型誤差結(jié)構(gòu)分析 采用似然函數(shù)法，分析基于5個不同方程的毛竹竹稈材積模型誤差結(jié)構(gòu)(Xiaoetal.， 2011； Dongetal.， 2014； 董利虎等， 2016)，確定應(yīng)當采用對數(shù)回歸或非線性回歸模型進行模型擬合。具體計算步驟如下：

第1步： 利用原始數(shù)據(jù)分別進行對數(shù)線性回歸和非線性回歸擬合，得到5個模型的參數(shù)估計值(ai)和方差(σ2)。計算5個模型各自的對數(shù)似然值(lnL)，通過下式計算各模型的赤池信息量準則(AICc)：

(12)

式中：k為模型參數(shù)個數(shù)；n為建模樣本數(shù)量； lnL為模型的對數(shù)似然值。

對數(shù)轉(zhuǎn)換線性回歸模型的赤池信息量準則稱為AICcln，非線性回歸模型的赤池信息量準則稱為AICcnorm。

第2步： 分別對比5個模型的2個赤池信息量準則，若ΔAICc(AICcnorm-AICcln)<-2，表明模型誤差項是相加的，應(yīng)基于原始數(shù)據(jù)進行非線性回歸擬合； 如果ΔAICc> 2，則說明模型誤差項是相乘的，應(yīng)進行對數(shù)轉(zhuǎn)換的線性回歸擬合。

1.2.3 模型評價與檢驗 模型擬合優(yōu)度和預(yù)估精度通過以下4個指標進行評價檢驗：

調(diào)整確定系數(shù)

(13)

(14)

(15)

(16)

2 結(jié)果與分析

2.1 竹稈材積模型誤差結(jié)構(gòu)分析

分別利用原始數(shù)據(jù)的非線性回歸和對數(shù)轉(zhuǎn)換的線性回歸擬合5個基礎(chǔ)材積模型[式(2)～(6)]，求得各模型的赤池信息量準則(AICcnorm和AICcln)，相比較得到5個模型對應(yīng)的ΔAICc(表2)。

表2 浙江省毛竹竹稈材積模型誤差結(jié)構(gòu)似然分析統(tǒng)計

由表2可知，通過對數(shù)轉(zhuǎn)換線性回歸擬合的模型與利用原始數(shù)據(jù)直接擬合的非線性模型獲得的AICc相比，5個模型的ΔAICc均大于2，說明本研究浙江省毛竹竹稈材積模型誤差結(jié)構(gòu)為乘積型誤差，5個模型都應(yīng)當采用對數(shù)轉(zhuǎn)換的線性回歸進行擬合分析。

2.2 竹稈材積模型擬合與評價檢驗

不分建模樣本和檢驗樣本而利用全部樣本建立材積模型，能充分利用其信息，可得到預(yù)估誤差最小的模型(Myers，1986； Kozaketal.， 2003； 曾偉生等， 2011c； 唐思嘉， 2017)。本研究利用全部樣本，采用對數(shù)轉(zhuǎn)換的線性回歸進行竹稈材積模型擬合分析，并通過擬合優(yōu)度和預(yù)估精度指標對模型進行評價和檢驗。

比較各模型的統(tǒng)計指標結(jié)果，模型M2最優(yōu)，其次是M5。從便于實際應(yīng)用角度，確定采用包含胸徑和胸高節(jié)長的毛竹竹稈材積模型M5。

表3 竹稈材積對數(shù)回歸模型擬合參數(shù)與統(tǒng)計指標①

①括號內(nèi)數(shù)值為各參數(shù)t檢驗值。Values in bracket arettest value for each parameter.

從對數(shù)回歸模型殘差分布(圖3)也可看出，各模型殘差隨著預(yù)測值增大基本呈均勻分布，不存在異方差。

以上檢驗是基于全部樣竹的平均值，并不能很好反映出不同大小(胸徑)毛竹竹稈材積的預(yù)估情況(Kozaketal.， 2003)。因此，本研究對5個對數(shù)回歸模型各徑階的平均偏差進行檢驗(圖4)?？梢姡?個模型的預(yù)估偏差隨徑階并未發(fā)生明顯的變化，均接近于0。從數(shù)值上看，各模型在4.0～5.9 cm和14.0～15.9 cm徑階時，偏差相比較其他徑階大，其中前者略有低估，后者略有高估。由圖4也可看出，各模型在不同徑階范圍的預(yù)估精度均較高。

圖3 竹稈材積對數(shù)回歸模型殘差分布

圖4 不同徑階竹稈材積對數(shù)模型平均偏差檢驗

2.3 對數(shù)模型校正

基于對數(shù)轉(zhuǎn)換的線性回歸模型預(yù)測的是期望材積對數(shù)值，而要獲得材積實際值(即算術(shù)尺度下的材積預(yù)估值)，則需要對對數(shù)模型預(yù)測值進行反對數(shù)轉(zhuǎn)換(Dongetal.， 2014)。一般認為，反對數(shù)轉(zhuǎn)換過程中會對材積等預(yù)測值產(chǎn)生系統(tǒng)上的低估，因此需要對其進行校正(Finney，1941； Baskerville，1974； Packardetal.， 2008； Cliffordetal.， 2013)。

表4所示為5個模型反對數(shù)轉(zhuǎn)換后的部分統(tǒng)計指標。可以看出，對對數(shù)模型進行校正后，模型擬合優(yōu)度和預(yù)估精度變化不顯著。Madgwick等(1975)研究指出，反對數(shù)轉(zhuǎn)換引入校正因子會高估預(yù)測值，認為當對數(shù)模型誤差項足夠小時，反對數(shù)轉(zhuǎn)換過程中無需進行校正； Dong等(2014)比較不同校正因子進行反對數(shù)轉(zhuǎn)換模型校正與不做校正的模型統(tǒng)計指標，也發(fā)現(xiàn)同樣規(guī)律。根據(jù)表4，利用本研究得到的對數(shù)模型進行反對數(shù)轉(zhuǎn)換得到材積預(yù)估值，從ME和MSE可以看出，模型本身具有較高的擬合優(yōu)度和預(yù)估性能，誤差較?。?而通過校正后，則會高估預(yù)測值，平均偏差和平均系統(tǒng)誤差均明顯增大。

表4 對數(shù)模型反對數(shù)轉(zhuǎn)換在原始尺度下的非線性模型校正評價檢驗①

3 討論

常用的生物模型如材積模型，一般采用冪函數(shù)等非線性函數(shù)的對數(shù)回歸擬合或直接非線性擬合的方法獲得相應(yīng)模型(Parresol， 2001； Zianisetal.， 2011； Dongetal.， 2014)，而對于采用何種方法進行模型擬合能獲得最佳效果依然存在爭論。近年來，林業(yè)和生態(tài)研究領(lǐng)域的一些學者也開始利用似然函數(shù)法對林木生物量等模型的誤差結(jié)構(gòu)進行分析，以確定應(yīng)采用何種建模方式進行建模(Xiaoetal.， 2011； Ballantyne， 2013； Laietal.， 2013； Dongetal.， 2014； 董利虎等， 2016)。

Dong等(2014)為了說明不同誤差結(jié)構(gòu)對模型擬合精度的影響，對落葉松根部生物量數(shù)據(jù)同時進行了基于冪函數(shù)形式的對數(shù)回歸模型和非線性模型的研建，為保證比較模型具有相同的響應(yīng)變量，其針對2種形式模型分別進行了原始尺度下的反對數(shù)轉(zhuǎn)換和對數(shù)轉(zhuǎn)換尺度下的對數(shù)轉(zhuǎn)換，即LR模型由線性回歸直接擬合，NLR*模型的相關(guān)統(tǒng)計指標基于直接非線性擬合的NLR模型預(yù)測值進行對數(shù)轉(zhuǎn)換重新計算獲取，利用對數(shù)尺度下的數(shù)據(jù)進行預(yù)測評價檢驗； 而在原始尺度下，NLR模型直接利用原始數(shù)據(jù)非線性擬合，LR*模型的相關(guān)統(tǒng)計指標則基于LR模型預(yù)測值進行反對數(shù)轉(zhuǎn)換重新計算獲取，利用的是算術(shù)尺度數(shù)據(jù)。本研究比較了2種尺度下模型M5各徑階的平均百分標準偏差(MPSE，%)，計算公式見參考文獻(曾偉生等， 2011)，結(jié)果如圖5所示。可以看出，在原始尺度下，NLR模型和LR*模型各徑階預(yù)估精度變化較小，分別在4.0～5.9 cm徑階和12.0～13.9 cm徑階時最高。在徑階4.0～9.9 cm之間，LR*模型的MPSE均小于NLR模型，但徑階大于10 cm時，2種模型精度值接近(圖5A)； 而在對數(shù)轉(zhuǎn)換尺度下，2種模型預(yù)估精度均呈隨徑階增大而升高的趨勢，且LR模型的MPSE在各徑階基本低于或接近NLR*模型(圖5B)。不難看出，當模型誤差結(jié)構(gòu)為乘積型時(對數(shù)尺度下)，利用對數(shù)轉(zhuǎn)換線性模型具有更佳的預(yù)估效果，偏差較小。

圖5 原始尺度(A)和對數(shù)轉(zhuǎn)換尺度(B)下線性模型和非線性模型各徑階的平均百分標準偏差

通過不同徑階竹稈材積對數(shù)模型平均偏差檢驗結(jié)果(圖4)與不同尺度下線性模型和非線性模型預(yù)估精度(圖5)的比較可以發(fā)現(xiàn)，各模型對于不同徑階毛竹竹稈材積的預(yù)估效果均較好，差異不大。但應(yīng)看到，在最小徑階時，模型預(yù)估精度易出現(xiàn)波動(圖5B)，這或許與該徑階樣本數(shù)據(jù)較少有關(guān)(4.0～5.9 cm徑階僅有12株)，后續(xù)應(yīng)對大徑階和小徑階毛竹開展更多研究。

許多關(guān)于對數(shù)回歸模型的研究均提出需對模型進行校正以減少或消除偏差的影響。本研究對比5個對數(shù)模型進行反對數(shù)轉(zhuǎn)換校正和未經(jīng)校正的統(tǒng)計指標(表4)發(fā)現(xiàn)，對數(shù)模型擬合優(yōu)度較高、預(yù)估誤差較小，進行反對數(shù)轉(zhuǎn)換校正后，會高估材積值，與Madgwick等(1975)的研究結(jié)果一致。故本研究在進行對數(shù)模型反對數(shù)轉(zhuǎn)換時未進行校正。

與二元材積模型相比，一元模型的擬合優(yōu)度和各徑階偏差等并無顯著差異(表3、圖4)，亦能滿足單株材積或單位面積竹林蓄積預(yù)估的精度要求，在生產(chǎn)上具有簡便和實用性； 但引入胸高節(jié)長和竹高變量后，所建立的二元竹稈材積模型具有更高的擬合優(yōu)度。在實踐中，由于胸高節(jié)長較竹高更易準確測量，因此包含胸高節(jié)長變量的毛竹材積或生物量模型具有進一步研究的應(yīng)用價值； 同時，有必要在更大范圍進行數(shù)據(jù)調(diào)查，提高模型精度，增強模型的適用性，并推廣應(yīng)用于我國毛竹林經(jīng)營實踐。

4 結(jié)論

本研究以毛竹為研究對象，采用排水法準確測定毛竹樣竹竹稈材積，以胸徑(D)、竹高(H)和胸高節(jié)長(L)為自變量，建立5個毛竹竹稈材積模型，通過模型誤差結(jié)構(gòu)和模型評價檢驗比較分析，得到以下主要結(jié)論：

1) 排水法是準確測量毛竹竹稈材積的有效手段，可以提高材積測量精度。

2) 似然函數(shù)法是進行模型誤差結(jié)構(gòu)分析比較與模型擬合方式選擇的較好方法。

3) 引入竹高(H)和胸高節(jié)長(L)變量后，基于胸徑-竹高(D-H)和胸徑-胸高節(jié)長(D-L)的二元竹稈材積模型評價檢驗指標有所提高，二元材積模型優(yōu)于一元模型。

4) 模型M5與M2相比，擬合優(yōu)度和預(yù)估精度略有提高，考慮實踐中變量測量的方便和準確性，基于胸徑-胸高節(jié)長的模型M5為預(yù)估毛竹竹稈材積的最優(yōu)模型，即V=0.191 2D2.114 9e-6.841 1/L。