基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法

1.本发明属于无创心脏电生理成像技术领域，具体涉及一种基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法。


背景技术：

2.目前社会上因心律失常引发的心血管疾病已经成为影响全世界人类健康甚至导致死亡的主要疾病之一，根据美国心脏协会和世界卫生组织统计，每年因心律失常造成的心脏疾病死亡的人数接近千万，而医疗水平如此发达的美国每年依然有接近38万死于冠心病，在美国每年死亡人数中的占比达到了惊人的17％左右；而根据中国心血管病报告中的内容显示，我国的形势同样严峻，现存的心血管疾病患者的人数数以亿计且其发病率和死亡率依旧在增加，目前已经超过癌症在内的其它所有疾病导致的死亡率。因此，发展心脏疾病的诊断技术具有重要意义。
3.由于心脏上的电位比较难以直接测量，因此临床上往往通过心电图这类方式对心脏异常情况做简单的初步诊断。医生通过对不同导联心电图的特征信息观察以及根据多年的临床经验可以大致推断出心脏疾病类型和模糊的病灶点位置。但是这仅仅是一个简单的定性评估，而在临床术中诊断治疗时，某些频发性的心律失常需要非常精确的病灶点位置信息，而仅根据体表心电图的判断显然不能完全满足临床上的需求。而侵入式的电位标测系统需要侵入病人心脏测量心脏表面电位，寻找确定需要消融的点，这会对身体造成损伤，有一定程度的风险；不仅如此，临床上借助这类系统往往也只是局部测量以及部分心脏重建，这不利于心脏电信号的完整呈现。
4.因此基于体表心电信号的深度信息，并结合病人的生理信息，简单且无创地获得心脏电位信息具有非常重要的临床意义。


技术实现要素：

5.为了解决心肌跨膜电位的无创重建问题，提出了一种基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法，本发明是采用以下技术方案来实现的：
6.本发明公开了一种基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法，具体步骤如下：
7.1)获得病人的体表高导联电极位置和躯干增强ct的平扫图像，所述的高导联数设置为256导联、120导联、64导联中的任意一种；
8.2)根据获得的躯干ct图像的胸腔部分平扫图像，得到病人的心脏-躯干几何模型；
9.3)基于上述所得的心脏-躯干几何模型，并根据电信号在心肌上的传导过程，得到体表高导联电位和心肌跨膜电位的映射关系φ＝hu，φ表示体表电位，u表示心肌跨膜电位，h表示传递矩阵；
10.4)在病人标记的体表电极位置获得高导联心电信号，并对信号进行滤波去噪，基
线校准，得到预处理后的体表电位信号；
11.5)根据步骤3)中得到的映射关系，基于传递矩阵h和步骤4)中预处理后的体表电位信号，心肌跨膜电位进行初始化得到低精度的初始解；
12.6)构建基于全局特征的快速迭代软阈值收缩算法网络，其包括梯度下降步骤、软阈值收缩步骤、全局特征提取步骤和动量积累步骤，基于步骤5)中初始化得到低精度的初始解，通过该网络的迭代重建得到高精度的心肌跨膜电位序列。
13.作为进一步地改进，本发明所述的步骤2)中：通过matlab软件对步骤1)中采集得到的多张躯干ct图像的胸腔部分中进行左右心室和心外膜轮廓的标记，结合3d slicer软件标记的体表电极位置，进一步得到心脏-躯干有限元网格模型。
14.作为进一步地改进，本发明所述的步骤3)中：通过边界元方法计算得到传递矩阵h。
15.作为进一步地改进，本发明所述的步骤5)中：初始化得到低精度的初始解的初始化方式为：
16.u0＝q
init
φ
[0017][0018]
其中，u
train
和φ
train
表示训练集中心肌跨膜电位和体表64导联心电数据，b表示正则化参数，i表示单位矩阵。
[0019]
作为进一步地改进，本发明所述的步骤6)中：心肌跨膜电位重建问题转化成以下的最小化优化问题：
[0020][0021]
其中，u是需要重建的心肌跨膜电位，λ表示正则化参数，ψ表示任意形式的转换算子，是graph算子，gradient算子，nonlocal算子中的任意一种，l1范数表示稀疏先验信息的约束形式；
[0022]
通过进一步迭代软阈值收缩算法的推导，该优化问题转化为下面两个子问题：
[0023][0024][0025]
通过动量积累步骤对两个子问题迭代过程进行加速，进一步讲该问题转化为：
[0026]rk
＝y
k-α
kht
(hy
k-φ)
[0027][0028]
[0029][0030]
其中，其中，tk是计算过程引入的中间变量，yk表示过去两次迭代的重建值u
k-1
和u
k-2
的线性组合，同时也是下一次迭代开始的输入。
[0031]
作为进一步地改进，本发明所述的步骤6)中：基于全局特征的快速迭代软阈值收缩算法网络的具体流程包括如下四个步骤：
[0032]
6.1梯度下降步骤：该步骤重复出现在网络中的每一个阶段中，对应的是传统fista中的梯度下降步骤：
[0033]rk
＝y
k-α
kht
(hy
k-φ)；
[0034]
其中yk表示上一个阶段的输出，这个步骤保留了由心脏-躯干模型得到的转换矩阵h，即保留了心脏相关的生理信息，对网络收敛起到了很关键的作用；
[0035]
6.2软阈值收缩步骤：该步骤中利用神经网络结构来替代从rk到uk的过程，该步骤的计算过程如下式：
[0036][0037]
其中，soft()表示软阈值收缩算子，γk表示正则化参数，和表示非线性转换算子，在该步骤中利用conv-bn-relu-conv的神经网络结构单元来替代算子而算子则是算子的正交对称形式；
[0038]
并在网络中引入残差结构，可以将其转化为残差的形式：
[0039][0040]
uk＝rk resk；
[0041]
其中，resk表示uk和rk之间的差值，即软阈值收缩步骤所学到的信息；
[0042]
6.3全局特征提取步骤：该步骤充分利用心肌跨膜电位的非局部先验信息并利用简单的网络结构更好获取特征图的全局信息，该步骤的计算过程为：
[0043]gk
＝ck(s(qk
t
)v) zk[0044][0045]
其中，和ck都是不同的卷积核大小为1*1的卷积层，s表示soft max激活函数，zk表示第k个阶段中经过了soft(fk(rk),γk)操作得到的高维特征图；
[0046]
加入全局特征提取步骤后，残差的计算进一步转化为：
[0047][0048]
6.4动量积累步骤：通过对前两步迭代的重建目标输出u
k-1
和u
k-2
的线性组合引入一个新的变量yk，作为下一次迭代的输入。我们直接将整体看成每个阶段中的一个可学习的参数，那么动量积累步骤在就可以写成简单的形式：
[0049]
y1＝u0[0050]yk 1
＝uk ρk(u
k-u
k-1
)。
[0051]
本发明的有益效果如下：
[0052]
本发明提供了一种基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法，同时结合传统的快速迭代软阈值收缩算法(fast iterative shrinkage-thresholding algorithm,fista)的收敛理论以及深度学习的灵活性泛化性优势，用神经网络结构来替代传统方法中需要人工不断调整的正则化参数和正则项形式，利用网络自动的学习心电数据中的普适性特征，并且保留了病人的生理信息和严谨的物理模型，算法的收敛性和重建的精度都有理论的支持。
[0053]
本发明中的梯度下降步骤保留了由心脏-躯干模型计算所得的转换矩阵，使得网络在训练过程中有了生理模型的约束从而更容易收敛；考虑到心肌跨膜电位的非局部相似性特征，本发明中利用注意力机制来提取全局特征，使得网络的感受野增大；本发明中利用神经网络来替代传统fista中的软阈值收缩步骤，从而替代掉传统方法中特定的正则化形式，同时省去了人工不断调整正则化参数的过程，增强了网络灵活性和泛化性；本发明在迭代网络中增加了一个可学习的动量积累模块来替代fista中的动量积累步骤，使得网络收敛速度更快。本发明针对心肌跨膜电位的各种特性做出的一系列创新点改进，使得最终能实现更好的重建精度和更快的重建速度。
附图说明
[0054]
图1为本发明心肌跨膜电位的重建流程图；
[0055]
图2为本发明基于全局特征的快速迭代软阈值收缩算法网络流程图；
[0056]
图3为本发明基于全局特征的快速迭代软阈值收缩算法网络中全局特征提取模块流程图；
[0057]
图4为本方法和现有tikhonov、tv方法的对于异位起搏情况中心脏不同节点的心肌跨膜电位重建结果对比图；
[0058]
图5为本方法和现有tikhonov、tv方法的对于心肌梗死情况中心脏不同节点的心肌跨膜电位重建结果对比图；
具体实施方式
[0059]
为了更为明确地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。
[0060]
如图1所示，本发明基于全局特征的快速迭代软阈值收缩算法网络的无创心肌跨膜电位重建方法，具体实施步骤如下：
[0061]
1)获得病人的体表高导联电极位置和躯干ct图像的胸腔部分平扫图像，所述的高导联数设置为256导联、120导联、64导联中的任意一种，具体实施过程中设置为64导联；使病人佩戴好高导联心电采集背心或者提前标记好高导联电极位置后，再进行躯干胸腔的增强ct扫描；
[0062]
2)根据获得的躯干ct图像的胸腔部分平扫图像，得到病人的心脏-躯干几何模型；通过matlab软件对步骤1)中采集得到的多张胸腔ct图像中进行左右心室和心外膜轮廓的标记，结合3d slicer软件标记的体表电极位置，进一步得到心脏-躯干有限元网格模型。
[0063]
3)基于上述所得的心脏-躯干几何模型，并根据电信号在心肌上的传导过程，得到
体表64导联电位和心肌跨膜电位的映射关系φ＝hu，φ表示体表电位，u表示心肌跨膜电位，h表示传递矩阵；基于边界元方法，计算得到传递矩阵h。
[0064]
4)在病人标记的体表电极位置获得64导联心电信号，并对信号进行滤波去噪，基线校准，得到预处理后的体表电位信号；通过高导联心电采集设备得到病人体表心电信号，然后对采集得到的心电信号进行滤波去噪，基线校准的预处理。
[0065]
5)根据步骤3)中得到的映射关系，基于传递矩阵h和步骤4)中预处理后的体表电位φ，对心肌跨膜电位进行初始化得到低精度的初始解；得到心肌跨膜电位低精度解的初始化方式为：
[0066]
u0＝q
init
φ
[0067][0068]
其中，u
train
和φ
train
表示训练集中心肌跨膜电位和体表高导联心电数据，b表示正则化参数，i表示单位矩阵。
[0069]
6)构建基于全局特征的快速迭代软阈值收缩算法网络，如图2所示，其包括梯度下降步骤、软阈值收缩步骤、全局特征提取步骤和动量积累步骤，其中全局特征提取步骤基于改进的自注意力机制来捕捉通道维度上远距离的相似性，从而同时获取不同通道之间的全局特征信息，如图3所示。基于步骤5)中初始化得到低精度的初始解，通过该网络的迭代重建得到高精度的心肌跨膜电位序列。心肌跨膜电位重建问题可以转化成以下的最小化优化问题：
[0070][0071]
其中，u是需要重建的心肌跨膜电位，λ表示正则化参数，ψ表示任意形式的转换算子，可以是graph算子，gradient算子，nonlocal算子中的任意一种，l1范数表示稀疏先验信息的约束形式。
[0072]
通过进一步迭代软阈值收缩算法的推导，该优化问题可以转化为下面两个子问题：
[0073][0074][0075]
通过动量积累步骤对两个子问题迭代过程进行加速，进一步讲该问题转化为：
[0076]rk
＝y
k-α
kht
(hy
k-φ)
[0077][0078]
[0079][0080]
其中，其中，tk是计算过程引入的中间变量，yk表示过去两次迭代的重建值u
k-1
和u
k-2
的线性组合，同时也是下一次迭代开始的输入。
[0081]
基于全局特征的快速迭代软阈值收缩算法网络的具体步骤包括：
[0082]
6.1梯度下降步骤：该步骤重复出现在网络中的每一个阶段中，对应的是传统fista中的梯度下降步骤，
[0083]rk
＝y
k-α
kht
(hy
k-φ)
[0084]
其中yk表示上一个阶段的输出，这个步骤保留了由心脏-躯干模型得到的转换矩阵h，即保留了心脏相关的生理信息，对网络收敛起到了很关键的作用。
[0085]
6.2软阈值收缩步骤：该步骤中利用神经网络结构来替代从rk到uk的过程，该步骤的计算过程如下式：
[0086][0087]
其中，soft()表示软阈值收缩算子，γk表示正则化参数，和表示非线性转换算子，在该步骤中利用conv-bn-relu-conv的神经网络结构单元来替代算子而算子则是算子的正交对称形式。
[0088]
并在网络中引入残差结构，可以将其转化为残差的形式：
[0089][0090]
uk＝rk resk[0091]
其中，resk表示uk和rk之间的差值，即软阈值收缩步骤所学到的信息。
[0092]
6.3全局特征提取步骤：该步骤充分利用心肌跨膜电位的非局部先验信息并利用简单的网络结构更好获取特征图的全局信息，该步骤的计算过程为：
[0093]gk
＝ck(s(qk
t
)v) zk[0094][0095]
其中，和ck都是不同的卷积核大小为1*1的卷积层，s表示softmax激活函数，zk表示第k个阶段中经过了soft(fk(rk),γk)操作得到的高维特征图。
[0096]
加入全局特征提取步骤后，残差的计算进一步转化为：
[0097][0098]
6.4动量积累步骤：通过对前两步迭代的重建目标输出u
k-1
和u
k-2
的线性组合引入一个新的变量yk，作为下一次迭代的输入。直接将整体看成每个阶段中的一个可学习的参数，那么动量积累步骤在就可以写成简单的形式：
[0099]
y1＝u0[0100]yk 1
＝uk ρk(u
k-u
k-1
)
[0101]
结合上述提到的梯度下降步骤、软阈值收缩步骤、全局特征提取步骤以及动量积
累步骤四部分，本发明提出的基于全局特征的快速迭代软阈值收缩算法网络如图2所示。图中每一个虚线框的网络块对应传统快速迭代软阈值收缩算法的一次迭代，除了最后输出最终心肌跨膜电位重建值的阶段外，每个阶段中都包含上述四个步骤。输入的体表电位先经过初始化模块得到初始值，然后输入后面多阶段的网络。下方放大细节中的conv层的卷积核大小为1*3，除了首尾两个用做编解码的卷积层外，非线性算子中的卷积层输入和输出的特征维度都为64，第一个用作编码的卷积层的输入维度是1，输出维度是64；反之最后用作解码的卷积层的输出维度还原到1，即重建的目标心肌跨膜电位的中间解。另外在图中下方细节方法部分，加入了自身映射的残差结构，这样网络不是直接学习，而是学两者之间的差值，这样的设计使得网络更容易训练。理论推导过程中的梯度下降步长以及正则化参数在网络中都看成可学习的参数，不再根据情况人工选择调整。
[0102]
本发明所设计的基于全局特征的快速迭代软阈值收缩算法网络所考虑的损失函数主要考虑以下三个方面：一是考虑最后输出的心肌跨膜电位重建值与真值的相近程度，这一部分采用均方误差来描述；二是在理论推导中要求两个非线性算子正交，才能保证整体网络的收敛；三是考虑到心肌跨膜电位的稀疏性特性，希望非线性算子是将中间值映射到稀疏域然后再进行后续的算子处理，采用了l1范数来描述其稀疏性损失。最终的损失函数是由上述三部分损失的线性叠加：
[0103]
l
total
＝β1l
mse
β2l
sym
β3l
spa
[0104][0105][0106][0107]
其中，l
mse
表示第一部分均方误差的损失，l
sym
表示第二部分正交性损失，l
spa
表示第三部分稀疏性损失，l
total
是三部分的加权组合，表示最终的损失，k表示整个网络的阶段数，表示权重系数，β
*
这里分别设定为β1＝1，β2＝0.01，β3＝0.001。
[0108]
初始化后的低精度解通过训练好后的网络各步骤迭代重建出高精度解。
[0109]
在具体的实施过程中，本发明整个算法在ubuntu系统中测试，其中cpu为intel(r)core(tm)i9-7900x cpu@3.30ghz 3.31ghz,主机内存为32gb ram，显卡型号为nvidia titan rtx，主要编程环境为tensorflow1.14和python 3.7。
[0110]
为了验证本方法对心肌跨膜电位的重建效果，讲上述提出的方法应用于异位起搏数据和心肌梗死数据的心肌跨膜电位时间分布重建，图4显示了本方法和现有tikhonov、tv方法在异位起搏时心脏不同节点的重建结果图对比，心肌跨膜电位波形快速去极化(即快速上升)的时间点表示心肌细胞刚开始激活发生动作电位，图中三行从上到下分别表示激活时间从早到晚的顺序，如图中箭头所示其中平滑的线表示真值，有波动的线表示重建值。从图中可以看出，tikhonov方法虽然能重建出tmp波形的大致轮廓，但是全程夹杂了大量的噪声。tv方法在最开始激活时间点上重建不够准确，如node_423的情况，甚至在后期容易出现一些伪动作电位，可以明显看出本发明提出的方法得到的结果与真值更相近，精度更高。图5显示了本方法和现有tikhonov、tv方法在心肌梗死时心脏不同节点的重建结果图对比，
图中三行自上而下分别是梗死核中心、梗死疤痕中间区域和梗死疤痕边缘。对于梗死核中心，由于心肌细胞几乎无法被激活，所以动作电位的幅值非常低，心肌跨膜电位的波形几乎是平的，这种情况传统方法无法重建出比较好的轮廓，而本方法可以比较好的表现出梗死的程度。对于梗死疤痕中间和边缘的节点，在tikhonov方法和tv方法的重建结果中存在大量的噪声或者多峰情况(多个动作电位)，而本方法提出的基于模型的深度学习方法始终有比较好的表现，说明在时间维度上，本方法有很好的去噪效果和鲁棒性。
[0111]
表1为心肌跨膜电位重建消融实验中的各模型重建值与真值之间的相对误差(re)和结构相似度(ssim)的定量分析统计表，可以看出本方法提出的所有优化方案都对相对误差和结构相似性两种评价指标有不同程度的正面改进，其中全局特征提取步骤和动量积累步骤对心肌跨膜电位重建的作用都非常关键。
[0112]
表1心肌跨膜电位重建消融实验中的各模型重建值与真值之间的re和ssim的定量分析统计表
[0113][0114]
表2为tikhonov、tv、ista-net和gfista-net四种不同方法的心肌跨膜电位重建值与真值之间的相对误差(re)和结构相似度(ssim)的定量对比实验结果表，可以看出相较于传统方法，本方法大大提升了重建效果。
[0115]
表2 tikhonov、tv、ista-net和gfista-net四种不同方法的心肌跨膜电位重建值与真值之间的re和ssim的定量对比实验结果表
[0116]
上述的对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。