CN106019189B - 一种基于自一致性的并行磁共振成像快速重构方法 - Google Patents

Abstract

本发明公开了一种基于自一致性的并行磁共振成像快速重构方法。本发明基于SPIRiT框架，针对含有JTV和JL1复合正则项的并行成像的重构问题，提出了一种笛卡尔采样的快速重构算法。解决现有技术并行磁共振成像技术虽然能够显著的减少扫描时间，但是现有的重构算法重构时间较长。该方法利用分裂Bregman技术（split Bregman）将原问题分解为多个易于求解的子问题。实验仿真结果表明，与现有的NLCG算法相比，新算法在保证重构质量的情况下，收敛速度有很大的提升。

Description

一种基于自一致性的并行磁共振成像快速重构方法

技术领域

本发明涉及一种基于自一致性的并行磁共振成像快速重构方法。

背景技术

磁共振成像(Magnetic Resonance Imaging,MRI)能够提供良好的人体软组织对比度，并且没有放射性，因此逐渐成为现代临床医学中不可缺少的成像工具。然而，受物理和生理因素的限制，采集磁共振信号的速度很慢。并行成像是一种提高采集速度的常用技术。SMASH和SENSE的提出，标志着并行成像成为切实可行的技术。并行成像使用对磁共振信号灵敏度不同的多个线圈同时采集磁共振信号。灵敏度信息的使用减少了用于重构的数据个数，从而提高了采集的速度。

过去十多年，产生了很多种并行成像重构算法。这些算法的不同之处在于：1)显式还是隐式使用灵敏度信息；2)恢复单一图像还是逐线圈恢复多幅图像。SMASH，SENSE，SPACE-RIP，kSPA属于采用显式灵敏度信息恢复单一图像的方法。PILS，GRAPPA，SPIRiT属于采用隐式灵敏度信息进行逐线圈恢复的方法。而AUTO-SMASH采用隐式灵敏度信息恢复单一图像，PARS则采用显式灵敏度信息进行逐线圈恢复。

SENSE是各种算法中应用最普遍的，它易于集成图像的先验知识。近年来，有研究人员将压缩感知与SENSE结合，增加正则项，取得了较好重构效果。当灵敏度已知时，SENSE能获得最优解。然而准确估计线圈的灵敏度往往非常困难。因此，不需显式使用线圈灵敏度信息的自校准方法就体现出了优势。

GRAPPA是一种应用广泛的逐线圈的自校准重构方法，由于不需要显式的使用灵敏度信息，从而避免了棘手的灵敏度估计。Lustig等在GRAPPA的基础上提出了一种更为有效的并行成像技术的理论框架——SPIRiT。和SENSE类似，SPIRiT将重构问题转化成一个逆问题来求解，易于集成图像的先验知识，而且能够重构任意k空间采样的数据。Murphy和Vasanawala等提出L1-SPIRiT模型，引入Joint L1正则项，并使用POCS和NLCG等算法来求解。然而POCS的重构质量欠佳，而NLCG算法的重构时间较长。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于自一致性的并行磁共振成像快速重构方法。并行磁共振成像技术能够显著的减少扫描时间，然而现有技术的重构算法重构时间较长。本发明基于SPIRiT框架，针对含有JTV和JL1复合正则项的并行成像的重构问题，提出了一种笛卡尔采样的快速重构算法。该算法利用分裂Bregman技术(splitBregman)将原问题分解为多个易于求解的子问题。实验仿真结果表明，与现有的NLCG算法相比，新算法在保证重构质量的情况下，收敛速度有很大的提升。

本发明的目的是通过以下技术方案来实现的：一种基于自一致性的并行磁共振成像快速重构方法，它包括以下步骤：

S0：初始化，令x⁰＝0，w⁰＝0，d⁰＝0，z⁰＝0，k＝0；

式中，为多线圈图像变量，每一列为一个按列堆栈的线圈图像N＝m×n，C为多线圈个数，m和n分别为单个线圈二维图像的行数和列数；k为循环变量；b_w、b_d和b_z为对偶变量，w、d和z为任意变量，w＝x，z＝Ψx，Ψ为逐线圈的小波变换， 其中D_n和D_m分别表示n×n和m×m的循环矩阵，所述的循环矩阵结构如下：

S1：分别计算w^k+1、d^k+1和z^k+1的值，其中：

式中，矩阵Γ＝λ(G-I)^T(G-I)+μI具有对角块结构，对该矩阵进行简单的重排序皆可产生块对角结构，由C×C块组成；μ为惩罚参数，λ为惩罚参数；

式中，μ₁表示惩罚参数，shrinkJ()表示联合一维收缩算子，计算公式如下：

式中，Ψ表示逐线圈的小波变换；repmat为matlab中的函数，功能为将第一个变量按照第二个参数的维度进行复制拓展，c表示循环索引；

式中，μ₂表示惩罚参数，shrink2J()表示联合二维收缩算子，计算公式如下：

式中，

S2：在循环变量c的值分别为1～C的情况下，计算

式中，为傅里叶变换矩阵，表示频率系数选择子，其每一行均从N×N的单位矩阵中抽出，M为欠采样的频率点数，M/N为欠采样率； 为采集到的频率系数，其每一列均为一个线圈图像的欠采样频率系数，α₁和α₂为正则化参数；

S3：分别对 和的值进行更新，其中：

S4：判断k值是否大于最大循环次数K，如果不是则对k值进行加1操作之后返回步骤S1，否则进入步骤S5；

S5：输出最终得到的x：

x＝x^k+1；

S6：针对每次循环重构处的各个线圈的图像数据使用SRSOS方法对各个线圈图像进行联合，得到最终的单幅重构图像：

本发明的有益效果是：本发明公布了一种基于自一致性的并行磁共振成像快速重构方法。并行磁共振成像技术能够显著的减少扫描时间，然而现有的重构算法重构时间较长。本发明基于SPIRiT框架，针对含有JTV和JL1复合正则项的并行成像的重构问题，提出了一种笛卡尔采样的快速重构算法。该算法利用分裂Bregman技术(split Bregman)将原问题分解为多个易于求解的子问题。实验仿真结果表明，与现有的NLCG算法相比，新算法在保证重构质量的情况下，收敛速度有很大的提升。

附图说明

图1为本发明方法流程图；

图2为T1加权的使用SPGR序列获取的全采样的高分辨率大脑成像数据(即data1)图；

图3为使用2维旋转回城序列获取的全采样的轴向脑成像数据(即data2)图；

图4为加速比约为6的泊松亚采样示意图；

图5为亚采样率为6的采用data1测试序列时本发明与NLCG性能比较；

图6为亚采样率为6的采用data2测试序列时本发明与NLCG性能比较；

图7为测试序列为data1时NLCG算法重构图像与原始图像的差值图像；

图8为测试序列为data1时本发明方法重构图像与原始图像的差值图像；

图9为测试序列为data2时NLCG算法重构图像与原始图像的差值图像；

图10为测试序列为data2时本发明方法重构图像与原始图像的差值图像。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案：

本发明是基于SPIRiT框架提出的一种高效的重构方法。

SPIRiT是一个基于自校准的并行成像重构方法，该方法逐线圈内插出亚采样时缺失的频率点。在SPIRiT中，一个内插核g_ij通过对频域中心全采样的数据(通常称为自校准信号)进行校准得到。若以x_i表示第i个线圈的整个频域数据，则基于校准的一致性准则可写成：

式中，C表示线圈数，“*”表示卷积操作。卷积核g_ij被称作SPIRiT核。

全部线圈的一致性准则可简化成矩阵形式：

x＝Gx (2)

式中，x为全部线圈的图像数据。G为内插操作子，能够作为频域的卷积或图像域的乘积来实现。

当然，除了基于校准的一致性，重构还需要与原采样的数据保持一致性，表示成矩阵形式：

其中为多线圈图像变量，每一列为一个按列堆栈的线圈图像N＝m×n，m和n分别为单个线圈二维图像的行数和列数。为傅里叶变换矩阵，表示频率系数选择子，其每一行均从N×N的单位矩阵中抽出，M为欠采样的频率点数，M/N为欠采样率； 为采集到的频率系数，其每一列均为一个线圈图像的欠采样频率系数。

由于多线圈图像在小波变换域和一阶有限差分变换域的相似性，因此我们同时引入Joint L1(JL1)正则项和Joint Total Variation(JTV)正则项。根据式(2)和式(3)的两个限制，并行成像的重构问题可转化成以下的优化问题：

其中α₁和α₂为正则化参数，Ψ表示逐线圈的小波变换，||x||_JTV按下式进行定义：

其中，D_hx和D_vx分别为对x进行行方向和列方向的一阶有限差分变换， D_m表示m×m的循环矩阵，结构如下：

为了推导方便，将行、列方向的一阶有限差分变换合并写作：即因此，x的一阶有限差分变换可简单表示成D_hvx。很明显有：

为了方便推导，将||x||_JTV写成以下形式：

因此，式(4)所示的优化问题可以写成：

为了求解式(8)所示的优化问题，应用二次罚函数技术，则可得：

其中λ为惩罚参数。

为了求解式(9)所示的优化问题，引入任意变量w＝x，z＝Ψx，d＝Dx，d有以下形式：

利用分裂Bregman技术，原问题式(9)可以通过以下迭代问题求解：

式(11)中，μ，μ₁和μ₂均为惩罚参数。对式(11)分别求解x,w,z,d，则得到以下子问题的求解：

式(15)的最优化条件为而矩阵

Γ＝λ(G-I)^T(G-I)+μI具有对角块(diagonal-block)结构，对该矩阵进行简单的重排序皆可产生块对角(block-diagonal)结构，由C×C块组成。Γ的逆出现在w的更新表达式中：

由于C比较小，因此可直接对Γ中的每个C×C块求逆来求得Γ^-1。

式(16)中，z的解由下式给出：

式中，μ₁表示惩罚参数，shrinkJ()表示联合一维收缩算子(联合一维软阈值法)，计算公式如下：

式中，repmat为matlab中的函数，功能为将第一个变量按照第二个参数的维度进行复制拓展。

式(17)中，d的解由下式给出：

式中，μ₂表示惩罚参数，shrink2J()表示联合二维收缩算子(联合二维软阈值法)，计算公式如下：

式中，

由于式(18)中各个线圈变量独立，故可以分离成单个线圈进行计算，即：

根据最优化条件，式(24)的解可通过下式得到：

由于Ψ为正交小波，故有Ψ^TΨ＝I；可被FFT对角化，因此该式可通过FFT快速求解：

既然所有子问题都能够高效求解，针对问题(9)的求解，本发明得到一种高效的重构算法——基于自一致性的含联合全变分的并行磁共振成像快速分裂Bregman算法(SplitBregman for SPIRiT based Parallel MR Imaging with JTV and JL1regularization,SB4SpMRI)，SB4SpMRI算法的具体流程如图1所示。在该流程中，C为线圈总数，c为相应的循环变量(即前文所述的线圈索引)；K为总迭代次数，k为相应的循环变量(即第k次迭代)。

算法1中，第2步为循环控制条件，当迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某个值时，停止循环。算法1中，每次循环重构出的是各线圈的图像数据，因此还需要使用SRSOS(Square Root Sum Of Squares)方法对各个线圈图像进行联合，才可得到最终的单幅重构图像。这个后处理过程可用SRSOS(x)表示，计算方法如式(27)所示：

本文所述算法使用matlab编程环境(Version R2008b,the MathWorks Inc.,Natick,MA)实现，并应用在亚采样的并行MR成像重构中。所有实验在配置为Intel Corei53230M@2.6GHz的CPU，4GB内存和Windows 764bit操作系统的笔记本上执行。

实验中所使用的素材有两组，都是真实的多线圈并行MR成像数据，如图2和3所示：

第一组为T1加权的使用SPGR序列获取的全采样的高分辨率大脑的成像数据(data1)。使用有8通道接收线圈的GE Signa-Excite 1.5T扫描系统扫描完成。扫描参数为：TE和TR分别为8ms和17.6ms，flip等于20°，BW＝6.94kHz。FOV为20×20×20cm，成像矩阵大小为200×200×200。

第二组是全采样的轴向脑成像数据(data2)。该成像数据使用2维旋转回声序列(2D spin echo sequence)在有8个线圈的GE 3T扫描设备(GE Healthcare，Waukesha，WI)上采样得到，TE和TR分别为11ms和700ms，矩阵尺寸为256×256，FOV为220×220mm。

本发明采用已广泛应用的泊松亚采样图案对全采样数据进行亚采样，从而得到不同加速倍数的亚采样数据。亚采样图案中，用于SPIRiT核自校准的中心全采样信号尺寸为24×24。图4为加速比约为6的采样图案的示意图(已考虑中心全采样信号，矩阵尺寸为256×256)。

此处使用SNR来定量评价重构图像的质量，SNR定义为：

式中，MSE表示参考图像与重构图像间的均方误差，Var表示参考图像的方差。

为了测试本发明的新算法的性能，我们将新算法(SB4SpMRI)与经典算法——非线性共轭梯度的算法(NLCG)进行比较，正则项均采用JTV和JL1复合正则项。实验中SPIRiT核大小均设置成5×5。对各个算法的正则化参数进行调优，以使算法达到最优性能。

图5和图6给出了亚采样率为6时新算法SB4SpMRI与NLCG的性能比较。当某种算法重构图像的SNR变化趋于稳定时，就认为该算法已经收敛。图5和图6分别为采用数据序列data1和data2时的实验结果。从图5可以看出，新算法SB4SpMRI能达到与NLCG算法相当的SNR，而收敛速度却大大的快于NLCG。而对于图6，也可以得到类似的结论。

图7和图8分别给出了测试序列为data1时，NLCG和SB4SpMRI算法重构图像与原始图像的差值图像。图9和图10分别给出了测试序列为data2时，NLCG和SB4SpMRI算法重构图像与原始图像的差值图像。需要指出的是，这里比较的算法的重构图像都是经过充分迭代之后得到的。从图7、图8、图9和图10可以看出，SB4SpMRI算法的重构误差与NLCG相当。

Claims (1)

1.一种基于自一致性的并行磁共振成像快速重构方法，其特征在于：它包括以下步骤：

S0：初始化，令x⁰＝0，w⁰＝0，d⁰＝0，z⁰＝0，k＝0；

式中，为多线圈图像变量，每一列为一个按列堆栈的线圈图像N＝m×n，C为多线圈个数，m和n分别为单个线圈二维图像的行数和列数；k为循环变量；b_w、b_d和b_z为对偶变量，w、d和z为辅助变量，w＝x，z＝Ψx，Ψ为逐线圈的小波变换，其中D_n和D_m分别表示n×n和m×m的循环矩阵，I_m和I_n分别为m×m和n×n的单位矩阵，所述的循环矩阵结构如下：

S1：分别计算w^k+1、d^k+1和z^k+1的值，其中：

式中，矩阵Γ＝λ(G-I)^T(G-I)+μI具有对角块结构，对该矩阵进行简单的重排序皆可产生块对角结构，由C×C块组成；其中，G为内插操作子，μ为惩罚参数，λ为惩罚参数；

式中，μ₁表示惩罚参数，shrinkJ()表示联合一维收缩算子，计算公式如下：

式中，下标c表示原变量的第c个线圈数据，Ψ表示逐线圈的小波变换；repmat为matlab中的函数，功能为将第一个变量按照第二个参数的维度进行复制拓展，c表示循环索引；

式中，μ₂表示惩罚参数，shrink2J()表示联合二维收缩算子，计算公式如下：

式中，其中，d_h＝D_hx，d_v＝D_vx；

S2：在循环变量c的值分别为1～C的情况下，计算x_c ^k+1：

式中，为傅里叶变换矩阵，表示频率系数选择子，其每一行均从N×N的单位矩阵中抽出，M为欠采样的频率点数，M/N为欠采样率；为采集到的频率系数，其每一列均为一个线圈图像的欠采样频率系数，α₁和α₂为正则化参数；

S3：分别对和的值进行更新，其中：

S4：判断k值是否大于最大循环次数K，如果不是则对k值进行加1操作之后返回步骤S1，否则进入步骤S5；

S5：输出最终得到的x：

x＝x^k+1；

S6：针对每次循环重构处的各个线圈的图像数据使用SRSOS方法对各个线圈图像进行联合，得到最终的单幅重构图像：