不同矯正位置對落葉松分位數(shù)削度方程預(yù)測精度的影響

辛士冬，何 培，姜立春

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

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

樹干形狀的變化可能會受立地質(zhì)量、林分密度、競爭和間伐等多種因素的影響。由于每一株樹木的干形都不同，在削度方程中引入樹干上部直徑進(jìn)行矯正，考慮單木之間以及林分之間的變異情況，這種矯正方式約束削度方程必須通過已知的樹干上部直徑位置，使削度方程的預(yù)測精度得到提高[5-6]。Czaplewski等[7]使用樹干上部直徑對Max-Burkhart分段削度方程進(jìn)行矯正，與未矯正的模型相比，其均方根誤差(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%測量的樹干上部直徑進(jìn)行矯正最好。

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

1.3 基礎(chǔ)削度方程

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

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

a4(b2-x)2I2+ε。

(1)

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

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

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

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

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

(2)

分位數(shù)回歸的優(yōu)勢在于可以提供任意分位點處解釋變量的變化趨勢，包括數(shù)據(jù)的上下邊界，誤差的分布也不需要任何的假設(shè)，并且不受數(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模塊進(jìn)行擬合[25]。

1.5 模型檢驗

(4)

1.6 模型評價指標(biāo)

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

2 結(jié)果與分析

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

基本模型利用SAS軟件中的PROC MODEL過程進(jìn)行參數(shù)估計，由于落葉松干形數(shù)據(jù)存在空間相關(guān)性，不滿足模型的基本假設(shè)條件，所以在參數(shù)估計過程中添加連續(xù)一階自回歸函數(shù)CAR(1)來消除干形數(shù)據(jù)的自相關(guān)性[27]。分位數(shù)回歸選取9個不同的分位點(τ=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9)，采用SAS軟件中的PROC NLP過程進(jìn)行參數(shù)估計，結(jié)果見表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ù)回歸的方法都可以進(jìn)行模擬，體現(xiàn)了分位數(shù)回歸方法較強的靈活性，可以更為全面地描述樹干形狀的變化趨勢。

2.3 模型預(yù)測精度分析

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

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

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

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

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

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

3 討 論

削度方程廣泛應(yīng)用于單木材積、林分蓄積的估算和森林碳匯的規(guī)劃，具有重要的實際應(yīng)用價值，其不僅是林分生長與收獲模型的重要組成部分，也是森林經(jīng)營管理中必不可少的手段，削度方程預(yù)測單木材積和林分材積是否精確直接決定了森林規(guī)劃系統(tǒng)方案的科學(xué)性。本研究利用一個樹干上部直徑對分位數(shù)分段削度方程進(jìn)行矯正，選取的矯正位置在高于胸徑位置并低于樹干上部拐點位置(即樹干中間拋物體轉(zhuǎn)變?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)計指標(biāo)可以發(fā)現(xiàn)，準(zhǔn)確選擇樹干上部直徑的矯正位置至關(guān)重要，不同的矯正位置對模型的精度影響較大。本研究的最佳矯正位置為樹干相對高40%處，這與Cao等[26]研究美國路易斯安納州的火炬松(P.taedaL.)所提出的最佳矯正位置并不一致，Cao等[26]提出5個分位數(shù)組合模型的預(yù)測精度最高，并且只研究對比了3個分位數(shù)組合(τ=0.1、0.5、0.9)和5個分位數(shù)組合(τ=0.1、0.3、0.5、0.7、0.9)模型。本次試驗考慮了多種分位數(shù)組合，并構(gòu)建了3、5、7、9個分位數(shù)的組合模型，但結(jié)果表明3個分位數(shù)組合模型(τ=0.3、0.5、0.7)的預(yù)測精度最高，該組合模型與未矯正的模型相比，MAB和MPB均下降13.5%，可以為該地區(qū)落葉松樹干干形提供相對精準(zhǔn)的預(yù)測。

4 結(jié) 論

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

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

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