1982—2022 年黄土高原归一化植被指数变化的时空特征

方德昆,闫小月,张学珍,罗海江

1.中国环境监测总站,国家环境保护环境监测质量控制重点实验室,北京 100012

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

植被覆盖是地表生态环境的重要组成部分,是保障生态质量的基础,也是实现“碳中和”的重要途径。 植被具有涵养水源、保持水土、固碳释氧和改善局地小气候的功能[1-2];此外,植被结构的改善有利于野生动植物种群的增加,对生物多样性的恢复起重要作用[3]。 因而,植被覆盖状况是区域生态环境质量评价的重要内容。

卫星遥感是区域植被覆盖监测的主要技术手段。 基于卫星遥感技术发展了一系列植被指数用于监测植被生长状况,如:归一化植被指数(Normalized Vegetation Index,NDVI)、比值植被指数(Ratio Vegetation Index,RVI)、增强型植被指数(Enhanced Vegetation Index,EVI)和土壤调节植被指数(Soil-Ad justed Vegetation Index,SAVI)等。其中,NDVI 物理意义清晰,数学计算过程简单,使用范围最广,成为我国生态环境质量评价的重要指标之一[4]。 与其他各指数相比,NDVI 具有产品完善、时间序列长等优点,尤其适合开展区域生态环境质量变化的长期研究。

黄土高原地处我国北方半干旱气候区,自21世纪初国家实施“退耕还林(草)工程”以来,黄土高原土地利用结构和格局发生了深刻变化。 特别是十八大以来,在习近平生态文明思想的指引下,生态恢复和生态建设逐步融入区域发展战略,山水林田湖草沙冰综合治理措施逐步实施,区域生态环境显著改善。 本文拟采用2 套卫星遥感的NDVI 数据产品(GIMMS 3g 和MODIS),通过数据融合的方法,建立一套长时间序列的NDVI 数据集,在此基础上揭示1982—2022 年黄土高原NDVI 变化的时空特征。

1 研究方法

本研究采用了GIMMS NDVI 3g 和MODIS NDVI 数据集。 GIMMS NDVI 3g 是美国航空航天局(NASA)全球监测与模型研究(GIMMS)项目发布的第三代基于NOAA 系列卫星AVHRR 传感器影像资料的NDVI 数据集,涵盖1981 年7 月至2015 年12 月的全球NDVI 数据,空间分辨率为1/12°(约8 km),时间分辨率为1 d。 数据生产过程中已处理了校准损失、轨道漂移、火山爆发等不利因素造成的干扰。 目前该数据集已被广泛应用于全球各地的长期植被覆盖变化研究[5-7]。

MODIS NDVI 数据集是NASA 地球观测系统(EOS)计划发布的全球NDVI 数据集,由地球遥测卫星Terra 和Aqua 上的中分辨率成像光谱仪(MODIS) 影像资料合成。 本研究采用编号为MOD13A2 v6 的数据产品,该产品使用基于一定合成周期的最大值合成法消除云层和气溶胶的干扰噪声,涵盖2000 年2 月至2023 年2 月的全球NDVI 数据,空间分辨率为1 km,时间分辨率为16 d。 目前该数据集已被广泛应用于植被动力学研究[8-9]。

由于GIMMS NDVI 3g 和MODIS NDVI 数据集的时间跨度均不足以单独覆盖目标研究时段,故需融合2 套数据重构1982—2022 年的NDVI时间序列。 现有研究表明,GIMMS NDVI 3g 和MODIS NDVI 之间存在高度相关性[10-11],故本研究采用方差匹配法[式(1)和式(2)]重构覆盖1982—2022 年且均一化的NDVI 时间序列。 方差匹配法的核心思想是通过调整其中一个数据集,使2 个时间序列重叠时段数据的平均值和标准差保持一致,从而确保重构序列前后2 段的方差可比。

式中:Mi表示i 年的MODIS NDVI;和分别表示重叠时段(即2003—2013 年)的MODIS NDVI和GIMMS NDVI 3g 的平均值; SM和SG分别表示重叠时段MODIS NDVI 和GIMMS NDVI 3g 的标准差;表示2003—2022 年MODIS NDVI 的平均值。

由于GIMMS NDVI 3g 的时间跨度比MODIS NDVI 更长,本研究以GIMMS NDVI 3g 为参考标准,调整MODIS NDVI 的平均值和标准差,使其在重叠时段与GIMMS 3g 的平均值和标准差一致。由于数据可用性原因,本研究采用2 个数据集的重叠时段2003—2013 年作为参考时段进行数据调整,之后通过组合1982—2002 年的GIMMS NDVI 3g 数据和调整后的2003—2022 年MODIS NDVI 数据,构建1982—2022 年均一化的NDVI时间序列。 此外,由于2 个数据集的空间分辨率不同,故在应用方差匹配法进行调整之前,需对MODIS NDVI 进行像元合并。 本研究采用面积平均法进行此升尺度操作,将MODIS NDVI 的空间分辨率由1 km 升至8 km,确保与GIMMS NDVI 3g 一致。

利用均一化的1982—2022 年NDVI 时间序列,首先分析逐像元的NDVI 变化趋势,为突出“退耕还林(草)工程”以及十八大以来NDVI 变化的空间特征,针对1982—2019 年、1982—2002年、2003—2022 年、2012—2022 年4 个时段分别进行分析,变化趋势采用最小二乘法进行拟合。采用土地利用类型变化情况解释黄土高原NDVI变化的原因,其中土地利用数据空间分辨率为1 km,来自中国多时期土地利用遥感监测数据集[12],包括1980、2000、2010、2020 年4 个时间点。 最后,以5 a 为一个时间段,对比分析了各时段之间各等级NDVI 的面积变化及主要NDVI 等值线的空间位置变化。

2 结果与讨论

2.1 区域NDV I 变化的时空特征及与土地利用变化的关系

1982—2022 年黄土高原暖季(5—9 月)NDVI时序变化如图1 所示。 黄土高原暖季NDVI 呈现出明显的增加趋势,增长率约为0.03/(10 a),暖季NDVI 在2018 年达到最大值(0.55),比最小值(0.40,出现在1982 年)高出约38%。 黄土高原暖季NDVI 具有明显的年际波动趋势,这与年际的降水不均匀、旱涝、气候变化有关联[13-14]。 其增加趋势表现出“先慢后快”的特征,具体可划分为2 个明显的时间段:1982—2002 年,暖季NDVI表现为缓慢的上升趋势,速率为0.01/(10 a),增幅仅为0.06;2003—2022 年,暖季NDVI 快速上升,速率高达0.06/(10 a), 增幅约为0.12。1982—2002 年,暖季NDVI 均值仅为0.44,但由于暖季NDVI 持续上升,均值在2012—2022 年达到较高水平,为0.52。 暖季NDVI 的增加反应了植被覆盖率的提高,表明近40 年来黄土高原地表持续绿化,2003—2022 年绿化效果更为明显,又以2012—2022 年尤为显著。

图1 1982—2022 年黄土高原暖季(5—9 月)NDVI 时序变化Fig.1 NDVI’s temporal tendency in the warm season (May-September) of the Loess Plateau from 1982 to 2022

黄土高原各时间段暖季NDVI 变化趋势的空间分布如图2 所示。 总体来看,近40 年来黄土高原86.4%区域的暖季NDVI 呈上升趋势,其余区域略有下降[图2(a)]。 其中黄土高原中部,即山西西部、陕西北部和宁夏中部暖季NDVI 上升速率较快,约为0.03/(10 a) ;这些地区主要为沙漠(毛乌素沙漠)和山区(吕梁山脉)。 暖季NDVI 呈下降趋势的区域面积较小且分布零散,最大下降速率约为-0.02/(10 a) 。 通过与同期土地利用类型的空间分布[图3(a)、3(d)]比较,可以发现暖季NDVI 上升速率较快的区域与同期耕地转换成林地、草地的空间范围基本一致;而在暖季NDVI 下降趋势的区域中,主要的土地利用变化则是林、草地转变为耕地及建设用地。

图2 黄土高原暖季NDVI 变化趋势空间分布Fig.2 Spatial distribution of NDV I’s temporal tendency in the warm season of the Loess Plateau

图3 黄土高原土地利用类型空间分布Fig.3 Spatial distribution of land use types in the Loess Plateau

以2002 年“退耕还林(草)工程”政策付诸实施为时间节点,之前(1982—2002 年)黄土高原暖季NDVI 变化趋势并不明显(图1)。 虽然70.9%的区域暖季NDVI 呈上升趋势,但上升速率较小,仅为0.01/(10 a) 。 暖季NDVI 呈下降趋势的区域面积较小,主要分布在山西中西部、河南与山西交界处以及陕西中部,且下降速率也较小,大多不超过-0.01/(10 a) 。 同期土地利用变化甚微[图3(a)、图3(b)]。 “退耕还林(草)工程”政策实施之后(2003—2022 年),黄土高原93.8%区域的暖季NDVI 呈上升趋势[图2(c)]。 其中,自黄土高原东北至西南方向出现一条醒目的“绿化带”,横穿内蒙古南部、山西北部、陕西中部、甘肃东南部和宁夏南部,暖季NDVI 上升速率平均高达0.04/(10 a) 。 同时期,黄土高原发生了大面积的“耕地范围缩减、林草范围扩增”的土地利用变化[图3(b)、图3(d)]。 可见,“退耕还林(草)工程”是近40 年黄土高原暖季NDVI 增加的重要原因。

十八大提出“加强生态文明建设”的总体要求以来,黄土高原暖季NDVI 的增加趋势愈加显著(图1)。 其中黄土高原东北部、西南部,即山西北部、陕西北部和甘肃东部上升速率较大,约为0.07/(10 a) ,其空间分布与2003—2022 年暖季NDVI 上升速率高值区基本吻合。 因此,2012—2022 年黄土高原暖季 NDVI 的显著上升是

2003—2022 年黄土高原地表绿化的主要成因。同时期,暖季NDVI 呈下降趋势的区域主要包括山西南部的运城盆地和洛阳盆地,平均下降速率约为-0.05/(10 a) 。 通过与同期土地利用变化的空间分布对比[图3(c)、图3(d)],可以发现暖季NDVI 增加区域与耕地向林、草地转换区域范围基本一致,暖季NDVI 降低区域则与城市、耕地扩张区域基本一致。 可见,耕地继续向林地、草地转换可能是这一时期暖季NDVI 呈上升趋势的主要成因;而城建用地和耕地范围的扩增,则可能是黄土高原东南部暖季NDVI 呈下降趋势的主要成因。

由上述结果可知,黄土高原暖季NDVI 变化与同期土地利用变化密切相关。 由于2003—2022 年暖季NDVI 增速较大,2012—2022 年增速尤为显著,因此选取黄土高原2000、2010、2020 年土地利用类型,分别对2003—2022 年、2012—2022 年暖季NDVI 增速大于0.05/(10 a) 区域的土地利用变化情况进行分析,旨在进一步认识暖季NDVI 快速增加背后的土地利用变化因素。 其中,2000—2020 年土地利用变化对应于2003—2022 年暖季NDVI 变化,2010—2020 年土地利用变化对应于2012—2022 年暖季NDVI 变化。 将土地利用类型转换面积由大至小排列,取前6 位,如图4、表1 所示。 总体来看,在暖季NDVI 增速大于0.05/(10 a)区域内,土地利用类型变化主要表现为耕地向草地、林地的转换。 2000—2020年,耕地转换为草地、耕地转换为林地、草地转换为林地的面积共计6.70×104km2,占同期土地利用变化总面积的43.9%,其中耕地转换为草地的面积百分比最大(为29.3%)。 2010—2020 年土地利用变化特征与2000—2020 年基本一致,均表现为“耕地面积减少、林草面积增加”。 可见,这种变化趋势是2 个时期暖季NDVI 快速增加的主要原因。 进一步观察2010—2020 年土地利用类型变化的空间分布情况可知[图4(b)],黄土高原东部暖季NDVI 增加的主要原因可能是草地向林地的转换,而其他区域则是耕地向草地的转换。

表1 2003—2022 年、2012—2022 年黄土高原暖季NDVI 增速大于0.05/(10 a) 区域的土地利用类型变化情况Table 1 Areas of each type of land use changes for regions with NDVI’s grow th rate greater than 0.05/(10 a) in the warm season of the Loess Plateau from 2003 to 2022 and from 2012 to 2022,respectively

图4 2003—2022 年、2012—2022 年黄土高原暖季NDV I 增速大于0.05/(10 a) 的区域及其土地利用类型变化Fig.4 Areas with NDVI’s growth rate greater than 0.05/(10 a) and their land use changes in the warm season of the Loess Plateau from 2003 to 2022 and from 2012 to 2022,respectively

2.2 典型暖季NDVI 等值线变化的主要特征

为了进一步探究近40 年来黄土高原暖季NDVI 的动态变化情况,以5(a)为一个时间段,对比分析了各时段各等级暖季NDVI 的面积变化及暖季NDVI 等值线的空间位置变化(图5)。 总体来看,以暖季NDVI 等值线0.4 为界,黄土高原西北区域等值线数值较小且分布稀疏,这是由于西北区域分布有沙漠(乌兰布和沙漠、库布齐沙漠和毛乌素沙漠)及高原(鄂尔多斯高原),土地利用类型以草地和荒漠为主,植被覆盖度较低。 而东南区域等值线数值较大且分布密集,这是由于该区域分布着山脉(吕梁山脉、子午岭山脉和六盘山脉)及盆地(临汾盆地、运城盆地和洛阳盆地),栽培植被、阔叶林广布,植被覆盖度较高。

图5 黄土高原暖季NDVI 等值线分布Fig.5 Distribution of NDVI contours in the warm season of the Loess Plateau

以5(a)为一个时间段,即1985—1989 年、1995—1999 年、2005—2009 年、2015—2019 年,通过前、后2 个时段的对比来看,暖季NDVI 等值线的空间位置发生了明显变化,主要表现为:低值等值线(0～0.4)由中部向西北方向推进(图5),导致暖季NDVI 小于0.4 的面积减少了11.81×104km2(表2);而高值(0.6 ～0.1)等值线高值由东南向中部推进,导致暖季NDVI 大于0.6 的面积增加了6.81×104km2。 谢佩君等[15-16]通过分析1987—2020 年的“植被绿度线”位置的变化发现,植被绿度线明显向西北方向移动,尤其是2005—2010 年北移幅度最大,在空间格局上,这一变化主要集中在黄土高原中部区域的黄土丘陵区,即陕西省北部榆林和延安一带。 这表明黄土高原中部地区暖季NDVI 持续增加,植被恢复效果最为明显。 此外,暖季NDVI 等值线空间位置变化、各值域面积变化在各时间段并不相同,即在“退耕还林(草)”前(1985—1989 年、1995—1999 年),暖季NDVI 等值线空间位置变化不明显、各值域面积变化较小;“退耕还林(草)”后(2005—2009 年、2015—2019 年),暖季NDVI 等值线空间位置变化较为明显,各值域面积变化较大,与前文分析结果一致。

表2 黄土高原暖季NDVI 各值域区间的面积Table 2 Areas for each range of NDVI in the warm season of the Loess Plateau

黄土高原植被覆盖度增加致使生态质量有所恢复转好,植被覆盖度和植被净初级生产力呈增加趋势,其中水源涵养服务持续转好的面积占比为15.46%[17]。 2000—2008 年,黄土高原生态系统净固碳能力从0.011 Pg 增至0.108 Pg,植被固碳量增加速率为每年9.4 g/m2[18]。 值得注意的是,暖季NDVI 在黄土高原南部、东部和青海东部一带呈缓慢下降趋势。 这可能是由于:黄土高原南部、东部人口较多,1990—2015 年人口增幅为29.83%[19];区域城镇化发展迅速,耕地分布较为广泛,受人类活动的影响,植被覆盖发生一定程度的退化[20]。 青海东部是黄土高原向青藏高原的过渡地带,干旱缺水的自然条件使农业生产较难达到可持续发展的必要条件[21]。 局部植被覆盖度减少对生态系统服务功能构成了一定威胁[18]。

3 结论

采用方差匹配方法对GIMMS NDVI 3g 和MODIS NDVI 进行了融合,重构形成了1982—2022 年黄土高原的NDVI 时间序列。 研究发现其间黄土高原暖季NDVI 呈显著上升趋势,其中1982—2002 年黄土高原暖季NDVI 出现小幅度的增加,2003—2019 年(尤其是十八大以来)黄土高原暖季 NDVI 上升趋势明显, 且增速较快。2003—2019 年的绿化区域主要分布在黄土高原中部,并向东北、西南方向延展,这些区域是开展“退耕还林(草)”的重点区域;同时,暖季NDVI等值线低值向西北方向明显推进。 说明近20 年的生态建设工程取得了明显成效,地表植被覆盖度大幅提升。 值得注意的是,最近20 年间黄土高原南部、东部和青海东部一带,由于城市扩展、气候变化等因素,NDVI 在部分区域稍有下降。