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

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

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

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

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

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

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

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

接下來，計算檢驗集Nt中的所有點對應(yīng)的預(yù)測誤差，將得到N/2個預(yù)測誤差．預(yù)測誤差的定義為由(6)式得到的下一時刻(m+1)的預(yù)測值和實際值的差值．

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

(3)假設(shè)檢驗

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

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

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

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

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

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

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

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

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

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

2 實例

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

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

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

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

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

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

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

計算得到的x1和x2的Lyapunov指數(shù)如表2所示．在估算這兩組信號的Lyapunov指數(shù)的過程中，相空間重構(gòu)的嵌入維數(shù)設(shè)為5，滯后點數(shù)設(shè)為7．從表2可以看出，x1信號有2個正的Lyapunov指數(shù)，即 1．1785、0．1846，信號x2有 3 個正的 Lyapunov指數(shù)，即 0．9674、0．1857 和 0．0184．所以，這兩組信號可以視為是混沌的，但具有不同的混沌階數(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ù)測方法對非平穩(wěn)信號進(jìn)行了比較分析．從定量分析研究的角度，所采用的代替數(shù)據(jù)法和Lyapunov指數(shù)估計方法對區(qū)分具有不同非線性性質(zhì)的非平穩(wěn)信號是有效的．

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

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

2 張賢達(dá)，保錚．非平穩(wěn)信號分析與處理．北京:國防工業(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)定性與分岔．動力學(xué)與控制學(xué)報，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 呂書鋒，胡宇達(dá)．正交各向異性疊層板的非線性主共振分析．動力學(xué)與控制學(xué)報，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 收到修改稿．

*國家自然科學(xué)基金資助項目(10972192)

E－mail:han1011@163．com