基于GEE云平台的榆林市植被覆盖度时空变化研究

白欣，贺炳彦，高婷，李雪，孙乔

（1.长安大学地质工程与测绘学院，西安 710000；2.中煤航测遥感集团有限公司，西安 710000；3.西安灵图信息科技有限公司，西安 710000；4.航天建筑设计研究院有限公司，北京 100162）

植被作为生态大循环的重要参与者，在全球能量流动和物质循环中起着重要作用［1］。通常用植被覆盖度（Fractional vegetation cover，FVC）表征地表植被覆盖情况，一般将植被覆盖度定义为植被冠层在地面上垂直投影面积与土地总面积的比值［2］。国内外学者常通过像元二分模型来估算植被覆盖度，多使用归一化植被指数来分解像元遥感信息，从而监测研究区的植被覆盖变化。张舒婷等［3］采用Landsat数据，利用像元二分法估算植被覆盖度，并监测了黄土高原丘陵沟壑区18年的植被覆盖变化；张志强等［4］采用像元二分模型、Hurst指数分析了2000—2019年黄河流域的植被覆盖度时空变化特征和发展趋势。

榆林市位于黄土高原腹地，水土流失严重，是中国典型的生态脆弱区。针对黄土高原风蚀沙化严重、土壤盐碱化等问题，从1999年开始，国家实施退耕还林还草措施，以期植被较快恢复，国内许多学者也开始关注榆林市退耕还林还草后的植被覆盖状况。杨波等［5］使用MODIS数据研究了榆林市自还林还草工程以来的植被变化和未来发展趋势，发现18年来榆林市归一化植被指数（Normalized difference vegetation index，NDVI）共增长了88.56%，未来有10%的增长潜力；马锋等［6］研究发现榆林市2000—2018年植被覆盖度呈波动上升的趋势，且南部黄土丘陵沟壑区的植被覆盖度增速高于北部风沙草滩区。但少有学者使用高分辨率数据对榆林植被覆盖进行长时序的监测，主要原因是植被覆盖度计算过程中高分辨率数据存储不便、计算量大、运行时间长等问题。

Google earth engine（GEE）平台的出现可以有效解决上述问题。GEE是一个基于云计算的“遥感大数据”分析和呈现平台，相较于单机地理信息处理平台而言，GEE不仅有大量的影像数据可供用户使用，还提供了高性能的集群服务器，实现海量影像数据在线可视化处理，大大提高了工作效率［7］。因此，本文基于GEE平台，针对榆林市的生态环境问题和还林还草工程，采用Landsat surface reflectance（陆地卫星地表反射率数据）和像元二分模型反演了2000—2019年榆林市植被覆盖度，分析植被覆盖时空变化特征，同时分析了气候因子（年平均气温、年降水量）和地形因子（高程、坡度）对植被覆盖变化的影响，最后通过Hurst指数预测了未来研究区植被覆盖的变化趋势，为评估榆林市还林还草工程提供科学依据。

1 研究区与数据

1.1 研究区概况

榆林市地处陕西省最北部、陕北黄土高原与毛乌素沙地南缘交界处，东隔黄河与山西相望，南与陕西省延安市接壤，西邻宁夏、甘肃，北连内蒙古，高程530～1 922 m（图1），面积约43 578 km2。研究区属温带大陆性季风气候，年均温度为10℃，年平均降水量为400 mm左右，气象灾害较多，立地条件较差［8］。以长城为界，北边为风沙草滩区，南部为黄土丘陵沟壑区。土壤多以风沙土、黄绵土为主，水土流失和风蚀沙化严重，是退耕还林还草工程的典型区域。

图1 榆林市位置

1.2 数据来源及预处理

本研究采用的遥感影像数据为Landsat surface reflectance，空间分辨率30 m、时间分辨率16 d，来源于GEE遥感计算云平台。首先，在GEE平台上调用“LANDSAT/LT05/C01/T1_SR”“LANDSAT/LE07/C01/T1_SR”“LANDSAT/LC08/C01/T1_SR”影像数据集，该产品已经进行了大气校正［9，10］。利用Landsat质量评估波段（QA）检测了像元的云、云影、雪等干扰因素，筛选出了3 042景低云量影像。依据GEE平台上的CFMASK算法去云处理后，在线计算归一化植被指数，采用最大值合成方法合成年度NDVI最大值影像，在此基础上选用像元二分模型估算植被覆盖度，最后采用qualityMosaic算法在线裁剪影像。

2000—2019年榆林市年平均气温和年降水量数据均来源于国家地球系统科学数据中心（http：//www.geodata.cn），空间分辨率为1 km。将该数据上传至GEE平台，用上述qualityMosaic算法在线裁剪，得到榆林市每年的年平均气温和年降水量分布数据。从地理空间数据云上获取坡度和高程数据，空间分辨率为30 m。

2 研究方法

2.1 植被覆盖度估算

植被指数通常用归一化植被指数值来表示，用来衡量植被生长状况和空间分布，是最常用的一种植被指数，其计算方法是近红外波段的反射值（NIR）和红外波段反射值（RED）之差与其之和的比值，计算公式如下：

采用像元二分模型估算植被覆盖度，计算公式如下：

式中，NDVIsoil和NDVIveg分别为纯裸土和纯植被像元的NDVI值，它们的理论值应为0和1，但受地表湿度、土壤类型、土壤颜色等的影响，其实际值会随空间和时间而变化［11］。因此，在年最大NDVI累积频率表上取5%和95%作为NDVIsoil和NDVIveg，然后代入式（2）中，计算榆林市2000—2019年植被覆盖度［12］。根据水利部2008年颁布施行的《土壤侵蚀分类分级标准》［13］，本研究将植被覆盖度分为5个等级：低植被覆盖度（0≤FVC＜30%）、中低植被覆盖度（30%≤FVC＜45%）、中 植 被 覆 盖 度（45%≤FVC＜60%）、中高植被覆盖度（60%≤FVC＜75%）、高植被覆盖度（75%≤FVC≤100%）。

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

为了分析研究区植被覆盖度在过去20年的变化趋势及其时空特征，采用Sen+Mann-Kendall法分析FVC长时间序列变化趋势。Sen斜率估计用于计算趋势值，其计算公式如下［14］：

式中，xi、xj分别代表第i、第j年的植被覆盖度；β是斜率，用于判断时间序列的上升或下降，当β小于0时表示其呈现下降趋势，相反表示上升。

Mann-Kendall属于非参数检验方法，与其他参数检验的方法相比，样本受异常值干扰小，不需要遵从一定的分布，更适合顺序变量［15］。其检验统计量S的计算方法如下［16］

本文所用时间序列长度n值为20，故趋势检验使用检验统计量Z，计算方法见式（6），取显著水平α=0.05，Z1-α/2=Z0.975=1.96。

2.3 相关性分析

偏相关分析是当两个变量与第三个变量相关时，控制第三个变量不变，只分析剩余两个变量之间的相关性。选取年平均气温和年降水量作为植被覆盖度的相关变量，分别计算其与植被覆盖度的偏相关系数，公式如下［18］：

式中，rxy为变量x、y的相关系数，xi和yi是第i年的值，xˉ和yˉ分别是两个变量的平均值。rxy*z为控制自变量Z不变时因变量x与自变量y的偏相关系数。显著性检验采用t检验法，统计量t的计算公式如下：

式中，rxy*z和rxz分别为偏相关系数和相关系数，n为年数，m为自变量个数。

2.4 Hurst指数

Hurst指数是描述自相似性和长期依赖性的有效方法，并在水文、气候、地质等领域广泛运用［19］。本文采用Hurst指数来描述研究区FVC未来的变化趋势，其具体计算过程参考严恩萍等［20］的方法。Hurst指数取值包括3种形式：当0.5≤Hurst时，表示时间序列具有长期持续性，未来的变化趋势与过去相同；当0≤Hrust＜0.5时，说明时间序列具有反持续性，未来的变化趋势与过去相反；当Hrust=0.5，表明时间序列是相互独立的，无相关性。

3 结果与分析

3.1 植被覆盖度空间分布特征

为了解榆林市植被覆盖度空间分布情况，利用2000—2019年的最大植被覆盖度，合成了榆林市20年平均植被覆盖度分布图，如图2a所示。榆林市植被覆盖度由东南向西北逐渐递减，这与榆林市北部风沙草滩区、南部黄土高原丘陵沟壑区的地貌有关。低植被覆盖度区域和中低植被覆盖度区域主要分布在北部风沙草滩区，该区域地处毛乌素沙地南缘地带，沙化风蚀严重，植被恢复困难。中植被覆盖度区域、中高植被覆盖度区域和高植被覆盖度区域主要分布在南部黄土丘陵沟壑区，该地区空气湿度、土壤肥力等自然环境较风沙区更好，加之退耕还林工程的实施，使得植被覆盖更高。统计各等级植被覆盖度所占面积比例，见图2b。榆林市中植被覆盖度所占面积比最大，占总面积的25.1%。

图2 2000—2019年榆林市平均FVC空间分布

3.2 植被覆盖度动态变化特征

图3表示研究期内榆林市植被覆盖度最大值的年际变化。榆林市植被覆盖度整体呈波动增长的趋势，2019年最大植被覆盖度为52.3%，较2000年增长了11.2个百分点，年增长率为1.37%（P＜0.05）。最大植被覆盖度在2001年为最低值，可能因为2001年退耕还林还草工程刚刚开始，耕地面积减少，而新种植的植被覆盖较低［21］。

图3 2000—2019年榆林市最大FVC年际变化

基于SEN+Mann-Kendall的植被覆盖度趋势分析结果表明（图4），榆林市8.3%的地区植被覆盖度呈严重退化趋势，19.6%的地区植被覆盖度呈轻微退化趋势，9.1%的地区植被覆盖度稳定不变，30.7%的地区植被覆盖度呈轻微改善趋势，32.3%的地区植被覆盖度有明显改善。植被覆盖度呈改善的地区集中在黄土丘陵沟壑区的子州县、吴堡县、绥德县、米脂县、清涧县，横山区、相比于风沙草滩区，这些地区降水量更多，土壤条件更好，植被恢复最为明显。呈退化趋势的地区主要分布在定边南部和靖边南部、榆阳区西部和神木西部，属于风沙区，说明这些地区沙地生态系统相对脆弱，植被生存条件不稳定［22］。

图4 2000—2019年榆林市植被覆盖度变化趋势分布

3.3 植被覆盖度变化与高程、坡度的分异特征

3.3.1 植被覆盖度变化与高程的分异特征由于不同高程地区的气温、空气湿度、土壤肥力等因素不同，因此对植被生长的影响程度也不同。从图5a可以看出，榆林市平均植被覆盖度随高程增加呈先上升后下降的趋势，平均植被覆盖度的最大值出现在高程730～930 m的区域内，值为0.62，这里位于各川川谷，地势平缓，土壤水分良好，适宜植被生长和植树种草；1 130～1 330 m高程区域的平均植被覆盖度最低，值为0.41。高程730～1 130 m范围内的植被覆盖度改善最为显著，之后随着高程增加，植被覆盖改善效果逐渐下降，植被覆盖退化渐为明显。1 530～1 922 m高程范围内的植被覆盖退化区域面积占比最大，该区域沟壑纵横，生态脆弱，水土流失严重，人工植树造林难度大，因此植被退化较为明显，可见大于1 530 m高程区域内的植被退化需引起进一步关注。3.3.2植被覆盖度变化与坡度的分异特征由图5b可知，榆林市植被覆盖度随坡度的增加呈缓慢上升的趋势，平均植被覆盖度在坡度为0°～5°的区域最低，为0.39，在坡度大于35°的区域最高，为0.59。一方面，坡度平缓的区域人类活动较为密集，坡度大的区域人类活动受限，植被生长受到的干扰减小，一定程度上导致植被覆盖随坡度增大而上升［23］；另一方面，退耕还林工程按照“先陡坡后缓坡”的原则实施，对于坡度大、水土流失严重的坡耕地优先造林，因此坡度大的地区植被覆盖度越大，植被覆盖度的改善面积所占比重也随之增大。坡度为0°～5°的区域内植被覆盖退化面积所占的比重最小，为31.8%，而改善的面积所占比重最大，为53.3%，可能是因为坡度平缓的地区植被覆盖受人为影响干扰最大。

图5 不同地形下植被覆盖度及时序变化趋势

3.4 植被覆盖度与气候因子相关分析

将榆林市2000—2019年植被覆盖度与年平均气温和年降水量进行偏相关分析和显著性检验，见图6。可以看出，植被覆盖度与降水的偏相关性更大，降水量更能决定植被的生长状况。植被覆盖度与年降水量呈显著相关（P＜0.05）的像元占榆林市总像元数的11.5%，其中，呈显著正相关的区域占9.0%，主要分布在中部的黄土丘陵沟壑区，该区域水土流失严重，植被生长对降水的依赖性较强。植被覆盖度与年平均气温的偏相关关系中，有7.2%的区域通过了显著性0.05的检验，呈显著正相关（P＜0.05）的像元占榆林市总像元数的5.4%，主要分布在府谷县、神木市、佳县等东部地区；仅有1.8%的地区植被覆盖度与年均气温呈显著负相关，主要位于定边县南部、靖边县西南部，这些地方海拔较高，降水量较少，温度升高导致植被表层土壤含水量减少，进而加剧植被根层干燥化，抑制植被的生长［24］。

图6 植被覆盖度与降水（a）和气温（b）的偏相关性

3.5 植被覆盖度变化的持续性分析

利用2000—2019年榆林市的植被覆盖度数据，计算Hurst指数，以预测植被覆盖度的未来变化趋势。榆林市Hurst指数均值为0.58，Hurst指数小于0.5的面积占19.1%，大于0.5的面积占80.9%，表明未来一段时期内榆林市植被覆盖度变化趋势与2000—2019年的变化趋势保持一致的地区面积占比为80.9%。结合前面的FVC变化趋势得到榆林市植被覆盖度变化的持续性特征空间分布图（图7），一共得到6种变化情形：持续性严重退化、持续性轻微退化、持续性稳定不变、持续性轻微改善、持续性明显改善、未来变化趋势不确定。

图7 基于Hurst指数的FVC变化趋势分布

植被覆盖度持续性轻微改善和明显改善的区域面积分别占21.9%、31.9%，主要分布在府谷县、佳县、米脂县、吴堡县、绥德县、清涧县、子州县、横山区、神木市东南部和榆阳区东南部；持续性稳定不变的区域占9.1%，主要分布在毛乌素沙地南缘等风沙草滩区和靖边县最南端；持续性轻微退化和持续性严重退化区域面积分别占14.0%、8.0%，主要分布在定边县南部、靖边县南部以及神木市西北部；未来变化趋势不确定的地区占15.1%，这部分地区需研究人员进一步关注和研究。同时，持续性退化的区域还需要研究人员重点关注其植被恢复受干扰的原因。

4 结论

1）榆林市植被覆盖度由西北向东南递增，北部风沙草滩区的植被覆盖度低，南部黄土丘陵沟壑区的植被覆盖度高；中植被覆盖度所占面积最大，约占榆林市总面积的25.1%。

2）2000—2019年榆林市植被覆盖度整体呈波动增长的趋势，20年共增长了11.2个百分点，年增长率为1.37%；研究区内有63.0%的区域植被覆盖得到了改善，其主要分布在土壤肥力、降水条件更好的黄土丘陵沟壑区；有27.9%的区域植被覆盖呈退化趋势。

3）榆林市植被覆盖度受地形因子和气候因子的影响，其高程、坡度分异特征明显，植被覆盖度随高程增加呈先上升后下降的趋势，随坡度的增加呈缓慢上升的趋势；在高程为730～1 130 m时植被改善区域面积占比最大，在坡度为0°～5°时植被退化区域面积占比最大；气候因子中，降水对榆林市植被覆盖的影响大于气温，植被覆盖度与降水呈显著相关的像元占比为11.5%，与气温呈显著相关的像元占比为5.4%。

4）榆林市Hurst指数均值为0.58，有80.9%的植被覆盖区域具有正向持续性，其中，持续改善的区域占53.8%，持续性稳定不变的区域占9.1%，持续退化的区域占22.0%，还有15.1%的区域未来变化趋势不确定，需研究人员继续关注。