在全球气候变化背景下，陆地生态系统碳循环对环境因子变化的响应是生态学研究的热点问题。作为陆地生态系统的主体，森林是一个巨大的碳库，其储存的碳量是大气的3倍，并且对全球气候变化极其敏感(王春林等，2006; Beer et al., 2010)。因此，准确测定和估算森林生态系统和大气间的碳交换量，对于评价区域碳源、汇功能以及深入研究森林生态系统对全球气候变化的反馈机制具有重要意义(王春林等，2006; Jha et al., 2013)。

目前，碳通量的观测主要有清单法、箱法和微气象法(王兴昌等，2015)，其中，微气象法中的涡度相关技术法可以在半小时尺度上直接、连续地观测植被群落和大气之间的气体交换通量及其环境控制因子，现已广泛用于各个生态系统的通量观测研究中，成为了一种通用的标准方法(Baldocchi et al., 2001; Law et al., 2002)。研究表明，森林生态系统CO₂通量易受温度、降水、光合有效辐射、林龄和树种组成等因子的影响, 且各因子之间交互作用对CO₂通量的影响也不可忽视(Saigusa et al., 2002; Barr et al., 2007; Allard et al., 2008; Grossiord et al., 2013; Fernández-Martínez et al., 2014; Xie et al., 2015)，因此，应用传统统计回归模型估测森林生态系统CO₂通量时，难以克服诸多因子之间复杂的复共线性问题，且对数据中存在的噪声和缺失数据敏感，极易造成模型精度降低(张雷等，2014; Lopatin et al., 2016)。此外，一些学者运用人工神经网络模型(ANN)拟合森林生态系统CO₂通量(He et al., 2006; Wen et al., 2014)，该方法虽然有强大的自身学习能力和较强的非线性映射能力，但也存在一些不可避免的问题，如对训练数据依赖性较高、容易出现过度拟合等(李旭青等，2014)。随机森林模型是由Breiman (2001)提出的数据挖掘算法，该算法属于“灰箱”模型，其结果具有可解释性，并且可以有效解决过度拟合和复共线性等问题，在生态学领域具有广阔应用前景。

东亚季风区主要分布在20°—40°N之间，该区域森林林型多样、植被类型丰富，具有很高的净CO₂吸收强度，每年的净生态系统生产力(NEP)达3.6 Mg C·hm^-2，并且森林生态系统NEP总量约0.72 Pg C·a^-1，占全球森林生态系统NEP的8%，在全球碳循环中扮演着重要角色(Yu et al., 2014;方精云等，2015)。毛竹(Phyllostachys edulis)作为该区域一种重要的林分类型，不仅具有可观的经济价值，而且具有巨大的固碳潜力，在减缓全球气候变化方面发挥了积极作用。在毛竹林生态系统固碳能力和碳分配方面，已有许多学者做了大量研究(Zhou et al., 2011; Tang et al., 2015)，但毛竹林碳循环形成机制、CO₂通量对环境因子的响应等方面还相对缺乏深入探讨。鉴于此，本研究基于浙江省安吉县山川乡2011—2014年毛竹林涡度相关数据，采用随机森林模型对毛竹林CO₂通量进行模拟，并分析光合有效辐射(PAR)、大气温度(T_A)、空气相对湿度(R_H)、空气CO₂浓度(C_CO2)、土壤热通量(Q_S)、饱和水汽压差(VPD)、生态系统有效能量(A_E)、风速(W_S)和土壤5 cm深处温度(T_S)9个环境因子对CO₂通量的影响，以期为亚热带典型森林碳循环模拟提供技术和理论支撑。

1 研究区概况

研究区位于浙江省安吉县山川乡毛竹林碳汇研究基地(119°40′25.7″ E，30°28′34.5″ N)，海拔380 m，主要土壤类型为黄壤土。气候属典型的海洋性季风气候，四季分明，雨热同期，年均降水1 400 mm，年均气温16.6 ℃。2010年，浙江农林大学在该基地内建立了40 m涡度相关通量观测塔，以塔为中心1 km×1 km范围内主要森林类型为毛竹林，平均胸径为12~18 cm，平均高为13~20 m。此外，还有少量针阔混交林、农田、道路和建筑。

2 研究方法 2.1 数据来源及处理

本研究采用2011—2014年毛竹林通量塔CO₂通量和气象原始数据，时间分辨率为半小时，其中通量数据由于降水及传感器自身原因会出现异常值和缺失值。因此，本研究首先利用EdiRe软件(http://www.geos.ed.ac.uk/abs/research/micromet/EdiRe)对通量数据进行处理，包括WPL空气密度校正、二次坐标旋转和虚温订正(于贵瑞等，2006)；然后进行冠层储存项校正和摩擦风速校正(于贵瑞等，2006)；最后，对通量数据进行质量控制和插补。数据质量控制过程主要有：1)剔除摩擦风速小于0.2 m·s^-1的数据；2)剔除超出仪器合理测量范围的观测数据，如水汽浓度超出0~40 g·m^-3、CO₂通量超出－2.0~2.0 mg·m^-3。得到的有效数据比例如表 1所示，可以看出，除了2013年有效数据比例偏低，外其余各年有效数据均在正常范围内(50%~83%)(Katul et al., 2001)。

对于缺失数据，白天缺失数据利用Michaelis-Menten方程插补(Falge et al., 2001)：

$ {\rm{NE}}{{\rm{E}}_{\rm{d}}} = {R_{\rm{e}}} - {P_{\rm{m}}} \times \alpha \times {\rm{PAR}}/(\alpha \times {\rm{PAR}} + {P_{\rm{m}}})。$

(1)

式中：NEE_d为白天缺失的CO₂通量数据；R_e为生态系统暗呼吸速率(mg·m^-2s^-1)；P_m为冠层最大净光合速率(mg·m^-2s^-1)；α为初始光能利用率。

夜间缺失数据利用基于温度响应的Lloyd & Taylor方程插补(Lloyd et al., 1994)：

$ {\rm{NE}}{{\rm{E}}_{\rm{n}}} = {\rm{R}}{{\rm{E}}_{{\rm{ref}}}}{\rm{exp}}[A(1/283.16 - {T_0} - 1/{T_t} - {T_0})]。$

(2)

式中：NEE_n为夜间缺失的CO₂通量数据；RE_ref为参考温度(10 ℃)下生态系统呼吸速率；A为温度敏感参数(本研究取309 K)；T₀为温度常数(本研究取-46 ℃)；T_t为土壤5 cm深处温度。

最后，对通量和气象数据求半小时观测值的月平均，得到本研究所用数据。

2.2 随机森林模型

随机森林模型是一种基于决策树的组合模型，其回归过程如图 1所示。首先，对训练数据集利用自助抽样(Bootstrap)形成N个bootstrap样本；其次，对每个bootstrap样本分别建立回归树模型，且该样本对应回归树模型的全部训练数据和测试数据；最后，以N个回归树结果的平均值作为因变量的预测值。此外，该模型迫使每个二叉分裂点仅考虑预测变量的一个子集，从而使得弱相关性变量获得更多参与回归树模型构建的机会，因此所得到的回归树模型具有更高的可信度(Breiman, 2001)。为获得优化的随机森林模型，本研究以2011—2012年通量观测数据为模型训练样本，对随机森林模型进行训练；采用2013年观测数据作为测试样本，对通过训练的模型进行测试；用通过测试后的模型对2014年毛竹林CO₂通量进行模拟，即2014年数据作为独立验证样本。随机森林模型构建通过R语言randomForest包实现。

2.2.1 随机森林模型输入环境变量

随机森林模型输入环境变量包括PAR、T_A、R_N、C_CO2、Q_S、VPD、A_E、W_S和T_S，其中A_E=R_N－Q_S，R_N为太阳净辐射(王兴昌等，2015)，其余8个环境变量均由通量塔观测获取。Spearman相关分析表明，这9个环境变量与毛竹林CO₂通量均极显著相关(P < 0.001)。

2.2.2 随机森林模型精度评价

采用均方根误差(RMSE)、决定系数(R²)和林氏调和系数(LCCC)3个指标对训练、测试和验证数据集的随机森林模型模拟结果进行精度检验：

$ {\rm{RMSE}} = \sqrt {\frac{1}{n}\sum\nolimits_{i = 1}^n {{{({p_i} - {o_i})}^2};} } $

(3)

$ {R^2} = 1 - \sum\nolimits_{i = 1}^n {{{({p_i} - {o_i})}^2}} /\sum\nolimits_{i = 1}^n {{{\left({{o_i} - {{\bar o}_i}} \right.}^2};} $

(4)

$ {\rm{LCCC}} = 2r{\sigma _{\rm{o}}}{\sigma _{\rm{p}}}/[{\sigma _{\rm{o}}}^2 + {\sigma _{\rm{p}}}^2 + \left({\bar o - \bar p} \right)]。$

(5)

式中：o_i和p_i分别为通量塔观测值(以下称观测值)和模型拟合值(以下称拟合值)；σ_o和σ_p分别为观测值和拟合值的标准差；r为观测值和拟合值之间的皮尔逊相关系数；${\bar o} $和$ {\bar p}$分别为观测值和拟合值的均值。LCCC表征观测值和拟合值回归拟合线靠近1:1回归线的程度，取值为0~1，值越大表明模型拟合效果越好(Lin, 1989)。

为进一步检验模型泛化能力，基于2014年通量塔观测数据，利用十折交叉验证法(10-fold cross-validation)评价已构建好的随机森林模型的预测误差。十折交叉验证法是将所有数据随机分成大小相近的10组，先以第1组数据作为验证数据，剩下9组数据用作拟合模型；再以第2组数据作为验证数据，剩下9组数据用作拟合模型。依此类推，重复10次，以10次测试误差的平均值作为最后精度评价标准。

2.3 环境因子对模型的重要性

随机森林模型中，在拟合或分类时会用到每个通过自助抽样形成的大约2/3 bootstrap样本数据，剩余1/3样本数据称为袋外样本数据(OOB)，用于评估模型整体误差以及变量的重要性。若X_j为输入变量之一，则计算变量X_j在第n株树的重要性VI_n时，首先利用自助抽样的样本数据建立回归树模型T_n，计算OOB的预测误差，然后随机置换变量X_j的观测值，再次建立回归树模型T_n′，计算OOB′ (随机置换变量观测后袋外数据)的预测误差，将2次袋外数据预测误差的差值进行标准化处理后，在所有回归树中取平均值即为变量X_j在第n株树的重要性VI_n(X_j) (Delgado-Baquerizo et al., 2016)，其计算公式如下：

$ \begin{array}{l} {\rm{V}}{{\rm{I}}_n}({X_j}) = \{ \sum\nolimits_{i = 1}^{{{\rm{N}}_{{\rm{OOB}}}}} {I[f({X_i}) = {\rm{ }}{f_n}({X_i})]} -\\ \sum\nolimits_{i = 1}^{{{\rm{N}}_{{\rm{OOB}}}}} {I[f({X_i}) = {f_n}(X{\prime _i})]} \} /{{\rm{N}}_{{\rm{OOB}}}}。\end{array} $

(6)

变量X_j在整个随机森林模型中的重要性得分VI(X_j)为：

$ {\rm{VI}}({X_j}) = \sum\nolimits_{n = 1}^{{N_{\rm{t}}}} {{\rm{V}}{{\rm{I}}_n}({X_j})/{N_{\rm{t}}}。} $

(7)

式中：N_OOB为袋外样本数；f(X_i)为袋外数据中第i个观测值；f_n(X_i)为在随机置换变量X_j的观测值前第n株树上袋外数据第i个观测值所对应的预测值；$ {{f_n}(X{\prime _i})}$为在随机置换变量X_j的观测值后第n株树上袋外数据第i个观测值所对应的预测值；$ I\left[ {f({X_i}) = {f_n}({X_i})} \right]$和$I\left[ {f({X_i}) = {f_n}(X{\prime _i})} \right] $为判别函数，当$f({X_i}) = {f_n}({X_i}) $或$f({X_i}) = {f_n}(X{\prime _i}) $时，取值为1，否则为0；N_t代表随机森林模型中树的数量。

本研究以输入参数对模型的重要性得分及其影响显著性为依据，评价环境因子对毛竹林CO₂通量的影响。环境因子重要性得分及显著性检验通过R语言rfPermute包实现。

3 结果与分析 3.1 模型参数优化结果

随机森林模型中有3个重要参数：一是随机森林模型中树的数量(N_t)，在计算过程中，N_t取值足够大时，模型误差将趋于平稳，默认取值为500；二是分裂的变量个数M，即每个树节点上通过随机特征确定用于分裂的变量数, 默认为自变量数的1/3；三是终端节点最小样本数(N_S)，本研究采用模型默认取值，5。目前，随机森林模型中M大多基于经验公式判定，本研究采用逐次筛选确定最佳M，即遍历所有M的取值(输入变量数1~9)，取模型误差最小时所对应的M。图 2显示，当M取值为4时，模型误差达到最小值，而N_t达到500以后模型误差趋于稳定，因此最终M、N_t和N_S的取值分别为4、500和5。

3.2 环境变量变化特征

图 3为2011—2014年模型输入变量变化特征。PAR、T_A、Q_S、A_E、W_S和T_S的日变化均为单峰型，在12:00—2:00达到极值。6月进入梅雨季，PAR开始下降，7月PAR又迅速增大，PAR极大值出现在2011年4月，为1 237.7 W·m^-2；Q_S夜间为负值，表明土壤为热源，白天为正值，表明为热汇；2013年夏季出现极端高温现象，Q_S、T_A和T_S的极大值均出现在7月，分别为37.5 W·m^-2、33.9 ℃和27.3 ℃。A_E在－59.84~482.9 W·m^-2之间；W_S在各年夏季(6—8月)与同年其他月份相比均显著下降(P < 0.05)；VPD主要受空气水热条件影响，在2013年波动最为剧烈，由于7、8月受极端高温胁迫影响，VPD显著高于同年其他月份(P < 0.05), 其最大值为23.4 hPa。R_H和C_CO2在日变化中均出现“谷峰交替”现象，由于夜间植被呼吸和土壤呼吸作用导致空气中CO₂和水汽积累，在清晨5:00左右达到最大值，而后进行光合作用，CO₂和水汽作为光合作用反应底物被消耗，其浓度逐渐下降并在14:30左右达到最小值，R_H变化范围为36.7%~98.6%。

3.3 CO₂通量模拟结果及评价

图 4为模型在训练和测试阶段的拟合结果，图 5为模型在训练和测试阶段的精度评价。可以看出，在2011—2012年模型训练阶段，R²=0.961 5，RMSE= 0.005 4 mg·m^-2s^-1，LCCC=0.980 1，模型拟合效果较好(图 5A)。而在模型测试阶段，R²下降至0.845 5，RMSE升高到0.437 7 mg·m^-2s^-1，LCCC下降至0.914 1(图 5B)，并在2月12:00—14:00出现了低估情况，其极值低估了31.5%；在7月夜间，拟合值出现了较大波动，没有很好地体现该时段CO₂通量真实变化趋势；在8和12月，CO₂吸收峰值出现了明显高估，分别高估了57.3%和106%(图 4)。但总的来说，随机森林模型能较好地模拟CO₂通量的数值大小和变化趋势。

为了验证随机森林模型的稳定性和泛化能力，首先，利用2014年独立验证数据对随机森林模型进行十折交叉验证，结果(表 2)表明，R²和LCCC较高，最大值分别为0.976 4和0.977 1，均值为0.970 5和0.974 1，而RMSE较小，为0.007 5~0.017 8 mg·m^-2s^-1，均值为0.011 6 mg·m^-2s^-1，这说明随机森林模型模拟精度高，在毛竹林CO₂通量模拟中具有很好的稳定性。

其次，在日尺度上，选取2014年每月1日作为典型日，利用已训练的随机森林模型进行模拟，结果(表 3)表明，实测值与拟合值之间误差较小，R²为0.783 3，RMSE为0.010 6 mg·m^-2s^-1，线性关系也较为明显，LCCC为0.850 5，这说明随机森林模型在日尺度上也具有良好的拟合效果。

3.4 环境因子重要性得分

图 6为计算得到的环境因子对CO₂通量影响的重要性得分。可以看出，环境因子对CO₂通量影响的重要性得分表现为PAR(63.332) > T_S(29.932) > R_H(25.839) > T_A(25.581) > C_CO2 (25.095) > VPD(24.123) > W_S(23.504) > A_E(19.323) > Q_S(18.502)，PAR在影响毛竹林CO₂通量变化方面发挥了决定性作用。进一步分析表明，在0.05显著水平下，PAR、T_S和VPD对毛竹林CO₂通量的影响较为显著，说明这3个因子是影响CO₂通量月尺度变化的重要因素。

4 讨论 4.1 样本容量和模型参数对随机森林模拟的影响

虽然随机森林模型对输入参数间的共线性关系不敏感，但是该模型需要足够的样本容量来保证结果的合理性(张雷等，2014)，因此模型精度除了与自身内部参数有关外，还与训练、测试和独立验证阶段的样本容量有关，当样本容量较多时，模型精度较高，反之则较低。本研究随机森林模型在训练、测试和独立验证阶段样本容量分别为1 152、576和576，样本容量较大。

模型误差与其内部参数M、N_t和N_S有关，M过小，即在分裂节点上选择过少的变量特征，会使模型最终精度降低；过大，则减小模型的随机性，也会影响模型的最终精度(刘敏等，2015)。因此，M对模型精度的影响不能忽略。然而，对于M的选择，大多研究只是基于经验公式，如Yang等(2016)、Rahmati等(2015)的研究中M取值都为变量数的1/3，本研究选取的输入变量为9个，若取变量数的1/3，即3，根据图 4，显然不是最佳M取值，因此，本研究采用确定M的方法更为科学和严谨。N_t对随机森林模型的性能和可解释性具有重要意义(刘敏等，2015)，Breiman(2001)早已证明当N_t足够大时，随机森林模型误差会逐渐收敛。根据图 2，本研究N_t取值为500时，模型误差已经收敛于一个极小值，达到稳定。随机森林模型对N_S的取值并不敏感(刘敏等，2015)，因此N_S对模型误差的影响不予考虑，根据表 2和表 3结果，说明本研究随机森林模型参数设置合理，达到模型泛化要求。

4.2 环境变量对随机森林模型模拟的影响

模型误差还与输入环境变量性质和通量数据质量有关，通量数据质量除易受降水、大气湍流及传感器自身噪声影响外，还受数据质量控制方法(数据剔除方法、数据插补方法)的影响，当缺失数据较多时，数据质量会下降，导致模型模拟值与实际值偏差大。本研究中9个输入变量为通量塔观测的气象因子，均属于非生物因素，而森林生态系统CO₂通量还受其他生物因素的显著影响，如叶片气孔导度、水分利用效率、叶面积指数等(Wilson et al., 2001;张军辉等，2006; Keenan et al., 2013)，生物因素在本研究并未充分考虑。此外，全球气候变化增加了极端气候出现的频率，利用模型拟合森林生态系统CO₂通量时，往往难以刻画这种短时、强烈的干扰对CO₂通量的影响(Campbell et al., 2009)。如在2013年夏季发生的极端高温干旱，使得模型拟合值没有很好地反映该年7、8月CO₂通量真实值的大小和变化趋势，直接导致较大误差的产生，在模型测试阶段出现了精度下降现象。所以未来CO₂通量拟合研究除应综合考虑生理因素和环境因素对森林生态系统碳通量的影响(Froelich et al., 2015)外，还应充分考虑极端气候的影响，增加模型物候参数，降低极端干扰对模型模拟的影响，以提高模型精度(Campbell et al., 2009; Richardson et al., 2012)。虽然模型误差来源多样，但随机森林模型总体上能够较好地拟合2011—2014年毛竹林CO₂通量的数值大小和变化趋势。

4.3 CO₂通量对主要环境因子的响应

PAR是影响CO₂通量的主导因子，而T_S和VPD是重要的控制因子(图 6)。PAR是绿色植物生命活动、有机物质合成的重要能量来源，在一定范围内，PAR增大，植物光合作用增强，使植被固碳能力也相应增强；当伴随高温、高VPD出现时，为防止叶片失水过快，叶片气孔收缩甚至关闭，导致气孔导度下降；进而使CO₂通量显著下降。尤其在夏季，由于高温干旱，毛竹林很容易出现“光合午休”现象(孙成等，2013)，这导致2013年7、8月CO₂通量比2011年同期分别小20.2%和40.1%，比2012年同期分别小15.7%和36.0%(图 5)。在不同时间尺度上，PAR和其他气象因子往往占据不同主导地位。在月尺度上，PAR对CO₂通量的主导作用也出现在其他森林类型中，如Gonsamo等(2015)对北美洲东北部2种成熟林(温带落叶林和针阔混交林)碳通量的研究表明，PAR是影响2种成熟林碳通量的关键因子；Xie等(2015)对北京奥林匹克国家森林公园碳通量的研究表明，PAR是影响该森林在2012和2013年生长季碳通量变化的主导因子。而在年际间或更长时间尺度上，温度对CO₂通量的影响比PAR更为显著(Yang et al., 2014)。此外，PAR和CO₂通量之间的关系尚存在不确定性，在Xie等(2015)和Li等(2015)的研究中，PAR与CO₂通量满足Michaelis-Menten关系式，而彭镇华等(2009)对安庆杨(Populus)树林生态系统碳通量与PAR关系的研究表明，在较低温度下，PAR和碳通量为线性关系。因此，准确揭示PAR对森林生态系统碳循环的影响，有利于更好认识植被对于全球气候变化的响应(Gonsamo et al., 2015)。本研究中，T_S表示土壤5 cm深处温度，地表通常包含大量有机质(微生物、细根等)，是异氧呼吸的主要来源，T_S主要控制土壤呼吸进而影响森林生态系统CO₂通量(Desai et al., 2005)。环境因子对CO₂通量的影响具有时间尺度性(Yang et al., 2014)，对于不同时间尺度下毛竹林CO₂通量对环境因子变化的响应需要进一步研究。

5 结论

本研究基于随机森林模型对毛竹林CO₂通量进行模拟，结果表明，与BP人工神经网络模型相比，随机森林模型能够以致高精度模拟毛竹林CO₂通量，说明随机森林模型在拟合毛竹林碳通量方面具有较好的适用性。在模型测试阶段，由于高温干旱使得模拟精度下降。十折交叉验证表明，随机森林模型内部3个参数(N_t、M和N_S)设置合理，模型误差主要来源于输入数据。环境因子重要性得分表明，在0.05显著水平下，光合有效辐射、土壤5 cm深处温度和饱合水汽压差在月尺度上对毛竹林CO₂通量变化的贡献显著，说明这3个因子在月尺度上对毛竹林CO₂通量具有控制作用。为深入揭示毛竹林碳循环机制，今后应深入研究不同时间尺度下毛竹林碳通量对环境因子的响应。本研究结果可为深入了解亚热带典型森林生态系统碳通量控制机制及为区域生态系统过程模型参数调整提供科学依据。