一种基于广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法

摘要

一种基于广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法，包括：(a)在双层伯格曼字典学习框架上融入带非凸函数p范数先验信息进行字典学习和系数稀疏，建立图像稀疏表示模型；(b)在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新学习字典和稀疏系数，特别地利用广义软阈值迭代方法求解非凸p范数先验信息的目标函数，更新稀疏系数；(c)在双层伯格曼字典学习外层迭代上进行图像更新，得到重建图像。本发明通过广义软阈值迭代方法求解带非凸p范数先验信息的目标函数，可以在更大范围内惩罚小系数且对大系数偏差更小，进一步地稀疏表示图像，在少的扫描测量下精确地重建图像，减少重建图像的伪影，恢复更多图像细节。

说明书

技术领域

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

背景技术

磁共振成像是一种利用图像进行医疗诊断的检查技术，它没有电离辐射伤害，可以直接扫描人体的各个体层图像。因此磁共振成像可以提供给医生清晰实用的人体内部结构图像。

但是，磁共振成像系统的主要缺点是成像速度慢，使得磁共振成像检查的适应症大为减少，比如不适合于运动性器官的检查和危重病人的检查等，再者对于噪动或丧失自制能力的患者，如果不使用镇静剂难以成像，儿科的某些应用同样受到限制。因此自磁共振成像出现以来，人们一直致力于成像速度的提高。此外，磁共振成像速度慢与扫描时间有关。扫描时间又与采样率成正比，减少扫描时间的同时会使采样率相应下降，重建的图像分辨率也会降低。

然而压缩感知和稀疏表示理论在磁共振图像上的应用可以使得在不改变目前磁共振硬件系统的条件下减少成像时间成为可能。压缩感知理论阐明，图像在某些特定的变换域中能够稀疏表示，称之图像的稀疏变换。奈奎斯特采样定理要求在信号的采集阶段以高于信号带宽的两倍采样率来获取信号，信号才能得到更好的重构，而压缩感知则对信号的带宽不再作要求，进一步延伸至信号的稀疏性，满足条件的信号则可在远少于采样率的情况下精确的重构信号。

以稀疏正则化为特征的压缩感知理论是减少扫描量而保证高效重建的一种有效途径。稀疏表达字典学习理论特点是设计最优的自适应字典，使得重叠的图像块在这些字典下是稀疏表示的。鉴于稀疏编码模型在图像处理各方面，特别是在图像恢复中的良好效果，设计鲁棒的数值算法使字典学习领域的一个极为重要的问题。

求解信号的稀疏表达问题，分成两种类型：直接求解l₀范数优化的贪婪算法和l₁范数近似最小化方法。直接求解l₀问题的数值结果很不稳定，为一个NP难问题，解决这一问题的方法通常为寻找次优解的贪婪算法。贪婪算法的优点是计算简单、重构效果好，但大部分贪婪算法的前提条件是已知待重构信号的稀疏度。l₁范数最小化是通过用l₁范数来近似l₀范数，将问题转化为凸优化问题，通过求解线性规划问题进行稀疏信号的重构，进一步得到等价解。

相比于l₁范数近似最小化方法，目前有研究表明使用p(0＜p＜1)范数能更好地稀疏表示，联合其它更好的变换可以得到较好的重建结果。Chartrand等将非凸的p(0＜p＜1)范数用于磁共振图像重建，其比l₁范数更近似的表示l₀范数，能更好地稀疏表示图像。其实验证明了基于p范数优化算法在重构信号效果方面优于l₁范数优化算法。Wong等将半局部先验算子引入到同伦l₀最小化中稀疏重建磁共振图像。Shi等结合同伦l₀最小化和非凸极大极小凹惩罚学习恢复自然图像。

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

$\left(\begin{array}{c}\underset{u,D,\Gamma }{\min }\{\underset{l}{\Sigma }{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2+v{\left|\right|F_{p}u-f\left|\right|}_2^2\}\\ s.t.{\left|\right|{\alpha }_l\left|\right|}_0\le T_0,\forall l\end{array}\right)$

此处，Γ＝[α₁,α₂,...,α_L]为所有图像块对应的稀疏稀疏矩阵。DLMRI直接通过正交匹配追踪方法解决l₀最小化问题，虽然这些数据的学习方法比以前预定义为基础字典的方法有很大的改善，但l₀最小化问题是非凸的NP难问题。模型中前一项用于图像块在自适应学习字典上的稀疏表示，后一项是在图像数据的拟合保真项。求解该模型的两步迭代方法，第一步是自适应字典学习；第二步是从高度欠采样的K-空间数据重建图像。

刘且根等人提出以双层伯格曼算法为主要框架的字典学习模型，双层伯格曼字典学习算法(TBMDU)用l₁范数来近似l₀范数，致力于用软阈值迭代算法解决l₁范数稀疏最小化问题。其外层迭代与数据保真度这一项有关，内层是与字典和基于图像块稀疏有关，内层的稀疏系数主要是l₁范数稀疏。改进的稀疏编码和字典更新应用在内层伯格曼迭代，如此使整个算法在较少次迭代后就能达到收敛。其算法模型如下：

$\left(\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\{\underset{D,\Gamma }{\min }\underset{l}{\Sigma }({\left|\right|{\alpha }_l\left|\right|}_1+\frac{\lambda }{2}{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2)+\frac{\mu }{2}{\left|\right|F_{p}u-f^i\left|\right|}_2^2\}\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right)$

式中D＝[d₁,d₂,…,d_J]∈C^M×J，Γ＝[α₁,α₂,…,α_L]∈C^J×L。λ代表在最优字典中图像块的稀疏水平。对于医学图像，λ的值可以凭经验来得到。

在磁共振成像中，变换后的大系数决定成像的高频成分，而系数的高频成分决定成像的细节部分和伪迹。现有技术的缺陷在于，贪婪算法解决l₀范数优化问题，是一个NP难问题，并且大部分贪婪算法需要知道重构信号的稀疏度，这在实际工程实施中是很难满足的。Avishankar提出的字典学习算法也是使用贪婪算法求解l₀最小化问题。刘且根等人提出的双层伯格曼字典学习算法用l₁范数来近似l₀范数求得等价解，但是其变换后小系数近似为零且大系数的偏差较大，会使重建图像解剖结构不清晰且出现混叠现象。因此业界需要一种更好的稀疏表示方法稀疏系数，以更好的描述解剖结构，重建更精确的磁共振图像。

发明内容

本发明的目的是提出一种基于广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法(GTBMDU)。

本发明实施例广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法中通过广义软阈值迭代方法求解带非凸p范数先验信息的目标函数，相对于凸1范数约束，使用非凸p范数约束可以在更大范围内惩罚小系数且对大系数偏差更小，更好的保留大系数而小系数近似为零，因此进一步的稀疏表示图像，并且可以在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。

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

步骤(a)：在双层伯格曼字典学习框架上融入带非凸函数p范数先验信息进行字典学习和系数稀疏，建立图像稀疏表示模型。

步骤(b)：在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新学习字典和稀疏系数，特别地利用广义软阈值迭代方法求解非凸p范数先验信息的目标函数，更新稀疏系数。

步骤(c)：在双层伯格曼字典学习外层迭代上进行图像更新，得到重建图像。

进一步地说，本发明所述步骤(a)为：在双层伯格曼迭代框架上通过带非凸p范数进行字典学习和稀疏编码表示，相对于凸1范数约束，使用非凸p范数约束可以在更大范围内惩罚小系数且对大系数偏差更小，更好地保留大系数而小系数近似为零，因此进一步地稀疏表示图像，得到图像稀疏表示模型。

进一步地说，本发明所述步骤(b)为：用增加辅助变量将步骤(a)中的无约束问题转化为约束问题求解，通过双层伯格曼字典学习算法内层迭代完成字典更新和系数稀疏，特别地通过广义软阈值迭代算法求解带非凸p范数先验信息约束的目标函数，以更新稀疏系数。

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

进一步地说，本发明所述步骤(a)为：在双层伯格曼字典学习迭代框架上融入非凸p范数约束稀疏系数，建立的图像稀疏模型为：

$\left(\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\{\underset{D,\Gamma }{\min }\underset{l}{\Sigma }(g_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2)+\frac{\mu }{2}{\left|\right|F_{p}u-f^i\left|\right|}_2^2\}\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right)$

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

进一步地说，本发明所述步骤(b)为：增加辅助变量，把非约束问题转化为约束问题，则图像稀疏模型的第一个子问题转化为：

$\left(\begin{array}{cc}\underset{D,{\alpha }_l}{\min }\left(g_{\tau ,p}{\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|z_l\left|\right|}_2^{2}s.t& z_l=D{\alpha }_l-R_{l}u,l=\mathrm{1,2},...,L\end{array}\right)$

进一步地说，本发明所述步骤(b)为:引入辅助变量后，采用分离双层伯格曼迭代技术，将图像模型变为：

$\left(\begin{array}{c}(D^{i+1},{\alpha }_l^{i+1},z_l^{i+1})=\arg \underset{D,{\alpha }_l,z_l}{\min }\underset{l}{\Sigma }{g}_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|z_l\left|\right|}_2^2+\frac{\beta }{2}{\left|\right|z_l+D{\alpha }_l-R_{l}u-y_l^i/\beta \left|\right|}_2^2\\ Y^{i+1}=Y^i+\beta (\mathrm{Ru}-z^{i+1}-D^{i+1}{\Gamma }^{i+1})\end{array}\right)$

其中λ表示在最优字典中图像块的稀疏水平。对于医学图像，λ的值可以凭经验来得到。

进一步地说，本发明所述步骤(b)为轮换地更新一个变量，同时固定其它变量：固定辅助变量z_l和系数α_l，通过最速下降法更新字典Dⁱ⁺¹；固定辅助字典D和系数α_l，通过最小化二次多项式更新固定辅助变量z_l和字典D，通过利用广义软阈值迭代算法求解非凸范数g_τ,p(α_l)的目标函数，更新稀疏系数

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

进一步地说，本发明所述步骤(c)为在双层伯格曼字典学习迭代框架上通过广义软阈值迭代算法求解带非凸p范数先验信息约束的目标函数，重建磁共振超欠采样K数据，得到最终成像结果。

本发明的技术方案具有以下的优点或有益效果：本发明实施例广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法中通过广义软阈值迭代方法求解带非凸l_p范数先验信息的目标函数。在磁共振成像中，图像经过傅里叶变换，在频域内的小系数决定图像的细节和重建图像的伪影。本技术方案发明方法相对于凸1范数约束，使用非凸p范数约束可以在更大范围内惩罚小系数且对大系数偏差更小，更好地保留大系数而小系数近似为零，因此进一步地稀疏表示图像，在更少的测量数据下更精确的重建图像，恢复更多的图像细节。

由于本发明的算法可以在更大范围内惩罚小系数且对大系数偏差更小，更好地保留大系数而小系数近似为零，因此进一步地稀疏表示图像，使得重构的图像达到令人满意的效果，因此能够快速、精确的重构磁共振图像。

附图说明

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

图2为同一采样率下不同采样轨迹模板的DLMRI和GTBMDU的重建结果。其中：(a)为原图，(b)为模拟螺旋采样轨迹模板，(c)为相位编码方向上的随机欠采样，(d)为模拟径向采样模板，其采样率为86％，(e)(f)(g)为DLMRI的在三种采样轨迹模板下的重建误差图，(h)(i)(j)为p＝0.5时GTBMDU的在三种采样轨迹下的重建误差图。

图3为模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法重建图像的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化情况。

图4为模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法重建图像的高频误差(HFEN)随欠采样因子(Downsampling Factor)的变化情况。

图5为10倍欠采样率的模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，TBMDU， GTBMDU三种算法10倍欠采样率下重建误差图。

图6为模拟径向采样轨迹下不同复杂程度的高斯白噪声下DLMRI，TBMDU，GTBMDU三种算法的重建结果。其中：(a)为原图，(b)为6.09倍的模拟径向采样轨迹，(c)(d)(e)为DLMRI、TBMDU和GTBMDU三种算法在σ＝14时的重建误差图。

图7为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，TBMDU，GTBMDU三种重建算法的峰值信噪比(PSNR)值，当σ＝2时，DLMRI重建图像的峰值信噪比(PSNR)为40.16dB、TBMDU重建图像的峰值信噪比(PSNR)为42.74dB和GTBMDU在p＝0.5重建图像的峰值信噪比(PSNR)为43.41dB。

图8为磁共振系统评估物理模型的比较。其中：(a)为全采样重建图像，(b)为相位编码方向上80％的随机欠采样模板；(c)(d)(e)分别为DLMRI、TBMDU和GTBMDU三种算法的重建图像，其中GTBMDU重建时p＝0.7；(f)为(a)(c)(d)(e)的放大图像。

图9为本技术方案所提出GTBMDU算法在复杂K空间数据下的噪声敏感度测试，其中噪声的标准偏差为σ＝30。(a)(b)(c)分别为DLMRI、TBMDU和GTBMDU三种算法在5倍欠采样率的重建图像，其中GTBMDU重建时p＝0.7，三种算法对应的重建图像的峰值信噪比分别为17.93dB、22.94dB和25.42dB；(d)(e)分别为(a)(b)(c)的部分放大图像。

图10为实值头部数据的重建结果的对比图。其中(a)为全采样重建图；(b)为欠采样80％的相位编码方向上的随机欠采样模板；(c)(d)为DLMRI和GTBMDU两种算法在欠采样80％的重建图像，其中GTBMDU重建时p＝0.7；(e)为(a)(c)(d)的扩展图。

具体实施方式

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

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

现参考附图1描述根据本发明的基于广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法。根据本发明的方法，本发明技术方案在双层伯格曼迭代框架字典学习算法上，利用广义软阈值迭代算法对非凸函数最优化求解的优势，得到更好的图像重建效果。本发明实施例广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法中通过广义软阈值方法求解非凸函数l_p，其在大范围内惩罚小系数而大系数的偏差更小，使得对大系数的影响较小而小系数近似为零，因此进一步地稀疏表示图像，并且可以在更少的测量K数据下更精确的重建图像，恢复更多的图像细节。本发明技术方案算法的重建结果可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述。

步骤(a)：在双层伯格曼字典学习框架上融入带非凸函数p范数先验信息进行字典学习和系数稀疏，建立图像稀疏表示模型：

$\left(\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\{\underset{D,\Gamma }{\min }\underset{l}{\Sigma }(g_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2)+\frac{\mu }{2}{\left|\right|F_{p}u-f^i\left|\right|}_2^2\}\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right)---\left(1\right)$

前两项作为磁共振图像字典稀疏和非凸函数p范数稀疏表示的正则项，如下式所示：

$J\left(u\right)=\underset{l}{\Sigma }(g_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2)---\left(2\right)$

(1)式中第三项保证重建结果与K空间欠采样数据保持匹配约束，λ表示在最优字典中图像块的稀疏水平，μ表示K数据拟合的权重。

在磁共振成像中，图像数据经过傅里叶变换到频率域得到一系列变换系数K数据，当使用凸优化l₁范数近似l₀范数时，会使一部分小系数近似为零，大系数会出现一定的偏差，由于这些部分缺失的系数使得重建的磁共振图像会出现拖尾现象，产生混叠效应，重建图像解剖结构不清晰，对比度不明显。因此在磁共振成像中需要尽可能的减少变换后系数的偏差，保留相应解剖结构的系数。

目前的算法中通过软阈值迭代算法求解凸优化的l₁范数，传统的软阈值算子定义为：

Shrink(x,τ)＝sgn(x)·max{0,|x|-τ}      (3) 

式中对每一个输入值都有一个对应的输出值。相应的l₁最小化问题函数||Ax-b||²₂+2τ||x||₁能够通过简单的软阈值迭代算法来求解。通过最近邻和梯度下降法技巧可以得到近似求解如下：

xⁿ⁺¹＝Shrink(xⁿ-τA^T(Ax-b)/γ,τ/γ)       (4) 

其中γ＞eig((A)^TA)。

广义软阈值迭代算法对逆问题的解决非常有优势，广义软阈值迭代算法在_p值减少的情况下对小系数惩罚多，同时当p＜1时对大系数影响较小，使得处理大系数的偏差更小。使用改进的广义软阈值迭代算法使非凸函数g_τ,p(x)最小化。类比的，广义软阈值算子定义为：

GShrink(x,τ,p)＝sgn(x)·max{0,|x|-τ|x|^p-1},0＜p≤1     (5) 

由上式(5)可知，阈值函数GShrink满足一个非凸最小化函数与软阈值迭代算法相似，函数可以通过广义软阈值迭代算法最小化，得到以下式子：

$x^{n+1}=\mathrm{Shrink}(x^n-\tau A^T(\mathrm{Ax}-b)/\gamma ,\frac{\tau }{\gamma }\frac{1}{{\left|x^n\right|}^{1-p}+ϵ}),---\left(6\right)$

其中，γ＞eig((A)^TA)，GISTA在p＝1时与软阈值迭代算法相等。

步骤(b)：在双层伯格曼迭字典学习内层迭代上利用增加辅助变量和轮换技术更新学习字典和稀疏系数，特别地利用广义软阈值迭代方法求解带非凸p范数先验信息的目标函数，更新稀疏系数。通过下式更新学习字典，对于模型等式(1)通过增加辅助变量把非约束问题转化为约束形式：

$\left(\begin{array}{cc}\underset{D,{\alpha }_l}{\min }\left(g_{\tau ,p}{\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|z_l\left|\right|}_2^{2}s.t& z_l=D{\alpha }_l-R_{l}u,l=\mathrm{1,2},...,L\end{array}\right)---\left(7\right)$

用增广拉格朗日方法解决上式(7)：

$\left(\begin{array}{c}(D^{i+1},{\alpha }_l^{i+1},z_l^{i+1})=\arg \underset{D,{\alpha }_l,z_l}{\min }\underset{l}{\Sigma }{g}_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|z_l\left|\right|}_2^2+\frac{\beta }{2}{\left|\right|z_l+D{\alpha }_l-R_{l}u-y_l^i/\beta \left|\right|}_2^2\\ Y^{i+1}=Y^i+\beta (\mathrm{Ru}-z^{i+1}-D^{i+1}{\Gamma }^{i+1})\end{array}\right)---\left(8\right)$

其中λ表示最优字典中图像块的稀疏水平。对于医学图像，λ的值可以凭经验来得到。求解二次多项式函数关于D的导数，得到如下的最速下降更新规则：

$D^{i+1}=D^i+\xi Y^{i+1}{\left({\Gamma }^{i+1}\right)}^T,d_q^{i+1}=d_q^{i+1}/{\left|\right|d_q^{i+1}\left|\right|}_2,q=\mathrm{1,2},...,Q.---\left(9\right)$

求解等式(8)通过轮换地更新一个变量，同时固定其它变量：固定辅助变量z_l和系数α_l，通过最速下降法更新字典Dⁱ⁺¹；固定辅助字典D和系数α_l，通过最小化二次多项式更新固定辅助变量z_l和字典D，通过利用广义软阈值迭代算法求解非凸函数g_τ,p(α_l)更新

由等式(8)的可知，在学习字典下需要更好的稀疏矩阵以得到更精确的重建图像。利用非凸函数g_τ,p(α_l)可以更近似的表示l₀范数，相应地仅需要更少的测量即可精确的重建图像。利用广义软阈值迭代算法的双层伯格曼字典学习模型求解等式(8)，得到更新系数:

${\alpha }_l^{i,m+1}=\arg \underset{{\alpha }_l}{\min }{\left|\right|D^i{\alpha }_l-b_l-y_l^i/\beta \left|\right|}_1^2+2\frac{\lambda +\beta }{\mathrm{\lambda \beta }}{g}_{\tau ,p}\left({\alpha }_l\right)---\left(10\right)$

函数满足广义软阈值函数，得到α_l的解：

${\alpha }_l^{i,m+1}=\mathrm{Shrink}({\alpha }_l^{i,m}+{\left(D^i\right)}^T{y}_l^{m+1}/\gamma T_0,\frac{1}{\gamma T_0}\frac{1}{{\left|{\alpha }_{q,l}^{i,m}\right|}^{1-p}+ϵ})---\left(11\right)$

其中T₀＝λβ/(λ+β)，b_l＝R_lu。此处可以在大范围内惩罚小系数且使变换后大系数的偏差更小。

步骤(c)：在双层伯格曼迭代外层迭代上通过最小二乘解析问题求解，进行图像更新，得到重建图像。由步骤S102可知，通过广义软阈值迭代算法使得大系数有更好的保真度，由此类推由经过傅里叶反变换得到更精确的图像数据，进而减少重建图像的伪迹，得到更精细的解剖结构图，更精确的重建图像。本技术方案提出的广义双层伯格曼迭代算法解决非凸型磁共振超欠采样K数据重建模型的外层伯格曼迭代如下：

$\left(\begin{array}{c}{u}^{i+1}=\arg \underset{u}{\min }\{\underset{D,\Gamma }{\min }\underset{l}{\Sigma }(g_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|D{\alpha }_l-R_{l}u\left|\right|}_2^2)+\frac{\mu }{2}{\left|\right|F_{p}u-f^i\left|\right|}_2^2\}\\ f^{i+1}=f^i+f-F_p{u}^{i+1}\end{array}\right)---\left(12\right)$

上式(12)为广义双层伯格曼迭代的外层迭代，通过求解最小二乘解析问题，更新图像uⁱ⁺¹，以得到最终的成像结果。图像u通过求解下面式子的最小化得到：

$\left(\begin{array}{c}{u}^{i+1}=\underset{u}{\arg \min }\underset{l}{\Sigma }\left((g_{\tau ,p}\left({\alpha }_l\right)+\frac{\lambda }{2}{\left|\right|z_l\left|\right|}_2^2+\frac{\mu }{2}{\left|\right|F_{p}u-f^i\left|\right|}_2^2)\right)\\ s.t.R_{l}u=D{\alpha }_l+z_l\end{array}\right)---\left(13\right)$

固定根据前面得到的D,α,z，得到更新的图像uⁱ⁺¹：

通过最小二乘解析上式得到：

$(\mu F_P^T{F}_P+\beta \underset{l}{\Sigma }{R}_l^T{R}_l)u^{i+1}=\mu F_P^T{f}^i+\beta \underset{l}{\Sigma }{R}_l^T(D^{i+1}{\alpha }_l^{i+1}+z_l^{i+1}-y_l^{i+1}/\beta )---\left(15\right)$

利用F^TF＝1_N，代入式子(15)得到：

$(\mathrm{\mu F}{F}_P^T{F}_P{F}^T+\mathrm{\beta F}\underset{l}{\Sigma }{R}_l^T{R}_l{F}^T)F{u}^{i+1}=\mathrm{\mu F}{F}_P^T{f}^i+\mathrm{\beta F}\underset{l}{\Sigma }{R}_l^T(D^{i+1}{\alpha }_l^{i+1}+z_l^{i+1}-y_l^{i+1}/\beta )---\left(16\right)$

矩阵FF_p^TF_pF^T是一个0和1的对角矩阵。用Γ⁰＝0代表D⁰，定义

$S_2=\frac{F\underset{l}{\Sigma }{R}_l^T(D^{i+1}{\alpha }_l^{i+1}+z_l^{i+1}-y_l^{i+1}/\beta )}{\omega }---\left(17\right)$

于是式子(16)可写成

$\mathrm{Fu}(k_x,k_y)=\left(\begin{array}{cc}{S}_2(k_x,k_y),& (k_x,k_y)\notin \Omega \\ [\mu S_1(k_x,k_y)+\beta S_2(k_x,k_y)]/(\mu +\beta ),& (k_x,k_y)\in \Omega \end{array}\right)---\left(18\right)$

综上所述在双层伯格曼字典学习迭代框架上通过广义软阈值迭代算法求解带非凸p范数先验信息约束的目标函数，重建磁共振超欠采样K数据，得到最终成像结果。

具体而言，本发明实施例广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法中通过广义软阈值方法求解带非凸p范数先验信息约束的目标函数，相对于1范数，非凸p范数在大范围内惩罚小系数而且使得在p＜1时大系数的偏差较小，进而对大系数的影响更小而小系数近似为零，因此进一步地稀疏表示图像。和常规的双层伯格曼迭代算法类似，其利用轮换求解技术更新图像u，可以在更少的测量数据下更精确地重建图像，恢复更多的图像细节。

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

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

(2):若||F_puⁱ-f||₂＞σ,则进行(3)-(12)，否则停止迭代

(3):若i＝1，i≤P则进行(4)-(11)，否则停止迭代

(4):若不满足迭代条件，则进行(5)-(6)

(5):Y^m+1＝T₀(-DⁱΓ^i,m+Ru^k,i+Cⁱ/β)

(6):${\alpha }_l^{i,m+1}=\mathrm{Shrink}({\alpha }_l^{i,m}+{\left(D^i\right)}^T{y}_l^{m+1}/\gamma T_0,\frac{1}{\gamma T_0}·\frac{1}{{\left|{\alpha }_{q,l}^{i,m}\right|}^{1-p}+ϵ})$

(7):Cⁱ⁺¹＝Y^m+1,Γ^i+1,0＝Γ^i,m+1,Dⁱ⁺¹＝Dⁱ+ξCⁱ⁺¹(Γ^i+1,0)^T,

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

(9):通过在频率域插值更新uⁱ⁺¹，

(10):$S_2=F\underset{l}{\Sigma }{R}_l^T(D^{i+1}{\alpha }_l^{i+1}+z_l^{i+1}-y_l^{i+1}/\beta )/\omega$

(11):$\mathrm{Fu}(k_x,k_y)=\left(\begin{array}{cc}{S}_2(k_x,k_y),& (k_x,k_y)\notin \Omega \\ [\mu S_1(k_x,k_y)+\beta S_2(k_x,k_y)]/(\mu +\beta ),& (k_x,k_y)\in \Omega \end{array}\right)$

(12):fⁱ⁺¹＝fⁱ+f-F_puⁱ⁺¹

本发明技术方案采用各种不同的欠采样因子和噪声对所提出方法的性能进行评估。本发明提出的GTBMDU方法与DLMRI和双层伯格曼字典学习方法TBMDU相比，DLMRI直接通过正交匹配追踪方法解决l₀最小化问题，TBMDU则致力于l₁解决稀疏最小化问题，GTBMDU方法则通过利用广义软阈值迭代算法求解非凸p范数先验信息的目标函数。实验过程各种参数的标准值分别设置如下：字典块尺寸字典的过完备性K＝1(相当于J＝36)，字典块的步长r＝1，因为所以总采样数据为L＝262144。参数设置为DLMRI的默认值。本发明技术方案中用到其他参数β＝0.0056,μ＝Nβ/M,ζ＝0.01,λ＝12,m＝3。为避免除零现象，ε设为5的初始值基础上每经过一次内部迭代则减少2％。重建图像的质量通过使用峰值信噪比(PSNR)和高频误差(HFEN)来衡量。

如下表格1列出DLMRI，TBMDU，GTBMDU三种算法在相同欠采样率下不同采样轨迹下的重建结果的峰值信噪比。可以看出，在不同的p值下，GTBMDU重建结果PSNR值不同。为平衡计算量，GTBMDU算法下当p＝0.5或者p＝0.7为较好的选择。

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

图2为同一采样率下不同采样轨迹模板的DLMRI和GTBMDU的重建结果。其中：(a)为原图，(b)为模拟螺旋采样轨迹模板，(c)为相位编码方向上的随机欠采样，(d)为模 拟径向采样模板，其采样率均为86％，(e)(f)(g)为DLMRI的在三种采样轨迹下的重建误差图，(h)(i)(j)为p＝0.5时GTBMDU的在三种采样轨迹下的重建误差图。

图3为模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法重建图像的峰值信噪比(PSNR)随欠采样因子(Downsampling Factor)的变化情况。从中可以看出TBMDU与GTBMDU两种算法重建的PSNR值之间的差距，GTBMDU算法有较好的重建性能，可以更好的重建出图像细节。

图4为模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法重建图像的高频误差(HFEN)随欠采样因子(Downsampling Factor)的变化情况。特别地，当欠采样因子增加到20倍，非凸优化的优势降低，三种算法都不能从较少的K数据采样中较好的重建图像。

图5为10倍欠采样率的模拟径向采样轨迹下DLMRI，TBMDU，GTBMDU三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，TBMDU，GTBMDU三种算法10倍欠采样率下重建误差图。

图6为模拟径向采样轨迹下不同复杂程度的高斯白噪声下DLMRI，TBMDU，GTBMDU三种算法的重建结果。其中：(a)为原图，(b)为6.09倍的模拟径向采样轨迹，(c)(d)(e)为DLMRI、TBMDU和GTBMDU三种算法在高斯白噪声标准偏差σ＝14时的重建误差图。可以看出，本发明技术方案的GTBMDU重建结果图可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述。

图7为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，TBMDU，GTBMDU三种重建算法的峰值信噪比(PSNR)值，当σ＝2时，DLMRI重建图像的峰值信噪比为40.16dB、TBMDU重建图像的峰值信噪比为42.74dB和GTBMDU在p＝0.5重建图像的峰值信噪比为43.41dB。

图8为磁共振系统评估物理模型的比较。其中：(a)为全采样重建图像，(b)为相位编码方向上80％的随机欠采样模板；(c)(d)(e)分别为DLMRI、TBMDU和GTBMDU三种算法的重建图像，其中GTBMDU重建时p＝0.7；(f)为(a)(c)(d)(e)的放大图像，可以看出GTBMDU算法在零值填充区域可以更好的消除混叠伪影。

图9为本技术方案所提出GTBMDU算法在复杂K空间数据下的噪声敏感度测试，其中噪声的标准偏差为σ＝30。(a)(b)(c)分别为DLMRI、TBMDU和GTBMDU三种算法在5倍欠采样率的重建图像，其中GTBMDU重建时p＝0.7，三种算法对应的重建图像的峰值信噪比分别为17.93dB、22.94dB和25.42dB；(d)(e)分别为(a)(b)(c)的部分放大图像。从三种重建图像可以看出DLMRI算法比其他两种方法有更严重的重建伪影，本实施例GTBMDU方法有更好的重建结果。

图10为实值头部数据的重建结果的对比图。其中(a)为全采样重建图；(b)为欠采样80％的相位编码方向上的随机欠采样模板；(c)(d)为DLMRI和GTBMDU两种算法在欠采样80％的重建图像，其中GTBMDU重建时p＝0.7；(e)为(a)(c)(d)的扩展图。可以看出GTBMDU在相位编码方向上比DLMRI呈现更少的混叠伪影，恢复的图像有更 好的保真度。

发明实施例广义双层伯格曼非凸型字典学习的磁共振超欠采样K数据成像方法中通过广义软阈值迭代方法求解带非凸p范数先验信息的目标函数，相对于凸1范数约束，使用非凸p范数约束可以在更大范围内惩罚小系数且对大系数偏差更小，更好地保留大系数而小系数近似为零，因此进一步地稀疏表示图像，可以在更少的测量数据下得到更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。本发明实施例算法重建结果可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述，恢复的图像有更好的保真度。

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