基于參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解的轉(zhuǎn)子故障診斷

唐貴基，周 翀，龐 彬，李楠楠

(華北電力大學(xué) 能源動力與機械工程學(xué)院，河北 保定 071003)

Keywords:rotor;fault diagnosis;time-varying filtering;empirical mode decomposition;parameter optimization;Hilbert transform(HT)

如何提取出轉(zhuǎn)子不平衡、碰磨、油膜渦動、裂紋等典型非線性、非平穩(wěn)故障信號的故障特征，以準(zhǔn)確診斷出轉(zhuǎn)子故障類型，是轉(zhuǎn)子系統(tǒng)故障診斷研究的熱點和難點[1]。

經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是一種自適應(yīng)的非線性、非平穩(wěn)信號分解方法[2]，在轉(zhuǎn)子故障診斷中得到了廣泛的應(yīng)用[3-4]。但是，該方法存在模式混疊的現(xiàn)象。為解決這一問題,相關(guān)學(xué)者提出了集合經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)[5]、噪聲輔助多元經(jīng)驗?zāi)B(tài)分解(Noise Assisted Multivariate Empirical Mode Decomposition,N-A MEMD)[6]等噪聲輔助方法,變分模態(tài)分解(Variational Mode Decomposition,VMD)[7]、經(jīng)驗小波變換(Empirical Wavelet Transform，EWT)[8]等基于濾波的方法。但噪聲輔助方法存在噪聲幅值、集合個數(shù)等參數(shù)選取困難，無法分離倍頻程范圍內(nèi)的分量等不足。經(jīng)驗?zāi)B(tài)分解本質(zhì)上是通過迭代對信號進行低通濾波以去除局部均值，從而分離出本征模態(tài)函數(shù)?，F(xiàn)有基于濾波的方法構(gòu)造截止頻率恒定的低通濾波器進行濾波，不適合局部均值隨時間改變的非平穩(wěn)信號的分析。

基于此，Li等[9]提出了時變?yōu)V波經(jīng)驗?zāi)B(tài)分解(Time Varying Filtering EMD,TVFEMD)法。該方法克服了經(jīng)驗?zāi)B(tài)分解過程中的模式混疊現(xiàn)象，提高了信號在噪聲和低采樣頻率下的魯棒性，并且能夠保持信號的時變特征，適合于非線性、非平穩(wěn)信號的分析。林近山等[10]使用時變?yōu)V波經(jīng)驗?zāi)B(tài)分解的方法診斷齒輪箱的故障，取得了良好的效果。但是TVFEMD使用過程中，需人為設(shè)定帶寬閾值和B樣條階數(shù)兩個參數(shù)，具有較大的主觀性，而且無法預(yù)知在所選參數(shù)組合下分解的有效性。因此，研究能夠?qū)崿F(xiàn)參數(shù)自動尋優(yōu)的時變?yōu)V波經(jīng)驗?zāi)B(tài)分解方法，是十分必要的。Zhang 等[11]針對軸承、齒輪故障提出了改進的時變?yōu)V波經(jīng)驗?zāi)B(tài)分解方法實現(xiàn)參數(shù)自動尋優(yōu)，但是衡量信號沖擊性的加權(quán)峭度指標(biāo)不適用于轉(zhuǎn)子故障信號。尋找新的衡量指標(biāo)并采用有效的尋優(yōu)算法以實現(xiàn)轉(zhuǎn)子故障信號TVFEMD過程的自動尋優(yōu)十分必要。

本文以Hilbert邊際譜部分頻帶能量比作為適應(yīng)度函數(shù)，采用粒子群算法進行參數(shù)自動尋優(yōu)，提出了參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解的方法，再將該方法和希爾伯特變換(Hilbert Transform,HT)相結(jié)合對轉(zhuǎn)子故障信號進行時頻分析，實現(xiàn)對轉(zhuǎn)子故障的診斷。

1 時變?yōu)V波經(jīng)驗?zāi)B(tài)分解原理

時變?yōu)V波經(jīng)驗?zāi)B(tài)分解本質(zhì)上是通過構(gòu)造截止頻率隨時間變化的低通濾波器，來完成EMD分解過程中的迭代去除均值操作，并以局部窄帶信號代替本征模態(tài)函數(shù)作為迭代停止條件。由于對于給定任意多分量信號x(t)，可以表示為一個雙分量信號

x(t)=A(t)ejφ(t)=a1(t)ejφ1(t)+a2(t)ejφ2(t)

(1)

因此，只需考慮雙分量信號的分解過程。對雙分量信號進行時變?yōu)V波經(jīng)驗?zāi)B(tài)分解的基本步驟如下：

步驟1對x(t)進行希爾伯特變換，得到復(fù)解析信號的幅值A(chǔ)(t)和相位φ(t)，對瞬時相位求導(dǎo)可得瞬時頻率φ′(t)。

(2)

步驟2確定幅值曲線A(t)的極小值和極大值所在的時刻{tmin},{tmax}以及幅值A(chǔ)({tmin}),A({tmax})。

步驟3對A(t)的極值點分別進行插值，所得曲線分別為β1(t)和β2(t)。

步驟4根據(jù)式(3)、式(4)計算出瞬時均值a1(t)和瞬時包絡(luò)a2(t)。

a1(t)=[β1(t)+β2(t)]/2

(3)

a2(t)=[β2(t)-β1(t)]/2

(4)

步驟5令

(5)

(6)

(7)

(1)定義信號x(t)的極大值點的時間序列ui,i=1,2,3,...;

如果

(8)

步驟7根據(jù)調(diào)整后的截止頻率重構(gòu)信號

(9)

將h(t)的極值點作為節(jié)點，將h(t)分成n段，每段步長為m。其中，n稱之為B樣條函數(shù)的階數(shù)。使用式(10)～式(13)進行B樣條插值逼近，逼近結(jié)果記為m(t)，代表局部均值曲線。

(10)

(11)

(12)

(13)

式中：[.]↓m為降采樣，即每隔m個點進行采樣；[.]↑m為過采樣，即在每兩個采樣點之間插入m個采樣點。

步驟8判斷是否滿足停止準(zhǔn)則θ(t)<ξ，其中ξ為給定帶寬閾值。如滿足，x(t)為IMF；如不滿足，令x(t)=x(t)-m(t)，重復(fù)步驟1～步驟7，直到滿足停止準(zhǔn)則為止，滿足停止準(zhǔn)則的x(t)為IMF。

(14)

(15)

(16)

2 參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解原理

TVFEMD執(zhí)行過程中，需人為指定帶寬閾值ξ和B樣條階數(shù)n兩個參數(shù)，存在較大盲目性和不確定性。為實現(xiàn)參數(shù)的自動尋優(yōu)，需首先定義目標(biāo)函數(shù)，確定尋優(yōu)算法。下面就本文所提Hilbert邊際譜部分頻帶能量比作為目標(biāo)函數(shù)的粒子群參數(shù)優(yōu)化方法進行介紹。

2.1 Hilbert邊際譜部分頻帶能量比

Hilbert 邊際譜以概率的形式表示了每個頻率上分布的整個時間序列累計的振幅或能量[12]。在邊際譜中，不同的轉(zhuǎn)子故障類型在轉(zhuǎn)頻及其分?jǐn)?shù)諧波、倍頻處表現(xiàn)出不同的特征，但故障特征頻率集中于一定的時間和頻率范圍內(nèi)[13]。參考文獻[14]提出的軸承故障包絡(luò)譜故障特征能量比指標(biāo)來衡量解調(diào)效果，本文提出Hilbert邊際譜部分頻帶能量比的指標(biāo)來定量衡量轉(zhuǎn)子故障信號TVFEMD的效果并作為粒子群算法尋優(yōu)的適應(yīng)度函數(shù)。該指標(biāo)定義為邊際譜0～5倍轉(zhuǎn)頻頻帶范圍內(nèi)能量與邊際譜總能量的比值，其計算公式為

(17)

2.2 粒子群算法

粒子群算法是一種具備全局尋優(yōu)能力的方法。X=(X1,X2,...,XM)代表M個粒子組成的種群，D維向量Xi=(xi1,xi2,...,xiD)表示第i個粒子的位置，向量Vi=(vi1,vi2,...,viD表示第i個粒子的速度，Qi=(qi1,qi2,...,qiD)為個體局部極值，G=(g1,g2,...,gD)為種群全局極值，各粒子通過式(18)和式(19)更新位置和速度

(18)

(19)

式中：ω為慣性權(quán)重;d=1,2,...,D;i=1,2,...,M;a為當(dāng)前迭代次數(shù);c1和c2為加速度因子;η為介于[0,1]的隨機數(shù)。

2.3 參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解流程

步驟1初始化粒子群各參數(shù)。

步驟2以ξ和n的參數(shù)組合作為位置坐標(biāo)，初始化粒子群的位置并隨機初始化粒子群的速度。

步驟3在不同粒子位置下對信號進行TVFEMD，計算各個位置下的Hilbert邊際譜部分頻帶能量比。

步驟4對比各個位置Hilbert邊際譜部分頻帶能量比大小，初始化個體局部極值和整體局部極值。

步驟5利用式(18)和式(19)更新粒子群的位置和速度。

步驟6迭代循環(huán)，轉(zhuǎn)至步驟3，直至迭代次數(shù)達到最大設(shè)定值后輸出最佳適應(yīng)度值和粒子位置。

3 診斷流程

基于參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解和希爾伯的轉(zhuǎn)子故障診斷方法基本步驟如下：

步驟1確定參數(shù)尋優(yōu)范圍，初始化粒子群各參數(shù)。

Zhang等推薦帶寬閾值0<ξ≤1。筆者經(jīng)多次實驗發(fā)現(xiàn)，隨著ξ的取值接近于0，計算耗時增加，而且邊際譜出現(xiàn)較多的特征頻率外的成分干擾。一般ξ取值達到約0.1時即可獲得較好的分解效果，邊際譜圖更加干凈，計算速度較快。因此，本文中帶寬閾值的尋優(yōu)范圍為0.1≤ξ≤1，B樣條階數(shù)的尋優(yōu)范圍在5≤n≤30，獲得的參數(shù)組合在實際工程應(yīng)用中已足夠達到較好的分解性能。

參考文獻[15-16]，并經(jīng)過多次數(shù)值實驗，本文中粒子群的參數(shù)設(shè)置見表1。

表1 粒子群算法參數(shù)設(shè)置Tab.1 Parameter setting of particle swarm optimization

步驟2使用粒子群尋優(yōu)的最佳參數(shù)組合進行TVFEMD，獲得若干IMF(Intrinsic Mode Function)。

步驟3對IMF進行HT，根據(jù)瞬時信息得到故障信號的Hilbert時頻圖和邊際譜。

步驟4根據(jù)時頻信息診斷轉(zhuǎn)子故障類型。

4 實驗驗證

采用如圖1所示的轉(zhuǎn)子實驗臺驗證本文所提方法的有效性和優(yōu)越性。該實驗臺由圓盤轉(zhuǎn)子、動壓軸承、調(diào)速電機控制器、預(yù)緊力支架、ZonicBook/618E型信號采集設(shè)備、后處理計算機等設(shè)備組成。

4.1 恒定轉(zhuǎn)速的轉(zhuǎn)子不平衡故障

在轉(zhuǎn)子圓盤轉(zhuǎn)子45°的標(biāo)記處安裝1 g的不平衡質(zhì)量來模擬轉(zhuǎn)子不平衡故障。轉(zhuǎn)子的轉(zhuǎn)速為3 000 r/min，采樣頻率為2 560 Hz。分別在水平和垂直兩個方向安裝傳感器采集信號，從水平方向的信號中任選連續(xù)的2 000點作為待分析信號。

圖1 Bently RK4轉(zhuǎn)子實驗臺和信號采集器Fig.1 Bently RK4 rotor test rig and signal collector

表2為在人為設(shè)定的不同參數(shù)組合下進行TVFEMD所得Hilbert邊際譜部分頻帶的能量比值。

表2 不平衡故障不同參數(shù)組合的Hilbert邊際譜部分頻帶能量比Tab.2 Partial spectrum energy ratio of Hilbert marginal spectrum with different parameter combinations for unbalanced faults

為對比分析不同參數(shù)取值時分解效果的差異，同時，驗證使用邊際譜部分頻帶能量比指標(biāo)作為適應(yīng)度函數(shù)衡量分解效果的有效性，現(xiàn)選取表2中能量比最大的參數(shù)組合1、居中的任意參數(shù)組合3以及最小的參數(shù)組合8，分別在三組參數(shù)組合下進行TVFEMD，得到Hilbert邊際譜如圖2所示。

圖2(a)中可見清晰的基頻及二倍頻成分，且集聚性好，幅值較大；圖2(b)雖然也存在清晰基頻及二倍頻成分，但與圖2(a)相比幅值較低；圖2(c)則完全無法分離出清晰的基頻及倍頻成分，分解的效果很差。可見隨著邊際譜部分頻帶能量比逐漸降低，其分解效果也隨之變差。由此驗證，使用Hilbert邊際譜部分頻帶能量比作為衡量分離效果的指標(biāo)是可行的。

同時，表2也表明不同的參數(shù)組合下進行TVFEDM所得Hilbert邊際譜的部分頻帶能量比差距較大，且無明顯規(guī)律性。因此，為達到良好的分解效果，適當(dāng)?shù)膮?shù)組合的選取十分必要。而人為選擇參數(shù)存在較大的主觀性和盲目性，使用優(yōu)化算法獲取最佳參數(shù)組合無疑具有重要的工程意義。

圖3為使用粒子群算法進行參數(shù)尋優(yōu)時Hilbert邊際譜部分頻帶能量比隨進化代數(shù)的變化曲線。粒子群進化到第5代即得到了信號最大的Hilbert邊際譜部分頻帶能量比0.011 8。尋優(yōu)所得最佳參數(shù)組合為ξ=0.31,n=13。在該參數(shù)組合下進行TVFEMD，結(jié)果如圖4所示。

圖2 不同參數(shù)組合下不平衡故障信號分解結(jié)果Fig.2 Decomposition results of unbalanced fault signals under different combinations of parameters

圖3 不平衡故障Hilbert 邊際譜部分頻帶能量比隨進化代數(shù)的變化Fig.3 Hilbert partial spectrum energy ratio of unbalanced fault changes with evolutionary algebra

圖4 不平衡故障信號使用本文方法分析所得Hilbert邊際譜和時頻圖Fig.4 The Hilbert marginal spectrum and time-frequency diagram of the unbalanced fault signal are analyzed using this method

從圖4可知，明顯的50 Hz的轉(zhuǎn)頻成分，以及與基頻相比微不足道的二倍頻成分，不平衡故障特征明顯。與圖2(a)相比，基頻的幅值更大，二倍頻成分進一步壓縮，分解的效果更好。圖4中的時頻圖，基頻集聚性強，不平衡特征明顯。使用粒子群算法尋優(yōu)的參數(shù)組合進行TVFEMD，能夠有效的診斷不平衡故障，實現(xiàn)良好的分解效果。

為驗證本文所提方法相比其他方法的優(yōu)越性，分別使用原始EMD方法、文獻[9]中的指定參數(shù)進行TVFEMD、文獻[11]所提參數(shù)尋優(yōu)方法對不平衡故障信號進行分析，結(jié)果如圖5所示。其中，文獻[9]采用ξ=0.1，n=28的參數(shù)組合；對實驗采集的轉(zhuǎn)子不平衡故障信號采用文獻[11]所提基于相關(guān)峭度指標(biāo)進行尋優(yōu)所得參數(shù)組合為ξ=0.37，n=20。

圖5 不同方法處理不平衡故障信號分析結(jié)果Fig.5 Analysis results of unbalanced fault signals processed by different methods

圖5(a)所示采用原始EMD方法分析結(jié)果中，僅能發(fā)現(xiàn)幅值極小的基頻附近成分。圖5(b)所示采用文獻[9]方法分析結(jié)果中，除基頻外還存在幅值較為突出的二倍頻成分。圖5(c)所示采用文獻[11]方法分析結(jié)果與圖4所示采用本文方法進行分析的結(jié)果相比，集聚性較差，且邊際譜基頻幅值較低?？傊?，采用其他三種現(xiàn)有方法，無法突出故障特征頻率，實現(xiàn)對故障類型的準(zhǔn)確判斷，而使用本文方法對故障信號進行處理故障特征頻率突出，能夠準(zhǔn)確對不平衡故障進行診斷，從而驗證了本文所提方法在轉(zhuǎn)子恒定轉(zhuǎn)速不平衡故障診斷中的優(yōu)越性。

4.2 變轉(zhuǎn)速的轉(zhuǎn)子油膜渦動故障

為驗證本文所提方法轉(zhuǎn)子變轉(zhuǎn)速工況下的適用性，采集實驗臺從1 400～2 000 r/min升速過程油膜渦動故障的振動信號，采樣頻率為1 280 Hz,采樣點數(shù)為3 000。圖6和圖7分別為升速油膜渦動故障的時域波形和頻譜。頻譜中存在嚴(yán)重的頻率模糊現(xiàn)象，無法準(zhǔn)確提取故障特征頻率。

圖6 升速油膜渦動故障的時域波形Fig.6 Time domain waveform of oil film whirling fault at rising speed

圖7 升速油膜渦動故障的頻譜Fig.7 Spectrum of oil film whirling fault at rising speed

圖8為使用本文所提的方法對故障信號進行分析所得Hilbert邊際譜部分頻帶能量比隨進化代數(shù)的變化曲線。其中，對于變轉(zhuǎn)速轉(zhuǎn)子故障信號，粒子群進化到第3代即可得到最大的幅值能量比0.037 45。尋優(yōu)所得最佳的參數(shù)組合ξ=0.12，n=20。在該參數(shù)組合下對故障信號進行分解，所得Hilbert時頻圖如圖9(a)所示。圖中能夠清晰提取轉(zhuǎn)頻及其二分頻成分，符合轉(zhuǎn)子油膜渦動故障的特征，且集聚性較好，具有良好的分解效果。

圖9(b)為采用原始EMD分析轉(zhuǎn)子故障信號的Hilbert時頻圖，譜圖中存在嚴(yán)重的頻率模糊現(xiàn)象，無法準(zhǔn)確提取故障特征頻率。圖9(c)為采用文獻[9]指定的參數(shù)組合ξ=0.1，n=28進行TVFEMD和HT所得的Hilbert時頻圖。雖然也能較為清晰提取轉(zhuǎn)頻及其二分頻，但特征頻率處能量明顯小于圖9(a)采用本文方法分析所得結(jié)果。采用文獻[11]所提出的相關(guān)峭度指標(biāo)進行尋優(yōu)，得到最優(yōu)參數(shù)組合ξ=0.42，n=12，在該參數(shù)組合下進行TVFEMD和HT，所得Hilbert時頻圖如圖9(d)所示。雖然能夠準(zhǔn)確提取轉(zhuǎn)子的二分頻，卻無轉(zhuǎn)頻成分，與頻譜所示的特征頻率不一致，由此造成故障特征信息的丟失。總之，本文所提方法能夠準(zhǔn)確提取非平穩(wěn)轉(zhuǎn)子故障信號的特征頻率，相比于其他三種方法具有明顯的優(yōu)越性。

圖8 油膜渦動故障Hilbert 邊際譜部分頻帶能量比隨進化代數(shù)的變化Fig.8 The variation of Hilbert marginal spectrum and partial band energy ratio of oil film whirling fault with evolution algebra

圖9 升速油膜渦動故障信號不同方法分析Hilbert時頻圖Fig.9 Analysis of Hilbert time spectrum by different methods of rising oil film whirl fault signal

5 結(jié) 論

本文提出了一種基于參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解和希爾伯特變換診斷轉(zhuǎn)子故障的方法。該方法以Hilbert邊際譜部分頻帶能量比來定量衡量TVFEMD的分解效果，并作為粒子群尋優(yōu)算法的適應(yīng)度函數(shù)實現(xiàn)參數(shù)尋優(yōu)；通過將分解所得的IMF進行HT得到的瞬時信息來診斷轉(zhuǎn)子故障類型。分別將該方法應(yīng)用于恒定轉(zhuǎn)速下的轉(zhuǎn)子不平衡和變轉(zhuǎn)速下的油膜渦動兩種典型轉(zhuǎn)子故障，驗證了該方法診斷轉(zhuǎn)子故障的有效性和優(yōu)越性，得到結(jié)論如下：

(1)使用Hilbert邊際譜部分頻帶能量比作為衡量TVFEMD分離效果的指標(biāo)是可行的。

(2)使用參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解結(jié)合希爾伯特變換對故障信號進行分析，能夠突出故障特征頻率成分，提高分解效果，從而有效的診斷非平穩(wěn)轉(zhuǎn)子故障。

(3)本文所提參數(shù)優(yōu)化時變?yōu)V波經(jīng)驗?zāi)B(tài)分解方法相比于原始經(jīng)驗?zāi)B(tài)分解和現(xiàn)有方法相比能夠克服模式混疊現(xiàn)象，更加準(zhǔn)確提取故障特征頻率，從而有效診斷轉(zhuǎn)子不平衡、油膜渦動等典型故障。