粟 軍
(西雙版納布龍州級自然保護區(qū)管護所,云南 景洪 666100)
利用一元線性回歸方程計算盜伐林木材積試驗
粟 軍
(西雙版納布龍州級自然保護區(qū)管護所,云南 景洪 666100)
盜伐林木案件中,如現(xiàn)地伐樁保存完好,但無采伐木,周圍有林相、結(jié)構(gòu)相同或相近林分的情況下,利用樣木的根徑和胸徑測量數(shù)據(jù),建立線性回歸方程計算出盜伐林木的回歸胸徑,再利用其樹高的測量數(shù)據(jù)建立樹高曲線對數(shù)回歸方程,計算出盜伐林木的回歸樹高,并利用統(tǒng)計學原理,對回歸胸徑進行相關系數(shù)檢驗,對回歸樹高進行F檢驗。試驗結(jié)果認為,利用該方法能使盜伐林木的蓄積量計算值更加接近實際值。
盜伐林木;胸徑回歸;樹高回歸;一元線性回歸;二元材積;相關系數(shù)檢驗;F檢驗;材積計算
在《云南省林地鑒定規(guī)范(暫行)》第十七條規(guī)定中,關于現(xiàn)地伐樁保存完好,但無采伐木,周圍有林相、結(jié)構(gòu)相同相近林分的,直接每木檢尺伐樁、鑒定樹種、查數(shù)年輪,統(tǒng)計株數(shù)。在有代表性的地段設標準地,建立胸徑與根徑回歸公式、胸徑與樹高的回歸公式或樹高曲線圖,用二元立木公式計算蓄積。根據(jù)上述規(guī)定,在盜伐林木案件中可以借鑒上述規(guī)定和方法計算林木蓄積。
通過實驗,驗證一元線性回歸分析方法計算現(xiàn)地伐樁材積的可行性,了解一元線性回歸分析方法和原理在計算現(xiàn)地伐樁材積中地應用,使理論與實踐相結(jié)合,并把此方法更加緊密地應用于計算現(xiàn)地伐樁材積中。
外業(yè)工具有測樹圍尺、測高器、皮尺、GPS儀、地形圖等。用測樹圍尺測量16棵被盜伐銀葉錐根徑(D5cm),在盜伐林木現(xiàn)場周圍對42棵坡向、林相、結(jié)構(gòu)、樹種相同的樣木的根徑(D5cm)、胸徑(D130cm)用測樹圍尺測量,用測高器測量樣木的樹高。盜伐林木根徑測量數(shù)據(jù)詳見表1,樣木測量數(shù)據(jù)詳見表2。內(nèi)業(yè)整理需電腦一臺,把外業(yè)測量的盜伐林木伐樁數(shù)據(jù)、樣木的根徑、胸徑、樹高數(shù)據(jù)錄入Excel文件中。
2.1.1 建立回歸模型
利用樣木的根徑和胸徑的測量數(shù)據(jù),建立線性回歸。預測對象“胸徑”作為因變量y, 將主要影響因素“地徑”作為自變量x,它們之間的關系用一元線性回歸模型表示為:y=a +bx。
參數(shù)計算公式為:
快速計算方法:應用Excel插入散點圖,添加線性趨勢線,Excel可自動計算出a、b參數(shù),用線性回歸方程推算被采伐木的胸徑(圖1)。 經(jīng)計算,a=-0.382 1,b=0.926 9。
2.1.2 相關系數(shù)檢驗
相關系數(shù)是描述兩個變量之間的線性相關關系的密切程度的數(shù)量指標,用R表示。當R=1時,變量x、y完全正相關;R在-1~1之間,當0﹤R﹤1時,變量x和y正相關,當-1﹤R﹤0時,為負相關;當R=0時,變量x和y沒有線性關系。通過應用Excel插入散點圖,添加線性趨勢線計算出R值后,查相關系數(shù)檢驗表,在自由度N-2(N為樣本個數(shù))和顯著性水平a(一般取a=0.05)下,若R大于臨界值,則變量x和y之間線性關系成立,否則兩個變量不存在線性關系。
在自由度(f=n-2, 即f=42-2 )和顯著性水平a=0.05的情況下,查“相關系數(shù)臨界值表”得值0.304 4。通過計算,相關系數(shù)R=0.963 ,大于臨界值0.304 4,且R值接近1 ,說明不但該地林木根徑x與胸徑y(tǒng)之間的線性關系成立,而且線性關系顯著,說明以該地被采伐林木根徑推算被采伐林木胸徑準確度高。利用回歸方程y=0.926 9x-0.382 1計算采伐木回歸胸徑結(jié)果詳見表2。
表1 樣木根徑、胸徑、樹高
表2 胸徑回歸統(tǒng)計
1)利用調(diào)查樣木的胸徑和樹高的測量數(shù)據(jù)作為樹高曲線的基礎數(shù)據(jù)(表1)。
2)剔除樣木異常數(shù)據(jù)。把測量數(shù)據(jù)按升序排列,然后利用這些數(shù)據(jù)繪成胸徑樹高曲線散點圖(圖2),利用這個散點圖用肉眼剔除不可靠數(shù)據(jù)(表3)。
3)根據(jù)樣木胸徑樹高的關系生成胸徑樹高曲線散點圖,利用胸徑樹高對數(shù)曲線回歸公式y(tǒng)=aln(x)+b計算回歸分析后的樹高,式中y為回歸樹高、a、b為參數(shù)、ln為自然對數(shù)、x為回歸胸徑。
a、b參數(shù)計算公式:把樣木胸徑轉(zhuǎn)換為自然對數(shù)值,即x=ln(x), 然后利用上述一元線性回歸模型y=a+bx及參數(shù)公式計算出對數(shù)回歸曲線公式中的a、b參數(shù)。快速計算方法:在Excel中插入散點圖,添加對數(shù)趨勢線自動計算出a、b參數(shù)a=7.501 5、b=-11.85 。
4)在Excel中計算出整化后的回歸胸徑自然對數(shù)值ln(x),把對應的自然對數(shù)代入y=7.501 5ln(x)-11.85中,計算出回歸后的樹高(表4)。
5)對回歸樹高進行F檢驗
F檢驗計算公式為:
式中:F為統(tǒng)計量,s2為方差。
在樣木中隨機抽取16棵樣木樹高與回歸樹高進行F檢驗,在Excel 中進行快速F檢驗(表5)。
圖1 根徑胸徑一元線性回歸散點圖Fig.1 Scatter diagram of unary linear regression between ground diameter and DBH
圖2 胸徑樹高對數(shù)曲線Fig.2 Logarithmic curve between DBH and tree height
序號胸徑/cm樹高/m157.514260.015352.721448.614556.415
表4 樹高回歸計算統(tǒng)計
表5 F-檢驗 雙樣本方差分析
在表6中可以看到,F(xiàn)統(tǒng)計量為1.932 095 82,F(xiàn)單尾臨界為2.403 447 07,F(xiàn)統(tǒng)計量
假設在盜伐林木現(xiàn)場只盜伐了1棵林木,是否也可以用上述方法計算立木材積呢?假設在盜伐現(xiàn)場有一棵銀葉錐伐倒木,根徑(D5cm)為56.5 cm,胸徑(D130cm)53.7 cm,樹高為17 m。對伐樁利用硬闊一元材積公式計算得材積0.822 m3,胸徑樹高利用南亞熱帶闊葉樹二元立木材積公式計算得材積1.661 m3。假設只有伐樁在現(xiàn)場,樹干已運出林區(qū),根據(jù)小樣本容量n≥k+1(k為解釋變量的數(shù)目)。利用在盜伐林木現(xiàn)場范圍選取6棵坡向、林相、結(jié)構(gòu)、樹種相同和根徑相似的樣木進行以該地被采伐林木根徑推算被采伐林木胸徑準確度高。利用樹高曲線對數(shù)回歸公式計算得樹高為18 m,利用回歸后的胸徑樹高用南亞熱帶闊葉樹二元立木材積公式計算得材積為1.711 m3。假設材積1.661 m3為實際值,利用一元材積計算得材積為0.822 m3,誤差值為51.5%,利用回歸胸徑樹高計算材積得1.711 m3,誤差值為3.01%,說明利用回歸后的胸徑樹高更接近實際值(表6)。
表6 小樣本求材積統(tǒng)計
在云南省西雙版納州景洪市勐龍鎮(zhèn)陸拉村委會松林新寨村小組集體林地內(nèi)(現(xiàn)場中心坐標為E 100° 32′ 1.6″, N 21° 35′ 28.8″),對盜伐的16棵銀葉錐伐樁和留在現(xiàn)場的伐倒木樹干利用測樹圍尺和皮尺對根徑、胸徑、樹高進行測量,并對周圍林相、結(jié)構(gòu)相同林分的42棵銀葉錐樣木利用測樹圍尺和測高器對樣木根徑、胸徑、樹高進行測量。對16棵盜伐的銀葉錐分別利用硬闊葉樹種根徑一元材積公式、南亞熱帶闊葉樹二元立木材積公式,以及利用回歸后的胸徑、樹高用南亞熱帶闊葉樹二元立木材積公式計算的蓄積進行對比(表7)。
經(jīng)計算,對16棵盜伐林木根徑數(shù)據(jù)利用硬闊葉樹種根徑一元材積公式計算得13.623 m3, 對胸徑樹高數(shù)據(jù)利用南亞熱帶闊葉樹二元立木材積公式計算得28.024 m3, 對回歸后的胸徑、樹高用南亞熱帶闊葉樹二元立木材積公式計算材積得27.395 m3。假設16棵盜伐林木以28.024 m3材積為實際值,那么,利用硬闊葉樹種根徑一元材積公式計算得的13.623 m3與實際值誤差為53.4%,利用回歸后的胸徑、樹高用南亞熱帶闊葉樹二元立木材積公式計算得的27.395 m3與實際值的誤差為2.3%。表7對比說明,對盜伐林木伐樁保存完好,但無采伐木,周圍有林相、結(jié)構(gòu)相同或相近林分的情況下,對回歸后的胸徑、樹高用二元立木材積公式計算得到的材積比根徑一元材積公式計算的材積更接近于實際值。
盜伐林木案件中,現(xiàn)地伐樁保存完好,但無采伐木,周圍有林相、結(jié)構(gòu)相同相近林分的情況下,利用樣木的根徑和胸徑的測量數(shù)據(jù),在Excel中建立根徑胸徑線性回歸散點圖及線性回歸公式計算出盜伐林木的回歸胸徑,利用樣木的胸徑和樹高的測量數(shù)據(jù),在Excel中建立樹高曲線散點圖及對數(shù)回歸公式,計算出盜伐林木的回歸樹高,并利用統(tǒng)計學原理,對回歸胸徑進行相關系數(shù)檢驗,對回歸樹高進行F檢驗。通過試驗,實現(xiàn)盜伐林木與樣木擬合后利用二元立木公式計算蓄積,使盜伐林木的蓄積量通過線性回歸更加接近實際值。
統(tǒng)計學一般經(jīng)驗認為,當樣本量n≥30或者至少n≥3(k+1)(k為解釋變量的數(shù)目)時,才能滿足模型估計的基本要求。通過實驗分析,對于盜伐林木在現(xiàn)場選擇與盜伐木樹種相同、根徑相近和林相、結(jié)構(gòu)相同或相近林分的立木做為樣木,樣木數(shù)量為n≥3(k+1)進行回歸分析,是可行的。另外,對于同一樣本數(shù)據(jù)和同一顯著性水平,相關系數(shù)R檢驗、F檢驗、t檢驗3種檢驗效果具有等價性。在有關一元線性回歸模型統(tǒng)計檢驗中,完全不必拘泥于選擇哪一種統(tǒng)計檢驗,任選3種統(tǒng)計檢驗中的一種均可以做出同樣的統(tǒng)計判斷。
表7 蓄積量對比
[1] 李潔明,祁新峨.統(tǒng)計學原理[M].上海:復旦大學出版社,2010:1-453.
[2] 龐云海.已采伐林木蓄積量測算的應用實例及分析[J].林業(yè)調(diào)查規(guī)劃,2006(5):13-17.
[3] 郭強,楊杰夫,施海波. Excel中自動完成方差齊性與非齊性t檢驗[J].湖南農(nóng)業(yè)大學學報(社會科學版.素質(zhì)教育研究),2007(1):18-23.
[4] 錢曉莉.回歸方程的檢驗標準[J].財貿(mào)研究,1994(2):66-69.
[5] 王天營. 一元線性回歸分析中三種檢驗的等價性研究[J].統(tǒng)計與決策,2011(3):8-11.
[6] 曹俊忠.一元線性回歸顯著性檢驗方法分析[J].西北紡織工學院學報,1988(3-4):78-82.
[7] 徐悅,陳昌華,蔣之富,等.天然赤松胸徑與樹高相關模型的研究[J].林業(yè)調(diào)查規(guī)劃,2008(3):56-58.
[8] 王曉林,郭斌.柞樹樹高與胸徑相關關系的研究[J].森林工程,2012(6):18-21.
[9] 伐樁現(xiàn)場勘驗和林木蓄積計算方法[EB/OL]. 紅安林業(yè)網(wǎng).
Timber Volume Calculation of Illegal Logging with Unary Linear Regression Equation
SU Jun
(Bulong State-Level Nature Reserve Administration of Xishuangbanna, Jinghong, Yunnan 666100, China)
For the illegal logging cases, in which stump was preserved without cutting wood and surrounded by stand with the same or similar forest form and structure, the linear regression equation was set up to calculate the regression DBH of the illegal logging based on the data of sample trees’ ground diameter and DBH, and the logarithmic curve regression equation was set up to calculate the regression height of the illegal logging based on the data of sample trees’ height. With the principle of Statistics, the related-coefficient test of regression DBH and F-test of regression tree height showed that stock volume of the illegal logging calculated by this method was precision.
illegal logging; DBH regression; tree height regression; unary linear regression; dual volume; related-coefficient test; F-test; timber volume calculation
2017-07-10;
2017-07-17.
粟 軍(1971-),男,云南勐臘人,工程師.從事3S技術(shù)應用、林地林木司法鑒定工作.
10.3969/j.issn.1671-3168.2017.05.003
S711;S758
A
1671-3168(2017)05-0014-06