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

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()

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()

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()

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

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

$$\begin{array}{rl}{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 2f_{m}t\right)\right]\left[\frac{1}{2}+\frac{1}{2}\cos \left(2\pi 2f_{z}t\right)\right]\end{array}$$

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

  我们选取213DE信号进行包络解调分析， 第125号数据对应的工况是转速为1797rpm、轴承内圈故障（故障尺寸为0.021inchs）。根据第一章节描述，内圈故障特征频率为$5.42\ast 1797÷60=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()