基于信号子空间主特征向量的心磁源重构

李浩然，朱俊杰

河南理工大学电气工程与自动化学院，河南焦作 454000

前言

心脏磁场研究的发展离不开超导量子干涉仪的成功研制，利用测量到的心脏磁场数据重构心脏的电流源，可视化心脏电流空间分布一直是心脏电生理活动逆向问题研究的重点和难点［1-5］。在心磁源重构的研究中有很多重构心脏电流源的方法，利用高精度的源重构算法重构心脏电流源，分析心脏电兴奋时电流分布的情况，可视化心脏电活动的传导过程，这一研究对心脏疾病的早期预防和诊断有重要的意义［6-7］。

波束形成器是一种常被用于心磁图和脑磁图源重构的成像方法。1997年，Van Veen等［8］在脑电研究的过程中利用线性约束最小方差波束形成（Minimum Variance Beamforming,MVB）方法对脑电信号源进行了定位；2001年，Sekihara 等［9］提出信号子空间特征投影的方法，并将该方法用于脑磁源定位的研究中；2006年，Kim 等［10］提出一种独立成分分析的预处理方法，利用该方法对预激综合征患者的传导旁路进行定位；2011年，王伟远等［11-13］针对QRS 波段心脏电活动变化快的特点，研究时间窗宽度、源移动速度以及信号噪声对源重构精度的影响；2014年，赵晨等［14］提出一种用于等效源重构的磁场极值信号法，该方法可以用于在非均匀介质的情况下求解心磁逆问题；2018年，周大方等［15］提出一种抑制空间滤波器噪声输出功率增益的波束形成方法，该方法具有很好的噪声抑制能力，提高了重建分布电流源偶极矩强度估计的源分辨率。

本研究提出一种改进的基于信号子空间主特征向量的波束形成（Improved Signal-subspace Principal Eigenvector based Beamforming,ISPEB）方法，此方法约束了空间滤波器输出的噪声功率和增益，补偿了信号主特征向量对噪声空间分布均匀性的影响。仿真实验中，通过SPEB 和ISPEB 方法对心磁仿真数据进行源重构，估计其对应的分布电流源空间谱。最后采用MVB、SPEB 和ISPEB 方法对2 例沿z 轴正向测量的健康受试的单周期61通道心脏磁场数据进行源重构成像，并比较不同两种方法的电活动成像结果。

1 波束形成方法

1.1 问题的提出

假设SQUID 系统在人体胸腔表面测量到的心脏磁场由n个呈网格状分布的电流源偶极子源产生［16-20］，则检测平面对应的电流源模型可以用如下线性方程表示:

其中，b(t)=[b1(t),b2(t),…,bc(t)]T是t时刻c通道阵列信号的列向量；q(t,rj)=q(t,rj)η(t,rj)(j= 1,2,…,n)表示n个分布电流源的偶极矩；rj=(xj,yj,zj)(j= 1,2,…,n)表示电流源空间的任意网格点位置；q(t,rj)=表示源偶极矩的强度；η(t,rj)=[ηX(t,rj),ηY(t,rj),ηZ(t,rj)]T是源偶极矩的方向；v(t)是t时刻测量的噪声向量。下文中，q(t,rj)、η(t,rj)、q(t,rj)、v(t)和b(t)将简写为qj、ηj、qj、v和b。

本研究使用空间滤波的方法重建心脏电流源，将心脏磁场测量数据b(t) 作为空间滤波器的输入，估计的分布源的偶极矩(t,rj)(j= 1,2,…,n)为输出，空间滤波可用加权的线性运算表示：

其中，W(t,rj)T是空间滤波的权矩阵，式（2）可简写为：

MVB方法的基本原理是先用空间滤波技术重建心脏的分布电流偶极子源，根据可描述电流源偶极矩平均强度的分布电流源空间谱估计，对心脏电流源成像。电流源偶极矩与该位置上的电流源强度有关。

带入式（3），得到空间滤波器的总输出功率为：

其中，tr[ ·] 表示矩阵的迹。假设式（1）中测量信号的噪声v(t)为高斯白噪声，满足E(v)= 0 与E(vvT)=σ02I，σ02是滤波器噪声输入功率，I是单位阵。将式（1）和E(v)= 0，E(vvT)=σ02I带入式（4），可以表示为：

其中，前一项为空间滤波器输出的所有电流源偶极矩qi对应功率和；后一项σ02tr[WjTWj]是滤波器输出噪声功率，其中tr[WjTWj]是输出噪声功率的增益。对式（5）开方，可得估计的电流源偶极矩平均强度空间谱：

式（6）中估计的电流源空间谱也可以估计为：

最小化空间滤波器输出的总功率，归一化空间滤波器输出的噪声空间谱强度，得到MVB 方法的空间滤波权矩阵：

SPEB 重建分布源方法的滤波权矩阵以MVB 滤波权矩阵为基础。将MVB 权矩阵Wj,MVB的列向量投影到信号子空间主特征向量x1为基的一维子空间，可得到新的矩阵：

其中，x1是对称阵E(bbT)最大特征值λ1(λ1∈(0,1))对应的主特征向量。满足x1Tx1= 1，属于信号子空间。将式（9）带入式（7），令E(bbT)=I，对MVB的权矩阵Wj,MVB使用信号子空间主特征向量投影策略后，新的噪声空间谱强度为：

用式（10）对应的参数标准化权矩阵Wj,1，归一化每个分布源位置的空间滤波器输出的噪声空间谱强度，使其重新分布均匀，归一化后的SPEB 空间滤波器权矩阵为：

1.2 ISPEB方法

本研究提出ISPEB方法，令空间滤波权矩阵为：

可以证明，其中的实对称矩阵V=E(bbT)(x1x1T)-1是一个对空间滤波器输出功率有影响的矩阵。由式（9）可知，因为最大特征值λ1相比其他特征值，受噪声影响最小，故可以采用该功率较小的滤波权矩阵来增强分布源重建与空间谱估计的抗噪能力。因此，ISPEB方法中输出噪声功率增益可以表示为：

对式（13）使用施瓦兹不等式可以得知，ISPEB方法可以更好地降低空间滤波器的输出噪声增益，将E(bbT)=I带入式（8）、（11）和（12）可知，ISPEB、SPEB 和MVB 的噪声空间谱强度都相等，等于1。综上所述，ISPEB方法可以约束噪声空间谱对分布源空间谱的影响和滤波器输出功率的增益，相对于MVB和SPEB 方法，ISPEB 方法可以提高估计空间谱的分辨率。

1.3 ISPEB方法的分析与比较

由于多电流源重建问题比较复杂，采用波束形成方法进行电流源的空间谱估计决定了电活动成像的分辨能力［21-24］。本研究比较了ISPEB和SPEB单电流源重建的源分辨率。

由式（2）可知，当单电流源偶极矩的方向已知后，空间滤波器的源偶极矩估计退化为源偶极矩的强度估计，其中，和wj是退化后空间滤波器的输出和权向量。由式（7）可知，任意位置的单电流源空间谱强度估计为假设给定单电流源S 在单源rs位置上的强度估计为。定义任意位置rj(j= 1,2,…,n) 上的点扩散函数为归一化后的估计空间谱可以表示为

由式（12）可知，单源重建时，采用的空间滤波权向量为：

由式（7）和（14）可得ISPEB方法的单源重构的空间谱强度估计：

其中，lj是电流源产生的磁场列向量。将式（15）进行化简，得到ISPEB的空间谱估计强度表示如下：

归一化后的点扩散函数为：

同理可得SPEB的点扩散函数：

比较式（17）和（18），有：

通常，由于所有测量通道上的理想磁场信号平均功率大于噪声平均功率，所以可以令α=根据施瓦兹不等式(ljT,f)2≤可知，由此，根据式（18）和（19）可得：

当rj=rs时，点扩散函数可以取到最大值1，但此时仍然小于1，所以永远小于1，即比小。当rj≠rs时，其他空间位置上，ISPEB 方法的点扩散函数仍然小于SPEB 方法的点扩散函数，点扩散函数可以反映空间其他位置的估计强度影响大小，取值小，说明单源对邻域的影响扩散小。定义空间谱估计的单源分辨率为ψj=将式（20）带入该式，得到：

因此，相比于SPEB 方法，ISPEB 方法有着更好的源分辨率。

2 仿真与实验

2.1 仿真数据和相关源重构

相关电流源时比较难分辨，本研究利用仿真的磁场数据，比较了MVB、SPEB 和ISPEB 方法估计相关电流源的能力。

在仿真实验中，采用双电流偶极子仿真产生61通道的磁场数据（图1），在3 组仿真数据中分别加入信噪比为10 dB 和20 dB 高斯白噪声，源模型均在检测平面下8 cm 处，检测平面范围为20 cm×20 cm。分别用3 种方法对双电源偶极子进行源重构，并分析3种源重构的精度，采用重构出的电流源极矩强度较大的位置与给定电流源的位置误差和重构电流源的分布面积作为评价源重构精度指标。

图2是3 种方法对双电流源偶极子重构得到的空间谱强度估计经过归一化处理后的等高线表示，图中的“×”表示给定的源位置。可以看出，在20 dB的信噪比下，SPEB 和ISEPB 方法明显比MVB 方法有着较好的源重构能力；在10 dB的信噪比下，ISPEB方法有着最好的源重构能力。在低信噪比情况下，ISPEB方法的源重构能力最强。

图3是SPEB 和ISPEB 方法对双电流源偶极子重构的噪声输出功率增益抑制比。当信噪比为20 dB 时，SPEB 方法的噪声输出功率增益抑制比为627～657，ISPEB 方法的噪声输出功率增益抑制比为1625～1 780；当信噪比为10 dB 时，SPEB 方法的噪声功率增益抑制比为635～657，ISPEB 方法的噪声输出功率抑制比为1 662～1 780。通过对比分析可知，当信号的信噪比从20 dB 下降到10 dB 时，ISPEB 方法的抑制能力更突出，因此，在多电流偶极子源重构中ISPEB 方法的噪声抑制能力和源分辨能力更好。

2.2 心脏电活动成像

本文使用的心磁图数据采集协议已得到当地伦理委员会批准，同时，受检者给出了书面知情同意书。对两名健康受试的61通道心脏磁场数据进行源重构分析。图4是两组健康者的心脏磁场数据图，在整个心脏活动中，心脏磁场信号的R峰时刻心脏电活动较为剧烈，此时心室内电活动比心房的电活动强，易于观察。分别用MVB、SPEB 和ISPEB 这3 种方法对两名健康受试的心脏磁场数据R 峰时刻进行源重构。

图5给出了两名健康受试实测心脏磁场数据的源重构结果，图中方框表示给定的源位置。可以看出，MVB方法重构出的心脏电流源分布效果最差，心脏磁场信号中的噪声严重影响了源重构的结果。在信噪比较低的情况下，SPEB 和ISPEB 方法都有较好的成像结果。在R峰时刻，健康受试的心房除极已结束，心室正处于除极期，心室内黄色表示的电活动强度比心房内红色表示的电活动强度高，心室与心房的源强度对比度增强了，心室电流源在心脏范围内的分布及其电活动成像的清晰度也改善了。MVB方法成像结果比较模糊，上述电活动特征不明显；而SPEB方法和ISPEB方法都有较好的清晰度。

3 结论

本研究在MVB 方法和SPEB 方法的基础上，在SPEB 方法的权矩阵上添加了一个约束矩阵，使其在信噪比较低的情况下有着较好的成像结果，ISPEB方法利用一种低迹半正定矩阵构造了一个滤波器权矩阵，可以有效地降低噪声功率增益，提高重建分布电流源偶极矩强度，即提高电流源空间谱估计的源分辨能力。实验结果和分析比较表明，当信噪比较低时，ISPEB 方法在3 种方法中有着最好的成像结果，当信噪比较低时，可以采用ISPEB方法进行成像。