王甜甜 王慧 朱艳春 王丽嘉†

1) (上海理工大学医疗器械与食品学院,上海 200093)

2) (中国联通医疗基地,官洲生命科学创新中心,广州 510000)

心血管疾病(cardiovascular diseases,CVDs)的高发病率和高死亡率已经严重影响了人类的生存质量.如何评估心脏功能、辅助临床CVDs 诊疗和预后评估,是一个迫切需要解决的问题.针对这个问题,本文在前期心脏电影磁共振(cardiac cine magnetic resonance,CCMR)图像左心肌分割的基础上,提出一种基于位移流U-Net(DispFlow_UNet)和生物力学变分自动编码器(variational autoencoder,VAE)的左心肌运动追踪方法:DispFlow_UNet_VAE.主要研究内容有:1) 搭建压缩激励残差U-net 网络精准分割左心肌,根据分割结果计算心室体积、心肌质量等,评估心脏整体功能;2) 根据DispFlow_UNet_VAE 估计CCMR 图像连续帧之间的左心室运动,结合左心肌分割掩膜得到左心肌密集位移场;3)利用模拟数据真实位移场、临床数据集对追踪结果进行对比和评估.结果表明,本文追踪算法具有较高的精度和泛化能力.

1 引言

据世界卫生组织(world health organization,WHO)估计,每年全球有1790 万人死于心血管疾病(cardiovascular disease,CVDs)[1],如果按照目前的趋势继续下去,2030 年全球将有2360 万人死于CVDs[2].因此,如何准确评估心脏功能对于疾病的早诊断、早治疗具有重要临床意义.为了评估心脏功能,使用了两类指标:全局指标[3−5]和局部指标[6,7].全局指标包括射血分数、心肌质量、心室容积等,在CVDs 诊断中具有重要意义.然而,该类指标不能用于局部分析,无法反映细微的室壁运动异常,因此,在整体评价心脏功能的同时,需通过探索心肌运动规律,进一步分析心肌局部功能变化.

目前,借助医学成像手段评估心脏结构及功能指标是临床诊断CVDs 的主要依据.其中,临床常用的技术为心动超声图和CCMR 成像.超声心动图已广泛应用于应变分析,但其结果的准确性取决于操作技巧,并受到狭窄声窗的限制,在评估胸骨后的心脏时效果欠佳[8].CCMR 成像由于无电离辐射、软组织对比度高、可多参数和多方位成像等优点而被广泛应用于临床[9].其中,基于稳态自由进动(steady state free precession,SSFP)序列的CCMR 成像可动态、精确地描绘心肌、心室及其内部的复杂结构,是无创性量化心脏功能的金标准技术[10].一个典型的基于SSFP 序列的CCMR 图像如图1 所示,代表由舒张末期(end diastole,ED)至收缩末期(end systole,ES),再至ED 的一系列相位组成的心脏运动周期,其中心肌在很大程度上是同质的,不能直接提供运动模式的信息[11],因此,需要通过追踪算法来获取心肌的运动信息,从而得到运动位移场,以辅助医生早诊断早治疗CVDs.

图1 基于SSFP 序列的CCMR 图像Fig.1.CCMR images based on SSFP sequence.

由于心脏的收缩是非刚性,心脏追踪通常在像素级水平上进行.在CCMR 图像中,一个相位的运动通常是相对参考相位来估计的,设t 相位的图像为I(x,y,t),参考相位的图像为I(x,y,tref).运动追踪的目的是找到映射 Fθ[12]:

其中,Fθ是参数为θ 的映射函数,Vx和 Vy分别是沿x 和y 方向的位移矢量.

运动跟踪方法通常可以根据 Fθ的不同进行分类:基于强度的方法[13]、光流法[14]、基于心膜轮廓的算法[15]、基于非刚性配准算法[16−18]和基于深度学习(deep learning,DL)的方法[19−23]等.这些方法在运动估计中呈现出不同的计算负担和精确度,其中传统的追踪算法缺乏完全自动化,需要加入符合心肌运动模式的约束,限制了复现性,而基于DL 的追踪算法已被证明具有高效性、精确性和可重复性.

研究表明,心脏运动估计可以重新表述为一个可学习的问题.Qin 等[19]提出一种结合分割和运动估计的多任务框架,学习的心脏运动场用于扭曲分割掩膜,并以半监督方式引导分割模块,结果表明,分割和运动估计性能都有所提高.Zheng 等[20]提出一种表观流网络,它是一种改进的U 型网,在训练中使用分割掩膜来改进运动估计.然而,DL追踪算法使用常见的机器学习技术(如所有可学习参数的L2 正则化[20,22]、运动估计的直接正则化(如平滑惩罚[19]、解剖感知[23]))作为约束项,缺乏应用特定的先验知识来指导优化.因此,本文提出一种DispFlow_UNet_VAE 左心肌追踪算法,并使用基于生物力学的先验知识来约束学习,其可以隐式地学习正则化行为以指导优化,提高算法的高效性及复现性.

2 方 法

2.1 基于压缩激励残差U-net 的左心肌分割方法

精准的心肌分割是运动追踪的基础.Ronneberger 等[24]提出的U 形网络(U-Net)凭借所需训练集数量少、分割模糊边界精度高等优点而广泛应用于左心肌分割中,此后多种U-Net 改进算法应运而生,本文采用王慧[25]提出的压缩激励残差U-net (squeeze and excitation residual U-net,SERU-net)用于全心肌的自动分割,如图2 所示.该网络融合了压缩激励(squeeze-and-excitation,SE)模块和残差模块,其中SE 模块使网络在提取特征时能够通过学习自动获取每个特征通道的重要程度,更好提取左心肌有效特征;残差模块有效抑制了梯度消失和训练过拟合问题,使得信息前后向传播更加顺畅,使网络学习效果更好,分割精度更高.

图2 压缩激励残差U-net 网络结构Fig.2.Squeeze-and-excitation residual U-shaped network.

图3 给出了测试集中某一例数据从基底到顶端的分割结果,其中黄色为专家手动分割的金标准,红色为SERU-net 分割结果.从图3 可以看出,SERU-net 能够更好适应左心肌形状变化,精确勾画出左心肌.

图3 SERU-net 左心肌分割结果Fig.3.Results of left myocardium segmentation by SERU-net.

2.2 运动追踪算法

深度学习方法通常可分为有监督训练策略与无监督训练策略.在医学图像中,存在大量未标记数据和少量已标记数据,半监督学习被认为是充分运用可获得数据进行自动训练的有效方法[26].因此本文采取半监督深度学习策略.图4 给出了DispFlow_UNet_VAE 框架,其 由位移流U-Net (2.2.1 节)和VAE 正则化器(2.2.2 节)两部分组成.首先,VAE 用于学习模拟变形的概率分布以捕捉潜在的生物力学特征;接着,将学习好的VAE 用作位移流U-Net 配准网络的正则化函数,这样做的优点是,使解空间正则化并且有助于预测生物力学上合理的运动特征,而不需要任何显示的惩罚项.

图4 DispFlow_UNet_VAE 运动追踪框架Fig.4.The motion tracking architecture of DispFlow_UNet_VAE.

2.2.1 DispFlow_UNet 网络

如图5 所示,DispFlow_UNet 是U-net 的变体.在给定同一切片上的目标相位和参考相位作为输入的情况下,DispFlow_UNet 产生两者之间像素级的位移矢量场(displacement vector field,DVF),结合分割掩膜,从位移场中提取心肌运动特征.

图5 DispFlow_UNet 网络框架Fig.5.The network architecture of DispFlow_UNet.

由于无法获取图像之间的真实密集位移场,因此网络通过时间跟踪空间特征,依靠时间强度变化作为自我监督.将两个输入图像在位置P=(x,y)处的像素强度记为Is(P)和It(P),DispFlow_UNet产生s(源)与t(目标)相位之间的位移图记为Ft,Ft(P)=(Ftx(P),Fty(P)),通过空间变化器将Ft进行图像重建,使得以下强度差异最小化:

2.2.2 VAE 正则化器

为了进一步确保DispFlow_UNet 产生的位移场符合心肌运动实际,本文提出学习模拟变形概率分布的VAE 正则化网络,以隐含地学习生物力学上的形变规律,如图6 所示.VAE[27,28]被用作一个基于学习的正则化器,以模拟生物力学似是而非的变形的概率分布,其编码器用于将变形映射到潜在变量,通过高斯分布进行正则化,并通过解码器解码为变形场.VAE 正则化器被训练于重建生物力学模拟变形的一阶导数,以消除刚性平移的任何影响.用 Φ=[u,v]∈R2×H×W表示二维(2D)空间中的变形,其中,u,v 分别表示沿x 和y 方向的位移.变形场的一阶梯度表示为

图6 VAE 网络Fig.6.The network architecture of VAE.

其中H 和W 代表空间的维度.VAE 的重构损失公式为

其中,∇Φ′表示输入 ∇Φ 的重建项,z 代表由VAE编码的潜在矢量,p(z) N(0,I) 代表先验的高斯分布,qθ代表由θ 参数化的编码器,DKL代表KL 散度(Kullback-Leibler divergence),α 是控制重建质量和潜在空间规则程度之间折衷的超参数.

VAE 损失 RVAE提供一个定量的度量指标,用于确定变形在生物力学上的可信性.靠近已学习的VAE 潜在流形的解将产生较低的 RVAE,而远离流形的解将给出较高的 RVAE.

2.2.3 损失函数

密集追踪网络DispFlow_UNet_VAE 的目标损失函数由图像强度差异 LIMG(Ft) 和VAE 重构损失 RVAE(∇Φ ⊙M) 两部分组成:

其中,M 是分割掩膜的二值图,仅用于调整感兴趣区域,⊙ 代表哈达玛积(Hadamard product).γ 是超参数,权衡了图像的相似性和变形的物理合理性.

3 实 验

3.1 实验数据

该方法在Automated Cardiac Diagnosis Challenge (ACDC)数据集[29]上进行训练、验证和预测;在临床数据集上进行模型泛化能力评估;在模拟MRI 数据集[30]上进行对比实验.实验数据的具体情况如表1 所列.

表1 实验数据Table 1.The experimental data.

ACDC 数据集包括100 名受试者(健康和患病病例)的短轴CCMR 图像,使用两种不同主磁场的磁共振扫描仪(1.5T Siemens Area 和3.0T Siemens Trio Tim)在六年内通过SSFP 序列收集得到.所有数据采集都是屏气进行的,以确保在视频中只能观察到心脏运动.每个受试者包含9—10层图像,相位数在12—35 之间变化.总体来说,该数据集有951 个,每个序列提供ED 和ES 相位的金标准,如图7 所示.此外,100 名受试者被平均分为5 类,每类20 名受试者,分别为正常病例(normal subjects,NOR)、伴有埂塞的收缩性心力衰竭(myocardial infarction,MINF)、扩张性心肌病(dilated cardiomyopathy,DCM)、肥厚型心肌病(hypertrophic cardiomyopathy,HCM)和右心室异常(right ventricle abnormality,ARV).

图7 ED (左)与ES (右)的原始图像及金标准Fig.7.Original image and its ground truth of ED (left) and ES (right).

75 例临床数据集在 GE 1.5T 磁共振扫描仪通过 SSFP 序列获得,所有数据均符合伦理要求.采集参数是:图像大小为256×256、视野大小(field of vision,FOV)为360 mm×360 mm、层厚为6—8 mm、层间距为2—4 mm、重复时间(repetition time,TR)为 3.5 ms、回波时间(echo time,TE)为1.4 ms、层数为6—10 层、相位数为20—28 个.

3.2 实验流程

实验将从数据分配、数据预处理、训练与验证三步骤进行.

3.2.1 数据分配

在数据驱动的机器学习中,常假设训练样本和测试样本来自同一个分布(内部分布),然而测试集中的外部分布通常会给出较差的模型泛化.本文按照留下一种疾病的方法在ACDC 数据集上进行了五重交叉验证,根据已知的疾病将内部分布和外部分布分开.由于不同疾病之间存在显著的心脏解剖和动力学差异,因此与其他4 种疾病相比,一种疾病类别可被视为外部分布.

对于内部分布中的受试者,将其分为训练集(80*80%=72 名)、验证集(80*10%=8 名)和预测集(80*10%=8 名),外部分布中的所有受试者(20 名)被用作测试.在训练过程中,排除了覆盖LV 腔的所有切片的顶部20%和底部20%,而选择了中间剩余的60%,该设计旨在进一步减少平面外运动的影响.表2 给出了详细的数据分配情况.

表2 ACDC 数据集分配情况Table 2.ACDC data set allocation.

3.2.2 数据预处理

为了避免周围结构对左心肌的影响,需要在原始图像上确定感兴趣区域(region of interest,ROI),所有训练图像被裁剪为96×96 大小,同时将强度归一化到[0,1]的范围,并通过随机旋转、平移和缩放等操作扩充数据,提高网络精确性及鲁棒性.

3.2.3 训练与验证

网络搭建采用Pytorch 框架,操作系统为Windows 10.通过(3)式中的损失函数来训练网络参数,选取Adam 作为优化器,学习率=0.0001,批量大小为16,训练次数epoch=300,参数 γ=0.001,通过验证集选择得出.利用训练好的模型对预测集进行预测,根据预测结果进行评估.

3.3 评价方法

从三个方面评估本文追踪方法的精确性和泛化能力.

1)利用几何指标值量化本文算法的追踪精度.一般来说,手动制作心脏运动的参考标准是不可能的,为了进行定量分析,通常利用分割掩膜作为独立的参考标准.应用训练好的密集追踪网络去产生FES,使用FES来扭曲ES 帧的分割金标准MES,得到 MES°WFES,再利用Dice 系数、Hausdroff 距离(HD)和心肌轮廓距离(myocardial contour distance,MCD)测量其与ED 帧心膜掩膜图之间的重叠程度,Dice 和HD 的公式见(4)式和(5)式.

其中,A 是变形得到的ED 帧心肌轮廓,G 是ED帧掩膜,D(.,.)为欧几里得距离.Dice 系数取值越高、MCD 和Hausdroff 距离取值越低代表追踪精度越高.

2)利用模拟数据位移场验证本文方法的可行性.计算预测位移场与模拟数据位移场之间的平均角度误差(average angle error,AAE)和平均终止点误差(average endpoint error,AEPE),测量两位移场之间的重合程度.

设图像某像素点P 的角度误差AE(P)是预测位移矢量 u0(P)=(x0,y0) 和位移真值矢量u1(P)=(x1,y1)在点P 二维空间的角度,则:

AAE 为计算区域中AE(P)的均值,其取值越小表明追踪结果越准确.

图像中某像素点P 的终止点误差EPE(P)表示两个位移终止点 u0(P)=(x0,y0)、u1(P)=(x1,y1)之间的欧氏距离:

AEPE 为计算区域中EPE(P)的均值,其取值越小表明追踪结果越准确.

3)利用临床数据集(即外部数据)验证本文模型的泛化能力.使用训练好的模型直接预测临床数据集,通过对比模型在内部数据和外部数据集上的性能评估模型的泛化能力.

4 结果与讨论

4.1 全局功能评估

根据SERU-net 方法分割左心肌,利用面积长度法(即:每一层面心室面积与层厚(包含层间距)的乘积)计算心室体积(EDV,ESV)、每搏输出量(SV)、射血分数(EF)、左心室质量(ED_LVM,ES_LVM)等指标,将指标值与正常范围值进行比较,可以辅助临床医师评估心脏整体功能.18 名受试者各功能指标数据的均值和标准差如表3 所列.

表3 左心室功能指参数 (均值±标准差)Table 3.Left ventricular function parameters (Mean±standard deviation).

4.2 追踪结果评估

基于目标损失函数((3)式)对密集追踪网络进行训练,网络训练 300 个epoch 后得到的训练集与验证集的损失曲线如图8 所示,可以看出,所提出的网络可以较好地拟合,验证集的损失率可以降到 0.02 左右.

图8 训练集(蓝)与验证集(橙)的损失曲线Fig.8.Loss curves of training set (blue) and verification set(orange).

利用预测位移场(图9 中DVF)将ES 帧掩膜重采样并与ED 帧金标准掩膜进行比较,结果如图9 所示,可以看出,扭曲后的图像(图9 中Warped_label_ES)与金标准(图9 中label_ED)的LVC 和LVM 高度重合,表明了左心肌运动追踪的精度.

图9 利用预测位移场将ES 扭曲至ED 的示例图Fig.9.Example diagram of warping ES to ED using the predicted displacement field.

此外,预测位移场可以辅助临床医生诊断不同类型的CVDs,以达到早诊断早治疗的目的.图10给出了几个例子,不难发现,网络追踪得到的位移场确实足够好,足以表征病理类别中典型病例的心脏运动.其中正常人(NOR)的位移场在整个LVM上有大致相同的振幅,表明MOR 病例LVM 的收缩和增厚是同步的;对于MINF 病例,可以看到LVM 上位移流不均匀,部分心肌节段收缩和增厚明显小于其他节段,这是MINF 病例的典型症状;对于HCM 病例,LVM 上的位移流过大,这意味着LVM 收缩和增厚过度;相反,DCM 病例LVM 上的位移流很小,这是因为DCM 的心脏通常没有足够的收缩和增厚.

图10 不同病例类型预测的位移图Fig.10.Predicted displacement fields of different case types.

4.3 对比结果

本文通过两个方面展开对比实验,首先将本文方法与其他追踪算法所得的结果进行对比,其次将预测得到的位移场与模拟数据位移场进行对比.

1)为了验证本文追踪方法的准确性,将本文方法与二维B 样条自由变形(free-form deformation,FFD)配准方法[31]、Qin 等[19]提出的基于深度学习的L2 范数正则化配准网络(DL+L2)、未引入VAE 正则化器的追踪方法(DispFlow_UNet)得到的结果进行对比.在5 重交叉验证的结果中,LVC 和LVM 获得的Dice 系数、MCD 和HD 的平均值(标准差)如表4,加粗字体表示结果更好,结果显示DispFlow_UNet 优于DL+L2,FFD 法,而DispFlow_UNet_VAE 相比于DispFlow_UNet进一步取得更好的指标值.同时,各方法计算得到LVM 的Dice 和MCD 箱形图如图11 所示,更加直观地看出本文方法取得最高的Dice 值和最低的MCD 值,表明追踪精度最优.

图11 本文方法与其他方法的DM (a)和MCD (b)指标箱形图Fig.11.Box chart of DM (a) and MCD (b) indicators of the method presented in this paper and other methods.

表4 不同追踪方法Dice 系数、MCD 和HD 的对比Table 4.Comparison of Dice coefficients,MCD and HD of different tracking methods.

2)将预测得到的左心肌位移矢量u0=(x0,y0)和位移真值矢量 u1=(x1,y1) 进行对比,结果如图12 所示,可以直观看出,两者在方向和幅值上都很相近.为了更加清晰地比较两者的相似性,图13 将预测位移矢量场(红色)与位移真值矢量场(绿色)叠加显示,放大区域很好地展示了两者的相似性.除了视觉上直观的对比,本文还计算了定量评估运动追踪精度的评价指标AEE 和AEPE.其中,AEE=9.79±7.43,AEPE=2.28±1.19 (均值±偏差)单位分别为度和像素,说明预测位移矢量和位移值真值矢量重合程度较高.

图12 预测位移矢量与位移真值矢量对比Fig.12.Comparison between the predicted displacement field (left) and the true displacement field (right).

图13 预测位移矢量(红色)与位移真值矢量(绿色)对比Fig.13.Predicted displacement field (left) and the true displacement field (right).

4.4 模型泛化能力评估

在对模型进行泛化性研究中,利用ACDC 数据集训练的模型对15 例临床数据集直接进行预测,测试结果如表5 所列,其中比较了Dice 和MCD,结果表明本文所提出的方法在ACDC 数据集(内部预测集)和临床数据集(外部预测集)上都取得了较满意的结果,表明该方法在未知区域的数据上具有很好的通用性,具备较高的配准精度,这可能是由于生物力学规范的好处,它强制生成的变形场在生物力学上是合理的,并且对域转移问题不敏感.

表5 临床数据集上的模型泛化性能Table 5.Model generalization performance on clinical datasets.

5 结论

本文提出一种基于位移流U-Net 和生物力学引导的VAE 左心肌运动追踪方法DisFlow_UNet_VAE.首先针对CCMR 图像进行了数据预处理,包括数据归一化、数据增强和感兴趣区域的提取;其次,搭建压缩激励残差U-Net 网络精准分割ED 和ES 相位的左心肌,为后续的定量计算服务.接着运用DisFlow_UNet_VAE 追踪心肌运动得到密集位移场,利用位移场将ES 相位的心肌掩膜变形至ED 相位,并计算其与ED 相位掩膜的重叠程度以评估追踪精度;最终对预测集及临床数据集进行预测,结果表明,DisFlow_UNet_VAE 能够精确地追踪左心肌运动,具有较强的模型泛化能力.密集位移场能够辅助医师诊断不同类型的CVDs,并为后续应变分析提供了很好的支持,以实现对心肌组织的活性以及局部心脏功能的评估.该研究将进一步将 2D 网络推广到 三维(3D) 网络,对 3D 数据进行追踪,更好地满足临床需求.