信号处理-基于希尔伯特解调(包络谱)的轴承故障诊断实战,通过python代码实现超详细讲解 轴承故障诊断实战-基于希尔伯特解调(包络谱)的轴承故障诊断实战,通过python代码实现超详细讲解


欢迎关注公众号《故障诊断与python学习》

代码位置:https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing
参考资料:
机械故障诊断及典型案例解析(第2版,时献江)
会议论文:Bearing Intelligent Fault Diagnosis Under Complex Working Condition Based on SK-ES-CNN,2021 Global Reliability and Prognostics and Health Management (PHM-Nanjing)

1、数据介绍


包括内圈、外圈、滚动体和正常数据,分别为一个。
1730_12k_0.007-InnerRace:内圈故障
1730_12k_0.007-OuterRace3:外圈故障
1730_12k_0.014-Ball:滚动体故障
1730_48k_Normal:正常数据

对CWRU轴承数据集不了解的同学见这里:
CWRU数据集介绍 第1期
CWRU数据集介绍 第2期
CWRU数据集介绍 第3期
CWRU数据集介绍 第3期

2、加载CWRU内圈故障数据

下面先以轴承内圈故障数据进行分析。原始数据为mat文件,是matlab文件,定义一个函数进行数据读取

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

data_acquision(FilePath)输入参数FilePath,输出一个1维array数据。下面进行演示

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-InnerRace.mat'
xt = data_acquision(file_path)
plt.figure(figsize=(12,3))
plt.plot(xt)
print(xt)

输出结果:

[ 0.22269856 0.09323776 -0.14651649 ... -0.36125573 0.31138814
 0.17055689]

3、希尔伯特解调(包络谱)分析

希尔伯特解调法,亦叫包络谱分析。

3.1希尔伯特黄变换

x ( t ) x(t) x(t)为一个实时域信号,其Hilbert变换定义为:
h ( t ) = 1 π ∫ − ∞ + ∞ x ( τ ) t − τ d τ = x ( t ) ∗ 1 π t h(t)=\frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{x(\tau)}{t-\tau} \mathrm{d} \tau=x(t) * \frac{1}{\pi t} h(t)=π1+tτx(τ)dτ=x(t)πt1
则原始信号 x ( t ) x(t) x(t)和它的Hilbert变换信号 h ( t ) h(t) h(t)可以构建一个新的解析信号 z ( t ) z(t) z(t):
z ( t ) = x ( t ) + j h ( t ) = a ( t ) e j φ t z(t)=x(t)+j h(t)=a(t) e^{j \varphi t} z(t)=x(t)+jh(t)=a(t)ejφt

# step1: 做希尔伯特变换
ht = fftpack.hilbert(xt)
print(ht)

输出结果:

[-0.02520403 -0.28707983 -0.00610516 ... 0.1100125 0.22821944
 -0.11203138]

此时输出的 h ( t ) h(t) h(t)是解析信号 a ( t ) a(t) a(t)的虚部系数

z ( t ) z(t) z(t)取模,得到其幅值 a ( t ) : a(t): a(t):

a ( t ) = ∣ z ( t ) ∣ = x 2 ( t ) + h 2 ( t ) {a(t)=|z(t)|=\sqrt{x^{2}(t)+h^{2}(t)}} a(t)=z(t)=x2(t)+h2(t)

注: a ( t ) a(t) a(t)即为包络信号,也叫解析信号

3.2获得包络信号

t = np.sqrt(ht**2+xt**2) # at = sqrt(xt^2 + ht^2)

接下来对包络信号做fft即为包络谱

3.3获得包络谱

sampling_rate = 12000
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.plot(freq, am)


观察发现:
(1)在频率为0Hz的地方幅值比较大
(2)在低频部分貌似看到一倍频和二倍频

3.4去直流分量

在0Hz的幅值比较高,使得其它频率幅值较低,不便观察。这种现象叫直流分量,去直流分量方法,y = y-mean(y)

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.plot(freq, am)

输出结果:

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.plot(freq, am)
plt.xlim(0,500)

输出结果:

4、计算故障特征频率

内圈故障特征频率: F B P F I = n f r 2 ( 1 + d D cos ⁡ α ) F_{\mathrm{BPFI}}=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \alpha\right) FBPFI=2nfr(1+Ddcosα)

外圈故障特征频率: F B P F O = n f r 2 ( 1 − d D cos ⁡ α ) F_{\mathrm{BPFO}}=\frac{n f_{r}}{2}\left(1-\frac{d}{D} \cos \alpha\right) FBPFO=2nfr(1Ddcosα)

滚动体故障特征频率: F B S F = D f r 2 d [ 1 − ( d D cos ⁡ α ) 2 ] F_{\mathrm{BSF}}=\frac{D f_{r}}{2 d}\left[1-\left(\frac{d}{D} \cos \alpha\right)^{2}\right] FBSF=2dDfr[1(Ddcosα)2]

n n n: 滚动体个数, f r f_{r} fr: 轴转速 d d d: 滚珠(子)直径 D D D: 轴承节径

轴承型号为:6205-2RSL JME SKF 深沟球滚珠轴承
d d d=7.94mm, D D D=39.04mm, α \alpha α=0, n n n=9

4.1定义一个轴承故障特征频率计算函数

为了方便,定义了一个轴承故障特征频率计算函数,只需输入参数 d d d, D D D, α \alpha α, n n n f r f_{r} fr即可

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

下面计算CWRU在转速1730rpm时的故障特征频率

bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('内圈故障特征频率',bpfi)
print('外圈故障特征频率',bpfo)
print('滚动体故障特征频率',bsf)
print(ftf)
print(fr)

5、理论故障特征频率与实际故障特征频率验证

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.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r') # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r') # 二倍频

输出结果:

红色为理论内圈故障特征频率,蓝线为实际故障特征频率。虽然没有完全重合,但这个是在允许范围内的。因此在实际情况中包络谱出现该情况,可以该轴承出现了内圈故障。

6、与fft进行对比分析

那如果不希尔伯特变换解调,直接fft,能够看到故障特征频率吗?下面进行对比分析

sampling_rate = 12000
am = np.fft.fft(xt) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值(此时的绝对值很大)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate) # 获取fft频率,此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.plot(freq, am)

输出结果:

可见原始内圈故障信号的fft高频较多

plt.plot(freq, am)
plt.xlim(0, 500)
plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r') # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r') # 二倍频


可见直接fft的话,故障特征频率幅值较低,被其它频率幅值掩盖了。反过来,希尔伯特解调可以更加方便观察故障特征频率低频。

7、封装包络谱函数

为了更加方便使用包络谱,这里封装了一个包络谱函数plt_envelope_spectrum

def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
 '''
 fun: 绘制包络谱图
 param data: 输入数据,1维array
 param fs: 采样频率
 param xlim: 图片横坐标xlim,default = None
 param vline: 图片垂直线,default = None
 '''
 #----去直流分量----#
 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)') # 纵坐标标签

7.1外圈故障数据测试

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-OuterRace3.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)

输出结果:

可见实际外圈故障特征频率比较明显

7.2滚动体故障数据测试分析

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.014-Ball.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)

可见实际滚动体故障特征频率不明显

8、总结

(1)包络谱能够检测出内圈、外圈故障,滚动体比较困难
(2)直接使用fft难以检测故障特征频率,故障特征频率易被高频掩盖

作者:故障诊断与python学习原文地址:https://blog.csdn.net/m0_47410750/article/details/126004829

%s 个评论

要回复文章请先登录注册