读研之掉进故障检测（一）—CWRU轴承数据集的三种已建立的诊断技术
`提示：本人经常断更，
文章目录
前言
关于怎么开始这方面的研究，其实本人一直在纠结。是该从时下的提出的模型入手，还是先从论文中使用的数据集入手？
因为研一开学的很晚，所以才刚刚开始这方面学习一个星期。
`提示：本章将针对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)
关于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
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)
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)
总结
仅做记录与参考。因为我的在下载数据的时候，已经不能从官网（https://csegroups.case.edu/bearingdatacenter/pages/welcome-case-western-reserve-university-bearing-data-center-website）下载数据了，因此我也不确定自己的结果是否正确。这里只提供论文中三种方法可能实现的方法。同时也没有按照论文中的相关参数进行设置，仅仅是为了方便对论文以及数据的理解做出的尝试。
6474
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
