張深逢　程曉萍　陳士釗　葉崗　宋云峰

（1.寧波舜宇智能科技有限公司　2.華北水利水電大學(xué)）

EMD故障診斷與激光測(cè)振技術(shù)的研究與應(yīng)用*

張深逢1程曉萍2陳士釗1葉崗1宋云峰1

（1.寧波舜宇智能科技有限公司2.華北水利水電大學(xué)）

介紹了經(jīng)驗(yàn)?zāi)B(tài)分解故障診斷方法，該方法統(tǒng)一了瞬時(shí)頻率的概念，產(chǎn)生時(shí)頻域分析方法—本征模態(tài)函數(shù)，可以突出局部數(shù)據(jù)特征，提取更準(zhǔn)確、更有效的原始信號(hào)特征信息，并經(jīng)過(guò)分解，提煉出有效時(shí)域信號(hào)，對(duì)其進(jìn)行Hilbert-Huang變換，實(shí)現(xiàn)信號(hào)在頻域中的再分析；提出的激光多普勒測(cè)振技術(shù)，具有抗干擾、高分辨率、高精度、非接觸式的振動(dòng)優(yōu)點(diǎn)。通過(guò)激光測(cè)振儀采集數(shù)控機(jī)床齒輪振動(dòng)信號(hào)，并借助經(jīng)驗(yàn)?zāi)B(tài)分解方法，診斷出軸承與軸瓦之間存在著頻率為33.3 Hz的周期性摩擦現(xiàn)象，從而證明了EMD能從大量的非線性、非平穩(wěn)信號(hào)中提取振動(dòng)特征信息和相關(guān)的模態(tài)參數(shù)，是一種非線性、非平穩(wěn)等時(shí)變信號(hào)處理方法。

經(jīng)驗(yàn)?zāi)B(tài)分解；時(shí)頻域聯(lián)合分析；故障診斷；激光多普勒

0　引言

在機(jī)械運(yùn)行狀態(tài)監(jiān)測(cè)、振動(dòng)分析和故障診斷過(guò)程中，存在大量的非線性、非平穩(wěn)信號(hào)。傅里葉變換（fast fourier transform，F(xiàn)FT）、時(shí)頻域分析（wigner-ville distrbution，WVD）和小波分析等信號(hào)分析方法，雖然給出時(shí)變性的描述，但對(duì)于出現(xiàn)瞬時(shí)頻率的局域性信號(hào)，則無(wú)法滿足功能需求。

美國(guó)宇航局Norden E. Huang于1998年提出經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，EMD）方法，它是一種時(shí)頻域分析方法，是指將一個(gè)復(fù)雜時(shí)間序列信號(hào)分解為有限個(gè)不同時(shí)間尺度的本征模態(tài)函數(shù)（intrinsic mode function，IMF）之和，每個(gè)本征模態(tài)函數(shù)所包含的頻率成分都與信號(hào)的分析頻率有關(guān)，而且隨著復(fù)雜時(shí)間序列信號(hào)的變化而變化。經(jīng)驗(yàn)?zāi)B(tài)分解方法的提出，統(tǒng)一了瞬時(shí)頻率概念，并且產(chǎn)生以時(shí)域基本信號(hào)為基礎(chǔ)的新時(shí)頻域分析方法—基本模式分量的概念，突出了數(shù)據(jù)局部特征，可以提取更準(zhǔn)確、更有效的原信號(hào)的特征信息。

目前，經(jīng)驗(yàn)?zāi)B(tài)分解應(yīng)用廣泛，如Huang N E 等人將EMD應(yīng)用于水波研究［1］、合成孔徑雷達(dá)圖像濾波，檢測(cè)沖擊波瞬時(shí)頻率的量值；Yue Huanyin等人將EMD應(yīng)用于地理學(xué)，設(shè)計(jì)了逐次逼近型模數(shù)轉(zhuǎn)換器的濾波器，提高了濾波效果［2］；Yu dejie等人將EMD運(yùn)用于滾子軸承的故障診斷中，檢測(cè)滾子軸承受力不均導(dǎo)致軸承和軸瓦周期性摩擦的沖擊故障［3-4］；Liu B等人將EMD應(yīng)用于齒輪箱的故障診斷中，檢測(cè)出齒輪嚙合出現(xiàn)周期性噪聲的原因是其中一個(gè)齒輪出現(xiàn)輕微裂紋［5］。

本文通過(guò)激光測(cè)振儀采集數(shù)控機(jī)床齒輪振動(dòng)信號(hào)，借助經(jīng)驗(yàn)?zāi)B(tài)分解方法，診斷出軸承與軸瓦之間存在著頻率為33.3 Hz的周期性摩擦現(xiàn)象，驗(yàn)證了EMD分解法在診斷摩擦故障方面的優(yōu)越性。

1　經(jīng)驗(yàn)?zāi)B(tài)分析理論

1.1Hilbert-Huang變換分析

經(jīng)驗(yàn)?zāi)B(tài)分解可以獲取原始信號(hào)的本征模態(tài)函數(shù)IMF，這些IMF是一系列相互關(guān)聯(lián)的時(shí)域信號(hào)，要從中獲取頻域中的信息，需要對(duì)本征模態(tài)進(jìn)行數(shù)據(jù)處理，Hilbert變換是一種最常用的方法。

Hilbert變換是信號(hào)分析中的重要工具，對(duì)任意給定信號(hào)x（ t），其Hilbert變換y（ t）定義為

其中，t是時(shí)間，τ是延遲時(shí)間。

由式（1）可知，Hilbert-Huang變換是x（ t）與1 πt的卷積，對(duì)應(yīng)頻域輸出如圖1所示。

圖1　Hilbert傳遞函數(shù)特性圖

由于任何信號(hào)都可以看作是正弦信號(hào)的疊加，因此，對(duì)于任意給定信號(hào)x（ t）可表示為

其中，a（ t）是信號(hào)x（ t）的調(diào)制信號(hào)；fs為載波高頻頻率；φ（t）為相位函數(shù)。

根據(jù)Hilbert-Huang變換原理，以x（ t）為實(shí)部，以x（ t）的Hilbert變換y（ t）為虛部，y（ t）亦可表達(dá)為y（ t）=a（ t）·sin［2π fst+φ（t）］，構(gòu)造解析函數(shù)z（ t）為

其中幅值函數(shù)

相位函數(shù)

a（ t）的絕對(duì)值|a（ t）|成為時(shí)間信號(hào)x（ t）的包絡(luò)線，用于信號(hào)的包絡(luò)線分析，也稱為被測(cè)信號(hào)的解調(diào)分析。這種分析方法成功運(yùn)用在齒輪、滾動(dòng)軸承的振動(dòng)分析中，并且得到了理想的解調(diào)分析結(jié)果。

對(duì)相位函數(shù)φ（t）求導(dǎo)，便可得到瞬時(shí)頻率

或

瞬時(shí)頻率是Hilbert變換的一個(gè)重要概念，瞬時(shí)頻率概念的提出衍生了本征模態(tài)函數(shù)。

1.2本征模態(tài)函數(shù)

經(jīng)過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解后的每個(gè)IMF分量，也可以構(gòu)造每個(gè)IMF的解析函數(shù)，進(jìn)而研究每個(gè)IMF的瞬時(shí)頻率，其算法

由此，可以分別得出每個(gè)IMF的幅值函數(shù)和相位函數(shù)。

幅值函數(shù)

瞬時(shí)相位函數(shù)

并由此可以求出瞬時(shí)頻率

或

這樣任意時(shí)間信號(hào)x（ t）可以表示為多個(gè)基本模式分量ci（ t）和余項(xiàng)rn（ t）之和。

其中rn（ t）是Hilbert變換的余項(xiàng)，也叫殘差函數(shù)，當(dāng)省略后，式（12）變?yōu)?/p>

式（13）稱為Hilbert幅值譜，另記作

式（14）精確描述了信號(hào)幅值在整個(gè)頻率段上的變化規(guī)律，因此，更確切地說(shuō)是完整信號(hào)能量時(shí)頻域圖譜。如果取式（14）中若干個(gè)IMF進(jìn)行局部分析，便可以得到Hilbert的局部圖譜。

式（15）精確地描述信號(hào)的幅值在所需要的頻率段上隨頻率和時(shí)間變化的規(guī)律。

由式（14）可以定義邊際譜h（ω）

可見(jiàn)，邊際譜h（ω）是時(shí)頻譜H（ω，t）對(duì)時(shí)間的積分，表達(dá)了頻率在整個(gè)頻域段上的能量貢獻(xiàn)程度，反映在概率上能量在整個(gè)頻率段上的積累。Hilbert變換實(shí)際是反映了在某個(gè)頻率上的幅值對(duì)整個(gè)頻率段上的權(quán)值貢獻(xiàn)，在邊際譜上某個(gè)頻率僅僅代表了該頻率振動(dòng)存在的概率大小，此振動(dòng)在Hilbert圖譜中發(fā)生的精確時(shí)間可用式（17）表示。

對(duì)比式（13）和式（17）可以看出，式（13）中ai（ t）、ωi（t）是時(shí)間函數(shù)的變量，式（17）中ai、ωi是常量。式（17）很好地描述了時(shí)間和頻率的定量關(guān)系，能夠?qū)崿F(xiàn)時(shí)變信號(hào)的完整分析。

對(duì)Hilbert時(shí)頻譜的平方再對(duì)頻率積分，可得到瞬時(shí)能量密度

由式（18）可以看出，對(duì)頻率積分后，只有自變量時(shí)間t，因此，瞬時(shí)能量密度IE（ t）是時(shí)間t的函數(shù)，反映能量隨時(shí)間t的波動(dòng)情況。以上基于經(jīng)驗(yàn)?zāi)B(tài)分解EMD的信號(hào)分析方法統(tǒng)稱為Hilbert-Huang變換（Hilbert-Huang Transformation，HHT）。

1.3分解篩選停止準(zhǔn)則

Hilbert-Huang變換用可變幅度（權(quán)值）和瞬時(shí)頻率分解信號(hào)，避免虛假諧波分量處理非平穩(wěn)、非線性信號(hào)存在的缺陷，一方面實(shí)現(xiàn)經(jīng)驗(yàn)?zāi)B(tài)分解將任意時(shí)間信號(hào)分解為本征模態(tài)函數(shù)加權(quán)之和；另一方面也實(shí)現(xiàn)用其求解每個(gè)本征模態(tài)函數(shù)的瞬時(shí)頻率，真正賦予具有局部時(shí)間尺度的本征模態(tài)函數(shù)和瞬時(shí)頻率實(shí)際的物理意義。Huang N E等研究發(fā)現(xiàn)，不是每一個(gè)時(shí)間信號(hào)x（ t）都具有本征模態(tài)函數(shù)，只有滿足以下2個(gè)基本條件，才能分解本征模態(tài)函數(shù)［6］。

1） 在整個(gè)時(shí)間段上，時(shí)間信號(hào)x（ t）極值點(diǎn)的數(shù)量Ne（含極大值和極小值）與過(guò)零點(diǎn)的數(shù)量Nz相等或差值為1，即

2） 對(duì)于任意給定的時(shí)間點(diǎn)ti∈（ta，tb），貫穿局部極大值的上包絡(luò)線fma（xt）和局部極小值的下包絡(luò)線fm（int）的均值為零，即

條件1類似于高斯平穩(wěn)過(guò)程的分布，條件2將全局限定變?yōu)榫植肯薅?，可避免波形不?duì)稱導(dǎo)致的瞬時(shí)頻率波動(dòng)，其實(shí)質(zhì)是用局部極大值和局部極小值的包絡(luò)線近似和代替，以使局部均值為零。鐘佑明等人在對(duì)本征模態(tài)函數(shù)的數(shù)學(xué)模型進(jìn)行詳細(xì)深入分析之后，論證了Hilbert建立的局部對(duì)稱性的必要性和用局部極大值極小值的包絡(luò)線近似和代替的合理性［7］。

然而，實(shí)際上很難滿足這2個(gè)條件，尤其是局部極大值、極小值點(diǎn)的均值為零的條件幾乎無(wú)法達(dá)到。而且，完全按照這2個(gè)條件分解本征模態(tài)函數(shù)需要多次篩選，計(jì)算量較大，也導(dǎo)致分解本征模態(tài)函數(shù)失去了其實(shí)際的物理意義，因此，必須確定篩選過(guò)程停止準(zhǔn)則。

篩選過(guò)程的停止準(zhǔn)則可用經(jīng)驗(yàn)法表示，通過(guò)限制2個(gè)連續(xù)準(zhǔn)本征模態(tài)函數(shù)處理結(jié)果的標(biāo)準(zhǔn)差Sd的大小來(lái)確定，如式（21）所示。

其中，i為第i個(gè)本征模態(tài)函數(shù)；T為信號(hào)的時(shí)間寬度；hi（k-1）（t）、hik（ t ）是篩選本征模態(tài)函數(shù)過(guò)程中2個(gè)連續(xù)的處理結(jié)果的時(shí)間信號(hào)；k為篩選的次數(shù)。

標(biāo)準(zhǔn)差Sd的值越小，分解得到的本征模態(tài)函數(shù)的線性和穩(wěn)定性越好，能夠分解的本征模態(tài)函數(shù)的數(shù)量也越多。經(jīng)驗(yàn)表明，當(dāng)Sd在0.2～0.3范圍時(shí)，不僅能保證得到的本征模態(tài)函數(shù)的線性性和穩(wěn)定性，而且能使所得本征模態(tài)函數(shù)具有實(shí)際的物理意義。

2　經(jīng)驗(yàn)?zāi)B(tài)分解算法過(guò)程

2.1算法分解流程

EMD貫穿了時(shí)域頻域分析，然而其有較大的局限性，如它要求信號(hào)必須滿足2個(gè)限定條件，才能分解出若干個(gè)IMF，但大多數(shù)信號(hào)并不能滿足這2個(gè)苛刻的條件。Huang N E等假設(shè)：任何信號(hào)都由一定內(nèi)在聯(lián)系的IMF組成，每個(gè)IMF，或非線性，各IMF之間相互關(guān)聯(lián)，形成組合信號(hào)。Huang N E等進(jìn)一步指出［8］，可用改進(jìn)的EMD算法，將所需本IMF先提取出來(lái)，再作進(jìn)一步分析，此過(guò)程也稱篩選過(guò)程，其實(shí)質(zhì)是基于數(shù)據(jù)的特征時(shí)間尺度來(lái)獲取各IMF。

基于IMF分解的定義，信號(hào)分解的最終目的是得到使瞬時(shí)頻率概念有意義的IMF，按照分解IMF的2個(gè)限定條件，其分解原理步驟如下［9］：

1） 把信號(hào)x（ t）作為待處理信號(hào)，找出該信號(hào)所有局部極大值點(diǎn)、極小值點(diǎn)；用3次樣條曲線連接局部極大、極小值點(diǎn)，得到信號(hào)x（ t）的上、下2條包絡(luò)線，則信號(hào)x（ t）的所有點(diǎn)均位于上、下包絡(luò)線之間；取上、下2條包絡(luò)線均值組成的有序序列為m（ t），如圖2所示。 “?！焙汀?”標(biāo)志分別代表信號(hào)x（ t）的極大、極小值點(diǎn)，上、下2條虛線表示極值點(diǎn)的上、下包絡(luò)線，點(diǎn)劃線表示上、下兩條包絡(luò)線的均值序列m（ t）。

圖2 EMD分解原理圖

2） 用信號(hào)x（ t）減去其上、下包絡(luò)線的均值m（ t），可以得到h1（ t）

核查h1（ t）是否滿足IMF的2個(gè)限定條件，如不能，則把h1（ t）當(dāng)待處理信號(hào)，重復(fù)以上操作，直到h1（ t）滿足IMF的條件要求，此時(shí)，記為

3） 提取信號(hào)x（ t）中的第1個(gè)c1（ t）后，從信號(hào)x（ t）中減去c1（ t），得到剩余時(shí)間r1（ t）。

4） 把r1（ t）當(dāng)作新的待處理信號(hào)“c1（ t）”，重復(fù)上述操作，依次可以得到n個(gè)IMF，分別記為c2（ t），c3（ t），···，cn（ t ），這個(gè)分解過(guò)程滿足停止準(zhǔn)則式（21）時(shí)停止，最后得到余項(xiàng)rn（ t）。這樣，將信號(hào)x（ t）分解為n個(gè)IMF和1個(gè)余項(xiàng)的和，即

歸納以上4個(gè)步驟，將Hilbert-Huang變換的算法以程序流程圖的形式繪出，如圖3所示。

圖3　Hilbert-Huang變換算法程序流程圖

然而，Hilbert-Huang變換亦有其缺點(diǎn)，基于EMD的Hilbert-Huang變換每次要迭代，信號(hào)包絡(luò)線都要用三次樣條插值前后各2個(gè)臨近點(diǎn)。但是，三次樣條插值要求信號(hào)臨近點(diǎn)是極值，但實(shí)際中信號(hào)臨近點(diǎn)不一定為極值，于是插值求包絡(luò)線時(shí)便把端點(diǎn)當(dāng)作極值點(diǎn)，從而引起差值失真，導(dǎo)致分解出的各個(gè)IMF在兩端點(diǎn)附近失去物理意義。當(dāng)只分解出1個(gè)IMF時(shí)，這種端點(diǎn)效應(yīng)影響較小，但當(dāng)分解出多IMF的復(fù)雜信號(hào)來(lái)說(shuō)，特別是需要作多次EMD分解，邊緣效應(yīng)會(huì)無(wú)限放大，嚴(yán)重淹沒(méi)了信號(hào)的特征信息，這種在分解過(guò)程中兩端點(diǎn)附近出現(xiàn)失真的現(xiàn)象就是“端點(diǎn)效應(yīng)”［10］。

2.2算法端點(diǎn)效應(yīng)

端點(diǎn)效應(yīng)的出現(xiàn)，嚴(yán)重限制了EMD的應(yīng)用，如對(duì)高頻信號(hào)來(lái)說(shuō)，失真程度不會(huì)很嚴(yán)重，分解的IMF誤差不太大；但對(duì)低頻信號(hào)來(lái)說(shuō)，三次樣條插值將不能決定靠近極值點(diǎn)以外的數(shù)據(jù)走向，端部的邊緣效應(yīng)會(huì)傳遞到信號(hào)的內(nèi)部，這樣分解的IMF被影響的程度和傳遞到信號(hào)內(nèi)部的長(zhǎng)度都受到了極大影響，更難以獲得有用的信息，甚至完全失去原始信號(hào)的絕大部分信息，或者獲得錯(cuò)誤很大的分析結(jié)果。因此，需要研究有效的端點(diǎn)效應(yīng)抑制方法［11-12］。

出現(xiàn)端點(diǎn)效應(yīng)的原因有：數(shù)據(jù)序列的長(zhǎng)度有限；三次樣條插值時(shí)需要用到前后各2個(gè)臨近點(diǎn)，而這2個(gè)端點(diǎn)不一定是極值點(diǎn)。因此，可以這2個(gè)原因?yàn)橹贮c(diǎn)分析研究抑制端點(diǎn)效應(yīng)的方法［13-15］：延長(zhǎng)信號(hào)序列的長(zhǎng)度或者在數(shù)據(jù)兩端增加極值點(diǎn)；采用其他的樣條插值函數(shù)；基于Hilbert-Huang變換理論。

通常采用以下幾種方法抑制端點(diǎn)效應(yīng)：

1） 直接對(duì)原始數(shù)據(jù)的端點(diǎn)為極值點(diǎn)進(jìn)行簡(jiǎn)單延拓；

2） 采用全局統(tǒng)計(jì)平均方法延拓極值；

3） 采用神經(jīng)網(wǎng)絡(luò)方法對(duì)數(shù)據(jù)延拓；

4） 采用局部統(tǒng)計(jì)平均方法延拓極值；

5） 在端點(diǎn)處數(shù)據(jù)“平衡位置”附加2條平行線段的方法延拓極值；

6） 采用多項(xiàng)式擬合的方法延拓極值點(diǎn)；

7） 基于AR模式的時(shí)間序列線性預(yù)測(cè)方法；

8） 鏡像延拓法延拓極值；

9） 具有更高預(yù)測(cè)精度的支持向量機(jī)（support vector machine，SVM）的雙邊延拓等。

以上9種抑制端點(diǎn)效應(yīng)的方法各自有其適用場(chǎng)合，也各有其優(yōu)勢(shì)側(cè)重點(diǎn)，現(xiàn)選取幾種較常用的方法總結(jié)歸納如下：

全局統(tǒng)計(jì)平均方法延拓極值是指分析待處理信號(hào)的極大值、極小值，定出兩端點(diǎn)附加極值點(diǎn)的位置和幅值，構(gòu)造新區(qū)間，并確保極大值、極小值的區(qū)間不小于當(dāng)前待分析信號(hào)的長(zhǎng)度。全局統(tǒng)計(jì)平均方法能在一定程度上抑制端點(diǎn)效應(yīng)，但是進(jìn)行全局統(tǒng)計(jì)時(shí)，如果信號(hào)中含有幾個(gè)幅值較大的極值時(shí)，會(huì)使兩端點(diǎn)附近的分析結(jié)果出現(xiàn)誤差和畸變，這種情況可采用局部統(tǒng)計(jì)方法延拓極值。局部統(tǒng)計(jì)方法延拓極值與此類似，不再贅述。全局統(tǒng)計(jì)平均方法延拓極值和局部統(tǒng)計(jì)方法延拓極值常常結(jié)合使用，以各展所長(zhǎng)，彌補(bǔ)缺陷。

采用神經(jīng)網(wǎng)絡(luò)方法對(duì)數(shù)據(jù)延拓，是基于神經(jīng)網(wǎng)絡(luò)對(duì)于集合復(fù)雜關(guān)系的強(qiáng)映射能力，其基本理論基礎(chǔ)是假定BP神經(jīng)網(wǎng)絡(luò)中隱含層單元可以根據(jù)需要自由設(shè)定，3層BP網(wǎng)絡(luò)可以實(shí)現(xiàn)任意精度的映射。其抑制端點(diǎn)效應(yīng)的具體步驟如下：

1） 以端點(diǎn)［x0，xn］為三次樣條的插值區(qū)間，選定兩端點(diǎn)作為起始點(diǎn)，向內(nèi)m個(gè)點(diǎn)為神經(jīng)網(wǎng)反向和前向外插訓(xùn)練樣本；

2） 建立反向和前向外插神經(jīng)網(wǎng)絡(luò)，該網(wǎng)絡(luò)結(jié)構(gòu)為帶有反饋的遞歸網(wǎng)絡(luò)結(jié)構(gòu)（RBP網(wǎng)絡(luò)），RBP網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示。

圖4　反向和前向外插值遞歸型RBP網(wǎng)絡(luò)結(jié)構(gòu)

3） 分別用反向和前向外插神經(jīng)網(wǎng)絡(luò)對(duì)本RBP網(wǎng)絡(luò)進(jìn)行訓(xùn)練，即用訓(xùn)練實(shí)現(xiàn)映射訓(xùn)練。同樣，其他映射反向訓(xùn)練過(guò)程類似前向訓(xùn)練。

用訓(xùn)練過(guò)RBP網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行反向和前向外插逐點(diǎn)計(jì)算，以前向訓(xùn)練為例，其過(guò)程用訓(xùn)練yn+1，然后將yn+1反饋到網(wǎng)絡(luò)輸入端，繼續(xù)用訓(xùn)練依次類推，預(yù)測(cè)步數(shù)在一定的范圍，以保證訓(xùn)練精度［16］。

經(jīng)過(guò)多年的分析研究，國(guó)內(nèi)眾多學(xué)者分別提出工程上比較實(shí)用的自適應(yīng)時(shí)變?yōu)V波法ATVFD［17-18］、極值域均值模式分解法EMMD［19］和計(jì)算精確度更高的極值域均值模式分解算法IEMMD［20］，這3種方法都是局部均值求解的實(shí)用化技術(shù)。下面闡述局部均值求解速度和精度都比較好的極值域均值模式分解算法IEMMD。

極值域均值模式分解算法IEMMD取消了相鄰極值點(diǎn)變化均勻的假設(shè)，求出所有局部極值點(diǎn)組成的極值點(diǎn)序列e（ ti），其中

其中，ti是第i個(gè)時(shí)間；ti+1是第i+1個(gè)時(shí)間。

由于極值點(diǎn)e（ ti）中的極大、極小值點(diǎn)呈間隔排列的，則局部極值序列mi（ tε1）與信號(hào)有唯一的交點(diǎn)，對(duì)應(yīng)時(shí)間為tε1；mi+1（tε2）與信號(hào)也有唯一交點(diǎn)，對(duì)應(yīng)時(shí)間tε2，如圖5所示。設(shè)mi（ tε）在待分析信號(hào)x（ t）中介于x（ tj）和x（ tj+1）之間，1≤j≤k-1。

圖5　信號(hào)、極值點(diǎn)與局部均值的關(guān)系圖

然后用均值序列mi（ tε1）和mi+1（tε2）以數(shù)學(xué)加權(quán)平均，求其在ti+1處的局部均值m（ ti+1），即

其中，k（ ti）和k（ ti+1）是通過(guò)梯形相似的幾何特性得到的加權(quán)系數(shù)

IEMMD既使用了信號(hào)局部數(shù)據(jù)，又用到中值定理求解mi（ tε）時(shí)對(duì)應(yīng)tε，這樣得到的局部均值較正確，瞬時(shí)頻率得到保證，分解出的IMF精度和時(shí)頻分辨率更高。

EMD分解方法的前提是所采集信號(hào)一定要準(zhǔn)確反映機(jī)械設(shè)備運(yùn)行狀態(tài)，因此，采集數(shù)據(jù)的質(zhì)量也是影響EMD分解方法的重要因素。本文采用激光高速測(cè)振的方法采集數(shù)據(jù)，該方法具有抗干擾、高分辨率、高精度、非接觸式的優(yōu)點(diǎn)。

3　EMD聯(lián)合激光測(cè)振技術(shù)的結(jié)果與驗(yàn)證

3.1激光測(cè)振技術(shù)

激光測(cè)振技術(shù)是一種非接觸的測(cè)量振動(dòng)技術(shù)，該測(cè)振技術(shù)具有抗干擾、高分辨率、高精度、非接觸的特性。激光多普勒測(cè)振儀可以應(yīng)用在許多接觸式測(cè)振方式無(wú)法測(cè)量的任務(wù)中，出色的頻率和相位響應(yīng)，可準(zhǔn)確地對(duì)各種物體的振動(dòng)、位移、速度及加速度等進(jìn)行測(cè)量。在滿足高精度、高速測(cè)量需求的同時(shí)，還可以彌補(bǔ)接觸式測(cè)量方法無(wú)法測(cè)量大幅度振動(dòng)的缺陷。以激光多普勒測(cè)振儀（Laser Doppler Vibrometer，LDV）為例，它主要由1臺(tái)高精度激光干涉儀和1臺(tái)信號(hào)處理器組成，高精度激光干涉儀內(nèi)的He-Ne激光器發(fā)出的偏振光（頻率為f0）由分光鏡分成2束：一路作為測(cè)量；一路用于參考。測(cè)量光通過(guò)聲光調(diào)制器具有一定頻移F，再被聚焦到被測(cè)物體表面，物體振動(dòng)引起頻移（f=2ν λ）。系統(tǒng)手機(jī)反射光并與參考光匯聚在傳感器上，2束光在傳感器表面形成干涉，干涉信號(hào)的頻率為F+f，攜帶了被測(cè)物體的振動(dòng)信息，信號(hào)處理器將頻移信號(hào)轉(zhuǎn)換為速度和位移信號(hào)［21］。激光測(cè)振系統(tǒng)原理圖如圖6所示。

圖6　激光測(cè)振系統(tǒng)原理圖

3.2EMD聯(lián)合激光測(cè)振技術(shù)驗(yàn)證分析

數(shù)控機(jī)床軸承、齒輪發(fā)生損傷時(shí)，在損失部位會(huì)產(chǎn)生沖擊脈沖激勵(lì)，出現(xiàn)振蕩衰減的脈沖響應(yīng)信號(hào)，該信號(hào)被軸承或齒輪的振動(dòng)頻率信號(hào)所調(diào)制，則對(duì)齒輪或軸承進(jìn)行故障診斷的理想方法是利用EMD的包絡(luò)技術(shù)，獲得損失部位對(duì)應(yīng)特征頻率，再利用EMD逐層抽取IMF函數(shù)，可以得到識(shí)別軸承或齒輪損失程度的結(jié)果。應(yīng)用EMD故障診斷方法，采用激光測(cè)振系統(tǒng)，對(duì)某型號(hào)數(shù)控機(jī)床的齒輪損失程度監(jiān)測(cè)診斷如下：采樣頻率為1000 Hz，采樣點(diǎn)數(shù)512，機(jī)床轉(zhuǎn)速2000 r/min，進(jìn)刀速度60 mm/min，單邊進(jìn)刀量2 mm，采用開(kāi)發(fā)的多通道數(shù)據(jù)采集診斷系統(tǒng)對(duì)機(jī)床某處齒輪進(jìn)行振動(dòng)數(shù)據(jù)采集，得到圖7所示振動(dòng)時(shí)域信號(hào)和圖8所示功率譜，圖7和圖8是在LabVIEW環(huán)境中編程實(shí)現(xiàn)的。

圖7　機(jī)床齒輪振動(dòng)信號(hào)的時(shí)域信號(hào)

圖8　機(jī)床齒輪振動(dòng)時(shí)域信號(hào)的頻譜圖

從圖7中可以看到，損傷的齒輪波形有些雜亂，但不太明顯，說(shuō)明齒輪損傷程度不大。但從圖8頻譜分析圖中可以看到，損傷齒輪存在工頻或噪聲成分，在33 Hz、66 Hz和158 Hz 3處存在著次譜峰、第3譜峰和高譜峰，由于機(jī)床轉(zhuǎn)速n=2000r/min，對(duì)應(yīng)頻率為f=n60=33.3 Hz 。由此看出，66 Hz譜峰是機(jī)床齒輪振動(dòng)二倍工頻，158 Hz的高譜峰與五倍工頻166 Hz相近，但不相等，說(shuō)明齒輪損傷影響著齒輪工頻，即振動(dòng)信號(hào)中存在頻率調(diào)制的現(xiàn)象，且調(diào)制頻率以一倍工頻33.3 Hz為基頻。這解釋了轉(zhuǎn)子周期性碰撞摩擦的原因：轉(zhuǎn)子每轉(zhuǎn)動(dòng)一周，摩擦一次，線速度減小一次，摩擦過(guò)后，線速度重新恢復(fù)為額定值，此即為工頻振動(dòng)分量的頻率調(diào)制現(xiàn)象。

運(yùn)用Matlab環(huán)境下EMD故障診斷算法繼續(xù)分析，得到Matlab下基于EMD分解的6個(gè)基本模式分量，對(duì)應(yīng)EMD公式中的c1（ t）、c2（ t）、c3（ t）、c4（ t ）、c5（ t）、c6（ t）以及最后殘差余項(xiàng)rn（ t），齒輪振動(dòng)信號(hào)EMD分解結(jié)果之殘差余項(xiàng)如圖9所示。由于數(shù)據(jù)量比較大，兩端點(diǎn)失真的現(xiàn)象從外部邊緣延伸至內(nèi)部數(shù)據(jù)的長(zhǎng)度相對(duì)較小，因此，選取信號(hào)內(nèi)部數(shù)據(jù)進(jìn)行EMD分析，端點(diǎn)抑制效應(yīng)的影響會(huì)很小，幾乎可以忽略不計(jì)。對(duì)圖7的振動(dòng)數(shù)據(jù)進(jìn)行EMD分析的前2個(gè)經(jīng)驗(yàn)?zāi)B(tài)分解IMF1和IMF2如圖10所示。

從圖10中可以看到，EMD分解結(jié)果的前2個(gè)IMF，即IMF1和IMF2的幅值較大，IMF1最大值可達(dá)4.8 μm，IMF2的最大值則為2.5μm，而且IMF1和IMF2所包含的頻率成分較多，適合對(duì)c1（ t）、c2（ t）作頻譜分析，以獲得更準(zhǔn)確的頻率特性信息。

而圖11中，EMD分解的后2個(gè)IMF，即IMF3和IMF4幅值為0.05μm，幅值較圖10小，縱坐標(biāo)最大幅值幾乎是圖10的1/100，而且從圖形看，IMF1和IMF2的波形與原始信號(hào)圖7相近，而圖11所示的結(jié)果低頻率區(qū)域（0～180 Hz）幾乎是水平的曲線。這也就是說(shuō)IMF1和IMF2包含的原始信號(hào)的特征信息較多，IMF3和IMF4包含的原始信息較少。

而圖12中，EMD分解的后2個(gè)IMF，即IMF5和IMF6縱坐標(biāo)最大值幅值為0.005μm，幅值更小，縱坐標(biāo)最大幅值幾乎是圖11的1/10，是圖10的1/1000，整條曲線在低頻率區(qū)域更趨于一條直線，包含的頻率成分更少，已經(jīng)不適合做頻率分析。因此，選取包含原始特征信息較多的IMF1和IMF2進(jìn)行分析。

圖9　齒輪振動(dòng)信號(hào)EMD分解結(jié)果之殘差余項(xiàng)

圖10　齒輪振動(dòng)信號(hào)EMD分解結(jié)果之IMF1和IMF2

圖11　齒輪振動(dòng)信號(hào)EMD分解結(jié)果之IMF3和IMF4

圖12　齒輪振動(dòng)信號(hào)EMD分解結(jié)果之IMF5和IMF6

從圖13所示的IMF1的瞬時(shí)頻率曲線的頻譜圖中可以看到，橫坐標(biāo)在頻率為33.3 Hz處譜峰最高，出現(xiàn)強(qiáng)烈振動(dòng)，在33.3 Hz的倍數(shù)頻出分別出現(xiàn)不同程度的譜峰，這與圖8所示的機(jī)床齒輪振動(dòng)時(shí)域信號(hào)的頻譜圖相吻合，驗(yàn)證了在轉(zhuǎn)子運(yùn)轉(zhuǎn)時(shí)，軸承與軸瓦之間存在著一定程度的摩擦現(xiàn)象，說(shuō)明了EMD分析方法的正確性，也使EMD分解法在診斷摩擦故障方面得到了突出表現(xiàn)。

圖13　齒輪振動(dòng)的IMF1的瞬時(shí)頻率曲線圖的頻譜圖

4　結(jié)論

實(shí)驗(yàn)結(jié)果表明，EMD能從大量的非線性、非平穩(wěn)信號(hào)中提取振動(dòng)特征信息和相關(guān)的模態(tài)參數(shù)，能夠定量地描述頻率和時(shí)間的關(guān)系，通過(guò)時(shí)頻域分析，把信號(hào)自適應(yīng)地分解到不同頻帶，分解出基本模式分量，再對(duì)基本模式分量進(jìn)行解調(diào)分析和包絡(luò)分析，實(shí)現(xiàn)對(duì)時(shí)變信號(hào)完整、準(zhǔn)確的分析，無(wú)需人為劃分頻帶，層層挖掘數(shù)據(jù)中的有效信息，提取出局部瞬時(shí)頻率和信號(hào)中的沖擊、脈沖、振蕩分量等非線性、非平穩(wěn)信息，在機(jī)械動(dòng)態(tài)分析和故障診斷中，尤其在軸承、齒輪的損傷診斷、調(diào)制解調(diào)振蕩信號(hào)以及診斷摩擦故障方面，有其獨(dú)特的優(yōu)勢(shì)，且效果明顯，是一種非常優(yōu)秀的非線性、非平穩(wěn)等時(shí)變信號(hào)處理方法。

［1］ Huang N E， Shen Zheng， Long S R. A new view of nonlinear water waves： The Hilbert spectrum［J］. Annu. Rev.Fluid Mech.，1999，31：417-457.

［2］ Yue Huanyin， Guo Huadong， Han Chunming. A SAR interferogram filter based on the empirical mode decomposition method［J］.Geosciences and Remote Sensing Symposium， 2001，5：2061-2063.

［3］ Yu Dejie， Cheng junsheng， Yang Yu. Application of EMD method and Hilbert spectrum to the fault diagnosis of roller bearings［J］.Mechanical System and Signal Processing 19（2005） ：259-270.

［4］ Huang N E， Shen Zheng， Long S R. The empirical mode decomposition and the Hilbert spectrum for non-stationary time series analysis ［J］. Proc.R.Soc.Lond，1988，454：903-995.

［5］ Liu B，Riemenschneider S， Xu Y. Gearbox fault diagnosis using empirical mode decomposition and Hilbert spectrum［J］.Mechanical System and Signal Processing 2006，20（3）：718-734.

［6］ Huang N E. A New method for nonlinear and non-stationary time series analysis and its application for civil infrastructurehealth monitoring. NASA Goddard Space Flight Center，Greenbelt， Maryland 20771，USA.

［7］ 鐘佑明，秦樹(shù)人，湯寶平.Hilber-Huang變換中的理論研究［J］.振動(dòng)與沖擊，2002，21（4）：13-17.

［8］ Huang N E. A new view of earthquake ground motion data：The Hilbert spectrum analysis［C］.Proc. Intl. workshop on annual commemoration of Chi-Chi Earthquake， 2000，Ⅱ：64-75.

［9］ Huang N E， Shen Zheng， Long S R. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［C］. Proceeding of the Royal Society of London ， 1998，A（454）：903-905.

［10］ Deng Yongjun， Wang Wei， Qian Chengchun， et al. Boundary processing technique in EMD method and Hilbert transform［J］.Chinese science Bulletin，2001，46（11）：257-263.

［11］ 譚冬梅，姚三，瞿偉廉.振動(dòng)模態(tài)的參數(shù)識(shí)別綜述［J］.華中科技大學(xué)學(xué)報(bào)：城市科學(xué)版，2002，19（3）：73-78.

［12］ 徐小力，許寶杰，韓秋實(shí).大型機(jī)械設(shè)備在線趨勢(shì)預(yù)測(cè)的預(yù)測(cè)模型及預(yù)測(cè)方法［C］.第十一屆全國(guó)設(shè)備監(jiān)測(cè)與診斷學(xué)術(shù)會(huì)議論文集.北京：中國(guó)宇航出版社，2003：80-84.

［13］ 胥永剛.機(jī)電設(shè)備監(jiān)測(cè)診斷時(shí)域新方法的應(yīng)用研究［D］.西安：西安交通大學(xué)，2003.

［14］ 王紅軍，付瑤.基于多項(xiàng)式擬合的EMD端點(diǎn)效應(yīng)處理方法研究［J］.機(jī)械設(shè)計(jì)與制造，2010（10）：197-199.

［15］ 許寶杰，張建民，徐小力，等.抑制EMD端點(diǎn)效應(yīng)方法的研究［J］.北京理工大學(xué)學(xué)報(bào)，2006，26（3）：196-200.

［16］ 許寶杰，徐小力，李建偉.用神經(jīng)網(wǎng)絡(luò)插值進(jìn)行EMD端點(diǎn)效應(yīng)抑制方法的研究［C］.第十二屆全國(guó)設(shè)備監(jiān)測(cè)與診斷學(xué)學(xué)術(shù)會(huì)議論文集.北京：機(jī)械工業(yè)出版社，2006：171-175.

［17］ 盧志茂，孫美玲，張春祥，等.基于極值域均值模式分解的語(yǔ)音增強(qiáng)方法［J］.系統(tǒng)工程與電子技術(shù)，2011，33（7）：1680-1684.

［18］ 崔寶珍.自適應(yīng)形態(tài)濾波與局域波分解理論及滾動(dòng)軸承故障診斷［D］.太原：中北大學(xué)，2013.

［19］ 余波.自適應(yīng)時(shí)頻方法及其在故障診斷方法中的應(yīng)用研究［D］.大連：大連理工大學(xué)，1998.

［20］ 蓋強(qiáng).局域波時(shí)頻分析方法的理論研究與應(yīng)用［D］.大連：大連理工大學(xué)，2001.

［21］ 葉聲華，楊學(xué)友，曲興華，等.激光及光電測(cè)試技術(shù)［C］.葉聲華科研團(tuán)隊(duì)論文集.北京：機(jī)械工業(yè)出版社，2013：14-21.

The Research and Application of the EMD Fault Diagnosis with Laser Vibration Technology

Zhang Shenfeng1Cheng Xiaoping2Chen Shizhao1Ye Gang1Song Yunfeng1
（1.Ningbo Shunyu Intelligent Measuring Instrument Co.， Ltd. 2. North China University of Water resources and Electric Power）

Experience mode decomposition fault diagnosis method is introduced， the method unifies the concept of instantaneous frequency and promethean combination of time domain and frequency domain， generates the time-frequency domain analysis method， the intrinsic mode function， can highlight local characteristics of data， extract more accurate， more effective original signal characteristic information， and through the decomposition， extract the time domain signal effectively， through the Hilbert-Huang transform， to realize signal's reanalysis in the frequency domain； Laser Doppler vibration technology， has advantages such as anti-jamming， high resolution， high accuracy， non-contact. The subject collects the gear vibration signals of nc machine tool through laser vibrometer， and using empirical mode decomposition method， diagnoses that there exists the periodic friction phenomenon at 33.3 Hz frequency between the bearing and bearing shell， which proves that the EMD can extract vibration characteristic information and the relevant modal parameters from a large number of nonlinear and non-stationary signals， and it is an excellent method in processing nonlinear， non-stationary and time-varying signals.

Empirical Mode Decomposition； Joint Time-Frequency Domain Analysis； Fault Diagnosis； Laser Doppler

張深逢，男，1988年生，碩士，工程師，主要研究方向：機(jī)械振動(dòng)狀態(tài)監(jiān)測(cè)及故障診斷技術(shù)、旋轉(zhuǎn)機(jī)械故障診斷技術(shù)以及信號(hào)分析處理等。E-mail： changgongzsf@163.com

國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專項(xiàng)（2013YQ470765）； 2015年 寧波市科技創(chuàng)新團(tuán)隊(duì)（2013B82005）