非线性时间序列的相空间重构技术研究_秦奕青

第20卷第11期 系

统 仿 真 学 报 V ol. 20 No. 11

2008年6月 Journal of System Simulation Jun.,2008

非线性时间序列的相空间重构技术研究

秦奕青

1,3

，蔡卫东

2,3

，杨炳儒

3

(1.北京信息科技大学 计算机学院, 北京 100192; 2. 济南大学 信息科学与工程学院, 济南 250022;

3. 北京科技大学 信息工程学院, 北京 100083)

摘 要：分析了混沌时间序列相空间重构中常用的C-C 方法所存在的四点不足，提出了改进的C-C-2方法。该方法改进了时间序列关联积分的计算方法和参数，利用混沌序列周期N 的概念，提出了通过寻找S cor (t ) 的第一个属于混沌序列周期N 的局部极小峰值，来确定最优延迟时间窗口的判断方式；并只寻找平均∆S 2(t ) 的第一个极小值来确定最优时间延迟，所得结果更合适、稳定，而且将原算法的抗噪能力由30%提高到80%。

关键词：相空间重构；关联积分；延迟时间窗口；非线性时间序列 中图分类号:TP18; O415.5 文献标识码:A 文章编号：1004-731X (2008) 11-2969-05

Research on Phase Space Reconstruction of Nonlinear Time Series

QIN Yi-qing, CAI Wei-dong, YANG Bing-ru

(1. Computer School, Beijing Information Science & Technology University, Beijing 100192, China;

2. School of Information Science and Engineering, University of Ji’nan, Ji’nan 250022, China;

3. School of Information Engineering, University of Science and Technology Beijing, Beijing 100083, China)

1,3

2,3

3

Abstract: A new method called C-C-2 is presented based on the analysis for the C-C method which is usually used to reconstruct the phase space of chaotic time series. First, the C-C-2 method improves the correlation integral algorithm of time series on the computing way and the parameters. Moreover, the method proposes a new way to determine the optimal delay time window by finding the first local minimum peak value of Scor (t) belonging to the period N of the chaotic time series based on the theory of the chaotic system period N. Finally, the method estimates the optimal delay time only by finding the first local minimum of the average ∆S 2(t). The experimental results shows that the C-C-2 method are more stable and more appropriate and also improves the robustness of C-C method from 30% to 80%.

Key words: phase space reconstruction; correlation integral; delay time window; nonlinear time series

引 言混沌时间序列分析与预测的基础是Takens 、Packard 等提出的状态空间的重构理论[1,2]，即把具有混沌特性的时间序列重建为一种低阶非线性动力学系统。通过相空间重构，可以找出混沌吸引子在隐藏区的演化规律，使现有的数据纳入某种可描述的框架之下，从而为时间序列的研究提供了一种崭新的方法和思路。相空间重构是非线性时间序列分析的重要步骤，重构的质量直接影响到模型的建立和预测。

相空间重构的具体方法是用原始系统中某个变量的延迟坐标来重构相空间，Takens [1]从数学上为其奠定了可靠的

基础。他的基本观点是：相空间重构法虽然是用一个变量在不同时刻的值构成相空间，但动力系统的一个变量的变化自即此变量随因此，重构的相空间的轨迹也反映系统状态的演化规律。

相空间重构的原理如下： 收稿日期：2007-03-02 修回日期：2007-09-15

基金项目：国家自然科学基金 (60675030)；山东省教育厅科技计划项目(J06G01)；济南大学科研基金项目(Y0614)。

作者简介：秦奕青(1969-), 女, 北京人, 博士生, 副教授, 研究方向为数据挖掘、建模与仿真；蔡卫东(1968-), 男, 山东济南人, 博士生, 副教授, 研究方向为混沌理论及应用、数据挖掘；杨炳儒(1943-), 男, 天津人, 教授, 博导, 研究方向为人工智能、数据挖掘和柔性建模。

设单变量的时间序列为{x (t i ), i =1,2,..., N }，其采样时间间隔△t 。通过时间延迟构成m 维向量：X (t i ) =(x (t i ), x (t i +τ),..., x (t i +(m −1) τ)), i =1,2,..., M ，其中m 为嵌入维

数，τ为时间延迟，X (t i ) 为m 维相空间中的相点，M 为相点个数，且M =N −(m −1) τ，集合{X (t i ), i =1,2,..., M }描述了系统在相空间中演化轨迹，这时可以在重构的m 维空间中研究系统的混沌行为。研究表明，只要m , τ选择合适，重构的相空间与原系统具有相同的拓扑性质。

在相空间重构中，嵌入维数m 和时间延迟τ的选取具有十分重要的意义，同时这种选取也是很困难的[3]。关于嵌入现在主要有以下两种观点[4-6]： 维数m 和时间延迟τ的选取，

第一种观点认为两者互不相关，即m 和τ的选取是独立进行的。一般情况下，持这种观点的学者，对于最佳嵌入维数m 的选取采用：首先计算出系统的分形维数d ，m 取大于或等于2d +1的整数，目前主要有以下几种方法：G-P 算法、特征值分解方法、替代邻点法等。最佳延迟时间τ的选取比

较复杂，方法也比较多，包括自相关函数法、互信息法等。

以上各种方法，或者是在预先假定一个合适的τ值的情况下求m ；或者是在假定一个合适的m 值的情况下求τ。当我们研究现实经济数据的时候，我们无法事先预知经济系统因此上的有关非线性特征，无从对m 或τ做出合适的假定，述方法都不适用。

第二种观点则认为m 和τ是相互关联的。因为现实中的

• 2969 •

2008年6月 系 统 仿 真 学 报 Jun., 2008

时间序列都是有限长且不可避免地受到各种噪声的影响。大量实验表明，m 和τ的关系与重构相空间的时间窗τW 密切相关(τW =(m −1) τ) ，对于特定的时间序列，其τW 相对固定，

令N →∞有

S 1(m , r , t ) =

1t

[C s (m , r , t ) −C s m (1,r , t )] (6) ∑t s =1

m 和τ的不恰当配对将直接影响重构后的相空间结构与原空间的等价关系，因此相应地产生了m 和τ的联合算法，如

如果时间序列x ={x i }独立同分布，那么对固定的m , t ，当N →∞时，对于所有的r ，均有S 1(m , r , t ) 恒等于零。但实际时间序列是有限长且元素间存在相关性，实际得到的结果一般不等于零。S 1(m , r , t ) ∼t 反映了时间序列的自相关特性，仿照求时间延迟的自相关法原理，最优时间延迟τd 可取

C-C 法、时间窗口法嵌入维、时间延迟自动算法等。

从操作的难易程度、计算量的大小方面来看，尤其是对小数据组而言，1999年由Kim H. S.等[7]提出应用关联积分同时估计出时间延迟τ和延迟时间窗口τW 的C-C 方法更为有效。C-C 方法虽然没有坚实的理论基础，但是在实际应用中有很好的效果，并且该方法相对简单，易于在计算机上实现，能够同时得到τd 和τW ，众多的实验证明该方法还有较好的抗噪声能力[8]。

本文针对C-C 方法在相空间重构中的不足，提出了一种基于C-C 方法的改进算法，来确定最优时间延迟τd 与最优延迟时间窗口τW 。该算法对关联积分的计算、最优延迟时间窗口、最优时间延迟的判断规则都进行了改进，使得最优延迟时间窗口τW 选取更可靠、准确，最优时间延迟τd 的选择更准确，而且还大大提高了原算法的鲁棒性。

S 1(m , r , t ) ∼t 的第一个零点。或者取S 1(m , r , t ) ∼t 对所有半径r 相互差别最小的时间点，此时表示重构相空间中的点最接近均匀分布，重构吸引子轨道在相空间完全展开。

选择最大和最小的两个半径r ，定义差量

∆S 1(m , t ) S 1(m , r i , t )}−min{1(m , r i , t )} (7)

∆S 1(m , t ) 度量了S 1(m , r , t ) ∼t 对所有半径r 的最大偏差。综上所述，最优时间延迟τd 可取S 1(m , r , t ) ∼t 的第一个零点或∆S 1(m , t ) ∼t 的第一个局部极小点。

根据BDS 统计结论可以得到和N , m , r 的合理估计，一

m =, r =k ×σ2, 般情况下取N =3000, 2,3,4,5

k =1,2,3,4, std(σ=x ) , (σ为时间序列x ={x i }的标准差), t =1,2,...,200，计算式(8)、式(9)

1(t ) =

1 C-C方法分析

1.1 C-C方法原理

考虑混沌时间序列x ={x i , i =1,2,..., N }，其采样时间间隔△t 。以m 为嵌入维数，τ为时间延迟，重构相空间X ={X i }，X i 为m 维相空间中的相点

T

144

∑∑S 1(m , r i , t ) (8) 16m =1i =1

∆1(t ) =

14

∑∆S 1(m , t ) (9) 4m =1

寻找1(t ) 的第一个零点或∆1(t ) 的第一个局部极小点即

X i =(x i , x i +τ,..., x i +(m −1) τ) , i =1,2,..., M (1) 为最优时间延迟τd 。另外，由于统计量式(6)采用分块平均的策

则嵌入时间序列的关联积分为[9] 当t =kT 时(k 为大于零的整数) ，略，对于周期为T 的时间序列，

2

C (m , N , r , t ) =∑θ(r −d ij ) , r >0 (2) S 1(t ) 与∆S 1(t ) 均为零。综合考虑S 1(t ) 与∆1(t ) ，定义指标

M (M −1) 1≤i

S 1_cor (t ) =∆S 1(t ) +S 1(t ) (10)

其中m 为嵌入维，N 是时间序列的数据个数，r 为计算中所

寻找S 1_cor (t ) 的全局最小点即可获得最优延迟时间窗口

取的搜索半径，t 为时间延迟，M =N −(m −1) t 表示m 维

τW ，即平均轨道周期的最优估计。

相空间中嵌入点数目，θ为Heaviside 函数：θ(x ) =0, if x

1.2 C-C方法仿真试验

仿真实验采用MATLAB 平台，考察了Lorenz 混沌系统

关联积分是个累积分布函数，表示相空间中任意两点之间距离小于r 的概率。

∞-范数表示。定义序列x ={x i }的检验统计量

m

(式(11))的x 分量，用四阶Runge-Kutta 法积分方程组，选择

初始值为[x , y , z ]=[−1,0,1]，参数值为[σ, r , b ]=[16.0,4.0,

S (m , N , r , t ) =C (m , N , r , t ) −C (1,N , r , t ) (3) 45.92]，积分步长h =0.01，积分区间为[0, 1000]。

式(3)式的计算过程为：将时间序列x 平均分解成t 个互不相交的子序列，t 为重构时间延迟，即

⎧x (1)={x 1, x t +1,..., x ⎢⎣N t ⎥⎦−t +1}⎪⎪⎪⎪试验过程中，我们选取了[53001, 56000]的3000个点，x (2)={x 2, x t +2,..., x ⎢⎣N ⎥⎦−t +2}⎪⎪ (4) ⎨⎪用C-C 方法对混沌系统x 分量的重构结果和Kim 在文[7]中... ... ... ⎪⎪⎪⎪的一样，如图1所示。 ⎪⎩x (t ) ={x t , x t +t ,..., x ⎢⎣N t ⎥⎦}

⎧dx dt =−σx +σy ⎪⎪⎪⎪⎨dy dt =−xz +rx −y (11)

⎪⎪⎪⎪⎩dz =xy −bz

1.3 C-C方法的不足 计算式(3)式定义的统计量采用分块平均的策略，即

1t

在进行系统仿真的同时，我们分别在不同的区间选取了S 1(m , N , r , t ) =∑[C s (m , N , r , t ) −C s m (1,N , r , t )] (5)

t s =1

3000个点，进行τ和τ的计算，结果见表1。

d

W

• 2970 •

2008年6月 秦奕青，等：非线性时间序列的相空间重构技术研究 Jun., 2008

图2 C-C方法分析的局部极小点与全局最小点

互矛盾的结论。因此本文认为，将S 1(t ) 的第一个零点视为最优时间延迟τd 是不合适的，只需考虑∆1(t ) 的第一个局部极小点作为最优时间延迟τd 。

图1 C-C方法分析Lorenz 混沌系统的x 分量 表1 C-C方法重构Lorenz 系统x 分量的计算结果 样本区间

嵌入维数m

时间延迟τd

10 10 10 11 11 11 11 10 11

嵌入窗口τW

191 132 184 104 137 137 84 152 101

10001--13000 21 20001--23000 15 30001--33000 20 40001--43000 11 50001--53000 14 60001--63000 14 70001--73000 9 80001--83000 17 90001--93000 11

针对以上C-C 方法的不足，本文提出了基于C-C 方法的改进的相空间重构方法。

2 改进的方法C-C-2

2.1 C-C-2方法计算原理

改进后的C-C-2方法的策略仍然根据原C-C 方法的基本策略：先定义关联积分，再构造统计量S (m , N , r , t ) ，依据

BDS 统计结论确定参数的合适取值范围。实际计算中利用S 1(m , N , r , t ) ∼t 的统计结论，实现最优时延τd 与最优延迟时间窗口τW 的估计。

因为非线性动力系统的一个变量的变化跟此变量与系统的其他变量的相互作用密切有关，改进后的C-C-2方法采用以下方法来改进关联积分：根据范数的等价性，利用2-范数代替原算法的∞-范数，其度量体现了各维的共同作用结果，而不是仅以其中绝对值最大的某一维作为标准，以更好的体现重构的相空间X 序列中各点的相关性质。

而且考虑到混沌时间序列在周期N 点的振荡，为了能寻找出更准确的混沌序列的固有属性，采取以下两个措施改进关联积分的参数r ：

通过深入分析，C-C 方法存在以下不足：

(1) 严格来讲混沌系统不存在周期，对于存在周期N 的低维混沌系统来讲，混沌时间序列在周期N 的点的数值均是存在着振荡的。式(2)定义的关联积分表示相空间中任意两对通常的周期函数统计的结果很点之间距离小于r 的概率，

明显，但对于混沌周期N 点的数值振荡，固定的r 值并不能有效体现出这种混沌的特性。

(2) 动力系统的一个变量的演化自然跟此变量与系统的其他变量的相互作用有关，作为重构的相空间X ={X i }序列中各点的相关性质度量值，式(2)定义的关联积分采用∞-范数进行计算，只考虑了m 维向量中度量最大的一维，显然是欠妥当的。

(1) 把搜索半径r 所依据的度量σ根据时间序列的变异程度进行适当扩大，目的是减少数值振荡所带来的干扰性，令σ=std(x ) *(1+cv /3) ，其中cv 是时间序列的离散变异系数(coefficient of variation)，cv =std(x ) mean(x ) ；

(3) 理想情况下S 1_cor (t ) 的全局最小点即是最优延迟时实际中存在若干个局部极小点与全局最小点在数间窗口τW ，

值上相当接近，干扰了全局最小点的判读；甚至最优延迟时间窗口τW 所对应的不是全局最小点，最终导致最优延迟时间窗口τW 的错误估计。如图2所示，所标记的几个点在不同的区域计算中都有可能是最优延迟时间窗口τW 。

(2) 将式(2)中的固定值r 改为和嵌入维数m 相关的

r (m ) =r *log(m +1) ，随着重构的维数增加适当地扩大搜索半径，目的也是为了减少振荡，特别是高维数据振荡产生的干扰。

其余的统计量均参照原算法定义进行计算。

通过大量的实验研究，发现在计算S (m , N , r , t ) 的结果分S 2_cor (t ) 存在着明显的具有周期N 的局部极小峰值，析图中，

(4) 实际计算中S 1(t ) 的第一个零点并不等于∆1(t ) 的第一个局部极小点。混沌系统虽然不存在周期，然而对于存在周期N 的低维混沌系统来讲，平均轨道周期T 是指混沌吸引子在永不重合而又彼此相似的相空间轨道上振荡的平均周期。对于平均轨道周期为T 的时间序列，当t =kT ，k =1,2,... 时，该零点很有可能既是S 1(t )

的第一个零点，也

而且这些时间点均为可能成为S 1_cor (t ) 最小值的局部极小值，因此我们给出了新的最优延迟时间窗口τW 的判断规则：原C-C 方法选择S 1_cor (t ) 在最优延迟时间窗口τW 的选择上，

的全局最小点，改进后的C-C-2方法主要考虑S 2_cor (t ) 的在计算结果图上明显的具有周期N 的第一个局部极小峰值，来确定最优延迟时间窗口τW ；在没有明显的具有周期N 的

是∆1(t ) 的零点，因是S 1_cor (t ) 的全局最小点，从而得到相

• 2971 •

2008年6月 系 统 仿 真 学 报 Jun., 2008

结果中，选择S 2_cor (t ) 的全局最小点来确定最优延迟时间窗口τW 。

另外，C-C 方法寻找1(t ) 的第一个零点或∆1(t ) 的第改进后的C-C-2方法一个局部极小点作为最优时间延迟τd ，

只用∆2(t ) 的第一个局部极小点作为最优时间延迟τd 。

的3000个点) 。

C-C-2方法的具体计算过程如下。 改进的嵌入时间序列的关联积分为

C 2(m , N , r (m ), t ) =

2

∑θ(r (m ) −d ij ) (12)

M (M −1) 1≤i

其中m 为嵌入维，N 是时间序列的数据个数，r (m ) 为计算m =1；r (m ) = 中所取的搜索半径，r (m ) >0，r (m ) =r , if

r *log(m +1), if 1m >，t 为时间延迟，M =N −(m −1) t 表

示m 维相空间中嵌入点数目，θ为Heaviside 函数：θ(x ) =0, if x

图3 原C-C 方法分析Lorenz 系统x 分量的结果图

2-范数。

S 2(m , N , r (m ), t ) , S 2(m , r (m ), t ) , ∆S 2(m , t ) , 2(t ) , ∆2(t ) 和S 2_cor (t ) 的定义均参照式(5)、式(6)、式(7)、式(8)、式(9)和式(10)，将r 改为r (m ) 。

分析计算结果时，寻找S 2_cor (t ) 的具有周期N 性质的第即平均轨一个局部极小峰值即可获得最优延迟时间窗口τW ，道周期的最优估计；在没有明显的具有周期N 的结果中，选择S 2_cor (t ) 的全局最小点来确定最优延迟时间窗口τW ；同时用∆2(t ) 的第一个局部极小点作为最优时间延迟τd 。

图4 C-C-2方法分析Lorenz 系统x 分量的结果图

2.2 C-C-2方法的仿真试验

为了验证算法的有效性和通用性，我们进行了大量的仿真试验，实验证明了我们改进的C-C-2方法具有很高的应用价值。

实验中，我们同时考察了Lorenz 混沌系统、Duffing 混

据图4，我们可以看到当t =46,92,138,183时，S 2_cor (t ) 均出现了明显的极小峰值，与图3相比较，这几个t 点S 1_cor (t ) 也均为局部极小值，因此我们可以得出一个非常重

要的结论：C-C-2方法可以通过S 2_cor (t ) 的明显的具有周期

沌系统(式(13))、Rossler 混沌系统(式(14))的x 分量。 N 的第一个局部极小峰值来确定最优延迟时间窗口τW 。

⎧dx dt =y ⎪表2说明了改进后的C-C-2算法的计算结果。我们同样⎪⎪⎪2

⎨dy dt =−δy −ax (1+x ) +f cos z (13) 采取了与原C-C 方法相同的测试范围。 ⎪⎪⎪=dz dt ω⎪⎩

表格2 C-C-2方法重构Lorenz 系统x 分量

⎧dx dt =−(y +z ) ⎪⎪样本区间 嵌入维数m 时间延迟τd 嵌入窗口τW ⎪⎪ (14) ⎨dy dt =x +d *y ⎪10001--130006 10 46 ⎪⎪dz dt e z (x f ) =+−⎪⎩20001--230006 10 46

Lorenz 系统测试条件同前，Duffing 系统的参数值[δ, a , f , ω]=[0.05,0.5,7.5,1]、Rossler 系统的参数值[d , e , f ]=[0.2,0.4,5.7]，分别用四阶Runge-Kutta 法积分方程组，选择初始值[x , y , z ]=[−1,0,1]，积分步长h =0.05，测试区间均选择[50001, 53000]的3000积分区间为[0, 5000]，个点。

为了体现出周期N 的性质，Duffing 系统、Rossler 系统的仿真试验中，令t =1,2,...,300。

通过对C-C-2方法得出的计算结果图进行分析，我们发现相空间重构图中普遍存在着具有周期N 规律的局部极小峰值，具体如图3, 4所示(测试样本区间选取[60001, 63000]

30001--3300040001--4300050001--5300060001--6300070001--7300080001--8300090001--93000

6 6 6 6 6 6 6

10 46 10 46 10 45 10 46 10 46 10 46 10 46

通过对表2分析，我们可以看到，C-C-2方法计算得出

C-C-2方法计算得出的的τd 与原C-C 方法得出的结果一样，

τW 却与原C-C 方法得出的结果不同，而是一个稳定的数值，因此通过计算得出的最佳嵌入维数也相应稳定和减小，而且与公认的Lorenz 关联维数的理论值(2.06)[10]所计算出的最佳

• 2972 •

2008年6月 秦奕青，等：非线性时间序列的相空间重构技术研究 Jun., 2008

嵌入维数相同，为混沌序列的预测、分析等后续工作降低了难度。

用C-C 方法和改进的C-C-2方法对Duffing 系统、Rossler 系统x 分量的重构结果分别见图5~8。结果同样表明C-C-2方法可以确定更准确、稳定的最优延迟时间窗口和最优时间延迟，说明了该方法具有较强的通用性。

延迟的选择更准确，而且大大增强了原算法的鲁棒性；利用

C-C-2方法分析的结果，为下一步混沌序列的预测、分析等工作奠定了更好的基础。

图8 C-C-2方法分析Rossler 系统x 分量的结果图

图6 C-C-2方法分析Duffing 系统x 分量的结果图

参考文献：

[1]

Takens F. Detecting Strange Attractors in Turbulence [C]// Dynamical

Systems and Turbulence, Lecture Notes in Mathematics. Berlin: Springer-Verlag, 1981, 898: 366-381.

Packard N H, Crutchfield J P, Farmer J D, et al. Geometry From a Time Series [J]. Physical Review Letters (S0031-9007), 1980, 45(9): 712-716.

Small M, Tse C K. Optimal Embedding Parameters: A Modeling Paradigm [J]. Physica D: Nonlinear Phenomena (S0167-2789), 2004, 194(3-4): 283-296.

李夕海, 刘代志, 张斌, 等. 基于重采样的混沌时间序列相空间重构研究[J]. 信号处理, 2006, 22(2): 248-251.

张雨, 任成龙. 确定重构相空间维数的方法[J]. 国防科技大学学报, 2005, 27(6): 101-105.

吕小青, 曹彪, 曾敏, 等. 确定延迟时间互信息法的一种算法[J]. 计算物理, 2006, 23(2): 184-188.

Kim H S, Eykholt R, Salas J D. Nonlinear Dynamics, Delay Times, and Embedding Windows [J]. Physica D (S0167-2789), 1999, 127(1-2): 48-60.

吕金虎, 陆君安, 陈士华. 混沌时间序列分析及其应用[M]. 武汉: 武汉大学出版社, 2002.

Grassberger P, Procaccia I. Measure the Strangeness of Strange Attractors [J]. Physica D (S0167-2789), 1983, 9(1-2): 189-208. 择[J]. 西安交通大学学报, 2004, 38(4): 335-338.

3 噪声效应

为了检验新算法的抗噪性能，我们在原Lorenz 时间序列上加上了Gaussian 白噪声进行测试。令y i =x i +ησεi ，其中x i 是原Lorenz 序列，σ为其标准差，εi 是均值为0、标准差为1的独立同分布的Gaussian 白噪声序列，η为噪声强度。选择η=0.1,0.2,...,0.9,1.0进行仿真试验，采取Kim 在发现新的C-C-2方法的抗噪性能由文[7]中使用的误差标准，原C-C 方法的30%提高到80%。

[2]

[3]

[4] [5] [6] [7]

4 结论

本文通过深入研究相空间重构的C-C 方法，针对该方法所存在的4点不足，提出了一种改进的C-C-2方法。新方法主要通过计算改进的关联积分，提出了新的最优延迟时间窗口的判断方式，并修改了原算法最优时间延迟的判断方式。仿真实验表明，该方法能明显体现出混沌序列周期N 的特性，对最优延迟时间窗口的选取更准确、稳定，对最优时间

[8] [9]

[10] 马红光, 李夕海, 王国华. 相空间重构中嵌入维和时间延迟的选

• 2973 •