Growth Models of Stand Volume and Biomass of Coniferous Species in Guangdong Province

广东省位于20°13′～25°31′ N和109°39′～117°19′ E，处于东亚季风区，包括了热带、中亚和南亚热带3种气候，年均气温22.3℃，月降水量约313 mm，年日照时数约1 746 h，太阳辐射总量约4 200～5 400 MJ·m⁻²，生长期水热充足，利于植被生长。地貌类型复杂多样，以低山、丘陵为主，地势北高南低。

建模数据：从广东省第5—9次连清数据中，选择马尾松（Pinus massoniana Lamb.）、湿地松（Pinus elliottii Engelm.）、杉木（Cunninghamia lanceolata (Lamb.) Hook.）3个针叶树种，筛选5次清查中皆存在且保留木大于5株的样地，共155块（马尾松林91块、湿地松林16块、杉木林48块），样地分布见图1。固定样地为0.0667 hm²，样木调查因子包括样木号、胸径，单株材积按广东省地方一元材积表^[21]计算获得，单株生物量按行业标准生物量模型^[22-24]获得。为避免株数差异造成样地蓄积量和生物量的差异，本研究中样地蓄积量和生物量为样地内5期保留木的蓄积量和生物量算数平均值，统计量如表1所示。

Figure 1.  Distribution map of fixed sample plots of modeling data

样地胸径采用林分平均胸径：

式中，$ {D}_{g} $为林分平均胸径，n为样木株数，d_i为样木胸径。

验证数据：选取2002年达到起测胸径、存在于4次连清中，且保留木大于5株的样地，共126块，其中马尾松59块、湿地松22块、杉木45块。各统计量如表2所示。

采用Mitscherlich理论生长方程，拟合5期林分平均胸径和相应时间间隔之间的关系；以Richards理论生长方程，拟合样地平均蓄积量和生物量与相应时间间隔之间的关系。

式中：$ {D}_{gij} $、$ {V}_{ij} $、$ {BIO}_{ij} $分别为第i个样地第j期的样地平均胸径、蓄积量和生物量，$ {age}_{i} $是第$ i $个样地的初期林龄，$ {ad}_{j} $为第$ j $期与初期（1997年）的时间间隔，$ a1、a2、a3、b1、b2、b3、c2、c3 $为参数。

基于上述胸径、蓄积量和生物量3个生长方程，分别以体现生长潜力的参数a分级和与生长速度有关的参数b分级、考虑竞争与否、分步或联合建模等8种组合，构建模型系（表3）。

对于林分水平模型来说, 与树冠有关的竞争指数能较好地反映树木间的竞争关系，因此选用基于交角且对象木与竞争木距离小于10 m的竞争指数$ ju10 $，详见惠刚盈等^[25]。由于林木竞争对单木胸径的生长具有较大的影响，但对林分平均胸径的影响并不算太大，为解释林木竞争对林分生长的影响，故将竞争指数引入了胸径生长方程：

式中：$ {b}_{0} $为竞争指数参数，$ {ju10}_{ij} $为第i个样地第j期的竞争指数（样地竞争指数为样地内所有保留木竞争指数的算数平均值），其余同上。

为了兼顾区域尺度上立地质量的多样性以提高模型在大尺度范围内的适用性，对参数a或b分级代替立地质量的差异，以参数a分级为例：

将式（2）作为立地质量划分等级的分级方程，已知林分平均胸径$ {D}_{gij} $和时间间隔$ {ad}_{j} $，以样地号和树种号构成的唯一编号作为参数a和age的哑变量，树种作为参数b的哑变量，拟合方程（6）。将参数${a1}_{i}$从低到高排序，并按照有序样本聚类分为6类。

式中： $ {a1}_{i}(i=\mathrm{1,2},\dots ,155) $是第i个样地的参数，$ {b1}_{k}(k=\mathrm{1,2},3 $时分别代表马尾松、湿地松和杉木)是第$ k $个树种的参数，其余同上。

以Model 5为例，基于式（2），以参数a区别树种，以参数b区别立地质量差异，以唯一编号作为初期林龄age（未知）的哑变量，根据5期林分平均胸径$ {D}_{gij} $和时间间隔$ {ad}_{j} $之间的关系来估计样地初期林龄$ {age}_{i} $，见式（7）。将$ {age}_{i} $作为式（3）和（4）的自变量，以参数c区别树种，其余参数处理同上，分别拟合样地平均蓄积量$ {V}_{ij} $和生物量$ {BIO}_{ij} $与林龄$ {age}_{i} $的关系，即式（8）和（9）。

式（7）～（9）中：$ {a1}_{k}、{a2}_{k}、{a3}_{k}\left(k=\mathrm{1,2},3\right) $分别为胸径、蓄积量和生物量生长方程中的第k个树种的参数，$ {b1}_{m}、{b2}_{m}、{b3}_{m}\left(m=\mathrm{1,2},\cdots,6\right) $分别为胸径、蓄积量和生物量生长方程中的第m个立地质量等级的参数，$ {c2}_{k}、{c3}_{k}\left(k=\mathrm{1,2},3\right) $分别为蓄积量和生物量生长方程中的第k个树种的参数，其余同上。

为了比较建模方法对生长模型拟合优度和估测精度的差异，除分步建模法外，本研究还采用了联合建模法，以Model 7为例。由于林龄未知，且3个方程中年龄参数$ {age}_{i} $是一致的，基于式（2）～（4），参数处理同分步建模法，根据林分平均胸径、蓄积量和生物量与时间间隔之间的关系，采用非线性误差变量联立方程组法同时拟合3个方程，即式（10）。

考虑各模型参数不同，为了检验8种模型系的拟合效果，采用以下4项模型评价指标：确定系数$ {R}^{2} $、估计值标准误差SEE、平均预估误差MPE、总体相对误差TRE，公式详见文献^[26]。

（1）确定立地质量等级和林龄（以参数b分级为例）

根据胸径生长方程反推林龄计算公式，如下：

式中，$ pag{e}_{ij} $为第i个样地第j期的样地林龄，当不引入竞争指数时$ {b}_{0}=0 $，其余同上。

由于验证样本所在样地的立地质量等级m和初期林龄皆未知，需先确定立地质量等级，再由式（11）求出林龄，才可估计蓄积量和生物量。当立地质量分级时，至少需要两期胸径才能确定样地的立地质量等级，薛春泉等^[27]发现长周期确定的等级比短周期稳定，因此本研究采用期初（2002年）和期末（2017年）2期胸径来确定样地的立地质量等级，确定方法如下：

式中，$ pag{e}_{i1}\mathrm{和}pag{e}_{i4} $为第i个样地期初和期末的样地林龄，该方法的直观解释就是使推算的2002年和2017年林龄间隔与实际林龄间隔（15年）相差最小的立地质量等级，即为该样地的立地质量等级。

根据式（12）确定样地的立地质量等级后可知参数$ {b}_{m} $，再由式（11）计算初期（2002年）林龄age，则2007年、2012年、2017年时样地的林龄分别由初期林龄相应加上5、10、15年。

（2）模型估测效果评价指标

由蓄积量和生物量生长方程估计样地平均蓄积量和生物量，样地总蓄积量和总生物量可由平均蓄积量和生物量乘以株数获得，按式（13）～（14）计算各树种4个时期下区域尺度蓄积量和生物量估计的总相对误差EV和EBIO。

式中：n为样地数，$ {N}_{i} $为第i个样地的保留木株数，$ P{V}_{i} $、$ {V}_{i} $分别为第i个样地平均蓄积量的估计值和真实值，$ PBI{O}_{i} $、$ BI{O}_{i} $分别为第i个样地平均生物量的估计值和真实值。

连年生长量函数$ Z\left(t\right) $为总生长过程曲线$ {y}\left({t}\right) $的一阶导数，为树木在某年的实际生长速度，即连年生长量随年龄t变化的函数。Richards生长方程的连年生长量方程为：

平均生长量函数$ \theta \left(t\right) $用来表示树木在某一时刻的平均生长速度，总生长过程曲线$ y\left(t\right) $除以年龄t则可得平均生长量依年龄变化的函数。Richards生长方程的平均生长量方程为：

表4列出了8种模型系的拟合评价指标。整体上看，8种模型系的3个生长方程拟合确定系数$ {R}^{2} $皆达0.980以上，总相对误差TRE在±0.5%以内，平均预估误差MPE在±1%以内，拟合效果良好，但各模型间的4项指标差异不大。当控制建模方法和立地质量分级的参数相同、仅比较竞争影响时，引入竞争指数对模型拟合优度的提升没有明显的作用，如Model 5和Model 6拟合指标间差异甚微、Model 2较Model 1的蓄积量和生物量生长方程的拟合效果甚至出现了一定程度的下降；类似地，可以发现参数b分级模型的拟合优度较参数a分级模型稍好；当控制竞争影响和分级参数相同时，相对于分步建模，采用含度量误差的联立方程组时，蓄积量和生物量生长方程的拟合效果有所提高，但是胸径生长方程的拟合效果稍差。

8种蓄积量和生物量生长方程的模型验证指标分别如图2和图3所示。总体而言，各模型的预测效果都较好，从图2可知，区域尺度蓄积量估计误差最大的仅有17.62%，说明蓄积量估计精度可达82.38%以上；从图3可知，区域尺度生物量估计误差最大的仅有12.51%，则生物量估计精度可达87.49%以上。8个模型中以基于参数b分级的分步模型（即Model 5和Model 6）预测效果最优，其对3个树种4期蓄积量和生物量的估计误差较其余模型稍小，且在不同树种和不同时期的表现都较为稳定，说明模型适用性较好。以Model 5为例，3个树种中对杉木的估计效果稍好，估计误差最大为−5.25%，而估计马尾松和湿地松时误差最大为10.00%和10.36%。8个模型对3个针叶树种的4期区域尺度蓄积量和生物量估计的误差，大都表现出中期误差高于首末2期，初期估计精度较高，但明显出现低估的现象，伴随着样地林龄的增大，中期出现高估且精度稍低，但到末期（2017年）时估计误差降低，精度有所提高。

Figure 2.  Regional-scale estimation errors of eight volume growth models

Figure 3.  Regional-scale estimation errors of eight biomass growth models

综合8个模型的拟合效果（表4）和验证效果（图2、图3），Model 5最优，Model 6次之。表5和表6列出了基于参数b分级的分步模型的模型参数估计值。由表5可知，马尾松胸径、蓄积量和生物量的生长潜力均为最大，杉木最小，其中马尾松蓄积量和生物量的生长极限值可达湿地松和杉木的2倍以上。由表6可知，立地质量等级越高，树木的生长速度越快，但是各等级下的林分蓄积量生长速度均为最快，其次为生物量，胸径生长速度最慢。

基于参数b分级不考虑竞争的分步建模（即Model 5）最优，为节约篇幅，选择样地立地质量等级占比较多的等级2、3和4，绘制Model 5下蓄积量、生物量的连年生长量和平均生长量的关系曲线图（图4）。蓄积量和生物量的连年生长量、平均生长量均显示出立地质量越好，生长量能达到的极限值越大，且所需年限更短，各相邻等级间的差异随着等级的增加而逐渐加大。对于不同的树种来说，同一立地质量等级下的生长量均为马尾松>湿地松>杉木，生长量达到最大时的林龄为杉木>湿地松>马尾松，各树种间生长量差异随着等级的增加而逐渐加大。立地质量等级越高，平均生长量和连年生长量曲线相交时的样地年龄越小。

Figure 4.  Annual increment and average increment from class 2 to class 4

处理多树种问题是建立林分模型的一个难点，通常是将组成树种进行分组或合并^[28]。曾伟生等^[2]建立了东北林区林分水平的三储量模型，筛选出10种主要森林类型并引入哑变量体现差异后，蓄积量模型、生物量模型和碳储量模型的确定系数分别从0.945、0.805、0.839提高到0.959、0.949、0.951。本研究对存在多个组成树种的样地按优势树种分属不同的样地，并构建哑变量区别广东省的3个主要针叶树种，从表5和图4可知，各树种间的生长潜力和生长量有较大的差异，表现为马尾松>湿地松>杉木，其中马尾松远超2倍左右，这是因为建模数据（表1）中马尾松保留木大多为天然林，而杉木的起源大多为人工林，Zeng等^[29]研究结果表明同一直径的单木地上生物量估计值在不同起源上存在一定差异，天然林的生物量大于人工林。此外，在同一立地质量下，不同树种间生长量的差异主要是由于树种特性的差异引起的，不同树种对同一区域的生长适应程度并不一致，马尾松可能更适宜广东的立地条件。各树种间区域尺度蓄积量和生物量的预测精度也有较大的差异，对杉木的估计效果较马尾松和湿地松稍好，可能原因是不同树种的建模样本量和检验样本量不一致，另一个可能的原因是不同的树种有不同最优的基础生长模型，如吴宏炜等^[30]研究发现，Schumacher模型是湿地松林分蓄积量模型的最优基础模型。对于估计误差，大都表现出中期高于首末两期的现象，一方面可能是由于采用首末两期胸径数据来确定验证样本所在样地的立地质量等级，相较于中期，推算的各样地立地质量等级更贴合期初和期末时的实际情况，因此用于林龄确定、估计蓄积量和生物量时，3个生长方程中反映立地质量差异的参数更为稳定，模型的预测效果会更好。另一方面随着林龄增大，蓄积量和生物量逐渐积累，公式（13）和（14）分母变大，这可能使后期的估计误差降低，估计效果更好。

非线性独立回归估计法的估计效果优于含度量误差变量联立方程组法，因为分步建模时，3个方程都是独立拟合给出参数的最优估计，而采用联立方程组法联立求解时，由于方程中林龄一致的约束，为保证方程之间的联系，得到更优化的参数估计，各方程较独立拟合时都有了不同的变化。从拟合评价指标（表4）可以看出，对比分步建模法，采用联立方程组法可以提高蓄积量和生物量生长方程的拟合效果，但是却降低了胸径生长方程的拟合优度，导致林分年龄的估计精度有所下降，使用验证数据集检验模型的适用性时，由于林分年龄未知，确定首期林龄成为关键的一步，所以胸径生长方程的拟合效果直接影响了验证时的估计精度。由于林分大部分为混交林，具有年龄不一、树种结构关系复杂等特点，因此探究林木间竞争作用对林分结构和生长变化的影响至关重要。一般而言，模型中的参数个数越多，拟合效果会更优，但是引入竞争指数后对于模型的拟合效果（表4）和估计精度（图2、3）没有明显的提高。虽然从理论上讲包含树木空间位置的竞争指数对预测林木生长更有效, 但Biging等^[31]研究发现与距离无关的竞争指数比与距离有关的竞争指数预测单木生长更好，臧灏等^[32]探究人工杉木林胸径生长对竞争的响应时也发现对比于与距离有关的单木竞争指标，采用与距离无关的单木竞争指标的生长模型估计精度更高。竞争指数表现的预测能力依赖于树种、立地条件等因素，这些因素的变化会导致不同类型的竞争指数表现不一致，且不具有决定性^[33-34]。本研究引入竞争指数对林分模型的适用性没有明显影响，可能与选取的竞争指数有关，也有可能和立地的限制有关，因为在较大范围内密度对林分蓄积量、生物量的影响不大^[35]。

模型的应用范围和模型的预测精度相关联，当模型应用于大尺度的森林资源测算时，若忽略立地质量的差异必然导致模型估计结果的不准确。立地质量是指在某一立地上既定森林或其他植被类型的生产潜力^[35]，既定立地上的不同树种和相同的树种在不同的立地上的立地质量都是有差异的。本研究为体现广东省样地间的立地质量多样性，分别对与生长潜力有关的参数a和与生长速度有关的参数b进行分级，结果表明以参数b分级体现立地质量的差异优于参数a分级。薛春泉等^[26]研究发现生长速度分级的胸径生长方程与地上生物量的异速生长方程构成的模型系，不仅拟合效果达到了0.95以上，区域尺度生物量的估计误差也下降到了8%以内，这与本研究参数b分级模型预测效果较好具有一致性。另外，林分蓄积量和生物量的生长量随着立地质量的提高而增加，说明立地质量等级越高，森林蓄积量和生物量密度越大。立地质量与森林生产力相关，不同的立地质量代表着不同的生产力，立地质量越高，林分的生长速率越快，达到最大生长量所需的年限也更短。赵菡等^[36]用优势木树高等级代替立地等级估计江西省不同立地质量等级的马尾松林生物量，研究发现不同立地等级的地上生物量均值估计结果随着立地质量等级的升高而增大，这与本研究结果一致。

本研究以广东省3种针叶优势树种为研究对象，利用1997—2017年的5期森林资源连续清查数据，基于林分平均胸径、蓄积量和生物量与相应调查时间间隔的关系建立林分生长模型，实现了对广东省蓄积量和生物量的准确估计。本研究结论如下：（1）基于与生长速度有关的参数b分级不考虑竞争的分步模型最优，引入竞争指数与否差异不大。（2）最优模型的区域尺度蓄积量和生物量估计误差最大分别为10.36%和10.22%，模型的适用性较好，4期的估计误差表现为中期高于首末两期。（3）不同树种中以马尾松的生长潜力最大，杉木最小，但是模型对杉木的估计效果较马尾松和湿地松稍好。（4）立地质量等级越高，生长量极值越大，达到最高峰所需年限也更短；同一立地质量下马尾松的生长量最大，其次为湿地松，杉木最小。