吳文壽, 李允公, 王 波, 李國萌, 石悅紅
(東北大學 機械工程與自動化學院,沈陽 110819)
?
基于聽覺模型和極值點概率密度的斷齒故障特征提取方法研究
吳文壽, 李允公, 王 波, 李國萌, 石悅紅
(東北大學 機械工程與自動化學院,沈陽 110819)
齒輪斷齒故障的重要特征是嚙合過程中在斷齒處產生碰撞與沖擊。考慮到人耳聽覺系統對于突發的瞬態聲信號具有本能的反應,為提取斷齒故障誘發的瞬態沖擊響應成分,提出一種基于聽覺模型和信號極值點概率密度的特征提取方法。該方法首先對信號進行GT帶通濾波、相位調整及極值點提取,然后計算各極值點的幅值概率密度,通過對其求導判斷各濾波通道中是否存在瞬態沖擊成分,繼而提取與之相關的極值點。同時,由于系統振動時會產生與斷齒沖擊無關的極值點,為準確提取斷齒沖擊,根據瞬態信號頻帶連續性和多頻段分布特點,設計了相應的提取方法。經實測信號驗證表明,所提方法能準確刻畫及提取斷齒故障特征,可以在含有多種類型的瞬態沖擊響應成分中提取出只由斷齒故障所誘發的沖擊成分,且提取結果精確度較高。
斷齒故障;聽覺模型;概率密度;瞬態信號;故障診斷;特征提取
齒輪傳動是一種最為常見的傳動形式。據統計,在機械設備故障或狀態不良的原因中,約60%[1]與齒輪失效有關,而輪齒斷裂在齒輪失效形式中所占比例又最高,約為41%[2]。斷齒使齒輪副不能正常嚙合,并產生周期性沖擊響應成分。斷齒故障的動力學特性十分豐富[3],但最本質、最能區別于其他齒輪故障特征的在于斷齒會誘發出周期性瞬態沖擊響應成分[4],且斷齒越嚴重,沖擊越明顯。因此,如果能提取斷齒故障的沖擊響應成分,并且能夠識別出沖擊響應出現的周期,則可對齒輪斷齒故障診斷提供有效依據。在瞬態沖擊成分提取方面,栗茂林等[5]提出一種對連續小波系數進行非線性化,實現沖擊成分特征提取。WANG等[6]提出一種基于Levenberg-Marquardt方法的瞬態模型和參數識別迭代提取方法,并最終應用于軸承和齒輪故障特征提取。嚴保康等[7]提出一種利用形態提升方法,通過放大微弱沖擊成分來實現微沖擊特征提取,并成功運用到滾動軸承早期微弱脈沖故障檢測中。近年以來,利用譜峭度[8]方法也受到學者深入研究,主要利用譜峭度對瞬態沖擊成分較為敏感,依此來發現沖擊成分所占據的頻段,以便確定濾波器的帶寬和中心頻率,此方法對于提取瞬態沖擊信號表現良好。
人耳聽覺系統對突發的瞬態聲音信號具有本能的反應[9],基于此,李允公等[10]提出一種基于極值點概率密度和聽覺模型的瞬態信號提取方法,該方法從瞬態信號導致的概率密度曲線的長拖尾現象出發,并結合聽覺模型提出相應的存在瞬態信號判斷方法。此方法提取的瞬態信號起始時間的準確性較高,且當背景信號或干擾噪聲較強時也可獲得滿意效果,但未考慮在對不同類型瞬態信號的區分。
考慮到即使在正常的齒輪系統中也有可能產生瞬態振動成分(如輕載時側隙過大),為提取只與斷齒有關的瞬態沖擊成分,本文基于聽覺模型和信號極值點幅值概率密度,提出了專門針對由斷齒故障所誘發的沖擊成分提取方法。經實驗驗證表明該方法能很好的提取僅與斷齒故障相關的沖擊響應成分,并且能夠識別沖擊響應出現的周期。
1.1 基本原理
本文所提方法的實現過程如圖1所示,首先利用Gammatone濾波器組對信號進行帶通濾波,為消除濾波時出現的相位不一致,對濾波之后的結果再采用逆序濾波,然后計算各濾波通道內幅值極大值點的概率密度及其導數。根據處理結果判斷是否存在瞬態沖擊成分,繼而篩選有效幅值極值點并進行區分提取,最終提取得到只含有斷齒故障的瞬態信號成分。

圖1 方法實現過程原理圖Fig.1 The principle diagram of the method implementation process
1.2 Gammatone濾波器
Gammatone濾波器組[11]是通過模擬人耳耳蝸基底膜的濾波特性所提出的一種帶通濾波方法。設振動信號為x(t),共有N個濾波器,則第i個濾波器的沖擊響應表達式為
h(i,t)=Bntn-1exp(-2πBt)cos(2πfit+φi)
(1)
式中:n為濾波器階數;φi為相位;fi為第i個濾波器的中心頻率。
參數B與fi的關系為:
B=1.019×ERB(fi)=1.019(24.7+0.108fi)
(2)
式中:ERB(fi)為Gammatone濾波器的等效矩形帶寬。
濾波器的時域波形如圖2所示。

圖2 濾波器時域波形Fig.2 The filter waveform in time domain
則帶通濾波為
y(i,t)=x(t)*h(i,t)
(3)
式中:*代表卷積。
為避免同一頻率成分在不同濾波通道中出現相位差異,在基于傅里葉變換進行卷積運算時忽略濾波器h(i,t)的相頻特性,實現逆序濾波,即
(4)
式中:F-1表示傅里葉逆變換,X(f)和H(i,f)分別為x(t)和h(i,t)的頻域表達。
隨后對逆序濾波之后的信號提取所有極大值點,其余幅值均置為零,得到y1(i,t)。
本文選擇200個Gammatone濾波器來實現帶通濾波,濾波器中心頻率通常按照對數均勻分布設置,圖3為不同中心頻率下濾波器的頻譜圖(圖中只繪出16個濾波器)。

圖3 濾波器頻譜圖Fig.3 The spectrum of filter
1.3 瞬態信號識別方法
若機械系統中存在碰撞與沖擊,一般情況下,其振動信號中會存在幅值明顯的瞬態沖擊響應成分,信號的概率密度曲線也必然會出現如圖4所示的長拖尾現象[8]。

圖4 瞬態信號所致概率密度曲線的長拖尾現象Fig.4 Heavy-tailed in probability densitycurve induced by transient signal
但是僅依據長拖尾部分難以判斷瞬態信號的幅值范圍,因此,可提取信號的極大值點,并計算極大值點的概率密度曲線p(i,g)。由于極值點幅值間往往具有較強的不連續性,所以在極值點的幅值概率密度曲線中的長拖尾部分會出現小幅值波動情況,而且這部分的幅值必然與瞬態沖擊成分有關,基于此特點,為實現對瞬態成分的自動識別,計算p(i,g)關于g的導數
(5)
式中:g表示逆序濾波之后信號幅值的極大值點。
利用d(i,g)曲線,按以下準則[10]判斷第i通道里是否含有瞬態沖擊成分:
① 在d(i,g)曲線經歷最小值之后,存在向上過零點,記該點的橫坐標為g0。
②ηq為一個小值,且
(6)
式中:θ是預先給定的閾值,本文取θ=0.1。
1.4 瞬態信號篩選提取
依據上述兩條判斷準則,可初步提取與瞬態沖擊成分相關的極值點,若d(i,g)滿足以上準則,則令:
(7)
式中:幅值為1的位置很可能為瞬態沖擊信號的極值點。
為使初步提取的極值點連續化,并形成波形,可將y2(i,t)和一方波信號做卷積得到y3(i,t),即
y3(i,t)=y2(i,t)*r(t)
(8)
式中:r(t)表示方波信號。
由于信號在經時頻分解之后,瞬態沖擊成分會占據一定的連續頻帶,且在某些頻率處達到峰值。為使沖擊成分辨識更加突出明顯,將y3(i,t)沿各通道方向累加求和得到θ1(t),即
(9)
此外,為明確沖擊成分的頻域分布情況,將y3(i,t)沿著時間方向累加求和得到θ2(i),即
(10)
繼而計算Ω(i,t),即
Ω(i,t)=θ1(t)θ2(i)
(11)
可知,若Ω(i,t)非零,則說明在第i通道、t時刻時存在瞬態沖擊成分。
為了進一步將y3(i,t)中無關的極值點刪除,使其達到簡煉的目的,繼而計算z(i,t),即
z(i,t)=y3(i,t)Ω(i,t)
(12)
1.5 斷齒瞬態信號區分提取
由于輪齒間嚙合拍打、存在側隙或軸承出現小故障,均會引起瞬態沖擊成分,但是其往往高頻部分能量較小,且占據的頻帶較窄,這也意味著會出現與斷齒故障無關的極值點。
而斷齒導致的沖擊響應成分往往能量較大,且會占據較寬的頻段。因此,為將淹沒在其中的斷齒沖擊響應成分準確提取,對z(i,t)進行分頻段統計,平均分為α個頻段依次累加分別得到s1(t)、s2(t)、s3(t),…,sα(t),即
(13)
(14)
(15)
?
(16)
式中:ν=N/α,N為濾波器個數。
將sk(t)(k=1,2,3,…,α)兩兩相乘得到alp(t),即
alp(t)=sl(t)sp(t)(l,p=1,2,3,…,αl≠p)
(17)
按圖5所示的流程原則比較有值區域和無值區域,通過選取sk(t)和alp(t)來獲取提取結果。

圖5 sk(t),alp(t)選擇流程圖Fig.5 Selection flow chart of sk(t) and alp(t)
圖5中
C(t)=count{sk(t)>0,k=1,2,3,…,α}
(18)
表示當在t時刻時,若第k個頻段上滿足sk(t)>0,則計為1,否則為零,各頻段依次將結果累加為C(t)。
繼而用alp(t)分別和含有瞬態沖擊成分的濾波信號相乘之后再累加求和得到Tq,即
(19)
最終,Tq(t)即為提取出的只含有由斷齒所引發的瞬態沖擊響應成分。
為驗證本文所提方法在齒輪斷齒故障特征提取中的有效性,搭建如圖6所示的二級斜齒圓柱齒輪減速器實驗臺及人為斷齒故障,其主要零部件參數為:電動機額定功率750 W,最大轉速1 400 r/min;減速器高低速級齒數分別為13,86,14,85;選用壓電式加速度傳感器拾取振動信號,考慮到沖擊成分會激發各階固有頻率,且其固有頻率往往較高,為滿足采樣定理且提高信號的還原度,設置采樣頻率為10 240 Hz。

圖6 減速器實驗臺和人為齒輪斷齒故障Fig.6 Redactor test rig and artificial gear fracture fault
測得減速器高速軸轉速為937 r/min,理論計算旋轉頻率fr=15.62 Hz,嚙合頻率fm=203.3 Hz,圖7為斷齒故障的振動加速度信號時域波形圖和頻譜圖。在時域波形圖中雖然可以明顯看出存在瞬態沖擊成分,但難以分辨出斷齒故障特征;頻譜圖中雖然存在旋轉頻率和嚙合頻率及其高次倍頻,但譜線過于密集,不能準確判斷出是何故障。

圖7 齒輪斷齒故障信號時域波形及其頻譜Fig.7 Time waveform and frequency spectrum of gear fracture fault
濾波器參數選擇如下:濾波器個數N=200;濾波器階數n=4;相位φi=0;中心頻率fi按照對數均勻分布設置。圖8是經過逆GT之后的時頻圖,從此圖中也難以辨識瞬態沖擊成分。

圖8 逆GT濾波結果Fig.8 Filtering results of the inverse GT

圖9 各通道過零點Fig.9 Zero-crossing of each channel

圖10 第196通道濾波結果Fig.10 The 196th channel’s filtering results
依據給定的兩條判斷準則判斷篩選后,各通道與其對應的過零點g0對應關系如圖9所示,其中幅值為零表示該通道無向上的過零點,即表明無瞬態沖擊成分。
以第196通道為例說明。圖10是振動信號經逆GT濾波之后的濾波結果及局部放大圖,在圖10(b)中所指幅值為9.509。
第196通道的概率密度曲線p(g)及其導數d(g)如圖11和圖12所示,從11(b)中看出p(g)的長拖尾部分出現小幅值波動,且在第三個波峰處幅值為9.469。在12(b)中,g0=9.469是d(196,g)經最小值之后逐漸上升且由負變正的第一個幅值點,并且絕大多數沖擊成分的極值點的幅值都會大于g0點,依此便可確定此通道的瞬態沖擊信號的幅值范圍。

圖11 第196通道概率密度曲線Fig.11 The 196th channel’s probability density curve

圖12 第196通道概率密度曲線的導數Fig.12 The 196th channel’s derivative of probability density
在10 (b)中幅值為9.509的極大值點與圖12 (b)的第一個過零點g0點相對應。繼而根據上述理論依據,計算θ1(t)、θ2(t)如圖13所示,從圖13 (a)中已明顯的看出周期性瞬態沖擊成分,但依然有其他沖擊成分。

圖13 各通道累加結果Fig.13 Accumulative results of each channel
4個不同頻段的sk(t)(k=1,2,3,4)分布情況如圖14所示。由圖可知低頻段s1(t)為零,中低頻段s2(t)已出現明顯的周期性,且有值區域均與高頻段s4(t)相對應。為此,經sk(t)兩兩相乘得到alp(t)如圖15所示,因s1(t)為零,圖中只給出a23(t)、a24(t)、a34(t)的分布情況。
經流程原則選擇篩選后,a24(t)=s2(t)s4(t)即為最終保留的篩選結果,最終提取波形如圖16(a)所示,從圖中可以看出其周期性非常明顯,每兩個沖擊之間的平均時間間隔為Δt=0.422 4 s,經計算得中間軸轉速r=2.360 6 r/s,亦即中間軸大齒輪每轉一圈耗時為0.423 6 s,和提取的瞬態沖擊時間間隔相一致,從而表明系統發生斷齒故障,且發生在中間軸的大齒輪上,與人為制造的斷齒位置相匹配。

圖14 不同頻段sk(t)圖Fig.14 sk(t) diagram in different frequency bands
圖16(b)是提取結果的局部放大圖。從圖中可以看出,在經一次嚙合過程中出現兩次沖擊(嚙入嚙出各一次),對于其余三個嚙合碰撞也出現同樣的結果,并且在圖16(c)所示的實測信號的對應位置波形振動趨勢和提取結果局部放大圖相一致。

圖15 alp(t)圖Fig.15 alp(t) diagram

圖16 最終信號提取結果Fig.16 Final signal extraction result
本文基于聽覺模型和極值點概率密度,提出了一種針對斷齒故障特征的提取方法,試驗驗證結果表明:
(1) 極值點概率密度經求導后,利用其是否存在向上過零點即可自適應的明確各濾波通道中是否存在瞬態成分及其對應的幅值范圍。
(2) 該方法能夠將斷齒故障所誘發的瞬態沖擊成分與其他無關沖擊相區分,只提取與斷齒有關的沖擊響應成分,且提取結果便于人和計算機識別,為設備斷齒故障檢測帶來便捷。
(3) 輪齒嚙合過程中,在斷齒處會產生兩次沖擊,可作為斷齒故障的又一特征。
[1] 唐貴基,龐爾軍,王曉龍. 基于EMD的齒輪箱齒輪故診斷的研究[J]. 機床與液壓,2013,41(13):188-190.
TANG Guiji,PANG Erjun,WANG Xiaolong. Research on gear fault diagnosis based on EMD[J]. Mchine Tool & Hydraulics,2013,41(13):188-190.
[2] 馮偉. 基于振動分析的齒輪斷齒故障研究[J]. 廣州航海:高等專科學校學報,2007,15(1):13-16.
FENG Wei. Investigation on gear fracture failure based on vibration analysis[J]. Journal of Guang Zhou Maritime College,2007,15(1):13-16.
[3] 馬銳,陳予恕. 齒輪傳動系統斷齒故障的機理研究[J]. 振動與沖擊,2013,32(21):47-51.
MA Rui,CHEN Yushu. Fault mechanism of a gear system with tooth broken[J]. Journal of Vibration and Shock,2013,32(21):47-51.
[4] JENA D P,SAHOO S,PANIGRAHI S N. Gear fault diagnosis using active noise cancellation and adaptive wavelet transform[J]. Measurement,2014,47:356-372.
[5] 栗茂林,梁霖,王孫安,等. 基于連續小波系數非線性流形 學習的沖擊特征提取方法[J].振動與沖擊,2012,31(1):106-111.
LI Maolin,LIANG Lin,WANG Sunan,et al. Mechanical impact feature extraction method based on nonlinear manifold learning of continuous wavelet coefficients[J]. Journal of Vibration and Shock,2012,31(1):106-111.
[6] WANG Shibin,CAI Gaigai,ZHU Zhongkui,et al. Transient signal analysis based on Levenberg-Marquardt method for fault feature extraction of rotating machines [J]. Mechanical Systems and Signal Processing,2015,54(55):16-40.
[7] 嚴保康,周鳳星. 一種基于形態提升的自適應軸承微沖擊提 取方法[J]. 振動與沖擊,2013,32(24):198-203.
YAN Baokang,ZHOU Fengxing. Adaptive weak impulse extraction method of rolling bearings based on morphological lifting wavelet[J]. Journal of Vibration and Shock,2013,32(24):198-203.
[8] GUO W,TSE P W,Djordjevich A. Faulty bearing singnal recovery from large noise using a hybrid method based on spectral kurtosis and ensemble empirical mode decomposition [J]. Measurement,2012,45:1308-1322.
[9] KAYA E M,ELHILALI M. A temporal saliency map for modeling auditoryattention [C]//46th Annual Conference on Information Sciences and Systems. Princeton:IEEE,2012:1-6.
[10] 李允公,張金萍,戴麗. 基于極值點概率密度和聽覺模型的瞬態信號提取方法研究[J]. 振動與沖擊,2015,34(21):37-53.
LI Yungong,ZHANG Jinping,DAI Li. A method extracting transient signals based on probability density of extreme points and on auditory model[J]. Journal of Vibration and Shock,2015,34(21):37-53.
[11] LI Yungong,ZHANG Jinping,DAI Li,et al. Auditory-model-based feature extraction method for mechanical faults diagnosis[J]. Chinese Journal of Mechanical Engineering,2010,21(3):391-397.
A method extracting fault features of gear teeth fractures based on an auditory modeland probability density of extreme points
WU Wenshou, LI Yungong, WANG Bo, LI Guomeng, SHI Yuehong
(School of Mechanical Engineering & Automation,Northeastern University,Shenyang 110819,China)
The collision and impact at gear teeth fractures are the most important features in meshing process of gear pairs. Considering a human auditory system has an instinctive response to sudden transient acoustic signals, in order to extract transient impulse response components induced by gear teeth fracture faults, a method extracting fault features of gear teeth fractures based on an auditory model and signal probability density of extreme points was proposed. Firstly, band-pass filtering with Gammatone filters , phase adjustment and extreme points extraction for signals were conducted, and then the amplitude probability densities of extreme points were calculated, their derivatives were used to judge if there are transient impact components in all filtered signals. Those extreme points with transient impact components were extracted. Meanwhile, the whole system vibration might produce extreme points being irrelevant to gear teeth fracture impacts. In order to accurately extract impacts of gear teeth fractures, according to transient signals’ frequency band continuity and multi-band distribution characteristics, an appropriate extraction method was designed. The real measured signals showed that the proposed method can accurately extract fault features of gear teeth fractures, and can extract impact components only induced by gear teeth fracture faults in a variety of transient impulse response components, and the extraction accuracy is higher.
gear teeth fracture fault; auditory model; probability density; transient signals; fault diagnosis; feature extraction
國家自然科學基金資助項目(51275080)
2015-12-31 修改稿收到日期:2016-03-04
吳文壽 男,碩士生,1991年生
李允公 男,博士,副教授,1976年生
TH17
A
10.13465/j.cnki.jvs.2016.19.017