信号处理-基于希尔伯特解调（包络谱）的轴承故障诊断实战，通过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
1
2
3
4
5
6
7
8
9
10
11

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)
1
2
3
4
5

输出结果：

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

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)
1
2
3

输出结果：

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

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

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

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

注：$a(t)$即为包络信号，也叫解析信号

3.2获得包络谱

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
2
3
4
5
6
7
8

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

3.3去直流分量

在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)
1
2
3
4
5
6
7
8
9

输出结果：

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)
1
2
3
4
5
6
7
8
9
10
11

输出结果：

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))*(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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

下面计算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)
1
2
3
4
5
6

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')  # 二倍频
1
2
3
4
5
6
7
8
9
10
11
12
13

输出结果：

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

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)
1
2
3
4
5
6
7
8

输出结果：

可见原始内圈故障信号的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')  # 二倍频
1
2
3
4

可见直接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)')   # 纵坐标标签
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

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)
1
2
3

输出结果：

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

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)
1
2
3

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

8、总结

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