徐亞軍,于德介,孫云嵩,趙 丹
(湖南大學汽車車身先進設計制造國家重點實驗室,湖南 長沙 410082)
滾動軸承是旋轉機械中應用最為廣泛的部件之一,各種旋轉機械輸出軸的支撐大多采用滾動軸承形式。滾動軸承的運行狀態對整臺機器的安全運行影響最大[1],其運行狀態的正常與否將直接影響到整個機組的性能,因此滾動軸承的狀態監測和故障診斷技術研究具有重要的工程應用價值。
滾動軸承發生故障時,會產生周期性的脈沖沖擊力,從而產生振動信號的調制現象,頻譜圖上表現為固有頻率兩側出現間隔均勻的調制邊頻帶[2]。傳統的軸承故障診斷方法通過對軸承振動信號進行解調診斷故障。常用的軟件解調方法有廣義檢波濾波解調、希爾伯特變換解調。在廣義檢波濾波解調分析中,由于取絕對值、檢波過程或平方過程都會使載波頻率有可能出現高次諧波而產生混頻效應,因此需要選擇合適的采樣頻率以避免這種混頻效應[3]。Hilbert解調由于Hilbert算子的加窗效應,使得解調結果出現非瞬時響應特性,即在解調出的調制信號兩端及有突變的中間部位產生調制,表現為幅值按指數衰減波動,從而使解調誤差增大[4]。
數學形態分析是基于積分幾何和隨機集的非線性信號分析方法。該方法通過形態變換將一個復雜信號分解為具有物理意義的各個部分,使其與背景剝離,同時保留信號主要的形狀特征[5]。數學形態學在考察信號時使用結構元素探針,通過結構元素探針在信號中不斷移動來提取有用信息進行特征分析和描述[6,7]。多尺度形態學的結構元素由信號產生,具有一定的適應性,能更有效地提取信號的沖擊特征[8~10]。對故障振動信號進行多尺度形態學操作,提取故障信號的沖擊特征,然后進行頻譜分析,能較好地達到解調的目的。與傳統的解調分析比較,由于算法只涉及加減運算,不需要對信號進行絕對值、Hilbert算子等運算,一方面可以減少由于算子運算產生的混頻效應、加窗效應等;另一方面降低了算法的復雜度。且不需要對運算結果進行低通濾波,無需預先選擇截止頻率,因此多尺度形態學解調是一種很好的解調分析方法。
根據軸承故障特征頻率公式[11],軸承故障特征頻率與軸承的相關尺寸和轉速相關。當轉速變化時,軸承的故障特征頻率會隨轉速一起變化,從而此時的故障特征頻率不再是某一固定值,而是隨時間變化的一條故障頻率曲線。傅里葉變換以及其他的解調方法需要信號滿足平穩化的要求,因此傳統的軸承故障診斷方法不適用于轉速變化情況下的軸承故障診斷。
Emmanuel J Candès等近年來提出的線調頻小波路徑追蹤算法能有效提取信號的瞬時頻率[12]。本文將線調頻小波路徑追蹤算法與多尺度形態學解調方法相結合,提出了基于線調頻小波路徑追蹤的階比多尺度形態學解調方法,并用于轉速大范圍波動情況下的滾動軸承故障診斷。當滾動軸承發生故障時,其載波頻率一般為外環的各階固有頻率,而調制頻率為產生疲勞剝落元件(內環、外環或滾動體)的通過頻率及其倍頻[13]。本文方法先用線調頻小波路徑追蹤算法估計轉速波動滾動軸承的故障特征頻率;再依據獲得的通過頻率對等時間間隔采樣的滾動軸承振動信號進行插值和重采樣,得到等角度采樣的角域平穩信號;最后用信號局部峰值方法確定多尺度形態學分析的結構元素,根據各結構元素對重采樣信號進行形態學操作,進而獲得操作結果平均值的階次譜以診斷滾動軸承故障。仿真信號和實驗信號分析結果表明,基于線調頻小波路徑追蹤的階比多尺度形態學解調的方法可以有效地應用于滾動軸承的故障診斷。

線調頻小波路徑追蹤算法采用的多尺度線調頻基元函數如下


式(1)定義的多尺度線性調頻基函數在動態分析時間段內的瞬時頻率為aμ+2bμt,在時頻面的表示如圖1所示。通過信號在多尺度線性調頻基函數上的逐段投影分析,計算獲得每個時間分析段I內的最大投影系數和對應的線調頻基元函數,該基元函數即為在時間分析段I中與分析信號最為相似的頻率成分,適合在小的動態分析時間段內逐段擬合頻率呈曲線變化的頻率成分。線調頻基函數的多尺度特性使得它具有了動態匹配分析信號的特性,而基函數中包含的調頻率信息則使得其適合分析頻率呈曲線變化的非平穩信號。

圖1 基元函數的時頻表示Fig.1 Time frequency representation of elementary function
當信號與多尺度線性調頻基函數越相似時,其投影系數也越大,基元函數的能量也越大,因此要求找到一種動態分析時間段連接方法,在該連接方法下,使連接的所有基元函數信號在整個分析時間內的總能量最大,即

∏覆蓋整個分析時間段,不能重疊,其對應的最大投影系數和基元函數分別為

∏的連接方法應保證在投影中使連接的基函數在整個分析時間段內的總能量最大,線調頻小波路徑追蹤算法提出的連接算法如下:
(1)初始化。以i為時間支持區序號,d(i)為第i個時間支持區之前分解信號的總能量,pre(i)為連接到第i個時間支持區的前置時間支持區序號,e(i)為第i個時間支持區最大投影系數對應的分解信號的能量,初始化時,置d(i),pre(i)=0;
(2)對于動態分析時間段集合{Ii,i∈Z}中的每一個元素Ii,查找出與其相鄰的所有下一個動態分析時間段集合{Ij},即{Ij}中所有元素的起始時間與Ii相鄰。如果

有

Π的連接方法可以保證在整個分析時間段內基元函數組合形成的信號與分析信號最為相似,而基元函數在動態分析時間支持區∏={I1,I2,…}內的瞬時頻率為aμi+2bμiti,ti∈Ii,對應的線性直線連接形成的頻率曲線則是對分解信號的瞬時頻率估計。
數學形態學有兩種基本的形態函數:腐蝕和膨脹,若待處理信號x(n)是采樣得到的一維多值信號,定義域為Df=0,1,2,…,N-1,g(n)為一維結構元素序列,定義域為Dg=0,1,2,…,M-1,其中N和M都是整數,且有N≥M。信號膨脹、腐蝕、開和閉運算分別定義為

式中m∈0,1,2,…,M-1。
數學形態學的開運算可用于濾除信號上方的峰值噪聲,去除信號邊緣的毛刺;閉運算可用于平滑或抑制信號下方的波谷噪聲,填補信號的漏洞和縫隙[16]。同時利用形態學開、閉運算可以構造AVG和DIF濾波器,構造方式分別為

式中 AVG濾波器能用于信號正負沖擊特征的平滑,相當于平滑濾波器,而DIF能用于提取信號的沖擊特征[17]。
令x和B分別代表一維信號和形態學的結構元素,若給定了一個二值形態學變換T,則基于T的多尺度形態學運算可定義為一族形態學變換{Tλ|λ>0,λ∈N},其中Tλ(x)=λT(x/λ),相似地,多尺度形態學膨脹、腐蝕操作定義為

式中λB=B⊕B⊕…⊕B(λ-1次),參照傳統的多尺度形態學在圖像處理中的應用,定義長度尺度λl和高度尺度λh,即λ=(λl,λh)。對長度尺度λl和高度尺度λh的搜索,本文采用局部極值自適應搜索方法,得到長度尺度λl的最大值和最小值分別為

因此長度尺度λl={λlmin,λlmin+1,…,λlmax},其中in為極值間隔,高度尺度λh定義為

式中j=0,1,2,…,λlmax-λlmin,其中pnmin,pnmax分別為信號的最大和最小值。β為尺度的幅值系數,(0<β<1),本文取β=1/3。利用λl和λh構建多尺度形態學的每一個尺度的結構元素λB為

式中λlB是將B進行λl-1次膨脹運算。
本文選取的形態學算子為差值算子DIF,結構元素B為三角型結構元素,B=[0,1,0]。則

多尺度形態學DIFλ操作能有效地提取信號的沖擊特征,若對滾動軸承故障信號進行多尺度形態學操作,提取故障信號的沖擊特征,然后進行頻譜分析,則能較好地達到解調的目的。但是頻譜分析是以分析穩態信號為前提的,而變轉速工況下拾取的軸承振動信號為非平穩信號。因此本文結合線調頻小波路徑追蹤方法和多尺度形態學分析,提出基于線調頻小波路徑追蹤的多尺度形態學解調方法。
信號等角度重采樣的具體步驟為:
(1)采用線調頻小波路徑追蹤方法分析滾動軸承振動信號,獲得變轉速滾動軸承振動信號的故障特征頻率曲線;
(2)用三階多項式對特征頻率曲線進行擬合,則有

(3)確定階次跟蹤的最大分析階Dmax;
(4)計算等角度重采樣的角度間隔Δθ,根據采樣定理,采樣率應不小于最大分析階次的兩倍,所以

(5)計算重采樣后數據的長度N

式中T為時域采樣的總時間,f(t)為頻率擬合函數;
(6)計算等角度重采樣的鍵相時標Tn

式中T0為時域采樣開始時間;
(7)等角度重采樣,根據所求出的鍵相時標Tn,利用Lagrange線性插值公式對振動信號進行插值,求出振動信號在角域里對應于鍵相時標Tn的幅值。對于給定的插值節點Tn及對應的幅值x(Tn),則Lagrange線性插值公式為

經以上步驟得到重采樣信號x(Tn)后,用信號局部極值自適應搜索方法得到重采樣信號x(Tn)的長度尺度和高度尺度,再利用式(17)即可構建多尺度形態學DIFλ操作的每一個尺度的結構元素λB。
根據式(18)分別計算出分析信號x(Tn)的各個尺度形態學分析結果,其反映了分析信號中含有的特定尺度的沖擊成分,在每一尺度的分析結果中也會含有隨機噪聲。對于特定的信號x(Tn),其沖擊成分應該出現在多個尺度中,因此將各個尺度的形態學分析得到的信號取平均值作為最終的多尺度形態學操作結果,這樣既突出信號中的沖擊成分,又有效地抑制了噪聲。最后對平均后的多尺度形態學操作結果做頻譜分析,便完成了對滾動軸承故障信號的階比多尺度形態學解調。
基于線調頻小波路徑追蹤的階比多尺度形態學解調過程可用圖2所示的原理框圖表示。

圖2 基于線調頻小波路徑追蹤的階比多尺度形態學解調原理框圖Fig.2 The schematic diagram of the order multi-scale morphology demodulation based on chirplet path pursuit algorithm
為驗證本文方法的有效性,對一軸承故障仿真信號進行分析。當軸承出現故障時,振動信號中會出現以軸承元件固有頻率為載波頻率的多頻率調制信號,通常表現為不斷重復的脈沖衰減信號,當軸承轉速變化時,脈沖出現的間隔即調制頻率會隨轉速變化。
取脈沖信號為

信號載波頻率為1 020Hz,衰減系數為-420。為模擬軸承振動信號被調制信號調制的現象,假定載波信號幅值被頻率f1調制,取調制頻率f1為

圖3為單個脈沖信號時域波形圖,通過構建repeat函數來實現每隔周期1/f1出現單個脈沖,即可組成周期脈沖信號。圖4為加噪仿真信號時域波形示意圖,信號采樣頻率為8 192Hz,采樣時長4s,信噪比10dB。對加噪仿真信號進行希爾伯特變換,得到包絡信號,因脈沖調制頻率最大為54Hz,所以將包絡信號采樣頻率由原來的8 192降為1 024Hz,對于脈沖調制頻率的提取仍能滿足采樣定理。圖5為用FFT變換得到的降采樣包絡信號頻譜圖,由于脈沖調制頻率隨時間變化,從該圖很難直接對調制頻率進行識別。圖6為仿真信號的調制頻率和擬合頻率時頻圖,圖中實線是仿真信號的調制頻率曲線,虛線為對包絡信號用線調頻小波路徑追蹤方法提取的調制頻率擬合曲線。從圖中不難看出擬合頻率曲線很好地匹配了信號的調制頻率曲線。將信號以線調頻小波路徑追蹤方法提取的調制頻率進行等角度重采樣,對重采樣信號進行多尺度形態學解調,再對解調信號做階次譜就得到圖7所示的仿真信號的多尺度形態學解調階次譜。圖中階次1,2,3,4恰好對應著模擬故障特征頻率及其倍頻,從而很好地驗證了本文方法的有效性。

圖3 脈沖信號時域波形圖Fig.3 The time domain waveform diagram of pulse signal

圖4 仿真信號時域波形圖Fig.4 The waveform of simulated signal

圖5 包絡信號頻譜圖Fig.5 Spectrogram of the envelope signal

圖6 仿真信號的調制頻率和擬合頻率時頻圖Fig.6 The time-frequency diagram of the simulation signal’s modulation frequency and fitting frequency

圖7 仿真信號的多尺度形態學解調階次譜Fig.7 The order spectrum of the simulation signal obtained by multi-scale morphology demodulation
圖8中實線是仿真信號的調制頻率曲線,虛線為仿真信號基于Wigner-Ville峰值跟蹤算法得到的擬合頻率曲線,與基于線調頻小波路徑追蹤方法得到的調制頻率曲線(圖6)相比,其精度明顯要低。
圖9為仿真信號基于Wigner-Ville的包絡階次譜,圖中出現了(2.705,77.08),(3.599,62.39)的干擾階次。可見在信噪比較低的情況下,基于線調頻小波路徑追蹤的階比多尺度形態學解調方法能得到更準確的階次譜。

圖8 仿真信號的調制頻率與 Wigner-Ville峰值跟蹤算法的擬合頻率Fig.8 The modulation frequency and the fitting frequency obtained by Wigner-Ville peak tracking algorithm of the simulation signal

圖9 基于Wigner-Ville峰值跟蹤的包絡階次譜Fig.9 The envelope order spectrum based on Wigner-Ville peak tracking algorithm
試驗軸承為6307E型滾動軸承,利用激光分別在外圈和內圈上切割寬0.15mm,深0.13mm的槽,模擬外圈和內圈故障,由于試驗條件的限制,未能在滾動體上設置故障。振動加速度傳感器固定在軸承座正上方,拾取軸承外圈和內圈故障振動信號時,轉速分別在809.09~1 280.60r/min之間和340.72~714.91r/min之間變化,采樣頻率8 192 Hz,采樣時長為4s。滾動軸承試驗臺如圖10所示。

圖10 滾動軸承試驗裝置Fig.10 The experiment device of rolling bearing

圖11 軸承外圈故障振動信號時域波形圖Fig.11 The time domain vibration signal of bearing outer race fault
圖11為軸承外圈故障振動信號的時域波形圖,從圖中可以看出存在明顯的沖擊脈沖,且脈沖間隔隨轉速變化。圖12中曲線分別為用轉速計實測的故障軸承所在軸轉頻曲線(曲線1)和根據軸承故障特征頻率計算公式計算獲得的軸承外圈故障特征頻率曲線(曲線2)[11],可以看出轉頻在13.48~21.34 Hz間變化,故障特征頻率隨轉速波動在41.26~65.31Hz之間變化。圖13為傳統的包絡信號頻譜圖,由于故障特征頻率隨轉速變化,無法準確判斷其峰值點屬于轉頻及其倍頻還是故障特征頻率及其倍頻。將包絡信號的采樣頻率由8 192Hz降為1 024 Hz,對于10倍故障特征頻率的提取仍能滿足采樣定理。當滾動軸承外圈發生故障時,其載波頻率一般為外環的各階固有頻率,而調制頻率為外環的通過頻率及其倍頻[13]。圖14中實線為圖12中計算得到的故障特征頻率曲線,虛線是降采樣的包絡信號經過線調頻小波路徑追蹤方法分解得到的擬合頻率曲線。從圖中可以看出其能很好地匹配根據轉頻計算得到的故障特征頻率曲線。以圖14中虛線所示頻率曲線對外圈故障信號進行等角度重采樣,對重采樣信號進行多尺度形態學解調,再對解調信號進行階次譜分析,得到外圈故障信號的階比多尺度形態學解調階次譜如圖15所示,圖中的峰值譜線對應外圈故障特征頻率及其倍頻的階次,從圖中能很清晰地辨別出故障特征頻率及其倍頻,表明軸承信號被外圈故障特征頻率及其倍頻成分調制了,這正是一般滾動軸承出現外圈故障時的振動信號特征,與實際情況相符。

圖12 軸轉動頻率與軸承外圈故障特征頻率曲線圖Fig.12 The Shaft rotation frequency and bearing outer race fault characteristic frequency

圖13 外圈故障包絡信號頻譜圖Fig.13 The envelope spectrum of bearing outer race fault’s vibration signal

圖14 外圈故障特征頻率曲線Fig.14 The characteristic frequency of the bearing outer race fault

圖15 外圈故障信號多尺度形態學解調階次譜Fig.15 The order spectrum of the bearing outer race fault signal obtained by multi-scale morphology demodulation

圖16 外圈故障信號的包絡階次譜Fig.16 The envelope order spectrum of the bearing outer race fault signal
同樣以圖14中虛線所示頻率曲線對外圈故障信號進行等角度重采樣,再對解調信號進行階次譜分析得到圖16所示的外圈故障信號的包絡階次譜。比較圖15與16,可以看出圖15中的階次更加清晰、明顯,既突出了沖擊成分,又有效地抑制了噪聲。從而可以驗證階比多尺度形態學對于軸承故障診斷的優越性與有效性。

圖17 軸承內圈故障振動信號時域波形圖Fig.17 The time domain vibration signal of bearing inner race fault

圖18 內圈故障信號多尺度形態學解調階次譜Fig.18 The order spectrum of the bearing inner race fault signal obtained by multi-scale morphology demodulation

圖19 內圈故障信號基于Wigner-Ville的包絡階次譜Fig.19 The envelope order spectrum of the bearing inner race fault signal based on Wigner-Ville
圖17是含有內圈故障的滾動軸承振動加速度信號。當滾動軸承內圈出現故障時,其載波頻率一般為外環的各階固有頻率,調制頻率為內環的通過頻率及其倍頻以及轉頻及其倍頻[13]。圖18為內圈故障振動信號基于線調頻小波路徑追蹤的多尺度形態學解調階次譜。階次譜中出現了1X1(0.992 5),2X1(1.975)的階次,這與故障特征頻率的1,2倍頻極為吻合。同時也出現了0.190 7,0.410 7等階次,這是振動信號被轉頻及其倍頻調制的結果。實際上,根據軸承故障特征頻率公式和軸承型號(6307E)[11],計算得到內圈的故障特征頻率是轉頻的4.937倍(則轉頻為故障特征頻率的1/4.937≈0.202),階次譜中出現的 X(0.190 7),2X(0.410 7),3X(0.603 3),4X(0.810 6),6X(1.210 8),7X(1.414 3),8X(1.602 1),9X(1.821 4)的階次為轉頻及其倍頻。
用基于Wigner-Ville峰值的階比跟蹤方法對內圈故障信號進行包絡階次譜分析,結果如圖19所示,其階次基本上不能區分出來,已淹沒在噪聲中。可見在內圈故障振動信號信噪比較低的情況下,基于線調頻小波路徑追蹤的階比多尺度形態學解調效果要明顯優于基于Wigner-Ville的階比跟蹤效果。
(1)線調頻小波路徑追蹤算法能在不額外安裝硬件的情況下得到軸承的故障特征頻率曲線,從而能將轉速波動的非平穩振動信號轉化為角域的平穩信號,降低了滾動軸承故障診斷的復雜度,節約了成本。用線調頻小波路徑追蹤算法獲得軸承故障特征頻率曲線的過程不會受到二次時頻交叉干擾項的影響,不需要通過帶通濾波器對其他頻率分量進行“掩蓋”,使得其估計精度高于峰值跟蹤法。
(2)多尺度形態學的尺度系數由信號產生,具有一定的自適應性,能夠提取振動信號中含有的特定尺度的沖擊成分,對各尺度形態學分析結果取平均能有效抑制噪聲,因此多尺度形態學階次分析比包絡階次分析更能有效地突出沖擊成分,抑制噪聲。
(3)仿真與應用實例表明,將線調頻小波路徑追蹤算法和多尺度形態學解調方法相結合,能很好對變轉速下的滾動軸承外圈與內圈故障進行識別。
[1] 鐘秉林,黃仁.機械故障診斷學[M].北京:機械工業出版社,2007.
Zhong Bing-lin,Huang Ren.Mechanical Fault Diagnosis[M].Beijing:China Machine Press,2007.
[2] 丁康,孔正國.振動調幅調頻信號的調制邊頻帶分析及其解調方法[J].振動與沖擊,2005,24(6):9—12.
Ding Kang,Kong Zheng-guo.Modulation principle of amplitude and frequency modulated signal and its demodulation procedure[J].Journal of Vibration and Shock,2005,24(6):9—12.
[3] 張帆,丁康.廣義檢波解調分析的三種方法及其局限性研究[J].振動工程學報,2002,15(2):243—248.
Zhang Fan,Ding Kang.Research on the three algorithms and limitations of generalized detection-filtering demodulation analysis[J].Journal of Vibration Engineering,2002,15(2):243—248.
[4] 劉紅星,陳濤,屈梁生,等.能量算子解調方法及其在機械信號解調中的應用[J].機械工程學報,1998,34(5):85—90.
Liu Hong-xing,Chen Tao,Qu Liang-sheng,et al.Energy operator demodulation method and its application in mechanical signal demodulation[J].Chinese Journal of Mechanical Engineering,1998,34(5):85—90.
[5] 龔煒,石青云,程民德.數字空間中的數學形態學——理論及應用[M].北京:科學出版社,1997.Gong Wei,Shi Qing-yun,Chen Min-de,Theory of Mathematical Morphology in Digital Space and Its Application[M].Beijing:Science Press,1997.
[6] NIKOLAOU N G,ANTONIADIS I A.Application of morphological operators as envelope extractors for impulsive-type periodic signals[J].Mechanical Systems and Signal Processing,2003,17(6):1 147—1 162.
[7] SUSANTA M,BHABATOSH C.A multiscale morphological approach to local contrast enhancement[J].Signal Processing,2000,80:685—696.
[8] Mukhopadhyay S,Chanda B.An edge preserving noise smoothing technique using multiscale morphology[J].Signal Processing,2002,82(4):527—544.
[9] Zhang Lijun,Xu Jinwu,Yang Jianhong,et a1.Multiscale morphology analysis and its application to fault diagnosis[J].Mechanical Systems and Signal Processing,2008,22:597—610.
[10]林京,屈梁生.基于連續小波變換的信號檢測技術與故障診斷[J].機械工程學報,2000,36(12):95—100.
Lin Jing,Qu Liang-sheng.Feature detection and fault diagnosis based on continous wavelet transform[J].Chinese Journal of Mechanical Engineering,2000,36(12):95—100.
[11]MALLAT S,ZHANG Z.Matching pursuit with timefrequency dictionaries[J].Signal Processing,1993,41(12):3 397—3 415.
[12]Candès E J,Charlton P R,Helgason H.Detecting highly oscillatory signals by chirplet path pursuit[J].Applied and Computational Harmonic Analysis,2008,24(1):14—40.
[13]PAN Minchun,CHIU Chunching.Investigation on improved Gabor order tracking technique and its applications[J].Journal of Sound and Vibration,2006,295(3-5):810—826.
[14]Candes E J,Charltion P R,Helgason Hannes.Detecting highly oscillatory signals by chirplet path pursuit[J].Applied and Computational Harmonic Analysis,2008,24(1):14—40.
[15]Qian S,Chen D.Joint Time-Frequency Analysis[M].Englewood Cliffs,Prentice Hall,1996.
[16]章立軍,徐金梧,陽建宏.自適應多尺度形態學分析及其在滾動軸承故障診斷中的應用[J].北京科技大學學報,2008,30(4):441—445.
Zhang Li-jun,Xu Jin-wu,Yang Jian-hong.Adaptive multiscale morphology analysis and its application in fault diagnosis of bearing[J].Journal of Beijing University of Science and Technology,2008,30(4):441—445.
[17]郝如江,盧文秀,褚福磊.滾動軸承故障信號的數學形態學提取方法[J].中國電機工程學報,2008,28(6):65—70.
Hao Ru-Jiang,Lu Wen-xiu,Chu Fu-lei.Mathematical morphology extracting method on roller bearing fault signals[J].Proceedings of the CSEE,2008,28(6):65—70.