1　引言

       磁共振成像(magnetic resonance imaging，MRI)作为具有前景的影像手段，已成为现代临床医学的"眼睛"，但存在扫描速度慢的瓶颈问题，限制了其在危重症疾病中的应用，如为避免运动伪影往往需要屏住呼吸，而这对重症患者来说几乎不可能完成。压缩感知(compressed sensing，CS)^[1,2]作为一种从源头上减少采样数据的新理论，被视为最具前景的快速MRI方法。2007年，随着压缩感知—磁共振成像(compressed sensing MRI，CS-MRI)的提出^[3]，CS-MRI因其在加快MRI扫描速度方面显露出巨大的潜能而备受关注，并迅速成为一个研究热点^[4,5]。然而，目前在CS-MRI临床试验中所获得的加速比远没有充分挖掘CS理论潜在的大幅降低采样数据的能力^[6]。究其原因，一方面，CS自身理论体系仍处在一个起步阶段，还存在许多有待解决的问题^[7,8]；另一方面，由于受制于MRI的物理约束，导致CS-MRI的感知矩阵(sensing matrix)与稀疏矩阵(sparse matrix)^[9,10]之间存在较强的相干性，犹如一道"屏障"横亘在CS与MRI之间，制约了CS的性能。这种因感知矩阵和稀疏矩阵的强相干性导致CS采样能力下降(即使信号本身非常稀疏)的现象被称为"相干性屏障"(coherence barrier)^[11]。为了解决"相干性屏障"问题，许多研究者从CS-MRI的编码模型和不相干采样轨迹设计方面提出了许多解决办法，如将具有普适不相干的随机编码替代傅里叶编码^[12]；利用Toeplitz随机矩阵来替代傅里叶编码^[13,14]；对成像对象预先调制的傅里叶编码方式^[15]；基于Noiselet的编码方式^[16]等。

       然而，在CS-MRI不相干评价机制方面，大多数情况下还是沿用CS理论的不相干评价指标，如零空间、约束等距(restricted isometry property，RIP)、不相干性(incoherence)等指标^[17]，这些指标是评价CS重建信号质量好坏的关键指标。对于CS-MRI而言，Lustig等^[3]首次提出了与CS Gram矩阵（测量矩阵与其转置的内积)对应的变换点扩散函数(transform point spread function，TPSF)概念，并将主瓣对旁瓣的比(类似于CS Gram矩阵的次大值)作为衡量不相干性的指标。文献[11,18,19,20]等相继指出，这种基于次大值的不相干度量指标是建立在测量矩阵在最坏情况的假设之上，因此，该指标是一个非常"紧"的指标，容易导致实际性能明显好于理论预测的结果。为此，线相干(line coherence)、块相干(block coherence)等渐进不相干(asymptotic incoherence)指标相继提出^{[21,22,23,24]}。这些新的不相干指标虽然较原来的指标"松弛"，但仍然是基于次大值的一维指标，忽略了其分布特征。此外，这些不相干性指标没有表示成MRI采样轨迹的显示函数，其作用更多地体现在"事后检验"测量矩阵的不相干程度上，很难形成指导设计不相干采样轨迹的方法。

       本文首先将CS-MRI的不相干性表示成TPSF矩阵形式，同时，通过给定稀疏矩阵为单位阵的情况下将TPSF表示成更能直接反映MRI扫描轨迹不相干特性的点扩散函数(point spread function，PSF)。通过分析PSF的数学定义给出了不同采样模式下不相干性对PSF形状的影响。

2　磁共振成像的数学模型

       磁共振重建图像与原始图像实际上可以表示成如下的卷积形式：

       式中为采样函数，它由一系列幅度为1的冲击函数组成，采样的位置由决定。也就是说轨迹的PSF函数定义为采样轨迹的傅里叶变换。磁共振成像是在有限的视野内获取信息，而且由于采样速度的因素也不能采集无限个数据点，因此，需要对磁共振成像的数学模型公式进行离散化得到：

       同样对PSF进行离散得到：

       式中。当采样间隔相等时，为常数；当采样间隔不相等时，是变化的，在不同位置权重不相同，因此需要对其做密度补偿，故将其称为密度补偿函数(density compensation function，DCF)。

点击查看大图

点击查看大图

点击查看大图

3　MRI不相干采样轨迹的评价指标

       就磁共振而言，衡量采样轨迹性能优劣的一个重要指标就是PSF，因为它能准确反映重建图像与采样轨迹的关系。对压缩感知而言，衡量采样矩阵性能的一个重要指标是采样矩阵列相干性。对于基于压缩感知的磁共振成像，该用什么指标来衡量轨迹的性能呢？Lustig给出了一个称之为"TPSF"的评价指标^[3]，其定义为：

       为傅立叶降采矩阵，为矩阵的逆操作，为稀疏矩阵，为矩阵的逆操作，ei、ej分别为第i、j个自然基底向量。当越小时，表示采样矩阵与稀疏矩阵越不相干，也就是说相应的采样轨迹越接近不相干采样。对比压缩感知理论的不相干性定义^[25]，将其表示成Gram矩阵形式如下：

       其不相干性为Gram矩阵的非对角线上的最大值为：

       对比TPSF和压缩感知的不相干定义可以发现，两者都是用来衡量采样矩阵与稀疏矩阵的相干性，其实质是一样的，只是在两个不同领域给出的名称不一样而已。然而，在MRI中，大多数情况下为了衡量采样轨迹的性能都会采用PSF这一指标，那么PSF是否也适用于CS-MRI？是否也能反映采样轨迹的相干采样程度？Lustig同样给出了在笛卡尔坐标系下PSF的定义^[3]。

       并把定义为其相干性。对比式(5)和式(8)，当图像本身是稀疏的，即，E为单位阵，这时TPSF与PSF就是一致的(需要强调的是必须是正交基)，即有：

点击查看大图

点击查看大图

点击查看大图

点击查看大图

点击查看大图

4　采样点位置对PSF形状的影响

       磁共振图像的重建过程实际上就是成像对象与采样轨迹的PSF的卷积，因此PSF的好坏与重建图像的质量密切相关。理想的情况下，PSF是一个冲击函数，这样可以精确重建原始图像。但事实上，在采集信号过程中，由于采样时间受限，通常都需要降采样，而这势必会导致PSF不再是一个理想的冲击函数，而是一个有一定宽度且有旁瓣存在的函数。这种PSF必将引起重建图像出现模糊、伪影等失真的情况。那么在降采样情况下，怎样来评价PSF的好坏呢？为此，有必要考察PSF的各个构成部分，重写PSF的离散形式如下：

       式中为采样点的位置，为对应采样点处的密度补偿函数。由式(10)可知，采样点的位置决定了PSF峰值出现的位置，而密度补偿函数则决定了PSF峰值的高度。理想的情况下PSF是一个冲击函数，也就是说PSF的峰值仅在原点处有最大峰值1，在其他地方都为0。如果PSF除了在原点处存在峰值外，在其他地方仍然有一个比较大的峰值(或者次大值)，那么反映在重建图像上会是什么样的情况呢？假设在笛卡尔坐标系下，重建图像的坐标(10，10)处为感兴趣点，当对应的PSF除了在原点处[对应到重建图像的坐标系就是目标点(10，10)]有一个最大峰值外，在坐标(90，90)处存在另一个峰值，则原始图像的坐标(100，100)处的像素值将会叠加在重建图像坐标(10，10)处。换句话说，原始图像坐标(90，90)处的信息将混迭进入重建后的图像坐标(10，10)处，PSF在坐标(90，90)处的峰值高度决定了混迭进入目标点位(10，10)的信号强度的大小，而出现峰值的位置即坐标(90，90)决定了原始图像的哪些点的信息将被混迭进入目标图像。

       为了考察不同采样模式的PSF，本文仅选取两个采样点作为研究对象。此外，为了简单起见，本节仅讨论在笛卡尔网格下的采样，且假设采样点关于中心原点是对称的，则PSF可以简化为：

       由式(11)可知，对于两个采样点的情况下，PSF为余弦函数的和。假定采样轨迹按以下3种方式行走：(1)；(2)；(3)。3种采样模式及其PSF如图1所示。由图1可知，不同采样位置的点，对应的PSF形状也完全不同。因此如何选取采样点的位置对快速重建高质量的磁共振图像非常重要。

点击查看大图

点击查看大图

点击查看大图

图1  采样点位置与PSF的关系。A~C为3种不同采样轨迹，D~F为其对应PSF
图2  4种笛卡尔采样模式。A~D分别为奈奎斯特降采样、行随机降采样、行列均随机降采样和随机降采样
Fig.1  The relationship between the positions of sampling points and PSFs. A—C are three different sampling tracks, D—F are the corresponding PSFs.
Fig.2  Four cartesian sampling modes. A—D are Nyquist down-sampling, random down-sampling in the row direction, all row and column directions, and randomly down-sampling.

5　仿真结果

       考虑以下4种二维采样模式：(1)行、列等间隔采样模式。对图像的行和列两个方向同时降采2倍，总的降采倍数为4，如图2A所示。(2)行随机采样模式。仅随机降采行，降采比为4，如图2B所示。(3)行、列均随机采样模式。同时随机降采行和列，每个方向降采2倍，总降采比为4，如图2C所示。(4)完全随机采样模式。以像素点为最小单位，随机降采，总降采比为4，如图2D所示。为了测试4种采样模式对PSF形状的影响，分别给出了如图3所示的测试图像。为了满足稀疏条件，该测试图像仅由像素值为1、0.8、0.7、0.6、0.5等5个非零像素点构成。由于重建图像等于原始图像与采样模式PSF的卷积，当图像信号仅为冲击函数时，其采样模式的PSF就是重建图像。为了简单起见，本文中的图像都是由一系列冲击函数构成，故对应采样模式的PSF就是对原图像的重建图像。图4、图5分别给出了测试图像在4种不同采样模式下的PSF。

       从图4、图5可以看出，测试图像在4种采样模式下，其PSF的形状各异。采样模式1下，其能量均匀地泄露到有限的4个位置，由于总能量一致，导致出现4个相同的峰值，从而无法区分主、旁瓣。这种PSF反映在重建图像中就是混迭伪影；采样模式2在列方向已经很难重建原始图像了，尽管其PSF的形状在某些行呈随机分布；采样模式3虽然可以重建原始图像，但其旁瓣幅度明显增加了，加大了重建的难度；采样模式4则非常清晰地重现了原始图像，而且其旁瓣幅度仍然非常小，旁瓣的分布也呈随机分布在整个平面。

点击查看大图

图3  原始图像。A、B分别为三维视图和一维视图
Fig.3  The original image. A, B are the overview of 3-D and 1-D.

点击查看大图

图4  原始图像在对应采样模式下的PSF。A~D分别表示图2A~D对应的4种采样模式的重建结果
Fig.4  The PSFs of the original image in the corresponding sampling modes. A—D are the corresponding reconstruction results in the Fig.2.

点击查看大图

图5  不同采样模式的一维PSF。A~D分别表示图4A~D对应的一维PSF
Fig.5  The 1-D PSF of corresponding sampling pattern. A—D are the corresponding 1D PSF in the Fig.4.

6　结论

       通过观察以上PSF的形状可以看出，当主瓣变窄时，会出现几个幅度相对大的旁瓣，但其他大部分旁瓣的幅度却并不高。也就是说，旁瓣本身也存在着幅度高低差异，呈不均匀分布。不均匀的旁瓣分布会导致重建图像呈现类噪声伪影，但这种伪影是良性的，现代信号处理有成熟的理论技术可以去除这些伪影。需要注意的是，不均匀的旁瓣分布并不是消除能量泄露，而是将原来均匀分布的旁瓣导致的折叠伪影转化为噪声伪影。综上可知，目前以主瓣宽度和旁瓣幅度来评价PSF好坏的指标存在优化的空间。如果在固定主瓣宽度的情况下，通过优化旁瓣的分布位置，使其随机地扩散开来，则其每个旁瓣的幅度均会下降(总能量守恒)，从而保证PSF既能满足图像的空间分辨率要求，又能减少混迭伪影的干扰。由于很难定量给出旁瓣的随机分布程度，且采样点的位置选择本身也受制于一些物理因素，因此如何建立基于这种思路的优化模型仍值得探讨。