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

Posted 2023-03-28 故障诊断与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.31138814
  0.17055689]

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

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

3.1希尔伯特黄变换

设$x(t)$为一个实时域信号，其Hilbert变换定义为：
$h(t)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x(\tau )}{t-\tau }\mathrm{d}\tau =x(t)\ast \frac{1}{\pi t}$
则原始信号$x(t)$和它的Hilbert变换信号$h(t)$可以构建一个新的解析信号$z(t)$:
$z(t)=x(t)+jh(t)=a(t)e^{j\phi t}$

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

输出结果：

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

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

对$z(t)$取模，得到其幅值$a(t):$

$a(t)=|z(t)|=\sqrt{x^2(t)+h^2(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_{\mathrm{BPFI}}=\frac{n{f}_r}{2}\left(1+\frac{d}{D}\cos \alpha \right)$

外圈故障特征频率：$F_{\mathrm{BPFO}}=\frac{n{f}_r}{2}\left(1-\frac{d}{D}\cos \alpha \right)$

滚动体故障特征频率：$F_{\mathrm{BSF}}=\frac{D{f}_r}{2d}\left[1-{\left(\frac{d}{D}\cos \alpha \right)}^2\right]$

$n$: 滚动体个数，$f_r$: 轴转速 $d$: 滚珠(子)直径 $D$: 轴承节径

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

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

为了方便，定义了一个轴承故障特征频率计算函数，只需输入参数$d$, $D$, $\alpha$, $n$和$f_r$即可

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))以上是关于信号处理-基于希尔伯特解调（包络谱）的轴承故障诊断实战，通过python代码实现超详细讲解的主要内容，如果未能解决你的问题，请参考以下文章
 
 故障诊断分析基于matlab小波变换外圈轴承故障诊断含Matlab源码 1678期
 故障诊断分析基于matlab小波包能量分析轴承故障诊断含Matlab源码 1620期
 轴承故障分解基于matlab ITD轴承故障信号分解含Matlab源码 1871期
 滚动轴承特征频率计算推导过程
 关于包络概念的理解
 数字信号基于matlab轴承故障分析（时频图+功率谱）含Matlab源码 002期