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

沈錢勇 湯孟平

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

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

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

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

1 材料與方法

1.1 數據來源

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

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

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

(1)

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

樣竹按胸徑(D)、竹高(H)、胸高節長(L)和竹稈材積(V)等主要指標統計情況如表1所示。可見，竹稈材積的變異系數最大，胸高節長的變異系數最小。

圖2 樣竹徑階分布

表1 毛竹樣竹實測數據統計

1.2 毛竹竹稈材積模型研建

1.2.1 竹稈材積模型選擇 立木材積方程一般包括4類，即組合變量方程、對數材積方程、Horner材積方程和廣義材積方程。通過初步篩選和比較，本研究基于廣義材積方程，選擇冪函數形式的一元模型及山本式和寺崎渡方程的4個二元模型進行分析比較。一元模型基于以胸徑為自變量的式(2)，二元模型基于以胸徑、竹高和胸高節長為自變量的式(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為胸高節長(cm)；ai為模型參數，i=0、1、2。

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

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

(7)

(8)

(9)

(10)

(11)

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

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

(12)

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

對數轉換線性回歸模型的赤池信息量準則稱為AICcln，非線性回歸模型的赤池信息量準則稱為AICcnorm。

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

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

調整確定系數

(13)

(14)

(15)

(16)

2 結果與分析

2.1 竹稈材積模型誤差結構分析

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

表2 浙江省毛竹竹稈材積模型誤差結構似然分析統計

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

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

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

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

表3 竹稈材積對數回歸模型擬合參數與統計指標①

①括號內數值為各參數t檢驗值。Values in bracket arettest value for each parameter.

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

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

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

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

2.3 對數模型校正

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

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

表4 對數模型反對數轉換在原始尺度下的非線性模型校正評價檢驗①

3 討論

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

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

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

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

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

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

4 結論

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

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

2) 似然函數法是進行模型誤差結構分析比較與模型擬合方式選擇的較好方法。

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

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