脑功能成像及在人文社会科学中的应用
上QQ阅读APP看书,第一时间看更新

3.2 脑电的信息反演技术

3.2.1 高分辨脑电唯一性定理

近年来,高分辨脑电成为脑电逆问题研究中的一个热门领域。它主要研究怎样从已知头皮表面的电位去推算皮层表面的电位分布。皮层表面的电位活动在通过颅骨区时,颅骨的低电导率会对皮层的脑电信号产生空间低通滤波的作用,从而导致头皮表面的电位信号变得模糊不清。因此,若能通过头皮表面观测到的微弱信号还原出皮层表面的电位分布,则可以获得高分辨率的脑电信息。正因为如此,一般称研究解决这种脑电计算的课题为高分辨脑电问题研究。

通常情况下,逆问题计算的解答并不唯一。这是由于辨识分析所获得的信息不足。那么高分辨脑电的计算是否存在唯一解呢?如果存在唯一解,就能为验证高分辨脑电算法的正确性提供有力的理论依据。

3.2.1.1 定理的表述

高分辨脑电唯一性定理的表述如下。

设头皮和皮层之间的区域为无源静态场,区域中的媒质是线性媒质,且没有电流从头皮流入空气中,如图3.5所示。当头皮表面的电位分布给定时,则皮层表面的电位分布是唯一确定的。

图3.5 头模型示意

3.2.1.2 相关的命题

无源场域Ω中的电位函数φ满足拉普拉斯方程 ▽2φ= 0。根据电磁场理论,若给定了场域的全部边界面上的电位或者场强法向分量,则拉普拉斯方程的解答是唯一确定的。全部边界面也可以被分割开,它的一部分给定电位值而另一部分给定场强法向分量,这样的解答仍然是唯一的。高分辨脑电问题中,无源场域Ω的全部边界是由头皮表面Ss和皮层表面Sc合成的,但给定的仅是头皮Ss处的边界值,缺少皮层Sc处的边界值,这是非全部边界条件的唯一性问题。然而,头皮Ss处还满足场强法向分量为零,由此提出以下两个命题。

命题1:若给定整个头皮表面Ss处电位和零场强法向分量这种双重边界值,拉普拉斯方程的解答唯一性定理是否成立?

命题2:实际头皮表面Ss处能测到的电位仅是局部。若给定局部头皮电位和整个头皮表面Ss零场强法向分量这种边界值,求解皮层表面Sc电位的解答是否唯一?

两个命题的主要区别是:命题1需要知道头皮表面(这里假设头皮表面是闭合面)所有点的电位,而命题2只需知道头皮表面局部的电位值。

Yamashita(1982)给出了心电问题中唯一性定理的证明。我们试图据此来证明高分辨脑电的唯一性,以证明上述两个命题。结果是前者可以证明,但后者难以证明。问题出在什么地方呢?

Yamashita(1982)证明唯一性定理时利用的是以下公式讨论命题2时,我们发现了若干这一命题的反例,现举一例说明如下。

式中:SHSL相当于脑电模型的皮层表面Sc; φ是所构造的格林函数;表示ψ沿面S法向的偏导数。利用这个公式以及皮层表面形状的不规则性,推出只有当皮层表面Sc处处电位为0时才有式(3.1)成立,因而证明皮层表面Sc处处电位为0,也就是皮层表面的电位唯一。但实际上,式(3.1)只能说明在皮层表面积分为0,即φ函数正交,而φ本身在皮层表面不一定处处为0。一般而言,式(3.1)的适用范围并不明确。因此,当采用这种证明方法

首先,将问题简化为二维的情况,并且假定模型是规则的几何形状。

图3.6 命题2的反例

如图3.6所示的单位正方形。假设在x=0的边上,在y∈[0,1]区间的电位和电位沿x正方向的偏导数满足

φ= 0,对任意的点{px, y)|x= 0, y∈[0,0.5]};

,对任意的点{px, y)|x= 0, y∈[0,1]};

φy+1)=φy),即周期性边界条件。

那么,根据命题2的条件可以推出电位φ的解答在整个正方形区域均为0。但我们找出了完全满足该边界条件的非零解,即

式中:U0为一任意常数。根据下列边界条件,利用无源场所满足的拉普拉斯方程可以求解该式。

φ= 0, ,对任意的点{px, y)|x= 0, y∈[0,0.5]};

φ=U0, ,对任意的点{px, y)|x= 0, y∈[0.5,1.0]};

③周期性边界条件:φy+1)=φy)。显然,它也满足只给定x=0, y∈[0,0.5]的无源场边值问题,因为

x= 0时,有

从而得到命题2的反例。

上面的反例一方面说明命题2不成立;另一方面,采用Yamashita(1982)的方法却能够“证明”命题2是成立的,这进一步说明该证明方法存在缺陷。因此,需要尝试通过新的方法来证明高分辨脑电唯一性定理(命题1)。

3.2.1.3 定理的证明

采用反证法证明命题1。设该命题的拉普拉斯方程有两组不同的解φ′φ″,令差场为φ=φ′-φ″,其边界值为

式中:下标s代表头皮(scalp), c代表皮层(cortical)。因为头皮上的电位已知,所以φSs= 0。故Ss为等位面,Jt Ss= 0, Jt代表切向电流密度;又因为Jn Ss=0, Jn为法向电流密度,因此J Ss=0,即头皮表面电流密度为0。因为所考察区域内为无源场,所以▽·J= 0,或▽2φ= 0,这里的φ, J均指差场中的变量。若能证明差场为零,即可证明唯一性。

由于生物体可视为准静态场,▽·E= 0,所以电场分布可用电位和等位面来描述。首先证明在所考察的场域Ω内部(不包括边界),不存在电位的极值点。以极大值点为例,极小值点的证明依此类推。

设点pΩ,且pSs, pSc。若点p为电位的极值点,例如极大值点,则可以包围p点作一个小的闭合面Sp,使Sp所限定的域含于Ω,且φp大于Sp面上任何点的电位。因此,面Sp上电力线的方向都是从面Sp发散出去的。对于线性媒质,J=σE,电流线和电力线的方向一致。所以,>0(本章均以外法线方向为正方向),而SpΩ, Sp包围的区域Ωp中处处有▽·J= 0。根据高斯通量定理,有

两个结果自相矛盾。同理,若p点电位为极小值,则,结果也相互矛盾。因此,在场域Ω内部不可能存在电位的极值点,p点只能处在边界上。(注:以上证明场域内部不存在电位的极值点也可以直接利用调和函数的最大值原理,即:如果函数φ在闭区域Ω+S是有定义的,而且连续,又假定φΩ内部满足方程式▽2φ=0,则函数φ必然在曲面S上达到最大值或最小值。)

现将场域分层剖分。定义参数,它表示Ss面上任一点到面Sc距离的最大值。在场域Ω中,包围面Sc作一闭合面S11,如图3.7所示,有S11Ω。令S10 =Ss, S12 =Sc,以及,使得D11=D12。这样,S11Ω分为Ω11Ω12两个场域。

图3.7 头模型剖分闭合面示意

下面继续递推。在Ω11Ω12中,包围S11S12分别作闭合面S21S23,有S21Ω11, S23Ω12, S21S23Ω11Ω12分别分为两个场域Ω21Ω22, Ω23Ω24,且D21=D22=D23=D24,这里令S22=S11, S24=S12

依此类推,设第n次剖分将Ω分为2 n个场域,分别为Ωn 1, Ωn 2, …, Ωn 2n,各个场域中Dn1 =Dn2 = …=Dn2n,各个闭合面分别设为Sn0, Sn1, Sn2, …, Sn2n,其中Sn0 =Ss, Sn2n =Sc

下面考察头皮上任一点p。对差场电位而言,头皮表面为等位面。在p点作切平面并取为xpy坐标平面,其外法线方向为z轴正方向,在ypz平面上头皮表面为apb曲线,如图3.8所示。在p点处,apb等位曲线的曲率半径为ry。同理,可定义在p点处,xpz平面上的相应曲率半径为rx。又因为头皮表面为等位面,p点邻域的拉普拉斯方程可以表示为

, R是表面主曲率半径的平均值,则式(3.8)可化为

解此方程,令

可得

图3.8 头皮表面附近的电位

由于z= 0点的Jz= 0,且φ= 0,所以C1, C2均为0。两个条件Jz= 0和φ= 0缺一不可。也就是说,在点p的充分小的邻域内满足Jz= 0和φ=0。因为p是头皮上的任意一点,按照上述分层剖分方法,由积分中值定理,式(3.11)中指数型函数的特点是不过零,由头皮表面的Jz= 0,可以推出系数C1为0,因此总可以找到n,使得在Ωn1域内,包括在最大值或最小值处,处处满足Jz= 0和φ= 0。因此,在Sn1面上有Jz= 0和φ= 0,表明Sn1为等位面。同理,可以推出在Sn2面上有Jz= 0和φ= 0, Sn2为等位面。依此类推,得到Sn2n面上有Jz= 0和φ= 0,即皮层表面的电位为0。命题1的唯一性得证。

上面的证明只能得出命题1的唯一性,却得不到命题2的唯一性。因为如果只知道头皮表面上的一部分电位,而头皮表面上其他区域的电位却是未知的,就无法形成一个闭合的层面,无法将邻域的性质推广为区域的性质,进而一层一层往里推。因此,上述证明方法推不出命题2。由此可见,定理的适用范围是已知头皮表面上所有点的电位,才能唯一确定皮层表面上各个点的电位。

3.2.1.4 小结

当已知整个头模型表面的电位分布,求皮层表面的电位分布时,根据唯一性定理,可以证明所得的解答是唯一确定的。也就是说,当采用不同的方法求解同一高分辨脑电问题时,解答必然是相同的。因此,高分辨脑电唯一性定理为算法的相互验证和实际的临床应用奠定了基础。本章进一步证明,由头皮表面上局部的电位分布求解皮层表面的电位分布,并不具有唯一性。

3.2.2 高分辨脑电

人脑电生理活动的对外表现形式之一是,在头皮表面记录到的头皮电位信号,即脑电图(electroencephalogram, EEG)。它已经在临床诊断与监护、病理和生理研究中获得了广泛的应用。脑电图具有的突出优点是其记录信号的时间分辨率很高,可以容易地做到毫秒级的水平,这是其他脑功能性成像手段,如正电子发射断层成像(positron emission tomography, PET)、功能性磁共振成像(fMRI)等目前无法达到的。但是,传统脑电图技术存在的严重缺陷是信号的空间分辨率太低,远不能满足现代神经科学研究对脑内电生理源进行定位的要求。造成脑电图空间分辨率不高的原因一方面是传统脑电图使用的电极数量太少(一般为16~20导),即信号的空间采样率太低,这可以通过增加电极数来克服;另一方面是颅骨的导电性能差,削弱了头皮电位的强度并使头皮电位的分布变得模糊(blurring),效果相当于对电位分布施加某种空间低通滤波,而这种空间模糊效应不可能通过增加电极来改善。据推算,颅骨低通滤波器的点扩散函数约为直径2.5cm的圆斑,所以电极数目增加到128~256导时就足够了,此时电极间的距离为1.95~2.25cm,已经接近或超过空间频率的Nyquist要求。

提高脑电图空间分辨率的根本途径是对脑内电活动源进行空间定位,这就是所谓的“脑电逆问题”。通常的方法是,把头皮电位等效地看成由皮层下一个或几个电流偶极子所产生,并通过非线性优化来确定偶极子位置,称为“等效偶极子定位”。但这类方法只能解决少数几个偶极子的定位问题,因此不适用于弥散型的电活动源定位。脑电逆问题研究近期的主要发展是把对孤立偶极子的定位扩展到对颅内电位、电流或偶极矩分布密度的计算,“高分辨脑电图”(high resolution EEG, HREEG)技术正是其中的主要形式之一。

3.2.2.1 高分辨脑电工作原理

从头皮电位求解皮层电位称为高分辨脑电。求解高分辨脑电需基于以下前提:大脑内部从皮层到头皮的区域(包括脑脊液、颅骨和头皮)是无源场,媒质是线性、分层均匀的。该区域的电磁过程可以用拉普拉斯方程来描述。利用边界条件,可以算出皮层电位。下面基于FVM,导出求解高分辨脑电的分层递推算法(layered recursive algorithm, LRA)。

首先引用N面体内的电流积分方程

其中▽φ可用φ来表示。根据场的无源性,式(3.13)等号的右端项I为0。因此,节点电流的连续性方程可以表示为

式中:N0代表三棱柱分层剖分网格中一层的节点数;n代表层数,由里到外逐渐增加;ai是矩阵的系数。

假设第n+1层和第n层的节点电位是已知的,将这些项移到式(3.14)等号的右边,则式(3.14)变为

式中的未知数是第n-1层的节点电位。式(3.15)也可以用矩阵形式表示为

据此,可由第n层和第n+1层电位递推出第n-1层电位。最外层的递推公式需要特殊处理,根据最外层的边界条件

可以得到最外层N和次外层N-1的节点电位之间的关系为

为形象起见,以二维四边形剖分网格为例加以说明,很容易将它推广到三维情况,甚至可以把二维情况简单地看作三维情况的截面表示,如图3.9所示。第N层代表头皮表面,第N-1层代表相邻第N层的剖分网格的次外层,其他各层依此类推。Pi代表剖分节点,实线的四边形代表体元。对任一节点作闭合面。以P5点为例,它的闭合面是B1B2B3B4。闭合面的每条边(注:对于二维情况,“边”和有限体积法的“面”表示同一几何实体)的电流通量用包含该边的体元上各节点的电位表示。例如,边B1C1和边B1D1的电流通量由节点P1P2P4P5的电位表示。但是,如果对头皮表面的节点作闭合面,例如对节点P8作闭合面B3B4B6B5,为了计算边C3B5B5B6B6C4的电流通量,需要在头皮表面以外添加虚拟的辅助层(第N+1层),用第N层和第N+1层上的节点电位表示其电流通量。我们已知头皮表面(第N层)的节点电位,但不知道第N+1层的节点电位,用这种方法计算高分辨脑电需要做一定的假设和近似。因此,利用边界条件,我们对于头皮表面的递推方程做特殊处理。在头皮表面的节点P8作新的闭合面C3B3B4C4,边C3B3B3B4B4C4的电流通量仍用第N层和第N-1层的节点电位来表示;由于头皮表面,边C3C4的电流通量为0。这样,头皮表面节点的电流通量只和第N层和第N-1层的节点电位有关,不需要再作辅助层(第N+1层)。

图3.9 头皮表面节点的闭合面

综上所述,最外层(第N层)的电位是测量值,是已知的,因此第N-1层的电位可以由式(3.18)推出;已知第N层和第N-1层的电位,可以进一步推出第N-2层的电位;依此类推,即可得到无源区的各层电位。这就是求解高分辨脑电的分层递推算法(LRA)。

注意,分层递推算法是综合考虑每层各节点间以及与上下相邻层节点间关系的基于层面的递推算法;而以往采用的是一点一点往里推的算法,它只考虑与该点相邻的各点之间的关系,并没有考虑到待求层面上各点之间的关系。我们将每层上的所有节点作为整体进行考虑,其准确度比逐点递推的准确度更高。这种分层的三棱柱剖分网格具有较小的矩阵规模,其系数矩阵的维数与一层的节点数相当,因此它的矩阵规模与三维问题中BEM的规模是一样的,而它的稀疏性和FVM相当。与正问题一样,采用稀疏矩阵的共轭梯度法计算,速度很快。同时,三棱柱剖分解决了六面体剖分的几何奇异性,在高分辨脑电问题中不存在几何奇异点。因此,LRA是处理高分辨脑电问题的一种理想方法。

3.2.2.2 高分辨脑电的模型

脑电正问题的计算可以分为两大类:解析解和数值解。

(1)解析解

对于球模型,可以采用解析解进行分析,下面简单介绍一下四层各向异性球模型的解析解。以下是电流偶极子源在多层球模型中产生的电位分布的解析解的推导公式。是系数因子。式(3.19)就是电流源偶极子在任意层球模型中产生的电位分布的解析解表达式。本章所用的四层球模型的解析解以及灵敏度的计算公式均由此式推得。

式中:φdip是待求点的电位;ε是介电常数;rθφ是待求点的球坐标;r0是偶极子球坐标的径向分量;Pn是勒让德级数;是一阶缔合勒让德级数;Mr是电流偶极子强度的径向分量;Mθ是电流偶极子强度的切向分量;

(2)数值解

各种球模型与头的实际结构形状相比都存在很大差异,用上述解析解会不可避免地带来误差。由于头是一个很不规则的几何形体,对于真实头模型,电位解析表达式的求解是不可能的,只能利用数值计算的方法解决边值问题。求解正问题的数值方法包括有限元法、边界元法、有限体积法、矩量法、模拟电荷法、有限差分法、虚拟媒质法等。下面首先对目前广泛使用的有限元法和边界元法做简单介绍,然后详细介绍有限体积法。

① 有限元法

有限元法(finite element method, FEM)是以变分原理和剖分插值为基础的数值计算方法。首先利用变分原理把所要求解的电磁场问题转化为等价的变分问题,然后利用剖分插值将变分问题离散化为普通多元函数的极值问题,最终归结为一组多元的代数方程组,解该方程组即可得到待求边值问题的数值解。

下面以二维静电场为例加以说明。它满足泊松方程

式中:φ是电位;ρ是电荷密度;ε是介电常数。将上面微分形式的定解问题转化为等价的变分问题,也就是泛函的极值问题。

式中:I是泛函;Ω是积分区域。对于Ω内无空间体电荷的情况,ρ=0,它实际上描述的是静电场中的汤姆逊定理:处于介质中一个固定的带电导体系统,其表面上电荷的分布,应使合成的静电场具有最小的静电能量。

在式(3.21)中,各单元中未知函数用该单元节点处的电位值展开,即有

式中:n0为单元的节点数;φi表示该单元中各节点i处的电位值;为单元ei点处的形状函数。

将式(3.22)代入式(3.21)以后,由于泛函中的ε, ρ都是事先给定的,形状函数也是事先选定的,所以泛函只是节点φi的函数。令泛函对各节点求导后等于零,即可得到一组联立代数方程,求解该代数方程组,即可得到各节点的电位值,这就是有限元法。其特点是对整个场域进行剖分插值,每个单元都认为是均匀的,而不同单元间的电磁特性可以不同。它适合对复杂场域进行剖分,比较全面灵活。Y.Yan应用有限元法分析了脑电正向问题,对三层球模型有限元解和解析解做了比较,并实际计算了真实头模型。

② 边界元法

有限元法是在整个场域设置待解量并进行离散化处理,而边界元法(boundary element method, BEM)则是在边界处设置待解量并进行离散。在媒质已确定的情况下,场函数是由域内的源与边界条件共同决定的。若域内的源已给定,则满足边界条件成为确定场函数的关键。边界条件反映了域外场或源对域内场的影响。边界元法首先将场的边值问题转化为边界积分方程问题,然后利用有限元离散技术形成边界量的代数方程组,解出边界量进而获得域内的场分布。该方法又分直接边界元法与间接边界元法。

直接边界元法的原理是位势理论。如静态电场问题,给定域内场源后,边界电位φ和与边界垂直的电场强度分量决定了域内任一点电位,并可得到其关系式。但是,边界上的φ通常并不全已知,它们也不是独立的。因此,需先找出边界上的φ的关系。例如,若已知边界上的φ,则可得到边界上的。两者均确定后,即可得出场域内任一点的场量。

间接边界元法的原理是场方程解答的唯一性定理。这时,未知量是分布在边界上的虚拟场源密度,如静电场问题为面电荷密度。寻求边界上场源密度,使满足给定的边界条件,则根据唯一性定理,解答必是正确的。对线性媒质进行离散化处理,通过集中量所产生的效应的叠加来表示边界上连续分布场源所产生的效应,从而形成代数方程组。解出虚拟场源后,便可求出场域中任一点的场量。

为形象起见,下面简单介绍要用到的边界元基本方程——位势方程。它可以表示为

式中:M0是待求电位的场点;MS; φ是电位;ε是介电常数;ρ是电荷密度;r是场点M0到dS的距离;rs是场源到场点M0的距离;θr和dS法向的夹角;Ω是场域;SΩ的边界;Sδ是场域Ω除去点M0以后的边界。这里只推导γ=1的情形。此时,式(3.23)可简化为

采用直接边界元法推导式(3.24)。由格林第二公式

将泊松方程

和格林函数v

代入式(3.25),整理得

式(3.27)中:δ为冲击函数,冲击点在Q,其满足如下性质

δP-Q)= ∞, P=Q

δP-Q)= 0, PQ

由于格林函数v在均匀媒质、无限大空间电磁场中的解答为,令u=φ,代入式(3.28)即可得式(3.24)。两种边界元法的特点包括:

其一,降低了问题求解的空间维数。三维问题可利用边界表面积分降维为二维问题,而二维问题则可降维为一维问题。代数方程组阶数降低,使许多问题得以简化。

其二,场强的计算准确度高。不必先解出电位分布,再求梯度以获得场强。可直接获得场强,避免了求梯度带来的误差。

其三,易于处理开域问题,即适用于场域外边界伸展至无穷远的情况。

然而,边界元法与有限元法相比存在以下不足之处。

其一,系数矩阵为非对称的满阵,求解大型代数方程组代价大。

其二,系数矩阵元素必须经数值积分处理,耗费机时较长。

其三,处理非线性媒质较为困难。

很多学者应用边界元法计算脑电正问题,其中王云华提出的“两步法”,有效地解决了分层媒质间电导率相差很大时的数值计算病态;Fuchs则改善了顶点格式的边界元法计算立体角的准确度。

③有限体积法

有限体积法(finite volume method, FVM)是从流体力学中发展起来的数值计算方法。它的中文译名有两个:有限体积法和有限体元法。为与有限元法从字面上相区分,避免误以为是有限元法应用于三维形体的计算,本章中均采用有限体积法这一中文译名。有限体积法最早引入到电磁场计算和生物医学工程中是1996年M.Rosenfeld发表在IEEE Transactions on Biomedical Engineering上的一篇文章。最早在国内介绍该领域FVM的是1997年尧德中教授发表在《国外医学·生物医学工程分册》上的一篇文章。下面具体介绍有限体积法。

假定人体组织中电阻抗的容性成分可以忽略,因此导体元问题是准静态和线性的。假定体积Ω和表面积S的电位分布由高斯通量定理给出

式中:σ是电导率;φ是电位;Jv是体元的电流密度。边界条件规定没有电流流入空气中,因此, n是边界的法线方向。

(a)离散化

物理区域被分成大量任意形状的六面体 —— 体元,这些分割可以产生合适的网格。两个相邻的体元如图3.10所示,每个元的表面根据曲线坐标的方向用ξ, η, ζ来表示。元的中心坐标为(i, j, k),8个顶点的坐标则相应地分别为i±1/2, j ±1/2, k ±1/2。元的表面方向由法向矢量Sl给定,这里l=ξ, η或者ζ。离散形式的未知数φi, j, k不像大多数离散方法定义在顶点上,而是定义在元的中心,以简化通量计算。

图3.10 体元和次级元示意

为得到准确的离散解,区域矢量由下式给出。

这里r是位置矢量。每个六面体的体积通过把它分成共用六面体的主对角线的三个金字塔形来计算。

(b)导体元方程的近似

在FVM中将积分守恒方程(3.29)离散化,得到

等号左边的每项均表示穿过相应表面的电流,等号右边表示单元中所有电流元的总和,V为体元的体积。用电流为零的边界条件可以消去边界处单元相应的电流流动项。

体元每个面的电位梯度由其积分定义

来计算,式中V为次级元的体积。这种方法优于微分形式。比起基于微分方程的有限差分法,积分形式的计算减少了网格平滑的限制。图3.10中的次级元可用于梯度的近似计算,它由两个体元的相邻各一半组成。下面的公式可以产生二阶精度的近似。

这里l表示次级元的一个面,φl是该面中心点的电位。例如,的近似公式为

因而,式(3.32)中对应的电流项为

类似地,可以推出用于计算其他各面电流通量的表达式。由此得出每个元涉及19个点的二阶精度的离散方程

有限体积法具有几个重要的特性。首先,离散FVM的数值解是保守的。在离散情况下,无论体元的形状如何安排,高斯通量定理在整个定义域仍然得到满足,因为它使用的是积分方程而不是差分方程。这一性质对于精确模拟电导率突变的问题非常重要。其次,FVM的计算量很小,因为体元离散化后的系数矩阵非常稀疏。FVM允许对相当精细的网格进行模拟,从而提高了给定计算结果的准确度。

头模型构造的研究是脑电问题研究中的一个热点,也是正问题计算的基础。通过建立的头模型,可以把脑内源的活动和头皮上测得的电位分布联系起来。建立头模型时,需要考虑头部不同组织的导电性,比如脑、头骨及脑皮层等。头模型建立的精确程度决定着脑电问题的定位精度。从以前完成的大部分工作来看,使用的头模型大致有球模型、椭球模型和真实头模型等几种形式。

球模型。最简单的球模型是均匀球模型,即用各向同性、电导率均匀的介质球来描述球的几何特性和电导率。当头模型为球时,模型中任意位置的偶极子所产生的电势在头皮上的分布都有解析解。也就是说,如果给出偶极子的位置、方向和极矩参数,就可以推导计算出模型表面任意一点电势的公式。假设模型是体积导体并且是线性的,则多个偶极子的电势分布可以通过简单的叠加来得出。均匀球模型会产生较大的定位误差,这是由于颅骨的低电导率不仅削弱了头皮电位的强度,而且使头皮电位变得模糊,其作用相当于对电位分布施加了低通的空间滤波。

三层同心球模型的使用使得定位准确性显著提高。在这个模型中,头被看作三个同心球,分别为头皮(scalp)、颅骨(skull)和大脑皮层(brain)三个同心球,又称之为SSB(scalp-skull-brain)头模型。为了更多地模拟头的组织,有的模型还包括了脑膜,这样就构成了四层同心球模型。

椭球模型。有些研究者采用椭球模型来研究脑电问题,希望通过采用扁长的或扁平的椭球模型,或者偏心球模型,来逼近真实头部形状,从而提高计算的精度。但研究结果表明,与采用球模型相比,采用椭球模型对计算的精度改进不大,故一般较少采用。

真实头模型。球模型被用于逆问题是由于它可以利用计算机快速求解,但球模型有两个缺点:球模型和人头部的真实形状有明显的差别,球模型只能作为对真实头的一个简单粗糙的近似;球模型的电导率一般假设为均匀的,与人头部的实际情况不相符。多层球模型克服了第二个缺点,但第一个缺点仍然存在。实际上,根据脑组织的结构构成可知,脑内每个点的电导率都是不一样的,即使它们属于同一组织。为了更好地体现人体头部的构造,人们引入了真实头模型。

真实头模型是根据人脑的实际结构(根据人脑的MRI、X-CT等成像手段或者解剖学知识)而得到的,它可以实际模拟人脑的结构,模型复杂。由于真实头模型更接近人脑的形状,与球模型和椭球模型相比,它的计算和定位精度显著提高,是目前头模型中研究的热点。真实头模型没有解析解,只能通过数值计算求解。在真实头模型上一般采用的数值计算有边界元素法(BEM)、有限元法(FEM;见图3.11)、有限差分法(FDM)等。

图3.11 有限元模型中的单元结构

为了更加精确地计算脑电问题,必须采用能够在复杂形状的头模型以及任意变化的体积电导率上进行仿真的数值解。而由于有限元计算可以出色地处理复杂的边界、形状和考虑各向异性,它已经是脑电计算中的主要方向。因此,在脑电建模中,如何建立贴近实际情况的有限元真实头模型是研究的热点问题。

为了得到真实头的有限元模型,首先需要有描述整个头部构造信息的数据源。通常,这些数据源可以通过对头部断层扫描,如CT或者核磁共振成像(MRI)等技术获得。而MRI图像由于可以较好地分辨出人脑内部的软组织,因此成了构建真实头有限元模型的首选数据源。一系列的二维MRI图像就相当于一组体数据,人们可以利用图像分割技术从每张二维图像中提取出各个组织,然后对每个组织进行三维重构,最后再将这些重构出来的组织的三维几何模型进行有限元剖分,获得真实头的有限元模型(见图3.12)。

图3.12 ANSYS中生成的3D头模型(头皮、颅骨、大脑三层组织)

对于真实头有限元模型的建立,研究人员已做了大量的尝试。Ziolkowski提出了快速二维有限元分割的方法。Hartmnna则提出了快速三维有限元模型建立的方法。在这些研究的基础上,还有一些可以用于有限元剖分的软件,如NETGEN35 I1等。

然而,由于人脑内部构造的复杂性和医学图像成像技术的特殊性,对MRI图像的分割以及对人脑组织的重构工作在具体实现上都还有很大的困难。特别是将数据量的限制考虑进去时,对脑灰质、脑白质等表面复杂组织的重构工作就变得更具有挑战性。

目前,边界元真实头模型在脑电研究中被广泛地应用,边界元模型的构造方法已经十分成熟,构造出的模型也相当稳定。

有限元真实头模型,一般是从头部的CT数据出发,经过信息提取、图像边缘检测、图像分割、三维重构重建及三维有限元体积剖分等步骤而完成的,其中涉及图像处理、模式识别、图像分割、三维图像重构和有限元体积剖分等多个学科的知识,建模方法繁琐耗时,而且建立的有限元模型还存在着一定的不稳定性,这些都极大地限制了有限元法在脑电脑磁研究中的进一步应用。

3.2.2.3 噪声对高分辨脑电的影响

实际采集的脑电数据中不可避免地会存在噪声,下面分析噪声对高分辨脑电的影响。以三层球模型为分析对象,如图3.13所示。

图3.13 噪声对高分辨脑电的影响

图3.13(A)表示头皮表面不加噪声时的皮层电位分布。

图3.13(B)表示头皮表面各点均加上0.001mV的微小电位变化量(它相当于头皮表面电位最大值的1%)所得到的高分辨脑电,它与图3.13(A)相比没有明显变化。

图3.13(C)表示头皮表面加上沿x轴正弦变化2%的空间低频信号(以球模型直径为一个空间周期)所得到的高分辨脑电,与图3.13(A)比较也没有明显变化。

图3.13(D)表示在球模型表面的一个点加上该点电位幅值2%的噪声,这时整个高分辨脑电图的分布发生了明显的改变,变化幅度最大的点达到不加噪声时原皮层电位的62%。

在图3.13(B)和(C)两种情况中,可以认为所加的噪声都是空间低频、慢变的噪声,而图3.13(D)的情况相当于施加了空间高频、快变的噪声。因此,可以得出如下结论:高分辨脑电对空间低频噪声不敏感,但空间高频噪声对高分辨脑电的影响很大。高分辨脑电相当于空间上的高通滤波器,它会增强高频成分,抑制低频成分。这个结论十分重要,也是很容易理解的。

3.2.3 脑电源定位问题

脑电正问题是由给定脑电活动的分布正向求解头皮电位。研究脑电正问题的目的在于,将正演结果与实际观测的头表数据(脑电图)或皮层表面电位数据(皮层电图)相比较,以达到辅助推断脑电发生源的目的。由此可见,脑电正问题是脑电逆问题的基础。

脑电正问题的研究包括脑模型和算法两个部分。真实头模型一般能够借助CT和MRI来构造,算法则因模型的不同而不同。真实头模型只能用数值方法求解,同时由于数值方法本身隐含的近似性,以及计算机真实头模型与实际头模型的差异,会产生一定的误差。此外,真实头模型的数值求解计算量大,且不便于分析因果关系,因此,各种近似头模型,尤其是球模型、多层球模型得到相当普遍的应用。

关于脑电逆问题的研究主要集中在两个方面:源的等效偶极子定位和不同的脑电源定位方法。

3.2.3.1 源的等效偶极子定位研究

偶极子源定位方法的目的是寻找产生头皮脑电的等效偶极子源(Mosher et al.,1992)。此方法中常用的模型有瞬时状态偶极子模型、时间—空间偶极子模型和由Pasual-Marqul等人提出的权重最小化方法中运用的三维网格模型。

(1)瞬时状态偶极子模型,实际上是一种空间模型,即只取某一时刻的各导联数据进行偶极子分析,一般取特殊时刻的数据,比如波幅最大或最小的时刻、在癫痫病中棘波出现的时刻等。这种模型一般用于寻找单个等效偶极子,对于两个或多个等效偶极子,很难得到稳定的有意义的解。运用此模型每隔一定的时间间隔取出各导联数据求出等效偶极子,可以进行偶极子跟踪。

(2)时间—空间偶极子模型(Christopher et al.,2000),即取一段时间观测记录。该模型又分为三类:移动偶极子、位置固定偶极子和位置与方向固定偶极子。位置与方向固定偶极子时空模型,即假定在观测的时间段内偶极子源固定,这段时间内头皮电位的变化仅仅源于偶极子源的强度变化,这样EEG就分为了空间成分(偶极子源的位置和方向)与时间成分(强度随时间而变化)。这三类模型最终都能归结为统一的表达方法,都能很好地利用观测数据的时间和空间信息。

(3)由Pasual-Marqui等人提出的权重最小化方法中运用的三维网格模型。该模型把10000个网格均匀分布在头内,每个网格点上有3个正交的单位强度偶极子源,通过不断减少最低电流密度区域的网格点的权重,来寻找最高电流密度区域。这种模型不必预先假定偶极子源的数目,但病态特性较重(Roberto et al.,1994)。

使用较多的算法有非线性最小二乘法、全局优化算法、子空间分解算法。

非线性最小二乘法即求函数FP)=‖V′-V‖取极小值的解(Salu et al.,1990),其中P为偶极子的参量,V′表示从N个头皮位置计算出的电位,F表示观测到的电位。此求逆过程就是FP)达到最小的优化过程。常用的优化方法有最速下降法、Marquardt法和Powell法。①最速下降法就是以FP)的负梯度方向作为搜索方向迭代寻优,每次迭代都先求得最优步长,当符合给定判据时,迭代结束。最速下降法具有线性的收敛速度,其特点是对开局有利,可以帮助我们很快接近最优点,但对收局不利,越接近最优点,最优步长越小,收敛越慢。②Marquardt法就是,如果初始量远离最优值点,则选较大的步长,按梯度法迭代到最优点附近,然后减少步长,按牛顿法寻找最优。这样既提高了算法的收敛速度,又防止了算法的发散。但相对梯度法而言,其计算量较大。③Powell法就是在导数的计算很困难时进行直接搜索的寻优方法,其特点是不必求导数,计算简单,但由于它利用的目标函数的信息少,所以迭代步骤较多。

全局优化算法借鉴物理、生物等学科的原理来克服局部最小问题,达到全局优化的目的。例如,模拟退火法就是借鉴热力学中的退火原理,即初始温度足够高,降温速度足够慢,能使模拟退火自动收敛到最优解。遗传算法是借鉴生物界自然选择和遗传机制,自适应地获取和积累搜索空间的知识,以控制搜索过程和达到局部最优解。

子空间分解算法被运用在多个偶极子源的定位中。首先用主成分分析(PCA)将脑电分解为信号子空间和噪声子空间,然后用独立分量分析(ICA)将信号分解,再由多信号分类方法(MUSIC)分别去定位(田源,2002)。MUSIC同时利用观测记录的时间和信息,用单个偶极子扫描整个大脑三维空间,得到多个偶极子的分布位置。子空间分解算法通过对多个偶极子分别定位,降低了搜索的复杂性。

3.2.3.2 不同的脑电源定位方法研究

在许多实际情况中,源的活动范围很广或其分布情况无法获取,因此,人们尝试重建可以了解神经活动的图像方法。理论研究(Cuffin,1985)表明,头皮电位分布可以被看作大脑皮层电位或源的分布,经头部组织(颅骨)空间卷积的结果,因此从头皮分布重建大脑皮层电位分布图像(或者源的分布),可利用空间解卷的方法实现。但要得到空间卷积阵,必须知道确定的源的分布。Nunez等人(贺士娟,2001)在平行平面(测量面与源分布面平行)的假设基础上,建立了空间解卷的理论公式,但其存在限制条件,即要求电极个数必须与源的离散化个数相等。

Freeman在同样的条件和限制下,由大脑嗅觉皮层电位分布对脑电源的深层分布进行了恢复(Gond,1986)。但是,在存在大面积皮层电活动的情况下,必须假设大量的源,此时解卷问题便成为一个严重的病态问题,从理论上是不可能的(Freeman,1980)。

Hjorth于1982年提出了不基于任何模型的空间解卷方法(Holder, 1987)。他利用多道脑电信号的空间相关阵来求得脑电的空间卷积矩阵,再用其对信号进行抽象的解卷来恢复。但是,从理论上讲,只有在源确定的情况下空间卷积矩阵才是确定的(Nicolas,1976)。因此,这一方法效果有限。Sidman等人于1988年提出了一种基于奇异值分解的皮层成像方法(Majer, 1991; Baker et al.,1991)。但这种方法只能对深部活动源在大脑皮层的表现进行比较好的恢复,对浅层源则无能为力。

皮层成像方法(Baker et al.,1991)先假定有一个半球形偶极子设置在一个各向同性的球模型中,即头模型由三层构成。最外层是单位半径的球,代表头皮;中间层是同心半球,代表大脑皮层,其电位是重建的;最内层也是同心半球,代表偶极子层。

He等人(He,1996; Wang & He,1998)认为,封闭的偶极子层能更好地等效代替分布于其中的源,他们采用球面偶极子层代替半球面偶极子层作为等效皮层源;为了克服脑壳低电导率的影响,采用同心三层非均质球作为头模型;为确保偶极子层包含皮层中的浅层脑电源,扩大偶极子层到皮层表面处,同时增加偶极子数目以确保偶极子密度足以反映其中的脑电源;通过仿真发现,采用偏心率在0.75~0.80,包含1200~1300个偶极子的偶极子层能取得较好的效果。为建立更为精确的头模型,Babiloni(1997)等人利用磁共振成像(MRI)建立具有真实几何形状的三层头模型,各层边界面由MRI通过三角剖分建立;采用包含364个等效径向电流偶极子的偶极子层作为源模型,偶极子层位于皮层4~10mm处,其形状与皮层表面相似。

实现皮层成像的方法还有很多。例如,He(1997)将头皮Laplacian值用于重建皮层电位,Srebro等人(1993)提出基于边界元的皮层成像方法, Gevins等人(Le & Gevins,1993; Gevins et al.,1994)将有限元方法用于真实形状、不均匀电导率头模型的皮层成像。詹望和杨福生(1999)提出建立在有限网络模型上的递归解卷算法,来逆推皮层电位。

3.2.3.3 脑内源定位的临床应用

脑电逆问题能消除颅骨等组织低电导率的影响,从而得到仅从头皮电位分布所看不到的信息(Baker et al.,1991)。在病灶定位方面,脑电图虽然可以反映脑部异常,但不能很好地进行定位;而偶极子定位则有明显的优势,它可用于癫痫定位(Schneider,1972; Sidman et al.,1982),也可用于局部脑血管病及肿瘤等(Watanabe & Sakai,1984; Lee et al.,1988)。根据对不同刺激的脑电诱发响应,可以利用偶极子进行功能定位(Maier et al., 1987)。通过追踪源随时间的变化情况,可以进行脑的高级神经活动(如思维)及心理学的分析研究。例如,据世界卫生组织(WHO)报告,全球癫痫患者约有5000万,3500万人未接受任何治疗或规范化治疗,癫痫在发展中国家存在很大的治疗缺口,百万患者中未接受治疗者占80%~94%。现在,对难以治愈而需要进行病灶切除手术的癫痫病人,在手术前往往会运用偶极子定位方法的脑电图记录并定出病患发作期间颅内痈源发放位置。在手术中以脑皮质电图(EcoG)及深电极图DEEG记录确定致痈区,对比偶极子定位方法,两者之间误差极小。偶极子定位方法的EEG具有无创性、准确性,相当于脑磁图(吴若秋等,2005)。EcoG能验证术前头皮EEG的可靠性,它是脑电生理活动的直接反映,由于去除了颅骨、硬膜的影响,灵敏度高,局灶容易识别,并可进行电刺激记录。