基于参数模型和分量分析的事件相关诱发电位单次提取

事件相关诱发电位单次提取

王荣昌　都思丹∃

(南京大学电子科学与工程系,南京　210093)

　　摘要　我们针对脑电事件相关电位(ER P)这种信噪比极低的信号检测问题,提出了两种ER P信号单次提取方法,能非常有效地同时去除自发脑电、眼动伪迹和工频噪声三种常见噪声。(1)首次对自发脑电、眼动伪迹和工频噪声这三种常见成分连同事件相关电位同时进行A RX建模,利用基于最小二乘(L S)的A RX算法进行参数辨识获得提取结果;(2)利用分量分析,采用Fast I CA算法进行事件相关电位的提取。明确指出I CA分解的一些重要分解特性及其内在机理,针对实际情况对Fast I CA算法进行了改进,实现了分解结果对ER P成分的自适应映射。数值仿真实验结果表明两种方法均有较高的信号分解提取能力。

关键词　事件相关电位　A RX模型　分量分析　Fast I CA算法

Extraction of Si ngle-tr i a l Even t-rela ted Poten ti a ls by M ean s of ARX M odel i ng and I ndependen t Com ponen t Ana lysis

W ang Rongchang　D u Sidan

(D ep art m ent of E lectronic S cience&E ng ineering,N anj ing U niversity,N anj ing　210093,Ch ina)

　　Abstract　T he p resen t paper focu sed on the ex tracti on of even t2related po ten tials on a single s w eep under ex trem ely low S N rati o.Tw o m ethods that can efficien tly remove spon taneou s EEG,ocu lar artifacts and pow er line in terference w ere p resen ted based on A RX modeling and independen t componen t analysis(I CA).T he fo rm er m ethod app lied A RX model to the m easu red compound signal that ex ten sively con tained the th ree k inds of o rdinary no ises m en ti oned above,and u sed A RX algo rithm fo r param etric iden tificati on.T he latter decompo sed the signal by m ean s of independen t componen t analysis.Besides,som e of I CA’s i m po rtan t decompo sing characters and its in trin sic cau sality w ere po in ted ou t defin itely.A cco rding to the p ractical situati on,som e modificati on on Fast I CA algo rithm w as also given,so as to i m p lem en t au to2adap tive m app ing of decompo sed resu lts to ER P componen t. T h rough si m u lati on,bo th the tw o w ays are p roved to be h igh ly capab le of signal ex tracti on and S N rati o i m p roving.

Key words　Even t2related po ten tials(ER P)　　A RX model　　Independen t componen t analysis(I CA)　　Fast I CA algo rithm

1　引　言

脑电中的事件相关电位(Even t2related po ten t2 ial,ER P)是研究人类认知心理的相当重要的一类微弱脑电信号,是现代研究人的认知功能最有效的途径之一。然而,从头皮检测得到的ER P信号往往只有几个微伏[6],且往往淹没在几十个微伏的自发脑电(E lectroencep halogram,EEG)背景噪声中。另外,不可避免地会混入幅度相当大的眼动伪迹,以及传统的50H z工频干扰。我们就如何从强大的随机背景噪声中提取出微弱的ER P阵发信号的新方法进

:128@..行了探索研究。

传统的提取ER P的方法被称作“迭加平均”,该方法建立在有用信号与噪声的统计以及噪声的零均值假设之上。然而由于眼动伪迹产生的不确定性以及50H z工频干扰的相对稳定性使得它们无法通过“迭加平均”的方法消除。另外,当涉及到记忆、认知、学习等应用领域时,每一次的诱发电位都有相应的改变和波动,而“迭加平均”却无法保留单次诱发信号中包含的信息。

近年来,各种现代信号处理的方法逐渐被应用于生物医学信号领域。其中参数模型辨识方法被运用到视觉诱发电位的提取[5]。20世纪90年代后期发展起来的分量分析()作为一种新兴的盲源

生物医学工程学杂志

J B i om ed Eng　2006;23(6)∶1222～1227　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

分离(B lind sou rce sep arati on ,B SS )技术,能够处理一组由统计信源线性组合得到的混合信号,最终得到各信源分量。I CA 分析在EEG 和脑磁信号(M egnetoencep halograp h ,M EG )的伪迹去除方面得到应用[3]

。我们也将I CA 应用于脑电中事件相关电位ER P 信号的提取方面的研究。

基于单次检测到的脑电信号可看作事件相关电位ER P 和三个“噪声”分量(包括自发脑电EEG ,眼动伪迹R EO G 和50H z 工频噪声)

之迭加的思想。从统计信号处理的角度来看,这四个成分可认为相互统计。基于上述普遍假设,分别利用A RX 模型参数辨识和分量分析(I CA )对ER P 信号进行单次提取。

2　理论与方法

2.1　ARX 参数模型

A RX 模型是基于自回归(A u to regressive ),引

入外来输入(exogenou s inp u t )的一种复杂模型。利

用A RX 模型对本问题建模基于这样一个普遍假设:检测所得的混合信号y (k )由几个加性成分构成,如图1所示。结合实测脑电信号可能存在的各种伪迹和干扰,本文特别引入了50H z 的工频干扰项,使得讨论更符合实际情况。

图1　脑电混合信号的ARX 模型

F ig 1　ARX model of measured co mpound signal

　　s (k ):事件相关电位ER P ,受靶刺激产生,波形持续时间约1s [6],具有一定形状的阵发信号。可看作一模板参考信号(由事先迭加平均所得)经过滤波的所得[4];n (k ):自发脑电背景噪声,EEG 。在短时间内可看作是白噪声经过一个A R 模型的输出[4,8];v (k ):眼动伪迹,R EO G ;p (k ):50H z 工频噪声。

以上的四个信号成分在频域上互相混叠,因此无法用带通滤波的手段进行分离。然而,注意到,这四种信号分别由不同的彼此无关信号源产生,因此可以看作具有统计性,可以利用信号的统计特性来分解。

　　模型的数学表达式:

y (k )=s (k )+n (k )+v (k )+p (k )

(1)

y (k )=-2p

i =1a i y (k -i )+

2q 1+d -1j =d

b j u (k -j )+

2q 2+∆-1

Σ=∆

c Σx (k -Σ)+

2q 3+Ρ-1

Κ=Ρ

d Κl (k -Κ)+

e (k )(2)

式中:p ,q 1,q 2和q 3分别为自回归(au to regressive )

和滑动平均(m oving average )模型系数;u (k )为ER P 参考信号,可通过先验的迭加平均信号获得;x (k )为眼动参考信号,从眼部电极直接得到;l (k )为

工频参考信号,直接用50H z 余弦信号表示;d ,∆和Ρ分别为单次信号相对于参考信号的时间延迟。　　转换到Z 变换域上,得到:

Y (z )=B 1(z )A (z )U (z )+B 2(z )

A (z )X (z )+

B 3(z )A (z )L (z )+1

A (z )E (z )

(3)

其中:

A (z )=1+2p

i =1

a i z -i

;B 1(z )=

2

q 1+d -1j =d

b j z

-j

;B 2(z )=

2q 2+∆-1

Σ=∆

c Σz -Σ

;B 3(z )=2q 3+∆-1

Κ=∆

d Κz -Κ。s (k )可看作一模板参考信号u (k )(由事先迭加平均所得)经过传递函数B 1(z )

A (z

)滤波的所得[4];同理,眼动

的伪迹v (k )由直接从眼电极测得的信号x (k )经过

传递函数B 2(z )

A (z

)滤波的所得;工频噪声p (k )由50H z

参考信号经过传递函数B 3(z )

A (z

)滤波表示。

根据A kaike 信息准则[9],在滤波效果(要求阶数高)和计算复杂度(阶数越低越好)之间取折中后,

选定滤波阶数p =12,q 1=q 2=q 3=8。

对比前面提到得数学表达式(1),可以直接检测得到的数据如下:一个模型输出y (k ),三个模型输入u (k )、x (k )和l (k )。于是,可以利用基于L S ,使估计均方误差最小化的A RX 算法来估计A RX 模型中的各模型参数,进而计算得到单次ER P 信号。2.2　分量分析

分量分析(Indep enden t com ponen t analy 2sis ,I CA )是近年来发展起来的一种新的盲源分离(B lind sou rce sep arati on ,B SS )方法,能够完全分离

多变量混合数据各成分。其分析处理的对象往往是非高斯信号,因此更具普遍适应性。

I CA 分量分析的使用有几个前提:(1)各源分量必须满足统计性;

(2)观测信号的个数必须大于等于待分解的独

3

221第6期　　　　　　　　　　　　王荣昌等。　基于参数模型和分量分析的事件相关诱发电位单次提取

假设有一组N维观测信号,记作x j(t),j=1,2,…,N;t=1,2,…,T。它由若干变量s i(t),i=1, 2,…,M经过加权混合产生。信号处理的目标是从这N个观测信号中分离出各个干净的s i(t)。

其中,任意x j(t)均可表示为:

x j(t)=a j1s1+a j2s2+…+a jM s M(4)　　表示成矩阵的形式:

x1(t)

　…x N(t)=

a11 (1)

　…　…　…

a N1…a N M

s1(t)

　…

s M(t)

(5)

　　简写为:

x=A・s(6)其中:A为由系数a ji构成的N×M维矩阵。上式中,唯一可知的是观测向量x。这就是I CA需要解决的盲源分离的问题:要估计解混矩阵W,使得当x通过W时的输出为y,并且满足y各个分量之间尽可能相互。

分解模型的数学表达用矩阵表示为:

y1(t)

　…y M(t)=

w11 (1)

　…　…　…

w M1…w M N

x1(t)

　…

x N(t)

(7)

y=W・x(8)

　　至此,多变量分析问题已经简化为对系数矩阵的求解问题。本文在求解此问题时采用了Fast I CA 算法[1]。

Fast I CA算法基于非高斯性(N ongau ssian ity)极大准则,而非高斯性的存在恰恰是分量分析方法必需的前提条件。中心极限定理表明:一组相互的随机变量线性组合后得到的混合变量必接近于高斯分布。换句话说,源信号的非高斯性比混合信号非高斯性强。于是,可以对分离结果的非高斯性(N ongau ssian ity)进行度量,当其达到最大时,可认为实现了最佳分离。大部分的I CA算法的设计准则也正是基于此。

Fast I CA算法的原理可描述为:使用某种目标函数(负信息熵,N egen trop y)来标度输出分量的非高斯性,然后找到某种迭代算法来使得目标函数最大化(M ax i m izati on)。

信息量公式H(y)=-f(y)logf(y)dy,表明“信号愈不确定(随机),包含信息量愈大”[3]。有结论表明:在方差相等情况下,所有随机变量中,高斯随机变量具有最大的熵。因此,可以用信息熵来度量非高斯性。为了理解和计算方便,这里引入“负”的概念,定义“负信息熵”,使其具有非负定性:高斯性的时候等于零,非高斯性的时候大于零;非高斯性越强,负信息熵越大。

　　负信息熵定义为:

J(y)=H(y gauss)-H(y)(9)其中:y gauss是和y具有相同的协方差矩阵的高斯随机信号。

负信息熵从统计意义上讲是非高斯性的最优估计,但是计算上比较困难。直接用定义计算需要估计信号的概率密度函数,这显然是不实际的。因此,通常利用式(10)估计负信息熵,即

J(y)∝{E[G(y)]-E[G(v)]}2(10)式中:y为零均值,单位方差的待估计信号;v为零均值,单位方差的高斯随机信号;G为任意非二次函数。一般选用上升比较慢的函数,有两个典型:

G1(u)=

1

a1

logco s(ha1u),G2(u)=-exp(-u2 2)。其中:a1取1～2之间的合适定值,一般取a1=1。

Fast I CA算法步骤:

(1)对观测信号去均值;

(2)利用EV D特征值分解法预白化观测信号,得到x

。令

k=0;

(3)令n=0,k=k+1;

(4)初始化解混向量w k;

(5)n=n+1;

(6)调整w k:g(u)=uexp(-u2 2),w k(n)= E{x g(w-(n-1)T x)}-E{g′(w(n-1)T x)}w(n-1);

(7)归一化w k:w k(n)=w k(n) ‖w k(n)‖;

(8)算法收敛否,如果未收敛则转到步骤(5),重新调整w k;如果收敛,执行下一步;

(9)分解矩阵所有行向量w k是否求解完毕(k =N)。如果未完毕,则转到步骤(3)继续分解下个行向量,如果求解完毕,则执行下一步;

(10)求分量:x=W x

ζ。

　　将分量分析应用于事件相关电位ER P的提取,有着物理意义明确,简洁、快速、高效的优点。其分析所用数据仅需混合信号(本文中即多导检测信号),更适合实时处理,同时也为脑机接口(BC I)方面的应用提供了思路。它不像前面提到的

4221　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　　　第23卷

参数模型方法,无需ER P 信号的事先迭加平均信息,无需眼部参考电极作为眼动伪迹的参考信号,因此也具有更强的适应性。

3　数值仿真

为了验证两种方法提取ER P 信号的效果,本文中用一组数值仿真实验来模拟真实检测综合脑电信号的ER P 提取。用于仿真的信号选取如下:s (k )(ER P )和v (k )(R EO G ),模拟各自真实信号的相应特征[4,6];n (k )(EEG )为仿真的12阶A R 模型输

出[4];p (k )(工频噪声):任意相位50H z 余弦信号。此四种信号的幅度按照实际信号的信噪比关系,调整到-30db 左右(按通常的情况考虑)[4]。由以上四种信号成分加性混合构造一仿真检测综合信号y (k ),如图2。图

3左显示了y (k )各成分。

图2　单次扫描仿真信号

F ig 2　Si m ulated

measured signal of si ngle -tr i al

图3　ARX 参数模型辨识提取ERP 结果

F ig 3　Effects of ERP extraction by mean s of ARX modeli ng

3.1　ARX 参数模型辨识

模型的参数估计思想如下:首先根据模型输出y (k ),模型输入u (k )、x (k )和l (k ),利用L S (L east Square )算法估计出模型各部分参数。包括自回归A (z )和滑动平均B 1(z ),B 2(z ),B 3(z )模型系数,这里

分别用a ,b 1,b 2,b 3表示之。然后,利用辨识得到的模型参数(即滤波系数),对参考信号进行滤波,得到结果。

图3显示了A RX 参数模型辨识分解脑电检测信号的效果。可以看到分解滤波之后得到的各信号成分与原仿真信号的成分非常相似:眼动信号的影响几乎被完全地再现出来;另外,有用信号ER P 从波形和时域上看,其基本信息也被保留了下来。

辨识的结果也很好地体现出该方法的一致性,即仿真得到的A RX 模型系数(A (

z )系数)只是轻微地与仿真所用的原系数有所不同,基本保持一致,如图4所示。

图4　ARX 模型参数辨识结果一致性比较

F ig 4　Con sistency of parameter iden tif ication by mean s of ARX

modeli ng

3.2　分量分解

3.2.1　利用Fast I CA 算法分解　同前所述,根据实

际从头皮测得的混合脑电信号的特征[4,6],认为其包含四个相互统计的分量:①事件相关电位ER P ;

②自发脑电EEG ;③眼动伪迹R EO G ;④50H z 工频噪声。式(1)在此处仍然成立,重新表述为:

　y (t )=s (t )+n (t )+v (t )+A sin (2Πf 0t +Η)(11)　　注意到这里工频噪声出现了一个随机相位的问题。与前面参数模型估计有所不同的是:I CA 分解

中工频信号的相位无法通过类似滤波的方法进行估计。为此将式(11)改写为:

y (t )=s (t )+n (t )+v (t )+asin (2Πf 0t )+bcos (2Πf 0t )

(12)

　　这样,就可以引入两路“工频参考信号”,p c (t )=co s (2Πf 0t )和p s (t )=sin (2Πf 0t )。这样,有理由认为观测信号是由目标信号源s 0(t ),n 0(t ),v 0(t ),sin (2Πf 0t )和co s (2Πf 0t )线性加权混合而成,因此其相应的I CA 线性混合模型可描述如下:

　　y j (t )=a j 1s 0(t )+a j 2n 0(t )+a j 3v 0(t )+

5

221第6期　　　　　　　　　　　　王荣昌等。　基于参数模型和分量分析的事件相关诱发电位单次提取

　　　　　　　　　　　　　j=1,2,3(13) r s(t)=asin(2Πf0t)=a44p s(t)

r c(t)=acos(2Πf0t)=a55p c(t)

(14)写成矩阵的形式:

y1(t)

y2(t)

y3(t) r s(t) r c(t)=

a11a12a13a14a15

a21a22a23a24a25

a31a32a33a34a35

000a440

0000a55

s0(t)

n0(t)

v0(t)

p s(t)

p c(t)

(15)

　　利用前面介绍的Fast I CA方法,求取解混矩阵W即可由下式得到各分量。

s0(t)

n0(t)

v0(t) p s(t) p c(t)=

w11w12w13w14w15

w21w22w23w24w25

000w440

0000w55

y1(t)

y2(t)

y3(t)

r s(t)

r c(t)

(16)

　　根据I CA前提②以及对ER P信号提取问题的抽象可知:必须获取

3导或3导以上(本例中设定3导信号)的同步头皮检测波形,同时提供两个0相位的正余弦信号。事实上,只需要从多个电极上同步记录检测波形即可。这里分别对A RX模型参数法仿真中用到的4种仿真成分通过随机加权的形式得到3个仿真的检测信号,以模拟从各不同电极得到的各自不同的信号。同时还必须确保各成分的比例关系符合真实脑电信号的特征。可参见图2,3导仿真检测信号的波形大致同此,且这3个检测信号互不相同。

I CA分解的结果如图5所示,从上至下依次为s0 (t),v0(t),p c(t),n0(t)和p s(t)。

3.2.2　分解向量的自适应映射　分量分离方法为了计算简便,归一化各分量为具有单位方差[1,2]。因此,这里得到的仅仅是假设的具有单位方差的分量,还需要进一步映射各分量到实际的混合信号成分。这就涉及到I CA分解的一些特别的性质:①无法确知分量的方差(能量),因为处理时均归一化为具有单位方差;②分解信号有时候会相差一个“负号”;③各分量分解结果的顺序随机不定。

　　前两点都可以解释为:各分量能量比例以及信号的正负都可以在计算混合矩阵时加入相应比例系数和正负号予以抵消;后一点可解释为:A中各列向量位置作相应的左右调整。

基于此,在分析了ER P分量间先验知识的前提下,提出了对Fast I CA算法的一种改进(事实上,如果研究时将ER P信号都建立在归一化为单位方差的情况之下,可以直接以这里得到的结果作为ER P 的提取信号以供后期分析之用)。

图5　FastI CA算法分解分量结果

从上至下依次为s0(t),v0(t),p c(t),n0(t)和p s(t)

F ig5　D eco mposition by mean s of FastI CA

from top to bo ttom:s0(t),v0(t),p c(t),n0(t)&p s(t)

　　分析脑电混合信号的各成分以及混合方程式(15)可知:①ER P信号的能量在各成分中最小,对应A中的元素应相对较小。②工频参考信号对应的A中的行向量其他元素理论值均应为0,参见式(15)。因此,根据上面两条准则,对混合矩阵A中的列向量进行恰当选取,就可以得到ER P信号对应的3个“比例”系数(分别对应着3导观测信号中ER P 成分的大小)。在处理过程中又将这三个ER P信号取平均,以进一步减小可能发生的某一个观测误差影响。采用这种方法,使得分解结果具有一定的自适应映射的能力:可随着每一次分解向量s的分量次序、正负,自适应地计算获得实际的ER P成分。如图6显示了映射结果。

图6中3条虚线分别是3个仿真检测信号中的ER P仿真成分,粗实线是提取结果。可以看到ER P 信号提取结果令人满意。其波形起伏以及诱发电位产生的时间信息基本被提取出来。信噪比得到了很大的改善。这将十分有利于后续对ER P成分的分析工作。

最后还有值得一提的是:I CA分量分解方法同时也可很方便地单独将脑电信号中的工频噪声成分剔除。传统的陷波器带宽不可能做的无限小,因

6221　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　　　第23卷此在使用时不免会损失有用信号的一部分重要细节,而I CA分解则完全不存在类似的问题。

图7显示了工频滤除的效果。图中的两条曲线显示的分别是I CA分解滤除前后的信号的功率谱密度(B u rg法[7,11])。很明显,原检测信号功率谱在50H z处有一非常明显的尖峰,说明原检测信号包含了比较大能量的50H z工频成分。滤波之后的信号功率谱曲线如实线所示,可见滤除效果较好,且能比较好地保留有用信号的信息成分

。

图6　事件相关电位ERP信号提取结果

F ig6　Effects of ERP extraction by mean s of FastI

CA

图7

　50Hz工频噪声滤除效果

F ig7　Results of Power-li ne i n terference re moval

4　结　论

事件相关电位ER P的动态变化情况对研究人

的认知过程具有重要意义,因此开发一种能够实现

单次提取的方法成为迫切需要。A RX参数模型辨识

和I CA分量分析这两种不同的方法进行事件相

关电位ER P信号提取,都取得了比较好的效果:能

够极大地保留ER P信号的有用成分(包括波形和出

现时间),同时可干净去除工频干扰以及眼动伪迹。

两种方法相比较,I CA分量分析作为一种与传

统信号分析手段较为不同的方法有一定优势:计算

速度快,所需的信号先验知识较少,无需提供ER P、

眼动参考信号,信号检测步骤简单等等。值得一提的

是,我们还对I CA分析方法做出了针对ER P信号提

取这个特定问题的一点改进,使得能够自动从分解

向量中获得所需要ER P信号成分,取得了较好的效

果。I CA分析将成为生物信号处理领域的一种有效

的分析工具。

参考文献

1H yvarinen A,O ja E.Independent component analysis:

algo rithm s and app licati ons.N eural N etwo rk s,2000;13(425)

∶411

2H yvarinen A.Fast and robust fixed2po int algo rithm s fo r

independent component analysis.IEEE T rans O n N eural

N etwo rk s,1999;10(3)∶626

3H yvarinen A,Karhunen J,O ja E.Independent component

analysis.N ew Yo rk:John W iley&Sons,Inc.2001:86288,

1052387,4072416

4Cerutti S,Ch iarenza G,L iberati D,et al.A param etric

m ethod of identificati on of single2trail event2related po tentials

in the brain.IEEE T rans B i om ed Eng,1988;35(9)∶701

5Cerutti S,Baselli G,L iberati D,et al.Single s w eep analysis

of visual evoked po tentials th rough a model of param etric

identificati on.B i o l Cybern,1987;56(223)∶111

6Pan YF.Evoked po tentials in m edical p ractice.2nd ed.

Beijing:Peop le’s M edical Publish ing House,2000:5792622

[潘映辐.临床诱发电位学.第二版.北京:人民卫生出版社,

2000:5792622]

7M ano lak is D G,Ingle V K,Kogon S M.Statistical and adap tive

signal P rocessing.Beijing:Publish ing House of E lectronics

Industry,2000:3762456,6562667[M ano lak is D G,Ingle V K,

Kogon S M.统计与自适应信号处理.北京:电子工业出版社,

2000:3762456,6562667]

8Sp reckelsen M V,B romm B.E sti m ati on of single2evoked

cerebral po tentials by m eans of param etric modeling and

kal m an filtering.IE E E T rans B io m ed E ng,1988;35(9)∶691

9A kaike H.Statistical p redicto r identificati on.A nn Inst S tatist

M ath,1970;22∶203

10Zhang DX.EEG signal p rocessing th rough w avelet and

independent component analysis.(D issertati on fo r M aster

D egree D ep t EE).A nhui U niversity,2002:15236[张道信.基

于小波和分量分析的脑电信号处理(硕士毕业论文).安

徽大学,2002:15236]

11W u XP,L i XH,Kong M,et al.H armonic esti m ati on and

removal based on independent component analysis.

T ransacti ons of Ch ina E lectro technical Society,2003;18(4)∶

56[吴小培,李晓辉,孔　敏等.基于分量分析的谐波估

计和消除.电子技术学报,2003;18(4)∶56]

12Zhang XD.M odern signal p rocessing.2nd ed.Beijing:

T singhua U niversity P ress,2003:692112,2152304[张贤达.

现代信号处理.第二版.北京:清华大学出版社,2003:692112,

2152304]

(收稿:2004204226　　修回:2004212201)

7221

第6期　　　　　　　　　　　　王荣昌等。　基于参数模型和分量分析的事件相关诱发电位单次提取