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

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

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

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

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

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

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

1 局部特征尺度分解

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

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

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

(1)

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

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

(2)

通常情況下，a取0.5，因此Ak=-xk，此時xk與Ak關于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]，因此需要將數據邊界進行延長，設A1(t0,x0)和AM(tM+1,xM+1)為延伸之后兩端的極值點，得到L1和LM。

3)將所有的Lk(k=1,…,M)用三次樣條曲線L(t)連接，該曲線即LCD的基線。設原始信號與基線之間的差值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)從原始數據x(t)中分離出來，殘差記作r1(t)，即

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

(6)

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

(7)

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

2 遞歸定量分析

2.1 遞歸圖的生成

遞歸圖的構造原理如下[18]：

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

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

(8)

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

(9)

3)計算遞歸值。

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

(10)

(11)

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

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

2.1.1 文獻[19]方法

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

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

(12)

定義2個判別準則，即

(13)

(14)

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

2.1.2 互信息法

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

對于離散變量X，Y，其聯合熵為

(15)

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

互信息可以定義為

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

(16)

設混沌時間序列s1,s2,…，時間延遲為d，嵌入維數為τ，重構相空間得

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

(17)

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

(18)

則系統對變量s的平均信息量為系統的熵H，即

(19)

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

(20)

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

設在時刻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)

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

(24)

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

對于一般情況，互信息為

(25)

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

2.2 遞歸定量特征的提取

定義遞歸圖中線結構和點密度中的特征遞歸率、確定性和層流性等參數來定量地描述時間序列的動態特性[3,13]。

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

(26)

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

(27)

(28)

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

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

(29)

(30)

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

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

3.1 算法流程

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

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

3.2 數據準備

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

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

表1 試驗數據描述

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

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

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

圖2 不同故障狀態軸承振動數據的LCD處理結果

表2 不同故障狀態軸承數據相空間重構參數

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

圖3 不同故障狀態軸承振動數據的遞歸圖

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

圖4 不同故障狀態軸承振動數據遞歸特征的三維散點圖

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

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

表3 不同故障程度軸承振動數據的相空間重構參數

圖5 不同內圈故障程度軸承振動數據的LCD處理結果

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

圖6 不同內圈故障程度軸承振動數據的遞歸圖

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

(31)

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

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

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

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

4 結束語

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