基于立地混合效應(yīng)的湖南丘陵平原區(qū)杉木多形立地指數(shù)模型研究

宋 爭，黃 朗，胡 松，張雄清，朱光玉*

(1. 中國林業(yè)科學研究院林業(yè)科技信息研究所，北京 100091；2. 中南林業(yè)科技大學林學院，湖南 長沙 410004；3. 中國林業(yè)科學研究院林業(yè)研究所，北京 100091)

立地質(zhì)量評價是對森林立地的宜林性或潛在生產(chǎn)力判斷或預(yù)測的基礎(chǔ)工作[1]，立地指數(shù)是一種最普遍的評價森林立地質(zhì)量的方法[2-3]，分析某樹種林分基準年齡時對應(yīng)的林分優(yōu)勢高可用于分析森林立地生產(chǎn)潛力的基礎(chǔ)[4]。除了根據(jù)優(yōu)勢高與年齡指標，還能夠通過樣地實測值和收獲表蓄積量估算獲得立地指數(shù)[5]。但由于立地條件的差異性或多樣性，基于優(yōu)勢高-年齡的導(dǎo)向曲線模擬精度往往較低，同一樹種在不同區(qū)域的立地指數(shù)估計結(jié)果隨立地質(zhì)量的不同而存在顯著差異。因此，很多學者在評價立地質(zhì)量時，優(yōu)先分析環(huán)境要素（氣候、植被、地形等）與立地指數(shù)的關(guān)系，得到包含環(huán)境因子的立地指數(shù)函數(shù)關(guān)系[6-7]。也有學者采用樹高-年齡分級啞變量的方法，用樣地林木平均高劃分等級代替立地因子評價立地質(zhì)量[8]。這些方法雖然提高了立地指數(shù)建模精度，但關(guān)于環(huán)境因子在模型中的表達仍需進一步探索。

目前，如何建立立地指數(shù)模型成為描述森林立地質(zhì)量的重點和中心問題。研究指出有單形和多形兩種關(guān)系，其中多形立地指數(shù)曲線結(jié)合優(yōu)勢木平均高和其他反映立地條件因子，共同評定立地質(zhì)量[9-11]。最常用的2種擬合方法為固定基準年齡的靜態(tài)方程和可變基準年齡的動態(tài)方程，前者依賴于選定的基準年齡，后者局限于嚴格的數(shù)據(jù)要求[12]。而混合效應(yīng)模型是一種最常用的處理局部變異性的建模方法[13-14]。對于不同立地類型的優(yōu)勢高生長，混合模型不僅能夠描述整體變化的平均趨勢，還能夠給出協(xié)方差和方差等相關(guān)數(shù)據(jù)以了解個體特點[15-16]。雖然在構(gòu)建立地指數(shù)模型時，相關(guān)文獻運用數(shù)量化方法I[17-18]，然而該方法前提是固定效應(yīng)，如果考慮隨機效應(yīng)，該方法不適用。

鑒于此，本研究以湖南丘陵平原區(qū)杉木（Cunninghamia lanceolata(Lamb.) Hook.）人工林為研究對象，首先利用數(shù)量化方法I確定影響杉木優(yōu)勢高生長的顯著性影響因子；然后以主導(dǎo)立地因子及其組合作隨機效應(yīng)，建立基于立地隨機效應(yīng)的杉木多形立地指數(shù)曲線模型，通過模型選優(yōu)對不同立地因子水平的立地質(zhì)量進行評價；最后以Kmeans聚類劃分立地類型組，解決復(fù)雜立地類型的模型應(yīng)用問題，以提高立地指數(shù)曲線模型的適用性和準確性。

1 材料與方法

1.1 研究區(qū)概況

湖南位于中亞熱帶，地貌大部分為丘陵，呈現(xiàn)大陸性季風濕潤氣候，全年平均氣溫17.5 ℃左右，平均降水量1200～1700 mm，日照時長1300～1800 h，雨水豐富且熱量充沛。研究區(qū)喬木樹種以杉木（Cunninghamia lanceolata(Lamb.)Hook）為主，其次為馬尾松（Pinus massonianaLamb.）與擬赤楊（Alniphyllum fortunei(Hemsl.)Makino）等。主要小喬木與灌木有冬青（Ilex chinensisSims）、木姜子（Litsea euosmaW.W.Smith）、山茶（Camellia japonicaL.）、杜莖山（Maesa japonica(Thunb.) Moritzi.）等。主要草本植物有狗脊蕨（Woodwardia japonica(L. f.) Sm.）、鐵芒萁（Dicranopteris linearis(Burm.) Underw.）等。

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

以湖南丘陵平原區(qū)杉木人工純林為研究對象，在湘鄉(xiāng)、懷化、會同、株洲、永州、臨湘、綏寧、桃源等26個縣，設(shè)置樣地（20 m × 30 m）共360塊，記錄相關(guān)立地因子，并對樣地內(nèi)活立木每木檢尺。樣地內(nèi)每0.01 hm2選取最粗優(yōu)勢木1株，計算優(yōu)勢木平均高(H)。相關(guān)林分調(diào)查因子統(tǒng)計詳見表1。

表1 杉木林分調(diào)查因子統(tǒng)計Table 1 Statistics of stand factors for Cunninghamia lanceolata

1.3 研究方法

1.3.1 立地分類方法 選取如表2所示的影響杉木優(yōu)勢高生長的6個立地因子，根據(jù)《中國森林立地》相關(guān)標準進行分級，將海拔每100 m劃分為一級，分析不同海拔的差異。

表2 立地因子等級劃分Table 2 The division of site factor grades

立地因子作自變量，各林分優(yōu)勢木平均高做因變量，運用數(shù)量化方法Ⅰ對杉木立地因子分別分析?；诜讲罘治鲋胁煌癙r>F”大小，篩選顯著因子并進一步得到主導(dǎo)因子。

1.3.2 非線性混合效應(yīng)模型 為更準確地評價杉木立地質(zhì)量，本研究考慮立地因子作隨機效應(yīng)對杉木立地指數(shù)模型的影響，隨機效應(yīng)包括了立地因子與其交互作用。非線性混合效應(yīng)模型的形式為[19-20]：

式中：Hij為 第i個立地類型下第j個樣地的平均優(yōu)勢高，Tij為第i個 立地類型下第j個 樣地的林齡，φij為參數(shù)向量，εij是誤差項，M是立地類型數(shù)量，ni為第i個立地類型的樣地數(shù)，β為（p×1）維固定效應(yīng)向量，bi為（q×1） 維隨機效應(yīng)向量，Aij,Bij為設(shè)計矩陣，D為隨機效應(yīng)的方差—協(xié)方差矩陣，Ri為立地類型內(nèi)的方差—協(xié)方差結(jié)構(gòu)，σ2為方差，Gi為描述方差異質(zhì)性的對角矩陣，Γi為描述隨機效應(yīng)自相關(guān)性的方差矩陣。

構(gòu)建非線性混合效應(yīng)模型，需要對固定與隨機效應(yīng)參數(shù)進行構(gòu)造。通常情況下對所有可能的參數(shù)考慮隨機效應(yīng)，但這會造成模型不收斂[21-22]。為避免模型不收斂問題，本研究選擇在1個參數(shù)上考慮立地的隨機效應(yīng)。其次，由于本研究中林分優(yōu)勢高與年齡不是重復(fù)調(diào)查數(shù)據(jù)，不存在樣地內(nèi)測量值之間的相關(guān)性。因此，假設(shè)誤差項εij服從正態(tài)分布。

最后，本模型只考慮單個隨機效應(yīng)參數(shù)來確定隨機效應(yīng)參數(shù)方差，因此假定隨機參數(shù)構(gòu)造類型為無結(jié)構(gòu)，并采用AIC、BIC、對數(shù)似然比（Loglikelihood）和確定系數(shù)R2選擇最優(yōu)混合效應(yīng)模型。AIC、BIC值越小，Log-likelihood和R2越大，模型的擬合效果越好。

1.4 候選基礎(chǔ)模型

基礎(chǔ)模型不僅要能夠很好描述林木優(yōu)勢高生長規(guī)律，還要具有生物學的可解釋性[23-24]。其中立地指數(shù)曲線基礎(chǔ)模型常為“S”形，本研究選擇常用的立地指數(shù)曲線方程（M1～M8），擬合杉木林分優(yōu)勢高（見表3）。

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

1.5 模型精度評價

利用確定系數(shù)（R2） 、預(yù)估精度（P）、平均絕對誤差（MAE）和均方根誤差（RMSE）進行模型評價。

式中：yi為 第i個 樣地林分優(yōu)勢高實測值，為第i個樣地林分優(yōu)勢高估計值，n為樣本數(shù)，tα為置信水平α=0.05時的t分布值，p為模型參數(shù)個數(shù)，為平均預(yù)估值。

2 結(jié)果與分析

2.1 因子顯著性分析

林分優(yōu)勢木平均高受多種立地因子綜合作用影響。根據(jù)優(yōu)勢木平均高與不同立地因子的相互關(guān)系，采用數(shù)量化方法Ⅰ對其進行顯著性檢驗，結(jié)果（表4）表明：對優(yōu)勢木平均高具有顯著影響的立地因子有海拔、坡度、坡向與土壤類型，其顯著性順序為土壤類型＞海拔＞坡向＞坡度。

表4 立地因子的顯著性檢驗Table 4 Significance testing of site factors

2.2 基礎(chǔ)模型擬合與選擇

利用 Forstat2.1軟件對候選基礎(chǔ)模型M1～M8進行估計，相關(guān)基礎(chǔ)模型參數(shù)估計值、確定系數(shù)（R2） 、預(yù)估精度（P）、平均絕對誤差（MAE）及均方根誤差（RMSE）如表5所示。

表5 候選基礎(chǔ)模型擬合結(jié)果Table 5 The fitting results of candidate basic models

從8個候選基礎(chǔ)模型（表5）擬合結(jié)果可知，這8個模型的擬合精度均較低（R2=0.4243～0.5644），也證明了區(qū)域性立地指數(shù)建模精度低的問題。其中，模型M4擬合效果最好，其確定系數(shù)（R2= 0.5644）與預(yù)估精度（P=0.9750）最大，平均絕對誤差（MAE=6.7626）及均方根誤差（RMSE=2.4100）最小。

2.3 非線性混合效應(yīng)模型模擬

以主導(dǎo)立地因子作隨機效應(yīng)，考慮立地單因素與多因素交互作用，構(gòu)建立地指數(shù)混合模型。結(jié)合基礎(chǔ)模型參數(shù)個數(shù)和隨機效應(yīng)因子組合類型，擬合所有參數(shù)及其組合形式，最優(yōu)隨機效應(yīng)立地指數(shù)模型通過AIC和BIC等評價指標獲?。ū?）。

計算結(jié)果顯示，除HB、PX、TL等單因子在其它參數(shù)組合的擬合效果較優(yōu)外，受水平數(shù)影響，其它隨機效應(yīng)組合均只在漸進參數(shù)a上時，模型收斂。因此，含立地隨機效應(yīng)的杉木多形立地指數(shù)模型表達式為：

式中：a，b，c為模型固定效應(yīng)參數(shù)；μi為隨機效應(yīng)參數(shù)；H為林分平均優(yōu)勢木高；Age為林分年齡。

擬合多種隨機效應(yīng)因子組合類型后的結(jié)果如表6所示，AIC、BIC相對于基礎(chǔ)模型M4都顯著下降，Log-likelihood和R2均大幅提高，這說明基于立地隨機效應(yīng)的混合模型明顯優(yōu)于基礎(chǔ)模型。此外，混合模型擬合效果的好壞，取決于立地因子顯著性的高低，立地因子越顯著，相應(yīng)隨機效應(yīng)組合類型的混合模型擬合結(jié)果越好。其中，海拔、坡度、坡向、土壤類型4個主導(dǎo)立地因子交互作用的混合模型AIC、BIC最低，Log-likelihood及R2最高。根據(jù)立地分類的原則與方法，4個主導(dǎo)立地因子排列組合可得到立地類型，綜上所述含立地類型的混合模型為最優(yōu)模型。

表6 基于立地效應(yīng)的混合模型擬合結(jié)果Table 6 The fitting results of mixed effects models considering site effect

2.4 模型應(yīng)用

依據(jù)海拔、坡度、坡向、土壤類型4個主導(dǎo)因子，將研究區(qū)劃分為78個立地類型，那么相應(yīng)的混合模型隨機效應(yīng)參數(shù)有78個，勢必會影響立地指數(shù)模型的應(yīng)用。為了簡化立地類型，將初始立地類型（ST）進行K-means聚類。參考森林立地分類的標準，確定系數(shù)≥0.99為聚類精度標準，合并得分值相近的立地類型為立地類型組（GST），使得本研究區(qū)的立地類型聚類為11種（表7）。

表7 立地類型聚類結(jié)果Table 7 Classification of cluster with site types

將聚類后的立地類型組（GST）作為隨機效應(yīng)構(gòu)建立地指數(shù)非線性混合模型，利用AIC、BIC、Log-likelihood及R2等4個指標進行模型選擇（表8），并與基礎(chǔ)模型（None）、初始立地類型（ST）模擬結(jié)果進行比較。如表8所示，考慮立地的隨機效應(yīng)后，模型確定系數(shù)從0.5644提高到0.8089～0.8117，其中基于立地類型組的混合效應(yīng)模型模擬精度最高（R2=0.8117，Loglikelihood=-704.9879）。因此，基于立地隨機效應(yīng)的立地指數(shù)曲線模型可以極大提高模型精度，同時利用K-means聚類劃分立地類型組的方法，不僅保證了預(yù)估精度，而且還簡化了立地類型，提高立地指數(shù)模型的實用性。

表8 模型精度評價Table 8 Model evaluation

利用該模型繪制湖南丘陵平原區(qū)杉木多形立地指數(shù)曲線（圖1）。圖中展示了不同立地類型組的立地指數(shù)曲線圖，各立地類型組的模型模擬效果均較好。

圖1 多形立地指數(shù)曲線Fig. 1 Polymorphic site index curve

3 討論

林分優(yōu)勢高生長受林分密度影響很小，是同齡林立地質(zhì)量評價采用最廣泛的指標，但關(guān)于林分優(yōu)勢高的定義及測定標準并不一致[25]。趙美麗等[26]對樟子松（Pinus elliottiiEngelm.）同齡林最高、最粗木平均優(yōu)勢高的相關(guān)性、穩(wěn)定性研究后認為，最粗6株木法能科學得到林分優(yōu)勢高，該方法具有簡易、穩(wěn)定的優(yōu)點，所以本研究以每100 m2選擇1株最粗木的方法計算平均優(yōu)勢高。

在森林立地分類與評價的因子中，地形屬間接因素，但綜合反映著氣候和土壤特性。本研究在進行顯著性研究中用數(shù)量化方法I，以分析不同立地因子對杉木優(yōu)勢高生長的影響，結(jié)果表明，海拔、土壤類型、坡向、坡度與杉木優(yōu)勢高生長有強相關(guān)性。不同立地因子對優(yōu)勢高生長的影響也不同。產(chǎn)生顯著性差異的原因有以下幾點：從海拔來說，南方丘陵平原地區(qū)熱量和降水量的垂直分異總的來說決定了研究對象的生長與分布，海拔為一個主要因素。從土壤類型與坡度來說，杉木對土壤濕度、養(yǎng)分、通氣性等要求較高，杉木生長由土壤類型的差異直接影響，坡度間接表現(xiàn)土壤特性。從坡向來說，不同坡向熱量光照差異，對杉木不足以引起其生長的很大差異，故坡向相比較而言，顯著性影響要低于土壤類型、海拔，但高于坡度的影響作用。王冬至等[27]在研究塞罕壩華北落葉松（Larix principisrupprechtiiMayr）人工林地位指數(shù)模型時，將林分優(yōu)勢高與不同立地因子進行相關(guān)性分析，發(fā)現(xiàn)林分優(yōu)勢高與立地因子顯著相關(guān)。

基于林分優(yōu)勢高-年齡關(guān)系的立地指數(shù)，已成為表征森林立地生產(chǎn)力的一個重要參數(shù)，應(yīng)用該方法的前提在于導(dǎo)向曲線具有較高的模擬精度。然而對于某一特定的樹種，導(dǎo)向曲線的形狀會因氣候、土壤、地形等的不同而有所變化，使用同一導(dǎo)向曲線模擬樹高生長必然會造成較大誤差[28]。而混合模型在估計總體的平均效應(yīng)和個體差異上具有獨特優(yōu)勢，可以通過多參數(shù)方法構(gòu)建立地指數(shù)模型[29]，有利于增加模型的評估精度。Wang等[30]利用啞變量與混合模型法分別模擬火炬松（Pinus taedaL.）樹高生長，最終表明混合模型方法準確性更好。因此，為了解釋不同立地水平的優(yōu)勢高生長差異，本研究采用單水平的非線性混合效應(yīng)模型進行立地隨機效應(yīng)模擬，確定系數(shù)（R2）從基礎(chǔ)模型擬合效果的0.4243～0.5644提高到0.5655～0.8089，以立地因子及其組合作隨機效應(yīng)的混合模型對比基礎(chǔ)模型優(yōu)勢顯著。

本研究以K-means聚類為基礎(chǔ)，將復(fù)雜的立地類型（ST）聚類劃分成立地類型組（GST），并在立地指數(shù)模型中考慮GST為隨機效應(yīng)，構(gòu)建多形立地指數(shù)模型（R2=0.8117）。此方法較原來基礎(chǔ)模型極大提高了模型預(yù)估精度，并且驗證了立地類型聚類構(gòu)造隨機效應(yīng)的方法具可行性。

4 結(jié)論

本研究采用非線性混合效應(yīng)模型方法，構(gòu)建了含立地隨機效應(yīng)的湖南丘陵平原區(qū)杉木多形立地指數(shù)曲線模型，利用Forstat求解模型參數(shù)，得出如下結(jié)論：海拔、坡度、坡向與土壤類型是影響湖南丘陵平原地區(qū)杉木優(yōu)勢高生長的顯著性因子，以顯著性因子及其組合作隨機效應(yīng)可以顯著提高立地指數(shù)建模精度，且包含立地類型的混合模型擬合精度最高，擬合精度的高低與主導(dǎo)立地因子的顯著性緊密相關(guān)；以K-means聚類將杉木78個初始立地類型劃分為11個立地類型組，包含立地類型組的混合模型在保留預(yù)估精度的前提下提高模型實用性。因此，基于立地隨機效應(yīng)的多形立地指數(shù)模型可以顯著提高建模精度，也為區(qū)域性復(fù)雜立地類型的立地質(zhì)量評價提供了思路與方法。