沈錢勇 湯孟平

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

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

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

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

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

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

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

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

(1)

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

樣竹按胸徑(D)、竹高(H)、胸高節(jié)長(zhǎng)(L)和竹稈材積(V)等主要指標(biāo)統(tǒng)計(jì)情況如表1所示?？梢?jiàn)，竹稈材積的變異系數(shù)最大，胸高節(jié)長(zhǎng)的變異系數(shù)最小。

圖2 樣竹徑階分布

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

1.2 毛竹竹稈材積模型研建

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

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é)長(zhǎng)(cm)；ai為模型參數(shù)，i=0、1、2。

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

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

(7)

(8)

(9)

(10)

(11)

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

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

(12)

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

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

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

1.2.3 模型評(píng)價(jià)與檢驗(yàn) 模型擬合優(yōu)度和預(yù)估精度通過(guò)以下4個(gè)指標(biāo)進(jìn)行評(píng)價(jià)檢驗(yàn)：

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

(13)

(14)

(15)

(16)

2 結(jié)果與分析

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

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

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

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

2.2 竹稈材積模型擬合與評(píng)價(jià)檢驗(yàn)

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

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

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

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

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

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

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

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

2.3 對(duì)數(shù)模型校正

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

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

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

3 討論

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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