并行成像技术与MRI加速方法

MRI加速技术概述

MRI成像的一个核心挑战是采集时间过长。传统的全采样方式需要逐行采集k空间数据，耗时较长，这在临床应用中会导致运动伪影等问题。因此，一个改进思路是，能否通过采集更少的k空间样本来缩短采集时间，同时仍然能够重建出清晰的图像？这就是欠采样（Undersampling）技术要解决的核心问题。

欠采样技术的分类

根据硬件配置的不同，欠采样技术可以分为两大类。第一类基于多接收线圈系统，包括SMS（Simultaneous Multi-Slice，同时多层激发）和并行成像（Parallel Imaging）技术。第二类基于单接收线圈系统，包括部分傅里叶采集（Partial Fourier）和压缩感知（Compressed Sensing）技术。

这些加速技术存在共同的限制：加速因子 R 通常不超过6，即 R \leq 6。更关键的是，信噪比（SNR）会随着加速因子 R 的增大而快速下降，因此所有欠采样方法都需要思考权衡问题。

并行成像与压缩感知的对比

并行成像和压缩感知代表了两种不同的加速重建框架。并行成像采用笛卡尔欠采样模式，依赖多接收线圈提供的空间编码冗余信息，通过线性重建算法来恢复完整图像。压缩感知则采用非笛卡尔的变密度采样模式（如SPARKLING轨迹），利用图像在某个变换域（字典）中的稀疏性先验，通过非线性重建算法来求解欠定问题。

并行成像的基本原理

多线圈信号融合

并行成像的核心思想是利用多个接收线圈同时采集信号。每个线圈由于其物理位置不同，具有不同的空间敏感度分布，这种敏感度差异为图像重建提供了额外的空间编码信息。当使用多个线圈元件（如4个线圈）进行采集时，每个线圈会产生一幅带有其敏感度加权的图像。将这些原始数据和临时图像进行适当组合，就能重建出最终的完整图像。

并行成像的优缺点

并行成像技术的优势：减少扫描时间，提高空间或时间分辨率 (在相同的扫描时间内可以采集更多数据)，限制几何畸变（更短的回波链减少了场不均匀性的累积效应）。

并行成像技术的缺点：信噪比降低（欠采样本身就减少了信号平均次数），线圈敏感度的非均匀性（因为会引入额外的图像伪影和噪声放大效应）

并行成像的两种主要方法

并行成像技术主要分为两种实现路径，它们的区别在于处理欠采样数据的域不同。

GRAPPA方法

GRAPPA（GeneRalized Autocalibrating Partially Parallel Acquisitions）方法在k空间域进行操作，由Griswold等人于2002年提出。当采用加速因子 R=2 进行欠采样时，k空间中相邻采集行的间距变为 2\Delta k_y，即只采集了一半的相位编码线。GRAPPA的策略是利用多线圈数据之间的相关性，通过线性组合来估计并填充那些未采集的k空间行，使行间距恢复为 \Delta k_y。填充完成后的完整k空间数据经过逆傅里叶变换（IFFT）即可得到完整视场（FOV）的无混叠图像。

SENSE方法

SENSE（SENSitivity Encoding）方法在图像域进行操作，由Pruessmann等人于1999年提出。与GRAPPA不同，SENSE首先对欠采样的k空间数据直接进行逆傅里叶变换。由于k空间欠采样导致采样间距增大，根据傅里叶变换的性质，图像域的视场会相应缩小为原来的 \text{FOV}/2（对于 R=2 的情况），产生混叠伪影，即图像的不同部分会折叠重叠在一起。SENSE方法随后利用各线圈的敏感度图信息，在图像域通过求解线性方程组来展开这些混叠，恢复出完整FOV的图像。

这两种方法本质上是等价的，都利用了线圈敏感度提供的空间编码信息来解决欠采样带来的信息缺失问题，只是操作的域不同：GRAPPA在k空间补全数据，SENSE在图像域分离混叠。

SENSE方法的数学原理

SENSE（Sensitivity Encoding，敏感度编码）是Pruessmann等人于1999年在Magnetic Resonance in Medicine期刊上发表的经典并行成像方法。该方法的核心思想是在图像域处理欠采样导致的混叠问题，利用多线圈敏感度信息将重叠的像素分离开来。

混叠现象的数学描述

当对k空间进行欠采样后，首先需要将数据变换回图像域。设 y 为欠采样的k空间数据，通过对采样位置 \Omega 上的数据进行逆傅里叶变换，得到混叠图像：

这里 \mathcal{F}_\Omega^* 表示在采样集合 \Omega 上的伴随傅里叶算子（即逆傅里叶变换）。由于k空间欠采样，图像域的视场缩小为 \text{FOV}/2（以加速因子 R=2 为例），导致原本位于不同空间位置的像素值在混叠图像中叠加到了同一个位置。

对于混叠图像中某一位置 r_1 处的像素值 d(r_1)，它实际上是原始完整图像中两个位置像素值的加权叠加：

其中:

• x(r_1) 是原始图像在位置 r_1 处的真实像素值
• x(r_1 + \frac{FOV}{2}) 是与 r_1 相距半个视场的另一位置的像素值
• 系数 \alpha 和 (1-\alpha) 取决于具体的采样模式

这个方程表明仅凭单个线圈的数据，无法从一个已知量 d(r_1) 中解出两个未知量 x(r_1) 和 x(r_1 + \frac{FOV}{2})，方程是欠定的。

多线圈展开的建模

两线圈情形

SENSE方法通过引入多个接收线圈来解决上述欠定问题。考虑使用2个线圈、加速因子 R=2 的情况。定义两个混叠位置的坐标：

每个线圈都有其特定的空间敏感度分布，用 C_l(r) 表示第 l 个线圈在位置 r 处的敏感度值。

线圈敏感度描述了该线圈对不同空间位置信号的接收效率，它取决于线圈的几何形状和物理位置。对于第1个线圈，其混叠图像在位置 r_1 处的像素值为：

线圈1在位置 r_1 处观测到的混叠值，等于真实图像在 r_1 位置的值乘以线圈1在该位置的敏感度，加上真实图像在 r_2 位置的值乘以线圈1在 r_2 位置的敏感度。

类似地，对于第2个线圈：

将这两个方程写成矩阵形式：

现在有2个方程、2个未知数，只要敏感度矩阵可逆（即两个线圈的敏感度分布线性无关），就可以唯一求解出 x(r_1) 和 x(r_2)。

一般情形的推广

将上述分析推广到 L 个线圈、加速因子为 R 的一般情况。此时，混叠图像中每个位置 r_1 实际上是 R 个不同位置像素值的叠加。这 R 个位置分别为：

对于第 L 个线圈，其混叠方程为：

将所有 L 个线圈的方程整合为矩阵形式：

这个方程可以用紧凑的向量形式表示。定义 \underline{r} = [r_1, \ldots, r_R] 为所有混叠位置的集合，则对于视场内的每一个位置 r_1，都有：

其中 d(r_1) 是 L \times 1 的观测向量，C(\underline{r}) 是 L \times R 的敏感度矩阵，x(\underline{r}) 是 R \times 1 的待求解像素值向量。

并行求解结构

SENSE重建的结构使其可以并行计算：对于图像中的每一个折叠位置 r_1，都需要独立求解一个小型线性系统。这些线性系统之间相互独立，可以并行计算，因此SENSE方法具有良好的计算效率。当线圈数 L 大于等于加速因子 R 时，方程组是超定或恰定的，可以通过最小二乘或直接求逆得到解；当 L < R 时，方程组欠定，需要额外的正则化约束。

SENSE重建的求解公式

前面建立了SENSE的线性方程组模型

现在需要求解这个方程组以获得原始图像的像素值。由于存在噪声，且线圈数 L 通常大于加速因子 R，方程组是超定的，需要采用加权最小二乘方法求解。

对于视场内的每一个位置 r_1，SENSE重建的解析解为：

这个公式是加权最小二乘估计的标准形式。其中:

• C^H 表示敏感度矩阵的共轭转置
• \Sigma 是噪声协方差矩阵，描述了不同线圈之间噪声的相关性

引入 \Sigma^{-1} 作为权重矩阵的目的是在求解时考虑噪声特性，给噪声较小的线圈数据更大的权重，从而获得最优的信噪比。

加速效果的直观展示

SENSE方法可以在保持空间分辨率不变的前提下减少扫描时间。以 1 \times 1 \times 1.1 \text{ mm}^3 分辨率的脑部成像为例：标准全采样需要9分14秒；使用 R=2 的并行采集，时间缩短为5分03秒；使用 R=4 的并行采集，时间进一步缩短为2分59秒。从重建图像可以看出，随着加速因子增大，虽然扫描时间大幅减少，但图像质量会有所下降，主要表现为噪声增加。

线圈敏感度图的估计

SENSE方法的实现依赖于准确的线圈敏感度图 C_\ell。敏感度图需要针对每个扫描对象单独测量，因为不同对象的几何形状和位置会影响线圈的有效敏感度分布。敏感度图可以在正式扫描之前单独采集，但这会引入运动或配准误差的问题——如果患者在敏感度校准扫描和正式扫描之间发生了移动，敏感度图就不再准确。由于线圈敏感度在空间上是平滑变化的，因此只需要低分辨率的图像就足够估计敏感度分布。

敏感度图估计流程

敏感度图的估计分为三个步骤。第一步是提取k空间中心的低频数据。对于每个线圈 \ell，将其完整k空间数据 y_\ell 与一个低通滤波器 M_{LF} 相乘，只保留中心区域的低频成分：

k空间中心包含了图像的低频信息，即整体的亮度和对比度分布，这正好对应于敏感度的平滑变化特性。

第二步是对每个线圈的低频数据进行逆傅里叶变换，得到低分辨率的线圈图像：

第三步是归一化处理。将每个线圈的低分辨率图像除以所有线圈图像的平方和开根号（即Sum-of-Squares组合图像的幅值），得到归一化的敏感度图：

这种归一化确保了敏感度图反映的是各线圈之间的相对敏感度差异，而不是绝对信号强度。归一化后，所有线圈在同一位置的敏感度平方和为1。

加速因子与信噪比的权衡

SNR损失公式

并行成像的加速是有代价的，主要体现为信噪比的下降。设 SNR_0 为不使用并行成像时的全采样信噪比，SNR_R 为使用加速因子 R 时的信噪比，两者之间的关系为：

这个公式包含两个导致SNR下降的因素。第一个是 \sqrt{R} 因子，这是由于欠采样本身造成的：采集的数据点减少为原来的 1/R，根据信号平均的原理，信噪比会下降 \sqrt{R} 倍，所有欠采样方法都无法避免此限制。第二个是几何因子 g（g-factor），这是SENSE方法特有的额外SNR损失，源于线圈敏感度矩阵的条件数。

几何因子的特性

几何因子 g 反映了线圈几何配置对重建质量的影响，与敏感度矩阵求逆过程中的误差放大有关。当敏感度矩阵接近奇异时，g 值会很大，导致严重的噪声放大。

g 因子具有以下特性：

• 当 R=1（即不加速）时，g=1，没有额外的SNR损失

• 当 R=2 时，g 通常在1.5到2之间

• g 因子会随着线圈通道数的增加而减小，因为更多的线圈提供了更丰富的空间编码信息，使敏感度矩阵的条件数改善

• g 因子也会随着主磁场强度 B_0 的增加而减小，因为更高场强下线圈敏感度的空间变化更加显著。

g因子的空间分布

g 因子不是一个常数，而是随空间位置变化的。这是因为对于图像中不同的折叠位置组合，敏感度矩阵 C(\underline{r}) 是不同的，因此矩阵的条件数也不同。

从g-map图像可以观察到，随着加速因子从1增加到4，g 因子的空间分布变得越来越不均匀，某些区域（通常是图像中心或线圈敏感度变化不大的区域）会出现很高的 g 值，导致这些区域的噪声被严重放大。对应的SENSE重建图像也显示出随着 R 增大，图像质量逐渐下降，噪声增加，特别是在 g 因子较高的区域。

SENSE方法特点总结

SENSE的处理流程是先对欠采样k空间数据进行逆傅里叶变换得到混叠图像，然后在图像域利用敏感度编码信息进行展开。这种处理顺序使得逆问题大大简化：整个重建问题在空间上是可分离的，即可以分解为针对每组 R 个混叠体素的独立小问题，其中 R 既是加速因子也是混叠因子。这意味着不需要求解一个涉及整幅图像所有像素的巨型方程组，而是并行求解大量相互独立的小型方程组，计算效率很高。

由于不同体素组对应的敏感度矩阵不同，矩阵的条件数也不同，因此噪声放大效应是空间变化的。这种空间变化的噪声放大用几何因子（g-factor）来量化。SENSE方法的一个前提条件是必须显式地知道或测量线圈敏感度图，这增加了方法的实现复杂度，也是其与GRAPPA方法的主要区别之一。

GRAPPA方法的基本原理

GRAPPA（Generalized Autocalibrating Partially Parallel Acquisitions，广义自校准部分并行采集）是Griswold等人于2002年在Magnetic Resonance in Medicine期刊上发表的并行成像方法。与SENSE在图像域处理混叠不同，GRAPPA直接在k空间域工作，通过估计和填充未采集的k空间数据点来实现图像重建。

k空间数据填充的思想

GRAPPA的核心策略是利用多线圈数据之间的冗余性，从已采集的k空间点（源点）来估计未采集的k空间点（目标点）。以加速因子 R=3 为例，每采集一行k空间数据就跳过两行，导致k空间中存在大量空白行。

GRAPPA的目标就是将这些空白位置的数据估计出来，使k空间恢复完整，然后通过标准的逆傅里叶变换即可得到无混叠的图像。

卷积核的几何结构

GRAPPA使用一个局部的卷积核来建立源点和目标点之间的线性关系。卷积核的几何形状定义了哪些源点参与估计某个目标点。例如，选择 [3,2] 的核几何意味着在 (k_x, k_y) 方向上分别取3个和2个采样点作为源点。这个核的形状决定了源点和目标点之间的空间关系：源点来自目标点周围的邻域，跨越多个已采集的k空间行。

卷积核权重建立了一个从所有线圈的源点到单个线圈中某个目标点的线性映射。在图像域中，每个线圈的图像等于真实图像乘以该线圈的敏感度图，而乘法在图像域对应于k空间中的卷积操作。因此，不同线圈的k空间数据之间存在卷积关系，这正是GRAPPA能够利用线圈间相关性的数学基础。

GRAPPA权重的矩阵表示

单个目标点的估计

以2个线圈、R=3 的情况为例来说明权重的具体含义。假设卷积核包含12个源点（来自两个线圈，每个线圈6个点），需要估计的目标点在两个线圈中分别记为 X（线圈1）和 Y（线圈2）。源点用编号1到12标识，其中1-6来自线圈1的已采集行，7-12来自线圈2的已采集行。

对于线圈1中的目标点 X，其估计值是所有源点的加权线性组合：

类似地，对于线圈2中的目标点 Y：

这里 w_{l,j} 表示第 l 个目标线圈使用第 j 个源点的权重系数。每个目标线圈都有自己的一组权重，用于从相同的源点集合中估计该线圈的缺失数据。

矩阵形式的表达

将所有权重组织成矩阵形式可以更清晰地表达GRAPPA的计算结构。定义源点向量 s_{\text{source}} 包含所有源点的k空间值，目标点向量 s_{\text{target}} 包含所有目标线圈在该位置的待估计值。对于第一个卷积核位置，权重矩阵 W^{(1)} 的结构为：

用紧凑的符号表示为：

在这个矩阵中，每一行对应一个目标线圈，每一列对应一个源点。矩阵的维度是（目标线圈数）\times（源点总数），在本例中为 2 \times 12。权重矩阵 W^{(1)} 的上标 (1) 表示这是针对第一个缺失位置（相对于已采集行的第一个空白行）的权重，对于 R=3 的情况，还需要另一个权重矩阵 W^{(2)} 来估计第二个缺失位置的数据。

多卷积核与平移不变性

第二个卷积核设置

对于加速因子 R=3 的情况，每采集一行就跳过两行，因此在两条已采集行之间存在两个不同位置的缺失点。第一个缺失位置（紧邻上方已采集行）使用权重矩阵 W^{(1)}，第二个缺失位置（紧邻下方已采集行）需要使用另一套权重矩阵 W^{(2)}。对于第二个卷积核位置，目标点记为 U（线圈1）和 V（线圈2），其估计公式为：

同样用紧凑形式表示为：

核类型数量与权重矩阵维度

对于加速因子 R，k空间中存在 R-1 种不同类型的缺失位置，因此需要 R-1 套不同的权重矩阵。每个缺失点都属于这 R-1 类中的某一类，通过判断其相对于最近已采集行的位置来确定应该使用哪套权重。

将所有核类型整合到一个统一的矩阵方程中，可以同时估计同一源点邻域内所有缺失位置的目标点。对于 L 个线圈、核大小 [k_x, k_y]、加速因子 R 的系统，完整的权重矩阵维度为 L \times (R-1) 行（对应所有目标线圈和所有核类型）乘以 L \times k_x \times k_y 列（对应所有源点）。整合后的矩阵方程为：

这种矩阵结构的优势在于可以一次性处理所有核类型，提高计算效率。

GRAPPA权重的校准

自校准信号区域

GRAPPA方法的一个核心优势是自校准（Autocalibrating）特性，即权重矩阵可以从数据本身学习得到，而不需要单独的校准扫描。这是通过在k空间中心额外采集一些全采样的校准数据来实现的，这个区域称为ACS（Auto-Calibration Signal）区域。

由于k空间中心包含图像的低频信息且信噪比较高，在此区域进行全采样可以提供可靠的训练数据来估计GRAPPA权重。

平移不变性假设

GRAPPA权重校准的数学基础是平移不变性假设：同一套权重可以应用于k空间中任意位置，只要源点和目标点之间的相对几何关系保持不变。这个假设来源于卷积的平移不变特性——线圈敏感度在图像域的乘法对应k空间中的卷积，而卷积核的权重与绝对位置无关。

基于平移不变性，可以在ACS区域内滑动卷积核，收集大量的源点-目标点配对作为训练样本。对于第一类核位置，设ACS区域内的源点矩阵为 S_{\text{source}}^{ACS}，对应的目标点向量为 s_{\text{target}}^{ACS}，则有：

这里下标 1,2 表示从线圈1和2的数据中进行校准。由于ACS区域提供了大量的训练样本，而待估计的权重参数数量相对较少，这构成了一个超定线性系统。通过最小二乘求解（伪逆），可以得到最优的权重估计：

其中 (\cdot)^\dagger 表示Moore-Penrose伪逆。平移不变性使得可以在大量源点-目标点对上拟合同一个线性系统，从而获得稳健的权重估计。

GRAPPA重建流程

以一个具体的重建实例来说明完整的GRAPPA流程。考虑 256 \times 256 的T1加权图像，使用 R=3 的加速因子，L=20 个线圈通道，ACS区域包含32条全采样的校准行。

第一步：权重校准

利用ACS区域的全采样数据，按照上述最小二乘方法估计 R-1=2 套权重矩阵 W^{(1)} 和 W^{(2)}。

第二步：k空间填充

使用校准得到的权重，对每个线圈的欠采样k空间数据进行填充。将卷积核在整个k空间上滑动，对每个缺失位置应用相应的权重矩阵，得到每个线圈的完整k空间数据。

第三步：逐线圈图像重建

对每个线圈填充后的k空间数据分别进行逆傅里叶变换（IFFT），得到每个线圈的复数图像 x_\ell。由于每个线圈的敏感度分布不同，这些图像会呈现不同的空间加权特性。

第四步：线圈组合

将所有线圈的图像组合成最终的重建图像。最常用的组合方法是Sum-of-Squares（SoS）：

这种组合方法不需要显式的敏感度图，通过取各线圈图像幅值的平方和再开根号，可以消除敏感度加权效应，得到均匀的最终图像。

GRAPPA与SENSE的对比

方法特性对比

SENSE和GRAPPA是并行成像的两种主流方法，它们在多个方面存在差异。从处理域来看，SENSE在图像域工作，而GRAPPA在k空间域工作。从前置条件来看，SENSE需要显式的线圈敏感度图，通常需要单独的校准扫描或从数据中估计；GRAPPA则需要ACS校准数据，这些数据可以嵌入到正式采集中。从重建目标来看，SENSE直接重建组合后的最终图像；GRAPPA重建每个线圈的k空间数据，再通过后续步骤组合成图像。

共同的SNR损失特性

尽管处理方式不同，GRAPPA和SENSE在信噪比性能上遵循相同的物理规律。两种方法都受到加速因子 R 和几何因子 g 的影响，SNR损失公式相同：

其中 SNR_1 是不加速时的参考信噪比。\sqrt{R} 因子来源于欠采样导致的信号平均次数减少，g 因子来源于线圈几何配置导致的噪声放大。从g-map图像可以看出，SENSE和GRAPPA的g因子空间分布模式相似，都呈现出空间变化的噪声放大特性，某些区域（尤其是图像中心）的g值明显高于边缘区域。

Breuer等人在2009年的研究（MRM 62(3): 739-46）详细分析了两种方法的等价性和g因子特性。实际应用中，GRAPPA因为不需要单独的敏感度校准扫描而更受欢迎，但SENSE在某些特定应用场景下可能具有计算优势。

基于模型的MR图像重建

前面介绍的SENSE和GRAPPA方法虽然原理不同，但都属于针对特定采样模式设计的专用重建算法。现在转向一种更通用的框架——基于模型的图像重建方法。这种方法从物理测量过程出发，建立完整的正向模型，然后通过求解逆问题来重建图像。这种框架的优势在于可以统一处理各种采样模式（笛卡尔或非笛卡尔），并且可以方便地引入各种先验信息和正则化约束。

并行MRI的编码矩阵

在并行MRI中，每个接收线圈采集的信号可以用连续形式的积分方程描述。对于第 \ell 个线圈，其在k空间轨迹 k(t) 位置处采集的信号为：

这个方程的物理含义是：k空间信号等于图像 x(r) 与线圈敏感度 C_\ell(r) 的乘积在整个视场上进行傅里叶变换的结果。指数项 e^{-2\iota\pi (k(t) \cdot r)} 是傅里叶变换的核函数，\iota 表示虚数单位。

将连续模型离散化以便于计算机处理。设视场大小为 \text{FOV} = \text{FOV}_x \times \text{FOV}_y = n_x \Delta x \times n_y \Delta y，图像被离散为 N = n_x \times n_y 个像素。将第 \ell 个线圈在所有 M 个采样时刻的k空间数据堆叠成向量：

离散化后的正向模型可以写成矩阵-向量乘法的形式：

其中各矩阵的定义如下。傅里叶编码矩阵 F 的元素为：

这是一个 M \times N 的复数矩阵，第 (i,j) 个元素表示第 j 个像素位置 r_j 对第 i 个k空间采样点 k(t_i) 的贡献。敏感度矩阵 C_\ell 是一个 N \times N 的对角矩阵：

对角线上的元素是第 \ell 个线圈在各像素位置的敏感度值。

并行MRI的完整测量模型

多线圈数据的整合

实际测量中还存在噪声，因此第 l 个线圈的测量数据可以写成：

其中 y_l \in \mathbb{C}^M 是第 l 个线圈测量到的含噪k空间数据，\varepsilon_l 是该线圈的测量噪声。将所有 L 个线圈的数据纵向堆叠，可以得到统一的矩阵方程：

这里 I_L 是 L \times L 的单位矩阵，\otimes 表示Kronecker积（克罗内克积）。(I_L \otimes F) 的作用是对每个线圈的数据独立应用相同的傅里叶编码。整个系统矩阵可以简记为：

模型中各变量的含义

在这个测量模型中，A \in \mathbb{C}^{ML \times N} 是系统矩阵，包含了傅里叶编码和线圈敏感度的联合作用；y \in \mathbb{C}^{ML} 是所有线圈测量到的k空间数据的堆叠向量；x \in \mathbb{C}^N 是待重建的潜在图像（真实图像）；\varepsilon 是复高斯白噪声。图像重建的目标就是从测量数据 y 中恢复出图像 x。这个模型框架还可以扩展以考虑其他物理效应，如弛豫效应和场不均匀性等。

普通最小二乘重建

超定系统的最小二乘解

当测量数据量大于等于未知像素数量时，即 ML \geq N，线性模型是超定的，可以使用普通最小二乘法求解。对于并行MRI，重建目标是从所有线圈的k空间测量数据中恢复原始图像 x。最小二乘估计定义为使所有线圈数据残差平方和最小的解：

其中 y_\ell \in \mathbb{C}^M 是第 \ell 个线圈测量到的k空间数据，F_\Omega 是在采样位置 \Omega 上的傅里叶编码矩阵，C_\ell 是第 \ell 个线圈的敏感度矩阵。敏感度矩阵是 N \times N 的对角矩阵，其对角元素为该线圈在各像素位置的敏感度值：

将所有线圈的数据和编码整合为统一的矩阵形式，系统矩阵可以写成 A = (I_L \otimes F)C，其中 I_L 是 L \times L 的单位矩阵，\otimes 表示Kronecker积。此时目标函数等价于：

最优性条件的推导

为求解上述优化问题，对目标函数关于 x 求梯度并令其为零。目标函数对 x 的梯度为：

其中 C_\ell^* 表示 C_\ell 的共轭转置（对于对角矩阵就是对角元素取共轭），F_\Omega^H 是傅里叶矩阵的共轭转置（伴随算子）。令梯度为零得到最优解 \hat{x} 满足的条件：

正规方程与解析解

将上述条件展开整理，得到正规方程：

这个方程的形式是标准形式

左边括号内的矩阵 \sum_{\ell=1}^{L} C_\ell^* F_\Omega^H F_\Omega C_\ell 是Gram矩阵，右边 \sum_{\ell=1}^{L} C_\ell^* F_\Omega^H y_\ell 是将各线圈数据通过伴随操作变换回图像域后的敏感度加权求和。当Gram矩阵可逆时，最小二乘解的解析形式为：

这个表达式说明重建过程包含两个步骤：首先对每个线圈的数据进行敏感度加权的傅里叶变换伴随操作 C_\ell^* F_\Omega^H y_\ell，然后将所有线圈的结果求和并通过Gram矩阵求逆进行去卷积。

全采样笛卡尔数据的简化

对于全采样的笛卡尔k空间数据，傅里叶矩阵 F 是酉矩阵（正交矩阵的复数推广），满足

即

利用这个性质，Gram矩阵中的 F^H F 项简化为单位矩阵的倍数，最小二乘解简化为：

这个结果被称为最优线圈组合方法。它先对每个线圈的k空间数据进行逆傅里叶变换 F^{-1}y_\ell 得到线圈图像，然后用敏感度的共轭 C_\ell^* 进行加权，最后通过敏感度Gram矩阵 \sum_{\ell=1}^{L} C_\ell^* C_\ell 的逆进行归一化。这正是Roemer等人在1990年提出的相控阵线圈组合方法的数学形式。

规则欠采样与SENSE的联系

对于规则欠采样的笛卡尔数据（即并行成像中常用的等间隔跳行采样），矩阵 F_\Omega^H F_\Omega 不再是单位矩阵的倍数，但具有简单的块结构。这种块结构反映了欠采样导致的图像混叠：空间上相距 \text{FOV}/R 的像素会在同一个混叠位置叠加。利用这种结构，大型矩阵方程可以分解为多个小型独立方程，每个方程对应一组混叠像素的分离问题。这正是SENSE重建方法的数学本质——它是基于模型重建框架在规则笛卡尔欠采样条件下的特例，通过敏感度编码信息将混叠的像素分离开来。Pruessmann等人1999年的SENSE论文正是基于这一数学框架推导的。

正则化重建方法

正则化的必要性

前面介绍的普通最小二乘方法在欠采样情况下会失效。当测量数据量小于未知像素数量时，即 ML < N，线性系统是欠定的，存在无穷多个解都能完美拟合数据。此外，即使系统不是严格欠定的，当采用非笛卡尔采样模式时，系统矩阵 A 的条件数可能很大，导致问题病态（ill-conditioned），微小的噪声会被放大成巨大的重建误差。正则化技术通过引入关于图像的先验知识来约束解空间，使问题变得适定。

二次正则化

最简单的正则化形式是二次正则化最小二乘（Regularized Least Squares, RLS），其优化目标为：

第一项 \frac{1}{2}\|Ax - y\|_2^2 是数据保真项，衡量重建图像通过正向模型生成的预测数据与实际测量数据之间的差异；第二项 \beta\|Tx\|_2^2 是正则化项，对图像施加平滑性约束。参数 \beta > 0 是正则化参数，控制数据保真和正则化之间的权衡——\beta 越大，解越平滑但可能偏离数据；\beta 越小，解越贴近数据但可能受噪声影响。矩阵 T 是 K \times N 的变换矩阵，常见选择是有限差分算子，此时 \|Tx\|_2^2 衡量图像梯度的能量，惩罚图像中的剧烈变化。

二次正则化的求解

由于目标函数是关于 x 的二次函数，其Hessian矩阵为 A'A + \beta T'T。共轭梯度（Conjugate Gradient, CG）算法非常适合求解这类问题，因为它只需要矩阵-向量乘法而不需要显式构造和存储完整矩阵。当 T 是有限差分算子时，Hessian矩阵 A'A + \beta T'T 近似具有Toeplitz结构（对于笛卡尔采样是精确的），这允许使用循环预条件技术加速CG的收敛。虽然RLS方法在历史上曾广泛使用，但由于二次惩罚的局限性，现在通常只作为更复杂算法的内层迭代步骤。

边缘保持正则化

二次正则化的缺陷

当 T 选择为有限差分算子时，二次正则化 \|Tx\|_2^2 会惩罚图像中所有的梯度，包括真实的边缘结构。这导致重建图像的边缘被模糊，细节丢失。医学图像中的组织边界、病灶轮廓等信息对诊断至关重要，因此需要一种能够保持边缘的正则化方法。

非二次势函数

边缘保持正则化通过将二次惩罚替换为非二次的势函数 \psi 来解决边缘模糊问题：

势函数 \psi 通常选择为凸且光滑的非二次函数。常见的选择包括：Huber函数，它对小梯度使用二次惩罚、对大梯度使用线性惩罚；Hyperbola函数；以及Fair势函数，其定义为：

其中 \delta 是控制从二次行为过渡到线性行为的阈值参数。这些势函数的共同特点是：对于小的梯度值（平坦区域的噪声），惩罚近似二次，起到平滑去噪的作用；对于大的梯度值（真实边缘），惩罚增长较慢（亚线性），不会过度抑制边缘。这种方法的理论根源可以追溯到贝叶斯方法中的马尔可夫随机场（Markov Random Fields）模型。

边缘保持正则化的优化算法

算法选择

由于边缘保持正则化的目标函数是非二次的，不能像RLS那样直接用线性方程组求解，需要使用迭代优化算法。常用的算法包括：非线性共轭梯度（Nonlinear CG）算法，它是线性CG向非二次目标函数的推广；3MG（majorize-minimize memory gradient）算法，采用优化-最小化策略；以及优化梯度法（Optimized Gradient Method, OGM）。

OGM是一种一阶优化方法，专门针对具有Lipschitz连续梯度的凸代价函数设计。它在理论上达到了这类问题的最优最坏情况复杂度，收敛速度界比经典的Nesterov快速梯度法快2倍。还可以结合线搜索技术进一步改进OGM的实际收敛速度。

重建示例

考虑一个具体的重建实验：k空间采样模式为中心区域全采样、外围欠采样的变密度模式。使用Fair势函数作为正则化项，参数 \delta = 0.1；变换矩阵 T 选择为有限差分算子，对应corner-rounded全变分（TV）正则化。

从重建结果可以看到，初始图像（直接逆傅里叶变换）由于欠采样产生了明显的混叠伪影，而经过边缘保持正则化重建后，伪影被有效抑制，同时边缘结构得到保持。最终的归一化均方根误差（NRMSE）达到1.55%。

从代价函数的收敛曲线可以观察到，OGM算法比CG算法收敛更快，在相同迭代次数下达到更低的代价函数值。

本讲内容总结

并行成像技术

并行成像通过笛卡尔欠采样来减少扫描时间，加速倍数等于欠采样因子。其加速能力受限于线圈数量和SNR损失（由g因子决定）。这类方法主要依赖线性重建算法，根据处理域的不同分为两类：SENSE在图像空间处理混叠，GRAPPA在k空间填充缺失数据。在高加速因子下，线圈敏感度的准确校准变得困难。并行成像技术在过去近20年里已经在临床应用中得到广泛使用，带来了革命性的影响。

基于模型的重建方法

当重建问题欠定时，需要引入正则化来约束解空间。二次正则化虽然计算简单，但会模糊图像边缘，现在已经不再推荐使用。边缘保持的凸正则化惩罚项能够在去噪的同时保持图像的边缘结构。针对这类非线性凸优化问题，存在多种高效的迭代算法，它们能够快速收敛到全局最优解。

压缩感知在MRI中的应用

压缩感知的历史背景与核心思想

压缩感知（Compressed Sensing, CS）是21世纪初由Donoho、Tao、Romberg、Candès等数学家提出的信号处理理论，它从根本上改变了我们对信号采样和重建的认识。传统的奈奎斯特采样定理要求采样频率至少是信号最高频率的两倍才能完美重建信号，而压缩感知理论表明：如果信号具有稀疏性或可压缩性，则可以用远少于奈奎斯特率的采样来精确恢复信号。

压缩感知的核心思想可以概括为：数据本身是稀疏的、可压缩的、冗余的，因此可以直接感知压缩后的信息，而不是先采集完整数据再压缩。这一思想建立在三个支柱之上：信号的稀疏性或可压缩性提供了信息冗余的来源；伪随机或非均匀采样模式确保欠采样产生的伪影呈现类噪声的非相干分布；非线性重建算法利用稀疏先验从欠采样数据中恢复原始信号。

压缩感知的基本概念

稀疏性与可压缩性

稀疏信号是指只有少量非零系数的信号。从直观上看，如果一个长度为 N 的向量 x 中只有 K 个元素非零（K \ll N），那么这个向量就是 K-稀疏的。稀疏信号的示意图显示大部分系数为零（黑色），只有少数系数非零（白色亮点）。

基（basis）是用于表示信号的一组基本元素。变换（transform）将信号从一个基变换到另一个基。常见的变换包括傅里叶变换 \mathcal{F}（对应傅里叶基）和小波变换 \psi（对应小波基）。同一个信号在不同基下的表示可能具有完全不同的稀疏程度。

可压缩性是指信号在某个特定的基下是稀疏的。真实世界的图像（包括MRI数据）虽然在原始像素域通常不是稀疏的，但在小波（wavelet）、曲波（curvelet）、脊波（ridgelet）等变换域中具有良好的可压缩性。

以经典的Lena图像为例：在傅里叶域，能量分布相对分散，在图像域，所有像素都有值，但在小波域，大部分系数接近零，只有少数系数具有显著幅值，这正是可压缩性的体现。

稀疏性的度量范数

衡量稀疏性需要合适的数学工具。三种常用的范数定义如下：

\ell_0 伪范数（pseudo-norm）直接计算非零系数的数量：

\ell_0 伪范数值越小，表示信号越稀疏。它是稀疏性最直接的度量，但在数学上不是真正的范数（不满足齐次性），且 \ell_0 最小化问题是NP难的组合优化问题。

\ell_1 范数是信号系数绝对值之和：

\ell_1 范数是 \ell_0 伪范数的凸松弛（convex relaxation），即稀疏性的凸代理。虽然它不直接计数非零元素，但最小化 \ell_1 范数往往能得到稀疏解，这一性质使得压缩感知问题变得计算可行。

\ell_2 范数是信号系数平方和的平方根：

\ell_2 范数衡量的是信号的能量。\ell_2 范数越小表示能量越低，但它与稀疏性没有直接关系。

范数单位球的几何直观

\ell_p 范数单位球的定义是满足 \|x\|_p = 1 的所有点构成的集合。不同 p 值对应的单位球形状差异很大，这种几何差异直接影响了正则化的效果。

在二维空间中：p=2 时单位球是圆形，p=1 时单位球是菱形（正方形旋转45度），p=0 时单位球退化为坐标轴上的点。在三维空间中：p=2 时单位球是球体，p=1 时单位球是八面体，p=0 时单位球是坐标轴上的六个点。

随着 p 从2减小到0，单位球变得越来越尖锐，在坐标轴方向上越来越突出。

\ell_1 正则化诱导稀疏解的几何解释

理解 \ell_1 正则化为何能产生稀疏解，需要从优化问题的几何角度来分析。考虑最小二乘问题加上正则化约束：

其中 f(x) 是数据拟合项（如 \|Ax-y\|_2^2），约束条件限制解必须落在 \ell_p 范数球内。

从几何上看，f(x) 的等高线是一系列嵌套的椭圆（对于二次函数），最优解 x^* 是使 f(x) 最小的等高线与约束区域（范数球）相切的点。

对于 \ell_2 正则化（p=2），约束区域是圆形。等高线椭圆与圆相切的点通常位于圆周的任意位置，不一定在坐标轴上，因此解的各分量一般都非零。

对于 \ell_1 正则化（p=1），约束区域是菱形。由于菱形的顶点位于坐标轴上，等高线椭圆与菱形相切时，切点很可能落在菱形的顶点处。顶点处的坐标有一个或多个分量为零，这正是稀疏解的特征。即使最优切点不在顶点，也更可能落在菱形的边上（某些分量为零），而不是内部。

这种几何直观解释了压缩感知为何选择 \ell_1 范数作为稀疏性的代理：它是最接近 \ell_0 的凸范数，既保留了诱导稀疏解的能力，又使得优化问题是凸的，可以高效求解全局最优。

压缩感知MRI的整体框架

压缩感知应用于MRI需要三个关键要素的配合：稀疏表示域、非相干测量和稀疏促进的非线性重建。

稀疏域与可压缩性

MR图像在小波域（wavelet domain）具有良好的可压缩性。设原始图像为 x，小波变换算子为 \Psi，则图像的小波系数向量为：

从示例图像可以看到，脑部MR图像经过小波变换后，大部分系数接近零（图中暗色区域），只有少数系数具有显著幅值（亮色区域），这正是稀疏性的体现。

测量域与欠采样

MRI的测量过程是在k空间（傅里叶域）进行采样。压缩感知MRI采用伪随机的欠采样模式，而不是规则的等间隔跳行采样。变密度采样（Variable Density Sampling）是一种常用的策略：在k空间中心（低频区域）采样密度高，在外围（高频区域）采样密度低。这种设计基于MR图像能量主要集中在低频的特点。

非相干性条件

压缩感知成功的关键条件之一是测量域与稀疏域之间的非相干性（incoherence）。在MRI中，测量变换是傅里叶变换 \mathcal{F}，稀疏化变换是小波变换 \Psi。傅里叶基和小波基之间具有天然的低相干性：傅里叶基函数是全局振荡的正弦波，而小波基函数是局部化的。这种非相干性确保了当在k空间欠采样时，产生的混叠伪影在小波域呈现类似噪声的非结构化分布，从而可以通过稀疏约束将其与真实信号分离。

优化重建模型

设 \mathcal{F}_K 表示在采样位置集合 K 上的傅里叶变换（欠采样傅里叶算子），y 为测量数据，\epsilon 为噪声，则测量模型为：

压缩感知MRI的重建通过求解以下优化问题实现：

这个目标函数包含两项：第一项 \frac{1}{2}\|\mathcal{F}_K x - y\|_2^2 是数据一致性项（Data Consistency），确保重建图像在已采样的k空间位置与测量数据一致；第二项 \lambda\|\Psi x\|_1 是稀疏性约束项（Sparsity Constraint），通过 \ell_1 范数促进小波系数的稀疏性。参数 \lambda > 0 控制两项之间的权衡。

压缩感知重建的迭代算法

交替优化策略

上述优化问题可以通过迭代方法求解，每次迭代包含两个步骤：数据一致性步骤和稀疏性强化步骤。

数据一致性步骤对数据保真项进行梯度下降：

其中 \eta 是步长参数，\mathcal{F}_K^* 是 \mathcal{F}_K 的伴随算子（即在采样位置填充数据、其余位置补零后进行逆傅里叶变换）。这一步计算当前估计图像在k空间的预测值 \mathcal{F}_K x^t 与实际测量值 y 的残差，然后将残差变换回图像域并用于更新图像估计。其物理意义是确保重建图像通过正向模型得到的k空间数据与实际采集的数据相匹配。

稀疏性强化步骤在小波域进行软阈值操作：

这一步首先将图像变换到小波域 \Psi x^{t+\frac{1}{2}}，对小波系数应用软阈值函数 \text{soft}_{\lambda\eta}(\cdot)，然后通过逆小波变换 \Psi^* 返回图像域。这个过程可以理解为在小波域进行去噪，将小幅值系数（很可能是噪声或伪影）压缩为零，保留大幅值系数（真实信号成分）。

阈值函数

软阈值（soft thresholding）函数的定义为：

对于幅值小于阈值 \lambda 的系数，软阈值将其置为零；对于幅值大于 \lambda 的系数，软阈值将其幅值减去 \lambda。软阈值函数是连续的，这使得整个算法的收敛性分析更加容易。

与之对比，硬阈值（hard thresholding）函数直接将小于阈值的系数置零，大于阈值的系数保持不变：

硬阈值虽然更直接地反映稀疏性，但函数不连续，在优化中可能导致收敛问题。\ell_1 正则化对应的近端算子恰好是软阈值函数，这是选择软阈值的数学原因。

迭代重建的流程图解

迭代重建算法可以用流程图清晰地展示。从观测到的k空间数据 y（只在采样位置 K 有值）开始，算法在图像域和变换域之间交替操作。

在数据一致性分支：首先计算当前图像估计 x^t 的k空间表示 \mathcal{F}_K x^t，然后与测量数据 y 相减得到k空间残差 \mathcal{F}_K x^t - y，再通过伴随算子 \mathcal{F}_K^* 将残差变换回图像域。这个残差图像显示了当前估计与数据不一致的区域。

在稀疏性强化分支：将更新后的图像 x^{t+\frac{1}{2}} 变换到小波域，得到小波系数 \Psi x^{t+\frac{1}{2}}。对这些系数应用软阈值操作 \text{soft}_{\lambda\eta}(\cdot)，小幅值系数被压缩。最后通过逆小波变换 \Psi^* 返回图像域，得到新的图像估计 x^{t+1}。

这两个步骤交替进行，直到收敛。每次迭代都同时改善数据一致性和稀疏性，最终得到既符合测量数据又具有稀疏小波表示的重建图像。

采样模式与非相干性

非相干性的本质

压缩感知理论要求测量矩阵与稀疏基之间具有非相干性。在MRI的语境下，这意味着傅里叶采样模式（加上欠采样）与小波稀疏表示之间应该是非相干的。

从示意图可以看到：图像 x 通过测量变换（傅里叶变换 \mathcal{F}）得到k空间数据，结合采样模式形成实际的测量；同时，图像通过稀疏化变换（小波变换 \Psi，通常是数据相关的）得到稀疏表示。非相干性条件要求这两个域之间的关联性低，即测量域加采样模式与稀疏域之间应该不相干。

直观理解非相干性：如果在k空间的某个位置欠采样，产生的混叠伪影在图像域会呈现某种模式。非相干性确保这种伪影模式在稀疏域（小波域）不会与真实信号的稀疏模式相似，而是分散成类似噪声的非结构化分布。这样，稀疏约束就能有效区分真实信号和欠采样伪影。

傅里叶-小波系统天然具有较好的非相干性，这是压缩感知MRI得以成功的物理基础之一。通过设计合适的欠采样模式（如伪随机采样、变密度采样），可以进一步优化非相干性，提高重建质量。

非相干性的量化定义

变换点扩散函数

非相干性可以通过变换点扩散函数（Transform Point Spread Function, TPSF）来定量描述。TPSF的计算过程如下：首先在稀疏域中取一个单点信号（即只有一个位置非零的稀疏向量）；然后将这个点变换到k空间，进行欠采样操作；最后将欠采样后的数据变换回稀疏域。TPSF就是这个过程的输出，它描述了一个稀疏域中的单点信号经过欠采样测量系统后，在稀疏域中的响应分布。

理想情况下，如果测量系统完美保持信息，TPSF应该仍然是一个单点。但由于欠采样，TPSF会出现旁瓣（sidelobe），表现为原始点位置之外的残余响应。这些旁瓣就是欠采样引入的伪影在稀疏域的表现。非相干性水平（coherence level）定义为TPSF中最大旁瓣的幅值——旁瓣越小，非相干性越好，压缩感知重建的效果越好。

这里需要强调的是，非相干性的要求是针对测量域和稀疏域之间的关系，而不一定涉及图像域。在MRI中，测量域是k空间（傅里叶域），稀疏域是小波域，图像域只是中间表示。

采样模式的比较

不同轨迹的特性

不同的k空间采样模式在非相干性和梯度效率方面有不同的表现。

从示例图像可以比较六种采样模式：2DFT（标准笛卡尔采样）、Random（随机采样）、Radial（径向采样）、Spiral（螺旋采样）、V.D. Spiral（变密度螺旋）和Perturbed Spiral（扰动螺旋）。

2DFT是传统的逐行扫描模式，欠采样时会产生规则的混叠伪影，非相干性较差。随机采样具有最好的非相干性，欠采样伪影呈现类噪声分布，但在实际MRI中难以实现，因为k空间位置由梯度积分决定，随机跳跃需要极快的梯度切换。径向和螺旋采样是实际可行的非笛卡尔轨迹，它们在欠采样时产生的伪影比笛卡尔采样更分散。变密度和扰动螺旋进一步优化了采样分布，提高了非相干性。

采样模式的设计需要平衡非相干性和梯度效率。非相干性越高，压缩感知重建效果越好；但过于复杂的轨迹可能超出梯度硬件的能力限制，或导致采集效率降低。变密度采样模式可以利用k空间能量分布的先验知识（中心能量高、外围能量低），在能量集中的区域增加采样密度。

变密度采样的理论基础

均匀随机采样的局限

考虑对MR图像使用5%的均匀随机采样进行压缩感知重建。结果显示重建质量很差，图像严重失真。这是因为MR图像的小波分解具有特殊的结构：低频的近似系数（approximation coefficients）包含图像的主要结构信息，能量集中；高频的细节系数（detail coefficients）是稀疏的，能量分散。

小波分解图中，左上角的近似系数对应图像的整体亮度和对比度，这部分信息在k空间中心（低频区域）。如果均匀随机采样，k空间中心的采样密度与外围相同，就无法充分捕获低频的近似信息。近似系数无法通过稀疏约束恢复，因为它们本身不是稀疏的，必须通过实际采样获得。相反，细节系数是稀疏的，可以通过稀疏约束从欠采样数据中恢复。

变密度采样的优势

变密度采样（Variable Density Sampling, VDS）根据k空间能量分布调整采样密度：在能量集中的中心区域（对应小波近似系数）增加采样密度，在能量分散的外围区域（对应小波细节系数）降低采样密度。这样既保证了低频信息的完整采集，又通过高频区域的欠采样实现加速，同时保持良好的非相干性。

使用5%的变密度随机采样进行压缩感知重建，结果显示图像质量明显优于均匀采样。

变密度采样的实际实现

MRI中的k空间轨迹

在MRI中，k空间位置由梯度场的时间积分决定：

其中 G(\tau) 是梯度波形。这意味着k空间轨迹是连续的曲线，不能任意跳跃到随机位置。

分段采集（segmented acquisition）将整个k空间轨迹分成多个片段（shots），每个射频激发后采集一段轨迹。扫描时间与片段数量成正比。常见的非笛卡尔轨迹包括径向（从中心向外的直线）和螺旋（从中心向外盘旋）。

梯度硬件约束

梯度系统存在物理限制：最大梯度幅值 G_{\max} < 40 \text{ mT/m}，最大梯度变化率（slew rate）S_{\max} < 200 \text{ T/m/s}。这些限制约束了k空间轨迹的速度和加速度，从而限制了可实现的采样模式。设计实用的变密度采样轨迹需要在满足硬件约束的前提下优化非相干性和采样效率。

压缩感知MRI的加速极限

实验研究结果

压缩感知在MRI中能够实现多大的加速？Lazarus等人（MRI 2018）针对T2*成像进行了系统的实验研究。他们使用结构相似性指数（Structural Similarity Index, SSIM）作为图像质量度量，SSIM的定义为：

其中 l、c、s 分别表示亮度、对比度和结构相似性分量，\alpha、\beta、\gamma 是权重参数。SSIM > 0.9通常被认为是可接受的图像质量。

实验结果显示，最大可实现的加速因子 R_{\max} 取决于图像尺寸和输入信噪比。对于256×256的图像，在高SNR条件下，可以达到约10倍的加速（SSIM > 0.9）。更大的图像尺寸允许更高的加速因子，因为图像越大，稀疏表示的冗余度越高。输入SNR越高，能够容忍的加速因子也越高。从右侧的曲线图可以看出，当图像尺寸达到1024时，在一定SNR条件下可以达到15倍以上的加速。

压缩感知与并行成像的结合

联合重建框架

压缩感知可以与并行成像技术结合，充分利用两者的优势。并行成像利用多线圈的空间编码冗余，压缩感知利用图像的稀疏性先验。联合重建的优化问题为：

其中 \mathcal{F}_\Omega 是在采样位置 \Omega 上的傅里叶变换（非均匀采样），S_\ell 是第 \ell 个线圈的敏感度图（对角矩阵形式），y_\ell 是第 \ell 个线圈的k空间测量数据，f(\Psi x) 是稀疏促进正则化项（如 \ell_1 范数或Group LASSO）。

实现流程

联合重建的实现涉及多个模块：非均匀采样提供k空间数据；敏感度图描述各线圈的空间响应；小波分解实现稀疏表示；稀疏促进算法（如 \ell_1 范数、Group LASSO）强化稀疏性约束；非线性迭代求解器（如FISTA、POGM）求解优化问题。

开源软件包如PySAP（Python Sparse data Analysis Package）和ModOpt提供了这些算法的实现，可以直接用于压缩感知MRI重建研究和应用。

本讲内容总结

压缩感知采集

压缩感知MRI目前通常以较低的加速率实现，这是因为欠采样性能受到多种因素限制。变密度采样不仅是经验上有效的策略，更是理论上必需的——它是突破相干性壁垒、实现高质量重建的关键。非笛卡尔轨迹（如SPARKLING）在采样效率和非相干性方面优于笛卡尔采样。压缩感知可以与并行成像技术结合使用，进一步提高加速能力。

压缩感知重建

重建图像的质量取决于稀疏化基的选择——图像在该基下越稀疏，重建效果越好。图像重建是迭代过程，计算量较大，需要在数据一致性和稀疏性之间反复优化。正则化参数（如 \lambda）需要仔细调节，以平衡数据保真和稀疏约束。目前已有多种开源软件包提供压缩感知MRI重建算法的实现，便于研究和应用。

SPARKLING轨迹优化方法

SPARKLING（Spreading Projection Algorithm for Rapid K-space sampLING，快速k空间采样的扩展投影算法）是一种基于优化的非笛卡尔k空间轨迹设计方法。该方法由Boyer、Lazarus、Chauffert等研究者在一系列工作中提出和发展，相关论文包括：Boyer等人2016年发表在SIAM Imaging Science上的采样方案生成方法；Lazarus等人2019年发表在Magnetic Resonance in Medicine上的用于T2*加权MRI加速的变密度k空间填充曲线；Lazarus等人2020年发表在NMR in Biomedicine上的3D变密度SPARKLING轨迹；以及Chaithya等人2022年发表在IEEE Transactions on Medical Imaging上的全3D SPARKLING轨迹优化方法。

轨迹设计的核心思想

SPARKLING的设计目标是生成满足MRI硬件约束的k空间采样轨迹，同时使采样点的分布尽可能接近给定的目标采样密度（通常是变密度分布）。

整个设计流程涉及多个要素：MRI图像通过傅里叶变换对应k空间；初始轨迹可以选择径向或螺旋等经典非笛卡尔模式；目标采样密度 \pi 通常是变密度分布（中心密度高、外围密度低）；扫描仪的梯度硬件存在幅值和变化率的限制。SPARKLING算法通过投影梯度下降（Projected Gradient Descent）方法，在满足硬件约束的前提下优化轨迹，使最终的采样点分布逼近目标密度。

距离最小化的优化目标

SPARKLING的数学表述是一个约束优化问题，目标是最小化实际采样密度与目标采样密度之间的距离。设 k = (k[i])_{1 \leq i \leq N} \in \mathcal{Q}_N 表示由 N 个采样点组成的轨迹，实际的采样密度可以表示为这些点位置上的狄拉克函数之和：

其中 \delta_{k[i]} 是位于 k[i] 处的狄拉克δ函数。优化问题的形式为：

这里 h 是一个平滑核函数，\star 表示卷积操作，\pi 是目标采样密度。通过与平滑核卷积，将离散的采样点分布和目标密度都转换为连续函数，从而可以计算它们之间的 \ell_2 距离。

硬件约束集的定义

约束集 \mathcal{Q}_N 包含了所有满足MRI硬件限制的合法轨迹。对于离散化的轨迹，约束通过有限差分来近似连续的导数条件：

其中轨迹的速度和加速度通过有限差分计算：

\ell_{2,\infty} 范数定义为所有时间点上向量2-范数的最大值：

速度约束 \|\dot{k}\|_{2,\infty} \leq \gamma G_{\max} 对应梯度幅值限制（G_{\max} 约为40 mT/m），加速度约束 \|\ddot{k}\|_{2,\infty} \leq \gamma S_{\max} 对应梯度变化率（slew rate）限制（S_{\max} 约为200 T/m/s）。\gamma 是旋磁比。

SPARKLING的优化驱动表述

目标函数的分解

SPARKLING的优化目标可以分解为两个竞争的项。最优轨迹 \hat{K} 通过最小化以下目标函数得到：

其中 F_N^a(K) 是吸引项（attraction term），F_N^r(K) 是排斥项（repulsion term）。

吸引项的定义为：

这一项衡量采样点与目标密度之间的匹配程度。H 是一个核函数，\pi(x) 是目标采样密度。吸引项的作用是将采样点拉向目标密度高的区域——当采样点靠近高密度区域时，积分值增大，因此最小化时会倾向于让采样点分布与目标密度一致。

排斥项的定义为：

这一项衡量采样点之间的相互作用。H(K[n] - K[m]) 度量第 n 个和第 m 个采样点之间的邻近程度。排斥项的作用是防止采样点聚集在一起——当两个采样点靠得很近时，H 值增大，因此最小化时会倾向于让采样点彼此远离，实现局部均匀分布。

吸引项和排斥项的竞争达到平衡时，采样点既遵循目标密度的整体分布（在密度高的地方点多），又保持局部的均匀性（相邻点之间不会过于靠近）。

SPARKLING的约束优化求解

多曲线轨迹表示

实际的SPARKLING轨迹由多条曲线组成，每条曲线对应一次射频激发后的k空间读出。设轨迹集合为 K = [k_i]_{i=1}^{N_c}，其中 N_c 是曲线（shots）的数量。优化问题变为：

约束集的完整定义

对于多曲线轨迹，约束集 \mathcal{Q} 对每条曲线 k_i（i \in \{1, 2, \ldots, N_c\}）施加以下条件：

前两个约束是梯度硬件限制，\alpha 和 \beta 分别是速度和加速度的上界。第三个约束 k_i[TE] = [0, 0, 0]^T 要求每条曲线在回波时间（Echo Time, TE）时刻经过k空间中心，这是MRI信号采集的物理要求——回波信号在 k=0 处最强。

投影梯度下降算法

SPARKLING使用投影梯度下降方法求解约束优化问题。迭代更新公式为：

每次迭代包含两个步骤：首先沿目标函数的负梯度方向更新轨迹 K^t - \eta^t \nabla F_N(K^t, \pi)，其中 \eta^t 是步长；然后通过投影算子 \Pi_\mathcal{Q} 将更新后的轨迹投影回约束集 \mathcal{Q}。

投影算子 \Pi_\mathcal{Q} 的作用是：如果梯度更新导致轨迹违反硬件约束（速度或加速度超限），投影操作会将轨迹调整到最接近的满足约束的位置。这确保了优化过程中始终保持轨迹的物理可行性。

从示意图可以看到优化过程的效果：初始的轨迹（如径向或螺旋）经过投影梯度下降迭代后，逐渐变形为复杂的曲线结构，既满足硬件约束（每条曲线平滑可行），又使整体采样点分布逼近目标密度（中心密集、外围稀疏）。最终的SPARKLING轨迹呈现出独特的扭曲螺旋形态，这种形态是优化算法在吸引、排斥和硬件约束三者之间取得平衡的结果。

SPARKLING的脉冲序列实现

T2*加权成像的序列设计

SPARKLING轨迹的实际应用通常基于T2加权的梯度回波（Gradient-Recalled Echo, GRE）序列，具体实现采用FLASH（Fast Low Angle SHot）序列结构。T2加权成像对磁敏感性差异敏感，能够提供丰富的组织对比度，特别适合检测静脉结构、出血和铁沉积等。

序列实现在7 Tesla高场强西门子扫描仪上进行。脉冲序列时序图显示了各通道的波形：RF通道发射小翻转角（\alpha = 25°）的激发脉冲；Slice通道施加层选梯度；RO（读出）和FE（频率编码）通道在ADC采集窗口期间执行SPARKLING轨迹对应的梯度波形；PE（相位编码）通道在传统笛卡尔序列中用于相位编码，但在SPARKLING中被整合到非笛卡尔轨迹设计中。ADC采集窗口的观测时间 T_{obs} = 30 ms，从时刻 t_0 开始，在回波时间TE结束。

序列参数配置

具体的序列参数包括：矩阵大小 N = 512；视场 \text{FOV} = 200 \times 200 \text{ mm}^2；重复时间 TR = 550 ms，采用交错采集模式；回波时间 TE = 30 ms；翻转角 \alpha = 25°；带宽 BW = 32.55 Hz/pixel；观测时间 T_{obs} = 30.72 ms；层厚3 mm；使用32通道接收线圈；梯度时间分辨率10 µs。

高场强（7T）相比低场强（1.5T）能够提供更高的信噪比和更好的空间分辨率，但同时也带来更强的磁敏感性伪影。T2*加权成像在7T下能够获得更丰富的磁敏感对比度，这对神经影像研究（如检测微小静脉结构）特别有价值。

2D SPARKLING的活体成像结果

在7 Tesla下进行的活体脑部成像实验验证了SPARKLING轨迹的优越性能。成像参数为 0.39 \times 0.39 \times 3 \text{ mm}^3 分辨率，采集11个层面。

参考图像使用全采样采集，扫描时间为4分42秒（TA = 4'42"）。三种加速采集方法都使用相同的加速因子 AF = 20，扫描时间缩短至14秒（TA = 14"）：In-out SPARKLING、In-out Radial（径向）和In-out Spiral（螺旋）。这里的In-out指的是轨迹从k空间中心向外再返回的双向读出模式。

从重建图像的对比可以看出，在相同的20倍加速条件下，SPARKLING轨迹提供了最接近参考图像的重建质量。放大的细节区域显示，径向和螺旋轨迹的重建图像存在明显的条纹伪影和细节模糊（图中箭头标注），而SPARKLING的重建保持了更清晰的解剖结构和更少的伪影。这验证了基于优化的变密度轨迹设计相比传统非笛卡尔轨迹的优势。

向3D SPARKLING的扩展

3D成像的优势与挑战

从2D向3D扩展SPARKLING具有信噪比方面的天然优势。信噪比与激发体积成正比：SNR \propto \text{Excited Volume}。2D多层采集中每层单独激发，激发体积小；而3D单slab采集一次激发整个体积，获得更高的SNR。

从示意图可以看到，2D多层采集的采样覆盖是平行的多个平面，而3D single-slab采集覆盖整个三维k空间体积。

然而，3D SPARKLING的优化面临严峻的计算挑战。优化目标函数仍然包含吸引项和排斥项：

约束集 \mathcal{Q} 的定义与2D情况类似，对每条曲线施加速度、加速度和过中心点约束：

计算复杂度分析

排斥项的计算涉及所有采样点对之间的相互作用，复杂度为 O(N^2)。投影操作 \Pi_\mathcal{Q} 对每条曲线独立进行，复杂度为 O(N_c)。整体的投影梯度下降迭代为：

对于2D SPARKLING，采样点数 N \approx 10^5，曲线数 16 \leq N_c \leq 64，优化时间约为1天（T_c = 1 day），这在实际应用中是可接受的。

对于3D SPARKLING，采样点数剧增至 N \approx 10^7，曲线数 4000 \leq N_c \leq 8000。如果直接使用朴素算法，由于 O(N^2) 的复杂度，优化时间将达到数月甚至数年（T_c = Months to Years），这在实际中是不可行的。

3D轨迹的初始化策略

为缓解计算负担，可以采用不同的初始化和分解策略。Stacks方法将3D k空间分成多个2D层面，分别优化后堆叠，优化时间约1天。Parcellations方法将k空间分成多个区块分别优化，优化时间约3-4天。这些方法虽然减少了计算时间，但由于分块边界处理的问题，可能无法获得全局最优的轨迹。

全3D SPARKLING的加速算法

N体问题与快速多极方法

排斥项的计算本质上是一个N体问题（N-Body problem）：计算N个粒子（采样点）之间的两两相互作用。直接计算的复杂度为 O(N^2)，对于 N = 10^7 量级的3D问题完全不可行。

快速多极方法（Fast Multipole Method, FMM）是解决N体问题的经典加速技术，可以将复杂度从 O(N^2) 降低到 O(N \log N)。FMM的核心思想是利用空间层次结构：对于距离较远的粒子群，可以用多极展开来近似它们的集体效应，而不需要计算每对粒子的精确相互作用。

FMM的层次分解

FMM使用树形数据结构（如八叉树）将空间递归分割。算法包含多种操作类型：P2M（Particle to Multipole）将粒子的影响聚合到多极展开；M2M（Multipole to Multipole）在树的上行阶段合并多极展开；M2L（Multipole to Local）将远场多极展开转换为局部展开；L2L（Local to Local）在树的下行阶段分发局部展开；L2P（Local to Particle）将局部展开应用到粒子；M2P（Multipole to Particle）处理中等距离的相互作用；P2P（Particle to Particle）处理近邻粒子的直接相互作用。

通过这种层次分解，远程相互作用通过低阶多极展开高效近似，只有近邻相互作用需要精确计算，从而实现整体的 O(N \log N) 复杂度。

GPU并行加速

除了算法层面的加速，SPARKLING优化还利用GPU进行并行计算。投影操作 \Pi_\mathcal{Q}(K) = [\Pi_\mathcal{Q}(k_i)]_{i=1}^{N_c} 可以对每条曲线独立进行，具有天然的并行性。FMM算法本身也具有良好的并行结构，适合GPU加速。

PBBFMM3D（Parallel Black-Box Fast Multipole Method）是Wang等人开发的GPU加速FMM实现。性能测试显示，对于 N \times 10^6 量级的问题，PBBFMM相比朴素GPU实现有数量级的速度提升。并行实现的对比显示，GPU实现相比单核CPU有约10倍的加速，相比16核CPU并行也有显著优势。

通过结合FMM算法和GPU并行化，全3D SPARKLING的优化时间可以从数月缩短到约4小时（T_c = 4 Hrs），使得高分辨率3D SPARKLING轨迹的实际应用成为可能。

轨迹质量的点扩散函数分析

点扩散函数的定义与意义

点扩散函数（Point Spread Function, PSF）是评估成像系统空间分辨率的标准工具。PSF定义为：对于一个假想的点源对象（理想的单点信号），通过MRI采集和重建后得到的图像。理想情况下，完美的成像系统应该将点源重建为一个点；但由于采样有限、k空间覆盖不完整等原因，实际重建的PSF会呈现一定的展宽和旁瓣。

PSF越尖锐（peakier），表示成像系统的空间分辨率越好。PSF的质量可以通过以下定量指标来衡量：半高全宽（Full-Width-at-Half-Maximum, FWHM）——PSF主峰宽度在峰值一半处的宽度，FWHM越小表示分辨率越高；峰值与旁瓣比（Peak to Side lobe level）——主峰幅值与最大旁瓣幅值之比，比值越高表示旁瓣抑制越好；峰值与噪声比（Peak to Noise level）——主峰幅值与背景噪声水平之比，比值越高表示信噪比越好。

Full 3D SPARKLING与SpSOS SPARKLING的对比

通过PSF分析可以定量比较不同轨迹设计方法的性能。SpSOS SPARKLING是一种基于堆叠（Stacks）策略的简化3D轨迹，将k空间分层分别优化2D SPARKLING再堆叠。Full 3D SPARKLING则是完全在三维空间中优化的轨迹。

从PSF图像可以看到，Full 3D SPARKLING在三个正交视图（y、z方向）中都呈现更尖锐的中心峰和更低的背景。SpSOS SPARKLING的PSF在某些方向上出现明显的条纹结构，这是分层优化边界效应的体现。

PSF剖面曲线的对比更加直观：Full 3D SPARKLING（蓝色曲线）具有更高的主峰、更低的旁瓣和更低的噪声底。SpSOS SPARKLING（青色曲线）虽然主峰位置相同，但旁瓣明显更高，噪声水平也更高。这说明完全3D优化的轨迹在成像质量上优于分层堆叠的简化方法。

非笛卡尔重建中的密度补偿

非笛卡尔采样模式（如SPARKLING）的k空间采样点分布是不均匀的。在变密度采样中，中心区域采样密度高，外围区域采样密度低。如果在重建时不考虑这种密度差异，会导致图像对比度失真——高密度区域的信息被过度加权。

密度补偿（Density Compensation）通过对k空间数据进行加权来校正采样密度的不均匀性。设 D_K 为采样密度图，迭代重建算法的数据一致性步骤需要修改为包含密度补偿：

在计算k空间残差 \mathcal{F}_K x^t - y 后，将残差与密度补偿因子结合，然后通过伴随算子变换回图像域。密度补偿的作用是给低密度区域的数据更高的权重，给高密度区域的数据较低的权重，从而平衡不同区域的贡献。稀疏性强化步骤保持不变：将图像变换到小波域，应用软阈值，再变换回图像域。

全3D SPARKLING在3T下的成像结果

10倍加速的前瞻性验证

在3 Tesla临床扫描仪上进行了全3D SPARKLING的前瞻性（prospective）验证实验。成像参数为600 µm各向同性分辨率的T2*加权成像，加速因子 AF = 10。

三种方法的对比结果：笛卡尔参考扫描需要15分13秒；Full 3D SPARKLING和SpSOS SPARKLING都只需要4分58秒，实现了约3倍的时间节省。从多个切面的重建图像可以看到，Full 3D SPARKLING的图像质量与笛卡尔参考最为接近，细节结构清晰可辨。SpSOS SPARKLING虽然也能获得可用的图像，但在某些区域（如放大框中显示的区域）存在轻微的伪影。

更高加速因子的探索

进一步测试了不同加速因子下Full 3D SPARKLING和SpSOS SPARKLING的性能。在600 µm各向同性分辨率的T2*加权成像条件下：AF = 10 时扫描时间为4分58秒；AF = 15 时扫描时间为3分22秒；AF = 20 时扫描时间为2分36秒；AF = 40 时扫描时间仅为1分16秒。

从重建图像可以观察到，随着加速因子增大，图像质量逐渐下降，但Full 3D SPARKLING始终优于SpSOS SPARKLING。即使在40倍加速下，Full 3D SPARKLING仍然保持了可识别的解剖结构，而SpSOS SPARKLING出现了更明显的噪声和伪影。作为对比，笛卡尔参考（图e）需要15分13秒才能获得高质量图像。

SPARKLING的应用领域与影响

磁敏感加权成像

T2* GRE（梯度回波）序列结合SPARKLING可以实现快速高分辨率的磁敏感加权成像（SWI）。对比显示，笛卡尔采集需要15分13秒，而SPARKLING仅需4分58秒即可获得相当质量的SWI图像。SWI对检测微出血、静脉结构、钙化等具有独特价值，扫描时间的缩短有助于提高临床可行性和患者舒适度。

功能磁共振成像

SPARKLING轨迹可以通过重复扫描应用于功能磁共振成像（fMRI）。fMRI需要快速采集大量时间序列数据以捕捉血氧水平依赖（BOLD）信号的动态变化。Amor等人（MRM 2023）研究了场不均匀性校正对3D-SPARKLING fMRI中BOLD敏感性的影响。

对比3D-SPARKLING和传统3D-EPI的激活检测结果（p < 0.001显著性水平），两种方法都能检测到感兴趣效应，但3D-SPARKLING在某些区域可能提供不同的敏感性分布。这为fMRI研究提供了新的采集策略选择。

钠离子成像

钠（^{23}Na）MRI由于信号极低、弛豫时间极短，对采集效率的要求极高。Baptista等人（Journal of Neuroradiology 2023）将SPARKLING应用于快速脑部钠离子成像。与传统的TPI（Twisted Projection Imaging）轨迹相比，SPARKLING能够在28秒的短扫描时间内获得可用的钠离子图像，显著提高了钠MRI的临床实用性。

血管成像与灌注成像

牛津大学FMRIB中心的Q. Shen和T. Okell等研究者将SPARKLING应用于血管成像（Angiography）和灌注成像（Perfusion imaging）。对比结果显示，SPARKLING相比现有最先进方法（SOTA）能够提供相当或更好的血管可视化效果，这对脑血管疾病的诊断具有潜在价值。

未来发展方向

SPARKLING技术正在向解剖成像领域扩展，特别是与MP2RAGE序列的结合。MP2RAGE是一种用于高分辨率T1加权解剖成像的序列，广泛应用于脑结构分析和分割。将SPARKLING与MP2RAGE结合有望在保持解剖图像质量的同时显著缩短扫描时间。

离共振效应的物理背景

B0场不均匀性的来源

MRI成像依赖于均匀的主磁场 B_0，但实际中由于多种因素，B_0 场会存在空间不均匀性。这些因素包括：磁体本身的制造和匀场限制；被扫描对象引入的磁敏感性差异，特别是空气-组织界面（如鼻窦、耳道附近）会产生显著的局部场畸变；以及涡流效应等。场不均匀性用 \Delta B_0(r) 表示，它是空间位置 r 的函数。

从场图（field map）可以看到，\Delta B_0 在大脑不同区域的分布范围可以从 -200 Hz 到 +200 Hz。这种场不均匀性在气-组织交界区域（如额叶底部、颞叶）尤其明显。

离共振对信号的影响

B_0 场不均匀性会导致局部共振频率偏移。定义离共振频率为：

其中 \gamma 是旋磁比。这个频率偏移会在信号采集过程中引入额外的相位累积。考虑离共振效应的连续信号方程为：

与理想的傅里叶编码相比，多了一个时间依赖的相位项 e^{-j\omega(r)t}。这个相位项随采集时间 t 累积，导致信号偏离理想的傅里叶关系。

离散化后的信号方程为：

笛卡尔与非笛卡尔采样的不同影响

离共振效应对笛卡尔和非笛卡尔采样的影响程度不同。在笛卡尔逐行采样中，每行的采集时间很短，相位累积有限，主要表现为图像的几何畸变和信号丢失（在离共振严重区域）。在非笛卡尔采样（如SPARKLING）中，单次读出时间较长，相位累积更显著，会产生更复杂的伪影模式，包括模糊和几何失真。

从示例图像可以看到，相同的模体在笛卡尔和SPARKLING轨迹下的重建结果：不考虑 \Delta B_0 时两者都有伪影，但SPARKLING的伪影更为明显；加入 \Delta B_0 校正后，两者的图像质量都得到改善，但这需要额外获取场图数据。

离共振校正的数学方法

超越傅里叶的信号模型

离共振效应使得信号不再是简单的傅里叶变换关系，需要更复杂的模型来描述。考虑第 m 个采样时刻的信号：

其中 f_n 是第 n 个体素的图像值，\omega_n = \omega(r_n) 是该体素处的离共振频率。离共振相位项 e^{-i\omega_n t_m} 使得正向模型不再是纯粹的傅里叶矩阵乘法，而是一个更复杂的非均匀算子。

基函数展开近似

直接处理上述模型的计算复杂度很高，因为需要为每个体素-时间对计算离共振相位。Sutton等人（IEEE TMI 2003）提出了一种基于基函数展开的近似方法。核心思想是将离共振相位项用少量基函数的线性组合来近似：

其中 b_{m,l} 是时间相关的系数，c_{l,n} 是空间相关的系数，L 是基函数的数量（通常远小于 M 和 N）。将这个近似代入信号方程：

这个形式的优势在于：内层求和是标准的非均匀傅里叶变换（NUFFT），可以用快速算法计算；外层只是 L 个这样的变换的加权求和。

基函数的联合估计

基函数系数 B = (b_{m,l}) 和 C = (c_{l,n}) 可以通过矩阵分解从场图数据中估计。定义离共振矩阵 E，其元素为：

这是一个 M \times N 的矩阵。基函数估计问题转化为低秩矩阵近似：

这个问题可以通过奇异值分解（SVD）或其他低秩分解方法求解。

磁敏感加权成像中的离共振校正

三种校正策略的对比

在实际应用中，离共振校正面临一个权衡：需要场图信息，但获取场图需要额外的扫描时间。以加速因子 AF = 20 的磁敏感加权成像（SWI）为例，对比三种策略的效果。

无校正时，图像存在明显的信号丢失和几何失真（红色箭头标注的区域）。使用额外采集的 \Delta B_0 场图进行校正，图像质量显著改善（绿色和黄色箭头标注的区域恢复），但需要额外约15分钟的场图扫描时间。使用从数据本身估计的 \Delta B_0 场图进行校正，同样能够获得良好的校正效果，无需额外扫描，但需要约8小时的计算时间来迭代估计场图。

Daval-Frérot等人（MRM 2022）发展了从非笛卡尔SWI数据中迭代估计静态场图的方法，实现了无需额外采集的离共振校正。

轨迹设计中的离共振考量

信号的时间依赖性

回顾SPARKLING的优化目标，吸引项和排斥项都只关注采样点在k空间中的空间分布，而没有考虑采样的时间顺序。然而，离共振效应与采样时刻 t 直接相关：相位累积 \omega(r)t 随时间增加。

更完整的信号模型需要考虑多种时间依赖效应。定义复合衰减参数：

其中 \alpha(r) 代表 T_2^* 弛豫效应（信号随时间指数衰减），\omega(r) 代表 B_0 不均匀性效应（相位随时间累积）。完整的信号方程为：

从轨迹的时间着色图可以看到，标准SPARKLING优化只保证空间分布均匀，但同一时刻采样的点（相同颜色）可能分布不均匀。轨迹开始时采样的点（蓝色）和结束时采样的点（黄色）在空间上的分布没有特殊约束。这意味着某些空间区域可能只在较晚时刻被采样，此时该区域的信号已经因 T_2^* 衰减和离共振相位累积而显著变化。

MORE-SPARKLING：引入时间权重

时间加权排斥项

为了在轨迹设计中考虑时间因素，Chaithya等人（MRM 2023）提出了MORE-SPARKLING（Minimized Off-Resonance Effects SPARKLING）方法。核心改进是在排斥项中引入时间权重。

标准SPARKLING的排斥项只考虑空间距离：

MORE-SPARKLING修改排斥项为包含时间权重的形式。对于2D情况：

对于3D情况，引入显式的时间权重因子：

其中 \tau \geq 0 是时间权重参数，t_i 和 t_j 分别是采样点 i 和 j 的采集时刻，p 是归一化因子。

时间权重的作用机制

指数权重因子 e^{\frac{|t_i - t_j| \tau}{p}} 的行为如下：当 t_i = t_j（两个采样点在同一时刻采集）时，权重因子等于1，排斥强度与标准SPARKLING相同；当 t_i \neq t_j（两个采样点在不同时刻采集）时，权重因子大于1，排斥强度增加。时间差 |t_i - t_j| 越大，权重因子越大，排斥越强。

这种设计的物理意义是：增加时间上相距较远的采样点之间的排斥力。这会促使优化算法生成这样的轨迹：在任意给定时刻，当时采样的点在k空间中分布更均匀。换句话说，轨迹不仅在整体上覆盖目标密度，而且在每个时间窗口内也有较好的k空间覆盖。

从优化后的轨迹时间着色图可以观察到，MORE-SPARKLING生成的轨迹在每个时间切片内的采样点分布更加均匀，相邻时刻的采样点在空间上有更好的交织，这有助于减轻离共振效应和 T_2^* 衰减对重建图像的不利影响。

MORE-SPARKLING的实验验证

活体成像结果对比

在活体脑部成像中验证了MORE-SPARKLING的效果。场图显示离共振频率分布范围约为 -400 Hz 到 +400 Hz，在额叶和颞叶区域存在显著的场不均匀性。

对比三种重建策略：使用标准SPARKLING需要额外的内部场图估计，重建时间约8小时；使用MORE-SPARKLING（\tau = 1.0 或 \tau = 1.2），无需额外的场图估计，重建时间仅约15分钟。从轴位和矢状位图像可以看到，MORE-SPARKLING在保持图像质量的同时大幅缩短了重建时间，绿色箭头标注的区域显示两种方法都能有效校正离共振伪影。

GoLF SPARKLING：低频网格化采样

变密度采样的根本问题

回顾变密度采样的原理：小波分解中的近似系数（低频成分）包含图像的主要结构信息，这部分无法通过稀疏约束恢复，必须在数据中直接采样获取；细节系数（高频成分）是稀疏的，可以通过稀疏约束从欠采样数据中恢复。因此，变密度采样在k空间中心（对应低频）的采样密度必须足够高。

标准SPARKLING虽然整体遵循变密度分布，但k空间中心的采样仍然是非笛卡尔的、不规则的。在某些情况下，中心区域的采样可能超过奈奎斯特率（过采样），但由于点分布不均匀，信息提取效率不是最优的。

低频笛卡尔网格化

GoLF SPARKLING（Gridding of Low Frequencies）的核心思想是：在k空间中心区域强制使用笛卡尔网格采样，这是对该区域最优的采样方式；在外围区域保持非笛卡尔SPARKLING轨迹，以获得更好的k空间覆盖效率。

这种混合设计带来额外的好处：中心的笛卡尔采样数据可以通过简单的逆傅里叶变换（IFFT）快速生成低分辨率图像，这个低分辨率图像可以用于估计线圈敏感度图，无需额外的校准扫描。

GoLF约束的数学表述

标准SPARKLING的约束集为：

GoLF SPARKLING在此基础上添加了低频网格化约束：

最后一个约束 k_i[n_i + l] = x_i^l 规定了每条轨迹在经过k空间中心区域时必须通过预定义的笛卡尔网格点 x_i^l。

从示意图可以看到，每条轨迹曲线在中心区域（绿色高亮部分）被约束到规则的网格位置上，而在外围区域则自由优化。

GoLF与MORE-SPARKLING的联合应用

多种方法的对比实验

在600 µm各向同性分辨率、加速因子 AF = 15、扫描时间3分20秒的条件下，对比了多种轨迹和重建方法。

笛卡尔参考图像（GRAPPA p4）作为质量基准。场图显示离共振频率分布范围约 -200 Hz 到 +200 Hz。标准SPARKLING（\tau = 0）的重建图像在某些区域存在伪影（红色箭头标注）。MORE-SPARKLING（\tau = 1.0）通过时间权重优化改善了部分伪影（绿色箭头标注），但仍有残余。GoLF-SPARKLING（\nu = 1）通过低频网格化也提供了改善，但效果与MORE-SPARKLING不同。GoLF + MORE-SPARKLING（\tau = 1.0, \nu = 1）结合了两种优化策略，获得了最佳的图像质量（橙色箭头标注区域显示最少伪影）。

作为对比，传统的Twisted Projection Imaging（TPI）和FLORET轨迹在相同条件下的表现不如优化的SPARKLING变体（红色和橙色箭头标注的区域显示更多伪影）。

2D+和3D非笛卡尔轨迹的系统分类

轨迹类型的维度

非笛卡尔k空间轨迹可以按照多个维度进行分类。基本2D轨迹类型包括：Radial（径向）——从中心向外的直线；Spiral（螺旋）——从中心向外盘旋的曲线；Cones（锥形）——具有锥面结构的轨迹；SPARKLING——优化生成的变密度曲线。

从2D扩展到3D的策略包括：Stacked（堆叠）——将2D轨迹在第三维度上平行堆叠；Rotated（旋转）——将2D轨迹绕某轴旋转形成3D覆盖；Conical（锥形）——将2D轨迹沿锥面排列。

完全3D的轨迹设计包括：3D Cones——三维锥形轨迹；TPI（Twisted Projection Imaging）——扭曲投影成像轨迹；3D SPARKLING——完全三维优化的SPARKLING；3D MORE-SPARKLING——带时间权重的三维SPARKLING。

轨迹性能的基准测试

视觉质量对比

对多种2D+和3D轨迹进行了系统的基准测试。

2D+轨迹（通过堆叠、旋转或锥形扩展到3D）的对比显示：Spiral轨迹在某些区域存在环形伪影（红框标注）；SPARKLING轨迹的伪影模式不同，且整体水平较低。3D轨迹的对比显示：3D Cones和TPI存在特定的伪影模式；3D SPARKLING的图像质量较好；3D MORE-SPARKLING在离共振严重区域（对照场图）表现最佳。

参考图像和 \Delta B_0 场图作为评估基准，场图显示的离共振范围约 -200 Hz 到 +200 Hz，主要集中在额叶底部区域。

定量性能指标

使用SSIM（结构相似性指数）和PSNR（峰值信噪比）作为定量评估指标。

2D+轨迹的结果显示，对于不同的扩展策略（Stacked、Rotated、Conical），SPARKLING在SSIM和PSNR上通常优于或接近其他基本轨迹类型。例如，Rotated Spiral达到SSIM = 0.886和PSNR = 27.82，而Rotated SPARKLING达到SSIM = 0.850和PSNR = 24.96。

3D轨迹的结果更加突出：3D Cones的SSIM = 0.827，PSNR = 23.10；TPI的SSIM = 0.870，PSNR = 25.46；3D SPARKLING的SSIM = 0.890，PSNR = 27.71；3D MORE-SPARKLING达到最高的SSIM = 0.907和PSNR = 28.93。

从箱线图可以看到，3D MORE-SPARKLING不仅在平均性能上最优，而且方差较小，表明其在不同数据集上的表现更加稳定。这验证了基于优化的轨迹设计（特别是考虑时间因素的MORE-SPARKLING）相比传统几何轨迹的显著优势。

加速解剖MRI成像

本节将前面介绍的各种技术整合应用于解剖MRI成像，展示如何将GRAPPA并行成像与GoLF SPARKLING轨迹相结合，实现高质量的快速解剖扫描。这种组合策略充分利用了笛卡尔采样在并行成像中的优势和非笛卡尔采样在k空间覆盖效率上的优势。

从k空间采样模式的示意图可以看到三个区域的划分：ACS（Auto-Calibration Signal）区域位于中心，采用全采样的笛卡尔网格用于GRAPPA校准和敏感度估计；GA（GRAPPA Accelerated）区域是中间环形区域，采用欠采样的笛卡尔网格，通过GRAPPA填充；NC（Non-Cartesian）区域是外围区域，采用SPARKLING非笛卡尔轨迹覆盖。

MPRAGE序列用于T1加权成像

MPRAGE的基本原理

MPRAGE（Magnetization Prepared Rapid Gradient Echo）是获取高分辨率T1加权解剖图像的标准序列。T1加权成像能够提供良好的灰质（Grey Matter, GM）、白质（White Matter, WM）和脑脊液（Cerebrospinal Fluid, CSF）之间的对比度，是脑结构分析、组织分割和形态学测量的基础。

MPRAGE序列的核心是反转恢复（Inversion Recovery）预备脉冲。首先施加一个非选择性的180°反转脉冲（Nonsel 180°），将所有组织的纵向磁化矢量 M_z 反转到负值。然后在反转时间（TI）期间，不同组织根据各自的T1弛豫时间以不同速率恢复。

从反转恢复曲线可以看到，白质（红色曲线）恢复最快，灰质（绿色曲线）次之，脑脊液（蓝色曲线）恢复最慢。

在恢复过程的特定时刻开始数据采集，此时不同组织的 M_z 值差异最大，从而产生最佳的T1对比度。数据采集使用一系列小翻转角 \alpha° 的梯度回波脉冲，快速采集k空间数据。时序图下方的彩色条带表示采集窗口内不同时刻对应的k空间位置——这与信号的T1恢复状态相关，需要在轨迹设计中考虑。

传统笛卡尔采样与SPARKLING

传统MPRAGE使用笛卡尔逐行采样，k空间数据呈现规则的网格模式。笛卡尔采样的优势是重建简单（直接FFT），且与GRAPPA等并行成像技术兼容性好。但笛卡尔采样的k空间覆盖效率有限，特别是在高分辨率3D成像中，扫描时间较长。

将SPARKLING轨迹应用于MPRAGE序列可以提高k空间覆盖效率。SPARKLING的变密度特性与MPRAGE的需求匹配：k空间中心的高密度采样确保了图像对比度的准确性，外围的欠采样通过压缩感知重建恢复细节信息。

3D SPARKLING轨迹的示意图显示了变密度的球形分布，中心区域采样点密集（暖色），外围区域稀疏（冷色）。

GoLF-SPARKLING用于MPRAGE序列

GoLF-SPARKLING特别适合MPRAGE应用，因为它在k空间中心强制使用笛卡尔网格采样。这带来多重好处：中心的笛卡尔数据可以直接用于GRAPPA校准；低分辨率图像可以快速生成用于敏感度估计；与传统并行成像流程兼容。

从GoLF-SPARKLING的3D轨迹图可以看到，轨迹在中心区域被约束到规则网格（绿色高亮的直线段），在外围区域则呈现优化后的曲线形态。这种混合设计在保持SPARKLING高效覆盖优势的同时，获得了笛卡尔采样在重建和校准方面的便利性。

将笛卡尔加速引入非笛卡尔MRI

混合采样策略

GoLF-SPARKLING的采样模式可以分解为三个组成部分。完整轨迹 K 包含非笛卡尔部分 K_{nc}（外围的SPARKLING曲线）和笛卡尔部分 K_c（中心的网格采样）。笛卡尔部分进一步可以在相位编码方向上进行GRAPPA欠采样。

从采样密度图可以直观理解不同策略的效果。无GRAPPA时（GS No GRAPPA），笛卡尔中心区域是全采样的，采样密度为1（奈奎斯特率）。使用GRAPPA=2加速后（GS GRAPPA=2），笛卡尔中心区域的采样密度降为0.5，但通过GRAPPA重建可以恢复完整数据。进一步增加外围非笛卡尔区域的欠采样因子（AF=5, AF=8），可以实现更高的整体加速。

两步重建流程

混合采样的重建分为两个步骤。第一步是GRAPPA重建：对中心笛卡尔区域的欠采样数据应用GRAPPA算法，填充缺失的k空间行，得到完整的中心k空间数据。第二步是压缩感知重建：将GRAPPA重建后的中心数据视为已采集数据，与外围非笛卡尔数据一起进行CS重建。

这种两步策略的优势在于：GRAPPA在中心笛卡尔区域工作最可靠（ACS数据充足、g因子低）；CS在外围稀疏区域更有效（高频细节本身是稀疏的）。两种方法各自发挥所长，实现协同加速。

GRAPPA在SPARKLING中的重要性

不同加速策略的对比

实验比较了多种加速策略在MPRAGE成像中的表现。笛卡尔参考图像（扫描时间9.2分钟）作为质量基准。GoLF-SPARKLING无GRAPPA加速（AF=5，扫描时间1.8分钟）的重建图像存在一定程度的噪声和模糊（红框标注区域）。GoLF-SPARKLING结合GRAPPA=2加速（AF=5，扫描时间1.8分钟），在相同扫描时间下图像质量明显改善（绿框标注区域更清晰）。GoLF-SPARKLING结合GRAPPA=2×2加速（AF=8，扫描时间1.1分钟），进一步缩短扫描时间的同时保持了可接受的图像质量（黄框标注区域）。

放大的细节对比更清楚地显示了GRAPPA加速的作用：无GRAPPA时的图像（b）存在明显的噪声纹理；加入GRAPPA=2后（c）纹理减少，细节更清晰；GRAPPA=2×2（d）虽然加速更高但质量仍可接受。

非笛卡尔加速的必要性探讨

笛卡尔与非笛卡尔的系统对比

一个自然的问题是：既然GRAPPA在笛卡尔采样上工作良好，是否还需要非笛卡尔加速？为回答这个问题，进行了系统的对比实验。

纯笛卡尔方案的结果：笛卡尔参考需要9分11秒；GRAPPA=2×2（AF=4）需要2分17秒；GRAPPA=4×2（AF=8）需要1分6秒；GRAPPA=3×3（AF=9）需要1分钟。随着加速因子增大，笛卡尔GRAPPA图像的噪声逐渐增加，AF=9时图像质量明显下降。

GoLF-SPARKLING方案的结果：GRAPPA=2（AF=5）需要1分48秒；GRAPPA=2×2（AF=8）需要1分6秒；GRAPPA=2×2（AF=10）仅需54秒。在相同或更高的加速因子下，GoLF-SPARKLING的图像质量优于纯笛卡尔GRAPPA。

非笛卡尔加速的优势

对比结果表明，在高加速因子下，非笛卡尔采样确实提供了额外的优势。这种优势来源于：SPARKLING的变密度采样更好地匹配了k空间能量分布；非笛卡尔轨迹的欠采样伪影更分散，更适合CS重建；GoLF设计使得GRAPPA在中心区域工作在最优条件下。

在AF=10、扫描时间54秒的条件下，GoLF-SPARKLING仍然保持了可用的图像质量，而纯笛卡尔方案在AF=9时已经出现明显的质量下降。这验证了将并行成像与非笛卡尔压缩感知相结合的策略在高加速应用中的价值。

加速功能磁共振成像采集

这一部分的核心目标是：在给定的扫描时间预算内，提升fMRI的空间分辨率。传统fMRI面临一个根本性的权衡问题——如果想要更高的空间分辨率，就需要更长的扫描时间；如果想要缩短扫描时间，就必须牺牲分辨率。

高场强fMRI的科学需求

BrainSync项目的第二工作包（WP2）致力于在11.7T超高场强下量化和建模学习与决策的神经机制。要实现这一科学目标，必须拥有一套性能卓越的fMRI序列。

从神经成像技术的覆盖范围来看，不同技术在空间和时间尺度上各有优势。在空间维度上，大脑结构从宏观到微观依次为脑回（Gyri）、皮层柱（Columns）、皮层层次（Layers）、神经元（Neurons）、树突（Dendrites）和突触（Synapses）。在时间维度上，神经活动发生在毫秒到小时的不同尺度。EEG和MEG具有极高的时间分辨率（毫秒级），但空间分辨率有限；光学成像和电生理技术可以达到细胞甚至突触级别的空间分辨率，但覆盖范围受限；fMRI则在空间覆盖和分辨率之间取得了较好的平衡，能够实现全脑成像。

要成功开展高场强fMRI研究，需要同时满足三个技术指标：亚毫米级的空间分辨率，以分辨皮层柱和皮层层次结构；小于0.5秒的时间分辨率，以捕捉血流动力学响应的动态变化；以及全脑覆盖能力，以研究分布式的神经网络。实现这些目标需要两个关键要素：一台性能卓越的MRI设备，以及一种创新的采集加速方法。

SPARKLING加速采集方案

即使拥有超高场强MRI提供的SNR优势，要在合理的扫描时间内实现高分辨率全脑fMRI仍然需要采集加速技术。本课程前面介绍的SPARKLING轨迹正是为此目的设计的解决方案。

针对fMRI应用，SPARKLING需要满足以下具体指标：空间分辨率达到1-1.5毫米（各向同性），以分辨皮层内的精细结构；时间分辨率优于0.5秒，以adequate采样血流动力学响应；同时实现全脑覆盖，以研究大尺度脑网络。SPARKLING通过其变密度非笛卡尔采样策略，结合压缩感知重建，能够在满足这些要求的同时实现显著的采集加速。

高分辨率fMRI的技术框架

fMRI数据处理的完整流程

高分辨率fMRI从数据采集到神经科学发现需要经历四个主要阶段。第一阶段是采集（Acquisition），MRI扫描仪按照设计的k空间轨迹采集原始数据，这些数据随时间累积形成一系列k空间体积。第二阶段是重建（Reconstruction），将欠采样的k空间数据通过压缩感知或并行成像算法转换为图像序列。第三阶段是统计分析（Statistical Analysis），对重建的图像时间序列进行统计建模，生成z-score激活图来识别与任务相关的脑区。第四阶段是神经科学解读（Neuroscience），将统计结果与解剖结构对应，得出关于脑功能的科学结论。

BOLD信号本质上是一个四维数据结构，包含三个空间维度和一个时间维度。每个时间点对应一个完整的三维脑部体积，整个fMRI实验产生数百个这样的体积。这种4D特性为采集加速提供了额外的自由度——不仅可以在空间维度上进行欠采样，还可以在时间维度上设计采样策略。

变密度采样与k空间覆盖

fMRI采集加速的核心策略是使用4D欠采样模式来高效获取数据。图像域与k空间通过傅里叶变换相互关联，MRI图像的低频信息（整体对比度和结构）集中在k空间中心，高频信息（边缘细节）分布在k空间外围。目标采样密度的设计遵循这一能量分布特性：中心区域采样密度高（图中亮区），外围区域采样密度低（图中暗区）。这种变密度设计确保了关键的低频信息被充分采样，而稀疏的高频信息可以通过压缩感知重建恢复。

BOLD信号的动态采样策略

时变欠采样的基本概念

传统方法使用静态采样模式，即每个时间帧采用完全相同的k空间轨迹。Comby等人2025年发表在Imaging Neuroscience上的工作提出了动态采样策略，利用SNAKE仿真工具生成模拟fMRI数据来验证这一方法。

时变欠采样使用螺旋堆叠（stack of spiral）采样方案，沿 k_z 维度进行欠采样，并且在不同时间帧之间改变采样位置。从全采样到不同加速因子 R=2、R=4、R=8 的对比可以看到，随着加速因子增大，每个时间帧内采集的 k_z 层数减少，但不同帧采集的层位置会有所不同。这种设计的关键在于：虽然单个时间帧的k空间是不完整的，但多个相邻时间帧的数据组合起来可以覆盖更完整的k空间。

静态采样与动态采样的对比

传统fMRI采用静态采样方式，每个时间帧独立采集完整（或固定欠采样模式）的k空间数据。

图中显示，在静态采样下，k空间位置随时间保持不变（绿色方块位置固定），完成一个完整体积需要2.4秒的时间分辨率。k空间中心位置 K=0 的采集时刻在每帧中相同。

4D动态fMRI采用时变采样策略，k空间采样位置随时间帧变化。从图中可以看到，不同时间点的k空间采集位置（绿色方块）在垂直方向上有所偏移，K=0 的采集时刻也随之变化。这种设计带来两个优势：第一是提高时间分辨率，由于每帧采集的数据量减少，单帧采集时间缩短，可以实现1.2秒甚至更短的时间分辨率；第二是支持全局重建，重建算法可以同时利用多个时间帧的k空间数据，将它们视为对同一动态对象的联合观测，从而更高效地利用所有采集数据。

动态采样的物理基础在于BOLD信号的时间连续性——相邻时间帧的脑部图像高度相似，信号变化相对缓慢。这种时间冗余性允许重建算法借用相邻帧的信息来填补当前帧k空间的缺失部分，实现比单帧独立重建更好的图像质量。

动态采样的实验验证

不同加速因子的性能对比

对动态采样策略在不同加速因子下的表现进行了系统评估。从k空间采样模式来看，全采样情况下所有 k_z 层都被采集；R=2 时采集一半的层；R=4 时采集四分之一的层；R=8 时仅采集八分之一的层。对应的z-score激活图显示了重建后进行统计分析得到的脑激活模式，z值范围从-11到+11，正值表示任务相关的激活，负值表示去激活。

从激活图可以观察到，R=4 加速因子（橙色框标注）表现出比全采样更强的激活信号。这个看似反直觉的结果可以解释为：动态采样结合全局重建能够更好地利用时间维度的信息冗余，在适当的加速因子下反而能够提升统计检测的敏感性。R=8 时加速过度，k空间采样过于稀疏，重建质量下降导致激活检测性能降低。

F2-score评估指标

使用F2-score作为激活检测性能的定量评估指标。F2-score是精确率和召回率的加权调和平均，对召回率给予更高权重，适合评估激活检测这类更关注避免漏检的任务。图中比较了两种重建方法：COLD重建（蓝色曲线）和MERGE GLOBAL + COLD重建（黄色曲线）。MERGE GLOBAL方法利用多个时间帧的k空间数据进行联合重建，而COLD是一种基于低秩约束的动态MRI重建算法。

曲线显示，结合全局重建的方法（MERGE GLOBAL + COLD）在各加速因子下都优于单独的COLD重建。两种方法都在 R=4 附近达到峰值性能，F2-score约为0.82-0.85。这个结果不仅实现了更高的时间分辨率，还验证了动态采样策略的有效性。目前正在使用MRI扫描仪采集的前瞻性活体数据对这些仿真结果进行验证。

高时间分辨率功能MRI的完整方案

EPI-SPARKLING混合轨迹

针对高时间分辨率fMRI应用，提出了EPI-SPARKLING混合采集方案。这种方案将传统EPI（Echo Planar Imaging）的快速采集特性与SPARKLING的优化k空间覆盖相结合。

无加速的EPI-SPARKLING轨迹在 (k_x, k_y) 平面内展示了完整的采样模式，每个 k_z 位置都采用优化的非笛卡尔轨迹。R=4 加速的EPI-SPARKLING则通过减少采集的 k_z 层数实现加速，同时在平面内保持高效的SPARKLING覆盖。这种设计由Pan等人开发，相关工作正在准备发表中。

深度学习重建

在 R=4 加速条件下，采用深度学习方法进行图像重建。重建得到的z-score激活图清晰显示了任务相关的脑区激活模式，z值分布从-11到+11。Comby等人在JOSS 2025和Imaging Neuroscience 2025上发表的工作为这一流程提供了软件工具和方法学基础。

重建的BOLD信号时间曲线展示了感兴趣区域（ROI）内信号随时间的变化。图中比较了原始BOLD信号（蓝色实线，经过中心化处理）和Savitzky-Golay平滑后的信号（粉色虚线）。原始信号均值约为0.0342，信号波动范围在±0.0010左右，体现了典型的BOLD响应幅度（约百分之几的信号变化）。

阶段性总结

技术发展现状

fMRI加速技术传统上依赖扫描-重复策略，即每个时间帧独立采集和重建。压缩感知理论为动态成像提供了新的框架，其核心是时变或动态采样——不同时间帧使用不同的k空间采样位置，使得联合重建能够更高效地利用数据。

这种动态采样策略在堆叠式2D采样模式中容易实现，包括螺旋堆叠和SPARKLING堆叠。通过改变不同时间帧在堆叠维度（k_z）上的采样位置，即可实现时变欠采样。EPI仍然是fMRI领域的主力序列，正如它在扩散加权成像（dwMRI）中的地位一样，其快速采集特性是实现高时间分辨率的基础。

对于解剖成像，GoLF SPARKLING与GRAPPA的结合已经展示了在一分钟内完成高质量解剖扫描的能力，为快速定位和配准提供了支持。

未来研究方向

图像重建方面的发展方向是利用学习先验来改进重建质量。NC-PDNet（非笛卡尔原始对偶网络）将深度学习与优化算法相结合；PnP（Plug-and-Play，即插即用）方法使用预训练的去噪网络作为正则化器；扩散模型（Diffusion models）则利用生成模型的强大表示能力来指导重建。这些方法都旨在超越传统的手工设计稀疏先验，学习更符合真实图像统计特性的先验知识。

临床转化方面，计划在Henri-Mondor医院进行部署，由Blanche Bapst医生进行放射学评估。这种从方法学研究到临床应用的转化是验证技术实用价值的关键步骤。

本讲内容总结

2D SPARKLING

2D SPARKLING轨迹特别适合T2*加权成像，这类成像使用长回波时间（TE）、长重复时间（TR）的梯度回波序列。在与其他非笛卡尔采样方案的对比中，2D SPARKLING表现出更好的图像质量。2D SPARKLING的一个局限是在z方向（层面选择方向）的分辨率有限，只能实现平面内的高分辨率。

3D SPARKLING

3D SPARKLING的实现存在不同的策略选择。堆叠策略（将2D SPARKLING分层堆叠）的表现优于局部3D优化方法。全局优化的完全3D轨迹又优于堆叠策略，能够获得最佳的点扩散函数特性。完全3D优化的计算需求较高，但通过快速多极方法和GPU加速，可以在可接受的时间内完成优化，从而实现真正的高分辨率三维成像。

MORE-SPARKLING

MORE-SPARKLING通过在排斥项中引入时间权重，对轨迹施加时间平滑性约束。这种设计的目的是抵消数据采集过程中由于 B_0 场不均匀性导致的离共振效应。时间权重使得在任意时刻采样的点在k空间中分布更均匀，减少了对额外场图校准扫描或耗时的内部场图估计的依赖。

GoLF SPARKLING与GRAPPA结合

GoLF SPARKLING通过在k空间中心强制使用笛卡尔网格采样，使采样密度在低频区域达到奈奎斯特率。这种设计对于非笛卡尔解剖MRI成像是必要的，因为它确保了良好的图像对比度。GoLF SPARKLING与GRAPPA并行成像的结合可以在一分钟内完成解剖扫描，实现了临床可行的超快速高分辨率成像。

阶段性总结与展望

技术发展脉络

MR加速技术经历了从并行成像到压缩感知再到优化轨迹设计的发展过程。传统的MR加速主要依赖GRAPPA等并行成像方法，利用多线圈的空间编码冗余实现加速。压缩感知理论揭示了变密度采样（VDS）是实现高倍加速的关键，稀疏性先验使得从欠采样数据中恢复图像成为可能。非笛卡尔成像范式通过SPARKLING等优化轨迹实现了变密度采样的高效实现，突破了笛卡尔采样的固有限制。

关键技术要点

GoLF-SPARKLING实现了笛卡尔采样与非笛卡尔成像的最优融合：中心区域采用笛卡尔网格达到奈奎斯特率，确保低频信息完整；外围区域采用优化的非笛卡尔轨迹实现高效加速。这种混合策略在保持图像对比度的同时最大化了加速潜力。

对于非笛卡尔解剖MRI成像，GoLF-SPARKLING的低频网格化设计是获得良好图像对比度的内在需求。GoLF SPARKLING与GRAPPA的结合已经展示了在一分钟内完成高质量解剖扫描的能力。

未来发展方向

下一步的研究方向包括：利用学习先验改进图像重建，如NC-PDNet（非笛卡尔原始对偶网络）、PnP（即插即用）方法和扩散模型等深度学习技术；以及在Henri-Mondor医院进行临床部署，由Blanche Bapst医生进行放射学评估，验证这些技术在实际临床环境中的可行性和诊断价值。