CN101650422A - 遥感植被指数时间序列数据去噪方法 - Google Patents

Abstract

本发明提供了一种遥感植被指数时间序列数据去噪方法，该方法包括：根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各等距离的数据点的数值并以此生成标准数据，其中，距离为遥感植被指数时间序列数据的合成周期；在标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；根据特征点数据对标准数据按照预定滤波规则进行滤波处理得到滤波数据，如果滤波数据满足预定判决条件，则将滤波数据作为去噪后的遥感植被指数时间序列数据。根据本发明，能够有效地去除遥感植被指数时间序列数据中的噪声，能够保留植被迅速变化的特征，并且处理速度快、效率高。

Description

遥感植被指数时间序列数据去噪方法

技术领域

本发明涉及遥感图像处理技术，具体地，涉及一种遥感植被指数时间序列数据去噪方法。

背景技术

遥感植被指数(Vegetation Index，简称为VI)时间序列数据记录了地表植被(以下简称为植被)的动态变化。但是受到云污染、大气变异等因素的影响，遥感植被指数时间序列数据往往存在大量噪声，这一噪声限制了该数据的进一步应用。目前，相关技术中已发展了很多方法来除去这些噪声(即遥感植被指数时间序列数据去噪方法)，例如：最佳指数斜率提取法(the Best Index Slope Extractionmethod，简称为BISE)、修正的BISE滤波方法(the Modified BISEFiltering)、傅里叶分析法(Fourier Analysis)、加权最小二乘法(theWeighted Least-Squares approach)、多项式最小二乘拟合法(thePolynomial Least Squares Operation method，简称为PoLeS)、地学统计法(Geostatistics)、非对称高斯函数拟合法(the AsymmetricGaussian Function Fitting)、双逻辑斯蒂函数拟合法(the DoubleLogistic Function-Fitting)、SG滤波法(Savitzky-Golay Filtering)、均值迭代滤波法(Mean-Value Iteration Filter，简称为MVI)等。

上述遥感植被指数时间序列数据去噪方法的工作原理一般分两步：首先，根据一定的策略进行噪声点识别；然后，去除识别出的噪声点。以下详细说明上述过程。

(一)识别噪声点

遥感植被指数时间序列数据的变化反映了植被的渐变过程，该变化表现出一定的趋势，如果该时间序列数据中的一些数据点偏离了该变化趋势(例如，一个遥感植被指数的值远高于或低于其周围的数据点值)，则认为这些数据点是噪声点。也即，遥感植被指数时间序列数据去噪方法都存在一个预设前提：遥感植被指数时间序列数据是一个缓慢变化的过程，而那些迅速变化的点被认为是噪声。

相关技术通常采取三种方法识别遥感植被指数时间序列数据的变化趋势、或者在识别变化趋势的同时判断噪声点：第一种，使用大滤波窗口对该时间序列数据进行平滑滤波或函数拟合，从而得到一条平滑的趋势时间序列，如SG滤波方法；第二种，根据一定的决策判断方法识别噪声点，即根据植被的生长规律来设置判断条件，如果数据点的值不满足该判断条件，则认为该数据点为噪声，例如：植被处于出苗至成熟期这段时间，其植被指数的值应该一直处于升高的趋势，如果在该时间段内某一植被指数值反而下降，则认为它是噪声点，如BISE方法；第三种，频率域滤波法，该方法认为遥感植被指数时间序列数据的高频部分为噪声数据，通过去除高频数据来去除噪声，即先将遥感植被指数时间序列数据转换到频率域，然后将高频部分(即噪声)抑制，将低频部分再反变换到空间域，从而得到遥感植被指数时间序列趋势，如傅里叶分析法，该方法在识别噪声点的同时去除噪声点。

(二)去除噪声点

通常采取两种方法去除噪声点：第一种是迭代滤波法，该方法用噪声点两边最近邻点的植被指数值的平均值来替换该噪声点的植被指数值，如均值迭代滤波法；第二种是最小二乘拟合法，该方法首先根据一定规则对遥感植被指数时间序列进行分段，然后对各段进行函数或多项式最小二乘拟合，将拟合的结果替换原来的遥感植被指数值，例如加权最小二乘法、多项式最小二乘拟合法、非对称高斯函数拟合法、双逻辑斯蒂函数拟合法、SG滤波法等。

但是，在相关技术中，遥感植被指数时间序列数据去噪方法至少存在以下三个问题。

第一，无法保留植被迅速变化的特征。相关技术中的去噪方法主要针对自然植被区，其预设前提均认为植被指数是一个缓慢变化的过程，因此那些迅速变化的点应该被当作噪声而去除。但是，对人工植被区(例如一年两季或一年三季的农业植被区)或受强烈干扰的自然植被区(例如发生大规模火灾或病虫害的天然森林)而言，这些快速变化的点则恰恰是地物变化本身的反映，它们作为遥感植被指数时间序列上的特征点，往往对进一步的应用具有重要作用，因此应该将他们保留而不是去除。

第二，去噪后的数据引入了新的误差。现有的遥感植被指数时间序列数据，其空间分辨率均比较低(一般在250m以上)，每个像元均可能包含多种不同物候的植被，因此各像元所反映的植被指数时间序列往往是不规划的曲线，很难用特定的数学函数进行刻画，所以相关技术中上述基于函数的去噪方法，如非对称高斯函数拟合法、双逻辑斯蒂函数拟合法等都会在去噪的时候引进一些新的错误数据；同样，基于傅里叶分析的方法也会存在类似问题，因为它们也是基于规则的正弦或余弦函数。

第三，去噪处理速度慢、效率低。上述大部分去噪方法由于存在最小二乘法拟合过程，里面涉及大量的循环迭代运算，降低了去噪处理的速度和效率，难以满足应用的要求。

由此可见，相关技术中遥感植被指数时间序列数据去噪方法存在无法保留植被迅速变化的特征、去噪后引入新的误差、去噪处理速度慢效率低这三方面的问题。

发明内容

有鉴于此，本发明提出了一种遥感植被指数时间序列数据去噪方法，用于解决相关技术中遥感植被指数时间序列数据去噪方法存在的无法保留植被迅速变化的特征、去噪后引入新的误差、去噪处理速度慢效率低的问题至少之一。

根据本发明的一个方面，提供了一种遥感植被指数时间序列数据去噪方法。

根据本发明的遥感植被指数时间序列数据去噪方法包括：根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各等距离的数据点的数值并以此生成标准数据，其中，距离为遥感植被指数时间序列数据的合成周期；在标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；根据特征点数据对标准数据按照预定滤波规则进行滤波处理得到滤波数据，如果滤波数据满足预定判决条件，则将滤波数据作为去噪后的遥感植被指数时间序列数据。

优选地，上述对时间轴上各等距离数据点进行插值的处理包括：根据质量评估文件，在遥感植被指数时间序列数据中选出各合成周期内质量最优的数据点，或者选出各合成周期内质量等级低于预定质量等级阈值的数据点，并以此生成最优数据；根据最优数据、合成日期文件依次对时间轴上各等距离数据点按照预定插值规则X_i计算其数值，其中， $X_i=\frac{Y_i-Y_{i-1}}{T_{y_i}-T_{y_{i-1}}}*(T_{x_i}-T_{y_{i-1}})+Y_{i-1},$ X_i为时间轴上各等距离数据点x_i的数值，Y_i为最优数据中数据点y_i的数值，

为合成日期文件中记录的y_i的合成日期，

为x_i在时间轴上的日期，i＝2，3，…，n，n为自然数。

优选地，上述预定查询规则包括：设置查询窗口 $f=\frac{r}{s},$ 其中，f为查询窗口的长度，r为地表植被生长周期的时长，s为合成周期，生长周期为植被从播种到收割的生长过程；对于分别以标准数据中的每个数据点为起点的多个查询窗口中，如果查询窗口中的极值数据点位于查询窗口的中间位置，则该极值数据点为极值点。

优选地，上述预定查询规则还包括：对于标准数据中第一个极值点之前的所有数据点、和最后一个极值点之后的所有数据点，如果其中的一个数据点和与该数据点最邻近的极值点的距离为至少半个查询窗口的长度、并且该数据点的数值和与该数据点最邻近的极值点的数值之间的差值大于第一预定阈值fet₁，则该数据点为极值点，并将查找到的极值点按照其合成日期的先后顺序存入特征点数据中；其中，fet₁小于或等于D_i，D_i与生长周期对应的遥感植被指数时间序列数据中的最小波峰的最大值和最小值的差值。

优选地，上述预定查询规则还包括：第一去伪处理：如果特征点数据中任意相邻的两个极值点均为最大值点，则删除其中数值较小的极值点，如果该相邻的两个极值点均为最小值点，则删除其中数值较大的极值点；第二去伪处理：如果特征点数据中任意相邻的两个极值点的差值小于第二预定阈值fet₂，则删除该相邻的两个极值点中的第一个极值点，并重新进行第一去伪处理；其中，fet₂小于或等于D_i。

优选地，上述查询窗口f取与

最接近的奇数为其值。

优选地fet₂大于或等于fet₁。

优选地，上述滤波处理包括：将标准数据作为第一数据，使用滤波器对第一数据进行卷积滤波并生成第二数据，其中，滤波器为 $(\frac{1}{k+2},\frac{k}{k+2},\frac{1}{k+2}),k=\mathrm{1,2,3},...,m,$ k为当前进行滤波处理的次数，m为自然数；用特征点数据中的所有数据点替换第二数据中相应位置的数据点，并将替换后的第二数据作为滤波数据；上述预定判决条件包括： $D_j^{k+1}{<\mathrm{fet}}_3,$ 如果滤波数据不满足预定判决条件，则将滤波数据作为第一数据，并对该第一数据进行滤波处理，直至滤波数据满足预定判决条件，其中， $D_j^{k+1}=|N_j^{K+1}-N_j^k|,$ N_j ^k是第k次滤波处理后滤波数据中数据点j的数值，D_j ^k+1是相邻两次滤波处理生成的两个滤波数据中相应位置的数据点之间的差值，fet₃为第三预定阈值，j＝1，2，3，…，n，k＝1，2，3，…，m，n为滤波数据中数据点的个数，n和m均为自然数；其中，fet₃小于或等于D_i，D_i为与生长周期对应的遥感植被指数时间序列数据中的最小波峰的最大值和最小值的差值。

优选地，上述预定判决条件还包括：对第一数据进行滤波处理的次数大于10次。

优选地，fet₃小于fet₁。

借助于本发明提供的技术方案，通过对时间轴上各等距离(即合成周期)的数据点插值，将遥感植被指数时间序列数据还原为等距离的标准数据，在标准数据中查找出所有极值点，并根据所有的极值点对标准数据进行滤波处理，能够有效地去除遥感植被指数时间序列数据中的噪声，并且能够保留遥感植被指数时间序列数据中植被迅速变化的特征。

在本发明的优选方案中，根据预定的查询规则能够在标准数据中可靠、全面、准确地查找出反映植被变化的极值点，并且根据全部极值点对标准数据进行滤波处理，能够有效、可靠地保留遥感植被指数时间序列数据中植被迅速变化的特征。

在本发明的优选方案中，根据预定的滤波器对数据进行卷积滤波，该滤波器不通过数学函数表达，不会在滤波处理的同时引入新的误差，并且该滤波器不进行最小二乘法拟合，从而能够加快处理速度快、提高处理效率。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是根据本发明实施例的遥感植被指数时间序列数据去噪方法的流程图；

图2是步骤S102的具体处理流程；

图3是步骤S104的具体处理流程；

图4是步骤S106的具体处理流程；

图5a为应用于本发明实施例的实验区中森林区域的NDVI时间序列数据；

图5b为应用于本发明实施例的实验区中一年两季作物冬油菜-夏水稻区域的NDVI时间序列数据；

图5c为应用于本发明实施例的实验区中一年两季作物冬小麦-夏水稻区域的NDVI时间序列数据；

图6是根据本发明实施例的遥感植被指数时间序列数据去噪方法具体应用的处理流程；

图7a是对试验区中森林区域的NDVI时间序列数据的处理效果对比图；

图7b是对试验区中一年两季作物冬油菜-夏水稻的NDVI时间序列数据的处理效果对比图；

图7c是对试验区中一年两季作物冬小麦-夏水稻的NDVI时间序列数据的处理效果对比图；

图8a是试验区的原始NDVI图像；

图8b是根据本发明实施例的遥感植被指数时间序列数据去噪方法的对试验区NDVI图像处理的效果图。

具体实施方式

以下结合附图对本发明的实施例进行说明，应当理解，此处所描述的实施例仅用于说明和解释本发明，并不用于限定本发明。

图1是根据本发明实施例的遥感植被指数时间序列数据去噪方法的流程图，如图1所示，根据本发明实施例的遥感植被指数时间序列数据去噪方法包括：

步骤S102，根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各等距离的数据点的数值并以此生成标准数据，其中，所述距离为遥感植被指数时间序列数据的合成周期；

步骤S104，在标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；

步骤S106，根据特征点数据对标准数据按照预定滤波规则进行滤波得到滤波数据，如果滤波数据满足预定判决条件，则将滤波数据作为去噪后的遥感植被指数时间序列数据。

下面对上述处理过程进行详细说明。

(一)步骤S102

在相关技术中，遥感植被指数时间序列数据产品往往包含质量评估文件和合成日期文件。其中，质量评估文件对遥感植被指数数据的每一个数据点都给出一个质量等级，合成日期文件对遥感植被指数数据的每一个数据点都给出了对应的合成日期。由于遥感植被指数时间序列数据在前期的预处理过程中采用了最大值合成技术，此技术会使最终的遥感植被指数数据的多个数据点在时间轴上表现为非等间隔(或非等距离)的数据点，而在后续的线性插值处理中会默认地将它们视为等距离的数据点来进行处理，这将会使遥感植被指数时间序列数据产生忽高忽低的抖动，这样很容易将这些抖动的数据点当作噪声来处理，而实际情况是这些数据点本身是正常的，只是因为在时间轴上被当作等间隔处理时造成了伪噪声。

例如：如果采用合成周期为16天的最大值合成遥感植被指数数据，对于时间序列上两个相邻的数据点，第一个数据点的合成时间是1月1日，第二个数据点的合成时间是1月31日，二者之间实际相差30天，而在默认情况下系统就会按16天来处理这两个数据，当它们的植被指数值保持不变时，其波形的斜率(植被指数值之差除以合成时间之差)就会明显增大，因为在分子(即植被指数值之差)不变的情况下分母(即合成时间之差)由30天减小到了16天，所以将这两个数据点在等间隔的时间序列上处理时就会表现为植被指数值突然升高的假象，这就造成了伪噪声，这样不能反映植被的真实变化，并且还会因为伪噪声造成新的误差。

针对这种问题，本发明实施例提供的遥感植被指数时间序列数据去噪方法对遥感植被指数时间序列数据进行预处理，将其还原为时间轴上多个等距离的数据点，这样在后续的处理过程中能够基于较为准确的数据进行去噪处理(即滤波处理)，从而能够避免上述类型的伪噪声以及由该种伪噪声带来的误差。

图2示出了步骤S102的具体处理流程，如图所示，步骤S102的具体处理过程包括：

步骤S1022，根据质量评估文件，在遥感植被指数时间序列数据中选出各合成周期内质量最优的数据点，或者选出各合成周期内质量等级低于预定质量等级阈值的数据点，并以此生成最优数据；需要说明的是，遥感植被指数的质量等级越高表示数据质量越差，质量等级越低表示数据质量越高；

步骤S1024，根据最优数据、合成日期文件依次对时间轴上各等距离数据点按照预定插值规则X_i计算其数值，其中， $X_i=\frac{Y_i-Y_{i-1}}{T_{y_i}-T_{y_{i-1}}}*(T_{x_i}-T_{y_{i-1}})+Y_{i-1},$ X_i为时间轴上各等距离数据点x_i的数值，Y_i为最优数据中数据点y_i的数值，为合成日期文件中记录的y_i的合成日期，

为x_i在时间轴上的日期，i＝1，2，3，…，n，n为自然数，上述距离为遥感植被指数时间序列数据的合成周期；

步骤S1026，根据上述各等距离的数据点生成标准数据。

通过上述处理，能够将遥感植被指数时间序列数据还原为等距离的标准数据，这样避免将遥感植被指数时间序列数据在时间轴上等间隔处理时造成的伪噪声。

(二)步骤S104

为了解决相关技术中无法保留植被迅速变化的特征的问题，本发明实施例将遥感植被指数时间序列数据中的反映植被迅速变化特征的极值点提取出来、生成特征点数据，并在后续处理中保留这些极值点，从而能够保留植被迅速变化的特征。

图3示出了步骤S104的具体处理流程，如图3所示，步骤S104的具体处理过程包括：

步骤S1042，设定查询窗口 $f=\frac{r}{s},$ 其中，f为查询窗口的长度，r为地表植被生长周期的时长，s为合成周期，该生长周期为植被从播种到收割的生长过程；优选地，查询窗口f取与

最接近的奇数为其值；

例如：中国北方的冬小麦一般是在10月底播种，在次年的3月份开始返青，6月上旬开始收割，这期间遥感植被指数时间序列数据会显示出两个波峰，一个波峰对应于播种至返青期，另一个则对应于返青期至收割期，因此冬小麦单个生长周期所覆盖的时长约为110天，所以查询窗口也应该覆盖110天左右(即r＝110)，对于合成周期为16天(即s＝16)的遥感植被指数数据来说，查询窗口f的大小应该为7；

步骤S1044，对于分别以标准数据中的每个数据点为起点的多个查询窗口中，如果查询窗口中的极值数据点位于查询窗口的中间位置，则该极值数据点为极值点，将查找到的所有极值点生成特征点数据；

另一方面，由于标准数据中的第一个查询窗口和最后一个查询窗口中的部分数据中也可能存在极值点，而这两部分数据无法用步骤S1044依次轮询到，这样需要对它们进行进一步的识别：

步骤S1046，对于标准数据中第一个极值点之前的所有数据点、和最后一个极值点之后的所有数据点，如果其中的一个数据点与和该数据点最邻近的极值点的距离为至少半个查询窗口的长度、并且该数据点的数值和与该数据点最邻近的极值点的数值之间的差值大于阈值fet₁，则该数据点为极值点，将查找到的极值点按照合成时间的先后顺序存入上述特征点数据，优选地，设置fet₁≤D_i，D_i为与生长周期对应的遥感植被指数时间序列数据中的最小波峰的最大值和最小值的差值；如上所述，植被生长周期的遥感植被指数时间序列数据可能会包括至少一个生长波峰；

在上述查找到的所有极值点中，可能会由于噪声而存在伪极值点，这些噪声可能是由于极端的大气环境造成的，例如降雪、长时间的阴雨天气，这样就需要根据极值点的特征和极值点间的关系将伪极值点去除；一方面，由于植被生长周期中的最大值和最小值是相间排列的，如果相邻的两个极值点均为最大值或最小值，则应该去除其中之一：

步骤S1048，依次判断特征点数据中任意相邻的两个极值点是否为最大值点和最小值点，如果均为最大值点，则处理进行到步骤S1050，如果均为最小值点，则处理进行到步骤S1052，否则，处理进行到步骤S1054；

步骤S1050，删除两个极值点中数值较小的极值点，处理进行到步骤S1048；

步骤S1052，删除两个极值点中数值较大的极值点，处理进行到步骤S1048；

另一方面，由于极大值和极小值代表了植被生长过程中两个截然不同的物候状况，这样任意相邻的两个极大值和极小值之间的差值应该大于一定的阈值，如果它们不满足这一条件，则这两个极值点中至少有一个数据点可能是伪数据点，优选地，可以设置fet₂≤D_i；

步骤S1054，依次判断特征点数据中任意相邻的两个所述极值点的差值是否小于阈值fet₂，如果判断为是，处理进行到步骤S1056，否则，处理进行到下述步骤S1060；优选地，还可以设置fet₁≤fet₂，并且设置fet₂接近于D_i，这样可以更为严格地筛除伪极值点；

步骤S1056，删除这两个极值点中的第一个极值点，处理进行到S1048；

通过上述处理，将标准数据中的极值点查找出来并生成特征点数据，能够记录植被迅速变化特征，为后续处理做准备。

(三)步骤S106

经过了上述步骤S102和步骤S104的预处理，可以基于较为准确的标准数据和特征点数据来进行滤波处理(即去噪处理)。

图4示出了步骤S106的具体处理流程，如图4所示，步骤S106的具体处理过程包括滤波处理和判决处理：

1、滤波处理

步骤S1060，将标准数据作为第一数据，并且设置k＝1，k为进行滤波处理的次数；

步骤S1062，使用滤波器

(或称为3点变权重滤波器)对第一数据进行卷积滤波，将滤波后的数据作为第二数据；

当k＝1时，上述滤波器为均值滤波，随着滤波处理次数的增加，k的值也越来越大，相应地，滤波器的中间值

的权重也越来越大，而两边相邻点的权重则越来越小；本发明实施例提供的这种变权重滤波策略，每次仅采用3个数据点进行变权重滤波，滤波幅度小，可以很好地保持原始遥感植被指数时间序列数据的曲线形状，并且计算量小、计算速度快，还可以促进算法收敛，从而可以提高运算效率；

但是，随着卷积滤波次数的增加，极值点的植被指数值将会被逐渐改变，这会导致最后的结果产生物候偏移，并且削弱遥感植被指数时间序列数据的形状特征，因此，极值点的数值应该在滤波数据中予以保留：

步骤S1064，用特征点数据中的所有数据点替换第二数据中相应位置的数据点，并将替换后的第二数据作为滤波数据；

经过了滤波处理(即步骤S1062的卷积滤波和步骤S1064的特征点替换)后的滤波数据中可能仍然会存在噪声，这种噪声通常是由于多云天气(云量大于10％)时的大气环境造成的，在滤波数据中这种噪声表现为幅度较小的数据抖动，如果这种数据抖动不满足滤波要求，可以对滤波数据进行再次滤波处理：

2、判决处理

步骤S1066，如果滤波数据满足预定判决条件则处理进行到步骤S1068，否则，处理进行到步骤S1070，该预定判决条件包括： $D_j^{k+1}<{\mathrm{fet}}_3,$ 其中， $D_j^{k+1}=|N_j^{K+1}-N_j^k|,$ N_j ^k是第k次滤波处理后滤波数据中数据点j的数值，D_j ^k+1是滤波数据和第一数据(或者为前后相邻两次滤波处理生成的两个滤波数据)中相应位置的数据点之间的差值，fet₃为阈值，j＝1，2，3，…，n，k＝1，2，3，…，m，n为滤波数据中数据点的个数，n和m均为自然数；其中，可以设置fet₃≤D_i；

滤波处理后产生的滤波数据会比第一数据更为光滑，这就意味着滤波数据(即新产生的遥感植被指数时间序列数据)含有更少的噪声；此外，滤波数据和第一数据之间的差异在有噪声的地方会比在没有噪声的地方更明显，如果该差异在一定的范围内时，即 $D_j^{k+1}<{\mathrm{fet}}_3,$ 则该噪声是可以接受的；即处理可以进行到步骤S1068，否则对该滤波数据继续进行滤波处理，即处理进行到步骤S1070；

此外，优选地，还可以设置0＜fet₃＜fet₁，fet₃的值可以根据对滤波曲线的具体要求来确定，fet₃的值越小，即对噪声的容忍度越小，经过滤波处理后的曲线越平滑，但是，需要进行的滤波处理的次数会越多，这样算法收敛会越慢，即处理速度会相对较慢、处理效率相对较低；与此相反，fet₃的值越大，即对噪声的容忍度越大，处理后的曲线会相对较粗糙，但算法收敛较快，即处理速度相对较快、处理效率相对较高；

步骤S1068，将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。

步骤S1070，如果k＞10，即进行滤波处理的次数大于10次，则处理进行到步骤S1072，否则，处理进行到步骤S1074；设置滤波处理次数的上限可以防止对极端异常的数据点进行滤波处理的死循环；

步骤S1072，将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。

步骤S1074，将滤波数据作为第一数据，并且，k＝k+1，处理转入步骤S1062。

通过上述处理，能够在对标准数据进行滤波的同时，保留植被迅速变化的特征，并且本发明实施例提供的滤波器能够有效地保持遥感植被指数时间序列数据的曲线形状，并且运算速度快、效率高，从而能够解决相关技术中遥感植被指数时间序列数据去噪方法存在的无法保留植被迅速变化的特征、去噪后引入新的误差、去噪处理速度慢效率低的问题。

下面对本发明实施例的具体应用进行说明。

应用本发明实施例将江苏省江宁区作为试验区，使用本实施例提供的方法对江宁区的MODIS(Moderate Resolution ImagingSpectrometer)归一化差值植被指数(即NDVI)进行处理。该NDVI的合成周期为16天，空间分辨率为250米，NDVI时间序列为2006年1月1日至2008年12月31日，3年共69个数据点(每年23个数据点)。图5a、图5b、图5c分别为应用于本发明实施例的实验区中3个像元点的NDVI时间序列数据，图5a为森林区域的NDVI时间序列数据，图5b为一年两季作物冬油菜-夏水稻的NDVI时间序列数据，图5c为一年两季作物冬小麦-夏水稻的NDVI时间序列数据。在图5a、图5b、图5c中，横轴原点时间为2006年1月1日，横轴的单位为天，纵轴为NDVI的数值；图中阴影区域的NDVI为2007年的时间序列数据。

图6示出了根据本发明实施例的遥感植被指数时间序列数据去噪方法具体应用的处理流程。

表一中示出了下列处理过程中涉及到的参数(质量等级阈值fet_q、fet₁、fet₂、fet₃)。

表一

如图6所示，根据本发明实施例的遥感植被指数时间序列数据去噪方法具体应用的处理流程包括如下步骤：

步骤S201，根据质量评估文件，从遥感植被指数时间序列数据中选出质量等级低于质量等级阈值的所有数据点，并将这些数据点作为最优数据，即Q_i＜fet_q，fet_q＝4，其中，Q_i为质量评估文件中记录的NDVI中数据点y_i的质量等级，fet_q为预定的质量等级阈值；

步骤S202，根据插值规则X_i对时间轴上各等距离数据点进行插值，并将这些等距离数据点作为标准数据，其中， $X_i=\frac{Y_i-Y_{i-1}}{T_{y_i}-T_{y_{i-1}}}*(T_{x_i}-T_{y_{i-1}})+Y_{i-1},$ X_i为时间轴上各等距离数据点x_i的数值，Y_i为最优数据中数据点y_i的数值，

为合成日期文件中记录的y_i的合成日期，

为x_i在时间轴上的日期，i＝1，2，3，…，n，n为自然数；

步骤S203，设置查询窗口 $f=\frac{r}{s},$ 对于该实验区的NDVI，主要农作物的最短生长波峰长度在r≈110天，对于合成周期s＝16，查询窗口 $f=\frac{r}{s}=\frac{110}{16}\&ap;7;$

步骤S204，依次在以标准数据中的每个数据点为起点的多个查询窗口中查找符合极值点(即特征点)判断条件的数据点，并将找到的所有极值点作为特征数据，该极值点判断条件为查询窗口中的极值数据点位于该查询窗口的中间位置；

步骤S205，查询标准数据两端可能存在的极值点，对于标准数据中第一个极值点之前的所有数据点、和最后一个极值点之后的所有数据点，如果其中的数据点与和该数据点最邻近的极值点的距离至少大于3.5(即至少半个查询窗口的长度)、并且该数据点的数值和与该数据点最邻近的极值点的数值之间的差值大于fet₁，则该数据点为极值点，将查找到的极值点按照合成时间的先后顺序存入上述特征点数据，其中，fet₁＝0.1；

步骤S206，判断特征点数据中任意相邻的两个极值点是否为最大值点和最小值点，如果均为最大值点，则删除其中数值较小的极值点，如果均为最小值点，则删除其中数值较大的极值点；

步骤S207a，如果特征点数据中任意相邻的两个极值点的差值小于阈值fet₂，处理进行到步骤S207b，否则，处理进行到步骤S208，其中，fet₂＝0.15；

步骤S207b，删除这两个极值点中的第一个极值点，并且，处理转入步骤S206；

步骤S208，将标准数据作为第一数据，并且令k＝1；

步骤S209，用滤波器对第一数据进行滤波，将滤波后的数据作为第二数据；

步骤S210，用特征点数据中的所有数据点替换第二数据中相应位置的数据点，并将替换后的第二数据作为滤波数据；

步骤S211，如果滤波数据中的所有数据点均满足预定判决条件： $D_j^{k+1}<{\mathrm{fet}}_3,$ 则处理进行到步骤S212，否则，处理进行到步骤S213，其中fet₃＝0.05；

步骤S212，将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。

步骤S213，判断滤波处理的次数是否大于10次，即k＞10，如果判断为是，则处理进行到步骤S214，否则，处理进行到步骤S215；

步骤S214，将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。

步骤S215，将滤波数据作为第一数据，并且，k＝k+1，处理进行到步骤S209。

经过上述处理后，能够在有效地去除试验区的遥感植被指数时间序列数据中存在的噪声，还能够有效地保存植被迅速变化的特征、提高处理速度和处理效率。下面对本发明实施例提供的方法的去噪效果进行详细说明。

(一)保存植被迅速变化的特征

图7a、图7b、图7c示出了本发明实施例提供的方法与相关技术的方法对试验区的遥感植被指数时间序列数据处理后波形的对比，相关技术中选择目前使用效果较好的非对称高斯函数拟合法(AG)、双逻辑斯蒂函数拟合法(DL)和SG滤波法，图7a、图7b、图7c显示了上述四种方法(即本发明实施例提供的方法和三种相关技术的方法)在实验区3个像元点上的去噪结果，这三个像元点分别为森林区域、一年两季作物冬油菜-夏水稻、一年两季作物冬小麦-夏水稻。在图7a、图7b、图7c中，横轴原点时间为2006年1月1日，横轴的单位为天，纵轴为NDVI的数值；图中阴影区域的NDVI为2007年的时间序列数据。

图7a显示了对试验区中森林区域的NDVI时间序列数据的处理效果对比，如图7a所示，自2006年至2008年，NDVI时间序列数据每年表现为一个波峰，四种方法的去噪结果类似，但本发明实施例提供的方法对2008年1月底至2月初的NDVI极小值保留得很好，NDVI在此时间会出现这种极小值，主要是2008年1月底至2月初中国南方经历了一次异常的冰冻天气过程，持续时间约一个月，大部分植被均遭到破坏，因此此时的NDVI值明显降低。

图7b显示了对试验区中一年两季作物冬油菜-夏水稻的NDVI时间序列数据的处理效果对比，如图7b所示，NDVI时间序列数据每年均会表现为两个波峰，前一个波峰对应于冬油菜，后一个波峰对应于夏水稻，而相关技术的三种去噪方法均把每年的两个波峰合成了一个没有实际意义的错误波峰，这样很容易让后续的用户将其判断为森林植被类型，而本发明实施例提供的新方法在除去噪声的同时却很好的保留了这两个波峰。

图7c显示了对试验区中一年两季作物冬小麦-夏水稻的NDVI时间序列数据的处理效果对比，图7c的效果与图7b的效果类似，这里不再赘述。

从图7a至图7c中可以看出，本发明实施例提供的方法对自然植被的去噪效果优于相关技术的上述三种方法，尤其是对农区植被来说，本发明实施例提供的方法能有效地保持遥感植被指数时间序列数据中的极值点(即特征点)。

图8b示出了根据本发明实施例的方法处理NDVI图像的效果图。选择实验区中一块大小为400×400像元的区域进行测试，图8a示出了试验区的原始NDVI图像，该图像为2008年6月上旬的一期图像。该实验区(江苏省江宁区)6月上旬还处于梅雨季节，NDVI的质量较差，造成NDVI值偏低并表现明显的马赛克现象，经过本发明的新方法进行滤波去噪后，如图8b所示，偏低的NDVI值基本得到恢复，而且马赛克现象也基本去除。

可见本发明实施例提供的方法能够有效地去除遥感植被指数时间序列数据中存在的噪声、并且能够可靠地保持植被迅速变化的特征。

(二)提高处理速度和处理效率

对于实验区的测试影像(400像元×400像元×69层)，本发明用交互数据语言(Interactive Data Language，简称为IDL)6.4编写的处理程序在桌面台式机(1个中央处理器即CPU，CPU主频为2.21GHz，内存为2GB随机存储器即RAM)上进行运算，所花费的时间是47s，根据这一时间推算，对于一景标准的MODIS影像(4800像元×4800像元)，如果时间序列为3年(69层)，其处理时间还不到2小时。对于SG方法来说，有研究表明(Chen et al.，2004)，SG方法处理一幅8849像元×5601像元×48层的图像，用桌面台式机(4个CPU，CPU主频为1.8GHz，内存为1GB RAM)花了22小时，如果折算成时间序列为3年的标准MODIS影像，它所花的时间约为14.7小时，这还没有考虑计算机性能上的差异。对于非对称高斯和双逻辑斯蒂函数拟合法来说，它们所花的时间比SG方法更长。

可见，本发明实施例提供的发明能够极大程度地提高遥感植被指数时间序列数据去噪的处理速度和处理效率。

综上所述，根据本发明实施例提供的技术方案，通过对时间轴上各等距离的数据点插值将遥感植被指数时间序列数据还原为等距离的标准数据，在标准数据中查找出所有极值点，并根据所有的极值点对标准数据进行滤波，能够有效地去除遥感植被指数时间序列数据中的噪声，并且能够保留遥感植被指数时间序列数据中植被迅速变化的特征；根据预定的查询规则能够在标准数据中可靠、全面、准确地查找出反映植被变化的极值点，并且根据全部极值点对标准数据进行滤波处理，能够有效、可靠地保留遥感植被指数时间序列数据中植被迅速变化的特征；根据预定的滤波器(即3点变权重滤波器)对数据进行卷积滤波，该滤波器不通过数学函数表达，不会在滤波处理的同时引入新的误差，并且该滤波器不进行最小二乘法拟合，从而能够加快处理速度快、提高处理效率。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种遥感植被指数时间序列数据去噪方法，其特征在于，包括：

根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各所述等距离的数据点的数值并以此生成标准数据，其中，所述距离为遥感植被指数时间序列数据的合成周期；

在所述标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；

根据所述特征点数据对所述标准数据按照预定滤波规则进行滤波处理得到滤波数据，如果所述滤波数据满足预定判决条件，则将所述滤波数据作为去噪后的遥感植被指数时间序列数据。

2.根据权利要求1所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述对时间轴上各等距离数据点进行插值包括：

根据所述质量评估文件，在所述遥感植被指数时间序列数据中选出各合成周期内质量最优的数据点，或者选出各合成周期内质量等级低于预定质量等级阈值的数据点，并以此生成最优数据；

根据所述最优数据、所述合成日期文件依次对所述时间轴

上各所述等距离数据点按照所述预定插值规则X_i计算其数值，

其中， $X_i=\frac{Y_i-Y_{i-1}}{T_{y_i}-T_{y_{i-1}}}*(T_{x_i}-T_{y_{i-1}})+Y_{i-1},$ X_i为所述时间轴上各所述等距离数据点x_i的数值，Y_i为所述最优数据中数据点y_i的数值，

为所述合成日期文件中记录的y_i的合成日期，

为x_i在所述时间轴上的日期，i＝1，2，3，…，n，n为自然数。

3.根据权利要求1所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述预定查询规则包括：

设置查询窗口 $f=\frac{r}{s},$ 其中，f为查询窗口的长度，r为地

表植被生长周期的时长，s为所述合成周期，所述生长周期为植被从播种到收割的生长过程；

对于分别以所述标准数据中的每个数据点为起点的多个所述查询窗口中，如果所述查询窗口中的极值数据点位于所述查询窗口的中间位置，则该极值数据点为所述极值点。

4.根据权利要求3所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述预定查询规则还包括：

对于所述标准数据中第一个所述极值点之前的所有数据点、和最后一个所述极值点之后的所有数据点，如果其中的一个数据点和与该数据点最邻近的所述极值点的距离为至少半个所述查询窗口的长度、并且该数据点的数值和与该数据点最邻近的所述极值点的数值之间的差值大于第一预定阈值fet₁，则该数据点为所述极值点，并将查找到的所述极值点按照其合成日期的先后顺序存入所述特征点数据中；

其中，所述fet₁小于或等于D_i，所述D_i为与所述生长周期对应的遥感植被指数时间序列数据中最小波峰的最大值和最小值的差值。

5.根据权利要求4所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述预定查询规则还包括：

第一去伪处理：如果所述特征点数据中任意相邻的两个所述极值点均为最大值点，则删除其中数值较小的极值点，如果该相邻的两个所述极值点均为最小值点，则删除其中数值较大的极值点；

第二去伪处理：如果所述特征点数据中任意相邻的两个所述极值点的差值小于第二预定阈值fet₂，则删除该相邻的两个所述极值点中的第一个所述极值点，并重新进行所述第一去伪处理；

其中，所述fet₂，小于或等于所述D_i。

6.根据权利要求3至5中任一项所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述查询窗口f取与最接近的奇数为其值。

7.根据权利要求5所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述fet₂大于或等于所述fet₁。

8.根据权利要求1至3中任一项所述的遥感植被指数时间序列数据去噪方法，其特征在于，

所述滤波处理包括：将所述标准数据作为第一数据，使用滤波器对所述第一数据进行卷积滤波并生成第二数据，其中，所述滤波器为 $(\frac{1}{k+2},\frac{k}{k+2},\frac{1}{k+2}),k=\mathrm{1,2,3},...,m,$ k为当前进行所述滤波处理的次数，m为自然数；用所述特征点数据中的所有数据点替换所述第二数据中相应位置的数据点，并将替换后的所述第二数据作为滤波数据；

所述预定判决条件包括： $D_j^{k+1}<{\mathrm{fet}}_3,$ 如果所述滤波数据不满足所述预定判决条件，则将所述滤波数据作为第一数据，并对该第一数据进行所述滤波处理，直至所述滤波数据满足所述预定判决条件，其中， $D_j^{k+1}=|N_j^{K+1}-N_j^k|,$ N_j ^k是第k次所述滤波处理后所述滤波数据中数据点j的数值，D_j ^k+1是相邻两次所述滤波处理生成的两个所述滤波数据中相应位置的数据点之间的差值，fet₃为第三预定阈值，j＝1，2，3，…，n，k＝1，2，3，…，m，n为所述滤波数据中数据点的个数，n和m均为自然数；

其中，所述fet₃小于或等于D_i，所述D_i为与所述生长周期对应的遥感植被指数时间序列数据中的最小波峰的最大值和最小值的差值。

9.根据权利要求8所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述预定判决条件还包括：

对所述第一数据进行所述滤波处理的次数大于10次。

10.根据权利要求8所述的遥感植被指数时间序列数据去噪方法，其特征在于，所述fet₃小于所述fet₁。

Images