CN101690659B - 脑电波分析方法 - Google Patents

Abstract

一种脑电波分析方法，应用了经典的时频域分析和主成份分析方法解决了脑电信号特征提取的难题，成功地提取了与人体紧张，疲劳与放松息息相关的时频域参数，并且将其映射到主成分空间中去，还运用了支持向量机在主成份空间中高效地分析非线性关系，提高了判读的精确性和有效性。

Description

脑电波分析方法

技术领域

本发明涉及人体健康工程、神经生理学、生物医学数字信号处理技术、模式提取与模式分类、软件工程，尤其涉及一种脑电波分析方法。

背景技术

过快的都市节奏及激烈的市场竞争给我们带来的除了城市的繁荣外，更多的是白领们工作压力的日益递增，由此带来的结果就是越来越多的人群处于亚健康状态。亚健康的主要症状应包括疲劳、失眠、头痛、肌肉和关节疼痛、容易反复感冒、胸闷气短、心悸、食欲不振等躯体方面的表现；抑郁、焦虑、恐惧、冷漠、孤独、嫉妒、精神不振、情绪低落等心理方面的表现；以及不能很好地承担相应的社会义务及责任、人际关系紧张、家庭关系不和睦、工作学习困难等社会适应方面的表现。

放松训练指的是在一个安静的环境中，按一定的要求去完成特定的动作程序，通过反复练习，使人学会有意识地控制自体的心理生理活动，以降低机体的唤醒水平，调整因紧张性刺激而引起的机体心理--生理功能的紊乱。

随着亚健康人群越来越多，心理危机成为重要杀手，疲劳等级判定和放松治疗意义重大，临床上面，比如脑瘫，植物人的康复治疗也对脑电波的正确判读具有迫切需求，而且，临床上医生和心理医生的人工诊断具有主观性，医生需要长期培训，积累丰富的临床经验才能胜任，且他们的工作量有时候很巨大，比如人工判读长时间记录的脑电图，病理情况也是因人而异，但是治疗方法可能只是单一性的。

脑力疲劳的定量判定-脑波研究：

脑力疲劳是指工作或学习过程中，由于长时间从事脑力劳动而导致人们作业机能衰退、脑力能力下降的现象。其主要表现为工作或学习效率低下，反应迟缓，记忆力下降，注意力不集中以及协调性能变差。在交通驾驶、航空航天活动、人-机系统监控等工作中，操控人员瞬间的注意力分散、反应迟缓或协调性不够，都可能导致极为严重的事故。因此，分析、研究脑力疲劳及其评价标准，就显得尤为重要。

目前主要采用主观评价方法和行为学描述方法对脑力疲劳程度进行评估。主观评价方法根据作业者的身体感受症状来判断疲劳的程度，主要通过问卷调查的形式来进行，但其评分标准不易统一、且易受主观因素的影响，无法对疲劳时心理生理状态进行客观评定。而行为学描述方法主要通过脑力任务的反应时间和正确率对脑力疲劳程度进行评价。不但与被试者的主动配合程度密切相关，而且还有较强的“学习效应”，必须事先对被试者进行较长时问的训练，以排除“学习效应”影响，因而无法实时动态检测脑力疲劳。

近年来，脑电信号的检测与分析技术取得了很大的进步。在脑力疲劳研究中，脑电现已成为最广泛的评定中枢神经系统变化的指标之一，被誉为监测疲劳的“金标准”。Murata等利用事件相关电位P300的幅度和潜伏期大小对脑力疲劳进行分析，发现随着脑力疲劳程度的增加，潜伏期延长，幅度降低；Brookhuis等使用脑能量参数(theta+alpha)/beta对驾驶员进行研究，发现随着工作时间的增加、疲劳状态的加重，该参数呈上升趋势；Jung等用脑电功率谱对大脑警觉性水平进行分析，发现脑电功率谱可以反映大脑警觉状态的波动情况。虽然先前研究已取得了一些有意义的成果，但还不是很成熟。

大脑的皮质活动与人体的放松，疲劳和紧张状态息息相关。由于脑电信号可以分解为4个不同的节律，这4个不同的节律会随着疲劳状态的变化而变化。通过脑电采集仪对采集到的脑波进行分析评估后采用相应的方案进行放松训练从而有效的缓减紧张疲劳的心理生理状态。为脑力疲劳研究开辟一条新途径。

模式识别方法在脑电数字信号处理中的应用：

自1924年，德国神经精神病学家汉斯·贝格尔(Hans Berger)开始了对人类脑电活动的研究，并于1929年发表了名为“人类脑电图”的第一篇研究报告后，人们对脑信号(electroencephalogram，EEG)无论是从研究手段还是研究方法上都取得了显著的进展。在研究手段方面，随着电子和通讯技术的发展，数字脑电图仪和无线脑电图仪相继问世，为脑电图的记录和分析提供了更多的便捷。

近年来，数字信号处理技术与模式识别方法蓬勃发展，在生物医学信号处理中的应用也越来越为人们所关注。从经典的时频域信号处理方法到现代的时频域、高阶累积量信号处理方法，从有维数灾难和易陷入局部极小的神经网络到适合高阶数据分析且全局最优的支持向量机，模式识别方法将继续向前发展，不断推进人类去揭开隐藏在大脑深处的奥秘。

在对脑电波的分析方面，存在以下难点：

难点一：脑电信号与人脑状态的特征提取。

如果说信号是携带有客观物体特征或状态的信息载体，则脑电信号就是携带有人脑特征或状态的信息载体。人脑是一个开放的、时变的和非线性的，主要是确定性特征的系统，因此，这样的系统产生的信号，也应该是时变的、非线性的，主要是确定性的。由于脑电信号在测量后会引入随机误差，并且信号本身存在随机干扰，这些都使信号随机化。随机信号处理技术的应用有助于提取包含在随机化了的信号中的确定成分(真实脑电信号)。

经典的信号处理技术针对的是线性的、非时变问题，而对于时变的、非线性问题，往往先提出一些近似合理的假设，即在相当短的时间内，认为脑电信号是线性的、非时变的。而针对非线性脑电信号处理技术也不断丰富起来，如时频分析技术、高阶谱分析技术、复杂度分析技术等。

难点二：对人脑的状态进行有效地等级判定。

由于受到脑电信号本身的时变性、非线性和随机性的影响，此外还受到个体与个体之间的差异性影响，对人脑状态的正确评定和划分成为了难点问题。

发明内容

本发明提供了一种脑电波分析方法，解决了背景技术中不能进行量化分析，实时性差，个体针对性不强的问题，应用了经典的时频域分析和主成份分析方法解决了脑电信号特征提取的难题，成功地提取了与人体紧张，疲劳与放松息息相关的时频域参数，并且将其映射到主成分空间中去，还运用了支持向量机在主成份空间中高效地分析非线性关系，提高了判读的精确性和有效性。

为了达到上述目的，本发明提供了一种脑电波分析方法，包含以下步骤：

步骤1、采集脑电数据；

8导脑电：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

1导皮肤阻抗SCR；

1导心电ECG；

1导水平眼电HEOG；

1导垂直眼电VEOG；

1导肌电EMG；

步骤2、对采集到的脑电数据进行时域分析；

步骤2.1、脑电回放(EEG review)，了解本次采集数据有无明显的伪差；引起产生伪差的原因主要有三个方面：

1、自身的生理现象和精神因素引起的伪差；

2、脑电记录仪本身带来的伪差；

3、空间电磁干扰产生的伪差；

步骤2.2、剔除明显漂移的脑电数据；

步骤2.3、观察眼电的方向与眼电伪迹的方向是否一致，以及眼电的幅值大小，以指导后续的去除眼电的步骤；

步骤2.4、观察心电导联和肌电导联的基本特征和幅值大小，目测对脑电影响的位置和程度，以指导后续的伪迹剔除的步骤；

步骤3、对脑电(EEG)进行预处理；

步骤3.1、去除生理伪迹；

去除眨眼伪迹(Blink)，心电伪迹(EKG)，肌电伪迹(EMG)；

采用相关法去除眨眼伪迹：

矫正的EEG＝原始EEG-b×EEG                        (1)

其中，b为传递系数，通过分别估算眼电导联和每导EEG的协方差来计算：

b＝cov(EOG，EEG)/var(EOG)

其中，cov为协方差“covariance”，var代表方差“variance”；

步骤3.2、去除技术伪迹；

去除由50周干扰(市用50Hz交流电混入脑电引起的在整个导联上都会出现50Hz的工频干扰)和阻抗过高引起的伪迹等；

选用数字滤波来去除；

步骤4、对脑波进行频域分析；

步骤4.1、先将采集得到的脑电时序列数据分段保存，经快速傅立叶变换(FFT)得到傅立叶成分；

N点有限长序列x(n)的离散傅立叶变换(DFT)为：

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{\mathrm{nk}},0\≤k\≤N-1.---\left(2\right)$

X(k)为变换得到的傅立叶成分，W_N来表示正交序列集中的基

即 $W_N^{\mathrm{nk}}=e^{-j\frac{2\mathrm{\π}}{N}\mathrm{nk}};$

步骤4.2、由傅立叶成分计算功率谱：

P(k)＝|X(k)|²                                    (3)

步骤4.3、由功率谱成份计算得出各个导联的脑电参数：

1、幅值[μv](符号为A_z(c))： $A_z\left(c\right)=4\sqrt{S_z\left(c\right)};$

2、百分比[％](符号为P_z(c))：P_z(c)＝S_z(c)/S_T(c)；

3、重心频率[Hz](符号为F_z(c))： $F_z\left(c\right)=(\underset{z}{\mathrm{\Σ}}p\left(f\left(c\right)\right)\×f)/\left(\underset{z}{\mathrm{\Σ}}p\right(f\left(c\right));$

其中，c代表各个脑电导联：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

z代表脑电的各个频率段：δ(0.5-4Hz)，θ(4-8Hz)，α(8-13Hz)，β(13-25Hz)和T(0.5-25Hz)；

p(f(c))是频率f(f取值0.5-25hz)处的脑电功率谱成份；

S_z(c)是在z频段范围内c导联处的脑电功率谱成份之和；

S_T(c)是T频段范围内(0.5～25Hz)c导联处的脑电功率谱成份之和；

步骤5、提取脑波特征，进行主成分分析；

主成分分析(principle component analysis，PCA)技术是依据原始数据的线性变换，由原始数据互相关矩阵的主要特征的大小来确定坐标变换和变量压缩，目的是在数据空间中找到一组m维的正交基，这组正交基最大可能地表示数据的方差和协方差，以便将数据从原始的l维空间映射到由这组正交基所构成的m维子空间上(m＜l)，从而达到降维和特征提取的目的。正是由于PCA技术原理清楚、计算简单、实现方便，因此被广泛的应用到信号分析的领域。

步骤5.1、选取脑电频域参数来构造原始数据S；

原始数据S的构成为：

S＝[s₁，s₂……s_p]^T                                        (4)

其中，x₁，s₂……s_p为p个脑波频率参数；

步骤5.2、进行KMO(Kaiser-Meyer-Olkin检验方法)检验，判断原始数据S是否适合做主成分分析；

KMO(Kaiser-Meyer-Olkin measure)，用于检验变量间的偏相关性是否足够小；

KMO计算公式如下：

$\mathrm{KMO}=\frac{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2}{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2+\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{a}_{\mathrm{ij}}^2}---\left(5\right)$

其中，r_ij ²为两变量间的简单相关系数，a_ij ²为两变量间的偏相关系数，KMO统计量取值在0至1之间，其值越大，做主成分分析的效果越好，当KMO＜0.5时，不适合做主成分分析；

步骤5.3、计算特征值和主成分；

设原始数据S＝[s₁，s₂……s_p]^T是p维随机向量，它的协方差阵的p个特征值为λ₁≥λ₂≥……≥λ_p≥0，相应的标准正交化的特征向量为α₁，α₂……α_p；

原始数据S的第i个主成分为：

$f_i={\mathrm{\α}}_i^{T}S---\left(6\right)$

特征值为：

且方差λ_i＝Var(f_i)，i＝1，2……p；                (7)

其中，Var( )为方差函数；

步骤5.4、根据特征值λ计算累计方差贡献率C，确定主成分个数q；

$C=\frac{{\mathrm{\Σ}}_{i=1}^q{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{i=1}^p{\mathrm{\λ}}_i}---\left(8\right)$

其中，p表示原始空间中特征的个数，q表示所取的主成分个数；通常，前q个主成分的累计贡献率C要求达到70％-90％；

步骤5.5、计算特征主成分的负荷矩阵：

主成分的负荷矩阵相应由步骤5.3中的特征向量α₁，α₂……α_p生成，前q个主成分的负荷矩阵表示为[α₁，α₂……α_q]；

步骤6、采用支持向量机进行脑波模式分类；

支持向量机是Vapnik在20世纪90年代中期提出的，它是一种建立在小样本理论基础上的通用学习机器。它不仅能够避免传统方法中出现的局部极小和过度拟合问题，而且通过引入核函数(kernel function)，使人们能够在高维空间中高效地分析非线性关系。支持向量机在分类中的优越性，不仅体现在分类正确率的提高，而且体现在训练速度的提升；

步骤6.1、构造q维的输入向量x＝(x₁，x₂……x_i……x_q)；

其中，q为步骤5.4中确定的主成分个数，x_i＝f_i是第i个主成分。

将步骤5中提取到的特征主成分代入上式；

步骤6.2、计算支持向量机的决策函数；

支持向量机的决策函数为：

$y=f\left(x\right)=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{a}_i^{*}{y}_{i}K(x_i,x)+b^{*}\},---\left(9\right)$

其中，sgn( )为符号函数，a_i ^*为最优解，K(x_i，x)表示内积计算，相当于把原来的特征空间变换到某一新的特征空间；设共有n个训练样本(xi，yi)，yi是训练集中样本的类别标号，y_i∈{+1，-1}，b°是新的特征空间中原点到分类面的距离，y＝f(x)是测试集中样本的类别标号，也就是最终判别出的状态；

步骤6.3、输出分类结果：y∈{+1，-1}，+1是疲劳状态，-1是放松状态。

本发明具有如下优点：

1、定量判定：利用脑电信号可以定量地测量和判读人的生理，心理和情绪状况；利用数字信号处理技术和特征提取算法能够提取与人的紧张，焦虑，抑郁情绪相对应的电生理参数，利用模式识别技术可以正确判定疲劳等级，区分疲劳和放松等不同的精神状态。

2、实时评价：分析脑电与人体精神状态的关系，实时地给出疲劳判定和放松效果的评价。

3、个体针对性：基于自身脑电的分析更加具有康复治疗的针对性，治疗方案可以因人而异；

4、广泛的适用人群：适用既面向正常人群，又面向临床上的病例，具有广泛的应用价值；

5、精确的判读结果：从医生定性的诊断治疗到脑电信号定量的测评，诊断和治疗，大大减少人工分析的主观性，提高工作效率和分析结果的一致性。

6、良好的通用性：适用于各种脑电采集仪采集的数据。

本发明应用了经典的时频域分析和主成份分析方法解决了脑电信号特征提取的难题。又应用了支持向量机的方法对脑力疲劳和放松效果给出了正确性很高的判读结果。

附图说明

图1是实施例中主成分分析中输入变量数目与特征值的关系图；

图2是实施例中支持向量机的示意图；

图3是实施例中训练集和测试集在主成分空间中的分布示意图。

具体实施方式

以下根据图1～图3，具体说明本发明的较佳实施例：

一种脑电波分析方法，包含以下步骤：

步骤1、采集脑电数据；

在25℃室温下安静独立房间中进行脑电数据采集；

参照国际10-20系统安放电极，双极导联，使用日本光电(Nihon-KodenEEG1100)进行采集，采样频率为200Hz，采集30Hz以下的脑电成分，共采集13导数据，包括：

8导脑电：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

1导皮肤阻抗SCR；

1导心电ECG；

1导水平眼电HEOG；

1导垂直眼电VEOG；

1导肌电EMG；

步骤2、对采集到的脑电数据进行时域分析；

步骤2.1、脑电回放，了解本次采集数据有无明显的伪差；

步骤2.2、剔除明显漂移的脑电数据；

步骤2.3、观察眼电的方向与眼电伪迹的方向是否一致，以及眼电的幅值大小，以指导后续的去除眼电的步骤；

步骤2.4、观察心电导联和肌电导联的基本特征和幅值大小，目测对脑电影响的位置和程度，以指导后续的伪迹剔除的步骤；

步骤3、对脑电(EEG)进行预处理；

步骤3.1、去除生理伪迹；

去除眨眼伪迹(Blink)，心电伪迹(EKG)，肌电伪迹(EMG)；

采用相关法去除眨眼伪迹：

矫正的EEG＝原始EEG-b×EEG                        (1)

其中，b为传递系数，可以通过分别估算眼电导联和每导EEG的协方差来计算：

b＝cov(EOG，EEG)/var(EOG)；

其中，cov为协方差“covariance”，var代表方差“variance”；

步骤3.2、去除技术伪迹；

去除由50周干扰和阻抗过高引起的伪迹等；

选用数字滤波来去除；

步骤4、对脑波进行频域分析；

步骤4.1、先将采集得到的脑电时序列数据分段保存(每段数据5.12秒，含1024个采样点)，经快速傅立叶变换(FFT)得到傅立叶成分，采用基于库利-图基型(频率抽选的以2为基底)的FFT算法；

N点有限长序列x(n)的离散傅立叶变换(DFT)为：

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{\mathrm{nk}},0\≤k\≤N-1.---\left(2\right)$

X(k)为变换得到的傅立叶成分，W_N来表示正交序列集中的基

即 $W_N^{\mathrm{nk}}=e^{-j\frac{2\mathrm{\π}}{N}\mathrm{nk}};$

步骤4.2、由傅立叶成分计算功率谱(周期图periodogram)：

P(k)＝|X(k)|²                            (3)

步骤4.3、由功率谱成份计算得出各个导联的脑电参数：

1、幅值[μv](符号为A_z(c))： $A_z\left(c\right)=4\sqrt{S_z\left(c\right)};$

2、百分比[％](符号为P_z(c))：P_z(c)＝S_z(c)/S_T(c)；

3、重心频率[Hz](符号为F_z(c))： $F_z\left(c\right)=(\underset{z}{\mathrm{\Σ}}p\left(f\left(c\right)\right)\×f)/\left(\underset{z}{\mathrm{\Σ}}p\right(f\left(c\right));$

其中，c代表各个脑电导联：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

z代表脑电的各个频率段：δ(0.5-4Hz)，θ(4-8Hz)，α(8-13Hz)，β(13-25Hz)和T(0.5-25Hz)；

p(f(c))是频率f(f取值0.5-25hz)处的脑电功率谱成份；

S_z(c)是在z频段范围内c导联处的脑电功率谱成份之和；

S_T(c)是T频段范围内(0.5～25Hz)c导联处的脑电功率谱成份之和；

步骤5、提取脑波特征，进行主成分分析；

步骤5.1、选取脑电频域参数来构造原始数据S；

选取脑电9个频域参数来构造S，它们包括前额中线的theta波，中央区的beta波，枕区的alpha波和beta波；

原始数据S＝[s₁，s₂……s₉]，其中s₁，s₂……s₉为9个脑波频率参数，它们分别是A_δ(Fz)，P_δ(Fz)，A_θ(Fz)，P_θ(Fz)，A_α(O1)，P_α(O1)，A_β(O1)，P_β(O1)，A_T(O1)；

测试集原始数据S如下表所示：

步骤5.2、进行KMO检验，判断原始数据S是否适合做主成分分析；

$\mathrm{KMO}=\frac{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2}{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2+\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{a}_{\mathrm{ij}}^2}---\left(5\right)$

其中，r_ij ²为两变量间的简单相关系数，a_ij ²为两变量间的偏相关系数；

这里KMO计算了9个脑波频率参数之间的偏相关性大小，KMO＝0.569＞0.5，适合做主成分分析；

步骤5.3、计算特征值和主成分；

设原始数据S＝[s₁，s₂……s₉]^T是9维随机向量，它的协方差阵的9个特征值为λ₁≥λ₂≥……≥λ₉≥0，相应的标准正交化的特征向量为α₁，α₂……α₉；

原始数据S的第i个主成分为：

$f_i={\mathrm{\α}}_i^{T}S---\left(6\right)$

特征值为：

且方差λ_i＝Var(f_i)，i＝1，2……9；                (7)

其中，Var( )为方差函数；

计算结果如图1所示，横轴为输入变量数目，纵轴为对应的特征值λ₁，λ₂……λ₉；

步骤5.4、根据特征向量值λ计算累计方差贡献率C，确定主成分个数q；

$C=\frac{{\mathrm{\Σ}}_{i=1}^q{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{i=1}^p{\mathrm{\λ}}_i}---\left(8\right)$

其中，p表示原始空间中特征的个数，q表示所取的主成分个数；通常，前q个主成分的累计贡献率C要求达到70％-90％；

计算结果如下表所示：

由上表中看出，前3个主成分的累计贡献已经达到了85.29％，所以在本实施例中，确定主成分个数q＝3，选取前3个特征主成分f1，f2，f3；

步骤5.5、计算特征主成分的负荷矩阵：

本实施例中，提取了前3个特征主成分f1，f2，f3，每一个主成分实际都是原始输入参数的一个线性组合，主成分的负荷矩阵见下表：

由上表的负荷矩阵，前3个主成分f1，f2，f3可以由输入变量s1-s9表示如下：

f1＝-0.634s1-0.867s2+0.385s3+0.402s4+0.319s5-0.958s6-0.592s7+0.964s8+0.783s9；

f2＝0.689s1+0.350s2+0.685s3+0.868s4+0.354s5-0.040s6-0.268s7+0.055s8-0.204s9；

f3＝-0.227s1+0.102s2+0.594s3+0.149s4-0.805s5-0.023s6+0.329s7-0.033s8+0.207s9；

前3个主成分的负荷矩阵表示为[α₁，α₂，α₃]；

其中：

α₁ ^T＝[-0.634，-0.867，0.385，0.402，0.319，0.958，-0.592，0.964，0.783]；

α₂ ^T＝0.689，0.350，0.685，0.868，0.354，-0.040，-0.268，0.055，-0.204]；

α₃ ^T＝[-0.227，0.102，0.594，0.149，-0.805，-0.023，0.329，-0.033，-0.207]；

步骤6、采用支持向量机进行脑波模式分类；

步骤6.1、构造输入向量x＝(x₁，x₂……x_i……x_q)，其中q为步骤5.4中确定的主成分个数，x_i＝f_i是第i个主成分。将步骤5中提取到的特征主成分代入上式；

步骤6.2：计算支持向量机的决策函数，如图2所示。

支持向量机的决策函数为：

$y=f\left(x\right)=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{a}_i^{*}{y}_{i}K(x_i,x)+b^{*}\},---\left(9\right)$

其中，sgn( )为符号函数，a_i ^*为最优解，K(x_i，x)表示内积计算，相当于把原来的特征空间变换到某一新的特征空间；设共有n个训练样本(xi，yi)，yi是训练集中样本的类别标号，y_i∈{+1，-1}，b°是新的特征空间中原点到分类面的距离，y＝f(x)是测试集中样本的类别标号，也就是最终判别出的状态；

选取分析数据样本，分为训练集和测试集；

训练集包括从7名被测试者身上记录所得的39段脑电时序列(每段5.12s，1024采样点，200Hz采样率)，其中21段记录状态为心算疲劳状态，18段记录状态为听音乐放松状态；

测试集包括从这7名相同的被测试者身上记录得到的26段脑电时序列(每段5.12s，1024采样点，200Hz采样率)，其中17段记录状态为心算疲劳状态，9段记录状态为听音乐放松状态；

训练集和测试集数据在三维主成分f1，f2，f3空间中的分布散点图如图3所示，可以看到在两种不同大脑状态下的训练集和测试集数据都是属于非线性可分的，而且测试集数据表现出与训练数据很相似的分布；

步骤6.3、输出分类结果：y∈{+1，-1}，+1是疲劳状态，-1是放松状态；下表为测试集数据的实际状态和支持向量机的判别分类结果：

上表的实验结果表明支持向量机完全正确地区分了测试集中的26段脑电时序列的两种脑波状态。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (3)

1.一种脑电波分析方法，其特征在于，包含以下步骤：

步骤1、采集脑电数据；

所述的步骤1中，采集如下数据：

8导脑电：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

1导皮肤阻抗SCR；

1导心电ECG；

1导水平眼电HEOG；

1导垂直眼电VEOG；

1导肌电EMG；

步骤2、对采集到的脑电数据进行时域分析；

步骤2.1、脑电回放，了解本次采集数据有无明显的伪差；

步骤2.2、剔除明显漂移的脑电数据；

步骤2.3、观察眼电的方向与眼电伪迹的方向是否一致，以及眼电的幅值大小，以指导后续的去除眼电的步骤；

步骤2.4、观察心电导联和肌电导联的基本特征和幅值大小，目测对脑电影响的位置和程度，以指导后续的伪迹剔除的步骤；

步骤3、对脑电EEG进行预处理；

步骤3.1、去除生理伪迹；

去除眨眼伪迹Blink，心电伪迹EKG，肌电伪迹EMG；

步骤3.2、去除技术伪迹；

去除由50周干扰和阻抗过高引起的伪迹；

步骤4、对脑波进行频域分析；

步骤4.1、先将采集得到的脑电时序列数据分段保存，经快速傅立叶变换得到傅立叶成分；

N点有限长序列x(n)的离散傅立叶变换(DFT)为：

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{\mathrm{nk}},$ 0≤k≤N-1.                                        (2)

X(k)为变换得到的傅立叶成分，W_N来表示正交序列集中的基

即 $W_N^{\mathrm{nk}}=e^{-j\frac{2\mathrm{\π}}{N}\mathrm{nk}};$

步骤4.2、由傅立叶成分计算功率谱：

P(k)＝|X(k)|²                                (3)

步骤4.3、由功率谱成份计算得出各个导联的脑电参数：

1、幅值： $A_z\left(c\right)=4\sqrt{S_z\left(c\right)};$

2、百分比：P_z(c)＝S_z(c)/S_T(c)；

3、重心频率： $F_z\left(c\right)=(\underset{z}{\mathrm{\Σ}}p\left(f\left(c\right)\right)\×f)/\left(\underset{z}{\mathrm{\Σ}}p\left(f(c\right)\right);$

其中，c代表各个脑电导联：Fp1-A1，F3-A1，O1-A1，Fp2-A2，F4-A2，O2-A2，Fz-Aav和Cz-Aav；

z代表脑电的各个频率段：δ频段范围是0.5-4Hz，θ频段范围是4-8Hz，α频段范围是8-13Hz，β频段范围是13-25Hz和T频段范围是0.5-25Hz；

p(f(c))是频率f处的脑电功率谱成份，f取值0.5-25Hz；

S_z(c)是在z频段范围内c导联处的脑电功率谱成份之和；

S_T(c)是T频段范围内c导联处的脑电功率谱成份之和；

步骤5、提取脑波特征，进行主成分分析；

步骤5.1、选取脑电频域参数来构造原始数据S；

原始数据S的构成为：

S＝[s₁，s₂……s_p]^T                                (4)

其中，s₁，s₂……s_p为p个脑波频率参数；

步骤5.2、进行KMO检验，判断原始数据S是否适合做主成分分析；

KMO计算公式如下：

$\mathrm{KMO}=\frac{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2}{\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{r}_{\mathrm{ij}}^2+\mathrm{\Σ}\underset{i\≠j}{\mathrm{\Σ}}{a}_{\mathrm{ij}}^2}---\left(5\right)$

其中，

为两变量间的简单相关系数，

为两变量间的偏相关系数，KMO统计量取值在0至1之间，其值越大，做主成分分析的效果越好，当KMO＜0.5时，不适合做主成分分析；

步骤5.3、计算特征值和主成分；

设原始数据S＝[s₁，s₂……s_p]^T是p维随机向量，它的协方差阵的p个特征值为λ₁≥λ₂≥……≥λ_p≥0，相应的标准正交化的特征向量为α₁，α₂……α_p；

原始数据S的第i个主成分为：

$f_i={\mathrm{\α}}_i^{T}S---\left(6\right)$

特征值为：

λ_i＝Var(f_i)，i＝1，2……p；                    (7)

其中，Var()为方差函数；

步骤5.4、根据特征值λ计算累计方差贡献率C，确定主成分个数q；

$C=\frac{{\mathrm{\Σ}}_{i=1}^q{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{i=1}^p{\mathrm{\λ}}_i}---\left(8\right)$

其中，p表示原始空间中特征的个数，q表示所取的主成分个数；前q个主成分的累计贡献率C要求达到70％-90％；

步骤5.5、计算特征主成分的负荷矩阵：

主成分的负荷矩阵相应由步骤5.3中的特征向量α₁，α₂……α_p生成，前q个主成分的负荷矩阵表示为[α₁，α₂……α_q]；

步骤6、采用支持向量机进行脑波模式分类；

步骤6.1、构造q维的输入向量x＝(x₁，x₂……x_i……x_q)；

其中，q为步骤5.4中确定的主成分个数，x_i＝f_i是第i个主成分；

将步骤5中提取到的特征主成分代入步骤6.1中的公式；

步骤6.2、计算支持向量机的决策函数；

支持向量机的决策函数为：

$y=f\left(x\right)=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{q}_i^{*}{y}_{i}K(x_i,x)+b^{*}\},---\left(9\right)$

其中，sgn()为符号函数，

为最优解，K(x_i，x)表示内积计算，相当于把原来的特征空间变换到某一新的特征空间；设共有n个训练样本(x_i，y_i)，y_i是训练集中样本的类别标号，y_i∈{+1，-1}，b^*是新的特征空间中原点到分类面的距离，y＝(x)是测试集中样本的类别标号，也就是最终判别出的状态；

步骤6.3、输出分类结果：y∈{+1，-1}，+1是疲劳状态，-1是放松状态。

2.如权利要求1所述的脑电波分析方法，其特征在于，所述的步骤3.1中，采用相关法去除眨眼伪迹：

矫正的EEG＝原始EEG-b×EEG            (1)

其中，b为传递系数，通过分别估算眼电导联和每导EEG的协方差来计算：

b＝cov(EOG，EEG)/var(EOG)

其中，cov为协方差“covariance”，var代表方差“variance”。

3.如权利要求1所述的脑电波分析方法，其特征在于，所述的步骤3.2中，选用数字滤波来去除技术伪迹。

Images