基于總變差去噪和快速譜相關(guān)的滾動軸承故障診斷

唐貴基， 田 甜， 龐 彬

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

作為旋轉(zhuǎn)機(jī)械的重要部件，滾動軸承在工業(yè)生產(chǎn)中占有重要地位。滾動軸承故障導(dǎo)致機(jī)器性能惡化，因此有必要準(zhǔn)確地檢測軸承的故障。滾動軸承的故障信號通常都是非平穩(wěn)信號[1]，采集到的滾動軸承故障信號中包含許多噪聲信號，常規(guī)信號分析方法難以有效提取出滾動軸承的故障特征。

循環(huán)平穩(wěn)隨機(jī)過程是一類特殊的非平穩(wěn)隨機(jī)過程。它的統(tǒng)計特性雖然是非平穩(wěn)的，但卻隨時間的變化呈現(xiàn)出周期性平穩(wěn)變化。循環(huán)平穩(wěn)信號的一個主要特性是：即使在頻譜重疊的情況下，它們也可以很容易地與其他干擾信號分離。機(jī)械信號循環(huán)譜分析的一個核心工具是譜相關(guān)(Spectral Correlation，SC)，它以雙譜圖的形式顯示信號中調(diào)制信號和載波信號的整個結(jié)構(gòu)[2-8]。雖然循環(huán)譜分析在狀態(tài)監(jiān)測中的能力已經(jīng)在一些研究工作中得到證實,然而在某些情況下，SC計算成本高，影響了它的應(yīng)用。Antoni等[9]提出了SC的替代方法——循環(huán)調(diào)制譜(Cyclic Modulation Spectrum，CMS)。CMS彌補了SC計算成本高的缺陷，并被證明在許多情況下，可以作為一個有效的診斷工具。然而CMS卻受不確定性原理的約束，只能檢測調(diào)制頻率低于頻率分辨率的周期模式。

Antoni等[10]提出了快速譜相關(guān)方法(Fast Spectral Correlation，F(xiàn)ast-SC)進(jìn)一步豐富了循環(huán)平穩(wěn)分析理論，此方法同時彌補了譜相關(guān)計算成本高和CMS受不確定性約束的缺陷。然而滾動軸承故障信號通常包含噪聲成分，筆者在研究中發(fā)現(xiàn)噪聲的存在會嚴(yán)重影響快速譜相關(guān)方法的分析效果。總變差去噪方法對噪聲十分敏感，可以很好地剔除信號中的噪聲成分，并保留良好的信號邊緣，近些年被成功引入到故障診斷領(lǐng)域[11-12]。本文利用其改善快速譜相關(guān)對于含強背景噪聲的滾動軸承故障振動信號的分析效果。首先利用總變差去噪對含噪的循環(huán)平穩(wěn)信號消除噪聲成分，然后通過對降噪信號進(jìn)行快速譜相關(guān)分析提取滾動軸承的故障特征頻率，最后應(yīng)用仿真與實驗分析驗證了此方法的有效性。

1 總變差去噪原理

TVD(Total Variation Denoising)是基于一個優(yōu)化問題的定義，最初由Rudin[13]提出，用于去除圖像中的噪聲。此后，由于TVD在保留給定信號的尖銳邊緣方面的優(yōu)勢，TVD在一維信號處理得到了廣泛的采用。

給定一個一維信號x(n)(0≤n≤N-1)，則信號x(n)的總變差可定義為

(1)

矩陣D的定義如下

(2)

則DDT為三對角陣

(3)

假設(shè)信號x(n)受到加性高斯白噪聲w(n)的污染，則含噪信號y(n)為

y(n)=x(n)+w(n)

(4)

最小化優(yōu)化目標(biāo)函數(shù)為

(5)

(6)

優(yōu)化最小化(Majorization-Minimization，MM)是解決難以直接求解優(yōu)化問題的有效方法。MM方法的思想是選擇一個接近F(x)且易于求解的Gk(x)。該方法產(chǎn)生一個序列xk，每個xk是通過最小化Gk-1(x)得到。參考文獻(xiàn)[14]，總結(jié)MM算法的步驟如下

步驟1 設(shè)k=0,初始化x0；

步驟2 選擇Gk(x)，使：

a)Gk(x)≥F(x),?x

b)Gk(xk)=F(xk)

c)Gk(x)為凸函數(shù)

步驟4k=k+1,重復(fù)步驟(2)。

(7)

式中：Λk∶=diag(|Dxk|)

將式(7)左右同乘λ，得

(8)

(9)

結(jié)合式(5)中給出的優(yōu)化目標(biāo)函數(shù)可以表示為

(10)

結(jié)合式(6)可得

(11)

2 快速譜相關(guān)原理

2.1 譜相關(guān)理論

SC可定義為瞬時自相關(guān)函數(shù)的雙離散傅里葉變換，也就是關(guān)于時間t的傅里葉級數(shù)以及關(guān)于時間滯后τ的傅里葉變換。

以時域信號x(tn)為例，其譜相關(guān)定義如下

(12)

式中:Rx(tn,τm)為x(tn)的瞬時自相關(guān)函數(shù);Rx(tn,τm)=E{x(tn)x(tn-τm)*}，tn=n/Fs,τm=m/Fs；Fs為采樣頻率；E為集總平均；τm為時延因子；*為共軛復(fù)數(shù)；α為循環(huán)頻率；f為頻率。

2.2 快速譜相關(guān)

SC的一個估計方法是平均循環(huán)周期圖法(Averaged Cyclic Periodogram，ACP)。平均循環(huán)周期圖的定義如下

(13)

ACP雖然可以同時吸收不同頻率下的頻譜分量，估計結(jié)果和SC接近，但是有著計算時間長的缺點。

為彌補ACP計算時間長，計算成本大的缺點，Antoni等提出了Fast-SC。它的原理如下

(1) 假設(shè)滿足f=fk=kΔf,令α=pΔf+δ，則|Xw(i,fk-α)|≈|Xw(i,fk-p)|；

(2) 利用一階泰勒展開式可得

綜上，

Xw(i,fk-α)=|Xw(i,fk-α)|ejφi(fk-α)?

(14)

將式(14)代入式(13)，得出掃描譜相關(guān)的定義式

(15)

式中：XSTFT(i，fk)為信號x(tn)的短時傅里葉變換；DFT(Discrete Fourier Transform)為離散傅里葉變換。

Sx(α,fk;p)在循環(huán)頻率區(qū)間[(p-1)Δf,(p+1)Δf]掃描SC，表明它對p的幾個值的聚合，從而在整個循環(huán)頻率范圍內(nèi)重建SC。p的最大值為

(16)

由此定義快速譜相關(guān)

(17)

其中核函數(shù)

(18)

由式(18)可知，當(dāng)Nw→∞，K→∞時，快速譜相關(guān)的結(jié)果近似于譜相關(guān)的結(jié)果。快速譜相關(guān)是譜相關(guān)的漸近收斂方差，具有類似ACP的統(tǒng)計性能，但計算量大大減少。因為快速譜相關(guān)易于并行計算，允許進(jìn)一步的潛在加速，因此計算速度快。離散傅里葉變換用零填充來估計Sx(α,fk;p)，從而提高循環(huán)頻率α的數(shù)值分辨率，以便降低柵欄效應(yīng)，所以比ACP的估計效果還要好。

一種基于快速譜相關(guān)的增強包絡(luò)譜(Enhanced Envelope Spectrum，EES)，定義為

(19)

快速譜相關(guān)算法最后得到一個快速譜相關(guān)圖和一個增強包絡(luò)譜圖，通過這兩個譜圖識別故障特征信息。

3 仿真信號分析

為驗證本文方法的有效性，對式(20)所示的仿真信號s(t)進(jìn)行模擬。

(20)

表1 外圈故障模型中的參數(shù)

為體現(xiàn)本文所選去噪方法的優(yōu)越性，將仿真信號總變差去噪效果和小波閾值去噪效果進(jìn)行了對比分析。小波閾值降噪采用MATLAB工具箱中的“wdpen”

圖1 仿真信號的時域波形

函數(shù)實現(xiàn)。小波閾值去噪是一種分頻段濾波方法，它的原理是首先利用小波把原信號分解成不同頻段的尺度系數(shù)，對小波系數(shù)進(jìn)行閾值處理再重構(gòu)，得到去噪后信號[16-18]。圖2為s(t)的波形圖；圖3為總變差去噪的時域圖；圖4(a)、圖4(b)、圖4(c)分別為db4小波、db5小波和db6小波閾值降噪的時域圖。通過圖3和圖4的對比可知，總變差去噪明顯優(yōu)于db4小波和db5小波的去噪效果。由表2的評價結(jié)果可知，總變差去噪的信噪比(Signal-Noise Ratio,SNR)最大、均方根誤差(Root Mean Squared Error,RMSE)最小，進(jìn)一步證明了總變差降噪的效果優(yōu)于db4小波、db5小波和db6小波閾值降噪的效果[19]。

圖2 s(t)的波形圖

圖3 總變差去噪后的波形圖

本文分別用三種方案對仿真信號進(jìn)行分析:方案一,對仿真信號直接進(jìn)行快速譜相關(guān)分析；方案二,首先用總變差去噪方法對信號s(t)進(jìn)行降噪，然后對去噪后的信號進(jìn)行快速譜相關(guān)分析。方案三,通過表2可知db6小波閾值比db4、db5降噪效果更好。因此首先用db6小波閾值去噪方法對信號s(t)進(jìn)行降噪，然后對去噪后的信號進(jìn)行快速譜相關(guān)分析。

(a) 小波閾值db4

(b) 小波閾值db5

(c) 小波閾值db6

Tab.2 Evaluation results of different evaluation methods for different denoising results

總變差去噪小波閾值db4小波閾值db5小波閾值db6RMSE1.56151.59751.60471.5681SNR8.72158.52338.48428.6846

圖5(a)為直接對s(t)進(jìn)行快速譜相關(guān)得出的譜相關(guān)圖，由此圖觀察到在外圈故障特征頻率處出現(xiàn)了能量集中；圖5(b)為其相應(yīng)的增強包絡(luò)譜圖，可識別出外圈故障特征頻率，但無法識別故障特征頻率的倍頻成分。圖6(a)為本文方法所得的譜相關(guān)圖，可看出在故障特征頻率及其倍頻處都出現(xiàn)了能量集中；圖6(b)為其相應(yīng)的增強包絡(luò)譜圖，可有效地識別故障特征頻率及其倍頻。總變差去噪和快速譜相關(guān)的結(jié)合使故障沖擊頻率更加明顯，效果很好。圖7(a)、圖7(b)分別為小波閾值去噪和快速譜相關(guān)結(jié)合的效果圖，圖7(a)中能量混雜，無法識別出故障特征頻率及其倍頻，圖7(b)只能識別出轉(zhuǎn)頻和二倍故障頻率。由此可見，總變差去噪和快速譜相關(guān)相結(jié)合的方法優(yōu)于直接使用快速譜相關(guān)方法，且優(yōu)于小波閾值去噪和快速譜相關(guān)相結(jié)合的方法。

4 實驗驗證

利用圖8所示的實驗平臺對本文所述方法進(jìn)行進(jìn)一步驗證。實驗平臺主要由電機(jī)、轉(zhuǎn)子、加載器及軸承組成。在滾動軸承內(nèi)圈模擬滾動軸承的局部損傷，圖9為內(nèi)圈故障圖。本文所分析的故障信號由電渦流傳感器測量主軸所得。電渦流采集振動信號的方式為非接觸測量方式，更加安全便捷。采樣頻率為12 800 Hz。表3為滾動軸承結(jié)構(gòu)參數(shù)，根據(jù)式(21)計算可得內(nèi)圈理論故障特征頻率fi為172 Hz。

(a) 直接快速譜相關(guān)

(b) 增強包絡(luò)譜

圖5 直接快速譜相關(guān)的結(jié)果

Fig.5 Results of using fast spectral correlation directly

(a) TVD降噪后快速譜相關(guān)

(b) 增強包絡(luò)譜

圖6 本文方法的分析結(jié)果

Fig.6 Results of analysis by proposed method

(a) 小波閾值降噪后快速譜相關(guān)

(b) 增強包絡(luò)譜

圖7 小波去噪與快速譜相關(guān)結(jié)合

Fig.7 Results of analysis of combination of the wavelet threshold and fast spectral correlation

圖8 實驗平臺

軸承型號節(jié)圓直徑/mm滾動體直徑/mm接觸角/(°)滾動體個數(shù)轉(zhuǎn)軸轉(zhuǎn)頻/HzN205397.501224

(21)

圖9 內(nèi)圈故障

式中：n為滾動體個數(shù)；fr為轉(zhuǎn)軸轉(zhuǎn)頻；d、D分別為滾動體直徑及節(jié)圓直徑；α為接觸角。

這里也給出了同小波閾值降噪的對比分析結(jié)果。圖10為實測信號的波形圖；圖11為總變差去噪的時域圖；圖12(a)、圖12(b)、圖12(c)分別為db4小波、db5小波和db6小波閾值降噪的時域圖。通過圖11和圖12的對比可知，總變差去噪明顯優(yōu)于db4小波、db5小波和db6小波的去噪效果。小波閾值降噪效果取決于小波基函數(shù)的選擇，當(dāng)基函數(shù)選取不合適的時候，降噪效果便會欠佳。由表4的評價結(jié)果可知，總變差去噪的SNR最大、RMSE最小，進(jìn)一步證明了總變差降噪的效果優(yōu)于db4小波、db5小波和db6小波閾值降噪的效果。

圖10 實測信號的波形

圖11 總變差去噪后的波形圖

(a) 小波閾值db4

(b) 小波閾值db5

(c) 小波閾值db6

圖12 小波閾值去噪后時域圖

Fig.12 Time-domain figure after the wavelet threshold denoising

表4 兩種評價方法對不同去噪結(jié)果的評價結(jié)果

Tab.4 Evaluation results of different evaluation methods for different denoising results

總變差去噪小波閾值db4小波閾值db5小波閾值db6RMSE39.951641.562440.107941.5681SNR21.791821.448521.757921.6846

首先利用快速譜相關(guān)直接對實測信號進(jìn)行處理，結(jié)果如圖13所示。圖13(a)為譜相關(guān)圖，可看出在故障特征頻率及其倍頻處有能量集中；圖13(b)為相應(yīng)的增強包絡(luò)譜，可識別出內(nèi)圈故障頻率和2倍頻，3～5倍頻則不夠明顯。下面利用本文方法對該組實測信號進(jìn)行處理。首先利用總變差去噪方法對信號進(jìn)行去噪，對去噪后的信號進(jìn)行快速譜相關(guān)，得到圖14所示的結(jié)果，圖14(a)為降噪信號的快速譜相關(guān)譜，此圖在故障特征頻率及其倍頻處都出現(xiàn)了明顯的能量集中；圖14(b)為其相應(yīng)的增強包絡(luò)譜，可看到在故障特征頻率及其倍頻處具有明顯峰值，并且辨識出了轉(zhuǎn)頻及內(nèi)圈故障特征頻率通轉(zhuǎn)頻的調(diào)制特征，可斷定軸承內(nèi)圈出現(xiàn)損傷。為體現(xiàn)本文方法的優(yōu)越性，通過小波閾值去噪和快速譜相關(guān)的結(jié)合效果與之作比較。由表4可知，小波db5信噪比大于db4、db6且均方根誤差小于db4、db6。因此先將實測信號進(jìn)行db5小波閾值去噪，然后將去噪后的信號進(jìn)行快速譜相關(guān)。圖15(a)為該方法的譜相關(guān)圖，可看到故障頻率及其倍頻處有能量集中；圖15(b)為相應(yīng)的增強包絡(luò)譜，只可看到轉(zhuǎn)頻處和故障頻率及其二倍頻處有明顯的峰值，相對于圖13(b)沒有明顯的改善，效果不如圖14(b)。以上結(jié)果說明，相比于直接快速譜相關(guān)和小波閾值去噪與快速譜相關(guān)相結(jié)合的方法，本文方法更好地提取了故障特征。

(a) 直接快速譜相關(guān)

(b) 增強包絡(luò)譜

圖13 直接快速譜相關(guān)分析結(jié)果

Fig.13 Results of using fast spectral correlation directly

(a) TVD降噪后快速譜相關(guān)

(b) 增強包絡(luò)譜

圖14 本文方法的分析結(jié)果

Fig.14 Results of analysis by proposed method

(a) 小波閾值降噪后快速譜相關(guān)

(b) 增強包絡(luò)譜

圖15 小波去噪與快速譜相關(guān)結(jié)合

Fig.15 Results of analysis of combination of the wavelet threshold and fast spectral correlation

5 結(jié) 論

快速譜相關(guān)是無偏無方差的譜相關(guān)估計方法，不僅克服了平均循環(huán)周期圖計算成本大的問題，還克服了循環(huán)調(diào)制譜受不確定原理約束的缺陷。在機(jī)械故障診斷領(lǐng)域具有巨大的應(yīng)用潛能。但該方法容易受到強背景噪聲的影響，因此本文將總變差去噪方法作為預(yù)處理方法，提出了總變差去噪結(jié)合快速譜相關(guān)的滾動軸承微弱故障診斷方法。利用總變差去噪提高振動信號的信噪比，再通過快速譜相關(guān)方法識別故障特征信息。仿真信號與實驗信號結(jié)果表明本文方法極大改善了快速譜相關(guān)方法的分析效果，促進(jìn)其在滾動軸承微弱故障診斷中的應(yīng)用。