基于能量重心法的列車軸承多普勒畸變故障聲信號校正診斷研究

張 翱， 胡 飛， 沈長青， 劉 方， 何清波， 孔凡讓

(中國科學技術大學 精密機械與精密儀器系，合肥 230027)

列車軸承故障是列車故障的主要類型，也是影響列車安全的最大根源之一[1]。美國的一項統計表明，每年大約有50起跟軸承有關的列車出軌事故發生，而且這個數據已經持續保持了近10年的時間[2]。因此，加強軸承的監測和診斷，及時了解和掌握軸承的工作狀態，可以盡量發揮軸承的工作潛力，避免或減少事故的發生，對列車的安全運行具有十分重要的意義。

軸承聲音檢測系統形成于20 世紀80 年代，系統組成如圖1 所示，采用軌邊監測麥克風陣列來采集軸承運行時發出的聲音，通過分析軸承聲音信號來判斷軸承的運行狀況。但由于麥克風放置位置與鐵軌的不可忽略的約一米的距離以及軸承聲源移動相對麥克風存在橫向速度，產生了不同于雷達等領域多普勒效應的畸變聲音信號，因此對其進行分析和校正是進行精確的軸承故障信號特征提取和診斷的前提。

圖1 高速列車軸承聲音檢測系統

20世紀90年代開始，Stojanovic等[3]提出了用鎖相環技術(PLL)進行多普勒聲音信號校正，隨后Johnson等[4]在此基礎上提出了鎖相環技術與DFE算法相結合的校正方法，該方法適用于通信級別的信號校正；楊殿閣等[5]提出非線性時間映射法，基于聲場中的運動學幾何關系，建立聲源與測量信號之間的非線性時間映射方法，采用運動聲源的聲源特征函數，在時域消除多普勒效應，但是該方法需要預測量列車行駛速度、初始位置和麥克風與鐵軌橫向距離等參數。近期Dybaa等[6-7]提出了局部擾動頻率的概念，基于Hilbert變換求解瞬時頻率對畸變信號進行重采樣糾正，但該頻域修正方法具有較多局限性，在特征頻率分布比較密集的情況下難以進行有效的帶通濾波，而且噪聲的存在和濾波器的缺陷造成瞬時頻率曲線的波動，從而帶來校正誤差。

本文提出了一種基于能量重心法的多普勒畸變校正方法，并將其應用在列車軸承道旁故障診斷中。首先由基于能量重心法的瞬時頻率(IF)估計(IFE)來獲取IF向量，在莫爾斯聲學理論的基礎上，對IF向量進行非線性擬合，進而得到重采樣時間間隔序列，再對原信號進行重采樣處理，即可實現多普勒畸變的消除。最后通過對仿真信號及列車軸承內外圈故障聲信號的分析處理，驗證了此種方法的有效性。

1 多普勒效應

當波源和觀察者有相對運動時，觀察者接受到的振動頻率與波源振動頻率不同的現象稱為多普勒效應。

在列車速度為亞聲速的情況下，考慮列車軸承聲源為單極子點聲源，并且傳播介質為理想流體，即不存在粘滯性，沒有能量損耗，其運動模型示意圖如圖2所示。根據莫爾斯聲學理論[8]，從波動方程和運動關系出發推導可以推到得出以下公式：

(1)

式中，P為麥克風處采集到的聲壓值，q為單極子點聲源的強度，t為運行時刻，R為發聲時刻聲源與麥克風之間的距離，c為聲音在介質中速度，θ為聲源與麥克風連線與聲源運動方向之間的夾角，V為聲源速度，M=V/c為馬赫數。式中第二項為小量，可以忽略不計。對于簡諧聲源q=q0sin(ω0t)有：

(2)

圖2 運動模型示意圖

在式(2)中，乘積符號左邊的部分決定信號的幅值，右邊的部分決定信號的相位。對相位進行求導即可得出頻率隨時間的變化。

(3)

式中，x為采集開始時刻時聲源所在位置與麥克風位置水平距離，r為聲源運動軌跡與麥克風的垂直距離，f0為運動聲源信號頻率，f為麥克風采集到的聲音信號頻率。由式(3)可以看出接收到的聲音信號頻率對比原信號頻率呈現非線性畸變，因此本文所提出的方法即是針對此種信號多普勒畸變現象進行消除。

2 基于瞬時頻率的多普勒畸變消除原理

信號的頻率非平穩性的還原通常使用重采樣方法，而建立一組重采樣時間序列就是該方法的核心。對于多普勒畸變信號，瞬時頻率與原信號頻率(假設原信號為單一頻率f0信號)存在如下關系[9]：

(4)

式中,n為周期內采樣點數,fs為原信號采樣頻率,fsi為畸變信號i點處的采樣頻率，fi為畸變信號i點處的瞬時頻率。

采樣時間間隔即是采樣頻率倒數，即dt=1/fs，代入上式可得:

const=fi×dti=f0×dt

i=1,…,N-1

(5)

式中，dti即為重采樣時間間隔，因此重采樣時間間隔序列可以得出dtrsp=[dt1dt2… dtN-1]，從該時間間隔可以進一步計算出重采樣時間點序列trsp=[0t1t2…tN-1]，重采樣時間點序列以畸變信號的起始點為起始點，因此trsp(1)=0，計算公式如下:

(6)

由于重采樣時間點超過畸變信號時間上限時，重采樣將失去意義，因此定義重采樣時間點序列上限為tM，即trsp=[0t1t2…tM]，其中M是不大于N-1的正整數，其值應滿足:

(7)

多普勒畸變信號的瞬時頻率f在任意區間[trsp(i)trsp(i+1)]內是連續變化的，所以采用某時刻的瞬時頻率將會帶來計算誤差，設存在一個較大的整數值K，在時間dti/K內可以認為瞬時頻率是恒定值，則:

(8)

對區間[trsp(i)trsp(i+1)]內的K個時段求和可得:

i=1,…,M-1

(9)

當K值足夠大時，上式可以用積分表示為:

i=1,…,M-1

(10)

由trsp(1)=0可進一步得到trsp(i) 的表達式為:

i=1,…,M-1

(11)

對上式進行求解即可得出重采樣時間點序列trsp=[0t1t2…tM]，最后通過三次樣條函數差值法對多普勒畸變信號x(t)進行重采樣，從而得出校正后信號y(t)。

y=[y(dt)y(2×dt) …y(M×dt)]=[x(trsp(1))x(trsp(2)) …x(trsp(M))]

(12)

3 基于能量重心法的瞬時頻率估計(IFE)

瞬時頻率( IF )是非平穩信號分析中的一個重要概念，Ville于1948年提出的瞬時頻率定義式是目前在學術界最為常用且得到了普遍認可。瞬時頻率估計(Instantaneous Frequency Estimation-IFE)方法主要分為相位法、時頻分布(Time-Frequency Distributions, TFD)法和希爾伯特黃變換(Hilbert-Huang Transform, HHT)[10]。相位法是利用信號的相位信息求取瞬時頻率，主要有相位差分法和相位建模法。基于時頻分布求取信號的瞬時頻率的方法有用時頻分布的一階矩求取瞬時頻率和峰值估計瞬時頻率兩種。HHT則是根據Hilbert變換求出實信號對應的解析信號，再利用解析信號求解瞬時頻率。目前提供的瞬時頻率估計方法大多只適用于單分量信號，對于多分量信號則先將其分解為單分量信號再進行計算。

對工程實際應用中的多分量信號，目前則主要是采用在STFT時頻面上進行峰值搜索來估計瞬時頻率。但是這種方法的估計精度依賴于頻率分辨率鄰近頻率成分的干擾也會對估計精度產生影響，存在時頻集聚性不是很理想的問題。另外，在求取多分量信號的IF時，一般通過時頻濾波將其他分量遮掩，提取某一分量的IF，依次遞推，增加了算法的計算量[11]。

離散頻譜校正是用來解決離散頻譜分析中時域截斷和頻域柵欄效應引起的誤差、精確地提取頻率成分的方法[12]。能量重心法則是根據對稱窗函數離散頻譜的能量重心特性推導出的一種離散頻譜校正方法[13]，對于Hanning窗、矩形窗、Hamming窗、Blackman窗、Blackm an-Harris窗等對稱窗函數而言，離散窗譜的能量重心都在原點或原點附近，利用這個特性可以進行頻率的校正，幅值的校正則可利用巴什瓦定理進行。這種方法在分析平穩或準平穩信號時，具有較高的分析精度。但列車軸承聲信號屬于非平穩信號，不能直接使用離散頻譜校正方法。因此，根據STFT算法原理，把短時間間隔的信號看成是準平穩信號，再通過疊代就能夠利用離散頻譜校正方法來校正整段信號的頻譜。本文通過加Hanning窗的能量重心法在頻域中對信號的頻譜進行瞬時頻率估計(IFE)，并通過離散頻譜校正方法來提高分析精度。

3.1 能量重心法的原理

Hanning窗的功率譜函數為:

(13)

對任意一確定值f，G(f)滿足下式:

n=∞

(14)

該式表明，Hanning窗離散頻譜的能量重心無窮逼近坐標原點。由于Hanning窗旁瓣的功率譜值很小，根據其能量重心的特性，用主瓣內功率譜值較大的幾條譜線可精確地求得主瓣的中心坐標。

對諧波信號其歸一化加Hanning窗的頻譜模函數的平方為:

(15)

相當于式(13)乘以系數A并平移到f=f0處，f0和A分別為分析諧波信號的歸一化頻率和幅值。

Hanning窗譜頻率校正原理如圖3所示，設Gi為功率譜第i條譜線值，Gk為主瓣內譜線最大值，k為幅值最大點對應的譜線號。

圖3 Hanning窗譜頻率校正

根據Hanning窗的能量重心特性有:

(16)

根據式(16)可以求得主瓣的中心:

(17)

設采樣頻率為fs，做譜點數為N，f0為主瓣中心，由式(17)就能得到能量校正法校正頻率的通用公式[14]:

n=∞

(18)

3.2 基于能量重心法進行IFE的流程

基于能量重心法進行IFE的流程如下:

②選定初始搜索頻率fsta頻率搜索范圍q,整周期采樣系數k, 根據fsta和k計算出整周期采樣點數N0；

③根據當前m，計算當前分段數據的起始點位置：p=(m-1)dm+1。在信號中截取M點數據{sig(i)},i∈[p,p+M-1]；

④以分析段中點為中心重新選取N0點作為新的分析段，并對其做加窗的DFT，求得搜索范圍q內最大頻率譜線位置以及其左右鄰近的n條譜線的功率值。其中，n的取值與加窗的類型有關，這里加Hanning窗，n=2；

⑤利用能量重心法校正公式(18)對搜索到的峰值頻率譜線進行頻譜校正，得到校正頻率f，令f為新的fsta，重新計算出整周期采樣所需的采樣點數N；

⑥如果N=N0或者達到最大迭代次數j時，f為分析段中點所對應的頻率估計值，否則fsta=f，轉步驟④；

⑦令m=m+1，返回步驟③，直到m=Num。

最終根據迭代求得的各段對應的IF，根據莫爾斯聲學理論進行非線性插值擬合后得到IF擬合值。

4 基于IFE的多普勒畸變信號校正方法

4.1 多普勒畸變信號校正流程

多普勒畸變信號校正的流程圖如圖4所示，具體校正方法如下：

①對原始信號進行預處理。軸承聲音信號中包含大量噪聲，因此在分析前首先對其進行濾波處理。由于信號中趨勢項的存在，會使時域中的相關分析或頻域中的頻率譜分析產生很大誤差[15]，因此還應進行去除趨勢項處理；

②采用STFT進行時頻分析。通過時頻譜，確定出步驟③所需求的初始搜索頻率fsta，頻率搜索范圍q,整周期采樣系數k；

③基于能量重心法的IFE。由3.2節所訴方法進行IFE，獲得各段數據對應的IF；

④非線性擬合。根據莫爾斯聲學理論進行非線性插值擬合后得到IF擬合值；

⑤多普勒畸變校正。由擬合后的IF計算出重采樣的時間間隔序列，對原始信號進行重采樣，以消除多普勒畸變；

⑥包絡譜分析。通過包絡譜解調出被調制的故障頻率，判斷出故障類型，由此驗證本方法在列車軸承故障診斷中的有效性。

圖4 多普勒畸變信號校正流程圖

4.2 仿真信號模型的建立

根據公式(3)及莫爾斯聲學理論，在此建立一個含有三個頻率的信號進行仿真分析，其中設置這三個頻率相近(f1=1 200 Hz,f2=1 300 Hz,f3=1 400 Hz)，由于經多普勒畸變后頻寬有重疊，不能夠簡單地通過帶通濾波器獲得其中任何一個頻率變化。采樣頻率設定為fs=50 kHz，這也是在后續實驗中對列車軸承信號進行聲音采集的采樣頻率。給定仿真參數x=5 m，r=2 m，c=340 m/s，以及v=20 m/s，仿真的原始信號時域圖如下所示。

圖5 帶多普勒畸變的仿真原始信號

4.3 仿真信號校正結果

在本文中，由多普勒效應引起的時域幅值變化不是重點內容，因此在后續的信號處理中僅僅針對頻域。由圖6(a)可以看出仿真的多普勒畸變原始信號頻率分布在1 000 Hz至1 600 Hz之間,其主要頻寬Δf=342 Hz。由圖6(b)原始信號的STFT圖，我們可以清楚的看到三個設定頻率的多普勒畸變現象。選取中間的頻率進行IFE，在時刻t=0時，峰值大概出現在1 380 Hz處，因此搜索算法中設定1 380 Hz為起始點，設定頻率搜索范圍q為20 Hz，再通過基于能量重心法瞬時頻率估計得到各段對應的IF。

圖6 校正前后仿真信號的頻譜對比

對得到的IF根據公式(12)進行非線性最小二乘法擬合，其IF擬合圖如圖7所示。擬合函數中有四個未知量，即是f0,r,v,x，通過擬合，我們即可獲得這些參數。將擬合所得參數與設定的參數進行比較，如表1所示。

圖7 中間頻率的IF

表1 設定參數與擬合后所得參數對照表

最后，根據第二節描述的多普勒畸變信號校正方法，結合非線性最小二乘法擬合所得到的IF曲線，求解得出重采樣時間點序列trsp，然后通過三次樣條函數差值法對畸變信號的重采樣，從而得出校正后信號。校正后信號的頻譜及STFT圖如圖6(c)和圖6(d)所示，可以明顯的看出，三個頻率均得到了很好的校正，而且從表1可以看出各參數的誤差均符合要求，從而驗證了該多普勒畸變校正方法的有效性。

5 道旁列車軸承聲音信號故障診斷

我國列車使用的軸承是單列向心短圓柱滾子軸承，主要使用的型號是NJ(P)3226X1，因此基于該型號軸承本項目組自行設計了一套實驗平臺(圖9(a))，用于獲取靜態列車軸承聲音信號。實驗中麥克風選用丹麥B&K公司的聲壓場麥克風4944-A，采集卡選用美國NI公司的PXI-4472動態信號采集模塊，采集箱選用美國NI公司的PXI-1033機箱。軸承故障采用線切割方式人為加工而成，軸承內外圈的切縫均為0.18 mm,如圖8所示。

圖8 列車軸承內外圈故障

圖9 多普勒畸變信號獲取

實驗中電機轉速設置為1 430 r/min，軸承加載負荷為3 t,采樣頻率為50 kHz，由列車軸承實驗平臺采集到的信號為不含多普勒畸變的靜態信號。如圖9(b)所示，汽車以108 km/h的速度沿直線勻速行駛，其上固定一音響，以播放上面采集到的靜態信號，麥克風安置于與汽車行駛軸向相距大約1.5 m處。

根據滾子軸承的運動關系可知，軸承外圈和內圈的故障頻率可以由以下公式得出：

(19)

(20)

其中,fr是軸的旋轉頻率，實驗中采用的是1 430 r/min,其他參數見實驗所用的列車軸承規格參數表即表2。從式(19)，式(20)我們可以得出理論的外圈故障頻率應為138.74 Hz，內圈故障頻率應為194.93 Hz。

表2 列車軸承NJ(P)3226X1規格參數

對實驗信號進行多普勒畸變校正的方法與上節所討論的仿真信號校正基本是一致的，不同的是，實驗信號比仿真信號包含更多的頻率分量以及大量噪聲，且是非平穩信號。因此在對其進行校正前，要首先進行去除趨勢項及濾波預處理。圖10為外圈故障信號，在圖12(a)中我們可以清楚的看到,主要的頻率分量集中在1 000 Hz到2 000 Hz之間，因此我們在這里只考慮3 000 Hz以下的頻率分量。圖11(a)為外圈故障信號的時頻分布圖，圖中我們可以看到兩個瞬時頻率分量，由此就可以確定初始搜索頻率fsta為1 350 Hz,選定搜索范圍q為20 Hz，采用上節提到的方法進行多普勒校正。

圖10 列車軸承外圈故障信號

圖11 瞬時頻率估計

在圖12(b)中可以看到校正前信號在故障頻率處有個頻寬為Δf的多普勒畸變，校正后的信號頻譜及包絡譜見圖12(c) 和(d)，圖中可以清楚的看到外圈故障頻率f0為138.9 Hz,多普勒畸變已經基本消除，并且結果與理論值138.74 Hz十分接近。

圖12 軸承外圈故障信號包絡譜分析

下面對內圈故障信號進行分析，其時域圖如圖13所示。根據圖15(a)，信號的主要頻率分布在1 200 Hz到2 200 Hz之間，因此采用6 000 Hz的低通濾波器進行濾波。在圖15(b)中，可以清楚的看到旋轉頻率fr，這是因為在低頻的時候，多普勒效應不是十分明顯，但是內圈故障頻率在理論值194.93 Hz附近沒有明顯的峰值。由圖14(a)可以大致確定初始搜索頻率fsta為1 350 Hz,選定搜索范圍q為20 Hz，通過校正之后，在圖15(d)中我們就可以清楚的看到故障頻率fi，但故障頻率被旋轉頻率調制了，因此有邊頻帶。此時fi為194.5 Hz，與理論值194.93 Hz也十分接近。

圖13 列車軸承內圈故障信號

圖14 瞬時頻率估計

圖15 軸承內圈故障信號包絡譜分析

6 結 論

本文提出了一種針對高速列車軸承聲音信號的多普勒畸變校正方法，并將其應用在列車軸承道旁故障診斷中。首先由基于能量重心法的瞬時頻率(IF)估計(IFE)來獲取IF向量，在莫爾斯聲學理論的基礎上，對IF向量進行非線性擬合，進而得到重采樣時間間隔序列，再對原信號進行重采樣處理，即可實現多普勒畸變的消除。通過對仿真信號進行分析處理后，由表1可以看出該方法的誤差在容許的范圍內。而通過對軸承信號進行多普勒畸變校正后能準確地判斷出軸承的故障類型，證明了該方法在軌邊列車軸承故障診斷中的應用是有效可行的。

與基于STFT峰值搜索及基于WVD峰值搜索的IFE相比，本文采用基于能量重心法的IFE有更高的精度，且分析精度受頻率分辨率影響小，同時在信號的頻域進行分析處理的方法，與以往在時域處理的方法相比，具有無需知道r、v和x這些參數的優點。但由于需要通過手動設置初始搜索頻率、搜索范圍和初始的整周期采樣系數，也就無法實現離線無人診斷操作。為彌補此局限性，筆者將進行進一步的研究改進。

參 考 文 獻

[1]Choe C H, 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, April, 1997, Vol. 3087, pp. 480-496.

[2]Irani F D, et al.先進道旁車輛狀態監視系統的開發和應用[J].國外鐵道車輛，2002, 39(2): 39-45.

Irani F D, et al. Development and deploymen t of advanced wayside condition monitoring systems[J].Foreign Rolling Stock, 2002, 39(2): 39-45.

[3]Stojanovic M, Catipovic J A, Proakis J G. Phase-coherent digital communications for underwater acoustic channels[J]. Oceanic Engineering, 1994, 19(1): 100-111.

[4]Johnson M, Freitag L, Stojanovic M. Improved Doppler tracking and correction for underwater acoustic communications[J]. IEEE International Conference on Acoustics, Speech, and Signal Processing, 1997, 1: 575-578

[5]楊殿閣，羅禹貢，李 兵，等. 基于時域多普勒修正的運動聲全息識別方法[J]. 物理學報，2010, 59(07): 4738-4747.

YANG Dian-ge, LUO Yu-gong, LI Bing , et al.Acoustic holography method for measuring moving sound source with correction for Doppler effect in time-domain[J]. Acta Physica Sinica, 2010: 59(07): 4738-4747.

[8]Morse P, Ingard K. Theoretical Acoustics[M]. Science Press, Beijing, 1986.

[9]Cline J E, Bilodeau J R. Acoustic Wayside Identification of Freight Car Roller Bearing Detects[C]. Proceedings of the 1998 ASME/IEEE Joint Railroad Conference, 1998. 79-83.

[10]賈繼德，孔凡讓，王建平等.基于瞬時頻率估計的內燃機信號階比分析[J].內燃機工程, 2005, 26(3): 15-21.

JIA Ji-de, KONG Fan-rang, WANG Jian-ping. Order analysis of internal combustion engine signal based on instantaneous frequency estimation[J]. Chinese Internal Combustion Engine Engineering, 2005, 26(3): 15-21.

[11]郭 渝，秦樹人.無轉速計旋轉機械升降速振動信號零相位階比跟蹤濾波[J].機械工程學報, 2004, 40(3): 51-54.

GUO Yu, QIN Shu-ren. Tacholess order tracking filtering for run-up or coast down vibration signal of rotating machinery based on zero-phase distortion digital filtering[J]. Chinese Journal of Mechanical Engineering, 2004, 40(3): 51-54.

[12]丁 康，張曉飛.頻譜校正理論的發展[J].振動工程學報, 2000, 13(1): 14-22.

DING Kang, ZHANG Xiao-fei. Advances in spectrum correction theory[J].Journal of Vibration Engineering, 2000, 13(1): 14-22.

[13]丁 康，江利旗.離散頻譜的能量重心校正法[J].振動工程學報, 2001, 14(3): 354-358.

DING Kang, JIANG Li-qi. Energy centrobaric correction method for discrete spectrum[J]. Journal of Vibration Engineering, 2001, 14(3): 354-358.

[14]段虎明，秦樹人，李 寧.離散頻譜的校正方法綜述[J].振動與沖擊, 2007, 26(11): 138-145.

DUAN Hu-ming, QIN Shu-ren, LI Ning.Review of correction methods for discrete spectrum[J]. Journal of Vibration and Shock, 2007, 26(11): 138-145.

[15]趙寶新，張保成，趙鵬飛,等.EMD在非平穩隨機信號消除趨勢項中的研究與應用[J].機械制造與自動化，2009, 38(5): 85-87.

ZHAO Bao-xin, ZHANG Bao-cheng, ZHAO Peng-fei, et al. Research on and application of EMD in eliminating trend Item of non-stationary random signals[J]. Machine Building and Automation, 2009, 38(5): 85-87.