基于K-L散度和TEO的滾動(dòng)軸承故障頻率識(shí)別方法

李一平,葉海天,曹 立,趙 琦,徐新勝

(中國(guó)計(jì)量大學(xué) 質(zhì)量與安全工程學(xué)院,浙江 杭州 310018)

滾動(dòng)軸承是各類旋轉(zhuǎn)機(jī)械的關(guān)鍵部件之一,廣泛應(yīng)用于航空航天、交通工具及機(jī)械設(shè)備等領(lǐng)域。由于復(fù)雜的工作環(huán)境和高強(qiáng)度的運(yùn)行,滾動(dòng)軸承內(nèi)圈、外圈及滾動(dòng)體等部件容易損壞失效,影響機(jī)械系統(tǒng)的正常運(yùn)行[1]。滾動(dòng)軸承發(fā)生故障時(shí)通常會(huì)產(chǎn)生沖擊、噪聲等表象,由此產(chǎn)生的振動(dòng)信號(hào)是一種有效的故障分析依據(jù),被廣泛應(yīng)用于機(jī)械故障診斷中。通常,滾動(dòng)軸承產(chǎn)生的振動(dòng)信號(hào)中包含了大量的運(yùn)行狀態(tài)信息,但振動(dòng)信號(hào)在傳遞過程中易受到環(huán)境噪聲和其他部件的干擾,因此,也會(huì)造成一些微弱的故障信號(hào)被淹沒,結(jié)果往往收集到的振動(dòng)信號(hào)具有非平穩(wěn)、非線性、非高斯性和低信噪比的特點(diǎn)[2]。于是,消除振動(dòng)信號(hào)的噪聲并有效提取故障特征,對(duì)后續(xù)滾動(dòng)軸承的故障診斷及預(yù)測(cè)性維修工作具有重要的意義。

近年來,專家學(xué)者對(duì)信號(hào)降噪問題,提出了如小波降噪[3]、自適應(yīng)濾波[3]、經(jīng)驗(yàn)?zāi)B(tài)分解[4]、頻域特征提取[5]等方法,并在信號(hào)特征提取、故障診斷、壽命預(yù)測(cè)等領(lǐng)域得到了廣泛應(yīng)用。但自適應(yīng)濾波和小波降噪技術(shù)的降噪效果過分依賴濾波器性能的優(yōu)劣[3],經(jīng)驗(yàn)?zāi)B(tài)分解雖然能將信號(hào)分解為一系列固有模態(tài)函數(shù)(IMF),但存在邊際效應(yīng)和模態(tài)混疊現(xiàn)象[4];同樣,頻域特征提取技術(shù)也存在計(jì)算復(fù)雜且降噪效果取決于信號(hào)的幅頻信息等問題[5]。因此,這些方法在工程實(shí)際應(yīng)用中存在一定的局限性。奇異值分解(Singular Value Decomposition,SVD)作為一種信號(hào)消噪方法近年來被廣泛的應(yīng)用于信號(hào)處理和統(tǒng)計(jì)學(xué)等領(lǐng)域[6]。該方法通過重構(gòu)矩陣及奇異值分解的形式,有效克服了原始信號(hào)信息丟失、計(jì)算復(fù)雜的問題,并通過選取有效的矩陣階次,可以達(dá)到較好的信號(hào)降噪效果。對(duì)奇異值分解(SVD)方法而言,確定有效的矩陣階次是保證信號(hào)降噪效果的關(guān)鍵。對(duì)此,錢征文[6]等根據(jù)信號(hào)快速傅里葉變換結(jié)果中主頻率個(gè)數(shù)來確定有效階次,采用不同頻率成分的幾組信號(hào)對(duì)該方法進(jìn)行了驗(yàn)證,但不適用強(qiáng)噪聲、多頻率的信號(hào)。周福成[7]等針對(duì)頻率切片小波變換在強(qiáng)背景噪聲條件下,故障特征識(shí)別能力不足的缺點(diǎn),提出根據(jù)奇異值差分譜最大值原則,確定有效階次進(jìn)行降噪處理,繼而利用頻率切片小波對(duì)降噪信號(hào)進(jìn)行了分析與應(yīng)用,而奇異值差分譜法在噪聲較強(qiáng)時(shí)無法保留全部有用信息,因此,實(shí)際應(yīng)用效果也存在不足。沈科宇[8]等通過S-T變換將時(shí)域信號(hào)向時(shí)頻域轉(zhuǎn)換,對(duì)時(shí)頻信號(hào)矩陣進(jìn)行奇異值分解,根據(jù)奇異值大小對(duì)信號(hào)降噪并重新構(gòu)造時(shí)域信號(hào),但該方法主觀性較強(qiáng),降噪效果及應(yīng)用場(chǎng)景也不理想。楊文獻(xiàn)[9]等通過分析信號(hào)信噪比與奇異熵之間的聯(lián)系,提出根據(jù)奇異熵增量漸進(jìn)特性來確定信號(hào)奇異值分解有效階次的方法,并采用仿真信號(hào)進(jìn)行了驗(yàn)證,但無法還原真實(shí)信號(hào)的狀態(tài)和特點(diǎn)。于是,針對(duì)信號(hào)數(shù)據(jù)的實(shí)際特點(diǎn)及信號(hào)降噪的應(yīng)用場(chǎng)景和目標(biāo),確定奇異值分解的有效階次,對(duì)信號(hào)有效降噪尤為重要。

此外,Teager能量算子(Teager Energy Operator,TEO)在20世紀(jì)90年代提出,由于它能追蹤并計(jì)算信號(hào)的瞬時(shí)能量,因此引起國(guó)內(nèi)外學(xué)者的關(guān)注并得到廣泛應(yīng)用[10]。王天金[11]等針對(duì)滾動(dòng)軸承故障振動(dòng)信號(hào)中的瞬態(tài)沖擊特點(diǎn),提出了基于TEO的頻譜分析方法,通過瞬時(shí)Teager能量的Fourier頻譜識(shí)別軸承的故障頻率。李輝[12]等利用TEO計(jì)算經(jīng)驗(yàn)?zāi)B(tài)分解后的齒輪箱軸承振動(dòng)信號(hào)分量,識(shí)別齒輪箱軸承的故障部位和類型。向玲[13]等提出了基于變分模態(tài)分解和1.5維Teager能量譜的滾動(dòng)軸承故障特征提取方法。

本文提出一種綜合K-L散度和TEO的滾動(dòng)軸承故障信號(hào)降噪及其頻率提取方法。首先通過SVD分解原始滾動(dòng)軸承振動(dòng)信號(hào)得到若干信號(hào)分量,為了克服傳統(tǒng)奇異值分解中確定有效階次主觀性強(qiáng)且降噪效果不顯著的缺點(diǎn),利用K-L散度計(jì)算原始信號(hào)與各信號(hào)分量的相關(guān)性來確定奇異值分解的有效階次,根據(jù)計(jì)算所得K-L散度值的差異性去除噪聲信號(hào),有效抑制噪聲干擾,完成信號(hào)降噪和重構(gòu)。由于滾動(dòng)軸承發(fā)生故障時(shí)會(huì)產(chǎn)生瞬態(tài)沖擊,TEO能夠及時(shí)跟蹤信號(hào)波形的變化且時(shí)間分辨率高,對(duì)重構(gòu)信號(hào)進(jìn)行TEO解調(diào),因而能增強(qiáng)重構(gòu)信號(hào)中微弱的故障頻率,提高故障頻率識(shí)別的準(zhǔn)確性。

1 奇異值分解

1.1 奇異值分解原理

假設(shè)任意實(shí)矩陣Am×n,秩為r,則A的奇異值分解為

A=USVT。

(1)

式(1)中,Um×m、Vn×n為標(biāo)準(zhǔn)化正交矩陣,Sm×n為奇異值對(duì)角矩陣,其對(duì)角線元素為矩陣A的非零奇異值。

矩陣A的秩為r,從式(1)中除去A的零奇異值,于是,得到A奇異值分解的精簡(jiǎn)形式:

(2)

式(2)中,σi為矩陣的非零奇異值,ui、vi分別為Um×m、Vn×n的第i個(gè)列向量。

1.2 奇異值分解降噪原理

奇異值分解(SVD)降噪是依據(jù)信號(hào)的信息量與奇異值的對(duì)應(yīng)關(guān)系,奇異值越大,包含信息量越多,奇異值越小,包含信息量越少,將包含信息量少的信號(hào)分量去除,即去除噪聲信號(hào)。由式(1)可知,將奇異值按降序排列,可認(rèn)為有效信號(hào)主要由數(shù)值較大的奇異值表示。通過設(shè)定閾值,將小于閾值的奇異值設(shè)為0,即可去除噪聲分量,保留原始信號(hào)的主要信息,從而完成信號(hào)降噪。

SVD降噪的步驟如下。

1) 一維信號(hào)轉(zhuǎn)化為矩陣

將一維信號(hào)進(jìn)行SVD分解的前提是將信號(hào)序列構(gòu)造為Hankel矩陣。對(duì)于一維振動(dòng)信號(hào)序列X(L)={x1,x2,…,xn},構(gòu)造Hankel矩陣如下:

(3)

式(3)中,A為m×n階Hankel矩陣A的行數(shù)和列數(shù)滿足L=m+n-1。

2) 對(duì)矩陣進(jìn)行奇異值分解

A=USVT。

(4)

3) 有效階次確定

奇異值分解以后,第個(gè)非零奇異值構(gòu)成的信號(hào)分量Xi=uieivi。通常,較大的奇異值構(gòu)成的信號(hào)分量包含原始信號(hào)更多的信息,通過設(shè)定閾值確定保留的奇異值個(gè)數(shù),將其余奇異值置為零,即確定奇異值分解的有效階次。

4) 重構(gòu)信號(hào)

在確定有效階次后,利用式(2)計(jì)算降噪后的信號(hào)矩陣X′,對(duì)X′的負(fù)對(duì)角線元素求均值,得到重構(gòu)信號(hào)x(i)=(x1,x2,…,xn),如式(5)所示:

(5)

式(5)中,α=max(1,i-m+1),β=min(n,1)。

重構(gòu)信號(hào)在保留原始振動(dòng)信號(hào)有用信息的同時(shí)進(jìn)行信號(hào)的降噪,并將振動(dòng)信號(hào)矩陣還原為一維振動(dòng)信號(hào),便于進(jìn)行后續(xù)的振動(dòng)信號(hào)解調(diào)和故障頻率提取。

2 基于K-L散度識(shí)別SVD噪聲信號(hào)

2.1 K-L散度基本理論

K-L(Kullback-Leible)散度是由庫(kù)爾貝克和萊伯勒提出并命名的,在信息論中得到廣泛應(yīng)用的信號(hào)相似性度量方法[14-16]。韓中合[14],王志杰[15],徐統(tǒng)[16]等人先后應(yīng)用K-L散度進(jìn)行軸承振動(dòng)信號(hào)的處理,并在故障信號(hào)處理上實(shí)現(xiàn)了應(yīng)用。為了更加準(zhǔn)確的區(qū)分有用信號(hào)和噪聲信號(hào),本文提出基于K-L散度識(shí)別SVD分解噪聲信號(hào)的方法。

K-L散度是量化兩個(gè)概率分布之間相似度的指標(biāo)。設(shè)兩個(gè)概率分布分別為p(x)、q(x),則q相對(duì)于p的K-L散度定義為:

(6)

D(p,q)=δ(p,q)+δ(q,p)。

(7)

K-L散度用來衡量?jī)蓚€(gè)概率分布的差異性,其物理意義是二者夾角的量度,因此,K-L散度值較大時(shí),表明兩個(gè)概率分布的差異較大;K-L散度值越小表示兩個(gè)概率分布的差異越小,即相似度越高。當(dāng)兩個(gè)概率分布完全相同時(shí),值為0。

2.2 K-L散度計(jì)算

K-L散度需要滿足的兩個(gè)條件:

1) 兩個(gè)概率分布中各自的概率和必須為1;

2) 若其中一個(gè)概率分布為0,則另一個(gè)概率分布也必須為0。

設(shè)兩信號(hào)分別為X={x1,x2,…,xn}和Y={y1,y2,…,yn},并且其概率分布為p(x)和q(y)。采用非參數(shù)估計(jì)法求解信號(hào)的概率分布,函數(shù)定義為

(8)

式(8)中,p(x)表示信號(hào)X的核密度估計(jì)函數(shù),k表示核函數(shù),h表示給定的正數(shù)(平滑參數(shù)或窗寬)。常用核密度函數(shù)是高斯核函數(shù),即

(9)

同理,可求得Y的概率分布q(y)。

將信號(hào)的概率密度函數(shù)分別代入式(6)、(7),即可求得信號(hào)X和信號(hào)Y的K-L散度值。

2.3 基于K-L散度識(shí)別SVD噪聲信號(hào)

利用K-L散度值來衡量經(jīng)過SVD分解后各個(gè)信號(hào)分量與原始信號(hào)的“關(guān)系遠(yuǎn)近”。K-L散度值越小,代表此分量與原信號(hào)的關(guān)系越近,是有用信號(hào);反之,K-L散度值越大,代表此分量與原信號(hào)關(guān)系越遠(yuǎn),是噪聲信號(hào)。具體步驟如圖1。

圖1 基于K-L散度識(shí)別SVD分解噪聲信號(hào)Figure 1 Recognition of SVD decomposition noise signal based on K-L divergence

3 基于K-L散度和TEO的滾動(dòng)軸承故障頻率提取

TEO是一種非線性差分算子,能夠增強(qiáng)信號(hào)的瞬時(shí)特性,有效提取信號(hào)中的瞬態(tài)沖擊頻率。設(shè)一組連續(xù)時(shí)間信號(hào)為x(t),TEO定義為φ,則:

(10)

滾動(dòng)軸承發(fā)生故障時(shí)會(huì)產(chǎn)生瞬態(tài)沖擊,TEO能夠及時(shí)跟蹤信號(hào)波形的變化且時(shí)間分辨率高,利用奇異值分解與K-L散度值降噪后的滾動(dòng)軸承故障信號(hào),經(jīng)過TEO解調(diào),能有效增強(qiáng)降噪信號(hào)的瞬態(tài)頻率,進(jìn)而有效提取軸承故障信號(hào)中的瞬態(tài)沖擊頻率。

基于K-L散度和TEO的滾動(dòng)軸承故障頻率提取方法步驟如圖2。

圖2 基于K-L散度和TEO的滾動(dòng)軸承故障頻率提取流程圖Figure 2 Flow chart of fault frequency extraction for rolling bearing based on K-L divergence and TEO

4 案例分析

為驗(yàn)證所提方法的有效性,本文采用美國(guó)西儲(chǔ)大學(xué)公開的軸承故障數(shù)據(jù)集進(jìn)行試驗(yàn)。該試驗(yàn)驅(qū)動(dòng)電機(jī)功率為2.2 kW、電機(jī)輸出轉(zhuǎn)速為1 730 r/min,采用電火花技術(shù)在軸承內(nèi)圈加工長(zhǎng)度為0.178 8 mm,深度為0.279 4 mm的單點(diǎn)故障,軸承信息如表1所示。

表1 軸承信息

依據(jù)軸承參數(shù),分別計(jì)算型號(hào)為SKF62052RS深溝球軸承的轉(zhuǎn)頻、內(nèi)圈故障頻率。

軸承旋轉(zhuǎn)頻率為

fo=V/60。

(11)

內(nèi)圈故障頻率計(jì)算公式為

fi=0.5fo(1+dcosα/D)z。

(12)

其中,V為軸旋轉(zhuǎn)速度,r/min。D為軸承節(jié)圓直徑,mm。d為滾動(dòng)體直徑,mm。α為接觸角,z為滾動(dòng)體數(shù)。

由式(11)、(12)計(jì)算可得滾動(dòng)軸承的旋轉(zhuǎn)頻率fo=29.16 Hz,內(nèi)圈故障頻率為f1=158.13 Hz。

本文以內(nèi)圈振動(dòng)信號(hào)為例,內(nèi)圈振動(dòng)信號(hào)的幅值圖如圖4。從圖中可以看出由于噪聲影響,無法直接識(shí)別和提取內(nèi)圈故障頻率。因此,為有效識(shí)別和提取被噪聲淹沒的故障頻率,需對(duì)內(nèi)圈振動(dòng)信號(hào)進(jìn)行降噪預(yù)處理。

根據(jù)式(3)及Hankel矩陣的構(gòu)成準(zhǔn)則,首先將滾動(dòng)軸承振動(dòng)信號(hào)構(gòu)成Hankel矩陣,采用SVD對(duì)信號(hào)矩陣進(jìn)行分解,奇異值點(diǎn)的分布如圖4所示。

圖3 原始信號(hào)幅值圖Figure 3 Amplitude diagram of original signal

圖4 奇異值分布圖Figure 4 Singular value distribution

將經(jīng)過SVD分解后的奇異值代入式(2),逐個(gè)計(jì)算信號(hào)分量與原始信號(hào)的K-L散度值,得到K-L散度值最小的信號(hào)分量,確定奇異值分解的有效階次。為了方便觀察,取前100個(gè)K-L散度值繪圖,歸一化的K-L散度值如圖5所示?？梢钥闯?第8個(gè)信號(hào)分量的K-L散度值最小,與原信號(hào)相似度最高,由式(5)計(jì)算得到降噪后的滾動(dòng)軸承振動(dòng)信號(hào),如圖6。

為驗(yàn)證本文方法的有效性,采用目前被廣發(fā)使用的有效階次確定方法——奇異值差分譜法和奇異值中值法作對(duì)比,并計(jì)算信號(hào)降噪指標(biāo)——信噪比(SNR)和均方根誤差(RMSE)。

奇異值差分譜法是將奇異值差分譜的最大突變作為有效階次,為了方便觀察,取前100個(gè)奇異差分值繪圖,如圖7所示。選取差分值最大值對(duì)應(yīng)的奇異值186.5作為閾值,選取有效階次為2,得到降噪后的滾動(dòng)軸承振動(dòng)信號(hào),如圖8。

圖5 歸一化K-L散度值Figure 5 Normalized K-L divergence

圖6 基于K-L-SVD方法降噪幅值圖Figure 6 Noise reduction amplitude image based on K-L-SVD method

圖7 奇異值差分譜Figure 7 Singular difference spectrum

圖8 奇異值差分譜法降噪幅值圖Figure 8 Singular value difference spectrum method for noise reduction amplitude diagram

奇異值中值法是將奇異值按從大到小的順序進(jìn)行排列,取奇異值的中值作為有效階次。假設(shè)奇異值個(gè)數(shù)為n,若n為偶數(shù),則中值的為第n/2個(gè)奇異值和第(n+2)/2個(gè)奇異值之和的平均數(shù),若為奇數(shù),則中值為第(n+1)/2個(gè)奇異值。由奇異值中值法降噪得到的信號(hào)如圖9。

圖9 奇異值中值法降噪幅值圖Figure 9 Singular value median method for noise reduction amplitude chart

信噪比(SNR)表示信號(hào)與噪聲的比值,信噪比的值越大則代表噪聲越小。均方根誤差(RMSE)表示計(jì)算真實(shí)值與預(yù)測(cè)值的比值,均方根誤差越小則代表降噪效果越好。計(jì)算上述三種方法所得降噪信號(hào)的SNR和RMSE如表2所示。從表中可以看出本文所提方法在增大SNR值的同時(shí)降低RMSE值,降噪效果更優(yōu)。

表2 三種方法的降噪指標(biāo)

由于滾動(dòng)軸承發(fā)生故障時(shí)會(huì)產(chǎn)生瞬態(tài)沖擊,且早期微弱的故障信號(hào)易被噪聲淹沒。利用TEO對(duì)降噪后的信號(hào)進(jìn)行解調(diào),增強(qiáng)信號(hào)中微弱的故障信號(hào),以便對(duì)滾動(dòng)軸承的故障頻率進(jìn)行準(zhǔn)確識(shí)別。經(jīng)過TEO解調(diào)后的滾動(dòng)軸承振動(dòng)信號(hào)的時(shí)頻圖如圖10。

圖10 TEO解調(diào)信號(hào)時(shí)頻圖Figure 10 Time-frequency diagram of TEO demodulation signal

從圖10中可以看出,經(jīng)過TEO解調(diào)后滾動(dòng)軸承的轉(zhuǎn)頻、內(nèi)圈故障特征頻率、內(nèi)圈故障特征頻率倍頻、內(nèi)圈故障特征頻率與轉(zhuǎn)頻之和的頻率得到增強(qiáng),說明該方法在滾動(dòng)軸承故障特征提取方面有很好的效果。

為了進(jìn)一步驗(yàn)證TEO解調(diào)對(duì)滾動(dòng)軸承故障頻率識(shí)別的有效性,選用快速傅里葉變換(Fast Fourier Transform,FFT)和包絡(luò)譜解調(diào)分別對(duì)振動(dòng)信號(hào)進(jìn)行分析,得到結(jié)果如圖11和圖12所示。

從圖中可以看出,振動(dòng)信號(hào)經(jīng)快速傅里葉變換后存在頻譜泄露的問題,使早期微弱故障信無法識(shí)別;包絡(luò)譜解調(diào)后雖然可以突出故障頻率,但存在頻譜混疊的問題,增加了故障頻率識(shí)別的復(fù)雜性。綜上所述,本文利用TEO對(duì)振動(dòng)信號(hào)進(jìn)行解調(diào)可以增強(qiáng)早期微弱故障信號(hào)并準(zhǔn)確識(shí)別故障頻率。

圖11 快速傅里葉變換信號(hào)時(shí)頻圖Figure 11 Time-frequency diagram of FFT

圖12 包絡(luò)譜解調(diào)信號(hào)時(shí)頻圖Figure 12 Time-frequency diagram of envelope spectrum demodulation signal

5 結(jié) 論

滾動(dòng)軸承振動(dòng)信號(hào)具有非平穩(wěn)、非線性、非高斯性和低信噪比的特點(diǎn),同時(shí),傳統(tǒng)的奇異值分解中閾值選取方法不能有效分離有效信號(hào)和噪聲信號(hào);然而利用K-L散度和TEO相結(jié)合的方法,能夠有效提取滾動(dòng)軸承故障頻率。

本文利用SVD與K-L散度對(duì)滾動(dòng)軸承振動(dòng)信號(hào)進(jìn)行降噪,克服了傳統(tǒng)奇異值分解確定有效階次方法中不確定性、主觀性強(qiáng)的缺點(diǎn),根據(jù)不同信號(hào)分量與原始信號(hào)K-L散度值的差異性,將噪聲信號(hào)從中分離出來,從而有效抑制噪聲的干擾,提高信號(hào)的信噪比。同時(shí),通過TEO對(duì)降噪信號(hào)解調(diào)分析,能夠增強(qiáng)信號(hào)中微弱的故障頻率,與其它信號(hào)解調(diào)方法相比,提取到的故障頻率更加全面、清晰。

本文采用美國(guó)西儲(chǔ)大學(xué)滾動(dòng)軸承數(shù)據(jù)集進(jìn)行了驗(yàn)證,為滾動(dòng)軸承的故障振動(dòng)信號(hào)分析提供了一種可行的辦法。