关键词 GP算法，无标度区，关联维，时间序列

1 引言

在科学、经济、工程等许多应用中都存在着在历史数据的基础上预测未来的问题。时间序列预测是预测领域内的一个重要研究方向，这是一种根据历史数据构造时间序列模型，再把模型外推来预测未来的一种方法。近年来，来自天文、水文、气象等领域如太阳黑子、径流量、降雨量等时间序列都被发现含有混沌特性。

要将混沌理论应用于一个具体的混沌时间序列预测中，首先进行混沌识别，然后借助于相空间重构理论，并应用神经网络等方法，在重构的相空间中分析该动力系统的运动规律，进行混沌预测研究。

相空间的重构理论，即把具有混沌特性的时间序列重建为一种低阶非线性动力学系统。通过相空间重构，可以找出隐藏在混沌吸引子中的演化规律，使现有的数据纳入某种可描述的框架之下。在重构相空间时，合适的延迟时间和嵌入维数的选取至关重要。1983年Grassberger和Procaccia建立了从时间数据序列计算系统奇异吸引子关联维数的算法，该方法同时可以对时间序列进行混沌识别。在GP算法计算过程中，通常需要确定无标度区，无标度区的准确性将直接影响分形维的准确性。

无标度区的确定通常有以下几种方法^[1]：(1)运用经验公式。由于公式本身含有分形维，因而需要一个迭代过程.另外这种方法还缺乏一定的客观标准。(2)用三段直线逼近，取中间一段为无标度区。这种方法计算量太大，且取中间一段也似乎存在缺乏一定的理论依据问题。(3)Yyokoya等提出的利用最大线性度法。这种方法在应用中也存在一定的局限性。(4)将拟合直线分成两部分，前一段称为纹理分维，后一段称为结构分维。这种细分的意义和可信度有待进一步的研究。

党建武等基于置信度和相关度两个评价指标提出了无标度区确定的分组递归计算机识别算法^[2]，但对于具体的系统，并没有给出无标度区确定中标度r始值的选择方法。

杨绍清等根据J.P Eckmann等^[3，4]的成果稍作变换得到了一种比较准确的在无标度区确定中标度r始值的选择方法^[1]。但是该方法需要预先求出构吸引子的直径R、2 阶熵K₂、延迟时间τ等量，并且有个致命的缺点，还需要事先知道嵌入维m的大小，所以不宜在GP算法中使用。

本研究基于GP算法，建立了求混沌时间序列关联维数的算法，在该算法中，可手动修改参数，计算简单，无需求其它量，非常容易得到无标度区的上、下限估计值，较好地解决了GP算法中无标度区的确定问题。

2 GP算法的描述

从一维时间序列直接计算关联维的一般方法是采用Grassberger和Procaccia在Whitney的“嵌入定理”和Packard的重建相空间理论基础上提出的算法^[5]，即GP算法。

考察m维相空间中的一对相点：

（1）

设它们之间的距离，即欧氏模为 r_ij（m）显然 r_ij（m）是相空间维数m的函数，即：

（2）

给定一临界距离r，距离小于r的点对数在所有点对中所占比例记为关联积分C(r，m)：

（3）

式中，N为总相点数，H(.)为Heaviside函数，定义如下：

（4）

C(r，m)是一个积累分布函数，它描写了相空间中的吸引子上两点之间距离小于r的概率，刻画了相对于相空间某参数点x(t )在r内的相点聚焦的程度，所以存在下面的标度关系：

（5）

则称D为关联维数，即

（6）

当相空间维数m达到饱和时，通过公式(6)计算所得的D不再变化，它就是动态系统吸引子的分维，所对应的m就是我们需要选取的最佳嵌入维数m。

传统的GP算法有很多弊病，其中很突出的一点就是运算工作量大，计算速度非常缓慢。朱家富等提出了一种改进算法^[6]，可以把公式(3)改写为：

（7）

其运算工作量约为N²/2，比传统GP算法的工作量减少了一半。

3 无标度区的确定

从公式（2） 很明显得到r不能大于 ，如果r大于 ，一切点对的距离都不会超过它，则C(r，m) =1。若r小于 ，则距离 都比r大，则C(r，m) = 0。所以，对于每个m，首先计算出所有的重构后两个点之间的距离 ，在matlab编程中可用d(i，j)存储其值，然后在计算出 和 （ r_ij（m）=0，不能取）。 这时可以用 作为理论上 的估计值，用作为 的估计值。 其中 ，K可在10到20之间取值，这样可以确定无标度区的近似取值范围，并且不增加时间复杂度，因为在后面计算 的时候，需要计算 的值，而此时直接调用前面计算好的 即可，此时不需要计

算，前面的计算发挥了作用，所以前面的计算 没有重复计算之嫌，确定了无标度区的近似取值范围后，然后可以用如下三种方法得到标量r的取值序列：（1）等距离间隔取值；（2）等比例递增取值；（3）等指数比例递增取值。确定r的取值序列后，再根据具体的每个m来手动调节

和 ，以便得到一个精确的无标度区。

在具体的计算中，对于一个时间序列，首先确定无标度区的近似取值范围，从上面三种方法中选取一种合适的方法得到标量r的取值序列（标量r的取值序列长度L可以自行按精度需要设定，一般取6~10即可），然后可以用最小二乘法计算r从 递增到 的斜率，并且对于每个m，观察lnr~ 曲线，如果发现出现平坦区，则可以再次手动调节 ， 的取值，缩小区间范围。比如采用方法1等距离间隔取值得到标量r的取值序列（L=9），如果对于某一个m，发现在lnr~ 曲线右边出现平坦区，见图1所示。

图1 出现平坦区的lnr~lnc(r,m)曲线

显然，此时第8个、第9个r(即r8、r9)取得过大， 和 之间并非是真实的无标度区间，并且这时用最小二乘计算出的D也不是准确的，所以这时可以把 调整为 ，然后再次按等值递增重新计算得到标量r的取值序列，再重新计算D，见图2所示。

图2 调整取值区间后的lnr~lnc(r,m)曲线

对于每个出现平坦区的m都做如此处理，在实验中发现，如果m=M要做如此变换，一般的m=M+1，m=M-1都要做如此处理。如果发现在曲线左边出现平坦区，那么用同样方法修改 即可，采用方法2、方法3得到标量r的取值序列，也可做类似处理。用该方法确定无标度区后，当m增大到一定值的时候，则看到D不再增大，意味着关联维数出现饱和现象，此时可用D作为关联维数的近似值，而此时的m则可当作嵌入维数。

4 实例仿真

文献[1]对Henon和Lorenz两种著名的混沌系统关联维数进行了计算，并与文献[4]的结果进行了比较，本文也采用这两种系统进行仿真，对于Henon系统和Lorenz 系统，分别取时间序列500和1000点，用于实现关联维数提取，计算值示于表1，表中参照值取自文献[1，4]。仿真结果表明本文方法能够准确提取混沌时间序列的关联维数，并且无需计算其它量。

表1 两个典型人工混沌时间序列无标度区和关联维计算表

5 结束语

文中提出的基于GP算法的关联维计算中无标度区的识别算法，避免了尺度选择盲目性，减少了计算冗余，计算效率高并且具有可以在不计算直径R、2 阶熵K₂条件下提取关联维的特点。无论实测时间序列无标度区跨度大小如何变化，均可以快速有效地确定混沌时间序列的无标度区和提取其关联维参数，从而提高处理效率。

参考文献

[1]杨绍清，章新华，肖明杰等.一种实用的混沌信号相关维的提取算法[J].电子学报，2000，28(1)：20~22

[2]党建武，施怡，黄建国.分形研究中无标度区的计算机识别[J].计算机工程与应用，2003，25(3)： 25-27

[3] J.P. Eckmann，D.Ruelle. Fundamental Limitations for Estimation Dimensions and Lyapunov Exponents in dynamical Systems[J]. Physica D，1992，56： 185-187

[4] Y.C. Lai. Effective Scaling Regime for computing the Correlation Dimension from Chaotic Time Series[J]. Physica D.，1998，115： 1-18

[5] Grassberger P，Procaccia I.Measuring the strangeness of strange attractors[J]. Physica D，1983，9： 189~208

[6]朱家富杨，杨浩，彭拥军.一种计算脑电信号相关维数的改进算法[J].西南师范大学学报(自然科学版)，2004，29(4)：623~626