CN105184755B

CN105184755B - 基于自一致性的含联合全变分的并行磁共振成像高质量重构方法

Info

Publication number: CN105184755B
Application number: CN201510676797.1A
Authority: CN
Inventors: 段继忠; 罗仁泽; 苏赋; 邓魁; 郑勉; 汪敏; 曹玉英
Original assignee: Southwest Petroleum University
Current assignee: Southwest Petroleum University
Priority date: 2015-10-16
Filing date: 2015-10-16
Publication date: 2017-12-05
Anticipated expiration: 2035-10-16
Also published as: CN105184755A

Abstract

本发明公开了一种基于自一致性的含联合全变分的并行磁共振成像高质量重构方法，本发明基于SPIRiT框架，针对含有JTV和JL1复合正则项的并行成像的重构问题，提出了一种高质量的重构算法。首先将有约束的重构问题，转化为无约束的最优化问题，再对数据保真项和自校正项进行简化，再采用算子分裂技术将简化后的重构问题转化成一个梯度计算问题和一个含有JTV和JL1复合正则项的去噪问题，复合正则项的去噪问题通过全新设计的基于Split Bregman技术的算法求解。最后再通过FISTA进行加速。本发明设计实验，比较了新算法与其他常用算法的重构性能。实验仿真表明，新算法的收敛速度与POCS算法相当，而重构图像的SNR有较大提升。

Description

基于自一致性的含联合全变分的并行磁共振成像高质量重构方法

技术领域

本发明涉及一种基于自一致性的含联合全变分的并行磁共振成像高质量重构方法。

背景技术

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

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

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

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

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于自一致性的含联合全变分的并行磁共振成像高质量重构方法。

本发明的目的是通过以下技术方案来实现的：基于自一致性的含联合全变分的并行磁共振成像高质量重构方法，其特征在于：它包括以下步骤：

S0：初始化，令u⁰＝D^Ty，z⁰＝0，t¹＝1，j＝1；

式中，x为全部线圈的频域数据，r表示空间位置索引；表示逐线圈的傅里叶变换，D和D_c分别表示选择原频率采样点和未采样点，D^T表示选择原采样频率点并将元采样频率点放回频率域的原位置，表示选择未被采样的频率点并将未被采样的频率点放回频率域的原位置，表示未被采样的频率点，y表示采集到的频率点，z表示中间变量，t表示与算法加速相关的中间变量，j表示循环变量；

S1：计算未采样的重构数据x_g，计算公式如下：

式中，x_g表示未采样的重构数据，b＝-(G-I)D^Ty，G为频域内插操作子；L为的梯度的Lipschitz常数；

S2：将未采样的重构数据x_g转换成步骤S3需要的含噪图像数据v，计算公式如下：

式中，

S3：计算出去噪后的图像数据u^j，包括以下子步骤：

S30：初始化，令

式中，表示多线圈图像变量，N＝m×n，C为多线圈个数，m和n分别为单个线圈二维图像的行数和列数；d＝D_hvx₁，z₁＝Ψx₁，d和z₁为任意变量；b_d和为对偶变量；其中D_n和D_m分别表示n×n和m×m的循环矩阵，所述的循环矩阵结构如下：

S31：令k的值为0，k表示循环变量；

S32：计算：

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

式中，c表示循环索引；

S32：计算：

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

式中，Ψ表示逐线圈的小波变换；

S33：令循环变量c(只在步骤S33～S35中应用)的值为1；

S34：计算：

式中，α₁和α₂为进行rescale处理的正则化参数；

S35：判断循环变量c的值是否大于C：若不大于则对循环变量c做加一操作之后返回步骤S34，否则进入步骤S36；

S36：计算：

S37：计算：

S38：判断k值的大小是否大于K：若不大于K则对k做加一操作之后返回步骤S32，否则进入步骤S39；K为循环次数；

S4：将步骤S34中得到的所有线圈的去噪后的图像数据的集合u^j转换成步骤S5的迭代算法需要的未采样重构数据z^j，计算公式如下：

S5：采用来自于FISTA的迭代算法进行差值加速，包括以下子步骤：

S51：更新t^j+1，计算公式如下：

S52：更新计算公式如下：

S6：判断是否满足条件，若满足条件则进入步骤S7，否则返回步骤S1；所述的条件为迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某值；

S7：计算最终重构的多线圈频域数据，计算公式如下：

S8：对步骤S7得到的各个线圈的频域数据进行傅里叶反变换，采用SRSOS对各个线圈图像进行联合，得到最终的单幅重构图像，计算公式如下：

本发明的有益效果是：本发明基于SPIRiT框架，针对含有JTV和JL1复合正则项的并行成像的重构问题，提出了一种高质量的重构算法。首先将有约束的重构问题，转化为无约束的最优化问题，再对数据保真项和自校正项进行简化，再采用算子分裂技术将简化后的重构问题转化成一个梯度计算问题和一个含有JTV和JL1复合正则项的去噪问题，复合正则项的去噪问题通过全新设计的基于Split Bregman技术的算法求解。最后再通过FISTA进行加速。本发明设计实验，比较了新算法与其他常用算法的重构性能。实验仿真表明，新算法的收敛速度与POCS算法相当，而重构图像的SNR有较大提升。

附图说明

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

图2为本发明步骤S3流程图；

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

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

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

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

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

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

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

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

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

具体实施方式

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

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

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

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

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

x＝Gx； (2)

式中，x为全部线圈的频域数据。G为频域内插操作子，使用整个频域数据与g_ij相卷积来完成内插。

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

y＝Dx； (3)

式中，D为选择操作子，从整个频域中选出采样的频率点，y为采集到的频率点。

由于不同线圈的图像在小波域具有相似性，因此引入联合稀疏性模型(GroupSparse)：

式中，c为线圈索引，r为空间位置索引。

根据公式(2)和公式(3)的两个限制，结合联合稀疏性模型公式(4)，并行成像的重构问题可转化成求解带联合稀疏性正则项(JL1)的优化问题：

式中，其中Ψ为逐线圈的小波变换，为逐线圈的傅里叶变换。Lustig等采用Projection Over Convex Sets(POCS)来求解公式(5)。

另外，Murphy和Vasanawala使用图像域校正，将重构问题写成：

式中，G为图像域内插操作子。Murphy和Vasanawala都是用POCS从图像域直接重构各个线圈图像。频域和图像域校正的POCS算法在本质上类似，因此具有非常相似的重构性能。

在考虑到多线圈图像梯度的相似性，我们在重构时引入Joint Total Variation(JTV)正则项。引入JTV之后，基于频域校正的并行成像的重构问题可转化为以下的优化问题：

式中，Ψ表示逐线圈的小波变换，表示逐线圈的傅里叶变换，G为频域内插操作子。

利用二次罚函数技术，将公式(5)变换为：

在很多情况下，保持已采样的频率系数不变是很有必要的。为此，我们提出了一种简化公式(8)的方法。

假设表示未被采样的频率点，D和D_c分别表示选择原频率采样点和未采样点，D^T表示选择原采样频率点并将其放回频率域的原位置，D_c ^T表示选择未被采样的频率点并将其放回频率域的原位置，则有x可以表示成因此，式(8)可以简化成：

式中，b＝-(G-I)D^Ty，正则化参数α₁和α₂已利用λ₂进行rescale处理。

引入算子分裂算法，公式(9)转化为：

式中，L为的梯度的Lipschitz常数。

根据x_g、的定义可知(不能写成因为和的维度不一样)：

令则将公式(11)改写为：

令则公式(13)可以改写为：

公式(14)描述的是一个含JTV和JL1复合正则项的去噪问题，下面进行求解。

含JTV和JL1复合正则项的去噪问题可以表示成(相当于对公式14的求解)：

式中，y为含噪声图像，为多线圈图像变量，每一列为一个按列堆栈的线圈图像，N＝m×n，C为多线圈个数，m和n分别为单个线圈二维图像的行数和列数。

||x||_JTV按下式进行定义：

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

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

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

于是，原Joint TV去噪问题，即公式(15)可写成：

为了求解公式(19)，引入任意变量d＝D_hvx，z＝Ψx，其中d有以下形式：

利用Split Bregman技术可得：

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

式中，shrink2J(d,T)表示联合二维收缩算子(联合二维软阈值法)，定义为：

其中，β₁表示惩罚参数。

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

式中，shrinkJ(z,T)表示联合一维收缩算子，计算公式如下：

式中，β₂表示惩罚参数。

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

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

由于Ψ^TΨ＝Ι，可被FFT对角化，因此公式(31)可通过FFT快速求解：

于是含JTV和JL1复合正则项的去噪问题可以通过步骤S3进行求解，流程如图2所示。具体地，C为线圈总数，c为相应的循环变量，K为总迭代次数，k为相应的循环变量(即第k次迭代)。

既然所有子问题都能够高效求解，针对公式(9)的求解，本发明提出一种高效的重构方法——基于自一致性的含联合全变分的并行磁共振成像快速算子分裂算法(FastOperator

Splitting Algorithm for SPIRiT-based Parallel Imaging Reconstructionwith Joint Total Variation Regularization，FOSJTV)，FOSJTV算法的具体流程如图1所示。

其中，步骤S6为循环控制条件，可以是迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某个值时，停止循环。步骤S3用于求解公式(14)所示的复合正则项的去噪问题，将含噪图像v进行去噪得到u^j。步骤S2和步骤S4为一对正反操作，步骤S2将未采样的重构数据x_g转换成步骤S3所需的含噪图像数据v，而步骤S4则将去噪后图像数据u^j转换成步骤S5所需的未采样重构数据z^j(原理如公式(12)-(14)之间的叙述)。步骤S5是插值加速方案，来自FISTA，拥有非常优秀的收敛性能。步骤S7为退出循环之后，计算最终重构的多线圈数据(频域)。

由于每次循环重构出的是各个线圈的频域数据，因此还需要对各个线圈频域数据分别进行傅里叶反变换，再使用SRSOS(Square Root Sum Of Squares)方法对各个线圈图像进行联合，才可得到最终的单幅重构图像。这个后处理过程可用SOSF(x)表示，计算方法如公式(33)所示。

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

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

第一组为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。

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

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

式中，MSE表示参考图像与重构图像间的均方误差，Var表示参考图像的方差。为了测试新算法的性能，我们将本发明算法(FOSJTV)与文献中最常用的算法(POCS)、基于非线性共轭梯度的算法(NLCG0)，以及基于简化非线性共轭梯度的算法(NLCG)进行比较，正则项均采用JTV和JL1复合正则项。实验中SPIRiT核大小均设置成5×5。对各个算法的正则化参数进行调优，以使算法达到最优性能。

图6和图7给出了亚采样率为6时新算法FOSJTV与POCS的性能比较。当某种算法重构图像的SNR变化趋于稳定时，就认为该算法已经收敛。图6和图7分别为采用数据序列data1和data2时的实验结果。

从图6可以看出，新算法FOSJTV与POCS算法收敛速度相当，而FOSJTV收敛时重构图像的SNR却大大高于POCS。既然POCS算法拥有非常简单的结构，极易实现，那么以上结论是可以预见的。而对于图7，也可以得到类似的结论。

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

Claims

1.基于自一致性的含联合全变分的并行磁共振成像高质量重构方法，其特征在于：它包括以下步骤：

S0：初始化，令u⁰＝D^Ty，z⁰＝0，t¹＝1，j＝1；

式中，x为全部线圈的频域数据，r表示空间位置索引；表示逐线圈的傅里叶变换，D和D_c分别表示选择原频率采样点和未采样点，D^T表示选择原采样频率点并将原采样频率点放回频率域的原位置，表示选择未被采样的频率点并将未被采样的频率点放回频率域的原位置，表示未被采样的频率点，y表示采集到的频率点，z表示中间变量，t表示与算法加速相关的中间变量，j表示循环变量；

S1：计算未采样的重构数据x_g，计算公式如下：

式中，

S3：计算出去噪后的图像数据u^j，包括以下子步骤：

S30：初始化，令d⁰＝0，z₁ ⁰＝0；

式中，x₁∈C^N×C表示多线圈图像变量，N＝m×n，C为多线圈个数，m和n分别为单个线圈二维图像的行数和列数；d＝D_hvx₁，z₁＝Ψx₁，d和z₁为辅助变量；b_d和为对偶变量；其中D_n和D_m分别表示n×n和m×m的循环矩阵，所述的循环矩阵结构如下：

S31：令k的值为0，k表示循环变量；

S32：计算：

<mrow> <mi>s</mi> <mi>h</mi> <mi>r</mi> <mi>i</mi> <mi>n</mi> <mi>k</mi> <mn>2</mn> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>D</mi> <mrow> <mi>h</mi> <mi>v</mi> </mrow> </msub> <msup> <msub> <mi>x</mi> <mn>1</mn> </msub> <mi>k</mi> </msup> <mo>+</mo> <msubsup> <mi>b</mi> <mi>d</mi> <mi>k</mi> </msubsup> <mo>,</mo> <mn>1</mn> <mo>/</mo> <msub> <mi>&beta;</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mrow> <mo>(</mo> <mo>|</mo> <msub> <mi>D</mi> <mrow> <mi>h</mi> <mi>v</mi> </mrow> </msub> <msup> <msub> <mi>x</mi> <mn>1</mn> </msub> <mi>k</mi> </msup> <mo>+</mo> <msubsup> <mi>b</mi> <mi>d</mi> <mi>k</mi> </msubsup> <mo>|</mo> <mo>-</mo> <mn>1</mn> <mo>/</mo> <msub> <mi>&beta;</mi> <mn>1</mn> </msub> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mfrac> <mrow> <msub> <mi>D</mi> <mrow> <mi>h</mi> <mi>v</mi> </mrow> </msub> <msup> <msub> <mi>x</mi> <mn>1</mn> </msub> <mi>k</mi> </msup> <mo>+</mo> <msubsup> <mi>b</mi> <mi>d</mi> <mi>k</mi> </msubsup> </mrow> <msub> <mi>s</mi> <mn>2</mn> </msub> </mfrac> <mo>;</mo> </mrow> 1

式中，c表示循环索引,

S32：计算：

<mrow> <mi>s</mi> <mi>h</mi> <mi>r</mi> <mi>i</mi> <mi>n</mi> <mi>k</mi> <mi>J</mi> <mrow> <mo>(</mo> <msubsup> <mi>&Psi;x</mi> <mn>1</mn> <mi>k</mi> </msubsup> <mo>+</mo> <msubsup> <mi>b</mi> <msub> <mi>z</mi> <mn>1</mn> </msub> <mi>k</mi> </msubsup> <mo>,</mo> <mn>1</mn> <mo>/</mo> <msub> <mi>&beta;</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mrow> <mo>(</mo> <mo>|</mo> <msubsup> <mi>&Psi;x</mi> <mn>1</mn> <mi>k</mi> </msubsup> <mo>+</mo> <msubsup> <mi>b</mi> <msub> <mi>z</mi> <mn>1</mn> </msub> <mi>k</mi> </msubsup> <mo>|</mo> <mo>-</mo> <mn>1</mn> <mo>/</mo> <msub> <mi>&beta;</mi> <mn>2</mn> </msub> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mfrac> <mrow> <msubsup> <mi>&Psi;x</mi> <mn>1</mn> <mi>k</mi> </msubsup> <mo>+</mo> <msubsup> <mi>b</mi> <msub> <mi>z</mi> <mn>1</mn> </msub> <mi>k</mi> </msubsup> </mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> </mfrac> <mo>;</mo> </mrow>

式中，Ψ表示逐线圈的小波变换；

S33：令循环变量c的值为1；

S34：计算：

式中，α₁和α₂为进行rescale处理的正则化参数；

S35：判断c的值是否大于C：若不大于则对c做加一操作之后返回步骤S34，否则进入步骤S36；

S36：计算：

S37：计算：

S51：更新t^j+1，计算公式如下：

S52：更新计算公式如下：

S7：计算最终重构的多线圈频域数据，计算公式如下：