立地因子對檸條生物量動態變化的影響

王 博，許 昊，金學娟，王曉華，鎖 嵐

（1.寧夏大學 農學院，寧夏銀川 750021；2.寧夏大學 經濟管理學院，寧夏銀川 750021）

森林生物量是指在一定時空范圍內森林生態系統中植物積累的干物質總量，是研究碳循環的重要指標。森林在全球碳循環中的作用和初級生產力的模擬等均與森林生物量有關，準確預測森林生物量是進行森林資源管理和綜合評價的重要內容[1]。森林各部分生物量的大小可側面反映當地氣候、地形條件變化等對森林生態系統碳分配模式的影響。森林生物量預測不僅對研究森林生態系統碳儲量和固碳能力意義重大，也可為改善當地生態環境提供參考。因此，森林生物量預測一直受到林業學者的高度關注。

我國干旱半干旱地區的主要植被類型為灌木[2]，精準獲取具有時效性的灌木生物量數據是研究其生態效益和可持續經營的基礎。直接獲取樹木生物量數據對樹木破壞大，不適合大規模應用[3]。生物量法是最直接且應用最廣泛的方法之一[4]，通常以數學模型描述植物生長過程，具有破壞少、使用簡便且相對精度較高的特點[5-6]。受自身生長和環境等因素長期影響，植物在生長發育過程中通過不斷優化自身資源分配提高其適應環境變化的能力。灌木生物量是反映灌木對環境變化長期適應的重要指標[7]。目前，國內關于灌木生物量模型的報道較多，對各器官生物量模型的建立主要以林分因子為自變量[8-11]，且多以靜態模型為主[12-16]。王新云等[17]以人工檸條（Caraganakorshinskii）分枝長度與基徑平方的乘積為自變量，建立不同林齡人工檸條灌叢生物量模型；曾偉生等[2]采用非線性誤差變量聯立方程組方法，建立檸條單株相容性地上生物量、地下生物量和根莖比模型，并在此基礎上建立群落水平的生物量模型。部分學者研究立地因子對植物生物量的影響[18-20]；灌木生物量在立地因子及其交互作用影響下變化規律的研究較少[7，21]。鞠靈等[7]基于立地因子及其交互作用進行檸條生物量異速生長方程研究，但仍以林木因子為自變量，且未考慮時間因素；許昊等[21]僅研究不同立地條件對檸條生物量分配格局的影響，無法預測檸條生物量。本研究以立地因子及其交互作用為自變量，建立灌木生物量動態預測模型，研究灌木生長變化規律。

檸條為豆科（Leguminosae）錦雞兒屬灌木，又名毛條、白檸條，具有抗旱、耐貧瘠等特點，在惡劣條件下可正常萌發，擁有較強的適應性[22]。檸條作為我國干旱半干旱地區主要的造林樹種，在防風固沙、水土保持等生態恢復方面發揮重要作用[22]。本研究以寧夏靈武市白芨灘國家級自然保護區的檸條灌木林為研究對象，基于非參數多元方差分析，探究立地條件對檸條生物量的影響，在構建生物量回歸模型時，以啞變量的形式引入立地因子，提高模型精度[23]，為進一步研究灌木林碳儲量提供科學依據。

1 材料與方法

1.1 研究區概況

研究區位于寧夏靈武市白芨灘國家級自然保護區（106°00′～106°21′E，37°46′～38°41′N），總面積約為8.18×104hm2，地處毛烏素沙地邊緣；保護區以檸條灌木林為主，其分布面積約為1.73×104hm2，占保護區總面積的21.15%；保護區還分布著沙拐棗（Calligonummongolicum）、北沙柳（Salixpsammophila）和細枝巖黃耆（Hedysarumscoparium）等灌木樹種[21]。

1.2 數據來源

結合森林資源二類清查數據和相關規劃造林數據，依據檸條在寧夏中部干旱風沙區的分布，在靈武市白芨灘國家級自然保護區內種植密度（5 000叢/公頃）一致的林地中選擇不同種植年限、不同立地條件（海拔、坡度和坡向）且具有代表性的檸條林分。2021—2022 年7—8 月，共設置26 塊225 m2的樣地（表1）。在每個樣地內測定所有檸條的地徑和株高，取平均值。在每個樣地內選擇3～5 叢生長中等、無病蟲害的樣株，測定年輪；測定完成后，采用整株收獲法獲取樣株地上部分和地下1 m 深度根系，共采集108 叢。將所有樣品帶回實驗室，在65 ℃烘箱內烘干至恒重[24]，稱量干重。記錄樣地海拔、坡度和坡向；將各個立地因子標準化后的數值由小到大均分為不同等級。海拔分為高海拔（H=1 391～1 476 m）、中海拔（M=1 304～1 390 m）和低海拔（L=1 218～1 304 m）。坡向分為陽坡（南坡，S）和陰坡（北坡，N）。依據《測樹學》[25]，坡度可劃分為0°～5°、6°～15°、16°～25°、26°～35°、36°～45°和>45°六個等級；本研究區坡度為6°～35°，可分為緩坡（HS=6°～15°）、斜坡（XS=16°～25°）和陡坡（DS=26°～35°）。

表1 樣地基本特征和生物量統計Tab.1 Basic characteristics and statistics on biomass in sample sites

1.3 研究方法

1.3.1 數據檢驗與非參數多元方差分析

調查海拔、坡向和坡度3個立地因子，比較不同立地條件下檸條生物量的變化。常規的參數方差分析需滿足方差齊次性和正態分布[26]；當數據不能滿足這兩個條件時，采用非參數方差分析。基于R語言，采用Shapiro-Wilk 函數對數據進行正態分布檢驗，采用Fligner-Killeen 函數對數據進行方差齊次性分析（表2）。兩種檢驗中，總生物量的P值均小于0.05，參數方差分析不適用本研究數據，需考慮非參數方差分析。

非參數方差分析結果顯示，海拔及其與坡度和坡向的交互作用對檸條生物量均影響不顯著，坡向、坡度及其交互作用對總生物量均影響極顯著（P<0.01）；在建立生物量模型時，應考慮坡向、坡度及其交互作用的影響（表3）。

表3 非參數方差分析結果Tab.3 Results of non-parametric analysis of variance

1.3.2 生物量模型構建與模型精度檢驗

建立森林生長與收獲動態預測數學模型，確定適宜的生長方程，是建立森林生物量模型的重要環節。從理論上講，理想的生長方程需滿足較高的適應性和精度，且該方程參數具有一定的生物學意義。Logistic[27]方程可準確描述林木生長和收獲。本研究采用R 語言對檸條總生物量進行擬合，計算公式為：

式中，WT為總生物量（kg）；a、b和c為模型參數；t為林齡。

模型參數個數相同時，可通過平均絕對誤差（Mean Absolute Error，MAE）、均方根誤差（Root Mean Squared Error，RMSE）和調整后的決定系數（adj-R2）檢驗模型精度，確定最優模型；模型參數個數不同時，可對不同模型進行方差分析，選擇精度最高的模型作為最優模型。MAE、RMSE 越接近0，adj-R2越接近1，模型精度越高，擬合效果越好[7，28-29]。

式中，yij為因變量實際觀測值；ij為因變量預測值；為因變量實際觀測值均值；n為樣本總量；ni為第i個樣地的實際觀測數；r為模型參數個數。

2 結果與分析

2.1 立地因子及其交互作用的啞變量設置

為體現立地因子對檸條生物量的影響，將立地因子作為啞變量引入基礎方程，改進檸條生物量基礎模型。利用R語言的dummyVars函數進行啞變量設置。綜合考慮坡向、坡度及其交互作用，共設置11 個啞變量；其中，單一因子啞變量5 個（坡向2 個、坡度3個）、交互作用啞變量6個（表4）。

表4 坡向、坡度及其交互作用的啞變量設置Tab.4 Dummy variables of slope aspect,slope gradient and interactions

2.2 總生物量模型與啞變量模型設置

采用Logistic 方程建立檸條總生物量與時間的回歸關系，不同類型啞變量參數個數不同，對應的方程結構也不同。

以坡度為啞變量：

以坡向為啞變量：

以坡向與坡度交互作用為啞變量：

基于Logistic 方程，選用或構建9 種方程結構，方程1 為基礎方程結構（未添加啞變量），方程2～9為引入啞變量后的方程結構（表5）。

表5 總生物量各啞變量方程表達式Tab.5 Equation expressions of total biomass for dummy variables

采用MAE、RMSE 和adj-R2為評價指標，對模型的擬合能力進行對比分析，剔除不能收斂的模型（表6）。與基礎模型相比，引入啞變量模型的MAE和RMSE 均較小，說明基礎模型引入啞變量能明顯提高模型精度。在引入啞變量坡度和坡向的方程（2、3和4及5、6和7）中，方程4和7的MAE 和RMSE均最小，決定系數均最大，兩個方程均將啞變量引入模型參數c處；在考慮交互作用的方程8 和9 中，兩者的MAE 和RMSE 均低于引入單因素啞變量的方程，決定系數均高于引入單因素啞變量的方程，其中方程9（啞變量引入模型參數b處）的決定系數最高（adj-R2=0.887 5）。引入單因素啞變量方程的決定系數均表現為參數c>參數b>參數a；將雙因素啞變量引入模型參數c處的方程擬合結果不收斂，且在單因素方程參數擬合中，參數b與參數c的擬合精度相差較小（方程4 的決定系數比方程3 高0.09%，方程7 的決定系數比方程6 高0.09%），最終選擇方程9為最優方程。

表6 總生物量模型擬合結果Tab.6 Results of model fitting for total biomass

方程9 的參數估計值顯著性分析結果顯示，參數a和添加坡向與坡度交互作用啞變量參數b的T檢驗結果均大于0.05；參數c的檢驗結果在0.01 水平上差異顯著（表7）。

表7 方程9參數顯著性分析Tab.7 Significance analysis on parameters in equation 9

2.3 啞變量模型擬合曲線分析

圖1為沒有添加啞變量時總生物量與林齡的擬合曲線圖。

圖1 總生物量變化趨勢Fig.1 Change trends of total biomass

當坡度啞變量引入模型參數a處時，總生物量表現為緩坡>陡坡>斜坡，說明坡度越低，總生物量越大（圖2a）；當坡度啞變量引入模型參數b處時，總生物量表現為陡坡>斜坡>緩坡，且隨林齡增加，陡坡總生物量增加趨勢更明顯，說明在陡坡條件下檸條達到最大生物量的時間比其他坡度條件下短（圖2b）；當坡度啞變量引入模型參數c處時，總生物量表現為緩坡>陡坡>斜坡，且隨林齡增加，緩坡和陡坡總生物量增加趨勢更明顯（圖2c）。坡度啞變量模型以方程4表現最優。

圖2 添加坡度啞變量的總生物量變化趨勢Fig.2 Change trends of total biomass with slope gradient as dummy variable

無論將坡向啞變量引入模型參數a、b和c任何一個位置，總生物量均表現為陰坡>陽坡，且隨林齡增加，兩者差別變大；說明坡向對總生物量擬合曲線中3 個生物學參數的影響較一致（圖3）。坡向啞變量模型以方程7表現最優。

圖3 添加坡向啞變量的總生物量變化趨勢Fig.3 Change trends of total biomass with slope aspect as dummy variable

當坡向與坡度交互作用啞變量引入模型參數a和b處時，總生物量變化趨勢幾乎相同，交互作用啞變量擬合曲線結果與單因素啞變量擬合曲線結果有一定關聯（圖4）。當坡度啞變量引入模型參數a處時，總生物量表現為緩坡>陡坡>斜坡（圖2a）；當交互作用啞變量引入模型參數a處時，陰坡條件下總生物量也表現為緩坡>陡坡>斜坡（圖4a）。無論將交互作用啞變量引入模型參數a、b任何一個位置，緩坡+陰坡的總生物量均最大，且與其他5 個交互作用啞變量擬合曲線差別較大；陡坡+陽坡的總生物量均最小。交互作用啞變量在林木最大收獲量和其達到最大收獲量的時間上均表現較優。坡向與坡度交互作用啞變量模型以方程9表現最優。

圖4 添加坡向與坡度交互作用啞變量的總生物量變化趨勢Fig.4 Change trends of total biomass with interaction of slope aspect and slope gradient as dummy variable

2.4 模型差異顯著性分析

利用R 語言的方程差異性分析功能對方程1（基礎模型）、方程4（單因素啞變量最優模型）和方程9（雙因素啞變量最優模型）進行分析（表8）。3個模型間均呈極顯著差異（P<0.001）；方程4 與方程1的差異最大（F=12.512）。

表8 模型預測效果比較Tab.8 Comparisons on effects of model prediction

2.5 模型檢驗

采用留一處法對最優模型進行檢驗（表9）。添加了啞變量方程的RMSE 和MAE 均值均小于基礎方程。方程9 的RMSE 均值最小（0.197 4），表示林齡每增加1 年，檸條總生物量僅變化0.197 4 個單位；MAE均值也最小（0.001 8），驗證結果較好。

表9 留一處法檢驗結果Tab.9 Results of leave-one-out test

2.6 研究區檸條生物量動態預測

根據2022年森林遙感影像，獲得研究區內檸條分布圖；采用基于立地因子影響的總生物量啞變量最優模型（方程9），對檸條生物量進行30 年內動態預測。2022 年，總生物量約為1 500 kg/hm2；2037年，總生物量約為2 500 kg/hm2，預測在未來15 年間增長約66.66%；2052年，總生物量約為4 500 kg/hm2，預測在未來30年間增長約200%（圖5）。

圖5 30年總生物量動態預測Fig.5 Dynamic prediction of total biomass after 30 years

3 討論與結論

坡度、坡向和海拔等立地因子與當地氣候、土壤環境等關聯密切。閔梓驍等[30]認為，陰坡土壤水分情況優于陽坡，沙棘（Hippophaerhamnoides）在陰坡和坡位較低條件下生長較好，與本研究結果相似。周振釗等[20]采用逐步回歸分析法研究立地因子對林下灌木生物量的影響，結果顯示林下灌木生物量與坡向和海拔均呈顯著相關；本研究結果顯示，海拔及其與坡度和坡向的交互作用對檸條生物量均影響不顯著，坡向及其與坡度的交互作用均影響極顯著。羅鵬飛等[31]研究立地因子對桉樹（Eucalyptusspp.）胸徑生長的影響，發現坡向和坡度對桉樹胸徑生長影響顯著，與本研究結果相似。候曉巍等[32]認為，陽坡條件下，祁連圓柏（Juniperusprzewalskii）生態系統碳密度略高于陰坡，與本研究結果相反；這是由于祁連圓柏屬喜光樹種，檸條屬高耗水樹種，陽坡光照條件優于陰坡，陰坡土壤水分條件優于陽坡。這些研究結果表明，立地條件不僅影響植被生長且對不同植被有不同影響。種群內部的競爭關系也對植被生長有一定影響，許昊等[33]考慮檸條種植行間距對其株高生長的影響，結果顯示種植行間距對檸條株高生長曲線的漸進值影響更大，種植行間距越大（密度越小），檸條株高越高。本研究為消除種植行間距的影響，選擇種植行間距相同的檸條林分開展研究。

本研究僅針對寧夏靈武市白芨灘國家級自然保護區的檸條灌木林進行研究，海拔分布范圍較小，這是海拔對于檸條生物量影響不顯著的重要原因。將坡度作為定量因子時，啞變量數量過多，各啞變量模型擬合圖間的差異難以辨別，無法得出理想結果；本研究將坡度劃分為與坡向相同的定性因子，將其分為3 個梯度，各梯度間差異較大，繪制的模型擬合曲線圖效果更佳。李盈等[19]采用賦值方法對立地因子進行劃分，優化灌木生物量異速生長方程，該方法可有效提高模型精度。立地因子數據在回歸模型中屬于定性因子，是不等距的，賦值的方式具有一定局限性。將定性因子直接加入回歸模型中，不滿足常規的回歸分析要求，也很難對結果進行解釋。隨著建模技術的深入發展，在回歸模型中以啞變量的形式引入定性因子能很好地解決以上問題。

本研究采用Logistic 方程構建回歸模型，擬合曲線呈S 形，結果顯示檸條總生物量在未來30 年的模擬過程中持續增大。這可能是因為試驗區常年日照時長較長，土壤含水量低，檸條生長緩慢，未在30 年內達到K值。模型擬合過程中，坡向與坡度交互作用能有效提高模型精度；在生物量回歸方程中，應結合當地條件適當考慮以立地因子交互作用為啞變量。本研究中，由于實際觀測林齡不夠全面，擬合出的模型精度均小于0.900，一般認為調整后的決定系數在0.900 以上的數學模型因果關系較強[34]，需在日后研究中繼續對靈武市白芨灘地區檸條灌木林進行更精確的觀測，以獲得更全面的數據。

海拔及其與坡度、坡向的交互作用對檸條生物量均影響不顯著，坡度和坡向及其交互作用對檸條生物量均影響極顯著。與基礎模型相比，引入啞變量明顯提高模型精度。剔除不收斂的模型，方程9表現最優（MAE=0.191 4、RMSE=0.278 9、adj-R2=0.887 5），擬合精度最高。

利益沖突：所有作者聲明無利益沖突。

作者貢獻聲明：王博負責研究計劃制定、外業調查、數據分析、文獻檢索和論文撰寫；許昊負責研究計劃指導和論文指導與修改；金學娟、王曉華和鎖嵐負責外業調查和數據整理。