CN109920017B - 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法 - Google Patents

Abstract

本发明涉及基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法，属于医学磁共振成像技术领域。本发明基于特征向量的迭代自一致性的并行成像重构ESPIRiT框架，提出了一种含联合总变分Lp伪范数正则项的ESPIRiT并行成像重构算法。首先，采用OS技术将重构问题分解为一个梯度计算问题和一个含稀疏正则项的去噪问题。其次，应用MM算法进行去噪。最后，使用FISTA进行加速。实验结果表明，提出的基于ESPIRiT重构模型使用LpJTV正则项的新算法与使用L1正则项的传统重构算法相比，提出的算法能够更有效地改善图像的重构质量。

Description

基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共 振成像重构方法

技术领域

本发明涉及基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法，属于磁共振成像技术领域。

背景技术

磁共振成像(Magnetic Resonance Imaging,MRI)是一种非侵入性和非电离方式的人体结构可视化的成像技术，与其他大多数成像方式相比，它在不同的软组织中提供了良好的对比度。然而，受到数据采样时间的限制，磁共振成像速度通常较慢。并行成像技术的提出，能够有效地减少扫描时间，加速磁共振成像。目前，主要有两种并行成像方法：(1)已知灵敏度信息的，每一个线圈均能准确估计灵敏度轮廓，通过估计的灵敏度信息重构图像。如：灵敏度编码(SENSitivity Encoding,SENSE)；然而，这一方法最主要的障碍是很难准确的测量接收线圈的灵敏度，微小的扰动误差会导致成像图像中出现可见的伪影，使重构结果不一致。(2)基于k空间局部内核的，灵敏度信息是来自于多通道的k空间邻域点的自校准信号。如：一般自校准部分并行采样(GeneRalized Autocalibrating PartiallyParallel Acquisitions,GRAPPA)和迭代自一致性的平行成像重构(iTerative Self-consistent Parallel Imaging Reconstruction,SPIRiT)。这一方法最主要的限制是要获取足够的自校准信号通常是不可行的，且对于高加速因子这一方法是无效的。

由于上述两种成像技术的不足，Uecker等提出了利用特征向量算子的迭代自一致性的并行成像重构(Iterative self-consistent parallel imaging reconstructionusing eigenvector maps,ESPIRiT)模型，将SENSE和GRAPPA限定为子空间来重构缺失的数据(SENSE通过使用预先计算的线圈灵敏度maps来实现线圈图像的重构，GRAPPA利用k空间的校准内核的滤波)。ESPIRiT通过间接的方式实现线圈灵敏度maps的准确估计。且应用k空间算子的主要特征向量来表示线圈灵敏度maps，这些灵敏度算子可以通过在图像空间中的使用奇异值分解来快速计算，从而能够从k空间中心的自校准区域中有效的估计线圈灵敏度。

另外，ESPIRiT还结合了L₁正则项约束模型，能有效的提高成像图像的质量。但是，由于动态腐蚀、化学位移和小视场(Field of view，FOV)等因素的影响，使得含L₁正则项的ESPIRiT重构方法的重构图像中出现一些少量的重叠伪影。因此，为了进一步改善图像的重构质量，本发明基于ESPIRiT重构模型，提出了含Lp伪范数联合总变差(Lp pseudo-normJoint Total Variation，LpJTV)正则项的ESPIRiT并行成像重构算法。

发明内容

本发明的目的在于克服现有技术的不足，提供基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法。

本发明采用的技术方案如下：基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法，包括以下步骤：

Step0：初始化，令w¹＝0，x¹＝0，t¹＝0，i＝1；

其中，w为含噪的图像数据，w¹表示含噪的图像数据的初始值；x表示重构的图像，x¹是重构图像的初始值；t为时间加速因子，t¹表示时间加速因子的初始值；i表示循环变量；

Step1：计算中间变量zⁱ⁺¹，计算公式如下：

式中，zⁱ⁺¹为第i+1次循环的中间变量，wⁱ为第i次循环的含噪的图像数据，

为系统编码矩阵，(·)^T表示矩阵的转置，

是在k空间位置中获取的算子，I_C是C×C的单位矩阵，

是由N×N的单位矩阵I_N的M行构成的亚采样矩阵，M＜N，M表示单位矩阵I_N中的M行，即采样数据的点数，N＝m×n表示二维图像的像素点个数，也表示所有k空间位置的总个数，m和n分别为单个线圈二维图像的行数和列数，C表示线圈图像的总个数，

表示Kronecker积，

是二维离散傅里叶算子，F_n分别是m点和n点傅里叶变换矩阵，

为灵敏度信息算子组矩阵，其中S_cj表示第c个线圈的第j组灵敏度信息算子的灵敏度，J表示灵敏度信息算子的总组数，c为线圈索引变量，j为灵敏度信息算子组数索引变量，

为亚采样数据，y_c表示第c个线圈的亚采样数据，L为

的梯度的Lipschitz常数；

Step2：利用MM算法将中间变量zⁱ⁺¹进行去噪，得到第i+1次循环的重构图像xⁱ⁺¹，计算公式如下：

其中，

表示待去噪的重构图像，θ_j表示第j组灵敏度信息算子的待去噪的重构图像，

表示Lp伪范数联合总变分正则项，

表示向量变量

的L₂范数的p次幂，r表示k空间位置索引变量，p是一个常数，0＜p＜1，

是一阶有限差分变换，

表示待去噪的重构图像θ的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，

和

分别为行和列方向的一阶有限差分变换，D_n和D_m分别为n×n和m×m的循环矩阵，其结构为

I_n和I_m分别为n×n和m×m的单位矩阵，λ＝α/L，λ表示正则化参数，α为平衡数据保真项和正则项的参数；

Step3：利用FISTA加速迭代法对步骤Step2中得到的重构图像xⁱ⁺¹进行加速；

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

Step5：输出重构图像x＝xⁱ⁺¹。

Step2的具体过程如下：

Step20：初始化，令θ⁰＝0，k＝1，

其中，θ⁰表示待去噪的重构图像的初始值，k表示迭代次数；

Step21：初始化，r＝1；

Step22：计算

式中，

表示第k次迭代的第r个k空间位置的总变分算子的辅助变量，

表示第k-1次迭代的待去噪的重构图像，

表示第k-1次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第k-1迭代的待去噪的重构图像θ^k-1的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，β表示惩罚项的参数；

Step23：判断r的值是否大于N；若不大于则对r做加一操作之后返回步骤Step22，否则进入步骤Step24；

Step24：初始化，j＝1；

Step25：计算：

式中，

表示第k次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第i+1次外循环的第j组灵敏度信息算子的中间变量，

表示第k次迭代的所有k空间位置的第j组灵敏度信息算子的总变分算子的辅助变量，F和F^-1分别为二维离散傅里叶变换和二维离散傅里叶逆变换，λ＝α/L表示正则化参数，α为平衡数据保真项和正则项的参数；

Step26：判断j的值是否大于J，若不大于，则对j做加一操作之后返回步骤Step25；否则，进入步骤Step27；

Step27：判断k值的大小是否大于等于最大迭代次数K，若不大于，则对k做加一操作之后返回步骤Step21；否则，进入步骤Step28；

Step28：输出xⁱ⁺¹＝θ^K，xⁱ⁺¹表示第i+1次循环的重构图像，θ^K是完成第K次迭代后的所有灵敏度信息算子的去噪后的重构的图像。

Step3的具体过程如下：

Step31：更新tⁱ⁺¹，计算公式如下：

其中，tⁱ、tⁱ⁺¹分别表示第i次和第i+1次循环的时间加速因子；

Step32：更新wⁱ⁺¹，计算公式如下：

其中，wⁱ⁺¹为第i+1循环的含噪的图像数据，xⁱ表示第i次循环的重构图像。

本发明的有益效果是：基于ESPIRiT重建模型，本发明提出了一种采用LpJTV正则化项的高质量并行成像重建算法。首先，我们将约束的非凸重构问题，转化为无约束的凸优化问题，再对拟合正则项和稀疏约束正则项进行简化，再采用OS技术将简化后的重构问题转化成一个梯度计算问题和一个去噪问题，其次，去噪问题应用MM算法进行求解。最后，再通过FISTA进行对新的算法进行加速。理论分析和成像实验表明，与传统的重构算法相比，我们提出的新的并行成像算法可以更有效地提高重构图像的质量。

附图说明

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

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

图3为使用了八通道的头部线圈进行完全采样的来获取受试者的体内人脑数据(即data1)图；

图4为使用FOV小于受试者的头部数据(即data2)图；

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

图6为数据集data1下，使用OSL1重建并行MR图像(具有6×加速度和24×24ACS线)；

图7为数据集data1下，使用OSLpJTV重建并行MR图像(具有6×加速度和24×24ACS线)；

图8为数据集data2下，使用OSL1重建并行MR图像(具有6×加速度和24×24ACS线)；

图9为数据集data2下，使用OSLpJTV重建并行MR图像(具有6×加速度和24×24ACS线)；

图10为数据集data1下，使用OSL1重建的误差图(具有6×加速度和24×24ACS线)；

图11为数据集data1下，使用OSLpJTV重建的误差图(具有6×加速度和24×24ACS线)；

图12为数据集data2下，使用OSL1重建的误差图(具有6×加速度和24×24ACS线)；

图13为数据集data2下，使用OSLpJTV重建的误差图(具有6×加速度和24×24ACS线)。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：如图1-2所示，在ESPIRiT模型中，图像每一个空域点r处的灵敏度信息算子s_r通常可通过求解以下问题得到：

式中，

是图像空间中的每个位置r的半正定矩阵值的卷积，s_r是

的特征值为“1”的特征向量，r表示k空间位置索引变量。

的定义为：

其中，

是针对向量化的多线圈图像进行二维(twodimensional，2D)傅里叶变换的矩阵，F_m和F_n分别是m点和n点傅里叶变换矩阵，I_C是C×C的单位矩阵，N＝m×n表示2D图像的像素点个数，也表示所有k空间位置的总个数，m和n分别为单个线圈2D图像的行数和列数,C表示线圈图像的总个数，

表示Kronecker积。其中

是具有矩阵值内核的卷积，

按下式定义：

其中

是R_r选择的每个k空间数据块的采样数，R_r表示k空间位置r的邻域的整个网格在所有线圈中选择的一个k空间块，(·)^H表示矩阵的共轭转置；算子

的计算(在所有的k空间位置r上)与V_||的每个内核是相关的，它可以表示为一组卷积。V_||是矩阵A的行空间跨度。矩阵A是通过在整个自动校准(Automatic calibration，AC)数据中滑动窗口来构造的校准矩阵。

Uecker等证明了校准矩阵A存在零空间向量(即：校准矩阵A存在冗余信息)，也就是说k空间块之间是相关的，可以通过k空间块之间的相关性来重构缺失数据。校准矩阵A的奇异值分解为：

A＝UΣV^H (4)

式(4)中矩阵V的列是A的行的基。在校准数据中所有的空间块都存在一个基，因此，可以将V分解为A的零空间跨度V_⊥和A的行空间跨度V_||；其中U，Σ，V是校准矩阵A的奇异值分解因子。

在一下情况下，采样误差会产生一些小于“1”特征值的特征向量，这表明信号分量不能用严格的SENSE模型来解释。

为此，Uecker等提出以下模型来获取多组灵敏度信息算子：

其中，J表示灵敏度信息算子的总组数，通常为1或2；l _j(r)是近似为“1”的常数；s_j(r)为空间位置为r的第j组灵敏度信息算子的灵敏度，j表示第j组灵敏度信息算子组数索引变量。将所有位置的s_j(r)重新排列可得对角矩阵

于是，可得基于ESPIRiT的重构模型：

其中，

是由N×N单位矩阵的M行构成的亚采样矩阵，相当于只选取k空间上相应位置的数据，M＜N，M表示单位矩阵I_N中的M行，即采样数据的点数；

表示2D傅里叶变换；S_cj表示第c个线圈图像的第j组灵敏度信息算子的灵敏度，c为线圈索引变量；

为第c个线圈的亚采样数据；

第j组灵敏度信息算子的重构的图像。

式(6)可以用最小二乘法求解。

Uecker利用正则项技术，将式(6)问题转化成如下优化问题：

其中Ψ是小波变换，||·||₁表示L₁-范数，

表示L₂-范数的平方。α为平衡数据保真项和正则项的参数。

为了表达的方便，将ESPIRiT的重构模型表示成更为紧凑的形式：

其中

是k空间位置获取的算子，

是2D傅里叶变换，

表示重构的图像，x_j表示第j组灵敏度信息算子的重构的图像，

为亚采样数据，y_c表示第c个线圈的亚采样数据，(·)^T表示矩阵的转置。而

表示灵敏度信息算子组矩阵，由所有S_cj构成：

为了提高重构质量，考虑x的J组灵敏度信息算子的Joint TV(JTV)正则项，再使用Lp伪范数(0＜p＜1)代替式(7)中L₁范数，则可得优化问题：

其中，

表示Lp伪范数联合总变分正则项，

表示向量变量

的L₂范数的p次幂，p是一个常数，0＜p＜1，

是一阶有限差分变换，

表示重构图像x的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，

和

分别为行和列方向的一阶有限差分变换，I_n和I_m分别为n×n和m×m的单位矩阵，D_n和D_m分别为n×n和m×m的循环矩阵，具有如下结构：

根据

式(10)可以重写为：

其中，

表示待去噪的重构图像，θ_j表示第j组灵敏度信息算子的待去噪的重构图像，θ和x规模相同。

表示待去噪的重构图像θ的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，

引入算子分裂技术(Operator Splitting，OS)，式(12)可以转换为：

式(13)中，zⁱ⁺¹表示第i+1次循环的中间变量，wⁱ为第i次循环的含噪的图像数据，i表示循环变量。L为

的梯度的Lipschitz常数，λ＝α/L表示正则项的参数。

式(14)中，xⁱ⁺¹表示第i+1次循环的重构图像，

使用了一种快速MM算法来求解(14)。考虑隐函数ψ的加性半二次优化，

可由下式表示：

其中，

和

是辅助向量，β是惩罚项的参数，ψ(s)是依赖于

的矩阵函数。式(15)可以通过如下个公式求解：

因此，令辅助变量

根据公式(15)可以将式(14)可以重写为：

其中

是总变分算子的辅助变量，u_r表示第r个k空间位置的总变分算子的辅助变量。问题(17)的解包含分别最小化θ和u。利用增强拉格朗日乘子法，可以将最小化问题(17)分解为两个子问题求解，如：

根据式(16)，可以得到(18)的解：

式中，

表示第k次迭代的第r个k空间位置的总变分算子的辅助变量。

表示第k-1次迭代的待去噪的重构图像，

表示第k-1次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第k-1次迭代的待去噪的重构图像θ^k-1的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，k表示迭代次数。

式(19)的求解公式如下：

其中，

是第k次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第i+1次循环的第j组灵敏度信息算子的中间变量，

表示第k次迭代的所有k空间位置的第j组灵敏度信息算子的总变分算子的辅助变量。

是一个在周期边带条件下的块循环矩阵，并且可以通过2D离散傅里叶变换对角化。因此，式(21)可以转化为：

式中，F和F^-1分别为2D离散傅里叶变换和2D离散傅里叶逆变换。

用MM算法经过K次迭代后(K为最大迭代次数)，得到了所有灵敏度信息算子的去噪后的重构图像θ^K，令xⁱ⁺¹＝θ^K，xⁱ⁺¹是第i+1循环的重构图像。利用FISTA迭代加速法对第i+1循环的重构图像xⁱ⁺¹进行加速，公式如下：

其中，t为时间加速因子，tⁱ、tⁱ⁺¹分别表示第i次和第i+1次循环的时间加速因子；w^{i +1}为第i+1次循环的含噪的图像数据，xⁱ表示第i次循环的重构的图像。

因此，针对公式(10)的求解，本发明提出一种有效的重构方法——一种基于特征向量的自一致性联合全变分Lp伪范数的并行磁共振成像重构算法(Operator Splittingmethod for solving the ESPIRiT Parallel MR Imaging Reconstruction problemwith the LpJTV regularization term，OSLpJTV)，具体流程如图1所示。

步骤如下：

Step0：初始化，令w¹＝0，x¹＝0，t¹＝0，i＝1；

其中，w¹表示含噪的图像数据的初始值；x¹是重构图像的初始值；t¹表示时间加速因子的初始值；

Step1：计算中间变量zⁱ⁺¹，计算公式如下：

Step2：利用MM算法将中间变量zⁱ⁺¹进行去噪，得到第i+1次循环的重构图像xⁱ⁺¹。

Step3：利用FISTA加速迭代法对步骤Step2中得到的重构图像xⁱ⁺¹进行加速，其子步骤如下：

Step31：更新tⁱ⁺¹，计算公式如下：

Step32：更新wⁱ⁺¹，计算公式如下：

Step4：判断是否达到停止准则，若满足条件则进入步骤Step5，否则返回步骤Step1；停止准则为迭代达到最大迭代次数或者当重构图像与前次迭代重构出的图像的相对误差小于某值时；

Step5：输出重构图像x＝xⁱ⁺¹。

其中，含LpTV正则项的去噪问题可以通过步骤Step2用MM算法(MajorizationMinimization method for the denoising problem with the LpJTV regularizationterm，MMLpJTV)进行求解，流程如图2所示。

步骤如下：

Step20：初始化，令θ⁰＝0，k＝1，

其中，θ⁰表示待去噪的重构图像的初始值；

Step21：初始化，r＝1；

Step22：计算

Step23：判断r的值是否大于N；若不大于则对r做加一操作之后返回步骤Step22，否则进入步骤Step24；

Step24：初始化，j＝1；

Step25：计算：

Step26：判断j的值是否大于J，若不大于，则对j做加一操作之后返回步骤Step25；否则，进入步骤Step27；

Step27：判断k值的大小是否大于等于最大迭代次数K，若不大于，则对k做加一操作之后返回步骤Step21；否则，进入步骤Step28；

Step28：输出xⁱ⁺¹＝θ^K，xⁱ⁺¹表示第i+1次循环的重构图像，θ^K是完成第K次迭代后的所有灵敏度信息算子的去噪后的重构的图像。

在下面的实验中，比较了提出的OSLpJTV算法(问题(10))和传统算法OSL1(问题(7))的算法性能。实验中的所有算法都是使用Matlab(MathWorks，Natick，MA)实现的。

为了测试所有比较的算法的性能，本发明使用了八通道的头部线圈进行完全采样的来获取受试者的体内人脑数据，即data1；使用FOV小于受试者的头部数据的data2；如图3和图4所示。为了生成测试数据集，使用加速因子为R×(不包括ACS线)的笛卡尔泊松亚采样，对全采样的数据集进行人工采样并用于仿真重建，如图5所示。在以下实验中，对于所有比较的算法，使用尺寸为24×24，且ESPIRiT内核尺寸为6×6的校准区域。

所有的实验都是在配置Intel Core i5-3230M@2.6GHz CPU，8GB内存和Windows10操作系统(64位)的笔记本电脑上进行的。比较算法的正则化参数是针对SNR最优进行调整的。在以下实验中，SNR和NRMSE用于定量评估重建图像的质量，分别定义为：

MSE表示通过使用来自全采样数据集的平方和的平方根(SRSoS)方法计算的重建图像和参考图像之间的均方误差，Var表示参考图像的方差。

其中x是参考图像，是

重构图像。

首先，对加速因子为6时，两个数据集下的各个算法的重构图像进行了视觉上的比较，如图6-图9所示。图6、图7展示了在数据集data1下，算法OSL1和算法OSLpJTV的重构图像。从图6可以看出，使用OSL1算法的重构图像有着明显的模糊伪影。图7使用LpJTV正则项的OSLpJTV算法的重构图像质量明显优于使用L1正则项的OSL1算法的重构质量。同样，在数据集data2下，图8和图9展示了两种算法的重构图像，从图8和图9中可以看出OSLpJTV算法重构的图像更清晰，并且能够有效的提高图像的重构质量。

为了进一步比较提出的新算法的重构性能，在图10-图13中展示了两种算法在数据集data1和data2下的误差图像(白点越多，误差越大)。从图10和图12可以明显看出，OSL1重构图像有较大的误差。而图11和图13能更有效地改善重构误差。综上所述，在两个测试集中，使用LpJTV正则项的OSLpJTV算法总是能够获得质量最好的重构图像，而OSL1重构图像质量较差。

表1对于data1和data2，加速因子为3-10的笛卡尔亚采样重构的两种比较方法的SNR和NRMSE

除了在视觉上比较重建结果外，还定量评估了这些方法在data1和data2的不同加速因子(3-10)下的笛卡尔亚采样重建MR图像的性能。表1给出了关于加速因子的重建结果的SNR和NRMSE。对于SNR，算法的值越高，重建性能越好。对于NRMSE，算法的值越低，重建性能越好。如表1所示，在加速因子3-10时，对于data1，提出的方法(OSLpJTV)与传统方法(OSL1)相比具有2dB的平均SNR增益。对于data2，提出的方法(OSLpJTV)与传统方法(OSL1)相比具有1.67dB的平均SNR增益。

同样，可以在表1中观察两种方法的NRMSE的相同对比。容易观察到，随着加速因子的增加，提出的方法重建的精度更低。与视觉比较类似，也可以得出同样的结论，与传统方法(OSL1)相比，所提出的方法(OSIRJTV)表现出更优越的性能。

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

Claims (3)

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

Step0：初始化，令w¹＝0，x¹＝0，t¹＝0，i＝1；

其中，w为含噪的图像数据，w¹表示含噪的图像数据的初始值；x表示重构的图像，x¹是重构图像的初始值；t为时间加速因子，t¹表示时间加速因子的初始值；i表示循环变量；

Step1：计算中间变量zⁱ⁺¹，计算公式如下：

式中，zⁱ⁺¹为第i+1次循环的中间变量，wⁱ为第i次循环的含噪的图像数据，

为系统编码矩阵，(·)^T表示矩阵的转置，

是在k空间位置中获取的算子，I_C是C×C的单位矩阵，

是由N×N的单位矩阵I_N的M行构成的亚采样矩阵，M＜N，M表示单位矩阵I_N中的M行，即采样数据的点数，N＝m×n表示二维图像的像素点个数，也表示所有k空间位置的总个数，m和n分别为单个线圈二维图像的行数和列数，C表示线圈图像的总个数，

表示Kronecker积，

是二维离散傅里叶算子，F_n分别是m点和n点傅里叶变换矩阵，

为灵敏度信息算子组矩阵，其中S_cj表示第c个线圈的第j组灵敏度信息算子的灵敏度，J表示灵敏度信息算子的总组数，c为线圈索引变量，j为灵敏度信息算子组数索引变量，

为亚采样数据，y_c表示第c个线圈的亚采样数据，L为

的梯度的Lipschitz常数；

Step2：利用MM算法将中间变量zⁱ⁺¹进行去噪，得到第i+1次循环的重构图像xⁱ⁺¹，计算公式如下：

其中，

表示待去噪的重构图像，θ_j表示第j组灵敏度信息算子的待去噪的重构图像，

表示Lp伪范数联合总变分正则项，

表示向量变量

的L₂范数的p次幂，r表示k空间位置索引变量，p是一个常数，0＜p＜1，

是一阶有限差分变换，

表示待去噪的重构图像θ的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，

和

分别为行和列方向的一阶有限差分变换，D_n和D_m分别为n×n和m×m的循环矩阵，其结构为

I_n和I_m分别为n×n和m×m的单位矩阵，λ＝α/L，λ表示正则化参数，α为平衡数据保真项和正则项的参数；

Step3：利用FISTA加速迭代法对步骤Step2中得到的重构图像xⁱ⁺¹进行加速；

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

Step5：输出重构图像x＝xⁱ⁺¹。

2.根据权利要求1所述的基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法，其特征在于，所述Step2的具体过程如下：

Step20：初始化，令θ⁰＝0，k＝1，

其中，θ⁰表示待去噪的重构图像的初始值，k表示迭代次数；

Step21：初始化，r＝1；

Step22：计算

式中，

表示第k次迭代的第r个k空间位置的总变分算子的辅助变量，

表示第k-1次迭代的待去噪的重构图像，

表示第k-1次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第k-1迭代的待去噪的重构图像θ^k-1的一阶有限差分变换，

是提取

的所有k空间位置为r的点组成的列向量，β表示惩罚项的参数；

Step23：判断r的值是否大于N；若不大于则对r做加一操作之后返回步骤Step22，否则进入步骤Step24；

Step24：初始化，j＝1；

Step25：计算：

式中，

表示第k次迭代的第j组灵敏度信息算子的待去噪的重构图像，

表示第i+1次外循环的第j组灵敏度信息算子的中间变量，

表示第k次迭代的所有k空间位置的第j组灵敏度信息算子的总变分算子的辅助变量，F和F^-1分别为二维离散傅里叶变换和二维离散傅里叶逆变换，λ＝α/L表示正则化参数，α为平衡数据保真项和正则项的参数；

Step26：判断j的值是否大于J，若不大于，则对j做加一操作之后返回步骤Step25；否则，进入步骤Step27；

Step27：判断k值的大小是否大于等于最大迭代次数K，若不大于，则对k做加一操作之后返回步骤Step21；否则，进入步骤Step28；

Step28：输出xⁱ⁺¹＝θ^K，xⁱ⁺¹表示第i+1次循环的重构图像，θ^K是完成第K次迭代后的所有灵敏度信息算子的去噪后的重构的图像。

3.根据权利要求1所述的基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法，其特征在于，所述Step3的具体过程如下：

Step31：更新tⁱ⁺¹，计算公式如下：

其中，tⁱ、tⁱ⁺¹分别表示第i次和第i+1次循环的时间加速因子；

Step32：更新wⁱ⁺¹，计算公式如下：

其中，wⁱ⁺¹为第i+1次循环的含噪的图像数据，xⁱ表示第i次循环的重构图像。

Images