基于不同立地質量評價指標的杉木大徑材林分樹高-胸徑模型

夏洪濤，郭曉斌，張 珍，田相林，郭福濤，孫帥超

（1.福建農林大學 林學院，福建 福州 350002；2.國家林業和草原局杉木工程技術研究中心，福建 福州 350002；3.福建省上杭白砂國有林場，福建 龍巖 364205；4.西北農林科技大學 林學院，陜西 楊凌 712100）

杉木Cunninghamialanceolata是我國主要的速生用材樹種之一，生長快、產量高，且材質優良。根據第九次全國森林資源清查結果，杉木人工林在我國的面積和蓄積均位居我國人工喬木林首位，分別達990萬hm2和7.55億m3[1]。杉木平均輪伐期一般在25～30 a，因此市場上的杉木多為中小徑材。近年來，隨著我國經濟的發展，市場對于高品質的大徑材杉木需求越來越大[2]。培育和經營大徑材杉木不僅能通過高價值木材獲取更多經濟效益，而且能大大增加杉木人工林的碳匯功能，有效提高其生態效益[3]。同時，輪伐期的延長還可以有效減緩傳統杉木人工林集約經營引起的地力衰退問題。因此，大徑材經營是當前我國杉木人工林經營的重要方向。

樹高和胸徑是林分結構特征中最重要的測樹因子，也是計算林分蓄積和生物量等關鍵指標的基本變量[4-7]。樹木胸徑可以快速、準確地測定，但樹高的測量成本相對較高，因此實際調查中往往按徑階測量部分樹高，并通過構建樹高-胸徑模型來估算缺失的樹高[8-11]。同時，樹高-胸徑模型也是鏈接不同維度林分結構的工具，是森林生長模型系統中的基本組分[12]，建立科學準確的樹高-胸徑模型是實現林分生長預估和結構合理調控的必要前提。

樹高和胸徑的關系除了取決于樹種特性外，還可能受到立地條件[13]、競爭[14-15]和營林措施[16]等因素的影響。其中，立地條件能夠反映出相同直徑樹木的不同樹高生長潛力，因此對樹高-胸徑關系的影響最為顯著。近年來，大量研究發現在樹高-胸徑模型中加入立地變量能有效提高模型的預測精度[17-19]。然而，關于模型中立地質量評價指標的選擇卻各有不同。目前，在同齡林研究中多采用地位指數作為立地質量指標[18]，而在異齡林中則多采用優勢高[19]。Molina-Valero等[20]在輻射松Pinus radiata人工林中對地位指數和立地形進行比較研究，發現二者在反映立地條件方面具有相同的效果。此外，地位級由于簡便易用且效果穩定，也是我國長期林業生長實踐中常用的立地質量評價指標。不同立地質量評價指標的構建原理不同，各有優點，但目前關于不同立地質量評價指標在樹高-胸徑模型中的效果評價還較為缺乏。

目前，學者們已開發出多種不同模型來描述樹高-胸徑關系[21-22]。但在利用森林調查數據建模時，應用傳統的最小二乘回歸技術經常違背隨機和獨立觀測的假設，并且沒有考慮到自相關的存在。為解決這一問題，混合效應模型被廣泛應用于樹高-胸徑模型的構建[23-25]。相較于固定效應模型，混合效應模型既能體現總體的平均變化趨勢，又能反映個體間的差異[26-29]。大量研究表明，混合效應模型在保證模型形式簡潔性的同時，也可顯著提高樹高-胸徑模型的預測精度[30-33]。

本研究以杉木大徑材經營林分為對象，構建基于不同立地質量評價指標的杉木樹高-胸徑混合效應模型，通過模型評價來探討不同立地質量評價指標對杉木大徑材樹高-胸徑模型的影響，以期為杉木大徑材林分的結構調控和經營措施制定提供理論依據。

1 研究區概況

研究地位于福建省龍巖市西南部的上杭白砂國有林場（116°30′～116°38′E，25°04′～25°15′N），屬亞熱帶氣候，雨量充沛。年平均氣溫20.1 ℃，年降水量約1 600 mm。該地區森林資源豐富，森林覆蓋率為96.6%，主要代表樹種有杉木、馬尾松Pinusmassoniana、柳杉Cryptomeria fortunei和福建柏Fokieniahodginsii等。

2 數據來源與處理

上杭白砂國有林場現存全國集中連片最大的杉木大徑材林分約164 hm2[34]，2021年7月采用分層系統抽樣方法設置了杉木大徑材固定樣地48塊，樣地規格為25.82 m×25.82 m。每塊樣地進行每木檢尺后選取5株最高的杉木為優勢木，計算優勢木平均樹高作為該樣地的優勢高。本研究共得到1 762株杉木的樹高-胸徑數據，其中優勢木240株。將所有數據按1∶1隨機劃分為2組，采用交叉驗證方法進行建模和檢驗。表1為各因子的統計結果。

表1 林分調查因子統計Table 1 Statistics of forest investigation factors

3 研究方法

3.1 樹高-胸徑基礎模型

用于描述樹高-胸徑關系的生長方程眾多，為了篩選出適合杉木大徑材的樹高-胸徑基礎模型，本研究利用全部樣地數據對Richards[35]、Mitcherlich[36]、Korf[37]、Vanclay[38]、Logistic[39]、Weibull[40]和Gompertz[41]7個生長方程進行擬合比較。其中，Richards方程、Logistic方程和Weibull方程均取得了最好的擬合效果（決定系數為0.840～0.841），由于Richards方程具有合理解析性和預測穩定性等優點，在生產實踐中具有廣泛的適用性[42-43]，故選取Richards方程作為基礎模型。其模型形式如下：

式中：H為林木樹高；D為林木胸徑；a為樹木生長的最大值參數，b為生長速率參數，c為同化作用冪指數有關的參數。

3.2 立地質量評價指標的計算

本研究選取優勢高、地位級、地位指數和立地形4個描述立地質量的指標，分別將其添加到杉木樹高-胸徑基礎模型中，評價不同立地質量評價指標對模型的影響。根據前人研究[25]，立地條件主要影響林分內的樹高最大生長潛力，故將立地質量評價指標添加在樹高最大值參數a上，模型形式如下：

式中：S為不同的立地質量評價指標，包括優勢高、地位級、地位指數和立地形；a、a0、b、c為模型參數。

3.2.1 優勢高

優勢高（Dominant height，HT）即優勢木的平均樹高，在很多林分中可直接作為評價立地質量的指標[44]。本研究計算每塊樣地中最高的5株杉木樹高平均值，將其作為該林分優勢高。

3.2.2 地位級

地位級（Site class，SC）是根據林分平均樹高和林分年齡來評定立地質量的指標，一般分為若干等級，地位級越高說明立地條件越好。本研究采用張博等[45]基于分位數回歸建立的福建杉木人工林地位級表來確定各林分對應的地位級。

3.2.3 地位指數

地位指數（Site index，SI）根據各不同林分在基準年齡時的優勢木樹高平均值作為立地質量高低的依據，是一種直觀的定量指標。本研究中各林分地位指數的計算采用前期根據福建省杉木人工林標準地構建的地位指數模型，導向曲線經測試選用Lundqvist-Korf方程：

式中：HT為林分優勢木平均高；A為林齡。由此得到杉木人工林地位指數方程：

式中：A0為基準年齡（25 a）。

3.2.4 立地形

立地形（Site form，SF）是根據基準胸徑下的優勢木平均樹高作為評價立地質量的依據。本研究對立地形的計算采用前期根據福建省杉木人工林標準地構建的杉木立地形模型，其導向曲線經測試選用Logistic方程：

式中：DT為優勢木的平均胸徑。由此得到如下立地形方程：

式中：D0為基準胸徑（25 cm）。

3.3 混合效應模型

在固定效應模型的基礎上，添加樣地水平的隨機效應，構建包含樣地效應的混合模型，其一般形式為：

式中：Hij和Dij分別為第i個樣地中第j株杉木的樹高和胸徑值；f(β,μi,Dij)為樹高與胸徑間的函數關系，β為固定效應參數向量，μi為隨機參數向量，服從期望為 0、方差為ψi的正態分布；εij為誤差項，服從期望為0、方差為Rij的正態分布；ψi和Rij分別為樣地間和樣地內的方差-協方差矩陣。對基礎模型（式1）中不同固定參數分別加入樣地效應，本研究對比測試了含1～3個隨機效應參數的模型形式。以固定參數a、b、c均添加隨機效應的混合模型為例，其一般表達式（式9），在此基礎上對包含不同立地質量評價指標的Richards方程（式2）中的固定參數添加樣地效應，得到其一般表達式（式10）：

式中：μ1、μ2、μ3和μ0為隨機效應參數；S為不同立地質量評價指標。本研究利用R軟件中的nlme包對不同的樹高-胸徑非線性混合效應模型進行擬合。

3.4 模型評價與檢驗

3.4.1 模型檢驗

為充分利用數據信息，本研究采取二折的交叉檢驗方法對各模型進行檢驗。其基本原理是將全部數據隨機分為2組，首先利用組1的樣本建模并用組2的樣本檢驗，然后利用組2的樣本建模并用組1的樣本檢驗，最后利用所有樣本的實測值與預測值之間的誤差來計算各評價指標。

3.4.2 模型評價指標

采用赤池信息準則（AIC）、貝葉斯信息準則（BIC）和-2倍對數似然值（-2LogL）作為模型擬合評價指標，3個指標值越小，說明模型擬合效果越好。采用擬合優度（R2）、平均絕對誤差（MAE）和均方根誤差（RMSE）作為模型檢驗指標，R2越大，MAE和RMSE越小，說明模型精度越高、預測效果越好。各指標計算公式如下[46]：

式中：L為估計模型的似然函數最大值；p為參數個數；n為樣本數；和Hi分別為第i株杉木的樹高預測值和觀測值。

4 結果與分析

4.1 杉木樹高-胸徑固定效應模型構建

通過在基礎模型中加入不同立地質量評價指標，構建出包含不同立地質量評價指標的杉木樹高-胸徑固定效應模型，表2為各模型的擬合結果。與基礎模型（M1）相比，增加不同立地質量評價指標后的模型擬合效果均有明顯提高。其中加入優勢高的固定效應模型（M2）擬合效果最好(AIC=7 228.198，BIC=7 255.387，-2LogL=7 218.198)，其次分別為加入地位指數（M4）和立地形（M5）的樹高-胸徑模型，而加入地位級的模型（M3）表現相對較差。

表2 含不同立地質量評價指標的樹高-胸徑模型擬合Table 2 Fitting results of height-diameter models with different site quality evaluation indexes

4.2 杉木樹高-胸徑混合效應模型構建

4.2.1 基礎混合模型

表3是在基礎模型（公式1）中加入不同參數組合樣地效應的擬合結果。含隨機效應的混合模型AIC、BIC和-2LogL值均小于基礎模型，說明其擬合結果均優于固定效應模型。當僅包含1個隨機參數時，混合模型中含μ1的擬合效果較好。當混合模型含2個隨機參數時，混合模型中含μ1、μ3的擬合效果最優。當基礎模型中包含3個隨機參數μ1、μ2、μ3時，AIC、BIC和-2LogL值均有所增大。因此，包含隨機參數μ1、μ3的杉木樹高-胸徑混合效應模型是最優的混合效應基礎模型。

表3 基礎混合效應模型擬合Table 3 Fitting results of mixed effect base model

4.2.2 含不同立地質量評價指標的混合模型

在公式2中加入不同參數組合的樣地隨機效應，分別篩選出包含不同立地質量評價指標的最優混合效應模型。只考慮1個隨機參數時，加入優勢高的混合模型中包含隨機參數μ3的擬合效果較優，而加入立地形、地位級和地位指數的混合模型均為包含隨機參數μ1的模型擬合效果較優?？紤]2個隨機參數時，加入不同立地質量評價指標的混合模型均為包含隨機參數μ1、μ3的模型擬合效果最好?？紤]3個或4個隨機參數時，模型擬合效果不能再提高。表4列出了包含不同立地質量評價指標的最優混合效應模型，其中以優勢高為立地質量評價指標且包含隨機參數μ1、μ3的混合模型（M7）取得了最好的擬合優度（AIC=7 121.858、BIC=7 159.923、-2LogL=7 107.858）。

表4 含不同立地質量評價指標的混合效應模型擬合?Table 4 Mixed effects model fitting with different site quality evaluation indexes

4.3 模型檢驗

表5為10個模型的檢驗結果。與基礎模型（M1）相比，加入不同立地質量評價指標或隨機效應后均能顯著提高模型的預測效果，混合效應模型（M6～M10）能比固定效應模型（M1～M5）取得更好的預測效果。在固定效應和混合效應模型中，加入優勢高的杉木樹高-胸徑模型都分別取得了最好的預測效果（M2和M7），其次為加入地位指數的杉木樹高-胸徑模型（M4和M9）。在所有模型中，以優勢高為立地質量評價指標的杉木樹高-胸徑混合效應模型（M7）的預測效果最優。

表5 不同模型檢驗結果Table 5 Validation results of different models

以4 cm為徑階距，對比杉木大徑材林分樹高-胸徑基礎模型（M1）和最優混合模型（M7）的不同徑階殘差（圖1）。從圖1可以看出，模型M7在不同徑階的殘差中值都在0軸附近，而基礎模型（M1）的殘差中值則會偏離0軸上下波動，說明加入優勢高的混合效應模型的殘差對稱性比基礎模型更好，預測結果更為穩定可靠。在12以上的徑階中，模型M7的殘差最大值、最小值和箱體均明顯小于基礎模型，說明其預測誤差的變異程度更小。表6為采用全部數據擬合得到的模型M7（最優混合模型）的參數估計值與隨機效應方差-協方差。

圖1 基礎模型與最優混合效應模型各徑階殘差對比Fig.1 Residual comparison of each diameter class between the base model and the optimal mixed effect model

表6 最優混合效應模型的參數估計Table 6 Parameter estimation of the optimal mixed effect model

圖2是利用建立的最優非線性混合效應模型對杉木大徑材林分的樹高預測結果。從圖2可以看出，除4徑階和48以上徑階由于樹木很少造成其樹高變異較小外，其余各徑階的樹高變異都很大，說明加入立地變量和樣地效應對杉木樹高的預測有較大影響。在徑階較大時，杉木樹高的整體變異性更高，此時立地條件對樹高生長的影響更顯著。

圖2 最優模型各徑階樹高預測Fig.2 Tree height prediction of each diameter class of the optimal model

5 討 論

準確預測樹高對于估算單木材積、評價立地質量、描述林分生長至關重要。目前，關于杉木樹高-胸徑模型的研究已有不少，但大多都是基于短輪伐期的中小徑材杉木林分[47-48]。本研究中的杉木大徑材林分平均年齡約40 a，此時杉木樹高生長趨于穩定，變化緩慢。因此，相較于其他杉木樹高-胸徑模型，本研究所構建的模型更適合于長輪伐期經營的杉木林分。同時，由于大徑材林分中杉木的樹高生長潛力已經充分發揮，不同立地條件所帶來的樹高生長差異也更為穩定，因而也更適用于不同立地質量評價指標的評價。

本研究通過比較7個常用的樹高-胸徑生長方程，確定Richards方程為杉木樹高-胸徑基礎模型。該模型的3個參數決定著樹高的上漸近值和曲線形狀變化，反映了樹高隨著胸徑生長逐漸增大后變緩的趨勢，使模型具有較好的預測功能和生物學意義。加入立地條件變量后，顯著提高了杉木樹高-胸徑模型的預測效果，這與一些學者的研究結果相一致[13,17-19]。地位級和地位指數主要用于同齡林，而優勢高和立地形由于不依賴林分年齡，可同時用于同齡林和異齡林[18-19]。過去對杉木人工林生長模型的研究主要以地位指數作為立地質量評價指標[47,50]，但本研究的結果表明，以優勢高為立地質量評價指標的杉木樹高-胸徑模型的擬合和預測效果均最好（MAE=1.596，RMSE=2.034，R2=0.871）。原因可能在于，優勢高是根據調查結果計算出的直接指標，而其他3個指標都是根據數據擬合加工后的間接指標，本研究中選取的杉木大徑材林分的樹高生長趨于穩定，直接使用優勢高可以避免其他指標折算過程中產生的誤差，因此能更好地反映立地質量。同時，本研究中立地形與地位指數的預測效果極為接近（表5），這與Molina-Valero等[20]的研究結果一致，說明特定情況下立地形可以作為地位指數的替代指標。由于兼具簡便易用和適用范圍廣的優點，優勢高和立地形在其他生長模型中的效果值得進一步研究。

與固定效應模型相比，混合效應模型能進一步提高杉木樹高-胸徑模型的精度，而且加入樣地效應后各立地質量評價指標的綜合表現排序并未發生變化（表5）。過去的大量研究已經證明混合效應模型在樹高-胸徑模擬中的優越性[26-28,31]，但對于隨機效應添加的參數卻有所不同。例如王冬至等[30]將樣地效應添加于截距參數上，臧顥等[51]在每個固定參數上都添加隨機效應，?z?elik等[25]通過比較發現不同樹種適合的隨機效應添加方式有所不同。本研究通過比較隨機效應添加在不同參數的擬合效果，均發現將隨機效應作用于Richards方程的μ1和μ3時模型擬合最優，這與陳浩等[28]建立的馬尾松樹高-胸徑模型的結果相一致。

本研究主要探究了立地條件對大徑材杉木樹高-胸徑關系的影響，但由于林分結構和生長動態的復雜性和多變性，杉木的樹高-胸徑生長關系還可能受到經營措施、競爭和氣候變化等因素的影響，從而導致建立的模型適用范圍受限。因此，未來研究還需進一步增加對其他影響因素的研究，并擴大研究范圍，以增強該模型的通用性，為杉木人工林可持續經營提供參考依據。

6 結 論

本研究以大徑材經營的杉木林分為對象，通過引入不同的立地質量評價指標，構建出杉木非線性混合效應樹高-胸徑模型，并對不同立地質量評價指標效果進行評價。研究結果表明，在杉木樹高-胸徑模型中加入立地質量評價指標和樣地隨機效應均能顯著提高模型擬合及預測精度，其中以優勢高為立地質量評價指標的非線性混合效應模型的表現效果最好。在4個不同立地質量評價指標中，優勢高的綜合表現最優，其次為地位指數和立地形，因此在一定情況下，優勢高和立地形可以替代地位指數用于反映杉木人工林立地質量。大徑材經營是當前我國杉木人工林經營的重要方向，本研究可為杉木大徑材林分的生長預測和林分結構調控提供科學依據。