基于时空Kriging改进的多云地区MODISNDVI时间序列重构方法

基于时空kriging改进的多云地区modis ndvi时间序列重构方法
技术领域
1.本发明属于遥感技术中植被指数时间序列重建领域，涉及一种基于时空kriging改进的ndvi时间序列重构方法。通过时空kriging插值补充多云地区缺失的ndvi信息，来重建高质量的ndvi时间序列。


背景技术：

2.归一化植被指数(ndvi)被定义为近红外波段与红光波段植被冠层反射率值之差和这两个波段数值之和的比值，它被广泛应用于研究植被分布情况、监测植被动态变化，以及预测环境变化带来的生态影响。长时间遥感观测获得的ndvi时间序列曲线，是反映植被生物学特征相随时间变化的最佳指示因子,也是季节变化和人为活动影响的重要指示器。目前ndvi时间序列产品主要由avhrr、modis、spot vgt、landsat tm/etm等卫星传感器提供，但由于受到恶劣大气、云、雪污染以及传感器误差的影响，ndvi观测数据通常会出现观测值降低或缺失等异常情况。因此，通过预处理对ndvi时间序列进行去噪和重构以提供更高质量的时序数据集，是植被指数研究的重要工作之一。
3.目前，国内外研究者已经提出了许多有效的方法来重建高质量的ndvi时间序列，主要分为纯时间域方法和时空方法两大类。其中纯时间域方法应用最广，如最大值合成法(mvc)、通过平滑函数拟合整个时间序列的ndvi曲线的非对称高斯函数拟合法(a-g)与双逻辑斯蒂拟合法(d-l)、通过滤波器滑动处理时间序列的局部窗口的s-g滤波、均值迭代滤波(mvi)。上述方法均能在一定条件下有效重构ndvi时间序列，但都存在各自的优缺点和局限性。例如，s-g滤波能够快速平滑时间序列并剔除异常值，在植被类型提取中表现出色，但其保真性较差，且存在边缘临界值错估的问题。a-g、d-l方法去噪能力良好，保真性强，但这两种方法参数设置十分复杂，不适用于有多个生长季的地区。此外这些纯时间域方法在数据出现连续缺失时表现较差，往往会出现错估。
4.而结合空间信息的时空方法主要包括预测缺失值的时空窗口法、利用土地覆盖信息的时空迭代法以及结合相似像元提供辅助信息的stsg方法。此类方法侧重于利用时空信息为时间序列重建提供参考，但另一方面由于忽略了长时间序列中植被的周期性特征，且无法实现时空信息的有效协同利用，这些时空方法在经常出现数据连续缺失的多云地区仍受到诸多限制。
5.kriging法是依据变异函数理论对区域内变化量进行线性无偏最优估计的方法，而时空kriging方法是在普通空间kriging变异函数的基础之上展开的时空扩展。该方法能够有效利用时空邻近的采样点进行插值，弥补了单一时刻有效样本不足带来的插值精度问题，在气温变化预测、pm2.5监测等多种领域都得到了成功应用。


技术实现要素：

6.针对以上方法提到的无法应对数据连续缺失以及无法有效利用时空信息等问题，
本方法通过引入结合质量信息的预插值以及时空kiging插值来分别提取长时间序列中的植被周期性特征以及邻近时空中的有效信息，可以有效解决多云地区ndvi时间序列中因云噪声所引起的连续数据缺失，同时在此基础之上开展的加权s-g迭代滤波能进一步消除噪声并平滑时间序列曲线，最终能够得到高质量的ndvi时间序列数据。
7.鉴于此，本发明采用的技术方案是：一种基于时空kriging改进的多云地区modis ndvi时间序列重构方法，包括以下步骤：
8.(1)对ndvi时间序列数据进行投影转换与图像裁剪预处理；采用mod13q1提供的pr数据和由vi qa提取出产品可用性vu数据转换为十进制后，可作为辅助质量信息标注有效、低质量和无效像元，为预插值处理提供参考依据。
9.(2)根据质量信息标注对时间序列进行预插值。
10.(3)将预插值后的ndvi时间序列数据集视为一个时空随机场，引入叠加窗口机制，借助空间探针对其附近搜索窗口内的所有像元做多时段筛选，得到满足拟合时间和空间变异函数所需数量的有效样本点。
11.(4)对数据进行stl时间序列分解以满足时空区域变化量二阶平稳的要求，将数据分解为趋势项t、季节项s和随机项r三部分并剔除季节项s，在此基础上开展变异函数构建与时空kriging插值，在插值过后加入季节项s得到最终的插值结果。
12.(5)进行基于质量权重的s-g迭代滤波，平滑时间序列曲线，得到高质量的ndvi时间序列曲线。
13.一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被执行时，可实现上述的基于时空kriging改进的多云地区modis ndvi时间序列重构方法。
14.本发明提出的一种基于时空kriging改进的重构方法。相比常规的时间域重构方法，引入了质量信息标注与预插值流程，能够有效从长时间序列中识别植被的周期性特征，使缺失数据得到初步恢复，同时又符合植被的生长变化趋势。
15.由于多云地区ndvi时间序列数据集存在严重的云噪声干扰与数据缺失问题，采用时空kriging插值对时序数据进行处理，可以有效利用缺失像元临近时空位置的其他未被云噪声干扰的正确像元的信息来补充缺失像元。通过实验发现，在预插值过后采用时空kriging插值能明显提高一些明显异常的低值，并保留了原始有效值，同时在空间信息影响下使某些连续值不再呈单纯的线性趋势，能更好地描述植被的生长变化。而时空kriging的插值精度与变异函数的构建直接相关，采用叠加窗口机制对样本点进行多时段有效筛选就是为了成功构建出一个最优的，最能描述随机场中各位置点之间的空间关系的变异函数，从而能够使用kriging方法对随机场中的待测点进行正确的预测与插值。
16.基于质量权重的s-g迭代滤波有别与传统的s-g滤波，采用了质量信息加权的滤波过程，能够为时间序列数据中不同质量的数据赋予不同的权重，更高质量的数据将占有更大的权重，而更低质量的数据占权重更小。在滤波过程中，高权重的数据会保留下来，低权重的数据会在迭代过程中更新，不断靠近整个时间序列曲线的上包络线，从而得到更加可靠、更加光滑的时间序列曲线。在时空kriging插值后会不可避免地出现一些误差，造成时间序列曲线出现波动和离群值的现象，而在插值后应用基于质量权重的s-g迭代滤波，可以修正这些误差，并使得整个时间序列曲线能够更好地拟合植物的真实生长变化趋势。
附图说明
17.为了更清楚地说明本发明的技术方案，下面将对实施例中所需要使用的附图作简单地介绍。
18.图1是本发明提出的基于时空kriging改进的多云地区modis ndvi时间序列重构方法技术路线图；
19.图2是结合质量信息的预插值流程中质量标注方法以及数据替换示意图；
20.图3是利用ndvi时间序列数据建立时空随机场示意图；
21.图4是利用时间切片与叠加窗口筛选机制筛选有效样本示意图。
具体实施方式
22.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、详细地描述。所描述的实施例仅仅是本发明的一部分实施例。
23.本发明拟解决实际应用中，多云地区modis ndvi时间序列存在数据连续缺失、常规方法无法有效利用时空信息等问题，提出一种基于时空kriging改进的重构方法。该方法基于以下2个假设：
①
空间上云噪声的存在不是连续的，且像元之间存在空间和时间自相关性，因此邻近时空未被云覆盖的像元可以对其他像元提供插值参考。
②
被云覆盖的像元在其他年份的同一日可能是清晰可用的，可以从长时间序列中提取出其他年份周期性特征相似的ndvi值作为参考。
24.图1所示的是本发明的具体技术路线，其中投影转换与图像裁剪为数据预处理步骤，pixel reliability、vi usefulness为modis提供的辅助数据，可以用于标注有效、低质量和无效像元。进一步的，结合质量信息的预插值首先需要对时间序列整体进行线性插值，以获得大致的植被生长趋势。以单个像元为例，在原始时间序列中的有效像元将在插值后保留其原始值。无效像元将以每年同一日期doy(the same day of years)为参考，经过测试，以无效像元为中心向两侧取出长度为5的窗口数据最能代表此像元周围的生长趋势，与其余年份doy有效值的窗口数据进行person相关性分析可以判断它们是否相似。若相关性系数大于0.9，说明与此doy生长趋势相似，无效值将被所有相似doy的平均值替换，若无相似doy，则无效值将由线性插值替换。低质量像元处理方法和无效像元相似，由于云雾噪声对ndvi的影响多为负偏差，因此只取高于原始值的结果替换。ndvi时序数据将通过结合质量信息的预插值初步恢复缺失的信息，具体质量标注方法与数据替换示意见图2。
25.预插值过后，将ndvi时间序列数据集视为一个时空随机场s
×
t，s表示空间域，t表示时间域。将z(s)定义为时空随机场s
×
t(s∈s)中空间点s处的空间探针，有z(s)＝{z(s,t1),z(s,t2),
…
,z(s,tn)}，t1...tn代表不同时刻t。z(s)可表示空间位置上任意一像元的全部时间序列；将z(t)定义为时空随机场s
×
t(t∈t)中时间t下的时间切片，有z(t)＝{z(s1,t),z(s2,t),
…
,z(sn,t)}，z(t)可表示某日期下实验区域内所有像元的ndvi观测值，时空随机场的建立如图3。为满足时空kriging插值有效样本点需求，本发明使用一种叠加窗口筛选机制，借助空间探针对其附近搜索窗口内的所有像元做多时段筛选，得到满足拟合时间和空间变异函数所需数量的有效样本点。具体方法为：
①
筛选记录某一时间切片中所有在预插值阶段标注为有效的像元；
②
筛选邻近的2个时间切片内的有效像元，将其在原时间切片中对应位置的参考值加入记录，经过叠加以后的有效时间切片可以保证拟合空间变异函
数所需的有效样本数，同时降低时间变异带来的误差。
③
检查记录中有效样本点占窗口内总像元数量的百分比，若低于20％则说明此区域受云雾噪声影响严重，此时间切片无效，需对窗口内像元进行均匀采样，并加入参考值。叠加窗口筛选机制示意图如图4。
26.筛选出有效数据后，还需要对数据进行stl时间序列分解以满足时空区域变化量二阶平稳的要求。时间序列一般可以分解为趋势项t、季节项s和随机项r三部分，其中随机项r满足二阶平稳。由于趋势项对时间序列的平稳性影响较小，因此在时序分解后选择剔除季节项s，剩余两部分之和则满足平稳性要求，可以开展时空kriging插值。
27.对ndvi时间序列数据进行时空kriging插值，将以邻近时空有效像元提供参考信息，插值补充缺失像元。时空kriging的插值公式定义如下：
[0028][0029]
其中(si,ti)代表某样本点i在时空随机场中的坐标，z(si,ti)表示该样本点的观测值。λi代表其权重系数，且需满足条件n为样本点总数，z
*
(s,t)代表时空点(s,t)处的估计值。同时假设时空区域变化量z(s,t)满足二阶平稳的条件下，时空变异函数公式为：
[0030][0031]
其中(s,t)代表样本点在时空随机场中的位置，h＝(hs,h
t
)代表相应的时空距离，γ(hs,h
t
)为变异函数，式中e、var表示数学期望和方差运算。
[0032]
本文发明使用一种常用的积和式模型来构建时空变异函数，其构建方法如下：
[0033]
γ(hs，h
t
)＝(k1c
t
(0)+k2)γs(hs)+(k1cs(0)k3)γ
t
(h
t
)-k1γs(hs)γ
t
(h
t
)
[0034][0035]
其中γ(hs,h
t
)、γs(hs)、γ
t
(h
t
)分别对应时空、空间、时间变异函数；c
st
(0,0)、cs(0)、c
t
(0)分别对应时空、空间、时间变异函数基台值。
[0036]
根据估计方差最小原则，可以根据积合式模型公式求出任意两点的时空变异函数值，通过在时空变异函数构建公式中引入拉格朗日乘除数可求得最优权重系数λi(i＝1,2,...,n)，带入插值公式即可对任意待估点z
*
(s,t)进行插值。完成时空kriging插值后还需要，加上之前去除的季节项才能得到最终结果。
[0037]
进一步的，为处理插值所带来的误差，更好地降低噪声，且使时序曲线能够正确地反映植被生长变化情况，将在插值的基础之上进行加权s-g迭代滤波。将不同质量的数据赋予不同的权重，质量可靠的数据将拥有更大的权重，可以作为控制点对植被的上包络线框架进行整体上的控制，通过多次迭代后拟合曲线将更为平滑，能够更有效地反映植被物候生长状态，且能够最大程度地保留有效信息。其滤波的基本公式为：
[0038][0039]
其中y
j+i
，表示原始时间序列值和滤波结果，wi表示相应的s-g滤波权重，n为滤波器长度，m为半窗口宽度。参考质量标注和两次插值过程其初始权重设置思路为：有效值始终是可信的，设置权重为1；无效值往往缺失严重，经过插值后误差较大，设置权重为0.6；低质量值只受到噪声一定程度的影响，权重设置为0.8。
[0040]
然后将滤波半窗口设置为9，多项式拟合阶数设置为6，迭代次数2次，结合表中权重wi对ndvi时间序列进行平滑滤波。对于wi＝1的点，两次迭代中值均不改变，其余值按下式更新权重：
[0041][0042]
其中wj代表更新后的权值，yj代表原始值，代表拟合过后的值。经过平滑滤波后，可以得到最终的高质量ndvi时序数据。
[0043]
在总共5种方法的噪声模拟实验中，本发明在各类植被类型上平均mae均为最低；在空间重建实验中rmse值、r2值、cc值均为最佳，整体重建效果较好；在不同植被类型测试点的时间序列重构实验中，本发明的拟合曲线也表现出更好的稳定性以及平滑度，说明该方法具有较强的数据保真性和去噪性能。
[0044]
以上这些实例应理解为仅用于解释本发明而不用于限制本发明的保护范围。依据本发明所作的各种改动或修改这类等效变化和修饰同样落入本发明权利要求所限定的范围。