退耕还草对锡林郭勒草原总初级生产力变化的影响

连浩宇,李愈哲

1 中国科学院地理科学与资源研究所陆地表层格局与模拟院重点实验室,北京 100101

2 中国人民大学公共管理学院,北京 100872

固碳是陆地生态系统一项重要的服务功能,地表植被通过光合作用吸收空气中的CO2并固定有机碳,可以减少温室效应、调节生态平衡、联系大气和陆地的物质能量交换[1]。定量评估陆地生态系统的固碳潜力是区域性碳汇管理和气候调节的基础[2],也是碳达峰碳中和战略背景下的热点问题[3]。其中,陆地植物总初级生产力(GPP)是全球最大的碳通量,描述了单位时间内绿色植物通过光合作用所固定的有机碳量,是定量化生态系统固碳潜力的关键参数[4]。GPP决定了进入陆地生态系统的初始物质和能量,可以反映植被整体的生长状况和固碳潜力[5];它受气候变化影响的同时,对气候系统具有反馈作用,可用于分析固碳过程带来的气候效应[6]。

随着近年遥感技术的发展,利用模型估算区域尺度GPP成为了一种重要且普遍的研究方法[7]。GPP估算模型包括植被光合作用模型(VPM)和植被光合作用和呼吸模型(VPRM)等光能利用率模型[8],回归树模型等统计模型[9]以及生物群落-生物地球化学循环模型(Biome-BGC Model)等过程模型[10]。其中,VI模型由Wu等在2010年建立,该模型基于植被指数和光合有效辐射(PAR)进行估算GPP,模型中的光合有效辐射分量(FPAR)与植被光能利用率(LUE)均通过植被指数进行量化,并将此方法应用到玉米地的GPP估算[11]。VI模型估算的GPP已经过大量的地面验证和地面校准工作,具有较高的数据精度,适宜做大尺度区域 GPP的时空变化研究。王克清等[12]结合中国通量观测研究联盟(ChinaFLUX)的台站数据,对中国典型植被类型的GPP进行了模拟,其中VI模型在锡林郭勒温性高原的模拟值与实测值相关系数R2达到0.86,模拟效果良好。

内蒙古草原属北方温性草原,是我国重要的防沙治沙屏障,具有涵养水源、净化环境作用。同时,内蒙古草原干燥多风沙、地表物质结构脆弱、土壤侵蚀强烈,生态环境极其脆弱[13—15]。本世纪初,全盟风蚀沙化面积12.6 km2,占全盟草原总面积的64%,其中强度风蚀沙化面积1.48 km2,植被覆盖率下降至27.2%;同时水土流失呈加剧趋势,全盟轻度以上的水土流失面积17.2 km2,占总土地面积的81.2%;其中浑善达克沙地从1949 年到1995 年沙漠化面积由0.57万 km2增加到3.05万 km2,平均每年以100 km2多的速度增加。

为了改善内蒙古草原的生态环境,我国政府开展了一系列生态修复工程,在锡林郭勒盟实施的退耕还林(还草)工程,在农牧交错带主要进行退耕还草[16]。随着“退耕还草工程”的开展,锡林郭勒草原的耕地面积明显减少[17],其中63.95%的耕地转化为草地[18]。截至2015年,全盟植被覆盖率上升至49.66%,产草量达到3.89 kg/hm2,重度沙漠化与中度沙漠化面积均有下降。

许多学者研究发现退耕还草可以显著提高草地植被数量,增加土壤养分和有机碳含量[19],增加土壤水分和碳固存[20],在退耕还草中期有显著的碳氮增汇效应[21]。但是目前关于GPP变化量与气候、土壤、地形等环境因子间响应机制的研究并不多见,对退耕还草后GPP变化的影响因素也缺乏控制变量的研究。基于此,研究利用VI模型对2010—2015年锡林郭勒草原退耕还草区域GPP变化量进行估算,并与同时期锡林郭勒草原各环境因子进行相关性分析,定量评估环境因子对退耕还草前后GPP变化量的影响。本研究有利于掌握退耕还草前后GPP的变化,明确退耕还草后GPP变化与环境因子的关系,为内蒙古草原退耕还草政策的制定提供科学基础和决策依据,为我国草原生态治理和区域碳中和提供建议。

1 材料与方法

1.1 研究区概况

研究区位于中国内蒙古自治区锡林郭勒盟(42°32′—46°41′N,111°59′—120°00′E)(图1),地势南高北低,南部为低山丘陵,北部为平缓的波状平原,总面积大致为14785 km2。研究区地处中温带半湿润到半干旱气候区,属中温带半干旱大陆性气候,年平均气温为0.6℃,年平均降水量294.9 mm。研究区内可利用草场面积138.56万 hm2、耕地17085 hm2,草原类型较为丰富,地跨草甸草原、典型草原和沙丘沙地草原,栗钙土为主要土壤类型[22]。本研究设立的两个涡度相关观测系统分别位于典型温性草原的草地生态系统(43°32′24″N,116°33′43″E)和农田生态系统(43°35′20″N,116°45′43″E)的中心地带,设立地块地形相对较平缓,草地植被建群种为大针茅(Stipagrandis);农田生态系统种植作物为小麦(Triticumaestivuml),种植区域不进行灌溉,在作物种植前施有机肥(表1)。

表1 草地和农田生态系统植被和土壤特征差异

图1 研究区概况图

锡林郭勒盟高程在760—1911 m,坡度主要在0°—2°,平均海拔在1000 m以上,最高峰是西乌珠穆沁旗内的古如格苏乌拉峰;地势南高北低,自西南向东北倾斜,东南处地势最高,东北处地势最低(图2)。西部和北部地形平坦,东南部多低山丘陵,盆地错落其间,形成广阔的高原草场;东北部和东部地区以低山丘陵、高平原和宽谷平原为主,是森林向草原的过渡地段。

图2 高程和坡度的空间分布

锡林郭勒盟属中温带干旱、半干旱大陆性季风气候,年平均降水量在99.58—531.29 mm,其中东部降水较多,大兴安岭余脉西坡及阴山余脉北坡局部地区在400 mm以上,西部降水较少,局部不足150 mm;锡盟年平均气温在-4.76—4.90℃,西南部气温较高,东北部气温较低,北部中蒙边境地区和灰腾梁一带年平均气温0℃以下(图3)。

图3 多年平均降水和气温的空间分布

锡林郭勒盟的土壤水分在0.07—0.20 m3/m3,其中东部土壤含水量较高,西部土壤含水量较低,与降水量由西向东递增有关;锡盟有机碳含量在0.20—5.98 g/kg,其中西南部有机碳含量较低,主要是荒漠草原,中部有机碳含量在10.00—20.00 g/kg,主要是典型草原,东部有机碳含量较高,平均在20.00 g/kg以上,主要为草甸草原,土壤有机碳含量与草原类型呈一致分布(图4)。

图4 土壤含水量与土壤有机碳的空间分布

1.2 数据来源与处理

1.2.1遥感数据

本研究的遥感数据来源于美国航空航天局(NASA)(https://ladsweb.modaps.eosdis.nasa.gov)提供的Modis/Terra卫星遥感产品MOD09A1,空间分辨率为500 m,时间分辨率为8 d,数据选取编号为h26v04的影像,产品数据下载后需要对遥感数据进行格式转换与投影变换。同时,本研究选择时域线性插值的方法对数据中的异常值进行修复,采用Savitzky-Golay滤波法对遥感数据进行平滑处理,表达式为:

(1)

1.2.2环境因子数据

本研究所采用的高程,气温,降水数据,均来自中国科学院地理科学与资源研究所的资源环境科学与数据中心(https://www.resdc.cn/Default.aspx)。其中,海拔高度(DEM)空间分布数据是基于的航天飞机雷达地形高程数据(SRTM V4.1)经整理拼接生成的90 m的栅格数据,采用WGS1984椭球投影。土地利用数据(LUCC)是在2010年土地利用遥感监测数据的基础上,基于Landsat遥感影像,通过人工目视解译生成,空间分辨率为1000 m。

本研究使用的坡度、土壤水分、土壤有机碳、土壤全氮、全磷、全钾含量数据均来自国家地球系统科学数据中心(http://www.geodata.cn/),其中坡度数据的空间分辨率为1000 m,将DEM原始数据通过投影转换、图幅拼接、数据重采样、数据裁剪,最后将投影拼接后的SRTM数据进行坡度计算得到。土壤数据的时间尺度为2010—2015年,空间分辨率为1000 m,通过Arcgis 10.7软件裁剪至研究区范围,利用栅格投影统一坐标系,最终得到锡林郭勒盟土壤数据集。

1.2.3地面观测数据

本研究通过在研究区内站点设置通量塔获得实际GPP数据,通量观测塔分别位于草地生态系统和农田生态系统的中心区域。利用中国通量观测研究联盟(ChinaFLUX)提供的程序对涡度相关系统观测的原始采样数据进行处理,包括剔除野点、密度效应校正(WPL)、坐标轴3次旋转变化修正等,得到采样间隔为30 min的GPP数据。研究区内草地和农田站点的能量闭合度分别为82.3%和81.6%,GPP观测数据质量较好,处于国际中上水平[23]。

1.2.4土地利用类型数据

本研究采用的土地利用类型数据来自中国科学院地理科学与资源研究所资源环境科学与数据中心(https://www.resdc.cn/Default.aspx),研究数据为2010年及2015年地利用现状遥感监测栅格数据,空间分辨率为1000 m。土地类型包括了耕地、林地、草地、水域、居民地和未利用土地6个一级类型以及25个二级类型。

1.3 数据分析与检验方法

1.3.1GPP计算

本文选用VI模型估算GPP。由于VI模型仅有植被指数(EVI)和光合有效辐射(PAR)两个模型驱动参数,结构较为简单。计算公式为:

GPP=(EVI×EVI×PAR)×m

(2)

式中,的参数m为光合转换系数,随植被类型和环境而变化,同时通过模型标定得到,是模拟GPP与站点实测GPP线性拟合曲线的斜率。

本研究所用PAR数据基于研究区及其周围气象站点的观测数据(观测数据来自国家气象局网站http://data.cma.cn/),利用ANUSPLIN插值软件,结合数字高程模型(数据来源于美国奋进号航天飞机的雷达地形测绘SRTM),经过空间插值获取得到。

本研究所用植被指数EVI由MODIS遥感数据的地表反射率产品计算得到的。其中,MODIS格网范围在H26V04,使用的遥感产品为地表反射率产品MOD09A1;EVI考虑了蓝光波段,增强了对植被覆盖度变化的敏感性,具体计算公式如下:

(3)

式中,ρNIR代表近红外(841—876 nm)波段,ρRED代表红光(620—570 nm)波段,ρBLUE为蓝光(459—479 nm)波段。

1.3.2残差分析法

本研究采用Evans和Geerken提出的残差分析法[24],逐像元量化退耕还草对GPP的影响。首先对锡林郭勒盟2010—2015年退耕还草的区域,基于1985—2010年的气温,降水和GPP数据,逐像元建立GPP与年平均降水量,年平均气温的回归模型。其中,1985—2010年的GPP数据由VI模型计算得到。随后将2015年平均降水量和平均气温代入模型,获得退耕还草区域2015年GPP预测值(GPPpre),所得GPP预测值只受气候变化影响。随后通过遥感观测退耕还草区域2015年PAR和EVI数据,利用VI模型,得到2015年的GPP实际值(GPPreal),GPP实际值受到气候变化与退耕还草的共同影响,整体上高于GPP预测值,具体结果如下(图5):

图5 GPP实际值与预测值

最后得到GPP实际值与GPP预测值之间的残差序列即为GPP变化量(ε)。GPP变化量可以反映不受气候影响,仅由退耕还草造成的GPP变化。残差分析表达式为:

ε=GPPreal-GPPpre

(4)

GPPpre=aT+bP+c

(5)

式中,GPPreal为2015年GPP的实际值,GPPpre为2015年GPP的预测值;a、b分别为GPP对气温和降水量的回归系数;c为回归常数项;T表示气温,℃;P表示降水,mm;ε>0表示退耕还草对GPP有增加作用,ε<0表示退耕还草对GPP有减少作用,ε=0表示退耕还草不改变GPP。

1.3.3相关性分析

本研究采用SPSS 26.0对GPP变化量与各环境因子进行Pearson相关分析,偏相关分析和复相关分析。涉及样本共计379个,均为研究区内管理方式发生转变的退耕还草涉及地块,每个像元点对应一组环境因子数据及GPP变化量。

2 结果

2.1 VI模型精度评价

本文通过比较锡林浩特市内草地站点和农田站点的实测GPP与理论GPP进行模型精度分析。总体上看,VI模型模拟的GPP值与涡度相关系统的测量值有较强的一致性,说明在8 d尺度上VI模型可以利用EVI与PAR较好地模拟出不同站点GPP的变化。

VI模型的模拟精度受不同土地利用方式的影响,2010年与2015年在草地区模型模拟的GPP值与实测GPP值在线性回归方程中R2为0.79,均方根误差(RMSE)为1.46 g C/m2,模拟值与观测值之间离散程度较低(图6)。农田生态系统区GPP模拟值与实测值的线性回归方程中R2为0.68,RMSE为2.78 g C/m2。检验结果说明VI模型对GPP的拟合情况较好,模拟GPP与实际GPP的误差较小,同时草地区的拟合精度高于农田区。

图6 VI模型站点模拟结果

草地站点模型模拟值与实测值之间线性拟合直线的斜率为0.79,很大部分模拟值与实测值的分布低于1∶1直线;农田站点GPP模拟结果与实测结果的线性拟合直线的斜率为0.68,模拟值与实测值的拟合线略低于1∶1直线,是VI模型未考虑LAI变化对植被的影响,造成了GPP值的低估[25]。

2.2 退耕还草区域GPP变化量

2010—2015年,锡林郭勒盟退耕还草区域的GPP总体稳定,退耕区域内GPP变化量平均值为0.47 g C m-2d-1。其中GPP增加的区域占研究区面积的67.20%,GPP平均增加1.98 g C m-2d-1,最多增加10.16;GPP降低的区域占研究区面积的32.80%,GPP平均减少3.39 g C m-2d-1,最多减少16.71 g C m-2d-1(表2)。

表2 GPP变化情况

2.3 退耕还草区域的空间布局

本研究通过比较锡林郭勒盟2010年及2015年土地利用类型数据,计算5年内各行政区耕地变为草地的面积。2010年—2015年锡林郭勒草原退耕还草区域共128.23 km2(表3),退耕还草区域主要集中在锡林郭勒盟南部,东北部也有零星分布,以东乌珠穆沁旗、西乌珠穆沁旗、锡林浩特市、正蓝旗、多伦县、正镶白旗和太仆寺旗为主。

表3 退耕还草面积

2.4 退耕还草GPP变化量的影响因素

为了探究决定退耕还草后GPP变化量的主要环境因子,用回归分析来确定退耕还草区域内GPP变化量与环境因子的关系。结果表明,气温、降水、高程、坡度均与GPP变化量呈线性负相关(P<0.01,表4);土壤含水量、土壤有机碳均与GPP变化量呈线性正相关(P<0.01,表4);土壤氮含量、土壤磷含量、土壤钾含量则与GPP变化量之间的关系不显著(P>0.05,表4)。

表4 GPP变化量与环境因子相关性

通过偏相关分析进行环境因子的进一步筛选发现,坡度和土壤含水量是影响 GPP变化量最主要的因素(P<0.01);高程和气温与GPP变化量也有显著偏相关,但偏相关性较弱(P<0.01,表5);降水量和土壤有机碳与GPP变化量无显著偏相关(P>0.05,表5)。以上述六种环境因子与GPP变化量建立线性回归模型,模型的复相关系数为0.652,达到极显著水平(P<0.01,表5),说明环境条件对退耕还草后GPP的变化量有显著影响。

表5 GPP变化量与环境因子偏相关性

3 讨论

3.1 退耕还草对GPP的影响

2010—2015年,锡林郭勒盟退耕还草区域的GPP总体稳定,GPP变化不显著。Jiang 等研究了近30 年来黄土高原典型黄土区植被生产力时空变异,发现其整体呈现上升趋势,但以2006为断点,在此前30%—40%的植被区域表现出植被生产力损失,随后逐渐恢复[26]。此外,Feng 等针对黄土高原退耕还草工程进行研究,发现植被生产力稳步上升,并在2016 年上升趋势趋强[27]。本研究表明5年内锡林郭勒草原退耕还草区域GPP变化不显著,部分地区GPP有所下降,退耕还草区域的GPP平均值仅略微上升。这一与其它区域有所不同的现象归咎于区域涉及退耕还草的地块相对更少,同时在温性草原退耕还草对于GPP的影响弱于其它区域和植被类型,即研究区开垦地块表现出的GPP未大幅低于其自然状态植被。这一现象是否具有更大时空范围的普遍性,后续研究可以通过分析更长时段的退耕还草措施持续时间、扩大研究区范围进行验证探讨。

此外,有研究表明退耕还林、草地治理等其它工程措施也可以不同程度提高GPP。任小玢[28]对宁夏草地变化的研究表明,耕地扩张等人为活动是导致草地退化的绝对主导因素,而草地治理,退耕还草等生态保护政策可以有效提高温性草原的植被生产力。茆杨等[29]对西南地区植被初级生产力时空变化的研究发现,退耕还林可以显著提高植被初级生产力,同时植被生产力的提高与累计造林面积存在显著关联。系列结果不仅暗示不同工程措施的GPP影响存在极大的区域性差异,还表明研究区GPP变化与工程措施所涉面积紧密相关,在退耕还草实施力度和规模较大的区域仍对区域GPP产生更大程度的改变。

3.2 环境因子对GPP变化量的影响

本研究发现土壤含水量和土壤有机碳均与GPP变化显著负相关(表4,P<0.05),气温、降水、高程和坡度均与GPP变化量显著正相关(表4,P<0.05)。研究发现水热条件较好的耕地,退耕后GPP呈上升趋势,与Liu等[30]针对黄土高原退耕还草工程的研究结论一致,即植被初级生产力并非在所有区域都呈现出生态恢复的状态,而只是在水热条件较好的区域增加。刘洋洋也发现退耕后植被生产力在不同类型以及不同地形下的草地变化趋势存在差异,降雨、人为活动是驱动其变化的主要因素[31];退耕后植被初级生产力变化在不同空间上的恢复状态、速率不同,在不同时间段的恢复状态、速率也不同[30]。这表明在未来退耕还草实施过程中,综合考虑实施地块的环境因子将有助于遴选出更具恢复价值和潜力的地块,在水热条件较好区域实施退耕还草能够更大程度增加区域的GPP。此外,不同植被类型下,GPP 对相同气候因子的响应程度存在差异,王大为等[32]研究表明,落叶阔叶林和常绿针叶林 GPP 主要受蒸发和相对湿度影响,草地主要受降水量影响,耕地主要受平均气温影响,灌木林、湿地和荒漠植被主要受水汽压影响,且越靠近干旱地区,水汽压对 GPP 的影响越强。

同时研究还发现,地势较高,坡度陡峭的区域退耕后GPP呈现上升趋势。这一现象与彭文英等发现渭河流域在25°的坡度上退耕还林效果最大化结论相似[33],说明在退耕实施中应优先退耕坡度较高的耕地。最后,土壤含水量较高的耕地退耕后GPP呈下降趋势。这是灌溉农田表现土壤含水量较高,显著高于自然退耕还草地和周边天然草地[34];相应地块退耕后土壤水分损失较大,从而造成退耕后草地GPP形成一定水分限制而下降[35]。由此可见,不恰当的实施退耕还草在部分区域和地块可能对草原的初级生产力和固碳潜力产生负面影响[36]。为更好实现退耕还草最优的生态效益,应合理选择退耕区域,减少退耕初期植被生产力损失,提高草原固碳潜力。

3.3 退耕还草对生态系统固碳潜力的影响

草地生态系统固碳主要通过光合作用将大气中的二氧化碳转化为有机碳,固定在植物体内或土壤中;而GPP反映了植物在生态系统水平上通过光合作用初始固定的有机碳总量,代表植物的光合生产能力,可以较好地量化生态系统的固碳潜力。锡林郭勒草原在退耕还草5年内,大部分地区GPP呈上升趋势,但GPP总量变化不显著(表2);毛绍娟等对日喀则河谷退耕还草工程的研究表明,退耕初期由于植物种单一,植被地上碳密度较低,随着退耕时期延长,植被碳密度有所增加,生态系统固碳潜力提高[37]。其中,2010—2015年退耕还草造成的GPP变化量有限,固碳潜力提高不明显,与没有建立有效的侵蚀控制机制,以及与退耕初期施肥的干扰作用有关[38]。此外,水热条件较好的耕地,退耕后固碳潜力提高显著(表4),与马南方等的研究结果类似,其结果显示退耕后固碳效益与气候密切相关,越湿润的地方退耕后固碳收益越显著[39];邹婧汝的研究也发现适宜的气温和降水有利于退耕后次生草地的有机物积累和植物萌发,加快植被对土壤养分的吸收,提高草地群落的多样性,对植被生产力和生态系统固碳有正向作用[40]。结合本研究的结果(表4),退耕区水热条件是影响草地退耕后固碳收益的最重要影响因素。在退耕还草措施以外,区域的整体GPP在研究期间还受到草原生态修复和补助奖励等其他非工程措施性政策因素的影响,其贡献和作用机还需更深入的研究加以剥离量化。

退耕还草通过播种优良牧草,提高植被高度、盖度和生物量,改善植物群落结构[41],增加了植物群落的光合速率和叶面积指数[42],有利于生态系统的碳固存[43];同时农地转变为多年生植被后,可以有效防控土壤侵蚀,改善土壤的物理、化学以及生物学质量,提高植被固碳效益[44]。研究表明,在<450mm降雨带(如锡林郭勒草原),土壤水分不足以支撑乔灌的正常生长,区域生态系统碳的增加多集中在草地生物量中[45];经过长期封育退耕,植被恢复后,枯落物积累量增多,加之土壤微生物和土壤酶的作用,退耕还草区域的土壤有机碳密度能显著高于周边原生草地[46]。关于退耕还草后固碳潜力随时间的变化,史利江等的研究表明从短期(<20 a)来看,退耕还林(草)的土壤固碳效应不明显,20 a后随着恢复年限的增加,植被的土壤碳氮增汇效应显著[21];许明祥等则认为以10 a为界,退耕还草的短期土壤碳增汇效应不明显,而10 a后土壤碳增汇效应逐渐明显[47]。相关研究均认为短期内退耕还草的固碳效益并未显著高于农田,与本研究结果相互印证,但为了进一步摸清区域长期退耕的效果,还需要未来对退耕后更长时段进行更深入分析与研究。

4 结论

本文以MODIS GPP数据为基础,系统研究了锡林郭勒草原2010—2015年退耕还草区域GPP的变化特征,并探讨了退耕还草后GPP变化与主要气候因子的关系,获得结论如下:

(1)2010—2015年,锡林郭勒盟退耕还草区域的GPP总体稳定,GPP变化量平均值为0.47 g C/m2,没有显著上升趋势。其中GPP增加的区域占研究区面积的67.20%,GPP降低的区域占研究区面积的32.80%;GPP变化量最大值为10.16 g C/m2,最小值为-16.71 g C/m2。

(2)锡林郭勒草原退耕还草5年内,退耕区域的土壤含水量和土壤有机碳与GPP变化量显著负相关,气温、降水、高程和坡度均与GPP变化量显著正相关,其中土壤含水量和坡度对GPP变化量的影响最大;可见退耕时应优先选择水热条件好,坡度较高低的耕地,减少退耕初期GPP的损失,增加固碳效益。