2013—2020年呼伦湖流域植被碳源/汇估算及影响因子分析

2022-12-06李朝晖李文静王增龙包萨茹窦华山王文林

生态与农村环境学报 2022年11期

李朝晖，单楠①，王琪，李文静，王增龙，包萨茹，窦华山，敖文，庞博，王文林

(1.生态环境部南京环境科学研究所，江苏南京 210042；2.国家环境保护呼伦湖湿地生态环境科学观测研究站，内蒙古呼伦贝尔 021000；3.呼伦贝尔市北方寒冷干旱地区内陆湖泊研究院，内蒙古呼伦贝尔 021000)

人类活动导致大气中CO2、CH4等温室气体浓度不断上升，逐渐影响人类生活，因此全球变暖问题受到越来越多的关注。植被净生态系统生产力(net ecosystem productivity，NEP)指单位时间、单位面积内植物将无机物合成为有机物的总量或固定的总能量与植被呼吸以及土壤异养呼吸(heterotrophic respiration，RH)之间的差值，表示植被生态系统碳的净吸收或净存储，通常用来定量表征植被生态系统碳源/汇能力。当NEP值>0时，表示生态系统为碳汇，反之为碳源[1-3]。碳汇(carbon sink)指森林和草地吸收大气中的CO2，从而降低大气CO2浓度的机制，碳汇在减缓气候变化方面有重要作用[4]。NEP不仅与陆地生态系统多种要素循环紧密相关，还对全球气候系统有重大影响，是表征陆地生态系统碳循环的一个重要指标，对于全球碳循环研究具有重要意义，因此准确监测其时空变化格局成为区域碳收支研究的关键问题之一[4]。

草地生态系统具有防风固沙、水源涵养、水土保持和生物多样性保护等重要生态功能，是中国陆地生态系统重要组成部分。草原碳汇作为生物固碳重要组成部分，对于减缓全球气候变暖具有重要意义。草原碳汇研究是保护草原、恢复退化草地的关键所在。同时，草原面积广阔，通过碳固定措施开展可持续管理，既可以实现草原保护，又可以实现草原碳汇价值增收。目前，已有学者结合遥感及地面观测数据估算了青海高原、青藏高原、三江源和内蒙古等区域植被净初级生产力(net primary productivity, NPP)和NEP，分析了不同区域植被生态系统碳源/汇时空格局变化及影响因素[1-2,5-7]。然而，目前针对呼伦湖流域碳源/汇的研究较少，区域碳汇功能时空格局仍不明晰。呼伦湖流域是北方防沙带和东北森林带交汇的重要节点，同时也是生物多样性重要保护区，东北亚鸟类迁徙的重要通道和集散地。但呼伦湖流域生态环境较为脆弱，生态系统植被变化对区域碳汇功能具有重要影响。基于30 m高分辨率遥感数据和气象再分析数据，采用光能利用率遥感模型和土壤呼吸模型估算呼伦湖流域不同生态系统植被NEP，并分析NEP时空分布格局与变化趋势，评估各生态系统碳库及其动态变化，揭示呼伦湖流域碳源/汇时空格局以及不同生态系统碳源/汇功能对气候变化的响应，研究结果可为呼伦湖流域的生态环境保护和修复、区域生态安全评估等提供支撑[8-9]。

1 研究区与数据

1.1 研究区概况

研究区域为呼伦湖流域，包括呼伦湖以及流域内主要河流(克鲁伦河和乌尔逊河)、草原和林地。流域属于中温带大陆性草原气候区，流域东南部为阿尔山区域。呼伦湖流域海拔在497～1 660 m之间，呈现中心区域海拔较低、呼伦湖及河流周围海拔逐渐升高的地势特征，流域生态系统类型以草原和湿地为主(图1)。

图1 呼伦湖流域地理位置以及DEM和地物覆盖类型

1.2 数据来源与处理

1.2.1Landsat 8遥感影像数据

Landsat 8卫星重访周期为16 d，成像宽幅为185 km×185 km，覆盖呼伦湖流域约需要6景。Google Earth Engine(GEE)是由Google公司、卡内基梅隆大学与美国地质调查局共同开发的用于处理卫星遥感图像和其他地球观测数据的云计算平台。该研究采用GEE平台提供的30 m空间分辨率Landsat 8地表反射率产品LANDSAT_LC08_C01_T1_SR。由于Landsat 8重访周期以及呼伦湖流域云、冰雪的限制，Landsat 8对研究区不能实现有效且完整的覆盖。已有研究[10]表明，对于高纬度且受云层影响较大区域，可以通过采用连续多月反射率中值(median)影像来有效解决。同时，由于7—9月夏季良好的水热条件为植被创造良好生长条件，植被以生殖生长为主，光合作用能力强[11]。故相比于其他月份，7—9月植被固碳能力较强，有助于呼伦湖流域植被高固碳能力期碳汇能力的量化和描述。因此采用2013—2020年呼伦湖流域7—9月反射率中值影像作为各年反射率数据。

1.2.2MODIS数据

MOD17A2是使用1 km范围内每日MODIS土地覆盖数据、fAPAP/LAI和表面气象数据来为陆地植被提供精确且定期测量的生态系统生产力估算产品。该产品基于光能利用率模型计算得到。基于GEE平台获取2013—2020年呼伦湖流域7—9月均值MOD17A2产品，并提取NPP(代码为PsnNet)作为NPP估算模型的验证数据。

1.2.3气象与辐射数据

气象与辐射数据包括月平均气温(temperature, T)、月总降水量(total precipitation, TP)和月总辐射量(surface solar radiation downwards, SSRD)，均来自ECMWF ERA5气候再分析数据集。ERA5数据集是利用物理定律将模型数据与来自世界各地的观测数据整合成一个全球完整且一致的数据集，提供了对过去气候的准确描述。采用ERA5月均值数据(ERA5-Land monthly averaged data from 1981 to present)，提取2013—2020年7—9月的2 m高度处气温、总降水量和地表下行太阳辐射均值分别作为气温、降水量和太阳辐射月均值数据。

1.2.4地表覆盖数据FROM-GLC

FROM-GLC为清华2017年发布的10 m分辨率全球地表覆盖制图数据集[12]。该数据集将全球地表覆盖分为10类，呼伦湖流域共有8种覆盖类型(图1)，其中，植被类型主要包括农田、森林、草地和湿地4类。以上述4类植被类型为基础，将流域非植被类型(包括水体、不透水面和裸土)区域和面积极小的灌木林归为其他。经统计，农田、森林、草地、湿地和其他覆盖类型面积分别为5 513.54、28 951.17、2 780 043.79、3 205.99和468 323.52 km2，其占流域面积比例分别为0.17%、0.88%、84.60%、0.10%和14.25%。地表覆盖数据可用来分析呼伦湖流域2013—2020年不同生态系统类型碳源/汇时空动态。

2 模型与方法

2.1 NEP估算模型

植被NEP为表征区域植被碳汇/源的重要指标，NEP可以表示为植被NPP与土壤异养呼吸消耗碳量的差值，计算公式为

NEP(x,t)=NPP(x,t)-RH(x,t)。

(1)

式(1)中，NEP(x,t)为像元x在时间段t内的植被净生态系统生产力(以C计)，g·m-2；NPP(x,t)为像元x在时间段t内的植被净初级生产力(以C计)，g·m-2；RH(x,t)为像元x在时间段t内的土壤微生物呼吸量(以C计)，g·m-2。

2.1.1光能利用率模型

采用CASA(Carnegie-Ames-Stanford approach)模型估算NPP[13]，其计算公式为

NPP(x,t)=APAR(x,t)×ε(x,t)。

(2)

式(2)中，APAR(x,t)为像元x在时间段t内吸收的光合有效辐射，MJ·m-2；ε(x,t)为像元x在时间段t内的实际光能利用率(以C计)，g·MJ-1。APAR(x,t)计算公式为

APAR(x,t)=SOL(x,t)×FPAR(x,t)×0.5。

(3)

式(3)中，SOL(x,t)为像元x在时间段t内的太阳入射总辐射，MJ·m-2；FPAR(x,t)为植被对入射光合有效辐射的吸收比例，其值采用归一化植被指数(NDVI)进行估算，并假设在时间段t内保持不变；0.5为植被所能利用的太阳有效辐射(波长范围为400～700 nm)占太阳总辐射的比例估算值。

ε(x,t)计算公式为

ε(x,t)=ε0×f1(x,t)×f2(x,t)×w(x,t)。

(4)

式(4)中，f1和f2为温度对光能利用率的2个影响系数，f1(x,t)和f2(x,t)表示温度对光能利用率的影响；w(x,t)为水分条件对光能利用率的影响；ε0为理想状态下植被最大光能利用率，其值参照朱文泉等[14-15]对全国植被NPP计算结果。f1(x,t)、f2(x,t)和w(x,t)计算公式为

f1(x,t)=0.8+0.02×Topt(x,t)-0.000 5×

Topt(x)2，

(5)

f2(x,t)=[1.184 ×(1+e0.3×[-Topt(x,t)-10+T(x,t)])]/

{1+e0.2×[Topt(x,t)-10-T(x,t)]}，

(6)

w(x,t)=0.5+0.5×E/Ep。

(7)

式(5)～(7)中，T(x,t)为平均气温，℃；Topt(x,t)为植被生长最适温度，℃；E和Ep分别为实际蒸散发量(W·m-2)和潜在蒸散发量(W·m-2)。

2.1.2土壤呼吸模型

RH(x,t)计算公式[16]为

lnRH(x,t)=0.22+0.87×lnRs(x,t)。

(8)

式(8)中，Rs(x,t)为土壤呼吸速率，其值采用土壤呼吸模型[17]进行估算。

Rs(x,t)=1.250×e[0.054 52×T(x,t)]×P(x,t)/

[4.259+P(x,t)]。

(9)

式(9)中，T(x,t)为像元x在时间段t内的平均气温，℃；P(x,t)为像元x在时间段t内的降水量，cm。

2.2 时空变化与气候因子分析

年际变化率采用NEP趋势斜率表示，NEP年际变化显著性采用NEP值与时间序列的相关关系表示。趋势斜率大于0表示NEP增加，反之表示NEP减少。2013—2020年NEP波动特征采用NEP变异系数CV表示。CV值越小，表示数据分布越集中，数据随时间的波动小，稳定性较高；反之，则表示数据随时间序列波动大，稳定性低。此外，温度、降水和辐射对NEP的影响采用像元尺度的偏相关分析方法得到。由于该研究有4个变量，因此采用二阶偏相关系数(partial correlation coefficient, PCC)进行分析。

3 结果与分析

3.1 呼伦湖流域典型气候因子与植被指数年际变化

如图2所示，2013—2020年呼伦湖流域气温总体呈现先升高后降低趋势，在2016年达到最高。2013—2020年流域降水量总体呈现先减少后增加趋势，在2016年达到最低。2014—2017年研究区月总辐射量明显高于平均值，表明该时期太阳入射辐射较强，且2014年为最大太阳入射辐射(图2)。2013—2020年流域NDVI值呈现先降低后升高趋势，在2016年达到最低值，这与降水量变化趋势大体一致(图2)。

3.2 呼伦湖流域植被碳源/汇时空变化特点

2013—2020年呼伦湖流域植被碳汇年际变化和时空变化见图3～4。在时间尺度上，碳汇总体呈现先减少再增大的年际变化特征(图3)，最大值和最小值分别为16.44和8.41 g·m-2·月-1(图3)，碳汇年际变化率为-0.12 g·m-2·月-1，植被固碳呈轻微下降趋势。这与研究区降水量变化趋势相似(图2)，与月平均气温变化趋势相反(图2)。在空间尺度上，2013—2020年呼伦湖流域碳汇高值区主要位于呼伦湖北部及西北部、克鲁伦河消落带、乌尔逊河消落带和东南部阿尔山森林区域，低值区主要位于呼伦湖西南部、克鲁伦河两岸、乌尔逊河两岸和阿尔山附近区域(图4)。

图2 2013—2020年7—9月呼伦湖流域平均气温(T)、降水量(TP)、月总辐射量(SSRD)和归一化植被指数(NDVI)的年际变化

R2为决定系数，k为斜率。

图4 呼伦湖流域2013—2020年7—9月和多年平均NEP空间分布

呼伦湖流域东南部为阿尔山森林区域，西北部为呼伦湖及其周边的湿地和草原。研究区内存在气候差异，温度和降水量空间分布不均匀。由图5可知，就NEP经向变化而言，呼伦湖流域碳汇功能呈现明显经度地带性，植被类型由西向东表现为由草地向森林转变，NEP值逐渐上升，碳汇增大；就NEP纬向变化而言，呼伦湖流域NEP随纬度增加呈现先减小后增大特征，其主要原因为东南部低纬度地区森林为碳汇高值区。

阴影表示1倍标准差。

基于像元尺度，采用一元线性回归分析法进一步分析2013—2020年碳汇年际变化情况，以期获得呼伦湖流域不同空间位置上植被碳循环在近8 a变化的区域差异性规律。如图6所示，2013—2020年呼伦湖流域碳汇整体上呈现相对稳定且轻微减少趋势。58.29%的区域植被NEP呈减少趋势，36.77%的区域植被NEP斜率在0～1 g·m-2·月-1之间，呈微弱增加趋势，而斜率高于1 g·m-2·月-1的区域面积仅占4.94%，且主要分布在东南部森林区域。从NEP变异系数计算结果可知，呼伦湖流域植被NEP年际变化变异系数主要集中在0.2～0.5之间，其区域面积占比为82.52%。呼伦湖流域大部分区域碳循环状态年际变化波动较小，总体上保持碳汇功能，但植被固碳能力总体呈下降趋势。

图6 2013—2020年呼伦湖流域NEP变化斜率空间分布和变异系数(CV)空间分布

3.3 不同植被类型的碳源/汇年际时序变化

研究区不同植被类型碳循环状态有较大差异(图7)。从NEP年际变化来看，农田生态系统平均NEP表现为升高—降低—升高—降低趋势，碳汇功能波动变化，NEP均值最低和最高分别为21.29和37.24 g·m-2·月-1；森林生态系统平均NEP值最高，碳汇功能最大，年际变化从2013年的23.88 g·m-2·月-1增大至2015年的最高值44.96 g·m-2·月-1，然后稳定在约40 g·m-2·月-1；草地生态系统碳汇功能显著低于农田和森林生态系统，2014年最高(16.60 g·m-2·月-1)，2016年最低(8.33 g·m-2·月-1)；湿地生态系统碳汇总体高于草地，低于农田和森林，最高值为17.96 g·m-2·月-1(图7)。对于碳汇总量，草地生态系统远高于其他生态系统碳汇总量。结果表明，2013—2020年呼伦湖流域农田和森林均为高固碳功能区，草地为低固碳功能区。流域内农田和森林面积较小，固碳总量较少，草地面积最大，碳汇总量较高。

2013—2020年研究区不同植被类型碳汇多年均值和总量以及频数分布见图7。从频数分布图〔图7(d)〕可以看出，农田生态系统NEP主要分布在20～40 g·m-2·月-1之间；森林生态系统NEP则主要分布在30～50 g·m-2·月-1之间，因此森林碳汇功能(38.10 g·m-2·月-1)最高；草地面积最广，其NEP频数分布的波峰范围为5～20 g·m-2·月-1，因此草地碳汇总量(37.00 Tg·月-1)最高；湿地生态系统NEP则主要分布在小于14 g·m-2·月-1区间，碳汇总量为0.04 Tg·月-1。总体上，农田、森林和草地碳汇功能和总量均为正数，表现为碳汇，碳汇能力由大到小依次为森林、农田和草地。流域内植被碳循环以碳汇为主，仅湿地生态系统少量植被表现为碳源〔图7(d)〕。

(a)不同植被类型的碳汇功能年际时序变化；(b)不同植被类型的碳汇总量年际时序变化；(c)不同植被类型的多年平均NEP与NEP总量；(d)不同植被类型的NEP频数分布。草地/50表示草地NEP总量值除以50。

3.4 呼伦湖流域植被碳源/碳汇对气候因子的响应

植被生长与气候变化密切相关，对2013—2020年呼伦湖流域NEP与温度、降水量和太阳入射辐射的偏相关系数逐像元进行计算(图8)。如图8所示，呼伦湖流域植被NEP与温度主要呈负相关，相关系数在-0.25～-1.0之间，部分区域呈正相关关系，如呼伦湖南部草原(相关系数在0～0.5之间)和东南部森林(相关系数在0.25～1.0之间)。如图8所示，呼伦湖流域植被NEP与降水量主要呈正相关，相关系数在0.25～1.0之间，这些区域分布在呼伦湖、克鲁伦河和乌尔逊河周边草原地区以及流域东南部森林地区。呼伦湖流域植被NEP与辐射的相关性较弱，相关系数主要在-0.5～0.5之间。

图8 2013—2020年呼伦湖流域NEP与温度(T)、降水量(TP)和太阳入射辐射(SSRD)变化的偏相关系数

为了分析不同植被类型NEP对气候因子的响应差异，统计不同植被类型NEP与气候因子的偏相关系数频数(图9)。由于草地像元数远多于其他3种植被类型，其频数分布曲线较平滑。如图9所示，农田NEP受温度影响最大，偏相关系数主要位于-0.25～-0.75和0.25～0.5之间；其受降水量和辐射影响较小，偏相关系数主要分布在-0.4～0.25之间。如图9所示，森林NEP与温度、降水量和辐射主要呈正相关，与各气候因子偏相关系数在0.25～1.0之间的像元数占比分别为74.03%、67.86%和66.23%。如图9所示，草地NEP对气候因子响应表现出明显差异，与温度呈负相关(偏相关系数在-1.0～-0.25之间)的像元数占草地像元总数的75.80%；与降水量呈正相关关系，偏相关系数分布在0.25～1.0之间的像元数占草地像元总数的58.91%；草地NEP与辐射相关性不明显，偏相关系数分布在大于0.5或小于0.5的像元数占比为78.15%。如图9所示，湿地生态系统NEP对温度和降水的响应不明显，其与温度和降水量呈正相关和负相关的像元数占比相差不大，其与辐射主要呈正相关，偏相关系数分布在0～0.5之间的像元数占比相对较高。

图9 不同植被类型NEP对温度、降水量和辐射的响应频数分布

从呼伦湖流域植被NEP对温度、降水量和辐射3种气候因子的偏相关系数分布(图10)来看，呼伦湖流域植被生长主要受温度和降水量影响。呼伦湖四周和乌尔逊河两岸主要为受降水量影响区域，而呼伦湖流域西南部、南部和阿尔山附近区域则主要为受温度影响区域。

图10 影响呼伦湖流域植被生长的气候主导因子分布

3.5 NEP估算模型精度

目前，已有研究中NEP估算模型的验证主要是基于植被地上生物量观测数据或研究区涡度相关通量观测数据[1-2,5,7]。由于缺乏碳通量实测值，该研究基于MODIS NPP产品对模型估算结果进行对比分析。如图11所示，2013—2020年呼伦湖流域估算的NPP与MODIS NPP间呈显著正相关，且相关性较高，皮尔逊相关系数分布在0.47～0.81之间。如图11所示，2013—2020年估算NPP与MODIS NPP之间均方根误差(RMSE)均小于7 g·m-2·月-1。总体来看，该研究得到的估算NPP变化趋势与MODIS NPP具有较好一致性。

图11 2013—2020年呼伦湖流域估算NPP与MODIS NPP皮尔逊相关系数与均方根误差(RMSE)

4 讨论与结论

4.1 讨论

4.1.1气候因子对不同植被类型NEP的影响

热量和水分是决定植被呈地带性分布的控制条件，不同植被类型对热量和水分变化的敏感性存在较大差异[2,18-19]。该研究结果表明森林生态系统NEP受温度、降水量和辐射3种气候因子共同影响，且3者影响程度无明显差异，均对森林植被生长有促进作用。农田生态系统NEP主要受温度变化控制，而草地生态系统受降水量和温度共同影响。该研究结果与以往相关研究结论[19-21]相似，均表明气温和降水量是植被生长主导因子，但其对不同植被类型的影响程度存在差异。然而大部分研究结果表明，存在对特定植被类型的影响程度高于其他气候因子的单个特定气候因子，笔者研究发现2种或3种气候因子呈耦合影响的结论与之存在差异。多种气候因子耦合影响可能是高纬度区域植被对气候因子的响应特征，较低的温度和较少的降水量导致植被未达到最佳光合功能。此外，呼伦湖对流域内植被碳汇也具有不可忽视的作用。呼伦湖流域水热变化受流域中呼伦湖的影响和调控，区域植被生长所需的气温、降水与呼伦湖湖水蒸发关系密切。此外，相比于草地、森林等生态系统，笔者研究发现呼伦湖流域湿地生态系统NEP年际波动较小，表明研究区湿地生态系统对气候因子响应不明显，这可能与湿地生态系统水热补充有关。

4.1.2NEP模拟结果比较与不确定性分析

采用相同模型估算不同区域相同植被类型NEP，其结果存在一定差异。如表1所示，笔者研究估算的2013—2020年呼伦湖流域植被NEP月均值为8.41～16.44 g·m-2·月-1，多年月均值为13.22 g·m-2·月-1。刘凤等[1]估算2000—2015年青藏高原多年NEP月均值为10.70 g·m-2·月-1，戴尔阜等[7]估算2001—2012年内蒙古草原生态系统NEP月均值在4.2～12.5 g·m-2·月-1之间，与笔者研究区域相近的呼伦贝尔平原草原区多年NEP月均值为9.00 g·m-2·月-1。上述研究结果均低于笔者研究，这可能是研究区域植被类型、气候条件和估算时段的差异导致。青藏高原以常绿针叶林、落叶阔叶林、高寒灌丛、高寒草原和草甸以及高寒荒漠等植被类型为主，平均海拔3 000 m以上，冬季漫长寒冷[1]。因此，青藏高原特殊地理位置和气候条件可能导致植被光合条件严峻，植被固碳能力小于呼伦湖流域。而戴尔阜等[7]研究中植被类型为草地，笔者研究的多种生态系统混合(包含森林等高固碳能力植被)与之存在差异，由此导致植被固碳较低。此外，笔者研究估算月份为7—9月，植被生长旺盛，这降低了其他植被生长状态低的月份对NEP月均值估算结果的影响，因此笔者研究结果高于其他研究。潘竟虎等[4]估算中国西北干旱区植被碳汇多年NEP月均值为49.33 g·m-2·月-1，远高于其他研究，这可能是因为其研究区域包括大兴安岭等森林资源丰富区域，碳汇功能强于其他生态系统，因此多年NEP均值高于以草地生态系统为主的其他研究模拟结果。此外，空间分辨率差异可能也是NEP估算结果存在一定差异的原因。

表1 笔者研究NEP模拟结果与其他研究模拟结果的比较

笔者研究模型估算结果有效描述了呼伦湖流域植被碳汇时空格局，但仍存在一定不确定性，主要体现在3个方面：(1)使用的数据来自遥感观测和模型模拟。由于缺少实测数据，在利用NDVI数据估算植被吸收光合有效辐射(FPAR)时可能存在误差。用于计算实际光能利用率的气温、降水量和太阳入射辐射等数据来自ECMWF模型气候再分析结果，因此估算所用的实际光能利用率也存在模拟误差。(2)光能利用率模型假设呼伦湖流域各生态系统植被最大光能利用率均为固定值，且保持一致。然而，已有研究[22-25]发现研究区植被类型、空间位置和水分胁迫等因素均会影响最大光能利用率。因此，笔者研究中最大光能利用率的估算是实际光能利用率的误差来源之一。(3)土壤呼吸受到多种因素影响，而笔者研究土壤呼吸模型对不同生态系统土壤类型、植被类型和土壤微生物均采用相同的模型参数，这可能影响土壤呼吸估算的准确性。