从NDVI序列检测植被“绿化”或“褐化”变化趋势的新方法

唐亮, 赵忠明, 唐娉

(1.海南热带海洋学院,三亚 572022; 2.中国科学院遥感与数字地球研究所,北京 100101)

0 引言

卫星数据是当今陆地生态系统监测的主要信息源。几十a的卫星时间序列为定量描述全球植被活动及其变化趋势提供了数据基础。植被作为陆地生物圈的主要组成部分，是气候变化中的关键因素[1]，而且人们认为气候变暖将强烈影响陆地生物圈[2]。植被状态常用于评估自然和农用土地的生产率和退化[3-5]，其变化趋势被认为与土地退化关联[6-8]。陆地生态系统的时间动态包括渐变、突然变化和没有变化。没有变化是很少见的。在几十a的时间框架内，一般认为突然的植被减少，可能由一些短期的过程导致，如火灾、作物收割或发生了其他灾害； 而降雨事件或雪盖的减少也会引起突然的植被变绿。植被渐变会发生在更长的时间内，体现植被对全球变化的适应过程，包括大洋振荡，或持续的气候变化，如年际降水减少或大气中二氧化碳浓度增加。整个序列的渐变是最不复杂的变化类型，这种变化也指单调的“绿化(greening)”或“褐化(browning)”。归一化植被指数(normalized difference vegetation index，NDVI)与植被生产率强相关[9-10]，NDVI的趋势可近似表达植被“绿化”或“褐化”的趋势[11-12]。对于长时间序列，植被活动及其趋势分析方法已经有了一些研究[13-14]，其中最常使用的长时间序列趋势分析方法包括参数化和非参数化方法，前者有最小二乘的一元回归变化斜率法； 后者包括Sen趋势度估计、Mann-Kendall(MK)方法及季节性MK方法等。Beurs等[13]指出，基于最小二乘回归得到的NDVI曲线趋势的有效性会受到影响，因为生物物理参数的时间序列在时间上有相关性，会导致最小二乘回归方法应用的假设条件遭到破坏，这些条件是： ①所有的函数值应当是彼此独立的； ②残差有随机性； ③残差是零均值的； ④残差的方差应当对所有时间点是一样的。

另一方面，Belle等[15]指出Sen趋势度估计、MK方法及季节性MK方法具有一个隐含的假设，即季节间趋势必须是一致的，否则有实例表明每个季节有明显的趋势，但整体统计却没有趋势。因此在趋势性分析时需同时对时间序列的周期性或季节性建模，以去除整体统计对季节变化趋势的掩盖。故而进行趋势性分析是当前时间序列趋势分析的主要内容之一。因地球绕太阳公转的周期性影响引起的气候、地下水等多种周期性影响及地形的影响，NDVI序列具有强烈的周期性或多种周期叠加效应，周期性影响是NDVI趋势性分析必须考虑的因素。时频分析是周期性分析的基本手段，其中Fourier分析和小波分析是最常使用的时频分析工具，但对时间序列进行Fourier分解和小波分解均尚不能直接得到时间序列的趋势分量。Verbesselt等[14]在Fourier分解基础上提出了季节—趋势分解模型，将具有周期性特征的时间序列曲线拟合分解为趋势项(一阶线性函数)、季节项(周期函数组合)及残差的叠加，从而将季节项和趋势项分离出来，直接得到分段表示的NDVI变化趋势，但因为是拟合方法，其趋势项的估计仍然采用了最小二乘估计。不同于小波分析，同样适合于非平稳信号的新的时频分析工具经验模态分解(empirical mode decomposition，EMD)，可将时间序列中不同频率和振幅的信息很好地分解出来，随着分解的分量中频率逐渐降低而振幅增大，可分解得到表征原序列整体趋势的分量，该分量是近似单调变化的序列，即通常所说的线性趋势项，具有非常明确的物理意义。

本文以2006—2015年间中国局部地区的NDVI数据为数据源，探求EMD方法结合MK显著性检验进行植被单调性变化分析的方法。实验证明该方法是一种有效的时间序列变化趋势分析方法。

1 NDVI变化趋势分析方法

1.1 EMD方法

EMD方法是Huang等[16]于1998年提出的一种时频信号分析方法，属于自适应局部时频分析方法，它能根据数据信号本身的特性将其分解为有限个本征模态函数(intrinsic mode function，IMF)的叠加，所分解出来的各IMF分量包含了原信号不同时间尺度的局部特征信号，特别适合非平稳信号的分析[17]。目前，EMD方法已经广泛应用于天气、地震和医疗等信号处理的各个领域。其基本思想是将一个不规则信号表示成若干IMF与单调残差函数相加的形式，一维信号x(t)的分解可表示为

(1)

式中：imfi(t)为第i个IMF；n为IMF的个数；rn(t)为残差函数。

EMD方法的本质是通过数据的时间尺度特征来获得本征波动模式，进行数据分解。可以形象地称这种分解过程为“筛选(sifting)”过程，用x(t)代表原始数据信号，其分解步骤如下：

1)求取原始数据信号的极大值点和极小值点。

2)分别由极大值点和极小值点利用三次样条插值函数拟合获得上包络线s1和下包络线s2。

3)计算上、下包络线的均值m1，m1=(s1+s2)/2，如图1所示。

图1 原始NDVI时间序列和上、下包络线及其均值

4)将原始数据序列x(t)减去包络线均值m1，得到一个新的数据序列h0，h0=x(t)-m1。

5)将新的序列重复步骤1)—4)进行迭代，直到满足迭代停止准则，获得第一个IMF，imf1=hk。迭代停止准则计算公式为

(2)

当sdT超过给定的门限阈值时，迭代停止，此时hk代表第k次迭代后的新序列。

图2为第1次迭代后的结果。

图2 第1次迭代后的上、下包络线及其均值

从图2中可以看出，包络线起伏较大，但总体均值趋于0值。图3为第4次迭代后的结果。

图3 第4次迭代后的上、下包络线及其均值

从图3中可以看出，包络线均值起伏不大，已基本上趋于0值。

6)从原始数据序列中减去获取的第一个IMF，得到残差r1，r1=x(t)-imf1。

7)如果r1不是单调函数，重复步骤1)—6)，获取imf2，imf3，…，直到r1为单调函数结束。

EMD方法分解结果如图4所示。从图4可以看出，通过EMD方法首先得到的是原始数据信号中的高频分量imf1，后续的imf2—imf6分量频率逐渐降低，残差为单调函数，表示平均趋势。由于数据长度的有限性，残差中的极值点数量逐渐减少，使得原始数据信号经过几次提取后可以得到有限个IMF分量和最后的残差分量RES。

图4 EMD方法分解结果

一般趋势可以用残差表示，也可以用残差叠加若干IMF得到。图5所示的曲线是根据中国植被区划(1980)图[18]选取的2条曲线及其EMD方法分解的残差曲线RES和imf6+RES的结果，一条曲线位置在广东省南部的热带季雨林、雨林区域； 另一条在黄土高原的暖温带草原区域。

(a) 热带季雨林、雨林区域 (b) 黄土高原暖温带草原区域

从图5可以看出，如果只是把残差RES看做平均趋势项，2006—2015年间NDVI相对呈现单调变化趋势。如果把imf6+RES看做趋势项，2006—2008年间NDVI会呈现分段变化，图5中2区域都呈先升后降的现象，图5(a)NDVI在2006—2012年间先升后降，在2012—2015年间缓慢上升； 图5(b)NDVI在2006—2014年间一直呈现上升趋势，在2014—2015年间呈下降态势，而平均趋势一个下降一个则是上升。

图5中的NDVI上升或下降趋势是目视判断的，是否呈现显著的单调增、单调降趋势还是没有显著变化，需要采用显著性检验方法进行检验。

1.2 MK检验方法

由于本文主要针对单调性的变化趋势进行研究，因此主要针对残余量，即平均趋势进行趋势检验，检验是显著单调增、单调降还是没有显著变化。检验方法采用MK趋势检验方法。MK检验方法最初由Mann[19]和Kendall[20]提出，用来检测水域降水的长期变化趋势和突变情况。在时间序列趋势分析中，许多学者应用MK检验方法分析气温、水质、降水和径流等要素时间序列趋势变化[21-24]。MK检验方法不受少数异常值的干扰，也不需要样本遵循一定的分布，适用于非正态分布的数据。在MK检验中，设一时间序列数据(x1，x2,...,xn)是n个独立的、随机变量同分布的样本； 检验的统计量S计算公式为

(3)

(4)

当n>10时，S近似正态分布，其均值为0，方差var(S)=n(n-1)(2n+5)/18，标准的正态系统变量ZMK计算公式为

(5)

这样，在双边的趋势检验中，在α显著水平上，如果|ZMK|>Z1-α/2，则原假设是不可接受的，即在α显著水平上，时间序列数据存在明显上升或下降趋势，否则接受原始时间序列无变化趋势的假设。对于统计量ZMK，当ZMK>0时为上升趋势； 当ZMK<0时为下降趋势。当显著性水平为α时，置信度为(1-α)×100%，本文要求95%的置信度，则ZMK=1.96。

2 实验与分析

2.1 实验区域及实验数据

为了检验本文方法的鲁棒性，根据气候和纬度的不同，在中国北部、西部和南部的区域任意固定大小的3个样区，通过实验结果与样区实际情况的对比，分析评价该方法的适用性。样区A位于榆林附近，是黄土高原与内蒙古高原的过渡区。北部是毛乌素沙漠南缘风沙草滩区，南部是黄土高原腹地，沟壑纵横，丘陵峁梁交错，海拔为1 000～1 200 m。植被分区属于温带草原区域和温带荒漠区域。气温四季变化明显。经过长期退耕还林还草的治理，该区域如今沙滩地植被茂盛； 海子(湖泊)星罗棋布。水土侵蚀逐步得到治理，水土流失得到初步控制，生态环境有较大改善。样区B位于川、青、藏三省区的结合部，大陆性季风高原气候，气温低，年均气温在0 ℃以下，日照时间长，昼夜温差大，年均降水量为400～700 mm，海拔在4 000 m以上，植被分区属于青藏高原高寒植被区域，包括山地寒温性针叶林区域及高寒灌丛草甸地带。样区C位于广州附近，属于东亚季风区，从北向南分别为中亚热带、南亚热带和热带气候，是中国光、热和水资源最丰富的地区之一。该地区年平均降水量在1 300～2 500 mm之间，降水充沛，空间分布基本上也呈南高北低的趋势。改革开放后经济发展迅速。

3个样区实验数据选取的是2006—2015年MODIS 16 d合成的NDVI数据。从http: //ladsweb.nascom.nasa.gov/data/search.html网站下载MODIS产品数据MOD13A2数据，即全球1 km空间分辨率植被指数16 d合成数据。MOD13A2数据16 d的组合，1 a有23个时相，共10 a的数据。

实验数据的预处理包括重投影、数据重排和数据重构。重投影即将正弦投影转换为经纬度投影； 数据重排即将按空间顺序组织的数据按照时间顺序进行了重排； 数据重构采用Chen等[25]的方法,核心是通过Savitzky-Golay滤波器的迭代滤波使NDVI不断逼近上包络来减弱云噪声和大气变化的影响。

2.2 结果与分析

每个样区的卫星影像及经过趋势提取后的NDVI中值和NDVI趋势如图6。图6中(c)，(f)和(i)为3样区NDVI趋势图，深灰表示10 a间植被没有显著变化，浅灰表示NDVI值小于0.1，即植被很少，橘色表示10 a间植被有变褐的趋势，绿色表示10 a间植被有变绿的趋势。由图6(a)—(c)可知，样区A植被从北到南NDVI从低到高分布，即该地区植被的覆盖程度从南向北变差，这与该地区的地貌地形和气候是相吻合的，但从变化趋势看，从南到北，植被总体上是呈现变绿的趋势。这种现象与这10 a间退耕还林、水土侵蚀得到治理、水土流失得到控制、生态环境改善相吻合。由图6(d)—(f)可知，样区B植被的分布情况相对复杂，因地形复杂，有常年积雪区域存在，相应植被变化的趋势也比较复杂。在NDVI值比较低的周围既有呈现“褐化”的趋势，也有呈现“绿化”的趋势，产生这种现象的原因需要后续进行更深入细致的分析。由图6(g)—(i)可知，样区C除城市区域外，植被覆盖度很高。这与该地区的气候吻合。从植被变化的趋势看，存在“褐化”的区域也存在“绿化”的区域。“褐化”的区域大多与城市扩展、经济开发相关联。

(d) 样区B 卫星影像 (e) 样区B NDVI中值 (f) 样区B NDVI趋势图

(g) 样区C 卫星影像 (h) 样区C NDVI中值 (i) 样区C NDVI趋势图

3 结论

本文提出了一种利用EMD方法结合 MK检验方法从NDVI序列检测植被“绿化”或“褐化”变化趋势的新方法，该方法在检测变化趋势时无需去除序列的周期性或季节性直接进行趋势分析，对序列不需要作任何隐含假设，且趋势检测过程中不需要使用最小二乘进行一元回归。本文通过对中国北部、西部和南部3个实验区域2006—2015年10 a间植被变化趋势的分析，表明本文方法是一种有效的时间序列变化趋势的分析方法，对生态长期变化研究有一定参考价值。

本文只研究了10 a间NDVI变化的总体平均趋势，事实上10 a间植被的变化不一定是单调变化的，可能会先上升后下降，也可能先下降后上升，也可能分段上升或分段下降，有些变化可以直接关联到变化的驱动力，有些变化还不能直接关联到驱动力，而仅凭NDVI单要素的时间序列变化趋势分析生态驱动机理尚有不确定性。后续的研究一方面将扩展到对全国的植被“绿化”、“褐化”的趋势进行研究； 另一方面将结合数据的时间尺度特征通过EMD方法分解的本征模态的叠加，分析NDVI在10 a间的波动模式，期望对生态变化的时间尺度特征有更清晰的了解，从而评价中国10 a间退耕还林、退耕还草所取得的实效，为陆地生态系统变化研究提供支撑。