基于虚拟通道ICA方法的诱发电位单导少次提取
2012-08-13王永轩邱天爽
王永轩 邱天爽 刘 蓉
(大连理工大学 电子信息与电气工程学部,大连 116024)
引言
诱发电位(evoked potential,EP)反映了中枢神经系统的感觉通路及皮层区域的神经电活动情况,是临床医学诊断神经系统损伤及病变的重要手段之一,也是实现脑-机接口(brain-computer interface,BCI)的重要脑电信号之一[1-2]。由头皮表面测得的EP信号是和自发脑电(EEG)等背景噪声同时记录到的,信噪频谱重叠,而且信噪比很低,所以提取难度较大。目前,临床医学通常使用的方法是叠加平均法[3-4],通过很多的刺激次数来得到足够的信噪比。然而,叠加平均法中太多的刺激次数会增大每次EP信号之间的差异,而且反复刺激容易引起神经疲劳甚至产生暴发伪迹,因此少次提取问题成为一个重要的研究热点和难点。此外,EP信号是由特定位置的大脑皮层细胞产生的反应,传导到头皮不同位置时相当于经过了不同的系统,很难用线性模型描述,从而 EP信号的单导提取也显得尤为重要。
独立分量分析(independent component analysis,ICA)是为解决盲源分离问题而产生的一种有效算法[5-8]。它从多路观测信号出发,对已知信息量很少的源信号进行估计,最后得到相互独立的源信号的近似值。ICA算法处理的对象是相互独立的源信号经未知线性组合而产生的一组混合信号,在实际使用时,一般假设观测信号路数大于或等于独立源的个数。显然在EP信号单导提取任务中无法利用多导信号提供的空间信息,因此不能直接使用ICA算法。为了实现这一目的,文献[9-10]均引入了一路噪声(EEG)作为虚拟通道,将单路观测信号扩展为两路,以满足ICA算法的应用条件。这个模型建立的前提必须是虚拟噪声通道与观测信号中的噪声具有很高的相关性,但实验表明EEG作为一种随机信号,不同时间段的EEG之间的相关性很弱,即虚拟噪声通道与观测信号中的噪声匹配度很低,因此这个模型很难消除EP信号中的噪声。文献[11,12]假设单路观测信号是由 EP信号和多个噪声分量叠加而成,引入的多路虚拟噪声通道分别代表了这些噪声分量。虽然模型中假设多路虚拟噪声通道是独立源信号,但是ICA算法要求输出信号最大程度符合独立性准则,所以引入的虚拟通道之间的实际相关程度会束缚ICA算法的效果。针对以上问题,本文在EP信号的单导少次提取中,将少次信噪混合信号作为虚拟观测通道,构建了新的虚拟通道模型,增加了ICA算法中解混矩阵的自由度,从而提升了ICA的应用性能。
1 独立分量分析算法
1.1 算法原理
ICA算法处理的信号是多路观测信号,将多路观测信号用向量表示为
式中,M表示观测信号路数,x的各个分量xm表示的是第m路观测信号,即
式中,N为采样点数。在ICA算法中把观测信号以及源信号都看作是随机变量,每个采样点是随机变量的不同取值,在算法中的各种统计量正是利用这些采样点来进行计算。多路观测信号的形成是多路未知源信号的未知线性混合,数学模型为
式中,s= [s1,s2,…,sM]T为未知源信号向量,A 为M×M维未知混合矩阵。ICA算法的应用前提是s的各个分量sm相互独立,且A是满秩的。在这个假设条件下根据观测信号x得到源信号的估计^s,因此^s是在独立性准则条件下的最佳估计。通过寻找一个解混矩阵BM×M,得到多路输出信号
如果y的各个分量ym在最大程度符上合独立性准则,那么就可以认为y是 s的一个最佳估计。显然如果解混矩阵B是混合矩阵A的逆矩阵,即
那么ICA算法将准确还原多路源信号。由于ICA算法假设源信号是统计独立的,因此需要给出独立性准则以及反映独立性准则的目标函数,当目标函数达到极值时,就可以认为输出信号是独立的,也就是源信号的估计。比较经典的算法包括FastICA算法和 Infomax算法等[13]。
1.2 虚拟通道
单次EP信号观测模型可表示为
式中,x1= [x1(1),x1(2),…,x1(N)]为单次含噪信号,s1为纯净EP信号,v1表示 EEG等噪声。为了应用ICA算法实现单导诱发电位提取,需要引入虚拟通道,将一路观测信号扩展为多路。文献[11-12]的方法是将刺激前记录的EEG按刺激间隔分段截取作为虚拟噪声通道,然后与单次含噪信号一起构成多路观测信号,从而建立如下模型
式中,A 为未知混合矩阵,s= [s1,v2,…,vM]T为独立源向量。虽然模型中假设多路虚拟噪声通道是独立源信号,但是ICA算法要求输出信号最大程度符合独立性准则,所以引入的虚拟通道之间的实际相关程度会束缚ICA算法的效果。事实上由于未知混合矩阵的元素给定值很多,所以解混矩阵的自由度会很低,这样ICA输出分量的独立性也会受到限制。
在本算法模型中,引入的是虚拟观测通道,也就是将刺激后记录的单导少次含噪信号形成多路观测信号,同时假设多路观测信号中的噪声均由有限个独立噪声模式组成,即
式中,uj,(j=2,3,…,M) 为独立噪声模式,aij为第j种噪声模式在第i路噪声中的权重。另外根据EP信号的锁时性[1],可假设单导相邻少次记录中 EP信号基本不变,这样虚拟通道模型可用矩阵表示为
式中,A 为未知混合矩阵,s= [s1,u2,u3,…,uM]T为独立源向量。在这个模型中未知混合矩阵假定元素很少,提高了解混矩阵的自由度,从而 ICA输出分量的独立程度得到了提高。因为模型中的每路观测信号之间(即每次含噪信号之间)的相关性不强,所以观测信号矩阵可以保证行满秩,也就是混合矩阵一定是满秩的。
2 实验与分析
2.1 ICA仿真实验
仿真实验选取4路独立源信号,分别代表 EP信号(通过叠加平均法得到的猫的体感诱发电位)、工频干扰(50 Hz正弦波)、自发脑电(高斯分布随机信号)和亚高斯噪声(周期方波)。4路源信号经过随机矩阵 A4×4混合得到4路观测信号,最后通过fastICA算法提取得到4路输出信号,源信号、观测信号和提取信号如图1所示。其中观测信号的混合矩阵 A4×4由 Matlab随机生成:
从图1中可以看到,源信号淹没在观测信号中,而提取后得到的信号与源信号非常相似,同时可以发现提取后得到的信号的顺序与源信号相比发生了变化,这是由 ICA的不确定性造成的。总体说来,对于符合ICA假设条件的源信号与混合方式得到的结果是比较理想的,即使在较低信噪比情况下提取结果与原始信号仍有很高的相关系数。图1中四路信号的相关系数分别为0.9994(EP)、-0.9999(正弦)、-0.9899(高斯)和0.9999(方波)。
2.2 虚拟通道ICA仿真实验
在虚拟通道ICA仿真实验中分别对两种模型进行了仿真。第1种是文献[11,12]的模型,其中第1路观测信号是EP信号与噪声的混合信号,信噪比为0 dB,另外3路为其他时间段记录的噪声,观测信号与提取信号如图2所示。第2种为所给出的模型,4路信号均为含噪信号,其中的EP信号是相同的,但噪声为不同时间段的记录,信噪比均为0 dB,观测信号与提取信号如图3所示。
从图2中可以看出4路输出信号中第一路为EP信号,但信噪比不但没有上升反而略有下降,这是因为4路观测信号中的噪声之间的相关程度很弱,彼此无法互相削弱抵消。其他3路输出信号与观测信号完全相同,这是由模型中混合矩阵的约束造成的,解混矩阵与混合矩阵有着完全相同的型式。而从图3中可以看到EP信号出现在输出信号中的第1路,信噪比得到了明显的提高。
2.3 视觉诱发电位的采集
图1 4路源信号、观测信号和提取信号。(a)EP信号;(b)正弦波;(c)高斯噪声;(d)周期方波;(e)~(h)相应的4路观测信号;(i)~(l)相应的4路提取信号Fig.1 Four channels of source signals,observed signals and extracted signals.(a)EP signal;(b)Sine wave;(c)Gaussian noise;(d)Cycle square wave;(e)~(h)Four channels of observed signals;(i)~(l)Four channels of extracted signals
图2 第1种虚拟通道模型的观测信号与提取信号。(a)~(d)4路观测信号;(e)~(h)相应的4路提取信号Fig.2 The observed and extracted signals with the first virtual-channel ICA model.(a) ~ (d)Four channels of observed signals;(e)~(h)Four channels of extracted signals
实验中使用的视觉诱发电位(visual evoked potential,VEP)的原始数据由美国 Neuroscan公司32导脑电设备Scan4.3采集得到。对12位受试者进行模式翻转视觉诱发电位(pattern reversal visual evoked potential,PRVEP)测试。在测试过程中进行多次刺激,经过机器自动识别与人工筛选去除了有明显伪迹的记录。使用传统的叠加平均法得到每导记录中的 VEP,图4为受试者 S1的 Cz,O1,Oz,O2四路导联的VEP波形。
从图4中可以看到 O1,Oz,O2等3个导联均在刺激后的100 ms左右出现了P100峰值,且每导记录中的VEP波形非常接近,而Cz导联没有出现明显的P100峰值,这反映了VEP主要出现在头皮枕骨处这一现象[1]。
2.4 多导联ICA的VEP提取
选择受试者 S1的 Cz,O1,Oz,O2等4个导联单次记录作为观测数据使用ICA算法进行处理,观测信号与提取信号如图5所示,由图5可见在4路提取信号中任何一路都没有明显出现期望的P100峰值。
图3 第2种虚拟通道模型的观测信号与提取信号。(a)~(d)4路观测信号;(e)~(h)相应的4路提取信号Fig.3 The observed and extracted signals with the second virtual-channel ICA model.(a) ~ (d)Four channels of observed signals;(e)~(h)Four channels of extracted signals
出现这种结果的原因在于,虽然多导脑电信号在形式上符合ICA的模型,但实际上并不真正符合。以两路观测信号为例进行分析,假设两路信号分别为xi=si+vi,i=1,2。对于相邻导联,其中的EP信号非常接近(见图4),因此观测信号可以简化为xi=s+vi(先不考虑EEG等噪声的相近程度)。
ICA算法的任何一路输出实际上就是多路输入的线性组合,即y=k1x1+k2x2,通过ICA算法可以得到在独立性准则条件下的加权系数k1和k2。然而无论k1和k2是通过什么准则确定的,y和s的最大相关系数是有限的。计算 y和s的相关系数 ry,然后通过对k1和k2求偏导可得最大相关系数为
式中,rv为v1和v2的相关系数,λ为信噪比。显然rv越小越好,而实验表明相邻导联中的EEG等噪声具有较大的正相关性(图5中 O1,Oz,O2这3个导联的波形比较接近),所以无论加权系数k1和k2如何选择,也就是ICA中的解混矩阵如何选择,都无法得到很好的结果。
2.5 虚拟通道ICA的VEP提取
图4 叠加平均法得到的VEP波形。(a)Cz导联;(b)O1导联;(c)Oz导联;(d)O2导联Fig.4 The VEP waveforms with ensemble average method.(a)The Cz channel;(b)The O1 channel;(c)The Oz channel;(d)The O2 channel
使用受试者S1的Oz导联的连续4次记录构建虚拟通道ICA模型,观测信号与提取信号如图6所示。图6中第3路提取信号中出现了P100峰值的轮廓,虽然与图4所示的叠加平均结果相比仍然有很大噪声,但原来完全淹没在噪声里的EP信号得到了提取,背景噪声得到了一定程度的去除。如果以图4所示的叠加平均结果作为纯净EP信号,经计算可知提取信号与观测信号相比信噪比提高了约为12 dB。此外实验也证明了随着虚拟通道数量(即刺激次数)的增加,EP信号提取效果也会更好。图7(b)为使用10路虚拟通道(连续10次记录)经ICA算法得到的结果,此时信噪比提高约为20 dB。
图5 使用Cz,O1,Oz,O2四导数据进行多导联 ICA的VEP提取。(a)Cz导联;(b)O1导联;(c)Oz导联;(d)O2导联;(e)~(h)相应的4路提取信号Fig.5 VEP extraction with the Cz,O1,Oz and O2 channels based on ICA method.(a)The Cz channel;(b)The O1 channel;(c)The Oz channel;(d)The O2 channel;(e)~(h)Four channels of extracted signals
图6 使用Oz导联连续4次数据进行虚拟通道ICA的VEP提取。(a)~(d)观测信号;(e)~(h)相应的提取信号Fig.6 VEP extraction with four times Oz channel observed signals based on virtual-channel ICA method.(a)-(d)The observed signals;(e)-(h)The extracted signals
图8为12位受试者分别在4路和10路虚拟通道ICA方法下得到的 Cz,O1,Oz,O2等 4个导联VEP相关系数的统计结果(均值与方差)。首先使用叠加平均法分别得到每位受试者 Cz,O1,Oz,O2各个导联的VEP作为标准信号。然后使用4路和10路虚拟通道ICA方法提取VEP,并计算与相应的标准信号的相关系数以及每类相关系数的均值与方差。通过统计分析可知在显著性水平α=0.01条件下4路和10路虚拟通道方法的 O1,Oz,O2导的VEP提取差别是显著的。此外,在实验中还观察到Cz导联中的VEP很微弱且成倒相现象,反映在统计结果中Cz导联的相关系数明显较小,这说明VEP并不主要出现在顶叶处。
图7 Oz导联VEP波形比较。(a)4路虚拟通道ICA的提取结果;(b)10路虚拟通道ICA的提取结果;(c)叠加平均法结果Fig.7 The extracted VEPs of Cz channel with different methods.(a)Four virtual channels ICA method;(b)Ten virtual channels ICA method;(c)The estimated result with ensemble average method.
图8 4路和10路虚拟通道 ICA方法下得到的 Cz,O1,Oz,O2四导联 VEP相关系数的统计结果。显著性差异由*表示(**P<0.01)Fig.8 The statistical results of estimated VEP correlation coefficients of the Cz,O1,Oz and O2 channels based on four and ten virtual channels ICA methods.Statistically significant differences denoted by asterisks(**P <0.01).
3 讨论和结论
ICA模型的解混矩阵相当于一个空间滤波器组,任何一路输出都是多路观测信号的加权组合。若多导观测信号中的EEG等噪声具有较强的正相关性,则任何一路输出都不可能减弱或消除噪声,所以很难得到期望的EP信号。如果能将观测信号在形式上组成ICA算法所需要的观测矩阵,而噪声之间具有较弱的相关性甚至负相关性,那么ICA是有可能达到目的的。这样空间信息可以来源于广义的信号空间而并不局限于物理空间(多传感器得到的多路观测信号)。因此,本研究针对单导EP信号提取问题,将单导多次观测信号以虚拟通道的观点构成多路观测信号,实现了EP信号的有效提取。仿真实验表明当源信号与混合方式符合ICA假设条件时可以得到很好的结果,即使在较低信噪比情况下提取结果与原始信号的相关系数仍然很高。但是在多导VEP提取中,并没得到预期的结果。这说明应用ICA对多导信号进行处理,其前提是多导观测信号必须符合ICA的假设条件。在多导脑电信号处理中,因为每导脑电信号都是EP与EEG的叠加,而ICA算法却要求它们都是来自某些独立源的某种线性组合,所以直接使用ICA方法处理多导脑电信号以实现信号分离是不适合的。而将单导观测信号以虚拟通道的观点构成多路观测信号以后,近似符合了ICA的应用条件,可以实现EP信号的有效提取。从空间滤波器角度分析可得出结论,基于虚拟通道的ICA算法实际上实现了独立性准则下的加权叠加平均效果。
[1]潘映辅.临床诱发电位学(第二版)[M].北京:人民卫生出版社,2000.
[2]Wang Yijun,Gao Xiaorong,Hong Bo,et al.Brain-computer interfaces based on visual evoked potentials [J]. IEEE Engineering in Medicine and Biology Magazine,2008,27(5):64-71.
[3]Davila CE,Mobin MS.Weighted averaging of evoked potentials[J].IEEE Trans Biomed Eng,1992,39:338-345.
[4]Bezerianos A,Laskaris N,Fotopoulos S,et al.Data dependent weighted averages for recording of evoked potential signals[J].Electroencephalography and Clinical Neurophysiology,1995,96:468-471.
[5]Aapo H,Erkki O.Independent component analysis:algorithms and applications[J].Neural Networks,2000,13(4-5):411-430.
[6]Zou Ling,Zhang Yingchun.Single-trial evoked potentials study by combining wavelet denoising and principal component analysis methods[J].Clin Neurophysiol,2010,27:17-24.
[7]Hu L,Zhang GZ,Hung YS,et al.Single-trial detection of somatosensory evoked potentials by probabilistic independent component analysis and wavelet filtering[J].Clin Neurophysiol,2011,122:1429-1439.
[8]Villa NC, James CJ. Independent component analysis for auditory evoked potentials and cochlear implant artifact estimation[J].IEEE Trans Biomed Eng,2011,58:348-354.
[9]师黎,钟丽辉,王端.基于虚拟通道ICA-WT大鼠视觉诱发电位少次提取[J].中国生物医学工程学报,2010,29(3):379-383.
[10]李洪,孙云莲.基于EMD虚拟通道的ICA算法在信号消噪中的应用[J].中国生物医学工程学报,2007,30(5):33-36.
[11]季忠,金涛,杨炯明,等.虚拟噪声通道在基于ICA消噪过程中的应用 [J]. 中国机械工程,2005,16(4):350-364.
[12]刘钦团,邱飞岳,李浩君.基于虚拟通道的 ICA的 P-VEP提取方法的研究[J].计算机工程与科学,2010,32(7):137-139.
[13]周宗潭,董国华,徐昕,等.独立分量分析[M].北京:电子工业出版社,2007.