基于GEE平台的1986—2021年黄土高原植被覆盖度时空演变及影响因素

赵安周，田新乐

河北工程大学矿业与测绘工程学院，河北 邯郸 056038

植被是陆地生态系统的第一性生产者和主要组成成分，不仅在全球物质循环和能量流动等方面具有关键的作用，而且对气候变化和人类活动响应敏感（李依璇等，2021；田智慧等，2022）。监测长时间序列植被的时空格局变化及其驱动机制对区域生态恢复及可持续发展具有重要的意义（赵安周等，2017；袁和第，2020）。

遥感技术具有宏观性、快速性、准确性以及能够进行动态监测等优势，目前已被广泛应用在长时间序列和不同空间尺度的植被动态监测中（Xu et al.，2022）。诸多遥感数据中，NDVI由于计算简单且与叶面积指数（Leaf Area Index，LAI）和净初级生产力（Net Primary Productivity，NPP）等植被指数密切相关而被广泛应用在区域及全球植被动态变化监测中（张宝庆等，2011；张园等，2020；黄栋等，2021）。如张园等（2020）利用 2001—2018年 MODIS NDVI数据，探究长白山自然保护区植被对温度和降水变化的响应。张宝庆等（2011）结合SPOT NDVI和GIMMS NDVI两种数据集对近30年黄土高原植被变化趋势及空间分布进行监测。黄栋等（2021）使用SPOT NDVI数据集，从区域和像元两个角度分析了环渤海地区植被空间动态演变趋势及其对气候和土地利用变化的响应。需要指出的是，NDVI会受到大气和土壤等背景信息的影响，不仅在高植被区域存在过饱和现象，而且难以区分低植被覆盖区的植被信息（杨嘉等，2008）。因此也有学者利用像元二分法模型估算植被覆盖度（Fraction of Vegetation Coverage，FVC）监测植被的时空演变格局及驱动因素（张家政等，2022）。目前的研究表明该指标可以更加准确表征植被状况，反映陆地生态系统的健康情况（易浪等，2014）。

黄土高原地处黄河中游，是构建国家“两屏三带”生态安全格局和落实黄河流域高质量发展的主体区域（李依璇等，2021）。长期以来，由于气候变化和不合理的人类活动，使该地区水土流失严重，是黄河主要的泥沙输入地（赵安周等，2017）。自1999年实施生态工程建设以来，该地区的植被得到了快速恢复，生态环境明显改善（袁和第，2020）。目前已有诸多研究利用NDVI、FVC和NPP等指标对该地区的植被时空演变格局及其影响因素进行了分析，易浪等（2014）、Han et al.（2022）和张家政等（2022）分别采用MODSI NDVI、SPOT NDVI和GIMMS NDVI等数据分析了黄土高原植被时空演变及影响因素，认为该地区植被呈显著增加的趋势，且对气候变化敏感。受限于本地计算资源和数据分发的方法，目前研究采用的 MODIS NDVI时间从2000年开始，空间分辨率为250 m、500 m或1 km；SPOT NDVI时间从1998年开始，其空间分辨率为1 km；而GIMMS数据的获取时间虽从1981年开始，但其空间分辨率仅为0.08333°（约8 km），这些数据都无法精确评估退耕还林前后黄土高原植被的动态演变。

谷歌地球引擎（Google Earth Engine）平台是基于云计算技术用于处理和分析大型地理空间数据集的大规模计算设施，为海量遥感数据的快速处理分析提供了前所未有的发展机遇（宁晓刚等，2022）。目前GEE已被广泛应用在大尺度的植被监测、土地利用分类以及碳循环等诸多领域。虽然也有研究利用 GEE平台对黄土高原的植被覆盖时空演变及其影响因素进行了分析（郭永强等，2019；常铮，2022），但较少考虑植被类型、地形因子等因素对FVC的影响，对退耕还林（草）前后其FVC的动态变化分析略显不足。

基于此，借助GEE平台，利用1986—2021年30 m空间分辨率的Landsat地表反射率数据计算长时间序列黄土高原NDVI和FVC，分析退耕还林（草）工程实施前后（1986—1999、2000—2021和1986—2021年）FVC的时空演变趋势，揭示地形因子对植被的影响。在此基础上，采用偏相关分析和残差分析等方法，定量探讨气候因子和人类活动对FVC的影响。研究结果可为黄土高原后续生态工程建设方案的制定提供科学依据和数据支撑。

1 研究区概况

黄土高原地处黄河中游（33°43′—41°16′N，100°54′—114°33′E），包含山西省、宁夏回族自治区的全部地区以及陕西省、甘肃省、内蒙古自治区、青海省和河南省的部分地区，总面积约 6.48×105km2（赵安周等，2017；张家政等，2022）。该地区地势由西北向东南倾斜，海拔差异超过3000 m（图1a）。气候类型属温带大陆性季风气候，冬季寒冷干旱，夏季炎热多雨，年平均气温 3.6—14.3 ℃，多年平均降水150—800 mm（易浪等，2014）。植被类型以草地、栽培植被和林地为主，其面积占研究区总面积的85%以上，中部区域为退耕还林（草）的重点区域（图 1b）（Zheng，2019；Zhao et al.，2022）。受地形地貌、气候变化以及不合理的人类活动，黄土高原是世界上水土流失最严重的地区之一（郭永强等，2019）。

图1 研究区位置和植被类型图Figure 1 Location and vegetation type of Loess Plateau

2 数据和方法

2.1 数据来源与预处理

2.1.1 遥感影像数据

1986—2021年 Landsat（5，7和 8）地表反射率数据集（Surface Reflectance，SR）来自于美国地质调查局，由 GEE平台获取（https://code.earthengine.google.com/），该数据产品已经过大气校正，可以消除大气散射、吸收以及反射引起的误差，其时间分辨率为16 d，空间分辨率为30 m×30 m。其中 1986—2011、2012、2013—2021年分别使用Landsat 5（19301景）、Landsat 7（786景）、Landsat 8（8800景）数据，黄土高原共计使用28887景数据（图2）。

图2 1986—2021年Landsat 5/7/8用于年合成采用的影像数量Figure 2 Number of scenes of Landsat 5/7/8 used for annual composite formation in 1986-2021

借助 GEE平台对黄土高原范围内遥感影像的QA质量波段做按位与运算，根据需求设定掩膜值去除受到云、阴影和雪等较大影响的像元以实现对像元的质量控制，进而达到去云的效果。在此基础上计算其NDVI，同时为进一步消除大气及传感器角度的影响，采用最大值合成法合成年 NDVI数据，获取1986—2021年黄土高原年最大NDVI数据集。

2.1.2 气象、植被类型和地形数据

1986—2021年黄土高原年平均气温和降水数据来源于中国科学院资源环境科学数据中心（http://www.resdc.cn），该数据是基于全国2400多个气象站点观测数据，利用 ANUSPLIN插值软件插值而成，空间分辨率为 1 km×1 km。太阳辐射（Solar Radiation，RAD）年合成数据来源于 ERA5再分析数据集（https://cds.climate.copernicus.eu/），其时间跨度为 1981—2021年，空间分辨率为0.1°×0.1°。

植被类型来自于中国科学院等单位编制的 1∶100万中国植被类型图（侯学煜，2001），依据该类型图将黄土高原分为栽培植被、林地、草地、灌丛、荒漠以及其他6种类型（图1b）。在计算FVC与降水、气温和RAD的偏相关系数以及残差分析前，利用GEE平台将FVC和RAD数据重采样为1 km×1 km。

数字高程模型（Digital Elevation Model，DEM）数据来自地理空间数据云平台（https://www.gscloud.cn/），该数据空间分辨率为30 m，利用该数据对黄土高原的坡向、高程、坡度3种地形因子进行提取。通过使用等间距分类法（李晶等，2021）对上述得到的地形因子进行重分类，统计不同地形因子的FVC等级占比及其均值变化。

黄土高原矢量边界数据来自于中国科学院资源环境科学数据中心（https://www.resdc.cn/），退耕还林重点区域的界线来源于Zheng et al.（2019）和Zhao et al.（2022）。

2.2 研究方法

2.2.1 像元二分模型

FVC的估算采用像元二分模型，该方法假定植被区的混合像元仅由植被和土壤构成，植被覆盖部分所占像元面积比例即为该像元的 FVC值（李苗苗等，2004；徐勇等，2022）。其计算公式如下：

式中：

FVC——植被覆盖度，其值介于0—1之间；

NDVIsoil和 NDVIveg——分别表示无植被覆盖和纯植被覆盖时其地表的NDVI值，其中，裸地像元值NDVIsoil和纯植被覆盖像元值NDVIveg的理论值应分别接近0和1；

RNI和 R——分别表示近红外和红光波段的反射率。参考已有的文献如李苗苗等（2004）和徐勇等（2022），本文取当年NDVI累计频率的5%和95%分别作为当年的NDVIsoil和NDVIveg。

2.2.2 趋势分析

采用非参数化趋势度（Sen）方法从像元尺度上计算黄土高原 FVC的变化趋势，该方法是一种稳健的非参数统计的趋势计算方法，计算效率高，且对测量误差和离群数据不敏感（Yu et al.，2020；Zhao et al.，2022），其具体计算公式如下（解晗等，2022）：

式中：

β——FVC的时间变化趋势，β＞0表示FVC呈现上升趋势，β＜0表示FVC呈现下降的趋势；

a和b——时间序数；

xa和xb——第a和b时间的FVC值。为进一步分析FVC变化的显著性，采用Mann-Kendall检验方法对其变化的显著性进行检验（徐建华，2014）。

2.2.3 偏相关分析

为探究气温、降水和RAD对植被生长的影响，对1986—2021年黄土高原FVC值与年平均气温、年降水和年RAD进行逐像元偏相关分析（Evans et al.，2004），其公式如下：

式中：

rxyz——变量 z固定后变量 x和 y的偏相关系数，即xy相关中剔除z的影响；

rxy——变量x与变量y的相关系数；

rxz——变量x与变量z的相关系数；

ryz——变量y与变量z的相关系数。

2.2.4 残差分析

采用残差分析区分气候变化和人类活动对FVC的影响，该方法是由Evans et al.（2004）提出。通过对降水、温度和 RAD进行逐像元多元线性回归分析，拟合得到FVC的预测值，将其视为气候因素对FVC的影响，进而对利用FVC真实值与预测值之差来区分气候变化和人类活动对 FVC变化的影响（李晶等，2021；杨灿等，2021；常铮等，2022），其公式如下：

式中：

x1、x2和x3——年平均气温、年降水和RAD 3个气候因子；

a、b和c——FVC与气候因子的回归系数；

d——回归常数项；

ε——残差值，ε的值域为[-1—1]，ε＞0和ε＜0的时候分别表示人类活动对FVC起促进和抑制作用，ε=0表示人类活动对FVC无影响；

FVCCC和FVCobs——FVC模拟值和实际值；

Trend(ε)——人类活动影响下FVC变化趋势；

Trend(FVCCC)——实际FVC变化趋势；

H——人类活动的贡献率。

3 结果分析

3.1 黄土高原FVC年际变化特征

1986—2021年黄土高原FVC呈波动上升的趋势，其值由 1986年的 0.3901上升到 2021年的0.5421，增速为 0.0044 a-1（P ＜0.01，图 3）。分时间段看，1986—1999年期间，其FVC值由0.3901增至 0.4369，增速为 0.0038 a-1（P＜0.01，图 3）；2000—2021期间，FVC值由0.4612增至0.5421，增速为0.0058 a-1（P＜0.01），增速明显快于前一阶段（图 3）。从不同的植被类型看（表 1），1986—2021年期间黄土高原各个类型植被FVC均呈显著上升趋势（P＜0.01），其中草地的上升速率最快（Trend=0.0055 a-1）。从不同时间段看，2000—2021年的耕地、林地、草地和灌丛的上升速率均高于1986—1999年。

表1 1986—1999、2000—2021和1986—2021年黄土高原不同植被类型的FVC变化趋势Table 1 Trends of FVC changes between different vegetation types in the Loess Plateau in 1986-2021, 1986-1999 and 2000-2021

图3 1986—1999、2000—2021和1986—2021年黄土高原FVC年际变化Figure3 The variation of average annual vegetation cover of the Loess Plateau in 1986-2021,1986-1999 and 2000-2021

3.2 黄土高原FVC空间格局演变

3.2.1 空间分布特征

将FVC划分为低覆盖度[0, 0.20]、中低覆盖度（0.20, 0.40]、中等覆盖度（0.40, 0.60]、中高覆盖度（0.60, 0.80]和高覆盖度（0.80, 1.00] 5个等级，并提取FVC值为0.60的分界线（李辉霞等，2011）。图4显示了 1986—1999、2000—2021和 1986—2021年黄土高原 FVC多年均值的空间分布。总体上，1986—2021年黄土高原FVC均值呈西北向东南递增的分布格局，低值区主要分布在北部的内蒙古自治区、宁夏回族自治区和甘肃省等地，这些地区多为草地和荒漠沙地；高值区域主要分布在东南部的山西和陕西省，这些地区的植被类型以灌溉农业、灌丛和林地为主（图4a）。与1986—1999年相比，2000—2021年FVC=0.60曲线呈向西北方向移动的趋势，表明生态工程建设以后，其FVC＞0.60的像元呈增多的态势（图4b和4c）。从1986—2021年的像元分布频率看，FVC处于中等覆盖度的像元占比最大（23.61%），其中1986—1999年黄土高原地区FVC以中等及以下覆盖度占比达67.20%，2000—2021年 FVC以中等及以下覆盖度占比减少至57.36%，表明生态工程建设以后黄土高原植被得到恢复（图4d）。

图4 1986—1999、2000—2021和1986—2021年黄土高原FVC均值空间分布及不同等级FVC像元占比Figure 4 Spatial distribution and the frequencies of FVC on the Loess Plateau in 1986-2021,1986-1999 and 2000-2021

3.2.2 空间趋势特征

图5显示了1986—2021、1986—1999和2000—2021年黄土高原像元尺度FVC的变化趋势。从图中可以看出，黄土高原FVC整体趋势以增加为主。1986—2021年 FVC呈增加和减小的面积分别占67.19%和32.81%，其中显著上升的面积占53.66%，主要位于黄土高原的中北部；显著退化的面积仅占8.68%，主要位于南部的西安等城市的周边以及西北部的内蒙古自治区和宁夏回族自治区等地（图5a和5d）。从不同时间段看，1986—1999年FVC值呈显著上升和显著减少的分别占18.38%和6.43%（图5b和5d）；到2000—2021年，FVC值呈显著上升的比例上升至48.12%，主要位于黄土高原中部地区（图5c和5d）。

图5 1986—1999、2000—2021 和1986—2021 年黄土高原FVC 空间趋势Figure 5 Spatial distribution of FVC trends on the Loess Plateau in 1986-2021, 1986-1999 and 2000-2021

3.3 FVC影响因素分析

3.3.1 FVC的地形效应

图6显示了1986—2021年黄土高原FVC随地形因子的变化情况。从图中可以看出，FVC值随高程的变化呈“下降—上升—下降”的变化趋势，最高值（0.7793）出现在 3.0—3.5 km，最低值（0.4285）出现在 1.0—1.5 km。FVC各等级面积占比随高程不同差异明显，其中低等级FVC和中低等级FVC面积占比随高程变化起伏较大，当海拔超过 4 km时，中高等级FVC和高等级FVC占比仅为2.41%，而中低等级FVC及低等级FVC的比例达到75.01%（图6a），主要原因是由于该海拔高度已远超林线，随着海拔升高气温和空气密度都会降低，不利于植被生长，使得植被覆盖情况较差（李辉霞等，2011；王晓蕾等，2022）。从不同坡向FVC的变化看，其均值在0.4571—0.4828之间，各等级FVC占比差异不大，表明坡向对黄土高原 FVC影响效果并不明显（图6b）。从不同的坡度等级看，FVC值随坡度的增加呈增加的趋势，其值从0°—5°的0.4282上升到 25°—45°的 0.7025（图 6c），具体而言，15°—25°中高等级 FVC和高等级 FVC比例最大（73.93%）；0°—5°的低等级FVC和中低等级FVC比例最大（42.82%），主要是由于城市等人口密集区域主要集中在坡度较小的地方，受城市扩张等人类活动的干扰，其低等级 FVC占比较高（马士彬等，2019）。

图6 黄土高原不同地形因子的FVC等级占比及其均值变化Figure 6 Proportion and FVC means of vegetation cover at different levels of terrain factors

3.3.2 气候变化对FVC的响应

除地形因子外，植被还会受到气候因子的影响（Feng et al.，2016；Liu et al.，2020；Gao et al.，2022；肖强等，2016；聂桐等，2022）。1986—2021年黄土高原FVC与年降水量、年平均气温和RAD的偏相关系数分别为0.239、0.093和-0.006，其中呈显著正相关（P＜0.05）像元比例分别为48.50%、22.51%和5.96%（图7）。FVC与年降水的偏相关系数整体上大于年均气温和RAD，表明降水是影响黄土高原植被变化的主要气候因素（图 7a）。黄土高原地处干旱半干旱区域，降水的增多会提高土壤含水量，促进植被生长（郭永强等，2019；张家政等，2022）。FVC与气温呈显著负相关的比例为6.15%，主要位于西北部的甘肃省、宁夏回族自治区以及内蒙古自治区（图7b）。此外，FVC和太阳辐射呈显著负相关的比例为7.81%，主要位于西北部的青海省、宁夏回族自治区、内蒙古自治区等地（图 7c），这些地区的太阳辐射过强，过度红外或紫外辐射不利于植被的生长（何亮，2021；孙高鹏等，2021）。

图7 1986—2021 年FVC 与气候因子偏相关系数的空间分布 Figure 7 Spatial distribution of partial correlation coefficients between FVC and climate factors in 1986-2021

3.3.3 人类活动对于FVC的影响

人类活动也是该区 FVC变化的重要驱动因素之一（Kou et al.，2021；张宝庆等，2021）。残差分析结果表明，人类活动是 FVC变化的主要贡献因子（80.92%），并将黄土高原FVC的残差趋势值分为显著负向影响（Trend＜0，P＜0.05）、轻微负向影响（Trend＜0，P＞0.05）、轻微正向影响（Trend＞0，P＞0.05）和显著正向影响（Trend＞0，P＜0.05）4 个等级。从图中可以看出，人类活动对黄土高原植被生长起正向和负向作用的像元比例分别为 73.20%和 26.80%，其中呈显著正向影响像元面积占比为26.99%，主要位于中部的陕西、甘肃等退耕还林（草）工程重点区域（图 8）。残差趋势为负值的区域主要分布在黄土高原东南部的西安、太原以及西北部的银川等省会城市附近，这些地区大多经济相对发达，人口密集，城市扩张较快，对植被生长起消极作用（图8）。总体上，人类活动对该地区植被呈积极影响，但应关注城市扩张等人类活动对植被的消极影响。

图8 1986—2021年黄土高原人类活动对FVC影响的空间分布Figure 8 Spatial distribution of human activities impacts on FVC in the Loess Plateau in 1986-2021

4 讨论

目前的研究表明近几十年黄土高原的植被呈改善的趋势，如Kou et al.（2021）利用SPOT NDVI和GIMMS NDVI分析了黄土高原1998—2018年黄土高原 FVC的时空演变及影响因素，结果表明近80%的区域呈变绿的趋势；Zhao et al.（2017）利用GIMMS NDVI数据分析了1982—2013年黄土高原生长季NDVI的时空演变，认为90%以上的区域呈变绿的趋势；李依璇等（2021）采用MODIS数据分析了 2000—2018年黄土高原植被动态变化，认为70%以上的区域被覆盖度呈增加趋势。上述研究存在差异的原因可能是数据源以及研究时间尺度的不同造成的。本研究结果表明，1986—2021年黄土高原FVC呈增加的面积占67.19%，2000—2021年黄土高原FVC的增加速率明显快于1986—1999年（图2和图5），尤其在生态工程重点建设区尤为显著（图5a），表明退耕还林（草）等生态工程建设对植被有积极影响，这与赵安周等（2017）、郭永强等（2019）的研究结果类似。本文采用GEE平台的Landsat地表反射率数据计算的FVC同时兼顾了高空间分辨率和长时间序列，不仅弥补了 MODIS、SPOT等数据时间跨度较短的不足，而且避免了GIMMS数据空间分辨率较低的问题，因此可以更加清晰地表征黄土高原，尤其是城市扩张区域植被的时空演变（图5）。黄土高原多为半干旱地区，水资源匮乏，但是水分作为光合作用的主要原料之一，降水的增加会促进植被的快速生长。1986—2021年黄土高原降水呈显著增加的趋势（Trend=1.975 mm·a-1，P=0.015，图 9），黄土高原FVC与降水呈正相关的面积比例达到 79.73%（图7a），进一步表明降雨的增多会促进该地区植被的生长。同时，由于生态系统的反馈机制，植被的快速增加会影响下垫面条件，进而会对陆面与大气间的交互作用机制产生影响，进一步增加局部降水（Gao et al.，2022；张宝庆等，2021）。此外黄土高原东南部的关中平原及北部的河套平原 FVC值较高，主要是由于这些地方主要为灌溉农业区，人为灌溉会对植被生长产生积极作用（Kou et al.，2021；金凯等，2020）。本研究结果也表明，在西安、太原、银川等城市周边的植被出现显著退化的趋势，主要是由于快速的城市扩张将草地、农田以及林地等转变成城市不透水表面，使得植被减少（Kou et al.，2021）。不同的地形条件下其FVC值不同主要是由于不同的海拔和坡度其土壤所含水分、热量条件以及矿物质含量与种类的不同，使得植被生长存在差异，进而产生空间异质性（王晓蕾等，2022）。

图9 1986—2021、1986—1999和2000—2021年黄土高原年降水量及趋势Figure 9 Trends of precipitation in the Loess Plateau in 1986-2021, 1986-1999 and 2000-2021

考虑到本地计算机的性能难以满足大空间尺度的数据处理，本研究利用 GEE平台对黄土高原的Landsat数据进行了处理得到了1986—2021年黄土高原逐年的 FVC数据，该平台可以在线高效处理大范围长时间序列的遥感数据，快速实现影像统计、趋势分析等方面的研究，对大区域环境下的植被变化进行实时监测，极大的缩短影像处理时间，提高工作效率。应该指出的是，虽然本文利用30 m空间分辨率的Landsat数据对1986—2021年黄土高原 FVC的时空演变格局及影响因素进行了分析，但仍存在一定的不确定性。首先，考虑到Landsat数据的可获取性，本研究仅对年尺度的 FVC变化进行了分析，并未分析其年内变化特征，已有研究表明植被对气候因子的响应具有季节性（阎世杰等，2019）。其次，FVC的变化对多种自然和人为要素具有较强的依赖性和敏感性，其变化会受到多因子间的交互作用，本研究虽采用残差法剥离了气候干扰，但人类活动影响因子多样，如何选择人类活动因子更详细地划分不同类型的人类活动（农业灌溉和城市扩张等）对FVC的影响是今后的研究重点。最后，虽然退耕还林（草）等生态工程提高了黄土高原的FVC，有效抑制入黄的泥沙量，但另一方面植被覆盖度的提高会增加该地区水资源的压力，进而带来土壤干化等问题（Wang et al.，2021），因此，如何因地制宜选择合适的植物进行植树造林值得进一步研究。

5 结论

（1）时间上，1986—2021年黄土高原FVC均值整体上呈波动上升趋势，2000—2021年的增加趋势快于1986—1999年。空间上，1986—2021年黄土高原多年 FVC均值总体呈从西北向东南递增的分布格局，其中呈显著上升面积主要集中在陕西省北部、山西省西部等地。

（2）对FVC变化进行地形效应响应分析可知，坡度和高程对FVC的影响较显著。FVC值随高程呈“下降—上升—下降”的变化趋势；FVC随坡度的增加呈增加的趋势。

（3）在黄土高原地区人类活动对 FVC变化的贡献程度高于气候的影响。在人类活动方面，73.20%的区域对 FVC变化起积极作用；在气候因子方面，年降水量对FVC变化影响最大。

GEE平台具有海量共享的遥感影像和地理数据，依托谷歌云平台和强大的后台处理器计算可处理大范围、长时间序列的遥感影像数据，极大的缩短影像处理时间，提高工作效率。因此本研究GEE平台的30 m空间分辨率的Landsat地表反射率数据计算了黄土高原地区 30 m分辨率、长时间序列（1986—2021年）的FVC数据，可以从更精细尺度分析黄土高原地区的植被覆盖变化及对气候变化和人类活动的响应，其研究结果可为黄土高原生态环境保护及未来退耕还林（草）政策的制定提供决策和建议。