周博，王瑤瑤，劉雙燕，文振華

(鄭州航空工業(yè)管理學院 a.民航學院；b.航空工程學院，鄭州 450046)

滾動軸承的工作狀況對旋轉(zhuǎn)設(shè)備的運轉(zhuǎn)起著至關(guān)重要的作用[1]，相關(guān)調(diào)查研究發(fā)現(xiàn)，滾動軸承引起的問題占旋轉(zhuǎn)設(shè)備故障的40%以上，因此，對滾動軸承進行故障診斷與健康評估的研究具有重要意義[2]，其關(guān)鍵是對傳感器所獲取信號進行處理并提取有效的故障特征。實際工作中，由于機械部件的負載、摩擦力、阻尼不斷變化，導致產(chǎn)生的振動信號常呈現(xiàn)出明顯的非線性非穩(wěn)態(tài)特點[3]。傳統(tǒng)基于傅里葉變換的方法只能處理線性和平穩(wěn)信號，因此對非平穩(wěn)信號的特征提取成為機械故障診斷領(lǐng)域的重要課題[4]。

國內(nèi)外研究人員針對非穩(wěn)態(tài)信號的故障特征提取做了許多研究，傅里葉變換、Wigner-Ville分布、小波變換、希爾伯特-黃變換(HHT)、局部均值分解(LMD)等時頻分析方法可以同時得到振動信號的時域信息和頻域信息，近年來在旋轉(zhuǎn)機械的故障診斷中應用廣泛，且診斷效果良好。然而，上述分析方法在使用中均存在一定的不足：短時傅里葉變換的時頻窗口是固定的[5]；小波變換母小波的選擇和分解層數(shù)不具備自適應性[6]；HHT實現(xiàn)了信號的自適應分解，但存在過包絡、欠包絡以及擬合精度差等問題[7]；LMD雖然能夠改善上述缺陷[8]，但其得到的分量往往摻雜著虛假頻率，處理過程中還可能會發(fā)生信號突變，而且該方法計算量大，運算速度慢。局部特征尺度分解(Local Characteristic-Scale Decomposition，LCD)具有較強的自適應性，能夠?qū)⒎瞧椒€(wěn)信號分解成為多個內(nèi)稟尺度分量(Intrinsic Scale Component，ISC)之和[9]，相對于前述時頻分析方法，LCD對信號的處理速度有很大提升，同時避免了過包絡和欠包絡的問題，而且降低了擬合誤差。

另外，在滾動軸承故障診斷和健康評估的過程中，提取有效的故障特征非常關(guān)鍵，同時也是信號分解的重要內(nèi)容。遞歸圖是基于相空間重構(gòu)理論分析時間序列的周期性、混沌性以及非平穩(wěn)性的重要方法，可以獲得信息量、相似性以及預測性等相關(guān)內(nèi)容，了解時間序列的內(nèi)在構(gòu)成[10-12]。然而遞歸圖只是一種定性分析方法，無定量判斷依據(jù)，主要依靠科研人員的經(jīng)驗進行結(jié)果判定，客觀性不足。因此，嘗試通過提取遞歸圖中的一些定量特性對非線性時間序列進行分析[13-14]，通過遞歸定量分析(Recurrence Quantification Analysis,RQA)對滾動軸承進行更準確的故障診斷與健康評估[15]。

綜上，將LCD和RQA進行融合，提出了基于LCD-RQA的滾動軸承故障診斷與健康評估方法：首先，利用LCD對傳感器獲取的滾動軸承振動信號進行自適應分解；然后，對分解的每條曲線進行相空間重構(gòu)并繪制成為遞歸圖，提取遞歸圖中的遞歸率、確定性和層流性等作為故障特征；最后，利用提取出來的遞歸定量特征進行滾動軸承的故障診斷與健康評估。

1 局部特征尺度分解

局部特征尺度分解適用于處理非線性和非穩(wěn)態(tài)信號，該過程將得到n個從高頻到低頻的內(nèi)稟尺度分量[16]，ISC分量必須滿足以下條件：

①在整個數(shù)據(jù)集上，任何相鄰2個極值點的符號互不相同。

②在整個數(shù)據(jù)集上，對于局部最大(最小)值點(tk,xk)，任意臨近的最大(最小)值(tk-1,xk-1)與(tk+1,xk+1)由一條直線相連，該直線可表示為

(1)

為保證ISC曲線的光滑性和對稱性，xk與Ak的比值為一個常數(shù)，即

Ak/xk=(a-1)/a；a∈(0,1)。

(2)

通常情況下，a取0.5，因此Ak=-xk，此時xk與Ak關(guān)于x軸對稱。

對于一個復雜的信號x(t)(t>0)，可按如下步驟進行分解得到ISC分量：

1)選取信號x(t)的極值點(tk,xk)。

2)利用(1)式計算Ak并計算對應的Lk,即

Lk=aAk+(1-a)xk；k=2,…,M-1，

(3)

由于k的取值范圍為從[2，M-1]，因此需要將數(shù)據(jù)邊界進行延長，設(shè)A1(t0,x0)和AM(tM+1,xM+1)為延伸之后兩端的極值點，得到L1和LM。

3)將所有的Lk(k=1,…,M)用三次樣條曲線L(t)連接，該曲線即LCD的基線。設(shè)原始信號與基線之間的差值h1(t)為

h1(t)=x(t)-L1(t)，

(4)

若h1(t)滿足條件①和②，則h1(t)作為第1個ISC分量；否則將h1(t)視為原始信號重復以上步驟，得

h11(t)=h1(t)-L11(t)，

(5)

若h11(t)依舊不滿足條件①和②，重復以上步驟k次，直到h1k(t)滿足條件成為第1個ISC分量。

4)將第1個ISC分量c1(t)從原始數(shù)據(jù)x(t)中分離出來，殘差記作r1(t)，即

x(t)-c1(t)=r1(t)。

(6)

5)繼續(xù)將殘差r1(t)作為原始信號進行處理，重復以上步驟直到最終的殘差rn(t)為常數(shù)或單調(diào)函數(shù)，以及極值點不超過3個的函數(shù)。則原始信號x(t)被分解成為多個ISC分量與殘差之和，即

(7)

式中：ci(t)為第i個ISC分量。

2 遞歸定量分析

2.1 遞歸圖的生成

遞歸圖的構(gòu)造原理如下[18]：

1)對于時間序列uk(k=1,2,…,N)，其采樣間隔為Δt，假設(shè)嵌入維數(shù)為m，延遲時間為τ，通過文獻[19]和互信息法選取恰當?shù)膍和τ，得到一個新的時間序列，重構(gòu)之后的相空間為

xi=(ui,ui+τ,…,ui+(m-1)τ)；i=1,2,…,N-(m-1)τ。

(8)

2)計算重構(gòu)后相空間中的點xi與點xj之間的距離，即

(9)

3)計算遞歸值。

R(i,j)=H(εi-Sij)，

(10)

(11)

式中:εi為截止距離；H(r)為Heaviside單位函數(shù)。

4)以i為橫坐標，j為縱坐標繪制遞歸圖R(i,j)，i，j分別為時間序列標號。

2.1.1 文獻[19]方法

文獻[19]在進行信號處理時可以清晰地將真實信號和噪聲分辨出來，且計算效率較高。

在不同嵌入維數(shù)條件下，相空間中各點的最鄰近點的距離變化值為

(12)

定義2個判別準則，即

(13)

(14)

當時間序列固定不變時存在嵌入維數(shù)，若d>d0時，E1(d)和E2(d)停止變化或者波動較小，則最小嵌入維數(shù)為d0。

2.1.2 互信息法

互信息函數(shù)方法在相空間重構(gòu)中運用較多，用于估計重構(gòu)相空間的延遲時間[20]。

對于離散變量X，Y，其聯(lián)合熵為

(15)

式中：pij為變量X在狀態(tài)i且變量Y在狀態(tài)j時出現(xiàn)的概率；m，n分別為X，Y的狀態(tài)數(shù)。

互信息可以定義為

I(X,Y)=H(X)+H(Y)-H(X,Y)。

(16)

設(shè)混沌時間序列s1,s2,…，時間延遲為d，嵌入維數(shù)為τ，重構(gòu)相空間得

X(t)={x0(t),x1(t),…,xn(t)}，

(17)

xn(t)=s(t+nτ)，

(18)

則系統(tǒng)對變量s的平均信息量為系統(tǒng)的熵H，即

(19)

記[s,q]=[x(t),x(t+τ)]，考慮一個總的耦合系統(tǒng)(S,Q)，假定s已知為si，則q的不定性為

(20)

式中：Pq|s(qj|si)為條件概率。

設(shè)在時刻t時x已知，則在x+τ時刻，x的平均不定性為

H(S,Q)-H(S)，

(21)

(22)

式中：H(S,Q)為孤立的q的不定性；H(Q|S)為已知s的q的不定性。所以s的已知減小q的不定性，互信息為

I(Q,S)=H(Q)-H(Q|S)=H(Q)+

H(S)-H(S,Q)=I(S,Q)。

(23)

互信息是聯(lián)合概率分布Psq的函數(shù)。從概率分布的直方圖中估計Psq通常采用以下方法：設(shè)在平面上點(s,q)處有一個大小為ΔsΔq的盒子，那么有

(24)

式中：Nsq，Ntotal分別為盒子中的點數(shù)和總點數(shù)。

對于一般情況，互信息為

(25)

若In(τ)是一個延遲時間重構(gòu)，相空間重構(gòu)的時間延遲則為In(τ)初次達到最小值時的時滯。

2.2 遞歸定量特征的提取

定義遞歸圖中線結(jié)構(gòu)和點密度中的特征遞歸率、確定性和層流性等參數(shù)來定量地描述時間序列的動態(tài)特性[3,13]。

對于相空間重構(gòu)后的N個向量，根據(jù)(10)式構(gòu)建遞歸圖后，遞歸率RRR的定義為

(26)

設(shè)P(l)，P(v)分別為遞歸圖中45°方向和垂直方向直線的長度分布，分別定義為

(27)

(28)

式中：Nl為長度l的45°方向直線的數(shù)量；Nv為長度v的垂直方向直線的數(shù)量；Nα為長度α的45°方向或垂直方向直線的數(shù)量；lmin，lmax為45°方向直線的最小長度(一般取2)和最大長度；vmin，vmax為垂直方向直線的最小長度(一般取2)和最大長度。

確定性RDET和層流性RLAM的定義為

(29)

(30)

遞歸率、確定性和層流性分別反映了某個特定狀態(tài)出現(xiàn)的概率以及系統(tǒng)的可預測性、間歇性和層次性，這些遞歸圖中提取的參數(shù)都是系統(tǒng)動力學特征的反映，可以作為特征參數(shù)用于滾動軸承的故障診斷與健康評估。

3 基于LCD-RQA的滾動軸承故障診斷與健康評估

3.1 算法流程

綜合上述理論分析，提出的基于LCD-RQA的滾動軸承故障診斷與健康評估方法的流程如圖1所示。

圖1 基于LCD-RQA的滾動軸承故障診斷與健康評估流程

3.2 數(shù)據(jù)準備

以美國Case Western大學所公開的滾動軸承故障診斷試驗數(shù)據(jù)[21]為例進行方法驗證。

故障診斷案例以驅(qū)動端SKF 6205-2RS型深溝球軸承為研究對象，分別選取1 797 r/min下正常狀態(tài)以及故障直徑為0.534 mm的內(nèi)圈故障、鋼球故障和外圈故障軸承所采集的數(shù)據(jù)，具體數(shù)據(jù)組成見表1，利用該數(shù)據(jù)驗證LCD-RQA滾動軸承故障診斷方法的可行性。

表1 試驗數(shù)據(jù)描述

健康評估的案例同樣來自驅(qū)動端軸承SKF 6205-2RS型深溝球軸承，分別選取1 797 r/min下正常狀態(tài)、內(nèi)圈滾道上故障直徑分別為0.178，0.356，0.534 mm的試驗數(shù)據(jù)，具體數(shù)據(jù)組成見表1，利用該數(shù)據(jù)集驗證LCD-RQA滾動軸承健康評估方法的可行性。

3.3 基于LCD-RQA的滾動軸承故障診斷

首先利用文獻[19]算法和互信息法計算相空間重構(gòu)的嵌入維數(shù)m和延遲時間τ，計算最小嵌入維數(shù)時，以E(d+1)-E(d)<0.01和E*(d+1)-E*(d)<0.01作為(13)式和(14)式的定量判定標準，具體結(jié)果見表2。然后利用LCD對4種故障狀態(tài)的軸承振動數(shù)據(jù)進行自適應分解，分解層數(shù)為5，分解后的各ISC分量如圖2所示。

圖2 不同故障狀態(tài)軸承振動數(shù)據(jù)的LCD處理結(jié)果

表2 不同故障狀態(tài)軸承數(shù)據(jù)相空間重構(gòu)參數(shù)

將獲得的ISC分量轉(zhuǎn)化成為遞歸圖，轉(zhuǎn)化時每次選取的ISC分量長度為1 000個數(shù)據(jù)點，且每次向后滑移100個數(shù)據(jù)點來形成多個遞歸圖。不同故障狀態(tài)軸承振動數(shù)據(jù)各ISC分量第1～1 000個數(shù)據(jù)點所轉(zhuǎn)化的遞歸圖如圖3所示，由圖可知，不同故障狀態(tài)軸承振動數(shù)據(jù)的遞歸圖表現(xiàn)出遞歸點的密度和水平垂直線的結(jié)構(gòu)均不相同。

圖3 不同故障狀態(tài)軸承振動數(shù)據(jù)的遞歸圖

為定量描述各故障狀態(tài)軸承振動數(shù)據(jù)的遞歸圖，分別計算各遞歸圖的遞歸率、確定性和層流性并將計算結(jié)果繪制在三維散點圖中，結(jié)果如圖4所示，由圖可知：軸承振動數(shù)據(jù)的遞歸特征具有很好的可分性，同一故障狀態(tài)的遞歸定量特征具有很好的聚集性，說明本方法可以很好地對滾動軸承進行故障診斷。

圖4 不同故障狀態(tài)軸承振動數(shù)據(jù)遞歸特征的三維散點圖

3.4 基于LCD-RQA的滾動軸承健康評估

同樣，利用文獻[19]算法和互信息法計算相空間重構(gòu)的嵌入維數(shù)m和延遲時間τ，計算最小嵌入維數(shù)時，以E(d+1)-E(d)<0.01和E*(d+1)-E*(d)<0.01作為(13)式和(14)式的定量判定標準，具體結(jié)果見表3。然后利用LCD對正常狀態(tài)及不同內(nèi)圈故障程度的軸承進行自適應分解，分解層數(shù)為5時所得各ISC分量如圖5所示。

表3 不同故障程度軸承振動數(shù)據(jù)的相空間重構(gòu)參數(shù)

圖5 不同內(nèi)圈故障程度軸承振動數(shù)據(jù)的LCD處理結(jié)果

將獲得的ISC分量轉(zhuǎn)化成為遞歸圖，轉(zhuǎn)化時每次選取的ISC分量長度為1 000個數(shù)據(jù)點，且每次向后滑移100個數(shù)據(jù)點來形成多個遞歸圖。不同內(nèi)圈故障程度軸承振動數(shù)據(jù)各ISC分量第1～1 000個數(shù)據(jù)點所轉(zhuǎn)化的遞歸圖如圖6所示，由圖可知遞歸點的密度和水平垂直線的結(jié)構(gòu)隨著軸承故障程度的增加而改變[3]。

圖6 不同內(nèi)圈故障程度軸承振動數(shù)據(jù)的遞歸圖

正常狀態(tài)循環(huán)計算300組，3種程度內(nèi)圈故障數(shù)據(jù)均循環(huán)計算200組，生成遞歸率、確定性和層流性這3類遞歸定量分析特征參數(shù)。將正常狀態(tài)的前100組特征作為健康評估的標準數(shù)據(jù)，其余數(shù)據(jù)的故障特征作為測試數(shù)據(jù)，計算測試數(shù)據(jù)遞歸定量特征與正常數(shù)據(jù)遞歸定量特征之間的歐式距離，并將所計算的歐式距離轉(zhuǎn)化成為置信度(CV)以表示滾動軸承的健康程度，即

(31)

式中：di為歐式距離；a為計算標準化參數(shù)。

正常數(shù)據(jù)與測試數(shù)據(jù)之間的歐式距離及CV值如圖7所示，由圖可知：測試數(shù)據(jù)正常樣本與標準數(shù)據(jù)之間的歐式距離非常小，在0附近波動，隨著故障程度的加深，各故障樣本與標準數(shù)據(jù)之間的歐式距離逐漸增大；將歐式距離曲線轉(zhuǎn)換成為CV值并歸一化到[0，1]，正常情況的CV值在1附近，隨著故障程度的加深，該滾動軸承的健康度在不斷下降。

圖7 不同故障程度軸承振動數(shù)據(jù)的歐式距離和健康度

根據(jù)經(jīng)驗將正常閾值設(shè)定為0.8，完全故障的閾值設(shè)定為0.4，根據(jù)故障注入的要求，故障程度1的滾動軸承故障程度較輕，該狀態(tài)下的健康度可接受，軸承可繼續(xù)使用；故障程度2的滾動軸承健康狀態(tài)已經(jīng)接近完全故障，而故障程度3的滾動軸承健康狀態(tài)低于設(shè)定的完全故障閾值，判斷為不可使用?？筛鶕?jù)該參數(shù)設(shè)定下的模型對后續(xù)設(shè)備進行健康評估。

4 結(jié)束語

針對滾動軸承振動信號非線性非平穩(wěn)的特點，提出了一種基于局部特征尺度分解及遞歸定量分析的滾動軸承故障診斷和健康評估方法。通過局部特征尺度分解對滾動軸承振動信號進行分解，將分解后的各ISC分量轉(zhuǎn)化為遞歸圖并計算其定量特征。試驗數(shù)據(jù)的三維散點圖可看出各故障狀態(tài)具有很好的可分性，健康評估中測試數(shù)據(jù)與標準數(shù)據(jù)的歐式距離及置信度則表明隨著故障程度的加深，滾動軸承健康度在逐漸下降。該方法對非線性振動信號具有良好的適用性，對滾動軸承的故障診斷和健康評估應用效果顯著。