相空间重构延迟时间互信息改进算法研究

蒋爱华，周璞，章艺，华宏星

(1.中国船舶重工集团公司第704研究所减振中心，上海 200031;2.上海交通大学机械系统振动国家重点实验室，上海 200240)

相空间重构延迟时间互信息改进算法研究

蒋爱华1，2，周璞1，章艺1，华宏星2

(1.中国船舶重工集团公司第704研究所减振中心，上海 200031;2.上海交通大学机械系统振动国家重点实验室，上海 200240)

针对改进互信息算法利于快速可靠获得时间序列相空间重构的延迟时间问题，通过等边缘分布2、4等分Lorenz时间序列构成平面分析Cellucci互信息算法缺陷;用大小顺序值代替原序列数值、判断新序列数值所在等边缘概率区间获得概率分布矩阵、修正概率分布矩阵最末行与列改进Cellucci互信息算法;以改进算法所得最佳延迟时间进行Lorenz时间序列相空间重构并以小数据量法得出其最大Lyapunov指数，对比雅可比矩阵法所得最大Lyapunov指数以确认改进算法的有效性。结果表明，时间序列长度不能整除划分区间数时Cellucci互信息算法会获得错误的最佳延迟时间;所提改进算法能消除Cellucci算法缺陷，且计算速度快于Fraser算法;数据序列长度较大时改进算法结果更稳定;由两种最大Lyapunov指数计算方法所得结果间误差较小，表明改进的互信息算法有效、可靠。

相空间重构;延迟时间;互信息;Cellucci算法;最大Lyapunov指数

对测试所得振动信号时间序列进行相空间重构，可在不清楚振动系统结构与影响参数时研究振动系统特性。因此在非线性振动系统分析中获得一定应用［1－2］。

互信息算法为相空间重构中确定延迟时间的常用方法。设{s(ti)}(i=1，2…N)为试验观测所得时间间隔为△t的一维性时间序列，该序列可重构为

式中:τ为延迟时间。

式中:H(Q)，H(S)，H(Q，S)分别为Q、S的信息熵及Q与S间互信息熵，通常以第一个互信息极小值点所在位置作为最佳延迟时间。

在计算两时间序列间互信息时常将两时间序列S与Q值对应到两相互垂直轴上，构造成SQ平面;并将每个序列中最大与最小值间区域按一定准则分别划分为若干区间si，qj，则该两序列数据对将对应于SQ平面中点。通过计算区间si中数据对点数，除以总数据对数N则可得各区间边缘分布概率。两时间序列区间相互重合区域为(si，qj)，计算该区域内数据对点数并除以总点数，则得重合区域内联合概率密度，从而可得S，Q序列间互信息。H(Q)，H(S)，H(Q，S)表达式为

式中:Pq(qi)，Ps(si)，Psq(si，qj)分别为Q在qi区域的边缘分布概率密度、S在si区域的边缘分布概率密度及S与Q在(si，qi)区域的联合概率密度。

SQ平面划分示意图见图1。Q，S间互信息化简为

图1 SQ平面划分示意图Fig.1 Skeleton map of SQ plane partition

1 互信息算法

1.1 SQ平面划分准则

划分SQ平面准则主要有等边缘概率分布与等间距两种(图1)。等边缘概率分布划分区间即使划分后区间具有相同数据对点数。第一个用互信息法确定最佳延迟时间Fraser算法即采用该方法［3］。Fraser算法按等边缘分布将SQ平面划分为4个网格，判断每个网格是否存在子结构或已稀疏。若无子结构或已稀疏则无需细分，若存在子结构则据边缘概率对含子结构区域进行等边缘分布划分，直到各区域内无子结构存在。判断是否存在子结构准则为数据对是否均匀分布于区域内。Fraser算法示意图见图2。由图2看出，区域R1(2)、R1(3)中数据点均匀分布无需进一步划分;区域R1(1)、R1(4)、R2(4，2)中存在子结构，则需进一步划分。

因此，Fraser算法计算互信息时子结构网格数按4的幂次方增加，且为确定子结构是否存在需将各区域按等边缘概率分布划分两层，判断各下层子结构网格中的数据点数。计算过程中需存储各层子结构概率分布，故Fraser算法耗时、耗空间。

图2 Fraser－Swinney算法示意图Fig.2 Skeleton map of Fraser－Swinney algorithm

等间距划分区间可使各划分区间距上下限差值相同。与Fraser算法相比该方法计算速度较且程序编译简单，但存在互信息对划分区间数极敏感，即划分区间太多则会有大量区域内数据点数为0或1，无法反映数据对的概率分布，并获得错的互信息。反之，若划分区间太少则各区域内存在子结构而获得互信息值不准确。此外，该方法所得结果不适用于对不同变量的时间序列需相同等间距划分区间数的多变量非线性时间序列重构。

图3 Rossler微分方程中三变量时间序列等间距划分下求互信息Fig.3 Mutual information of time series from the three variables ofRossler with equal－distance elements

以Rossler数值大小分布不均时间序列为例，采用相同等区间划分数划分三变量的时间序列，会得到错误的最佳延迟时间。图3(a)为运用四阶龙格－库塔所得10000点Rossler离散序列三维图，起始点为(－1，0，1)，时间步长为0.05。图3(b)、(c)分别为Rossler离散序列中三变量的数值分布柱状图、等间距40等分数值区间时延迟点数1～200的互信息值图。可见由变量Z时间序列所得最佳延迟时间不准确。

Rossler方程为

式中:d=0.2，e=0.4，f=5.7均为系数。

大量研究集中于等间距划分区间数确定。对总长度为N的序列，Mosteller等［4］提出划分个区间，Bendat等［5］提出划分1.87(n－1)0.4个区间，Rissanen［6］由算法的随机复杂性通过对划分区间进行理论推导认为，拥有最小随机复杂性的划分区间数即为最佳划分区间。随机复杂性公式为

式中:F(m)为随机复杂性，m为划分区间数;R为序列中最大最小值之差值;Ni(i=1，2…)为划分后各区间数据点数，即

由于N一般较大，式(5)、(6)中存在大数阶乘会在计算中溢出，故将式(4)后两项变换为

由于复杂度本身的计算为耗时耗空间过程，从而可抵消等间距划分节约时间与空间优点，目前尚无较好方法用于确定最佳等间距划分区间数。Cellucci等［7］提出的基于统计的等边缘划分区间互信息算法，不仅能满足多变量时间序列重构要求，且计算速度更快。

1.2 Cellucci互信息算法

互信息算法［7］建立于两序列S，Q间统计独立假设，此时S，Q形成的数据点会在SQ平面均匀分布(图1)，由式(2)知S，Q间互信息将为零。此时SQ平面内任意区域点数为已知，设S轴第i区间内点数为OS(i)，Q轴第j区间内点数为OQ(i)，则区域(i，j)内点数ESQ(i，j)为

采用等边缘概率划分区间，且S，Q用相同划分区间数NE，边缘分布密度PS(i)，PQ(j)变为PS(i)= PQ(j)=1/NE，则式(7)变为

此时若已知ESQ(i，j)，则可得NE，取ESQ(i，j)≥5，则式(8)变为

取NE为满足以上要求的最大整数。当N为NE的倍数时，将S，Q轴分别按等边缘分布概率划分为NE个区间，则式(2)可变为

当N不是NE的倍数时则按要求确定NE，即按NE划分SQ平面后，所有区域内实际点数满足OSQ(i，j)≥1;80%区域内满足OSQ(i，j)≥5。判别是否满足第二要求所用均匀分布卡方检验完成，卡方检验公式为

卡方检验中所用自由度ν=(NE－1)2，则数据对在SQ平面内均匀分布概率为

1.3 Cellucci算法缺陷

该算法缺陷为:当时间序列S或Q长度不能整除划分区间数NE时其计算的互信息出现错误。以Lorenz系统中x变量序列为例，Lorenz系统的微分方程为

式中:σ=10，b=8/3，r=28。

用4阶龙格－库塔获得100 000点离散序列，起始点为(8.331，13.291，18.063)，时间步长0.01，其中前4 096个数据点三维图见图4，该过程用Tisean3.0.0［8］完成。

图4 Lorenz吸引子三维图Fig.4 Three－dimension diagram of 4096 Lorenz data pairs whose original

S，Q长度取4 096，按等边缘分布对其划分2等份，划分后SQ平面各区域内数据点数用矩阵表示。延时点数1、10、100、200的两时间序列等份后各矩阵为

按等边缘分布对S，Q分别划分4等份，延时点数1、10、100、200的两时间序列等份后各矩阵为

由各矩阵看出，对SQ平面按等边缘分布2等份时数据对点分布矩阵中无零元素，需进一步划分确认是否存在子结构或零元素;对SQ平面按等边缘分布4等份时矩阵中出现零元素，对SQ平面进一步细分时对应的新矩阵中元素亦为零，与算法［7］中要求所有区域内数据对数大于1不符。即数据对数非划分区间整数倍时，用该算法对SQ平面等边缘概率划分将会仅二等份，致最佳延时不准确。4 096(非NE倍数)与4 500 (NE倍数)个数据点按算法［7］计算的Lorenz与Rossler延迟点数1～200的互信息见图5。

图5 不同长度数据序列Cellucci法所求互信息Fig.5 Mutual information fromCellucci's algorithm with different series lengths

此外，其它长度的Lorenz、Rossler序列亦存在相同情况，且Cellucci算法中用于确定等边缘分布划分区间数的另个判据较复杂，程序不易实现。

2 互信息改进算法

2.1 最末边缘区间分布点数处理

在Cellucci算法基础上本文提出改进算法。以满足式(9)最大整数NE为等边缘分布划分S或Q轴区间数，每个区间内数据点数取小于或等于N/NE的最大整数N1，以N1由小到大依次划分S或Q轴直至结束。若N/NE为整数，则划分后各边缘区间内数据点数完全相同;若N/NE不为整数，则划分后S，Q轴最末个边缘区间内分布点数少于其它区间。此时按最后一个区间分布点数与前各区间分布点数的比例进行修正，即设N2为N整除NE时所得余数在最后一个边缘区间内各区域分布点数得出后分别乘以N1/N2，得各区域内新的分布点数，在计算互信息时最后一个区间点总数则按N1计算，此时的概率分布矩阵用C表示。

2.2 概率分布矩阵计算

由式(2)知，互信息值大小取决于划分SQ平面后各区间的概率分布矩阵。在等边缘概率划分区间前提下分布概率矩阵只与各元素在其序列中大小顺序有关，与序列中各元素具体数值大小无关，因此在不改变各数值大小顺序前提下对两序列中数值进行任意变换，可简化SQ平面内各区域的分布概率判断。判断各区域内分布点数时所用方法为:①对S，Q两序列分别排序，并用S，Q各元素在其序列中的大小顺序值代替实际值，分别用两向量A、B表示，无论S、Q数值范围大小，A、B所含数据均为1～N的整数，判断各区域内分布点数则可直接用A、B完成;②因已知每个边缘分布区间所含点数为N1，则可依次判断数据对(A(i)，B(i))(i=1，2…N)所在区域，并将区域概率分布矩阵C中元素C(m，n)加1，其中m为1+A(i)/N1最大整数，n为1+B(i)/N1最大整数。

获得概率分布矩阵C后用第一种处理方法时式(2)变为

3 改进算法验证

为确认算法的有效性需对最佳延迟时间进行验证，用相空间重构的非线性不变量，即利用Lyapunov指数、关联维数等判断重构效果［9－10］。本文用最大Lyapunov指数作为判断依据。该指数计算有两条途径，即①已知系统微分方程或映射关系，利用该指数的定义或系统微分方程的雅可比矩阵特征值在一段时间内的平均值计算获得Lyapunov指数谱［11］，其中最大值即为最大Lyapunov指数;②已知系统试验的时间序列，则重构后数据可由小数据量［12］、Wolf［13］、BBA［14］等方法及改进算法获得Lyapunov指数，其中小数据量法与Wolf法所得为最大Lyapunov指数，BBA法可获得所有Lyapunov指数谱。

对Lorenz系统，由于微分方程已知，由①可得理论的最大Lyapunov指数。Lorenz系统雅可比矩阵为

通过4阶龙格－库塔法获得该系统中不同变量的时间序列，由时间序列可由②获得另个Lyapunov指数，所得指数值与理论指数值接近程度则可反映重构效果的好坏，即所选延迟时间与嵌入维数的好坏。

图6 不同数据长度的互信息Fig.6 Mutual information by different pointslengthes

分别用Lorenz系统中变量x的前1 024、2 048、4 096、8 192、16 384、24 576、32 768、65 536个数据点进行互信息计算。据改进算法所得1～200延迟点数时互信息见图6。计算过程中各数据点数对应参数见表1。由表1看出，用于计算互信息的序列长度较短时计算结果与稳定值有一定偏差，而随采用数据长度增加，用以上方法所得第一最小互信息点即最佳延迟时间均趋于稳定;用于计算互信息的序列长度大于等于16 384后最佳延迟点数均为17。

表1 不同数据长度计算互信息时参数Tab.1 Parameters appearing in the process of mutual information calculation by series with different lengths

以上计算互信息算法均在Matlab6.5上实现，所用计算机CPU频率3 GHz，内存2 GB。由数据长度计算时间知，改进方法的计算时间远小于Fraser算法，相同运行环境下用Fraser算法计算长4 096序列互信息费时已超200 s。

表2 两种方法所得最大Lyapunov指数Tab.2 Maximal Lyapunov exponents gained from Jacobi matrix and phase space constructed by time series with different time delay

据Takens定理，采用大于2倍于分形维数的嵌入维数构造相空间时重构的相空间与原系统非线性不变量［15］相同。已知Lorenz序列的分形维数为2.07［16］，故本文嵌入维数取5，通过产生的100 000点时间序列对相空间进行重构，由重构数据通过小数据量法计算出最大Lyapunov指数，该计算过程可直接用Tisean 3.0.0完成。用Lorenz系统方程的雅可比矩阵通过QR分解方法获得理论最大Lyapunov指数。两种方法计算值见表2。由表2看出，延迟时间取17时小数据量法所得最大Lyapunov指数与理论值相差0.276%，与最佳延迟时间18相差仅1点。

4 结论

(1)时间序列长度为等边缘划分区间数的整数倍时Cellucci互信息算法有效;反之，时间序列长度不为划分区间的整数倍时Cellucci互信息算法会获得错误的最佳延迟时间。

(2)对两时间序列进行排序，可用序列中各数值在序列中大小顺序值代替原数值实现时间序列的数值转换;判断序列各数值所在等边缘概率区间方法可得概率分布矩阵，计算速度快，程序易编译。

(3)数据序列长度较大时用改进的最佳延迟时间计算方法结果更稳定;改进互信息算法计算时间远小于Fraser互信息算法。

(4)用改进互信息算法的最佳延迟时间所得Lyapunov指数与理论值的相对误差较小，表明该算法有效可靠，可消除Cellucci互信息算法缺陷。

(5)本文互信息算法构造中均设每个区域内平均分布5个数据点，虽所得最佳延迟时间近似值较准确，但该假设是否合理尚需理论推导;是否有其它值作为每个区域内平均分布点数获得更准确结果有待验证。

［1］He Qing－bo，Du Ru－xu，Kong Fan－rang.Phase space feature based on independent component analysis for machine health diagnosis［J］.Journal of Vibration and Acoustics ASME，2012(2):10－14.

［2］Zhang Hong－wei，Li You－rong，Xiao Han，et al.Bearing fault diagnosis based on weighted phase space reconstruction［C］. Digital Manufacturing and Automation，IEEE International Conference，2010:315－318.

［3］Fraser A M，Swinney H L.Independent coordinates for strange attractors from mutual information［J］.Physical Review A，1986，33(2):1134－1140.

［4］Mosteller F，Tukey J W.Data analysis and regression［M］. Addison－Wesley，1977.

［5］Bendat J B，Piersol A G.Measurement and analysis of random Data［M］.New York:John Wiley，1966.

［6］Rissanen Y.Stochastic complexity in statistical inquiry［C］. World Scientific，Singapore，1992.

［7］Cellucci C J，Albano A M，Rapp P E.Statistical validation of mutual Information calculations:comparisons of alternative numerical algorithms［J］.Physical Review E，2005(71):1－14.

［8］Hegger R，Kantz H，Schreiber T.Practical implementation of onlinear time series methods［C］.The Tisean package，CHAOS，1999:413－435.

［9］Kantz H，Schreiber T.Nonlinear time series analysis［M］. Cambridge:Ambridge University Press，2004.

［10］Edward O.Chaos in dynamic system［M］.Cambridge: Cambridge University Press，2002.

［11］Kawabe T.Indicator of chaos based on the riemannian geometric approach［J］.Physical Review E，2005(71):1－4.

［12］Rosenstein M T，Collins J J，De Luca C J.A practical method for calculating largest Lyapunov exponents from small data sets［J］.Physical Review D，1993(65):117.

［13］Wolf A，Swift J B，Swinney H L，et al.Determining Lyapunov exponents from a time series［J］.Physical Review D，1985(16):285－317.

［14］Brown R，Bryant P，Abarbanel H D I.Computing the Lyapunov spectrum of a dynamical system from observed time series［J］.Physical Review A，1991(43):2787－2906.

［15］Takens F.Detecting strange attractor in turbulence［J］. Lecture Notes in Mathematics，1981(898):361－381.

［16］Viswanath D.The fractal property of the Lorenz attractor［J］.Physical Review D:Nonlinear Phenomena，2004(190): 115－128.

Improved mutual information algorithm for phase space reconstruction

JIANG Ai－hua1，2，ZHOU Pu1，ZHANG Yi1，HUA Hong－xing2
(1.704 rsearch Institution，China Shipbuilding Industry Corporation，Shanghai 200031，China;
2.State Key Laboratory of Mechanical System and Vibration，Shang Hai JiaoTong University，Shanghai 200240，China)

The mutual information algorithm was improved for gaining rapidly and reliably the time delay in phase space reconstrution of time series.The defect of Cellucci＇s mutual information algorithm was analyzed based on respectively partitioning the plane，constructed by a pair of Lorenz series with the same size，into four or sixteen grids with equal distribution probability in elements on each axis.The improved mutual information algorithms was then promoted based on the original probability matrix that shows the distribution of points corresponding to data pairs of Lorenz series on the plane via the process of sorting the two series，replacing each numerical value by its order number in its own series so as to judge in which data set the element is located and revising the last column and row of the matrix.Finally，after reconstructing the phase space with the optimal time delay，the comparison between the maximal Lyapunov exponent calculated by Rosenstein＇s algorithm from time series and that gained by Jaccobi matrix from Lorenz equation was used to confirm the validity of the new mutual information algorithm.The results show that Cellucci＇s mutual information algorithm may lead to wrong optimal time delay when the series size is not a multiple of elements.The new algorithm，whose result is steadier when large numbers of data pairs are used，can not only eliminate the default of Cellucci＇s algorithm but also is faster than Fraser＇s algorithm.Besides，the lesser difference between the maximal Lyapunov exponents calculated by the two algorithms shows that the new mutual information algorithm is available and feasible.

phase space reconstruction;time delay;mutual information;maximal Lyapunov exponent

023

A

10.13465/j.cnki.jvs.2015.02.014

2013－06－02修改稿收到日期:2013－09－10

蒋爱华男，博士生，1980年11月生

华红星男，教授，博士生导师，1954年生邮箱:hhx@sjtu.edu.cn