混沌振動Lyapunov特征指數(shù)重構(gòu)演化矩陣算法

劉樹勇， 位秀雷，方遠(yuǎn)， 楊慶超

(海軍工程大學(xué) 動力工程學(xué)院，湖北 武漢 430033)

?

混沌振動Lyapunov特征指數(shù)重構(gòu)演化矩陣算法

劉樹勇， 位秀雷，方遠(yuǎn)， 楊慶超

(海軍工程大學(xué) 動力工程學(xué)院，湖北 武漢 430033)

針對Lyapunov指數(shù)計算過程中，鄰近點(diǎn)對演化難以有效刻畫混沌軌道在相空間不同方向的動態(tài)特性問題，本文研究了混沌重構(gòu)相空間中映射矩陣的構(gòu)造方法，推導(dǎo)了混沌振動實(shí)測時間序列重構(gòu)相空間中Lyapunov特征指數(shù)演化矩陣算法表達(dá)式，分析了軌道演化過程中切空間級聯(lián)矩陣對參考點(diǎn)鄰域的拉伸折疊作用，揭示其特征向量方向是鄰域變形的方向、奇異值大小對應(yīng)特征方向上鄰域拉伸和壓縮的大小。構(gòu)造了包含映射值域和定義域的增廣矩陣，通過微擾法分析了計算過程的誤差。計算結(jié)果表明，應(yīng)用重構(gòu)演化矩陣算法能夠有效提取混沌時間序列最大Lyapunov指數(shù)。設(shè)計了非線性混沌振動實(shí)驗(yàn)裝置，證實(shí)了該方法的有效性。

演化矩陣；切空間；投影算子；混沌；重構(gòu)相空間；Lyapunov指數(shù)；混沌振動

混沌是機(jī)械、生物、物理、化學(xué)乃至社會、經(jīng)濟(jì)、大氣等復(fù)雜系統(tǒng)中存在的普遍現(xiàn)象[1]。Lyapunov指數(shù)是實(shí)現(xiàn)混沌識別的重要指標(biāo)。國內(nèi)外學(xué)者對Lyapunov指數(shù)的計算進(jìn)行了大量的研究。主要集中系統(tǒng)數(shù)學(xué)模型已知的情況下 Lyapunov指數(shù)譜[2-3]的計算以及直接從時間序列中計算Lyapunov指數(shù)。其中，WOLF等提出的時間序列最大Lyapunov指數(shù)的軌道追蹤法對于計算該指數(shù)具有開創(chuàng)性的意義[4]。在該方法中，通過面積和體積的演化可以計算時間序列的所有Lyapunov指數(shù)譜。但該方法對軌道追蹤有嚴(yán)格的要求，如：當(dāng)兩軌線的距離大于閾值時，需要重新搜索一個最近鄰點(diǎn)，同時距離參考軌線的距離必須最小、與被置換向量方向之間的夾角盡可能小等。隨后，KANTZ等提出了穩(wěn)健的Lyapunov指數(shù)算法，通過追蹤?quán)徲虻难莼^程來進(jìn)行計算指數(shù)值[5]；ROSENSTEIN等根據(jù)實(shí)測數(shù)據(jù)的特點(diǎn)提出了小數(shù)據(jù)量算法，其優(yōu)點(diǎn)是在數(shù)據(jù)量較少時也可以實(shí)現(xiàn)對指數(shù)的計算[6]。隨后，國內(nèi)學(xué)者對小數(shù)據(jù)量最大Lyapunov指數(shù)計算參數(shù)選擇提出了改進(jìn)方法[7-8]，并探討了多變量含噪聲混沌時間序列Lyapunov指數(shù)的計算問題[9]。文獻(xiàn)[10]提出了基于奇異值的Lyapunov指數(shù)計算方法，文獻(xiàn)[11]則研究了基于智能網(wǎng)絡(luò)的Lyapunov指數(shù)計算。但針對映射矩陣的重構(gòu)技術(shù)和演化特征還缺乏深入研究，沒有從本質(zhì)上闡述映射過程中的物理含義。同時，實(shí)現(xiàn)時間序列Lyapunov指數(shù)的快速計算也是非常具有挑戰(zhàn)性課題[12-13]。本文對演化矩陣重構(gòu)方法進(jìn)行了研究，并給出仿真和實(shí)驗(yàn)計算結(jié)果。

1 Lyapunov特征指數(shù)重構(gòu)演化矩陣算法

1.1 相空間重構(gòu)

1.2 演化矩陣的計算

圖1 演化矩陣對鄰域的變形作用Fig.1 The effect of evolution matrix on the neighborhood deformation

其中，φn是φ的第n次迭代。應(yīng)用鏈?zhǔn)揭?guī)則可以沿著系統(tǒng)的軌道估計出每一個步長的切映射，并得到n個切映射的乘積：

(1)

(2)

(3)

對于時間序列相空間重構(gòu)吸引子而言，追蹤重構(gòu)軌道上參考點(diǎn)的鄰近點(diǎn)或鄰域的演化是計算Lyapunov指數(shù)譜的基本思想，如圖2所示。但在計算過程中，切空間映射矩陣的計算和投影算子的重構(gòu)規(guī)律還需要深入研究。

圖2 重構(gòu)軌道切空間的演化Fig.2 The evolution of reconstruction track tangent space

假設(shè)Δxi是在xi處切向量的近似，根據(jù)式(3)可得

(4)

{Δxi(jk)=xjk-xi|k=1，2，…，Nb}

(5)

由這些位移向量形成的矩陣記為鄰域矩陣Bxi∈RNb×m。由于參考點(diǎn)和鄰近點(diǎn)都隨著相空間中的軌線不斷演化，由此可以計算出經(jīng)過一次演化后的位移向量：

{ΔXi+1(jk)=Xjk+1-Xi+1|k=1，2，…，Nb}

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式(7)中的逆問題可以采取如下形式來描述：

(15)

(16)

(17)

根據(jù)上面的分析可知，這種方法只需要采用m′×m′矩陣，而不需要m×m矩陣，使計算量減少。將式(15)改寫成如下形式:

(18)

1.3 誤差分析

(19)

(20)

(21)

(22)

(23)

沿著相空間軌道計算出切映射矩陣后，進(jìn)行

QR分解(QxjRxj)，得到Lyapunov指數(shù)譜的計算

(24)

式中:i=1，2，…，m；N是切映射矩陣演化次數(shù)。

2 Hénon混沌系統(tǒng)算例

以Hénon離散混沌系統(tǒng)為研究對象，其中，參數(shù)a=1.4，b=0.3，取3 000個點(diǎn)進(jìn)行分析。

重構(gòu)相空間過程，其延遲時間選擇為1，得到其重構(gòu)吸引子，如圖3(a)所示。計算得到系統(tǒng)在重構(gòu)相空間中的Lyapunov指數(shù)譜曲線，如圖3(b)所示，在演化到400個步長時，曲線達(dá)到穩(wěn)定。根據(jù)計算，得到不同嵌入維數(shù)時的Lyapunov指數(shù)譜，最大Lyapunov指數(shù)與標(biāo)準(zhǔn)值(0.408)的誤差百分比如表1所示。從表中可知，該離散系統(tǒng)在嵌入維數(shù)合適時，計算誤差很小。在嵌入維數(shù)增加過大時，會導(dǎo)致計算結(jié)果誤差迅速大，需要通過相空間重構(gòu)選擇合適的維數(shù)。

圖3 Hénon系統(tǒng)仿真計算Fig.3 Simulation of Hénon

Table 1 Lyapunov exponents of different embedding dimension

嵌入維數(shù)Lyapunov指數(shù)最大Lyapunov指數(shù)與標(biāo)準(zhǔn)值[4]的誤差/%20.4035,-1.33211.1030.3776,-0.7245,-0.1309-7.4540.4564,-1.1943,0.1008,-2.216311.86

3 Lyapunov指數(shù)計算實(shí)驗(yàn)研究

為了驗(yàn)證本文所提方法的效果，設(shè)計了單端磁吸式雙勢阱混沌振動試驗(yàn)裝置，如圖4所示，并對產(chǎn)生的混沌振動信號進(jìn)行分析。該裝置由如下部分組成：激振器、支座、板簧、磁鐵、質(zhì)量塊、外部固定框架、固定板簧裝置。各部件的具體功能為：1)激振器：提供穩(wěn)定外部激勵，保持整個系統(tǒng)持續(xù)穩(wěn)定振動；2)支座：固定激振器，為了保證裝置在水平方向上不受其他外力影響，需嚴(yán)格保持激振器水平放置；3)板簧：為整個裝置的主體元件，為保證板簧具有良好的恢復(fù)力，剛度不易太硬，并且板簧的寬度不宜太小，這樣可以減少其做水平擺動之外的運(yùn)動；4)磁鐵(方形磁鐵)：提供板簧端部吸引力，固定磁鐵的支座可以左右水平移動，進(jìn)而調(diào)節(jié)磁力作用距離；5)質(zhì)量塊：安裝在板簧自由端，受到磁力作用，它和磁鐵一起形成了裝置的雙勢阱，為保證板簧不產(chǎn)生水平擺動之外的運(yùn)動，其質(zhì)量不宜過大；6)外部固定框架：框架需滿足兩個要求，不能是磁性材料，以免對板簧運(yùn)動形成外部干擾；質(zhì)量不宜過大，否則難以保持整個框架的水平性，板簧的水平擺動會受到外部力的干擾；7)板簧固定裝置：在兩個夾緊螺母和板簧之間設(shè)置一個方形鐵片，將螺母和板簧之間的點(diǎn)接觸擴(kuò)展成面接觸，從而使板簧固定更牢靠。

圖4 單端磁吸式混沌振動裝置Fig.4 Two-well potential magnetic leaf spring chaotic vibration experimental rig

實(shí)驗(yàn)激勵頻率范圍為5～25 Hz，采樣頻率為2 KHz，數(shù)據(jù)采集時長為5 s。調(diào)節(jié)功率放大器增益為0.96 V、激勵頻率為15 Hz時，采集數(shù)據(jù)并選擇數(shù)據(jù)長度為3 000的信號段進(jìn)行分析。由偽鄰點(diǎn)嵌入維算法和第一最小互信息法[5]計算重構(gòu)相空間參數(shù)中嵌入維數(shù)為3，延遲時間為10 s，得到的重構(gòu)吸引子如圖5(a)、(b)所示。

圖5 實(shí)測信號的重構(gòu)吸引子Fig.5 Reconstruction attractor of the measured signal

從圖5中可以看出吸引子不同于極限環(huán)和二維環(huán)面等平凡吸引子。為了計算信號的Lyapunov指數(shù)譜，對軌道上的鄰域進(jìn)行分析，并重構(gòu)演化矩陣。計算得到Lyapunov指數(shù)譜為(0.562 80，-0.012 5， -0.659 4)，如圖6所示。由于最大Lyapunov指數(shù)大于零，因此表明所設(shè)計的系統(tǒng)處于混沌振動狀態(tài)。

圖6 重構(gòu)吸引子上的鄰域演化Fig.6 Neighborhood evolution of reconstruction attractor

4 結(jié)論

1)研究了重構(gòu)相空間中映射矩陣的鄰域構(gòu)造法，推導(dǎo)了混沌振動實(shí)測時間序列Lyapunov特征指數(shù)算法表達(dá)式；

2)分析了切空間級聯(lián)矩陣對參考點(diǎn)鄰域的拉伸折疊作用，其奇異值的大小度量了初始距離的拉伸程度，其拉伸折疊的方向?qū)?yīng)于演化矩陣的特征向量方向，矩陣奇異值漲落是導(dǎo)致鄰域發(fā)生幾何形變的內(nèi)在動力；

3)設(shè)計了混沌振動試驗(yàn)臺，并有效提取了實(shí)測信號的Lyapunov指數(shù)譜。數(shù)值計算結(jié)果和實(shí)測數(shù)據(jù)分析結(jié)果證明了該重構(gòu)演化矩陣算法是切實(shí)可行的。

[1]LORENZ E N. A study of the predictability of a 28-variable atmospheric model[J]. Tellus， 1965， 17: 321-333.

[2]GREENE J M，KIM J S. The calculation of Lyapunov spectra[J]. Physica D， 1987， 24(1-3)： 213-225.

[3]BENETTIN G， GALGANI L， GIORGILLI A， et al. Lyapunov characteristic exponents for smooth dynamical systems; a method for computing all of them. Part 1: theory[J]. Meccanica， 1980， 15(1)： 9-20.

[4]WOLF A，SWIFT J B，SWINNEY H L，et al．Determining Lyapunov exponent from a time series[J]．Physica D，1985，16(3)： 285-317．

[5]KANTZ H， SCHREIBER T．Nonlinear time series analysis [M]．Cambridge：Cambridge Univ．Press， 1997： 65-74．

[6]ROSENTEIN M T，COLLINS J J，CARLO J D．A practical method for calculating largest Lyapunov exponents from small data sets [J]．Physical D， 1993， 65(1-2)：117-134．

[7]沈紅兵，索洪敏．一種計算最大Lyapunov指數(shù)的改進(jìn)小數(shù)據(jù)量法及其實(shí)現(xiàn)[J]．遵義師范學(xué)院學(xué)報， 2009， 11(6)： 73-75．

SHEN Hongbing， SUO Hongmin. A kind of calculating the improved small datum of the largest lyapunov exponent and its realization [J]. Journal of zunyi normal college， 2009， 11(6)： 73-75．

[8]楊永鋒，仵敏娟，高喆，等．小數(shù)據(jù)量法計算最大Lyapunov指數(shù)的參數(shù)選擇[J]．振動、測試與診斷， 2012， 32(3)： 371-374．

YANG Yongfeng ， WU Minjuan ， GAO Zhe， et al. Parameters selection for calculating largest lyapunov exponent from small data sets[J]. Journal of vibration， measurement & diagnosis， 2012， 32(3)： 371-374．

[9]陶慧，李瑩，馬小平．含噪聲多變量混沌時間序列的最大Lyapunov指數(shù)計算[J]．河南理工大學(xué)學(xué)報：自然科學(xué)版， 2014， 33(6)： 770-775．

TAO Hui， LI Ying ， MA Xiaoping. Calculation for the largest Lyaponuv exponent of multivariate chaotic time-series with noise[J]. Journal of Henan Polytechnic University：Natural Science， 2014， 33(6)： 770-775．

[10]張曉丹，李志萍，張麗麗．一類基于奇異值分解的Lyapunov指數(shù)計算方法[J]．北京科技大學(xué)學(xué)報， 2005， 27(3)： 371-374．

ZHANG Xiaodan， LI Zhiping， ZHANG Lili. A method based on singular value decomposition for computation of Lyapunov exponent[J]. Journal of University of Science and Technology Beijing， 2005， 27(3)： 371-374．

[11]李冬梅，王正鷗．基于RBF網(wǎng)絡(luò)的動力系統(tǒng)Lyapunov指數(shù)的計算方法[J]．信息與控制， 2004， 33(5)： 523-526．

LI Dongmei， WANG Zhengou. An algorithm for computing lyapunov exponents of a dynamical system based on RBF neural networks[J]. Information and control， 2004， 33(5)： 523-526．

[12]楊愛波，王基，劉樹勇，等．基于空間柵格法的最大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 electronicasinica， 2012， 40(9)： 1871-1875．

[13]劉樹勇，楊慶超，位秀雷，等．鄰近點(diǎn)快速搜索方法在混沌識別中的應(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 method to chaos identification[J]. J. Huazhong Univ. of Sci.& Tech：Natural Science Edition， 2012， 40(11)： 89-92．

[14]蔣愛華，周璞，章藝，等．相空間重構(gòu)延遲時間互信息改進(jìn)算法研究[J]．振動與沖擊， 2015， 14(2)： 71-74．

JIANG Aihua， HOU Pu， HANG Yi， et al. Improved mutual information algorithm for phase space reconstruction[J]. Journal of vibration and shock， 2015， 14(2)： 71-74．

[15]龍海輝，張佃中．基于等概率符號分析方法計算互信息確定延遲時間[J]．計算物理， 2010， 27(3)： 468-473．

LONG Haihui， ZHANG Dianzhong. Determination of delay time with mutual information by equal-probability symbolic analysis method[J]. Chinese journal of computational physics， 2010， 27(3)： 468-473．

[16]OSELEDEC V I. A multiplicative ergodic theorem Liapunov characteristic numbers for dynamical systems[J]. Trans. moscow math. soc.， 1968， 19(2)： 197-231.

[17]ULRICH P. Identification of true and spurious Lyapunov exponents from time series[J]. International journal of bifurcation and chaos， 1992， 2(1)： 155-165.

本文引用格式：

劉樹勇， 位秀雷，方遠(yuǎn)，等.混沌振動Lyapunov特征指數(shù)重構(gòu)演化矩陣算法[J]. 哈爾濱工程大學(xué)學(xué)報， 2017， 38(5)： 746-751.

LIU Shuyong，WEI Xiulei，F(xiàn)ang Yuan，et al.Study on the Lyapunov exponent algorithm based on the reconstruction matrices evolution[J]. Journal of Harbin Engineering University， 2017， 38(5)： 746-751.

A reconstruction evolution matrices algorithm based on Lyapunov characteristic exponent of chaotic vibration

LIU Shuyong， WEI Xiulei， FANG Yuan， YANG Qingchao

(College of Power Engineering，Naval University of Engineering，Wuhan 430033，China)

The reconstruction map matrix of chaotic phase space is studied since the dynamical characteristics of the chaotic orbit at different directions are discovered by the nearest orbit ineffectively， and the formulation of Lyapunov exponent based on evolution matrices in reconstructed phase space is derived. It explains that in the process of orbit evolution， the stretch and fold of the neighbor is the rotation operation of the cascade matrices and the fluctuation of the singular value of the reconstructed tangent space matrices. The SVD matrix of the reconstructed tangent space plays important role on the system dynamical behavior characteristic. The eigenvector direction is the deformation direction of the neighbor， and the stretch extent corresponds to the eigenvalue. The augmented matrix including range and definition domains is constructed， and the computation error is analyzed with perturbation. Numerical results show that the presented algorithm is useful to extract the maximal Lyapunov exponent of the chaotic time series. The experimental rig is designed and the Lyapunov exponent of the measured signal is computed with the method effectively.

evolution matrix; tangent space; projection operator; chaos; reconstructed phase space； Lyapunov exponent; chaotic vibration

2015-12-24.

日期：2017-04-27.

國家自然科學(xué)基金項(xiàng)目(5157942) ；國家自然科學(xué)基金青年基金項(xiàng)目(51509253)

劉樹勇(1975-)， 男， 副教授; 位秀雷(1988-)， 男， 博士研究生.

位秀雷， E-mail:wxlcln@163.com.

10.11990/jheu.201512081

TN911

A

1006-7043(2017)05-0746-06

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170427.1644.184.html