爱分享

电脑基础 · 2023年4月19日

基于小波时频图和2D-CNN的滚动轴承故障检测

目录

一、研究思路

 1、基于小波时频图和CNN的滚轴故障诊断方法的研究思路如下：

二、数据集介绍与数据处理

  1、数据集介绍

2、数据集分割与合并

3、数据集分析

三、小波时频图导出

四、CNN网络的构建和测试

1、CNN网络构建

 2、训练参数的优化

4、训练结果

一、研究思路

 1、基于小波时频图和CNN的滚轴故障诊断方法的研究思路如下：

（1）信号采集系统通过轴承试验台采集滚动轴承在不同状态下的振动信号，并通过数据分割与合并  构造信号样本集。

（2）对样本集中的振动信号进行CWT，再生成于信号相对应的时频图。

（3）建立CNN模型，选择一定数量的特征图作为训练样本，对CNN进行训练。

（4）将特征图样本集中的剩余样本作为测试样本，对训练好的CNN进行测试。

二、数据集介绍与数据处理

  1、数据集介绍

       振动数据来源于美国凯斯西储大学的开发数据集。轴承试验装置的示意图如图1-1所示。平台组成：一个1.5KW（2马力）的电动机（图左侧）；一个扭矩传感器/ 译码器（图中间连接处）；一个功率测试计（图右侧）；驱动端轴承为SKF6205 ，采样频率为12Khz和48Khz；电子控制器(图中没显示) 。轴承试验装置的示意图如图1-1所示。平台组成：一个1.5KW（2马力）的电动机（图左侧）；一个扭矩传感器/ 译码器（图中间连接处）；一个功率测试计（图右侧）；驱动端轴承为SKF6205 ，采样频率为12Khz和48Khz；电子控制器(图中没显示) 。

基于小波时频图和2D-CNN的滚动轴承故障检测

图1-1 轴承试验装置示意图

DE - drive end accelerometer data 驱动端加速度数据；

FE - fan end accelerometer data 风扇端加速度数据；

BA - base accelerometer data 基座加速度数据（正常）；

time - time series data 时间序列数据；

RPM- rpm during testing 转每分钟，除以60为旋转频率；

B -滚动体故障；IR – 内圈故障；OR –外圈故障；

驱动端和风扇端轴承外圈的损伤点分别放置在3点钟、6点钟、12点钟三个不同位置。

数据文件为Matlab格式。每个文件包含风扇和驱动端振动数据，以及电机转速。

      12k_Fan_End_OR007@6_3_297.mat这个是12KHZ采样频率下驱动端外圈故障、@6是外圈的损伤点在6点钟位置、3代表电机载荷模式、297是编号。数据集由12kHz采样率下电机端轴承的故障数据（DE）60个、12kHz采样率下风扇端轴承的故障数据（FE）45个、48kHz采样率下电机端轴承的故障数据（DE）52个，还有正常运转的轴承数据4个组成。

基于小波时频图和2D-CNN的滚动轴承故障检测

图1-2  12k驱动端故障数据列表

       如图1-2所示，左边第一列是4种故障尺寸，由于第4种故障缺少外圈故障，所以一般选取前3种尺寸。这样的话内圈、滚动体、外圈3种故障就有3*3也就是9种故障分类，再加上正常状态就是最常见的10分类方式。如果不考虑故障尺寸进行4分类的比较少见，主要是这样的话难度降低不少，结果提升的空间不大，毕竟10分类也没有很麻烦，而且10分类可以直接根据文件来创建数据集，不用进行数据集的合并。外圈故障分为3种，这里我只取取6点钟的，因为不同方向的差异时间上并不大，进行区分也没有那么多必要性。

2、数据集分割与合并

      将采集到的振动信号进行分割，每个样本包含1024个采样点。每种工况下正常状态的样本为350个，不同故障类型下，不同损伤大小及不同工况的样本为117个，最后得到正常状态和三种故障状态的样本各1400个。

表2-1是轴承不同状态的信号样本列表。

表2-1  轴承不同状态的信号样本列表

损伤尺
寸/inch

负载/
hp

电机转速/
（r· min-')

正常
信号

内圈
故障

滚动体
故障

外圈
故障

0

0

1797

350

——

——

——

1

1772

350

——

——

——

2

1750

350

——

——

——

3

1730

350

——

——

——

0.007

0

1797

——

117

117

117

1

1772

——

117

117

117

2

1750

——

117

117

117

3

1730

——

117

117

117

0.014

0

1797

——

117

117

117

1

1772

——

117

117

117

2

1750

——

117

117

117

3

1730

——

117

117

117

0.021

0

1797

——

116

116

116

1

1772

——

116

116

116

2

1750

——

116

116

116

3

1730

——

116

116

116

总计

　

——

1400

1400

1400

1400

       表2-1是我选取论文里的信号样本。这里3个损伤尺寸、3种故障类型。一共3*3=9种故障类型，加上正常，一共10种。4种负载状态，4*10=40。这还不算外圈故障里面损伤点在不同位置的情况。

       这里分类选择可以选择4分类。即正常信号、滚动体故障、内圈故障、外圈故障。

       从分类讲，我找相关的论文大多数是用10分类。因为4分类的指标普遍都太好了，做的人确实也较少，但是实际操作起来也不会简单多少。从工程上讲，区分故障种类的同时我们肯定更希望知道故障的大致严重程度，所以我又做了一个10分类的程序。

表2-2  不同故障样本的数量及标签配置

信号类型

训练样本

测试样本

样本标签

正常信号

900

500

[1 0 0 0 ]

内圈故障

900

500

[0 1 0 0 ]

滚动体故障

900

500

[0 0 1 0 ]

外圈故障

900

500

[0 0 0 1 ]

总计

3 600

2 000

　

基于小波时频图和2D-CNN的滚动轴承故障检测

图2-1 不同故障样本的数量及标签配置

这里选取出的数据集一共是5600个（训练：3600，测试：2000），每一类是1400个（训练：900，测试：500），训练集和测试集的比例大约是0.6428：0.3572。

3、数据集分析

时域分析：在相关论文中，以时域分析作为输入的较多。时域分析中处理起来也很简单，编程容易。

频域分析：频域数据的优点就是提取了频域信息，特征明显。因为是研究信号，还是频域图更加直观，传统的方法里故障频率是至关重  要的指标，而且我们数字信号处理的重点也是频域的内容，但是频域数据丢失了时间特征。

时频图分析：通过小波变换或者短时傅里叶变换形成时频图，将图片作为输入。这里的输入信号是图像。缺点就是数据量大，处理起来  有点慢。这篇论文使用的是连续小波时频图，我用程序来对比一下不同数据的差别。

基于小波时频图和2D-CNN的滚动轴承故障检测

图2-2 不同故障信号进行时域和频域分析对比

 基于小波时频图和2D-CNN的滚动轴承故障检测

图2-3 不同故障信号进行时频图对比

采用cmor3－3小波基对轴承样本数据集中的信号进行CWT，生成时频图。以图2-2的第一列4个信号为例，其时频图如图2-3第一列所示。

matlab的数据分析部分代码

%% 加载原始数据
load 0/12k_Drive_End_B007_0_118;      %读取数据集
a1=X118_DE_time';
load 0/12k_Drive_End_IR007_0_105;
a2=X105_DE_time';
load 0/12k_Drive_End_OR007@6_0_130;
a3=X130_DE_time';
load 0/normal_0_97;
a0=X097_DE_time';%10
L=1024;                                %取数据集里的样本长度为L=864
data0=[]
start_point=randi(length(a0)-L);     %随机取一个起点
end_point=start_point+L-1;
s0=[data0 ;a0(start_point:end_point)];
data1=[]
start_point=randi(length(a1)-L);     %随机取一个起点
end_point=start_point+L-1;
s1=[data1 ;a1(start_point:end_point)];
data2=[]
start_point=randi(length(a2)-L);     %随机取一个起点
end_point=start_point+L-1;
s2=[data2 ;a2(start_point:end_point)];
data3=[]
start_point=randi(length(a3)-L);     %随机取一个起点
end_point=start_point+L-1;
s3=[data3 ;a3(start_point:end_point)];
fs=48000;
y0=fft(s0,length(s0));
f0=fs*(0:length(s0)-1)/length(s0);
y0=fft(s0,length(s0));
t0=0:1/fs:length(s0)/fs;
y1=fft(s1,length(s1));
f1=fs*(0:length(s1)-1)/length(s1);
y1=fft(s1,length(s1));
t1=0:1/fs:length(s1)/fs;
y2=fft(s2,length(s2));
f2=fs*(0:length(s2)-1)/length(s2);
y2=fft(s2,length(s2));
t2=0:1/fs:length(s2)/fs;
y3=fft(s3,length(s3));
f3=fs*(0:length(s3)-1)/length(s3);
y3=fft(s3,length(s3));
t3=0:1/fs:length(s3)/fs;
%绘制时域图
%--------------------------------------------
subplot(421);
plot(s0);
title('正常信号时域图');xlabel('时间/S');ylabel('幅度');
%绘制频谱响应图
%--------------------------------------------
subplot(422);
plot(f0,abs(y0));title('正常信号频谱');xlabel('频率/HZ');ylabel('幅度/dB');
%--------------------------------------------
subplot(423);
plot(s1);
title('B故障信号时域图');xlabel('时间/S');ylabel('幅度');
%绘制频谱响应图
%--------------------------------------------
subplot(424);
plot(f1,abs(y1));title('B0故障信号频谱');xlabel('频率/HZ');ylabel('幅度/dB');
%--------------------------------------------
subplot(425);
plot(s2);
title('IR故障信号时域图');xlabel('时间/S');ylabel('幅度');
%绘制频谱响应图
%--------------------------------------------
subplot(426);
plot(f2,abs(y2));title('IR故障信号频谱');xlabel('频率/HZ');ylabel('幅度/dB');
%--------------------------------------------
subplot(427);
plot(s3);
title('OR故障信号时域图');xlabel('时间/S');ylabel('幅度');
%绘制频谱响应图
%--------------------------------------------
subplot(428);
plot(f3,abs(y3));title('OR故障信号频谱');xlabel('频率/HZ');ylabel('幅度/dB');
N=1024;
fs=12000;
t=0:1/fs:N/fs;
%% 各种变换时频图
wavename='cmor3-3';%cmor是复Morlet小波，其中3－3表示Fb－Fc，Fb是带宽参数，Fc是小波中心频率。
totalscal=256;
Fc=centfrq(wavename); % 小波的中心频率
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
coefs0=cwt(s0,scals,wavename); % 求连续小波系数
coefs1=cwt(s1,scals,wavename); % 求连续小波系数
coefs2=cwt(s2,scals,wavename); % 求连续小波系数
coefs3=cwt(s3,scals,wavename); % 求连续小波系数
figure
subplot(4,3,1)
imagesc(t,f,abs(coefs0)/max(max(abs(coefs0))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('正常小波时频图');
subplot(4,3,4)
imagesc(t,f,abs(coefs1)/max(max(abs(coefs1))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('B故障小波时频图');
subplot(4,3,7)
imagesc(t,f,abs(coefs2)/max(max(abs(coefs2))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('IR故障小波时频图');
subplot(4,3,10)
imagesc(t,f,abs(coefs3)/max(max(abs(coefs3))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('OR故障小波时频图');
% 短时傅里叶变换时频图
% spectrogram(s,256,250,256,fs);
% 时频分析工具箱里的短时傅里叶变换
f = 0:fs/2;
[tfr,t,f] = tfrstft(s0');
tfr = tfr(1:floor(length(s0)/2), :);
subplot(4,3,2)
imagesc(t, f, abs(tfr)/max(max(abs(tfr))));%得到一个彩色图像
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('正常短时傅里叶变换时频图');
%st
[st_fre,st_times,st_frequencies] = st(s0,4);
subplot(4,3,3)
imagesc(t,f,abs(st_fre)/max(max(abs(st_fre))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('正常广义S变换时频图');
%%%%%%%%%%%%%%%%%%%%
f = 0:fs/2;
[tfr,t,f] = tfrstft(s1');
tfr = tfr(1:floor(length(s1)/2), :);
subplot(4,3,5)
imagesc(t, f, abs(tfr)/max(max(abs(tfr))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('B故障短时傅里叶变换时频图');
%st
[st_fre,st_times,st_frequencies] = st(s1,4);
subplot(4,3,6)
imagesc(t,f,abs(st_fre)/max(max(abs(st_fre))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('B故障广义S变换时频图');
%%%%%%%%%%%%%%%%%%%%
f = 0:fs/2;
[tfr,t,f] = tfrstft(s2');
tfr = tfr(1:floor(length(s2)/2), :);
subplot(4,3,8)
imagesc(t, f, abs(tfr)/max(max(abs(tfr))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('IR故障短时傅里叶变换时频图');
%st
[st_fre,st_times,st_frequencies] = st(s2,4);
subplot(4,3,9)
imagesc(t,f,abs(st_fre)/max(max(abs(st_fre))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('IR故障广义S变换时频图');
%%%%%%%%%%%%%%%%%%%%
f = 0:fs/2;
[tfr,t,f] = tfrstft(s3');
tfr = tfr(1:floor(length(s3)/2), :);
subplot(4,3,11)
imagesc(t, f, abs(tfr)/max(max(abs(tfr))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('OR故障短时傅里叶变换时频图');
%st
[st_fre,st_times,st_frequencies] = st(s3,4);
subplot(4,3,12)
imagesc(t,f,abs(st_fre)/max(max(abs(st_fre))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');ylabel('频率 f/Hz');title('OR故障广义S变换时频图');

三、小波时频图导出

把训练集用连续小波变换转换成时频图，对训练集数据进行小波时频提取，并保存在对应的文件夹下   

小波时频图_mat/test_img/“个数”-“类型”.jpg

小波时频图_mat/test_img/“个数”-“类型”.jpg

这里要注意我用MATLAB导出时因为样本有5600个，所以速度很慢，用Python导出时速度就会快一些，但是Python是从0开始的，不管是“个数”还是“类型”，而MATLAB是从1开始。

python的部分代码

import pywt
import matplotlib.pyplot as plt
import numpy as np
from data_preprocess import prepro
path = '0'
x_train, y_train, x_test, y_test = prepro(
    d_path=path,
    length=1024,
    number=100,
    normal=True,
    rate=[0.7,  0.3],
    enc=False, enc_step=28)
for i in range(0, len(x_train)):
    N = 1024
    fs = 12000
    t = np.linspace(0, N/fs, N, endpoint=False)
    wavename = 'cmor3-3'
    totalscal = 256
    fc = pywt.central_frequency(wavename)
    cparam = 2 * fc * totalscal
    scales = cparam / np.arange(totalscal, 1, -1)
    [cwtmatr, frequencies] = pywt.cwt(x_train[i], scales, wavename, 1.0 / fs)
    plt.contourf(t, frequencies, abs(cwtmatr))
    plt.axis('off')
    plt.gcf().set_size_inches(875/96, 656/96)
    plt.gca().xaxis.set_major_locator(plt.NullLocator())
    plt.gca().yaxis.set_major_locator(plt.NullLocator())
    plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
    plt.margins(0, 0)
    plt.savefig('./小波时频图 py/train_img/' + str(i) + '-' + str(y_train[i]) + '.jpg')
    plt.clf()  #避免内存溢出
    plt.close() #释放内存
for i in range(0, len(x_test)):
    N = 1024
    fs = 12000
    t = np.linspace(0, N/fs, N, endpoint=False)
    wavename = 'cmor3-3'
    totalscal = 256
    fc = pywt.central_frequency(wavename)
    cparam = 2 * fc * totalscal
    scales = cparam / np.arange(totalscal, 1, -1)
    [cwtmatr, frequencies] = pywt.cwt(x_test[i], scales, wavename, 1.0 / fs)
    plt.contourf(t, frequencies, abs(cwtmatr))
    plt.axis('off')
    plt.gcf().set_size_inches(875/96, 656/96)  #设置图片大小以英尺为单位
    plt.gca().xaxis.set_major_locator(plt.NullLocator())
    plt.gca().yaxis.set_major_locator(plt.NullLocator())
    plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)  #去除边缘
    plt.margins(0, 0)
    plt.savefig('./小波时频图 py/test_img/' + str(i) + '-' + str(y_test[i]) + '.jpg')
    plt.clf()  #避免内存溢出
    plt.close() #释放内存

四、CNN网络的构建和测试

1、CNN网络构建

CNN网络一共有5个层级结构：输入层、卷积层、激励层、池化层、全连接FC层。论文中的采样层即池化层。

1.输入层：与传统神经网络/机器学习一样，模型需要输入的进行预处理操作。即第二部分数据处理。

2.卷积层：找出图像中的某些特征与过滤器(Filter)进行过滤，过滤的过程就是求卷积的过程。这里过滤器就是卷积核。

3.激励层：所谓激励，实际上是对卷积层的输出结果做一次非线性映射。Sigmoid函数、Tanh函数、ReLU、Leaky ReLU、ELU、Maxout。

4.池化层：池化（Pooling）也称为欠采样或下采样。主要用于特征降维，压缩数据和参数的数量，减小过拟合，同时提高模型的容错性。

5.输出层（全连接层）：经过前面若干次卷积+激励+池化后，终于来到了输出层，模型会将学到的一个高质量的特征图片全连接层。其实在全连接层之前，如果神经元数目过大，学习能力强，有可能出现过拟合。

这里我们通过不同图片的命名来区分标签：x_1.jpg是正常状态，x_2.jpg是内圈故障，x_3.jpg是滚动体故障，x_4.jpg是外圈故障。

输入特征图的大小为28×28（同时用64×64对比）; 隐层由两层卷积层和两层采样层交替组成，第一层和第二层卷积层的卷积核个数分别为6和12，卷积核的大小取为5 ×5，激活函数选择 sigmoid 函数; 采样层（池化层）的采样方式选择均值采样，即对特征图求其p × p区域内的平均值，区域大小取2×2且区域不重叠; 分类器选择softmax 分类器。

网络结构图3-1。

基于小波时频图和2D-CNN的滚动轴承故障检测

图3-1  CNN网络结构

 2、训练参数的优化

采用学习率自适应调整、以训练误差作ρ为停止条件的方法对网络作进一步优化，具体方法如下: 初始学习率设定为0.009，最小学习率设定为0.6，每迭代五次就调整一次学习率，学习率的衰减因子设定为0.97。

如果我们直接用默认的学习率（0.001）和ReLU激活函数效果会更好，收敛的更快。

4、训练结果

基于小波时频图和2D-CNN的滚动轴承故障检测

基于小波时频图和2D-CNN的滚动轴承故障检测

图3-2  训练过程的准确率和误差曲线

这里我们发现第二个图的收敛比第一个图慢，所以第一个图的网络结构参数更优，第一个图输入层28*28 ，第一层和第二层卷积层的卷积核个数分别为6 和12，卷积核的大小取为 5 × 5，激活函数选择 sigmoid 函数。

基于小波时频图和2D-CNN的滚动轴承故障检测

图3-3  测试混淆矩阵

代码下载链接：

matlab版：🍞正在为您运送作品详情

python版：🍞正在为您运送作品详情

标签：

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注