立地因子对柠条生物量动态变化的影响

王 博，许 昊，金学娟，王晓华，锁 岚

（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），拟合精度最高。

利益冲突：所有作者声明无利益冲突。

作者贡献声明：王博负责研究计划制定、外业调查、数据分析、文献检索和论文撰写；许昊负责研究计划指导和论文指导与修改；金学娟、王晓华和锁岚负责外业调查和数据整理。