萬書亭，張 雄，龐 彬，豆龍江

華北電力大學(xué) 機(jī)械工程系，河北 保定 071003)

0 引言

滾動軸承是機(jī)械設(shè)備中應(yīng)用最廣泛的零部件之一，其運(yùn)行狀況影響著整個系統(tǒng)的工作狀態(tài)，滾動軸承的故障識別與診斷對于保證機(jī)械設(shè)備的安全可靠運(yùn)行具有重要意義。實(shí)際工況下，滾動軸承受運(yùn)行狀態(tài)及工作環(huán)境影響，其故障沖擊信號常常與多種噪源干擾融合，且存在轉(zhuǎn)頻調(diào)制現(xiàn)象。從調(diào)制信號中將故障沖擊信號解調(diào)出來一直是機(jī)械故障領(lǐng)域的難點(diǎn)問題。近年來，故障沖擊特征提取方法被廣泛研究，從傳統(tǒng)的遞歸式經(jīng)驗?zāi)B(tài)分解(EMD)，到局部均值分解，再到變分模態(tài)分解(VMD)，以及各類分解方法的改進(jìn)形式[1-5]，信號的分解能力及自適應(yīng)性逐漸增強(qiáng)，對沖擊成分的提取愈加有效。Antoni等系統(tǒng)地定義了譜峭度，先后提出了基于短時傅里葉變換的譜峭度法及譜峭度圖的快速算法，并將其應(yīng)用到旋轉(zhuǎn)機(jī)械的監(jiān)測和故障診斷中。最大相關(guān)峭度解卷積[6-10]是一種可以在低信噪比的情況下，通過提升信號峭度值達(dá)到降噪目的，從而提取微弱沖擊成分的方法。無論是濾波方法、信號分解方法還是解卷積方法，其對信號共振頻帶的依賴均很大，且其通過抑制噪聲干擾增強(qiáng)沖擊特征，并未針對性地對工頻干擾采取措施，而實(shí)際工況中故障沖擊特征往往被工頻干擾所淹沒。

陷波具有快速濾除工頻載波信息的能力，與其他濾波器不同，陷波的參數(shù)設(shè)計取決于工頻轉(zhuǎn)速信息，而不依賴于共振頻帶。但是由于陷波具有窄帶濾波特性，其對中心頻率及帶寬參數(shù)變化較為敏感。而人工參數(shù)選擇的主觀性較大，且多參數(shù)選取工作量大。因此，本文提出一種粒子群參數(shù)優(yōu)化方案以有效地解決陷波過程中中心頻率及帶寬的參數(shù)選擇問題。

1 基本理論

1.1 自適應(yīng)陷波器

1.1.1 陷波原理

在信號分析中，信號會被工頻干擾所淹沒，通過窄帶濾波可以將這些工頻干擾濾除。陷波的目的就是通過構(gòu)成窄帶濾波消除單個的正弦干擾信號，且其他頻率成分都能通過。

對于二階數(shù)字陷波，假設(shè)要濾除的頻率為ω0，濾波器的幅度特性在ω=ω0處為0，在其他頻率上接近常數(shù)。零點(diǎn)z=e±jω0使陷波器的幅度特性在ω=±ω0處為0。為使離開零點(diǎn)后，幅度迅速上升到某個常數(shù)，將2個極點(diǎn)放在非?？拷泓c(diǎn)的地方，極點(diǎn)為z=re±jω0，傳遞函數(shù)為：

(1)

其中，r為半徑，與陷波器的帶寬有關(guān)，0≤r≤1；b0為增益系數(shù)。設(shè)陷波器在 -3 dB處的頻帶帶寬定義為Bw，則Bw和增益系數(shù)b0為：

(2)

當(dāng)使用歸一化頻率時，Bw=1-r。不同的r和帶寬Bw，對應(yīng)不同的 -3 dB處的頻帶和陷波器深度，如圖1所示。

圖1 陷波器幅值響應(yīng)Fig.1 Amplitude response of trap filter

當(dāng)r較小而Bw較大時，ω=±ω0附近的頻率分量會受到顯著的影響；而當(dāng)r→1時，陷波器帶寬Bw減小，同時陷波深度也減小。所以設(shè)計陷波器時需要根據(jù)實(shí)際要求選擇r值。

1.1.2 粒子群優(yōu)化陷波器參數(shù)

粒子群優(yōu)化算法[11]具有良好的全局尋優(yōu)能力，因此采用粒子群優(yōu)化算法實(shí)現(xiàn)陷波中心頻率及帶寬這2個參數(shù)的自適應(yīng)選取。假設(shè)D維空間中，種群X=(X1，X2，…，XM)包含M個粒子，其中第i個粒子Xi=(xi1，xi2，…，xiD)表示1個D維向量，代表第i個粒子在D維搜索空間中的位置。第i個粒子的速度為Vi=(vi1，vi2，…，viD)，Pi=(pi1，pi2，…，piD)為個體局部均值，G=(g1，g2，…，gD)為種群全局極值，各粒子通過Pi和G迭代更新自身速度和位置：

(3)

其中，ω為慣性權(quán)重；d=1，2，…，D；i=1，2，…，M；α為當(dāng)前迭代次數(shù)；c1和c2為加速度因子；η為[0，1]范圍內(nèi)的隨機(jī)數(shù)。

利用粒子群優(yōu)化算法尋優(yōu)時，需要確定一個適應(yīng)度函數(shù)，粒子每次更新位置時需要計算當(dāng)前位置對應(yīng)的函數(shù)值，通過對比進(jìn)行更新。信號中沖擊成分的比重影響峭度指標(biāo)的大小，信號包含的沖擊成分越多[12]，峭度值越大，因此本文將進(jìn)行陷波處理后信號的峭度值K作為適應(yīng)度函數(shù)，表達(dá)式如下：

(4)

其中，x為振動信號；μ為信號x的均值；E(·)表示取均值運(yùn)算；δ為信號x的標(biāo)準(zhǔn)差。

利用粒子群參數(shù)尋優(yōu)的具體步驟如下：

a. 初始化粒子群優(yōu)化算法的各項參數(shù)，參數(shù)值如表1所示；

b. 初始化粒子種群，以參數(shù)組合[ω0，Bw]作為粒子的位置，隨機(jī)初始化各粒子的位置與移動速度；

c. 計算每個粒子的位置對應(yīng)的適應(yīng)度值;

d. 對比適應(yīng)度值的大小并更新種群全局極值；

e. 更新粒子的位置和速度；

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

表1 粒子群優(yōu)化算法各項參數(shù)Table 1 Parameters of particle swarm optimization algorithm

1.2 自適應(yīng)陷波流程圖

自適應(yīng)陷波流程圖如圖2所示。

圖2 自適應(yīng)陷波流程圖Fig.2 Flowchart of adaptive trap

2 仿真信號分析

為驗證自適應(yīng)陷波方法能夠有效從調(diào)制信號中濾除工頻干擾，還原故障沖擊特征，采用文獻(xiàn)[13]中的滾動軸承內(nèi)圈故障模型進(jìn)行分析，數(shù)學(xué)模型如式(5)所示。

(5)

其中，τi為第i次沖擊相對于平均周期T的微小波動；Ai為以1/fr為周期的幅值調(diào)制，fr為軸承所在工作軸的轉(zhuǎn)頻，fr=20 Hz；CA=0；A0=2；h(t)為指數(shù)衰減脈沖；B為系統(tǒng)的衰減系數(shù)；內(nèi)圈故障通過頻率為130 Hz；fn=3 kHz為系統(tǒng)固有頻率；φA、φω為初相位。設(shè)置采樣頻率fs=8 192 Hz，取8 192點(diǎn)數(shù)據(jù)分析。沖擊信號的時域波形如圖3所示。

圖3 沖擊信號時域波形Fig.3 Time-domain waveform of impulse signal

在原始故障仿真模型中混入fc=20 Hz的工頻載波信號C(t)=5 sin(2πfct)，調(diào)制信號的時域波形如圖4所示，信號中原有的沖擊特征被工頻載波信號淹沒。

圖4 調(diào)制信號時域波形Fig.4 Time-domain waveform of modulating signal

陷波器的參數(shù)設(shè)計取決于工頻轉(zhuǎn)速信息，而不依賴于共振頻帶，但是陷波器具有窄帶濾波特性，所以對中心頻率及帶寬參數(shù)的變化較為敏感?？紤]到人工參數(shù)選擇的主觀性較大，且多參數(shù)選取工作量大，本文提出自適應(yīng)陷波方法，利用粒子群多參數(shù)尋優(yōu)能力，根據(jù)最大峭度準(zhǔn)則自適應(yīng)選取陷波參數(shù)。

初始化參數(shù)，設(shè)定中心頻率搜索范圍為[19，21] Hz，帶寬系數(shù)搜索范圍為[0，1]。圖5為峭度值隨進(jìn)化代數(shù)變化的曲線，峭度值8.66出現(xiàn)在第15代進(jìn)化種群中，此時陷波中心頻率ω0=20 Hz、帶寬Bw=0.2。所設(shè)計的陷波器幅值響應(yīng)曲線如圖6所示，工頻陷波后，其波形如圖7所示，從圖7中可以看出，調(diào)制信號的工頻信息被濾除，沖擊特性被有效解調(diào)出來。

圖5 峭度值隨進(jìn)化代數(shù)變化的曲線Fig.5 Value of kurtosis varying with evolutionary algebra

圖6 陷波器幅值響應(yīng)圖Fig.6 Amplitude response of trap filter

圖7 工頻陷波后的時域波形Fig.7 Time-domain waveform of power frequency trap

為驗證粒子群優(yōu)化參數(shù)選擇的合理性，將中心頻率及帶寬修改為[20 Hz，0.1]、[20 Hz，0.6]、[21 Hz，0.2]3種組合，在3種組合參數(shù)下陷波的時域波形見圖8。對比圖3、圖7、圖8可看出：圖8(a)陷波后波形邊界翹曲嚴(yán)重；圖8(b)陷波過度，沖擊能量衰減；圖8(c)中心線波動，存在微弱正弦工頻影響。

圖8 不同陷波參數(shù)下的時域波形Fig.8 Time-domain waveforms of different trap parameters

進(jìn)一步驗證，以波形匹配方差S2為評價指標(biāo)表征不同參數(shù)下陷波后的沖擊特征與原始信號的匹配程度，波形匹配方差數(shù)值越小則波形匹配程度越高，其結(jié)果見表2。

(6)

其中，yi為陷波后信號序列；xi為原信號序列；n為采樣點(diǎn)數(shù)。

表2 不同參數(shù)陷波后的波形匹配方差Table 2 Values of S2 under different trap parameters

3 實(shí)驗信號分析驗證

3.1 微弱故障實(shí)驗分析

為了驗證本文方法對軸承輕微故障診斷的有效性，采用美國Case Western Reserve大學(xué)的滾動軸承數(shù)據(jù)，軸承型號JEMSKF6023-2RS。故障源是滾動體表面通過電火花加工的直徑分別為0.177 8 mm、0.355 6 mm、0.533 4 mm的凹坑。選用最輕微的0.177 8 mm故障數(shù)據(jù)進(jìn)行分析，采樣頻率為12 kHz，軸的轉(zhuǎn)速為1 730 r/min。滾動軸承參數(shù)尺寸如表3所示，滾動軸承發(fā)生外圈故障、內(nèi)圈故障和滾動體故障時的故障特征頻率分別為88、143、115 Hz。

表3 滾動軸承參數(shù)Table 3 Parameters of rolling bearing

取8 192點(diǎn)數(shù)據(jù)進(jìn)行分析，圖9為滾動軸承故障信號的時域波形，直接對其進(jìn)行包絡(luò)分析，結(jié)果如圖10所示。由圖10可見，包絡(luò)譜中成分復(fù)雜，滾動軸承的滾動體微弱故障信息被噪聲頻率淹沒。

圖9 滾動軸承故障信號時域波形Fig.9 Time-domain waveform of rolling bearing’s fault signal

圖10 滾動軸承故障信號包絡(luò)譜Fig.10 Spectrum envelope of rolling bearing’s fault signal

利用粒子群多參數(shù)尋優(yōu)特性設(shè)計自適應(yīng)陷波器，初始化參數(shù)，根據(jù)軸轉(zhuǎn)速(1 730 r/min)設(shè)置中心頻率搜索范圍為[26，30]Hz、帶寬范圍為[0，1]，根據(jù)最大峭度準(zhǔn)則確定陷波中心頻率及帶寬。圖11為峭度值隨進(jìn)化代數(shù)變化的曲線，最大峭度出現(xiàn)在第16代進(jìn)化種群中，此時陷波中心頻率ω0=29.3 Hz(此時實(shí)驗轉(zhuǎn)頻為28.8Hz，考慮為實(shí)驗轉(zhuǎn)速設(shè)置誤差引起)，帶寬Bw=0.6。所設(shè)計的陷波器幅值響應(yīng)曲線如圖12所示。工頻陷波后，其波形如圖13所示，從圖中可以看出，調(diào)制信號工頻信息被濾除，沖擊特性被解調(diào)出來。參數(shù)選擇不當(dāng)(ω0=28 Hz，Bw=0.1)時的陷波后信號的包絡(luò)譜如圖14所示。對比圖13、14可知，當(dāng)參數(shù)選擇不當(dāng)時，陷波未得到理想效果。

圖11 峭度值隨進(jìn)化代數(shù)變化的關(guān)系曲線Fig.11 Value of kurtosis varying with evolutionary algebra

圖12 陷波器幅值響應(yīng)圖Fig.12 Amplitude response of trap filter

圖13 陷波后信號的包絡(luò)譜(ω0=26 Hz，Bw=0.2)Fig.13 Spectrum envelope of signal after trap，when ω0=26 Hz and Bw=0.2

圖14 陷波后信號的包絡(luò)譜(ω0=28 Hz，Bw=0.1)Fig.14 Spectrum envelope of signal after trap，when ω0=28 Hz and Bw=0.1

以陷波作為前置工頻降噪處理方法，結(jié)合集合經(jīng)驗?zāi)B(tài)分解(EEMD)及VMD方法對滾動軸承故障信號進(jìn)行分析。

EEMD算法是Wu[14]等在EMD算法的基礎(chǔ)上，利用高斯白噪聲頻率均勻分布統(tǒng)計特性，在原信號中混入適當(dāng)標(biāo)準(zhǔn)差和疊加次數(shù)的高斯白噪聲，以消除原信號間斷現(xiàn)象的一種自適應(yīng)分解方法；VMD[15]方法可以將1個信號分解成為預(yù)定個數(shù)的準(zhǔn)正交子模態(tài)，每個模態(tài)分量具有頻率中心和有限帶寬，將信號分解問題轉(zhuǎn)換到變分模型中進(jìn)行分解，在各分量之和等于原始輸入信號的約束下，尋求各分量的聚集帶寬之和最小，其實(shí)質(zhì)是多個自適應(yīng)維納濾波器組，具有類似于小波包的頻率分割特性。

對原始信號進(jìn)行EEMD，通過相關(guān)系數(shù)準(zhǔn)則[16]選取最大的2組本征模態(tài)函數(shù)(IMF)分量進(jìn)行信號重構(gòu)，對重構(gòu)信號進(jìn)行包絡(luò)譜分析，結(jié)果如圖15(a)所示。從圖15(a)中可以看出，滾動軸承故障特征頻率115.7 Hz被提取出來，但其譜線并不突出，相較于圖13所示的陷波處理結(jié)果，其干擾譜線較多，且倍頻成分仍未能準(zhǔn)確提取到。進(jìn)而對原始信號進(jìn)行工頻陷波，對陷波后信號進(jìn)行EEMD后重構(gòu)降噪，結(jié)果如圖15(b)所示。從圖15(b)中可以看出，故障特征譜線較突出，干擾譜線得到有效抑制，且能尋找到2倍頻(230 Hz)。

圖15 滾動軸承故障信號的EEMD和陷波-EEMD包絡(luò)譜Fig.15 Spectrum envelopes of rolling bearling’s fault signal after EEMD and EEMD-trap

同理，以陷波為原始前置處理，以VMD后重構(gòu)為后置降噪處理(VMD的分解個數(shù)K采用中心頻率相離準(zhǔn)則[17]選取，懲罰因子默認(rèn)為2 000)，VMD信號包絡(luò)譜及陷波-VMD包絡(luò)譜如圖16所示，可見陷波對于增強(qiáng)沖擊特征有一定的意義。

圖16 滾動軸承故障信號的VMD和陷波-VMD包絡(luò)譜Fig.16 Spectrum envelopes of rolling bearing’s fault signal after VMD and VMD-trap

3.2 復(fù)合故障試驗分析

為了進(jìn)一步驗證自適應(yīng)陷波方法對軸承故障特征頻率信息提取的有效性，在QPZZ軸承故障模擬試驗臺上完成了滾動軸承故障模擬實(shí)驗，試驗平臺的結(jié)構(gòu)見附錄中的圖A1。試驗采用6205E軸承，使用線切割機(jī)在滾動軸承內(nèi)圈上加工出深1.5 mm、寬0.2 mm的凹槽模擬軸承外圈、內(nèi)圈故障，如附錄中的圖A2所示。采用數(shù)據(jù)采集卡由安裝在軸承座上的加速度傳感器采集振動信號，采樣頻率為12 800 Hz，電機(jī)轉(zhuǎn)速為1 468 r/min，滾動軸承參數(shù)尺寸如表4所示。

表4 滾動軸承參數(shù)Table 4 Parameters of rolling bearing

軸承內(nèi)圈故障特征頻率fi、外圈故障特征頻率fo分別為：

(7)

其中，Z為滾動體個數(shù)；D為節(jié)圓直徑；d為滾珠直徑；α為接觸角；N為電機(jī)轉(zhuǎn)速。

圖17 時域波形Fig.17 Time-domain waveforms

圖18 不同陷波參數(shù)下的時域波形Fig.18 Time-domain waveforms under different trap parameters

對時域信號進(jìn)行采集，取8 192點(diǎn)數(shù)據(jù)分析，未處理前的轉(zhuǎn)頻調(diào)制波形及自適應(yīng)陷波(ω0=50 Hz，Bw=0.2)后的沖擊特征時域波形如圖17所示，由圖可見，自適應(yīng)陷波方法可以將正弦載波信息濾除，解調(diào)出間隔明顯的沖擊成分。圖18為不同陷波參數(shù)下的時域波形。通過對比圖17、18可以發(fā)現(xiàn)，陷波參數(shù)選擇不當(dāng)將導(dǎo)致陷波過度，沖擊能量降低(ω0=50 Hz、Bw=0.6情況下)，亦或是陷波不徹底，殘留趨勢項(ω0=56 Hz、Bw=0.2情況下)。

進(jìn)一步進(jìn)行包絡(luò)譜分析，分析結(jié)果如圖19所示。由圖19可見，利用陷波可以在原始信號包絡(luò)譜圖中提取到外圈故障特征頻率fo及其倍頻2fo分別為87.5 Hz、176.6 Hz(由于計算機(jī)進(jìn)行包絡(luò)譜分析時存在頻譜分辨率問題，與理論計算公式計算所得的結(jié)果會有誤差)，較微弱的內(nèi)圈故障信息被淹沒。自適應(yīng)陷波后的沖擊信號的包絡(luò)譜可以進(jìn)一步還原內(nèi)圈故障的沖擊特性fi為131.3 Hz(由于計算機(jī)進(jìn)行包絡(luò)譜分析時存在頻譜分辨率問題，與理論計算公式計算所得的結(jié)果會有誤差)。可見，工頻干擾嚴(yán)重的環(huán)境下，陷波對沖擊特征增強(qiáng)效果明顯，尤其在復(fù)合故障源情況下，可以避免包絡(luò)解調(diào)不足造成的故障情況誤判。若進(jìn)一步針對外圈故障特征頻率進(jìn)行陷波處理，可以從復(fù)合故障譜線中提取到明顯的內(nèi)圈故障頻率及其倍頻，如圖20所示。

圖19 滾動軸承內(nèi)、外圈復(fù)合故障信號包絡(luò)譜Fig.19 Spectrum envelopes of outer and inner ring compound fault signal of rolling bearing

圖20 內(nèi)圈沖擊特性包絡(luò)譜Fig.20 Spectrum envelope of inner impact characteristic

4 結(jié)論

針對軸承故障診斷中調(diào)制信息易被工頻載波信號淹沒的問題，將工頻陷波理論引入到軸承故障診斷領(lǐng)域。陷波消除工頻干擾的參數(shù)選取只取決于轉(zhuǎn)速信息，不依賴于共振頻帶選擇，具有很好的自適應(yīng)性。通過粒子群優(yōu)化算法優(yōu)化陷波理論中的中心頻率及帶寬選取方法，去除了人為因素影響，提高了陷波效率及準(zhǔn)確性。仿真分析結(jié)果表明，利用陷波可以從工頻調(diào)制信號中解調(diào)出原始沖擊信息，對于工頻干擾較大的工況具有明顯的工頻去噪效果。實(shí)驗驗證結(jié)果表明，利用工頻陷波輔助傳統(tǒng)信號分解降噪算法，能夠準(zhǔn)確提取淹沒于噪聲干擾中的沖擊信息。

附錄見本刊網(wǎng)絡(luò)版(http:∥www.epae.cn)。