CN112617798B - 一种基于Lp范数联合全变分的并行磁共振成像重构方法 - Google Patents

Abstract

本发明涉及一种基于Lp范数联合全变分的并行磁共振成像重构方法，属于磁共振成像技术领域。为了提高重构图像的质量，本发明结合了Lp范数联合全变分正则项，提出了一种基于交替方向乘子法求解的局部k空间邻域建模算法。在人体的脑部和膝盖数据集上进行了实验，结果表明提出的算法与并行局部k空间邻域建模算法和结合了全变分的并行局部k空间邻域建模算法相比，提出的算法能更好地保留重构图像的边缘和细节信息，重构效果更好。

Description

一种基于Lp范数联合全变分的并行磁共振成像重构方法

技术领域

本发明涉及一种基于Lp范数联合全变分的并行磁共振成像重构方法，属于磁共振成像技术领域。

背景技术

磁共振成像(Magnetic resonance imaging,MRI)是一种利用核磁共振原理进行的医学成像技术，具有无创伤、无电离辐射、高分辨率和任意方向断层扫描等优点，是当今医学诊断和治疗监测的重要手段之一。尽管MRI技术已经广泛应用到临床医疗中，但由于成像时间较长极大地限制了MRI技术的发展。

压缩感知技术(Compressed Sensing,CS)是近年来新兴的一种图像处理技术，突破了传统奈奎斯特采样频率的限制，在不满足香农定理的情况下，从少量的测量数据中恢复原始信号。因此将CS应用到MRI中，可以大幅度地减少图像的扫描时间，提高成像速度。

目前，常用的MRI加速方法还有并行成像。并行磁共振成像技术以多通道相控阵线圈同时采集磁共振信号，利用各个线圈的空间灵敏度差异来编码空间信息，通过减少相位编码有效提高了MRI成像速度。一类并行成像方法是已知灵敏度信息的，根据测量的线圈灵敏度信息来重构图像。如：SMASH，SENSE和SPACE RIP。然而，这些方法最主要的障碍是很难准确的测量接收线圈的灵敏度，微小的扰动误差会导致成像图像中出现可见的伪影，使重构结果不一致。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于Lp范数联合全变分的并行磁共振成像重构方法，基于PLORAKS模型，结合了Lp范数联合全变分正则项(joint TotalVariation,JTV)，实验仿真结果表明，提出的新算法LpJTV-PLORAKS在保证MRI成像速度的同时，提高了重构图像的质量。

本发明采用的技术方案是：一种基于Lp范数联合全变分的并行磁共振成像重构方法，包括以下步骤：

S0：初始化，令：f¹＝0，V¹＝0，B_V ¹＝0，i＝1；

其中，f＝[f₁,...,f_j,...,f_L]表示重构的多线圈k空间数据，L表示接收线圈的总个数，f¹表示待重构k空间数据的初始值，V＝[V₁,...,V_z,...,V_N]为Dx对应的辅助变量，x＝F^-1f，F^-1为逐线圈二维傅里叶反变换，x＝[x₁,...,x_j,...,x_L]表示重构的多线圈图像，Dx表示多线圈图像x的一阶有限差分变换，z表示k空间位置索引变量，

表示提取Dx的空间位置为z的所有点组成的列向量，式中，

和

分别为行方向和列方向的一阶有限差分变换，I_m和I_n分别为m×m和n×n的单位矩阵，D_m和D_n是在周期性边界条件下的一阶有限差分算子矩阵，分别为m×m和n×n的循环矩阵，其结构为

是辅助变量V对应的对偶变量，i是循环变量，表示数据的第i次迭代，上标1表示第1次迭代；

S1：计算x＝F^-1f，并计算辅助变量V_z ⁱ⁺¹，计算公式如下：

其中，V_z ⁱ⁺¹表示第i+1次迭代的第z个k空间位置的多线圈图像x的一阶有限差分变换(Dxⁱ⁺¹)_z对应的辅助变量，(Dxⁱ)_z表示第i次迭代的k空间位置为z的多线圈图像x的一阶有限差分变换，

表示第i次迭代辅助变量

对应的对偶变量，(||g||)^p-2表示

的l₂范数的p-2次幂且p∈(0,1)，δ是惩罚项的参数；

S2：判断z的值是否大于N；若是，则进入步骤S3；否则，对z进行加一操作后返回S1；

S3：计算待重构的多线圈图像的k空间数据fⁱ⁺¹，使用共轭梯度法求解下述问题得到fⁱ⁺¹：

其中，fⁱ是第i次迭代的多线圈图像的k空间数据，

是N×N的单位阵I_N的M行组成的采样矩阵，

是M×1的欠采样数据，λ为正则化参数，α为JTV正则化参数，δ是惩罚项的参数，Vⁱ表示第i次迭代的多线圈图像x的一阶有限差分变换Dxⁱ对应的辅助变量，Dxⁱ表示第i次迭代的多线圈图像x的一阶有限差分变换，B_V ⁱ表示第i次迭代辅助变量V对应的对偶变量，

是一个构造高维矩阵的线性算子，能够将k空间数据f_j(n_h,n_y)构造成如下矩阵：

f_j表示第j个线圈图像的k空间数据，矩阵R_S(f_j)是近似低秩的，其中，

矩阵

分别表示为：

k＝1,...,K，m＝1,...,N_R，

和

分别代表

的实部和虚部，上标～表示该数据是频域数据，

表示K个不同的k空间位置，(n_h,n_y)满足-N_h+R≤n_h≤N_h-R，-N_y+R≤n_y≤N_y-R，N_h和N_y表示k空间度量区域的正整数，R为傅里叶截断半径，当

时，

使用

表示集合

中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下表m表示第m个像素点，N_R表示Λ_R的基数，J_r(·)是非凸正则化函数，假设

则J_r(X)的表达式为：

S4：计算对偶变量B_V ⁱ⁺¹，计算公式如下：

B_V ⁱ⁺¹＝B_V ⁱ+(DF^-1fⁱ⁺¹-Vⁱ⁺¹)；

其中，B_V ⁱ⁺¹表示第i+1次迭代辅助变量V对应的对偶变量，DF^-1fⁱ⁺¹表示第i+1次迭代的多线圈图像F^-1f的一阶有限差分变换；

S5：判断是否达到停止准则，若满足条件则进入步骤S6，否则对循环变量i做加一操作之后返回步骤S1；停止准则为迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某值时；

S6：输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成重构图像，计算公式如下：

本发明的有益效果是：基于LORAKS模型，本发明提出了一种基于Lp范数联合全变分正则项的并行成像重建算法，首先，使用ADMM算法将磁共振图像重构问题转换成一个梯度计算问题和一个去噪问题，其次分别求解。理论分析和成像实验表明，与传统的重构算法相比，本发明提出的新的并行成像算法可以更有效地提高重构图像的质量。

附图说明

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

图2为8通道的膝盖数据(即data1)图；

图3为加速因子为6的泊松亚采样示意图；

图4-6为数据集data1下，分别使用P-LORAKS，JTV-PLORAKS，LpJTV-PLORAKS重建并行MR图像(具有6倍加速度和24×24中心校准区)；

图7-9为数据集data1下，分别使用P-LORAKS，JTV-PLORAKS，LpJTV-PLORAKS重建的误差图(具有6倍加速度和24×24中心校准区)。

具体实施方式

下面结合附图和具体实施例，进一步详细描述本发明的技术方案。

实施例1：如图1-9所示，本发明是基于PLORAKS框架提出的一种高效的重构方法。

并行磁共振成像技术是一种有效加快磁共振成像(Magnetic resonanceimaging,MRI)速度的方法，它使用多个接收线圈同时采集感应信号，从而有效地减少MRI的扫描时间。局部k空间邻域模型是最近提出的一种k空间低秩约束重构模型，它利用图像的线性位移不变性，将图像的k空间数据映射到高维矩阵中来解决图像重构问题。

在PLORAKS模型中，当线圈图像在k空间中具有线性依赖关系或者缓慢变化的相位时，全采样的k空间数据可以映射到一个低秩矩阵，这就使得利用低秩正则化来重构欠采样数据或者噪声数据成为可能，本发明用下面这个式子表示MRI的采样过程：

Af＝d (1)

其中，f＝[f₁,...,f_j,...,f_L]，

为第j个线圈图像的k空间数据，

是欠采样矩阵，

表示欠采样数据。

S矩阵可通过对k空间数据f_j(n_h,n_y)选取K个不同的像素点进行线性运算的线性算子

来构造，具体形式如下：

R_PS(f)＝[R_S(f₁),...,R_S(f_j),...,R_S(f_L)] (3)

矩阵R_S(f_j)是近似低秩的，

其中，矩阵

分别表示为：

k＝1,...,K，m＝1,...,N_R，

和

分别代表

的实部和虚部。

表示K个不同的k空间位置，(n_h,n_y)满足-N_h+R≤n_h≤N_h-R，-N_y+R≤n_y≤N_y-R，N_h和N_y表示k空间度量区域的正整数，R为傅里叶截断半径，当

时，

使用

表示集合

中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下表m表示第m个像素点，N_R表示Λ_R的基数。

P-LORAKS重建问题可以用下面的公式表示：

式中，

是每个线圈全采样k空间数据f_j的集合；

是每个线圈欠采样数据d_j的集合。

其中σ_k是矩阵X的第k个奇异值，I_Q-r是一个(Q-r)×(Q-r)其中，λ为正则化参数，J_r(·)是非凸正则化函数，它迫使矩阵参数的秩小于等于r，r为一个整数，假设

则：

其中，σ_k是矩阵X的第k个奇异值，I_Q-r是一个(Q-r)×(Q-r)的单位矩阵。T是X的最优低秩近似，G是一个零填充矩阵。

PLORAKS模型为并行成像提供了一种新的重构方法。尽管PLORAKS指出了结合了JTV正则项的PLORAKS模型重构算法可有效改善重构质量，但仍有较大的改进余地。为了进一步提高重构质量，本发明将lp范数联合全变分正则项与PLORAKS模型相结合，提出LpJTV-PLORAKS算法，得到如下最优化问题：

其中，α为JTV正则项的参数，

表示l_P范数联合全变分正则项，

表示

的l₂范数的p次幂，

表示提取Dx的空间位置为z的所有点组成的列向量，表达式为：

式中

和

分别为行方向和列方向的一阶有限差分变换，D_m和D_n是在周期性边界条件下的一阶有限差分算子矩阵，分别为m×m和n×n的循环矩阵，其结构为

由上可得，问题(10)还可以转化为：

考虑隐函数

的加性半二次优化，

可表示为：

式中，

和

是辅助变量，δ是惩罚项的参数，ψ(w)是依赖于

的矩阵函数。式(13)可以利用下式来求解：

使用增广拉格朗日乘子法引入辅助变量

和对偶变量

即

则式(12)可以转化为：

使用ADMM算法，可以将最小化问题(15)分解为如下子问题：

B_V ⁱ⁺¹＝B_V ⁱ+(DF^-1fⁱ⁺¹-Vⁱ⁺¹) (18)

其中，子问题(18)是对偶变量的更新方法。

子问题(17)可以根据式(14)求解得：

子问题(16)可以通过共轭梯度法来求解fⁱ⁺¹，且在周期性边界条件下，D^TD可以通过二维快速傅里叶变换实现对角化。

具体流程如图1所示，其中步骤如下：

S0：初始化，令：f¹＝0，V¹＝0，B_V ¹＝0，i＝1；

其中，f¹表示重构的多线圈k空间数据的初始值，V¹表示多线圈图像x的一阶有限差分变换对应的辅助变量的初始值，B_V ¹表示V¹对应的对偶变量的初始值；

S1：计算辅助变量V_z ⁱ⁺¹，计算公式如下：

S2：判断z值是否大于N；若是，则进入步骤S3；否则，对z进行加一操作后返回S1；

S3：计算待重构的多线圈图像的k空间数据fⁱ⁺¹，使用共轭梯度法求解下述问题得到fⁱ⁺¹：

S4：计算对偶变量B_V ⁱ⁺¹，计算公式如下：

B_V ⁱ⁺¹＝B_V ⁱ+(DF^-1fⁱ⁺¹-Vⁱ⁺¹)；

S5：判断是否达到停止准则，若满足条件则进入步骤S6，否则对循环变量i做加一操作之后返回步骤S1；停止准则为迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某值时；

S6：输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成重构图像，计算公式如下：

下面结合具体实验数据，说明本发明的效果。

为了测试提出的LpJTV-PLORAKS算法的重构性能，本发明分别与P-LORAKS算法和JTV-PLORAKS算法进行了比较。所有算法均用MATLAB实现。

开展仿真实验所用的数据为一个人体膝盖数据集P1_imk_100，命名为data1(如图2)，尺寸为320×320。为了生成测试数据集，使用加速因子为6的笛卡尔泊松亚采样，对全采样的数据集进行人工采样并用于仿真重建，如图3所示。在以下实验中，对于所有比较的算法，使用尺寸为24×24校准区域。

所有实验均在配置为Intel core i7-8700@3.20GHz的GPU，16GB的内存和Ubuntu18.04.5操作系统的台式机上进行。所有实验参数均根据信噪比最佳进行调优。本发明使用信噪比(Signal Noise Ratio,SNR)定量评估重建图像的质量。SNR的数值越高表明图像重构质量越好。SNR的计算方法如下：

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

首先，本发明对数据集data1，在加速因子为6时各个算法的重构图像进行了视觉上的比较，图4-6展示了在数据集data1上，算法P-LORAKS、算法JTV-PLORAKS和算法LpJTV-PLORAKS的重构图像，在图7-9中展示了三种算法在数据集data1下的误差图像(误差图颜色越深，表示误差越小)。从图4-6可看出，P-LORAKS算法的重构图像在纹理细节方面比较模糊，整个重构图像存在轻微的伪影，且误差点较为分散；JTV-PLORAKS算法与PLORAKS算法相比，重构图像的中间区域比较清晰，图像细节得以保留，整个图像的误差点明显减少；LpJTV-PLORAKS算法与以上两种算法相比，更好地保留了中间区域的纹理细节和边缘的轮廓信息，且误差图的颜色也是明显变深，表明了重构图像与参考图像的一致性较高。

为了更加客观地评价本发明所提算法的重构性能，本发明还使用SNR定量地评估在加速因子为6时，以上三种算法在数据集data1上的重构性能。LpJTV-PLORAKS算法SNR值比TV-PLORAKS算法和PLORAKS算法分别提高了1.34dB和3.11dB。

以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (1)

1.一种基于Lp范数联合全变分的并行磁共振成像重构方法，其特征在于：包括以下步骤：

S0：初始化，令：f¹＝0，V¹＝0，B_V ¹＝0，i＝1；

其中，f＝[f₁,...,f_j,...,f_L]表示重构的多线圈k空间数据，L表示接收线圈的总个数，f¹表示待重构k空间数据的初始值，V＝[V₁,...,V_z,...,V_N]为Dx对应的辅助变量，x＝F^-1f，F^-1为逐线圈二维傅里叶反变换，x＝[x₁,...,x_j,...,x_L]表示重构的多线圈图像，Dx表示多线圈图像x的一阶有限差分变换，z表示k空间位置索引变量，

表示提取Dx的空间位置为z的所有点组成的列向量，式中，

和

分别为行方向和列方向的一阶有限差分变换，I_m和I_n分别为m×m和n×n的单位矩阵，D_m和D_n是在周期性边界条件下的一阶有限差分算子矩阵，分别为m×m和n×n的循环矩阵，其结构为

是辅助变量V对应的对偶变量，i是循环变量，表示数据的第i次迭代，上标1表示第1次迭代；

S1：计算x＝F^-1f，并计算辅助变量V_z ⁱ⁺¹，计算公式如下：

其中，V_z ⁱ⁺¹表示第i+1次迭代的第z个k空间位置的多线圈图像x的一阶有限差分变换(Dxⁱ⁺¹)_z对应的辅助变量，(Dxⁱ)_z表示第i次迭代的k空间位置为z的多线圈图像x的一阶有限差分变换，

表示第i次迭代辅助变量

对应的对偶变量，(||g||)^p-2表示

的l₂范数的p-2次幂且p∈(0,1)，δ是惩罚项的参数；

S2：判断z的值是否大于N；若是，则进入步骤S3；否则，对z进行加一操作后返回S1；

S3：计算待重构的多线圈图像的k空间数据fⁱ⁺¹，使用共轭梯度法求解下述问题得到f^{i +1}：

其中，fⁱ是第i次迭代的多线圈图像的k空间数据，

是N×N的单位阵I_N的M行组成的采样矩阵，

是M×1的欠采样数据，λ为正则化参数，α为JTV正则化参数，δ是惩罚项的参数，Vⁱ表示第i次迭代的多线圈图像x的一阶有限差分变换Dxⁱ对应的辅助变量，Dxⁱ表示第i次迭代的多线圈图像x的一阶有限差分变换，B_V ⁱ表示第i次迭代辅助变量V对应的对偶变量，R_S(·):

是一个构造高维矩阵的线性算子，能够将k空间数据f_j(n_h,n_y)构造成如下矩阵：

f_j表示第j个线圈图像的k空间数据，矩阵R_S(f_j)是近似低秩的，其中，

矩阵

分别表示为：

k＝1,...,K，m＝1,...,N_R，

和

分别代表

的实部和虚部，上标～表示该数据是频域数据，

表示K个不同的k空间位置，(n_h,n_y)满足-N_h+R≤n_h≤N_h-R，-N_y+R≤n_y≤N_y-R，N_h和N_y表示k空间度量区域的正整数，R为傅里叶截断半径，当

时，

使用

表示集合

中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下标m表示第m个像素点，N_R表示Λ_R的基数，J_r(·)是非凸正则化函数，假设

则J_r(X)的表达式为：

σ_k是矩阵X的第k个奇异值，T是X的最优低秩近似；

S4：计算对偶变量B_V ⁱ⁺¹，计算公式如下：

B_V ⁱ⁺¹＝B_V ⁱ+(DF^-1fⁱ⁺¹-Vⁱ⁺¹)；

其中，B_V ⁱ⁺¹表示第i+1次迭代辅助变量V对应的对偶变量，DF^-1fⁱ⁺¹表示第i+1次迭代的多线圈图像F^-1f的一阶有限差分变换；

S5：判断是否达到停止准则，若满足条件则进入步骤S6，否则对循环变量i做加一操作之后返回步骤S1；停止准则为迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某值时；

S6：输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成重构图像，计算公式如下：

Images