本发明涉及地理信息与遥感技术领域，尤其涉及一种时序的ndvi数据序列的优化方法。
背景技术：
ndvi(normalizeddifferencevegetationindex)归一化差分植被指数，也称为生物量指标变化。ndvi是植被生长状态及植被覆盖度的最佳指示因子。ndvi能够准确的反映出地表植被的活力状态和植被季相变化特征。基于这一特性，ndvi现已被广泛的应用于各种变化的监测中，如城市扩张监测、森林扰动监测、植被长势监测、气候变化监测、开采扰动监测等。原始的时序ndvi数据是通过对卫星遥感器获取的遥感数据进行辐射定标、大气校正等计算得到的。遥感器获取遥感数据的过程容易受到云、气溶胶、太阳高度角和地物双向性反射等一系列不良因素的影响，导致遥感器获取的遥感数据中包含了一定量的噪声数据。这也导致在此基础上计算得出的时序ndvi数据中存在一定量的噪声数据。当上述时序ndvi数据用于变化检测时，噪声数据会在一定程度上影响检测结果的真实性，导致检测结果不能体现被检测对象的真实变化，使人们根据检测结果做出错误的判断。
常见的时序ndvi数据序列的优化方法，如采用最佳坡度系数截取法bise优化时序ndvi数据序列，虽能很大程度地去除时序ndvi数据序列中的噪声数据，但优化后的时序ndvi数据序列仍然具有较高的离散性，尤其是在连续噪声处bise表现不佳；采用savitzky-golay滤波法优化时序ndvi数据序列，在去除时序ndvi数据序列中的噪声数据的同时，还可以降低优化后的时序ndvi数据序列中数据的离散性，但优化后的时序ndvi数据序列与真实的数据序列的时相与幅度存在一定误差，目前该误差无法消除。
综上所述，如何去除时序ndvi数据序列中的这些噪声数据，并尽可能的保证时序ndvi数据的真实性，是目前我们面临的首要问题。
技术实现要素：
为解决上述问题，本发明提供了一种时序ndvi数据序列的优化方法，包括如下步骤：
获取初始的第一时序ndvi数据序列；在限制条件下，采用最佳坡度系数截取法bise对第一时序ndvi数据序列进行优化，得到第二时序ndvi数据序列；采用小波变换法对第二时序ndvi数据序列进行优化，得到目标时序ndvi数据序列。
优选地，获取初始的第一时序ndvi数据序列步骤，包括：获取原始卫星遥感数据；对原始卫星遥感数据进行辐射定标；对定标后的数据进行大气校正；对大气校正后的数据进行配准和几何校正,并对计算得出的ndvi数据按时间先后顺序重新组合，得到第一时序ndvi数据序列。
优选地，特定的限制条件的表达式为：
其中xi-1、xi、xi+1为在第一时序ndvi数据序列中选取的3个连续的数据,i为大于等于2的整数。
优选地，在限制条件下，采用最佳坡度系数截取法bise对第一时序ndvi数据序列进行优化,得到第二时序ndvi数据序列步骤，包括:
数据采集：在第一时序ndvi数据序列中依次选取3个的连续的数据：xi-1、xi、xi+1；
数据预处理：根据公式
数据优化：设定一个阈值ti，当斜率slop(i-1,i)的绝对值|slop(i-1,i)|和斜率slop(i+1,i)的绝对值|slop(i+1,i)|同时大于阈值ti，且满足
优化数据的整理：按时序先后逐一选取第一时序ndvi数据序列中的数据，并对所选取的数据进行优化，优化后的数据按时序由先到后重新组合，得到第二时序ndvi数据序列。
优选地，采用小波变换法对第二时序ndvi数据序列进行优化，得到目标时序ndvi数据序列步骤，包括：采用预先设定的小波对第二时序ndvi数据序列进行n层分解计算，n为预先设定的分解层数，n≥2；采用预先设定的阈值模式y，对第1层到第n层的每一层高频系数进行阈值量化处理；根据小波分解的第n层的低频系数和经过量化处理后的第1层到第n层的高频系数，对第二时序ndvi数据序列进行小波重构，得到目标信号波；从目标信号波中提取ndvi数据，将提取到的ndvi数据按时间先后顺序重新组合得到目标时序ndvi数据序列。
优选地，预先设定的小波为db7小波；分解层数n具体为：2；特定阈值模式y设置为：极大极小阈值准则：minimaxi，且设置为软阈值：soft。
采用本发明的方法优化初始的时序ndvi数据序列，能有效的去除初始的时序ndvi数据序列中的噪声数据，同时保持了时序ndvi数据的真实性，使得ndvi能够更为准确的体现被监测目标的真实变化信息。
附图说明
图1为本发明实施例提供的一种时序ndvi数据序列的优化方法流程图；
图2为本发明获取初始的第一时序ndvi数据序列步骤的流程图；
图3为本发明在限制条件下，采用最佳坡度系数截取法bise对所述第一时序ndvi数据序列进行优化,得到第二时序ndvi数据序列步骤的流程图；
图4为采用小波变换法对所述第二时序ndvi数据序列进行优化，得到目标时序ndvi数据序列步骤的流程图；
图5为检测区域一：中国内蒙宝日希勒镇
图6为检测区域二：美国阿拉巴拉契亚马丁镇
图7为第一样点的数据优化对比图；
图8为第二样点的数据优化对比图；
图9为第三样点的数据优化对比图；
图10为第四样点的数据优化对比图；
图11为采用最佳坡度系数截取法bise优化后的时序ndvi数据序列得到的变化检测结果；
图12为采用savitzky-golay滤波法优化后的时序ndvi数据序列得到的变化检测结果；
图13为采用本发明所记载的新方法优化后的时序ndvi数据序列得到的变化检测结果；
图14为采用原始时序ndvi数据序列得到的变化检测结果。
具体实施方式
下面结合附图和实施例，对本发明的技术方案做进一步的详细描述。
图1为本发明实施例提供的一种时序ndvi数据序列的优化方法流程图。如图1所示，时序ndvi数据序列的优化方法包括以下步骤：
s101，获取初始的第一时序ndvi数据序列。
具体地，该步骤具体的实施过程，如图2所示，包括：首先，获取原始卫星遥感数据；然后，对该原始卫星遥感数据进行辐射定标；紧接着，对定标后的数据进行大气校正；然后，对大气校正后的数据进行配准和几何校正；最后，根据配准和几何校正后的遥感数据计算遥感数据的ndvi数据并将计算得出的ndvi数据按时间先后顺序重新组合，得到第一时序ndvi数据序列。
s102，在限制条件下，采用最佳坡度系数截取法bise对第一时序ndvi数据序列进行优化,得到第二时序ndvi数据序列。
具体地，限制条件的表达式为：
该步骤具体的实施过程，如图3所示，包括：
首先，进行数据采集，即：在第一时序ndvi数据序列中依次选取3个的连续的数据：xi-1、xi、xi+1；
然后，对采集到的数据进行预处理：根据公式
紧接着，对采集到的数据进行优化：设定一个阈值ti，当斜率slop(i-1,i)的绝对值|slop(i-1,i)|和斜率slop(i+1,i)的绝对值|slop(i+1,i)|同时大于阈值ti,且满足限制条件公式：
最后，优化后的数据整理：按时序先后逐一选取第一时序ndvi数据序列中的数据，并对所选取的数据进行优化，优化后的数据按原时序由先到后重新组合，得到第二时序ndvi数据序列。
s103，采用小波变换法对第二时序ndvi数据序列进行优化，得到目标时序ndvi数据序列。
具体地，该步骤具体过程，如图4所示，包括：
首先，采用预先设定的小波对第二时序ndvi数据序列进行n层分解计算，n为预先设定的分解层数，n≥2；其次，采用预先设定的阈值模式y，对第1层到第n层的每一层高频系数进行阈值量化处理；再次，根据小波分解的第n层的低频系数和经过量化处理后的第1层到第n层的高频系数，对第二时序ndvi数据序列进行小波重构，得到目标信号波；最后，从目标信号波中提取ndvi数据,将提取到的ndvi数据按时间先后顺序重新组合得到目标时序ndvi数据序列。
优选的，预先设定的小波为db7小波；分解层数n具体为：2；特定阈值模式y设置为：极大极小阈值准则：minimaxi，且设置为软阈值：soft。可以利用数据仿真软件matlab实现本方法。
为验证新方法的可靠性，实施例一选取了代表不同植被类型的两个区域作为目标检测区域。所述两个区域为：植被覆盖以草地为主的中国内蒙宝日希勒镇(如图5所示)和植被覆盖以森林为主的美国阿拉巴拉契亚马丁镇(如图6所示)。
分别在两个区域选取200个样点，将每个样点所对应的时序ndvi数据序列作为原始数据序列。
分别采用最佳坡度系数截取法bise、savitzky-golay滤波法以及本发明记载的新方法对原始数据序列进行优化处理。
随机抽取第一、第二、第三、第四，4个样点，对经过上述三种方法优化处理后的数据进行目视比较，如图7-10所示。
其中，图7中的原始数据为第一样点对应的时序ndvi数据序列。图7中的图1-1为采用最佳坡度系数截取法bise优化处理后的数据与原始数据的对比图；图7中的图1-2为采用savitzky-golay滤波法优化处理后的数据与原始数据的对比图；图7中的图1-3为采用本发明记载的新方法优化处理后的数据与原始数据的对比图。
相应地，图8中的原始数据为第二样点对应的时序ndvi数据序列。图8中的图2-1为采用最佳坡度系数截取法bise优化处理后的数据与原始数据的对比图；图8中的图2-2为采用savitzky-golay滤波法优化处理后的数据与原始数据的对比图；图8中的图2-3为采用本发明记载的新方法优化处理后的数据与原始数据的对比图。
相应地，图9中的原始数据为第三样点对应的时序ndvi数据序列。图9中的图3-1为采用最佳坡度系数截取法bise优化处理后的数据与原始数据的对比图；图9中的图3-2为采用savitzky-golay滤波法优化处理后的数据与原始数据的对比图；图9中的图3-3为采用本发明记载的新方法优化处理后的数据与原始数据的对比图。
相应地，图10中的原始数据为第四样点对应的时序ndvi数据序列。图10中的图4-1为采用最佳坡度系数截取法bise优化处理后的数据与原始数据的对比图；图10中的图4-2为采用savitzky-golay滤波法优化处理后的数据与原始数据的对比图；图10中的图4-3为采用本发明记载的新方法优化处理后的数据与原始数据的对比图。
分别对比图7中的图1-1、1-2、1-3，图8中的图2-1、2-2、2-3，图9中的图3-1、3-2、3-3，图10中的图4-1、4-2、4-3；宏观上可以看出采用最佳坡度系数截取法bise优化后的数据仍然保持着较高的离散型，最佳坡度系数截取法bise在连续噪声数据处的优化效果不佳；savitzky-golay滤波法能很好地去除原始数据序列中的噪声数据，但是在数据陡降时，savitzky-golay滤波法会改变陡降处数据的时相，导致优化后的数据不准确。采用本发明实施例所记载的新方法，可以有效的消除掉时序ndvi数据序列中的噪声数据，同时保持了时序ndvi数据的真实性。
实施例二：通过随机选取样本的方式，确定研究区域的ndvi植被阈值；
记录从有植被到无植被变化的年份，以该年作为样例区域的扰动年份；
随机选取两百个样本点，对不同方法得到时序ndvi数据序列及原始时序ndvi数据序列的变化检测结果进行验证。
对不同方法得到的优化数据与原始数据的变化检测结果进行精度评价和结果分析。
图11为采用最佳坡度系数截取法bise优化后的时序ndvi数据序列得到的变化检测结果；
图12为采用savitzky-golay滤波法优化后的时序ndvi数据序列得到的变化检测结果；
图13为采用本发明所记载的新方法优化后的时序ndvi数据序列得到的变化检测结果；
图14为采用原始时序ndvi数据序列得到的变化检测结果；
通过对上述变化检测结果进行精度评价可以得出以下精度检测结果：
表1精度检测结果
从表1中可以看出三种滤波方法都可以提高变化检测的真实度，其本发明实施例记载的新方法在变化检测中表现最佳，将整体精度提高57％。
采用本发明实施例所记载的新方法，可以有效的消除掉时序ndvi数据序列中的噪声数据，同时保持了时序ndvi数据的真实性。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。