王君杰，夏宛琦，姜立春

（東北林業(yè)大學 林學院，黑龍江 哈爾濱 150040）

樹高-胸徑模型在計算樹干材積、出材率表編制、生長和收獲模型、森林資源調(diào)查、立地質(zhì)量評價、生物量和碳儲量研究等方面有著廣泛的應用[1-5]，目前，樹高-胸徑模型的模型形式主要有非線性模型和混合效應模型。非線性模型采用普通最小二乘法（OLS）得到參數(shù)估計[6-7]，但只能得到樹高的平均預測。當數(shù)據(jù)來自較大區(qū)域，即立地條件不同時，總體使用最小二乘法對一些區(qū)域可能產(chǎn)生較大誤差。此外，最小二乘法假定模型誤差服從獨立正態(tài)分布，當違背這種假設時，得到的參數(shù)估計不可靠?；旌闲Ｐ屯瑫r考慮固定效應參數(shù)和隨機效應參數(shù)，使用固定效應參數(shù)可以得到樹高的平均預測，如果能夠計算隨機效應參數(shù)值就可以進行個體預測[8-9]。個體預測需要進行二次抽樣校正隨機參數(shù)。二次抽樣對于小尺度是可行的，但對于大尺度區(qū)域進行二次抽樣會大大增加調(diào)查成本，且會增大樹高測量的誤差，在實際應用時限制了混合模型個體預測的應用，而混合模型總體平均預測效果有時不如傳統(tǒng)的最小二乘法[10-11]。此外混合效應模型應用數(shù)據(jù)條件必須是分組數(shù)據(jù)或縱向數(shù)據(jù)，如果數(shù)據(jù)沒有這樣的信息，如結(jié)合森林采伐獲得的伐倒木實測樹高和胸徑數(shù)據(jù)，沒有樣地信息，無法應用混合效應模型進行分析，因此不能進行混合效應樹高-胸徑模型的構(gòu)建。

分位數(shù)回歸（Quantile regression）由Koenker和Bassett[12]提出，它可以給出一組用來描述響應變量的分位數(shù)回歸曲線。不同回歸曲線具有各種形狀，可以靈活地反映數(shù)據(jù)的變化規(guī)律。近年來，分位數(shù)回歸方法在林業(yè)上已經(jīng)用于構(gòu)建和研究自稀疏邊界線模型[13]、胸徑百分位模型[14]、胸徑生長模型[15]、削度方程[16]和最大密度線模型[17]等。

啞變量（Dummy variable）模型是將一個分類變量引入基本模型中，僅擬合一個模型就可以得到不同分類變量模型的參數(shù)估計值。在林業(yè)中常用于不同區(qū)域模型的構(gòu)建，并與混合效應模型等方法共同使用，提高模型精度[9,18]。已有很多學者對啞變量方法與線性混合效應模型[19]、非線性混合效應模型[20-21]等其他方法進行比較，但與分位數(shù)回歸模型的比較還未見報道。

本研究以大興安嶺4個區(qū)域的興安落葉松Larix gmelinii為研究對象，采用分位數(shù)回歸方法和啞變量方法構(gòu)建各區(qū)域的樹高-胸徑模型，確定各區(qū)域樹高-胸徑模型的回歸系數(shù)，然后對分位數(shù)回歸模型、啞變量模型和傳統(tǒng)的基本模型進行詳細的對比分析，為林業(yè)樹高-胸徑模型的構(gòu)建提供新的思路和方法。

1 研究區(qū)概況

大興安嶺地區(qū)是我國天然林和重點國有林的主要分布地區(qū)之一，屬大陸寒溫帶季風性針葉林氣候，平均海拔573 m，年平均氣溫-2～4 ℃，年降水量450 mm左右。物種資源豐富，主要樹種為興安落葉松，其次是樟子松Pinus sylvestrislvar.mongouca、云杉Picea koraiensis、白樺Berula platypgylla等。本區(qū)屬高緯度山地，地勢西高東低，主要土壤類型為棕色針葉林土、暗棕壤、暗色草甸土和沼澤土等。地貌基本屬于低山、丘陵地貌，地形起伏不大。全區(qū)河谷寬闊，江河網(wǎng)布，水資源豐富，存在大量的沼澤地和草甸地貌，年徑流量達149億m3以上。

根據(jù)張萬儒等[22]的方法，將研究地區(qū)劃分為4個區(qū)域：漠河、圖強和阿木爾劃分為區(qū)域1（西北）；塔河、十八站和韓家園劃分為區(qū)域2（東北）；松嶺、加格達奇劃分為區(qū)域3（中部）；呼中、新林林場劃分為區(qū)域4（南部）。分區(qū)情況見圖1。在4個區(qū)域?qū)崪y了2 411株興安落葉松伐倒木的胸徑和樹高。為了使數(shù)據(jù)具有代表性，每個區(qū)域樹高和胸徑數(shù)據(jù)都盡量接近正態(tài)分布。采用隨機的方法將每個區(qū)域的數(shù)據(jù)分成建模數(shù)據(jù)（80%）和檢驗數(shù)據(jù)（20%），即建模樣本數(shù)據(jù)共1 933對，檢驗樣本數(shù)據(jù)共478對。各區(qū)域的建模樣本和檢驗樣本數(shù)據(jù)統(tǒng)計量基本情況見表1。

2 研究方法

2.1 基本模型

圖1 研究區(qū)域分區(qū)情況Fig.1 Partition of the study area

選取國內(nèi)外林業(yè)上廣泛應用的12個樹高-胸 徑 方 程：Weibull方 程、Schumacher方 程、Richards方程、Korf方程、Loetsh方程、Curtis方程、Mitscherlich方程、唐守正方程、Sibbesen方程、Wykoff方程、Hossfeld方程、Ratkowsky方程[3,23]。初步分析表明Richards方程表現(xiàn)了較好的擬合效果，因此本研究選取Richards方程作為基本模型，模型形式如下：

式中：H為樹高；D為胸徑；a、b、c為模型參數(shù)；ε為誤差項。

2.2 分位數(shù)回歸

分位數(shù)回歸（Quantile regression）利用自變量的多個分位數(shù)得到因變量條件分布的相應分位數(shù)方程，可以捕捉分布的尾部特征，對條件均值以外的預測非常有效[24]。當自變量對不同部分的因變量的分布產(chǎn)生不同的影響時，例如出現(xiàn)左偏或右偏的情況時，更加全面地刻畫分布的特征，從而得到全面的分析[25]。

表1 興安落葉松樣本統(tǒng)計量Table 1 Summary statistics for sample trees of Larix gmelinii

基于Richards方程，其9個分位數(shù)（τ=0.1, 0.2,0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9）樹高-胸徑模型形式如下：

式中：Hτ是第τ個分位數(shù)的樹高預測值；aτ、bτ、cτ是第τ個分位數(shù)模型的參數(shù)，其他參數(shù)定義同上。

分位數(shù)模型中第τ個分位數(shù)模型的參數(shù)估計可以通過最小化殘差絕對值的非對稱損失函數(shù)來估計[13]，如式（3）所示：

不同分位數(shù)模型擬合具體采用SAS軟件中的非線性規(guī)劃（NLP）模塊計算不同分位數(shù)τ上的樹高-胸徑方程的參數(shù)估計。

2.3 啞變量模型

啞變量（Dummy variable）模型是在基礎(chǔ)模型中引入一個定性變量，將不同區(qū)域的生長模型整合成一個模型，進行參數(shù)估計時僅對啞變量模型進行擬合即可得到各區(qū)域的參數(shù)估計值，既減少了工作量，又解決了不同區(qū)域單獨建立模型所造成的整體和分量之間不相容的問題[26]。啞變量法不同于最小二乘法將所有參數(shù)都視為全局參數(shù)，也不同于混合模型方法將模型參數(shù)視為適用于總體的固定參數(shù)和適用于個體的隨機參數(shù)的和，它把個體參數(shù)看作是固定的，但在不同的個體之間卻是不同的[20]，Wang等[27]稱之為固定個體效應。

Richards方程在不同區(qū)域都表現(xiàn)了顯著的不同，因此本研究也構(gòu)建了不同區(qū)域的啞變量樹高-胸徑模型，基于Richards方程的啞變量模型如下：

式中：a1、a2、a3、b1、b2、b3、c1、c2、c3是與區(qū)域有關(guān)的參數(shù)值；S1、S2、S3為啞變量，其他參數(shù)定義同上。當區(qū)域為區(qū)域1時，S1=1，而S2、S3為0，依次類推。當區(qū)域為區(qū)域4時，S1、S2、S3都為0。

采用非線性額外平方和法比較不同區(qū)域間的差異[28]，如果沒有顯著不同，則各區(qū)域之間可以共用一套參數(shù)。該方法需要擬合無限制模型（Full model）和限制模型（Reduced model）。無限制模型對不同區(qū)域有不同的一組參數(shù)值，限制模型對所有區(qū)域有一組相同的參數(shù)值。在本文中，無限制模型即加入啞變量的模型，限制模型即基本模型（1）。采用F檢驗進行比較，F(xiàn)值計算如下：

式中：SSEF為無限制模型的誤差平方和；SSER為限制模型的誤差平方和；dfF為無限制模型的自由度；dfR為限制模型的自由度。

2.4 模型評價和檢驗指標

擬合結(jié)果采用平均絕對誤差（MAE）、均方根誤差（RMSE）、確定系數(shù)（R2）和赤池信息量（AIC）、貝葉斯信息量（BIC）來對比不同模型。檢驗結(jié)果采用平均絕對誤差（MAE）、平均預測誤差百分比（MPE）、平均絕對百分比誤差（MAPE）、均方根百分比誤差（RMSPE）進行檢驗。

3 結(jié)果與分析

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

選擇9個不同的分位點（τ=0.1, 0.2, 0.3, 0.4, 0.5,0.6, 0.7, 0.8, 0.9），利用SAS統(tǒng)計軟件的非線性規(guī)劃PROC NLP模塊對分位數(shù)模型（2）進行非線性分位數(shù)回歸擬合。利用SAS軟件中的PROC NLIN過程對啞變量模型（4）和基本模型（1）進行擬合，得到各方法模型的參數(shù)估計值，結(jié)果見表2。由表2可以看出，基于不同分位數(shù)的參數(shù)估計各不相同，整體有一定的趨勢走向。隨著分位點τ的增大，漸近線參數(shù)a逐漸增加，斜率參數(shù)b先增加后減小，形狀參數(shù)c逐漸變小，這反映出在不同分位數(shù)處的樹高曲線的變化規(guī)律?；趩∽兞糠椒ǖ母鲄^(qū)域參數(shù)也各不相同，且不同于基本模型求出的參數(shù)。

表3顯示了整個大興安嶺地區(qū)和4個劃分區(qū)域之間樹高-胸徑模型差異的測試結(jié)果。式（5）用來檢驗使用全部數(shù)據(jù)的限制模型是否足以滿足4個劃分區(qū)域的需要。從表3可以看出，使用全部數(shù)據(jù)的限制模型不足以描述4個劃分區(qū)域的樹高-胸徑關(guān)系（F=49.32，P=0）。然后利用4個區(qū)域兩兩進行檢驗：將4個區(qū)域的數(shù)據(jù)兩兩進行合并，分別構(gòu)建成對區(qū)域的含有啞變量的無限制模型（啞變量個數(shù)都為1），對成對區(qū)域的無限制模型和限制模型進行擬合，并利用式（5）計算F值及相應的P值。除區(qū)域2和區(qū)域4外，其他5個成對區(qū)域P值都顯著，表明其他5個成對區(qū)域有顯著不同；而區(qū)域2和區(qū)域4差異不顯著（F=1.27，P=0.284 4），說明區(qū)域2和區(qū)域4的樹木長勢沒有明顯區(qū)別，可以利用同一個模型進行擬合。

表2 基于不同方法的樹高-胸徑模型參數(shù)估計?Table 2 Parameter estimates of height-diameter models based on different methods

表3 不同區(qū)域模型對比F檢驗Table 3 F-test for models for different regions

為了更直觀地分析分位數(shù)模型模擬樹高曲線的趨勢，利用擬合數(shù)據(jù)繪制不同分位數(shù)的樹高曲線（圖2），其中A為基本模型，B為中位數(shù)模型（τ=0.5），C為3個分位數(shù)組合（τ=0.1,0.5,0.9），D為5個分位數(shù)組合（τ=0.3,0.4,0.5,0.6,0.7），E為7個分位數(shù)組合（τ=0.1,0.2,0.3,0.5,0.7,0.8,0.9），F(xiàn)為9個分位數(shù)組合（τ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9）。由圖2A—B可知，中位數(shù)曲線與基本模型曲線十分接近，說明擬合數(shù)據(jù)近似正態(tài)分布。而從圖2B—F可以看出，各分位數(shù)曲線都在樹高觀測值中穿過，使得觀測值插在兩個相鄰的分位數(shù)曲線中。無論是數(shù)據(jù)密集的小徑階，還是數(shù)據(jù)分布稀疏的大徑階，分位數(shù)回歸方法可以清楚地描述樹高曲線的變化規(guī)律，這也說明不同分位數(shù)曲線能夠表達不同區(qū)域不同立地條件的樹高曲線變化。

圖2 不同分位數(shù)興安落葉松樹高-胸徑曲線模擬Fig.2 Height-diameter models simulation of Larix gmelinii based on base model, median model and different quantile groups

3.2 模型擬合結(jié)果總體評價

一組QR模型可以用來模擬樹高-胸徑關(guān)系，并檢驗不同大小的樹的分位數(shù)之間的差異[9]。因此可以根據(jù)擬合統(tǒng)計量選擇出各區(qū)域所對應的最優(yōu)分位數(shù)模型，進而比較各區(qū)域樹木生長的差異。將表2中的參數(shù)估計值代入Richards方程中，計算各樹高-胸徑模型的擬合統(tǒng)計量，結(jié)果見表4。由表4可看出，各個區(qū)域都可以根據(jù)擬合統(tǒng)計量從其9個分位數(shù)模型中選出一個最優(yōu)的分位數(shù)模型，即在所有分位數(shù)模型中擁有最小的MAE、RMSE、AIC、BIC值和最大的R2值，4個區(qū)域所對應的最優(yōu)分位數(shù)模型分別為τ=0.7、τ=0.3、τ=0.5和τ=0.3。各區(qū)最優(yōu)分位數(shù)模型與各區(qū)啞變量模型擬合結(jié)果差異不大，除區(qū)域3外，這2種方法均優(yōu)于基本模型。區(qū)域3最優(yōu)分位數(shù)模型為中位數(shù)模型，即τ=0.5，且基本模型的擬合效果也較好，說明區(qū)域3的樹高曲線與總體擬合數(shù)據(jù)的樹高曲線趨勢比較接近，適用于代表平均水平的基本模型。區(qū)域3的基本模型與其最優(yōu)分位數(shù)模型（τ=0.5）的AIC、BIC值相差很?。éIC=3.116、ΔBIC=3.115），幾乎可以認為是同一模型。而區(qū)域1、2、4如果僅用基本模型進行估計誤差會較大。區(qū)域1的最優(yōu)分位數(shù)模型為τ=0.7，說明區(qū)域1樹木長勢分布偏高，立地條件較好。區(qū)域2、4的最優(yōu)分位數(shù)模型都為τ=0.3，說明這2個區(qū)域的樹木長勢分布相似，這與F檢驗時得到的結(jié)論一致，即區(qū)域2和區(qū)域4的樹木長勢沒有明顯區(qū)別，可以利用同一個模型進行擬合。且與其他2個區(qū)域相比，區(qū)域2、4都出現(xiàn)長勢偏低的趨勢，也說明這2個區(qū)域的立地條件相對較差。

表4 基于不同方法的樹高-胸徑模型擬合統(tǒng)計量?Table 4 Goodness-of-fit statistics of height-diameter models based on different methods

3.3 模型檢驗

為了更好地評價樹高-胸徑模型的預測能力，模型檢驗至關(guān)重要。利用檢驗數(shù)據(jù)對3種方法構(gòu)建的樹高-胸徑模型進行檢驗，結(jié)果見表5。就本文中構(gòu)建的3種模型來看，區(qū)域1、3、4的最優(yōu)分位數(shù)回歸模型都要優(yōu)于啞變量模型，且這2種模型都優(yōu)于基本模型。雖然區(qū)域2的啞變量模型要優(yōu)于最優(yōu)分位數(shù)模型，但區(qū)域2的啞變量模型沒有通過正態(tài)性檢驗（P=0.028 6），則區(qū)域2的最優(yōu)模型仍然為τ=0.3時的分位數(shù)模型。因此，4個區(qū)域的最優(yōu)樹高-胸徑模型如下：

區(qū)域1：

區(qū)域2和區(qū)域4：

區(qū)域3：

表5 不同方法的樹高-胸徑模型獨立性檢驗?Table 5 Validation statistics for height-diameter models based on different methods

3.4 不同樹高-胸徑模型應用模擬

為了更直觀地比較3個模型的預測效果，利用檢驗數(shù)據(jù)進行不同模型的樹高曲線模擬。圖3給出了各區(qū)域?qū)淖顑?yōu)分位數(shù)模型、啞變量模型、基本模型的興安落葉松樹高模擬曲線。由圖3可以看出，在4個區(qū)域中，最優(yōu)分位數(shù)模型與啞變量模型在各區(qū)域擬合曲線十分接近，趨勢也幾乎相同，都反映了數(shù)據(jù)的平均變化趨勢。基本模型在區(qū)域1、2和4偏差較大，而在區(qū)域3與最優(yōu)分位數(shù)模型和啞變量模型模擬曲線比較接近，這也說明區(qū)域3的數(shù)據(jù)分布近似于總體數(shù)據(jù)的分布。

4 結(jié)論與討論

圖3 不同方法的樹高-胸徑曲線預測模擬Fig.3 The simulation of height-diameter models based on different methods

準確的樹高預測對于估算樹木材積、評價立地指數(shù)、描述林分的發(fā)展和演替至關(guān)重要[28]。建立大尺度范圍的林木模型是林業(yè)建模中的一項重要工作，但隨著方程通用性的提高，用于局部地區(qū)的誤差也會增大[29]，因此迫切地需要了解不同小區(qū)域之間的差異。本研究按照地理方位，將大興安嶺地區(qū)劃分為西北（區(qū)域1）、東北（區(qū)域2）、中部（區(qū)域3）和南部（區(qū)域4）4個區(qū)域，引入啞變量建立大興安嶺地區(qū)4個區(qū)域的無限制模型，并與整體區(qū)域的限制模型進行F檢驗，發(fā)現(xiàn)限制模型不足以描述4個區(qū)域的樹高-胸徑關(guān)系。對兩兩區(qū)域之間的差異進行比較，發(fā)現(xiàn)東北（區(qū)域2）和南部（區(qū)域4）之間沒有明顯差別，可以共用同一模型。其他5對區(qū)域都有顯著不同，說明大興安嶺東北部和南部的樹木長勢十分相似，當?shù)氐沫h(huán)境條件，如氣候、土壤和水源等可能非常相似，需要進一步調(diào)查研究。

樹高-胸徑模型擬合時最常使用的方法就是最小二乘法。賴巧玲等[30]發(fā)現(xiàn)最小二乘法對異常數(shù)據(jù)比較敏感，因而當樣本量不夠大時，最小二乘法得到的參數(shù)估計不夠穩(wěn)定。Calama等[2]也發(fā)現(xiàn)，當數(shù)據(jù)有相關(guān)性時使用普通最小二乘回歸法會導致置信區(qū)間的偏差。與最小二乘法相比，分位數(shù)方法具有很多優(yōu)勢[31]：適合異方差模型，不需要對模型干擾項做假定，不易受到異常值的影響，能更全面地刻畫條件分布的特征等。Zang等[9]構(gòu)建了落葉松人工林的樹高-胸徑基本模型和分位數(shù)模型，但并沒有對兩者進行比較。本研究利用擬合數(shù)據(jù)進行9個分位點（τ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6,0.7, 0.8, 0.9）的分位數(shù)模型的擬合，發(fā)現(xiàn)無論是數(shù)據(jù)密集與否，分位數(shù)回歸都可以清楚地描述樹高曲線的變化規(guī)律。4個區(qū)域各自的最優(yōu)分位數(shù)模型和啞變量模型沒有顯著差異，可以根據(jù)實際情況進行選擇：若數(shù)據(jù)量足夠大，滿足獨立、正態(tài)、等方差的假設，這兩種方法都可以使用；若數(shù)據(jù)不滿足基本假設，或數(shù)據(jù)量較少時，建議使用分位數(shù)模型。盡管區(qū)域2的擬合統(tǒng)計量顯示啞變量模型優(yōu)于最優(yōu)分位數(shù)模型，但區(qū)域2數(shù)據(jù)不符合正態(tài)分布，因此認為區(qū)域2的最優(yōu)模型為分位數(shù)模型。此外，根據(jù)各區(qū)域得到的最優(yōu)分位數(shù)模型的分位點的不同，也可以得到不同區(qū)域間樹木長勢的信息。區(qū)域1的最優(yōu)分位數(shù)模型所對應的分位數(shù)是τ=0.7，則相對于其他3個區(qū)域，當樹木具有相同的胸徑時，區(qū)域1所對應的樹高偏高，說明區(qū)域1的立地條件、水源、光照比較適合樹木的生長。區(qū)域2和區(qū)域4的最優(yōu)分位數(shù)模型所對應的分位數(shù)都是τ=0.3，說明2個區(qū)域的樹木長勢相似，環(huán)境條件可能十分相似，這一點可以從區(qū)域2和區(qū)域4的限制模型和非限制模型的比較中得到相同的結(jié)論。區(qū)域3的最優(yōu)分位數(shù)模型所對應的分位數(shù)是τ=0.5，且擬合效果與基本模型相似，說明區(qū)域3的數(shù)據(jù)分布與總體數(shù)據(jù)的分布十分相似，可以利用基本模型直接擬合。

本研究采用林業(yè)上廣泛應用的Richards生長方程構(gòu)建了大興安嶺不同區(qū)域的興安落葉松樹高-胸徑模型，對比了分位數(shù)模型、啞變量模型和基本模型的擬合效果和預測能力。結(jié)果表明：各個區(qū)域都有其對應的最優(yōu)分位數(shù)模型，且各區(qū)域最優(yōu)分位數(shù)模型與啞變量模型擬合結(jié)果差異不大，并且都優(yōu)于基本模型。從預測能力來看，最優(yōu)分位數(shù)模型要略優(yōu)于啞變量模型，且都優(yōu)于基本模型。利用啞變量模型和分位數(shù)模型可以了解不同區(qū)域之間的差異，為森林資源管理者和規(guī)劃者提供更準確的樹木生長發(fā)育信息。在進行方法選擇時，可以根據(jù)數(shù)據(jù)特征進行選擇。本研究所使用的是伐倒木實測數(shù)據(jù)，由于樣本數(shù)量有限，只比較了樹高-胸徑模型的平均預測，隨著數(shù)據(jù)的積累，可以將分位數(shù)模型、組合分位數(shù)模型和混合效應模型進行比較。