CN110148215B - 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法 - Google Patents

Abstract

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，属于四维磁共振图像重建领域，为了解决现有重建方法存为的重建图像伪影严重以及重建时间较长的问题。包括：将欠采样K空间数据进行三维逆傅里叶变换得到初始重建图像；建立图像重建模型；引入辅助变量和正则化参数；获得三个子问题；初始化内、外循环最大迭代次数和收敛阈值；在内循环中迭代求解重建图像；判断重建图像是否满足收敛条件，满足则退出内循环，否则返回内循环继续迭代；未达到外循环最大迭代次数则更新正则化系数后再进入内循环迭代，否则终止，获得最终重建图像。实验表明本发明能有效加速磁共振成像且能较好地恢复图像细节，并获得具有临床诊断意义的重建图像。

Description

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重 建方法

技术领域

本发明属于磁共振成像领域的四维磁共振图像重建领域，具体涉及一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法。

背景技术

四维磁共振成像技术(Four-dimensional Magnetic Resonance Imaging,4DMRI)可以同时展现受检组织器官的三维解剖结构信息和动态变化信息。例如，四维心脏磁共振(4D CMR)成像可提供可视化三维心脏解剖结构和心脏瓣膜运动信息，为心血管疾病的临床诊断以及房室结构和心肌功能评价提供帮助；四维放射治疗(4D Radiation Therapy,4DRT)借助四维磁共振图像可实现对肿瘤的四维跟踪，避免对健康组织有所损害。

四维磁共振数据采集时需要获取包含几乎连续时刻上的三维数据以致所需扫描时间非常长，在长时间扫描中，患者不自主的生理运动导致图像产生运动伪影，给临床诊断带来困难。同时对一些患心脏疾病或哮喘病的患者及儿童患者，无法进行长时间扫描。所以加速成像过程成为了四维磁共振成像中亟待解决的问题。加速四维磁共振成像的关键一方面在于对硬件条件的改善和快速扫描序列的研究，一方面在于图像重建方法的研究。受磁共振扫描仪器物理特性及患者生理情况的限制，前者加速磁共振成像的发展空间有限，所以对图像重建算法的研究是解决问题的关键。目前大部分重建算法都仅在空间域或时间域进行约束，并未充分利用数据中的先验信息，因此在重建图像质量等方面存在一定缺陷。

发明内容

本发明为了解决现有重建方法图像细节的重建效果差以及重建时间较长的问题，充分挖掘了四维磁共振图像在时空域上的先验信息，提出了一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，达到了高质量重建图像细节和加速四维磁共振成像的目的。

为实现上述目的，本发明所采取的技术方案是：

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，所述四维包括三维空间域和一维时间域；所述方法包含以下步骤：

(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据b；

(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X⁽⁰⁾＝F^-1b，F表示傅里叶变换操作；

(3)结合四维磁共振图像的平滑特性和局部化低秩特性建立如下图像重建模型：

式中，X^*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项

为数据保真项，b表示测量数据；第二项/>

表示对空间域及时间域的平滑约束，采用l₁范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||_*表示核范数；其中P_i为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用/>

表示分组构建局部化低秩矩阵的操作算子；λ₁、λ₂分别表示平衡第二项、第三项的约束项系数；

(4)引入辅助变量T和W_i，将式(1)表示的重建模型改写为如下易于求解的形式：

利用拉格朗日乘子法，将其写为如下无约束形式：

其中β₁和β₂为正则化参数；

(5)本发明采用内外双重循环的模式；初始化外循环最大迭代次数，内循环最大的迭代次数和收敛阈值；将上述式(3)所示问题分为X、T和W_i三个子问题，在内循环中采用交替的方法迭代求解X、T和W_i三个子问题，在外循环中更新β₁和β₂；

得到第k次迭代结果为：

(6)判断本次迭代重建图像X^(k)是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环；

(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β₁和β₂，重置内循环迭代次数为k＝1且将本次重建图像作为X⁽⁰⁾，返回步骤(5)。

进一步地，在步骤(1)中，所述参考四维磁共振图像大小为：D1×D2×D3×N_t；D1、D2、D3分别表示三个空间维度尺寸，N_t表示时间维的大小；

所述欠采样模板为径向SoS采样模板(采样模板采用如图2所示的SoS方式)，该采样模板在k_x-k_y平面进行径向采样，在k_z方向进行笛卡尔采样；数据采集过程中先采集所有k_x-k_y平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°(黄金角度)，重复此过程；所需采集的辐条数量根据欠采样倍数确定；k_x、k_y、k_z表示笛卡尔坐系中K空间的三个坐标轴。

进一步地，在步骤(2)中，只在空间维度上进行逆傅里叶变换。

进一步地，在步骤(3)中，在建立的基于平滑约束和局部低秩约束的四维磁共振图像重建模型中，平滑约束项中D(X)定义为如下加权梯度的形式：

其中，由于四维磁共振图像在时间域和空间域上的平滑性不同，且表现为时间域的平滑性强于空间域，所以引入了参数γ增强了在时间维度上的惩罚，γ>1；

x、y、z、t表示在笛卡尔坐标系下的四维图像坐标。

进一步地，在步骤(3)中(按照如图3所示的方法构建局部化低秩矩阵)，

小图像块的划分方式：设置空间域每个分块尺寸为M×N×Q和扫描步长S₁,S₂,S₃，其中1<S₁<M,1<S₂<N,1<S₃<Q，将每一时刻图像按相同的分块尺寸和扫描步长划分为若干个部分重叠的小图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到图像块数量，由于对三维图像的分块方式是相同的，所以每时刻的分块数量P也相同的；将所有时刻对应位置上大小为M×N×Q小图像块视为一组，共P组，

分组构建局部化低秩矩阵的操作算子

具体为：将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×t矩阵，获得P组低秩矩阵；由于局部像素之间的相关性，所以每组矩阵都具有低秩性。

进一步地，在步骤(5)中，为了使式(3)与式(1)所示问题等价，在外循环中，β₁与β₂的值需要不断增大，达到10⁶数量级。

进一步地，在步骤(5)中，X、T和W_i三个子问题解的获得过程为：

通过保持其中两个变量不变，更新另一个变量的方式来交替求解式(3)，因此包括以下三个步骤：

步骤5-1，固定X和W_i，求解T，得到如下T子问题：

利用软阈值法求解得：

步骤5-2，固定X和T，求解W_i，得到如下W_i子问题:

上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：

考虑到每小块在计算上是独立的，且对

进行奇异值分解有：

则上式用软阈值法求解，得：/>

步骤5-3，固定W_i和T，求解W_i，得到如下X子问题：

分析上式的第三项，其是对各分块矩阵差值的平方和，所以可以利用权重矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即

则式(13)等价于如下形式：

显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解；

(2A^TA+λ₁β₁D^TD+λ₂β₂)X＝2A^Tb+λ₁β₁D^T(T)+λ₂β₂W (15)

同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：

由式(9)得到迭代形式的公式(4)，由式(12)得到迭代形式的公式(5)，由式(16)得到迭代形式的公式(6)。

进一步地，步骤(6)中包含以下操作步骤：

步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代；

步骤6-2，当连续两次迭代的相对误差小于允许范围时，算法收敛，退出内循环，收敛条件为：

|ε^(k-1)-ε^(k)|<0.005ε^(k) (17)

ε表示重建误差。

进一步地，步骤(7)中包含以下操作步骤：

步骤7-1，判断是否达到最大外循环迭代次数，是则退出外循环，终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β₁和β₂，如式(18)和式(19)所示：

β₁＝β₁·β_inc1 (18)

β₂＝β₂·β_inc2 (19)

且为了获得与原问题式(1)相同的解，在外循环中，β₁与β₂的值需要不断增大，达到10⁶数量级，β_inc1、β_inc2分别表示β₁、β₂的更新倍数；

步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X⁽⁰⁾，并返回步骤(5)。

本发明的有益效果是：

本发明研究在K空间数据高倍欠采样的情况下，利用少量采集数据重建高分辨率四维磁共振图像的方法。

本发明充分挖掘四维磁共振图像在时空域上的先验信息，建立并求解重建模型，包括以下步骤：利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据；将欠采样四维K空间数据在时间帧上分别进行三维逆傅里叶变换得到初始重建图像；构建局部化低秩矩阵，建立图像重建模型；引入辅助变量和正则化参数；将问题分为三个子问题，采用内外双重循坏模式，在内循环中迭代求解获得重建图像；判断重建图像是否满足收敛条件及是否达到最大内循环迭代次数，不满足收敛条件并且未达到最大内循环迭代次数则返回内循环继续迭代，否则退出内循环，进入外循环；未达到外循环最大迭代次数则更新正则化系数，并将本次重建图像作为内循环迭代初始图像再进入内循环迭代，否则算法终止并获得最终重建图像。

实验表明本发明具有能有效加速磁共振成像、较好地重建图像细节及能较好地克服成像伪影的优点，具有临床诊断意义。所提出的算法计算效率较高，CPU运行时间较短，是一种高效的算法；同时利用本发明提出的算法，可以较好地重建磁共振图像，较完整地保留图像中的细节区域，克服成像伪影。实验结果表明，与传统方法相比，本发明能够在较高加速倍数下获得高质量重建图像，同时提高了重建速度，实现了加速磁共振图像重建的目的。实验表明本发明能有效加速磁共振成像且能较好地恢复图像细节，并获得具有临床诊断意义的重建图像。

本发明的技术要点为：

(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据；

(2)将欠采样K空间四维数据在每一时间帧上进行三维逆傅里叶变换得到初始重建图像；

(3)设置分块尺寸和扫描步长，按相同的分块尺寸和扫描步长将四维磁共振图像每一时间帧上的三维图像划分为若干个部分重叠的图像块，同时获得权重系数矩阵，并将图像块构建成局部化低秩矩阵；结合图像空间域和时间域平滑先验信息和局部化低秩特性建立图像重建模型；

(4)引入辅助变量和正则化参数β₁和β₂；

(5)将问题分为易求解的三个子问题；初始化外循环最大迭代次数、内循环的最大迭代次数和收敛阈值；采用内外双重循环模式，在内循环中迭代求解获得重建图像；

(6)判断重建图像是否满足收敛条件及是否达到内循环最大迭代次数，满足收敛条件或达到内循环最大迭代次数则退出内循环，进入外循环，否则返回步骤(5)继续迭代；

(7)判断是否达到外循环最大迭代次数，若达到则获得最终重建图像，否则在外循环更新正则化系数β₁和β₂后，将本次重建图像作为内循环迭代初始图像返回步骤(5)。

附图说明

图1为本发明方法流程图(基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法)；

图2为径向SoS采样模板示意图；

图3为局部化低秩矩阵构建方法示意图；

图4为6倍欠采样下重建效果图；

图5为8倍欠采样下重建效果图；

图6为12倍欠采样下重建效果图；

图7为16倍欠采样下重建效果图。

具体实施方法

下面结合附图1至7对本发明进行详细说明。

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

(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据，用b表示；

(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X⁽⁰⁾＝F^-1b，F表示傅里叶变换操作；

(3)利用四维磁共振图像空间域、时间域平滑特性和局部化低秩特性建立如下图像重建模型；

式中，X^*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项

为数据保真项，b表示测量数据；第二项/>

表示对空间域及时间域的平滑约束，采用l₁范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||_*表示核范数；其中P_i为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用/>

表示分组构建局部化低秩矩阵的操作算子。λ₁、λ₂分别表示平衡第二项、第三项的约束项系数；

(4)引入辅助变量T和W_i，将式(1)变为如下易求解方程：

利用拉格朗日乘子法，将其写为如下无约束形式：

其中β₁和β₂为正则化参数。

(5)初始化外循环最大迭代次数Initer，内循环最大的迭代次数Outiter和收敛阈值ε；在内循环中对变量X，T和W_i进行交替求解，求得第k次迭代结果为：

/>

式(5)中U_i ^(k-1)与V_i ^(k-1)表示对

进行奇异值分解所得的奇异值向量，如：

(6)判断本次迭代重建图像X是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环，进入外循环；

(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β₁和β₂，重置内循环迭代次数为k＝1且将本次重建图像作为X⁽⁰⁾，返回步骤(5)。

步骤(1)中欠采样模板使用如图2所示的径向SoS采样模板，该采样模板在k_x-k_y平面进行径向采样，在k_z方向进行笛卡尔采样。数据采集过程中先采集所有k_x-k_y平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°(黄金角度)，重复此过程。所需采集的辐条数量根据欠采样倍数确定，最终获得的采集数据记为测量数据b，欠采样倍数ACC为全采样K空间数据量与欠采样K空间数据量的比值。

步骤(1)中，参考四维磁共振图像大小为D1×D2×D3×N_t(3D空间域+1D时间域)步骤(3)具体包括：

步骤3-1，设置分块尺寸M×N×Q和扫描步长S₁,S₂,S₃，其中1<S₁<M,1<S₂<N,1<S₃<Q，将各时刻上的三维图像划分成若干重叠的图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到的图像块数量；

步骤3-2，按照如图3所示的方法构建局部化低秩矩阵。所有时刻对应位置上的大小为M×N×Q图像块视为一组，一共P组，

再将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×T矩阵，得到P组矩阵。由于局部像素之间的相关性，所以每组矩阵都具有低秩性；

步骤3-3，建立基于平滑约束和局部低秩约束的四维磁共振图像重建模型

步骤(3)平滑约束项中D(X)定义为如下加权梯度的形式：

本发明中，步骤(4)引入了辅助变量T和W_i，获得了易求解的目标函数，并引入正则化参数β₁和β₂将其写成了无约束的形式。

本发明采用内外双重循环模式，在步骤(5)中初始化外循环的最大迭代次数，内循环的最大迭代次数和收敛阈值，并在内循环中迭代交替求解X，T和W_i。具体是通过保持其中两个变量不变，更新另一个变量的方式来交替求解，求解过程具体有以下三个步骤：

步骤5-1，固定X和W_i，求解T，得到如下T子问题:

利用软阈值法求解得：

步骤5-2，固定X和T，求解W_i，得到如下W_i子问题:

上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：

考虑到每小块在计算上是独立的，且对

进行奇异值分解有：

则上式用软阈值法求解，得：

步骤5-3，固定W_i和T，求解W_i，得到如下X子问题:

分析上式的第三项，其是对各分块矩阵差值的平方和，所以利用权重系数矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即

则上式等价于如下形式：

显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解。

(2A^TA+λ₁β₁D^TD+λ₂β₂)X＝2A^Tb+λ₁β₁D^T(T)+λ₂β₂W (14)

同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：

步骤(6)具体包括：

步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代。

步骤6-2，当连续两次迭代的相对误差小于允许范围时，退出内循环。收敛条件为：

|ε^(k-1)-ε^(k)|<0.005ε^(k) (16)

ε表示重建误差

步骤(7)中具体包括：

步骤7-1，判断是否达到最大外循环迭代次数，是则终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β₁和β₂，如式(18)和式(19)所示：

β₁＝β₁gβ_inc1 (17)

β₂＝β₂gβ_inc2 (18)

为了获得与原问题式(1)相同的解，在外循环中，β₁与β₂的值需要不断增大，最终达到10⁶数量级。

步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X⁽⁰⁾，返回步骤(5)。

本发明将通过以下实验进行有效性的说明

1.实验环境：

实验室用台式机参数：CPU为Inter(R)Core(TM)i5-4590，主频为3.30GHz，内存为16G，操作***为Win10 64位***，实验软件为MatlabR2016b。

2.实验结果及结果分析

本发明实验所用的图像数据为四维腹部磁共振图像，实验的欠采样倍数分别为6，8，12，16倍。为了验证本方法的有效性，本发明将利用信号误差比(SER)与结构性相似度(SSIM)对重建效果进行定量评估，并与四维全变分方法(4D TV)进行对比。同时，为了说明本方法的可行性，本方法重建图像将与初始重建图像进行比较。

选择一幅四维腹部磁共振图像作为参考图像，大小为：D1×D2×D3×N_t，将其每个像素值减去最小像素值后，再除以最大与最小像素值之差，得到归一化到[0，1]后的四维参考磁共振图像。之后按照如图1所示流程图进行重建实验。

本发明进行了欠采样倍数分别为6，8，12，16倍的四组实验，对实验结果进行定量评估如表1及表2所示。重建图像如图4到图7所示。

选取四维图像一帧的一个二维层面图像进行实验效果展示，图4展示了6倍欠采样下重建效果，图4(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图4(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图5展示了8倍欠采样下重建效果，图5(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图5(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图6展示了12倍欠采样下重建效果，图6(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图6(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图7展示了16倍欠采样下重建效果，图7(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图7(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。

实例一：本实例欠采样倍数为6倍，初始重建图像SER：15.74dB；SSIM：0.8137；4DTV重建图像SER：24.93dB；SSIM：0.9902；本发明重建图像SER：30.20dB，SSIM：0.9952；

实例二：本实例欠采样倍数为8倍，如表2所示初始重建图像SER：14.00dB；SSIM：0.7645；4D TV重建图像SER：22.94dB；SSIM：0.9854；本发明重建图像SER：27.88dB；SSIM：0.9922；

实例三：本实例欠采样倍数为12倍，初始重建图像SER：11.45dB；SSIM：0.6758；4DTV重建图像SER：20.52dB；SSIM：0.9737；本发明重建图像SER：25.92dB；SSIM：0.9898；

实例四：本实例欠采样倍数为16倍，初始重建图像SER：9.85dB；SSIM：0.6284；4DTV重建图像SER：18.95dB；SSIM：0.9665；本发明重建图像SER：24.12dB；SSIM：0.9859。表1重建图像信号误差比(dB)

表2重建图像结构性相似度

实验结果表明在相同的欠采样倍数下，利用本发明中提出的方法重建出的四维磁共振图像具有较高的信号误差比与结构性相似度，并且从重建图像上看，本发明方法对细节的重建效果好。相比于四维全变分方法，本发明所提方法可以在高倍欠采样的情况下获得更好的重建结果，重建图像细节更为清晰，与参考图像最为接近。

Claims (6)

1.一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，所述四维包括三维空间域和一维时间域；其特征在于，所述方法包含以下步骤：

(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据b；

(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X⁽⁰⁾＝F^-1b，F表示傅里叶变换操作；

(3)结合四维磁共振图像的平滑特性和局部化低秩特性建立如下图像重建模型：

式中，X^*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项

为数据保真项，b表示测量数据；第二项/>

表示对空间域及时间域的平滑约束，采用₁范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||_*表示核范数；其中P_i为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用/>

表示分组构建局部化低秩矩阵的操作算子；λ₁、λ₂分别表示平衡第二项、第三项的约束项系数；

(4)引入辅助变量T和W_i，将式(1)表示的重建模型改写为如下易于求解的形式：

利用拉格朗日乘子法，将其写为如下无约束形式：

其中β₁和β₂为正则化参数；

(5)采用内外双重循环的模式；初始化外循环最大迭代次数，内循环最大的迭代次数和收敛阈值；将上述式(3)所示问题分为X、T和W_i三个子问题，在内循环中采用交替的方法迭代求解X、T和W_i三个子问题，在外循环中更新β₁和β₂；

得到第k次迭代结果为：

(6)判断本次迭代重建图像X^(k)是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环；

(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β₁和β₂，重置内循环迭代次数为k＝1且将本次重建图像作为X⁽⁰⁾，返回步骤(5)；

在步骤(1)中，所述参考四维磁共振图像大小为：D1×D2×D3×N_t；

D1、D2、D3分别表示三个空间维度尺寸，N_t表示时间维的大小；

所述欠采样模板为径向SoS采样模板，该采样模板在k_x-k_y平面进行径向采样，在k_z方向进行笛卡尔采样；数据采集过程中先采集所有k_x-k_y平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°，重复此过程；所需采集的辐条数量根据欠采样倍数确定；k_x、k_y、k_z表示笛卡尔坐系中K空间的三个坐标轴；

在步骤(3)中，在建立的基于平滑约束和局部低秩约束的四维磁共振图像重建模型中，平滑约束项中D(X)定义为如下加权梯度的形式：

其中，由于四维磁共振图像在时间域和空间域上的平滑性不同，且表现为时间域的平滑性强于空间域，所以引入了参数γ增强了在时间维度上的惩罚，γ>1；

x、y、z、t表示在笛卡尔坐标系下的四维图像坐标；

在步骤(3)中，小图像块的划分方式：设置空间域每个分块尺寸为M×N×Q和扫描步长S₁,S₂,S₃，其中1<S₁<M,1<S₂<N,1<S₃<Q，将每一时刻图像按相同的分块尺寸和扫描步长划分为若干个部分重叠的小图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到图像块数量，由于对三维图像的分块方式是相同的，所以每时刻的分块数量P也相同的；将所有时刻对应位置上大小为M×N×Q小图像块视为一组，共P组，

分组构建局部化低秩矩阵的操作算子

具体为：将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×t矩阵，获得P组低秩矩阵。

2.根据权利要求1所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(2)中，只在空间维度上进行逆傅里叶变换。

3.根据权利要求2所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(5)中，为了使式(3)与式(1)所示问题等价，在外循环中，β₁与β₂的值需要不断增大，达到10⁶数量级。

4.根据权利要求3所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(5)中，X、T和W_i三个子问题解的获得过程为：

通过保持其中两个变量不变，更新另一个变量的方式来交替求解式(3)，因此包括以下三个步骤：

步骤5-1，固定X和W_i，求解T，得到如下T子问题：

利用软阈值法求解得：

步骤5-2，固定X和T，求解W_i，得到如下W_i子问题:

上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：

考虑到每小块在计算上是独立的，且对

进行奇异值分解有：/>

则上式用软阈值法求解，得：

步骤5-3，固定W_i和T，求解W_i，得到如下X子问题：

分析上式的第三项，其是对各分块矩阵差值的平方和，所以可以利用权重矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即

则式(13)等价于如下形式：

显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解；

(2A^TA+λ₁β₁D^TD+λ₂β₂)X＝2A^Tb+λ₁β₁D^T(T)+λ₂β₂W (15)

同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：

由式(9)得到迭代形式的公式(4)，由式(12)得到迭代形式的公式(5)，由式(16)得到迭代形式的公式(6)。

5.根据权利要求4所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，

步骤(6)中包含以下操作步骤：

步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代；

步骤6-2，当连续两次迭代的相对误差小于允许范围时，算法收敛，退出内循环，收敛条件为：

|ε^(k-1)-ε^(k)|＜0.005ε^(k) (17)

ε表示重建误差。

6.根据权利要求5所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，

步骤(7)中包含以下操作步骤：

步骤7-1，判断是否达到最大外循环迭代次数，是则退出外循环，终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β₁和β₂，如式(18)和式(19)所示：

β₁＝β₁·β_inc1 (18)

β₂＝β₂·β_inc2 (19)

且为了获得与原问题式(1)相同的解，在外循环中，β₁与β₂的值需要不断增大，达到10⁶数量级，β_inc1、β_inc2分别表示β₁、β₂的更新倍数；

步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X⁽⁰⁾，并返回步骤(5)。

Images