Crazy Climber算法與重采樣技術在消除多普勒效應及列車軸承診斷中的應用

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

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

引 言

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

然而當列車高速運行經過ADBD時，其診斷的有效性將有所降低，其中一個重要的原因就是列車與麥克風陣列之間的高速相對運動引起的多普勒效應，多普勒效應會在時域和頻域上對采集到的信號產生較大的畸變和頻移，這無疑將嚴重影響診斷效果[7]。因此為了提高列車高速運行中麥克風采集到信號的可用性，進而提高列車軸承狀態監測和故障診斷的有效性，多普勒效應的消除勢在必行。

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

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

文中首先介紹了聲學運動學模型、時域重采樣方法和基于STFT-CC的瞬時頻率估計方法；然后針對提出的方法進行了信號仿真；最后將此方法應用于實際故障軸承的診斷中。仿真和實驗結果均表明本文方法的可行性與有效性。

1 診斷方法和原理

1.1 診斷方法

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

圖1 消除信號多普勒效應流程圖

1.2 時域重采樣原理

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

列車經過麥克風時的情況如圖2所示，這里假定列車以恒定速度沿筆直鐵軌運行，僅考慮單一聲源和單一麥克風的情況。

圖2 列車經過麥克風示意圖

由圖2可知，在聲源發聲時刻與麥克風接收時刻之間存在一段延時，隨著聲源和麥克風之間距離的縮短，其延時時間將逐漸變短，這意味著信號在時間軸上被壓縮了，因此信號頻率提高了，反之當聲源遠離麥克風時，信號在時間軸上被拉伸，因此信號頻率降低，這就是運動聲源測量中多普勒效應的產生原因。

圖2中點A即初始位置，此處時間為0。假設在時間為ts時，聲源發出的信號幅值為x(ts)，此信號經過傳播在時間tr時到達麥克風，因此有下式成立

tr=ts+Δt

(1)

式中 Δt即聲音信號在空氣中的傳播時間，可由下式得到

(2)

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

(3)

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

重采樣的過程通過對原信號進行三次樣條插值完成，如圖3所示。圖中tm為麥克風實際采樣時間序列，tp為重采樣時間序列，xp為原信號幅值序列，xr為校正信號的幅值序列。

圖3 重采樣過程示意圖

時域重采樣的具體步驟如下：

(1)以聲源在A點出發時刻為0時刻，校正麥克風采集信號時間原點，即tm=tr+R0/c，其中R0為出發時刻聲源與麥克風之間的距離，c為聲速。

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

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

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

1.3 基于STFT-CC的瞬時頻率估計

定義一個頻率調制信號為

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

(4)

則信號的瞬時頻率定義為信號相位的一階導數ω(t)=φ′(t)。信號的瞬時頻率是一個非常重要的參量，其估計的準確程度關系到對其他信號參量的求取[10]。

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

獲取信號瞬時頻率的一般方法是求取信號時頻分布的譜峰值，然而該方法不能同時提取多條瞬時頻率曲線，且在高噪聲和復雜工況環境下，時頻分布譜峰檢測法并不能給出理想的瞬時頻率估計[15]。

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

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

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

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

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

(5)

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

當系統溫度的值低于某個設定的閾值時，這個迭代移動過程結束。

當整個迭代過程結束后，需要對網格上的每個點進行度量，度量的方式存在兩種。第一種稱之為均值度量，在某一時刻，考慮每個移動點相當于質量為1/N的質點，則可以度量該網格點在t時刻的度量值為

(6)

第二種度量方式稱之為加權度量，該方式考慮到了將原始的時頻分布平面上的數據用來對上式的結果進行加權，期望不僅能提取到脊所在的正確的位置，而且脊上參數的變化能更加貼近原來的時頻平面上參數的變化規律。該網格點在t時刻的加權度量值為

(7)

因為每個可移動點的運動是一個隨機的過程，因此上述測量值也是一個隨機量，為了對每個網格點進行最終的度量，需要用這個隨機量的均值來代表度量值，因此有

(8)

式中T為總時間長度。基于上述算法的結果，將度量值相對突出的這些點連接起來實現對所有IF的提取。

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

(9)

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

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

2 仿真信號處理實例

2.1 構造仿真信號

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

(10)

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

f=f0·

(11)

式中Y為鐵軌與麥克風的垂直距離，S為聲源起始點與麥克風之間的水平距離。

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

圖4 仿真信號

圖5 原信號與重采樣信號

圖6 仿真信號瞬時頻率與非線性擬合曲線

2.2 仿真結果

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

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

3 實驗驗證

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

圖7 內外圈故障

圖8 實驗平臺和實驗過程

被測軸承的具體參數如表2所示，對于軸承內外圈單點故障特征頻率可有下式獲得

(12)

式中fr為軸承旋轉頻率，即1 430 r/min，由此可計算得到外圈故障特征頻率為138.7 Hz，內圈故障特征頻率為194.8 Hz。

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

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

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

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

圖10 外圈故障軸承信號頻譜分析

圖11 外圈故障信號瞬時頻率估計

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

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

圖13 內圈故障軸承信號頻譜分析

圖14 內圈故障信號瞬時頻率估計

4 結 論

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

文中以3個頻率彼此相近的信號對本文方法進行了仿真，并在實際測試平臺上對外圈故障和內圈故障的軸承分別進行了實驗，仿真和實驗結果表明了本文方法在去除多普勒效應和道旁軸承故障診斷方面的有效性，展現了本文方法相比現有方法的優勢所在。

參考文獻：

[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變換提取信號瞬時特征參數的問題研究 [J]. 電訊技術,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] 趙義正,楊景曙.基于分數階Fourier變換的瞬時頻率估計方法 [J]. 安徽大學學報(自然科學版),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] 王磊,郝士琦,戎雁. 瞬時頻率的Prony方法提取及MATLAB實現 [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] 趙曉平,趙秀莉,侯榮濤,等. 一種新的旋轉機械升降速階段振動信號的瞬時頻率估計算法 [J]. 機械工程學報, 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] 張曉冬,王橋,吳樂南. 用于信號特征提取和重建的脊提取算法 [J].電子與信息學報, 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]. 武漢大學學報(自然科學版), 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. 理論聲學(下冊)[M]. 楊訓仁,呂如榆,戴根華,譯. 北京:科學出版社, 1986:822—850.