目录

一、研究思路

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；电子控制器(图中没显示) 。

图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个组成。

图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

　

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

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

3、数据集分析

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

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

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

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

图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';%10L=1024;                                %取数据集里的样本长度为L=864data0=[]
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 prepropath = '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 = 1024fs = 12000t = np.linspace(0, N/fs, N, endpoint=False)wavename = 'cmor3-3'totalscal = 256fc = pywt.central_frequency(wavename)cparam = 2 * fc * totalscalscales = 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 = 1024fs = 12000t = np.linspace(0, N/fs, N, endpoint=False)wavename = 'cmor3-3'totalscal = 256fc = pywt.central_frequency(wavename)cparam = 2 * fc * totalscalscales = 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。

图3-1  CNN网络结构

2、训练参数的优化

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

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

4、训练结果

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

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

图3-3  测试混淆矩阵

注意matlab版本要有深度学习工具箱，我的是2022a版本，python要用tensorflow2.6

代码下载链接：

matlab版：

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

  1. 基于FFT频谱与小波时频图的双流CNN轴承故障诊断模型

    在博客https://blog.csdn.net/qq_41043389/article/details/106150980里,我们提出了采用小波时频图作为轴承信号的故障特征数据,即首先利用提取各样本 ...

  2. Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”

    Python实现"EMD\EEMD\VMD+Hilbert时频图"与"CWT小波时频图"   信号处理中常需要分析时域统计量.频率成分,但不平稳信号的时域波形往 ...

  3. 融合DE 端和FE端数据，利用小波变换生成时频图，再分别利用DCNN、KNN和DNN进行对比实验(python代码)

    1.数据集介绍: 试验台如图所示,试验台左侧有电动机,中间有扭矩收集器,右侧有动力测试仪,控制电子设备在图中没有显示.SKF6203轴承使用16通道数据采集卡采集轴承的振动数据,并在驱动端部分(DE) ...

  4. 基于小波精英解学习和多角度搜索的新型阴阳平衡优化算法

    文章目录 一.理论基础 1.基本阴阳平衡优化算法 (1)基于档案集的解更新阶段 (2)基于超球体的解更新阶段 (3)阴阳平衡优化(YYPO)算法伪代码 2.基于小波精英解学习和多角度搜索的新型阴阳平衡 ...

  5. matlab模式识别提取特征向量,一种基于小波特征向量提取的手机检测方法与流程...

    本发明涉及到手机检测领域,尤其涉及到一种基于小波特征向量提取的手机检测方法. 背景技术: 随着保密要求的不断提高,很多场合严禁携带手机.录音笔.录像机等电子产品,亟需一种设备可以检测出该类电子产品.目 ...

  6. 基于小波的图像边缘检测,小波变换边缘检测原理

    1.什么是"小波神经网络"?能干什么用呀 小波神经网络(Wavelet Neural Network, WNN)是在小波分析研究获得突破的基础上提出的一种人工神经网络.它是基于小波 ...

  7. 002基于小波的多类癫痫类型分类系统-2021

    WAVELET-BASED MULTI-CLASS SEIZURE TYPE CLASSIFICATION SYSTEM ABSTRACT 癫痫是最常见的脑部疾病之一,影响着全球1%以上的人口.它的特 ...

  8. bp神经网络预测未来五年数据_基于小波神经网络的数据中心KPI预测

    随着软件和微服务的发展,智能运维越来越受到人们的重视.在大量的运维数据里,最不可忽视的就是各种关键性能指标数据(Key Performance Indicators,KPI),它们在数学上都可以被表达 ...

  9. 基于小波图像去噪的MATLAB实现

    基于小波图像去噪的MATLAB实现 一.课题背景 数字图像处理(Digital Image Processing,DIP)是指用计算机辅助技术对图像信号进行处理的过程.数字图像处理最早出现于 20世纪 ...

最新文章

  1. 《Greenplum企业应用实战》一2.3　畅游Greenplum
  2. 在网页中怎样打印网页中的一部分（比如打印网页中的一个表格）
  3. CES线下展回归在即：飞行汽车外骨骼智能戒指吸足眼球
  4. Golang 协程goroutine的调度模型-MPG模式
  5. 三维位姿：***图像特征-特征提取-姿态估计
  6. c# mysql app.config_60. C# -- 读取 appconfig文件配置数据库连接的方法
  7. apue.3e环境配置
  8. 微信小程序或微信网页里关注公众号
  9. 适合程序员的护眼显示器——大上科技Paperlike系列电子墨水显示器
  10. 自己写了个磁力链搜索引擎
  11. pci串口驱动安装失败_PCI并口卡驱动安装不上
  12. 2019开源BI软件排行榜
  13. 从小白开始教你怎样在Eclipse中使用Git(番外) - 各种图标的含义
  14. TPH-YOLOv5: Improved YOLOv5 Based on Transformer Prediction Head forObject Detection on Drone-captur
  15. JQuery Ajax 参数含有特殊字符
  16. 最新引流脚本之窃语漂流瓶引流脚本，如何使用窃语脚本
  17. 可以写一个表白代码吗
  18. 《等着我吧，我会回来》 苏·西蒙诺夫
  19. 建造者（Builder）模式
  20. 投资入门第 3 步：技术分析法(常用技术分析)

热门文章

  1. 阿里开源性能测试神器，性能监控分析工具Arthas
  2. Camtasia快捷键大全
  3. CDH安装失败，如何重新安装
  4. 【BZOJ】1617: [Usaco2008 Mar]River Crossing渡河问题（dp）
  5. CentOS 9 开局配置
  6. 禁止访问 (403)CSRF验证失败. 请求被中断.
  7. 小偏方！不看，走宝.
  8. 三种循环的流程图画法总结 （转载）
  9. java int.tryparse_【转载】 C#中使用int.TryParse方法将字符串转换为整型Int类型
  10. php的构造函数理解