20世纪80年代末以来，全球气候发生显著变化，引发了一系列生态问题；在这种背景下，森林的碳汇作用越来越受到重视(Nabuurs et al., 2003；Richards et al., 2004；张颖等，2015；Zou et al., 2015)。作为森林碳汇计量的基础，迄今为止，国内外建立了数以千计的生物量模型(Ter-Mikaelian et al., 1997；Zianis et al., 2005；曾伟生，2011；Dong et al., 2014；Ozdemir et al., 2019)，我国还发布了13个主要树种的立木生物量模型和碳计量参数的行业标准(国家林业局，2014a；2016a)。这些模型基本均为异速生长方程(相对生长方程)，是解决某一时间节点树木生物量与自身其他测定因子(直径、树高)之间的关系模型，未涉及时间和环境因子(Rubner，1883；罗红等，2011；汪珍川等，2015；Cao et al., 2019)，模型精度虽高，但模型的应用范围与建模样本大小和抽样范围有关，在估算森林碳汇(源)时，至少需要2次调查，且只能获得调查期间的碳汇变动情况。理论生长方程作为描述一种有机体或一个种群大小随年龄变化的模型，可以反映某些生物生长的规律性(Zeide，1993；Rohner et al., 2018)，其特点是逻辑性强，参数具有生物学意义解释，并可从理论上对尚未观察的事实进行预测，属于机理模型范畴(张建国等，2003；2004；Sharma et al., 2004)。生物量理论生长方程，以年龄为自变量，可以通过1次调查获得一个时间系列多个节点的生物量估算结果，为森林生物量动态评估提供了极大便利，如果时间系列以1年为间隔，则可获得碳储量的年增长量(年碳汇量)；但是由于时间的一维性，模型数据获取比较困难，大范围应用时，空间异质性和复杂性常导致模型精度不理想，因此生物量生长方程研究少有报道，影响了模型在实践中的应用。

木荷(Schima superba)为山茶科(Theaceae)木荷属(Schima)大乔木，是一种优良的绿化、珍贵用材树种，同时又具备较好的耐火、抗火和难燃特性(郑万钧，1997)。近年来，广东省人工造林大量营造木荷，第7次清查时全省木荷样木为2 249株，第8次清查时增长到7 465株，第9次清查时达10 464株。本研究利用90株木荷实测生物量、胸径和树高等因子，对其中所有16 cm径阶以上的大树和12 cm径阶中采用均匀抽样方法抽取的树木共40株进行树干解析，获得相应基础数据，构建独立的生物量异速生长方程和兼容的生物量理论生长方程，分别在单木、样地和区域尺度上对广东省第7—9次清查中的木荷生物量进行估算和动态预测，以期为广东省木荷林碳汇计量提供模型和方法，同时为其他树种碳汇计量提供可借鉴的方法学支持。

1 数据 1.1 建模数据

木荷生物量建模数据源自广东省21个地、市，采用全株收获法于2013年获得，共90株。起测径阶2 cm，其中，胸径小于10 cm按2 cm径阶整化，胸径10~20 cm按4 cm径阶整化，胸径20 cm以上按6 cm径阶整化，最大径阶38 cm，共10个径阶，每个径阶选取6株样木(国家林业局，2014b)；另外30株按2012年第8次广东省连续清查固定样地中木荷的实际数量和大小在广东全省进行分配。各径阶样木的分配状况如下：2 cm径阶6株，4 cm径阶6株，6 cm径阶14株，8 cm径阶16株，12 cm径阶14株，16 cm径阶8株，20 cm径阶8株，26 cm径阶6株，32 cm径阶6株，38 cm径阶6株。根据建模要求，按胸径和树高2项主要因子进行控制，分不同径阶和树高选取标准样木。样木选择没有发生断梢、分叉且生长正常的树木，采用单株伐倒法取样。全部样木实测胸径、树高和冠幅等因子，分别称量树枝、树叶、干材和干皮鲜质量，然后分别抽取样品带回实验室，85 ℃恒温下烘干至恒重，根据各样品鲜质量和干质量的关系计算各组分含水率，再推算出样木各部分干物质量，各部分相加得到地上总生物量(国家林业局，2014b)。

树干解析时，0号盘取自0.3 m，1号盘取自1.3 m(胸高处)，2号盘取自3.3 m，继续以2 m为长度分段，以此类推直至尾端长度不足2 m但大于1 m则加以1 m分段，各取一圆盘，树干解析的取样数量根据伐倒木长度确定，依据树干北向在每个圆盘上进行北向标记(孟宪宇，2006)。内业对经刨光和标记南北向的圆盘样品，采用砂轮将每个圆盘工作面磨平后扫描，使用WinDENDRO软件进行年轮数计数和直径生长量测量。由于0号盘受根部扩张影响大，断面积形状与圆形相差较大，直径年生长量不稳定，故本研究利用1号盘获得的解析木年龄(胸高处年龄)和胸径生长量建立生长方程。解析木生长过程中各年龄的地上生物量，用该年龄时区分求积获得的材积乘以本株树木地上总生物量和材积之比推算，其余没有进行树干解析的树木胸高处年龄，用生长锥在胸高处钻取木芯，通过查数年轮获得。建模样本和解析木样本统计信息见表 1。

1.2 预测数据

选择广东省最近3期(2007、2012和2017年)森林资源连续清查中固定样地均存在共计1 053株木荷作为预测数据，样木分布在119个样地中，样地大多为阔叶混交林，只考虑木荷，平均每块样地约9株，每块样地分布有1~70株木荷样木。表 2给出了样木在单木、样地(林分)和区域(省级)尺度上的胸径统计量，其中单木平均值为算术平均值，样地尺度和区域尺度平均值为平方平均值。

2 方法 2.1 模型构造

以胸径为自变量，构建生物量异速生长方程[式(1)]；选择Richards方程建立胸径和胸高处年龄的理论生长模型，通过联立方程组构建理论生长方程与异速生长方程相结合的生物量估算模型系(Richards，1959)[式(2)]：

$ \mathrm{AGB}=a_{0} D^{b_{0}}; $

(1)

$ \left\{\begin{array}{l}{D=a\left(1-\mathrm{e}^{-b t}\right)^{c}} \\ {\mathrm{AGB}=a_{0} D^{b_{0}}}\end{array}\right.。$

(2)

式中：AGB为地上生物量(kg)；D为胸径(cm)，为了统一，利用树干解析资料中的树皮系数，方程中的胸径均已转换为带皮胸径；t为胸高处年龄(a)；a、b、c、a₀和b₀为参数。

由于抽样范围较广，不同样木之间生长速度差异很大，所以样木按生长速度(快、中、慢)分为3个等级，式(2)可变形为：

$ \left\{\begin{array}{l}{D=a\left(1-\mathrm{e}^{-b_{i} t}\right)^{c}} \\ {\mathrm{AGB}=a_{0} D^{b_{0}}}\end{array}, i=1, 2, 3\right.。$

(3)

式中：i为生长速度等级，1表示快，2表示中等，3表示慢，每株样木(非解析木)只能属于一个等级。

对每株解析木，约束整个生长过程具有相同的等级，具体方法和分级过程详见文献(李海奎等，2011；Li et al., 2013)。以参数b区分生长速度等级，主要原因有3点：1)在Richards模型中，b与生长速度有关；2) b变异系数最大；3)如果用形状参数c反映生长速度可能会比较小(前期数据预处理发现)。

2.2 年龄确定

由式(2)中的第1个方程可知：

$ t=\frac{-\ln \left[1-(D / a)^{1 / c}\right]}{b}。$

(4)

当生长速度分级时，式(4)变为：

$ t_{i}=\frac{-\ln \left[1-(D / a)^{1 / c}\right]}{b_{i}}, i=1, 2, 3。$

(5)

式中：t_i为生长速度等级为i时的推算年龄。

由胸径推算年龄，如果树木生长速度没有分级，只要知道胸径，就可以通过式(4)直接获得年龄；当分级时，则需要至少2期(间隔期明确)胸径才能获得树木生长速度等级，此时再用式(4)获得年龄。生长速度等级确定方法如下：

$\begin{array}{c} \Delta {t_i} = \frac{{ - \ln \left[ {1 - {{\left({{D_1}/a} \right)}^{1/c}}} \right]}}{{{b_i}}} - \frac{{ - \ln \left[ {1 - {{\left({{D_0}/a} \right)}^{1/c}}} \right]}}{{{b_i}}}, \\ \min \left| {\Delta {t_i} - \Delta t} \right|。\end{array}$

(6)

式中：Δt_i为生长速度等级为i时的推算年龄间隔；D₁为当期胸径；D₀为前期胸径；Δt为实际年龄间隔。式(6)的直观解释就是使推算年龄间隔与实际年龄间隔最小的生长速度等级。

如果只有1期胸径数据，同时知道年龄，对每株树木，已知参数a和c，按式(4)推算出参数b的估计值$\hat{b} $，按$\min \left|\hat{b}-b_{i}\right| $也可获得树木生长速度等级。

2.3 生物量预测

由于预测数据只有胸径没有年龄，故在树木生长速度分级情况下，需要2期数据才能确定1株树木的胸高处年龄。本研究将由前期胸径推算后期生物量的方法称为预测估计，如利用2007和2012年2期数据确定树木生长速度分级，然后估计得2017年生物量。具体方法如下：

1) 由式(6)获得生长速度分级，假设为2，如果2007年的年龄为t₀，则2012年的年龄为t₀+5；

2) 2017年的年龄为t₀+10；

3) 2017年的胸径为a[1-e^{-b₂(t₀+10})]^c；

4) 利用式(3)的第2个方程获得2017年和等级为2的生物量。

类似地，本研究将由后期胸径推算前期生物量的方法称为回溯估计，前后2期胸径估计中间年度生物量的方法称为插值估计。具体步骤和预测估计类似，只是年龄的确定不同。

某一期总生物量的估算按以下3种方式进行：

$ \mathrm{PBio}_{\mathrm{T}}=\sum\limits_{i=1}^{n} \mathrm{AGB}_{i}; $

(7)

${\rm{PBi}}{{\rm{o}}_{\rm{T}}} = \sum\limits_{i = 1}^m {{{\overline {{\rm{AGB}}} }_{si}}} {n_i};$

(8)

$ \mathrm{PBio}_{\mathrm{T}}=\overline{\mathrm{AGB}}_{r i} n。$

(9)

式中：PBio_T为估算的某一期总生物量；AGB_i为第i株树的生物量；n为样木总株数；AGB_si为第i块样地平均木的生物量；n_i为第i块样地样木株数；m为样地数；AGB_ri为区域内所有样木平均木的生物量。

式(7)—(9)分别为基于单木、林分和区域水平的总生物量估算方法，其中平均木的计算均采用平方平均值。

2.4 评价指标

异速生长方程[式(1)]采用非线性回归拟合，理论生长方程与异速生长方程相结合的生物量估算模型系[式(2)和式(3)]采用非线性联立方程组进行联合估计，拟合采用权函数消除异方差(Tang et al., 2001；2002)。为了避免利用反函数[式(4)和式(5)]推算年龄时可能出现的超出定义域现象，考虑木荷的生物学特性，本研究固定参数a为100。所有拟合和计算在Forstat软件(唐守正等，2009)和R软件(Henningsen et al., 2007)上完成，使用确定系数(R²)和均方根误差(RMSE)作为模型精度评价指标(曾伟生等，2011)：

${R^2} = 1 - \sum\limits_{i = 1}^n {{{\left({{y_i} - {{\hat y}_i}} \right)}^2}} /\sum\limits_{i = 1}^n {{{\left({{y_i} - \hat y} \right)}^2}} ;$

(10)

$ \mathrm{RMSE}=\sqrt{\sum\limits_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2} / n}。$

(11)

式中：y_i为第i个样本实测值；$\hat{y}_{i} $为第i个样本估计值；y为平均实测值；n为样本数。

生物量预测时，以存量估计误差(P_B)和增量估计误差(P_Δ)作为评价指标：

$ P_{\mathrm{B}}=\frac{\left(\mathrm{PBio}_{\mathrm{T}}-\mathrm{Bio}_{\mathrm{T}}\right)}{\mathrm{Bio}_{\mathrm{T}}} \times 100 \%; $

(12)

$ P_{\Delta}=\frac{\left(\Delta \mathrm{PBio}_{\mathrm{T}}-\Delta \mathrm{Bio}_{\mathrm{T}}\right)}{\Delta \mathrm{Bio}_{\mathrm{T}}} \times 100 \%。$

(13)

式中：PBio_T为某期总生物量估计值；Bio_T为某期总生物量真值，通常情况下是未知的，以异速生长方程[式(1)]单木尺度的生物量合计代替；ΔPBio_T为估计的2期总生物量增加值；ΔBio_T为2期总生物量增加值的真值。

为了比较拟合效果，采用林业行业标准中木荷异速生长方程参数估算的总生物量作为对照(国家林业局，2016b)。

3 结果与分析

表 3所示为模型拟合的参数及评价指标，可以看出，独立的生物量异速生长方程中，胸径对地上生物量的解释率达0.88以上；生长速度不分级的模型系，虽然年龄对胸径的解释率达0.78以上，但均方根误差仍超过4 cm，造成预测胸径对地上生物量的解释率只有0.53，均方根误差高达80 kg以上；生长速度分级的模型系，理论生长方程中年龄对胸径的解释率达0.95以上，比不分级提高0.166 3，均方根误差下降到1.97 cm，降低2.16 cm以上，预测胸径对地上生物量的解释率提高到近0.82，提高近0.30，均方根误差下降到51 kg左右，下降30 kg以上。

表 4列出了生长速度分级和不分级情况下，2007、2012和2017年3期木荷总生物量的估计误差。不分级时，生物量估计值由1期数据即可以获得，如-5.83%就是由2012年的调查数据在时间上减去5年，在单木水平上获得的2007年生物量估计误差。分级时，首先要用2期数据获得速度分级，然后再由1期数据获得不同时期的估计值，如-0.41%就是由2012和2017年2期数据获得速度分级后，在单木水平上由2012年数据推算的2007年生物量估计误差。表 4显示，不分级时，各期生物量的估计误差变动很大，从-46.31%到77.45%，在不同时期的估计误差中，2012年估计2007年的误差最小，其次是2017年估计2012年，2个时期均低估生物量10%左右；2007年估计2017年的误差最大，3种尺度均高估，幅度在68%以上，其次是2007年估计2012年，均高估50%以上。前期估计后期，生物量均高估，后期估计前期，生物量均低估，高估幅度大于低估幅度。间隔期越长，估计误差越大，如单木水平上，2007年估计2012年的误差为51.08%，估计2017年的误差为68.82%。在估计尺度上，除了单木水平上2017年估计2007年没有结果外，均呈相同规律，即单木误差小于林分误差、林分误差小于区域误差，且水平间幅度基本相当，单木水平和区域尺度间的误差不大于10%。生长速度分级时，各期生物量的估计误差显著减小，除了在区域尺度上2017年估计2007年的误差为-16.13%外，其余不同时期的估计误差均不超过10%。前期估计后期，生物量均高估，后期估计前期，生物量均低估，与不分级情况一致，但高估幅度小于低估幅度，与不分级情况相反。间隔期长短对估计误差影响不大，特别是在单木水平上几乎无差别。在估计尺度上，单木水平误差基本不超过3%，单木误差小于林分误差、林分误差小于区域误差，前期估计后期，单木水平和区域尺度间的误差小于5%；后期估计前期，除2017年估计2007年外，单木水平和区域尺度间的误差均小于8%。

表 5所示为不同间隔期木荷生物量增量的估计误差。不分级时，各间隔期生物量增量的估计误差均为正值，说明估计值普遍偏大，从32.57%到115.45%，变化幅度较大。后期估计前期的误差较小，如在林分水平上，由2017年估计2007和2012年间的生物量增量，误差为32.57%(P_Δ2007-2012)，前期估计后期的误差较大(P_Δ2012-2017)，在3种尺度上均大于100%，而10年间隔的估计误差(P_Δ2007-2017)处于中等，主要是因为由2012年估计2007和2017年的间隔都是5年，各期生物量总量较为准确，所以增量估计误差不是最大。在估计尺度上，单木误差小于林分误差、林分误差小于区域误差，且水平间幅度基本相当，单木水平和区域尺度间的误差介于10%~15%之间。生长速度分级时，各间隔期生物量增量的估计误差均显著减小，除了在林分水平和区域尺度上前期估计后期的误差(P_Δ2012-2017)大于10%外，其余估计误差均小于10%。前期估计后期及由中间时间估计前后2期，生物量均高估，后期估计前期，在3种尺度上估计误差符号不一，但绝对误差量区域最小，单木最大，与其余2种估计类型不一致。与真值接近程度上，3种间隔的增量估计在单木水平最为稳定，均不超过10%，其次是林分水平，区域尺度普遍高估。从单木到林分再到区域，增量估计值不断增加，但3种间隔增量估计误差的差异比较稳定，均在8%左右。

4 讨论

生长速度分级，可以显著改善模型拟合精度(表 3)，极大降低预测数据估计误差，对生物量存量的估计(表 4)及对增量的预测都有很好效果(表 5)。如何确定生长速度等级是关键问题，在应用范围较广、立地条件存在较大差异时尤为重要。本研究使用森林资源连续清查3期固定样地样木数据，任意2期用于确定胸径生长速度等级，1期用于预测，此时预测时期可能在2期之前、之后或中间，预测数据起点可以是确定胸径生长速度等级的2期中任意一期。在确定生长速度等级后进行增量预测时，间隔期的两端均为预测值，另一期为预测数据起点，如利用2007和2012年数据确定样木生长速度等级，以2007年为预测起始点，预测2012和2017年的总生物量，二者之差即为2012—2017年的生物量增量。表 6列出了不同间隔期确定的生长速度等级一致性比较，可以看出，单木水平和样地水平，前期(2007—2012年，表 6中下标1)和后期(2012—2017年，表 6中下标2)确定的生长速度等级相同比例均在70%以上，5年间隔期(表 6中下标1和2)和10年间隔期(表 6中下标3)确定的生长速度等级相同比例均在80%以上，这说明，周期较长的等级确定比短周期更为稳定。3期确定的生长速度等级完全相同比例也平均达70%以上，单木水平和样地水平的比例基本相同，表明尺度对生长速度等级的确定基本没有影响。

碳汇造林时，由于是人工林，可以准确知道树木年龄，而式(4)中年龄为胸高处年龄，需要进行适当转换，将树木实际年龄转换为胸高处年龄，这样只要在某个年度测定1期树木胸径，就可以推算树木胸径的生长速度参数，与各等级生长速度参数比较获得具体的生长速度等级，然后就可以预测未来任何年度的生物量，从而获得任意年度间的碳汇量。根据本研究的解析木数据，在生长速度等级为1、2或3时，树木达到1.3 m(胸高)需要的时间约为1、2或3年，此时胸高处年龄则分别为树木实际年龄减去1、2或3年。

决定估计误差的另一个重要因子是生物量真值。本研究使用的生物量真值由独立的异速生长方程获得，数据集与建立异速生长方程和理论生长方程模型系的数据相同。一般来讲，独立的异速生物量模型精度较高，主要是因为模型的胸径、树高等变量均为生长的结果，包含立地因子，模型的拟合没有条件限制(袁亚湘等，1997)，本研究中只有胸径1个变量的模型，就解释了生物量变化的近90%。分级后，胸径和年龄的关系也显著改善，决定系数达0.95以上，从而预测胸径和地上生物量的关系也达0.80以上，与独立异速生物量模型的决定系数仅相差0.06(表 3)。木荷作为我国南方主要的常绿阔叶树种之一，国家林业局于2016年发布了其立木生物量模型和碳计量参数的行业标准(LY/T 2660—2016)，建模样木共150株，涉及云南、福建、广西和广东等7个省(区)，但广东的样木只有10株，一元模型为异速生长方程，决定系数为0.956 4。表 7给出了不同模型估计的3期总生物量对比，可以看出，即使模型决定系数相差较大，但估计的总生物量非常相近，如在单木尺度上，2007年总生物量3种模型估计结果基本相等，2017年也仅相差5%左右。在尺度转换上，3种模型具有相同规律，即单木水平较大，区域尺度较小，林分水平位于中间，但理论生长方程和异速生长方程相结合的模型系最为稳定，3期均不超过2%，行业标准变动最大，3期变动幅度在3.5%~5.0%之间，独立模型处于中间，在2%左右变动。分析各时期总生物量发现，随着年龄增大，区域尺度和单木水平的总生物量差别不断减小，3种模型表现出相同规律性，使模型间的差别越来越小，如在2017年，行业标准与联合估计的模型系，单木水平的差异(以联合估计的生物量为基数)为5.7%，区域尺度下降到3.3%。

5 结论

1) 生长速度分级的胸径生长方程与地上生物量的异速生长方程构成的模型系，可以极大提高模型精度，减小预测估计误差。

2) 生长速度不分级时，仅利用胸径或年龄数据，就可以预测未来的生物量动态，推算以前的生物量；生长速度分级时，利用2期胸径数据或1期胸径和年龄数据，能够更为准确地预测未来的生物量动态。

3) 不同尺度上，对生物量的动态预测误差不同，单木水平误差最小，区域尺度误差最大，但二者差异不大，对存量的估计误差不超过6%，对增量的估计误差不超过8%。