CN102697493B - 一种快速的脑电信号中眼电伪迹自动识别和去除的方法 - Google Patents

Abstract

本发明提出了一种快速的脑电信号中眼电伪迹自动识别和去除的方法，属于生物信息技术领域，主要应用于脑电信号采集与预处理的过程中。具体包括：将采集得到的多导脑电信号和眼电信号进行离散小波变换，获取多尺度的小波系数；将串接小波系数作为独立分量分析的输入，利用基于负熵判据的FastICA算法实现独立成分的快速获取；通过夹角余弦法识别出眼电伪迹后，将该独立成分置零，并经过ICA逆变换将其余成分投影返回到原信号各个电极；最后通过反演小波变换得到去除眼电伪迹的脑电信号。本发明解决了ICA方法应用于含噪脑电信号中分离效果差，收敛速度慢的问题，实现了从脑电中快速自动去除眼电伪迹的功能。

Description

一种快速的脑电信号中眼电伪迹自动识别和去除的方法

技术领域：

本发明属于脑电信号(Electroencephalography，EEG)采集与预处理技术领域。具体涉及一种基于离散小波变换(Discrete Wavelet Transform,DWT)与独立分量分析(IndependentComponentAnalysis，ICA)的快速的脑电信号中眼电伪迹自动去除方法。

背景技术

脑电信号是一种反映大脑活动的生物电信号，在研究人脑功能和疾病的过程中发挥着越来越大的作用。然而脑电信号非常微弱，具有很高的时变敏感性，在采集时极易受到无关噪声的干扰，从而形成了各种EEG伪迹。这些伪迹给脑电信号的分析处理带来了很大的困难。眼电(Electro-oculogram，EOG)是EEG信号中的一种最主要的干扰成分，它会随机地出现在脑电信号中，且幅度一般较大，导致采集的EEG信号产生明显畸变，形成眼电伪迹(Ocular Artifact,OA)。在临床脑电检查中，医生通常舍弃含有EOG干扰的脑电数据段，寻找较干净的脑电信号进行观测分析。但是舍弃数据段可能导致重要信息丢失。因此，眼动干扰的消除一直是脑电信号预处理中重要的研究内容。

近年来，研究者们提出了多种脑电信号中眼电伪迹去除的方法，常用的有回归方法、主成分分析PCA、小波变换和独立分量分析ICA等。

其中，回归方法主要用于眨眼或眼球的运动引起的眼电伪迹，依赖于建立一个正确的回归导联，但由于眼电信号和脑电信号的激活扩散都具有双向性(bi-directionality)，因此回归方法去除伪迹时会错误地去除某些脑电信号。

Jung和Berger等提出了用PCA方法进行眼电伪迹的去除，在被试验者完成眼动和眨眼任务时记录EEG和EOG信号，再计算出这些信号的主成分，作为眼动和眨眼伪迹的主成分。然后从混合信号去除该成分，得到校正后的信号。研究表明PCA在效果上显著优于回归方法，然而，PCA不能完全从EEG中分离与它的波形相似的电位的噪声。

小波变换是傅立叶变换的新发展，小波变换系数能反映信号在时域及频域的局部信息。因而，小波分析作为一种时频分析方法，在生物医学信号处理方面有着广阔的应用前景，特别适合像EEG这类非平稳信号的处理。小波门限法去噪是基于小波变换多分辨率分析的一种方法。由于信号和噪声经过小波变换后的统计特性不同，从而在多尺度分析中呈现出不同的传播行为。选取一个合适的阈值，并用此阈值对小波分解后得到的各层小波空间里的细节系数进行截断，而保持尺度空间里的逼近系数不变，然后再进行逆变换，即可得到去噪后的信号。小波门限法去噪要求信号和噪声的频带不能混叠，但是EEG和EOG伪迹的频带相混叠，因此去噪效果不太理想。

独立分量分析是由盲源信号分离技术发展来的多道信号处理方法，能够取得比较理想的去噪效果，被广泛应用于去除EEG信号中的眼电伪迹。ICA的思想来自于中心极限定理：一组均值和方差为同一数量级的随机变量，共同作用的结果必接近高斯分布。因此若要分离一组由相互统计独立的信源经线性组合而产生的混和信号时，只要对分离结果的非高斯性进行度量，当其非高斯性达到最大时，可以认为混合信号实现了最佳分离。ICA的模型可以用下式表示：

c=A·s    (1)

r=W·c    (2)

其基本思想可描述如下：设c(t)=[c₁(t),c₂(t),…,c_n(t)]^T是n维的观测信号，s(t)=[s₁(t),s₂(t),…,s_n(t)]^T是产生观测信号的n个相互统计独立的源信号，且观测信号c(t)是源信号s(t)经过一个未知矩阵A∈R^n×n线性混和而产生的，如公式(1)所示。在混和矩阵A和源信号s(t)未知的情况下，仅利用观测信号c(t)和源信号统计独立的假设，寻找一个线性变换分离矩阵W∈R^n×n，希望输出信号r(t)=Wc(t)=WAs(t)尽可能逼近真实的源信号s(t)，如公式(2)所示。

理论上认为脑电信号和其他干扰源所产生的干扰信号都是由相互独立的信源产生的，从而适用ICA方法来除去伪迹，提取出有用的脑电信号。ICA方法是PCA的一种延伸，是在所有统计意义下的去相关，因此比PCA更具有优越性。但是ICA去除伪迹也存在一些需要探讨和解决的问题，这些问题极大限制了其在脑电信号实际在线预处理中的广泛应用。

首先，由于经典的ICA模型没有考虑噪声的存在，如果用ICA算法对含高斯噪声的观测信号直接进行独立分量分析，有时会产生较大的误差。因此该方法要求观测信号中最多只有一个高斯型信源，而脑电信号的产生机理复杂，且在实际采集过程中除了眼电的影响，还受到各种各样其余噪声的干扰。这些噪声不仅使ICA算法的分离效果严重变坏，而且需要通过迭代多次才能得到ICA算法的分离矩阵，计算量大且速度很慢。Jafari于2003年提出基于小波变换的ICA方法来提取胎儿心电(FECG)，近年来有研究将该方法应用于图像处理，事件相关电位的提取等。这种方法可以较好得克服ICA模型中混有噪声的问题，具有很好的抗噪能力。但目前这种方法还没有应用于去除脑电信号中的眼电伪迹。

其次，因为ICA分离信号顺序的不确定性，判断独立成分是否为伪迹是个困难的问题，传统采用的方法是结合视觉观测的半自动方法，但这种方法难以应用于实时在线地处理脑电信号。

发明内容

针对独立分量分析ICA在去除脑电中的眼电伪迹时，因受到多种噪声的干扰而影响算法的分离效果和速度，以及无法自动识别眼电伪迹这两个不足，本发明提出了一种基于离散小波变换与独立分量分析的脑电信号中眼电伪迹的自动去除方法，即DWT-ICA，在去除眼迹的过程中完全无须人工干预，该方法的抗噪能力很强，收敛速度快且实时性好，同时极大地提高了脑电信号的信噪比，为脑电信号的在线预处理提供了新的思路。

本发明为解决所述技术问题，采用如下的技术方案：

首先将实际采集得到的多导脑电信号和一导眼电参考信号进行离散小波变换，获取多尺度的小波变换系数，将各导小波变换系数作为独立分量分析的输入，利用基于负熵判据的快速ICA（FastICA）算法实现独立成分的快速获取，通过夹角余弦法准则识别出眼电伪迹后，将该独立成分置零，并经过ICA逆变换将其余成分进行投影映射，最后通过反演离散小波变换进行信号重构即可得到去除眼电伪迹的脑电信号。具体步骤如下：

(1)采集n+1导信号x(t)；

首先按照国际标准10-20系统来安放电极，通过电极帽采集n导脑电信号，同时在眼睛附近放置一个电极，采集1导眼电参考信号，共采集了n+1导信号x(t)=[x₁(t),x₂(t)，…,x_n+1(t)]^T∈R^(n+1)×N。

其中，x_i(t)为采集到的任意导信号，i代表每导信号对应的序号，i为正整数且i∈[1，n+1]，t为每导信号的样本点，t为正整数且t∈[1，N]，N为样本点总数，R为实数集，设第l导为眼电信号；

眼电信号EOG是EEG信号中的一种最主要的干扰成分，它会随机地出现在脑电信号中，且幅度一般较大，导致采集到的EEG信号产生明显畸变，形成眼电伪迹。眼电伪迹给脑电信号的分析处理带来了很大的困难。因此，眼电伪迹的消除是脑电信号预处理中必不可少的研究内容。该步骤中在采集脑电信号的同时，采集了一导眼电参考信号，以便于在后续步骤中利用有效方法自动快速地识别眼电伪迹。

(2)计算采集到的n+1导信号x(t)的离散小波变换系数矢量

具体包括以下步骤；

2.1计算采集到的任意导信号x_i(t)的第j层分解的逼近系数分量

与细节系数分量

选择小波基函数，并采用Mallat算法对采集到的任意导信号x_i(t)进行L层分解，得到采集到的任意导信号x_i(t)的第j层分解的逼近系数分量

与细节系数分量

其中，j代表分解尺度，j∈[1，L]且为整数；

2.2计算采集到的任意导信号x_i(t)的离散小波变换系数矢量

将采集到的任意导信号x_i(t)的逼近系数分量及细节系数分量进行串接，得到x_i(t)的离散小波变换系数矢量

其结构如下：

${\tilde{x}}_i\left(t\right)={[a_i^L,d_i^L,d_i^{L-1},\·\·\·,d_i^1]}^T\∈R^{M\×1}---\left(3\right)$

其中，R为实数集，M为每导信号小波系数的样本点总数；

2.3计算n+1导信号x(t)的离散小波变换系数矢量

即

$\tilde{x}\left(t\right)={[{\tilde{x}}_1\left(t\right),{\tilde{x}}_2\left(t\right),\·\·\·,{\tilde{x}}_{n+1}\left(t\right)]}^T\∈R^{(n+1)\×M}---\left(4\right)$

小波变换是在傅立叶变换的基础上发展起来的一种时频分析方法，在生物医学信号处理方面有着广阔的应用前景，特别适合像EEG这类非平稳信号的处理。由于离散小波变换的计算速度快，适合于在线分析，因此在本发明中采用了离散小波变换方法。实践证明，就概率密度函数而言，离散小波变换系数比原始信号的超高斯性更强，因而在小波域对离散小波变换系数进行独立分量分析ICA，将使迭代的收敛速度更快，抗噪能力更强。这正是在独立分量分析算法前对信号进行离散小波变换的目的。

为了简化小波变换算法并加快运算的速度，Mallat根据图像分解与重构塔式算法的思想，提出了小波分析与重构的快速算法。该算法在多分辨率分析的基础上，给出了如何通过双通道滤波器组实现信号的离散小波变换及反变换，为离散小波变换在工程中的应用提供了非常便捷的手段。

典型的三层多分辨率小波分解树如图3所示。图中，A_j为j尺度下的逼近系数分量，D_j为j尺度下的细节系数分量，分解满足表达式：

f(t)=A₁+D₁+D₂+D₃。由此可以推断，L层小波分解满足表达式：

$f\left(t\right)=A_L+\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{D}_j---\left(5\right)$

式中，L代表最大分解尺度。

(3)把步骤(2)中得到的离散小波变换系数矢量

作为独立分量分析算法所需要的多导输入，独立分量分析算法采用基于负熵判据的FastICA算法，FastICA算法的逼近算法选用紧缩逼近法，依次提取单个权值向量w_i(μ)(i为正整数且i∈[1，n+1])，并在提取每个权值向量w_i(μ)前利用格拉姆-施密特正交化分解剔除所有已提取过的权值向量，通过不断迭代得到n+1个权值向量，从而计算分离矩阵W∈R^(n+1)×(n+1)；

独立分量分析ICA方法的关键问题是建立一个能够度量分离结果独立性的目标函数及其相应的分离算法。本发明中选择负熵最大作为目标函数，并采用快速ICA（FastICA）算法作为相应的分离算法。与其它分离算法相比，FastICA性能较好，具有很多优点，主要有以下几个方面：①该算法的收敛速度至少是二次方的，意味着FastICA具有非常快的收敛速度。②该算法和基于梯度的算法不同，其不需要选择步长参数，使用方便。③FastICA算法继承了神经网络算法并行、分布的优点，计算简单，且所需的内存空间较小。鉴于此，本发明采用的是FastICA算法。

FastICA算法有两种逼近方法：紧缩逼近法和均衡逼近法。本发明中采用的是紧缩逼近法。该逼近方法是分别对分离矩阵W的每一列进行更新，通过系统学习找出一个方向，即W中的一个权值向量w_i(μ)，使得信号在该方向上的投影具有最大的非高斯性。因此紧缩逼近法每次只迭代得到一个权值向量w_i(μ)，同时为了保证每次提取出来的都是尚未提取过的权值向量，本发明采用格拉姆-施密特正交化分解将已提取过的权值向量去掉，以解决权值向量相同的问题。

具体提取单个权值向量w_i(μ)的步骤如下：

3.1对输入信号

进行预处理，得到新的白噪声化信号z(t)，计算公式如下：

$z\left(t\right)=V\hat{x}\left(t\right)---\left(6\right)$

其中：V∈R^(n+1)×(n+1)是

的白化矩阵，满足E(VV^T)＝I， $\hat{x}\left(t\right)=\tilde{x}\left(t\right)-E\left\{\tilde{x}\left(t\right)\right\}$

3.2计算w_i(μ+1)＝E{z(t)g[w_i ^T(μ)z(t)]}-E{g′[w_i ^T(μ)z(t)]，}w_i(μ)

其中，μ为迭代次数，g为一个任意非二次型函数的导数，取任意初始向量w_i(0)∈R^(n+1)×1，要求它满足||w_i(0)||＝1，i为正整数且i∈[1，n+1]；

3.3归一化： $w_i(\mathrm{\μ}+1)=\frac{w_i(\mathrm{\μ}+1)}{\left|\right|w_i(\mathrm{\μ}+1)\left|\right|};$

3.4判断w_i(μ+1)是否收敛，如果没有收敛返回步骤3.3，否则迭代结束，输出最终w_i(μ+1)作为得到的权值向量w_i(μ)；

将每次得到的w_i(μ)组成分离矩阵W，即

W＝[w₁(μ)，w₂(μ)，…，w_n+1(μ)]^T∈R^(n+1)×(n+1)；    (7)

(4)对离散小波变换系数矢量进行变换，计算独立成分矩阵y(t)，计算公式如下：

$y\left(t\right)=W\tilde{x}\left(t\right)={[y_1\left(t\right),y_2\left(t\right),\·\·\·,y_{n+1}\left(t\right)]}^T---\left(8\right)$

其中，y_i(t)为第i个独立成分，y_i(t)是由M个元素组成的列向量；

由于步骤(3)与步骤(4)在小波域对离散小波变换系数进行独立分量分析，因此迭代求取独立成分的收敛速度很快，也具有非常好的分离效果，克服了其余噪声的干扰。

(5)分别计算每个独立成分y_i(t)与第l导眼电参考信号x_l(t)的离散小波变换系数矢量

的夹角余弦值cosθ_i，其取值范围为[-1，1]，计算公式如下：

$\cos {\mathrm{\θ}}_i=\frac{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{y}_{\mathrm{iq}}{\tilde{x}}_{\mathrm{lq}}}{\sqrt{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{y_{\mathrm{iq}}}^2}\sqrt{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{{\tilde{x}}_{\mathrm{lq}}}^2}}---\left(9\right)$

其中，M为每导信号小波系数的样本点总数，y_iq是列向量y_i(t)中的元素，

是由M个元素组成的列向量，是列向量

的元素，q代表元素的序号；

因为独立分量分析ICA是以源信号彼此相互独立为前提，实际情况中得到的是观测信号，无法得到真正的源信号，独立分量分析ICA通过构造对分离结果独立性度量的目标函数，采用分离算法实现独立成分的提取。但提取的各独立成分幅值的大小、正负性及顺序均具有不确定性。因此如何自动判别独立成分是否为伪迹是个困难的问题。

几何中采用夹角余弦来衡量两个模式向量的相似度，机器学习中常借用这一概念来衡量样本向量之间的差异，本发明中提出用该方法来衡量在小波域进行ICA分解得到的各导独立成分与眼电参考信号的离散小波系数的差异，从而自动识别眼电伪迹。

(6)剔除独立成分矩阵y(t)中的眼电伪迹成分，得到无眼电伪迹成分的独立成分矩阵

对步骤(5)中得到的各个独立成分y_i(t)与第l导眼电信号x_l(t)的离散小波变换系数矢量

的夹角余弦值cosθ_i的绝对值进行排序，最大的cosθ_i|值所对应的独立成分y_i(t)即为眼电伪迹成分，将独立成分矩阵y(t)中的眼电伪迹成分y_i(t)置零，其余独立成分不变，此时信号即为

(7)经过独立分量分析ICA逆变换将进行投影映射，得到小波变换系数u(t)，即:

$u\left(t\right)=W^{-1}\tilde{y}\left(t\right)={[u_1\left(t\right),u_2\left(t\right),\·\·\·,u_{n+1}\left(t\right)]}^T---\left(10\right)$

ICA算法提取得到的各独立成分幅值的大小、正负性及顺序均具有不确定性，而该步骤经过独立分量分析ICA逆变换可将独立成分进行投影映射，使其与源信号幅值的大小、正负以及顺序均具有一致性。

(8)选择与步骤2.1中相同的小波基函数，采用Mallat塔式重构算法对步骤7中得到的各导小波变换系数u_i(t),i∈[1，n+1]且i≠l,重构脑电信号,即得到去除眼电伪迹后的n导脑电信号。

本发明的有益效果是：

(1)本发明提出了一种基于离散小波变换与独立分量分析的脑电信号中眼电伪迹的去除方法，即DWT-ICA，该方法具有较强的抗噪能力，有效缩短了ICA的迭代过程，因而收敛速度更快，实时性更好，同时极大地提高了脑电信号的信噪比。

(2)本发明针对ICA分离信号顺序的不确定性，提出了使用夹角余弦法来识别眼电伪迹，将度量信号间的相似度转换为计算向量间的相似度，该方法具有很好的分辨能力，能够实现眼电伪迹的自动识别，完全无需人工的干预，且精确度很高，这为脑电信号的在线预处理提供了新的思路。

附图说明：

图1基于离散小波变换与ICA方法去除眼电伪迹的流程框图；

图2国际标准10-20电极系统的放置示意图；

图3三层小波分解树结构图；

图4采集的脑电及眼电信号，其中，(a)为C3导脑电信号，(b)为C4导脑电信号，(c)为眼电信号；

图5去除眼电伪迹后的脑电信号，其中，(a)为C3导去除眼迹的脑电信号，

(b)为C4导去除眼迹的脑电信号；

图6国际标准数据库中纯净的脑电信号，其中，(a)为C3导脑电信号，(b)为Cz导脑电信号，(c)为C4导脑电信号；

图7混入眼电伪迹及高斯白噪声的脑电信号，其中，(a)为C3导含噪脑电信号，(b)为Cz导含噪脑电信号，(c)为C4导含噪脑电信号；

图8ICA方法去除眼电伪迹后的脑电信号，其中，(a)为C3导去除眼迹的脑电信号，(b)为Cz导去除眼迹的脑电信号，(c)为C4导去除眼迹的脑电信号；

图9DWT-ICA方法去除眼电伪迹后的脑电信号，其中，(a)为C3导去除眼迹的脑电信号，(b)为Cz导去除眼迹的脑电信号，(c)为C4导去除眼迹的脑电信号；

图10DWT-ICA方法去除眼电伪迹前后信号的功率谱对比，其中，(a)为C3导脑电信号的功率谱对比，(b)为Cz导脑电信号的功率谱对比，(c)为C4导脑电信号的功率谱对比；

图4至图10均使用matlab仿真环境下的plot函数实现。

具体实施方式

以下结合附图对本发明的具体实施做进一步说明。本实施例是在matlab的仿真环境下进行的。

1.参见图2，按照国际标准导联10-20系统，在C3,C4处安放电极，采集了2导脑电信号x₁(t)与x₂(t)，另外还在F7处同步采集了一导眼电信号x₃(t)，共采集了3导信号x(t)=[x₁(t),x₂(t),x₃(t)]^T，参考电极放置在A1和A2。信号的采样频率为250Hz，模拟滤波范围为0.1~100Hz。每次EEG记录时间为10s，采样点数为2500。

参见图4，为采集的2导脑电信号及1导眼电信号波形，可见眼电信号对各导脑电均产生了不同程度的影响，形成了眼电伪迹。

2.利用matlab环境下离散小波变换工具箱中的wavedec函数对每导信号x₁(t)，x₂(t)，x₂(t)进行离散小波变换，选择小波基‘sym4’，此处进行3层分解，得到信号x₁(t)的逼近系数分量

与细节系数分量

得到信号x₂(t)的逼近系数分量

与细节系数分量

得到信号x₃(t)的逼近系数分量

与细节系数分量

将每一导信号经过3层分解后得到的逼近系数分量及细节系数分量进行串接，构成一个一维矢量，其结构即为 ${\tilde{x}}_1\left(t\right)={[a_1^3,d_1^3,d_1^2,d_1^1]}^T\∈R^{M\×1},$ ${\tilde{x}}_2\left(t\right)={[a_2^3,d_2^3,d_2^2,d_2^1]}^T\∈R^{M\×1},$ ${\tilde{x}}_3\left(t\right)={[a_3^3,d_3^3,d_3^2,d_3^1]}^T\∈R^{M\×1},$ 其中，R为实数集，M为每导信号小波系数的样本点总数，由选择的小波基和分解层数来决定，此处为2501。

离散小波变换系数矢量 $\tilde{x}\left(t\right)={[{\tilde{x}}_1\left(t\right),{\tilde{x}}_2\left(t\right),{\tilde{x}}_3\left(t\right)]}^T.$

3.把作为ICA方法的输入，利用matlab环境下FastICA工具箱中的fastica函数，通过不断迭代得到分离矩阵W，此处设定迭代的精度为O.0001。当权值向量达到此精度时，就结束迭代过程，此时的权值即为最终的权值向量。

经过迭代，此处求出的分离矩阵为：

$W=\left[\begin{array}{ccc}0.0088& 0.0063& -0.0205\\ 0.0731& -001743& 0.0085\\ -0.2373& 0.0731& 0.0196\end{array}\right]$

再由

求取独立成分y(t)＝[y₁(t)，y₂(t)，y₃(t)]^T。

4.分别计算每导独立成分y₁(t)，y₂(t)，y₃(t)与第3导眼电参考信号的小波变换系数

的夹角余弦值cosθ₁，cosθ₂，cosθ₃。将这3个夹角余弦值的绝对值进行排序，最大的|cosθ_i|(i∈[1，3])值所对应的独立成分即为眼电伪迹成分，将该成分y_i(t)置零，其余独立成分不变，即得到无眼电伪迹成分的独立成分矩阵

5.经过ICA逆变换将

投影进行投影映射，得到小波变换系数u(t)，即：

$u\left(t\right)=W^{-1}\tilde{y}\left(t\right)={[u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)]}^T---\left(11\right)$

这时已在小波域中去除了眼电伪迹。

6.选择小波基函数‘sym4’，采用离散小波变换工具箱中的waverec函数对上一步骤中去除了眼电伪迹的小波变换系数u₁(t)与u₂(t)重构脑电信号，如图5所示，可以看出此时脑电信号中的眼电伪迹已经被很好的去除了，信噪比得到了极大的提高。

为了进一步验证本发明去除眼电伪迹的有益效果，下面将采用国际标准数据库中纯净的脑电信号来构造混有眼电伪迹和其它噪声的采集数据，并从去噪效果以及时间方面来定量评估该方法的效果。

采集数据的构造：

脑电数据来源于“BCI Competition 2003”竞赛数据库，该实验采用Agcl电极，数据由差分电极从国际标准的10-20导联系统的C3、Cz、C4（图2）三个通道获得，脑电信号的采样频率是128Hz，每次采集实验持续9s，共采集了140组脑电数据，本发明选取了其中一组纯净的脑电数据，如图6所示。将采集的眼电信号按一定比例混合到C3、Cz、C4导脑电中，并加入5分贝瓦(dBW)的高斯白噪声来模拟其他噪声的影响，以验证DWT-ICA方法的抗噪能力。各导信号的波形如图7所示，其信噪比(signal-to-noise ratio,SNR)分别为-3.96dB，-2.73dB，-3.80dB。SNR的计算公式如下：

$\mathrm{SNR}=10\lg \{\underset{k=1}{\overset{C}{\mathrm{\Σ}}}{h}^2\left(k\right)/\underset{k=1}{\overset{C}{\mathrm{\Σ}}}{[h\left(t\right)-\hat{h}\left(t\right)]}^2\}---\left(12\right)$

其中，h(k)为某个电极上不含眼电伪迹的纯净脑电信号，

为目前重建的脑

电信号，k为样本点，C为样本点总数。信噪比的值越大意味着重建的EEG信号越接近真实的EEG信号。

分别采用ICA方法来去除眼电伪迹，去除眼电后的脑电信号波形如图8所示，以及采用本发明提出的DWT-ICA方法来去除眼电伪迹，去除眼电后的脑电信号波形如图9所示。其中，本发明提出的DWT-ICA方法中的离散小波变换选用小波基‘sym8’，进行3层分解，ICA算法和DWT-ICA方法迭代过程中的迭代精度均设为0.0001。两种方法去除脑电信号中眼电伪迹后的信噪比以及所用的时间如表1所示：

表1ICA与DWT-ICA去除眼电的信噪比及所耗时间

从表中可以看出DWT-ICA方法与ICA相比，去噪效果较好，同时DWT-ICA在时间上大大减少了，可见本发明提出的眼电伪迹自动去除方法具有很强的抗噪能力，收敛速度快且实时性好，同时极大地提高了脑电信号的信噪比，为脑电信号的在线预处理提供了新的思路。

为进一步从能量角度观察去噪的效果，分别对原纯净的脑电信号，含有眼电伪迹的EEG以及去除眼电伪迹的EEG进行AR模型功率谱估计，如图10所示，眼电伪迹的影响一般在低频段，经过DWT-ICA方法的去噪后，眼电伪迹几乎被全部剔除了，与原纯净脑电的功率谱吻合较好，即脑电的能量得到了很好的恢复，从而再次说明本研究提出的DWT-ICA方法应用与EOG的去除的有效性及正确性。

Claims (1)

1.一种快速的脑电信号中眼电伪迹自动识别和去除的方法，其特征在于，包括以下步骤：

(1)采集n+1导信号x(t)；

首先按照国际标准10-20系统来安放电极，通过电极帽采集n导脑电信号，同时采集1导眼电信号，共采集了n+1导信号x(t)=[x₁(t),x₂(t),…,x_n+1(t)]^T∈R^(n+1)×N，

其中，x_i(t)为采集到的任意导信号，i代表每导信号对应的序号，i为正整数且i∈[1,n+1]，t为每导信号的样本点且t为正整数，N为样本点总数，R为实数集，设第l导为眼电信号；

(2)计算采集到的n+1导信号x(t)的离散小波变换系数矢量

具体包括以下步骤；

2.1选择小波基函数，并采用Mallat算法对采集到的任意导信号x_i(t)进行L层分解，得到采集到的任意导信号x_i(t)的第j层分解的逼近系数分量

与细节系数分量

其中，j代表分解尺度，j∈[1，L]且为整数；

2.2将采集到的任意导信号x_i(t)的逼近系数分量及细节系数分量进行串接，得到任意导信号x_i(t)的离散小波变换系数矢量

其结构如下：

${\tilde{x}}_i\left(t\right)=[a_i^L,d_i^L,d_i^{L-1},...,d_i^1{]}^T\∈R^{M\×1}---\left(1\right)$

其中，R为实数集，M为每导信号小波系数的样本点总数；

2.3计算n+1导信号x(t)的离散小波变换系数矢量

即

$\tilde{x}\left(t\right)=[{\tilde{x}}_1\left(t\right),{\tilde{x}}_2\left(t\right),...,{\tilde{x}}_{n+1}\left(t\right){]}^T\∈R^{(n+1)\×M}---\left(2\right)$

其中

为步骤2.2得到的离散小波变换系数矢量；

(3)计算分离矩阵W；

把步骤(2)中得到的离散小波变换系数矢量

作为独立分量分析算法所需要的多导输入，独立分量分析算法采用基于负熵判据的FastICA算法，FastICA算法的逼近算法选用紧缩逼近法，依次提取单个权值向量w_i(μ),i为正整数且i∈[1,n+1]，并在提取每个权值向量w_i(μ)前利用格拉姆-施密特正交化分解剔除所有已提取过的权值向量，通过不断迭代得到n+1个权值向量，从而计算分离矩阵W，计算公式如下：

W=[w₁(μ),w₂(μ),…,w_n+1(μ)]^T∈R^(n+1)×(n+1)    (3)

（4）对离散小波变换系数矢量进行变换，计算独立成分矩阵y(t)，计算公式如下：

$y\left(t\right)=W\tilde{x}\left(t\right)=[y_1\left(t\right),y_2\left(t\right),...,y_{n+1}\left(t\right){]}^T---\left(4\right)$

其中，y_i(t)表示第i个独立成分，y_i(t)是由M个元素组成的列向量；

（5）分别计算每个独立成分y_i(t)与第l导眼电信号x_l(t)的离散小波变换系数矢量

的夹角余弦值cosθ_i，其取值范围为[-1,1]，计算公式如下:

$\cos {\mathrm{\θ}}_i=\frac{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{y}_{\mathrm{iq}}{\tilde{x}}_{\mathrm{lq}}}{\sqrt{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{y_{\mathrm{iq}}}^2}\sqrt{\underset{q=1}{\overset{M}{\mathrm{\Σ}}}{{\tilde{x}}_{\mathrm{lq}}}^2}}---\left(5\right)$

其中，M为每导信号小波系数的样本点总数，y_iq是列向量y_i(t)中的元素，

是由M个元素组成的列向量，

是列向量

的元素，q代表元素的序号；

（6）剔除独立成分矩阵y(t)中的眼电伪迹成分，得到无眼电伪迹成分的独立成分矩阵

对步骤（5）中得到的各个独立成分y_i(t)与第l导眼电信号x_l(t)的离散小波变换系数矢量

的夹角余弦值cosθ_i的绝对值进行排序，最大的|cosθ_i|值所对应的独立成分y_i(t)即为眼电伪迹成分，将独立成分矩阵y(t)中的眼电伪迹成分y_i(t)置零，其余独立成分不变，即得到无眼电伪迹成分的独立成分矩阵

（7）经过独立分量分析ICA逆变换将

进行投影映射，得到小波变换系数u(t)，即:

$u\left(t\right)=W^{-1}\tilde{y}\left(t\right)=[u_1\left(t\right),u_2\left(t\right),...,u_{n+1}\left(t\right){]}^T---\left(6\right)$

其中，u_i(t)表示各导小波变换系数，i∈[1,n+1]；

（8）选择与步骤2.1中相同的小波基函数，采用Mallat塔式重构算法对步骤7中得到的各导小波变换系数u_i(t)重构脑电信号,其中i∈[1,n+1]且i≠l，即得到去除眼电伪迹后的n导脑电信号。

Images