Computer Engineering and Applications 计算机工程与应用
2017,53(13)
9
公共空间模式算法结合经验模式分解的 EEG 特征提取
张学军 1,2,黄婉露 1,黄丽亚 1,2,成谢锋 1,2
ZHANG Xuejun1,2, HUANG Wanlu1, HUANG Liya1,2, CHENG Xiefeng1,2
1. 南京邮电大学 电子科学与工程学院,南京 210023
2. 江苏省射频集成与微组装工程实验室,南京 210023
1.School of Electronic Science and Engineering, Nanjing University of Posts and Telecommunications, Nanjing 210023, China
2.Jiangsu Province Engineering Lab of RF Integration & Micropackage, Nanjing 210023, China
ZHANG Xuejun, HUANG Wanlu, HUANG Liya, et al. EEG signals feature extraction combined with empirical
mode decomposition and common spatial pattern. Computer Engineering and Applications, 2017, 53(13):9-15.
Abstract: Normal Common Spatial Pattern(CSP)method is restricted to the abundant input channels and lacking fre-
quency information. This paper puts forward an improved CSP method combined with the Empirical Mode Decomposition
(EMD-CSP)to achieving feature vector by a different component choice. Firstly, the EMD method is proposed to decom-
pose the EEG signal into a set of stationary time series called Intrinsic Mode Functions(IMF). Secondly, these IMFs are
analyzed with the band-power to detect the valuable IMFs with characteristics of sensorimotor rhythms(5~28 Hz), and
then the improved CSP filter is attached to the feature extraction of screening IMFs. Finally, once the feature vector is
built, the classification of MI is performed using Support Vector Machine(SVM). The results obtained show that the
EMD-CSP allow the most reliable features and that the accurate classification rate obtained is 92% which confirms the
feasibility and availability of this method.
Key words: Electroencephalogram(EEG); Empirical Mode Decomposition(EMD); Common Spatial Pattern(CSP)
摘 要:常规的公共空间模式分解方法需要大量的输入通道、缺乏频域信息,发展受到限制。为了克服以上缺点,将
经验模式分解(Empirical Mode Decomposition,EMD)和公共空间模式算法结合,改变 CSP 滤波器成分选择方式,提
出 EMD-CSP 算法来获取特征向量。该算法对预处理后的信号进行经验模式(EMD)分解 ,得到固有模态函数
(Intrinsic Mode Functions,IMFs),观察并计算每个 IMF 分量的能量谱,筛选有效的 IMF 频段(5~28 Hz),使用改进
的 CSP 滤波器进行滤波获取特征,最后使用支持向量机(Support Vector Machine,SVM)进行分类。分类结果得到 9
位受试的想象运动平均分类正确率为 92%,证实了该算法的可行性与有效性。
关键字:脑电信号 ;经验模式分解 ;公共空间模式分解
文献标志码:A 中图分类号:R318;TP39
doi:10.3778/j.issn.1002-8331.1703-0199
1 引言
脑机接口是一种基于脑电信号实现人脑与计算机
或其他设备之间通信和控制的接口[1],通过大脑活动产
生的脑电信号(Electroencephalogram,EEG)为人类和环
境之间提供了新的通信和控制渠道,拓宽了人类的控制
能力,具有广泛的应用前景。该领域的发展不仅帮助瘫
痪病人使用计算机、神经假体、机器臂等电子设备,也
实现了包括:运动恢复、通信、环境控制甚至娱乐等其他
功能[2]。
目前基于 EEG的脑机接口主要包含 5个主要部分:获
取信号、预处理、特征提取、特征分类和接口设备控制。
在以上步骤中:因为有效地在降维特征空间展现了输入
信号特征并适用于进一步标识不同想象运动脑电信号
的判别信息,特征提取在 BCI 研究界受到广泛的关注[2]。
基金项目:国家自然科学基金(No.61271334)。
作者简介:张学军(1969—),男,博士,教授,主要研究方向为智能信息处理、认知网络频谱感知、无线射频识别技术等,E-mail:
xjzhang@njupt.edu.cn;黄婉露(1993—),女,硕士研究生,主要研究领域为智能信息处理;黄丽亚(1972—),女,博士,教
授;成谢锋(1956—),男,教授,博士生导师。
收稿日期:2017-03-13 修回日期:2017-05-03 文章编号:1002-8331(2017)13-0009-07
计算机工程与应用www.ceaj.org
10
2017,53(13)
Computer Engineering and Applications 计算机工程与应用
作 为 基 于 运 动 想 象 的 脑 机 接 口(Motor Imagery-
Brain Computer Interface,MI-BCI)系统的信号源,脑电
图(EEG)和脑皮层电图(Electrocorticography,ECoG),
EEG 信号是一种弱的、非线性、非平稳并且随时间变化
的信号。因此提出一种有效的特征提取方法是改善识
别精度的关键。目前,时间-频率分析作为处理智能信
号一种潜在的强大方法被广泛应用于脑信号的研究。
传统的时间 -频率方法包括:短时傅里叶变换(Short-
Time Fourier Transform,STFT)、小 波 变 换(Wavelet
Transform,WT)、Winger-Ville 分布,但是这些方法的本
质都是基于傅里叶变换,根据海森堡不确定性原理,该
方法不可能同时得到时间-频率的良好分辨率。近年,
希尔伯特黄变换( Hilbert-Huang Transform,HHT)作为
另一种时间-频率分析法已经变得越来越流行,其可以
有效地分析非线性和非平稳信号。原信号经过经验模
式分解(EMD)被分解为一系列固有模态函数(IMFs),
随后对每个固有模态函数进行希尔伯特黄变换。HHT
不涉及海森堡不确定性原理可以获得时域和频域的高
分辨率。目前它被广泛应用于许多信号处理领域,如雷
达探测、地震信号和生物医学信号等[3]。
再者,由于 EEG 信号低空间分辨率,EEG 信号构成
的 BCI 系统需要进行有效的空间滤波,从而确保从受试
的相关脑域中提取特征信息。在这一方面,常用的算法
有:共域空间分解(CSP)、独立主成分分析(Independent
Component Analysis,ICA)和共域空间谱模式(CSSP)、
滤 波 器 CSP(Filter Bank Common Spatial Pattern,
FBCSP)、判别滤波 CSP(Discriminant Filtering Common
Spatial Pattern,DFBCSP)等 多 种 CSP 的 修 改 版 本 。
Koles 于 1991 年提出使用基于小波变换的 CSP 来提取
EEG 信号中的非正常成分[4],紧接着文献[5]首次将该方
法运用到运动想象 EEG 信号分类中,文献[6]采用 CSP
方法对多导 EEG 信号进行了特征提取,文献[7]在此基
础上提出的 sub-CSP 算法对运动想象脑电信号的分类
识别进行了进一步研究。
然而,传统的 CSP 需要大量的输入通道,缺乏频率
信息 [1]。本文提出 EMD 结合 CSP 的特征提取算法,在
CSP 的基础上加入 EMD 的频域信息,为更好地进行想
象运动判别提供了一个可选的途径。主要分几个步骤:
(1)将滤波后的信号经过经验模式分解得到各固有模态
函数(IMF);(2)通过计算每个固有模态函数的能量谱
筛选在 μ 节律和 β 节律(5~28 Hz)的振荡模态,将多个
IMF 看成新的多通道信号;(3)对多通道信号通过公共
空间模式分解,通过公共空间模式算法滤波后数据在不
同类别状态下的能量差异选取空间滤波器的构造成分,
构造新的空间滤波器,获取特征向量组;(4)采用支持向
量机(SVM)分类器对特征进行分类。
文章第 2 章对本文使用的 EMD 分解、CSP 特征提
取、频域能量分析、支持向量机等算法进行了简要的介
绍。第 3 章为数据集描述。第 4 章介绍了数据分析结
果:和以往只采用 CSP 进行特征提取相比较,不仅减少
了电极数,经过 EMD 分解后,再用改进的 CSP 滤波器提
取特征的方法对实验中 9 位受试的左右手想象运动均
获得了较好的分类准确度,9 位受试所有实验平均分类
正确率在 92%,最高达 93.8%。
2 方法
2.1 经验模式分解
经验模式分解[3]作为一种自适应的纯数字导向的分
解方法,主要过程是将信号分解为一系列完全的几乎正
交 的 独 立 成 分 ,通 常 把 该 成 分 称 为 固 有 模 态 函 数
(IMFs),一个固有模态函数代表嵌入在信号中的一种
简单振荡模式。所以,EMD 分解的基函数由信号本身
决定,具有很好的自适应性并且能够精确的表征信号[8]。
对于实信号 x(t) ,给出如下 EMD 算法:
(1)判断每个 x(t) 的局部极值,用三次样条曲线进
行曲线拟合,局部极大值形成上包络 emax(t) ,局部极小
值形成下包络 emin(t) 。
2
(1)
emin(t) + emax(t)
(2)求 emax(t) 和 emin(t) 的均值:
m(t) =
(3)计算输入信号 x(t) 和 m(t) 的差值:
c(t) = x(t) - m(t)
(2)
如果 c(t) 不能满足 IMF 的定义截至条件,重复上述
过程(1)~(3),否则,提取 c(t) 作为固有模态函数,剩余
量 r(t) 计算如下:
r(t) = x(t) - c(t)
(3)
(4)剩余量作为一个新的数据经过相同的筛选过程
以获得下一个更低频率的固有模态函数,直到剩余函数
r(t) 为一个单调函数或者仅有一个极致时,分解过程停
止。假设原始信号 x(t) 被分解为 n 个固有模态函数和
一个剩余函数量 r(t) ,可以得重构信号:
x(t) = ∑
n ci(t) + rn(t)
i = 1
(4)
2.2 公共空间模式提取
CSP 主要思想是通过线性变换得到投影矩阵,将多
通道 EEG 数据投影到低维子空间,投影矩阵每一行由信
道的权重因子组成,这种变幻可以使信号矩阵的方差最
大化。CSP方法是基于两类的协方差矩阵同时对角化[8-9]。
下面通过一个实验例子来详细说明这个算法 [9]。
XL 和 XR 表示预处理过的 EEG 矩阵(下标 L 和 R 分别
代表左右想象运动和右手想象运动),N × T 维,N 为
通道数,T 为采样点数。标准化的空间协方差表示为:
RH =
XL X T
L
trace(XL X T
L )
(5)
计算机工程与应用www.ceaj.org
张学军,黄婉露,黄丽亚,等:公共空间模式算法结合经验模式分解的 EEG 特征提取
2017,53(13)
11
RF =
XR X T
R
trace(XR X T
R )
0
0
RR = U0∑U T
对复合协方差矩阵进行对角化分解:
RL + - ---
R = - ---
白化矩阵:
P = ∑-1/2U T
对平均协方差矩阵进行变形:
SL = P
SR = P
SL 和 SR 有共同的特征向量,满足下式:
SL = B∑
SR = B∑
∑
+∑
- ---
RL P T
- ---
RR P T
BT
R
= I
BT
L
(6)
(7)
(8)
(9)
(10)
(11)
(12)
L
R
(13)
SL 最大特征值对应的特征向量对应 SR 最小特征值,反
之亦然。对应白化 EEG 的最大特征值的特征向量进行
变换,可以获得两个信号矩阵的最优分离方差。投影矩
阵 W 表示为:
W = U T P
(14)
此处 W 即空间滤波器,原始的 EEG 信号 x(t) 经过 W 可
以被投影得到新的数据集 Z0 ,如式(13)所示:
Z0 = x( )t ∙W
构建矩阵 Z = [
2 行和后 m 行构成,特征向量 f =[
ö
ø÷
(15)
z1,z2,…,z2m ∈ R N × 2m 由 Z0 的前 m
]
T ∈
f1,f2,…,f2m
]
,i = 1,2,…,2m
(16)
m < T
æ
èç
R2m × 1 即定义为:
)zi
var(
∑
2m var(
)zj
fi =
j = 1
2.3 频域能量分析
脑电信号很多特征表现为频域变化,傅里叶变换可
以实现信号从时域到频域的转化,进而分析脑电信号,
具体变换过程如下:
x( )j ω
X ( )k = ∑
(17)
k - 1
N
)
)
(
j - 1 (
N
j = 1
其中,x( j) 为原始信号;X(k) 为傅里叶变换后的频率信
号;ω N = e(
-2πi N 为 N 阶傅里叶级数。
)
傅里叶变换后的信号以复数形式表示,为了便于分
析,做共轭相乘变换,如式(16)所示:
X ( )k ′ = X ( )k × -
变换后不会改变信号的能量变化,便于分析[10]。
-- -----
X ( )k ′ 1 000
(18)
2.4 支持向量机(SVM)
支持向量机(SVM)[11]是基于统计学习理论的一种
机器学习方法,它通过适当的非线性映射将输入向量映
射到一个高维的特征空间,使得数据(属于两类)总能被
一个超平面分割。所谓最优分类面就是要求分类面不
但能将两类数据正确分开,而且使分类间隔最大[12]。
对于非线性问题,只需要将它通过非线性变换转化
为另一个空间中的线性问题,然后再构造最优分类面。
核参数 γ 和误差惩罚因子 C 是影响 SVM 性能的主要参
数;γ 的取值影响空间变换后的数据分布,惩罚因子 C
则决定了支持向量机的收敛速度及推广能力。因此,对
γ 和 C 的选择很大程度上影响了脑电信号的识别率。
本文采用网格化交叉验证方法进行 γ 和 C 最优参
数的选择,将训练集作为原始的数据集,在一定范围内
改变核函数和惩罚因子的值,运用交叉验证方法进行分
类,选择分类准确率最高的 γ 和 C 作为最佳参数。
3 数据集描述
实 验 数 据 来 自 BCI Competition 2008 Data Sets
2b 数据[13]。数据包含 9 个被试的脑电 EEG 数据。被试
都是右利手,有正常或者矫正正常视力。所有被试坐在
扶手椅上,注视距离眼部 1 米的屏幕[13]。
如图 1所示,每个被试采集数据过程分为 5段,实验采
集所有 5 段数据但仅对后 3 段数据进行记录反馈。每
段包含一些组(run),每段开始,记录大约 5 分钟估计眼电
EOG 影响。5 分钟具体为:(1)2 分钟睁眼(注视屏幕上固
定的“+”号);(2)1分钟闭眼;(3)1分钟眼动;(4)1分钟休息。
EOG
EOG
EOG
Eyes open
Eyes closed
Movement
Run 1
…
Run N
图 1 EEG 信号采集流程图
实验分 3 个电极记录(C3、Cz 和 C4)的脑电信号,
采样频率为 250 Hz。数字信号通过 0.5~100 Hz 的带通
滤波,50 Hz 时提供陷波滤波器。电极在各个被试的头
部位置如图 2 所示。电极在 Fz 的位置设置为 EEG 接
地。除了 EEG 的通道,用 3 个电极收录眼电数据(如图
2,左边电极作为参考量)。
#2
#3
#1
图 2 想象运动人脑电极分布图
每一个被试在两周内的不同天,在不记录反馈的情
况下,参加了 2 期(session)只有屏幕显示的测试。每一
期包括 6 组(run),一组包含 10 条(trial),两类动作想
象。一共可以获得每组 20 条(trial)数据,一期共有 120
条数据。总计可获得每一个人的每一个动作图像类别
的 120 次重复的数据。在首个运动想象训练之前,被试
执行和想象不同的身体部位的动作,选择想象效果最好
的动作(比如压一个球或者拉闸)。
如图 3 所示,屏幕显示范式包括两类:一种是左手
计算机工程与应用www.ceaj.org
12
2017,53(13)
Computer Engineering and Applications 计算机工程与应用
(类别 1)运动想象;一种右手(类别 2)运动想象。每一
实验从固定的“+”开始,加短期声音提醒发音(1 kHz,
70 ms)。若干秒后,一个视觉提示以箭头的形式在 3~
4.25 s 期间呈现,提示时长 1.25 s;之后被试想象相应的
手臂动作持续 4 s;每一条实验完成后至少休息 1.5 s;加
入 1 s 内任意时间来打断被试的适应性。
Screening
beep
Imagination of
left hand movement
Imagination of
right hand movement
Fixation cross Cue Imagery Period Pause
0
8
9
图 3 实验界面和单次实验时间序列图
2
3
4
5
6
7
1
V
μ
/
幅
振
30
20
10
0
-10
-20
0
1
2
3
4
时间/s
5
6
7
8
(a)受试 B01C3 电极第一个 trial 滤波前 EEG 信号
V
μ
/
幅
振
30
20
10
0
-10
-20
0
1
2
3
4
时间/s
5
6
7
8
time in s
(b)受试 B01C3 电极第一个 trial 滤波后 EEG 信号
图 5 受试 B01C3 电极第一个 trial 滤波前后对比图
4 结果
4.1 预处理
本文采用 9 位受试的 EEG 数据。首先,对原始信号
进行基线校正;其次,使用独立分量分析(Independent
Component Analysis,ICA)去 除 原 始 信 号 中 的 眼 电 伪
迹,利用 Matlab 工具箱 EEGLAB 进行 ICA 脑电伪迹去
除[14]。图 4 是用 EEGLAB 对受试 01 的脑电信号进行 ICA
处理后,得到的 C3、C4 以及 Cz 通道的独立成分图。
2
1
3
+
(a)ICA 处理后 C3、C4、Cz 通道独立成分能量分布脑
-
GDF file
地形图(1∶Cz,2∶C3,3∶C4)
Channel 3
Continous Data
100
s
l
a
i
r
T
1000
500
1500
2000
2500
3000
Frames
222
111
0
-111
-222
Channel 4
Continous Data
100
s
l
a
i
r
T
500
1000
1500
2000
2500
3000
Frames
222
111
0
-111
-222
Activity power spectrum
20
10
0
-10
5
10 15 20 25 30 35 40 45 50
g
l
·
0
1
)
z
H
/
2
r
e
w
o
P
V
μ
(
g
l
·
0
1
)
z
H
/
2
r
e
w
o
P
V
μ
(
30
20
10
0
-10
Activity power spectrum
5
10 15 20 25 30 35 40 45 50
Frequency/Hz
Frequency/Hz
(b)处理后 C3、C4 通道信号独立源成分功率
谱图(左:C3 ,右:C4)
图 4
ICA 处理后 C3 、C4 、Cz 通道独立成分图
下文主要对和运动想象相关的频率段进行分析,即
μ 节律和 β 节律,所以随后进行 5~28 Hz 的带通滤波,如
图 5、图 6 所示,为去伪迹后,滤波前后的对比图。
4.2 EMD 分解与频段筛选
V
μ
/
幅
振
30
20
10
0
-10
-20
0
1
2
3
4
时间/s
5
6
7
8
(a)受试 B01C4 电极第一个 trial 滤波前 EEG 信号
V
μ
/
幅
振
30
20
10
0
-10
-20
0
1
2
3
4
时间/s
5
6
7
8
(b)受试 B01C4 电极第一个 trial 滤波后 EEG 信号
图 6 受试 B01C4 电极第一个 trial 滤波前后对比图
一个 trial(0~8 s)的 EEG 信号进行 EMD 分解,如图 7 所
示,图(a)为左手想象运动信号 EMD 分解结果,图(b)为
右手想象运动信号 EMD 分解结果。
每一个 IMF 分量代表一种振荡模式,包含了相应的
频域信息。选取前 8 阶 IMF 分量,求能量谱。图 8 和图 9
分别为左手进行想象运动 C3 通道、C4 通道前 8 阶 IMF
分量的信号能谱图;图 10 和图 11 分别为右手进行想象
运动 C3、C4 通道前 8 阶 IMF 分量的信号能谱图。可以
看出,在 5~28 Hz 内,IMF1、IMF2、IMF3 的能量较为集
中 IMF4~IMF8 则可能为低频高能量的肌电信号。由
此,选取 IMF1~IMF3 做进一步 CSP 特征提取。
4.3 改进和 CSP 滤波
将一次实验的 C3、C4 通道前 3 阶 IMF 分量合并,构
成 一 个 6 × 2 000 的 矩 阵 X i(i = L 表 示 想 象 左 手 运 动 ,
i = R 表示想象右手运动),其中,6 为 IMF 个数,可看作
通道数,2 000 为一次实验的采样点个数,即窗口长度。
整个实验过程包含 120 组实验,共得到 120 组向量矩阵,
分为 40 组测试向量矩阵和 80 组训练向量矩阵,分别进
行公共空间模式分解。具体过程如下:
选取一位受试,以受试 01T 为例。对该受试滤波后
(1)使用 2.2 节中所述的 CSP 特征提取算法进行计
计算机工程与应用www.ceaj.org
张学军,黄婉露,黄丽亚,等:公共空间模式算法结合经验模式分解的 EEG 特征提取
2017,53(13)
13
Empirical Mode Decomposition,ϕ/2π = 0
Empirical Mode Decomposition,ϕ/2π = 0
l
a
n
g
i
S
1
f
m
i
2
f
m
i
3
f
m
i
4
f
m
i
5
f
m
i
6
f
m
i
7
f
m
i
i
8
f
m
0.04
0.02
0
.
s
e
r
l
a
n
g
i
S
1
f
m
i
2
f
m
i
3
f
m
i
4
f
m
i
5
f
m
i
6
f
m
i
7
f
m
i
i
8
f
m
0
-0.02
-0.04
-0.06
0
.
s
e
r
l
a
n
g
i
S
1
f
m
i
2
f
m
i
3
f
m
i
4
f
m
i
5
f
m
i
6
f
m
i
7
f
m
i
8
f
m
i
i
9
f
m
0.04
0.02
.
s
e
r
1
2
3
4
5
6
7
8
0
1
2
3
4
5
6
7
8
时间/s
时间/s
(a)受试 B01 左手想象运动信号 EMD 分解结果(左边∶C3 电极,右边∶C4 电极)
Empirical Mode Decomposition,ϕ/2π = 0
Empirical Mode Decomposition,ϕ/2π = 0
l
a
n
g
i
S
1
f
m
i
2
f
m
i
3
f
m
i
4
f
m
i
5
f
m
i
6
f
m
i
7
f
m
i
8
f
m
i
i
9
f
m
0.06
0.04
0.02
0-0.02
.
s
e
r
1
2
3
4
时间/s
5
6
7
8
0
1
2
3
4
时间/s
5
6
7
8
(b)受试 B01 右手想象运动信号 EMD 分解结果(左边:C3 电极,右边:C4 电极)
图 7 受试 B01 左右手想象运动信号 EMD 分解图
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 0.4
0.2
0
μ
/
y
g
r
e
n
e
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF1 frequency/Hz
V 0.2
0.1
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF3 frequency/Hz
V 0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF5 frequency/Hz
10 15 20 25 30
5
IMF2 frequency/Hz
10 15 20 25 30
5
IMF4 frequency/Hz
10 15 20 25 30
5
IMF6 frequency/Hz
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 0.4
0.2
0
μ
/
y
g
r
e
n
e
V 0.2
0.1
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF1 frequency/Hz
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF3 frequency/Hz
V 0.4
0.2
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF5 frequency/Hz
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF2 frequency/Hz
10 15 20 25 30
5
IMF4 frequency/Hz
10 15 20 25 30
5
IMF6 frequency/Hz
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
V 0.2
0.1
0
μ
/
y
g
r
e
n
e
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF7 frequency/Hz
5
IMF8 frequency/Hz
图 8 左手进行想象运动 C3 通道前 8 阶 IMF 分量能谱图
10 15 20 25 30
10 15 20 25 30
5
IMF7 frequency/Hz
5
IMF8 frequency/Hz
图 9 左手进行想象运动 C4 通道前 8 阶 IMF 分量能谱图
10 15 20 25 30
算,选取所有的特征向量组成空间滤波器 W 。
(2)记原始输入的80组6×2 000的矩阵数据为Data1,
对 Data1 用上述滤波器 W 进行滤波,滤波后的数据为
Data2。
(3)计算 Data2 中 6 个通道( C3 电极前 3 阶 IMF 分
量和 C4 电极前 3 阶 IMF 分量)数据的能量,选取每组训
练集中不同类别之间能量差异最明显的前 3 对分量对
应的特征向量,构成新的滤波器 W ′ 。
计算机工程与应用www.ceaj.org
14
2017,53(13)
Computer Engineering and Applications 计算机工程与应用
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
V 0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF1 frequency/Hz
V 0.4
0.2
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF3 frequency/Hz
V 0.2
0.1
0
μ
/
y
g
r
e
n
e
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF5 frequency/Hz
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF7 frequency/Hz
10 15 20 25 30
5
IMF2 frequency/Hz
10 15 20 25 30
5
IMF4 frequency/Hz
10 15 20 25 30
5
IMF6 frequency/Hz
10 15 20 25 30
5
IMF8 frequency/Hz
图 10 右手进行想象运动 C3 通道前 8 阶 IMF 分量能谱图
V 0.4
0.2
0
μ
/
y
g
r
e
n
e
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF1 frequency/Hz
10 15 20 25 30
5
IMF3 frequency/Hz
V 1.0
0.5
0
μ
/
y
g
r
e
n
e
V 0.10
0.05
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF2 frequency/Hz
10 15 20 25 30
5
IMF4 frequency/Hz
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF5 frequency/Hz
V 0.06
0.04
0.02
0
μ
/
y
g
r
e
n
e
10 15 20 25 30
5
IMF7 frequency/Hz
10 15 20 25 30
5
IMF6 frequency/Hz
10 15 20 25 30
5
IMF8 frequency/Hz
图 11 右手进行想象运动 C4 通道前 8 阶 IMF 分量能谱图
(4)40 组测试数据作为 Data3,使用改进后的滤波器
W ′ 进行滤波,按照公式(16)进行特征提取。
(5)最后通过 SVM 进行分类。
经过对 9 位受试的 27 段 session 的数据分析发现,在
5~30 Hz 频段中,每组训练集中总存在不同类别之间能
量相差明显的前 3 组分量,如表 1 所示。以第一位受试
单次实验训练矩阵的各阶 IMF 能量为例,C3 通道的
IMF1 分量,C4 通道的 IMF2、IMF3 分量的两个类别能
量差异明显,故可以选取此 3 个分量对应的特征向量添
加到新的滤波器中。
表 1 受试 B01 单次实验用于构造新 CSP 滤波器的训练矩阵
的各阶 IMF 能量最大值
能量类别
想象左手运动
最大能量/μV
想象右手运动
最大能量/μV
IMF1
IMF2
IMF3
C3
C4
C3
C4
C3
C4
0.68
0.65
0.44
0.57
0.52
0.20
5.80×10-17
0.50
0.61
1.18×10-17
0.58
2.50×10-18
4.4 结果与比较
将每个 trial 的 C3、C4 通道前三阶 IMF 分量合并构
成观察向量,经过改进后的 CSP 进行空域滤波并提取特
征,为进一步验证该特征的有效性,图 12(a)、(b)、(c)为
受试 01 左右手想象运动时,单次 trial 的脑电信号 EMD
分解后前 3 阶 IMF 分量经过 CSP 空域滤波后特征的脑
地形图。由图可见,在进行两类想象运动时,C3、C4 电
极能量相差较大,同时,想象左手运动和时 C4 电极明显
比 C3 电极活跃,想象右手运动时则刚好相反,表现为
C3 电极比 C4 活跃,这与运动想象中时间按相关同步/
去同步特性刚好相符。由此说明,该特征可以有效表征
出“源”的空间位置信息,可以作为分类特征。
3 000
2 000
1 000
0
-1 000
-2 000
-3 000
(a)想象运动第 1 阶 IMF 经 CSP 滤波后特征能量的脑地
2 000
1 500
1 000
500
0
-500
-1 000
-1 500
-2 000
形图(左图:想象左手运动,右图:想象右手运动)
2 000
1 500
1 000
500
0
-500
-1 000
-1 500
-2 000
(b)想象运动第 2 阶 IMF 经 CSP 滤波后特征能量的脑地
3 000
2 000
1 000
0
-1 000
-2 000
-3 000
形图(左图:想象左手运动,右图:想象右手运动)
2 000
1 500
1 000
500
0
-500
-1 000
-1 500
-2 000
(c)想象运动第 3 阶 IMF 经 CSP 滤波后特征能量的脑地
3 000
2 000
1 000
0
-1 000
-2 000
-3 000
形图(左图:想象左手运动,右图:想象右手运动)
图 12 想象运动前 3 阶 IMF 经 CSP 滤波后特征能量的脑地形图
根据 CSP 特征提取的信号,采用 SVM 分类器进行
分类。得到图 13 的 9 位受试的分类结果,首先,图 13(a)、
(b)给出了 9 位受试单次实验(包含 120 组左右手想象运
动)的分类正确率。由图可见,除了个别较低的情况外,
分类正确率基本在 85%~95%之间浮动,且 90%~95%区
间内居多。如图 13(c)所示,9 位受试想像运动的所有实
验平均分类准确度在 88%~94%之间浮动,其中,受试 07
达到最高 93.8%;9 位受试所有实验平均分类准确度为
92%,可见本文方法具有较为有效且稳定的分类正确率。
其次,将 EMD-CSP 方法与其他方法比较,如表 2 所
示,列出了第三届 BCI 竞赛中成绩较好的 3 组,可见采
用相同数据集的情况下,本文提出的方法的分类准确度
高于其他 3 组。表 3 给出了 9 位受试的所有实验的分类
计算机工程与应用www.ceaj.orgFc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -2000-1500-1000-5000500100015002000Fc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -3000-2000-10000100020003000Fc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -2000-1500-1000-5000500100015002000Fc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -3000-2000-10000100020003000Fc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -2000-1500-1000-5000500100015002000Fc5 Fc3 Fc1 Fcz Fc2 Fc4 Fc6 C5 C3 C1 Cz C2 C4 C6 Cp5 Cp3 Cp1 Cpz Cp2 Cp4 Cp6 Fp1 Fpz Fp2 Af7 Af3 Afz Af4 Af8 F7 F5 F3 F1 Fz F2 F4 F6 F8 Ft7 Ft8 T7 T8 T9 T10 Tp7 Tp8 P7 P5 P3 P1 Pz P2 P4 P6 P8 Po7 Po3 Poz Po4 Po8 O1 Oz O2 Iz -3000-2000-10000100020003000
张学军,黄婉露,黄丽亚,等:公共空间模式算法结合经验模式分解的 EEG 特征提取
2017,53(13)
15
表 2 第三届 BCI 竞赛成绩和本文方法结果的比较
姓名
特征提取方法
所用通道
Wei Qingguo[15]
CSSD+waveformmean+FDA
AR+spectralpower+wavelet coefficients
Hammon[16]
Blankertz[17]
Offset+spectral power
Not Given
Logistic regression
SVM
Author in this paper
EMD+CSP
B01
B02
B03
B04
10
20
30
实验序数
40
50
60
(a)B01~B04 受试单次实验分类正确率
B05
B06
B07
B08
B09
10
20
30
实验序数
40
50
60
(b)B05~B09 受试单次实验分类正确率
10
20
30
实验序数
40
50
60
100
95
90
85
80
75
70
0
100
95
90
85
80
75
70
0
100
95
90
85
80
75
70
0
/
%
率
确
正
类
分
/
%
率
确
正
类
分
/
%
率
确
正
类
分
均
平
(c)9 位受试的所有实验平均分类正确率
图 13
9 位受试分类正确率
正确率以及建模时间。本文仅采用了 6 个通道,去除 3
个眼电通道,仅提取 C3、C4 通道的信号进行分析,在保
证一定分类准确度的情况下,大大减少了通道数。另
外,如图 14 所示,9 位受试的所有实验的平均分类响应
时间在 0.17~0.21 s 之间,为便携式 BCI 系统的信号在线
采集提供了可行性。
表 3 B01~B09 位受试平均分类正确率和建模时间
受试编号
平均分类正确率/%
建模时间/s
B01
B02
B03
B04
B05
B06
B07
B08
B09
总平均
92.8
93.5
90.0
93.4
91.4
90.4
93.8
89.7
92.3
92.0
0.19
0.18
0.19
0.17
0.19
0.21
0.19
0.21
0.18
0.19
分类方法
SVM
Regularized logistic regression
64
33
结果/%
91
87
86
92
6
s
/
间
时
应
响
均
平
0.22
0.21
0.20
0.19
0.18
0.17
0.16
1
2
3
4
5
编号
6
7
8
9
图 14
9 位受试所有实验平均响应时间
5 讨论与总结
针对公共空间模式分解算法需要大量通道信号和
缺乏频域信息的缺点进行改进,本文提出 EMD-CSP 算
法。先将信号进行经验模式分解得到多个固有模态函
数,选取符合 μ 节律和 β 节律的振荡模式组合成多通道
信息簇,基于由公共空间模式滤波后的信号能量选取成
分构造新的空间滤波器提取特征,得到特征向量集,经
支持向量机分类后得到所有 9 位受试平均分类正确率
为 92%,其中,最高的受试达到 93.8%。改进的 CSP 滤
波成分选择方法在保证分类准确度的基础上大大减少
了实际使用的通道数量,为便携式脑机接口设计提供了
简单可行的方法。
参考文献:
[1] Yang Banghua,Li Huarong,Wang Qian,et al.Subject-based
feature extraction by using fisher WPD-CSP in brain-
computer interfaces[J].Computer Methods and Programs in
Biomedicine,2016,129:21-28.
[2] Aydemir O.Common spatial pattern-based feature extraction
from the best time segment of BCI data[J].Turkish Journal
of Electrical Engineering&Computer Sciences,2016,24:
3976-3986.
[3] Li Mingai,Cui Yan,Hao Dongmei,et al.An adaptive fea-
ture extraction method in BCI- bsed rehabilitation[J].Journal
of Intelligent & Fuzzy Systems,2015,28:525-535.
[4] Koles Z J.The quantitative extinction and topographic
mapping of the abnomnal components in the clinical
EEG[J].Electroencephalography and Clinical Neurophysiology,
1991,79(6):440-447.
[5] Muller Gerking J,Pfurtscheller G,Flyvbjerg H.Dasigning
optimal spatial filters for single- trial EEG classification
task[J].Clinical Neurophysiology,1999,
in a movement
110(5):787-798.
(下转 54 页)
计算机工程与应用www.ceaj.org0102030405060707580859095100试验序数单次试验分类正确率(%)123456789859095100被试编号1234567890.160.170.180.190.20.210.22被试编号