辛士冬，何 培，姜立春

(東北林業(yè)大學林學院，森林生態(tài)系統(tǒng)可持續(xù)經(jīng)營教育部重點實驗室，黑龍江 哈爾濱 150040)

削度方程可以描繪樹干輪廓，預測樹干任意高度處的直徑，預測樹干任意區(qū)間的材積以及商品材積，也是森林經(jīng)營管理和規(guī)劃、林分生長與收獲的重要工具[1-3]。曾偉生等[4]指出未來材種出材率表的編制將采用削度方程來進行。

樹干形狀的變化可能會受立地質量、林分密度、競爭和間伐等多種因素的影響。由于每一株樹木的干形都不同，在削度方程中引入樹干上部直徑進行矯正，考慮單木之間以及林分之間的變異情況，這種矯正方式約束削度方程必須通過已知的樹干上部直徑位置，使削度方程的預測精度得到提高[5-6]。Czaplewski等[7]使用樹干上部直徑對Max-Burkhart分段削度方程進行矯正，與未矯正的模型相比，其均方根誤差(RMSE)降低了10%～25%。Cao[8]研究發(fā)現(xiàn)利用胸高以上50%高度處的直徑矯正模型最為合理。Kozak[9]利用可變指數(shù)模型提出了使用樹干上部直徑矯正的最佳高度區(qū)間為胸高以上40%～50%的高度；然而，Gómez-García等[10]使用Kozak的可變指數(shù)削度方程發(fā)現(xiàn)利用樹干8 m高度處直徑提供的矯正可以使模型的RMSE從1.55 cm降低到1.34 cm。Sharma等[11]基于黑云杉(Piceamariana)和北美短葉松(Pinusmariana)提出以樹干總高度的34%～38%測量的樹干上部直徑進行矯正最好。

分位數(shù)回歸方法被提出以來，已經(jīng)被廣泛使用于描述響應變量與解釋變量的關系[12]。分位數(shù)回歸估計可以量化響應變量的整個條件分布，并對響應變量不同分位數(shù)點的協(xié)變量效應進行評估，從而為分析數(shù)據(jù)提供了更豐富的信息[13]。落葉松(Larixgmelinii)是中國東北、內蒙古林區(qū)的主要樹種，分布于溫帶及寒溫帶，是東北地區(qū)主要三大針葉用材林樹種之一，所以對天然落葉松干形精準預測的研究具有重要意義[14]。而目前利用分位數(shù)回歸以及系統(tǒng)研究不同矯正位置對樹木干形的影響還鮮見報道，鑒于此，本研究利用非線性分位數(shù)方法建立落葉松分段削度方程，并以樹干相對高20%、30%、40%、50%、60%、70%處及胸徑到樹尖的中間位置(50%*)的直徑對削度方程進行矯正，分析使用不同矯正位置對削度方程預測精度的影響，為落葉松天然林干形的精準預測提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究地點位于大興安嶺松嶺林業(yè)局壯志林場(123°29′～125°11′E，50°9′～51°23′N)，其整個自然實體由伊勒呼里山綿延南伸的兩條低山丘陵組成，西北高、東南低。海拔為400～700 m。地處寒溫帶，屬大陸性氣候。最高氣溫35.5 ℃，最低氣溫-43.8 ℃，年平均氣溫-3 ℃。年降雨量415～500 mm，多集中在7、8月，占全年降水量的48.9%，年無霜期為100～110 d，植物生長期為90～100 d，森林植被由東南向西北明顯過渡。該林場主要樹種包括興安落葉松、白樺(Betulaplatyphylla)、樟子松(Pinussylvestrisvar.mongholica)、山楊(Populusdavidiana)等。

1.2 樣本數(shù)據(jù)的采集

本研究所用的數(shù)據(jù)來源于該林場不同林分的落葉松天然林樣地，樹木被伐后，測量胸徑、樹高及樹干相對高(0、2%、4%、6%、8%、10%、15%、20%、30%、40%、50%、60%、70%、80%、90%)15個位置的直徑。樣地采集樣本214株，剔除2株明顯異常的數(shù)據(jù)，用于建模的樣本數(shù)量為159株，用于檢驗的樣本數(shù)量53株，共計212株，將其以75%和25%隨機分成建模數(shù)據(jù)和獨立檢驗數(shù)據(jù)，落葉松各樣本統(tǒng)計量見表1。

表1 落葉松天然林各樣木調查因子統(tǒng)計量

1.3 基礎削度方程

目前，削度方程在國內外成功應用主要為可變指數(shù)和分段削度方程[15-19]。分段削度方程不僅可以靈活地適應樹干復雜的幾何變化形狀，還可以直接積分計算材積和商品材積，這也是分段削度方程在應用時優(yōu)于可變指數(shù)削度方程的主要方面[20]。因此本研究采用Max and Burkhart分段削度方程作為基礎模型[20]，此方程由3個不同的二次函數(shù)和2個連接點組成，能夠適應樹干上部、中部、下部的幾何形狀，其相應的幾何體為錐體、拋物線體、凹曲線體，對方程進行直接積分可以得到材積方程[21]。具體方程形式如下：

y=a1(x-1)+a2(x2-1)+a3(b1-x)2I1+

a4(b2-x)2I2+ε。

(1)

式中：y=d2/D2，d為樹干高h處的直徑，D為胸徑；x=h/H，H為全樹高，h為從地面到所需位置的高度；當x≤b1時，I1=1；當x>b1時，I1=0; 當x≤b2時,I2=1；當x>b2時，I2=0；a1、a2、a3、a4為待定參數(shù)；b1、b2為樹干下部和上部拐點處的相對高度；ε為模型的誤差項。

1.4 分位數(shù)回歸模型

分位數(shù)回歸是一種估計因變量完整的條件分布和評估不同分位數(shù)預測變量影響的方法，把模型(1)表示成分位數(shù)削度方程如下：

yτ=a1τ(x-1)+a2τ(x2-1)+a3τ(b1τ-x)2I1+

a4τ(b2τ-x)2I2+ετ。

(2)

分位數(shù)回歸的優(yōu)勢在于可以提供任意分位點處解釋變量的變化趨勢，包括數(shù)據(jù)的上下邊界，誤差的分布也不需要任何的假設，并且不受數(shù)據(jù)中異方差的影響，是一種穩(wěn)健的回歸建模和推斷方式[22-24]。分位數(shù)回歸的參數(shù)估計具體是使損失函數(shù)Lmin完成最小值的線性規(guī)劃得以實現(xiàn)：

(3)

式中：τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9，具體采用SAS軟件的 PROC NLP模塊進行擬合[25]。

1.5 模型檢驗

(4)

1.6 模型評價指標

模型的評價采用平均誤差[MAB，式中記為σ(MAB)]、相對誤差[MPB，式中記為σ(MPB)]來對比不同矯正位置及不同檢驗方式模型的預測精度。它們相應的數(shù)學表達式為：

2 結果與分析

2.1 基本模型與分位數(shù)模型的參數(shù)估計

基本模型利用SAS軟件中的PROC MODEL過程進行參數(shù)估計，由于落葉松干形數(shù)據(jù)存在空間相關性，不滿足模型的基本假設條件，所以在參數(shù)估計過程中添加連續(xù)一階自回歸函數(shù)CAR(1)來消除干形數(shù)據(jù)的自相關性[27]。分位數(shù)回歸選取9個不同的分位點(τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)，采用SAS軟件中的PROC NLP過程進行參數(shù)估計，結果見表2?？梢钥闯?，不同分位數(shù)模型的參數(shù)估計值都有所波動，并且Max等[21]分段削度方程在9個分位點下都可以得到參數(shù)估計，說明分位數(shù)回歸方法的靈活性較強，可以提供不同分位點下樹干曲線的變化趨勢，為干形數(shù)據(jù)分析提供了更豐富的信息。

表2 落葉松削度方程的參數(shù)估計值

2.2 落葉松樹干曲線模擬

為了避免分位數(shù)曲線出現(xiàn)不能描述樹干曲線的異常情況，本次試驗利用擬合數(shù)據(jù)和表2的參數(shù)估計值，并以基本模型、中位數(shù)模型、基于3個分位數(shù)組合(τ=0.1、0.5、0.9)、5個分位數(shù)組合(τ=0.1、0.3、0.5、0.7、0.9)、7個分位數(shù)組合(τ=0.1、0.2、0.3、0.5、0.7、0.8、0.9)和9個分位數(shù)組合(τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)為例，模擬了落葉松的分段削度方程在不同分位點的樹干曲線(圖1)?；灸Ｐ?圖1A)和中位數(shù)(τ=0.5)模型(圖1B)只提供了數(shù)據(jù)中心位置的變化趨勢。

圖1 落葉松樹干曲線模擬Fig.1 Stem taper curve simulation of Larix gmelinii

從圖1也可以看出，隨著分位點個數(shù)的增加，基于分位數(shù)回歸的削度方程可以涵蓋的數(shù)據(jù)范圍更大，提供的信息也更加豐富。無論是在數(shù)據(jù)的中間還是數(shù)據(jù)的邊緣，利用分位數(shù)回歸的方法都可以進行模擬，體現(xiàn)了分位數(shù)回歸方法較強的靈活性，可以更為全面地描述樹干形狀的變化趨勢。

2.3 模型預測精度分析

2.3.1 傳統(tǒng)預測檢驗

基于表2中各模型的參數(shù)估計值和53株落葉松獨立檢驗樣本(表1)，利用SAS軟件計算基本模型及各分位數(shù)削度方程的平均誤差(MAB)和相對誤差(MPB)見表3。從表3可以發(fā)現(xiàn)，分位數(shù)(τ=0.5)模型的MAB和MPB略低于基本模型，更適合于預測落葉松樹干曲線的變化趨勢。

表3 落葉松削度方程的獨立性檢驗

2.3.2 不同矯正位置對分位數(shù)組合模型預測精度的影響

使用不同矯正位置的分位數(shù)組合模型，采用Cao等[26]提出的檢驗方法進行檢驗，矯正位置的選取考慮了樹干的形狀變化規(guī)律，其高于胸徑、低于樹干最上部的拐點(即從拋物體轉變?yōu)殄F體的連接點)，所以本次試驗選取樹干相對高20%、30%、40%、50%、60%、70%的直徑和胸徑到樹尖的中間位置(50%*)的直徑對模型進行矯正?；诒?的參數(shù)估計值和公式(4)，采用SAS軟件編寫程序計算不同分位數(shù)組合模型及不同矯正位置的平均誤差(MAB)和相對誤差(MPB)，結果見表4?？梢钥闯?，利用相對高20%和70%處的直徑進行矯正時，各分位數(shù)組合模型的平均誤差和相對誤差高于基本模型的統(tǒng)計量，說明矯正位置選取胸徑和樹干上部拐點附近不會提高模型的預測精度。使用樹干相對高30%～60%處直徑進行矯正的大多數(shù)分位數(shù)組合模型的統(tǒng)計量均優(yōu)于未矯正的基本模型，而且使用樹干不同高度的上部直徑進行矯正時，最優(yōu)的分位數(shù)組合模型并不一致，當使用相對高20%、30%、40%處的直徑和胸徑到樹尖的中間位置(50%*)直徑對分位數(shù)組合模型進行矯正時，其最優(yōu)模型為3個分位數(shù)的組合模型(τ=0.3、0.5、0.7)，而當使用相對高50%、60%、70%的直徑矯正時，其最優(yōu)模型為5個分位數(shù)的組合模型(τ=0.3、0.4、0.5、0.6、0.7)。從表4也可以看出，并不是使用分位點越多的分位數(shù)組合模型預測精度越高，原因在于使用不同分位點的組合，樹干上矯正位置高度h處的直徑測量值位于不同的分位點之間，通過Cao等[26]提出的插值法生成不同的預測曲線，進而得到不同的預測值。本研究的最佳矯正位置為樹干相對高40%處，其各分位數(shù)組合模型的MAB和MPB最低，最優(yōu)模型為3個分數(shù)的組合模型(τ=0.3、0.5、0.7)，與未矯正的基本模型相比，MAB和MPB分別均下降13.5%?？傮w分位數(shù)組合模型的預測精度順序為40% >50%*>50%>60%> 30%> 20%> 70%。

表4 基于不同矯正位置落葉松削度方程的統(tǒng)計量

3 討 論

削度方程廣泛應用于單木材積、林分蓄積的估算和森林碳匯的規(guī)劃，具有重要的實際應用價值，其不僅是林分生長與收獲模型的重要組成部分，也是森林經(jīng)營管理中必不可少的手段，削度方程預測單木材積和林分材積是否精確直接決定了森林規(guī)劃系統(tǒng)方案的科學性。本研究利用一個樹干上部直徑對分位數(shù)分段削度方程進行矯正，選取的矯正位置在高于胸徑位置并低于樹干上部拐點位置(即樹干中間拋物體轉變?yōu)闃涓身敳垮F體的直徑)。最低矯正位置的選擇是由數(shù)據(jù)中最小樹高決定的，本研究數(shù)據(jù)中胸徑高度為1.3 m，最小樹高為7.1 m，其胸徑高度占總樹高的18.3%。最高的矯正位置要低于樹干上部拐點，主要由模型的參數(shù)估計值決定，基本模型的樹干上部拐點參數(shù)為0.789 2，基于不同分位數(shù)的樹干上部拐點參數(shù)范圍為0.791 3～0.862 9，所以本研究矯正位置選取樹干相對高度的20%～70%(以10%為間隔)。選取胸徑到樹尖的中間(50%*)處的樹干上部直徑作為矯正位置，是因為有研究者指出這一矯正位置為最佳矯正位置[26]。

從不同矯正位置的落葉松各分位數(shù)組合模型的統(tǒng)計指標可以發(fā)現(xiàn)，準確選擇樹干上部直徑的矯正位置至關重要，不同的矯正位置對模型的精度影響較大。本研究的最佳矯正位置為樹干相對高40%處，這與Cao等[26]研究美國路易斯安納州的火炬松(P.taedaL.)所提出的最佳矯正位置并不一致，Cao等[26]提出5個分位數(shù)組合模型的預測精度最高，并且只研究對比了3個分位數(shù)組合(τ=0.1、0.5、0.9)和5個分位數(shù)組合(τ=0.1、0.3、0.5、0.7、0.9)模型。本次試驗考慮了多種分位數(shù)組合，并構建了3、5、7、9個分位數(shù)的組合模型，但結果表明3個分位數(shù)組合模型(τ=0.3、0.5、0.7)的預測精度最高，該組合模型與未矯正的模型相比，MAB和MPB均下降13.5%，可以為該地區(qū)落葉松樹干干形提供相對精準的預測。

4 結 論

1)基于中位數(shù)回歸(τ=0.5)的削度模型預測精度略優(yōu)于非線性回歸的削度模型。

2)在削度模型中引入樹干相對高40%處的直徑并采用3個分位數(shù)的組合模型(τ=0.3、0.5、0.7)進行矯正，預測能力明顯得以提高。

3)本研究選用Max and Burkhart削度方程，是因為該方程具有明顯的幾何生物學意義以及能直接積分計算各種材積[20]，當選用不同的削度方程進行矯正時，分位數(shù)組合以及矯正位置都可能發(fā)生變化。