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

专利名称：遥感植被指数时间序列数据去噪方法
技术领域：
本发明涉及遥感图像处理技术，具体地，涉及一种遥感植被指
凄t时间序列凝:据去噪方法。
背景技术：
遥感档L:帔指数(Vegetation Index,简称为VI)时间序列数据记 录了地表植被(以下简称为植被)的动态变化。但是受到云污染、 大气变异等因素的影响，遥感植被指数时间序列数据往往存在大量 噪声，这一噪声限制了该数据的进一步应用。目前，相关技术中已 发展了很多方法来除去这些噪声(即遥感植被指数时间序列数据去 p桑方法)，例如最佳指翁:杀牛率4是取法(the Best Index Slope Extraction method,筒称为BISE )、 f务正的BISE滤波方法(the Modified BISE Filtering )、 4專里叶分冲斤法(Fourier Analysis )、力口冲又最小二乘法(the ^Veighted Least-Squares approach )、 多项式最小二乘4以合法(the Polynomial Least Squares Operation method,简称为PoLeS )、地学统 i十法(Geostatistics )、 非只十^尔高#Jt函凄t4以合法(the Asymmetric Gaussian Function Fitting),只又逻4辱^j"蒂函凄史4以合》去(the Double Logistic Function-Fitting )、 SG滤波法(Savitzky—Golay Filtering )、 均值迭代滤波法(Mean-Value Iteration Filter,简称为MVI)等。
上述遥感相^皮指数时间序列凄t据去p喿方法的工作原理一4殳分两 步首先，根据一定的策略进行噪声点识别；然后，去除识别出的 f喿声点。以下详细"i兌明上述过程。(一) 识别，噪声点
遥感植被指数时间序列数据的变化反映了植被的渐变过程，该 变化表现出 一定的趋势，如果该时间序列数据中的 一些数据点偏离 了该变化趋势(例如， 一个遥感植一皮指数的值远高于或低于其周围 的数据点值)，则i人为这些数据点是噪声点。也即，遥感賴3皮指数时
间序列凄t据去p喿方法都存在一个预-没前提遥感档^皮指凄t时间序列 数据是一个緩慢变化的过程，而那些迅速变化的点被认为是噪声。
相关技术通常采取三种方法识别遥感植被指数时间序列数据的 变化趋势、或者在识别变化趋势的同时判断噪声点第一种，使用 大滤波窗口对该时间序列凄t才居进4亍平滑滤波或函凄t拟合，乂人而纟寻到 一条平滑的趋势时间序列，如SG滤波方法；第二种，根据一定的 决策判断方法识别噪声点，即根据植被的生长规律来设置判断条件， 如果数据点的值不满足该判断条件，则认为该数据点为噪声，例如 植被处于出苗至成熟期这段时间，其植被指数的值应该一直处于升 高的趋势，如果在该时间段内某一植被指数值反而下降，则认为它 是p喿声点，如BISE方法；第三种，频率i或滤波法，该方法i^为遥 感植被指数时间序列数据的高频部分为噪声数据，通过去除高频数 据来去除噪声，即先将遥感植被指数时间序列数据转换到频率域， 然后将高频部分(即噪声)抑制，将低频部分再反变换到空间域， 从而得到遥感植被指数时间序列趋势，如傅里叶分析法，该方法在 识别噪声点的同时去除噪声点。
(二) 去除噪声点
通常采取两种方法去除噪声点第一种是迭代滤波法，该方法 用噪声点两边最近邻点的植被指数值的平均值来替换该噪声点的植 被指数值，如均值迭代滤波法；第二种是最小二乘拟合法，该方法 首先根据一定规则对遥感植被指数时间序列进行分段，然后对各段进行函数或多项式最小二乘拟合，将拟合的结果替换原来的遥感植 被指数值，例如加^又最小二乘法、多项式最小二乘拟合法、非对称 高斯函数拟合法、双逻辑斯蒂函数拟合法、SG滤波法等。
但是，在相关4支术中，遥感核j皮指凄t时间序列数据去噪方法至 少存在以下三个问题。
第一，无法保留植被迅速变化的特征。相关技术中的去噪方法 主要针对自然植被区，其预设前提均认为植被指数是一个緩慢变化 的过程，因此那些迅速变化的点应该被当作噪声而去除。但是，对 人工植被区(例如一年两季或一年三季的农业才直被区)或受强烈干 扰的自然植被区(例如发生大规模火灾或病虫害的天然森林)而言， 这些快速变化的点则恰恰是地物变化本身的反映，它们作为遥感植
净皮指^:时间序列上的特4i点，往往对进一步的应用具有重要作用， 因此应该将他们^呆留而不是去除。
第二，去噪后的数据引入了新的误差。现有的遥感植被指数时 间序列数据，其空间分辨率均比较低(一般在250m以上)，每个像 元均可能包含多种不同物候的植被，因此各像元所反映的植被指数 时间序列往往是不失见划的曲线，4艮难用特定的数学函数进行刻画， 所以相关技术中上述基于函数的去噪方法，如非对称高斯函数拟合 法、双逻辑斯蒂函数拟合法等都会在去噪的时候引进一些新的错误 数据；同样，基于傅里叶分析的方法也会存在类似问题，因为它们 也是基于规则的正弦或余弦函数。
第三，去噪处理速度慢、效率低。上述大部分去噪方法由于存 在最小二乘法拟合过程，里面涉及大量的循环迭代运算，降{氐了去 噪处理的速度和效率，难以满足应用的要求。由此可见，相关技术中遥感植被指数时间序列数据去噪方法存 在无法保留植被迅速变化的特征、去噪后引入新的误差、去噪处理 速度慢效率低这三方面的问题。

发明内容
有鉴于此，本发明提出了一种遥感植被指数时间序列数据去噪 方法，用于解决相关技术中遥感植被指数时间序列数据去噪方法存 在的无法保留档^皮迅速变化的特征、去噪后引入新的误差、去噪处 理速度慢效率低的问题至少之一。
根据本发明的一个方面，提供了 一种遥感植被指数时间序列数 才居去噪方法。
才艮据本发明的遥感档」故指数时间序列凄t据去噪方法包括才艮据 遥感植被指数时间序列数椐、合成日期文件、和质量评估文件对时 间轴上各等距离的数据点按照预定插值规则进行插值，得到各等距 离的数据点的数值并以此生成标准数据，其中，距离为遥感植被指 数时间序列数据的合成周期；在标准数据中按照预定查询规则查找 出所有极值点并以此生成特征点数据；根据特征点数据对标准数据 按照预定滤波规则进行滤波处理得到滤波数据，如果滤波数据满足 预定判决条件，则将滤波数据作为去噪后的遥感植被指数时间序列 数据。
优选地，上述对时间轴上各等距离数据点进行插值的处理包括 根据质量评估文件，在遥感植被指数时间序列数据中选出各合成周 期内质量最优的数据点，或者选出各合成周期内质量等级低于预定 质量等级阈值的数据点，并以此生成最优数据；根据最优数据、合 成日期文件依次对时间轴上各等距离数据点按照预定插值规则z,计
算其数值，其中，z,^^*(t; -rv )+l， z,为时间轴上各等距离
Y, 一 乂-i数据点x,的数值，y;为最优数据中数据点y,.的数值，7;为合成日期文 件中记录的y,的合成日期，7;为x,在时间轴上的日期，^=1，2,3广'，"， n为自然凄t。
优选地，上述预定查询规则包括设置查询窗口/=二，其中，
/为查询窗口的长度，r为地表才直被生长周期的时长，s为合成周期，
生长周期为柏3皮从播种到收割的生长过程；对于分别以标准数据中 的每个数据点为起点的多个查询窗口中，如果查询窗口中的极值数 据点位于查询窗口的中间位置，则该极值数据点为极值点。
优选地，上述预定查询规则还包括对于标准数据中第一个极 值点之前的所有数据点、和最后一个极值点之后的所有数据点，如 果其中的一个数据点和与该数据点最邻近的极值点的距离为至少半 个查询窗口的长度、并且该数据点的数值和与该数据点最邻近的极 值点的数值之间的差值大于第 一预定阈值，则该数据点为极值 点，并将查找到的极值点按照其合成日期的先后顺序存入特征点数 据中；其中，/a小于或等于A， A.为与生长周期对应的遥感档^皮指 数时间序列邀:据中的最小波峰的最大值和最小值的差值。
Y尤选地，上述预定查询^见则还包4舌第一去伪处理如果特4正 点数据中任意相邻的两个极值点均为最大值点，则删除其中数值较 小的4及值点，如果该相邻的两个才及值点均为最小值点，则删除其中 数值较大的极值点；第二去伪处理如果特征点数据中任意相邻的 两个极值点的差值小于第二预定阈值/ef2,则删除该相邻的两个极值 点中的第一个极值点，并重新进行第一去伪处理；其中，/^2小于或 等于A。
优选地，上述查询窗口/耳又与二最接近的奇数为其值。j尤选:t也大于或等于优选地，上述滤波处理包括将标准lt据作为第一ft据， -使用 滤波器对第一数据进行巻积滤波并生成第二数据，其中，滤波器为(l,l,l)，A: -1，2，3广'，附，*为当前进4亍滤波处理的次凄t, m为自 A: + 2 "2 "2然数；用特征点数据中的所有数据点替换第二数据中相应位置的数 据点，并将替换后的第二数据作为滤波数据；上述预定判决条件包 括D)+1</"3,如果滤波数据不满足预定判决条件，则将滤波数据作为第一数据，并对该第一数据进行滤波处理，直至滤波数据满足预定判决条件，其中，"+1=|wf+1-《|,《是第A:次滤波处理后滤波 数据中数据点的数值，^"是相邻两次滤波处理生成的两个滤波凄t据中相应位置的数据点之间的差值，为第三预定阈值， /=1，2,3，"'," ， A:=l,2，3，"'，m ，"为^t^皮lt才居中凄t才居,存、^^个凄t,"和mi勾为自然数；其中，/^小于或等于A， A为与生长周期对应的遥感植 被指数时间序列数据中的最小波峰的最大值和最d 、值的差值。优选地，上述预定判决条件还包^l舌对第一数据进行滤波处理 的次lt大于10次。优选地，'J 、于/e"。借助于本发明提供的技术方案，通过对时间轴上各等距离(即 合成周期)的数据点插值，将遥感植-故指数时间序列数据还原为等 距离的标准数据，在标准数据中查找出所有极值点，并根据所有的 极值点对标准数据进行滤波处理，能够有效地去除遥感植被指数时 间序列数据中的噪声，并且能够保留遥感植被指数时间序列数据中 才直被迅速变化的特征。在本发明的优选方案中，根据预定的查询规则能够在标准数据 中可靠、全面、准确地查找出反映植被变化的极值点，并且根据全 部才及值点对标准:数据进4于滤波处理，能够有效、可靠地保留遥感植 被指数时间序列数据中植被迅速变化的特征。在本发明的优选方案中，根据预定的滤波器对数据进行巻积滤 波，该滤波器不通过数学函数表达，不会在滤波处理的同时引入新 的误差，并且该滤波器不进行最小二乘法拟合，从而能够加快处理 速度快、提高处理效率。本发明的其它特征和优点将在随后的说明书中阐述，并且，部 分地从说明书中变得显而易见，或者通过实施本发明而了解。本发 明的目的和其他优点可通过在所写的"i兌明书、权利要求书、以及附 图中所特别指出的结构来实现和获得。


附图用来提供对本发明的进一步理解，并且构成说明书的一部 分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中图1是根据本发明实施例的遥感植被指数时间序列数据去噪方 法的流程图；图2是步骤S102的具体处理流程；图3是步-骤S104的具体处理流程；图4是步骤S106的具体处理流程；序列数据；图5b为应用于本发明实施例的实验区中一年两季作物冬油菜 —夏水^舀区i或的NDVI时间序列数据；图5c为应用于本发明实施例的实验区中一年两季作物冬小麦 一夏水稻区域的NDVI时间序列凝:据；图6是才艮据本发明实施例的遥感植被指数时间序列数据去噪方 法具体应用的处理流程；图7a是对试验区中森林区域的NDVI时间序列数据的处理效果 对比图；图7b是对试验区中一年两季作物冬油菜一夏水稻的NDVI时 间序列凄t据的处理效果对比图；图7c是对试验区中一年两季作物冬小麦一夏水稻的NDVI时间 序列凝:据的处理效果对比图；图8a是试验区的原始NDVI图像；图8b是根据本发明实施例的遥感植被指数时间序列数据去噪 方法的对试-验区NDVI图〗象处理的效果图。
具体实施方式
以下结合附图对本发明的实施例进行:说明，应当理解，此处所 描述的实施例仅用于说明和解释本发明，并不用于限定本发明。图1是根据本发明实施例的遥感植被指数时间序列数据去噪方 法的流程图，如图1所示，根据本发明实施例的遥感植被指数时间 序列数据去噪方法包括
步骤S102,根据遥感植被指数时间序列数据、合成日期文件、 和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进 行插值，得到各等距离的数据点的数值并以此生成标准数据，其中， 所述距离为遥感才直被指数时间序列凝:据的合成周期；
步骤S104,在标准数据中按照预定查询规则查找出所有极值点 并以此生成特^正点数据；
步骤S106,根据特征点数据对标准数据按照预定滤波规则进行 滤波得到滤波凄t据，如果滤波数据满足预定判决条件，则将滤波数 据作为去噪后的遥感植被指数时间序列数据。
下面对上述处理过禾呈进4亍详细i兌明。
( 一 )步-骤S102
在相关技术中，遥感植被指数时间序列数据产品往往包含质量 评估文件和合成日期文件。其中，质量评估文件对遥感才直纟皮指数凄t 据的每一个数据点都给出 一个质量等级，合成日期文件对遥感植被
指数数据的每一个数据点都给出了对应的合成日期。由于遥感植被 指数时间序列凄t据在前期的预处理过程中采用了最大值合成技术， 此技术会使最终的遥感植被指数数据的多个数据点在时间轴上表现 为非等间隔(或非等距离)的数据点，而在后续的线性插值处理中 会默认地将它们视为等距离的数据点来进行处理，这将会使遥感植 被指数时间序列数据产生忽高忽低的抖动，这样很容易将这些抖动 的数据点当作噪声来处理，而实际情况是这些数据点本身是正常的， 只是因为在时间轴上被当作等间隔处理时造成了伪噪声。例如如果采用合成周期为16天的最大值合成遥感植被指数数 据，对于时间序列上两个相邻的数据点，第一个数据点的合成时间 是1月1日，第二个数据点的合成时间是1月31日，二者之间实际 相差30天，而在默i人情况下系统就会4要16天来处理这两个凄t据， 当它们的植被指数值保持不变时，其波形的斜率(植被指数值之差 除以合成时间之差)就会明显增大，因为在分子(即植被指数值之 差)不变的情况下分母(即合成时间之差)由30天减小到了 16天， 所以将这两个数据点在等间隔的时间序列上处理时就会表现为植被 指数值突然升高的假象，这就造成了伪噪声，这样不能反映植被的 真实变化，并且还会因为伪噪声造成新的误差。
针对这种问题，本发明实施例提供的遥感植被指数时间序列数 据去噪方法对遥感植被指数时间序列数据进行预处理，将其还原为 时间轴上多个等距离的数据点，这样在后续的处理过程中能够基于 较为准确的数据进行去噪处理(即滤波处理)，从而能够避免上述类 型的伪噪声以及由该种伪噪声带来的误差。
图2示出了步骤S102的具体处理流程，如图所示，步骤S102 的具体处理过程包括
步骤S1022,根据质量评估文件，在遥感植被指数时间序列数 据中选出各合成周期内质量最优的数据点，或者选出M成周期内 质量等级低于预定质量等级阈值的数据点，并以此生成最优数据； 需要说明的是，遥感植被指数的质量等级越高表示数据质量越差， 质量等级越低表示数据质量越高；
步-骤S1024,才艮据最优^t据、合成日期文件依次对时间轴上各 等距离数据点按照预定插值规则Z,.计算其数值，其中，
X, = "，1 *( ; -1;,》+ L,, Z,.为时间轴上各等距离数据点x,.的数值，^为最优数据中数据点y,.的数值，r),为合成日期文件中记录的y,的合 成日期，7;为;c,在时间轴上的日期，/=1,2，3广',"，"为自然数，上述
距离为遥感植被指数时间序列数据的合成周期；
步骤S1026, 4艮据上述各等距离的数据点生成标准数据。
通过上述处理，能够将遥感植被指数时间序列数据还原为等距 离的标准数据，这样避免将遥感植被指数时间序列数据在时间轴上 等间隔处理时造成的伪噪声。
(二)步骤S104
发明实施例将遥感植被指数时间序列数据中的反映植被迅速变化特 征的极值点提取出来、生成特征点数据，并在后续处理中保留这些 极值点，从而能够保留植被迅速变化的特征。
图3示出了步骤S104的具体处理流程，如图3所示，步骤S104 的具体处理过程包4舌
步艰《S1042，《殳定查询窗口/=二，其中，/为查询窗口的长度，
r为地表植被生长周期的时长，s为合成周期，该生长周期为植被从
播种到收割的生长过程；优选地，查询窗口 /取与二最接近的奇数为 其值；
例如中国北方的冬小麦一般是在10月底播种，在次年的3 月份开始返青，6月上旬开始收割，这期间遥感植被指数时间序列 数据会显示出两个波峰， 一个波峰对应于播种至返青期，另一个则 对应于返青期至收割期，因此冬d 、麦单个生长周期所覆盖的时长约 为110天，所以查询窗口也应该^隻盖110天左右(即A"-llO ),对于合成周期为16天(即5=16 )的遥感植被指数数据来说，查询窗口 /
的大小应该为7;
步骤S1044，对于分别以标准数据中的每个数据点为起点的多 个查询窗口中，如果查询窗口中的招J直数据点位于查询窗口的中间 位置，则该极值数据点为极值点，将查找到的所有才及值点生成特征 点数据；
另 一方面，由于标准凝:才居中的第一个查询窗口和最后一个查询 窗口中的部分数据中也可能存在极值点，而这两部分数据无法用步 骤S1044依次4仑询到，这样需要对它们进行进一步的识别
步骤S1046,对于标准数据中第一个极值点之前的所有数据点、 和最后一个核 (t点之后的所有lt据点，如果其中的一个数据点与和 该数据点最邻近的极值点的距离为至少半个查询窗口的长度、并且 该数据点的数值和与该数据点最邻近的极值点的数值之间的差值大 于阈值/A，则该数据点为极值点，将查找到的极值点按照合成时间 的先后顺序存入上述特征点数据，优选地，设置/w^A, A为与生 长周期对应的遥感植被指数时间序列数据中的最小波峰的最大值和 最小值的差值；如上所述，植被生长周期的遥感植被指数时间序列 数据可能会包括至少一个生长波峰；
在上述查找到的所有极值点中，可能会由于噪声而存在伪极值 点，这些噪声可能是由于极端的大气环境造成的，例如降雪、长时 间的阴雨天气，这样就需要才艮据极值点的特征和极值点间的关系将 伪极值点去除； 一方面，由于植被生长周期中的最大值和最小值是 相间排列的，如果相邻的两个极值点均为最大值或最小值，则应该 去除其中之一步骤S1048,依次判断特征点数据中任意相邻的两个极值点是
否为最大值点和最小值点，如果均为最大值点，则处理进行到步骤
S1050，如果均为最小值点，则处理进行到步骤S1052，否则，处理 进行到步骤S1054;
步骤S1050,删除两个极值点中数值较小的极值点，处理进行 到步骤S1048;
步骤S1052，删除两个极值点中数值较大的极值点，处理进行 到步骤S1048;
另 一方面，由于极大值和极小值代表了植被生长过程中两个截 然不同的物候状况，这样任意相邻的两个极大值和极小值之间的差 值应该大于一定的阈值，如果它们不满足这一条件，则这两个极值 点中至少有一个数据点可能是伪数据点，优选地，可以设置/^sD,.;
步骤S1054，依次判断特征点数据中任意相邻的两个所述极值 点的差值是否小于阈值/e"，如果判断为是，处理进行到步骤S1056, 否则，处理进行到下述步骤S1060;优选地，还可以i殳置/^s/^, 并且设置2接近于A ,这样可以更为严格地筛除伪极值点；
步骤S1056，删除这两个极值点中的第一个极值点，处理进行 到S1048;
通过上述处理，将标准数据中的核 f直点查找出来并生成特征点 数据，能够记录植被迅速变化特征，为后续处理做准备。
(三)步-骤S106
经过了上述步-骤S102和步骤S104的预处理，可以基于4交为准 确的标准凄t据和特4正点ft据来进行滤波处理(即去噪处理)。
18图4示出了步骤S106的具体处理流程，如图4所示，步骤S106 的具体处理过禾呈包4舌滤波处理和判决处理
1、滤波处5里
步骤S1060,将标准lt据作为第一凄t据，并且设置;t-l, A:为进 行滤波处理的次凄t;
步骤S1062，使用滤波器(l,l,l)(或称为3点变权重
/: + 2 A: + 2 /: + 2
滤波器)对第一数据进行巻积滤波，将滤波后的数据作为第二数据；
当"i时，上述滤波器为均值滤波，随着滤波处理次数的增加， )t的值也越来越大，相应地，滤波器的中间值(即l)的权重也
越来越大，而两边相邻点(即l)的权重则越来越小；本发明实
施例提供的这种变权重滤波策略，每次仅采用3个数据点进行变权 重滤波，滤波幅度小，可以4艮好地^呆持原始遥感才直被指凄t时间序列 数据的曲线形状，并且计算量小、计算速度快，还可以促进算法收 敛，乂人而可以提高运算效率；
但是，随着巻积滤波次数的增加，极值点的植被指数值将会被 逐渐改变，这会导致最后的结果产生物候偏移，并且削弱遥感植#皮 指数时间序列数据的形状特征，因此，极值点的数值应该在滤波数 据中予以保留
步骤S1064，用特征点数据中的所有数据点替换第二数据中相 应位置的数据点，并将替换后的第二数据作为滤波数据；
经过了滤波处理(即步骤S1062的巻积滤波和步骤S1064的特 征点替换)后的滤波数据中可能仍然会存在噪声，这种噪声通常是由于多云天气(云量大于10%)时的大气环境造成的，在滤波数据 中这种噪声表现为幅度较小的数据抖动，如果这种数据抖动不满足
滤波要求，可以对滤波教:据进4亍再次滤波处理 2、判决处理
步骤S1066,如果滤波数据满足预定判决条件则处理进行到步 骤S1068,否则，处理进行到步骤S1070,该预定判决条件包括 z^+1</&3,其中，￡>)+1=^，+1-jv)I，《是第a:次滤波处理后滤波^:据
中数据点y的数值，D"是滤波数据和第一数据(或者为前后相邻两
次滤波处理生成的两个滤波数据)中相应位置的数据点之间的差值， /"3为阈值，/=1A3，".,"， /fe=l，2，3，...，m ,"为滤波凄史据中凄丈据点的个f史， "和附均为自然^t;其中，可以i殳置/"—Z),.;
滤波处理后产生的滤波数据会比第 一数据更为光滑，这就意味 着滤波数据(即新产生的遥感植被指数时间序列数据)含有更少的 噪声；此外，滤波数据和第一数据之间的差异在有噪声的地方会比 在没有噪声的地方更明显，如果该差异在一定的范围内时，即 ￡>)"</"3,则该噪声是可以接受的；即处理可以进行到步骤S1068，
否则对该滤波数据继续进行滤波处理，即处理进行到步骤S1070;
此夕卜，优选地，还可以设置0</"3 </^， /^的值可以才艮据对滤 波曲线的具体要求来确定，/"3的值越小，即对噪声的容忍度越小，
经过滤波处理后的曲线越平滑，但是，需要进行的滤波处理的次数 会越多，这样算法收敛会越慢，即处理速度会相对较慢、处理效率 相对较低；与此相反，/^的值越大，即对噪声的容忍度越大，处理 后的曲线会相对较粗糙，但算法收敛较快，即处理速度相对较快、 处J里凌文率相只于4交高；步骤S1068,将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。
步-骤S1070,如果*>10，即进《亍滤波处理的次^:大于10次，则处理进行到步骤S1072,否则，处理进行到步骤S1074;设置滤波处理次凌t的上限可以防止对纟及端异常的数据点进4于滤波处理的死循环；
步骤S1072，将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。
步骤S1074,将滤波数据作为第一数据，并且，A: = ifc+l，处理转入步骤S1062。
通过上述处理，能够在对标准数据进行滤波的同时，保留植被迅速变化的特征，并且本发明实施例提供的滤波器能够有效地保持
遥感植:故指数时间序列数据的曲线形状，并且运算速度快、效率高，从而能够解决相关技术中遥感植被指数时间序列数据去噪方法存在的无法保留植被迅速变化的特征、去噪后引入新的误差、去噪处理速度慢效率低的问题。
下面对本发明实施例的具体应用进行说明。
应用本发明实施例将江苏省江宁区作为试验区，使用本实施例
提供的方法对江宁区的MODIS (Moderate Resolution ImagingSpectrometer )归一化差值植被指数(即NDVI)进行处理。该NDVI的合成周期为16天，空间分辨率为250米，NDVI时间序列为2006年1月1日至2008年12月31日，3年共69个数据点(每年23个数据点)。图5a、图5b、图5c分别为应用于本发明实施例的实验区中3个l象元点的NDVI时间序列数据，图5a为^M木区i或的NDVI时间序列数据，图5b为一年两季作物冬油菜一夏水稻的NDVI时间序列数据，图5c为一年两季作物冬小麦一夏水稻的NDVI时间序列数据。在图5a、图5b、图5c中，横轴原点时间为2006年1月1日，横轴的单位为天，纵轴为NDVI的数值；图中阴影区域的NDVI为2007年的时间序列凄t据。
图6示出了根据本发明实施例的遥感植被指数时间序列数据去噪方法具体应用的处理流程。
表一中示出了下列处理过程中涉及到的参数(质量等级阈值/ 、、 /"2 、 )。
表一
类型作用步骤取值
质量等级阈值/e^识别时间序列数据中质量较优的数据点S2014
第一阈^f直/A补充时间序列数据两端可能存在的极值点S2050.1
第二阈4直/&2删.除伪扭zft点S2070.15
第三阈值/"3在滤波处理中判别数据点是否为噪声S2110.05
如图6所示，根据本发明实施例的遥感植被指数时间序列数据去噪方法具体应用的处理流程包4舌如下步-豫
步骤S201，根据质量评估文件，从遥感植被指数时间序列数据中选出质量等级^^于质量等级阈值的所有数据点，并将这些数据点作为最优数据，即2'<7 ，/e~=4，其中，e,为质量评估文件中记录的NDVI中数据点y,'的质量等级，y 为预定的质量等级阈值；
步骤S202，根据插值规则《对时间轴上各等距离数据点进行插值，并将这些等距离数据点作为标准数据，其中，
z, =^— 7;)《
Tv.—I'' ' ， A为时间轴上各等距离数据点^的数值，^为最优数据中数据点y'的数值，、为合成日期文件中记录的凡的合
成日期，T、为"'在时间轴上的日期，Z-l，2，3，…，w，"为自然凄t;
/=二
步骤S203,设置查询窗口 s，对于该实验区的NDVI,主要农作物的最短生长波峰长度在r-110天，对于合成周期"16,查询
窗口S 16 ;
步骤S204,依次在以标准数据中的每个数据点为起点的多个查询窗口中查找符合极值点(即特征点)判断条件的数据点，并将找到的所有极值点作为特征数据，该极值点判断条件为查询窗口中的才及^直凄t才居点^立于该查询窗口的中间《立置；
步骤S205，查询标准数据两端可能存在的极值点，对于标准数据中第一个极值点之前的所有数据点、和最后一个才及值点之后的所有凄t据点，如果其中的凄t据点与和该凄t据点最邻近的极值点的距离至少大于3.5 (即至少半个查询窗口的长度)、并且该凄t据点的数值
和与该数据点最邻近的才及值点的数值之间的差值大于，则该数据点为极值点，将查找到的极值点按照合成时间的先后顺序存入上述特征点数据，其中，/A=G-1;
步骤S206,判断特征点数据中任意相邻的两个极值点是否为最大值点和最小值点，如果均为最大值点，则删除其中数值较小的极值点，如果均为最小值点，则删除其中数值较大的极值点；
步骤S207a，如果特征点数据中任意相邻的两个极值点的差值
小于阈值/仏-，处理进^f于到步骤S207b,否则，处理进行到步骤S208，其中，M = 0.15;
23步骤S207b,删除这两个极值点中的第一个极值点，并且，处理转入步-骤S206;
步-骤S208,将标准数据作为第一l史据，并且令*=1;
步骤S209,用滤波器& + 2、 + 2、 + 2对第一数据进行滤波，将滤波后的数据作为第二数据；
步骤S210,用特征点数据中的所有数据点替换第二数据中相应位置的数据点，并将替换后的第二数据作为滤波数据；
步骤S211，如果滤波数据中的所有数据点均满足预定判决条件"「1 < /e",则处理进行到步骤S212，否则，处理进行到步骤S213,
其中/"3-ao5;
步骤S212,将滤波数据作为去噪后的遥感植被指数时间序列数据，至此，处理结束。
步-骤S213，判断滤波处理的次凄t是否大于IO次，即*>10，如果判断为是，则处理进行到步骤S214，否则，处理进行到步骤S215;
步4繁S214，将滤波数据作为去噪后的遥感植;陂指数时间序列数」據，至》匕，处理结束。
步-骤S215，将滤波数据作为第一凄t据，并且，*="1,处理进4亍到步-骤S209。
经过上述处理后，能够在有效地去除试—验区的遥感ti^皮指凄史时间序列数据中存在的噪声，还能够有效地保存植被迅速变化的特征、效率。下面对本发明实施例^是供的方法的去噪效果进行详细i兌明。
(一)保存植;故迅速变化的特征
图7a、图7b、图7c示出了本发明实施例提供的方法与相关技术的方法对试-验区的遥感才直一皮指凄t时间序列数据处理后波形的对比，相关技术中选择目前使用效果较好的非对称高斯函数拟合法(AG)、双逻辑斯蒂函数拟合法(DL)和SG滤波法，图7a、图7b、
技术的方法)在实验区3个像元点上的去噪结果，这三个像元点分别为森林区域、 一年两季作物冬油菜一夏水稻、 一年两季作物冬小麦一夏水稻。在图7a、图7b、图7c中，4黄轴原点时间为2006年1月1日，横轴的单位为天，纵轴为NDVI的数值；图中阴影区域的NDVI为2007年的时间序列数据。
岁文果只于比，力口图7a所示，自2006年至2008年，NDVI时间序列凄史据每年表现为一个波峰，四种方法的去噪结果类似，但本发明实施例提供的方法对2008年1月底至2月初的NDVI才及小值保留得很好，NDVI在此时间会出现这种极小值，主要是2008年1月底至2月初中国南方经历了一次异常的冰冻天气过程，持续时间约一个月，大部分植被均遭到破坏，因此此时的NDVI值明显降低。
图7b显示了刈-试—睑区中一年两季作物冬油菜一夏水稻的NDVI时间序列4改据的处理效果对比，如图7b所示，NDVI时间序列数据每年均会表现为两个波峰，前一个波峰对应于冬油菜，后一个波峰对应于夏水稻，而相关技术的三种去噪方法均把每年的两个波峰合成了 一个没有实际意义的错误波峰，这样4艮容易让后续的用户将其
25判断为森林植一皮类型，而本发明实施例才是供的新方法在除去P桑声的同时却很好的保留了这两个波峰。
图7c显示了对试—睑区中一年两季作物冬小麦一夏水稻的NDVI时间序列数摔的处理效果对比，图7c的效果与图7b的效果类似，这里不再赘述。
从图7a至图7c中可以看出，本发明实施例提供的方法对自然植被的去噪效果优于相关技术的上述三种方法，尤其是对农区植被来说，本发明实施例提供的方法能有效地保持遥感植被指数时间序列数据中的极值点(即特征点)。
图8b示出了才艮据本发明实施例的方法处理NDVI图^f象的效果图。选择实验区中一块大小为400x400像元的区域进行测试，图8a示出了试验区的原始NDVI图像，该图像为2008年6月上旬的一期图像。该实验区(江苏省江宁区)6月上旬还处于梅雨季节，NDVI的质量较差，造成NDVI值偏低并表现明显的马赛克现象，经过本发明的新方法进行滤波去噪后，如图8b所示，偏低的NDVI值基本得到恢复，而且马赛克现象也基本去除。
可见本发明实施例提供的方法能够有效地去除遥感植;陂指数时间序列数据中存在的噪声、并且能够可靠地保持植被迅速变化的特征。
(二)提高处理速度和处理效率
对于实验区的测试影像(400像元x 400像元x 69层)，本发明用交互凄t据语言(Interactive Data Language,简称为IDL) 6.4编写的处理禾呈序在桌面台式木L( 1个中央处理器即CPU,CPU主频为2.21GHz，内存为2 GB随机存储器即RAM)上进行运算，所花费的时间是47s,根据这一时间推算，对于一景标准的MODIS影像(4800^象元x 480(H象元)，如果时间序列为3年(69层)，其处理时间还不 到2小时。只于于SG方法来"i兌，有研究表明(Chen et al" 2004), SG 方法处理一幅8849 <象元x 5601像元x 48层的图<象，用桌面台式枳^ (4个CPU，CPU主频为1.8 GHz，内存为1 GB RAM )花了 22小时， 如果折算成时间序列为3年的标准MODIS影像，它所花的时间约 为14.7小时，这还没有考虑计算机性能上的差异。对于非对称高斯 和双逻辑斯蒂函数拟合法来i兌，它们所花的时间比SG方法更长。
可见，本发明实施例提供的发明能够极大程度地提高遥感植被 指凄t时间序列数据去噪的处理速度和处理效率。
综上所述，根据本发明实施例提供的技术方案，通过对时间轴 上各等距离的数据点插值将遥感植被指数时间序列数据还原为等距 离的标准数据，在标准数据中查找出所有极值点，并根据所有的极 值点对标准数据进4亍滤波，能够有效地去除遥感档::被指数时间序列 数据中的噪声，并且能够保留遥感植被指数时间序列数据中植被迅 速变化的特征；才艮据预定的查询规则能够在标准数据中可靠、全面、 准确地查找出反映植被变化的才及值点，并且根据全部极值点对标准 数据进行滤波处理，能够有效、可靠地保留遥感植被指数时间序列 数据中植被迅速变化的特征；根据预定的滤波器(即3点变权重滤 波器)对数据进行巻积滤波，该滤波器不通过数学函数表达，不会 在滤波处理的同时引入新的误差，并且该滤波器不进行最小二乘法 拟合，从而能够加快处理速度快、提高处理效率。
以上^f又为本发明的优选实施例而已，并不用于限制本发明，对 于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本 发明的精神和原则之内，所作的任何修改、等同替换、改进等，均 应包含在本发明的保护范围之内。
权利要求
1.一种遥感植被指数时间序列数据去噪方法，其特征在于，包括根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各所述等距离的数据点的数值并以此生成标准数据，其中，所述距离为遥感植被指数时间序列数据的合成周期；在所述标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；根据所述特征点数据对所述标准数据按照预定滤波规则进行滤波处理得到滤波数据，如果所述滤波数据满足预定判决条件，则将所述滤波数据作为去噪后的遥感植被指数时间序列数据。
2. 根据权利要求1所述的遥感植被指数时间序列数据去噪方法， 其特征在于，所述对时间轴上各等距离lt据点进行插值包括才艮据所述质量评估文件，在所述遥感植;陂指数时间序列数 据中选出各合成周期内质量最优的数据点，或者选出各合成周 期内质量等级^f氐于预定质量等级阈值的lt据点，并以此生成最 优数据；根据所述最优数据、所述合成日期文件依次对所述时间轴 上各所述等距离数据点按照所述预定插值规则计算其数值，其中，z = H！ *(r -r ) + ^ z,.为所述时间轴上各所述等距离数据点x,的数值，^为所述最优数据中数据点y,的数值，7;为所述合成日期文件中记录的凡.的合成日期，7;为x,.在所述时 间轴上的曰期，！'-1^2，3,…，"，Aj为自然^L
3. 根据权利要求1所述的遥感植被指数时间序列数据去噪方法，其特4i在于，所述预定查询井见则包^舌:没置查询窗口/=匸，其中，/为查询窗口的长度，r为地表植被生长周期的时长，s为所述合成周期，所述生长周期为 植被从播种到收割的生长过程；对于分别以所述标准数据中的每个数据点为起点的多个 所述查询窗口中，如果所述查询窗口中的极值数据点位于所述 查询窗口的中间位置，则该才及值数据点为所述招J直点。
4. 根据权利要求3所述的遥感植被指数时间序列数据去噪方法， 其特4i在于，所述预定查询^见则还包4舌对于所述标准凄t据中第 一个所述招 f直点之前的所有凄t据 点、和最后一个所述才及值点之后的所有凄t据点，如果其中的一 个数据点和与该数据点最邻近的所述极值点的距离为至少半 个所述查询窗口的长度、并且该#:据点的数值和与该#:据点最 邻近的所述极值点的数值之间的差值大于第一预定阈值凤， 则该凝:据点为所述极值点，并将查找到的所述极值点按照其合 成曰期的先后顺序存入所述特4正点数据中；其中，所述/^小于或等于A，所述A.为与所述生长周期 对应的遥感植被指数时间序列数据中最小波峰的最大值和最 小值的差值。
5. 4艮据4又利要求4所述的遥感才直4皮指数时间序列数据去噪方法， 其特4i在于，所述预定查询^见则还包括第一去伪处理如果所述特征点数据中任意相邻的两个所 述极值点均为最大值点，则删除其中数值较小的极值点，如果 该相邻的两个所述极值点均为最小值点，则删除其中数值较大 的招 f直点；第二去伪处理如果所述特4正点凄G居中任意相邻的两个所 述极值点的差值小于第二预定阈值/"2,则删除该相邻的两个 所述招J直点中的第 一个所述才及值点，并重新进4亍所述第 一去伪 处理；其中，所述/^2小于或等于所述化。
6. 根据权利要求3至5中任一项所述的遥感植被指数时间序列数 据去噪方法，其特征在于，所述查询窗口/取与二最接近的奇数为其值。
7. 根据权利要求5所述的遥感植被指数时间序列数据去噪方法， 其特征在于，所述/^大于或等于所述/a。
8. 根据权利要求1至3中任一项所述的遥感植被指数时间序列数 据去噪方法，其特征在于，所述滤波处理包括将所述标准数据作为第一数据，使用 滤波器对所述第一数据进行巻积滤波并生成第二数据，其中，所述滤波器为(l，l，l),A:"，2，3，…，m ， A:为当前进行所述 fc + 2 A: + 2 A: + 2滤波处理的次数，m为自然数；用所述特征点数据中的所有数据点替换所述第二数据中相应位置的#:据点，并将替换后的所 述第二凄t据作为滤波数据；所述预定判决条件包括村+1</"3,如果所述滤波数据不 满足所述预定判决条件，则将所述滤波数据作为第一数据，并 对该第一数据进4于所述滤波处理，直至所述滤波lt据满足所述 预定判决条件，其中，z^+1=|jvf+1-《I，《是第;t次所述滤波处理后所述滤波凄t据中数据点的数值，化"1是相邻两次所述滤 波处理生成的两个所述滤波数据中相应位置的凄t据点之间的差^直，/^3为第三予贞定阈<直，/ =1,2,3，...,"， A: =l，2，3，..',m ,"为所 述滤波数据中数据点的个数，az和"!均为自然数；其中，所述/A小于或等于D,.,所述Z),为与所述生长周期 对应的遥感植被指数时间序列数据中的最小波峰的最大值和 最小值的差值。
9. 根据权利要求8所述的遥感植被指数时间序列数据去噪方法， 其特征在于，所述预定判决条件还包括对所述第一数据进行所述滤波处理的次数大于10次。
10. 根据权利要求8所述的遥感植被指数时间序列数据去噪方法， 其特4i在于，所述/^小于所述/K。
全文摘要
本发明提供了一种遥感植被指数时间序列数据去噪方法，该方法包括根据遥感植被指数时间序列数据、合成日期文件、和质量评估文件对时间轴上各等距离的数据点按照预定插值规则进行插值，得到各等距离的数据点的数值并以此生成标准数据，其中，距离为遥感植被指数时间序列数据的合成周期；在标准数据中按照预定查询规则查找出所有极值点并以此生成特征点数据；根据特征点数据对标准数据按照预定滤波规则进行滤波处理得到滤波数据，如果滤波数据满足预定判决条件，则将滤波数据作为去噪后的遥感植被指数时间序列数据。根据本发明，能够有效地去除遥感植被指数时间序列数据中的噪声，能够保留植被迅速变化的特征，并且处理速度快、效率高。
文档编号G01S7/48GK101650422SQ20091017713
公开日2010年2月17日 申请日期2009年9月27日 优先权日2009年9月27日
发明者朱文泉 申请人:北京师范大学