王 基, 楊琪斌,2, 劉樹勇, 位秀雷
(1.海軍工程大學(xué) 動力工程學(xué)院,武漢 430033; 2.國家海洋技術(shù)中心漳州基地籌建辦公室,北京 100018)
信息熵和HQ準(zhǔn)則在最大Lyapunov指數(shù)計算中的應(yīng)用
王 基1, 楊琪斌1,2, 劉樹勇1, 位秀雷1
(1.海軍工程大學(xué) 動力工程學(xué)院,武漢 430033; 2.國家海洋技術(shù)中心漳州基地籌建辦公室,北京 100018)
最大Lyapunov指數(shù)是判斷時間序列是否為混沌的一個重要判據(jù),目前應(yīng)用比較廣泛的是小數(shù)據(jù)量法。將信息熵和HQ準(zhǔn)則應(yīng)用在最大Lyapunov指數(shù)的算法中,改進(jìn)了小數(shù)據(jù)量法。信息熵優(yōu)化了相空間重構(gòu)參數(shù),克服了獨立求解重構(gòu)參數(shù)的不足;利用HQ準(zhǔn)則確定鄰近點個數(shù)增加了計算時的精度。仿真實驗表明該改進(jìn)的小數(shù)據(jù)量法在計算最大Lyapunov時具有良好的準(zhǔn)確性,對噪聲具有良好的魯棒性。
信息熵;HQ準(zhǔn)則;小數(shù)據(jù)量法;Lyapunov指數(shù)
混沌已經(jīng)應(yīng)用于許多領(lǐng)域,想要利用混沌,必須要對系統(tǒng)進(jìn)行混沌識別。一般來說,一個動力學(xué)系統(tǒng)的最大Lyapunov指數(shù)大于零時,系統(tǒng)處于混沌狀態(tài)[1]。求得實測時間序列的Lyapunov指數(shù)對于故障信號實時診斷具有重要的意義。在Lyapunov指數(shù)的計算過程中,存在很多問題,比如計算復(fù)雜、計算精度不夠高、無法得到實測時間序列的動力學(xué)系統(tǒng)的數(shù)學(xué)表達(dá)式等。WOLF等[2]提出的軌道跟蹤法具有開創(chuàng)性的意義。為后續(xù)的各種算法的出現(xiàn)打下了堅實的基礎(chǔ)。但是軌道跟蹤法由于容易受到參數(shù)的影響,計算精確度比較差。ROSENSTEIN等[3]在Wolf方法基礎(chǔ)上,提出了計算最大Lyapunov指數(shù)的小數(shù)據(jù)量法。小數(shù)據(jù)量法改善了計算精度,但是在參數(shù)的選取上仍然存在諸多不足。蔣愛華等[4]使用了改進(jìn)的互信息法計算時間延遲,提高了計算最大Lyapunov指數(shù)的速度。楊愛波等[5]在利用小數(shù)據(jù)量法計算最大Lyapunov指數(shù)時,使用空間柵格法選取最近鄰點,大大提高了計算速度,但是對于噪聲魯棒性不佳。劉樹勇等[6]在鄰近點搜索時應(yīng)用了kd樹算法,提高了鄰近點搜索效率,加快了計算速度。楊永鋒等[7]使用加權(quán)平均計算平均周期,并用最大無波動區(qū)間作為計算最大Lyapunov指數(shù)的擬合區(qū)域,具有便捷性,易于實現(xiàn)。李彬彬[8]將若干最優(yōu)的時間延遲點對應(yīng)的最大Lyapunov值求均值,過程簡單,但是誤差較大。在所有的算法當(dāng)中,小數(shù)據(jù)量法由于可以實現(xiàn)對于不完全數(shù)據(jù)的計算,使用的最多。但是在確定相空間重構(gòu)的參數(shù)和鄰近點這兩個重要的環(huán)節(jié)上還有許多的不足亟待解決。例如確定重構(gòu)參數(shù)時缺少整體性,鄰近點個數(shù)主要是靠經(jīng)驗主觀確定等。
在有關(guān)算法中,主要采用TAKENS[9]的嵌入定理進(jìn)行相空間重構(gòu),目前對嵌入維數(shù)m和時間延遲τ這兩個參數(shù)的選取主要是把嵌入維數(shù)和時間延遲分別單獨求解,但是這種方法因為不能夠很好的保持原動力系統(tǒng)整體的特性,所以確定的相空間并不一定最佳。本文提出一種新的相空間重構(gòu)方法,此方法采用信息熵模型來確定嵌入維數(shù)m和時間延遲τ,用遺傳算法對建立在高維空間的信息熵模型進(jìn)行求解,從而實現(xiàn)了對重要的重構(gòu)參數(shù)的優(yōu)化。這種方法不僅保證了兩個重構(gòu)參數(shù)的相互聯(lián)系性,擴充了兩個重構(gòu)參數(shù)的整體性關(guān)系,還可以在重構(gòu)之后保持原有的動力學(xué)關(guān)系。鄰近點的個數(shù)的選?。亨徑c數(shù)量太少會導(dǎo)致計算精度差;數(shù)量太多則會使計算變得繁瑣。在計算中通常使用最多的是固定鄰近點個數(shù)法和固定鄰域半徑法,但是它們都存在著明顯的不足,缺乏足夠的說服力。本文利用HQ(Hannan-Quinn)準(zhǔn)則[10]來實現(xiàn)鄰近點個數(shù)的選取,避免了引入質(zhì)量差的鄰近點和偽鄰近點引起的不利影響,有效增加了計算精度。
1.1 m和τ的信息熵模型
設(shè)定兩個變量為X={x1,x2,…,xn}和Y={y1,y2,…,yn},變量的先驗概率為{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è)混沌時間序列為x(1),x(2),…,x(n),…,則一定有合適的嵌入維數(shù)m和時間延遲τ的相空間X(n)=(x(n),x(n+τ),…,x(n+(m-1)τ))∈Rm,(n=1,2,…),使得重構(gòu)相空間與原混沌系統(tǒng)具有等價關(guān)系。即存在一個映射F:Rm→Rm能將原混沌系統(tǒng)復(fù)原出來,相空間點的軌跡表達(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)行化簡后為
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)致無法確定時間序列未來某時刻的值。
由式(7)可得,在選取合適的嵌入維數(shù)m和時間延遲τ下,f可以反映原系統(tǒng)運動模式。首先,得出f中m和τ具有的一般熵關(guān)系;其次,用神經(jīng)網(wǎng)絡(luò)逼近f。
式(7)還說明了,在知道x(n),x(n+τ),…,x(n+(m-1)τ)后能夠確定x(n+mτ),即只有在知道m(xù)個時刻的值x(n),x(n+τ),…,x(n+(m-1)τ)后,才可以徹底消除未來某時刻x(n+mτ)的不定性。故嵌入維數(shù)m和延遲時間τ是必須緊密聯(lián)系才能消除未來值的不確定性。信息熵能夠刻畫不定性,所以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ù)是一個有關(guān)m和τ的表達(dá)式,雖然利用傳統(tǒng)的優(yōu)化算法可以求解式(8),但是由于這個熵函數(shù)非常復(fù)雜,使用傳統(tǒng)算法的可操作性不高,所以采用遺傳算法進(jìn)行求解。
遺傳算法(GA)是一種非數(shù)值優(yōu)化算法,使用它求解優(yōu)化問題只需目標(biāo)函數(shù)就可以進(jìn)行優(yōu)化問題求解,并且沒有傳統(tǒng)的優(yōu)化算法的缺點。
遺傳算法求解的算法如下:
1)編碼:參數(shù)m和τ為非負(fù)整數(shù),采用二進(jìn)制編碼。
2)初始群體的確定:參照實際問題反復(fù)試驗隨機產(chǎn)生群體規(guī)模為Q[30,80]。
4)選擇算子:使用比例選擇算子。即個體在下一代群體中的個數(shù)由該個體的適應(yīng)值在種群總的適應(yīng)值中的比例來決定。
5)交叉算子:使用兩點交叉算子。交叉概率pc為0.7 6)變異算子:使用基本位變異算子。變異概率pm:0.01 7)終止條件:最大迭代次數(shù)T<100。 鄰近點的個數(shù)是混沌特征指數(shù)計算過程中的另一個重要參數(shù),如何確定最優(yōu)個數(shù)值得研究。本文利用HQ準(zhǔn)則[10]來計算得到最優(yōu)鄰近點個數(shù)值。 經(jīng)過研究發(fā)現(xiàn),對于建立擬合模型,需要考慮模型的復(fù)雜度和擬合效果。充足數(shù)量的模型參數(shù)可以保證擬合精度;參數(shù)過多則使復(fù)雜度增大,甚至導(dǎo)致過擬合現(xiàn)象的發(fā)生。 赤池弘次提出赤池信息準(zhǔn)則,簡稱AIC準(zhǔn)則[11-12],用于確定ARMA(p,q)模型的獨立參數(shù)個數(shù): AIC(p,q)=lnσ2+2(p+q+1)/N (10) 式中:σ2是擬合方差,N是擬合數(shù)據(jù)個數(shù)。AIC準(zhǔn)則的意義在于平衡模型的擬合精度和復(fù)雜度。但AIC準(zhǔn)則仍存在著一定的局限性。 HANNAN等對赤池弘次的AIC準(zhǔn)則進(jìn)行了完善和優(yōu)化,然后提出了HQ準(zhǔn)則[10](Hannan-Quinn定階準(zhǔn)則): (11) 式中:D是表征權(quán)重的常數(shù),取D>2;σ2需要根據(jù)不同的需要對具體定義進(jìn)行修改。擬合模型的精度和復(fù)雜度之間的最佳平衡點是式(11)取得最小值時。 先對鄰近點個數(shù)K設(shè)定一個比較大的取值范圍K∈[Kmin,Kmax]。分別對每個K值對應(yīng)的HQ準(zhǔn)則值進(jìn)行計算,公式如下: (12) 得到一系列準(zhǔn)則值后,式(12)的最小值對應(yīng)的K值就是鄰近點個數(shù)的最優(yōu)選擇。 其中,歸一化均方誤差: σ2= (13) 為了檢驗本文中基于信息熵和HQ準(zhǔn)則改進(jìn)的小數(shù)據(jù)量法在計算Lyapunov指數(shù)時的準(zhǔn)確性和對噪聲的魯棒性,分別對Lorenz系統(tǒng)和含噪聲的Henon系統(tǒng)的最大Lyapunov指數(shù)采用不同的方法進(jìn)行計算。 (1)Lorenz系統(tǒng) Lorenz方程 (14) 式中:σ=10,r=28,b=8/3,令x(0)=1,y(0)=0,z(0)=1。采用四階Runge-Kutta法對Lorenz方程進(jìn)行求解,取步長為0.02,得到x的7 500個點的數(shù)據(jù)集,去除暫態(tài)過程的前面5 000個點,最后得變量x的一個2 500個點的時間序列,并將其歸一化到[0,1]區(qū)間,原始數(shù)據(jù)選前2 000個點。 通過第2節(jié)所介紹的方法得到的最優(yōu)嵌入?yún)?shù):m=9,τ=4。 通過第3節(jié)中所介紹方法,得到每一個K值對應(yīng)的HQ準(zhǔ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時,HQ準(zhǔn)則值最小,因此最佳鄰近點個數(shù)可取為11。 圖2中橫坐標(biāo)為演化步數(shù),縱坐標(biāo)為演化距離,虛線斜率為最大Lyapunov指數(shù)的理論值。圖中曲線1,2,3分別為本文方法、文獻(xiàn)[5]方法和小數(shù)量法的演化曲線。采用不同方法計算Lorenz系統(tǒng)的最大Lyapunov指數(shù)和誤差,如表1。從表中可以看出,采用本文算法計算最大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,求得混沌時間序列,去除前面2 000個點,取后邊2 000點。 通過第2節(jié)中所介紹的方法得到的最優(yōu)嵌入?yún)?shù):m=2,τ=1。 通過第3節(jié)中所介紹方法,得到每一個K值對應(yīng)的HQ準(zhǔn)則值,見圖3。 圖3 HQ準(zhǔn)則值與K值的關(guān)系Fig.3 The relationship of HQ values and K values 從圖3看出當(dāng)K=5時,HQ準(zhǔn)則值最小,因此最佳鄰近點個數(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%的噪聲,采用不同的方法計算不同噪聲水平下Henon時間序列的最大Lyapunov指數(shù)和誤差,如表2。仿真結(jié)果表明,計算精度隨著噪聲水平的增大而降低。但是本文方法可以計算含噪聲的Henon時間序列的最大Lyapunov指數(shù),并且計算準(zhǔn)確度比其它方法高,說明該方法對噪聲具有很好的魯棒性。 表2 不同方法計算最大Lyapunov指數(shù)的對比 本文利用信息熵和HQ準(zhǔn)則對小數(shù)據(jù)量法進(jìn)行改進(jìn),用來計算混沌序列最大Lyapunov指數(shù)。利用信息熵來確定相空間重構(gòu)的參數(shù),保證了兩個重構(gòu)參數(shù)之間的相互聯(lián)系,擴充了兩個重構(gòu)參數(shù)的整體性關(guān)系,實現(xiàn)了對重構(gòu)參數(shù)的優(yōu)化求解,還可以在重構(gòu)之后保持原有的動力學(xué)關(guān)系;利用HQ準(zhǔn)則來確定鄰近點的個數(shù),可以消除引入質(zhì)量差的鄰近點和偽鄰近點對計算的不利影響,增加了計算準(zhǔn)確性。仿真實驗的結(jié)果表明,將信息熵和HQ準(zhǔn)則應(yīng)用于小數(shù)據(jù)量法來計算最大Lyapunov指數(shù)是可行的,并且具有良好的準(zhǔn)確性;對混沌系統(tǒng)加入噪聲之后算出的Lyapunov指數(shù)準(zhǔn)確度仍然很高,說明本文方法對噪聲具有良好的魯棒性。 [ 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 ] 蔣愛華,周璞,章藝,等.相空間重構(gòu)延遲時間互信息改進(jìn)算法研究[J].振動與沖擊,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 ] 楊愛波,王基,劉樹勇,等.基于空間柵格法的最大Lyapunov指數(shù)算法研究[J].電子學(xué)報,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 ] 劉樹勇,楊慶超,位秀雷,等.鄰近點快速搜索方法在混沌識別中的應(yīng)用[J].華中科技大學(xué)學(xué)報(自然科學(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ù)量法計算最大Lyapunov指數(shù)的參數(shù)選擇[J].振動、測試與診斷,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 ] 李彬彬.非線性心音時間序列的最大Lyapunov指數(shù)[J]. 上海電機學(xué)院學(xué)報,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) 國家自然科學(xué)基金(51179197);海洋工程國家重點實驗室(上海交通大學(xué))開放課題(1009) 2015-08-17 修改稿收到日期:2016-01-03 王基 男,副教授,1964年6月生 楊琪斌 男,助理工程師,1991年8月生 O322 A 10.13465/j.cnki.jvs.2017.01.0192 鄰近點個數(shù)的選擇
3 仿真實驗
4 結(jié) 論