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

(哈爾濱理工大學(xué)機(jī)械動(dòng)力工程學(xué)院 哈爾濱，150080)

引 言

當(dāng)機(jī)械設(shè)備運(yùn)行時(shí)，由于現(xiàn)場(chǎng)環(huán)境及工況變化等因素的影響會(huì)導(dǎo)致信號(hào)中會(huì)摻雜大量的噪聲，同時(shí)也會(huì)引起動(dòng)態(tài)信號(hào)出現(xiàn)非平穩(wěn)性。這些非平穩(wěn)信號(hào)的統(tǒng)計(jì)量是時(shí)變函數(shù)，需要通過時(shí)頻分析方法處理信號(hào)，常用的有Wigner-Ville分布、經(jīng)驗(yàn)?zāi)B(tài)分解和小波變換[1-3]等，但是以上方法存在欠包絡(luò)、過包絡(luò)和模式混疊等問題。

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

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

綜上所述，針對(duì)傳統(tǒng)EWT存在的問題，筆者提出MBCV-EWT與奇異值差分譜相結(jié)合的信號(hào)降噪方法，通過自適應(yīng)的對(duì)信號(hào)頻率進(jìn)行劃分克服模式混疊等問題，利用評(píng)價(jià)指標(biāo)篩選最優(yōu)調(diào)幅-調(diào)頻分量并根據(jù)最大奇異值突變點(diǎn)實(shí)現(xiàn)信號(hào)的降噪。

1 最大類間方差-經(jīng)驗(yàn)小波變換分解

經(jīng)驗(yàn)小波變換的實(shí)質(zhì)是對(duì)信號(hào)的頻譜進(jìn)行劃分，并相應(yīng)的在每個(gè)區(qū)間構(gòu)造帶通濾波器，實(shí)現(xiàn)對(duì)信號(hào)中不同的調(diào)幅-調(diào)頻成分的提取，為后續(xù)降噪提供基礎(chǔ)。

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

最大類間方差(maximum between-cluster variance，簡(jiǎn)稱MBCV)[20-22]通過計(jì)算目標(biāo)與背景的類間方差并將最大值作為圖像劃分閾值的標(biāo)準(zhǔn)。對(duì)于給定圖像，假設(shè)像素?cái)?shù)為N，灰度取[0,L-1] ，則

(1)

其中：ni為素?cái)?shù)；pi為像素點(diǎn)出現(xiàn)的概率，pi=ni/N,i=0,1,…,L-1。

圖像根據(jù)閾值K分為目標(biāo)c0和非目標(biāo)c1兩部分，分別為灰度值在[0,k]和[k+1,L-1]中的像素構(gòu)成，圖像的均值表示為

(2)

c0與c1的均值為

(3a)

(3b)

類間方差表示為

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

(4)

為了提高M(jìn)BCV的自適應(yīng)性，引入P作為循環(huán)結(jié)束的判定準(zhǔn)則確定閾值的數(shù)量，P的取值范圍為[0，1] ，即

(5)

P越大，表明類間的差別越大，當(dāng)P近似于1時(shí)，類間方差為最大值，可以確定閾值數(shù)，自適應(yīng)的區(qū)間就被劃分出來。

1.2 構(gòu)建區(qū)間內(nèi)的帶通濾波器

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

(6)

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

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

(7)

經(jīng)驗(yàn)小波變換的細(xì)節(jié)系數(shù)為

(8)

經(jīng)驗(yàn)小波變換的近似系數(shù)為

(9)

得到經(jīng)驗(yàn)?zāi)B(tài)函數(shù)為

(10a)

(10b)

1.3 調(diào)幅-調(diào)頻分量篩選指標(biāo)

有量綱分析的故障診斷主要是從時(shí)域方面進(jìn)行分析，主要包含均值、均方根值和方差3種幅域參數(shù)，有量綱指標(biāo)對(duì)故障程度反映較為敏感，通常用于中度至重度損傷的故障診斷與預(yù)測(cè)，受外部環(huán)境影響較大，且易受復(fù)合故障的影響。

無量綱指標(biāo)是根據(jù)有量綱指標(biāo)的比值轉(zhuǎn)化而來，反應(yīng)頻率密度函數(shù)的形狀，對(duì)于時(shí)間序列信號(hào)X={xi},i=1,2,…,n，筆者應(yīng)用無量綱幅域參數(shù)，例如：波形指標(biāo)Sf、峰值指標(biāo)Cf、脈沖指標(biāo)If、裕度指標(biāo)CLf和峭度指標(biāo)Kr。

無量綱的優(yōu)點(diǎn)可歸納為：a.完全具有反映故障特征的能力；b.幾乎與振動(dòng)信號(hào)絕對(duì)水平無關(guān)；c.對(duì)不同類型故障存在不同敏感性；d.對(duì)復(fù)合并發(fā)故障不敏感；e.基本不受工況、載荷和轉(zhuǎn)速等變化的影響。

2 基于奇異值差分譜的調(diào)幅-調(diào)頻分量降噪

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

構(gòu)造Hankel矩陣為

(11)

其中：1

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

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

A=USVT

(12)

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

從本質(zhì)上來說，奇異值分解就是將調(diào)幅調(diào)頻分量分解成一組成分的簡(jiǎn)單線性疊加，成分之間不存在相位差，有用信息與噪聲分別儲(chǔ)存在不同的成分上，通過篩選成分進(jìn)行重構(gòu)可以實(shí)現(xiàn)有用信息的保留和噪聲成分的去除。

集合B=[b1,b2,…,bq-1]為奇異值差分譜序列，其中，bi為相鄰奇異值的差值。選擇最大突變點(diǎn)對(duì)調(diào)幅調(diào)頻分量進(jìn)行重構(gòu)，根據(jù)矩陣特點(diǎn)將其按行與列收尾連接并進(jìn)行逆變換，實(shí)現(xiàn)調(diào)幅調(diào)頻分量的重構(gòu)，重構(gòu)后的信號(hào)即為降噪后的有用信號(hào)。

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

定義奇異值差分譜為

bi=σi-σi+1

(13)

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

3 仿真與實(shí)驗(yàn)

3.1 仿真驗(yàn)證

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

(14)

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

圖2 原始信號(hào)時(shí)域和頻域圖Fig.2 Time domain, frequency domain diagram of the original signal

圖3 加噪信號(hào)時(shí)域和頻域圖Fig.3 Time domain and frequency domain diagram of the noisy signal

對(duì)于傳統(tǒng)EWT方法，頻譜的劃分必須計(jì)算出每個(gè)邊界的具體位置。首先，求出頻譜幅值的極大值，如圖4所示。

圖4 幅值極大值點(diǎn)Fig.4 Amplitude maximum point

根據(jù)Mm+α(M1-Mm)

圖5 EWT自適應(yīng)頻譜劃分Fig.5 EWT adaptive spectrum division

采用MBCV對(duì)頻譜進(jìn)行處理，頻譜劃分如圖6所示。可以看出，MBCV對(duì)于頻譜的劃分與傳統(tǒng)方法有所區(qū)別，在5 Hz和20 Hz，20 Hz和100 Hz之間有且僅有一條邊界，說明此方法可以實(shí)現(xiàn)準(zhǔn)確的劃分。

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

此外，MBCV方法并不需要計(jì)算準(zhǔn)確的邊界，根據(jù)劃分區(qū)間的個(gè)數(shù)，MBCV方法可以自適應(yīng)地實(shí)現(xiàn)邊界的確定。如圖7所示，當(dāng)設(shè)定只劃分3個(gè)區(qū)域的時(shí)候， MBCV分解后正是信號(hào)所含有的3個(gè)主頻。

對(duì)于較為復(fù)雜的信號(hào)，仿真信號(hào)表達(dá)式為

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

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

(15)

如圖8所示，對(duì)該信號(hào)進(jìn)行MBCV-EWT分解，可以看出， 第1個(gè)區(qū)間內(nèi)的成分即為正弦信號(hào)， 第2個(gè)區(qū)間為調(diào)幅調(diào)頻信號(hào)，第3個(gè)區(qū)間為調(diào)幅信號(hào)。如圖9所示，當(dāng)正弦信號(hào)的頻率改為125 Hz時(shí)，正弦信號(hào)的主頻在調(diào)幅調(diào)頻信號(hào)頻率之間，MBCV-EWT可以將該頻率與其他頻率區(qū)分開。

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

圖8 復(fù)雜信號(hào)的MBCV-EWT分解Fig.8 MBCV-EWT decomposition of complex signals

圖9 正弦信號(hào)改變后分解圖Fig.9 Decomposed graph after sinusoidal signal changed

采用MBCV-EWT， 完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(complete ensemble emprirical mode decomposition,簡(jiǎn)稱CEEMD)和雙層小波變換對(duì)信號(hào)進(jìn)行分解，如圖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

通過對(duì)比發(fā)現(xiàn)，MBCV-EWT分解后的信號(hào)與小波變換、CEEMD分解相比得到的分量更少，每個(gè)調(diào)幅-調(diào)頻分量中僅含有一個(gè)主頻成分，不同的分量之間不存在模態(tài)混疊現(xiàn)象，更不存在固定小波后再分解導(dǎo)致的窗固定問題。從Matlab計(jì)算速度來看，對(duì)同一個(gè)信號(hào)MBCV-EWT的速度要優(yōu)于其他幾種。

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

對(duì)于EWT分解后的3個(gè)調(diào)幅-調(diào)頻分量，以第1個(gè)調(diào)幅-調(diào)頻分量為例，其時(shí)域波形圖和頻譜圖如圖14所示。

圖13 脈沖指標(biāo)隨時(shí)間變化曲線圖Fig.13 Pulse indicator over time curve

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

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

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

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

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

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

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

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

為驗(yàn)證降噪方法的有效性，將奇異值差分譜降噪分別與小波閾值降噪、CEEMD強(qiáng)制降噪進(jìn)行對(duì)比，如圖19所示。可以看出，奇異值差分譜降噪后的信號(hào)更接近原始信號(hào)。

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

圖20 滾動(dòng)軸承振動(dòng)測(cè)試試驗(yàn)臺(tái)Fig.20 Rolling bearing vibration test rig

采用3種方法降噪后，分別計(jì)算重構(gòu)后信號(hào)與原始信號(hào)的相關(guān)性，相關(guān)性越高，說明降噪方法對(duì)信號(hào)的還原度越高。MBCV-EWT與奇異值差分譜降噪后信號(hào)與原始信號(hào)的相關(guān)系數(shù)為0.936 1；CEEMD的相關(guān)系數(shù)為0.438 2；小波變換的相關(guān)系數(shù)為0.488 5。可見，MBCV-EWT與奇異值差分譜相結(jié)合的降噪方法明顯優(yōu)于前兩種。

3.2 實(shí)驗(yàn)驗(yàn)證

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

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

表1 具體參數(shù)Tab.1 Specific parameters

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

滾動(dòng)軸承外圈故障頻率計(jì)算公式為

(16)

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

圖22 信號(hào)對(duì)比Fig.22 Signal comparison chart

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

圖23 EWT自適應(yīng)頻譜劃分Fig.23 EWT adaptive spectrum division

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

4 結(jié) 論

1) MBCV-EWT相對(duì)于傳統(tǒng)信號(hào)分解方法，能夠自適應(yīng)地將信號(hào)頻譜劃分成區(qū)間，克服了模式混疊，每個(gè)分解得到的調(diào)幅-調(diào)頻分量只對(duì)應(yīng)一個(gè)頻率，分解準(zhǔn)確性高。

2) 通過評(píng)價(jià)指標(biāo)對(duì)各調(diào)幅-調(diào)頻分量進(jìn)行篩選，可以得到相關(guān)性最高的一組，有效減少后續(xù)降噪的計(jì)算量。

3) 基于奇異值差分譜的降噪能夠?qū)⒄{(diào)幅-調(diào)頻分量中的主要頻率識(shí)別出來，抑制多余噪聲，降噪效果明顯。

4) 通過仿真和實(shí)驗(yàn)數(shù)據(jù)的驗(yàn)證，該方法能夠有效地對(duì)振動(dòng)信號(hào)進(jìn)行自適應(yīng)分解和降噪，得到具體的故障頻率，實(shí)現(xiàn)故障位置的識(shí)別或預(yù)測(cè)。