基于DS-VMD 及相關峭度的滾動軸承故障診斷?

石文杰， 黃 鑫， 溫廣瑞，2， 張志芬

（1.西安交通大學機械工程學院 西安，710049） （2.新疆大學機械工程學院 烏魯木齊，830047）

引 言

滾動軸承是旋轉(zhuǎn)機械中應用最為廣泛的零部件之一，復雜的工作環(huán)境和長時間的交變載荷極易使軸承產(chǎn)生故障。對軸承運行狀態(tài)進行實時監(jiān)測是保證高端機械裝備正常運行的必要環(huán)節(jié)，也是提高機械設備運行安全性、穩(wěn)定性和可靠性的重要手段。因此，對于滾動軸承的狀態(tài)監(jiān)測和故障診斷具有重要意義。

常用的滾動軸承故障診斷方法有快速峭度圖、最小熵解卷積和包絡譜分析等。鐘先友等［1］將迭代濾和快速峭度圖相結合用于滾動軸承的微弱故障特征提取。祝小彥等［2］提出了多點最優(yōu)調(diào)整的最小熵解卷積與Teager 能量算子相結合的滾動軸承故障診斷方法，并對信號進行包絡分析實現(xiàn)軸承的故障診斷。陳祥龍等［3］提出平方包絡譜相關峭度新指標，結合Morlet 小波濾波和粒子群優(yōu)化算法，實現(xiàn)滾動軸承故障診斷與識別。

近年來，自適應信號分析方法在軸承故障診斷領域得到了廣泛應用。徐可等［4］提出一種自適應波形匹配的延拓方法，對EMD 存在的端點效應進行改進，并結合優(yōu)化的支持向量機實現(xiàn)軸承故障診斷。栗蘊琦等［5］采用經(jīng)驗小波變換對信號進行自適應分解，對各分量進行奇異值分解重構，有效提取出軸承故障特征。王普等［6］將信號進行局部均值分解，通過相關系數(shù)選擇合適的分量并獲得新的時間序列，構建多尺度熵能量進行故障特征提取。Dragomiretskiy 等［7］提出的VMD 方法是一種常用的自適應信號分析方法。Tang 等［8］結合VMD 和盲源分離方法，通過獨立分量分析對軸承復合故障進行診斷。Jiang 等［9］提出了一種能量波動譜來突出潛在模態(tài)初始中心頻率，并構造了一種基于中心頻率引導的VMD 優(yōu)化策略，通過自適應細化平衡參數(shù)來提取最優(yōu)模式以進行旋轉(zhuǎn)機械故障診斷。Yan 等［10］提出了一種改進的尺度空間VMD 方法，通過應用尺度空間表示方法獲取待分析信號的先驗知識以確定VMD 所需的參數(shù)。相關研究表明，VMD 方法的分解效果受分解模態(tài)數(shù)和懲罰因子參數(shù)的影響較大，分解參數(shù)的確定通常依據(jù)經(jīng)驗或者待處理信號的先驗知識，而在工程實際中，振動信號先驗知識不易獲得，VMD 參數(shù)選取不準確會導致信號處理結果出現(xiàn)偏差。

針對上述問題，筆者提出一種基于微分搜索優(yōu)化的VMD 參數(shù)自適應尋優(yōu)算法，使用相關峭度指標對分解后的各本征模態(tài)函數(shù)進行評價及加權重構，并對重構信號進行包絡譜分析以提取軸承故障特征。仿真和實驗結果表明，該方法可以有效提取滾動軸承故障特征。

1 基本原理

1.1 變分模態(tài)分解

VMD 作為一種非遞歸信號分解方法，其本質(zhì)是一種具有窄帶特性的自適應維納濾波器組。VMD 將實信號分解為指定數(shù)量的有限帶寬的IMF，并使每個分量的估計帶寬之和最小，其約束變分問題模型為

其中：uk為分解后的各個IMF；ωk為每個IMF 的中心頻率。

引入二次懲罰因子α和拉格朗日乘法算子λ，使約束變分問題轉(zhuǎn)變?yōu)榉羌s束變分問題

其中：f為原始信號。

在VMD 方法中，分解模態(tài)數(shù)K和二次懲罰因子α對分解結果有較大影響。K值決定了信號分解后的模態(tài)數(shù)，α值影響分解后各IMF 的頻帶寬度。α值越大，則分解后得到的IMF 分量所處的頻帶分布范圍越窄。二者的組合對VMD 方法的影響如表1所示。可見，K值的選擇對于能否準確確定故障特征頻率所處的共振頻帶非常重要，而α值的正確選擇保證了VMD 算法進行信號重構時的精度［11］。分解結果的偏差最終會導致振動信號故障信息的缺失。

表1 K 值和α 值對VMD 方法的影響Tab.1 Influence of K and α on VMD

1.2 微分搜索

微分搜索算法由Civicioglu［12］首次提出。Wen等［13］證明該算法具有靈活性和穩(wěn)定性，且對高維多模態(tài)函數(shù)的求解更有效。DS 算法模擬了自然界中動物種群因自然資源容量的變化而發(fā)生的季節(jié)性遷徙行為。DS 算法中，問題的隨機解組成一個種群，種群會逐步遷徙到問題的最優(yōu)值。在遷徙中種群會比較暫時停留位置是否優(yōu)于當前位置，如果優(yōu)于當前位置，種群的成員立刻定居在這個位置臨時停留并且在下一次迭代中從這個位置繼續(xù)遷徙，否則將本次迭代的暫時停留位置丟棄，保持原位置不變。

在DS 算法中，參與遷徙的所有個體（xi，i=1，2，…，N）組成一個種群（Sg，g=1，2，…，m），其中每個個體包含的元素（xi，j，j=1，2，…，D）等于問題的維數(shù)。其中：N為個體總數(shù)；m為總的迭代次數(shù)；D為問題的維數(shù)。

個體中元素的初始位置定義為

其中：xi，j為種群的第i個個體的第j維元素；r為0～1之間滿足均勻分布的隨機數(shù)；Uj為搜索空間第j維的上限；Lj為搜索空間第j維的下限。

在DS 算法中，尋找一個暫時停留位置的機制可以描述為一個類布朗隨機行走運動模型。種群向目標移動和個體元素位置的變化大小受比例值的控制，個體元素參加暫時停留位置的搜索過程由一個隨機過程決定。該隨機過程表示為

其中：d表示種群移動目標，為一個N×D維的矩陣；比例值sscale為一個滿足伽馬分布的隨機數(shù)，通過該隨機過程產(chǎn)生暫時停留位置；Ssite為一個N×D維的矩陣；S為種群當前所在的位置，為一個N×D的矩陣；mmap為個體元素的選擇策略，為一個N×D的［0，1］隨機整數(shù)矩陣；0 表示元素沒被選中，1 表示元素被選中；N為個體數(shù)量；D為問題維數(shù)。

通過式（6）可以看出，移動目標并不等于暫時停留位置，種群中的個體會選擇一些元素按一定的比例向移動目標方向產(chǎn)生暫時停留位置，選擇策略和比例值都是按一定規(guī)律隨機產(chǎn)生。如果一個暫時停留位置元素因為某些原因超出尋優(yōu)空間的界限，則該元素隨機轉(zhuǎn)移到區(qū)域中另一個位置。

1.3 相關峭度

相關峭度是一種基于峭度衍生的統(tǒng)計量。峭度是信號尖峰的量度，反應隨機變量的分布特性，是檢測旋轉(zhuǎn)部件故障沖擊性的重要指標［14］。在軸承正常運行過程中，其振動信號的幅值分布接近于正態(tài)分布，因此峭度值約等于3。若軸承發(fā)生故障，其振動信號幅值分布將偏離正態(tài)分布，峭度值也會隨之發(fā)生變化。因此，峭度值可以用來區(qū)分滾動軸承的正常狀態(tài)與故障狀態(tài)［15］，其計算公式為

其中：x為離散信號為信號均值；N為采樣長度；σ為信號的標準差。

由于峭度指標未考慮軸承故障沖擊的周期性，因此基于峭度對軸承運行狀態(tài)進行判斷的方法易受到異常值的影響，魯棒性較差。相關峭度充分考慮滾動軸承故障信號的周期性，能更加穩(wěn)定、準確地提取滾動軸承故障特征［16］。相關峭度的表達式為

其中：T為解卷積周期，其設定與故障沖擊周期相關，代表信號相鄰兩個沖擊脈沖間的數(shù)據(jù)點數(shù)；M為移位數(shù)，其影響信號解卷積處理后所提取沖擊脈沖的個數(shù)。

需要注意的是，當M=1，T=0 時，相關峭度等同于峭度。與峭度相比，相關峭度更適合于衡量滾動軸承故障的周期性沖擊成分。相關峭度值越大，代表所感興趣的周期性沖擊成分在信號中所占的比重越大，軸承的故障特征越明顯。

2 基于DS-VMD 及相關峭度的信號重構

針對滾動軸承故障振動沖擊信號的周期性特性，筆者提出了一種基于DS-VMD 及相關峭度的信號重構方法，其處理流程如圖1 所示。具體的處理步驟如下。

1）設定VMD 中參數(shù)α和K的搜索范圍。

圖1 所提出方法的處理流程Fig.1 Flowchart for the proposed method

2）設定DS 算法的適應度函數(shù)。利用DS 算法搜尋最優(yōu)的分解參數(shù)時，需要設定一個適應度函數(shù)。Shannon 信息熵作為一種評價信號稀疏性的標準，其值的大小反映了信號的不確定程度，值越大，則信號的不確定性越大［17］。因此，將信號分解后的包絡信號序列pj的熵值作為適應度值。對于給定信號s(t)，其包絡熵值Ep為

包絡熵值越小，說明信號稀疏度越高，包含的周期性成分越多。因此，令分解后IMF 的最大包絡熵值最小，可以保證各分量都具有較好的稀疏性，即DS 算法的適應度函數(shù)為

其中：Epk為第k個IMF的包絡熵值。

將局部min{Ep}作為尋優(yōu)過程中的適應度值，最終目標是求得全局min{Ep}。

3）通過DS 算法得到VMD 參數(shù)最優(yōu)組合。

4）根據(jù)步驟3 得到的參數(shù)結果對時域振動信號進行VMD。

5）根據(jù)軸承故障特征頻率，計算各IMF 的相關峭度。

6）根據(jù)相關峭度值的大小，賦予各IMF 相應的權重并重新構造時域振動信號。假設通過DS 算法確定的VMD 分解模態(tài)數(shù)為K，各層信號分別為x1，x2，…，xk，各模態(tài)的相關峭度值分別為ck1，ck2，…，ckk，各IMF 對應的重構系數(shù)分別為c1，c2，…，ck，則重新構造的信號xnew為

7）對新構造的信號做包絡譜分析，提取相應的軸承故障特征頻率。

3 數(shù)據(jù)驗證分析

3.1 仿真信號分析

構建模擬振動信號來驗證所提出方法的有效性。模擬信號x(t)由齒輪嚙合分量xg(t)、軸承故障分量xb(t)和噪聲分量ε組成

每部分信號的構造方法為

當轉(zhuǎn)頻fr設定為80 Hz 時，故障特征頻率fp設定為50 Hz，軸承部件的共振頻率fn設定為8 kHz，阻尼比ζ設定為0.01，齒數(shù)Z設定為30。

添加高斯白噪聲使仿真信號的信噪比為?5 dB。信號的采樣頻率為20480 Hz，所分析的數(shù)據(jù)長度為20480 個點。仿真信號的時域波形如圖2 所示，發(fā)現(xiàn)軸承故障的周期性沖擊特征被湮沒在噪聲中。

圖3 為仿真信號快速傅里葉變換（fast Fourier transform，簡稱FFT）頻譜及包絡譜。從圖3（a）可看出，仿真信號x(t)的傅里葉譜、齒輪振動分量和軸承振動分量可以很容易地在2.4 kHz，4.8 kHz 和8 kHz的頻譜范圍內(nèi)找到。圖3（b）為仿真信號包絡譜，從中可以清晰地觀察到80 Hz 的轉(zhuǎn)頻。雖然圖中50 Hz的故障特征頻率及其倍頻也可以被反映出來，但其幅值分量明顯小于轉(zhuǎn)頻分量，不易直接從包絡譜圖中讀取。

圖2 仿真信號的時域波形Fig.2 The bearing fault simulated signal

圖3 仿真信號FFT 頻譜及包絡譜Fig.3 FFT spectrum and envelope spectrum of simulated signal

故障特征比（fault feature ratio，簡稱FFR）可以作為定量指標來對軸承某一感興趣的故障頻率成分進行選擇。FFR 值越大，說明所選擇的頻率分量中所包含的周期性沖擊信息越豐富。故障特征比Rf定義［18］為

其中：f為所感興趣的故障特征頻率；S為時域信號包絡譜的幅值總和；S(kf)為故障特征頻率各倍頻對應的包絡頻譜幅值。

圖4 為仿真信號DS-VMD 分解結果。將仿真信號x(t)輸入到所提出方法的模型中，通過DS 算法確定VMD 的參數(shù)組合。設定二次懲罰因子α和分解模態(tài)數(shù)K的搜索范圍分別為［1000，5000］和［2，10］。 DS 算 法 得 到 的 全 局 最 優(yōu) 解 為α=2653，K=5。

根據(jù)DS 得到的參數(shù)對信號進行VMD，利用相關峭度指標對各IMF 進行加權重構，處理后信號包絡譜如圖5 所示，從圖中可以清晰地看到軸承故障特征頻率及其各倍頻。與原始振動信號相比，重構信號的FFR 值由0.065 提升至0.164，提升率達到152.3%，更加說明了所提出方法對于提取滾動軸承故障特征的有效性。

圖4 仿真信號DS-VMD 分解結果Fig.4 DS-VMD result of simulated signal

圖5 處理后信號包絡譜Fig.5 Envelope spectrum after processing

為了驗證所提方法的優(yōu)越性，分別采用EMD方法和傳統(tǒng)VMD 方法進行對比。VMD 中模態(tài)數(shù)K取5，α取500。圖6 為仿真信號EMD 結果，仿真信號被分解為13 個有效模態(tài)。

EMD 與VMD 分解后通過相關峭度指標加權的分析結果如圖7 所示?？梢钥闯?，2 種方法得到的結果中，轉(zhuǎn)頻幅值均較大，且在故障特征頻率及其倍頻處噪聲較多，均無法明顯提取軸承故障特征。

圖6 仿真信號EMD 結果Fig.6 Simulated signal EMD result

表2 為3 種方法得到的包絡譜FFR 值?？梢?，所提出方法的FFR 值最大，說明該方法相較于傳統(tǒng)方法具有明顯的優(yōu)越性。

圖7 EMD 與VMD 分析結果Fig.7 Result of EMD and VMD

表2 不同方法得到的包絡譜FFR 值Tab.2 FFR of envelope spectrum obtained by different methods

3.2 實驗信號分析

為了進一步驗證所提方法的有效性，在SpectraQuest 公司設計的風力渦輪機傳動系診斷模擬器上進行了實驗測試，如圖8 所示。該平臺由2.24 kW驅(qū)動電機、兩級變速箱、行星齒輪箱和磁粉制動器組成。將有外圈缺陷的滾動軸承安裝在齒輪箱的輸入軸上，軸承的外圈故障特征頻率是輸入軸轉(zhuǎn)動頻率的3.048 倍。

直流電機的轉(zhuǎn)速設定為1500 r/min，輸入軸的轉(zhuǎn)速為25 Hz，故障軸承的外圈缺陷特征頻率為76.2 Hz，采樣頻率為20480 Hz，采樣數(shù)為20480。圖9 為振動信號時域波形及包絡譜。從圖9（a）可以看出，軸承故障的周期性沖擊被軸的振動和齒輪嚙合所湮沒。圖9（b）所示的包絡譜也無法得出相應的故障特征頻率信息。

圖8 風力渦輪機傳動系診斷模擬器Fig.8 Wind turbine drive train diagnostic simulator

圖9 振動信號時域波形及包絡譜Fig.9 Time domain waveform and envelope spectrum of vibration signal

將實驗信號輸入到所提出方法的模型中，通過DS 算法得到VMD 的參數(shù)α=2926，K=3，其分解結果如圖10 所示。

重構信號時域波形及包絡譜如圖11 所示?？梢钥闯?，相比于原始信號包絡譜，重構信號包絡譜中軸承故障特征頻率及其各倍頻均可以更清晰地被觀察到，且故障特征頻率的幅值大于轉(zhuǎn)頻幅值。與原始信號相比，重構信號的FFR 值由0.067 提升至0.128，提升率高達91%。

圖10 實驗信號DS-VMD 分解結果Fig.10 DS-VMD result of vibration signal

圖11 重構信號時域波形及包絡譜Fig.11 Time domain waveform and envelope spectrum of reconstructed signal

從圖12 所示的DS 算法迭代過程可以看出，在運算進行到第8 次時，算法已經(jīng)收斂。這說明將DS算法運用在軸承實際振動信號中，可以較快地收斂到較優(yōu)的結果，且算法效率較高。

圖12 DS 算法迭代過程Fig.12 Iteration process of DS algorithm

為了說明所提方法的優(yōu)越性，選擇EMD 方法及人為確定參數(shù)的傳統(tǒng)VMD 方法進行對比。VMD 方法中K在未知的情況下取5，α取500。2 種方法分析得到的包絡譜如圖13 所示。圖13 所示的包絡譜中，低頻處譜線干擾噪聲均較多，無法準確識別軸承故障特征頻率。

表3 不同方法所得包絡譜的FFR 值。所提方法的FFR 值最大，說明該方法相較于傳統(tǒng)方法更優(yōu)越。

圖13 兩種不同方法得到的包絡譜Fig.13 Envelope spectrum obtained by two different methods

為了進一步說明所提出方法的優(yōu)勢，使用快速峭度圖對實驗信號進行分析?？焖偾投葓D是軸承故障診斷中的常用方法，文獻［19-20］提出譜峭度概念，正式給出了其數(shù)學定義，并首次將其應用在滾動軸承的故障診斷中。Antoni［21］簡化了譜峭度的計算，給出了快速峭度圖方法，并給出沖擊分量所處的頻帶位置。筆者利用快速峭度圖法對仿真信號進行分析，以表明所提方法的優(yōu)越性。

表3 不同方法所得包絡譜的FFR 值Tab.3 The FFR of envelope spectrum obtained by different methods

圖14 為振動信號快速峭度圖及濾波結果。從圖14（b）可以看出，濾波后的時域信號包絡在0.1 s處有一明顯的瞬態(tài)沖擊，但該沖擊并不具有周期性。這說明該時刻的沖擊成分對應的峭度值最大，導致最終結果無法反映出軸承故障特征。筆者提出的方法由于考慮到滾動軸承故障沖擊的周期性，可實現(xiàn)軸承故障特征的準確提取。

圖14 振動信號快速峭度圖及濾波結果Fig.14 Fast kurtogram of vibration signal and results

4 結 論

1）采用DS 算法自適應確定VMD 參數(shù)最優(yōu)組合，避免了在先驗知識不足的情況下設置參數(shù)所導致的分解結果不理想的情況。仿真和實驗信號均證明了該方法的有效性。

2）采用相關峭度指標對VMD 分解后的各IMF進行加權重構。相比于原始振動信號，重構信號同時包含轉(zhuǎn)頻信息及故障特征信息。與傳統(tǒng)濾波方法相比，所提出方法可以再識別軸承故障特征方面表現(xiàn)更好。