李兆飛，任小洪，黃臣程

(1.四川理工學(xué)院 自動化與電子信息學(xué)院，四川 自貢 643000；2.人工智能四川省重點實驗室，四川 自貢 643000)

研究滾動軸承振動的非線性演化動態(tài)特性有助于認識其內(nèi)在變化規(guī)律，并利用這些規(guī)律對其進行診斷和監(jiān)測，具有一定的普遍意義。隨著對軸承振動非線性機理的認識逐漸加深，采用關(guān)聯(lián)維數(shù)[1]、最大Lyapunov指數(shù)[2-3]、Kolmogorov熵及相關(guān)復(fù)雜性測度[4-6]等方法對振動的非線性混沌特性進行判斷，很好地解釋了軸承振動現(xiàn)象的復(fù)雜運動現(xiàn)象，判定軸承振動是一種在有限相空間軌道永不重復(fù)、形態(tài)復(fù)雜的運動。然而，這些特征量只能表現(xiàn)軸承振動的某方面特性，而且受信號長度及噪聲影響，這些特征量的計算結(jié)果會有所偏差，導(dǎo)致判斷錯誤。因此，嘗試采用遞歸圖、CLY方法和功率譜3種方法對軸承不同振動狀態(tài)實測信號的混沌特性進行試驗分析。

如前所述，目前對軸承振動的混沌特性研究多是基于混沌的辨別，但對其振動超混沌特性的分析還未見報道，分析軸承振動的超混沌特性有助于研究振動在相空間多個不同方向的復(fù)雜演化規(guī)律。通常通過計算實測滾動軸承振動信號的Lyapunov指數(shù)譜判斷其超混沌特性。Lyapunov指數(shù)譜的計算往往要通過重構(gòu)相空間，之后再判定未知系統(tǒng)方程狀態(tài)。判斷未知系統(tǒng)方程通常采用人工神經(jīng)網(wǎng)絡(luò)[7-9]，存在訓(xùn)練時間較長、易陷入局部極小點及過擬合，且確定網(wǎng)絡(luò)結(jié)構(gòu)困難的問題[10]。因此，采用訓(xùn)練耗時較小且較少出現(xiàn)過擬合現(xiàn)象的最小二乘支持向量機[11](Least Squares Support Vector Maehine，LS-SVM)判別系統(tǒng)方程計算信號Lyapunov指數(shù)譜的方法，分析軸承振動的超混沌特性。

1 軸承振動信號采集及相空間重構(gòu)

1.1 軸承振動信號的采集

采用Case West Reserve University滾動軸承實驗室實測的內(nèi)、外圈和滾動體的不同故障程度與轉(zhuǎn)速振動數(shù)據(jù)庫，該數(shù)據(jù)庫采用深溝球軸承6205-2RS JEM SKF，電動機空載，采樣頻率fs=12 kHz，轉(zhuǎn)速為1 797 r/min，軸承故障由電火花加工機在球軸承內(nèi)、外圈及鋼球上模擬損傷性故障，故障直徑分別為0.18，0.36，0.54和0.72 mm(文中選取內(nèi)圈模擬點蝕故障并在后續(xù)分析中分別對應(yīng)簡化表示為f1～f4)，故障深度為0.28 mm。正常狀態(tài)和內(nèi)圈故障振動信號的時域波形如圖1所示，并分別取1 024個點研究不同故障類型及不同程度故障軸承振動的混沌特性。

(b)內(nèi)圈故障

1.2 相空間重構(gòu)技術(shù)

相空間重構(gòu)技術(shù)是混沌特性研究的前提。混沌吸引子是混沌系統(tǒng)的特征之一，體現(xiàn)了混沌系統(tǒng)的運動規(guī)律，相空間重構(gòu)目的就是為了在高維空間中恢復(fù)混沌吸引子。系統(tǒng)運動是有n個變量的動力系統(tǒng)，即一組n個變量的一階微分方程

(1)

由連續(xù)變量坐標(biāo)x(t)及其(n-1)階導(dǎo)數(shù)構(gòu)成的n維相空間就表示了系統(tǒng)隨時間演變的狀態(tài)空間，即

X(t)=[x(t),x1(t),…,xn-1(t)]T。

(2)

Ruelle用離散序列x(t)及其(n-1)時滯位移將其代替，構(gòu)成一個新的n維嵌入相空間，即

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

,(3)

式中：τ為時延；m為重構(gòu)相空間的維數(shù)，m≥2d+1(d為狀態(tài)空間的關(guān)聯(lián)維數(shù)，也稱拓撲維)。對xn(n=1,2,…,N)，取定τ，其重構(gòu)相空間為

Xl= [xl,xl+τ,…,xl+(m-1)τ]T=

(4)

式中：l=N-(m-1)τ。求取τ使原序列時延后能作為獨立坐標(biāo)，可根據(jù)平均位移[14]、自相關(guān)[13]、互信息[12]及(去偏)復(fù)自相關(guān)技術(shù)等計算。好的重構(gòu)相空間是使重構(gòu)后的吸引子與系統(tǒng)真正的吸引子拓撲等價。m是嵌入維，可根據(jù)預(yù)測誤差最小[16]、幾何不變量[14]、虛假鄰近點[17]及其改進后的Cao方法[18]計算。確定時延與嵌入維有2種爭議：一種認為可先求得時延再選擇嵌入維；另一種則認為時延與嵌入維是相關(guān)的。研究者提出較常用的C-C[19]及其改進法[20-21]能同時計算時延與嵌入窗。因為不能得到混沌系統(tǒng)先驗知識，這些方法均有一定的主觀性[22]。

分析中采用自相關(guān)法[13]求τ，根據(jù)xn與xn+τ有些獨立又不完全不相關(guān)的特性，為使其能在重構(gòu)相空間中作為獨立坐標(biāo)處理，計算自相關(guān)函數(shù)為

(5)

則C(τ)首次下降到初值的1-1/e時的延時為最優(yōu)τ，不同嵌入維下正常及內(nèi)圈故障振動信號的自相關(guān)函數(shù)分別如圖2所示，從圖中可以看出，正常狀態(tài)最佳延時為3，不同內(nèi)圈故障狀態(tài)的最佳延時均為2。

(b)內(nèi)圈故障

再用虛假鄰近點法[17]計算m，先選定m，求得Xl的最近鄰點Xη(l)，即

‖Xl-Xη(l)‖m=min{‖Xl-Xj‖:j=1,2,…,l,l≠j}，

(6)

將m增加到m+1，當(dāng)滿足

(7)

Xη(l)即為Xl的虛假鄰近點。(7)式說明當(dāng)相空間從m維演變到m+1維時，2個相點相差很大距離是因為高維混沌吸引子中2個不相鄰的點投影到低維軌道上時變成相鄰的2個臨近虛假點造成的。求取最佳嵌入維數(shù)m時一般取RT=10[17]，統(tǒng)計每個m時其虛假鄰近點占所有重構(gòu)向量數(shù)的比例P(m)，直到P(m)不再隨m的增加而減小或小于某個值(如l%或5%)時，吸引子被認為完全展開。

在不同嵌入維下正常及故障狀態(tài)統(tǒng)計的P(m)如圖3所示。從圖中可以看出，當(dāng)m=5和m=4時，正常狀態(tài)及內(nèi)圈故障狀態(tài)的P(m)不再明顯減少，因此正常狀態(tài)最佳嵌入維取5，故障狀態(tài)最佳嵌入維均取4。

(a)正常狀態(tài)

根據(jù)嵌入維和延時，采用相空間重構(gòu)技術(shù)繪制正常和內(nèi)圈故障振動信號的混沌吸引子如圖4所示。從圖中可以看出，混沌吸引子在有限的空間不斷的纏繞，呈現(xiàn)一定自相似的狀態(tài)。

圖4 滾動軸承振動信號混沌吸引子

2 滾動軸承振動的混沌特性分析

首先采用遞歸圖對實測軸承不同振動狀態(tài)的相關(guān)性和確定性進行檢驗，然后使用CLY方法判定軸承振動的混沌特性，最后結(jié)合功率譜對振動進行頻譜分析。

2.1 遞歸圖

為分析序列周期及非周期成分的相對大小，揭示動力系統(tǒng)內(nèi)部相似性的先驗知識，文獻[23]提出了遞歸圖(Recurrence Plot，RP)方法。周期系統(tǒng)信號被噪聲污染時可能表現(xiàn)出非周期性，這表明系統(tǒng)不會準(zhǔn)確地回到以前的狀態(tài)變量，但可能返回很接近以前某時刻的狀態(tài)，因此可用來表現(xiàn)周期及非周期的程度。由(4)式采用嵌入技術(shù)把信號重構(gòu)為向量Xl(l=1,2,…,N-(m-1)τ)。然后求Xl的任意2個時刻i,j(i,j=1,2,…,N-(m-1)τ)的距離，即

dij=‖Xi-Xj‖ 。

(8)

如果dij≤r，以(i,j)為縱坐標(biāo)做一點，即繪制出遞歸圖，表現(xiàn)重構(gòu)的軌線重復(fù)或遞歸其自身的信息，表示了系統(tǒng)的時間關(guān)聯(lián)情況。試驗中，計算正常狀態(tài)及內(nèi)圈故障狀態(tài)信號時分別取r=2和r=1，得到的遞歸圖如圖5所示。

圖5 滾動軸承振動信號的遞歸圖

與文獻[24]中周期序列、混沌序列及隨機序列的遞歸圖進行對比分析發(fā)現(xiàn)：軸承正常狀態(tài)的遞歸圖存在較明顯的平行于對角線的線段，但帶狀線不是連續(xù)的，說明正常振動為非嚴格周期狀態(tài)，存在隨機性因素(部分為噪聲干擾)；內(nèi)圈故障狀態(tài)遞歸圖也存在較明顯的平行于對角線的線段，與混沌系統(tǒng)類似，表明其也具有周期成分，但介于周期與隨機之間，振動為混沌特征，且故障狀態(tài)越大，遞歸圖與白噪聲信號越相似；但軸承正常和內(nèi)圈故障狀態(tài)周期性均不明顯，由于繪制遞歸圖受噪聲和r的影響，其周期和非周期性大小只可作為定性的判斷，還需要參考其他方法進行綜合的分析判定。

2.2 CLY方法

為判別信號的混沌特性，以虛假鄰域方法為基礎(chǔ)，文獻[25]提出了CLY算法：由(4)式將xn(n=1,2,…,N)嵌入到m維的相空間，對選定的m，由 (6) 式計算每個Xl(m)(l=1,2,…,N-(m-1)τ)的最近鄰點Xη(l)(m)，得到平均預(yù)測誤差

。(9)

定義函數(shù)E(m)=E*(m+1)/E*(m)。對選擇的任意m，通常隨機信號的E(m)=1或在其上下波動；而對混沌系統(tǒng)，隨m的增大穩(wěn)定在1附近。故根據(jù)m-E(m)曲線可簡單的判別信號的混沌特性。正常及內(nèi)圈故障狀態(tài)振動信號的m-E(m)曲線如圖6所示，從圖中可以看出，當(dāng)m≥5時E(m)趨于1，可簡單判別出不同振動狀態(tài)都存在混沌特性。

2.3 功率譜

對xn(n=1,2,…,N)加上周期條件xn=xn+j，求其自相關(guān)函數(shù)(離散卷積)

(10)

然后對Cn進行離散Fourier變換，計算Fourier變換系數(shù)[26]

(11)

圖6 滾動軸承振動信號的m-E(m)曲線

圖7 不同信號的功率譜

正常和內(nèi)圈故障狀態(tài)振動信號的功率譜如圖8所示，從圖中可以看出：正常狀態(tài)振動信號有較明顯的獨立尖峰譜線，表明出較強的周期性，而存在連續(xù)的譜線，又表明其也具有混沌性，部分譜線分布類似于隨機信號譜線，則對應(yīng)于其中由噪聲引起的隨機成分；內(nèi)圈故障狀態(tài)振動信號的功率譜類似于Lorenz序列的功率譜圖，其功率譜線不同于周期和隨機信號的離散譜，表現(xiàn)為寬帶連續(xù)譜，還表現(xiàn)出在高頻段隨頻率指數(shù)衰減的特性，因此故障狀態(tài)存在很強的混沌性，但周期性并不明顯。綜合分析可知，軸承不同狀態(tài)振動信號的頻率分布特性存在顯著的差異，正常狀態(tài)具有明顯的周期性，而內(nèi)圈故障狀態(tài)具有明顯的混沌特性，與遞歸圖判定結(jié)論基本一致。

圖8 滾動軸承振動信號的功率譜

3 滾動軸承振動的超混沌特性分析

Lyapunov指數(shù)譜λi(i=1,2,…,m-1)是表示相空間各方向軌道旁的膨脹及收縮的平均度量，由系統(tǒng)長時間演變決定，反映了初始時刻2個無限靠近的點隨時間分離的情況，每個λi都是相空間不同方向上相對運動局部變形的平均，不同特征的運動與λi的關(guān)系見表1。

混沌吸引子的局部不穩(wěn)定和某方向上的指數(shù)發(fā)散或膨脹相對應(yīng)，有2個或2個以上正的λi系統(tǒng)就為超混沌，否則系統(tǒng)為非混沌?，F(xiàn)采用LS-SVM判別系統(tǒng)方程方法得到軸承振動的Lyapunov指數(shù)譜[27-28]，研究滾動軸承不同狀態(tài)振動的超混沌特性。

表1 不同運動與Lyapunov指數(shù)的關(guān)系

3.1 LS-SVM算法

由(4)式重構(gòu)xn(n=1,2,…,N)在m維歐式空間的一條軌道Xl(l=1,2,…,N-(m-1)τ)，則Xl能再現(xiàn)原系統(tǒng)的動力學(xué)特性，設(shè)Xl到Xl+τ的映射關(guān)系為G[6]，即

(12)

顯然Xi+τ的前m-1個數(shù)值即為Xi的后m-1個數(shù)值，只要確定Xl→Xl+mτ的關(guān)系F即可確定G，設(shè)

xl+mτ=F(Xl) 。

(13)

1)LS-SVM辨識未知系統(tǒng)

采用LS-SVM辨識F。先由原始信號xn(n=1,2,…,N)構(gòu)造訓(xùn)練數(shù)據(jù)集(Xl,xj+τ)，其中l(wèi)=1,2,…,N-mτ。利用高維特征空間的線性函數(shù)擬合樣本集[29]

xl+mτ=ωTφ(Xl)+b，

(14)

則原空間中的約束條件和優(yōu)化問題可表達為

s.t.xj+mτ=ωTφ(Xl)+b+el，

(15)

式中：el為誤差，其組成的誤差向量為e，偏差為b，權(quán)重為C。用Lagrange乘子al有

b+ej-xj+mτ}。

(16)

由KKT條件并消除el和ω，得到線性方程組

(17)

式中：Y=[x1+mτ,x2+mτ,…,xN]T；α=[α1,α2,…,αN-mτ]T；I=[1,1,…,1]T；E為N-mτ階單位陣；l=(ki,j)(N-mτ)×(N-mτ)，采用徑向基函數(shù)得到K的每個元素

(18)

解方程(17)得到α和b，則LS-SVM回歸函數(shù)逼近 (12) 式中的F為

(19)

2)計算Lyapunov指數(shù)譜

由 (12) 式得到

ΔXl+τ=DG(Xl)·ΔXl，

(20)

(21)

式中：DG(Xn)為映射G的Jacobian矩陣；ΔXl表示重構(gòu)的向量Xl在l時刻的微小變化。通過DG(Xl)矩陣，在l+τ時刻重構(gòu)相空間的向量Xl+τ值的微小變化中將被表現(xiàn)出來。繼續(xù)下去，該作用會累積到l+rτ(r∈N)時刻重構(gòu)相空間的向量Xl+rτ取值的變化，其關(guān)系為

ΔXl+rτ=DG(Xl+(r-1)τ)…DG(Xl)·ΔXl。(22)

(23)

按大小排列λi(i=1,2,…,m)。因為 (20) 式中矩陣會產(chǎn)生分數(shù)冪和指數(shù)，所以，ΔXl+τ往往是病態(tài)陣，不能求得其全部精確特征值。可利用長乘積矩陣的分解方法對 (23) 式求解。首先定義

Ar=DG(Xl+(r-1)τ)…DG(Xl)=A(r)A(r-1)…A(1)

，(24)

對此長乘積矩陣，遞歸計算

A(i)Q(i-1)=Q(i)R(i),i=1,2…r，

(25)

式中：Q(i)為正交矩陣；R(i)為上三角矩陣；Q(0)為m階單位陣。按 (25) 式方法將QR分解r次，則Ar的QR分解為Ar=Q(r)R(r)R(r-1) …R(1)，故 (22) 式的本征值可以通過下式求出

(26)

進而求得系統(tǒng)的Lyapunov指數(shù)譜為

(27)

3.2 計算結(jié)果分析

采用前文給出的實測滾動軸承正常和內(nèi)圈故障狀態(tài)的振動信號，分析其振動的超混沌特性。計算同樣采用相關(guān)法和虛假鄰近點法分別確定延時τ和嵌入維m。LS-SVM計算參數(shù)取為：C=1 000，a=1，計算出的Lyapunov指數(shù)譜見表2。

表2 不同狀態(tài)振動的Lyapunov指數(shù)譜

由表可知，軸承正常及內(nèi)圈故障狀態(tài)振動均為超混沌，但正常狀態(tài)振動的超混沌特征不明顯。當(dāng)故障增大時，軸承超混沌特性程度也增強，并且表現(xiàn)的更加明顯。

此外，對軸承外圈、滾動體及保持架不同故障狀態(tài)進行了相同的混沌及超混沌試驗分析，均得到同樣的結(jié)果。因此軸承振動的混沌程度與故障大小存在一定的對應(yīng)關(guān)系，正常狀態(tài)振動比較平穩(wěn)有序，較大故障的振動狀態(tài)會增加混亂程度。

4 結(jié)束語

應(yīng)用4種方法，從多個方面對實測軸承正常及故障狀態(tài)振動的非線性混沌特性進行了判定和分析，研究表明：軸承正常及故障狀態(tài)振動的混沌特征具有明顯的差異性，研究軸承振動的非線性混沌特性要結(jié)合不同的振動狀態(tài)。上述研究結(jié)果有助于理解軸承振動的復(fù)雜演化規(guī)律，為今后使用其振動信號混沌特性提取故障特征及進行故障診斷奠定一定的基礎(chǔ)。