一种肺部磁共振动态成像的方法

技术领域

本发明涉及磁共振成像领域，尤其涉及一种肺部磁共振动态成像的方法。

背景技术

随着当代技术的发展，磁共振成像(Magnetic Resonance Imaging,MRI)技术已经展示出检查诸如肺栓塞、肺结节和新型冠状肺炎等肺部疾病的能力。与计算机断层扫描(Computed Tomography,CT)技术相比，肺部磁共振可以在无电离辐射的情况下提供良好的软组织对比度和肺部解剖结构信息，更加适用于儿童和其他需要长期随访检查的患者。然而，目前大部分肺部磁共振成像技术是借助呼吸门控来采集成像数据，这种方法扫描效率低，且只能获得单一呼吸状态的静态图像。动态磁共振成像是一项很有前景的技术，这种技术在进行肺部结构检查时扫描效率高，同时借助动态图像可以评估肺部通气功能。然而，由于肺部质子密度低、信号衰减快(T2*短)、呼吸运动以及呼吸过程中肺实质信号强度的改变等综合因素，使得动态肺部MRI具有一定难度。

为了解决肺部质子密度低、信号衰减快的问题，在肺部磁共振成像时一般选用超短回波时间(Ultrashort Time ofEcho,UTE)序列采集信号。UTE序列在射频激发后立即施加读出梯度并开始采集K空间信号，因此能够以亚毫秒级别的回波时间(echo time,TE)捕捉肺实质的短T2*信号。近期的一些研究表明，使用三维径向UTE序列可以检测一些肺部疾病，为CT扫描提供了一个可能的替代方案。尽管借助UTE序列较好解决了肺部磁共振信号采集的问题，肺部UTE成像依然需要准确的呼吸监测，达到减少图像中的呼吸伪影的问题。

在MR数据采集过程中，监测呼吸运动的方法有许多种。临床上最常用的方法是借助外部设备，如呼吸绑带，该方法通过间接测量腹部拉伸来估计呼吸运动。然而，这种方法操作过程繁琐，并且估计结果不太准确。另外一种方法是自导航(self navigation)，该方法是从采集到的K空间数据中直接估计呼吸信号，不需要额外的呼吸监测设备。在这种方法中，自导航信号可以是K空间中心点的数据，也可以是沿上下方向的导航回波(Superior-Inferior navigator，SI navigator)。此外，获取高时间分辨率、低空间分辨率的图像序列也可用于自导航。

利用检测到的呼吸信号，有多种运动校正或补偿重建方法可以用来重建出动态图像。目前最常用的方法是运动分解(motion resolved)重建，在该方法中，所有采集到的数据被分入不同的运动状态，并使用基于图像序列时空相关性的压缩感知(CompressedSensing,CS)算法重建出动态图像。这种方法已成功应用于腹部成像、心脏电影成像、动态对比度增强MRI和动态肺部UTE成像，代表性的有XD-GRASP(eXtra-Dimensional Golden-angle Radial Sparse Parallel)和XD-UTE(eXtra-Dimensional UTE)。然而，在运动位移较大的情况下，这种方法的重建结果容易产生模糊伪影。为了减少大的运动位移对运动分解重建的影响，一些研究提出了基于运动补偿(Motion Compensation,MoCo)的重建方法，并已成功应用于动态心脏MRI。近期，这些MoCo重建方法的一个变种，称为iMoCo(iterativeMoCo)UTE，已成功应用于肺部UTE MRI。在iMoCo中，使用配准得到的运动场(motion field)将参考状态下的单帧图像转换到其他运动状态来模拟呼吸运动，重建得到参考状态下的静态图像，而不是动态图像序列。然而，由于基于CS的运动分解重建和MoCo重建方法没有充分考虑肺实质信号强度与肺部磁共振动态图像的时间域相关性的变化关系，导致重建图像质量降低。

因此，本领域的技术人员致力于开发一种肺部磁共振动态成像的方法，利用运动补偿增加动态图像序列中图像间的时间相关性，减少重建伪影，并应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高重建质量。

发明内容

与人体其他部位的动态成像相比，肺部磁共振动态成像具有如下特征：在呼吸气末端(End Expiration,EE)到吸气末端(End Inspiration,EI)的过程中，肺泡密度会降低，从而导致动态图像序列中肺实质区域的信号强度呈现下降趋势。因此，肺部磁共振动态图像的时间域相关性会随着肺实质信号强度的变化而降低。然而，现有已知的基于运动校正或补充实现肺部磁共振动态成像的方法(例如基于图像序列时空相关性的压缩感知(Compressed Sensing,CS)的重建方法，以及基于iMoCo(iterative MoCo)UTE的运动补偿重建方法)没有充分考虑到这种信号强度的变化，导致重建图像质量降低。

本发明提供了一种肺部磁共振动态成像方法，通过运动补偿重建算法来实现动态成像，包括三维径向UTE采集，数据采集过程中的呼吸状态监测，运动分解重建和运动场估计，以及一种新型动态MRI重建算法。该算法主要是利用运动补偿来增加动态图像序列中图像间的时间相关性，减少重建伪影；并应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高重建质量。

为实现上述目的，本发明提供了一种肺部磁共振动态成像方法，包括：

步骤S10采用超短回波时间序列，采集自由呼吸状态下的肺部磁共振信号；

步骤S20获取呼吸曲线；

步骤S30利用所述呼吸曲线，将所述步骤S10中采集的数据进行运动分解重建，得到肺部运动分解图像；

步骤S40将所述肺部运动分解图像进行运动场估计；

步骤S50将所述步骤S40得到的结果进行运动状态加权的运动补偿重建，得到肺部磁共振动态图像；

步骤S60将所述肺部磁共振动态图像进行通气量图估计，得到肺部的通气量。

进一步地，在所述步骤S10中，所述超短回波时间序列的采样轨迹为交错式的螺旋叶序形式。

进一步地，所述步骤S10包括：在采样过程中，将每个叶序相较于前一个叶序的Z轴旋转一个黄金角度137.51°。

进一步地，所述步骤S10还包括：在数据采集过程中，监测呼吸状态；其中，对所述呼吸状态的监测可以使用自导航方式，在每个交错叶序采集之前，先采集一条沿着导航方向的回波信号用作导航信号；或者，对所述呼吸状态的监测也可采用外部呼吸监测设备，比如使用呼吸绑带，记录采集过程中的呼吸状态。

进一步地，所述步骤S20包括：

使用所述自导航方式监测呼吸状态时，对所述导航信号进行一维傅里叶反变换，得到所述肺部磁共振信号沿导航方向的各线圈通道的投影分布图；

从所述各线圈通道的所述投影分布图中估计出呼吸信号；

从所述呼吸信号中选出真实可信的呼吸信号；

将所述真实可信的呼吸信号的平均值作为最终估计出的呼吸曲线；或者

使用所述外部呼吸监测设备监测所述呼吸状态时，从所述外部呼吸监测设备直接获取记录到的所述呼吸曲线。

进一步地，所述步骤S30包括：

基于图像序列时空相关性的压缩感知算法重建出所述肺部运动分解图像。

进一步地，所述基于图像序列时空相关性的压缩感知算法的重建用如下公式表示：

其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，TV

进一步地，利用共轭梯度法求解所述公式1。

进一步地，所述步骤S40包括：

根据所述肺部运动分解图像，将每个运动状态的图像配准到其他运动状态的图像，得到所述每个运动状态到其他所有运动状态的运动场。

进一步地，基于Demons的三维非刚性图像配准技术进行图像配准。

进一步地，所述步骤S50包括：

利用如下公式进行运动补偿重建：

其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，λ

进一步地，所述Ope()的表达式如下：

其中，N是所述待求解的肺部运动分解图像X的总的运动状态数，x

进一步地，所述运动状态加权系数的表达式如下：

其中，c表示衰减系数。

进一步地，利用交错方向乘子法求解所述公式2。

进一步地，重复执行所述步骤S40和所述步骤S50。

进一步地，重复执行所述步骤S40和所述步骤S50包括：

将所述步骤S40得到的结果进行所述运动状态加权的运动补偿重建后，得到肺部磁共振动态图像序列；

借助质量更高的图像，将所述肺部磁共振动态图像序列重复执行所述步骤S40；

重复以上两个步骤直至达到迭代停止条件。

进一步地，所述步骤S60包括：

基于雅各比行列式估计所述通气量。

进一步地，所述步骤S60包括：

将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；

计算每个所述运动场的雅各比行列式，以获得从所述中间运动状态到所述每个运动状态的肺部体积的百分比变化；

将呼吸周期内的所述雅各比行列式的曲线拟合到所述肺部体积的演化函数，然后估计所述通气量，其中，所述演化函数如下：

其中，v

进一步地，所述步骤S60包括：

根据肺实质信号强度的时域变化估计所述通气量。

进一步地，所述步骤S60包括：

将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；

将呼吸周期内每个像素的信号强度值随时间变化的值拟合到肺部密度的演化函数，然后估计所述通气量，其中，所述演化函数如下：

其中，i

本发明提供的肺部磁共振动态成像方法，具有以下有益的技术效果：

1、利用SI导航的UTE序列采集肺部磁共振信号，无需外部呼吸检测设备的肺部自由呼吸扫描；

2、以及基于交错式螺旋叶序形式采集肺部磁共振信号，有助于减少动态图像重建时的伪影；

3、利用运动补偿来增加动态图像序列中图像间的时间相关性，减少了重建伪影；

4、应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高了重建质量；

5、本申请在欠采样下能够获取高质量肺图像和通气量定量图像，缩短了扫描时间；且能提供肺部结构与通气功能的三维可视化。

以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。

附图说明

图1是本发明一个较佳实施例提供的肺部磁共振动态成像方法的流程图；

图2是本发明一个较佳实施例提供的带SI导航功能的UTE序列的数据采集方式图示；

图3是本发明一个较佳实施例提供的从SI导航信号中提取呼吸曲线的流程示意图；

图4是本发明一个较佳实施例提供的呼吸曲线图；

图5是本发明一个较佳实施例提供的运动分解重建示意图；

图6是本发明一个较佳实施例提供的运动场估计方式的图示；

图7是本发明一个较佳实施例提供的基于运动状态加权的运动补偿重建算法的框架图；

图8是图7的流程图；

图9是本发明一个较佳实施例提供的估计通气量的示意图；

图10是本发明一个较佳实施例提供的离体猪肺成像测试结果；

图11是本发明一个较佳实施例提供的离体猪肺通气量图估计结果；

图12是本发明一个较佳实施例提供的人体肺部自由呼吸动态成像测试结果；

图13是本发明一个较佳实施例提供的人肺通气量图估计结果。

具体实施方式

以下参考说明书附图介绍本发明的多个优选实施例，使其技术内容更加清楚和便于理解。本发明可以通过许多不同形式的实施例来得以体现，本发明的保护范围并非仅限于文中提到的实施例。

如图1所示，本发明提供的肺部磁共振动态成像方法，包括如下步骤：

S10、数据采集：可以利用带导航功能的超短回波时间(UTE)序列来进行采集。例如，利用采用带上下方向(SI)导航功能的UTE序列，在自由呼吸状态下采集肺部磁共振信号。其中UTE序列的采样轨迹为交错式的螺旋叶序形式，如图2所示，在采样过程中，采样方式为：将每个叶序相较于前一个叶序的Z轴旋转了一个黄金角度137.51°，黄金角度(GoldenAngle,GA)的定义为

S20、呼吸信号获取：借助步骤S10中的SI导航信号，估计采集过程中的呼吸曲线。如图3所示，对步骤S10所得各线圈通道SI导航信号(图3中的第一列)沿Kz方向进行一维傅里叶反变换，得到每个线圈通道中人体信号(即所采集的肺部磁共振信号)沿着SI导航方向的投影分布图(图3中的第二列)。采用主成分分析的算法从各线圈通道的投影分布图中估计出呼吸信号，具体做法是对各线圈通道的投影分布图进行特征值分解，将最大特征值对应的时域特征向量作为各线圈通道的呼吸信号(图3中的第三列)。在得到从各线圈通道中估计出的呼吸信号后，使用线圈聚类(coil clustering)的方法从所有的呼吸信号中选出真实可信的一部分，线圈聚类的具体做法是将所有线圈通道的呼吸曲线两两组合，计算两条曲线间的相关性，相关系数高于0.95的两条曲线就作为真实可信的呼吸曲线。如图3中第三列所示，线圈1和n的呼吸曲线相关性高，而线圈1和N的相关性低，线圈1和n的呼吸曲线就是真实可信的。最后，将这部分真实可信的呼吸曲线的平均值作为最终估计出的呼吸曲线(图3中的第四列)。图4为从某次扫描得到的SI导航信号中提取出的呼吸曲线。

在一些实施方式中，当使用外部呼吸监测设备监测呼吸状态时，从该外部呼吸监测设备直接获取记录的呼吸曲线。

S30、运动分解重建：借助步骤S20所得图4中的呼吸曲线，将步骤S10中采集到的成像数据按呼吸位移大小分入N个不同的运动状态，将同一个呼吸状态下采集到的数据分入同一个运动状态，这样共得到N个不同运动状态下的K空间数据，如图5第二列所示。此时，由于各个运动状态的K空间数据是高度欠采样的，直接对其进行非均匀快速傅里叶逆变换(inverse Non Uniform Fast Fourier Transform,iNUFFT)得到的图像有很严重的欠采样伪影，如图5第三列所示欠采样伪影图像。因此，这里使用基于图像序列时空相关性的压缩感知算法重建出运动分解图像。该重建可以用公式表示为如下优化问题：

其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道UTE序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换(NUFFT)算子。TV

S40、运动场估计：如图6所示，根据步骤S30所得运动分解图像，将每个运动状态的图像配准到其它(N-1)个运动状态的图像，得到每个运动状态到其它所有运动状态共计N×(N-1)个运动场，用于之后的图像重建。此步骤进行图像配准时可选用基于Demons的三维非刚性图像配准技术(基于Demons的三维非刚性图像配准技术可以采用现有技术，例如参考文献“Thirion JP.Image matching as a diffusion process:an analogy withMaxwell's demons.Med Image Anal.1998；2:243-60”)，也可选用其它常用的三维非刚性图像配准方法。

S50、运动状态加权的运动补偿重建，得到动态图像：如图7所示，使用步骤S30所得分好运动状态的成像数据与步骤S40所得运动场，进行运动补偿重建，得到肺部磁共振动态图像。运动补偿重建可以用公式表示为如下优化问题：

其中，以上公式的数据一致性项(第一项)中的X，d，D，S，F，空间域稀疏项(第二项)中的λ

其中N是X总的运动状态数，x

对于某一运动状态k，系数α

在一些较佳的实施方式中，可以通过重复迭代的方式提升图像重建质量。如图7和图8所示，重复执行步骤S40和步骤S50。具体地，步骤S50可以进一步包括：S501：运动状态加权的运动补偿重建，即利用公式2、3、4得到动态图像：S502；在得到动态图像序列后，借助质量更高的该图像重新估计运动场(执行步骤S40)，如S503所示。通过将该重建与运动场估计这两个步骤迭代执行，可以增加运动场估计精度并提升图像重建质量。迭代停止的条件可以根据实际需求设定，例如，设定具体的次数，或者设定图像质量达到预设条件。

S60、通气量图估计：根据步骤S50所得肺部动态图像，估计肺部通气量图。在呼吸过程中，肺泡的体积会发生改变，从动态图像可以估计出这种改变导致的运动场，运动场的雅各比行列式代表肺泡体积的百分比变化(参见图9)，用这种方式可以估计出肺部的通气量值；此外，呼吸过程中肺泡体积的改变也导致肺泡密度的改变，反映为磁共振图像中肺实质信号强度的时域变化，借助这种信号强度变化的程度大小也可以估计出肺部的通气量值，如图9所示。

基于雅各比行列式估计通气量的方法：如图9所示，将动态图像配准到一个选取的中间运动状态(如运动状态N/2)，得到每个运动状态到中间状态的运动场，其中，中间运动状态可以是第1个运动状态到第N个运动状态中的任意一个，较佳地，选取第N/2个运动状态。然后计算每个运动场的雅各比行列式(Jacobian determinant)，以获得从中间状态到每个运动状态肺部体积的百分比变化，雅各比行列式的计算方式如下(公式5)：

其中(u

其中，v

基于肺部图像信号强度估计通气量的方法：与基于雅各比行列式的方法类似，如图9所示，将动态图像配准到一个中间运动状态，然后将呼吸周期内每个像素的信号强度值随时间变化的值拟合到肺部密度也即肺部磁共振信号强度的演化函数(公式7)：

其中，i

本发明较佳实施例在3.0特斯拉的医用磁共振成像系统(uMR790,上海联影医疗)上进行了离体猪肺和人肺自由呼吸成像测试，证实了本发明前述较佳实施例进行肺部动态成像的可行性。

图10为离体猪肺成像测试结果。为了得到离体猪肺在不同呼吸状态下的全采样数据，用打气筒向猪肺中打入不同气量的空气：0毫升、200毫升、400毫升、600毫升、800毫升，然后对各个状态的猪肺做静态扫描，得到全采样数据。扫描参数如下：视场(field ofview)＝32×14×32cm

图10中，第一行是从全采样的K空间数据中重建出的作为参考的真实图像(GroundTruth)，第二行是第一行图像中方框区域的放大图。第三行是基于运动分解重建(XD-UTE)得到的作为参考的真实图像，第四行是第三行图像的误差图。第五行是基于iMoCo得到的作为参考的真实图像，第六行是第五行的误差图。第七行是基于本发明提供的方法得到的真实图像，第八行是第七行的误差图。从呼气末端到吸气末端，由于肺泡密度降低，可以看到肺实质的信号强度明显下降。在XD-UTE图像(第三行)中，一些气管结构和猪肺边缘存在图像模糊，在误差图(第四行)中用绿色箭头表示。在iMoCo的结果(第五行)中，这种图像模糊大大减少。然而，在iMoCo重建的图像中，肺实质信号强度的时域变化大大减少，在iMoCo的误差图(第六行)中可以看到由此导致的明显的重建误差。总的来说，与iMoCo相比，使用本发明的重建算法能够更好地保留了肺实质信号强度的时域变化，并且与XD-UTE相比，随机噪声减少，气管结构更清晰，肺部边缘更清晰。

图11是借助图10中动态猪肺图像，用基于雅各比行列式的方法估计出的猪肺通气量图。第一行分别显示了借助真实图像(Ground Truth)、XD-UTE、iMoCo和本发明重建算法图像估计出的通气量图。第二行显示由不同图像得出的通气量图与Ground Truth之间的误差。与另外两种方法相比，本发明重建算法的通气量图估计误差更小，如红色箭头所示。

图12为人肺自由呼吸成像测试结果。人体肺部成像扫描参数如下：视场(field ofview)＝40×28×30cm

图12中，第一行是对高度欠采样的不同运动状态数据直接进行NUFFT变换得到的图像。其他几行是NUFFT、XD-UTE、iMoCo和使用本重建算法的结果放大图，放大区域如第一行中的方框所示。在NUFFT图像中，存在严重的欠采样伪影，这些伪影在XD-UTE图像中有所减少。然而，在XD-UTE图像中，一些肺部结构和膈肌区域存在图像模糊，如绿色箭头所示。在iMoCo图像中，这些残留噪声和图像模糊进一步减少，但iMoCo的图像看起来过于平滑，特别是位于红圈区域的一些肺实质结构模糊不清。相较于iMoCo，使用本发明重建算法的结果中肺实质的对比度更高，肺部结构和膈肌的清晰度与iMoCo的相当。

图13是借助图12中肺部动态图像，采用基于雅各比行列式的方法估计出的人肺通气量图。与图11中展示的猪肺通气量图评估结果一致，使用本发明重建算法的人肺通气量图要比XD-UTE更加真实。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。