基于EMD-DTW的土壤水分观测数据异常值检测方法与流程

文档序号:26645694发布日期:2021-09-15 03:07
基于EMD-DTW的土壤水分观测数据异常值检测方法与流程
基于emd

dtw的土壤水分观测数据异常值检测方法
技术领域
1.本发明涉及土壤水分观测数据质量检测技术领域,更具体地说,它涉及基于emd

dtw的土壤水分观测数据异常值检测方法。


背景技术:

2.土壤水分观测数据与其他气象观测资料一样,在投入使用之前,应当先经过数据质量的检测,剔除和订正异常(可疑)数据,保证数据的可靠性。
3.目前对于土壤水分观测数据异常值检测的研究成果,通常是利用土壤中水分因子与水文常数之间的关系理论,并参考地面气象要素质量控制中的空间一致性、时间一致性和内部一致性等方法,以单位时间内数据变化率计算为主要研究目标,通过临近观测站之间数据变化趋势对比,人工设定变化率阈值范围以查找异常数据。但是,由于土壤水分受土壤质地、水文环境、气象条件和人类活动等多重因素影响,其变化规律具有极强的个性化特征,因而在现有条件下,无论是利用传统的统计回归方法还是采用临近站参考方法,都会造成数据特征值的丢失和变化率误判,也就无法准确实现异常数据检测。
4.因此,本发明旨在设计提供基于emd

dtw的土壤水分观测数据异常值检测方法,方法采用观测站点数据自检的方式,以自身数据序列为参照,有效提取土壤水分变化特征,以解决上述问题。


技术实现要素:

5.本发明的目的是针对上述技术问题,提供基于emd

dtw的土壤水分观测数据异常值检测方法,本发明以土壤水分观测站点自身数据序列为参照,通过提取土壤水分观测数据变化特性实现异常值检测,检测过程中不需本观测站点以外的他站数据作为参考,具有较高普适性。
6.本发明的上述技术目的是通过以下技术方案得以实现的:基于emd

dtw的土壤水分观测数据异常值检测方法,利用emd分解和重构操作,将土壤水分时间序列数据特有变化特征提取出来,对提取出来的变化特征进行dtw相似度的匹配和序列长度对齐,并通过计算变化幅度差异来确定异常点的位置,以实现土壤水分异常值检测,具体包括以下步骤:
7.s1、将某观测站可能存在异常值的一定时间段内连续的土壤水分观测数据作为待检序列x,长度为t,则待检序列为x=[x1,x2...,x
t
];将该站已经通过质量检验的土壤水分观测数据序列作为参照序列y,长度为s,则y=[y1,y2,...,y
s
];
[0008]
s2、采用emd算法对步骤s1中的待检序列x和参照序列y进行emd分解,分别得到本征模态函数重构序列imfs(x)和imfs(y);并对imfs(x)和imfs(y)进行首尾延展,分别形成序列c=[0,imfs(x),0]和序列q=[0,imfs(y),0],序列c的长度为m,m=t+2,序列q的长度为n,n=s+2;
[0009]
s3、采用动态时间规整算法dtw对步骤s2中的序列c和序列q进行dtw相似度匹配和序列长度对齐,并输出dtw对齐后的序列c’和序列q’,和
k为对齐后序列的长度;
[0010]
s4、计算序列c’和序列q’的幅值差序列d’,即:
[0011][0012]
其中,d’=(c
’‑
q’)的长度为k,diff(c,q)表示c和q围绕原点(纵坐标0值线)的振幅值的差;
[0013]
s5、设threshold为异常判定门限值,且threshold>0,并遍历幅值差序列中的各个时序点:
[0014][0015]
若则将该时序点标记为超限点,通过该超限点的下标i定位并输出检测待检序列为x的异常点x
i
‑1。
[0016]
进一步地,步骤s2中,利用emd算法将待检序列x和参照序列y分解出来的imf分量各自累加,分别获得的本征模态函数重构序列imfs(x)和imfs(y)。
[0017]
进一步地,步骤s4中diff(c,q)的计算规则为:
[0018][0019]
进一步地,步骤s5中所述的异常判定门限值threshold为非负值,所述异常判定门限值threshold的取值方法为:
[0020]
1)设待检序列x的均值t为待检序列x的长度,则待检序列x的标准差为:
[0021]
2)设参照序列y的均值s为参照序列y的长度,则参照序列y的标准差为:
[0022]
3)设即以序列x和序列y标准差之和平均值的3倍作为异常判定门限值。
[0023]
综上所述,本发明具有以下有益效果:
[0024]
1、本发明的方法全流程保留土壤水分量纲信息,无需归一化等预处理,数据展现始终在时间序列坐标系下,过程清晰;
[0025]
2、本发明的方法以土壤水分观测站点自我检查为实现条件,从站点数据自身提取
变化特征,无需他站参照和限制,也无需建立复杂的多因素分析模型,因而该方法具有较高的普适性。
附图说明
[0026]
图1是本发明实施例中的流程图;
[0027]
图2是本发明实施例中dtw行走路径示意图;
[0028]
图3是本发明实施例中土壤体积含水量时间序列图;
[0029]
图4是本发明实施例中时间序列的emd分解与重构图(a.待检序列x的emd分解b.参照序列y的emd分解c.重构和首尾延展形成序列c和q);
[0030]
图5是本发明实施例中dtw对齐和超限点图(a.dtw行走路径b.对齐序列c.超限点检查);
[0031]
图6是本发明实施例中土壤体积含水量异常点结果图。
具体实施方式
[0032]
以下结合附图1

6对本发明作进一步详细说明。
[0033]
实施例:基于emd

dtw的土壤水分观测数据异常值检测方法,利用emd分解和重构操作,将土壤水分时间序列数据特有变化特征提取出来,对提取出来的变化特征进行dtw相似度的匹配和序列长度对齐,并通过计算变化幅度差异来确定异常点的位置,以实现土壤水分异常值检测,如图1所示,具体包括以下步骤:
[0034]
s1、将某观测站可能存在异常值的一定时间段内连续的土壤水分观测数据作为待检序列x,长度为t,则待检序列为x=[x1,x2...,x
t
];将该站已经通过质量检验的土壤水分观测数据序列作为参照序列y,长度为s,则y=[y1,y2,...,y
s
];
[0035]
s2、采用emd算法对步骤s1中的待检序列x和参照序列y进行emd分解,分别得到本征模态函数重构序列imfs(x)和imfs(y)。并对imfs(x)和imfs(y)进行首尾延展,分别形成序列c=[0,imfs(x),0]和序列q=[0,imfs(y),0],序列c的长度为m,m=t+2,序列q的长度为n,n=s+2;
[0036]
s3、采用动态时间规整算法dtw对步骤s2中的序列c和序列q进行dtw相似度匹配和序列长度对齐,并输出dtw对齐后的序列c’和序列q’,和k为对齐后序列的长度;
[0037]
s4、计算序列c’和序列q’的幅值差序列d’,即:
[0038][0039]
其中,d’=(c
’‑
q’)的长度为k,diff(c,q)表示c和q围绕原点(纵坐标0值线)的振幅值的差;
[0040]
s5、设threshold为异常判定门限值,且threshold>0,并遍历幅值差序列中的各个时序点:
[0041]
[0042]
若则将该时序点标记为超限点,通过该超限点的下标i定位并输出检测待检序列为x的异常点x
i
‑1。
[0043]
其中,步骤s2中,emd(empirical mode decomposition)分解算法将原始时间序列逐级分解为若干个本征模态函数(intrinsic mode function,imf)和一个余项(residual,res),其分解过程为筛分(sifting),步骤如下:
[0044]
step 1设待分解的时间序列为x(t),长度为t,t=1,...,t;
[0045]
step 2初始化,令r
i
‑1(t)=x(t),i=1;
[0046]
step 3获得第i阶的imf:
[0047]
a、令h
k
‑1(t)=r
i
‑1(t),k=1;
[0048]
b、判断出h
k
‑1(t)的所有极大值和极小值点;
[0049]
c、通过三次样条插值分别对极大值和极小值点进行拟合,求出上下包络线e
+
(t)和e

(t),并计算出上下包络线之间的均值m
k
(t)及h
k
(t)=h
k
‑1(t)

m
k
(t);
[0050]
d、判断是否不大于给定的门限值(通常取值0.2),若不大于,则c
i
(t)=h
k
(t),进行step 4;否则,令k=k+1,并转到b。
[0051]
step 4计算余量r
i
(t)=r
i
‑1(t)

c
i
(t),并判断r
i
(t)是否为单调函数或者常数,若是,则res=r
i
(t),分解结束;若否,则令i=i+1,并转到step 3的a。
[0052]
时间序列x(t)经过分解后,可以表示为各imf分量和余项res之和,即
[0053][0054]
式(1)中,m为分解出来的imf的阶数。设
[0055][0056]
式(2)中,imfs为各阶imf分量累加形成的本征模态函数重构序列。
[0057]
步骤s2中利用emd算法将待检序列x和参照序列y分解出来的imf分量各自累加,分别获得的本征模态函数重构序列imfs(x)和imfs(y)。
[0058]
步骤s3中利用dtw(dynamic time warping)动态时间规整算法,是把时间规整和距离测度计算结合起来的一种非线性归正计算方法,其能够实现不同长度序列的相似度匹配,并且具有较好的鲁棒性。
[0059]
dtw算法的基本思路:给定时间序列c=[c1,c2,

,c
m
]和q=[q1,q2,

,q
n
],长度分别为m和n,则可以构造一个n
×
m的矩阵网格,此时从矩阵的左下角w1=(1,1)出发,每次行走一步,方向只能为自己的右、上或右上三个方向之一的连续格点,一直到右上角w
k
=(m,n)结束,其中k为行走步数,满足max(m,n)≤k≤m+n

1。
[0060]
在上述行走原则下,可获得多条行走路径,通过公式:
[0061][0062]
求得行走距离代价最小的那条路径,即为最优化路径。而原时间序列c被非线性缩放为c

,q被非线性缩放为q

,c

和q

的长度相同,都为k。式(3)中,d(w
k
)是w
k
=(i,j)时c
i
和q
j
之间的距离,通常用欧式距离计算。为直观表述上述算法思想,我们以c=[c1,c2,c3,c4,c5],q=[q1,q2,q3,q4,q5,q6]两个不等长的时间序列举例,他们的长度分别为5和6。
[0063]
假设,经过dtw动态规整计算得出的最优化路径为图2所示,可以看出从矩阵网格
左下角走到右上角,最终经过8步行走到达终点。行走完毕后,原c序列被非线性缩放为原q序列被非线性缩放为c

和q

的长度都为8,这样就实现了两个序列的最优化对齐。同时,缩放后序列c

和q

中的各个点完整保留位置信息,上标代表在新序列的位置,下标代表其所在原序列中的位置。
[0064]
步骤s4中diff(c,q)表示c和q围绕原点(纵坐标0值线)的振幅值的差,计算规则为:
[0065][0066]
步骤s5中所述的异常判定门限值为非负值,所述异常判定门限值的取值方法为:
[0067]
1)设待检序列x的均值t为待检序列x的长度,则待检序列x的标准差为:
[0068]
2)设参照序列y的均值s为参照序列y的长度,则参照序列y的标准差为:
[0069]
3)设:
[0070][0071]
即以序列x和序列y标准差之和平均值的3倍作为异常判定门限值。
[0072]
以下通过运行实例对本发明做进一步说明:
[0073]
选择山东省d8815郓城北站2020年雨季时段的10cm深度土壤体积含水量观测数据(数据单位:%)为例,对上文实施方式进行演示和验证。
[0074]
为方便展示,以该站8月单月观测数据为待检序列x,以6月和7月两月的历史观测数据为参照序列y,观测频次为每小时1次,并假设参照序列y已经通过质量检查。序列原始统计信息如表1所示:
[0075]
表1序列统计信息表
[0076][0077]
此外,如图3所示,其为绘制的该站的10cm深度土壤体积含水量时间序列曲线图,从图3中可以看出,体积含水量观测值在所示时段内起伏较大,待检序列x在8月15日之后有
一个短时急降和急升的变化过程。
[0078]
按照本发明的方法的步骤流程,将待检序列x和参照序列y分别进行emd分解,由图4可以看出,待检序列x被分解为imf1~imf5共计5个imf(如图4a所示),参照序列y被分解为imf1~imf8共计8个imf(如图4b)所示,将各自imf累加重构并做首尾延展后,形成新的时间序列c和q(如图4c所示)。
[0079]
将序列c和序列q进行dtw对齐,行走路径如图5a所示,行走完成后,形成对齐序列c

和q

(如图5b)。按照上述公式(4)的计算规则,我们可以计算出c

和q

序列的幅值差序列d’=(c
’‑
q’)(如图5c所示),并根据表1中x序列和y序列的标准差按照上述公式(5)获得门限值threshold=3*(3.08+6.29)/2=14.055。由图5c可以看出,序列d’=(c
’‑
q’)中部分时序点的幅值差大于门限值threshold,记录这些点为超限点,并根据超限点时序号位置信息,定位输出原始待检序列x中异常点的位置,对应信息如表2所示:
[0080]
表2超限点与异常点的对应信息
[0081][0082]
因此,由图6中所示,可以标记出d8815郓城北站2020年8月10cm深度土壤体积含水量(数据单位:%)序列中的异常点,验证了8月15日之后的急降与急升过程属于观测异常,流程结束。
[0083]
在本发明的上述实施例中,本发明的方法首先利用时间序列emd分解操作,将土壤水分观测数据特有变化特征提取出来,并对提取出来的变化特征进行了dtw相似度的匹配和序列长度对齐,最后通过计算变化幅度差异来确定异常点的位置。本发明的方法的整个流程能够保留土壤水分量纲信息,无需归一化等预处理,数据展现始终在时序坐标系下,因而计算过程清晰;同时,本发明以观测站点自我检查为实现条件,从站点数据自身提取个性变化特征,无需选择其他参考站,也无需考虑土壤水文物理和天气条件因素,因而该方法具有较高的普适性。
[0084]
本具体实施例仅仅是对本发明的解释,其并不是对本发明的限制,本领域技术人员在阅读完本说明书后可以根据需要对本实施例做出没有创造性贡献的修改,但只要在本发明的权利要求范围内都受到专利法的保护。
再多了解一些
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1