基于图像微结构频数分析的时间序列复杂度测算方法与流程

本发明涉及一种时间序列复杂度测算方法。特别是涉及一种基于图像微结构频数分析的时间序列复杂度测算方法。

背景技术：

复杂性科学兴起于20世纪80年代，是系统科学发展的新阶段，也是当代科学发展的前沿领域之一。复杂性科学的发展，不仅引发了自然科学界的变革，而且也日益渗透到哲学、人文社会科学领域。随着大数据的流行，基于数据驱动的方法成为研究复杂系统的重要手段，而时间序列分析成为探究和解读复杂系统内在运行机理的强有力工具，相关成果已广泛应用于物理学、经济学、气象学、语言学和信息科学等众多领域。复杂性测度是(非线性)时间序列最重要的动力学指标之一，众多科研人员致力于该问题的研究并提出了一系列有益的解决方法，目前应用最广的一类方法是基于熵概念的方法，其中经典的方法包括K-S熵、近似熵、排列熵和功率谱熵等。这些方法各具特点，例如，K-S熵描述了动力学系统轨道随时间演化信息的产生率，用来度量系统运动的混乱或无序程度。近似熵是从不同维数条件下新模式产生的概率大小测算时间序列的复杂度。排列熵是依据各种排列模式出现的相对频率测算时间序列的随机程度。功率谱熵是通过分析信号的频谱结构测算时间序列的不确定性及复杂程度。基于熵概念方法的优势在于其深刻的物理背景及方法设计的多样性，研究人员可以依据不同的信息侧重，构造丰富多样的熵方法。研究发现，现有方法主要局限在时域和频域两个角度设计熵方法，因此复杂性分析存在一定的片面性。

技术实现要素：

本发明所要解决的技术问题是，提供一种将复杂的一维时间序列映射为二维图像进行处理的基于图像微结构频数分析的时间序列复杂度测算方法。

本发明所采用的技术方案是：一种基于图像微结构频数分析的时间序列复杂度测算方法，包括如下步骤：

1)构造信号的递归矩阵；

2)绘制灰度图，根据第i时刻和第j时刻的递归状态的递归矩阵绘制灰度图像；

3)对灰度图像I(x,y)，利用不同尺度的高斯核函数进行滤波，得到不同高斯灰度图像，形成高斯金字塔，对不同尺度的高斯核函数与灰度图像I(x,y)进行卷积运算得到图像的尺度空间C(x,y,R)：

其中x,y为灰度图像坐标，G(x,y,R)代表高斯核函数，代表卷积操作，R为尺度大小；

4)在尺度空间中初步确定特征点位置和尺度；

5)利用尺度空间函数的泰勒函数二次展开式进行最小二乘拟合，利用拟合曲面的极值来去除不稳定特征点，

其中H为Hessian矩阵，计算公式为：

6)对剩余特征点进行聚类，使各个剩余特征点与所属类代表点的相似度之和最大，用公式可以表示为：

其中i是第i个样本点，b_i指第i个样本点的类代表点，w(i,b_i)是指第i个样本点与样本点的类代表点之间的欧几里德距离的负值，当一个类代表点k同时又是其他类别中的点时，这种聚类结果是错误的并需要被剔除掉，故此时δ_k(b)＝-∞，其余情况δ_k(b)＝0；

7)为了得到有效的聚类中心，通过计算不断改变影响度M(i,k)和有效度N(i,k)的值，其中M(i,k)反映的是k点是否适合作为i点的聚类中心，而N(i,k)则反映i点是否选择k点作为i点聚类中心，计算公式如下：

8)将聚类结果进行信息度量，过程如下：设得到特征点数量为Num，聚类结果为s类，其中第i类包含的点数为n，i＝1,....,s，则第i类出现的概率为：

P(i)＝n/Num (7)

定义时间序列的微结构递归熵为：

9)使用近似熵和排列熵计算不同信号的复杂度，并与微结构递归熵的结果进行对比分析。

步骤1)所述的构造信号的递归矩阵是：

给定时间序列x₁.x₂,...,x_L，利用延迟坐标法在嵌入空间上重构嵌入向量：

X_k＝(x_k.x_k+τ,...,x_k+(m-1)τ) (9)

式中k＝1,2,...,L-(m-1)τ，τ为延迟时间，m为嵌入维数；在相空间中，第i时刻和第j时刻的递归状态的递归矩阵表示：

Y_i,j＝ε-||X_i-X_j||，i,j＝1,...,N (10)

式中，N是状态向量X_i的个数，ε为判别距离，||·||表示L₂范数，作递归图时，对递归矩阵二值化，即

步骤4)包括对相邻尺度的两个高斯灰度图像相减得到高斯差分尺度空间，对尺度空间每个点与相邻位置和相邻尺度的所有点进行比较，得到局部极值和极值对应尺度，

D(x,y,R)＝(G(xykR)-G(xyR))*I(xy) (12)

其中k常数，D(x,y,R)为尺度空间函数。

本发明的基于图像微结构频数分析的时间序列复杂度测算方法，将复杂的一维时间序列映射为二维图像，然后借助强大的图像分析工具统计图像局部微结构的频数，从而设计新的时间序列复杂性测算的熵方法。本发明的方法更进一步丰富和拓展了熵概念，并为新熵方法的设计提供有益的参考。本发明主要优点及特色体现在如下几个方面：

1)本发明从时间序列的图像微结构这一全新的视角设计新的时间序列复杂性测算方法，因此可以挖掘出蕴藏在复杂时间序列中时域和频域方法无法获取的新信息。

2)本发明提出的熵方法将复杂的一维时间序列映射为二维图像，然后借助强大的图像分析工具统计图像局部微结构的频数，从而测算时间序列的复杂性。这种时间序列图模式分析思路更进一步丰富和拓展了熵概念，并为新的熵方法的设计提供有益的参考。

3)本发明提出的基于图像微结构频数分析的时间序列复杂度测算方法是一种全新的普适性算法，可广泛应用于生物学、经济学、信息科学等多个领域。

附图说明

图1是本发明基于图像微结构频数分析的时间序列复杂度测算方法的流程图；

图2a是正弦信号的递归图；

图2b是p＝0.2的混合信号的递归图；

图2c是p＝0.5的混合信号的递归图；

图2d是p＝0.7的混合信号的递归图；

图2e是高斯白噪声信号的递归图；

图2f是Lorenz信号的递归图；

图3a是正弦信号的灰度图；

图3b是p＝0.2的混合信号的灰度图；

图3c是p＝0.5的混合信号的灰度图；

图3d是p＝0.7的混合信号的灰度图；

图3e是高斯白噪声信号的灰度图；

图3f是Lorenz信号的灰度图；

图4a是正弦信号的递特征点聚类结果图；

图4b是p＝0.2的混合信号的特征点聚类结果图；

图4c是p＝0.5的混合信号的特征点聚类结果图；

图4d是p＝0.7的混合信号的特征点聚类结果图；

图4e是高斯白噪声信号的特征点聚类结果图；

图4f是Lorenz信号的特征点聚类结果图。

具体实施方式

下面结合实施例和附图对本发明的基于图像微结构频数分析的时间序列复杂度测算方法做出详细说明。

本发明的基于图像微结构频数分析的时间序列复杂度测算方法，将从时间序列的图像微结构这一全新的视角设计新的时间序列复杂性测算方法，该方法可挖掘出蕴藏在复杂时间序列中时域和频域方法无法获取的新信息。

如图1所示，本发明的基于图像微结构频数分析的时间序列复杂度测算方法，包括如下步骤：

1)构造信号的递归矩阵；具体是：

给定时间序列x₁.x₂,...,x_L，利用延迟坐标法在嵌入空间上重构嵌入向量：

X_k＝(x_k.x_k+τ,...,x_k+(m-1)τ) (1)

式中k＝1,2,...,L-(m-1)τ，τ为延迟时间，m为嵌入维数；在相空间中，第i时刻和第j时刻的递归状态的递归矩阵表示：

Y_i,j＝ε-||X_i-X_j||，i,j＝1,...,N (2)

式中，N是状态向量X_i的个数，ε为判别距离，||·||表示L₂范数，作递归图时，对递归矩阵二值化，即

2)绘制灰度图，根据第i时刻和第j时刻的递归状态的递归矩阵绘制灰度图像；

3)对灰度图像I(x,y)，利用不同尺度的高斯核函数进行滤波，得到不同高斯灰度图像，形成高斯金字塔，对不同尺度的高斯核函数与灰度图像I(x,y)进行卷积运算得到图像的尺度空间C(x,y,R)：

其中x,y为灰度图像坐标，G(x,y,R)代表高斯核函数，代表卷积操作，R为尺度大小；

4)在尺度空间中初步确定特征点位置和尺度；包括对相邻尺度的两个高斯灰度图像相减得到高斯差分尺度空间(DoG)，对尺度空间每个点与相邻位置和相邻尺度的所有点进行比较，得到局部极值和极值对应尺度，

D(x,y,R)＝(G(x y kR)-G(x y R))*I(x y) (6)

其中k常数，D(x,y,R)为尺度空间函数。

5)由于高斯差分尺度空间会产生较强的边缘效应，会导致极值点对比度偏低，为了保证结果的可靠性，要对高斯差分尺度空间中检测到的局部极值点进一步筛选才能选为最终的特征点。利用尺度空间函数的泰勒函数二次展开式进行最小二乘拟合，利用拟合曲面的极值来去除不稳定特征点，

其中H为Hessian矩阵，计算公式为：

6)对剩余特征点进行聚类，使各个剩余特征点与所属类代表点的相似度之和最大，相似度值储存在矩阵W中，用公式可以表示为：

其中i是第i个样本点，b_i指第i个样本点的类代表点，w(i,b_i)是指第i个样本点与样本点的类代表点之间的欧几里德距离的负值，当一个类代表点k同时又是其他类别中的点时，这种聚类结果是错误的并需要被剔除掉，故此时δ_k(b)＝-∞，其余情况δ_k(b)＝0；

7)为了得到有效的聚类中心，通过计算不断改变影响度M(i,k)和有效度N(i,k)的值，其中M(i,k)反映的是k点是否适合作为i点的聚类中心，而N(i,k)则反映i点是否选择k点作为i点聚类中心，计算公式如下：

8)将聚类结果进行信息度量，过程如下：设得到特征点数量为Num，聚类结果为s类，其中第i类包含的点数为n，i＝1,....,s，则第i类出现的概率为

P(i)＝n/Num (11)

定义时间序列的微结构递归熵为：

9)使用近似熵和排列熵计算不同信号的复杂度，并与微结构递归熵的结果进行对比分析。

为了考察微结构递归熵对不同随机信号的识别能力，选取下面序列长度为800点的不同信号进行复杂性计算：

(1)正弦信号y₁＝sin(x)，采样间隔为π/32。

(2)正弦信号与白噪声信号的混合序列，y＝y₁+py₂，其中y₁为正弦信号，y₂为白噪声序列，p为随机成分混入比例，分别取p＝0.2,0.5,0.7。

(3)由matlab产生的实高斯白噪声序列。

(4)Lorenz信号，采用的Lorenz映射方程为：

固定δ＝10，b＝8/3，初始值为(2,2,20)，每隔0.01s取一个点。舍弃前1000个点后取800个点生成序列。

复杂性计算如下：

1)构造不同信号的递归矩阵。

给定时间序列x₁.x₂,...,x_L，利用延迟坐标法在嵌入空间上重构嵌入向量：

X_k＝(x_k.x_k+τ,...,x_k+(m-1)τ) (13)

式中k＝1,2,...,L-(m-1)τ，τ为延迟时间，m为嵌入维数。通常在相空间中，第i时刻和第j时刻的递归状态使用递归矩阵表示：

Y_i,j＝ε-||X_i-X_j||，i,j＝1,...,N (14)

式中，N是状态向量X_i的个数，ε为判别距离，实验中ε＝0.25*std(x)，||·||表示L₂范数。得到不同信号的递归图(如图2a～图2f所示)。

2)对递归矩阵二值化，即

根据灰度矩阵Y_i,j¹做出不同信号的灰度图，如图3a～图3f所示。

3)通过不同的高斯卷积核与灰度图像进行卷积得到多尺度图像，

其中G(x,y,R)代表可变尺度高斯函数，C(x,y,R)代表图像的尺度空间。代表卷积运算，R为尺度因子，I(x,y)代表输入的灰度图像。

4)在尺度空间中检测可能的极值点。对相邻的比例因子为k的尺度高斯函数差分，并与输入的灰度图像进行卷积运算获得差分高斯尺度空间的极值。

D(x,y,R)＝(G(x y kR)-G(x yR))*I(x y) (18)

5)为了得到特征点，每个样本点与其周边的26个点进行比较，若该点是这些点中的最大值或者最小值，那么该点就是可能的极值点。同时为了保证结果可靠性，去除对比度较低的、边界的极值点，及对噪声敏感的极值点。

其中H为Hessian矩阵，计算公式为：

6)对特征点进行聚类。计算所有特征点之间的相似度值，构成相似度矩阵W；

其中i是第i个样本点，b_i指第i个样本点的类代表点，w(i,b_i)是指第i个样本点与其类代表点之间的负的欧几里德距离，当某个类代表点k同时又是其他类别中的点时，这种聚类结果错误并应被排除，故此时δ_k(b)＝-∞，其余情况δ_k(b)＝0

7)在迭代过程中，计算每一次迭代的影响度M(i,k)和有效度N(i,k)，其中M(i,k)说明的是k点是否能作为i点的聚类中心，而N(i,k)则反映i点是否选择k点作为其聚类中心。

最后得到不同信号的特征点聚类结果图，如图4a～图4f所示。

8)根据聚类结果计算不同信号的微结构递归熵H(i)：

其中P(i)表示不同类别出现的概率。

9)使用近似熵和排列熵计算不同信号的复杂度。

表1给出了不同信号的微结构递归熵、近似熵和排列熵的计算结果。测试结果表明，本发明可有效区分不同随机程度的信号，描述了时间序列的演化结构特征。

表1