CN108229760B

CN108229760B - 面向灾害预测的不等间距时间序列异常趋势分析方法

Info

Publication number: CN108229760B
Application number: CN201810129484.8A
Authority: CN
Inventors: 李勤勇; 吕江花; 宋建功; 杜建海; 马世龙
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2018-02-08
Filing date: 2018-02-08
Publication date: 2021-12-03
Anticipated expiration: 2038-02-08
Also published as: CN108229760A

Abstract

本发明公开了一种面向灾害预测的不等间距时间序列异常趋势分析方法，属于灾害预测技术领域。首先获取所要研究的灾害的设定指标的动态观测数据，确定不等间距时间序列的最优间距；然后去除动态观测数据的噪声成分；最后对得到的动态观测数据进行异常趋势检测；本发明提出不等间距时间序列最优间距的概念及计算方法，计算简便，利用该最优间距将不等间距时间序列转换为等间距时间序列；使用本发明提出的方法可以快速地识别出不等间距时间序列中的异常趋势，辅助相关专家进行灾害预测。

Description

面向灾害预测的不等间距时间序列异常趋势分析方法

技术领域

本发明属于灾害预测技术领域，具体是一种面向灾害预测的不等间距时间序列异常趋势分析方法。

背景技术

自然灾害是全球共同面临的重大威胁，突发性强、破坏力大是它们的共同特点，难以准确预测并避免，往往造成重大的人口伤亡和经济损失，带来巨大而深远的社会影响。灾害预测是防灾减灾、灾害应急事业的重要组成部分和未来发展方向。相关研究和实践经验表明，有效地灾害预测能够在较大程度上避免自然灾害造成的重大损失：例如1975年2月发生在我国辽宁省海城、营口地区的里氏7.0级地震，由于中国科学家对该次地震进行了准确预测并及时发布了短临预报，灾区伤亡人口共计18308人，仅占灾区总人口的0.22％，死亡1328人，仅占灾区总人口的0.02％。

目前，基于数据分析的方法是灾害预测的科学手段之一。基于数据分析的灾害预测方法的基本假设是长时间持续观测的目标数据中隐藏着灾害发生的前兆特征，其基本原理是力图从目标数据的一般变化趋势中，通过若干数据分析的方法和技术，发现异常变化趋势，然后结合自然灾害发生的一般规律和特点进行灾害预测。长期以来，我国相关防灾减灾部门针对灾害预测建立了规模宏大的数据观测体系，获取了十分丰富的动态观测数据，为利用数据分析方法进行灾害预测提供了宝贵数据资源。

灾害预测领域中的动态观测数据具有观测持续时间长、数据量大的一般特点。与此同时，灾害预测领域中的动态观测数据还具有不同于其他领域观测数据的特点：一是从数据形式上看，动态观测数据具有强烈的不等间距性，为一种不等间距的时间序列；二是从数据内容上看，动态观测数据在观测过程中普遍受到不同程度地自然环境影响，因此数据中往往存在不同程度的噪声成分。动态观测序列具有强烈的不等间距性，因此传统时间序列分析方法难以直接对其进行处理，其存在的噪声成分又在一定程度上掩盖了真实变化趋势。因此，灾害预测领域动态观测数据的上述两个特点使得从动态观测数据中提取异常变化趋势进而进行灾害预测变得异常困难。

发明内容

为了解决上述问题，本发明提出一组面向灾害预测的不等间距时间序列异常趋势分析方法，为基于数据分析的灾害预测提供一种新颖的分析思路和手段。由于本发明方法是针对不等间距时间序列这一动态观测数据的本质特征提出的，因此能够应用于灾害预测的不同领域，具有较广泛的适用性。

本发明提供的面向灾害预测的不等间距时间序列异常趋势分析方法，实现步骤包括：

步骤1，获取所要研究的灾害的设定指标的动态观测数据；动态观测数据为不等间距时间序列；

步骤2，确定不等间距时间序列的最优间距φ_b，目的是使得不等间距时间序列在φ_b为间距分隔时，能在最大程度上保证每个间距内都存在观测值；

步骤3，去除动态观测数据的噪声成分，设得到的动态观测数据表示为不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，v_i为在t_i时刻的观测值；

步骤4，对步骤3得到的动态观测数据进行异常趋势检测；

获取NTS的速率趋势序列

求取设定长时间段的均值A：

求

的均方差σ_Δ：

服从正态分布N(A,σ_Δ)，速率

超出方差限2σ_Δ、2.5σ_Δ、3σ_Δ的概率分别为5％、1％、0.3％，如果超出界限则可能属于某种因素造成的异常趋势变化。

相对于现有技术，本发明优点与积极效果在于：

(1)灾害预测领域的动态观测序列多为不等间距时间序列，经典时间序列的相关处理方法难以直接应用到不等间距时间序列上，本发明专利在明确不等间距时间序列定义的基础上，提出不等间距时间序列最优间距的概念及计算方法，利用该最优间距将不等间距时间序列转换为等间距时间序列；

(2)传统去除时间序列中噪声成分的方法大多计算复杂，且不适合于灾害预测领域的特殊应用需求，本发明专利提出一种计算较为简便的适用于灾害预测领域的不等间距时间序列成分组成模型，并应用该模型配合相关方法去除了灾害预测领域不等间距时间序列中的噪声成分；

(3)传统时间序列异常趋势识别方法都是基于等间距时间序列的，且不适合于灾害预测领域的特殊需求，本发明专利提出一种面向灾害预测的不等间距时间序列异常趋势识别方法，使用本专利提出的方法可以快速地识别出不等间距时间序列中的异常趋势，辅助相关专家进行灾害预测。

附图说明

图1为本发明对动态观测序列最优间距的计算流程图；

图2为本发明对异常趋势检测算法中主要参数的示意图；

图3为本发明的异常趋势检测方法的执行流程图；

图4为测项汤家坪-水准-AB原始数据图；

图5为汤家坪-水准-AB的线性趋势图；

图6为汤家坪-水准-AB消趋势、消周期数据图；

图7为利用本发明异常趋势检测方法对汤家坪-水准-AB的异常趋势提取的效果图。

具体实施方式

下面结合附图和实施例来说明本发明的技术方案。

针对灾害预测的动态观测数据具有强烈的不等间距性，本发明提出一种面向灾害预测的不等间距时间序列异常趋势分析方法。具体包括：第一，明确不等间距时间序列的定义，在此基础上提出不等间距时间序列最优间距的概念，给出计算不等间距时间序列最优间距的方法；第二：提出一种灾害预测领域动态观测数据成分组成模型，提出去除动态观测数据中长期趋势、周期趋势等噪声成分的具体方法和步骤；第三，确定动态观测序列中的异常趋势判定标准，给出一种基于速率变化的不等间距时间序列异常趋势检验方法。由于本发明面向灾害预测领域，因此在本发明中将不加区分地使用“不等间距时间序列”和“动态观测序列”。下面首先对本发明中用到的基本概念进行定义和说明，然后对上述方法进行详细说明。

首先，说明本发明中关于不等间距时间序列的相关基本概念。

定义1：时间序列TS；

时间序列，或称等间距时间序列，是指将同一统计指标的数值按其发生的时间先后顺序排列而成的序列，表示为一系列有序的(时间,数值)键值对，即时间序列TS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，其中，对i＝1,...,n，t_i+1-t_i＝δ，常数δ称为间距，v_i称为时间序列在t_i时刻的观测值；n为该时间序列的长度，记为|TS|。

定义2：子序列TS_sub；

给定一个长度为n的时间序列TS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，TS的子序列TS_sub定义为一个长度为m(m＜n)的、从TS中连续获取的时间序列，即TS_sub＝TS_p+m-1，1＜p＜n-m+1。定义时间序列TS上的操作TS||·||为获取TS的子序列操作，·表示子序列应符合的条件或规则。

定义3：不等间距时间序列NTS；

给定序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，其中，对i＝1,...,n，定义t_i+1-t_i＝φ_i，φ_i称为元素(t_i,v_i)与(t_i+1,v_i+1)之间的间距，称该序列为不等间距时间序列，n为NTS的长度，记为|NTS|。

定义4：时间间隔Δ；

给定长度为n的不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，i＝1,...,n，

且i＜j，则NTS中观测值(t_j,v_j)以及(t_i,v_i)之间的时间间隔表示为

其中，

表示时间单位T选取函数。

定义5：时间窗口Window；

时间窗口Window定义为一个四元组，表示为

其中w_s表示该时间窗口的开始时间，w_e表示该时间窗口的结束时间，w_l表示该时间窗口的长度，简称窗长，

w_s,w_e均为时间变量，窗长

代表自然数。窗长的大小由窗口开始时间w_s、窗口结束时间w_e以及时间单位选取函数共同确定。

定义6：数据匹配Match；

给定不等间距时间序列NTS，其间距集合φ_NTS＝(φ₁,...,φ_n-1)，其最优间距为φ_b，以及时间间隔Δ，设观测值(t_l,v_l)、(t_p,v_p)∈NTS，且l＞p，观测值(t_p,v_p)与观测值(t_l,v_l)配对，当且仅当t_l-Δ-φ_b/2＜t_p＜t_l-Δ+φ_b/2且t_p＝min{dist(t_l-Δ,t_q)}，t_q∈NTS||(t_l-Δ-φ_b/2,t_l-Δ+φ_b/2)||。dist(t_l-Δ,t_q)表示t_q和t_l-Δ之间的距离。NTS||(t_l-Δ-φ_b/2,t_l-Δ+φ_b/2)||为获取NTS的子序列。

根据不等间距时间序列的定义(上述定义3)，长度为n的不等间距时间序列中各相邻元素的间距为Φ＝(φ₁,...,φ_n-1)。为了最大程度上保证后续异常趋势识别算法的有效性，本发明按照如下思路确定不等间距时间序列的最优间距：寻找一个尽量小的自然数

使得不等间距时间序列在以此为间距分隔时，能够在最大程度上保证每个间距内都存在观测值。最优间距取值较小保证了后续异常趋势识别算法在进行数据匹配时的识别准确度，尽量多的间距内存在观测值则保证了后续异常趋势识别算法进行异常识别时的准确度。基于上述思路，本发明提出不等间距时间序列NTS在指定间距下的不缺数率的概念，基于不缺数率计算不等间距时间序列的最优间距。

设不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，间距

NTS以φ为间距时的不缺数率p_φ如下所示：

上述表示中，其中k_φ表示NTS以φ为步长时NTS中存在观测值的间距个数，q_φ为NTS以φ为步长时NTS中全部间距的个数。不缺数率p_φ刻画了不等间距时间序列以φ为间距时的数据缺失程度。

不等间距时间序列最优间距的计算步骤如图1所示，现对其进行详细描述：给定不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，开始日期start，结束日期end以及候选间距集合Φ＝{φ₁,φ₂,...,φ_m}，NTS的最优间距φ_b详细计算步骤如下：

步骤1.1，对每一个候选间距φ_i∈Φ进行处理，如果当前候选间距不为空，即读取成功，则继续执行步骤1.2；否则转步骤1.9；

步骤1.2，计算以φ_i为间距时NTS的不缺数率，设置日期指针ρ←start；初始化变量k_φ←0、q_φ←0；“←”表示赋值；

步骤1.3，如果ρ小于等于结束日期end，q_φ←q_φ+1，转步骤1.4，否则转步骤1.1；

步骤1.4，如果ρ+φ_i超过结束日期end，转步骤1.7；否则转步骤1.5；

步骤1.5，q_φ←q_φ+1，获取NTS中(ρ,ρ+φ_i)范围内的子序列NTS_sub＝NTS||(ρ,ρ+φ_i)||；如果NTS||(ρ,ρ+φ_i)||中存在观测值，则k_φ←k_φ+1，然后转步骤1.8执行；否则，转步骤1.6执行；

步骤1.6，日期指针以φ_i为步长向后滑动ρ←ρ+φ_i，转步骤1.3；

步骤1.7，q_φ←q_φ+1，获取NTS中(ρ,end)范围内的子序列NTS||(ρ,end)||，如果NTS||(ρ,end)||中存在观测值，则k_φ←k_φ+1，然后转步骤1.8执行；否则，转步骤1.6执行；

步骤1.8，计算NTS以φ_i为间距时的不缺数率

转步骤1.6；

步骤1.9，计算NTS的最佳间距

其中，argmin、argmax分别为求取使目标函数取得最小值或最大值时的变量值。

上面

是求取不缺数率最大的候选间距的i值，

是对所获得的i值选取其中的最小值，即从Φ中选取不缺数率最大且i值最小的φ_i作为NTS的最佳间距φ_b。

其次，本发明提出了一种灾害预测领域动态观测数据成分组成模型，并提出去除动态观测数据中长期趋势、周期趋势等噪声成分的方法。

根据经典时间序列分析理论中的Cramer分解定理，考虑面向灾害预测领域中动态观测数据的观测规律和特点，本发明将灾害预测领域不等间距时间序列的各种趋势变化归结为长期趋势、周期趋势、异常趋势和随机波动4大因素的综合影响，使用加法模型表示如下：

NTS＝T_t+S_t+A_t+I_t

上述表示中：

(1)T_t表示长期趋势(Trend)：由于灾害预测领域动态观测序列的观测时间跨度一般很长，可达十几年甚至几十年，因此观测序列在一定程度上受到长趋势(递增、递减等)的影响；

(2)S_t表示周期趋势(Season)：由于灾害预测领域动态观测序列的观测时间跨度一般很长，期间容易受到各种环境因素(如温度、降水等)的影响，因此，观测序列中往往存在一定的周期变化规律；

(3)A_t表示异常趋势(Abnormal)：这是面向灾害预测的时间序列分析中最关注的成分，序列中出现的异常点或异常波动往往预示着相关观测对象发生了异常现象，有可能是灾害发生的前兆；

(4)I_t表示随机波动(Immediate)：除了上述因素外，时间序列还受到各种其他因素的综合影响，这些影响导致序列呈现出一定的随机波动趋势。

按照上述建模设定，在对灾害预测领域的动态观测序列进行进一步分析以前，应首先将序列中的长期趋势成分以及周期趋势去除。下面分别描述去除长期趋势和周期趋势的方法。

设不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，如果该序列的长期趋势表现为线性特征，则考虑使用线性模型来对该序列进行拟合，将{t₁,...,t_n}作为自变量，{v₁,...,v_n}作为因变量，建立序列值随时间变化的回归模型，从而找到该序列中的线性成分，并将其去除。模型具体表示为：

上述模型中，{I_t}为随机波动，T_t＝a+bt就是消除随机波动影响之后该序列的长期趋势，{v_t-T_t}_t＝1,...,n为消除长程趋势后的序列。E()、Var()分别表示求均值和方差。

给定不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>，其最优间距为φ_b，去除NTS中周期趋势成分的步骤如下：

步骤2.1，使用φ_b计算周期个数及每个周期内的时期数，m＝φ_b，n＝12/m；其中，n为周期数，m为一个周期内的时期数；

步骤2.2，计算每个周期内的所有时期的平均值

v_ik表示第k个周期的第i个时期的观测值。

步骤2.3，对NTS中的每一个观测值(t_i,v_i)∈NTS，减去其对应的对应周期内的平均值

即得到去除周期趋势变化后的序列。

最后，在确定动态观测序列的最优间距φ_b以及去除序列中的噪声成分后，接下来对序列中存在的异常趋势进行提取，本发明提出一种基于速率变化的异常趋势检测方法。之所以选择速率作为表征序列异常趋势的物理量，一方面是因为速率的计算简单，另一方面，考虑到速率作为一种带有普遍意义的物理量，能够在较广泛的领域中进行应用。本发明提出的动态观测序列中异常趋势的检测方法，速率计算是其中的关键步骤，为了得到更准确的计算结果，本方法使用滑动窗口技术对得到的速率序列进行低通滤波。

本发明提供的面向灾害预测的不等间距时间序列异常趋势分析方法，首先获取所要研究的灾害的设定指标的动态观测数据，再对动态观测数据去除噪声，然后按照图3所示来对动态观测数据进行分析检测。设去噪后的不等间距时间序列NTS＝<(t₁,v₁),(t₂,v₂),…,(t_i,v_i),...,(t_n,v_n)>、开始时间start、结束时间end、窗长wlen、步长slen，时间间隔Δ，其最优间距为φ_b，假设NTS已去除噪声成分。如图3所示，本发明对动态观测序列中异常趋势检测的步骤具体如下：

步骤3.1，读取不等间距时间序列NTS；

步骤3.2，计算不等间距时间序列NTS的最优间隔φ_b；

步骤3.3，设置计算起始时间点ρ←start+wlen+Δ；

步骤3.4，如果ρ+wlen≤end，则转步骤3.5，否则转步骤3.15；

步骤3.5，生成时间窗口window←[ρ,ρ+wlen]；

步骤3.6，获取时间窗口中的子序列NTS_win＝NTS||[ρ,ρ+wlen]||；

步骤3.7，如果子序列NTS||[ρ,ρ+wlen]||中元素个数大于零，转步骤3.8；否则转步骤3.14；

步骤3.8，获取子序列NTS||[ρ,ρ+wlen]||中的元素(t_i,v_i)，如果(t_i,v_i)不为空，转步骤3.9；否则转3.13；

步骤3.9，利用最优间距φ_b和时间间隔Δ对(t_i,v_i)进行数据匹配，构造时间匹配区间L；

L←[t_i-Δ-φ_b/2,t_i-Δ+φ_b/2]；

步骤3.10，查找NTS中位于L中的子序列NTS_L，如果NTS_L中元素个数不为零，转步骤3.11；否则，更新i＝i+1，然后转步骤3.8；

步骤3.11，确定NTS_L中与(t_i,v_i)时间匹配的观测值(t_j,v_j)，匹配规则为

步骤3.12，计算(t_i,v_i)与(t_j,v_j)之间的速率r_i,j＝(v_i-v_j)/(t_i-t_j)*365，同时i＝i+1，转步骤3.8；

步骤3.13，计算窗口window内的平均速率(w_e,R_w)，R_w＝(∑r_i,j)/|NTS||window|||；

步骤3.14，滑动窗口ρ←ρ+slen，转步骤3.4；

步骤3.15，生成时间窗口window←[ρ,end]；

步骤3.16，获取时间窗口中的子序列NTS_win＝NTS||[ρ,end]||；

步骤3.17，如果子序列NTS||[ρ,end]||中元素个数大于零，转步骤3.18；否则转步骤3.24；

步骤3.18，获取子序列NTS||[ρ,end]||中的元素(t_i,v_i)，如果(t_i,v_i)不为空，转步骤3.19；否则转步骤3.23；

步骤3.19，利用最优间距φ_b和时间间隔Δ对(t_i,v_i)进行数据匹配，构造时间匹配区间L；

L←[t_i-Δ-φ_b/2,t_i-Δ+φ_b/2]；

步骤3.20，查找NTS中位于L中的子序列NTS_L，如果NTS_L中元素个数不为零，继续执行步骤3.21；否则，i＝i+1，转步骤3.18；

步骤3.21，确定NTS_L中与(t_i,v_i)时间匹配的观测值(t_j,v_j)，匹配规则为

步骤3.22，计算(t_i,v_i)与(t_j,v_j)之间的速率r_i,j＝(v_i-v_j)/(t_i-t_j)*365，同时i＝i+1，转步骤3.18；

步骤3.23，计算窗口window内的平均速率(w_e,R_w)，R_w＝(∑r_i,j)/|NTS||window|||；

步骤3.24，返回得到的NTS速率序列。

上述方法的处理对象为不等间距时间序列NTS，其最优间距根据本发明如图1所示最优间距计算方法确定为φ_b。由于实践中往往对目标序列部分时间段的观测值感兴趣，因此使用开始时间start和结束时间end获取NTS的一段子序列进行处理；使用滑动窗口技术对速率趋势序列进行平滑处理，窗长wlen和滑动步长slen用于控制窗口滑动过程中的窗口大小以及步长大小；时间间隔Δ用于确定进行速率计算的观测值之间的日期距离。上述变量的概念及其意义如图2所示，其中w₁,w₂...,w_n表示依次产生的滑动窗口。

在对动态观测序列的速率趋势计算完成后就可以进行异常趋势的提取了。

设

为速率趋势变化值，取一定长时间段的均值A：

并求

的均方差σ_Δ：

经过对大量数据进行χ²检验得知，

的基本上服从正态分布N(A,σ_Δ)。故速率

超出方差限2σ_Δ、2.5σ_Δ、3σ_Δ的概率分别为5％、1％、0.3％，如果超出界限而又查无其他原因，则可能属于某种因素造成的异常趋势变化。

下面以地震预测领域的跨断层流动形变资料分析为例，说明本发明方法的有效性和适用性。

地震来自断层的快速错动或是岩体的突然破裂，是一种以力学运动为主体的自然现象。实践表明，利用地壳形变观测结果来描述地壳的形变及其演化过程，对研究现今的地壳运动规律和提取与地震有关的前兆信息是有效的。下面以我国四川省汤家坪形变观测站的水准-AB观测结果为例，说明本发明方法的有效性和适用性。

汤家坪场地位于四川省宁南县汤家坪附近的黑水河北东岸(102.46N，27.0E)。地质构造为则木河断裂带紧靠南东的尾端地段，此地断裂走向略向东偏转，存在多组断裂和次级分支断裂。水准标志为铜质，基线墩为钢筋混凝土强制归心嵌铁标志，标石埋设深度为1.5米。该场地建成于1980年，始测于1981年9月，已经连续累积了30年的观测资料，观测精度较稳定，测站高差中每千米偶然中误差平均为±0.2mm，观测资料质量良好，该场地水准-AB测项截止至2011年9月的观测数据如下所示，原始数据时序图如图4所示。

表1汤家坪-水准-AB动态观测序列

从表1中可以看出，该测项始测于1981年9月18日，截止至2013年9月15日，共计160组观测值，最小值为(200+-7-15,-3764.70)，最大值为(1993-5-16,-3737.72),平均值A＝-3751.905，标准差σ＝7.32，最大值与平均值之比为0.99，说明在目前的观测区间内，最大值小于其平均值，因此难以判定该场站是否出现异常趋势。下面使用本发明方法对该测项的异常趋势进行提取。

按发明内容中所述的方法，第1步是确定动态观测序列的最优观测周期。给定备选最佳间距为Φ＝{1,2,3,6,12}，单位为月，利用发明内容中动态观测序列最优间距的计算方法，计算得到各备选间距的不缺数率分别为：p₁＝50.35％，p₂＝98.04％，p₃＝100％，p₆＝100％，p₁₂＝100％，可以看出，当φ_b＝2时观测数据可以取到除100％外最高的不缺数率，因此将这组数据的最佳间距定为2(月)。

第2步，去除该动态观测序列中的长程趋势和周期趋势。该测项的线性趋势如图5所示，从图中可以看出，该测项具有明显的递减趋势，对其进行线性显著性分析，在显著性水准α＝0.05的条件下，检验值p_α＝2.2e-16＜＜α，因此可以进一步认定该序列具有明显的线性趋势；利用发明内容中的动态观测序列噪声成分去除方法去除该测项中存在的噪声成分后，结果如图6所示。从该图中已经可以看出，经过上述处理后，该测项在2013年左右存在超过其应在范围的异常变化趋势。为了更加全面的分析，下面再利用本发明提出的异常趋势检验方法对该测项可能存在的异常趋势进行进一步验证。

第3步，使用异常趋势检验算法对去除噪声成分后的动态观测序列进行处理，得到动态观测序列的变化速率。本例中异常趋势检验算法的参数设置如表2所示。按照该参数设置对目标观测序列进行处理后，异常趋势检验算法得到的断层速率变化趋势特征值为：平均速率为0.58mm/a，均方差3.87mm/a，最大速率17.14mm/a，最小速率-5.23mm/a，异常趋势检验算法得到的处理结果如图7所示。

表2异常趋势检验算法参数取值

参数名	参数值
		开始时间start	1985-1-1
结束时间end	2013-9-30
		窗长w	6
步长s	3
		时间间隔Δ	12
最佳间距φ<sub>b</sub>	2

第4步，利用异常判定范围确定异常趋势。按照本发明的异常趋势判定方法，汤家坪-水准-AB的速率变化异常趋势判定标准为8.2mm/a，超过该数值的观测值都可以被判断为异常趋势。应用异常趋势检验算法后，汤家坪-水准-AB在芦山地震发生前的异常变化趋势被充分提取出来：该测项的速率变化在2010～2012年之间快速增加，2012至芦山地震发生前又发生快速下降，其最大值17.14mm/a超过异常判定预置8.2mm/a 2.1倍，因此可以判定该测项在芦山地震发生前已经存在异常变化。

在原始数据、预处理数据显示出不明显地异常趋势变化的情况下，本发明提出的方法能够从目标数据中提取出显著的异常变化趋势，说明本发明对灾害预测领域动态观测序列异常趋势识别的有效性和适用性。

Claims

1.一种面向灾害预测的不等间距时间序列异常趋势分析方法，其特征在于，实现步骤如下：

确定最优间距φ_b的方法如下：设候选间距集合Φ＝{φ₁,φ₂,...,φ_m}，设开始日期为start，结束日期为end，执行下面步骤：

步骤1.1，读取候选间距φ_i∈Φ，判断是否读取成功，若是继续步骤1.2，否则转步骤1.9；

步骤1.2，设置日期指针ρ←start；初始化变量k_φ←0、q_φ←0；“←”表示赋值；

步骤1.3，如果ρ小于等于结束日期end，变量q_φ自增1，继续步骤1.4；否则，转步骤1.1执行；

步骤1.4，如果ρ+φ_i超过结束日期end，转步骤1.7；否则，继续执行步骤1.5；

步骤1.5，变量q_φ自增1，获取NTS中(ρ,ρ+φ_i)范围内的子序列NTS||(ρ,ρ+φ_i)||；如果NTS||(ρ,ρ+φ_i)||中存在观测值，则变量k_φ自增1，然后转步骤1.8执行；否则，继续执行步骤1.6；NTS为不等间距时间序列；

步骤1.6，日期指针以φ_i为步长向后滑动，更新ρ增加φ_i，转步骤1.3执行；

步骤1.7，变量q_φ自增1，获取NTS中(ρ,end)范围内的子序列NTS||(ρ,end)||，如果NTS||(ρ,end)||中存在观测值，则变量k_φ自增1，然后转步骤1.8执行；否则，转步骤1.6执行；

步骤1.8，计算NTS以φ_i为间距时的不缺数率

然后转步骤1.6；

步骤1.9，从Φ中选取不缺数率最大且i值最小的φ_i作为NTS的最佳间距φ_b，计算如下：

其中，argmin、argmax分别为求取使目标函数取得最小值或最大值时的变量值；

步骤4，对步骤3得到的动态观测数据进行异常趋势检测；

获取NTS的速率趋势序列

求取设定长时间段的均值A：

求

的均方差σ_Δ：

服从正态分布N(A,σ_Δ)，速率

2.根据权利要求1所述的方法，其特征在于，所述的步骤3中，将动态观测数据表示为长期趋势T_t、周期趋势S_t、异常趋势A_t和随机波动I_t之和；去除动态观测数据的噪声成分，包括去掉其中的长期趋势和周期趋势成分；

(1)如果长期趋势表现为线性特征，将动态观测数据的不等间距时间序列中的时间作为自变量，观测值作为因变量，建立序列值随时间变化的回归模型，然后从中找到序列中的线性成分去除；

(2)去除周期趋势成分的方法是：

首先，使用φ_b计算周期个数及每个周期内的时期数，m＝φ_b，n＝12/m；其中，n为周期数，m为一个周期内的时期数；

其次，计算每个周期内的所有时期的平均值

v_ik表示第k个周期的第i个时期的观测值；

最后，对动态观测数据中的每一个观测值，减去其所对应的周期内的平均值

去除周期趋势成分。

3.根据权利要求1所述的方法，其特征在于，所述的步骤4中，获取NTS的速率趋势序列，具体方法如下：

步骤3.1，设置计算起始时间点ρ←start+wlen+Δ；start为开始时间，wlen为设定的窗长，Δ为设定的时间间隔；

步骤3.2，如果ρ+wlen≤end，执行步骤3.3，否则转步骤3.13；end为结束时间；

步骤3.3，生成时间窗口window←[ρ,ρ+wlen]；

步骤3.4，获取NTS位于时间窗口window中的子序列NTS_win；

步骤3.5，如果子序列NTS_win中元素个数大于零，执行步骤3.6；否则转步骤3.12；

步骤3.6，获取子序列NTS_win中的元素(t_i,v_i)，如果(t_i,v_i)不为空，执行步骤3.7；否则转步骤3.11；

步骤3.7，利用最优间距φ_b和时间间隔Δ对(t_i,v_i)进行数据匹配，构造时间匹配区间L；

L←[t_i-Δ-φ_b/2,t_i-Δ+φ_b/2]；

步骤3.8，查找NTS中位于L中的子序列NTS_L，如果NTS_L中元素个数不为零，执行步骤3.9；否则转步骤3.6；

步骤3.9，确定NTS_L中与(t_i,v_i)时间匹配的观测值(t_j,v_j)，匹配规则为

|t_j-t_i|表示t_j-t_i的时间长度；

步骤3.10，计算(t_i,v_i)与(t_j,v_j)之间的速率r_i,j＝(v_i-v_j)/(t_i-t_j)*365，同时i自增1，转步骤3.6；

步骤3.11，计算窗口window内的平均速率(w_e,R_w)，R_w＝(∑r_i,j)/|NTS||window|||；w_e表示时间窗口的结束时间；

步骤3.12，更新滑动窗口，ρ←ρ+slen，转步骤3.2；slen为设置的步长；

步骤3.13，生成时间窗口window←[ρ,end]；

步骤3.14，获取时间窗口window中的子序列NTS_win；

步骤3.15，如果子序列NTS_win中元素个数大于零，转步骤3.16；否则转步骤3.22；

步骤3.16，获取子序列NTS_win中的元素(t_i,v_i)，如果(t_i,v_i)不为空，转步骤3.17；否则转步骤3.21；

步骤3.17，利用最优间距φ_b和时间间隔Δ对(t_i,v_i)进行数据匹配，构造时间匹配区间L；

L←[t_i-Δ-φ_b/2,t_i-Δ+φ_b/2]；

步骤3.18，查找NTS中位于L中的子序列NTS_L，如果NTS_L中元素个数不为零，继续执行步骤3.19；否则，i自增1，转步骤3.16；

步骤3.19，确定NTS_L中与(t_i,v_i)时间匹配的观测值(t_j,v_j)，匹配规则为

步骤3.20，计算(t_i,v_i)与(t_j,v_j)之间的速率r_i,j，同时i自增1，转步骤3.16；

步骤3.21，计算窗口window内的平均速率(w_e,R_w)；

步骤3.22，返回得到的NTS速率趋势序列。

4.根据权利要求3所述的方法，其特征在于，所述的数据匹配的定义为：设观测值(t_l,v_l)、(t_p,v_p)∈NTS，且l＞p，观测值(t_p,v_p)与观测值(t_l,v_l)配对，当且仅当t_l-Δ-φ_b/2＜t_p＜t_l-Δ+φ_b/2且t_p＝min{dist(t_l-Δ,t_q)}，t_q∈NTS||(t_l-Δ-φ_b/2,t_l-Δ+φ_b/2)||，其中，dist(t_l-Δ,t_q)表示t_q和t_l-Δ之间的距离。