基于过程模型的泰国热带森林碳通量动态模拟解析

廖立国，谭正洪，蒋 龙，符 妙，金 艳，刘应帅，章 杰

（海南大学 生态与环境学院，海口 570228）

森林生态系统作为陆地生态系统主体，在维持全球碳平衡、对抗气候突变等方面发挥着重要作用[1]。森林的生物量占陆地生态系统的80%，其碳储量约占全球陆地碳储量的一半，是大气碳储量的3倍[2−4]，在调节大气CO2浓度方面起到了重要的碳汇作用[5]。森林生态系统碳循环、碳固持、碳交换等研究已成为生态学研究的热点问题，因此，对森林生态系统进行的碳通量动态研究具有十分重要的科学意义。

大气CO2浓度的变化是导致全球变暖的重要因子，对大气CO2浓度长期监测可为人类制定碳排放及碳贸易政策提供理论基础[6]；通过大气CO2浓度的监测数据反演区域及全球的碳源−汇分布[7−8]；利用自由大气增浓实验（FACE）研究生态系统对CO2施肥效应的响应机制[9]。森林的生长代谢与大气CO2浓度变化密切相关。森林植被的光合作用与呼吸作用会影响CO2浓度[10]，而林冠层的CO2浓度的时空变化又对森林生态系统的生产力产生重要影响[11]。目前，与森林内部碳循环过程密切相关且作为涡度通量观测重要补充的森林生态系统CO2浓度的研究尚不多见[12]。陈步峰等[13]在海南尖峰岭热带雨林观测到林内CO2浓度最高，谭正洪等[12]得到近地层CO2浓度日尺度变化曲线为“双峰”的观测结果；焦振等[14]在帽儿山温带阔叶林观测到CO2浓度在日尺度上呈现“单峰”，且垂直方向上CO2浓度随高度增加而降低；陈晓峰等[15]利用CO2廓线系统在安吉毛竹林生长季的监测结果显示，CO2浓度在7月达到最低值，林冠层与林冠上层垂直梯度变化趋势相似。NPP模型对生态系统净初级生产力（Net Primary Productivity，NPP）模拟可以反应生产力的变化趋势。GROSSO等[16]利用超过5 600个站点的全球生态观测网数据集合（包括年均NPP、土地覆盖等级、降水和温度等），部分修正了Miami模型湿润区NPP高估的问题，提出了NCEAS（National Center for Ecological Analysis and Synthesis）模型。近年来，NCEAS模型已在全球/区域生态定量评估中验证了该模型的可靠性和准确性[17−18]。董丹等[19]利用改进CASA模型估算中国西南喀斯特地区1999—2003年的NPP总体呈现增长的趋势。李传华等[20]基于Biome-BGC模型，利用参数本地化对五道梁地区草地生态系统1961—2015年净初级生产力的情景模拟表明研究区的NPP呈显著上升的趋势。年尺度上的NPP模拟研究已取得较明确的结果，但由于植物生长对降水存在的滞后效应[21−22]使得模型对月尺度上的净初级生产力模拟具有不确定性，模拟结果与当月的降水及温度相关性低[23]。

基于生态系统过程模型的模拟研究已成为描述森林生态系统碳动态的重要研究手段。王媛等[24]基于EPPML模型千烟洲马尾松人工林的碳循环过程进行的模拟分析得到的结果与实测值十分吻合；韩其飞等[25]基于Biome-BGG模型对天山北坡的针叶林长达50年的模拟结果显示气温的升高促进了生态系统的NPP； PAULICK等[26]利用FORMIND模型模拟了两种地形不同的热带森林净生态系统交换量（Net Ecosystem Exchange，NEE），发现在演替早期的NEE最大，在演替后期的NEE接近零。然而对于可以深入了解森林碳循环过程的CO2通量动态模拟研究目前鲜有报道。热带森林通常指的是南北纬回归线之间具有典型天然热带植被的森林[27]，作为地球系统碳循环的组成部分，在缓解全球气候变化中扮演着重要角色，近年来，众多专家学者从热带森林土壤呼吸、天然林再生过程中的碳积累速率、森林恢复对地面碳密度复原等方面开展了研究[28−30]，但对于热带森林的CO2通量在森林碳循环过程的动态变化关注较少。鉴于此，笔者基于泰国SKR环境研究中心站点2001—2003年热带森林数据集，结合以往的研究结果[16−18]，选择NCEAS模型进行改进，对月尺度NPP进行预测，耦合Fick第一定律表征森林的CO2扩散过程，同时考虑温度变化对凋落物分解速率的影响，构建热带森林碳通量动态模型，旨在模拟热带雨林CO2通量月动态循环变化，并对模型中参数进行敏感性分析，确定各参数对森林碳通量影响程度，为进一步完善热带森林碳循环模拟研究提供技术及理论基础。

1 材料与方法

1.1 研究区概况泰国位于亚洲中南半岛的中南部，地形以山地与高原为主，地势北高南低，根据地理特征分为4个自然区域。北部丛林密布，东北部属于湄公河流域，为高原地带，干季少雨，雨季泥泞。中部为湄南河冲积平原，南部是狭长的丘陵地带。全国森林总面积为1 440万hm2，森林覆盖率达25%。

研究区位于泰国呵叻高原西北部Sakaerat环境研究中心（地理坐标14°29′32.5″N，101°54′58.7″E），Sakaerat研究区（以下简称为SKR）占地面积约为78 km2，地形较平缓，坡度不超过4°，属于热带草原气候，年平均降水量在1 200～1 300 mm，年均温24 ℃左右，干湿季分明。干季为11月至翌年3月，炎热干燥，月降水少于40 mm；雨季为4−10月，降雨充沛[31]。研究区主要包含29.5 km2的干旱常绿林区域和12.2 km2的干旱龙脑香科森林区域。植被类型是热带季节性常绿林，森林冠层高度为35～40 m。海拔高度约543 m，主要土壤类型是由砂岩母质演化而来的浅石质盐土。

1.2 数据与方法

1.2.1 数据来源与处理收集、整理SKR环境研究中心站点通量塔2001−03−01—2003−12−31综合气象要素和CO2浓度数据集（时间分辨率为30 min）。气象观测系统依4个不同垂向梯度（距离地面：5 、20 、35 、45 m）架设，同步测定温度和CO2浓度（图1）。翻斗式雨量计被固定安置于通量塔顶（47 m，远高于冠层高度），以消除冠层截留的影响。数据处理选取线性内插法对少量的缺测数据（缺测数据占比小于总体样本的5%）进行插补，以确保原始数据的完整性。最终，将逐半小时序列数据融合为逐月数据集，用于构建森林碳通量动态模型。

图1 涡度通量塔监测示意图Fig. 1 Schematic diagram of eddy flux tower for monitoring

1.2.2 GLUE算法1992年，BEVEN等[32]提出普适似然函数不确定性评估方法，即GLUE（Generalized Likelihood Uncertainty Estimation）算法，对模型的参数进行识别和率定。在预先设置好的模型参数范围内运用蒙特卡洛方法随机获取参数组合，代入模型中进行模拟，计算出模拟值与实测值的似然函数值，再计算函数值权重，得到参数组的似然值；根据预先设置好的似然值阈值，筛选出参数组似然值并重新归一化，分析参数不确定性[33]。

1.2.3 模型构建与评价本研究将森林生态系统作为一个整体进行模拟研究，关注影响林中CO2通量的各个组分，包括森林向大气边界层的扩散过程、植被的光合及呼吸作用、凋落物分解，并对各组分进行分析，将其在模型中重现，模拟月尺度下森林内部CO2通量的动态变化。

通过数据处理，将CO2浓度整理为通量数据，以林冠以下的3层CO2通量均值作为森林内部的CO2通量，以林冠上层梯度作为大气边界层CO2通量，运用Fick第一定律，模拟林内向大气的CO2扩散通量。净初级生产力是由植被以CO2为底物进行光合作用所产生的有机质总量中扣除植被自养呼吸后的剩余部分，因此可以表征光合与呼吸作用对CO2通量的影响。应用此种关系，改进后的NCEAS模型可利用SKR监测站的降水数据对月尺度森林NPP进行模拟。利用NCEAS模型模拟的月净初级生产力对凋落物质量进行约束，同时考虑温度变化对凋落物分解速率的影响，对森林凋落物的分解进行模拟。

根据图2对模型形式化，森林CO2通量变化量（dCO2）计算公式为：

图2 热带森林碳通量模型基本构成Fig. 2 Basic structure of carbon flux model of tropical forest

式中，Fd是凋落物分解CO2通量，NPPmonth为月净初级生产力，Fc森林内部与边界层之间的CO2扩散通量。

模型的循环迭代步骤为：1）设置迭代元素i，i值为6～36。2）判断i是否等于6，若i=6，则第一个月的CO2通量也就是6月的CO2通量数据等于5月的实测值加上6月CO2通量变化量dCO2。3）若i≠6，则当月CO2通量数据等于上一个月的通量模拟值加上当月的通量变化量dCO2。直到i=36，模型的迭代结束，否则，i=i+1，继续进行迭代计算。

运用GLUE算法，设置3个参数，在预先确定的参数范围内随机选取参数组的模拟（随机模拟次数为100 000次），似然函数选取均方根误差（Root Mean Squared Error，RMSE），进行参数敏感性分析，确定模型的最优参数组合，将最优参数组合的CO2通量模拟值与实测值比较，对模型优度评价。

1.2.4 模型各组分计算（1）引入植被对降水滞后响应的月尺度NPP模型优化。

NCEAS模型的NPP的计算公式为：

公式（2）为裸地计算公式，式中，ANPP单位为gC·（m2·a）−1。

公式（3）为非裸地计算公式，式中，ANPP单位为gC·（m2·a）−1。

式中，F（MAP）单位为mm。

式中，F（MAT）单位为℃。

ANPP为地上NPP，MAP是年降水量，MAT是年均温度，比较MAP及MAT模拟结果，取最低值为NPP。根据图3所示，F（MAP）显著低于F（MAT），因此，选取MAP对NPP进行模拟。模型可模拟年尺度的NPP，而本研究使用月尺度的NPP对森林的碳通量动态变化进行模拟，因此根据研究内容，结合研究区域气象状况，提出1种新方法，对NCEAS模型改进，从而模拟月尺度的NPP。

图3 2002和2003年SKR站热带森林NPP的2种模拟结果 Fig. 3 Simulation results of Tropical Forest NPP at SKR station in 2002 and 2003 NCEAS

由于植被生长对降水存在明显的滞后响应[34−35]。结合NCEAS模型的运算过程，将前期不同时间跨度月降水量（考虑前1−12 月）均值作为当月降水量，并假设全年的月降水量都相同，则每个月的净初级生产力也应当一致，将当月降水量乘以12，获得年降水量，代入到NCEAS模型中，获得在全年的月降水量都相同的假设情景下年NPP，再除以12，得到月NPP。以此类推，在假设全年每个月降水量都相同的情景下，利用NCEAS模型，模拟当月的NPP，计算公式为：

式中，NPPmonth为月净初级生产力，单位为 gC·m−2·月，PPT为不同时间跨度月降水量均值，单位为mm，统计1−12月不同尺度的NPPmonth，得到各月尺度模拟下的全年净初级生产力NPPm，与NCEAS模型利用年降水量计算的NPP对比，确定最优月尺度。

（2）CO2扩散通量Fc。模型假定森林内部与边界层之间的物质交换为平衡状态，主要受浓度扩散控制[36]，因此其CO2扩散通量Fc可以通过Fick第一定律进行估算，计算公式为：

式中，Fc是森林内部与边界层之间的CO2扩散通量，单位为 gC·m−2·月；k为森林内部与边界层之间的扩散系数，单位为 gC·m−2·月，用于表征扩散速度；dC/dx是边界层（45 m）与森林内部(5、20、35 m均值)的CO2浓度梯度。

（3）凋落物分解CO2通量Fd。模型使用一种简单，容易参数化的方程表达凋落物每月分解产生的CO2通量，计算公式为：

式中，Fd是每月凋落物分解产生的CO2通量，单位为 gC·m−2·月；M0是每月的凋落物质量，单位为 gC·m−2·月；rinst是由TOMCZYK[37]提出的考虑到温度变化对估计温度敏感性潜在影响的每日分解速率 d−1；Days表示各月包含的天数。

式中，rref是在当温度参考值Tref在288 K的分解速率0.020 4 d−1[38]。kB是玻尔兹曼常数（8.617×10−5eV·K−1）。Ea为活化能（eV）。

1.2.5 模型的参数取值范围碳通量模型需要率定的参数共有3个，见表1。

表1 各参数内涵及取值范围Tab. 1 Connotation and value range of each parameter

根据参数的物理特性，确定了k、Ea、M03个参数的初始取值范围，M0的取值范围通过最大月NPP进行约束。

1.2.6 模型优度评价选用皮尔森相关系数（Pearson Correlation Coefficient，r）与均方根误差（RMSE）对模型结果与实际测量值进行检验，r表示两个变量之间线性相关关系，在模型模拟精度分析中，用r表示测量值与模拟值的相关关系。RMSE是用于衡量观测值同模拟值之间的偏差，其值越低表明模型的预估结果越接近于实际测量值，即模型的拟合能力越高。

式中，CO2_sim是模型模拟值，CO2_ob是实际观测值，是模拟值的平均值，是观测值的平均值，模拟值与观测值之间的协方差与标准差的商为样本皮尔森相关系数r。

2 结果与分析

2.1 生态系统气象观测值由图4可见，研究区降水量的季节分布差异明显，降水主要集中分布在4−10月，11月至翌年3月为干季。2002年4−10月的降水量为1 720.76 mm，占全年降水量的98.9%，降水量最多的月份为5月，高达887.36 mm；2003年4−10月的降水量为1 109.31 mm，占全年降水量的96.2%，降水量最多的月份为5月，为355.66 mm。全年降水量与温度变化趋势呈现明显的“单峰”，雨热同期较为明显，是典型的热带草原气候，温度在4−5月达到高峰期，2002年年均温度为23.38 ℃，4月平均温度最高（5.65 ℃），1月平均温度最低（20.31）；2003年年均温度为22.88 ℃，4月平均温度最高（25.85 ℃），12月平均温度最低（18.54 ℃）。2002与2003年的降水与温度季节变化趋势基本一致。

图4 2001年3月至2003年12月的降水量及温度变化趋势Fig. 4 The variation trend of precipitation and temperature in the SKR monitoring station from March 2001 to December 2003

2.2 CO2浓度变化特征2002—2003年4个梯度的CO2浓度连续监测结果（图5）表明，SKR研究区的常绿林的CO2浓度特征呈现1个近U型的曲线，夏秋季低，冬春季高，2002年5 m层CO2浓度最 低 值 为6月 的378.71 μmol·mol−1，最 高 值 为12月的390.93 μmol·mol−1，2003年最底层5 m高CO2浓度最低值为8月的387.21 μmol·mol−1，最高值为12月的398.52 μmol·mol−1。垂直方向上，从最底层5 m的CO2浓度值到最高层45 m的CO2浓度值依次降低。

图5 2002年和2003年CO2浓度月变化趋势Fig. 5 The monthly variation trend of CO2 concentration in 2002 and 2003

2.3 最优月尺度计算假设全年月降水量都相同，考虑降雨的滞后效应，分析前期不同时间跨度的月尺度降水所模拟的年净初级生产力模拟结果（NPPm）与以全年降水量的年净初级生产力模拟结果（NPP）。

由图6-a可见，在2002年，1−8月尺度的降水量均值的NPPm存在较大差异，呈现线性增加的趋势，表明降水量对年净初级生产力具有显著影响。8−12月尺度的降水量均值的NPPm差异较小。比较2种年净初级生产力果，接近NPP 562.33 gC·(m2·a）−1的是前6个月降水量均值的NPPm，为531.82 gC·(m2·a）−1，前7个月降水量均值的NPPm，为531.82 gC·(m2·a）−1。2003年降水量相对2002年较低（图4），以不同时间跨度月尺度（1−12月）降水量均值作为当月降水量的NPPm呈现线性增加的趋势（图6-b），植被的生产力受到水分的限制。对比2种年净初级生产力模拟结果，接近NPP 438.82 gC·（m2·a）−1的是前6个月降水量均值的NPPm，为416.75 gC·（m2·a）−1，以上结果表明，每个月的净初级生产力除了与当月降水量相关，还与前期的降水量相关。

图6 2002年和2003年基于不同月尺度的NPP模拟结果Fig. 6 NPP simulation results based on different monthly scales in 2002 and 2003

2.4 模型拟合结果与检验经过GLUE算法的循环模拟后得到拟合效果最佳的CO2通量模拟值与实测值比较，结果见图7，模拟值与实测值的趋势具有较高的一致性，r=0.69，RMSE=14.93 gC·m−2·月。

图7 2001年6月至2003年12月森林内部CO2通量模拟值与实测值对比Fig. 7 Comparison of simulated and measured CO2 concentration in the forest from June 2001 to December 2003

从全年尺度来看，模型能准确模拟出森林在雨季与干季的CO2通量变化特征，雨季时的高温与丰富降水有利于植被的光合作用，森林的CO2通量降低。干季时的温度较低、辐射不足、降水稀少，不利于植物的生长，森林CO2通量达到全年的最高值。2002年，CO2通量模拟值与实测值的峰值均出现在12月，实测值为510.11 gC·m−2·月，模拟值为524.78 gC·m−2·月；2003年，CO2通量模拟值与实测值的峰值出现在12月，实测值为526.08 gC·m−2·月，模拟值为537.57 gC·m−2·月。根据模拟结果来看，实测数据与模拟数据的相关性较强，数值大小接近，总体而言，模型能够较好地模拟泰国热带雨林的CO2通量数值大小及变化趋势。

2.5 参数敏感性分析选取100 000组参数进行不确定性分析，比较不同参数在模型的敏感性大小（图8）。结果表明，随着取值的不同，M0和k两个参数的拟合优度存在较为明显的差别，表明参数在模型中较为敏感。其中参数k最为敏感，在参数分布区间可看出较为明显的极值区间（k取值0.4时）；M0在区间0～13.9呈现下降趋势，但在M0>13.9的区间内都具有较低似然值取值；而参数Ea在整个参数分布区都有低似然值取值，表明Ea不敏感。

图8 模型的3个参数与似然函数关系散点图Fig. 8 Scatter plot of the relationship between three parameters of the model and likelihood function

3 讨 论

本研究对NCEAS模型进行了改进，并基于构建的热带森林碳通量动态模型，应用泰国SKR监测站的降水量及温度数据，CO2涡度通量数据，对森林内部的CO2通量变化状况进行模拟，并与监测站的实际观测结果比较，得到的结果分析如下：

（1）统计结果表明，泰国SKR监测站热带森林干季与雨季的降水量差异较明显，雨季降水量占全年的90%以上，显著高于干季，并且月均温度也高于干季；CO2浓度垂直梯度格局上受土壤呼吸与冠层碳代谢等因子影响，在垂直方向上随高度的增加而降低，造成该现象的原因是由于夏秋季节气温的上升促进了森林生态系统中木本植物的生长[14]，增强了植物群落的光合作用，导致CO2浓度的下降。不同高度间CO2浓度变化趋势无显著差异

（2）NCEAS模型一般通过利用气候和土地覆盖的数据来估计全球年尺度上的NPP，本研究尝试改进NCEAS年尺度NPP模型，对月尺度上的森林净初级生产力进行模拟，并与原模型年尺度的模拟结果进行比较。结果表明，模型在使用降水驱动时，NPP值明显低于气温驱动下的NPP值，说明降水是限制该区域植被生长最重要因子，符合利用NCEAS模型估算全球NPP，得到低纬度地区NPP的变化主要受降水量驱动的结论[30]。基于植物滞后响应的特点，从植被对降水的滞后时间跨度着手，并且考虑了模型适用范围存在的差异。假设全年月降水量都相同，利用NCEAS模型采用降水量计算NPP，通过与年降水量模拟结果对比确定了前6个月降水均值为影响当月生产力的最优时间尺度，提高了NCEAS模型在时间尺度上的通用性与准确性，表明了对植物的干旱状况评估时，应当考虑不同时间尺度下降水积累效应及水分的亏缺对植被的影响[39]，此结果也与崔林丽等[40]分析中国东部植物NDVI对降水响应特征，得到植被对降水变化最大响应的滞后期随着纬度降低由北到南逐渐延长的结果相符。

（3）构建的模型对SKR森林的CO2通量的模拟结果与实际测量值的对比说明新模型在拟合热带森林在碳通量方面具有较好的适用性、可靠性。参数敏感性分析结果表明k是影响热带森林CO2通量变化的主要因子，当通过生态系统过程模拟森林碳通量变化时，要着重考虑扩散过程，同时耦合其他生态系统功能。模型的模拟值与实测值之间相关性好，相关系数达0.69。模拟的CO2通量数据与实测CO2通量数据的季节性差异基本一致，均呈现干季高、雨季低，与同在热带地区的海南尖峰岭通量观测结果一致[13]。在研究区同一经度带（101°E）的西双版纳热带季雨林的研究中，张一平等[41]发现林冠−大气间的日均值呈现明显的季节性变化特征，干季的碳通量绝对值高于雨季，但在更高的纬度地区，赵仲辉等[42]在亚热带杉木人工林生态系统的开展的研究结果与张元媛等[43]在寒带峨眉冷杉成熟林生态系统的研究结果却呈现CO2通量冬季低、夏季高的现象，与本研究的结果相反，这可能是因为热带地区生长旺季冠层的CO2浓度由于强烈的光合作用而变得很低[12]，导致在不同温度带的观测出现了不同的结果。通量数据质量是影响模型误差的重要因子，并且易受降水、大气湍流及传感器自身噪声影响[44]。2002年雨季的CO2通量模拟值出现了较大的波动，雨季的气象状况对模拟结果存在明显影响，没有很好体现雨季的CO2通量变化特征。此外，模型的误差与参数的选择有重要关系，碳通量动态模型根据生态系统过程输入了3个参数，以RMSE作为似然函数，对参数敏感性进行分析，发现模型对参数Ea和M0不敏感，而当k取值0.4时，取得极值，意味扩散系数k影响热带森林CO2通量变化的主要因子，并说明若只受扩散作用影响，森林生态系统的CO2每个月将会有40%扩散到大气中，经过2.5个月后扩散完成。表明了在代谢旺盛的热带森林，向外界扩散的CO2是森林CO2通量变化的主要原因，这与基于SVAT模型在印度尼西亚苏拉威热带森林开展的冠层与大气的CO2交换研究结果相似[45]。

尽管构建的模型模拟结果与实际观测值的趋势具有较高的一致性，但与通量观测结果相比较还存在一定的偏差，原因可能在于：1）模型所用数据只来自单一气象站点，结果缺少普适性，后续工作需要整合更多站点数据并进行模拟，以提高模型的适用范围；2）通量塔的最高层的CO2气体分析仪只有45 m，作为边界层的高度有所不足，但更高区域的CO2通量数据难以获取；3）构建的碳通量模型中只设置了3个参数，并未使用遥感NPP数据对模型进行验证；4）森林生态系统CO2通量还受其他因素的影响，比如土壤呼吸。虽然存在以上问题，但模型总体上能够较好地模拟SKR监测站热带森林CO2通量的数值大小和变化趋势。因此，在未来的研究中，为进一步提高模型的精确性，减小误差，深入揭示森林碳通量的动态变化过程，还需选取不同似然函数对模型参数率定，增加土壤呼吸参数。研究结果可为典型热带森林生态系统碳通量动态模拟及森林碳循环研究提供科学依据与理论基础。