CN104688220A - 一种去除脑电信号中眼电伪迹的方法 - Google Patents

Abstract

一种去除脑电信号中眼电伪迹的方法，结合ICA和MEMD这两种方法，MEMD从脑电信号中完全提取出眼电伪迹，并且保留了部分脑电信息，这样去除眼电的时候，避免了部分脑电信息的损失；ICA从脑电信号中分离出来眼电伪迹，通过计算独立分量的四阶累积量判定伪迹，然后去除眼电；本发明的优势用于眼电伪迹的去除，不仅能够自动识别眼电伪迹并将其去除，同时保留了大量脑电信息，为脑电去噪提供了一种新思路；该方法不仅能够有效地去除眼电伪迹，而且保留了大量有用的脑电信息，通过与单独的ICA方法进行对比，结果表明IMEMD方法去噪效果更好，能够进一步提高信号信噪比和减少信号均方误差。

Description

一种去除脑电信号中眼电伪迹的方法

技术领域

本发明涉及脑电信号预处理方法，特别涉及一种去除脑电信号中眼电伪迹的方法，主要应用于大脑工作记忆、脑电信号特征提取以及辅助疾病诊断治疗等等。

背景技术

脑电信号是由大脑神经细胞产生的，反映大脑活动的生物电信号。由于脑电具有容易获取，非侵入性，很高的时间分辨率的特点，所以脑电在科学研究和疾病诊断方面发挥着越来越重要的作用。但是脑电信号非常微弱，幅值很小，而且随机性强，在采集信号的过程中受到各种噪声的干扰，比如：眼电，肌电，心电，工频等等。因此，采集到的脑电信号中包含多种伪迹。眼电作为常见的一种伪迹，幅值远远大于脑电的幅值，而且能量主要集中在低频段，这严重影响了脑电基本节律波中的δ波(0.5-4Hz)和θ波(4-8Hz)。为了消除眼电伪迹的影响，在临床上，医生通常将含有眼电干扰的脑电数据段舍弃掉，这样可能导致一些重要脑电信息的丢失。显然，这种方法是不可取的。所以，眼电伪迹的去除，一直是脑电信号预处理的非常重要的研究内容。

传统的滤波器在消除工频干扰和其他高频伪迹有很好的效果，可是眼电和脑电信号频谱有部分重叠，所以在去掉眼电伪迹的同时势必 造成脑电信息的损失。自回归方法也一直用于去除眼电伪迹，其基本要求需要一个好的回归导联(例如眼电参考导联)。因为眼电和脑电导联交叉影响，头皮电极采集到的眼电和脑电信号不可能是纯净信号，这样会过高估计眼电伪迹，导致校正后的脑电信号丢失一些有用的脑电信息。小波变换是由傅里叶变换发展而来，具有时频局部化和多分辨率特性，很适合对微弱的脑电信号进行去噪。但是，在去噪之前，需要大量的实验去选择适当的小波基函数和分解层数，不仅消耗大量时间，也增加了计算复杂度。独立分量分析(independent component analysis,ICA)方法是九十年代发展起来的一项新的信号分解技术，用来从混合信号中提取具有统计独立性的成分。由于ICA方法分解出来的各个独立分量之间是相互独立的，所以其在信号处理领域得到很好的应用。因为眼电和脑电是由不同的信号源产生，相互独立，这样就可以利用ICA把眼电从脑电信号中分离出来，从而把眼电伪迹消除。但是在识别伪迹的时候，通常是根据脑电波形与地形图去判定，直接去除掉与眼电相关的独立分量会损失部分脑电信息，并且该方法非常耗时，容易使人疲劳，不适合实时处理脑电信号。多元经验模式分解(multivariate empirical mode decomposition,MEMD)方法是最近这些年发展起来的，是一种非平稳信号分解方法，该方法事先不需要选定基函数，其最大特点是依据数据本身的时域信息进行时域分解，得到的多元本征模式函数(multiple multivariate intrinsic mode functions,MIMFs)通常个数是有限和平稳的，而且是具有实际意义的窄带信号。因此该方法也被逐渐应用到脑电伪迹去除方面。但是，直接将与眼电 伪迹相关的MIMFs去掉，可能导致损失大量的脑电信息，显然这是不可取的。

发明内容

为了解决以上方法存在的问题，本发明的目的在于提供一种去除脑电信号中眼电伪迹的方法，结合ICA和MEMD这两种方法，MEMD从脑电信号中完全提取出眼电伪迹，并且保留了部分脑电信息，这样去除眼电的时候，避免了部分脑电信息的损失；ICA可以从脑电信号中分离出来眼电伪迹，通过计算独立分量的四阶累积量判定伪迹，然后去除眼电；本发明结合这两种方法的优势用于眼电伪迹的去除，不仅能够自动识别眼电伪迹并将其去除，同时保留了大量脑电信息，为脑电去噪提供了一种新思路。

为了达到上述目的，本发明的技术方案为：

一种去除脑电信号中眼电伪迹的方法，包括以下步骤：

步骤一、对含眼电伪迹的脑电信号X(t)进行MEMD处理，设含眼电伪迹的脑电信号X(t)＝[x₁(t),x₂(t),…,x_n(t)]^T，其中n表示脑电信号的通道个数，t表示时间，T表示矩阵的转置，脑电信号经MEMD分解后产生m个MIMFs即其中每个MIMF_i均为一个n通道的信号，且各通道的信号都处于同一频段内，同时分解产生的全部MIMF_i(i＝1,2,…,m)的频率范围是按照从高到低的顺序排列的；

步骤二、计算每个MIMF(t)的功率谱以及眼电参考信号v(t)的功率谱，选取功率集中频率范围处于参考眼电的功率集中频率范围之内的MIMFs进行重构，具体为：若第j个MIMF满足功率集中频率范围 处于参考眼电的功率频率范围之内的要求，则将第j个MIMF直到最后一个MIMF进行重构即 $U\left(t\right)={[u_1\left(t\right),u_2\left(t\right),...,u_n\left(t\right)]}^T={\mathrm{\Σ}}_{i=j}^m{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，U(t)表示已经完全包含了眼电伪迹的信号分量，同时保留剩下的MIMF_k(k＝1,2,…,j-1)用于后续处理；

步骤三、对重构的信号U(t)进行ICA处理，按照统计独立的原则找到估计分离矩阵，然后将重构信号分解出n个独立分量即C(t)＝[c₁(t),c₂(t),...,c_n(t)]^T＝W*U(t)，其中c₁(t)_,c₂(t)_,…_,c_n(t)代表n个独立分量，W表示分离矩阵，t表示时间，T表示矩阵的转置，然后计算每个独立分量的四阶累积量，将四阶累积量大于阈值的独立分量判定为伪迹并置零，阈值设定为1.5，其它独立分量保持不变；

步骤四、信号的重构，对步骤三处理后的独立分量进行ICA逆变换即H(t)＝[h₁(t),h₂(t),...,h_n(t)]^T＝W^-1*C(t)，其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，W^-1为W的逆矩阵，C(t)中含眼电伪迹的独立分量已经置零，通过逆变换得到去除眼电伪迹的脑电信号H(t)，再与步骤二保留的MIMFs进行MEMD逆变换即  $Y\left(t\right)={[y_1\left(t\right),y_2\left(t\right),...,y_n\left(t\right)]}^T=H\left(t\right)+{\mathrm{\Σ}}_{i=j}^{j-1}{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，最终得到去除眼电伪迹同时保留了大部分信息的脑电信号Y(t)。

本发明的优势在于：本发明提出了一种基于独立分量分析与多元经验模式分解的脑电信号中眼电伪迹去除的方法即IMEMD方法，该方法不仅能够有效地去除眼电伪迹，而且保留了大量有用的脑电信 息，通过与单独的ICA方法进行对比，结果表明IMEMD方法去噪效果更好，能够进一步提高信号信噪比和减少信号均方误差。

附图说明

图1是本发明的流程图。

图2(a)是仿真脑电信号；图2(b)是仿真眼电信号；图2(c)是仿真混合信号,由脑电信号和眼电信号混合而成，眼电参考信号为眼电信号幅值最大通道，仿真信号都是在MATLAB环境下产生。

图3是混合信号经MEMD处理，得到12个MIMF，图中显示了前4个MIMF的波形图，图3(a)～(d)分别表示MIMF₁，MIMF₂，MIMF₃，MIMF₄的波形图。

图4是MIMF₁，MIMF₂，MIMF₃，MIMF₄以及眼电的功率谱图，其中黑色线条表示眼电功率谱。

图5(a)是由MIMF₄～MIMF₁₂重构的信号,图5(b)是MIMFs重构的信号经ICA处理产生的独立分量，图5(c)是ICA逆变换，图5(d)是MEMD逆变换。

图6是含有眼电伪迹的真实脑电信号以及眼电参考导联，其中最后一个波形图表示眼电参考信号。

具体实施方式

下面结合附图对本发明做详细叙述。

参照图1，一种去除脑电信号中眼电伪迹的方法，包括以下步骤：

步骤一、对含眼电伪迹的脑电信号X(t)进行MEMD处理，相对于 经验模式分解(empirical mode decomposition,EMD)方法而言，MEMD能够对多通道脑电信号同时进行处理，可以使产生的各个MIMF中所有通道的信号都处于同一频段之内，以便进行同频带的后续处理，设含眼电伪迹的脑电信号X(t)＝[x₁(t),x₂(t),…,x_n(t)]^T，其中n表示脑电信号的通道个数，t表示时间，T表示矩阵的转置，脑电信号经MEMD分解后产生m个MIMFs即其中每个MIMF_i均为一个n通道的信号，且各通道的信号都处于同一频段内，同时分解产生的全部MIMF_i(i＝1,2,…,m)的频率范围是按照从高到低的顺序排列的；

步骤二、计算每个MIMF(t)的功率谱以及眼电参考信号v(t)的功率谱，选取功率集中频率范围处于参考眼电的功率集中频率范围之内的MIMFs进行重构，具体为：若第j个MIMF满足功率集中频率范围处于参考眼电的功率频率范围之内的要求，则将第j个MIMF直到最后一个MIMF进行重构即 $U\left(t\right)={[u_1\left(t\right),u_2\left(t\right),...,u_n\left(t\right)]}^T={\mathrm{\Σ}}_{i=j}^m{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，U(t)表示已经完全包含了眼电伪迹的信号分量，同时保留剩下的MIMF_k(k＝1,2,…,j-1)用于后续处理；

步骤三、对重构的信号U(t)进行ICA处理，按照统计独立的原则找到估计分离矩阵，然后将重构信号分解出n个独立分量即C(t)＝[c₁(t),c₂(t),...,c_n(t)]^T＝W*U(t)，其中c₁(t)_,c₂(t)_,…_,c_n(t)代表n个独立分量，W表示分离矩阵，t表示时间，T表示矩阵的转置，然后计算每个独立分量的四阶累积量，四阶累积量俗称峭度，用来判断信 号的高斯性，若信号具有超高斯特性即信号有尖峰分布特点，比如眼电、心电信号等，则峭度值就大于零，将四阶累积量大于阈值的独立分量判定为伪迹并置零，阈值设定为1.5，其它独立分量保持不变；

步骤四、信号的重构，对步骤三处理后的独立分量进行ICA逆变换即H(t)＝[h₁(t),h₂(t),...,h_n(t)]^T＝W^-1*C(t)，其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，W^-1为W的逆矩阵，C(t)中含眼电伪迹的独立分量已经置零，通过逆变换得到去除眼电伪迹的脑电信号H(t)，再与步骤二保留的MIMFs进行MEMD逆变换即  $Y\left(t\right)={[y_1\left(t\right),y_2\left(t\right),...,y_n\left(t\right)]}^T=H\left(t\right)+{\mathrm{\Σ}}_{i=j}^{j-1}{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，最终得到去除眼电伪迹同时保留了大部分信息的脑电信号Y(t)。

下面分别对仿真信号和真实信号中的眼电伪迹进行去除，在仿真信号去噪的过程中，同时与单独ICA方法伪迹去除结果进行对比。

一、仿真信号中眼电伪迹的去除

步骤1.构造仿真脑电信号A(i)＝[a₁(i),a₂(i),a₃(i),a₄(i),a₅(i)]^T(i＝1,2,…,1600),和仿真眼电信号B(i)＝[b₁(i),b₂(i),b₃(i),b₄(i),b₅(i)]^T(i＝1,2,…,1600),其中i表示信号采样点序号，T表示矩阵的转置，信号通道个数都为5，信号采样率200Hz，数据长度8秒，脑电信号和眼电信号波形如图2(a)～(b)所示。将脑电信号和眼电信号按照以下方式混合：

C(i)＝[c₁(i),c₂(i),c₃(i),c₄(i),c₅(i)]^T＝A(i)+λB(i),i∈{1,2,...,1600}  (1) 

其中i表示信号采样点序号，T表示矩阵的转置。通过调整参数λ，可以得到不同信噪比下的混合信号。仿真混合信号参见图2(c)，可以 看到混合信号中含有明显的眼电伪迹，相对于脑电信号，眼电的幅值很大。然后对附图中构造的仿真混合信号进行MEMD处理，通过处理将信号分解成频率从高到低排列的12个MIMFs即  $C\left(i\right)={[c_1\left(i\right),c_2\left(i\right),c_3\left(i\right),c_4\left(i\right),c_5\left(i\right)]}^T={\mathrm{\Σ}}_{j=1}^{12}{\mathrm{MIMF}}_j\left(i\right)(i=\mathrm{1,2},...,1600),$ 其中i表示信号采样点序号，T表示矩阵的转置。参见图3，画出了前4个MIMF的波形图，明显可看出MIMF₁幅值很大，功率集中的频率最高，而MIMF₂，MIMF₃以及MIMF₄幅值逐渐减小，并且功率集中的频段也同样降低。

步骤2.使用Welch方法计算步骤1中的每一个MIMF以及参考眼电信号的功率谱，这里的眼电信号是选取B(i)中幅值最大的眼电通道作为参考，如图4所示,其中黑色线条表示眼电功率谱。选取功率集中频率范围在眼电功率集中频率范围的MIMFs进行重构，由图可看到，眼电的功率主要集中在3Hz以下，MIMF₁的功率集中在20Hz以上，同样MIMF₂和MIMF₃功率集中的频段也不在3Hz以下，而MIMF₄有部分功率集中的频段在眼电功率集中频段范围内，所以MIMF₄满足要求，由于MIMFs是按频率从高到低排列的，所以排序在MIMF₄后面的MIMFs也都满足要求。接着将满足要求的MIMFs进行重构即  $U\left(i\right)={[u_1\left(i\right),u_2\left(i\right),u_3\left(i\right),u_4\left(i\right),u_5\left(i\right)]}^T={\mathrm{\Σ}}_{j=1}^{12}{\mathrm{MIMF}}_j\left(i\right)(i=\mathrm{1,2},...,1600),$ 其中i表示信号采样点序号，T表示矩阵的转置。重构信号U(i)如图5(a)，可看到重构信号含有明显的眼电伪迹，同时保留前三个MIMF_j(j＝1,2,3)，通过与图2(c)中的混合信号对比，可发现重构信号已经将眼电伪迹完全提取出来。

步骤3.使用ICA方法对步骤2重构的信号U(i)进行处理即D(i)＝[d₁(i),d₂(i),d₃(i),d₄(i),d₅(i)]^T＝W*U(i)(i＝1_,2,…,1600),其中d₁(i),d₂(i),d₃(i),d₄(i),d₅(i)代表5个独立分量，i表示信号采样点序号，W表示分离矩阵，T表示矩阵的转置。独立分量D(i)参见图5(b)。然后计算每个独立分量的四阶累积量，经计算它们的四阶累积量分别为1.79，-0.45，-0.44，-0.46，-0.44，将四阶累积量大于阈值的分量判定为伪迹并置零(本发明设定的阈值为1.5)，所以将第一个独立分量置零，其他的独立分量保持不变。

步骤4.对步骤3处理后的独立分量进行ICA逆变换即H(i)＝[h₁(i),h₂(i),h₃(i),h₄(i),h₅(i)]^T＝W^-1*D(i)(i＝1_,2,…,1600),其中i表示信号采样点序号，W^-1为W的逆矩阵，T表示矩阵的转置。经过逆变换得到的脑电信号H(i)如图5(c)，可以看到眼电伪迹已经被完全去除掉，再将与信号H(i)步骤2保留的MIMFs进行MEMD逆变换即  $Y\left(i\right)={[y_1\left(i\right),y_2\left(i\right),y_3\left(i\right),y_4\left(i\right),y_5\left(i\right)]}^T={\mathrm{\Σ}}_{j=1}^3{\mathrm{MIMF}}_j\left(i\right)+H\left(i\right)(i=\mathrm{1,2},...,1600),$ 其中i表示信号采样点序号，T表示矩阵的转置。将没有参加ICA处理的脑电信息重构回去，减少了脑电信息的损失，最终得到去除眼电伪迹的脑电信号Y(i)如图5(d)所示。

为了定量评价本发明去噪效果好坏，在仿真实验中，与单独ICA方法处理结果进行比较。使用信噪比(SNR)和均方误差(MSE)作为指标。信噪比的公式如下：

$\mathrm{SNR}=10*\lg ({\mathrm{\Σ}}_{i=1}^n{x}_i^2/{\mathrm{\Σ}}_{i=1}^n{y}_i^2)---\left(2\right)$

均方误差的表达式如下：

$\mathrm{MSE}=\frac{1}{n}{\mathrm{\Σ}}_{i=1}^n{(x_i-y_i)}^2---\left(3\right)$

根据公式(1)得到四种信噪比下90组脑电信号，然后计算脑电信号去噪后所有通道SNR和MSE的均值和标准差，结果如表1和表2所示。

表1不同信噪比下，IMEMD和ICA去噪后信噪比的均值和标准差

表2不同信噪比下，IMEMD和ICA去噪后均方误差的均值和标准差

从表1和表2可以看出IMEMD方法在提高信号信噪比和降低均方误差都要优于ICA方法，说明前者不仅能够有效地去除掉眼电伪迹，相对于后者，损失的脑电信息更少。通过作统计分析，两者具有显著性差异(p<0.05)。

二：真实脑电信号中眼电伪迹的去除

为了进一步验证本发明的有效性，采集了三个被试脑电信号以及眼电参考导联信号，脑电采样率500Hz。每个被试选取了含有眼电伪迹明显的脑电信号30段，本发明对伪迹进行去除，然后计算去噪前后每个通道脑电信号与眼电参考信号的互相关系数Coe。一段受到眼电干扰的脑电信号，可看到越靠近眼部的导联受到的干扰越强，远离眼部的导联受到伪迹干扰逐渐减弱，其中最后一个通道是眼电参考，眼电幅值比脑电信号幅值大得多。使用IMEMD方法对脑电信号去除伪迹，可发现眼电伪迹几乎完全去除。处理结果如表3所示。

表3各脑电通道去噪前后与眼电参考信号的互相关系数Coe

从表3可看出，相对于去噪前的互相关系数，去噪后的互相关系数明显降低，通过作统计分析，去噪前后的互相关系数有显著性差异(p<0.05)，表明本发明能够彻底地将眼电伪迹从脑电信号中去除。

Claims (1)

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

步骤一、对含眼电伪迹的脑电信号X(t)进行MEMD处理，设含眼电伪迹的脑电信号X(t)＝[x₁(t),x₂(t),…,x_n(t)]^T，其中n表示脑电信号的通道个数，t表示时间，T表示矩阵的转置，脑电信号经MEMD分解后产生m个MIMFs即其中每个MIMF_i均为一个n通道的信号，且各通道的信号都处于同一频段内，同时分解产生的全部MIMF_i(i＝1,2,…,m)的频率范围是按照从高到低的顺序排列的；

步骤二、计算每个MIMF(t)的功率谱以及眼电参考信号v(t)的功率谱，选取功率集中频率范围处于参考眼电的功率集中频率范围之内的MIMFs进行重构，具体为：若第j个MIMF满足功率集中频率范围处于参考眼电的功率频率范围之内的要求，则将第j个MIMF直到最后一个MIMF进行重构即 $U\left(t\right)=[u_1\left(t\right),u_2\left(t\right),...,u_n\left(t\right){]}^T={\mathrm{\&Sigma;}}_{i=j}^m{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，U(t)表示已经完全包含了眼电伪迹的信号分量，同时保留剩下的MIMF_k(k＝1,2,…,j-1)用于后续处理；

步骤三、对重构的信号U(t)进行ICA处理，按照统计独立的原则找到估计分离矩阵，然后将重构信号分解出n个独立分量即C(t)＝[c₁(t),c₂(t),...,c_n(t)]^T＝W*U(t)，其中c₁(t),c₂(t),…,c_n(t)代表n个独立分量，W表示分离矩阵，t表示时间，T表示矩阵的转置，然后计算每个独立分量的四阶累积量，将四阶累积量大于阈值的独立分量判定为伪迹并置零，阈值设定为1.5，其它独立分量保持不变；

步骤四、信号的重构，对步骤三处理后的独立分量进行ICA逆变换即H(t)＝[h₁(t),h₂(t),...,h_n(t)]^T＝W^-1*C(t)，其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，W^-1为W的逆矩阵，C(t)中含眼电伪迹的独立分量已经置零，通过逆变换得到去除眼电伪迹的脑电信号H(t)，再与步骤二保留的MIMFs进行MEMD逆变换即 $Y\left(t\right)=[y_1\left(t\right),y_2\left(t\right),...,y_n\left(t\right){]}^T=H\left(t\right)+{\mathrm{\&Sigma;}}_{i=1}^{j-1}{\mathrm{MIMF}}_i\left(t\right),$ 其中n表示信号的通道个数，t表示时间，T表示矩阵的转置，最终得到去除眼电伪迹同时保留了大部分信息的脑电信号Y(t)。