非笛卡尔并行磁共振成像数据的自适应约束重建新算法

周爱珍 梅颖洁 王 聪 莫纪江 陈武凡 冯衍秋

（南方医科大学生物医学工程学院，广州 510515）

引言

并行磁共振成像（PMRI）概念的提出，突破了传统MR成像时间受到磁场梯度以及射频场硬件性能约束的限制。多线圈部分采集方式是：使用多个射频接收线圈同时接收感应信号，利用各个线圈的空间敏感度差异来编码空间信息，降低成像所需的梯度编码步数，提高数据获取速度。现有的并行磁共振成像方法主要有3类［1，3］：第一类是基于图像域的成像方法，如敏感度编码方法（SENSE）和基于局部化敏感度分布的并行成像方法（PILS）；第二类是基于k空间域的成像方法，运用比较成熟的是GRAPPA方法，另外还有空间谐波并行采集技术（SMASH），以及在该算法上改进出的AUTO-SMASH和VD-AUTO-SMASH方法；第三类是同时基于k空间域和图像域的重建方法，如SPACE-RIP。在已知线圈敏感度的情况下，SENSE重建方法可以得到相对较好的图像。在SENSE中，重建问题是解一系列的线性方程［2］，即

式中，E为编码矩阵，；sγ为第γ个线圈的线圈敏感度；kκ为k空间第κ个采样位置；rρ为第ρ个体素的位置；f为待求的全视野图像；d为由各个线圈采集的部分k空间数据形成的矢量。

通常采用最小化残差的L2范数的平方方法求最优解。但是，方程的解非常依赖于线圈形状和采样轨迹，特别是当欠采样因子比较大时，方程往往是病态的，重建结果的信噪比（SNR）会降低。为了提高图像信噪比，通常在重建方程中添加惩罚项。TV约束和Tikhonov约束是两种常用的约束方式［4］，PM模型也是较好的约束方式。

SENSE重建TV约束模型为freg=

SENSE重建Tikhonov约束模型［2］为freg=，其中，ξ为正则化参数。和TV约束方式不同，L2范数惩罚项值的大小和作用域的整个趋势都有关。当图像存在“跳变”时，L2范数值会很大，约束最小化会使“跳变”被约束，这样能很好地去除平滑区域内的“跳变”噪声点，但也同时会导致图像边缘信息被模糊。

综上可知，TV模型能保护边缘信息，但会产生“阶梯效应”；Tikhonov约束可以在较大程度上平滑噪声及避免“阶梯效应”，但是会导致图像边缘模糊［16］。为了获得更好的约束方式，考虑将这两种平滑性度量结合起来，形成某种“混合”测度。笔者提出，由先验图像的梯度特征决定惩罚项，对重建方程约束方式进行自适应调整，并借鉴PM模型的思想，进一步调整约束方式，以便在梯度值较大的边缘区惩罚项使用TV约束方式，在梯度值较小的平滑区使用Tikhonov约束方式。这样，可以在保留图像细节的同时抑制噪声，增加图像信噪比。

1 自适应约束SENSE重建理论

对采样数据进行网格化处理，然后利用平方和重建方法［18］进行图像重建，将得到的结果作为先验图像。根据先验图像的梯度特征决定惩罚项［14，17］，从而对多通道非笛卡尔轨迹采样数据进行SENSE重建，重建模型为：

式中，p（x′）=）是关于先验图像m0的梯度幅度下降函数，且，。也就是说，在梯度较大的区域p取1，在梯度较小的区域p取2。

这里，定义p=2-M，且M定义为［5］

式中，m0为先验图像，G为高斯卷积核，λ为正则化参数，x′为像素位置，为高斯平滑后的初始图像的梯度模值的直方图中元素个数最多的值，*代表卷积运算。

从式（3）可以推出，在梯度值较小区域，表达式分母中的分子项值较小，M接近于0；在梯度值较大区域，表达式分母中分子项值较大，M接近于1。根据前述基于L1范数和L2范数正则化重建方法的优缺点，为了得到较好的结果，应该在梯度较大的边缘区域采用TV约束方式（2），在平滑区域采用Tikhonov约束方式。模型在梯度较大区域p=1，在平滑区域p=2，符合要求。

另外，由于梯度模值较小的部分多是较平滑的区域，梯度模值较大的部分往往是边缘区，如果在这两种区域里也采用自适应约束，则扩散速度较慢。笔者借鉴PM模型的思想，进一步调整约束方式，根据图像的不同设定两个边界值，对于梯度模值小于较小边界值的部分直接采用Tikhonov约束方式，梯度模值大于较大边界值的部分直接采用TV约束方式，其余部分采用自适应约束方式，即

式中：α1、α2为设定的两个边界值，可以根据图像效果选择；p的值可以在迭代过程中根据新得到的中间结果更新，也可以由初始图像给出后不再更新，比较简单的方法是中间过程不再进行p的更新。

2 方法

对于多通道笛卡尔坐标系采样数据，SENSE可以对混叠图像逐点处理，求其对应系统矩阵的伪逆，得出原始图像中发生混叠的像素值。但是，对于非笛卡尔坐标系采样数据，由于无法直接估计出混叠像素点之间的相互作用，所以不能使用上述方法求解。笔者采用非线性共轭梯度下降方法迭代解方程［6，13］，流程如图1所示。

图1 重建过程流程Fig.1 Schematic diagram of reconstruction process

此过程即将敏感度加权后的图像首先进行傅里叶变换，然后按非笛卡尔采样轨迹进行k空间数据重采样，见图2（a）。另外，有

此过程先将数据网格化至笛卡尔坐标系，然后进行傅里叶逆变换至图像域得到混叠图像，再对混叠图像进行敏感度共轭加权求和，见图2（b）。所以，E*并不是简单的E的逆［7，17］。

在上述两个过程中，涉及网格化和逆网格化。网格化即利用卷积核将非笛卡尔轨迹数据插值到笛卡尔坐标系，逆网格化过程即将笛卡尔轨迹数据插值到非笛卡尔采样轨迹。卷积核通常采用Kaiser-Bessel窗［8-12］，有

图2 图像域与k空间域数据变换过程。（a）图像域数据至k空间域非笛卡尔轨迹数据的变换过程；（b）k空间域非笛卡尔轨迹数据至图像域数据的变换过程Fig.2 Transformation of data between image domain and k-space.（a）transformation of data from image domain to non-Cartesian trajectory k-space；（b）transformation of data from non-Cartesian trajectoryk-space to image domain

另外，由于受k空间采样模式的影响，得到的图像往往会变模糊，如对螺旋采样轨迹数据来说会产生螺旋环状边叶。所以，要将得到的图像进行去极化，即除以卷积核的傅里叶变换。这样做尽管会使图像产生一些阴影，但是它能很好地抑制由于数据和采样函数卷积产生的边叶。去极化核近似为

卷积核的宽度和采样密度相关，比较宽的核在成像中往往有平滑信息的作用，这样也可以平滑掉一部分噪声，比较窄的核有锐化图像的作用，可以更好地保留小的细节信息，但是不利于消除噪声。在本实验中，取L=4，α=2。

3 实验结果与分析

实验中采用的数据来自2010 ISMRM会议，为人体动静脉畸形瘤动脉注射的X线照射结果。数据采样轨迹为可变密度螺旋轨迹，数据采集方式为8通道并行采集，每个通道采集200条螺旋，每条螺旋上采集2 000个数据点，线性欠采样螺旋按golden angle角度旋转，在13个TR后，在k空间边缘的欠采样因子近似为15。各个线圈的敏感度值由轴向层计算得来，是已经提供的数据。在每个通道采集的数据中，都加入了独立噪声。图像大小设定为512像素×512像素，本实验采用了其中的50个（序号61～110）个螺旋作为采样数据，欠采样倍数2.6。分别用平方和（sum-of-squares，SOS）重建方法、无约束SENSE重建方法、TV约束SENSE重建方法，以及所提出的自适应约束方法进行了重建，结果如图3～图5所示。

从实验结果中可以得出，同平方和重建方法、无约束SENSE重建方法相比，约束重建方法得到的结果边缘更加清晰，与传统的约束重建方法相比，本研究提出的算法同时保持了去噪及保护细节的功能。图4和图5显示，无约束重建结果中比较模糊的小细节在本重建算法中得到保护，平滑区域的噪声也得到抑制。但是，也可以从图中看到，部分小细节模糊，这主要是由于惩罚项约束方式p值的大小受先验图像（见图3（a））梯度的影响。如果先验图像中某些部位图像质量过差，那么相应部位的p值会受到影响，从而可能导致相应部位重建图像质量降低。在本实验中，先验图像在图4所示部位含有的噪声较多，而在图5所示部位的含噪声较少，p值受到先验图像中噪声影响，从而图4（d）中的部分小细节恢复得不是很好，但是图5（d）中的小细节普遍恢复得很好。

图3 4种方法重建结果。（a）平方和重建结果；（b）无约束SENSE重建结果；（c）TV约束SENSE重建结果；（d）自适应约束SENSE重建结果Fig.3 Reconstruction results of four methods.（a）the result of SOS method；（b）the result of SENSE；（c）the result of SENSE with TV constraint；（d）the result of SENSE with adaptive constraint

图4 图3（a）中实线矩形框内对应部分的放大图。（a）平方和重建结果部分放大；（b）无约束SENSE重建结果部分放大；（c）TV约束SENSE重建结果部分放大；（d）自适应约束SENSE重建结果部分放大Fig.4 Enlargement of solid rectangular ROI in Fig.3（a）.（a）the enlargement of SOS method；（b）the enlargement of SENSE；（c）the enlargement of SENSE with TV constraint；（d）the enlargement of SENSE with adaptive constraint

图5 图3（a）中虚线矩形框内对应部分的放大图。（a）平方和重建结果部分放大；（b）无约束SENSE重建结果部分放大；（c）TV约束SENSE重建结果部分放大；（d）自适应约束SENSE重建结果部分放大Fig.5 Enlargement of dotted rectangular ROI in Fig.3（a）.（a）the enlargement of SOS method；（b）the enlargement of SENSE；（c）the enlargement of SENSE with TV constraint；（d）the enlargement of SENSE with adaptive constraint

4 结论

在多通道磁共振并行成像中，非笛卡尔轨迹部分采样k空间数据重建仍是一个难点，在重建中需要考虑噪声和伪影的去除以及计算量问题。傅里叶变换和部分采样多通道非笛卡尔轨迹并行成像技术能大幅度地减少梯度编码步数，提高成像扫描速度，相比单线圈成像技术能更加有效地缩短数据获取时间，从而达到优化成像的目的。但是，由于部分采样数据的不完整性及敏感度信息获取的不准确性，导致单纯的并行重建出现伪影和噪声。笔者在本文中提出利用先验图像的梯度特征决定惩罚方式的自适应约束方法，第一部分介绍了常用的多通道非笛卡尔采样部分数据的各种SENSE约束重建方法，第二部分综合各种约束重建方式的特点提出了自适应算法。

实验结果显示，相比原来的方法，所提出的方法可以有效地抑制噪声、保留细节。但是，由于重建方法需要多次迭代完成，所以重建时间上有一定的延长。如何进一步的优化算法，是今后研究的一个重点。

［1］Liu B，King K.Regularized sensitivity encoding（SENSE）reconstruction using bregman iterations［J］.Magnetic Resonance in Medicine，2009，61（1）：145-152.

［2］Pruessmann KP，Weiger.M.SENSE：Sensitivity Encoding for fast MRI［J］.Magnetic Resonance in Medicine，1999，42（5）：952-962.

［3］陈武凡.并行磁共振成像的回顾、现状与发展前景［J］.中国生物医学工程学报，2005，24（6）：649-654.

［4］Vijayakumar S，Duensing R，Huang F.G-factor and gradient weighted denoising with edge restoration（g-denoiser）for SENSE reconstruction MR images［C］//Christophe J，Marin O，eds.2008 5th IEEE International Symposium on Biomedical Imaging：from Nano to Macro.Paris：IEEE，2008：472-475.

［5］Perona P，Malik J.Scale space and edge detection using anisotropic diffusion［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12（7）：629-639.

［6］Lustig M，Donoho D，Pauly JM.Sparse MRI：the application of compressed sensing for rapid MR imaging［J］.Magnetic Resonance in Medicine，2007，58（6）：1182-1195.

［7］Pruessmann KP，Weiger M.Advances in sensitivity encoding with arbitrary k-space trajectories［J］.Magnetic Resonance in Medicine，2001，46（4）：638-651.

［8］Rasche V，Proksa R，Sinkus R，et al.Resampling of data between arbitrary grids using convolution interpolation［J］.IEEE Trans Med Imaging，1999，18（5）：385-392.

［9］Jackson JI，Meyer CH，Nishimura DG.Selection of a convolution function for fourier inversion using gridding［J］.IEEE Trans Med Imaging，1991，10（3）：473-478.

［10］Osulllivan J.A fast sinc function gridding algorithm for fourier inversion in computed tomography［J］.IEEE Trans Med Imaging，1985，4（4）：200-207.

［11］Çukur T，Santos JM，Nishimura DG，et al.Varying kernel-extent gridding reconstruction for under sampled variable-density spirals［J］.Magnetic Resonance in Medicine，2008，59（1）：196-201.

［12］Hoge RD，Kwan RKS，Pike GB.Density compensation functions for spiral MRI［J］.MRM，1997，38（1）：117-128.

［13］Beatty PJ，Nishimura DG.Rapid gridding reconstruction with aminimal oversampling ratio［J］.IEEE Trans Med Imaging，2005，24（6）：799-808.

［14］Osher S，Burger M，Goldfarb D.An iterative regularization method for total variation-based image restoration［J］.Multiscale Model Simul，2005，4（2）：460-489.

［15］Huang F，Chen YM，Duensing GR.Application of partial differential equation-based in painting on sensitivity maps［J］.Magnetic Resonance in Medicine，2005，53（2）：388-397.

［16］Rudin LI，Osher S，Fatemi E.Nonlinear total variation based noise removal algorithms［J］.Physica，1992，60（1-4）：259-268.

［17］Block KT，Uecker M，Frahm J.Under sampled radial MRI with multiple coils.Iterative image reconstruction using a total variation constraint［J］.Magnetic Resonance in Medicine，2007，57（6）：1086-1098.

［18］Larsson EG，Erdogmus D，Yan R，et al.SNR-optimality of sum of-squares reconstruction for phased-array magnetic resonance imaging［J］.Journal of Magnetic Resonance，2003，163（1）：121-123.