龔?fù)稹≡瑫暂x　王細(xì)洋

1.華中科技大學(xué)，武漢，430074　　2.南昌航空大學(xué)，南昌，330063

?

最小熵反褶積的數(shù)學(xué)形態(tài)法在滾動(dòng)軸承故障特征提取中的應(yīng)用

龔?fù)?，2袁曉輝1王細(xì)洋2

1.華中科技大學(xué)，武漢，4300742.南昌航空大學(xué)，南昌，330063

針對(duì)強(qiáng)噪聲背景下滾動(dòng)軸承故障特征提取，提出了基于最小熵反褶積的數(shù)學(xué)形態(tài)法。該方法先應(yīng)用最小熵反褶積算法加強(qiáng)信號(hào)中的沖擊特性，再利用數(shù)學(xué)形態(tài)法進(jìn)行故障特征提取，其中選取具有雙向脈沖提取能力的DIF濾波器作為形態(tài)算子，并以峭度值作為結(jié)構(gòu)元素長(zhǎng)度選取依據(jù)。仿真信號(hào)和滾動(dòng)軸承的內(nèi)外故障實(shí)例分析表明該方法具有較好的特征提取效果。通過對(duì)比發(fā)現(xiàn):最小熵反褶積算法能夠增大信號(hào)中峭度值，有效加強(qiáng)信號(hào)脈沖特性。

最小熵反褶積；數(shù)學(xué)形態(tài)；故障診斷；滾動(dòng)軸承

0　引言

滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械中最常用的通用零部件之一，同時(shí)也是最易損的零件。根據(jù)統(tǒng)計(jì)，旋轉(zhuǎn)機(jī)械的故障中有30%是由軸承故障引起的，因此它的工作狀態(tài)對(duì)機(jī)器運(yùn)行具有巨大的影響[1]。當(dāng)故障點(diǎn)出現(xiàn)在軸承的內(nèi)外圈及滾動(dòng)體上時(shí)，周期性脈沖會(huì)產(chǎn)生。而且這些脈沖信號(hào)的頻寬范圍常常覆蓋了軸承和機(jī)械系統(tǒng)固有振動(dòng)頻率，易引起軸承甚至機(jī)械系統(tǒng)的共振，導(dǎo)致機(jī)器振動(dòng)和產(chǎn)生噪聲，甚至損壞。因此這些沖擊脈沖的成功提取成為軸承故障診斷的一種有效方法。

對(duì)于這些沖擊脈沖性的振動(dòng)信號(hào)分析有很多方法，如局部特征分解、經(jīng)驗(yàn)?zāi)B(tài)分解和小波變換等[2-6]。相比前述方法，一種非線性的數(shù)學(xué)形態(tài)技術(shù)近來得到了廣泛關(guān)注。文獻(xiàn)[7]應(yīng)用形態(tài)閉算子作為形態(tài)算子，以0.6倍的脈沖周期作為結(jié)構(gòu)元素長(zhǎng)度，成功提取出了軸承內(nèi)外圈故障。但這種方法在選取結(jié)構(gòu)元素長(zhǎng)度時(shí)需要依賴經(jīng)驗(yàn)知識(shí)。為了克服這種缺陷，Dong等[8]提出了一種基于信噪比的改進(jìn)形態(tài)分析方法。然而，在強(qiáng)噪聲背景下形態(tài)方法所能提取出的故障形態(tài)特征，即脈沖沖擊特性是非常有限的。因此在形態(tài)方法處理前，希望通過一種方法能夠簡(jiǎn)化信號(hào)，突出信號(hào)中的沖擊特性，這樣形態(tài)方法就能提取出更多的故障形態(tài)特征，便于故障診斷。因此，基于上述考慮本文提出了一種基于最小熵反褶積的數(shù)學(xué)形態(tài)法(mathematical morphology based on minimum entropy deconvolution, MMBMED)，并應(yīng)用于故障軸承的振動(dòng)信號(hào)特征提取中。該方法通過最小熵反褶積(minimum eneropy deconvolution, MED)方法進(jìn)行預(yù)處理，而后應(yīng)用數(shù)學(xué)形態(tài)方法提取故障特征，選取具有雙向脈沖提取能力的DIF濾波器作為形態(tài)算子，結(jié)構(gòu)元素長(zhǎng)度的確定通過峭度作為選取原則。仿真信號(hào)和滾動(dòng)軸承內(nèi)外圈故障分析表明，基于最小熵反褶積的數(shù)學(xué)形態(tài)方法具有較好的故障特征提取效果。

1　相關(guān)理論

1.1數(shù)學(xué)形態(tài)法

數(shù)學(xué)形態(tài)學(xué)的基本思想是用結(jié)構(gòu)簡(jiǎn)單、緊湊的結(jié)構(gòu)元素去量度和提取圖像中的對(duì)應(yīng)形狀，以達(dá)到對(duì)圖像分析和識(shí)別的目的。形態(tài)學(xué)的非線性特性可以簡(jiǎn)化圖像，在保持圖像本質(zhì)基礎(chǔ)上除去噪聲部分，最初是作為圖像處理工具的，后來才發(fā)展到一維信號(hào)處理。

數(shù)學(xué)形態(tài)學(xué)中包含4個(gè)基本形態(tài)運(yùn)算算子,即膨脹算子、腐蝕算子、開算子和關(guān)算子。它們的定義分別如下：

設(shè)f(n)為一維多值信號(hào)，s(n)為結(jié)構(gòu)元素。H、K分別為f(n)和s(n)的定義域，H=(0,1，…，n-1)，K=(0，1，…，m-1)，且m

(f⊕s)(n)=min(f(n-m)+s(n))

(1)

m∈{0,1,…，M-1}

(fΘs)(n)=min(f(n+m)-s(n))

(2)

m∈{0,1,…M-1}

根據(jù)上述定義，形態(tài)開運(yùn)算和形態(tài)閉運(yùn)算的定義分別為

(3)

(f?s)(n)=(f⊕sΘs)(n)

(4)

Nikolaou等[7]通過分析發(fā)現(xiàn)4種形態(tài)算子對(duì)脈沖信號(hào)的作用不同。如腐蝕運(yùn)算對(duì)負(fù)脈沖具有平滑作用，會(huì)減少脈沖信號(hào)中的正脈沖部分。膨脹運(yùn)算具有和腐蝕運(yùn)算相反的效果。閉算子能夠保留信號(hào)中的正脈沖信號(hào)，開算子能夠保留信號(hào)中的負(fù)脈沖信號(hào)。顯然，實(shí)際應(yīng)用中很難得到正負(fù)脈沖的先驗(yàn)知識(shí)。因此，常常用4個(gè)基本算子的算術(shù)組合作為新的形態(tài)算子。例如，形態(tài)閉和開的差分運(yùn)算得到的DIF(differencefilter)算子，其數(shù)學(xué)定義如下式所示：

(5)

該算子具有同時(shí)提取正負(fù)脈沖信號(hào)的特性，因此具有較好的脈沖提取能力。

1.2最小熵反褶積算法

Wiggins[9]首次提出了MED的概念，并應(yīng)用在地震信號(hào)處理當(dāng)中。求取MED是指尋求一個(gè)反轉(zhuǎn)濾波器的最優(yōu)系數(shù)集合，使得所恢復(fù)的輸出信號(hào)具有最大峭度。因?yàn)榍投仁且粋€(gè)對(duì)故障敏感的因子，反映了信號(hào)中的峰值(即脈沖屬性)，增大信號(hào)峭度值有助于軸承中的故障診斷。

圖1說明了MED算法原理。信號(hào)e(n)通過結(jié)構(gòu)濾波器后的輸出和噪聲v(n)共同作為一個(gè)反轉(zhuǎn)濾波器 f 的輸入，輸出為y(n)。為了處理后有更好的效果，我們希望反轉(zhuǎn)濾波器f的輸出y(n)與最初的輸入信號(hào)e(n)盡可能相似。雖然系統(tǒng)輸入e(n)是未知的，但希望這些就是脈沖序列。Endo等[10]總結(jié)了MED方法具體過程：

圖1　MED算法原理

(1)最小熵算法為長(zhǎng)度L的反轉(zhuǎn)濾波器 f找到一個(gè)系數(shù)集合，優(yōu)化輸出信號(hào)y(n)的峭度：

(6)

其中，[·]表示信號(hào)為離散信號(hào)。

(2)反轉(zhuǎn)濾波器的卷積表達(dá)式為

(7)

(3)當(dāng)目標(biāo)函數(shù)(式(6))出現(xiàn)最大峭度值時(shí)，即

(8)

(4)對(duì)式(7)兩邊求導(dǎo)可得

(9)

(5)由式(7)和式(9)可將式(8)寫成：

(10)

式(10)寫成矩陣形式如下：

b=Af

(11)

式中，A為輸出信號(hào)y(n)的自相關(guān)矩陣；b為濾波器的輸出立方與輸入信號(hào)的互相關(guān)加權(quán)和；f為反轉(zhuǎn)濾波器系數(shù)矩陣。

由式(11)通過迭代計(jì)算可以計(jì)算出反轉(zhuǎn)濾波器系數(shù)：

f=A-1b

(12)

由以上分析可以歸納出最小熵遞歸迭代算法步驟如下：

(1)初始化反轉(zhuǎn)濾波器的系數(shù)，設(shè)定初始狀態(tài)f(0)中元素全為1；

(2)通過輸出信號(hào)x(0)和反轉(zhuǎn)濾波器初始值f(0)(按式(7))計(jì)算y(0)；

(3)通過y(0)和x(0)(按式(11))計(jì)算列矢量b(1)(b(1)為第一次計(jì)算b的結(jié)果)；

(4)計(jì)算得到濾波器系數(shù)矩陣f(1)=A-1b(1)；

(5)計(jì)算預(yù)設(shè)定誤差。如果誤差小于設(shè)定值，結(jié)束計(jì)算。否則，返回步驟(3)。

2　信號(hào)仿真

仿真信號(hào)由諧波信號(hào)、脈沖信號(hào)和噪聲三個(gè)信號(hào)復(fù)合構(gòu)成，其數(shù)學(xué)表達(dá)式如下：

x(t)=x1+x2+x3

(13)

其中，諧波信號(hào)x1=cos(60t)+cos(130t)；x2為周期為32 Hz的指數(shù)衰減序列，該信號(hào)在一個(gè)周期內(nèi)的衰減函數(shù)為e-10tsin(10t)；x3為均值為0，標(biāo)準(zhǔn)偏差為1的高斯白噪聲。采樣頻率fs為2048 Hz,采樣時(shí)間為1 s。圖2為仿真信號(hào)的時(shí)域圖和頻譜圖。在圖2b中脈沖特征頻率32 Hz、64 Hz等成分由于諧波信號(hào)和噪聲抑制，沒有被提取出來，因此采用建議方法進(jìn)行分析。新方法處理過程為：①先采用MED方法預(yù)處理；②通過數(shù)學(xué)形態(tài)法提取特征，以式(5)中的DIF算子作為形態(tài)算子，選擇扁平型高度為0的結(jié)構(gòu)元素。對(duì)于其長(zhǎng)度，本文以脈沖重復(fù)周期T的0.1倍為間隔，構(gòu)造9個(gè)不同長(zhǎng)度的結(jié)構(gòu)元素，即 {0.1T，0.2T，0.3T，0.4T，0.5T，0.6T，0.7T，0.8T，0.9T}，分別計(jì)算9個(gè)不同長(zhǎng)度的峭度值，最后選取最大峭度的長(zhǎng)度作為結(jié)構(gòu)元素長(zhǎng)度。

(a)時(shí)域圖

(b)頻譜圖圖2　仿真信號(hào)

仿真信號(hào)經(jīng)過MED處理后，繼續(xù)進(jìn)行數(shù)學(xué)形態(tài)分析。計(jì)算9個(gè)不同長(zhǎng)度的結(jié)構(gòu)元素，發(fā)現(xiàn)長(zhǎng)度為0.1T時(shí)的峭度值7.7879為最大，因此長(zhǎng)度選擇0.1T，相應(yīng)結(jié)果的頻譜圖見圖3a。從圖3a中發(fā)現(xiàn)脈沖特征頻率32 Hz，以及它的諧波成分(64～160 Hz) 都被提取出來了，諧波成分和噪聲得到了較好的抑制。特征提取效果說明最小熵反褶積和數(shù)學(xué)形態(tài)法相結(jié)合能夠有效提取出脈沖特征成分，較好地抑制噪聲。

為了說明MED方法的作用，仿真信號(hào)將直接運(yùn)用新方法的第二步進(jìn)行處理與對(duì)比。計(jì)算發(fā)現(xiàn)在0.1T時(shí)的峭度值只有4.2503，其結(jié)果頻譜圖見圖3b。圖中只發(fā)現(xiàn)了幾個(gè)脈沖特征頻率(64 Hz,96 Hz,128 Hz和160 Hz)。與圖3a對(duì)比可知：仿真信號(hào)未經(jīng)過MED預(yù)處理，結(jié)構(gòu)元素長(zhǎng)度在0.1T時(shí)峭度值從4.2503增大到7.7879，峭度值明顯增大。圖3b雖然提取出了幾個(gè)脈沖特征頻率，但它們周圍的噪聲較大，而且脈沖特征頻率32 Hz沒有提取出來。因此通過這個(gè)對(duì)比說明MED算法能夠優(yōu)化信號(hào)的峭度，加強(qiáng)信號(hào)的脈沖特性。

(a)有MED處理

(b)無MED處理圖3　仿真信號(hào)頻譜圖對(duì)比

3　實(shí)例分析

3.1滾動(dòng)軸承外圈故障分析

實(shí)例分析的數(shù)據(jù)來自于美國凱西西楚大學(xué)軸承試驗(yàn)中心。實(shí)驗(yàn)中外圈故障測(cè)試軸承為SKF6205-2RSJEMSKF深溝球軸承。采樣頻率為12kHz，采樣點(diǎn)數(shù)為12288。轉(zhuǎn)速為1750r/min，由計(jì)算可得外圈故障通過頻率(theballpassingfrequencyouterrace,BPFO)約為105.1Hz。外圈故障原始振動(dòng)信號(hào)如圖4a所示。在圖4b中BPFO等特征成分沒有顯現(xiàn)，因此采用本文方法處理。

(a)時(shí)域圖

(b)頻譜圖圖4　外圈故障原始信號(hào)

外圈故障信號(hào)經(jīng)MED處理后，計(jì)算發(fā)現(xiàn)當(dāng)結(jié)構(gòu)元素長(zhǎng)度為0.9T時(shí)峭度值達(dá)到最大，為8.2639，因此結(jié)構(gòu)元素長(zhǎng)度為0.9T，結(jié)果如圖5所示。圖5b中清楚提取出了BPFO和它的幾個(gè)諧波成分(2BPFO，3BPFO，4BPFO，5BPFO和6BPFO)。此外轉(zhuǎn)速fr也被找到。上述分析表明新方法有效診斷出了滾動(dòng)軸承外圈故障。

(a)時(shí)域圖

(b)頻譜圖圖5　外圈故障的新方法分析結(jié)果

3.2滾動(dòng)軸承內(nèi)圈故障分析

為了進(jìn)一步說明新方法對(duì)軸承內(nèi)圈故障的有效性，以一個(gè)內(nèi)圈故障為例進(jìn)行分析。測(cè)試軸承型號(hào)為6203-2RSJEMSKF。采樣頻率為12kHz，采樣點(diǎn)數(shù)為12 288，轉(zhuǎn)速為1772r/min。計(jì)算可知外圈故障通過頻率(theballpassingfrequencyinnerrace,BPFI)約為146Hz。內(nèi)圈故障原始振動(dòng)信號(hào)如圖6a所示，其頻普?qǐng)D如圖6b所示。

(a)時(shí)域圖

(b)頻譜圖圖6　內(nèi)圈故障原始信號(hào)

新方法處理后計(jì)算發(fā)現(xiàn)結(jié)構(gòu)元素長(zhǎng)度為0.6T時(shí)峭度值最大，為14.7514，因此結(jié)構(gòu)元素長(zhǎng)度為0.6T，相應(yīng)結(jié)果如圖7所示。在圖7b頻譜圖中，內(nèi)圈故障特征頻率BPFI成為了主要成分，它的三個(gè)諧波成分(2BPFI,3BPFI和4BPFI)也能夠識(shí)別出來，同時(shí)還能發(fā)現(xiàn)轉(zhuǎn)速 fr和它的2倍頻成分。此外特征頻率BPFI及它的諧波兩側(cè)還清晰找到了豐富的邊頻帶，即頻率成分為146κ±30i(其中，i=1,2;κ=1,2,3,4)。因此，這些結(jié)果說明新方法對(duì)滾動(dòng)軸承內(nèi)圈故障的特征提取也是有效的。

(a)時(shí)域圖

(b)頻譜圖圖7　內(nèi)圈故障的新方法分析結(jié)果

3.3有無最小熵反褶積方法的軸承內(nèi)外圈故障對(duì)比

與仿真信號(hào)中的對(duì)比相同，軸承內(nèi)外圈故障原始振動(dòng)信號(hào)直接采用數(shù)學(xué)形態(tài)法處理，外圈和內(nèi)圈故障時(shí)結(jié)構(gòu)元素長(zhǎng)度仍然是0.9T(峭度值2.9685)和0.6T(峭度值6.0521)，相應(yīng)頻譜圖分別為圖8a和圖8b。圖8a外圈故障中只有BPFO和2BPFO被提取出來。圖8b內(nèi)圈故障中，也只有BPFI和2BPFI被識(shí)別，而且特征頻率兩側(cè)的邊頻帶也未發(fā)現(xiàn)。與圖5b和圖7b相比：經(jīng)過MED方法處理后，外圈故障時(shí)峭度值從2.9685增大到8.2639，所提取BPFO的諧波成分更多(從BPFO～6BPFO)。內(nèi)圈故障時(shí)峭度值從6.0521增大到14.7514，所提取BPFI的諧波成分更多(從BPFI～4BPFO)，以及它們兩側(cè)邊頻帶也較為明顯。因此，事實(shí)說明MED算法的峭度值增大效果明顯，有效加強(qiáng)了軸承內(nèi)圈故障信號(hào)的脈沖特征，特征提取效果更好。

(a)外圈無MED分析

(b)內(nèi)圈無MED分析圖8　軸承外圈和內(nèi)圈故障時(shí)的頻譜圖對(duì)比

4　結(jié)語

為了提取出強(qiáng)噪聲背景下滾動(dòng)軸承內(nèi)外圈故障的特征頻率，提出了基于最小熵反褶積的數(shù)學(xué)形態(tài)法。新方法利用MED方法進(jìn)行預(yù)處理，以提高原始信號(hào)峭度值，加強(qiáng)故障的沖擊特性。然后運(yùn)用數(shù)學(xué)形態(tài)方法提取故障特征。選取了DIF濾波器作為形態(tài)算子，結(jié)構(gòu)元素長(zhǎng)度以9個(gè)不同長(zhǎng)度中峭度值作為選取依據(jù)。通過仿真信號(hào)和滾動(dòng)軸承內(nèi)外圈故障分析說明新方法能夠提取出豐富的故障特征。通過對(duì)比發(fā)現(xiàn)最小熵反褶積使得信號(hào)峭度值增大明顯，軸承故障信號(hào)沖擊特征得到了加強(qiáng)，故障特征提取效果明顯。

[1]McFaddenPD,SmithJD.VibrationMonitoringofRollingElementBearingsbytheHigh-frequencyResonanceTechnique—AReview[J].TribologyInternational, 1984, 17:3-10.

[2]GeorgoulasG,LoutasT,StyliosCD,etal.BearingFaultDetectionBasedonHybridEnsembleDetectorandEmpiricalModeDecomposition[J].MechanicalSystemsandSignalProcessing, 2013, 41:510-525.

[3]GuoW,TsePW,ANovelSignalCompressionMethodBasedonOptimalEnsembleEmpiricalModeDecompositionforBearingVibrationSignals[J].JournalofSoundandVibration, 2013, 332:423-441.

[4]QiuH,LeeJ,LinJ,etal,WaveletFilter-basedWeakSignatureDetectionMethodandItsApplicationonRollingElementBearingPrognostics[J].JournalofSoundandVibration, 2006,289:1066-1090.

[5]楊宇, 程軍圣. 基于LCD降噪和VPMCD的滾動(dòng)軸承故障診斷方法[J]. 中國機(jī)械工程, 2013, 24(24):3338-3344.

YangYu,ChengJunsheng.ARollingBearingFaultDiagnosisMethodBasedonLCDDenoisingandVPMCD[J].ChinaMechanicalEngineering, 2013, 24(24):3338-3344.

[7]NikolaouNG,AntoniadisIA.ApplicationofMorphologicalOperatorsasEnvelopeExtractorsforImpulsive-typePerodicSignals[J].MechanicalSystemsandSignalProcessing, 2003,17:1147-1162.

[8]DongY,LiaoM,ZhangX,etal.FaultsDiagnosisofRollingElementBearingsBasedonModifiedMorphologicalMethod[J].MechanicalSystemsandSignalProcessing, 2011,25:1276-1286.

[9]WigginsRA.MinimumEntropyDeconvolution[J].Geoexploration, 1978,16:21-35.

[10]EndoH,RandallR.EnhancementofAutoregressiveModelBasedGearToothFaultDetectionTechniquebytheUseofMinimumEntropyDeconvolutionFilter[J].MechanicalSystemsandSignalProcessing, 2007,21:906-919.

(編輯王艷麗)

Applications of Mathematical Morphology Method to Fault Feature Extraction of Rolling Bearings Based on Minimum Entropy Deconvolution

Gong Tingkai1,2Yuan Xiaohui1Wang Xiyang2

1.Huazhong University of Science and Technology, Wuhan,430074 2.Nanchang Hangkong University, Nanchang, 330063

Aiming at the extractions of fault features of rolling bearings under the strong noise background, a novel method, called mathematical morphology based on minimum entropy deconvolution (MMBMED) was proposed herein. In this method, MED was first introduced to enhance the impact behaviors of the signals. And the fault characteristics of defective bearings were extracted by MM method, where the DIF filter with catching the bidirectional pulses was adopted as morphological operator, and kurtosis was employed as the criterion to the length selection of structural elements. The effectiveness of the new method was validated by both of simulation signals and vibration signals of rolling bearings with the outer and inner race faults. Finally, the results disclose that MED is effective to increase signal kurtosis and to strengthen the impulsive characteristics.

minimum entropy deconvolution(MED); mathematical morphology(MM); fault diagnosis; rolling bearing

2015-07-13

國家自然科學(xué)基金資助項(xiàng)目(51379080，AA201406168)

TH 165；TN911

10.3969/j.issn.1004-132X.2016.18.009

龔?fù)?，男?1977年生。華中科技大學(xué)數(shù)字化水電工程學(xué)院博士，南昌航空大學(xué)飛行器工程學(xué)院講師。主要研究方向?yàn)樾D(zhuǎn)機(jī)械故障診斷及信號(hào)處理。袁曉輝，男， 1971年生。華中科技大學(xué)數(shù)字化水電工程學(xué)院教授、博士研究生導(dǎo)師。王細(xì)洋，男，1967年生。南昌航空大學(xué)飛行器工程學(xué)院教授。