基于B樣條本征時間尺度分解和對角切片譜的軸承故障診斷

鐘先友，曾良才，趙春華，陳保家

滾動軸承是機械設(shè)備中最常用的、最易損傷的零部件之一，其狀態(tài)監(jiān)測和故障診斷一直受重視。滾動軸承發(fā)生故障時，其振動信號通常表現(xiàn)出非線性、非平穩(wěn)特征。在目前的非平穩(wěn)信號處理方法中，常用于滾動軸承故障診斷的方法有小波變換［1-2］、經(jīng)驗?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)［3-5］和局部均值分解(Local Mean Decomposition，LMD)［6-7］。但是這些方法都存在著一定的局限性，其中小波變換需要選擇基函數(shù)和分解層數(shù)，因而對信號的處理缺乏自適應(yīng)性;EMD存在過包絡(luò)、欠包絡(luò)、模態(tài)混淆、端點效應(yīng)等問題。LMD也存在迭代計算量大、模態(tài)混淆、端點效應(yīng)等問題。

Frei等［8］提出了一種自適應(yīng)時頻分析方法—本征時間尺度分解方法(ITD)，并將這種方法應(yīng)用于生物醫(yī)學信號處理中，取得了較好的效果。ITD方法能夠自適應(yīng)地將一個復雜信號分解為若干個瞬時頻率具有物理意義的PR分量之和。與EMD相比，ITD在計算效率方面有明顯優(yōu)勢，適合在線分解。林近山［9］將ITD算法成功地應(yīng)用于齒輪箱的故障診斷中。但是ITD方法中對基線的定義是基于信號本身的線性變換，因此從第二個分量開始，有明顯的信號失真。鄭進德和楊宇等［10-11］采用三次樣條插值對ITD方法進行改進并應(yīng)用到軸承和轉(zhuǎn)子的故障診斷中。但三次樣條擬合曲線時會產(chǎn)生的過包絡(luò)、欠包絡(luò)問題［12］，文獻［12-13］采用B樣條插值分別改進EMD和LMD算法，并驗證了B樣條插值的優(yōu)越性。因此，本文提出了B樣條插值改進的本征時間尺度分解方法(BITD)。

Teager能量算子適合檢測軸承振動信號中的沖擊成分，對信號的瞬時變化具有良好自適應(yīng)能力，而且計算復雜性低，計算效率高。文獻［4，14］運用Teager能量算子成功地識別出滾動軸承的故障特征頻率。

包絡(luò)分析技術(shù)是軸承故障振動信號分析的最有效方法之一。在進行包絡(luò)分析之前，為了消除噪聲的干擾，須先對振動信號進行高通濾波，以消除低頻噪聲的干擾，但高頻濾波的中心頻率和帶寬往往難以選擇，濾波參數(shù)的選擇對分析結(jié)果影響很大［15］。滾動軸承發(fā)生故障時產(chǎn)生的周期性沖擊引起軸承系統(tǒng)的高頻固有振動，故障特征被調(diào)制到高頻段，采用BITD方法將軸承振動信號分解為若干個PR分量，這些分量由高頻到低頻依次被分解出來，選取分解結(jié)果的前幾個分量重構(gòu)就可能突出故障信號的高頻段，從而克服了包絡(luò)分析需要預(yù)先確定濾波器中心頻率和帶寬的難題。

實際軸承振動信號中往往混入了大量的噪聲，使得原始故障特征信息與噪聲混淆而不易提取，對角切片譜具有抑制高斯白噪聲的特性，同時可以識別軸承振動信號中的二次相位耦合成分［16］。

基于以上原因，本文提出了基于B樣條的本征時間尺度分解、Teager能量算子解調(diào)和對角切片譜相結(jié)合的軸承故障診斷方法，成功地提取出滾動軸承的故障特征。

1 ITD方法

ITD方法能夠自適應(yīng)地將一個復雜信號分解為若干個相互獨立的合理旋轉(zhuǎn)分量和一個趨勢項之和。設(shè)Xt是待分析的原信號，分解前先定義一個基線提取算子L，使得從原始信號中去掉該基線后剩下的余量信號成為一個合理旋轉(zhuǎn)分量。一次分解的表達式為［8］:

(1)確定信號Xt的極值點Xk及對應(yīng)的時刻τk(k=1，2，…M，M為所有極值點的個數(shù))，并計算

式中:0 ＜a＜1，一般地，a=0.5。

(2)定義信號的分段線性基線提取算子如下:

(3)將基線信號Lt作為原始信號，重復上述步驟，直到基線信號為一單調(diào)函數(shù)或常函數(shù)。原始信號被分解為:

考察式(5)所示信號

式中:x1(t)為調(diào)頻信號，x2(t)和x3(t)為兩個余弦信號，時域波形如圖1。對x(t)別進行ITD分解，得到的分解結(jié)果如圖2所示，其中PR1、PR2和PR3分別為ITD分解的前三個分量，R為殘余分量。

圖1 仿真信號的時域波形Fig.1 The time domain of simulated signal

圖2 仿真信號的ITD分解結(jié)果Fig.2 ITD decomposition results of the simulation signal

從圖2中可以看出，PR1分量、PR2分量和PR3分量分別對應(yīng)于仿真信號x(t)的三個分量x1(t)、x2(t)和x3(t)，第二個分量和第三個分量出現(xiàn)了失真，這是由于ITD方法中是以原始信號任意兩個相鄰的極值點為跨度對信號進行分段線性變換來構(gòu)造基線信號，這導致第二個分量和第三個分量信號波形出現(xiàn)了毛刺而失真。因此，本文對ITD方法進行改進，用B樣條插值來代替ITD方法中的線性變換。

2 BITD方法

BITD方法基本的分解過程如下:

(1)確定原始信號Xt所有的局部極值點，方法與ITD相同，通過式(2)和(3)和計算各基線的控制點Xk。

(2)對序列端點采用鏡像對稱延拓方法進行延拓，得到左右兩端點極值(τ0，X0)和(τM+1，XM+1)，令 k 分別等于0和M，按式(2)求出L1與LM的值，然后對所有Lk用B樣條函數(shù)來擬合，得到基線信號L1。

(3)將L1從原始信號中分離出來，得到P1，若P1是一個PR分量，則P1作為信號Xt的第一個分量，否則將P1作為原始信號重復上述步驟，循環(huán)k次，直到得到Pk是一個PR分量，Pk即為信號Xt的第一個PR分量PR1，將PR1信號中分離出來，得到一個新的信號r1。

圖3 仿真信號的BITD分解結(jié)果Fig.3 BITD decomposition results of the simulation signal

(4)再將r1作為原始信號重復上面的步驟，得到Xt的第二個滿足條件的PR2。重復循環(huán)n次，得到信號Xt的n個滿足PR條件的分量，直到rn為一單調(diào)函數(shù)或常函數(shù)，這樣就把Xt分解為n個PR分量和一個單調(diào)或常函數(shù)之和，即:

用BITD方法對圖1仿真信號x(t)進行分解，得到的分解結(jié)果如圖3所示，從圖中可以看出，BITD方法對第二個分量和第三個分量的分解也取得較好的效果。

文獻［12］采用B樣條插值改進EMD算法，簡稱BEMD，現(xiàn)對于圖1所示的信號，分別用EMD、BEMD、ITD和BITD進行分解，比較四者的分解效果，在同一臺電腦上，各運行20次，取平均值，得到三種方法分解的時間用t來表示，用r1和r2表示分解得到的第一個分量和第二個分量與真實信號的相關(guān)系數(shù)。結(jié)果與表1所示。

表1 四種方法分解效果比較Tab.1 Comparison of four methods of decomposition

從表1可以看出，ITD分解速度最快，但分解所得的分量與真實分量的相關(guān)性最小，而BITD分解速度比EMD和BEMD快，且分解所得的分量與真實分量的相關(guān)性最大，相比另外三種分解方法具有一定的優(yōu)勢。

3 Teager能量算子的基本原理

Teager能量算子計算簡單，計算效率高。在機械故障診斷中，能量算子解調(diào)方法被用于處理振動信號來提取故障信息［4，14］。連續(xù)信號 x(t)的 Teager能量算子可定義為:

式中:x(t)為測得振動信號，x'(t)和x″(t)分別為信號x(t)的一階和二階導數(shù)。

離散信號x(n)的Teager能量算子可定義為:

對于離散時間信號，能量算子只需要三個樣本數(shù)據(jù)就可以計算任意時刻n處的信號源能量。

4 對角切片譜的基本原理

軸承發(fā)生故障時，系統(tǒng)會表現(xiàn)出一定的非線性，最常見的表現(xiàn)形式就是二次相位耦合。對于這種非線性耦合現(xiàn)象，僅用基于二階統(tǒng)計量的分析方法如自相關(guān)、功率譜來處理，很難將故障特征提取出來，這主要是由于二階統(tǒng)計量不提供任何相位信息。對角切片譜能檢測二次相位耦合特征，并且當信號中混有高斯噪聲時，理論上可被對角切片譜完全抑制掉。

對于隨機變量 x(t)，它的三階累積量 c3x(τ1，τ2)(τ1、τ2為時間延遲)的對角切片為 c3x(τ，τ)，定義該對角切片的傅里葉變換為隨機變量x(t)的對角切片譜:

5 基于BITD、能量算子和對角切片譜的診斷方法

將基于BITD、Teager能量算子和對角切片譜的診斷方法應(yīng)用到滾動軸承診斷中，主要包括以下步驟:

(1)對信號x(t)進行BITD分解，得到若干個PR1，PR2…PRn分量;

(2)計算各PR分量與原信號的相關(guān)系數(shù);

(3)計算各PR分量的樣本熵;

(4)取相關(guān)系數(shù)和樣本熵值都較大的PR分量重構(gòu)信號，對重構(gòu)信號進行能量算子解調(diào)并求對角切片譜。重構(gòu)的分量個數(shù)要通過分析來取值，本文取相關(guān)系數(shù)大于0.3和樣本熵值大于1的PR分量重構(gòu)信號。

6 軸承故障仿真信號分析

根據(jù)文獻［17］建立滾動軸承元件發(fā)生單點局部損傷時傳感器所采集到的信號模型為:

其中:m(t)是沖擊幅值，是幅值調(diào)制函數(shù);f1是第一調(diào)制頻率，等于軸的轉(zhuǎn)頻或滾動體的公轉(zhuǎn)頻率;T為故障特征周期;f2是軸承座-傳感器系統(tǒng)的某一固有頻率，即載波頻率;c為沖擊信號衰減指數(shù);U(t)為單位階躍函數(shù);n(t)為噪聲。其中取固有頻率f2為2.5 kHz，阻尼系數(shù)為c=0.1，n(t)=0，T=1/200，m(t)=1。根據(jù)式(10)和式(11)可得到外圈故障的時域波形和頻譜圖如圖4。向信號中添加Gauss白噪噪聲，使得信噪比為-10 dB，得到外圈故障的時域波形和頻譜圖如圖5和圖6。

圖5 加噪后軸承故障模擬信號的時域圖Fig.5 The time domain of simulate bearing failure signal with noise

采用BITD方法對該故障仿真信號進行分解，得到8個PR分量和一個殘余分量，前4個分量如圖7所示，計算前3個PR分量與原信號的相關(guān)系數(shù)和PR分量的樣本熵，如表2所示。

圖6 加噪后軸承故障模擬信號的頻譜圖Fig.6 The frequency domain of simulate bearing failure signal with noise

圖7 BITD方法對圖5中軸承信號的分解結(jié)果Fig.7 The decomposed results of bearing signal shown in Fig.5 by BITD

表2 PR分量與原信號的相關(guān)系數(shù)及樣本熵Tab.2 The correlation coefficients between PR andoriginal signal，and the PR sample entropy

從表2看出，前兩個分量相對較大，對前兩分量進行重構(gòu)，得到分量PR12，計算PR12分量與原信號的相關(guān)系數(shù)和樣本熵，其值分別為0.82和1.32，PR12相比單個PR分量提高了相關(guān)系數(shù)和樣本熵值，保留了原始信號更多的沖擊特征信息，故對PR12分量進行能量算子解調(diào)后求對角切片譜，結(jié)果如圖8。從圖8中可以清楚地看出，在故障特征一倍頻(200 Hz)及其倍頻處有明顯的譜線，表明本文提出的方法可以將故障頻率成分成功地提取出來。

為了進行對比分析，對原信號進行Hilbert變換求包絡(luò)譜，如圖9所示，在包絡(luò)譜中可以找到故障特征頻率一倍頻(200 Hz)，但故障特征頻率二倍頻和三倍頻無法識別，說明本文所提方法效果更好。

對PR12分量進行Hilbert變換求包絡(luò)譜，如圖10所示，在包絡(luò)譜中可以找到故障特征頻率一倍頻(200 Hz)和故障特征的三倍頻，對比圖9和圖10可知，可以看出，經(jīng)BITD分解后作包絡(luò)譜分析，可以進一步削弱部分干擾成分，提高信噪比，而與圖8比較，可以看出本文提出的方法分析結(jié)果更為準確，信噪比更高。

圖8 PR12分量的對角切片譜Fig.8 Diagonal slice spectrum of PR12 component

圖9 圖5中軸承信號的包絡(luò)譜Fig.9 The envelope spectra of bearing signal shown in Fig.5

圖10 PR12分量的包絡(luò)譜Fig.10 The envelope spectra of PR12 component

7 發(fā)電機滾動軸承故障診斷實例

圖11是某公司850 kW發(fā)電機軸承振動信號的時域波形及頻譜。軸承故障為外圈損傷，軸承型號6326，采樣頻率25 600 Hz。轉(zhuǎn)頻為24.33 Hz，軸承外圈故障特征頻率為85.8 Hz。從圖11可以看出，時域波形比較復雜，難以分辨出信號的具體特征。在頻譜圖中，故障信號的低頻特征淹沒在背景噪聲中，無法識別故障頻率及其諧波。

圖11 軸承振動信號時域圖和頻譜Fig.11 The time domain and frequency domain of bearing vibration signal

采用BITD方法對該故障仿真信號進行分解，得到10個PR分量和一個殘余分量，前3個分量如圖12所示，計算PR與原信號的相關(guān)系數(shù)和PR的樣本熵，如表3所示。

圖12 BITD方法對軸承信號的分解結(jié)果Fig.12 The decomposed results of bearing signal by BITD

表3 PR分量與原信號的相關(guān)系數(shù)及樣本熵Tab.3 The correlation coefficients between PR and original signal，and the PR sample entropy

圖13 PR12分量的對角切片譜Fig.13 Diagonal slice spectrum of PR12 component

圖14 PR12分量的包絡(luò)譜Fig.14 The envelope spectra of PR12 component

圖15 重構(gòu)信號的對角切片譜Fig.15 diagonal slice spectrum of the reconstruction signal

從表3看出，前兩個分量相對較大，對前兩分量進行重構(gòu)，得到分量PR12，計算PR12分量與原信號的相關(guān)系數(shù)和樣本熵，其值分別為0.961和2.103，PR12相比單個PR分量提高了相關(guān)系數(shù)和樣本熵值，故對PR12分量進行能量算子解調(diào)后求對角切片譜，結(jié)果如圖13。從圖中可以清楚地看出，在故障特征一倍頻、二倍頻和三倍頻處有明顯的譜線，表明本文所提出的方法可以將故障頻率成分成功提取出來。

對PR12分量進行Hilbert變換求包絡(luò)譜，如圖14所示，在包絡(luò)譜中可以找到與故障特征頻率一倍頻(85.8 Hz)相近的頻率和故障特征的二倍頻，但一倍頻處存在干擾成分，且故障特征的三倍頻無法識別，對比圖13和圖14可知，可以看出本文所提出的方法的分析結(jié)果更為準確，信噪比更高。

對振動信號進行EMD分解，求各分量與原信號的相關(guān)系數(shù)及樣本熵，前兩個分量IMF1和IMF2與原信號的相關(guān)系數(shù)大于0.3，樣本熵值大于1，故對前兩分量進行重構(gòu)，對重構(gòu)信號進行能量算子解調(diào)后求對角切片譜，結(jié)果如圖15。從圖15中可以看出，在故障特征一倍頻(85.8 Hz)、二倍頻和三倍頻附近處有較明顯的譜線，但頻率值沒有圖13中的準確，且圖13中二倍頻更為明顯，表明本文所提出的方法識別效果更好。

8 結(jié)論

(1)本文所提出的BITD方法是一種自適應(yīng)時頻分析方法，可以解決ITD方法分解信號產(chǎn)生的波形失真問題;

(2)比較對原信號直接做包絡(luò)譜分析和對信號進行BITD分解后作包絡(luò)譜分析的效果，結(jié)果表明BITD分解可以削弱部分干擾成分，提高信噪比，但與本文所提出的方法相比，基于BITD、Teager能量算子和對角切片譜的分析結(jié)果更為準確，信噪比更高;

(3)BITD方法計算效率高，能量算子解調(diào)和對角切片譜，快速，易于實現(xiàn)，仿真信號與軸承故障診斷工程實例的分析表明，本文所提出的方法能有效地提高信噪比，突出故障特征，提高故障診斷的準確性，具有良好的應(yīng)用前景。

［1］王曉冬，何正嘉，訾艷陽.滾動軸承故障診斷的多小波譜峭度方法［J］.西安交通大學學報，2010，44(3):77-8 WANG Xiao-dong，HE Zheng-jia，ZI Yan-yang.Spectral kurtosis of multiwavelet for fault diagnosis of rolling bearing［J］.Journal of Xi'an Jiaotong University，2010，44(3):77-81.

［2］彭志科，何永勇，褚福磊.小波尺度譜在振動信號分析中的應(yīng)用研究［J］.機械工程學報，2002，38(3):122-126.PENG Zhi-ke，HE Yong-yong，CHU Fu-lei.Using wavelet scalogram for vibration signals anaylsis［J］.Chinese Journal of Mechanical Engineering，2002，38(3):122-126.

［3］湯寶平，蔣永華，張詳春.基于形態(tài)奇異值分解和經(jīng)驗?zāi)B(tài)分解的滾動軸承故障特征提取方法［J］.機械工程學報，2010，46(5):37-42.TANG Bao-ping，JIANG Yong-hua， ZHANG Xiang-chun.Feature extraction method of rolling bearing fault based on singular value decomposition-morphology filter and empirical mode decomposition［J］.Chinese Journal of Mechanical Engineering，2010，46(5):37-42.

［4］李 輝，鄭海起，楊紹普.基于EMD和Teager能量算子的軸承故障診斷研究［J］.振動與沖擊，2008，27(10):15-22.LI Hui， ZHENG Hai-qi， YANG Shao-pu. Bearing fault diagnosis based on EMD and teager kaiser energy operator［J］.Journal of Vibration and Shock，2008，27(10):15-22.

［5］蘇文勝，王奉濤，張志新，等.EMD降噪和譜峭度方法在滾動軸承早期故障診斷中的應(yīng)用［J］.振動與沖擊，2010，29(3):18-21.SU Wen-sheng，WANG Feng-tao，ZHANG Zhi-xin，et al.Application of EMD denoising and spectral kurtosis in early fault diagnosis of rolling element bearings［J］.Journal of Vibration and Shock，2010，29(3):18-21.

［6］程軍圣，楊 怡，楊 宇.基于LMD的能量算子解調(diào)機械故障診斷方法［J］.振動、測試與診斷，2012，32(6):915-918.CHENG Jun-sheng，YANG Yi，Yang Yu.Energy operator demodulation mechanical fault diagnosis method based the LMD ［J］.Journal of Vibration，Measurement ＆ Diagnosis，2012，32(6):915-918.

［7］楊 宇，王歡歡，程軍圣，等.基于LMD的包絡(luò)譜特征值在滾動軸承故障診斷中的應(yīng)用［J］.航空動力學報，2012，27(5):1153-1158.YANG Yu，WANG Huan-huan，CHENG Jun-sheng，et al.Application of envelope spectrum characteristics based on LMD to roller bearing fault diagnosis［J］.Journal of Aerospace Power，2012，27(5):1153-1158.

［8］Frei M G，Osorio I.Intrinsic time-scale decomposition:analysis and real-time filtering of non-stationary signals［J］.Proceedings of the Royal Society，2007，463:321-342.

［9］林近山.基于本征時間尺度分解算法的齒輪箱故障診斷［J］.機械傳動，2011，35(9):51-53.LIN Jin-shan.Fault diagnosis of gear box based on Intrinsic time-scale decomposition algorithm［J］.Journal of Mechanical Transmission，2011，35(9):51-53.

［10］鄭近德，程軍圣，楊 宇.基于改進的ITD和模糊熵的滾動軸承故障診斷方法［J］.中國機械工程，2012，23(19):2372-2377.ZHENG Jin-de，CHENG Jun-sheng，YANG Yu.A rolling bearing fault diagnosis method based on improved ITD and fuzzy entropy ［J］.Chinese Mechanical Engineering，2012，23(19):2372-2377.

［11］楊 宇，王歡歡，程軍圣.基于ITD改進算法和關(guān)聯(lián)維數(shù)的轉(zhuǎn)子故障診斷方法［J］.振動與沖擊，2012，31(23):67-70.YANG Yu，WANG Huan-huan，CHENG Jun-sheng.A rotor fault diagnosis method based on ITD improved algorithm and correlation dimension ［J］.Journal of Vibration and Shock，2012，31(23):67-70.

［12］邱綿浩，劉 箐，叢 華.基于B樣條插值曲線直接篩選的EMD及其在機械振動信號處理中的應(yīng)用［J］.裝甲兵工程學院學報，2007，21(3):29-33.QIU Mian-hao，LIU jing，CONG Hua.Research of direct sifting EMD based on cubic B-Spline interpolation curve and its application in processing mechanical vibrating signals［J］.Journal of Academy of Armored Force Engineering，2007，21(3):29-33.

［13］王明達，張來斌，梁 偉，等.基于B樣條插值的局部均值分解方法研究［J］.振動與沖擊，2010，29(11):73-77.WANG Ming-da，ZHANG Lai-bin，LIANG Wei，et al.Local mean decomposition method based on B-spline interpolation［J］.Journal of Vibration and Shock，2010，29(11):73-77.

［14］王天金，馮志鵬，郝如江，等.基于Teager能量算子的滾動軸承故障診斷研究［J］.振動與沖擊，2012，31(2):1-5.WANG Tian-jin，F(xiàn)ENGZhi-peng，HAORu-jiang，et al.Fault diagnosis of rolling element bearings based on Teager energy operator［J］.Journal of Vibration and Shock，2012，31(2):1-5.

［15］高 強，杜小山，范 虹，等.滾動軸承故障的EMD診斷方法研究［J］.振動工程學報，2007，20(1):15-18.GAO Qiang，DU Xiao-shan，F(xiàn)AN Hong，et al.An empirical mode decomposition based method for rolling bearing fault diagnosis［J］.Journal of Vibration Engineering，2007，20(1):15-18.

［16］ Kachenoura A，Albera L，Bellanger J J，et al.Nonminimum phase identification based on higher order spectrum slices［J］.IEEE Transaction on Signal Processing，2008，56(5):1821-1829.

［17］Mc Fadden PD，Smith JD.Model for the vibration produced by a single point defect in a rolling element bearing［J］.Journal of Sound and Vibration，1984，96(1):69-82.