Python軸承故障診斷之經(jīng)驗(yàn)?zāi)B(tài)分解EMD原理介紹

1 經(jīng)驗(yàn)?zāi)B(tài)分解EMD原理介紹

1.1 EMD概述

經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）方法是一種自適應(yīng)信號(hào)時(shí)頻處理方法，特別適用于非線性、非平穩(wěn)信號(hào)的分析處理[1]。其本質(zhì)是一種對(duì)信號(hào)進(jìn)行分解的方法,將信號(hào)分解為各個(gè)相互獨(dú)立的成分的疊加,依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征來(lái)進(jìn)行信號(hào)分解，具備自適應(yīng)性。

EMD的優(yōu)點(diǎn)在于它是一種自適應(yīng)的、數(shù)據(jù)驅(qū)動(dòng)的分解方法，不需要預(yù)先假設(shè)信號(hào)的分布或結(jié)構(gòu)。這使得它適用于處理各種類型的信號(hào)，包括非線性和非平穩(wěn)信號(hào)。

EMD 認(rèn)為任何一個(gè)復(fù)雜序列都是由多個(gè)單頻率信號(hào)疊加而成，因此可以分解成若干個(gè) 本征模態(tài)函數(shù)（Intrinsic Mode Functions, IMF），IMF 的各個(gè)分量即代表了原始信號(hào)中的各頻 率分量，并按照從高頻到低頻的順序依次排列，這也是 IMF 的物理含義[2]。

1.2 本征模態(tài)函數(shù)IMF

本征模態(tài)函數(shù)（Intrinsic Mode Functions, IMF）就是原始信號(hào)被 EMD 分解之后得到的各層信號(hào)分量。任何信號(hào)都可以拆分成若干個(gè) IMF 之和。IMF 有兩個(gè)假設(shè)條件：

• 在整個(gè)數(shù)據(jù)段內(nèi)，極值點(diǎn)的個(gè)數(shù)和過(guò)零點(diǎn)的個(gè)數(shù)必須相等或相差最多不能超過(guò)一 個(gè)；
• 在任意時(shí)刻，由局部極大值點(diǎn)形成的上包絡(luò)線和由局部極小值點(diǎn)形成的下包絡(luò)線 的平均值為零，即上、下包絡(luò)線相對(duì)于時(shí)間軸局部對(duì)稱。

對(duì)于上述的條件理解如下：

第一，圖線要反復(fù)跨越 x 軸，比如：

反復(fù)跨越 x 軸

而不能像下面這樣某次穿過(guò)零點(diǎn)后出現(xiàn)多個(gè)極點(diǎn)：

某次穿過(guò)零點(diǎn)后出現(xiàn)多個(gè)極點(diǎn)

第二，上下包絡(luò)線要對(duì)稱，比如：

上下包絡(luò)線對(duì)稱

而不能像如下這樣，上下包絡(luò)線不對(duì)稱

1.3 EMD 分解的基本假設(shè)

• 信號(hào)至少有兩個(gè)極值點(diǎn)：一個(gè)極大值點(diǎn)和一個(gè)極小值點(diǎn)；
• 特征時(shí)間尺度由極值之間的時(shí)間間隔定義；
• 如果數(shù)據(jù)完全沒有極值，但只包含拐點(diǎn)，那么可以一次或多次劃分來(lái)揭示極值點(diǎn)，最終的結(jié)果可以通過(guò)積分得到

2 EMD分解的基本原理和步驟

EMD的分解過(guò)程是一個(gè)迭代的過(guò)程。首先，對(duì)原始信號(hào)進(jìn)行極值點(diǎn)的提取，然后通過(guò)連接極值點(diǎn)的均值得到第一輪的近似IMF（也叫做“本征模態(tài)”）。接下來(lái)，將這個(gè)近似IMF從原始信號(hào)中減去，得到一個(gè)新的信號(hào)，然后對(duì)這個(gè)新信號(hào)再次進(jìn)行極值點(diǎn)提取和均值連接，得到第二輪的近似IMF。如此往復(fù)，直到得到的近似IMF滿足某種停止準(zhǔn)則。

對(duì)于原始信號(hào) X(t)

第一步，極值點(diǎn)提?。?/p>

從待分解的信號(hào)中識(shí)別局部極值點(diǎn)，包括局部極大值和局部極小值。極值點(diǎn)是信號(hào)中的局部特征，能夠幫助刻畫信號(hào)的振蕩特性。

第二步，構(gòu)建上下包絡(luò)線：

通過(guò)連接相鄰的局部極大值和局部極小值，構(gòu)建信號(hào)的上包絡(luò)線和下包絡(luò)線。上包絡(luò)線 U(t) 由局部極大值連接而成，下包絡(luò)線 L(t) 由局部極小值連接而成。包絡(luò)線用于描述信號(hào)的振蕩范圍。

第三步，提取均值函數(shù)：

計(jì)算上包絡(luò)線和下包絡(luò)線的平均值，得到均值函數(shù) m1。將原始信號(hào)減去均值函數(shù)，得到一維信號(hào) h1。

m1 = ( U(t) + L(t) ) / 2

h1 = X(t) - m1

第四步，迭代分解：

對(duì)減去均值函數(shù)后的一維信號(hào) h1，重復(fù)步驟1-3的過(guò)程，直到得到的剩余信號(hào)為“單調(diào)信號(hào)”，或者滿足IMF的兩個(gè)假設(shè)條件。迭代k次的IMF為

hk = h(k-1) - mk

第五步，確定本征模態(tài)函數(shù)（IMF）：

在每一次迭代中，通過(guò)極值點(diǎn)提取、構(gòu)建包絡(luò)線等步驟，得到的剩余信號(hào)被稱為一個(gè)本征模態(tài)函數(shù)（IMF）。IMF具有局部特征，并且代表了信號(hào)在不同尺度上的振蕩模式。使用上述方法得到的第一個(gè)IMF記為c1, 然后將c1從原始信號(hào)中分離，得到

r1 = X(t) - c1

由于r1仍然包含大量信息，因此將r1作為新的原始信號(hào)，再通過(guò)步驟1-4的分析，可以得到IMF2，以此類推，得到

r1 - c2 = r2,... ..., r(n-1) -cn = rn

當(dāng)cn或rn小于某一設(shè)定值，或者得到的剩余信號(hào)為“單調(diào)信號(hào)”，無(wú)法提取更多的IMF時(shí)，迭代終止，得到最終的分解結(jié)果為：

第六步，重構(gòu)信號(hào)：

將得到的IMF函數(shù)進(jìn)行逐個(gè)提取，直到無(wú)法再得到新的IMF為止。最終得到的IMF函數(shù)可以被看作是信號(hào)在不同時(shí)間尺度上的振蕩模式，它們的組合可以重構(gòu)原始信號(hào)。

這些基本步驟構(gòu)成了EMD方法的核心流程，通過(guò)這些步驟，EMD可以將復(fù)雜的信號(hào)分解成不同尺度和頻率的本征模態(tài)函數(shù)，從而揭示信號(hào)的局部特征和振蕩模式。

3 基于Python的EMD實(shí)現(xiàn)

在 Python 中，使用 PyEMD 庫(kù)來(lái)實(shí)現(xiàn)經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）

2.1 代碼示例

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD


# 生成一個(gè)示例信號(hào)
t = np.linspace(0, 1, 1000)
s = np.sin(11*2*np.pi*t*t) + 6*t*t


# 創(chuàng)建 EMD 對(duì)象
emd = EMD()


# 對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解
IMFs = emd(s)


# 繪制原始信號(hào)和每個(gè)本征模態(tài)函數(shù)（IMF）
plt.figure(figsize=(15,10))
plt.subplot(len(IMFs)+1, 1, 1)
plt.plot(t, s, 'r')
plt.title("Original signal")


for num, imf in enumerate(IMFs):
    plt.subplot(len(IMFs)+1, 1, num+2)
    plt.plot(t, imf)
    plt.title("IMF "+str(num+1))


plt.show()

2.2 軸承故障數(shù)據(jù)的分解

選擇 0.021英寸滾珠故障信號(hào)數(shù)據(jù)來(lái)做EMD分解

2.2.1 凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)的加載

第一步，導(dǎo)入包，讀取數(shù)據(jù)

import numpy as np
from scipy.io import loadmat
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rc("font", family='Microsoft YaHei')


# 讀取MAT文件
data = loadmat('21_2.mat') # 0.021英寸 滾珠
# 注意，讀取出來(lái)的data是字典格式，可以通過(guò)函數(shù)type(data)查看。

第二步，數(shù)據(jù)集中讀取 驅(qū)動(dòng)端加速度數(shù)據(jù)，取一個(gè)長(zhǎng)度為1024的信號(hào)進(jìn)行后續(xù)觀察和實(shí)驗(yàn)

# DE - drive end accelerometer data 驅(qū)動(dòng)端加速度數(shù)據(jù)
data_list = data['X222_DE_time'].reshape(-1)
# 劃窗取值（大多數(shù)窗口大小為1024）
data_list = data_list3[0:1024]
#  進(jìn)行數(shù)據(jù)可視化
plt.figure(figsize=(20,10))
plt.plot(data_list)
plt.title("滾珠")
plt.show()

2.2.2 滾珠故障信號(hào)EMD分解

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD


t = np.linspace(0, 1, 1024)
data = np.array(data_list)
# 創(chuàng)建 EMD 對(duì)象
emd = EMD()


# 對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解
IMFs = emd(data)


# 繪制原始信號(hào)和每個(gè)本征模態(tài)函數(shù)（IMF）
plt.figure(figsize=(15,10))
plt.subplot(len(IMFs)+1, 1, 1)
plt.plot(t, data, 'r')
plt.title("Original signal", fontsize=10)


for num, imf in enumerate(IMFs):
    plt.subplot(len(IMFs)+1, 1, num+2)
    plt.plot(t, imf)
    plt.title("IMF "+str(num+1), fontsize=10)
    # 增加第一排圖和第二排圖之間的垂直間距
plt.subplots_adjust(hspace=0.4, wspace=0.2)
plt.show()

2.3 信號(hào)分量的處理

通過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）得到了信號(hào)的分量，可以進(jìn)行許多不同的分析和處理操作，以下是一些常見的對(duì)分量的利用方向：

（1）信號(hào)重構(gòu)：將分解得到的各個(gè)本征模態(tài)函數(shù)（IMF）相加，可以重構(gòu)原始信號(hào)。這可以用于驗(yàn)證分解的效果，或者用于信號(hào)的重建和恢復(fù)。

（2）去噪：對(duì)于復(fù)雜的信號(hào)，可能存在噪聲或干擾成分。通過(guò)分析各個(gè)IMF的頻率和振幅，可以識(shí)別和去除信號(hào)中的噪聲成分。

（3）頻率分析：分析每個(gè)IMF的頻率成分，可以幫助理解信號(hào)在不同頻率上的振蕩特性，從而揭示信號(hào)的頻域特征。

（4）特征提取：每個(gè)IMF代表了信號(hào)的局部特征和振蕩模式，可以用于提取信號(hào)的特征，并進(jìn)一步應(yīng)用于機(jī)器學(xué)習(xí)或模式識(shí)別任務(wù)中。

（5）信號(hào)預(yù)測(cè)：通過(guò)對(duì)分解得到的各個(gè)IMF進(jìn)行分析，可以探索信號(hào)的未來(lái)趨勢(shì)和發(fā)展模式，從而用于信號(hào)的預(yù)測(cè)和預(yù)測(cè)建模。

（6）模式識(shí)別：分析每個(gè)IMF的時(shí)域和頻域特征，可以幫助對(duì)信號(hào)進(jìn)行模式識(shí)別和分類，用于識(shí)別信號(hào)中的不同模式和特征。

（7）異常檢測(cè)：通過(guò)分析每個(gè)IMF的振幅和頻率特征，可以用于探測(cè)信號(hào)中的異?；蛲话l(fā)事件，從而用于異常檢測(cè)和故障診斷。

在得到了信號(hào)的分量之后，可以根據(jù)具體的應(yīng)用需求選擇合適的分析和處理方法，以實(shí)現(xiàn)對(duì)信號(hào)的深入理解、特征提取和應(yīng)用。對(duì)于后續(xù)的研究，主要利用IMF分類來(lái)對(duì)故障信號(hào)做模式識(shí)別，即故障分類。