CN106446548A - 一种基于小波分析的动力学结构突变检测方法 - Google Patents

Abstract

本发明公开了一种基于小波分析的系统动力学结构突变的检测方法，依据描述系统状态的时间序列，通过滑动移除技术从原序列中选取子序列，采用小波变换系数估计子序列的标度指数，依据各子序列标度指数趋势判断系统动力学结构发生突变点或突变区间，并对突变区间进行方差贡献分析以验证结果的可靠性。本发明的突变检测方法不仅能够对平稳序列进行检测，而且适用于非平稳序列，能够快速、有效和准确地判断系统动力学结构的突变情况。

Description

一种基于小波分析的动力学结构突变检测方法

技术领域

本发明涉及物理学领域，具体而言涉及一种动力学结构突变的检测方法。

背景技术

标度指数作为表征系统内部本质特征的一个重要参数，在地球物理学时间序列相关动力学进程的研究中有着重要的意义。如基于标度指数和滑动移除技术结合而提出的动力学结构突变检测方法：滑动移除去趋势波动分析、滑动移除重标极差分析、滑动移除重标方差分析等，这些方法相对于传统统计方法如滑动t检验、克拉默法、Yamamoto信噪比法等不仅能够有效地检测系统动力学结构突变的位置，而且能够对系统动力学结构突变前后的状态进行表征，极大地丰富了系统动力学结构突变检测理论。其中，标度指数的快速、准确的计算对于系统动力学结构突变的检测尤为重要。目前标度指数的计算方法主要有重标极差分析法、去趋势波动分析、小波分析法、重标方差法等。其中，重标极差分析法是最经典的标度指数分析方法，但在长序列的分析中易受短期相关性和周期的影响，导致其结果精度较差；重标方差法分析法是重标极差分析法的一种改进，具有较强的稳健性，但在进行大数据量的分析中计算时间较长；去趋势波动分析方法可以有效地消除序列中的各阶趋势成分，能够很好地处理具有噪声及趋势的非平稳信号；通过小波分析法来计算标度指数是近些年来的热点，由于小波分析是在序列尺度和时间域上进行，是对序列长相关特性的表征，通过小波变换系数对序列的标度指数进行估计，能够更好地克服长序列的周期、短期相关性等问题，且具有计算速度快、收敛性好、精度高等特点，适合非平稳数据的分析。

实际中由于地球物理气候系统是一个复杂、动态的非线性系统，在采用基于标度指数理论的方法对系统动力学结构突变进行检测时，周期性、短期相关性等非平稳特征会对标度指数计算精度产生一定的影响。其次，对于大数据量的分析，由于算法本身的原因造成计算时间的太长则会给系统动力学结构突变分析带来不便。因此，如何快速、有效、准确地提取序列的标度指数是进行相关时间序列动力学结构突变分析的一个关键性问题。鉴于小波分析法在标度指数计算过程中的快速、准确和有效性，发明一种新的动力学结构突变检测方法——滑动移除小波分析法，该方法是基于数据的移除对于具有相同动力学属性的相关序列标度指数的估算几乎没有影响的这一特征而提出。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种系统动力学结构突变的检测方法，通过滑动移除技术从原序列中选取子序列，采用小波变换系数估计子序列的标度指数，依据子序列标度指数趋势判断系统动力学结构发生突变点或突变区间，并对突变区间进行验证，能够快速、有效和准确地判断系统动力学结构的突变情况。

技术方案：本发明所述的滑动移除小波分析法包括以下步骤：

(1)定义滑动移除窗口长度和滑动步长

假设描述系统状态的时间序列x(t)的数据量为N，定义滑动移除窗口长度为h，滑动步长为l。

(2)定义子序列

以步骤(1)定义的滑动移除窗口长度h从x(t)的第t·l(t＝1，2，…，int(N/l))个数据开始连续移除h个数据，将剩余N-h个数据连接形成子序列，其中int表示取整。

(3)计算子序列的标度指数

以步骤(2)定义的子序列进行Mallat一维分解，计算小波系数d_x(j，k)，其中j＝1，2，…，J；k＝1，2，…，2^-jn；(j为尺度参数，J为小波分解层数，k为位置参数，n＝N-h为子序列长度)；

由小波系数d_x(j，k)计算中间参量η_j，s_j(j＝1，2，…，J)：

s_j＝(n ln²2)/2^j+1

式中n_j＝2^-jn。

计算子序列的标度指数γ的小波估计值H_w(j₁，j₂)，1≤j₁≤j₂≤J：

重复以上(2)-(3)步骤，得到M＝int(N/l)个标度指数γ值。

(4)检测系统动力学结构的突变

绘制步骤(3)得到的标度指数γ值序列随时间变化的曲线，根据曲线的变化趋势初步判断系统的动力学结构的突变点或突变区间，对于突变区间进一步地采用标度指数方差贡献值进行验证。

(5)计算标度指数γ值的方差贡献

将步骤(3)计算得出的标度指数γ值序列定义为z(i)(i＝1，2，…，M)序列，利用如下公式计算每个标度指数γ的方差贡献C(i)值，即：

式中，表示整个标度指数序列z(i)的均值。

(6)验证动力学结构突变检测结果

绘制步骤(5)中C(i)序列随时间变化的曲线，定义方差贡献阈值为3倍的z(i)序列标准差值S，基于方差贡献C(i)值是否超过阈值S来验证系统动力学结构是否发生区间突变。

进一步，前述步骤1中，滑动移除窗口长度h至少包含2个数据点以确保移除后新序列标度指数γ值计算的差异性。

进一步，前述步骤1中，滑动步长l应小于等于滑动移除窗口长度h，以确保不会遗漏数据，取l＝h为宜。

进一步，前述步骤6中，标准差值S：

本发明与现有技术相比，其有益效果是：

(1)本发明将小波标度指数计算法与滑动移除技术相结合，不仅能够快速、准确地计算子序列的标度指数，从而进行系统动力学结构突变分析，而且能够更好地克服大数据量序列标度指数估计中的周期长、误差大等问题，适用于大数据量序列突变分析。

(2)本发明不仅能够检测出系统的动力学结构突变点或突变区间，而且能够对突变区间的正确性进行验证，准确地描述系统突变前后的动力学结构状态。

附图说明

图1a为两种具有不同动力学稳定结构的系统的理想时间序列示意图。

图1b为采用本发明的方法对图1a的两种具有不同动力学稳定结构系统的突变检测结果示意图。

图2a为一种具有动力学稳定结构和随机系统的理想时间序列示意图。

图2b为采用本发明的方法对图2a的具有动力学稳定结构和随机系统的序列突变检测结果示意图。

图3a为一种动力学稳定结构系统内的随机系统的理想时间序列示意图。

图3b为采用本发明的方法对图3a的突变初步诊断示意图。

图3c为采用本发明的方法对图3a的突变检测结果验证示意图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

记观测时间序列集为D＝{x(t)，t＝1，2，…，N}，其中N为序列总长度，滑动移除窗口为h∈N，滑动步长l∈N。用数学公式可以把滑动移除窗口后M个子序列y(t)描述如下：

式中τ＝1，2，…M，

对每个子序列进行Mallat一维分解，计算各子序列的小波系数d_x(j，k)，其中j＝1，2，…，J；k＝1，2，…，2^-jn。(j为尺度参数，J为小波分解层数，k为位置参数，n为子序列长度)

由小波系数计算中间参量η_j，s_j(j＝1，2，…，J)：

s_j＝(n ln²2)/2^j+1 (3)

式中n_j＝2^-jn。

计算各子序列的标度指数γ的小波估计值H_w(j₁，j₂)，1≤j₁≤j₂≤J：

绘制各子序列标度指数γ值随时间变化的曲线，并基于该曲线的变化情况判断系统动力学结构的突变点或突变区间，对于突变区间进一步地采用以下标度指数方差贡献值进行验证。

将计算得出的各子序列标度指数γ值序列定义为z(i)(i＝1，2，…，M)序列，利用如下公式计算每个标度指数γ的方差贡献C(i)值：

式中表示整个标度指数序列z(i)的均值。

绘制C(i)序列随时间变化的曲线，定义方差贡献阈值为3倍的z(i)序列标准差值S，即：

基于子序列标度指数的方差贡献C(i)值是否超过阈值S来验证系统动力学结构是否发生区间突变。

在上述计算中，滑动步长l应小于等于滑动移除窗口长度h，以确保不会遗漏数据，经过多次试验，取l＝h为宜。滑动移除窗口长度h的选择则需要从大到小逐渐缩小范围，直到计算结果稳定为止。

根据本发明的上述实施例，针对一种突变发生在两种具有不同动力学稳定结构的系统之间的情况，构造如下一理想时间序列G(t)：

图1(a)给出了序列G(t)随时间变化的情况。显然，系统在t＝1000处发生了突变，由一种稳定的动力学结构突变为另一种稳定的动力学结构。图1(b)为采用本发明方法突变检测结果，取滑动步长、滑动移除窗口长度为l＝h＝2。从图1(b)可以看到，在t＝1000处，标度指数γ值发生了一次明显的突变，其演化过程和序列G(t)完全一致。

根据本发明的上述实施例，针对一种突变发生在一种具有动力学稳定结构和随机系统之间的情况，待分析时间序列长度为2000个，其中前1000个数据点由混沌系统Logistic虫口模型产生，后1000个数据由符合正态分布的随机数产生。Logistic虫口模型如下：

x_n+1＝ux_n(1-x_n)，x∈[0，1] (8)

式中x为虫口数状态变量，大于0小于1，x_n记为第n代虫口数，x_n+1为第n+1代虫口数。当3.569945672＜u＜4.0，系统进入混沌状态。选取虫口初值x₀＝0.8，控制参数u＝3.8。图2(a)为采用该方法产生的序列随时间变化情况。图2(b)为采用本发明方法突变检测的初步结果，取滑动步长、滑动移除窗口长度为l＝h＝25。显然t＝1000处发生了一次明显的突变，系统由一种稳定的状态突变为一种不稳定的状态，即Logistic映射所产生的数据序列变化幅度相对平稳，而由随机数据生成的数据其变化幅度相对较大。

根据本发明的上述实施例，针对另一种突变发生在动力学稳定结构系统之内情况，待分析时间序列长度为1000个，由混沌系统Logistic虫口模型产生，其中第300-330个数据点由均匀分布的随机数模拟代替。图3(a)为采用本发明方法所产生的数据序列随时间变化情况，取滑动步长、滑动移除窗口长度为l＝h＝30。图3(b)为采用本发明方法突变检测的初步结果。显然，标度指数γ值在区间[301，330]处发生了突变，脱离了稳定的状态。进一步验证该区间突变的可靠性，对标度指数γ序列进行方差贡献分析，图3(c)为给出了方差贡献分析的结果。显然，除了在区间[301，330]内标度指数γ计算的方差贡献超过了三倍方差阈值，在其他区域内方差贡献基本接近于0值，可以判定序列在区间[301，330]发生了突变。

如上，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (4)

1.一种系统动力学结构突变的检测方法，其特征在于，包括以下步骤：

(1)定义滑动移除窗口长度和滑动步长

假设描述系统状态的时间序列x(t)的数据量为N，定义滑动移除窗口长度为h，滑动步长为l。

(2)定义子序列

以步骤(1)定义的滑动移除窗口长度h从x(t)的第t·l(t＝1，2，…，int(N/l))个数据开始连续移除h个数据，将剩余N-h个数据连接形成子序列，其中int表示取整。

(3)计算子序列的标度指数

以步骤(2)定义的子序列进行Mallat一维小波分解，计算其小波系数d_x(j，k)，其中j＝1，2，…，J；k＝1,2，…，2^-jn；(j为尺度参数，J为小波分解层数，k为位置参数，n＝N-h为子序列长度)；

由小波系数d_x(j，k)计算中间参量η_j，s_j(j＝1，2，…，J)：

s_j＝(nln²2)/2^j+1

式中n_j＝2^-jn。

计算子序列的标度指数γ的小波估计值H_w(j₁，j₂)，1≤j₁≤j₂≤J：

(4)检测系统动力学结构的突变

重复以上(2)-(3)步骤，得到M＝int(N/l)个标度指数γ值。绘制标度指数γ值随时间变化的曲线，根据曲线的变化趋势初步判断系统的动力学结构的突变点或突变区间，对于突变区间进一步地采用标度指数方差贡献值进行验证。

(5)计算标度指数γ值的方差贡献

将步骤(3)计算得出的标度指数γ值序列定义为z(i)(i＝1，2，…，M)序列，利用如下公式计算每个标度指数γ的方差贡献C(i)值，即：

式中，表示整个标度指数序列z(i)的均值。

(6)验证动力学结构突变检测结果

绘制步骤(5)中C(i)序列随时间变化的曲线，定义方差贡献阈值为3倍的z(i)序列标准差值S，基于方差贡献C(i)值是否超过阈值S来验证系统动力学结构是否发生区间突变。

2.根据权利要求1所述的系统动力学结构突变的检测方法，其特征在于步骤1中，滑动移除窗口长度h至少包含2个数据点以确保移除后新序列标度指数γ值计算的差异性。

3.根据权利要求1所述的系统动力学结构突变的检测方法，其特征在于步骤1中，滑动步长l应小于等于滑动移除窗口长度h，以确保不会遗漏数据，取l＝h为宜。

4.根据权利要求1所述的系统动力学结构突变的检测方法，其特征在于步骤(6)中，标准差值S：

。