基于时间尺度局部Hurst指数的时空标度分析方法与流程

本发明涉及一种标度分析方法，特别是涉及一种基于时间尺度局部Hurst指数的时空标度分析方法，属于数据分析技术领域。

背景技术：

自然过程或受多因素影响的复杂系统生成的序列过程通常体现为无序、非稳态、随机和非线性，传统的统计方法已经无法表征这种复杂的波动过程。近年来，研究发现复杂的时间结构的震荡行为是无尺度特性的，而自相似性是这种无特征尺度行为的必然性质。Hurst指数作为判断时间序列遵从随机行走还是分数有偏随机游走过程的指标，被广泛用来表征具有不同相关性的标度行为。当0<H<0.5时，波动为反持续波动，当0.5<H<1时，为长相关性波动。

目前，估计Hurst指数的方法主要有重标极差估计法(R/S)、小波变换极大模值法(WTMM)、波动分析法(FA)、频谱分析法、广义Hurst指数估计法(GHE)、去趋势波动分析法(DFA)、多重去趋势波动分析(MF-DFA)等(见：文献1Kantelhardt J W.Fractal and Multifractal Time Series[J].Physics,2008:463-487)，这些方法被成熟地广泛应用于时间序列的标度行为或多重标度行为的估计上，为判别、比较和模拟动力系统提供了识别参数。但是，这些方法得到的Hurst指数仅仅表征的是整体的标度行为，也就是依靠估计的一个或有限的Hurst指数对生成序列的动力系统进行行为特征估计。然而，实际情况下在各种内在或外界因素的影响下，标度行为是随着时间不断发生变化的，而这些方法无法估计具体某一时刻或具体某一事件影响下的局部Hurst指数。

基于以上原因，局部Hurst指数的估计(见：文献2Ihlen E A.Introduction to Multifractal Detrended Fluctuation Analysis in Matlab[J].Frontiers in Physiology,2012,3(141):141)为我们提供了对复杂动力系统进行局部分析和深入理解的必要认识，弥补了传统估计方法结果的不足。但据我们所知，局部Hurst指数是由局部波动所估计的，单一的局部Hurst指数对标度行为特征没有明显意义上的体现。

因此，需要建立一种依托局部Hurst指数反映整体行为的分析技术以解决完善上述方法所存在的问题。

技术实现要素：

本发明的主要目的在于，克服现有技术中的不足，提供一种基于时间尺度局部Hurst指数的时空标度分析方法，从统计角度以时间尺度局部Hurst指数的分布构造分形结构，通过统计特征值量化不同空间和时间上标度结构发生的转变，从而反映整体行为。

为了达到上述目的，本发明所采用的技术方案是：

一种基于时间尺度局部Hurst指数的时空标度分析方法，包括以下步骤：

1)选择研究对象，获取目标数据，建立相应的时间序列；

2)检验时间序列中的目标时间序列是否含周期成分、是否呈随机行走形态；

若是，则去除目标时间序列的周期成分，生成对应的波动增量时间序列；

若否，即目标时间序列不含周期成分、且呈噪音形态，则该目标时间序列为波动增量时间序列，直接执行下一步；

3)计算波动增量时间序列的局部Hurst指数，构成H_t序列；

4)对不同空间环境背景的H_t序列进行频数统计，计算基于不同空间下的概率密度分布；

5)进行空间上的标度行为分析；

5-1)对基于不同空间下的概率密度分布进行高斯拟合，确定原分布和拟合分布的空间尺度下的统计特征值；

5-2)通过空间尺度下的统计特征值，考察并量化标度行为随空间因素的变化；

6)进行时间上的标度行为分析；

6-1)按照特征时间对H_t序列进行区间划分，确定每个时间区间段的H_t序列的基于不同时间下的概率密度分布，并确定时间尺度下的统计特征值；

6-2)通过时间尺度下的统计特征值，考察并量化标度行为随时间因素的变化。

本发明进一步设置为：所述步骤1)中的时间序列，是将待考察的现象或已测得的数据整理为连续的时间序列u，u＝u₁,u₂,…,u_i,…,u_N,i∈[1,N]；其中，N为自然数。

本发明进一步设置为：所述步骤2)中的去除目标时间序列的周期成分，生成对应的波动增量时间序列，具体为，

令u_i为时间序列u中的目标时间序列，v_i为去除u_i的周期成分后的时间序列，则其中，<>为均值运算；

令x_i为v_i的波动增量时间序列，则x_i＝v_i+1-v_i。

本发明进一步设置为：所述步骤3)中的计算波动增量时间序列的局部Hurst指数，具体为，

3-1)将波动增量时间序列x_i转变为随机行走结构序列Y(k)，

其中，<x>为波动增量时间序列均值；

3-2)确定计算局部Hurst指数的区间长度s，保证s尺度大小的区间内囊括我们要分析的特征波动，使得计算的局部Hurst指数具有表征意义，同时保证为后面的局部Hurst指数的统计分析提供足够的支持；

对Y(k)序列进行等距重叠划分，要求根据确定的区间长度s，保证重叠区间的区间中值是连续的；

3-3)每个区间内，用p阶多项式进行局部拟合，v＝1,…,N_v，N_v＝N-s，p取1或大于1的阶数；

3-4)得出s尺度下每个区间的局部波动均方根RMS，

其中，为第v局部区间的p阶多项式拟合，p取1或大于1的阶数；

3-5)以q取0时、对数坐标下定义的直线F_q(s)～s^H(0)定义一条回归线，斜率H(0)等于0阶Hurst指数；

根据确定的区间尺度s值，确定回归线上相应的波动函数值F₀(s)；

在区间尺度为s值情况下，残差波动均方根resRMS(v)为，

resRMS(v)＝logF₀(s)-log(RMS(v))

即为波动函数值F₀(s)与局部波动均方根RMS{s}(v)在对数坐标下的差；

局部Hurst指数H_t由下式求得，

其中，maxL为随机行走结构序列Y(k)的长度。

本发明进一步设置为：所述步骤3-5)中的0阶Hurst指数H(0)，F_q(s)～s^H(0)，具体为，

3-5-1)确定区间尺度s，保证s的大小不超过总序列长度的十分之一；

3-5-2)将波动增量时间序列x_i转变为随机行走结构序列Y(k)，

根据确定的区间尺度s对随机行走结构序列Y(k)进行等距非重叠划分；

3-5-3)对每个区间进行去除趋势计算，计算随机行走结构序列Y(k)每个区间内相对局部趋势偏差的均方根

其中，为第v局部区间的p阶多项式拟合，p取1或大于1的阶数；

3-5-4)改变区间尺度s的大小，重复步骤3-5-1)至步骤3-5-3)，由求得0阶Hurst指数H(0)，

其中，N_s＝N/s。

本发明进一步设置为：所述步骤5)中的高斯拟合，是用高斯分布拟合H_t序列的分布；

其高斯分布公式为其中，σ为标准差，h为拟合变量，H₀为均值。

本发明进一步设置为：所述步骤6)中的确定时间尺度下的统计特征值，具体是通过主导Hurst指数，以概率密度分布的最大值对应的Hurst指数确定的。

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

本发明的时空标度分析方法，以时间尺度局部Hurst指数为基础构造分形结构，可以对局部具体时刻的标度行为或具体某外界因素的影响进行识别，避免了传统方法得到的Hurst指数所估计的标度指数只能表征整体空间标度行为的缺点；并利用统计分布的方式构造分形结构，简单直接地表示了时间序列的奇异性，统计特征值被用来表征整体的标度行为，提供了直接可比较的方法；同时，结合时间尺度局部Hurst指数和统计分布，拓展了标度分析的内容，建立了空间行为变化分析和时间变化分析，丰富了识别系统的途径，为多领域的相关性分析提供了系统的分析方法，具有重要的理论和工程意义。

上述内容仅是本发明技术方案的概述，为了更清楚的了解本发明的技术手段，下面结合附图对本发明作进一步的描述。

附图说明

图1为本发明基于时间尺度局部Hurst指数的时空标度分析方法的流程图；

图2为三个不同测量点d1、d2、d3的地下水波动数据及计算得出的H_t序列；

图2(a)为测量点d1的地下水波动数据及计算得出的H_t序列；

图2(b)为测量点d2的地下水波动数据及计算得出的H_t序列；

图2(c)为测量点d3的地下水波动数据及计算得出的H_t序列；

图3为不同测量点的局部指数H_t序列的概率密度分布与拟合曲线；

图4为不同测量点的局部指数H_t序列的概率密度分布逐年变化；

图5为不同测量点的主导Hurst指数的逐年变化及趋势线。

具体实施方式

下面结合说明书附图，对本发明作进一步的说明。

本发明提供一种基于时间尺度局部Hurst指数的时空标度分析方法，如图1所示，包括以下步骤：

1)选择研究对象，获取目标数据，建立相应的时间序列。

其中的时间序列，是将待考察的现象或已测得的数据整理为连续的时间序列u，u＝u₁,u₂,…,u_i,…,u_N,i∈[1,N]；其中，N为自然数。

2)检验时间序列中的目标时间序列是否含周期成分、是否呈随机行走形态。

若是，则去除目标时间序列的周期成分，生成对应的波动增量时间序列；具体为，

令u_i为时间序列u中的目标时间序列，v_i为去除u_i的周期成分后的时间序列，则其中，<>为均值运算；

令x_i为v_i的波动增量时间序列，则x_i＝v_i+1-v_i。

若否，即目标时间序列不含周期成分、且呈噪音形态，则该目标时间序列为波动增量时间序列，直接执行下一步。

3)计算波动增量时间序列的局部Hurst指数，以去除趋势波动分析(DFA)和多重去除趋势波动分析(MF-DFA)算法为基础，构成H_t序列。

3-1)将波动增量时间序列x_i转变为随机行走结构序列Y(k)，

其中，<x>为波动增量时间序列均值；

3-2)确定计算局部Hurst指数的区间长度s，保证s尺度大小的区间内囊括我们要分析的特征波动，使得计算的局部Hurst指数具有表征意义，同时保证为后面的局部Hurst指数的统计分析提供足够的支持；

对Y(k)序列进行等距重叠划分，要求根据确定的区间长度s，保证重叠区间的区间中值是连续的；

3-3)每个区间内，用p阶多项式进行局部拟合，v＝1,…,N_v，N_v＝N-s，p取1(DFA1)或大于1的阶数(DFA2，DFA3，……)；

3-4)对每个区间进行去除趋势计算，得出s尺度下每个区间的局部波动均方根RMS，

其中，为第v局部区间的p阶多项式拟合，p取1或大于1的阶数；

3-5)以q取0时、对数坐标下定义的直线F_q(s)～s^H(0)定义一条回归线，斜率H(0)等于0阶Hurst指数；

根据确定的区间尺度s值，确定回归线上相应的波动函数值F₀(s)；

在区间尺度为s值情况下，残差波动均方根resRMS(v)为，

resRMS(v)＝logF₀(s)-log(RMS(v))

即为波动函数值F₀(s)与局部波动均方根RMS{s}(v)在对数坐标下的差；

局部Hurst指数H_t由下式求得，

其中，maxL为随机行走结构序列Y(k)的长度。

步骤3-5)中的0阶Hurst指数H(0)，F_q(s)～s^H(0)，具体为，

3-5-1)确定区间尺度s，保证s的大小不超过总序列长度的十分之一；

3-5-2)将波动增量时间序列x_i转变为随机行走结构序列Y(k)，

根据确定的区间尺度s对随机行走结构序列Y(k)进行等距非重叠划分；

3-5-3)对每个区间进行去除趋势计算，计算随机行走结构序列Y(k)每个区间内相对局部趋势偏差的均方根

其中，为第v局部区间的p阶多项式拟合，p取1(DFA1)或大于1的阶数(DFA2，DFA3，……)；

3-5-4)改变区间尺度s的大小，重复步骤3-5-1)至步骤3-5-3)，由求得0阶Hurst指数H(0)，

其中，N_s＝N/s。

4)对不同空间环境背景的H_t序列进行频数统计，计算基于不同空间下的概率密度分布。

5)进行空间上的标度行为分析。

5-1)对基于不同空间下的概率密度分布进行高斯拟合，确定原分布和拟合分布的空间尺度下的统计特征值；

其中的高斯拟合，是用高斯分布拟合H_t序列的分布；

其高斯分布公式为其中，σ为标准差，h为拟合变量，H₀为局部指数序列的均值。

5-2)通过空间尺度下的统计特征值，考察并量化标度行为随空间因素的变化。

6)进行时间上的标度行为分析。

6-1)按照特征时间对H_t序列进行区间划分，确定每个时间区间段的H_t序列的基于不同时间下的概率密度分布，并确定时间尺度下的统计特征值；

其中的确定时间尺度下的统计特征值，具体是通过主导Hurst指数以概率密度分布的最大值对应的Hurst指数确定的。

6-2)通过时间尺度下的统计特征值，考察并量化标度行为随时间因素的变化。

实施例：

实施例数据采用来自美国国家地质调查局(USGS)按日测量的地下水波动数据，如图2所示，该数据位于密西西比河下游流域的哈奇河流域，收集了三个不同数据测量点水井连续五年测得的地下水波动数据。其中，图2(a)为d1数据测量点位于哈奇河下游近河岸测量点的地下水波动数据，图2(b)为d2数据测量点位于哈奇河下游远离河岸测量点的地下水波动数据，图2(c)为d3数据测量点位于哈奇河上游测量点的地下水波动数据，三个测量点所在的水文环境具有各不相同的特征。

根据示例数据，局部区间尺度s大小取为30天，计算得出相应的H_t序列，如图2所示。

对三个测量点的H_t序列的概率分布进行高斯分布拟合，如图3所示。考虑统计参数主导Hurst指数H评估整体标度行为，均值H₀，标准差σ为高斯分布拟合参量，得到三个测量点相应的水文环境下的局部Hurst指数的统计特征值，见表1。

表1

由表1可得，由下游远离河岸测量点d2到下游近河岸测量点d1，再到上游测量点d3，三个测量点的H_t序列分布向右移动；这意味着波动相关性在逐渐加强，同时分布的范围变大，意味着标度行为的复杂性在变大，意味着上游流域复杂的水文环境。进一步地，d1测量点的主导Hurst指数说明其波动呈现的是长相关稳态波动；d2测量点主导Hurst指数说明其波动呈反持续波动；d3测量点主导Hurst指数说明其波动呈长相关性非稳态波动。而Hurst指数则说明局部Hurst指数分布存在正偏斜。

对每个波动时间序列，按年进行区域划分，即划分为5个区间，对每一年的局部Hurst指数进行数理统计，如图4所示，则每个测量点逐年的标度统计结构变化，可以对比整体序列分布图3得出，图4中的背景灰色分布即为得出的五年整体序列分布。

从图4可以发现，d1测量点的标度行为变化是逐年波动持续性增强；d2测量点的整体行为趋于稳定，但分布范围逐渐变大，说明它的标度行为复杂程度在变大；d3测量点的变化呈现无规律性。

将逐年分布的主导Hurst指数及其趋势线按年画出，如图5所示，可以清楚地表示出按时间的变化趋势。

本发明的创新点在于，从统计角度以时间尺度局部Hurst指数的分布构造分形结构，通过统计特征值量化不同空间和时间上标度结构发生的转变，从而反映整体行为。

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。