滚动轴承故障检测---基于希尔伯特变换解调(python实现)
一、基本原理
请转向这里查看
二、方法实现
2.1 加载故障数据
下面先以轴承外圈故障数据进行分析。原始数据为mat文件,是matlab文件,定义一个函数进行数据读取
from scipy.io import loadmat
import matplotlib.pyplot as plt
import math
import numpy as np
import scipy
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
def data_acquision(file_path):
"""
fun: 从cwru mat文件读取加速度数据
param file_path: mat文件绝对路径
return accl_data: 加速度数据,array类型
"""
file = loadmat(file_path) # 加载mat数据
file_keys = file.keys()
for key in file_keys:
if 'DE' in key:
files = file[key].ravel()
return files
if __name__ == "__main__":
# step1: 获取测试数据
path = r'D:\360安全浏览器下载\轴承故障检测\数据源\西储大学\CWRU轴承数据\cwru_data\0HP\12k_Drive_End_OR007@6_0_130.mat'
data = data_acquision(path) #获取原始数据
plt.figure(figsize=(12,3))
plt.title('原始数据')
plt.plot(data)
print(data)
2.2 包络谱分析
# step2: 截取数据 符合实际检测要求
xt = data[:2000]
plt.figure(figsize=(12,3))
plt.title('截取数据')
plt.plot(xt)
print(xt)
# step3: 做希尔伯特变换
ht = scipy.fftpack.hilbert(xt)
plt.figure(figsize=(12,3))
plt.title('希尔伯特变换')
plt.plot(ht)
print(ht)
# step4: 获取包络信号
# z(t)=x(t)+jh(t)=a(t)ejφt
# 对z(t)取模,得到其幅值a(t). a(t)即为包络信号,也叫解析信号
at = np.sqrt(ht**2+xt**2) # at = sqrt(xt^2 + ht^2) 实部
plt.figure(figsize=(12,3))
plt.title('包络信号')
plt.plot(at)
print(at)
# step5: 获取包络谱
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
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 / sampling_rate) # 获取fft频率,此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.figure(figsize=(12,3))
plt.title('包络谱')
plt.plot(freq, am)
plt.xlim(0,500)
2.3 特征频率计算
2.3.1 轴承故障的特征频率
内圈故障特征频率:
外圈故障特征频率:
滚动体故障特征频率:
保持架故障特涨频率:
其中, : 滚动体个数;轴转速;: 滚珠(子)直径;: 轴承节径
2.3.2 轴承故障特征频率计算函数
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
滚珠数量
fr: float(r/min)
转速
d: float(mm)
滚珠直径
D: float(mm)
轴承直径
alpha: float(°)
接触角
fr::float(r/min)
转速
Returns
-------
BPFI: float(Hz)
内圈故障频率
BPFO: float(Hz)
外圈故障频率
BSF: float(Hz)
球体故障频率
FTF: float(Hz)
保持架频率
'''
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.3.3 轴承特征频率计算
轴承型号为:6205-2RSL JME SKF 深沟球滚珠轴承
下面计算CWRU在转速1797rpm时的故障特征频率
# step 6
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1797)
print('内圈故障特征频率',bpfi)
print('外圈故障特征频率',bpfo)
print('滚动体故障特征频率',bsf)
print('保持架故障特征频率',ftf)
print('当前转速',fr)
2.4 理论与实际验证
plt.vlines(x=bpfo, ymin=0, ymax=0.2, colors='r') # 一倍频
plt.vlines(x=bpfo*2, ymin=0, ymax=0.2, colors='r') # 二倍频