基于高阶全变分正则化的快速迭代磁共振图像重建方法与流程

本发明涉及磁共振成像技术领域，具体涉及一种压缩感知理论下的基于高阶全变分正则化模型的增广拉格朗日快速迭代磁共振图像重建方法。

背景技术：

磁共振成像因其非创伤性、无电离辐射等优势，在医学检测中发挥着重要作用。然而，过长的扫描时间等问题阻碍着磁共振成像的发展。压缩感知理论利用图像的稀疏性，在采集部分k空间样本的情况下，通过非线性算法重建磁共振图像，实现加速磁共振成像这一目标。

在基于压缩感知理论的磁共振图像重建算法中，常用方法是正则化方法，即将图像重建问题归纳为一个求解代价方程的问题，代价方程中包含数据一致性项和对所求解问题进行约束的正则项；再通过求解此代价方程获得重建图像。

目前使用广泛的一种正则化方法为全变分方法，这一方法因其能在平滑图像的同时有效地保留图像的纹理和边缘等细节信息而被广泛应用。然而，这一方法仍然存在着许多问题和局限性，例如全变分正则化比较适合具有分块常量特性的图像处理问题，而在重建具有较多细节的磁共振图像时常会产生片状伪影等。这主要是由于全变分方法利用图像的一阶导数，因此对于平滑区域、细节区域等重建效果不佳。为了改善这一状况，研究人员提出了许多扩展方法，相比于全变分方法来说，这些方法虽然重建效果有所提高，但大多以增加计算复杂度为代价，降低了计算效率。因此，设计一种计算效率高、耗时短、且能克服全变分产生的伪影问题的磁共振图像重建方法是目前该领域中的一个重要目标。

技术实现要素：

本发明的目的是针对目前磁共振图像重建方法的不足，提供一种可以同时提高重建图像质量及计算效率的基于高阶全变分正则化的快速迭代磁共振图像重建方法。

本发明为解决上述技术问题采取的技术方案是：

(1)利用预先设置的欠采样模板获取磁共振图像的部分k空间数据；

(2)利用获取的部分k空间数据和高阶全变分正则化法建立磁共振图像重建模型；

(3)对部分k空间数据直接进行傅里叶逆变换，得到空间域预估磁共振图像作为初始重建图像；

(4)引入辅助变量和增广拉格朗日乘子以便重建模型的快速迭代求解；

(5)计算初始重建图像中每个像素的高阶方向导数，利用软阈值收缩方法更新高阶方向导数，然后用更新的高阶方向导数获得本次迭代得到的磁共振重建图像；

(6)判断当前重建图像结果是否满足收敛条件，若是则获得最终重建磁共振图像，否则进入步骤(7)；

(7)增加迭代参数取值，以当前迭代步骤中更新的磁共振图像为初始重建图像，返回步骤(5)继续进行循环迭代操作。

上述步骤(2)中重建模型的建立，具体采用如下方式进行：

若将D_u，n定义为n阶微分算子，f为目标重建图像，u∈S表示单位圆S上任一方向的向量，则图像在某一方向上的n阶导数可表示为如下向量形式：

其中，D_u，nf表示图像在单位向量u＝(u_x，u_y)方向上的n阶方向导数，s_n(u)为n阶导数下的单位向量分量，为图像的n阶偏导。实验发现，方向导数阶数过高会导致过度平滑现象，因此本发明所使用的主要为二阶方向导数，即：

其中，

因此，可以根据n阶方向导数，建立高阶全变分正则项即基于高阶方向导数的图像L₁-L₁范数：

利用高阶全变分正则项和欠采样获取的k空间数据，建立基于高阶全变分正则化的压缩感知图像重建模型：

其中A为磁共振图像傅里叶变换后k空间欠采样操作算子，b为获得的k空间欠采样数据，f为待重建图像，D_u，nf为图像的每一像素点在单位圆上任意方向的n阶方向导数，表示L₁范数，λ为用于平衡数据一致性项||Af-b||²和正则项的正则化参数，C(f)为需要求解的代价方程。式(5)的目的即为通过最小化代价方程C(f)，获得重建出的磁共振结果图像

上述步骤(4)具体包括：

由于图像重建模型优化问题(5)中D_u，nf的L₁范数不是处处可微，为了便于求解，将L₁范数用Huber方程进行估计，Huber方程的形式如下：

利用Huber方程估计后的代价方程形式变为：

上式中，β为迭代参数，在每次迭代中，逐渐增加β的取值，当β→∞时，C_β(f)＝C(f)。

上述步骤(4)具体包括：

对于式(7)的求解，我们引入辅助变量z_u和可进一步加快求解过程收敛速度的拉格朗日乘子Λ_u，得到的C_β(f)的逼近形式如下：

至此，原始图像重建优化问题(5)便转换成为一个较易求解的优化问题，即：

上述步骤(5)具体包括：

分别固定方程中的z_u和D_u，nf，可利用交替最小化求解方程，获得重建图像。

首先，假定D_u，nf为固定常数，以z_u为未知数求解(9)方程，具体可表示为：

利用软阈值方法求解上式可得：

其次，假定z_u为固定常数，以D_u，nf为未知数求解(9)方程，忽略式中与z_u相关项可得如下方程：

最小化上式可通过直接求一阶导数来得到结果，因此我们可以得出关于f的二次方程(12)最小时满足下式：

利用导数算子的旋转可操控性质，可以得到和的简单表示形式：

其中，q_n＝∫_Ss_n(u)z_udu，r_n＝∫_Ss_n(u)Λ_udu，E_n为n阶偏导算子，即本发明中的二阶偏导算子即为E₂f＝[f_xx，f_xy，f_yy]^T；

在求解z_u的过程中，为了能够获得较好的图像重建结果，需在所有的方向上对z_u进行考虑，这使得当z_u在某一方向上数量级比较大时所需的计算存储容量会变得很大。但在f的求解中可发现，实际只需要计算z_u在s(u)上的投影q_n(u)＝∫_Ss(u)z_udu，这样的简化便使得计算所需的存储代价减少了很多，因此可以大幅提高计算效率。

根据上述推导，式(13)等价于：

在离散傅里叶域，二阶算子能够简单地表达为如下形式：

其中，为离散傅里叶变换算子。当对磁共振图像的欠采样在笛卡尔坐标系下进行时，可将A算子表示为其中S为频域采样算子，表示对图像的傅里叶变换和欠采样两个步骤。因此，式(17)可以写成：

由此，可以得到f的解析表达式，即本次迭代中的重建图像结果：

当采样不在非笛卡尔坐标系下进行，也可简单地利用共轭梯度法进行最小化计算，求解(17)。

在每一次迭代计算后需要将引入的拉格朗日乘子Λ_u进行如下更新：

其中为本次迭代所使用的拉格朗日乘子，为用于下一次迭代中的拉格朗日乘子，D_u，nf和z_u分别为当前迭代所使用初始图像的高阶方向导数和z_u的解。

上述步骤(6)具体包括：

判断当前重建图像是否满足收敛条件的方法为计算相邻两次重建图像的L₂范数误差，即当误差小于预定阈值时，停止迭代，获得最终重建磁共振图像，否则继续循环迭代。

上述步骤(7)具体包括：

每次迭代中增加迭代参数β的取值：

β＝β₀·β_inc (22)

其中β₀为本次迭代中所使用的参数值，β_inc为每次循环所增加的倍数。

本发明的有益效果是：

本发明所提出的算法计算效率较高，CPU运行时间较短，是一种高效的算法；同时利用本发明提出的算法，可以较好地重建磁共振图像，较完整地保留图像中的细节区域。实验表明，与全变分方法、图像高阶导数Laplacian方法、小波方法等方法相比，本发明能够在较高加速系数下获得较高质量的重建图像，同时提高了重建速度，从而实现了加速磁共振图像重建的目的

附图说明

图1为本发明方法流程图(高阶全变分正则化磁共振图像重建算法流程图)。

图2为仿真实验所用参考磁共振图像；图中：(a)为侧向脑部图像，(b)为轴向脑部图像，(c)为血管造影图像，(d)为腕部图像。

图3为利用不同方法对图2-(a)侧向脑部磁共振图像在不同欠采样参数和噪声水平下重建图像对比(磁共振图像侧向脑部图像重建结果)；图中：

(a)为侧向脑部图像放大图(原磁共振图像)，(b)为三阶变分重建结果图(4.35倍欠采样，添加40dB噪声；信噪比SNR＝21.01dB)，(c)为曲波变换重建结果图(4.35倍欠采样，添加40dB噪声；SNR＝17.25dB)，(d)为拉普拉斯重建结果图(4.35倍欠采样，添加40dB噪声；SNR＝16.19dB)，

(e)为全变分重建结果图(SNR＝23.55dB；4.35倍欠采样，添加40dB噪声)，(f)为二阶变分重建结果图(SNR＝24.35dB；4.35倍欠采样，添加40dB噪声)，(g)为各项同性二阶变分图(SNR＝16.19dB；4.35倍欠采样，添加40dB噪声)，(h)为Lysaker重建结果图(SNR＝16.19dB；4.35倍欠采样，添加40dB噪声)，

(i)为全变分重建结果图(SNR＝28.40dB；2倍欠采样，添加40dB噪声)，(j)为二阶变分重建结果图(SNR＝30.15dB；2倍欠采样，添加40dB噪声)，(k)为各项同性二阶变分图(SNR＝29.90dB；2倍欠采样，添加40dB噪声)，(l)为Lysaker重建结果图(SNR＝26.89dB；2倍欠采样，添加40dB噪声)。

图4为利用不同方法对图2-(d)腕部磁共振图像在不同欠采样参数和噪声水平下重建图像对比(欠采样和添加噪声下的腕部图像重建结果图)；图中；

(a)为原图像，(b)为全变分正则化重建结果图(3倍欠采样，添加噪声35dB，重建结果SNR＝25.34dB)，(c)全变分正则化误差图(3倍欠采样，添加噪声35dB，重建结SNR＝25.34dB)，

(d)为原图像放大图，(e)为二阶变分正则化结果图(3倍欠采样，添加噪声35dB，重建结果SNR＝25.65dB)，(f)为二阶变分正则化误差图(3倍欠采样，添加噪声35dB，重建结果SNR＝25.65dB)。

图5为利用不同方法对图2-(c)血管造影磁共振图像在不同欠采样参数和噪声水平下重建图像对比(欠采样和添加噪声下的血管造影图像重建结果)；图中：

(a)为原图像，(b)为全变分方法重建结果图(2.86倍欠采样，添加噪声30dB，重建结果SNR＝29.24dB)，(c)为全变分方法误差图(2.86倍欠采样，添加噪声30dB，重建结果SNR＝29.24dB)，

(d)为原图像放大图，(e)为二阶全变分方法重建结果图(2.86倍欠采样，添加噪声30dB，重建结果SNR＝31.30dB)，(f)为二阶全变分方法误差图(2.86倍欠采样，添加噪声30dB，重建结果SNR＝31.30dB)。

图6为本发明方法的收敛性能及运算效率分析对比图；图中：(a)为不同β值所对应的收敛趋势曲线图，(b)为本发明所用算法的CPU计算时间(图中，短的浅色曲线为对尺寸为256×256的图像进行重建；长的深色曲线为对尺寸为512×512的图像进行重建)。

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

如图1所示，本发明的具体实施步骤如下：

(1)利用预先设置的欠采样模板获取部分k空间数据，为了验证本发明的效果，采用了四组参考磁共振图像，如图2所示，分别为侧向脑部磁共振图像(a)，轴向脑部磁共振图像(b)，血管造影磁共振图像(c)，腕部磁共振图像(d)，对参考图像进行傅里叶变换，采集原始k空间数据，采集到的欠采样k空间数据表示为b＝Af+n，其中A为对磁共振图像进行傅里叶变换后在k空间进行欠采样操作算子，n为实际采样中可能存在的加性噪声，b为获得的k空间欠采样数据，f为待重建图像；

对测量数据b直接进行补零傅里叶逆变换得到空间域初始重建图像其中表示对测量数据做傅里叶逆变换；

(2)利用获取的部分k空间数据和高阶全变分正则化法建立磁共振图像重建模型；

(3)对部分k空间数据直接进行傅里叶逆变换，得到空间域预估磁共振图像作为初始重建图像；

(4)引入辅助变量和增广拉格朗日乘子以便重建模型的快速迭代求解；

(5)计算初始重建图像中每个像素的高阶方向导数，利用软阈值收缩方法更新高阶方向导数，然后用更新的高阶方向导数获得本次迭代得到的磁共振重建图像；

(6)判断当前重建图像结果是否满足收敛条件，若是则获得最终重建磁共振图像，否则进入步骤(7)；

(7)增加迭代参数取值，以当前迭代步骤中更新的磁共振图像为初始重建图像，返回步骤(5)继续进行循环迭代操作。

上述步骤(2)中重建模型的建立，具体采用如下方式进行：

将图像在某一方向上的n阶方向导数表示为向量形式：

其中f为目标图像，D_u，n为n阶微分算子，D_u，nf表示图像在单位向量u＝(u_x，u_y)方向上的n阶方向导数，s_n(u)为n阶导数下的单位向量分量，为n阶图像偏导。实验发现方向导数阶数过高会导致过度平滑现象，因此本发明所使用的主要为二阶方向导数，即：

其中，

因此，可以根据n阶方向导数，建立高阶全变分正则项即基于高阶方向导数的图像L₁-L₁范数：

本发明中所使用的是二阶全变分正则项，即：

利用高阶全变分正则项和欠采样获取的k空间数据，建立基于高阶全变分正则化的压缩感知图像重建模型：

其中A为磁共振图像傅里叶变换后k空间欠采样操作算子，b为获得的k空间欠采样数据，f为待重建图像，D_u，nf为图像的每一像素点在单位圆上任意方向的n阶方向导数，λ为用于平衡数据一致性项||Af-b||²和正则项的正则化参数，C(f)为需要求解的代价方程。式(6)的目的即为通过最小化代价方程C(f)，获得重建出的磁共振图像

上述步骤(4)具体包括：

由于图像重建模型优化问题(6)中L₁范数方程不是处处可微，为了便于求解，将L₁范数用Huber方程进行估计，Huber方程的形式如下：

利用Huber方程估计后的代价方程形式变为：

上式中，当β→∞时，C_β(f)＝C(f)。

对于式(8)的求解，我们引入辅助变量z_u和可进一步加快求解过程收敛速度的拉格朗日乘子Λ_u，得到的C_β(f)的逼近形式如下：

至此，原始图像重建优化问题(6)便转换成为一个较易求解的优化问题，即：

上述步骤(5)具体包括：

分别固定方程中的z_u和D_u，nf，可利用交替最小化求解方程，获得重建图像。

首先，假定D_u，nf为固定常数，以z_u为未知数求解(10)方程，具体可表示为：

利用软阈值方法求解上式可得：

其次，假定z_u为固定常数，以D_u，nf为未知数求解(10)方程，忽略式中与z_u相关项可得如下方程：

最小化上式可通过直接求一阶导数来得到结果，因此我们可以得出关于f的二次方程(13)最小时满足下式：

利用导数算子的旋转可操控性质，可以得到和的简单表示形式：

其中，q_n＝∫_Ss_n(u)z_udu，r_n＝∫_Ss_n(u)Λ_udu，E_n为n阶偏导算子，即本发明中的二阶偏导算子即为E₂f＝[f_xx，f_xy，f_yy]^T；

在求解z_u的过程中，为了能够获得较好的图像重建结果，需在所有的方向上对z_u进行考虑，这使得当z_u在某一方向上数量级比较大时所需的计算存储容量会变得很大。但在f的求解中可发现，实际只需要计算z_u在s(u)上的投影q_n(u)＝∫_Ss(u)z_udu，这样的简化便使得计算所需的存储代价减少了很多，因此可以大幅提高计算效率。

根据上述推导，式(14)等价于：

在离散傅里叶域，二阶算子能够简单地表达为如下形式：

其中，为离散傅里叶变换算子。当对磁共振图像的欠采样在笛卡尔坐标系下进行时，可将A算子表示为其中S为频域采样算子，表示对图像的傅里叶变换和欠采样两个步骤。因此，式(18)可以写成：

因此，可以得到f的解析表达式，即本次迭代中的重建图像结果：

当采样不在非笛卡尔坐标系下进行，也可简单地利用共轭梯度法进行最小化计算，求解(18)。

在每一次迭代计算后需要将引入的拉格朗日乘子Λ_u进行如下更新：

其中为本次迭代所使用的拉格朗日乘子，为用于下一次迭代中的拉格朗日乘子，D_u，nf和z_u分别为当前迭代所使用初始图像的高阶方向导数和z_u的解。

上述步骤(6)具体包括：

判断当前重建图像是否满足收敛条件的方法为计算相邻两次重建图像的L₂范数误差，即当误差小于预定阈值时，停止迭代，获得最终重建磁共振图像，否则继续循环迭代。

上述步骤(7)具体包括：

每次迭代中增加迭代参数β的取值：

β＝β₀·β_inc (23)

其中β₀为本次迭代中所使用的参数值，β_inc为每次循环所增加的倍数。

为了定量分析本发明中所进行实验的结果，本发明采用SNR和CPU运行时间两个指标对结果进行分析：

其中f_ref为参考图像，为重建出的磁共振结果图像。

图3为利用图2(a)侧向脑部磁共振图像，在不同欠采样参数下基于不同方法重建出的结果图像对比。其中，(b)-(h)为在4.35倍欠采样参数下，分别使用三阶全变分、曲波变换、拉普拉斯高阶方法、全变分、二阶全变分、各向同性二阶全变分和Lysaker高阶方法这七种方法对其进行重建的效果图。可以看出，三阶全变分由于图像导数阶数较高，因此造成了图像过度模糊，信噪比较低；全变分导致了重建结果出现片状伪影；曲波变换和拉普拉斯高阶方法在重建图像中产生曲线状和圆形伪影，因此重建效果不佳。综合考虑，二阶全变分对图像重建效果最好。(i)-(l)为2倍欠采样参数下分别用全变分、二阶全变分、各向同性二阶全变分和Lysaker高阶方法这四种方法的重建结果图，结合箭头所指位置进行比对，二阶全变分优于其他方法。

图4和图5分别为对图2(d)血管造影磁共振图像和2(c)腕部磁共振图像进行3倍欠采样同时添加35dB噪声，对图2-(c)进行2.86倍欠采样同时添加30dB噪声后全变分和二阶全变分方法重建结果比对。图4(c)、4(f)、5(c)、5(f)是对原图像和重建结果图像做差值得出的误差图，从中可以看出，本文所提出的二阶全变分方法无论在定性方面或者定量方面都优于全变分正则化方法。

在图6(a)中，通过在迭代过程中改变β参数来观察其对收敛速度的影响，两红颜色折线代表的是β初始迭代值为15，在每一次迭代中分别增大2倍和5倍的误差收敛趋势；两蓝颜色折线是初始迭代值为50，在每一次迭代中分别增大2倍和5倍的误差收敛趋势。从此图中可看出，β的初始值以及每次迭代的增量大小并不会对收敛速度产生特别明显的影响。因此，只要满足β值随每次迭代而增加，都能在一定迭代次数内获得较好的重建结果。图6(b)为本发明所使用算法的CPU运行时间。其中浅颜色与深颜色曲线分别为使用本发明方法对尺寸为256×256和512×512的磁共振图像进行重建的信噪比改善趋势，从对应坐标中可看出，本发明所使用的算法CPU运行时间非常短，是一种高效的算法。

表1为部分实验数据，若以结果图的SNR为衡量重建质量的标准，可看出高阶变分方法较其他方法效果更好。表1为利用不同方法对图2中四幅参考磁共振图像在不同欠采样在不同欠采样参数和噪声水平下重建图像的信噪比值列表。

表1在不同采样参数和噪声水平下利用不同方法重建磁共振图像的SNR列表