2020年华为杯数学建模竞赛C题论文和代码

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

    摘           要：

脑电信号的识别和分类是脑机接口技术中非常重要的一环，使用者无需通过复杂训练就可以获得较高的识别准确率，具有稳定的锁时性和高时间精度特性。本文使用基于监督学习的随机森林，SVM算法，半监督学习S3VM 等算法，以较高的分类准确度，完成了对P300脑电信号的分析和判别。

针对问题一，我们设计了0.1~20Hz，阻带增益-80dB的46阶带通巴特沃斯IIR 滤波器对冗余数据做预处理，以是否出现P300脑电信号作为输入标签。考虑到发生P300事件和不发生的样本数不平衡，采用随机排列抽取法生成训练集，并建立了基于随机森林算法的监督学习分类模型。在python环境下搭建算法，预测5轮测试数据，采取“投票”机制得到不同受试者的10个测试结果(共48个目标字符)， 目标字符识别准确率为92%，相比KNN 模型84%的准确率，识别效果更好。

针对问题二，在问题一数据处理的基础上，从闪烁时刻起记录有效训练数据，并生成一个维度为7500×20的矩阵，同样地对P300 电位做标签化处理，得到一个7500维的向量。将通道选择问题转换成回归模型，进一步引入“ spike and slab”先验分布，建立贝叶斯稀疏回归模型。使用python3.8实现EP算法，推导出模型后验分布近似分布，依据支撑向量后验分布参数选择出最优通道组合，基于通道选择结果，再次采用问题一的算法和模型，找到了对所有被试者都适用的一组最优通道组合：Fz，F3，F4，C3，Cz，C4，CP3,CP4,CP5,CP6,P3,P4, P7。

针对问题三，在问题一数据处理和问题二通道选择的基础上，选择10个字符的数据作为训练集(共12×5×10×5=3000个样本), 其中前 30%作为有标签数据集,后70%作为无标签数据集，剩下作为验证集(共12×5×2×5=600个样本) 。建立半监督学习的S3VM 分类模型，并引入正负样本比例因子r排除正负样本数目不平衡的问题。模型通过带标记样本预训练分类器，用拉格朗日启发式方法对无标签样本进行标签预测，扩充训练集。通过模拟退火的最优化方法进行迭代，并在训练时对初始无标签样本和有标签样本赋予不同的权重。最终在验证集上的平均准确率达到75.3%，在测试数据中得到的结果基本吻合。

针对问题四，我们根据附件2所给的精确无噪声的样本标签数据，以及四个分区 Alpha, Beta, Theta, Delta所占能量强度不同的特点, 分别建立随机森林, 支持向量机(SVM) ，以及LightGBM 三种监督学习算法训练睡眠预测模型。在测试集和训练集的选取

上，采用“交叉验证法”，克服因为训练集测试集不同分割方法引起模型性能不稳定。测试选取不同比例训练集数据对模型预测精度影响。综合比较SVM，随机森林和LGBM算法预测精度和时间复杂度，得出，LGBM算法得到的精度为89%，且算法测试运行时间较短，其综合运行性能最优。

综上所述，本文建立的脑电信号分析和判别模型地实现了P300脑-机接口对目标字符的识别，最优通道选择，解决了带标签样本较难获得的问题。基于睡眠脑电信号频谱能量不同特点，建立模型实现了对睡眠状态的预测，有望成为评估睡眠质量、诊断和治疗睡眠相关疾病的重要辅助工具。

关键词：P300脑机接口； 监督学习； 随机森林算法； 通道选择； 贝叶斯稀疏回归模型；EP; spike and slab; 半监督学习S3VM 算法

即有效信号仅会发生在目标所在行和列闪烁后的125个采样点内，而目标信号闪烁之前人脑并不会产生 P300 信号，因此采样闪烁之前的信号是没有意义的，不仅增加了数据冗余，还会干扰随后预测模型的精度。当然，我们不仅要采集目标字符出现后125个采样点这些有效数据，也需要非目标字符所在行或列闪烁后的125个采样点后的信号，这个数据作为人体不出现 P300 信号时候的正常脑电波信号，以此作为负的训练样本。每个行和列的闪烁时间为80ms，闪烁后间隔80ms进行下一次闪烁，因此数据选取窗口的滑动点数大小应为 250Hz*80+80ms=40,数据的选取如下图所示：

4.2.3 数据标签的制作

预处理的另外一个重要的问题就是标签的制作。使用4.2.1以及4.2.2数据预处理及冗余去除过程后，我们仅仅需要判断当一个行或列闪烁时所采集的信号，判断 P300 信号是否出现，这里P300信号只有两种状态，有和无，所以将脑电信号的输入标签化处理：

区间E中出现波峰特征

区间E中没有波峰特征

其中，P₄，用于判断是否出现P300 电位， E∈ti,jti,j+500,ti,j表示第i个目标字符在第j轮的起始采样时间。

4.2.4基于随机森林的分类问题算法

随机森林( Random Forest，RF) 算法是一种组合分类器算法，包含许多独立的决策树( Decision Tree, DT) [6]。决策树算法以树状结构表示数据分类的结果，是一种典型的数据分类器。决策树是随机森林的核心部分，依据条件熵和信息增益保证决策树的左右生成，其基本原理为[7]：

(1)采用 Bagging的方法从样本数N的数据集D中随机地选取不同变量数，用以形成众多的随机决策树并决定每棵树的分类节点。

(2)计算数据集D的信息熵，以及各个变量Xᵢ和D之间的条件熵。信息熵H(D)体现了数据发生的不确定性，信息熵的值越高，表明越混乱，分类效果不好。条件熵H(D|X₄)是在数据集D发生的前提下，X；变量所带来的熵值变化。用公式可表示为：

4.3模型的求解

4.3.1平衡训练数据种类

使用4.1预处理数据训练模型时，存在一个非常严重的问题。每个被试者进行的每个字符的每轮实验中，每个行和列都会闪烁一次，即一共闪烁12次，但在每轮实验中，出现目标字符的行和列一共只有2个。因此，经过预处理后标签为1的数据，即出现 P300 波形的数据仅有16.7%，其余84%的数据的标签均为0，即没有P300信号。用这样一个标签极其不平衡的数据训练的模型，最后模型必然出现预测结果大部分为0的情况。模型只要预测结果都为0，训练精度就能达到 84%。

为了解决不同标签数据量不平衡问题，我们采用随机排列抽取法来降低标签为0 数据量，保留标签为1数据量的方法。具体做法如下：每个被试者在进行同一个字符的5轮实验中，每轮实验目标字符所在行和列闪烁时所采集的数据，即标签为1的数据保留。剩下的行或列闪烁时，在每轮实验均分别取两个不重复的行或列作为有效训练数据。比如：第一个被试者在注视字符‘B’时，5轮实验中，每次行1或者列8闪烁时的数据保留，然后取其采样点后 125个数据作为有效数据，经过滤波，标签为1。如预处理部分所述，共得到10组标签为1的有效信号。剩下共10个行或者列闪烁时分别在不同组取，如2，3行闪烁时采集的标签为0的信号取与实验1，行4，5闪烁时候，经过预处理后标签为0的数据取于实验 2，以此类推。其余字符闪烁时也做同样的处理。这样每个实验者，在进行每个不同字符的实验时候，均能得到 10个标签为0 和标签为1 的有效训练数据。它们数量相等，不会出现直接将所有实验数据做同样处理后出现的不平衡情况。

另外，除有效字符所在位置外，余下的 10 个行和列在 5 个实验中的取法共有 C102C82⋅⋅⋅C42种方法，并且不同的字符进行不同取法进行组合，会产生大量不同的数据集。为了消除不同取法所带来的随机影响，每次采用不同取法整合出标签0的数据，与标签为1 的数据合成一个有效训练数据训练模型，将每种取法均实验一次显然是不可行，也是没有必要的。我们采取下列这样方法：将实验字符所在行或者列外的其他10个行或列随机重排两次，每次依次从排列的行或列中选择2个，这两个行和列的数据按顺序采集于实验1，2，…。比如1号被试者在进行字符‘B’实验时，将‘B’所在行或列去除，剩余10个行或列的一次随机排列为：

9, 5, 3, 7, 2, 4, 12, 10, 11

这些行或者列闪烁后采集数据标签一定为0，其中列9行5 闪烁时的数据采集于实验1，行3列7采集于实验2，以此类推。这样的取法每个被试者每个数据采集2次。所以每个被试者均可以得到两组有效训练数据。

4.3.2实验结果

基于 Python3.8搭建出随机森林字符预测模型，使用上述的数据处理得到训练数据训练模型，对测试集数据做同样的预处理后，送入训练出的模型进行预测计算。最为理想的数据预测情况应该是每个被试者的每个字符，行和列均出现一个P300信号，取出行和列，即可定位出预测字符。但模型不可能具有 100%精确度，有可能在其他位置也预测出 P300信号，我们采用“少数服从多数”的投票方法进行结果的提取。每个字符的5组测试数据中，模型预测出出现P300信号位置，就给该位置加1，即投给该位置一票，5组预测实验结束后，分别取得票最高的行和列作为最终结果，并且以此定位字符作为预测的最终结果。五个被试者的“投票”表格情况如下：

表 4-1 被测者 S1 P300信号预测位置统计表

表4-2 被测者 S2 P300信号预测位置统计表

表 4-3 被测者 S3 P300信号预测位置统计表

表4-4 被测者 S4 P300信号预测位置统计表

表 4-5 被测者 S5 P300信号预测位置统计表

依据每个字符的5 轮实验中，模型预测的各位置出现 P300 信号次数，分别取出现次数最高的行和列序号，以此定位出需要识别的字符，其结果如下：

表4-6待识别字符统计表

表4-10 随机森林与KNN对比表

5. 问题二的模型建立与求解

5.1问题二的分析

题目给出的原始的20个通道采集的脑电数据非常庞大，其中参杂着很多冗余信息，不仅会带来系统的复杂性，增大系统数据传输和处理的压力，甚至会造成模型的分类和识别精度的下降，因此我们将剔除影响较小甚至没有影响的通道以选取最优的通道组合。传统的通道选择方法分为两类，一类是经验选择法，另外一种是机器学习法。经验选择法需要扎实的生物学相关的知识和前期大量的实验来论证各个通道的比例权重，而且其选择的算法不能满足每个被试者的特点，会造成重要数据丢失，影响实验结果。机器学习算法则是通过挖掘数据内在的特点，能够根据不同的被试者自动调整各个通道的采集数据的重要性。所以我们选择了机器学习算法来解决通道选择问题。

考虑到脑电信号子集的空间稀疏性， 目前已经存在多种稀疏贝叶斯回归方法，比如对所有脑电通道施加L₁范数约束，选择出具有空间稀疏性的最优通道子集 LASSO 方法[8], 对通道不同组合施加L₁范数约束的组GROUP-LASSO法[9]。本文以“ spike and slab”先验建立稀疏贝叶斯回归模型，以此训练出每个被试者的最佳通道组合。

5.2模型建立

不失一般性，我们以第一个被试者为例来推导我们的模型。首先，将某个采样点20个通道采集的数据记作d₁，d₄为一个长度为20的向量，即为给的表格数据中的一行。将每行或每列闪烁时刻起N个采样数据作为有效数据，每个被试者均进行5轮实验，每轮实验行和列一共闪烁12次，所有每个被试者有60N有效数据。和问题一保持一致，仍然选取N=125。将被试者所有的有效数据组合成一个矩阵，并且用字母D表示， D=d1d2d60kT.D是维度为7500×20的矩阵，在压缩感知领域也称为“字典”。同样地，我们也将P300信号进行标签化处理，记作y，y是维度为7500的向量。如果相应的采样时刻出现P300信号，那么y对应的位置为1，否则为0。20个通道的权重向量为x。那么权系数求解问题就转换成一个线性回归模型：

y= Dx+ε                                             (7)

为了方便模型的推导，将字典D的维度记作P×K.ε为具有精度为σ的高斯零均值的随机变量。为了得到关于x的稀疏解, 我们引入“ spike and slab”稀疏促进函数[5]。

                           x|z-∏j=1K1-zjδxj+zjNxj|0τ0                   (8)

其中x₁是K维权重系数向量x的第i个分量，他们是K个满足均值为0，方差为τ₀且互相独立高斯变量。z称为支撑向量，是K个i. i. d伯努利分布的二维随机变量。 zᵢ=1意味着 xᵢ≠0,，进一步意味着通道i 是需要选择的，是对睡眠状态判定有着至关重要的通道； zᵢ=0意味着 xᵢ=0则意味这条通道是多余的。δ(·)代表冲激函数。支撑向量满足的分布为

1.  function Hd = bandpass

2.  %BANDPASS Returns a discrete- time filter object.

3.  % Butterworth Bandpass filter designed using FDESIGN. BANDPASS.

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  = ' stopband';                     % Band to match exactly

16.

17. % Construct an FDESIGN object and call its BUTTER method.

18. h  = fdesign. bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, …

19.                                     Astop2, Fs);

20. Hd = design(h, ' butter', 'MatchExactly', match);

21.

22.% [EOF]

2.随机森林预测模型 ( python 3.8) :

1.  from sklearn. ensemble import RandomForestClassifier

2.  from sklearn. model   selection import train   test   split

3.  from sklearn. metrics import accuracy   score

4.  from sklearn. preprocessing import StandardScaler

5.  import nump y as np

6.  import scipy. io as sio

7.

8. #一号被试预测模型

9. #数据导入

10. file = sio. loadmat(' train. mat')

11. s1   tem = file['S1  train']

12.

13. f ile1 = sio. loadmat(' trainevent. mat')

14. s1   eve = file1['S1  train   event']

15.

16. file2 = sio. loadmat('testEVENT. mat')

17. s1   test   tem = file2['S1  test']

61. scaler = StandardScaler() # 标准化转换

62. scaler. fit( predict   data)  # 训练标准化对象

63. predict   data= scaler. transform( predict   data)

64.

65. res1 = clf. predict( predict   data)

66. vote = np. zeros((10, 13))

67. index = 0

68. for i in range(0,20,2): #列

69.              for j in range(66): #行

70.                  h = s1   test   eve[j][i]

71.        ifh<100:

72.                         if res1[ index]:

73.                               vote[i//2,h] +=1

74.                         index+=1

3.基于EP的稀疏贝叶斯回归代码实现(python3.8)

1.  import nump y as np

2.

3.  import matplotlib. pyplot as plt

4.  import scipy. io as sio

5.  def reconstruct(A, y, P, K, sigma02, nvi2, tau0, alpha):

6.             ATA, ATy = np. dot(A. T, A), np. dot(A. T, y)

7.             J, H = ATA / sigma02, ATy / sigma02

8.             # Initialize Global Approximation

9.              mi, sigmai, nvi = np. zeros(A. shape[1]), np. ones(A. shape[1]), 0.5 *

10. np. ones(A. shape[1])

11.             mi   old = mi

12.             nvi   old = nvi

13.

14.            # Initialize Site Approximation

15.            mi2 sigmai2  inv, sigmai2  inv, nvi1 = np. zeros(A. shape[1]), 1e-6 *

16. np. ones(A. shape[1]), 0.5 * np. ones(K)

17.             for ittc s in range( max   itt):

18.                  # compute cavity

19.                   mi   sigmai   inv   bar, sigmai   inv   bar = mi / sigmai - mi2 sigmai2  inv, 1 /

20. sigmai - sigmai2  inv

21.                   nvi   bar = nvi2

22.

23.                  # positive cavity variance

24.                   index = sigmai   inv   bar < 0

25.                   update   index = sigmai   inv   bar > 0

26.                   sigmai   inv   bar[ index] = 1

27.

28.                  # Compute matching moments