2016, Vol. 37
2.大连海洋大学 信息工程学院，辽宁 大连 116023;
3.东北大学 机械工程与自动化学院，辽宁 沈阳 110819
2.College of Information Engineering, Dalian Ocean University, Dalian 116023, China;
3.School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
在脑电(electroencephalogram，EEG)信号的测量与分析中，头皮处得到的观测EEG是由大脑皮层产生的源信号经组织液、头骨等容积导体传输到头皮处产生的混合信号，并在传输期间混有工频干扰和高斯噪声等外部信源.描述信号系统混合过程的模型有瞬时混合和卷积混合两种，相应的信号盲分离方法为瞬时盲分离和盲反卷积.考虑到EEG传输的时延性及多向性，卷积混合模型更符合其物理意义，因此，EEG盲分离应采用盲反卷积.
目前大多数研究中的EEG盲分离步骤均采用瞬时盲分离[1-3]，采用盲反卷积的研究较少.Abutaleb等[4]使用随机微积分的方法对观测EEG信号进行盲反卷积，取得了较好效果，但该方法只能针对3个观测信号分离出2个源信号，即假设观测信号仅有2个源信号混合而成.Dyrholm等[5]提出一种基于极大似然估计的盲反卷积算法，对观测EEG进行盲分离.该方法可以分离出多个源信号，但需要针对所有观测EEG信号首先进行Informax ICA，然后在得到的分量中选取若干个分量进行盲反卷积; 其不足在于Informax ICA步骤会造成数据损失，且分量挑选是否合适直接影响最终分离结果.
针对上述情况，本文提出一种基于盲反卷积的EEG盲分离方法，结合脑电源成分的独立性确定代价函数，并采用共轭梯度法进行迭代寻优.利用EEG仿真实验进行方法验证，结果表明，本文方法可以有效进行观测EEG盲分离.
1 盲反卷积算法 1.1 EEG的混合模型描述设m维观测EEG信号xi(i=1, 2, …, m)是由n维源信号sj(j=1, 2, …, n)卷积混合而成，则混合模型可以表示如下：
|
(1) |
式中：t为采样点; hij(·)为第j个源信号传输到第i个观测点的卷积混合过程，用FIR滤波器表示; l为FIR滤波器阶数.通常将卷积混合模型转化成瞬时混合模型进行求解，具体转化方法如下[6]：
设置长度为w的时间窗，令mw≥n(w+l-1)，则窗内第i个观测EEG为xi(t)=[xi(t), xi(t-1), …, xi(t-w+1)]T，定义矢量x(t)=[x1T(t), x2T(t), …, xmT(t)]T，则公式(1)可重写为
|
(2) |
式中：s(t)=[s1T(t), s2T(t), …, snT(t)]T为加窗处理后的EEG源信号矢量，sj(t)=[sj(t), sj(t-1), …, sj(t-w-l+2)]T; A=(Aij)(i=1, 2, …, m, j=1, 2, …, n)是一个分块矩阵，其中每个分块均是w×(w+l-1)维Toeplitz矩阵：
|
(3) |
由此，式(1)卷积混合模型就转化为式(2)瞬时混合模型.
1.2 代价函数脑电源信号之间是互相独立的，且短时间内是平稳的，其自相关矩阵的非对角块部分等于零，符合JBD原理.所以，基于JBD确定代价函数，具体方法如下[7].
由于构成矢量s(t)的EEG源信号互相独立，且同一源信号在不同时延情况下是相关的，因此，s(t)的自相关矩阵Rs(τ)=E{s(t) sH(t-τ)}=bdiag{[Rs1(τ), Rs2(τ), …, Rsn(τ)]}为分块对角矩阵，其中Rsj(τ)=E{ sj(t)sjH(t-τ)}.对式(2)求自相关，可得
|
(4) |
设矩阵W为A的逆矩阵，即EEG解混矩阵，可得
|
(5) |
根据联合块对角化原理，取多个时延τq, q=1, …, Q，得到多个Rx(τq)，同时进行分块对角化，令WRx(τq)WH的非对角块部分尽量等于零，从而得到如下代价函数：
|
(6) |
式中：‖·‖F为矩阵的Frobenius范数; offbdiag(·)表示矩阵的非对角块部分.求得对W的估计
采用共轭梯度法对式(6)进行迭代寻优：
|
(7) |
|
(8) |
|
(9) |
式中：αk为步长因子，采用强Wolfe不精确线性搜索方法产生，gk=∇J(Wk)，计算公式为[8]
|
(10) |
具体迭代步骤如下：
①白化观测EEG信号.对EEG的零时延自相关矩阵Rx(0)进行特征值分解，得到Rx(0)=VDVH.式中，V为特征向量矩阵，D为特征值矩阵，D中所有特征值以降序排列.构造白化矩阵Λ=D
②给出解混矩阵的初始值W0并对其进行标准化，同时给出算法的终止阈值ε(ε > 0)，令k=0.
③依据式(10)计算g0=∇J(W0)，令d0=-g0.
④采用强Wolfe不精确线性搜索方法计算αk，依据式(7)计算Wk+1，并进行标准化Wk+1=Wk+1/‖Wk+1‖F.
⑤如果‖Wk+1-Wk‖F≤ε，算法停止，输出
⑥令ndim为W的维数，如果k=ndim，则令W0=Wk+1，转到步骤③.
⑦依据式(10)计算gk+1, 依据式(8)，式(9)计算dk+1，令k=k+1，转到步骤④.
1.4 收敛性分析由于在共轭梯度法寻优中，采用强Wolfe不精确线性搜索产生αk时，方向dk肯定是下降方向[9]，即J(Wk+1)-J(Wk)≤0.同时，由式(6)可以看出，代价函数的值是基于多个矩阵Frobenius范数的平方和计算的，因此该代价函数是非负的，存在一个下界0.则对于
给出3个时域信号s1，s2，s3，作为n=3的三维源信号.s1为Cz导联1s时间段静息EEG信号，采集设备使用Neuroscan 4.3系统，采样频率1 000 Hz，记录带宽0.1~100 Hz，皮肤阻抗≤5 kΩ.s2为频率50 Hz的正弦信号，用来模拟工频干扰，s3表示测量过程中混入的高斯噪声.
设观测EEG的维数m=4，传输路径FIR滤波器阶数l=3，滤波器系数随机生成，具体如下：
|
(1) |
其中，Hij=aij+bijz-1+cijz-2.
|
图 1 源信号 Fig.1 Source signals (a)—EEG信号; (b)—工频干扰; (c)—噪声. |
|
图 2 混合信号 Fig.2 Mixed signals (a)—混合信号1; (b)—混合信号2;(c)—混合信号3; (d)—混合信号4. |
取时间窗长度w=6，时延数Q=20，阈值ε=0.01，使用本文方法和文献[5]方法得到的对源信号的估计分别如图 3和图 4所示.
|
图 3 本文方法得到的分离信号 Fig.3 Signals obtained with the method proposed (a)—EEG信号估计; (b)—工频干扰估计; (c)—噪声估计. |
|
图 4 文献[5]方法得到的分离信号 Fig.4 Signals obtained with the method in document[5] (a)—工频干扰估计; (b)—EEG信号估计; (c)—噪声估计. |
为了验证方法的有效性，采用相关系数衡量实验结果，结果如表 1所示.可以看出，与文献[5]相比，本文方法得到的分离信号与源信号之间具有更大的相关程度，说明了本文方法的有效性.
|
表 1 两种方法得到的相关系数对比 Table 1 Comparison of correlation coefficients |
盲分离算法应用于生理信号的有效性与算法假设条件和生理信号之间的适应性有关.本文考虑源信号传输过程中的时延性和多向性，用卷积混合模型描述观测EEG，直接进行盲反卷积，得到较好的分离效果.而文献[5]则需要结合使用瞬时盲分离和盲反卷积，这可能是导致丢失部分有用信息、分离效果受到影响的原因.
4 结语本文提出一种基于盲反卷积的EEG盲分离方法，并通过仿真实验进行验证.结果表明，该方法可有效进行EEG盲分离，为针对EEG信号和其他生理信号的盲分离处理提供了借鉴和思路.
| [1] |
Sardouie S H, Shamsollahi M B, Albera L, et al.
Denoising of ictal EEG data using semi-blind source separation methods based on time-frequency priors[J]. IEEE Journal of Biomedical and Health Informatics , 2015, 19 (3) : 839–847.
DOI:10.1109/JBHI.2014.2336797 ( 0)
|
| [2] |
Selvan S, George S, Balakrishnan R.
Range-based ICA using a nonsmooth quasi-Newton optimizer for electroence phalographic source localization in focal epilepsy[J]. Neural Computation , 2015, 27 (3) : 628–671.
DOI:10.1162/NECO_a_00700 (0)
|
| [3] |
Hamaneh M B, Chitravas N, Kaiboriboon K, et al.
Automated removal of EKG artifact from EEG data using independent component analysis and continuous wavelet transformation[J]. IEEE Transactions on Biomedical Engineering , 2014, 61 (6) : 1634–1641.
DOI:10.1109/TBME.2013.2295173 (0)
|
| [4] |
Abutaleb A, Fawzy A, Sayed K.Blind deconvolution of EEG signals using the stochastic calculus[C]// Proceedings of 2012 Cairo International Biomedical Engineering Conference.Giza, 2012:175-178.
(0)
|
| [5] |
Dyrholm M, Makeig S, Hansen L K.
Model selection for convolutive ICA with an application to spatiotemporal analysis of EEG[J]. Neural Computation , 2007, 19 (4) : 934–955.
DOI:10.1162/neco.2007.19.4.934 (0)
|
| [6] |
Bousbia-Salah H, Belouchrani A, Abed-Meraim K.Blind separation of convolutive mixtures using joint block diagonalization[C]// Proceedings of the Sixth International Symposium on Signal Processing and Its Applications.Kuala Lumpur, 2001:13-16.
(0)
|
| [7] |
Ghennioui H, Fadaili E M, Thirion-Moreau N, et al.
A nonunitary joint block diagonalization algorithm for blind separation of convolutive mixtures of sources[J]. IEEE Signal Processing Letters , 2007, 14 (11) : 860–863.
DOI:10.1109/LSP.2007.903273 (0)
|
| [8] |
Ghennioui H, Thirion-Moreau N, Moreau E, et al.
Gradient-based joint block diagonalization algorithms:application to blind separation of FIR convolutive mixtures[J]. Signal Processing , 2010, 90 (6) : 1836–1849.
DOI:10.1016/j.sigpro.2009.12.002 (0)
|
| [9] |
孙文瑜, 徐成贤, 朱德通.
最优化方法[M]. 北京: 高等教育出版社, 2010 : 149 -158.
( Sun Wen-yu, Xu Cheng-xian, Zhu De-tong. The optimization method[M]. Beijing: Higher Education Press, 2010 : 149 -158. ) ( 0)
|