本期给大家推荐南卫理公会大学公开的齿轮箱数据集,含有三种齿轮故障,并且公布齿轮啮合频率。有写论文需要的同学们赶紧用起来吧!
数据集是获取齿轮箱实验台(见图1和图2)的径向振动信号,测试了三个不同的齿轮健康状态:健康齿轮、1个齿轮出现缺损和3个齿轮出现磨损。当使用数据集时,请引用下列出版物。
数据集已在MATLAB中记录,并以.mat格式存储。在MATLAB中通过加载命令打开文件。所有测量的通道都以伏特为单位。信号进行去趋势处理,以消除加速度计中的偏置。
实验中包括一个模拟数字控制器(研华PCI-1710,12位,100kHz)采样频率为10 kHz的加速度计(模拟装置,ADXL210JQC),用于测量振动信号。采集器是用实时MATLAB工作站进行的。
加速度计配置:Vref=5V,ADXL210 灵敏度=100mV/g,AD单位数=Vref/4096=5/4096=1.22mV,分辨率=AD单位数/ADXL210的灵敏度=1.22mV/100 mV/g=0.0122g=12.2mg。
## 导入包
from matplotlib import pyplot as plt
from matplotlib import rcParams
import numpy as np
import pandas as pd
import os
config = {
"font.family": 'serif', # 衬线字体
"font.size": 14, # 相当于小四大小
"font.serif": ['SimSun'], # 宋体
"mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体,和Times New Roman差别不大
'axes.unicode_minus': False # 处理负号,即-号
}
rcParams.update(config)
##========绘制时域信号图========##
def plt_time_domain(arr, fs=1600, ylabel='Amp($m/s^2$)', title='原始数据时域图', img_save_path=None, x_vline=None, y_hline=None):
"""
:fun: 绘制时域图模板
:param arr: 输入一维数组数据
:param fs: 采样频率
:param ylabel: y轴标签
:param title: 图标题
:return: None
"""
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
plt.rcParams['axes.unicode_minus'] = False # 显示负号
font = {'family': 'Times New Roman', 'size': '20', 'color': '0.5', 'weight': 'bold'}
plt.figure(figsize=(12,4))
length = len(arr)
t = np.linspace(0, length/fs, length)
plt.plot(t, arr, c='g')
plt.xlabel('t(s)')
plt.ylabel(ylabel)
plt.title(title)
if x_vline:
plt.vlines(x=x_vline, ymin=np.min(arr), ymax=np.max(arr), linestyle='--', colors='r')
if y_hline:
plt.hlines(y=0.2, xmin=np.min(t), xmax=np.max(t), linestyle=':', colors='y')
#===保存图片====#
if img_save_path:
plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight')
plt.show()
##========绘制频域信号图========##
def plt_fft_img(arr, fs, ylabel='Amp(mg)', title='频域图', img_save_path=None, vline=None, hline=None, xlim=None):
"""
:fun: 绘制频域图模板
:param arr: 输入一维时域数组数据
:param fs: 采样频率
:param ylabel: y轴标签
:param title: 图标题
:return: None
"""
# 计算频域幅值
length = len(arr)
t = np.linspace(0, length/fs, length)
arr = arr - np.mean(arr)
fft_result = np.fft.fft(arr)
fft_freq= np.fft.fftfreq(len(arr), d=t[1]-t[0]) # FFT频率
fft_amp= 2*np.abs(fft_result)/len(t) # FFT幅值
# 绘制频域图
plt.figure(figsize=(12,4))
plt.title(title)
plt.plot(fft_freq[0: int(len(t)/2)], fft_amp[0: int(len(t)/2)], label='Frequency Spectrum', color='b')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅值')
plt.legend()
if vline:
plt.vlines(x=vline, ymin=np.min(fft_amp), ymax=np.max(fft_amp), linestyle='--', colors='r')
if hline:
plt.hlines(y=hline, xmin=np.min(fft_freq), xmax=np.max(fft_freq), linestyle=':', colors='y')
#===保存图片====#
if img_save_path:
plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight')
if xlim: # 图片横坐标是否设置xlim
plt.xlim(xlim)
plt.tight_layout()
plt.show()
def plt_envelope_spectrum(data, fs, ylabel='Amp(mg)', title='包络谱图', img_save_path=None, vline=None, hline=None, xlim=None):
'''
fun: 绘制包络谱图
param data: 输入数据,1维array
param fs: 采样频率
param xlim: 图片横坐标xlim,default = None
param vline: 图片垂直线,default = None
'''
from scipy import fftpack
#=========做希尔伯特变换=======#
xt = data
ht = fftpack.hilbert(xt)
at = np.sqrt(xt**2+ht**2) # 获得解析信号at = sqrt(xt^2 + ht^2)
at = at - np.mean(at) # 去直流分量
fft_amp = np.fft.fft(at) # 对解析信号at做fft变换获得幅值
fft_amp = np.abs(fft_amp) # 对幅值求绝对值(此时的绝对值很大)
fft_amp = fft_amp/len(fft_amp)*2
fft_amp = fft_amp[0: int(len(fft_amp)/2)] # 取正频率幅值
fft_freq = np.fft.fftfreq(len(at), d=1 / fs) # 获取fft频率,此时包括正频率和负频率
fft_freq = fft_freq[0:int(len(fft_freq)/2)] # 获取正频率
# 绘制包络谱图
plt.figure(figsize=(12,4))
plt.title(title)
plt.plot(fft_freq, fft_amp, color='b')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅值')
if vline:
plt.vlines(x=vline, ymin=np.min(fft_amp), ymax=np.max(fft_amp), linestyle='--', colors='r')
if hline:
plt.hlines(y=hline, xmin=np.min(fft_freq), xmax=np.max(fft_freq), linestyle=':', colors='y')
#===保存图片====#
if img_save_path:
plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight')
if xlim: # 图片横坐标是否设置xlim
plt.xlim(xlim)
plt.tight_layout()
plt.show()
def calculate_kurtosis(signal):
"""
计算信号的峭度
:param signal: 输入信号
:return: 峭度值
"""
x = signal
X_rms2 = np.sum(x**2)/len(x) # 3.均方值
X_rms = np.sqrt(X_rms2) # 4.均方根值(有效值)
X_beta = np.mean( np.power(x, 4) ) # 12.峭度
X_kf = X_beta/X_rms ** 4 # 18.峭度指标
return X_kf
有3个文件,从上至下分别是齿轮缺损、健康、齿轮磨损。
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
data = scio.loadmat(file_path)
data
可知振动数据在'acc'里,且振动数据是单轴的。
def data_read(file_path):
"""
:fun: 读取csv数据
:param file_path: 文件路径
:return df:
"""
data = scio.loadmat(file_path)
acc_arr = data['acc'].flatten()
return acc_arr
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
acc_arr
2.2.2 时域分析
fs = 10000
##=======健康=====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_no_fault_full_load_01_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'no_fault'
plt_time_domain(acc_arr, fs=fs, title=file_name)
kurtosis_value = calculate_kurtosis(acc_arr)
print(f"Kurtosis of the signal: {kurtosis_value}")
##=====断齿====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'chipped_tooth'
plt_time_domain(acc_arr, fs=fs, title=file_name)
kurtosis_value = calculate_kurtosis(acc_arr)
print(f"Kurtosis of the signal: {kurtosis_value}")
##=====齿轮磨损====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_three_worn_teeth_full_load_13_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'worn_teeth'
plt_time_domain(acc_arr, fs=fs, title=file_name)
kurtosis_value = calculate_kurtosis(acc_arr)
print(f"Kurtosis of the signal: {kurtosis_value}")
Kurtosis of the signal: 2.474412899848076
Kurtosis of the signal: 3.0782704766646565
Kurtosis of the signal: 3.0312783300070545
2.2.3 频域分析
fs = 10000
##=======健康=====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_no_fault_full_load_01_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'no_fault'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[365*2, 365*3, 365*4, 365*5, 365*6, 365*7, 365*8, 365*9], xlim=None)
##=====断齿====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'chipped_tooth'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[365, 365*2, 365*3, 365*4, 365*5, 365*6, 365*7, 365*8, 365*9])
##=====齿轮磨损====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_three_worn_teeth_full_load_13_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'worn_teeth'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[365, 365*2, 365*3, 365*4, 365*5, 365*6, 365*7, 365*8, 365*9])
可见明显的四个倍频,主频为齿轮啮合频率。下面在300-500Hz范围内放大看主频及边频。
fs = 10000
##=======健康=====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_no_fault_full_load_01_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'no_fault'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[355, 365,365*2], xlim=[300,500])
##=====断齿====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'chipped_tooth'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[355, 365, 365*2], xlim=[300,500])
##=====齿轮磨损====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_three_worn_teeth_full_load_13_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'worn_teeth'
plt_fft_img(acc_arr, fs=fs, title=file_name, vline=[355, 365, 365*2], xlim=[300,500])
分析:
齿轮故障频谱分析,主要是看边频幅值是否升高,可以看到故障齿轮的边频明显增多。
2.2.4 包络谱分析
fs = 10000
##=======健康=====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_no_fault_full_load_01_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'no_fault'
plt_envelope_spectrum(acc_arr, fs=fs, vline=[365, 365*2], xlim=[0, 1000], title=file_name)
##=====断齿====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'chipped_tooth'
plt_envelope_spectrum(acc_arr, fs=fs, vline=[365, 365*2], xlim=[0, 1000], title=file_name)
##=====齿轮磨损====##
file_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_three_worn_teeth_full_load_13_December_2009_10kHz_pos1.mat'
acc_arr = data_read(file_path)
fs = 10000
file_name = r'worn_teeth'
plt_envelope_spectrum(acc_arr, fs=fs, vline=[365, 365*2], xlim=[0, 1000], title=file_name)
分析:
可以看到齿轮啮合频率,这是因为啮合频率被转频调制了。
编辑:李正平
校核:陈凯歌、赵栓栓、董浩杰、曹希铭、赵学功、白亮、陈少华
该文资料搜集自网络