基于不完全Cholesky分解相關(guān)熵雙譜的軸承故障診斷

李 輝, 郝如江

(1.天津職業(yè)技術(shù)師范大學(xué) 機械工程學(xué)院，天津 300222; 2.石家莊鐵道大學(xué) 機械工程學(xué)院，石家莊 050043)

滾動軸承是機械傳動中最常用的零件，滾動軸承的健康狀態(tài)是保證機械傳動正常運行的重要因素，因此，檢測滾動軸承的健康狀態(tài)至關(guān)重要，基于振動信號的滾動軸承故障檢測與診斷是最常用的方法之一[1-2]。當滾動軸承內(nèi)圈、外圈或滾動體產(chǎn)生表面點蝕、疲勞脫落、局部裂紋等故障時，由于故障部位對與其接觸的軸承其它部件產(chǎn)生沖擊作用，從而產(chǎn)生周期性的瞬態(tài)脈沖沖擊現(xiàn)象，這種瞬態(tài)沖擊會激勵機電設(shè)備的固有振動，并與軸承的回轉(zhuǎn)頻率產(chǎn)生調(diào)制，使得振動信號產(chǎn)生復(fù)雜的幅值、相位調(diào)制現(xiàn)象，不僅具有非線性、非高斯信號特征，而且往往包含大量的背景噪聲，將軸承微弱故障特征淹沒，這些因素加大了軸承故障診斷的難度。軸承故障振動信號的非線性、非高斯特征使傳統(tǒng)基于信號二階統(tǒng)計量的方法，如相關(guān)分析、功率譜分析等，性能衰退，甚至失效[3-4]。因此，針對強背景噪聲下滾動軸承故障特征難以有效提取的問題，許多學(xué)者進行了研究。王宏超等[5]對二階循環(huán)統(tǒng)計量的譜相關(guān)進行了改進，提出一種新的時頻分析方法。鄢小安等[6]提出了一種基于改進奇異譜分解的形態(tài)學(xué)解調(diào)方法，并用于軸承故障診斷。胥永剛等[7]提出了改進奇異譜分解方法，并成功應(yīng)用于軸承故障診斷。胡超凡等[8]提出了基于張量分解的多維降噪技術(shù)，用于滾動軸承復(fù)合故障多維降噪以及特征提取。在眾多的機械設(shè)備故障診斷方法中，基于信號高階統(tǒng)計量的信號處理方法，由于在非線性、非高斯信號處理方面有著獨特優(yōu)勢，得到了更廣泛的研究和應(yīng)用，如雙譜[9-11]、切片雙譜[12]、倒雙譜[13]、循環(huán)雙譜[14-16]、高階譜[17]分析等，已廣泛用于機械設(shè)備的故障診斷，取得了良好效果。

軸承故障振動信號呈周期性瞬態(tài)沖擊和復(fù)雜調(diào)制的特點，具有明顯的非線性、非平穩(wěn)、非高斯特點。雙譜是分析高斯信號的有效工具，理論上高斯噪聲的雙譜為零, 但非高斯噪聲的雙譜并不為零，傳統(tǒng)的雙譜不能有效消除非高斯噪聲的干擾，因此會造成雙譜分析性能衰退，甚至失效。相關(guān)熵是處理高斯、非高斯噪聲的有效方法，已在雷達和通信信號檢測、信號濾波、波達方向估計和時延估計等方面得到應(yīng)用，取得了良好效果[18-20]。盡管相關(guān)熵方法在通信領(lǐng)域的應(yīng)用已經(jīng)展開，但在機電設(shè)備故障診斷領(lǐng)域的應(yīng)用才剛剛起步[21-22]。本文針對傳統(tǒng)雙譜難以有效處理非高斯噪聲干擾以及相關(guān)熵計算量大的問題，提出了一種基于相關(guān)熵和雙譜的軸承故障診斷方法，該方法首先利用核函數(shù)和不完全Cholesky分解(incomplete Cholesky decomposition, ICD)算法計算核矩陣的低秩分解—下三角矩陣G；其次，利用Gini指數(shù)(Gini index，GI)選取下三角矩陣G的主分量，由下三角矩陣G的主分量計算信號的核矩陣；根據(jù)核矩陣計算信號的相關(guān)熵；最后，計算相關(guān)熵的雙譜，根據(jù)相關(guān)熵的雙譜特征識別軸承故障，并利用仿真信號和電機軸承故障試驗信號驗證了該方法的有效性和可靠性。

1 相關(guān)熵和ICD方法簡介

1.1 相關(guān)熵

對于實信號x(t)，時變自相關(guān)熵可定義為

(1)

式中:Ε(·)為期望算子;τ為時延;κσ(·)通常采用高斯核函數(shù)。κσ(·)表達式為

(2)

(3)

式中，m為時延的采樣點數(shù)，m=0,1,2,…,N-1。

1.2 不完全Cholesky矩陣分解

(4)

式(4)稱為矩陣κ的不完全Cholesky分解。當d<

1.3 基于Gini指數(shù)的矩陣降維

GI最初用于判斷收入分配公平程度的指標，是一個比例數(shù)值，在0和1之間，GI越小，表示信號的幅值越均衡，反之，表示信號的幅值變化越大[24]。因此，可利用GI有效提取軸承故障振動信號中周期性的脈沖沖擊。

向量z的GI定義為

(5)

z(1)≤z(2)≤…≤z(N)

(6)

(7)

(8)

2 基于相關(guān)熵的雙譜估計

設(shè){x(n)}是零均值的平穩(wěn)隨機過程，其三階累積量Rxx(τ1,τ2)可表示為

Rxx(τ1,τ2)=E[x(n)x(n+τ1)x(n+τ2)]

(9)

式中，τ1、τ2為時延的采樣點數(shù)。

信號{x(n)}的雙譜可表示為三階累積量Rxx(τ1,τ2)的傅里葉變換，即

(10)

式中，ω為圓頻率，且|ω1|≤π,|ω2|≤π, |ω1+ω2|≤π。

(11)

(12)

從式(12)可以看出基于相關(guān)熵的雙譜是圓頻率ω1和ω2的函數(shù)。

3 基于相關(guān)熵雙譜的故障診斷步驟

在實際工程應(yīng)用中，根據(jù)以下步驟進行基于相關(guān)熵雙譜的軸承故障診斷：

(5) 根據(jù)雙譜的頻譜結(jié)構(gòu)識別軸承故障。

基于相關(guān)熵雙譜的軸承故障診斷流程圖如圖1所示。

圖1 相關(guān)熵雙譜故障診斷流程圖Fig.1 Flowchart of correntropy bispectrum based fault diagnosis

4 相關(guān)熵雙譜軸承故障診斷仿真

滾動軸承內(nèi)圈、外圈或滾動體產(chǎn)生表面點蝕、疲勞脫落、局部裂紋等是滾動軸承的典型故障。滾動軸承點蝕故障振動響應(yīng)信號模型可以表示為[25]

sin[2πfbi(t-ti)+θbi]

(13)

x(t)=x1(t)+n1(t)+n2(t)

(14)

式中:Ai為瞬態(tài)脈沖幅值;Ci為阻尼衰減因子;ti為沖擊持續(xù)的時間;θbi為初始相位;fbi電動機系統(tǒng)的共振頻率；n1(t)為零均值高斯噪聲;n2(t)為脈沖噪聲。

函數(shù)Θ(t-ti)用來指定沖擊發(fā)生的時間，可用下式定義

(15)

利用信號仿真x(t)驗證相關(guān)熵的降噪性能和本文提出方法的有效性。設(shè)軸承外圈故障特征頻率為fouter=120 Hz，電動機系統(tǒng)的固有振動頻率fb=3 000 Hz，瞬態(tài)沖擊振幅Ai=4.5，采樣頻率fs=20 000 Hz，信號采樣點數(shù)n=4 000。

圖2(a)和圖2(b)分別為仿真信號x1(t)的時域波形和FFT，可以清晰看到信號x1(t)在頻域內(nèi)圍繞軸承系統(tǒng)共振頻率fb形成邊頻帶簇，邊頻帶的間隔為軸承外圈故障特征頻率120 Hz。

(a) 時域圖

(b) 頻域圖圖2 仿真信號x1(t)及其FFTFig.2 Simulative signal x1(t) and its FFT

當軸承外圈發(fā)生故障時，在外圈故障雙譜平面內(nèi)，會產(chǎn)生以系統(tǒng)共振頻率為中心、軸承外圈故障特征頻率為間隔的譜峰簇。圖3為仿真信號x1(t)相關(guān)熵(σ=3)的雙譜，從圖3的雙譜平面內(nèi)可以清晰看到：圍繞系統(tǒng)共振頻率fb，在(±fb,0)、(0,±fb)、(fb,-fb)以及(-fb,fb)等位置附近存在顯著的譜峰簇，相鄰譜峰之間的間隔為軸承外圈故障特征頻率fouter=120 Hz，這種頻譜特征表明了軸承外圈的故障特征，與理論分析一致。

圖3 仿真信號x1(t)的雙譜Fig.3 Bi-spectrum of simulative signal x1(t)

在仿真信號x1(t)中加入零均值高斯噪聲n1(t)，形成信噪比為SNR=-10 dB的染噪信號，之后再隨機加入幾個幅值不同的脈沖信號，以模擬非高斯脈沖噪聲，圖4(a)和圖4(b)分別為仿真信號x2(t)的時域波形和FFT，由于仿真信號x2(t)完全被強高斯和非高斯噪聲淹沒，因此從圖4(a)已完全看不出信號幅值的變化規(guī)律。在圖4(b)中，除在系統(tǒng)共振頻率3 000 Hz位置出現(xiàn)比較大的峰值外，已無法識別在共振頻率兩側(cè)的邊頻帶簇。

(a) 時域圖

(b) 頻域圖圖4 仿真信號x2(t)及其FFTFig.4 Simulative signal x2(t) and its FFT

圖5 矩陣各列的Gini指數(shù)Fig.5 Gini index of each column for matrix

圖6為信號x2(t)的相關(guān)熵 (σ=3)。從圖6可以看出：仿真信號x2(t)相關(guān)熵的時域波形呈現(xiàn)出比較明顯的周期性幅值調(diào)制特征，在圖4(a)中，不但幅值較大的非高斯脈沖噪聲已被完全消除，而且強高斯噪聲也得到了很好地抑制，但相關(guān)熵仍然保留了信號中的周期瞬態(tài)沖擊成分。

圖6 仿真信號x2(t)的相關(guān)熵Fig.6 Correntropy of simulative signal x2(t)

圖7和圖8為仿真信號x2(t)相關(guān)熵的雙譜。從圖7和圖8可以看出，在強高斯噪聲和非高斯噪聲的影響下，基于相關(guān)熵的雙譜仍然能夠準確識別軸承外圈的故障特征，在(±fb,0)、(0,±fb)、(fb,-fb)以及(-fb,fb)等位置附近存在顯著的譜峰簇，這些獨立的譜峰仍然清晰可辨，相鄰譜峰之間的間隔為軸承外圈故障特征頻率fouter=120 Hz，由此可見，在強高斯和非高斯噪聲干擾下，基于相關(guān)熵的雙譜仍然能準確提取軸承外圈的故障特征。

圖7 仿真信號x2(t)相關(guān)熵的雙譜Fig.7 Correntropy based bi-spectrum of signal x2(t)

圖8 仿真信號x2(t)相關(guān)熵的雙譜Fig.8 Correntropy based bi-spectrum of signal x2(t)

為彰顯相關(guān)熵的降噪能力，將基于相關(guān)熵的雙譜與傳統(tǒng)雙譜和文獻[10]中基于小波閾值降噪雙譜進行了對比。圖9和圖10給出了利用直接法計算的仿真信號x2(t)的雙譜，圖11和圖12為采用小波閾值降噪后信號的雙譜，采用5層sym8小波分解。從圖9、圖10和圖11、圖12可以看出，仿真信號x2(t)的傳統(tǒng)雙譜和小波閾值降噪雙譜已完全被強高斯和非高斯噪聲淹沒，完全不能分辨軸承故障特征頻率。對比圖8和圖10、圖12，圖8中基于相關(guān)熵的雙譜噪聲方差很小，而傳統(tǒng)雙譜(圖10)和基于小波閾值降噪雙譜(圖12)的噪聲方差仍然很大，軸承故障頻率已完全被強噪聲淹沒，不能有效識別。由此可見，盡管雙譜對高斯噪聲具有一定的抑制能力，但當信號中含有很強的高斯噪聲和非高斯噪聲時，雙譜也難以有效提取軸承外圈的故障特征。而小波閾值降噪方法，在處理強高斯噪聲和非高斯噪聲干擾時，也難以取得理想的效果。

圖9 仿真信號x2(t)的雙譜Fig.9 Bi-spectrum of simulative signal x2(t)

圖10 仿真信號x2(t)的雙譜Fig.10 Bi-spectrum of simulative signal x2(t)

圖11 仿真信號x2(t)小波閾值降噪雙譜Fig.11 Wavelet de-noising based bi-spectrum of x2(t)

圖12 仿真信號x2(t)小波閾值降噪雙譜Fig.12 Wavelet de-noising based bi-spectrum of x2(t)

通過上述仿真信號可知：由于相關(guān)熵能有效抑制高斯噪聲和脈沖噪聲，因此基于相關(guān)熵的雙譜分析具有從強高斯噪聲和非高斯噪聲背景中提取軸承故障特征的能力，基于相關(guān)熵的雙譜分析為高斯、非高斯噪聲的處理提供了一種嶄新的魯棒性解決方法。

5 電機軸承故障診斷

采用美國凱斯西儲大學(xué)(Case Western Reserve University，CWRU)軸承數(shù)據(jù)中心的網(wǎng)站上公布的試驗數(shù)據(jù)[26]，CWRU軸承數(shù)據(jù)集采集系統(tǒng)如圖13所示。試驗軸承型號為：深溝球軸承6205-2RS JEM SKF，采樣頻率fs=12 000 Hz，電機負載為空載，電機轉(zhuǎn)速1 797 r/min(fr=29.95 Hz)，軸承內(nèi)圈故障振動信號數(shù)據(jù)記錄號為105DE，軸承外圈故障振動信號數(shù)據(jù)記錄號為130DE。計算得到滾動軸承內(nèi)圈、外圈故障特征頻率如下：

finner=162.185 Hz，fouter=107.365 Hz

圖13 CWRU軸承數(shù)據(jù)采集系統(tǒng)Fig.13 Bearing data collection system of CWRU

5.1 軸承內(nèi)圈故障診斷

圖14為滾動軸承內(nèi)圈故障振動信號及其FFT。在圖14(a)中存在明顯的幅值調(diào)制現(xiàn)象，但根據(jù)時域波形還不能識別軸承故障。在圖14(b)中，頻譜幅值最大值處的頻率fb=3 585 Hz，反映了電機系統(tǒng)的共振頻率，但在fb的周圍沒有出現(xiàn)清晰的以軸承內(nèi)圈故障特征頻率為間隔的邊頻帶簇，因此根據(jù)頻域圖不能有效識別軸承的故障類型。

(a) 時域圖

(b) 頻域圖圖14 軸承內(nèi)圈故障振動信號的時域圖和頻譜圖Fig.14 Bearing inner race fault vibration signal and its FFT

圖15 矩陣各列的Gini指數(shù)Fig.15 Gini index of each column for matrix

圖16 軸承內(nèi)圈故障振動信號的相關(guān)熵Fig.16 Correntropy of bearing inner race fault vibration signal

圖17和圖18為軸承內(nèi)圈故障振動信號相關(guān)熵的雙譜，從圖17和圖18可以看出，雙譜圖中的譜峰呈離散點狀分布，幾乎沒有受到噪聲的影響。在低頻范圍，沿f1=0和f2=0軸的兩側(cè)，約在2 720 Hz處存在明顯的譜峰簇，譜峰的頻率等于軸承內(nèi)圈故障特征頻率finner，譜峰之間的間隔為軸承內(nèi)圈故障特征頻率。在高頻范圍，在f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz、f1=±3 585 Hz等10條直線的交點附近，存在明顯的譜峰簇且譜峰之間的間隔為軸承內(nèi)圈故障特征頻率，這種頻譜結(jié)構(gòu)表明了軸承的內(nèi)圈故障特征，同時也說明電動機系統(tǒng)在1 310 Hz、2 720 Hz以及3 585 Hz頻率處具有較大的振動能量，與圖14(b)FFT頻譜圖的結(jié)果一致。

圖17 軸承內(nèi)圈故障振動信號相關(guān)熵的雙譜Fig.17 Correntropy based bi-spectrum of bearing inner race fault vibration signal

圖18 軸承內(nèi)圈故障信號相關(guān)熵的雙譜Fig.18 Correntropy based bi-spectrum of bearing inner race fault vibration signal

圖19和圖20給出了利用直接法計算的軸承內(nèi)圈故障振動信號的雙譜，圖21和圖22為采用小波閾值降噪后信號的雙譜。在圖19、圖20中，因噪聲的影響，可模糊看出f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz、f1=±3 585 Hz等10條頻譜直線，但卻不能識別軸承內(nèi)圈故障特征頻率。在圖21、圖22中，因小波濾波可視為帶通濾波，經(jīng)過小波閾值濾波后，去掉了信號的高頻部分，只保留了信號中的中、低頻部分，可模糊看出f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz等8條頻譜直線，同樣也不能識別軸承內(nèi)圈故障特征頻率。在上述頻譜直線的交點處出現(xiàn)譜峰，表明軸承內(nèi)圈故障振動信號具有很強的非線性耦合，但軸承內(nèi)圈故障振動信號的雙譜和小波閾值降噪后雙譜已完全被強噪聲淹沒，因此，難以準確分辨軸承內(nèi)圈故障特征頻率及其高次諧波。

圖19 軸承內(nèi)圈故障信號的傳統(tǒng)雙譜Fig.19 Bi-spectrum of bearing inner race fault vibration signal

圖20 軸承內(nèi)圈故障信號的傳統(tǒng)雙譜Fig.20 Bi-spectrum of bearing inner race fault vibration signal

圖21 軸承內(nèi)圈故障信號小波閾值降噪后的雙譜Fig.21 Wavelet de-noising based bi-spectrum of bearing inner race fault vibration signal

圖22 軸承內(nèi)圈故障信號小波閾值降噪后的雙譜Fig.22 Wavelet de-noising based bi-spectrum of bearing inner race fault vibration signal

5.2 軸承外圈故障診斷

圖23為滾動軸承外圈故障時采集的振動信號及其FFT。在圖23(a)中存在明顯的幅值調(diào)制現(xiàn)象，在圖23(b)中的頻率等于3 445 Hz處存在顯著的峰值，但根據(jù)時域波形及其FFT還不能準確識別軸承故障類型。

(a) 時域圖

(b) 頻域圖圖23 軸承外圈故障振動信號的時域圖和頻譜圖Fig.23 Bearing outer race fault vibration signal and its FFT

圖24 矩陣各列的Gini指數(shù)Fig.24 Gini index of each column for matrix

圖25 軸承外圈故障振動信號的相關(guān)熵Fig.25 Correntropy of bearing outer race fault vibration signal

圖26和圖27分別為軸承外圈故障振動信號相關(guān)熵(σ=3)的輪廓圖和三維圖。從圖26和圖27可以清楚地看到：在f1=±3 445 Hz、f2=±3 445 Hz等4條頻譜直線的交點附近，存在6個明顯的譜峰簇，且譜峰之間的間隔為軸承外圈故障特征頻率，這種頻譜結(jié)構(gòu)表明了軸承的外圈故障特征，同時也說明電動機系統(tǒng)在3 445 Hz頻率處具有較大的振動能量，與圖23(b)FFT頻譜圖的結(jié)果一致。

圖26 軸承外圈故障振動信號相關(guān)熵的雙譜Fig.26 Correntropy based bi-spectrum of bearing outer race fault vibration signal

圖27 軸承外圈故障信號相關(guān)熵的雙譜Fig.27 Correntropy based bi-spectrum of bearing outer race fault vibration signal

為了與本文提出的方法進行對比，圖28和圖29給出了利用直接法計算的軸承外圈故障振動信號的雙譜，圖30和圖31為采用小波閾值降噪后信號的雙譜。在圖28和圖29中，因噪聲的影響，可模糊看出f1=±3 445 Hz、f2=±3 445 Hz等頻譜線，但卻不能清晰識別以軸承外圈故障特征頻率為間隔的邊頻帶簇。在圖30和圖31中，小波閾值降噪后信號的雙譜也難以正確分辨軸承外圈故障特征頻率。

圖28 軸承外圈故障信號的傳統(tǒng)雙譜Fig.28 Bi-spectrum of bearing outer race fault vibration signal

圖29 軸承外圈故障信號的傳統(tǒng)雙譜Fig.29 Bi-spectrum of bearing outer race fault vibration signal

圖30 軸承外圈故障信號小波閾值降噪后的雙譜Fig.30 Wavelet de-noising based bi-spectrum of bearing outer race fault vibration signal

圖31 軸承外圈故障信號小波閾值降噪后的雙譜Fig.31 Wavelet de-noising based bi-spectrum of bearing outer race fault vibration signal

通過上述分析可以看出：基于相關(guān)熵的雙譜具有很強的從強噪聲背景中提取軸承故障特征的能力，在雙譜圖中能清晰地刻畫軸承內(nèi)圈、外圈故障特征，提高了軸承故障診斷的可靠性和準確性，為從強高斯和非高斯噪聲環(huán)境中提取軸承故障特征的有效方法。

6 結(jié) 論

傳統(tǒng)雙譜具有較強的抑制高斯噪聲的能力，但卻不能有效抑制非高斯噪聲的干擾。相關(guān)熵不僅能有效抑制高斯噪聲，而且能有效抑制非高斯噪聲，因此，相關(guān)熵為高斯、非高斯噪聲的處理提供了一種嶄新的魯棒性方法。綜合利用了相關(guān)熵和雙譜的優(yōu)點，提出了基于相關(guān)熵雙譜分析的軸承故障診斷方法。首先利用不完全Cholesky分解算法，計算核矩陣的低秩下三角矩陣，再利用Gini指數(shù)選擇下三角矩陣的主分量，不僅降低了運算量，提高了核矩陣的計算速度，而且突出了軸承故障的瞬態(tài)沖擊特征，提高了軸承故障診斷的準確性和可靠性。仿真和試驗結(jié)果表明：基于相關(guān)熵和雙譜分析技術(shù)，能有效提取淹沒在強高斯和非高斯噪聲環(huán)境中的微弱信號，提高了信噪比，為軸承故障診斷的有效方法。