基于多尺度形態濾波和遞歸求差的沖擊特征自適應分離方法

2024-03-19 07:18湯明軍

振動與沖擊 2024年5期

和丹,權偉,湯明軍,劉暉

(1.西安工程大學機電工程學院,西安 710048; 2.西安市現代智能紡織裝備重點實驗室,西安 710600;3.西安交通大學機械工程學院,西安 710049)

周期性沖擊特征提取和辨識是旋轉機械故障診斷中的關鍵環節。但在實際工業現場,由于振動傳輸路徑和多振動源耦合等因素影響,測得的振動信號中不僅包含與故障相關的弱周期性沖擊,還包含信號傳遞中強背景噪聲[1]。因此,從復雜振動信號中準確提取和分離出周期性沖擊特征對機械設備故障診斷具有極其重要的意義。

針對周期性沖擊特征提取和分離,國內外研究者開展了大量研究工作,可分為三類:第一類是利用周期和沖擊特性,結合數字濾波理論構造相應濾波器。如:最大相關峭度解卷積(maximum correlation kurtosis deconvolution,MCKD)[2]、多點最優最小熵解卷積(multipoint optimal minimum entropy deconvolution,MOMEDA)[3]、最大二階循環平穩盲解卷積(maximum second order cyclostationary blind deconvolution,CYCBD)[4]和譜峭度[5]等。這類方法通常能夠最大程度強化沖擊特征,但容易受隨機沖擊干擾影響,且在復合故障診斷時需要設計多個濾波器進行,適合使用場景為無隨機沖擊干擾影響,共振帶和故障頻率未知的故障沖擊識別[6]。第二類是利用周期性沖擊波形特征和內積理論設計相應濾波器,如:匹配追蹤[7],拉布拉斯小波[8]等。這類方法通常需要設計與故障匹配的原子庫,且計算量大,優點是能最大限度還原故障信號,適合故障特征已經情況下的故障信號提取[9]。第三類是利用通用型分解方法開展沖擊特征提取研究。如:經驗模態分解(empirical mode decomposition,EMD)[10]、局域均值分解(local mean decomposition,LMD)[11]、變分模態分解(variational modal decomposition,VMD)[12]等。該類方法可將沖擊特征分解到多個模式分量,通常需要結合其他指標進一步篩選故障特征,優點是不限于沖擊特征分離,適合故障未知及多種故障類型的故障信號分離[13]。上述方法雖然在故障沖擊特征提取中取得一定效果,但適用范圍不盡相同。

形態濾波是一種非線性濾波方法,因其具有較強的非線性信號處理和信號形態細節特征保持能力成為故障診斷中研究熱點。該方法利用結構元素對信號進行濾波,在保留特征信息的同時消除噪聲成分,其效果主要依賴結構元素類型,結構元素尺度和形態學算子等。近些年研究主要關注如何自適應選擇最優結構元素以及如何設計與結構元素匹配的形態學算子兩方面。Li等[14]針對結構元素尺度對形態濾波結果不確定性問題,依據沖擊信號特征,構造自適應結構元素,試驗結果表明自適應結構元素較直線型,三角形元素在特征提取和計算效率方面效果更優。Chen等[15]為了提高形態濾波算法計算性能,將具有消除噪聲的對角切片譜(diagonal slice spectrum,DSS)算法與形態濾波器結合,有效提升了特征提取能力。Duan等[16]針對隨機沖擊干擾下的弱沖擊特征提取問題,將形態學濾波融入最小熵解卷積方法(minimum entropy deconvolution,MED)方法,成功消除路徑中隨機沖擊干擾,準確提取出微弱故障特征。Wang等[17]為了提高形態濾波算法性能,利用多尺度思想構造多尺度增強最優梯度積濾波(enhanced optimal gradient product filter,EOGPF)器,并結合特征頻率強度系數(characteristic frequency intensity coefficient,CFIC)指標實現了最優功率譜的選取。Guo等[18]為了提高形態濾波算法計算性能,利用迭代思想設計一種迭代形態差分積小波流程,有效提高了故障診斷精度。上述研究雖然在微弱沖擊特征提取和抗隨機沖擊干擾方面取得一定成果,然而多故障耦合中周期性沖擊特征的準確提取和分離并未完全解決。

考慮到形態學濾波具有較好的泛化能力,本文將開展基于多尺度形態學濾波的沖擊特征自適應提取和分離技術研究,旨在實現弱故障提取和復合故障分離。首先,篩選適合沖擊特征分離的形態學濾波算子,并分析其濾波特性;然后,利用形態學算子、加權諧噪比、SOSO(strengthen operate subtract operate)特征增強理論和迭代求差思想構造一種多尺度故障特征提取和分離方法;最后,利用仿真和試驗信號驗證該方法有效性。對比CYCBD和譜峭度方法,本文所提出的方法無需預設參數,且不受隨機沖擊影響,不僅對單一弱沖擊提取有效,而且可以實現復合沖擊有效分離。

1 基本原理

1.1 形態學濾波

形態學濾波由于其計算效率高并能處理非線性信號,被廣泛應用于旋轉機械的故障診斷上?；镜男螒B學算子是由腐蝕、膨脹、開、閉算子組成。假設原始信號f(n)是一維離散信號定義為F=(0,1,…,N-1),g(m)也是一維離散信號定義為G=(0,1,…,M-1),且(N≥M),g(m)是結構元素,則數學形態學的四種基本運算定義如下:

膨脹算子 (f⊕g)(n)=max{f(n-m)+g(m)}

{1≤n≤N;1≤m≤M}

(1)

腐蝕算子 (fΘg)(n)=min{f(n+m)-g(m)}

{1≤n≤N;1≤m≤M}

(2)

開運算 (f°g)(n)=(fΘg⊕g)(n)

(3)

閉運算 (f·g)(n)=(f⊕gΘg)(n)

(4)

式中,⊕Θ°·分別表示膨脹、腐蝕、開、閉運算。

通過四種基本的形態學算子的級聯組合可以構造出復雜的組合形態學算子,如下:

開閉算子FOC=(f°g·g)(n)

(5)

閉開算子FCO=(f·g°g)(n)

(6)

閉-膨脹算子FCD=(f·g⊕g)(n)

(7)

膨脹-閉算子FDC=(f⊕g·g)(n)

(8)

腐蝕-開算子FEO=(fΘg°g)(n)

(9)

開-腐蝕算子FOE=(f°gΘg)(n)

(10)

形態梯度算子(morphological gradient,MG)

MG=(f⊕g)(n)-(fΘg)(n)

(11)

形態差分算子(difference filter,DIF)

DIF=(f·g)(n)-(f°g)(n)

(12)

閉-開平均算子(average operator,AVG)

AVG=((f·g)(n)+(f°g)(n))/2

(13)

平均組合差分算子(average combination difference,ACDIF)

ACDIF=(FCD(f(n))+FDC(f(n)))/2-

FEO(f(n))

(14)

組合形態濾波-帽變換算子(combination morphological filter-hat transform,CMFH)

CMFH=f(n)-(FCO(f(n))+FOC(f(n)))/2

(15)

1.2 CMFH特征分析

在四種基本形態學算子中,膨脹算子和閉運算能夠保留信號的正脈沖,抑制負脈沖,相反,腐蝕算子和開運算保留負脈沖,移除正脈沖[19]。AVG為實現正負脈沖特征同步提取,將其進行級聯組合得到FOC、FCO、FCD、FDC、FEO、FOE算子,但該類算子會導致信號產生偏移。為了解決偏移問題,前人又提出改進算子MG、DIF、AVG、ACDIF和CMFH[20]。

為篩選出最優沖擊特征提取算子,利用本文仿真信號2、西儲大學軸承外圈故障信號(#311,采樣頻率為12 kHz,損傷尺寸為0.014英寸,負載為2馬力,轉速為1 750 r/min)和渥太華軸承內圈故障信號[21](#I-A-1,采樣頻率為20 kHz,采樣時長10 s,運行轉速從12.5 Hz增加到27.8 Hz)分析形態學算子(AVG、ACDIF、DIF、MG、CMFH)沖擊特征提取和分離特性。引入能量幅值(energy amplitude,EA)作為沖擊特征提取能力評價標準,EA值越大,沖擊特征提取能力越強。公式如下:

(16)

圖1分別展示了三組不同信號下形態學算子(AVG、ACDIF、DIF、MG、CMFH)在尺度為1～60時EA值的變化情況,圖1(a)、圖1(b)和圖1(c)中CMFH算子的EA值在各個尺度下都遠大于其他算子,表明CMFH算子對沖擊特征提取能力更佳。為掌握CMFH算子濾波特性,對其進行幅頻響應分析,如圖2所示。

(a)

圖2 CMFH幅頻響應

從圖2可以看出,隨著尺度的減小,截止頻率逐漸增大,即保留了更多的高頻分量。這表明CMFH濾波器可以實現將故障特征按照尺度變化規律進行逐級分離。

1.3 SOSO增強技術

為了更大限度地提高故障信號中沖擊特征強度,利用SOSO[22]增強技術提高整體算法去噪性能。具體流程,如圖3所示。

圖3 SOSO強化信號特征流程圖

圖3中:Sk(t)為第k次濾波信號;S0(t)為測試信號,F(·)為濾波算法;k為循環濾波次數。上述增強算法可表示為SOSO(Sk(t))。在本文中,濾波算法采用多尺度形態濾波,k取20。具體見本文第2章流程圖步驟1～3。

2 沖擊特征自適應分離方法

本文結合多尺度形態濾波器、SOSO強化理論以及迭代求差思想提出了沖擊信號分離法,可實現故障信號中沖擊特征自適應分離和提取。流程如圖4所示。

圖4 本文方法分析流程圖

具體實現過程如下:

步驟1輸入測試信號fk(t),同時初始化參數k=1;

步驟2利用結構元素(structuring elements,SE)、CMFH算子和λmax構造多尺度形態濾波器。直線型結構元素可以避免信號幅值修改和減少運算量,故選用直線型結構元素[23]。依據式(17)確定濾波器尺度λ范圍為1～λmax

(17)

式中:τmax為測試信號fk(t)自相關函數中除零位置外的最大極值位置;fs為測試信號fk(t)的采樣頻率。

步驟3利用多尺度形態濾波器對測試信號fk(t)進行處理得到多組濾波信號f(k,1)(t)、f(k,2)(t)、…、f(k,λmax)(t),并結合峭度指標(kurtosis,KR)和包絡諧噪比[24](envelope harmonic to noise ratio,EHNR)構造加權諧噪比(weighted harmonic to noise ratio,WHNR)指標,該指標能夠有效指示沖擊特征和沖擊特征周期性。利用WHNR指標作為濾波信號擇優篩選標準,得到最佳沖擊信號分量IMF(k,est)(t),WHNR指標計算公式如下

(18)

式中:μ為濾波信號f(k,λ)(t)的平均值;σ為濾波信號f(k,λ)(t)的標準差;E(f(k,λ)(t))是對濾波信號f(k,λ)(t)包絡處理的時間序列;rf為E(f(k,λ)(t))的自相關函數;τmax是rf的最大位置;rf(τmax)為最大位置對應的幅值;rf(0)為包絡總能量。

步驟4利用SOSO增強技術結合多尺度形態濾波算法對得到的沖擊信號IMF(k,est)(t)進行特征增強,進而得到沖擊信號MIMFk(t);

步驟5判斷沖擊信號MIMFk(t)中ER(energy ratio)是否小于能量比閾值ε,如果是,則得到z組沖擊信號分量,其中z≤k,完成分離。如果否,則利用遞歸求差實現測試信號分離,遞歸求差公式為fk+1(t)=fk(t)-MIMFk(t),fk(t)為測試信號第k次循環后信號。循環終止條件為ER<ε,ER見式(19)。經過大量試驗信號驗證,閾值ε取0.001。

(19)

式中:EMIMFk(t)為MIMFk(t)信號幅值平方和;Efk(t)為測試信號幅值平方和。

3 仿真信號分析

3.1 滾動軸承內外圈復合故障仿真分析

為了驗證本文所提方法在提取和分離沖擊信號方面的優越性,模擬滾動軸承內外圈復合故障仿真信號,并添加諧波干擾、隨機沖擊以及強背景噪聲進行分析。具體仿真信號描述為

(20)

式中:f(t)為故障仿真信號1;y1(t)為軸承外圈發生故障的模擬信號,其中fo為軸承外圈故障頻率;y2(t)為軸承內圈發生故障的模擬信號,其中fi為軸承內圈故障頻率,fn為系統固有頻率;y3(t)為諧波干擾信號,兩個諧波頻率fh1、fh2分別為30 Hz、50 Hz;y4(t)為隨機沖擊干擾信號,其中隨機沖擊的發生時刻分別為tr1=0.1 s,tr2=0.2 s,tr3=0.4 s,δ(t)為信噪比為-4 dB高斯白噪聲。仿真信號采樣頻率fs=10 240 Hz,仿真時長為1 s,其參數具體數值如表1所示。

表1 仿真信號參數

滾動軸承內外圈復合故障仿真分析選用仿真信號1。仿真信號成分如圖5所示。圖5(a)、圖5(b)分別為軸承外圈沖擊和軸承內圈沖擊信號,圖5(c)為諧波干擾信號,圖5(d)為隨機沖擊干擾信號。圖6(a)、圖6(b)分別為復合故障仿真信號1時域圖和包絡譜。在時域波形中可以看出故障沖擊信號已經被噪聲完全掩蓋,在包絡譜中發現諧波耦合頻率(fh2-fh1=20 Hz)和內圈故障特征頻率基頻fi,但軸承外圈故障特征頻率fo在干擾信號下無法識別,表明包絡譜分析方法失效。

(a) 外圈沖擊特征

(a) 時域圖

利用不添加SOSO增強的文中方法對仿真信號1進行分析,結果如圖7所示。得到兩組信號模式分量。其中MIMF1中包絡譜發現了外圈故障基頻fo、內圈故障fi和2fi,MIMF2分量中可以找到內圈故障基頻fi、內圈故障fo和2fo,可以看出該方法未能有效分離出內外圈故障信號。利用本文方法對仿真信號1進行分析,結果如圖8所示。得到3組故障信號模式分量,從MIMF1時域圖中可以看到清晰的沖擊特征,包絡譜上可以清晰地找到軸承內圈故障特征頻率基頻及其倍頻,且譜線突出;同樣,從MIMF2信號包絡譜中可以找到外圈故障特征頻率基頻及其倍頻;MIMF3信號的時域圖出現了三組明顯的高幅值隨機沖擊信號,包絡譜中沒有發現明顯的故障信息。對比圖7和圖8可知,包含SOSO增強算法的本文方法可以有效增強內外圈故障的局部特征,并去除隨機沖擊、諧波以及高斯白噪聲干擾,成功分離出軸承復合故障。這驗證了本文方法在復合故障分離中優勢。

(a)

3.2 滾動軸承單一故障仿真分析

為了進一步驗證本文方法在提取單一沖擊故障特征方面的有效性,選用仿真信號1中的y1(t)、y2(t)、y3(t)和δ(t)分量構造仿真信號2,如圖9所示。在圖9時域波形中可以看出,故障沖擊信號已被噪聲完全掩蓋,包絡譜中無法有效識別出軸承外圈故障特征頻率fo及其倍頻。

(a) 時域圖

利用不添加SOSO增強的文中方法對仿真信號2進行分析,結果如圖10所示。從圖10中MIMF1分量包絡譜可以發現軸承外圈故障頻率fo以及2～4倍頻;MIMF2分量中也可以找到外圈故障基頻fo,可以看出該方法可以實現軸承外圈故障信號提取。利用本文方法對仿真信號2進行分析,結果如圖11所示。MIMF1分量時域圖中呈現明顯的沖擊成分,包絡譜中可以看出軸承外圈故障特征頻率基頻及其2～8倍頻,可以判斷軸承發生了外圈故障。MIMF2和MIMF3信號時域圖和包絡譜沒有發現明顯的故障特征信息。這表明本文方法可以有效提取單一沖擊故障特征。對比圖10和圖11可知,SOSO增強方法有效去除了噪聲對外圈故障特征提取的干擾,增強了軸承外圈沖擊故障特征的局部特征,使得故障特征頻率在譜線上更加明顯。

(a)

4 試驗數據驗證

以YQ-365型異步牽引電機中圓柱滾子軸承為研究對象,振動數據由圖12中的傳感器b獲得。圓柱滾子軸承安裝在電機內部的主軸上。電機結構簡圖如圖13所示。

圖12 YQ-365異步牽引電機

圖13 電機內部結構

4.1 圓柱滾子軸承內圈故障

在進行異步牽引電機的滾子承內圈故障試驗中,采用線切割方式在輸出軸軸承(型號:NU214 C4)上加工內圈故障尺寸為4×20 mm的細槽,用于模擬內圈微小故障。以12.8 kHz采樣頻率采集振動信號。圖14為軸承內圈故障的實物圖。電機轉速為600 r/min,根據軸承型號軸承故障頻率計算公式得到軸承內圈故障特征頻率fi為80.56 Hz。圖15(a)為傳感器所采集到的振動波形,直接對該信號進行包絡譜分析如圖15(b)所示,從包絡譜可以找到轉動頻率fr及其二倍頻,但無法判別是否存在軸承故障。

圖14 軸承內圈損壞

(a) 時域圖

采用本文所提方法對實測信號1進行分析,結果如圖16所示,從圖中可以看出,MIMF1信號時域圖中隨機沖擊成分被成功剔除,強化了沖擊性特征,包絡譜中找出軸承內圈故障特征頻率fi和轉頻fr,以及轉頻fr的2～6倍頻。MIMF2信號包絡譜中只找到轉動頻率的基頻fr?？梢耘卸ㄔ撦S承內圈發生故障。

(a)

利用CYCBD[25]、匹配追蹤[26]、自適應噪聲完全集合經驗模態分解[27](complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)三類不同沖擊特征提取方法對實測信號1進行處理,并與本文方法處理結果對比分析。三個對比方法關鍵參數設置如下:

(1) 利用粒子群優化算法預先確定CYCBD算法參數,最終將CYCBD濾波器的長度和特征頻率值分別設置為312和80 Hz。

(2) 分析實測信號1,利用Gabor解析字典,設置字典參數范圍分別為:波形尺度參數Sf為[0.1∶0.005∶5],時移參數Uf為[0.1∶0.01∶10],頻率參數βf為[1∶1∶6 400]。

(3) 利用CEEMDAN算法對實測信號1分解產生14個IMF分量,通過各模式能量含量和白噪聲特性確定最佳的濾波信號。

三個對比方法分析結果如圖17～圖19所示。圖17展示CYCBD方法分析結果,在時域圖中可以看到兩個十分明顯的隨機沖擊特征,但包絡譜中未能找到有關內圈故障的特征成分;圖18展示匹配追蹤分析結果,時域圖中既包含隨機沖擊信號也包含周期性的沖擊特征,包絡譜出現轉頻fr及其2倍頻,如果想去除隨機沖擊,需要進一步優化原子庫和匹配規則;圖19展示CEEMDAN算法對實測信號1處理結果,從時域圖中可以找到明顯的沖擊特征,隨機沖擊干擾得到有效去除,但包絡譜只能發現轉頻fr及其2倍頻,無法判別軸承是否出現故障;通過對三種方法結果分析得出,CYCBD和匹配追蹤方法受隨機沖擊干擾較大,CEEMDAN分解方法消除隨機沖擊干擾影響,但受背景噪聲干擾較大,無法有效提取出軸承故障特征,相比之下,本文方法提取故障特征效果更佳。

(a) 時域圖

4.2 圓柱滾子軸承內圈和滾動體故障

同樣在YQ-365異步牽引電機上進行圓柱滾子軸承內圈和滾動體復合故障試驗,采用線切割方式在輸出軸承(型號:NU214 C4)上加工滾動體故障尺寸為20×0.25 mm,內圈故障尺寸為8×0.25 mm的細槽,用于模擬內圈和滾動體故障。實物如圖20所示。采樣頻率為12.8 kHz,數據長度為102 400個點,電機轉速為600 r/min。軸承故障特征頻率如表2所示。

表2 軸承特征參數

(a) 內圈損壞

本次試驗采集得到實測信號2,其時域波形和包絡譜如圖21所示。由于受到背景噪聲干擾,時域圖中無法找到較明顯的故障沖擊,包絡譜中只發現保持架的故障頻率fc及其二倍頻,但未發現滾子軸承內圈和滾動體故障特征。表明包絡譜分析方法失效。

(a) 時域圖

應用本文所提方法對實測信號2進行分析得到圖22所示結果。在MIMF1時域信號中沖擊特征明顯,噪聲成分有所降低,包絡譜中可以找到軸承內圈故障特征頻率fi和調制頻率fr,以及頻率fr的2～4倍頻。MIMF2信號中背景噪聲得到有效去除,包絡譜中發現滾動體故障特征頻率fb及其2倍頻,以及調制頻率fc及其2、3倍頻。在MIMF3信號包絡譜只找到保持架基頻fc及其2倍頻。根據分析結果可判斷出軸承滾動體和內圈同時發生故障,診斷結果與實際情況完全相符。

(a)

為驗證本文方法優越性,利用快速譜峭度[28]、匹配追蹤[29]和雙約束非負矩陣分解[30]三類不同沖擊特征提取方法對實測信號2進行分析,圖23為快速譜峭度分析結果。取圖23中2個較明顯的共振頻帶(白色實線框),共振頻帶中心頻率分別取1 860 Hz和5 980 Hz,對應帶寬分別為332 Hz、286 Hz。利用共振解調分別得到2個共振帶進行包絡譜分析,結果如圖24所示。從譜圖中可以看出該方法可以提取出轉動頻率,但并未找到軸承的故障特征頻率。與圖22結果比較,快速譜峭度方法僅僅增強了軸承轉頻fr,不能識別出軸承故障特征。利用匹配追蹤方法對實測信號2進行分離處理,得到兩組分離信號如圖25所示。其中分離信號1時域圖中出現明顯的沖擊特征,包絡譜中可以識別出滾動體故障特征頻率fb、調制頻率fc及其2、3倍頻,分離信號2包絡譜中僅找到調制頻率fc及其3倍頻。與圖22結果比較分析得知,匹配追蹤方法在分析實測信號2中的復合故障時效果不佳,可以提取出滾動體故障特征,但不能識別出軸承內圈故障特征。利用雙約束非負矩陣分解方法對實測信號2進行處理,結果如圖26所示。分離信號1包絡譜中可以找到調制頻率fc的二倍頻,分離信號2包絡譜中也同樣只找到調制頻率fc的二倍頻,分析可知雙約束非負矩陣分解方法在分析實測信號2中的復合故障時效果不佳。經過上述分析可知,本文方法在處理復合故障特征分離問題上有自身優勢。