王 超， 孔凡讓， 胡 飛， 劉 方

(中國科學(xué)技術(shù)大學(xué)精密機(jī)械與精密儀器系， 安徽 合肥 230022)

引 言

隨著現(xiàn)代鐵路運(yùn)輸?shù)目焖侔l(fā)展和不斷提速，其安全性問題變得日益突出。軸承是火車上必不可少的重要零部件之一，而大多數(shù)旋轉(zhuǎn)機(jī)械的失效都是由軸承故障和缺陷引起的[1～3],因此建立健全火車軸承的狀態(tài)監(jiān)測(cè)與故障診斷系統(tǒng)變得尤為重要。20世紀(jì)80年代，列車聲音檢測(cè)系統(tǒng)(ADBD)技術(shù)的發(fā)展在預(yù)報(bào)和診斷軸承失效和過熱方面取得良好的效果[4]。ADBD通過在道旁設(shè)置固定麥克風(fēng)陣列，當(dāng)火車經(jīng)過時(shí)麥克風(fēng)陣列采集聲音信號(hào)并傳輸至信號(hào)處理系統(tǒng)，由此實(shí)現(xiàn)火車軸承的監(jiān)測(cè)與診斷，其相比于機(jī)載監(jiān)視系統(tǒng)(OBM)[5]和熱軸承檢測(cè)系統(tǒng)(HBD)[6]有著經(jīng)濟(jì)實(shí)用、適應(yīng)性強(qiáng)等優(yōu)點(diǎn)[3]。

然而當(dāng)列車高速運(yùn)行經(jīng)過ADBD時(shí)，其診斷的有效性將有所降低，其中一個(gè)重要的原因就是列車與麥克風(fēng)陣列之間的高速相對(duì)運(yùn)動(dòng)引起的多普勒效應(yīng)，多普勒效應(yīng)會(huì)在時(shí)域和頻域上對(duì)采集到的信號(hào)產(chǎn)生較大的畸變和頻移，這無疑將嚴(yán)重影響診斷效果[7]。因此為了提高列車高速運(yùn)行中麥克風(fēng)采集到信號(hào)的可用性，進(jìn)而提高列車軸承狀態(tài)監(jiān)測(cè)和故障診斷的有效性，多普勒效應(yīng)的消除勢(shì)在必行。

近年來，國內(nèi)外許多學(xué)者針對(duì)多普勒效應(yīng)的消除做了大量研究并取得相應(yīng)的成果。楊殿閣等建立了基于測(cè)量面、輻射面和聲全息面的運(yùn)動(dòng)學(xué)幾何關(guān)系，提出了聲源與測(cè)量信號(hào)之間的非線性時(shí)間映射方法，從而消除了多普勒效應(yīng)[8]。但這種方法需要預(yù)先測(cè)量傳感器與鐵軌的垂直距離以及列車速度等參量，使得其在實(shí)際使用中受到了限制；Dybala提出了一種基于希爾伯特變換的面向干擾的動(dòng)態(tài)信號(hào)重采樣方法以消除道旁監(jiān)測(cè)系統(tǒng)受到的多普勒效應(yīng)的影響[7,9]，然而這種方法在頻域處理時(shí)只能包含單一的頻率譜線，這與實(shí)際使用環(huán)境中存在的寬帶頻率成分相矛盾；此外希爾伯特變換的端點(diǎn)效應(yīng)也對(duì)此方法的有效性產(chǎn)生嚴(yán)重的影響。

針對(duì)上述兩種方法存在的問題，提出一種基于短時(shí)傅里葉變換和Crazy Climber算法(STFT-CC)的瞬時(shí)頻率估計(jì)方法，應(yīng)用此方法得到擬合的瞬時(shí)頻率曲線，結(jié)合莫爾斯理論得到信號(hào)重采樣所需的各項(xiàng)參數(shù)，并針對(duì)多普勒效應(yīng)進(jìn)行運(yùn)動(dòng)學(xué)建模，確定時(shí)域重采樣時(shí)間間隔，進(jìn)而消除信號(hào)中多普勒效應(yīng)的影響。此方法能夠有效提取重采樣所需各項(xiàng)參數(shù)，從而避免了對(duì)未知量的預(yù)測(cè)量。

文中首先介紹了聲學(xué)運(yùn)動(dòng)學(xué)模型、時(shí)域重采樣方法和基于STFT-CC的瞬時(shí)頻率估計(jì)方法；然后針對(duì)提出的方法進(jìn)行了信號(hào)仿真；最后將此方法應(yīng)用于實(shí)際故障軸承的診斷中。仿真和實(shí)驗(yàn)結(jié)果均表明本文方法的可行性與有效性。

1 診斷方法和原理

1.1 診斷方法

消除列車軸承多普勒效應(yīng)并進(jìn)行故障診斷的方法主要分為五個(gè)步驟，即基于STFT的時(shí)頻分析，基于CC算法瞬時(shí)頻率估計(jì)，對(duì)估計(jì)得到離散時(shí)頻點(diǎn)進(jìn)行非線性擬合，基于時(shí)域重采樣技術(shù)消除多普勒效應(yīng)，包絡(luò)分析得出診斷結(jié)果。圖1給出了本文提出方法的流程圖。首先對(duì)多普勒畸變信號(hào)進(jìn)行STFT時(shí)頻分析，從時(shí)頻譜圖上觀察瞬時(shí)頻率變化，估計(jì)CC算法閾值，然后利用CC算法對(duì)時(shí)頻譜圖上的各條瞬時(shí)頻率脊線進(jìn)行估計(jì)，并根據(jù)聲學(xué)理論進(jìn)行非線性最小二乘法擬合，得到連續(xù)的瞬時(shí)頻率曲線，對(duì)擬合得到的連續(xù)瞬時(shí)頻率做時(shí)域重采樣，即可對(duì)原信號(hào)中的多普勒效應(yīng)進(jìn)行消除。最后，對(duì)去除多普勒效應(yīng)的信號(hào)進(jìn)行包絡(luò)分析，分析包絡(luò)譜圖后可得到故障診斷結(jié)果。

圖1 消除信號(hào)多普勒效應(yīng)流程圖

1.2 時(shí)域重采樣原理

消除多普勒效應(yīng)最常用的方法是用不同的采樣頻率對(duì)原信號(hào)進(jìn)行重采樣，進(jìn)而還原出不失真的信號(hào)，其中的關(guān)鍵就是確定重采樣的時(shí)間序列。最近，Dybala提出了一種基于希爾伯特變換的面向干擾的動(dòng)態(tài)信號(hào)重采樣方法以消除多普勒效應(yīng)[7]，這種方法需要已知被測(cè)對(duì)象的特征頻率，而且其在頻域處理時(shí)只包含有單一譜線。時(shí)域重采樣方法與之相比具有明顯的優(yōu)勢(shì)，由于其完全基于聲學(xué)運(yùn)動(dòng)學(xué)模型，因此時(shí)域重采樣方法更加易于理解和應(yīng)用。

列車經(jīng)過麥克風(fēng)時(shí)的情況如圖2所示，這里假定列車以恒定速度沿筆直鐵軌運(yùn)行，僅考慮單一聲源和單一麥克風(fēng)的情況。

圖2 列車經(jīng)過麥克風(fēng)示意圖

由圖2可知，在聲源發(fā)聲時(shí)刻與麥克風(fēng)接收時(shí)刻之間存在一段延時(shí)，隨著聲源和麥克風(fēng)之間距離的縮短，其延時(shí)時(shí)間將逐漸變短，這意味著信號(hào)在時(shí)間軸上被壓縮了，因此信號(hào)頻率提高了，反之當(dāng)聲源遠(yuǎn)離麥克風(fēng)時(shí)，信號(hào)在時(shí)間軸上被拉伸，因此信號(hào)頻率降低，這就是運(yùn)動(dòng)聲源測(cè)量中多普勒效應(yīng)的產(chǎn)生原因。

圖2中點(diǎn)A即初始位置，此處時(shí)間為0。假設(shè)在時(shí)間為ts時(shí)，聲源發(fā)出的信號(hào)幅值為x(ts)，此信號(hào)經(jīng)過傳播在時(shí)間tr時(shí)到達(dá)麥克風(fēng)，因此有下式成立

tr=ts+Δt

(1)

式中 Δt即聲音信號(hào)在空氣中的傳播時(shí)間，可由下式得到

(2)

式中c為空氣中的聲速(設(shè)為c=340 m/s)。R是聲源與麥克風(fēng)在發(fā)聲時(shí)刻ts時(shí)的距離，將式(2)代入式(1)可得

(3)

式(3)建立了發(fā)聲時(shí)刻ts和麥克風(fēng)接收采樣時(shí)刻tr之間的關(guān)系，因此重采樣時(shí)間序列可通過式(3)的計(jì)算得到，式中參量Y,S和v可從1.3節(jié)中得到。

重采樣的過程通過對(duì)原信號(hào)進(jìn)行三次樣條插值完成，如圖3所示。圖中tm為麥克風(fēng)實(shí)際采樣時(shí)間序列，tp為重采樣時(shí)間序列，xp為原信號(hào)幅值序列，xr為校正信號(hào)的幅值序列。

圖3 重采樣過程示意圖

時(shí)域重采樣的具體步驟如下：

(1)以聲源在A點(diǎn)出發(fā)時(shí)刻為0時(shí)刻，校正麥克風(fēng)采集信號(hào)時(shí)間原點(diǎn)，即tm=tr+R0/c，其中R0為出發(fā)時(shí)刻聲源與麥克風(fēng)之間的距離，c為聲速。

(2)由式(3)計(jì)算得到tp，將tm代入式(3)中的ts，得到的tr即為tp。

(3)對(duì)離散的xp(tm)進(jìn)行三次樣條擬合，對(duì)擬合后的曲線利用tp插值。

(4)所得信號(hào)xr(tp)即為還原信號(hào)。

1.3 基于STFT-CC的瞬時(shí)頻率估計(jì)

定義一個(gè)頻率調(diào)制信號(hào)為

x(t)=Aexp(jφ(t))

(4)

則信號(hào)的瞬時(shí)頻率定義為信號(hào)相位的一階導(dǎo)數(shù)ω(t)=φ′(t)。信號(hào)的瞬時(shí)頻率是一個(gè)非常重要的參量，其估計(jì)的準(zhǔn)確程度關(guān)系到對(duì)其他信號(hào)參量的求取[10]。

自從Ville于1948年提出瞬時(shí)頻率的概念以來，發(fā)展出許多計(jì)算和估計(jì)信號(hào)瞬時(shí)頻率的方法，例如希爾伯特變換、分?jǐn)?shù)階傅里葉變換和Prony法等[11～14]，然而這些方法大都因頻率估計(jì)效果不理想而在實(shí)際使用中受到限制。而基于STFT的瞬時(shí)頻率估計(jì)以其線性時(shí)頻變換、無交叉項(xiàng)等優(yōu)點(diǎn)得到廣泛應(yīng)用，因此這里使用STFT獲得信號(hào)的時(shí)頻分布。此外，基于STFT的時(shí)頻分布在信號(hào)瞬時(shí)頻率變化較慢的情況下能夠在后續(xù)的信號(hào)處理過程中表現(xiàn)出較為滿意的效果。

獲取信號(hào)瞬時(shí)頻率的一般方法是求取信號(hào)時(shí)頻分布的譜峰值，然而該方法不能同時(shí)提取多條瞬時(shí)頻率曲線，且在高噪聲和復(fù)雜工況環(huán)境下，時(shí)頻分布譜峰檢測(cè)法并不能給出理想的瞬時(shí)頻率估計(jì)[15]。

Crazy Climber 算法是由Carmona 于1999年提出的[16]，該算法在一定程度上能夠同時(shí)提取多個(gè)瞬時(shí)頻率。在時(shí)頻譜圖上，信號(hào)的能量都集中在一些被稱為脊線的地方，通過隨機(jī)分布在該時(shí)頻譜圖上的質(zhì)點(diǎn)受到脊線的吸引而建立的馬爾可夫鏈。在提取脊的過程中，該算法將一系列按一定規(guī)則隨機(jī)運(yùn)動(dòng)的點(diǎn)視為一種密度分布，所有的點(diǎn)在經(jīng)過一定程度的移動(dòng)以后形成整體上的密度分布圖，然后將密度最為集中的區(qū)域指定為脊線的軌跡。

Crazy Climber 算法提取多脊線的核心思想是這樣的[17]：在起始時(shí)刻時(shí)，一系列的質(zhì)點(diǎn)(可稱之為Climbers)隨機(jī)的安排在時(shí)頻分布圖上，每一個(gè)質(zhì)點(diǎn)可認(rèn)為是分布密度的度量；然后每一個(gè)質(zhì)點(diǎn)在時(shí)頻面上按照一定的規(guī)則進(jìn)行移動(dòng)，而這個(gè)規(guī)則取決于該質(zhì)點(diǎn)以及該點(diǎn)下一時(shí)刻可能在時(shí)頻分布的非零矩陣M的值；經(jīng)過移動(dòng)后，質(zhì)點(diǎn)逐漸被時(shí)頻分布圖上脊所在的位置吸引而逐漸聚集，就像在爬山一樣。類似于模擬退火方法[18]，整個(gè)系統(tǒng)也是有一個(gè)設(shè)定的初始溫度，并且溫度逐漸降低，隨之各個(gè)可以移動(dòng)的點(diǎn)逐漸穩(wěn)定下來，聚集在脊所在的位置上。

假設(shè)時(shí)頻分布矩陣M的大小是B×U，M(j,k)表示在時(shí)頻分布上位置為(j,k)的值，其中j=1,2,…,B為水平方向，k=1,2,…,U為豎直方向。所有質(zhì)點(diǎn)獨(dú)立地依據(jù)規(guī)則移動(dòng)，該規(guī)則描述如下：

在初始時(shí)刻t=0，質(zhì)點(diǎn)總數(shù)設(shè)為N，它們平均分布在網(wǎng)格上Γ={1,2,…,B}×{1,2,…,U}。質(zhì)點(diǎn)的初始位置用Xα(0)表示，其中α=1,2,…,N作為質(zhì)點(diǎn)的標(biāo)記。另外，初始溫度設(shè)為T(0)。

在某時(shí)刻t時(shí)，某一個(gè)質(zhì)點(diǎn)停留在位置上Xα(t)=(j,k)，那么它的下一個(gè)位置Xα(t+1)=(j′,k′)由兩個(gè)條件確定。首先是水平方向，如果2≤j≤B-1，那么j′=j-1或者j'=j+1各占50%的概率，即該點(diǎn)以相同的概率向左或向右移動(dòng)一格；如果該質(zhì)點(diǎn)剛好處于時(shí)頻分布的邊際位置，那么該質(zhì)點(diǎn)須向邊際的反方向移動(dòng)，即如果j=1，質(zhì)點(diǎn)向右移動(dòng)一格j′=j+1，反之如果j=B，質(zhì)點(diǎn)向左移動(dòng)一格j′=j-1；因此在水平方向移動(dòng)后，該質(zhì)點(diǎn)的位置變?yōu)?j′,k)。然后是豎直方向，與水平方向一樣質(zhì)量向上或向下移動(dòng)的概率相同(邊際情況處理與水平方向相同)，但是它也可能不做移動(dòng)而停留原位置，如果M(j′,k′)>M(j′,k)，質(zhì)點(diǎn)須移動(dòng)，即Xα(t+1)=(j′,k′)；反之，這個(gè)移動(dòng)按照概率

(5)

執(zhí)行，按照概率(1-pt)不執(zhí)行。同時(shí)，系統(tǒng)溫度將更新為T(t+1)。

當(dāng)系統(tǒng)溫度的值低于某個(gè)設(shè)定的閾值時(shí)，這個(gè)迭代移動(dòng)過程結(jié)束。

當(dāng)整個(gè)迭代過程結(jié)束后，需要對(duì)網(wǎng)格上的每個(gè)點(diǎn)進(jìn)行度量，度量的方式存在兩種。第一種稱之為均值度量，在某一時(shí)刻，考慮每個(gè)移動(dòng)點(diǎn)相當(dāng)于質(zhì)量為1/N的質(zhì)點(diǎn)，則可以度量該網(wǎng)格點(diǎn)在t時(shí)刻的度量值為

(6)

第二種度量方式稱之為加權(quán)度量，該方式考慮到了將原始的時(shí)頻分布平面上的數(shù)據(jù)用來對(duì)上式的結(jié)果進(jìn)行加權(quán)，期望不僅能提取到脊所在的正確的位置，而且脊上參數(shù)的變化能更加貼近原來的時(shí)頻平面上參數(shù)的變化規(guī)律。該網(wǎng)格點(diǎn)在t時(shí)刻的加權(quán)度量值為

(7)

因?yàn)槊總€(gè)可移動(dòng)點(diǎn)的運(yùn)動(dòng)是一個(gè)隨機(jī)的過程，因此上述測(cè)量值也是一個(gè)隨機(jī)量，為了對(duì)每個(gè)網(wǎng)格點(diǎn)進(jìn)行最終的度量，需要用這個(gè)隨機(jī)量的均值來代表度量值，因此有

(8)

式中T為總時(shí)間長(zhǎng)度?；谏鲜鏊惴ǖ慕Y(jié)果，將度量值相對(duì)突出的這些點(diǎn)連接起來實(shí)現(xiàn)對(duì)所有IF的提取。

考慮到噪聲和脊線邊緣的影響，需要設(shè)定一個(gè)閾值用于將脊線從低度量值中區(qū)分出來。

(9)

式中SH表示設(shè)定的閾值，由于一般情況下噪聲和脊線邊緣的測(cè)度值遠(yuǎn)小于脊線處的測(cè)度值，因此可以使用整個(gè)分布測(cè)度值的均值或是均值的倍數(shù)作為閾值。經(jīng)過上述處理后的脊線點(diǎn)已經(jīng)清晰顯示出來，考慮瞬時(shí)頻率是沿著時(shí)間軸方向緩慢變化的曲線，從網(wǎng)格左端時(shí)刻開始沿著時(shí)間軸提取與前一刻脊點(diǎn)頻率坐標(biāo)相近的脊線點(diǎn)作為下一時(shí)刻脊點(diǎn)頻率值，并將二者連接成一條脊線，重復(fù)這一過程直到所有的脊線均被找到。

另外一個(gè)值得注意的是，由于一些具有較高能量的隨機(jī)噪聲可能在脊線提取過程中被凸顯出來，形成局部極值點(diǎn)孤立的存在于時(shí)頻矩陣中，但是由于隨機(jī)噪聲脊線的長(zhǎng)度較信號(hào)脊線的長(zhǎng)度要短很多，因此在提取IF時(shí)須設(shè)定一個(gè)最短長(zhǎng)度，小于該長(zhǎng)度的脊線認(rèn)為是隨機(jī)噪聲，從而消除這些局部高能量隨機(jī)噪聲對(duì)脊線提取的影響。

2 仿真信號(hào)處理實(shí)例

2.1 構(gòu)造仿真信號(hào)

圖2給出了信號(hào)發(fā)生與接收的示意圖，在列車速度為亞聲速的情況下，考慮列車軸承聲源為單極子點(diǎn)聲源，并且傳播介質(zhì)為理想流體即不存在粘滯性，沒有能量損耗。根據(jù)莫爾斯聲學(xué)理論[19]，從波動(dòng)方程和運(yùn)動(dòng)關(guān)系出發(fā)，略去影響甚微的高階項(xiàng)后可以推導(dǎo)得出以下公式

(10)

式中P為麥克風(fēng)處采集到的聲壓，q=q0·sin(2πf0t)為單極子點(diǎn)聲源質(zhì)量總流率，q′=?q/?t，R為發(fā)聲時(shí)刻聲源與麥克風(fēng)之間的距離，θ為發(fā)聲時(shí)刻聲源與麥克風(fēng)連線與聲源運(yùn)動(dòng)方向之間的夾角，M=v/c為聲源運(yùn)動(dòng)的馬赫數(shù)。由式(10)可推導(dǎo)出信號(hào)頻率的表達(dá)式[9]

f=f0·

(11)

式中Y為鐵軌與麥克風(fēng)的垂直距離，S為聲源起始點(diǎn)與麥克風(fēng)之間的水平距離。

仿真中采用3個(gè)頻率彼此非常接近的信號(hào)(f1=1 000 Hz,f2=1 100 Hz,f3=1 200 Hz)以保證其不會(huì)被帶通濾波器所提取。預(yù)先給定的參量為：S=8 m,Y=2 m,c=340 m/s和v=20 m/s，仿真信號(hào)的時(shí)域波形如圖4所示。

圖4 仿真信號(hào)

圖5 原信號(hào)與重采樣信號(hào)

圖6 仿真信號(hào)瞬時(shí)頻率與非線性擬合曲線

2.2 仿真結(jié)果

對(duì)仿真信號(hào)以4 000 Hz頻率進(jìn)行采樣，圖5(a)給出了原信號(hào)的頻譜圖，從中可以明顯看到由多普勒效應(yīng)引起的頻移。圖5(b)給出的STFT時(shí)頻分布圖也明顯受到了多普勒效應(yīng)的影響。 對(duì)仿真信號(hào)進(jìn)行CC算法處理，以3倍整體測(cè)度值的均值作為處理閾值，處理得到的離散瞬時(shí)頻率及其非線性擬合曲線如圖6所示。選擇f2=1 100 Hz的畸變信號(hào)，根據(jù)擬合曲線并結(jié)合式(11)可以得到f0,Y,S和v四個(gè)參量，其值在表1中給出。將得到的參量代入式(3)可得到重采樣時(shí)間序列tr，由此對(duì)原信號(hào)進(jìn)行重采樣，經(jīng)過重采樣后信號(hào)的頻譜和STFT譜圖分別如圖5(c)和5(d)所示，可以看到經(jīng)過本文方法處理后的信號(hào)已經(jīng)完全消除了多普勒效應(yīng)的影響。由此可以看出，在處理包含多種頻率成分的多普勒畸變信號(hào)時(shí)，本文方法無需預(yù)測(cè)量其他參數(shù)，即可得出正確的診斷結(jié)果，相比楊和Dybala的方法具有明顯的優(yōu)勢(shì)[7～9]。

表1 非線性曲線擬合參數(shù)與給定參數(shù)的比較

3 實(shí)驗(yàn)驗(yàn)證

本實(shí)驗(yàn)平臺(tái)由作者所在團(tuán)隊(duì)自主設(shè)計(jì)開發(fā)，被測(cè)軸承型號(hào)為目前列車上廣泛使用的NJ(P)3226X1，聲音采集采用B&K公司生產(chǎn)的4944-A型麥克風(fēng)，對(duì)被測(cè)軸承的內(nèi)圈和外圈分別進(jìn)行線切割，寬度均為0.18 mm，如圖7所示。實(shí)驗(yàn)分兩組進(jìn)行，分別采集內(nèi)圈故障軸承和外圈故障軸承的聲音信號(hào)，圖8(a)為故障軸承測(cè)試平臺(tái)，軸承由電機(jī)驅(qū)動(dòng)(轉(zhuǎn)速為1 430 r/min)，在一定負(fù)載作用下(負(fù)載為3 t)產(chǎn)生聲音信號(hào)。使用NI公司開發(fā)的數(shù)據(jù)采集系統(tǒng)(DAS)采集數(shù)據(jù)信號(hào)(采樣頻率為50 kHz)。圖8(b)為實(shí)驗(yàn)場(chǎng)景，測(cè)試平臺(tái)和揚(yáng)聲器被固定在移動(dòng)的汽車上，汽車沿直線以恒定速度行駛，產(chǎn)生的聲音信號(hào)被固定的麥克風(fēng)采集。

圖7 內(nèi)外圈故障

圖8 實(shí)驗(yàn)平臺(tái)和實(shí)驗(yàn)過程

被測(cè)軸承的具體參數(shù)如表2所示，對(duì)于軸承內(nèi)外圈單點(diǎn)故障特征頻率可有下式獲得

(12)

式中fr為軸承旋轉(zhuǎn)頻率，即1 430 r/min，由此可計(jì)算得到外圈故障特征頻率為138.7 Hz，內(nèi)圈故障特征頻率為194.8 Hz。

分析實(shí)驗(yàn)信號(hào)和消除多普勒效應(yīng)的具體步驟如圖1所示。然而實(shí)際實(shí)驗(yàn)環(huán)境遠(yuǎn)比仿真情況復(fù)雜，實(shí)驗(yàn)環(huán)境更易參雜各種頻率的噪聲，因此這里僅選擇頻率成分相對(duì)較高的部分進(jìn)行分析。此外實(shí)驗(yàn)中的采樣頻率非常高，如果直接對(duì)采樣所得數(shù)據(jù)進(jìn)行分析，顯然超出了現(xiàn)有計(jì)算機(jī)的能力范圍。對(duì)圖9中采集到的軸承外圈故障信號(hào)進(jìn)行頻譜分析，從圖10(a)中可看到，外圈故障軸承的頻率成分在1 000～2 000 Hz部分相對(duì)較高，因此這一部分是大家感興趣的部分，所以選擇2 500 Hz以下的頻率成分進(jìn)行分析，對(duì)信號(hào)低通濾波后以5 000 Hz進(jìn)行降采樣。為了更清楚地觀察信號(hào)的時(shí)頻分布，圖11(a)給出了頻率為800～1 600 Hz之間的STFT時(shí)頻分布，從圖中可以看到存在2條瞬時(shí)頻率曲線，經(jīng)過STFT-CC方法處理后的瞬時(shí)頻率估計(jì)曲線如圖11(b)所示，由非線性擬合得到的4個(gè)未知參量分別為f0=1 251 Hz;v=30.26 m/s;Y=2 m以及S=6 m，由此可以根據(jù)1.2節(jié)中介紹的時(shí)域重采樣原理確定重采樣時(shí)間序列，對(duì)原信號(hào)重采樣后即可消除多普勒效應(yīng)，圖11(c)給出了去除多普勒效應(yīng)后的信號(hào)STFT譜圖，對(duì)比圖11(a)，可以看到由多普勒效應(yīng)引起的頻移已得到較大的改善。

表2 NJ(P)3226X1軸承參數(shù)

圖9 列車軸承外圈故障信號(hào)

重采樣后信號(hào)的頻譜圖和包絡(luò)譜圖如圖10(c)和(d)所示，圖中得到的軸承外圈故障特征頻率為fo=138.3 Hz，與理論值138.7 Hz非常接近，而在原信號(hào)包絡(luò)譜圖10(b)中卻由于多普勒效應(yīng)的存在而不能得到正確的診斷結(jié)果。因此本實(shí)驗(yàn)驗(yàn)證了本文方法去除信號(hào)多普勒效應(yīng)的有效性。

圖10 外圈故障軸承信號(hào)頻譜分析

圖11 外圈故障信號(hào)瞬時(shí)頻率估計(jì)

對(duì)于內(nèi)圈故障軸承信號(hào)的處理與外圈故障軸承信號(hào)的處理極為相似，這里就直接給出分析結(jié)果而不再贅述。圖(12)為采集到的內(nèi)圈故障軸承的時(shí)域信號(hào)。從圖13(a)的頻譜圖中可以確定感興趣的頻帶范圍：1 300～2 100 Hz，與之對(duì)應(yīng)的瞬時(shí)頻率估計(jì)如圖14所示。圖13(b)為原信號(hào)的包絡(luò)譜圖，從中僅能得到軸承的旋轉(zhuǎn)頻率fr，而故障信息由于多普勒效應(yīng)引起的頻帶展寬而被湮沒了，因此對(duì)信號(hào)進(jìn)行包絡(luò)分析不能給出有效的故障診斷信息。而在重采樣后信號(hào)的包絡(luò)譜圖13(d)中，由于多普勒效應(yīng)的明顯降低，使得故障信息集中地體現(xiàn)出來，可以明顯地看到內(nèi)圈故障特征頻率fi=194 Hz，其與理論值194.8 Hz很接近。同時(shí)從圖14(c)中也可以看出信號(hào)的多普勒效應(yīng)得到了有效的去除。因此軸承內(nèi)圈故障的實(shí)驗(yàn)再次表明了本文方法的有效性；同時(shí)實(shí)驗(yàn)的結(jié)果也再次顯示了本文方法相比現(xiàn)有的兩種多普勒校正方法在適用范圍和可操作性上的明顯優(yōu)勢(shì)。

圖12 列車軸承內(nèi)圈故障信號(hào)

圖13 內(nèi)圈故障軸承信號(hào)頻譜分析

圖14 內(nèi)圈故障信號(hào)瞬時(shí)頻率估計(jì)

4 結(jié) 論

本文提出了一種基于STFT-CC時(shí)頻估計(jì)和時(shí)域重采樣去除多普勒效應(yīng)的方法，以此提高列車道旁軸承故障診斷的診斷效果。文中首先介紹了時(shí)域重采樣理論和STFT-CC算法。STFT-CC可以通過非線性擬合的方式同時(shí)獲取多個(gè)信號(hào)的時(shí)頻估計(jì)曲線，這很好地解決了帶通濾波器不能同時(shí)獲得多條瞬時(shí)頻率曲線的問題，避免了希爾伯特變換方法中存在端點(diǎn)效應(yīng)的缺點(diǎn)；最后根據(jù)獲得的時(shí)頻曲線推導(dǎo)出曲線的對(duì)應(yīng)參數(shù)，將其應(yīng)用于時(shí)域重采樣原理，可以得到重采樣時(shí)間序列。因此時(shí)域重采樣的過程不需要對(duì)其他參量進(jìn)行預(yù)測(cè)量，大大提高了診斷的可行性和穩(wěn)定性。

文中以3個(gè)頻率彼此相近的信號(hào)對(duì)本文方法進(jìn)行了仿真，并在實(shí)際測(cè)試平臺(tái)上對(duì)外圈故障和內(nèi)圈故障的軸承分別進(jìn)行了實(shí)驗(yàn)，仿真和實(shí)驗(yàn)結(jié)果表明了本文方法在去除多普勒效應(yīng)和道旁軸承故障診斷方面的有效性，展現(xiàn)了本文方法相比現(xiàn)有方法的優(yōu)勢(shì)所在。

參考文獻(xiàn)：

[1] Xavier C, Fabrice B, Jean P D. Early detection of fatigue damage on rolling element bearings using adapted wavelet [J]. Trans. ASME Journal of Vibration and Acoustics, 2007, 129(4):495—506.

[2] Lei Yaguo, He Zhengjia, Zi Yanyang. Application of a novel hybrid intelligent method to compound fault diagnosis of locomotive roller bearings [J]. Trans. ASME Journal of Vibration and Acoustics, 2008, 130 (3):1—6.

[3] Yan Ruqiang, Gao R. Rotary machine health diagnosis based on empirical mode decomposition [J]. Trans. ASME Journal of Vibration and Acoustics, 2008, 130(2):021007.

[4] Choe H C, Wan Y, Chan A K. Neural pattern identification of railroad wheel-bearing faults from audible acoustic signals:comparison of FFT, CWT and DWT features [J]. SPIE Proceedings on Wavelet Applications, 1997, 3 087:480—496.

[5] Sneed W H, Smith R L. On-board real-time bearing defects detection and monitoring[A]. Proceedings of the 1998 ASME/IEEE Joint Railroad Conference[C]. Philadelphia,1998:149—153.

[6] Barke D, Chiu W K. Structural health monitoring in the railway industry:a review [J]. Structural Health Monitoring, 2005, 4:81—93.

[8] Yang Diange, Wang Ziteng, Li Bing, et al. Quantitative measurement of pass-by noise radiated by vehicles running at high speeds [J]. Journal of Sound and Vibration, 2011, 330:1 352—1 364.

[10] Djurovic′ I. Viterbi algorithm for chirp-rate and instantaneous frequency estimation [J]. Signal Processing, 2011, 91:1 308—1 313.

[11] 張旻,程家興,樊甫華,等.利用Hilbert變換提取信號(hào)瞬時(shí)特征參數(shù)的問題研究 [J]. 電訊技術(shù),2003,4:44—48.Zhang Min, Cheng Jiaxing, Fan Puhua, et al. Study on the problems in extracting instantaneous characters of signals based on Hilbert transform [J]. Telecommunication Engineering, 2003,4:44—48.

[12] Wang Yanxue, He Zhengjia, Zi Yanyang. A comparative study on the local mean decomposition and empirical mode decomposition and their applications to rotating machinery health diagnosis [J]. Trans. ASME Journal of Vibration and Acoustics, 2010, 132(2):021010.

[13] 趙義正,楊景曙.基于分?jǐn)?shù)階Fourier變換的瞬時(shí)頻率估計(jì)方法 [J]. 安徽大學(xué)學(xué)報(bào)(自然科學(xué)版),2005, 29(1):44—49.Zhao Yizheng, Yang Jingshu. A method of instantaneous frequency estimation based on fractional Fourier transformation [J]. Journal of Anhui University Natural Science Edition, 2005, 29(1):44—49.

[14] 王磊,郝士琦,戎雁. 瞬時(shí)頻率的Prony方法提取及MATLAB實(shí)現(xiàn) [J]. 計(jì)算機(jī)仿真, 2008,25(2):303—305,309.Wang Lei, Hao Shiqi, Rong Yan. Prony method for instantaneous frequency detection and its implementation based on matlab[J]. Computer Simulation,2008,25(2):303—305.

[15] 趙曉平,趙秀莉,侯榮濤,等. 一種新的旋轉(zhuǎn)機(jī)械升降速階段振動(dòng)信號(hào)的瞬時(shí)頻率估計(jì)算法 [J]. 機(jī)械工程學(xué)報(bào), 2011, 47(7):103—108.Zhao Xiaoping, Zhao Xiuli, Hou Rongtao, et al. A new method for instantaneous frequency estimation of run-up or run-down vibration signal for rotating machinery [J]. Journal of Mechanical Engineering, 2011, 47(7)：103—108.

[16] Carmona R A, Hwang W L, Torresani B. Multiridge detection and time-frequency reconstruction [J]. Signal Processing, 1999, 47(2):480—492.

[17] 張曉冬,王橋,吳樂南. 用于信號(hào)特征提取和重建的脊提取算法 [J].電子與信息學(xué)報(bào), 2003, 25(7):878—883.Zhang Xiaodong, Wang Qiao, Wu Lenan. Characterization extraction and reconstruction of signal by the ridges of continuous wavelet transform [J]. Journal of Electronics and Information Technology, 2003, 25(7):878—883.

[18] 丁立新,康立山,陳毓屏,等. 演化計(jì)算研究進(jìn)展 [J]. 武漢大學(xué)學(xué)報(bào)(自然科學(xué)版), 1998,44(5):561—568.Ding Lixin, Kang Lishan, Chen Yuping, et al. Research development of evolutionary computation [J]. Journal of Wuhan University (Natural Science Edition), 1998, 44(5):561—568.

[19] Morse P M, Ingard K U. 理論聲學(xué)(下冊(cè))[M]. 楊訓(xùn)仁,呂如榆,戴根華,譯. 北京:科學(xué)出版社, 1986:822—850.