python实战故障诊断之CWRU数据集（二）：异常数据剔除及包络解调初步探索

作者：繁依Fanyi0 | 2024-03-07 01:01:11

踩

包络解调

文章目录

1. 概述
2. 异常数据探索
3. 正常信号的平方包络解调分析

1. 概述

在完成了CWRU数据的初步探索后，我们需要注意到，该数据库中存在较多异常数据，具体异常表现形式包括：电噪声干扰、驱动端与风扇端传感器信号混淆以及分段采集信号整合。本章节将首先带领大家观察并剔除这些异常信号，此外对于正常信号，我们采用平方解包络的方法对信号进行初步探索。

2. 异常数据探索

2.1. 电噪声干扰

通过对时域信号进行观察，可发现数据集中177.mat文件存在电噪声，具体表现形式如图2.1所示,具体代码如下。
在这里插入图片描述

图2.1 177DE时域信号


from scipy.io import loadmat
import numpy as np
from scipy import signal, fftpack, stats
import matplotlib.pyplot as plt

# 从.mat文件中解析数据
def get_data(source, fs):
    data_DE = []
    data_FE = []
    data_Speed = []
    
    for i, key in enumerate(source.keys()):
        if i == 3:
            data_DE = source[key].flatten()
        elif i == 4:
            data_FE = source[key].flatten()
        elif i == 5:
            data_Speed = source[key].flatten()
    t = np.arange(len(data_DE))/fs
    return t, data_DE, data_FE, data_Speed

source = loadmat('.//data//48k Drive End Bearing Fault Data/177.mat')
fs = 48000
t, data_DE, data_FE, data_Speed = get_data(source, fs)


# 展示时域信号及频域信号
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t, data_DE, linewidth=0.4)
plt.xlim((0.4, 0.6))
plt.xlabel('time/s')
plt.ylabel('magnitude/(m/(s^2))')
plt.title('data_DE')
plt.subplot(1, 2, 2)
plt.plot(t, data_FE, linewidth=0.4)
plt.xlim((0.4, 0.6))
plt.xlabel('time/s')
plt.ylabel('magnitude/(m/(s^2))')
plt.title('data_FE')
plt.show()

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
43

2.2. 驱动端与风扇端传感器信号混淆

通过对时域信号进行观察，可发现数据集中189、201、213、226、238文件，均存在驱动端与风扇端传感器信号混淆的现象，具体表现为两种信号为同一种信号，仅存在1.0154的比值，如图2.2所示，具体代码如下。

在这里插入图片描述

图2.2 177DE时域信号

paths = [189, 201, 213, 226, 238]
for _path in paths:
    file_path = './/data//48k Drive End Bearing Fault Data/' + str(_path) +  '.mat'
    source = loadmat(file_path)
    fs = 48000
    t, data_DE, data_FE, data_Speed = get_data(source, fs)


    # 展示时域信号及频域信号
    plt.figure(figsize=(10, 4))
    plt.subplot(1, 2, 1)
    plt.plot(t, data_DE, linewidth=0.4)
    plt.xlim((0.4, 0.6))
    plt.xlabel('time/s')
    plt.ylabel('magnitude/(m/(s^2))')
    plt.title('data_DE')
    plt.subplot(1, 2, 2)
    plt.plot(t, data_FE, linewidth=0.4)
    plt.xlim((0.4, 0.6))
    plt.xlabel('time/s')
    plt.ylabel('magnitude/(m/(s^2))')
    plt.title('data_FE')
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

2.3. 分段采集信号整合

通过对时域信号进行观察，可发现数据集中191, 228, 229文件，均存在分段采集信号整合的现象，具体表现为信号中间由明显间断，如图2.3所示，具体代码如下。

在这里插入图片描述

图2.3 191DE时域信号

paths_1 = [191, 228, 229]
for _path in paths_1:
    file_path = './/data//48k Drive End Bearing Fault Data/' + str(_path) +  '.mat'
    source = loadmat(file_path)
    fs = 48000
    t, data_DE, data_FE, data_Speed = get_data(source, fs)

    # 展示时域信号及频域信号
    plt.figure(figsize=(10, 4))
    plt.subplot(1, 2, 1)
    plt.plot(t, data_DE, linewidth=0.4)
    # plt.xlim((0.4, 0.6))
    plt.xlabel('time/s')
    plt.ylabel('magnitude/(m/(s^2))')
    plt.title('data_DE')
    plt.subplot(1, 2, 2)
    plt.plot(t, data_FE, linewidth=0.4)
    # plt.xlim((0.4, 0.6))
    plt.xlabel('time/s')
    plt.ylabel('magnitude/(m/(s^2))')
    plt.title('data_FE')
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

3. 正常信号的平方包络解调分析

平方解调的步骤为：求取信号平方，低通滤波。对信号取平方可得：

\begin{aligned} z_{s} (t) & = x_{m}^{2} (t) = A_{m}^{2} {[1 + α_{m} \cos (2 π f_{m} t)]}^{2} {[\cos (2 π f_{z} t)]}^{2} \\ = A_{m}^{2} [\frac{1 + α_{m}^{2}}{2} + 2 α_{m} \cos (2 π f_{m} t) + \frac{α_{m}^{2}}{2} \cos (2 π 2 f_{m} t)] [\frac{1}{2} + \frac{1}{2} \cos (2 π 2 f_{z} t)] \end{aligned}

$\begin{aligned} z_{s}(t) &=x_{m}^{2}(t)=A_{m}^{2}\left[1+\alpha_{m} \cos \left(2 \pi f_{m} t\right)\right]^{2}\left[\cos \left(2 \pi f_{z} t\right)\right]^{2} \\ &=A_{m}^{2}\left[\frac{1+\alpha_{m}^{2}}{2}+2 \alpha_{m} \cos \left(2 \pi f_{m} t\right)+\frac{\alpha_{m}^{2}}{2} \cos \left(2 \pi 2 f_{m} t\right)\right]\left[\frac{1}{2}+\frac{1}{2} \cos \left(2 \pi 2 f_{z} t\right)\right] \end{aligned}$

z_{s} (t) = x_{m}^{2} (t) = A_{m}^{2} [1 + α_{m} cos (2 π f_{m} t)]^{2} [cos (2 π f_{z} t)]^{2} = A_{m}^{2} [\frac{1 + α _{m}^{2}}{2} + 2 α_{m} cos (2 π f_{m} t) + \frac{α _{m}^{2}}{2} cos (2 π 2 f_{m} t)] [\frac{1}{2} + \frac{1}{2} cos (2 π 2 f_{z} t)]

根据频域卷积定理，频谱中频率成分包括含有调制频率 $f_{\mathrm{m}}$ 及其 2 倍频谐波的低频部分，以及包含有以 2 倍载波频率 $f_{x}$ 为中心频率，以调制频率 $f_{m}$ 及其 2 倍频谐波频率为边带间隔频率的高频调制部分。

我们选取213DE信号进行包络解调分析，第125号数据对应的工况是转速为1797rpm、轴承内圈故障（故障尺寸为0.021inchs）。根据第一章节描述，内圈故障特征频率为 $5.42*1797\div60=162.3Hz$ ,通过包络解调也正好可得到该频率，证明了我们算法的准确性。

在这里插入图片描述

图3.1 213包络解调结果

source = loadmat('.//data//48k Drive End Bearing Fault Data/213.mat')
fs = 48000
t, data_DE, data_FE, data_Speed = get_data(source, fs)

# 信号解包络
sig_envelope = data_DE[:4*fs] ** 2
# 对信号进行频域低通滤波(100Hz)
low_pass = 400
# 求取包络信号的FFT变换结果
sig_fft = fftpack.fft(sig_envelope)
# 求取包络信号的FFT变换频率
sig_fre = fftpack.fftfreq(len(sig_fft), d=1/fs)
# 低通滤波器
sig_fft[abs(sig_fre)>low_pass] = 0
# 转换为滤波后时域信号
sig_filter = fftpack.ifft(sig_fft).real
# 求取包络曲线的FFT结果
sig_fre, sig_mag = signal.welch(sig_filter, fs=fs, nfft=4*fs, 
                                nperseg=4*fs, scaling = 'spectrum')
sig_mag = np.sqrt(2 * sig_mag)

# 展示结果
plt.figure(figsize=(10, 4))
plt.plot(sig_fre, sig_mag)
plt.xlabel('frequence(Hz)')
plt.ylabel('magnitude')
plt.vlines(sig_fre[np.argmax(sig_mag)], 0, np.max(sig_mag), 'r', 'dashdot')
plt.text(sig_fre[np.argmax(sig_mag)], np.max(sig_mag), 
         str(sig_fre[np.argmax(sig_mag)]) + 'Hz', c = 'r')
plt.xlim((120, 180))
plt.ylim((0, 0.25))
plt.title('envelope signal-frequence')
plt.show()
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

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/繁依Fanyi0/article/detail/202608