CN104574456A - 一种基于图正则化稀疏编码的磁共振超欠采样k数据成像方法 - Google Patents

Abstract

一种基于图正则化稀疏编码的磁共振超欠采样K数据成像方法，包括以下步骤：(a)：在双层伯格曼迭代框架上进行图正则化稀疏编码表示，得到图像稀疏模型；(b)：引入辅助变量和轮换求解的技术，在双层伯格曼迭代的内层迭代上更新学习字典和稀疏系数；(c)：利用部分超欠采样的K数据约束，在双层伯格曼迭代的外层迭代上进行图像更新，以得到成像结果。本发明通过自适应字典学习引入图正则化稀疏编码方法，建立邻近图来编码局部结构数据以及挖掘其在几何数据方面的约束，使得图像数据可以更好的稀疏表示；另外本发明可以处理局部几何特征更复杂的图像，能有效的捕获局部的图像结构，恢复出更多的图像细节，得到的图像结果具有更好的保真度。

Description

一种基于图正则化稀疏编码的磁共振超欠采样K数据成像方法

技术领域

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

背景技术

磁共振成像是一种重要的医学诊断手段，尤其在缺乏离子化的情况下为临床医生提供了重要的解剖结构。尽管，磁共振成像使高分辨率图像很好的区分表示软组织，但是其成像速度受到物理和生理上的限制。成像速度慢是磁共振成像***的主要缺点，使得磁共振成像检查的适应症大为减少，再者它不适合于运动性器官的检查和危重病人的检查，增加扫描持续时间可能会带来一些生理运动伪影等。因此自磁共振成像出现以来，人们一直致力于成像速度和成像质量的提高。

磁共振成像速度慢与扫描时间有关。扫描时间又与采样率成正比，减少扫描时间的同时会使采样率相应下降，重建的图像分辨率也会降低。为了保证图像的质量，需要的增加图像的先验知识，这类方法被称为正则化方法。近年来发展起来的以稀疏表示为特征的压缩感知理论是减少扫描量而保证高效重建的一种有效途径。

图像的稀疏表示具有坚实的生物背景，它最早源于“有效编码假说”。若信号在某个基越稀疏，那么所需要的采样量则越少。因此，压缩感知理论中的一个重要的问题就是稀疏基的选取。稀疏表示的好处就在于，非零系数揭示了信号与图像的内在结构与本质属性，同时非零系数具有显式的物理意义。

传统的压缩感知磁共振成像通常采用预定义字典，其可能不能很好的稀疏表达重建图像。全变差更倾向于卡通式的图像，这种图像是分段常数。因此图像处理、信息传输和计算机视觉等领域一直在寻求信号与图像的稀疏而简洁的表示方式。在图像处理应用中，小波变换、傅里叶变换、离散余弦变换和变差分等，都是典型的固定模式的稀疏变换算子，但是这些算子并不能充分利用处理对象的特征，他们可以看作是固定的字典。

鉴于变差分模型有很好的保存图像边缘的能力，其在磁共振成像邻域也得到广泛的应用。但是在欠采样大的情况下可能会产生块状效应，将该模型直接用于磁共振成像领域，图像的重质量会受到影响，Yang等人提出在TV模型上添加其它的稀疏约束项来提高重建图像质量。其模型如下：

$\underset{u}{\min }\{{\mathrm{\μ}}_1{\left|\right|u\left|\right|}_{\mathrm{TV}}+{\mathrm{\μ}}_2{\left|\right|\mathrm{\ψu}\left|\right|}_1+\frac{1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}$

式中，ψ表示稀疏约束项，μ₁,μ₂＞0，用于权衡两个前两个正则项和第三个保真项。

目前除了基于已知固定的稀疏变换的正则化来提高稀疏性，基于图像块的字典学习的图像稀疏方法也得到广泛的研究。很多研究表明，基于字典学习的稀疏编码比固定字典的稀疏编码更具有优势，在自适应的动态字典下，信号能够得到更为稀疏的表达，在很多应用中得到了非常好的效果。鉴于稀疏编码模型在图像处理各方面，特别是在图像恢复中的良好效果，设计鲁棒的数值算法是字典学习领域的一个极为重要的问题。该理论的特点是设计最优的自适应字典，使得重叠的图像块在这些字典下是更优的稀疏表示。

当前典型的字典学习方法为一种基于奇异值分解的算法K-SVD，在训练过程中迭代多次，每次计算一次SVD分解。K-SVD方法训练字典分为稀疏编码和字典更新两部分，稀疏编码时固定字典D，用匹配追踪(MP)或者正交匹配追踪(OMP)等算法迭代求解信号在字典上的稀疏系数，然后根据求得的稀疏系数再用SVD分解更新字典中的每一列(即字典的每一个原子)，如此反复得到最优化的字典。但是K-SVD的缺陷就是易于陷入局部解，即算法的有效性严重依赖于训练样本集的类型或初始字典。

Avishankar等人提出了一种两步迭代方法，将字典学习模型用到K空间欠采样的磁共振图像重建上，提出字典学习重建磁共振图像(DLMRI)模型：

$\begin{array}{c}\underset{u,D,\mathrm{\Γ}}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|{\mathrm{D\α}}_l-R_{l}u\left|\right|}_2^2+v{\left|\right|F_{p}u-f\left|\right|}_2^2\}\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,\∀l\end{array}$

此处，Γ＝[α₁,α₂,...,α_L]为所有图像块对应的稀疏稀疏矩阵。前一项用于图像块在自适应学习字典上的稀疏表示，后一项是在图像数据的拟合保真项。求解该模型的两步迭代方法，第一步是自适应字典学习；第二步是从高度欠采样的K-空间数据重建图像。虽然这些数据的学习方法比以前预定义为基础字典的方法有很大的改善，但大多数方法没有考虑K空间数据几何结构的信息，会导致图像细节的丢失。

刘且根等人提出以双层伯格曼迭代算法为主要框架的字典学习模型，例如最近提出基于双层伯格曼字典学习(TBMDU)，外层迭代与数据保真度这一项有关，内层是与字典和基于图像块稀疏有关。改进的稀疏编码和字典更新应用在内层伯格曼迭代，这样做可以使整个算法在较少次迭代后就能达到收敛。其算法模型如下：

$\left\{\begin{array}{c}{u}^{k+1}=\arg \underset{u}{\min }\{\underset{D,\mathrm{\Γ}}{\min }\underset{i}{\mathrm{\Σ}}({\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|{\mathrm{D\α}}_i-R_{i}u\left|\right|}_2^2)+\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2\}\\ f^{k+1}=f^k+f-F_p{u}^{k+1}\end{array}\right.$

式中D＝[d₁,d₂,...,d_J]∈C^M×J，Γ＝[α₁,α₂,...,α_L]∈C^J×L。λ代表在最优字典中图像块的稀疏水平。对于医学图像，λ的值可以凭经验来得到。J＝K·M，其中K测量字典过完备的程度。双层伯格曼字典学习克服了对初值敏感和计算量大的两个难题。但希望在此基础上仍能成倍地提高扫描速度，并实现K数据欠采样的高质量成像。

现有技术的缺陷在于，变差分正则化及小波稀疏约束等固定基变换，并不能理想的稀疏表达所有的图像；Avishankar和刘且根等人提出的字典学习非固定基学习方法，仅仅考虑图像的稀疏表示，没有考虑图像块空间数据所蕴含的几何结构信息，这可能会导致细节的丢失。因此，业界需要一种能够快速，精细的重构磁共振图像几何结构特征和图像细节的算法。

发明内容

本发明的目的是提出一种基于图正则化稀疏编码的磁共振超欠采样K数据成像方法(GSCMRI)。

本发明要解决的技术问题是通过图正则化稀疏编码方法，利用其在几何特征数据约束方面的优势，建立邻近图来编码局部结构数据，可以更好的描述结构信息和实现磁共振图像块的重建识别，解决现有的磁共振成像方法不能理想地稀疏表示图像结构的缺陷以及快速精确地重建磁共振图像。

本发明是通过以下技术方案步骤实现：

步骤(a)：在双层伯格曼迭代框架上进行图正则化稀疏编码表示，建立图像稀疏模型。

步骤(b)：引入辅助变量和轮换求解的技术，在双层伯格曼迭代的内层迭代上更新学习字典和稀疏系数。

步骤(c)：利用部分超欠采样的K数据约束，在双层伯格曼迭代的外层迭代上进行图像更新，以得到成像结果。

进一步地说，本发明所述步骤(a)为：在双层伯格曼迭代框架上进行图正则化稀疏编码表示，即对图像的系数进行稀疏表示的基础上，将图像数据的局部几何结构约束通过对系数图正则化以得到更好的稀疏表示，得到图像稀疏模型。

进一步地说，本发明所述步骤(b)为：引入辅助变量将步骤(a)中的无约束问题转化为约束问题求解，通过伯格曼迭代算法完成字典更新以及通过软阈值迭代更新稀疏系数。

进一步地说，本发明所述步骤(c)为：在步骤(b)完成字典学习和稀疏系数后，通过伯格曼外层迭代，求解最小二乘的解析问题，进行图像更新，得到图像重建结果。

进一步地说，所述步骤(a)为：在双层伯格曼迭代框架的稀疏模型上融入图正则化稀疏编码过程，建立的图像稀疏模型为：

$\left\{\begin{array}{c}{u}^{k+1}=\arg \underset{u}{\min }\{\underset{D,\mathrm{\Γ}}{\min }\underset{i}{\mathrm{\Σ}}({\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|{\mathrm{D\α}}_i-R_{i}u\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right))+\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2\}\\ f^{k+1}=f^k+f-F_p{u}^{k+1}\end{array}\right.$

其中，模型中第一个子问题的前三项为图像在字典上进行稀疏和图正则化表示，u表示磁共振图像，D表示学习字典，α_l表示图像块稀疏系数，R_l表示图像块选取矩阵，Γ表示图正则化稀疏表达系数，F_p表示部分傅立叶变换，f表示K数据；第四项保证重建结果与K空间欠采样数据保持匹配约束，λ，η和μ分别表示学习字典稀疏系数，图正则化编码系数和K数据拟合的权重。

进一步地说，本发明步骤(b)引入辅助变量，将图像模型的第一个子问题转化为：

$\begin{array}{c}\underset{D,\mathrm{\Γ}}{\min }\frac{\mathrm{\λ}}{2}{\left|\right|Z\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\underset{i}{\overset{M}{\mathrm{\Σ}}}{\left|\right|{\mathrm{\α}}_i\left|\right|}_1\\ s.t.Z={\mathrm{D\α}}_i-R_{i}u,{\left|\right|d_j\left|\right|}^2\≤1,j=\mathrm{1,2},...,J\end{array}$

图正则化约束在流形假设下可以有效地利用几何数据。其假设两个数据点x_i,x_j在数据分布的内在几何结构上相近，它们的表示系数α_i和α_j在新字典中也相近。给定的数据X，最近邻图G是由M个顶点构造的，每个顶点代表X中的一个数据点，W表示G的权重矩阵。如果x_i是在x_j的K接近邻域或者x_j是在x_i的K接近邻域，W_ij＝1；否则W_ij＝0。此外，x_i的程度得定义为和C＝diag(c₁,c₂,...,c_M)，L＝C-W。映射权重图G来稀疏表达系数Γ，选择合理映射的准则最小化下列目标函数：

$\frac{1}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{({\mathrm{\α}}_i-{\mathrm{\α}}_j)}^2{W}_{\mathrm{ij}}=\mathrm{Tr}\left({\mathrm{\ΓL\Γ}}^T\right)$

引入辅助变量后，采用双层伯格曼技术，将模型变更为：

$\begin{array}{c}(D^{k+1},{\mathrm{\α}}_i^{k+1},z_i^{k+1})=\arg \underset{D,{\mathrm{\α}}_i,z_i}{\min }{\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|z_i\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\frac{\mathrm{\β}}{2}{\left|\right|z_i+{\mathrm{D\α}}_i-R_{i}u-\frac{y_i^k}{\mathrm{\β}}\left|\right|}_2^2\\ Y^{k+1}=Y^k+\mathrm{\β}(\mathrm{Ru}-Z^{k+1}-D^{k+1}{\mathrm{\Γ}}^{k+1})\end{array}$

进一步地说，本发明步骤(b)所述的轮换求解：轮换地更新一个变量，同时固定其它变量：固定辅助变量z_i和系数α_i，通过最速下降法更新字典D^k+1，经奇异值分解逐列更新字典的每一列，最小化近似误差；固定辅助字典D和系数α_i，通过最小化二次多项式更新固定辅助变量z_i和字典D^k，通过利用软阈值迭代算法更新

进一步地说，本发明步骤(c)所述的图像更新是：对于磁共振图像u的子问题，它在双层伯格曼迭代过程的每一次内层迭代之后进行更新。即在外层的伯格曼迭代中，磁共振图像u通过求解模型最小化得到：

$u^{k+1}=\arg \underset{u}{\min }\underset{i}{\mathrm{\Σ}}\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2+\frac{\mathrm{\β}}{2}{\left|\right|D^{k+1}{\mathrm{\α}}_i^{k+1}-R_{i}u+z_i^{k+1}-\frac{y_i^{k+1}}{\mathrm{\β}}\left|\right|}_2^2$

进一步地说，本发明步骤(c)所述的图像更新是：使用所述的图正则化稀疏编码表示重构未采样的K空间数据，得到最终成像结果。

本发明的技术方案具有以下的优点或有益效果：本发明通过自适应字典学习引入图正则化稀疏编码方法，建立邻近图来编码局部结构数据以及挖掘其在几何数据方面的约束，使得图像数据可以更好的稀疏表示；另外本发明可以处理局部几何特征更复杂的图像，能有效的捕获局部的图像结构，恢复出更多的图像细节，得到的图像结果具有更好的保真度。

由于本技术方案能够有效的捕获局部的图像结构，达到快速、精确的重构磁共振图像。因此本发明的基于图正则化稀疏编码的磁共振超欠采样K数据成像方法使得重构的图像达到令人满意的效果。

附图说明

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

图2为模拟径向欠采样轨迹时DLMRI，TBMDU和GSCMRI三种算法下的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化。

图3为模拟径向欠采样轨迹时DLMRI、TBMDU和GSCMRI三种算法下的高频误差(HFEN)随着模拟径向欠采样轨迹欠采样因子(Downsampling Factor)的变化。

图4为模拟径向欠采样轨迹下的DLMRI，TBMDU和GSCMRI三种算法的重建效果的对比结果。(a)为原图，(b)(c)(d)和(e)(f)(g)分别是DLMRI、TBMDU和GSCMRI三种方法在8倍模拟径向轨迹欠采样率下的重建结果图和重建误差图。

图5为7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种算法在二维随机欠采样轨迹，笛卡尔采样轨迹和相位编码方向上随机欠采样K空间数据的中间部分三种采样欠采样轨迹下的峰值信噪比(PSNR)随迭代次数(Interation Number)变化的曲线图。图5中9根曲线依次分别表示：

二维随机欠采样轨迹下采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(Random GSCMRI)；

二维随机欠采样轨迹下采用双层伯格曼字典学习算法(Random TBMDU)；

二维随机欠采样轨迹下采用字典学习重建磁共振图像算法(Random DLMRI)；

笛卡尔采样轨迹下采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(Cartesian GSCMRI)；

笛卡尔采样轨迹下采用双层伯格曼字典学习算法(Cartesian TBMDU)；

笛卡尔采样轨迹下采用字典学习重建磁共振图像算法(Cartesian DLMRI)；

相位编码方向上随机欠采样K空间数据采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(LowResolution GSCMRI)；

相位编码方向上随机欠采样K空间数据采用双层伯格曼字典学习算法(LowResolution TBMDU)；

相位编码方向上随机欠采样K空间数据采用字典学习重建磁共振图像算法(LowResolution DLMRI)。

图6为7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种算法在二维随机欠采样轨迹，笛卡尔采样轨迹和相位编码方向上随机欠采样K空间数据的中间部分三种欠采样轨迹下高频误差(HFEN)随迭代次数(Interation Number)变化的曲线图。图6中9根曲线的意义同图5。

图7(a)为原图，(b)(c)(d)和(e)(f)(g)分别为二维随机欠采样轨迹在7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种方法下重建结果图和误差图。

图8为DLMRI、TBMDU和GTBMDU三种重建算法在不同标准偏差(Standard-Deviation)的高斯白噪声下的重建峰值信噪比值(PSNR)。

图9(a)为原图，(b)为二维随机欠采样模板，(c)(d)(e)和(f)(g)(h)分别是DLMRI、TBMDU和GSCMRI三种算法在σ＝5标准差的高斯白噪声时的重建结果图及误差图。

图10为4倍模拟径向欠采样轨迹下DLMRI、TBMDU和GSCMRI三种算法对水模数据的测试情况。其中(a)为原图，(b)为4倍欠采样的模拟径向轨迹欠采样轨迹模板，(c)(d)(e)分别为DLMRI、TBMDU和GSCMRI三种算法在该条件下的重建结果图，(f)(g)为相应部分结果的放大输出图。

图11为5倍欠采样率的笛卡尔采样轨迹下DLMRI，TBMDU和GSCMRI三种算法的重建效果，其中(a)为原图，(b)为5倍欠采样的相位编码方向上笛卡尔采样轨迹模板，(c)(d)(e)和(f)(g)(h)分别为DLMRI、TBMDU和GSCMRI三种算法在5倍欠采样下的重建结果图和残差图。

图12为GSCMRI算法下峰值信噪比(PSNR)随图正则化拉普拉斯算子权重η值的变化情况。

图13为GSCMRI算法下高频误差(HFEN)随图正则化拉普拉斯算子权重η值的变化情况。

图14(a)为原图，(b)(c)(d)分别为GSCMRI算法在η＝10^-1,10^-3,10^-5的重建结果图，(e)(f)(g)为相应的残差图。

具体实施方式

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

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

根据本发明的方法，本发明的技术方案通过图正则化稀疏编码模型并融入双层伯格曼迭代框架得到最终的成像结果。具体地说本发明实施例通过图正则化稀疏编码约束局部几何结构，可以更加精确地稀疏表示图像，处理几何特征更加复杂的图像，并恢复出更多的图像细节。现参考图1描述本发明的基于图正则化稀疏编码的磁共振超欠采样K数据成像方法。

步骤(a)：在双层伯格曼迭代框架上进行图正则化稀疏编码表示，建立图像稀疏模型。

基于先验信息进行约束的磁共振成像模型为以下目标函数：

$\begin{array}{cc}\underset{u}{\min }J\left(u\right)& s.t.{\left|\right|F_p-f\left|\right|}_2\≤{\mathrm{\ϵ}}_{\mathrm{\σ}}\end{array}---\left(1\right)$

此处代表需要重建图像，表示欠采样的傅里叶测量。通过部分采样的傅里叶编码矩阵映射图像u到频域空间f得到F_pu≈f。

本技术方案提出通过图正则化稀疏编码模型并融入双层的伯格曼迭代框架建立图像稀疏模型，为图像在学习字典上和图正则化上的稀疏编码表示，建立新的图像稀疏模型如下：

$\left\{\begin{array}{c}{u}^{k+1}=\arg \underset{u}{\min }\{\underset{D,\mathrm{\Γ}}{\min }\underset{i}{\mathrm{\Σ}}({\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|{\mathrm{D\α}}_i-R_{i}u\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right))+\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2\}\\ f^{k+1}=f^k+f-F_p{u}^{k+1}\end{array}\right.---\left(2\right)$

其中和λ均衡图像块的稀疏和更新字典的近似误差。对于许多的自然和医学图像，λ值可以凭经验来得到。第一个子问题的前三项利用这部分中的更新变量D，Γ来进行学习字典更新、稀疏系数编码。第四项保证重建结果与K空间欠采样数据保持匹配约束，λ，η和μ分别表示字典学习稀疏系数，图正则化编码系数和K数据拟合的权重。

其中，图正则化稀疏编码的目标函数由三部分线性组合构成，数据拟合经验项：拉普拉斯算子正则化项：ηΤr(ΓLΓ^Τ)；基于L₁范数的补偿项：其图像模型如下：

$\begin{array}{c}\underset{D,\mathrm{\Γ}}{\min }\frac{\mathrm{\λ}}{2}{\left|\right|X-\mathrm{D\Γ}\left|\right|}_F^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\left|\right|{\mathrm{\α}}_i\left|\right|}_1\\ s.t.{\left|\right|d_j\left|\right|}^2\≤1,\∀j=\mathrm{1,2},...,J\end{array}---\left(3\right)$

其中参数λ表示字典学习稀疏系数，η决定图表正则化拉普拉斯算子的权重。

图正则化约束在流形假设下可以有效地利用几何数据。其假设两个数据点x_i,x_j在数据分布的内在几何结构上相近，它们的表示系数α_i和α_j在新字典中也相近。给定的数据X，最近邻图G是由M个顶点构造的，每个顶点代表X中的一个数据点，W表示G的权重矩阵。如果x_i是在x_j的K接近邻域或者x_j是在x_i的K接近邻域，W_ij＝1；否则W_ij＝0。此外，x_i的程度得定义为和C＝diag(c₁,c₂,...,c_M)，L＝C-W。映射权重图G来稀疏表达系数Γ，选择合理映射的准则最小化下列目标函数：

$\frac{1}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{({\mathrm{\α}}_i-{\mathrm{\α}}_j)}^2{W}_{\mathrm{ij}}=\mathrm{Tr}\left({\mathrm{\ΓL\Γ}}^T\right)---\left(4\right)$

步骤(b)：引入辅助变量和轮换求解的技术，在双层伯格曼迭代的内层迭代上更新学习字典和稀疏系数。

定义X＝Ru，通过增加辅助变量，把无约束问题转化为约束问题，即

$\begin{array}{c}\underset{D,\mathrm{\Γ}}{\min }\frac{\mathrm{\λ}}{2}{\left|\right|Z\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\underset{i}{\overset{M}{\mathrm{\Σ}}}{\left|\right|{\mathrm{\α}}_i\left|\right|}_1\\ s.t.Z{=X-\mathrm{D\Γ},\left|\right|d_j\left|\right|}^2\≤1,j=\mathrm{1,2},...,J\end{array}---\left(5\right)$

技术方案中通过双层伯格曼迭代求解约束问题等式(5)，将其转化为一系列无约束问题，其中目标函数是由原目标函数加上约束的惩罚项组成，即

$(D^{k+1},{\mathrm{\α}}_i^{k+1},z_i^{k+1})=\arg \underset{D,{\mathrm{\α}}_i,z_i}{\min }{L}_{\mathrm{\β}}(D,{\mathrm{\α}}_i,z_i,y_i)---\left(6\right)$

其中 $\begin{array}{c}{L_{\mathrm{\β}}(D,{\mathrm{\α}}_i,z_i,y_i)=\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|z_i\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\frac{\mathrm{\β}}{2}{\left|\right|z_i+{\mathrm{D\α}}_i-R_{i}u-\frac{y_i^k}{\mathrm{\β}}\left|\right|}_2^2\\ Y^{k+1}=Y^k+\mathrm{\β}(\mathrm{Ru}-Z^{k+1}-D^{k+1}{\mathrm{\Γ}}^{k+1})\end{array}---\left(7\right)$

其中为乘变量，

求解二次多项式函数关于D的导数，得到如下的最速下降更新规则：

$D^{k+1}=D^k+{\mathrm{\ξY}}^{k+1}{\left({\mathrm{\Γ}}^{k+1}\right)}^T,d_i^k=\frac{d_i^k}{{\left|\right|d_i^{k+1}\left|\right|}_2},i=\mathrm{1,2},...,J---\left(8\right)$

经过最速下降迭代后，再将字典的每一列向量进行归一化处理使得D^k+1每一列的范数为1。本技术方案的字典更新的突出优势是每一次迭代可以看成是一次字典的加细操作，若把每一步看成是一个尺度，则字典的更新是从低尺度到高尺度的精细过程。即将“小尺度”(Y^k+1(Γ^k+1)^T)的原型加到一个字典中，该方式可以有效的避免陷入局部解。

本技术方案中采用的双层伯格曼迭代算法，使D^kΓ^k收敛到X直到最优化问题的约束满足为止。相应地，变量的D^kΓ^k稀疏度在初始迭代时非常大，随着迭代的进行逐渐的减少。

如图1所示，在步骤S102中，进行图正则化稀疏编码求解。为最小化二次多项式来表示Z，获得最优解如下：

$Z=\frac{\mathrm{\λX}+\mathrm{\β}(D^k\mathrm{\Γ}-Y^k/\mathrm{\β})}{\mathrm{\λ}+\mathrm{\β}}---\left(9\right)$

将最优解的(9)式的Z代入到等式(6)中，得到

${\mathrm{\α}}_i^{k+1}=\arg \underset{{\mathrm{\α}}_i}{\min }\frac{\mathrm{\λ\β}}{2(\mathrm{\λ}+\mathrm{\β})}{\left|\right|D^k{\mathrm{\α}}_i-x_i-y_i^k/\mathrm{\β}\left|\right|}_F^2+\mathrm{\η}{L}_{\mathrm{ii}}{\mathrm{\α}}_i^T{\mathrm{\α}}_i+{\mathrm{\α}}_i^T{h}_i+{\left|\right|{\mathrm{\α}}_i\left|\right|}_1---\left(10\right)$

式中再利用软阈值迭代算法求解即：

$\begin{array}{c}{\mathrm{\α}}_i^{m+1}=\arg \underset{{\mathrm{\α}}_i}{\min }\{\mathrm{\γ}{\left|\right|{\mathrm{\α}}_i-[{\mathrm{\α}}_i^m+\frac{-2\mathrm{\η}{L}_{\mathrm{ii}}{\mathrm{\α}}_i^m-h_i+{\left(D^k\right)}^T{y}_i^m}{2\mathrm{\γ}}]\left|\right|}_2^2+{\left|\right|{\mathrm{\α}}_i\left|\right|}_1\}\\ =\mathrm{shrink}({\mathrm{\α}}_i^m+\frac{-2\mathrm{\η}{L}_{\mathrm{ii}}{\mathrm{\α}}_i^m-h_i^m+{\left(D^k\right)}^T{y}_i^m}{2\mathrm{\γ}},\frac{1}{2\mathrm{\γ}})\end{array}---\left(11\right)$

其中， $\mathrm{\γ}\≥\left[\frac{\mathrm{\λ\β}}{2(\mathrm{\λ}+\mathrm{\β})}\right]\mathrm{eig}\left({\left(D^k\right)}^T{D}^k\right)+\mathrm{\α}{L}_{\mathrm{ii}}.$

步骤(c)：利用部分超欠采样的K数据约束，在双层伯格曼迭代框架的外层迭代上进行图像更新，以得到成像结果。

本技术方案提出的图正则化稀疏编码模型的外层伯格曼迭代如下：

$\left\{\begin{array}{c}{u}^{k+1}=\arg \underset{u}{\min }\{\underset{D,\mathrm{\Γ}}{\min }\underset{i}{\mathrm{\Σ}}({\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|{\mathrm{D\α}}_i-R_{i}u\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right))+\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2\}\\ f^{k+1}=f^k+f-F_p{u}^{k+1}\end{array}\right.---\left(12\right)$

方程(12)为外层伯格曼迭代，求解最小二乘解析问题，更新u^k+1，得到最终成像结果。

如图1所示，对于图像u的子问题，它在双层伯格曼迭代过程的每一次内部迭代之后进行更新。即在外层的伯格曼迭代中，图像u通过求解最小二乘解析问题得到：

$u^{k+1}=\arg \underset{u}{\min }\underset{i}{\mathrm{\Σ}}\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2+\frac{\mathrm{\β}}{2}{\left|\right|D^{k+1}{\mathrm{\α}}_i^{k+1}-R_{i}u+z_i^{k+1}-\frac{y_i^{k+1}}{\mathrm{\β}}\left|\right|}_2^2---\left(13\right)$

图像u的最小二乘解析如下：

$(\mathrm{\μ}{F}_p^T{F}_p+\mathrm{\β}\underset{i}{\mathrm{\Σ}}{R}_i^T{R}_i)u^{k+1}=\mathrm{\μ}{F}_p^T{f}^k+\mathrm{\β}\underset{i}{\mathrm{\Σ}}{R}_i^T(D^{k+1}{\mathrm{\α}}_i^{k+1}+z_i^{k+1}-y_i^{k+1}/\mathrm{\β})---\left(14\right)$

定义：

$S_2=\frac{F\underset{i}{\mathrm{\Σ}}{R}_i^T(D^{k+1}{\mathrm{\α}}_i^{k+1}+z_i^{k+1}-y_i^{k+1}/\mathrm{\β})}{\mathrm{\ω}}---\left(15\right)$

其中定义为全傅里叶编码矩阵，标准化如下F^ΤF＝1_N，Fu代表全K空间数据，Ω表示被采样数据的子集。其中平均每块图像结果并且变换到傅里叶域，则:

$\mathrm{Fu}(k_x,k_y)=\left\{\begin{array}{cc}{S}_2(k_x,k_y),& (k_x,k_y)\∉\mathrm{\Ω}\\ \frac{\mathrm{\μ}{S}_1(k_x,k_y)+\mathrm{\β}{S}_2(k_x,k_y)}{\mathrm{\μ}+\mathrm{\β}},& (k_x,k_y)\∈\mathrm{\Ω}\end{array}\right.---\left(16\right)$

具体而言，本技术方案提出的引入自适应字典的更新可以比预先构造的字典更好的稀疏表示图像以及图正则化稀疏编码方法可以有效的捕获图像的局部几何结构。通过图正则化局部结构约束的引入从而可以更精确的稀疏表示图像，并恢复出更多的图像细节，使得重建图像具有更好的保真度。

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

(1):初始化:Γ⁰＝0，D⁰，f⁰＝f

(2):若||F_pu^k-f||₂＞σ则进行(3)-(12)；否则迭代停止

(3):若j≤P则进行(4)-(11)；否则迭代停止

(4):当不满足迭代停止条件则进行(5)-(6)

(5): $Y^{m+1}=\frac{\mathrm{\λ\β}(X-D^m{\mathrm{\Γ}}^{m+1}+Y^m/\mathrm{\β})}{\mathrm{\λ}+\mathrm{\β}}$

(6): $\begin{array}{c}{\mathrm{\α}}_i^{j,m+1}=\mathrm{shrink}({\mathrm{\α}}_i^{j,m}+\frac{-2\mathrm{\η}{L}_{\mathrm{ii}}{\mathrm{\α}}_i^{j,m}-h_i^m+{\left(D^j\right)}^T{y}_i^{m+1}}{2\mathrm{\γ}},\frac{1}{2\mathrm{\γ}})\end{array}$

(7):C^j+1＝Y^m+1,Γ^j+1,0＝Γ^j,m+1,D^j+1＝D^j+ξC^j+1(Γ^j+1,0)^Τ

(8): $d_q^{j+1}=d_q^{j+1}/{\left|\right|d_q^{j+1}\left|\right|}_2,q=\mathrm{1,2},...,Q$

(9):通过频率插值更新u^k+1

(10): $S_2=\frac{F\underset{i}{\mathrm{\Σ}}{R}_i^T(D^{k+1}{\mathrm{\α}}_i^{k+1}+z_i^{k+1}-y_i^{k+1}/\mathrm{\β})}{\mathrm{\ω}}$

(11): $\mathrm{Fu}(k_x,k_y)=\left\{\begin{array}{cc}{S}_2(k_x,k_y),& (k_x,k_y)\∉\mathrm{\Ω}\\ \frac{\mathrm{\μ}{S}_1(k_x,k_y)+\mathrm{\β}{S}_2(k_x,k_y)}{\mathrm{\μ}+\mathrm{\β}},& (k_x,k_y)\∈\mathrm{\Ω}\end{array}\right.$

(12):f^k+1＝f^k+f-F_pu^k+1

本发明实施例中通过实值图像数据实验和复数值图像数据实验评价本技术算法的性能。本发明技术方案的采样方案包括二维随机采样，一维随机相位编码的笛卡尔采样和模拟径向欠采样。实值实验使用体内的MR扫描图像，大小为512×512。附图5和附图6复数值图像数据实验的图像大小为256×256和512×512。各参数的标准值设置和TBMDU方法中一样，图像块大小为过完备字典中K＝1(对应于J＝36)，重叠块r＝1(对应于J＝36和数据采样率为L＝267289时512×512的图像大小对的环境)，对于其它特定的参数，DLMRI和TBMDU方法设置为默认值。本发明实施例中实验实现过程，选取参数η＝10^-3。此外，引入峰值信噪比(PSNR)和高频误差(HFEN)来量化重建图像质量。

图2为模拟径向欠采样轨迹时DLMRI，TBMDU和GSCMRI三种算法下的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化。可以看出在不同的欠采样因子下GSCMRI的重建效果都相对另外两种重建算法较好。

图3为模拟径向欠采样轨迹时DLMRI、TBMDU和GSCMRI三种算法下的高频误差(HFEN)随着模拟径向欠采样轨迹欠采样因子(Downsampling Factor)的变化。可以看出在不同的欠采样因子下GSCMRI的重建效果都相对另外两种重建算法较好。

图4为模拟径向欠采样轨迹下的DLMRI，TBMDU和GSCMRI三种算法的重建效果的对比结果。(a)为原图，(b)(c)(d)和(e)(f)(g)分别是DLMRI、TBMDU和GSCMRI三种方法在8倍模拟径向轨迹欠采样率下的重建结果图和重建误差图。可以看出在同样的采样率下，GSCMRI在不同的采样轨迹下的重建效果都优于其它两种重建算法。

图5为7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种算法在二维随机欠采样轨迹，笛卡尔采样轨迹和相位编码方向上随机欠采样K空间数据的中间部分三种采样欠采样轨迹下的峰值信噪比(PSNR)随迭代次数(Interation Number)变化的曲线图。可以看出在同样的采样率下，GSCMRI在不同的采样轨迹下的重建效果都优于其它两种重建算法。图5中9根曲线依次分别表示：

二维随机欠采样轨迹下采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(Random GSCMRI)；

二维随机欠采样轨迹下采用双层伯格曼字典学习算法(Random TBMDU)；

二维随机欠采样轨迹下采用字典学习重建磁共振图像算法(Random DLMRI)；

笛卡尔采样轨迹下采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(Cartesian GSCMRI)；

笛卡尔采样轨迹下采用双层伯格曼字典学习算法(Cartesian TBMDU)；

笛卡尔采样轨迹下采用字典学习重建磁共振图像算法(Cartesian DLMRI)；

相位编码方向上随机欠采样K空间数据采用图正则化稀疏编码的磁共振超欠采样K数据成像算法(LowResolution GSCMRI)；

相位编码方向上随机欠采样K空间数据采用双层伯格曼字典学习算法(LowResolution TBMDU)；

相位编码方向上随机欠采样K空间数据采用字典学习重建磁共振图像算法(LowResolution DLMRI)。

图6为7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种算法在二维随机欠采样轨迹，笛卡尔采样轨迹和相位编码方向上随机欠采样K空间数据的中间部分三种欠采样轨迹下高频误差(HFEN)随迭代次数(Interation Number)变化的曲线图。可以看出在同样的采样率下，GSCMRI在不同的采样轨迹下的重建效果都优于其它两种重建算法。图6中9根曲线的意义同图5。

图7(a)为原图，(b)(c)(d)和(e)(f)(g)分别为二维随机欠采样轨迹在7.11倍欠采样率下DLMRI，TBMDU和GSCMRI三种方法下重建结果图和误差图。

图8为DLMRI、TBMDU和GTBMDU三种重建算法在不同标准差(Standard-Deviation)的高斯白噪声下的重建峰值信噪比值(PSNR)。

图9(a)为原图，(b)为二维随机欠采样模板，(c)(d)(e)和(f)(g)(h)分别是DLMRI、TBMDU和GSCMRI三种算法在σ＝5标准差的高斯白噪声时的重建结果图及误差图。可以看出，GSCMRI重建算法与参照图像误差最小，TBMDU算法次之。

图10为4倍模拟径向欠采样轨迹下DLMRI、TBMDU和GSCMRI三种算法对水模数据的测试情况。其中(a)为原图，(b)为4倍欠采样的模拟径向轨迹欠采样轨迹模板，(c)(d)(e)分别为DLMRI、TBMDU和GSCMRI三种算法在该条件下的重建结果图，(f)(g)为相应部分结果的放大输出图。从重建结果来看，GSCMRI方法实现更高分辨率重建。

图11为5倍欠采样率的笛卡尔采样轨迹下DLMRI，TBMDU和GSCMRI三种算法的重建效果，其中(a)为原图，(b)为5倍欠采样的相位编码方向上笛卡尔采样轨迹模板，(c)(d)(e)和(f)(g)(h)分别为DLMRI、TBMDU和GSCMRI三种算法在5倍欠采样下的重建结果图和残差图。从误差图可以看出GSCMRI方法的重建效果和参照图像非常接近。

图12为GSCMRI算法下峰值信噪比(PSNR)随图正则化拉普拉斯算子权重η值的变化情况。可以看出，η＝10^-3时，本发明算法有相对较好的重建效果。

图13为GSCMRI算法下高频误差(HFEN)随图正则化拉普拉斯算子权重η值的变化情况。可以看出，η＝10^-3时，本发明算法有相对较好的重建效果。

图14(a)为原图，(b)(c)(d)分别为GSCMRI算法在η＝10^-1,10^-3,10^-5的重建结果图，(e)(f)(g)为相应的残差图。

本发明实施例通过图正则化稀疏编码的引入约束局部几何结构，利用其在几何数据方面的优势，处理几何结构特征更复杂的图像，恢复出更多的图像细节，从而得到更精确的重建图像，使重建图像具有更好的保真度。

Claims (9)

1.一种基于图正则化稀疏编码的磁共振超欠采样K数据成像方法，其特征是包括以下步骤：

步骤(a)：在双层伯格曼迭代框架上进行图正则化稀疏编码表示，得到图像稀疏模型；

步骤(b)：利用引入辅助变量和轮换求解的技术，在双层伯格曼迭代的内层迭代上更新学习字典和稀疏系数；

步骤(c)：利用部分超欠采样的K数据约束，在双层伯格曼迭代的外层迭代上进行图像更新，以得到成像结果。

2.根据权利要求1所述的成像方法，其特征是步骤(a)所述在双层伯格曼迭代框架上进行图正则化稀疏编码表示，是对图像的系数进行稀疏表示的基础上，将图像数据的局部几何结构约束通过对系数图正则化以得到更好的稀疏表示，得到图像稀疏模型。

3.根据权利要求1所述的成像方法，其特征是在步骤(b)中引入辅助变量将步骤(a)中的无约束问题转化为约束问题求解，通过双层伯格曼内层迭代算法完成字典更新以及通过软阈值迭代更新稀疏系数。

4.根据权利要求1所述的成像方法，其特征是步骤(c)在步骤(b)完成字典学习和稀疏系数后，通过双层伯格曼外层迭代，求解最小二乘的解析问题，进行图像更新，得到图像重建结果。

5.根据权利要求2所述的成像方法，其特征是在双层伯格曼迭代框架的稀疏模型上融入图正则化稀疏编码过程，建立的图像稀疏模型为：

$\left\{\begin{array}{c}{u}^{k+1}=\arg \underset{u}{\min }\{\underset{D,\mathrm{\Γ}}{\min }\underset{i}{\mathrm{\Σ}}({\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|D{\mathrm{\α}}_i-R_{i}u\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right))+\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2\}\\ f^{k+1}=f^k+f-F_p{u}^{k+1}\end{array}\right.$

其中，模型中第一个子问题的前三项为图像在字典上进行稀疏和图正则化表示，u表示磁共振图像，D表示学习字典，α_l表示图像块稀疏系数，R_l表示图像块选取矩阵，Γ表示图正则化稀疏表达系数，F_p表示部分傅立叶变换，f表示K数据；第四项保证重建结果与K空间欠采样数据保持匹配约束，λ，η和μ分别表示学习字典稀疏系数，图正则化稀疏系数和K数据拟合权重。

6.根据权利要求5所述的成像方法，其特征是引入辅助变量，将图像模型的第一个子问题转化为：

$\underset{D,\mathrm{\Γ}}{\min }\frac{\mathrm{\λ}}{2}{\left|\right|Z\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\underset{i}{\overset{M}{\mathrm{\Σ}}}{\left|\right|{\mathrm{\α}}_i\left|\right|}_1$

s.t.Z＝Dα_i-R_iu,||d_j||²≤1,j＝1,2,...,J

7.根据权利要求6所述的成像方法，其特征是引入辅助变量后，采用双层伯格曼技术，将图像模型变更为：

$(D^{k+1},{\mathrm{\α}}_i^{k+1},z_i^{k+1})=\arg \underset{D,{\mathrm{\α}}_i,z_i}{\min }{\left|\right|{\mathrm{\α}}_i\left|\right|}_1+\frac{\mathrm{\λ}}{2}{\left|\right|z_i\left|\right|}_2^2+\mathrm{\ηTr}\left({\mathrm{\ΓL\Γ}}^T\right)+\frac{\mathrm{\β}}{2}{\left|\right|z_i+D{\mathrm{\α}}_i-R_{i}u-\frac{y_i^k}{\mathrm{\β}}\left|\right|}_2^2$

Y^k+1＝Y^k+β(Ru-Z^k+1-D^k+1Γ^k+1)

8.根据权利要求7所述的成像方法，其特征是所述的图像模型的求解是轮换地更新一个变量，同时固定其它变量；固定辅助变量z_i和系数α_i，通过最速下降法更新字典D^k+1；固定辅助字典D和系数α_i，通过最小化二次多项式更新固定辅助变量z_i和字典D^k，通过利用软阈值迭代算法更新

9.根据权利要求6、7或8所述的成像方法，其特征是对于磁共振图像u的子问题，它在双层伯格曼迭代过程的每一次内层迭代之后进行更新，在外层的伯格曼迭代中，磁共振图像u通过求解模型最小化得到：

$u^{k+1}=\arg \underset{u}{\min }\underset{i}{\mathrm{\Σ}}\frac{\mathrm{\μ}}{2}{\left|\right|F_{p}u-f^k\left|\right|}_2^2+\frac{\mathrm{\β}}{2}{\left|\right|D^{k+1}{\mathrm{\α}}_i^{k+1}-R_{i}u+z_i^{k+1}-\frac{y_i^{k+1}}{\mathrm{\β}}\left|\right|}_2^2$