2019, Vol. 40
2. 东北大学 医学与生物信息工程学院, 辽宁 沈阳 110169
2. School of Medicine & Biological Information Engineering, Northeastern University, Shenyang 110169, China
3D显示带来的强烈感和沉浸感是2D影像无法比拟的, 3D时代已经到来.然而, 观看立体显示会给观看者带来一些不良反应:头痛、恶心和视觉疲劳等.因此, 评估观看者舒适与否在立体显示技术的发展中尤为重要.
立体深度运动可引发频繁的眼动效应, 眼动效应和视觉不舒适有密切联系.Yano等[1]发现, 立体深度运动是导致立体显示不舒适的重要原因之一.Jung等[2]利用人眼注意力机制, 提出了可衡量立体深度运动的速度与舒适度关系的模型.Speranza等[3]通过主观实验得到了视差变化率严重影响视觉舒适度这一结论.Shen等[4]通过脑电实验得到了与立体深度匀变速运动有关的三种ERP成分, 进而得到人们能够感知的运动速度阈值.
脑电(electroencephalography, EEG)是大脑内细胞活动在脑皮层的综合反映, 反应人们生理心理活动[3], 有非入侵、实时性等优势.但EEG极其微弱且多噪, 因此提取与任务相关的有用信息至关重要.
频域和空域分析法常用于分析EEG:小波变换、傅里叶变换、熵及共空间模式(common spatial pattern, CSP)[5]等.Huang等[6]提出一种自适应的时频分析方法——经验模式分解(empirical mode decomposition, EMD), 该方法处理单通道数据且模式混叠现象严重.Rehman等[7]提出了多变量经验模式分解(multivariate EMD, MEMD), 可处理多维数据, 避免分解出的本征模态函数(intrinsic mode function, IMF)数量频率不匹配等问题, 但仍存在模式混叠问题.Rehman等[8]基于模式混叠问题提出了噪声辅助经验模式分解(noise-assisted MEMD, NA-MEMD), 能更好地解决这些问题.
选择包含有用信息的IMF分量也是应用模式分解方法中的重要一步.文献[9]根据研究对象相关频率的先验信息对IMF进行选取.韩笑等[10]提出了一种不依赖于先验知识的方法, 基于互信息构建了一种能够衡量信号与IMF分量之间相似性的有效因子, 根据相似性决定哪些IMF是有用成分.
本文提出一种基于PMEMD-2JSD-CSP模型的特征提取和分类方法.首先对EEG进行通道筛选, 采用部分噪声辅助多变量经验模式分解(partial noise assisted multivariate EMD, PNA-MEMD)对EEG进行分解得到多尺度IMF, 使用基于詹森-香农散度(Jensen-Shannon distance, JSD)的有效因子和本文提出的筛选准则筛选IMF有用成分.最后, 利用CSP算法提取其特征, 支持向量机(support vector machines, SVM)对所提特征进行分类, 结果验证了该方法的有效性和采集数据的可分性.
1 实验方案通过主观实验找出运动场景中目标运动加速度与舒适度的关系, 得到可用于脑电实验的准确材料.
1.1 主观实验实验素材由Maya制作.刺激场景中运动目标为深灰色球体, 在深度方向上做周期性匀加速运动, 加速度a=0, 1, 2, 3, 4, 5, 6 m/s2, 平均速度为2 m/s, 单次运动时长为0.6 s.运动模式见表 1.整个实验包括4个子部分, 每部分之间有3 min休息.每部分包括70个试次, 试次时长为6 s, 其中7种运动模式各10次, 且采用随机呈现的方式, 每一试次呈现后被试进行舒适度判断[11], 判断之后进入下一试次.
|
表 1 主观实验的运动模式 Table 1 Movement patterns of subjective experiment |
模式M1~M7条件下主观判断分数:4.75, 4.52, 3.98, 3.19, 2.41, 1.64, 1.23.可以看出立体视频做深度运动及加速度超过4 m/s2时, 将引起被试的不舒适感; 加速度低于2 m/s2时, 则为舒适.因此选择加速度为1, 1.5 m/s2的运动模式作为脑电实验中标签为舒适类别的实验素材, 选择加速度为5.5, 6 m/s2的运动模式作为标签为不舒适类别的实验素材.
1.2 脑电实验实验包括10个子部分, 子部分之间休息2 min.每个子部分包括20个试次, 试次时长为4.2 s, 子部分中运动模式数量均匀且随机呈现, 试次之间有2 s停顿.实验流程如图 1所示.
|
图 1 实验流程 Fig.1 Flow diagram of experiments |
脑电数据采集采用国际10-20系统, 嵌入式34导脑电帽与Neuroscan系统.参考电极使用左耳电极A1, 采样频率为1 kHz.对EEG进行预处理:1~30 Hz的带通滤波及利用独立成分分析(ICA).本文构建的EEG数据库包括200个试次(训练集:150, 测试集:50), 每个试次时长4.2 s, EEG信号结构为30×4 200.
EEG信号包含了丰富的频率, 其中δ(0.5~4 Hz), θ(4~8 Hz), α(8~13 Hz), β(13~30 Hz), γ(30~44 Hz)与认知相关.在视觉认知领域, 一般只对θ, α, β三个波段进行研究.
2 算法本文提出的方法如图 2所示.首先对EEG进行导联选择, 采用PNA-MEMD方法得到多尺度IMF, 通过所提出的基于JSD有效因子的筛选准则完成IMF的筛选, 按权重叠加得到重构信号, 利用CSP得到重构信号的特征, 送入SVM进行分类.
|
图 2 算法流程 Fig.2 Process of the proposed algorithm |
EMD适用于非平稳非线性信号的时频分析, 可将信号分解为多尺度IMF和残差信号.PNA-MEMD利用高频带限噪声代替NA-MEMD中的白噪声, 分离了干扰成分且对残留信号的分解不会造成干扰, 由于添加的白噪声去掉了很多无用频率成分, 故IMF成分数量相比NA-MEMD明显减少[12].
PNA-MEMD算法的具体步骤如下:
1) 生成m通道的随机高斯白噪声信号W(t)={w1(t), w2(t), …, wm(t)}, W(t)中信号长度均为T;
2) 对m通道白噪声进行MEMD[7]分解, 分解后得到的第二个IMF分量IMF2作为PNA-MEMD叠加的噪声, 记为S(t)={s1(t), s2(t), …, sm(t)};
3) 将步骤2)中产生的m通道噪声S(t)与N通道EEG输入X(t)={x1(t), x2(t), …, xN(t)}结合, 生成信号Y(t), 其通道数为p=m+N;
4) 对Y(t)进行MEMD分解.
2.2 基于2JSD的IMF自适应选择使用基于JS散度的有效因子[13], 结合本文提出的筛选准则, 完成IMF的自适应选择.
JS散度和相对熵用于衡量相同事件空间里两个概率分布的差异情况, JS散度是相对熵对称化变形.设两个随机离散变量X, Y, 概率分布分别为P, Q.相对熵和JS散度公式为
|
(1) |
|
(2) |
式中:D(P||Q)和Djs(P||Q)分别表示X, Y之间的相对熵和JS散度; X, Y分别有n种不同取值, 每一个取值对应的概率为pi及qi, i=1, 2, …, n.
EEG信号经过分解具有多尺度IMF成分, 在每个尺度上考虑两种JS散度值:EEG与其多尺度IMF之间的JS散度值; EEG的IMF与噪声的IMF之间的JS散度值[13].m通道噪声与N通道EEG结合后的信号为Y(t)={x1(t), x2(t), …, xN(t), s1(t), s2(t), …, sm(t)}.分解后j尺度的IMF表示为YFj(t)={xf1j(t), xf2j(t), …, xfNj(t), sf1j(t), sf2j(t), …, sfmj(t)}, 信号长度均为L.
有效指标由两部分构成:
1) Dijx代表EEG i通道的j尺度IMF分量xfij(t)和EEG的i通道的原始信号xi(t)之间的JS散度值为
|
(3) |
2) Dikjxs代表了EEGi通道的j尺度IMF分量xfij(t)和噪声k通道的j尺度IMF分量sfkj(t)之间的JS散度值:
|
(4) |
有效因子为
|
(5) |
式中:αj为最终用于筛选IMF成分的有效指标, ej越小, esj越大, 则αj越小, 对应的j尺度IMF分量包含的有用信息越多.
本文提出的自适应筛选准则具体步骤如下:
1) 原始EEG信号经过自适应通道选择及模式分解, 得到各通道相同数量的IMF分量;
2) 去掉频带过高的IMF分量之后, 将IMF分量分为两个集合, 集合H包含大量α, β频段, 集合L包含少量α, θ及更低频段;
3) 根据式(3)~(5)计算IMF分量的有效指标αj, 在集合H和L中分别进行升序排序;
4) 依据最佳分类准确率, IMF选择个数为3.取集合H升序排序后前两个IMF分量和集合L升序排序后第一个IMF分量;
5) 对3个IMF分量按权重相加得重构信号.
本文的筛选准则不仅考虑了α, β频段, 也考虑了θ频段等低频频段对于任务的轻微影响.该筛选准则结合排序的机制, 对于不同被试可更加清晰地得知任务对其EEG信号的何种频带影响更大.
3 结果与分析 3.1 通道选择结果多通道EEG信号会带来噪声和冗余, 故通道优化对于后续信号处理是非常必要的.本文利用基于CSP权重的l1范数方法进行通道选择[14], 得到被试的通道排序, 选择前n个通道组合, 通过CSP和SVM得到30个分类率, 选择使分类率最大的通道组合进行后续处理.被试的通道组合结果见表 2.
|
|
表 2 不同被试的最优通道组合及对应的最大分类率 Table 2 Optimal channel combination and corresponding maximum classification rate of different subjects |
对于5个被试, 分别选择了11, 14, 14, 14, 11个导联.表 2中A1表示通道选择后的分类率, A2表示原始分类率, 5名被试的A1均大于相应的A2.图 3为被试平均的频域r2脑电地形图, r2值探究两种标签下各通道的频带频谱变化.r2值越高, 表示对应通道下两类试次相应的频段越可分.表 2中出现了2, 3, 4次的导联基本均在图 3中变化明显, 这证明了基于CSP权重选择通道的可行性.
|
图 3 被试在8~18 Hz平均频域r2的脑电地形图 Fig.3 Average topographic map of r2 EEG between 8~18 Hz |
使用4通道高斯白噪声, 长度为4 200点, 均值为0, 方差为输入信号方差的10%.将白噪声进行MEMD分解, 得到的4通道IMF2作为独立噪声与EEG结合得到混合信号, 使用MEMD对其进行分解.模态分解会分解出时间尺度小的高频信号, 某些IMF存在着大量超过30 Hz的无用信号, 计算平均频率用于剔除无用IMF.设时域信号为x(t), x(t)的频谱为F(f), 则x(t)的平均频率MF为
|
(6) |
对5名被试每个尺度的平均IMF求平均频率, IMF1和IMF2的MF值分别达到了27.297 5, 29.575 1 Hz, 故后续分析中剔除IMF1, IMF2.
为了探究IMF分量包含哪些频带, 需要观察IMF分量的边际谱.希尔伯特-黄变换(Hilbert-Huang transform, HHT)是一种用于处理非平稳、非线性信号的方法.利用模态分解方法分解出原信号的IMF, 再对IMF进行HHT得到其瞬时频率和瞬时幅度, 进而得到IMF分量的边际谱.图 4为被试S1的IMF3~6的平均边际谱.观察被试S1分解后的IMF3~6边际谱, IMF3几乎包含整个β频段; IMF4包含部分α和β频段; IMF5包含α和θ频段; IMF6包含θ及更低频段.由于视觉认知任务对EEG信号的α, β频段影响较大, 对θ频段影响较小[15-16], 故进行信号重构时要考虑包含α, β, θ频段的成分.
|
图 4 被试S1分解后的IMF3~6边际谱 Fig.4 Marginal spectrum of decomposed IMF3~6 for one subject S1 (a)—舒适试次；(b)—不舒适试次. |
视觉认知任务对3种频段的影响有着大致相似的走向, 但不同被试对于任务有着不同的感知, 故应依据被试自身数据来选择IMF成分用于信号重构.本文采用如下机制:对于每一试次的EEG信号进行PNA-MEMD分解得到多尺度IMF, 分别对IMF3~5和IMF6~10计算有效因子并升序排序得到序列D3~5和D6~10, 将序列D3~5中的前两阶IMF分量IMFh1, IMFh2, 和D6~10中的第一阶IMF分量IMFl1叠加作为重构信号, 公式为
|
(7) |
重构信号不仅包含α和β频段, 也包含了部分θ或更低频段, 由于θ频段对EEG信号的影响很小, 故引入权重机制.本文提出的筛选机制具有两个优点:考虑θ频带对EEG信号的轻微影响; 权重机制的引入减弱了偶尔引入的无用IMF对重构结果的影响.为方便显示, 只列出被试S1的IMF3~9的有效因子平均值及标准差, 如表 3所示.
|
|
表 3 有效因子的值(均值±标准差) Table 3 Effective factor value(mean±SD) |
CSP算法使用线性变换将多通道时域信号投影到低维空间的子空间中, 使信号的某类别方差最大化且另一类别方差最小化[5].SVM能较好地解决高维数、小样本等问题.因此本文采用CSP处理重构信号, SVM对所得特征进行分类.分类时, 选用RBF核, 使用网格搜索选择其最优参数c和r, 训练样本和测试样本按150:50的比例随机选取, 对测试集进行1 000次测试取平均值为最后结果.
本文提出的方法记作2JSD, 手动筛选IMF3~5和IMF3~4并叠加重构的方法分别记作IMF3+4+5和IMF3+4[9], 三种方法使用相同的数据和参数，分类结果见表 4.对于被试S2, IMF3+4+5取得最高分类率, 对于其他被试, 2JSD均取得了最高分类结果, 这表明个体差异影响实验结果.2JSD的平均分类率为69.324%, 为三种方法中的最高, 表明立体视频的匀加速深度运动舒适与否的可分性.
|
|
表 4 分类结果比较 Table 4 Comparison of classification results |
本文探究立体深度匀加速运动对于视觉舒适度的影响, 构建PMEMD-2JSD-CSP模型用于提取EEG信号特征.对EEG信号进行通道选择及PNA-MEMD分解后, 采用基于JSD的有效因子及本文提出的筛选准则对分解得到的IMF进行两次不同范围的筛选, 并在重构信号时结合权重机制, 得到了较为满意的结果.除了被试S2, 其余被试均在此模型下得到最优分类率, 最高值为73.16%, 平均值为69.324%, 此结果说明立体视频做匀加速深度运动舒适与否具有可分性.
| [1] |
Yano S, Emoto M, Mitsuhashi T.
Two factors in visual fatigue caused by stereoscopic HDTV images[J]. Displays, 2004, 25(4): 141–150.
DOI:10.1016/j.displa.2004.09.002 |
| [2] |
Jung Y J, Lee S I, Sohn H, et al.
Visual comfort assessment metric based on salient object motion[J]. Journal of Electronic Imaging, 2012, 21(1): 011008.
DOI:10.1117/1.JEI.21.1.011008 |
| [3] |
Speranza F, Tam W J, Renaud R, et al.Effect of disparity and motion on visual comfort of stereoscopic images[C]//Proceedings of SPIE.San Jose, 2006: 6055.
|
| [4] |
Shen L L, Du W R, Wang C.
Event-related potentials measurement of perception to 3D motion in depth[J]. China Communications, 2015, 12(5): 93–96.
DOI:10.1109/cc.2015.7112047 |
| [5] |
Park C, Looney D, Rehman N U, et al.
Classification of motor imagery BCI using multivariate empirical mode decomposition[J]. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2013, 21(1): 10–22.
DOI:10.1109/TNSRE.2012.2229296 |
| [6] |
Huang N E, Zheng S, Steven R, et al.
The empirical mode decomposition and the Hilbert spectrum for nonlinear non-stationary time series analysis[J]. Proceedings of the Royal Society A—Mathematical Physical and Engineering Sciences, 1998, 454(1971): 903–995.
DOI:10.1098/rspa.1998.0193 |
| [7] |
Rehman N, Mandic D P.
Multivariate empirical mode decomposition[J]. Proceedings of the Royal Society A-Mathematical Physical and Engineering Sciences, 2010, 466(2117): 1291–1302.
DOI:10.1098/rspa.2009.0502 |
| [8] |
Rehman N, Mandic D P.
Filter bank property of multivariate empirical mode decomposition[J]. IEEE Transactions on Signal Processing, 2011, 59(5): 2421–2426.
DOI:10.1109/TSP.2011.2106779 |
| [9] |
Gaur P, Pachori R B, Wang H, et al.A multivariate empirical mode decomposition based filtering for subject independent BCI[C]//The 27th Irish Signals and Systems Conference(ISSC).Derry: IEEE, 2016: 1-7.
|
| [10] |
韩笑, 佘青山, 高云园, 等.
基于NA-MEMD和互信息的脑电特征提取方法[J]. 传感技术学报, 2016, 29(8): 1140–1148.
( Han Xiao, She Qing-shan, Gao Yun-yuan, et al. Feature extraction of EEG based on NA-MEMD and mutual imformatin[J]. Chinese Journal of Sensors and Actuators, 2016, 29(8): 1140–1148. DOI:10.3969/j.issn.1004-1699.2016.08.003 ) |
| [11] |
Du S P, Masia B, Hu S M, et al.
A metric of visual comfort for stereoscopic motion[J]. ACM Transactions on Graphics, 2013, 32(6): 1–9.
|
| [12] |
Huang W P, Zeng J, Wang Z Y, et al.
Partial noise assisted multivariate EMD:an improved noise assisted method for multivariate signals decomposition[J]. Biomedical Signal Processing and Control, 2017, 36: 205–220.
DOI:10.1016/j.bspc.2017.04.003 |
| [13] |
She Q S, Ma Y L, Meng M, et al.
Noise assisted MEMD based relevant IMFs identification and EEG classification[J]. Journal of Central South University, 2017, 24(3): 599–608.
DOI:10.1007/s11771-017-3461-5 |
| [14] |
Li M A, Ma J Y, Jia S M.Optimal combination of channels selection based on common spatial pattern algorithm[C]//IEEE International Conference on Mechatronics and Automation.Beijing: IEEE, 2011: 295-300.
|
| [15] |
Chen C X, Li K, Wu Q Y, et al.
EEG-based detection and evaluation of fatigue caused by watching 3DTV[J]. Displays, 2013, 34(2): 81–88.
DOI:10.1016/j.displa.2013.01.002 |
| [16] |
Malik A S, Khairuddin R N, Amin H U, et al.
EEG based evaluation of stereoscopic 3D displays for viewer discomfort[J]. Biomedical Engineering Online, 2015, 14(12): 21.
|