相空間重構(gòu)延遲時間互信息改進算法研究

蔣愛華，周璞，章藝，華宏星

(1.中國船舶重工集團公司第704研究所減振中心，上海 200031;2.上海交通大學(xué)機械系統(tǒng)振動國家重點實驗室，上海 200240)

相空間重構(gòu)延遲時間互信息改進算法研究

蔣愛華1，2，周璞1，章藝1，華宏星2

(1.中國船舶重工集團公司第704研究所減振中心，上海 200031;2.上海交通大學(xué)機械系統(tǒng)振動國家重點實驗室，上海 200240)

針對改進互信息算法利于快速可靠獲得時間序列相空間重構(gòu)的延遲時間問題，通過等邊緣分布2、4等分Lorenz時間序列構(gòu)成平面分析Cellucci互信息算法缺陷;用大小順序值代替原序列數(shù)值、判斷新序列數(shù)值所在等邊緣概率區(qū)間獲得概率分布矩陣、修正概率分布矩陣最末行與列改進Cellucci互信息算法;以改進算法所得最佳延遲時間進行Lorenz時間序列相空間重構(gòu)并以小數(shù)據(jù)量法得出其最大Lyapunov指數(shù)，對比雅可比矩陣法所得最大Lyapunov指數(shù)以確認(rèn)改進算法的有效性。結(jié)果表明，時間序列長度不能整除劃分區(qū)間數(shù)時Cellucci互信息算法會獲得錯誤的最佳延遲時間;所提改進算法能消除Cellucci算法缺陷，且計算速度快于Fraser算法;數(shù)據(jù)序列長度較大時改進算法結(jié)果更穩(wěn)定;由兩種最大Lyapunov指數(shù)計算方法所得結(jié)果間誤差較小，表明改進的互信息算法有效、可靠。

相空間重構(gòu);延遲時間;互信息;Cellucci算法;最大Lyapunov指數(shù)

對測試所得振動信號時間序列進行相空間重構(gòu)，可在不清楚振動系統(tǒng)結(jié)構(gòu)與影響參數(shù)時研究振動系統(tǒng)特性。因此在非線性振動系統(tǒng)分析中獲得一定應(yīng)用［1－2］。

互信息算法為相空間重構(gòu)中確定延遲時間的常用方法。設(shè){s(ti)}(i=1，2…N)為試驗觀測所得時間間隔為△t的一維性時間序列，該序列可重構(gòu)為

式中:τ為延遲時間。

式中:H(Q)，H(S)，H(Q，S)分別為Q、S的信息熵及Q與S間互信息熵，通常以第一個互信息極小值點所在位置作為最佳延遲時間。

在計算兩時間序列間互信息時常將兩時間序列S與Q值對應(yīng)到兩相互垂直軸上，構(gòu)造成SQ平面;并將每個序列中最大與最小值間區(qū)域按一定準(zhǔn)則分別劃分為若干區(qū)間si，qj，則該兩序列數(shù)據(jù)對將對應(yīng)于SQ平面中點。通過計算區(qū)間si中數(shù)據(jù)對點數(shù)，除以總數(shù)據(jù)對數(shù)N則可得各區(qū)間邊緣分布概率。兩時間序列區(qū)間相互重合區(qū)域為(si，qj)，計算該區(qū)域內(nèi)數(shù)據(jù)對點數(shù)并除以總點數(shù)，則得重合區(qū)域內(nèi)聯(lián)合概率密度，從而可得S，Q序列間互信息。H(Q)，H(S)，H(Q，S)表達式為

式中:Pq(qi)，Ps(si)，Psq(si，qj)分別為Q在qi區(qū)域的邊緣分布概率密度、S在si區(qū)域的邊緣分布概率密度及S與Q在(si，qi)區(qū)域的聯(lián)合概率密度。

SQ平面劃分示意圖見圖1。Q，S間互信息化簡為

圖1 SQ平面劃分示意圖Fig.1 Skeleton map of SQ plane partition

1 互信息算法

1.1 SQ平面劃分準(zhǔn)則

劃分SQ平面準(zhǔn)則主要有等邊緣概率分布與等間距兩種(圖1)。等邊緣概率分布劃分區(qū)間即使劃分后區(qū)間具有相同數(shù)據(jù)對點數(shù)。第一個用互信息法確定最佳延遲時間Fraser算法即采用該方法［3］。Fraser算法按等邊緣分布將SQ平面劃分為4個網(wǎng)格，判斷每個網(wǎng)格是否存在子結(jié)構(gòu)或已稀疏。若無子結(jié)構(gòu)或已稀疏則無需細(xì)分，若存在子結(jié)構(gòu)則據(jù)邊緣概率對含子結(jié)構(gòu)區(qū)域進行等邊緣分布劃分，直到各區(qū)域內(nèi)無子結(jié)構(gòu)存在。判斷是否存在子結(jié)構(gòu)準(zhǔn)則為數(shù)據(jù)對是否均勻分布于區(qū)域內(nèi)。Fraser算法示意圖見圖2。由圖2看出，區(qū)域R1(2)、R1(3)中數(shù)據(jù)點均勻分布無需進一步劃分;區(qū)域R1(1)、R1(4)、R2(4，2)中存在子結(jié)構(gòu)，則需進一步劃分。

因此，F(xiàn)raser算法計算互信息時子結(jié)構(gòu)網(wǎng)格數(shù)按4的冪次方增加，且為確定子結(jié)構(gòu)是否存在需將各區(qū)域按等邊緣概率分布劃分兩層，判斷各下層子結(jié)構(gòu)網(wǎng)格中的數(shù)據(jù)點數(shù)。計算過程中需存儲各層子結(jié)構(gòu)概率分布，故Fraser算法耗時、耗空間。

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

等間距劃分區(qū)間可使各劃分區(qū)間距上下限差值相同。與Fraser算法相比該方法計算速度較且程序編譯簡單，但存在互信息對劃分區(qū)間數(shù)極敏感，即劃分區(qū)間太多則會有大量區(qū)域內(nèi)數(shù)據(jù)點數(shù)為0或1，無法反映數(shù)據(jù)對的概率分布，并獲得錯的互信息。反之，若劃分區(qū)間太少則各區(qū)域內(nèi)存在子結(jié)構(gòu)而獲得互信息值不準(zhǔn)確。此外，該方法所得結(jié)果不適用于對不同變量的時間序列需相同等間距劃分區(qū)間數(shù)的多變量非線性時間序列重構(gòu)。

圖3 Rossler微分方程中三變量時間序列等間距劃分下求互信息Fig.3 Mutual information of time series from the three variables ofRossler with equal－distance elements

以Rossler數(shù)值大小分布不均時間序列為例，采用相同等區(qū)間劃分?jǐn)?shù)劃分三變量的時間序列，會得到錯誤的最佳延遲時間。圖3(a)為運用四階龍格－庫塔所得10000點Rossler離散序列三維圖，起始點為(－1，0，1)，時間步長為0.05。圖3(b)、(c)分別為Rossler離散序列中三變量的數(shù)值分布柱狀圖、等間距40等分?jǐn)?shù)值區(qū)間時延遲點數(shù)1～200的互信息值圖。可見由變量Z時間序列所得最佳延遲時間不準(zhǔn)確。

Rossler方程為

式中:d=0.2，e=0.4，f=5.7均為系數(shù)。

大量研究集中于等間距劃分區(qū)間數(shù)確定。對總長度為N的序列，Mosteller等［4］提出劃分個區(qū)間，Bendat等［5］提出劃分1.87(n－1)0.4個區(qū)間，Rissanen［6］由算法的隨機復(fù)雜性通過對劃分區(qū)間進行理論推導(dǎo)認(rèn)為，擁有最小隨機復(fù)雜性的劃分區(qū)間數(shù)即為最佳劃分區(qū)間。隨機復(fù)雜性公式為

式中:F(m)為隨機復(fù)雜性，m為劃分區(qū)間數(shù);R為序列中最大最小值之差值;Ni(i=1，2…)為劃分后各區(qū)間數(shù)據(jù)點數(shù)，即

由于N一般較大，式(5)、(6)中存在大數(shù)階乘會在計算中溢出，故將式(4)后兩項變換為

由于復(fù)雜度本身的計算為耗時耗空間過程，從而可抵消等間距劃分節(jié)約時間與空間優(yōu)點，目前尚無較好方法用于確定最佳等間距劃分區(qū)間數(shù)。Cellucci等［7］提出的基于統(tǒng)計的等邊緣劃分區(qū)間互信息算法，不僅能滿足多變量時間序列重構(gòu)要求，且計算速度更快。

1.2 Cellucci互信息算法

互信息算法［7］建立于兩序列S，Q間統(tǒng)計獨立假設(shè)，此時S，Q形成的數(shù)據(jù)點會在SQ平面均勻分布(圖1)，由式(2)知S，Q間互信息將為零。此時SQ平面內(nèi)任意區(qū)域點數(shù)為已知，設(shè)S軸第i區(qū)間內(nèi)點數(shù)為OS(i)，Q軸第j區(qū)間內(nèi)點數(shù)為OQ(i)，則區(qū)域(i，j)內(nèi)點數(shù)ESQ(i，j)為

采用等邊緣概率劃分區(qū)間，且S，Q用相同劃分區(qū)間數(shù)NE，邊緣分布密度PS(i)，PQ(j)變?yōu)镻S(i)= PQ(j)=1/NE，則式(7)變?yōu)?/p>

此時若已知ESQ(i，j)，則可得NE，取ESQ(i，j)≥5，則式(8)變?yōu)?/p>

取NE為滿足以上要求的最大整數(shù)。當(dāng)N為NE的倍數(shù)時，將S，Q軸分別按等邊緣分布概率劃分為NE個區(qū)間，則式(2)可變?yōu)?/p>

當(dāng)N不是NE的倍數(shù)時則按要求確定NE，即按NE劃分SQ平面后，所有區(qū)域內(nèi)實際點數(shù)滿足OSQ(i，j)≥1;80%區(qū)域內(nèi)滿足OSQ(i，j)≥5。判別是否滿足第二要求所用均勻分布卡方檢驗完成，卡方檢驗公式為

卡方檢驗中所用自由度ν=(NE－1)2，則數(shù)據(jù)對在SQ平面內(nèi)均勻分布概率為

1.3 Cellucci算法缺陷

該算法缺陷為:當(dāng)時間序列S或Q長度不能整除劃分區(qū)間數(shù)NE時其計算的互信息出現(xiàn)錯誤。以Lorenz系統(tǒng)中x變量序列為例，Lorenz系統(tǒng)的微分方程為

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

用4階龍格－庫塔獲得100 000點離散序列，起始點為(8.331，13.291，18.063)，時間步長0.01，其中前4 096個數(shù)據(jù)點三維圖見圖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平面各區(qū)域內(nèi)數(shù)據(jù)點數(shù)用矩陣表示。延時點數(shù)1、10、100、200的兩時間序列等份后各矩陣為

按等邊緣分布對S，Q分別劃分4等份，延時點數(shù)1、10、100、200的兩時間序列等份后各矩陣為

由各矩陣看出，對SQ平面按等邊緣分布2等份時數(shù)據(jù)對點分布矩陣中無零元素，需進一步劃分確認(rèn)是否存在子結(jié)構(gòu)或零元素;對SQ平面按等邊緣分布4等份時矩陣中出現(xiàn)零元素，對SQ平面進一步細(xì)分時對應(yīng)的新矩陣中元素亦為零，與算法［7］中要求所有區(qū)域內(nèi)數(shù)據(jù)對數(shù)大于1不符。即數(shù)據(jù)對數(shù)非劃分區(qū)間整數(shù)倍時，用該算法對SQ平面等邊緣概率劃分將會僅二等份，致最佳延時不準(zhǔn)確。4 096(非NE倍數(shù))與4 500 (NE倍數(shù))個數(shù)據(jù)點按算法［7］計算的Lorenz與Rossler延遲點數(shù)1～200的互信息見圖5。

圖5 不同長度數(shù)據(jù)序列Cellucci法所求互信息Fig.5 Mutual information fromCellucci's algorithm with different series lengths

此外，其它長度的Lorenz、Rossler序列亦存在相同情況，且Cellucci算法中用于確定等邊緣分布劃分區(qū)間數(shù)的另個判據(jù)較復(fù)雜，程序不易實現(xiàn)。

2 互信息改進算法

2.1 最末邊緣區(qū)間分布點數(shù)處理

在Cellucci算法基礎(chǔ)上本文提出改進算法。以滿足式(9)最大整數(shù)NE為等邊緣分布劃分S或Q軸區(qū)間數(shù)，每個區(qū)間內(nèi)數(shù)據(jù)點數(shù)取小于或等于N/NE的最大整數(shù)N1，以N1由小到大依次劃分S或Q軸直至結(jié)束。若N/NE為整數(shù)，則劃分后各邊緣區(qū)間內(nèi)數(shù)據(jù)點數(shù)完全相同;若N/NE不為整數(shù)，則劃分后S，Q軸最末個邊緣區(qū)間內(nèi)分布點數(shù)少于其它區(qū)間。此時按最后一個區(qū)間分布點數(shù)與前各區(qū)間分布點數(shù)的比例進行修正，即設(shè)N2為N整除NE時所得余數(shù)在最后一個邊緣區(qū)間內(nèi)各區(qū)域分布點數(shù)得出后分別乘以N1/N2，得各區(qū)域內(nèi)新的分布點數(shù)，在計算互信息時最后一個區(qū)間點總數(shù)則按N1計算，此時的概率分布矩陣用C表示。

2.2 概率分布矩陣計算

由式(2)知，互信息值大小取決于劃分SQ平面后各區(qū)間的概率分布矩陣。在等邊緣概率劃分區(qū)間前提下分布概率矩陣只與各元素在其序列中大小順序有關(guān)，與序列中各元素具體數(shù)值大小無關(guān)，因此在不改變各數(shù)值大小順序前提下對兩序列中數(shù)值進行任意變換，可簡化SQ平面內(nèi)各區(qū)域的分布概率判斷。判斷各區(qū)域內(nèi)分布點數(shù)時所用方法為:①對S，Q兩序列分別排序，并用S，Q各元素在其序列中的大小順序值代替實際值，分別用兩向量A、B表示，無論S、Q數(shù)值范圍大小，A、B所含數(shù)據(jù)均為1～N的整數(shù)，判斷各區(qū)域內(nèi)分布點數(shù)則可直接用A、B完成;②因已知每個邊緣分布區(qū)間所含點數(shù)為N1，則可依次判斷數(shù)據(jù)對(A(i)，B(i))(i=1，2…N)所在區(qū)域，并將區(qū)域概率分布矩陣C中元素C(m，n)加1，其中m為1+A(i)/N1最大整數(shù)，n為1+B(i)/N1最大整數(shù)。

獲得概率分布矩陣C后用第一種處理方法時式(2)變?yōu)?/p>

3 改進算法驗證

為確認(rèn)算法的有效性需對最佳延遲時間進行驗證，用相空間重構(gòu)的非線性不變量，即利用Lyapunov指數(shù)、關(guān)聯(lián)維數(shù)等判斷重構(gòu)效果［9－10］。本文用最大Lyapunov指數(shù)作為判斷依據(jù)。該指數(shù)計算有兩條途徑，即①已知系統(tǒng)微分方程或映射關(guān)系，利用該指數(shù)的定義或系統(tǒng)微分方程的雅可比矩陣特征值在一段時間內(nèi)的平均值計算獲得Lyapunov指數(shù)譜［11］，其中最大值即為最大Lyapunov指數(shù);②已知系統(tǒng)試驗的時間序列，則重構(gòu)后數(shù)據(jù)可由小數(shù)據(jù)量［12］、Wolf［13］、BBA［14］等方法及改進算法獲得Lyapunov指數(shù)，其中小數(shù)據(jù)量法與Wolf法所得為最大Lyapunov指數(shù)，BBA法可獲得所有Lyapunov指數(shù)譜。

對Lorenz系統(tǒng)，由于微分方程已知，由①可得理論的最大Lyapunov指數(shù)。Lorenz系統(tǒng)雅可比矩陣為

通過4階龍格－庫塔法獲得該系統(tǒng)中不同變量的時間序列，由時間序列可由②獲得另個Lyapunov指數(shù)，所得指數(shù)值與理論指數(shù)值接近程度則可反映重構(gòu)效果的好壞，即所選延遲時間與嵌入維數(shù)的好壞。

圖6 不同數(shù)據(jù)長度的互信息Fig.6 Mutual information by different pointslengthes

分別用Lorenz系統(tǒng)中變量x的前1 024、2 048、4 096、8 192、16 384、24 576、32 768、65 536個數(shù)據(jù)點進行互信息計算。據(jù)改進算法所得1～200延遲點數(shù)時互信息見圖6。計算過程中各數(shù)據(jù)點數(shù)對應(yīng)參數(shù)見表1。由表1看出，用于計算互信息的序列長度較短時計算結(jié)果與穩(wěn)定值有一定偏差，而隨采用數(shù)據(jù)長度增加，用以上方法所得第一最小互信息點即最佳延遲時間均趨于穩(wěn)定;用于計算互信息的序列長度大于等于16 384后最佳延遲點數(shù)均為17。

表1 不同數(shù)據(jù)長度計算互信息時參數(shù)Tab.1 Parameters appearing in the process of mutual information calculation by series with different lengths

以上計算互信息算法均在Matlab6.5上實現(xiàn)，所用計算機CPU頻率3 GHz，內(nèi)存2 GB。由數(shù)據(jù)長度計算時間知，改進方法的計算時間遠(yuǎn)小于Fraser算法，相同運行環(huán)境下用Fraser算法計算長4 096序列互信息費時已超200 s。

表2 兩種方法所得最大Lyapunov指數(shù)Tab.2 Maximal Lyapunov exponents gained from Jacobi matrix and phase space constructed by time series with different time delay

據(jù)Takens定理，采用大于2倍于分形維數(shù)的嵌入維數(shù)構(gòu)造相空間時重構(gòu)的相空間與原系統(tǒng)非線性不變量［15］相同。已知Lorenz序列的分形維數(shù)為2.07［16］，故本文嵌入維數(shù)取5，通過產(chǎn)生的100 000點時間序列對相空間進行重構(gòu)，由重構(gòu)數(shù)據(jù)通過小數(shù)據(jù)量法計算出最大Lyapunov指數(shù)，該計算過程可直接用Tisean 3.0.0完成。用Lorenz系統(tǒng)方程的雅可比矩陣通過QR分解方法獲得理論最大Lyapunov指數(shù)。兩種方法計算值見表2。由表2看出，延遲時間取17時小數(shù)據(jù)量法所得最大Lyapunov指數(shù)與理論值相差0.276%，與最佳延遲時間18相差僅1點。

4 結(jié)論

(1)時間序列長度為等邊緣劃分區(qū)間數(shù)的整數(shù)倍時Cellucci互信息算法有效;反之，時間序列長度不為劃分區(qū)間的整數(shù)倍時Cellucci互信息算法會獲得錯誤的最佳延遲時間。

(2)對兩時間序列進行排序，可用序列中各數(shù)值在序列中大小順序值代替原數(shù)值實現(xiàn)時間序列的數(shù)值轉(zhuǎn)換;判斷序列各數(shù)值所在等邊緣概率區(qū)間方法可得概率分布矩陣，計算速度快，程序易編譯。

(3)數(shù)據(jù)序列長度較大時用改進的最佳延遲時間計算方法結(jié)果更穩(wěn)定;改進互信息算法計算時間遠(yuǎn)小于Fraser互信息算法。

(4)用改進互信息算法的最佳延遲時間所得Lyapunov指數(shù)與理論值的相對誤差較小，表明該算法有效可靠，可消除Cellucci互信息算法缺陷。

(5)本文互信息算法構(gòu)造中均設(shè)每個區(qū)域內(nèi)平均分布5個數(shù)據(jù)點，雖所得最佳延遲時間近似值較準(zhǔn)確，但該假設(shè)是否合理尚需理論推導(dǎo);是否有其它值作為每個區(qū)域內(nèi)平均分布點數(shù)獲得更準(zhǔn)確結(jié)果有待驗證。

［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－h(huán)ua1，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月生

華紅星男，教授，博士生導(dǎo)師，1954年生郵箱:hhx@sjtu.edu.cn