李　民（上海久創(chuàng)建設(shè)管理有限公司　上海）

基于小波分析的電機(jī)軸承診斷Matlab程序設(shè)計(jì)

李民
（上海久創(chuàng)建設(shè)管理有限公司上海）

采用Matlab編程快速地在計(jì)算機(jī)上實(shí)現(xiàn)基于小波分析的電機(jī)滾動(dòng)軸承故障診斷。利用Daubechies小波對(duì)外圈或內(nèi)圈故障軸承的振動(dòng)信號(hào)進(jìn)行小波分解?；谛〔ǚ纸庀禂?shù)對(duì)含有故障特征頻率的第一層細(xì)節(jié)信號(hào)進(jìn)行小波重構(gòu)并提取其Hilbert包絡(luò)譜，從中檢測(cè)出故障特征頻率，據(jù)此判斷故障類型。設(shè)計(jì)了該算法的Matlab程序?；诿绹?guó)凱斯西儲(chǔ)大學(xué)提供的電機(jī)軸承振動(dòng)加速度信號(hào)采用上述算法進(jìn)行實(shí)驗(yàn)的結(jié)果表明，該方法可準(zhǔn)確地分辨出軸承故障的類型，非常適用于軸承的早期故障診斷。

軸承故障診斷小波分析Hilbert包絡(luò)譜Matlab

滾動(dòng)軸承是支承旋轉(zhuǎn)軸的關(guān)鍵部件，容易發(fā)生損壞，是電動(dòng)機(jī)中常見(jiàn)故障之一［1］。軸承故障診斷方法很多，其中振動(dòng)分析法是一種較實(shí)用的方法，其優(yōu)點(diǎn)在于振動(dòng)發(fā)自軸承本身，無(wú)需另加信號(hào)源；測(cè)試方便；可發(fā)現(xiàn)軸承的早期故障。當(dāng)軸承表面出現(xiàn)局部損傷時(shí)，在受載運(yùn)轉(zhuǎn)過(guò)程中軸承的其他零件會(huì)周期地撞擊損傷點(diǎn)，產(chǎn)生的沖擊力激勵(lì)軸承座及其支撐結(jié)構(gòu)，形成一系列由沖擊激勵(lì)產(chǎn)生的減幅振蕩，減幅振蕩發(fā)生的頻率即故障特征頻率［2］，根據(jù)該頻率可判斷發(fā)生故障的部位。然而故障信號(hào)在軸和軸上多種零部件振動(dòng)的干擾下往往被淹沒(méi)，無(wú)法直接利用頻譜分析檢測(cè)其特征頻率。小波變換是一種時(shí)頻分析方法，可進(jìn)行多分辨率分析，即，將信號(hào)分解成若干層次的細(xì)節(jié)信號(hào)及概貌信號(hào)。對(duì)軸承振動(dòng)信號(hào)進(jìn)行小波變換，提取其中具有故障特征的細(xì)節(jié)信號(hào)進(jìn)行重構(gòu)；對(duì)重構(gòu)信號(hào)做Hilbert包絡(luò)譜分析，從中檢測(cè)出軸承的故障特征頻率，據(jù)此判斷故障類型。

利用Matlab軟件編程快速地實(shí)現(xiàn)了基于小波分析的電機(jī)滾動(dòng)軸承外圈和內(nèi)圈故障診斷，并采用美國(guó)凱斯西儲(chǔ)大學(xué)提供的電機(jī)軸承振動(dòng)加速度信號(hào)進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果驗(yàn)證了該方法的有效性。

一、基于小波分析的電機(jī)軸承故障診斷

1.電機(jī)滾動(dòng)軸承故障的特征頻率

滾動(dòng)軸承由外圈、內(nèi)圈、滾動(dòng)體和保持架組成，工作時(shí)外圈與軸承座或機(jī)殼相聯(lián)接，固定或相對(duì)固定，內(nèi)圈與機(jī)械傳動(dòng)軸相聯(lián)接，隨軸一起轉(zhuǎn)動(dòng)。當(dāng)滾動(dòng)軸承表面發(fā)生損傷故障，如內(nèi)圈、滾動(dòng)體或外圈出現(xiàn)點(diǎn)蝕、裂紋或剝落等，根據(jù)不同的損傷部位，按以下公式分別計(jì)算軸承故障的特征頻率［3］（以外圈和內(nèi)圈剝落一點(diǎn)為例），見(jiàn)式（1）和式（2）。

2.基于小波分析的軸承故障診斷算法

一般采用加速度傳感器在軸承座上檢測(cè)滾動(dòng)軸承的振動(dòng)信號(hào)。若軸承表面出現(xiàn)局部損傷，在受載運(yùn)轉(zhuǎn)時(shí)軸承其他零件會(huì)周期地撞擊損傷點(diǎn)產(chǎn)生低頻的沖擊信號(hào)，其頻率即故障特征頻率，但檢測(cè)該頻率主要會(huì)遇到2個(gè)問(wèn)題：（a）沖擊信號(hào)的寬頻帶性質(zhì)會(huì)激起軸承結(jié)構(gòu)及傳感器本身在各自固有頻率上發(fā)生諧振，故軸承振動(dòng)信號(hào)中還含有故障特征頻率的高次諧波分量。（b）由于軸向間隙的存在，沖擊信號(hào)還要對(duì)軸承的高頻固有振動(dòng)信號(hào)進(jìn)行調(diào)制［4］。導(dǎo)致故障信號(hào)被其他振動(dòng)所干擾而無(wú)法直接通過(guò)頻譜分析檢測(cè)出故障特征頻率。本文分別采用小波分析與Hilbert包絡(luò)譜分析解決上述兩個(gè)問(wèn)題。

（1）小波分析提取含故障特征頻率的細(xì)節(jié)信號(hào)。小波是一種均值為零，很快衰減的瞬時(shí)振蕩函數(shù)，小波分析是一種時(shí)頻分析方法，它利用一系列伸縮和平移的小波函數(shù)對(duì)信號(hào)進(jìn)行展開(kāi)，該過(guò)程等效于用一系列不同頻帶的高通和低通濾波器將信號(hào)分解成若干層次的高頻細(xì)節(jié)信號(hào)及低頻概貌信號(hào)（即信號(hào)的主體輪廓），可對(duì)信號(hào)進(jìn)行多分辨率分解，被譽(yù)為“數(shù)學(xué)顯微鏡”。

小波分析算法的步驟包括分解與重構(gòu)，為在計(jì)算機(jī)上實(shí)現(xiàn)小波分析，根據(jù)二進(jìn)離散小波變換的快速算法Mallat算法進(jìn)行計(jì)算［5］，小波變換公式為式（3）。

式中cj，k——第j級(jí)小波分解所得低頻系數(shù)，設(shè)c0，k為原信號(hào)x（k）

n——小波分解的級(jí)數(shù)

dj，k——第j級(jí)小波分解所得高頻系數(shù)

h0（k）——離散尺度序列，是一低通濾波器

h1（k）——離散小波序列，是一高通濾波器

不同類型的小波，如Daubechies小波、Haar小波、墨西哥草帽小波等，濾波系數(shù)h0（k）與h1（k）均不相同。序列｛dn，k，dn-1，k，…，d1，k，c1，k｝是x（k）的二進(jìn)離散小波變換。利用小波分解系數(shù)重構(gòu)原信號(hào)的公式為式（4）。

根據(jù)公式（3）（4）對(duì)軸承振動(dòng)信號(hào)進(jìn)行小波分解與重構(gòu)可獲得其各層概貌信號(hào)Cj，k及細(xì)節(jié)信號(hào)Dj，k，其中幅度最大的細(xì)節(jié)信號(hào)中包含軸承故障的特征頻率。

（2）Hilbert變換包絡(luò)譜檢測(cè)軸承的故障特征頻率。含有軸承故障特征頻率的細(xì)節(jié)信號(hào)是一種調(diào)幅信號(hào)，它是故障信號(hào)對(duì)軸承的高頻固有振動(dòng)進(jìn)行幅度調(diào)制形成的，設(shè)其為式（5）。

式中A（t）——故障信號(hào)

fm——軸承固有振動(dòng)頻率

故f（t）的Hilbert變換可看成是f（t）通過(guò)一個(gè)幅度為1的全通濾波器輸出，其正頻率成分做-90°相移，負(fù)頻率成分做+90°相移。則調(diào)幅信號(hào)f（t）的Hilbert變換為式（8）。

設(shè)f（t）的解析信號(hào)為：s（t）=f（t）+jf?（t），則有式（9）。

因此可利用Hilbert變換提取f（t）的包絡(luò)，即故障信號(hào)A（t），再用傅里葉變換對(duì)其進(jìn)行功率譜分析，功率譜中幅度最大處的頻率即故障特征頻率。

二、軸承故障診斷的Matlab程序設(shè)計(jì)

Matlab（Matrix Laboratory矩陣實(shí)驗(yàn)室）是一款由美國(guó)MathWorks公司出品的數(shù)學(xué)軟件，包含許多專業(yè)領(lǐng)域的算法工具箱，編程十分方便。軸承故障診斷的實(shí)驗(yàn)數(shù)據(jù)來(lái)源于美國(guó)凱斯西儲(chǔ)大學(xué)（Case Western Reserve University）的軸承試驗(yàn)臺(tái)測(cè)試數(shù)據(jù)［6］，試驗(yàn)軸承支承電機(jī)主軸，試驗(yàn)軸承為6205-2RS JEM SKF，在其內(nèi)圈、外圈和滾動(dòng)體上應(yīng)用電火花技術(shù)分別加工了直徑0.18778 mm，深0.2794 mm的凹坑，模擬單點(diǎn)損傷。在電機(jī)主軸驅(qū)動(dòng)端軸承座上對(duì)應(yīng)軸承的正下方（即軸承負(fù)荷區(qū)）設(shè)置加速度傳感器測(cè)試軸承振動(dòng)信號(hào)，采樣頻率fs=12 000 Hz，試驗(yàn)中電機(jī)負(fù)荷為0，轉(zhuǎn)速n=1797 r/min，轉(zhuǎn)動(dòng)頻率f=n/60=29.95 Hz，滾動(dòng)體個(gè)數(shù)Z=9，接觸角α=0°，滾動(dòng)體直徑d=0.3126 mm，軸承的節(jié)徑E=1.537 mm，根據(jù)公式（1）和（2）計(jì)算得：外圈和內(nèi)圈故障特征頻率的理論值分別為f0=107.36 Hz，f1=162.19 Hz。設(shè)振動(dòng)信號(hào)采樣點(diǎn)數(shù)N= 12 000，由于使用快速傅里葉變換對(duì)信號(hào)做功率譜，其點(diǎn)數(shù)nfft應(yīng)＞N，且是2的指數(shù)次方，故取nfft=214=16 384；采用Daubechies小波“db2”對(duì)軸承振動(dòng)信號(hào)進(jìn)行3級(jí)小波分解與重構(gòu)，濾波系數(shù)為h0（k）=｛-0.1294，0.2241，0.8365，0.4830｝，h1（k）=｛-0.4830，0.8365，-0.2241，-0.1294｝，對(duì)第1層細(xì)節(jié)信號(hào)d1進(jìn)行Hilbert包絡(luò)譜分析。軸承故障診斷matlab程序及注釋（百分號(hào)后）如下：

clear%清除內(nèi)存

load 156.mat%調(diào)入故障軸承振動(dòng)信號(hào)的數(shù)據(jù)文件，存入sig

sig=X156_DE_time；

fs=12000；N=12000；Ts=1/fs；sig=sig（1：N）；%設(shè)置取樣頻率fs，取樣點(diǎn)數(shù)N

t=0：Ts：（N-1）*Ts；%時(shí)間軸t

sig=（sig-mean（sig））/std（sig，1）；%對(duì)sig進(jìn)行歸一化

subplot（211）；plot（t，sig）；%繪制sig波形

xlabel（＇時(shí)間t/s＇）；ylabel（＇振動(dòng)加速度/V＇）；

nfft=16384；

S=psd（sig，nfft）；%對(duì)sig做功率譜

subplot（212）；plot（（0：nfft/2-1）/nfft*fs，S（1：nfft/2））；%繪制功率譜

xlabel（＇頻率f/Hz＇）；ylabel（＇功率譜P/W＇）

［c，l］=wavedec（sig，4，＇db2＇）；%利用db2對(duì)sig進(jìn)行3級(jí)小波分解

c3=wrcoef（＇a＇，c，l，＇db2＇，3）；d3=wrcoef（＇d＇，c，l，＇db2＇，3）；%重構(gòu)第1-3層細(xì)節(jié)d1～d3和第3層概貌c3

d2=wrcoef（＇d＇，c，l，＇db2＇，2）；d1=wrcoef（＇d＇，c，l，＇db2＇，1）；figure；

subplot（414）；plot（t，c3）；ylabel（＇c3＇）；%繪制c3

subplot（413）；plot（t，d3）；ylabel（'d3'）；%繪制d3

subplot（412）；plot（t，d2）；ylabel（＇d2＇）；%繪制d2

subplot（411）；plot（t，d1）；ylabel（＇d1＇）；%繪制d1

y=hilbert（d1）；%對(duì)d1進(jìn)行Hilbert變換，得y

ydata=abs（y）；%ydata=|y|

ydata=ydata-mean（ydata）；%對(duì)ydata去均值（目的是去除幅度較大的直流分量）

P=psd（ydata，nfft）；%ydata的功率譜為P

figure；

plot（（0：nfft/2-1）/nfft*fs，P（1：nfft/2））；xlabel（＇頻率f/Hz＇）；%繪出d1的Hilbert包絡(luò)譜

P=P（1：nfft/2）；［M，f1］=max（P）；f1=f1*fs/nfft-1%故障頻率f1為包絡(luò)譜中幅度最大處的頻率

三、實(shí)驗(yàn)結(jié)果分析

首先對(duì)軸承振動(dòng)信號(hào)156.mat運(yùn)行上述Matlab程序，結(jié)果如圖1～圖3所示。圖1為該信號(hào)的歸一化時(shí)域波形及功率譜，顯然從功率譜中難以檢測(cè)故障頻率，圖2為對(duì)振動(dòng)信號(hào)做3級(jí)小波分解與重構(gòu)所得第1～3層細(xì)節(jié)信號(hào)d1～d3和第3層概貌信號(hào)c3，對(duì)整體幅度較大的細(xì)節(jié)信號(hào)d1做Hlibert包絡(luò)譜（見(jiàn)圖3（a）），其幅度最大處的頻率f1=107.3984 Hz，因此軸承故障類型為外圈故障；再用同樣方法對(duì)振動(dòng)信號(hào)105.mat進(jìn)行處理，得d1的包絡(luò)譜如圖3（b）所示，f1=161.5977 Hz，該軸承故障類型為內(nèi)圈故障。與理論值相比，利用db2小波檢測(cè)外圈和內(nèi)圈故障特征頻率的絕對(duì)誤差分別為：e0=0.0384 Hz和e1=0.5923 Hz，誤差較小，不影響對(duì)故障類型的識(shí)別；Matlab編程語(yǔ)句短小精悍，可快速地在計(jì)算機(jī)上實(shí)現(xiàn)故障診斷算法。

圖1　軸承振動(dòng)信號(hào)的時(shí)域波形及功率譜

圖2　軸承振動(dòng)信號(hào)的db2小波分析圖

圖3　第1層細(xì)節(jié)信號(hào)d1的包絡(luò)譜圖

四、結(jié)論

基于Matlab軟件編程，通過(guò)對(duì)電機(jī)滾動(dòng)軸承的振動(dòng)加速度信號(hào)進(jìn)行小波分解與重構(gòu)，提取其第1層細(xì)節(jié)信號(hào)做Hilbert包絡(luò)譜分析，快速地實(shí)現(xiàn)了對(duì)軸承外圈和內(nèi)圈（剝落1點(diǎn)）的故障特征頻率檢測(cè)，從而準(zhǔn)確地判別出故障類型。由于在表面損傷故障的初期，軸承元件的損傷點(diǎn)往往只有1點(diǎn)，故小波分析是一種有效的電機(jī)軸承早期故障診斷方法。

1 褚福磊，彭志科，馮志鵬，李志農(nóng).機(jī)械故障診斷中的現(xiàn)代信號(hào)處理方法［M］.北京：科學(xué)出版社.2009：159-160

2 葛哲學(xué)，陳仲生.MATLAB時(shí)頻分析技術(shù)及其應(yīng)用［M］.北京：人民郵電出版社.2006：164

3 沈玉娣，劉雄.機(jī)械故障診斷FORTRAN源程序匯編［M］.西安：西安交通大學(xué)出版社1990：170-171

4 劉圓，張永建，李元寶.改進(jìn)小波能量分析法在軸承故障診斷中的應(yīng)用［J］.煤礦機(jī)械.2012，33（5）：254-255

5 戴虹，錢(qián)晉武，張震，章亞男，沈林勇.GRNN在肌電預(yù)測(cè)踝關(guān)節(jié)運(yùn)動(dòng)中的應(yīng)用［J］.儀器儀表學(xué)報(bào).2013，34（4）：845-852

6 邱阿端.電機(jī)滾動(dòng)軸承故障診斷的新方法.中小型電機(jī)［J］.1996，23 （5）：47-50

7 美國(guó)凱斯西儲(chǔ)大學(xué)電機(jī)滾動(dòng)軸承振動(dòng)信號(hào)數(shù)據(jù)集［EB/OL］.http：// csegroups.case.edu/bearing data center/pages/12k-drive-end-bearingfault-data

〔編輯利文〕

TH113.1

B