基于線性分位數(shù)組合的興安落葉松冠幅預測

王君杰， 姜立春

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

樹冠大小作為樹木的基本屬性，與樹冠表面積和體積呈顯著正相關，影響樹木的光合和呼吸作用[1-2]，進而影響樹木活力[3]、林木競爭[4]、林分健康[5]及生物圈中的碳循環(huán)[6]，是預測森林動態(tài)的重要指標[7]。冠幅(crown width)是判斷樹冠大小的主要指標，常作為重要變量引入生物量模型[8]、削度方程[9]和樹干材積[10]等單木模型中，此外冠幅大小影響林冠的光截獲量[11]，從而影響林下的生物多樣性和野生動植物分布[12]。冠幅用途廣泛，但測量費時費力，因此林業(yè)上經(jīng)常采用相關模型對冠幅進行預測[2]。冠幅大小受立地條件、林分密度和林分中各種隨機因素的影響[13]，通常將胸徑和其他單木或林分變量引入到基礎冠幅模型中，以降低模型預測誤差[1-2,7,12-13]。

單木冠幅數(shù)據(jù)通常來自樣地的每木檢尺或重復測量，常呈層次嵌套結(jié)構[14]，這導致構建冠幅模型時存在異方差和自相關等問題。雖然可以通過加權回歸和引入自相關函數(shù)來矯正這些問題，但會增加模型的擬合難度，有時矯正效果并不理想。分位數(shù)回歸(quantile regression)[15]對數(shù)據(jù)以及模型誤差項的假設條件比較寬松，保證了參數(shù)估計的穩(wěn)健性和可靠性，同時也可以評估不同分位數(shù)對響應變量的影響[16-17]。近年來，分位數(shù)回歸在林業(yè)中得到一定的應用，但這些大多都是研究自變量對因變量在某個特定分位數(shù)水平上影響的邊際效應或平均效果，如模擬自稀疏邊界線[16]、直徑分布[17]、樹高曲線[18]、潛在最大冠幅[19]、削度方程[20]等。雖然一些學者使用不同分位數(shù)組合對胸徑生長[21]和削度方程[22-23]進行預測，但由于方法的限制，都只探討了一個抽樣數(shù)量對模型預測精度的影響。雖然?z?elik等[24]提出使用多個抽樣數(shù)據(jù)和分位數(shù)組合對樹高曲線進行預測，但沒有深入研究不同抽取方案 (隨機抽樣、選取最大樹、平均木和最小樹)對預測結(jié)果的影響。

興安落葉松(Larixgmelinii)是大興安嶺地區(qū)的主要樹種，林木蓄積豐富，被廣泛用于木材生產(chǎn)和造林，同時在區(qū)域碳儲存和碳循環(huán)中起著重要的作用。因此常常需要利用模型對興安落葉松的冠幅進行準確預測。目前綜合應用分位數(shù)組合和不同抽取方案對興安落葉松冠幅模型的研究還未見報道。因此，本研究使用傳統(tǒng)方法構建大興安嶺落葉松基礎和多元冠幅模型，評估影響冠幅變化的主要因素；利用分位數(shù)回歸和分位數(shù)組合對基礎和多元模型進行擬合和預測，對比三分位數(shù)、五分位數(shù)、七分位數(shù)和九分位數(shù)組合的預測能力，評價傳統(tǒng)回歸、分位數(shù)回歸和分位數(shù)組合的預測效果；評估模型實際應用時的抽取方案，以期為野外實際應用提供指導。

1 材料與方法

1.1 研究地概況

研究地點位于黑龍江省大興安嶺新林區(qū)(123°41′～125°25′E，51°20′～52°10′N)。新林區(qū)地處大興安嶺山脈東坡，伊勒呼里山北坡。氣候寒冷濕潤，氣溫年較差和日較差較大，年平均氣溫為-2.6 ℃。降水多集中于7—8月，年降水量為513.9 mm。8月下旬開始出現(xiàn)初霜，無霜期平均為90 d左右。風速較小，年平均風速一般為2～3 m/s。新林區(qū)水系為外流流域，河谷密集，年平均秒流量20.95 m3/s。植被屬于寒溫帶針葉林區(qū)的大興安嶺山地寒溫針葉林帶，森林覆蓋率達84%，但植物種類相對貧乏。主要喬木樹種有興安落葉松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongouca)、白樺(Betulaplatyphylla)、山楊(Populusdavidiana)等。

1.2 數(shù)據(jù)來源及統(tǒng)計

研究用數(shù)據(jù)為2012年和2017年調(diào)查的65塊落葉松天然林樣地(各樣地落葉松斷面積比例≥65%)的3 452株落葉松。樣地涵蓋廣泛的林分結(jié)構和立地條件，面積大小取決于林分密度，范圍為0.02～0.10 hm2。實測每個樣地胸徑大于5 cm的所有樹木的胸徑(DBH，記為DBH)、樹高(HT)、枝下高(HCB)，并記錄每株樹的相對位置。測量東、南、西、北4個方向的最大樹冠半徑，得到東西、南北兩個方向的樹冠寬度，冠幅(CW，記為WC)計算為兩個方向樹冠寬度的算術平均值[1]。除上述變量，還計算了每株樹的冠長率(RC)、高徑比(RHD)、樣地平均樹高(Hmean，m)、優(yōu)勢木平均高(HDOM，m)、優(yōu)勢木平均胸徑(DDOM，cm)、每公頃斷面積(AB，m2/hm2)、每公頃株數(shù)(N，株/hm2)、相對間距指數(shù)(SR)、樣地內(nèi)所有大于對象木胸徑的立木斷面積總和(BAL,m2)、樣地內(nèi)落葉松斷面積(Glarch,m2)。所有單木及林分變量統(tǒng)計信息見表1。

表1 落葉松天然林各樣木調(diào)查因子數(shù)據(jù)Table 1 Deseriptive statistics for Larix gmelinii sample trees

1.3 基礎冠幅模型選擇

根據(jù)興安落葉松冠幅和胸徑的關系，即冠幅隨胸徑的增大而增大，本研究選用了6個生長模型：線性、邏輯斯蒂(Logistic)、冪函數(shù)、理查德(Richards)和威布爾(Weibull)函數(shù)來模擬落葉松冠幅和胸徑的關系。具體模型形式見表2。

表2 候選冠幅基礎模型Table 2 Candidate basic crown width models

1.4 多元冠幅模型構建

選用胸徑作為唯一預測因子的基礎冠幅模型不能反映其他單木因子和林分變量對冠幅的影響[1]，因此本研究考慮引入其他影響冠幅的因子，如樹木大小、立地條件和林分競爭等變量[2,7,13-14]?；谧顑?yōu)基礎冠幅模型構建含有其他變量的多元冠幅模型，采用多元逐步回歸技術，同時考慮冠幅與其他自變量的相關關系，篩選模型中的其他自變量，包括表1中提到的變量及變量的變形形式，如對數(shù)轉(zhuǎn)換和開方等，同時考慮多變量的多重共線性，剔除方差膨脹因子(VIF)大于5的自變量[25]。

1.5 分位數(shù)回歸

分位數(shù)回歸可以對數(shù)據(jù)的任意分位點進行響應變量的估計。第τ個分位點的參數(shù)估計可以通過最小化殘差絕對值的非對稱損失函數(shù)來估計[16]。如式(1)所示：

(1)

1.6 分位數(shù)組合應用方法

本研究主要是應用分位數(shù)組合對冠幅模型進行預測，測試了7種分位數(shù)組合：三分位數(shù)組合 (τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)、五分位數(shù)組合 (τ=0.1, 0.3, 0.5, 0.7, 0.9和τ=0.3, 0.4, 0.5, 0.6, 0.7)、七分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 0.9和τ=0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.9)、九分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)[21,26]。

(2)

1.7 分位數(shù)組合抽取方案

應用分位數(shù)組合時需要使用抽樣樹木的數(shù)據(jù)對模型進行樣地局部化校準，從而預測未抽樣樹木的冠幅。采用以下4種抽取方案得到抽樣數(shù)據(jù)：① 每個樣地隨機抽取1～9株樹。為了得到可靠的結(jié)果，此過程重復100次，然后計算平均值。② 選取每個樣地胸徑最大的1～9株樹。③ 選取每個樣地與該樣地平均胸徑最接近的1～9株樹。④ 選取每個樣地胸徑最小的1～9株樹。

1.8 模型評價和檢驗指標

本研究使用K折交叉驗證(K-fold cross validation，K=5)對模型進行檢驗[27]，將數(shù)據(jù)隨機分為5組，每組包含13塊樣地。使用4組數(shù)據(jù)估計模型參數(shù)，對剩余1組數(shù)據(jù)預測，該過程重復5次，計算各模型評價指標。使用平均誤差[ME，式中記為σ(ME)]、均方根誤差[RMSE，式中記為σ(RMSE)]、和決定系數(shù)(R2)選擇最優(yōu)基礎模型。使用平均絕對誤差[MAE，記為σ(MAE)]和均方根誤差對各模型擬合優(yōu)度進行評價。使用平均百分比誤差[MPE，式中記為σ(MPE)]、平均絕對百分比誤差[MAPE，式中記為σ(MAPE)]、均方根誤差和決定系數(shù)(R2)對模型K折交叉驗證結(jié)果進行評價。各評價指標相應的數(shù)學表達式為：

(3)

(4)

R2=1-

(5)

(6)

(7)

(8)

使用R軟件的lm過程擬合基礎和多元冠幅模型，用quantreg包擬合線性分位數(shù)回歸(τ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)，使用SAS的Proc NLP過程計算組合分位數(shù)的插值系數(shù)。

2 結(jié)果與分析

2.1 基礎、多元和分位數(shù)冠幅模型確定

利用R軟件的lm和nls過程分別擬合候選基礎模型的線性和非線性形式，擬合結(jié)果見表3。

表3 候選冠幅基礎模型擬合統(tǒng)計量Table 3 Fitting statistics with candidate basic crown width models

由表3可以看出，6個模型評價指標都很相近。線性模型參數(shù)估計穩(wěn)定，不受初始值影響，且通常可以很好地擬合冠幅[7,25]，因此選擇線性模型作為最優(yōu)基礎冠幅模型。模型形式如下：

WCij=β0+β1DBHij+εij。

(9)

式中：WCij、DBHij分別為第i個樣地第j株樹的冠幅和胸徑，β0、β1為模型參數(shù)，εij為模型誤差項。

基于最優(yōu)基礎冠幅模型，采用逐步回歸技術構建多元冠幅模型，最終確定模型包括的變量為：單木枝下高HCB(樹木大小)、樣地平均高Hmean(立地質(zhì)量)和樣地內(nèi)落葉松斷面積Glarch(競爭)。具體模型形式如下：

WCij=β0+β1DBHij+β2HCBij+β3Hmeani+β4Glarch,i+εij。

(10)

式中：HCBij為第i個樣地第j株樹的枝下高；Hmeani和Glarch,i分別為第i個樣地的平均高和落葉松斷面積；β0、β1、β2、β3、β4為模型參數(shù)。

將基礎和多元冠幅模型分別表示為分位數(shù)回歸模型，模型形式如下：

(11)

(12)

2.2 模型的參數(shù)估計

利用R軟件的lm過程和quantreg包分別對基礎和多元冠幅模型進行線性和分位數(shù)回歸擬合。線性回歸和分位數(shù)回歸的參數(shù)估計值和擬合統(tǒng)計量見表4，所有參數(shù)均極顯著(P<0.001)。多元冠幅模型(10)的參數(shù)β1和β3的估計值為正數(shù)，β2和β4的參數(shù)估計值為負數(shù)，說明冠幅與胸徑和樣地平均高成正比，與枝下高和落葉松斷面積成反比。與基礎模型線性回歸相比，多元冠幅模型線性回歸的MAE和RMSE分別降低6.9%和6.3%?；A和多元分位數(shù)模型在不同分位點的參數(shù)估計值都顯著。分位數(shù)回歸的高分位點和低分位點的誤差MAE和RMSE變化較大，離中位數(shù)越遠，擬合效果越差。兩個中位回歸(τ=0.5)擬合效果最好，且評價指標與相應模型的線性回歸差異不大，可選為最優(yōu)分位數(shù)回歸模型。

表4 線性回歸和分位數(shù)回歸的參數(shù)估計和擬合統(tǒng)計量Table 4 Parameter estimation and fitting statistics for linear regression and quantile regression at nine quantiles

表4(續(xù))

為了直觀地分析冠幅隨胸徑的變化，圖1繪制了線性回歸和基于9個分位點(τ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)的基礎冠幅模型和多元冠幅模型的模擬曲線(圖1)?？梢钥闯?，分位數(shù)回歸可以提供不同分位點的估計，產(chǎn)生了一簇冠幅預測直線，可以體現(xiàn)整個分布的信息。中位數(shù)回歸與線性回歸預測幾乎重合，較其他分位數(shù)回歸能更好地擬合數(shù)據(jù)的平均趨勢。

圖1 線性回歸和基于9個分位點的分位數(shù)回歸的基礎和多元模型的模擬曲線Fig.1 Basic and multivariate model simulation of linear regression and quantile regression based on nine quantiles

2.3 模型比較及組合分位數(shù)抽取方案和數(shù)量的確定

計算線性和分位數(shù)回歸各分位數(shù)的預測統(tǒng)計量比較簡單，就是把參數(shù)估計值直接代入相應模型得到預測值，然后計算各模型的統(tǒng)計量。然而，分位數(shù)組合的統(tǒng)計量計算需要確定：分位數(shù)組合、抽取方式、抽取數(shù)量。本研究分別對比了三分位數(shù)組合 (τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)、五分位數(shù)組合 (τ=0.1, 0.3, 0.5, 0.7, 0.9和τ=0.3, 0.4, 0.5, 0.6, 0.7)、七分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 0.9和τ=0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.9)和九分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)，采用4種抽取方案(隨機抽樣、選取最大樹、平均木、最小樹)和9種抽取數(shù)量(1～9株)對冠幅進行預測，分位數(shù)組合預測誤差RMSE見圖2。為了比較分位數(shù)組合與中位數(shù)回歸(最優(yōu)分位數(shù))的預測效果，中位數(shù)回歸模型預測誤差用橫線表示?？梢钥闯隼梅治粩?shù)組合和抽樣數(shù)據(jù)對模型進行校準可以提高模型的預測能力，且抽取方案和數(shù)量的不同導致分位數(shù)組合的預測能力產(chǎn)生明顯差異。7種分位數(shù)組合方式在相同抽取方案和樣本大小下的差異都很小，特別是對于隨機抽樣和選取平均木。對于選取最大樹和最小樹，7種分位數(shù)組合方式略有差異，大多數(shù)情況下，七分位數(shù)組合和九分位數(shù)組合的RMSE都比三分位數(shù)和五分位數(shù)組合大?？傮w上看，三分位數(shù)組合的預測能力較好，因此本研究選用三分位數(shù)組合(τ=0.1,0.5,0.9和τ=0.3, 0.5, 0.7)對冠幅進行預測。

圖2 模型預測誤差比較Fig.2 Comparison of the prediction errors of the models

隨著抽取數(shù)量的增加，模型的RMSE越來越小，呈現(xiàn)出先陡后緩的趨勢，對于4種抽取方案，選取最小樹的預測能力不理想，基礎和多元分位數(shù)組合的抽樣數(shù)量分別要超過8株(圖2G)和6株(圖2H)才能優(yōu)于中位數(shù)回歸。從RMSE值來看，其他3種抽取方案的預測能力排序均為：選取最大樹>隨機抽樣>選取平均木。對于基礎和多元分位數(shù)組合，樣本大小分別超過3株和4株時，RMSE即可小于中位數(shù)回歸。當使用6株最大樹進行校準時，與中位數(shù)回歸相比，基礎和多元分位數(shù)組合的RMSE分別降低了16.3%和11.6%。雖然選擇更多數(shù)量的最大樹可以進一步降低RMSE，但RMSE的降低速率減緩，且會增加更多采樣成本。每個樣地測量6株樣本可以達到采樣成本和預測精度之間的平衡，因此，對于基礎和多元分位數(shù)組合，建議抽取6株最大樹方案。

2.4 模型K折交叉驗證評價

分別對基礎和多元冠幅模型的線性回歸、中位數(shù)回歸(最優(yōu)分位數(shù)回歸)和三分位數(shù)組合(τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)進行K折交叉驗證，分位數(shù)組合的抽取方案采用選取最大樹6株，評價結(jié)果見表5?？梢钥闯龆嘣Ｐ偷腗PE、MAPE、RMSE和R2都優(yōu)于相應的基礎模型。相對于基礎模型，多元模型的線性和中位數(shù)回歸的評價指標有明顯改善，多元分位數(shù)組合的改善不明顯?；A模型中，線性回歸與中位數(shù)回歸的差異不大，2種三分位數(shù)組合(τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)的檢驗結(jié)果相似。多元模型中，中位數(shù)回歸略優(yōu)于線性回歸，分位數(shù)組合τ=0.3, 0.5, 0.7的預測能力略優(yōu)于τ=0.1, 0.5, 0.9。無論是基礎還是多元模型，分位數(shù)組合的檢驗結(jié)果都明顯優(yōu)于線性回歸和中位數(shù)回歸，具有更好的預測能力。

表5 線性回歸、中位數(shù)回歸和三分位數(shù)組合的K折交叉驗證統(tǒng)計量Table 5 K-fold cross validation statistics of linear regression, median quantile regression and three quantiles groups

3 討 論

本研究基于大興安嶺東部興安落葉松天然林數(shù)據(jù)，構建了基礎和多元冠幅模型，并利用線性回歸和分位數(shù)回歸對冠幅模型進行擬合，利用分位數(shù)組合對冠幅進行預測，同時考慮不同抽取方案和數(shù)量對預測結(jié)果的影響。研究表明，線性模型可較好地描述樹木冠幅與胸徑的關系，這與雷相東等[25]和Bechtold[28]的研究結(jié)果一致。加入單木枝下高HCB(樹木大小)、樣地平均高Hmean(立地質(zhì)量)和樣地內(nèi)落葉松斷面積Glarch(競爭)的多元冠幅模型可以明顯提高模型的擬合和預測能力。冠幅隨著胸徑和樣地平均高的增大而增大，隨著枝下高和落葉松斷面積的增大而減小。Fu等[2]和Sharma等[29-30]分別對杉木(Cunninghamialanceolata)、挪威云杉(Piceaabies)和歐洲山毛櫸(Fagussylvatica)的冠幅模型進行研究，得到了類似的結(jié)果，即立地質(zhì)量越好，枝下高和競爭越小，則樹木冠幅越大。樣地平均高是描述林分立地質(zhì)量的指標[29]，樣地平均高越大，林分的生長條件越好，則冠幅越大[2,29-30]。單木枝下高與樹冠大小(如冠長率)密切相關，當樹冠底部的枝條死亡時，枝下高增大，冠幅變小[2,29-30]。樹冠的生長對林分密度非常敏感。樣地落葉松斷面積是一種與密度相關的競爭指標[5]，對于同一林分，落葉松斷面積越大，對同一空間資源利用的競爭越大，則冠幅越小[30]。

分位數(shù)回歸能靈活地反映冠幅的分布，中位數(shù)回歸與其他分位點相比擬合效果更好[18]，且與線性回歸基本一致。Bohora等[21]利用分位數(shù)組合對胸徑生長進行預測時，發(fā)現(xiàn)三分位數(shù)(τ=0.1, 0.5, 0.9)、五分位數(shù)(τ=0.1, 0.3, 0.5, 0.7, 0.9)和九分位數(shù)組合都可以提高預測精度，但九分位數(shù)組合會導致預測誤差增大。Cao等[22]、?z?elik等[23-24]和馬巖巖等[26]分別對削度方程和樹高曲線進行研究，發(fā)現(xiàn)五分位數(shù)組合略優(yōu)于三分位數(shù)組合，但都沒有考慮抽取方案對分位數(shù)組合的影響。本研究將分位數(shù)組合形式擴展至7種：三分位數(shù)組合 (τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)、五分位數(shù)組合 (τ=0.1, 0.3, 0.5, 0.7, 0.9和τ=0.3, 0.4, 0.5, 0.6, 0.7)、七分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 0.9和τ=0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.9)和九分位數(shù)組合(τ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)，同時考慮不同抽取方案的影響，發(fā)現(xiàn)7種分位數(shù)組合都可以提高冠幅模型的預測精度，但各分位數(shù)組合間的差異會隨著抽取方案的不同產(chǎn)生變化。對于選取最大樹和最小樹產(chǎn)生的抽樣數(shù)據(jù)，7種分位數(shù)組合產(chǎn)生較大差異，而對于隨機抽樣和選取平均木差異不大?？傮w上，三分位數(shù)組合 (τ=0.1, 0.5, 0.9和τ=0.3, 0.5, 0.7)的預測能力較好且穩(wěn)定，增加分位數(shù)組合個數(shù)不能使預測結(jié)果有實質(zhì)性的提高。

分位數(shù)組合應用時使用1個已知觀測的抽樣僅適合部分單木模型預測[21-23,26]，對于其他單木模型，使用1個觀測值反而會增大預測誤差[24]，本研究對冠幅模型預測的研究也得到相同的結(jié)果。?z?elik等[24]利用多個抽樣數(shù)據(jù)和分位數(shù)組合對樹高曲線進行預測，但沒有對抽取方案進行比較。本研究使用4種抽取方案(隨機抽樣、選取最大樹、平均木和最小樹)對分位數(shù)組合的預測能力進行比較，發(fā)現(xiàn)抽取方案的不同會使預測結(jié)果產(chǎn)生明顯差異。對于4種抽取方案，分位數(shù)組合的預測效果都會隨著抽樣數(shù)量的增加而提高，但選取最小樹的效果不理想，這可能是因為小樹的代表性較差，不能充分反映樣地的信息[31]。對于基礎和多元分位數(shù)組合，選取最大樹的預測能力最好。選取最大樹的另一個好處是，可以對材積較大的樹木進行更精確的校準[32]。在平衡了模型預測精度和抽樣成本之后，建議應用分位數(shù)組合進行預測時每塊樣地抽樣數(shù)量為6株。當無法獲得抽樣數(shù)據(jù)時，可以選擇線性回歸或中位數(shù)回歸對冠幅進行預測，多元模型的中位數(shù)回歸略優(yōu)于線性回歸。