魏春雨 楊威 王梓卉敏 韓清鵬

(遼寧科技大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，鞍山 114051)

基于非線性混沌理論的非平穩(wěn)信號(hào)的比較分析*

魏春雨 楊威 王梓卉敏 韓清鵬?

(遼寧科技大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，鞍山 114051)

本文研究采用基于非線性混沌理論的兩種非線性參數(shù)估計(jì)方法(代替數(shù)據(jù)法和Lyapunov指數(shù)估計(jì)法)對(duì)非平穩(wěn)信號(hào)進(jìn)行分析．首先對(duì)上述兩種非線性方法的具體算法進(jìn)行介紹，然后對(duì)兩組本質(zhì)不同的非平穩(wěn)振動(dòng)信號(hào)進(jìn)行對(duì)比分析．這兩組信號(hào)是通過(guò)測(cè)試具有不同非線性約束邊界條件的薄壁構(gòu)件獲得．分析結(jié)果表明，在時(shí)域波形上直觀相似的非平穩(wěn)信號(hào)，用上述非線性混沌分析的方法可以有效地加以定量區(qū)分．

非線性混沌理論， 非平穩(wěn)信號(hào)， 代替數(shù)據(jù)法， Lyapunov指數(shù)

引言

人們通?？紤]振動(dòng)信號(hào)是否具有平穩(wěn)或非平穩(wěn)的性質(zhì)，而線性或非線性特性則是對(duì)其源系統(tǒng)而言［1］．在工程中，大多數(shù)系統(tǒng)所測(cè)得的信號(hào)是非平穩(wěn)的，而系統(tǒng)是非線性的．區(qū)分非平穩(wěn)信號(hào)在現(xiàn)代信號(hào)處理中，非平穩(wěn)信號(hào)的處理是最為活躍、發(fā)展最為迅速的方向之一［2］，在通訊、雷達(dá)、自動(dòng)控制、模式識(shí)別、機(jī)械振動(dòng)和生物醫(yī)學(xué)工程等領(lǐng)域廣泛應(yīng)用．通過(guò)分析信號(hào)的時(shí)變特征，構(gòu)造合適的時(shí)頻分布并進(jìn)行恰當(dāng)?shù)奶幚?，達(dá)到不同的信號(hào)處理目的．近年來(lái)，人們采用多種時(shí)頻分析方法對(duì)非線性和非平穩(wěn)信號(hào)進(jìn)行分析，包括小波分析、Wagner－Ville分布等［3－6］．此外，還有許多非線性分析的方法也加以應(yīng)用［7－9］．但是，從非線性角度，基于混沌理論的一些方法，如Lyapunov指數(shù)和分?jǐn)?shù)維，對(duì)于具有非線性特性的非平穩(wěn)信號(hào)分析更為有效［10，11］．另外，對(duì)信號(hào)進(jìn)行混沌辨識(shí)，并將非線性信號(hào)從隨機(jī)噪聲中區(qū)分出來(lái)往往是比較困難的．這是因?yàn)榛煦缧盘?hào)和隨機(jī)噪聲信號(hào)通常都具有相似的寬頻特征．目前區(qū)分混沌和噪聲的常用方法主要有兩類:維數(shù)估計(jì)方法和非線性預(yù)測(cè)法．第一類方法是基于隨機(jī)只能存在無(wú)限維吸引子的原理，而有限維吸引子意味著混沌，如在 Grassberger－Procaccia算法中所指出的那樣［10］．第二類方法則是根據(jù)時(shí)間序列的非線性預(yù)測(cè)原理、混沌具有不同于隨機(jī)數(shù)據(jù)的短期預(yù)測(cè)性的性質(zhì)來(lái)區(qū)分混沌和噪聲信號(hào)．也就是對(duì)于一個(gè)混沌時(shí)間序列，它的預(yù)測(cè)值和實(shí)際值之間的相關(guān)系數(shù)會(huì)隨預(yù)測(cè)時(shí)間的增加而減弱;而對(duì)于隨機(jī)時(shí)間序列，這種相關(guān)性不會(huì)因預(yù)測(cè)時(shí)間的增加而變化．

本文針對(duì)由金屬薄壁構(gòu)件實(shí)測(cè)得到的具有不同非線性約束邊界條件、不同性質(zhì)的非平穩(wěn)振動(dòng)信號(hào)，采用兩種不同的非線性預(yù)測(cè)方法進(jìn)行對(duì)比分析，即代替數(shù)據(jù)法和Lyapunov指數(shù)方法．其中代替數(shù)據(jù)法主要參考了文獻(xiàn)［12］．這種方法可以通過(guò)對(duì)比混沌信號(hào)和噪聲信號(hào)(包括白噪聲和有色噪聲)的預(yù)測(cè)誤差分布及其所對(duì)應(yīng)的代替數(shù)據(jù)集的預(yù)測(cè)誤差分布的差異來(lái)實(shí)現(xiàn)．Lyapunov指數(shù)是某過(guò)程在其相平面內(nèi)相鄰軌線的平均發(fā)散速率的量化定義．正的Lyapunov指數(shù)(一個(gè)或多個(gè))是進(jìn)行混沌辨識(shí)的重要指標(biāo)．時(shí)間序列的Lyapunov指數(shù)的估算方法有多種，見(jiàn)文獻(xiàn)［13－17］．

1 兩種非線性預(yù)測(cè)方法的算法原理

1．1 基于代替數(shù)據(jù)法的混沌辨識(shí)

(1)信號(hào)的代替數(shù)據(jù)集的生成

信號(hào)所對(duì)應(yīng)的代替數(shù)據(jù)集可以由基于時(shí)間序列的 Gauss隨機(jī)過(guò)程假設(shè)得到．本文通過(guò)對(duì)非Gauss過(guò)程有效的非線性直方圖變換方法來(lái)實(shí)現(xiàn)．對(duì)于一個(gè)原始信號(hào)，一般生成包含128個(gè)不同時(shí)間

接下來(lái)，通過(guò)將上述復(fù)數(shù)乘以eiφ使其相位角在每個(gè)頻率上隨機(jī)化，其中φ是歸一化的、在［0，2π］區(qū)間內(nèi)變化的隨機(jī)量，得到新的X'(k)．對(duì)其進(jìn)行逆傅立葉變換，可以得到Gauss型代替數(shù)據(jù)序列，也就是得到具有原始信號(hào)相同的幅值分布形式的代替數(shù)據(jù)序列x'(n)，如下式序列的代替數(shù)據(jù)集．

進(jìn)行直方圖變換時(shí)，首先生成一個(gè)與原給定時(shí)間序列長(zhǎng)度相同的Gauss隨機(jī)數(shù)集合，然后對(duì)這個(gè)Gauss數(shù)據(jù)集的順序進(jìn)行重排．新的時(shí)間序列應(yīng)是具有Gauss概率密度函數(shù)分布的、與原序列相對(duì)應(yīng)的非線性尺度變換的結(jié)果．對(duì)于原始信號(hào)是Gauss隨機(jī)的情況，變換后的序列也具有Gauss分布．第一步是，將時(shí)間序列x(n)進(jìn)行傅立葉變換，得到

由于上式的逆傅立葉變換是實(shí)數(shù)，其相位角具有對(duì)稱性，即φ(k)=－φ(N－k)．

(2)統(tǒng)計(jì)量的計(jì)算

由于低維混沌意味著其系統(tǒng)在短期內(nèi)可以視為是確定性的，而隨機(jī)過(guò)程與此不同．可以將預(yù)測(cè)誤差ε取為統(tǒng)計(jì)量．首先利用狀態(tài)變量x將時(shí)間序列x(n)進(jìn)行相空間重構(gòu)，采用時(shí)間滯后的嵌入法，即

其中嵌入維數(shù)d應(yīng)該滿足d≥2D+1，D是系統(tǒng)的真實(shí)吸引子維數(shù)，在實(shí)際計(jì)算過(guò)程中往往需要依靠經(jīng)驗(yàn)選?。⒁獾竭@里的n小于數(shù)據(jù)長(zhǎng)度N，即1≤n≤N．時(shí)間滯后點(diǎn)τ的選擇有時(shí)具有一定隨意性，在這里簡(jiǎn)單地取作1．

下面將時(shí)間序列的數(shù)據(jù)點(diǎn)集分成長(zhǎng)度相等(Nf=Nt)的擬合集和檢驗(yàn)集兩部分．在擬合集中，尋找與當(dāng)前點(diǎn)歐幾里德距離最為相鄰的k個(gè)狀態(tài)點(diǎn)．這k個(gè)當(dāng)前時(shí)刻為m的狀態(tài)點(diǎn)與對(duì)應(yīng)的下一時(shí)刻m+1的狀態(tài)點(diǎn)組成如下k個(gè)點(diǎn)對(duì):

采用如(5)式所示的預(yù)測(cè)公式對(duì)這k個(gè)點(diǎn)對(duì)進(jìn)行擬合

其中，a和b是擬合系數(shù)，ˉx是x的預(yù)測(cè)點(diǎn)．

接下來(lái)，計(jì)算檢驗(yàn)集Nt中的所有點(diǎn)對(duì)應(yīng)的預(yù)測(cè)誤差，將得到N/2個(gè)預(yù)測(cè)誤差．預(yù)測(cè)誤差的定義為由(6)式得到的下一時(shí)刻(m+1)的預(yù)測(cè)值和實(shí)際值的差值．

進(jìn)行混沌識(shí)別的統(tǒng)計(jì)量取為上述全部預(yù)測(cè)誤差的平均絕對(duì)誤差(MAE)．

(3)假設(shè)檢驗(yàn)

對(duì)于生成的所有代替數(shù)據(jù)集，根據(jù)如下公式(7)來(lái)計(jì)算顯著度:

其中，QD是由原始信號(hào)時(shí)間序列計(jì)算得到的統(tǒng)計(jì)量值(MAE值)，us和σs分別是由生成的128個(gè)代替數(shù)據(jù)序列計(jì)算得到的統(tǒng)計(jì)量值的均值和方差．

計(jì)算得到的χ值可以用于分析原始信號(hào)數(shù)據(jù)和代替數(shù)據(jù)的差異．如果χ是一個(gè)較小的數(shù)，這意味著原始信號(hào)和它的所有代替數(shù)據(jù)集具有相同的性質(zhì)，因此隨機(jī)假設(shè)可以接受，也就是原始信號(hào)是隨機(jī)的．相反，如果χ值較大，可以認(rèn)為代替數(shù)據(jù)序列與原始信號(hào)有較大的差別，拒絕隨機(jī)假設(shè)．

更進(jìn)一步，為了辨識(shí)原始數(shù)據(jù)序列是隨機(jī)的還是混沌的，定義如式(8)所示的置信判據(jù)．拒絕隨機(jī)假設(shè)的最大概率也就是相應(yīng)的顯著度P定義為

其中erf()為數(shù)據(jù)序列所有元素的誤差函數(shù)，其定義如式(9)所示．

根據(jù)經(jīng)驗(yàn)，如果計(jì)算得到的概率P值小于0．05，原始信號(hào)數(shù)據(jù)將顯著地不同于它的代替數(shù)據(jù)集，這時(shí)可以拒絕隨機(jī)假設(shè)，認(rèn)為原始信號(hào)在95%置信度下是混沌的．如果P值大于0．05，則認(rèn)為原始信號(hào)是隨機(jī)的(95%置信度)．當(dāng)然，這里的臨界值0．05可以根據(jù)實(shí)際情況的不同而不同．

1．2 Lyapunov指數(shù)的估計(jì)算法

Lyapunov指數(shù)的估計(jì)算法也是基于非線性預(yù)測(cè)理論的．進(jìn)行相空間重構(gòu)后，考慮兩條具有不同初始點(diǎn)的相鄰軌線L1和L2，它們的初始點(diǎn)分別是x0和z0．這兩個(gè)初始點(diǎn)之間的距離d0=|z0－x0|．經(jīng)過(guò)時(shí)間Δt后，x0和z0將沿著各自的軌線到達(dá)x1和y1，這時(shí)新的距離d1=|y1－x1|．在x1和y1之間選擇一個(gè)新點(diǎn)z1，并設(shè)d0=|z1－x1|．點(diǎn)x1和z1分別位于L1和L3之上．再過(guò)時(shí)間Δt，以x1和z1為起始點(diǎn)的軌線L1和L3將到達(dá)它們對(duì)應(yīng)的新點(diǎn)x2和y2．這樣經(jīng)過(guò)p次重復(fù)，將得到di=|yi－xi|(i=1，2，…，p)．Lyapunov指數(shù)的計(jì)算公式如式(10)所示．

按數(shù)值大小進(jìn)行重排，得到該時(shí)間序列的Lyapunov指數(shù)如下

因?yàn)閜是一個(gè)較大的整數(shù)，由此得到的Lyapunov指數(shù)是每個(gè)相點(diǎn)在其軌線上以指數(shù)形式發(fā)散的統(tǒng)計(jì)平均值．在上述Lyapunov指數(shù)中，一個(gè)或多個(gè)Lyapunov指數(shù)可能都是正的．根據(jù)非線性理論，正的Lyapunov指數(shù)意味著該時(shí)間序列(信號(hào))是混沌的．

2 實(shí)例

對(duì)某薄壁構(gòu)件進(jìn)行寬頻隨機(jī)激振試驗(yàn)．外激勵(lì)為500Hz范圍內(nèi)的標(biāo)準(zhǔn)白噪聲，對(duì)構(gòu)件的振動(dòng)響應(yīng)進(jìn)行測(cè)試．由于是兩種性質(zhì)不同的非線性邊界條件，得到的構(gòu)件振動(dòng)響應(yīng)呈現(xiàn)出明顯的非平穩(wěn)振動(dòng)信號(hào)特征，其中應(yīng)該蘊(yùn)含著不同的特征．

圖1所示為測(cè)試得到的兩組振動(dòng)響應(yīng)時(shí)間信號(hào)，分別記為x1和x2．它們所對(duì)應(yīng)的功率譜密度(PSD)如圖2所示．從圖(2)只能看出，這兩個(gè)信號(hào)都具有0－500Hz有限頻帶分布的特征，形狀上僅有稍微區(qū)別．因此，僅從圖1和圖2人們很難區(qū)分這兩個(gè)信號(hào)．

首先從非線性定性分析的角度對(duì)這兩組信號(hào)進(jìn)行比較．可以利用連續(xù)峰值來(lái)繪制這兩個(gè)信號(hào)的偽Poincare映射圖．設(shè)是一個(gè)時(shí)間序列的峰值集合，則它的偽Poincare映射圖是指對(duì)應(yīng)繪出的點(diǎn)圖．如果考慮采樣誤差，信號(hào)的周期點(diǎn)會(huì)在偽Poincare映射圖表現(xiàn)為一個(gè)較小的區(qū)域．而偽Poincare映射圖中出現(xiàn)分散的點(diǎn)區(qū)域時(shí)，根據(jù)非線性混沌理論，表明存在不規(guī)則或奇怪吸引子．對(duì)于這兩個(gè)信號(hào)x1和x2，它們的偽Poincare映射圖如圖3所示．從圖3(a)和(b)可以看出，信號(hào)x1和x2的偽Poincare映射圖在形貌上還是有區(qū)別的．

圖1 非平穩(wěn)時(shí)間序列信號(hào)x1和x2Fig．1 Non － stationary time series of x1and x2

圖2 x1和x2的功率譜密度Fig．2 Power spectral densities of x1and x2

根據(jù)前文所介紹的非線性混沌分析理論，用代替數(shù)據(jù)法計(jì)算得到的關(guān)于x1和x2的概率P值見(jiàn)表1．如1．1節(jié)所言，如果P值小于0．05，則拒絕隨機(jī)假設(shè)，信號(hào)的混沌特性得到辨識(shí)(95%的置信度)．因此信號(hào)x1和x2都可以認(rèn)為是混沌時(shí)間序列．另外，從表1中的其它特征值來(lái)看，如QD，us和χ，也可以對(duì)這兩組信號(hào)加以定量地對(duì)比區(qū)分．

圖3 x1和x2偽Poincare映射圖Fig．3 Pseudo Poincare mapping portraits of x1and x2

計(jì)算得到的x1和x2的Lyapunov指數(shù)如表2所示．在估算這兩組信號(hào)的Lyapunov指數(shù)的過(guò)程中，相空間重構(gòu)的嵌入維數(shù)設(shè)為5，滯后點(diǎn)數(shù)設(shè)為7．從表2可以看出，x1信號(hào)有2個(gè)正的Lyapunov指數(shù)，即 1．1785、0．1846，信號(hào)x2有 3 個(gè)正的 Lyapunov指數(shù)，即 0．9674、0．1857 和 0．0184．所以，這兩組信號(hào)可以視為是混沌的，但具有不同的混沌階數(shù)．

表1 代替數(shù)據(jù)法得到的x1和x2的特征參數(shù)數(shù)值Table 1 Characteristic parameter values of x1and x2by surrogate data method

表2 估算得到的x1和x2的Lyapunov指數(shù)Table 2 Lyapunov index of estimated x1and x2

3 結(jié)論

本文利用基于混沌理論的非線性預(yù)測(cè)方法對(duì)非平穩(wěn)信號(hào)進(jìn)行了比較分析．從定量分析研究的角度，所采用的代替數(shù)據(jù)法和Lyapunov指數(shù)估計(jì)方法對(duì)區(qū)分具有不同非線性性質(zhì)的非平穩(wěn)信號(hào)是有效的．

對(duì)于兩個(gè)典型的非平穩(wěn)振動(dòng)信號(hào)x1和x2，代替數(shù)據(jù)法所得到的特征參數(shù)的數(shù)值是互不相同的，估算得到的最大Lyapunov指數(shù)也是正的．根據(jù)代替數(shù)據(jù)法中的概率值的大小可以看出，x1比x2具有更明顯的混沌特征，但是x1有2個(gè)正的Lyapunov指數(shù)，而x2有3個(gè)正的Lyapunov指數(shù)．

1 Kantz H，Schreiber T．Nonlinear time series analysis．Cambridge:Cambridge University Press，1999

2 張賢達(dá)，保錚．非平穩(wěn)信號(hào)分析與處理．北京:國(guó)防工業(yè)出版社，1998(Zhang X D，Bao Z．Nonstationary signal analysis and processing．Beijing:National Defence Industry Press，1998(in Chinese))

3 He Zhengjia，Meng Qingfeng，Zhao Juyuan，et al．Timefrequency(scale)analysis and diagnosis for nonstationary dynamic signal of machinery．International Journal of Plant Engineering and Management，1996，1(1):40 ～47

4 Meltzer G S，Ivanov Y Y．Fault detection in gear drives with non－stationary rotational speed－ Part I:the time－frequency approach．Mechanical Systems and Signal Processing，2003，17(5):1033 ～1047

5 Wang Wenyi，Harrap M J．Condition monitoring of rolling element bearings by using cone kernel time－frequency distribution．In:Processings of the SPIE Conference on Measurement Technology and Intelligent Instrument．Bellingham，WA，1993:290～298

6 Katvonik V，Stankovic L．Instantaneous frequency estimation using the Wigner distribution with varying and datadriven widow length．IEEE Transactions on Signal Processing，1998，46(9):2351 ～2325

7 Casdagli M．Nonlinear prediction of chaotic time series．Physica D，1989，35(3):335～356

8 Engbert R．Testing for nonlinearity:the role of surrogate data．Chaos，Solitons＆Fractals，2002，13(1):79 ～84

9 Banbrook M，Ushaw G，McLaughlin S．Lyapunov exponents from a time series:a noise－robust extraction algorithm．Chaos，Solitons＆Fractals，1996，7(7):973 ～976

10 Theiler J，Eubank S，Longtin A，et al．Testing for nonlin－earity in time series:the method of surrogate data．Physica D:Nonlinear Phenomena，1992，58(1－4):77 ～94

11 Yonemoto K，Yanagawa T．Estimating the Lyapunov exponent from chaotic time series with dynamic noise．Statistical Methodology，2007，4(4):461 ～480

12 Kennel M，Isabelle S．Method to distinguish possible chaos from coloured noise and to determine embedding parameters．Physical Review A，1992，46(6):3111 ～3118

13 Kantz H．A robust method to estimate the maximal Lyapunov exponent of a time series．Physics Letters A，1994，185(1):77～87

14 Wolf A，Swift J V，Swinney H L，et al．Determining Lyapunov exponents from a time series．Physica D:Nonlinear Phenomena，1985，16(3):285 ～317

15 Rosenstein M T，Collins J J，Luca C J．A practical method for calculating largest Lyapunov exponents from small data sets．Physica D:Nonlinear Phenomena，1993，65(1－2):117～134

16 周良強(qiáng)，陳芳啟，葛鵬．一類大展弦比機(jī)翼非線性模型的穩(wěn)定性與分岔．動(dòng)力學(xué)與控制學(xué)報(bào)，2012，10(3):222 ～227(Zhou L Q，Chen Q F，Ge P．Stability and bifurcation for a nonlinear model of high－aspect－ratio wings．Journal of Dynamics and Control，2012，10(3):222～27(in Chinese))

17 呂書(shū)鋒，胡宇達(dá)．正交各向異性疊層板的非線性主共振分析．動(dòng)力學(xué)與控制學(xué)報(bào)，2009，7(1):35～38(Lu S F，Hu Y D．Nonlinear principal resonance of orthotropic laminated plates．Journal of Dynamics and Control，2009，7(1):35～38(in Chinese))

*The project supported by the National Natural Science Foundation of China(10972192)

? Corresponding author E－mail:han1011@163．com

COMPARATIVE ANALYSES OF NON-STATIONARY SIGNALS BASED ON NONLINEAR CHAOTIC THEORIES*

Wei Chunyu Yan WeiWang Zihuimin Han Qingpeng?
(College of Mechanical Engineering，Liaoning Science and Technology University，Anshan114051，China)

In the paper，two nonlinear estimation methods based on nonlinear chaotic theory，surrogate data method and Lyapunov exponents，are used to distinguish the difference of non－stationary signals．After brief introduction of the corresponding algorithms，two typical different non－stationary signals measured from a thin－plate structure with different nonlinear constraining boundaries are taken to compare by using the above two methods respectively．The obtained results demonstrate that the apparently similar signals are distinguished effectively in quantitative way with applying above nonlinear chaotic analyses．

nonlinear chaotic theory， non－stationary signals， surrogate data method， Lyapunov exponents

9 October 2012，

1 July 2013．

10．6052/1672－6553－2014－014

2012－10－09 收到第 1 稿，2013－07－01 收到修改稿．

*國(guó)家自然科學(xué)基金資助項(xiàng)目(10972192)

E－mail:han1011@163．com