基于分位數(shù)回歸的杉木人工林地位級劃分方法研究

張博，陳科屹，周來，Sajjad Saeed，張雅馨，孫玉軍＊

(1.北京林業(yè)大學(xué)森林資源和環(huán)境管理國家林業(yè)和草原局重點實驗室，北京 100083；2.中國林業(yè)科學(xué)研究院林業(yè)科技信息研究所，

北京 100091；3.中國林業(yè)科學(xué)研究院森林生態(tài)環(huán)境與保護研究所，北京 100091)

森林立地質(zhì)量是林地更新、樹種選擇、地力維持、生產(chǎn)力評估和經(jīng)營管理等林業(yè)工作和研究的基礎(chǔ)[1-2]。作為評價林分生長類型和林地生產(chǎn)力的重要依據(jù)，立地質(zhì)量評價對研究森林生長收獲規(guī)律、預(yù)估林地生產(chǎn)力和制訂相應(yīng)營林措施具有重要的指導(dǎo)意義[3-5]。目前，利用林分生長量的數(shù)據(jù)來定量評價立地質(zhì)量的方法包括：地位級法、地位指數(shù)法和立地形法等[1,6]。3種方法分別根據(jù)林分平均樹高與林分平均年齡、優(yōu)勢木平均高與標準年齡以及優(yōu)勢木平均高與基準胸徑的關(guān)系來反映立地條件對樹高生長的影響。在森林經(jīng)理經(jīng)營調(diào)查中，由于林分平均高和林分平均年齡是必備調(diào)查因子，地位級法被國內(nèi)外學(xué)者廣泛關(guān)注[7-9]。

傳統(tǒng)的地位級表編制方法基于生長模型、誤差調(diào)整和簡單的圖形解釋雖然能客觀反映總體的立地水平及其差異，但采用均值回歸模型和標準差調(diào)整法編制地位級表對建模數(shù)據(jù)本身不具有描述性，林分地位級的估測結(jié)果可能出現(xiàn)不同程度的高估或低估現(xiàn)象，或者現(xiàn)實林地的生長狀況可能根本達不到相應(yīng)的林分條件平均高[10]。因此，本研究提出一種基于分位數(shù)回歸模型（Quantile Regression Model）構(gòu)建地位級表的方法。分位數(shù)回歸是一種估計因變量或特定分位數(shù)函數(shù)的完全條件分布的方法，由Koenker等人[11]提出，是統(tǒng)計學(xué)和計量經(jīng)濟學(xué)常用的回歸分析方法之一[12]。一般的均值回歸估計方法僅針對協(xié)變量的條件均值或中心效應(yīng)[13]，而分位數(shù)回歸方法可靈活地描述相應(yīng)條件分位數(shù)自變量和因變量之間關(guān)系并得到的回歸曲線。目前，分位數(shù)回歸模型已被應(yīng)用于多個林業(yè)研究中，如：森林資源量的估算[14]，林分密度[15]，直徑分布預(yù)測[16]，直徑增長[17]，削度[18]和森林的地上生物量[19]等。

根據(jù)分位數(shù)模型可得出各分位數(shù)對應(yīng)的條件估計值和不易受到極端值影響的特點，該模型可用于描述樹高生長過程與立地質(zhì)量的關(guān)系[20-21]，但目前該模型在立地質(zhì)量評價研究中鮮有報道。基于此，本研究以福建省三明市將樂國有林場的杉木（Cunninghamia lanceolata(Lamb.)Hook.）純林為例，構(gòu)建基于分位數(shù)回歸模型的地位級表，并與標準差調(diào)整法進行比較，對將樂縣國有林場立地質(zhì)量進行評價與分析，為進一步提高地位級分級策略的效率和立地質(zhì)量評價的準確性提供理論依據(jù)和參考。

1 數(shù)據(jù)與研究區(qū)概況

數(shù)據(jù)來源于福建省三明市將樂國有林場2012 至2017年期間調(diào)查的418 塊杉木純林小班調(diào)查樣地數(shù)據(jù)，主要分布于將樂縣南口鄉(xiāng)、完全鄉(xiāng)、黃潭鎮(zhèn)、白蓮鎮(zhèn)、水南鎮(zhèn)、余坊鄉(xiāng)、光明鄉(xiāng)和古鏞鎮(zhèn)。樣地均為面積0.06 hm2的方形樣地。主要調(diào)查內(nèi)容包括各樣地地理坐標、坡度、坡向、坡位、海拔等地形因子，土壤類型、土壤厚度、腐殖質(zhì)層厚度等土壤因子，每木檢尺測定樣地內(nèi)每株樹木的胸徑、樹高、冠幅等，通過計算得出各樣地的林分平均年齡、林分平均胸徑、林分平均樹高等林分因子。樣地各齡級的信息統(tǒng)計見表1。

表1 各齡級樣地基本信息Table 1 Summary of basic information statistics of sample plots for each age-class

2 研究方法

2.1 基礎(chǔ)模型選擇

基于已有的研究成果選擇了7個常用的生長模型（表2）作為基礎(chǔ)模型[22-25]。

為了比較模型的擬合優(yōu)度，用赤池信息準則（AIC）、貝葉斯信息準則（BIC）、對數(shù)似然值（logLik）、均方根誤差（RMSE）、決定系數(shù)（R2）和平均絕對誤差（MAE），選擇擬合模型，所有計算均使用R（Version 3.6.1）軟件進行。

2.2 地位級表的編制（標準差調(diào)整法）

以導(dǎo)向曲線為基礎(chǔ)，按標準年齡時樹高值和地位級距（C），采用標準差調(diào)整法，可形成地位曲線簇（即樹高生長曲線簇）[1]。杉木在25 a 左右樹高生長區(qū)域穩(wěn)定，且杉木在20～30 a 時達到數(shù)量、經(jīng)濟成熟齡[26]。因此，本研究以林分樹高生長量趨于穩(wěn)定、杉木達到成熟齡確定基準年齡（A0）為25 a。

（1）擬合各齡級樹高標準差方程

根據(jù)各齡級樹高標準差（SH）與齡級平均年齡（Ai），利用SH=a+b×lg（Ai）式擬合齡級樹高標準差方程。將各齡級代入，計算出各齡級樹高標準差理論值（SA）。本研究方程為：

（2）導(dǎo)算地位級表

通常在基準年齡（A0）時，由于導(dǎo)向曲線的理論樹高值可能不是地位級數(shù)值，因此需要根據(jù)基準年齡時的樹高（H0）與標準差理論值（S0）的大小進行調(diào)整，公式如下：

式中，Hij為第i 齡級第j 地位級調(diào)整后的樹高；Hik為第i 齡級的導(dǎo)向曲線樹高；H0j為基準年齡時第j 地位級的樹高；H0k為基準年齡時導(dǎo)向曲線樹高；為基準年齡所在齡級樹高標準差理論值；為第i 齡級樹高標準差理論值。

以調(diào)整后的導(dǎo)向曲線為準，按地位級距C 逐齡級導(dǎo)算出各地位級曲線上的樹高值，其余地位級的調(diào)整系數(shù)Kj 為：

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

分位數(shù)回歸模型基于表2 中的線性和非線性模型來預(yù)測第τ 分位數(shù)的樹高模型：

與最小二乘方法相比，分位數(shù)回歸模型的參數(shù)通過最小化分位回歸領(lǐng)域的損失函數(shù)（或稱為檢驗函數(shù)）獲得[27]。

在本研究中，分位數(shù)回歸模型的構(gòu)建通過使用R 語言中的“quantreg”包來完成[28]。

2.4 地位級分級和評價

計算每個樣地林分平均高與各分級曲線樹高預(yù)測值的差值平方和（或差值的絕對值），以確定每個樣地的地位級，具有最小殘差平方和的模型即為該樣地所屬的立地類型，從而確定該林分目前的地位級。分別采用傳統(tǒng)方法與分位數(shù)回歸模型方法統(tǒng)計出每個林分的地位級，對分級結(jié)果進行比較和差異分析。

3 結(jié)果與分析

3.1 導(dǎo)向曲線的擬合

根據(jù)數(shù)據(jù)資料整理，采用最小二乘法和非線性擬合技術(shù)擬合7個基礎(chǔ)模型，其擬合結(jié)果如表3。依據(jù)AIC、BIC、RMSE 和MAE 最小，logLik 和R2值最大的原則，選出最優(yōu)模型Mod.4（Logistic）為最優(yōu)導(dǎo)向曲線方程：

表3 導(dǎo)向曲線模型擬合結(jié)果匯總Table3 Fitting statistics for Guide curve growth models

3.2 地位級表

（1）基準年齡及地位級距

根據(jù)基準年齡確定條件，本研究中杉木的標準年齡為25 a。杉木達到基準年齡時基準樹高H0為14.24 m，樹高變動范圍為6.8～19.3 m。根據(jù)將樂地區(qū)杉木的編表資料，以及樹高、胸徑的絕對變動幅度和經(jīng)營水平，確定地位級距C 為2 m，即地位級H0j分別為6 至20 的8個地位級。

（2）地位級表的編制

標準年齡A0（25 a）代入導(dǎo)向曲線方程得到樹高理論值H0k（14.24 m）并計算調(diào)整系數(shù)Kj和各相應(yīng)齡級樹高值繪制地位級表（表4）。根據(jù)立地級表和導(dǎo)向曲線方程擬合8個地位級的生長趨勢模型（圖1），模型參數(shù)及統(tǒng)計量如表5。

表4 杉木人工林地位級及相應(yīng)樹高Table 4 Site class and corresponding tree height of Chinese fir

表5 地位級模型參數(shù)及統(tǒng)計量Table 5 Parameter estimation and fitting statistics for site class models

圖1 傳統(tǒng)方法的地位級分布和建模數(shù)據(jù)散點分布Fig.1 Distribution of site classes by traditional method and scatter plot of modeling data

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

地位等級通常分成5～7 級，為便于與傳統(tǒng)地位級表編制方法比較，本研究分位數(shù)回歸模型的地位級表也分為8個地位級。根據(jù)數(shù)據(jù)的分布和導(dǎo)向曲線方程（Logistic 模型）選擇了8個分位數(shù)（0.01、0.05、0.15、0.30、0.70、0.85、0.95、0.99）（圖2）。其中，分位數(shù)0.01 和0.99 接近于地位級的下限和上限，可作為地位級Ⅷ和Ⅰ，分位數(shù)0.05 和0.95是研究區(qū)地位級的最前5%和最后5%水平，分別記作地位級Ⅶ和Ⅱ，分位數(shù)0.15、0.30、0.70、0.85則依據(jù)數(shù)據(jù)的分布狀況以及保持分位數(shù)回歸曲線簇形狀的相對均勻，分別記作地位級Ⅵ、Ⅴ、Ⅳ和Ⅲ，模型參數(shù)統(tǒng)計結(jié)果（表6）。

表6 分位數(shù)回歸模型參數(shù)及統(tǒng)計量Table 6 Parameter estimation and fitting statistics for quantile regression models of site class

圖2 分位數(shù)回歸方法的地位級分布和建模數(shù)據(jù)散點分布Fig.2 Distribution of site classes by quantile regression method and scatter plot of modeling data

3.4 兩種分級方法結(jié)果比較

從兩種方法的地位級分級過程看，傳統(tǒng)方法地位級分級曲線規(guī)則，能夠反映出林分地位級的普遍規(guī)律，得到的地位級表示基準年齡時的林分平均高。而分位數(shù)回歸的8個分位數(shù)點的選擇是根據(jù)建模數(shù)據(jù)分布來確定的，地位級分級結(jié)果反映出的是不同特定條件下的變化規(guī)律。因此，擬合的分級曲線可能存在過于接近的現(xiàn)象，基準年齡較小時無法直接判斷地位級分級狀況，造成對地位級分級結(jié)果模糊或?qū)ζ毡橐?guī)律的理解不足的現(xiàn)象。

利用各模型的最小殘差平方和來判斷418 塊樣地的地位級狀況，從兩種方法的分級結(jié)果（表7）可以看出，傳統(tǒng)地位級表中大多數(shù)樣地地位級分布于10 至16 地位級，與分位數(shù)回歸模型中多數(shù)樣地分級分布于Ⅵ至Ⅲ地位級結(jié)果基本一致。對兩種分級結(jié)果進行差異顯著性檢驗，分位數(shù)回歸模型的地位級分級效果與傳統(tǒng)方法沒有顯著差異（F值0.130 3，P值0.719 3），表明分位數(shù)回歸模型能夠較為準確的描述研究區(qū)內(nèi)林分的地位級的分布特征。兩種方法分級的主要差異在于對最大和最小地位級的分級表現(xiàn)，分位數(shù)回歸方法將更多的樣地劃分到最大和最小地位級中，反映了研究區(qū)各齡級林分平均高的完整分布狀況。總體而言，傳統(tǒng)方法的擬合曲線表現(xiàn)出的是平均狀態(tài)和預(yù)設(shè)的變化規(guī)律，而分位數(shù)回歸模型利用數(shù)據(jù)分布特征所反映出的是特定狀態(tài)和條件對應(yīng)的變化規(guī)律，其劃分的地位級結(jié)果在很大程度上取決于建模樣本的結(jié)構(gòu)和數(shù)據(jù)質(zhì)量。

表7 傳統(tǒng)方法（左）和分位數(shù)回歸模型（右）對樣地地位級分級的結(jié)果Table 7 The result of site classes grouped by traditional method (left) and quantile regression model (right)

4 結(jié)論

本研究提出一種基于分位數(shù)回歸模型的立地質(zhì)量分級和評價方法，以福建省三明市將樂縣國有林場小班調(diào)查樣地數(shù)據(jù)構(gòu)建基于分位數(shù)回歸模型的地位級分級模型，并與傳統(tǒng)地位級劃分方法相比較。結(jié)果表明，基于不同分位數(shù)擬合曲線簇，可準確評估林地地位級水平；根據(jù)林分平均高與各分級曲線樹高預(yù)測值的差值平方和（或差值絕對值）最小的原則，可迅速確定該林分的地位級和生長類型，實現(xiàn)快速、準確的立地質(zhì)量評價，為精準評估森林生產(chǎn)力和進一步提升森林質(zhì)量評價效率提供方法和依據(jù)。

我國森林資源清查數(shù)據(jù)中主要記錄平均木的樹高，采用林分平均樹高構(gòu)建地位級表可以使分位數(shù)回歸方法與森林清查數(shù)據(jù)兼容。資源清查的小班調(diào)查數(shù)據(jù)雖然沒有標準地調(diào)查數(shù)據(jù)準確，但其數(shù)據(jù)量大且覆蓋廣，能夠從一定程度上反映立地質(zhì)量對林分生產(chǎn)力的影響[7]。分位數(shù)回歸模型可描述、分類、預(yù)測和驗證小班數(shù)據(jù)林分生長與地位級之間的關(guān)聯(lián)性，基于導(dǎo)向生長模型的分位數(shù)回歸曲線簇能夠更直觀的反映出不同地位級下杉木樹高的變化軌跡，從整體上描述了各小班林地生產(chǎn)力，并有效地簡化了地位級劃分的測算工作。因此，結(jié)合分位數(shù)回歸方法和森林資源連續(xù)清查數(shù)據(jù)構(gòu)建林分尺度和區(qū)域尺度的地位級分級體系，可以實現(xiàn)對數(shù)據(jù)的充分利用。但是，由于分位數(shù)回歸模型受數(shù)據(jù)分布的影響，地位級的劃分結(jié)果很大程度上取決于建模樣本的結(jié)構(gòu)和數(shù)據(jù)質(zhì)量。今后的研究中需要考慮地位級劃分和分位點選擇的聯(lián)系，以進一步提高分位數(shù)回歸模型在生產(chǎn)力評價中的適用性。