数学建模社区-数学中国

标题: C题相关资源脑电信号分析判别 [打印本页]

作者: 1047521767 时间: 2021-10-15 11:28
标题: C题相关资源脑电信号分析判别
脑电信号按其产生的方式可分为诱发脑电信号和自发脑电信号。诱发脑电信号是通过某种外界刺激使大脑产生电位变化从而形成的脑电活动；自发脑电信号是指在没有外界特殊刺激下，大脑自发产生的脑电活动。

P300电信号：是诱发脑电信号的一种，在小概率刺激发生后300毫秒范围左右出现的一个正向波峰，如下图所示。

[color=rgba(0, 0, 0, 0.5)]它不受刺激物理特性影响，与知觉或认知心理活动有关，与注意、记忆、智能等加工过程密切相关。被测者观察字符矩阵(图1-a)，并记下左上角的字符，然后按顺序闪烁字符矩阵行列(图1-b)的同时，记录被测者的P300脑电信号。

[color=rgba(0, 0, 0, 0.5)]睡眠脑电：是自发脑电信号的一种，能够反映身体状态的自身变化，也是用来诊断和治疗相关疾病的重要依据。根据睡眠状态进一步分为清醒期、睡眠I期、睡眠II期、深度睡眠和快速眼动期，脑电信号如图2所示。

[color=rgba(0, 0, 0, 0.5)]P300脑电数据由5位被测者的20次视觉刺激实验的脑电数据组成，每次实验的字符循环闪烁5轮，记录头顶20个位置的脑电信号，脑电信号记录位置如下图所示。

[color=rgba(0, 0, 0, 0.5)]睡眠脑电数据3000各睡眠脑电样本，由状态标签、Alpha、Beta、Theta和Delta波在“8-13Hz”，“14-25Hz”，“4-7Hz”和“0.5-4Hz”频率范围内的能量占比。
本文主要解决以下3个问题：
1、对P300脑电数据去除冗余信息降维。
2、对不同视觉刺激下的脑电波进行分类，判断闪烁的行或列中是否包含目标字符，最终确定目标字符。
3、对睡眠脑电信号进行分类，根据不同波对应频率范围内的能量比，确定睡眠状态。

数据预处理
建立数据字典
P300脑机接口数据包含训练用的12个已知字符的脑电数据，和测试用的8个未知字符的脑电数据。将训练数据和测试数据，以字符为键脑电数据为值，组成字典，如下所示。
字典 “Sx_train_data_dict” “Sx_test_data_dict” “Sx_train_event_dict” “Sx_test_event_dict”
键目标字符待测字符目标字符待测字符
值字符对应的P300数据表字符对应的P300数据表字符对应的事件标签表字符对应的事件标签表

生成数据索引表计算索引

每个数据表首行首列为字符序号num(见图2-b，从左上到右下依次为1~36)，计算目标字符在矩阵上的行列位置row和col。

row=num//6+1

c o l = n u m % 7 + 1 col = num \% 7 +1
col=num%7+1
创建索引表

针对训练集和测试集的不同字符对应的事件标签表，去除标识符和结束标签所在的行，生成与P300数据表对应的索引表，STrait(t)表示一次闪烁刺激记录的脑电数据：第1列：字符矩阵某次STrait(t)发生位置，记录行号或列号。

第2列：字符矩阵某次STrait(t)开始时对应P300数据表的位置索引

第3列：字符矩阵某次STrait(t)结束时对应P300数据表的位置索引

第4列：字符矩阵某次STrait(t)记录的样本数

第5列：1（闪烁中包含目标字符）或0（闪烁中不包含目标字符）

数据规整

由于每次STrait(t)的P300数据采样情况不同，导致每次采样数据量存在波动。对所有STrait(t)的数据采样量进行统计分析，采样量最小值为3，最大值为46，平均数为39，下四分位数为41，中位数为41，上四分位数为43。因此，删除采样数据量小于40的索引记录。

数据切片

索引表中每行数据记录了每个STrait(t)的闪烁的行/列、开始位置索引、结束位置索引、采样样本量和包含目标字符标志。逐行遍历索引表，将之对应的P300数据表切片，获得每次STrait(t)的P300数据。

问题一：数据去冗余

每个被试者的脑电信号，对于每个目标字符闪烁事件记录的脑电信号共20个通道，大约3000多条数据。原始脑电数据量较大，数据已经在正负类比例上产生较大差距，正例较少导致训练结果一般。另一方面，原数据包含的信息过多，部分通道间的信息有较强相关性，造成数据冗余，对分类产生一定的干扰。

度量指标

皮尔逊相关系数

皮尔逊相关系数是最早由统计学家卡尔·皮尔逊设计的统计指标，是研究变量之间线性相关程度的量，一般用字母 r 表示。

相关系数是按积差方法计算，同样以两变量与各自平均值的离差为基础，通过两个离差相乘来反映两变量之间相关程度；着重研究线性的单相关系数。依据相关现象之间的不同特征，其统计指标的名称有所不同，将反映两变量间线性相关关系的统计指标称为相关系数。

两变量间的皮尔逊相关系数可通过以下公式计算：

Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y方差。对于变量X和Y，当相关系数为0时，X和Y两变量无关系；当X的值增大（减小），Y值增大（减小），两个变量为正相关，相关系数在0.00与1.00之间；当X的值增大（减小），Y值减小（增大），两个变量为负相关，相关系数在-1.00与0.00之间。相关系数的绝对值越大，相关性越强，相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。通常情况下通过以下取值范围判断变量的相关强度：0.8-1.0 为极强相关；0.6-0.8为强相关；0.4-0.6为中等程度相关；0.2-0.4为弱相关；0.0-0.2为极弱相关或无相关。

皮尔逊相关系数实现过程
对于每个表格循环读取每行脑电数据，并保存为最终二维数组。接着将每个通道作为一组向量，求取两两向量间的皮尔逊相关系数，生成大小为[20,20]的相关系数矩阵，可将其存为表格以便观察，如图5所示，是第5位受试者训练数据相关系数矩阵。行列索引表示通道标识符，为该通道的数据向量，矩阵具有沿对角线的对称性，所以只观察下三角区即可。图示中的数据说明，不同通道向量之间的相关性不同，有正相关有负相关，有相关性较强的，也有无关的。相关系数越高，说明这两个通道的脑波的变化特征越相似，相关性极强的可只用一个通道表征这种特征，另一个通道是冗余的。

主成分分析
主成分分析（Principal Component Analysis，PCA），是一种统计方法。通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的综合变量，转换后的这组变量叫主成分，主成分数量较原变量尽量少。在研究多变量任务或者做数据分析处理时，变量个数太多会增加任务复杂性，希望变量个数较少而得到的信息较多。图6为主成分分析原理示意图。

主成分分析实现过程
对于脑电训练数据的数组，与每个通道的均值做差，然后计算协方差矩阵，并求解协方差矩阵的特征值和特征向量。将特征值从大到小进行排序，作为一个特征数组保存，以便之后使用。按照顺序，越靠前特征值越大，包含的信息量就越多，越靠后则信息量越少。

结果融合
对于每位受试者，遍历训练数据中所有目标字符的子表，对其对应数组计算相关系数矩阵和PCA的特征排序数组。对于相关系数矩阵，每位受试者设置不同的阈值，系数大于此值的，将其行列索引即通道标识符进行保存。对于每个子表，对筛选出的相关通道对，比较其在特征排序数组中的前后位置，若比较靠后的通道其索引>10，则说明此通道包含的信息较少，且其信息少于另一通道，则将此通道删除，另一通道保留，依次循环所有通道对。

上述过程，得到特定受试者的12个子表的12个通道组合，将所有组合合成一维数组，计算1-20通道标识符出现的次数，若出现次数>10，则启用此通道，最后得到筛选出的个人最优通道组合。

实验中输出的相关数据如下，图7为受试者1的12个子表PCA特征排序结果，可以看出从13通道以后排列各不相同，说明不同目标字符实验中的脑波也有一定的变化，与实际情况相符。根据受试者生理特征的独特性，相关系数的选择阈值设置也不同。阈值越高，说明对数据相关性或负相关性的要求越高，达到相关性要求的通道对就越少，自然从特征排序中删除的通道数目减少，反之在组合中留下的通道越多。但是由于选取通道个数有一定限制，我们将各受试者组合通道数控制在14和16之间，因此对于不同受试者选定了不同阈值。或者，根据全局最优通道数的要求，可以调整。

所有被试者通道融合
在得出各受试者的最有通道组合后，采取同理的融合方式，将5个通道组合数组合成一维数组，对通道标识符计数，若大于3次，则启用此通道。基于上述结果，全局最优通道组合标识符为下图所示，共包含13个通道，对应的通道名称为Fz、F3、F4、Cz、C3、C4、T7、T8、CP3、CP4、CP5、Pz、P3。有了此结果，可以通过给脑电数据降维，删除冗余通道，留下对应的主特征，进一步提高分类准确率。

问题二：P300脑电信号分类
实验流程
以经过去冗余处理后选取得13个通道数据作为分类输入，通过多种模型得性能对比，选取其中最优的模型，实验流程如下图所示。

数据去噪——滤波
脑电信号时非平稳得，对分类视觉刺激的准确率有很大影响。为提高模型识别准确率与稳定性，我们考虑对脑电信号进行滤波处理。经验模态分解(empirical mode decomposition，EMD) 算法[1]是一种将随机信号分解成多个内涵模态分量（IMF），将非平稳信号转化为平稳信号的算法。在该问题中，我们对给出的非平稳脑电信号序列，依次找出其极大值点与极小值点，通过样条插值，拟合出上下包络曲线，根据上下包络曲线求出平均包络曲线，将原信号数据减去该平均包络曲线便得出分解的第一个IMF分量，将上述过程进行N次迭代，得到N个IMF分量，直至拟合曲线无负局部极小值与正局部极大值，停止迭代，公式如下所示：

其中，x(t)代表原非稳定随机信号，Ci (t)表示经过第i次EMD分解之后的IMF分量，rn(t)代表原信号被分解成一系列IMF分量后的剩余部分。
在解决问题1时，我们对附件1中的信号序列进行EMD分解，结果如下图所示，其中imf1~7为从原始序列中分离出来的频率依次从高至低的局部特征信号，res为经过7次imf分量分离后的剩余信号分量。为减少高低频噪声对后续模型训练的影响，我们将原始信号分别剔除高低频分量作为后续的训练数据。

模型效果比较

输入预处理完毕的Train，Test数据集。并将其划分为10等分，训练集占2成，轮次为4。本文采用共计6个模型，进行分类识别比较，实验如下所示：

通过实验可以发现，具有较好精度性能的模型为Random Forest、Decision TreeGini、fastDTW-KNN。其中Decision TreeGini计算时间最短，但测试集精度最低；SVM算法训练集精度良好，但测试集精度较低且耗时较长；fastDTW-KNN性能相对稳定。

问题三：睡眠脑电数据分类
提供3000个睡眠脑电特征样本及其标签数据，第一列为已知标签，第二至五列为从原始时序列中计算得到的特征参数，依次包括“Alpha”，“Beta”，“Theta”，“Delta”，分别对应了脑电信号在“8-13Hz”，“14-25Hz”，“4-7Hz”和“0.5-4Hz”频率范围内的能量占比，特征参数单位为百分比。其中，清醒期共记录633条，快速眼动期共599条，睡眠期（I期）共记录512条，睡眠期（II期）共记录604条，深睡眠期共记录602条。

本问题旨在使用尽可能少的训练样本数据的基础上实现相对较高的预测准确率，为此，我们将训练集数据占比从0.7、0.5、0.3、0.1，依次输入模型训练并测试，最终确定数据占比和模型选择。

模型性能对比
候选实验模型在不同数据集比例下的准确率如下所示。

通过对比实验可知，在不断减少训练集占比的情况下，两候选模型普遍出现过拟合现象。通过模型对验证集和测试集分类的准确率结果中可知，随机森林模型相对于SVM具有更加稳定的预测性能。

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)