基于局部均值分解和切片雙譜的滾動軸承故障診斷研究

唐貴基，王曉龍

局部均值分解(Local Mean Decomposition，LMD)是在經驗模態分解(Empirical Mode Decomposition，EMD)理論基礎上發展起來的一種較新的時頻分析方法，該方法同樣具有自適應性，它可以將一個復雜的多分量復合信號分解成多個乘積函數(Product Function，PF)分量和一個余量，得到的PF分量能真實的反應原信號的信息［1］。相對于EMD方法來說，LMD方法受端點效應影響較小，并且可以避免EMD算法自身存在的問題，用于機械故障診斷，具有獨特優勢。文獻［2－3］中分別運用LMD方法對齒輪故障和轉子裂紋故障進行診斷，獲得了令人滿意的診斷結果。

雙譜是高階譜中運算最簡單，應用范圍最廣的分析方法，理論上可以完全抑制高斯有色噪聲，提高參數估計的準確性，近些年不少學者對雙譜進行了深入研究，并將雙譜分析方法應用于故障診斷領域。趙慧敏等［4］利用雙譜對柴油發動機曲軸軸承故障信號進行分析，段向陽等［5］則采用切片雙譜對離心泵故障信號進行分析，均成功識別出故障特征頻率。

滾動軸承是旋轉機械系統中的重要組成零部件，其工作狀態正常與否決定了整個系統性能的好壞。當其內圈、外圈或滾動體出現局部損傷或缺陷時，將會產生周期性沖擊振動，輕則使設備產生噪音、振動異常，重則損壞設備，因此對滾動軸承進行故障診斷研究具有重要意義。針對這一問題，本文提出一種結合LMD算法和切片雙譜的診斷新方法，診斷實例證明，本方法能夠滿足實際工程需要。

1 LMD基本原理及算法改進

1.1 LMD算法簡介

對于給定的信號x(t)，LMD算法流程為［6］:

(1)首先確定原給定信號x(t)的所有局部極值點(極大值點和極小值點均包括在內)ni(i=1，2，…，M)，其中M為極值點總數。

(2)計算相鄰極值點均值mi和包絡估計值ai:

(3)將相鄰極值點均值mi和包絡估計值ai分別用直線連接后再利用移動平均法對其進行平滑處理，可得到局部均值函數m11(t)和包絡估計函數a11(t)。

(4)將局部均值函數m11(t)從原給定信號x(t)中分離出來，并利用包絡估計函數a11(t)對其進行歸一化處理，公式表述如下:

如果s11(t)的包絡估計函數a12(t)≠1，則將s11(t)作為原始數據重復步驟(1)和步驟(2)n次，直到它的包絡估計函數a1(n+1)(t)=1，此時，s1n(t)是一個純調頻信號。將迭代過程中生成的所有包絡估計函數相乘便可得到包絡信號a1(t)，公式表示如下:

(5)將包絡信號a1(t)和純調頻信號s1n(t)相乘可得到第一個乘積函數PF1(t)，即:

PF1(t)是一個單分量調幅－調頻信號，包絡信號a1(t)就是它的瞬時幅值，而它的瞬時頻率f1(t)則可通過純調頻信號s1n(t)求取，求取公式為:

(6)從原給定信號x(t)中減去第一個乘積函數PF1(t)可獲得一個新的信號u1(t)，再將u1(t)作為原始數據重復上述步驟(1)—(5)k次，直到uk(t)為一個單調函數或常量則停止迭代。至此，原給定信號x(t)被分解成k個PF分量和1個余量uk(t)，信號x(t)可表示成k個乘積函數分量和1個余量相加的形式，公式表述如下:

1.2 LMD算法改進

原始信號經過LMD算法自適應分解后獲得的PF分量是按頻率由高到低的順序依次排列的。本文認為，LMD分解出的能量過小的乘積函數分量受噪聲干擾嚴重，很難準確從中提取出信號特征信息，對于實際故障診斷基本沒有太大幫助，并且分解過多無意義的PF分量無疑會增加運算量。為減少整個分解過程耗時，加快運算速度，本文提出一個基于能量表述的迭代控制終止條件，具體公式表示如下:

其中:x(i)為原始信號序列，pf(i)為當前分解獲得的乘積函數分量序列(i=1，2，…，m)，m 為序列點總數，θ為評價指標。當評價指標θ＜0.1時，認為當前分解獲得的PF分量能量占原始信號能量小于10%，則可以停止迭代，余下分量均計為殘余分量。

2 峭度篩選準則

峭度是反映信號分布特性的無量綱統計參量，可以描述信號中沖擊成分所占比重的大小，表達式為:

其中:σ和μ分別為信號的標準差和均值，E(t)表示變量t的數學期望。當滾動軸承正常運轉時，振動信號的幅值分布接近正態分布，峭度值約為3，隨著故障的出現和發展，信號中的沖擊成分增多，幅值分布將偏離正態分布，峭度值也隨之變大［7］。

原始信號經過局部均值分解后會得到多個乘積函數分量，各分量都或多或少的攜帶了與軸承故障特征相關的信息，而PF分量的峭度值越大，則說明該分量中的沖擊成分比重越大，越能更好的表現出故障特征。本文認為，峭度指標對于PF分量的篩選具有一定的指導性作用，并在此基礎上提出峭度篩選準則，通過LMD分解結果的峭度值對比，直接篩選出峭度值最大的PF分量，并對其做進一步分析，將更有利于故障特征頻率信息的提取。

3 切片雙譜的定義及性質

零均值、平穩隨機信號x(n)的三階累積量c(τ1，τ2)定義為［8］:

其中，E{·}代表統計均值。

則定義三階累積量c(τ1，τ2)的二維傅里葉變換為雙譜 B(ω1，ω2)，即:

但是雙譜的計算量龐大，不利于過長數據的實時分析，為簡化計算過程，可將雙譜投影到一維頻率空間，通過計算雙譜切片來實現減少運算量的目的。

取 τ1=τ2=τ，則三階累積量的對角切片 c1D(τ)可定義為［9］:

定義三階累積量對角切片c1D(τ)的一維傅里葉變換為切片雙譜B1D(ω)，即:

切片雙譜雖然形式上與功率譜相似，但卻不像功率譜那樣具有明確的物理意義，可將切片雙譜理解為歪度在頻域的分解［10］。當滾動軸承發生局部損傷時，振動信號為多分量調幅調頻信號，并且會偏離高斯分布，解調后的信號中包含了故障特征頻率及其倍頻、轉頻及其倍頻等幅值調制成分，這些調制成分的頻率及相位是互相關聯的，即存在二次相位耦合關系，而切片雙譜由于同時包含了信號的幅值和相位信息，從根本上彌補了功率譜的不足，能夠有效剔除非二次相位耦合諧波項，準確檢測出故障信號的非線性耦合特征，并且切片譜對高斯白噪聲具有良好的抑制作用，將其應用于軸承故障的診斷，分析結果將更加可靠。

4 仿真算例驗證

利用仿真算例來驗證LMD算法的信號分解能力以及切片雙譜的噪聲抑制和非二次相位耦合諧波剔除能力，設模擬信號為:

其中:

w(n)為高斯白噪聲，信噪比SNR為－4 dB;分析點數N=20 000，采樣頻率fs=20 000 Hz，模擬信號及其LMD運算結果如圖1所示。為驗證本文提出的新迭代停止條件在運算量控制上的優勢，將原停止條件分解結果一并列出，進行比較，結果如圖2所示。

圖1 模擬信號及其LMD結果Fig.1 Simulation signal and LMD results

圖2 原停止條件的LMD結果Fig.2 LMD results of original stopping condition

新迭代停止條件LMD運算后將模擬信號分解成2個PF分量和一個余量，而原迭代控制條件LMD運算后共獲得7個PF分量和一個余量，但幾個能量相對較小的低頻分量對于診斷過程并沒有實際意義，反而增加了分解運算量，而新條件獲得的分量個數則相對較少，一定程度上可以減小分解工作量。通過PF1和PF2的幅值譜(如圖3所示)可以看出，原信號中的兩個單分量幅值調制信號被準確分解出來，并且各分量的幅值譜中，調制邊帶隱約可見。圖4和圖5分別為PF1、PF2分量包絡信號的幅值譜和切片譜，下面對幅值譜和切片譜對比結果進行分析。

圖3 PF1和PF2分量的幅值譜Fig.3 Amplitude spectrum of PF1 and PF2 component

圖4 PF1分量包絡信號的幅值譜和切片雙譜Fig.4 Amplitude spectrum and slice bispectrum of PF1 component envelope signal

圖5 PF2分量包絡信號的幅值譜和切片雙譜Fig.5 Amplitude spectrum and slice bispectrum of PF2 component envelope signal

圖4 中，PF1分量包絡信號的切片譜與幅值譜相比，僅在 100 Hz、250 Hz、350 Hz處譜線峰值明顯，而650 Hz處的峰值譜線則消失。圖5中，PF2分量包絡信號的切片譜與幅值譜相比，100 Hz頻率成分消失，15 Hz、40 Hz、55 Hz頻率成分卻依然存在。上述對比說明，PF1 包絡信號中，100 Hz、250 Hz、350 Hz這三個分量存在二次相位耦合，同理，PF2包絡信號中的15 Hz、40 Hz、55 Hz頻率成分也滿足耦合關系，而不滿足非線性耦合關系的650 Hz和100 Hz頻率成分則在切片譜中被剔除，并且與幅值譜對比后發現，切片譜中白噪聲成分幅值相對較小，從而說明切片譜對噪聲具有一定的抑制作用。通過仿真算例，驗證了LMD算法的多分量調制信號分解能力以及切片雙譜的高斯白噪聲抑制特性和非線性耦合特征檢測能力。

5 滾動軸承故障診斷方法流程

基于LMD和切片雙譜的滾動軸承故障診斷方法主要分為以下幾個步驟:

(1)首先對故障信號進行LMD分解，得到一組按照頻率由高到低順序排列的PF分量。

(2)分別計算出各PF分量的峭度值。

(3)篩選出峭度值最大的PF分量，對其包絡信號(瞬時幅值)做進一步的切片雙譜分析。

(4)將軸承故障特征頻率理論計算值與切片雙譜中峰值明顯的譜線進行對比，從而判斷故障類型。

6 滾動軸承故障診斷實例

實測信號為試驗臺模擬的N205滾動軸承內圈、外圈故障信號，實驗臺的結構如圖6所示。采樣頻率fs=3 200 Hz，采樣點數M=4 096，軸承滾動體個數Z=12，轉軸轉速N=1 440 r/min，軸承節圓直徑D=39 mm，滾動體直徑d=7.5 mm，壓力角α=0°，通過計算可得，

圖6 實驗平臺Fig.6 Experimental platform

6.1 滾動軸承外圈故障診斷

外圈故障信號時域波形及幅值譜如圖7所示，圖7(a)中的時域波形雜亂無章，基本無法辨識出信號特征。圖7(b)中三根譜線幅值較明顯，經過分析，僅第二根譜線較接近理論計算的外圈故障特征頻率的二倍頻，而第一、三條譜線均與特征頻率無關，僅靠這一有限信息，無法對故障類型做出精確診斷，為此，利用本文提出的方法對其進行診斷。首先對故障信號進行LMD運算，分解結果如圖8所示。

圖7 外圈故障信號的時域波形及幅值譜Fig.7 Time domain waveform and amplitude spectrum of outer ring fault signal

圖8中新迭代控制條件將故障信號分解成4個PF分量和1個余量，分別計算出故障信號LMD運算獲得的4個PF分量的峭度值，結果如表1所示。為提取外圈故障特征頻率，根據篩選準則，選擇峭度值最大的PF1分量的包絡信號做進一步切片雙譜分析，結果如圖9(a)所示。

圖8 新迭代停止條件的LMD結果Fig.8 LMD results of new iteration stopping condition

表1 各PF分量的峭度值Tab.1 Kurtosis value of PF components

圖9 PF1分量包絡信號的切片雙譜和幅值譜Fig.9 Slice bispectrum and amplitude spectrum of PF1 component envelope signal

當軸承外圈發生故障時，由于外圈固定，因此故障點所處位置及受載荷大小均保持不變，故障信號是以系統固有頻率為高頻載波，外圈故障特征頻率及其倍頻為調制頻率的調制信號［11］，因此包絡信號主要成分即為故障特征頻率及其倍頻。圖9(a)切片雙譜中，117.18 Hz和232.81 Hz兩處譜線峰值突出，分別對應理論計算的軸承外圈故障特征頻率(116.3 Hz)及其二倍頻(232.6 Hz)，綜上所述，可判定滾動軸承外圈存在局部缺陷，理論分析與實際情況相符。

圖9(b)為包絡信號的幅值譜，幅值譜中，116.63 Hz處譜線峰值較高，對應外圈故障特征頻率，但是與切片雙譜對比后可以發現，幅值譜中背景噪聲較嚴重，干擾譜線較多，特征頻率的二倍頻被噪聲淹沒，診斷效果明顯沒有切片譜好。

6.2 滾動軸承內圈故障診斷

內圈故障信號的時域波形及幅值譜如圖10所示，雖然圖10(a)時域波形出現明顯沖擊，但是波形較復雜，僅通過時域波形無法了解故障信息。圖10(b)中，由于故障沖擊的作用，系統的固有頻率已被激起，并且邊頻成分豐富，調制現象明顯，但低頻段譜線雜亂，干擾較多，故障特征頻率譜線并不明顯，利用LMD算法對信號進行分解，結果如圖11所示。

圖10 內圈故障信號的時域波形及幅值譜Fig.10 Time domain waveform and amplitude spectrum of inner ring fault signal

圖11 新迭代停止條件的LMD結果Fig.11 LMD results of new iteration stopping condition

在圖11中，新迭代控制條件將內圈故障信號分解成4個PF分量和1個余量，LMD分解獲得的4個PF分量的峭度值如表2所示，由于PF2的峭度值最大，因此選擇PF2分量的包絡信號做切片雙譜分析，結果如圖12(a)所示。

表2 各PF分量的峭度值Tab.2 Kurtosis value of PF components

圖12 PF2分量包絡信號的切片雙譜和幅值譜Fig.12 Slice bispectrum and amplitude spectrum of PF2 component envelope signal

軸承內圈發生故障時，由于內圈隨轉軸一起旋轉，因此故障點與滾動體接觸位置及受載大小均呈周期性變化，故障信號是以系統固有頻率為載波，轉頻及內圈故障特征頻率為調制頻率的調制信號，并且轉頻對特征頻率也具有調制作用［11］，因此包絡信號中主要成分為轉頻、內圈故障特征頻率及其邊帶頻率。圖12(a)切片雙譜中有四根譜線幅值明顯，其中 23.44 Hz和46.88 Hz分別對應轉頻(24 Hz)及其二倍頻(48 Hz)，171.87 Hz和148.43 Hz則分別對應內圈故障特征頻率理論計算值(171.7 Hz)及其轉頻調制邊帶(147.7 Hz)，結合切片譜檢測故障耦合特征的性質可知，轉頻與內圈故障特征頻率通過二次相位耦合產生了特征頻率左側的邊頻成分，綜合分析后可斷定軸承內圈存在局部損傷。通過對圖12(b)幅值譜的分析，同樣能識別出故障特征頻率(171.82 Hz)及其調制邊帶(148.39 Hz)譜線，但是切片譜對噪聲的抑制效果更好，切片譜中的故障特征頻率及其相關譜線與幅值譜相比也更加明顯，內圈故障診斷實例對比再次驗證了切片雙譜在噪聲抑制和故障非線性耦合特征檢測上的優勢。

7 結論

LMD算法適合處理多分量調幅—調頻信號，而切片雙譜作為高階譜的簡化算法，又是分析非平穩、非高斯信號的強有力工具，本文將局部均值分解方法與切片雙譜分析方法相結合，并以峭度指標做為PF分量的篩選準則，提出了基于LMD和切片雙譜的滾動軸承故障診斷方法。運用該方法對故障信號進行分析，比單獨使用LMD診斷方法效果更明顯，并通過滾動軸承內圈、外圈故障診斷實例，驗證了新方法的可行性和準確性。

［1］程軍圣，楊 怡，張 亢，等.基于局部均值分解的循環頻率和能量譜在齒輪故障診斷中的應用［J］.振動工程學報，2011，24(1):78 －83.CHENG Jun-sheng，YANG Yi，ZHANG Kang，et al.Application of cycle frequency and energy spectrum based on local mean decomposition to gear fault diagnosis［J］.Journal of Vibration Engineering，2011，24(1):78－83.

［2］何 田，林意洲，郜普剛，等.局部均值分解在齒輪故障診斷中的應用研究［J］.振動與沖擊，2011，30(6):196－201.HE Tian，LIN Yi-zhou，GAO Pu-gang，et al.Application of local mean decomposition in gear fault diagnosis［J］.Journal of Vibration and Shock，2011，30(6):196 －201.

［3］任達千，楊世錫，吳昭同，等.LMD時頻分析方法的端點效應在旋轉機械故障診斷中的影響［J］.中國機械工程，2012，23(8):951 －956.REN Da-qian，YANG Shi-xi，WU Zhao-tong，et al.Research on end effect of LMD based time-frequency analysis in rotating machinery fault diagnosis［J］.China Mechanical Engineering，2012，23(8):951－956.

［4］趙慧敏，夏超英，肖云魁，等.柴油發動機曲軸軸承振動信號的雙譜分析［J］.振動、測試與診斷，2009，29(1):14－18.ZHAO Hui-min，XIA Chao-ying，XIAO Yun-kui，et al.Bispectrum analysis for vibration data of crankshaft bearing in diesel engine［J］.Journal of Vibration，Measurement ＆Diagnosis，2009，29(1):14 －18.

［5］段向陽，王永生，蘇永生，等.切片雙譜分析在離心泵故障診斷中的應用［J］.振動、測試與診斷，2010，30(5):581－584.DUAN Xiang-yang，WANG Yong-sheng，SU Yong-sheng，et al.Application of slice bispectrum analysis to fault diagnosis of centrifugal pump［J］.Journal of Vibration，Measurement ＆Diagnosis，2010，30(5):581 －584.

［6］Smith J S.The local mean decomposition and its application to EEG perception data［J］.Journal of the Royal Society Interface，2005，2(5):443 －454.

［7］蘇文勝，王奉濤，張志新，等.EMD降噪和譜峭度法在滾動軸承早期故障診斷中的應用［J］.振動與沖擊，2010，29(3):18－21.SU Wen-sheng，WANG Feng-tao，ZHANG Zhi-xin，et al.Application of EMD denoising and spectral kurtosis inearly fault diagnosis of rolling element bearings［J］.Journal of Vibration and Shock，2010，29(3):18 －21.

［8］ Nikias C L，Raghuveer M R.Bispectrum estimation:A digital signal processing framework［J］.Proceedings of IEEE，1987，75(7):869－891.

［9］楊江天，陳家驥，曾子平.雙譜分析及其在機械診斷中的應用［J］.中國機械工程，2000，11(4):424 －426.YANG Jiang-tian，CHEN Jia-ji，ZENG Zi-ping.Bispectral analysis and its application in machinery diagnosis［J］.China Mechanical Engineering，2000，11(4):424－426.

［10］ Rivola A，White P R.Detection of rolling element bement bearing damage by statistical vibration analysis［J］.Journal of Sound and Vibration，1998，216(6):889 －910.

［11］胡愛軍，馬萬里，唐貴基.基于集成經驗模態分解和峭度準則的滾動軸承故障特征提取方法［J］.中國電機工程學報，2012，32(11):106 －111.HU Ai-jun，MA Wan-li，TANG Gui-ji.Rolling bearing fault feature extraction method based on ensemble empirical mode decomposition and kurtosis criterion［J］.Proceedings of the CSEE，2012，32(11):106－111.