基于GEE 的大青山国家级自然保护区NDVI变化及影响因素分析

艾丽亚， 王永芳,2,3， 郭恩亮,2， 银 山,4， 顾锡羚

（1.内蒙古师范大学地理科学学院，内蒙古 呼和浩特 010022；2.内蒙古自治区蒙古高原灾害与生态安全重点实验室，内蒙古 呼和浩特 010022；3.蒙古高原气候变化与区域响应高校重点实验室，内蒙古 呼和浩特 010022；4.内蒙古自治区遥感与地理信息系统重点实验室，内蒙古 呼和浩特 010022）

自然保护区建设是我国生态文明建设的重要组成部分，在维持生物多样性、涵养水源和固碳等方面起着关键作用［1］。尤其是近些年来，我国在国家级自然保护区建设和管理层面上取得了巨大发展，截至2018 年，已建立国家级自然保护区474 处，对区域自然资源合理利用和生态环境的恢复与保护方面均起到了积极作用［2］。由此，自然保护区生态环境效益和保护成效评估研究也成为诸多学者所关注的热点问题之一［3-4］。

植被动态变化研究是区域生态环境效益评估的有效手段。归一化植被指数（Normalized difference vegetation index，NDVI）作为反映地表植被生长状况的重要指标，能够定性、定量评价植被生长和覆盖状况，得到了广泛应用［5-7］。因高分辨率遥感影像使用门槛较高，已有研究大多采用中、低分辨率影像数据开展长时序的植被NDVI 变化研究［8-9］，但此举在小区域尺度往往无法满足精度上的需求。Landsat是一种常见的中分辨率遥感影像，具有时间跨度长、多波段等优势，在免费开放的影像中具有较高的空间分辨率，在植被NDVI 获取方面具有较高的实用价值。然而，因该影像的处理过程复杂，较难实现长时间序列的植被监测工作［10］。谷歌地球引擎（Google Earth Engine，GEE）云平台，兼顾了数据收集下载和强大的空间分析处理能力，能够支持长时间序列的遥感影像获取和分析，实现了数据获取、处理、分析与应用的一体化，为长时序、中高分辨率植被遥感信息提取提供了新的技术手段［11-12］。

阴山山脉是中国北方重要的自然和生态界山，其中大青山地处阴山山脉中段，是生态多样性最丰富的区域之一。该保护区具有涵养水源、水土保持、防风固沙、固碳释氧等生态功能，也是黄河上中游重要的水源补给区，在调节我国北方水分平衡和水资源供给中起着重要作用［13］。然而，气候干旱、土地贫瘠，加之人为过度开发利用导致该区生态环境极其脆弱［14］。为了扭转此局面，政府部门于2000年在内蒙古自治区成立了自然保护区，在2008年设立了国家级自然保护区，并实施了一系列生态环境保护工程。但目前有关该区生态环境效益评估还未得到有效开展。因此，本文借助GEE 平台，采用Landsat 影像对比分析了大青山国家级自然保护区建立前后的NDVI 时空演变特征，并定量评估了气候因子和人类活动对植被变化的影响。研究成果可以探清大青山国家级自然保护区成立对该区带来的生态效益，同时可为保护区未来生态环境管理决策提供基础信息与技术支撑。

1 研究区概况

内蒙古大青山国家级自然保护区位于土默特川平原北部的呼和浩特市、包头市和乌兰察布市境内（109°47′～112°17′E，40°34′～41°14′N），总面积约3888.70 km2，东西长度为217 km，南北平均宽度为18 km。保护区属中低山地，山区海拔为1012～2319 m，北坡平缓南坡陡峭（图1a）。该区属于温带半干旱大陆性季风气候，年均气温3～4 ℃，年降水量400～500 mm。保护区北坡易受蒙古高原干燥气流的影响，气候寒冷、干燥；南坡由于山地的阻挡，气候温暖、湿润。该保护区是山地森林、灌丛、山地草原和重要水源涵养地为主的综合性山地森林类保护区［13］。本文采用的保护区界线以国务院批准的文件为准［15］，保护区按照不同功能，可分为核心区、缓冲区和实验区（图1b）［16］。

图1 研究区概况Fig.1 Overview of the study area

2 数据与方法

2.1 数据来源

本文所用的遥感数据是Landsat-5 TM、Landsat-7 ETM+和Landsat-8 OLI 时间序列数据。选择1995—2020 年4—10 月植被生长季云量低于10%的影像（图2），利用GEE 平台对数据进行格式转换、数据拼接与裁剪，并计算NDVI值［17］。由于Landsat-8 卫星传感器和Landsat-5、7 卫星传感器之间有明显差异，根据Roy等［18］提供的截距和偏移调整Landsat-8 各个波段的反射率，实现Landsat-5、7、8 协同，并以最大值合成法得到1995—2020年NDVI年值数据集，空间分辨率为30 m［19］。该数据集NDVI 最大值均集中于7 月，最大值在0.89～1.00 之间，年际变化可对比性较强。

图2 Landsat影像具体日期Fig.2 Specific date of the Landsat images

本文采用的气象数据为中国地面气温降水日值0.5°×0.5°格点数据集（V2.0）。选取研究区1995—2020 年4—10 月生长季的9 个像元数据，利用Python对数据进行反距离权重插值，并通过平均、求和获取了研究区气温、降水量年值数据集，空间分辨率为30 m。

土地利用数据源自武汉大学杨杰、黄昕教授基于GEE、利用Landsat 数据提取的1985—2020 年中国土地覆盖数据集，空间分辨率为30 m［20］。本文提取了1995—2020年的草地和林地面积，用于表征研究区生态环境工程实施成效。

2.2 研究方法

2.2.1 Sen氏斜率与Mann-Kendall非参数检验Sen氏斜率是非参估算线性趋势的算法［21］，本研究用于获取1995—2008 年和2008—2020 年2 个时间段植被NDVI的变化趋势，其计算公式为：

式中：β为趋势度；xi和xj分别为某一像元在i年和j年的值。使用趋势度β判断时间序列趋势的升降，当β＞0 时，表示植被NDVI 呈增加趋势；当β＜0时，则呈减少趋势。

本文采用Mann-Kendall方法在0.05置信水平上对NDVI的变化趋势进行显著性检验［22］。

2.2.2 相关性分析本文采用皮尔逊相关系数测定NDVI与气温和降水量之间的相关关系［23-24］，计算公式为：

式中：rxy为相关系数；n为研究时段年数；x、y为相关分析的2 个变量；xi、yi分别为其样本值；xˉ为气候因子的平均值；yˉ为NDVI 的平均值。本文同样在0.05置信水平上对相关性进行显著性检验［25］。

2.2.3 残差分析本文假设NDVI 仅受人类活动和气候条件影响，并利用残差分析法评估人类活动对NDVI的影响［26］。计算公式为：

式中：y为因变量NDVI 值；x1、x2为自变量；a为 常数；b、c为回归系数；σ为残差值；NDVI实际和NDVI预测分别为实际NDVI 值和预测NDVI 值。残差趋势为正表示人类活动对NDVI 起促进作用，为负表示人类活动对NDVI起抑制作用［27］。

2.2.4 Lindeman-Merenda-Gold模型相关性分析和残差分析可以反映气候因子与人类活动对NDVI的影响，但无法量化因子对NDVI 变化的贡献率［28］。因此本文通过多元线性回归模型的方差分解量化各变量对NDVI 变化的相对重要性。并选择计算最为密集的Lindeman-Merenda-Gold（LMG）模型［29］，区分多元线性回归中不同相关回归量的贡献率，由R语言包“relaimpo”实现。贡献率越大说明变量对植被变化的影响越大，反之则相反［30］。选取的气候因子包括气温和降水量；人类活动由林地和草地面积表征，即人工造林、封山育林、退耕还林还草和生态移民等生态环境工程所带来的成效。

3 结果与分析

3.1 NDVI时空变化特征

3.1.1 年际变化特征如图3所示，1995—2008年保护区植被NDVI 斜率为-0.002 a-1，显著性检验Z值为-0.77，说明NDVI 整体呈非显著下降趋势。从NDVI 值的变化看出，2001—2003 年NDVI 出现持续增加，且增速较快；2008—2020 年研究区植被NDVI呈上升趋势，斜率为0.011 a-1，Z值为2.20，表明该时间段内NDVI 呈显著上升趋势。整体而言，保护区成立以前NDVI 呈非显著下降趋势。2008 年国家级自然保护区成立后NDVI 呈显著上升趋势，说明国家级保护区的成立促进了研究区植被恢复。

图3 NDVI年际变化特征Fig.3 Characteristics of interannual variation of NDVI

3.1.2 空间变化特征1995—2008 年保护区NDVI减少区域共占总面积的69.04%，其中显著减少区域占7.78%，主要分布于保护区东西延伸的山脉以北地区（图4a～b、表1）。NDVI 增加区域占总面积的19.90%，其中显著增加区域占1.55%，主要分布于保护区南部沿线地区。说明在国家级保护区成立前NDVI 变化以非显著减少为主导，呈现“北部减小、南部增长”的空间变化特征。分区方面，1995—2008 年核心区、缓冲区和实验区的非显著减少区域占比分别为60.26%、64.97%和60.29%，显著减少区域占比分别为8.14%、6.64%和8.04%，说明3个分区植被NDVI的变化程度较为接近。

表1 不同时期NDVI变化趋势统计Tab.1 NDVI variation trend statistics in different time periods

图4 NDVI变化趋势及显著性检验空间分布Fig.4 Spatial distributions of NDVI variation trend and significance test

2008—2020 年NDVI 减少区域占总面积的2.26%，其中显著减少区域仅占0.10%，分布于保护区中心地带。而NDVI 增加区域占总面积的94.98%，其中显著增加区域占比为38.47%，在保护区东、西部地区较为明显（图4c～d）。说明国家级自然保护区成立后植被NDVI 整体呈现出向好的特征。核心区、缓冲区和实验区的显著增加区域分别占其总面积的30.35%、36.19%和43.88%。说明在2008 年自然保护区成立之后，实验区的植被恢复高于其他功能区。

3.2 气候因子对NDVI的影响

3.2.1 气候因子的年际变化特征1995—2008 年保护区降水量的变化速率为-4.737 mm·a-1，整体呈非显著下降趋势（图5a）；气温则以0.053 ℃·a-1的速率呈非显著上升趋势（图5b），说明该时期保护区气候较为干旱，不利于植被生长与恢复。其中2001—2003 年降水量持续增加、气温持续降低、气候相对湿润为该时间段内NDVI的增长提供了有利条件。

图5 气候因子年际变化特征Fig.5 Characteristics of interannual variation of climatic factors

2008—2020年降水量变化速率为3.059 mm·a-1，整体呈非显著上升趋势。气温上升速率由0.053 ℃·a-1减缓为0.019 ℃·a-1，说明气候暖干化趋势有所缓解。尤其是降水量的增加对植被生长、恢复起到促进作用，使该时期NDVI增加。

3.2.2 气候因子与植被NDVI的相关关系由图6和表2 可知，2 个时间段内保护区降水量与NDVI 的相关性在空间上有所不同。1995—2008 年期间，NDVI 与降水量的相关系数以正值为主，其中非显著正相关区域占79.06%。显著正相关区域占研究区总面积的18.54%，主要集中于保护区西部和东北部，在中部地区也有零星分布（图6a～b），说明在这些地区降水量对植被的影响较大。核心区、缓冲区和实验区的显著正相关区域分别为15.86%、16.90%和20.69%，说明降水量对实验区的影响大于核心区和缓冲区。

表2 NDVI与降水量之间的相关性统计Tab.2 Statistic for correlation between NDVI and precipitation

图6 NDVI与降水量之间的相关系数及显著性检验空间分布Fig.6 Spatial distributions of correlation coefficient and significance test between NDVI and precipitation

保护区成立后，NDVI 与降水量呈非显著正相关区域占总面积的72.08%。显著正相关区域占比为24.57%，主要集中于保护区西部，在南部沿线和东北部地区也有体现（图6c～d），说明降水量对NDVI 的影响范围在扩大。并且，这些地区植被NDVI在该时间段内显著增加（图6d），也表明了降水量的增加是植被恢复的主要原因之一。核心区、缓冲区和实验区的显著正相关区域分别为19.13%、19.58%和29.61%，说明降水量对实验区的影响仍大于核心区和缓冲区。

如图7 和表3 所示，1995—2008 年保护区气温与NDVI 以负相关为主，其中非显著负相关区域占研究区总面积的56.17%。显著负相关区域占比为41.75%，在研究区东北、中部和西部地区均有体现（图7a～b），说明气温的升高导致这些区域植被NDVI的减少。分区方面，核心区、缓冲区和实验区的显著负相关区域占比分别为41.01%、41.47% 和42.27%，说明气温对3个功能区的影响较为一致。

表3 NDVI与气温之间的相关性统计Tab.3 Statistic for correlation between NDVI and temperature

图7 NDVI与气温之间的相关系数及显著性检验空间分布Fig.7 Spatial distributions of correlation coefficient and significance test between NDVI and temperature

2008—2020 年气温与NDVI 以非显著正相关为主，占保护区的84.46%，显著正相关区域仅占1.58%（图7c～d），说明该时间段内气温对NDVI 无显著影响。核心区、缓冲区和实验区显著正相关区域分别为2.40%、2.15%、0.89%，表明气温对实验区NDVI的影响最小。

3.3 人类活动对植被NDVI的影响

根据国家相关规定，保护区内的核心区和缓冲区禁止人类活动。因此，本文以实验区为例，评估人类活动对NDVI 的影响。图8 和表4 所示，1995—2008 年残差趋势为抑制NDVI 的区域占实验区总面积的33.84%，主要分布于北部沿线地区，均为NDVI下降或显著下降地区，说明不合理的人类活动导致实验区北部NDVI 的下降；残差趋势为促进NDVI 的区域占实验区总面积的51.57%，主要分布在中部与南部（图8a），与NDVI 呈增加趋势的地区分布较为一致（图4a），说明正向的人类活动促进了该区植被NDVI的增加。

表4 实验区NDVI残差趋势统计Tab.4 Statistics of NDVI residual trends in the experimental area

图8 实验区NDVI残差趋势空间分布Fig.8 Spatial distributions of NDVI residual trends in the experimental area

2008—2020 年实验区残差趋势为抑制NDVI 的区域占比仅为2.79%；而残差趋势为促进NDVI的区域明显增加，覆盖实验区94.62%的面积，与图4b 中NDVI 呈增加趋势的区域基本重叠（图8b）。说明2008年国家级自然保护区成立以后人类活动对NDVI 的影响以促进作用为主，并驱动了对实验区植被的恢复。

3.4 气候因子和人类活动对NDVI变化的贡献率

为了进一步量化气候因子和人类活动对NDVI变化的影响，同时也为了找出保护区成立前后植被NDVI 变化的主导因素，本文利用LMG 模型计算了气温、降水量、林地面积和草地面积对NDVI 的贡献率。结果显示，1995—2008年降水量对实验区NDVI变化的贡献率为32.00%，气温的贡献率为58.42%，林地和草地面积的贡献率不足10%（图9a），说明气候暖干化趋势是该时间段实验区NDVI 减少的主要影响因素；2008—2020 年林地面积对实验区NDVI变化的贡献率为34.28%，草地面积的贡献率为43.93%，降水量和气温对NDVI 的贡献率分别为19.02%和2.78%（图9b）。说明大青山国家级自然保护区成立后，正向的人类活动为实验区NDVI 变化的主要影响因素，也说明实验区内实施的生态环境保护工程取得了显著效益。

图9 气候因子和人类活动对实验区NDVI变化的贡献率Fig.9 Contribution of climate factors and human activities to NDVI change in the experimental area

4 讨论

本文借助GEE云平台，利用1995—2020年Landsat遥感影像，提取了研究区30 m 分辨率植被NDVI 数据集，对大青山国家级自然保护区长时序植被动态变化进行了监测。结果显示在国家级自然保护区成立后，NDVI 呈增加趋势的地区占研究区总面积的94.98%，植被生态环境大幅好转，这与胡尔查等［31］和王佳新等［16］的研究结果较为一致。因选用了Landsat 遥感影像，且借助了GEE 集成运算环境，本文获取了时间跨度长、空间分辨率较高的植被监测结果，更为精确地刻画了植被NDVI 的时空演变特征，丰富了长时序、小区域尺度植被监测技术手段。

本文研究结果显示气候的暖干化是国家级自然保护区成立前NDVI 减少的主要原因，这与陈淑君等［7］的研究结果较为一致。气温升高、降水减少容易导致气象干旱的发生，从而使得大气、土壤水分亏缺，加剧植被的干旱胁迫，抑制植被生长［32］。本文还发现大青山国家级自然保护区成立后，正向的人类活动是实验区NDVI 变化的主要影响因素。据不完全统计，保护区已累计完成人工造林6.45×105hm2，完成封育面积1.32×106hm2，这些工程措施使得大青山国家级自然保护区森林覆盖率大面积提高［33］。另外，以往的研究大多采用土地利用类型转移面积，如耕地转化为林地和草地面积［34-35］表征生态工程实施成效。而本文直接采用草地和林地面积予以表征，原因在于研究区耕地面积较少［36］，耕地转为林草地面积仅能表征小部分生态环境工程效益。

本研究聚焦于国家级自然保护区设立对生态环境产生的正向作用，研究结果也证明了生态工程实施后植被向好发展的事实。但是由于缺乏放牧强度、砍伐树木等表征负向人类活动的时空数据，本研究未能将这些指标应用到LMG 模型的运算。也因缺乏保护区生态环境工程实施的准确信息，未能紧密结合生态工程实施的具体时间、方位等开展深层次的人类活动对植被NDVI 变化的驱动机制研究。此外，本文只讨论了气温和降水量两种气候因子对NDVI 的影响。内蒙古大青山国家级自然保护区作为半干旱地区，植被对气候因子的响应关系是一个复杂的过程，植被变化还会受到潜在蒸散等其他气候因子的影响［37］。因此，在未来将考虑多种气候因子对植被NDVI的综合影响。

5 结论

（1）1995—2008年保护区成立之前，植被NDVI以-0.002 a-1的速率减少，减少区域共占研究区总面积的69.04%。2008—2020 年保护区成立以后植被NDVI 以0.011 a-1的速率增加，增加区域占总面积的94.98%，说明保护区的建立对植被的好转起到明显促进作用。

（2）1995—2008 年研究区气候呈暖干化趋势。NDVI 与气温呈显著负相关区域面积占比为41.75%，与降水量呈显著正相关区域面积占比为18.54%。2008—2020 年NDVI 与降水量、气温均以正相关为主，显著正相关区域分别为24.57%和1.58%，说明在该时间段内降水量对植被NDVI 的影响有所增加，而气温的影响则显著减弱。

（3）1995—2008 年降水量和气温对实验区NDVI 变化的贡献率分别为32.00%和58.42%，说明气候因子是实验区NDVI 变化的主要影响因素。2008—2020 年实验区林地面积和草地面积对NDVI 变化的贡献率分别为34.28%和43.93%，说明正向的人类的活动是植被恢复的主导因素。