森林是陆地生态系统的主体，在维护全球气候系统、调节全球碳平衡、减缓大气温室气体浓度上升等方面具有不可替代的作用，世界各国越来越重视对森林生物量的估计、监测和评价。以往在大尺度监测森林生物量时，更多关注的是总生物量及地上和地下生物量，但是在研究森林生态系统生产力、过程或机制模型、能量和营养物质分布时，则需将总生物量分解成不同的组成部分(Waring et al., 1998)，如将地上生物量分解成树干和树冠、树干分解成干材和干皮、树冠分解成树枝和树叶等；而且，在进行森林生态系统分析时，总物质必须满足分解与合成相一致(分量加起来等于总量)。因此，构建各分项单木生物量相容性方程系统一直为世界生态学家和林学家所重视。到目前为止，国内外已提出了不同方法构建生物量相容性(或可加性)方程系统，其中最常见的有非线性似然无关回归方法(Parresol，2001；Bi et al., 2004；符利勇等，2014；Dong et al., 2015)、比例平差法(唐守正等，2000；符利勇等，2014)和线性或非线性联合估计方法等(唐守正等，2000；2015；符利勇等，2014；Dong et al., 2015；Zou et al., 2015)。

马尾松(Pinus massoniana)是我国亚热带地区分布面积最广的针叶树种，生长快、适应性强、耐干旱瘠薄，为重要的工业原料造林树种之一，广泛用于制浆造纸、建筑、松香等产业，在速生丰产用材林基地建设中具有重要地位和作用(丁贵杰等，2006)。迄今，关于马尾松生物量建模国内学者已做了大量研究，符利勇等(2011a)基于混合效应模型分析地域对马尾松生物量的影响，结果发现混合效应模型不仅可提高模型预测精度，而且模型参数还能有效定量分析地域对生物量的影响；符利勇等(2014)利用比例平差法和非线性联立方程组法建立不同起源地上生物量以及干材、干皮、树枝和树叶各分项生物量相容的通用性模型，得出非线性联立方程组法所构建的各分项生物量模型精度较高；Fu等(2017a)基于1950—2000气候数据，采用非线性混合效应模型构建气候敏感的非线性混合效应地上总生物量模型，结果表明生长季总降雨量、生长季平均温度、最湿季平均温度和最湿季总降雨量对马尾松地上总生物量具有显著影响，与未考虑气候因子的传统生物量模型相比，气候敏感的地上生物量模型可明显改进模型预测精度。

IPCC(2007)研究表明，近100年来，全球地表平均气温上升了0.74 ℃，近50年(1956—2005)的增温趋势尤为明显。国家气候变化评估报告显示，近100年来我国平均温度上升了0.5~0.8 ℃，预计到2100年地表平均温度将升高3.9~6.0 ℃(国家气候变化评估报告编辑委员会，2007)。气候变化除了温度升高外，降水的时空格局也会发生变化(IPCC，2007)，因此，气候变化已成为影响森林生物量的主要因子(Zeng et al., 2017；Fu et al., 2017a)。Fu等(2017a)分析了气候变化对马尾松生物量的影响，但是其研究对象为地上总生物量，对于气候变化与各分项生物量，如干材、干皮、树枝和树叶生物量之间的关系，当前尚未开展研究。鉴于此，本研究在Fu等(2016；2017a)所构建生物量模型的基础上，基于覆盖我国南方整个马尾松分布区域的生物量单木数据，采用非线性联立方程组法构建气候敏感的马尾松生物量相容性方程系统，分析气候因子对马尾松各分项生物量的影响，以期为森林碳汇监测和森林可持续经营提供技术支撑。

1 数据 1.1 生物量数据

生物量数据为我国南方马尾松立木地上生物量实测数据，共150株样木，采集时间为2009年6—9月，采集地点为马尾松主要分布区的江苏、浙江、安徽、福建、江西、湖南、广东、广西和贵州9省(区)，覆盖我国亚热带地区(包括北亚热带、中亚热带和南亚热带)。样木数按2、4、6、8、12、16、20、26、32和38 cm以上共10个径阶均匀分配(除26、32 cm径阶分别为14、16株外，其他径阶均为15株)，每个径阶的样木数按树高级从低到高也尽量均匀分配，在大尺度范围内具有广泛代表性。全部样木实测胸径、地径和冠幅，伐倒后，测量其树干长度(树高)和活树冠长度(冠长)，分干材、干皮、树枝和树叶称鲜质量，并分别抽取样品带回实验室，85 ℃恒温下烘干至恒重，根据样品鲜质量和干质量分别推算出样木各部分干质量并汇总得到地上部分干质量。150株样木生物量数据统计见表 1，具体数据调查、测定方法、各径阶和气候带的抽样百分比分布以及地上生物量(AGB)与直径(D)和树高(H)的散点分布见符利勇等(2011a)。

1.2 气候数据

气候数据从WorldClim数据库(www.worldclim.org)下载得到(Hijmans et al., 2005)。每块样地下载19个用于构建气候敏感的生物量相容性方程系统的候选因子，分别为平均生长季温度(mean growing season temperature，MGT)、年平均温度(annual mean temperature，AMT)、温度月平均变化幅度(平均每个月的最高气温-最低气温) [mean diurnal range(mean of monthly maximum temperature-minimum temperature)，MDR]、温度季节性(temperature seasonality，TS)、最热月的最高气温(maximum temperature of the warmest month，MTWM)、最冷月的最低气温(minimum temperature of the coldest month，MTCM)、温度年变化幅度(temperature annual range，TAR)、最湿季平均温度(mean temperature of the wettest quarter，TWQ)、最干季平均温度(mean temperature of the driest quarter，TDQ)、最热季平均温度(mean temperature of the warmest quarter，MTQ)、最冷季平均温度(mean temperature of the coldest quarter，MTCQ)、生长季总降雨量(total growing season precipitation，TGP)、年降雨量(Annual precipitation，AP)、最湿月降雨量(precipitation of the wettest month，PWM)、最干月降雨量(precipitation of the driest month，PDM)、最湿季降雨量(precipitation of the wettest quarter，PWQ)、最干季降雨量(precipitation of the driest quarter，PDQ)、最热季降雨量(precipitation of the hottest quarter，PHQ)、最冷季降雨量(precipitation of the coldest quarter，PCQ)，分辨率均为1 km×1 km。19个候选因子按照属性可分为与温度有关的变量(11个)和与降雨量有关的变量(8个)，每个气候因子为1950—2000年的平均值。马尾松生长季由气温确定，即每年从连续5日的日平均气温大于6 ℃开始，到当年8月第1个日气温低于10.5 ℃结束(Feng et al., 2011)。

2 方法

Fu等(2016)提出总量控制联立方程组法且各分项生物量基础模型为二元异速方程对应的马尾松生物量相容性方程系统预测效果最好，本研究在Fu等(2016)所构建生物量方程系统的基础上进一步考虑气候因子，重点分析气候变化对各分项生物量的影响。

2.1 总量控制联立方程组

总量控制联立方程组是指由总量直接平差分配给干材、干皮、树枝和树叶，从而保证各分量之和等于总量(唐守正等, 2000；Fu et al., 2016)。模型表达式如下：

$ \left\{ \begin{align} & {{w}_{1i}}=\frac{{{f}_{1}}\left( {{x}_{1i}},{{\phi }_{1}} \right)}{{{f}_{1}}\left( {{x}_{1i}},{{\phi }_{1}} \right)+{{f}_{2}}\left( {{x}_{2i}},{{\phi }_{2}} \right)+{{f}_{3}}\left( {{x}_{3i}},{{\phi }_{3}} \right)+{{f}_{4}}\left( {{x}_{4i}},{{\phi }_{4}} \right)}{{\dot{W}}_{0i}}; \\ \\ & {{w}_{2i}}=\frac{{{f}_{2}}\left( {{x}_{2i}},{{\phi }_{2}} \right)}{{{f}_{1}}\left( {{x}_{1i}},{{\phi }_{1}} \right)+{{f}_{2}}\left( {{x}_{2i}},{{\phi }_{2}} \right)+{{f}_{3}}\left( {{x}_{3i}},{{\phi }_{3}} \right)+{{f}_{4}}\left( {{x}_{4i}},{{\phi }_{4}} \right)}{{\dot{W}}_{0i}}; \\ & {{w}_{3i}}=\frac{{{f}_{3}}\left( {{x}_{3i}},{{\phi }_{3}} \right)}{{{f}_{1}}\left( {{x}_{1i}},{{\phi }_{1}} \right)|+{{f}_{2}}\left( {{x}_{2i}},{{\phi }_{2}} \right)+{{f}_{3}}\left( {{x}_{3i}},{{\phi }_{3}} \right)+{{f}_{4}}\left( {{x}_{4i}},{{\phi }_{4}} \right)}{{\dot{W}}_{0i}}; \\ & {{w}_{4i}}=\frac{{{f}_{4}}\left( {{x}_{4i}},{{\phi }_{4}} \right)}{{{f}_{1}}\left( {{x}_{1i}},{{\phi }_{1}} \right)+{{f}_{2}}\left( {{x}_{2i}},{{\phi }_{2}} \right)+{{f}_{3}}\left( {{x}_{3i}},{{\phi }_{3}} \right)+{{f}_{4}}\left( {{x}_{4i}},{{\phi }_{4}} \right)}{{\dot{W}}_{0i}}; \\ &{{Y}_{i}}={{y}_{i}}+{{e}_{i}},{{Y}_{i}}={{\left( {{W}_{1i}},{{W}_{2i}},{{W}_{3i}},{{W}_{4i}} \right)}^{\text{T}}},{{y}_{i}}={{\left( {{w}_{1i}},{{w}_{2i}},{{w}_{3i}},{{w}_{4i}} \right)}^{\text{T}}}; \\ & E\left( {{e}_{i}} \right)=0,\text{var}\left( {{e}_{i}} \right)={{\sigma }^{2}}\psi ,i=1,\cdots ,\text{ }N。\end{align} \right. $

(1)

式中：$\hat{W}_{0 i}$为第i(i=1, …, N)株样木地上总生物量估计值，由独立的地上总生物量方程计算得到；W_1i、W_2i、W_3i和W_4i分别为第i株样木对应的干材、干皮、树枝和树叶生物量观测值，假定各观测值含有度量误差(唐守正等，2015)；w_0i，w_1i、w_2i、w_3i和w_4i分别为W_0i、W_1i、W_2i、W_3i和W_4i对应的真实值；f₁(x_1i, ϕ₁)、f₂(x_2i, ϕ₂)、f₃(x_3i, ϕ₃)和f₄(x_4i, ϕ₄)分别为干材、干皮、树枝和树叶生物量的基础模型，x_1i、x_2i、x_3i和x_4i分别为各基础模型对应的自变量，ϕ₁、ϕ₂、ϕ₃和ϕ₄分别为各基础模型对应的待估参数；e_i为Y_i=(W_1i, W_2i, W_3i, W_4i)^T对应的残差向量，假定服从期望为0、方差为σ²ψ的正态分布，本研究中ψ假定为4×4维的无结构类型误差结构矩阵。

Fu等(2016)从D、H、D₀、Age、CH和CW 6个变量中选取不同变量构建一元、二元和三元生物量模型，并采用加权最小二乘回归法消除模型中存在的异方差。基于独立形式的一元、二元和三元生物量模型，利用非线性似然无关回归法、比例平差法和非线性联立方程组法构建相应的相容性或可加性生物量方程系统。通过比较，由模型(1)所构建的二元总量控制生物量联立方程组预测精度较高，即f_j(x_ji, ϕ_j)=a_jD_i^b_jH_i^c_j(j=1, 2, 3, 4)，式中：a_j、b_j和c_j(j=1, 2, 3, 4)分别为各分项生物量方程对应的参数；D_i和H_i分别为第i(i=1, …, N)株样木对应的胸高直径(cm)和树高(m)。方程系统表达式如下：

$ \left\{ \begin{align} & {w_{1i}} = \frac{{{a_1}D_i^{{b_1}}H_i^{{c_1}}}}{{{a_1}D_i^{{b_1}}H_i^{{e_1}} + {a_2}D_i^{{b_2}}H_i^{{c_2}} + {a_3}D_i^{{b_3}}H_i^{{e_3}} + {a_4}D_i^{{b_4}}H_i^{{c_4}}}}{{\dot W}_{0i}};\\ & {w_{2i}} = \frac{{{a_2}D_i^{{b_2}}H_i^{{e_2}}}}{{{a_1}D_i^{{b_1}}H_i^{{e_1}} + {a_2}D_i^{{b_2}}H_i^{{c_2}} + {a_3}D_i^{{b_3}}H_i^{{c_3}} + {a_4}D_i^{{b_4}}H_i^{{c_4}}}}{{\dot W}_{0i}};\\ & {w_{3i}} = \frac{{{a_3}D_i^{{b_3}}H_i^{{e_3}}}}{{{a_1}D_i^{{b_1}}H_i^{{c_1}} + {a_2}D_i^{{b_2}}H_i^{{c_2}} + {a_3}D_i^{{b_3}}{H^3} + {a_4}D_i^{{b_4}}H_i^{{c_4}}}}{{\dot W}_{0i}};\\ & {w_{4i}} = \frac{{{a_4}D_i^{{b_4}}H_i^{{e_4}}}}{{{a_1}D_i^{{b_1}}H_i^{{c_1}} + {a_2}D_i^{{b_2}}H_i^{{c_2}} + {a_3}D_i^{{b_3}}H_i^{{e_3}} + {a_4}D_i^{{b_4}}H_i^{{c_4}}}}{{\dot W}_{0i}};\\ & {{Y_i} = {y_i} + {e_i}, {Y_i} = {{\left({{W_{1i}}, {W_{2i}}, {W_{3i}}, {W_{4i}}} \right)}^{\rm{T}}};}\\ & {{y_i} = {{\left({{w_{1i}}, {w_{2i}}, {w_{3i}}, {w_{4i}}} \right)}^{\rm{T}}}。} \end{align} \right. $

(2)

2.2 气候敏感的马尾松生物量相容性方程系统

模型(2)中各分项(干材、干皮、树枝和树叶)生物量方程均为异速方程：

$ W_{j i}=f_{j}\left(D_{i}, H_{i}, \phi_{j}\right)=a_{j} D_{i}^{b_{j}} H_{i}^{\epsilon_{j}}+\varepsilon_{j i}, j=1, 2, 3, 4。$

(3)

式中：W_ji为第i(i=1, …, N)株样木对应的第j分项生物量；ε_ji为误差项。

为了构建不同起源(人工林和天然林)马尾松通用的生物量方程系统，创建一个哑变量S用于区别不同起源之间的差异效应，即：S=1表示人工林；S=0表示天然林。与Fu等(2017a)相似，为区别亚热带各气候亚区(包括北亚热带、中亚热带和南亚热带)对生物量影响的差异程度，构建以下哑变量：T₁=1、T₂=0表示南亚热带区；T₁=0、T₂=1表示中亚热带区；T₁=0、T₂=0表示北亚热带区。假定所有哑变量作用在参数a_j上，各分项生物量模型变为：

$ \begin{array}{l} {W_{ji}} = {f_j}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\phi _j}} \right) = \\ \begin{array}{*{20}{c}} {\left({{a_j} + {L_j}{S_i} + {k_{1j}}{T_{1i}} + {k_{2j}}{T_{2i}}} \right)D_i^{{b_j}}H_i^{{\varepsilon _j}} + {\varepsilon _{ji}}}\\ {j = 1, 2, 3, 4。} \end{array} \end{array} $

(4)

式中：L_j、k_1j和k_2j为模型参数。

为了定量描述气候变化对马尾松各分项生物量的影响，本研究在模型(4)的基础上进一步考虑温度和湿度有关的气候因子，利用逐步回归法分析各气候因子的不同组合对各分项生物量的影响以及对模型精度的贡献程度(唐守正等，2008)，在显著水平α=0.05情形下，如果变量显著影响模型精度将被用于各分项气候敏感的生物量模型构建。同时，为了避免变量冗余，本研究利用方差膨胀因子剔除共线性变量。假定F_j为第j分项生物量所选的气候因子向量，各分项气候敏感的生物量模型为W_j=f_j(D_i, H_i, S_i, T_i, F_ji, ϕ_j)，用f_j(D_i, H_i, S_i, T_i, F_ji, ϕ_j)替换模型(1)中的f_j(x_ji, ϕ_j)，即为最终构建的气候敏感的马尾松生物量相容性方程系统。非线性联立方程组及其参数估计方法详细介绍见唐守正等(2008)。需要指出的是，地上马尾松$\dot{W}_{0 i}$的估计本研究直接选用Fu等(2017a)所构建的马尾松气候敏感的地上生物量模型，表达式为：

$ \begin{array}{l} {W_{0i}} = \exp \left({{\beta _0} + L \times {S_i} + {k_1}{T_{1i}} + {k_2}{T_{2i}}} \right)\\ D_i^{{\beta _1}}{\rm{TGP}}_i^{{\beta _2}}{\rm{MGT}}_i^{{\beta _3}}{\rm{TWQ}}_i^{{\beta _4}}{\rm{PWQ}}_i^{{\beta _{\rm{s}}}} + \varepsilon 。\end{array} $

(5)

式中：W_0i为第i株样木对应的地上生物量(i=1, …, N)；β₀~β₅、L、k₁和k₂为模型参数；TGP_i、MGT_i、TWQ_i和PWQ_i分别为第i株样木对应的生长季总降雨量、生长季平均温度、最湿季平均温度和最湿季降雨量。

2.3 模型评价

利用平均残差(e)、残差方差(δ)、均方根误差(RMSE)和总相对误差(TRE)对模型进行评价(Fu et al., 2017a；2017b)：

$ \bar e = \sum {{e_i}} /N = \sum {\left({{y_i} - {{\hat{y}}_i}} \right)} /N; $

(6)

$ \delta = \sum {{{\left({{e_i} - \bar e} \right)}^2}} /(N - 1); $

(7)

$ {\rm{RMSE}} = \sqrt {{{\bar e}^2} + \delta } ; $

(8)

$ {\rm{TRE}} = \sum_{i = 1}^N {\left({{y_i} - \mathop {{{\hat{y}}_i}}\limits^. } \right)} /\sum_{i = 1}^N {{{\hat{y}}_i}} \times 100\% 。$

(9)

式中：y_i和$\hat{y}_{i}$分别为第i(i=1, …, N)株单木生物量各分项实测值(kg)和估计值(kg)；N为总观测数。

鉴于试验数据较少，为了更准确评价所构建的模型，本研究选用一阶交叉验证法，即在每次计算时删除1株马尾松，剩余数据用于估计模型参数，得到模型参数估计量之后预测所删除的马尾松各分项生物量，共计算150次。逐步回归计算在ForStat2.2(唐守正等，2008)上实现，其他所有计算在R3.4.0(R Core Team, 2012)上实现。

3 结果与分析 3.1 基础模型

基于逐步回归方法和方差膨胀因子检验，最终确定与各分项生物量相关性最强的前4个气候因子。对于干材，气候因子分别为温度月平均变化幅度(MDR)、温度季节性(TS)、最冷季平均温度(MTCQ)和生长季总降雨量(TGP)；对于干皮，气候因子分别为年平均温度(AMT)、最干季平均温度(TDQ)、最热季降雨量(PHQ)和平均生长季温度(MGT)；对于树枝，气候因子分别为年降雨量(AP)、最冷季降雨量(PCQ)、平均生长季温度(MGT)和生长季总降雨量(TGP)；对于树叶，气候因子分别为最湿季平均温度(TWQ)、最热月的最高气温(MTWM)、年降雨量(AP)和平均生长季温度(MGT)。将所选中的气候因子加入到各分项生物量模型(4)中，得到最终用于构建气候敏感的马尾松生物量相容性方程系统的各分项生物量基础模型：

$ \begin{array}{l} {W_{1i}} = {f_1}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{MD}}{{\rm{R}}_i}, {\rm{T}}{{\rm{S}}_i}, {\rm{MTC}}{{\rm{Q}}_i}, {\rm{TG}}{{\rm{P}}_i}, {\phi _1}} \right) + {\varepsilon _{1i}} = \\ \left({{a_1} + {L_1}{S_i} + {k_{11}}{T_{1i}} + {k_{21}}{T_{2i}} + {\beta _{11}}\ln MD{R_i} + } \right.\\ {\beta _{21}}\ln {\rm{T}}{{\rm{S}}_i} + {\beta _{31}}\ln {\rm{MTC}}{{\rm{Q}}_i} + \\ {\beta _{41}}\ln {\rm{TG}}{{\rm{P}}_i})D_i^{{b_1}}H_i^{{c_1}} + {\varepsilon _{1i}}; \end{array} $

(10)

$ \begin{array}{l} {W_{2i}} = {f_2}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{AM}}{{\rm{T}}_i}, {\rm{TD}}{{\rm{Q}}_i}, {\rm{PW}}{{\rm{Q}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\phi _2}} \right) + {\varepsilon _{2i}} = \\ \left({{a_2} + {L_2}{S_i} + {k_{12}}{T_{1i}} + {k_{22}}{T_{2i}} + {\beta _{12}}\ln {\rm{AM}}{{\rm{T}}_i} + } \right.\\ {\beta _{22}}\ln {\rm{TD}}{{\rm{Q}}_i} + {\beta _{32}}\ln {\rm{PW}}{{\rm{Q}}_i} + \\ {\beta _{42}}\ln {\rm{MG}}{{\rm{T}}_i}|)D_i^{{b_2}}H_i^{{c_2}} + {\varepsilon _{2i}}; \end{array} $

(11)

$ \begin{array}{l} {W_{3i}} = {f_3}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{A}}{{\rm{P}}_i}, {\rm{PC}}{{\rm{Q}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\rm{TG}}{{\rm{P}}_i}, {\phi _3}} \right) + {\varepsilon _{3i}} = \\ \left({{a_3} + {L_3}{S_i} + {k_{13}}{T_{1i}} + {k_{23}}{T_{2i}} + {\beta _{13}}\ln A{P_i} + } \right.\\ {\beta _{23}}\ln {\rm{PC}}{{\rm{Q}}_i} + {\beta _{33}}\ln {\rm{MG}}{{\rm{T}}_i} + \\ {\beta _{43}}\ln {\rm{TG}}{{\rm{P}}_i})D_i^{{b_3}}H_i^{{c_3}} + {\varepsilon _{3i}}; \end{array} $

(12)

$ \begin{array}{l} {W_{4i}} = {f_4}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{TW}}{{\rm{Q}}_i}, {\rm{MTW}}{{\rm{M}}_i}, {\rm{A}}{{\rm{P}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\phi _4}} \right) + {\varepsilon _{4i}} = \\ \left({{a_4} + {L_4}{S_i} + {k_{14}}{T_{1i}} + {k_{24}}{T_{2i}} + {\beta _{14}}\ln {\rm{TW}}{Q_i} + } \right.\\ {\beta _{24}}\ln {\rm{MT}}{{\rm{W}}_i} + {\beta _{34}}\ln {\rm{A}}{{\rm{P}}_i}| + \\ {\beta _{44}}\ln {\rm{MG}}{{\rm{T}}_i})D_i^{{b_4}}H_i^{{c_4}} + {\varepsilon _{4i}}。\end{array} $

(13)

式中：β₁₁~β₄₁、β₁₂~β₄₂和β₁₃~β₄₃、β₁₄~β₄₄分别为各分项生物量模型参数。

利用全部试验数据对各分项对应的生物量基础模型(3)、考虑起源和气候亚区但不考虑气候因子的基础模型(4)以及同时考虑起源、气候亚区和气候因子的基础模型(10)~(13)进行拟合，相应的参数估计值见表 2，通过T检验可知，各模型参数显著不等于0(α=0.05)。干材、干皮、树枝和树叶对应的拟合指标决定系数(R²)，对于模型(3)分别为0.95、0.92、0.87和0.61，对于模型(4)分别为0.95、0.93、0.89和0.64，对于模型(10)~(13)分别为0.96、0.95、0.90和0.60。这表明模型(3)、(4)和(10)~(13)都有较好的拟合精度，当模型进一步考虑起源和气候亚区以及同时考虑起源、气候亚区和气候因子时，拟合精度进一步提高。

各分项所对应基础模型的一阶交叉检验评价指标见表 3。与传统基础模型(3)相比，当考虑林分起源和气候亚区(南亚热带、中亚热带和北亚热带)对各分项生物量的影响时，能提高模型预测精度，如与模型(3)相比，模型(4)对应干材、干皮、树枝和树叶的TRE分别下降7.04%、18.90%、19.68%和6.47%。在模型(4)的基础上，当进一步考虑气候因子时，各分项生物量模型预测精度进一步提高，如与模型(4)相比，模型(10)~(13)对应干材、干皮、树枝和树叶的RMSE分别下降23.51%、11.85%、4.72%和7.78%，TRE分别下降42.51%、22.85%、9.60%和19.18%。这进一步表明气候因子对干材、干皮、树枝和树叶生物量的影响较大。

3.2 气候敏感的马尾松生物量相容性方程系统

将模型(10)~(13)代入模型(1)中，得到气候敏感的马尾松生物量相容性方程系统，表达式如下：

$ \left\{ \begin{align} & {w_{1i}} = \frac{{{f_1}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{MD}}{{\rm{R}}_i}, {\rm{T}}{{\rm{S}}_i}, {\rm{MTC}}{{\rm{Q}}_i}, {\rm{TG}}{{\rm{P}}_i}, {\phi _1}} \right)}}{T}{{\hat W}_{0i}};\\ & {w_{2i}} = \frac{{{f_2}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{AM}}{{\rm{T}}_i}, {\rm{TD}}{{\rm{Q}}_i}, {\rm{PW}}{{\rm{Q}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\phi _2}} \right)}}{T}{{\hat W}_{0i}};\\ & {w_{3i}} = \frac{{{f_3}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{A}}{{\rm{P}}_i}, {\rm{PC}}{{\rm{Q}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\rm{TG}}{{\rm{P}}_i}, {\phi _3}} \right)}}{T}{{\hat W}_{0i}};\\ & {w_{4i}} = \frac{{{f_4}\left({{D_i}, {H_i}, {S_i}, {T_i}, {\rm{TW}}{{\rm{Q}}_i}, {\rm{MTW}}{{\rm{M}}_i}, {\rm{A}}{{\rm{P}}_i}, {\rm{MG}}{{\rm{T}}_i}, {\phi _4}} \right)}}{T}{{\hat W}_{0i}};\\ & {Y_i} = {y_i} + {e_i}, {Y_i} = {\left({{W_{1i}}, {W_{2i}}, {W_{3i}}, {W_{4i}}} \right)^{\rm{T}}};\\ & {y_i} = \left({{w_{1i}}, {w_{2i}}, {w_{3i}}, {w_{4i}}} \right)^{\rm{T}}。\end{align} \right. $

(14)

式中：T=f₁(D_i, H_i, S_i, T_i, MDR_i, TS_i, MTCQ_i, TGP_i, ϕ₁)+f₂(D_i, H_i, S_i, T_i, AMT_i, TDQ_i, PWQ_i, MGT_i, ϕ₂)+f₃(D_i, H_i, S_i, T_i, AP_i, PCQ_i, MGT_i, TGP_i, ϕ₃)+f₄(D_i, H_i, S_i, T_i, TWQ_i, MTWM_i, AP_i, MGT_i, ϕ₄)。

模型(14)不仅考虑了干材、干皮、树枝和树叶各分项生物量之间的相关性，同时也可保证各分项生物量之和等于地上总生物量。利用全部试验数据对模型(14)进行拟合，得到模型参数估计值见表 4，通过T检验可知，各模型参数显著不等于0(α=0.05)。干材、干皮、树枝和树叶对应的R²分别为0.96、0.95、0.90和0.59，表明对于干材、干皮、树枝，模型(14)具有较高的拟合精度，而对于树叶，模型(14)的拟合精度相对较低。

3.3 模型评价

表 5列出了一阶交叉验证法对应模型(14)的评价指标。模型(14)对应各分项生物量的平均残差非常小，尤其是干皮和树叶接近0。与传统基础模型(3)相比，模型(14)对应干材、干皮、树枝和树叶的RMSE分别下降26.09%、20.54%、13.77%和10.01%，TRE分别下降46.52%、37.42%、27.20%和24.40%，这表明当考虑林分起源、气候亚区(南亚热带、中亚热带和北亚热带)和气候因子对各分项生物量的影响时，能提高模型预测精度。与模型(4)相比，模型(14)对应干材、干皮、树枝和树叶的RMSE分别下降23.37%、12.37%、4.60%和7.85%，TRE分别下降42.47%、22.83%、9.37%和19.18%，这进一步表明气候因子对干材、干皮、树枝和树叶生物量的影响较大。与模型(10)~(13)相比，模型(14)对应各分项生物量的预测精度略有下降，但模型(14)能保证所预测的干材、干皮、树枝和树叶生物量之和等于总生物量。图 1为模型(14)对应各分项生物量残差分布，对于干材、树皮、树枝和树叶尽管存在异方差，但是不明显，进一步说明模型(14)预测各分项生物量效果较好。

4 讨论

随着全球气候变暖，树木各分项(如干材、干皮、树枝和树叶)的生长都随气候变化而变化，因此，在现有生物量建模方法(Parresol，1999；2001；Bi et al., 2004；Zeng et al., 2011；2012；2014；Zeng, 2015；Fu et al., 2012；2014；Dong et al., 2015；Zou et al., 2015；曾伟生等，2010a；2010b；符利勇等；2011a；2011b；2014)基础上研究气候变化下的各分项生物量可加性模型显得更为重要(Fu et al., 2017a；Zeng et al., 2017)。本研究假定气候因子作用在异速方程常数项上来分析气候变量对各分项生物量的影响，实际应用中还可利用气候因子对D和H所对应的参数进行参数化。与传统未考虑气候因子的基础模型相比，本研究构建的气候敏感的马尾松生物量相容性方程系统能明显提高模型预测精度，表明考虑气候因子对各分项生物量的影响对减少或避免生物量模型的不确定性有着重要意义。

本研究分析了温度和降雨量2类因子对马尾松各分项生物量的影响，从选中的气候因子可知，平均生长季温度对各分项(尤其是干皮、树枝和树叶)生物量影响最大，这是因为在亚热带地区，马尾松所处生长季正好属于雨季，雨水充沛，因此限制马尾松生长的最主要因子是温度。关于各分项生物量，干材与最冷季平均温度呈反比，是因为每年最冷季处于马尾松休眠季，平均温度越高，马尾松自身消耗的碳越多，而自身却不生长，因此使得干材与其呈负相关。干皮与最干季平均温度呈负相关，是因为最干季本身缺少水分，如果温度升高将促进水分蒸发，因此阻碍了树皮的生长。树叶与最热月的最高气温呈反比，与其他温度呈正比，这主要是因为最热月中阻碍树叶生长的因子是水分，温度越高，水分蒸发越快，从而抑制了树叶的生长。总体而言，与降雨量相比，马尾松各分项生物量受温度的影响较大。

本研究重点是如何构建气候敏感的马尾松生物量相容性方程系统，因此气候数据只应用历史数据(1950—2000年)而没有对不同未来气候模式下马尾松各分项生物量进行预测。需注意的是，本研究构建的模型没有考虑气候因子对马尾松直径和树高的影响，因此当利用模型(14)分析未来气候变化对马尾松各分项生物量的影响时，需给定马尾松直径和树高，即分析未来气候变化对相同直径和树高的马尾松的影响。如要分析同一样木未来不同时间点气候变化对其生物量的影响，还需在模型(14)基础上进一步考虑直径和树高受气候因子的影响。此外，WorldClim数据库中只能下载1950—2000年的气候数据，2001—2009年的气候数据不能下载，客观上讲，缺失2001—2009年的气候数据可能会对模型结果有一定影响，但由于本研究是应用近50年气候的平均值来构建生物量模型，并且短期几年气候变化不会很明显，因此该影响不会很大。除此之外，与Fu等(2017a)相类似，本研究侧重建模方法，对于气候变化对马尾松生物量的影响还需进一步详细研究。

本研究构建的气候敏感的马尾松生物量相容性方程系统只适用于模型参数不受随机因素影响的试验条件(唐守正等，2015)，但现实中，由于每株树木自身生长状况、遗传因素、所处立地和生境条件不同，各分项生物量与林木因子之间的关系可能存在较大的随机差异。对于既能反映总体平均变化情形又能反映个体间随机差异的各分项生物量相容性方程系统，即考虑随机干扰的生物量相容性方程系统迄今尚未报道，因此如何应用非线性混合效应模型联立方程组方法构建各分项气候敏感的生物量相容性方程系统是后续需解决的问题。尽管一些学者已利用度量误差联立方程组构建了生物量相容性系统(曾伟生等，2010b；Zeng et al., 2012)，但是现有方法大多假定方程系统中各分项生物量含有度量误差，而其他变量都未含有度量误差(实质上也是传统的非线性联立方程组)。现实中，尤其在林业领域，许多林分或树木因子实际测量时都存在一定程度的度量误差。度量误差模型理论指出，当模型中因变量和自变量都存在度量误差时，采用传统方法估计的参数为有偏估计，即使样本数量非常大也无法减少估计偏差(统计学中称为没有“相合性”)(Fu et al., 2016)。前期分析中，本研究假定树高和直径存在度量误差，即利用含度量误差的非线性联立方程组构建气候敏感的生物量相容性方程系统，但是计算不收敛。Fu等(2016)利用非线性联立方程组构建马尾松各分项生物量相容性方程时发现，当直径和树高都考虑度量误差时能提高模型预测精度，但是提高的幅度非常小。尽管本研究尚未考虑生长状况、遗传因子等对各分项生物量的随机干扰以及直径和树高中的度量误差对各分项生物量的影响，但总体而言所提出的气候敏感的马尾松生物量相容性方程系统能有效分析气候因子对马尾松生物量的影响，可为我国生物量以及碳汇监测和预估提供理论支撑。

5 结论

针对我国南方亚热带马尾松树种，本研究构建了一个不同林分起源(人工林和天然林)、不同亚热带亚区气候带(南亚热带、中亚热带和北亚热带)通用的气候敏感的马尾松生物量相容性方程系统。为了揭示气候变化对马尾松各分项(干材、干皮、树枝和树叶)生物量的影响，模型考虑了12个气候因子，其中，温度月平均变化幅度(平均每个月的最高气温-最低气温)、温度季节性、最冷季平均温度和生长季总降雨量对干材生物量影响显著；年平均温度、最干季平均温度、最热季降雨量和平均生长季温度对干皮生物量影响显著；年降雨量、最冷季降雨量、平均生长季温度和生长季总降雨量对树枝生物量影响显著；最湿季平均温度、最热月的最高气温、年降雨量和平均生长季温度对树叶生物量影响显著。而且，模型还考虑了3个气候带亚区对马尾松各分项生物量的影响差异。结果表明，与传统未考虑气候因子的各分项生物量模型相比，本研究构建的气候敏感的马尾松生物量相容性方程系统预测精度明显提高。相同直径条件下，在最湿季，平均气温过高或最湿季降雨量过大将减少生物量，但在生长季，大量降雨或温度较高会增加单木生物量，且模型能定量描述不同气候带亚区生物量差异。总之，本研究构建的气候敏感的马尾松生物量相容性方程系统能有效分析气候因子对生物量的影响，可应用于其他树种的生物量预估。