宋 爭,黃 朗,胡 松,張雄清,朱光玉*
(1. 中國林業科學研究院林業科技信息研究所,北京 100091;2. 中南林業科技大學林學院,湖南 長沙 410004;3. 中國林業科學研究院林業研究所,北京 100091)
立地質量評價是對森林立地的宜林性或潛在生產力判斷或預測的基礎工作[1],立地指數是一種最普遍的評價森林立地質量的方法[2-3],分析某樹種林分基準年齡時對應的林分優勢高可用于分析森林立地生產潛力的基礎[4]。除了根據優勢高與年齡指標,還能夠通過樣地實測值和收獲表蓄積量估算獲得立地指數[5]。但由于立地條件的差異性或多樣性,基于優勢高-年齡的導向曲線模擬精度往往較低,同一樹種在不同區域的立地指數估計結果隨立地質量的不同而存在顯著差異。因此,很多學者在評價立地質量時,優先分析環境要素(氣候、植被、地形等)與立地指數的關系,得到包含環境因子的立地指數函數關系[6-7]。也有學者采用樹高-年齡分級啞變量的方法,用樣地林木平均高劃分等級代替立地因子評價立地質量[8]。這些方法雖然提高了立地指數建模精度,但關于環境因子在模型中的表達仍需進一步探索。
目前,如何建立立地指數模型成為描述森林立地質量的重點和中心問題。研究指出有單形和多形兩種關系,其中多形立地指數曲線結合優勢木平均高和其他反映立地條件因子,共同評定立地質量[9-11]。最常用的2種擬合方法為固定基準年齡的靜態方程和可變基準年齡的動態方程,前者依賴于選定的基準年齡,后者局限于嚴格的數據要求[12]。而混合效應模型是一種最常用的處理局部變異性的建模方法[13-14]。對于不同立地類型的優勢高生長,混合模型不僅能夠描述整體變化的平均趨勢,還能夠給出協方差和方差等相關數據以了解個體特點[15-16]。雖然在構建立地指數模型時,相關文獻運用數量化方法I[17-18],然而該方法前提是固定效應,如果考慮隨機效應,該方法不適用。
鑒于此,本研究以湖南丘陵平原區杉木(Cunninghamia lanceolata(Lamb.) Hook.)人工林為研究對象,首先利用數量化方法I確定影響杉木優勢高生長的顯著性影響因子;然后以主導立地因子及其組合作隨機效應,建立基于立地隨機效應的杉木多形立地指數曲線模型,通過模型選優對不同立地因子水平的立地質量進行評價;最后以Kmeans聚類劃分立地類型組,解決復雜立地類型的模型應用問題,以提高立地指數曲線模型的適用性和準確性。
湖南位于中亞熱帶,地貌大部分為丘陵,呈現大陸性季風濕潤氣候,全年平均氣溫17.5 ℃左右,平均降水量1200~1700 mm,日照時長1300~1800 h,雨水豐富且熱量充沛。研究區喬木樹種以杉木(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.)等。
以湖南丘陵平原區杉木人工純林為研究對象,在湘鄉、懷化、會同、株洲、永州、臨湘、綏寧、桃源等26個縣,設置樣地(20 m × 30 m)共360塊,記錄相關立地因子,并對樣地內活立木每木檢尺。樣地內每0.01 hm2選取最粗優勢木1株,計算優勢木平均高(H)。相關林分調查因子統計詳見表1。

表1 杉木林分調查因子統計Table 1 Statistics of stand factors for Cunninghamia lanceolata
1.3.1 立地分類方法 選取如表2所示的影響杉木優勢高生長的6個立地因子,根據《中國森林立地》相關標準進行分級,將海拔每100 m劃分為一級,分析不同海拔的差異。

表2 立地因子等級劃分Table 2 The division of site factor grades
立地因子作自變量,各林分優勢木平均高做因變量,運用數量化方法Ⅰ對杉木立地因子分別分析。基于方差分析中不同“Pr>F”大小,篩選顯著因子并進一步得到主導因子。
1.3.2 非線性混合效應模型 為更準確地評價杉木立地質量,本研究考慮立地因子作隨機效應對杉木立地指數模型的影響,隨機效應包括了立地因子與其交互作用。非線性混合效應模型的形式為[19-20]:

式中:Hij為 第i個立地類型下第j個樣地的平均優勢高,Tij為第i個 立地類型下第j個 樣地的林齡,φij為參數向量,εij是誤差項,M是立地類型數量,ni為第i個立地類型的樣地數,β為(p×1)維固定效應向量,bi為(q×1) 維隨機效應向量,Aij,Bij為設計矩陣,D為隨機效應的方差—協方差矩陣,Ri為立地類型內的方差—協方差結構,σ2為方差,Gi為描述方差異質性的對角矩陣,Γi為描述隨機效應自相關性的方差矩陣。
構建非線性混合效應模型,需要對固定與隨機效應參數進行構造。通常情況下對所有可能的參數考慮隨機效應,但這會造成模型不收斂[21-22]。為避免模型不收斂問題,本研究選擇在1個參數上考慮立地的隨機效應。其次,由于本研究中林分優勢高與年齡不是重復調查數據,不存在樣地內測量值之間的相關性。因此,假設誤差項εij服從正態分布。
最后,本模型只考慮單個隨機效應參數來確定隨機效應參數方差,因此假定隨機參數構造類型為無結構,并采用AIC、BIC、對數似然比(Loglikelihood)和確定系數R2選擇最優混合效應模型。AIC、BIC值越小,Log-likelihood和R2越大,模型的擬合效果越好。
基礎模型不僅要能夠很好描述林木優勢高生長規律,還要具有生物學的可解釋性[23-24]。其中立地指數曲線基礎模型常為“S”形,本研究選擇常用的立地指數曲線方程(M1~M8),擬合杉木林分優勢高(見表3)。

表3 候選基礎模型Table 3 Candidate basic models
利用確定系數(R2) 、預估精度(P)、平均絕對誤差(MAE)和均方根誤差(RMSE)進行模型評價。

式中:yi為 第i個 樣地林分優勢高實測值,為第i個樣地林分優勢高估計值,n為樣本數,tα為置信水平α=0.05時的t分布值,p為模型參數個數,為平均預估值。
林分優勢木平均高受多種立地因子綜合作用影響。根據優勢木平均高與不同立地因子的相互關系,采用數量化方法Ⅰ對其進行顯著性檢驗,結果(表4)表明:對優勢木平均高具有顯著影響的立地因子有海拔、坡度、坡向與土壤類型,其顯著性順序為土壤類型>海拔>坡向>坡度。

表4 立地因子的顯著性檢驗Table 4 Significance testing of site factors
利用 Forstat2.1軟件對候選基礎模型M1~M8進行估計,相關基礎模型參數估計值、確定系數(R2) 、預估精度(P)、平均絕對誤差(MAE)及均方根誤差(RMSE)如表5所示。

表5 候選基礎模型擬合結果Table 5 The fitting results of candidate basic models
從8個候選基礎模型(表5)擬合結果可知,這8個模型的擬合精度均較低(R2=0.4243~0.5644),也證明了區域性立地指數建模精度低的問題。其中,模型M4擬合效果最好,其確定系數(R2= 0.5644)與預估精度(P=0.9750)最大,平均絕對誤差(MAE=6.7626)及均方根誤差(RMSE=2.4100)最小。
以主導立地因子作隨機效應,考慮立地單因素與多因素交互作用,構建立地指數混合模型。結合基礎模型參數個數和隨機效應因子組合類型,擬合所有參數及其組合形式,最優隨機效應立地指數模型通過AIC和BIC等評價指標獲取(表6)。
計算結果顯示,除HB、PX、TL等單因子在其它參數組合的擬合效果較優外,受水平數影響,其它隨機效應組合均只在漸進參數a上時,模型收斂。因此,含立地隨機效應的杉木多形立地指數模型表達式為:

式中:a,b,c為模型固定效應參數;μi為隨機效應參數;H為林分平均優勢木高;Age為林分年齡。
擬合多種隨機效應因子組合類型后的結果如表6所示,AIC、BIC相對于基礎模型M4都顯著下降,Log-likelihood和R2均大幅提高,這說明基于立地隨機效應的混合模型明顯優于基礎模型。此外,混合模型擬合效果的好壞,取決于立地因子顯著性的高低,立地因子越顯著,相應隨機效應組合類型的混合模型擬合結果越好。其中,海拔、坡度、坡向、土壤類型4個主導立地因子交互作用的混合模型AIC、BIC最低,Log-likelihood及R2最高。根據立地分類的原則與方法,4個主導立地因子排列組合可得到立地類型,綜上所述含立地類型的混合模型為最優模型。

表6 基于立地效應的混合模型擬合結果Table 6 The fitting results of mixed effects models considering site effect
依據海拔、坡度、坡向、土壤類型4個主導因子,將研究區劃分為78個立地類型,那么相應的混合模型隨機效應參數有78個,勢必會影響立地指數模型的應用。為了簡化立地類型,將初始立地類型(ST)進行K-means聚類。參考森林立地分類的標準,確定系數≥0.99為聚類精度標準,合并得分值相近的立地類型為立地類型組(GST),使得本研究區的立地類型聚類為11種(表7)。

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

表8 模型精度評價Table 8 Model evaluation
利用該模型繪制湖南丘陵平原區杉木多形立地指數曲線(圖1)。圖中展示了不同立地類型組的立地指數曲線圖,各立地類型組的模型模擬效果均較好。

圖1 多形立地指數曲線Fig. 1 Polymorphic site index curve
林分優勢高生長受林分密度影響很小,是同齡林立地質量評價采用最廣泛的指標,但關于林分優勢高的定義及測定標準并不一致[25]。趙美麗等[26]對樟子松(Pinus elliottiiEngelm.)同齡林最高、最粗木平均優勢高的相關性、穩定性研究后認為,最粗6株木法能科學得到林分優勢高,該方法具有簡易、穩定的優點,所以本研究以每100 m2選擇1株最粗木的方法計算平均優勢高。
在森林立地分類與評價的因子中,地形屬間接因素,但綜合反映著氣候和土壤特性。本研究在進行顯著性研究中用數量化方法I,以分析不同立地因子對杉木優勢高生長的影響,結果表明,海拔、土壤類型、坡向、坡度與杉木優勢高生長有強相關性。不同立地因子對優勢高生長的影響也不同。產生顯著性差異的原因有以下幾點:從海拔來說,南方丘陵平原地區熱量和降水量的垂直分異總的來說決定了研究對象的生長與分布,海拔為一個主要因素。從土壤類型與坡度來說,杉木對土壤濕度、養分、通氣性等要求較高,杉木生長由土壤類型的差異直接影響,坡度間接表現土壤特性。從坡向來說,不同坡向熱量光照差異,對杉木不足以引起其生長的很大差異,故坡向相比較而言,顯著性影響要低于土壤類型、海拔,但高于坡度的影響作用。王冬至等[27]在研究塞罕壩華北落葉松(Larix principisrupprechtiiMayr)人工林地位指數模型時,將林分優勢高與不同立地因子進行相關性分析,發現林分優勢高與立地因子顯著相關。
基于林分優勢高-年齡關系的立地指數,已成為表征森林立地生產力的一個重要參數,應用該方法的前提在于導向曲線具有較高的模擬精度。然而對于某一特定的樹種,導向曲線的形狀會因氣候、土壤、地形等的不同而有所變化,使用同一導向曲線模擬樹高生長必然會造成較大誤差[28]。而混合模型在估計總體的平均效應和個體差異上具有獨特優勢,可以通過多參數方法構建立地指數模型[29],有利于增加模型的評估精度。Wang等[30]利用啞變量與混合模型法分別模擬火炬松(Pinus taedaL.)樹高生長,最終表明混合模型方法準確性更好。因此,為了解釋不同立地水平的優勢高生長差異,本研究采用單水平的非線性混合效應模型進行立地隨機效應模擬,確定系數(R2)從基礎模型擬合效果的0.4243~0.5644提高到0.5655~0.8089,以立地因子及其組合作隨機效應的混合模型對比基礎模型優勢顯著。
本研究以K-means聚類為基礎,將復雜的立地類型(ST)聚類劃分成立地類型組(GST),并在立地指數模型中考慮GST為隨機效應,構建多形立地指數模型(R2=0.8117)。此方法較原來基礎模型極大提高了模型預估精度,并且驗證了立地類型聚類構造隨機效應的方法具可行性。
本研究采用非線性混合效應模型方法,構建了含立地隨機效應的湖南丘陵平原區杉木多形立地指數曲線模型,利用Forstat求解模型參數,得出如下結論:海拔、坡度、坡向與土壤類型是影響湖南丘陵平原地區杉木優勢高生長的顯著性因子,以顯著性因子及其組合作隨機效應可以顯著提高立地指數建模精度,且包含立地類型的混合模型擬合精度最高,擬合精度的高低與主導立地因子的顯著性緊密相關;以K-means聚類將杉木78個初始立地類型劃分為11個立地類型組,包含立地類型組的混合模型在保留預估精度的前提下提高模型實用性。因此,基于立地隨機效應的多形立地指數模型可以顯著提高建模精度,也為區域性復雜立地類型的立地質量評價提供了思路與方法。