基于长时间序列NDVI的黄土高原延河流域及其沟壑区植被覆盖变化分析

贾云飞, 李云飞, 范天程, 曾 竞, 赵建林

(长安大学 地质工程与测绘学院, 西安 710054)

黄土高原位于我国中部偏北部，属半干旱大陆性季风气候区，气候变化敏感，生态环境脆弱[1]，是中国水土流失的严重区域之一。严重的水土流失会影响社会经济的发展，高泥沙输移会威胁黄河下游地区[2]。相关研究表明黄河近90%的输沙量来源于黄土高原[3]。

为了控制水土流失和高输沙量，过去50 a中国在黄土高原实施了一系列水土保持措施：其中包括坡面区域的大量梯田建设和生态修复项目，以及沟壑区域大量淤地坝的建设[4]。其中最为显著的生态修复项目是自1999年开展的“退耕还林(草)”。相关研究表明自“退耕还林”政策实施以来，黄土高原植被覆盖度显著提高，植被状况明显好转；郭敏杰等[5]研究表明1982—2006年黄土高原植被覆盖度呈现稳定的增加趋势，增速为0.075%/a；黑哲[6]研究表明2000—2019年超过90%的区域植被NDVI呈增加趋势，其中71.56%的地区呈显著增加；郭永强等[7]研究表明实施“退耕还林”(2000—2015)后，植被覆盖率年均增幅为0.59%，年植被覆盖率从1987年的41.78%提高到2015年的53.23%。胡春宏等[8]研究表明新中国成立后，经过治理，入黄沙量由1919—1959年16亿t/a锐减至2000—2018年约2.5亿t/a。黄土高原具有典型的丘陵沟壑地貌，而侵蚀过程和泥沙主要来源于沟壑区与坡面。其中坡面以面蚀、细沟侵蚀以及耕作侵蚀为主，而沟壑区主要以沟蚀和重力侵蚀为主，研究表明在黄土高原流域中，沟壑间侵蚀产沙大于坡面产沙。Zhao等[9]研究表明当沟壑区的密度大于30%时，沟壑区对流域泥沙的贡献大于75%。以往的研究表明坡面的植被覆盖情况有了明显的改善，但坡面的植被覆盖情况改善是否会诱发坡面间的沟壑区植被覆盖增加，目前还尚未有研究表明。

本文以黄土高原退耕还林重点区域延河流域为研究对象，基于Google Earth(GE)平台采用系统样本法和人工交汇法提取该流域的沟壑区与坡面地貌；基于长时间序列NDVI数据，采用趋势分析法、Pettitt突变点分析、残差分析法以及地学统计方法，在分析延河流域过去20 a(2000—2019年)整体的植被覆盖变化情况以及降雨量和人类活动对NDVI变化的影响的基础上，分析黄土高原坡面和沟壑区的植被覆盖变化情况，探索坡面与沟壑区植被变化的差异性及相关性。

1 研究区域概况

本文的研究区域是延河流域，延河是黄河的一级支流，全长286.9 km。延河流域(108°45′—110°28′E，36°23′—37°17′N)位于黄土高原中部，流域面积7 725 km2，地势西北高、东南低。延河流域属于黄土丘陵沟壑区，地表破碎，沟壑密度在2.1～4.6 km/km2，水土流失严重。流域的植被区划属于森林草原地带，主要是人工种植的次生植被和干旱草本植物。该流域属暖温带半湿润地区，年平均降水量约为520 mm，降水主要集中在6—9月份，约占全年降雨量的75%。

2 数据来源与研究方法

2.1 数据来源

2.1.1 NDVI数据 本文使用归一化植被指数(Normalized Difference Vegetation Index，NDVI)来分析植被覆盖变化。NDVI与植被的初级生产力、叶覆盖、生物量具有很好的相关性，已被广泛用于量化植被的生长趋势和过程[10]。本文使用的NDVI数据为MODIS 16 d合成的NDVI产品MOD13Q1，空间分辨率为250 m，时间跨度为2000—2019年，由Google Earth Engine平台提供，使用最大值合成法合成影像，得到延河流域各年的NDVI最大值栅格数据集，并将NDVI像元初始值转化到-1～1[11]。

2.1.2 降雨数据 为分析降雨对延河流域NDVI植被覆盖变化的贡献，本文采用1 km分辨率的降雨数据集，该数据来源于中国科学院资源环境科学与数据中心的中国1980年以来逐年年降水量空间插值数据集，该数据集是基于全国2 400多个气象站点日观测数据，通过整理、计算和空间插值处理生成。本文使用了该数据集2000—2019年的降雨数据。

2.1.3 沟壑区样本提取与栅格化 为了区分沟壑区和坡面植被覆盖变化情况，本文首先在延河流域提取沟壑区样本，提取过程如下：首先，在延河流域系统均匀地分布272个3 km×3 km正方形样本区域，该样本覆盖整个延河流域。然后基于Google Earth平台，人工勾绘出每个样本区内的沟壑区域，并假设样本区内除沟壑区域外就是非沟壑区域即坡面区域。在此基础上，对提取的矢量沟壑区和坡面进行栅格化。栅格化过程如下：首先创建与NDVI像元一致的缓冲区，基于缓冲区对沟壑区矢量与坡面区矢量进行分割，使得分割后的每一个区域能够正好覆盖所对应的像元，然后计算得到分割后的沟壑区和坡面区内像元对应的矢量面积，再结合像元面积得到每个像元的沟壑区面积占比与坡面面积占比。具体流程见图1。

注：A为272个样本区中的某一样本区；B为基于Google Earth人工勾绘的样本区内沟壑区域；C，D为沟壑区矢量和坡面区矢量分割结果；E，F为沟壑区和坡面区像元面积占比。

由于本文采用的NDVI数据的空间分辨率为250 m，故栅格化后的沟壑区和坡面的空间分辨率与NDVI像元一致，栅格化后的像元值为对应栅格的沟壑区面积占比(GD)，计算公式如下：

(1)

式中：GD为每个像元值，即沟壑区面积占比；Sgully为每个像元中人工勾绘沟壑的面积(m2)。

坡面面积占比(SD)计算公式如下：

(2)

式中：SD为每个像元值，即坡面面积占比；Sslope为每个像元中非沟壑区即坡面的面积(m2)。

2.2 研究方法

2.2.1 延河流域沟壑区和坡面NDVI 为分析延河流域沟壑区和坡面NDVI变化趋势，本文分别计算了272个样本区2000—2019年的样本区(NDVImean)，坡面(NDVIslope)和沟壑区(NDVIgully)的平均NDVI年最大值。其中，样本区平均NDVI等于样本区内各像元的平均值。沟壑区NDVI采用面积加权法计算，计算公式如下：

(3)

式中：i=1,2,…,272；NDVIgully(i)为第i个样本区的沟壑区平均NDVI最大值；NDVIij为第i个样本区第j个像元的年NDVI最大值；GDij为第i个样本区第j个像元对应的沟壑区面积占比。

坡面NDVI同样采用面积加权法计算，计算公式如下：

(4)

式中：i=1,2,…,272；NDVIslope(i)为第i个样本区的坡面平均NDVI最大值；NDVIij为第i个样本区第j个像元的年NDVI最大值；SDij为第i个样本区第j个像元对应的坡面面积占比。

在此基础上，本文采用ANOVA分析方法，分析延河流域沟壑区NDVI(NDVIgully)和坡面区NDVI(NDVIslope)是否具有显著差异(以p<0.05检验其显著性)。以及采用线性回归方程和R2分析样本区中两者相关性即：

NDVIgully=a+b×NDVIslope

(5)

式中：a和b为线性回归系数，如果b>1则表明样本区中沟壑区NDVI大于坡面NDVI，反之沟壑区NDVI小于坡面NDVI。

2.2.2 趋势分析法 本文采用一元线性趋势分析法分析整个延河流域20 a的NDVI变化趋势。基于像元尺度拟合NDVI变化，并整合每个像元的时变特征和区域空间变化，以反映研究区域NDVI的20 a时空格局演变[12]。一元线性趋势线斜率的计算公式为：

(6)

式中：n为研究数据的时间长度，本文为20 a；NDVIj为整个延河流域第j年的NDVI最大值；k为像元回归方程趋势线斜率，若k>0，表示研究时段内NDVI变化呈增加趋势，反之，则呈减少趋势[13]。为了更好地判断整个研究区域植被覆盖的动态变化趋势，根据k值的大小，可定义7个NDVI趋势变化水平[14]：重度退化、中度退化、轻度退化、基本不变、轻度改善、中度改善以及明显改善(表1)。

表1 NDVI趋势变化水平

2.2.3 突变点分析 为了检验延河流域过去20 a植被覆盖变化是否连续，本文采用Pettitt[15]突变检验方法检验延河植被覆盖变化是否存在显著突变点(p<0.05)。Pettitt突变检验是一种非参数检验法，其假设样本容量中存在突变点t，采用Mann-Whitney的统计量Ut,n检验样本序列突变点t前后两个子样本，x1,…,xt和xt+1,…,xn二者累积分布是否存在显著差异。

(7)

式中：xt，xi分别为第t年和第i年的NDVI相邻时间段斜率样本；n为数据系列长度；Ut,n为将根据第一个样本序列超过第二个样本序列次数的统计组成新的序列。

Pettitt法原假设H0为序列不存在突变点。若t时刻满足：

(8)

则t点处为突变点。

同时计算统计量：

(9)

若概率p≤0.05，则认为检测出的突变点在统计意义上是显著的。

2.2.4 NDVI变化驱动力分析 由于延河流域过去20 a人类活动频繁，特别是退耕还林项目的实施，同时降雨也会对区域植被覆盖产生影响。因此，本文采用由Evans等[16]提出的残差分析方法研究降雨量和人类活动对延河流域NDVI变化的影响及相对贡献。该方法的主要步骤为：首先基于2000—2019年NDVI和降雨量时间序列数据，以NDVI为因变量，以降雨量为自变量，建立一元线性回归模型，计算模型中的指标参数(斜率和截距)；基于降雨量与回归模型的参数，计算得到NDVI的预测值(NDVIPRE)，表示降雨量对NDVI变化的影响；最后计算NDVI观测值与NDVIPRE的差值，即残差(NDVIHA)，表示人类活动对NDVI变化的影响。具体计算公式如下：

NDVIPRE=c+d×PRE

(10)

NDVIHA=NDVIOBS-NDVIPRE

(11)

式中：NDVIPRE为回归模型的NDVI预测值，即降雨量对NDVI变化的贡献；PRE为年降雨量(mm/a)；NDVIHA为残差，即人类活动对NDVI变化的贡献；NDVIOBS为基于NDVI栅格数据集的观测值；c，d为模型参数。

3 结果与分析

3.1 延河流域NDVI时空变化特征

(1) 2000—2019年延河流域年平均NDVI整体呈现增加趋势，在研究时段内，延河流域年平均NDVI由2000年的0.415增加到2018年的0.750，其增长率达到44.60%，且年平均趋势率为1.30%(p<0.001)，相关系数R2=0.544，表明2000—2019年延河流域的植被覆盖状况改善较为明显。基于Pettitt′s突变点分析表明延河流域NDVI突变年份在2008年(p<0.001)，因此将延河流域NDVI变化整体趋势分为两个阶段：2000—2008年和2009—2019年。2000—2008年延河流域年平均NDVI增加显著，平均趋势率为2.00%(p<0.001)，相关系数R2=0.343；2009—2019年延河流域年平均NDVI增加趋势减缓，平均趋势率为0.70%(p<0.001)，相关系数R2=0.105(图2)。

(2) 2000—2019年延河流域NDVI变化趋势具有一定的空间异质性，延河流域大面积植被覆盖恢复良好，年平均NDVI呈增加和减小的区域分别占98.90%，0.80%。其中呈现明显改善(k>0.009)的区域面积占比约为81.27%，约6 314.20 km2；其次是中度改善的区域面积占比约为14.44%，约1 121.84 km2，主要分布在西北部、西南部及东南部；重度退化(k<-0.009)的面积占比约为0.31%，约24.42 km2(表2)，主要分布在中部地区，该地区为延安市区。

图2 2000-2019年延河流域NDVI年际变化及分段变化

通过对比可以发现2000—2008年延河流域植被覆盖状况增加趋势更加明显，整个延河流域植被覆盖情况都有所改善，退化的区域面积可以忽略不计；而2009—2019年植被覆盖状况增加趋势明显减缓，且具有很大的空间异质性；除西北部、中部以及东南部地区均出现不同水平的退化趋势外，其他地区的增加趋势也有所减缓(图3)。2000—2008年延河流域年平均NDVI呈增加和减小趋势的区域面积分别占总面积的98.57%，0.90%。其中呈明显改善(k>0.009)的区域面积占比约为88.85%，约6 903.16 km2；其次是中度改善的区域面积占比约为7.41%，约575.86 km2；退化的面积占比都在0.50%以下。2009—2019年延河流域年平均NDVI变化趋势具有很大的空间异质性；年平均NDVI呈增加和减小趋势的区域面积分别占总面积的87.13%，7.33%。其中明显改善的区域面积占比约为33.70%，约2 617.78 km2；其次是中度改善的区域面积占比约为34.01%，约2 642.71 km2；再次是轻度改善的区域面积占比约为19.42%，约1 508.60 km2；基本不变的区域面积约占5.54%，约430.51 km2；轻度退化的区域面积约占4.65%，约361.20 km2(表2)。

图3 不同年份延河流域NDVI变化趋势空间分布对比

表2 2000-2019及分段时期延河流域趋势变化分布

3.2 NDVI变化驱动力分析

延河流域NDVI变化驱动力分析表明2000—2019年延河流域降雨量对NDVI变化的影响一直保持稳定，平均趋势率为0.037%/a(p<0.05)，且相关性较弱(R2=0.001)；而2000—2019年人类活动对NDVI变化的影响呈现明显增加的趋势，平均趋势率为1.29%/a(p<0.001)，相关性较强(R2=0.755)。这在一定程度上说明人类活动对延河流域植被覆盖的改善产生了积极影响。人类活动对延河流域NDVI变化的影响也呈现分段特征，2000—2008年人类活动对延河流域NDVI变化的贡献明显大于2009—2019年(图4)。

3.3 延河流域沟壑区与坡面植被覆盖变化特征

基于272个样本区内平均沟壑区NDVI和坡面NDVI分析表明，随着坡面NDVI的增加，沟壑区NDVI也呈现明显增加的趋势，沟壑区与坡面的线性趋势率为1.02(p<0.001)，即沟壑区植被覆盖要略大于坡面植被覆盖，两者相关系数R2=0.97(图5)。

然而，两者相关性存在年际差异，且与延河流域整体NDVI变化趋势存在关联。如图6所示，基于每年样本区中沟壑区NDVI和坡面NDVI线性相关分析表明，2000—2008年，各样本区中坡面NDVI与沟壑区NDVI的线性趋势率均≥1，且相关系数R2均>0.90，即当整体植被覆盖呈现显著增加时，沟壑区植被覆盖要好于坡面区域植被覆盖，且两者相关性较强；而2009—2019年，各样本区中坡面NDVI与沟壑区NDVI的线性趋势率均<1，即沟壑区植被覆盖小于坡面区域植被覆盖，且相关系数R2除了2010年、2011年、2015年外，均小于0.90。

注：NDVIP，NDVIH分别为降水、人类活动对NDVI变化的贡献。

图5 坡面NDVI与沟壑区NDVI线性相关

非独立样本t检验的结果表明沟壑区NDVI与坡面NDVI存在明显差异(p<0.05)。2000—2019年沟壑区NDVI与坡面NDVI差异的均值为0.012(p<0.001)；2000—2008年沟壑区NDVI与坡面NDVI差异的均值为0.009(p<0.001)；2009—2019年沟壑区NDVI与坡面NDVI差异的均值为0.014(p<0.001)。

4 讨论与结论

4.1 讨 论

本文研究表明，2000—2019年延河流域整体植被覆盖改善较快，但具有一定的空间异质性，中部地区城市建设导致植被覆盖减少，其他地区都植被覆盖都呈现增加趋势。延河流域植被覆盖变化具有分段特征，2000—2008年植被覆盖恢复速度较快，2009—2019年植被覆盖恢复速度有所减缓，恢复趋势趋于平稳。退耕还林还草分为1999年起实施的前一轮退耕还林还草和2014年起实施的新一轮退耕还林还草。1999—2001年是退耕还林还草的试点示范阶段，2002年全国全面启动退耕还林还草工程，2008年起开始巩固各工程省区的退耕还林成果，为扩大退耕还林、退牧还草规模，2014年开始了新一轮退耕还林还草[17-18]。与此同时降雨和人类活动对NDVI变化归因分析表明，降雨量对延河流域NDVI变化的贡献保持稳定，而人类活动对延河流域NDVI的贡献有着明显增加的趋势。1999年来开展的如“退耕还林”等一系列的生态工程是2000—2019延河流域NDVI明显增加的主要原因。

自退耕还林项目实施以来，延河流域丘陵沟壑区植被恢复态势较为显著，坡面植被覆盖状况也得到了明显改善[19]。基于人工提取的大量沟壑区和坡面样本，本文初步分析了退耕还林项目对延河流域沟壑区和坡面区植被覆盖影响。初步结果表明样本区中沟壑区与坡面区的植被覆盖具有同步性，两者相关性较强，这表明当坡面植被覆盖增加时，也间接增加了沟壑区植被覆盖。这与前人研究一致，相关研究表明沟壑区由于弃耕时间长且人为扰动较小，形成了比较稳定的植物群落，因此在坡面植被覆盖改善时，沟壑区的植被覆盖也随之改善，且沟壑植被覆盖优于坡面植被覆盖。然而，坡面上的植被恢复与坡度、水土保持政策密切相关，近年来也得到了极大的改善[20]。本文研究表明，延河流域坡面NDVI与沟壑区NDVI具有较强的相关性，且同样具有分段特征，坡面NDVI增加的趋势越快，沟壑区NDVI增加的趋势也越快，且相关性较强。

图6 2000-2019各年坡面NDVI与沟壑区NDVI的线性关系

4.2 结 论

(1) 2000—2019年延河流域植被NDVI增加显著，平均趋势率为1.30%/a(p<0.001)，其中，中部地区是城市用地NDVI呈减小趋势；其他地区NDVI都呈现增加趋势，北部相比南部NDVI增加较明显。2000—2008年延河流域植被NDVI增加显著，平均趋势率为2.00%/a(p<0.001)，且整个延河流域植被覆盖都呈现改善状况；2009—2019年植被NDVI增加趋势减缓，平均趋势率为0.70%/a(p<0.001)，且西北部、中部和东南部植被覆盖呈现轻度退化趋势。

(2) 降雨和人类活动对延河流域NDVI变化的贡献具有明显的差异性，2000—2019年降雨对NDVI变化的贡献保持稳定，且相关系数R2=0.001，相关性较弱；而人类活动对NDVI变化的贡献较大，人类活动对延河流域NDVI变化的贡献的平均趋势率为1.30%/a(p<0.001)，相关系数R2=0.755，存在较强的相关性。

(3) 2000—2019年延河流域沟壑区域NDVI随着坡面NDVI的增加而增加，年平均趋势率为1.02(p<0.001)，且相关系数R2=0.967，相关性较强。坡面NDVI与沟壑区NDVI的线性关系也呈现分段趋势，2000—2008各年坡面NDVI与坡面NDVI的平均趋势率均>1，且相关性较强，2009—2019各年的平均趋势率均<1，且相关性较2008年以前有所减弱。