滚动轴承故障诊断(附MATLAB程序)

第二组实验

轴承故障数据：

Test2.mat 数据打开后应采用X105_DE_time作为分析数据，其他可作为参考，转速1797rpm

轴承型号：6205-2RS JEM SKF, 深沟球轴承

采样频率：12k Hz

1、确定轴承各项参数并计算各部件的故障特征频率

通过以上原始数据可知次轴承的参数为：

轴承转速r=1797r/min；滚珠个数n=9；滚动体直径d=7.938mm；

轴承节径D=39mm；：滚动体接触角α=0

由以上数据计算滚动轴承不同部件故障的特征频率为：

外圈故障频率f1=r/60 * 1/2 * n(1-d/D *cosα)=107.34Hz

内圈故障频率f2=r/60 * 1/2 * n(1+d/D *cosα)=162.21Hz

滚动体故障频率f3=r/60*1/2*D/d*[1-(d/D)^2* cos^2(α)]=70.53Hz 保持架外圈故障频率f4=r/60 * 1/2 * (1-d/D *cosα)=11.92Hz

2.对轴承故障数据进行时域波形分析

将轴承数据Test2.mat导入MATLAB中直接做FFT分析得到时域图如下：

并求得时域信号的各项特征：

（1）有效值：0.2909；（2）峰值：1.5256；（3）峰值因子：5.2441；（4）峭度：5.2793；（5）脉冲因子：7.2884；（6）裕度因子：9.1083：

3.包络谱分析

对信号做EMD 模态分解，分解得到的每一个IMF 信号分别和原信号做相关分析，找出相关系数较大的IMF 分量并对此IMF 分量进行Hilbert 变换。

s i g n a l

Empirical Mode Decomposition

i m f 1

i m f 2

i m f 3

i m f 4

i m f 5

i m f 6

i m f 7

i m f 8

r e s .

由图中可以看出经过EMD 分解后得到的９个IMF 分量和一个残余量。IMF 分量分别和原信号做相关分析后得出相关系数如下：

由上表得：IMF1的相关系数明显最大，所以选用IMF1做Hilbert 包络谱分析。所得Hilbert 包络谱图如下：

对包络谱图中幅值较大区域局部放大得到下图

由以上包络图的局部放大图中可以看出包络图中前三个峰值最大也最明显，三个峰值频率由小到大排列分别为58.59Hz、105.5Hz、164.1Hz。把这三个频率数值和前文计算所得的理论值进行比较可知：频率值最大为164.1Hz和内圈的故障理论计算特征频率f2=162.21Hz 相近，说明此轴承的故障发生在轴承的内圈。

clc

程序1：原始信号时域分析及小波去噪处理

clear all

z=importdata('C:\Users\wangkun\Desktop\轴承诊断\test2.mat');

x1=z.X105_DE_time(1:4096);

clear z;

N=4096;

fs=12000;

n=0:N-1;

t=n/fs;

f=n*fs/N;

figure(1);

plot(t,x1);

xlabel('t');

ylabel('幅值');

title('原信号时域图')

%小波去噪

[thr,sorh,keepapp]=ddencmp('den','wv',x1);

xd=wdencmp('gbl',x1,'db3',2,thr,sorh,keepapp);

figure(2);

plot(t,xd);

xlabel('t');

ylabel('幅值');

title('小波去噪后时域图')

程序2：EMD分解及Hilbert包络

clc

clear all

z=load('C:\Users\wangkun\Desktop\轴承诊断\test2.mat');

x=z.X105_DE_time(1:1024);

N=1024;

fs=12000;

n=0:N-1;

f=n*fs/N;

lag=N;

n=0:N-1;

t=n/fs;

imf=emd(x);

[m,n]=size(imf); %imf为一m*n阶矩阵，m是imf分量，n为数据点

emd_visu(x,1:length(x),imf,m); %实信号的信号重构及emd结果显示函数for i=1:m

a(i)=kurtosis(imf(i,:));%峭度

b(i)=mean(imf(i,:)); %均值;

c(i)=var(imf(i,:)); %方差;

d(i)=std(imf(i,:)); %均方值

e(i)=std(imf(i,:)).^0.5; %均方根值

f(i)=skewness(imf(i,:)); %计算偏度

end

[k,c]=max(a); %k为峭度最大值，c为最大元素在数组中的位置

[r,lags]=xcorr(x,lag,'unbiased'); %计算序列的自相关函数

for i=1:m

[R,lags]=xcorr(imf(i,:),lag,'unbiased'); %计算序列的自相关函数

a=corrcoef(R(1:N/2),r(1:N/2)); %相关系数矩阵【对称】，主对角元素为1 xg(i)=abs(a(1,2)); %相关系数

end

[R,C]=max(xg); %R为最大值，C为最大元素在数组中的位置

figure(4);

y = hilbert(imf(C,:));

a = abs(y);%包络

b=fft(a);

mag1=abs(b);

mag=mag1*2/N;

f1=(0:N-1)*fs/N;

plot(f1(1:N/2),mag(1:N/2));

%set(gca,'xlim',[0,.400]);

title('包络');

xlabel('频率');

ylabel('幅值');