CN107219481B - 一种通用的三维欠采样轨迹设计方法 - Google Patents

Abstract

本发明涉及磁共振成像，更具体地涉及一种通用的三维欠采样轨迹设计方法。提供了一种通用的三维欠采样轨迹设计方法，所述方法包括以下步骤：i)生成径向或螺旋轨迹线；ii)网格化轨迹空间；iii)从n＝1开始，依次将轨迹线上的轨迹点映射到笛卡尔网格点上；以及iv)重复步骤iii)直到所有轨迹点映射到笛卡尔网格点上。所述方法可以灵活地根据重建算法、成像序列以及图像对比的要求，直接计算出K空间的填充顺序，并用于序列设计。同时，所述方法满足重建算法的要求，能通过CS和PI联合算法、半傅立叶重建算法等重建图像。本发明提出的方法可以提高成像质量，可用于医学成像等。

Description

一种通用的三维欠采样轨迹设计方法

技术领域

本发明涉及磁共振成像，更具体地涉及一种通用的三维欠采样轨迹设计方法。

背景技术

磁共振成像(MRI)是一种强大的医用成像模式，它没有电离辐射，而且能提供多种图像对比度，可用于人体解剖结构、生理功能、血流和代谢信息的成像。然而，与其他成像模式如计算机断层扫描(CT)相比，MRI的成像速度较慢，检查一次约为几分钟至十几分钟，这极大的限制了图像的时间和空间分辨率。而且由于成像过程中需保持静止状态，成像时间过长增加了病人的不舒适感，并且容易受到呼吸运动的影响，产生运动伪影。

为了形成MR图像，需要对信号进行空间定位。三维磁共振成像中，每个体素的信号的位置坐标由相位编码、频率编码以及层编码确定。频率编码梯度场(G_fE)通过使该方向不同位置的自旋有不同的频率来对信号进行一个方向上的空间定位，一般是在读出数据(采集数据)时接通。相位编码梯度场(G_PE)改变不同位置的回波信号的相位，使质子轻微散相。梯度场打开后，由于不同位置的自旋所处的局部磁场不一样，因而进动频率有差异，一段时间后，形成一定的相位关系。梯度场关闭以后，自旋保持相位关系继续以相同的进动频率进动，相位与相位编码梯度场方向的坐标关系一一对应，因而可以确定一个维度的坐标。层编码(G_SS)为另一方向的相位编码，作用原理与前述相位编码一致。但只有一次编码很难兼顾图像对信噪比和空间分辨率的需求，因而需要进行多个幅度不同的相位编码及层编码。每个TR周期内进行层编码，相位编码和频率编码后的信号包含了位置信息。下一TR只改变层编码梯度场和相位编码梯度场幅度。总体采集时间＝重复时间(TR)×相位编码线×层编码线。由于频率编码(FE)方向欠采并不影响总的采集时间，因而三维采集的加速主要发生在相位编码(PE)方向和层编码(SS)方向上。采集轨迹由梯度场的切换来控制，但梯度场受到硬件和电生理的约束，所以采样轨迹需要设计的尽量平滑。

随着磁共振扫描设备硬件的发展，以及新的采集轨迹和重建方法组合的提出，MRI成像所需时间已经大幅度降低。20世纪90年代末，并行成像的出现打破了磁共振梯度和脉冲序列对MRI扫描速度的限制。并行成像技术利用阵列接收线圈的灵敏度信息减少相位编码条数从而加速了MRI数据的采集。并行成像已成为目前的主流技术，在临床上广泛应用，几乎所有的常规序列都能使用并行成像加速。并行成像主要应用的算法是灵敏度编码(SENSE)算法和广义自动校准部分并行采集(GRAPPA)算法。但是，并行成像加速倍数受接收线圈限制，理论上加速倍数不能超过阵列线圈个数。实际应用时，加速倍数通常为两倍。

压缩感知(CS)是近年来提出的一种新的采样理论，可以突破奈奎斯特采样定理的限制，较好的恢复出原始信号。压缩感知磁共振成像(CSMRI)利用k空间信息冗余特性，经过稀疏变换、不相干欠采样、非线性重建三个部分，实现部分k空间数据重建图像的过程，从而加快成像速度。根据CS理论，要完全恢复欠采图像需要满足下列三个条件：(1)MRI图像是可稀疏的；(2)采样域(k空间)与稀疏变换域是高度不相干的；(3)受稀疏变换约束的非线性重建方法。因为图像具有分段光滑的特性，存在大量的冗余信息，所以满足第一个条件。稀疏变换域一般采用小波变换、有限差分或者离散余弦变换，具体使用哪种变换域需要根据信号本身特性来选择。而稀疏程度决定着欠采样数据量的数目。稀疏变换域是确定的，为了满足不相干特性，因而需要在k空间进行随机采样，并且采样越随机，CS重建的效果越好。

由于并行成像和压缩感知利用信息不同，可以将二者结合起来获得更高的加速倍数，这就要求采样模板的采样尽量随机并且分布较为均匀。现有方法，网格化泊松盘(Poisson-disk)分布的三维笛卡尔(Cartesian)采样模板很好地平衡了均匀性和不相干的需求。泊松盘采样是一种随机的采样过程，通过限定相邻采样间的最小距离，以免在k空间出现采集点局部聚集的情况，使得分布更加均匀。然而，信号的能量在k空间中是不均匀的，图像大量信息集中在k空间的低频区域(也即k空间中心)。因而在k空间需要对采样点间的距离进行约束，并且进行变密度数据采集(中心的采集需要更密集，往k空间高频区域或k空间边缘方向越来越稀疏)或者中心区域全采。

网格化泊松盘分布的三维笛卡尔模板通常在相位编码和层编码方向采用变密度Poisson-disk分布。基本思想是在现有采集点(每个点代表其在频率编码(kx)方向为一条信号读出线，下同)的周围生成数个点，并检查它们是否满足与其它现有采集点的限定距离，满足则添加为下一个采集点。在新生成的点被接纳为采样点之前，我们必须检查该点与以前生成的点是否离得太近。依照一定的选择标准将所有潜在采集点都进行判断，最终得到采样模板。通过对最小距离的设计(如将距中心距离的平方根作为最小距离)，可以得到变密度的Poisson-disk采样。

但是，这种随机采样仍可能会出现采集点与采集点之间的距离过大的问题，这将使CS和并行成像(PI)联合重建算法不稳定。并且该采样技术只生成了采样模板，但无法直接和序列中K空间的采集顺序相结合。

发明内容

本发明的解决方案的目的是解决前面强调的问题。

从前述可知数据采集与梯度场有关，由梯度场的切换控制采样轨迹。现有技术(例如，泊松盘(Poisson-disk)采样)仅生成了采样模板，但无法直接和序列中K空间的采集顺序相结合。结合笛卡尔采样及非笛卡尔采样的各自的优势，我们提出了一种通用的三维笛卡尔采样轨迹设计方法。所述方法可以灵活地根据重建算法、成像序列以及图像对比的要求，直接计算出K空间的填充顺序，并用于序列设计。同时，所述方法满足重建算法的要求，能通过CS和PI联合算法、半傅立叶重建算法等重建图像。本发明提出的方法可以提高成像质量，可用于医学成像等。

因此，根据本发明的解决方案，提供了一种通用的三维欠采样轨迹设计方法。本发明采用允许变密度采样的类径向或类螺旋路径的旋转笛卡尔采样轨迹。

本发明提供了一种通用的三维欠采样轨迹设计方法，其可以包括以下步骤：i)生成径向或螺旋轨迹线；ii)网格化轨迹空间；iii)从n＝1开始，依次将轨迹线上的轨迹点映射到笛卡尔网格点上；以及iv)重复步骤iii)直到所有轨迹点映射到笛卡尔网格点上。

在一些实施方式中，其中步骤i)中所述的径向或螺旋轨迹线可以表示为：

其中y_n(r_m)和z_n(r_m)为第n条轨迹线的空间坐标，r_m为轨迹点半径，是黄金角，κ是相对角速度，N是轨迹线数，M为轨迹线上的轨迹点个数，k空间总采集点数为N×M。并且其中轨迹点半径的绝对值|r_m|为轨迹点距轨迹填充中心(y＝0,z＝0)的距离，当r_m为负时表示轨迹点在填充中心的左侧，当r_m为正时表示轨迹点在填充中心的右侧。

在一些实施方式中，

在一些实施方式中，

在一些实施方式中，r_m∈[0,1]。

在一些实施方式中，r_m∈[-a,1]或者r_m∈[a,-1]，其中a＝0.5或0.75。

在一些实施方式中，r_m∈[-1,1]。

在一些实施方式中，r_m+1＝r_m+Δr,m＝[0…M-1]，其中r_m，r_m+1为同一轨迹线上相邻轨迹点半径，Δr为增量；并且，当Δr>0时，从中心向外围采样；而当Δr<0时，从外围向中心采样。

在一些实施方式中，当Δr＝常数时，采样为均匀采样；而当Δr＝变量时，采样为变密度采样。

在一些实施方式中，可以通过参数κ控制轨迹线类型，当κ＝0时轨迹线类型为径向轨迹线；κ>0时轨迹线类型为螺旋轨迹线，并且κ越大，轨迹旋转速度越快。

在一些实施方式中，在所述步骤ii)中，可以设k_v方向采样线数为N_y，k_z方向采样线数为N_z，加速倍数为

在一些实施方式中，在所述步骤iii)中，所述第n条轨迹上的轨迹点坐标为(y_n(r_m),z_n(r_m))，可以从最临近当前路径上点的区域开始搜索潜在的笛卡尔网格点(k_y，k_z)：

a.若最近邻的笛卡尔网格点(k_y，k_z)未被采集过，则直接将y_n(r_m),z_n(r_m))映射到(k_y，k_z)，然后进行下一个轨迹点的映射搜索；

b.若最近邻的笛卡尔网格点(k_y，k_z)已经采集，则在以(k_y，k_z)为中心，设定D(y,z)为搜索范围矩阵，在大小为Sr×Sr的范围内进行搜索，选择所述范围内以下代价函数最小的网格点：

其中是潜在网格点(k_y，k_z)的方位角，θ_y,z是当前轨迹点的方位角，λ是权重因子(0<λ<1)，D(y,z)是在轨迹点(y,z)上的搜索范围矩阵，若在设定大小内未搜索到合适网格点，则需扩大Sr×Sr继续搜索直到轨迹点被成功映射到网格点上。

在一些实施方式中，可以通过约束(k_y，k_z)的搜索范围，避免网格点的重复采样，也防止网格点落在太远的区域。

本领域技术人员在阅读整个说明书和权利要求书时将理解本发明的这些优点和其它优点。

附图说明

为了更好地理解本发明，现在仅通过非限制性示例的方式参照所附附图描述其优选实施方式，其中：

图1示出了根据本发明的方法的步骤的流程图。

图2示出了根据本发明的方法的黄金角类螺旋三维笛卡尔采样随时间的采样过程。

图3示出了根据本发明的方法的黄金角类螺旋轨迹的三维笛卡尔采样。

图4示出了根据本发明的方法的黄金角类径向轨迹的三维笛卡尔采样。

图5a、5b分别示出了根据泊松盘采样方法获得的采样模板与根据本发明的方法的采样模板。

图6示出了使用三种模板采样(从左至右分别为：全采集方法、本发明提出的方法、泊松盘采样方法)后重建得到的图像(上)及与参考图像的误差(下)。

具体实施方式

下面结合附图对本发明的具体实施例进行说明。在下文所描述的本发明的具体实施例中，为了能更好地理解本发明而描述了一些很具体的技术特征，但显而易见的是，对于本领域的技术人员来说，并不是所有的这些技术特征都是实现本发明的必要技术特征。下文所描述的本发明的一些具体实施例只是本发明的一些示例性的具体实施例，其不应被视为对本发明的限制。另外，为了避免使本发明变得难以理解，对于一些公知的技术没有进行描述。

非笛卡尔如径向、螺旋轨迹等欠采样本身与稀疏域高度不相干，并且比笛卡尔采样对运动更不敏感。但受到磁共振硬件等的限制，非笛卡尔轨迹成像质量不稳定。

另外，相邻采集轨迹通过黄金角旋转可以避免采集相同数据点，并覆盖k空间的大部分信息。黄金角最早发现于植物学科，科学家们研究发现新叶子会尽量避免与上一片叶子重叠，以获得最大的光照，并发现叶子偏转的角度为360°除以黄金分割比率g＝(1+√5)/2＝1.618。发现这个最优问题与CT、PET等的投影角度的最优排序是类似的。随着径向和螺旋轨迹等非笛卡尔采样的出现，Winkelmann将此黄金角引入磁共振成像的非笛卡尔采样模板中。黄金角采样确保了信号在k空间中分布较为均匀，在动态成像中可以灵活选择时间分辨率。

结合黄金角和非笛卡尔轨迹的优势，我们在三维笛卡尔采样上通过控制类似径向或者类似螺旋轨迹的相对角速度和半径，以及黄金角为方位角的旋转方式设计采样轨迹，并通过距离和邻域点约束将其映射到三维笛卡尔上的采样点。

在图1中，示出了根据本发明的方法的步骤的流程图。变密度三维笛卡尔采样模板的设计过程如下：

i)首先，生成径向或螺旋轨迹线：

其中y_n(r_m)和z_n(r_m)为第n条轨迹线的空间坐标，r_m为轨迹点半径，是黄金角，κ是相对角速度，N是轨迹线数，M为轨迹线上的轨迹点个数。其中轨迹点半径的绝对值|r_m|为轨迹点距轨迹填充中心(y＝0,z＝0)的距离，当r_m为负时表示轨迹点在填充中心的左侧，当r_m为正时表示轨迹点在填充中心的右侧。参数r_m根据用户需求定义，一般范围设定在-1到1之间。调节r_m可以控制轨迹线的方向、密度的变化。令r_m+1＝r_m+Δr,m＝[0…M-1],其中r_m，r_m+1为同一轨迹线上相邻轨迹点半径，Δr为增量。当Δr＝C(常数)时，采样为均匀采样；当Δr＝var(变量)时，采样为变密度采样。通过参数κ控制轨迹类型为径向或者螺旋轨迹，当κ＝0时为径向轨迹；κ>0时为螺旋轨迹,κ越大，轨迹旋转速度越快。k空间总采集点数为N×M。

ii)网格化轨迹空间：

设k_y方向采样线数为N_y，k_z方向采样线数为N_z，加速倍数为

iii)从n＝1开始，依次将轨迹线上的轨迹点映射到笛卡尔网格点上。轨迹上的轨迹点坐标为(y_n(r_m),z_n(r_m))，从最临近当前路径上点的区域开始搜索潜在的笛卡尔网格点(k_y，k_z)：

1)若最近邻的笛卡尔网格点(k_y，k_z)未被采集过，则直接将(y_n(r_m),z_n(r_m))映射到(k_y，k_z)，进行下一个轨迹点的映射搜索。

2)若最近邻的笛卡尔网格点(k_y，k_z)已经采集，则在以(k_y，k_z)为中心的，设定D(y,z)为搜索范围矩阵，大小为Sr×Sr的范围(最初可以设置成5×5)内进行搜索，选择该范围内以下代价函数[2]最小的网格点：

其中，是潜在网格点(k_y，k_z)的方位角，θ_y,z是当前轨迹点的方位角，λ是权重因子(0<λ<1)，D(y,z)是在轨迹点(y,z)上的搜索范围矩阵，若在设定大小内未搜索到合适网格点，则需扩大Sr×Sr继续搜索直到轨迹点被成功映射到网格点上。通过约束(k_y，k_z)的搜索范围，避免网格点的重复采样，也防止网格点落在太远的区域。

iv)重复步骤iii)直到所有轨迹点映射到笛卡尔网格点上。

下面介绍几种具体的参数设置，作为实施例：

a)类螺旋轨迹采样

基本参数：κ>0，N×M根据用户需求设置

轨迹点半径：r_m∈[0,1]

采样模式：Δr＝C(均匀)或者Δr＝var(变密度)

采集顺序：r₀＝0，r_M-1＝1，Δr>0(从中心向外填充),或者

r₀＝1，r_M-1＝0，Δr<0(从外围向中心填充)。

b)类径向轨迹采样

(1)由中心向外四周发散

基本参数：κ＝0且,N×M根据用户需求设置

轨迹点半径：r_m∈[0,1]

采集顺序选择同上，

采样模式选择同上。

(2)横贯中心

基本参数：κ＝0,N×M根据用户需求设置

轨迹点半径：r_m∈[-1,1]

采集顺序：r₀＝0，r_M-1＝1，Δr>0(从k空间一端的外围穿过中心到另一端的外围),或者

r₀＝1，r_M-1＝0，Δr<0(从k空间另一端的外围穿过中心到一端的外围)，

采样模式选择同上。

c)半傅立叶采样模板

基本参数：κ＝0,N×M根据用户需求设置

轨迹点半径：r_m∈[-a,1]或者r_m∈[a,-1]，其中a根据采集数据的多少设定，如1/2傅立叶采集则a＝0.5，3/4傅立叶采集则a＝0.75等。

采集顺序选择同上，

采样模式选择同上。

根据图像对比度需求设定采集顺序，若想获得T1加权图像则采样从外围向中心填充，若想获得T2加权图像则采样从中心向外填充。为迎合k空间信号的能量是不均匀的，即图像大量信息集中在k空间的低频区域(也即k空间中心)，Δr应设置成随距离k空间中心位置越远，Δr越大。

现在参考图2-4，图2示出了根据本发明的方法的黄金角类螺旋三维笛卡尔采样随时间的采样过程，其为黄金角类螺旋轨迹线的形成示意图，将轨迹线上的轨迹点映射为笛卡尔网格点后即可得到图5b，将上述a)中的相关参数代入式[1]即可生成图2。图3示出了根据本发明的方法的黄金角类螺旋轨迹的三维笛卡尔采样，其中为了强调黄金角只显示了两条螺旋轨迹线。图4示出了根据本发明的方法的黄金角类径向轨迹的三维笛卡尔采样，其用黄金角类径向轨迹来说明采样点从径向轨迹线映射到笛卡尔网格上的过程。如图4所示，点(灰点以及已经映射的黑点)代表的是笛卡尔网格，虚线代表径向采样，黑点表示的是点已经从径向轨迹线上映射到了笛卡尔网格上，而映射是搜索周围区域未被填充的笛卡尔网格点利用代价函数 作为约束条件来找到对应的映射位置的。

本发明提供的通用采样模板并不影响成像时间，并且直接将采样模板和采样轨迹结合在一起，无需另外设计采样轨迹。参考本发明中提供的几种参数设置可以直接根据所需的图像对比度设置参数。而泊松盘(Poisson-disk)采样轨迹的核心是泊松盘采样模板，实际运用时需要另外设计采样轨迹，而采样轨迹的设计方法因人而异，没有一个统一的标准。所以使用泊松盘采样轨迹生成所需图像对比度需另外设计采样轨迹与该采样模板一同使用。

现在对使用本发明提供的采样模板与使用泊松盘采样模板得到k空间数据采用同一重建方法得到的图像可以进行图像质量的比较。图像质量比较的实验步骤：1)全采一幅图像，取得该图像的k空间数据，直接重建作为参考图像；2)分别用泊松盘分布模板及我们提出的通用模板对k空间数据仿真欠采样过程。3)对上述模板得到的新k空间数据用同一种重建方法重建图像。

现在参考图5a-5b，图5a-5b分别示出了根据泊松盘采样方法获得的采样模板与根据本发明的方法的采样模板，其中本发明提出的采样模板在相同参数条件下加速倍数(加速倍数等于实际采样点数比上全采的点数)要大于泊松盘分布采样模板。

图5a-5b所示生成图像大小的采样模板结果如图6所示。图6示出了使用三种模板采样(从左至右分别为：全采集方法、本发明提出的方法、泊松盘采样方法)后重建得到的图像(上)及与参考图像的误差(下)。从图6中可以看出，本发明提出的采样模板重建得到的图像与参考图像几乎一致，而泊松盘分布采样重建后得到的图像细节丢失严重，如大部分灰质信号丢失(见白框内区域)。从整个成像及其误差图像上看，我们提出的采样模板要显然优于泊松盘分布采样模板。

本发明提出的设计方法可以灵活地根据重建算法、成像序列以及图像对比的要求，直接计算出K空间的填充顺序，并用于序列设计。同时，所述方法满足重建算法的要求，能通过CS和PI联合算法、半傅立叶重建算法等重建图像。本发明提出的方法可以提高成像质量，可用于医学成像等。

尽管已经根据优选的实施方案对本发明进行了说明，但是存在落入本发明范围之内的改动、置换以及各种替代等同方案。还应当注意的是，存在多种实现本发明的工艺的可选方式。因此，意在将所附的权利要求书解释为包含落在本发明的主旨和范围之内的所有这些改动、置换以及各种替代等同方案。

Claims (12)

1.一种通用的三维欠采样轨迹设计方法，其特征在于包括以下步骤：

i)生成径向或螺旋轨迹线；

ii)网格化轨迹空间；

iii)从n＝1开始，依次将所述轨迹线上的轨迹点映射到笛卡尔网格点上；以及

iv)重复步骤iii)直到所有轨迹点映射到笛卡尔网格点上；

其中，在所述步骤iii)中，第n条轨迹线上的轨迹点坐标为(y_n(r_m)，z_n(r_m))，从最临近当前路径上点的区域开始搜索潜在的笛卡尔网格点(k_y，k_z)：

a.若最近邻的笛卡尔网格点(k_y，k_z)未被采集过，则直接将(y_n(r_m)，z_n(r_m))映射到(k_y，k_z)，然后进行下一个轨迹点的映射搜索；

b.若最近邻的笛卡尔网格点(k_y，k_z)已经采集，则在以(k_y，k_z)为中心，设定D(y，z)为搜索范围矩阵，在大小为Sr×Sr的范围内进行搜索，选择所述范围内以下代价函数最小的网格点：

其中是潜在网格点(k_y，k_z)的方位角，θ_y，z是当前轨迹点的方位角，λ是权重因子，其中0＜λ＜1，D(y，z)是在轨迹点(y，z)上的搜索范围矩阵，若在设定大小内未搜索到合适的网格点，则需扩大Sr×Sr继续搜索直到所述轨迹点被成功映射到网格点上。

2.根据权利要求1所述的方法，其特征在于，其中步骤i)中所述的径向或螺旋轨迹线为：

其中y_n(r_m)和z_n(r_m)为第n条轨迹线的空间坐标，r_m为轨迹点半径，是黄金角，κ是相对角速度，N是轨迹线数，M为轨迹线上的轨迹点个数，k空间总采集点数为N×M；并且

其中所述轨迹点半径的绝对值|r_m|为轨迹点距轨迹填充中心(y＝0，z＝0)的距离，当r_m为负时表示轨迹点在填充中心的左侧，当r_m为正时表示轨迹点在填充中心的右侧。

3.根据权利要求2所述的方法，其特征在于：

4.根据权利要求2所述的方法，其特征在于：

5.根据权利要求2所述的方法，其特征在于：r_m∈[0，1]。

6.根据权利要求2所述的方法，其特征在于：r_m∈[-a，1]或者r_m∈[a，-1]，其中a＝0.5或0.75。

7.根据权利要求2所述的方法，其特征在于：r_m∈[-1，1]。

8.根据权利要求2所述的方法，其特征在于：

r_m+1＝r_m+Δr，m＝[0...M-1]，其中r_m，r_m+1为同一轨迹线上相邻轨迹点半径，Δr为增量；并且，

当Δr＞0时，从中心向外围采样；而

当Δr＜0时，从外围向中心采样。

9.根据权利要求8所述的方法，其特征在于：

当Δr＝常数时，采样为均匀采样；而

当Δr＝变量时，采样为变密度采样。

10.根据权利要求2所述的方法，其特征在于：

通过参数κ控制轨迹线类型，当κ＝0时轨迹线类型为径向轨迹线；κ＞0时轨迹线类型为螺旋轨迹线，并且κ越大，轨迹旋转速度越快。

11.根据权利要求1所述的方法，其特征在于：

在所述步骤ii)中，设k_y方向采样线数为N_y，k_z方向采样线数为N_z，加速倍数为

12.根据权利要求1所述的方法，其特征在于：

通过约束(k_y，k_z)的搜索范围，避免网格点的重复采样，也防止网格点落在太远的区域。