王 博，許 昊，金學(xué)娟，王曉華，鎖 嵐

（1.寧夏大學(xué) 農(nóng)學(xué)院，寧夏銀川 750021；2.寧夏大學(xué) 經(jīng)濟(jì)管理學(xué)院，寧夏銀川 750021）

森林生物量是指在一定時空范圍內(nèi)森林生態(tài)系統(tǒng)中植物積累的干物質(zhì)總量，是研究碳循環(huán)的重要指標(biāo)。森林在全球碳循環(huán)中的作用和初級生產(chǎn)力的模擬等均與森林生物量有關(guān)，準(zhǔn)確預(yù)測森林生物量是進(jìn)行森林資源管理和綜合評價的重要內(nèi)容[1]。森林各部分生物量的大小可側(cè)面反映當(dāng)?shù)貧夂?、地形條件變化等對森林生態(tài)系統(tǒng)碳分配模式的影響。森林生物量預(yù)測不僅對研究森林生態(tài)系統(tǒng)碳儲量和固碳能力意義重大，也可為改善當(dāng)?shù)厣鷳B(tài)環(huán)境提供參考。因此，森林生物量預(yù)測一直受到林業(yè)學(xué)者的高度關(guān)注。

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

1.3 研究方法

1.3.1 數(shù)據(jù)檢驗與非參數(shù)多元方差分析

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

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

表3 非參數(shù)方差分析結(jié)果Tab.3 Results of non-parametric analysis of variance

1.3.2 生物量模型構(gòu)建與模型精度檢驗

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

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

模型參數(shù)個數(shù)相同時，可通過平均絕對誤差（Mean Absolute Error，MAE）、均方根誤差（Root Mean Squared Error，RMSE）和調(diào)整后的決定系數(shù)（adj-R2）檢驗?zāi)Ｐ途?，確定最優(yōu)模型；模型參數(shù)個數(shù)不同時，可對不同模型進(jìn)行方差分析，選擇精度最高的模型作為最優(yōu)模型。MAE、RMSE 越接近0，adj-R2越接近1，模型精度越高，擬合效果越好[7，28-29]。

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

2 結(jié)果與分析

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

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

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

2.2 總生物量模型與啞變量模型設(shè)置

采用Logistic 方程建立檸條總生物量與時間的回歸關(guān)系，不同類型啞變量參數(shù)個數(shù)不同，對應(yīng)的方程結(jié)構(gòu)也不同。

以坡度為啞變量：

以坡向為啞變量：

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

基于Logistic 方程，選用或構(gòu)建9 種方程結(jié)構(gòu)，方程1 為基礎(chǔ)方程結(jié)構(gòu)（未添加啞變量），方程2～9為引入啞變量后的方程結(jié)構(gòu)（表5）。

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

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

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

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

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

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

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

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

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

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

無論將坡向啞變量引入模型參數(shù)a、b和c任何一個位置，總生物量均表現(xiàn)為陰坡>陽坡，且隨林齡增加，兩者差別變大；說明坡向?qū)偵锪繑M合曲線中3 個生物學(xué)參數(shù)的影響較一致（圖3）。坡向啞變量模型以方程7表現(xiàn)最優(yōu)。

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

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

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

2.4 模型差異顯著性分析

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

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

2.5 模型檢驗

采用留一處法對最優(yōu)模型進(jìn)行檢驗（表9）。添加了啞變量方程的RMSE 和MAE 均值均小于基礎(chǔ)方程。方程9 的RMSE 均值最?。?.197 4），表示林齡每增加1 年，檸條總生物量僅變化0.197 4 個單位；MAE均值也最?。?.001 8），驗證結(jié)果較好。

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

2.6 研究區(qū)檸條生物量動態(tài)預(yù)測

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

圖5 30年總生物量動態(tài)預(yù)測Fig.5 Dynamic prediction of total biomass after 30 years

3 討論與結(jié)論

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

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

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

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

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

作者貢獻(xiàn)聲明：王博負(fù)責(zé)研究計劃制定、外業(yè)調(diào)查、數(shù)據(jù)分析、文獻(xiàn)檢索和論文撰寫；許昊負(fù)責(zé)研究計劃指導(dǎo)和論文指導(dǎo)與修改；金學(xué)娟、王曉華和鎖嵐負(fù)責(zé)外業(yè)調(diào)查和數(shù)據(jù)整理。