基于VMD-SWT滾動軸承故障診斷方法研究＊

周進群，劉義亞

（1.深南電路股份有限公司，深圳 518000；2.無錫深南電路股份有限公司，江蘇 無錫 214000）

前言

隨著社會的不斷發展進步，工業設備日漸趨于智能化、大型化，這就給故障診斷領域提出了嚴峻的挑戰。滾動軸承在所有機械設備中占據著舉足親重的地位，是機械設備中最重要的結構之一，一旦損壞，輕則導致巨大的經濟損失，重則危及人生安全戰[1]。滾動軸承的工作機理以及所處的工作環境，使得故障信號的信噪比極低，而早期故障又十分微弱，故而導致故障特征難以提取[2]。現階段較為成熟的分析算法如Hilbert包絡、短時傅里葉、小波變換、Wigner-Ville分布等都有各自的不足。如短時傅里葉窗長度固定、小波變換會造成能量泄漏、Wigner-Ville分布存在交叉項等。2014年，Konstantin Dragomiretskiy等在研究維納濾波的基礎上，參考經驗模式分解(Empirical Mode Decomposition ，EMD)提出了變分模態分解(Variational Mode Decomposition，VMD)[3]。該方法通過在變分框架內求解最優變分模型來求取模態分量。WangYanxue等對其進行研究，并使用VMD對轉子碰模故障進行分析診斷，取得了良好的效果[4]錯誤！未找到引用源。。TripathyRK等深入研究VMD并將其引入醫學心電圖研究中，提高了診斷的準確率[5]。李志農等人將其應用于機械的故障診斷中，驗證了該方法能夠揭示出碰模故障數據的頻率結構，區分碰模故障的嚴重程度[6]。馬增強等將Teager能量算子引入VMD中，對滾動軸承故障進行精確診斷，取得了良好效果[7]。但VMD的結果依賴于預先給定的模態數K，K的選取決定了最終的分解結果，K過小，則無法有效反應信號的特征信息，而K過大，則會導致頻率混疊、過分解等現象，因此K的選擇尤其重要。

2011年，Ingrid Daubechies等提出了同步壓縮變換(Syn-chrosqueezing transform，SST)[8]，同時結合小波變換所構成的同步壓縮小波變換(Synchrosqueezing Wavelet Transform，SWT)方法能夠通過對小波變換后的小波系數進行重組，將相同頻率下的尺度相加，從而抑制能量泄露，提高時頻分辨率。Gaurav Thakur等對其進行深入研究并將其引入古氣候分析中，取得理想效果[9]。Roberto HHerrera等將此方法用于地震波檢測分析中，提取出較為清晰的波形特征[10]。SWT能夠有效處理頻率較低及緩慢變換的信號中的噪聲，但當信號處于強噪聲中時，該方法難以達到理想效果。而在實際工程中，滾動軸承的工作環境及工作機理會造成提取的信號中夾雜著較大的噪聲，導致無法有效提取信號的時頻特征。

峭度是描述隨機變量分布特性的數值統計量，反應信號沖擊成分的大小。本文將最大峭度指標引入 VMD，選擇出最優分解模態數及懲罰因子 α，再利用峭度指標選擇有效VMD分解分量從而從噪聲信號中提取出有效信息，最后使用 SWT提取信號中的時頻特征，從而有效判斷滾動軸承的運行狀態。相比于EMD等方法來說，VMD具有嚴謹的理論支撐，且是根據中心頻率選取模態數，提高了計算效率，能夠將原始信號進行分解從而獲得K個含有不同頻率成分的模態分量；SWT方法雖然能夠提高時頻譜的可讀性，但當噪聲較大時仍然無提取出有效特征，而工業現場會由于環境、滾動軸承工作機理等會造成故障振動信號中包含很大的噪聲，導致信噪比極低，將VMD引入SWT中，能夠將原始信號進行分解，從而提取出含有有效頻率成分的模態分量，提高信噪比，再使用 SWT進行分析，從而有效提取信號的時頻特征。

1 變分模態分解VMD

VMD方法參考了 EMD算法，借助了內稟尺度分量（Intrinsic Scale Component，ISC）的概念，重新定義為一個調頻調幅信號[11]。

式中，φk(t)非遞減，；包絡線非負，Ak(t)≥0；瞬時頻率；Ak(t)及ωk(t)相對于φk(t)來說是緩變的。因此，在區間內，uk(t)可看作是幅值為 Ak(t)且頻率為ωk(t)的諧波信號。

可以通過如下步驟構造出變分模型：1）借助Hilbert變換獲得各模態函數uk(t)的解析信號，從而獲得信號的單邊頻譜；2）通過指數修正，將每個圍繞各自估算的中心頻率的模態函數調制到響應基頻帶；3）通過高斯平滑解調信號獲得每段帶寬，即L2范數梯度的平方根。受約束的變分模型為：

引入懲罰因子α構造增廣Lagrange函數從而求取上述變分模型的最優解，Lagrange構造函數如下：

式中：α為懲罰因子；λ為Lagrange乘子。

將上述Lagrange函數從時域變換到頻域，求解相應的極值，可得到uk和ωk的頻域表達式。

再利用交替方向乘子算法(Alternate Dirsplimi Method of Multipliers，ADMM)求取變分模型的最優解，從而實現VMD分解。算法具體步驟如下：

由上述步驟即可將信號分解成K個模態分量。

2 同步壓縮小波變換

峭度(Kurtosis)是反映振動信號分布特性的數值統計量，是無量綱參數，對信號的沖擊十分敏感，適用于表面損傷類故障，尤其是早期故障的診斷[12]：

在軸承正常運轉狀態下，其振動信號的幅值分布接近正態分布，當故障開始出現時，信號幅值的分布會偏離正態分布，峭度值也隨之增大。峭度指標的絕對值越大，說明軸承偏離其正常狀態，故障越嚴重。

同步壓縮小波變換(SWT)能夠在一定程度上有效抑制噪聲，改善小波變換的能量泄露，提高時頻分辨率[13]。

利用CWT對信號x(t)進行處理，可定義為：

此時便可將任意頻率 ωl周圍上的值壓縮到 ωl上，從而獲取同步壓縮變換的值，提高時頻分辨率。同步壓縮變換可定義為：

同步壓縮小波變換的逆變換(ISWT)定義為：

其中

改進SWT的算法(VSWT)流程如下：

（1）初始化VMD參數：分解層數K=2(K∈[2,15])，懲罰因子 α=10(α∈[10,2000])；

（2）計算各模態分量的峭度值，選出最大值，令K=K+1，重新計算峭度值，以此類推算出所有K情況下各模態分量的峭度值，選出所有K情況下峭度的最大值進行比較，最大的峭度值所在的K即為最優值；α=α+10；

（3）同理計算出最優α值(α步長為10)；

（4）將計算出的最優K和α帶入VMD中，使用VMD將信號進行分解；

（5）計算各模態峭度值，選出含有最大峭度值的模態分量；

（6）使用SWT對選出的模態分量進行分析。

本文提出使用最大峭度法提取出VMD的最優參數，將VMD引入SWT算法中構成一種新的分析方法(VSWT)，首先使用峭度指標判斷出VMD最優參數，再將提取出的參數帶入VMD中對原始信號進行分解，最后領用峭度指標提取有效模態分量后使用 SWT進行處理分析，提取出信號的時頻特征，從而對滾動軸承運轉狀態進行判斷。

3 仿真分析

通過構造一組仿真故障信號ft進行分析以便驗證提出方法的有效性。該信號構造軸承故障的仿真信號，其中軸承的固有頻率f=3000Hz，位移常數y0=5，阻尼系數ζ＝0.1，沖擊故障發生的周期為 0.01s，采樣頻率 fs=20kHz，采樣點數N=4096：

對信號加入信噪比(SNR)為-5的高斯白噪聲，圖1a和1b分別給出了原始信號和加噪信號幅值譜。首先設定VMD的初始化模態數為 K=2，懲罰因子α=10，對信號進行 VMD分解，計算各模態的峭度值，選擇最大的峭度值作為K=2時的峭度，對K進行疊加直至K=15停止，比較各模態數下的最大峭度值，圖1c給出了各模態數下的最大峭度值，從圖中可以看出K=7時峭度最大，因此模態數K定為7。同理對懲罰因子進行優化，圖1d給出了各懲罰因子下的最大峭度值，從圖中可以看出當 α=20時峭度值最大。由上述步驟確定VMD的模態數K=7，懲罰因子α=20，使用VMD對信號進行分解，圖1d給出了VMD的分解后各模態分量的幅值譜。計算各模態分量的峭度值，第一層峭度最大，達到了4.9201，因此選擇第一層模態分量為有效分量，使用 SWT對其進行處理，圖2給出了SWT處理結果。

圖1 VMD分解Fig.1 Variational Mode Decomposition

圖2(a)給出了原始信號的SWT處理結果，可以看出2000至4000Hz之間，以3000為中心頻率有一條清晰的能量集中帶，圖 2(b)給出了加噪信號的 SWT處理結果，反應出特征頻率已經完全淹沒在噪聲之中，圖 2(c)給出了對加噪信號使用VSWT處理的結果，從圖中可以反映出在3000附近出現了很強的能量集中帶，其他位置的噪聲已基本濾除，充分說明了本文所提出方法的可行性與有效性。

圖2 SWT處理結果Fig.2 Processing results by SWT

圖3給出了使用從圖3中可以反映出除邊緣部分有少量差別以外，其余部分高度吻合，反映了VSWT具有較高的重構能力。

圖3 原信號與重構信號對比Fig.3 The original signal and the signals reconstructed by IVSWT

綜上所述，可以判斷出VSWT具有良好的抗噪能力以及較高的信號提取精度，能夠在強噪聲背景下有效提取信號時頻特征，提高時頻可讀性，同時能夠有效的重構信號。

4 實驗和驗證

4.1 實驗說明

為驗證改進 SWT的工程實用性，本文以深南電路鉆孔工序的鉆床為實驗平臺，通過更換故障軸承采集滾動軸承正常振動信號、外圈故障信號、內圈故障信號、滾動體故障信號，分別應用 Hilbert包絡、SWT、VSWT分析故障信號，通過提取的時頻特征判斷滾動軸承的故障類型。

通常滾動軸承常見的故障類型為外圈故障、內圈故障以及滾動體故障。其故障特征頻率計算公式為：

內圈故障特征頻率

外圈故障特征頻率

滾動體故障特征頻率

式中，fr為轉頻；d為滾珠直徑；D為節圓直接； 為接觸角；Z為滾珠數。

圖4 軸承故障類型Fig.4 The type of bearing fault

表1 軸承參數Tab.1 Parameter of bearing

實驗平臺采用PCB MA352A60型加速度傳感器，采集垂直方向的振動信號數據，實驗滾動軸承具體參數見表 1。該實驗轉速設定為1000r/min，采樣頻率fs為50 kHz。以工程實際中滾動軸承故障現象為參照，利用線切割加工技術，在試驗臺滾動軸承外圈加工0.3×0.05(寬×深)微小凹痕模擬外圈故障；在試驗臺滾動軸承內圈加工0.3×0.05(寬×深)微小凹痕模擬內圈故障；在試驗臺滾動軸承內圈加工0.3×0.05(寬×深)微小凹痕模擬滾動體故障，具體現象如圖4所示。由式(18)～式(20)可求得其外圈、內圈和滾動體的故障特征頻率分別為88.64Hz、128.03Hz以及44.32Hz。

4.2 信號處理及分析

對采集的故障信號進行處理，圖5給出了滾動軸承正常振動信號、外圈故障振動信號、內圈故障振動信號以及滾動體故障振動信號的時域波形圖。

圖5 軸承振動信號Fig.5 The bearing vibration signal

為驗證改進 SWT方法的有效性與實用性，本文使用不同方法處理滾動軸承故障數據，并對處理結果進行比較。

圖6給出了外圈故障信號處理結果。使用Hilbert包絡對故障故障信號進行處理，結果如圖6a所示，含有噪聲較大，特征頻率完全被噪聲淹沒，從圖中找不出任何特征。圖 6b為 SWT處理后的時頻譜，同樣，從圖中找不到任何特征。使用本文提出的改進SWT(VSWT)進行處理，首先確定VMD分解的模態數和懲罰因子，令K=2，α=10，使用VMD進行分解，求取個模態的峭度值，選取最大的值作為K=2時的最大峭度值，保持α不變，令K=K+1(K∈[2,15])，同樣求出此時的最大峭度值，以此類推，求出不同K值下的最大峭度值，選擇所有最大峭度值中的最大峭度值所在的模態數K，此時K即為最優解，同理求取最優懲罰因子 =10( ∈[10,2000])，從圖6d～圖6e中可以看出當K=3，α=180時，峭度最大，分別為62.72和64.36，確定VMD參數后使用VMD對信號進行分解，經計算第3層峭度值最大，達到了64.3655，故選擇第3模態作為SWT輸入信號，使用SWT進行處理，結果如圖6c所示，從圖中可以看出，在89.27Hz、180.03Hz、265Hz、355.09Hz處都出現了清晰的能量集中帶，與外圈故障特征頻率的1倍頻88.64Hz、2倍頻177.28Hz、3倍頻265.92Hz、4倍頻354.56Hz較為接近，由此可以判斷滾動軸承外圈出現了故障，與實驗結果一致。

圖6 外圈故障信號處理結果Fig.6 Processing results of outer ring fault signal

針對內圈故障信號，使用上述的同樣方法進行分析。圖7a和圖7b分別給出了Hilbert包絡后的包絡譜及SWT后的時頻譜，由于采集過程中噪聲較大以及其他因素干擾，從圖無法找出任何特征。使用VSWT進行處理，首先確定VMD分解的模態數和懲罰因子，令K=2，α=10，使用VMD進行分解，求取個模態的峭度值，選取最大的值作為K=2時的最大峭度值，保持α不變，令K=K+1(K∈[2,15])，同樣求出此時的最大峭度值，以此類推，求出不同K值下的最大峭度值，選擇所有最大峭度值中的最大峭度值所在的模態數K，此時K即為最優解，同理求取最優懲罰因子 =140( ∈[10,2000])，從圖7d～圖7e中可以看出當K=9，α=140時，峭度最大，分別為67.93和72.78，確定VMD參數后使用VMD對信號進行分解，經計算第9層峭度值最大，達到了72.7765，故選擇第9模態作為SWT輸入信號，使用SWT進行處理，結果如圖 7(c)所示，從圖中可以看出，雖然周圍仍存在噪聲，但在130.27Hz、255.89Hz處都出現了清晰的能量集中帶，與外圈故障特征頻率的1倍頻128.03Hz、2倍頻256.06Hz、較為接近，由此可以判斷滾動軸承內圈出現了故障，與實驗結果一致。

圖7 內圈故障信號處理結果Fig.7 Processing results of inner ring fault signal

針對滾動體故障振動信號，使用同樣方法進行分析。圖8(a)和圖8(b)分別給出了Hilbert包絡后的包絡譜及SWT后的時頻譜，同樣，由于噪聲較大以及滾動體早期微弱故障較難提取等因素，從圖無法找出任何特征。使用VSWT進行處理，同理，從圖8(d)-8(e)中可以看出當K=3，α=550時，峭度最大，分別為245.5和292，確定VMD參數后使用VMD對信號進行分解，經計算第3層峭度值最大，達到了292.0223，故選擇第3模態作為SWT輸入信號，使用SWT進行處理，結果如圖8(c)所示，從圖中可以看出，雖然周圍仍存在噪聲，但在45.07Hz、90.28Hz處都出現了較為清晰的能量集中帶，與滾動體故障特征頻率的1倍頻44.32Hz、2倍頻88.64Hz、較為接近，由此可以判斷滾動軸承滾動體出現了故障，與實驗結果一致。

圖8 滾動體故障信號處理結果Fig.8 Processing results of rolling body fault signal

經實驗驗證，使用常見分析手法，如Hilbert包絡、SWT等提取強噪聲背景下的滾動軸承早期微弱故障振動信號已較為困難。VSWT通過使用峭度指標優化VMD的模態數K及懲罰因子α，使用VMD對信號進行分解后根據最大峭度原則選擇含有有效信息最多的IMF，使用SWT對提取出的最優模態分量進行處理從而提取有效時頻特征。相比 SWT而言，VSWT能夠很好的處理強噪聲背景下的故障信號，優勢顯著。

5 結論

本文研究同步壓縮小波變換(SWT)結合引入最大峭度原則的變分模態分解(VMD)所形成的 VSWT方法能夠將信號分解成K個模態分量，不同的模態分量包含不同的頻率特征，根據最大峭度原則選擇最優IMF作為SWT的輸入信號，即利用CWT對有效IMF進行處理，再使用SST對CWT得到的系數進行壓縮，最后實驗結果表明，該方法能夠有效抑制噪聲，即使在較強的噪聲背景下，也能夠從信號中提取出有效特征頻率，同時VSWT也具有較高的時頻分辨率。