含間伐強度啞變量的關帝山油松天然林單木斷面積生長模型構建

李學輝，白 寧，郭晉平，馬建峰，張蕓香

(1.山西農業大學林學院，山西 太谷 030801；2.關帝山國有林管理局，山西 文水 032100)

森林生長模型可以用來預估林分的生長和收獲量，單木生長模型可以預測各單株林木的生長狀況和生長潛力[1-4]，構建特定森林樹種的單木生長模型對于區域森林經營具有重要意義[5]。在生長模型中用胸高斷面積反映林木生長狀況具有較高的準確性和穩定性，是構建林分和單木模型的合適指標[2，6]。根據競爭指標是否含有對象木與競爭木之間相對位置的信息，將單木模型分為距離有關單木模型和距離無關單木模型[7-9]。距離無關單木生長模型會導致相同大小林木若干年后仍為同樣大小林木的結果，也無法反映林木在林分間伐前后競爭壓力變化[10]。距離有關單木模型可以反映林木在林分中所處不同小生境的差異，并能較準確地反映競爭木間伐對主林木的影響[11]。

單木模型的建模方法包括潛在生長量修正法、回歸估計法和生長分析法。潛在生長量修正法需要確定無競爭壓力的疏開木，利用單木競爭指標對每株林木的潛在生長函數進行調整和修正，由于疏開木難以確定，因此該方法使用較少。回歸估計法利用多元回歸建立林木生長量與林木大小、林木競爭狀態和所處立地條件等因子之間的回歸方程，該方法建立的單木模型適應性差，預估能力依賴建模樣本數據，且方程參數沒有生物學意義。生長分析法以理論生長方程為基礎模型，通過分析參數與單木競爭之間的關系構造單木模型，該方法不依賴疏開木的生長，選擇合理的理論生長模型，可以獲得良好的預測效果[12-14]。

在影響林木生長發育的因子中，不僅有許多定量因子，還有一些定性因子需要在建模中加以考慮，而構建含啞變量的模型可以有效表示定性因子[15-16]，因此，在單木模型中加入立地類型、林分類型、競爭類型等定性因子作為啞變量可以有效提高模型預測效果[17]。撫育間伐是基本的森林經營措施，顯著影響林木胸徑、斷面積和材積生長，間伐強度是關鍵影響因素[18]。但將間伐強度作為啞變量加入到單木模型中的研究還鮮見報道。

油松(PinustabulaeformisCarr)為松科松屬針葉常綠喬木，是我國特有樹種，是黃土高原地區最主要的造林樹種之一，油松天然林不僅在山西集中分布，在遼寧、河北、北京、內蒙古、陜西等地區也有廣泛分布，在構建區域尺度景觀多樣性、保持水土和維持森林生態系統的穩定性等方面發揮著重要作用，具有較高的研究價值[19]。因此，通過選擇合理的理論生長方程和參數，提高油松單木模型精度是當前在構建單木生長模型中亟待解決的問題。

本研究以山西關帝山林區油松天然林為研究對象，采用生長分析法的建模技術，將間伐強度作為啞變量引入模型，建立距離有關單木斷面積生長模型，探究單木生長模型中加入間伐強度啞變量的建模途徑和效果，構建含間伐強度啞變量的最優單木斷面積生長模型，為油松天然林的經營管理提供科學依據。

1 研究區概況

研究區位于呂梁山中段關帝山林區(111°21′～111°37′E，37°45′～37°59′N)，是呂梁山斷裂隆起的背斜構造山地，地形破碎，山體陡峻，溝壑縱橫。林區氣候受季風影響和控制，屬于暖溫帶大陸性山地氣候，夏季濕潤多雨，冬季寒冷干燥，年均溫4.2°C，年均降水量600～820 mm，年均相對濕度70.9%，年無霜期100～130 d，雪期和凍土期6個月，森林覆蓋率71%。

孝文山林場位于關帝山林區中部，地處111°24′～112°37′E，37°41′～37°54′N，南北長約2 km，東西寬約17 km，海拔 1 345～2 659 m。主要森林類型建群樹種有油松、華北落葉松(LarixprincipisvupprechtiiMayr)、青杄(PiceawilsoniiMast)、白杄(PiceameyeriRehd.et Wils)、白樺(BetulaplatyphyllaSuk)、遼東櫟(QuercusliaotungensisMayr)等；主要林下灌木樹種有三椏繡線菊(SpiraeatrilobataL)、黃刺玫(RosaxanthinaLinn.)、山刺玫(RosadavuricaPall.)、胡枝子(LespedezabicolorTurcz.)、水栒子(CotoneastermultiflorusBge.)等；主要林下草本植物有披針苔草(CarexlanceolataBoott)、蕨類(Pteridium)、地榆(SanguisorbaofficinalisLinn.)、金蓮花(TrolliuschinensisBunge)等。

2 研究方法

2.1 樣地設置與指標測定

于2019年7月對關帝山林區孝文山林場五十溝的9塊固定樣地進行復測，復測項目包括：號牌編碼、胸徑、單木位置和生長狀況。經復測后根據林分內間伐(2013年)的林木斷面積與林分總斷面積之比，確定間伐強度，分為重度間伐、中度間伐、輕度間伐3類。并在每個樣地內確定3株優勢木，計算各樣地優勢木平均高。樣地調查間隔期為5年。樣地概況見表1。

表1 樣地概況

2.2 對象木和競爭木確定

將樣地內所有胸徑≥5 cm的活立木作為對象木，靠近樣地邊界的對象木采用擴大緩沖區的辦法確定競爭木[20]，緩沖區設在樣地外圍8 m范圍內。9塊樣地的對象木共360株。

競爭木確定采用固定半徑法[21]。為確定合理的競爭木測定范圍，選取7株對象木，按競爭指數計算方法(公式1)，以1～12 m共計12個固定半徑，分別計算其競爭指數，繪制固定半徑與競爭指數增量的關系圖(圖1)，確定本研究競爭木的固定半徑為6.4 m。

圖1 不同半徑的競爭指數增量Fig.1 Increment of competition index on different radius

2.3 樣地調查取樣和對象木胸高斷面積測算

對樣地內樹高≥1.3 m的活立木進行每木編號和檢尺，測定其胸徑、樹高(圖帕斯TurPulse 200激光測距儀)和冠幅，并測定其相對位置，繪制林木定位圖，用生長錐鉆取對象木樹芯，確定其年齡。由胸徑計算出每株立木的胸高斷面積。

2.4 競爭指數計算

從立木定位圖測出對象木與各競爭木之間的距離，采用公式(2)(HEGYI公式)，計算出每株對象木的競爭指數：

(1)

式中：CI為某株對象木的競爭指數，Di為樣地中第i株對象木的胸徑，Dj為該對象木第j株競爭木的胸徑，dij為該對象木i與其第j株競爭木之間的距離。

3 模型構建

3.1 基礎模型構建

選用Logistic方程、Korf方程、Richards方程、Weibull方程，考慮立地質量、競爭指數、林木年齡3個因素。其中，立地質量用立地指數表達，本研究用樣地優勢木平均高代替地位指數。定義參數A是林木斷面積的最大值[22]，受立地質量的影響；參數c與生長速度相關[23]，受競爭壓力的影響。構造出4個單木斷面積生長的基礎模型表達式(表2)。

表2 基礎模型表達式

3.2 間伐啞變量設置

將間伐強度作為啞變量引入單木斷面積生長模型，用δ(x,i)表示[24]。其公式為：

(2)

3.3 模型參數求解和精度檢驗

運用R語言進行數據分析，采用最小二乘法求解模型參數。

為檢驗模型精度，以決定系數R2和預估精度P值選擇模型，并對所選模型進行F檢驗，判斷所選模型回歸效果是否顯著；選擇總相對誤差(TRE)、平均絕對誤差(MAE)、均方根誤差(RMSE)作為模型精度檢驗指標，比較基礎模型和含啞變量模型的模擬效果。

(3)

(4)

P-1)

(5)

(6)

(7)

(8)

4 結果分析

4.1 基礎模型擬合和精度檢驗

以關帝山360株油松為基礎，利用基礎模型表達式對樹齡、優勢木平均高及競爭指數進行擬合，結果見表3。

表3 基礎模型參數與模型精度擬合

4種基礎模型的R2為0.46～0.52，P值為0.89～0.92。各基礎模型擬合度以Logistic模型最高，為 0.514 5；P值以Korf模型最高，為 0.911 7。

4.2 含間伐強度啞變量模型擬合與檢驗

將間伐強度分成3級，將間伐強度啞變量加入到模型的不同參數及其組合上，共構建了28種參數組合，利用360株對象木的生長數據進行參數估計，剔除不收斂參數組合，以參數顯著性t檢驗和P值大小(以P≤0.05為標準)，剔除參數不符合的模型，再根據模型決定系數R2確定斷面積生長模型的最優啞變量參數組合，進行模型選優。模型R2和P值見表4。

表4 含間伐強度啞變量模型精度擬合

含間伐強度啞變量的模型與基礎模型相比P值和R2均有顯著提高。含間伐強度啞變量的W模型R2總體來說大于其他含間伐強度啞變量模型，含間伐強度啞變量的K模型P值均大于其他含間伐強度啞變量模型。Logistic模型將間伐強度啞變量加到a、a2參數效果最優(L5)；Weibull模型將間伐強度啞變量加到a、a1、a2參數效果最優(W7)；Korf模型將間伐強度啞變量加到a1、a2參數效果最優(K6)；Richards模型將間伐強度啞變量加到a、a1、a2參數效果最優(R7)。

4.3 基礎模型與含間伐強度啞變量模型的比較

為進一步確定基礎模型與含間伐強度啞變量模型的精度差異，對基礎模型和各含間伐強度啞變量最優模型的總相對誤差(TRE)、平均絕對誤差(MAE)及均方根誤差(RMSE)3個精度檢驗指標進行對比分析(表5)。

表5 基礎模型與含間伐強度啞變量模型擬合效果對比

從表5模型模擬結果可以看出，含間伐強度啞變量模型的TRE、MAE、RMSE檢驗指標值均小于基礎模型的檢驗指標值。

為更直觀地反映基礎模型與含間伐強度啞變量模型擬合效果差異，建立斷面積殘差分布圖(圖2)和斷面積實測值與預測值相關關系圖(圖3)。

由圖2、圖3可以看出，含間伐強度啞變量模型要優于基礎模型，其擬合效果明顯優于基礎模型。說明含間伐強度啞變量斷面積生長模型模擬效果要優于基礎模型的模擬效果，且含間伐強度啞變量的Weibull模型最優。利用F檢驗對模型進行檢驗，檢驗結果F=57.04>F0.95(9，337)=1.91，說明回歸效果顯著。最優模型為：

圖2 含間伐強度啞變量模型殘差Fig.2 Residual Diagram of dummy variable model with thinning intensity

圖3 含間伐強度啞變量預測值與實測值相關關系Fig.3 Correlation between predicted value and measured value with dummy variable of thinning intensity

g=(a+c1K1+c2K2)SI[1-exp-(a1+c3K1+c4K2)t(a2+c5K1+c6K2)CI]

5 討論與結論

5.1 討論

使用固定半徑法確定競爭木時，通常依據樹種的不同來研究確定一個合適的競爭半徑。本實驗將油松競爭半徑確定為6.4 m。隨著目標樹近自然經營在國內的發展，通過目標樹競爭距離構建目標樹競爭強度與目標樹生長關系尤為重要，油松的斷面積受競爭強度影響顯著，通過目標樹與競爭木之間距離的改變，目標樹空間分布格局發生變化，林木的空間分布格局直接或間接影響林木生長及林下空間合理分配[25]。隨著林木的生長發育導致樣地立地條件改變、林木年齡增長、人為采伐等，目標樹競爭壓力必然會發生變化，在固定半徑內的目標樹競爭壓力變化規律還有待進一步研究。因此通過對目標樹競爭距離的研究，可以較好地為森林近自然經營和森林撫育間伐提供更加科學的依據。

目前并沒有一個單一的競爭指數或單一類別的指數是普遍優越的，指標的表現因森林類型和立地條件而異[26]。Hegyi簡單競爭指數是目前應用比較廣泛的一類指標，其從胸徑與距離的角度構建了林木之間的競爭關系。有學者研究認為利用樹冠重疊指數可以較好地反映競爭關系，因此在進一步研究林木競爭指標中可以綜合考慮指標與胸徑、樹冠、樹高之間的相關關系，從而研究林木水平和垂直空間的綜合競爭狀況。

不同間伐強度的樹種結構、徑階結構及對象木、目標樹所受的競爭壓力均表現出差異性，導致林木生長發育的差異，采用林分尺度上的間伐強度對林木生長發育影響同樣顯著，但對于單木尺度的間伐強度對油松斷面積的影響本實驗也做了相關分析，結果并不理想，相關問題有待進一步研究。

含間伐強度啞變量模型可以很好地融合于油松單木模型，提高斷面積建模精度，擴大模型的相容性。由于此次研究林木數據有限，缺乏有效的適用性檢驗，以后的研究中可以在增加樣本數量的基礎上，研究更多不同類型的單木斷面積生長模型，以增加模型的適用范圍。

5.2 結論

采用4種理論生長方程建立關帝山油松與距離有關的單木斷面積生長基礎模型，各基礎模型預測精度均在90%以上，預測效果良好。

在考慮林分間伐效應的基礎上，將間伐強度作為啞變量引入基礎模型，構建含間伐強度啞變量與距離有關的油松單木斷面積生長模型，引入啞變量后，各模型預測精度顯著提高，其中Weibull模型引入間伐強度啞變量模型擬合度和預測精度最優，分別為0.606 1、99%，TRE、MAE及RMSE檢驗指標值最小，含間伐強度啞變量的Weibull模型可作為關帝山油松最優單木斷面積生長模型，啞變量最優組合方式為a、a1、a2。

本研究建立的模型對關帝山地區油松天然林是適用的，可以模擬其單木的生長，從而為油松林分的經營管理決策提供參考依據。