CN103902802A

CN103902802A - 一种顾及空间信息的植被指数时间序列数据重建方法

Info

Publication number: CN103902802A
Application number: CN201210589550.2A
Authority: CN
Inventors: 梁守真; 陈劲松; 李洪忠; 李晓明; 崔晓伟; 王贺
Original assignee: Shenzhen Institute of Advanced Technology of CAS
Current assignee: Shenzhen Institute of Advanced Technology of CAS
Priority date: 2012-12-29
Filing date: 2012-12-29
Publication date: 2014-07-02
Anticipated expiration: 2032-12-29
Also published as: CN103902802B

Abstract

本发明提供一种顾及空间信息的植被指数时间序列数据重建方法，包括S1、NDVI数据的数值转换；S2、建立冰/雪和云掩膜文件；S3、替换冰/雪覆盖像元NDVI值；S4、生成年际NDVI时间序列数据；S5、土地覆盖数据预处理；S6、对NDVI时间序列数据进行空间域处理；S7、对NDVI时间序列数据进行时间域处理。本发明的植被指数时间序列数据重建方法，充分利用冰/雪和云掩膜文件，从时间维和空间维两个角度NDVI对时间序列进行重建，通过上述方法不仅可以有效去除云、冰/雪对NDVI数据的影响，并且重建后的时间序列在最大程度上保留了植被的变化信息，弥补了当前技术方法的不足。

Description

一种顾及空间信息的植被指数时间序列数据重建方法

【技术领域】

本发明涉及植被遥感领域，尤其涉及一种顾及空间信息的植被指数时间序列数据重建方法。

【背景技术】

在植被遥感领域，通常采用植被指数来对地表植被进行研究。由于卫星传感器所获取的归一化植被指数NDVI(Normalized Difference Vegetation Index)时间序列数据能够精确地反映地表植被的活力状态和植被季相变化特征，NDVI成为当前最常用的一种植被指数，并已在全球、洲际以及地区范围的植被覆盖制图、陆地生态系统动态监测和模拟方面发挥了重要作用。

目前常用的NDVI时间序列数据大多是利用多天的数据应用最大值合成法MVC(Maximum Value Composite)或者限制观测角的最大值合成法CV-MVC(Constrained-view Angle Maximum Value Composite)生成。最大值合成法可以消除一定的云、气溶胶、水汽的影响，但是由于算法本身的缺陷和环境影响的复杂性，生成的NDVI数据仍然受残余云、大气以及地表二向性反射等因素的影响。时间序列数据整体波动较大，相邻值高低变化没有规律，曲线季节变化趋势不明显，使得NDVI时间序列数据无法进行各种趋势分析和信息提取。因此，为更真实地反映植被的动态，还需对遥感植被指数时间序列数据进行重建，最大程度上去除或降低数据的噪声水平，提高数据质量。

目前已经发展的多种NDVI时间序列数据的重建方法，包括最佳指数斜率提取法BISE(Best Index Slope Extraction)、时间窗口线性内插法TWO(TemporalWindows Operation)、均值迭代滤波法MVI(Mean-value Iteration Filter)、Savizky-Golay法、时间序列谐波分析法HANTS(Harmonic ANalysis of TimeSeries)、非对称高斯函数拟合法(Asymmetric Gaussian Funtion-fitting)、双Logistic(Double Logistic Funtion-fitting)函数拟合法等。

上述时序数据重建方法都基于以下两个假设：①NDVI的时序变化对应着植被的生长和衰老过程；②植被生长是一个渐进的变化过程，时序数据中的突变点不符合植被的生长规律，应作为噪声点。最佳指数斜率提取法是采用滑动窗口识别并替代序列中的噪声，其是在一个设定的滑动窗口内，从第一个点向前搜索，如果下一点的值高于（第一个点）起始点的值，则接受该点；当遇到下降点时，如果在该滑动窗口内不存在高于该低值和先前高值之差的20%的点，该点被接受，否则认为是噪声，对噪声点采用线性拟和的方法补充。时间窗口线性内插法是在设定的窗口内从始点寻找离它最近的比它大的值作为下一个窗口的起始点，如果窗口内没有找到比开始点更大的值，那么选择离它最近的除它之外的最大值作为下一个起点，当前点和下一点之间的时间点被认为是噪声影响点，对其进行线性插值。均值迭代滤波法，首先计算时间序列某一点与其前后相邻两个点的平均值的差值，如果该差值超过设定的阈值，认为是噪声点，用前后两点的平均值替换该点的值，反复迭代进行，直到NDVI序列中不再存在满足条件的点。Savizky-Golay法是用一定长度的窗口和待处理数据作卷积，对待处理的数据作加权多项式拟和，拟合的目标是求得最小均方根误差，而一些远离大多数点的边沿点不参与拟合。时间序列谐波分析法又称为傅立叶变换法，该方法首先将时间序列信号进行傅立叶变换，之后选择谐波组分进行反傅立叶变换生成新的时间序列，检验出偏离时间序列拟合曲线较远的数据，并将它们去掉，然后用傅立叶特征分量计算出的时间插值去填充。双Logistic函数拟合法和非对称高斯拟合法相似，都是利用固定的函数形式用优化算法去拟合NDVI时间序列曲线，然后将拟合曲线的点作为重建的时间序列。

总体上，现有技术方法大致可分为三类，第一类是基于移动窗口的处理方法，如最佳指数斜率提取法、时间窗口线性内插法、Savizky-Golay法，均以移动窗口作为处理单元，利用一定的规则检测噪声点，对窗口内数据进行线性或非线性拟合，并用拟合值代替噪声影响的点；第二类是谐波分析法或傅立叶变换法，该方法是利用傅立叶变换将时间序列从时间域转换到频率域，通过选择低频信息、去除高频组分来重建时间序列；第三类是模型拟合法，该方法假设植被的生长过程符合某种数学模型，首先利用时间序列数据解算模型的各个参数，然后模型的拟合值来代替原始时间序列的值，双Logistic函数拟合法和非对称高斯拟合法都是属于该种方法。

尽管上述方法能够实现NDVI时间序列的重建，并在许多方面得到了运用，但是仍然存在一些问题：

首先，仅从时间维处理，缺乏对空间信息的考虑。上述NDVI时间序列重建方法均是从时间维的角度出发，单纯从时间序列数据本身的相关性考虑，使用前后多期的数据拟合曲线或平均值来替代受噪声影响像元的NDVI值，很少考虑地物在空间上的相关性，从空间角度来处理数据。上述技术方法在处理NDVI时间序列数据时，假设植被生长是一个缓慢的渐变过程，环境胁迫或人为干扰导致的迅速变化被认为是噪声而被处理掉，但这些信息往往是很重要的，应该保留而不应处理掉。在NDVI图像上，相邻像元在空间上往往具有明显的依赖性或者相关性。如能结合像元在空间上的相关性信息，在一定程度上可以保留这些突变点，避免重要信息的丢失。因此，在重建时间序列时，需要考虑空间信息。

其次，引入新的误差，不利于定量分析。时间序列重建的目的是消除噪声的影响，但是目前现有技术方法在消除噪声影响时也会导致部分高质量像元的值发生改变，从而引入新的误差。上述技术方法大都是在分析美国国家海洋和大气局(National Oceanic and Atmospheric Administration,NOAA)提供的NDVI时间序列数据的基础上发展起来的，由于该时序数据没有相应的质量控制信息，所以构建的时间序列重建方法通常需要先检测噪声，然后再对噪声点进行处理。因此，相关技术在时间序列数据重建时引入的误差主要来自两个方面，一是噪声检测错误而导致的误差，这是由于混合像元的存在和植被生长过程的复杂性，根据假设条件检测到的像元未必是受到噪声的影响；二是拟合值替代导致的误差。尽管基于函数拟合和谐波分析的方法不需要检测噪声像元，但是该方法采用拟合值代替时间序列中所有点的NDVI值，导致非噪声影响像元的值也发生改变。最佳指数斜率提取法和时间窗口线性内插法可以保证非噪声点像元的值不变，但是简单的线性插值得到的结果往往精度较低。对某些应用来说，上述方法对时间序列的处理结果能满足要求，但是在定量分析NDVI的时空分布和变化时可能会产生较大误差，尤其是在受云雨影响严重的中国南方地区。

【发明内容】

本发明要解决的技术问题在于克服现有时间序列重建方法的不足，构建一种能将时间和空间信息有效结合并保持高质量像元点值不变的时间序列数据重建方法。

为此，本发明提供一种顾及空间信息的植被指数时间序列数据重建方法，包括以下步骤，

S1、NDVI数据的数值转换；

S2、建立冰/雪和云掩膜文件；

S3、替换冰/雪覆盖像元NDVI值；

S4、生成年际NDVI时间序列数据；

S5、土地覆盖数据预处理；

S6、对NDVI时间序列数据进行空间域处理；

S7、对NDVI时间序列数据进行时间域处理。

进一步地，所述步骤S1具体为将每个合成期文件的NDVI数据值转换到[-1,1]。

进一步地，所述步骤S2具体为从每个合成期的质量评估文件中提取冰/雪和云覆盖像元，建立每个合成期的冰/雪和云掩膜文件。

进一步地，所述步骤S3具体为利用每个合成期的冰/雪和云掩膜文件，对冰/雪和云掩膜文件相应的NDVI图像进行掩膜，并分别计算每个合成期的无云和无冰/雪覆盖像元的多年均值，用均值代替每个合成期冰/雪覆盖像元NDVI值。

进一步地，所述步骤S4具体为利用初始化NDVI数据，生成年际的NDVI时间序列数据。

进一步地，所述步骤S5预处理具体为投影转换，使土地覆盖数据NDVI数据具有相同的投影。

进一步地，土地覆盖数据为矢量数据，并且步骤S5的所述预处理还包括矢量数据向栅格数据的转换。

进一步地，所述步骤S6具体为基于土地覆盖数据及像元之间的相关性,对云覆盖的像元进行空域处理。

进一步地，所述步骤S7具体为采用样条插值法对时间序列中云覆盖的像元进行插值计算。

本发明空间信息的植被指数时间序列数据重建方法，充分利用冰/雪和云掩膜文件，从时间维和空间维两个角度NDVI对时间序列进行重建，通过上述方法不仅可以有效去除云、冰/雪对NDVI数据的影响，并且重建后的时间序列在最大程度上保留了植被的变化信息，弥补了当前技术方法的不足。

【附图说明】

图1示出本发明的指数时间序列重建方法的流程图。

图2示出本发明的植被指数数据原始数据和结果数据图。

图3示出本发明的植被指数数据原始时间序列和结果时间序列图。

图4示出本发明的植被指数数据重建方法的比较图。

【具体实施方式】

下面结合附图和具体实施例对本发明作进一步详细说明。

如图1所示，本发明的顾及空间信息的植被指数时间序列数据重建方法，步骤如下：

S1、NDVI数据的数值转换。具体为将每个合成期文件的NDVI数据值转换到[-1,1]。由于现有的NDVI数据大都通过线性转换将其数值转换为整数型或字节型进行存储，以节省存储空间，所以其数据范围并非真实的NDVI值的范围[-1,1]，在应用时有必要对其数值进行反转换。

S2、建立冰/雪和云掩膜文件。具体为从每个合成期的质量评估文件中提取冰/雪和云覆盖像元，建立每个合成期的冰/雪和云掩膜文件。

由于每个合成的NDVI数据都有一个对应的质量评估文件，该文件记录了NDVI数据中每个像元的质量，在NDVI数据中需要处理的像元主要是冰/雪和云覆盖的像元，因此从每个合成期的质量评估文件中提取冰/雪和云覆盖像元，建立每个合成期的冰/雪和云掩膜文件。

S3、替换冰/雪覆盖像元NDVI值。

在中国，冰/雪覆盖的像元主要出现在冬季的北方和海拔较高的地区，并且冬季大部分植被都停止生长，因此NDVI在年际之间的差异较小。因此可以利用每个合成期多年高质量像元的均值来代替合成期冰/雪覆盖像元的NDVI。

具体为利用每个合成期的冰/雪和云掩膜文件，对冰/雪和云掩膜文件相应的NDVI图像进行掩膜，并分别计算每个合成期的无云和无冰/雪覆盖像元的多年均值，用均值代替每个合成期冰/雪覆盖像元NDVI值。

S4、生成年际NDVI时间序列数据。具体为利用初始化NDVI数据，生成年际的NDVI时间序列数据。

S5、土地覆盖数据预处理。预处理具体为投影转换，使土地覆盖数据NDVI数据具有相同的投影。若土地覆盖数据为矢量数据，上述预处理还包括矢量数据向栅格数据的转换。

S6、对NDVI时间序列数据进行空间域处理。具体为基于土地覆盖数据及像元之间的相关性,对云覆盖的像元进行空域处理。

在每个合成期的NDV图像上建立一个3×3的移动窗口，如果NDVI窗口中心的像元为云覆盖像元，则寻找周围临近的8个像元中是否存在无云覆盖的像元。

此处有三种情况：①、不存在满足条件的像元，则不做任何处理；②、如果存在一个满足条件的像元，则根据土地覆盖数据判断满足条件像元的土地覆盖类型是否与中心像元一致：如果一致，则用该像元的NDVI值替代云覆盖的像元NDVI值；如果不一致，则不做任何处理；③、如果存在多个满足条件的像元，则根据土地覆盖数据判断哪些像元与中心像元的土地覆盖类型一致：如果只有一个像元与中心像元的一致，则用其NDVI替代中心像元的NDVI；如果有多个像元的土地覆盖类型与中心像元一致，首先根据它们各自的NDVI时间序列，选择相应时期的无云覆盖像元分别统计它们与中心像元的相关性，选择相关程度最高的临近像元的NDVI值来替代中心云覆盖像元的NDVI值。

计算相关性，采用相关系数作为判定标准，其计算公式为：

R_{xy} = \frac{Σ_{i = 1}^{n} (x_{i} - \overset{&OverBar;}{x}) (y_{i} - \overset{&OverBar;}{y})}{\sqrt{Σ_{i = 1}^{n} {(x_{i} - \overset{&OverBar;}{x})}^{2}} \times \sqrt{Σ_{i = 1}^{n} {(y_{i} - \overset{&OverBar;}{y})}^{2}}}

式中，R_xy表示待处理像元和临近像元之间的相关系数，

和

分别表示待处理像元和其临近像元时间序列中不受云影响的点的平均值。

S7、对NDVI时间序列数据进行时间域处理。具体为采用样条插值法对时间序列中云覆盖的像元进行插值计算。

采用样条插值法对时间序列中云覆盖的像元进行插值计算。采用样条插值法的好处在于使用低阶的多项式样条而能产生和高阶的多项式插值类似的效果，插值误差较小，并且避免了使用高阶多项式所出现的龙格现象。为避免在时间序列两侧进行连续插值的时出现异常值，在待处理的时间序列两侧补充一定长度的时间序列，可以选择前一年和后一年的部分时间或全年的时间序列。

本发明的顾及空间信息的植被指数时间序列数据重建方法，充分利用冰/雪和云掩膜文件，从时间维和空间维两个角度NDVI对时间序列进行重建，通过上述方法不仅可以有效去除云、冰/雪对NDVI数据的影响，并且重建后的时间序列在最大程度上保留了植被的变化信息，弥补了当前技术方法的不足。

实施例1

本实施例进一步具体说明本发明提供的顾及空间信息的植被指数时间序列数据重建方法。

数据采自辽宁省MODIS（Moderate Resolution Imaging Spectroradiometer）16天合成的NDVI数据和土地覆盖数据，共获取2003-2007年5年的NDVI数据，其中2005年的NDVI时间序列数据作为待重建的案例数据。由于MODIS的NDVI数据是同质量评估数据以及其他数据集成在同一个数据集中，并且其投影方式为正弦投影(SIN)，因此在进行数据重建前先利用MODIS重投影工具MRT(MODIS Reprojection Tool)将NDVI数据及其质量评估数据从数据集中提取并转为地理坐标。

步骤S1、MODIS NDVI数据的数值转换。

MODIS NDVI数据的原始有效数据范围为[-2000,10000]，非常见的[-1,1]，由于该数据的缩放因子为0.0001，所以其转换公式为

NDVI=0.0001*DN

式中，DN为原始数据的NDVI值。通过该公式即可将数据的范围转换到[-1,1]。

步骤S2、建立冰/雪和云覆盖像元掩膜。

MODIS的质量评估数据记录了NDVI数据中每个像元的质量。在MODIS质量评估数据-像元可靠度图像中，冰/雪和云像元的灰度值分别为2和3，利用两个标识分别建立每个合成期的冰/雪和云掩膜文件。

步骤S3、替换冰/雪覆盖像元NDVI值。

利用每个合成期的冰/雪和云掩膜文件对与其相对应的NDVI数据进行掩膜，掩膜值可以定义为[-1,1]范围之外的任何值，在本实例中冰/雪覆盖区域的值设为-10，云覆盖区域的掩膜值设为-20；之后计算2003-2007年每个合成期非冰/雪和云覆盖像元的平均值，用计算的均值替代相应合成期冰/雪覆盖像元的NDVI值。

步骤S4、生成年际NDVI时间序列数据。

将每年的NDVI数据进行层堆栈(Layer Stacking），生成每年的NDVI时间序列数据。MODIS NDVI数据采用16天合成，每年有23个合成期的数据，这样，在时间序列数据中每个像元有23个值。

步骤S5、土地覆盖数据预处理。

由于获取的土地覆盖数据为矢量数据，首先需将矢量数据转化为栅格数据，空间分辨率与NDVI数据一致，然后对其投影进行转换，以便相同地物在NDVI图像和土地覆盖数据上具有相同的坐标。

步骤S6、对NDVI数据的进行空间域处理。

计算相关性，采用相关系数作为判定标准，其计算公式为：

R_{xy} = \frac{Σ_{i = 1}^{n} (x_{i} - \overset{&OverBar;}{x}) (y_{i} - \overset{&OverBar;}{y})}{\sqrt{Σ_{i = 1}^{n} {(x_{i} - \overset{&OverBar;}{x})}^{2}} \times \sqrt{Σ_{i = 1}^{n} {(y_{i} - \overset{&OverBar;}{y})}^{2}}}

式中，R_xy表示待处理像元和临近像元之间的相关系数，

和

步骤S7、对NDVI时间序列数据进行时间域处理。

在待处理的2005年的时间序列两侧补充2004年和2006年的NDVI时间序列数据，采用样条插值法对2005年时间序列中云覆盖的像元进行插值计算。

如图2所示，经过本方法重建后的NDVI数据已经消除了云和冰雪影响，数据质量大幅提高；如图3所示，为研究区某个像元的重建前后的时间序列，从图中可以看出，由于受冰雪或云影响而导致的NDVI偏低的像元已经得到纠正，重建的时间序列能更好地反应植被随时间的变化特征；并且相比于目前常用的其他方法，例如Savizky-Golay法、时间序列谐波分析法、非对称高斯函数拟合法和双Logistic函数拟合法，如图4所示，本发明的重建方法能更好地保留NDVI时间序列的细节特征。

以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。

Claims

1.一种顾及空间信息的植被指数时间序列数据重建方法，其特征在于，包括以下步骤，

S1、NDVI数据的数值转换；

S2、建立冰/雪和云掩膜文件；

S3、替换冰/雪覆盖像元NDVI值；

S4、生成年际NDVI时间序列数据；

S5、土地覆盖数据预处理；

S6、对NDVI时间序列数据进行空间域处理；

S7、对NDVI时间序列数据进行时间域处理。

2.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，步骤S1具体为将每个合成期文件的NDVI数据值转换到[-1,1]。

3.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S2具体为从每个合成期的质量评估文件中提取冰/雪和云覆盖像元，建立每个合成期的冰/雪和云掩膜文件。

4.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S3具体为利用每个合成期的冰/雪和云掩膜文件，对冰/雪和云掩膜文件相应的NDVI图像进行掩膜，并分别计算每个合成期的无云和无冰/雪覆盖像元的多年均值，用均值代替每个合成期冰/雪覆盖像元NDVI值。

5.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S4具体为利用初始化NDVI数据，生成年际的NDVI时间序列数据。

6.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S5预处理具体为投影转换，使土地覆盖数据NDVI数据具有相同的投影。

7.根据权利要求6所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述土地覆盖数据为矢量数据，并且步骤S5的所述预处理还包括矢量数据向栅格数据的转换。

8.根据权利要求6所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S6具体为基于土地覆盖数据及像元之间的相关性,对云覆盖的像元进行空域处理。

9.根据权利要求1所述的顾及空间信息的植被指数时间序列数据重建方法，其特征在于，所述步骤S7具体为采用样条插值法对时间序列中云覆盖的像元进行插值计算。