基于多時頻曲線提取廣義特征的變轉速軸承故障診斷

肖 飛， 張宏立， 馬 萍， 王 聰

(新疆大學 電氣工程學院， 烏魯木齊 830000)

滾動軸承是旋轉機械中的重要部件，由于復雜的工作環境和運行狀態(急劇變轉速、變負載等)，使其在實際生產環境中極易受到損壞。根據相關資料：旋轉機械中的超過30%的機械故障是由軸承故障引起的[1]。而且，在實際生產過程中，機械設備常在變轉速工況下運行。因此研究變轉速工況下的滾動軸承故障診斷具有重要的生產實踐意義。

變轉速工況下，滾動軸承振動信號呈現的時變非平穩特性使得故障特征難以提取。近年來，國內外許多學者針對變轉速工況下的滾動軸承故障診斷開展了大量研究[2-5]。較為典型的方法是基于階次跟蹤的軸承故障診斷技術，其核心思想是將角域上的信號重采樣得到循環平穩的信號，硬件階次跟蹤和計算階次跟蹤都要求設備安裝轉速計來獲取轉速，這增加了空間和經濟成本[6]。因此，國內外學者開始研究無轉速計下的階次跟蹤方法，并提出基于瞬時頻率估計的轉速計算方法[7-8]。此外，階次跟蹤方法還存在一些局限性，其在信號重采樣的過程中存在幅值誤差、包絡畸變以及計算復雜等問題[9-10]。

Olhede等[11]提出的廣義解調方法，是一種可以將非線性非平穩信號轉換為線性平穩信號的方法。劉東東等[12-16]將其應用到旋轉機械的故障診斷中。由于原始廣義解調算法只能對原始信號中的單一分量進行單次解調變換，難以處理多分量信號。因此有學者對廣義解調進行了改進，提出了迭代廣義解調算法[17]。迭代廣義解調算法存在一些不足，由于其對振動信號循環多次進行解調變換，導致結果出現頻譜混疊[18]。

Huang等[19]提出了一種基于快速路徑優化的多時頻曲線提取(multiple time-frequency curve extraction, MTFCE)方法，能夠將包絡時頻圖中的瞬時故障特征頻率(instantaneous fault characteristic frequency, IFCF)和瞬時轉頻(instantaneous shaft rotational frequency, ISRF)曲線提取出來，利用IFCFs與ISRF的比值來描述提取的曲線之間的關系。因此，可通過比值與分量信號特征系數的比較來實現分量信號的特征提取。在變轉速振動信號中，由于低轉速下時頻圖中轉頻和故障沖擊成分被淹沒在噪聲中，導致曲線之間比值誤差較大，使得無法通過比值對故障特征進行有效提取。

綜上所述，為解決變轉速滾動軸承故障信號由于時變非平穩特性、故障沖擊幅值淹沒在噪聲中，導致故障診斷困難的問題，提出一種新的基于多時頻曲線提取廣義特征的變轉速軸承故障診斷方法。利用MTFCE從包絡時頻譜中提取IFCFs和ISRF，基于假設思想，根據廣義解調理論計算變轉速軸承發生故障時的廣義特征指標，再構建量化診斷模型進行故障診斷。該方法相較于迭代廣義解調方法無需對信號進行解調變換，避免了多次解調導致的頻譜混疊問題，同時無需利用曲線的比值進行診斷，提高了診斷的準確率。仿真信號和實例分析結果證明了所提方法的魯棒性和有效性。

1 算法部分

1.1 多時頻曲線提取算法

多時頻曲線提取算法是一種基于快速路徑優化的時頻曲線提取算法，能夠將時頻圖中的IFCFs和ISRF曲線提取出來。

假設信號x(t)經過短時傅里葉變換后的時頻表達(time-frequency representation, TFR)為X(τ,f)，其中τ為時間變量，f為頻率變量。快速路徑優化算法可從X(τ,f)中提取出T-F曲線fp(τ)。定義TFR中時間τn處的峰值表示為Np(τn)，第m個峰值對應頻率為vm(τn)，第m個峰值對應的TFR幅值為Qm(τn)。

(1)

如果X(τ,f)的時間跨度為[τ1,τ2,…,τN]，路徑優化可描述為

(2)

其中mc(τN)確定在時間τn處提取的峰值，F[]為優化選擇支持函數，{m1,…,mN}表示沿時間跨度的峰值數數列。文獻[20]提出了一種只依賴于具有支持函數的有限個前向點的快速路徑優化算法，算法具體步驟如下所示

(3)

其中：

(4)

(5)

m[]=perc0.5[],IQR[]=perc0.75[]-perc0.25[]

(6)

其中fd(τn-1)是τn-1處的候選峰值的頻率，fd是在時間[τ1,…,τn-1]的一系列候選峰值的頻率。Δfd是fd的導數。m[]表示一系列數的中位數，IQR表示四分位距。percp表示序列的第p個分位數。λ1和λ2是懲罰因子，可以取1。權函數w1()和w2()分別用于抑制脊線頻率值和其導數的非典型變化。

快速路徑優化問題可表示為

forn=1,…,N,m=1,…,Np(τn)

andk=1,…,Np(τn-1)

(7)

其中，q(m,τn)用來表示前一時刻τn-1的峰值應與當前時刻τn的峰值相聯系。U(m,τn)是有助于優化的中間向量。

將最后計算結果U(m,τn)最大的峰值作為提取曲線的終點，然后基于q(m,τn)的結果向前追蹤峰值提取出整條曲線。

基于快速路徑的多時頻曲線則是在提取出一條時頻曲線后，將該條曲線從時頻圖中刪去，再重復應用快速路徑優化算法提取下一條時頻曲線，直到提取曲線數量達到預設的值。

1.2 基于假設思想的廣義特征指標計算

變轉速下，滾動軸承的IFCFs在包絡時頻圖中表現出明顯的峰值，且IFCFs與ISRF存在固定的比例關系，這個比例關系與軸承的故障特征系數相關。

廣義解調可通過選擇合適相位函數可以將時頻圖中特定的曲線轉換成平行于時間坐標軸的直線，通過頻譜圖可看到與該曲線對應的突出峰值。該峰值與軸承故障聯系緊密。

因此，基于假設思想，假設發生故障時的該峰值為變轉速軸承故障的廣義特征指標。可利用廣義解調理論定義廣義特征指標。

廣義解調理論的基本思想是對信號x(t)進行廣義傅里葉變換

(8)

式中，s0(t)表示隨時間變化的實值函數。

變轉速滾動軸承發生外圈、內圈和滾珠故障時對應的故障特征系數Fo,Fi和Fb計算如下所示

(9)

(10)

(11)

式中：fo,fi和fb分別表示外圈、內圈和滾珠故障特征頻率；n為滾動軸承滾珠個數；fr為軸承轉頻；D為滾動軸承節圓直徑；d代表滾珠直徑；α為接觸角。

假設預設轉頻方程為

y=c0+c1t+c2t2+…+cmtn

(12)

由廣義解調變換理論，提出解調變換后的故障廣義特征指標如表1所示。

表1 廣義特征指標計算Tab.1 Generalized feature calculation

1.3 量化診斷模型建立

由廣義解調算法可知信號x(t)的時頻分布是一條由f(t)確定的曲線。而MTFCE算法可從時頻圖中提取出f(t)。因此可以不通過解調變換直接從時頻圖中獲取解調后的f0，再由f0進行故障診斷。

由MTFCE從時頻圖中提取出ISRF和IFCFs曲線，截取平滑段對其進行多項式擬合得到方程如下所示

yISRF=a0+a1t+a2t2+…+antn

(13)

yIFCF(k)=bk0+bk1t+bk2t2+…+bkntn

(14)

n=1,2,…，N且a0

由于故障特征頻率趨勢線在包絡信號的時頻圖中表現為具有特定分布規律的曲線。MTFCE可以從時頻圖中直接提取出轉頻曲線。再經由多項式擬合獲取轉頻方程，不需要額外安裝轉速計獲取轉速信息。擬合轉頻yISRF和預設轉頻y吻合時則表示該方法估計轉頻的有效性。

由于時頻分辨率的限制，MTFCE提取的能夠完整表征故障特征的曲線一般為3～4條。假設提取曲線為ISRF和IFCF1和IFCF2。結合表1計算故障廣義特征頻率，可構建量化診斷模型如圖1所示。

圖1 量化診斷模型Fig.1 Quantitative diagnostic model

將b10,b20與假設故障廣義特征對比確定故障類型。為了方便故障診斷以及算法對比，繪制以bl0為頻率的固定幅值頻譜圖結合量化診斷模型來進行故障診斷。

2 故障診斷

基于上述分析，本文提出一種基于多時頻曲線提取廣義特征的變轉速軸承故障診斷方法。該方法流程圖如圖2所示。具體算法步驟如下所示：

圖2 故障診斷方法流程圖Fig.2 Flow chart of fault diagnosis method

(1) 利用快速譜峭度算法計算最優帶通濾波參數，并對原始振動信號進行濾波。

(2) 對濾波后的共振帶信號進行短時傅里葉變換，并采用多時頻曲線提取算法提取IFCFs和ISRF。

(3) 截取平滑曲線擬合IFCFs和ISRF曲線的多項式方程，再由擬合后的ISRF曲線參數結合表1計算廣義特征。

(4) 由擬合后的IFCFs曲線參數繪制固定初始值頻譜圖，再結合圖1量化診斷模型進行故障診斷。

3 仿真分析

為驗證提出算法有效性，構造了升速條件下軸承故障仿真信號[21]。

xb(t)={1+A′(t)sin[2πf(t)·t]}·

μ(t-tm)

(15)

式中：A′(t) 為轉頻的瞬時幅值；f(t)為軸承轉頻；N為故障信號的沖擊個數；Am為第m個故障沖擊幅值；β為軸承衰減系數；ωr為故障的共振頻率；μ(t)為單位階躍函數；tm為第m個故障沖擊發生時刻。其計算公式如下所示

(16)

其中τ為滾珠滑移引起的軸承故障沖擊間隔誤差，其取值通常為0.01～0.02。仿真信號具體參數如表2所示。

表2 故障軸承仿真信號參數Tab.2 Parameters of simulated faulty bearing signal

根據上述仿真模型得到的故障軸承振動信號如圖3所示，對仿真信號進行快速譜峭度濾波得到共振帶信號，濾波參數由圖4的中心頻率和帶寬確定。然后利用短時傅立葉變換得到共振帶信號的包絡時頻圖如圖5所示。可以從時頻圖中得到清晰的IFCFs和ISRF，而且可以發現初始階段，轉頻較小時，故障特征的相關故障沖擊幅值被噪聲淹沒。

圖3 故障滾動軸承仿真信號Fig.3 Faulty rolling bearing simulation signal

圖4 仿真信號Kurtogram圖Fig.4 Kurtogram of simulated signal

圖5 仿真信號包絡時頻圖Fig.5 Envelope time-frequency diagram of simulated signal

利用多時頻曲線提取方法提取ISRF和IFCFs曲線，截取平滑的T-F (time-frequency)曲線如圖6所示，圖中存在ISRF以及IFCFs以及IFCF1與ISRF的耦合曲線。對T-F曲線進行多項式擬合，擬合后的曲線參數如表3所示。從表中看出提取出的轉頻曲線的擬合方程f0=-1.228t2+12.884+23.589與預設的轉頻曲線的擬合方程基本一致。

圖6 仿真信號截取的T-F曲線Fig.6 T-F curve after intercepting of simulated signal

表3 擬合曲線參數Tab.3 Parameters of the fitted curve

假設內圈故障系數Fi=5.4，外圈故障系數Fo=3.5，滾珠故障系數Fb=2.3。利用廣義解調理論結合表1計算的仿真信號的廣義特征指標如表4所示。由擬合曲線參數繪制固定幅值的初始值頻譜圖如圖7所示，圖中外圈故障特征頻率fo及其諧波2fo被有效提取，可以確定軸承存在外圈故障。

表4 仿真信號廣義特征指標計算Tab.4 Calculation of generalized features of simulated signals

圖7 仿真信號固定幅值的頻譜圖Fig.7 Spectrogram of the simulated signal with fixed amplitude

為了驗證本文所提方法的有效性，利用迭代廣義解調算法，根據故障系數所確定的相位函數對信號循環進行解調變換。所得頻譜圖如圖8所示。圖中可以看出，頻率峰值與軸承外圈故障特征頻率吻合，但從圈出的橢圓部分可知，迭代廣義解調算法存在頻譜混疊現象，會造成故障誤診。本文方法相對于迭代廣義解調算法表現更豐富故障相關信息，例如轉頻和轉頻與故障特征的耦合。體現了基于多時頻曲線提取廣義特征的變轉速軸承故障診斷方法的有效性和優越性。

圖8 迭代廣義解調頻譜圖Fig.8 Spectrogram of iterative generalized demodulation

4 實例分析

利用加拿大渥太華大學采集到的變轉速滾動軸承的振動信號對算法的有效性進行進一步的驗證。采用型號為MFS-PK5M的機械故障模擬器所采集的變轉速工況下軸承數據進行驗證，該數據采集頻率為200 000 Hz，采集時間10 s，測試軸承的相關參數如表5所示，試驗平臺如圖9所示，實例振動信號如圖10所示。

表5 滾動軸承參數Tab.5 Rolling bearing parameters

圖9 試驗平臺Fig.9 Test set-up

圖10 實例振動信號Fig.10 Real faulty bearing signal

由滾動軸承參數可以計算出，內圈故障系數Fi=5.43，外圈故障系數Fo=3.57，滾珠故障系數Fb=2.32。利用快速譜峭度濾波算法對振動信號濾波得到共振帶信號，濾波參數由圖11的中心頻率和帶寬確定。對共振帶信號的包絡信號進行短時傅立葉變換得到時頻圖如圖12所示。

圖11 實例信號Kurtogram圖Fig.11 Kurtogram of example signal

圖12 實例信號包絡時頻圖Fig.12 Envelope time-frequency diagram of example signal

利用MTFCE提取時頻曲線，圖13是截取平滑段以后的T-F曲線，對曲線進行多項式擬合。擬合后的數據如表6所示，轉頻方程為1.500t+12.568。由于該數據集提供了轉速信息。圖14顯示了轉頻對比，從圖中可以看出擬合后的轉頻曲線與實際測量轉頻誤差很小，而從時頻圖中提取的轉頻與測量轉頻相差很大。

圖13 實例信號截取以后的T-F曲線Fig.13 T-F curve after interception of example signal

表6 擬合曲線參數Tab.6 Parameters of the fitted curve

圖14 轉頻對比圖Fig.14 RPM comparison chart

由擬合后的ISRF和IFCFs繪制固定幅值的頻譜圖如圖15所示。再根據表1計算的實例信號的廣義特征如表7所示。結合量化診斷模型可以確定故障類型為內圈故障。圖16為迭代廣義解調頻譜圖，可以看出低頻部分由于頻譜混疊導致會存在故障的誤診。

圖15 實例信號固定幅值頻譜圖Fig.15 Spectrogram of fixed amplitude of example signal

圖16 實例信號迭代廣義解調頻譜圖Fig.16 Example signal spectrogram of IGD

5 結 論

針對低變轉速下滾動軸承故障振動信號呈現時變非平穩性，被噪聲干擾無法有效診斷的問題，提出一種基于多時頻曲線提取廣義特征的變轉速滾動軸承故障診斷方法，通過仿真和實例分析驗證了所提方法的有效性。得出了以下結論：

(1) 相對于迭代廣義解調算法，本文方法無需對信號進行后續的解調變換，避免了頻譜混疊問題。

(2) 通過基于假設思想的廣義特征計算和量化診斷模型的構建，所提方法可對變轉速工況下滾動軸承故障特征進行提取，有效的實現了變轉速滾動軸承故障診斷，且魯棒性較高。

(3) 由于海森堡不確定原理，短時傅里葉變換不能同時得到高的時間分辨率和頻率分辨率。而MTFCE從時頻圖中提取IFCFs和ISRF，時頻分辨率越高提取效果越好。后續將會提高時頻分辨率上做研究，進一步提高曲線提取的準確性。