P波间期的心脏电流源重建及电活动磁成像*

周大方 蒋式勤† 赵晨 Peter van Leeuwen

1)(同济大学电子与信息工程学院,上海 201804)

2)(维藤/黑尔德克大学健康学院,德国维藤 D-58448)

1 引 言

1970年,Cohen等[1]利用超导量子干涉仪(superconducting quantum interference device,SQUID)在人体胸腔表面无接触、无创地测量到了心脏磁场图(magnetocardiogram,MCG).心脏磁场十分微弱(约10–12Tesla),远低于地球磁场(约10–6Tesla).SQUID的优点是[2,3]: 心磁测量无需外加激励,无需接触式电极,受试个体完全无损,操作方便.尤其是测量磁信号的时间分辨率高(≦1ms).目前,SQUID测量MCG信号的多通道系统已从早期的4通道发展到300通道.利用心磁图无创获取心脏功能信息的研究也取得了实质的进展[4−7].虽然心磁图的研究比心电图(electrocardiogram,ECG)推迟了80多年,但是,心脏电/磁信号同生共源,心磁图研究可以作为心电图的一种补充.MCG对心脏的切向电流(tangential current)和再入电流(reentrant current)较为敏感,可提供心脏的功能信息[3,8].

重建磁场电流源需要求解场–源逆问题[9,10].近20年脑磁(magnetoencephalogram,MEG)[11−13]和脑电(electroencephalogram,EEG)[14,15]的研究中发展了一种基于空间滤波器的重建分布电流源的方法.该方法给研究对象划分网格,并在每个网格交点上构造一个空间滤波器.然后,重建滤波器位置上电流源的偶极矩,获取分布源空间谱的信息.研究表明,重建分布源的精度与空间滤波器的输出噪声有关[16,17].虽然心磁测量是在屏蔽室内完成的,但是,在信号测量、数据处理和心脏电流源重建的环节中,仍需采用抑制噪声的方法[18−22],尤其是当P波间期信噪比(signal-to-noise ratio,SNR)较低(≤ 30 dB)的时候.

最小方差波束成形(minimum variance beamforming,MVB)是一种改进的空间滤波方法[11].为了提高重建电流源的精度,MVB采用了空间滤波器的输出总功率最小化和噪声空间谱强度归一化,以及在线调整权矩阵中测量磁信号的自适应技术[23,24].因此,MVB重建电流源的精度比非自适应的最小范数空间滤波(minimum norm spatial filtering,MN)方法高[25,26].Sekihara等[11,13]还提出了一种可降低MVB空间滤波器输出噪声功率的信号子空间投影方法,并将其用于脑功能的磁源定位.脑磁源定位的方法需要假设脑的有效电流源数目少于磁场测量通道的数目,以及综合导联场矩阵(composite lead-field matrix)是列满秩的.然而,由于心脏结构复杂且有动态变化,一般情况下这两个条件不能同时满足[17].为此,我们曾经提出了一种抑制空间滤波器输出噪声功率增益(suppressing spatial filter output noise-power gain,SONG)的波束成形,可用于R波峰时刻心脏电活动成像[27].

对心脏P波间期电活动无创成像的难点是P波心脏磁信号比R波峰弱,信噪比较低.针对这一问题,本文提出了一种基于空间滤波器的可提高分布源空间谱估计强度对比度(improving intensity contrast of distributed source spatial spectrum estimation,IIC)的波束成形.因为分布源模型中导联矩阵表示传感器阵列对一个源所在位置的测量灵敏度,所以,IIC方法在空间滤波器的加权矩阵中引入导联场矩阵,可以使滤波器输出估计对磁场电流源及其分布比较敏感,从而改进分布源强度的对比度.再通过设置源强度阈值,提取每个时刻重建源中偶极矩强度极大的电流源,就可以消除其他位置上相对弱的重建分布源与超出心脏范围的伪源,提高P波间期电流源重建的精度.文中给出了IIC方法单源重建的理论分析与仿真试验,并比较了该方法与MVB和SONG电流源重建方法的性能.文中还用2个健康人的P波段61通道心脏磁场测量数据重建心脏电流源,并分析了IIC方法的成像结果.

2 方 法

2.1 磁场分布源的模型及其空间谱估计

假设SQUID系统在人体胸腔表面测量到的心脏磁场由n个呈网格分布的电流源产生.测量面上第k通道的心磁测量信号用标量bk(t)表示,则t时刻c通道阵列信号的列向量可用b(t)=[b1(t),b2(t),...,bc(t)]T表示.磁场与分布源的关系模型可用如下线性方程表示[13,16]:

心脏电活动是一个随机过程[3,17],故文中假设电流源的偶极矩是一个随机变量,简记为个分布源的空间谱可表示为

本文利用空间滤波器方法重建产生磁场的电流源.已知空间滤波器的输入是测量信号向量输出是需要计算的电流源偶极矩估计它们的关系可用线性加权运算表示[13,24]:

2.2 IIC波束成形方法

本文在MVB[12]和SONG[27]的基础上,提出一种IIC波束成形方法.目的是通过提高分布源强度的对比度,减小P波间期信号噪声的影响.

2.2.1 滤波器的加权矩阵与比较

MVB和SONG的滤波权矩阵分别为[12,27]:

为了改进分布源强度的对比度,本文在此基础上,构造了IIC的滤波权矩阵:

根据实对称矩阵有关定理[28,29]可以证明,(5)式中的实对称阵U是一个“特征值不大于1,矩阵的迹小于其阶数”的低迹的半正定阵.且有

如果空间滤波器的输入v是高斯白噪声,则分布源空间谱的强度估计由(7)式可见,当空间滤波器的输入v是高斯白噪声时,IIC方法降低了空间滤波器的输出噪声功率增益和输出噪声功率且降噪性能比SONG和MVB好.

2.2.2 空间谱估计强度对比度的分析

由于多源空间谱估计分析的困难,本文分析了用IIC方法的单电流源空间谱估计.首先定义波束成形方法的点扩散函数为任意空间位置上对电流源强度估计的归一化.假设任意位置上分布源和上给定单电流源s的强度估计分别为并定义源强度的对比度与源强度的对比度成反比.

和

建立科学有效的评价机制。打破传统的把考试分数作为唯一标准来衡量学生的教学评价方法，采用非标准化考试。从注重学习知识的考评向运用知识的考量转变，主要考量学生运用知识分析解决问题的能力；从注重最终结果的考核向学生成长过程考核，主要考量课堂参与度、线上学习次数和时长、团队汇报展示、团队作业、随堂实验等，同时客观记录并科学评价学生成长历程，学生的创新实验、技术研发、发表论文、研究课题、获得专利、竞赛成绩和自主创业等都要记录在学生成长档案；由教师单一评价主体向自评、互评、他评和教师评价等多元评价主体转变，进行综合评价。

可见,这3种波束成形方法中,IIC的单源空间谱估计强度的对比度最高,即

2.2.3 局部强度极大电流源提取

本文在改进空间谱估计对比度的基础上,采用了设置源强度阈值,提取每个时刻重建分布源中局部强度极大电流源的方法.这是因为在重建磁场分布源时,需要消除环境噪声和计算误差引起的大量分布的重建弱源与重建的伪源(可能超出心脏范围).本文假设源空间中每个分布源的最小邻域中包括26个其他分布源,图1中它们的位置分别用红色和蓝色表示.并用每个时刻空间谱估计的最大强度作为提取局部极大源的阈值的基数,其中比例系数可根据经验确定.这样,改进源强度对比度后,可以通过设定阈值,消除部分重建的伪源与大量较弱的重建分布源.最后,利用那些局部强度极大的电流源来研究心脏的电活动.

图1 分布源位置及指定的26个邻域位置,分别用红色和蓝色表示Fig.1.The position of the source and 26 specified adjacent positions are shown in red and blue colors,respectively.

3 心脏电活动的仿真与成像

电活动的仿真与成像是指用仿真的磁场数据重建电流源,并记录电流源在空间中的位移,给出源位移随时间变化的图像.目的是检验用局部强度极大电流源近似磁场电流源,以及重建电流源在空间中的移动模拟电活动成像的效果.本文将根据重建电流源中伪源的个数、重建源的误差、以及成像结果中电活动的方向,分析比较IIC,SONG,MVB以及一种非空间滤波算法即信赖域反射(trust region reflective,TRR)方法的仿真与成像的结果.TRR是一种单电流源重建方法[30,31],它通过最小化磁场残差的范数和迭代更新,求解等效的单电流源[32].

3.1 仿真数据

假设包含心脏的人体躯干是沿坐标系Z轴水平分层的导体(horizontally layered conductor,HLC),心脏–躯干模型如图2(a)所示.层间电导率之差与该处静电势的乘积被称作二次电流密度,其方向与沿Z方向测量的心脏磁场平行[17].根据毕奥萨法尔定律,导体中二次电流密度产生的磁场可以忽略[10].对图2(a)中包含心脏在内的导体

划分网格,由(1)式可得,间距为1 cm的网格交点将构成一个n=12168的分布源空间.

3.1.1 磁场导联场矩阵

参照Magnes 1300C生物磁强计系统[2],假设仿真的1 kHz磁场数据是在胸腔表面用61通道环状阵列传感器沿Z方向测量的.如图2(a)所示,是通道的磁信号测量位置,故(1)式中的磁场导联矩阵可计算如下[17]:

图2 (a)水平分层导体.其中黑色直线表示分层躯干的边界,黄色椭圆体表示心脏;(b),(c)测量通道1和61的磁场导联,(ζ=X,Y,Z)曲线.其中红、绿和蓝色分别表示X,Y和Z方向导联的曲线Fig.2.(a)The horizontally layered conductor(HLC),where the black line indicates the boundary of HLC,and the yellow spheroid the heart;(b),(c)The X,Y and Z lead-field based plots of channels 1 and 61 are expressed in red,green and blue,respectively.

3.1.2 仿真的61通道磁场数据

首先,用已知偶极矩强度的单或双电流源产生模拟心脏的磁场.为简单起见,假定该单或双电流源通过了源空间中的直线参考路径,用于仿真沿参考路径的电活动.因为P波间期心脏电活动的强度先增后减,所以用半周期正弦函数近似描述该单或双电流源偶极矩强度的变化[17]:

其中rp是参考路径上某个分布源的位置.τrp是电流源出现在位置rp的时刻.ω是弧度频率,正弦函数的周期T=2π/ω.根据P波间期人体心脏的心房和心室肌中电活动的最大速度约1 m/s[33],令直线参考路径上相邻位置间电流源移动的时间∆t=τrp+1−τrp=0.01s,以及移动速度(moving velocity)MV=|rp+1−rp|/∆t≈1m/s.其中rp与rp+1表示直线参考路径上相邻位置的前后.当τrp≦t<τrp+1时,电流源出现在位置rp,其偶极矩强度q(t,rp)按照半周期正弦函数先增后减.当t=τrp+1时,该电流源到达位置rp+1,其偶极矩强度q(t,rp+1)将先增加后减小至0.

然后,在给定单或双电流源偶极矩和计算相应磁场导联矩阵的基础上,61通道的仿真磁场数据可用下式计算:

其中b(t)=[b1(t),b2(t),...,b61(t)]T是61通道磁场测量信号.q(t,rp)=q(t,rp)η(t,rp)表示单或双电流源的偶极矩.Lp是单或双电流源的磁场导联矩阵.v为高斯白噪声,根据文中P波间期心磁测量信号的信噪比给定.

如第2节所述,利用61通道的仿真磁场数据求解逆问题,即可得到分布源的空间谱.再从中提取局部偶极矩强度极大的电流源作为重建电流源,并给出其沿给定参考路径位移的图像.

3.2 仿真结果1

假设参考路径1是一条起止点坐标为Ps=(7,5,6)cm和Pe=(–3,–5,16)cm的空间直线.如图3(a)—图3(e)中黑色圆圈所示,该路径上相邻分布源的间距为单电流源通过参考路径上相邻位置的时间速度图3(b)—图3(e)给出了SNR=20 dB时,用上述4种方法得到的重建电流源和电活动成像结果.图中可见,用IIC每个时刻重建的电流源沿着参考路径1上的黑圈移动.黑圈中的绿色由浅到深表示时间变化.在该参考路径以外,IIC出现的伪源最少,SONG与TRR的结果次之.

表1中比较了上述4种方法的仿真结果中的分布电流源估计误差,分别考虑了无噪声以及信噪比SNR为30与20 dB的情况.表1中分别表示X,Y,Z方向分布电流源估计的根均方误差和总误差.无噪声时,IIC,SONG和TRR的源估计误差较小.MVB的误差相对较大.SNR=30 dB时,相比其他方法,IIC可以降低误差E.SNR=20 dB时,IIC的单源估计误差比其他方法小.如图3所示,IIC能够仿真单源沿参考路径的电活动,伪源数量最少,单源重建的精度相对其他方法要好.

3.3 仿真结果2

假设参考路径2包括两条起止点坐标分别为[Ps1,Pe1]=[(7,5,6),(7,–5,6)]cm和[Ps2,Pe2]=[(–3,5,16),(–3,–5,16)]cm的空间直线.如图4(a)—图4(e)中黑圈所示,这两条路径上相邻分布源的间距均为1 cm.仿真中,令两个电流源分别通过各自参考路径上相邻位置的时间速度其中一个电流源比另一个电流源提前0.005 s开始移动.图4给出了SNR=20 dB时,4种方法的仿真结果和电活动成像.图中可见,用IIC每个时刻重建的电流源沿参考路径2中的黑圈移动.黑圈中的绿色由浅到深表示时间变化.在该参考路径以外,IIC出现的伪源最少,明显优于其他3种方法.

表2比较了上述3种空间滤波方法仿真结果中的分布电流源估计误差,分别考虑了无噪声以及信噪比SNR为30与20 dB的情况.TRR是一种非空间滤波的等效单源重建方法,表中未列出.从表2可见,信噪比对双电流源估计的误差有影响.表中用IIC方法同时估计两个电流源的误差较小,而用SONG和MVB的估计误差相对较大.SNR=20 dB时,IIC的双电流源估计误差最大值E=2.44cm,比其他方法小.当SNR=30 dB时,IIC的双电流源估计误差相对SNR=20 dB时小.如图4所示,IIC亦能够仿真双电流源沿参考路径的电活动,伪源数量较少,其双电流源重建的精度相对其他方法要高.

图3 (a)黑色圆圈○和箭头表示沿直线方向的参考路径1;(b)−(e)是SNR=20 dB时,IIC,SONG,MVB和TRR方法的源重建结果.黑圈中的浅绿色到深绿色表示重建电流源位置随时间的变化Fig.3.(a)The black circles ○ and arrow indicate the reference path 1 along the straight line;(b)−(e)The results of the current source reconstruction using IIC,SONG,MVB and TRR methods when SNR is 20 dB.The green points at different color levels denote the reconstructed source locations at different times.

4 心脏电活动的磁成像

本文利用P波间期采集的磁场信号,分析了健康人A和B的心脏电活动及其成像结果.这些数据是在磁屏蔽室内用Magnes 1300C 61通道生物磁强计系统沿Z方向测量的.信号采样频率1 kHz[2].假设P波间期心磁测量信号及其噪声的平均功率为Ps和Pv,则信噪比 SNR=10lg[Ps/Pv].计算可得健康人A和B测量信号的P波间期信噪比 SNR∈[20,30]dB.

图5中给出了健康人A和B的P波峰值(Ppeak)时刻心磁信号的源空间谱估计强度等高线图.经归一化处理后的等高线值用色标表示.图中比较了用IIC,SONG和MVB三种空间滤波方法重建电流源的结果.为了显示人体心脏及其房室的相对位置,图中给出了一个健康人的心脏核磁共振影像(magnetic resonance imaging,MRI)冠状位(coronal)和水平位(transverse)视图,将其作为图5中的XY与XZ面视图,并分别与A和B的测量坐标系配准.图5(c)和图5(d)中A和B的心房内,用绿色三角形表示了冠状位与水平位面的交点(0,−9,10)和(0,−6,10)cm.

表1 参考路径1对应的单源估计误差Table 1.The error of single source estimation correspond to the reference path 1.

图5(c)和图5(d)中给出了Ppeak时刻重建电流源的结果.IIC图中,黄红色标和紧密环绕的等高线表示较强电流源均集中在A和B的心房.健康人心房内的电活动较心室内强.这与Ppeak时刻健康人的心房正在除极,而心室尚未开始除极有关[33,34].然而,SONG的等高线扩散到了心脏以外,不像IIC那样集中,说明SONG的空间谱强度估计的对比度相对较低,影响了重建电流源的精度.MVB的等高线图显示,心脏内外有很多无规律和强弱不同的伪源.重建电流源误差相对较大是MVB电活动成像结果不理想的主要原因.

图6(c)和图6(d)是健康人A和B的P波前半段(Ponset—Ppeak,Ponset是P波起点)电活动成像结果.图中用圆点表示重建源的位置,其中颜色由浅红到深红表示重建源位置随时间的变化.图6(c)和图6(d)中IIC的结果显示,P波前期健康人右心房除极时电活动从右心房的右上方向左移动.这与右心房的右上方靠近窦房结,右心房的左侧靠近右心室有关[33].由于P波后半段(Ppeak—Pend,Pend是P波终点)的磁信号比较复杂,心房内电活动成像的特征不明显,故文中没有图示.图6(c)和图6(d)中,SONG与IIC的电活动成像结果类似,A和B的心房外各有少量重建的伪源.图6(c)中IIC和SONG显示,虽然心脏P波前期磁信号呈单调上升,但是Ponset—Ppeak间期电活动的方向有变化.从起始时刻Ponset到结束时刻Ppeak,重建源移动方向如图中绿色箭头所示.MVB的成像结果比较模糊,所以A和B的心房内电活动特征不明显.TRR是一种重建等效单源的方法,所以图6(c)和图6(d)中显示,等效单电流源是从右上向左下移动的[32].实验研究结果显示,IIC图6中P波间期健康人A和B心脏内分别有14个和26个极大电流源.相比其他方法超出心脏范围的伪源较少.IIC和SONG的成像结果显示了健康人右心房中电活动的方向.

5 讨 论

用MVB,SONG和IIC研究源重建的结果表明,MVB采用的空间滤波器噪声空间谱强度归一化方法,SONG采用的抑制空间滤波器输出噪声功率增益方法,以及IIC采用的上述两种方法,可以不同程度地抑制噪声的影响,改进重建电流源的精度.表1、表2和图3、图4的仿真结果表明,IIC方法优于MVB和SONG.由(14)式和(15)式可知,IIC空间谱估计的源强度对比度比MVB源强度对比度大.

本文采用在滤波权矩阵中引入导联场矩阵的方法,研究了有多个电流源的情况.图3和图4的仿真结果说明,IIC可以较好地消除重建源中的弱源与伪源,改进了2个源重建的精度.然而,当多个电流源的偶极矩方向变化时,多源重建的精度还需要深入研究.

表2 参考路径2对应的双电流源估计误差Table 2.The error of two sources estimation correspond to the reference path 2.

图4 (a)黑色圆圈○和箭头表示沿直线方向的参考路径2;(b)−(e)SNR=20 dB时,IIC,SONG,MVB和TRR方法的源重建结果.黑圈中的浅绿色到深绿色表示重建电流源位置随时间的变化Fig.4.(a)The black circles ○ and arrow indicate the reference path 2 along the straight line;(b)−(e)The results of the current source reconstruction using IIC,SONG,MVB and TRR methods when SNR is 20 dB.The green points at different color levels denote the reconstructed source locations at different times.

IIC采用提高空间谱估计的源强度对比度和提取偶极矩极大源的方法去除重建的伪源.从图3—图6的源重建和电活动成像结果可见,这种方法比其他方法要好.仿真和实测成像的结果中,根据经验,重建源强度的阈值分别等于当前时刻电流源空间谱的最大强度估计的0.01%和1%.图6中用IIC时,健康人A的瞬时极大源数最多是2个,健康人B的最多是4个.

图5中2个健康人的IIC的电活动成像结果显示,参照MRI中人体心脏与其左右房室的位置,用黄色表示的Ppeak时刻较强电流源均集中在他们的右心房.图6中红色圆点表示Ponset—Ppeak间期强度极大源也集中在右心房.他们右心房的电活动均比左心房强.这与右心房靠近窦房结,即心脏电兴奋的起搏点,以及健康人右心房的电兴奋早于左心房相符合[33].图6(c)中IIC和SONG的电活动成像结果显示,虽然Ponset—Ppeak间期心磁信号单调上升,但是,这时电活动的方向有变化.可能这与参考文献[33]指出的“心脏肌纤维的走向呈螺旋状,以及在纤维的走向上,肌电阻率较低,肌电活动较强”有关.图6中未能给出左心房电活动的成像结果,可能与心脏磁场只有Z轴方向的信号测量有关.

我们还用IIC方法分析了2个病人的数据[4].他们心脏左、右冠状动脉严重狭窄,P波间期心房电活动与健康人明显不同.可能是心脏冠脉严重狭窄引起心肌缺血,他们的心房电活动相对较弱,未显示出心房中电活动的方向.因此,文中没有给出病人的电活动成像结果,有关问题还需要深入研究.

6 结 论

本文针对心脏P波间期电流源重建及心脏电活动成像的问题,提出了一种可提高分布源空间谱估计强度对比度(IIC)的波束成形方法.该方法在空间滤波权矩阵中引入导联场矩阵,可以改进分布源空间谱强度的对比度,有降低滤波器输出噪声功率增益的作用.再根据经验设置源强度阈值,提取重建分布源中偶极矩强度极大的电流源,可以去除重建源中相对较弱的分布源和心脏范围以外的伪源.可以提高P波间期电流源重建的精度.文中对4种方法的理论分析,源重建仿真和电活动成像结果的比较可见,IIC方法将有助于P波间期的心脏电流源重建和电活动的成像.对有局部心肌缺血的病人,心脏电活动成像的问题有待深入研究.

感谢中国科学院上海微系统与信息技术研究所的张懿教授,谢晓明教授和孔祥燕教授及张树林博士对本研究的支持和技术交流.