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

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

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

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

Olhede等[11]提出的廣義解調(diào)方法，是一種可以將非線性非平穩(wěn)信號(hào)轉(zhuǎn)換為線性平穩(wěn)信號(hào)的方法。劉東東等[12-16]將其應(yīng)用到旋轉(zhuǎn)機(jī)械的故障診斷中。由于原始廣義解調(diào)算法只能對(duì)原始信號(hào)中的單一分量進(jìn)行單次解調(diào)變換，難以處理多分量信號(hào)。因此有學(xué)者對(duì)廣義解調(diào)進(jìn)行了改進(jìn)，提出了迭代廣義解調(diào)算法[17]。迭代廣義解調(diào)算法存在一些不足，由于其對(duì)振動(dòng)信號(hào)循環(huán)多次進(jìn)行解調(diào)變換，導(dǎo)致結(jié)果出現(xiàn)頻譜混疊[18]。

Huang等[19]提出了一種基于快速路徑優(yōu)化的多時(shí)頻曲線提取(multiple time-frequency curve extraction, MTFCE)方法，能夠?qū)j(luò)時(shí)頻圖中的瞬時(shí)故障特征頻率(instantaneous fault characteristic frequency, IFCF)和瞬時(shí)轉(zhuǎn)頻(instantaneous shaft rotational frequency, ISRF)曲線提取出來(lái)，利用IFCFs與ISRF的比值來(lái)描述提取的曲線之間的關(guān)系。因此，可通過(guò)比值與分量信號(hào)特征系數(shù)的比較來(lái)實(shí)現(xiàn)分量信號(hào)的特征提取。在變轉(zhuǎn)速振動(dòng)信號(hào)中，由于低轉(zhuǎn)速下時(shí)頻圖中轉(zhuǎn)頻和故障沖擊成分被淹沒(méi)在噪聲中，導(dǎo)致曲線之間比值誤差較大，使得無(wú)法通過(guò)比值對(duì)故障特征進(jìn)行有效提取。

綜上所述，為解決變轉(zhuǎn)速滾動(dòng)軸承故障信號(hào)由于時(shí)變非平穩(wěn)特性、故障沖擊幅值淹沒(méi)在噪聲中，導(dǎo)致故障診斷困難的問(wèn)題，提出一種新的基于多時(shí)頻曲線提取廣義特征的變轉(zhuǎn)速軸承故障診斷方法。利用MTFCE從包絡(luò)時(shí)頻譜中提取IFCFs和ISRF，基于假設(shè)思想，根據(jù)廣義解調(diào)理論計(jì)算變轉(zhuǎn)速軸承發(fā)生故障時(shí)的廣義特征指標(biāo)，再構(gòu)建量化診斷模型進(jìn)行故障診斷。該方法相較于迭代廣義解調(diào)方法無(wú)需對(duì)信號(hào)進(jìn)行解調(diào)變換，避免了多次解調(diào)導(dǎo)致的頻譜混疊問(wèn)題，同時(shí)無(wú)需利用曲線的比值進(jìn)行診斷，提高了診斷的準(zhǔn)確率。仿真信號(hào)和實(shí)例分析結(jié)果證明了所提方法的魯棒性和有效性。

1 算法部分

1.1 多時(shí)頻曲線提取算法

多時(shí)頻曲線提取算法是一種基于快速路徑優(yōu)化的時(shí)頻曲線提取算法，能夠?qū)r(shí)頻圖中的IFCFs和ISRF曲線提取出來(lái)。

假設(shè)信號(hào)x(t)經(jīng)過(guò)短時(shí)傅里葉變換后的時(shí)頻表達(dá)(time-frequency representation, TFR)為X(τ,f)，其中τ為時(shí)間變量，f為頻率變量。快速路徑優(yōu)化算法可從X(τ,f)中提取出T-F曲線fp(τ)。定義TFR中時(shí)間τn處的峰值表示為Np(τn)，第m個(gè)峰值對(duì)應(yīng)頻率為vm(τn)，第m個(gè)峰值對(duì)應(yīng)的TFR幅值為Qm(τn)。

(1)

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

(2)

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

(3)

其中：

(4)

(5)

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

(6)

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

快速路徑優(yōu)化問(wèn)題可表示為

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

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

(7)

其中，q(m,τn)用來(lái)表示前一時(shí)刻τn-1的峰值應(yīng)與當(dāng)前時(shí)刻τn的峰值相聯(lián)系。U(m,τn)是有助于優(yōu)化的中間向量。

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

基于快速路徑的多時(shí)頻曲線則是在提取出一條時(shí)頻曲線后，將該條曲線從時(shí)頻圖中刪去，再重復(fù)應(yīng)用快速路徑優(yōu)化算法提取下一條時(shí)頻曲線，直到提取曲線數(shù)量達(dá)到預(yù)設(shè)的值。

1.2 基于假設(shè)思想的廣義特征指標(biāo)計(jì)算

變轉(zhuǎn)速下，滾動(dòng)軸承的IFCFs在包絡(luò)時(shí)頻圖中表現(xiàn)出明顯的峰值，且IFCFs與ISRF存在固定的比例關(guān)系，這個(gè)比例關(guān)系與軸承的故障特征系數(shù)相關(guān)。

廣義解調(diào)可通過(guò)選擇合適相位函數(shù)可以將時(shí)頻圖中特定的曲線轉(zhuǎn)換成平行于時(shí)間坐標(biāo)軸的直線，通過(guò)頻譜圖可看到與該曲線對(duì)應(yīng)的突出峰值。該峰值與軸承故障聯(lián)系緊密。

因此，基于假設(shè)思想，假設(shè)發(fā)生故障時(shí)的該峰值為變轉(zhuǎn)速軸承故障的廣義特征指標(biāo)?？衫脧V義解調(diào)理論定義廣義特征指標(biāo)。

廣義解調(diào)理論的基本思想是對(duì)信號(hào)x(t)進(jìn)行廣義傅里葉變換

(8)

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

變轉(zhuǎn)速滾動(dòng)軸承發(fā)生外圈、內(nèi)圈和滾珠故障時(shí)對(duì)應(yīng)的故障特征系數(shù)Fo,Fi和Fb計(jì)算如下所示

(9)

(10)

(11)

式中：fo,fi和fb分別表示外圈、內(nèi)圈和滾珠故障特征頻率；n為滾動(dòng)軸承滾珠個(gè)數(shù)；fr為軸承轉(zhuǎn)頻；D為滾動(dòng)軸承節(jié)圓直徑；d代表滾珠直徑；α為接觸角。

假設(shè)預(yù)設(shè)轉(zhuǎn)頻方程為

y=c0+c1t+c2t2+…+cmtn

(12)

由廣義解調(diào)變換理論，提出解調(diào)變換后的故障廣義特征指標(biāo)如表1所示。

表1 廣義特征指標(biāo)計(jì)算Tab.1 Generalized feature calculation

1.3 量化診斷模型建立

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

由MTFCE從時(shí)頻圖中提取出ISRF和IFCFs曲線，截取平滑段對(duì)其進(jìn)行多項(xiàng)式擬合得到方程如下所示

yISRF=a0+a1t+a2t2+…+antn

(13)

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

(14)

n=1,2,…，N且a0

由于故障特征頻率趨勢(shì)線在包絡(luò)信號(hào)的時(shí)頻圖中表現(xiàn)為具有特定分布規(guī)律的曲線。MTFCE可以從時(shí)頻圖中直接提取出轉(zhuǎn)頻曲線。再經(jīng)由多項(xiàng)式擬合獲取轉(zhuǎn)頻方程，不需要額外安裝轉(zhuǎn)速計(jì)獲取轉(zhuǎn)速信息。擬合轉(zhuǎn)頻yISRF和預(yù)設(shè)轉(zhuǎn)頻y吻合時(shí)則表示該方法估計(jì)轉(zhuǎn)頻的有效性。

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

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

將b10,b20與假設(shè)故障廣義特征對(duì)比確定故障類型。為了方便故障診斷以及算法對(duì)比，繪制以bl0為頻率的固定幅值頻譜圖結(jié)合量化診斷模型來(lái)進(jìn)行故障診斷。

2 故障診斷

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

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

(1) 利用快速譜峭度算法計(jì)算最優(yōu)帶通濾波參數(shù)，并對(duì)原始振動(dòng)信號(hào)進(jìn)行濾波。

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

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

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

3 仿真分析

為驗(yàn)證提出算法有效性，構(gòu)造了升速條件下軸承故障仿真信號(hào)[21]。

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

μ(t-tm)

(15)

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

(16)

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

表2 故障軸承仿真信號(hào)參數(shù)Tab.2 Parameters of simulated faulty bearing signal

根據(jù)上述仿真模型得到的故障軸承振動(dòng)信號(hào)如圖3所示，對(duì)仿真信號(hào)進(jìn)行快速譜峭度濾波得到共振帶信號(hào)，濾波參數(shù)由圖4的中心頻率和帶寬確定。然后利用短時(shí)傅立葉變換得到共振帶信號(hào)的包絡(luò)時(shí)頻圖如圖5所示。可以從時(shí)頻圖中得到清晰的IFCFs和ISRF，而且可以發(fā)現(xiàn)初始階段，轉(zhuǎn)頻較小時(shí)，故障特征的相關(guān)故障沖擊幅值被噪聲淹沒(méi)。

圖3 故障滾動(dòng)軸承仿真信號(hào)Fig.3 Faulty rolling bearing simulation signal

圖4 仿真信號(hào)Kurtogram圖Fig.4 Kurtogram of simulated signal

圖5 仿真信號(hào)包絡(luò)時(shí)頻圖Fig.5 Envelope time-frequency diagram of simulated signal

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

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

表3 擬合曲線參數(shù)Tab.3 Parameters of the fitted curve

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

表4 仿真信號(hào)廣義特征指標(biāo)計(jì)算Tab.4 Calculation of generalized features of simulated signals

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

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

圖8 迭代廣義解調(diào)頻譜圖Fig.8 Spectrogram of iterative generalized demodulation

4 實(shí)例分析

利用加拿大渥太華大學(xué)采集到的變轉(zhuǎn)速滾動(dòng)軸承的振動(dòng)信號(hào)對(duì)算法的有效性進(jìn)行進(jìn)一步的驗(yàn)證。采用型號(hào)為MFS-PK5M的機(jī)械故障模擬器所采集的變轉(zhuǎn)速工況下軸承數(shù)據(jù)進(jìn)行驗(yàn)證，該數(shù)據(jù)采集頻率為200 000 Hz，采集時(shí)間10 s，測(cè)試軸承的相關(guān)參數(shù)如表5所示，試驗(yàn)平臺(tái)如圖9所示，實(shí)例振動(dòng)信號(hào)如圖10所示。

表5 滾動(dòng)軸承參數(shù)Tab.5 Rolling bearing parameters

圖9 試驗(yàn)平臺(tái)Fig.9 Test set-up

圖10 實(shí)例振動(dòng)信號(hào)Fig.10 Real faulty bearing signal

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

圖11 實(shí)例信號(hào)Kurtogram圖Fig.11 Kurtogram of example signal

圖12 實(shí)例信號(hào)包絡(luò)時(shí)頻圖Fig.12 Envelope time-frequency diagram of example signal

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

圖13 實(shí)例信號(hào)截取以后的T-F曲線Fig.13 T-F curve after interception of example signal

表6 擬合曲線參數(shù)Tab.6 Parameters of the fitted curve

圖14 轉(zhuǎn)頻對(duì)比圖Fig.14 RPM comparison chart

由擬合后的ISRF和IFCFs繪制固定幅值的頻譜圖如圖15所示。再根據(jù)表1計(jì)算的實(shí)例信號(hào)的廣義特征如表7所示。結(jié)合量化診斷模型可以確定故障類型為內(nèi)圈故障。圖16為迭代廣義解調(diào)頻譜圖，可以看出低頻部分由于頻譜混疊導(dǎo)致會(huì)存在故障的誤診。

圖15 實(shí)例信號(hào)固定幅值頻譜圖Fig.15 Spectrogram of fixed amplitude of example signal

圖16 實(shí)例信號(hào)迭代廣義解調(diào)頻譜圖Fig.16 Example signal spectrogram of IGD

5 結(jié) 論

針對(duì)低變轉(zhuǎn)速下滾動(dòng)軸承故障振動(dòng)信號(hào)呈現(xiàn)時(shí)變非平穩(wěn)性，被噪聲干擾無(wú)法有效診斷的問(wèn)題，提出一種基于多時(shí)頻曲線提取廣義特征的變轉(zhuǎn)速滾動(dòng)軸承故障診斷方法，通過(guò)仿真和實(shí)例分析驗(yàn)證了所提方法的有效性。得出了以下結(jié)論：

(1) 相對(duì)于迭代廣義解調(diào)算法，本文方法無(wú)需對(duì)信號(hào)進(jìn)行后續(xù)的解調(diào)變換，避免了頻譜混疊問(wèn)題。

(2) 通過(guò)基于假設(shè)思想的廣義特征計(jì)算和量化診斷模型的構(gòu)建，所提方法可對(duì)變轉(zhuǎn)速工況下滾動(dòng)軸承故障特征進(jìn)行提取，有效的實(shí)現(xiàn)了變轉(zhuǎn)速滾動(dòng)軸承故障診斷，且魯棒性較高。

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