基于平方包絡譜負熵準則的軸承早期復合故障特征提取方法

陳 鵬， 趙小強,2,3

(1.蘭州理工大學 電氣工程與信息工程學院, 蘭州 730050；2.甘肅省工業過程先進控制重點實驗室，蘭州 730050；3.蘭州理工大學 國家級電氣與控制工程實驗教學中心，蘭州 730050)

在石油化工、有色冶金等領域的生產裝置中，鼓風機、壓縮機、風機和泵等旋轉設備是整個生產裝置的“心臟”，也被稱為“關鍵機組”。這些關鍵機組具有轉速高、功率大、無備機和維修周期長的特點，如一旦發生故障而停機檢修，則會對整個企業生產造成牽一發而動全身的重大影響，而軸承是關鍵機組在長期運行過程中最易發生故障的零部件，如果不對其早期故障進行及時監測和診斷，則會引發嚴重故障而導致停機檢查和維修，從而造成嚴重經濟損失。而在早期故障中，早期復合故障由于具有多信號耦合、信號微弱和強噪聲干擾的特點，相對早期單故障和成熟復合故障而言更加難以辨識[1]。因此，研究能夠提取軸承早期復合故障特征的有效方法一直是故障診斷領域的難點[2]。

目前，以時頻分析為主的信號分解方法在軸承的降噪和故障診斷中得到了廣泛應用[3]。如Huang等[4]提出了經驗模態分解(empirical mode decomposition,EMD)方法，但是EMD方法存在對噪聲敏感、模態混疊和端點效應的問題，因此有學者們針對EMD的問題，提出了集總經驗模態分解(ensemble empirical mode decomposition，EEMD)[5]、互補集合經驗模態分解(coplementary ensemble empirical mode decomposition,CEEMD)[6]。同時，學者們提出了其他一些自適應分解方法，如優化集合局部均值分解(ensemble local mean decomposition,ELMD)[7]、局部特征尺度分解(local characteristic-scale decomposition,LCD)[8]和固有時間尺度分解(intrinsic time scale decomposition,ITD)[9]等。但這些方法都是基于遞歸分解的原理，不可避免地存在一定的模態混疊和端點效應[10]。而Dragomiretskiy等[11]于2014年提出了一種新的信號自適應分解方法—變分模態分解(variational mode decomposition，VMD)，近年來被廣泛應用于旋轉機械的振動信號處理。研究結果表明，VMD可以在一定程度上改善模態混疊和端點效應[12]，但是該方法在早期復合故障診斷領域的效果欠佳[13]。對此，Apostolidis等[14]于2016年在群濾波(swarm filter, SWF)的基礎上，提出了一種新的信號智能分解方法—群分解(swarm decomposition, SWD)，它通過迭代濾波，可將多個頻率相近的復合信號自適應地分解為若干單一模態的振蕩分量，具有較高的頻率區分能力，在解決多模式混合信號的分解方面具有突出的優勢。李娟等[15]將SWD分解方法與平均差值算子方法相結合實現了齒輪箱的復合故障診斷。但是，未對閾值參數值PSth和Tth設置，導致很難達到預期效果，故建立能夠提取復合故障特征的準則，獲得到合理閾值參數，從而使 SWD將不同故障特征分解在不同的分量下，對于提高SWD的分解能力具有重要意義。

基于以上分析，本文構建基于平方包絡譜負熵的優化準則，采用在全局、局部搜索能力和穩定性方面較優的改進蝗蟲優化算法(improved grasshopper optimization algorithm, IGOA)對SWD中的閾值參數進行優化，實現對軸承早期復合故障特征的提取。

1 SWD方法

SWD作為一種智能分解方法，在參數化SWF的基礎上，可以從復合信號中迭代提取出主振蕩分量(oscillatory component ,OC)[14]。SWF的基本原理是將輸入信號x[n]作為群體的捕食路徑，輸出群體軌跡。在群濾波的模型中，獵物在第n步的位置表示為pprey[n]，M為群的數量，群的第i個個體在第n步的位置和速度分別為pi[n]和vi[n]。群個體在獵物過程受兩種作用力的控制，第一種作用力是由自身的驅動力產生，第二種作用力是由相互之間的吸引或排斥而產生。分別定義為式(1)和式(2)

(1)

(2)

式中，f()函數定義為

(3)

式中：sgn()為符號函數；ln()自然對數；d|為兩個群個體之間距離的絕對值；dcr為2個群個體之間相互不影響的臨界距離，即f(dcr)=0。一般來說，dcr可以用來控制群個體的分布，其大小為輸入信號的均方根值。

為跟蹤獵物，群個體必須通過改變自己的位置來更新其狀態，其第i個個體在第n步的速度和位置可表示為

vi[n]=vi[n-1]+δ(Fdr-i[n]+Fco-i[n])

(4)

pi[n]=pi[n-1]+δvi[n]

(5)

式中，δ為控制群個體的靈活性。

SWF的輸出為群的軌跡，可以通過式(6)表示

(6)

式中：β為尺度因子，一般取β=0.005；群參數β和M為控制群行為的關鍵決定因素。對于不同的信號，通過式(7)來選擇δ和M的值

(7)

式中：Yδ,M(k)和S(k)為序列Yδ,M(n)和S(n)傅里葉變換的幅值；Yδ,M(n)為群濾波在參數為δ和M時的輸出；S(n)為一個預設的非平穩信號，該信號可以是包含多個分量的復合信號。參數δ和M與各M頻率之間的關系為

(8)

(9)

通過對SWF的迭代執行，不斷分解最新的振蕩分量，直到當輸入信號中不包含任何足夠能量的振蕩模式時算法終止，即當連續兩次迭代的偏差小于閾值Tth時，SWD迭代終止。在每一次迭代中，選取能量譜密度幅值較大的頻帶作為提取新模式的候選分量。

由于SG(savitzky-golay)濾波器具有平滑能量譜的作用，為了提高計算效率，在找到最高峰值之前對分解信號進行SG濾波具有重要意義。同時，為了減少搜索區域，預先設置一個合適的峰值閾值PSth，其最優頻率ωm可以通過下式估計

(10)

(11)

式中，S為使用SG濾波xit[n]信號后的傅里葉變換。

通過對SWD分解原理的分析，SWD分解算法的具體步驟總結如下：

步驟1輸入信號x[n]并初始化閾值PSth和Tth；

步驟2離散化xit[n]←x[n]，xit[n]為離散數據序列，it為序列號。首先，設置it=0，然后指定y0[n]←xit[n]；

步驟3平滑y0[n]的能量譜，估計最優頻帶；

步驟4初始化群參數和δ；

步驟5當i=1時，通過SWF進行群濾波，yi[n]←FSWF(xit[n],M,δ)。

在SWD的分解過程中，閾值PSth和Tth對于SWF的分解過程起著關鍵作用，其大小決定了是否能夠合理分解。如果設置不合理，則會導致分解時間增加或分解分量增多。

2 提出的早期復合故障診斷方法

2.1 建立平方包絡譜負熵優化準則

當滾動軸承發生故障時產生循環瞬態沖擊，并引起振動信號瞬時能量的波動，而平方包絡方差可以用來檢測信號時域內的瞬時能量波動[16]。對于長度為N的信號x[n]，其平方包絡SSE(x[n])可表示為

(12)

但是平方包絡方差不能檢測周期性沖擊，而故障產生的沖擊往往是周期性的，對于周期性沖擊可以通過平方包絡譜來檢測，平方包絡譜定義為[17]

SSES(f)=TDFT(SSE(x[n]))

(13)

信息熵由于可以衡量信號的不規則性而得到廣泛應用，其也可以作為用來測量信號周期性沖擊強度的指標。因此，根據信息熵定義平方包絡譜熵為

(14)

式中，〈〉為取均值。由于平方包絡譜熵隨著信號的周期性沖擊越強，其值越小。為了使其變化性同周期性沖擊強度變化一致，定義平方包絡譜負熵為[18]

(15)

在智能算法優化參數的過程中，建立一個簡單、魯棒的優化準則對于提高SWD的分解能力非常重要[19]。滾動軸承振動信號的平方包絡譜在故障特征頻率及其多次諧頻處均有峰值，即故障特征頻率在平方包絡譜中是周期分布的。因此本文根據以上分析，提出以振動信號的平方包絡譜負熵為準則，進行自適應SWD分解，將分解信號分量的平均平方包絡負熵最大對應的閾值參數作為自適應SWD的分解參數，即

(16)

式中，k為分解后模態分量的個數。

2.2 蝗蟲優化算法及改進

Saremi等[20]提出了蝗蟲優化算法(grasshopper optimization algorithm, GOA) 。該算法由于原理簡單、收斂速度較快以及平衡全局和局部搜索能力強而得到一定應用，其數學模型為

Xi=r1Si+r2Gi+r3Ai

(17)

式中：Xi為第i只蝗蟲的位置；Si為在群體中第i只蝗蟲與其他蝗蟲的相互影響；Gi為第i只蝗蟲的重力引力；Ai為第i只蝗蟲所受風力影響；r1，r2和r3為[0,1]的隨機數，表示Si,Gi和Ai受不同環境影響的系數。

(18)

式中：dij=xj-xi|表示第i只蝗蟲和第j只蝗蟲之間的距離；dij=xj-xi|/dij為第i只蝗蟲和第j只蝗蟲之間的單位向量；s(r)群體中蝗蟲相互影響力，其函數表達式為

(19)

式中：f為吸引力常數；l為吸引力長度尺度參數。在蝗蟲群體覓食的過程中，形成了舒適區、排斥區和吸引區三種互區域，三種區域的范圍通過s(r)的值來界定：當s(r)>0時，蝗蟲間相互吸引，稱r的取值范圍為吸引區；當s(r)<0時，蝗蟲間相互排斥，稱r的取值范圍為排斥區；當s(r)=0時，蝗蟲之間既不相互吸引也不相互排斥，稱r的取值為舒適距離。但當r的取值過大時，s(r)的取值同樣接近于0，而此時r的取值并非舒適距離。此外，f和l二者的取值情況可以控制吸引區、排斥區及舒適區的分布情況，且二者取值通常較小，本文中取f=1.5，l=0.5，式(17)中G分量的表達式為

Gi=-geg

(20)

式中：g為重力常數；eg為朝向地心的單位向量。A分量的計算公式為

Ai=-uew

(21)

式中：u為漂移常數；ew為風向的單位向量。在式(17)中代入S，G和A的值，得到

(22)

式中：ubd為d維度變量的上限；lbd為d維度變量的下限；Td為目標中d維度變量的值，最外側的c用于平衡全局與局部尋優能力，內側的c用于控制舒適域、排斥域及吸引域的伸縮，c的值按線性遞減，其表達式為

(23)

式中：cmax為最大值；cmin為最小值；l為當前迭代次數；L為最大迭代次數。

在此，通過改進的圓混沌映射來代替式(23)來更新參數c的變化，從而提高GOA算法的全局搜索能力和穩定性，其表達式為

(24)

2.3 早期復合故障診斷實現步驟

通過平方包絡譜負熵準則，建立IGOA優化SWD的復合早期故障診斷方法如圖1所示。其具體步驟如下：

圖1 提出的早期復合故障診斷流程Fig.1 The flow of early composite fault diagnosis

步驟1輸入采集的軸承振動信號和初始化GOA的參數，蝗蟲的種群N=20、最大迭代次數L=10和閾值參數ubd=0.99，lbd=0.01，根據設定參數產生初始化蝗蟲群體；

步驟2計算輸入信號的能量譜，估計最優頻帶，初始化群參數M和δ；

步驟3進行SWD迭代分解，并計算分解模態分量的適應度值；

步驟4如果l≥L，迭代準止，否則轉入步驟3，且l=l+1，更新當前最優的個體；

步驟5保存輸出最優閾值參數，構建優化群分解(optimized swarm decomposition，OSWD)方法；

步驟6利用OSWD方法分解振動信號，獲得所有模態分量；

步驟7對分解后的分量進行包絡譜分析，得到故障特征頻率。

3 試驗驗證

3.1 仿真試驗

為了驗證所提方法的準確性和有效性，構建一組周期性脈沖序列信號并添加高斯白噪聲，其仿真信號模型為

y(t)=s(t)+n(t)

(25)

式中：s(t)為模擬故障沖擊的周期性沖擊分量；固有頻率fn=4 000 Hz；采樣點數N=4 096；采樣頻率為20 kHz；n(t)為高斯白噪聲。其中：外圈位移常數A0=3；阻尼系數g=0.05；故障特征頻率為100 Hz；產生的模擬外圈故障振動信號如圖2所示。內圈位移常數A0=2.5，阻尼系數g=0.1，故障特征頻率為125 Hz，產生內圈模擬故障振動信號如圖3所示。通過內外圈模擬信號并添加一組隨機擾動噪聲復合而成的故障信號，如圖4所示。與內外圈單故障信號相比，復合故障信號中反映內外圈故障的周期性沖擊基本被淹沒。對圖4所示的復合故障仿真信號直接進行包絡分析得到如圖5所示的結果。從圖5可知，外圈故障及倍頻較為突出，內圈故障特征較微弱，容易被忽略，可見直接進行包絡譜分析實現不了內外圈的復合故障診斷。

圖2 內圈故障時域波形Fig.2 Time domain waveform of inner ring fault

圖3 外圈故障時域波形Fig.3 Time domain waveform of outer ring fault

圖4 復合故障信號時域波形Fig.4 Time domain waveform of composite fault signal

圖5 復合故障信號包絡譜Fig.5 Envelope spectrum of composite fault signal

在此，設置蝗蟲的種群數N=20、最大迭代次數L=10以及閾值參數ubd=0.99，lbd=0.01。通過提出的OSWD方法對內外圈復合故障仿真信號進行故障特征提取，首先經過OSWD分解后得到的分量如圖6所示。從圖6可知，故障信號被分解為2個分量，其次對分解后的分量進行包絡譜分析如圖7所示。由圖7可知，在第1個分量中可以清晰地得到內外圈故障特征頻率的一倍頻和二倍頻；在第2個分量中同樣可以觀察到內外圈故障特征頻率的一倍頻和二倍頻。與理論故障頻率進行對比，易判斷出內外圈發生了故障。因此，OSWD對內外圈復合故障特征的分離是有效的。

圖6 OSWD分解后OC分量Fig.6 OSWD decomposition component

圖7 OSWD分解后OC分量包絡譜Fig.7 Envelope spectrum of OSWD decomposition component

為了驗證本文所提方法的有效性，設置SWD中2個閾值參數為默認值0.2，進行內外圈復合故障信號的分解和特征提取。首先，對復合故障仿真信號進行分解，分解后的分量和包絡譜如圖8和圖9所示。由圖可知，復合故障仿真信號被分解為5個分量，在第2個分量中得到了外圈故障特征頻率的一倍頻和內圈故障特征頻率的一倍頻與二倍頻；在第3個分量中得到了內外圈故障特征頻率的一倍頻和二倍頻。因此，通過SWD方法也可以判斷出內圈產生故障，但是SWD的閾值參數是手動選擇而不具備自適應性，同時分解得到了3個冗余分解分量，這也必然導致分解過程中計算時間的增加。

圖8 SWD分解后OC分量Fig.8 SWD decomposition component

圖9 SWD分解后OC分量包絡譜Fig.9 Envelope spectrum of SWD decomposition component

為了進一步驗證OSWD的優越性，本文繼續采用當前在模態混疊和端點效應問題方面較有優勢的VMD分解方法對復合故障仿真信號進行分解，設置參數k=6，α=8 000，其中k一般認為值為5～6時可以實現較好的分解，α為VMD中的默認值。分解得到的分量和分量的包絡譜如圖10和11所示。由圖可知，在第2個和第3個分量的包絡譜中得到了內外圈故障特征頻率的一倍頻和二倍頻。但是選取不同的[k,α]組合參數，在VMD分解中得到的結果不同，此結果為本文經過反復選擇參數測試得到的VMD分解提取效果較優的結果。

圖10 VMD分解后分量Fig.10 VMD decomposition component

圖11 VMD分解后分量包絡譜Fig.11 Envelope spectrum of VMD decomposition component

通過上述三種方法對仿真信號進行復合故障特征提取分析，可以得知OSWD在復合故障診斷方面具有一定的優越性，能夠實現復合故障的診斷，有必要進行進一步的應用研究。

3.2 工程驗證

將所提方法應用于某石油化工企業生產裝置中關鍵機組風機軸承的早期復合故障診斷，驗證該方法的有效性和優越性。在風機運行過程中，其型號為SKF6030C3的風機驅動端軸承內外圈發生早期磨損，如圖12所示。從圖12可知，其故障程度還比較輕微，屬于早期故障，相比內圈，外圈發生了程度更輕微。振動信號通過安裝在驅動端的水平和垂直加速度傳感器獲得，其中電機轉速為2 985 r/min，信號的采樣頻率為25.6 kHz，數據采樣點為16 384。電機驅動風機運行的系統結構圖如圖13所示。該軸承不同部件的理論故障特征頻率，如表1所示。

圖12 內外圈早期復合故障Fig.12 Early composite failure of inner and outer rings

圖13 關鍵機組風機運行系統結構圖Fig.13 The system structure of key unit operation

表1 軸承故障特征頻率

風機振動信號時域波形如圖14所示。從圖14可知，不太規律的周期性沖擊，但僅通過時域波形很難判斷出內外圈故障。對其直接進行包絡譜分析如圖15所示。從圖15可知，僅有內圈故障特征頻率的一倍頻比較突出而外圈故障可能是因為內圈故障發生后才產生而比較微弱，因此故障特征頻率難以發現。

圖14 振動信號時域波形Fig.14 Time domain waveform of vibration signal

圖15 振動信號包絡譜Fig.15 Envelope spectrum of vibration signal

設置IGOA的初始化參數與仿真案例中的參數相同，通過IGOA優化SWD得到的最佳參數為[0.950 55,0.202 66]。應用提出的OSWD方法對風機軸承的振動信號進行分解得到的分量，如圖16所示。從圖16可知，振動信號被分解為2個分量，然后對2個分量進行包絡譜分析得到的結果，如圖17所示。在分解后的第1個分量的包絡譜中較為突出的得到外圈故障特征頻率。同時，在第2個分量的包絡譜中可以明顯的看到轉頻、轉頻倍頻和外圈故障特征頻率。由此可見，通過本文提出的方法可以實現關鍵機組軸承內外圈早期復合故障的特征提取和診斷。

圖16 OSWD分解后OC分量Fig.16 OSWD decomposition component

圖17 OSWD分解后OC分量包絡譜Fig.17 Envelope spectrum of OSWD decomposition component

為驗證本文所提方法的有效性和優越性，采用SWD分解方法對關鍵機組軸承振動信號進行分析，按原始文獻設置閾值參數為默認值0.2，振動信號被SWD分解為4個分量如圖18所示。對其進行包絡譜分析如圖19所示。由圖19可知，在第1個分量的包絡譜中獲得了外圈故障特征頻率，在第2、第3、第4個分量中都得到了內圈故障特征頻率和轉頻的倍頻，實質上在第2個分量已經得到了內圈故障特征頻率，但是因為未通過建立的平方包絡負熵準則優化來實現自適應分解導致分解后產生冗余分量，增加了分解時間。

圖18 SWD分解后OC分量Fig.18 SWD decomposition component

圖19 SWD分解后OC分量包絡譜Fig.19 Envelope spectrum of SWD decomposition component

為了進一步對比分析，設置VMD參數k=6，α=8 000，繼續采用VMD分解方法對采集的風機振動信號進行分析，得到分解后的分量和分量包絡譜，分別如圖20和圖21所示。由圖20可知，在分解后第2個分量中得到了外圈的故障特征頻率，而內圈故障特征頻率沒有在任何分量中有效提取，因此通過VMD方法無法判斷出內圈故障，故無法實現復合故障特征的有效提取。

圖20 VMD分解后分量Fig.20 VWD decomposition component

圖21 VMD分解后包絡譜Fig.21 Envelope spectrum of VMD decomposition component

通過以上試驗方法的對比分析可以得到，通過本文提出的OSWD方法對于關鍵機組軸承早期復合故障振動信號進行自適應分解，能夠實現早期復合故障特征的有效提取，解決了早期復合故障信號由于微弱、多信號相互耦合以及強噪聲干擾導致的故障特征難以提取的問題，實現了滾動軸承早期復合故障的準確診斷。

4 結 論

本文提出了基于平方包絡譜負熵準則的軸承早期復合故障特征提取方法，主要結論如下：

(1)針對軸承早期復合故障診斷中故障特征難以提取的問題，提出了基于OSWD的復合故障特征提取方法。

(2)建立了能夠檢測周期性故障沖擊能量的平方包絡譜負熵。

(3)就閾值參數在SWD在分解過程中對SWD分解效果影響較大的問題，通過建立的平方包絡譜負熵準則，應用IGOA優化SWD實現了閾值參數組合優化。

(4)通過與SWD、VMD對比分析表明，該方法在復合故障診斷方面的性能更加優越，具備較好的應用性。