都灵理工大学滚动轴承振动数据集(含轴承、变转速)
本期给大家推荐都林理工大学公开的轴承数据集，含有不同转速、不同故障程度。有写论文需要的同学们赶紧用起来吧！
目录
摘要
1 实验台描述
本节的各部分主要介绍试验台、传感器和采集系统。目的是准确、全面地说明在各种条件下以及在不同故障轴承的情况下所进行的试验台和测量。
1.1 实验台
试验台如图1a所示。基本上包括在一个高速主轴，用于驱动轴的旋转。
主轴的轴承(其特性未知)是润滑的，通过液体(乙二醇/水)制冷回路进行冷却。主轴的本体固定在一个单一且极其刚性的支架上，该支架位于一个非常巨大的钢制底板上。主轴的速度是通过变频器的控制面板设置的，但不能主动控制：主轴不仅没有键相传感器或转速计来检测其实际速度，而且也没有反馈给变频器的控制器。其直接后果是，轴的实际转速总是低于理想转速，并且随着加载的增加，其转速差值增加。
两个相同的滚子轴承的外圈通过同样的轴承座进行支撑(位置B1和B3)，这些轴承的内圈连接到一个非常短且带有滴漏的空心轴上，该轴特别设计用于最高达35000rpm的转速(图1c)。
轴最初是一个完整的齿轮箱的一部分，并携带一个正齿轮，用于驱动旋转。由于施加的扭矩，直齿轮产生径向和切向接触力，这种接触力作用在一对滚子轴承上。在试验台上，由直齿轮产生的径向力被第三个滚子轴承施加的负载所取代(B2，图1b/c所示)。轴承B2的外圈与一个精密加载部件相连，其运动方向与轴正交：当加载部件通过螺母的旋转被拉动时，两个平行的弹簧(图1a中绿色)被压缩起来并产生所需负载。静态力传感器能够测量产生的力，其方向是纯径向的。
专门为这种高速航空应用制造的三种轴承的主要几何参数列于表1。通过向空心轴注入油来保证轴承B1，B2和B3的润滑；然后通过小径向管将离心油(稳定状态工作速度至少为6000rpm)输送到内圈并随后输送到滚子。
表1 滚子轴承主要参数
1.2 传感器
试验台在其工作寿命期间数据收集一直主要集中在轴承支架和主轴上的加速度数据。第3节描述的文件包含图1b中结构的两个最重要点A1和A2的加速度，分别位于故障轴承B1的支座上和轴承B2的支承上。加速度计为三轴IEPE型，频率范用为1-12000Hz(振幅±5%，相位±10°)，标称谐振频率为55kHz，标称灵敏度为1mV/ms-2.
中心轴承上的径向力通过力传感器进行测量，该传感器的灵敏度为-0.499 mV/N，图2b已经通过重复的负载循环测量，示出了几乎为零的滞后循环。
最后，尽可能地将K型热电偶和分辨率为0.1°的专用数字温度计放置在故障轴承的外圈附近。温度已在"耐久性"测试期间手动记录(第3.2节)并从上午和每次测试开始前约20°C升至当天结束时约60°C-70°C
1.3 数据采集
用OROS公司生产的OR38信号分析仪实现了数字化数据采集，其输入通道精度为：相位±0.02°，振幅±0.02dB，频率±0.005%。
模拟-数字转换采用24位Δ-Σ转换器进行，在所有通道上均同步；每个通道的范围分别设置在最小(±17 mV)和最大(±40V)之间，以避免信道饱和，同时达到最佳幅值分辨率。
加速度计的灵敏度被设置在 OR38 中，以便附录A中所列的文件包含以m/s2为单位的加速度时间历史。六个记录通道对应于放置在点A1和A2的加速度计的输出(图1b)；表2指定了每个通道的加速度方向。
表2 测得的加速度的方向(见图2b)
每个数据集具有相同的持续时间T和相同的采样频率fs。采样频率不得除以信道的数目，因此每个时域加速度由T × fs个样本点组成。
T和fs的值在第3节中给出，对轴承进行的两次测试是不同的。
2 测试说明
本文说明了两个不同的实验阶段。在第一个实验中，加速度是在不同故障程度、不同速度和不同负载下运行的；第二个实验中是单个故障轴承在恒速和负载下接受长期(约330小时)测试。
2.1 不同速度和负载测试
每个轴承，从0A到6A，都经历了相同的测试过程：
在最低转速(100 Hz)和无负载下短暂运行，以检查安装是否正确； 施加静载荷：首先是1000N，然后是1400N，最后是1800N； 将轴的速度从0Hz增加到500Hz，步长为100 Hz; 一旦达到轴的稳定速度，立即测量加速度。
整个过程耗时约30分钟，由于逆变器的功率有限，在较高的负载条件下无法达到较高的转速。表4给出速度-载荷组合的列表。
表4 速度-载荷组合的列表
在采样频率fs=51200Hz的条件下，收集了6个通道(表2)的时间，时间长度为数据被记录在文件中，文件名的格式如下：CnA_fff_vvv_m.mat。
● C：根文件名，通用于所有文件；
● n：从0到6的整数值，表示缺陷的种类，例如0A、1A、…、6A(表3)；
● fff：从100到500的整数值，表示轴的标称速度(赫兹)；
● vvv：与称重传感器电压(mV)相对应的整数值，表示所施加的负载；
● M：整数值，表示是否重复测量(m=2)或不重复(m=1)；
● .mat：matlab文件扩展名。
例如，C4A_100_702_1.mat包含由轴承4A(滚子故障)在标称转速为100Hz(即6000rpm)下产生的加速度信号的第一次测试，该加速度信号在约1407N的负载下产生，该负荷对应于负荷单元702mV的输出。(灵敏度0.499mV/N)。
每个文件都包含一个与文件名称相同的矩阵(除了.m扩展名)，该矩阵有512000行(时间示例)和6列(每个通道一个)。所有文件名的列表在表A1中给出。
表A1 八个不同故障程度的轴承的文件名，从0A到6A
2.2 加速寿命测试
对4A轴承(滚子故障)进行了加速寿命测试，以监测缺陷的演变，最初是锥形压痕，最大直径约450μm(图3a)。必须强调的是，在开始这项试验之前，轴承已经安装在完整的变速箱上，在各种负载和速度条件下运行，由于保密协议，本文没有记录。在这第一次测试结束时(约19小时)的状况在图3b中可以看出压痕的边界不像以前那样有规律，并且从滚子表面上除去了少量的材料。在加速寿命试验70、144和332小时后，还检查了缺陷的演变情况，并用图3所示的图片加以说明，遗憾的是，这些照片并不那么统一，因为它们是由不同的人用不同的显微镜拍摄的，但趋势是明显的：随着时间的推移，较大的部分材料被移除，同时，凹痕区域周围的小“坑”被碾平了。
值得注意的是，在开始耐久性试验之前，润滑油已经更换。新油的粘度与以前的油几平相同(40°时15 cSt，ASTM D 445)，但它不是专门为高速应用而设计的。例如，其产生泡沫的阻力不符合所需的航空标准。与原油相反，它具有无神经毒性的优点，这是我们实验室实验的强制性要求。
加速寿命评估的每次测量都是在相同的条件下进行的：额定速度：300Hz(即18000转/分钟)，额定负荷：1800N。
定期记录两个三轴加速度计的时间历程，每次测试持续8秒，采样频率为102.4kHz。
温度的指示由安装在轴承外圈附近尽可能靠近的热电偶发出。每天，在测试开始和打开主轴之前，试验台和轴承的温度约为18-20°C(室温)；由于质量结构和润滑油的热惯性(约15L)，轴承加热持续了数小时。如附录A表A2所示，只有大约三小时后才能达到几乎稳定的温度，因此每个工作日的第一份文件所对应的温度比最终值低5-6°℃。
所有文件名的格式为：E4A_nn.mat
E4A：根文件名，通用于所有文件; nnn：007到226之间的整数值，表示采集的数量。
例如，E4A_083.mat包含测量到一天结束，也就是运行7.5小时后记录的加速度信号。
每个文件都包含一个与文件名称相同的矩阵(除了.m扩展名)，该矩阵有819200行(时间示例)和6列(每个通道一个)。
附录A表A2列出了总共66个文件的清单，并列出了相应的温度和获取时间，这些时间是早上开关主轴的时间。
3 时频域分析
3.1.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 pltplt.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: # 图片横坐标是否设置xlimplt.xlim(0, xlim)plt.tight_layout()plt.show()
def plt_stft_img(arr, fs, ylabel='Amp(mg)', title='stft时频域图', img_save_path=None, vline=None, hline=None, xlim=None):""":fun: 绘制stft时频域图模板:param arr: 输入一维时域数组数据:param fs: 采样频率:param ylabel: y轴标签:param title: 图标题:return: None"""import scipy.signal as signalimport numpy as npimport matplotlib.pyplot as pltf, t, nd = signal.stft(arr, fs=fs, window='hann', nperseg=128, noverlap=64,nfft=None,detrend=False, return_onesided=True, boundary='odd', padded=False, axis=-1)# fs:时间序列的采样频率, nperseg:每个段的长度，默认为256(2^n) noverlap:段之间重叠的点数。如果没有则noverlap=nperseg/2#window ：字符串或元组或数组，可选需要使用的窗。# #如果window是一个字符串或元组，则传递给它window是数组类型，直接以其为窗，其长度必须是nperseg。# 常用的窗函数有boxcar，triang，hamming， hann等，默认为Hann窗。#nfft ：int，可选。如果需要零填充FFT，则为使用FFT的长度。如果为 None，则FFT长度为nperseg。默认为无# detrend ：str或function或False，可选# 指定如何去除每个段的趋势。如果类型参数传递给False，则不进行去除趋势。默认为False。# return_onesided ：bool，可选# 如果为True，则返回实际数据的单侧频谱。如果 False返回双侧频谱。默认为 True。请注意，对于复杂数据，始终返回双侧频谱。# boundary ：str或None，可选# 指定输入信号是否在两端扩展，以及如何生成新值，以使第一个窗口段在第一个输入点上居中。# 这具有当所采用的窗函数从零开始时能够重建第一输入点的益处。# 有效选项是['even', 'odd', 'constant', 'zeros', None].# 默认为‘zeros’,对于补零操作[1, 2, 3, 4]变成[0, 1, 2, 3, 4, 0] 当nperseg=3.# padded：bool，可选# 指定输入信号在末尾是否填充零以使信号精确地拟合为整数个窗口段，以便所有信号都包含在输出中。默认为True。# 填充发生在边界扩展之后，如果边界不是None，则填充为True，默认情况下也是如此。# axis ：int，可选# 绘制STFT时频域图plt.figure(figsize=(12,4))plt.pcolormesh(t, f, np.abs(nd), vmin = np.min(np.abs(nd)), vmax = np.max(np.abs(nd)))plt.title(title)plt.xlabel('时间（t）')plt.ylabel('频率 (Hz)')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: # 图片横坐标是否设置xlimplt.xlim(0, xlim)plt.tight_layout()plt.show()
def data_read(file_path):""":fun: 读取xls数据:param file_path: 文件路径:return df:"""file_name = file_path.split('/')[-1].split('.')[0]arr = scio.loadmat(file_path)[file_name] # 含有6列数据df = pd.DataFrame(arr)df.columns = ['A1_x', 'A1_y', 'A1_z', 'A2_x', 'A2_y', 'A2_z']return dffile_path = r'E:/03-公开数据集/都灵理工大学变转速轴承数据集/都灵理工大学数据集/VariableSpeedAndLoad/C4A_400_496_1.mat'df = data_read(file_path)acc_arr= df['A1_z']df
3.1.2 绘制时域图
fs = 51200file_name = r'C4A_400_496_1.mat'##=====绘制时域数据====##time_img_save_path = file_path.replace('.mat', '时域_.png')plt_time_domain(acc_arr, fs=fs, title=file_name, img_save_path=time_img_save_path)
3.1.3 绘制频域图
fs = 51200file_name = r'C4A_400_496_1.mat'##=====绘制频域数据====##fft_img_save_path = file_path.replace('.mat', '频域_.png')plt_fft_img(acc_arr, fs=fs, title=file_name, img_save_path=fft_img_save_path)
3.1.4 绘制时频域图
fs = 51200file_name = r'C4A_400_496_1.mat'##=====绘制STFT时频域数据====##stft_img_save_path = file_path.replace('.mat', '时频域_.png')plt_stft_img(acc_arr, fs=fs, img_save_path=stft_img_save_path)
3.2 包络谱分析
3.2.1 计算故障特征频率
def bearing_fault_freq_cal(n, d, D, alpha, fr=None):'''基本描述：计算滚动轴承的故障特征频率详细描述：n, fr, d, D, alphareturn C_bpfi, C_bpfo, C_bsf, C_ftf, fr外圈 滚针 保持架 转速Parameters----------n: integerThe number of roller elementfr: float(r/min)Rotational speedd: float(mm)roller element diameterD: float(mm)pitch diameter of bearingalpha: float(°)contact anglefr:：float(r/min)rotational speedReturns-------BPFI: float(Hz)Inner race-way fault frequencyBPFO: float(Hz)Outer race-way fault frequencyBSF: float(Hz)Ball fault frequencyFTF: float(Hz)Cage frequency'''C_bpfi = n*(1/2)*(1+d/D*np.math.cos(alpha))C_bpfo = n*(1/2)*(1-(d/D)*np.math.cos(alpha))C_bsf = D*(1/(2*d))*(1-np.square(d/D*np.math.cos(alpha)))C_ftf = (1/2)*(1-(d/D)*np.math.cos(alpha))if fr!=None:return C_bpfi*fr/60, C_bpfo*fr/60, C_bsf*fr/60, C_ftf*fr/60, fr/60else:return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
D = 40.5d = 9alpha = 0n = 10fr = 24000 # 转速C_bpfi, C_bpfo, C_bsf, C_ftf, fr = bearing_fault_freq_cal(D=D, d=d, alpha=alpha, fr=fr, n=n)print(C_bpfi, C_bpfo, C_bsf, C_ftf, fr)
>>>输出2444.444 1555.555 855.555 155.555 400.0
3.2.2 不同故障类型包络谱
def plt_envelope_spectrum(data, fs, ylabel='Amp(mg)', title='包络谱图', img_save_path=None, vline=None, hline=None, xlim=None):'''fun: 绘制包络谱图param data: 输入数据，1维arrayparam fs: 采样频率param xlim: 图片横坐标xlim，default = Noneparam vline: 图片垂直线，default = None'''from scipy import fftpack#=========做希尔伯特变换=======#xt = dataht = 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)*2fft_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: # 图片横坐标是否设置xlimplt.xlim(0, xlim)plt.tight_layout()plt.show(
file_path = r'E:/03-公开数据集/都灵理工大学变转速轴承数据集/都灵理工大学数据集/VariableSpeedAndLoad/C4A_400_496_1.mat'df = data_read(file_path)acc_arr= df['A1_x']plt_envelope_spectrum(acc_arr, fs=fs, xlim=3000, vline=[C_bsf, C_bsf*2])
选择滚动体故障。可以看到明显的倍频，但是由于实际转速比额定转速低，未能匹配上对应的理论特征频率
file_path = r'E:/03-公开数据集/都灵理工大学变转速轴承数据集/都灵理工大学数据集/VariableSpeedAndLoad/C1A_400_702_2.mat'df = data_read(file_path)acc_arr= df['A1_x']plt_envelope_spectrum(acc_arr, fs=fs, xlim=6000, vline=[C_bpfi, C_bpfi*2])
选择内圈故障。可以看到明显的倍频，但是由于实际转速比额定转速低，未能匹配上对应的理论特征频率。
3.2.3 加速寿命测试数据集包络谱
fs = 102400D = 40.5d = 9alpha = 0n = 10fr = 18000C_bpfi, C_bpfo, C_bsf, C_ftf, fr = bearing_fault_freq_cal(D=D, d=d, alpha=alpha, fr=fr, n=n)print(C_bpfi, C_bpfo, C_bsf, C_ftf, fr)
>>>输出结果1833.333 1166.666 641.666 116.667 300.
file_path = r'E:/03-公开数据集/都灵理工大学变转速轴承数据集/都灵理工大学数据集/EnduranceTest/E4A_007.mat'df = data_read(file_path)acc_arr= df['A1_x']plt_envelope_spectrum(acc_arr, fs=fs, xlim=3000, vline=[C_bsf, C_bsf*2])
file_path = r'E:/03-公开数据集/都灵理工大学变转速轴承数据集/都灵理工大学数据集/EnduranceTest/E4A266.mat'df = data_read(file_path)acc_arr= df['A1_x']plt_envelope_spectrum(acc_arr, fs=fs, xlim=3000, vline=[C_bsf, C_bsf*2]
编辑：李正平
校核：陈凯歌、赵栓栓、董浩杰、曹希铭、赵学功、白亮、陈少华
