CN104091172A - 一种运动想象脑电信号的特征提取方法 - Google Patents

Abstract

本发明涉及一种运动想象脑电信号的特征提取方法。所述方法首先将采集到的脑电信号进行预处理，然后对每导信号进行经验模态分解(EMD)，得到多阶的固有模态函数(IMF)信号，接着选取相同阶数的IMF信号作为新的信号，通过共同空间模式算法(CSP)获得空间滤波器，提取脑电信号的特征，并输入到分类器中进行分类，并根据分类准确率选取EMD和CSP中参数的最优值，最终获得最优参数下的脑电特征。本发明采用基于EMD与CSP对运动想象脑电信号进行特征提取，可以根据每个人脑电信号的特点，将信号自适应地分解成多个IMF信号，只需较少的电极就可以进行脑电信号的特征提取，并且从而较大程度地提高了脑电信号的分类准确率。

Description

一种运动想象脑电信号的特征提取方法

技术领域

本发明属于脑电信号处理技术领域，具体涉及脑-机接口(Brain-Computer Interface,BCI)系统中对运动想象脑电信号特征的提取方法，采用经验模态分解和共同空间模式相结合的方法对运动想象脑电信号进行特征提取。

背景技术

人口老龄化以及交通事故造成了大量脊髓病变或损伤的患者，其大脑的动作指令无法通过正常的体内通路传达到肌肉，从而丧失了肢体的运动能力。由于现代医学的进步，这些患者可以继续在轮椅和床上生存很长时间，但失去工作能力，生活难以自理，不仅患者十分痛苦，也给家人和社会带来负担。

BCI建立了人与计算机之间的一个交流和控制通道，是一种通过脑电信号实现人脑与计算机或其他电子设备间通讯和控制的系统，它不依赖外周神经和肌肉组织等这些常用的大脑输出通道。通过脑-机接口技术可以帮助脊髓受损的患者，实现对外部环境和设备的控制，提高生活自理能力。

BCI系统结构如图1所示，其基本过程为：脑电采集装置首先从大脑采集脑电信号并转换成数字信号，然后经过信号处理和模式识别模块对信号依次进行预处理、特征提取和模式分类，最后通过控制器输出控制信号，驱动外部装置做相应的动作。

临床医学研究表明，当人在做单侧肢体运动想象时，对侧大脑皮层相应区域的μ节律(8-12Hz)和β节律(18-23Hz)的振幅会较未做运动前有明显减小，这种现象称为事件相关去同步(event-related desynchronizations,ERD)；与此同时，同侧大脑皮层相应区域的μ节律和β节律的振幅会较未做运动前有明显增大，称为事件相关同步(event-related synchronizations,ERS)。脑电的这种节律性差异，可用于实现BCI技术。

共同空间模式(Common Spatial Pattern,CSP)方法被认为是提取运动想象脑电信号特征最为有效的方法之一。CSP算法利用代数上矩阵同时对角化的理论，寻找一组空间滤波器，使得在这组滤波器的作用下，一类信号的方差达到极大，另一类信号的方差达到极小，从而达到提取特征的目的。但是CSP算法存在两方面的不足：第一，CSP算法对大量电极的脑电信号效果显著，这限制了其在便携式BCI系统中的应用；第二，在执行运动想象任务时，每个受试者发生ERS/ERD的频段具有个体差异性，而CSP算法在提取脑电特征时，没有考虑受试者的个体差异性及脑电信号的频率特性，因此会造成分类准确率不高。

发明内容

针对现有CSP方法两方面的不足，本发明提出了一种基于经验模态分解(Empirical ModeDecomposition,EMD)与共同空间模式的运动想象脑电信号特征提取的方法，即EMD-CSP方法。此方法只需较少的电极就可以进行脑电信号的特征提取，并且可以根据每个人脑电信号的特点，将信号自适应地分解成多个固有模态函数(Intrinsic mode function,IMF)信号，从而较大程度地提高了脑电信号的分类准确率。

本发明采用的技术方案为：首先，将采集到的脑电信号进行预处理，然后对每导信号进行EMD，得到多阶的IMF信号，接着选取相同个数的IMF信号构成新的信号，通过CSP获得空间滤波器，提取出脑电信号的特征，将其输入到分类器中进行分类，并根据分类准确率选取EMD和CSP中参数的最优值，最终提取出最优参数下的脑电特征。

本发明方法的具体实现过程如下：

步骤1，信号采集及预处理。

首先通过脑电采集装置采集n导脑电信号，并将采集到的信号通过有限脉冲响应(FiniteImpulse Response,FIR)滤波器进行8-30Hz带通滤波，滤波后信号为x(t)＝[x₁(t),x₂(t),…,x_n(t)]∈R^N×n。其中N为样本点总数，n为脑电导联数目，x_i(t)为第i导滤波后的脑电信号，i＝1,2,…,n，t＝{1,2,…,N}。

步骤2，对步骤1得到的脑电信号x_i(t)(i＝1,2…,n)进行经验模态分解。

步骤2.1，确定x_i(t)的所有极值点，采用三次样条法对极大值点和极小值点进行拟合，得到上包络线e_imax(t)和下包络线e_imin(t)。

步骤2.2，计算上、下包络线的平均值，公式如下：

$m_i\left(t\right)=\frac{e_{i\max }\left(t\right)+e_{i\min }\left(t\right)}{2}---\left(1\right)$

步骤2.3，计算剩余信号r_i(t)。

令：

c_i(t)＝x_i(t)-m_i(t)        (2)

若c_i(t)不满足IMF整个信号段内极值点和过零点数量相等或至多相差1的筛选停止条件，则将c_i(t)取代x_i(t)重复步骤2.1～2.3；否则，取c_i(t)为分离出的一阶IMF分量，并按照下式计算剩余信号r_i(t)：

r_i(t)＝x_i(t)-c_i(t)        (3)

步骤2.4，令x_i(t)＝r_i(t)，重复步骤2.1～2.3的筛选过程，直到r_i(t)的极值点个数小于等于2时，终止经验模态分解过程。

在EMD过程结束时，假设x_i(t)被分解成k阶IMF分量和最终剩余信号r_ik(t)之和，即：

$x_i\left(t\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{c}_{\mathrm{ij}}\left(t\right)+r_{\mathrm{ik}}\left(t\right)---\left(4\right)$

其中，c_ij(t)表示第j次筛选得到的IMF分量。取前g阶IMF分量c_ij(t)(j＝1,2,…,g)构成其中g≤k，g的值可以根据分类准确率进行选取。从而，获得 $\hat{x}\left(t\right)={[{\hat{x}}_1\left(t\right),{\hat{x}}_2\left(t\right),...,{\hat{x}}_n\left(t\right)]}^T\∈R^{E\×N},$ 其中E＝g·n。

步骤3，采用共同空间模式方法提取脑电信号特征。

对两类运动想象A和B分别进行T_A、T_B次实验，T_A、T_B为正整数。CSP计算过程如下：

步骤3.1，计算混合空间的协方差。

首先，计算两类运动想象信号每次实验的协方差，公式如下：

$C=\frac{\hat{x}\left(t\right)\hat{x}{\left(t\right)}^T}{\mathrm{trace}\left\{\hat{x}\left(t\right)\hat{x}{\left(t\right)}^T\right\}}\∈R^{E\×E}---\left(5\right)$

其中，为矩阵的迹，即矩阵的对角线元素之和。

然后，分别计算两类运动想象的平均协方差：

$C_A=\frac{1}{T_A}\underset{i=1}{\overset{T_A}{\mathrm{\Σ}}}{C}_{A,i}---\left(6\right)$

$C_B=\frac{1}{T_B}\underset{i=1}{\overset{T_B}{\mathrm{\Σ}}}{C}_{B,i}---\left(7\right)$

其中，C_A,i、C_B,i分别表示运动想象A和B的第i次实验的协方差。

进而求得混合空间的协方差：

C_M＝C_A+C_B        (8)

步骤3.2，对混合空间协方差进行特征值分解，公式如下：

$C_M=U_M{\mathrm{\Λ}}_M{U}_M^T---\left(9\right)$

其中，U_M为特征向量矩阵，Λ_M为特征值对角矩阵。

步骤3.3，进行白化处理。

对Λ_M进行降序排序得到Λ_Md，并对U_M做同样的行列变换得到U_Md。令对C_A、C_B分别进行白化处理，如式如下：

S_A＝PC_AP^T        (10)

S_B＝PC_BP^T        (11)

利用S_A、S_B具有相同的特征向量的特点，经特征值分解后可得：

S_A＝BΛ_AB^T        (12)

S_B＝BΛ_BB^T        (13)

其中，B为S_A与S_B的共同特征向量，Λ_A、Λ_B分别为S_A和S_B的特征对角矩阵，且Λ_A+Λ_B＝I，I为单位矩阵。

因此求得空间滤波器矩阵为：

W＝B^TP        (14)

对进行W滤波得： $Z_0=W\·\hat{x}\left(t\right)\∈R^{E\×N}.$

步骤3.4，求特征向量f。

提取Z₀的前m行和后m行(m≤E/2)，构成Z＝[z₁,z₂,…,z_2m]^T∈R^2m×N，然后进行特征提取，计算公式如下：

$f_i=\frac{\mathrm{var}\left(z_i\right)}{\underset{j=1}{\overset{2m}{\mathrm{\Σ}}}\mathrm{var}\left(z_j\right)}---\left(15\right)$

其中，var(·)表示方差，i＝1,2,…,2m，则特征向量f＝[f₁,f₂,…,f_2m]^T。

步骤4，根据分类结果优化参数。

首先将脑电训练数据和测试数据按照步骤1、2、3提取特征向量，分别对分类器进行训练和测试，然后根据分类准确率选取步骤2.4中g和步骤3.4中m的最优值，即选取分类准确率最高时的g和m值。

与现有技术相比，本发明具有以下优点：

(1)本发明针对传统CSP方法提取运动想象脑电特征时需要大量导联的不足，先对脑电信号进行自适应经验模态分解，将一导信号变为多导信号，然后再输入到CSP中进行特征提取。实验表明，在只选用8导导联的情况下，本发明达到了95％的分类准确率，相比传统CSP方法相同条件下80％的分类准确率提高了很多，成功解决了传统CSP在导联较少情况下分类准确率低的问题，为便携式BCI系统的应用打下了坚实的基础。

(2)本发明针对CSP方法提取脑电特征时无法反映受试者的个体差异性和频率特性的不足，采用具有自适应特点的信号处理方法，即经验模态分解，将脑电信号从高频到低频自适应地分解为多阶IMF信号，并且根据分类准确率选取最优的g值(即前g阶IMF)和m值(即CSP方法中提取前m行和后m行特征)，从而获得与受试者运动想象最相关的IMF信号。实验表明，在采用相同的数据集情况下，比BCI Competition III第一名提高了4个百分点，充分说明本发明较大程度上提高了分类准确率。

附图说明

图1为BCI系统基本结构示意图；

图2为本发明所涉及方法的流程框图；

图3为8×8电极阵列植入位置示意图；

图4为实验时序示意图；

图5为任意一导脑电信号的经验模态分解结果图。

具体实施方式

本实施例是在Matlab2014a的仿真环境下进行的。

实施例采用BCI Competition III的Data Set I数据。如图3所示，BCI Competition III的Data Set I数据是由植入大脑右侧运动皮层的一个8×8(共64导)的电极阵列采集的，实验执行想象左手小指和想象舌头两种运动想象任务。信号的采样频率为1000Hz，经过0.016-300Hz的带通滤波。

实验时序图如图4所示，每次实验持续7s。0-1s时，屏幕出现一个十字光标；1-5s时，屏幕出现提示图片以提示测试者想象左手小指或舌头运动；5-7s屏幕处于白屏状态，受试者休息；随后进行下一次实验。Data Set I数据集为1.5-4.5s间共3000个采样点，包括278组训练数据(两种运动想象各139组)和100组测试数据(两种运动想象各50组)。

(1)数据预处理及导联选取

首先将Data Set I数据集的采样频率由1000Hz降为100Hz，将采样点数由3000降至300，然后对数据采用FIR带通滤波器，进行8-30Hz带通滤波。

为了应用于便携式BCI系统，需要减少导联的数目，并且由于ERS/ERD主要通过脑电信号的功率谱进行区分，故本例通过计算每导信号平均功率谱的Fisher Score来选择最优的8导脑电信号。Fisher Score的计算公式如下：

$F_i=\frac{{|P_{1,i}-P_{2,i}|}^2}{\mathrm{var}\left(P_{1,i}\right)+\mathrm{var}\left(P_{2,i}\right)}---\left(16\right)$

其中，P_k,i表示执行第k类运动想象时第i导脑电信号的平均功率谱，var(·)表示方差，F_i表示第i导脑电信号的Fisher Score。将64导脑电信号的Fisher Score按顺序从大到小排列，取前8导联，包括38、39、30、29、31、40、46和12。通过上述方法，得到降维后的训练数据集和测试数据集分别为X_Tr∈R^8×300×278和X_Te∈R^8×300×100。

(2)对优选导联的脑电信号进行经验模态分解

任取一导脑电信号x_i(t)∈R^1×300(i＝1,2,…,n)进行EMD，用三次样条法对极大值点和极小值点进行拟合，得到上包络线下包络线，进而求得第1阶IMF分量c_i1，然后判断剩余信号r_i1是否满足停止条件，若不满足，则重复上述步骤，最终得到c_i1、c_i2、c_i3、c_i4、c_i5、c_i6共6阶IMF分量和剩余信号r_i6，如图5所示。

取x_i(t)经过EMD后的第1阶到第g阶IMF，得同理，对其余脑电信号做上述处理，得到维数有效扩展后的训练数据集和测试数据集分别为X和X从而有利于CSP算法的运用。

(3)共同空间模式法提取脑电特征

令和分别为想象左手小指和舌头的脑电数据集。首先计算和的平均协方差，分别为C_f和C_t，进而求得混合空间的协方差矩阵C_M，对其进行特征值分解，求得特征向量矩阵U_M和特征值对角矩阵Λ_M；然后分别对C_f和C_t白化处理，得到相应的白化矩阵S_f和S_t，经过特征值分解后，得到共同的特征向量B，求得空间滤波器矩阵W；最后将通过W滤波得Z₀，取Z₀的第m行和最后m行构成Z，并进行特征提取。最终，得到训练特征集和测试特征集分别为F_Tr∈R^2m×278和F_Te∈R^2m×100。

(4)利用LibSVM进行分类并选取最优值

LibSVM是台湾大学林智仁(Lin Chih-Jen)副教授等开发设计的一个简单、易于使用和快速有效的支持向量机(Support Vector Machine,SVM)模式识别的软件包，该软件对SVM所涉及的参数调节相对较少。本实施例中选择径向基函数作为LibSVM分类器的核函数。

由于核参数γ、误差惩罚因子C以及EMD中g值和CSP特征提取中m值对LibSVM分类器的分类效果有很大的影响，故本实施例中采用遍历搜索法寻找最优的γ、C、g和m值。

第一步，γ和C以2^s变化，s的步长为1，s的取值范围为[-15,15]，即γ＝2^-15,2^-14,…,2¹⁵和C＝2^-15,2^-14,…,2¹⁵，g的取值范围为[1,5]，即g＝1,2,…,5，m的取值范围为[1,5g]，即m＝1,2,…,5g。取γ、C、g和m在取值范围内任意组合，以F_Tr为训练数据，F_Te为测试数据，选取分类准确率最高时γ、C、g和m的值。采用本发明EMD-CSP方法时，当γ＝2⁴、C＝2⁴、g＝2和m＝7时，分类准确率为94％。

第二步，为了更精确地选取γ和C的值，在第一步的基础上，s的步长改为0.1，即γ＝2^3.1,2^3.2,…,2^4.9和C＝2^3.1,2^3.2,…,2^4.9，g和m的取值范围不变，再次进行遍历搜索。当γ＝2^4.1、C＝2^4.3、g＝2和m＝7时，得到最高分类准确率为95％。

采用传统CSP方法按上述步骤进行特征提取，并用LibSVM进行分类，当γ＝2^2.4、C＝2^2.1和m＝2时，最高分类准确率为80％，可见本发明方法较CSP方法在分类准确率上有很大的提高。

表1是针对相同的Data Set I脑电数据集，采用本发明EMD-CSP方法及BCI CompetitionIII中获得前两名的识别方法、导联数目和相应的分类准确率。

表1本发明方法与脑机接口大赛前两名的分类准确率比较

特征提取方法	分类方法	导联数目	分类准确率
				CSSD+均值+FDA	SVM	64	91％
AR+功率谱+小波系数	正则化逻辑回归	64	87％
				EMD-CSP	SVM	8	95％

从表1中可知，本发明方法的分类准确率为95％，比BCI Competition III中第一名提高了4个百分点，并且本发明只选用了8导脑电信号进行分析，远小于竞赛第一名所选的64导脑电信号，充分说明本发明方法可以在选用较少导联情况下，取得较高的分类准确率。

经验模态分解可以根据信号特点进行自适应分解的性质以及共同空间模式算法具有排除相同任务成分提取不同任务成分的特点，使得EMD-CSP方法用于脑电特征提取时克服了共同空间模式方法需要较多导联和不能反映受试者个体差异性和频率特性的缺点，较大程度地提高了分类准确率，使得便携式BCI系统应用于实际中成为可能。

Claims (1)

1.一种运动想象脑电信号的特征提取方法，其特征在于包括以下步骤：

步骤1，信号采集及预处理；

首先通过脑电采集装置采集n导脑电信号，并将采集到的信号通过有限脉冲响应滤波器进行8～30Hz带通滤波，滤波后信号为x(t)＝[x₁(t),x₂(t),…,x_n(t)]∈R^N×n；其中N为样本点总数，n为脑电导联数目，x_i(t)为第i导滤波后的脑电信号，i＝1,2,…,n，t＝{1,2,…,N}；

步骤2，对步骤1得到的脑电信号x_i(t)(i＝1,2,…,n)进行经验模态分解；

步骤2.1，确定x_i(t)的所有极值点，采用三次样条法对极大值点和极小值点进行拟合，得到上包络线e_imax(t)和下包络线e_imin(t)；

步骤2.2，计算上、下包络线的平均值，公式如下：

$m_i\left(t\right)=\frac{e_{i\max }\left(t\right)+e_{i\min }\left(t\right)}{2}$

步骤2.3，计算剩余信号r_i(t)；

令：

c_i(t)＝x_i(t)-m_i(t)

若c_i(t)不满足固有模态函数IMF整个信号段内极值点和过零点数量相等或至多相差1的筛选停止条件，则将c_i(t)取代x_i(t)重复步骤2.1～2.3；否则，取c_i(t)为分离出的一阶IMF分量，并按照下式计算剩余信号r_i(t)：

r_i(t)＝x_i(t)-c_i(t)

步骤2.4，令x_i(t)＝r_i(t)，重复步骤2.1～2.3的筛选过程，直到r_i(t)的极值点个数小于等于2时，终止经验模态分解过程；

在经验模态分解过程结束时，假设x_i(t)被分解成k阶IMF分量和最终剩余信号r_ik(t)之和，即：

$x_i\left(t\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{c}_{\mathrm{ij}}\left(t\right)+r_{\mathrm{ik}}\left(t\right)$

其中，c_ij(t)表示第j次筛选得到的IMF分量；取前g阶IMF分量c_ij(t)(j＝1,2,…,g)构成其中g≤k，g的值可以根据分类准确率进行选取；从而，获得 $\hat{x}\left(t\right)={[{\hat{x}}_1\left(t\right),{\hat{x}}_2\left(t\right),...,{\hat{x}}_n\left(t\right)]}^T\∈R^{E\×N},$ 其中E＝g·n；

步骤3，采用共同空间模式方法提取脑电信号特征；

对两类运动想象A和B分别进行T_A、T_B次实验，T_A、T_B为正整数；

步骤3.1，计算混合空间的协方差；

首先，计算两类运动想象信号每次实验的协方差，公式如下：

$C=\frac{\hat{x}\left(t\right)\hat{x}{\left(t\right)}^T}{\mathrm{trace}\left\{\hat{x}\left(t\right)\hat{x}{\left(t\right)}^T\right\}}\∈R^{E\×E}$

其中，为矩阵的迹，即矩阵的对角线元素之和；

然后，分别计算两类运动想象的平均协方差：

$C_A=\frac{1}{T_A}\underset{i=1}{\overset{T_A}{\mathrm{\Σ}}}{C}_{A,i}$

$C_B=\frac{1}{T_B}\underset{i=1}{\overset{T_B}{\mathrm{\Σ}}}{C}_{B,i}$

其中，C_A,i、C_B,i分别表示运动想象A和B的第i次实验的协方差；

进而求得混合空间的协方差：

C_M＝C_A+C_B

步骤3.2，对混合空间协方差进行特征值分解，公式如下：

$C_M=U_M{\mathrm{\Λ}}_M{U}_M^T$

其中，U_M为特征向量矩阵，Λ_M为特征值对角矩阵；

步骤3.3，进行白化处理；

对Λ_M进行降序排序得到Λ_Md，并对U_M做同样的行列变换得到U_Md；令对C_A、C_B分别进行白化处理，如式如下：

S_A＝PC_AP^T

S_B＝PC_BP^T

利用S_A、S_B具有相同的特征向量的特点，经特征值分解后可得：

S_A＝BΛ_AB^T

S_B＝BΛ_BB^T

其中，B为S_A与S_B的共同特征向量，Λ_A、Λ_B分别为S_A和S_B的特征对角矩阵，且Λ_A+Λ_B＝I，I为单位矩阵；

因此求得空间滤波器矩阵为：

W＝B^TP

对进行W滤波得： $Z_0=W\·\hat{x}\left(t\right)\∈R^{E\×N};$

步骤3.4，求特征向量f；

提取Z₀的前m行和后m行(m≤E/2)，构成Z＝[z₁,z₂,…,z_2m]^T∈R^2m×N，然后进行特征提取，计算公式如下：

$f_i=\frac{\mathrm{var}\left(z_i\right)}{\underset{j=1}{\overset{2m}{\mathrm{\Σ}}}\mathrm{var}\left(z_j\right)}$

其中，var(·)表示方差，i＝1,2,…,2m，则特征向量f＝[f₁,f₂,…,f_2m]^T；

步骤4，根据分类结果优化参数；

首先将脑电训练数据和测试数据按照步骤1、2、3提取特征向量，并依据特征训练集和测试集分别对分类器进行训练和测试，然后根据分类准确率确定步骤2.4中g和步骤3.4中m的最优值，即选取分类准确率最高时的g和m值。