CN109791617A - 低秩建模和并行成像的实时相位对比血流mri - Google Patents

Abstract

提出了一种新型的基于模型的成像方法，以实现稀疏采样的高分辨率实时相位对比磁共振成像(PC‑MRI)。所述实时相位对比血流磁共振成像方法包括：采集实时PC‑MRI数据，包括训练数据和成像数据；进行基于低秩模型的图像重建，包括基于采集的训练数据和成像数据的时间子空间和空间子空间的估计；基于重建的实时血流图像，计算速度图；以及基于计算出的速度图，进行血流量化分析。提出的方法在没有心电图(ECG)门控和呼吸控制的情况下能够实现2D高分辨率实时PC‑MRI，并首次实现3D成像。

Description

低秩建模和并行成像的实时相位对比血流MRI

技术领域

本发明主要涉及磁共振成像(MRI)，具体涉及稀疏采样的高分辨率实时相位对比磁共振成像(PC-MRI)。

背景技术

在过去的几十年中，相位对比磁共振成像(PC-MRI)已经发展成为一种对于心脏及大血管的血流动力学的定量及可视化的很有效的一种工具[1-5]。它推进了对于各种心血管疾病的理解和诊断，例如动脉粥样硬化[6]，动脉瘤[7]和动静脉畸形[8]。传统的PC-MRI[9,10]采用具有呼吸控制的心电图(ECG)同步电影采集来从多个心动周期中采集数据，进而重建出时空平均流速图。虽然这种方法在生物医学研究和临床实践中有很广泛的应用，但是它存在很多局限。例如它通常需要周期性或者准周期性心脏运动来确保有效的数据采集；不规则心脏运动常常导致采集时间的延长。另外，由于它的基本假设，这个方法只能得到多个心动周期的平均流动信息，不能够解决不规则的心脏运动，例如心律不齐所导致的不同心动周期血流变化的问题。捕获生理和/或病理的血流变化一直是PC-MRI研究的很重要的目标之一[11-14]。

无需ECG门控和呼吸控制的实时PC-MRI[15-17]是解决这些限制的很有前景的方向。然而，它需要更高的成像速度，对数据采集和图像重建都提出重大的挑战。已有很多技术用于推进在一个方向进行速度编码的空间二维实时PC-MRI。例如，先进的采集方法，像平面回波[18]、径向[19]、螺旋[20-23]采集方案都已经被用于实时PC-MRI。另外，一些从降采数据集中重建图像的加速方法已被提出。例如，伴随着并行成像的出现，敏感度编码(SENSE)[24]和广义自动校准部分并行采集(GRAPPA)[25]已经被应用于实时PC-MRI[26-30]。最近，使用正则化非线性反演[33]的基于模型的重建[31,32]被发展，通过联合重建质子密度图，相位图和一系列线圈敏感度，实现了空间分辨率1.5mm和时间分辨率25.6ms的2D实时血流成像。尽管已有很多方法用于实现2D实时PC-MRI，然而3D实时PC-MRI还没有任何报道。考虑到能够对复杂的血流流动模式和血管的几何结构进行成像，实现3D实时PC-MRI是很有意义的。然而，3D实时血流成像通常涉及到更具挑战性的在空间分辨率、时间分辨率、成像时间和更大计算量之间的权衡问题。

在本实施例中，我们对于稀疏采集的实时PC-MRI提出一种新的基于模型的方法，它基于对新型低秩模型与并行成像的集成。所述方法实现了无需ECG门控和呼吸控制的实时PC-MRI，并且很好地解决了常规的电影方法无法获得的每次心跳之间的血流变化问题。我们还将这种成像方法进行拓展，通过进一步稀疏建模的整合，实现在三个方向进行血流编码的3D实时PC-MRI。更重要的是，我们第一次验证了3D实时PC-MRI的可行性。所述方法的有效性已由体内实验在2D和3D实时PC-MRI中进行了系统地评估。

发明内容

本发明的一个方面提供了一种实时相位对比血流磁共振成像(MRI)方法，包括：采集实时相位对比MRI(PC-MRI)数据，包括训练数据和成像数据；基于采集的训练数据和成像数据进行基于低秩模型的图像重建；基于重建的实时血流图像计算速度图；根据计算出的速度图进行血流量化分析。

在根据本公开的方法中，进行基于低秩模型的图像重建还包括：针对训练数据进行时间插值；针对插值后的训练数据进行时间子空间估计；针对采集的成像数据进行线圈灵敏度的估计；基于估计出的时间子空间、线圈灵敏度和采集的成像数据进行空间子空间估计。

根据本公开的方法可以被用于2D和3D实时PC-MRI。

附图说明

附图中相同的附图标记指代贯穿各个部分的相同或功能上类似的元素，附图标记和下面的详细描述将一起被并入说明书，形成说明书的一部分。这些图示用于进一步说明本发明的具体实现结果，同时帮助更好地解释本发明的原理和优点。

图1图示根据本公开的实施例的提出的(k,t)-空间的采样机制；

图2图示根据本公开的实施例的提出的实时PC-MRI方法的流程图；

图3示出根据本公开的实施例的对两个健康受试者的2D实时血流成像和2D电影血流成像的对比；

图4示出根据本公开的实施例的一名健康受试者来自2D实时血流成像的重建速度波形；

图5示出2D实时PC-MRI和2D电影PC-MRI之间的一致性评价(Bland-Altman)分析；

图6示出根据本公开的实施例的一名23岁心律不齐患者的2D实时PC-MRI；

图7示出根据本公开的实施例的一名72岁心律不齐患者的2D实时PC-MRI；

图8示出根据本公开的实施例的一名健康受试者分别由传统3D电影血流成像和提出的3D实时血流成像得到的速度图；

图9示出根据本公开的实施例的一名健康受试者的3D实时血流成像的流量波形；

图10示出3D实时PC-MRI和3D电影PC-MRI之间峰值速度和每搏输出量之间的一致性评价(Bland-Altman)分析；

图11示出根据本公开的实施例的一名心率不齐患者由3D实时PC-MRI得到的重建流量波形图和迹线可视化图。

技术人员将会认识到图中的元素出于简化和清晰的目的而被图示，且没有必要按比例绘制。比如，图中一些元素相比于其他元素，维度可能被放大，这是为了能帮助更好地理解本发明的实施例。

已经在附图中通过常规符号在适当处呈现了方法和设备组件，仅示出了与理解本发明的实施例相关的具体细节，以便不会因为可受益于此描述的本领域普通技术人员所显而易见的细节而妨碍对本公开的理解。

具体实施方式

缩略词

PC-MRI：相位对比磁共振成像(phase-contrast magnetic resonance imaging)；

ECG：心电图(electro-cardiogram)；

SENSE：灵敏度编码(sensitivity encoding)；

GRAPPA：广义自动校准部分并行采集(generalized autocalibrating partiallyparallel acquisitions)；

VENC：编码速度(encoding velocity)；

AAo：升主动脉(ascending aorta)；

Dao：降主动脉(descending aorta)。

在本发明具体实施方式的最后部分，列出了所有的参考文献。请注意所有列出和提到的参考文献全部被包含在本发明中用作参考。

实施例：实时PC-MRI

理论

在忽略读出梯度时间范围内血液流动的情况下，我们可以对相位对比磁共振实时成像过程进行如下建模：

d_v,i(k,t)＝∫S_i(r)ρ_v(r,t)e^-j2πk·rdr+η_v,i(k,t) (1)

其中ρ_v(r,t)代表动态图像，当v＝1时与血流补偿图像序列相关联，当v＝2,…,N_v时与血流编码图像序列相关联，S_i(r)表示第i个接收线圈的敏感度分布，i＝1,2,…,N_c，d_v,i(k,t)和η_v,i(k,t)分别代表从(k,t)-空间测量的数据以及测量噪声。另外，v代表不同的血流编码(v＝1：血流补偿数据；v＝2,…,N_v：每个血流编码方向的血流编码数据)，k是k-空间的采样位置，r是相应的空间位置，t是时间。这里的目标是从欠采样数据{d_v,i(k,t)}中重建ρ_v(r,t)，然后计算速度图其中Δφ(r,t)＝∠ρ_v(r,t)-∠ρ₁(r,t)表示血流编码图像序列和血流补偿图像序列的相位差，VENC表示在扫描前预先设定的编码速度。而在实时PC-MRI中，由于没有ECG门控的数据共享，所采集得到的(k,t)-空间数据通常处于非常高欠采样状态。对{d_v,i(k,t)}进行直接求逆的运算会导致非常严重的混叠伪影，进而得到不准确的速度测量值。

这里我们引入基于低秩模型的重建方法及并行成像来解决这个问题。为了方便表述，我们采用离散的图像模型，每一个血流图像序列可以被表示为在空间和时间上分布的卡索拉蒂(Casorati)矩阵[34]，即：

与心脏动态成像中的应用[35-37]相似，由于随时间连续变化的动态图像序列都具有很强的时间空间相关性，于是这里每一个C_v都符合一个低秩近似。进一步，由于血流编码本身的特质，在不同的血流图像序列之间同样具有很强的时间空间相关性。为了利用这样一个性质，我们引入了下面的联合卡索拉蒂(Casorati)矩阵：

使其符合低秩结构，即rank(C)≤L，其中L是指矩阵C的秩，并且满足L≤min(M,N)。存在很多方法可以施加低秩约束[34,36,38,39]。这里，我们通过矩阵因子分解来利用一个显式的秩约束，即C＝UV，其中且在这个低秩表达中，U的列和V的行分别是支撑C的空间子空间和时间子空间。

接下来，我们进一步用公式表示低秩约束重建问题。首先，注意到，通过矩阵向量化标注，我们可以将公式(1)改写成：

d_i＝Ω(F_sS_iC)+n_i (4)

其中，d_i表示测量得到的数据，Ω表示稀疏采样运算符，F_s表示空间傅里叶变换矩阵，S_i和n_i表示敏感度分布以及相应的测量噪声。通过施加低秩约束，图像重建问题可用以下方程进行表示：

该问题是一个非凸优化问题，目前已有一些现有的算法[40,41]可以用来解决该问题。

图像重建问题可以被进一步简化。将早期的心脏成像工作[34,36,37]进行拓展，我们可以采用特殊的数据采集机制，利用获取的训练数据预估计出时间子空间V。我们采用依次交错采集k-空间中心相位编码数据行与k-空间外周相位编码数据行分别获取训练数据和成像数据的方式，如图1(a)所示。其中，训练数据由k-空间中心采集得到，而成像数据由在剩余的(k,t)-空间进行随机采集得到。采用如此数据采集策略，两组数据可以为低秩模型提供两部分相互之间互为补充的信息：训练数据具有非常高的时间分辨率，成像数据具有非常高的空间分辨率。从训练数据中，利用主成分分析法估计得到时间子空间V[34,42]。然后利用分析得到的信号动态变化特征以及从剩余k-空间采集到的成像数据估计得到空间子空间U。此外，我们还进行了时间维度的数据插值以将训练数据与成像数据这两组数据的采集时间点进行匹配，如图1(b)中所示，这里我们主要是将训练数据线性插值到与成像数据采集相同的时间点上。利用这样的数据采集策略，我们最终重建得到的实时血流图像可以获得的时间分辨率为2×N_v×TR。另外需要注意的是，这里的敏感度分布S_i是由流动参考序列的(k,t)-空间在时间维度上进行平均之后的数据进行估计得到的。

在前面估计出之后，我们可以通过解如下的凸优化问题来确定U：

通过预先将时域子空间估计出来之后，原先非凸的低秩矩阵恢复问题简化为了一个简单的最小二乘问题。而在求解得到之后，由前文定义的联合卡索拉蒂(Casorati)矩阵即可通过方程进行重建得到，继而可以得到每一个血流图像序列并估计出相应的流速数值。其中U/V/C是指真实值，是指相应的重建得到的值。

公式(6)可充分用于2D实时血流成像。然而，对于3D实时血流成像，则会涉及到由于高度欠采样数据导致的超大规模的运算量问题。这里，考虑到联合卡索拉蒂(Casorati)矩阵C在空间域的稀疏表示，进一步利用稀疏约束求解空间子空间U[36,43]。这样，空间子空间的重建问题可以用以下方程表示：

其中F_t表示时间傅里叶变换矩阵，λ是正则化参数，vec(·)运算符是将矩阵按列连成向量。

总结上述提出方法的流程图展示在图2中。流程图主要包括以下三部分：数据采集、图像重建和后处理。

实现

我们通过在体实验对所提出的2D和3D实时PC-MRI方法进行系统的评估。实验在Philips Achieva 3.0T TX全身磁共振双源扫描仪(Achieva,Philips Medical System,Best,The Netherlands)上进行的，并使用了原装的32通道心血管线圈。我们在原装的基于梯度回波(GRE)的脉冲序列中进行编程，实现了如图1(a)中所示的数据采集策略的在线实时采集方式。在此，未使用任何ECG门控或呼吸控制来辅助数据采集重建。而用于与之进行对照的传统电影PC-MRI则使用了厂家原装的基于GRE的脉冲序列，并使用了回溯式的ECG门控。

在2D实时PC-MRI实验中，我们征集了10位无任何心血管疾病征兆的健康志愿者(7名男性，年龄：22-29岁，中位值：25岁)进行磁共振扫描。此外，我们还对2位心率不齐的病人(2名男性，分别为23岁和72岁)进行了磁共振血流测量，进而对所提出的血流实时成像的方法在心率不齐时血流测量的可行性进行验证。该实验已获清华大学伦理审查委员会批准，并且所有受试人员均签写了知情同意书。在本实验中，相位对比磁共振电影成像和实时成像均在与升主动脉和降主动脉的管腔垂直的一个平面内进行扫描，在扫描过程中，受试者均处于自由呼吸的状态，这里均采用沿头脚的一个方向进行速度编码。在相位对比磁共振电影成像序列扫描中，我们采用了回溯式心电同步信号触发的方式进行数据采集，并根据每位受试者各自的心跳周期进行扫描前的心率设置。同时，我们连续采集三组相位对比磁共振电影成像的数据然后进行平均以减少由于呼吸运动产生的伪影。对于相位对比磁共振电影成像和实时成像，我们应用了以下的成像参数：视野＝240mm×225mm，矩阵大小＝132×124，空间分辨率＝1.80mm×1.80mm，层厚＝5mm，TR/TE＝4.5/2.8ms，翻转角＝10°，VENC＝200cm/s。其中，对于所提出的相位对比磁共振实时成像来说，时间分辨率为4×TR＝18ms，而对于电影成像来说，在采集28个心跳时相的条件下时间分辨率为36ms。两组成像序列的采集时间均在94s附近。

此外，我们还进行了3D实时血流成像实验。我们征集了5位健康志愿者和一名心律不齐患者。作为对比，我们还利用原装的2x SENSE[24]的成像序列进行3D电影成像数据采集。实验中，电影成像和实时成像的目标扫描部位均为包含整个主动脉，且在整个扫描过程中，受试者均处于自由呼吸的状态。具体的成像参数：视野＝180mm×256mm×43mm(FH/RL/AP)，空间分辨率＝2.40mm×2.40mm×2.40mm，矩阵大小＝76×108×18，TR/TE＝4.4/2.6ms，翻转角＝5°，VENC＝200/150/150cm/s(FH/RL/AP)，时间分辨率为8×TR＝35.2ms。

对于相位对比磁共振血流电影成像，血流补偿和血流编码图像序列分别直接由满采样的(k,t)-空间数据进行二维傅里叶变换进而重建得到。对于所提出的相位对比磁共振血流实时成像，我们按照图2中的流程图进行相应的数据处理。具体而言，我们首先利用采集得到的训练数据进行时间维度的线性插值，使其与近邻的相位编码行具有相同的采样时间点，进而由插值后的训练数据计算得到时间域子空间V。然后我们对所有采集得到的(k,t)-空间数据进行时间维方向的平均，根据平均值后的k-空间数据进行傅里叶变换进而估算出线圈敏感度分布图S_i。在将血流补偿和血流编码序列对应采集到的数据进行排列形成联合卡索拉蒂(Casorati)矩阵之后，基于前面计算得到的S_i和V，我们可以分别根据公式(6)和公式(7)计算分别得到2D和3D实时成像空间域的子空间U。为了进一步提高运算效率，可以首先对采集到的训练数据和成像数据进行通道压缩(例如[44])。在做完图像重建之后，我们另外进行了相位矫正[45]来修正由于涡流导致的相位偏差。最后可以将得到的流速分布图用于定量的流量分析。

对于健康受试者的体内实验，我们评估了从电影方法测得的流速与所提出方法测得的流速之间的一致程度。具体而言，我们对从两种方法得到的峰值速度和每搏输出量进行了一致性评价(Bland-Altman)分析以及配对学生t检验(Student's t-test)。在这里，峰值速度被定义为一个心动周期内的最大速度，每搏输出量是升主动脉内一个心动周期内流速的积分。对于心律失常患者的实验，我们以外部ECG记录的心脏运动为参考来评估所提出的方法捕获流速变化的有效性。

2D实时PC-MRI的实验结果

图3示出了根据本公开的实施例对两名健康志愿者进行2D实时血流成像和2D电影血流成像的对比。此处我们主要比较了由两种成像方法重建得到的在收缩期和舒张期时的模值图像和速度图。从该图中可以看出，由所提出的相位对比磁共振实时成像方法得到的模值图像和速度图的重建质量与电影成像方法相当。虽然这两种方法都可以较好地重建出血管结构，但实时成像方法比传统的电影成像方法具有更好地减少运动伪影的性能。

此外，我们对升主动脉(AAo)和降主动脉(DAo)中分别提取了管腔以内的区域作为感兴趣区域，并分析该区域内所有像素的平均血流速度。图4示出了根据本公开的实施例对一名健康受试者的重建速度波形。具体来说，图4(a)和(b)展示了10个心动周期中升主动脉(AAo)和降主动脉(DAo)的血液流速波形。从这两个图中可以清楚地看到，所提出的方法可以很好地捕捉不同心跳周期的周期性血流变化。我们进一步地分析比较了由所提出的血流实时成像方法得到的血流速度波形图与由传统的血流电影成像得到的速度波形图之间的关联性。我们对于由所提出的血流实时成像方法得到的连续30个心跳周期的血流速度波形图进行平均，得到一个综合的心跳周期的速度波形图，然后将其与由传统的血流电影成像得到的血流速度波形图进行比较。升主动脉和降主动脉的比较结果分别见于图4(c)和(d)。从结果中可明显看出由血流实时成像方法得到的平均之后的血流速度波形图与由电影成像得到的血流速度波形图在整体随时间变化的趋势上几乎完全重合，尤其是根据这两种方法对升主动脉和降主动脉两个感兴趣区域内峰值流速的估计也是非常相近。

我们进一步对所有10名健康志愿者用相位对比磁共振电影成像和实时成像这两种方法得到的结果进行了统计分析。图5(a)和(b)分别展示了对用两种方法得到的峰值流速和每搏输出量进行一致性评价(Bland-Altman)分析的结果。实时成像方法中的峰值流速和每搏输出量是30个连续心动周期的平均值。从结果中可以看出由所提出的血流实时成像方法得到的血流测量值与传统的血流电影成像方法得到的测量值表现出了非常好的一致性。此外，我们还对两种方法进行了配对学生t检验(Student's t-test)的统计分析。这两种方法在峰值流速的测量和在每搏输出量方面的测量的相关系数分别为0.94(P<0.001)和0.90(P＝0.0002)。这一统计结果进一步证实了这两种方法存在较强的相关性。

图6示出了根据本公开的实施例对一名23岁心律不齐的受试者(具有轻度心律失常)进行2D实时PC-MRI的结果。图6(a)展示了心电门控采集的升主动脉与降主动脉的心电信号随时间变化的波形图，图6(b)展示了在心律不齐期间截取的三个峰值时间点对应的由实时成像方法重建得到的模值图像和速度图。与设想的结果相同，所提出的血流实时成像方法能够很好地重建出不同心跳周期的血流速度的相对变化。具体来说，由图6(b)所示，所提出的血流实时成像方法可以很好地捕捉到在心律不齐阶段出现的一个骤降的血流速度变化。可以注意到，这种类型的血流动态变化的情况并不能够被传统的血流电影成像方法测量得到。再者，由该图还可以看出，根据所提出的血流实时成像方法得到的血流速度随时间变化的波形图与由ECG心电记录仪在数据采集期间同步记录得到的心电信号随时间变化的波形图有非常好的相关性。这里我们还给出了由所提出的方法重建得到的结果中选取三个时间点所进行的模值图像和速度图，此外，我们展示了来自所提出的方法的三个快照图像。显然，速度图证实了心律失常期内剧烈的血流变化。

图7示出了对一名72岁心律不齐的受试者进行2D实时PC-MRI的结果。图7(a)展示了升主动脉与降主动脉的速度波形。所提出的方法很好地捕捉不规则的血流变化，这比上一名患者更加显著。此外我们在图7(b)中展示了重建得到的模值图像和速度图。

3D实时PC-MRI的实验结果

图8示出了对一名健康志愿者分别利用3D血流电影成像和提出的实时成像这两种成像方法得到的重建速度图的结果。图8(a)和(b)分别展示的是收缩期和舒张期。可以看出，提出的实时成像方法提供了与电影成像方法可比较的成像质量。

图9示出了对一名健康志愿者利用提出的3D实时血流成像方法得到的重建血流波形。在图9(a)中标记出的升主动脉和降主动脉的血流波形分别展示在图9(b)和(c)。可以看出，所提出的方法能够显示血流量变化。我们进一步将来自3D实时成像方法的流动波形在10个心动周期内平均为一个合成流动波形，并将其与来自电影成像方法获取的波形进行比较。如图9(d)和(e)所示，所提出方法的平均流量波形与从电影成像方法获得的平均流量波形很好地相关，这是非常理想的。

此外，图10中示出了对健康受试者的峰值速度和每搏输出量之间的一致性评价(Bland-Altman)分析。具体地，图10(a)中示出的峰值速度和图10(b)中示出的每搏输出量的一致性评价(Bland-Altman)分析被用来对5名健康志愿者(三名男性两名女性，平均年龄21岁)用传统的3D相位对比磁共振电影成像和新提出的3D实时成像这两种方法得到的结果进行了比较。实时成像方法中的峰值流速和每搏输出量是10个连续心动周期的平均值。这些结果进一步证实了这两种方法之间的一致性。

图11示出了对一名心律不齐患者的3D实时PC-MRI的结果。图11(a)和(b)分别示出了升主动脉和降主动脉的重建血流量波形。图11(c)展示了一个正常心跳周期内的时间帧(A和B)和一个异常心跳周期内的时间点(C和D)对应的3D迹线可视化结果图。由该图可见，所提出的三维磁共振血流实时成像方法可以很好地解决每次心跳之间拍摄的病理变化问题，特别是它可以明显地捕捉到在异常心跳周期中出现的血流状态变化。而这是传统电影成像方法无法检测到的。

讨论

在本实施例中，我们介绍了一种新的实时血流成像方法并系统地证明了其在体内实验中的有效性。在这里，值得重申所提出的方法的关键特征。首先，它可以用作传统电影成像方法的可行替代方案，因为它为健康受试者提供了可比较的图像质量和流量信息。其次，它能够解析每次心跳之间的生理和/或病理血流变化，这是传统电影方法无法获得的。而这些信息通常在临床上很重要(例如，用于评估心律失常)。

与其他基于模型的图像重建算法一样，这里提出的方法涉及模型选择(即，对秩L的选择)。通常，对于L的选择需要在模型表征能力、数据样本的大小(即采集时间)和信噪比(SNR)等多方面进行权衡[36]。在本研究中，我们通过手动选择L来权衡上述因素，并发现其具有较为稳定的良好的重建结果。在今后的研究中，其他原理模型选择方法(如[46,47])值得进一步的分析和探讨。

我们提出的2D实时血流成像计算效率比较高。在64GB RAM、3.47GHz CPU的工作站上，对在体数据(94s的实时采集)进行图像重建运算大约只需10分钟。但是，对于3D实时血流成像，对在体数据(20分钟的实时采集)进行图像重建运算将会超过一小时。为了增强实用性，通过图形处理单元的实现可以提高计算效率。在未来的研究中这一研究方向值得继续深入探索。

本项工作主要集中在新型实时血流成像技术的开发上，该技术应作为我们后续临床研究的基础。鉴于所提出的方法很好地解决了每次心跳之间的流量变化问题，它可以为具有明显心律不齐的患者提供更多关于血流动力学的信息。在未来的工作中，我们计划针对各种潜在临床应用(例如心房颤动，房性早收缩或先天性心脏病)进行系统研究。

结论

本实施例提出了一种不用ECG门控和呼吸控制的高分辨率实时PC-MRI的基于模型的新方法，并且首次实现3D实时PC-MRI。它具有新颖的低秩模型结合并行成像，它们共同实现了从高度欠采样(k,t)-空间实时PC-MRI空间数据的高质量重建。该方法的有效性和效用性已经在针对2D和3D实时PC-MRI的体内实验中被证明。我们期待该提议方法将增强实时PC-MRI在各种临床应用中的实用性。

在前面的说明书中，已经描述了具体的实施例。然而，本领域的普通技术人员应该理解，在不脱离所附的权利要求书中提出的本发明的范围的情况下，可以进行各种修改和变化。因此，说明书和附图是说明性的而非限制性的，并且所有符合上述条件的修改旨在被包括在本教导的范围内。

益处、优点、问题的解决方案以及任何可能导致益处、优点或解决方案发生或变得更加明显的任何要素均不被解释为任何或所有权利要求的关键、必需或基本的特征或要素。本发明仅由所附权利要求限定，包括在本申请待审查期间所作的任何修改以及授权发布的那些权利要求的所有等同物。

此外，在该文件中，诸如第一和第二、顶部和底部等的关系术语仅用于将一个实体或动作与另一个实体或动作区分开，而不一定要求或暗示它们之间的任何实际的关系或顺序。术语“包括”、“包含”、“具有”、“含有”等关联词旨在非排他性地包含类似元素，由此，过程、方法、物品或装置等不仅包括原有的元素，还包括未明确列出的其他相似元素。如果没有更多的约束条件，那么之前冠以“包含……”、“具有……”、“包括……”的元素不排除在包括、具有、包含该元素的过程、方法、物品或装置中额外等同元素的存在。术语“a”和“an”被定义为一个或多个，除非在此另有明确说明。术语“大体上”、“基本上”、“大约”、“约”或任何其他形式，被定义为本领域普通技术人员所理解的“接近于”，且在一个非限制实施例中，该术语被定义为在10％内，在另一实施例中被定义为在5％内，在又一实施例中被定义为在1％内，还有一个实施例被定义为在0.5％内。本文所使用的术语“耦合”被定义为连接，尽管不一定直接且不一定机械地连接。以一种方式“配置”的设备或结构至少以该方式配置，但也可以以未列出的方式配置。

应该理解，一些实施例可以包括一个或多个通用或专用处理器(或“处理设备”)，例如微处理器、数字信号处理器、定制处理器和现场可编程门阵列(FPGA)以及单独的存储程序指令(包括软件和硬件)，来控制一个或多个处理器以结合某些非处理器电路实现本文描述的方法和/或设备的一些、大部分或全部功能。或者，可以通过没有存储程序指令的状态机或者在一个或多个专用集成电路(ASIC)中实现一些或全部功能，其中每个功能或某些功能的一些组合被实现为定制逻辑。当然，也可以使用这两种方法的组合。

此外，可以将实施例实现为其上存储有计算机可读代码的计算机可读存储介质，用于对计算机(例如，包括处理器)进行编程以执行如本文所述和要求保护的方法。这样的计算机可读存储介质的示例包括但不限于硬盘、CD-ROM、光存储设备、磁存储设备、ROM(只读存储器)、PROM(可编程只读存储器)、EPROM(可擦除可编程只读存储器)、EEPROM(电可擦除可编程只读存储器)和闪存。此外，在本文公开的概念和原理的指导下，预期普通技术人员尽管可能花费由例如可用时间、当前技术和经济考虑所激励的大量努力和许多设计选择，将易于能够以最少的试验来产生这些软件指令和程序及IC。

提供摘要以使得读者快速了解公开技术的本质。提交时的想法是，它不会被用来解释或限制权利要求的范围或含义。另外，在前面的具体实施方式中，可以看出，出于简化本公开技术的目的，各种特征在各种实施例中被组合在一起。本公开的方法不应被解释为反映所要求保护的实施例需要比每个权利要求中明确记载的特征更多的意图，而是如所附权利要求所反映的，发明主题在于少于单个公开实施例的所有特征。因此，所附权利要求在此被并入到具体实施方式中，其中每个权利要求本身作为单独要求保护的主题。

资金支持

本项工作部分地得到了“十三五规划”期间国家重点研发计划(2016YFC1301601)和美国国立卫生研究院(NIH-RO1-EB013695)的支持。

参考文献

[1]Van Dijk,P.:Direct cardiac NMR imaging of heart wall and bloodflow velocity.J Comput Assist Tomogr 8(3),429–436(1984)

[2]Nayler,G.,Firmin,D.,Longmore,D.:Blood flow imaging by cinemagnetic resonance.J Comput Assist Tomogr 10(5),715–722(1986)

[3]Gatehouse,P.D.,Keegan,J.,Crowe,L.A.,Masood,S.,Mohiaddin,R.H.,Kreitner,K.-F.,Firmin,D.N.:Applications of phase-contrast flow and velocityimaging in cardiovascular MRI.Eur Radiol 15(10),2172–2184(2005)

[4]Markl,M.,Kilner,P.J.,Ebbers,T.:Comprehensive 4D velocity mappingof the heart and great vessels by cardiovascular magnetic resonance.JCardiovasc Magn Reson 13(7),1–22(2011)

[5]Nayak,K.S.,Nielsen,J.-F.,Bernstein,M.A.,Markl,M.,D.Gatehouse,P.,M.Botnar,R.,Saloner,D.,Lorenz,C.,Wen,H.,S.Hu,B.,Epstein,F.H.,N.Oshinski,J.,Raman,S.V.:Cardiovascular magnetic resonance phase contrast imaging.JCardiovasc Magn Reson 17(1),1(2015)

[6]Markl,M.,Frydrychowicz,A.,Kozerke,S.,Hope,M.,Wieben,O.:4D flowMRI.J Magn Reson Imaging 36(5),1015–1036(2012)

[7]Hope,T.A.,Hope,M.D.,Purcell,D.D.,von Morze,C.,Vigneron,D.B.,Alley,M.T.,Dillon,W.P.:Evaluation of intracranial stenoses and aneurysms withaccelerated 4D flow.Magn Reson Imaging 28(1),41–46(2010)

[8]Ansari,S.,Schnell,S.,Carroll,T.,Vakil,P.,Hurley,M.,Wu,C.,Carr,J.,Bendok,B.,Batjer,H.,Markl,M.:Intracranial 4D flow MRI:Toward individualizedassessment of arteriovenous malformation hemodynamics and treatment-inducedchanges.Am J Neuroradiol34(10),1922–1928(2013)

[9]Lenz,G.W.,Haacke,E.M.,White,R.D.:Retrospective cardiac gating:Areview of technical aspects and future directions.Magn Reson Imaging 7(5),445–455(1989)

[10]Dyverfeldt,P.,Bissell,M.,Barker,A.J.,Bolger,A.F.,Carlh¨all,C.-J.,Ebbers,T.,Francios,C.J.,Frydrychowicz,A.,Geiger,J.,Giese,D.,Hope,M.D.,Kilner,P.J.,Kozerke,S.,Myerson,S.,Neubauer,S.,Wieben,O.,Markl,M.:4D flowcardiovascular magnetic resonance consensus statement.J Cardiovasc Magn Reson17(1),1–19(2015)

[11]Finn,J.P.,Nael,K.,Deshpande,V.,Ratib,O.,Laub,G.:Cardiac MRimaging:State of the technology.Radiology 241(2),338–354(2006)

[12]Thavendiranathan,P.,Verhaert,D.,Walls,M.C.,Bender,J.A.,Rajagopalan,S.,Chung,Y.-C.,Simonetti,O.P.,Raman,S.V.:Simultaneous right andleft heart real-time,free-breathing CMR flow quantification identifiesconstrictive physiology.J Am Coll Cardiol:Cardiovasc Imaging 5(1),15–24(2012)

[13]Kowallick,J.,Joseph,A.,Unterberg-Buchwald,C.,Fasshauer,M.,VanWijk,K.,Merboldt,K.,Voit,D.,Frahm,J.,Lotz,J.,Sohns,J.:Real-time phase-contrast flow MRI of the ascending aorta and superior vena cava as a functionof intrathoracic pressure(Valsalva manoeuvre).Br J Radiol 87(1042),20140401(2014)

[14]H.,Barth,P.,Gieseke,J.,Müller,K.,Burchert,W.,Esdorn,H.,Kececioglu,D.,Beerbaum,P.,Laser,K.T.:Impact of respiration on stroke volumesin paediatric controls and in patients after Fontan procedure assessed by MRreal-time phase-velocity mapping.Eur Heart J Cardiovasc Imaging 16(2),198–209(2015)

[15]Riederer,S.J.,Wright,R.C.,Ehman,R.L.,Rossman,P.J.,Holsinger-Bampton,A.E.,Hangiandreou,N.J.,Grimm,R.C.:Real-time interactive color flow MRimaging.Radiology181(1),33–39(1991)

[16]Liu,C.-Y.,Varadarajan,P.,Pohost,G.M.,Nayak,K.S.:Real-time color-flow MRI at 3 T using variable-density spiral phase contrast.Magn ResonImaging 26(5),661–666(2008)

[17]Joseph,A.A.,Merboldt,K.-D.,Voit,D.,Zhang,S.,Uecker,M.,Lotz,J.,Frahm,J.:Real-time phase-contrast MRI of cardiovascular blood flow usingundersampled radial fast low-angle shot and nonlinear inversereconstruction.NMR Biomed 25(7),917–924(2012)

[18]Eichenberger,A.,Schwitter,J.,McKinnon,G.,Debatin,J.,VonSchulthess,G.:Phase-contrast echo-planar MR imaging:Real-time quantificationof flow and velocity patterns in the thoracic vessels induced by Valsalva’smaneuver.J Magn Reson Imaging 5(6),648–655(1995)

[19]Shankaranarayanan,A.,Simonetti,O.P.,Laub,G.,Lewin,J.S.,Duerk,J.L.:Segmented k-space and real-time cardiac cine MR imaging with radialtajectories.Radiology 221(3),827–836(2001)

[20]Gatehouse,P.D.,Firmin,D.N.,Collins,S.,Longmore,D.B.:Real timeblood flow imaging by spiral scan phase velocity mapping.Magn Reson Med 31(5),504–512(1994)

[21]Nayak,K.S.,Pauly,J.M.,Kerr,A.B.,Hu,B.S.,Nishimura,D.G.:Real-timecolor flow MRI.Magn Reson Med 43(2),251–258(2000)

[22]Steeden,J.A.,Atkinson,D.,Taylor,A.M.,Muthurangu,V.:Assessingvascular response to exercise using a combination of real-time spiral phasecontrast MR and noninvasive blood pressure measurements.J Magn Reson Imaging31(4),997–1003(2010)

[23]Kowalik,G.T.,Steeden,J.A.,Pandya,B.,Odille,F.,Atkinson,D.,Taylor,A.,Muthurangu,V.:Real-time flow with fast GPU reconstruction for continuousassessment of cardiac output.J Magn Reson Imaging 36(6),1477–1482(2012)

[24]Pruessmann,K.P.,Weiger,M.,Scheidegger,M.B.,Boesiger,P.:SENSE:Sensitivity encoding for fast MRI.Magn Reson Med 42(5),952–962(1999)

[25]Griswold,M.A.,Jakob,P.M.,Heidemann,R.M.,Nittka,M.,Jellus,V.,Wang,J.,Kiefer,B.,Haase,A.:Generalized autocalibrating partially parallelacquisitions(GRAPPA).Magn Reson Med 47(6),1202–1210(2002)

[26]Weiger,M.,Pruessmann,K.P.,Boesiger,P.:Cardiac real-time imagingusing SENSE.Magn Reson Med 43(2),177–184(2000)

[27]Hoogeveen,R.,Leone,B.,van der Brink,J.:Real-time quantitativeflow using EPI and SENSE.In:Proceedings of the 9th Annual Meeting of theInternational Society for Magnetic Resonance in Medicine:21-27 April 2001；Glasgow,Scotland,vol.9,p.114(2001)

[28]H.,Gieseke,J.,Barth,P.,Hoogeveen,R.,Esdorn,H.,A.,Meyer,H.,Beerbaum,P.:Flow volume and shunt quantification inpediatric congenital heart disease by real-time magnetic resonance velocitymapping.Circulation 109(16),1987–1993(2004)

[29]Wintersperger,B.J.,Nikolaou,K.,Dietrich,O.,Rieber,J.,Nittka,M.,Reiser,M.F.,Schoenberg,S.O.:Single breath-hold real-time cine MR imaging:Improved temporal resolution using generalized autocalibrating partiallyparallel acquisition(GRAPPA)algorithm.Eur Radiol 13(8),1931–1936(2003)

[30]Jones,A.,Steeden,J.A.,Pruessner,J.C.,Deanfield,J.E.,Taylor,A.M.,Muthurangu,V.:Detailed assessment of the hemodynamic response to psychosocialstress using real-time MRI.J Magn Reson Imaging 33(2),448–454(2011)

[31]Joseph,A.,Kowallick,J.T.,Merboldt,K.-D.,Voit,D.,Schaetz,S.,Zhang,S.,Sohns,J.M.,Lotz,J.,Frahm,J.:Real-time flow MRI of the aorta at aresolution of 40 msec.J Magn Reson Imaging 40(1),206–213(2014)

[32]Tan,Z.,Roeloffs,V.,Voit,D.,Joseph,A.A.,Untenberger,M.,Merboldt,K.D.,Frahm,J.:Model-based reconstruction for real-time phase-contrast flowMRI:Improved spatiotemporal accuracy.Magn Reson Med(2016).doi:10.1002/mrm.26192

[33]Uecker,M.,Hohage,T.,Block,K.T.,Frahm,J.:Image reconstruction byregularized nonlinear inversion-joint estimation of coil sensitivities andimage content.Magn Reson Med60(3),674–682(2008)

[34]Liang,Z.-P.:Spatiotemporal imaging with partially separablefunctions.In:Proceedings of IEEE International Symposium on BiomedicalImaging:April 2007；Washington,DC,USA,pp.988–991(2007)

[35]Zhao,B.,Haldar,J.P.,Brinegar,C.,Liang,Z.-P.:Low rank matrixrecovery for real-time cardiac MRI.In:Proceedings of IEEE InternationalSymposium on Biomedical Imaging:April2010；Rotterdam,The Netherlands,pp.996–999(2010)

[36]Zhao,B.,Haldar,J.P.,Christodoulou,A.G.,Liang,Z.-P.:Imagereconstruction from highly undersampled-space data with joint partialseparability and sparsity constraints.IEEE Trans Med Imaging 31(9),1809–1820(2012)

[37]Christodoulou,A.G.,Zhang,H.,Zhao,B.,Hitchens,T.K.,Ho,C.,Liang,Z.-P.:High-resolution cardiovascular MRI by integrating parallel imaging withlow-rank and sparse modeling.IEEE Trans Biomed Eng 60(11),3083–3092(2013)

[38]Haldar,J.P.,Liang,Z.-P.:Spatiotemporal imaging with partiallyseparable functions:A matrix recovery approach.In:Proceedings of IEEEInternational Symposium on Biomedical Imaging:April 2010；Rotterdam,TheNetherlands,pp.716–719(2010)

[39]Lingala,S.G.,Hu,Y.,DiBella,E.,Jacob,M.:Accelerated dynamic MRIexploiting sparsity and low-rank structure:k-t SLR.IEEE Trans Med Imaging 30(5),1042–1054(2011)

[40]Haldar,J.P.,Hernando,D.:Rank-constrained solutions to linearmatrix equations using power factorization.IEEE Signal Proc Let 16(7),584–587(2009)

[41]Shen,Y.,Wen,Z.,Zhang,Y.:Augmented lagrangian alternatingdirection method for matrix separation based on low-rank factorization.OptimMethod Softw 29(2),239–263(2014)

[42]Gupta,A.S.,Liang,Z.-P.:Dynamic imaging by temporal modeling withprincipal component analysis.In:Proceedings of the 9th Annual Meeting of theInternational Society for Magnetic Resonance in Medicine:21-27 April 2001；Glasgow,Scotland,UK,p.10(2001)

[43]Zhao,B.,Haldar,J.P.,Christodoulou,A.G.,Liang,Z.P.:Furtherdevelopment of image reconstruction from highly undersampled(k,t)-space datawith joint partial separability and sparsity constraints.In:Proceedings ofIEEE International Symposium on Biomedical Imaging:April 2011；Chicago,USA,pp.1593–1596(2011)

[44]Zhang,T.,Pauly,J.M.,Vasanawala,S.S.,Lustig,M.:Coil compressionfor accelerated imaging with Cartesian sampling.Magn Reson Med 69(2),571–582(2013)

[45]Walker,P.G.,Cranney,G.B.,Scheidegger,M.B.,Waseleski,G.,Pohost,G.M.,Yoganathan,A.P.:Semiautomated method for noise reduction and backgroundphase error correction in MR phase velocity data.J Magn Reson Imaging 3(3),521–530(1993)

[46]Stoica,P.,Selen,Y.:Model-order selection:a review of informationcriterion rules.IEEE Signal Process Mag 21(4),36–47(2004)

[47]Haldar,J.P.:Constrained imaging:denoising and sparse sampling.PhDthesis,University of Illinois at Urbana-Champaign,Electrical&ComputerEngineering Department(2011)

Claims (5)

1.一种实时相位对比磁共振成像(MRI)方法，包括：

采集实时相位对比MRI(PC-MRI)数据，其包括训练数据和成像数据；

基于采集的训练数据和成像数据，进行基于低秩模型的图像重建；

基于重建的实时血流图像，计算速度图；以及

基于计算出的速度图，进行血流量化分析。

2.如权利要求1所述的方法，其中，进行基于低秩模型的图像重建进一步包括：

针对训练数据进行时间插值；以及

针对插值后的训练数据进行时间子空间的估计。

3.如权利要求2所述的方法，其中，进行基于低秩模型的图像重建进一步包括：

针对采集的成像数据进行线圈灵敏度的估计；以及

基于估计的时间子空间、估计的线圈灵敏度和采集的成像数据，进行空间子空间估计。

4.如权利要求1所述的方法，其中，进行基于低秩模型的图像重建进一步包括：

将每一个血流图像序列表示为在空间和时间上分布的卡索拉蒂(Casorati)矩阵：

其中，ρ_v(r,t)代表动态图像，当v＝1时与血流补偿图像序列相关联，当v＝2,…,N_v时与血流编码图像序列相关联，

引入如下的联合卡索拉蒂矩阵：

使其符合低秩结构rank(C)≤L，其中L是指矩阵C的秩，并且满足L≤min(M,N)；

通过矩阵因子分解而利用显式的秩约束：

C＝UV，

其中且并且其中，U的列和V的行分别是支撑C的空间子空间和时间子空间；

通过解以下方程来进行对2D实时PC-MRI的空间子空间估计：

其中，d_i表示测量的数据，Ω表示稀疏采样运算符，F_s表示空间傅里叶变换矩阵，S_i表示第i个接收线圈的敏感度分布，i＝1,2,…,N_c；

在解后，联合卡索拉蒂矩阵可重建为：

继而得到每一个血流图像序列并估计血流速度，其中U/V/C是指真实量，而是指对应的重建结果。

5.如权利要求1所述的方法，其中，进行基于低秩模型的图像重建进一步包括：

将每一个血流图像序列表示为在空间和时间上分布的卡索拉蒂(Casorati)矩阵：

其中，ρ_v(r,t)代表动态图像，当v＝1时与血流补偿图像序列相关联，当v＝2,…,N_v时与血流编码图像序列相关联，

引入如下的联合卡索拉蒂矩阵：

使其符合低秩结构rank(C)≤L，其中L是指矩阵C的秩，并且满足L≤min(M,N)；

通过矩阵因子分解而利用显式的秩约束：

C＝UV，

其中且并且其中，U的列和V的行分别是支撑C的空间子空间和时间子空间；

通过解以下方程来进行对3D实时PC-MRI的空间子空间估计：

其中，d_i表示测量的数据，Ω表示稀疏采样运算符，F_s和F_t分别表示空间和时间傅里叶变换矩阵，S_i表示第i个接收线圈的敏感度分布，i＝1,2,…,N_c，λ是正则化参数，vec(·)运算符是将矩阵按列连成向量；

在解后，联合卡索拉蒂矩阵可重建为：

继而得到每一个血流图像序列并估计血流速度，其中U/V/C是指真实量，而是指对应的重建结果。