唐貴基,周 翀,龐 彬,李楠楠
(華北電力大學(xué) 能源動(dòng)力與機(jī)械工程學(xué)院,河北 保定 071003)
Keywords:rotor;fault diagnosis;time-varying filtering;empirical mode decomposition;parameter optimization;Hilbert transform(HT)
如何提取出轉(zhuǎn)子不平衡、碰磨、油膜渦動(dòng)、裂紋等典型非線性、非平穩(wěn)故障信號(hào)的故障特征,以準(zhǔn)確診斷出轉(zhuǎn)子故障類型,是轉(zhuǎn)子系統(tǒng)故障診斷研究的熱點(diǎn)和難點(diǎn)[1]。
經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是一種自適應(yīng)的非線性、非平穩(wěn)信號(hào)分解方法[2],在轉(zhuǎn)子故障診斷中得到了廣泛的應(yīng)用[3-4]。但是,該方法存在模式混疊的現(xiàn)象。為解決這一問(wèn)題,相關(guān)學(xué)者提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)[5]、噪聲輔助多元經(jīng)驗(yàn)?zāi)B(tài)分解(Noise Assisted Multivariate Empirical Mode Decomposition,N-A MEMD)[6]等噪聲輔助方法,變分模態(tài)分解(Variational Mode Decomposition,VMD)[7]、經(jīng)驗(yàn)小波變換(Empirical Wavelet Transform,EWT)[8]等基于濾波的方法。但噪聲輔助方法存在噪聲幅值、集合個(gè)數(shù)等參數(shù)選取困難,無(wú)法分離倍頻程范圍內(nèi)的分量等不足。經(jīng)驗(yàn)?zāi)B(tài)分解本質(zhì)上是通過(guò)迭代對(duì)信號(hào)進(jìn)行低通濾波以去除局部均值,從而分離出本征模態(tài)函數(shù)?,F(xiàn)有基于濾波的方法構(gòu)造截止頻率恒定的低通濾波器進(jìn)行濾波,不適合局部均值隨時(shí)間改變的非平穩(wěn)信號(hào)的分析。
基于此,Li等[9]提出了時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解(Time Varying Filtering EMD,TVFEMD)法。該方法克服了經(jīng)驗(yàn)?zāi)B(tài)分解過(guò)程中的模式混疊現(xiàn)象,提高了信號(hào)在噪聲和低采樣頻率下的魯棒性,并且能夠保持信號(hào)的時(shí)變特征,適合于非線性、非平穩(wěn)信號(hào)的分析。林近山等[10]使用時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解的方法診斷齒輪箱的故障,取得了良好的效果。但是TVFEMD使用過(guò)程中,需人為設(shè)定帶寬閾值和B樣條階數(shù)兩個(gè)參數(shù),具有較大的主觀性,而且無(wú)法預(yù)知在所選參數(shù)組合下分解的有效性。因此,研究能夠?qū)崿F(xiàn)參數(shù)自動(dòng)尋優(yōu)的時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解方法,是十分必要的。Zhang 等[11]針對(duì)軸承、齒輪故障提出了改進(jìn)的時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解方法實(shí)現(xiàn)參數(shù)自動(dòng)尋優(yōu),但是衡量信號(hào)沖擊性的加權(quán)峭度指標(biāo)不適用于轉(zhuǎn)子故障信號(hào)。尋找新的衡量指標(biāo)并采用有效的尋優(yōu)算法以實(shí)現(xiàn)轉(zhuǎn)子故障信號(hào)TVFEMD過(guò)程的自動(dòng)尋優(yōu)十分必要。
本文以Hilbert邊際譜部分頻帶能量比作為適應(yīng)度函數(shù),采用粒子群算法進(jìn)行參數(shù)自動(dòng)尋優(yōu),提出了參數(shù)優(yōu)化時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解的方法,再將該方法和希爾伯特變換(Hilbert Transform,HT)相結(jié)合對(duì)轉(zhuǎn)子故障信號(hào)進(jìn)行時(shí)頻分析,實(shí)現(xiàn)對(duì)轉(zhuǎn)子故障的診斷。
時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解本質(zhì)上是通過(guò)構(gòu)造截止頻率隨時(shí)間變化的低通濾波器,來(lái)完成EMD分解過(guò)程中的迭代去除均值操作,并以局部窄帶信號(hào)代替本征模態(tài)函數(shù)作為迭代停止條件。由于對(duì)于給定任意多分量信號(hào)x(t),可以表示為一個(gè)雙分量信號(hào)
x(t)=A(t)ejφ(t)=a1(t)ejφ1(t)+a2(t)ejφ2(t)
(1)
因此,只需考慮雙分量信號(hào)的分解過(guò)程。對(duì)雙分量信號(hào)進(jìn)行時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解的基本步驟如下:
步驟1對(duì)x(t)進(jìn)行希爾伯特變換,得到復(fù)解析信號(hào)的幅值A(chǔ)(t)和相位φ(t),對(duì)瞬時(shí)相位求導(dǎo)可得瞬時(shí)頻率φ′(t)。
(2)
步驟2確定幅值曲線A(t)的極小值和極大值所在的時(shí)刻{tmin},{tmax}以及幅值A(chǔ)({tmin}),A({tmax})。
步驟3對(duì)A(t)的極值點(diǎn)分別進(jìn)行插值,所得曲線分別為β1(t)和β2(t)。
步驟4根據(jù)式(3)、式(4)計(jì)算出瞬時(shí)均值a1(t)和瞬時(shí)包絡(luò)a2(t)。
a1(t)=[β1(t)+β2(t)]/2
(3)
a2(t)=[β2(t)-β1(t)]/2
(4)
步驟5令
(5)
(6)
(7)
(1)定義信號(hào)x(t)的極大值點(diǎn)的時(shí)間序列ui,i=1,2,3,...;
如果
(8)
步驟7根據(jù)調(diào)整后的截止頻率重構(gòu)信號(hào)
(9)
將h(t)的極值點(diǎn)作為節(jié)點(diǎn),將h(t)分成n段,每段步長(zhǎng)為m。其中,n稱之為B樣條函數(shù)的階數(shù)。使用式(10)~式(13)進(jìn)行B樣條插值逼近,逼近結(jié)果記為m(t),代表局部均值曲線。
(10)
(11)
(12)
(13)
式中:[.]↓m為降采樣,即每隔m個(gè)點(diǎn)進(jìn)行采樣;[.]↑m為過(guò)采樣,即在每?jī)蓚€(gè)采樣點(diǎn)之間插入m個(gè)采樣點(diǎn)。
步驟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)
TVFEMD執(zhí)行過(guò)程中,需人為指定帶寬閾值ξ和B樣條階數(shù)n兩個(gè)參數(shù),存在較大盲目性和不確定性。為實(shí)現(xiàn)參數(shù)的自動(dòng)尋優(yōu),需首先定義目標(biāo)函數(shù),確定尋優(yōu)算法。下面就本文所提Hilbert邊際譜部分頻帶能量比作為目標(biāo)函數(shù)的粒子群參數(shù)優(yōu)化方法進(jìn)行介紹。
Hilbert 邊際譜以概率的形式表示了每個(gè)頻率上分布的整個(gè)時(shí)間序列累計(jì)的振幅或能量[12]。在邊際譜中,不同的轉(zhuǎn)子故障類型在轉(zhuǎn)頻及其分?jǐn)?shù)諧波、倍頻處表現(xiàn)出不同的特征,但故障特征頻率集中于一定的時(shí)間和頻率范圍內(nèi)[13]。參考文獻(xiàn)[14]提出的軸承故障包絡(luò)譜故障特征能量比指標(biāo)來(lái)衡量解調(diào)效果,本文提出Hilbert邊際譜部分頻帶能量比的指標(biāo)來(lái)定量衡量轉(zhuǎn)子故障信號(hào)TVFEMD的效果并作為粒子群算法尋優(yōu)的適應(yīng)度函數(shù)。該指標(biāo)定義為邊際譜0~5倍轉(zhuǎn)頻頻帶范圍內(nèi)能量與邊際譜總能量的比值,其計(jì)算公式為
(17)
粒子群算法是一種具備全局尋優(yōu)能力的方法。X=(X1,X2,...,XM)代表M個(gè)粒子組成的種群,D維向量Xi=(xi1,xi2,...,xiD)表示第i個(gè)粒子的位置,向量Vi=(vi1,vi2,...,viD表示第i個(gè)粒子的速度,Qi=(qi1,qi2,...,qiD)為個(gè)體局部極值,G=(g1,g2,...,gD)為種群全局極值,各粒子通過(guò)式(18)和式(19)更新位置和速度
(18)
(19)
式中:ω為慣性權(quán)重;d=1,2,...,D;i=1,2,...,M;a為當(dāng)前迭代次數(shù);c1和c2為加速度因子;η為介于[0,1]的隨機(jī)數(shù)。
步驟1初始化粒子群各參數(shù)。
步驟2以ξ和n的參數(shù)組合作為位置坐標(biāo),初始化粒子群的位置并隨機(jī)初始化粒子群的速度。
步驟3在不同粒子位置下對(duì)信號(hào)進(jìn)行TVFEMD,計(jì)算各個(gè)位置下的Hilbert邊際譜部分頻帶能量比。
步驟4對(duì)比各個(gè)位置Hilbert邊際譜部分頻帶能量比大小,初始化個(gè)體局部極值和整體局部極值。
步驟5利用式(18)和式(19)更新粒子群的位置和速度。
步驟6迭代循環(huán),轉(zhuǎn)至步驟3,直至迭代次數(shù)達(dá)到最大設(shè)定值后輸出最佳適應(yīng)度值和粒子位置。
基于參數(shù)優(yōu)化時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解和希爾伯的轉(zhuǎn)子故障診斷方法基本步驟如下:
步驟1確定參數(shù)尋優(yōu)范圍,初始化粒子群各參數(shù)。
Zhang等推薦帶寬閾值0<ξ≤1。筆者經(jīng)多次實(shí)驗(yàn)發(fā)現(xiàn),隨著ξ的取值接近于0,計(jì)算耗時(shí)增加,而且邊際譜出現(xiàn)較多的特征頻率外的成分干擾。一般ξ取值達(dá)到約0.1時(shí)即可獲得較好的分解效果,邊際譜圖更加干凈,計(jì)算速度較快。因此,本文中帶寬閾值的尋優(yōu)范圍為0.1≤ξ≤1,B樣條階數(shù)的尋優(yōu)范圍在5≤n≤30,獲得的參數(shù)組合在實(shí)際工程應(yīng)用中已足夠達(dá)到較好的分解性能。
參考文獻(xiàn)[15-16],并經(jīng)過(guò)多次數(shù)值實(shí)驗(yàn),本文中粒子群的參數(shù)設(shè)置見(jiàn)表1。
表1 粒子群算法參數(shù)設(shè)置Tab.1 Parameter setting of particle swarm optimization
步驟2使用粒子群尋優(yōu)的最佳參數(shù)組合進(jìn)行TVFEMD,獲得若干IMF(Intrinsic Mode Function)。
步驟3對(duì)IMF進(jìn)行HT,根據(jù)瞬時(shí)信息得到故障信號(hào)的Hilbert時(shí)頻圖和邊際譜。
步驟4根據(jù)時(shí)頻信息診斷轉(zhuǎn)子故障類型。
采用如圖1所示的轉(zhuǎn)子實(shí)驗(yàn)臺(tái)驗(yàn)證本文所提方法的有效性和優(yōu)越性。該實(shí)驗(yàn)臺(tái)由圓盤轉(zhuǎn)子、動(dòng)壓軸承、調(diào)速電機(jī)控制器、預(yù)緊力支架、ZonicBook/618E型信號(hào)采集設(shè)備、后處理計(jì)算機(jī)等設(shè)備組成。
在轉(zhuǎn)子圓盤轉(zhuǎn)子45°的標(biāo)記處安裝1 g的不平衡質(zhì)量來(lái)模擬轉(zhuǎn)子不平衡故障。轉(zhuǎn)子的轉(zhuǎn)速為3 000 r/min,采樣頻率為2 560 Hz。分別在水平和垂直兩個(gè)方向安裝傳感器采集信號(hào),從水平方向的信號(hào)中任選連續(xù)的2 000點(diǎn)作為待分析信號(hào)。
圖1 Bently RK4轉(zhuǎn)子實(shí)驗(yàn)臺(tái)和信號(hào)采集器Fig.1 Bently RK4 rotor test rig and signal collector
表2為在人為設(shè)定的不同參數(shù)組合下進(jìn)行TVFEMD所得Hilbert邊際譜部分頻帶的能量比值。
表2 不平衡故障不同參數(shù)組合的Hilbert邊際譜部分頻帶能量比Tab.2 Partial spectrum energy ratio of Hilbert marginal spectrum with different parameter combinations for unbalanced faults
為對(duì)比分析不同參數(shù)取值時(shí)分解效果的差異,同時(shí),驗(yàn)證使用邊際譜部分頻帶能量比指標(biāo)作為適應(yīng)度函數(shù)衡量分解效果的有效性,現(xiàn)選取表2中能量比最大的參數(shù)組合1、居中的任意參數(shù)組合3以及最小的參數(shù)組合8,分別在三組參數(shù)組合下進(jìn)行TVFEMD,得到Hilbert邊際譜如圖2所示。
圖2(a)中可見(jiàn)清晰的基頻及二倍頻成分,且集聚性好,幅值較大;圖2(b)雖然也存在清晰基頻及二倍頻成分,但與圖2(a)相比幅值較低;圖2(c)則完全無(wú)法分離出清晰的基頻及倍頻成分,分解的效果很差。可見(jiàn)隨著邊際譜部分頻帶能量比逐漸降低,其分解效果也隨之變差。由此驗(yàn)證,使用Hilbert邊際譜部分頻帶能量比作為衡量分離效果的指標(biāo)是可行的。
同時(shí),表2也表明不同的參數(shù)組合下進(jìn)行TVFEDM所得Hilbert邊際譜的部分頻帶能量比差距較大,且無(wú)明顯規(guī)律性。因此,為達(dá)到良好的分解效果,適當(dāng)?shù)膮?shù)組合的選取十分必要。而人為選擇參數(shù)存在較大的主觀性和盲目性,使用優(yōu)化算法獲取最佳參數(shù)組合無(wú)疑具有重要的工程意義。
圖3為使用粒子群算法進(jìn)行參數(shù)尋優(yōu)時(shí)Hilbert邊際譜部分頻帶能量比隨進(jìn)化代數(shù)的變化曲線。粒子群進(jìn)化到第5代即得到了信號(hào)最大的Hilbert邊際譜部分頻帶能量比0.011 8。尋優(yōu)所得最佳參數(shù)組合為ξ=0.31,n=13。在該參數(shù)組合下進(jìn)行TVFEMD,結(jié)果如圖4所示。
圖2 不同參數(shù)組合下不平衡故障信號(hào)分解結(jié)果Fig.2 Decomposition results of unbalanced fault signals under different combinations of parameters
圖3 不平衡故障Hilbert 邊際譜部分頻帶能量比隨進(jìn)化代數(shù)的變化Fig.3 Hilbert partial spectrum energy ratio of unbalanced fault changes with evolutionary algebra
圖4 不平衡故障信號(hào)使用本文方法分析所得Hilbert邊際譜和時(shí)頻圖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)相比,基頻的幅值更大,二倍頻成分進(jìn)一步壓縮,分解的效果更好。圖4中的時(shí)頻圖,基頻集聚性強(qiáng),不平衡特征明顯。使用粒子群算法尋優(yōu)的參數(shù)組合進(jìn)行TVFEMD,能夠有效的診斷不平衡故障,實(shí)現(xiàn)良好的分解效果。
為驗(yàn)證本文所提方法相比其他方法的優(yōu)越性,分別使用原始EMD方法、文獻(xiàn)[9]中的指定參數(shù)進(jìn)行TVFEMD、文獻(xiàn)[11]所提參數(shù)尋優(yōu)方法對(duì)不平衡故障信號(hào)進(jìn)行分析,結(jié)果如圖5所示。其中,文獻(xiàn)[9]采用ξ=0.1,n=28的參數(shù)組合;對(duì)實(shí)驗(yàn)采集的轉(zhuǎn)子不平衡故障信號(hào)采用文獻(xiàn)[11]所提基于相關(guān)峭度指標(biāo)進(jìn)行尋優(yōu)所得參數(shù)組合為ξ=0.37,n=20。
圖5 不同方法處理不平衡故障信號(hào)分析結(jié)果Fig.5 Analysis results of unbalanced fault signals processed by different methods
圖5(a)所示采用原始EMD方法分析結(jié)果中,僅能發(fā)現(xiàn)幅值極小的基頻附近成分。圖5(b)所示采用文獻(xiàn)[9]方法分析結(jié)果中,除基頻外還存在幅值較為突出的二倍頻成分。圖5(c)所示采用文獻(xiàn)[11]方法分析結(jié)果與圖4所示采用本文方法進(jìn)行分析的結(jié)果相比,集聚性較差,且邊際譜基頻幅值較低??傊捎闷渌N現(xiàn)有方法,無(wú)法突出故障特征頻率,實(shí)現(xiàn)對(duì)故障類型的準(zhǔn)確判斷,而使用本文方法對(duì)故障信號(hào)進(jìn)行處理故障特征頻率突出,能夠準(zhǔn)確對(duì)不平衡故障進(jìn)行診斷,從而驗(yàn)證了本文所提方法在轉(zhuǎn)子恒定轉(zhuǎn)速不平衡故障診斷中的優(yōu)越性。
為驗(yàn)證本文所提方法轉(zhuǎn)子變轉(zhuǎn)速工況下的適用性,采集實(shí)驗(yàn)臺(tái)從1 400~2 000 r/min升速過(guò)程油膜渦動(dòng)故障的振動(dòng)信號(hào),采樣頻率為1 280 Hz,采樣點(diǎn)數(shù)為3 000。圖6和圖7分別為升速油膜渦動(dòng)故障的時(shí)域波形和頻譜。頻譜中存在嚴(yán)重的頻率模糊現(xiàn)象,無(wú)法準(zhǔn)確提取故障特征頻率。
圖6 升速油膜渦動(dòng)故障的時(shí)域波形Fig.6 Time domain waveform of oil film whirling fault at rising speed
圖7 升速油膜渦動(dòng)故障的頻譜Fig.7 Spectrum of oil film whirling fault at rising speed
圖8為使用本文所提的方法對(duì)故障信號(hào)進(jìn)行分析所得Hilbert邊際譜部分頻帶能量比隨進(jìn)化代數(shù)的變化曲線。其中,對(duì)于變轉(zhuǎn)速轉(zhuǎn)子故障信號(hào),粒子群進(jìn)化到第3代即可得到最大的幅值能量比0.037 45。尋優(yōu)所得最佳的參數(shù)組合ξ=0.12,n=20。在該參數(shù)組合下對(duì)故障信號(hào)進(jìn)行分解,所得Hilbert時(shí)頻圖如圖9(a)所示。圖中能夠清晰提取轉(zhuǎn)頻及其二分頻成分,符合轉(zhuǎn)子油膜渦動(dòng)故障的特征,且集聚性較好,具有良好的分解效果。
圖9(b)為采用原始EMD分析轉(zhuǎn)子故障信號(hào)的Hilbert時(shí)頻圖,譜圖中存在嚴(yán)重的頻率模糊現(xiàn)象,無(wú)法準(zhǔn)確提取故障特征頻率。圖9(c)為采用文獻(xiàn)[9]指定的參數(shù)組合ξ=0.1,n=28進(jìn)行TVFEMD和HT所得的Hilbert時(shí)頻圖。雖然也能較為清晰提取轉(zhuǎn)頻及其二分頻,但特征頻率處能量明顯小于圖9(a)采用本文方法分析所得結(jié)果。采用文獻(xiàn)[11]所提出的相關(guān)峭度指標(biāo)進(jìn)行尋優(yōu),得到最優(yōu)參數(shù)組合ξ=0.42,n=12,在該參數(shù)組合下進(jìn)行TVFEMD和HT,所得Hilbert時(shí)頻圖如圖9(d)所示。雖然能夠準(zhǔn)確提取轉(zhuǎn)子的二分頻,卻無(wú)轉(zhuǎn)頻成分,與頻譜所示的特征頻率不一致,由此造成故障特征信息的丟失??傊?,本文所提方法能夠準(zhǔn)確提取非平穩(wěn)轉(zhuǎn)子故障信號(hào)的特征頻率,相比于其他三種方法具有明顯的優(yōu)越性。
圖8 油膜渦動(dòng)故障Hilbert 邊際譜部分頻帶能量比隨進(jìn)化代數(shù)的變化Fig.8 The variation of Hilbert marginal spectrum and partial band energy ratio of oil film whirling fault with evolution algebra
圖9 升速油膜渦動(dòng)故障信號(hào)不同方法分析Hilbert時(shí)頻圖Fig.9 Analysis of Hilbert time spectrum by different methods of rising oil film whirl fault signal
本文提出了一種基于參數(shù)優(yōu)化時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解和希爾伯特變換診斷轉(zhuǎn)子故障的方法。該方法以Hilbert邊際譜部分頻帶能量比來(lái)定量衡量TVFEMD的分解效果,并作為粒子群尋優(yōu)算法的適應(yīng)度函數(shù)實(shí)現(xiàn)參數(shù)尋優(yōu);通過(guò)將分解所得的IMF進(jìn)行HT得到的瞬時(shí)信息來(lái)診斷轉(zhuǎn)子故障類型。分別將該方法應(yīng)用于恒定轉(zhuǎn)速下的轉(zhuǎn)子不平衡和變轉(zhuǎn)速下的油膜渦動(dòng)兩種典型轉(zhuǎn)子故障,驗(yàn)證了該方法診斷轉(zhuǎn)子故障的有效性和優(yōu)越性,得到結(jié)論如下:
(1)使用Hilbert邊際譜部分頻帶能量比作為衡量TVFEMD分離效果的指標(biāo)是可行的。
(2)使用參數(shù)優(yōu)化時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)合希爾伯特變換對(duì)故障信號(hào)進(jìn)行分析,能夠突出故障特征頻率成分,提高分解效果,從而有效的診斷非平穩(wěn)轉(zhuǎn)子故障。
(3)本文所提參數(shù)優(yōu)化時(shí)變?yōu)V波經(jīng)驗(yàn)?zāi)B(tài)分解方法相比于原始經(jīng)驗(yàn)?zāi)B(tài)分解和現(xiàn)有方法相比能夠克服模式混疊現(xiàn)象,更加準(zhǔn)確提取故障特征頻率,從而有效診斷轉(zhuǎn)子不平衡、油膜渦動(dòng)等典型故障。