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

希尔伯特解调（包络谱）python代码实战及详细讲解，在CWRU数据上验证

1、数据介绍
2、加载CWRU内圈故障数据
3、希尔伯特解调(包络谱)分析
- 3.1希尔伯特黄变换
- 3.2获得包络信号
- 3.3获得包络谱
- 3.4去直流分量
4、计算故障特征频率
- 4.1定义一个轴承故障特征频率计算函数
5、理论故障特征频率与实际故障特征频率验证
6、与fft进行对比分析
7、封装包络谱函数
- 7.1外圈故障数据测试
- 7.2滚动体故障数据测试分析
8、总结

欢迎关注公众号《故障诊断与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.311388140.17055689]

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

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

3.1希尔伯特黄变换

设x(t)x(t)x(t)为一个实时域信号，其Hilbert变换定义为：
h(t)=1π∫−∞+∞x(τ)t−τdτ=x(t)∗1πth(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)+jh(t)=a(t)ejφtz(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)∣=x2(t)+h2(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、计算故障特征频率

内圈故障特征频率：FBPFI=nfr2(1+dDcos⁡α)F_{\mathrm{BPFI}}=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \alpha\right)FBPFI=2nfr(1+Ddcosα)

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

滚动体故障特征频率：FBSF=Dfr2d[1−(dDcos⁡α)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]

nnn: 滚动体个数，frf_{r}fr: 轴转速 ddd: 滚珠(子)直径 DDD: 轴承节径

轴承型号为：6205-2RSL JME SKF 深沟球滚珠轴承
ddd=7.94mm, DDD=39.04mm, α\alphaα=0, nnn=9

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

为了方便，定义了一个轴承故障特征频率计算函数，只需输入参数ddd, DDD, α\alphaα, nnn和frf_{r}fr即可

def bearing_fault_freq_cal(n, d, D, alpha, fr=None):'''基本描述：计算滚动轴承的故障特征频率详细描述：输入4个参数 n, fr, d, D, alphareturn C_bpfi, C_bpfo, C_bsf, C_ftf,  fr内圈    外圈    滚针   保持架  转速Parameters----------n: integerThe number of roller elementfr: float(r/min)Rotational speedd: float(mm)roller element diameterD: float(mm)pitch diameter of bearingalpha: float(°)contact anglefr:：float(r/min)rotational speedReturns-------BPFI: float(Hz)Inner race-way fault frequencyBPFO: float(Hz)Outer race-way fault frequencyBSF: float(Hz)Ball fault frequencyFTF: 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/60else: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维arrayparam fs: 采样频率param xlim: 图片横坐标xlim，default = Noneparam vline: 图片垂直线，default = None'''#----去直流分量----#data = data - np.mean(data)#----做希尔伯特变换----#xt = dataht = 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)*2am = 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: # 图片横坐标是否设置xlimplt.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难以检测故障特征频率，故障特征频率易被高频掩盖