ABSTRACT

Long time series remote sensing data set is widely applied to global or regional environmental change, vegetation dynamic change, land cover change and plant biophysical parameters inversion. It is affected by clouds, aerosols, solar elevation and other factors, which make the data has a lot of noise. This paper is based on the MODIS-NDVI data, the Heihe River Basin, in the MATLAB programming environment, using time series fitting software TIMESAT3.1, using Savitzky-Golay (S-G) filter, asymmetric Gauss function (AG) and dual logic curve fitting (DL) method for the reconstruction of NDVI time series data of Heihe River Basin in 2012. The results show that the TIMESAT3.1 software can effectively denoise the NDVI time series, and the three methods have achieved good results on the time series curve. However, the unsymmetrical Gauss function and the double logical curve fitting results have the overfitting problem of "warping" in the beginning of the growing season of vegetation. In contrast, the S-G filter has a good fitting effect. The data source provides a good basis for the vegetation detection in the Heihe River Basin.

Keywords:Time Series Data, MODIS NDVI, Denoising, TIMESAT3.1, Heihe Basin

黑河流域时间序列植被指数曲线的平滑去噪对比研究

和萍，高书鹏

云南师范大学旅游与地理科学学院，云南 昆明

收稿日期：2018年2月1日；录用日期：2018年2月15日；发布日期：2018年2月23日

摘 要

长时间序列遥感数据集被广泛应用于全球或区域环境变化、植被动态变化、土地覆盖变化和植物生物物理参数反演等研究，受云、气溶胶、太阳高度角等因素的影响使数据存在很多的噪声。影响了数据分析和应用的效果。本文以MODIS-NDVI数据为基础，黑河流域为研究区，在MATLAB编程环境下，利用时间序列拟合软件TIMESAT3.1，使用Savitzky-Golay (S-G)滤波、非对称高斯函数(AG)和双逻辑曲线拟合(DL)方法对黑河流域2012年NDVI时间序列数据进行了重构。结果表明，TIMESAT3.1软件能有效的对NDVI时间序列进行平滑去噪，三种方法在时间序列曲线上均取得了较好的效果。但是不对称高斯函数和双逻辑曲线拟合结果在植被的生长季起始阶段均存在“翘起”的过度拟合问题。相对而言，S-G滤波拟合效果较好，NDVI曲线在保持原有基本形状的基础上更加有效的揭示所蕴含的物候周期性变化规律，该数据源对黑河流域的植被检测提供了良好的基础。

关键词 :时间序列数据，MODIS NDVI，去噪，TIMESAT3.1，黑河流域

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

陆地生态系统中的植被在全球物质与能量循坏中起着重要作用，植被通常是遥感观测和记录的第一表层，它影响着能量、生态、气候、水文、地貌、土壤等各种地理环境要素，是遥感工作者重点关注的对象之一 [1] [2] [3] 。NDVI曲线是NDVI时间序列数据构成的反映植被生物学特征相随时间变化的最佳指示因子，也是季节变化和人为活动影响的重要指示器 [2] 。植被指数时间序列分析使得研究区域乃至全球范围的物候现象(如返青、生长区间、生长衰退期)成为可能 [4] [5] 。基于NOAA/AVHRR、SPOT/VEGETATION以及MODIS等高时间分辨率的卫星传感器得到的植被指数之间序列资料已经在植被动态、变化监测、宏观植被覆盖分类和植物生物物理参数反演方面得到了广泛的应用 [6] [7] 。理论上，由于植被冠层随时间变化幅度较小，该曲线是一条连续平滑的曲线 [8] 。然而，星载传感器在数据采集和处理过程中受到云覆盖(是最主要的影响因素)、大气干扰、双向反射(太阳入射角、传感器观测角)、数据传输错误等各种因素的干扰，使得植被指数时序曲线波动大并出现许多噪声 [9] 。计算得到的NDVI也会受到很大的影响，在NDVI曲线中总是会有明显的突升或突降 [10] ，使数据的质量和应用效果受到严重影响 [11] ，阻碍了对数据的进一步分析利用并可能导致错误的结论 [12] [13] 。因此，在进行各种趋势分析和信息提取之前，有必要先对时间序列植被指数进行去噪和平滑处理，即时序植被指数重构 [14] 。

尽管目前大多植被指数产品都采用了最大值合成法来部分消除一些干扰的因素 [15] [16] ，如PAL数据集、GIMMS NDVI数据集、SPOT GVT产品以及MODIS时序产品，但是这种方法只能在一定程度上减小由于大气条件的影响而引起的噪声，由其他因素引起的噪声仍然存在 [17] 。为此，国内外众多研究已经开发使用了多种时序植被指数去噪方法，以实现去云处理、噪声去除和重建平滑植被指数曲线，这些拟合方法可以在大体上分为三类：1) 时间域上的处理，包括最佳指数斜率提取法 [8] [18] [19] (BISE)、中值迭代滤波法(MIF) [10] 、及频率域上的处理，如傅里叶变化 [20] [21] ，不对称函数拟合法，如不对称高斯函数拟合、双逻辑曲线函数拟 [17] 合与基于加权最小二乘回归的S-G滤波(Savitzky-Golay)算法 [20] 。上述时序数据的去噪方法广泛应用在遥感时序数据的去噪处理中，但由于各种方法构造的方式和使用的数学模型不同，在方法使用和去噪结果也表现出差别。不同地区适用的去噪方法也不一样，宋春桥等(2011) [22] 对藏北地区2007~2009年MODIS数据进行去噪对比研究，得出的结论为AG与D-L拟合的结果均优于Savitzky-Golay滤波法；草云锋等(2010年) [23] 对辽宁省长白山自然保护区为研究区，进行去噪对比，发现非对称高斯算法(AG)对原始高质量数据保真性最高。由此看来不同的研究区适用的去噪方法是不一样的，对黑河流域地区适用的拟合算法的研究还鲜见。本文基于黑河流域2012年时序NDVI数据，用非对称高斯函数拟合、双Logistic曲线拟合和Savitzky-Golay滤波3种常见的算法对时序植被指数进行去噪对比，基于该方法重构的高质量MODIS NDVI时间序列数据，为进一步利用该数据进行黑河流域生态、环境监测提供了较高的基础，同时可为对黑河流域位研究区的学者在进行陆地生态环境各方面研究中去噪预处理的方法选择提供参考。

2. 研究区域与数据

黑河流域(97˚21E~101˚44E，98˚36N~100˚48N)发源于青藏高原东北部，全长821 km，流域总面积约14.29万km²，是我国西北地区第二大内陆流域，位于河西走廊中部，远离海洋，周围高山环绕，气候干燥，降水稀少而集中，日照充足，太阳辐射强烈，昼夜温差大。流域内垂直高度变化大，最高海拔5573米，最低海拔1131米，海拔高差达4442米，东西和南北气候差异显著，且上、中、下游的地形地貌特征各不相同。上游以高山、冰川、森林、草甸草原为主，植被覆盖度较高。中、下游分布着大面积戈壁和沙漠。研究区如图1所示。

3. 数据及预处理

本论文采用的数据是NASA数据中心(https://wist.echo.nasa.gov/api)免费下载的MODIS数据产品，8天合成的500 m分辨率MOD09A1地表反射率产品，时间范围为2012年1月到2012年12月，共有44期影像。该数据集是通过对原始的MODIS数据进行大气校正、辐射校正、几何校正和拉伸后合成的数据，使用MRT (MODIS Reprojection Tools)软件对数据进行格式转换和投影转换，数据投影为Albers 投影，地理坐标为WGS-84，并完成图像的拼接和裁剪。最后计算研究区的归一化植被指数：

$\text{NDVI}=\frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+{\rho }_{red}}$ (1)

式(1)中，NDVI为归一化植被指数， ${\rho }_{nir}$ 为近红外通道反射率，MODIS对应第2通道； ${\rho }_{red}$ 为红光通道反射率，对应MODIS第1通道。需要对计算出来的NDVI值进行质量控制，通过ENVI的命令行批处理获得研究区Albers标准投影Geotiff格式的NDVI数据以及相应的质量控制数据，去掉有云的数据。计算时间序列NDVI值。研究区得到的原始的时序NDVI曲线如图2所示，可以看到，曲线有强烈的上下波动。

4. 算法原理和实现

4.1. 算法模型

4.1.1. S-G滤波算法

最早于1964年由Savizky和Goaly [24] 提出，又称最小二乘法或数据平滑多项式滤波器 [8] ，可以有效抑制时间序列数据中所包含的随机涨落。这种方法重建的NDVI时间序列能够清晰描述序列的长期变

图1. 研究区地理位置

图2. 原始数据时序曲线

化趋势以及局部的突变信息，且不受时间、空间尺度和传感器的限制 [9] 。基于S-G滤波原理，NDVI时间序列数据的S-G滤波过程可由下式描述：

$Y_j^{*}={\displaystyle \sum _{i=-m}^{i=m}\frac{C_i{Y}_{i+j}}{2m+1}}$ (2)

式中， $Y_j^{*}$ 为合成序列数据， $Y_{j+i}$ 代表原始序列数据， $C_i$ 为滤波系数。N为滑动窗口所包含的数据点( $2m+1$ ) [25] [26] 。

4.1.2. 非对称高斯函数

非对称高斯滤波(AG)由Jönsson和Eklundh [27] 于2002年提出，该算法使用分段高斯函数组合来模拟植被生长季物候，一个组合代表依 $\alpha$ 次植被盛衰过程，最后通过平滑连接各高斯拟合曲线，实现时间序列重建 [23] 。其主要过程大致可以分为区间提取、局部拟合和整体连接三个步骤 [9] 。先提取原始时序数据曲线中的谷值和峰值，采用高斯函数分别拟合曲线的左右部分。局部拟合函数为：

$f\left(t\right)=f\left(t,c_1,c_2,a_1,\cdots ,a_5\right)=c_1+c_{2}g\left(t;a_1,\cdots ,a_5\right)$ (3)

其中 $g\left(t;a_1,\cdots ,a_5\right)$ 为高斯函数，式中 $c_1和c_2$ 控制曲线的基准和幅度， $a_1$ 决定峰值和谷值的位置， $a_4、a_5$ 和 $a_2、a_3$ 分别控制曲线左、右部分的宽度和陡峭度。整体拟合函数为：

$F\left(t\right)=\{\begin{array}{c}\alpha \left(t\right)f_L\left(t\right)+\left(1-\alpha \left(t\right)\right)f_L\left(t\right),t_L<t<t_C\\ \beta \left(t\right)f_C\left(t\right)+\left(1-\alpha \left(t\right)\right)f_R\left(t\right),t_L<t<t_R\end{array}$ (4)

其中 $g\left[t_L,t_R\right]$ 是时序数据中待拟合部分的变化区间， $f_L\left(t\right)$ 、 $f_C\left(t\right)$ 和 $f_R\left(t\right)$ 分别代表 $g\left[t_L,t_R\right]$ 区间内左边谷值、中间峰值及右边谷值对应的局部拟合函数， $\alpha \left(t\right)$ 和 $\beta \left(t\right)$ 为介于0和1之间的剪切系数。通过整体拟合函数将局部拟合函数连接起来是该方法的关键之一，这种从局部到整体的拟合策略避免了整体数据对局部拟合的干扰，拟合后的曲线更加接近真实情况。

4.1.3. 双逻辑函数滤波

双重逻辑函数滤波(Double Logistic, DL)是由Beck [28] 于2006年提出的一种新的算法，与非对称高斯滤波函数类似，双逻辑函数滤波也是一种局部拟合函数，它们具有相同的基本公式。Beck在其文章中指出，相比于基于傅立叶变换的滤波算法，双重逻辑函数滤波对植被生长季的估计更加准确，特别对高纬度地区此算法更加适用。基本方程式为：

$g\left(t;x_1,\cdots ,x_4\right)=\frac{1}{1+\exp \left(\frac{x_1-t}{x_2}\right)}-\frac{1}{1+\exp \left(\frac{x_3-t}{x_4}\right)}$ (5)

式中t为生长季内某一时期；𝑥₁, 𝑥₂, 𝑥₃, 𝑥₄为曲线形状与位置参数，𝑥₁决定了曲线左边拐点的位置，𝑥₂确定了左侧曲线的变化率；𝑥₃确定了曲线右边拐点位置，𝑥₄决定了右侧曲线的变化率。设置𝑥₁, 𝑥₂, 𝑥₃, 𝑥₄不同的取值则拟合出不同的曲线位置与形状。

4.2. 算法实现

重建过程在Timesat3.1 [29] 中完成(图3)，TIMESAT是一个软件包用于分析时间序列卫星传感器数据，能够调查的卫星的时间序列数据的季节性和它们与植被的动态特性，在时域保存关于短期和长期植被变化的重要信息 [29] 。软件集成了目前已被国内外学者广泛应用的非对称高斯函数拟合、双逻辑函数拟合以及Savizky-Golay滤波三种方法，并取得较好的效果 [30] [31] 。由FORTRAN语言编译为最初的版本，之后由MATLAB两种语言编译成相对高的版本，在功能上有了较大的提高 [17] 。由于TIMESAT软件在数据平滑处理中要求数据具备两个以上的完整的生长周期，本文数据仅为一个周期，人为的将一个生长季分为两个。图4是TIMESAT3.1参数设置窗口，在设置中将映像方式设置为1 = images files，将image file type设置为8 bit，窗口大小设置为10，取值范围为−1~1，噪声去阈值(spike)设置为2，拟合峰值(alptitude)参数为1，迭代一年的数据行有337，列有377。实验中多次尝试设置不同的参数，但是结果变化不大。

5. 结果与分析

图5中表示的是在MATLAB环境下，利用TIMESAT3.1软件用三种方法对黑河流域时序NDVI曲线实现去噪的结果。三种拟合方法得到的结果如图5的(a)、(b)、(c)所示，从处理的结果看，各种拟合方法均可在一定程度上消除部分噪声的影响，实现了NDVI时序数据的去噪，能有效的对NDVI时间序列进行平滑去噪，可使NDVI曲线在保持原有基本形状的基础上更加有效的揭示所蕴含的物候周期性变化规律，突出了MODIS时间序列的优点，并对物候参数提取有很大的帮助。分析结果，(a)和(b)表示的非对称高斯函数和双逻辑函数的拟合结果在植被的生长季起始阶段存在“翘起”的过度拟合问题，相对而

图3. TIMESAT界面

图4. TIMESAT3.1参数设置窗口

图5. 时间序列曲线去躁前后对比((a)：非对称高斯函数；(b)：双逻辑回归函数；(c)：双S-G滤波)

表1. 非对称高斯函数重构前后残差

表2. 双逻辑函数重构前后残差

表3. 双S-G滤波重构前后残差

言，S-G滤波拟合效果较好，更逼近原始数据，S-G滤波曲线变化趋势相对稳定。对S-G滤波的结果进行第二次S-G滤波，称为双S-G滤波。

如表1、表2、表3分别表示的是时间序列三种滤波方法在重构前后NDVI比较，从表示可以看到三种滤波方式中，双S-G滤波的残差是最小的，为0.00027，其次为双逻辑函数拟合法，残差为0.000723，最后为非对称高斯函数拟合法，残差为0.01975。

6. 结论

本文以MODIS-NDVI数据为基础，综合运用MRT (MODIS Reprojection Tools)、Timesat3.1和ENVI IDL软件中的信息提取、分析等技术，运用相关方法将黑河流域的时序植被指数进行平滑去噪的研究。研究内容概括为以下几个部分，以黑河流域为研究区域，采用MRT对数据进行数据重采样(包括投影转换以及格式转换)、数据提取、图像拼接，计算NDVI值。本文使用ENVI对数据质量控制，去掉有云的数据，进行空间插值。在TIMESAT3.1中用双逻辑回归函数、非对称高斯函数、S-G滤波的三种方法对黑河流域NDVI曲线平滑去噪的实现，对结果进行分析，反映植被变化趋势。结果表明，三种方法在一定程度上均实现了去噪过程，但在黑河流域地区，S-G滤波的方法更适合用于数据重建，非对称高斯函数和双逻辑回归函数均存在起点拟合过度的问题，起点处存在“翘起”过度拟合的现象。双S-G滤波去噪能够有效的去噪，提高了数据质量。不同的地区适用的去噪方法因区域、气候、植被等原因适用的去噪方法均不同的。因此，分为不同植被类型不同区域进行比较研究有待进一步尝试。

文章引用

和 萍,高书鹏. 黑河流域时间序列植被指数曲线的平滑去噪对比研究
Study on the Smooth Denoising Contrast of the Vegetation Index Curve of the Time Series of the Heihe River Basin[J]. 地球科学前沿, 2018, 08(01): 32-41. http://dx.doi.org/10.12677/AG.2018.81004

参考文献 (References)