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

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

1 引 言

奇異譜通過把時間序列信號向嵌入維數(shù)主軸分解，分解的每一個主成分代表著信號的特征成分[1]。信號主成分的大小分布反映出信號能量在對應(yīng)嵌入維數(shù)主軸上的大小分布情況，主成分各個狀態(tài)變量反映在整個系統(tǒng)中所占能量的相對關(guān)系[2]。因此，奇異譜技術(shù)不僅可以對動力系統(tǒng)的時間序列進行降噪，還可以作為動力系統(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ù)。

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

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

本文方法在滾動軸承早期故障識別中應(yīng)用，首先通過最大Lyapunov指數(shù)[7,8]判別機械故障數(shù)據(jù)具有混沌特性；然后，利用本文改進的混沌相空間重構(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ù)。

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

(6)

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

(7)

3 改進的Cao算法

3.1 Cao算法存在的問題

當(dāng)該時間序列為一串隨機數(shù)字時，隨著m的增加，E1(m)永遠不會達到飽和值。但是在實際計算中，當(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ù)的改進方法

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

針對這一問題，本文利用補充準(zhǔn)則E2(m)進行聯(lián)合判斷，同時對CAO算法提出了改進，給出了一種改進的嵌入維數(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)合整體不確定度達到極大[13]。

為了保證重構(gòu)坐標(biāo)之間最大限度地相互獨立，應(yīng)使X和Xτ聯(lián)合整體的不確定性達到最大。由信息理論可知，聯(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)格或者進行網(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)合熵

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

(14)

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

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

按照上式進行標(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ù)上面所介紹的粗粒化符號方法將其編碼成能捕獲有用定量信息的特殊的十進制數(shù)序列Lx(n)和Lxτ(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進行數(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)是由改進后Cao方法求得E1的最佳嵌入維數(shù)，圖3(d)是由改進后Cao方法求得E2的最佳嵌入維數(shù)，從圖3(b)中可以看出E1、E2都在m大于3后看似不再發(fā)生變化，實際m在3之后有微小波動，經(jīng)過本文提出的改進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)的雙圈拓撲結(jié)構(gòu)。從而驗證了改進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是一個正定對稱矩陣，因此其特征值非負，即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，分別對每條序列進行相空間重構(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，然后分別對每個序列進行相空間重構(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)理論對其進行分析。

圖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，負載為0.735 N·m，下面實驗中的振動信號都是采取時間為0.2 s，共2 400個采樣點。

采用混沌奇異譜對驅(qū)動端軸承10種運行狀態(tài)進行故障診斷，每種運行狀態(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，進而產(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)后以圖形的形式進行了有效區(qū)分，線型區(qū)分明顯，具有很好的準(zhǔn)確性及可觀性。

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

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