基于MODIS NDVI数据的陕西省植被绿度时空变化及人类活动影响

殷崎栋,柳彩霞,田 野

1 陕西生态环境规划设计院有限公司, 西安 710075 2 中国科学院空天信息创新研究院遥感科学国家重点实验室, 北京 100101 3 中国科科学院城市环境研究所城市环境与健康重点实验室, 厦门 361021

陆地植被是大气、岩石圈、土壤圈、水圈和生物圈相互作用的产物。它是连接气候变化、人类活动和生态系统的一个重要环节。气候和人类活动是控制和影响植被空间分布及其变化的基本驱动力[1-2],辨别气候和人类活动影响区域对于理解和管理地表植被非常重要,特别是在干旱和半干旱地区,因为在这些地区,年降雨量相对较低,年降雨量波动性较大。在辨识人类活动造成的影响程度时,气候的变异性很难将植被中的自然变化与直接人类活动引起的自然变化分开[3- 6]。鉴于土地退化或者绿化工程可能在几十年内逐步发生,对植被的观测时间序列必须一致,因此遥感技术成为监测并分析植被变化的一项可靠技术。遥感观测的植被指数,由于具备统一的观测基准和长时间序列特征,可以定量描述植被覆盖和植被生产力的变化。利用归一化差异植被指数(NDVI)和气候因子的关系,可以研究人为和自然因素引起的植被变化[6- 8]。RESTREND算法是一种从NDVI趋势中消除气候影响的方法[9],基于气候数据和植被指数进行土地退化的检测[3, 9]。早期的RESTREND算法通过计算年最大NDVI、生态系统生产力和降水之间的线性回归来控制影响植被的气候要素的变化[4, 9],然后从线性回归关系中预测的NDVI与观测NDVI之间的差异,得到NDVI残差,之后对残差进行线性趋势分析,该分析结果被认为是人类活动引起的地表植被变化。然而,由于降水量和植被指数之间的线性关系存在不稳定性。当在时间序列的中间发生退化时,这种线性关系可能会发生改变,导致不可靠的结果[5]。因此,检测植被变化是否引起植被与降水之间的线性关系(也称为断点)至关重要。遥感产品在时间序列上的断点检测已应用于土地覆盖制图[10- 12]和森林管理[13-14]。已有研究表明,当在植被与降水之间的关系中断时检测断点时,可以改善植被变化检测精度[3,5,15]。时间序列分割和残差趋势分析(Time Series Segmentation and Residual Trend, TSS-RESTREND)是由Burrell, Evans[6]提出,加入了断点检测的RESTREND算法,成功用于检测澳大利亚全国的植被变化情况[6,8],能够改进对变化区域和变化方向的检测。

陕西省地处中国大陆腹地,是水土流失、沙化等环境问题最严峻的地区之一。自1999年起,陕西省相继启动并实施了生态环境建设综合治理工程、天然林资源保护工程、退耕还林工程、重点防护林工程、水土保持工程和天然草场恢复与建设工程等一批重点生态建设工程,使陕西省特别是陕北地区植被状况有了明显好转,抑制了陕北毛乌素沙漠南扩的势头[16]。最新研究表明,近10多年来,陕西省部分区域的植被生长发生了重大变化。如栾金凯等[17]利用MODIS/Terra NDVI时间序列遥感影像(2000—2016年)和温度、降水气候数据,分析了17年来陕西省榆林市的植被指数空间分布变化规律,并且利用复直线回归法揭示了气象因素对植被生长和演化的影响,榆林市一半以上的区域人类活动对植被生长起到了促进作用。邓晨晖等[18]基于MODIS/Terra NDVI时间序列遥感影像(2000—2015年)、DEM和气象数据,采取趋势分析、多元回归残差法、偏最小二乘回归法,反映了秦岭地区植被覆盖度变化特征,探究其对气候变化和人类活动的双重响应机制。在区域尺度上,更早的研究集中在关中[19]、陕北[20- 24]等地。在全省尺度上,韩红珠等[25]利用MODIS数据提取了2001—2016年间的植被物候时空特征,结果表明退耕还林工程效益显著[26-27]。目前,还没有相关研究报道利用公里空间分辨率数据从宏观尺度分析陕西全省的植被绿度时空变化特征及其人类活动影响区域分析,人类活动对陕西全省植被的影响程度和范围还缺乏全局了解。本文从全省尺度,细分为陕北、关中平原、陕南地区,利用可靠的TSS-RESTREND算法解释人类活动对陕西省NDVI(2001—2018年)的贡献影响程度,算法剥离了气候变化本身对NDVI的贡献影响,对学者、管理者摸清陕西全省植被绿度的长时间序列变化具有重要意义。同时,2019年是陕西省全面实施退耕还林生态建设工程20周年,利用产品一致性较好的遥感时间序列产品,全面分析近20年来全省植被变化情况和客观评价生态建设工程成效,可为后续生态环境管理提供可靠的本底数据。

1 研究区概况

图1 研究区位置意图Fig.1 Location map of the study area

陕西省地理范围为105°29′—11°15′E, 31°42′—39°35′ N之间(图1),地处东部湿润地区和西部干旱区的交界地带,省内气候类型多样,各地的年均气温在7—6℃,属于典型的大陆性季风气候区。省内地形地貌类型丰富,地势特点是南北高,中间低,由西向东倾斜,年降水量受山地地形影响显著。自北向南可粗略划分为三大自然区域,分别为陕北黄土高原区、关中平原区和陕南秦巴山地区。陕北黄土高原区包括榆林地区和延安地区,海拔900—1600 m,地势较高,地处干旱、半干旱地区,年均气温7—11℃,年降水量约400—600 mm,地区分配不均,季节变化大。关中平原区由渭河干流以及两岸支流日久冲积、沉积形成,主要包含西安市、铜川市、宝鸡市、咸阳市和渭南市,地势相对较低,海拔在300—800 m之间,中部较为平坦,年均气温11—13℃,年降水量约500—700 mm,为温带半湿润气候。陕南秦巴山地处于秦岭和大巴山系,海拔200—1200 m,主要包括汉中市、安康市、商洛市,水热条件相对较好,年均气温14—15℃,年降水量约700—900 mm。同时,陕西省地处黄河中游、长江上中游地区,属于国家生态环境建设重点区域。

2 数据和方法

2.1 数据来源与预处理

本文使用的NDVI数据是2001年1月到2018年12月的MODIS NDVI (MOD13A2)数据,出自美国国家航空航天局NASA发布的EOS/MODIS数据产品,空间分辨率为1 km×1 km,时间分辨率为16 d。MODIS/NDVI产品经过了水、云、重气溶胶等预处理,保证了数据质量,被广泛应用于区域植被覆盖变化研究。本文数据处理在谷歌地球数据云平台 (Google Earth Engine) 中处理,NDVI数据获取ID为MODIS/006/MOD13A2,利用数据集的质量控制字段(SummaryQA),选取所有最优观测值(SummaryQA等于0)的最大值作为月NDVI。

降雨数据为1999—2018年陕西省及其周边省份33个气象台站的月数据,来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn/)。利用薄板样条空间插值算法[28-29],得到1999年1月至2018年12月的1 km空间分辨率的逐月降水栅格数据。

2.2 生长季NDVI和降水观测时间序列统计

计算观测异常值。以生长季(6、7、8月)降水异常值计算为例,利用生长季每个月的降雨插值图,计算如下变量：

(1)计算每一年生长季降水均值,即3个月(6、7、8月)的平均值(mm)；

(2)计算2001—2018年生长季降水均值(mm)；

(3)用(1)除以(2)得到异常百分比(%)

同样地,计算生长季NDVI异常百分比。为了规避观测数据缺失带来的误差,我们没有计算6、7、8月的平均值,而只用8月的NDVI来代表生长季的NDVI。

2.3 TSS-RESTREND方法

标准RESTREND是由Evans and Geerken[9]开发并经Wessels, Prince[4]等人进一步修改。算法基于像素的分析,用于揭示气候要素和人类活动导致的植被生态系统的退化过程及其规律[3, 30]。利用最小二乘回归(Ordinary Least Squares Regression, OLSR),计算植被净初级生产力(Net Primary Productivity, NPP)与降水量的关系,即植被—降水关系(Vegetation Precipitation Relationship, VPR),其中NPP用最佳生长季节NDVI(NDVImax)表示[5-6],降水量用最优累积降水量表示。根据OLSR模型得到的预测NDVImax与每次观测NDVImax的差值称为VPR—残差[9]。

yi=β0+β1xi

(1)

式中,yi是VPR—残差,xi是年份,β0是截距,β1是斜率。

一般情况下,认为VPR—残差的趋势与降水无关,是地表土地变化的表征[9]。利用标准RESTREND算法,必须满足以下3个条件[6]：

1)VPR关系必须是正向(斜率>0)且显著的。对显著的推荐值为：在P<0.05水平上,R2>0.3。

2)VPR—残差在时间序列上是单调函数[31]。

3)VPR与时间高度相关,即VPR在时间序列上具有可比性,例如在研究区间内生态系统没有发生重大结构性变化[5]。

在我们的研究中,最佳生长季NDVI发生在夏季,最优累积则是通过逐像素计算累积周期(1—12个月)和偏移期(1—3个月)的组合与NDVImax的相关系数,相关系数最高的降水累积量为最优累积降水。标准RESTREND算法有一定的局限性,当变化的速度和方向在时间序列中发生变化时,无法识别变化趋势[32]。

TSS-RESTREND (Time series segmentation and residual trend) 能克服上述问题,最先被Burrell, Evans[6]提出并用于澳大利亚全境的植被变化分析,该方法的主要目的是首先利用BFAST(Breaks for Additive Season and Trend)算法[33-34]检测出NDVI时间序列的断点。不像之前的研究应用[13, 35],在TSS-RESTREND中的BFAST算法不直接用于NDVI时间序列,而是用于NDVI—降水关系的残差时间序列上,在此应用中关闭BFAST的季节成分(Seasonal Component).

TSS-RESTREND首先在研究时间段内利用BFAST算法找出NDVI时间序列的不一致年份。在每一个子时间段内,利用Chow检验准则判定不一致年份对植被NDVI的影响程度,当F检验在α=0.05水平时,BFAST检测的不一致年份并没有改变植被的结构,BFAST检测结果被驳回,对检测出来的时间段进行前后合并。TSS-RESREND有四种情形：

(1)当VPR (α=0.05) 和 VPR—残差 (α=0.05)显著相关, 符合标准RESTREND算法。

(2)当在VPR—残差中检测到不一致年份,Chow检验同样要再用于VPR。当VPR—残差(α=0.05)显著而VPR (α=0.05)不显著时,采用分段RESTREND算法(segmented RESTREND):

yi=β0+β1xi+β2zi+β3xizi

(2)

式中,xi为年份,zi为哑变量(0 或者1).β0为截距,β1为斜率,β2为断点处的抵消项,β3为断点处的斜率。

(3)当一个像元在VPR中有显著断点时,可能预示着生态系统的结构性改变。在这种情况下,假定降雨对断点两侧的影响是一致的假设就不合理,因此,时间序列NDVI需要分段分析。为了使两侧的数据有可比性,可以对降雨进行归一化得到一个标准项:

(3)

式中,zi为标准项,xi为观测值,μ为时间序列均值,σ为标准差。利用这个标准项,可以建立如下多变量回归关系：

yi=β0+β1xi+β2zi+β3xizii=2001,…,2018

(4)

式中,xi为式(3)计算的降水标准项,zi为哑变量(0 或者1)。β0为截距,β1为斜率,β2为断点处的抵消项,β3为断点处的斜率。对一个像素而言,像元绿度的变化需要包含残差的变化。

表1 NDVI变化分类

根据线性拟合得到的F检验的显著性水平和斜率，可以将NDVI变化趋势分为增加、下降和基本不变(No Significant Change, NSC)；其中增加细分为剧烈增加(Dramatic Increase,I1)(α<0．01)， 明显增加(Obvious Increase, I2)(0．01≤α<0．025)， 中等增加(Moderate Increase, I3) (0．025≤α<0．05)和轻微增加(Slight Increase, INC) (0．05≤α<0．1)。下降细分为剧烈下降(Dramatic Decrease, D1)(α<0．01)， 明显下降(Obvious Decrease, D2)(0．01≤α<0．025)， 中等下降(Moderate Decrease, D3)(0．025≤α<0．05) 和轻微下降(Slight Decrease, DNC)(0．05≤α<0．1)

(4)不符合上述条件的像素被认为不适宜用TSS-RESTREND算法。

根据上述回归关系中的β1和显著性水平α[7, 36],根据NDVI残差与时间序列的线性拟合,我们利用F检验的α值将NDVI的变化趋势分为9个等级(表1)。NDVI下降的趋势等级分为4个等级,分别为D1 (α<0.01), D2 (0.01≤α<0.025), D3 (0.025≤α<0.05) 和 DNC (0.05≤α<0.1),NDVI增加的趋势等级分为4个等级,分别为I1 (α<0.01), I2 (0.01≤α<0.025), I3 (0.025≤α<0.05)和 INC (0.05≤α<0.1)。NSC为显著性水平不通过(α>0.1),表示NDVI无变化。

3 结果与讨论

3.1 2001—2018年生长季NDVI时空分布与统计分析

图2展示了陕西省生长季最大NDVI在2001—2018年间的空间分布。从图中可以看出,NDVI大小随纬度变化形成南北差异,随地形起伏变化形成区域差异,总体NDVI的空间分布特征与陕西省不同的地形地貌空间(图1)具有较好的一致性,NDVI时间序列南高北低,关中平原包含少量低值区。

时间维度上,全省NDVI平均值的年增长率为0.006/年,如图3所示。由于2011年研究区域生长季受云影响较大,缺失值较多,因此统计时去除2011年的NDVI值。该增长率超过了三北防护林(0.0007/a)[37]、黄河上游(0.0023/a)[38]和三江源地区(0.0001/a)[39],这说明退耕还林工程对陕西省的影响非常显著。截至2018年,相比18年来的平均值(图4),77.29%的区域大于均值。其中,陕北的榆林、延安大于均值的区域较大,分别为97.52%和89.03%(表2),秦巴山区次之,为73.91%。特别是2012年之后, NDVI高值向北逐年推进趋势明显(图2)。

图2 2001—2018年陕西省年最大归一化植被指数Fig.2 Maximum NDVI(Normalized Difference Vegetation Index) in Shanxi in August from 2001 to 2018

表2 NDVI与18年均值(2001—2018年)比较

图3 2001—2018年陕西省生长季年NDVI曲线Fig.3 Variations in annual maximum NDVI and abnormal (%) in Shaanxi in August from 2001 to 2018

图4 2001—2018年陕西省生长季NDVI异常百分比Fig.4 Maximum NDVI abnormal (%) in Shaanxi in August from 2001 to 2018

3.2 人类活动造成的NDVI变化趋势与分级

利用TSS-RESTREND算法计算得到NDVI变化量如图5所示。2001—2018年间,陕西省植被绿度受剧烈的人类活动影响,NDVI呈增加态势,NDVI增加的区域达71.77%,而陕北地区的增加量明显大于关中平原区和陕南秦巴山地,其中陕北的榆林NDVI增加区域为72.11%,延安为86.44%,均超过了全省平均水平,这是由于陕西省北部NDVI增长与人工造林面积密切相关,有关研究报道达96%以上,说明陕西省北部实施的退耕还林生态工程效果显著,据统计,2000—2014年间,陕北地区山地和沙地造林面积达238.6×102km2[40]。

NDVI减少区域主要分布在城市周边,如西安市、咸阳市、汉中市、安康市。如图所示。关中平原NDVI减少量区域面积最大,占该区域的6.04%,秦巴山地次之(1.20%)。

图6 2001—2018年陕西省绿度分级Fig.6 Greenness level determined by TSS-RESTREND in Shaanxi from 2001 to 2018 根据线性拟合得到的F 检验的显著性水平和斜率，可以将NDVI 变化趋势分为增加、下降和基本不变(No Significant Change, NSC)。其中增加细分为剧烈增加(Dramatic Increase,I1)(α<0．01)， 明显增加(Obvious Increase, I2)(0．01≤α<0．025)， 中等增加(Moderate Increase, I3) (0．025≤α<0．05)和轻微增加(Slight Increase, INC) (0．05≤α<0．1)。下降细分为剧烈下降(Dramatic Decrease, D1)(α<0．01)， 明显下降(Obvious Decrease, D2)(0．01≤α<0．025)， 中等下降(Moderate Decrease, D3)(0．025≤α<0．05) 和轻微下降(Slight Decrease, DNC)(0．05≤α<0．1)

(3)按照表1所示的绿度分级规则,可以将全省人类活动引起的NDVI变化量定性分为增长(INC,I1,I2,I3)、减少(D1,D2,D3)和无明显变化(NSC)8类,总体上陕西全省呈变绿趋势(图6)。榆林市和延安市的变绿区域明显多于关中平原和秦巴山地(图6,图7),在空间上呈连续分布,而关中平原和秦巴山地的变绿区域大部分呈点状分布。表3所示的TSS-RESTREND检测出的剧烈增长区域,延安为55.46%(11935.9km2),榆林为34.34% (13825.7 km2),而陕南为41.03% (14746.8 km2),这说明处于湿润气候区的陕南地区也有显著变绿趋势,与已有研究结果一致[18, 41]。在这里需要说明是,当NDVI和降雨量关系不符合TSS-RESTREND算法适用条件时,该像素不做任何处理。因此,表3的面积统计未包含算法不适用像元。虽然陕西省位于东部湿润和西部干旱过渡的半干旱地区,降雨要素是影响植被生长变化的首要气候要素,但是局部区域,如陕南秦巴山地,气温和水热因素也是抑制植被生长的重要要素。虽然TSS-RESTREND算法可以在算法层面剔除算法不适用的像元,但是人类活动对陕西南部植被的影响结果还须进一步论证。TSS-RESTREND算法内置的BFAST时间序列分析方法可以检测出VPR关系变化年份,从而可以扩大RESTREND算法的应用区域,并且可以反映NDVI变化带来的生态累积效应。图7展示了自2001—2018年,TSS-RESTREND检测的VPR断点年份。整体上,断点年份发生在2008年之后,这说明自2008年后,原有的植被NPP与降水的关系被打破,常规的RESTREND等趋势分析算法将不适用。断点年份发生在2008年之前的区域集中在关中平原西安市周边以及秦巴山区的安康市南部。需要说明的是,本文研究区间虽然是在2001—2018年,但是检测的断点年份却在2004—2016年,当绿度变化发生在研究时间段的前3年或者后3年期间时,TSS-RESTREND算法无法检测出断点年份,这是由于当生态系统结构状态的改变需要几年的适应周期,反映在NDVI上,显示出一定的滞后性[36]。

表3 绿度分级统计

图7 断点年份检测Fig.7 Break year detection in Shaanxi from 2001 to 2018

4 结论

本文分析了自退耕还林生态工程实施以来的2001—2018年,陕西省NDVI的时空变化规律,利用TSS-RESTREND算法剥离了气候要素(降雨)对植被NDVI的影响,得到了人类活动对植被变化的影响程度和区域,主要结论如下：

(1) 2001—2018年间,陕西省NDVI呈显著增加,全省平均增加速率为0.006/a;

(2) 相比18年来的平均值,77.29%的区域大于均值。其中,陕北的榆林、延安大于均值的区域较大,分别为97.52%和89.03%,秦巴山区次之,为73.91%。2012年之后,NDVI高值向北逐年推进趋势明显。

(3)在只考虑降雨对植被的影响前提下,人类活动对植被的变化影响巨大,全省NDVI增加的区域达71.77%,而陕北地区的增加量明显大于关中平原区和陕南秦巴山地,其中陕北的榆林NDVI增加区域为72.11%,延安为86.44%,均超过了全省平均水平。

(4)总体上陕西全省呈变绿趋势。榆林市和延安市的变绿区域明显多于关中平原和秦巴山地,延安和榆林的剧烈增长区域分别为55.46%和34.34%,而陕南为41.03%,说明处于湿润气候区的陕南地区也有显著变绿趋势。