CN108388866A - 一种植被单调变化趋势检测方法及相关装置 - Google Patents
一种植被单调变化趋势检测方法及相关装置 Download PDFInfo
- Publication number
- CN108388866A CN108388866A CN201810161681.8A CN201810161681A CN108388866A CN 108388866 A CN108388866 A CN 108388866A CN 201810161681 A CN201810161681 A CN 201810161681A CN 108388866 A CN108388866 A CN 108388866A
- Authority
- CN
- China
- Prior art keywords
- time sequence
- trend
- pixel
- ndvi
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000012544 monitoring process Methods 0.000 title abstract 3
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 19
- 238000003860 storage Methods 0.000 claims abstract description 11
- 230000002159 abnormal effect Effects 0.000 claims description 41
- 238000001514 detection method Methods 0.000 claims description 28
- 230000008859 change Effects 0.000 claims description 27
- 238000001914 filtration Methods 0.000 claims description 22
- 238000012360 testing method Methods 0.000 claims description 21
- 238000012545 processing Methods 0.000 claims description 13
- 238000004590 computer program Methods 0.000 claims description 10
- 238000012800 visualization Methods 0.000 claims description 8
- 238000007689 inspection Methods 0.000 claims description 5
- 238000010998 test method Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 230000000630 rising effect Effects 0.000 claims description 3
- 230000003247 decreasing effect Effects 0.000 abstract description 3
- 230000000694 effects Effects 0.000 abstract description 3
- 230000000875 corresponding effect Effects 0.000 description 14
- 230000008569 process Effects 0.000 description 11
- 230000001932 seasonal effect Effects 0.000 description 9
- 230000006870 function Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000007423 decrease Effects 0.000 description 4
- 238000009826 distribution Methods 0.000 description 4
- 230000001174 ascending effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000015556 catabolic process Effects 0.000 description 2
- 238000006731 degradation reaction Methods 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 238000011497 Univariate linear regression Methods 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000002146 bilateral effect Effects 0.000 description 1
- 229910002092 carbon dioxide Inorganic materials 0.000 description 1
- 239000001569 carbon dioxide Substances 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000003306 harvesting Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/188—Vegetation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/2433—Single-class perspective, e.g. one-against-all classification; Novelty detection; Outlier detection
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Multimedia (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种植被单调变化趋势检测方法，利用EMD分解方法对每个像元重构后的NDVI时间序列进行分解得到趋势分量，进而对趋势分量进行单调性检验得到每个像元的趋势结果。不依赖一阶回归这种必须用最小二乘求解的技巧，且不管季节的趋势如何或局部的趋势如何，一定能获得整个序列单调增或单调减的变化趋势。本发明还提供了一种植被单调变化趋势检测系统、装置及计算机可读存储介质，同样可以实现上述技术效果。
Description
技术领域
本发明涉及遥感数据数据分析领域，更具体地说，涉及一种植被单调变化趋势检测方法、系统、装置及计算机可读存储介质。
背景技术
遥感长时间序列数据提供了关于全球地表动态变化以及变化趋势等重要信息，这些信息的有效提取可以为生产、社会发展、环境保护等各项重要决策提供依据。特别是植被的长期变化信息是气候变化中的关键因素，植被作为陆地生物圈的主要组成部分，植被变化趋势被认为与土地退化关联，植被状态常用于评估自然和农用土地的生产率和退化。
在几十年的时间框架内，植被的突然减少，一般认为是由一些短期过程导致的，如火灾、农作物收获或灾害等，植被的突然增加则被认为可能由降雨事件或雪盖的减少引起；植被的渐变则认为是体现了植被对全球变化的适应过程，如大洋振荡、持续的气候变化、年际降雨减少或大气中二氧化碳浓度增加等。植被的渐变，单调地变得更绿或更不绿，这种单调变化的趋势也称为“绿化”或“褐化”。
遥感数据中表征植被信息的特征是植被指数(NDVI，Normalized DifferenceVegetation Index)。植被“绿化”或“褐化”的趋势可以通过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(Empirical ModeDecomposition，经验模态分解)分解，将重构NDVI时间序列分解为有限个本征模函数(IMF，intrinsic mode function)和残差的叠加，每个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的具体计算方法如下：
其中，表示第k+1次迭代后生成的新的时间序列 表示第二时间序列w(yi)表示第i个数据点的质量函数权值，w(yi)在所有野点定义为零，在正常点定义为1。
下面对本发明实施例提供的一种植被单调变化趋势检测系统进行介绍，下文描述的一种植被单调变化趋势检测系统与上文描述的一种植被单调变化趋势检测方法可以相互参照。
图5为本发明实施例提供的一种植被单调变化趋势检测系统的结构框图，参照图5，本发明实施例提供的一种植被单调变化趋势检测系统装置可以包括：
提取模块401，用于提取NDVI时间序列图像中每个像元的时间序列；
重构模块402，用于对每个像元的所述时间序列进行重构得到每个像元的重构NDVI时间序列；
分解模块403，用于对每个像元的所述重构NDVI时间序列进行EMD分解，得到对应每个像元的所述重构NDVI时间序列的趋势分量；
检验模块404，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
本实施例的一种植被单调变化趋势检测系统用于实现前述的一种植被单调变化趋势检测方法，因此植被单调变化趋势检测系统中的具体实施方式可见前文中的植被单调变化趋势检测方法的实施例部分，例如，提取模块401，重构模块402，分解模块403，检验模块404，分别用于实现上述植被单调变化趋势检测方法中步骤S101，S102，S103和S104，所以，其具体实施方式可以参照相应的各个部分实施例的描述，在此不再赘述。
下面对本发明实施例提供的一种植被单调变化趋势检测装置进行介绍，下文描述的一种植被单调变化趋势检测装置与上文描述的一种植被单调变化趋势检测方法可以相互参照。
本发明实施例提供的一种植被单调变化趋势检测装置具体包括：
存储器，用于存储计算机程序；
处理器，用于执行所述计算机程序时实现如上述实施例所述植被单调变化趋势检测方法的步骤。
下面对本发明实施例提供的一种计算机可读存储介质进行介绍，下文描述的一种计算机可读存储介质与上文描述的一种植被单调变化趋势检测方法可以相互参照。
本发明实施例提供的一种计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如上述实施例所述植被单调变化趋势检测方法步骤。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
Claims (10)
1.一种植被单调变化趋势检测方法，其特征在于，包括：
提取NDVI时间序列图像中每个像元的时间序列；
对每个像元的所述时间序列进行重构得到每个像元的重构NDVI时间序列；
对每个像元的所述重构NDVI时间序列进行EMD分解，得到对应每个像元的所述重构NDVI时间序列的趋势分量；
对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
2.根据权利要求1所述的方法，其特征在于，所述对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果之后，还包括：
将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图。
3.根据权利要求2所述的方法，其特征在于，所述将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图，包括：
将每个像元的上升趋势使用绿色调表示绿化，下降趋势使用黄色调表示褐化，利用色调饱和度表示趋势结果的趋势显著程度。
4.根据权利要求1所述的方法，其特征在于，所述对每个像元的所述时间序列进行重构得到每个像元的重构NDVI时间序列，包括：
S301，确定每个像元的所述时间序列中的异常点；
S302，将每个像元的时间序列中的每个所述异常点更换为正常点，得到每个像元更新后的NDVI时间序列，作为第一时间序列；
S303，利用第一S-G滤波获得的缓变曲线与每个所述第一时间序列进行对比，确定每个第一时间序列中在相同时间点低于所述缓变曲线值的点作为异常值，并将每个所述异常值更换为所述缓变曲线中与所述异常值对应的时间点的值，得到每个更新后的NDVI时间序列，作为第二时间序列；
S304，利用第二S-G滤波对每个所述第二时间序列进行滤波得到滤波后的第二时间序列；
S305，利用每个所述第一时间序列与对应的每个所述滤波后的第二时间序列计算对应每个第一时间序列与第二时间序列的拟合残差指数，将本次迭代后的残差指数小于上次迭代后的残差指数的时间序列作为第一时间序列，返回S303；将本次迭代后的残差指数不小于上次迭代后的残差指数的时间序列作为重构NDVI时间序列。
5.根据权利要求4所述的方法，其特征在于，所述将每个像元的时间序列中的每个所述异常点更换为正常点，得到每个像元更新后的NDVI时间序列，作为第一时间序列，包括：
确定每个像元的时间序列中标记有云的第一数据点；
确定每个像元的时间序列中标记无云，且与相邻数据点的数值相差超过预设阈值的第二数据点；将所述第一数据点与所述第二数据点作为异常数据点；
判断所述异常数据点的相邻点是否为异常数据点；
若是，则将所述异常数据点的数据值更新为年内标记无云的点的数据值，或将所述异常数据点的数据值更新为其他年份同一时期的标记无云的数据点的值，得到更新后的NDVI时间序列，作为第一时间序列；
若否，则将所述异常数据点的值更新为相邻点的数据值，得到更新后的NDVI时间序列，作为第一时间序列。
6.根据权利要求1所述的方法，其特征在于，所述对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果，包括：
利用Mann-Kendall检验方法对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
7.一种植被单调变化趋势检测系统，其特征在于，包括：
提取模块，用于提取NDVI时间序列图像中每个像元的时间序列；
重构模块，用于对每个像元的所述时间序列进行重构得到每个像元的重构NDVI时间序列；
分解模块，用于对每个像元的所述重构NDVI时间序列进行EMD分解，得到对应每个像元的所述重构NDVI时间序列的趋势分量；
检验模块，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果。
8.根据权利要求7所述的系统，其特征在于，还包括：
可视化模块，用于对每个所述趋势分量进行单调性检验，得到每个像元的趋势结果之后，将所述每个像元的趋势结果进行可视化处理得到每个像元的单调性趋势图。
9.一种植被单调变化趋势检测装置，其特征在于，包括：
存储器，用于存储计算机程序；
处理器，用于执行所述计算机程序时实现如权利要求1至6任一项所述植被单调变化趋势检测方法的步骤。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至6任一项所述植被单调变化趋势检测方法步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810161681.8A CN108388866B (zh) | 2018-02-27 | 2018-02-27 | 一种植被单调变化趋势检测方法及相关装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810161681.8A CN108388866B (zh) | 2018-02-27 | 2018-02-27 | 一种植被单调变化趋势检测方法及相关装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108388866A true CN108388866A (zh) | 2018-08-10 |
| CN108388866B CN108388866B (zh) | 2021-08-24 |
Family
ID=63068706
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810161681.8A Active CN108388866B (zh) | 2018-02-27 | 2018-02-27 | 一种植被单调变化趋势检测方法及相关装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108388866B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109670467A (zh) * | 2018-12-25 | 2019-04-23 | 河南大学 | 一种基于sg滤波的森林变化快速识别方法 |
| WO2021022656A1 (zh) * | 2019-08-02 | 2021-02-11 | 湖南联智桥隧技术有限公司 | 一种数据处理方法 |
| CN113003033A (zh) * | 2021-02-19 | 2021-06-22 | 南京机电职业技术学院 | 基于StEMD_VGG的智能垃圾分类抓取机械手臂及控制方法 |
| CN113111799A (zh) * | 2021-04-19 | 2021-07-13 | 北华航天工业学院 | 一种基于集合经验模态分解的耕地土壤肥力水平监测方法 |
| CN113553549A (zh) * | 2021-07-26 | 2021-10-26 | 中国科学院西北生态环境资源研究院 | 一种植被覆盖度反演方法、装置、电子设备及存储介质 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102073869A (zh) * | 2010-12-27 | 2011-05-25 | 中国农业大学 | 基于点位数据和遥感影像数据的区域耕地质量监测方法 |
| CN103778241A (zh) * | 2014-02-10 | 2014-05-07 | 中国科学院南京地理与湖泊研究所 | 一种大尺度植被退化区域遥感识别方法 |
| US9292747B2 (en) * | 2013-03-15 | 2016-03-22 | The Boeing Company | Methods and systems for automatic and semi-automatic geometric and geographic feature extraction |
| CN106803059A (zh) * | 2016-12-02 | 2017-06-06 | 浙江工业大学 | 一种遥感植被指数时间序列森林监测方法 |
-
2018
- 2018-02-27 CN CN201810161681.8A patent/CN108388866B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102073869A (zh) * | 2010-12-27 | 2011-05-25 | 中国农业大学 | 基于点位数据和遥感影像数据的区域耕地质量监测方法 |
| US9292747B2 (en) * | 2013-03-15 | 2016-03-22 | The Boeing Company | Methods and systems for automatic and semi-automatic geometric and geographic feature extraction |
| CN103778241A (zh) * | 2014-02-10 | 2014-05-07 | 中国科学院南京地理与湖泊研究所 | 一种大尺度植被退化区域遥感识别方法 |
| CN106803059A (zh) * | 2016-12-02 | 2017-06-06 | 浙江工业大学 | 一种遥感植被指数时间序列森林监测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 陈鹏: "基于SPOT-VEG的内蒙古草原植被覆盖时空动态变化监测研究", 《中国优秀硕士学位论文全文数据库 农业科技辑》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109670467A (zh) * | 2018-12-25 | 2019-04-23 | 河南大学 | 一种基于sg滤波的森林变化快速识别方法 |
| WO2021022656A1 (zh) * | 2019-08-02 | 2021-02-11 | 湖南联智桥隧技术有限公司 | 一种数据处理方法 |
| CN113003033A (zh) * | 2021-02-19 | 2021-06-22 | 南京机电职业技术学院 | 基于StEMD_VGG的智能垃圾分类抓取机械手臂及控制方法 |
| CN113003033B (zh) * | 2021-02-19 | 2022-06-07 | 南京机电职业技术学院 | 基于StEMD_VGG的智能垃圾分类抓取机械手臂及控制方法 |
| CN113111799A (zh) * | 2021-04-19 | 2021-07-13 | 北华航天工业学院 | 一种基于集合经验模态分解的耕地土壤肥力水平监测方法 |
| CN113111799B (zh) * | 2021-04-19 | 2024-01-30 | 北华航天工业学院 | 一种基于集合经验模态分解的耕地土壤肥力水平监测方法 |
| CN113553549A (zh) * | 2021-07-26 | 2021-10-26 | 中国科学院西北生态环境资源研究院 | 一种植被覆盖度反演方法、装置、电子设备及存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108388866B (zh) | 2021-08-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108388866B (zh) | 一种植被单调变化趋势检测方法及相关装置 | |
| Zhu et al. | A changing-weight filter method for reconstructing a high-quality NDVI time series to preserve the integrity of vegetation phenology | |
| Adarsh et al. | Trend analysis of rainfall in four meteorological subdivisions of southern India using nonparametric methods and discrete wavelet transforms. | |
| Liu et al. | Global evaluation of gap-filling approaches for seasonal NDVI with considering vegetation growth trajectory, protection of key point, noise resistance and curve stability | |
| Ivatt et al. | Improving the prediction of an atmospheric chemistry transport model using gradient-boosted regression trees | |
| Wang et al. | Evaluating climate field reconstruction techniques using improved emulations of real-world conditions | |
| Chen et al. | Monitoring of rice cropping intensity in the upper Mekong Delta, Vietnam using time-series MODIS data | |
| CN108564002B (zh) | 一种遥感影像时间序列变化检测方法及系统 | |
| Araghi et al. | Associations between large-scale climate oscillations and land surface phenology in Iran | |
| CN116468958B (zh) | 通信铁塔安全检测方法及系统 | |
| Sur | An a-contrario approach to quasi-periodic noise removal | |
| Tang et al. | SURE-based optimum-length SG filter to reconstruct NDVI time series iteratively with outliers removal | |
| Deshcherevskii et al. | Iterative algorithm for time series decomposition into trend and seasonality: Testing using the example of CO2 concentrations in the atmosphere | |
| CN115019188B (zh) | 基于多源遥感数据的洲滩植被量动态解译方法及装置 | |
| CN113837042B (zh) | 一种基于局部均值分解的小波阈值自动化变形监测信号去噪方法 | |
| Dirkson et al. | Calibration of subseasonal sea‐ice forecasts using ensemble model output statistics and observational uncertainty | |
| CN104268838B (zh) | 一种面向超光谱数据库的小波去噪算法 | |
| Furby et al. | Continental scale land cover change monitoring in Australia using Landsat imagery | |
| Hong et al. | Haze removal for new generation optical sensors | |
| CN107506779B (zh) | 一种植物茎干含水量的估算方法及系统 | |
| Ahmad et al. | Haze effects on satellite remote sensing imagery and their corrections | |
| CN113031074B (zh) | 一种测井曲线消噪方法、装置、设备及存储介质 | |
| US11721097B2 (en) | Field image correlation differential change detection and alerting system | |
| CN114841899B (zh) | 时空频联合紧致编码去除红外图像横纹的方法及红外设备 | |
| US20210247297A1 (en) | Systems and Methods for Converting Satellite Images to Surface Reflectance Using Scene Statistics |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20230330 Address after: Hainan Tropical Ocean College National University Science Park Entrepreneurship Incubation Base, No.1 Yucai Road, Jiyang District, Sanya City, Hainan Province, 572022 Patentee after: Sanya huidao Technology Co.,Ltd. Address before: Hainan Institute of Tropical Oceanography, No.1 Yucai Road, Jiyang District, Sanya City, Hainan Province, 572022 Patentee before: HAINAN TROPICAL OCEAN University |