简介：针对2020年研究生数学建模竞赛C题「脑电波分析」的备赛团队与个人，这是一份集代码、数据与文档于一体的完整资料包，围绕面向康复工程的脑电信号分析和判别模型展开，覆盖数据预处理、特征提取、建模与结果输出等常用环节。包内共263个文件，大小118.59MB，以Python源码及编译产物（.py/.pyc）为主，配合25个Excel表格、20个XML配置、18个TXT说明和大量PNG可视化图表，另有Word版题目分析文档，便于按模块查阅。资源涉及P300脑机接口等典型应用场景，并包含多组实验数据与训练脚本，可支撑从信号导入、特征计算到模型评估的完整实践。目前已有43人学习下载。对于备赛团队或个人，可直接参考或复用其中的算法脚本与数据组织方式，结合图表和表格理解脑电信号处理与分类建模流程，有效减少从零搭建代码框架的时间。
1. 2020年研究生数学建模C题脑电波分析：一份数据包背后的完整技术链路
2020年研究生数学建模竞赛C题给出的附件1-P300脑机接口数据，表面上是几个矩阵和一个docx，实际却是一整套脑电信号分析流程的浓缩样本。很多参赛队拿到数据后，最先想到的是跑一个现成的SVM，但忽略了从原始EEG到训练特征之间还有滤波、分段、去伪迹、特征降维这一长串工序。任何一个参数选得不对，最终准确率都会掉到随机水平。这篇文章直接拆解这个题目的数据处理链路，用mne和scikit-learn把每一步怎么落地、参数怎么定、坑在哪里讲清楚，适合准备数模竞赛的团队，也适合刚接触EEG分类的算法工程师。
2. P300脑电数据的读取、滤波与epoch切割：预处理全流程
2.1 附件里的数据结构和采样率确认
解压zip后，除了docx题目描述，最常见的是以eml或mat格式存储的EEG原始数据，也可能是已经按刺激编号切好的npy或csv文件。无论哪种形式，第一步都是确认数据格式和采样率。以numpy矩阵为例，通常形状是(n_times, n_channels)，其中每一行是一个采样点，列对应Fz、Cz、Pz、Oz等电极。如果文件中没有明确标注采样率，可以从相邻时间戳的差里估算：假设时间戳单位是秒，采样间隔为0.004s，那么采样率就是250Hz。这个数字直接决定了后面的滤波截止频率和特征窗口长度。
我一般会用pandas读入原始矩阵，打印shape和前五行，同时检查是否包含NaN。脑电数据偶尔会因为电极接触不良出现整列全为0的通道，这时候需要直接剔除该通道，否则后续ICA和分类特征都会引入无效信息。确认结构之后，把原始矩阵转型成mne能处理的连续信号对象。
2.2 从numpy矩阵构造mne的Raw对象
mne是Python生态中专门处理脑磁图/脑电图数据的库，内置了滤波、重参考、分段、独立成分分析等功能。即使原始数据不是标准的EDF/BDF格式，也可以手动构造Raw对象，代码很直观：
import numpy as np import mne # data: 形状 (n_times, n_channels) sfreq = 250.0 # 第2.1节确认的采样率 ch_names = ['Fz', 'Cz', 'Pz', 'Oz'][:data.shape[1]] # 根据实际通道数截取 info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types='eeg') raw = mne.io.RawArray(data.T, info) # 0.5~30Hz带通滤波，消除基线漂移和高频噪声 raw_filtered = raw.copy().filter(l_freq=0.5, h_freq=30.0, fir_design='firwin', verbose=False)这段代码把原始数组包装成了连续记录的Raw对象。RawArray要求输入形状是(n_channels, n_times)，所以要转置。带通下限0.5Hz用于去除电极缓慢漂移；上限30Hz是因为P300波的频率能量集中在低频段，30Hz以上的肌电噪声需要过滤掉。fir_design='firwin'指定窗函数法设计有限冲激响应滤波器，这种滤波器具有线性相位，不会像IIR滤波器那样导致波形相位扭曲，对后续P300的峰值潜时计算很重要。
如果数据采集时没有做好屏蔽，50Hz工频干扰会在频谱上形成尖峰，还要再加一个陷波滤波器。中国电网工频是50Hz，美国是60Hz，具体看数据来源。raw.notch_filter(freqs=[50.0])可以去除该频率附近的能量。注意如果带通上限正好高于工频，陷波是必须的；如果数据采集系统内置了硬件滤波，也可以省略这一步。
2.3 分段、基线校正、伪迹剔除与验证
P300实验通常是连续呈现一串刺激，每个刺激有一个时间戳和对应的目标/非目标标签。预处理的核心是把连续信号切成一个个以刺激为锚点的epoch。窗口一般取刺激前0.2秒到刺激后0.8秒，刺激前0.2秒作为基线。具体操作如下：
events = np.array([(t0, 0, code) for t0, code in zip(event_times, event_codes)]) epochs = mne.Epochs(raw_notch, events, tmin=-0.2, tmax=0.8, baseline=(-0.2, 0.0), preload=True, event_id={'target': 1, 'nontarget': 2}, picks='eeg', verbose=False) epochs.drop_bad(reject=dict(eeg=100e-6))baseline=(-0.2, 0.0)表示用刺激前200ms的平均幅值作为零点，排除直流漂移。event_id将原始事件码映射为可读标签。reject里设置的100e-6单位是伏特，也就是100微伏，超过这个幅值的epoch大概率包含眨眼、头动等伪迹，直接丢弃。这个阈值并非固定，如果丢弃比例超过30%，说明阈值太严格，或某个通道噪声过大，需要检查通道波形。
预处理完成后，立刻用叠加平均法验证P300是否真的存在。两种trial分别平均，画在同一个坐标系里：
epochs['target'].average().plot(time_unit='s', spatial_colors=True) epochs['nontarget'].average().plot(time_unit='s', spatial_colors=True)如果目标刺激的平均波形在300-500ms处明显比非目标更正，说明数据和预处理都没问题。如果波形高噪声或两条曲线重叠，回头检查滤波参数和分段窗口。这个验证步骤几乎不用花时间，却能够提前发现80%的预处理错误。
下面的表格汇总了预处理阶段最关键的参数以及调整原则，方便直接抄作业。
| 参数 | 常见取值 | 调整依据 |
|---|---|---|
| 带通下限 | 0.1~1Hz | 如果基线漂移严重，设为0.5Hz以上 |
| 带通上限 | 30~40Hz | 上限越高噪声越多；只关注P300就用30Hz |
| 陷波频率 | 50Hz或60Hz | 查看频谱峰位置 |
| epoch窗口 | -0.2 ~ 0.8s | 要覆盖刺激前基线和P300晚期成分 |
| 伪迹阈值 | 100μV | 按通道绝对值统计，放到5%~20%的trial |
| 基线校正 | 刺激前0.2s | 若基线噪声明显，可用全epoch平均做校正 |
3. 特征工程：从时域窗口到小波分解的完整实现
3.1 为什么不能把原始时间点直接送入分类器
一个epoch如果取1秒、4个通道，在250Hz采样率下就有1000个原始数值点。直接用这些点作为SVM输入，特征维度是样本数的几十倍，分类器很容易把噪声记住，在验证集上产生过拟合。此外，P300信号只占据一段特定时间窗，其它时间点基本属于无关背景，特征工程的作用就是把这部分有用信息提取出来，压缩成几十维的小向量。
数学建模竞赛还要求方法可解释，评审希望看到你用了什么特征、为什么这个特征有效。所以特征工程不仅是技术手段，也是最终报告的核心素材。
3.2 时域特征提取函数与参数解析
下面这个函数提取了三个层面的时域特征：整个epoch的统计量、P300时间窗的平均值、窗口内峰值和潜时。输入epochs_data的形状是(n_epochs, n_channels, n_times)，输出是训练矩阵。
def extract_time_features(epochs_data, sfreq=250.0): """ 从epoch数据中提取时域特征 返回 (n_epochs, n_features) """ t = np.arange(epochs_data.shape[2]) / sfreq # 时间轴 n_epochs, n_ch, n_times = epochs_data.shape # 全epoch统计量 mean_all = epochs_data.mean(axis=2) # (n, ch) var_all = epochs_data.var(axis=2) # 反映总能量 ptp_all = np.ptp(epochs_data, axis=2) # 最大最小值差 # 关键时间窗 win1 = (t >= 0.20) & (t <= 0.30) # 早期成分 win2 = (t >= 0.30) & (t <= 0.50) # P300主窗 win3 = (t >= 0.50) & (t <= 0.70) # 晚期正成分 seg2 = epochs_data[:, :, win2] feat_p300 = seg2.mean(axis=2) peak_val = seg2.max(axis=2) peak_lat = seg2.argmax(axis=2) / sfreq + 0.30 w1 = epochs_data[:, :, win1].mean(axis=2) w3 = epochs_data[:, :, win3].mean(axis=2) feats = np.hstack([mean_all, var_all, ptp_all, w1, feat_p300, w3, peak_val, peak_lat]) return feats特征维度等于n_channels * 8：整体均值、方差、峰峰值、早期窗、P300主窗、晚期窗、峰值幅值、峰潜时。这里的feat_p300是P300窗口的平均幅值，是分类贡献最大的特征。peak_lat代表该通道在P300窗口内达到最大值的时刻，目标刺激的P300潜时通常稳定在300-450ms之间，而非目标刺激的最大值往往更早，所以这个特征也很有区分度。
注意np.hstack是将所有通道的特征按顺序排成一行。如果四个通道，最终得到32维特征。这个维度对数百个trial的数据集来说是安全的，不至于过拟合，同时保留了足够的判别信息。
3.3 频域能量特征：小波分解的用法与收益
时域特征容易受基线漂移影响，频域特征则更稳定。P300在频域上表现为低频段能量增强，尤其是delta波段（0.5-4Hz）和theta波段（4-7Hz）。离散小波变换（DWT）可以在不同频带上给出时频能量，我常用pywt这个库。
import pywt def wavelet_energy(epoch_signal, wavelet='db4', level=4): """ epoch_signal: 单个通道、一维时间序列 返回: 各层归一化能量 """ coeffs = pywt.wavedec(epoch_signal, wavelet, level=level) energy = np.array([np.sum(np.square(c)) for c in coeffs]) return energy / np.sum(energy)db4小波形状与P300波形接近，分解后第3层和第4层系数对应低频分量。归一化能量消除了个体幅值差异，让不同受试者的特征可比。实际操作中，我会将每个通道的4个能量值追加到时域特征后面，使单通道特征数变为12，总数达到48维。加入能量特征后，SVM和LDA通常会有2~4%的AUC提升。代价是计算时间变长，因为wavedec需要对每个epoch每个通道分别跑一遍。如果trial数量超过5万，可以先截取0.2~0.7s片段再分解，能省掉将近一半的时间。
3.4 特征选择与PCA降维的取舍
特征并非越多越好。下表列出常用特征维度以及使用建议，帮助快速决定保留哪些。
| 特征类型 | 特征维度 | 优势 | 劣势 | 建议 |
|---|---|---|---|---|
| 全时间点原始值 | 200*通道 | 信息最完整 | 维度高、噪声大 | 不推荐 |
| 时域统计量 | 3*通道 | 计算快、稳定 | 丢失波形细节 | 保留均值与方差 |
| 分窗均值 | 3*通道 | 对齐P300时间 | 窗口选择敏感性高 | 核心特征 |
| 峰值与潜时 | 2*通道 | 直接捕获P300 | 受滤波影响大 | 保留 |
| 小波能量 | 4*通道 | 频域鲁棒性好 | 计算较慢 | 视计算资源 |
| PCA降维结果 | 10-20维 | 去冗余 | 可解释性下降 | 与特征联合使用 |
在特征高度相关时，PCA可以压缩维度。但PCA必须在训练集上拟合，然后用这个拟合好的变换去转换测试集，而不能在整体数据上先做PCA再切分，否则会把测试集分布信息泄漏到训练过程中，造成验证精度虚高。如果你只是想提高模型分数，我更推荐先保留解释性强的原始特征，把PCA交给分类器前的Pipeline处理，而不是预先做一次全局降维。
4. 分类建模：LDA基线与RBF-SVM调优
4.1 标签构造和类别不平衡问题
P300实验的原始事件码通常有两类：目标刺激和非目标刺激。目标刺激占比一般在20%左右，属于典型的类别不平衡问题。如果直接用准确率评估模型，即使全部预测为非目标，也能得到80%的准确率，但这显然没有意义。因此标签构造时，需要将目标类指定为正类，评估指标以ROC-AUC或F1为主，而不是准确率。
另外，要检查event_id是否与真实含义对应。有些数据包的标签用1表示目标、0表示非目标，有些用1和2，还有的用字符串。构造标签向量时务必输出前10个样本的人工核对。
4.2 LDA基线：协方差收缩的稳定性论证
线性判别分析是P300分类的传统方法，优势在于无需调参、训练快。当特征维度接近样本数时，协方差矩阵估计会不稳定，所以建议使用收缩估计器：
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.3, stratify=y, random_state=42 ) lda = LinearDiscriminantAnalysis(solver='lsqr', shrinkage='auto') lda.fit(X_train, y_train) pred = lda.predict(X_test) print(classification_report(y_test, pred, target_names=['NonTarget', 'Target']))solver='lsqr'表示用最小二乘法求解，该实现配合shrinkage='auto'可以自动估计收缩强度，即使特征数量接近样本量也不会出现奇异矩阵。默认的svd求解器虽然不需要计算协方差，但无法使用收缩机制。stratify=y让训练集和测试集的正负样本比例与原始数据一致，避免切分后少数类别在一侧缺失。
LDA作为基线的意义重大：如果LDA的AUC已经很接近1，说明特征工程已经足够好，没有必要上复杂的深度学习模型。如果LDA的AUC低于0.7，问题大概率在特征或预处理，而不是分类器。
4.3 使用Pipeline + GridSearchCV做SVM
SVM的RBF核可以学习非线性边界，在P300数据上通常能比LDA高几个百分点。但RBF核有两个超参数需要调：正则化系数C和核宽度gamma。正确的调参方式是把标准化、分类器、网格搜索串成一个Pipeline：
from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline pipe = Pipeline([ ('scaler', StandardScaler()), ('svm', SVC(kernel='rbf', class_weight='balanced', probability=True)) ]) param_grid = { 'svm__C': [0.1, 1, 10], 'svm__gamma': [0.001, 0.01, 0.1] } grid = GridSearchCV(pipe, param_grid, cv=5, scoring='roc_auc', n_jobs=-1) grid.fit(X_train, y_train) print("Best C, gamma:", grid.best_params_)StandardScaler至关重要，因为特征中既有微伏级别的幅值，又有秒级别的潜时，量纲差异会导致核距离计算被大数值特征支配。class_weight='balanced'会自动加大少数类的错误惩罚，缓解目标刺激比例过低的问题。网格搜索的候选集合不要设得太大，三个C值、三个gamma值足够在AUC上找到优质区域。如果训练很多轮后分数仍无改善，改回LDA可能是更理性的选择。
GridSearchCV内置的5折交叉验证完全在训练集内部执行，网格搜索每次只接触训练数据，测试集始终保留到最终评估，不存在调参泄漏。
4.4 混淆矩阵与误判分析
分类报告和数据即将结束时，务必输出混淆矩阵用于后续诊断：
from sklearn.metrics import confusion_matrix best_svm = grid.best_estimator_ pred_test = best_svm.predict(X_test) cm = confusion_matrix(y_test, pred_test) print(cm) # [[TN, FP], [FN, TP]]假如目标刺激被判为非目标的次数很多（假阴性偏高），说明P300特征强度不够，需要回到特征提取阶段，检查P300窗口是否应该从300ms改为250-450ms，或者增加单trial的小波能量特征。如果非目标被误判为目标（假阳性偏高），一般是残留眨眼伪迹带有很大的正向偏转，可以尝试将伪迹阈值调严到80μV，或者加入ICA去除眼电成分后再重新提取特征。
下面这张表总结了SVM调参时常见参数组合及其表现倾向，供实际训练时对照。
| C | gamma | 常见表现 | 使用时机 |
|---|---|---|---|
| 0.1 | 0.001 | 欠拟合，AUC偏低 | 特征噪声大 |
| 1 | 0.01 | 平衡 | 默认起始点 |
| 10 | 0.1 | 过拟合，训练集AUC极高 | 特征很干净、样本量大 |
5. 进阶技巧：交叉验证、抵抗数据泄露与报告呈现
5.1 分层K折与GroupKFold的使用时机
单次train_test_split的结果受随机种子影响太大，不能代表模型真实水平。更可靠的评估是分层K折交叉验证，让每一折的正负比例保持一致。如果数据中存在一个受试者多段连续记录，还要用GroupKFold按受试者或block分组，防止同一受试者的数据同时出现在训练和验证集，导致分数虚高。实际使用中，我习惯在外层用StratifiedKFold固定5折，把SVM的超参数事先通过内部网格搜索确定，再用这5折报告均值和方差。
5.2 数据泄露的三种常见来源
数据泄露是数模竞赛分数崩盘的最常见原因。第一，在完整数据上拟合标准化或PCA，然后再划分训练测试集，导致测试集信息参与了训练。正确做法是把标准化和PCA放进Pipeline，让它只在训练集上估计参数。第二，ICA去伪迹时用全部数据估计独立成分，然后再划分数据。ICA分解矩阵本身带有全局分布信息，应该先用训练集估计，应用到测试集时只是做矩阵变换。第三，用滑动窗口生成样本时，相邻窗口有重叠，随机切分会把同一个原始信号片段分到两边，需要用GroupKFold按刺激编号分组。
5.3 在竞赛报告中如何呈现模型结果
报告不要只给一个最终准确率。评审更看重特征与脑电生理机制的可解释性。建议在方法部分明确写出“采用Pz通道，300-500ms窗口的平均幅值作为特征之一，因为P300在此区域最突出”，并在结果部分画ROC曲线，标注AUC值。再加上一句“目标刺激与非目标刺激在该特征上差异显著（配对t检验，p<0.01）”，整个模型的合理性一下就立住了。
最后一个实用技巧：把数据读取、滤波、特征提取、分类分别封装成独立脚本，并在报告附录列出运行顺序。这样即使评审想复现，也能一步步执行，而不是在一大段代码里找入口。这个习惯对于比赛和后续项目开发都同样重要。
本文还有配套的精品资源，点击获取
