陆地生态系统是一个巨大的碳库且对气候变化极其敏感，在全球碳循环中起到重要的作用^[1]。森林生态系统是陆地生态系统的重要组成部分，其每年碳固定量约占陆地生态系统植被碳固定量的2/3^[2]，通过物理、化学和生物过程调节地表和大气之间碳循环^[3]。近年来，气候变化变得异常，全球降水区域差异性明显以及亚热带和热带区域发生干旱的频率和强度在不断增加^{[1, 4]}。另一方面，在太阳活动等自然外强迫作用下，全球增温开始出现减缓甚至停滞的趋势。在此复杂的气候变化背景下，森林生态系统碳收支的时空变化趋势及其与环境和生物因子之间的关联性存在极大的不确定性。

环境和生物因子是森林生态系统与大气之间碳通量的主要驱动因子^{[5, 6]}，不同时间尺度上，对森林生态系统碳通量变化起决定性的因素不同。在秒到小时尺度，大气湍流是碳通量变化的主要因素；在小时到月尺度，气象因子决定碳通量的变化，如气温、湿度和太阳辐射；在月到年尺度，季节性天气变化以及植被物候变化是碳通量变化的主要影响因素^[7]。目前，已有大量研究分析了森林生态系统总初级生产力(GPP)及净初级生产力(NPP)年际变化与环境和生物因素的关联性，如与年均温度正相关关系^{[8, 9]}、年均温度负相关关系^[10]、降水量负相关关系^[10]、降水量正相关关系^{[9, 11, 12]}、夏季日均最大气温负相关关系^[13]、降水量和温度协同影响^[14]、厄尔尼诺与南方涛动的年代际变化影响^[15]。叶面积指数和温度分别控制着长白山阔叶红松林生态系统GPP的季节动态和年际差别^[16]。研究表明森林生态系统GPP和NPP年际变化对环境因子的响应方式存在不一致性，没有特定的机理来解释GPP和NPP年际变化^{[17, 18]}。

竹林生态系统是森林生态系统的重要组成部分。毛竹(Phyllostachys edulis)是中国的主要竹种，其面积约占中国竹林面积的70%，占世界毛竹林总面积的80%^[19]。与其他森林类型相比，毛竹在生长和经营管理方面具有特殊性，是竹产区林农经济收入主要来源^{[20, 21]}。国内外学者开展了大量毛竹林生态系统碳储量研究，结果表明毛竹林具有巨大的固碳和蒸腾潜力^{[22, 23, 24, 25, 26, 27]}，在固碳减排和缓解气候变化方面起到重要的贡献^[28]。安吉县毛竹林地上碳储量在1986—2008年期间呈增加趋势^[29]。福建省毛竹林生态系统碳储量在1985—2006 年期间也呈上升趋势，年降雨量与碳储量具有正相关关系，年均温与碳储量呈负相关关系^[30]。但是针对毛竹林生态系统总初级生产力的长时间动态变化及其与环境和生物因子的关联性分析还很缺乏。

本研究基于涡度相关技术获取的毛竹林生态系统GPP和MODIS数据以及MERRA(The Modern Era Retrospective-analysis for Research and Applications)再分析数据，采用结合涡度相关技术的光能利用率模型(Eddy Covariance Light Use Efficiency，EC-LUE)估算2004—2011年安吉县毛竹林生态系统GPP，分析安吉县毛竹林生态系统GPP变化趋势及其驱动因素，为了解和预测未来气候变化对毛竹林生态系统碳通量的影响提供理论依据。

安吉县位于浙江省西北部(119°14′—119°52′E，30°23′—30°53′N)，全县面积1886.45 km²。气候属亚热带海洋性季风气候，总特征为光照充足、气候温和、雨量充沛、四季分明，年均降水量1400 mm，年均气温15.6 ℃。森林覆盖率达到69.4%，拥有山林13.2×10⁴ hm²，其中竹林面积6.33×10⁴ hm²。毛竹林面积为4.99×10⁴ hm²，占林地总面积的37.8%，为全国著名的“中国竹乡”。安吉县毛竹林具有明显的大小年交替特征，逢奇数年东部区域为大年，西部区域为小年，逢偶数年则相反^[20]。毛竹林通量观测塔位于安吉县东南部的山川与天荒坪两个乡镇交界处。以观测塔为中心1000 m范围内主要森林类型为毛竹林，其平均冠层高度为11 m，林分平均密度为每公顷3235株。此外，还有少量针阔混交林、农田、城镇和道路等地类。通量观测塔介绍详见文献^[31]。

借助EdiRe(http://www.geos.ed. ac.uk/abs/ research/ micromet/ EdiRe/)和Matlab 7.0软件对2011年通量塔观测的原始数据进行数据处理，主要包括以下3个方面^[32]：通量数据校正、缺失数据插补和组分分解。通量数据校正包括二次坐标轴旋转^[33]、空气密度脉动影响校正^[34]、储存项校正^[16]、异常值剔除^{[35, 36]}和摩擦风速校正。采用基于夜间数据的非线性回归法(Lloyd & Taylor和Michaelis-Menten模型)进行缺失数据插补和组分分解^[37]。

采用通量塔观测数据计算的GPP对EC-LUE模型进行参数化并估算2004—2011年安吉县毛竹林GPP空间分布。EC-LUE 模型属于光能利用率模型，它的一个独特之处在于采用蒸发系数(EF)作为湿度对最大光能利用率进行调节，将EF引入模型后使碳通量和水通量得到耦合，提高了模型预测精度^[23]。EC-LUE模型的输入变量主要包括归一化植被指数(NDVI)、光合有效辐射(PAR)、温度(T)和波文比^[38]。Xu等^[27]参照光响应曲线模型将GPP饱和效应引入EC-LUE模型，提高了预测精度，其模型表达式见式(1)。

点击浏览公式

(1)环境因子 从MERRA再分析数据集(http://gmao.gsfc.nasa.gov/ research/merra/)中获取2004—2011年环境因子，包括太阳入射短波辐射、净向下短波辐射、净向下长波辐射、比湿、10 m处的气温和大气压。大气压的空间分辨率为1.25°×1.25°，其它数据空间分辨率为纬度0.5°×经度0.667°，时间分辨率为天(d)。

(2)生物因子 从美国橡树岭国家实验室(http://daac.ornl.gov/ MODIS/modis.shtml)获取2004—2011年1 km空间分辨率的8 d合成MODIS叶面积指数(Leaf area index,LAI)(MOD15A2)和250 m空间分辨率16 d合成的 MODIS NDVI(MOD13Q1)产品数据。根据质量控制值来判断LAI/NDVI数据是否可靠。采用上边缘平滑方法对MODIS LAI/NDVI 产品做平滑处理，修正不可靠数据^[39]。

采用2008年Landsat TM(Thematic Mapper)影像band1、band3、band4、band5 和NDVI作为特征变量，基于后向传播神经网络提取安吉县毛竹林信息，分类总体精度和Kappa系数分别为91.8%和0.89。根据样地调查实测数据精度评价得出毛竹林的分类精度为82.19%，达到应用的要求。为了与MODIS数据空间尺度相一致，需要对30 m空间分辨率的分类图进行尺度上推处理。采用移动窗口平均法将30 m × 30 m分辨率分别上推至250 m × 250 m(MOD13Q1)和1 km × 1 km(MOD15A2)分辨率。首先，将TM影像分类图中毛竹林设为1，其他地类为0。其次，分别计算8×8像元窗口(250 m × 250 m)和33 × 33像元窗口(1 km × 1 km)内毛竹林面积百分比。最后，为了确保尺度上推后的像元内主要地类为毛竹林，将毛竹林面积百分比大于80%的像元设定为1，其他为0。

采用EC-LUE模型模拟2004—2011年安吉县毛竹林生态系统GPP空间分布(图 1)。GPP在空间分布上表现为偶数年(如2004)东部区域(小年)高于(P<0.01)西部区域(大年)，奇数年(如2005)西部区域(小年)稍高于(P=0.26)东部区域(大年)(图 2)，即小年毛竹林生态系统GPP稍高于大年GPP。2004—2011年，安吉县毛竹林年日均GPP分别为6.15、6.19、5.89、6.12、5.67、5.95、5.95和5.82 g C m^-2 W^-1。单样本t检验表明安吉毛竹林年均 GPP在年际间不存在显著差异(P=0.66)。根据毛竹林GPP空间分布统计出安吉县东西部毛竹林年日均GPP，用于后续GPP年际变化及其驱动因素分析。

图1 基于EC-LUE模型的安吉县毛竹林GPP空间分布图 Fig.1 Spatial distribution map of GPP for Moso bamboo in Anji County based on EC-LUE model

图选项

图2 东部和西部毛竹林年日均GPP对比图 Fig.2 Comparison of annual average daily GPP between eastern and western areas

图选项

根据2004—2011年日均GPP与LAI、T和PAR的相关性分析，年尺度GPP与温度的相关性最高，但未达到0.05显著性水平(表 1)，温度是驱动年尺度GPP变化的主要因素。2004—2011年日均GPP呈下降趋势，东部、西部和整个安吉县毛竹林年日均GPP变化速率分别为-0.064、-0.033和-0.045 g C m^-2 W^-1，东部毛竹林年日均GPP变化大于西部(图 3)。从年均温度和GPP变化趋势分析得出年均温度持续降低是导致年尺度GPP下降的主要驱动因素(图 3)。

图3 2004—2011年日均GPP和温度变化趋势及其关联性 Fig.3 Change trends of annual average daily GPP and temperature and their relationship from 2004 to 2011

图选项

GPP年际变化(年际差值)与LAI、T和PAR年际变化的相关性分析表明，安吉县东部和西部毛竹林LAI年际变化是GPP年际变化的主要驱动因素，整个安吉县温度和LAI年际变化对GPP年际变化起主要驱动作用，但均未达到0.05显著性水平(表 2)。东部和西部毛竹林GPP年际变化的变异系数分别为54%和70%，西部毛竹林GPP年际变化幅度大于东部，这主要是由于西部环境因素和生物因素的驱动效应具有正相关关系，两者相互促进使得GPP年际变化幅度变大；然而东部环境因素和生物因素的驱动效应为负相关关系，两者之间相互抵消使得GPP年际变化幅度变小(图 4)。

图4 安吉县东部和西部LAI年际变化与温度关系图 Fig.4 Relationships between interannual variability of LAI and temperature in eastern and western Anji county

图选项

安吉县毛竹林生态系统年尺度GPP在空间分布上东部和西部具有差异性，这主要是由于东部和西部毛竹林具有相反的大小年交替规律以及大小年毛竹林生理生态特征具有明显差异引起的。大年和小年毛竹林主要特征差异为：(1)采伐经营措施通常在大年进行，导致LAI降低；(2)小年毛竹林的竹叶为1龄叶，而大年除新竹竹叶外，其余为老龄叶。竹叶的光合作用能力随年龄增加而降低^[40]，新叶的光合速率是老叶的3倍以上^[41]，因此小年竹叶具有更强的光合能力；(3)小年竹林的叶片在冬季仍保持鲜绿，能储存大量的营养物质，通过根系为次年新竹生长提供营养^{[40, 42, 43]}。从MODIS LAI产品数据分析中也证实小年的LAI确实要稍大于大年的LAI。综上所述，小年毛竹林GPP高于大年的主要原因是小年具有更高的有效叶面积和叶片光合能力。

年均温度是2004—2011年安吉县毛竹林年日均GPP呈下降趋势的主要驱动因素。温度是影响植被光合效率三大环境因素之一。毛竹林光合效率的最适合温度在20—30 ℃之间，随着温度降低，其光合效率呈下降趋势^{[44, 45]}。政府间气候变化专门委员会第五次报告指出近15年来(1998—2012)全球温度上升出现减缓和停顿，主要原因是海洋表层温度上升的减弱和太阳活动的显著减弱等自然因素，该结论与根据MERRA数据得出的2004—2011年均温度呈下降趋势相一致。因此，该研究得出年均温度下降导致毛竹林年日均GPP下降的结论是有依据的。已有研究表明森林生态系统GPP随年均温度呈线性增加，而呼吸碳通量呈指数上升^[8]。年均温度降低时，呼吸碳通量和GPP同时降低，对净生态系统碳固定量的影响取决于两者的综合效应。因此，2004—2011年安吉县毛竹林GPP呈下降趋势并不能表明近年来毛竹林生态系统净固碳量的下降。

LAI是毛竹林GPP年际变化的主要驱动因素，该结论与其他森林类型一致^[46]。东部和西部区域LAI年际变化对GPP年际变化的贡献要高于温度，主要原因是大小年生理生态特征差异引起年际间LAI具有明显的规则变化规律。对于整个安吉县LAI和温度年际变化对GPP年际变化起到相当的贡献，这可能是整个安吉县东部和西部大小年相反均衡了LAI年际变化对GPP年际变化的影响。LAI和温度分别控制着长白山阔叶红松林生态系统GPP的季节动态和年际差别^[16]。LAI和温度对毛竹林GPP年际变化的作用方向决定了GPP年际变化波动幅度，当两者作用方向一致(同时为正或负效应)时，GPP年际变化变异性较大。东部和西部大面积的大年或小年毛竹林中夹杂着小部分大小年相反的毛竹林^[20]，另一方面，东部与西部大小年分割界线难以确定，本研究根据调查样地新竹量来确定东部和西部界线难免会造成大年竹和小年竹混在一起，使研究结果存在不确定性。自然干扰在很大程度上影响GPP年际变化^[47]，因此2008年的冰雪灾害以及2011年春季干旱对GPP年际变化也造成一定的影响^{[27, 48]}。此外，由于MERRA数据的空间分辨率较低，难以反映小尺度区域环境因子变化，对分析结果也带来不确定性。

(1)毛竹林具有大小年交替的生理特征，由于大小年之间叶面积指数、叶片光合参数及其光合效率的差异造成了大小年之间GPP存在差异，同时也是GPP年际变化的主要驱动因素，因此深入了解该生理生态规律对毛竹林生态系统固碳量动态变化的影响具有重要意义。

(2)近15年来，由于自然因素原因，全球气温呈减缓甚至下降的趋势，使得安吉县毛竹林生态系统GPP呈下降趋势，该结论为解释毛竹林固碳能力和竹材产量的变化以及了解短期气候变化对GPP影响提供了参考依据。

(3)GPP年际变化的波动幅度受环境和生物因素共同影响，当两者作用方向相同时，GPP年际变化的波动幅度越大，该结论为解释毛竹林生态系统碳收支年际变化的波动性提供了理论基础。在今后研究中，还需进一步定量分析毛竹林大小年生理生态特征对其固碳能力影响以及环境和生物因素对GPP年际变化的贡献。