基于SSA-VMD-MCKD的強背景噪聲環境下滾動軸承故障診斷

任 良， 甄龍信， 趙 云， 董前程， 張云鵬

(燕山大學 河北省特種運載裝備重點實驗室，河北 秦皇島 066004)

滾動軸承是機械設備的關鍵部件，應用非常廣泛，實際生產中其產生的故障往往會造成非常嚴重的后果[1]，其工作環境通常比較惡劣，故障信號特征往往被背景噪聲和機械轉動信號覆蓋，不易提取[2]，因此，需要發展在強背景噪聲干擾下的故障診斷方法。

滾動軸承故障信號通常呈現非線性、非平穩的特點且表現為周期性沖擊特征。隨著信號處理技術的發展，滾動軸承故障診斷的研究取得一定成效。文獻[3]利用經驗小波分解(empirical wavelet transfor, EWT)提取故障特征，再通過獨立分量分析(independent component analysis, ICA)算法降噪并成功進行軸承故障診斷。文獻[4]利用小波分解與FastICA算法對軸承信號聯合降噪，然而這種經典降噪算法存在小波基難以選擇的問題，且實際操作中FastICA算法用于軸承信號降噪的效果一般。經驗模態分解(empirical mode decomposition, EMD)是由Huang等[5]提出的一種自適應分解的時頻分析方法，但其存在模態混疊的缺陷，為減弱這種缺陷，文獻[6]提出集成經驗模態分解(ensemble empirical mode decomposition, EEMD)。文獻[7]先采用最大相關峭度解卷積(maximum correlated kurtosis deconvolution，MCKD)算法降噪，再用EEMD分解，成功提取柔性薄壁軸承故障特征信號。文獻[8]提出變分模態分解(variational mode decomposition，VMD)算法，相比EMD缺乏嚴格的理論推導，VMD具有嚴格的數學理論支持且可避免模態混疊的產生。文獻[9]將VMD成功應用于提取軸承故障信號。然而VMD的二次懲罰因子alpha和分解層數K等參數取值對VMD分解結果影響較大，故文獻[10]采用相對熵最小化作為目標函數對VMD進行參數優化，優化思路為先設定一個alpha尋最優K值，再利用最優的K值尋最優alpha值，該算法存在難以找到全局最優參數的弊端。

文獻[11]提出MCKD算法，該算法通過構建和尋找最優濾波器以突出信號中的連續周期脈沖信號，但濾波器長度L、移位數M等參數對MCKD算法的效果影響較大，文獻[12]采用粒子群算法對MCKD和VMD進行參數優化，并成功提取了滾動軸承早期微弱故障特征。文獻[13]提出一種處理非線性非平穩信號的新方法——奇異譜分析(singular spectrum analysis，SSA)，通過對信號軌跡矩陣的奇異值分解，將信號重構成不同趨勢的分量，奇異值分解可有效將信號和噪聲分量分開，故該算法可用于信號降噪[14-15]。

為充分利用SSA、VMD、MCKD在降噪和特征提取方面的優勢，本文提出將三者結合用于強背景噪聲環境下的滾動軸承早期故障診斷，為提高診斷準確性，以包絡熵極小值為目標函數，利用鯨魚優化算法(whale optimization algorithm，WOA)優化VMD和MCKD參數。

1 相關理論

1.1 SSA算法

SSA是一種處理非線性非平穩時間序列信號的方法，可實現信號的去噪、動態重構和特征提取。本文利用SSA進行信號去噪，其核心在于利用奇異值分解原理將特征信號與噪聲分離，主要步驟為如下。

步驟1輸入一維時域信號[s1,s2,s3,…,sn]。

步驟2確定合適窗口長度l,得到如下軌跡矩陣

(1)

步驟3進行奇異值分解，將S分解為如下形式

S=αΣVT

(2)

式中：α為l×l的矩陣;Σ為l×(n-l+1)的矩陣;V為(n-l+1)×(n-l+1)的矩陣。

計算軌跡矩陣的協方差矩陣

C=SST

(3)

(4)

步驟4分解所得的l個分量代表著不同的趨勢成分，本文通過求l個分量與源信號的互相關系數進行篩選重構。

1.2 VMD算法

VMD算法將信號x(t)分解成K個互不相關的稀疏子信uk(t)，具體分解步驟如下。

步驟1首先定義本征模態函數(intrinsic mode function，IMF)

uk(t)=Ak(t)cos[φk(t)],k=1,2,…

(5)

步驟2將uk(t)進行Hilbert變換得到單邊譜

(6)

式中:δ(t)為狄利雷克分布函數;*為卷積。

步驟3將信號頻譜移頻至基帶上

(7)

步驟4基于梯度二范數的平方估計各個uk(t)的帶寬，保證所有IMF帶寬之和最小

(8)

L({uk},{wk},λ)=alpha·

(9)

式中：alpha為二次懲罰因子;η(t)為拉格朗日乘子。

步驟5引進交替方向乘子算法對式(8)、式(9)進行迭代尋優，可得K個IMF，如式(10)所示

(10)

由式(5)可知K影響瞬時頻率的預估，不恰當的K取值將導致模態混疊的產生，由式(9)和式(10)可知alpha值會影響IMF的帶寬，故參數alpha和K的值對VMD分解效果影響較大[16]。

1.3 MCKD算法

MCKD算法能有效提取周期性沖擊脈沖分量，算法的本質是在構建最合適的濾波器用于突出被噪聲掩蓋的特征信號。假設產生周期性沖擊脈沖信號s=[s1,s2,…,sn],信號在傳遞過程中會混有各類噪聲。觀測信號的表達式為

x=h·s+n

(11)

式中：x為觀測信號;s為原始周期脈沖信號;h為系統傳遞函數;n(t)為噪聲。

構建濾波器f=[f1,f2,…,fL]，以相關峭度最大為目標函數，尋找最優濾波器，用于分離出原始周期脈沖信號，先忽略噪聲，則表達式為

(12)

式中,L為濾波器長度。

相關峭度的表達式為

(13)

式中：T為沖擊信號周期;M為移位數。

為了求得最大的相關峭度，對式(13)進行求導即

(14)

從式(14)中解得最優濾波器f

(15)

其中

算法的流程如下。

步驟1輸入參數L,M和T。

步驟3求得經濾波器輸出的信號s。

步驟4由s計算出β和Ψm。

步驟5更新濾波器f′。

步驟6若信號相關峭度值CKM(T)達到最大值，則停止迭代，否則重復步驟3～步驟5，直至滿足CKM(T)最大為止并輸入此時對應的s值。

顯然，參數L，T和M的取值對MCKD算法結果影響較大，M的取值大于7時會使分解精度降低，一般取1～7。T可通過式(16)確定

(16)

式中：fs為采樣頻率，fi為故障特征率。

1.4 基于WOA算法的VMD和MCKD參數優化

本文利用WOA[17]尋找VMD的參數alpha和K以及MCKD算法的參數L和M的最優值，均以包絡熵極小值為目標函數，包絡熵越小則代表信號中含特征信號越多。

信號x(i)的包絡熵Ep的計算如式(17)所示

(17)

式中：a(i)為原信號x(i)經Hilbert解調后的包絡信號;ε(i)為a(i)的歸一化形式;N為信號x(i)的長度。

鯨魚群體狩獵時，鯨魚個體有兩種捕食策略:一種是直接包圍獵物;另一種是環形游動并產生汽泡形成氣泡網驅趕獵物至氣泡網中心，算法流程如下。

步驟1鯨魚種群個體數量、位置、進化代數等參數初始化，第i個鯨魚個體位置為

Xi=r·(ub-lb)+lb

(18)

式中：r∈[0 1];Xi的取值范圍為[lb,ub],lb為待尋優參數的下邊界，ub為上邊界。

步驟2p<0.5且|A|<1時，WOA按照式(19)搜索

(19)

式中：r1,r2,p為[0 1]的隨機數;t為迭代次數;tmax為最大迭代次數;X(t)為當前解位置;X*(t)為當前最優解位置;A,C為系數。

步驟3當p<0.5且|A|≥1時，WOA根據式(20)進行隨機搜索迭代更新

(20)

式中，Xrand為隨機選取鯨魚個體位置。

步驟4當p≥0.5時，WOA根據式(21)進行螺旋收縮方式迭代更新

(21)

式中：D為個體與獵物之間的距離;b為螺線常數;l∈[-1,1]且為隨機數。

步驟5若t>tmax則迭代停止，輸出最優尋優結果，若不滿足則返回步驟2。

WOA優化VMD和MCKD算法時，目標函數O為包絡熵Ep極小值，即作為WOA算法的適應度函數。表達式如下

o=min{Ep}

(22)

參數優化的表達式如下

(23)

式中,lb和ub為待優化的參數的取值邊界。

具體流程如圖1所示。

圖1 WOA優化VMD，MCKD流程圖

2 故障診斷流程

本文提出采用SSA-VMD-MCKD方法提取強背景噪聲干擾下的滾動軸承微弱故障信號特征并進行故障診斷，步驟如下。

步驟1先將信號經SSA分解降噪，采用時域互相關準則將分解產生的不同趨勢的子信號篩選重構。時域互相關準則是衡量兩個時域信號相關程度的有效指標，計算公式如下

(24)

步驟2重構信號經WOA優化的VMD分解提取故障特征，WOA優化時設定的參數大小為：種群規模為100，最大迭代次數為10，alpha的尋優范圍為[0,10 000]，K的尋優范圍為[2,10]。采用峭度準則，將分解產生的IMF篩選重構，通常軸承正常運轉時信號的峭度值接近3，大于3時則認為含有較多的沖擊成分，故本文選擇值大于3的IMF進行重構，峭度值計算公式如下

(25)

式中：d為信號長度;u為信號均值;σ為標準差。

步驟3對重構信號包絡解調，確定最大幅值所對的特征頻率視為fi，根據式(16)計算MCKD的參數T，重構信號經WOA優化的MCKD分析，以加強故障特征從而確定故障。WOA優化時設定的參數為：種群規模為100，最大迭代次數為10，L的尋優范圍為[20,500]，M的尋優范圍為[1,7]。

步驟4對經MCKD處理的信號進行包絡解調，診斷故障。

總結以上診斷步驟，繪制故障診斷流程圖如圖2所示。

圖2 基于參數優化的SSA-VMD-MCKD的故障診斷步驟

3 仿真信號分析

3.1 單一軸承仿真信號故障診斷

為驗證本文所提方法的有效性，構建單一軸承周期性沖擊故障仿真信號，表達式如下

(26)

式中：s(t)為周期性沖擊故障信號分量;fn為固有頻率;A0為位移常數;g為阻尼系數;n(t)為高斯白噪聲信號。

取fn=1 600，A0=0.8，g=0.12，沖擊故障的時間間隔T為0.006 25 s,故特征頻率fi=1/T=160 Hz,采樣頻率為12 kHz,采樣點數為7 800。未添加高斯白噪聲的故障信號的時域圖，如圖3所示。由圖3可看出明顯的周期沖擊特征，采用短時傅里葉分析，畫出信號的時頻圖，如圖4所示，可知圖4不存在零星的噪點，表明不含噪聲。沖擊故障的包絡譜圖，如圖5所示。由圖5可清晰看出故障特征頻率fi及其倍頻。

為模擬強背景噪聲環境，添加-14 dB的高斯白噪聲干擾。添加-14 dB噪聲后的故障信號的時域圖，如圖6(a)所示。與圖3相比，可知信號特征已被噪聲嚴重污染。強背景噪聲軸承信號的時頻圖，如圖6(b)所示。可以看到圖中存在大量噪點，故障信號的頻譜圖和包絡譜圖，如圖6(c)、圖6(d)所示。顯然，信號特征頻率及其倍頻已被強背景噪聲掩蓋無法識別。

(a) 添加-14 dB噪聲的沖擊信號時域圖

采用本文所提SSA-VMD-MCKD方法對仿真故障信號進行處理，先將信號進行SSA分解，窗口長度l值越大，算法運行的時間越長且不合適的取值會破壞故障特征，經多次驗證窗口長度l=13時算法的效率較好，即將信號分解成13個不同趨勢的子信號，計算每個子信號與原信號的時域互相關系數，結果如圖7所示。通常相關系數大于0.5表示與原信號相似性較高，為避免故障特征信號丟失，本文設置閾值為0.5，故選取子信號2、信號4、信號5、信號7、信號8、信號11進行信號重構。再利用WOA優化算法尋找VMD的alpha和K的最優值，WOA優化算法的種群規模設為100，最大迭代次數為10，參數優化的目標函數即適應度函數為包絡熵極小值，alpha的尋優范圍為[0,10 000],K的尋優范圍為[2，10],計算機參數為Intel(R) Core(TM) i5-10300H CPU@2.50 GHz,運行內存為16 GB,利用MATLAB2014a處理，尋得最優參數組合為[12,4],故設定VMD分解的alpha=12，K=4，對重構信號利用VMD分解以提取故障特征，分解所得IMF的峭度值如圖8所示。由圖8可知IMF1和IMF4峭度值大于3，故選取這兩個信號進行重構。

圖7 SSA分解各子信號的互相關系數統計圖

圖8 IMF的峭度統計圖

重構信號的時域圖和頻譜圖如圖9(a)、圖9(b)所示，與圖6(a)、圖6(c)對比可知經SSA-VMD分解后信號中的毛刺和噪聲頻帶減少很多，表明降噪效果明顯，為了更直觀看出降噪效果，作出時頻圖，如圖9(c)所示，對比圖6(b)可知時頻圖上的中高頻噪點消除很多，即大部分中高頻噪聲被濾除，但仍存在少量噪聲。經參數優化的VMD分解后的重構信號的包絡譜圖，如圖9(d)所示。選擇圖中最大峰值處頻率作為特征頻率fi，故fi=160，根據式(16)計算出MCKD算法的參數T=75；采用WOA優化算法優化MCKD的參數L和M，WOA算法的種群數設為100，最大迭代次數為10，參數優化的目標函數即適應度函數為包絡熵值極小值，參數L的尋優范圍為[20,500],M的范圍為[1,7],尋得最優參數為[231,4],故設定MCKD算法的參數L=231，M=4，T=75。再將VMD分解后的重構信號經MCKD算法分析以加強故障特征，最后對信號進行包絡解調，診斷故障。最終所得信號的包絡譜圖，如圖10所示。顯然，圖中特征譜線十分突出，可輕易識別故障的特征頻率fi及其2倍～6倍頻，從而證實了所提診斷方法的有效性。

(a) VMD分解后重構信號時域圖

圖10 本文方法處理后信號包絡譜圖

3.2 驗證所提方法的合理性與必要性

為了驗證采用WOA算法優化MCKD的參數L和M的合理性，現改變L的參數，保持M不變觀察結果，即隨機令L=142，M值不變，經MCKD算法分析的信號包絡譜圖，如圖11所示。從圖11可以看出,特征頻率及其倍頻周圍存在較多譜線干擾，且5倍頻未成功提取。再保持L值不變，改變M值，即隨機取M=1，L值不變，經MCKD算法分析的信號包絡譜圖，如圖12所示。圖12已無法發現故障特征頻率及其倍頻，表明本文采用WOA優化MCKD的參數L和M結果的合理性。

Sugita等已證明VMD與MCKD結合的合理性，為證明強背景噪聲下SSA與VMD-MCKD結合的必要性，現將仿真信號直接通過WOA優化的VMD和MCKD算法進行分析，為了使對比試驗更具對比性及說服力，WOA優化的種群規模、最大迭代次數以及VMD和MCKD的參數尋優范圍均與之前設置相同。經參數優化VMD-MCKD算法分析的信號包絡譜圖，如圖13所示。顯然，圖13突出的譜線頻率為146 Hz，并非故障特征頻率，表明強背景噪聲下，以包絡熵為目標函數，直接通過參數優化的VMD-MCKD算法分析，故障特征提取效果不佳，故間接證明了本文采用SSA結合參數優化的VMD-MCKD的必要性，強背景噪聲干擾下，本文采用SSA分解的目的是降低噪聲以便于特征提取。

4 實測信號分析

4.1 實測軸承信號故障診斷

采用凱斯西儲大學公開的滾動軸承故障模擬平臺試驗數據進行分析。試驗平臺如圖14所示，電機驅動端安裝SKF6205-2RS型深溝球軸承，其參數如表1所示。采用電火花加工技術設置故障，電機轉速為1 772 r/min,采樣頻率為12 kHz，本文選取7 800點進行分析。電機的轉頻fr及內圈故障特征頻率fi的計算公式如下

圖14 滾動軸承故障模擬試驗平臺

(27)

式中,n為電動機轉速。

依據表1參數，由式(27)可計算故障特征頻率的理論值為fi=159.93 Hz。

表1 滾動軸承結構參數

為了模擬強背景噪聲環境，將試驗信號中再添加-5 dB的高斯白噪聲。信號的時域圖、頻譜圖和包絡譜圖，如圖15(a)～圖15(c)所示。由圖可知故障特征無法辨別且特征頻率及其倍頻已完全被轉頻及其倍頻以及噪聲干擾頻帶掩蓋無法進行故障診斷。故障信號的時頻圖，如圖15(d)所示。可觀察到圖中噪點密布,表明含噪較多。

(a) 實測信號時域圖

采用本文所提方法，先將信號經SSA分解，窗口長度l設為13，分解所得子信號與原故障信號的時域互相關系數如圖16所示。取相關系數大于0.5的子信號重構，即取子信號1、信號2、信號3、信號4、信號6、信號7進行信號重構。采用WOA優化VMD參數alpha和K，計算機參數與3.1節相同，WOA的種群規模設為100，最大迭代次數為10，參數優化的目標函數即適應度函數為包絡熵極小值，alpha的尋優范圍[0,10 000],K的尋優范圍為[2,10],最優參數組合為[2 017,5]，故設置alpha=2 017,K=5,將重構信號通過VMD分解，分解所得的各個IMF的峭度值如圖17所示，可知選取IMF1,IMF2,IMF3,IMF5進行信號重構。

圖16 SSA分解各子信號的互相關系數統計圖

圖17 IMF的峭度統計圖

重構信號的時域圖、時頻圖，如圖18(a)、圖18(b)所示。圖中可觀察到沖擊信號特征且時頻圖中噪點大幅減少，表明所提方法降噪效果明顯。重構信號的包絡譜圖，如圖18(c)所示。取圖中最大幅值處所對的頻率為特征頻率fi,其值為160，根據式(16)計算出MCKD算法的參數T的值，可得T=75，利用WOA優化MCKD，WOA算法的種群規模設為100，最大迭代次數為10，L的尋優范圍為[20,500],M的尋優范圍[1,7],尋優耗時為1 726.306 3 s，最終所得最優參數為[399，4],將重構信號經MCKD算法，按照最優參數設置，L=399，M=4。經參數優化的MCKD算法處理所得信號的包絡圖，如圖18(d)所示。圖18(d)故障的特征頻率及其倍頻顯而易見，故可確定滾動軸承故障位置在軸承內圈。

(a) VMD分解重構信號時域圖

4.2 與EWT-ICA方法結果比較

呂躍剛等[3]對故障信號先進行EWT分解，計算分解所得IMF的峭度值和時域互相關系數，選取峭度值及相關系數值都較大的IMF作為快速獨立分量分析(FastICA)算法的觀測信號，其余信號則作為該算法的虛擬噪聲通道信號，經FastICA算法解混后，再將特征信號包絡解調，此方法在弱背景噪聲環境下，故障診斷效果較好。為驗證本文方法的優越性，采用呂躍剛等的方法對添加-5 dB高斯白噪聲的實測故障信號進行分析，并將結果與本文所提方法對比。首先，將信號通過EWT分解采用“locmax”準則分解為5個IMF,每個信號的相關系數及峭度值如表2所示，選取IMF1,IMF2,IMF5作為觀測信號，然后將剩余IMF作為虛擬噪聲通道信號輸入FastICA算法，最后將EWT-ICA算法分解所得源信號包絡解調，包絡譜如圖19所示。可知圖19只出現特征頻率160 Hz，其倍頻并未出現，表明EWT-ICA算法在強背景噪聲環境下無法完全提取故障特征相關信息，間接表明本文方法的優越性。

表2 各IMF的峭度值和相關系數

圖19 EWT-ICA方法所得包絡譜圖

5 結 論

(1) 本文提出基于SSA-VMD-MCKD算法的滾動軸承故障診斷方法，能在強背景噪聲環境下有效提取軸承故障特征，準確診斷軸承故障。

(2) SSA算法能有效抑制噪聲干擾，MCKD算法能有效加強周期沖擊特征，強背景噪聲環境下，其他條件相同時，基于SSA-VMD-MCKD算法的診斷結果優于VMD-MCKD算法。

(3) VMD中的alpha、K參數及MCKD的L,M參數的取值對軸承故障診斷結果影響較大，采用WOA算法優化這4個參數能有效提高故障診斷效果。

(4) 強背景噪聲環境下，基于SSA-VMD-MCKD算法的滾動軸承故障診斷結果優于EWT-ICA算法。