基于滚动轴承加速度信号的故障诊断

一、介绍：改编自一次实验报告，数据为滚动轴承的一维加速度信号，长度为20w+，原码见文末，对于报告中出现的错误欢迎指正

二、实验名称：基于加速度信号的故障诊断实验

三、实验目的：

基于有不同类型损坏（内圈，外圈）的3/4英寸滚动轴承的振动加速度信号，采用人工智能算法与时频域分析的思路完成对信号的处理分类及故障诊断。

四、实验原理：

信号处理是故障诊断领域的一种重要方法，可以通过对系统产生的信号进行处理和分析，确定故障原因与位置，在机械系统损坏早期发出预警以得到及时维修，从而大幅降低财产损失并提高工作效率。常用的方法包括时域分析，频域分析与时频域分析（小波分析）以及近年来流行的机器学习方法。

随机森林算法是一种集成学习方法，它基于多个决策树模型进行分类。在训练过程中，随机森林算法通过随机选择特征和样本，构建多个决策树模型，并将它们集成为一个强分类器。在分类时，随机森林算法通过多数表决的方法，选取最终的分类结果。

傅里叶变换是一种经典的数字信号处理方法，可以在频域图上更好的观察信号的组成成分，从而甄别特殊信号实现故障诊断的目的。

小波变换是一种将信号分解成不同频率子带的方法，其核心思想是采用滤波器组将信号分解成一组不同频率的子信号。小波变换具有良好的时频局部性，因此适用于对信号进行分析和处理。

五、实验内容

5.1数据预处理

我们采集了三种不同类型的加速度数据，并对其进行了标注，包括正常、内圈故障和外圈故障，并将采集的数据进行了滤波和采样处理。在预处理过程中，由于不同数据长度并不一致，我们对数据进行了缺失值处理。考虑到机械结构运转产生的噪声，我们又利用滑移平均的方法进行了异常值处理，使得数据整体更平滑，以保证数据的质量和可靠性。下图即为10s内原始信号（蓝色）与处理后信号（红色）。

5.2随机森林算法对信号的分类尝试

在课堂上我们了解到，与传统的决策树算法相比，随机森林算法可以有效避免过拟合的问题，同时具有较高的准确性和可靠性。随机森林算法的优势在于能够处理高维度数据和大量数据，对于离散型和连续型变量的分类能力也很强。而我们的振动信号便是一个长度逾20w的一维数据，虽不满足高维，但也可做一次大胆的尝试。

随机森林算法的实现过程相对简单，主要分为两个步骤：构建决策树和随机森林集成。在构建决策树的过程中，随机森林算法会随机选择部分特征进行划分，避免了对所有特征的过分依赖，从而减少了过拟合的可能。在随机森林集成的过程中，随机森林算法会通过投票或取平均值等方法，对所有决策树的结果进行集成，从而得到最终的分类结果。

将三组不同类型的轴承数据打上对应类型标签，分为训练集与测试集（8:2）以进行有监督的机器学习，在训练过程中，我们依据数据大致特征手动设置了树的数量、最大深度、最小叶子节点数等参数。在测试过程中，又使用新的数据集对模型进行了验证，并计算了分类准确率和召回率等指标。

但是结果并不如人意，对于此三分类问题，模型平均准确率仅有37%—40%，详细数据如下：

我们又查询了其它资料，将降噪处理从原先的滑移窗口改为减去高斯白噪声以更好的保留不同类型信号特征，更改数据集的规模，都没有取得更好的方式。从数据集提升准确率看似行不通，我们又尝试从森林构造的角度进行分析，学习到了grid-search方法，详细介绍见6.1。

因此，我们认为，一维振动加速度数据集的问题可能在于其特征维度较少，而随机森林算法通常在处理高维度数据时表现更好。对于具有较少特征的数据集，可能会出现过拟合或者欠拟合的情况，导致模型的性能不佳。

5.3 小波分析判定故障类型

小波变换的基本过程可以分为两个步骤：分解和重构。由于在前述已经进行了数据填充与平滑处理，故不再赘述，下面介绍我们进行小波分析的步骤。

选择小波基函数：根据需要分析的信号类型和分析目的选择合适的小波基函数。

Daubechies小波具有以下特点：在时域是有限支撑的，即长度有限。在频域在处有N阶零点。和它的整数位移正交归一，即小波函数可以由所谓“尺度函数”求出来。尺度函数为低通函数，长度有限，支撑域在的范围内。因此常常应用在工程领域当中，我们这里利用db4小波，其形状与morlet小波相似。

进行小波变换：将信号在小波基函数下进行分解，得到多个子信号与每个子信号的能量分布。

图序数为偶数时，则表示经过低通滤波系数 g(k) 分解得到的低频成分信号，反之为奇数时，表示经过高通滤波系数 h(k) 分解得到的高频成分信号，按照上述方式分解，信号在第 i 层小波包分解后，将得到 2^i ( 2的 i 次方 ) 个特征信号，每个特征信号与相应频率段相匹配。

分解结果见下：

Figure2：normal circumstance wavelet result

Figure3：inner circumstance wavelet result

Figure4：outer circumstance wavelet result

可以看到，当轴承正常工作时，其小波能量谱较为集中，第一子波几乎占据了整体能量的百分之九十。面对不同的故障类型所进行的分类效果也是非常出色的：内圈故障反映在图上则为高频信号突出，外圈故障反映在图上则为低频信号明显。我们又重载了平行试验的其他数据集，均很好地完成了分类。

在进行小波分解的过程时，为了获得更精细的频域信息和更准确的故障诊断结果，我们进行了多层分解。多层分解就是将分解得到的低频部分再次进行小波分解，重复进行直至满足分解的需求。查阅资料我们得知，为确定小波分解的层数，可以考虑信号的特点：不同类型的信号具有不同的特点，例如振动信号和电气信号的频率范围和能量分布等不同，需要根据信号特点确定分解层数。

这里我们分解四阶得到16个子波，充分展示了不同信号的差异。

5.4 傅里叶变换提取特征频率

前述工作已经完成了对陌生信号的处理，包括数据处理与类别判断。而故障诊断的时频域分析更大价值在于精细捕捉出现故障的频率值。这里我们想先采用（快速）傅里叶变换得到其频域图，对不同情况的频域进行分析，从而抓住特征谱线判定故障频率。

从图我们可以清楚的看出，在基频以及其对应的倍频位置，均出现了局部极大值，这也是原始信号的频率。此外，由轴承说明书记载，内圈故障特征频率在基频的4.95倍处，外圈则为3.048倍处，

Normal，图示高峰为基频与倍频

Inner，注意到5倍频附近有很多尖峰，且566hz处高度明显高于normal状况下幅值

Outer，注意到3倍频处信号得到加强（但并非很明显）

5.5 希尔伯特变换（包络谱）提取特征频率

由于实验步骤的失误，用傅里叶变换提取特征谱线进行比对的方法无法进行，需要再尝试其它方案。

查阅文献我们得知，传统的故障诊断有经验公式法，利用一些理论推导的公式即可计算出现相应故障的特征频率，这里我们以内圈故障为例。

希尔伯特变换的主要特点是可以将一个实数信号转换为一个复数信号，并且可以提取出信号的振动成分，其它一些特点包括：

1、具有线性性质，可以用于处理任何形式的线性信号。

2、对高频噪声有很好的抑制作用，因为它只保留了信号中的振动成分。

3、可以应用于非平稳信号和周期性信号的分析，因为它不受信号周期的限制。

4、具有相位信息，可以用于分析信号的相位特性。

希尔伯特变换可以将一个实数信号x（t）转换成一个复数信号，具体的转换方法是通过对＄x（t）＄进行卷积得到一个称为希尔伯特核的函数h（t），然后将它和x（t）进行卷积得到XH（t）。

XH（t）是一个复数信号，它的实部是原信号x（t），而虚部则是x（t）的希尔伯特变换后的结果。可以证明，XH（t）的振动部分是由x（t）中所有频率高于0的正弦和余弦分量构成的。因此，通过提取XH（t）的振动部分，就可以分析x（t）中的所有频率高于0的成分。

希尔伯特变换在信号处理中有广泛的应用，特别是在振动分析、滤波、调制等方面。

可以看到，变换后得到的单个频域分析图与理论值有着非常完美的契合。

六、扩展思考与学习

6.1 随机森林算法调优与失败分析

随机森林算法是一种强大的集成学习算法，其在很多数据集上都能够取得良好的效果。然而，在特定的数据集上，随机森林算法的效果可能会不如其他算法。

从适用范围来看，随机森林算法适用于具有以下特点的数据集：

1、数据集具有较多的特征维度。随机森林算法通过随机选择特征和样本来构建决策树，可以有效地处理高维数据集。

2、数据集具有复杂的特征之间的关系。随机森林算法可以通过随机选择特征和样本来减小特征之间的相关性，提高模型的泛化能力。

3、数据集中存在噪声或缺失值。随机森林算法可以通过随机选择特征和样本来降低模型的过拟合风险，提高模型的鲁棒性。

Grid Search是一种自动化地搜索最佳参数组合的方法，它通过穷举搜索所有可能的参数组合来找到最佳的模型参数。

定义参数空间：首先，需要确定需要调整的参数以及它们的取值范围，将它们组合成一个参数空间。

定义评估函数：其次，需要定义一个评估函数，这个函数用来评估每个参数组合的性能，并返回性能的评价指标（如准确率、F1-score等）。

执行搜索：然后，执行Grid Search，搜索所有可能的参数组合，并使用评估函数计算每个参数组合的性能。

寻找最佳参数：最后，寻找具有最佳性能评价指标的参数组合，可以将其作为最终的模型参数。

我们在有限范围内想要穷尽288个可能的组合，但无奈硬件条件不足，计算机运行四小时仍未得到优化结果，因此只得在这里留下伏笔，期待以后可以有机会得到优化的答案。

6.2 一些神经网络的应用

我们在上学期的实践中已经尝试了各种神经网络在信号处理的应用，但当时原理含糊不清，只是套壳使用，由于目前我仍对一些参数的调节不是十分清楚，因此本次实践并没有套用现有的板子来尝试。不过在查阅资料的过程中，我了解到目前第三代神经网络SNN已经初露头角，得到了大批学者的青睐。在西安交通大学孙闯教授的一篇论文中提到，“Spiking neural network(SNN),also known as the third-generation neural network,mimics the dynamics of the biological brain and is more powerful for processing time-series information.As we know,vibration data is a typical time-series data,SNN would have stronger feature extraction potential forit.”我们研读了本篇文章，由于时间较短未成功复现其巧妙的数据处理方法，只是大致一窥，了解了SNN的强大之处，也相信在未来会有更多与其打交道的机会。

6.3 反思总结

本次实验我们体悟了故障诊断的全流程，包括数据采集、平台搭建、数据分析所有环节，与上学期项目课中直接给数据复现论文不同，在亲自实践的过程中才真正体悟到要想做出好看的结果，单单有优秀的分析方法是远远不够的，同样重要的是高质量的数据。在方法的海洋中挑取所谓名称高级的以美化结果是非常简单的，但一个好的工程师应当清晰知道每个环节需要注意的点，博采众法之长为基础，更要以审慎的目光注视数据，挖掘数据价值，这也是我们需要进一步提高的。

七、原码

"""RF
@author: infinitesimal
"""from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn import metrics
from sklearn.metrics import plot_confusion_matrix
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pddata=pd.read_csv
labels =[]
data['normal']=10*(data['normal']-np.mean(data['normal']))
#data['normal'] = moving_average(interval =data['normal'], window_size = 3)
data['inner']=10*(data['inner']-np.mean(data['inner']))
#data['inner']= moving_average(interval =data['inner'], window_size = 3)
data['outer']=10*(data['outer']-np.mean(data['outer']))
#data['outer']= moving_average(interval =data['outer'], window_size = 3)
for i in range(600000):if i <200000:labels.append('1')if 200000<=i <400000:labels.append('2')if 400000<=i <600000:labels.append('3')
num_example = data.shape[0]
data=data.values.reshape(600000,1,order='f')
arr = np.arange(600000)
np.random.shuffle(arr)
new_data=[]
new_labels=[]
for i in range(600000):k=arr[i]new_data.append(data[k])new_labels.append(labels[k])target_name=['1','2','3']
x_train, x_test, y_train, y_test = train_test_split(new_data, new_labels, test_size=0.1, random_state=2)# 创建一个随机森林分类器的实例
randomforest = RandomForestClassifier(n_estimators =100, n_jobs = -1,random_state =50,max_features = "auto", min_samples_leaf = 50)# 利用训练集样本对分类器模型进行训练
randomforest.fit(x_train, y_train)
expected = y_test  # 测试样本的期望输出
predicted = randomforest.predict(x_test)  # 测试样本预测from sklearn.model_selection import GridSearchCV
def evaluate(model, test_features, test_labels):predictions = model.predict(test_features)errors = abs(predictions - test_labels)mape = 100 * np.mean(errors / test_labels)accuracy = 100 - mapeprint('Model Performance')print('Average Error: {:0.4f} degrees.'.format(np.mean(errors)))print('Accuracy = {:0.2f}%.'.format(accuracy))return accuracy
# Create the parameter grid based on the results of random search
param_grid = {'bootstrap': [True],'max_depth': [80, 90, 100, 110],'max_features': [2, 3],'min_samples_leaf': [3, 4, 5],'min_samples_split': [8, 10, 12],'n_estimators': [100, 200, 300, 1000]
}
# Create a based model
rf = RandomForestClassifier()
# Instantiate the grid search model
grid_search = GridSearchCV(estimator = rf, param_grid = param_grid, cv = 3, n_jobs = -1, verbose = 2)# 输出最佳参数组合和最佳得分
# Fit the grid search to the data
grid_search.fit(x_train, y_train)
grid_search.best_params_print("\n")
best_grid = grid_search.best_estimator_
grid_accuracy = evaluate(best_grid, x_train, y_train)print("\n")
print('Improvement of {:0.2f}%.'.format( 100 * (grid_accuracy - 0.39) /0.39))plt.figure(figsize=(8,8),dpi=800)
plot_confusion_matrix(randomforest,x_test,y_test,display_labels=target_name,cmap="Blues",normalize="true",
)
accuracy = metrics.accuracy_score(y_test, predicted)
print("Accuracy: %.2f%%" % (accuracy * 100.0))
plt.title("ECG classification results")
plt.show()

#希尔伯特
import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
import pandas as pd
from scipy import signal, fftpack, stats
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号xt=pd.read_csv
def moving_average(interval, window_size):window = np.ones(int(window_size)) / float(window_size)return np.convolve(interval, window, 'same')fig = plt.figure(figsize=(30, 10),dpi=500)
t = np.linspace(start =0, stop=len(xt), num =len(xt))
average=xt.mean()
y_av = moving_average(interval =xt, window_size = 10)
plt.plot(t, xt, "b.-", t, y_av, "r.-")# step1: 做希尔伯特变换
ht = fftpack.hilbert(xt)
print(ht)t = np.sqrt(ht**2+xt**2)
at =np.sqrt(xt**2 + ht**2)sampling_rate =20000
at = at - np.mean(at)
am = np.fft.fft(at)
am = np.abs(am)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)
freq = freq[0:int(len(freq)/2)]
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.03, colors='r')
plt.vlines(x=156.13*2, ymin=0, ymax=0.03, colors='r')

"""
wavelet@author: infinitesimal
"""import numpy as np
import matplotlib.pyplot as plt
import os
from sklearn import preprocessing
import pywt
import pywt.data
import pandas as pd
import math
import scipy
def get_rms(records):"""均方根值 反映的是有效值而不是平均值"""return math.sqrt(sum([x ** 2 for x in records]) / len(records))data =pd.read_csv
tr =data
wp = pywt.WaveletPacket(data=tr, wavelet='db4',mode='symmetric',maxlevel=4)
#根据频段频率（freq）进行排序
print([node.path for node in wp.get_level(1, 'freq')])
print([node.path for node in wp.get_level(2, 'freq')])
print([node.path for node in wp.get_level(3, 'freq')])
print([node.path for node in wp.get_level(4, 'freq')])aaa = wp['aaa'].data
print(pywt.families())
def wpd_plt(signal,n):#wpd分解wp = pywt.WaveletPacket(data=signal, wavelet='db1',mode='symmetric',maxlevel=n)#计算每一个节点的系数，存在map中，key为'aa'等，value为列表map = {}map[1] = signalfor row in range(1,n+1):lev = []for i in [node.path for node in wp.get_level(row, 'freq')]:map[i] = wp[i].data#作图plt.figure(figsize=(24, 10))plt.subplot(n+1,1,1) #绘制第一个图plt.plot(map[1], color="skyblue")for i in range(2,n+2):level_num = pow(2,i-1)  #从第二行图开始，计算上一行图的2的幂次方#获取每一层分解的node：比如第三层['aaa', 'aad', 'add', 'ada', 'dda', 'ddd', 'dad', 'daa']re = [node.path for node in wp.get_level(i-1, 'freq')]  for j in range(1,level_num+1):plt.subplot(n+1,level_num,level_num*(i-1)+j)plt.plot(map[re[j-1]], color="skyblue") #列表从0开始wpd_plt(tr, 4)
n = 4
re = []  #第n层所有节点的分解系数
for i in [node.path for node in wp.get_level(n, 'freq')]:re.append(wp[i].data)
#第n层能量特征
energy = []
for i in re:energy.append(pow(np.linalg.norm(i,ord=None),2))
#for i in energy:
#    print(i)
# 创建一个点数为 8 x 6 的窗口, 并设置分辨率为 80像素/每英寸
plt.figure(figsize=(20, 7), dpi=80)
# 再创建一个规格为 1 x 1 的子图
# plt.subplot(1, 1, 1)
# 柱子总数
N =16
values = energy
# 包含每个柱子下标的序列
index = np.arange(N)
# 柱子的宽度
width = 0.45
# 绘制柱状图, 每根柱子的颜色为紫罗兰色
p2 = plt.bar(index, values, width, label="num", color="blue")
# 设置横轴标签
plt.xlabel('clusters')
# 设置纵轴标签
plt.ylabel('number of reviews')
# 添加标题
plt.title('Cluster Distribution')
# 添加纵横轴的刻度
plt.xticks(index, ('1', '2', '3', '4', '5', '6', '7', '8','9', '10', '11', '12', '13', '14', '15', '16'))
# plt.yticks(np.arange(0, 10000, 10))
# 添加图例
plt.legend(loc="upper right")
plt.show()for i in range(0,16):a=get_rms(re[i])print(a)

参考文献:

[1] Adele Cutler, D. Richard Cutler & John R. Stevens Bioinformatics ,Random Forests. 24 (18) pp. 2010–2014 (2011).

[2] Z Zhao, T Li, J Wu, C Sun, S Wang, R Yan, X Chen,Deep learning algorithms for rotating machinery intelligent diagnosis: An open source benchmark study.ISA transactions 107, 224-255

基于滚动轴承加速度信号的故障诊断相关推荐

【故障诊断分析】基于matlab小波变换外圈轴承故障诊断【含Matlab源码 1678期】
一.简介在滚动轴承的故障诊断时,传统的频谱分析法通常采用共振解调技术 , 具有良好的效果 ,但当内圈 . 滚动体或多点故障时,解调谱线却很难分辨故障类型小波包是小波理论在信号处理应用领域的又一重大发 ...
基于matlab的信号与系统实例,华南理工大学信号与系统实验基于Matlab的信号处理实例...
第2讲基于Matlab的信号处理实例实验内容(1) (1) 读取给定的3D加速度信号文件,绘出信号波形: 程序源代码: function sy2 fid = fopen('run 100m_TROUS ...
基于matlab的电机故障诊断,基于BP神经网络的电机故障诊断的研究(MATLAB,程序)
基于BP神经网络的电机故障诊断的研究(MATLAB,程序)(课题申报表,任务书,开题报告,中期检查表,外文翻译,论文25600字,程序,答辩PPT) 摘要随着我国经济发展的不断提高,电气化的应用越 ...
基于UDS的汽车通信故障诊断机制与处理策略
基于UDS的汽车通信故障诊断机制与处理策略 2015-06-10 22:12:34| 分类: AUTOSAR_|举报|字号订阅下载LOFTER 我的照片书 | 摘要:阐述一种诊断控制单元之间通 ...
matlab方波假频现象分析,基于MATLAB 的信号时域采样及频率混叠现象分析
龙源期刊网 http://www.doczj.com/doc/7d04774eae45b307e87101f69e3143323968f58a.html 基于MATLAB 的信号时域采样及频率混叠现象 ...
基于matlab的信号去噪,基于matlab的信号去噪研究.doc
基于matlab的信号去噪研究.doc 学生毕业设计(论文) 课题名称基于 MATLAB 的信号去噪研究姓名陈文学号 051220206 院系物电系专业电子信息工程 ...
基于matlab的信号合成与分解,基于matlab的信号合成与分解
基于matlab的信号合成与分解 - I - 摘要为了便于进行周期信号的分析与处理,常要把复杂的周期信号进行分解,即将周期信号分解为正余弦等此类基本信号的线性组合,通过对这些基本信号单元在时域和 ...
基于matlab的回波,基于MATLAB回波信号产生与消除.doc
基于MATLAB回波信号产生与消除摘要 MATLAB可以进行矩阵运算.绘制函数和数据.实现算法.创建用户界面.连接其他编程语言的程序等,主要应用于工程计算.控制设计.信号处理与通讯.图像处理.信号 ...
matlab已知ft求f2t,基于Matlab的信号与系统实验指导2
基于Matlab 的信号与系统实验指导实验一连续时间信号在Matlab 中的表示一.实验目的 1.学会运用Matlab 表示常用连续时间信号的方法 2.观察并熟悉这些信号的波形和特性二.实验原 ...

基于滚动轴承加速度信号的故障诊断

基于滚动轴承加速度信号的故障诊断相关推荐

最新文章

热门文章