2020年全国研究生数学建模竞赛华为杯C题面向康复工程的脑电信号分析和判别模型求解全过程文档及程序

2020年全国研究生数学建模竞赛华为杯

C题 面向康复工程的脑电信号分析和判别模型

原题再现：

  背景和意义
  大脑是人体中高级神经活动的中枢，拥有着数以亿计的神经元，并通过相互连接来传递和处理人体信息。脑电信号按其产生的方式可分为诱发脑电信号和自发脑电信号。诱发脑电信号是通过某种外界刺激使大脑产生电位变化从而形成的脑电活动；自发脑电信号是指在没有外界特殊刺激下，大脑自发产生的脑电活动。
  （1）诱发脑电信号（P300脑-机接口）
  在日常生活中，人的大脑控制着感知、思维、运动及语言等功能，且以外围神经为媒介向身体各部分发出指令。因此，当外围神经受损或肌肉受损时，大脑发出指令的传输通路便会受阻，人体将无法正常完成大脑指令的输出，也就失去了与外界交流和控制的能力。研究发现，在外围神经失去作用的情况下，人的大脑依旧可以正常运行，而且其发出指令的部分信息可以通过一些路径表征出来。脑-机接口技术旨在不依赖正常的由外围神经或肌肉组织组成的输出通路的通讯系统，实现大脑与外部辅助设备之间的交流沟通。
  P300事件相关电位是诱发脑电信号的一种，在小概率刺激发生后300毫秒范围左右出现的一个正向波峰（相对基线来说呈现向上趋势的波）。由于个体间的差异性，P300的发生时间也有所不同，图1表示的是在刺激发生后450毫秒左右的P300波形。P300电位作为一种内源性成分，它不受刺激物理特性影响，与知觉或认知心理活动有关，与注意、记忆、智能等加工过程密切相关。基于P300的脑-机接口优点是使用者无需通过复杂训练就可以获得较高的识别准确率，具有稳定的锁时性和高时间精度特性。

  （2）自发脑电信号（睡眠脑电）
  睡眠是身体休整积蓄能量的重要环节，睡眠质量对人的身心状态也有着重大影响。如何提高睡眠质量，减少睡眠相关疾病对健康的影响，日益受到广泛关注。睡眠过程中采集的脑电信号，属于自发型的脑电信号。自发型的睡眠脑电信号能够反映身体状态的自身变化，也是用来诊断和治疗相关疾病的重要依据。
  睡眠过程是一个动态变化的复杂过程。在国际睡眠分期的判读标准R&K中，对睡眠过程中的不同状态给出了划分：除去清醒期以外，睡眠周期是由两种睡眠状态交替循环，分别是非快速眼动期和快速眼动期；在非快速眼动期中，根据睡眠状态由浅入深的逐步变化，又进一步分为睡眠I期，睡眠II期，睡眠III期和睡眠IV期；睡眠III期和睡眠IV期又可合并为深睡眠期。图2给出了不同睡眠分期对应的脑电信号时序列，自上而下依次为清醒期、睡眠I期、睡眠II期、深睡眠和快速眼动期。从图2中可以观察到，脑电信号在不同睡眠分期所呈现的特点有所不同。基于脑电信号进行自动分期，能够减轻专家医师的人工负担，也是评估睡眠质量、诊断和治疗睡眠相关疾病的重要辅助工具。

  课题任务
  本赛题包含2个附件（数据文件），四个课题任务。具体说明如下，
附件1：P300脑机接口实验数据
  提供了5个健康成年被试（S1-S5）的P300脑机接口实验数据，平均年龄为20岁。在实验的过程中，要求每一位被试（被测试者）集中注意力。P300脑机接口实验的设计如下：每位被试能够观察到一个由36个字符组成的字符矩阵，如图3所示，字符矩阵以行或列为单位（共6行6列）。每轮实验的设计流程：首先，提示被试注视“目标字符”，例如在图3的字符矩阵上方，出现的灰色字符“A”；其次，进入字符矩阵的闪烁模式，每次以随机的顺序闪烁字符矩阵的一行或一列，闪烁时长为80毫秒，间隔为80毫秒；最后，当所有行和列均闪烁一次后，则结束一轮实验。在被试注视“目标字符”的过程中，当目标字符所在行或列闪烁时，脑电信号中会出现P300电位；而当其他行和列闪烁时，则不会出现P300电位。上述实验流程为1轮，共重复5轮。

  每位被试的P300脑电数据包含有4个文件，具体说明如下，
  train_data：训练用数据；
  train_event：训练数据的事件标签；
  test_data：测试用数据；
  test_event：测试数据的事件标签。
  训练用数据包括12个已知目标字符的数据（char01-char12），测试用数据包括10个待识别目标字符的数据（char13-char22）。每个字符矩阵闪烁实验中，脑电数据表格包含有20列（每列表示1个记录通道，记录通道依次进行编号，表1为记录通道的标识符，图5对应了记录通道的位置），脑电数据表格的行表示样本点数据，采样频率为250Hz。信号采集设备设置了参考电极和接地电极，即记录通道的信号为作用电极与参考电极之间的差值。

  训练数据中的标签文件同样是以子表形式与实验数据相对应，子表的名称为“charXX(Y)”，XX对应相应字符的序列号，Y表示实际的目标字符。子表的内容包含了两列，第一列表示标签，第二列为采样点序号。每轮实验的起始标签为目标字符对应的标识符（字符矩阵中36个字符的标识符详见表2，如“101”表示“A”），接下来为闪烁的行或列的标识符（详见图6，如“2”表示第2行，“9”表示第3列），一轮实验的结束标签为“100”。在训练数据的事件标签文件中，第一行给出了目标字符的标识符和对应的采样点序号，接下来是随机闪烁的行和列的标识符和对应的采样点序号，每轮实验以“100”标识符结束，共重复5次；
  测试数据中的标签文件同样是以子表形式与实验数据相对应，子表的名称为“charXX”，XX对应相应字符的序列号。在测试数据的事件标签文件中，第一行给出了待识别目标字符的标识符，统一表示为“666”，需要通过对脑电信号进行分析后，得到出现P300电位的行和列，并判断得到目标字符的识别结果。

  附件2：睡眠脑电数据
  提供3000个睡眠脑电特征样本及其标签，取自不同的健康成年人整夜睡眠过程。第一列为“已知标签”，用数字形式来表示不同的睡眠分期：清醒期（6），快速眼动期（5），睡眠I期（4），睡眠II期（3），深睡眠期（2）；第二至五列为从原始时序列中计算得到的特征参数，依次包括“Alpha”，“Beta”，“Theta”，“Delta”，分别对应了脑电信号在“8-13Hz”，“14-25Hz”，“4-7Hz”和“0.5-4Hz”频率范围内的能量占比，特征参数单位为百分比。
  根据以上附件所给出的数据来源和实验数据，请研究以下问题：
  问题一：在脑-机接口系统中既要考虑目标的分类准确率，同时又要保证一定的信息传输速率。请根据附件1所给数据，设计或采用一个方法，在尽可能使用较少轮次（要求轮次数小于等于5）的测试数据的情况下，找出附件1中5个被试测试集中的10个待识别目标，并给出具体的分类识别过程，可与几种方法进行对比，来说明设计方法的合理性。
  问题二：由于采集的原始脑电数据量较大，这样的信号势必包含较多的冗余信息。根据图5和表1，在20个脑电信号采集通道中，无关或冗余的通道数据不仅会增加系统的复杂度，且影响分类识别的准确率和性能。请分析附件1所给数据，并设计一个通道选择算法，给出针对每个被试的、更有利于分类的通道名称组合（要求通道组合的数量小于20大于等于10，每个被试所选的通道可以不相同，具体的通道名称见图5和表1）。基于通道选择的结果，进一步分析对于所有被试都较适用的一组最优通道名称组合，并给出具体分析过程。为了方便参赛者对最优通道组合进行选择，赛题给出了测试数据（char13-char17）的结果，它们的字符分别是：M、F、5、2、I。
  问题三：在P300脑-机接口系统中，往往需要花费很长时间获取有标签样本来训练模型。为了减少训练时间，请根据附件1所给数据，选择适量的样本作为有标签样本，其余训练样本作为无标签样本，在问题二所得一组最优通道组合的基础上，设计一种学习的方法，并利用问题二的测试数据（char13-char17）检验方法的有效性，同时利用所设计的学习方法找出测试集中的其余待识别目标（char18-char22）。
  问题四：根据附件2中所给出的特征样本，请设计一个睡眠分期预测模型，在尽可能少的训练样本的基础上，得到相对较高的预测准确率，给出训练数据和测试数据的选取方式和分配比例，说明具体的分类识别过程，并结合分类性能指标对预测的效果进行分析。
注意：在研究问题一时，不能将问题二提供的5个测试数据的结果作为已知条件，若使用该条件，则被判为作弊处理。

整体求解过程概述(摘要)

  脑电信号的识别和分类是脑机接口技术中非常重要的一环，使用者无需通过复杂训练就可以获得较高的识别准确率，具有稳定的锁时性和高时间精度特性。本文使用基于监督学习的随机森林，SVM算法，半监督学习S3VM等算法，以较高的分类准确度，完成了对P300脑电信号的分析和判别。
  针对问题一，我们设计了0.1~20Hz，阻带增益-80dB的46阶带通巴特沃斯IIR滤波器对冗余数据做预处理，以是否出现P300脑电信号作为输入标签。考虑到发生P300事件和不发生的样本数不平衡，采用随机排列抽取法生成训练集，并建立了基于随机森林算法的监督学习分类模型。在python环境下搭建算法，预测5轮测试数据，采取“投票”机制得到不同受试者的10个测试结果（共48个目标字符），目标字符识别准确率为92%，相比KNN模型84%的准确率，识别效果更好。
  针对问题二，在问题一数据处理的基础上，从闪烁时刻起记录有效训练数据，并生成一个维度为的矩阵，同样地对P300电位做标签化处理，得到一个7500维的向量。将通道选择问题转换成回归模型，进一步引入“spike and slab”先验分布，建立贝叶斯稀疏回归模型。使用python3.8实现EP算法，推导出模型后验分布近似分布，依据支撑向量后验分布参数选择出最优通道组合，基于通道选择结果，再次采用问题一的算法和模型，找到了对所有被试者都适用的一组最优通道组合：Fz, F3, F4, C3, Cz, C4, CP3, CP4, CP5, CP6, P3, P4, P7。
  针对问题三，在问题一数据处理和问题二通道选择的基础上，选择10个字符的数据作为训练集（共个样本），其中前30%作为有标签数据集，后70%作为无标签数据集，剩下作为验证集（共个样本）。建立半监督学习的S3VM分类模型，并引入正负样本比例因子r排除正负样本数目不平衡的问题。模型通过带标记样本预训练分类器，用拉格朗日启发式方法对无标签样本进行标签预测，扩充训练集。通过模拟退火的最优化方法进行迭代，并在训练时对初始无标签样本和有标签样本赋予不同的权重。最终在验证集上的平均准确率达到75.3%，在测试数据中得到的结果基本吻合。
  针对问题四，我们根据附件2所给的精确无噪声的样本标签数据，以及四个分区Alpha，Beta，Theta，Delta 所占能量强度不同的特点，分别建立随机森林，支持向量机（SVM）,以及LightGBM三种监督学习算法训练睡眠预测模型。在测试集和训练集的选取上，采用“交叉验证法”，克服因为训练集测试集不同分割方法引起模型性能不稳定。测试选取不同比例训练集数据对模型预测精度影响。综合比较SVM，随机森林和LGBM算法预测精度和时间复杂度，得出，LGBM算法得到的精度为89%，且算法测试运行时间较短，其综合运行性能最优。
  综上所述，本文建立的脑电信号分析和判别模型地实现了P300脑-机接口对目标字符的识别，最优通道选择，解决了带标签样本较难获得的问题。基于睡眠脑电信号频谱能量不同特点，建立模型实现了对睡眠状态的预测，有望成为评估睡眠质量、诊断和治疗睡眠相关疾病的重要辅助工具。

模型假设：

  （1）假设五名被测试成年人身体健康，无精神疾病史、脑部疾病史且视力正常；
  （2）假设设备确定能够采集到P300电位
  （3）假设两次目标刺激产生的正电位波峰不会重合；
  （4）假设个体间的差异不会超过设定值；

问题分析(部分)：

  从问题的描述可知，在小概率刺激发生后300毫秒范围左右P300 将会出现一个正向的波峰（根据P300的时域特征）。考虑到个体间的差异性，我们首先假设P300的波峰特征出现在目标符号所在行列闪现后的300毫秒左右的可偏差范围内。本题在预处理不同通道采集信号的基础上，结合是否在规定区间内出现P300脑电信号波峰特征为特征向量（分类标准），考虑到训练数据种类的不平衡，需要平衡P300和非P300电位的数据样本，并以此形成训练集并构建决策树，训练随机森林算法。使用算法预测 P300 信号出现的行和列，通过行列位置最终确定出S1,S4,S5的10个和S2,S3的9个待识别字符。
  题目给出的原始的20个通道采集的脑电数据非常庞大，其中参杂着很多冗余信息，不仅会带来系统的复杂性，增大系统数据传输和处理的压力，甚至会造成模型的分类和识别精度的下降，因此我们将剔除影响较小甚至没有影响的通道以选取最优的通道组合。传统的通道选择方法分为两类，一类是经验选择法，另外一种是机器学习法。经验选择法需要扎实的生物学相关的知识和前期大量的实验来论证各个通道的比例权重，而且其选择的算法不能满足每个被试者的特点，会造成重要数据丢失，影响实验结果。机器学习算法则是通过挖掘数据内在的特点，能够根据不同的被试者自动调整各个通道的采集数据的重要性。所以我们选择了机器学习算法来解决通道选择问题。 考虑到脑电信号子集的空间稀疏性，目前已经存在多种稀疏贝叶斯回归方法，比如对所有脑电通道施加范数约束，选择出具有空间稀疏性的最优通道子集LASSO方法，对通道不同组合施加范数约束的组GROUP-LASSO法[9]。本文以“spike and slab”先验建立稀疏贝叶斯回归模型，以此训练出每个被试者的最佳通道组合。

论文缩略图(全部)：

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码(部分 切勿抄袭)：

1. function Hd = bandpass   
2. %BANDPASS Returns a discrete-time filter object.   
3. % Butterworth Bandpass filter designed using FDESIGN.BANDPASS.   
4.    
5. % All frequency values are in Hz.   
6. Fs = 250;  % Sampling Frequency   
7.    
8. Fstop1 = 0.1;         % First Stopband Frequency   
9. Fpass1 = 0.2;         % First Passband Frequency   
10. Fpass2 = 20;          % Second Passband Frequency   
11. Fstop2 = 30;          % Second Stopband Frequency   
12. Astop1 = 80;          % First Stopband Attenuation (dB)   
13. Apass  = 1;           % Passband Ripple (dB)   
14. Astop2 = 80;          % Second Stopband Attenuation (dB)   
15. match  =