CN110393522B - 一种基于图总变分约束的无创心脏电生理反演方法 - Google Patents

Abstract

本发明公开了一种基于图总变分约束的无创心脏电生理反演方法，基于图总变分约束进行了无创心脏电生理反演，图结构的建立考虑了空间不相邻节点之间相似性，充分利用了心肌跨膜电位序列分布的底层特征，同时考虑了电兴奋传播的时间相关性，提供了高精度的心肌跨膜电位序列分布。本发明通过最终得到的心肌跨膜电位的时空分布，可以准确的观察到心肌梗死时的梗死区域位置及边界、异位起搏发生时的起搏点以及其它电活动异常的区域，在临床的诊断及治疗上具有重要的参考意义。本发明可直接同时重建激活的整个序列，而非按照单个时间点的求解。

Description

一种基于图总变分约束的无创心脏电生理反演方法

技术领域

本发明属于心脏电生理反演技术领域，具体涉及一种基于图总变分约束的无创心脏电生理反演方法。

背景技术

心率失常是威胁人类健康的最严重疾病之一，它可发生在心肌梗死之后的疤痕区域，也可自行发病。梗死的边缘区存在异质性，部分存活的细胞和缺血细胞间形成折返回路，极易造成异位起搏，引起心律失常。梗死愈合后，纤维化的细胞形成慢传导区域也会对心室的正常兴奋传递造成影响，边缘区域可能继续扩散影响正常心肌细胞，直至引发心力衰竭。因此，准确的定位梗死的位置与大小，描绘其边缘区，定位异位起搏点，可以揭示患者心律失常的风险，以预防更严重的心脏疾病。

临床上应用较多的识别病变位置的方法为侵入性导管起搏指导法，将导管放置在心内膜中并进行起搏，以模仿所寻求的搏动，然后将预先诊断为异位起搏的心电图的模式与导管起搏产生的心电图进行比较，相似程度最高的即为目标病变点，以此指导后续的消融手术；但这种方法手术时间长，对病人造成的压力较大，对医生的经验要求较高。因此，用无创的方式精准的定位病变区域，可以减少手术的难度与时间，对医生与病人来说都是迫切的需求，是一个十分值得研究的内容。

通过对人体表面测量的电位的数学处理来推断心脏的电信息，我们称之为心电逆问题，这是一个较为困难的问题，它具有物理和数学上的不适定性。由于测量数据的维度(体表电极的数量，通常为12/64/128导联)远小于未知量的维度(心脏上的节点数，约为10³数量级)，得到的解往往是不唯一的。常见的解决方法有基于2范数的吉洪诺夫正则化、基于1范数的总变分正则化、截断最小二乘法等，这些方法以心脏跨膜电位的空间平滑特性作为约束，得到了一定精度的解；但它们大都孤立的研究了某一时间点上整个心脏的电活动，而没有考虑其时间相关性，因此精度不高，且求解一个心动周期内心脏的兴奋传递情况所需的时间随序列长度线性增加，十分耗时。

因此，选取合适的先验条件作为约束，将电位分布的时间相关性与空间特性结合起来，以提高解的精度，降低计算复杂度，是一个十分值得研究的内容。

发明内容

为了解决心脏电生理的反演问题，本发明提供了一种基于图总变分约束的无创心脏电生理反演方法，以心脏节点上的跨膜电位序列分布的相关性构建图结构，增强相同状态心脏节点间的相似性，区分不同状态的节点。

一种基于图总变分约束的无创心脏电生理反演方法，包括如下步骤：

(1)采集病人的64导联体表电位序列以及胸腔CT扫描图像；

(2)根据胸腔CT扫描图像分别建立躯干和心脏的三维网格模型，并在同一坐标系下进行配准；

(3)结合心肌各向异性传导信息，得到心肌跨膜电位与体表电位间的映射关系Φ＝HU，Φ为体表电位序列，U为心肌跨膜电位序列，H为两者间的转换矩阵；

(4)根据上述映射关系，在已知体表电位序列Φ的基础上利用二阶吉洪诺夫正则化求解心肌跨膜电位序列关于电生理反演问题低精度的初始解

(5)根据初始解

并考虑各心脏节点状态的相似程度，建立图结构G＝(V,E,W)，其中V为图结构G的顶点集即所有心脏节点的集合，E为图结构G的边集合即记录了心脏节点间的相似程度，W为权重矩阵即为每条边赋予不同的权重值；

(6)根据上述图结构G根据下式计算心肌跨膜电位序列U的图总变分

其中：i为集合V中的任一心脏节点，N_i为心脏节点i所有邻接节点的集合，j为集合N_i中的任一心脏节点，w(i,j)为连接心脏节点i和j的边的权重值，u_i和u_j分别为心脏节点i和j上的跨膜电位序列,|| ||₁表示1范数；

(7)根据图总变分

建立以下目标函数，并对其进行优化求解得到高精度的心肌跨膜电位分布序列U；

其中：μ为正则化参数(用以权衡数据保真项与图总变分正则化项),|| ||_F表示F范数。

进一步地，所述步骤(1)中使病人穿戴64导联电极设备进行CT扫描，扫描范围需包含所有导联的位置。

进一步地，所述步骤(1)中得到病人64导联体表电位序列后需对其进行去噪以及拉平预处理。

进一步地，所述步骤(2)中三维展示包含64导联位置的胸腔CT扫描图像，标记64个导联电极在躯干上的位置，进而通过德洛内三角剖分法构建躯干三维网格模型。

进一步地，所述步骤(2)中截取心脏短轴方向的约20张CT切片图像，范围从右心室流出道高度至心尖，在每张CT切片图像上标记左心室、右心室与心外膜的轮廓，进而通过德洛内三角剖分法构建心脏三维网格模型。

进一步地，所述步骤(3)中的转换矩阵H根据bidomain假设并结合心肌的各向异性传导信息及特性通过计算得到。

进一步地，所述步骤(4)中通过以下算式求解心肌跨膜电位序列U关于电生理反演问题的低精度解；

其中：L为拉普拉斯矩阵，λ为正则化系数，^T表示转置。

进一步地，所述步骤(5)的具体实现方式如下：

5.1以各心脏节点作为图中顶点，并利用初始解

为各顶点赋予初值即各心脏节点对应的初始跨膜电位序列，从而得到顶点集V；

5.2对于集合V中的任一心脏节点i，根据初始跨膜电位序列计算心脏节点i与其他所有心脏节点之间的欧氏距离；

5.3根据欧氏距离利用K邻近搜索算法确定与心脏节点i最邻近的k个节点，组成集合N_i，k为大于1的自然数；

5.4将连接心脏节点i与最邻近k个节点的k条边保存至边集合E中；

5.5为不同的边赋予不同的权重值，从而得到权重矩阵W，其中：

其中：σ为图结构G中所有边的平均长度，|| ||₂表示2范数。

进一步地，所述步骤(7)中采用前后原始对偶法交替求解目标函数中可导的2范数数据保真项和不可导的1范数正则化项，直至收敛得到唯一的最优解。

本发明基于图总变分约束进行了无创心脏电生理反演，图结构的建立考虑了空间不相邻节点之间相似性，充分利用了心肌跨膜电位序列分布的底层特征，同时考虑了电兴奋传播的时间相关性，提供了高精度的心肌跨膜电位序列分布。本发明通过最终得到的心肌跨膜电位的时空分布，可以准确的观察到心肌梗死时的梗死区域位置及边界、异位起搏发生时的起搏点以及其它电活动异常的区域，在临床的诊断及治疗上具有重要的参考意义。本发明可直接同时重建激活的整个序列，而非按照单个时间点的求解。

附图说明

图1为本发明基于图总变分约束无创心脏电生理反演方法的流程示意图。

图2为不同噪声扰动下正则化参数μ的最优值测试结果序列。

图3(a)为心肌梗死时心肌跨膜电位分布的真值图。

图3(b)本发明重建心肌梗死时心肌跨膜电位分布结果图。

图3(c)为二阶吉洪诺夫方法重建心肌梗死时心肌跨膜电位分布结果图。

图4为真实室性早搏病人的第37导联心电图。

图5为本发明重建真实室性早搏病人心肌跨膜电位分布结果图。

图6为室性早搏病人腔内Ensite3000诊断图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，本发明基于图总变分约束的无创心脏电生理反演方法，包括如下步骤：

(1)病人穿戴64导联电极设备进行CT扫描，扫描范围需包含所有导联的位置，得到含有电极位置的胸腔CT图像，对得到的64导联体表电位序列进行去噪与拉平预处理，得到体表电位序列Φ。

(2)3D展示包含64导联位置的CT图像，标记64个导联电极在躯干上的位置，由德洛内三角剖分法则构建躯干三维网格模型。

截取心脏短轴方向的约20张CT切片，范围从右心室流出道高度到心尖，在每张切片上标记左心室，右心室与心外膜的轮廓，由德洛内三角剖分法则构建心脏三维网格模型，并在同一坐标系内将心脏和躯干模型配准。

根据bidomain假设，考虑到心肌的各向异性特性，计算得到转换矩阵H，结合各向异性的传导信息，得到心肌跨膜电位与体表电位间的映射关系：

Φ＝HU

其中：Φ为体表电位序列，U为心肌跨膜电位序列，H为两者间的转换矩阵。

(3)根据映射关系Φ＝HU，利用二阶吉洪诺夫正则化，得到已知体表电位序列Φ求解心肌跨膜电位序列U的电生理反演问题的低精度的解：

其中：L为拉普拉斯矩阵。

(4)根据二阶吉洪诺夫正则化得到的初始低精度的解，考虑各个心脏节点状态的相似程度，建立图结构G＝(V，E，W)，其中V为图G的顶点集，即所有心脏节点的集合；E为图G的边的集合，记录了各个节点间的相似程度；W为权重矩阵，为不同的边赋予不同的权重。

构建图结构G＝(V，E，W)的具体过程如下：

4.1根据二阶吉洪诺夫正则化得到的初始低精度的解，以各个节点为顶点，以各节点上心肌跨膜电位序列为该顶点的值，得到顶点集V。

4.2计算各个心脏节点两两之间的距离：

l＝||u_i-u_j||₂

其中：l为以节点i，j上的值为元素的两个向量u_i，u_j之间的欧式距离，两个节点越相似则距离越近。

4.3利用K邻近搜索算法寻找与每个节点最接近的k个节点，保存在该节点的邻接点集合中。

4.4连接各邻接点对之间的边保存在E中。

4.5为不同的边赋予不同的权重，得到权重矩阵W，保证相似程度越高的邻接点对间的边权重越大，按照热核权重的方式定义：

其中：σ为各边的平均长度。

(5)由图结构计算图总变分：

其中：n表示代表心脏的节点总数，n_i是节点i的所有邻接节点的集合，u_i，u_j为心脏节点i，j上的值，即该节点上的心机跨膜电位序列。

(6)由上述图总变分得到目标优化问题：

其中：μ为正则化参数，用以权衡数据保真项与图总变分正则化项，μ应随测量受噪声扰动的程度增大而增大，其测试结果如图2所示。

(7)用前后原始对偶法求解上述目标优化问题，令：

交替求解梯度下降和对偶问题：

(U_n+1,V_n+1)＝(U_n,V_n)+γ((P_n,Q_n)-(U_n,V_n))

迭代至算法收敛，得到最优解U_N，即为高精度的心肌跨膜电位分布序列。

在具体的实施过程中，本发明整个算法在Windows 10(64-bit)系统中测试，其中CPU为Inter(R)Core(TM)i7-8700(3.2GHz)，主机内存为16GB RAM，显卡型号为NVIDIAGeForce GTX 1070；在编程中，采用Matlab R2016a平台来进行编程。

为了验证本发明方法在识别心脏病变部位时的准确性与优越性，我们通过上述方法对仿真心肌梗死和真实室性早搏的病例进行实验。图3(a)～图3(c)分别显示了心肌梗死时心肌跨膜电位分布的真值、本发明方法以及二阶吉洪诺夫方法的重建结果；浅色部位为梗死部位，噪声扰动为20dB，本发明方法在准确勾勒梗死部位位置与边界上有明显优势。图4为真实室性早搏病人的心电图(第37导联)，图5为本发明方法重建的该病人心肌跨膜电位分布，箭头所指为异位起搏点，位于右心室流出道前间隔侧，与图6所示的病人腔内Ensite3000诊断图所示的异位起搏点一致。

上述的对具体实施方式的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述具体实施方式做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (1)

1.一种基于图总变分约束的无创心脏电生理反演方法，包括如下步骤：

(1)采集病人的64导联体表电位序列以及胸腔CT扫描图像，具体地：使病人穿戴64导联电极设备进行CT扫描，扫描范围需包含所有导联的位置，得到病人64导联体表电位序列后需对其进行去噪以及拉平预处理；

(2)根据胸腔CT扫描图像分别建立躯干和心脏的三维网格模型，并在同一坐标系下进行配准，具体地：三维展示包含64导联位置的胸腔CT扫描图像，标记64个导联电极在躯干上的位置，截取心脏短轴方向的约20张CT切片图像，范围从右心室流出道高度至心尖，在每张CT切片图像上标记左心室、右心室与心外膜的轮廓，进而通过德洛内三角剖分法构建心脏三维网格模型；

(3)结合心肌各向异性传导信息，得到心肌跨膜电位与体表电位间的映射关系Φ＝HU，Φ为体表电位序列，U为心肌跨膜电位序列，H为两者间的转换矩阵，其根据bidomain假设并结合心肌的各向异性传导信息及特性通过计算得到；

(4)根据上述映射关系，在已知体表电位序列Φ的基础上利用二阶吉洪诺夫正则化求解心肌跨膜电位序列关于电生理反演问题低精度的初始解

具体计算表达式如下：

其中：L为拉普拉斯矩阵，λ为正则化系数，^T表示转置；

(5)根据初始解

并考虑各心脏节点状态的相似程度，建立图结构G＝(V,E,W)，具体实现方式如下；其中V为图结构G的顶点集即所有心脏节点的集合，E为图结构G的边集合即记录了心脏节点间的相似程度，W为权重矩阵即为每条边赋予不同的权重值；

5.1以各心脏节点作为图中顶点，并利用初始解

为各顶点赋予初值即各心脏节点对应的初始跨膜电位序列，从而得到顶点集V；

5.2对于集合V中的任一心脏节点i，根据初始跨膜电位序列计算心脏节点i与其他所有心脏节点之间的欧氏距离；

5.3根据欧氏距离利用K邻近搜索算法确定与心脏节点i最邻近的k个节点，组成集合N_i，k为大于1的自然数；

5.4将连接心脏节点i与最邻近k个节点的k条边保存至边集合E中；

5.5为不同的边赋予不同的权重值，从而得到权重矩阵W，其中：

其中：σ为图结构G中所有边的平均长度，|| ||₂表示2范数；

(6)根据上述图结构G根据下式计算心肌跨膜电位序列U的图总变分

其中：i为集合V中的任一心脏节点，N_i为心脏节点i所有邻接节点的集合，j为集合N_i中的任一心脏节点，w(i,j)为连接心脏节点i和j的边的权重值，u_i和u_j分别为心脏节点i和j上的跨膜电位序列,|| ||₁表示1范数；

(7)根据图总变分

建立以下目标函数，并对其进行优化求解得到高精度的心肌跨膜电位分布序列U，具体则采用前后原始对偶法交替求解目标函数中可导的2范数数据保真项和不可导的1范数正则化项，直至收敛得到唯一的最优解；

其中：μ为正则化参数,|| ||_F表示F范数。

Images