引言

磁共振成像（magnetic resonance imaging，MRI）技术是利用特定原子核（如氢核）的核磁共振（nuclear magnetic resonance，NMR）现象，通过定量检测物质中特定原子核共振过程中射频能量吸收和发射的规律，来显示物体内部结构及特征的一种成像技术.MRI扫描仪根据主磁场强度的不同，可分为低场、中场、高场、超高场等.一般情况下，低场≤0.3 T、0.3 T≤中场≤1.0 T、1.0 T≤高场≤3 T、超高场≥3 T.高场和超高场磁共振设备磁源多采用超导磁体产生高强磁场，采集到的信号幅值强、信噪比高.但超导磁体对设备环境要求高，且购买、安装和维护费用非常昂贵^[1].中场磁共振设备磁源采用超导磁体、永磁体或电磁，采集到的信号信噪比较高，一般用于医用MRI等.低场磁共振设备的磁源一般采用永磁体或电磁（室温下可使用），采集到的信号幅值很弱、信噪比极低.但是与中高场相比，低场磁共振设备具有便携、廉价、制造维护应用成本低，测量方便、快速、无损、无副作用等优点，已被广泛应用在食品化工、在线无损检测、聚合物材料、灾害防治、生物医药、能源勘探、石油测井、岩心分析等各种领域^[2]，并逐渐成为一种必要的分析检查手段.近年来，随着电磁技术和计算机技术的发展，低场磁共振成像（low-field MRI，LF MRI）也迅速发展起来，高性能的低场开放型永磁MRI系统开始受到人们的青睐^[3].本文旨在研究低场磁共振图像重建算法，改善其图像质量、提高其成像速率.

压缩感知（compressive sensing，CS）是一个新的信号采样和压缩理论.它利用信号在特定变换域上的稀疏性及可压缩性，能在远小于Nyquist采样率的条件下，从少量的信号中通过非线性重建算法恢复出原始信号^[4].传统的MRI技术需要很长的扫描时间，重建图像的时间较长.而将CS理论应用于MRI（CS-MRI）能够极大地减少傅里叶变换域的采样数据，降低扫描时间；且CS-MRI在不显著降低图像质量的情况下克服了MRI成像速度慢的缺点；验证了从少量稀疏测量数据重建出图像的可能性^[5].本文提出的低场CS-MRI重建技术借鉴了现已发展成熟的中高场CS-MRI技术，吸收其优点，并规避其缺点.

对于CS-MRI，有两个关键点需要进一步研究.第一：稀疏变换，基于CS理论的MRI技术能从欠采样的稀疏图像数据中重建出原始图像，而磁共振图像本身不具有稀疏性，因此需要将图像在稀疏域进行稀疏表示.传统CS-MRI常用的稀疏基有小波变换、全变分（total variation，TV）等，然而它们有时不能提供足够的稀疏表示.如：小波变换（wavelet transform，WT）有以下缺点：具有转移敏感性、缺乏方向性、缺乏相位信息，且不能稀疏表示二维信号的一维奇异点，很难捕捉到规则的图像轮廓^[6].然而，通过改善可以弥补WT的缺点，如文献[7]中使用双树小波变换代替WT作为稀疏基，得到了好的稀疏效果.第二：重建算法，由于CS-MRI的重建过程是一个病态问题，一般使用正则约束项去求得病态问题的稳定解.例如：文献[8]提出利用共轭梯度算法去解TV和小波稀疏联合的正则化目标函数问题.文献[9]提出了一种用来解决与文献[5]同样问题的部分傅里叶系数重建（reconstruction partial Fourier，RecPF）算法.文献[10]提出一种将TV和CS应用于MRI的算法（TVCMRI），它采用一种算子分裂算法求解重建问题.文献[11]提出的快速复合分裂算法（fast composite splitting algorithm，FCSA）将正则化问题分裂为两个简单的子问题，每个子问题再通过快速收缩算法求解^{[11, 12]}，使正则化问题更有效的得到解决.最近，一种先进的小波树稀疏MRI重建算法（wavelet tree sparsity MRI，WaTMRI）^[13]也已经被提出，它在TV和小波稀疏联合正则项基础上增加了小波树稀疏正则项.与先前的算法相比，WaTMRI算法使测量的复杂性由O(K+Klogn)减少到O(K+logn) ^[14].其中K为信号稀疏度，n为原始信号的大小.

基于以上理论，本文提出了一种双树小波变换与小波树稀疏联合的低场CS-MRI算法，称为DT-WaTMRI，在提高成像速度的同时重建出高分辨率的磁共振图像.以下部分是对双树小波变换与小波树稀疏联合的低场CS-MRI算法的详细阐述.

1 理论与算法 1.1 CS-MRI重建模型

CS的基本问题是从欠采样的线性观测向量$y=\mathit{\Phi }x$中恢复出信号x，其中$y\in {{R}^{m}}$、$x\in {{R}^{n}}$、$\mathit{\Phi }\in {{R}^{m\times n}}$，$m < n$（Φ是测量矩阵，R表示实数，n、m分别表示信号x、y的维度大小）；从数学角度求解以上欠定线性系统有无穷多解.然而，根据CS理论，在稀疏的假设下信号能通过(1)式所示的优化问题被恢复.

$ \underset{x}{\mathop{\text{min}}}\, {{\left\| x \right\|}_{0}}\ \ \ \text{s}\text{.t}.\ \ \ y=\mathit{\Phi }x $

(1)

其中${{\left\| x \right\|}_{0}}$是l₀范数，表示向量$\left\| x \right\|$中值不为0的个数.由于最小l₀范数的优化问题是NP-hard问题，难以解决，Chrétien等人^[15]证明可以利用凸的l₁来近似非凸l₀范数，如(2)式所示.

$ \underset{x}{\mathop{\text{min}}}\, {{\left\| x \right\|}_{1}}\ \ \ \text{s}\text{.t}\text{.}\ \ \ y=\mathit{\Phi }x $

(2)

(2) 式是一个l₁范数求解最小化的线性规划（linear programming, LP）问题，其中${{\left\| x \right\|}_{1}}$是变量x的绝对值之和，对其进行求解即可恢复出原始信号.

CS-MRI就是将欠采样得到的少量K空间数据，进行某种稀疏变换，然后再结合(2)式的l₁范数凸优化问题，最终重构出图像.CS-MRI的数学模型^[8]可表示为(3)式

$ \begin{align} & \text{minimize}\ \ \ {{\left\| \mathit{\Psi }x \right\|}_{1}} \\ & \text{s}\text{.t}\text{.}\ \ \ \ \ \ \ \ \ \ \ \ {{\left\| {{F}_{S}}x-y \right\|}_{2}}<\varepsilon \\ \end{align} $

(3)

其中x是待重建的磁共振图像；Ψ表示线性稀疏变换；F_S表示欠采样傅里叶变换；y是从MRI扫描仪欠采样测量得到的K空间数据；ε控制重建数据和测量数据的保真度，代表测量数据的噪声水平^[16].利用拉格朗日乘子法将(3)式转换为无约束优化问题形式，如(4)式所示

$ x=\underset{x}{\mathop{\text{arg}\ \text{min}}}\, \left\{ \left\| {{F}_{S}}x-y \right\|_{2}^{2}+\lambda {{\left\| \mathit{\Psi }x \right\|}_{1}} \right\} $

(4)

其中λ为正则化参数，对(4)式进行求解即可重建出磁共振图像.近年来提出的各种CS-MRI方法都是基于改变(4)式的稀疏变换基或加入不同的正则项演化而来.目前，求解(4)式的CS-MRI算法多采用TV和小波稀疏约束项线性结合的方式来重建磁共振图像，如(5)式所示.

$ x=\underset{x}{\mathop{\text{arg}\ \text{min}}}\, \left\{ \frac{1}{2}\left\| {{F}_{S}}x-y \right\|_{2}^{2}+\alpha {{\left\| x \right\|}_{TV}}+\mathit{\beta }{{\left\| \mathit{\Psi }x \right\|}_{1}} \right\} $

(5)

其中α和β是两个正则参数；Ψ是小波变换；${{\left\| x \right\|}_{TV}}={{\Sigma }_{i}}{{\Sigma }_{j}}\sqrt{{{({{\nabla }_{1}}{{x}_{ij}})}^{2}}+{{({{\nabla }_{2}}{{x}_{ij}})}^{2}}}$，${{\nabla }_{1}}$和${{\nabla }_{2}}$分别为图像水平方向和垂直方向的梯度算子，i、j表示磁共振图像在x、y方向像素位置.

1.2 双树小波变换与小波树稀疏联合的CS-MRI算法 1.2.1 双树小波变换理论

稀疏变换在CS理论中扮演着至关重要的角色，CS-MRI重建的结果很大程度上取决于稀疏变换.一般来说，信号越稀疏重建所需的测量值越少.除此之外，当存在测量不足和/或噪声时，重构图像中伪影的消除也依赖于稀疏转换^[17].双树小波变换（dual-tree complex wavelet transform, DT-CWT）是离散WT的增强形式，它能通过使用复值小波和尺度函数来弥补离散WT的缺点^[6]. DT-CWT在二维或高维信号中具有移位不变性和定向选择性，能够很好地解决传统离散WT带来的在CS-MRI重建图像时产生边缘伪影、混淆现象及不具方向选择性等缺点.这使得DT-CWT成为CS-MRI中有效的稀疏变换工具.

1.2.2 小波树稀疏理论

自然数据（信号或图像）的小波系数的大部分系数值近似为0，只有少量的系数有较大值，具有很好的稀疏性.如：二维图像的小波系数服从四叉树结构，其在粗尺度上的系数被称为根节点，在细尺度上的系数被称为叶节点，每一根节点下面跟着四个细尺度上的子节点^[18]，具有很高的稀疏性.图 1所示为图像的小波四叉树结构.

小波树结构为CS-MRI提供了良好的技术优势^[19].根据凸优化理论，结构稀疏可以通过凸优化算法求解；原始信号的先验信息，如组或图结构能被利用，还可以进一步降低采样率.当前提出的一些算法已经开始利用小波系数树结构来改进CS重构，这些算法很好的利用了树结构的理论，如：树结构的母-子关系，母系数有一个大的或小的值，那么它的子系数值往往与其保持一致.以上先验知识若能被完全利用，仅仅$O(K+\text{log}n)$的观测数据就能恢复出满足树稀疏结构的原始数据，相比CS理论中标准K稀疏数据需要$O(K+K\text{log }n)$的采样数据才足够恢复出长度为n的原始数据，采样量减少了一半；当n无限大时，利用树结构的优势就更加显著了.文献[14]也验证了小波树结构能帮助提高成像速度.

因此，本文提出了一种双树小波变换、小波树和梯度稀疏联合的新的稀疏变换模型，并将其应用于LF MRI，以便在提高其成像速度的同时，改善其图像质量.该模型中的稀疏变换双树小波变换和小波树稀疏被建模为一个重叠组正则项、TV作为惩罚项用来抑制噪声.TV惩罚项和l₁范数约束项实际上是作为树形结构模型的补充，使新提出的稀疏变换模型对实际低场磁共振图像更具鲁棒性.实现以上条件的目标函数如(6)式所示：

$ x=\underset{x}{\mathop{\text{arg}\ \text{min}}}\, \left\{ \frac{1}{2}\left\| {{F}_{S}}x-y \right\|_{2}^{2}+\alpha {{\left\| x \right\|}_{TV}}+\beta \left( {{\left\| \mathit{\Phi }x \right\|}_{1}}+\sum\limits_{g\in G}{{{\left\| \mathit{\Phi }{{x}_{g}} \right\|}_{2}}} \right) \right\} $

(6)

式中Φ表示双树复小波变换系数，G表示所有母-子系数组的集合，g为G的一个组，x_g是属于母-子系数组的所有数据.

1.3 双树小波变换与小波树联合的CS-MRI重建算法求解

现有的CS-MRI模型主要是求解最小二乘保真项、TV与小波稀疏正则项线性组合的目标函数最小化的问题.最早使用的是CG下降算法求解该问题.而TVCMRI和RecPF算法分别用作一种算子分裂算法和变量分裂算法来求解该问题；FCSA算法将原始问题分解为两个简单的子问题并用快速迭代收缩阈值算法（fast iterative shrinkage threshold algorithm, FISTA）分别求解它们.上述算法都是CS-MRI最常用的算法，但是它们都没有利用树状结构这一先验知识来增强其算法性能.文献[12]提出的WaTMRI算法能求解树结构、小波稀疏和TV联合的凸优化CS-MRI模型问题.并且该算法收敛速度快，每次迭代仅花费$O(n\ \text{log}\ n)$时间；而且在实际图像应用中，优于其它CS-MRI算法和几种一般的基于树的算法.因此本文也采用WaTMRI算法来求解双树小波变换、小波树和梯度稀疏联合的CS-MRI模型，称为DT-WaTMRI算法.

DT-WaTMRI算法与先前CS-MRI的重建算法的主要不同是将稀疏变换基由小波变换改变为双树小波变换.这个算法不能用凸优化理论有效地求解.因此引入变量z约束x重叠结构，使(6)式变为非重叠凸优化问题，令$G\mathit{\Phi }x=z$，则(6)式衍变为下式：

$ x=\underset{x, z}{\mathop{\text{arg}\ \text{min}}}\, \left\{ \frac{1}{2}\left\| {{F}_{S}}x-y \right\|_{2}^{2}+\alpha {{\left\| x \right\|}_{TV}}+ \\ \beta \left( {{\left\| \mathit{\Phi }x \right\|}_{1}}+\sum\limits_{i=1}^{s}{{{\left\| {{z}_{gi}} \right\|}_{2}}} \right)+\frac{\lambda }{2}\left\| z-G\mathit{\Phi }x \right\|_{2}^{2} \right\} $

(7)

对于x子问题，我们可以联合(7)式中第一个和最后一个二次项惩罚项来求解，而z子问题则通过分组软阈值来解决，剩下部分与FCSA算法类似，可通过迭代方案有效解决.z子问题如(8)式所示：

$ {{z}_{gi}}=\underset{{{z}_{gi}}}{\mathop{\text{arg}\ \text{min}}}\, \left\{ \beta {{\left\| {{z}_{gi}} \right\|}_{2}}+\frac{\lambda }{2}\left\| {{z}_{gi}}-{{(G\mathit{\Phi }x)}_{gi}} \right\|_{2}^{2} \right\}, \, \ i=1, 2, \cdots , s $

(8)

其中gi表示第i组，s表示全部组数.z子问题有一个封闭式的软阈值解决方案如(9)式所示：

$ {{z}_{gi}}=\text{max}\left( {{\left\| ri \right\|}_{2}}-\frac{\beta }{\lambda }, 0 \right)\frac{ri}{{{\left\| ri \right\|}_{2}}}, \ \, i=1, 2, \cdots , s $

(9)

式中$ri={{(G\Phi x)}_{gi}}$.令$f(x)=\frac{1}{2}\left\| {{F}_{S}}x-y \right\|_{2}^{2}+\frac{\lambda }{2}\left\| z-G\mathit{\Phi }x \right\|_{2}^{2}$，这是一个有着莱布尼兹常数（Lipschitz, L_f）的凸平滑函数；再令${{g}_{1}}(x)=\alpha {{\left\| x \right\|}_{TV}}$，${{g}_{2}}(x)=\beta {{\left\| \mathit{\Phi }x \right\|}_{1}}$，${{g}_{1}}(x)$和${{g}_{2}}(x)$是凸的非平滑函数.为方便起见，将(9)式表示为(10)式：

$ z=shrinkgroup\left( G\mathit{\Phi }x, \frac{\beta }{\lambda } \right) $

(10)

本文提出的DT-WaTMRI算法的主要步骤为：

Input: $\rho {\rm{ = }}{\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle {{L_f}}$}},\alpha ,\beta ,\lambda ,{t^1} = 1,{r^1} = {x^0}$

For  k=1 to N do

$z = shrinkgroup(G\mathit{\Phi }{x^{k - 1}}, \beta /\lambda)$

    ${x_g} = {r^k} - \rho \nabla f({r^k})$

    ${x_1} = pro{x_\rho }(2\alpha ||x|{|_{TV}})({x_g})$

    ${x_2} = pro{x_\rho }(2\mathit{\beta }||\mathit{\Phi }x|{|_1})({x_g})$

    ${x^k} = \left({{x_1} + {x_2}} \right)/2$

    ${t^{k{\rm{ + }}1}}{\rm{ = }}\left[{{\rm{1 + }}\sqrt {1 + 4{{({t^k})}^2}} } \right]{\rm{/2}}$

    ${r^{k + 1}} = {x^k} + \left({{t^k} - 1/{t^{k + 1}}} \right)\left({{x^k} - {x^{k - 1}}} \right)$

End for

算法中${L_f} > 0$，N表示最大迭代次数，x表示待重建的图像，r迭代时存放相应值的中间变量，t表示适当的迭代步长；$\nabla f({r^k}) = {F_S}^T({F_S}{r^k} - y) + \lambda {\mathit{\Phi }^T}{G^T}({\mathit{\Phi }^T}{G^T}{r^k} - z)$，表示函数f在点${r^k}$处的梯度，其中${F_S}^T$表示部分逆傅里叶转换.对任意标量$\rho > 0$，近端映射（proximal map，prox）归一化连续凸函数$g(x)$的定义如(11)式所示：

$ pro{x_\rho }(g)(x) = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_u \left\{ {g(u) + \frac{1}{{2\rho }}{{\left\| {u - x} \right\|}^2}} \right\} $

(11)

$ project\left( {x, [l, u]} \right) = \left\{ \begin{array}{l} x, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} l \le x \le u\\ l, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \le l\\ u, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \ge u \end{array} \right. $

(12)

(12) 式将图像的像素值规范化到[l, u]范围内，且u≥l≥0，若不执行该操作，则重建图像由于存在负数像素数值而出现伪影.

2 材料和方法

为验证所提出模型的良好性能，我们将DT-WaTMRI与FCSA和WaTMRI算法的重建结果进行比较.所有的算法都在处理器为Intel(R)Core(TM)i5-4460 CPU @3.20GHz的Dell Inspiron3847的Matlab2016b环境中运行.

2.1 数据采集

实验测试数据分为三组，分别为Matlab自带体模仿真数据，以及使用MesoMR23-040H-1的LF MRI扫描仪（苏州纽迈分析仪器有限公司）采集的哈密瓜及小鼠腹部的大小为256*128的全K空间数据.MRI扫描仪的磁场强度为0.5 T，脉冲序列选用自旋回波（spin echo，SE）序列，数据采集时实验参数分别设置如下：序列重复时间（repetition time, TR）为2 000 ms、回波时间（echo time, TE）为100 ms、扫描次数（number of scan, NS）为4.根据设定参数在LF MRI扫描仪上进行回波扫描得到所需数据.

2.2 预处理工作

基于CS理论的MRI技术包括图像的稀疏表示、观测矩阵设计和算法重构三个方面.CS算法仅需采集Nyquist采样定理要求数据量的30%，甚至更少的数据重建原始信号.而磁共振图像的数据空间是K空间，K空间中存储在不同位置的磁共振信号对图像的对比度和空间分辨率的贡献有所差异：K空间的中心部分，信号强度大、空间频率低，形成图像对比度的权重大；K空间的外围部分，信号强度小、空间频率高，形成图像空间分辨率的权重大.因此CS-MRI对K空间数据欠采样选用变密度的随机测量矩阵，在K空间中心部分的采样权重多于K空间外围部分，即能加快采样速率又能保证成像质量.随机测量矩阵能很大概率的满足约束等距性质（restricted isometry property, RIP），RIP性质是CS理论完整恢复出原始信号的保证.但是随机测量矩阵硬件实现困难，所以当前实际应用中测量矩阵一般选用如图 2所示（白色区域表示采样，黑色区域表示未采样）的伪随机确定性测量矩阵.由K空间的冗余特性采用图 2所示的任意一种变密度的采样模式，都能缩短采样时间、提高成像速度、减少运动伪影.但是，不同采样脉冲序列产生不同的K空间模型，而不同的K空间模型需选用不同的随机变密度采样模式.由于实验所用采样脉冲序列为多回波SE序列，得到的K空间为笛卡尔单向平行轨迹，因此采用图 2(a)所示的一维变密度采样方式.

2.3 实验及评价指标

我们定义随机变密度欠采样矩阵的欠采样率为m/n（m和n分别为K空间的行、列大小），为方便起见，实验中所有算法的欠采样率均设为30%.利用DT-WaTMRI算法求解重建目标函数时，使其中的稀疏约束采用两级分解的Daubechies小波，正则化参数λ、β、γ分别设为0.1、0.1、0.025，最大迭代次数设置为200.

在评价重建图像质量时，信噪比是一项重要的技术指标，但信噪比高的图像还不能确保将两个图像相邻区域或结构有效地区分开来，图像还必须在特定的组织和周围组织间有足够的对比度.因此，本文对实验中不同算法的重建性能评估采用均方根误差（root mean square error, RMSE）、峰值信噪比（peak signal to noise ratio, PSNR）、方差（variance, Var）等指标来衡量.其公式分别如(13)、(14)、(15)式所示.

$ RMSE = \sqrt {\frac{1}{{MN}}\sum\nolimits_{i = 0}^M {\sum\nolimits_{j = 0}^N {{{[x(i, j) - \hat x(i, j)]}^2}} } } $

(13)

$ PSNR = 20{\rm{lg}}\frac{{{x_{{\rm{max}}}}}}{{RMSE}} $

(14)

$ Var=\frac{\sum\nolimits_{i=0}^{M}{\sum\nolimits_{j=0}^{N}{{{[\hat{x}(i,j)-\bar{\hat{x}}(i,j)]}^{2}}}}}{\sum\nolimits_{i=0}^{M}{\sum\nolimits_{j=0}^{N}{\hat{x}(i,j)}}} $

(15)

公式(13)、(14)、(15)中x(i, j)表示全采样重建图像，$\hat x(i, j)$表示欠采样重建图像，$\bar{\hat{x}}(i,j)$表示欠采样重建图像平均灰度值，M、N分别代表磁共振图像在x、y方向像素点个数，i、j分别代表磁共振图像在x、y方向像素点位置.

3 结果和讨论

图 3~5所示为不同数据的重建结果，图 3(a)、4(a)、5(a)分别为不同样品的全采样重建图像、变密度采样矩阵、欠采样零填充重建图像；其中全采样图像作为欠采样重建图像的评价标准，变密度采样矩阵表示变密度欠采样K空间数据的采样方式，欠采样零填充重建图像是对欠采样的K空间数据不使用任何算法直接进行傅里叶变换得到的重建图像，零填充重建图像也作为不同重建算法欠采样重建图像的对比.图 3(b)、4(b)、5(b)分别为使用FCSA、WaTMRI、DT-WaTMRI算法对体模数据、哈密瓜横断面数据、小鼠腹部横断面图像数据的重建结果. 图 5中小鼠腹部横断面数据重建图像太小主要是因为对整个小鼠样本进行扫描时，受限于LF MRI扫描仪放置样品的口径，在前期处理中设定了成像视野（field of view, FOV），后续不易改动，但是不影响图像分辨率和后续的说明与结论.图 6所示为使用FCSA、WaTMRI、DT-WaTMRI三种算法对体模数据、哈密瓜数据、小鼠腹部横断面数据的的重建结果的客观评价；图 6(a)~(c)是对体模数据的重建结果的客观评价，图 6(d)~(f)是对哈密瓜横断面数据的重建结果的客观评价，图 6(g)~(i)是对小鼠腹部横断面图像数据的重建结果的客观评价；图 6(a)、(d)、(g)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像与全采样重建图像的RMSE随迭代次数的变化情况，(b)、(e)、(h)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像与全采样重建图像的PSNR随迭代次数的变化情况，(c)、(f)、(i)为使用FCSA、WaTMRI、DT-WaTMRI三种算法欠采样重建图像的Var随迭代次数的变化情况.

对CS-MRI来说重建时间也是评价算法性能的一项重要指标，因此本文也对不同算法的重建时间做了评价.表 1列举了不同算法对不同样本的重建时间（单位为s）.

图 3、4、5所示为不同样本磁共振图像的重建结果，图 3(a)、4(a)、5(a)中的零填充重建图像是不同样本的全采样K空间数据经一维变密度随机欠采样后得到欠采样的K空间数据，对欠采样的K空间中缺失的数据简单填零后直接傅里叶变换的重建结果，由结果可知零填充重建所得到的重建图像伪影比较严重.图 3(b)、4(b)、5(b)是利用不同的CS-MRI重建算法对不同样本的欠采样数据进行重建的结果；由重建结果可知本文提出的DT-WaTMRI算法与FCSA算法和WaTMRI算法相比并没有显著改善图像的重建质量，而且使用本文算法对不同样本数据的重建结果还出现了类似伪影的竖条纹，尤其在体模数据中表现明显，这可能是因为体模数据是不含任何噪声的仿真数据，当用随机变密度矩阵对其进行欠采样，相当于原始无噪声数据矩阵与变密度矩阵相乘，造成了类似于条状伪影的数据缺失；也可能是本文算法还有需要改善的地方.但是本文算法与FCSA算法和WaTMRI算法相比，确实提高了重建图像的对比度和亮度，由图 3也可以看出本文算法重建出图像的清晰度程度实际上是与原图最接近的，这点可由重建图像的评价指标Var验证.

对重建出的结果图像可以通过视觉观察来判断其优劣，然而视觉观察往往带有选择性和主观性.为了提高说服力，图 6采用客观评价指标如RMSE、PSNR、Var等，来判断不同算法的重建结果的优劣.RMSE是用来计算全采样重建图像和欠采样重建图像的均方值，然后通过均方值的大小来确定重建图像的失真程度.对不含噪声的体模数据来说，DT-WaTMRI算法的RMSE小于FCSA和WaTMRI算法；对含有噪声的哈密瓜和小鼠腹部低场磁共振图像来说，DT-WaTMRI算法的RMSE略微小于FCSA和WaTMRI.PSNR是评价画质的客观测量法，值越大图像质量越好.对不含噪声的体模数据来说，DT-WaTMRI算法的PSNR大于FCSA和WaTMRI算法；对含有噪声的哈密瓜和小鼠腹部低场磁共振图像来说，DT-WaTMRI算法的PSNR略大于FCSA和WaTMRI算法.Var反应重建图像高频部分的大小，图像的对比度越大，方差就越大.对不含噪声的体模数据来说，DT-WaTMRI算法重建图像的Var大于FCSA和WaTMRI算法；对含有噪声的哈密瓜低场磁共振图像来说，DT-WaTMRI算法重建图像的Var值与FCSA算法相近；对含有噪声的小鼠腹部低场磁共振图像来说，DT-WaTMRI算法重建图像的Var稍大于WaTMRI算法.

表 1中对所有实验样本来说，稀疏变换采用传统小波变换的FCSA算法重建出磁共振图像所用时间比稀疏变换采用传统小波变换与结构稀疏联合的WaTMRI算法少，这是由于算法中引入了小波树结构增加了计算冗余性.而稀疏变换采用双树小波与小波树结构联合的DT-WaTMRI算法重建出磁共振图像所用时间比WaTMRI算法少，在保证重建质量的同时提高了成像效率.因此结合理论和实验结果可知本文提出的DT-WaTMRI算法改进了现有的重建算法，在一定程度上提高了磁共振图像的重建质量.

4 结论

在本文中，我们提出了一种基于双树小波变换与小波树结构联合的DT-WaTMRI重建算法，并用WaTMRI算法加上一定的数据一致性约束对其进行求解.实验结果表明，该算法与其它算法在PSNR、RMSE和Var三个评价指标上，具有较好的性能；且与WaTMRI树稀疏算法相比，该算法重建计算消耗时间更少，适用于低场磁共振图像重建.