基于參數(shù)自尋優(yōu)變分模態(tài)分解的信號(hào)降噪方法

何成兵, 車其祥, 徐振華, 于慶彬, 董玉亮, 程睿翔

(1.華北電力大學(xué) 能源動(dòng)力與機(jī)械工程學(xué)院,北京 102206;2.國(guó)網(wǎng)福建省電力公司電力科學(xué)研究院,福州 350007;3.國(guó)網(wǎng)山東省電力公司電力科學(xué)研究院,濟(jì)南 250000)

滾動(dòng)軸承在各類旋轉(zhuǎn)機(jī)械設(shè)備中應(yīng)用廣泛,同時(shí)也是設(shè)備中最易損壞部件之一。統(tǒng)計(jì)數(shù)據(jù)表明,在使用滾動(dòng)軸承的旋轉(zhuǎn)機(jī)械中,約30%的機(jī)械故障是由滾動(dòng)軸承引起的。滾動(dòng)軸承振動(dòng)信號(hào)具有沖擊性、非線性、非平穩(wěn)的特性,信號(hào)微弱、調(diào)制性強(qiáng),而且由于旋轉(zhuǎn)機(jī)械設(shè)備運(yùn)行工況復(fù)雜,背景噪聲污染嚴(yán)重,導(dǎo)致滾動(dòng)軸承故障信號(hào)常常淹沒(méi)于噪聲之中[1]。因此,對(duì)滾動(dòng)軸承故障信號(hào)進(jìn)行降噪處理,獲得有效故障特征,對(duì)于后續(xù)的故障信號(hào)分析與故障診斷具有重要的意義。

目前常用的降噪方法有小波變換(wavelet transform, WT)、經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)、集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)以及變分模態(tài)分解(variational mode decomposition, VMD)等。WT在信號(hào)降噪處理時(shí)的缺點(diǎn)是小波基函數(shù)和小波閾值函數(shù)的選擇不具有唯一性,受主觀性影響較大,很難保證降噪效果最優(yōu)。為克服WT的缺點(diǎn),曹玲玲等[2]結(jié)合快速譜峭度和帶通濾波以及Hilbert包絡(luò)分析,提出了一種新的改進(jìn)小波閾值函數(shù)進(jìn)行信號(hào)降噪;Malla等[3]利用復(fù)Morlet小波分析,得到滾動(dòng)軸承故障的相位圖和振幅圖,有效提取出了軸承故障特征。

EMD方法在信號(hào)分解過(guò)程中存在模態(tài)混疊、端點(diǎn)效應(yīng)及產(chǎn)生虛假分量等缺陷,使得信號(hào)中特征頻率與噪聲的分離程度不足[4]。EEMD對(duì)EMD進(jìn)行了優(yōu)化改進(jìn),李思琦等[5]利用EEMD對(duì)信號(hào)進(jìn)行分解,提出了一種基于EEMD與卷積神經(jīng)網(wǎng)絡(luò)相結(jié)合的滾動(dòng)軸承故障診斷方法;Yin等[6]利用改進(jìn)EEMD將原始信號(hào)中隱藏的固有噪聲分離,提出了一種基于改進(jìn)EEMD和自適應(yīng)閾值去噪的滾動(dòng)軸承弱故障特征提取新方法。EEMD可以在一定程度上抑制模態(tài)混疊現(xiàn)象,但也存在計(jì)算量大、頻譜剖分效果不理想的缺點(diǎn)。

VMD是一種新的非遞歸信號(hào)分解方法,它克服了EMD等算法存在的模態(tài)混疊和頻率效應(yīng)等缺點(diǎn),具有很好的噪聲魯棒性和降噪效果,得到了廣泛應(yīng)有。Wang等[7]采用VMD進(jìn)行碰摩信號(hào)的等效濾波特性分析,并與EWT、EMD和EEMD方法進(jìn)行了比較,證明VMD在提取瞬態(tài)沖擊方面更有效;朵慕社等[8-10]提出一種基于改進(jìn)VMD降噪和卷積神經(jīng)網(wǎng)絡(luò)的軸承故障智能診斷方法,其VMD分解層數(shù)通過(guò)排列熵閾值法確定,結(jié)合峭度準(zhǔn)則和互相關(guān)準(zhǔn)則,進(jìn)行IMF分量選取與信號(hào)重構(gòu),取得了較為理想的降噪效果。

在利用VMD進(jìn)行信號(hào)降噪處理時(shí),需提前設(shè)定模態(tài)數(shù)K和二次懲罰因子α等參數(shù),由于K和α的取值對(duì)降噪效果有很大的影響,它需要自適應(yīng)獲得最佳值,而不是人為主觀設(shè)定。目前主要采用以下三種途徑來(lái)獲得K和α的優(yōu)化參數(shù)[11]:一是根據(jù)先驗(yàn)知識(shí)或中心頻率觀察法選取,如Qiao等[12]對(duì)原始信號(hào)進(jìn)行VMD分解時(shí),根據(jù)每個(gè)IMF分量的中心頻率確定K值,這種方法適應(yīng)性較差,且不能保證VMD分解精度;二是通過(guò)評(píng)價(jià)指標(biāo)來(lái)選取,如Zhang等[13]通過(guò)能量和相關(guān)系數(shù)確定K值,這種方法選取的K值并不適用于所有的信號(hào);三是采用元啟發(fā)式算法進(jìn)行參數(shù)自適應(yīng)尋優(yōu),如人工魚群算法(artificial fish swarm, AFS)[14]、灰狼算法(grey wolf, GW)[15]、蚱蜢算法(grasshopper optimization, GO)[16]、鯨魚算法(whale optimization, WO)[17]、蝙蝠算法(bat optimization, BO)[18]和粒子群算法(particle swarm optimization, PSO)[19]等,這些優(yōu)化算法可搜尋到優(yōu)化參數(shù),但一般需要大量的迭代計(jì)算,同時(shí),由于各實(shí)際應(yīng)有領(lǐng)域的信號(hào)與噪聲特點(diǎn)不同,這些算法在應(yīng)用于具體領(lǐng)域時(shí),如何設(shè)定恰當(dāng)?shù)倪m應(yīng)度函數(shù)、如何提高算法搜索能力以及避免陷入局部最優(yōu)解,是急需解決的問(wèn)題。

針對(duì)上述問(wèn)題,本文提出了一種基于參數(shù)自尋優(yōu)變分模態(tài)分解的信號(hào)降噪方法。首先提出了一種改進(jìn)粒子群算法(improved particle swarm optimization, IPSO),以實(shí)現(xiàn)VMD最優(yōu)模態(tài)數(shù)K和二次懲罰因子α的自適應(yīng)尋優(yōu),該算法結(jié)合滾動(dòng)軸承故障信號(hào)特點(diǎn),建立了模態(tài)復(fù)合熵新指標(biāo),并以之作為適應(yīng)度函數(shù),同時(shí)設(shè)計(jì)了慣性權(quán)重值隨搜索進(jìn)程先大后小的計(jì)算公式、邊界粒子以及粒子群優(yōu)化處理原則,既提高了IPSO算法搜索能力,又避免了算法早熟收斂、陷入局部最優(yōu)解;然后基于最優(yōu)K和α,對(duì)原始信號(hào)進(jìn)行VMD分解,獲得K個(gè)IMF分量;利用相關(guān)系數(shù)篩選法,進(jìn)行IMF分量的有效模態(tài)和含噪模態(tài)識(shí)別,并利用小波閾值去噪方法對(duì)含噪模態(tài)進(jìn)行去噪處理;最后將有效模態(tài)與去噪后的含噪模態(tài)進(jìn)行重構(gòu),實(shí)現(xiàn)信號(hào)降噪。數(shù)值仿真和試驗(yàn)數(shù)據(jù)分析表明本文所提方法降噪效果明顯,有利于滾動(dòng)軸承故障特征的提取。

1 相關(guān)理論基礎(chǔ)

1.1 VMD算法

VMD將一個(gè)時(shí)域信號(hào)f(t)分解為K個(gè)具有調(diào)頻-調(diào)幅特性和稀疏特性的本征模態(tài)函數(shù)uk(t)

uk(t)=Ak(t)cos(φk(t))

(1)

各IMF中心頻率為ωk(t),為使每個(gè)模態(tài)的估計(jì)帶寬之和最小,建立如下約束變分模型

(2)

為求取式(2)最優(yōu)解,引入擴(kuò)展Lagrange函數(shù),將約束性變分問(wèn)題變換為非約束性變分問(wèn)題,其表達(dá)式為

(3)

式中,α為二次懲罰因子,用于保證信號(hào)的重構(gòu)精度,α值越大,各模態(tài)分量的頻率帶寬就越小。

利用交替方向乘子算法求取式(3)的鞍點(diǎn),即求取約束變分模型的最優(yōu)解,其中模態(tài)分量uk及中心頻率ωk分別為

(4)

(5)

VMD算法實(shí)現(xiàn)過(guò)程如下:

(6)

步驟4 重復(fù)步驟2、步驟3,直到滿足式(7)所示的迭代停止條件,結(jié)束循環(huán),得到K個(gè)變分模態(tài)分量。

(7)

式中,ε為求解精度。

1.2 PSO算法

PSO算法模擬鳥群的捕食行為,將每只鳥視為粒子,作為基本運(yùn)算單位,其優(yōu)化問(wèn)題的實(shí)質(zhì)是通過(guò)粒子的運(yùn)動(dòng)方向和速度確定最優(yōu)解在空間中的位置。其基本原理為:設(shè)搜索域?yàn)镈維,共有m個(gè)粒子構(gòu)成群體,第i個(gè)粒子位置表示為向量xi=(xi1,xi2,…,xiD),飛行速度表示為向量vi=(vi1,vi2,…,viD),xi位置變化就是解的軌跡,用適應(yīng)度函數(shù)Q(xi)來(lái)權(quán)衡粒子位置的好壞。第i個(gè)粒子所經(jīng)歷的個(gè)體最優(yōu)位置記為pi=(pi1,pi2,…,piD) ,搜索域中整個(gè)粒子群全局最優(yōu)位置記為pg=(pg1,pg2,…,pgD)。各粒子的速度和位置通過(guò)下式更新

vid(k+1)=ωvid(k)+c1r1(pid(k)-xid(k))+

c2r2(pgd(k)-xid(k))

(8)

xid(k+1)=xid(k)+vid(k+1)

(9)

式中:k為粒子群迭代次數(shù);i=[1,m]為粒子數(shù);d=[1,D]為空間維數(shù);r1和r2為相互獨(dú)立的0～1之間的隨機(jī)數(shù);c1和c2分別為粒子向個(gè)體和全局最優(yōu)位置方向的學(xué)習(xí)因子;vid(k)和xid(k)分別為第i個(gè)粒子在d維空間內(nèi)k次迭代的速度和位置;pid(k)和pgd(k)分別為第i個(gè)粒子在d維空間內(nèi)k次迭代的個(gè)體和全局最優(yōu)位置;ω為慣性權(quán)重。

PSO算法實(shí)現(xiàn)過(guò)程如下:

步驟1 算法初始化,設(shè)定粒子群的維度D、群體規(guī)模m、最大迭代次數(shù)Np、粒子位置范圍、速度范圍。

步驟2 計(jì)算每個(gè)粒子的適應(yīng)度函數(shù)值,確定粒子個(gè)體最優(yōu)位置pi=(pi1,pi2,…,piD)和粒子群全局最優(yōu)位置pg=(pg1,pg2,…,pgD)。

步驟3 根據(jù)式(8)、式(9)更新調(diào)整粒子的速度和位置。

步驟4 迭代終止條件一般選為最大迭代次數(shù)Np或粒子群最優(yōu)位置小于設(shè)定的最小閾值,判斷是否滿足迭代終止條件,不滿足則轉(zhuǎn)到步驟2,滿足則結(jié)束。

1.3 小波閾值去噪算法

小波閾值去噪是建立在小波變換基礎(chǔ)上的算法,其基本原理為:對(duì)含噪信號(hào)進(jìn)行小波變換,根據(jù)信號(hào)具有真實(shí)信號(hào)通常為低頻、噪聲信號(hào)多為高頻的特點(diǎn),對(duì)各頻帶小波系數(shù)與小波閾值進(jìn)行比較,甄別出有效信號(hào)與噪聲信號(hào),并通過(guò)小波重構(gòu)得到去噪信號(hào)。

小波閾值去噪算法實(shí)現(xiàn)過(guò)程如下:

步驟1 設(shè)定小波分解層數(shù),選定小波基函數(shù),構(gòu)造小波閾值函數(shù)。

步驟2 使用選定的小波基函數(shù)對(duì)含噪信號(hào)進(jìn)行Mallat算法處理,得到分解層數(shù)對(duì)應(yīng)的小波高頻系數(shù)和低頻系數(shù);

步驟3 將小波系數(shù)與小波閾值進(jìn)行比較,如果小波系數(shù)大于小波閾值,對(duì)應(yīng)數(shù)據(jù)為有效信號(hào),則該小波系數(shù)保留;反之,為噪聲信號(hào),該小波系數(shù)舍棄。

步驟4 對(duì)保留的小波系數(shù)重構(gòu),得到去噪信號(hào)。

2 改進(jìn)粒子群算法(IPSO)

本文提出了一種改進(jìn)粒子群算法(improved particle swarm optimization,IPSO)進(jìn)行VMD參數(shù)自適應(yīng)尋優(yōu),以確定VMD最優(yōu)的K和α。該算法的特點(diǎn)在于:從滾動(dòng)軸承故障信號(hào)特點(diǎn)出發(fā),提出了模態(tài)復(fù)合熵新指標(biāo),并以之作為適應(yīng)度函數(shù);設(shè)計(jì)了慣性權(quán)重值隨搜索進(jìn)程先大后小的計(jì)算公式、邊界粒子以及粒子群優(yōu)化處理原則,提高了IPSO算法搜索能力,也避免了算法早熟收斂、陷入局部最優(yōu)解。

IPSO算法流程如圖1所示,具體實(shí)現(xiàn)過(guò)程如下:

圖1 IPSO流程圖Fig.1 The flow diagram of IPSO

步驟1 算法初始化。設(shè)定粒子群維度數(shù)D為 2,對(duì)應(yīng)VMD的模態(tài)數(shù)K和二次懲罰因子α;在2維可行域中隨機(jī)選取m個(gè)粒子,對(duì)粒子位置、飛行速度隨機(jī)初始化;第i個(gè)粒子在2維空間的位置表示為向量xi=(xi1,xi2),飛行速度表示為向量vi=(vi1,vi2)。

步驟2 針對(duì)某一粒子的2維位置組合,對(duì)原始信號(hào)做VMD運(yùn)算,得到K個(gè)模態(tài)分量{uk},計(jì)算每一粒子的適應(yīng)度函數(shù)值Q(uki(1)),對(duì)各粒子進(jìn)行評(píng)價(jià),選取適應(yīng)度函數(shù)值最小者為最優(yōu),存儲(chǔ)對(duì)應(yīng)粒子的個(gè)體最優(yōu)位置pi=(pi1,pi2) 以及整個(gè)粒子群全局最優(yōu)位置pg=(pg1,pg2)。

通過(guò)對(duì){uk}分析發(fā)現(xiàn):如果{uk}中包含噪聲越多,則信號(hào)復(fù)雜性越高,對(duì)應(yīng)的模態(tài)復(fù)合熵越大,反之,則模態(tài)復(fù)合熵越小。因此,以{uk}的模態(tài)復(fù)合熵Ec作為適應(yīng)度函數(shù)Q(uk),且以模態(tài)復(fù)合熵最小作為適應(yīng)值,適應(yīng)度函數(shù)如式(10)所示

Q(uk)=Ec=β1·Eenergy+β2·Eenvelope

(10)

式中:Eenergy為模態(tài)能量熵;Eenvelope為模態(tài)包絡(luò)熵;β1,β2為二者的權(quán)重系數(shù),β1+β2=1,對(duì)于如滾動(dòng)軸承故障這樣具有明顯沖擊特征的信號(hào),取β1<β2,否則取β1≥β2。

模態(tài)能量熵Eenergy的計(jì)算式如下

(11)

模態(tài)包絡(luò)熵Eenvelope的計(jì)算式如下

(12)

式中,H為信號(hào)uk的Hilbert變換。

步驟3 根據(jù)迭代公式(8)更新各粒子速度,根據(jù)式(9)更新各粒子位置。

式(8)中慣性權(quán)重ω對(duì)算法搜索結(jié)果影響較大,其取值大則種群粒子搜索能力強(qiáng),可探索較大的區(qū)域,取值小則可精細(xì)搜索目標(biāo),局部搜索能力強(qiáng)。為提高算法搜索能力,同時(shí)避免陷入局部最優(yōu)解,采用ω權(quán)重值隨搜索進(jìn)程先大后小的原則動(dòng)態(tài)設(shè)定

(13)

式中:ωmax,ωmin為慣性權(quán)重的最大值和最小值,取ωmax=0.9,ωmin=0.4;t為粒子迭代次數(shù);Np為最大迭代次數(shù);Dec為衰減因子,其值為小于1的正數(shù)。

粒子位置更新后,有可能出現(xiàn)粒子的某一維度超出設(shè)定范圍的情況,此時(shí)粒子容易在邊界附近聚集,導(dǎo)致算法早熟收斂。因此,對(duì)超出邊界范圍的粒子進(jìn)行邊界優(yōu)化處理:將該維度的位置值設(shè)置為邊界值;該維度的速度乘以0～1之間的隨機(jī)數(shù),并將其方向設(shè)為反方向,以此緩解因粒子在邊界聚集而引起的算法早熟收斂。計(jì)算公式如下

(14)

式中:xid,vid是第i個(gè)粒子第d維的位置和速度;xmaxd,xmind是粒子第d維的上限和下限值。

步驟4 針對(duì)位置更新后的粒子,對(duì)原始信號(hào)做VMD運(yùn)算,計(jì)算每一粒子新的適應(yīng)度函數(shù)值Q(uki(k+1)),并更新粒子個(gè)體最優(yōu)位置pi(k+1)及粒子群全局最優(yōu)位置pg(k+1)。粒子i的個(gè)體最優(yōu)位置更新公式如下

pi(k+1)=

(15)

粒子群全局最優(yōu)位置更新公式如下

(16)

步驟5 為提升算法性能,提高算法收斂速度,對(duì)粒子群進(jìn)行優(yōu)化處理:按照適應(yīng)度函數(shù)值對(duì)所有粒子進(jìn)行排序,并將粒子中適應(yīng)度較差的一半用另一半粒子的速度和位置代替,但不改變?cè)辛Ｗ拥膫€(gè)體最優(yōu)位置pi(k+1)和粒子群全局最優(yōu)位置pg(k+1)。

步驟6 判斷當(dāng)前迭代次數(shù)是否達(dá)到最大迭代次數(shù)Np,如果未達(dá)到,返回步驟3進(jìn)行下一次迭代;如果已達(dá)到,則迭代終止,輸出最優(yōu)粒子2維位置,即獲得最優(yōu)模態(tài)數(shù)K和二次懲罰因子α。

3 本文所提降噪方法流程

本文提出的基于參數(shù)自尋優(yōu)變分模態(tài)分解與小波閾值處理的信號(hào)降噪方法流程如圖2所示,具體實(shí)現(xiàn)步驟為:

圖2 本文降噪方法流程圖Fig.2 The flow diagram of signal denoising method in this paper

步驟1:采集獲取原始信號(hào)x(t);

步驟2:以模態(tài)復(fù)合熵作為適應(yīng)度函數(shù),采用IPSO算法進(jìn)行VMD參數(shù)自適應(yīng)尋優(yōu),確定VMD最優(yōu)模態(tài)數(shù)K和二次懲罰因子α;

步驟3:基于最優(yōu)模態(tài)數(shù)K和二次懲罰因子α,對(duì)含有噪聲的原始信號(hào)進(jìn)行VMD處理,得到K個(gè)本征模態(tài)分量{uk};

步驟4:利用相關(guān)系數(shù)篩選法,進(jìn)行模態(tài)分量{uk}的有效模態(tài)和含噪模態(tài)分類,利用小波閾值去噪方法對(duì)含噪模態(tài)進(jìn)行去噪處理。具體流程如下:

(1) 利用式(17)計(jì)算各模態(tài)分量{uk}與原始信號(hào)x(t)之間的相關(guān)系數(shù)r,相關(guān)系數(shù)越大,表示兩信號(hào)之間的關(guān)聯(lián)性越強(qiáng)

(17)

利用下式計(jì)算相關(guān)系數(shù)閾值r0

(18)

式中,rmax為各模態(tài)分量相關(guān)系數(shù)最大值。

(2) 各模態(tài)分量相關(guān)系數(shù)ri分別與閾值r0進(jìn)行比較,如果ri≥r0,對(duì)應(yīng)的模態(tài)分量歸于有效模態(tài)組{uv1,uv2,…,uvp};如果ri

(3) 當(dāng)含噪模態(tài)組{un1,un2,…,unq}只有1個(gè)模態(tài)分量時(shí),直接去掉該模態(tài)分量;有兩個(gè)及以上模態(tài)分量時(shí),因相關(guān)系數(shù)最小的模態(tài)分量unq與原始信號(hào)關(guān)聯(lián)性最弱,通常為噪聲干擾信號(hào),將其先去除,含噪模態(tài)組變?yōu)閧un1,un2,…,un(q-1)},針對(duì)該含噪模態(tài)組按1.3節(jié)方法進(jìn)行小波閾值去噪處理,其中,設(shè)定小波基函數(shù)為db8、小波分解層數(shù)為5層、小波閾值公式如下

(19)

式中:N為信號(hào)長(zhǎng)度;σ為信號(hào)噪聲方差,其計(jì)算式為

σ=median(|Wj,k|)/0.674 5

(20)

式中,Wj,k為信號(hào)分解后得到的小波系數(shù)。

(21)

4 仿真信號(hào)降噪分析

4.1 降噪評(píng)價(jià)指標(biāo)

本文采用的信號(hào)降噪結(jié)果評(píng)價(jià)指標(biāo)包括:信噪比SNR、能量百分比ESN、均方根誤差RMSE以及波形相似系數(shù)NCC。

信噪比SNR計(jì)算式如下

(22)

式中:x(t)為含噪信號(hào);x′(t)為降噪后信號(hào);N為數(shù)據(jù)點(diǎn)數(shù)。SNR值越大表示降噪效果越好。

能量百分比ESN計(jì)算式如下

(23)

式中:E′為降噪后信號(hào)能量;E為含噪信號(hào)能量,ESN值越大表示降噪效果越好。

均方根誤差RMSE計(jì)算式如下

(24)

RMSE值越小表示降噪效果越好。

波形相似系數(shù)NCC計(jì)算式如下

(25)

NCC值越大表示降噪效果越好。

4.2 仿真信號(hào)構(gòu)建

基于滾動(dòng)軸承外圈故障機(jī)理與信號(hào)特征,構(gòu)建了式(26)所示仿真信號(hào)[20]

(26)

式中:x(t)為軸承外圈故障含噪信號(hào);y(t)為不含噪信號(hào);x1(t)為共振干擾脈沖信號(hào),設(shè)fr為軸轉(zhuǎn)頻,取fr=18 Hz,則相應(yīng)的干擾脈沖間隔周期T1=1/fr=0.056 s;x2(t)為外圈故障脈沖信號(hào),設(shè)f0為外圈故障頻率,取f0=107 Hz,則相應(yīng)的故障脈沖間隔周期T2=1/f0=0.009 3 s;s(t)為諧波信號(hào),其幅值取B0=0.4 m/s2;n(t)為高斯白噪聲信號(hào);Ai為調(diào)頻信號(hào),其頻率為fr,幅值取B1=0.6 m/s2;C為衰減系數(shù),取C=800;fn為系統(tǒng)固有頻率,取fn=2 300 Hz。仿真信號(hào)采樣頻率設(shè)為12 000 Hz,采樣點(diǎn)數(shù)為16 384個(gè)。

圖3為不含噪仿真信號(hào)y(t)的時(shí)域波形、幅值譜和包絡(luò)譜圖,圖4為含噪仿真信號(hào)x(t)的時(shí)域波形、幅值譜和包絡(luò)譜圖。由圖3可知,不含噪聲仿真信號(hào)外圈故障特征明顯,尤其是包絡(luò)圖上,軸轉(zhuǎn)頻fr(18 Hz)、2倍頻2fr(36 Hz)、3倍頻3fr(54 Hz),外圈故障頻率f0(107 Hz)、2倍頻2f0(214 Hz)、3倍頻3f0(321 Hz)、4倍頻4f0(428 Hz),以及以故障特征頻率f0為中心、以轉(zhuǎn)頻fr為邊帶的各種調(diào)制頻率成分f0+fr(125 Hz)、f0-fr(89 Hz)、2f0+fr(232 Hz)、2f0-fr(196 Hz)等清晰可見,易于故障判定。由圖4可知,由于疊加了較強(qiáng)的高斯噪聲,外圈故障特征幾乎被淹沒(méi),體現(xiàn)在圖4(c)的包絡(luò)譜圖上,幾乎無(wú)法提取出外圈故障特征頻率。

(a) 時(shí)域波形曲線

(b) FFT頻譜圖

(c) 包絡(luò)譜圖圖3 不含噪仿真信號(hào)y(t)Fig.3 Simulated signal y(t) without noise

(a) 時(shí)域波形曲線

(b) FFT頻譜圖

(c) 包絡(luò)譜圖圖4 含噪仿真信號(hào)x(t)Fig.4 Simulated signal x(t) with noise

4.3 仿真信號(hào)分析

針對(duì)含噪仿真信號(hào)x(t)分別進(jìn)行EMD和本文方法的降噪處理。圖5(a)為EMD分解得到的每個(gè)IMF分量時(shí)域波形,圖5(b)為每個(gè)IMF分量對(duì)應(yīng)的FFT頻譜圖。由圖5可見,EMD分解結(jié)果中存在較為嚴(yán)重的模態(tài)混疊現(xiàn)象,各模態(tài)分量的頻段分離效果較差,以IMF2為例,系統(tǒng)固有頻率fn及其邊頻帶和其他中低頻率成分混疊在一起,不能有效提取有用信息;IMF2～I(xiàn)MF5,相鄰的IMF之間均產(chǎn)生了模態(tài)混疊,每個(gè)IMF均包含了低頻成分。與此同時(shí),EMD分解過(guò)程中產(chǎn)生了虛假分量,如IMF3和IMF4中對(duì)應(yīng)的頻率成分并不存在。

(a) IMF時(shí)域波形曲線

(b) IMF頻譜圖圖5 含噪仿真信號(hào)EMD分解Fig.5 EMD decomposition of simulation signal with noisy

利用本文提出的方法,首先進(jìn)行基于IPSO算法的K和α參數(shù)尋優(yōu),K和α的尋優(yōu)范圍設(shè)置為[3,10]和[200,6 000],獲得的最優(yōu)參數(shù)值K=5,α=2 182;基于該最優(yōu)參數(shù)進(jìn)行VMD分解,得到的各IMF分量時(shí)域波形如圖6(a)所示,各IMF分量對(duì)應(yīng)的FFT頻譜如圖6(b)所示。由圖6可見,本文方法可有效去除模態(tài)混疊現(xiàn)象,含噪信號(hào)經(jīng)分解后實(shí)現(xiàn)了信號(hào)頻域及各個(gè)IMF分量的自適應(yīng)剖分,每個(gè)IMF都圍繞在某一中心頻率處,軸轉(zhuǎn)頻成分被分解到IMF1,系統(tǒng)固有頻率fn被分解到IMF3,而系統(tǒng)固有頻率fn邊頻帶則被分解到IMF2和IMF4上,噪聲成分被分解到IMF5中;同時(shí)該方法也避免了虛假模態(tài)的產(chǎn)生。由此可見,相較于EMD算法,本文方法具有更好的分解效果與噪聲魯棒性。

(a) IMF時(shí)域波形曲線

(b) IMF頻譜圖圖6 含噪仿真信號(hào)VMD分解Fig.6 VMD decomposition of simulation signal with noisy

分別針對(duì)EMD和本文方法分解得到的各IMF分量與含噪仿真信號(hào)x(t)進(jìn)行互相關(guān)計(jì)算,得到的互相關(guān)系數(shù)r如表1所示。

表1 各IMF分量與x(t)的互相關(guān)系數(shù)riTab.1 Cross correlation coefficient ri between each IMF component and x(t)

根據(jù)式(18)確定相關(guān)系數(shù)閾值r0,然后對(duì)表1中滿足系數(shù)閾值條件的IMF 進(jìn)行信號(hào)重構(gòu)。經(jīng)計(jì)算,EMD分解的各模態(tài)函數(shù)分量r0=0.16,本文方法分解的各模態(tài)函數(shù)分量r0=0.17。因此,結(jié)合表1,EMD分析結(jié)果中選取IMF1、IMF6,其他各模態(tài)函數(shù)分量舍棄,以此重構(gòu)信號(hào);本文方法分析結(jié)果中將IMF1、IMF2、IMF3、IMF4歸為有效模態(tài)組,IMF5歸為含噪模態(tài)組,由于含噪模態(tài)組只有1個(gè)分量,不用進(jìn)行小波閾值處理,直接舍棄,以有效模態(tài)組4個(gè)分量重構(gòu)信號(hào)。最終得到2組降噪后信號(hào),分別示于圖7和圖8。圖7(a)是EMD降噪信號(hào)時(shí)域波形,圖8(a)是本文方法降噪信號(hào)時(shí)域波形,與圖3(a)不含噪仿真信號(hào)y(t)的時(shí)域波形相比,可定性看出EMD降噪后仍存在一定程度上的噪聲干擾,本文方法降噪信號(hào)波形比EMD降噪更接近不含噪仿真信號(hào)。圖7(b)和圖8(b)是兩組降噪后信號(hào)的包絡(luò)譜圖,與圖3(c)不含噪仿真信號(hào)y(t)的包絡(luò)圖均極為接近,故障頻率f0及其各種調(diào)制頻率等特征頻率清晰可見。

(a) 時(shí)域波形曲線

(b) 包絡(luò)譜圖圖7 EMD方法降噪結(jié)果Fig.7 Noise reduction results of EMD method

(a) 時(shí)域波形曲線

(b) 包絡(luò)譜圖圖8 本文方法降噪結(jié)果Fig.8 Noise reduction results of this paper method

應(yīng)用式(22)～式(25),針對(duì)降噪后信號(hào)進(jìn)行降噪指標(biāo)計(jì)算,結(jié)果如表2所示(表中還列出了EEMD評(píng)價(jià)結(jié)果)。可以觀察到本文方法降噪結(jié)果信噪比更高、能量百分比更大、均方根誤差更小、平滑度更接近1,各項(xiàng)評(píng)價(jià)指標(biāo)全面好于EMD和EEMD方法,說(shuō)明本文方法降噪效果明顯優(yōu)于EMD和EEMD降噪方法。

表2 降噪評(píng)價(jià)指標(biāo)結(jié)果Tab.2 Evaluation index results of noise reduction

5 試驗(yàn)信號(hào)降噪分析

采用美國(guó)凱斯西儲(chǔ)大學(xué)的滾動(dòng)軸承故障試驗(yàn)數(shù)據(jù)驗(yàn)證本文降噪方法的有效性,試驗(yàn)系統(tǒng)如圖9所示。

圖9 滾動(dòng)軸承試驗(yàn)臺(tái)Fig.9 Rolling bearing test bench

實(shí)驗(yàn)臺(tái)驅(qū)動(dòng)端軸承型號(hào)為6205-2RS JEM SKF深溝球軸承,技術(shù)參數(shù)和規(guī)格信息如表3所示。

表3 滾動(dòng)軸承技術(shù)參數(shù)和規(guī)格信息Tab.3 Technical parameters and specifications of rolling bearing

采用電火加工技術(shù)分別在軸承內(nèi)、外圈布置單點(diǎn)故障,選擇驅(qū)動(dòng)端軸承外圈故障數(shù)據(jù),故障直徑為0.533 4 mm,深度為0.279 4 mm。該外圈故障對(duì)應(yīng)電機(jī)轉(zhuǎn)速為1 721 r/min(即轉(zhuǎn)頻fr為28.68 Hz),采樣頻率為12 000 Hz,采樣點(diǎn)數(shù)為16 384。結(jié)合表3中軸承各參數(shù),通過(guò)下式計(jì)算得到軸承外圈故障頻率f0為102.81 Hz。

(27)

試驗(yàn)采集到的原始故障信號(hào)相關(guān)信息如圖10所示,其中圖10(a)是其時(shí)域波形曲線,圖10(b)是其FFT分析幅值譜圖,圖10(c)是其包絡(luò)譜圖。由圖可知,故障時(shí)域波形沖擊特征較為明顯,但也存在大量噪聲;幅值譜圖上譜線主要分布在4個(gè)頻段區(qū)域:高頻區(qū)(中心頻率約3 400 Hz)、次高頻區(qū)(中心頻率約2 800 Hz)、中頻區(qū)(中心頻率約1 300 Hz)和低頻區(qū)(中心頻率約600 Hz),信號(hào)能量主要集中在高頻和次高頻區(qū),低頻區(qū)的能量較小,而軸承的外圈故障特征頻率幾乎無(wú)法提取;包絡(luò)圖上能觀察到軸轉(zhuǎn)頻和故障特征頻率,但由于噪聲干擾,圖上干擾譜線較多,故障特征頻率的調(diào)制頻帶等故障相關(guān)頻率成分均很難觀察到。

(a) 時(shí)域波形曲線

(b) FFT頻譜圖

(c) 包絡(luò)譜圖圖10 試驗(yàn)原始信號(hào)圖形Fig.10 Experimental original signal graph

利用本文方法,通過(guò)IPSO算法獲得VMD最優(yōu)參數(shù)值K=7,α=2 684,進(jìn)一步對(duì)原始故障信號(hào)進(jìn)行VMD分解,得到的各IMF分量時(shí)域波形如圖11(a)所示,對(duì)應(yīng)的FFT幅值譜如圖11(b)所示。將各IMF分量與原始信號(hào)進(jìn)行互相關(guān)計(jì)算,得到的互相關(guān)系數(shù)如表4所示。

表4 各IMF分量與原始信號(hào)的互相關(guān)系數(shù)Tab.4 Cross correlation coefficient between each IMF component and original signal

(a) IMF時(shí)域波形曲線

(b) IMF頻譜圖圖11 試驗(yàn)原始信號(hào)VMD分解Fig.11 VMD decomposition of test original signal

結(jié)合圖11和表4可以看出含噪信號(hào)經(jīng)VMD分解后實(shí)現(xiàn)了信號(hào)頻域及各個(gè)IMF分量的自適應(yīng)剖分,每個(gè)IMF都圍繞在4個(gè)頻段區(qū)域的某一中心頻率處,低頻段區(qū)被分解到IMF1和IMF2,二者有一定的頻率重疊,其互相關(guān)系數(shù)值也較小,中頻段區(qū)被分解到IMF3,次高頻段區(qū)被分解到IMF4,高頻段區(qū)由于能量占比最大,它被分別分解到IMF5、IMF6和IMF7,其互相關(guān)系數(shù)值也最大。

由式(18)可確定相關(guān)系數(shù)閾值r0=0.18,結(jié)合表4,將IMF4、IMF5、IMF6、IMF7歸為有效模態(tài)組,IMF1、IMF2、IMF3歸為含噪模態(tài)組,其中IMF2相關(guān)系數(shù)最小,直接舍棄,對(duì)含噪模態(tài)組的IMF1和IMF3進(jìn)行小波閾值處理,并與有效模態(tài)組的4個(gè)IMF分量一起重構(gòu)信號(hào),得到降噪信號(hào),其時(shí)域波形、FFT幅值譜、包絡(luò)譜示于圖12。與試驗(yàn)原始信號(hào)圖10相比,降噪后時(shí)域波形沖擊特征更明顯,幅值譜能量分布基本不變,噪聲干擾更少,譜線更清晰;降噪效果在包絡(luò)譜圖12(c)上體現(xiàn)最明顯,其譜圖較為干凈,干擾成分較少,可清晰發(fā)現(xiàn)軸轉(zhuǎn)頻fr(約29 Hz)、2倍轉(zhuǎn)頻2fr(58 Hz)成分,也可發(fā)現(xiàn)明顯的約103 Hz故障區(qū)域,對(duì)應(yīng)軸承外圈故障頻率f0,同時(shí)存在故障頻率的2倍頻2f0(206 Hz)、3倍頻3f0(309 Hz)成分,以及以故障特征頻率f0為中心以轉(zhuǎn)頻fr為邊帶的調(diào)制頻率成分f0+fr(132 Hz)和f0-fr(74 Hz),這與外圈故障特征頻率被轉(zhuǎn)頻所調(diào)制的特性相符。說(shuō)明本文方法有很好的降噪效果,據(jù)此提取的故障特征可用于軸承外圈故障的判定。

(a) 時(shí)域波形曲線

(b) FFT頻譜圖

(c) 包絡(luò)譜圖圖12 降噪后信號(hào)圖形Fig.12 Signal graph after noise reduction

6 結(jié) 論

本文針對(duì)滾動(dòng)軸承強(qiáng)背景噪聲故障信號(hào),提出了一種基于參數(shù)自尋優(yōu)變分模態(tài)分解的信號(hào)降噪方法。為實(shí)現(xiàn)VMD模態(tài)數(shù)K和二次懲罰因子α的自適應(yīng)尋優(yōu),提出了一種改進(jìn)粒子群算法(IPSO),通過(guò)構(gòu)建適合滾動(dòng)軸承故障信號(hào)特點(diǎn)的適應(yīng)度函數(shù)、慣性權(quán)重公式、邊界粒子以及粒子群優(yōu)化處理方法,提高了VMD參數(shù)自尋優(yōu)搜索能力;基于最優(yōu)K和α,對(duì)含噪信號(hào)進(jìn)行VMD分解,利用相關(guān)系數(shù)篩選法,對(duì)K個(gè)IMF分量進(jìn)行有效模態(tài)和含噪模態(tài)識(shí)別,對(duì)后者進(jìn)行小波閾值去噪,并與前者一起進(jìn)行信號(hào)重構(gòu),實(shí)現(xiàn)信號(hào)降噪。從數(shù)值仿真和試驗(yàn)數(shù)據(jù)角度,將本文所提方法與EMD和EEMD降噪方法進(jìn)行對(duì)比分析,從信噪比、能量百分比、均方根誤差以及波形相似系數(shù)4個(gè)評(píng)價(jià)指標(biāo),定量說(shuō)明本文方法在滾動(dòng)軸承故障信號(hào)降噪中具有更好的降噪效果,是一種有效的降噪方法。