基于自适应字典学习的动态磁共振并行重建

王 悦 蒋慧敏 汪 洋,2*

1(台州学院 浙江 台州 318000)2(河海大学 江苏 南京 210098)

0 引 言

磁共振成像(MRI)技术具有无创伤、无辐射、分辨率高和可多维成像等优点，被广泛应用于临床医学各个系统，是继CT以后的又一重要临床检测方法。但MR成像速度慢是其一大缺点，尤其是动态磁共振成像(dynamic magnetic resonance imaging，dMRI)，需要在较短时间内获得高时空分辨率的MRI图像序列，目前是医学界的一个难题[1]。过长的扫描时间加上病人的器官运动(如呼吸、吞咽等)，会导致成像模糊，同时也无法满足动态实时成像和功能成像的需求。在k空间进行降采样是提高成像速度的一种方法，但如果直接从k空间逆傅里叶变换重建图像，根据奈奎斯特采样定理，会导致重建图像产生混叠效应。动态磁共振成像数据在时空域具有很强的稀疏特性，使得压缩感知(compressive sensing，CS)技术被广泛应用于MR图像重建。CS理论指出，如果一个信号在变换域是稀疏的，且变换基和测量矩阵是不相关的(又称有限等距性质，RIP)，则采用CS方法可以从降采样(远小于奈奎斯特采样率条件下)的数据样本当中，通过非线性重建算法完美重建该信号。

dMRI重建方法可以分为在线模式和离线模式。采用离线模式重建时，需要在重建之前获得整个dMRI序列的采样数据。常见的离线重建方法有运动校正[2-3]、字典学习[4-6]、低秩近似[7-8]等。这些方法充分利用整个dMRI序列稀疏特性进行高精度重建，其缺点是重建速度较慢，并且需要等待较长的扫描时间。采用在线方法重建时，每一帧的重建仅仅跟之前的帧有关，可以边扫描边重建，节省了等待时间，但由于缺乏整个序列的完整信息，重建精度无法保障，同时对重建速度也提出了更高的要求。

在线重建有两种常用的方案：串行和并行。串行方案通常利用图像或变换域中相邻帧之间的稀疏性，这也是大多数现有在线方法中常用的策略。然而，这些方法往往会导致累积误差。Chen等[9]采用了一种新的动态全变差(dTV)的并行重建方案来解决误差累积问题。但该方法每次只能利用两帧之间的稀疏作为重建的先验知识，导致比离线方法精度低。

本文提出了一种基于并行自适应字典学习的实时动态磁共振重建方法，将原本在离线模式下运行的相对较慢的字典学习算法应用到在线模式当中来。如图1所示，以高精度采样的第一帧为参考，实现对任意n个相邻帧MR图像的实时在线重建。以三维图像小块作为重建对象，采用K-SVD和batch-OMP两种算法优化重建速度。实验结果表明，该算法在重建精度和速度上都具有一定优势。

图1 并行自适应字典学习dMRI重建原理(以3帧为例)

1 压缩感知MRI理论

MRI采集k域(傅里叶变换域)信号，再通过逆傅里叶变换重建空域信号。傅里叶变换是一种线性变换，要求k域信号数必须等于图像域的像素数。2007年Lustig等[10]提出了压缩感知磁共振成像(CSMRI)的概念，以远低于奈奎斯特频率的采样频率对k空间信号进行采集，从而极大降低成像时间，并通过非线性重建算法完美重建原始图像。

给定一个dMRI图像序列表示为x，k-t空间的欠采数据为y，压缩感知dMRI重建问题可以归结为如下l0范数最小化问题。

(1)

CSMRI需要满足3个条件：

(1) 稀疏矩阵Ψ，使得原始MR图像在稀疏域中仅有少量非零系数；

(2) 矩阵FuΨ-1满足RIP条件，一般采用随机或伪随机采样矩阵，降低采样数据的相关性；

(3) 非线性重建算法，能快速准确重建原始图像。

2 自适应字典学习并行重建模型

2.1 算法总体架构

本文设计一套基于并行自适应字典学习的实时动态磁共振重建方法，算法流程如图2所示。

图2 流程框图

2.2 字典学习和稀疏编码

过完备字典在信号的稀疏表示当中有着广泛的应用。各种基于图像库训练的非自适应图像重建方法因其适应性较差、重建效果一般等缺点逐渐被基于自适应字典学习的图像重建算法所取代。这里的自适应是指由压缩感知系统采集到的测量值作为图像观测样本数据来训练字典。由于字典学习的计算量会随着字典尺寸的增加而快速增加，我们通常将待处理的图像分为重叠的小块，来代替整幅图像进行稀疏表示。字典学习问题可描述为：

(2)

式中：D为过完备字典；Ri为重叠取块算子；αi为x的稀疏表示系数；Γ={α1,α2，…,αI}为稀疏稀疏αi的集合；T0为稀疏度阈值常数。式(2)中第一项和稀疏约束条件确保了过完备字典对每个图像块的最优稀疏逼近。第二项为数据保真项，v是一个与k空间采样时叠加的高斯白噪声的标准差σ有关(v=λ/σ)的参数。

由于dMRI动态扫描时，每一帧的主要器官相对位置近似不变，我们选取序列的第一帧图像为高精度采样(一般取大于50%的采样率)，为后续帧的重建提供高精度参考，这可由磁共振机的一次预扫描完成。第一帧与其余任意相邻的n-1帧图像组成一系列n帧的子序列图像Xs(j)，以3帧为例，如图1所示，Xs(j)=[x1,x2j,x2j+1]，j=1,2,…，x2j表示原dMRI序列x中的第2j帧图像。只要相邻两帧采样完毕，即可进行实时并行重建，无需等待后续帧采样完成。每个子序列的重建与其他子序列不相关。

式(2)可通过以下过程进行优化求解：

字典学习过程，固定子序列图像Xs，将问题转化为求解式(3)中D和ɑ最优解的子问题。

(3)

(4)

(5)

算法1Batch-OMP算法

1: 输入：α0=DTx,0=xTx,G=DTD,重建误差

2: 输出：稀疏系数γ满足x≈Dγ

3: Init: SetI:=(),L1:=[1],δ0:=0,n:=1

4:whilen-1>do

7:ifn>1then

8:w:=Solve forw{Ln-1w=gIn}

10:endif

13:βn=GIncn

14:αn:=α0-βn

17:n:=n+1

18:endwhile

表1 算法1参数说明

算法1中，上标n为迭代次数，D为初始化DCT字典，x为字典学习图像小块对应的列向量，α为稀疏系数。这里G在每次迭代当中保持不变，可以在迭代之前预先计算好，以节省重建时间。在每次迭代过程当中只更新α当中的一列(|αk|取值最大的那一列)。算法1采用丘拉斯基分解方法(8-9)，将正定矩阵GIn,In分解为下三角矩阵L和其转置的乘积的分解，减少了矩阵求逆过程的迭代运算量。最后采用误差约束模型(15-16)循环更新D和α。

2.3 图像重构

图像重构过程，固定过完备字典和稀疏系数，将问题转化为求解式(5)中待重建子序列Xs最优解的子问题。

(6)

式(6)是一个普通的最小二乘法问题，对唯一的变量xs求导并令其等于0得到：

(7)

直接对式(7)进行求解，计算量相当大，因为要对一个P×P的矩阵进行求逆才可得到重建结果，这使得求解此问题的时间复杂度达到O(p3)。我们可将运算由图像空间变换到傅里叶空间进行，对式(7)两边进行傅里叶变换得到：

(8)

(9)

3 实 验

所有的重建实验都在一台带有3.0 GHz Intel i7 CPU的笔记本电脑上使用MATLAB 2013a进行仿真。选取一组心肌灌注dMRI序列(分辨率为192×192×30帧)用于评估算法的性能。伪射线采样矩阵用来对测试图像在k-t空间进行欠采样，模拟磁共振机加速成像。如图3所示，第一帧采样率为50%，其他帧采样率为15%。

(a) 第一帧采样矩阵(b) 其他帧采样矩阵图3 采样矩阵

以3帧构成一个子序列为例，进行并行重建。图像质量的客观评价方法采用峰值信噪比(Peak Signal to Noise Rate,PSNR)，单位为dB。PSNR计算公式如下：

PSNR(x,y)=20lg(Imax/RMSE(x,y))

(10)

式中：x和y分别为原始图像和重建图像；Imax为图像像素最大值，一般为了减少矩阵运算量，要对图像进行归一化操作，则Imax=1。RMSE为两幅图像的均方根误差。

实验参数设置：为了快速并行重建dMRI序列，本文算法将三维图像块大小设置为最小的3×3×3，经过实验测试，重建图像质量与图像块尺寸关系不大，但是重建速度会随着图像块尺寸的增加而迅速下降。字典大小设为27×108，每次在线学习的原子个数为5 000。迭代次数默认设置为20次，迭代结束条件设置为两次迭代的PSNR值误差在0.02 dB以内。对比的算法有离线模式的k-t SLR[8]和在线模式的TV[9]、DTV[9]三种。三个对比实验算法的参数均按其最优参数设置。

图4给出了本文算法重建一个dMRI子序列(以3帧为例)的过程。经过14次迭代，耗时2.8秒，就将k空间欠采样后零填充的低信噪比图像(31.1 dB)，重建为高信噪比的图像(36.3 dB)。图5给出了重建前后的图像视觉对比，可见，重建后的图像已经完全去除了欠采样造成的伪影。

图4 迭代次数与PSNR关系

(a) 零填充图像 (b) 重建后图像图5 重建前后图像对比

图6对比了无采样噪声条件下各种算法的重建效果，以列为单位，第一行为各自的重建图像，第二行为重建图像与原始图像的误差图，为凸显效果，我们将误差放大5倍进行比较。显然本文算法的误差最小，尤其是关键的心脏部位，视觉效果最为逼真。

图6 各算法重建MR图像和误差图像比较

在有采样噪声情况下，假设k空间叠加的高斯白噪声标准差为σ，式(9)中λ=q/σ的取值就与重建效果息息相关。图7给出了参数q与重建图像的RMSE之间的关系。采样信噪比分别取23 dB(强噪声环境，σ=0.05)、37 dB(中等噪声环境，σ=0.01)和57 dB(低噪声环境，σ=0.001)。可见，q=0.005是考虑了在不同的采样噪声环境下的一个最优选择。尤其是在强噪声环境下，q=0.005时重建效果最佳。

图7 重建均方根误差与q的关系

本文算法适用于任意相邻的n帧子序列构成的并行重建，表2给出了子序列帧数n与重建PSNR值和算法重建时间之间的关系(20%的采样率条件下)。可见，以3帧作为一个子序列重建单元重建的精度(PSNR值)是最高的，而重建的时间又是最少的。并且各种参数设置下，重建精度都在可控的误差范围内，体现了本文算法的稳定性。

表2 重建子序列帧数与重建效果之间关系

本文算法与三个对比算法的重建PSNR值对比结果如表3所示。在不同的采样率下，本文算法均是最佳的。尤其在15%的低采样率条件下，相比其他算法有0.4 dB以上的性能提升。以3帧为一个dMRI子序列为例，本文算法平均每1.4 s即可重建一帧图像，重建速度远远快于离线模式(k-t SLR)，但慢于两种采用TV方法的算法。

表3 不同采样率下重建图像的PSNR值

4 结 语

快速成像一直是MRI创新的重点，从最早期扫描一个序列要花费一个多小时，到如今最新的压缩感知技术已然可以实现10倍以上加速的超快速成像。本文使用一种基于自适应字典学习的动态磁共振并行重建算法，实现了对dMRI图像的实时高精度重建。与传统的离线模式k-t SLR算法和在线模式DTV算法相比，本文算法在不同的采样率下均有0.2 dB以上的性能提升，为实时在线重建dMRI提供了一种解决方案，在磁共振成像领域具有广阔的应用前景。