滾動(dòng)軸承振動(dòng)的非線性超混沌特性研究

李兆飛，任小洪，黃臣程

(1.四川理工學(xué)院 自動(dòng)化與電子信息學(xué)院，四川 自貢 643000；2.人工智能四川省重點(diǎn)實(shí)驗(yàn)室，四川 自貢 643000)

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

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

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

1.1 軸承振動(dòng)信號(hào)的采集

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

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

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

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

(1)

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

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

(2)

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

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

,(3)

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

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

(4)

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

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

(5)

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

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

再用虛假鄰近點(diǎn)法[17]計(jì)算m，先選定m，求得Xl的最近鄰點(diǎn)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的虛假鄰近點(diǎn)。(7)式說(shuō)明當(dāng)相空間從m維演變到m+1維時(shí)，2個(gè)相點(diǎn)相差很大距離是因?yàn)楦呔S混沌吸引子中2個(gè)不相鄰的點(diǎn)投影到低維軌道上時(shí)變成相鄰的2個(gè)臨近虛假點(diǎn)造成的。求取最佳嵌入維數(shù)m時(shí)一般取RT=10[17]，統(tǒng)計(jì)每個(gè)m時(shí)其虛假鄰近點(diǎn)占所有重構(gòu)向量數(shù)的比例P(m)，直到P(m)不再隨m的增加而減小或小于某個(gè)值(如l%或5%)時(shí)，吸引子被認(rèn)為完全展開(kāi)。

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

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

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

圖4 滾動(dòng)軸承振動(dòng)信號(hào)混沌吸引子

2 滾動(dòng)軸承振動(dòng)的混沌特性分析

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

2.1 遞歸圖

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

dij=‖Xi-Xj‖ 。

(8)

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

圖5 滾動(dòng)軸承振動(dòng)信號(hào)的遞歸圖

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

2.2 CLY方法

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

。(9)

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

2.3 功率譜

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

(10)

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

(11)

圖6 滾動(dòng)軸承振動(dòng)信號(hào)的m-E(m)曲線

圖7 不同信號(hào)的功率譜

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

圖8 滾動(dòng)軸承振動(dòng)信號(hào)的功率譜

3 滾動(dòng)軸承振動(dòng)的超混沌特性分析

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

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

表1 不同運(yùn)動(dòng)與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)的動(dòng)力學(xué)特性，設(shè)Xl到Xl+τ的映射關(guān)系為G[6]，即

(12)

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

xl+mτ=F(Xl) 。

(13)

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

采用LS-SVM辨識(shí)F。先由原始信號(hào)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)化問(wèn)題可表達(dá)為

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的每個(gè)元素

(18)

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

(19)

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

由 (12) 式得到

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

(20)

(21)

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

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

(23)

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

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

，(24)

對(duì)此長(zhǎng)乘積矩陣，遞歸計(jì)算

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) 式的本征值可以通過(guò)下式求出

(26)

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

(27)

3.2 計(jì)算結(jié)果分析

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

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

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

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

4 結(jié)束語(yǔ)

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