一种自适应高欠采样超极化气体肺部动态MRI重建方法转让专利
申请号 : CN201810993635.4
文献号 : CN109247939B
文献日 : 2020-06-19
发明人 : 周欣 , 肖洒 , 邓鹤 , 段曹辉 , 孙献平 , 叶朝辉
申请人 : 中国科学院武汉物理与数学研究所
摘要 :
本发明公开了一种自适应高欠采样超极化气体肺部动态MRI重建方法。首先根据已有的超极化气体肺部动态像的K空间数据优化各帧的K空间欠采样率；根据优化的K空间欠采样率，生成欠采样轨迹进行肺部吸气动态成像，获得欠采样K空间数据；然后利用欠采样K空间数据和重建目标函数重建出超极化气体肺部动态磁共振图像。本方法可以获得具有更高时间分辨率及更丰富细节信息的超极化气体肺部动态磁共振图像。
权利要求 :
1.一种自适应高欠采样超极化气体肺部动态MRI重建方法，其特征在于，包括以下步骤：
步骤1，获取已有的全采样超极化气体肺部动态像的M帧K空间数据d1，设置欠采样率SR的最小值TH；
步骤2、对第1帧K空间数据～第M帧K空间数据根据各自对应的第1帧空间欠采样率～第M帧K空间欠采样率进行欠采样获得第1帧欠采样K空间数据～第M帧欠采样K空间数据，对第
1帧欠采样K空间数据～第M帧欠采样K空间数据分别通过CS重建算法进行重建并评估重建误差，修改第k帧K空间数据对应的第k帧K空间欠采样率直至对应的重建误差小于或等于设定的重建误差的上确界，进而获得优化后的第k帧K空间欠采样率SRk，并根据优化后的第k帧K空间欠采样率SRk，确定优化后的第k帧K空间欠采样轨迹，其中k＝1,…,M；
步骤3，根据优化后的第1帧K空间欠采样轨迹～第M帧K空间欠采样轨迹进行肺部吸气动态成像，获得欠采样K空间数据d2；
步骤4，构建重建目标函数为：
其中，E为图像到欠采样K空间数据d2的变换矩阵，L为图像的低秩部分，S为图像的稀疏部分，d2为步骤3所获得的欠采样K空间数据，Ψ为计算相邻图像的差值的操作符，T为整体动态图像的稀疏变换矩阵，Φ为单帧动态图像的稀疏变换矩阵，为两个矩阵进行点乘的运算符，λG、λL、λS、λn分别为权衡动态图像连续性、图像低秩性、图像整体稀疏性和各帧图像稀疏性的正则化参数；
步骤5，根据步骤3所获得的欠采样K空间数据d2及步骤4所构建的目标函数，采用迭代软阈值算法重建图像。
2.根据权利要求1所述的一种自适应高欠采样超极化气体肺部动态MRI重建方法，其特征在于，所述的步骤2具体包括以下步骤：
步骤2.1、初始化第k帧K空间的欠采样率＝TH，其中k＝1,…,M；
步骤2.2、根据当前的第k帧K空间欠采样率生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，再根据概率密度矩阵生成第k帧K空间欠采样轨迹，利用第k帧K空间欠采样轨迹对第k帧K空间数据进行欠采样获得第k帧欠采样K空间数据；
步骤2.3、采用CS重建算法对第k帧欠采样K空间数据进行重建，通过计算本步骤获得的重建图像的平均绝对误差或均方误差或均方根误差评估重建误差；
步骤2.4、设置重建误差的上确界α，若重建误差小于或等于重建误差的上确界α，设置当前第k帧K空间欠采样率为优化后的第k帧K空间欠采样率，若重建误差大于重建误差的上确界α，当前的第k帧K空间欠采样率增加一个固定步长h返回步骤2.2；
步骤2.5、根据优化后的第k帧K空间欠采样率，确定优化后的第k帧K空间欠采样轨迹，具体为：以优化的第k帧K空间欠采样率生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，然后，根据概率密度矩阵生成第k帧K空间欠采样轨迹。
3.根据权利要求1所述的一种自适应高欠采样超极化气体肺部动态MRI重建方法，其特征在于，所述的步骤4中，λn∈{λ1～λM}，λ1+λ2+...+λM＝λv，λv为设定值，λ1SR1＝λ2SR2＝...＝λMSRM。
说明书 :
一种自适应高欠采样超极化气体肺部动态MRI重建方法
技术领域
[0001] 本发明涉及肺部磁共振成像(Magnetic resonance imaging,MRI)技术、欠采样和压缩感知(Compressed sensing,CS)理论等领域，具体涉及一种自适应高欠采样超极化气体肺部动态MRI重建方法。
背景技术
[0002] 超极化气体肺部MRI是一种新兴的肺部成像方法，能够对肺部通气功能(如气-气交换、气-血交换)进行有效评估[M.S.Albert et al.Nature,1994,370:199-201.]，在肺部重大疾病[如慢性阻塞性肺疾病(Chronic obstructive pulmonary disease，COPD)、哮喘等]早期诊断方面具有很大的潜力[F.C.Horn et al.Radiology,2017,284:854-861.]。
[0003] 常用的超极化气体肺部MRI是在屏气状态下对肺部进行成像，从而获得静态图像，如通气像，扩散加权像等。与这种静态成像方法相比，肺部超极化气体动态成像不仅能够动态观测超极化气体在肺部的流动/扩散过程，获取肺部在不同通气及扩张状态下的结构/生理信息，还有助于多尺度分析肺部的各级气管及肺泡的通气功能。这对存在肺部通气缺陷的肺部疾病(如COPD、哮喘、肺部炎症)的诊断和预后评估具有潜在的临床应用价值。
[0004] 由于超极化气体的核自旋极化度的衰减具有不可恢复性(随时间和激发次数呈指数衰减)，因此，成像速度制约着超极化气体肺部动态成像的发展。有研究者直接在吸气过程中对志愿者进行动态成像，包括采用radial采样轨迹及spiral采样轨迹进行成像。但是这些方法均具有较大伪影，图像信噪比较低。CS是一种不需更改MRI硬件系统就可以在不降低图像质量的情况下，将成像速度提高若干倍的快速成像技术。因此，也有研究者结合CS理论进行动态成像[S.Xiao et al.J.Magn.Reson.,2018,290:29-37.]，并在6.67s内获得15帧的超极化气体肺部动态图像，空间分辨率达到3mm。
[0005] 然而，目前已有的CS方法一般在动态成像过程中保持恒定的欠采样率，难以实现肺部气体流动/扩散的精确动态观测。当超极化气体进入肺部时，依次经过主气管，主支气管，次级支气管，细支气管，终末小支气管，呼吸性细支气管，肺泡小管，肺泡等结构。其中，主气管到细支气管的管径较大，气体流动阻力小，流速快，且结构简单，包含细节信息少(在稀疏变换域中具有高稀疏性)。根据CS理论，可以采用高倍欠采样(即低欠采样率)进行动态成像，从而提高其时间分辨率。而终末小支气管及呼吸性支气管的管径较小，气体流动阻力较大，流速较慢，并且这部分区域气管分支开始显著增多，细节信息增加(在稀疏变换域中稀疏性有所下降)。这时应当逐步降低欠采样倍数，从而实现对通气过程的有效观测(即空间分辨率有所提高)。当气体进入到肺泡小管和肺泡囊时，由于其直径极小，气体已经无法流动，只存在扩散运动[J.C.Woods et al.J.Magn.Reson.,2018,In press.]。但是，受限于超极化气体MRI的空间分辨率，气体在这部分区域的扩散对肺部细节信息的增加无显著影响，因此稳定的欠采样倍数可以满足通气功能有效观测。
[0006] 根据肺部的结构特征获得自适应不同欠采样率条件下的K空间数据后，需要有效的重建算法重建出高质量的动态图像。有研究者依据动态图像的低秩性，稀疏性和连续性(气体流入效应)能够在3倍恒定欠采样的情况下有效重建动态图像[S.Xiao et al.J.Magn.Reson.,2018,290:29-37.]。然而，这些方法无法满足自适应高欠采样K空间数据重建需求，导致重建结果具有较明显的伪影。
发明内容
[0007] 本发明的目的是针对现有技术存在的不足，提出一种自适应高欠采样超极化气体肺部动态MRI重建方法。通过分析重建误差，自适应地调整动态成像过程中所需的欠采样率，即在吸气过程早期，采用高倍欠采样采集K空间数据，在吸气过程中期，逐步降低欠采样倍数，在吸气过程末期，保持稳定的欠采样倍数。之后，考虑到肺部图像的稀疏性会随时间进行有规律的变化，在重建过程中针对该变化构建对应的目标函数，避免细节信息丢失或伪影、噪声难以消除等情况，从而获得高质量动态图像。
[0008] 为了实现上述目的，本发明采用了以下技术措施：
[0009] 一种自适应高欠采样超极化气体肺部动态MRI重建方法，包括以下步骤：
[0010] 步骤1，获取已有的全采样超极化气体肺部动态像的M帧K空间数据d1，设置欠采样率SR的最小值TH(一般TH≥10％)。其中，SR最小值TH的设置标准为目前已有的方法中能够准确从欠采样数据中重建超极化图像的最低欠采样率。M为进行一次全采样超极化肺部动态成像所获得动态图像总帧数。
[0011] 步骤2、对第1帧K空间数据～第M帧K空间数据根据各自对应的第1帧空间欠采样率～第M帧K空间欠采样率进行欠采样获得第1帧欠采样K空间数据～第M帧欠采样K空间数据，对第1帧欠采样K空间数据～第M帧欠采样K空间数据分别通过CS重建算法进行重建并评估重建误差，修改第k帧K空间数据对应的第k帧K空间欠采样率直至对应的重建误差小于或等于重建误差的上确界，获得第k帧K空间数据对应的优化后的第k帧K空间欠采样率SRk，并根据优化后的第k帧K空间欠采样率SRk，确定优化后的第k帧K空间欠采样轨迹，其中k＝1,…,M。
[0012] 步骤2具体包括以下步骤：
[0013] 步骤2.1、初始化第k帧K空间的欠采样率＝TH，其中k＝1,…,M；
[0014] 步骤2.2、根据当前的第k帧K空间欠采样率生成生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，再根据概率密度矩阵生成第k帧K空间欠采样轨迹，利用第k帧K空间欠采样轨迹对第k帧K空间数据进行欠采样获得第k帧欠采样K空间数据；
[0015] 步骤2.3、采用CS重建算法(包括非线性共轭梯度下降算法，迭代软阈值算法，交替方向乘子法等)对第k帧欠采样K空间数据进行重建，这些方法通过线性迭代优化，能够有效的从欠采样数据中重建图像。之后，通过计算本步骤获得的重建图像的平均绝对误差(Mean Absolute Error，MAE)或均方误差(Mean Square Error，MSE)或均方根误差(Root Mean Square Error，RMSE)评估重建误差，这些是评价重建质量时所常用的评价参数，通过定量评估两种图像之间的差异；
[0016] 步骤2.4、设置重建误差的上确界α(一般0≤α≤0.1)，并判断重建误差是否小于或等于重建误差的上确界α，若重建误差小于或等于重建误差的上确界α，设置当前第k帧K空间欠采样率为优化后的第k帧K空间欠采样率，若重建误差大于重建误差的上确界α，当前的第k帧K空间欠采样率增加一个固定步长h返回步骤2.2，理论上，增加固定步长会减小相应的重建误差；
[0017] 步骤2.5、根据优化后的第k帧K空间欠采样率，确定优化后的第k帧K空间欠采样轨迹。具体方法为：以优化的第k帧K空间欠采样率生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，然后，根据概率密度矩阵生成第k帧K空间欠采样轨迹。
[0018] 步骤3，根据优化后的第1帧K空间欠采样轨迹～第M帧K空间欠采样轨迹进行肺部吸气动态成像，获得欠采样K空间数据d2。
[0019] 步骤4，构建重建目标函数为：
[0020]
[0021] 其中，E为图像到欠采样K空间数据d2的变换矩阵，L为图像的低秩部分，S为图像的稀疏部分，d2为步骤3所获得的欠采样K空间数据，Ψ为计算相邻图像的差值的操作符，T为整体动态图像的稀疏变换矩阵，Φ为单帧动态图像的稀疏变换矩阵，ο为两个矩阵进行点乘的运算符，λG、λL、λS、λn分别为权衡动态图像连续性、图像低秩性、图像整体稀疏性和各帧图像稀疏性的正则化参数。其中，λG、λL和λS按照常用的迭代软阈值算法中方法来设置，目标该函数中最后一项所对应的正则化参数λ1，...，λM设置需要满足两个限定条件：(a)λ1+λ2+...+λM＝λv，λv为设定值，λv按照常用的迭代软阈值算法中方法来设置；(b)λ1SR1＝λ2SR2＝...＝λMSRM。目前已有重建方法考虑了动态图像的低秩性，稀疏性和连续性，但忽略了各帧图像稀疏性的差异性。由于获得各帧图像时超极化气体在肺部分布范围不同，稀疏性也将发生有规律的变化。然而，图像稀疏性与加权系数的乘积对重建结果具有较大影响。当乘积过大时，欠采样伪影和噪声等可以得到有效抑制，但同时会丢失图像的一些细节信息；反之，当乘积过小时，能够避免细节信息的丢失，但存在较明显的欠采样伪影和噪声。本方法通过合理设置各帧图像稀疏性的加权系数，能够避免细节信息丢失，同时有效抑制欠采样伪影和噪声。
[0022] 步骤5，根据步骤3所获得的欠采样K空间数据d2及步骤4所构建的目标函数，采用迭代软阈值(Iterative soft thresholding)算法重建图像，该算法具有较好的收敛性，较易得到步骤5目标函数中L及S的最优解。
[0023] 本发明相对于现有技术存在以下优点：该方法经过自适应调整欠采样率，能够在图像稀疏性较高的时间段进行高倍欠采样，显著增加成像速度；在图像稀疏性较低时降低欠采样倍数，能较好的保持图像细节信息。同时，通过合理设置每一帧图像稀疏性的加权系数保证图像质量，进而获得具有更低重建误差及更多细节信息的图像。
附图说明
[0024] 图1为实施例1中的步骤5获得的重建图像；其中，(a)～(f)分别为实施例1步骤5获得的重建图像中的第1～5帧重建图像及第15帧重建图像；
[0025] 图2为实施例1中的步骤5获得的重建图像的比对图；其中，(a)～(d)分别为实施例1中根据15帧K空间数据d1获得的第15帧的全采样图像、L+S重建图像、L+S+G重建图像，实施例1步骤5获得的重建图像中的第15帧重建图像；
[0026] 表1为实施例1中采用不同方法获得的动态图像与全采样图像之间的重建误差(以MAE为评价参数)。
具体实施方式
[0027] 为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
[0028] 实施例1：
[0029] 本实施例中超极化气体为129Xe。
[0030] 一种自适应高欠采样超极化气体肺部动态MRI重建方法，包括以下步骤：
[0031] 步骤1，获取已有的全采样超极化气体肺部动态像的M帧K空间数据d1，设置欠采样率SR的最小值TH＝10％，本实施例中，M为15。
[0032] 步骤2、对第1帧K空间数据～第15帧K空间数据根据各自对应的第1帧K空间欠采样率～第15帧K空间欠采样率进行欠采样获得第1帧欠采样K空间数据～第15帧欠采样K空间数据，对第1帧欠采样K空间数据～第15帧欠采样K空间数据分别通过CS重建算法进行重建并评估重建误差，修改第k帧K空间数据对应的第k帧K空间欠采样率直至对应的重建误差小于或等于重建误差的上确界，获得第k帧K空间数据对应的优化后的第k帧K空间欠采样率，并根据优化后的第k帧K空间欠采样率SRk，确定优化后的第k帧K空间欠采样轨迹，其中k＝1,…,M。得到的优化后的第1帧K空间欠采样率～第15帧K空间欠采样率分别为：0.10 0.10
0.10 0.10 0.12 0.20 0.28 0.30 0.30 0.32 0.34 0.34 0.34 0.34 0.34。
0.10 0.10 0.12 0.20 0.28 0.30 0.30 0.32 0.34 0.34 0.34 0.34 0.34。
[0033] 步骤2具体包括以下步骤：
[0034] 步骤2.1、初始化第k帧K空间欠采样率＝TH＝10％；
[0035] 步骤2.2、根据当前的第k帧K空间欠采样率生成生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，再根据概率密度矩阵生成第k帧K空间欠采样轨迹，利用第k帧K空间欠采样轨迹对第k帧K空间数据进行欠采样获得第k帧欠采样K空间数据；
[0036] 步骤2.3、采用CS重建算法中的非线性共轭梯度下降算法(还可采用迭代软阈值算法，交替方向乘子法等方法)对第k帧欠采样K空间数据进行重建，这些方法通过线性迭代优化，能够有效的从欠采样数据中重建图像。之后，通过计算本步骤获得的重建图像的平均绝对误差(Mean Absolute Error，MAE)评估重建误差(还可通过均方误差Mean Square Error，MSE或均方根误差Root Mean Square Error，RMSE进行评估重建误差)；
[0037] 步骤2.4、设置重建误差(MAE)的上确界α＝0.01，并判断重建误差是否小于或等于0.01，若是，设置当前第k帧K空间欠采样率为优化后的第k帧K空间欠采样率，若否，当前的第k帧K空间欠采样率增加一个固定步长h＝0.02返回步骤2.2，理论上，增加固定步长会减小相应的重建误差；
[0038] 步骤2.5、根据优化后的第k帧K空间欠采样率，确定优化后的第k帧K空间欠采样轨迹。具体方法为：以优化的第k帧K空间欠采样率生成概率密度从中心行向边缘行逐渐降低的概率密度矩阵，然后，根据概率密度矩阵生成第k帧K空间欠采样轨迹。
[0039] 步骤3，在1.5T MRI谱仪上，设置其它的成像参数：TR＝10.5ms，TE＝5.8ms，采样矩阵为64×64，成像视野大小为384×384mm2，定角激发(角度为7°)。根据步骤2得到的优化后的第1帧K空间欠采样轨迹～第15帧K空间欠采样轨迹进行肺部吸气动态成像，获得欠采样K空间数据d2。
[0040] 步骤4，构建重建目标函数为：
[0041]
[0042] 其中E为图像到欠采样K空间数据d2的变换矩阵，L为图像的低秩部分，S为图像的稀疏部分，d2为步骤3所获得的欠采样K空间数据，Ψ为计算相邻图像的差值的操作符，T为整体动态图像的稀疏变换矩阵，Φ为单帧动态图像的稀疏变换矩阵，ο为两个矩阵进行点乘的运算符，λG、λL、λS、λn分别为权衡动态图像连续性、图像低秩性、图像整体稀疏性和各帧图像稀疏性的正则化参数。其中，λG、λL和λS按照常用的迭代软阈值算法中方法设置为0.005、0.01、0.01，目标该函数中最后一项所对应的正则化参数λ1，...，λ15设置需要满足两个限定条件：(a)令λ1+λ2+...+λ15＝λv，λv按照常用的迭代软阈值算法中方法设置为0.01；(b)λ1SR1＝λ2SR2＝...＝λ15SR15。因此获得的λ1，...，λ15为1.0e-04×[0.2762，0.2762，0.2762，
0.2762，0.3315，0.5525，0.7735，0.8287，0.8287，0.8840，0.9392，0.9392，0.9392，0.9392，
0.9392]。
0.2762，0.3315，0.5525，0.7735，0.8287，0.8287，0.8840，0.9392，0.9392，0.9392，0.9392，
0.9392]。
[0043] 步骤5，根据步骤3所获得的欠采样K空间数据d2及步骤4所构建的目标函数，采用迭代软阈值(Iterative soft thresholding)算法重建图像，最终得到符合步骤5中目标函数的L和S，重建的图像为abs(L+S)，abs为求绝对值的运算符。
[0044] 在图1中显示了在1.5T MRI谱仪上经过步骤1～5处理后最终获得的人体肺部超极化气体129Xe动态成像第1～5帧及最后一帧的图像，可以有效的观测超极化气体在肺部流动\扩散的动态过程。同时，将现有的两种重建方法(L+S，L+S+G)所获得的重建结果与本方法的结果进行了比较。其中，L+S方法将动态图像视为低秩和稀疏两部分数据的和，通过迭代软阈值算法进行求解，在动态图像重建中具有较好效果[R.Otazo et al.Magn.Reson.Med.,2015,73:1125-1136.]；L+S+G方法在L+S方法的基础上引入气体流入效应，能有效抑制超极化气体肺部动态成像过程中的噪声和伪影[S.Xiao et
al.J.Magn.Reson.,2018,290:29-37.]。图2比较了本方法(自适应高欠采样)与L+S、L+S+G方法在3倍欠采样条件下的重建结果，可以看到本方法在更高欠采样倍数的情况下获得了更好的重建结果。表1为采用L+S、L+S+G和本方法所获得重建结果的重建误差(以MAE为评价参数)，可以看到与另两种方法相比，本方法能够显著降低欠采样所带来的重建误差。
al.J.Magn.Reson.,2018,290:29-37.]。图2比较了本方法(自适应高欠采样)与L+S、L+S+G方法在3倍欠采样条件下的重建结果，可以看到本方法在更高欠采样倍数的情况下获得了更好的重建结果。表1为采用L+S、L+S+G和本方法所获得重建结果的重建误差(以MAE为评价参数)，可以看到与另两种方法相比，本方法能够显著降低欠采样所带来的重建误差。
[0045] 表1
[0046]
[0047] 本文中所描述的具体实施方式仅仅是对本发明精神作举例说明。本发明中超极化气体不限于实施例中的129Xe，也包括3He、83Kr等。
[0048] 本发明所属技术领域的技术人员可以对所描述的具体实施方式做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。