本发明涉及遥感数据数据分析领域，更具体地说，涉及一种植被单调变化趋势检测方法、系统、装置及计算机可读存储介质。
背景技术：
遥感长时间序列数据提供了关于全球地表动态变化以及变化趋势等重要信息，这些信息的有效提取可以为生产、社会发展、环境保护等各项重要决策提供依据。特别是植被的长期变化信息是气候变化中的关键因素，植被作为陆地生物圈的主要组成部分，植被变化趋势被认为与土地退化关联，植被状态常用于评估自然和农用土地的生产率和退化。
在几十年的时间框架内，植被的突然减少，一般认为是由一些短期过程导致的，如火灾、农作物收获或灾害等，植被的突然增加则被认为可能由降雨事件或雪盖的减少引起；植被的渐变则认为是体现了植被对全球变化的适应过程，如大洋振荡、持续的气候变化、年际降雨减少或大气中二氧化碳浓度增加等。植被的渐变，单调地变得更绿或更不绿，这种单调变化的趋势也称为“绿化”或“褐化”。
遥感数据中表征植被信息的特征是植被指数(ndvi，normalizeddifferencevegetationindex)。植被“绿化”或“褐化”的趋势可以通过ndvi的变化趋势表达。因此植被“绿化”或“褐化”信息可通过检测ndvi单调变化的趋势获得。一般的时间序列趋势检测常用的方法有一元回归变化斜率法、sen趋势度估计方法、曼肯德尔(mann-kendall)检验方法及季节性mann-kendall方法等。
ndvi序列不同于一般时间序列的特性有两点：强烈的季节性特点和数据在时间的有相关性特点。数据在时间上的相关性特点会使得由一元回归变化斜率法得到的ndvi趋势的有效性受到影响，因为这种相关性会破坏一元回归方法的假设条件。一元回归使用的假设条件通常是：回归变量是彼此独立的；回归残差有随机性，是零均值的；残差的方差对所有时间点是基本一样的。而sen趋势度估计方法、mann-kendall检验方法及季节性mann-kendall方法应用时也有一个隐含的假设，即每年每个季节或者每月变化的趋势必须是一致的，或者都是上升的，或者都是下降的，否则尽管每个季节有明显的趋势，但整体上却没有趋势。这一假设对具有强烈季节性特点的ndvi也是很难满足的。
因此，如何不依赖季节趋势、回归技巧获得整个事件序列趋势，是本领域技术人员需要解决的问题。
技术实现要素：
本发明的目的在于提供一种植被单调变化趋势检测方法、系统、装置及计算机可读存储介质，以不依赖季节趋势、回归技巧获得整个事件序列趋势。
为实现上述目的，本发明实施例提供了如下技术方案：
一种植被单调变化趋势检测方法，包括：
提取ndvi时间序列图像中每个像元的时间序列；
对所述每个像元的时间序列进行重构得到每个像元的重构ndvi时间序列；
对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量；
对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
其中，所述对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果之后，还包括：
将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图。
其中，所述将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图，包括：
将每个像元的上升趋势使用绿色调表示绿化，下降趋势使用黄色调表示褐化，利用色调饱和度表示趋势结果的趋势显著程度。
其中，所述对所述每个像元的时间序列进行重构得到每个像元的重构ndvi时间序列，包括：
s301确定每个像元的时间序列中的异常点；
s302，将每个像元的时间序列中的每个所述异常点更换为正常点，得到每个像元更新后的ndvi时间序列，作为第一时间序列；
s303，利用第一s-g滤波获得的缓变曲线与每个所述第一时间序列进行对比，确定每个第一时间序列中在相同时间点低于所述缓变曲线值的点作为异常值，并将每个所述异常值更换为所述缓变曲线中与所述异常值对应的时间点的值，得到每个更新后的ndvi时间序列，作为第二时间序列；
s304，利用第二s-g滤波对每个所述第二时间序列进行滤波；
s305，利用每个所述第一时间序列与对应的每个所述第二时间序列计算对应每个第一时间序列与第二时间序列的拟合残差指数，将本次迭代后的残差指数小于上次迭代后的残差指数的每个时间序列作为第一时间序列，返回s303；将本次迭代后的残差指数不小于上次迭代后的残差指数的每个时间序列作为重构ndvi时间序列。
其中，所述将每个像元的时间序列中的每个所述异常点更换为正常点，得到每个像元更新后的ndvi时间序列，作为第一时间序列，包括：
确定每个像元的时间序列中标记有云的第一数据点；
确定每个像元的时间序列中标记无云，且与相邻数据点的数值相差超过预设阈值的第二数据点；将所述第一数据点与所述第二数据点作为异常数据点；
判断所述异常数据点的相邻点是否为异常数据点；
若是，则将所述异常数据点的数据值更新为年内标记无云的点的数据值，或将所述异常数据点的数据值更新为其他年份同一时期的标记无云的数据点的值，得到更新后的ndvi时间序列，作为第一时间序列；
若否，则将所述异常数据点的值更新为相邻点的数据值，得到更新后的ndvi时间序列，作为第一时间序列。
其中，所述对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果，包括：
利用mann-kendall检验方法对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
为解决上述技术问题，本发明还提供了一种植被单调变化趋势检测系统，包括：
提取模块，用于提取ndvi时间序列图像中每个像元的时间序列；
重构模块，用于对所述每个像元的时间序列进行重构得到每个像元的重构ndvi时间序列；
分解模块，用于对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量；
检验模块，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
其中，还包括：
可视化模块，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果之后，将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图。
本发明还提供了一种植被单调变化趋势检测装置，包括：
存储器，用于存储计算机程序；
处理器，用于执行所述计算机程序时实现如权利要求1至6任一项所述植被单调变化趋势检测方法的步骤。
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至6任一项所述植被单调变化趋势检测方法步骤。
通过以上方案可知，本发明提供的一种植被单调变化趋势检测方法，包括：提取ndvi时间序列图像中每个像元的时间序列；对每个像元的所述时间序列进行重构得到每个像元的重构ndvi时间序列；对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量；对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
由此可见，本发明实施例提供的一种植被单调变化趋势检测方法，利用emd分解方法对每个像元重构后的ndvi时间序列进行分解得到趋势分量，进而对趋势分量进行单调性检验得到每个像元的趋势结果。不依赖一阶回归这种必须用最小二乘求解的技巧，且不管季节的趋势如何或局部的趋势如何，一定能获得整个序列单调增或单调减的变化趋势。本发明还提供了一种植被单调变化趋势检测系统、装置及计算机可读存储介质，同样可以实现上述技术效果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例公开的一种植被单调变化趋势检测方法流程图；
图2为本发明实施例公开的一种emd分解结果示意图；
图3为本发明实施例公开的一种具体的植被单调变化趋势检测方法流程图；
图4为本发明实施例公开的一种具体的emd分解流程图；
图5为本发明实施例公开的一种植被单调变化趋势检测系统结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例公开了一种植被单调变化趋势检测方法、系统、装置及计算机可读存储介质，以不依赖季节趋势、回归技巧获得整个事件序列趋势。
参见图1，本发明实施例提供的一种植被单调变化趋势检测方法，具体包括：
s101，提取ndvi时间序列图像中每个像元的时间序列。
具体地，首先在ndvi时间序列图像中提取出每个像元的时间序列，以对每个像元的时间序列进行重构。
s102，对每个像元的所述时间序列进行重构得到每个像元的重构ndvi时间序列。
在本方案中，在提取了每个像元的时间序列之后需要对其进行数据重构。数据重构的目的有两个，一个是滤波，从而消除或削弱因云或大气影响产生的随机噪声；另一个目的是通过重构使时间序列的时间维的采样间隔一致。
具体地，对时间序列进行重构可以包括两个部分，即野点处理和迭代逼近ndvi曲线上包络的过程。需要说明的是，在多光谱图像中被云覆盖像元的ndvi值会出现异常低的值，这些值即是野点；另外，大气效应也会引起ndvi值的异常变化，出现野点。野点处理就是尽可能滤除野点并相应的替换给出更合理的数值。而迭代逼近上包络则可以利用滤波方法进一步滤除异常值，以及可以滤除随机噪声。对于野点处理与迭代逼近上包络的具体步骤将在下述实施例做具体介绍，此处不再赘述。
s103，对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量。
具体地，对每个像元的重构后的ndvi时间序列均进行emd(empiricalmodedecomposition，经验模态分解)分解，将重构ndvi时间序列分解为有限个本征模函数(imf，intrinsicmodefunction)和残差的叠加，每个imf分量包含了每个原信号即重构ndvi时间序列的不同时间尺度的局部特征信号，而残差分量就是趋势项。
需要说明的是，emd是huang等人于1998年提出的一种不同于小波分析的新的时频信号分析方法，该方法属于自适应局部时频分析方法，它能根据数据信号本身的特性将其分解为有限个本征模函数imf和残差的叠加，每个imf分量包含了原信号的不同时间尺度的局部特征信号，非常适合非平稳信号的分析。自1998年以来，emd方法已经广泛应用于天气、地震、医疗等信号处理的各个领域。
emd的基本思想是将一个不规则信号表示成若干imf与单调残差函数相加的形式。对于自变量表示时间变化的一维信号，该单调残差函数就是趋势分量。一维信号x(t)的emd分解可表示为：
emd方法的本质是通过数据的时间尺度特征来获得本征波动模式，进行数据分解。这种分解过程可以形象地称之为“筛选(sifting)”过程，用x(t)代表原始数据信号，其分解步骤如下：
(1)求取原始数据信号即重构ndvi时间序列的极值点，包括极大值点和极小值点；
(2)分别由极大值点和极小值点利用三次样条插值函数拟合获得上包络线s1和下包络线s2；
(3)计算上下包络线的均值m1，m1＝(s1+s2)/2；
(4)将原数据序列即重构ndvi时间序列x(t)减去包络平均m1，得到一个新的数据序列h0，h0＝x(t)-m1；
(5)将新的序列重复步骤(1)到(4)进行迭代，直到满足迭代停止准则，获得第一个imf，imf1＝hk。
迭代停止准则通过如下公式计算：
当sdt超过给定的门限阈值时，迭代停止，此时hk代表第k次迭代后的新的序列。
参见图2，图2是emd分解结果的图示，从上到下分别是原始ndvi曲线，本征函数从imf1到imf6，以及残差也即趋势分量res。
s104，对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
需要说明的是，单调性检验是检验趋势是显著单调增、单调降还是没有显著变化。对每个趋势分量进行单调性检验就可以得到对应每个像元的趋势结果。
作为优选的可以利用mann-kendall(mk)趋势检验方法对趋势分量进行单调性检验，mann-kendall(mk)检验方法，最初由mann和kendall提出，用来检测水域降水的长期变化趋势和突变情况，后来广泛用于各类时间序列趋势分析中。mk检验方法不受少数异常值的干扰，也不需要样本遵循一定的分布，适用于非正态分布的数据。在mk检验中，设一时间序列数据(x1,x2,...,xn)是n个独立的、随机变量同分布的样本；检验的统计量s如下公式：
其中，
当n>10时，s近似正太分布，其均值为0，方差var(s)＝n(n-1)(2n+5)/18，标准的正态系统变量zmk通过以下公式进行计算：
在双边的趋势检验中，在α显著水平上，如果|zmk|＞z1-α/2则原假设是不可接受的，即在α显著水平上，时间序列数据存在明显的上升或下降趋势，否则接受原始时间序列无趋势的假设。统计量zmk大于0时是上升趋势；zmk小于0时是下降趋势。当显著性水平是α时，置信度是(1-α)100％。本发明要求95％的置信度，即zmk>1.96序列上升趋势明显，zmk<1.96序列下降趋势明显。
由此可见，本发明实施例提供的一种植被单调变化趋势检测方法，利用emd分解方法对每个像元重构后的ndvi时间序列进行分解得到趋势分量，进而对趋势分量进行单调性检验得到每个像元的趋势结果。不依赖一阶回归这种必须用最小二乘求解的技巧，且不管季节的趋势如何或局部的趋势如何，一定能获得整个序列单调增或单调减的变化趋势。
下面对本发明实施例提供的一种具体的植被单调变化趋势检测方法进行介绍，下文描述的一种具体的植被单调变化趋势检测方法与上文描述的实施例可以相互参照。
参考图3，本发明实施例提供的一种具体的植被单调变化趋势检测方法，具体包括：
s201，提取ndvi时间序列图像中每个像元的时间序列。
s202，对所述每个像元的时间序列进行重构得到每个像元的重构ndvi时间序列。
s203，对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量。
s204，对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
其中s201，s202，s203，s204可以参考上述实施例s101至s104，此处不再赘述。
s205，将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图。
具体的，为了更直观的了解每个像元的趋势结果，在本方案中得到趋势结果后对每个像元的趋势结果进行可视化处理，得到每个像元的单调性趋势图。
需要说明的是，成图时上升趋势与下降趋势可以用不同的色调标识，如上升趋势用绿色调表示“绿化”，下降趋势可以用黄色调表示“褐化”，趋势的显著程度可以利用色调的饱和度大小表示。
由此可见，在得到趋势结果后，对其进行可视化，并利用不同色调标识不同趋势，可以更直观的显示每个像元的趋势结果。
本发明实施例提供了一种具体的植被单调变化趋势检测方法，区别于上述实施例，本发明实施例对上述实施例的s102做了进一步的限定和说明，其他步骤内容与上述实施例大致相同，具体可以参考上述实施例，此处不再赘述。参考图4，s103具体包括：
s301，确定每个像元的时间序列中的异常点。
s302，将每个像元的时间序列中的每个所述异常点更换为正常点，得到每个像元更新后的ndvi时间序列，作为第一时间序列。
具体地，对数据重构的操作包括两个过程，即野点处理和迭代逼近ndvi曲线上包络的过程。野点处理是需要确定每个像元的时间序列中的异常点，然后将异常点更换为正常点。
需要说明的是，ndvi的野点都是明显低于正常ndvi值的点，因此需要滤除ndvi的异常值。需要处理的野点可以分为两类：一类是已经在ndvi数据的质量标记符中标记的野点，也就是已经确定为是野点的值；另一类则没有标明是否是野点，需要进一步的判断。
在进行野点处理时首先需要确定ndvi值中的野点，对于标记有云的数据点则将其直接确定为野点。对于没有标记的数据点，则需要进一步判断其是否为野点。判断过程具体为，将每个标记无云的数据点与相邻的数据点的值进行比较，判断其差值是否大于预设阈值，例如0.2，若大于，则说明为异常点，将其确定为野点。需要说明的是，在一个周期内，周期一般为20天，一个数据点的ndvi的变化不可能超过0.2，因此将判断的阈值设为0.2。通过判断确定的野点可以利用与该点相邻的点的值作线性插值，将该点的值进行替换。
对于有标记的点，例如标记有云的数据点，以相邻的标记无云的数据值代替。需要说明的是，对于标记有云的数据点，如果其相邻的数据点的值也标记有云时，则用年内标记无云的数据值代替。需要说明的是，如果有多年数据，则可以将多年数据中，同一时期的一个数据值代替。
s303，利用第一s-g滤波获得的缓变曲线与每个所述第一时间序列进行对比，确定每个第一时间序列中在相同时间点低于所述缓变曲线值的点作为异常值，并将每个所述异常值更换为所述缓变曲线中与所述异常值对应的时间点的值，得到每个更新后的ndvi时间序列，作为第二时间序列。
需要说明的是，s-g(savitzky-golay)滤波是一种局部(窗口或带宽)多项式逼近的方法，s-g滤波器不同带宽和多项式阶数的组合表示的不同含义，宽带低阶的s-g滤波得到的结果是比窄带高阶的s-g滤波更平滑。本发明在迭代过程中，两次使用s-g滤波。第一次使用低阶宽带的第一s-g滤波以便检出疑似异常值，将异常值用新的值替代后，再次使用高阶窄带的第二s-g滤波方法平滑滤除随机噪声且尽可能拟合数据。
迭代逼近ndvi曲线上包络的过程中，首先要进一步检测异常值。具体地，利用第一s-g滤波获得缓变曲线，同一时间低于缓变曲线的点即确定为意思的异常值，将异常值用新的值替代后得到更新后的ndvi时间序列，即第二时间序列。作为优选的，第一s-g滤波使用带宽为8、多项式p＝2的低阶带宽的s-g滤波。
需要说明的是，一般认为缓慢变化曲线体现了每年植被的循环，而缓慢变化的过程，云或者恶劣的大气条件造成的大多数噪声点的值应该低于缓慢变化的曲线。
当有高于缓变曲线的点的值时，则需要将这些异常值进行替换。由于要迭代逼近上包络线，因此将缓变曲线与第一时间序列进行比较，二者取同一时间点较大的值，将所有异常值更换后得到新的序列，即第二时间序列。
s304，利用第二s-g滤波对每个所述第二时间序列进行滤波得到滤波后的第二时间序列。
在本步骤中使用高阶窄带的第二s-g滤波方法平滑滤除随机噪声且尽可能拟合数据。作为优选的，第二s-g滤波选择的多项式阶数为6，带宽是4。
s305，利用每个所述第一时间序列与对应的每个所述滤波后的第二时间序列计算对应每个第一时间序列与第二时间序列的拟合残差指数，将本次迭代后的残差指数小于上次迭代后的残差指数的时间序列作为第一时间序列，返回s303；将本次迭代后的残差指数不小于上次迭代后的残差指数的时间序列作为重构ndvi时间序列。
每次迭代后，即执行s303至s304后，都需要计算拟合残差指数，并且将本次迭代的拟合残差指数与上一次迭代结算的拟合残差指数进行比较，当拟合残差指数不再下降时，就不再进行下一次迭代，将残差指数没有下降的时间序列作为重构ndvi时间序列，而残差指数相比上次迭代有下降的时间序列则返回s303，继续迭代。
其中，拟合残差指数为fk,fk的具体计算方法如下：
其中，
下面对本发明实施例提供的一种植被单调变化趋势检测系统进行介绍，下文描述的一种植被单调变化趋势检测系统与上文描述的一种植被单调变化趋势检测方法可以相互参照。
图5为本发明实施例提供的一种植被单调变化趋势检测系统的结构框图，参照图5，本发明实施例提供的一种植被单调变化趋势检测系统装置可以包括：
提取模块401，用于提取ndvi时间序列图像中每个像元的时间序列；
重构模块402，用于对每个像元的所述时间序列进行重构得到每个像元的重构ndvi时间序列；
分解模块403，用于对每个像元的所述重构ndvi时间序列进行emd分解，得到对应每个像元的所述重构ndvi时间序列的趋势分量；
检验模块404，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
本实施例的一种植被单调变化趋势检测系统用于实现前述的一种植被单调变化趋势检测方法，因此植被单调变化趋势检测系统中的具体实施方式可见前文中的植被单调变化趋势检测方法的实施例部分，例如，提取模块401，重构模块402，分解模块403，检验模块404，分别用于实现上述植被单调变化趋势检测方法中步骤s101，s102，s103和s104，所以，其具体实施方式可以参照相应的各个部分实施例的描述，在此不再赘述。
下面对本发明实施例提供的一种植被单调变化趋势检测装置进行介绍，下文描述的一种植被单调变化趋势检测装置与上文描述的一种植被单调变化趋势检测方法可以相互参照。
本发明实施例提供的一种植被单调变化趋势检测装置具体包括：
存储器，用于存储计算机程序；
处理器，用于执行所述计算机程序时实现如上述实施例所述植被单调变化趋势检测方法的步骤。
下面对本发明实施例提供的一种计算机可读存储介质进行介绍，下文描述的一种计算机可读存储介质与上文描述的一种植被单调变化趋势检测方法可以相互参照。
本发明实施例提供的一种计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如上述实施例所述植被单调变化趋势检测方法步骤。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。