python实战故障诊断之CWRU数据集（三）：信号预白化处理-倒谱预白化（CEP pre-whitening）

1. 概述

  在轴承故障中，故障信号通常较为微弱，很可能被设备周期性运转产生的强周期性信号所淹没，导致无法准确识别故障，这时我们需要用到信号预白化方法来提取故障信息。信号预白化的思路是通过抑制振动信号中的强周期部分，提升振动信号的冲击特性，从而提取出故障信息，目前常用的信号预白化方法有两种：

• 通过线性预测方法实现， 如基于信号自回归模型分离振动信号的确定性分量，留下的残余信号包含白噪声和轴承损伤引起的非平稳冲击;
• 通过倒谱编辑实现信号的预白化。

  本章主要介绍基于倒频谱的信号预白化方法。

2. 倒频谱预白化原理

2.1 倒频谱

  倒谱定义为“对数谱的逆傅里叶变换”。 而复倒谱定义为“复对数谱的逆傅里叶变换”。 显然复倒谱的优点是与时间信号可逆， 这样使得在倒谱域的信号处理成为可能，但是需要将相位展开成频率的连续函数。倒频谱的定义如下：
C ( τ ) = F − 1 [ ln ⁡ ( X ( f ) ) ] X ( f ) = ℑ [ x ( t ) ] = A ( f ) exp ⁡ ( j φ ( f ) ) (2-1)

$$\begin{gathered} C(\tau)=F^{-1}[\ln (X(f))] \\ X(f)=\Im[x(t)]=A(f) \exp (\mathrm{j} \varphi(f)) \end{gathered}$$ \tag{2-1} C(τ)=F−1[ln(X(f))]X(f)=ℑ[x(t)]=A(f)exp(jφ(f))​(2-1)
  写成幅值和相位的形式：
$$\ln (X(f))=\ln (A(f))+j\phi (f)\text{\hspace{1em}(2-2)}$$
  $X(f)$ 是复数时, $C(\tau )$ 表示的是复倒谱。

2.2 基于倒频谱的预白化处理

  基于倒频谱的振动信号预白化框图如图 2.1 所示。将除零倒频率处以外的倒谱分量均置零，此时对应的对数谱变成均匀分布，均值是原始对数幅值的平均值。倒频率零处的值确定对数谱的均值，因此给残余信号提供了准确的尺度。

图2.1 倒频谱预白化流程图

  预白化操作通过对整个实倒频谱设置一个零值(0除外)，然后逆变换回频域，得到的信号与原始信号的相位重新组合，最后逆变换到时域。这相当于围绕确定性周期信号进行一系列削减操作，结果几乎完全消除了它们对信号的影响。另一方面，与轴承损坏相关的信号，并不是严格周期性的，在倒频谱的绝对值中不会出现任何强峰值，也不会受到削减操作的影响。这个过程也可以用一种非常简单的方式来实现，避免向倒谱域的变换，只需要将傅里叶变换后的信号除以它的绝对值然后再变换回时域。
$$x_{cpw}=IFFT\left\{\frac{FT(x)}{|FT(x)|}\right\}\text{\hspace{1em}(2-3)}$$

3. 代码实现

  我们选取273FA信号进行包络解调分析， 第273号数据对应的工况是转速为1730rpm、风扇端轴承内圈故障（故障尺寸为0.021inchs）。根据第一章节描述，内圈故障特征频率为$4.95\times 1730÷60=142.7Hz$。通过预白化分析最终得到故障频率为143.5Hz(轴承运转过程中可能有滑移，导致频率漂移)，顺利获得故障频率。

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_BA = []
    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:
            if len(source.keys()) == 6:
                data_Speed = source[key].flatten()
            else:
                data_BA = source[key].flatten()
        elif i == 6:
            data_Speed = source[key].flatten()
                
    t = np.arange(len(data_DE))/fs
    return t, data_DE, data_FE, data_BA, data_Speed


def cep_prewhite(sig):

    # 求取信号的倒频谱
    sig_fft = fftpack.fft(sig)
    sig_fft_log = np.log(sig_fft)
    sig_rceps = fftpack.ifft(abs(sig_fft_log))
    # 对倒频谱中除0外全部置零
    sig_rceps[1:] = np.zeros(len(sig_rceps[1:]))
    # 反求原始信号（保留原始信号相位）
    _res = np.exp(fftpack.fft(sig_rceps) + 1j * sig_fft_log.imag)
    result = fftpack.ifft(_res).real
    '''
    # 等效简便算法
    sig_fft = fftpack.fft(sig)
    result = fftpack.ifft(sig_fft/abs(sig_fft)).real
    '''
    return result

# 获取待处理数据
source = loadmat('.//data/12k Fan End Bearing Fault Data/273.mat')
# 设定采样频率
fs = 12000
# 解析数据
t, data_DE, data_FE, data_BA, data_Speed = get_data(source, fs)
sig = data_FE[:4*fs]
# 信号预白化处理(去除边界点)
sig_prewhite = cep_prewhite(sig)[100:-100]

plt.figure(figsize=(10, 4))
plt.subplots_adjust(hspace=0.6)
plt.subplot(2, 1, 1)
plt.plot(np.arange(len(sig))/fs, sig)
plt.xlabel('time/s')
plt.ylabel('magnitude/(m/(s^2))')
plt.title('signal_orignal')
plt.xlim((0, 0.05))
plt.subplot(2, 1, 2)
plt.plot(np.arange(len(sig_prewhite))/fs, sig_prewhite)
plt.xlabel('time/s')
plt.ylabel('magnitude/(m/(s^2))')
plt.title('signal_prewhite')
plt.xlim((0, 0.05))
plt.show()

# 信号解包络
sig_envelope = sig_prewhite ** 2
# 对信号进行频域低通滤波(200Hz)
low_pass = 200
# 求取包络信号的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=2*fs, 
                                nperseg=2*fs, scaling = 'spectrum')
sig_mag = np.sqrt(2 * sig_mag)

# 展示包络谱频谱特征
plt.figure(figsize=(10, 4))
plt.plot(sig_fre, sig_mag, c='black')
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((0, 200))
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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

图2.2 预白化前后对比

图2.3 预白化后解包络结果