黄璐^1,2, 王宏³
1.东北大学 中荷生物医学与信息工程学院，辽宁 沈阳 110167;
2.大连海洋大学 信息工程学院，辽宁 大连 116023;
3.东北大学 机械工程与自动化学院，辽宁 沈阳 110819
收稿日期: 2015-05-18
基金项目: 国家自然科学基金资助项目(61071057).
作者简介: 黄璐(1979-), 女, 辽宁开原人, 东北大学博士研究生;
王宏(1960-), 女, 辽宁沈阳人, 东北大学教授, 博士生导师。

摘要: 采用卷积混合模型描述观测脑电信号(electroencephalogram，EEG)，提出一种基于盲反卷积的EEG盲分离方法.结合EEG源成分的独立性确定代价函数，并采用共轭梯度法进行迭代寻优.针对EEG仿真实验数据进行方法验证，采用分离信号与源信号之间的相关系数作为验证指标.实验结果表明，本文方法可以较好地实现EEG盲分离，为EEG信号处理和其他生理信号处理分析提供理论和方法借鉴.
关键词：脑电信号信号处理盲分离盲反卷积
Blind Separation of EEG Based on Blind Deconvolution
HUANG Lu^1,2, WANG Hong³
1.Sino-Dutch Biomedical and Information Engineering School, Northeastern University, Shenyang 110167, China;
2.College of Information Engineering, Dalian Ocean University, Dalian 116023, China;
3.School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Hong, E-mail:hongwang@mail.neu.edu.cn
Abstract: The convolution mixture model was adopted to describe the as-observed EEG signals, and then a method for the blind separation of EEG based on blind deconvolution was put forward. The cost function was established based on the independence of EEG sources, and iteration was carried out using conjugate gradient method. The verification was implemented with simulation experiment, adopting the correlation coefficients between separated signals and source signals as the verification indexes. Experimental results show that the method proposed can achieve blind separation of EEG successfully, providing a theoretical and practical reference for the processing of EEG and other physiological signals.
Key Words: physiological signalelectroencephalogramsignal processingblind separationblind deconvolution
在脑电(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信号x_i(i=1, 2, …, m)是由n维源信号s_j(j=1, 2, …, n)卷积混合而成，则混合模型可以表示如下：

(1)

式中：t为采样点; h_ij(·)为第j个源信号传输到第i个观测点的卷积混合过程，用FIR滤波器表示; l为FIR滤波器阶数.通常将卷积混合模型转化成瞬时混合模型进行求解，具体转化方法如下^[6]：
设置长度为w的时间窗，令mw≥n(w+l-1)，则窗内第i个观测EEG为x_i(t)=[x_i(t), x_i(t-1), …, x_i(t-w+1)]^T，定义矢量x(t)=[x₁^T(t), x₂^T(t), …, x_m^T(t)]^T，则公式(1)可重写为

(2)

式中：s(t)=[s₁^T(t), s₂^T(t), …, s_n^T(t)]^T为加窗处理后的EEG源信号矢量，s_j(t)=[s_j(t), s_j(t-1), …, s_j(t-w-l+2)]^T; A=(A_ij)(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)的自相关矩阵R_s(τ)=E{s(t) s^H(t-τ)}=bdiag{[R_s1(τ), R_s2(τ), …, R_sn(τ)]}为分块对角矩阵，其中R_sj(τ)=E{ s_j(t)s_j^H(t-τ)}.对式(2)求自相关，可得

(4)

设矩阵W为A的逆矩阵，即EEG解混矩阵，可得

(5)

根据联合块对角化原理，取多个时延τ_q, q=1, …, Q，得到多个R_x(τ_q)，同时进行分块对角化，令WR_x(τ_q)W^H的非对角块部分尽量等于零，从而得到如下代价函数：

(6)

式中：‖·‖_F为矩阵的Frobenius范数; offbdiag(·)表示矩阵的非对角块部分.求得对W的估计后，通过计算(t)=x(t)即可求得对EEG源信号的估计.
1.3 迭代算法采用共轭梯度法对式(6)进行迭代寻优：

(7)

(8)

(9)

式中：α_k为步长因子，采用强Wolfe不精确线性搜索方法产生，g_k=?J(W_k)，计算公式为^[8]

(10)

具体迭代步骤如下：
①白化观测EEG信号.对EEG的零时延自相关矩阵R_x(0)进行特征值分解，得到R_x(0)=VDV^H.式中，V为特征向量矩阵，D为特征值矩阵，D中所有特征值以降序排列.构造白化矩阵Λ=DV^H，以此矩阵对观测EEG信号x(t)进行白化.
②给出解混矩阵的初始值W₀并对其进行标准化，同时给出算法的终止阈值ε(ε > 0)，令k=0.
③依据式(10)计算g₀=?J(W₀)，令d₀=-g₀.
④采用强Wolfe不精确线性搜索方法计算α_k，依据式(7)计算W_k+1，并进行标准化W_k+1=W_k+1/‖W_k+1‖_F.
⑤如果‖W_k+1-W_k‖_F≤ε，算法停止，输出=W_k+1.
⑥令n_dim为W的维数，如果k=n_dim，则令W₀=W_k+1，转到步骤③.
⑦依据式(10)计算g_k+1, 依据式(8)，式(9)计算d_k+1，令k=k+1，转到步骤④.
1.4 收敛性分析由于在共轭梯度法寻优中，采用强Wolfe不精确线性搜索产生α_k时，方向d_k肯定是下降方向^[9]，即J(W_k+1)-J(W_k)≤0.同时，由式(6)可以看出，代价函数的值是基于多个矩阵Frobenius范数的平方和计算的，因此该代价函数是非负的，存在一个下界0.则对于ε > 0，总存在K，使得k > K时，|J(W_k+1)－J(W_k)| < ε，因此，算法收敛.
2 仿真实验数据给出3个时域信号s₁，s₂，s₃，作为n=3的三维源信号.s₁为Cz导联1s时间段静息EEG信号，采集设备使用Neuroscan 4.3系统，采样频率1?000 Hz，记录带宽0.1~100 Hz，皮肤阻抗≤5 kΩ.s₂为频率50 Hz的正弦信号，用来模拟工频干扰，s₃表示测量过程中混入的高斯噪声.
设观测EEG的维数m=4，传输路径FIR滤波器阶数l=3，滤波器系数随机生成，具体如下：

(1)

其中，H_ij=a_ij+b_ijz^-1+c_ijz^-2.
源信号及混合信号分别如图 1和图 2所示.
图 1(Fig. 1)

图 1 源信号Fig.1 Source signals (a)—EEG信号; (b)—工频干扰; (c)—噪声.

图 2(Fig. 2)

图 2 混合信号Fig.2 Mixed signals (a)—混合信号1; (b)—混合信号2;(c)—混合信号3; (d)—混合信号4.

3 实验结果分析取时间窗长度w=6，时延数Q=20，阈值ε=0.01，使用本文方法和文献[5]方法得到的对源信号的估计分别如图 3和图 4所示.
图 3(Fig. 3)

图 3 本文方法得到的分离信号Fig.3 Signals obtained with the method proposed (a)—EEG信号估计; (b)—工频干扰估计; (c)—噪声估计.

图 4(Fig. 4)

图 4 文献[5]方法得到的分离信号Fig.4 Signals obtained with the method in document[5] (a)—工频干扰估计; (b)—EEG信号估计; (c)—噪声估计.

为了验证方法的有效性，采用相关系数衡量实验结果，结果如表 1所示.可以看出，与文献[5]相比，本文方法得到的分离信号与源信号之间具有更大的相关程度，说明了本文方法的有效性.
表 1(Table 1)

表 1 两种方法得到的相关系数对比Table 1 Comparison of correlation coefficients 信号类别 本文方法 文献[5]方法 EEG信号 0.971 3 0.939 5 工频干扰 0.986 0 0.952 4 噪声 0.945 1 0.920 8 平均 0.967 5 0.937 6

表 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)