王 基, 楊琪斌，2, 劉樹勇, 位秀雷

(1.海軍工程大學(xué) 動(dòng)力工程學(xué)院，武漢 430033； 2.國(guó)家海洋技術(shù)中心漳州基地籌建辦公室，北京 100018)

信息熵和HQ準(zhǔn)則在最大Lyapunov指數(shù)計(jì)算中的應(yīng)用

王 基1, 楊琪斌1，2, 劉樹勇1, 位秀雷1

(1.海軍工程大學(xué) 動(dòng)力工程學(xué)院，武漢 430033； 2.國(guó)家海洋技術(shù)中心漳州基地籌建辦公室，北京 100018)

最大Lyapunov指數(shù)是判斷時(shí)間序列是否為混沌的一個(gè)重要判據(jù)，目前應(yīng)用比較廣泛的是小數(shù)據(jù)量法。將信息熵和HQ準(zhǔn)則應(yīng)用在最大Lyapunov指數(shù)的算法中，改進(jìn)了小數(shù)據(jù)量法。信息熵優(yōu)化了相空間重構(gòu)參數(shù)，克服了獨(dú)立求解重構(gòu)參數(shù)的不足；利用HQ準(zhǔn)則確定鄰近點(diǎn)個(gè)數(shù)增加了計(jì)算時(shí)的精度。仿真實(shí)驗(yàn)表明該改進(jìn)的小數(shù)據(jù)量法在計(jì)算最大Lyapunov時(shí)具有良好的準(zhǔn)確性，對(duì)噪聲具有良好的魯棒性。

信息熵；HQ準(zhǔn)則；小數(shù)據(jù)量法；Lyapunov指數(shù)

混沌已經(jīng)應(yīng)用于許多領(lǐng)域，想要利用混沌，必須要對(duì)系統(tǒng)進(jìn)行混沌識(shí)別。一般來(lái)說(shuō)，一個(gè)動(dòng)力學(xué)系統(tǒng)的最大Lyapunov指數(shù)大于零時(shí)，系統(tǒng)處于混沌狀態(tài)[1]。求得實(shí)測(cè)時(shí)間序列的Lyapunov指數(shù)對(duì)于故障信號(hào)實(shí)時(shí)診斷具有重要的意義。在Lyapunov指數(shù)的計(jì)算過(guò)程中，存在很多問(wèn)題，比如計(jì)算復(fù)雜、計(jì)算精度不夠高、無(wú)法得到實(shí)測(cè)時(shí)間序列的動(dòng)力學(xué)系統(tǒng)的數(shù)學(xué)表達(dá)式等。WOLF等[2]提出的軌道跟蹤法具有開創(chuàng)性的意義。為后續(xù)的各種算法的出現(xiàn)打下了堅(jiān)實(shí)的基礎(chǔ)。但是軌道跟蹤法由于容易受到參數(shù)的影響，計(jì)算精確度比較差。ROSENSTEIN等[3]在Wolf方法基礎(chǔ)上，提出了計(jì)算最大Lyapunov指數(shù)的小數(shù)據(jù)量法。小數(shù)據(jù)量法改善了計(jì)算精度，但是在參數(shù)的選取上仍然存在諸多不足。蔣愛(ài)華等[4]使用了改進(jìn)的互信息法計(jì)算時(shí)間延遲，提高了計(jì)算最大Lyapunov指數(shù)的速度。楊愛(ài)波等[5]在利用小數(shù)據(jù)量法計(jì)算最大Lyapunov指數(shù)時(shí)，使用空間柵格法選取最近鄰點(diǎn)，大大提高了計(jì)算速度，但是對(duì)于噪聲魯棒性不佳。劉樹勇等[6]在鄰近點(diǎn)搜索時(shí)應(yīng)用了kd樹算法，提高了鄰近點(diǎn)搜索效率，加快了計(jì)算速度。楊永鋒等[7]使用加權(quán)平均計(jì)算平均周期，并用最大無(wú)波動(dòng)區(qū)間作為計(jì)算最大Lyapunov指數(shù)的擬合區(qū)域，具有便捷性，易于實(shí)現(xiàn)。李彬彬[8]將若干最優(yōu)的時(shí)間延遲點(diǎn)對(duì)應(yīng)的最大Lyapunov值求均值，過(guò)程簡(jiǎn)單，但是誤差較大。在所有的算法當(dāng)中，小數(shù)據(jù)量法由于可以實(shí)現(xiàn)對(duì)于不完全數(shù)據(jù)的計(jì)算，使用的最多。但是在確定相空間重構(gòu)的參數(shù)和鄰近點(diǎn)這兩個(gè)重要的環(huán)節(jié)上還有許多的不足亟待解決。例如確定重構(gòu)參數(shù)時(shí)缺少整體性，鄰近點(diǎn)個(gè)數(shù)主要是靠經(jīng)驗(yàn)主觀確定等。

在有關(guān)算法中，主要采用TAKENS[9]的嵌入定理進(jìn)行相空間重構(gòu)，目前對(duì)嵌入維數(shù)m和時(shí)間延遲τ這兩個(gè)參數(shù)的選取主要是把嵌入維數(shù)和時(shí)間延遲分別單獨(dú)求解，但是這種方法因?yàn)椴荒軌蚝芎玫谋３衷瓌?dòng)力系統(tǒng)整體的特性，所以確定的相空間并不一定最佳。本文提出一種新的相空間重構(gòu)方法，此方法采用信息熵模型來(lái)確定嵌入維數(shù)m和時(shí)間延遲τ，用遺傳算法對(duì)建立在高維空間的信息熵模型進(jìn)行求解，從而實(shí)現(xiàn)了對(duì)重要的重構(gòu)參數(shù)的優(yōu)化。這種方法不僅保證了兩個(gè)重構(gòu)參數(shù)的相互聯(lián)系性，擴(kuò)充了兩個(gè)重構(gòu)參數(shù)的整體性關(guān)系，還可以在重構(gòu)之后保持原有的動(dòng)力學(xué)關(guān)系。鄰近點(diǎn)的個(gè)數(shù)的選?。亨徑c(diǎn)數(shù)量太少會(huì)導(dǎo)致計(jì)算精度差；數(shù)量太多則會(huì)使計(jì)算變得繁瑣。在計(jì)算中通常使用最多的是固定鄰近點(diǎn)個(gè)數(shù)法和固定鄰域半徑法，但是它們都存在著明顯的不足，缺乏足夠的說(shuō)服力。本文利用HQ(Hannan-Quinn)準(zhǔn)則[10]來(lái)實(shí)現(xiàn)鄰近點(diǎn)個(gè)數(shù)的選取，避免了引入質(zhì)量差的鄰近點(diǎn)和偽鄰近點(diǎn)引起的不利影響，有效增加了計(jì)算精度。

1 基于信息熵的相空間重構(gòu)

1.1 m和τ的信息熵模型

設(shè)定兩個(gè)變量為X={x1,x2,…,xn}和Y={y1,y2,…,yn}，變量的先驗(yàn)概率為{p(xi)}i=1,2,…,n和{p(yi)}i=1,2,…,k，可將信息熵定義：

(1)

此定義描述了變量X的不定性。

類似地，聯(lián)合熵定義為：

(2)

式中:p(xi,yi)是聯(lián)合概率。

(3)

(4)

1.2 相空間重構(gòu)的參數(shù)模型建立

設(shè)混沌時(shí)間序列為x(1),x(2),…,x(n),…，則一定有合適的嵌入維數(shù)m和時(shí)間延遲τ的相空間X(n)=(x(n),x(n+τ),…,x(n+(m-1)τ))∈Rm,(n=1,2,…)，使得重構(gòu)相空間與原混沌系統(tǒng)具有等價(jià)關(guān)系。即存在一個(gè)映射F:Rm→Rm能將原混沌系統(tǒng)復(fù)原出來(lái)，相空間點(diǎn)的軌跡表達(dá)式

X(n+τ)=F(X(n)),n=1,2,…

(5)

其中X(n+τ)=(x(n+τ),x(n+2τ),…,x(n+mτ))

式(5)的分量形式為

x(n+jτ)=fj(x(n),x(n+τ),…,x(n+(j-1)τ))，
j=1,2,…,m，n=1,2,…

(6)

將式(6)進(jìn)行化簡(jiǎn)后為

x(n+mτ)=f(x(n),x(n+τ),…,
x(n+(m-1)τ))，n=1,2,…

(7)

混沌系統(tǒng)的復(fù)雜性導(dǎo)致很難直接得到f的解析式；混沌系統(tǒng)高度的非線性則導(dǎo)致無(wú)法確定時(shí)間序列未來(lái)某時(shí)刻的值。

由式(7)可得，在選取合適的嵌入維數(shù)m和時(shí)間延遲τ下，f可以反映原系統(tǒng)運(yùn)動(dòng)模式。首先，得出f中m和τ具有的一般熵關(guān)系；其次，用神經(jīng)網(wǎng)絡(luò)逼近f。

式(7)還說(shuō)明了，在知道x(n),x(n+τ),…,x(n+(m-1)τ)后能夠確定x(n+mτ)，即只有在知道m(xù)個(gè)時(shí)刻的值x(n),x(n+τ),…,x(n+(m-1)τ)后，才可以徹底消除未來(lái)某時(shí)刻x(n+mτ)的不定性。故嵌入維數(shù)m和延遲時(shí)間τ是必須緊密聯(lián)系才能消除未來(lái)值的不確定性。信息熵能夠刻畫不定性，所以m和τ的熵關(guān)系可以建立。記

由以上討論，得到求m和τ的優(yōu)化模型：

目標(biāo)函數(shù)

(8)

約束條件：m和τ為非負(fù)整數(shù)。

根據(jù)式(4)，條件熵變換成聯(lián)合熵為：

minH(m,τ)=
H(X1,X2,…,Xm,Xm+1)-H(X1,…,Xm)

(9)

1.3 求解相空間重構(gòu)參數(shù)的信息熵模型

式(8)的目標(biāo)函數(shù)中的熵函數(shù)是一個(gè)有關(guān)m和τ的表達(dá)式，雖然利用傳統(tǒng)的優(yōu)化算法可以求解式(8)，但是由于這個(gè)熵函數(shù)非常復(fù)雜，使用傳統(tǒng)算法的可操作性不高，所以采用遺傳算法進(jìn)行求解。

遺傳算法(GA)是一種非數(shù)值優(yōu)化算法，使用它求解優(yōu)化問(wèn)題只需目標(biāo)函數(shù)就可以進(jìn)行優(yōu)化問(wèn)題求解，并且沒(méi)有傳統(tǒng)的優(yōu)化算法的缺點(diǎn)。

遺傳算法求解的算法如下：

1)編碼：參數(shù)m和τ為非負(fù)整數(shù)，采用二進(jìn)制編碼。

2)初始群體的確定：參照實(shí)際問(wèn)題反復(fù)試驗(yàn)隨機(jī)產(chǎn)生群體規(guī)模為Q[30，80]。

4)選擇算子：使用比例選擇算子。即個(gè)體在下一代群體中的個(gè)數(shù)由該個(gè)體的適應(yīng)值在種群總的適應(yīng)值中的比例來(lái)決定。

5)交叉算子：使用兩點(diǎn)交叉算子。交叉概率pc為0.7

6)變異算子：使用基本位變異算子。變異概率pm：0.01

7)終止條件：最大迭代次數(shù)T<100。

2 鄰近點(diǎn)個(gè)數(shù)的選擇

鄰近點(diǎn)的個(gè)數(shù)是混沌特征指數(shù)計(jì)算過(guò)程中的另一個(gè)重要參數(shù)，如何確定最優(yōu)個(gè)數(shù)值得研究。本文利用HQ準(zhǔn)則[10]來(lái)計(jì)算得到最優(yōu)鄰近點(diǎn)個(gè)數(shù)值。

經(jīng)過(guò)研究發(fā)現(xiàn)，對(duì)于建立擬合模型，需要考慮模型的復(fù)雜度和擬合效果。充足數(shù)量的模型參數(shù)可以保證擬合精度；參數(shù)過(guò)多則使復(fù)雜度增大，甚至導(dǎo)致過(guò)擬合現(xiàn)象的發(fā)生。

赤池弘次提出赤池信息準(zhǔn)則，簡(jiǎn)稱AIC準(zhǔn)則[11-12]，用于確定ARMA(p，q)模型的獨(dú)立參數(shù)個(gè)數(shù)：

AIC(p,q)=lnσ2+2(p+q+1)/N

(10)

式中:σ2是擬合方差，N是擬合數(shù)據(jù)個(gè)數(shù)。AIC準(zhǔn)則的意義在于平衡模型的擬合精度和復(fù)雜度。但AIC準(zhǔn)則仍存在著一定的局限性。

HANNAN等對(duì)赤池弘次的AIC準(zhǔn)則進(jìn)行了完善和優(yōu)化，然后提出了HQ準(zhǔn)則[10](Hannan-Quinn定階準(zhǔn)則)：

(11)

式中:D是表征權(quán)重的常數(shù)，取D>2；σ2需要根據(jù)不同的需要對(duì)具體定義進(jìn)行修改。擬合模型的精度和復(fù)雜度之間的最佳平衡點(diǎn)是式(11)取得最小值時(shí)。

先對(duì)鄰近點(diǎn)個(gè)數(shù)K設(shè)定一個(gè)比較大的取值范圍K∈[Kmin,Kmax]。分別對(duì)每個(gè)K值對(duì)應(yīng)的HQ準(zhǔn)則值進(jìn)行計(jì)算，公式如下：

(12)

得到一系列準(zhǔn)則值后，式(12)的最小值對(duì)應(yīng)的K值就是鄰近點(diǎn)個(gè)數(shù)的最優(yōu)選擇。

其中，歸一化均方誤差：

σ2=

(13)

3 仿真實(shí)驗(yàn)

為了檢驗(yàn)本文中基于信息熵和HQ準(zhǔn)則改進(jìn)的小數(shù)據(jù)量法在計(jì)算Lyapunov指數(shù)時(shí)的準(zhǔn)確性和對(duì)噪聲的魯棒性，分別對(duì)Lorenz系統(tǒng)和含噪聲的Henon系統(tǒng)的最大Lyapunov指數(shù)采用不同的方法進(jìn)行計(jì)算。

(1)Lorenz系統(tǒng)

Lorenz方程

(14)

式中:σ=10，r=28，b=8/3，令x(0)=1，y(0)=0，z(0)=1。采用四階Runge-Kutta法對(duì)Lorenz方程進(jìn)行求解，取步長(zhǎng)為0.02，得到x的7 500個(gè)點(diǎn)的數(shù)據(jù)集，去除暫態(tài)過(guò)程的前面5 000個(gè)點(diǎn)，最后得變量x的一個(gè)2 500個(gè)點(diǎn)的時(shí)間序列，并將其歸一化到[0，1]區(qū)間，原始數(shù)據(jù)選前2 000個(gè)點(diǎn)。

通過(guò)第2節(jié)所介紹的方法得到的最優(yōu)嵌入?yún)?shù):m=9，τ=4。

通過(guò)第3節(jié)中所介紹方法，得到每一個(gè)K值對(duì)應(yīng)的HQ準(zhǔn)則值，見(jiàn)圖1。

圖1 HQ準(zhǔn)則值與K值的關(guān)系Fig.1the relationship of HQ values and K values based on HQ rules

從圖1看出當(dāng)K=11時(shí)，HQ準(zhǔn)則值最小，因此最佳鄰近點(diǎn)個(gè)數(shù)可取為11。

圖2中橫坐標(biāo)為演化步數(shù)，縱坐標(biāo)為演化距離，虛線斜率為最大Lyapunov指數(shù)的理論值。圖中曲線1,2,3分別為本文方法、文獻(xiàn)[5]方法和小數(shù)量法的演化曲線。采用不同方法計(jì)算Lorenz系統(tǒng)的最大Lyapunov指數(shù)和誤差，如表1。從表中可以看出，采用本文算法計(jì)算最大Lyapunov指數(shù)，準(zhǔn)確性最高。

圖2 演化曲線圖Fig.2 Evolution curve

方法LE誤差小數(shù)據(jù)量法2.89776.8文獻(xiàn)[5]方法2.90326.6本文方法3.03732.3Standard[13]3.1096

(2)含噪聲的Henon系統(tǒng)

Henon映射：

(15)

式中:a=1.4，b=0.3。令y1=y2=0.5，求得混沌時(shí)間序列，去除前面2 000個(gè)點(diǎn)，取后邊2 000點(diǎn)。

通過(guò)第2節(jié)中所介紹的方法得到的最優(yōu)嵌入?yún)?shù)：m=2，τ=1。

通過(guò)第3節(jié)中所介紹方法，得到每一個(gè)K值對(duì)應(yīng)的HQ準(zhǔn)則值，見(jiàn)圖3。

圖3 HQ準(zhǔn)則值與K值的關(guān)系Fig.3 The relationship of HQ values and K values

從圖3看出當(dāng)K=5時(shí)，HQ準(zhǔn)則值最小，因此最佳鄰近點(diǎn)個(gè)數(shù)可取為5。

圖4是10%噪聲水平下演化圖的局部放大圖。橫坐標(biāo)為演化步數(shù), 縱坐標(biāo)為演化距離，虛線斜率為最大Lyapunov指數(shù)的理論值。圖中曲線1,2,3分別為本文方法、文獻(xiàn)[6]方法和小數(shù)據(jù)量法的演化曲線。

圖4 演化曲線圖Fig.4 Evolution curve

分別向Henon系統(tǒng)中加入10%、30%和50%的噪聲，采用不同的方法計(jì)算不同噪聲水平下Henon時(shí)間序列的最大Lyapunov指數(shù)和誤差，如表2。仿真結(jié)果表明，計(jì)算精度隨著噪聲水平的增大而降低。但是本文方法可以計(jì)算含噪聲的Henon時(shí)間序列的最大Lyapunov指數(shù)，并且計(jì)算準(zhǔn)確度比其它方法高，說(shuō)明該方法對(duì)噪聲具有很好的魯棒性。

表2 不同方法計(jì)算最大Lyapunov指數(shù)的對(duì)比

4 結(jié) 論

本文利用信息熵和HQ準(zhǔn)則對(duì)小數(shù)據(jù)量法進(jìn)行改進(jìn)，用來(lái)計(jì)算混沌序列最大Lyapunov指數(shù)。利用信息熵來(lái)確定相空間重構(gòu)的參數(shù)，保證了兩個(gè)重構(gòu)參數(shù)之間的相互聯(lián)系，擴(kuò)充了兩個(gè)重構(gòu)參數(shù)的整體性關(guān)系，實(shí)現(xiàn)了對(duì)重構(gòu)參數(shù)的優(yōu)化求解，還可以在重構(gòu)之后保持原有的動(dòng)力學(xué)關(guān)系；利用HQ準(zhǔn)則來(lái)確定鄰近點(diǎn)的個(gè)數(shù)，可以消除引入質(zhì)量差的鄰近點(diǎn)和偽鄰近點(diǎn)對(duì)計(jì)算的不利影響，增加了計(jì)算準(zhǔn)確性。仿真實(shí)驗(yàn)的結(jié)果表明，將信息熵和HQ準(zhǔn)則應(yīng)用于小數(shù)據(jù)量法來(lái)計(jì)算最大Lyapunov指數(shù)是可行的，并且具有良好的準(zhǔn)確性；對(duì)混沌系統(tǒng)加入噪聲之后算出的Lyapunov指數(shù)準(zhǔn)確度仍然很高，說(shuō)明本文方法對(duì)噪聲具有良好的魯棒性。

[ 1 ] GENCAY R，DECHERT DAVIS W. An algorithm for the n-dimensional unknown dynamical system [J]. Physica D，1992，59：142-157．

[ 2 ] WOLF A，SWIFT J B，SWINNEY H L，et al. Determining Lyapunov exponents from a time series [J]. Physica D，1985，16：285-317.

[ 3 ] ROSENSTEIN M T ，COLLINS J J， DE LUCA C J. A practical method for calculating largest Lyapunov exponents from small data sets[J]. Physica D，1993，65：117-134.

[ 4 ] 蔣愛(ài)華，周璞，章藝,等.相空間重構(gòu)延遲時(shí)間互信息改進(jìn)算法研究[J].振動(dòng)與沖擊，2015，14(2)：71-74. JIANG Aihua，ZHOU Pu，ZHANG Yi. Improved mutual information algorithm for phase space reconstruction[J]. Journal of Vibration and Shock，2015，14(2)：71-74.

[ 5 ] 楊愛(ài)波，王基，劉樹勇,等.基于空間柵格法的最大Lyapunov指數(shù)算法研究[J].電子學(xué)報(bào)，2012，40(9)：1871-1875. YANG Aibo，WANG Ji，LIU Shuyong,et al. An algorithm for computing the largest Lyapunov exponent based on space grid method[J]. Acta Electronica Sinica，2012，40(9) ：1871-1875.

[ 6 ] 劉樹勇，楊慶超，位秀雷,等.鄰近點(diǎn)快速搜索方法在混沌識(shí)別中的應(yīng)用[J].華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，2012,40(11)：89-92. LIU Shuyong，YANG Qingchao，WEI Xiulei,et al.The application of fast searching nearest points methodto chaos identification[J].J.Huazhong Univ. of Sci.&Tech.(Natural Science Edition)，2012,40(11)：89-92.

[ 7 ] 楊永鋒，仵敏娟，高喆,等.小數(shù)據(jù)量法計(jì)算最大Lyapunov指數(shù)的參數(shù)選擇[J].振動(dòng)、測(cè)試與診斷，2012，32(3)：371-374. YANG Yongfeng, WU Minjuan, GAO Zhe,et al. Parameter selection of maximum Lyapunov exponent for small data volume method[J]. Journal of Vibration,Measurement & Diagnosis，2012,32(3)：371-374.

[ 8 ] 李彬彬.非線性心音時(shí)間序列的最大Lyapunov指數(shù)[J]. 上海電機(jī)學(xué)院學(xué)報(bào)，2011，14(1)：17-20. LI Binbin. Largest lyapunov exponents of nonlinear heartbeat time series[J]. Journal of Shanghai Dianji University，2011，14(1)：17-20.

[ 9 ] TAKENS F. Dynamical systems and turbulence[M].Berlin：SpringVerlag,1981,366.

[10] HANNAN E J，QUINN B G.The determination of the order of an autoregression [J].Journal of the Royal Statistical Society SeriesB(Methodological)，1979:190-195.

[11] AKAIKE H. Autoregressive model fitting for control [J]. Annals of the Institute of Statistical Mathematics，1971，23(1)： 163-180.

[12] AKAIKE H. A new look at the statistical model identification [J]. Automatic Control，IEEE Transactions on，1974，19(6)：716-723.

[13] GAO Jianbo，ZHENG Zheming. Local exponential divergence plot and optimal embedding of a chaotic time series[J]. Physics Letters A,1993(181)：153-158.

Application of information entropy and HQ rule in estimating largest Lyapunov exponent

WANG Ji1, YANG Qibin1,2, LIU Shuyong1, WEI Xiulei1

(1. College of Power Engineering，Naval Uniwersity of Engineering，Wuhan 430033，China;2. The Zhangzhou Base Preparation Office of National Ocean Technology Center, Beijing 100018, China)

The largest Lyapunov exponent is an essential criterion to judge if a time series is chaos or not. The small-data method is widely used in chaotic characteristic extraction at present. Here, the information entropy and HQ rule were applied in estimating the largest Lyapunov exponent to improve the small-data method. The information entropy was applied to optimize parameters of phase space reconstruction, and disadvantages of traditional algorithms were overcome clearly. The computational accuracy of LE was improved greatly by using the HQ rule to calculate the number of neighbouring points. Simulation results showed that the improved small-data method here has good performances in estimating the largest Lyapunov exponent, and the algorithm is robust to noise.

information entropy；HQ rule；small-data method；lyapunov exponent (LE)

國(guó)家自然科學(xué)基金(51179197);海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(上海交通大學(xué))開放課題(1009)

2015-08-17 修改稿收到日期：2016-01-03

王基 男，副教授，1964年6月生

楊琪斌 男，助理工程師，1991年8月生

O322

A

10.13465/j.cnki.jvs.2017.01.019