李書樂，馬潔

(1. 北京信息科技大學(xué) 自動化學(xué)院，北京 100192；2. 北京信息科技大學(xué) 機(jī)電學(xué)院，北京 100192)

滾動軸承是機(jī)械設(shè)備中的重要部件,也是極易出現(xiàn)故障的部件。由于振動信號易于監(jiān)測和包含特征信息豐富等優(yōu)點(diǎn),70%的故障診斷都是通過振動分析實(shí)現(xiàn)的[1-2],短時(shí)傅里葉變換、小波分析和經(jīng)驗(yàn)?zāi)Ｊ椒纸獾葌鹘y(tǒng)的時(shí)頻分析方法可以有效地處理振動信號,如唐廣先等[3]用基于短時(shí)傅里葉變換等方法對滾動軸承進(jìn)行故障診斷。李恒等[4]提出的基于短時(shí)傅里葉變換和卷積神經(jīng)網(wǎng)絡(luò)的故障診斷方法實(shí)現(xiàn)了端到端的故障模式識別。鄭勇[5]用基于小波和能量特征提取的方法進(jìn)行旋轉(zhuǎn)機(jī)械故障診斷。張琛等[6]實(shí)現(xiàn)了基于小波包時(shí)延相關(guān)解調(diào)的滾動軸承故障診斷。Chhaya G等[7]用基于經(jīng)驗(yàn)?zāi)B(tài)分解和Hjorth參數(shù)的方法進(jìn)行滾動軸承故障診斷。

2014年Dragomiretskiy和Zosso提出了一種新的時(shí)頻分析方法——變分模態(tài)分解(Variational mode decomposition,VMD),克服了傳統(tǒng)時(shí)頻分析方法和經(jīng)驗(yàn)?zāi)B(tài)分解等自適應(yīng)信號分解方法的缺陷,廣泛應(yīng)用于故障診斷領(lǐng)域,如張俊等[8]用基于PSO-VMD-MCKD的方法對風(fēng)機(jī)軸承進(jìn)行故障診斷。

然而,這些基于振動測試的滾動軸承故障診斷是在假設(shè)轉(zhuǎn)速恒定的條件下進(jìn)行的,實(shí)際環(huán)境中的滾動軸承通常是在變轉(zhuǎn)速工況下運(yùn)行的,在設(shè)備的啟、停過程,電機(jī)電壓不穩(wěn)定,轉(zhuǎn)子負(fù)載變化時(shí)都會引起滾動軸承轉(zhuǎn)速的波動,故其故障振動信號為時(shí)變非平穩(wěn)信號,其故障特征頻率是變化的,而現(xiàn)有的時(shí)頻分析方法在分析時(shí)變多分量信號時(shí)存在模態(tài)混疊現(xiàn)象,不能有效處理這類信號,再加上在實(shí)際運(yùn)行環(huán)境中故障源產(chǎn)生強(qiáng)大背景噪聲的同時(shí),其他的機(jī)械設(shè)備也會產(chǎn)生環(huán)境噪聲,所以故障信號極易淹沒在噪聲中,而傳統(tǒng)的故障診斷方法往往會造成漏診斷現(xiàn)象,本文應(yīng)用基于廣義變分模態(tài)分解(GVMD)和分?jǐn)?shù)階傅里葉變換(FRFT)的方法處理時(shí)變非平穩(wěn)信號,實(shí)現(xiàn)對變轉(zhuǎn)速工況下滾動軸承故障的精細(xì)診斷。

1 廣義變分模態(tài)分解

1.1 廣義解調(diào)

廣義解調(diào)是一種新型的信號處理工具,尤其適用于多分量非平穩(wěn)信號的分離。廣義解調(diào)的核心在于對信號進(jìn)行廣義傅里葉變換,將呈曲線分布的時(shí)頻特性曲線轉(zhuǎn)換為平行于時(shí)間軸的直線,從而避免信號分量在時(shí)頻域的重疊[9]。

對于信號s(t),其廣義傅里葉變換定義為

(1)

式中：d(t)為一個(gè)實(shí)值函數(shù),exp[-j2πd(t)]為核函數(shù)。式(1)等同于對s(t)exp[-j2πd(t)]進(jìn)行標(biāo)準(zhǔn)傅里葉變換,同理定義逆廣義傅里葉變換為

(2)

由此可得,如果令SG(f)=δ(f-f0),則有s(t)=exp{j2π[f0t+d(t)]},所以,假設(shè)一個(gè)時(shí)頻分布為f0+d′(t)的信號s(t),只要找到一個(gè)近似d(t)的相位函數(shù)對原始信號進(jìn)行廣義解調(diào),得到的解調(diào)函數(shù)s(t)exp[-j2πd(t)]的時(shí)頻分布將是一條平行于時(shí)間軸的直線f(t)=f0。

1.2 變分模態(tài)分解

變分模態(tài)分解(VMD)是自適應(yīng)求解約束變分問題的方法,通過循環(huán)迭代來求解不同的頻率中心及帶寬,實(shí)現(xiàn)對信號分解[10]。

VMD是一個(gè)含約束的最優(yōu)化問題

(3)

式中：δ(t)為脈沖函數(shù);uk(t)為分解得到的k個(gè)本征模態(tài)函數(shù)(IMF)分量;ωk為第k個(gè)分量的中心頻率;f為原始信號。

為了求解上述的變分問題,引入拉格朗日算子λ(t)和二次懲罰因子α,把約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題,其表達(dá)式為

(4)

通過乘法算子交替方向法(Alternate direction method of multipliers, ADMM)交替更新的求解過程如下:

2) 循環(huán)條件n=n+1;

3) 滿足ω≥0時(shí),更新uk；

(5)

4) 更新ωk(t)；

(6)

5) 更新λ；

(7)

式中τ為噪聲的容限參數(shù)。

6) 循環(huán)步驟2)～步驟5),直到滿足迭代條件(式(8))時(shí)停止迭代。

(8)

1.3 廣義變分模態(tài)分解算法流程

廣義變分模態(tài)分解是通過將時(shí)頻特性呈曲線變化的故障特征信號轉(zhuǎn)化為時(shí)頻特性近似線性變化的信號以實(shí)現(xiàn)變轉(zhuǎn)速下軸承的故障診斷,適用于時(shí)變多分量信號的處理,具體步驟如下[11]:

1) 通過小波變換預(yù)估計(jì)信號x(t)中包含的模態(tài)個(gè)數(shù)K。

2) 對實(shí)信號x(t)進(jìn)行希爾伯特變換,得到解析信號y(t)=x(t)+jH[x(t)],j為虛數(shù)單位。

3) 對實(shí)信號x(t)進(jìn)行希爾伯特黃變換,求得信號的相位函數(shù)dk(t),k=1,2,…,K。

4) 選取步驟3)中沖擊性最強(qiáng)分量的相位函數(shù)dk(t),對y(t)進(jìn)行廣義傅里葉變換,即Z(t)=y(t)·exp[-j2πdk(t)]。

2 分?jǐn)?shù)階傅里葉變換

2.1 基于分?jǐn)?shù)階傅里葉變換的濾波原理

滾動軸承的故障特征頻率是非平穩(wěn)信號,且實(shí)際工況下,常常伴隨強(qiáng)噪聲信號的干擾,分?jǐn)?shù)階傅里葉變換(FRFT)是處理非平穩(wěn)信號的有效方法,是一種將信號從時(shí)域變換到分?jǐn)?shù)階域的精細(xì)濾波方法,由于非平穩(wěn)信號在某個(gè)分?jǐn)?shù)階域具有能量聚集的特性,表現(xiàn)出明顯的峰值,而噪聲信號不具備該特性且其能量會均勻分布在任何的分?jǐn)?shù)階域上,故只要選擇好合適的階次,就能實(shí)現(xiàn)精細(xì)地濾波。如文獻(xiàn)[12]將分?jǐn)?shù)階傅里葉變換應(yīng)用在軸承的故障診斷中。

定義在u′域的函數(shù)f(u′)的p階分?jǐn)?shù)階Fourier[13]是

(9)

如圖1所示,非平穩(wěn)信號的一個(gè)分量的時(shí)頻分布與時(shí)間軸的夾角為β,對該分量進(jìn)行分?jǐn)?shù)階傅里葉變換,當(dāng)其旋轉(zhuǎn)角度α與β正交時(shí),在u0點(diǎn)處則會表現(xiàn)為一個(gè)沖擊函數(shù),而噪聲信號將均勻地分布在整個(gè)時(shí)域平面。以分?jǐn)?shù)階域聚焦點(diǎn)u0為中心做窄帶濾波,再做FRFT逆變換就可以得到濾波后的時(shí)域信號。

圖1 FRFT提取信號分量的物理意義圖

由圖1可以得出,最佳的旋轉(zhuǎn)角度α、階數(shù)p與調(diào)頻率fm的關(guān)系為

(10)

2.2 分?jǐn)?shù)階傅里葉變換的濾波算法

FRFT對含有噪聲的信號濾波過程一般是建立在二維搜索基礎(chǔ)上的,通過自由參數(shù)的變換改變旋轉(zhuǎn)角度α,在二維平面(α,u)中進(jìn)行二維搜索得到參數(shù)的估計(jì)值[14]。計(jì)算步驟如下:

1) 對故障信號采用離散采樣法FRFT,通過峰值搜索的方法得到最大值對應(yīng)的α0和u0;

2) 對信號進(jìn)行p0(p0=2α0/π)階分?jǐn)?shù)階傅里葉變換,經(jīng)過變換所得到的信號為

Xp0(u)=Fp0(u)+Np0(u)

(11)

3) 在u域?qū)π纬傻募夥暹M(jìn)行遮隔處理,得

(12)

該理想帶通濾波器以u0為濾波中心,選擇合適的帶寬即可濾除掉大部分噪聲。

4) 對濾波后的信號進(jìn)分?jǐn)?shù)階傅里葉逆變換,得到了消除噪聲的信號。

3 仿真信號分析

3.1 GVMD仿真分析

為了驗(yàn)證GVMD算法可以有效處理時(shí)變多模態(tài)的信號,對復(fù)合信號x(t)=x1(t)+x2(t)+x3(t)進(jìn)行GVMD,其中：

x1(t)=3.5exp(-t)cos[2π(16t3-4t2+25t)] (13)

x2(t)=4cos[2π(22t3-8t2+50t)]

(14)

x3(t)=(3+2t)cos[2π(23t3-12t2+80t)]

(15)

t∈[0,1],其時(shí)域波形如圖2所示。

圖2 x(t)及3個(gè)成分的時(shí)域波形

首先,對原始信號進(jìn)行小波變換確定其預(yù)分解模態(tài)個(gè)數(shù),變換后的時(shí)頻圖如圖3所示。

圖3 小波變換時(shí)頻圖

由圖3可清晰的看出原始信號被分解為3個(gè)分量,即K=3,這與所給仿真信號的分量個(gè)數(shù)相符,說明了小波變換可以用來確定預(yù)分解模態(tài)的個(gè)數(shù)。首先,對原始信號進(jìn)行希爾伯特黃變換得到相位函數(shù),選取沖擊性最強(qiáng)分量的相位函數(shù)d(t)=22t3-10t2對解析信號進(jìn)行廣義傅里葉變換。其次,對廣義傅里葉變換后的信號進(jìn)行VMD分解,懲罰因子a=2 000,分解后信號的時(shí)域圖如圖4a)所示。由此可以看出分解后的分量c1、c2、c3分別對應(yīng)于真實(shí)分量x1(t)、x2(t)、x3(t),虛線為真實(shí)分量。再次,為了對比GVMD算法的有效性,采用VMD算法對x(t)進(jìn)行處理,其分解結(jié)果如圖4b)所示。

圖4 x(t)的2種分解結(jié)果

因此,由圖4可以看出,GVMD算法的分解結(jié)果與真實(shí)分量非常接近,證明該算法可以有效地處理包含時(shí)變多模態(tài)的信號,而VMD算法并不適用。

3.2 FRFT仿真分析

為了驗(yàn)證FRFT算法可以有效降低其他分量和噪聲的干擾,對帶噪聲的仿真信號x(t)=2exp(jπkt2+j2πft)+n(t)進(jìn)行FRFT濾波。初始頻率f為8 Hz,調(diào)頻率k為8 Hz,采樣點(diǎn)數(shù)N為1 025,采樣頻率fs為512 Hz,n(t)為高斯白噪聲,信噪比為-15 dB,仿真信號的時(shí)域圖如圖5所示。分?jǐn)?shù)階傅里葉變換采用采樣型算法[15]進(jìn)行分析,從p=0開始,以p=0.01為步長,得到以p為x軸,以采樣點(diǎn)數(shù)u為y軸,以FRFT幅值為z軸的分?jǐn)?shù)階傅里葉域譜,如圖6所示,由圖可以看出在p=1.02處出現(xiàn)了峰值,峰值為73.77,再以u=529為中心濾波,對濾波后的信號逆變換得到的信號時(shí)域圖如圖7所示,與圖3仿真信號的時(shí)域圖對比可知FRFT有良好的的濾波效果。

圖5 仿真信號時(shí)域圖

圖6 分?jǐn)?shù)階傅里葉變換后的故障診斷分析圖

圖7 濾波后信號的時(shí)域圖

4 變轉(zhuǎn)速工況下滾動軸承故障特征提取的試驗(yàn)驗(yàn)證

4.1 試驗(yàn)數(shù)據(jù)介紹

在本節(jié)中,將用試驗(yàn)采集到的故障數(shù)據(jù)證明該方法的有效性。試驗(yàn)數(shù)據(jù)來自渥太華大學(xué)實(shí)驗(yàn)室的旋轉(zhuǎn)機(jī)械故障模擬器MKF-PK5M。試驗(yàn)臺如圖8所示[16],主軸由一個(gè)交流逆變器控制的電動機(jī)驅(qū)動,兩端分別由兩個(gè)軸承支撐。兩個(gè)軸承中間裝有負(fù)載,故障軸承位于左側(cè)。在試驗(yàn)臺上安裝轉(zhuǎn)速計(jì)和加速度計(jì),分別采集軸轉(zhuǎn)速和振動信號。數(shù)據(jù)經(jīng)信號調(diào)節(jié)器(PCB PIEZOTRONIC 482C15)處理后,送入NI數(shù)據(jù)采集模塊(NI USB-6212 BNC),由計(jì)算機(jī)軟件Labview記錄。本文采用轉(zhuǎn)速從1 092 r/min變化至2 100 r/min時(shí)的滾動軸承內(nèi)圈故障數(shù)據(jù),表1所示為所用軸承的具體參數(shù)情況。

圖8 旋轉(zhuǎn)機(jī)械故障試驗(yàn)臺

表1 內(nèi)圈故障軸承參數(shù)

4.2 GVMD-FRFT滾動軸承故障特征提取算法

本文所提出的 GVMD-FRFT算法不僅可以有效地抑制噪聲干擾,而且可以解決變工況情況下信號處理困難的問題,能有效地對信號進(jìn)行分離,提取包含故障特征最完整的信號分量,再通過對該分量的包絡(luò)分析提取出故障特征頻率,滾動軸承故障特征提取的具體操作流程如圖9所示。

圖9 滾動軸承特征提取流程圖

4.3 滾動軸承內(nèi)圈故障信號分析

選取滾動軸承內(nèi)圈故障數(shù)據(jù)進(jìn)行GVMD分解,通過對原始數(shù)據(jù)進(jìn)行小波變換得到預(yù)分解分量個(gè)數(shù),小波變換后的時(shí)頻圖如圖10所示,雖然小波變換得到的時(shí)頻譜分辨率低,但可以看出原始信號包含了4個(gè)單分量,故K=4。對數(shù)據(jù)進(jìn)行希爾伯特黃變換得到相位函數(shù),選取沖擊性最強(qiáng)分量的相位函數(shù)d(t)=47 540t+51對解析信號進(jìn)行廣義傅里葉變換,再對變換后的信號進(jìn)行VMD分解。原始故障信號的時(shí)域圖如圖11所示,GVMD分解后的分量如圖12與圖13所示。

圖10 小波變換時(shí)頻圖

圖11 滾動軸承內(nèi)圈故障數(shù)據(jù)時(shí)域圖

圖12 GVMD分解后各分量時(shí)域圖

圖13 GVMD分解后各分量頻譜圖

由相關(guān)系數(shù)法則可知,相關(guān)系數(shù)越大的分量與原始信號的故障特征越接近,分解后得到各IMF分量的相關(guān)系數(shù)如表2所示。

表2 IMF分量的相關(guān)系數(shù)

由表2可知,IMF1的相關(guān)系數(shù)最大,因此,選用IMF1作為下一步分析的目標(biāo)信號。

令p∈[1.4,1.6],以Δp=0.001為步長,利用離散采樣法FRFT算法可得到分?jǐn)?shù)階傅里葉變換結(jié)果如圖14所示,利用峰值搜索方法可得出峰值,即p=1.452,u0=9 416,故可設(shè)計(jì)濾波器為Mp(u)=[zeros(1,u0-1),ones(1,1),zeros(16384-u0)],該濾波器僅允許u0處的信號通過。

圖14 分?jǐn)?shù)階傅里葉變換后的結(jié)果圖

含噪聲信號的分?jǐn)?shù)域與幅值波形如圖15所示,經(jīng)過分?jǐn)?shù)階域?yàn)V波后信號的分?jǐn)?shù)域與幅值波形如圖16所示,由此可以明顯看出,濾波前的信號包含很多噪聲干擾,而濾波后的信號無噪聲干擾,這說明FRFT實(shí)現(xiàn)了對信號的降噪處理。

圖15 含噪聲的故障振動信號分?jǐn)?shù)域與幅值關(guān)系

圖16 濾波后信號的分?jǐn)?shù)域與幅值關(guān)系

最后,對該信號做分?jǐn)?shù)階傅里葉逆變換,得到如圖17所示的時(shí)域波形圖。

圖17 濾波后的信號時(shí)域圖

濾波后對故障的時(shí)域信號進(jìn)行1.5維包絡(luò)解調(diào)如圖18所示,故障頻率值為122.1 Hz。所用實(shí)驗(yàn)數(shù)據(jù)采集時(shí)間很短,持續(xù)了5.46 s,轉(zhuǎn)速從1 092 r/min變化至2 100 r/min,波動范圍很小,故以1 350 r/min為參考標(biāo)準(zhǔn)計(jì)算得出的理論故障頻率值為122.1 Hz,與實(shí)際值保持一致,成功地提取出了故障特征頻率。

圖18 1.5維包絡(luò)譜圖

5 結(jié)論

基于 GVMD 和FRFT的故障特征提取方法的準(zhǔn)確性和優(yōu)越性主要表現(xiàn)在:

1) 針對變工況下滾動軸承故障特征頻率為時(shí)變非平穩(wěn)信號提出了廣義變分模態(tài)分解的方法,該方法將時(shí)頻曲線呈非線性分布的變工況振動故障信號解調(diào)為時(shí)頻曲線近似線性分布的信號,有效抑制了模態(tài)混疊,從而提取出有效的故障特征分量。

2) 針對強(qiáng)大背景噪聲下存在故障特征提取困難,以及傳統(tǒng)的降噪方法由于經(jīng)驗(yàn)設(shè)置導(dǎo)致結(jié)果不準(zhǔn)確的問題,提出了分?jǐn)?shù)階傅里葉變換濾波的方法,該方法克服了信號交叉項(xiàng)的干擾,能在分?jǐn)?shù)階域上對噪聲信號進(jìn)行分離,可以實(shí)現(xiàn)精細(xì)的濾波。