python实战故障诊断之CWRU数据集（一）：数据集初识

1. 概述

  在完成了振动信号处理的基础篇与高级篇后，不少读者建议采用真实数据对各类振动信号算法进行剖析。经过几个月的搜集与积累，笔者决定从公开的凯斯西储大学(CWRU)轴承数据入手，以新南威尔士大学Robert B. Randall教授数十年多篇论文为切入点，继续与大家探索隐藏在振动信号中的奥秘。

  本系列第一篇文章，是对CWRU数据集进行一个简介，同时为大家提供了相关数据的下载地址和数据初识。

2. CWRU数据集简介

  基于振动信号的滚动轴承故障诊断经过多年的发展，已经成为一个较为成熟的领域，但目前相关研究人员仍在不断开发新的故障诊断算法，提升相关技术。而对于新的故障诊断算法，若没有公认的基准，很难合理有效地评估其性能。在过去的十年中，来自凯斯西储大学(CWRU)轴承数据中心的数据已经成为测试这些算法的参考标准。数据下载地址

2.1. 试验设施简介

  凯斯西储大学(CWRU)轴承数据中心试验台的基本布局如图2.1所示。它由2马力伺服电机直驱，驱动轴尾部安装了测功仪，用于模拟负载。同时在驱动轴上还安装了扭矩测量装置（扭矩传感器与编码器），可将采集到的扭矩信息反馈给控制系统，从而调整测功仪的功率，进而控制负载。

图2.1 实验设施示意图

2.2. 试验数据简介

  在测试中，使用电火花加工(EDM)，在电机的驱动端和风扇端轴承(SKF深沟球轴承:6205-2RS JEM和6203-2RS JEM)上植入直径为0.007 - 0.028英寸(0.18 - 0.71 mm)的故障，故障分别被设置在滚动体、内圈和外圈上。每种故障轴承均被以同样的状态安装在试验台上，然后在0到3马力的电机负载下(大约电机转速为1797到1720 rpm)以恒速运行。

  在每次测试中，测量驱动端轴承壳体(DE)垂直方向上的加速度，在一些测试中也测量风扇端轴承壳体(FE)和电机支撑底板(BA)垂直方向上的加速度。某些测试使用的采样率为12 kHz，其他测试使用的采样率为48 kHz。通过测试共计获得161组数据，各数据对应的具体工况如下图所示：

图2.2 数据集简介

3. 轴承数据初步探索

  提取125DE数据，结果如图3.1（相关代码如下），信号时域部分较为平稳（信号峰度为3.04），但其频谱图中存在较多特征线谱，特别是在11-14 kHz波段。

图3.1 125DE信号

from scipy.io import loadmat
import numpy as np
from scipy import signal, fftpack, stats
import matplotlib.pyplot as plt

# 从.mat文件中解析数据
def get_data(source, fs):
    for i, key in enumerate(source.keys()):
        if i == 3:
            data_DE = source[key].flatten()
        elif i == 4:
            data_FE = source[key].flatten()
        elif i == 5:
            data_Speed = source[key].flatten()
    t = np.arange(len(data_DE))/fs
    return t, data_DE, data_FE, data_Speed

# 导入数据
source = loadmat('.//data//48k Drive End Bearing Fault Data/125.mat')
fs = 48000
t, data_DE, data_FE, data_Speed = get_data(source, fs)

# 计算信号偏度
sig_kurtosis = stats.kurtosis(data_DE, fisher=False)

# 计算FFT结果
sig_fre, sig_mag = signal.welch(data_DE, fs=fs, nperseg=fs, noverlap=2/3, nfft=fs, scaling = 'spectrum')
sig_mag = 120 + 20 * np.log10(np.sqrt(2 * sig_mag))
# 展示时域信号及频域信号
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t, data_DE, linewidth=0.4)
plt.xlim(0, 10)
plt.ylim(-0.6, 0.6)
plt.xlabel('time/s')
plt.ylabel('magnitude/(m/(s^2))')
plt.subplot(1, 2, 2)
plt.plot(sig_fre, sig_mag, linewidth=0.4)
plt.xlim(0, 20000)
plt.ylim(0, 90)
plt.xlabel('frequence/Hz')
plt.ylabel('magnitude/dB')
plt.show()

  经过前人研究，得出结论，信号中的中高频离散线谱成分主要来源于机械部分以及电磁部分。

  首先介绍机械部分中高频线谱，从图3.2中（相关代码如下）可以发现，其在11-14 kHz中存在较多以轴频为间隔的线谱，其产生的原因可能是机械松动，机械松动会造成非常尖锐的冲击，其在频谱中表现为机械固有频率为载波频率，周围存在大量的轴谐频边带谱。

图3.2 125DE信号（11-14kHz）

# 展示11-14 kHz
plt.figure(figsize=(10, 4))
plt.plot(sig_fre, sig_mag, linewidth=0.4)
plt.xlim(11000, 14000)
plt.ylim(0, 90)
plt.xlabel('frequence/Hz')
plt.ylabel('magnitude/dB')
plt.show()

  紧接着，我们介绍电磁部分中高频线谱，从图3.3中（相关代码如下）可以发现，其在4-9 kHz中存在较多以近似轴频为间隔的线谱，其可能来源于电磁测功仪，因其工作原理，其在负载时存在速度差，会造成控制频率（即轴频）滑移。

图3.3 125DE信号（4-9kHz）

# 展示4-9 kHz
plt.figure(figsize=(10, 4))
plt.plot(sig_fre, sig_mag, linewidth=0.4)
plt.xlim(4000, 9000)
plt.ylim(0, 90)
plt.xlabel('frequence/Hz')
plt.ylabel('magnitude/dB')
plt.show()

4. 轴承的故障特征频率探索

  为了顺利完成故障诊断的任务，我们还需深入理解轴承的故障特征，滚动轴承的个故障诊断频率具体推导过程可借鉴滚动轴承的运动学(特征频率与阶次)，我们这里将结果直接列出。

  滚动体通过外圈频率：
$$\begin{array}{rl}& BPFO=\frac{n{f}_r}{2}\left(1-\frac{d}{D}\cos \varphi \right)\end{array}\text{\hspace{1em}(4-1)}$$

  滚动体通过内圈频率：
$$\begin{array}{rl}& BPFI=\frac{n{f}_r}{2}\left(1+\frac{d}{D}\cos \varphi \right)\end{array}\text{\hspace{1em}(4-2)}$$

  保持架旋转频率：
$$\begin{array}{rl}& FTF=\frac{f_r}{2}\left(1-\frac{d}{D}\cos \varphi \right)\end{array}\text{\hspace{1em}(4-3)}$$

  滚动体自传频率：
$$\begin{array}{rl}& BSF=\frac{D{f}_r}{2d}\left(1-{\left[\frac{d}{D}\cos \varphi \right]}^2\right)\end{array}\text{\hspace{1em}(4-4)}$$

  其中$f_r$为驱动轴旋转频率，$n$为滚动体个数，$d$为滚动体直径，$D$为滚动体中心所在圆直径，$\varphi$为滚动体接触角。我们将各频率进行代码化，具体代码如下：

import numpy as np
def get_fre(n, fr, D, d, phi):
    BPFO = (n*fr)/2*(1-d/D*np.cos(phi))
    BPFI = (n*fr)/2*(1+d/D*np.cos(phi))
    FTF = fr/2*(1-d/D*np.cos(phi))
    BSF = (D*fr)/(2*d)*(1-(d/D*np.cos(phi))**2)    
    return BPFO, BPFI, FTF, BSF

  下面我们对CWRU数据集中所用的滚动轴承的特征频率进行分析，驱动端和风扇端用到的轴承具体参数（单位为：inches）如下表所示：

表4.1 驱动端轴承参数(n=9)

内圈直径	外圈直径	滚动体直径	节圆直径
0.9843	2.0472	0.3126	1.537

表4.2 风扇端轴承参数(n=8)

内圈直径	外圈直径	滚动体直径	节圆直径
0.6693	1.5748	0.2656	1.122

  深沟球轴承接触角为$0^{。}$,带入各计算参数可得：

表4.3 各轴承特征频率阶次

轴承名称	BPFO	BPFI	FTF	BSF
驱动端轴承	3.58	5.42	0.40	2.36
风扇端轴承	3.05	4.95	0.38	1.99