EEG信号分析实战用Python解码癫痫脑电波的数字指纹脑电图EEG就像大脑的“心电图”它记录着神经元群体此起彼伏的电活动交响曲。对于神经科学研究者、生物医学工程师或是任何对脑机接口、神经疾病辅助诊断感兴趣的开发者而言能够亲手处理和分析这些微伏级的信号无异于掌握了一门与大脑直接对话的语言。癫痫作为脑电活动异常同步化的典型代表其EEG中蕴含的“尖波”、“棘波”等特征波形正是我们试图从复杂背景噪声中捕捉的关键信号。本文不会重复教科书上的医学理论而是直接切入实战带你用Python搭建一套从原始EEG数据到特征识别的完整分析流水线。我们将聚焦于如何利用MNE-Python、PyEEG等开源工具包将抽象的脑电波形转化为可计算、可分析的数字特征一步步揭开癫痫脑电波的神秘面纱。1. 环境搭建与数据准备构建你的EEG分析工作台工欲善其事必先利其器。在开始任何信号分析之前一个稳定、高效的Python环境是基石。我强烈建议使用conda来管理你的环境它能很好地解决科学计算包之间复杂的依赖关系。# 创建一个名为eeg_analysis的独立环境并指定Python版本 conda create -n eeg_analysis python3.9 conda activate eeg_analysis # 安装核心数据分析库 conda install numpy scipy pandas matplotlib seaborn # 安装EEG处理的核心工具包MNE-Python pip install mne # 安装用于特征提取的PyEEG注意PyEEG可能需从源码安装或寻找替代 pip install pyEEG注意PyEEG库在某些系统上可能安装不顺利它是一个较早的用于提取多种非线性动力学特征的库。如果遇到问题可以考虑使用AntroPy或PyEntropy作为替代它们提供了更现代、维护更好的熵特征计算实现。准备好环境后下一步是获取数据。公开的癫痫EEG数据集是我们学习和验证算法的宝贵资源。其中最著名的是波恩大学癫痫EEG数据集。这个数据集结构清晰非常适合入门。import mne import numpy as np import pandas as pd import matplotlib.pyplot as plt # 假设我们已经将波恩数据集下载到本地并解压为文本文件 # 这里以加载一个片段为例实际中你需要根据数据格式调整 # 波恩数据集通常为单通道文本文件采样频率为173.61 Hz data_path ./data/epilepsy_eeg/ file_name Z001.txt # 使用numpy加载原始数据 raw_data np.loadtxt(data_path file_name) # 创建MNE的RawArray对象需要指定采样频率和信息结构 sfreq 173.61 # 波恩数据集的采样频率 info mne.create_info(ch_names[EEG], sfreqsfreq, ch_types[eeg]) raw mne.io.RawArray(raw_data.reshape(1, -1), info) # 简单查看数据基本信息 print(raw) print(f数据时长{raw.times[-1]:.2f} 秒) print(f数据点数{len(raw.times)})对于更复杂、多通道的标准EEG数据如.edf格式MNE提供了直接的读取函数# 加载一个标准的EDF格式EEG文件 edf_file ./data/sample_epilepsy.edf raw_edf mne.io.read_raw_edf(edf_file, preloadTrue) raw_edf.plot(duration10, scalingsauto) # 可视化前10秒数据下表对比了初学者常用的几个癫痫EEG公开数据集帮助您根据需求选择数据集名称主要特点数据格式适用场景波恩大学数据集单通道五类状态健康睁/闭眼癫痫病灶/非病灶区发作期结构简单。文本文件(.txt)算法验证、特征提取入门、二分类/多分类实验。CHB-MIT Scalp EEG多通道23通道长程头皮EEG包含发作标注临床价值高。EDF格式发作检测算法开发、长时间序列分析、多通道分析。TUH EEG Corpus目前最大的公开EEG数据库之一包含大量异常含癫痫记录数据多样。EDF格式开发鲁棒性强的模型、处理真实世界噪声。Kaggle癫痫预测数据集来自竞赛包含狗和人的颅内EEG旨在预测癫痫发作。Mat文件(.mat)发作预测模型、时序预测任务。2. EEG数据预处理从原始信号到干净波形原始的EEG信号几乎总是被各种噪声和伪迹所污染比如工频干扰50/60Hz、眼动EOG、肌电EMG以及基线漂移。不经过预处理的信号就像被雪花覆盖的电视屏幕重要的特征被完全掩盖。预处理的目标就是最大限度地保留大脑活动信号剔除这些干扰。滤波是预处理的第一步也是最关键的一步。EEG分析通常关注特定的频带如Delta: 0.5-4 Hz, Theta: 4-8 Hz, Alpha: 8-13 Hz, Beta: 13-30 Hz, Gamma: 30-45 Hz。癫痫的棘波、尖波通常包含较宽的频率成分但主要集中在Beta和Gamma频段。同时我们必须滤除无所不在的工频干扰。# 对之前加载的raw_edf数据进行滤波 raw_filtered raw_edf.copy() # 1. 首先进行高通滤波去除低频漂移例如低于0.5Hz的缓慢波动 raw_filtered.filter(l_freq0.5, h_freqNone) # 高通滤波 # 2. 进行低通滤波去除高频噪声并抑制工频干扰的谐波 raw_filtered.filter(l_freqNone, h_freq40.0) # 低通滤波到40Hz # 3. 专门针对50Hz工频干扰进行陷波滤波 raw_filtered.notch_filter(freqs50.0) # 对比滤波前后的功率谱密度 raw_edf.compute_psd(fmax50).plot(pickseeg, averageTrue, amplitudeFalse) raw_filtered.compute_psd(fmax50).plot(pickseeg, averageTrue, amplitudeFalse)除了滤波坏道检测与插值、伪迹剔除也是重要环节。眼电伪迹可以通过独立成分分析ICA来识别和去除。# 使用ICA去除眼电伪迹需要先进行高通滤波以利于ICA计算 raw_for_ica raw_filtered.copy().filter(l_freq1.0, h_freqNone) ica mne.preprocessing.ICA(n_components20, random_state97, max_iter800) ica.fit(raw_for_ica) # 自动寻找与眼电相关的成分 eog_indices, eog_scores ica.find_bads_eog(raw_for_ica, ch_name[Fp1, Fp2]) print(f检测到的眼电成分索引{eog_indices}) # 可视化这些成分的拓扑图和时间过程确认后将其剔除 ica.plot_properties(raw_for_ica, pickseog_indices[:2]) raw_cleaned raw_filtered.copy() ica.apply(raw_cleaned, excludeeog_indices)预处理流程可以总结为以下核心步骤但顺序和参数需根据具体数据调整数据导入与检查查看原始波形识别明显坏段。滤波高通如0.5Hz、低通如40Hz、陷波50/60Hz。重参考将参考电极转换为平均参考或其他生理学上更合理的参考方式。坏道处理检测并插值或剔除噪声过大的通道。伪迹剔除使用ICA或SSP等方法去除眼电、心电等伪迹。分段如果是事件相关电位分析需要根据事件标记进行分段。3. 特征工程提取癫痫脑电波的数字“指纹”经过预处理的干净EEG信号我们需要从中提取能够表征癫痫活动的特征。这些特征可以分为时域特征、频域特征和时频域特征。癫痫的棘波持续时间70ms和尖波持续时间70-200ms在时域上表现为快速的瞬态尖峰在频域上则会有特定的能量分布变化。时域特征直接刻画波形的形态。对于癫痫检测以下特征非常有用import numpy as np from scipy import stats, signal def extract_time_features(signal_segment): 从一个EEG信号片段中提取时域特征。 signal_segment: 一维numpy数组。 features {} # 1. 统计特征 features[mean] np.mean(signal_segment) features[std] np.std(signal_segment) # 标准差反映波动大小 features[variance] np.var(signal_segment) features[skewness] stats.skew(signal_segment) # 偏度衡量分布不对称性 features[kurtosis] stats.kurtosis(signal_segment) # 峰度衡量分布尖锐度 # 2. 基于幅值的特征对尖峰敏感 features[max_amplitude] np.max(np.abs(signal_segment)) features[line_length] np.sum(np.abs(np.diff(signal_segment))) # 线长对高频变化敏感 # 3. Hjorth参数刻画信号复杂性 diff1 np.diff(signal_segment) diff2 np.diff(diff1) activity np.var(signal_segment) mobility np.sqrt(np.var(diff1) / activity) complexity np.sqrt(np.var(diff2) / np.var(diff1)) / mobility features[hjorth_activity] activity features[hjorth_mobility] mobility features[hjorth_complexity] complexity return features # 应用将一段长信号分割成重叠的窗口并提取特征 def sliding_window_feature_extraction(raw_signal, window_sec2, overlap_ratio0.5, sfreq256): 滑动窗口特征提取。 window_points int(window_sec * sfreq) step_points int(window_points * (1 - overlap_ratio)) num_windows (len(raw_signal) - window_points) // step_points 1 feature_list [] for i in range(num_windows): start i * step_points end start window_points segment raw_signal[start:end] feat extract_time_features(segment) feature_list.append(feat) return pd.DataFrame(feature_list)频域特征揭示了信号在不同频率上的能量分布。癫痫发作间期或发作期特定频带如Gamma的能量可能会显著升高。def extract_frequency_features(signal_segment, sfreq): 提取频域特征。 features {} # 计算功率谱密度 freqs, psd signal.welch(signal_segment, fssfreq, npersegmin(256, len(signal_segment))) # 定义经典频带边界 band_def { delta: (0.5, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30), gamma: (30, 45) } total_power np.sum(psd) for band, (low, high) in band_def.items(): idx_band np.logical_and(freqs low, freqs high) band_power np.sum(psd[idx_band]) features[f{band}_power] band_power features[f{band}_relative_power] band_power / total_power if total_power 0 else 0 # 计算频带内功率谱的中值频率 if band_power 0: cumsum np.cumsum(psd[idx_band]) median_freq_idx np.where(cumsum band_power / 2)[0][0] features[f{band}_median_freq] freqs[idx_band][median_freq_idx] else: features[f{band}_median_freq] 0 # 谱熵 - 衡量频谱的平坦程度 psd_norm psd / np.sum(psd) spectral_entropy -np.sum(psd_norm * np.log2(psd_norm 1e-12)) features[spectral_entropy] spectral_entropy return features非线性动力学特征近年来备受关注它们试图刻画大脑作为一个复杂系统的特性。例如近似熵或样本熵可以量化信号的规律性癫痫发作时信号可能变得更规则或更不规则取决于发作类型。# 使用AntroPy库计算样本熵安装pip install antropy import antropy as ant def extract_nonlinear_features(signal_segment): features {} features[permutation_entropy] ant.permutation_entropy(signal_segment, order3, delay1, normalizeTrue) features[sample_entropy] ant.sample_entropy(signal_segment, order2, metricchebyshev) features[hjorth_fractal_dim] ant.hjorth_fd(signal_segment) # 分形维数 # 注意计算这些特征需要足够长的数据段且参数选择影响结果 return features在实际项目中我们通常会将数十甚至上百个特征组合成一个高维特征向量用于后续的机器学习模型。4. 构建癫痫检测模型从特征到分类提取了大量特征后我们面临两个核心问题1哪些特征最有效2用什么模型进行分类这里我们构建一个简单的机器学习流水线用于区分“发作间期”和“正常”或“发作期”的EEG片段。首先我们需要准备带标签的数据。以波恩数据集为例我们可以将SetA健康睁眼和SetB健康闭眼作为正常类将SetC、SetD发作间期和SetE发作期作为异常癫痫类。对每个数据文件进行滑动窗口分割并为每个窗口提取上一节所述的特征同时打上标签。import os from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import classification_report, confusion_matrix, roc_auc_score # 假设我们已经有一个函数 get_features_and_labels(data_folder)它返回特征DataFrame X 和标签数组 y X, y get_features_and_labels(./data/bonn_processed/) # 划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42, stratifyy) # 标准化特征对于基于距离的模型如SVM和神经网络至关重要 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test) # 使用随机森林进行训练和评估 rf_clf RandomForestClassifier(n_estimators100, max_depth10, random_state42, class_weightbalanced) rf_clf.fit(X_train_scaled, y_train) y_pred rf_clf.predict(X_test_scaled) y_pred_proba rf_clf.predict_proba(X_test_scaled)[:, 1] print(分类报告) print(classification_report(y_test, y_pred, target_names[正常, 癫痫])) print(fROC-AUC分数{roc_auc_score(y_test, y_pred_proba):.3f}) # 查看特征重要性辅助特征选择 feature_importance pd.DataFrame({ feature: X.columns, importance: rf_clf.feature_importances_ }).sort_values(importance, ascendingFalse) print(\nTop 10重要特征) print(feature_importance.head(10))随机森林给出的特征重要性排名能让我们直观地看到哪些特征对区分癫痫活动贡献最大。可能是gamma_relative_power、line_length或sample_entropy。基于此我们可以进行递归特征消除或选择重要性高于阈值的特征构建一个更精简、更不易过拟合的模型。对于更复杂的任务如发作期检测在长程EEG中定位发作起始点问题变成了时序上的二分类。我们可以使用滑动窗口提取特征但模型需要能够捕捉时序依赖关系。这时循环神经网络RNN或一维卷积神经网络1D-CNN更为合适。1D-CNN能够自动学习局部时间模式如棘波形状而无需手动设计时域特征。# 使用Keras构建一个简单的1D-CNN模型示例用于窗口分类 import tensorflow as tf from tensorflow import keras from tensorflow.keras import layers # 假设输入是原始信号窗口形状为 (样本数, 时间步长, 通道数) # 这里以单通道256Hz下2秒窗口为例 (512个时间点) model keras.Sequential([ layers.Input(shape(512, 1)), layers.Conv1D(filters32, kernel_size5, activationrelu, paddingsame), layers.MaxPooling1D(pool_size2), layers.Conv1D(filters64, kernel_size3, activationrelu, paddingsame), layers.GlobalAveragePooling1D(), # 替代Flatten参数更少 layers.Dense(32, activationrelu), layers.Dropout(0.3), layers.Dense(1, activationsigmoid) ]) model.compile(optimizeradam, lossbinary_crossentropy, metrics[accuracy]) model.summary() # 训练模型需要准备相应格式的数据 # history model.fit(X_train_3d, y_train, validation_split0.2, epochs30, batch_size32)提示深度学习模型需要大量的数据才能避免过拟合。在EEG领域数据增强技术如添加轻微的高斯噪声、随机缩放、时间扭曲等是增加数据多样性的有效手段。5. 结果可视化与解释让模型决策变得透明模型训练好后我们不能只满足于一个准确率数字。尤其是在医疗辅助诊断领域模型的可解释性至关重要。我们需要知道模型是依据什么做出“癫痫”判断的。对于随机森林或逻辑回归等模型特征重要性本身是一种全局解释。我们还可以使用SHAP或LIME等工具对单个预测样本进行解释展示每个特征对该次预测的贡献度。# 使用SHAP解释随机森林模型安装pip install shap import shap # 创建一个SHAP解释器 explainer shap.TreeExplainer(rf_clf) shap_values explainer.shap_values(X_test_scaled) # 可视化单个样本的解释 sample_idx 0 shap.force_plot(explainer.expected_value[1], shap_values[1][sample_idx, :], X_test.iloc[sample_idx, :]) # 汇总图显示所有样本的特征影响 shap.summary_plot(shap_values[1], X_test, plot_typebar)对于CNN模型可视化其学到的特征也很有意义。我们可以通过梯度加权类激活映射Grad-CAM的变体适用于1D信号来查看在做出“癫痫”判断时模型更关注输入信号窗口中的哪些时间点。# 这是一个概念性代码展示如何为1D-CNN实现简单的激活可视化 def visualize_cnn_activation(model, sample_signal): 为单个样本生成卷积层的激活图。 # 获取目标卷积层的输出和模型输出 conv_layer_model keras.Model(inputsmodel.inputs, outputs[model.layers[-4].output, model.output]) with tf.GradientTape() as tape: conv_output, predictions conv_layer_model(tf.expand_dims(sample_signal, axis0)) loss predictions[:, 0] # 假设是二分类取正类的输出 # 计算梯度 grads tape.gradient(loss, conv_output)[0] # 计算每个特征图的权重全局平均池化梯度 weights tf.reduce_mean(grads, axis0) # 生成类激活图 cam tf.reduce_sum(weights * conv_output[0], axis-1).numpy() # 上采样CAM到原始信号长度如果需要 # ... return cam # 然后可以将CAM叠加在原始信号图上看哪些区域被高亮最后将模型预测结果与原始EEG信号在时间轴上对齐展示是最直观的评估方式。我们可以绘制一段长程EEG信号并在其下方以颜色条或标记点的形式显示模型输出的概率或分类结果直观地观察模型是否在癫痫发作开始和结束时做出了正确的反应。我在处理一个长程头皮EEG项目时就发现模型有时会将强烈的肌电伪迹误判为癫痫发作。通过这种对齐可视化我们迅速定位了问题并通过在预处理中加强肌电伪迹的剔除以及增加对肌电伪迹有区分度的特征如特定高频段的能量比显著提升了模型的鲁棒性。这个过程让我深刻体会到EEG分析不仅仅是调包和跑模型更是一个需要结合信号处理先验知识、不断观察、迭代和调试的综合性工程。