黑龙江省森林大气CO2施肥效应的时空变化趋势1)

莫祝坤 李明泽 王斌

(东北林业大学，哈尔滨，150040)

自工业革命以来，由于煤炭、石油等化石燃料焚烧和人类活动对森林资源的破坏，大量的温室气体(CO2、CH4、O3等)被排放到大气中。其中，大气CO2作为最主要的温室气体，其体积分数已从工业前期的227×10-6上升到目前的411×10-6[1]。尤其在20世纪60年代以后，大气CO2上升速率明显加快。有研究显示，在过去的60 a中，陆地生态系统的碳汇量已占全球碳汇的一半以上，且越来越多的证据表明，北半球中、高纬度森林生态系统起着巨大的碳汇作用[2-5]。一些学者将部分陆地碳汇的增长归因于大气CO2体积分数的快速上升所引起的植被生产力的提升，并将这一影响称之为大气CO2的“施肥效应”[6-7]。

理论上，大气CO2体积分数上升对植被生长的施肥效应可以体现在直接施肥和间接施肥2个方面[8]：首先，外界环境中更高CO2体积分数将直接增加绿色植被叶片内的CO2体积分数以及CO2与O2的比值，这有利于叶片内化学酶进行羧化作用，进而提高植被的光合作用强度并抑制呼吸作用，即增长植被的净光合速率，从而促进植被的生长[9]，这就是CO2的直接施肥效应。另一方面，大气环境中不断增长的CO2体积分数可能会降低叶片气孔的开放程度，从而减少由于蒸腾作用所引起的水分的流失，增长植被体内的水分利用率。因此，这种CO2提升的间接影响可以延长季节性干旱生态系统的生长期，从而可能刺激植物体内生物量的积累[10]，这称为CO2的间接施肥效应。最近几十年，为了证实现实环境中大气CO2施肥效应假设的真实性及其时空变化趋势，一系列相关检测方法被提出，并在不同时空尺度上得以应用[11-13]。然而，这些研究并没有得到一致性的结论，现阶段，关于大气CO2是否对植被生长存在施肥效应仍存在较大的争论和不确定性[14-16]。

自然条件下的树木生长轮能够提供树木长期的生长记录，是记录环境信息的天然档案。Cook[17]最早将任何一年树木的径向增长分解为年龄因素，气候变化因素，CO2施肥效应因素，林分内部和外部干扰因素影响的集合。为了探究环境对树木的影响因素，首先必须消除林龄相关的生长趋势，并尽可能地保留类似大气CO2施肥效应这种低频信号。在树木生长轮分析方法中，国内外学者以往都是通过传统的保守去趋势法[18]，区域曲线标准化法[19]，断面积修正法[20]和径级法[21]等去除年龄趋势后再探究气候对树木径向生长的影响。目前，已有很多学者利用自然条件下的树木生长轮数据来检测CO2的施肥信号，并提出了不同观点[13，22]。其中许多学者已经在不同区域利用自然生长的树木生长轮数据通过各种方法证实了大气CO2施肥效应这一假设，如Wang et al.[23]通过对比加拿大白云杉(Piceaglauca)在不同CO2体积分数下的生长轮指数发现，处于相同生长阶段的树木在高CO2体积分数下的生长速率要明显高于工业前期CO2体积分数相对较低时的生长速率，以此说明大气CO2施肥效应的存在。如Soule et al.[24]通过对生长轮数据进行分析发现，在1950年以后大气CO2体积分数升高的条件下，树木的径向增长速率显著增加，特别是在干旱年份，且一般在水资源最有限的地点增幅最大。然而，还有一部分研究则通过对树木生长轮的分析认为，大气CO2体积分数的增加对树木的生长不存在明显的施肥效应，或者现有的检测手段很难检测到施肥信号，如Arrigo et al.[25]通过对美国北部白云杉生长轮数据进行检测，认为大气CO2可能需要达到一定的阈值才会对树木生长产生施肥效应，而目前的大气CO2体积分数远远没有达到这个阈值，树木径向的快速增长是因为温度的上升而并非大气CO2的施肥作用。目前,利用自然条件下的树木生长轮检测大气CO2施肥效应研究尚无定论，且各种检测途径仍存在一定的争议和不足之处。相比于传统的去趋势方法，广义加性混合模型法(GAMMs)是分析树木年龄等因子影响的一个更加有效工具[26]。该方法是在断面积修正法的基础上，将树木每年断面积增量看作为胸径、生理年龄、树高等生理或立地条件指标的加和函数，通过平滑函数来拟合断面积增量随年龄、胸径大小等因子变化的趋势。另有研究表明，相比于传统的保守去趋势法，GAMMs模型能够更多地保留类似大气CO2施肥这种低频的信号[13]。

除了利用树木生长轮数据来检测大气CO2的施肥效应外，基于光的利用效率(ELU)模型和遥感卫星影像所驱动的总初级产生力(PGP)数据集也是一种分析植被生长对气候变化响应的有效工具[27-28]。然而，这类数据很少被用于大气CO2施肥效应的检测研究中，这是因为目前大部分ELU模型主要由温度、降水和太阳辐射等因素驱动，而对大气CO2的影响仅考虑了与叶面积密切相关的光合有效辐射观测值的变化来影响植被生产力，却忽略了大气CO2体积分数上升可以使光的利用效率增加，这就导致了所估算的PGP明显低估了大气CO2的施肥效应[29]。最近，一种经过通量站点观测数据修正后的ELU模型所估算得到的全球PGP长期数据集(EC-LUE GPP)很好地弥补了上述问题[30]，该数据集通过一个尺度因子有效整合了大气CO2体积分数、太阳辐射分量和水蒸气压差对植被生长的影响，这使得利用该数据集检测大气CO2施肥效应的时空变化趋势成为了可能。如Wang et al.[28]利用修正后的EC-LUE GPP数据分析发现，CO2体积分数的上升以及大气氮沉降等都会对植被的光合作用产生显著的影响，但最近几十年来全球范围内大气CO2施肥对植被光合作用的影响呈下降趋势。

综上所述，大气CO2施肥效应的检测一直是全球碳循环研究的热点问题，但目前有关大气CO2施肥效应的检测及其时空变化趋势的研究结果不一致，尚存在较大的不确定性，且各种检测途径仍然存在一定的不足之处。黑龙江省森林一直是全球碳循环研究的重点区域，该区域是否存在大气CO2的显著施肥效应将直接影响到中国乃至全球陆地碳收支的预测和评价。然而，目前针对该区域大气CO2施肥效应的研究却鲜有报道[31]。因此，本研究以黑龙江省森林为研究对象，旨在利用树木生长轮观测资料和遥感数据协同检测黑龙江省森林大气CO2施肥效应并进一步分析其时空变化趋势，进而对大气CO2的施肥效应假设的真实性做出准确评价以及为明确黑龙江省森林在全球碳循环中所扮演的角色提供科学的依据。

1 数据来源

本研究的黑龙江省森林林区是全国重点林区之一，全省森林面积2 145万hm2，约占全国的1/10，森林蓄积量为20.52亿m3，森林覆盖率为46.14%。

1.1 生长轮数据

树干作为树木最大的生物质碳库，其碳库的年际变化可能对大气CO2的上升存在响应。本研究选取黑龙江省森林中最具有代表性的树种红松(Pinuskoraiensis)的生长轮样芯作为地面观测数据集，根据第九次全国森林调查报告，该树种在研究区内蓄积量占比较大，而且分布广泛，能够较准确地记录研究区内森林的环境信息。同时基于国际树木生长轮库标准，在野外随机选取受干扰较小且生长健康的树木，利用生长锥在树干1.3 m处沿水平坡方向分别取两根样芯。将采集好的样芯按照树木生长轮样本的基本程序进行预处理，包括风干、打磨等。生长轮样品经预处理后，在双筒显微镜下进行目视交叉定年，然后用VELMAX生长轮仪(型号A60)通过手摇每个生长轮的位移在显微镜下判断的方式进行生长轮宽度的测量，测量精度为0.001 mm。对于生长轮不清晰的样芯采用浸湿的方法以更精确地辨认生长轮界限，测量后用COFECHA程序[32](版本6.06)对交叉定年和宽度测量的结果进行检验，根据检验报告的提示，剔除和修正不准确的定年，以确保每个生长轮有准确的日历年代，最终收集了9个地区288棵树的生长轮样芯(共576根样芯)。各生长轮样芯的采集地区基本信息如表1所示。

表1 生长轮采样地区基本信息

1.2 GPP数据集

本研究所采用的EC-LUE GPP数据(https://doi.org/10.6084/m9.figshare.8942336.v1)是以全球涡度相关碳通量观测资料为基础并对光能利用率模型修正后得到的，该数据模型的驱动变量仅为：归一化植被指数、光合有效辐射、空气温度和波文比。利用全球范围近百个涡度相关站点的测量数据验证表明，该数据可以解释超过70%的GPP数据变化，能够很好的刻画GPP的季节和年际变化，其模拟能力超过了MODIS-GPP产品，同时该数据集在传统光能利用率模型中考虑了大气CO2施肥效应对光能利用率的影响。

1.3 气象及大气数据

每月气象数据采用CRU(Climatic Research Unit)气候数据集，版本为4.0.4。包括每月平均温度(Tm)，每月总降水(Pre)、每月真实水蒸气压(PVA)和水蒸气压差(DVP)(https://data.ceda.ac.uk/badc/cru/data/cru_ts)。其中每月水蒸气压差(DVP)是根据每月饱和水蒸气压(PSV)与PVA的差值计算而来，每月DVP的计算方法如公式(1)～(4)所示[33]。

DVP=PSV-PVA；

(1)

PSV=6.112fwe17.67Tm/Tm+243.5；

(2)

fw=1+7×10-4+3.46×10-6Pmst；

(3)

(4)

公式(3)与(4)中，Pmst为海平面气压，Pmsl为平均海平面气压(1 013.25 hPa)，Z为海拔高度。

大气CO2体积分数数据来源于美国夏威夷的冒纳罗亚天文台(Mauna Loa observatory, Hawaii)(https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html)。

全球大气氮和磷沉降数据(Ndep、Pdep)来源于全球气溶胶化学模型(LMDZ-INCA)的模拟数据(https://figshare.com/s/898e5d6c4f0982674ae3) ，该数据集能较好地评估氮和磷沉降的变化对森林碳吸收的影响[34]。

2 研究方法

2.1 线性混合效应模型检测大气CO2施肥效应

为了探索气候对树木径向生长的影响，首先必须消除林龄相关的生长趋势，并尽可能地保留类似大气CO2施肥效应这种低频信号。因此，本研究采用广义加性混合模型去除树木自身年龄以及胸径大小对树木径向生长的影响，首先需要将每年测得的树木生长轮宽度转换为树干断面积增量(IBA)，IBA通常表现为在树木幼龄期呈指数形式增长，树木成熟时会逐渐趋于稳定状态[26]。

(5)

公式(5)中：IBA为树干断面积增量，Rt和Rt-1分别为从树木生长到第t年和第t-1年的树干生长轮宽度累积之和。由于树干断面积增量呈偏态分布，需要对树干断面积增量进行对数转换，再通过GAMMs模型拟合每个研究样地上每棵树木年龄和胸径对IBA的影响。在GAMMs模型中，以胸径(HDB)和树木形成层年龄(Age)为固定因子，以树木号为随机因子，对每一个红松样地进行单独建模分析，模型形式如下：

ln(IBAi,t+1)=s(HDBi,t)+s(Agei,t)+Zi,tBi,t+vi,t+εi,t。

(6)

公式(6)中，i代表树木号，t代表年份，BZ为随机效应，v为具有一阶自回归结构的误差项，ε为残差项，s为平滑函数，采用3次样条平滑函数来进行平滑处理。

在去除树木年龄与胸径大小的影响后，使用线性混合效应模型分析了1960—2015年间，每个研究样地上树木生长轮剩余生长信号中气候因子以及大气CO2体积分数对树木径向生长的影响。在线性混合效应模型中，以每年大气CO2体积分数、树木生长前1年秋季平均气温、树木生长当年春季平均气温、夏季平均气温、夏季总降水以及寒冷季节降雪量为固定因子[13]，以树木号为随机因子，对每一个样地进行单独建模分析，模型形式如公式(7)所示：

(7)

公式(7)中，IBAIres为IBA的残差，代表去除树木年龄和胸径大小等影响后的树木剩余生长信号，An为模型的第n个参数，Clim为气候变量，k为气候变量个数，v为具有一阶自回归结构的误差项，ε为残差项。在线性混合效应模型中，根据修正后的赤池信息量准则，对不同气候因子可以生成的所有潜在模型进行排序。在所有排序模型中，计算每个模型的权重值，累计权重在95%以内的模型被认为是最优模型[13]，之后将选择出来的最优模型进行模型平均。本研究通过参数t检验来说明大气CO2体积分数以及各气候因子对树木剩余生长信号的影响。如果通过统计检验证明模型中CO2变量对树木剩余生长信号(BAIres)具有显著的积极影响，则可说明大气CO2施肥效应真实存在。

2.2 随机森林模型定量大气CO2施肥效应的时空变化趋势

本研究采用随机森林算法来进行历史观测变量PGP的重建并估算植被生产力的变化对大气CO2体积分数增加的响应，即大气CO2施肥效应的强度(β)[33]。研究以修正后的每年生长季(5—9月)的EC-LUE GPP作为观测变量，以每年大气CO2平均体积分数、生长季的平均温度和总降水、水蒸气压差、大气氮沉降以及磷沉降等作为影响因子，分别重建黑龙江省森林每个像元上的PGP数据，研究时段为1982—2015年，数据空间分辨率均重采样到0.5°×0.5°。研究设计了两种大气CO2体积分数变化情境模式，分别进行历史观测变量PGP的重建，第一种是所有驱动因子随时间真实变化，第二种则保持大气CO2体积分数在1982年的水平不变，其他驱动因子随时间真实变化[27，33]。之后利用两种情境下随机森林算法模拟得到的估计值的差值来拟合β值，具体公式如下：

ΔPGP=β×Δφ(CO2)+ε。

(8)

公式(8)中，ΔPGP代表利用随机森林算法所模拟得到的观测变量在不同CO2体积分数变化情境下的估计值的差值，Δφ(CO2)为两种情境下大气CO2体积分数的差值，ε为误差项。之后对所估算的β值进行归一化处理，通过设置一个15 a移动窗口来检测黑龙江省1982—2015年间每个像元β值的变化趋势(例如:1982—1996，1983—1997，…，2001—2015)，并采用非参数趋势检验(Man-Kendall检验)估计β时间上的变化趋势和统计学显著性[27]。

3 结果与分析

3.1 树木生长轮途径结果与分析

本研究通过GAMMs模型拟合了树干断面积增量随树木年龄以及胸径大小的变化趋势，表2给出了9个红松样地GAMMs模型的拟合结果。通过对各个样地GAMMs模型进行非参数项显著性检验发现P值均小于0.01，说明树木年龄以及树木胸径大小对IBA均有显著的影响。同时，各样地GAMMs模型的修正决定系数介于0.73～0.85之间，表明GAMMs模型能够很好地拟合树木年龄和树木胸径大小对IBA的影响。

本研究利用线性混合效应模型来检测树木生长轮剩余生长信号(BAIres)中气候和大气CO2的影响。从表3中可以看出，在所有红松样地中有8个样地表明大气CO2体积分数上升对树木生长呈现出积极影响，其中7个样地表明有显著的积极影响(P<0.05)，剩余1个样地则显示出大气CO2体积分数上升对树木生长不存在显著的影响(P<0.05)。同时，结果还表明大多数研究样地，春季平均气温以及夏季总降水的上升会对树木径向生长产生显著的积极影响(P<0.05)，而夏季高温以及寒冷季节降雪量增加会限制树木的径向生长。

3.2 EC-LUE GPP数据结果与分析

本研究采用修正后的EC-LUE GPP数据来检测β的时空变化趋势，通过对每年生长季的EC-LUE GPP数据进行分析发现(图1)黑龙江省森林1982—2015年间，生长季EC-LUE GPP呈显著的上升趋势(P<0.05)，从1982年892 g·m-2·a-1上升到2015年的978 g·m-2·a-1，年均上升速率为1.04 g·m-2，整体上升了约9.6%。

表2 广义加性混合模型拟合结果

表3 线性混合模型参数系数

1.04表示趋势变化率，*表示P<0.05，灰色阴影区域表示±1标准差，红色虚线为EC-LUE GPP线性拟合。

本研究利用随机森林模型通过控制大气CO2是否变化重建了2种情境下的PGP的年际变化，进而计算二者的差值，并以15 a为窗口逐像元检测1982—2015年黑龙江省森林的β时空变化趋势。在空间上(图2)，黑龙江省森林有88.1%区域存在显著的大气CO2施肥效应，但大多数区域β呈现显著下降趋势(P<0.05)。其中小兴安岭西北部下降最为明显，下降趋势超过3%，而大兴安岭北部和长白山南部的少部分地区出现上升趋势，且大部分地区介于0～3%。

黑点标表示P<0.05。

在1982—2015年，整个黑龙江林区β的均值为27.3%±9.0%，并在整体上呈现出显著的下降趋势(P<0.05)，年均下降速率约为1.1%(图3)。通过对黑龙江省整个林区归一化后的β值分析发现，在1996—2010年以前，β的最大下降速率达到3.5%，在这之后β的下降趋势开始逐渐变缓，并呈现出轻微的上升趋势，年均上升速率约为0.3%。

灰色阴影区域表示±1标准差，红色虚线表示β值的线性拟合。

4 讨论与结论

遥感数据在时间上和空间上均具有较好的连续性，被广泛应用于研究全球大尺度植被与气候的交互作用，为植被与气候关系的研究提供了足够时间长度、一致性和连续性的数据基础[35]。本研究采用一种经过涡度相关数据修正后的光能利用率模型所得到的PGP数据集来检测β值的时空变化趋势，结果表明在1982—2015年间，大气CO2体积分数上升对黑龙江省森林PGP有显著地促进作用，说明大气CO2体积分数上升对森林植被总体的生产力存在施肥效应。然而β值却存在显著的下降趋势(P<0.01)，下降速率约为每年1.1%。1982—2015年期间全球β值均值为16.1%±11.5%，而全球β值以0.92%的速度呈显著下降趋势[28]。一项长期CO2体积分数测量研究发现，PGP对高纬度地区CO2的响应呈下降趋势，陆地植被碳汇对CO2体积分数增加的响应逐渐减少[36]。此外，水蒸气压差的增加也可能导致植物生产力的下降，并抵消部分CO2带来的积极影响，导致了β值的下降[30]。不过β值的下降趋势在1996—2010年开始逐渐变缓，这可能与1998年以后我国开始实行“天保”工程有关。退耕还林和植树造林等政策使得幼、中龄林面积得以增加，有研究表明[37-38]，幼、中龄林树木对于大气CO2体积分数的上升更加敏感，进而缓解了β值的下降趋势。

本研究同时基于红松天然林树木生长轮数据，采用GAMMs模型去除树木自身生理因素后，利用线性混合效应模型分析了大气CO2体积分数对树木生长轮剩余生长信号的影响。通过对树木生长轮进行分析发现，在大多数样地中存在显著的大气CO2施肥效应，这与Chen et al.[39]利用剩余树木生长轮指数观察到的结果一致，Wang et al.[31]同样使用树木生长轮分析发现，中国东北地区兴安落叶松存在显著的大气CO2施肥效应。而本研究样地中的桶子沟原始森林没有检测出大气CO2施肥效应，这可能是由于该地区年均温度较低(年均温度为-2.4 ℃，在所有生长轮采集地区中温度最低)，导致大气CO2施肥效应的信号比较微弱，因此难以通过生长轮数据进行检测。同时也有一些研究表明[40-41]，低温可能会降低羧化酶的活性，进而抵消大气CO2体积分数上升对植被生长的影响。如Camarero et al.[26]利用GAMMs方法分析欧洲山松的施肥效应时，只在潮湿地区的老树中发现了显著的大气CO2施肥效应，这可能是由于研究区域处于高寒山区，降水的匮乏以及温度的限制导致了树木径向生长受到影响。此外，季节性温度和水分之间的相互作用也会影响树木生长对大气CO2体积分数上升的响应[42]。本研究结果也表明，夏季高温会对树木生长有负面影响。夏季高温会加速植被体内水分的蒸腾和土壤水分的蒸发作用，而冬季降雪量的增加则可能会导致积雪不易融化，缩短植被生长季，这一结果与在我国东北森林中的一些研究结果相似[31]。而春季温度上升和充足的降水则有利于植被新芽的生长，进而促进植被的光合作用[43]。

本研究通过自然生长的树木生长轮以及遥感产品清晰地检测到了大气CO2施肥效应的存在，但也存在一些不足之处：(1)收集的树种种类较少，只包含了红松一个针叶树种。不同的树种对大气CO2体积分数提升的敏感度以及生物量体内碳分配的方式有所区别，可能会做出不同的响应；(2)生长轮采样点较少，只包含了9个红松采样地区，不能遍布整个黑龙江森林研究区中；(3)没有加入林分竞争因子以及在样本选择时刻意选取大树和成熟树木进行检测，也可能会对结果产生影响。因此，未来有必要从更多地点收集多个树种的生长轮样芯和采样地区的立地条件信息来检测大气CO2体积分数上升对树木径向生长的影响。

1982—2015年间，黑龙江省森林中存在显著的大气CO2施肥效应，但这一影响在研究区域内的大多数地方呈显著的下降趋势，尤其在小兴安岭西北部；而黑龙江省森林整体β值的变化趋势比较一致，年均以1.1%速率逐渐下降。本研究揭示了植被生产力对大气CO2体积分数上升的响应规律，为评估黑龙江省森林在全球碳循环中扮演的角色提供重要的科学依据。