第100+32步 ChatGPT学习：时间序列EMD分解

基于Python 3.9版本演示

一、写在前面

之前我们介绍过时间序列的季节性分解。

最近又学到了好几种骚操作分解，且可以用这些分解优化时间序列预测性能。

首先，我们来学一学经验模态分解（Empirical Mode Decomposition，EMD）。

二、经验模态（EMD）分解

时间序列的经验模态分解（Empirical Mode Decomposition，EMD）是一种数据驱动的信号处理技术，用于分解非线性、非平稳信号。EMD技术由诺贝尔奖获得者黄锷（Norden E. Huang）和他的同事们在1998年提出，它特别适用于分析自然界中复杂的时间序列数据。

EMD的基本思想是将一个复杂的信号分解成若干个具有不同时间尺度的内禀模态函数（Intrinsic Mode Functions，IMFs）和一个残余信号。这些IMFs代表信号中的不同振荡模式，而残余信号则反映了信号的整体趋势。

（1）EMD的优点

1）自适应性：EMD是一种完全自适应的方法，不需要任何预设参数，能够处理非线性、非平稳信号。。

2）直观性：分解得到的IMFs通常具有物理意义，便于理解和解释信号中的不同成分。

3）广泛应用：EMD已在气象学、地震学、生物医学信号处理、机械故障诊断等领域得到了广泛应用。

（2）EMD的缺点

1）边界效应：由于边界处数据不足，包络线构建过程中可能出现误差，影响分解结果。

2）模式混叠：在某些情况下，分解得到的IMFs可能会包含多个不同频率的成分，导致模态混叠问题。

3）计算复杂度：EMD的计算过程较为复杂，尤其是对长时间序列数据，计算时间较长。

三、EMD代码Pyhton实现

下面，我使用的是之前分享过的肺结核的数据做演示：

Python代码：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

# 读取数据
file_path = 'pone.0277314.s006.xlsx'
data = pd.read_excel(file_path)

# 提取时间和PTB病例数
time_series = data['Time']
ptb_cases = data['PTB cases']

# 将时间转换为数值形式
time_numeric = np.arange(len(time_series))

def get_envelope_mean(signal):
    """计算信号的上包络线和下包络线的均值"""
    maxima = np.where(np.r_[True, signal[1:] > signal[:-1]] & np.r_[signal[:-1] > signal[1:], True])[0]
    minima = np.where(np.r_[True, signal[1:] < signal[:-1]] & np.r_[signal[:-1] < signal[1:], True])[0]
    
    if len(maxima) < 2 or len(minima) < 2:
        return np.zeros_like(signal)
    
    upper_env = CubicSpline(maxima, signal[maxima])(time_numeric)
    lower_env = CubicSpline(minima, signal[minima])(time_numeric)
    
    return (upper_env + lower_env) / 2

def sift(signal, max_iter=1000, tol=1e-6):
    """对信号进行sifting操作，提取IMF"""
    h = signal
    for _ in range(max_iter):
        m = get_envelope_mean(h)
        h1 = h - m
        
        if np.mean(np.abs(h - h1)) < tol:
            break
        h = h1
    
    return h

def emd(signal, max_imfs=10):
    """进行EMD分解"""
    residual = signal
    imfs = []
    for _ in range(max_imfs):
        imf = sift(residual)
        imfs.append(imf)
        residual = residual - imf
        
        if np.all(np.abs(residual) < 1e-6):
            break
    
    return np.array(imfs), residual

# 执行EMD分解
imfs, residual = emd(ptb_cases.values)

# 绘制分解结果
num_imfs = imfs.shape[0]
plt.figure(figsize=(12, 9))
for i in range(num_imfs):
    plt.subplot(num_imfs + 1, 1, i + 1)
    plt.plot(time_series, imfs[i], label=f'IMF {i + 1}')
    plt.legend()

# 绘制剩余信号
plt.subplot(num_imfs + 1, 1, num_imfs + 1)
plt.plot(time_series, residual, label='Residual')
plt.legend()
plt.tight_layout()
plt.show()

输出：