黃 靚，盧 侃，易 烜，朱晉梅，王琢玙，劉文劍

（1.湖南省青羊湖國有林場，湖南 寧鄉(xiāng) 410627； 2.湖南省林業(yè)局，湖南 長沙 410007）

樹高和胸徑是森林資源調(diào)查過程中最重要的兩個變量和常用的測量指標(biāo)，它們在描述林分結(jié)構(gòu)、評價林地生產(chǎn)潛力以及估算林分蓄積量、生物量的過程中，發(fā)揮了至關(guān)重要的作用［1］。胸徑通常可用圍尺等測量工具直接測量，而樹高的測量在較為郁閉的林分中受視覺影響較大，導(dǎo)致測量誤差較大［2］。盡管采用高精度的測量儀器（如激光、超聲波測高器等）能減小樹高測量誤差，但對于大區(qū)域、大面積測量來說仍然費時費力，測量成本過高［3］。因此，估計林分或林木高度的最常用方法是測量其子樣本，以子樣本為基礎(chǔ)，構(gòu)建一個樹高-胸徑模型來估計它們的高度。

傳統(tǒng)的樹高-胸徑模型構(gòu)建往往采用普通最小二乘法，但最小二乘法通常要考慮樣本的獨立性假設(shè)，即獨立、等方差、正態(tài)分布［4-5］。然而，林分生長數(shù)據(jù)一般具有層次結(jié)構(gòu)、重復(fù)測量等特點，表現(xiàn)為數(shù)據(jù)之間存在異方差和自相關(guān)，不滿足普通最小二乘法擬合回歸模型的獨立性假設(shè)［5］。同時，林木的生長通常受到樹種本身特性以及復(fù)雜的生境條件影響，如區(qū)域、樣地、林分條件、立地條件以及干擾等，很多學(xué)者在構(gòu)建樹高-胸徑模型時，通過在模型中引入這些除胸徑以外的變量來提高樹高-胸徑模型的精度與適用性［1，6-9］?；旌闲?yīng)模型作為處理多層次重復(fù)測量數(shù)據(jù)的有力工具，已廣泛應(yīng)用于樹高-胸徑模型等林業(yè)生長模型的構(gòu)建，它既可以基于固定效應(yīng)參數(shù)預(yù)估樹高的平均變化趨勢，又能夠以隨機效應(yīng)參數(shù)來反映不同區(qū)域、樣地以及立地條件間的樹高曲線變化。Crecente-Campo等［10］以西班牙北部藍桉（Eucalyptus globulus）林分為研究對象，構(gòu)建了廣義的樹高-胸徑混合效應(yīng)模型，發(fā)現(xiàn)這些模型解釋了83%的觀測變異，且平均誤差小于2.5 m。樊偉等［8］考慮到樣地水平的隨機效應(yīng)，構(gòu)建了大別山地區(qū)杉木（Cunninghamia lanceolata）樹高-胸徑混合效應(yīng)模型，發(fā)現(xiàn)混合效應(yīng)模型擬合效果優(yōu)于基礎(chǔ)模型以及含林分變量的模型。Ciceu等［1］基于雷泰扎特國家公園178個固定樣地歐洲云杉（Picea abies）異齡林?jǐn)?shù)據(jù)，考慮了描述林分結(jié)構(gòu)、物種混交和競爭的23個林分變量，利用其中3個林分預(yù)測因子構(gòu)建了一個廣義的樹高-胸徑混合效應(yīng)模型，能夠較為精確地預(yù)測該地區(qū)挪威云杉的樹高。

馬尾松（Pinusmassoniana）是我國的主要用材樹種之一，其分布范圍廣泛，其資源在我國森林資源總量中占有重要比重。根據(jù)第九次全國森林資源清查結(jié)果，馬尾松林面積和蓄積分別占我國喬木林面積和蓄積的4.47%和3.67%，在水源涵養(yǎng)、碳增匯以及用材方面有著至關(guān)重要的地位［11］。為精確估算湘西土家族苗族自治州馬尾松林樹高，本研究以現(xiàn)有調(diào)查的湘西土家族苗族自治州8個縣（市）馬尾松林調(diào)查數(shù)據(jù)為基礎(chǔ)，擬構(gòu)建含區(qū)域、齡組以及立地條件中顯著影響因子的樹高-胸徑混合效應(yīng)模型，為找到適用于湘西土家族苗族自治州馬尾松林樹高-胸徑模型提供理論依據(jù)。

1 研究區(qū)概況

研究區(qū)位于湘西土家族苗族自治州，地理坐標(biāo)為109°10＇-110°23＇E、27°44＇-29°38＇N，地處云貴高原東側(cè)的武陵山區(qū)，地勢東南低、西北高，海拔97.1～1 736 m。屬于亞熱帶季風(fēng)濕潤氣候，冬季比較寒冷，夏季溫高濕重，春夏之交陰濕多雨，年平均溫度在16.0～17.0℃之間，年平均降水量為1 284.2～1 416.9mm，降雨主要集中在每年4-9月。主要土壤有黃紅壤、山地黃壤、山地黃棕壤、草甸土、石灰土、酸性紫色土、堿性紫色土，全州大部分土壤呈微酸性至中性，過酸或過堿性土壤不多。森林覆蓋率為70.24%，森林資源豐富，但林種較為單一，山地植被的豐富度不夠，以杉、松為主的針葉純林面積占喬木林面積的56%，生態(tài)功能較強的闊葉林、混交林、復(fù)層林只占32%。

2 材料與方法

2.1 數(shù)據(jù)來源

數(shù)據(jù)來源于湘西土家族苗族自治州2019年7-8月調(diào)查的馬尾松人工林?jǐn)?shù)據(jù)。數(shù)據(jù)采集于湘西土家族苗族自治州的保靖縣、鳳凰縣、古丈縣、花垣縣、吉首市、龍山縣、瀘溪縣和永順縣8個縣（市），共77個樣地，樣地大小為25 m×25 m。在樣地中選取3株標(biāo)準(zhǔn)木，用激光測高測距儀測量樹高，用圍尺測量胸徑，最后取各樣地3株標(biāo)準(zhǔn)木樹高、胸徑的平均值分別作為林分平均高、林分平均胸徑。采用三倍標(biāo)準(zhǔn)差法剔除異常數(shù)據(jù)后為71株馬尾松林分平均木。除了林分平均高和平均胸徑外，同時還調(diào)查并記錄了各樣地馬尾松林的齡組、海拔、坡度、坡向、坡位及林分密度等。由于所調(diào)查的樣地林分密度為156 0～163 0株·hm-2，較為接近，對研究區(qū)樣地馬尾松生長影響不大，故未考慮林分密度對馬尾松的影響。各縣（市）樣本的樹高和胸徑數(shù)據(jù)見表1。

表1 樣木基本情況Tab.1 Basic information of sam ple wood

2.2 研究方法

2.2.1 立地因子等級劃分

參照《中國森林立地》劃分標(biāo)準(zhǔn)，對調(diào)查的海拔、坡度、坡位、坡向等4個立地因子進行分級［12］。其中，為了區(qū)分不同海拔梯度的影響，將海拔每200 m劃分為一級，具體劃分標(biāo)準(zhǔn)見表2。

表2 立地因子等級劃分Tab.2 Classification of site factors

2.2.2 顯著性影響因子選取

基于R 4.2.3軟件，采用方差分析、多重比較Scheffe法［13］并采用箱線圖將結(jié)果可視化，對湘西8個不同區(qū)域、4個不同齡組（這里將成熟林與過熟林歸為一組）、不同立地因子（包括不同海拔、坡度、坡向以及坡位）間的樹高是否存在顯著差異進行分析，從而篩選出對湘西土家族苗族自治州馬尾松樹高生長影響顯著的因子。

2.2.3 基礎(chǔ)模型的篩選

選取6組基本非線性模型（3個2參數(shù)模型和3個3參數(shù)模型）作為基礎(chǔ)模型（表3），比較它們的模型精度，選取最優(yōu)基礎(chǔ)模型進而構(gòu)建馬尾松樹高-胸徑混合效應(yīng)模型。

表3 候選基礎(chǔ)模型Tab.3 Candidate basic models

2.2.4 混合效應(yīng)模型構(gòu)建

大量研究表明樹高-胸徑模型的構(gòu)建過程中，樹高除了與胸徑呈明顯數(shù)量關(guān)系外，其他變量也會影響樹高曲線的變化，包括立地條件、林分因子、競爭以及經(jīng)營措施等［20-22］。本研究在最優(yōu)基礎(chǔ)模型的基礎(chǔ)上，考慮到區(qū)域、齡組以及4個不同立地因子對樹高曲線的影響，選擇其中的顯著性影響因子，擬構(gòu)建包含這些顯著性影響因子及其組合隨機效應(yīng)的馬尾松樹高-胸徑混合效應(yīng)模型，進而從中擇出最優(yōu)模型。非線性混合效應(yīng)模型（nonlinear mixed-effectmodels，NLMEMs）是依據(jù)回歸函數(shù)依賴于效應(yīng)參數(shù)非線性關(guān)系而建立的數(shù)學(xué)模型，其中，效應(yīng)參數(shù)包括固定效應(yīng)參數(shù)和隨機效應(yīng)參數(shù)，這里假設(shè)區(qū)域因子顯著，以其作為隨機效應(yīng)為例，其一般形式如下［23］：

式中：hi與Di分別為第i個區(qū)域的林分平均高和林分平均胸徑，εi為誤差向量，β、μi分別為固定效應(yīng)參數(shù)和隨機效應(yīng)參數(shù)。

混合效應(yīng)模型在參數(shù)效應(yīng)和隨機效應(yīng)的方差協(xié)方差結(jié)構(gòu)確定后，有兩個問題必須要考慮：一是模型的異方差，二是模型的自相關(guān)性［24］。由于本研究的數(shù)據(jù)不是重復(fù)觀察數(shù)據(jù)，所以不需要考慮自相關(guān)性。

2.2.5 模型評價與檢驗

使用R 4.2.3軟件的nls函數(shù)、nlme包的nlme函數(shù)分別對上述基礎(chǔ)模型以及混合效應(yīng)模型進行擬合，求解模型參數(shù)、計算模型精度評價與檢驗指標(biāo)。以赤池信息量準(zhǔn)則 （Akaike information criterion，AIC）及貝葉斯信息準(zhǔn)則 （Bayesian information criterion，BIC）選擇基礎(chǔ)模型 （即AIC、BIC值越小，模型擬合效果越好），同時對模型進行F檢驗，判斷模型回歸效果是否顯著。采用留一交叉驗證法［25-26］對模型進行檢驗，在此過程中，假設(shè)樣本量為n，將單個觀測值用于驗證集，其余n-1個觀測值作為訓(xùn)練集；然后，對排除的觀測值進行預(yù)測，并在n-1個觀測上重復(fù)n次此過程，產(chǎn)生n個測試誤差估計；基于n個測試誤差估計值計算確定系數(shù) （R2）、均方根誤差（root mean square error，RMSE）、總相對誤差（total relative error，TRE）。R2越接近于1，RMSE值、TRE的絕對值越小，模型精度越高。根據(jù)3個指標(biāo)計算結(jié)果比較基礎(chǔ)模型與樹高-胸徑混合效應(yīng)模型的擬合效果。各指標(biāo)計算公式如下：

式（2）～（6）中：CAI為赤池信息量準(zhǔn)則值，CBI為貝葉斯信息準(zhǔn)則值，R2為確定系數(shù)，ERMS為均方根誤差，ETR為總相對誤差，n為樣本數(shù)，l表示模型極大似然函數(shù)值，p為模型中參數(shù)的個數(shù)，hi為第i樣本樹高實測值，是第i樣本樹高預(yù)估值，為實測樹高的平均值。

3 結(jié)果與分析

3.1 影響因子篩選及樹高差異特征分析

由方差分析可知：區(qū)域、齡組以及4個立地因子（海拔、坡度、坡向、坡位）中，區(qū)域、齡組、坡向為研究區(qū)內(nèi)馬尾松樹高的顯著性影響因子（P＜0.05）（表4）。通過Scheffe法多重比較以及箱線圖可視化（圖1）可得：研究區(qū)8個不同縣（市）內(nèi)的馬尾松樹高生長差異明顯，其中古丈縣的馬尾松整體長勢最好，而花垣縣、龍山縣馬尾松的樹高總體偏低；從齡組來看，馬尾松樹高隨幼齡林、中齡林、近熟林以及成熟林、過熟林趨勢增加，其中馬尾松成熟林、過熟林明顯高于幼齡林、中齡林；不同坡向間馬尾松樹高生長差異顯著，研究區(qū)半陽坡馬尾松樹高整體高于其他三個坡向的；研究區(qū)內(nèi)不同海拔、坡度、坡位馬尾松樹高整體差異均不顯著。

圖1 不同區(qū)域、齡組以及立地條件馬尾松林分平均高Fig.1 Stand mean height of Pinusmassoniana in different regions，age groups and site conditions

表4 影響因子的顯著性檢驗Tab.4 Significance test of influence factors

3.2 基礎(chǔ)模型的擬合與篩選

使用R語言的nls函數(shù)對6個樹高-胸徑基礎(chǔ)模型進行擬合，得到模型擬合參數(shù)，并使用R2、AIC、BIC對基礎(chǔ)模型進行篩選，具體結(jié)果如表5所示。由表5可知，6個基礎(chǔ)模型的R2相差不大，為0.493 7～0.510 6，其中3參數(shù)模型M4的R2最大，但其3個參數(shù)a、b、c中僅參數(shù)a顯著（具有統(tǒng)計學(xué)意義）。而AIC、BIC值最小的模型為2參數(shù)模型 M2 （Schumacher）， R2＝0.509 7、AIC＝347.636 4、BIC＝354.424 4，同時其參數(shù)a、b均極顯著，具有統(tǒng)計學(xué)意義。綜合可得，選擇M2為最優(yōu)基礎(chǔ)模型。

3.3 混合效應(yīng)模型構(gòu)建

使用R語言中nlme包nlme函數(shù)擬合樹高-胸徑混合效應(yīng)模型。根據(jù)3.1的顯著性因子篩選結(jié)果，本研究考慮的隨機效應(yīng)因子為區(qū)域（AR）、齡組（AG）、坡向（AP）。在最優(yōu)基礎(chǔ)模型M2中加入這些隨機變量，擬合結(jié)果如表6所示。

表6 混合效應(yīng)模型參數(shù)估計及模型精度評價Tab.6 Parameter estimation and model accuracy evaluation ofm ixed-effectsmodels

根據(jù)R2、RMSE、TRE可以看出，混合效應(yīng)模型的擬合效果均優(yōu)于最優(yōu)基礎(chǔ)模型。只含一個隨機效應(yīng)因子時，M2.3擬合效果最好；當(dāng)包含兩個隨機效應(yīng)因子時，M2.12擬合效果最好；當(dāng)包含三個隨機效應(yīng)因子時，M2.21擬合效果最好。同時，可以看出M2.21的擬合效果優(yōu)于其他任意一個模型，其模型精度最高，R2＝0.845 3，與最優(yōu)基礎(chǔ)模型M2相比精度提升了65.84%；RMSE、TRE分別為1.730 0、 -0.060 7，RMSE 比 M2 的降低了39.01%，TRE絕對值比M2的降低了73.52%。

一般地，可以通過比較擬合模型的殘差分布圖很好地判斷誤差的異方差性。本研究繪制2個不同模型的殘差圖用以檢驗擬合模型是否存在異方差情況（圖2）。從圖2可以看出，最優(yōu)基礎(chǔ)模型M2、最優(yōu)混合效應(yīng)模型M2.21的殘差均表現(xiàn)為隨機分布的趨勢，不存在明顯異方差。同時可以發(fā)現(xiàn)，最優(yōu)混合效應(yīng)模型M2.21殘差分布更加收斂，擬合精度高于最優(yōu)基礎(chǔ)模型M2。

圖2 最優(yōu)基礎(chǔ)模型M 2與最優(yōu)混合效應(yīng)模型M 2.21的殘差圖Fig.2 Residual plot of optimal basic model M 2 and optimalm ixed-effectsmodel M 2.21

4 結(jié)論與討論

4.1 結(jié)論

本研究以湘西土家族苗族自治州8個縣（市）采集的馬尾松樣本數(shù)據(jù)為基礎(chǔ)，首先利用方差分析篩選出研究區(qū)馬尾松樹高生長顯著性影響因子，同時通過多重比較Scheffe法以及箱線圖可視化方法對不同區(qū)域、齡組以及立地條件的樹高差異特征分析；然后采用最小二乘法擬合6個候選基礎(chǔ)模型，選擇最優(yōu)基礎(chǔ)模型；最后構(gòu)建含顯著性影響因子及其組合的混合效應(yīng)模型，并比較最優(yōu)基礎(chǔ)模型與混合效應(yīng)模型以及混合效應(yīng)模型之間的擬合精度、擬合誤差、模型復(fù)雜程度，篩選出適用于湘西土家族苗族自治州的最優(yōu)樹高-胸徑模型。主要結(jié)論為：

（1）區(qū)域、齡組、坡向為研究區(qū)內(nèi)馬尾松樹高的顯著性影響因子。研究區(qū)8個不同縣（市）內(nèi)的馬尾松樹高生長差異明顯，其中古丈縣的馬尾松整體長勢最好，而花垣縣和龍山縣的馬尾松樹高總體偏低；從齡組來看，馬尾松樹高隨幼齡林、中齡林、近熟林以及成熟林、過熟林趨勢增加，其中馬尾松成熟林、過熟林明顯高于幼齡林、中齡林；不同坡向間馬尾松樹高生長差異顯著，研究區(qū)半陽坡馬尾松樹高整體高于其他三個坡向的；研究區(qū)內(nèi)不同海拔、坡度、坡位馬尾松樹高生長整體差異并不顯著。

（2）6組候選基礎(chǔ)模型中，M2（Schumacher）模型最優(yōu)，其AIC、BIC值最小，分別為347.636 4、354.424 4，R2＝0.509 7，同時其參數(shù)a、b均極顯著，具有統(tǒng)計學(xué)意義。

（3）在21個混合效應(yīng)模型中，最優(yōu)混合效應(yīng)模型為M2.21，即區(qū)域、齡組、坡向作為隨機效應(yīng)同時作用在最優(yōu)基礎(chǔ)模型M2的參數(shù)a、b時的混合效應(yīng)模型，其模型精度最高，R2＝0.845 3，與最優(yōu)基礎(chǔ)模型M2相比精度提升了65.84%；RMSE、TRE分別為1.730 0、-0.060 7，與M2相比，RMSE降低了39.01%，TRE的絕對值降低了73.52%。

4.2 討論

馬尾松林是湘西土家族苗族自治州主要用材林之一，其在生產(chǎn)經(jīng)濟用材、涵養(yǎng)水源、維持碳平衡等方面有著關(guān)鍵作用。建立準(zhǔn)確有效的馬尾松林樹高-胸徑模型對指導(dǎo)湘西土家族苗族自治州林業(yè)生產(chǎn)和實踐具有重要的現(xiàn)實意義。以往的研究表明，不同地區(qū)、不同林分類型以及不同生境條件下的林分生長規(guī)律不一。Zeide等［27］研究發(fā)現(xiàn)，林分密度是影響樹木的胸徑-樹高關(guān)系的主要因素。Meng等［28］研究表明，包含環(huán)境變量的樹高-胸徑模型能夠更好地預(yù)測不同區(qū)域內(nèi)的加州扭葉松（Pinus contorta var.latifolia）樹高。本研究顯示，區(qū)域、齡組以及坡向是影響湘西土家族苗族自治州馬尾松林樹高的主導(dǎo)因子，不同區(qū)域的綜合立地條件、微氣候、土壤理化性質(zhì)不一，從而導(dǎo)致馬尾松樹高生長具有差異；不同的齡組本身反映了馬尾松林生長的時間差；而不同坡向則影響了馬尾松林生長的光照時間與強度，從而使其生長產(chǎn)生差異。實際上，許多研究表明林分密度是影響樹木生長的關(guān)鍵因子之一，陳浩等［7］構(gòu)建了貴州省東北部地區(qū)馬尾松人工林樹高-胸徑非線性混合效應(yīng)模型，發(fā)現(xiàn)將密度、優(yōu)勢木平均高、胸高斷面積以不同個數(shù)及組合形式加入基礎(chǔ)模型，同時考慮樣地水平隨機效應(yīng)的混合效應(yīng)模型預(yù)測精度高、適用性強。本研究各樣地林分密度相近，尚未考慮林分密度對研究區(qū)內(nèi)馬尾松生長的影響，在后續(xù)的研究中，可補充調(diào)查不同林分密度的馬尾松人工林樣地，探究林分密度對湘西土家族苗族自治州馬尾松人工林生長的作用機制。

同時，不同基礎(chǔ)模型對不同地區(qū)、不同林分類型林分以及不同生境林分的適用性也不一樣。Xu等［6］在構(gòu)建中國東南地區(qū)杉木單木樹高-胸徑模型時發(fā)現(xiàn)，5個代表性基礎(chǔ)模型中，Korf模型最優(yōu)；而王君杰等［29］研究發(fā)現(xiàn)，12個基礎(chǔ)模型中，Richards方程更適用于大興安嶺的落葉松（Larix gmelinii）樹高-胸徑模型的構(gòu)建。本研究選取了6組非線性樹高-胸徑基礎(chǔ)模型，其中Schumacher模型（M2）最適用于湘西土家族苗族自治州馬尾松林樹高-胸徑模型的擬合。

混合效應(yīng)模型具有很好的兼容性，能夠?qū)?fù)雜的立地、林分以及干擾等因子在模型中表達出來，通過這種方法構(gòu)建的模型往往優(yōu)于傳統(tǒng)最小二乘法構(gòu)建的回歸模型，符利勇等［5］、段光爽等［30］以及Bronisz等［31］證明了這一點?？紤]到湘西土家族苗族自治州地處云貴高原東側(cè)的武陵山區(qū)，地形地貌復(fù)雜多樣，以及調(diào)查的各齡組、坡向馬尾松林樹高差異顯著，本研究構(gòu)建了包含區(qū)域、齡組、坡向隨機效應(yīng)的馬尾松樹高-胸徑混合效應(yīng)模型，結(jié)果顯示混合效應(yīng)模型擬合效果均優(yōu)于基礎(chǔ)模型。另外，本研究尚未采用分位數(shù)回歸、貝葉斯分層模型以及機器學(xué)習(xí)等方法來構(gòu)建湘西土家族苗族自治州馬尾松林樹高-胸徑模型，未來可比較這些方法擬合模型的優(yōu)劣性，進一步為湘西土家族苗族自治州馬尾松林林業(yè)生產(chǎn)實踐提供理論基礎(chǔ)。