CN106056647A - 一种基于卷积稀疏双层迭代学习的磁共振快速成像方法 - Google Patents

Abstract

一种基于卷积稀疏双层迭代学习的磁共振快速成像方法，包括以下步骤：在双层伯格曼字典迭代学习框架上融入卷积字典训练和稀疏系数更新，建立图像稀疏表示模型；在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新卷积字典和稀疏系数，特别地利用软阈值迭代法求解稀疏系数的子问题，更新稀疏系数；在双层伯格曼外层迭代上进行频域数据更新和图像更新，得到重建图像。本发明根据模型从自然图像库中获取卷积字典，该卷积字典能够有效地表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征，因此可以在更少的测量下更精确的重建磁共振图像，减少重建图像的伪影，恢复更多的图像细节。

Description

一种基于卷积稀疏双层迭代学习的磁共振快速成像方法

技术领域

本发明属于医学成像领域，尤其涉及磁共振成像。

背景技术

磁共振成像(MRI)作为一项较新的医学影像诊断技术，虽然具有无辐射和无电离、高的软组织对比分辨力，无骨伪影干扰，不用对比剂即可进行血流成像，其多参数成像便于对照比较并可获得多方位成像等优点。但是，磁共振成像系统的主要缺点是成像速度慢。扫描时间较长不适合检查运动性器官和危重病人，尤其对于噪动或丧失自制能力的患者，如不使用镇静剂，难以成像。因此如何减少扫描时间而又得到较好的重建图像成为一个至关重要的问题。

稀疏表示与图像处理中的许多问题紧密相关，它已广泛应用于图像处理和计算机视觉的各个领域，如图像压缩、图像去噪、图像去模糊、图像修复、场景分类与模式识别等。在过去几十年中，国内外学者关于稀疏表示理论的研究主要是针对稀疏表示模型和字典设计及其应用两方面，并取得了一系列的研究成果。由于图像是一种复杂的高维数据结构，现有的稀疏表示算法还存在一些缺点和不足之处。

稀疏表示模型中的一个关键问题是过完备字典的设计。字典中原子表征了图像信息的几何结构特征，原子特征种类越丰富越能够对图像形成最优稀疏表示。目前过完备字典的获取方式有两种：利用数学工具生成字典和基于训练样本学习字典。基于数学调和分析生成的字典有：傅立叶字典、小波包字典、各向同性Gabor字典、各向异性高斯混合字典、Gabor感知多成份字典等。傅立叶字典能够有效表示平滑信号，然而现实中绝大部信号具有奇异点，是非平滑的。随着小波理论的发展，学者们提出了许多小波变换工具，例如脊波、轮廓波、曲线波等。这些小波分析工具能够有效匹配图像中的方向、边缘、平滑特征，但却不能表示纹理类复杂几何结构特征。自然图像统计特性表明，理想的表征图像特征的原子应具有以下性能：多分辨率、时频局部化、多方向性、各向异性。

为获取理想特征原子，更为精确、自适应的匹配图像特征，研究者们提出从自然图像中通过学习构建字典，主要方法有：最优方向法、K-SVD、CNDL-FOCUSS、在线学习。MOD算法学习到的字典原子特征与简单细胞感受野特性极为相似，具有类Gabor特性，但该算法在字典更新过程中涉及逆矩阵问题使得计算复杂度高。而Michal Aharon提出的K-SVD算法利用奇异值分解方式代替逆矩阵逐个更新原子，在很大程度上简化了计算复杂度。从数据样本中通过学习方式获取的字典能够有效提取到表征图像局部几何结构的特征，但是该类算法一般基于图像块训练，忽略了块之间的空间结构，不能完全有效刻画图像各个层面的结构特征。1998年，Lecun等人在文献中首次提出卷积神经网络模型，并采用监督学习方式利用反向传播算法训练模型参数，分层提取图像的初级特征、中级特征、高级特征。该模型学习到的滤波器具有平移不变性，与细胞的感受野等级特性相对应。

在卷积模型中，把图像看成是一组特征响应与一组滤波器卷积求和，并对特征响应采用l₁范数进行稀疏约束，建立如下目标函数：

$\underset{z}{min}\left|\right|x-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2,s.t.\underset{k=1}{\overset{K}{\Σ}}|\left|z_k\right|{|}_0\≤\mathrm{\ϵ}$

其中，x表示自然图像；d_k表示第k个滤波器，所有的滤波器构成滤波器组d，称之为卷积字典；z_k表示自然图像与第k个滤波器的卷积结果即特征响应；K为卷积字典大小，即滤波器个数；ε表示所有特征响应的非零系数个数，体现特征响应的稀疏度。

发明内容

本发明目的是提出一种基于卷积稀疏双层迭代学习的磁共振快速成像方法(CSCMRI)。

本发明所述的一种基于卷积稀疏双层迭代学习的磁共振快速成像方法，从磁共振图像库中训练获取卷积字典，该字典能够有效表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征。再在卷积字典张成的特征空间下，利用增广拉格朗日和交替方向法的卷积实现快速磁共振图像重建，进一步地稀疏表示图像，并且在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。

本发明通过以下技术方案实现。

本发明所述的一种基于卷积稀疏双层迭代学习的磁共振快速成像方法，包括以下步骤：

步骤(a)：在双层伯格曼字典迭代学习框架上融入卷积字典训练和稀疏系数更新，建立图像稀疏表示模型。

步骤(b)：在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新卷积字典和稀疏系数，特别地利用软阈值迭代法求解稀疏系数的子问题，更新稀疏系数。

步骤(c)：在双层伯格曼外层迭代上进行频域数据更新和图像更新，得到重建图像。

进一步地说，本发明所述步骤(a)为：在双层伯格曼字典学习迭代框架上利用卷积字典学习更新字典，建立的图像稀疏模型：

$\left\{\begin{array}{c}{u}^{i+1}=\underset{u}{\arg \min }\left\{\underset{d,z}{\min }\frac{1}{2}\right||u-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+v_1\left|\right|F_{p}u-f^i|{|}^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|z_k\right|{|}_1+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(d_k\right)\}\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right.$

其中，模型中第一项为图像进行字典学习和稀疏系数表示的正则项，第二项保证重建结果与K空间欠采样数据保持匹配约束，第三项表示用L1范数约束稀疏系数，第四项表示对字典d_k的约束，β表示在最优字典中图像块的稀疏水平，v₁表示K数据拟合的权重，F_p表示部分傅里叶变换，f表示相应的频域数据。

利用分离变量法，对上述非约束问题分离出分别与d，z，u有关的项。

进一步地说，本发明所述步骤(b)为：

1)其中分离出的与d有相关的子问题为：

$\underset{d}{min}\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(d_k\right)$

通过增加辅助变量t₁,t₂，把非约束问题转化为约束问题，则关于的_d子问题转化为：

$min\frac{1}{2}\left|\right|u-t_1|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(t_{2,k}\right)$

$\begin{array}{cc}s.t.& t_1=\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k,t_2=d\end{array}$

增加辅助变量t₁,t₂后，将约束问题变为非约束问题：

$\underset{d^{i,j+1},t_1^{j+1},t_2^{j+1}}{min}\frac{1}{2}\left|\right|u-t_1|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(t_{2,k}\right)+\frac{{\mathrm{\μ}}_1}{2}||t_1-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k^{i,j}+{\mathrm{\λ}}_1^j|{|}_2^2+\frac{{\mathrm{\μ}}_2}{2}\left|\right|t_2-d+{\mathrm{\λ}}_2^j|{|}_2^2$

步骤(b)中与d有相关的子问题为轮换地更新一个变量，同时固定其它变量固定辅助变量t₁,t₂，通过最小化二次多项式更新卷积字典d；固定卷积字典d和辅助变量t₂，通过最小化二次多项式更新辅助变量t₁；固定卷积字典d和辅助变量t₁，通过最小化二次多项式更新辅助变量t₂。

2)分离出的与z有相关的子问题为：

$\underset{z}{min}\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|z_k\right|{|}_1$

通过增加辅助变量t₃,t₄，把与z有相关的非约束问题转化为约束问题，则关于的z子问题转化为：

$\underset{z}{min}\frac{1}{2}\left|\right|u-t_3|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|t_{4,k}\right|{|}_1$

$\begin{array}{cc}s.t.& t_3=\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k,t_4=z\end{array}$

增加辅助变量t₃,t₄后，再将约束问题变为非约束问题：

$\underset{z^{i,j+1},t_3^{j+1},t_4^{j+1}}{\min }\frac{1}{2}{\left|\right|u-t_3\left|\right|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|t_{4,k}\left|\right|}_1+\frac{{\mathrm{\μ}}_3}{2}{\left|\right|t_3-\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{d}_k^{i,j+1}*z_k+{\mathrm{\λ}}_3^j\left|\right|}_2^2+\frac{{\mathrm{\μ}}_4}{2}{\left|\right|t_4-z+{\mathrm{\λ}}_4^j\left|\right|}_2^2$

步骤(b)中与z有相关的子问题为通过轮换地更新一个变量，同时固定其它变量：固定辅助变量t₃,t₄，通过最小化二次多项式更新稀疏系数z；固定稀疏系数z和辅助变量t₄，通过最小化二次多项式更新辅助变量t₃；固定稀疏系数z和辅助变量t₃，通过最小化二次多项式更新辅助变量t₄。

进一步地说，本发明所述步骤(c)与u有相关的子问题，它在双层伯格曼字典学习迭代过程的每一次内部迭代之后进行更新。

$\left\{\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+v_1||F_{p}u-f^i|{|}^2\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right.$

本发明的技术方案具有以下的优点或有益效果：本发明基于卷积稀疏编码学习的磁共振快速成像方法，按模型训练方法从自然图像库中获取卷积字典，该字典能够有效表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征。再在卷积字典张成的特征空间下，利用增广拉格朗日和交替方向法实现快速磁共振图像重建。进一步的稀疏表示图像，并且可以在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。

由于本发明的算法可以有效表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征，因此进一步地使得重构的图像达到令人满意的效果，因此能够快速、精确的重构磁共振图像。

附图说明

图1为本发明算法步骤的流程图。

图2为模拟径向采样轨迹下CSCMRI算法的重建图像。(a)为模拟径向采样轨迹；(b)为训练得到的卷积字典；(c)为重建图像。

图3为TV、DLMRI和CSCMRI三种算法重建图像的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化情况。

图4为TV、DLMRI和CSCMRI三种算法重建图像的高频误差(HFEN)随欠采样因子(Downsampling Factor)的变化情况。

图5为10倍欠采样率的模拟径向采样轨迹下TV、DLMRI和CSCMRI三种算法的重建性能分析情况。其中(a)、(b)、(c)分别为TV、DLMRI和CSCMRI三种算法10倍欠采样率下重建图像；(d)、(e)、(f)分别为TV、DLMRI和CSCMRI三种算法10倍欠采样率下重建图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。此处所描述的具体实施例仅用于解释本发明技术方案，并不限于本发明。

参见本发明实施例的附图，下文将更详细地描述本发明。

现参考附图1描述根据本发明的本发明实施例基于卷积稀疏编码学习的磁共振快速成像方法。根据本发明的方法，本发明技术方案在双层伯格曼迭代框架卷积字典学习算法上，利用增广拉格朗日和交替方向法实现快速磁共振图像重建，得到更好的图像重建效果。本发明技术方案算法的重建结果可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述。

步骤(a)：在双层伯格曼字典学习框架上融入图像域u的卷积稀疏编码，建立图像稀疏表示模型：

$\begin{array}{c}\underset{u,d,z}{\min }\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+v_1||F_{p}u-f^i|{|}^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}\left|\right|z_k|{|}_1\\ \begin{array}{cc}s.t.\left|\right|d_k|{|}_2^2\≤1& \∀k\∈\{1,\mathrm{...},k\}\end{array}\end{array}---\left(1\right)$

其中，u表示待重建磁共振图像；d_k表示第k个滤波器，所有的滤波器构成滤波器组d，称之为卷积字典；z_k表示自然图像与第k个滤波器的卷积结果即特征响应；K为卷积字典大小，即滤波器个数；表示所有特征响应的非零系数个数，体现特征响应的稀疏度。

将式(1)约束问题变为非约束问题:

$\underset{u,d,z}{min}\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+v_1||F_{p}u-f^i|{|}^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}\left|\right|z_k|{|}_1+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(d_k\right)---\left(2\right)$

通过分离变量法，分别分离出与d,z,u有关的项。

步骤(b)：对于公式(2)，为求解子问题d，首先，通过分离变量法分离出与d有关的项，在双层伯格曼内层迭代上利用增加辅助变量和轮换技术更新学习卷积滤波器字典d：

$\underset{d}{min}\frac{1}{2}\left|\right|u^i-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(d_k\right)---\left(3\right)$

用增广拉格朗日方法解决上式(3)，增加辅助变量t₁,t₂，使得子问题更易求解，得:

$\begin{array}{c}\{d^{i,j+1},t_1^{j+1},t_2^{j+1}\}=\arg \min \frac{1}{2}\left|\right|u^i-t_1|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(t_{2,k}\right)\\ \begin{array}{cc}s.t.& t_1=\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k^{i,j},t_2=d\end{array}\end{array}---\left(4\right)$

将约束问题变为非约束问题：

$\underset{d^{i,j+1},t_1^{j+1},t_2^{j+1}}{min}\frac{1}{2}\left|\right|u^i-t_1|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(t_{2,k}\right)+\frac{{\mathrm{\μ}}_1}{2}||t_1-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k^{i,j}+{\mathrm{\λ}}_1^j|{|}_2^2+\frac{{\mathrm{\μ}}_2}{2}\left|\right|t_2-d+{\mathrm{\λ}}_2^j|{|}_2^2---\left(5\right)$

通过轮换地更新一个变量求解等式(5)，同时固定其它变量：固定辅助变量t₁,t₂，通过最小化二次多项式更新卷积字典d_k；固定卷积字典d和辅助变量t₂，通过最小化二次多项式更新辅助变量t₁；固定卷积字典d和辅助变量t₂，通过最小化二次多项式更新辅助变量t₂。

对于公式(5)中的卷积字典d，利用分离变量法，分离出式(5)中与d有关的项:

$d^{i,j+1}=\arg \underset{d}{min}\frac{{\mathrm{\μ}}_1}{2}\left|\right|t_1^j-\underset{k=1}{\overset{K}{\Σ}}{d}_k*z_k^{i,j}+{\mathrm{\λ}}_1^j|{|}_2^2+\frac{{\mathrm{\μ}}_2}{2}||t_2^j-d+{\mathrm{\λ}}_2^j|{|}_2^2---\left(6\right)$

利用帕赛瓦定理求解，将公式(6)转换到频域中，得到

${\hat{d}}^{i,j+1}=\arg \underset{d}{min}\frac{{\mathrm{\μ}}_1}{2}\left|\right|{\hat{t}}_1^j-\underset{k=1}{\overset{K}{\Σ}}{\hat{z}}_k^{i,j}{\hat{d}}_k+{\widehat{\mathrm{\λ}}}_1^j|{|}_2^2+\frac{{\mathrm{\μ}}_2}{2}||{\hat{t}}_2^j-\hat{d}+{\widehat{\mathrm{\λ}}}_2^j|{|}_2^2---\left(7\right)$

通过最小化二次多项式更新卷积字典d：

$-{\mathrm{\μ}}_1{\left({\hat{z}}_k^{i,j}\right)}^T({\hat{t}}_1^j-\underset{k=1}{\overset{K}{\Σ}}{\hat{z}}_k^{i,j}{\hat{d}}_k+{\widehat{\mathrm{\λ}}}_1^j)-{\mathrm{\μ}}_2({\hat{t}}_2^j-\hat{d}+{\widehat{\mathrm{\λ}}}_2^j)=0---\left(8\right)$

${\hat{d}}^{i,j+1}={({\left({\hat{z}}_k^{i,j}\right)}^T{\hat{z}}_k^{i,j}+\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}I)}^{-1}\[{\left({\hat{z}}_k^{i,j}\right)}^T({\hat{t}}_1^j+{\widehat{\mathrm{\λ}}}_1^j)+\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}({\hat{t}}_2^j+{\widehat{\mathrm{\λ}}}_2^j)\]---\left(9\right)$

再对进行傅里叶反变换，将变换回到时域：

$d^{i,j+1}=F^{-1}\{{\left({\left({\hat{z}}_k^{i,j}\right)}^T{\hat{z}}_k^{i,j}+\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}I\right)}^{-1}\[{\left({\hat{z}}_k^{i,j}\right)}^T({\hat{t}}_1^j+{\widehat{\mathrm{\λ}}}_1^j)+\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}({\hat{t}}_2^j+{\widehat{\mathrm{\λ}}}_2^j)\]\}---\left(10\right)$

对于公式(5)中的辅助变量t₁，利用分离变量法，分离出其中与t₁有关的项:

$\underset{t_1}{min}\frac{1}{2}\left|\right|u^i-t_1|{|}_2^2+\frac{{\mathrm{\μ}}_1}{2}||t_1-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j}+{\mathrm{\λ}}_1^j|{|}_2^2---\left(11\right)$

求解公式(11)最小二乘问题，求得：

$t_1^{j+1}=\frac{1}{{\mathrm{\μ}}_1+1}\[{\mathrm{\μ}}_1(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j}-{\mathrm{\λ}}_1^j)+u^i\]---\left(12\right)$

对于公式(5)中的辅助变量t₂，利用分离变量法，分离出式(5)中与t₂有关的项:

$t_2^{j+1}=\underset{t_2}{\arg \min }\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(t_{2,k}\right)+\frac{{\mathrm{\μ}}_2}{2}\left|\right|t_2-d^{i,j+1}+{\mathrm{\λ}}_2^j|{|}_2^2---\left(13\right)$

求解公式(13)最小二乘问题，求得

$t_2^{j+1}=\left\{\begin{array}{cc}\frac{d_k^{i,j+1}-{\mathrm{\λ}}_2^j}{\left|\right|d_k^{i,j+1}-{\mathrm{\λ}}_2^j|{|}_2},& \left|\right|d_k^{i,j+1}-{\mathrm{\λ}}_2^j|{|}_2^2\≥1\\ d_k^{i,j+1}-{\mathrm{\λ}}_2^j,& else\end{array}\right.---\left(14\right)$

对于公式(2)，在双层伯格曼内层迭代上利用增加辅助变量和轮换技术更新稀疏系数z。为求解子问题z，首先，通过分离变量法分离出与z有关的项：

$\underset{z}{\min }\frac{1}{2}\left|\right|u^i-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|z_k\right|{|}_1---\left(15\right)$

使用增广拉格朗日方法解决上式(15)，增加辅助变量t₃,t₄，使得求解稀疏系数z的子问题更易求解:

$\begin{array}{c}\underset{z}{\min }\frac{1}{2}\left|\right|u^i-t_3|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|t_{4,k}\right|{|}_1\\ \begin{array}{cc}s.t.& t_3=\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k,t_4=z\end{array}\end{array}---\left(16\right)$

将约束问题变为非约束问题：

$\underset{z^{i,j+1},t_3^{j+1},t_4^{j+1}}{min}\frac{1}{2}\left|\right|u^i-t_3|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|t_{4,k}\right|{|}_1+\frac{{\mathrm{\μ}}_3}{2}\left|\right|t_3-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k+{\mathrm{\λ}}_3^j|{|}_2^2+\frac{{\mathrm{\μ}}_4}{2}||t_4-z+{\mathrm{\λ}}_4^j|{|}_2^2---\left(17\right)$

通过轮换地更新一个变量求解等式(17)，同时固定其它变量：固定辅助变量t₃,t₄，通过最小化二次多项式更新稀疏系数z；固定稀疏系数z和辅助变量t₄，通过最小化二次多项式更新辅助变量t₃；固定稀疏系数z和辅助变量t₃，通过最小化二次多项式更新辅助变量t₄。

对于公式(17)中的稀疏系数z，利用分离变量法，分离出式(17)中与z有关的项:

$z_k^{i,j+1}=\arg \underset{z}{min}\frac{{\mathrm{\μ}}_3}{2}\left|\right|t_3^j-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j}+{\mathrm{\λ}}_3^j|{|}_2^2+\frac{{\mathrm{\μ}}_4}{2}||t_4^j-z+{\mathrm{\λ}}_4^j|{|}_2^2---\left(18\right)$

利用帕赛瓦定理求解，将公式(18)转换到频域中，得到

${\hat{z}}_k^{i,j+1}=\arg \underset{z}{min}\frac{{\mathrm{\μ}}_3}{2}\left|\right|{\hat{t}}_3^j-\underset{k=1}{\overset{K}{\Σ}}{\hat{z}}_k^{i,j}{\hat{d}}_k^{i,j+1}+{\widehat{\mathrm{\λ}}}_3^j|{|}_2^2+\frac{{\mathrm{\μ}}_4}{2}||{\hat{t}}_4^j-\hat{z}+{\widehat{\mathrm{\λ}}}_4^j|{|}_2^2---\left(19\right)$

通过最小化二次多项式更新稀疏系数z：

$-{\mathrm{\μ}}_3{\left({\hat{d}}_k^{i,j+1}\right)}^T({\hat{t}}_3^j+{\widehat{\mathrm{\λ}}}_3^j)+{\mathrm{\μ}}_3{\left({\hat{d}}_k^{i,j+1}\right)}^T\underset{k=1}{\overset{K}{\Σ}}{\hat{z}}_k^{i,j}{\hat{d}}_k^{i,j+1}+{\mathrm{\μ}}_4\hat{z}-{\mathrm{\μ}}_4({\hat{t}}_4^j+{\widehat{\mathrm{\λ}}}_4^j)=0---\left(20\right)$

${\hat{z}}_k^{i,j+1}={({\left({\hat{d}}_k^{i,j+1}\right)}^T{\hat{d}}_k^{i,j+1}+\frac{{\mathrm{\μ}}_4}{{\mathrm{\μ}}_3}I)}^{-1}\[{\left({\hat{d}}_k^{i,j+1}\right)}^T({\hat{t}}_3^j+{\widehat{\mathrm{\λ}}}_3^j)+\frac{{\mathrm{\μ}}_4}{{\mathrm{\μ}}_3}({\hat{t}}_4^j+{\widehat{\mathrm{\λ}}}_4^j)\]---\left(21\right)$

再对进行傅里叶反变换，将变换回到时域：

$z_k^{i,j+1}=F^{-1}\{{({\left({\hat{d}}_k^{i,j+1}\right)}^T{\hat{d}}_k^{i,j+1}+\frac{{\mathrm{\μ}}_4}{{\mathrm{\μ}}_3}I)}^{-1}\[{\left({\hat{d}}_k^{i,j+1}\right)}^T({\hat{t}}_3^j+{\widehat{\mathrm{\λ}}}_3^j)+\frac{{\mathrm{\μ}}_4}{{\mathrm{\μ}}_3}({\hat{t}}_4^j+{\widehat{\mathrm{\λ}}}_4^j)\]\}---\left(22\right)$

对于公式(18)中的辅助变量t₃，利用分离变量法，分离出式(18)中与t₃有关的项:

$\underset{t_3^{j+1}}{min}\frac{1}{2}\left|\right|u^i-t_3|{|}_2^2+\frac{{\mathrm{\μ}}_3}{2}||t_3-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}+{\mathrm{\λ}}_3^j|{|}_2^2---\left(23\right)$

求解公式(23)最小二乘问题，求得：

$t_3^{j+1}=\frac{1}{1+{\mathrm{\μ}}_3}\left({\mathrm{\μ}}_3\right(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j}-{\mathrm{\λ}}_3^j)+u^i)---\left(24\right)$

对于公式(18)中的辅助变量t₄，利用分离变量法，分离出式(18)中与t₄有关的项:

$t_4^{j+1}=\underset{t_4}{\arg \min }\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}\left|\right|t_{4,k}|{|}_1+\frac{{\mathrm{\μ}}_4}{2}||t_4-z_k^{i,j+1}+{\mathrm{\λ}}_4^j|{|}_2^2---\left(25\right)$

使用软阈值迭代算法(ISTA)求解公式(25)，得:

同时更新拉格朗日乘子λ₁,λ₂,λ₃,λ₄

${\mathrm{\λ}}_1^{j+1}={\mathrm{\λ}}_1^j-(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}-t_1^{j+1}),{\mathrm{\λ}}_2^{j+1}={\mathrm{\λ}}_2^j-(d^{i,j+1}-t_2^{j+1})---\left(27\right)$

${\mathrm{\λ}}_3^{j+1}={\mathrm{\λ}}_3^j-(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}-t_3^{j+1}),{\mathrm{\λ}}_4^{j+1}={\mathrm{\λ}}_4^j-(z^{i,j+1}-t_4^{j+1})---\left(28\right)$

步骤(c)：外层迭代中：求子问题u，与u有关的子问题

$\left\{\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\frac{1}{2}\left|\right|u-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}|{|}_2^2+v_1||F_{p}u-f^i|{|}^2\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right.---\left(29\right)$

上式为双层伯格曼迭代的外层迭代，通过求解最小二乘解析问题，并且左右两边同时反傅里叶变换，得到更新的图像uⁱ⁺¹：

$u^{i+1}=F^{-1}\left(\frac{F\[t_5+2v_1{F}_p^T{f}^i\]}{(F+2v_1{\mathrm{FF}}_p^T{F}^T{\mathrm{FF}}_p{F}^T)}\right)---\left(30\right)$

综上所述在双层伯格曼字典学习迭代框架上通过卷积稀疏编码学习的磁共振快速成像方法，得到最终成像结果。

具体而言，本发明实施例基于卷积稀疏编码学习的磁共振快速成像方法，按模型训练方法从自然图像库中获取卷积字典，该字典能够有效表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征。再在卷积字典张成的特征空间下，利用增广拉格朗日和交替方向法的卷积稀疏逼近实现快速磁共振图像重建。进一步的稀疏表示图像，并且可以在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。

综上所述，本发明的实施例完整的CSCMRI算法可归纳如下：

(1)初始化：

$d^{\left(0\right)}=0,z^{\left(0\right)}=0,{\mathrm{\λ}}_k^{\left(0\right)}=0(k=1,2,3,4),u^{\left(0\right)}=F_p^{T}f,f^0=f$

(2)内层迭代更新卷积字典d：

$\underset{d}{min}\frac{1}{2}\left|\right|u^i-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j}*z_k^{i,j}|{|}_2^2+\underset{k=1}{\overset{K}{\Σ}}{\mathrm{ind}}_C\left(d_k^{i,j}\right)$

(3)内层迭代更新稀疏系数z：

$\underset{z}{min}\frac{1}{2}\left|\right|u^i-\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j}|{|}_2^2+\mathrm{\β}\underset{k=1}{\overset{K}{\Σ}}|\left|z_k^{i,j}\right|{|}_1$

(4)更新拉格朗日乘子：

${\mathrm{\λ}}_1^{j+1}={\mathrm{\λ}}_1^j-(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}-t_1^{j+1});{\mathrm{\λ}}_2^{j+1}={\mathrm{\λ}}_2^j-(d^{i,j+1}-t_2^{j+1})$

${\mathrm{\λ}}_3^{j+1}={\mathrm{\λ}}_3^j-(\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}-t_3^{j+1});{\mathrm{\λ}}_4^{j+1}={\mathrm{\λ}}_4^j-(z^{i,j+1}-t_4^{j+1})$

(5)外层迭代更新图像u：

$u^{i+1}=F^{-1}\left(\frac{F\[\underset{k=1}{\overset{K}{\Σ}}{d}_k^{i,j+1}*z_k^{i,j+1}+2v_1{F}_p^T{f}^i\]}{(F+2v_1{\mathrm{FF}}_p^T{F}^T{\mathrm{FF}}_p{F}^T)}\right)$

(6)更新f：

fⁱ⁺¹＝fⁱ+f-F_puⁱ⁺¹

(7)结束。

本发明技术方案采用各种不同的欠采样因子对所提出方法的性能进行评估。本发明提出的CSCMRI方法与DLMRI和双层伯格曼字典学习方法TV相比，DLMRI直接通过正交匹配追踪方法解决l₀最小化问题，并且是基于块的字典稀疏表示方法，CSCMRI方法则从磁共振图像库中训练获取卷积字典，该卷积字典能够有效表示图像中的许多几何结构特征。实验过程各种参数的标准值分别设置如下：卷积核大小为11×11，字典的过完备性(相当于J＝36)，字典块的步长r＝1，因为所以总采样数据为L＝262144。参数设置为DLMRI的默认值。重建图像的质量通过使用峰值信噪比(PSNR)和高频误差(HFEN)来衡量。

图1为示出本发明算法步骤的流程图；

图2为模拟径向采样轨迹下CSCMRI算法的重建图像。(a)为模拟径向采样轨迹；(b)为训练得到的卷积字典；(c)为重建图像。

图3为下TV，DLMRI和CSCMRI三种算法重建图像的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化情况。

图4为下TV，DLMRI和CSCMRI三种算法重建图像的高频误差(HFEN)随欠采样因子(Downsampling Factor)的变化情况。

图5为10倍欠采样率的模拟径向采样轨迹下TV、DLMRI和CSCMRI三种算法的重建性能分析情况。其中，(a)、(b)、(c)分别为TV、DLMRI和CSCMRI三种算法10倍欠采样率下重建图像；(d)、(e)、(f)分别为TV、DLMRI和CSCMRI三种算法10倍欠采样率下重建图像。

本发明实施例基于卷积稀疏编码学习的磁共振快速成像方法，按模型训练方法从自然图像库中获取卷积字典，该字典能够有效表示图像中的许多几何结构特征，如平滑、方向、边缘、斑块、纹理等特征。再在卷积字典张成的特征空间下，利用增广拉格朗日和交替方向法的卷积稀疏逼近实现快速磁共振图像重建。进一步的稀疏表示图像，并且可以在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。

因本技术领域的技术人员应理解，本发明可以以许多其他具体形式实现而不脱离本发明的精神或范围。尽管业已描述了本发明的实施例，应理解本发明不应限制为这些实施例，本技术领域的技术人员可如所附权利要求书界定的本发明精神和范围之内作出变化和修改。

Claims (4)

1.一种基于卷积稀疏双层迭代学习的磁共振快速成像方法，其特征是包括以下步骤：

步骤(a)：在双层伯格曼字典迭代学习框架上融入卷积字典训练和稀疏系数更新，建立图像稀疏表示模型；

步骤(b)：在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新卷积字典和稀疏系数，特别地利用软阈值迭代法求解稀疏系数的子问题，更新稀疏系数；

步骤(c)：在双层伯格曼外层迭代上进行频域数据更新和图像更新，得到重建图像。

2.根据权利要求1所述的基于卷积稀疏双层迭代学习的磁共振快速成像方法，其特征是所述步骤(a)为：

在双层伯格曼字典学习迭代框架上利用卷积字典学习更新字典，建立的图像稀疏模型：

其中，模型中第一项为图像进行字典学习和稀疏系数表示的正则项，第二项保证重建结果与K空间欠采样数据保持匹配约束，第三项表示用L1范数约束稀疏系数，第四项表示对字典d_k的约束，β表示在最优字典中图像块的稀疏水平，v₁表示K数据拟合的权重，F_p表示部分傅里叶变换，f表示相应的频域数据；

利用分离变量法，对上述非约束问题分离出分别与d，z，u有关的项。

3.根据权利要求1所述的基于卷积稀疏双层迭代学习的磁共振快速成像方法，其特征是所述步骤(b)为：

1)其中分离出的与d有相关的子问题为：

通过增加辅助变量t₁,t₂，把非约束问题转化为约束问题，则关于的d子问题转化为：

增加辅助变量t₁,t₂后，将约束问题变为非约束问题：

与d有相关的子问题为轮换地更新一个变量，同时固定其它变量固定辅助变量t₁,t₂，通过最小化二次多项式更新卷积字典d；固定卷积字典d和辅助变量t₂，通过最 小化二次多项式更新辅助变量t₁；固定卷积字典d和辅助变量t₁，通过最小化二次多项式更新辅助变量t₂；

2)分离出的与z有相关的子问题为：

通过增加辅助变量t₃,t₄，把与z有相关的非约束问题转化为约束问题，则关于的z子问题转化为：

增加辅助变量t₃,t₄后，再将约束问题变为非约束问题：

与z有相关的子问题为通过轮换地更新一个变量，同时固定其它变量：固定辅助变量t₃,t₄，通过最小化二次多项式更新稀疏系数z；固定稀疏系数z和辅助变量t₄，通过最小化二次多项式更新辅助变量t₃；固定稀疏系数z和辅助变量t₃，通过最小化二次多项式更新辅助变量t₄。

4.根据权利要求1所述的基于卷积稀疏双层迭代学习的磁共振快速成像方法，其特征是所述步骤(c)与u有相关的子问题，它在双层伯格曼字典学习迭代过程的每一次内部迭代之后进行更新：

。