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

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

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

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

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

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

本研究以杉木大徑材經(jīng)營林分為對象，構(gòu)建基于不同立地質(zhì)量評價指標(biāo)的杉木樹高-胸徑混合效應(yīng)模型，通過模型評價來探討不同立地質(zhì)量評價指標(biāo)對杉木大徑材樹高-胸徑模型的影響，以期為杉木大徑材林分的結(jié)構(gòu)調(diào)控和經(jīng)營措施制定提供理論依據(jù)。

1 研究區(qū)概況

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

2 數(shù)據(jù)來源與處理

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

表1 林分調(diào)查因子統(tǒng)計Table 1 Statistics of forest investigation factors

3 研究方法

3.1 樹高-胸徑基礎(chǔ)模型

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

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

3.2 立地質(zhì)量評價指標(biāo)的計算

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

式中：S為不同的立地質(zhì)量評價指標(biāo)，包括優(yōu)勢高、地位級、地位指數(shù)和立地形；a、a0、b、c為模型參數(shù)。

3.2.1 優(yōu)勢高

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

3.2.2 地位級

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

3.2.3 地位指數(shù)

地位指數(shù)（Site index，SI）根據(jù)各不同林分在基準(zhǔn)年齡時的優(yōu)勢木樹高平均值作為立地質(zhì)量高低的依據(jù)，是一種直觀的定量指標(biāo)。本研究中各林分地位指數(shù)的計算采用前期根據(jù)福建省杉木人工林標(biāo)準(zhǔn)地構(gòu)建的地位指數(shù)模型，導(dǎo)向曲線經(jīng)測試選用Lundqvist-Korf方程：

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

式中：A0為基準(zhǔn)年齡（25 a）。

3.2.4 立地形

立地形（Site form，SF）是根據(jù)基準(zhǔn)胸徑下的優(yōu)勢木平均樹高作為評價立地質(zhì)量的依據(jù)。本研究對立地形的計算采用前期根據(jù)福建省杉木人工林標(biāo)準(zhǔn)地構(gòu)建的杉木立地形模型，其導(dǎo)向曲線經(jīng)測試選用Logistic方程：

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

式中：D0為基準(zhǔn)胸徑（25 cm）。

3.3 混合效應(yīng)模型

在固定效應(yīng)模型的基礎(chǔ)上，添加樣地水平的隨機(jī)效應(yīng)，構(gòu)建包含樣地效應(yīng)的混合模型，其一般形式為：

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

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

3.4 模型評價與檢驗

3.4.1 模型檢驗

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

3.4.2 模型評價指標(biāo)

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

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

4 結(jié)果與分析

4.1 杉木樹高-胸徑固定效應(yīng)模型構(gòu)建

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

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

4.2 杉木樹高-胸徑混合效應(yīng)模型構(gòu)建

4.2.1 基礎(chǔ)混合模型

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

表3 基礎(chǔ)混合效應(yīng)模型擬合Table 3 Fitting results of mixed effect base model

4.2.2 含不同立地質(zhì)量評價指標(biāo)的混合模型

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

表4 含不同立地質(zhì)量評價指標(biāo)的混合效應(yīng)模型擬合?Table 4 Mixed effects model fitting with different site quality evaluation indexes

4.3 模型檢驗

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

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

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

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

表6 最優(yōu)混合效應(yīng)模型的參數(shù)估計Table 6 Parameter estimation of the optimal mixed effect model

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

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

5 討 論

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

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

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

本研究主要探究了立地條件對大徑材杉木樹高-胸徑關(guān)系的影響，但由于林分結(jié)構(gòu)和生長動態(tài)的復(fù)雜性和多變性，杉木的樹高-胸徑生長關(guān)系還可能受到經(jīng)營措施、競爭和氣候變化等因素的影響，從而導(dǎo)致建立的模型適用范圍受限。因此，未來研究還需進(jìn)一步增加對其他影響因素的研究，并擴(kuò)大研究范圍，以增強(qiáng)該模型的通用性，為杉木人工林可持續(xù)經(jīng)營提供參考依據(jù)。

6 結(jié) 論

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