《最大李指数计算方法.doc》由会员分享，可在线阅读，更多相关《最大李指数计算方法.doc（9页珍藏版）》请在金锄头文库上搜索。
1、物理学报ACTA PHYSICA SINICA2000Vol.49No.4P.636-640一种最大李雅普诺夫指数估计的稳健算法杨绍清章新华赵长安最大李雅普诺夫指数是诊断和描述动态系统混沌的重要参数.在深入研究相空间重构技术和轨道跟踪法的基础上，提出了一种从标量混沌时间序列中估计最大李雅普诺夫指数的新算法.该算法能够克服现有算法的不足，主要有以下三个优点：1）很高的精度；2）几乎不受噪声的影响；3）所需的计算时间和存贮空间小, 能进行在线计算.PACC： 0545A ROBUST METHOD FOR ESTIMATING THE LARGESTLYAPUNOV EXPONENTYANG SH
2、AO-QING(Harbin Institute of Technology,Harbin150001,China)ZHANG XIN-HUA(Dalian Naval Academy,Dalian116018,China)ZHAO CHANG-AN(Harbin Institute of Technology,Harbin150001,China)ABSTRACTThe largest Lyapunov exponent is an important parameter of detecting and characterizing chaos produced from a dynami
3、cal system. In this paper, based on the technology of phase space reconstruction and the methods of trajectory tracing, a new algorithm is proposed for estimating the largest Lyapunov exponent from a scalar chaotic time series. This method, which can overcome the deficiencies of the existing methods
4、, has three main advantages: (1) It has highly accurate results; (2) It is little affected by noise; (3) It only needs a little time of computation and small space of memory and can calculate the largest Lyapunov exponent on line.1引言近几年来，混沌信号的诊断及其特性的描述已经广泛地应用于时间序列的分析中1.在诊断和描述混沌信号时，最大李雅普诺夫指数(1)不仅是一个很
5、重要的不变量，而且是判断混沌存在的一个重要依据.因此，利用1去诊断和描述混沌仍然是一种主要方法.其他有些方法(如K-熵法等2)与这种方法没有本质的区别.目前估计1的方法主要有两种3：一种是分析法(analytic approach)，一种是轨道跟踪法(trajectory tracing method).前者是用一个函数（如局部多项式或神经网络等）来建立系统模型，然后估计系统的雅可比矩阵，进而求取1；后者是直接从1的定义出发跟踪系统的两条轨道，获取1.由于轨道跟踪法不像分析法那样易受系统拓扑结构的影响，因而受到许多学者的高度重视.自从1985年Wolf4提出轨道跟踪法以来，这种方法得到了较大的
6、改进.最有代表性的是Rosenstein等5的工作,他们给出了一种比较优化的算法: 首先绘制 相对j。t的曲线，然后通过计算曲线中直线段的斜率来获得1.这里dk(j)是第k对最近邻沿各自的轨道经j步后的距离，表示对所有k的平均，而t为对时间序列的采样周期.这种算法由于充分地利用了能够利用的数据，并对它们进行了某种意义上的平均，因而结果比较稳定，但却大大地增加了计算量.总之，目前轨道跟踪法仍然存在以下三个问题：1）精度不高；2）受噪声影响大；3）计算量大、所需存贮空间大、收敛速度慢、很难应用.在应用轨道跟踪法时，首先要重构相空间.自从1981年Takens6提出了嵌入定理以来，相空间重构技术得到
7、了很大的发展.但是仍然存在着计算复杂，整体性差等缺点.本文在Kugiumtzis等工作的基础上715，提出了一套简单易行的相空间重构方法，并针对轨道跟踪法的现有问题提出了一种新的计算1的方法.在仿真试验中, 本算法的计算结果表现出了很强的鲁棒性.本方法计算量小、所需存贮空间小、收敛速度较快、能实现对1的在线计算，因而是一种很实用的好方法. 2本文的改进算法本文算法的主要思想是：在计算最大李雅普诺夫指数时，不能忽视混沌吸引子中的收缩作用（the action of contraction）.但至今，混沌吸引子中的这一重要作用并没有得到高度的重视，因而也未能在现有的算法中体现出来.混沌吸引子是耗散
8、系统的专有特征，是扩散（expansion）与收缩两种作用互相牵制的结果.为了维持混沌吸引子的存在，总的收缩率必然要超过总的扩散率.这样，即使在系统具有几个正的李雅普诺夫指数时，如果对整个李雅普诺夫指数谱进行求和，所得结果仍然是负值.因此，在运用轨道跟踪时，不能用对某一方向或所有方向上的李雅普诺夫指数进行平均的方法来求取1，而应该每次都在最大方向上估计1，然后再对所有估计值进行平均求取1.2.1相空间重构重构的相空间是否贴近实际未知系统，对系统特性的进一步刻划具有较大影响.对于一N点标量时间序列v(t0+k。t):k=0,1,N-1, 可以用Takens嵌入定理去重构相空间RmXi=(x(ti
9、), x(ti+p。t), x(ti+(m-1)p。t)i=1,2,3,Mx(tr)=v(t0+(r-1)。t)r=1,2,N. 这里Xi是重构相空间Rm中M点重构轨迹中的第i个点，而M=N-(m-1)p ，m是嵌入维，= p。t是时间延迟, w=(m-1)是时间窗，其中t为采样周期.可见，在重构相空间时只要选择m，和w中的任意两个参数即可，另一个参数可由w=(m-1)直接求得.研究中发现，可选择原始时间序列的平均峰值时间(mtbp)作为w .另外，w的大小有一个范围，w必须满足:w minww max，否则重构将会变坏，其中w min可用w min=（12）1/2(x)/(dx/dt)来估计
10、8，其中(。)为均方差，对Lorenz-1（见表1）来说，w min0.67 (N=4000时的结果).不难看出mtbp=4(x)/(dx/dt) ，因此mtbp是始终大于w min的.对Lorenz-1 系统来说，本文选择w=0.78，此结果与上述这些分析是完全一致的.在选择w以后，还要对或m作出选择.由于的大小与信号自身的相关程度有关，相关程度越低，就应越小，因此对带有噪声的信号重构时，不能太大.另一方面，随着m的增大，也应减小911，当m足够大时，就可以选得很小.因此，本文选择最大的m而得到最小的，即=t，p=1.本文同时注意到Schouten等12在估计带有噪声的吸引子分形维时也采用了
11、这种用法.综上所述 ，本文采用下列方法来重构相空间：1) 求取原始时间序列的mtbp，让w=mtbp（可以让w稍大于mtbp）.带有噪声的时间序列可经滤波后再求取mtbp;2）选择p=1，即=t;3） m=w/+1.2.2计算1假设重构轨迹中的第i个点Xi的最近邻为Xip，即满足：Xip-Xi=minj Xj-Xi,。为欧氏范数.在这里需要对两点的时间间隔作出规定，假设w/t|j-i|Rw/t，1Rw/t可以保证Xi和Xip是在不同的轨道上.而|j-i|Rw/t可以保证一次计算1所需的点数很少，这样不仅节省了大量的计算时间，而且大大地减少了计算时所需的存贮空间.R可以选择为一个常数，用以避免最
12、近邻之间的距离为零，不可太大和太小.R太小会使最近邻之间的距离受到噪声的影响很大，而R太大又会使两个最近邻完全失去相关性，同样可能会出现距离为零的情况.为了能在扩散的最大方向上计算1(i)，可用下式来求1(i): (1)这里需要附加两个条件由于时间相隔小于w时的两点可以认为是在同一条轨道上，所以条件（i）保证了由最近邻演化后的两点仍然分别保持在各自的轨道上，其中a也为常数.然而虽然有了条件（i）的限制，但还不能保证Xip+jhXi+jh与XipXi的方向是基本一致的，用方向上相差比较大的Xip+jhXi+jh与XipXi来计算1(i)同样是没有意义的.条件（ii）就是为了防止在有较大扭曲的轨道
13、中计算1(i)，m越大扭曲的程度越大，因而,限制的角度应相应增大,也即限制的条件应适当放松，但最大不能超过80.条件（ii）使(1)式对1(i)的计算更加合理.由此可见，利用(1)式及其两个附加条件计算的1(i)真正代表了系统局部最大的李雅普诺夫指数，其中参数b与嵌入维m成正比,它们可以用一简单的线性关系来表示. 最后，可得表1本文算法的计算结果系统t1理论值SNR1计算值|误差|%LogisticXi+1=4xi(1-xi)10.693-500.398342.5-40 0.4064 41.4 -20 0.3960 42.9 -150.4046 41.6 -10 0.4030 41.8 -5
14、0.3947 43.0 0 0.4160 40.0 5 0.4210 39.2 10 0.4378 36.8 1000.7060 1.9 Hnonxi+1=1-1.4x2i+yiyi=0.3xi10.418-50 0.3340 20.1 -40 0.3225 22.8 -20 0.3260 22.0 -15 0.3310 20.8 -10 0.3300 21.1 -5 0.3180 24.9 0 0.3220 23.0 5 0.3046 27.1 10 0.3054 26.9 100 0.3960 5.3 Lorenz-1=-(8/3)x+yz=-10y+10z=-yx+28y-z0.012.06-50 2.9440 42.9 -40 2.6612 29.2 -20 2.8115 36.5 -15 2.4082 16.9 -10 2.2550 9.5 -5 1.9628 4.7 0 2.0560 0.2 5 1.7986 12.7 10 1.7600 14.6 100 2.0690 0.4 Lorenz-2=16.0(y-x)=x(45.92-z)-y=xy-4.0z0.011.50-50 0.8496 43.4 -40 0.8602 42.7 -20 0.81
