混沌时间序列预测方法中的嵌入延迟确认

嵌入延迟τ的选取方法较复杂。实际的观测序列，因为存在噪声干扰和估计误差，如果τ太小，相空间轨迹会向同一位置挤压，信息不易显露，产生冗余误差；如果τ太大，会导致某一时刻的动力学性态与后一时刻的动力学性态变化剧烈，使得简单的几何对象表现得很复杂，动力系统信号失真，产生不相关误差。仍以Lorenz吸引子为例，画出其在二维相空间中，不同延迟时间下不同的演化轨迹如图2.8所示，可以看出奇异吸引子形状的变化。

图2.8　不同延迟时间下Lorenz吸引子的二维相空间重构

（a）τ=2；（b）τ=8；（c）τ=18；（d）τ=40

现有的延迟时间选择方法一般基于两个准则：

（1）相空间扩展法[34，35]，目的是使重构相空间轨迹从相空间的主对角线尽可能的扩展，但又不出现重叠，如平均位移法、摆动量法等。

（2）序列相关法[36～38]，即降低X（t）内元素的相关性，同时保持X（t）中包含的原动力学系统的信息，如自相关法，互信息量法，高阶相关法等。

准确求取时间延迟的方法应具有如下优点：

（1）概念清晰，在相空间几何学中有明确的意义。

（2）数学表达简单易算。

（3）对数据量要求不大，较小数据量求解结果的可靠性仍较高。

（4）对含噪数据集合计算结果的鲁棒性较高。

（5）易于扩展到高维等。

下面介绍几种常用的方法。

1.自相关法

自相关法为也称自相关函数法，是一种比较成熟的求时间延迟τ的方法，此方法主要是提取时间序列间的线性相关性。一般来说，对于一个混沌时间序列，可以先写出其自相关函数，然后做出自相关函数关于时间t的图像。根据实验测得的结果，当自相关函数下降到初始值的（1—1/e）倍时，所得到时间t即为重构相空间的时间延迟τ。

对于连续变量x（t），其自相关函数C（τ）定义为

式中：τ为时间的移动值，表示两时刻t和t+τ运动过程的相互关联或相近似的程度。

当x（t）的幅值一定时，C（τ）越大，则意味着x（t）与x（t+τ）关联越大。当τ变小时，两个时刻的间隔也变小，两个运动过程的关联程度变大；反之，当τ变大时，两个运动过程的关联程度变小，最后趋近于0。

对于离散混沌时间序列x（1），x（2），…，x（t），…，序列的时间跨度为jτ的自相关函数为：

由此可固定j，做出自相关函数关于时间τ（即取τ=1，2，…）的函数图像，则自相关函数下降到初始值的（1—1/e）倍时，所得到的时间τ即是重构相空间的时间延迟τ。

2.改进自相关法

自相关法是较成熟的延迟时间求取方法，但其只能提取序列间线性相关性。而且，根据自相关法得到的τ，可分别让X（t）和X（t+τ），X（t+τ）和X（t+2τ）之间最小相关，但X（t）和X（t+2τ）之间却仍然可能存在较强的相关性，很难推广到高维。Luis[39]发展了自相关法的思想，引入一个非线性相关函数，使延迟时间τ的计算度量变成两部分：

（1）采用线性相关函数。(www.xing528.com)

检测状态之间的线性相关性。

（2）采用非线性相关函数。

检测状态之间的非线性相关性。设τx和τx2分别对应Φxx（τ）和（τ）的第一个极小值，则

该方法同时检测状态之间的线性和非线性关系，（τ）可以检测到Φxx（τ）检测不到的时间序列中的变化，一般会得到比单独检测线性相关性更短的延迟时间。

3.互信息法

互信息函数是两个随机变量间一般性（包括线性和非线性）随机关联的度量，是非线性的分析方法，在延迟时间的选取上要优于自相关函数法。互信息法[40～42]利用计算互信息函数的第一极小值来确定最佳延迟时间τ。

令﹛x（i），i=1，2，…，N﹜表示一组信号，设点x（i）的概率密度为Px[x（i）]，它将信号映射为概率。令﹛y（j），j=1，2，…，N﹜，表示另一组信号，在y（j）处的概率密度为Py[y（j）]。两组信号同时测量得x（i）、y（j）的概率表示为Pxy[x（i），y（j）]，又称作联合概率。对于两组信号﹛x（i），y（j）﹜，给定x（i）的一个测量值，预测y（j）的平均信息量为互信息函数

其中，H（x）是信号﹛x（i）﹜的熵，表示对指定系统的N 个x（i）测量得到的平均信息量。H（y）定义与H（x）类似，H（x，y）是联合熵。

利用实验数据计算互信息函数，其主要工作是计算Pxy[x（i），y（j）]、Px[x（i）]和Py[y（j）]。通常采取划分网格的方法，将变量组成的样本空间划分为若干网格，然后通过统计各个格子中的点数来求出其概率值。格子的划分可以有不同的方法，Fraser等使用了等概率划分空间网格的方法，由于此法网格数很难确定，在实际的计算中并不实用。这里使用等间距划分网格的方法[41]计算互信息，将样本空间分别按x（i）和y（j）组成的可能状况值等间隔地划分若干组。显然，落入各个盒中的样本点数将与变量x（i）和y（j）的概率分布存在着一定的联系。

画出[x（i），y（j）]在二维平面上的重构图形，在图形上对欲求I（x，y）值的点集区域（吸引子）画出矩形框Δ。设Δx是Δ在x方向的长度，Δy是Δ在y方向的长度，然后把Δ在等间距的划分成小格子ε，记x方向小格子的长度为εx，x方向小格子的数目为M，所以有Δx=Mεx。同理，记y方向小格子的长度为εy，y方向小格子的数目为M′，所以有Δy=M′εy。在Δ中小格子的总数为M×M′。设（a，b）是所划格子的起始点，对采样点（x，y）作下述判断：

若满足a≤x（i）≤a+Δx，b≤y（j）≤b+Δy，i，j=1，2，…，N，则点[x（i），y（j）]在Δ中，否则该点不在Δ中。

如果点[x（i），y（j）]在Δ中，再判断是否满足（k—1）εx≤x（i）—a＜kεx，（l—1）εy≤y（j）—b＜lεy，k=1，2，…，M，l=1，2，…，M′，即

如果也满足，表明[x（i），y（j）]也在小格子εk，l中，此时进行一次记录，直到将所有的数据点全部搜寻一遍，记录落入标号为（k，l）的小格子中的所有数据点Nxy，同时记录落入范围在k—1到k内的点数Nx，和从l—1到l内的点数Ny，即可得到

式中：Ntotal为全部采样点数。

将计算出的Pxy、Px和Py代入式（2.86）～式（2.88）进行计算，可以求出在给定时间序列下的互信息值I。

在时间序列重构问题中，通常要求的是x（t+τ）值对x（t）值的依赖性。因此，可以令（x，y）=[x（t），x（t+τ）]。为了得到互信息I与延迟时间τ的函数关系，先取τ值接近于0，做出重构图，在这个重构图上进行取框、划分和计算，得到一个I值；然后将τ值逐渐增加，重复上述过程，最后得到了所需的结果。

这种通过计算互信息函数第一极小值来确定延迟时间的方法效果较好，但计算方法过于复杂，无法避免大量的计算和复杂的空间划分要求。同时，得到的互信息值与所分区间的个数有关，这是互信息函数方法的不足之处。