经验模态分解

经验模态分解(Empirical Mode Decomposition, EMD)是依据数据自身的时间尺度特征来进行信号分解，无需预先设定任何基函数，是一种时频域信号处理方式。EMD在处理非平稳及非线性数据上具有明显的优势，适合分析非线性非平稳的信号序列，具有较高的信噪比。
平稳信号：分布参数（均值，方差，协方差等）或者分布规律随着时间不发生变化
非平稳信号：分布参数（均值，方差，协方差等）或者分布规律随着时间发生变化

该方法的关键是经验模式分解，使负载信号分解为有限个本征信号（Intrinsic Mode Function, IMF），分解出来的各个IMF分量包含了原信号的不同时间尺度的局部特征信息。通俗理解就是，EMD是一台机器，将一对混杂的硬币投进去，其会自动按照1元、5毛、1毛、5分、1分分好。

EMD原理

在物理上，如果瞬时频率有意义，那么函数必须是对称的，局部均值为零，并且具有相同的过零点和极值点数目。任何信号都是由若干本征模函数组成，一个本征模函数必须满足以下两个条件：

  • 函数在整个时间范围内，局部极值点和过零点的数目必须相等或最多相差一个
  • 在任意时刻点，局部最大值的包络（上包络线）和局部最小值的包络（下包络线）平均必须为零，即上下包络线相对于时间轴局部对称

画图解释：

  1. 图线要反复跨越x轴

    img

    在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个

    而不能像这样某次穿过零点后出现多个极点：

    img

极点数目偏多

  1. 包络线要对称：

    img

    包络线对称

    而不能像这样：

    img

    包络线不对称

对于上述第二条说明：他把经典的全局性要求修改为局部性要求，使瞬时频率不再受不对称波形所形成的不必要的波动所影响，实际上，这个条件应为“数据的局部均值是零”。但是对于非平稳数据来说，计算局部均值涉及到“局部时间尺度”的概念，而这是很难定义的，因此，在第二个条件中使用了局部极大值和局部极小值包络的平均为零来代替，是信号的波形局部对称。

EMD将输入信号分解为几个本征模函数和一个残差组成，即由下列公式组成：

$$
I ( n ) = \sum _ { m = 1 } ^ { M } \operatorname { I M F _ { m } ( n ) + \operatorname { Res } _{m} ( n ) }
$$
其中I(n)表示输入信号，IMFm表示Mth的本征模函数，ResM(n)表示残差

EMD分解过程

提取IMF的过程称为筛选，筛选的过程如下：

emd分解

  1. 标出局部极值点
  2. 通过三次样条插值（cubic spline line）连接极大值点构成上包络线（upper envelope）,连接极小值点构成下包络线（lower envelope）
  3. 求上下包络线的均值m1
  4. 用输入信号减去上下包络线均值

$$
X ( t ) - m _ { 1 } = h _ { 1 }
$$

​ 上述过程的一次迭代不能保证h1是本征模函数（IMF），需要重复上述过程，直到h1是本征模函数（IMF）

停机准则

停机准则决定了一个本征模函数（IMF）筛选过程执行的数目，有如下停机准则

  • 标准偏差（Standard Deviation, SD）
    $$
    S D _ { k } = \sum _ { t = 0 } ^ { T } \frac { | h _ { k - 1 } ( t ) - h _ { k } ( t ) | ^ { 2 } } { h _ { k - 1 } ^ { 2 } (t)}
    $$
    当SD的值小于给定的阈值时，筛选过程停止

  • S Number准则

    定义为过零点和极值点相等或者至多差为1的连续筛选数目。一个S-Number被提前设置，只有当S次连续筛选后，每一次过零点和极值点保持相同（相等或者至多差1），筛选过程才停止

  • 阈值方法

    阈值方法设置两个阈值，确保全局小的扰动同时考虑局部大的偏移

选择停机准则后，第一个IMF（c1）可以获得，c1为包含输入信号最大频率的成分（component）,之后分离c1
$$
X ( t ) - c _ { 1 } = r _ { 1 }
$$
利用r1作为输入，获得其他的本征模函数

EMD限制（Limitations）

EMD的主要优点有如下：

  • EMD具有数据驱动的自适应性，能分析非线性非平稳信号，不受Heisenberg测不准原理制约等优点。
  • EMD在非线性非平稳信号分析中具有显著优势。与传统分析技术相比，EMD无需选择基函数，其分解基于信号本身极值点的分布。

EMD的主要缺点有如下：

  • 末端效应

    末端效应发生在信号的开始和结尾，因为在信号开始之前和结尾之后没有样本点被考虑。大多数情况下，末端点并不是信号的极值，但是在执行EMD的过程中，极值包络线会在末端点发散（diverge），导致错误，进而扭曲了IMF在末端点的波形，而且这种错误在EMD分解过程中会累积

  • 模态混叠问题

    模态混叠问题发上在EMD的执行过程中。出现下列情况之一就称为模态混合

    • 在同一个IMF分量中，存在尺度分布范围很宽却又各不相同的信号
    • 在不同的IMF分量中，存在着尺度相近的信号

    image-20211227120152132

    模态混叠问题使得特征提取、模型训练、模式识别变得困难，IMF失去了单一特征尺度的特征。集成经验模态分解(Ensemble empirical mode decomposition, EEMD)被提出用来解决模态混叠问题

EMD和其他方法对比

结果对比

集成经验模态分解（EEMD）

为了改善测量的准确性，集成平均是有效的方法（也就是多次测量取平均值）。信号极值点影响IMF,若分布不均匀时会出现模态混叠，白噪声的频谱均匀分布，白噪声是的信号会自动分布到合适的参考尺度上。由于零均值噪声的特性，噪音经过多次的平均计算后会相互抵消，这样集成均值的计算结果与原始信号的差值随着集成平均的次数增加而减少

EEMD的分解

  1. 给原始信号添加白噪声序列
  2. 分解带有白噪声序列的输入信号，得到IMFs（第一次测量得到一系列的IMFs）
  3. 重复第一和第二步，每次添加不同的白噪声序列（执行多次测量）
  4. 获取相关IMFs集成的均值作为最后的结果（理解为多次测量取平均）

添加的白噪声的两个属性

  • 添加的白噪声导致所有时间尺度上机制分布的相对均匀分布（The added white noise leads to relatively even distribution of extrema distribution on all timescales）
  • 通过集成平均，添加的噪声会被移除（噪声的均值为0）

集成经验模态分解（EEMD）优缺点

优点：

  • 该算法利用EMD滤波器组行为及白噪声频谱均匀分布的统计特性，使Sifting过程信号极值点分布更趋匀称，有效抑制由间歇性高频分量等因素造成的模态混叠。

缺点：

  • 在EEMD中，每个加噪信号hi(t)独立地被分解，使得每个hi(t)分解后可能产生不同数量的IMF，导致集合平均时IMF分量对齐困难。
  • 此外，添加的白噪声幅值和迭代次数依靠人为经验设置，当数值设置不当时，无法克服模态混叠。
  • 集总平均次数一般在几百次以上，非常耗时。虽然增加集合平均次数可降低重构误差，但这是以增加计算成本为代价，且有限次数的集合平均并不能完全消除白噪声，导致算法重构误差大，分解完备性差。

代码实现EMD

  • 安装EMD包

    1
    pip install EMD-signal

  • 测试代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    values = dataset['request'].tolist()
    S = np.array(values)
    emd = EMD()
    # emd.emd(S,max_imf=3)
    emd.emd(S)
    imfs, res = emd.get_imfs_and_residue()

    # 绘制 IMF
    vis = Visualisation()
    vis.plot_imfs(imfs=imfs, residue=res, include_residue=True)
    # 绘制并显示所有提供的IMF的瞬时频率
    # vis.plot_instant_freq(imfs=imfs)
    vis.show()

代码实现EEMD

  • 测试代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    # Prepare and run EEMD
    eemd = EEMD(trials=50)
    eemd.noise_seed(12345)
    E_IMFs = eemd.eemd(S)
    imfNo = E_IMFs.shape[0]
    # Plot results in a grid
    c = np.floor(np.sqrt(imfNo + 1))
    r = np.ceil((imfNo + 1) / c)
    plt.ioff()
    plt.subplot(r, c, 1)
    plt.plot(S, "r")
    plt.title("Original signal")

    for num in range(imfNo):
        plt.subplot(r, c, num + 2)
        plt.plot(E_IMFs[num], "g")
        # plt.xlim((tMin, tMax))
        plt.title("Imf " + str(num + 1))
    plt.show()

    image-20211227132458310

参考：