脑电中眼电伪迹的自动识别与去除

脑电中眼电伪迹的自动识别与去除李明爱;刘帆【摘要】目的为改善脑电中眼电伪迹的去除效果,基于脑电的非平稳性和模糊特点,提出一种将离散小波变换与二阶盲辨识相结合，并以模糊熵为眼电伪迹判别准则的眼电伪迹去除方法.方法首先,采用离散小波变换对含噪的脑电信号进行多分辩分析,获得平稳性更好的多尺度小波系数;进而，选择同层的小波系数构成小波系数矩阵，并基于二阶盲辨识对其盲源分离,得到源信号的估计;进一步以模糊熵为判别依据，实现眼电伪迹的自动判别与剔除.实验数据采用BCICompetitionIV公开数据库，使用信噪比、相关系数及均方根误差等常用伪迹判别指标进行衡量.结果本文方法相对于常用的眼电伪迹去除方法在多个性能指标上均取得最大值.结论本文提出的眼电伪迹去除方法，实现了眼电伪迹的自动精确判断与剔除，并表现出很好的稳定性.【期刊名称】《北京生物医学工程》【年(卷)，期】2018(037)006【总页数】7页(P559-565)【关键词】脑电信号;眼电伪迹;离散小波变换;二阶盲辨识;模糊熵【作者】李明爱浏帆【作者单位】北京工业大学信息学部北京100124;计算智能与智能系统北京市重点实验室北京100124;北京工业大学信息学部北京100124【正文语种】中文【中图分类】R318.040引言脑电信号(electroencephalogram,EEG)是大脑中的电生理活动在大脑皮质或头皮表面神经细胞的整体反映，其蕴含丰富的生理、心理及病理信息，在神经科学、临床医学、康复工程等研究领域的作用日益显著。但EEG信号微弱，在采集过程中不可避免地受到眼电(electrooculogram,EOG)、心电(electrocardiogram,ECG)、肌电(electromyogram,EMG)等噪声干扰，特别是眼电信号的巨大幅值使得眼电伪迹(oculararitifial,OA)对EEG信号产生严重影响。如何有效去除眼电伪迹，获得纯净的脑电信号具有重要的研究价值[1-2]。脑电信号的伪迹去除方法主要有两类，小波变换法[3](wavelettransform,WT)和盲源分离法[4-5](blindsourcesseparation,BSS)。小波变换法利用WT具有良好的时频特性和自适应特点，将含眼电伪迹的EEG源信号进行多尺度分解，通过将与眼电伪迹频带吻合的部分多尺度信号置零实现伪迹去除。由于EEG和EOG存在频带重叠现象，致使伪迹去除效果并不理想。盲源分离法利用BSS将含噪声的EEG观测信号分解成多个源信号分量，通过判断并剔除与眼电伪迹相似或相关的部分信号分量，还原出纯净的脑电信号。常用的BSS包括：(1)主成分分析(principlecomponentanalysis,PCA)[6]O该方法假设源信号是非相关的，基于协方差矩阵对观测信号进行正交分解，得到一组线性无关的变量，即主成分。当眼电信号与脑电信号成分具有相当的幅值时，不能彻底去除眼迹成分，无法得到理想的去噪效果。(2)独立成分分析(independentcomponentanalysis,ICA)[7]OICA假设源信号是统计独立的，且至多一个服从高斯分布，基于信号的二阶及高阶统计信息进行信号分解，从混合信号中恢复出原始信号的各统计独立分量，计算量较大，且EEG难以满足假设条件，去噪效果欠佳。FastICA是一种典型的ICA的快速实现算法，计算效率得到很好提升，但易受噪声干扰，且存在迭代运算不收敛现象。⑶二阶盲辨识(secondorderblindidentification,SOBI)［8-9］。该方法假设源信号为不相关的平稳信号，利用观测信号的二阶统计量实现源信号的盲分离，不仅减少了计算量，而且增强了方法的鲁棒性。脑电信号的统计特性与各种BSS方法假设条件的吻合程度是影响基于BSS的伪迹去除方法性能的重要因素。为了进一步改善眼迹去除效果，利用各种方法的优势，将两种及以上去噪方法相结合是很好的选择［10］。文献［11］提出离散小波变换(discretewavelettransform,DWT)和FastICA相结合的伪迹去除方法，即DWICA，它将EEG信号进行小波变换，获得比原始信号超高斯性更强、峰度更大的小波系数，进而对其进行独立成分分析，使收敛速度和抗噪能力同时得以提高。影响眼迹去除效果的另一重要因素是伪迹判断准则。基于WT的眼迹去除方法中，频率范围是衡量眼电伪迹成分的标准，由于不同受试者的眼电频率范围存在差异，且无法精确估算，与多尺度信号间频率匹配时还会有误差，导致眼电伪迹剔除不准确。基于BSS的眼迹去除方法中，通过计算眼电参考信号与信号分量间的相关系数或夹角余弦进行眼电伪迹判断。相关系数与夹角余弦准则均具有较强的分辨能力，且计算速度快，但要求脑电信号的数据序列足够长。文献［12］根据眼电伪迹相对脑电信号的复杂度更低，因而具有更小的熵值，采用样本熵来识别眼电伪迹。然而，样本熵的相似度量公式为Heaviside二值函数，其突变性使得样本熵值缺乏连续性，对阈值的取值非常敏感，其微弱变化，可能导致样本熵值的突变，造成伪迹的错误判断。本文提出将离散小波变换和二阶盲辨识算法相结合的脑电中眼电伪迹的去除方法。利用离散小波变换对含眼电噪声的非平稳脑电信号进行多尺度分析，以放大突变的眼电伪迹信号；对同一尺度下具有更好平稳性的小波系数［13］进行二阶盲辨识，并引入模糊熵作为眼电伪迹的判别依据，实现眼电伪迹的自动判别与去除。1基本原理1.1二阶盲辨识二阶盲辨识是一种基于二阶统计量的盲源分离方法。设信号线性混合模型如式(1)所示：x(t)=As(t)+n(t)(1)式中：s(t)=[s1(t),s2(t),...，sn(t)]T代表n个未知源信号；T代表矩阵的转置；x(t)=[x1(t),x2(t),...,xm(t)]T为通过传感器获得的m个观测信号；AuRmn表示线性混合矩阵，n(t)=[n1(t),n2(t),...,nm(t)]T表示噪声信号。假设源信号s(t)为多变量平稳信号，混合矩阵A为列满秩矩阵，n(t)为独立于源信号的零均值加性白噪声。则二阶盲辨识方法通过对一批协方差矩阵进行联合近似对角化来求解分离矩阵，以获得源信号s(t)的一个估计y(t)，主要步骤如下[8,14]:计算白化矩阵W,将观测信号x(t)白化处理，使得z(t)的协方差矩阵是单位矩阵，除去各分量之间的相关性。z(t)=Wx(t)(2)对于固定延时TU｛Tj|j=12...,k｝，计算白化数据的采样协方差矩阵R(t):R(T)=E[z(t+T)zT(t)]=ARz(T)AT⑶对所有的R(Tj),采用联合近似对角化法，求得正交归一化矩阵V:VTR(Tj)V=Dj(4)式中：｛Dj｝表示一组对角矩阵。计算混合矩阵A=W#V，求得去除相关后的源信号的估计y(t):y(t)=A-1x(t)=VTWx(t)⑸1.2离散小波变换小波变换是时间、频率的局部化分析方法，通过伸缩平移实现对信号的逐步细化，反映信号在时域和频域的局部特征。因具有多分辨分析和对信号的自适应特点，小波变换特别适合于处理EEG这类非平稳信号。在实际中，为方便计算，减少运算量，通常将处理的信号离散化，此即为离散小波变换［15］。在多分辨分析基础上，Mallat提出离散小波快速算法，即Mallat算法［16］。利用Mallat算法可以根据不同尺度进行多级分解，把源信号分解成细节分量和逼近分量，如式(6)所示。(6)式中：an,k、dn,k分别代表第n层的逼近分量和细节分量。Mallat塔式重构算法是分解算法的逆过程，如式(7)所示，当n逐级递减至1时，完成对信号的重构。⑺1.3模糊熵模糊熵(fuzzyentropy，FE)［17］是一种时间序列复杂度的测量方法。熵值越小，序列越简单；熵值越大，序列越复杂。针对眼电伪迹频率低、复杂性低，其模糊熵值相对脑电信号更小的特点，通过计算FE可实现伪迹判别，计算式如式(8)。(8)式中：em(n,r)由相似度计算得到；m为嵌入维数；n为指数边界的梯度；r为相似容限。具体过程详见文献[16]。2基于二阶盲辨识和模糊熵的眼电伪迹自动去除方法本文提出一种眼电伪迹自动去除方法，将二阶盲辨识方法与离散小波变换相结合，实现测量信号在小波域的盲源分离，并采用模糊熵进行眼电伪迹判别，具体步骤如下。2.1基于sym8小波基函数，对原始信号进行离散小波变换分解假设x(t)=[x1(t),x2(t),...,xm(t)]T表示通过电极采集的m导信号，xi(t)为第i导脑电信号，i=12...,m-1,xm(t)为一导参考眼电信号。采用sym8小波母函数，对m导测量信号分别进行3层离散小波分解，并将得到的逼近分量和细节分量构成小波系数向量Cx,i(t)=[Ai3Di3Di2Di,1](i=12...,m),进而构造小波系数矩阵Cx(t)=[Cx,1(t),Cx,2(t),...,Cx,m(t)]T。2.2基于SOBI算法的小波系数矩阵盲源分离采用二阶盲辨识算法，对小波系数矩阵进行盲源分离，求得其源信号的估计，记为Cs(t)=[Cs,1(t),Cs,2(t),...,Cs,m(t)]T。2.3基于模糊熵的眼电伪迹判别与去除计算各个源估计信号Cs,i(t)的模糊熵值FE(Cs,i(t))(i=12...,m),对其进行排序，其中具有最小模糊熵值的源估计信号即自动判定为眼电伪迹。将其置零，其余源估计信号保持不变，以达到伪迹剔除的目的。2.4基于SOBI投影逆变换假设眼电伪迹去除后的源估计信号记为则基于SOBI算法对其投影变换后的信号为2.5小波的逆变换基于式(7)所示Mallat塔式重构算法，对信号进行小波逆变换，还原出伪迹去除后的头皮脑电信号xc(t)。3实验研究与结果分析3.1实验数据3.1.1数据来源实验数据来自公开的BCICompetitionIV的DataSets2b数据库。采用Ag/AgCl电极，电极位置服从国际标准的10~20导联系统分布［18］。数据记录9名右利手、视觉良好的健康受试者。每个受试者的实验数据采集包含3个阶段，如图1所示。(1)空载阶段，包括2min睁眼和1min闭眼。(2)眼电采集阶段80s，包括眨眼、转动眼球、上下眼动和左右眼动各15s，受试者根据屏幕上的提示，依次进行。每种眼电伪迹采集完成时，随后休息5s进入下一阶段。(3)脑电采集阶段，采集左、右手运动想像EEG各60次共计120组数据。选择C3、Cz和C4,3个电极位置记录脑电数据(Cz为C3和C4的参考电极)，采样频率为250Hz,并将采集的信号通过0.5~100Hz的带通滤波及50Hz的陷波滤波器滤除工频干扰。采集实验时序如图1所示。图1实验采集时序Figure1Thetimingdiagramofcollectionexperiment3.1.2实验数据构造为了精确衡量脑电信号中眼电伪迹的去除效果，这里将根据脑电信号和眼电信号的双向激活扩散特性［19］，利用纯净的眼电信号和脑电信号构造含噪声的眼电信号和脑电信号。模型表达式如下：EOGraw(t)=EOGclean(t)+n1EEGclean(t)(9)EEGraw(t)=EEGclean(t)+n2EOGclean(t)+WGN(10)式中:EOGraw(t)、EEGraw(t)分别代表实际采集到的含噪声的眼电信号和脑电信号;EOGclean(t).EEGclean(t)代表纯净的眼电信号和脑电信号;WGN为高斯白噪声；系数n1和n2为干扰影响系数，n1u(0,0.5］,n2u(0,0.2］,干扰系数均随机产生，以保证实验结果的可信度。同时，为模拟肌电信号、心电信号等其他噪声对脑电信号的干扰，在C3，Cz，C4三导脑电信号中加入10dBW的高斯白噪声。3.2眼电伪迹去除性能评价指标3.2.1均方根误差均方根误差(rootmeansquarederror，RMSE)是一种衡量伪迹去除效果最常用的性能指标。假设xclean(t)和xc(t)分别表示某一导联上的纯净脑电信号和伪迹去除后的脑电信号，则均方根误差计算式为：(11)式中:N为采样点个数。RMSE的值越小，说明伪迹去除后的EEG越接近纯净的EEG，即去除效果越好。3.2.2信噪比第二种常用的伪迹去除效果评价指标为信噪比(signalnoiserate，SNR)，其计算式为：(12)xclean(t)和xc(t)式中标量含义同上。SNR的值越大，意味着xclean(t)和xc(t)越接近，说明去噪效果越好。3.2.3相关系数相关系数(correlationcoefficient,CC)是用以反映变量之间相关关系密切程度的统计指标，是用于评价伪迹去除效果的又一常用指标，计算式为：(13)式中：分别表示信号xclean(t)和xc(t)的均值。CC数值越大，表示xclean(t)和xc(t)之间相似度越高，即去噪效果越理想。3.3实验结果与分析3.3.1眼电伪迹去除的实验研究基于数据库中某受试者(subject1)的水平眼电(horizontalEOG,H-EOG)和运动想像脑电实验数据，将本文方法与常用眼电伪迹去除方法进行对比实验研究，并使用多种性能指标进行性能评价［9］。为克服实验数据的偶然性对实验结果的影响，依式(9)和式(10)构造含噪眼电和运动想像脑电的实验数据时，相互干扰系数n1、n2分别随机产生10组。针对每组干扰系数，从120组运动想像脑电数据中任选一组与水平眼电分别截取等长的750个采样点，构造含噪的运动想像脑电和眼电信号，并构成4导源信号进行实验研究。采用WT、PCA、ICA等多种方法去除眼电伪迹，计算伪迹去除后的脑电信号与纯净脑电信号的10组实验均方根误差的平均值如表1所示，WT方法去除水平眼电伪迹时效果最差，这是由于WT方法无法针对EEG信号与EOG信号的频带混叠问题进行伪迹有效去除；基于PCA、ICA、SOBI的水平眼电伪迹去除效果得到明显改善，且ICA和SOBI的效果较为接近并优于PCA,这说明脑电源信号和观测信号能够较好地满足盲源分离的原理和条件，但眼电伪迹去除效果还取决于脑电源信号的统计特点和分布特性；基于DWICA的去除效果进一步得到提高，这是由于小波系数相对于EEG源信号超高斯性更强；本文方法(DWSOFE)去除水平眼电伪迹效果最佳，在C3、Cz、C4中任意导联上均方根误差都取得相对最小值，这主要是由于该方法利用小波变换对非平稳的MI-EEG信号进行多分辨分析，获得了具有较好平稳性的小波系数，使之更加吻合SOBI假设条件，有效增强了对于EEG信号的适应性。表1多种方法去除H-EOG后的平均RMSE比较(单位：^V)Table1AverageRMSEcomparisonofmultiplemethodsafterH-EOGremoval(unit:pV)导联WTPCAICASOBIDWICADWSOFEC320.4111.797.247.367.176.40Cz26.0413.978.298.958.237.41C421.2912.007.607.076.976.08多种方法用于去除水平眼电在各导联上的平均信噪比和相关系数如图2和图3所示。从图2可知，WT方法获得的平均信噪比最低，ICA和PCA方法获得的平均信噪比依次得到极大提高，SOBI和DWICA的效果进一步得以改善，且效果相当，而本文方法得到最高平均信噪比。另外，从图3可见，WT方法在三导联上的平均相关系数指标均取得最小值，PCA和SOBI得到的平均相关系数依次得以明显提升，ICA和DWICA以微弱优势得以进一步提高，DWSOFE获得最高平均相关系数。充分展现了DWSOFE在以相关系数为伪迹去除性能指标时的绝对优势和更好的结果稳定性。综上可见，采用多种指标衡量水平眼电伪迹去除效果时，不同方法的表现并不具有完全一致性。例如，当采用RMSE和CC性能指标时，ICA效果好于PCA，而当采用SNR指标时，PCA却优于ICA，这与脑电源信号的统计特点以及分布特性密切相关。但对3种性能指标而言，本文方法均获得最佳效果，足以证明本文方法去除水平眼电伪迹的有效性和稳定性。图2多种方法去除H-EOG后的平均SNR比较Figure2AverageSNRcomparisonofmultiplemethodsafterH-EOGremoval3.3.2眼电伪迹判别方法对伪迹去除效果影响的实验研究表2眼电伪迹判别方法对伪迹去除效果的影响Table2ComparisonofOAremovaleffectfordifferentmethods评价指标RMSE/pVSNR/dBCC判别方法C3CzC4C3CzC4C3CzC4相关系数9.17610.6499.11511.75812.89412.3270.9220.9140.924夹角余弦7.0427.8306.74513.23614.50613.8830.9640.9690.969样本熵7.0387.8266.74013.24014.51213.8890.9650.9690.970模糊熵6.4007.4106.08014.09014.79014.7400.9720.9740.976图3多种方法去除H-EOG后的平均CC比较Figure3AverageCCcomparisonofmultiplemethodsafterH-EOGremoval脑电信号与眼电信号的复杂度不同，且脑电信号具有模糊特性，为此，本文通过计算模糊熵值来达到伪迹辨别的目的。采用3.3.1节中构造的含噪声的眼电数据和运动想像脑电数据，依据本文提出的SOBI和WT相结合的伪迹去除方法，仅改变眼电伪迹的判别准则，实现模糊熵和相关系数、夹角余弦、样本熵等常用伪迹判别准则的对比研究，基于3种评价指标的平均实验结果如表2所示。表2表明，以相关系数进行眼电伪迹效果相对较差，夹角余弦与样本熵的实验结果有明显改善，且二者效果较为接近，而模糊熵则对3种评价指标均取得最好的效果。究其原因，在模糊熵的相似度量计算式中采用指数函数，使熵值变化更加平稳平滑，不易造成伪迹判别错误，因而有利于提高伪迹识别的准确性。4讨论与结论本文提出一种二阶盲辨识结合离散小波变换的眼电伪迹去除方法，并引入模糊熵实现伪迹自动判别的同时，有效增强了伪迹判别的准确性和伪迹去除的精度。基于国际公开脑电和眼电数据库进行实验研究，结果表明，采用均方根误差、信噪比和相关系数等多项性能指标时，本文方法在去除水平眼电伪迹时均有明显优势，对脑电信号的在线预处理具有重要意义。在未来的工作中，将加强本文方法用于去除多个受试者垂直眼电、眼睛转动等多种眼电伪迹的适用性研究，并从统计学的角度验证本文方法的显著性。参考文献【相关文献】BurgerC,vandenHeeverDJ.RemovalofEOGartefactsbycombiningwaveletneuralnetworkandindependentcomponentanalysis[J].BiomedicalSignalProcessingandControl,2015,15:67-79.QnarS,AcirN.AnovelsystemforautomaticremovalofocularartefactsinEEGbyusingoutlierdetectionmethodsandindependentcomponentanalysis[J].ExpertSystemswithApplications,2017,68:36-44.KrishnaveniV,JayaramanS,AravindS,etal.AutomaticidentificationandremovalofocularartifactsfromEEGusingwavelettransform[J].MeasurementScienceReview,2008,6(4):45-57.OosugiN,KitajoK,HasegawaN,etal.AnewmethodforquantifyingtheperformanceofEEGblindsourceseparationalgorithmsbyreferencingasimultaneouslyrecordedECoGsignal[J].NeuralNetworkstheofficialjournaloftheInternationalNeuralNetworkSociety,2017,93:1-6.FitzgibbonSP,PowersDM,PopeKJ,etal.RemovalofEEGnoiseandartifactusingblindsourceseparation[J].JournalofClinicalNeurophysiology:officialpublicationoftheAmericanElectroencephalographicSociety,2007,24(3):232-243.CasarottoS,BianchiAM,CeruttiS,etal.Principalcomponentanalysisforreductionofocularartefactsinevent-relatedpotentialsofnormalanddyslexicchildren[J].ClinicalNeurophysiology,2004,115(3):609-619.KongW,ZhouZ,HuS,etal.AutomaticanddirectidentificationofblinkcomponentsfromscalpEEG[J].Sensors,2013,13(8):10783-10801.BelouchraniA,Abed-MeraimK,CardosoJF,etal.Ablindsourceseparationtechniqueusingsecond-orderstatistics[J].IEEETransactionsonSignalProcessing,1997,45(2):434-444.BridwellDA,RachakondaS,SilvaRF,etal.Spatiospectraldecompositionofmultisubjecteeg:evaluatingblindsourceseparationalgorithmsonrealandrealisticsimulateddata[J].BrainTopography,2016,31(1):47-61.UriguenJA,Garcia-ZapirainB.EEGartifactremoval-state-of-the-artandguidelines[J].JournalofNeuralEngineering,2015,12(3):031001.LiM,CuiY,YangJ.AutomaticremovalofocularartifactfromEEGwithDWTandICAmethod[J].AppliedMathematics&InformationSciences,2013,7(2):809-816.MowlaMR,NgSC,ZilanyMSA,etal.Artifacts-matchedblindsourceseparationandwavelettransformformultichannelEEGdenoising[J].BiomedicalSignalProcessingandControl,2015,22:111-118.黎志勇，李宁.基于小波的非平稳时间序列预测方法研究[J].计算机工程与应用,2014,50(10):38-43.LiZY,LiN.Researchonnon-stationarytimeseriespredictionbasedonwavelet[J].ComputerEngineeringandApplications,2014,50(10):38-43.BrewickPT,SmythAW.Increasingtheefficiencyandefficac