一种基于自编码器的心肌缺血位置的预测方法
【专利摘要】本发明公开了一种基于自编码器的心肌缺血位置的预测方法，该预测方法结合深度学习，对心肌缺血的运动场信息进行稀疏编码，能够对心肌缺血位置进行准确的预测；与临床上现有的心肌缺血诊断手段相比，本发明不需要给病人注射造影剂进行延迟增强造影，而是通过普通磁共振电影序列数据就可以对心肌缺血的位置进行准确的预测。
【专利说明】
一种基于自编码器的心肌缺血位置的预测方法
技术领域
[0001] 本发明属于心脏运动分析技术领域，具体涉及一种基于自编码器的心肌缺血位置 的预测方法。
【背景技术】
[0002] 心脏可以看作一个弹性体，正常的组织和病变组织往往具有不同的运动特征，因 此对心肌运动进行定量分析，对诊断心肌缺血具有重要意义。心脏发生形变的程度可以用 弹性力学中的应变来描述，根据应变的定义，需要先求出心肌的位移场，再计算其偏导数。 因此对心肌应变场的分析就转化为计算心肌的运动场，而且在计算过程中需要注意对噪声 的抑制，因为位移场的噪声在求取应变的时候会由于导数而放大，这对成像和计算的精确 度有了更高的要求。
[0003] 对心脏运动的定量分析有很多种方法，研究人员在早期把许多自然图像处理的方 法和理论引入心脏运动的定量分析，包括光流场(Optical Flow)以及通过基于曲率的匹配 和基于标记的方法来寻找图像序列中的观测点之间对应关系等。追踪人工植入标记的方法 可以作为分析心脏运动的金标准，但是这样侵入式的方法可能影响心脏原有的运动规律， 也有可能给病人带来额外的健康风险。结合现代医学影像模式，比如核磁共振标记等方法 在分析心脏运动领域可以克服前面提出的若干问题，具有安全、准确等优点。
[0004] 然而目前临床上诊断心肌缺血以核磁共振延迟对比增强（Delay Contrast Enhancement)影像作为金标准，为了获得延迟对比增强影像需要病人注射造影剂，这样的 做法增加了诊断心肌缺血的成本和安全风险。因此，如何从现有的核磁共振电影序列出发， 对心肌缺血的特征进行表达，进一步给出有效的诊断和缺血位置预测成为一个十分有意义 的研究问题。
【发明内容】
[0005] 针对现有技术所存在的上述技术问题，本发明提供了一种基于自编码器的心肌缺 血位置的预测方法，通过深度神经网络学习表达心肌缺血的运动特征信息，最终以磁共振 电影序列数据为依据，对心肌缺血的位置进行有效的预测。
[0006] -种基于自编码器的心肌缺血位置的预测方法，包括如下步骤：
[0007] (1)采集同一心动周期内关于心肌缺血的多帧磁共振电影序列影像作为一组训练 样本，依此对连续多个心动周期进行采集对应得到多组训练样本；
[0008] (2)对所述的训练样本进行特征提取，以得到左心内膜边界的运动场信息；
[0009] (3)对于任一组训练样本，根据其对应的磁共振延迟对比增强影像通过人工判断 对该影像中左心室缺血位置进行标记，得到该训练样本对应心肌缺血位置的标签图像；
[0010] (4)构建由多个自编码器累加而成的神经网络模型，进而使训练样本对应的运动 场信息作为输入量，对应的标签图像作为输出量对该神经网络模型进行训练，并最终确立 得到心肌缺血位置的预测模型；
[0011] (5)根据步骤（1)采集同一心动周期内关于心肌缺血的多帧磁共振电影序列影像 作为一组测试数据;然后根据步骤(2)对测试数据进行特征提取，得到测试数据关于左心内 膜边界的运动场信息；最后将测试数据的运动场信息输入上述预测模型中，从而输出得到 对应关于心肌缺血位置的标签图像。
[0012] 所述的步骤(2)中在特征提取之前，先对训练样本中各帧磁共振电影序列影像依 次进行高斯平滑和各向异性扩散滤波的预处理;能够实现对训练集数据中噪声的抑制，各 向异性扩散滤波在去除噪声的同时还能够较好地保持图像的边缘和重要细节信息。
[0013] 所述的步骤(2)中对训练样本进行特征提取的过程依次包括心内膜边界提取和轨 迹追踪测量。
[0014] 所述心内膜边界提取的过程为:首先，对于训练样本中任一帧磁共振电影序列影 像，通过人工在该影像上粗略地勾勒出左心室的心内膜边界曲线，以此作为初始状态;然 后，采用梯度矢量流计算图像灰度在两个正交方向上的梯度作为外力场，计算心内膜边界 曲线的弹性能量和弯曲能量作为内力场;最后，利用Snake算法使心内膜边界曲线在外力场 和内力场的作用下达到平衡状态，平衡状态下左心室的心内膜边界曲线即作为该帧磁共振 电影序列影像的心内膜边界提取结果。
[0015] 所述轨迹追踪测量的方法为:在完成对训练样本的心内膜边界提取后，采用基于 形状特征的追踪算法，确定同一心动周期内左心室的心内膜边界曲线上每一边缘像素点在 各帧磁共振电影序列影像中的对应关系，从而得到左心内膜边界的运动场信息；左心内膜 边界的运动场信息由左心室的心内膜边界曲线上所有边缘像素点的运动轨迹组成，任一边 缘像素点的运动轨迹即为该点在各帧磁共振电影序列影像中的坐标信息所组成的列向量。
[0016] 所述的自编码器由输入层、隐藏层和输出层组成;其中，前一个自编码器的隐藏层 为后一个自编码器的输入层;对于任一自编码器，其隐藏层的神经元个数比输入层的神经 元个数少。
[0017] 所述自编码器的函数模型如下：
[0018] h = 〇(wt+b)
[0019] z = 〇(w，h+b，)
[0020] 其中：t、h和z分别为自编码器的输入层、隐藏层和输出层，w和b均为输入层与隐藏 层之间的模型参数，w'和b'均为隐藏层与输出层之间的模型参数， 〇(s)为神经元函数且
,s为神经元函数0(s)的自变量。
[0021 ]所述的步骤(4)中对神经网络模型进行训练的具体方法如下：
[0022] 对于神经网络模型中的第一个自编码器，以训练样本的输入量作为该自编码器的 输入层，使该自编码器输出层与输入层的损失函数L最小为目标，通过梯度下降法求解出该 自编码器输入层与隐藏层之间以及隐藏层与输出层之间的模型参数；
[0023] 对于神经网络模型中除第一个和最后一个以外的任一自编码器，以前一个自编码 器的隐藏层作为该自编码器的输入层，使该自编码器输出层与输入层的损失函数L最小为 目标，通过梯度下降法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的 模型参数；
[0024] 对于神经网络模型中的最后一个自编码器，以前一个自编码器的隐藏层作为该自 编码器的输入层，使训练样本的输出量与该自编码器输入层的损失函数L'最小为目标，通 过反向传播法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的模型参 数。
[0025]所述的损失函数L和L'的表达式如下：
[0028]其中：x为训练样本的输出量，N为训练样本的个数，KL(h| |p)表示h与P之间的相对 熵即KL散度，P为对应h的期望平均激活函数，λ和β均为权重系数。
[0029] 本发明心肌缺血位置的预测方法结合深度学习，对心肌缺血的运动场信息进行稀 疏编码，能够对心肌缺血位置进行准确的预测;与临床上现有的心肌缺血诊断手段相比，本 发明不需要给病人注射造影剂进行延迟增强造影，而是通过普通磁共振电影序列数据就可 以对心肌缺血的位置进行准确的预测。
【附图说明】
[0030] 图1为本发明基于自编码器神经网络模型预测心肌缺血位置的框架示意图。
[0031] 图2(a)为心脏短轴方向的磁共振延迟对比增强影像。
[0032] 图2(b)为左心室在分成4扇区时预测的缺血位置的标签图像。
[0033] 图2(c)为左心室在分成8扇区时预测的缺血位置的标签图像。
[0034]图2(d)为左心室在分成12扇区时预测的缺血位置的标签图像。
【具体实施方式】
[0035] 为了更为明确地描述本发明，下面结合附图及【具体实施方式】对本发明的技术方案 进行详细说明。
[0036] 本发明基于自编码器的心肌缺血位置的预测方法，具体实施步骤如下：
[0037] S1.获取训练样本。
[0038] 采集心肌缺血病人的磁共振电影序列数据和磁共振延迟对比增强数据，要求同一 心动周期内采集不少于5帧图像，并以若干个心动周期的上述影像作为训练集，其中磁共振 电影序列数据用于提供心脏运动场信息，延迟对比增强数据用于提供训练所需的心肌缺血 位置标签。
[0039] S2.对训练样本进行滤波。
[0040] 对训练集中每个心动周期的数据进行预处理，包括高斯平滑和各向异性扩散滤 波，实现对训练集数据中噪声的抑制，各向异性扩散滤波在去除噪声的同时还能够较好地 保持图像的边缘和重要细节信息。
[0041] S3.提取心内膜边界。
[0042] 采用梯度矢量流(Gradient vector flow)计算边缘图的外力场，利用Snake算法 对左心室进行分割，从而实现了心内膜边界提取。具体地:对于心脏磁共振电影序列的某一 帧短轴图像，先通过人工标记的方法粗略地画出左心室的心内膜边界曲线，以此作为初始 状态;然后计算图像灰度在两个正交方向上的梯度作为外力场(也称为GVF场），计算边界曲 线的弹性能量和弯曲能量作为内力场，最终心内膜边界曲线在外力场和内力场的作用下达 到平衡状态，停留在图像灰度梯度值较大的地方，同时保持曲线的连续性和平滑性。
[0043] S4.获取心肌的运动场信息和缺血位置标签。
[0044] 采用基于形状特征的追踪算法，确定了同一心动周期内相邻帧图像上边界点之间 的对应关系，从而实现了心肌运动轨迹的追踪测量。具体地:将一个心动周期里心内膜边界 点的运动场写成一个列向量，再加上相应的标示心肌缺血位置的标签，便组成了一个训练 样本，追踪多个心动周期的心肌运动轨迹，进而组建得到多组训练样本;每组训练样本包括 左心内膜边界的运动场信息和相应的缺血位置标签，运动场信息包括通过估计得到的左心 室心内膜轮廓上若干点在同一心动周期内的对应关系，缺血位置标签P为通过相应的磁共 振延迟对比增强影像判断得到，其中P=1表示延迟对比增强影像中该位置出现心肌缺血症 状，P = 0表示影像中该位置没有明显的缺血症状。
[0045] S5.构建基于自编码器的神经网络模型。
[0046] 构建由多个自编码器叠加而成的深度神经网络模型，如图1所示；自编码器基本结 构由输入层、隐藏层和输出层组成;前一个自编码器的隐藏层可以作为下一个自编码器的 输入层;而且对于任意一个自编码器，其隐藏层的节点个数比输入层的节点个数少。选取若 干种不同的隐藏层节点个数组合，设计测试实验确定隐藏层节点数，需要综合考虑预测精 确度和程序运行时间两方面的因素。自编码器的函数模型如下：
[0047] h = 〇(wt+b)
[0048] z = 〇(w，h+b，)
[0049] 其中：t、h和z分别为自编码器的输入层、隐藏层和输出层，w和b均为输入层与隐藏 层之间的模型参数，w'和b'均为隐藏层与输出层之间的模型参数， 〇(s)为神经元函数且
,s为神经元函数0(s)的自变量。
[0050] S6.将心肌运动场信息和相应标签用于训练由堆叠自动编码器构成的深度神经网 络结构，得到每一层的结构参数，包括权重矩阵W和偏置b等。
[0051 ]首先将训练样本的输入量作为该自编码器的输入层，以该自编码器输出层与输入 层的损失函数L最小为目标，通过梯度下降法求解出该自编码器输入层与隐藏层之间以及 隐藏层与输出层之间的结构参数；
[0052] 对于神经网络模型中除第一个和最后一个以外的任一自编码器，以前一个自编码 器的隐藏层作为该自编码器的输入层，使该自编码器输出层与输入层的损失函数L最小为 目标，通过梯度下降法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的 结构参数；
[0053] 对于神经网络模型中的最后一个自编码器，以前一个自编码器的隐藏层作为该自 编码器的输入层，使训练样本的输出量与该自编码器输入层的损失函数L最小为目标，通过 反向传播法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的模型参数。
[0054] 损失函数L定义如下：
[0056]其中：pj表示第1层中第j个单位的平均激活函数，P是对应单元的期望平均激活函 数，N表示训练集中样本的个数。最后一项为KL散度，定义如下，
[0058] 损失函数中的第一项是输出值和输入值的均方根误差，在重建过程中要使其尽可 能地小;第二项是正则项，使得权重矩阵极可能接近零;第三项起稀疏约束作用，最终隐藏 层中大部分单位激活值为〇,将不起作用，从而确保输入层与隐藏层、隐藏层和输出层之间 的稀疏连接。
[0059] S7.将测试数据输入模型，得到心肌缺血位置标签。
[0060] 根据步骤S1采集同一心动周期内的多帧磁共振电影序列数据作为测试数据;然后 根据步骤S2对测试集数据进行预处理，根据步骤S3和S4通过追踪同一个心动周期的心肌运 动轨迹得到心肌的运动场信息；最后将心肌的运动场信息输入上述神经网络模型中，从而 输出得到预测心肌缺血位置的标签图像。
[0061] 将预测得到的心肌缺血位置标签图像与磁共振延迟对比增强影像真值对比如图2 所示，可见采用本发明预测得到的缺血位置与真实缺血位置非常接近。
[0062]以下我们通过实验来验证本发明的可靠性，实验的运行环境为：8G内存，CHJ为 intel i5,主频3.47GHz。结合表1的测试实验结果，在该实验中采用两个隐含层，第一个隐 含层有100个节点，第二个隐含层有30个节点，该节点数组合在预测缺血位置的准确率方面 优于其它9种组合方案，并且在训练时间上也可以接受。需要注意的是，在测试实验中采用 三个隐含层时，也可以达到相同的最高准确率，但是训练时间显著增长，因此在同等条件 下，我们选择设置两个隐含层，并将该组合方案应用到后续所有的实验中。
[0063] 表 1
[0064]
[0066]通过前边的测试实验，可以确定深度网络的隐含层节点数和隐含层数量，然后将 该网络在经过数据预处理后的训练集上训练，先后经过预训练和基于反向传播算法的调优 两个过程，最终确定深度网络的其它结构参数比如权重系数和偏置等。
[0067] 为了排除数据特异性对深度神经网络参数训练的影响，更加精确地反映所训练的 深度神经网络性能，实验中采取了五折交叉验证(five-fold cross-validation)框架。在 五折交叉验证框架下，我们将得到的数据集随机地分成5份，然后每次将其中一份数据用作 测试集，另外的四份数据用于训练深度神经网络得到心肌缺血位置预测模型。这一过程先 后重复5次，从而保证每一部分数据都被测试过，最后将预测的平均准确率作为衡量预测模 型精确度的指标。表2显示了训练得到的深度神经网络在五折交叉验证框架下，先后5次实 验中将左心室划分为不同数量的扇区时对心肌缺血位置的预测准确率。
[0068] 表 2
[0069]
[0070] 上述的对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发 明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的 一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例， 本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护 范围之内。
【主权项】
1. 一种基于自编码器的心肌缺血位置的预测方法，包括如下步骤： (1) 采集同一心动周期内关于心肌缺血的多帧磁共振电影序列影像作为一组训练样 本，依此对连续多个心动周期进行采集对应得到多组训练样本； (2) 对所述的训练样本进行特征提取，以得到左心内膜边界的运动场信息； (3) 对于任一组训练样本，根据其对应的磁共振延迟对比增强影像通过人工判断对该 影像中左心室缺血位置进行标记，得到该训练样本对应心肌缺血位置的标签图像； (4) 构建由多个自编码器累加而成的神经网络模型，进而使训练样本对应的运动场信 息作为输入量，对应的标签图像作为输出量对该神经网络模型进行训练，并最终确立得到 心肌缺血位置的预测模型； (5) 根据步骤（1)采集同一心动周期内关于心肌缺血的多帧磁共振电影序列影像作为 一组测试数据;然后根据步骤(2)对测试数据进行特征提取，得到测试数据关于左心内膜边 界的运动场信息；最后将测试数据的运动场信息输入上述预测模型中，从而输出得到对应 关于心肌缺血位置的标签图像。2. 根据权利要求1所述的预测方法，其特征在于:所述的步骤(2)中在特征提取之前，先 对训练样本中各帧磁共振电影序列影像依次进行高斯平滑和各向异性扩散滤波的预处理。3. 根据权利要求1所述的预测方法，其特征在于:所述的步骤(2)中对训练样本进行特 征提取的过程依次包括心内膜边界提取和轨迹追踪测量。4. 根据权利要求3所述的预测方法，其特征在于:所述心内膜边界提取的过程为:首先， 对于训练样本中任一帧磁共振电影序列影像，通过人工在该影像上粗略地勾勒出左心室的 心内膜边界曲线，以此作为初始状态;然后，采用梯度矢量流计算图像灰度在两个正交方向 上的梯度作为外力场，计算心内膜边界曲线的弹性能量和弯曲能量作为内力场;最后，利用 Snake算法使心内膜边界曲线在外力场和内力场的作用下达到平衡状态，平衡状态下左心 室的心内膜边界曲线即作为该帧磁共振电影序列影像的心内膜边界提取结果。5. 根据权利要求4所述的预测方法，其特征在于:所述轨迹追踪测量的方法为:在完成 对训练样本的心内膜边界提取后，采用基于形状特征的追踪算法，确定同一心动周期内左 心室的心内膜边界曲线上每一边缘像素点在各帧磁共振电影序列影像中的对应关系，从而 得到左心内膜边界的运动场信息;左心内膜边界的运动场信息由左心室的心内膜边界曲线 上所有边缘像素点的运动轨迹组成，任一边缘像素点的运动轨迹即为该点在各帧磁共振电 影序列影像中的坐标信息所组成的列向量。6. 根据权利要求1所述的预测方法，其特征在于:所述的自编码器由输入层、隐藏层和 输出层组成;其中，前一个自编码器的隐藏层为后一个自编码器的输入层;对于任一自编码 器，其隐藏层的神经元个数比输入层的神经元个数少。7. 根据权利要求6所述的预测方法，其特征在于:所述自编码器的函数模型如下： h = 〇(wt+b) z = 〇(w，h+b，） 其中：t、h和Z分别为自编码器的输入层、隐藏层和输出层，W和b均为输入层与隐藏层之 间的模型参数，W'和b'均为隐藏层与输出层之间的模型参数，〇(s)为神经元函数且 CF(ty)^ 1 ，S为神经元函数〇(s)的自变量。 w 1 + e8. 根据权利要求7所述的预测方法，其特征在于:所述的步骤(4)中对神经网络模型进 行训练的具体方法如下： 对于神经网络模型中的第一个自编码器，以训练样本的输入量作为该自编码器的输入 层，使该自编码器输出层与输入层的损失函数L最小为目标，通过梯度下降法求解出该自编 码器输入层与隐藏层之间以及隐藏层与输出层之间的模型参数； 对于神经网络模型中除第一个和最后一个以外的任一自编码器，以前一个自编码器的 隐藏层作为该自编码器的输入层，使该自编码器输出层与输入层的损失函数L最小为目标， 通过梯度下降法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的模型 参数； 对于神经网络模型中的最后一个自编码器，以前一个自编码器的隐藏层作为该自编码 器的输入层，使训练样本的输出量与该自编码器输入层的损失函数L'最小为目标，通过反 向传播法求解出该自编码器输入层与隐藏层之间以及隐藏层与输出层之间的模型参数。9. 根据权利要求8所述的预测方法，其特征在于:所述的损失函数L和L '的表达式如下：其中：X为训练样本的输出量，N为训练样本的个数，KL(h | | P)表示h与P之间的相对熵即 KL散度，P为对应h的期望平均激活函数，λ和β均为权重系数。
【文档编号】G06F19/00GK105868572SQ201610257836
【公开日】2016年8月17日
【申请日】2016年4月22日
【发明人】刘华锋, 马骁勇
【申请人】浙江大学