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

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

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

北京 100091；3.中國林業科學研究院森林生態環境與保護研究所，北京 100091)

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

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

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

1 數據與研究區概況

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

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

2 研究方法

2.1 基礎模型選擇

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

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

2.2 地位級表的編制（標準差調整法）

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

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

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

（2）導算地位級表

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

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

以調整后的導向曲線為準，按地位級距C 逐齡級導算出各地位級曲線上的樹高值，其余地位級的調整系數Kj 為：

2.3 分位數回歸模型

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

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

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

2.4 地位級分級和評價

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

3 結果與分析

3.1 導向曲線的擬合

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

表3 導向曲線模型擬合結果匯總Table3 Fitting statistics for Guide curve growth models

3.2 地位級表

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

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

（2）地位級表的編制

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

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

表5 地位級模型參數及統計量Table 5 Parameter estimation and fitting statistics for site class models

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

3.3 分位數回歸模型

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

表6 分位數回歸模型參數及統計量Table 6 Parameter estimation and fitting statistics for quantile regression models of site class

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

3.4 兩種分級方法結果比較

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

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

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

4 結論

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

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