基于VMD-SWT滾動(dòng)軸承故障診斷方法研究＊

周進(jìn)群，劉義亞

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

前言

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

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

峭度是描述隨機(jī)變量分布特性的數(shù)值統(tǒng)計(jì)量，反應(yīng)信號(hào)沖擊成分的大小。本文將最大峭度指標(biāo)引入 VMD，選擇出最優(yōu)分解模態(tài)數(shù)及懲罰因子 α，再利用峭度指標(biāo)選擇有效VMD分解分量從而從噪聲信號(hào)中提取出有效信息，最后使用 SWT提取信號(hào)中的時(shí)頻特征，從而有效判斷滾動(dòng)軸承的運(yùn)行狀態(tài)。相比于EMD等方法來(lái)說(shuō)，VMD具有嚴(yán)謹(jǐn)?shù)睦碚撝?，且是根?jù)中心頻率選取模態(tài)數(shù)，提高了計(jì)算效率，能夠?qū)⒃夹盘?hào)進(jìn)行分解從而獲得K個(gè)含有不同頻率成分的模態(tài)分量；SWT方法雖然能夠提高時(shí)頻譜的可讀性，但當(dāng)噪聲較大時(shí)仍然無(wú)提取出有效特征，而工業(yè)現(xiàn)場(chǎng)會(huì)由于環(huán)境、滾動(dòng)軸承工作機(jī)理等會(huì)造成故障振動(dòng)信號(hào)中包含很大的噪聲，導(dǎo)致信噪比極低，將VMD引入SWT中，能夠?qū)⒃夹盘?hào)進(jìn)行分解，從而提取出含有有效頻率成分的模態(tài)分量，提高信噪比，再使用 SWT進(jìn)行分析，從而有效提取信號(hào)的時(shí)頻特征。

1 變分模態(tài)分解VMD

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

式中，φk(t)非遞減，；包絡(luò)線非負(fù)，Ak(t)≥0；瞬時(shí)頻率；Ak(t)及ωk(t)相對(duì)于φk(t)來(lái)說(shuō)是緩變的。因此，在區(qū)間內(nèi)，uk(t)可看作是幅值為 Ak(t)且頻率為ωk(t)的諧波信號(hào)。

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

引入懲罰因子α構(gòu)造增廣Lagrange函數(shù)從而求取上述變分模型的最優(yōu)解，Lagrange構(gòu)造函數(shù)如下：

式中：α為懲罰因子；λ為L(zhǎng)agrange乘子。

將上述Lagrange函數(shù)從時(shí)域變換到頻域，求解相應(yīng)的極值，可得到uk和ωk的頻域表達(dá)式。

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

由上述步驟即可將信號(hào)分解成K個(gè)模態(tài)分量。

2 同步壓縮小波變換

峭度(Kurtosis)是反映振動(dòng)信號(hào)分布特性的數(shù)值統(tǒng)計(jì)量，是無(wú)量綱參數(shù)，對(duì)信號(hào)的沖擊十分敏感，適用于表面損傷類故障，尤其是早期故障的診斷[12]：

在軸承正常運(yùn)轉(zhuǎn)狀態(tài)下，其振動(dòng)信號(hào)的幅值分布接近正態(tài)分布，當(dāng)故障開(kāi)始出現(xiàn)時(shí)，信號(hào)幅值的分布會(huì)偏離正態(tài)分布，峭度值也隨之增大。峭度指標(biāo)的絕對(duì)值越大，說(shuō)明軸承偏離其正常狀態(tài)，故障越嚴(yán)重。

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

利用CWT對(duì)信號(hào)x(t)進(jìn)行處理，可定義為：

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

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

其中

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

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

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

（3）同理計(jì)算出最優(yōu)α值(α步長(zhǎng)為10)；

（4）將計(jì)算出的最優(yōu)K和α帶入VMD中，使用VMD將信號(hào)進(jìn)行分解；

（5）計(jì)算各模態(tài)峭度值，選出含有最大峭度值的模態(tài)分量；

（6）使用SWT對(duì)選出的模態(tài)分量進(jìn)行分析。

本文提出使用最大峭度法提取出VMD的最優(yōu)參數(shù)，將VMD引入SWT算法中構(gòu)成一種新的分析方法(VSWT)，首先使用峭度指標(biāo)判斷出VMD最優(yōu)參數(shù)，再將提取出的參數(shù)帶入VMD中對(duì)原始信號(hào)進(jìn)行分解，最后領(lǐng)用峭度指標(biāo)提取有效模態(tài)分量后使用 SWT進(jìn)行處理分析，提取出信號(hào)的時(shí)頻特征，從而對(duì)滾動(dòng)軸承運(yùn)轉(zhuǎn)狀態(tài)進(jìn)行判斷。

3 仿真分析

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

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

圖1 VMD分解Fig.1 Variational Mode Decomposition

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

圖2 SWT處理結(jié)果Fig.2 Processing results by SWT

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

圖3 原信號(hào)與重構(gòu)信號(hào)對(duì)比Fig.3 The original signal and the signals reconstructed by IVSWT

綜上所述，可以判斷出VSWT具有良好的抗噪能力以及較高的信號(hào)提取精度，能夠在強(qiáng)噪聲背景下有效提取信號(hào)時(shí)頻特征，提高時(shí)頻可讀性，同時(shí)能夠有效的重構(gòu)信號(hào)。

4 實(shí)驗(yàn)和驗(yàn)證

4.1 實(shí)驗(yàn)說(shuō)明

為驗(yàn)證改進(jìn) SWT的工程實(shí)用性，本文以深南電路鉆孔工序的鉆床為實(shí)驗(yàn)平臺(tái)，通過(guò)更換故障軸承采集滾動(dòng)軸承正常振動(dòng)信號(hào)、外圈故障信號(hào)、內(nèi)圈故障信號(hào)、滾動(dòng)體故障信號(hào)，分別應(yīng)用 Hilbert包絡(luò)、SWT、VSWT分析故障信號(hào)，通過(guò)提取的時(shí)頻特征判斷滾動(dòng)軸承的故障類型。

通常滾動(dòng)軸承常見(jiàn)的故障類型為外圈故障、內(nèi)圈故障以及滾動(dòng)體故障。其故障特征頻率計(jì)算公式為：

內(nèi)圈故障特征頻率

外圈故障特征頻率

滾動(dòng)體故障特征頻率

式中，fr為轉(zhuǎn)頻；d為滾珠直徑；D為節(jié)圓直接； 為接觸角；Z為滾珠數(shù)。

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

表1 軸承參數(shù)Tab.1 Parameter of bearing

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

4.2 信號(hào)處理及分析

對(duì)采集的故障信號(hào)進(jìn)行處理，圖5給出了滾動(dòng)軸承正常振動(dòng)信號(hào)、外圈故障振動(dòng)信號(hào)、內(nèi)圈故障振動(dòng)信號(hào)以及滾動(dòng)體故障振動(dòng)信號(hào)的時(shí)域波形圖。

圖5 軸承振動(dòng)信號(hào)Fig.5 The bearing vibration signal

為驗(yàn)證改進(jìn) SWT方法的有效性與實(shí)用性，本文使用不同方法處理滾動(dòng)軸承故障數(shù)據(jù)，并對(duì)處理結(jié)果進(jìn)行比較。

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

圖6 外圈故障信號(hào)處理結(jié)果Fig.6 Processing results of outer ring fault signal

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

圖7 內(nèi)圈故障信號(hào)處理結(jié)果Fig.7 Processing results of inner ring fault signal

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

圖8 滾動(dòng)體故障信號(hào)處理結(jié)果Fig.8 Processing results of rolling body fault signal

經(jīng)實(shí)驗(yàn)驗(yàn)證，使用常見(jiàn)分析手法，如Hilbert包絡(luò)、SWT等提取強(qiáng)噪聲背景下的滾動(dòng)軸承早期微弱故障振動(dòng)信號(hào)已較為困難。VSWT通過(guò)使用峭度指標(biāo)優(yōu)化VMD的模態(tài)數(shù)K及懲罰因子α，使用VMD對(duì)信號(hào)進(jìn)行分解后根據(jù)最大峭度原則選擇含有有效信息最多的IMF，使用SWT對(duì)提取出的最優(yōu)模態(tài)分量進(jìn)行處理從而提取有效時(shí)頻特征。相比 SWT而言，VSWT能夠很好的處理強(qiáng)噪聲背景下的故障信號(hào)，優(yōu)勢(shì)顯著。

5 結(jié)論

本文研究同步壓縮小波變換(SWT)結(jié)合引入最大峭度原則的變分模態(tài)分解(VMD)所形成的 VSWT方法能夠?qū)⑿盘?hào)分解成K個(gè)模態(tài)分量，不同的模態(tài)分量包含不同的頻率特征，根據(jù)最大峭度原則選擇最優(yōu)IMF作為SWT的輸入信號(hào)，即利用CWT對(duì)有效IMF進(jìn)行處理，再使用SST對(duì)CWT得到的系數(shù)進(jìn)行壓縮，最后實(shí)驗(yàn)結(jié)果表明，該方法能夠有效抑制噪聲，即使在較強(qiáng)的噪聲背景下，也能夠從信號(hào)中提取出有效特征頻率，同時(shí)VSWT也具有較高的時(shí)頻分辨率。