CN104248437A - 动态磁共振图像成像方法和系统 - Google Patents

Abstract

本发明提供了一种动态磁共振图像成像方法和系统，用以在常规的PS方法的基础上对空间基进行一范数的稀疏约束，从而更有效地减少了图像伪影和重建噪声，大大提高图像质量，有效实现高时空分辨的动态磁共振图像的重建。其方法包括：获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；根据所述导航数据，获得所述时间基函数；基于所述动态成像重建模型，利用所述时间基函数和动态图像数据，采用最小二乘法估计所述空间基函数；利用所述时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

Description

动态磁共振图像成像方法和系统

技术领域

本发明涉及动态磁共振图像成像技术，特别是涉及一种动态磁共振图像成像方法和系统。

背景技术

基于磁共振成像的无辐射，高软组织对比及多平面成像使得它成为目前最为流行的临床诊断工具之一。然而，由于奈奎斯特采样理论(即Nyquist采样定理)及硬件系统等条件限制，在获取动态MRI信号时，往往只能获取有限的一部分K空间信号。比如像心脏成像，动态对比度增强的肝脏成像，血流成像等，当用这些有限采集到的数据直接去傅里叶重建就会产生很大的运动伪影，因此就需要根据动态磁共振数据在时间维度上的低秩和稀疏特性，对其进行限制约束，去恢复出未采集数据来实现高分辨动态磁共振(MRI)成像。

对于像心脏，血流这种动态MR(Magnetic Resonance)成像需要在很短的时间内获取一个层面的切片才能避免图像出现严重的运动伪影，但是由于MR成像的本身的物理因素限制，只能在短的时间内扫描一个层面的某些行，所以在这种高度稀疏采样的数据中恢复出高时空分辨的动态图像，是目前研究动态磁共振成像的一个很重要的方向。在过去的几十年，很多理论算法被提出去解决这种稀疏重建问题，这些大致可以分为两类：一是稀疏约束方法，比如压缩感知(Compressed Sensing:CS)和它的各种变化模型,压缩感知重建是利用随机测量矩阵可把一个稀疏的高维信号随机投影到低维空间上，这样的随机投影包含了待重建信号的足够信息，这样利用压缩传感理论就可以设计随机测量矩阵在K空间中稀疏采样，只须采集少量数据就可以得到高质量的MR图像。二是依据时空相关的部分可分离(partial separability:PS)函数方法以及它的各种改进，部分可分离函数模型是通过将动态磁共振信号分解为与空间和时间相关的函数参数进行成像，这两者是相互独立的，这样就克服了传统磁共振成像时间和空间分辨率此消彼长的制约关系，还能同时实现高时间和高空间分辨率成像，非常适用于心脏，血流等动态成像研究。这两种方法各有优缺点。

虽然部分可分离函数算法和压缩传感理论能够证明磁共振成像在稀疏采样条件下可以实现高质量图像重建，从理论上突破了Nyquist采样定理的限制，但是在磁共振成像的临床应用中还存在一定的困难。由于这两种稀疏成像模型侧重点各有不同，部分可分离函数算法侧重于K-t空间数据的稀疏性，并且常规的部分可分离函数算法算法无法准确计算高度欠采样的图像信号，重建图像往往会伴随一定的噪声；而压缩传感理论则侧重于K空间数据的稀疏性，并且模型参数较多，鲁棒性不高等缺点。

基于现有技术中存在的问题，有待进一步地改进。

发明内容

基于此，有必要针对现有技术中的问题，提供一种动态磁共振图像成像方法和系统，用以在常规的PS方法的基础上对空间基进行一范数的稀疏约束，从而更有效地减少了图像伪影和重建噪声，大大提高图像质量，有效实现高时空分辨的动态磁共振图像的重建。

本发明提供的一种动态磁共振图像成像方法，其包括：

获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；

基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

根据所述导航数据，获得所述时间基函数；

基于所述动态成像重建模型，利用所述时间基函数和动态图像数据，采用最小二乘法估计所述空间基函数，在所述利用最小二乘法进行估计时，对所述空间基函数进行L₁范数的约束，获得所述空间基函数；

利用所述时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

在其中一个实施例中，所述采用最小二乘法估计所述空间基函数时，最小二乘法估计的优化判据表示为下述公式(1)：

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵；λ是正则化参数；d是所述动态图像数据；表示所述空间基形变矩阵的估计值；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；

根据所述优化判据获得对所述空间基形变矩阵的估计值，利用所述空间傅里叶变化矩阵求解获得所述空间基函数。

在其中一个实施例中，所述空间基形变矩阵为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶。

在其中一个实施例中，在根据所述优化判据获得对所述空间基形变矩阵的估计值时，引入下述公式(2)表示的Huber函数：

其中，t表示Huber函数的输入，α表示控制参数；

将所述空间基形变矩阵U_s输入所述Huber函数后，去替换所述公式(1)优化判据中λ||U_s||₁部分内的所述空间基形变矩阵U_s，然后求解所述优化判据，获得所述空间基形变矩阵的估计值。

在其中一个实施例中，所述根据所述优化判据获得对所述空间基形变矩阵的估计值的过程包括：

步骤10，调用按照下述公式(3)构建的对所述空间基形变矩阵进行最小二乘法估计的优化判据：

$\underset{\{U_s,G\}}{\min }{\left|\right|d-\mathrm{\Ω}\left(F_s{U}_s{V}_t\right)\left|\right|}_2^2+\frac{\mathrm{\λ}}{2\mathrm{\α}}{\left|\right|U_s-G\left|\right|}_F^2+\mathrm{\λ}{\left|\right|\mathrm{vec}\left(G\right)\left|\right|}_1$    公式(3)

其中，vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，所述辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；α表示Huber函数的控制参数；|| ||_F表示取F-范数；

步骤20，基于所述控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解所述公式(3)中的辅助矩阵G和空间基形变矩阵，直至所述辅助矩阵G和空间基形变矩阵的中间量均收敛；

步骤30，逐渐减小所述控制参数α的值，基于每一次调整所述控制参数α的调整值执行所述步骤20的交替最小化求解过程，直至将所述控制参数α的值减小到零，将最后一次迭代计算获得的所述空间基形变矩阵的中间量作为所述空间基形变矩阵的估计值。

在其中一个实施例中，所述步骤20的过程包括以下步骤：

步骤21，获取控制参数α的值；

步骤22，基于所述空间基形变矩阵的初始值，按照下述公式(4)计算所述辅助矩阵G的中间量；

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||_？表示取F-范数；vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；

步骤23，基于所述辅助矩阵G的中间量，按照下述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时所述空间基形变矩阵的初始值：

其中，表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；

步骤24，判断所述辅助矩阵G的中间量和所述空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续执行步骤22至步骤23；若是，则将最后一次迭代计算获得的所述空间基形变矩阵的中间量，作为基于所述控制参数α的另一调整值、执行所述步骤20时，所采用的所述空间基形变矩阵的初始值。

在其中一个实施例中，所述步骤22的计算过程包括以下步骤：

首先，提取第l次迭代时获得的所述空间基形变矩阵的中间量中各个元素的值；

然后，判断所述空间基形变矩阵的中间量的元素值的绝对值与所述控制参数α值的大小，当所述空间基形变矩阵的中间量的元素值的绝对值小于所述控制参数α的值时，在所述辅助矩阵G中的相应位置赋值为零；当所述空间基形变矩阵的中间量的元素值的绝对值大于等于所述控制参数α的值时，在所述辅助矩阵G中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\α})$     公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的所述空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示所述Huber函数中的控制参数；

最后，按照所述赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的所述辅助矩阵G的中间量。

在其中一个实施例中，所述步骤23中按照所述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量的过程包括以下步骤：

首先，调用按照公式(7)构建的关联所述空间基形变矩阵的中间量和辅助矩阵的中间量G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\λ}}{2\mathrm{\α}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\λ}}{2\mathrm{\α}}{G}^{(l+1)}$    公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表示欠采的图像数据矩阵；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；是单位算子，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；

然后，利用共轭梯度算法迭代优化求解所述公式(7)中的所述空间基形变矩阵的中间量

基于上述方法，本发明还提供了一种动态磁共振图像成像系统，其包括：

数据提取模块，用于获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；

模型构建模块，用于基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

时间基提取模块，用于根据所述导航数据，获得所述时间基函数；

空间基提取模块，用于基于所述动态成像重建模型，利用所述时间基函数和动态图像数据，采用最小二乘法估计所述空间基函数，在所述利用最小二乘法进行估计时，对所述空间基函数进行L₁范数的约束，获得所述空间基函数；及

图像恢复重建模块，用于利用所述时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

在其中一个实施例中，所述空间基提取模块包括：

模型构建单元，用于调用按照下述公式(1)构建的所述最小二乘法估计的优化判据：

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵；λ是正则化参数；d是所述动态图像数据；表示所述空间基形变矩阵的估计值；U_s表示所述空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；

最优解计算单元，用于根据所述优化判据获得对所述空间基形变矩阵的估计值；及

空间基求解单元，用于利用所述空间傅里叶变化矩阵求解获得所述空间基函数。

在其中一个实施例中，所述模型构建单元中引入一Huber函数将所述优化判据变更为下述公式(3)：

$\underset{\{U_s,G\}}{\min }{\left|\right|d-\mathrm{\Ω}\left(F_s{U}_s{V}_t\right)\left|\right|}_2^2+\frac{\mathrm{\λ}}{2\mathrm{\α}}{\left|\right|U_s-G\left|\right|}_F^2+\mathrm{\λ}{\left|\right|\mathrm{vec}\left(G\right)\left|\right|}_1$ 公式(3)

其中，vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，所述辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；α表示Huber函数的控制参数；|| ||_F表示取F-范数；V_t表示时间基函数；

所述最优解计算单元包括：

交替最小化求解单元，用于基于所述控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解所述公式(3)中的辅助矩阵G和空间基形变矩阵，直至所述辅助矩阵G和空间基形变矩阵的中间量均收敛；和

控制参数调整单元，用于逐渐减小所述控制参数α的值，基于每一次调整所述控制参数α的调整值调用执行所述交替最小化求解单元，直至将所述控制参数α的值减小到零，将最后一次迭代计算获得的所述空间基形变矩阵的中间量作为所述空间基形变矩阵的估计值。

在其中一个实施例中，所述交替最小化求解单元包括：

控制参数获取单元，用于获取控制参数α的值；

第一计算单元，用于基于所述空间基形变矩阵的初始值，按照下述公式(4)计算所述辅助矩阵G的中间量；

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||_F表示取F-范数；vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；α表示Huber函数的控制参数；

第二计算单元，用于基于所述辅助矩阵G的中间量，按照下述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时所述空间基形变矩阵的初始值：

其中，表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；α表示Huber函数的控制参数；及

第一判断单元，用于判断所述辅助矩阵G的中间量和所述空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续调用所述第一计算单元和第二计算单元；若是，则将最后一次迭代计算获得的所述空间基形变矩阵的中间量，作为基于所述控制参数α的另一调整值、调用所述交替最小化求解单元时所采用的所述空间基形变矩阵的初始值。

在其中一个实施例中，所述第一计算单元包括：

元素值提取单元，用于提取第l次迭代时获得的所述空间基形变矩阵的中间量中各个元素的值；

第二判断单元，用于判断所述空间基形变矩阵的中间量的元素值的绝对值与所述控制参数α值的大小，当所述空间基形变矩阵的中间量的元素值的绝对值小于所述控制参数α的值时，在所述辅助矩阵G中的相应位置赋值为零；当所述空间基形变矩阵的中间量的元素值的绝对值大于等于所述控制参数α的值时，在所述辅助矩阵G中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\α})$      公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的所述空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示所述Huber函数中的控制参数；及

输出单元，用于按照所述赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的所述辅助矩阵G的中间量。

在其中一个实施例中，所述第二计算单元包括：

优解模型构建单元，用于调用按照如下公式(7)所构建的关联和G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\λ}}{2\mathrm{\α}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\λ}}{2\mathrm{\α}}{G}^{(l+1)}$    公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表示欠采的图像数据矩阵；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；是单位算子，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；和

梯度计算单元，用于利用共轭梯度算法迭代优化求解所述公式(7)中的所述空间基形变矩阵的中间量

本发明提出了一种高分辨的动态磁共振成像方法，即将基于PS的稀疏采样方法和低秩矩阵的稀疏约束技术相结合的方法，且提供了一种新的优化模型，对原有的PS成像算法做了新的改进，提高重建质量，避免了拟合不准确而产生的伪影以及对噪声重建。

附图说明

图1为本发明动态磁共振图像成像方法的流程示意图；

图2为本发明动态磁共振图像成像方法中获得空间基形变矩阵的估计值的流程示意图；

图3为本发明动态磁共振图像成像方法的一最优实施例流程图；

图4为本发明动态磁共振图像成像系统的结构示意图；

图5为本发明动态磁共振图像成像系统的一最优实施例结构示意图；

图6为针对PS图像数据帧数为20的原始数据进行图像重建的结果对比图；

图7为针对PS图像数据帧数为14的原始数据进行图像重建的结果对比图。

具体实施方式

基于动态磁共振图像成像技术，本发明基于现有的部分可分离函数模型进行改进，利用低秩矩阵恢复的原理将磁共振信号分离成有限的时间和空间基函数，并通过结合时间基函数和测量的动态图像数据，对空间基函数采用L1范数限制约束避免像TV(total variation)算子过度的平滑图像，从而避免磁共振成像中的时间和空间分辨此消彼长的制约关系，避免传统方法中利用K空间和时间的相关性来提高对丢失数据的估计时对丢失未采数据拟合不准确的问题，还可以降低对噪声对图像的影响，减少了图像伪影和重建噪声，大大提高图像质量。以下将结合附图详细说明本发明的各个具体实施方式。

如图1所示，本实施例提供了一种动态磁共振图像成像方法，其包括：

步骤100，获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据，这里的导航数据是部分可分离函数方法中采集的一部分数据，起到时间导航的作用；而动态图像数据是欠采样的图像数据；

步骤200，基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

步骤300，根据导航数据，获得时间基函数；本步骤中可以利用PS采样模式采集数据，用奇异值分解方法得到时间基函数。

步骤400，基于动态成像重建模型，利用时间基函数和动态图像数据，采用最小二乘法估计空间基函数，在利用最小二乘法进行估计时，对所述空间基函数进行L1范数的约束，获得空间基函数；

步骤500，利用时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

在本实施例中，步骤100所基于的部分可分离函数方法的采样模式具体原理如下所示。步骤100的采样模式主要是基于PS理论模型，即在运动部位(例如，心脏，冠状动脉等)的磁共振动态成像中，实际接收的信号S(k,t)和图像函数ρ(x,t)的关系为下述公式(8)所示。

$S(k,t)={\∫}_{-\∞}^{+\∞}\mathrm{\ρ}(x,t)e^{-i2\mathrm{\πk}\·x}\mathrm{dx}+\mathrm{\η}(k,t)$      公式(8)

其中，ρ(x,t)为理想的图像函数，S(k,t)是理想的k空间数据，η(k,t)是测量噪声。在磁共振动态成像中，目标就是尽量高分辨率地恢复出ρ(x,t)信号。然而对于心脏等运动部位经过一段时间的信号采集，空间中某一平面数据将包含来自运动物体各个不同部位的信息，而且如果采集的时间越长，来自不同部位的混叠信息越多，如果对这些信息不加任何处理直接进行图像重建，则重建结果将发生模糊和运动伪影。基于上述原理采集的动态图像数据实际是部分图像数据，需要通过以下步骤的图像重建来实现完整的图像数据ρ(x,t)信号，并同时降低重建结果中的伪影和噪声。

本实施例中，步骤200主要是基于PS成像算法，通过该算法根据图像函数ρ(x,t)的空间变化和时间变化性质，利用部分可分离函数的性质，将步骤100中实际采集到的信号S(k,t)分离成L阶的时间基函数和空间基函数，那么，S(k,t)的L阶PS模型可定义为下述公式(9)。

$S(k,t)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{u}_l\left(\stackrel{\→}{k}\right)v_l\left(t\right)$       公式(9)

其中，L是模阶，可以取很小的值。和分别是L阶的空间基函数和L阶的时间基函数。k是个向量。

S(k,t)可以表示成一个低阶的Casorati矩阵，如下述公式(10)。

其中，S表示Casorati矩阵；N是K空间的所有像素点数，M是K空间的时间点数，秩(S)≤L，L是动态成像重建模型中空间基函数与时间基函数的模阶，那么就可以把Casorati矩阵S简化写成如下公式(11)：

S＝U_kV_t       公式(11)

其中，空间基函数时间基函数

基于上述原理，通过步骤200构建了如上述公式(11)所示的包含空间基函数与时间基函数的动态成像重建模型，假设(k,t)空间是稀疏采样的，那么Casorati矩阵S就有很多缺失项，这里有很多种方法通过“阶”限制可以恢复矩阵S，其中的一种常用的方法就是通过导航数据得到时间基函数V_t，然后利用最小二乘法估计得到空间基函数U_k，如下述公式(12)所表示的最小二乘法估计的优化判据。

其中，d是测量的图像数据，Ω是空间稀疏采样因子，|| ||₂表示取2-范数，表示空间基函数U_k的估计值。

由于高度欠采的Casorati矩阵S是病态的，直接运用最小二乘很难进行准确的估计，因此可以对其进行稀疏约束，对空间基函数U_k和时间基函数V_t同时进行估计。而本实施例中优选的方案是根据PS模型中采集到的导航数据，先获得时间基函数V_t，然后基于动态成像重建模型(基于上述公式(11))，用已知获得得时间基函数V_t，联合PS模型采集的动态图像数据去约束，采用最小二乘法估计出空间基函数U_k。

本实施例中，为了避免直接利用最小二乘法的公式(12)的优化判据进行估计时，会将部分噪声一同重建入图像，而增加了重建后动态图像的伪影和噪声的问题，上述步骤400中在通过最小二乘法拟合对图像重建时，对空间基函数进行L1范数的约束可以在在不平滑图像信噪比的前提下，大大提高重建图像的质量。

本实施例的上述步骤500中利用时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换的过程主要是依据上述公式(8)进行逆向运算过程，从而恢复出完整的图像信号ρ(x,t)。

基于上述实施例，本实施例中步骤400中采用最小二乘法估计上述空间基函数时，具体可参见下述公式(1)所示的最小二乘法估计的优化判据。

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵，此矩阵是行为N、列为L的矩阵；λ是正则化参数；d是动态图像数据；表示空间基形变矩阵的估计值；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示动态成像重建模型中空间基函数与时间基函数的模阶；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；V_t表示时间基函数。利用最小二乘法进行估计时，用空间基形变矩阵的L₁范数乘以一正则化参数去约束最小二乘法估计的优化判据，可以让数据更便于计算或获得更加泛化的结果，提高计算机的运算速度，降低计算的复杂度。

本实施例根据上述公式(1)的优化判据获得对空间基形变矩阵的估计值，然后基于U_s＝F_s ^-1U_k利用空间傅里叶变化矩阵求解获得相应的空间基函数。基于上述(1)对空间基函数采用最小二乘拟合估计时，在高阶重建的情况下，可以降低对噪声进行重建，能够提高PS模型的拟合精度，尤其对于高阶的PS模型以及有限的数据集的情况下能得到更准确地实现数据拟合。在本实施例中，联合采用部分可分离函数方法和空间基的L1范数约束限制，可以有效地减少了图像伪影和重建噪声。本实施例引入了一个空间傅里叶变化矩阵F_s，通过引入的空间傅里叶变化矩阵F_s对步骤200分解得到的空间基函数进行形变，获得空间基形变矩阵U_s，其为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示上述动态成像重建模型中空间基函数与时间基函数的模阶；从而将利用最小二乘法拟合估计中对空间基函数的估计转化为对空间基形变矩阵U_s的估计，并结合L₁范数去约束最小二乘法估计的优化判据，从而联合了部分可分离函数和空间基的L₁范数约束限制，来有效地减少图像伪影和重建噪声，大大提高图像质量，有效实现高时空分辨的动态磁共振图像的重建。

基于上述实施例，为了降低上述公式(1)的运算复杂度，在根据上述公式(1)所示的优化判据获得对空间基形变矩阵的估计值时，引入下述公式(2)表示的Huber函数：

其中，t表示Huber函数的输入，α表示控制参数；表示Huber函数。

将空间基形变矩阵U_s输入Huber函数后，去替换上述公式(1)优化判据中λ||U_s||₁部分内的空间基形变矩阵U_s，具体参见以下公式(2-1)。

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；α表示Huber函数中的控制参数；其他元素的解释参见上述内容。

然后求解上述公式(2-1)的优化判据，获得空间基形变矩阵的估计值α是一个控制着公式(2-1)与公式(1)的最大近似程度的参数，当α趋近于0时，公式(2-1)的解就接近于公式(1)。因为上式公式(1)中的L₁范数约束项是不可微的，因此引入可微的Huber函数来降低上述公式(1)的运算复杂度。基于公式(2)所示的Huber函数来求解上述公式(1)的优化判据，可以采用很多方法，以下将提供一种最优的求解方式。

基于上述实施例，本实施例提供了一种基于半二次(half-quadratic)最小化方法来求解上述公式(1)的优化判据的实现方法。如图2所示，上述根据优化判据获得对空间基形变矩阵的估计值的过程可以包括：

步骤10，调用按照下述公式(3)构建的对空间基形变矩阵进行最小二乘法估计的优化判据：

其中，vec(G)表示将辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；α表示Huber函数的控制参数；|| ||_F表示取F-范数；V_t表示时间基函数。

步骤20，基于控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解公式(3)中的辅助矩阵G和空间基形变矩阵U_s，直至辅助矩阵G和空间基形变矩阵U_s的中间量均收敛；

步骤30，逐渐减小控制参数α的值，基于每一次调整控制参数α的调整值执行步骤20的交替最小化求解过程，直至将控制参数α的值减小到零，将最后一次迭代计算获得的空间基形变矩阵的中间量作为空间基形变矩阵的估计值通过步骤30逐步减小控制参数α的值可以实现将公式(3)的计算结果更加接近所需要对空间基形变矩阵的真实估计值。

基于上述实施例，本实施例中，上述步骤20的过程优选采用以下方式，如图3所示，其包括以下步骤：

步骤21，获取控制参数α的值，迭代过程源于控制参数α的一个初始值，可以取控制参数α可能的最大值，然后在后续过程中逐渐减小此控制参数α；

步骤22，基于空间基形变矩阵的初始值，按照下述公式(4)计算辅助矩阵G的中间量；

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||_F表示取F-范数；vec(G)表示将辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的辅助矩阵G的中间量；V_t表示时间基函数。在第一次进行迭代计算时，空间基形变矩阵的初始值为一零矩阵。

步骤23，基于辅助矩阵G的中间量，按照下述公式(5)计算此次迭代过程中的空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时空间基形变矩阵的输入：

其中，表示单位算子，N表示K空间的所有像素点数，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的空间基形变矩阵的中间量；Vt表示时间基函数。

步骤24，判断辅助矩阵G的中间量和空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续执行步骤22至步骤23；若是，则将最后一次迭代计算获得的空间基形变矩阵的中间量，作为基于控制参数α的另一调整值、执行步骤20时，所采用的空间基形变矩阵的初始值。

基于上述实施例执行的步骤20，则上述步骤30的过程优选采用如下方式，如图3所示，其包括：

在执行完步骤24之后，若辅助矩阵G的中间量和空间基形变矩阵的中间量均同时收敛，则执行步骤31，判断控制参数α的值是否为零，

若否，则执行步骤32：逐渐减小调整控制参数α的值，获得另一个控制参数α的调整值，然后返回至步骤21，开始基于每一次调整控制参数α的调整值执行步骤20的交替最小化求解过程；

若是，则视为已将控制参数α的值减小到零，执行步骤33：将最后一次迭代计算获得的空间基形变矩阵的中间量作为空间基形变矩阵的估计值

图3提供了本发明基于半二次(half-quadratic)最小化方法来求解上述公式(1)的优化判据的最佳实施方式。

基于上述最佳实施方式，上述步骤22中对辅助矩阵G中间量的计算过程可以优选包括以下步骤：

首先，提取第l次迭代时获得的空间基形变矩阵的中间量中各个元素的值；

然后，判断空间基形变矩阵的中间量的元素值的绝对值与控制参数α值的大小，当空间基形变矩阵的中间量的元素值的绝对值小于控制参数α的值时，则在辅助矩阵G中的相应位置赋值为零；当空间基形变矩阵的中间量的元素值的绝对值大于等于控制参数α的值时，则在辅助矩阵G中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\α})$     公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示Huber函数中的控制参数；

最后，按照赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的辅助矩阵G的中间量。

以上过程主要是对上述公式(4)的最优求解过程，利用求导等常规方法可以将上述公式(4)的形式看作是看作以下函数表示的内容：

$G^{(l+1)}=S\left(U_s^{\left(l\right)}\right)$       公式(6-1)

其中是一个软阈值算子，表示第l次迭代时获得的空间基形变矩阵的中间量，然后结合上述公式(2)的Huber函数可以获得以下公式(6-2)所示的变形形式。

$S{\left(Q\right)}_{n,l1}=\left\{\begin{array}{cc}0& \mathrm{if}\left|Q_{n,l1}\right|<\mathrm{\&alpha;}\\ \frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\&alpha;})& \mathrm{if}\left|Q_{n,l1}\right|\&GreaterEqual;\mathrm{\&alpha;}\end{array}\right.$      公式(6-2)

根据上述公式(6-2)的内容即可实现对第l+1次迭代时获得的辅助矩阵G的中间量G^(l+1)的最优求解，这样做可以简化计算过程，提高计算机运算速度。

基于上述最佳实施方式，上述步骤步骤23中按照公式(5)计算此次迭代过程中的空间基形变矩阵的中间量的过程可以优选包括以下步骤：

首先，调用按照如下公式(7)所构建的关联空间基形变矩阵的中间量和辅助矩阵的中间量G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{G}^{(l+1)}$       公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表欠采的图像数据矩阵；表示第l+1次迭代时获得的空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的辅助矩阵G的中间量；是单位算子，L表动态成像重建模型中空间基函数与时间基函数的模阶。

然后，利用共轭梯度算法迭代优化求解上述公式(7)中的第l+1次迭代时获得的空间基形变矩阵的中间量

在上述实施例中，在公式(7)中只有L是未知的，因此可以进行求解，由于该公式(7)的方程系数是对称矩阵，所以可以运用共轭梯度算法迭代优化求解，可以得到唯一的局部最优解从而得到最终约束求解的空间基函数U_k，最后通过多次α循环获得最终最优解。在调用上述公式(7)计算时，优选先基于上述公式(7)进行傅里叶变换的到时间基U_k，然后再基于U_s＝F_s ^-1U_k逆运算得到U_s。

基于上述一种动态磁共振图像成像方法，本实施例还提供了一种动态磁共振图像成像系统，如图4所示，其系统包括：

数据提取模块601，用于获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；

模型构建模块602，用于基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

时间基提取模块603，用于根据导航数据，获得时间基函数；

空间基提取模块604，用于基于动态成像重建模型，利用时间基函数和动态图像数据，采用最小二乘法估计空间基函数，在利用最小二乘法进行估计时，对空间基函数进行L₁范数的约束，获得空间基函数；及

图像恢复重建模块605，用于利用时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

基于上述实施例，本实施例中，如图4所示，上述空间基提取模块604优选包括：

模型构建单元614，用于调用按照下述公式(1)所构建的最小二乘法估计的优化判据：

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵；λ是正则化参数；d是动态图像数据；表示空间基形变矩阵的估计值；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示动态成像重建模型中空间基函数与时间基函数的模阶；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；

最优解计算单元624，用于根据优化判据获得对空间基形变矩阵的估计值；及

空间基求解单元634，用于利用空间傅里叶变化矩阵求解获得空间基函数。

基于上述实施例，本实施例中，如图5所示，上述模型构建单元614中引入一Huber函数将优化判据变更为下述公式(3)：

$\underset{\{U_s,G\}}{\min }{\left|\right|d-\mathrm{\&Omega;}\left(F_s{U}_s{V}_t\right)\left|\right|}_2^2+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{\left|\right|U_s-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|\mathrm{vec}\left(G\right)\left|\right|}_1$     公式(3)

其中，vec(G)表示将辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；|| ||₂表示取2-范数；|| ||₁表示取L₁范数；α表示Huber函数的控制参数；|| ||_F表示取F-范数；V_t表示时间基函数。

则上述最优解计算单元624优选包括：

交替最小化求解单元701，用于基于控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解公式(3)中的辅助矩阵G和空间基形变矩阵，直至辅助矩阵G和空间基形变矩阵的中间量均收敛；和

控制参数调整单元702，用于逐渐减小控制参数α的值，基于每一次调整控制参数α的调整值调用执行交替最小化求解单元，直至将控制参数α的值减小到零，将最后一次迭代计算获得的空间基形变矩阵的中间量作为空间基形变矩阵的估计值。

基于上述实施例，本实施例中，如图5所示，上述交替最小化求解单元701优选包括：

控制参数获取单元711，用于获取控制参数α的值；

第一计算单元721，用于基于空间基形变矩阵的初始值，按照下述公式(4)计算辅助矩阵G的中间量；

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；|| ||_F表示取F-范数；vec(G)表示将辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的辅助矩阵G的中间量；α表示Huber函数的控制参数。这里在进行第一次迭代计算时空间基形变矩阵的初始值取零矩阵。

第二计算单元731，用于基于辅助矩阵G的中间量，按照下述公式(5)计算此次迭代过程中的空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时空间基形变矩阵的初始值：

其中，表示单位算子，N表示K空间的所有像素点数，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的空间基形变矩阵的中间量；α表示Huber函数的控制参数。

和第一判断单元741，用于判断辅助矩阵G的中间量和空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续调用第一计算单元721和第二计算单元731；若是，则将最后一次迭代计算获得的空间基形变矩阵的中间量，作为基于控制参数α的另一调整值、调用交替最小化求解单元时，所采用的空间基形变矩阵的初始值。

基于上述实施例，本实施例中，如图5所示，上述第一计算单元721优选包括：

元素值提取单元71，用于提取第l次迭代时获得的空间基形变矩阵的中间量中各个元素的值；

第二判断单元72，用于判断空间基形变矩阵的中间量的元素值的绝对值与控制参数α值的大小，当空间基形变矩阵的中间量的元素值的绝对值小于控制参数α的值时，在辅助矩阵G中的相应位置赋值为零；当空间基形变矩阵的中间量的元素值的绝对值大于等于控制参数α的值时，在辅助矩阵G中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\&alpha;})$       公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示Huber函数中的控制参数；及

输出单元73，用于按照赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的辅助矩阵G的中间量。

基于上述实施例，本实施例中，如图5所示，上述第二计算单元731优选包括：

优解模型构建单元74，用于调用按照如下公式(7)所构建的关联和G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{G}^{(l+1)}$        公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表欠采的图像数据矩阵；表示第l+1次迭代时获得的空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的辅助矩阵G的中间量；是单位算子，L表示动态成像重建模型中空间基函数与时间基函数的模阶。

和梯度计算单元75，用于利用共轭梯度算法迭代优化求解上述公式(7)中的第l+1次迭代时获得的空间基形变矩阵的中间量

基于上述实施例，如图5所示，本实施例中的控制参数调整单元702可以包括以下单元：

第三判断单元712，用于在辅助矩阵G的中间量和空间基形变矩阵的中间量均同时收敛时，判断控制参数α的值是否为零，

参数输出单元722，用于当控制参数α的值不为零时逐渐减小调整控制参数α的值，获得另一个控制参数α的调整值，然后基于每一次调整控制参数α的调整值调用上述交替最小化求解单元；及

结果输出单元732，用于当控制参数α的值减小到零时，将最后一次迭代计算获得的空间基形变矩阵的中间量作为空间基形变矩阵的估计值。

以上动态磁共振图像成像系统的各个模块所实现的功能基于动态磁共振图像成像方法的各个步骤，针对上述系统中的各个功能和公式的解释参见上述动态磁共振图像成像方法中的相关步骤和公式的解释，在此不作累述。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个非易失性计算机可读存储介质(如ROM、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的系统结构和方法。

基于上述各个实施例，本发明提出了一种高分辨的动态磁共振成像方法，即将基于PS的稀疏采样方法和低秩矩阵的稀疏约束技术相结合的方法，且提供了一种新的优化模型，对原有的PS成像算法做了新的改进，提高重建质量，避免了拟合不准确而产生的伪影以及对噪声重建。本发明具有以下优点：对传统的PS进行了改进，采用稀疏约束，引入L1范数约束，把空间基的求解转化成一个线性问题求解，把存在很多个解的病态问题转化成能获得唯一解的二次函数求解。这样就能把原始的PS模型采用最小二乘拟合不准确的问题，而把噪声进行重建的问题得到解决而获得更准确的重建图像，从而获得高质量的动态重建图像。如图6和图7所示，本发明基于图3和图5的最优实施例，利用采集的大脑原始数据进行了验证，图6(a)针对PS图像数据帧数为20的原始数据、利用原始PS理论进行的重建图像，图6(b)针对PS图像数据帧数为20的原始数据、利用本发明的方法和系统所获得的重建图像；图7(a)针对PS图像数据帧数为14的原始数据、利用原始PS理论进行的重建图像，图7(b)针对PS图像数据帧数为14的原始数据、利用本发明的方法和系统所获得的重建图像。从图6和图7来看，对采集的原始数据进行本发明提出的重建方法进行离线重建，其结果确实在较高的模阶L下能够有效的抑制噪声，减少图像伪影，实现高时空分辨率的成像。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (14)

1.一种动态磁共振图像成像方法，其包括：

获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；

基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

根据所述导航数据，获得所述时间基函数；

基于所述动态成像重建模型，利用所述时间基函数和动态图像数据，采用最小二乘法估计所述空间基函数，在所述利用最小二乘法进行估计时，对所述空间基函数进行L₁范数的约束，获得所述空间基函数；

利用所述时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

2.根据权利要求1所述的动态磁共振图像成像方法，其特征在于，所述采用最小二乘法估计所述空间基函数时，最小二乘法估计的优化判据表示为下述公式(1)：

   公式(1)

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵；λ是正则化参数；d是所述动态图像数据；表示空间基形变矩阵的估计值；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_kU_k表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；||||₂表示取2-范数；||||₁表示取L₁范数；V_t表示时间基函数；

根据所述优化判据获得对所述空间基形变矩阵的估计值，利用所述空间傅里叶变化矩阵求解获得所述空间基函数。

3.根据权利要求2所述的动态磁共振图像成像方法，其特征在于，所述空间基形变矩阵为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶。

4.根据权利要求2所述的动态磁共振图像成像方法，其特征在于，在根据所述优化判据获得对所述空间基形变矩阵的估计值时，引入下述公式(2)表示的Huber函数：

   公式(2)

其中，t表示Huber函数的输入，α表示控制参数；

将所述空间基形变矩阵U_s输入所述Huber函数后，去替换所述公式(1)优化判据中λ||U_s||₁部分内的所述空间基形变矩阵U_s，然后求解所述优化判据，获得所述空间基形变矩阵的估计值。

5.根据权利要求4所述的动态磁共振图像成像方法，其特征在于，所述根据所述优化判据获得对所述空间基形变矩阵的估计值的过程包括：

步骤10，调用按照下述公式(3)构建的对所述空间基形变矩阵进行最小二乘法估计的优化判据：

$\underset{\{U_s,G\}}{\min }{\left|\right|d-\mathrm{\&Omega;}\left(F_s{U}_s{V}_t\right)\left|\right|}_2^2+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{\left|\right|U_s-G\left|\right|}_F^2+{\left|\right|\mathrm{vec}\left(G\right)\left|\right|}_1$    公式(3)

其中，vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，所述辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；||||₂表示取2-范数；||||₁表示取L₁范数；α表示Huber函数的控制参数；||||_F表示取F-范数；V_t表示时间基函数；

步骤20，基于所述控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解所述公式(3)中的辅助矩阵和空间基形变矩阵，直至所述辅助矩阵和空间基形变矩阵的中间量均收敛；

步骤30，逐渐减小所述控制参数α的值，基于每一次调整所述控制参数α的调整值执行所述步骤20的交替最小化求解过程，直至将所述控制参数α的值减小到零，将最后一次迭代计算获得的所述空间基形变矩阵的中间量作为所述空间基形变矩阵的估计值。

6.根据权利要求5所述的动态磁共振图像成像方法，其特征在于，所述步骤20的过程包括以下步骤：

步骤21，获取控制参数α的值；

步骤22，基于所述空间基形变矩阵的初始值，按照下述公式(4)计算所述辅助矩阵G的中间量；

   公式(4)

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；||||_F表示取F-范数；vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；

步骤23，基于所述辅助矩阵的中间量，按照下述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时所述空间基形变矩阵的初始值：

   公式(5)

其中，表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；V_t表示时间基函数；

步骤24，判断所述辅助矩阵的中间量和所述空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续执行步骤22至步骤23；若是，则将最后一次迭代计算获得的所述空间基形变矩阵的中间量，作为基于所述控制参数α的另一调整值、执行所述步骤20时，所采用的所述空间基形变矩阵的初始值。

7.根据权利要求6所述的动态磁共振图像成像方法，其特征在于，所述步骤22的计算过程包括以下步骤：

首先，提取第l次迭代时获得的所述空间基形变矩阵的中间量中各个元素的值；

然后，判断所述空间基形变矩阵的中间量的元素值的绝对值与所述控制参数α值的大小，当所述空间基形变矩阵的中间量的元素值的绝对值小于所述控制参数α的值时，在所述辅助矩阵中的相应位置赋值为零；当所述空间基形变矩阵的中间量的元素值的绝对值大于等于所述控制参数α的值时，在所述辅助矩阵中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\&alpha;})$    公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的所述空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示所述Huber函数中的控制参数；

最后，按照所述赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的所述辅助矩阵的中间量。

8.根据权利要求6所述的动态磁共振图像成像方法，其特征在于，所述步骤23中按照所述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量的过程包括以下步骤：

首先，调用按照公式(7)构建的关联所述空间基形变矩阵的中间量和辅助矩阵的中间量G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{G}^{(l+1)}$    公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表示欠采的图像数据矩阵；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；是单位算子，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；

然后，利用共轭梯度算法迭代优化求解所述公式(7)中的所述空间基形变矩阵的中间量

9.一种动态磁共振图像成像系统，其特征在于，所述系统包括：

数据提取模块，用于获取基于部分可分离函数方法的采样模式所采集的K-t空间导航数据和动态图像数据；

模型构建模块，用于基于部分可分离函数方法构建包含空间基函数与时间基函数的动态成像重建模型；

时间基提取模块，用于根据所述导航数据，获得所述时间基函数；

空间基提取模块，用于基于所述动态成像重建模型，利用所述时间基函数和动态图像数据，采用最小二乘法估计所述空间基函数，在所述利用最小二乘法进行估计时，对所述空间基函数进行L₁范数的约束，获得所述空间基函数；及

图像恢复重建模块，用于利用所述时间基函数和空间基函数对K空间的动态图像数据进行插值恢复，再进行傅里叶变换，从而获得重建后的动态磁共振图像。

10.根据权利要求9所述的动态磁共振图像成像系统，其特征在于，所述空间基提取模块包括：

模型构建单元，用于调用按照下述公式(1)构建的所述最小二乘法估计的优化判据：

   公式(1)

其中，Ω是采样矩阵；F_s是一空间傅里叶变化矩阵；λ是正则化参数；d是所述动态图像数据；表示所述空间基形变矩阵的估计值；U_s表示所述空间基形变矩阵，即表示空间基函数；表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；||||₂表示取2-范数；||||₁表示取L₁范数；

最优解计算单元，用于根据所述优化判据获得对所述空间基形变矩阵的估计值；及

空间基求解单元，用于利用所述空间傅里叶变化矩阵求解获得所述空间基函数。

11.根据权利要求10所述的动态磁共振图像成像系统，其特征在于，所述模型构建单元中引入一Huber函数将所述优化判据变更为下述公式(3)：

$\underset{\{U_s,G\}}{\min }{\left|\right|d-\mathrm{\&Omega;}\left(F_s{U}_s{V}_t\right)\left|\right|}_2^2+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{\left|\right|U_s-G\left|\right|}_F^2+{\left|\right|\mathrm{vec}\left(G\right)\left|\right|}_1$ 公式(3)

其中，vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；G表示辅助矩阵，所述辅助矩阵G为行为N、列为L的矩阵，其中，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；F_s是空间傅里叶变化矩阵；λ是正则化参数；U_s表示空间基形变矩阵，即U_s＝F_s ^-1U_k，U_k表示空间基函数；||||₂表示取2-范数；||||₁表示取L₁范数；α表示Huber函数的控制参数；||||_F表示取F-范数；V_t表示时间基函数；

所述最优解计算单元包括：

交替最小化求解单元，用于基于所述控制参数α的一固定值，通过多次迭代的交替最小化求解过程求解所述公式(3)中的辅助矩阵和空间基形变矩阵，直至所述辅助矩阵和空间基形变矩阵的中间量均收敛；和

控制参数调整单元，用于逐渐减小所述控制参数α的值，基于每一次调整所述控制参数α的调整值调用执行所述交替最小化求解单元，直至将所述控制参数α的值减小到零，将最后一次迭代计算获得的所述空间基形变矩阵的中间量作为所述空间基形变矩阵的估计值。

12.根据权利要求11所述的动态磁共振图像成像系统，其特征在于，所述交替最小化求解单元包括：

控制参数获取单元，用于获取控制参数α的值；

第一计算单元，用于基于所述空间基形变矩阵的初始值，按照下述公式(4)计算所述辅助矩阵G的中间量；

   公式(4)

其中，表示软阈值算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；||||_F表示取F-范数；vec(G)表示将所述辅助矩阵G的列向量拉直转为多维向量；l表示迭代次数，表示第l次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；α表示Huber函数的控制参数；

第二计算单元，用于基于所述辅助矩阵的中间量，按照下述公式(5)计算此次迭代过程中的所述空间基形变矩阵的中间量，并将该值作为下一次执行步骤22时所述空间基形变矩阵的初始值：

   公式(5)

其中，表示单位算子，N表示K空间的所有像素点数，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；α表示Huber函数的控制参数；及

第一判断单元，用于判断所述辅助矩阵的中间量和所述空间基形变矩阵的中间量是否均同时收敛；若否则将迭代次数加一，继续调用所述第一计算单元和第二计算单元；若是，则将最后一次迭代计算获得的所述空间基形变矩阵的中间量，作为基于所述控制参数α的另一调整值、调用所述交替最小化求解单元时所采用的所述空间基形变矩阵的初始值。

13.根据权利要求12所述的动态磁共振图像成像系统，其特征在于，所述第一计算单元包括：

元素值提取单元，用于提取第l次迭代时获得的所述空间基形变矩阵的中间量中各个元素的值；

第二判断单元，用于判断所述空间基形变矩阵的中间量的元素值的绝对值与所述控制参数α值的大小，当所述空间基形变矩阵的中间量的元素值的绝对值小于所述控制参数α的值时，在所述辅助矩阵中的相应位置赋值为零；当所述空间基形变矩阵的中间量的元素值的绝对值大于等于所述控制参数α的值时，在所述辅助矩阵中的相应位置赋值为下述公式(6)所表示的值：

$G_{n,l1}^{(l+1)}=\frac{Q_{n,l1}}{\left|Q_{n,l1}\right|}\left(\right|Q_{n,l1}|-\mathrm{\&alpha;})$    公式(6)

其中，n＝1,...,N，N表示K空间的所有像素点数；l1＝1,...,L，L表示动态成像重建模型中空间基函数与时间基函数的模阶；表示第l+1次迭代时获得的所述辅助矩阵G的中间量中矩阵位置在(n,l1)处的值；Q_n,l1表示第l次迭代时获得的所述空间基形变矩阵的中间量中矩阵位置在(n,l1)的值；α表示所述Huber函数中的控制参数；及

输出单元，用于按照所述赋值后的结果形成N行、L列的矩阵，作为第l+1次迭代时获得的所述辅助矩阵的中间量。

14.根据权利要求12所述的动态磁共振图像成像系统，其特征在于，所述第二计算单元包括：

优解模型构建单元，用于调用按照如下公式(7)所构建的关联和G^(l+1)的最优解方程：

$(A^{*}A+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}I)U_s^{(l+1)}=A^{*}\left(d\right)+\frac{\mathrm{\&lambda;}}{2\mathrm{\&alpha;}}{G}^{(l+1)}$    公式(7)

其中，λ是正则化参数；α表示Huber函数的控制参数；A表示求导获得的矩阵Ω(F_sU_sV_t)，A^*表示矩阵A的共轭矩阵；A^*(d)中的d表示欠采的图像数据矩阵；表示第l+1次迭代时获得的所述空间基形变矩阵的中间量；G^(l+1)表示第l+1次迭代时获得的所述辅助矩阵G的中间量；是单位算子，L表示所述动态成像重建模型中空间基函数与时间基函数的模阶；和

梯度计算单元，用于利用共轭梯度算法迭代优化求解所述公式(7)中的所述空间基形变矩阵的中间量