沖擊特征受控極小化通用稀疏表示及其在機(jī)械故障診斷中的應(yīng)用

賀王鵬,訾艷陽,陳彬強(qiáng)

(1.西安交通大學(xué)機(jī)械制造系統(tǒng)工程國家重點(diǎn)實(shí)驗室, 710049, 西安;2.美國紐約大學(xué) 工程學(xué)院, 11201, 美國紐約;3.廈門大學(xué)航空航天學(xué)院, 361005, 福建廈門)

?

沖擊特征受控極小化通用稀疏表示及其在機(jī)械故障診斷中的應(yīng)用

賀王鵬1,2,訾艷陽1,陳彬強(qiáng)3

(1.西安交通大學(xué)機(jī)械制造系統(tǒng)工程國家重點(diǎn)實(shí)驗室, 710049, 西安;2.美國紐約大學(xué) 工程學(xué)院, 11201, 美國紐約;3.廈門大學(xué)航空航天學(xué)院, 361005, 福建廈門)

針對強(qiáng)背景噪聲下機(jī)械故障微弱暫態(tài)特征表示和有效提取的難題,提出了通用的稀疏優(yōu)化特征提取算法。算法針對含噪聲沖擊性特征提取問題設(shè)計了稀疏優(yōu)化表征函數(shù),該函數(shù)融合了沖擊特征的保真度與懲罰函數(shù)因子,考慮了正則化參數(shù)以適應(yīng)不同工程背景下各分析因子的實(shí)際影響,實(shí)現(xiàn)處理結(jié)果稀疏性極大化。同時,引入受控極小化方法對設(shè)計的表征函數(shù)進(jìn)行轉(zhuǎn)化,分解成一系列凸優(yōu)化分析問題。提出了針對離散信號的有限差分式數(shù)值迭代算法,驗證了其快速收斂性和數(shù)值穩(wěn)定性,提出的算法對機(jī)械故障診斷的數(shù)字采樣信號具有普遍適用性。將所提出的算法應(yīng)用于實(shí)驗室環(huán)境下的軸承故障特征識別中,無論是低噪聲還是低信噪比白噪聲環(huán)境下,振動信號中的沖擊特征都得到了顯著增強(qiáng),在Hilbert包絡(luò)譜中的故障特征頻率及其高次諧波比能量中占優(yōu)。所提出的算法還應(yīng)用于電力機(jī)車走行部輪對的故障診斷中,在高強(qiáng)度的工程有色噪聲環(huán)境下精確提取了其中的沖擊衰減成分,在時域和頻域診斷結(jié)果中都得到了準(zhǔn)確的驗證,指導(dǎo)了診斷實(shí)踐。

稀疏表示;凸優(yōu)化問題;受控極小化;特征提取;故障診斷

我國正在大力發(fā)展高端裝備制造產(chǎn)業(yè),高端裝備的核心價值主要體現(xiàn)在其運(yùn)行過程。各行業(yè)中的許多重大關(guān)鍵機(jī)械設(shè)備長期運(yùn)行在重載、腐蝕、高溫等復(fù)雜惡劣的工況下,核心零部件將不可避免地發(fā)生各種程度的機(jī)械故障。信號處理和特征提取技術(shù)是機(jī)械設(shè)備狀態(tài)監(jiān)測和故障診斷的核心技術(shù)。部件故障的早期微弱故障特征提取的一個重要障礙是強(qiáng)背景噪聲。

稀疏表示是近年來信號處理領(lǐng)域的熱點(diǎn)之一[1],對特定信號進(jìn)行稀疏表示的意義就是在給定的超完備字典中用盡可能少的原子來表示該信號,獲得信號更為簡潔的表示方式,以便更容易獲取隱藏在信號中的有用信息[2]。工程中廣泛應(yīng)用的小波變換是基于內(nèi)積原理的特征波形基函數(shù)信號分解,與故障特征最相似或局部最相似的小波基函數(shù)可以最佳地稀疏表征隱藏在混合信號中的故障信息[3-4],從時頻分布特性的差別來說,小波變換可分為Mallat分形、Selesnick分形和Zi分形等種類[5]。目前的降噪方法往往針對信號中的高階可導(dǎo)趨勢項、分段光滑趨勢項,以及沖擊性信號出現(xiàn)時間的時刻特征進(jìn)行了較好的處理,而在強(qiáng)噪聲背景下較完整地提取出能量微弱的振蕩衰減沖擊性特征依然具有很大的挑戰(zhàn)。

本文針對強(qiáng)背景噪聲中微弱故障特征提取的難題,構(gòu)造了新的稀疏優(yōu)化目標(biāo)函數(shù),并且引入受控極小化(MM)方法進(jìn)行推導(dǎo)求解。利用受控極小化方法對構(gòu)造的機(jī)械故障特征提取優(yōu)化問題進(jìn)行求解,同時充分利用帶狀稀疏矩陣的特點(diǎn),提出了具有高計算效率的迭代收斂算法流程[6]。將所提出的稀疏優(yōu)化特征提取算法以滾動軸承的故障診斷為應(yīng)用對象,分別分析了軸承實(shí)驗臺與電力機(jī)車的實(shí)際案例,稀疏有效地提取出了微弱的故障特征,驗證了方法的有效性。

1 問題描述和模型構(gòu)造

1.1 符號定義

本文提出的稀疏優(yōu)化理論是基于l1范數(shù)正則化進(jìn)行構(gòu)建的,l1范數(shù)定義為

(1)

式中:s={s(n)|n=0,1,…,N-1}是定義在RN空間中的一個離散時間序列。序列s的一階有限差分Δ(1)(s,m)和二階有限差分Δ(2)(s,m)分別為

(2)

Δ(2)(s,m)=s(m+1)-2s(m)+s(m-1)

(3)

對于離散時間序列,定義如下兩個差分矩陣

(4)

(5)

1.2 問題描述和模型構(gòu)造

令w、x、y都是RN的離散時間序列,并且滿足

(6)

式中:w為服從高斯分布的白噪聲序列;x為隱藏在背景噪聲中的機(jī)械故障特征;y為觀測到的含噪聲信號序列。

為了從觀測數(shù)據(jù)y中提取稀疏故障特征x,構(gòu)造如下的凸優(yōu)化問題

(7)

式中:F(x)為目標(biāo)函數(shù),由兩部分組成,第1部分為數(shù)據(jù)保真度度量算子,用于度量x={x(n)}和y={y(n)}的差異,第2部分為正則項(懲罰函數(shù));參數(shù)λ是正則化參數(shù),用于調(diào)節(jié)組成目標(biāo)函數(shù)的兩部分在優(yōu)化中所占的權(quán)重,其取值應(yīng)滿足λ>0,取值設(shè)定一般與信號中噪聲的標(biāo)準(zhǔn)差成比例。式(7)中的稀疏矩陣G應(yīng)該根據(jù)待提取的有用特征信號x的特點(diǎn)具體設(shè)定,如x的系數(shù)本身具有稀疏性,則G可取為單位矩陣,此時式(7)處理信號時相當(dāng)于軟閾值收縮降噪[7];如x是逐段恒定的信號,則G可取為式(4)中定義的一階微分矩陣操作算子,此時式(7)就是經(jīng)典的總變分(全變分)降噪方法[8]。

2 受控極小化方法

受控極小化適用于解決那些直接求解比較困難的優(yōu)化問題。受控極小化方法并不直接極小化式(7)中的目標(biāo)函數(shù)F(x),相反地,它求解一系列的優(yōu)化問題,記為Qk(x,xk),k=0,1,2,…。這種思路的出發(fā)點(diǎn)就是求解Qk(x,xk)比直接求解F(x)要簡單。受控極小化方法通過極小化Qk(x,xk)獲得相應(yīng)的序列值xk。

MM方法中要求Qk(x,xk)是凸函數(shù),并且是目標(biāo)函數(shù)F(x)的受控優(yōu)化算子,即滿足如下約束

圖1所示為極小化單變量函數(shù)的MM方法的執(zhí)行過程,圖中x0為初始值。由圖1可觀察到Q0(x,x0)在x0處與F(x)在此處的取值相等,其余值大于F(x)在相應(yīng)處的值。第1次迭代時,通過極小化Q0(x,x0)得到下一次迭代的初始值x1,即Q1(x,x1)在x1處與F(x)在此處的取值相等。通過類似的反復(fù)迭代,最終收斂于目標(biāo)函數(shù)F(x)的最小值。以此類推,可以得到MM方法用于求解多變量函數(shù)的極小化。由于MM方法使用簡單、計算高效,因此在多變量的優(yōu)化求解中其作用更加突出。

圖1 極小化單變量函數(shù)的MM方法的執(zhí)行過程

綜上所述,利用MM方法求解目標(biāo)函數(shù)F(x)可以概括為以下步驟:

(1)令k:=0,初始化x0;

(2)構(gòu)造合適的Qk(x,xk)使其滿足式(8)和式(9)中的約束;

(3)記Qk(x,xk)的極小值為xk+1,即

(10)

(4)令k:=k+1,轉(zhuǎn)步驟(2)。

3 算法推導(dǎo)

采用MM方法求解式(7)中定義的優(yōu)化問題,一種通用的方法是在構(gòu)造便于優(yōu)化的Qk(x,xk)(即受控優(yōu)化算子)時利用二次函數(shù)代替式(7)中的正則項。此時,目標(biāo)函數(shù)F(x)的優(yōu)化問題就轉(zhuǎn)化為二次項優(yōu)化問題,可以通過求解一系列的線性方程組得到解決。

對于x的l1范數(shù),構(gòu)造滿足式(8)和式(9)中的約束條件的函數(shù)如下

(11)

式(11)可以更簡便地表示如下

(12)

式中:Λk為對角矩陣,定義為Λk:=diag(|xk|)。

由Gx代換式(12)中的變量x,得到

(13)

此時,對角矩陣Λk:=diag(|Gxk|)。

利用式(13),可以得到式(7)中目標(biāo)函數(shù)F(x)的Majorizer,即

(14)

式中:對角矩陣Λk:=diag(|Gxk|)。采用受控極小化方法,通過極小化式(12)中的Qk(x,xk)得到xk+1,其解析解為

(15)

式(15)中的對角矩陣Λk里的部分元素可能會隨著Gxk收縮到0而趨于無窮大。因此,需要對式(15)進(jìn)行適當(dāng)?shù)母膶?使其避免產(chǎn)生上述問題。利用矩陣求逆引理[9],式(15)等價于

(16)

式(16)中需要對線性方程組進(jìn)行求解,通常這樣的求解在y點(diǎn)數(shù)較多時計算量會比較大。然而,稀疏矩陣G通常選取為帶狀矩陣,則式(16)中需要求逆的部分也是帶狀矩陣。此時式(16)中的線性系統(tǒng)可以進(jìn)行高效的求解。算法的具體計算步驟歸納如下。

輸入:y∈RN,λ,G;

初始化:x=y;

步驟1:通過式(16)計算xk+1;

步驟2:令xk:=xk+1,轉(zhuǎn)步驟1,直到收斂;

步驟3:輸出x。

4 仿真驗證

為了檢驗本文提出的稀疏優(yōu)化特征提取算法在周期性沖擊故障特征提取中的有效性和實(shí)用性?，F(xiàn)模擬了一組周期性沖擊信號,其中的沖擊單元為

(17)

令采樣頻率fs=1 kHz,信號采樣長度L=2 048。為了模擬工程實(shí)踐中普遍存在的背景噪聲,對原信號加入了高斯白噪聲,信噪比為-3 dB。兩信號的波形如圖2所示。

(a)周期性沖擊信號

(b)含噪聲的觀測信號圖2 仿真信號波形圖

采用本文提出的稀疏優(yōu)化特征提取算法對該含噪信號進(jìn)行分析,其中正則化參數(shù)設(shè)置為λ=0.66,矩陣G選取為式(5)中定義的二階差分矩陣D(2),算法迭代次數(shù)為20,分析結(jié)果及算法收斂特性分別如圖3和圖4所示。在圖3中可以很清晰地發(fā)現(xiàn),以0.2 s(對應(yīng)的頻率為5 Hz)為時間間隔的周期性沖擊單元,與仿真沖擊信號相一致。因此,本文提出的稀疏優(yōu)化特征提取算法可以有效地提取微弱沖擊故障特征,同時對噪聲也有很強(qiáng)的抑制作用。從圖4中可以觀察到,優(yōu)化的目標(biāo)函數(shù)隨著算法迭代次數(shù)的增加而逐步收斂于最優(yōu)解。

圖3 提取出的波形

圖4 目標(biāo)函數(shù)隨迭代次數(shù)的變化

5 工程應(yīng)用

滾動軸承是關(guān)鍵設(shè)備中應(yīng)用最為廣泛的一種通用機(jī)械部件,同時也是一種易損部件,其一旦出現(xiàn)故障,不僅會影響系統(tǒng)的正常運(yùn)行,而且嚴(yán)重時還會導(dǎo)致重大的惡性事故。因此,為了確保設(shè)備運(yùn)行安全,對具有重要用途的軸承進(jìn)行有效的故障診斷是非常必要的[10]。為驗證所提出的稀疏優(yōu)化特征提取方法在機(jī)械微弱故障特征提取中的有效性,本文將其應(yīng)用于實(shí)測的軸承故障振動信號分析中。

5.1 實(shí)驗背景說明

采用美國凱斯西儲大學(xué)軸承數(shù)據(jù)中心的滾動軸承故障模擬實(shí)驗臺數(shù)據(jù)進(jìn)行實(shí)驗驗證分析。實(shí)驗平臺如圖5所示,包括1臺電動機(jī),1個扭矩傳感器/譯碼器,1個功率測試計,還有電子控制器(圖中沒顯示)。被測試軸承支承電機(jī)軸,驅(qū)動端滾動軸承型號為SKF 6205,使用電火花加工技術(shù)在軸承上布置了單點(diǎn)故障。實(shí)驗中轉(zhuǎn)速為1 730 r/min,采樣頻率為12 kHz。表1中列出了軸承各部件的故障頻率倍數(shù),計算得到軸承外圈故障特征頻率為103.4 Hz,內(nèi)圈故障特征頻率為156.1 Hz。

圖5 軸承測試實(shí)驗平臺結(jié)構(gòu)示意圖

故障頻率倍數(shù)內(nèi)圈外圈保持架組滾動體54152358480398347135

圖6所示為在電機(jī)殼體的驅(qū)動端上采集的外圈故障振動信號的時域波形圖。采用本文提出的稀疏優(yōu)化特征提取算法對該含噪信號進(jìn)行分析,設(shè)置正則化參數(shù)λ=0.56,矩陣G選取為單位矩陣,算法迭代次數(shù)為20,分析結(jié)果如圖7所示。在提取得到的特征波形圖中可以很清晰地發(fā)現(xiàn)周期性沖擊單元,與外圈故障特征頻率相一致。

圖6 外圈故障軸承振動加速度信號的時域波形

圖7 外圈故障軸承振動信號特征提取波形

圖8 內(nèi)圈故障軸承振動加速度信號的時域波形

(a)提取出來的波型

(b)包絡(luò)譜圖9 內(nèi)圈故障軸承振動信號特征提取波形及其包絡(luò)譜

圖8所示為在電機(jī)殼體的驅(qū)動端上采集的內(nèi)圈故障振動信號的時域波形圖。采用本文提出的稀疏優(yōu)化特征提取算法進(jìn)行分析,設(shè)置正則化參數(shù)λ=0.46,其余參數(shù)設(shè)置與外圈故障分析一致,分析結(jié)果及其包絡(luò)譜如圖9所示。在圖9a中可以很清晰地發(fā)現(xiàn)周期性沖擊單元與內(nèi)圈故障特征頻率相一致。從圖9b中可以清晰地觀察到存在內(nèi)圈故障特征頻率及其倍頻的譜峰(圖中紅色的點(diǎn)劃線),同時周圍還伴隨著轉(zhuǎn)頻及其諧波成分的邊頻帶(圖中綠色的虛線)。分析結(jié)果與軸承內(nèi)圈存在局部故障時的振動信號特征相一致[11]。

將本文所提出的算法應(yīng)用于同型號同轉(zhuǎn)速下正常滾動軸承的信號分析中,振動信號時域波形如圖10所示,可以觀察到正常軸承振動信號振幅明顯低于故障軸承振動信號。采用本文提出的算法對該振動信號進(jìn)行分析,設(shè)置正則化參數(shù)λ=0.21,矩陣G選取為單位矩陣,算法迭代次數(shù)為20,分析結(jié)果如圖11所示。本文分析結(jié)果顯示,所提出的算法可以將正常軸承振動信號中大量隨機(jī)噪聲濾除,消噪效果較好,同時還不會產(chǎn)生偽故障沖擊特征。

圖10 正常軸承振動加速度信號的時域波形

圖11 正常軸承振動信號特征提取波形

5.2 稀疏分析結(jié)果與對比

將本文提出的算法應(yīng)用于某電力機(jī)車的滾動軸承的外圈輕微損傷故障診斷。軸承置于測試機(jī)車軸承的某實(shí)驗裝置上,軸承型號為52732QT,具體參數(shù)見表2,其中接觸角為0。實(shí)驗中采樣頻率為12.8 kHz,軸承轉(zhuǎn)速為481 r/min,數(shù)據(jù)長度為4 096,通過計算得到軸承外圈故障特征頻率為57.8 Hz。

表2 測試滾動軸承具體參數(shù)

利用加速度傳感器采集的振動數(shù)據(jù)如圖12所示。時域信號中可以觀察到噪聲比較大,看不到明顯的沖擊信號。利用本文提出的算法對該含噪聲信號進(jìn)行分析,設(shè)置正則化參數(shù)λ=0.4,矩陣G選取為式(5)中定義的二階差分矩陣D(2),算法迭代次數(shù)為20,提取特征波形及其包絡(luò)譜如圖13所示。從圖13a中可以很清晰地發(fā)現(xiàn)間隔約為17.1 ms的周期性沖擊單元,其對應(yīng)的頻率為58.5 Hz,與計算得到的外圈故障特征頻率57.8 Hz相一致。圖13b中存在很明顯的外圈故障特征頻率57.8 Hz及其倍頻成分。

圖12 故障軸承振動加速度信號的時域波形

(a)提取出來的波形

(b)包絡(luò)譜圖13 故障軸承振動信號特征提取波形及其包絡(luò)譜

6 結(jié) 論

(1)本文針對機(jī)械故障診斷中微弱故障特征提取的難題,構(gòu)造了有效的稀疏優(yōu)化目標(biāo)函數(shù),該目標(biāo)函數(shù)中的正則項(懲罰函數(shù))可根據(jù)待分析信號特征靈活調(diào)節(jié),實(shí)現(xiàn)分析結(jié)果的最大稀疏化。

(2)引入受控極小化方法進(jìn)行求解所構(gòu)造的目標(biāo)函數(shù)。利用受控極小化方法對構(gòu)造的機(jī)械故障特征提取優(yōu)化問題進(jìn)行求解,同時充分地利用帶狀稀疏矩陣的特點(diǎn),提出了高效的迭代收斂算法,便于在工程中推廣應(yīng)用。

(3)本文提出的稀疏優(yōu)化模型及算法可以拓展到其他變換域,處理變換得到的表征信號的系數(shù),如短時傅里葉變換系數(shù)和小波變換系數(shù);該稀疏優(yōu)化算法也可以作為后處理工具與現(xiàn)有的方法進(jìn)行結(jié)合,如基于小波的收縮降噪法。

(4)將所提出的稀疏優(yōu)化特征提取算法以滾動軸承的故障診斷為應(yīng)用對象,分別分析了軸承實(shí)驗臺與電力機(jī)車的實(shí)際案例,從振動信號中稀疏有效地提取出了微弱故障特征,驗證了算法的有效性。

[1] LEI Y G, LIN J, ZUO M J, et al. Condition monitoring and fault diagnosis of planetary gearboxes: a review [J]. Measurement, 2014, 48: 292-305.

[2] MALLAT S. A wavelet tour of signal processing: the sparse way [M]. New York, USA: Academic Press, 2008: 292-305.

[3] HE Wangpeng, ZI Yanyang, CHEN Binqiang, et al. Automatic fault feature extraction of mechanical anomaly on induction motor bearing using ensemble super-wavelet transform [J]. Mechanical Systems and Signal Processing, 2015, 54: 457-480.

[4] HE Wangpeng, ZI Yanyang, CHEN Binqiang, et al. TunableQ-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis [J]. Science China Technological Sciences, 2013, 56(8): 1956-1965.

[5] CHEN Binqiang, ZHANG Zhousuo, ZI Yanyang, et al. Detecting of transient vibration signatures using an improved fast spatial-spectral ensemble kurtosis kurtogram and its applications to mechanical signature analysis of short duration data from rotating machinery [J]. Mechanical Systems and Signal Processing, 2013, 40(1): 1-37

[6] FIGUEIREDO M A, BIOUCAS-DIAS J M, NOWAK R D. Majorization-minimization algorithms for wavelet-based image restoration [J]. IEEE Transactions on Image Processing, 2007, 16(12): 2980-2991.

[7] AFONSO M V, BIOUCAS-DIAS J M, FIGUEIREDO M A. Fast image recovery using variable splitting and constrained optimization [J]. IEEE Transactions on Image Processing, 2010, 19(9): 2345-2356.

[8] CONDAT L. A direct algorithm for 1D total variation denoising [J]. IEEE Signal Processing Letters, 2013, 20(11): 1054-1057.

[9] TYLAVSKY D J, SOHIE G R L. Generalization of the matrix inversion lemma [J]. Proceedings of the IEEE, 1986, 74(7): 1050-1052.

[10]何正嘉, 訾艷陽, 陳雪峰, 等. 內(nèi)積變換原理與機(jī)械故障診斷 [J]. 振動工程學(xué)報, 2007, 20(5): 528-533. HE Zhengjia, ZI Yanyang, CHEN Xuefeng, et al, Transform principle of inner product for fault diagnosis [J]. Journal of Vibration Engineering, 2007, 20(5): 528-533.

[11]RANDALL R B, JEROME A. Rolling element bearing diagnostic: a tutorial [J]. Mechanical Systems and Signal Processing, 2011, 25(2): 485-520.

(編輯 杜秀杰)

Majorization Minimization Oriented Sparse Optimization Method for Feature Extraction Technique in Machinery Fault Diagnosis

HE Wangpeng1,2,ZI Yanyang1,CHEN Binqiang3

(1. State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an 710049, China; 2. Polytechnic School of Engineering, New York University, New York 11201, USA; 3. School of Aerospace Engineering, Xiamen University, Xiamen, Fujian 361005, China)

To deal with the effective representation and extraction of incipient transient features of mechanical fault, a general sparsity based identification approach is proposed. This approach designs a sparsity optimization function that integrates impulsive feature preserving factor and penalty function factor, and takes the regularization parameter into consideration such as to address the actual factor weights in different situations. The majorization minimization is introduced to simplify the designed function into a series of convex optimization problems. A finite difference based numerical iterative method is developed for the proposed approach, and its fast convergence and numerical stability are illustrated. The proposed approach is versatile to digital signal processing of mechanical fault detecting practices, and is applied to bearing fault identifications in lab. It is shown that no matter in high or low noise backgrounds, the impulsive components are significantly enhanced, which can be verified in the dominant energy ratios of characteristic frequencies in the Hilbert envelope spectrum. This approach is also utilized to conduct bearing fault diagnosis of traction part of electrical locomotive, and the impulsive features masked by heavy colored noises are effectively detected in time domain and spectral domain.

sparse representation; convex optimization; majorization minimization; feature extraction; fault diagnosis

2015-08-11。 作者簡介:賀王鵬(1989—),男,博士生;訾艷陽(通信作者),男,教授,博士生導(dǎo)師。 基金項目:國家自然科學(xué)基金資助項目(51275384);福建省重大科技項目(2014H6025);福建省高端裝備制造協(xié)同創(chuàng)新中心資助項目。

時間:2015-12-30

10.7652/xjtuxb201604015

TH17

A

0253-987X(2016)04-0094-06

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20151230.1814.004.html