读研之掉进故障检测（一）—CWRU轴承数据集的三种已建立的诊断技术_西储大学轴承数据集进行快速傅里叶变换

读研之掉进故障检测（一）—CWRU轴承数据集的三种已建立的诊断技术

---

`提示：本人经常断更，

---

前言

关于怎么开始这方面的研究，其实本人一直在纠结。是该从时下的提出的模型入手，还是先从论文中使用的数据集入手？
因为研一开学的很晚，所以才刚刚开始这方面学习一个星期。
1
2

---

`提示：本章将针对Wade A. Smith在2015年发表的关于CWRU轴承论文进行了一些讨论。主要是文中提到的三种方案的实现。
具体的论文解读，可以看该博客。介绍的详细的

一、凯斯西储大学（CWRU）轴承数据

该数据应用十分广泛，凯斯西储大学（CWRU）轴承数据中心提供的数据集已成为轴承诊断领域的标准参考，2004年至2015年初发表在《机械系统和信号处理》上的41篇使用CWRU数据的论文。

论文通过将三种已建立的诊断技术应用于整个CWRU数据集来提供这样的基准。所有方法都使用平方包络频谱（即平方包络的频谱）作为最终诊断工具，但是在获取包络信号之前使用了不同的预处理步骤

二、三种方案

1、基于包络分析的方法

为什么使用包络分析法？
包络分析是一种广泛应用于故障检测领域的方法。它通过提取信号的振幅包络，在频域表示了信号的能量分布情况。包络分析对于检测低频振动特征和轴承等机械部件的故障非常有效，因为这些故障通常表现为周期性脉冲或冲击。
关于包络分析代码，需要感谢该文章。博主比较详细的介绍怎么使用包络分析进行故障检测。

那么，由于源代码是.ipynb。需要.py代码，可以如下：

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fftpack, stats
from matplotlib import rcParams
from BFF_cal import bearing_fault_freq_cal #可以直接按原博客的直接写成函数放在代码中
config = {
    "font.family": 'serif', # 衬线字体
    "font.size": 10, # 相当于小四大小
    "font.serif": ['SimSun'], # 宋体
    "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False # 处理负号，即-号
}
rcParams.update(config)

def data_acquision(FilePath):
    """
    fun: 从cwru mat文件读取加速度数据
    param file_path: mat文件绝对路径
    return accl_data: 加速度数据，array类型
    """
    data = scio.loadmat(file_path)  # 加载mat数据
    data_key_list = list(data.keys())  # mat文件为字典类型，获取字典所有的键并转换为list类型
    accl_key = data_key_list[3]  # 获取'X108_DE_time'
    accl_data = data[accl_key].flatten()  # 获取'X108_DE_time'所对应的值，即为振动加速度信号,并将二维数组展成一维数组
    return accl_data

#包络谱分析
def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
    '''
    fun: 绘制包络谱图
    param data: 输入数据，1维array
    param fs: 采样频率
    '''
    #----去直流分量----#
    data = data - np.mean(data)
    #----做希尔伯特变换----#
    xt = data
    ht = fftpack.hilbert(xt)
    at = np.sqrt(xt**2+ht**2)   # 获得解析信号at = sqrt(xt^2 + ht^2)
    am = np.fft.fft(at)         # 对解析信号at做fft变换获得幅值
    am = np.abs(am)             # 对幅值求绝对值（此时的绝对值很大）
    am = am/len(am)*2
    am = am[0: int(len(am)/2)]  # 取正频率幅值
    freq = np.fft.fftfreq(len(at), d=1 / fs)  # 获取fft频率，此时包括正频率和负频率
    freq = freq[0:int(len(freq)/2)]  # 获取正频率
    plt.plot(freq, am)
    if vline:  # 是否绘制垂直线
        plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r')  # 高度y 0-0.2，颜色红色
    if xlim: # 图片横坐标是否设置xlim
        plt.xlim(0, xlim)
    plt.xlabel('freq(Hz)')    # 横坐标标签
    plt.ylabel('amp(m/s2)')   # 纵坐标标签
    plt.show()



file_path = r'F:\pycharm_dataset\cwru\12k Fan End Bearing Fault Data\270.mat'
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfi)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

关于BFF_cal.py文件代码如下

import numpy as np

def bearing_fault_freq_cal(n, d, D, alpha, fr=None):
    '''
    基本描述：
        计算滚动轴承的故障特征频率
    详细描述：
        输入4个参数 n, fr, d, D, alpha
    return C_bpfi, C_bpfo, C_bsf, C_ftf,  fr
           内圈    外圈    滚针   保持架  转速

    Parameters
    ----------
    n: integer
        The number of roller element
    fr: float(r/min)
        Rotational speed
    d: float(mm)
        roller element diameter
    D: float(mm)
        pitch diameter of bearing
    alpha: float(°)
        contact angle
    fr:：float(r/min)
        rotational speed
    Returns
    -------
    BPFI: float(Hz)
        Inner race-way fault frequency
    BPFO: float(Hz)
        Outer race-way fault frequency
    BSF: float(Hz)
        Ball fault frequency
    FTF: 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/60
    else:
        return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

2、倒谱预白化

该方案的描述如下：
1、倒谱预白化，将所有频率分量设置为相同的幅度
2、全带宽信号的包络分析（平方包络谱）

为什么使用倒谱预白化？
倒谱白化有助于使信号中的频率成分具有相同的幅度，从而更容易检测到与故障相关的特征。通过消除由不同频率幅度引起的频谱偏差，该方法可以增强原始信号中可能被掩盖的故障特征

代码如下（示例）：

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fftpack, stats
from matplotlib import rcParams
from BFF_cal import bearing_fault_freq_cal

config = {
    "font.family": 'serif',  # 衬线字体
    "font.size": 10,  # 相当于小四大小
    "font.serif": ['SimSun'],  # 宋体
    "mathtext.fontset": 'stix',  # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False  # 处理负号，即-号
}
rcParams.update(config)


def data_acquisition(file_path):
    """
    fun: 从cwru mat文件读取加速度数据
    param file_path: mat文件绝对路径
    return accl_data: 加速度数据，array类型
    """
    data = scio.loadmat(file_path)  # 加载mat数据
    data_key_list = list(data.keys())  # mat文件为字典类型，获取字典所有的键并转换为list类型
    accl_key = data_key_list[3]  # 获取'X108_DE_time'
    accl_data = data[accl_key].flatten()  # 获取'X108_DE_time'所对应的值，即为振动加速度信号，并将二维数组展成一维数组
    return accl_data


# 包络谱分析
def plt_envelope_spectrum(data, fs, xlim=None, vline=None):
    '''
    fun: 绘制包络谱图
    param data: 输入数据，1维array
    param fs: 采样频率
    '''
    # ----去直流分量----#
    data = data - np.mean(data)

    '''
    Cepstrum预白化：首先，将输入数据进行傅里叶变换，然后取其幅度的对数。这相当于将信号从时域转换到倒谱（cepstrum）域。
    通过执行这一步骤，我们可以将信号的频谱信息转移到倒谱域，其中低频成分位于倒谱的高频部分，高频成分位于倒谱的低频部分。
    设置前几个倒谱系数为零：在进行倒谱操作后，通常会出现一个主要峰值，该峰值对应于信号的主要周期或频率。
    为了消除主要峰值对信号的影响，代码将前几个倒谱系数设置为零。这样做可以减小主要峰值的幅度，从而减小对信号频谱的影响
    '''
    # ----Cepstrum prewhitening----#
    cepstrum = np.fft.ifft(np.log(np.abs(np.fft.fft(data))))
    cepstrum[:10] = 0  

    # ----做希尔伯特变换----#
    xt = np.real(np.fft.fft(cepstrum))  
    ht = fftpack.hilbert(xt)
    at = np.sqrt(xt ** 2 + ht ** 2)  # 获得解析信号at = sqrt(xt^2 + ht^2)
    am = np.fft.fft(at)  # 对解析信号at做fft变换获得幅值
    am = np.abs(am)  # 对幅值求绝对值（此时的绝对值很大）
    am = am / len(am) * 2
    am = am[0: int(len(am) / 2)]  # 取正频率幅值
    freq = np.fft.fftfreq(len(at), d=1 / fs)  # 获取fft频率，此时包括正频率和负频率
    freq = freq[0:int(len(freq) / 2)]  # 获取正频率
    plt.plot(freq, am)
    if vline:  # 是否绘制垂直线
        plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r')  # 高度y 0-0.2，颜色红色
    if xlim:  # 图片横坐标是否设置xlim
        plt.xlim(0, xlim)
    plt.xlabel('freq(Hz)')  # 横坐标标签
    plt.ylabel('amp(m/s2)')  # 纵坐标标签
    plt.show()


file_path = r'F:\pycharm_dataset\cwru\12k Fan End Bearing Fault Data\270.mat'
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
data = data_acquisition(file_path)
plt_envelope_spectrum(data=data, fs=12000, xlim=300, vline=bpfi)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

3、基准方法（benchmark method）*（关于谱峰度带通滤波存在一定的错误，待修改24.1.07）

该方案的描述如下：
1、离散/随机分离（DRS）以去除确定性（离散频率）分量
2、谱峰度确定最冲带，其次进行带通滤波
3、带通滤波信号的包络分析（平方包络谱）

为什么？
该方法通过将信号分解为离散分量和随机分量，将信号中的确定性（离散频率）部分与随机噪声部分进行分离。这种方法对于去除背景噪声、基线漂移等干扰成分非常有用，使得后续的故障特征提取更加准确。

代码如下（示例）：

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fftpack, stats
from matplotlib import rcParams
from BFF_cal import bearing_fault_freq_cal
'''
添加了drs函数实现离散/随机分离（Discrete/Random Separation）操作，用于去除确定性（离散频率）成分。该函数将输入数据根据设定的阈值将其划分为随机分量和离散分量。

添加了spectral_kurtosis函数实现谱峭度（Spectral Kurtosis）和带通滤波操作。该函数计算输入数据的峭度和谱峭度，然后根据谱峭度的值确定最具冲击性的频带，并进行带通滤波。

在主程序中按照步骤调用相应的函数，首先进行离散/随机分离操作，然后进行谱峭度和带通滤波操作。最后，对滤波后的数据进行包络分析。
'''
config = {
    "font.family": 'serif',  # 衬线字体
    "font.size": 10,  # 相当于小四大小
    "font.serif": ['SimSun'],  # 宋体
    "mathtext.fontset": 'stix',  # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False  # 处理负号，即-号
}
rcParams.update(config)


def data_acquisition(file_path):
    """
    fun: 从cwru mat文件读取加速度数据
    param file_path: mat文件绝对路径
    return accl_data: 加速度数据，array类型
    """
    data = scio.loadmat(file_path)  # 加载mat数据
    data_key_list = list(data.keys())  # mat文件为字典类型，获取字典所有的键并转换为list类型
    accl_key = data_key_list[3]  # 获取'X108_DE_time'
    accl_data = data[accl_key].flatten()  # 获取'X108_DE_time'所对应的值，即为振动加速度信号，并将二维数组展成一维数组
    return accl_data


#离散/随机分离（DRS）
def drs(data):
    """
    fun: 离散/随机分离（DRS）以去除确定性分量
    param data: 输入数据，1维array
    return random_comp: 随机分量，1维array
           discrete_comp: 离散分量，1维array
    """
    mean_val = np.mean(data)
    std_val = np.std(data)
    threshold = 3 * std_val
    random_comp = np.where(np.abs(data - mean_val) <= threshold, data, 0)
    discrete_comp = np.where(np.abs(data - mean_val) > threshold, data, 0)
    return random_comp, discrete_comp


#谱峰度与带通滤波
def spectral_kurtosis(data, fs):
    """
    fun: Spectral kurtosis to determine the most impulsive band, followed by bandpass filtering
    param data: 输入数据，1维array
    param fs: 采样频率
    return filtered_data: 经过带通滤波后的数据，1维array
    """
    skewness = stats.skew(data)
    kurtosis = stats.kurtosis(data)
    freq, psd = signal.welch(data, fs=fs, nperseg=len(data))
    spectral_kurt = np.sum(psd**4 * freq) / np.sum(psd**2 * freq)**2
    if spectral_kurt < 100:
        lower_cutoff = 0.5 * freq[np.argmax(psd)]
        upper_cutoff = min(8 * freq[np.argmax(psd)], fs / 2 - 1)
        b, a = signal.butter(4, [lower_cutoff, upper_cutoff], btype='band', fs=fs)
        filtered_data = signal.filtfilt(b, a, data)
        return filtered_data
    else:
        return data


# 包络谱分析
def plt_envelope_spectrum(data, fs, xlim=None, vline=None):
    '''
    fun: 绘制包络谱图
    param data: 输入数据，1维array
    param fs: 采样频率
    '''
    # ----去直流分量----#
    data = data - np.mean(data)

    # ----做希尔伯特变换----#
    xt = data
    ht = fftpack.hilbert(xt)
    at = np.sqrt(xt ** 2 + ht ** 2)  # 获得解析信号at = sqrt(xt^2 + ht^2)
    am = np.fft.fft(at)  # 对解析信号at做fft变换获得幅值
    am = np.abs(am)  # 对幅值求绝对值（此时的绝对值很大）
    am = am / len(am) * 2
    am = am[0: int(len(am) / 2)]  # 取正频率幅值
    freq = np.fft.fftfreq(len(at), d=1 / fs)  # 获取fft频率，此时包括正频率和负频率
    freq = freq[0:int(len(freq) / 2)]  # 获取正频率
    plt.plot(freq, am)
    if vline:  # 是否绘制垂直线
        plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r')  # 高度y 0-0.2，颜色红色
    if xlim:  # 图片横坐标是否设置xlim
        plt.xlim(0, xlim)
    plt.xlabel('freq(Hz)')  # 横坐标标签
    plt.ylabel('amp(m/s2)')  # 纵坐标标签
    plt.show()


file_path = r'F:\pycharm_dataset\cwru\12k Fan End Bearing Fault Data\270.mat'
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
data = data_acquisition(file_path)

# Step 1: Discrete/Random Separation (DRS)
random_comp, discrete_comp = drs(data)

# Step 2: Spectral Kurtosis and Bandpass Filtering
filtered_data = spectral_kurtosis(discrete_comp, fs=12000)

# Step 3: Envelope Analysis
plt_envelope_spectrum(data=filtered_data, fs=12000, xlim=300, vline=bpfi)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

---

总结

仅做记录与参考。因为我的在下载数据的时候，已经不能从官网（https://csegroups.case.edu/bearingdatacenter/pages/welcome-case-western-reserve-university-bearing-data-center-website）下载数据了，因此我也不确定自己的结果是否正确。这里只提供论文中三种方法可能实现的方法。同时也没有按照论文中的相关参数进行设置，仅仅是为了方便对论文以及数据的理解做出的尝试。
1