从脑电波到股票预测：变分模态分解（VMD）在Python里的3个实战应用

从脑电波到股票预测：变分模态分解（VMD）在Python里的3个实战应用

变分模态分解（VMD）作为一种自适应信号处理方法，近年来在多个领域展现出强大的交叉应用潜力。不同于传统傅里叶变换或小波分析，VMD能够将复杂信号分解为一系列具有稀疏性的本征模态函数（IMF），每个IMF都具有明确的中心频率和有限带宽。这种特性使得VMD特别适合处理非平稳、非线性的复杂信号，从生物医学到金融分析，都能看到它的身影。

本文将聚焦三个典型场景：脑电信号分析、机械故障诊断和股票价格预测。每个案例都将展示如何通过Python实现VMD分解，并解释关键参数的选择逻辑。我们不会停留在理论层面，而是直接进入实战环节，通过代码和可视化结果，让你直观感受VMD在不同领域的应用价值。

1. 脑电信号分析：提取特定节律成分

脑电信号（EEG）是研究大脑活动的重要窗口，但原始EEG信号往往包含多种节律成分（如α波、β波等）以及各种噪声。传统滤波方法难以准确分离这些成分，而VMD提供了一种更智能的解决方案。

1.1 数据准备与预处理

典型的EEG信号采样频率在250-1000Hz之间。我们使用BCI竞赛数据集中的一段睁眼静息态EEG作为示例：

import numpy as np from scipy.io import loadmat import matplotlib.pyplot as plt # 加载EEG数据 eeg_data = loadmat('eeg_sample.mat')['data'][0] fs = 256 # 采样频率256Hz t = np.arange(len(eeg_data)) / fs # 绘制原始信号 plt.figure(figsize=(12,4)) plt.plot(t, eeg_data) plt.xlabel('Time (s)') plt.ylabel('Amplitude (μV)') plt.title('Raw EEG Signal') plt.show()

提示：EEG信号通常需要进行去趋势和带通滤波预处理，但VMD对预处理的要求相对较低，这是其优势之一。

1.2 VMD参数选择与分解

针对EEG信号，我们关注4-30Hz范围内的主要节律（θ、α、β波），因此设置K=4个模态：

from vmdpy import VMD alpha = 2000 # 带宽约束 tau = 0 # 噪声容限 K = 4 # IMF数量 DC = 0 # 不提取直流分量 init = 1 # 均匀初始化频率 tol = 1e-7 # 收敛容差 u, u_hat, omega = VMD(eeg_data, alpha, tau, K, DC, init, tol)

关键参数说明：

• alpha：控制带宽，值越大带宽越小。EEG信号通常选择2000-3000
• K：根据目标节律数量设置，通常4-6个足够覆盖主要脑电节律

1.3 结果分析与可视化

# 绘制各IMF及其频谱 plt.figure(figsize=(12,8)) for i in range(K): plt.subplot(K, 2, 2*i+1) plt.plot(t, u[i], label=f'IMF {i+1}') plt.legend() plt.subplot(K, 2, 2*i+2) freqs = np.linspace(0, fs/2, len(u[i])//2) fft_vals = np.abs(np.fft.fft(u[i]))[:len(u[i])//2] plt.plot(freqs, fft_vals) plt.xlim(0, 40) plt.xlabel('Frequency (Hz)') plt.tight_layout() plt.show()

典型结果会显示不同IMF对应不同频段：

• IMF1：δ波（0.5-4Hz）
• IMF2：θ波（4-8Hz）
• IMF3：α波（8-13Hz）
• IMF4：β波（13-30Hz）

这种分解可以帮助研究者单独分析特定脑电节律，而不受其他成分干扰。

2. 机械故障诊断：轴承振动信号分析

旋转机械的故障诊断依赖于对振动信号的特征提取。当轴承出现故障时，会产生特定的冲击特征频率，但这些频率往往被强噪声和机械运转的基本频率所掩盖。

2.1 轴承故障信号特点

典型的轴承故障信号包含：

• 基本旋转频率（通常较低）
• 故障特征频率（由损伤点周期性撞击产生）
• 宽带噪声

我们使用凯斯西储大学轴承数据中心的数据作为示例：

import pandas as pd # 加载轴承数据 bearing_data = pd.read_csv('bearing_fault.csv').values.flatten() fs = 12000 # 采样频率12kHz t = np.arange(len(bearing_data)) / fs # 绘制时域波形 plt.figure(figsize=(12,4)) plt.plot(t, bearing_data) plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.title('Bearing Vibration Signal') plt.show()

2.2 VMD参数优化

对于机械振动信号，我们需要更精细的频带划分：

alpha = 1500 # 稍小的alpha以获得更宽频带 tau = 0.1 # 允许一定噪声 K = 5 # 增加IMF数量 DC = 0 init = 1 tol = 1e-6 # 更严格的容差 u, u_hat, omega = VMD(bearing_data, alpha, tau, K, DC, init, tol)

2.3 故障特征提取

通过包络谱分析提取故障特征：

from scipy.signal import hilbert plt.figure(figsize=(12,8)) for i in range(K): # 计算包络谱 analytic_signal = hilbert(u[i]) env = np.abs(analytic_signal) env_fft = np.abs(np.fft.fft(env))[:len(env)//2] freqs = np.linspace(0, fs/2, len(env)//2) plt.subplot(K,1,i+1) plt.plot(freqs, env_fft) plt.xlim(0, 1000) # 聚焦在故障特征频率范围 plt.ylabel(f'IMF {i+1}') plt.xlabel('Frequency (Hz)') plt.tight_layout() plt.show()

通常会在某个IMF的包络谱中清晰地看到故障特征频率（如轴承外圈故障频率BPFO），这为故障诊断提供了直接依据。

3. 金融时间序列分析：股票价格预测

金融时间序列具有高度非线性和非平稳特性，传统时间序列分析方法往往难以捕捉其复杂动态。VMD可以将价格序列分解为不同时间尺度的成分，分别建模后再集成预测。

3.1 股票数据特点

以苹果公司(AAPL)的日收盘价为例：

import yfinance as yf # 获取股票数据 data = yf.download('AAPL', start='2020-01-01', end='2023-12-31') prices = data['Close'].values dates = data.index # 绘制价格曲线 plt.figure(figsize=(12,4)) plt.plot(dates, prices) plt.title('AAPL Daily Close Price') plt.ylabel('Price ($)') plt.grid(True) plt.show()

3.2 VMD分解策略

金融时间序列分解需要特别注意：

• 趋势成分的提取
• 季节性/周期性成分的分离
• 噪声的处理

# 归一化处理 prices_norm = (prices - np.mean(prices)) / np.std(prices) alpha = 1000 # 较小的alpha适应金融数据的宽频特性 tau = 0.01 # 允许一定重构误差 K = 3 # 通常分解为趋势、周期和噪声 DC = 1 # 提取直流分量作为趋势 init = 1 tol = 1e-5 u, u_hat, omega = VMD(prices_norm, alpha, tau, K, DC, init, tol)

3.3 预测模型构建

对分解后的各IMF分别建立预测模型：

from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import mean_squared_error # 准备训练数据 X = [] for i in range(len(prices)-10): X.append(prices_norm[i:i+10]) # 使用前10天数据预测第11天 X = np.array(X) # 对每个IMF训练预测模型 models = [] for k in range(K): y = u[k][10:] # 对应的目标值 model = RandomForestRegressor(n_estimators=100) model.fit(X, y) models.append(model) # 集成预测 test_input = prices_norm[-10:].reshape(1,-1) preds = [model.predict(test_input)[0] for model in models] final_pred = np.sum(preds) * np.std(prices) + np.mean(prices)

这种分解-集成的方法通常比直接对原始序列建模获得更好的预测性能，因为不同成分可能遵循不同的动态规律。

4. VMD参数选择经验总结

通过上述三个案例，我们可以总结一些VMD参数选择的实用经验：

其他实用技巧：

• 初始化策略：对于周期性明显的信号，可以使用init=2并指定初始频率
• 收敛判断：可以监控omega的变化，确保中心频率已稳定
• 带宽控制：如果模态混叠严重，适当增加alpha；如果信号丢失，则减小alpha

# 监控中心频率变化的实用函数 def plot_omega(omega): plt.figure(figsize=(10,6)) for i in range(omega.shape[1]): plt.plot(omega[:,i], label=f'IMF {i+1}') plt.xlabel('Iteration') plt.ylabel('Center Frequency') plt.legend() plt.show()

在实际项目中，建议先用少量数据测试不同参数组合，观察分解效果后再扩展到全数据集。记住，没有放之四海而皆准的最优参数，理解自己数据的特点才是关键。