MBCV-EWT和奇異值差分譜的滾動軸承信號降噪方法*

王亞萍， 崔 巍， 葛江華， 許 迪， 李云飛

(哈爾濱理工大學機械動力工程學院 哈爾濱，150080)

引 言

當機械設備運行時，由于現場環境及工況變化等因素的影響會導致信號中會摻雜大量的噪聲，同時也會引起動態信號出現非平穩性。這些非平穩信號的統計量是時變函數，需要通過時頻分析方法處理信號，常用的有Wigner-Ville分布、經驗模態分解和小波變換[1-3]等，但是以上方法存在欠包絡、過包絡和模式混疊等問題。

為了克服傳統信號分解中的問題，Gilles[4]提出經驗小波變換(empirical wavelet transform，簡稱EWT)方法，可以自適應劃分初始信號的頻譜，同時通過相對應的帶通濾波器在各自的劃分區間內構造正交帶通濾波器組，以此提取具有緊支撐的調幅-調頻分量。Aneesh等[5]分別用變分模態分解與EWT對電能質量進行分類并做對比。Thirumala等[6]將EWT用于電能質量指標(power-quality indices，簡稱PQIS)的估計。為了驗證EWT與傳統方法的差異，國內學者將EWT應用于轉子，并與經驗模態分解(empirical mode decomposition，簡稱EMD)、集合經驗模態分解(ensemble empirical mode decomposition，簡稱EEMD)等方法進行對比，證明了EWT的有效性[7-10]。EWT方法中最核心的部分是信號頻譜的劃分，劃分后的區間影響后續信號分解的效果，文獻[11]對經典劃分方法進行了闡述。祝文穎等[12]提出了一種單分量個數的估算方法用于頻譜劃分，給出了敏感信號分量選取方法，仍存在EWT區間定位模糊的問題。劃分后的調幅-調頻分量較多，缺乏尋優過程，需要建立一個篩選指標來挑選最優的分量，而無量綱指標基本不受工況、載荷和轉速等變化影響的特點為此提供了新的思路[13]。

為了進一步實現故障特征的提取，需要先對噪聲干擾進行抑制。Jumah等[14]提出基于小波變換系數取閾值的方法，該方法對去除一維高斯白噪聲具有較好效果，但不能改進小波變換理論不足的缺點。Imaouchen等[15]將互補集合經驗模態分解用于強噪聲背景下信號的處理實現強制降噪，模態混疊問題依舊存在。奇異值差分譜可以識別出奇異值的最大突變點位置，通過對突變點前一系列成分的重構可以實現信號的降噪，該方法在處理分解后的分量時效果突出[16-18]。胥永剛等[19]提出雙樹復小波的信號分解方法，通過奇異值差分譜實現單分量的降噪，但其計算復雜，很難應用到實際生產生活中。

綜上所述，針對傳統EWT存在的問題，筆者提出MBCV-EWT與奇異值差分譜相結合的信號降噪方法，通過自適應的對信號頻率進行劃分克服模式混疊等問題，利用評價指標篩選最優調幅-調頻分量并根據最大奇異值突變點實現信號的降噪。

1 最大類間方差-經驗小波變換分解

經驗小波變換的實質是對信號的頻譜進行劃分，并相應的在每個區間構造帶通濾波器，實現對信號中不同的調幅-調頻成分的提取，為后續降噪提供基礎。

1.1 多域多類別故障特征提取

最大類間方差(maximum between-cluster variance，簡稱MBCV)[20-22]通過計算目標與背景的類間方差并將最大值作為圖像劃分閾值的標準。對于給定圖像，假設像素數為N，灰度取[0,L-1] ，則

(1)

其中：ni為素數；pi為像素點出現的概率，pi=ni/N,i=0,1,…,L-1。

圖像根據閾值K分為目標c0和非目標c1兩部分，分別為灰度值在[0,k]和[k+1,L-1]中的像素構成，圖像的均值表示為

(2)

c0與c1的均值為

(3a)

(3b)

類間方差表示為

w1(c1-ck)2=w0w1(c0-c1)2

(4)

為了提高MBCV的自適應性，引入P作為循環結束的判定準則確定閾值的數量，P的取值范圍為[0，1] ，即

(5)

P越大，表明類間的差別越大，當P近似于1時，類間方差為最大值，可以確定閾值數，自適應的區間就被劃分出來。

1.2 構建區間內的帶通濾波器

確定了每個區間的邊界之后，模仿小波變換的方式在區間上構造帶通濾波器，根據2π的周期性只研究區間[0，2π]，可以得到

(6)

在Tn不重疊的情況下也同樣適用于τn，即

τn+τn+1<ωn+1-ωn?γωn+γωn+1<

(7)

經驗小波變換的細節系數為

(8)

經驗小波變換的近似系數為

(9)

得到經驗模態函數為

(10a)

(10b)

1.3 調幅-調頻分量篩選指標

有量綱分析的故障診斷主要是從時域方面進行分析，主要包含均值、均方根值和方差3種幅域參數，有量綱指標對故障程度反映較為敏感，通常用于中度至重度損傷的故障診斷與預測，受外部環境影響較大，且易受復合故障的影響。

無量綱指標是根據有量綱指標的比值轉化而來，反應頻率密度函數的形狀，對于時間序列信號X={xi},i=1,2,…,n，筆者應用無量綱幅域參數，例如：波形指標Sf、峰值指標Cf、脈沖指標If、裕度指標CLf和峭度指標Kr。

無量綱的優點可歸納為：a.完全具有反映故障特征的能力；b.幾乎與振動信號絕對水平無關；c.對不同類型故障存在不同敏感性；d.對復合并發故障不敏感；e.基本不受工況、載荷和轉速等變化的影響。

2 基于奇異值差分譜的調幅-調頻分量降噪

MBCV-EWT分解后得到的分量通過脈沖指標選取最優的調幅-調頻分量，并引入奇異值差分譜降噪方法。令信號矩陣為X=[x(1),x(2),…,x(N)]，對其進行奇異值分解及奇異值差分譜去噪。

構造Hankel矩陣為

(11)

其中：1

令m=N-n+1，A∈Rm×n，則矩陣A為Hankel矩陣，即重構吸引子軌道矩陣?？梢钥闯? Hankel矩陣相鄰兩行的矢量存在一定聯系，后一行滯后一個位置，有用信號矩陣的奇異值特點為數值較大的集中在前端，矩陣秩對應點處差距巨大，之后的奇異值趨近于零。當矩陣中存在零奇異值，則該矩陣必然不滿秩，即為奇異矩陣，其誤差矩陣E的范數為零，而噪聲矩陣為滿秩均值，行矢量互不相關。

對于吸引子軌道矩陣A∈Rm×n，無論矩陣的行和列是否相關，必然存在正交矩陣U=(u1,u2,L,um)和V=(v1,v2,L,vn)，滿足

A=USVT

(12)

其中：U和V分別為矩陣A的左右奇異陣；S=(diag(σ1,σ2,L,σq),0)或其轉置，這由mn決定；A∈Sm×n，0為零矩陣；q=min(m,n)，σ1≥σ2≥…≥σq≥0，為矩陣A的奇異值。

從本質上來說，奇異值分解就是將調幅調頻分量分解成一組成分的簡單線性疊加，成分之間不存在相位差，有用信息與噪聲分別儲存在不同的成分上，通過篩選成分進行重構可以實現有用信息的保留和噪聲成分的去除。

集合B=[b1,b2,…,bq-1]為奇異值差分譜序列，其中，bi為相鄰奇異值的差值。選擇最大突變點對調幅調頻分量進行重構，根據矩陣特點將其按行與列收尾連接并進行逆變換，實現調幅調頻分量的重構，重構后的信號即為降噪后的有用信號。

圖1 降噪流程圖Fig.1 Noise reduction flowchart

定義奇異值差分譜為

bi=σi-σi+1

(13)

可以看出，集合中必然存在一個最大的峰值bk滿足最大突變點的要求。最大突變點的本質是源于不同信號奇異值差別較大，有用信號通常處于前k個奇異值成分，而噪聲在之后的成分中，因此最大突變點的位置就是重構調幅調頻分量所需的分量數。降噪過程如圖1所示。

3 仿真與實驗

3.1 仿真驗證

仿真實驗驗證了該方法的有效性，仿真信號分別由頻率為5 Hz的信號x1、20 Hz的信號x2、100 Hz的信號x3和高斯白噪聲n組成，信號的采樣頻率fs=1 024 Hz。仿真信號表達式為

(14)

原始信號幅值圖及其頻譜圖，噪聲信號的幅值圖及其頻譜圖分別如圖2，3所示。

圖2 原始信號時域和頻域圖Fig.2 Time domain, frequency domain diagram of the original signal

圖3 加噪信號時域和頻域圖Fig.3 Time domain and frequency domain diagram of the noisy signal

對于傳統EWT方法，頻譜的劃分必須計算出每個邊界的具體位置。首先，求出頻譜幅值的極大值，如圖4所示。

圖4 幅值極大值點Fig.4 Amplitude maximum point

根據Mm+α(M1-Mm)

圖5 EWT自適應頻譜劃分Fig.5 EWT adaptive spectrum division

采用MBCV對頻譜進行處理，頻譜劃分如圖6所示?？梢钥闯觯琈BCV對于頻譜的劃分與傳統方法有所區別，在5 Hz和20 Hz，20 Hz和100 Hz之間有且僅有一條邊界，說明此方法可以實現準確的劃分。

圖6 MBCV-EWT頻譜劃分Fig.6 MBCV-EWT spectrum division

此外，MBCV方法并不需要計算準確的邊界，根據劃分區間的個數，MBCV方法可以自適應地實現邊界的確定。如圖7所示，當設定只劃分3個區域的時候， MBCV分解后正是信號所含有的3個主頻。

對于較為復雜的信號，仿真信號表達式為

x=sin(20πt)sin(100πt+cos(20πt))+

1.5cos(4πt)cos(200πt)+sin(10πt)

(15)

如圖8所示，對該信號進行MBCV-EWT分解，可以看出， 第1個區間內的成分即為正弦信號， 第2個區間為調幅調頻信號，第3個區間為調幅信號。如圖9所示，當正弦信號的頻率改為125 Hz時，正弦信號的主頻在調幅調頻信號頻率之間，MBCV-EWT可以將該頻率與其他頻率區分開。

圖7 N=3時的MBCV-EWT分解Fig.7 MBCV-EWT decomposition at N=3

圖8 復雜信號的MBCV-EWT分解Fig.8 MBCV-EWT decomposition of complex signals

圖9 正弦信號改變后分解圖Fig.9 Decomposed graph after sinusoidal signal changed

采用MBCV-EWT， 完備總體經驗模態分解(complete ensemble emprirical mode decomposition,簡稱CEEMD)和雙層小波變換對信號進行分解，如圖10～12所示。

圖10 MBCV-EWT分解圖Fig.10 MBCV-EWT exploded view

圖11 CEEMD分解圖Fig.11 CEEMD exploded view

圖12 雙層小波變換Fig.12 Double layer wavelet transform

通過對比發現，MBCV-EWT分解后的信號與小波變換、CEEMD分解相比得到的分量更少，每個調幅-調頻分量中僅含有一個主頻成分，不同的分量之間不存在模態混疊現象，更不存在固定小波后再分解導致的窗固定問題。從Matlab計算速度來看，對同一個信號MBCV-EWT的速度要優于其他幾種。

評價指標部分需要采用真實的全壽命周期實驗，這里采用辛辛那提大學的實驗數據。試驗臺上安裝4個雙列滾動軸承，實驗時間為2003-10-22T 12∶06∶24～2003-11-25T 23∶59∶56，共進行了34 d，每10 min采集一次，在實驗的最后階段軸承3發生了內圈故障。通過時間變化對每個指標的表現進行計算，如圖13所示，脈沖指標前26 d表現平穩，從第26 d起，兩個指標的幅值都產生明顯上升，說明對早期故障非常敏感，同時在失效之前都能保持較為明顯的波動。

對于EWT分解后的3個調幅-調頻分量，以第1個調幅-調頻分量為例，其時域波形圖和頻譜圖如圖14所示。

圖13 脈沖指標隨時間變化曲線圖Fig.13 Pulse indicator over time curve

圖14 AM-FM1時域、頻域圖Fig.14 AM-FM1 time domain, frequency domain diagram

以第1個分量為例構建Hankel矩陣，畫出奇異值分解圖15(a)和差分譜15(b)，取前50個點繪制在一起如圖15(c)所示。

圖15 奇異值差分譜降噪過程Fig.15 Singular value difference spectrum de-noising process

從圖15(c)可以看到，最大的峰值發生在第2個點，表明奇異值序列在此位置發生了最大的突變。取前兩個分量進行重構，如圖16所示。

圖16 AM-FM1降噪時域、頻譜圖Fig.16 Time domain, frequency domain diagram of the AM-FM1 after de-noising

同理可以得到AM-FM2和AM-FM3的降噪時域、頻譜圖，如圖17,18所示。

圖17 AM-FM2降噪時域、頻譜圖Fig.17 Time domain, frequency domain diagram of the AM-FM2 after de-noising

圖18 AM-FM3降噪時域、頻譜圖Fig.18 Time domain, frequency domain diagram of the AM-FM3 after de-noising

為驗證降噪方法的有效性，將奇異值差分譜降噪分別與小波閾值降噪、CEEMD強制降噪進行對比，如圖19所示?？梢钥闯?，奇異值差分譜降噪后的信號更接近原始信號。

圖19 降噪效果對比Fig.19 Contrast of de-noising effect

圖20 滾動軸承振動測試試驗臺Fig.20 Rolling bearing vibration test rig

采用3種方法降噪后，分別計算重構后信號與原始信號的相關性，相關性越高，說明降噪方法對信號的還原度越高。MBCV-EWT與奇異值差分譜降噪后信號與原始信號的相關系數為0.936 1；CEEMD的相關系數為0.438 2；小波變換的相關系數為0.488 5。可見，MBCV-EWT與奇異值差分譜相結合的降噪方法明顯優于前兩種。

3.2 實驗驗證

本論文主要研究早期故障階段，采用如圖20所示的試驗臺模擬實驗來驗證筆者提出的微弱信號檢測、信號分解、降噪與特征提取方法的有效性。滾動軸承振動測試試驗臺主要由SGM7J-04AFC6S伺服電機、YMC122A100加速度傳感器、POD-0.6 kg磁粉制動器、GFC-40X66梅花聯軸器和底座等連接件與緊固件組成。

該實驗用CoCo-80動態信號分析儀采集信號數據，3204ATN雙列角接觸球軸承，參數如表1所示。根據滾動軸承常見故障位置與故障類型，實驗主要模擬外圈裂痕故障。外圈裂痕情況如圖21所示，圖中圓圈即為裂痕位置，裂痕為電火花線切割的0.31 mm的輕度故障。轉速n=1 kr/min，采樣頻率f=400 Hz，通過計算得到外圈的故障頻率約為 52.98 Hz。

表1 具體參數Tab.1 Specific parameters

圖21 裂痕故障程度實物圖Fig.21 Physical map of crack faults

滾動軸承外圈故障頻率計算公式為

(16)

數據傅里葉變換后，裂痕故障信號與原始信號如圖22所示。可以看出，裂痕故障對信號的影響比較明顯，信號的轉頻幅值較低。從信號頻譜中可以看到，在50～60 Hz的范圍內有一個波峰，其值與故障頻率接近，同時根據式(16)可以求得該故障頻率所對應的即為軸承外圈故障。

圖22 信號對比Fig.22 Signal comparison chart

采用MBCV-EWT對信號進行分解，如圖23所示。故障頻率在第3個調幅-調頻分量中，并對分解得到的調幅-調頻分量進行奇異值差分譜降噪，通過圖24可以看出降噪后的信號具有較明顯的周期成分。

圖23 EWT自適應頻譜劃分Fig.23 EWT adaptive spectrum division

圖24 降噪后AM-FM3頻譜圖Fig.24 AM-FM3 frequency spectrum after de-noising

4 結 論

1) MBCV-EWT相對于傳統信號分解方法，能夠自適應地將信號頻譜劃分成區間，克服了模式混疊，每個分解得到的調幅-調頻分量只對應一個頻率，分解準確性高。

2) 通過評價指標對各調幅-調頻分量進行篩選，可以得到相關性最高的一組，有效減少后續降噪的計算量。

3) 基于奇異值差分譜的降噪能夠將調幅-調頻分量中的主要頻率識別出來，抑制多余噪聲，降噪效果明顯。

4) 通過仿真和實驗數據的驗證，該方法能夠有效地對振動信號進行自適應分解和降噪，得到具體的故障頻率，實現故障位置的識別或預測。