李 波, 胡哿郗, 石劍鈞, 劉恒暢, 洪 濤

(1. 電子科技大學(xué)航空航天學(xué)院, 四川 成都 611731; 2.飛行器集群智能感知與協(xié)同控制四川省重點(diǎn)實(shí)驗(yàn)室, 四川 成都 611731)

0 引 言

滾動(dòng)軸承因其摩擦阻力小、功率消耗小和機(jī)械效率高等特點(diǎn),被廣泛應(yīng)用于航空發(fā)動(dòng)機(jī)[1-2]、鐵路運(yùn)輸設(shè)備[3-4]以及大型工業(yè)機(jī)械設(shè)備[5-6],承擔(dān)動(dòng)力傳遞和動(dòng)力轉(zhuǎn)換的任務(wù),是機(jī)械發(fā)生故障的主要故障源之一[7-8]。由于滾動(dòng)軸承的健康狀態(tài)直接關(guān)系設(shè)備的安全穩(wěn)定運(yùn)行,因此及時(shí)發(fā)現(xiàn)關(guān)鍵滾動(dòng)軸承中存在的故障,對(duì)確保機(jī)械安全穩(wěn)定運(yùn)行具有重大意義[9]。通常情況下,機(jī)械設(shè)備運(yùn)行狀態(tài)的變化可通過其動(dòng)態(tài)信號(hào)(例如振動(dòng)信號(hào))進(jìn)行有效反映,利用信號(hào)處理技術(shù)分離出其中隱藏的故障特征信息。然而,其動(dòng)態(tài)信號(hào)中蘊(yùn)含了復(fù)雜機(jī)械結(jié)構(gòu)中各個(gè)部件產(chǎn)生的信號(hào)成分,而且不可避免地存在背景噪聲,導(dǎo)致其中的故障特征信息被干擾成分淹沒。同時(shí),滾動(dòng)軸承復(fù)雜的工況使得動(dòng)態(tài)信號(hào)表現(xiàn)出調(diào)頻、調(diào)幅等復(fù)雜的時(shí)變特性,提高了提取其故障特征信息的難度[10-11]。由于上述軸承工作環(huán)境的特點(diǎn),如何在復(fù)雜噪聲環(huán)境下有效提取滾動(dòng)軸承故障特征變得非常關(guān)鍵。

迄今為止,學(xué)者們提出了多種信號(hào)處理技術(shù)提取滾動(dòng)軸承振動(dòng)信號(hào)中的故障特征信息,而最有效的則是基于時(shí)頻域的信號(hào)處理技術(shù)。研究初期,小波變換被應(yīng)用于滾動(dòng)軸承的故障特征提取中,實(shí)現(xiàn)了對(duì)信號(hào)時(shí)頻變化特性的分析,但其變換效果極其依賴經(jīng)驗(yàn)選取的小波基函數(shù)[12]。為此,Huang等[13]首次提出了一種自適應(yīng)時(shí)頻處理的經(jīng)驗(yàn)?zāi)B(tài)分解方法,合理解決了小波變換存在的“基選擇”經(jīng)驗(yàn)依賴問題,但其同時(shí)也存在端點(diǎn)效應(yīng)和模態(tài)混疊的缺點(diǎn)。在此基礎(chǔ)上,Wu等[14]提出了集成經(jīng)驗(yàn)?zāi)B(tài)分解,有效地解決了經(jīng)驗(yàn)?zāi)B(tài)分解存在的模態(tài)混疊現(xiàn)象,但該方法存在一個(gè)白噪聲選擇問題,易出現(xiàn)分解效果不佳的情況。進(jìn)而,Dragomiretskiy等[15]提出了變分模態(tài)分解方法,相比于集成經(jīng)驗(yàn)?zāi)B(tài)方法具有更好的信號(hào)局部性能和更優(yōu)的噪聲抑制能力,但也存在分解過程中參數(shù)難以設(shè)定的缺點(diǎn)。

為了解決變分模態(tài)分解中參數(shù)難以確定的問題,學(xué)者們陸續(xù)開展了相應(yīng)研究。王奉濤等[16]提出通過分解信號(hào)的能量差值來確定分解層數(shù)。畢鳳榮等[17]提出通過計(jì)算分解后各個(gè)本征模態(tài)函數(shù)(intrinsic mode function,IMF)的中心頻率差值來確定分解層數(shù)。李華等[18]提出通過最大峭度原則優(yōu)化分解層數(shù)。然而,上述研究只優(yōu)化了變分模態(tài)分解(variational mode decomposition, VMD)分解層數(shù)和懲罰因子中的一個(gè),沒有全面考慮到這兩個(gè)參數(shù)之間的相互作用對(duì)分解效果的影響。為此,程軍圣等[19]以正交低峰值為優(yōu)化目標(biāo)函數(shù),利用螢火蟲算法對(duì)分解層數(shù)與懲罰因子進(jìn)行聯(lián)合優(yōu)化。唐貴基等[20]提出以信號(hào)包絡(luò)熵作為粒子群算法的優(yōu)化目標(biāo)函數(shù)對(duì)參數(shù)進(jìn)行優(yōu)化。Xu等[21]引入Teager能量算子并結(jié)合最大峭度系數(shù)和平均峭度系數(shù)優(yōu)化分解層數(shù)和懲罰因子。雖然上述研究實(shí)現(xiàn)了VMD分解層數(shù)與懲罰因子的自適應(yīng)設(shè)定,解決了需手動(dòng)設(shè)置參數(shù)的問題,但沒有考慮IMF各個(gè)分量對(duì)應(yīng)的懲罰因子αk對(duì)VMD結(jié)果的影響。因此,鄭圓等[22]利用智能算法確定了分解層數(shù)和多懲罰因子αk,然而其在建立優(yōu)化目標(biāo)函數(shù)時(shí),只考慮了分解后信號(hào)與原始信號(hào)的時(shí)域差異,無法有效地反映分解后信號(hào)包含原始信號(hào)的信息完整度。由此可見,當(dāng)前利用VMD提取滾動(dòng)軸承的故障特征時(shí),其預(yù)先設(shè)置懲罰因子具有不確定性的問題并未得到合理解決,仍存在較大提升空間。

因此,本文提出一種基于多懲罰因子優(yōu)化VMD的滾動(dòng)軸承故障特征提取方法。首先,構(gòu)建了融合皮爾遜相關(guān)系數(shù)、時(shí)頻譜相似度以及最大峭度值的優(yōu)化目標(biāo),然后,通過最大峭度準(zhǔn)則優(yōu)化其分解層數(shù),再利用灰狼優(yōu)化(grey wolf optimization,GWO)算法獲取最佳分解層數(shù)以及最佳分解模態(tài)函數(shù),保證了分解后模態(tài)函數(shù)最大限度地包含原始信號(hào)的故障特征信息,提升了滾動(dòng)軸承振動(dòng)信號(hào)經(jīng)提取后故障特征的顯著性。最后,利用公開數(shù)據(jù)集分析驗(yàn)證,通過對(duì)比多種不同方法,得出本文所提方法的故障特征提取效果優(yōu)于其他同類方法,證明了該方法的有效性。

1 算法基本原理

1.1 VMD算法分解原理

VMD非遞歸式地將具有多個(gè)頻率成分的實(shí)信號(hào)分解為多個(gè)本征模態(tài)函數(shù)uk,并求得每個(gè)uk的單邊譜,同時(shí)對(duì)每個(gè)uk進(jìn)行調(diào)頻,最后添加約束條件得到一個(gè)帶約束的變分模型[23]:

(1)

式中:ωk是假設(shè)分解得到的各階模態(tài)分量的頻率中心;uk是第k個(gè)IMF模態(tài)分量;f是原始信號(hào);δ(t)為一脈沖函數(shù)。

通過拉格朗日函數(shù)將上述問題轉(zhuǎn)化為無約束變分問題[15]:

(2)

式中:α為懲罰因子;λ為拉格朗日乘子。再利用交替乘子方向算法求解,更新各IMF分量的中心頻率以及帶寬[15]:

(3)

(4)

1.2 GWO算法原理

GWO算法將灰狼種群捕食活動(dòng)中的追蹤、接近、圍捕以及進(jìn)攻行為模擬運(yùn)用到算法搜索尋優(yōu)過程中,其種群將個(gè)體分為4個(gè)等級(jí):α、β、γ和ω。通過前3種最好的個(gè)體α、β、γ引導(dǎo)群體向獵物方向(即最優(yōu)解)靠近,該算法原理比較簡單。相比其他算法,在運(yùn)算過程中需要調(diào)整的參數(shù)較少且比較容易實(shí)現(xiàn),且收斂速度較快,收斂結(jié)果較準(zhǔn)[23]。

灰狼種群追捕獵物的公式為

D=|C·Xp(t)-X(t)|

(5)

X(t+1)=Xp(t)-A·D

(6)

式中:C和A為系數(shù);Xp(t)為獵物位置;X(t)為灰狼位置。C和A的更新公式為

C=2r2

(7)

A=2ar1-a

(8)

式中:r1和r2為隨機(jī)系數(shù),其模值范圍為[0,1];a為收斂因子,隨迭代次數(shù)從2到0線性遞減。再由式(5)和式(6)求出3條頭狼的位置Xα(t+1)、Xβ(t+1)、Xγ(t+1),然后通過3條頭狼的位置確定種群中其他灰狼的位置:

(9)

GWO算法根據(jù)當(dāng)前解空間中前3個(gè)最優(yōu)解搜索全局最優(yōu)解,同時(shí)在搜索過程充分借助輔助系數(shù)改變某個(gè)體權(quán)重,因此算法陷入局部最優(yōu)的概率較低。

2 基于GWO-VMD的滾動(dòng)軸承故障特征提取方法

2.1 基于皮爾遜相關(guān)函數(shù)、最大峭度和時(shí)頻譜相似度的優(yōu)化目標(biāo)函數(shù)

原始信號(hào)分解后的信號(hào)不能丟失原始信號(hào)的關(guān)鍵信息特征,且分解后的信號(hào)應(yīng)更明顯地體現(xiàn)出其內(nèi)部的特征,否則做此分解操作并無意義。因此,量化分解效果的有效性應(yīng)有兩個(gè)準(zhǔn)則[24-26]:一是分解后的信號(hào)應(yīng)盡量包含原始信號(hào)完整信息;二是分解后應(yīng)更明顯地突出故障特征。

基于上述兩個(gè)準(zhǔn)則,本研究提出3個(gè)優(yōu)化目標(biāo):皮爾遜相關(guān)系數(shù)、時(shí)頻譜相似度以及最大峭度。

(1) 皮爾遜相關(guān)系數(shù)

皮爾遜相關(guān)系數(shù)可衡量兩種信號(hào)的線性相關(guān)程度,絕對(duì)值越大代表線性相關(guān)的程度越強(qiáng)。

首先將信號(hào)進(jìn)行VMD,然后通過得到的模態(tài)分量將信號(hào)重構(gòu),得到重構(gòu)信號(hào)X(t),Xi為重構(gòu)信號(hào)X(t)某一個(gè)樣本點(diǎn)。設(shè)原始信號(hào)為S(t),Si為原始信號(hào)S(t)某一個(gè)樣本點(diǎn),則可計(jì)算X(t)和S(t)的皮爾遜相關(guān)系數(shù)如下:

(10)

(2) 時(shí)頻譜相似度

時(shí)頻譜能夠描述兩種信號(hào)在時(shí)域上其頻率成分的變化情況,恰好解決了非平穩(wěn)信號(hào)在時(shí)域上其頻率不斷變化的問題,通過時(shí)頻譜相似度可以有效地反映兩種非平穩(wěn)信號(hào)的相關(guān)性。

利用希爾伯特黃變換求得重構(gòu)信號(hào)時(shí)頻譜Xt-f和原始信號(hào)時(shí)頻譜St-f,其時(shí)頻譜相似度為

(11)

(3) 最大峭度

峭度是一種無量綱參數(shù)。通常情況下滾動(dòng)軸承振動(dòng)信號(hào)峭度值為3左右,信號(hào)幅值呈正態(tài)分布,當(dāng)出現(xiàn)故障時(shí),振動(dòng)信號(hào)峭度值增大且信號(hào)幅值偏離正態(tài)分布。峭度值越大,表明該滾動(dòng)軸承故障越嚴(yán)重。

當(dāng)原始信號(hào)分解后,其各個(gè)模態(tài)分量中的最大峭度越大,說明分解后的信號(hào)故障特征變現(xiàn)越明顯,因此計(jì)算分解后各個(gè)本征模態(tài)分量的峭度,取最大的峭度:

(12)

兩個(gè)時(shí)域信號(hào)的皮爾遜相關(guān)系數(shù)只能表征兩個(gè)信號(hào)時(shí)域上的相關(guān)性,且在此背景下,皮爾遜相關(guān)系數(shù)只有正向判定而無反向判定,即數(shù)值高表征兩個(gè)信號(hào)高度相似,但數(shù)值低并不能表征兩個(gè)信號(hào)毫無相關(guān)性。因此,求得以上3個(gè)指標(biāo)后,引入閾值融合皮爾遜相關(guān)系數(shù)、時(shí)頻譜相似度和最大峭度,融合流程圖如圖1所示。

圖1 融合3個(gè)目標(biāo)流程圖

如圖1所示,當(dāng)兩種信號(hào)的皮爾遜相關(guān)系數(shù)大于某閾值時(shí),則可充分表明重構(gòu)信號(hào)包含原始信號(hào)的完整信息特征,因此融合最大峭度系數(shù),提高其對(duì)故障特征的敏感性,不必再進(jìn)行大量的計(jì)算求得時(shí)頻譜相似度。當(dāng)皮爾遜相關(guān)系數(shù)小于閾值,此時(shí)皮爾遜相關(guān)系數(shù)失效,必須求得時(shí)頻譜相似度來標(biāo)定重構(gòu)信號(hào)與原始信號(hào)的相關(guān)性。

2.2 基于GWO-VMD的滾動(dòng)軸承故障特征提取方法

在VMD算法中懲罰因子α是一個(gè)重要的參數(shù),其決定著VMD后各IMF的帶寬。通常,滾動(dòng)軸承的振動(dòng)信號(hào)包含諧波分量與脈沖分量,諧波信號(hào)頻帶帶寬較小,因此VMD時(shí)需要較大的α,有利于提取其特征;而故障脈沖信號(hào)相反,其時(shí)域較短,頻帶帶寬較大,因此VMD時(shí)需要較小的α,有利于提取其特征。

由此可見,單懲罰因子VMD無法同時(shí)提取滾動(dòng)軸承振動(dòng)信號(hào)包含的諧波特征與故障脈沖特征。因此,本文首先利用最大峭度準(zhǔn)則優(yōu)化VMD的分解層數(shù)K,然后通過GWO算法,以上文的優(yōu)化目標(biāo)函數(shù)為適應(yīng)度,為每個(gè)IMF分量搜尋最優(yōu)的懲罰因子,信號(hào)經(jīng)分解后,對(duì)各個(gè)本征模態(tài)分量作包絡(luò)譜分析,提取其故障頻率特征。算法具體流程圖如圖2所示。

算法具體步驟如下所示。

步驟 1首先需要確定分解層數(shù)K的值,設(shè)置K從2遞增到n,每次迭代中對(duì)振動(dòng)信號(hào)作VMD得到K個(gè)模態(tài)分量,分別求得每個(gè)本征模態(tài)分量的峭度值,并取最大值,取第一個(gè)達(dá)到極大值時(shí)對(duì)應(yīng)的K。

步驟 2設(shè)置GWO初始參數(shù)a,A、C的值。

步驟 4將每只灰狼的位置Xi作為懲罰因子[α1,α2,…,αK]對(duì)信號(hào)進(jìn)行VMD,根據(jù)式(10)～式(12)計(jì)算得出適應(yīng)度值,并記錄最優(yōu)的3個(gè)頭狼個(gè)體位置Xα、Xβ和Xγ。

步驟 5通過式(9)更新其他灰狼個(gè)體的位置,并通過式(7)和式(8)更新參數(shù)A、C和a的值。

步驟 6循環(huán)執(zhí)行步驟4和步驟5,直至迭代完畢。

步驟 7選取適應(yīng)度值最大的個(gè)體作為優(yōu)化后的懲罰因子[α1,α2,…,αK]以及最佳分解層數(shù)K作為VMD參數(shù)對(duì)信號(hào)作分解。

步驟 8對(duì)各個(gè)IMF做包絡(luò)解調(diào),提取故障特征信息。

3 實(shí)驗(yàn)與分析

3.1 數(shù)據(jù)源選擇

本實(shí)驗(yàn)數(shù)據(jù)來源于某電氣工程實(shí)驗(yàn)室的滾動(dòng)軸承數(shù)據(jù),是相關(guān)研究學(xué)者普遍采用的數(shù)據(jù)集。該實(shí)驗(yàn)采用一個(gè)參數(shù)可調(diào)節(jié)的電機(jī)作為負(fù)載,負(fù)載可為1 hp、2 hp以及3 hp。實(shí)驗(yàn)中的故障類型為表面加工的單點(diǎn)故障。實(shí)驗(yàn)中被測試的軸承有兩種,其中SKF型號(hào)軸承的故障直徑包含0.007 mils、0.014 mils、0.021 mils,NTN型號(hào)軸承的故障直徑包含0.028 mils以及0.040 mils。本研究選擇了故障直徑為0.021 mils的SKF型號(hào)軸承振動(dòng)信號(hào)數(shù)據(jù)來驗(yàn)證所提方法的有效性,其轉(zhuǎn)速為1 772 r/min,采樣頻率為12 kHz。

3.2 實(shí)驗(yàn)過程與結(jié)果分析

首先,根據(jù)模型算法求得最佳分解層數(shù),固定分解因子,計(jì)算出不同分解層數(shù)對(duì)應(yīng)的最大峭度值,最大峭度與分解層數(shù)的關(guān)系圖如圖3所示。由圖3可以觀察到,最大峭度值對(duì)應(yīng)的分解層數(shù)為4,因此最佳分解層數(shù)為4層。求出最佳分解層數(shù)后,運(yùn)用GWO-VMD算法,為每層模態(tài)求得最佳的懲罰因子。

圖3 最大峭度與分解層數(shù)的關(guān)系

本實(shí)驗(yàn)還將本算法與基于參數(shù)優(yōu)化的VMD(parameter optimized VMD,PO-VMD)、基于蝗蟲優(yōu)化算法的VMD(grasshopper optimization algorithm-VMD, GOA-VMD)兩種算法求解最優(yōu)值的收斂情況進(jìn)行對(duì)比,3種算法的收斂迭代折線圖如圖4所示。從圖4中可觀察到,GWO-VMD算法收斂效果最好,在47代左右基本收斂達(dá)到最優(yōu)值;PO -VMD算法雖在前期收斂較快,但最終在70代左右才基本收斂達(dá)到最優(yōu)值;GOA-VMD算法的效果則完全不如GWO -VMD算法。表1為各算法達(dá)到最優(yōu)值的收斂情況,可觀察到雖然PO-VMD算法可達(dá)到最快的收斂速度,但其整體性能遠(yuǎn)遠(yuǎn)不如GWO-VMD算法,而GOA-VMD算法各項(xiàng)指標(biāo)都不如GWO-VMD算法,因此可表明GWO-VMD算法的效率性強(qiáng)于其他算法。由算法求得4層對(duì)應(yīng)的懲罰因子分別為[5 224, 4 110, 1 860, 540]。

表1 3種算法收斂情況

將上述求得的最佳分解層數(shù)和懲罰因子代入,對(duì)信號(hào)作VMD,并通過分解后的各個(gè)模態(tài)分量重構(gòu)原始振動(dòng)信號(hào)數(shù)據(jù)如圖5和圖6所示。由圖5和圖6可觀察到,原始信號(hào)與重構(gòu)信號(hào)的時(shí)域波形及其相似。對(duì)這兩種信號(hào)作包絡(luò)解調(diào),如圖7所示,可觀察到原始信號(hào)與重構(gòu)信號(hào)的包絡(luò)譜圖極其相似,重構(gòu)信號(hào)包絡(luò)譜圖僅僅在某些頻率上的幅值有些許損失,因此表明該重構(gòu)信號(hào)包含原始信號(hào)的完整特征,證明了該方法的有效性。

圖5 原始振動(dòng)信號(hào)波形圖

圖6 重構(gòu)信號(hào)波形圖

圖7 原始信號(hào)與重構(gòu)信號(hào)包絡(luò)譜

再進(jìn)一步對(duì)信號(hào)作分解,并求出各個(gè)IMF分量的頻譜圖。GWO-VMD算法分解后各個(gè)IMF的時(shí)域波形圖以及對(duì)應(yīng)的頻譜圖如圖8所示。可觀察到,原始信號(hào)通過分解后,分離出了各個(gè)頻段部分的分量,分解效果較好。PO -VMD算法和VMD算法分解各分量時(shí)域波形和頻譜分別如圖9和圖10所示?？梢杂^察到,圖9中IMF2的頻譜圖與IMF4頻譜圖處有重疊,并且在低頻處丟失了信號(hào)應(yīng)有的頻率特征。圖10也有相似的特點(diǎn),圖10中IMF1的頻譜圖與IMF4的頻譜圖處有重疊,也丟失了信號(hào)本該有的低頻特征,這證明了本文方法的有效性以及優(yōu)越性。

圖9 PO-VMD各分量時(shí)域波形和頻譜

圖10 VMD各分量時(shí)域波形、頻譜

然后,再對(duì)3種算法分解得到的IMF分量作包絡(luò)解調(diào),得到各個(gè)IMF分量的包絡(luò)譜圖。GWO-VMD后各IMF的包絡(luò)圖如圖11所示?？梢钥吹?在故障頻率處以及其2～5倍頻處有明顯的譜峰,說明提取到了其故障特征頻率,證明本文提出方法的有效性,并且譜圖在故障頻率及其倍頻處比較干凈,干擾較少。

圖11 GWO-VMD后4個(gè)IMF的包絡(luò)譜

PO-VMD后的各IMF包絡(luò)譜圖如圖12所示。從圖12可以觀察到,通過PO -VMD后,雖在IMF4包絡(luò)譜圖中在故障特征頻率處以及其2～3倍頻處有譜峰,但是其譜峰不是很明顯,受周圍噪聲頻率干擾會(huì)比較大,并且在4倍頻以及4倍頻以上處基本看不到譜峰,抑制效果過于嚴(yán)重導(dǎo)致所需要的故障特征也被抑制,無法有效地表征信號(hào)中蘊(yùn)含的故障特征信息。

圖12 PO-VMD后各IMF的包絡(luò)譜

VMD后的各IMF包絡(luò)譜圖如圖13所示。從圖13可觀察到,VMD后各個(gè)IMF都出現(xiàn)了較大的噪聲,分解過程中噪聲抑制不明顯,并且分解后頻率成分幾乎集中在低頻中,該頻段噪聲干擾很大,無法有效地觀察到故障頻率特征。

圖13 VMD后各IMF的包絡(luò)譜

上述分析結(jié)果表明了本文方法提取效果更好,能夠更有效地提取出故障特征頻率。表2中列出了3種方法提取的特征頻率幅值情況,并做了歸一化方便對(duì)比。同時(shí),從表2也可以觀察到本文所提的GWO-VMD方法提取的特征頻率幅值遠(yuǎn)大于其他兩種對(duì)比方法,進(jìn)一步驗(yàn)證了本方法的有效性和優(yōu)越性。

表2 3種算法提取特征頻率幅值

4 結(jié) 論

本文提出了一種基于多懲罰因子優(yōu)化VMD的滾動(dòng)軸承故障特征提取方法,該方法以VMD理論為基礎(chǔ),以提升故障特征顯著性和包含信息完整性為目標(biāo),構(gòu)建融合皮爾遜相關(guān)系數(shù)、時(shí)頻譜相似度以及最大峭度值的優(yōu)化目標(biāo)函數(shù),并利用GWO算法自適應(yīng)地搜索VMD的最佳分解層數(shù)和每個(gè)IMF分量對(duì)應(yīng)的懲罰因子,解決了VMD需根據(jù)經(jīng)驗(yàn)手動(dòng)設(shè)置參數(shù)導(dǎo)致分解效果不佳的問題,實(shí)現(xiàn)了滾動(dòng)軸承故障特征的有效提取。實(shí)驗(yàn)結(jié)果分析表明,GWO-VMD與PO -VMD和VMD等算法相比,可以在噪聲環(huán)境下有效地提取軸承故障特征。