本期给大家推荐南卫理公会大学公开的齿轮箱数据集，含有三种齿轮故障，并且公布齿轮啮合频率。有写论文需要的同学们赶紧用起来吧！目录1 数据集描述2 时、频域分析 2.1 导入包及函数 2.2 SMU齿轮故障分析 2.2.1 数据读取 2.2.2 时域分析 2.2.3 频域分析 2.2.4 包络谱分析1 数据集描述数据集是获取齿轮箱实验台(见图1和图2)的径向振动信号，测试了三个不同的齿轮健康状态：健康齿轮、1个齿轮出现缺损和3个齿轮出现磨损。当使用数据集时，请引用下列出版物。[1] A.H. Zamanian and A. Ohadi, “Gear fault diagnosis based on gaussiancorrelation of vibrations signals and wavelet coeffcients,” Applied SofComputing,vol.1l,no.8, pp.4807-4819,2011.[2]A.H. Zamanian and A. Ohadi,“Gearbox fault detection through PSOexact wavelet analysis and SVM classifer,, in 18th Annual InternationalConference on Mechanical Engineering-ISME2010, 2010. 图1 实验装置图图2 实验装置示意图故障齿轮如图3所示，图4演示了加速度计安装方向和位置。数据集可以从下面的链接下载：https://goo.gl/TorZJq。小齿轮(主动轮，15个齿)和大齿轮(从动轮，110个齿)，在变速箱中提供7.33的减速比。如图3所示，变速箱在小齿轮三种不同的健康状态下进行了测试：健康，1个齿轮缺损，3个齿轮磨损。小齿轮以1420rpm的额定转速旋转。在额定转速下，齿轮啮合频率(GMF)为365Hz；然而，精确的GMF可以通过观察不同齿轮健康状态下所获取信号FFT中的靠近365Hz的主频来计算。在每一种情况下，加速度计采集持续时间为10秒。 图3 齿轮磨损(左)和断齿(右) 图4 加速度计的方向和位置数据集已在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。2 时、频域分析以下程序适合在jupyter notebook编辑器上运行。 2.1 导入包及函数## 导入包from matplotlib import pyplot as pltfrom matplotlib import rcParamsimport numpy as npimport pandas as pdimport osconfig = { "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_kf2.2 SMU齿轮故障分析有3个文件，从上至下分别是齿轮缺损、健康、齿轮磨损。2.2.1 读取数据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_arrfile_path = r'E:\03-公开数据集\南卫理公会大学齿轮数据集(SMU)\数据/Gearbox_a_chipped_tooth_full_load_03_December_2009_10kHz_pos1.mat'acc_arr = data_read(file_path)acc_arr2.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 = 10000file_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 = 10000file_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 = 10000file_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.474412899848076Kurtosis of the signal: 3.0782704766646565Kurtosis 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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_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 = 10000file_name = r'worn_teeth'plt_envelope_spectrum(acc_arr, fs=fs, vline=[365, 365*2], xlim=[0, 1000], title=file_name)分析：可以看到齿轮啮合频率，这是因为啮合频率被转频调制了。编辑：李正平校核：陈凯歌、赵栓栓、董浩杰、曹希铭、赵学功、白亮、陈少华该文资料搜集自网络来源：故障诊断与python学习
