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

白银时代_

已于 2022-07-26 17:00:22 修改

阅读量8.5k

点赞数 3

分类专栏：深入浅出学习凯斯西储大学(CWRU)轴承数据集文章标签： python 人工智能开发语言

于 2022-07-20 22:38:15 首次发布

本文链接：https://blog.csdn.net/pengchen_wuhan/article/details/125902107

版权

深入浅出学习凯斯西储大学(CWRU)轴承数据集专栏收录该内容

5 篇文章

订阅专栏

本文介绍了CWRU轴承数据集，包括试验设施和数据简介，详细探讨了轴承数据的初步探索，特别是125DE信号的时域和频谱特性，揭示了机械和电磁部分的中高频线谱。通过对特征频率的分析，为滚动轴承的故障诊断提供基础。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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{aligned} &B P F O=\frac{n f_{r}}{2}\left(1-\frac{d}{D} \cos \phi\right) \end{aligned}\tag{4-1}$

滚动体通过内圈频率：
$\begin{aligned} &B P F I=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \phi\right) \end{aligned}\tag{4-2}$

保持架旋转频率：
$\begin{aligned} &F T F=\frac{f_{r}}{2}\left(1-\frac{d}{D} \cos \phi\right) \end{aligned}\tag{4-3}$

滚动体自传频率：
$\begin{aligned} &B S F=\frac{D f_{r}}{2 d}\left(1-\left[\frac{d}{D} \cos \phi\right]^{2}\right) \end{aligned}\tag{4-4}$

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

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