基于迭代廣義解調算法的變轉速滾動軸承復合故障特征的提取

趙德尊， 李建勇, 2， 程衛東， 溫偉剛

(1.北京交通大學 機械與電子控制工程學院，北京 100044; 2.載運工具先進制造與測控技術教育部重點實驗室, 北京 100044)

滾動軸承是旋轉機械中的重要部件，由于復雜的工作環境(高溫、重載)及運行狀態(急劇變轉速、變負載)，滾動軸承也是旋轉機械中最易受損的零部件之一。據有關資料統計：在使用滾動軸承的旋轉機械中有超過30%的機械故障是故障軸承引起的[1]，而且實際工況中，故障點往往不止一處，某些故障常常會誘發其位置的故障發生，進而出現多個故障并存的狀態[2]。另外，變轉速模式普遍存在于旋轉機械的實際運行中[3-4]。該條件下，由故障引起的沖擊間隔將隨轉速變化發生改變，使得頻率譜中出現頻率模糊現象，無法有效的完成故障特征的識別。因此，對變轉速條件下滾動軸承復合故障特征提取的研究具有重要的意義。

近年來，許多學者致力于滾動軸承復合故障提取的研究，并取得了很多成果。明安波等[5]將小波分解與頻譜自相關相結合，通過對單通滾動軸承道復合故障信號的分解，有效的實現了相互交錯、疊加在同一頻帶的內、外圈故障特征的分離和提取。馬新娜等[6]將EMD算法與自適應陷波器系統相結合在一定程度上滿足了滾動軸承復合故障信號分離和故障診斷的要求。Chen等[7]提出了基于改進的自適應沉余提升多小波和希爾伯特解調方法相結合的滾動軸承復合故障特征提取方法。仿真和實驗對比結果表明了該算法相對于傳統算法的優越性；Zhang等[8]則提出了基于偽時頻分析與Dopplerlet濾波器相結合的滾動軸承道旁聲音信號復合故障特征提取方法。該方法可以有效的利用道旁聲音信號實現滾動軸承的復合故障診斷。然而上述方法都受制于轉速恒定或小范圍波動這一前提條件，對于轉速大范圍變化的故障軸承信號無能為力。

針對變轉速工作模式對故障軸承信號產生的影響，計算階比分析是公認的最為有效的算法[9-10]。該算法的核心是獲得相對于參考軸的恒定角增量采樣，將時域的非周期信號轉化為角域周期信號以消除轉速變化對振動信號的影響。然而，如果利用計算階比分析算法對復合故障軸承信號進行處理，由于復合故障相互干擾、彼此疊加，將使得包絡階比譜中代表故障特征的峰并不突出以及具有較多的干擾成分，因此容易出現漏判或誤判。另外，階比分析方法存在理論上的缺陷，即重采樣過程中大量二次或高次方程的求解使得計算效率較低[11]，以及包絡畸變等原因容易產生誤差[12]，進而影響故障診斷的結果。因此有必要提出新的方法以消除變轉速工作模式對故障軸承振動信號的影響。

廣義解調算法是Olhede等[13]提出的一種可以將時頻分布是傾斜非平穩的信號轉換成時頻分布是線性的且平行于時間軸的平穩信號的一種分析方法，并將其成功的應用于語音信號。近年來，基于廣義解調算法適用于調幅-調頻信號的特點，許多學者將其應用于機械設備的故障診斷[14-16]。然而，廣義解調算法只能對信號時頻譜中的單一曲線進行識別與重置，難以處理時頻圖中含有非平行多曲線的多分量信號。因此有學者對廣義解調算法進行了改進和優化，提出了迭代廣義解調算法，并將其用于行星齒輪箱的故障診斷[17-18]。

綜上所述，本文提出了基于迭代廣義解調算法的變轉速滾動軸承復合故障特征提取方法。該方法的核心是利用迭代廣義解調算法將復合故障軸承信號的包絡時頻譜中與故障特征有關的多條特定時頻成分轉換成平行于時間軸的直線，使其滿足快速傅里葉變換的要求，進而通過識別頻率譜中突出峰的位置確定目標軸承的故障數量及故障位置。基于迭代廣義解調算法的變轉速滾動軸承復合故障特征提取方法的優點主要包括以下兩個方面：首先，舍棄計算階比分析方法，直接通過迭代廣義解調算法消除轉速變化對滾動軸承故障特征的影響，為變轉速條件下滾動軸承復合故障特征的提取提供了新思路；其次，該方法通過相位函數只對故障軸承信號中的特定成分進行分析，使得頻譜中不存在明顯的峰值干擾，易于故障特征的提取和故障點位置的判別。

1 算法部分

1.1 復合故障軸承信號的時頻特性

當軸承某一位置出現表面損傷類故障時，故障點與其對應的配合表面之間因為碰撞會在相應的振動信號中產生一個高幅值并且快速衰減的沖擊。隨著滾動軸承的不斷運行，上述沖擊將會以固定的時間間隔重復出現，其對應的重復頻率即為故障特征頻率。故障特征頻率公式可參考文獻[10]。根據故障特征頻率計算公式得知滾動軸承內、外圈及滾動體對應的故障特征頻率與軸承轉頻具有固定的比例關系且只與目標軸承幾何參數有關。將該比例關系定義為故障特征系數，其公式形式如下：

Fo=fo/fr

(1)

Fi=fi/fr

(2)

Fb=fb/fr

(3)

式中：Fo、Fi和Fb分別表示目標軸承內圈、外圈及滾珠故障時對應的故障特征系數;fo、fi和fb分別表示目標軸承的外圈、內圈及滾珠故障時對應的故障特征頻率；fr為目標軸承轉頻。

圖1 復合故障軸承仿真信號的包絡時頻圖Fig.1 Envelope time-frequency representation of the multi-fault bearing signal

故障點與其對應的配合表面之間因碰撞而產生的故障沖擊將進一步引起機械結構的高頻共振。因此，故障軸承振動信號可以看作以機械結構的共振頻率為載頻以故障特征系數為調頻的幅值解調信號。利用Hilbert變換對故障軸承振動信號進行解調，解調信號的時頻譜中將會出現成特定規律變化的時頻曲線。上述時頻曲線即為故障特征頻率趨勢線及其倍頻。當軸承出現多個故障時，其包絡信號的時頻圖中同樣出現多組故障特征頻率趨勢線及其倍頻，如圖1所示(信號為仿真信號，具體參數見仿真部分)。圖中1和2分別表示故障點1和故障點2所對應的故障特征頻率趨勢線，1-2、1-3以及2-2，2-3分別表示兩處故障點故障特征頻率的諧波。由于轉速曲線成線性變化，根據公式(1)～(3)，故障特征頻率也以同樣的變化趨勢成線性變化。這也是造成以快速傅里葉變換為核心的頻譜分析方法失效，產生頻率模糊現象的原因。

根據以上分析，如果能夠利用轉頻信息以及目標軸承的故障特征系數對其包絡時頻圖中的故障特征頻率進行識別和重置，使其成為平行與時間軸的直線，即可利用頻譜分析對故障特征頻率進行提取和識別。進而確定目標軸承的故障數量以及類型。

1.2 迭代廣義解調算法

廣義解調可以通過預設的相位函數把信號時頻分布中的特定曲線成分轉換成線性的、平行于時間軸的直線。其本質為廣義傅里葉變換。對于任意單分量信號x(t)，其廣義傅里葉變換的定義為：

(4)

式中：s0(t)表示隨時間t變化的實值函數，實際上是對x(t)e-2jπs0 (t )做標準的傅里葉變換。

實際工況中，多分量信號普遍存在于旋轉機械中。而原始廣義解調算法只能對原始信號中的單一分量進行單次解調。為有效的解調多分量信號中的多個特定成分，有學者對廣義解調進行了優化，提出了迭代廣義解調算法。迭代廣義解調算法的核心通過改變預設的相位函數，對原始信號進行反復的廣義解調。其具體步驟如下：

(1)通過傳統的時頻分析方法如短時傅里葉變換確定原始信號x(t)中的特定時頻成分；

(2)利用數據擬合算法估計原始信號中每一成分xi(t)的瞬時頻率fi(t)；

(4)對原始信號x(t)進行Hilbert變換，獲得解析信號y(t)=x(t)+jH[x(t)]，其中H[x(t)]是x(t)的Hilbert變換；

(6)對解調信號d(t)進行Hilbert變換得到新的解析信號z(t)=d(t)+jH[d(t)]，其中H[d(t)]是d(t)的Hilbert變換。

1.3 算法流程圖

變轉速條件下，滾動軸承復合故障信號的包絡時頻圖中存在成規律分布且能夠代表軸承故障特征的時頻曲線。上述時頻曲線與轉頻之間有固定的比例關系，比例系數即為故障特征系數。因此本文利用故障軸承信號的同步轉頻信息以及目標軸承的故障特征系數確定迭代廣義解調算法所用到的相位函數，結合迭代廣義解調算法可以將多分量信號中特定時頻曲線轉換成平行于時間軸的直線這一優點，提出了基于迭代廣義解調算法的變轉速滾動軸承復合故障特征提取方法。該方法的流程圖如圖2所示，具體步驟如下：

(1)利用同步測取的目標軸承轉頻信號獲取轉頻曲線；

(2)根據轉頻曲線方程以及目標軸承的故障特征系數計算迭代廣義解調算法所需的相位函數及相位點；

(3)通過Hilbert變換獲取目標軸承的包絡信號；

(4)對目標軸承的包絡信號進行迭代廣義解調得到解調信號；

(5)通過快速傅里葉變換計算解調信號的頻率譜；

(6)判斷頻率譜中突出峰的橫坐標值與相位點的關系完成故障診斷。

圖2 算法流程圖Fig.2 Flowchart of the proposed method

2 仿真驗證

本節構造了變轉速工作模式下滾動軸承復合故障仿真信號對本文算法的有效性進行驗證。變轉速條件下滾動軸承單一故障仿真信號的構造公式如下：

(5)

式中：N為信號的長度；Am=a·tm+b代表第m個沖擊的幅值；β為結構的衰減系數；wr表示軸承故障激起的共振頻率；μ(t)為單位階躍函數；tm表示第m個沖擊出現的時間，計算公式如下：

(6)

式中：m=2,3,…,N，f(t)=6t+38表示軸承轉頻隨時間變化的規律；τ代表由滾動體滑移帶來的故障沖擊間隔之間的誤差，其取值一般為0.01～0.02；n代表軸承每轉出現的故障沖擊數。

根據公式(5)和(6)構造變轉速條件下滾動軸承復合故障仿真信號：

xbearing(t)=x1(t)+x2(t)+n(t)

(7)

式中：x1(t)是故障點1引起的沖擊脈沖序列，x2(t)是故障點2引起的沖擊脈沖序列，n(t)為高斯白噪聲。仿真模型的其他參數見表1。

表1 變轉速條件下滾動軸承復合故障仿真模型參數

根據上述仿真模型以及參數確定的變轉速條件下

滾動軸承復合故障引起的沖擊響應(局部)如圖3所示。根據圖3可以看出，兩個故障點產生的沖擊相互交叉在一起，有的沖擊甚至重疊在一起。仿真信號整體時域波形圖如圖4所示。

根據仿真信號的頻率曲線方程f(t)=6t+38以及故障特征系數確定迭代廣義解調算法所需要的相位函數以及相位點如表2所示。為保證及驗證故障診斷結果的可靠性，結合故障軸承信號包絡時頻圖中的故障特征頻率曲線擁有多條諧波，且各個諧波相對于故障特征頻譜分別具有N(N=2,3,…)倍的倍數關系這一特點，本文增加了對故障特征頻率趨勢線的2倍及3倍諧波的處理。

圖3 滾動軸承復合故障沖擊脈沖仿真示意圖Fig.3 Multi-fault impulses of simulated signal

圖4 變轉速條件滾動軸承復合故障仿真信號Fig.4 Simulated multi-fault vibration signal

為對比分析，隨機選取了一個故障特征系數F3=2.5，以假設該故障軸承還包含第3個故障點。該故障點對應的故障特征頻率方程、相位函數以及相位點如表2所示。

表2 故障軸承仿真信號各個故障對應的相位函數及相位點

根據表2中故障點1和2分別對應的3個相位函數以及假設的故障點3對應的相位函數逐次對仿真信號的包絡信號進行廣義解調，得到解調信號后對解調信號進行FFT，計算的頻率譜如圖5所示。最后利用表中的各個故障對應的相位點以及頻率譜中的突出峰橫坐標完成故障診斷。

圖5中，以圓圈標注的突出峰橫坐標分別為117.5、235和352.5。三組數據分別與故障點1對應的相位點數值相吻合。而圖中以箭頭標注的3個突出峰對應的橫坐標分別為163、326.5和489，分別近似于故障點2所對應的3個相位點。因此可以斷定仿真信號中存在兩種故障。

另外，圖5中在假設的故障點3所對應的相位點95 Hz處不存在明顯且獨立的突出峰，以此可以說明該仿真信號中不存在所假設的故障點3。

綜上分析，本文提出的基于迭代廣義解調算法的滾動軸承復合故障診斷方法可以有效的確定目標軸承的故障點的數量和位置。

圖5 解調信號的頻譜圖Fig.5 Frequency spectrum of demodulated signal

3 實驗驗證

本節利用滾動軸承振動試驗臺上的實測數據對本文算法進行驗證。圖6為實驗裝置。其中加速度傳感器安裝在目標軸承附近以準確測取其振動信號。轉速計安裝在軸端用于測量軸承轉速。采集裝置為 YE6231 采集卡及其配套的軟件。對滾動軸承進行電火花切割模擬其外圈和內圈裂紋復合故障，其中故障程度分別為輕度故障(寬：0.2 mm,深：0.4 mm)、中度故障(寬：0.4 mm,深：0.6 mm)和重度故障(寬：0.6 mm,深：0.8 mm)。本文選取內圈輕度、外圈中度故障軸承作為實驗軸承如圖7所示。目標軸承的型號、幾何參數以及內外圈和滾動體對應的故障特征系數分別見表3。

圖8為測取包含內外圈故障的目標軸承在減速條件下的振動信號，其中采樣頻率為24 000 Hz，時長為3 s。圖9中左下角為同步測取轉速脈沖信號(局部)，從中也可以看出隨著轉速的減小，脈沖逐漸稀疏。根據轉速脈沖信號計算得到的轉速曲線如圖9所示，其擬合方程為f(t)=-0.7t2-5.8t+52.9。

圖6 滾動軸承試驗臺Fig.6 Experimental setup of rolling element bearing

圖7 復合故障軸承Fig.7 Target bearing with outer race defect and inner race defect

參數數值軸承型號6000滾動體數n7滾動體直徑d/mm4.8節圓直徑D/mm17.65接觸角α0外圈故障特征系數2.5內圈故障特征系數4.4滾珠故障特征系數1.7

圖8 減速條件下復合故障軸承實測信號Fig.8 Measured signal of multi-fault bearing under deceleration

根據轉頻曲線方程、目標軸承的故障特征系數以及公式(1)～(3)計算得到內、外圈和滾動體對應的故障特征頻率趨勢線方程。對上述方程分別乘以N(N=2, 3)即可得到內、外圈和滾動體故障特征頻率趨勢線的二倍及三倍頻趨勢線方程。上述故障特征頻率趨勢線方程及其2、3倍頻方程以及分別對應的相位函數、相位點分別見表4。

分別根據表4中的9個相位函數對實測信號的包絡信號逐次廣義解調。解調信號的頻率譜如圖10所示。圖中含有成規律分別且相對獨立的突出峰。以圓圈標注的突出峰的橫坐標分別為132.2、264.8和396.2。

這與目標軸承外圈故障對應的3個相位點相吻合，即可說明該軸承含有外圈故障。

圖9 局部轉速脈沖信號及轉頻曲線Fig.9 Rotational frequency curve and partial rotational speed signal

故障位置故障特征系數故障特征頻率方程相位函數相位點外圈Fo=2.5fo(t)=-1.75t2-14.5t+132.3so(t)=-0.56t3-7t2fpo=132.32fo(t)=-3.5t2-29t+264.62so(t)=-1.12t3-14t2fpo=264.63fo(t)=-5.52t2-43.5t+396.93so(t)=-1.68t3-21t2fpo=396.9內圈Fi=4.4fi(t)=-3.1t2-25.5t+232.8si(t)=-1.03t3-12.8t2fpi=232.82fi(t)=-6.2t2-51t+465.62si(t)=-2.06t3-25.6t22fpi=465.63fi(t)=-9.3t2-76.5t+698.43si(t)=-3.09t3-38.4t23fpi=698.4滾珠Fb=1.7fb(t)=-1.19t2-9.9t+90sb(t)=-0.4t3-4.95t2fpb=902fb(t)=-2.38t2-19.8t+1802sb(t)=-0.8t3-9.9t22fpb=1803fb(t)=-3.57t2-29.7t+2703sb(t)=-1.2t3-14.85t23fpb=270

圖10中以箭頭標注的3個突出峰的橫坐標分別為234.2、468.2和701.1。后兩個突出峰橫坐標分別是第一個突出峰橫坐標的2倍和3倍的同時，三個突出的橫標值與表4中目標軸承內圈故障所對應的3個相位點相吻合。因此可以判斷出該目標軸承內圈存在缺陷。

圖10 實測信號解調后的頻率譜Fig.10 Frequency spectrum of demodulated signal

最后，對目標軸承滾珠的健康狀況進行判別。如果滾珠產生表面損傷類故障，圖10中橫標分別為90、180和270處將會出現明顯的突出峰，尤其是在頻率為90 Hz處的突出峰最高。然而解調信號頻率譜中的上述3個位置并沒有出現突出峰(由于一階故障特征頻率的峰值具有較高的幅值優勢，因此，圖10中只給出了橫坐標為90 Hz附近的放大圖)。因此可以斷定該目標軸承不存在滾珠故障。

為了進一步揭示新方法的優點，本文利用計算階比分析方法對實測故障軸承信號進行分析作為對比。通過角域重采樣和包絡分析得到階比譜如圖11所示。圖中圓圈標注的峰代表目標軸承外圈故障時對應的階比，箭頭標注的峰代表內圈故障時對應的階比。通過圖10和圖11對比得知：其一，包絡階比譜中代表故障階比的峰聚集性不高，尤其在高階比處表現更為明顯；其二，包絡階比譜中的干擾峰較多，且目標峰幅值優勢不明顯，易造成誤判或漏判。因此本文方法獲得的頻譜圖更簡潔，峰值的辨識度更高，也更為有效的提取故障信息，定位故障點位置。

圖11 實測信號的包絡階比譜Fig.11 Envelope order spectrum of demodulated signal

4 結 論

本文提出了基于迭代廣義解調算法的變轉速滾動軸承復合故障特征提取方法。該方法的優勢主要表現在以下幾個方面：

(1) 利用迭代廣義解調算法可以將多分量信號中的多條特定的時頻曲線轉換成平行于時間軸的直線這一優勢消除了轉速變化對滾動軸承復合故障特征的影響。

(2)利用滾動軸承復合故障信號的時頻特性，即故障信號的包絡時頻譜中存在明顯的且與轉頻成固定規律變化能明確表征軸承故障類型的時頻曲線，轉頻曲線以及故障特征系數確定相位點用以判斷軸承的故障數量及位置。

(3)該算法無需角域重采樣即可完成復合故障特征的提取。另外，信號的處理過程只針對信號中表征軸承故障的特定成分，免除了其他成分的干擾，為變轉速條件下滾動軸承復合故障特征的提取提供了新思路。

需要指出的是本文方法僅適用于目標軸承某一部件上只包含一處故障點的復合故障情況，例如內圈、外圈和滾珠上都存在一處故障點。

[ 2 ] 王曉冬,何正嘉, 訾艷陽.小波自適應構造方法及滾動軸承復合故障診斷研究[J]. 振動工程學報, 2010, 23(4): 438-444.

WANG Xiaodong, HE Zhengjia, ZI Yanyang. Adaptive construction of multiwavelet and research on composite fault diagnosis of rolling bearing[J]. Journal of Vibration Engineering, 2010, 23(4): 438-444.

[ 3 ] LUO Jiesi, YU Dejie, LIANG Ming. Application of multi-scale chirplet path pursuit and fractional Fourier transform for gear fault detection in speed up and speed down processes [J]. Journal of Sound and Vibration, 2012, 331:4971-4986.

[ 4 ] 趙德尊, 李建勇, 程衛東. 變轉速及齒輪噪源干擾下基于IDMM與EMD的滾動軸承故障診斷方法[J]. 振動與沖擊, 2016, 35(10): 101-119.

ZHAO Dezun, LI Jianyong, CHENG Weidong. Method for rolling element bearing fault diagnosis based on IDMM and EMD under time-varying rotational speed and gear noise [J]. Journal of Vibration and Shock, 2016, 35(10): 101-119.

[ 5 ] 明安波,褚福磊,張煒. 滾動軸承復合故障特征分離的小波-頻譜自相關方法[J]. 機械工程學報，2013,49(3): 80-87.

MING Anbo, CHU Fulei, ZHANG Wei. Compound fault features separation of rolling element bearing based on the wavelet decomposition and spectrum auto-correlation [J]. Journal of Mechanical Engineering, 2013, 49(3): 80-87.

[ 6 ] 馬新娜, 楊紹普. 滾動軸承復合故障診斷的自適應方法研究[J]. 振動與沖擊, 2016, 35(10):145-150.

MA Xinna, YANG Shaopu. Adaptive compound fault diagnosis of rolling bearings [J]. Journal of Vibration and Shock, 2016, 35(10): 145-150.

[ 7 ] CHEN Jinglong, ZI Yanyang, HE Zhengjia, et al. Compound faults detection of rotating machinery using improved adaptive redundant lifting multiwavelet [J]. Mechanical Systems and Signal Processing, 2013, 38: 36-54.

[ 8 ] ZHANG Haibin, LU Siliang, HE Qingbo, et al. Multi-bearing defect detection with trackside acoustic signal based on a pseudo time-frequency analysis and Dopplerlet filter [J]. Mechanical Systems and Signal Processing, 2016, 70/71: 176-200.

[ 9 ] FYFE K R, MUNCK E D S. Analysis of computed order tracking [J]. Mechanical System and Signal Processing, 1997, 11(2): 187-205.

[10] RANDALL R B, ANTONI J. Rolling element bearing diagnostics-a tutorial [J]. Mechanical Systems and Signal Processing, 2011, 25 (2): 485-520.

[11] SAAVEDRA P N, RODRIGUEZ C G. Accurate assessment of computed order tracking [J]. Shock and Vibration, 2006, 13 (1): 13-21.

[12] CHENG Weidong, ROBERT X G, WANG Jingjiang, et al. Envelope deformation in computed order tracking and error in order analysis [J]. Mechanical System and Signal Processing, 2014, 48(1/2): 92-102.

[13] OLHEHE S, WALDEN A T. A generalized demodulation approach to time-frequency projections for multicomponent signals [J]. Proceedings of the Royal Society A, 2005, 461(2059): 2159-2179.

[14] CHENG Junsheng, YANG Yu, YU Dejie. Application of the improved generalized demodulation time-frequency analysis method to multi-component signal decomposition [J]. Signal Processing, 2009, 89: 1205-1215.

[15] 黃椒治, 林慧斌, 丁康. 基于廣義解調平滑能量分離算法的瞬時頻率估計[J]. 振動工程學報, 2014, 27(2): 281-288.

HUANG Jiaozhi, LIN Huibin, DING Kang, Instantaneous frequency estimation based on generalized demodulation smoothed energy separation algorithm [J]. Journal of Vibration Engineering, 2014, 27(2): 281-288.

[16] 皮維, 于德介, 彭富強. 基于多尺度線調頻基稀疏信號分解的廣義解調方法及其在齒輪故障診斷中的應用[J]. 機械工程學報, 2010, 46(15): 59-64.

PI Wei, YU Dejie, PENG Fuqiang. Generalized demodulation method based on multi-scale chirplet and sparse signal decomposition and its application to gear fault diagnosis [J]. Journal of Mechanical Engineering, 2010, 46(15): 59-64.

[17] FENG Zhipeng, CHU Fulei, ZUO Mingliang. Time-frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized demodulation [J]. Journal of Sound and Vibration, 2011, 330: 1225-1243.

[18] CHEN Xiaowang, FENG Zhipeng. Iterative generalized time-frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions [J]. Mechanical Systems and Signal Processing, 2016, 80: 429-444.