曾 群，殷守敬，田礼乔，陈晓玲

（1.华中师范大学学报编辑部，武汉430079；2.华中师范大学城市与环境科学学院，武汉430079；3.环境保护部卫星环境应用中心，北京100094；4.武汉大学测绘遥感信息工程国家重点实验室，武汉430079）

通过遥感技术对地球表面的长时间、连续观测，已经形成了海量的时序遥感数据，为全球地表过程和环境变化研究提供了重要数据源.归一化差分植被指数（NDVI）作为一个重要遥感参数，针对NOAA／AVHRR（GIMMS、Pathfinder）、MODIS、SPOT VEGETATION等传感器都有专门的NDVI产品，形成了长达12～21年的数据积累，在全球和区域环境变化、地表覆盖动态监测、植被生长过程分析等领域得到广泛应用［1－2］.

大区域覆盖、长时间跨度的NDVI数据集由于受到成像天气条件、气候变化、传感器衰减等因素的影响包含了大量噪声，尽管数据集经过了最大值合成、云检测等处理，但是数据产品仍然不同程度上受到噪声的影响，使得对数据的分析可能会导致错误的结论［2］，大大限制了其在全球变化定量分析研究中的应用.因此，对NDVI时间序列数据集进行去噪处理，剔除噪声对影像质量的影响，显得尤其重要.

针对NDVI时序数据去噪问题，目前已发展出许多算法，对各类算法优缺点也有过系统研究［3－5］.MVC（Maximum Value Composition）法［6］操作简单但残差保留过多，主要应用于NDVI数据集初级产品的生产和对定量化分析要求不高的研究；BISE（Best Index Slope Extraction）法［7］主要用于去除NDVI突降值，但阈值设置受人为影响较大；非线性高斯函数拟合法［8］适用于周期及趋势分析，但算法复杂不适合大数据量运算，异常值较多时会造成局部变化趋势错判；傅里叶变换［9］会对NDVI曲线过平滑造成变化信息的丢失，并且对NDVI数据中的伪极值比较敏感造成结果产生偏移；S－G（Savitzky－Golay）滤波［10］能够清晰描述NDVI长期变化趋势，同时保留局部的突变信息，该方法关键因素在于滑动窗口的大小和平滑多项式的阶数，其取值需要对NDVI数据分析后确定.

NDVI作为一种典型的地理数据，具备了地物时空相关的基本特性.现有的NDVI时序数据重建方法，往往只注重了其时间连续性，只对时域噪声进行了处理，而忽略了空间域噪声的去除.本文从时空相关性出发，提出一种综合利用S－G滤波、异常值检测和空间域中值滤波，适用于长时间序列NDVI数据的噪声去除方法.

1 基于时空相关的NDVI时序重建方法

1.1 基于S－G滤波的时域噪声去除

通过对前述NDVI时域噪声处理方法的对比，从NDVI数据的趋势变化分析和定量化研究等需求出发，综合考虑各种方法的噪声去除效果、数据保真性、参数自动化程度、计算效率等因素，选取S－G滤波方法进行时域噪声的去除.

基于S－G滤波的NDVI时序数据处理方法从提出［11］后被广泛应用［12－13］.该方法能够利用NDVI长期变化趋势（半年或一年）有效修正异常值，

随着迭代次数的增加，Fk呈抛物线形状变化，先是逐渐减小达到最小值而后逐渐增大.因此，Fk达到最小值的判断条件可以定义如下：

1.2 年序列异常值滤波

部分区域，如南部高海拔山区，由于长期云覆盖等条件影响在NDVI序列中连续出现较多的异常值，使得拟合出来的NDVI变化趋势偏差较大，S－G滤波后也难以达到理想的结果.对此，本文提出年际异常值滤波法，利用相邻年份同一时相的NDVI数据对S－G滤波后的值进行异常值判断和调整.具体如下：

如果待检测点，第t年NDVI值NDVIt如果满足：

即，比前1年和后1年的NDVI值NDVIt－1和NDVIt＋1都低，并且超过阈值，则判定为异常值，用前后年的均值代替，T1、T2表示判定阈值，表示NDVI的年际正常波动幅度，是一个大于零的小数，一般取0.2.

1.3 基于超限中值平滑滤波的空间域噪声去除

空间域噪声去除方法主要是借助地物的空间相关性，根据“距离越近的地物其相似性也较高”这一原理，利用邻域像素信息对影像进行平滑以去除噪声点.一般的空间平滑滤波在去噪的同时，往往造成边界模糊和信息失真.对于变化检测等定量分析研究，类型过渡的边界区域往往是研究的重点，在数据预处理中，要在最大程度上保留像素的原始值.基于该需求，结合现有的空间滤波方法，提出超限中值平滑滤波方法（Anomaly based Median Smoothing Filter，AMSF），即只是对超限像素进行中值平滑处理.

由于地物分布的空间异质性，用固定的阈值判定像素值是否超限是不合理的.在此我们借助基于统计的异常点检测方法判断像素是否超限.方法如下：

设有数组｛x｝，首先计算｛d（x）｝＝｛｜x0－m（x）｜，…，｜x1－m（x）｜，…，｜xi－m（x）｜，…，｜xn－m（x）｜｝，m（x）为｛x｝的中值，并计算数组｛d（x）｝的中值Md和中值绝对偏差MAD：MAD＝1.4826·Md.

然后计算异常度L：

一般情况下，如果有L＞3，认为该点为异常点.

中值滤波具有较好的椒盐噪声去除效果，并且具有较好的保边缘效果.在此，对异常点取窗口像素中值替代中心像素值.这样既可以有效抑制原始影像中的噪声，同时可以保证中心像素具备较好的同质性和连续性.

2 数据准备与实验

2.1 数据准备

选取江西省1998年4月1日～2009年4月1日共11年的SPOT VGT NDVI旬最大值合成产品S10.该数据以8位整型存储，在产品中自带了云标记产品，描述了成像质量及各像素成像状况（晴空、阴影或云等）.

对数据的预处理步骤包括：1）通过线性变换转化为反射率值；2）利用云标记产品，识别出云覆盖或阴影区域；3）利用前后时相无云数据线性插值，公式如下：同时可以体现NDVI的短期波动规律.为避免出现局部拟合不足的同时其它部分局部拟合过度的问题，对长时间序列NDVI数据，本文中以一年为一个周期进行分段拟合.具体步骤如下：

1）长期变化趋势拟合：拟合多项式表示NDVI的长期变化趋势，取较大的窗口宽度和较小的多项式系数可以获取平滑的拟合曲线.对不同的NDVI数据类型，需要根据数据时间分辨率等因素对参数进行设置.

2）NDVI序列初始化：低于长期变化趋势的数据点判定为噪音.利用长期变化趋势线拟合出的新值，取代NDVI序列中的噪声点，从而产生新的NDVI时序数列.新曲线会接近NDVI数据的上包络线.

3）NDVI序列迭代平滑：在迭代过程中，仍然利用S－G滤波来拟合长期变化趋势.但是重复迭代时，窗口宽度取值相对较小，多项式系数取值相对较大，以尽量体现出短期的变化趋势.令第k次拟合的拟合效果系数为Fk：NDV＊t表示时间t处插值后的NDVI值，Δt1、Δt2表示该像素点上后一时相和前一时相中第一个无云日期相隔天数.

2.2 时序重建

时序重建技术路线设计如图1所示.

图1 NDVI时序数据重建技术路线Fig.1 The technical route NDVI time series

2.2.1 S－G滤波 首先利用S－G滤波，对预处理后的NDVI数据以1年为一个时间段分段进行序列重建.S－G滤波参数设置时，滑动窗口越大，被平滑的峰谷值也就越多；平滑多项式次数越低，滤波结果会越平滑，但是也有可能会保留异常值，次数越高，就越容易去掉异常值，但是也有可能会因为拟合过度而出现更多的噪声.

对长期变化趋势模拟窗口半径宽度取9，即取6个月的开窗大小，多项式次数取2；短期变化趋势模拟窗口半径宽度取5，即开窗大小为3 mon，多项式次数取4.

2.2.2 年际异常滤波 对连续长时间受云覆盖等因素影响的区域，会产生连续异常值无法滤除.因此，利用历年同一时相的NDVI数据组成年序列影像（如1998年～2008年4月1日影像组成长度为11的序列），对每个像素点利用异常值检测，并对异常值利用公式（5）基于年序列数组进行线性插值补值.2.2.3空间平滑滤波 利用AMSF空间滤波（窗口大小取5×5）对数据逐波段平滑.

3 结果与讨论

如图2为山区林地区域样本点预处理前后NDVI曲线.对云标记已经标出的3个噪声数据点A、M、N，通过线性插值进行了修复，但是对于明显的噪声B点，未能检测出来.并且由于临近噪声点B，A点的修复效果也明显受到影响.

图2 样本点NDVI时序数据云插值预处理结果Fig.2 The cloud interpolation of sample points’NDVI time series

如图3所示，是对样本点预处理后的NDVI序列先后利用S－G滤波和年序列异常滤波方法处理后的结果.

从图3可以看出，通过S－G滤波，对数据进行了很好的平滑作用.不仅对噪声点A、B，原NDVI曲线的突降值（如E点）等大多数噪声点也通过多项式拟合检测出来，并进行了修复，修复后的曲线更加连续，能更真实地反映植被变化情况.

利用年序列异常滤波，对C、D时间段内的低异常值成功检测出来并进行了修复.5～6月是林地植被茂盛季，该时间段内出现了连续低值并且在之后NDVI值又恢复到峰值，很有可能是由于连续的多云天气造成的（可以排除发生地表覆盖类型变化、病虫害、火灾等情况），连续异常值的出现使得利用S－G滤波局部趋势拟合时出现偏差，未能有效修复.而通过年际同时相影像的异常分析，该时间段内数据均被检测为异常低值，并且通过插值达到了对其进行有效修复的目的.

图3 S－G滤波和年序列异常滤波效果图Fig.3 The effect picture of Savitzky－Golay filter and anomaly removal algorithm

如图4所示，经过前述时域处理的影像，通过在空间域中的分析，可以看出数据中存在空间离散点噪声.利用超限中值滤波，可以将空间上不连续的异常值噪声检测出来并修复.从图中可以看出，滤波后的影像中大部分区域像素值没有改变，仅对散点噪声进行了有效滤除，同时很好地保留了边缘部分信息.

图4 超限像素中值平滑滤波效果图左图为滤波前NDVI，右图为滤波后NDVIFig.4 The effect picture of median smoothing filter of overrun pixel

4 结语

本文从地理数据的时空相关性出发，建立了一种长时间序列NDVI数据重建方法，并利用11年时间序列的SPOT VGT NDVI旬合成数据进行了实验.首先利用分段S－G滤波对长时间序列NDVI数据时域噪声去除；基于年际异常信息判断，对长时间云覆盖等因素造成的连续异常值进行修复；然后利用超限中值滤波滤除了空间域离散噪声.结果表明，该方法可以有效地去除时空域噪声，为全球变化定量分析研究提供高质量数据集.

［1］ Wang L，Chen J，Gong P，et al.Land cover change detection with a cross－correlogram spectral matching algorithm［J］.International Journal of Remote Sensing，2009，30（12）：3259－3273.

［2］ Bradley B A，Jacob R W，Hermance J F，et al.A curve fitting procedure to derive inter－annual phonologies from time series of noisy satellite NDVI data［J］.Remote Sensing of Environment，2007，106：137－145.

［3］ Jennifer N H，Gregory J M.Noise reduction of NDVI time series：An empirical comparison of selected techniques［J］.Remote Sensing of Environment，2009，113：248–258.

［4］ 李 儒，张 霞，刘 波，等.遥感时间序列数据滤波重建算法发展综述.遥感学报，2009，13（2）：335－341.

［5］ 顾 娟，李 新，黄春林.NDV I时间序列数据集重建方法述评［J］.遥感技术与应用，2006，21（4）：391－395.

［6］ Holben B N.Characteristic of maximum value composite images for temporal AVHRR data［J］.International Journal of Remote Sensing，1986，7（11）：1417－1434.

［7］ Viovy N，Arino O，Belward A S.The best index slope extraction（BISE）：a method for reducing noise in NDVI timeseries［J］.International Journal of Remote Sensing，1992，13（8）：1585－1590.

［8］ Jensson P，Eklundh L.Seasonality extraction by function fitting to time－Series of satellite sensor data［J］.IEEE Transaction on Geoscience and Remote Sensing，2002，40（8）：1824－1832.

［9］ Roerink G，Menenti M，Verhoef W.Reconstructing cloud free NDVI composites using Fourier analysis of time series［J］.International Journal of Remote Sensing，2000，21（9）：1911－1917.

［10］ Savitzky A，Golay M J E.Smoothing and differentiation of data by simplified least squares procedures［J］.Analytical Chemistry，1964，36：1627－1639.

［11］ Chen J，Jonsson P，Tamura M，et al.2004.A simple method for reconstructing a high－quality NDVI time－series data set based on the Savitzky－Golay filter［J］.Remote Sensing of Environment，91（34）：332－344.

［12］ 边金虎，李爱农，宋孟强，等.MODIS植被指数时间序列Savitzky－Golay滤波算法重构［J］.遥感学报，2010，14（4）：725－741.

［13］ 宫 攀，陈仲新，唐华俊，等.基于MODIS温度／植被指数的东北地区土地覆盖分类［J］.农业工程学报，2006，22（9）：94－99.