全球气候变化与陆地生态系统响应是国际地圈生物圈计划的核心研究重要内容, 得到了国际科学界和国际社会的高度关注^[1]。植被作为陆地生态系统的主体, 是连接土壤、大气和水分等要素的自然纽带, 影响着地球系统的能量平衡, 在气候、水文、土壤等方面扮演关键角色, 陆地生态系统的任何变化必然在植被类型、数量或质量方面有所响应^[2]。在全球气候变化的大环境下, 植被变化备受关注, 在地球表面的能量交换、碳循环和区域人类活动中起到重要作用^[3], 会影响全球碳循环和气候系统的稳定性, 对研究全球气候变化背景下的陆地生态系统具有重要的意义, 地表植被覆盖的研究也是全球变化科学研究的热点领域^[4]。同时, 植被覆盖变化也是生态环境变化的直接结果, 它在很大程度上代表了生态系统功能及生态环境质量状况, 也是评价生态环境质量的重要参数, 对生态环境时空变化具有一定的指示作用^[5]。

21世纪以来, 随着全球气候变化和人类活动加剧, 中国土地退化、水土流失、沙漠化等一系列生态环境问题频发, 自然生态系统服务能力严重下降, 生态环境问题已成为人们最为关注和迫切需要解决的热点问题^[6—7]。为此, 生态保护恢复成为保证我国经济可持续发展的必然选择, 已经上升为国家战略, 成为实现生态文明建设的重要途径^[8]。我国开展了一系列生态保护修复工程, 这些工程的实施对生态保护修复起到了一定的促进作用, 生态系统恶化的趋势得到遏制^[9]。2020年, 国家发展改革委、自然资源部联合印发了《全国重要生态系统保护和修复重大工程总体规划(2021—2035年)》^[10], 对增强中国生态服务功能、保障国家生态安全具有重要的意义。退化生态系统的修复通常离不开植被的再恢复过程, 植被的修复深刻记录了人类活动的烙印^[11]。对大尺度、长时期的植被变化进行监测已成为生态保护修复区域生态监测和成效评估的重要内容, 并可为人类制定合理的土地利用及生态环境保护策略等提供科学依据^[12—13]。《全国重要生态系统保护和修复重大工程总体规划(2021—2035年)》中涉及的重要生态保护修复工程所在区域的自然植被覆盖度本底和时空变化特征如何, 目前尚不清晰, 因此亟需对这些区域在大规模生态恢复措施背景下的植被变化进行定量监测与评估, 有助于全面准确掌握我国生态系统保护和修复“家底”, 为中国生态环境保护和建设提供决策支持。

如何量化全国尺度植被变化、评估生态修复成效, 成为当前陆地表层生态系统研究领域共同的科学问题, 也是当前生态学研究的热点之一^[14]。绿色植被的物理和功能特性, 如植被结构、覆盖率和生产力等植被参数常作为评估植被恢复成效的重要指标^[15]。其中, 植被覆盖度(FVC)是反映地表植被群落生长态势的重要指标和描述生态系统的重要基础数据, 它与生长状况、生物量和光合作用强度密切相关, 能够在较大时空尺度上客观反映研究区植被覆盖信息, 被认为是检测区域生态修复有效指标, 对区域生态系统环境变化有着重要指示作用^[16], 同时其能够有效降低无植被区域光谱特征带来的不确定性, 提升分析精度^[9]。国内外学者关于生态保护修复区域植被覆盖度的研究无论在全球、全国还是区域尺度上都取得了显著进展, 许多学者已经在各尺度区域、多种时间范围进行探索^[17—26], 这些研究对于认识中国各区域植被变化的特征及驱动因素等具有重要影响。

卫星图像是监测大规模植被动态变化的重要数据来源之一, 也是生态保护修复成效评估中广泛使用的工具^[27]。针对大范围、长时间序列植被覆盖变化, 通常使用遥感手段来监测。因此本文基于2000—2020年的MODIS-NDVI遥感数据获取了植被覆盖度(FVC)的时空变化信息, 采用增幅变化分析、Theil-Sen Median趋势分析、稳定性分析方法, 对《全国重要生态系统保护和修复重大工程总体规划(2021—2035年)》涉及的青藏高原生态屏障区、黄河重点生态区、长江重点生态区、东北森林带、北方防沙带、南方丘陵山地带六大区域的自然植被覆盖度本底和时空变化特征进行了研究, 拟研究的具体问题如下：① 2000—2020年六大区域植被覆盖特征及不同植被类型覆盖的时空格局如何变化。②气候因子主要是降水、温度与植被覆盖度的时空变化相关性如何。研究以期为未来中国的生态保护和修复工程实施选区提供基础本底数据参考, 对分析中国重要生态保护修复区域对气候变化的响应以及生态与环境变化的趋势具有重要参考意义。

1 研究区概况与数据来源 1.1 研究区概况

本文研究区域选取了青藏高原生态屏障区、黄河重点生态区、长江重点生态区、东北森林带、北方防沙带、南方丘陵山地带六大重点陆域区域作为研究范围(图 1), 总面积约为7900×10³km², 涉及全国27省(区、市)。其中, 青藏高原生态屏障区以高原地貌为主, 海拔多在3000m以上, 属青藏高寒气候, 年均降水量大多400mm以下, 自东南向西北水热条件呈现由暖湿向寒旱过渡的特征；黄河重点生态区大部分位于干旱、半干旱地带, 地貌以山地、丘陵、高塬为主, 下游地区以平原为主, 本区域属大陆性季风气候, 丰水年和干旱年降水量相差2—5倍, 年降水量在150—750mm, 时间和空间分布十分不均；长江重点生态区大部分属典型的亚热带季风湿润气候, 年平均降水量500—1400mm, 具有复杂的地质构造和多样的地貌类型, 包括高原、山地、盆地、丘陵、平原等类型；东北森林带地貌类型多样, 分布着大兴安岭、小兴安岭、长白山地、松嫩平原和三江平原, 温带季风气候显著, 自南向北地跨中温带和寒温带, 降水量400—1000mm, 是我国沼泽湿地最丰富、最集中的区域；北方防沙带属干旱、半干旱地区, 沙化土地广布, 温带大陆性气候显著, 大部分地区年降水量在400mm以下；植被稀疏, 以草原、灌木、荒漠为主, 土地沙化、次生盐渍化严重, 是我国生态环境最脆弱的地区之一；南方丘陵山地带以丘陵地貌为主, 低山、盆地, 南岭山地横贯东西, 属于热带、亚热带季风气候, 雨热同季, 年平均降水量1000—2500mm^[10]。

1.2 数据来源 1.2.1 遥感和气候数据

遥感数据是基于美国国家航天局MOD13A1产品数据集(http://ladsweb.modaps.eosdis.nasa.gov)采用年平均值计算得到的年度NDVI时序数据, 空间分辨率为500m, 时间分辨率为16d, 时间尺度为2000年至2020年。由于MODIS/NDVI产品经过水、云等处理, 保证了质量, 同时拥有较高的空间分辨率, 被广泛应用于区域植被覆盖变化研究。基于Google Earth Engine(GEE)平台, 利用像元二分模型计算得到年平均植被覆盖度数据(FVC), 通过工程区域的矢量边界裁剪得到2000—2020年六大区域FVC数据, 空间坐标系为Albers_Conic_Equal_Area。

气候数据主要为2000—2020年空间分辨率为1km的温度和降水量TerraClimate数据集(来源于 https://climate.northwestknowledge.net/TERRACLIMATE/index), 为保证数据的一致性, 同样将温度和降水量数据采用平均值合成法生成年值数据。

1.2.2 矢量数据

矢量边界数据来源于2020年国家发展改革委和自然资源部发布的《全国重要生态系统保护和修复重大工程总体规划(2021—2035年)》中七大工程区域规划示意图, 主要选取了青藏高原生态屏障区、黄河重点生态区、长江重点生态区、东北森林带、北方防沙带、南方丘陵山地带六大区域示意图进行下载^[10], 利用ArcGIS软件并结合全国行政边界底图数据进行空间位置校正, 然后根据中国行政县域进行区域范围选取得到六大区域矢量数据, 空间坐标系为Albers_Conic_Equal_Area。

2 研究方法 2.1 植被覆盖度估算

利用遥感技术手段获得植被覆盖度的方法主要有经验模型法、植被指数法和像元分解模型法。像元分解模型中的像元二分模型具有计算模型简单可靠、数据参数通用易得、反演精度较高的优点, 因此采用像元二分模型来反演研究区域的植被覆盖度。根据研究需要, 本文基于GEE平台, 利用归一化植被指数通过像元二分模型方法^[27]估算植被覆盖度, 计算公式如下：

式中, FVC代表像元的植被覆盖度; NDVI为该像元上的归一化植被指数值; NDVI_veg为影像上植被部分对应的NDVI值; NDVI_soil为影像上裸地部分对应的NDVI值。参考植被覆盖度的遥感估算方法^[28], 并基于研究区域内的土地利用数据, 根据ENVI计算像元NDVI的累计分布情况, 截取置信度99.5%和0.5%分别作为NDVI最大值和最小值, 即NDVI_veg和NDVI_soil。利用年度合成的NDVI数据, 通过以上公式计算得到2000—2020年研究区域的植被覆盖度数据。参考中华人民共和国水利部2008年颁布的《土壤侵蚀分类分级标准》, 并结合相关研究^[29—30], 将植被覆盖度(FVC)划分为5个等级: ＜30%(低植被覆盖度)、30%—45%(中低植被覆盖度)、45%—60%(中植被覆盖度)、60%—75%(中高植被覆盖度)、＞75%(高植被覆盖度)。

2.2 变化趋势分析

为研究区域植被生长状况的变化及变化速率在空间上的差异, 本文结合R语言, 利用Theil-Sen Median方法^[31]来计算研究区域2000—2020年植被覆盖度的变化趋势。Theil-Sen Median方法又称为Sen斜率估计, 是一种稳健的非参数统计的趋势计算方法。该方法计算效率高, 对于测量误差和利群数据不敏感, 适用于长时间序列数据的趋势分析。其计算公式为：

式中, Median()代表中值, β值表示其变化增减趋势, β>0表示其变化为增加趋势, β＜0表示其变化为减少趋势。通过Mann\\Kendall统计检验法对变化趋势的显著性进行检验^[32], 以P值表示, 显著性仅代表趋势性变化可置信程度的高低。P＜0.01为极显著变化, 0.01＜P＜0.05为显著变化, P＞0.05为不显著变化。结合以上趋势增减和显著性变化, 本文将植被覆盖度变化程度划分为五个等级：极显著增加(β>0, P＜0.01)、显著增加(β>0, 0.01＜P＜0.05)、不显著变化(P＞0.05)、显著减少(β＜0, 0.01＜P＜0.05)、极显著减少(β＜0, P＜0.01)。

2.3 稳定性分析

引入变异系数来研究长时间序列植被覆盖度的稳定性^[33]。变异系数的大小反映了数据系列的波动程度, 变异系数越小, 植被稳定性越好；反之, 变异系数越大, 植被稳定性越差。其计算公式为：

式中, Cv表示变异系数, SD为标准差, M表示研究时段多年平均植被覆盖度, 按照变异系数的划分等级, 将变异系数划分为弱变异(Cv＜= 0.1), 中等变异(0.1＜Cv＜= 1)以及强变异(Cv>1)。

2.4 相关性分析

相关性分析是在剔除控制变量的影响下, 分析指定变量之间是否存在显著的相关性。通过计算年均FVC与对应年均温度和年降水量的相关系数来反映FVC对不同气候因子的响应, 其计算公式为：

式中, n为监测时间段的累计年数, x_i、y_i分别为相关分析的2个变量, x、y为变量样本的均值。相关系数取值在[-1, 1]之间, 相关系数的绝对值越大, 说明植被覆盖度与气候因子的相关性越高。

3 结果分析 3.1 植被覆盖度增幅变化分析

2020年, 六大区域平均植被覆盖度为27.66%, 其中, 青藏高原生态屏障区、北方防沙带、南方丘陵山地带、长江重点生态区、东北森林带和黄河重点生态区分别为13.63%、13.42%、91.70%、76.99%、64.47%和41.74%。2000—2020年, 六大区域平均植被覆盖度增幅为21.5%, 比全国范围平均植被覆盖度增幅21.0%略高。六大区域平均植被覆盖度增幅差异明显：黄河重点生态区、北方防沙带和青藏高原生态屏障区增幅分别为54.4%、34.6%和21.8%, 是全国平均值的1.04-2.59倍, 黄河重点生态区增幅最高。长江重点生态区、东北森林带和南方丘陵山地带增幅分别为18.0%、13.7%和12.9%, 均低于平均值, 这3个区域2000年植被覆盖本底状况较好, 增幅变化相对不明显, 南方丘陵山地带增幅最低(图 2)。

从六大区域各等级植被覆盖度面积变化来看(见表 1), 青藏高原生态屏障区年均变化率为0.01/a, 以低植被覆盖度减少87732 km²和中植被覆盖度增加57918 km²变化为主。黄河重点生态区年均变化率为0.022/a, 以低植被覆盖度减少82352 km²和高植被覆盖度增加65808 km²变化为主。长江重点生态区年均变化率为0.008/a, 以中高植被覆盖度减少190843 km²和高植被覆盖度增加286512 km²变化为主。东北森林带年均变化率为0.007/a, 以中植被覆盖度减少131177 km²和中高、高植被覆盖度分别增加81726 km²、41470 km²变化为主。北方防沙带年均变化率为0.015/a, 低植被覆盖度减少136220 km², 中低、中、中高、高植被覆盖度面积分别增加32539 km²、25475 km²、67606 km²、10600 km², 单一变化不明显。南方丘陵山地带年均变化率为0.006/a, 以中高植被覆盖度减少62675 km²和高植被覆盖度增加81187 km²变化为主。综上分析, 黄河重点生态区年均变化最高, 并且主要为低植被覆盖度和高植被覆盖度的相互变化, 与该地区原本植被覆盖度较低, 在一系列生态工程实施下, 植被覆盖呈显著上升趋势^[34]变化基本一致。

3.2 植被覆盖度变化趋势分析

2000—2020年, 六大区域植被覆盖度变化总体呈增加趋势(图 3—8)。黄河重点生态区、东北森林带、南方丘陵山地带和长江重点生态区增加区域面积占比均超过60%, 北方防沙带和青藏高原生态屏障区植被增加区域面积占比相对较低, 分别为24.9%和17.6%, 同时这两个区域植被减少区域面积占比相对较高, 原因在于这两个区域范围面积较大, 且主要位于干旱、半干旱气候区, 植被变化受气候和人类活动影响较明显(图 9)。

21年间, 六大区域自然植被覆盖度显著增加和极显著增加区域面积明显高于显著减少和极显著减少区域(见表 2)。(1)青藏高原生态屏障区增加区域主要分布在该屏障区东部, 其中祁连山脉、三江源东部等区域极显著增加；减少区域主要以散点状分布于该区域东北部, 其中党河南山周边极显著减少较明显(图 10)。(2)黄河重点生态区增加区域基本遍布整个生态区；减少区域主要分布于该区域西北部和东部局部地区, 其中腾格里沙漠、河口区等区域较显著(图 11)。(3)长江重点生态区增加区域主要分布在中部和东部区域；减少区域零散分布在横断山脉南麓、洞庭湖及鄱阳湖周边区域(图 12)。(4)北方防沙带增加区域主要分布在河北北部及太行山东麓、内蒙古高原东侧；减少区域主要以零散点状分布于北方农牧交错生态脆弱区、腾格里沙漠南侧(图 14)。(5)东北森林带增加区域主要分布在长白山、小兴安岭地带、大兴安岭中北部；减少区域主要分布于三江平原和大兴安岭伊勒呼里山西侧(图 13)。(6)南方丘陵山地带增加区域基本遍布整个丘陵山地带；减少区域主要分布于金华江、郁江等河流沿岸周边(图 15)。

3.3 植被覆盖度稳定性分析

经统计分析(见表 3), 黄河重点生态区、长江重点生态区、东北森林带和南方丘陵山地带4个区域95%以上的范围植被覆盖度稳定性处于弱变异和中等变异, 而青藏高原生态屏障区和北方防沙带80%以上的范围植被覆盖度稳定性处于弱变异和中等变异, 较前4个区域范围相对有减少。总体来说, 2000—2020年六大区域自然植被覆盖度相对稳定。

从强变异空间位置来看(图 16—17), 主要分布在：(1)青藏高原生态屏障区的雅鲁藏布江中下游复合侵蚀脆弱区和中西部草原荒漠区, 远离城市建成区, 说明植被覆盖度波动变化的主要原因为自然因素, 受气候、地质灾害影响较大；(2)黄河重点生态区的半干旱、半湿润气候过渡区与北方防沙带的阿尔泰山西南侧以及阴山北部, 农耕、畜牧活动较多, 植被覆盖度变化的主要原因为人为因素。

4 结论与讨论 4.1 讨论

(1) 六大区域植被时空特征变化

从时间尺度上来看, 2000—2020年, 随着生态保护修复工程如退耕还林还草、天然林保护工程、石漠化治理等工程的陆续启动, 大大改善了六大区域植被生长环境并促进植被生长, 因此六大区域FVC整体呈上升趋势, 植被改善较为明显, 与赵安周等^[34]、郑朝菊等^[35]、曹亚楠等^[36]、阿多等^[37]的研究结果较一致。

从空间分布特征来看, 六大区域FVC变化的空间分布特征呈现由东南至西北递减的趋势, 表明已实施的生态修复工程在中国西南地区^[35]、西北地区^[36]、华北平原^[37]、黄土高原^[34]、长江流域^[38]等地区取得了成效。但仍有一些生态脆弱区域如青藏高原生态屏障区与北方防沙带受到人类活动的负面影响, 原有森林和草地面积减少, 生态系统部分退化^[39], 植被覆盖度显著退化面积较小, 这与郭建晓等^[40]与刘利民等^[41]的研究结果一致。

(2) 植被覆盖度变化与气温、降水的相关性

分别将六大区域植被覆盖度变化值与年均温度和年降水量进行相关分析(图 22—27), 结果显示青藏高原生态屏障区、长江重点生态区和南方丘陵山地带FVC与年降水量相关性不显著(P＞0.05), 但与年均温度显著相关(P＜0.05), 这些区域植被的生长主要受气温的限制, 植被覆盖度对温度的变化较敏感, 因此植被覆盖度的变化与年均温度的相关性更高, 这与相关学者的研究结果一致^[42—44]。北方防沙带植被覆盖度变化与年降水量显著相关(P＜0.05), 与年均温度相关性不显著(P＞0.05), 该区域处于干旱半干旱区域, 气候干燥, 降水稀少, 植被对降水变化的响应更明显, 这与沈贝贝等^[45]学者的研究结果一致。黄河重点生态区和东北森林带植被覆盖度变化与年均温度和年降水量相关性均不显著(P＞0.05), 这两个区域植被覆盖度的变化可能更受20世纪末以来中国在黄土高原等地实施的退耕还林和封山育林等一系列生态工程或其他人类活动的影响^[46—47]。

(3) 不足及下一步方向

本文在分析上, 采用的植被覆盖度数据的空间分辨率为500m, 虽然已胜任研究区大范围的植被覆盖度变化监测, 但低空间分辨率对于植被覆盖度的反映比较宏观, 目前缺少微观尺度的考量, 未来可以利用更高空间分辨率的数据获取更精准的变化信息, 以期能落实在主体实施工程的县域单元上。本文对FVC与对应年均温度和年降水量的相关性进行了分析, 但缺乏深入研究不同气候带、不同物候等地理环境背景及不同生态修复措施下的植被覆盖变化规律, 下一步将会深入分析这些区域植被覆盖度变化对地理环境、生态工程措施等的响应机制, 探讨自然因素和人为因素对生态修复背景下植被覆盖变化的贡献。

4.2 结论

本文基于研究区域内2000—2020年500m植被覆盖度数据, 从增幅变化、趋势分析、稳定性三方面分析了青藏高原生态屏障区、黄河重点生态区、长江重点生态区、东北森林带、北方防沙带、南方丘陵山地带等六大区域的自然植被覆盖度时空变化特征, 以期为后续保护修复工程实施选址提供基础支撑。

(1) 2000—2020年, 六大区域平均植被覆盖度增幅为21.5%, 比全国范围平均植被覆盖度增幅21.0%略高。黄河重点生态区、北方防沙带和青藏高原生态屏障区增幅分别为54.4%、34.6%和21.8%, 是全国平均值的1.04-2.59倍。长江重点生态区、东北森林带和南方丘陵山地带增幅分别为18.0%、13.7%和12.9%, 均低于平均值, 这3个区域2000年植被覆盖本底状况较好, 增幅变化相对不明显, 我国西北方较南方改善更明显。

(2) 2000—2020年, 六大区域植被覆盖度变化总体呈增加趋势。黄河重点生态区、东北森林带、南方丘陵山地带和长江重点生态区增加区域面积占比均超过60%, 北方防沙带和青藏高原生态屏障区植被增加区域面积占比相对较低, 分别为24.9%和17.6%, 植被变化受人类活动影响较明显。

(3) 六大区域自然植被覆盖度变化相对稳定。黄河重点生态区、长江重点生态区、东北森林带和南方丘陵山地带4个区域95%以上的范围植被覆盖度稳定性处于弱变异和中等变异, 而青藏高原生态屏障区和北方防沙带80%以上的范围植被覆盖度稳定性处于弱变异和中等变异, 较前4个区域范围相对有减少。强变异区域主要位于青藏高原生态屏障区、北方防沙带和黄河重点生态区的生态环境敏感脆弱区。

(4) 不同修复区域植被覆盖度变化与降水和温度的相关性有差异。青藏高原生态屏障区、长江重点生态区和南方丘陵山地带FVC与年降水量相关性不显著, 但与年均温度显著相关；北方防沙带植被覆盖度变化与年降水量显著相关, 与年均温度相关性不显著；黄河重点生态区和东北森林带植被覆盖度变化与年均温度和年降水量相关性均不显著。