基于相空間重構(gòu)改進(jìn)算法的混沌奇異譜分析及應(yīng)用

張立國， 劉 婉， 張淑清， 劉海濤， 董 偉， 宋姍姍

(燕山大學(xué)電氣工程學(xué)院，河北秦皇島066004)

1 引 言

奇異譜通過把時間序列信號向嵌入維數(shù)主軸分解，分解的每一個主成分代表著信號的特征成分[1]。信號主成分的大小分布反映出信號能量在對應(yīng)嵌入維數(shù)主軸上的大小分布情況，主成分各個狀態(tài)變量反映在整個系統(tǒng)中所占能量的相對關(guān)系[2]。因此，奇異譜技術(shù)不僅可以對動力系統(tǒng)的時間序列進(jìn)行降噪，還可以作為動力系統(tǒng)特征提取的有效方法[3]?；煦缙娈愖V通過相空間重構(gòu)，得到各主分量的圖像譜值分布，從而提取混沌信號的特征信息。

相空間重構(gòu)作為混沌時間序列處理中的重要課題，其最佳嵌入維數(shù)m和延遲時間τ的選取至關(guān)重要。如果m取值過小，吸引子會發(fā)生折疊甚至在某些地方會出現(xiàn)自相交，重構(gòu)吸引子的幾何形狀和原始狀態(tài)吸引子可能完全不同；如果m取值過大，當(dāng)維數(shù)m大于最小嵌入維數(shù)的時候，幾何結(jié)構(gòu)已經(jīng)被完全打開，此時這些幾何不變量與嵌入的維數(shù)無關(guān)；如果τ取值太小，那么在相空間中各向量的分量間幾乎不包括新信息，從而低估關(guān)聯(lián)維數(shù)。相反，τ取值太大，相空間重構(gòu)的相關(guān)有效信息會遺漏，這樣會高估關(guān)聯(lián)維數(shù)[4]。

m的確定方法很多，如飽和關(guān)聯(lián)維數(shù)法、鄰近點維數(shù)法、虛假鄰近點法[5]。這些方法確定最佳嵌入維的原理大致相同，都是依據(jù)隨嵌入維數(shù)m升高而逐步收斂的情況，但是沒有給出如何確定m最終值的依據(jù)。

本文提出基于改進(jìn)Cao算法[6]對m和τ的求取改進(jìn)方法，解決了常用求取方法中所存在的計算方法復(fù)雜、計算量龐大、計算結(jié)果不精確等問題，并且可以快速有效地提取混沌動力系統(tǒng)中有用定量信息，重構(gòu)系統(tǒng)的相空間。通過對Lorenz典型混沌系統(tǒng)進(jìn)行數(shù)值仿真試驗，結(jié)果表明，這兩種方法結(jié)合能夠有效地重構(gòu)原系統(tǒng)的相空間。

在求取m和τ基礎(chǔ)上，對混沌奇異譜特征量化特征進(jìn)行分析，通過實驗驗證奇異譜分析技術(shù)不僅具有很強的穩(wěn)定性以及抗噪性能，還能對不同信號進(jìn)行明顯區(qū)分。

本文方法在滾動軸承早期故障識別中應(yīng)用，首先通過最大Lyapunov指數(shù)[7,8]判別機械故障數(shù)據(jù)具有混沌特性；然后，利用本文改進(jìn)的混沌相空間重構(gòu)方法得到最佳延遲時間和嵌入維數(shù)；最后得到混沌奇異譜，通過混沌奇異譜提取故障特征，將不同故障識別出來。為機械故障早期診斷提供一種新的有效途徑。

2 基于CAO算法的混沌相空間重構(gòu)

設(shè)x1,x2,…,xN是一個時間序列，并且其重構(gòu)相空間向量[9]：

Ym(i)=[x(i),x(i+τ),…,x(i+(m-1)τ)]
i=1,2,…,N-(m-1)τ

(1)

(2)

定義：

(3)

(4)

E(m)僅取決于嵌入維數(shù)m和延遲時間τ，為了研究E(m)從m到m+1的變化，定義

E1(m)=E(m+1)/E(m)

(5)

從而發(fā)現(xiàn)當(dāng)m大于某個值m0時，E1(m)停止變化。如果吸引子可以從m維相空間重構(gòu)中獲得，那么m0就是尋找的最小嵌入維數(shù)。

在進(jìn)行E*(m)數(shù)值驗證之前，需要定義另一個用于區(qū)分確定性信號和隨機信號E2(m)的量：

(6)

E2(m)=E*(m+1)/E*(m)

(7)

3 改進(jìn)的Cao算法

3.1 Cao算法存在的問題

當(dāng)該時間序列為一串隨機數(shù)字時，隨著m的增加，E1(m)永遠(yuǎn)不會達(dá)到飽和值。但是在實際計算中，當(dāng)m足夠大時，很難去判斷E1(m)緩慢增加還是停止改變。事實上，由于可觀測的數(shù)據(jù)樣本是有限的，即使時間序列是隨機的，也有可能E1(m)在某個維度是停止變化的。

為了解決這個問題，增加了E2(m)這個參考量。對于隨機數(shù)據(jù)，未來值對于過去值是無關(guān)的，因此，這種情況下，對于任何嵌入維數(shù)，E2(m)將趨近于1。然而，對于確定的時間序列，E2(m)是與m有關(guān)的，即對于所有的m，E2(m)不可能為某一定值。那么在不同的嵌入維數(shù)下，必存在一些E2(m)≠1的值。E1(m)用于計算確定時間序列的最小嵌入維數(shù)，E2(m)則是在實際應(yīng)用中判斷有序數(shù)列是在緩慢變化還是已經(jīng)趨向于穩(wěn)定。E1(m)和E2(m)的值在m大于某一特定值m0時停止變化。m0即為最佳嵌入維數(shù)。

3.2 確定嵌入維數(shù)的改進(jìn)方法

嵌入維數(shù)m的確定是依據(jù)E1(m)停止變化為標(biāo)準(zhǔn)的，這個標(biāo)準(zhǔn)并沒有給出判斷停止變化的準(zhǔn)則，只是依靠主觀判斷，實際上，時間序列E1(m)經(jīng)常是有起伏的，很少出現(xiàn)嚴(yán)格意義上的停止變化，所以，這給m的確定帶來了困難。

針對這一問題，本文利用補充準(zhǔn)則E2(m)進(jìn)行聯(lián)合判斷，同時對CAO算法提出了改進(jìn)，給出了一種改進(jìn)的嵌入維數(shù)的穩(wěn)定性準(zhǔn)則[10]，計算步驟如下：

1) 計算Δi：

Δi=|E1(i)-E1(i+1)|, 1≤i≤N-1

(8)

式中：|E1(i)-E1(i+1)|表示計算E1(i)與E1(i+1)的差值的絕對值，結(jié)果用Δi來表示。

2) 根據(jù)E1(i)的波動情況選取一個閾值e，找到第一個Δi

Δi=max(Δi);n≤i≤N-1

(9)

3) 重新設(shè)置：

(10)

4) 取j≤i≤N-2，當(dāng)滿足Δi>Δi+1，Δi+1>Δi+2，Δi

同理，當(dāng)用E2(i)、E2(i+1)換成計算式(8)中的E1(i)、E1(i+1)時，可計算出由式(5)、式(6)確定的關(guān)于E2(i)的m值。

3.3 基于符號分析的極大聯(lián)合熵延遲時間求取方法

3.3.1 最佳延遲時間的極大聯(lián)合熵準(zhǔn)則

考慮時間序列x(n)及其延遲時間序列xτ(n)=x(n+τ)，n=1,2,…,N。根據(jù)互信息函數(shù)的遞推公式[11]，兩組序列的互信息可表示為

I(X,Xτ)=H(X)-H(X|Xτ)=H(X)+
H(Xτ)-H(X|Xτ)=I(Xτ,X)

(11)

式中：X代表時間序列X(n)；Xτ代表其延遲序列X(n+τ)；H(X)是孤立的X的不定性，H(X|Xτ)是已知Xτ的X的不定性，所以Xτ的已知減少了X的不定性，則I(X,Xτ)的第一極小值處的τ即為最佳延遲時間。并且為了計算I(X,Xτ)，F(xiàn)raser A M等[12]提出了復(fù)雜的劃分網(wǎng)格的方法。

互信息I(X,Xτ)表征X和Xτ的相關(guān)程度，I(X,Xτ)取極小值時X和Xτ的相關(guān)度也極小，即X和Xτ的聯(lián)合整體不確定度達(dá)到極大[13]。

為了保證重構(gòu)坐標(biāo)之間最大限度地相互獨立，應(yīng)使X和Xτ聯(lián)合整體的不確定性達(dá)到最大。由信息理論可知，聯(lián)合熵H(X|Xτ)是X和Xτ的聯(lián)合整體的不確定性的度量，H(X|Xτ)越大，則X和Xτ聯(lián)合整體的不確定性也越大。由此推斷，聯(lián)合熵H(X|Xτ)的第一個極大值點即為相空間重構(gòu)的最佳延遲時間點[14]。

H(X,Xτ)=H(X)+H(Xτ)-I(X,Xτ)

=Const-I(X,Xτ)

(12)

由此可見，H(X,Xτ)和I(X,Xτ)呈近似相反的變化規(guī)律。互信息的極小值點即為聯(lián)合熵的極大值點在理論上得到驗證。

因此,通過復(fù)雜的劃分網(wǎng)格或者進(jìn)行網(wǎng)格標(biāo)記求取互信息第一極小值點，從而確定相空間重構(gòu)最佳延遲時間可以轉(zhuǎn)換為求取聯(lián)合熵H(X,Xτ)的第一極大值點。聯(lián)合熵的計算公式[15]可以表示為

(13)

P(xi,xtj)為變量xi、xtj的聯(lián)合概率分布，由此可見，求取不同延遲時間下的聯(lián)合熵[16]，找出聯(lián)合熵的第一極大值點所對應(yīng)的τ，即為所求的最佳延遲時間點。

3.3.2 符號分析法求取極大聯(lián)合熵

二進(jìn)制劃分是最簡單的一種劃分規(guī)則，只需給定一個閾值P0，時間序列中大于此閾值P0的取1，否則取0。閾值P0的選取有均值、零值等。則時間序列X(n)最終被轉(zhuǎn)換成符號序列{S(n)}，所得符號序列表征了2種數(shù)據(jù)元模式，即

(14)

時間序列信號經(jīng)二進(jìn)制符號化規(guī)則進(jìn)行處理，其具體過程描述如圖1所示。

圖1 二進(jìn)制符號化規(guī)則Fig.1 Binary symbolization rule

按照上式進(jìn)行標(biāo)記和辨識，則離散的符號替代連續(xù)數(shù)據(jù)的原始時間序列，見圖2。

圖2 符號化時間序列圖Fig.2 Symbolized time series diagram

符號數(shù)d可以通過使符號熵最大化來尋找，為了方便起見，將混沌序列的臨界點(d+1)的個數(shù)置為11，即d=10。且分割長度L取2。

對時間序列X(n)及其延遲時間序列X(n+τ)，n=1,2,…,N，根據(jù)上面所介紹的粗?；柗椒▽⑵渚幋a成能捕獲有用定量信息的特殊的十進(jìn)制數(shù)序列Lx(n)和Lxτ(n)，則各個特殊十進(jìn)制數(shù)出現(xiàn)的頻率為時間序列分析的指標(biāo)，即為聯(lián)合概率P(Lx,Lxτ)，則符號分析法求取的聯(lián)合熵公式(13)可以改寫為

(15)

4 數(shù)值驗證

為了驗證這2種方法結(jié)合求取最佳延遲時間和最佳嵌入維數(shù)的準(zhǔn)確性和優(yōu)越性，對常見的典型混沌時間序列Lorenz進(jìn)行數(shù)值仿真,求取這兩種典型混沌系統(tǒng)的最佳延遲時間和最佳嵌入維數(shù)并畫出其吸引子的重構(gòu)。

Lorenz系統(tǒng)的數(shù)值仿真實驗如圖3所示。

圖3 Lorenz系統(tǒng)(σ=16, b=4, r=45.92)的數(shù)值仿真實驗Fig.3 Numerical simulation experiment of Lorenz system (σ=16, b=4, r=45.92)

該系統(tǒng)方程可描述[17]為

(16)

式中：x、y、z分別是三維坐標(biāo)；σ=16；b=4；r=45.92。

選取Lorenz系統(tǒng)參數(shù)確定情況下該系統(tǒng)為混沌系統(tǒng)，用Runge-Kutta法求解方程(16)，步長h=0.01，取變量x為研究對象，去除前8 000個暫態(tài)點，得到一個7 000個點的時間序列。

通過本文所提出的符號分析法求取該系統(tǒng)時間序列聯(lián)合熵H(Lx,Lxτ)的極大值點，求取時間序列的有效延遲時間的圖形如圖3(a)所示。從圖中可以看出，聯(lián)合熵H(Lx,Lxτ)在τ=11時取得第一個所對應(yīng)的極值點，故可得由聯(lián)合熵第一極大值點法求得的相空間重構(gòu)最佳延遲時間；試驗過程中也大大地簡化了其計算量。圖3(b)是由Cao方法求得的E1&E2的最佳嵌入維數(shù)，圖3(c)是由改進(jìn)后Cao方法求得E1的最佳嵌入維數(shù)，圖3(d)是由改進(jìn)后Cao方法求得E2的最佳嵌入維數(shù)，從圖3(b)中可以看出E1、E2都在m大于3后看似不再發(fā)生變化，實際m在3之后有微小波動，經(jīng)過本文提出的改進(jìn)Cao算法，取e=0.1,如圖3(c)確定出E1的最佳嵌入維數(shù)為m=12，如圖(d)E2的最佳嵌入維數(shù)為m=18，所以該時間序列的最佳嵌入維數(shù)為18；圖 3(e)是m=12時LorenzX相三維重構(gòu)吸引子圖；圖3(f)是m=12時的重構(gòu)吸引子圖；圖3(g)是m=18時LorenzX相三維重構(gòu)吸引子圖；圖3(h)是m=18時的重構(gòu)吸引子圖。從圖中可以看出m=12時Lorenz混沌時間序列在X相的三維重構(gòu)吸引子圖和平面的重構(gòu)吸引子圖都無法完全展開，而m=18時重構(gòu)圖能夠很好的反映出Lorenz系統(tǒng)的雙圈拓?fù)浣Y(jié)構(gòu)。從而驗證了改進(jìn)Cao算法確定最佳嵌入維數(shù)的有效性和優(yōu)越性。

5 混沌奇異譜特性分析

5.1 混沌奇異譜特征量化分析

由嵌入定理Takens可知，對于混沌時間序列x(t)，可以得出m維相空間的重構(gòu)吸引子。于是可以得到軌道矩陣：

(17)

式中：N=n-m-1為相點個數(shù)，每一個相點坐標(biāo)為

X(j)=[x(j),x(j+τ),…,x(j+(m-1)τ)]

(18)

設(shè)混沌動力系統(tǒng)的m個狀態(tài)變量相互獨立且相互正交，狀態(tài)空間分布在m維嵌入空間向量的主分量方向上，則存在變換矩陣A，使得YTTT的協(xié)方差矩陣為對角陣，即

(19)

式中:CX和CY分別是和的延時-協(xié)變矩陣，Sj(j=1,2,…,m)為CY特征值。將特征值由大到小排列，即S1≥S2≥…≥Sm，記

(20)

稱S1≥S2≥…≥Sm為系統(tǒng)的奇異譜，它表示各個狀態(tài)變量在整個系統(tǒng)中所占能量的相對關(guān)系。式(20)中對于延時-協(xié)變矩陣CX定義為

(21)

式中：c(j)表示延時為時x(n)的協(xié)方差，其計算公式如下：

(22)

由協(xié)方差公式可知，矩陣CX是一個正定對稱矩陣，因此其特征值非負(fù)，即e1≥e2≥…≥em>0，這些特征向量Ek稱為經(jīng)驗正交函數(shù)(EOF)，第k個主成分(PC)定義為混沌時間序列x(t)在第k個經(jīng)驗正交函數(shù)Ek上的正交投影系數(shù)。從頻譜上來看，一個EOF對應(yīng)著一個自適應(yīng)滑動平均濾波器，每一個PC代表著信號的特征成分。奇異譜的意義在于：通過它把時間序列信號向嵌入維數(shù)主軸分解，信號主成分的大小分布反映出信號能量在對應(yīng)嵌入維數(shù)主軸上的大小分布情況。

5.2 奇異譜的特性驗證

針對Lorenz系統(tǒng)，取m=18。實驗采用10對分量不同時刻的時間序列，序列點數(shù)為5 000，分別對每條序列進(jìn)行相空間重構(gòu)，計算奇異譜，將10對奇異譜對應(yīng)的主成分分量放在同一個譜圖中，實驗結(jié)果如圖4所示。

圖4 Lorenz系統(tǒng)x分量的10對奇異譜圖Fig.4 10 pairs of singular spectra of the x component of the Lorenz system

以x分量時間序列{xn}(n=1,2,…,N,N為時間序列長度)為研究對象。向時間序列{xn}添加10 db噪聲，采用m=18，τ=11。取x(t)在不同時刻的10對時間序列，每個序列點數(shù)為5 000，然后分別對每個序列進(jìn)行相空間重構(gòu)，計算奇異譜，將10對奇異譜對應(yīng)的主成分分量放在同一個譜圖中，實驗結(jié)果如圖5所示。

圖5 Lorenz系統(tǒng)x分量的10對含10 db噪聲的奇異譜圖Fig.5 10 pairs of singular spectra with 10 db noise for the x component of the Lorenz system

接著再對Lorenz系統(tǒng)中x方向的一條時間序列，點數(shù)為5 000，分別加入10 db～20 db噪聲。采用m=18，τ=11。計算其奇異譜，得到奇異譜圖如圖6所示。

圖6 Lorenz系統(tǒng)x分量的不同噪聲含量的奇異譜圖Fig.6 Singular spectrum of different noise content of the x component of the Lorenz system

從圖4、圖5和圖6可以看出，奇異譜分析不僅具有很強的穩(wěn)定性及抗噪性能，還能將不同信號明顯區(qū)分。一方面對于含有噪聲的信號可以成功提取有用信息，另一方面可以將不同的有用信息有效區(qū)分。

6 混沌奇異譜分析在滾動軸承故障診斷中的應(yīng)用

6.1 軸承振動信號的混沌判別

首先，對4種軸承狀態(tài)振動信號選用10對不同時間段的樣本，時間序列的長度為2 400；然后，用本文方法分別計算樣本序列的嵌入維數(shù)和延遲時間，將它們作為小數(shù)據(jù)法的輸入，從而估計出樣本序列的最大Lyapunov指數(shù)。

圖7為不同狀態(tài)的最大Lyapunov指數(shù)，觀察可知不同狀態(tài)下振動信號的時間序列的最大Lyapunov指數(shù)均大于0，表明振動信號具有混沌特性，可以利用混沌相關(guān)理論對其進(jìn)行分析。

圖7 不同狀態(tài)下的最大Lyapunov指數(shù)Fig.7 Maximum “Lyapunov” index in different states

6.2 軸承振動信號的奇異譜

從圖7中可以看出，不同時間段的外圈故障的最大Lyapunov指數(shù)波動比較大，正常狀態(tài)和內(nèi)圈狀態(tài)下時間序列的最大Lyapunov指數(shù)交織在一起，僅利用設(shè)置臨界值的方法很難直接將不同狀態(tài)的信號分開。故需要有效的特征提取方法提取故障特征。

針對4種軸承類型下的振動信號，每種選擇50個樣本，樣本長度均為2 400，采用本文提出的混沌奇異譜方法。將數(shù)據(jù)代入嵌入維數(shù)的Cao方法和基于符號分析的極大聯(lián)合熵求延遲時間的方法，取所有數(shù)據(jù)的最大嵌入維數(shù)和延遲時間，計算奇異譜，最終相空間重構(gòu)過程中選用m=15，τ=6，每種類型隨機選取4個序列，繪制4種振動信號的奇異譜圖如圖8所示。

圖8 不同損傷位置的奇異譜圖Fig.8 Singular spectrum of different damage locations

從圖8中看出，不同損傷位置對應(yīng)奇異譜圖能量分布明顯不一致，但相同故障類型的奇異譜圖分布較為一致，上下波動不大，各個主成分的譜值基本一致，說明同一故障類型的軸承振動信號具有比較穩(wěn)定的奇異譜。所以，混沌奇異譜可以代表不同損傷位置振動信號的內(nèi)在動力學(xué)特征。

為了確定上述實驗效果，選擇電機驅(qū)動端振動傳感器采集的正常狀態(tài)、不同位置損傷、不同損傷直徑等10種狀態(tài)下的振動信號，采樣頻率為12 kHz，轉(zhuǎn)速為1 772 r/min，負(fù)載為0.735 N·m，下面實驗中的振動信號都是采取時間為0.2 s，共2 400個采樣點。

采用混沌奇異譜對驅(qū)動端軸承10種運行狀態(tài)進(jìn)行故障診斷，每種運行狀態(tài)選50組無重疊數(shù)據(jù)，每組數(shù)據(jù)的長度為2 400。首先計算10組原始振動信號數(shù)據(jù)的混沌奇異譜，將數(shù)據(jù)代入求嵌入維數(shù)的Cao方法和基于符號分析的極大聯(lián)合熵求延遲時間的方法，取所有數(shù)據(jù)的最大嵌入維數(shù)和延遲時間，最終相空間重構(gòu)過程中嵌入維數(shù)選用m=20，時間延遲選擇τ=6，進(jìn)而產(chǎn)生500×20的矩陣，為不同運行狀態(tài)的數(shù)據(jù)貼上標(biāo)簽，狀態(tài)類別依次為1，2，3，4，5，6，7，8，9和10。每種運行狀態(tài)得到50個特征樣本，其中30個作為訓(xùn)練集，剩余20個作為測試集，訓(xùn)練集樣本都是采用的有標(biāo)簽樣本。圖9是10種狀態(tài)的混沌奇異譜特征向量圖。

圖9 10種狀態(tài)的混沌奇異譜特征向量Fig.9 Chaotic singular spectral eigenvectors of 10 states

從圖9中看出，不同運行狀態(tài)下的振動信號的混沌奇異譜特征都是前期譜值較大，后期譜值較小。同時，不同運行狀態(tài)的混沌奇異譜特征前期譜值相差很大，后期譜值卻相差不大，譜值的大小代表能量的分布，說明運行狀態(tài)發(fā)生改變后系統(tǒng)的能量分布也發(fā)生了改變。

由此可見，不同的故障類型和不同的故障尺寸，均可通過奇異譜分析提取有效的特征，并且有效地區(qū)分，說明奇異譜分析可以有效地應(yīng)用于機械故障診斷。

7 結(jié) 論

(1) 混沌奇異譜作為一個有效的特征提取方法，能夠?qū)⑿盘栔械牟煌煞滞ㄟ^相空間重構(gòu)后以圖形的形式進(jìn)行了有效區(qū)分，線型區(qū)分明顯，具有很好的準(zhǔn)確性及可觀性。

(2) 針對混沌奇異譜分析時間序列的延遲時間和嵌入維數(shù)存在的不足，利用基于符號分析的極大聯(lián)合熵方法，采用Cao改進(jìn)算法確定嵌入m，提高了準(zhǔn)確度，兩者準(zhǔn)確實現(xiàn)相空間重構(gòu)，取得了很好的準(zhǔn)確性、穩(wěn)定性及實效性。

(3) 研究表明，機械故障信號具有混沌的特性，本文方法已在滾動軸承早期故障識別中應(yīng)用，為機械故障早期診斷提供一種新的有效途徑。