基于GEE的诺木洪洪积扇植被时空变化特征、成因及趋势分析

姚金玺， 张志， 张焜

(1.中国地质大学(武汉)地球物理与空间信息学院，武汉 430074； 2.青海省青藏高原北部地质过程与矿产资源重点实验室，西宁 810300)

0 引言

植被是陆地生态系统中的重要组分。理解植被对气候变化的响应对认识与预测未来生态系统的演化具有重要意义[1]，是干旱、半干旱区生态环境监测的内容[2]。近几十年来，受到全球气候变暖的影响，青海省诺木洪地区气候向暖湿化发展，径流有所增加。尽管对诺木洪河的管理得到加强，但由于农业用水的增加，对下游生态用水、植被发育产生一定的影响。

利用时间序列植被指数监测陆地植被生长变化规律并探讨气候要素的驱动作用已成为植被-气候相互作用研究中的重要方向之一[3-4]。气候变化和人类活动可能导致植被出现新的变化，需要对其重新认识[5]。基于传统数据平台和低空间分辨率影像分析大区域植被变化和驱动因素已开展较多研究[6-9]。基于遥感和地理信息系统(geographic information system,GIS)技术,徐浩杰等[6]和王林林等[7]利用生长季MODIS数据分析植被季节时空变化及其驱动因子，发现柴达木盆地植被覆盖呈逐步改善趋势; 张斯琦等[8]和李艳丽等[9]利用MODIS数据分析植被覆盖度及人类活动与植被演化、径流改变之间的关系。但前人研究多利用低空间分辨率遥感数据监测大区域植被宏观变化，并采用传统方法下载、存储与处理遥感影像，研究成果精度与效率低。

Google Earth Engine(GEE)平台中数据包括近40 a的全球卫星影像，且提供超过800 种功能函数[10]。基于云计算的 GEE 平台可在线高效处理大规模地理数据集，有效地解决了遥感大数据中的处理难题[11-12]。

基于GEE平台，利用2000—2017年MODIS增强植被指数(enhanced vegetation index，EVI)数据和Landsat归一化差异植被指数(normal difference vegetation index，NDVI)数据，分析长时间序列植被变化特征，重点研究不同年代诺木洪洪积扇上枸杞种植区和盐碱化区植被变化特征以及气候对植被变化的影响。在此基础上分析人为因素在植被时序演化中所起到的作用，并利用重标极差分析方法研究诺木洪洪积扇植被变化趋势的可延续性，以期为诺木洪河水资源综合管理和可持续发展决策提供参考。

1 研究区概况及数据源

研究区主体位于青海省海西蒙古族藏族自治州都兰县诺木洪地区，地理坐标为N36°10′～36°35′，E96°15′～96°40′，海拔为2 767～3 191 m(图1)，属高原干旱大陆性气候，干燥少雨，年均气温为1.2～4.3 ℃，降水量为17.8～177.5 mm，集中在6—9月，表现出雨热同季[13]。区内主要植被类型为梭梭、柽柳和芦苇等植物，覆盖度不高，在晚期洪积扇扇缘及河流沿岸等含水量高的地区，植被长势良好。

图1 诺木洪洪积扇位置示意图

本文MODIS，Landsat5与Landsat8卫星数据均经过了辐射定标和大气校正，使用的EVI数据是250 m空间分辨率的MODIS卫星16 d合成的EVI产品MOD13Q1[14]，NDVI数据是基于Landsat5与Landsat8影像计算所得。通过GEE算法去云得到高质量的影像数据后还需进行空间和时间双维度的数据筛选，空间筛选利用clip或clipToCollection函数，时间跨度为2000—2017年。预处理过程和数值计算等操作在GEE平台上由JavaScript语言编程实现。

气象数据为诺木洪气象观测站的年均温度及总降水量，时间跨度为2000—2017年，其来源于中国气象科学数据共享服务网(http: //cdc.cma.gov.cn)。

2 研究方法

研究主要内容包括： 遥感影像的收集与预处理，获取降水量和温度数据； 计算植被覆盖特征因子和自然环境因子； 变化特征的分析及结果验证； 植被变化影响因素分析和未来变化趋势分析4个部分。技术路线如图2所示。

图2 技术路线

2.1 最大化合成与最大值计算法

对2000—2017年研究区Landsat卫星数据进行计算NDVI值，并筛选每个像素的NDVI最大值，予以合成，针对MODIS卫星EVI产品进行筛选最大的EVI值。利用2组数据以及不同的计算方法来表达整个区域当年植被生长最繁茂的时期，分析该时间段内月际、年际植被变化特征。

2.2 Pearson相关分析与一阶偏相关分析法

对温度、降水量和各个时间点的最大EVI值数据做相关分析与偏相关分析，由此解释温度、降水量与最大化EVI数据的相关程度[15]。2个变量之间Pearson相关系数ρXY定义为这2个变量的协方差与二者标准差积的商，即

(1)

式中：X和Y为2个变量；cov(X,Y)为X与Y的协方差；σXσY为X与Y的标准差之积；E(X)为X的数学期望。

当相关系数处于[0.8,1]时，二者高度相关； 处于[0.5,0.8)时，二者中度相关； 处于[0.3,0.5)时，二者低度相关； 小于0.3时则为弱度相关。

在3个变量中，任意2个变量的偏相关系数是在排除其余一个变量影响后计算得到的，称为一阶偏相关系数，公式为：

(2)

式中：rij·h为排除变量Xh影响后的Xi和Xj的偏相关系数；rij为变量Xi与Xj的简单相关系数；rih为变量Xi与Xh的简单相关系数；rjh为变量Xj与Xh的简单相关系数。

偏相关系数检验的零假设为： 总体中2个变量间的偏相关系数为0。使用t检验方法，公式为：

(3)

式中：r为相应的偏相关系数；n为样本观测数；k为可控制变量的数目；n-k-2为自由度。当t>t0.05(n-k-2)或p<0.05时，拒绝原假设。

2.3 Hurst指数

Hurst指数常用于定量描述时间序列变化趋势的可延续性[16]，可以更好地分析植被的年际变化特征。Hurst指数由英国水文学家 Hurst[17]提出，现已应用于地质、遥感和水文等领域中。

本文中利用重标极差分析方法计算Hurst指数，分析研究区域内未来短期的植被变化趋势。Hurst指数值H总体处于[0,1]区间之内，将此区间进行划分代表着不同的含义。当0≤H<0.5时，表明时间序列具有长期相关性，但将来的总体趋势和过去相反，即反持续[18]； 当H≥0.5时，表明时间序列具有长期相关的特征，即此过程具有持续性，且H越接近1，持续性越强。

3 结果与分析

3.1 植被时序变化特征分析

3.1.1 植被覆盖度

利用像元二分模型结合NDVI数据获取研究区2000年和2017年的植被覆盖度(图3)。参考相关植被覆盖度划分的文献[19],结合该研究区的气象数据、植被生长状态以及实地勘察情况，综合分析可知，在戈壁以及盐碱化严重的区域，往往只有少量梭梭树以及柽柳，覆盖度比较低，而在枸杞种植园地由于人为的作用，植被发育较好。因此，将植被覆盖度以0.5为节点，处于[0,0.3)，[0.3,0.5)，[0.5,0.8)和[0.8,1]分别划分为极低植被覆盖度、低等植被覆盖度、中等植被覆盖度与高等植被覆盖度。在整个区域18 a的时间内，极低植被覆盖度由66.04%上升到67.71%，低等植被覆盖度由12.84%上升到13.63%，中等植被覆盖度由12.33%下降到11.28%，高等植被覆盖度由8.79%下降到7.38%，总体上来看植被变化并不大，整体植被发育比较弱。

(a) 2000年(b) 2017年

3.1.2 植被指数

MODIS卫星和Landsat卫星的空间分辨率分别为250 m和30 m，时间分辨率分别为1 d和16 d。采用不同空间、时间分辨率的卫星产品，在不同空间、时间尺度上分析植被的变化特征，以对比分析使结果更具代表性。

利用最大化合成法得到的NDVI年均值和最大EVI年均值来表达每年植被最茂盛的时期(图4)。最大化合成NDVI年均值由0.029上升到0.054，增幅为0.025，平均年增速为0.208%(k=0.002 19,R2=0.806 73)； 最大EVI年均值从0.633上升到0.771，增幅为0.138，平均年增速为1.15%(k=0.007 73,R2=0.816 56)。结果显示在2000—2008年内2组数据均处于相对平稳状态，在此时间段内植被变化不大，在2008—2017年内2组数据显示的植被变化趋势处于大幅度的上升状态。

(a) 最大化合成NDVI年均值(b) 最大EVI年均值

利用时间分辨率为1 d的MODIS数据分析月际植被变化特征，最大EVI月均值的结果高值集中在每年5—10月(图5(a))，以此判断植被比较茂盛的时间段是在5—10月份，植被增长最快的时间约为每年6—7月。

(a) 最大EVI月均值(b) 月平均降水量和温度

植被变化时间特征与降水量、温度在每年的月际分布基本一致，平均降水量和温度在1月份出现最小值(0.3 mm，-8.4 ℃)，在6—7月份出现最大值(17.1 mm，19.0 ℃)。从5月份开始，温度逐渐上升，降水量也有所增加(图5(b))，且最大EVI月均值对降水、温度的变化亦有所响应，为了进一步反映它们之间的关联，则引入相关与偏相关分析方法。植被生长季开始前温度多在0 ℃以下，温度升高有利于植物萌发。在植物生长季期间，平均温度约为10 ℃，热量充足且降水量的增加有利于提高土壤的含水量，进而促进植被发育[20]。

3.2 植被变化影响因素分析

基于18 a间诺木洪洪积扇的MODIS卫星16 d合成的EVI产品MOD13Q1，1 a内会有比较多的合成时间点，对每2个合成时间点间的温度数据取平均，降水量取总和，初步形成对应时间点的变量集合(时间点也称为样本点，根据时间点筛选自然因素数据更有利于研究其与植被变化之间的关系)。为了使每个月的数据点具有代表性，对每个月的最大化EVI和温度数据取平均值，降水数据取累计便得到需要进行分析的最大化EVI、温度与降水变量集合(3个数据集合的变量个数均为411个)。

对最大化EVI值的平均数据集分别与平均温度数据集、降水数据集进行相关分析，结果显示EVI数据集与温度数据集之间的Pearson相关系数是0.839，显著性0.000<0.01,说明二者之间高度相关，即温度的变化在植被变化过程中起到很大的作用。EVI数据集与降水数据集之间的Pearson相关系数是0.457，显著性0.000<0.01,说明二者之间具有低度相关性，则该研究区域的植被对温度的敏感度要大于对降水的敏感度。

对3个变量集合进行一阶偏相关分析，EVI数据集与温度、降水数据集的一阶偏相关系数分别为0.816和0.327，显著性0.000<0.01。2个一阶偏相关系数均大于0，即EVI数据与温度、降水均呈现正相关的相关关系，且EVI数据和温度、降水之间的正相关具有显著性。

综合分析结果可知： 在2000—2017年间，EVI数据与温度之间具有显著的强正相关性，与降水量则是呈现出弱正相关关系，说明温度与降水是影响植被发育的重要因素，且影响该区域植被变化的主要因素为温度，而降水对植被变化的影响效果比较弱。原因是诺木洪洪积扇降水总体稀少，年降水量仅17.8～177.5 mm，植被多为农作物与人工防护林，生产用水多依靠地下水与远程高山融雪[21]。

3.3 枸杞种植区和盐碱化区植被变化与制约关系

3.3.1 重点区植被覆盖度变化趋势分析

诺木洪洪积扇植被在18 a间内整体变化并不明显，但枸杞园地和盐碱化区变化则是相对突出，故在研究区内圈定2个重点研究区——枸杞种植区和盐碱化区(图1斜线区域)，划分依据是以研究区2000年和2017年植被覆盖度等级图为基础，结合水系支流数量和方向，将诺木洪河左侧支流(哈西瓦河)、青藏公路和新扇扇缘部分划为良好绿洲——枸杞种植区。水是影响盐碱化区植被发育好坏的最重要的因素，一般河流尾部是植被变化相对剧烈的位置[22]，故将早期洪积扇的扇缘位置划为盐碱化区。

根据像元二分模型对枸杞种植区和盐碱化区近18 a内的植被覆盖度进行计算和拟合(图6)。

(a) 枸杞种植区(b) 盐碱化区

枸杞种植区的极低植被覆盖度等级在时间段内有很大程度降低(k=-14.41)，中等植被覆盖度有所降低，但变化程度不大(k=-2.02)； 相反高等植被覆盖度上升较多(k=10.90)，低等植被覆盖度也有所上升(k=5.58)。综合分析枸杞种植区极低覆盖度降低而高等植被覆盖度上升现象，即在2000—2017年内枸杞种植区植被发育较好。而盐碱化区的各个植被覆盖度等级变化均不大，高等植被覆盖度有降低的趋势(k=-0.45)，低等和中等植被覆盖度虽有略微上升趋势，结合覆盖度等级之间转变关系及该区域主要受高等植被覆盖度影响分析，盐碱化区域在近18 a内植被是有所退化的。

3.3.2 重点区植被覆盖度变化趋势验证

利用GEE平台在枸杞种植区和盐碱化区均匀生成验证点，计算近18 a的NDVI进行一元线性回归以获取NDVI值和时间回归方程的斜率值k，如图7所示。枸杞种植区的k值较大，即植被改善速率较快，而样本点10，15，23和55斜率值为负值，说明有些地块植被出现衰减，考虑其是不适合农作物生长的地块，通过人为因素的干预导致植被衰减。盐碱化区k值比较小，多数值处于[0,5]，说明盐碱化区域植被在此时间段内并没有很大的改善。相反，样本点1，8，16，17，37，45，51和54斜率值出现负值，说明盐碱化区域部分植被退化与消亡。

图7 2000—2017年间诺木洪洪积扇枸杞种植区和盐碱化区NDVI变化k值

综合分析18 a整个区域内植被逐渐发育，但枸杞种植区和盐碱化区植被改善速率差异较大，枸杞种植区植被变化速率k均值为13.23，下游植被变化速率均值为3.47。故考虑是由于枸杞种植区域植被受到人为作用的干扰[23]，导致用水量逐渐增大，导致盐碱化区域水资源减少，进而使得盐碱化区植被改善速率远低于枸杞种植区，反而有部分区域植被逐渐退化。

3.4 不同空间尺度的植被未来变化趋势分析

3.4.1 研究区植被未来变化趋势

利用Hurst指数对研究区未来变化趋势分析(图8)，研究区像元的H主要集中0.5～0.8区间之内，H≥0.5的像元占比为95.89%，H<0.5的像元占比为4.11%。故在未来一段时间内该区域植被变化与18 a来植被变化具有极大的一致性，洪积扇扇缘位置的持续性要低于扇中和扇根位置，原因是扇缘有发育较好的植被，植被发育受多方面因素影响，而扇中和扇根主要是戈壁和岩体，植被发育很低。

(a) Hurst指数分布(b) Hurst指数分布折线图

3.4.2 重点区植被未来变化趋势

枸杞种植区植被变化具有强持续性(图9(a))，H均值为0.759，其中：H≥0.5，即持续性比重为99.169%，H<0.5，即反持续性比重为0.764%，即表示未来一段时间内植被变化与18 a间植被变化趋势有很大程度的一致性。

(a) 枸杞种植区Hurst指数分布 (b) 盐碱化区Hurst指数分布

盐碱化区植被变化持续性相对于枸杞种植区要弱，但仍具有强持续性(图9(b))，H均值为0.624，H≥0.5，即持续性比重为95.788%，H<0.5，即反持续性比重为3.605%，即表示未来一段时间内盐碱化区植被变化与18 a间的植被变化趋势有较大程度一致性。

分析可知盐碱化区的持续性相对于枸杞种植区要低，由于盐碱化区除了受到温度和降水等气象因素的影响之外，更重要的还受上游枸杞地种植区植被变化影响，上游人为因素干扰使得上游植被进一步发育时，则会导致下游区域植被持续退化。

4 结论

通过对研究区与重点区进行综合分析其植被变化特征、植被变化成因以及未来趋势，可以得到以下结论：

1)2000—2017年间，诺木洪洪积扇植被向好的方向转变，但改善速率不大，主要原因是受到气象因子的影响，且温度的影响要远大于降水的影响。

2)枸杞种植区受到人为作用使得植被持续发育，枸杞种植区的植被增长对盐碱化区植被有着一定的制约效应，且在未来一段时间内会一直存在。

3)未来枸杞种植区与盐碱化区均具有明显的强持续性。

本文是依托于自然因素对植被变化的影响上来反映人类活动的作用，虽然采用18 a间的遥感和气候数据进行分析，时间上满足长时间序列分析的要求，但研究区内地表植被的变化十分复杂，物候信息在植被变化中也十分重要，若增加每个时间段内植被物候信息，分析结果精度将会进一步提高，也可以更好地说明人类活动的作用。