基于模型的高级图像重建方法

本讲属于深度学习时代计算MRI课程的第六讲，主要内容是从基于模型的方法过渡到基于学习的图像重建方案。本节首先介绍稀疏正则化器及其相关优化算法。

稀疏性模型的两种形式

在压缩感知MRI重建中，稀疏性假设是核心。实现稀疏正则化有两种不同的建模方式：综合模型和分析模型。

综合模型

综合模型假设待重建的图像 x 可以表示为某个字典或基矩阵与稀疏系数的乘积，即：

其中 B 是一个 N \times K 的矩阵，通常称为字典或基矩阵。这里 N 是图像的维度（像素数），K 是字典的原子数。当 K > N 时，B 是一个宽矩阵，称为过完备字典，意味着用于表示信号的基向量数量超过了信号本身的维度，这样可以提供更灵活的稀疏表示。系数向量 z \in \mathbb{C}^K 被假设是稀疏的，即大部分元素为零或接近零。为了促进稀疏性，使用 \ell_1 范数 \|z\|_1 作为正则项。

分析模型

分析模型则从另一个角度出发，假设对图像 x 施加某个变换 T 后得到的结果是稀疏的，即 Tx 是稀疏的。这里 T 是一个 K \times N 的变换矩阵，通常是高矩阵，即 K > N。典型的例子包括有限差分算子（对应全变分正则化）或小波变换。相应的正则项为 \|Tx\|_1。这种分析模型在FDA批准的临床压缩感知MRI方法中被广泛采用。

两种模型的关系与权衡

当 B = T^{-1} 时，综合模型和分析模型在数学上是等价的。然而在实际应用中，B 和 T 通常都不是方阵，因此无法简单地互相求逆。从优化角度看，综合模型中对 \|z\|_1 的优化通常比分析模型中对 \|Tx\|_1 的优化更容易处理，因为前者直接对稀疏系数施加约束。但代价是综合模型的优化变量维度为 K，当使用过完备字典时 K > N，变量维度更高。

综合模型的优化问题

LASSO问题的形式

基于综合模型的压缩感知MRI重建，其典型的代价函数形式为：

这里 A 是MRI的测量矩阵（包含傅里叶欠采样和线圈灵敏度等），y 是采集到的k空间数据，\beta > 0 是正则化参数，用于平衡数据保真项和稀疏正则项。第一项 \frac{1}{2}\|ABz - y\|_2^2 是数据保真项，衡量重建结果与观测数据的一致性；第二项 \beta\|z\|_1 是稀疏正则项，鼓励系数向量 \hat{z} 具有稀疏结构。

\ell_1 范数是 \ell_0 计数测度的凸松弛。\ell_0 范数直接计算非零元素的个数，但它是非凸的，导致优化问题是NP难的。\ell_1 范数作为最紧的凸松弛，在一定条件下可以恢复与 \ell_0 相同的稀疏解，同时保持问题的凸性，使得全局最优解可以被高效求得。这个优化问题在统计学中被称为LASSO问题。

近端梯度方法

求解上述LASSO问题需要专门的优化算法，因为 \ell_1 范数虽然是凸的，但在零点处不可微。近端方法是处理这类问题的主流框架。

ISTA算法

迭代软阈值算法（Iterative Soft Thresholding Algorithm, ISTA）是最经典的近端梯度方法，也称为近端梯度方法（Proximal Gradient Method, PGM）或近端前向-后向分裂法。其迭代格式为：

这个迭代可以理解为两步：首先沿着数据保真项的负梯度方向进行梯度下降，得到一个中间结果；然后对这个中间结果应用软阈值操作，实现对 \ell_1 正则项的近端映射。

软阈值函数的定义为：

这个操作将绝对值小于阈值 c 的分量置零，将绝对值大于 c 的分量向零收缩 c 的距离。这正是 \ell_1 范数的近端算子的解析形式。

式中 D = \text{diag}\{d\} 是一个对角矩阵，需要满足条件 D \succeq B'A'AB，即 D - B'A'AB 是半正定的。这个条件保证了算法的收敛性。标准ISTA对应特殊情况 D = LI，其中 L = \|B'A'AB\|_2 是矩阵 B'A'AB 的谱范数（最大特征值）。使用对角矩阵 D 而非标量倍数的单位矩阵是对标准ISTA的推广，可以根据不同坐标方向的曲率选择不同的步长。

ISTA的收敛速度为 O(1/k)，即目标函数值与最优值的差距以迭代次数的倒数衰减，这个速度相对较慢。

FISTA算法

快速迭代软阈值算法（Fast Iterative Soft Thresholding Algorithm, FISTA）通过引入动量项加速收敛，也称为快速近端梯度方法（Fast Proximal Gradient Method, FPGM）。FISTA在每次迭代中不仅使用当前点的信息，还利用前一次迭代点的信息构造一个外推点，然后在外推点处进行梯度下降和近端操作。这种加速技术源于Nesterov的加速梯度方法，将收敛速度从 O(1/k) 提升到 O(1/k^2)，即平方加速。

近端优化梯度方法

POGM算法框架

考虑更一般的复合优化问题：

其中 f(x) 是光滑函数（具有Lipschitz连续梯度），g(x) 是近端友好函数，即其近端算子有解析解或可以高效计算。在MRI重建问题中，f(x) 对应数据保真项，g(x) 对应稀疏正则项。

近端优化梯度方法（Proximal Optimized Gradient Method, POGM）是近年提出的改进算法。与FISTA相比，POGM在最坏情况下的收敛界改善了约2倍。从实现复杂度来看，ISTA、FISTA和POGM三者相近，每次迭代的计算量相当，主要由梯度计算 \nabla f 和近端算子 \text{prox}_g 的求解决定。

在实际应用中，POGM的收敛速度经验上优于FISTA，特别是当结合自适应重启策略时效果更为显著。自适应重启是指当检测到目标函数值不再下降或出现振荡时，重置动量项，重新开始加速过程，这可以避免加速方法在接近最优解时的振荡现象。这一领域仍是活跃的研究方向，不断有新的算法改进被提出。

FISTA与POGM的迭代格式

算法初始化与参数更新

FISTA和POGM算法共享相同的初始化设置：令 w_0 = x_0，\theta_0 = 1。随后对 k = 1, 2, \ldots, N 进行迭代。

动量参数 \theta_k 的更新规则根据迭代阶段有所不同：

当 k < N 时使用第一个公式，这对应标准的FISTA动量参数更新；当 k = N（最后一次迭代）时使用第二个公式，这是POGM的特殊处理，通过在最后一步使用更大的动量来进一步优化收敛性能。步长参数 \gamma_k 定义为：

其中 L 是目标函数光滑部分梯度的Lipschitz常数。

迭代主体步骤

每次迭代包含三个计算步骤。首先计算梯度下降点：

这一步沿着光滑项 f 的负梯度方向移动，步长为 1/L。

然后计算动量组合点：

这个表达式是POGM区别于FISTA的核心。它结合了三个历史信息：当前梯度步与上一梯度步的差异 (w_k - w_{k-1})、当前梯度步与上一迭代点的差异 (w_k - x_{k-1})、以及上一动量点与迭代点的差异 (z_{k-1} - x_{k-1})。这种多重动量组合使得POGM能够更有效地利用历史信息加速收敛。

最后执行近端映射：

近端算子的作用是在动量点 z_k 附近找到一个点，使其既接近 z_k，又能使非光滑正则项 g(x) 较小。参数 \gamma_k 控制正则化的强度。

POGM方法适用于最小化 f(x) + g(x) 形式的问题，要求 f 是凸函数且具有 L-Lipschitz连续梯度，g 是凸函数。结合自适应重启策略可以进一步提升实际性能。

算法性能对比

演示实验使用正交离散小波变换（ODWT）作为稀疏基 B，对欠采样k空间数据进行重建。k空间采样模式显示中心区域密集采样、外围稀疏采样的变密度策略。重建结果显示POGM配合ODWT能够达到1.71%的归一化均方根误差（NRMSE）。

从收敛曲线可以观察到：代价函数值方面，POGM（绿色）下降最快，FISTA（蓝色）次之，ISTA（红色）最慢；NRMSE曲线同样反映出POGM的优势，在相同迭代次数下达到更低的重建误差。三种算法最终收敛到相同的解，但POGM所需迭代次数明显更少。

低秩加稀疏重建模型

动态MRI的结构假设

动态MRI（如fMRI）数据具有特殊的时空结构，可以利用两个互补的先验假设：沿时间轴的稀疏性和空间上的低秩性。时间稀疏性意味着在某个变换域（如时间差分或傅里叶变换）中，信号随时间的变化是稀疏的，这符合压缩感知的框架。空间低秩性则反映了动态图像序列中各帧之间存在强相关性，将所有帧排列成矩阵后，该矩阵的秩远小于其维度。

L+S分解优化问题

基于上述假设，将动态图像 \hat{X} 分解为低秩分量 L 和稀疏分量 S 的和，优化问题形式为：

第一项是数据一致性项，A 是包含傅里叶变换和线圈灵敏度的测量算子，y 是采集的k空间数据。第二项 \|L\|_* 是核范数正则项，核范数定义为矩阵奇异值的 \ell_1 范数，它是矩阵秩的凸松弛，用于促进低秩结构。第三项 \|\mathcal{T}S\|_1 是稀疏正则项，\mathcal{T} 是某种稀疏变换（如时间有限差分），用于促进稀疏分量在变换域中的稀疏性。

正则化参数 \lambda_L 和 \lambda_S 都是正数，需要通过经验调节来平衡三项之间的权重。\lambda_L 过大会过度压制低秩分量的能量，\lambda_S 过大则会丢失稀疏分量中的细节信息。

SNAKE-fMRI模拟器

SNAKE是Simulator from Neuro-Activation to K-space Evaluation的缩写，是一个用于fMRI数据模拟的软件工具。它可以从神经激活模式出发，模拟完整的fMRI信号采集过程，生成k空间数据。

模拟器的输入包括大脑解剖图像和激活区域的定义。输出包括期望的BOLD响应信号和实际的激活事件时间序列。从演示图中可以看到，BOLD信号强度（蓝色曲线）呈现周期性变化，与激活事件（橙色阶梯曲线）的开关模式相对应。这种block设计范式（20秒开-20秒关）是fMRI实验中常用的刺激模式。

fMRI采集轨迹设计

二维与三维采集模式

fMRI数据采集可以采用2D或3D的k空间轨迹。

2D EPI SPARKLING模式在 k_x-k_y 平面内设计优化的非笛卡尔轨迹，但由于只在二维平面内运动，设计自由度受限。为了获得完整的三维k空间覆盖，需要在不同的 k_z 位置重复采集多个2D切片。

3D采集模式则在整个三维k空间中设计轨迹。演示图显示了加速倍数 R = 3.75 的3D螺旋轨迹，在中心位置采集7个切片。3D模式的优势在于沿 k_z 轴有更大的设计自由度，可以通过调整不同帧在z方向的采样位置来优化整体覆盖。右侧散点图显示了前100帧中采集螺旋的z值分布，体现了这种沿z轴的变化策略。

EPI SPARKLING采集参数

采集实验在7 Tesla超高场强下进行。2D Slices模式和3D Stacked模式的主要参数对比如下：

空间分辨率均为3mm各向同性。视场（FOV）方面，2D模式为 (182, 217, 3) mm，3D模式为 (182, 217, 182) mm，后者覆盖更大的z方向范围。

时间参数方面，2D模式的TR为2000ms，体积TR为2s，采集180帧共360秒；3D模式的TR仅为50ms，体积TR为0.8s，采集450帧同样持续360秒。3D模式通过更短的TR实现了更高的时间分辨率。

回波时间TE均为25ms，翻转角2D模式为73°、3D模式为10°。实验范式采用20秒开-20秒关的block设计。线圈数量2D模式使用8通道，3D模式使用单通道。每次观测时间均为30ms。2D模式采集单层切片，3D模式同时采集16层切片。信噪比设置为10000，视为无限大（理想情况）。

低秩加稀疏重建的FISTA算法

算法输入与参数说明

该算法选用FISTA是因为其快速收敛性和高精度。算法需要以下输入：y 是多线圈欠采样k空间数据；A 是信号采集的编码矩阵，包含傅里叶变换和线圈灵敏度信息；\lambda_L 是低秩正则化参数；\lambda_S 是稀疏正则化参数；\gamma 是梯度步长。算法输出 X 是包含所有重建帧的矩阵。

算法的设计基于不相干性假设：低秩分量 L 和稀疏分量 S 之间应当不相干，同时采集空间与这两个分量之间也应保持不相干。这种不相干性保证了分解的唯一性和稳定性。

初始化步骤

初始化设置为：\chi_0 = M_0 = A^*y，即用测量数据的伴随变换结果作为初始估计；\tilde{\chi}_0 = \chi_0；L_0 = 0；S_0 = 0；\theta_0 = 1。这里 A^* 是 A 的伴随算子（共轭转置）。

迭代主体

对于 k = 1 到 N 的每次迭代，依次执行四个步骤：

低秩更新采用奇异值阈值化：

这里 \text{SV}_{\lambda_L}(\cdot) 表示对矩阵进行奇异值分解，然后对奇异值应用软阈值操作。具体地，若 M = U\Sigma V^* 是奇异值分解，则 \text{SV}_{\lambda}(M) = U \cdot \text{soft}(\Sigma, \lambda) \cdot V^*，其中软阈值作用在对角矩阵 \Sigma 的每个奇异值上。这个操作是核范数的近端算子，用于促进矩阵的低秩结构。

稀疏更新采用变换域软阈值化：

首先将残差 M_{k-1} - \hat{L}_{k-1} 变换到稀疏域 \mathcal{T}，然后应用软阈值，最后通过逆变换 \mathcal{T}^* 返回图像域。这实现了在变换域中促进稀疏性。

动量更新采用FISTA加速格式：

第一个公式更新动量参数，第二个公式计算外推点，将当前点沿着与上一步相同的方向进行外推，这是FISTA加速的核心。

数据一致性步骤：

这一步沿着数据保真项的负梯度方向更新，确保重建结果与观测数据保持一致。步长 \gamma 控制梯度下降的幅度。

重建结果的演变

从实验结果可以观察迭代过程中图像质量的改善。2D模式下，第1次迭代的重建图像模糊且存在明显伪影，经过50次迭代后对比度显著提升，解剖结构更加清晰。

3D模式同样显示出类似的改善趋势，50次迭代后采集伪影明显减少，图像质量大幅提高。

BOLD信号的时间分辨率分析

信号重建质量对比

BOLD信号是fMRI的核心测量量，反映神经活动引起的血氧水平变化。2D和3D采集模式重建的BOLD信号时间曲线显示了不同的特征。

在第1次迭代时，2D模式（左上）和3D模式（右上）的重建信号都存在较大噪声，橙色的重建信号与青色的期望信号之间有明显偏差。经过50次迭代后，两种模式都能较好地恢复周期性的BOLD信号模式。2D模式（左下）的时间范围约为0-140秒，3D模式（右下）的时间范围为0-350秒，后者由于更高的时间分辨率能够采集更长时间的信号。

重建信号呈现清晰的周期性波动，这与实验设计的block范式相对应。每个周期包含信号上升、平台和下降阶段，反映了神经激活开启和关闭时BOLD响应的血流动力学特性。

统计激活图

统计激活图通过z分数统计量来识别大脑中对刺激有显著响应的区域。z分数衡量每个体素的信号变化相对于噪声水平的强度，正值表示激活，负值表示抑制。

2D模式（左图）和3D模式（右图）的激活图均使用0.3作为z分数阈值来筛选显著激活的区域。颜色条范围从-11到11，红色/黄色区域表示正激活（z > 0.3），蓝色区域表示负激活（z < -0.3）。两种模式都能检测到视觉皮层区域的激活，红色轮廓标记了预期的激活区域边界。3D模式由于覆盖更大的空间范围，可以同时显示多个脑切片的激活分布。

分析模型的稀疏正则化

优化问题形式

分析模型的典型优化问题为：

与综合模型不同，这里直接对图像 x 进行优化，而非对稀疏系数。稀疏变换算子 T 作用在图像上，常见选择包括：小波变换，将图像分解为多尺度多方向的系数；有限差分算子，计算相邻像素的差值，对应全变分（Total Variation, TV）正则化；或者同时使用小波和TV的组合。

FDA批准的临床压缩感知MRI方法据推测与这种分析形式相关。然而，分析优化问题比综合形式更难求解，原因在于变换矩阵 T 出现在 \ell_1 范数内部，导致近端算子没有简单的解析解。

近端梯度方法的困难

对分析正则化问题应用近端梯度方法（PGM），迭代格式为：

第一步是标准的梯度下降，L = \|A\|_2^2 是数据保真项梯度的Lipschitz常数。第二步需要计算带有变换 T 的 \ell_1 范数的近端算子。

与综合模型中 \ell_1 范数近端算子有解析的软阈值解不同，\|Tx\|_1 的近端算子没有简单的闭式解。求解这个近端子问题本身需要使用内层迭代方法，通常涉及对偶问题的构造和求解。这构成了分析正则化的主要缺点：每次外层迭代都需要运行一个内层优化算法，计算成本显著增加。

由于这个原因，PGM、FPGM（即FISTA）和POGM等近端方法对于分析正则化问题的吸引力下降。在实际应用中，需要权衡分析模型更自然的建模方式与其更高的计算复杂度。一些替代方法包括ADMM（交替方向乘子法）等分裂算法，它们可以更高效地处理这类问题。

分析正则化的近似求解方法

分析正则化问题中变换矩阵 T 出现在 \ell_1 范数内部，导致近端算子难以计算。为绕过这一困难，有几种近似或替代策略。

角平滑方法

角平滑的核心思想是用光滑函数近似绝对值函数在零点处的尖角。具体地，用以下光滑近似替代绝对值：

其中 \epsilon > 0 是一个小的正数。当 |z| 远大于 \sqrt{\epsilon} 时，\sqrt{|z|^2 + \epsilon} \approx |z|；当 |z| 接近零时，该函数在零点附近变得光滑而非尖锐。这种近似与边缘保持正则化（edge-preserving regularization）的思想相似，在图像处理中用于保护边缘的同时去除噪声。

使用角平滑后，近端算子不再具有产生精确零值的阈值效应，因此无法真正诱导稀疏性。重建结果中的小系数会被压缩但不会精确为零，这与原始 \ell_1 正则化的稀疏促进特性有所不同。

惩罚方法

另一种策略是将分析正则化问题重新表述为：

这里引入了一个辅助的最小化问题来定义正则项 R_\alpha(x)。R_\alpha(x) 实际上等于Huber函数 \psi(Tx, \alpha)，Huber函数是绝对值函数的光滑近似，在原点附近是二次的，在远离原点处是线性的。这种方法相当于对原问题进行了角平滑处理，需要选择参数 \alpha 来控制平滑程度。

迭代重加权最小二乘方法（FOCUSS）也属于类似的角平滑类方法。

变量分裂方法

等价约束问题的构造

变量分裂是处理分析正则化的另一条重要途径。通过引入辅助变量 z，将原问题中的 \|Tx\|_1 替换为等价的约束优化问题：

这个问题与原问题完全等价：当约束 z = Tx 满足时，目标函数就退化为原来的形式。引入辅助变量的好处在于将变换 T 从 \ell_1 范数中分离出来，使得对 z 的优化子问题可以用简单的软阈值求解。

求解这类约束优化问题的算法包括：分裂Bregman算法、增广拉格朗日（AL）方法、交替方向乘子法（ADMM）以及Douglas-Rachford分裂方法。这些方法在数学上紧密相关，ADMM可以视为增广拉格朗日方法的一种特殊实现。

增广拉格朗日函数

对于约束问题，构造增广拉格朗日函数：

这里 \gamma \in \mathbb{C}^K 是拉格朗日乘子，用于强制约束 Tx = z；\mu > 0 是增广拉格朗日惩罚参数，它影响算法的收敛速度但不影响最终解 \hat{x}。增广项 \frac{\mu}{2}\|Tx - z\|_2^2 是对约束违反的二次惩罚，使得子问题更容易求解且改善收敛性。

引入缩放对偶变量 \eta \triangleq (1/\mu)\gamma，增广拉格朗日函数可以写成更简洁的形式：

这种缩放形式的优势在于将乘子和增广项合并为一个平方项，简化了后续的推导和实现。

增广拉格朗日方法在原始变量 x、z 的下降更新和对偶变量 \eta 的上升更新之间交替进行。

增广拉格朗日的迭代更新

各变量的更新规则

对 z 的更新：固定 x 和 \eta，对 z 最小化增广拉格朗日，得到简单的软阈值操作：

这是因为关于 z 的子问题形如 \min_z \beta\|z\|_1 + \frac{\mu}{2}\|z - (Tx_k + \eta_k)\|_2^2，其解正是对 Tx_k + \eta_k 应用软阈值，阈值为 \beta/\mu。

对 x 的更新：固定 z 和 \eta，对 x 最小化，这是一个二次问题，解为：

这个线性系统需要求解 (A'A + \mu T'T) 的逆与某个向量的乘积。直接求逆计算量大，实际中常用预条件共轭梯度（PCG）等迭代方法求解。

对 \eta 的更新：对偶变量沿着约束违反的方向上升：

这个更新使得当约束 Tx = z 被违反时，乘子会增长，从而在后续迭代中对约束违反施加更大的惩罚。

参数 \mu 的选择影响收敛速度，存在自适应调节 \mu 的方法。并行ADMM更新也是可行的，可以同时更新多个变量以加速计算。

并行MRI中的AL/ADMM

测量矩阵的结构

在并行MRI中，测量矩阵具有特殊结构 A = F_LC，其中 F_L \triangleq I_L \otimes F，\otimes 表示Kronecker积，I_L 是 L \times L 单位矩阵（L 是线圈数），F 是傅里叶变换矩阵。C_l 是第 l 个线圈的灵敏度对角矩阵，C 是所有线圈灵敏度的组合。

关键的矩阵性质是 F'F 的结构：对于笛卡尔采样，F'F 是循环矩阵；对于非笛卡尔采样，F'F 是Toeplitz矩阵。这些特殊结构可以利用FFT或DCT进行快速计算。

替代变量分裂策略

为了获得更简单的子问题更新，可以引入更多辅助变量：

这种多重分裂将原问题分解为多个简单的子问题。各子问题的更新特点如下：z 的更新是软阈值操作；x 的更新涉及对角矩阵 C'C，可以逐元素计算；v 的更新涉及 T'T，当 T 是循环或Toeplitz矩阵时可用FFT或DCT快速计算；u 的更新涉及 F_L'F_L，同样具有循环或Toeplitz结构。

这种方法的主要缺点是需要调节更多的增广拉格朗日惩罚参数，每个约束对应一个参数。条件数准则等方法可以帮助选择这些参数。

基于图像块的正则化

从像素到图像块

全变分正则化使用有限差分算子 T，正则项为 R(x) = \|Tx\|_1。有限差分操作等价于使用 2 \times 1 大小的图像块（patch）：水平差分用 [-1, 1] 滤波器，垂直差分用相应的垂直排列。这种最小的块只考虑相邻两个像素的差异。

使用更大的图像块可以提供更多的上下文信息，有助于区分信号和噪声。直观地说，更大的块能够捕捉更复杂的局部结构，而不仅仅是简单的梯度信息。这与卷积神经网络（CNN）的思想相似，CNN使用学习得到的卷积核在不同大小的感受野上提取特征。

基于图像块的正则化器可以采用综合模型或分析模型的形式。综合模型中，每个图像块被表示为字典原子的稀疏线性组合；分析模型中，对每个图像块应用分析变换后要求结果稀疏。这种方法在自然图像处理中取得了很好的效果，因为自然图像中的图像块往往具有低秩或稀疏的结构。

基于图像块的字典稀疏模型

核心假设与数学表示

基于图像块的稀疏模型假设：如果 x 是一幅合理的图像，那么从中提取的每个图像块都可以用字典原子的稀疏线性组合来近似表示：

其中 P_p 是图像块提取算子，从图像 x 中提取第 p 个图像块（共 P 个块）；D 是图像块字典，通常是过完备的；z_p 是对应于第 p 个图像块的稀疏系数向量。这是综合模型在图像块层面的应用。

从示意图可以看到，CT图像中不同位置的图像块被提取出来，每个块呈现不同的局部结构特征。有的块主要是均匀区域，有的块包含边缘或纹理。字典的作用是提供一组基本模式（原子），使得各种类型的图像块都能用少量原子的组合来表示。

图像块综合模型的正则化

优化问题形式

图像块综合模型使用字典原子的稀疏线性组合来表示每个图像块，即 P_p x \approx D z_p。其中 P_p \in \{0, 1\}^{d \times N} 是提取算子，从 N 像素的图像 x 中提取第 p 个 d 像素的图像块；D \in \mathbb{C}^{d \times J} 是包含 J 个原子的图像块字典；z_p \in \mathbb{C}^J 是第 p 个图像块的稀疏系数向量。

基于这种模型的正则化器定义为：

正则项 R(x) 的含义是：对于给定的图像 x，寻找最优的稀疏系数集合 \{z_p\}，使得每个图像块 P_p x 与其字典表示 D z_p 的误差最小，同时系数 z_p 保持稀疏。第一项 \frac{1}{2}\|P_p x - D z_p\|_2^2 是表示误差，第二项 \alpha\|z_p\|_1 是稀疏惩罚。

优化算法

这个双层优化问题通常用交替最小化算法求解。对 x 的更新：固定所有 \{z_p\}，关于 x 的子问题是二次的，可以用预条件共轭梯度（PCG）方法求解。对 z_p 的更新：固定 x，每个图像块的稀疏编码子问题可以独立求解，这是标准的LASSO问题，可以用POGM等近端方法求解。

这种方法的一个计算挑战是存储需求：基本实现需要存储所有 P 个图像块的 J 维系数向量，总共 JN 个系数值（假设图像块有重叠，P 与 N 同阶）。

图像块分析模型的正则化

分析模型的假设与形式

图像块分析模型假设对每个图像块施加稀疏变换后的结果是稀疏的，即 T P_p x 趋向于稀疏。这里 P_p \in \{0, 1\}^{d \times N} 仍是图像块提取算子，T \in \mathbb{C}^{K \times d} 是图像块稀疏化变换，典型选择如离散余弦变换（DCT）。

对应的正则化器形式为：

与综合模型相比，这里变换 T 作用在图像块上，辅助变量 z_p 近似变换后的系数 T P_p x，而非直接作为图像块的表示系数。

优化与实现

同样采用交替最小化算法。对 x 的更新是二次问题，对应的Hessian矩阵为 A'A + \beta \sum_p P_p' T' T P_p。当变换 T 是酉矩阵（如DCT）时，T'T = I，Hessian简化为 A'A + \beta \sum_p P_p' P_p，计算更加简单。

对 z_p 的更新仅需简单的软阈值操作：

这比综合模型中的稀疏编码更容易，因为软阈值有解析解，无需迭代求解。

为避免存储所有 KN 个系数值 \{z_p\}，可以使用累加器技巧。定义：

这样在更新 x 时直接累加各图像块的贡献，无需显式存储中间系数。

稀疏正则化方法总结

图像模型与代价函数

本节介绍的图像模型和对应的代价函数可以从多个维度分类。图像稀疏综合模型假设图像可以表示为基向量（如小波基）的稀疏线性组合。图像变换/分析稀疏模型假设对图像施加某种变换（如有限差分对应全变分、小波变换等）后结果是稀疏的。图像块字典稀疏模型在图像块层面应用综合模型，每个块用字典原子稀疏表示。图像块变换稀疏模型在图像块层面应用分析模型。卷积稀疏模型使用卷积结构的字典或变换，可以视为全局一致的图像块模型。

算法组件

求解这些优化问题的算法主要基于以下组件。基于梯度的方法如共轭梯度（CG），用于求解光滑的二次子问题。近端操作包括软阈值和硬阈值，分别对应 \ell_1 和 \ell_0 范数的近端算子。交替最小化在多个变量之间轮流优化，每次固定其他变量只优化一个。对偶方法包括增广拉格朗日（AL）、ADMM和原始-对偶算法，通过引入对偶变量将困难的约束问题转化为更易求解的形式。

自适应正则化与字典学习

稀疏模型的关键组件

不同类型的稀疏模型需要不同的关键组件：基于图像块的综合模型需要字典 D；基于图像块的分析模型需要稀疏化变换 T；卷积模型需要滤波器组 \{h_k\}。这些组件的选择直接影响重建质量。

组件获取的方式

获取这些组件有多种途径。手工设计的数学模型如离散余弦变换（DCT）、小波变换等，具有良好的数学性质和快速算法，但可能无法最优地适应特定类型的图像。群体自适应方法从高质量的代表性训练数据（如完全采样的参考图像）中学习模型，学到的字典或变换可以捕捉特定成像模态或解剖结构的统计特性。患者自适应方法针对每个患者的数据，联合优化图像 x、稀疏系数 Z 以及模型组件（D 或 T 或 \{h_k\}），这种方法可以最好地适应个体差异，但计算量更大。还可以考虑群体自适应和患者自适应方法的混合策略。

对于自适应方案，必须对学习得到的组件施加约束，否则优化问题可能是病态的。例如，字典的列需要归一化以避免与稀疏系数的尺度模糊。

自适应正则化的设计选择

三个维度的选择

自适应正则化的设计涉及三个主要维度的选择。数据维度决定模型的适应范围：群体自适应方法使用外部训练数据学习通用模型；患者自适应方法针对当前重建任务学习专用模型，特别适用于动态MRI等场景，因为同一序列的不同帧之间有强相关性。

空间结构维度决定如何建模图像的局部特性：基于图像块的模型将图像分割为重叠的小块分别处理；卷积模型使用全局一致的滤波器在整幅图像上滑动。

正则化器形式维度决定稀疏性假设的数学表达：综合（字典）方法假设图像块是字典原子的稀疏组合；分析（稀疏化变换）方法假设变换后的系数是稀疏的。

这三个维度的不同组合产生了多种具体方法，实际选择需要根据具体应用场景、可用数据和计算资源来权衡。

本课程关注的方法

在众多可能的组合中，本课程重点关注的方法选择是：数据方面采用患者自适应方法，空间结构方面采用基于图像块的模型，正则化器形式方面采用综合（字典）方法。这种组合的优势在于可以充分利用当前患者数据的特性，通过学习适应性字典来捕捉个体特异的图像结构，同时图像块模型提供了足够的灵活性来处理局部变化。分析方法（稀疏化变换）作为对比也会讨论。

字典学习MRI重建

字典盲重建问题

字典盲MR图像重建是指在重建过程中同时学习适应性字典，优化问题的形式为：

其中正则项 R(x) 定义为：

这里 P_m 是提取第 m 个图像块的算子，共有 M 个图像块；D 是待学习的字典，约束在可行集 \mathcal{D} 内（通常要求列归一化）；z_m 是第 m 个图像块的稀疏系数；Z = [z_1 \ldots z_M] 是所有系数的矩阵。

这里使用 \ell_0 范数 \|z_m\|_0 而非 \ell_1 范数作为稀疏惩罚。\ell_0 范数直接计算非零元素个数，虽然非凸，但在自适应字典学习场景下，由于整个问题本身已经是非凸的（字典 D 和系数 Z 的乘积形式），使用 \ell_0 反而更直接，且实验表明可以获得更高的信噪比。

交替嵌套最小化算法

由于优化问题涉及三个变量 x、D 和 Z，采用交替最小化策略依次更新。

固定 x 和 D，更新 Z：每个图像块的稀疏系数 z_m 可以独立更新。由于使用 \ell_0 范数，最优解通过硬阈值操作获得：保留绝对值最大的若干系数，其余置零。Z 矩阵的每一行（对应同一字典原子在所有图像块中的系数）可以顺序更新。

固定 x 和 Z，更新 D：使用SOUP-DIL（Sum of Outer Products Dictionary Learning）算法。这个算法利用了字典学习问题的特殊结构，可以高效求解。

固定 Z 和 D，更新 x：关于 x 的子问题是二次的，目标函数形如 \frac{1}{2}\|Ax - y\|_2^2 + \frac{\beta}{2}\sum_m \|P_m x - D z_m\|_2^2。对于单线圈笛卡尔MRI，利用傅里叶变换的性质可以用FFT高效求解。对于非笛卡尔采样或并行MRI，需要使用共轭梯度（CG）迭代方法。

整个问题是非凸的，但算法具有单调下降性，即每次迭代目标函数值不增。在一定条件下存在收敛性理论保证。

二维字典学习MRI实验结果

实验设置与可视化

实验使用 6 \times 6 大小的图像块，字典维度为 D \in \mathbb{C}^{36 \times 144}，即36维的图像块空间中有144个原子，是4倍过完备的字典。初始字典 D_0 可以选择DCT矩阵或随机矩阵。

从重建结果可以看到：完全采样的参考图像清晰显示了三个圆形结构；零填充重建（直接对欠采样k空间补零后逆变换）存在明显的混叠伪影；SOUP-DILLO-MRI方法的重建结果伪影大幅减少，图像质量接近完全采样。

采样模式显示了2.5倍加速下的笛卡尔欠采样，呈现规则的条纹图案。初始字典显示为规则的DCT基函数图案。学习得到的字典（分别显示实部和虚部）呈现出更丰富的结构，包含了边缘、纹理等适应于具体图像的特征。

收敛性能与稀疏惩罚比较

收敛曲线显示PSNR（峰值信噪比，相对于完全采样图像）随迭代次数的变化。SOUP-DILLO MRI（红色实线）使用 \ell_0 范数惩罚，SOUP-DILLI MRI（蓝色虚线）使用 \ell_1 范数惩罚。两种方法都在约40次迭代后基本收敛，最终PSNR约为37dB。使用 \ell_0 范数的SOUP-DILLO比使用 \ell_1 范数的SOUP-DILLI获得略高的SNR。由于自适应字典学习本身就是非凸问题，使用非凸的 \ell_0 范数并不会增加本质困难，反而可能带来性能提升。

多图像多方法性能对比

在多种测试图像上比较不同方法的PSNR性能。测试图像包括：(a) 膝关节、(b) 脑部、(c) 仿真模型、(d) 心脏、(e) 脊柱、(f) 脑部、(g) 另一脑部图像。采样方式包括笛卡尔采样和2D随机采样，加速倍数从2.5x到7x不等。

对比方法包括：0-fill（零填充）、Sparse MRI、PANO、DLMRI（基于K-SVD的字典学习）、SOUP-DILLI（\ell_1惩罚）和SOUP-DILLO（\ell_0惩罚）。表格数据显示，在大多数情况下SOUP-DILLO获得最高或接近最高的PSNR。例如图像(b)在2.5x笛卡尔采样下，零填充PSNR为27.7dB，Sparse MRI为31.6dB，而SOUP-DILLO达到42.3dB，提升显著。

误差图比较

通过误差图可以更直观地比较不同方法的重建质量。误差图显示重建图像与完全采样参考图像之间的绝对差值，颜色越深表示误差越小。DLMRI方法的误差图显示存在明显的残余结构；PANO方法误差有所减少但仍可见；FDLCP方法存在水平条纹状伪影；SOUP-DILLO方法的误差图最为均匀，整体误差最小。

总结而言，基于自适应字典学习的2D静态MR重建方法可以从欠采样数据获得高质量图像。SOUP-DIL算法比传统的K-SVD方法（如DLMRI使用的）更快速，同时保持收敛性保证。

深度学习MRI重建

本节进入课程的核心内容——将深度学习应用于MRI图像重建。这一领域近年来发展迅速，多位研究者做出了重要贡献，包括在大规模公开数据集上对重建网络进行基准测试、开发用于非笛卡尔MRI的密度补偿展开网络（NC-PDNet）、以及基于深度学习的离共振校正方法等。

深度学习的预期收益

MRI成像流程包括数据采集和图像重建两个主要耗时环节。传统MR扫描的采集时间和重建时间都较长。引入压缩感知、非笛卡尔成像和SPARKLING等加速采集技术后，采集时间显著缩短，但迭代重建算法的计算负担使得重建时间反而增加。深度神经网络的引入有望同时缩短采集时间和重建时间：通过学习从欠采样数据到高质量图像的直接映射，可以在保持或提升图像质量的同时大幅加速重建过程。

深度学习方法的分类框架

基于优化的分类

MRI重建的基本问题可以表述为：给定观测模型 y = Ax + \epsilon，求解优化问题：

其中 f(x) 是数据保真项，g(x) 是正则项。深度学习方法可以按照与这一优化框架的关系进行分类。

单域方法（Single domain）直接学习从测量数据 A^H y 到重建图像 \hat{x} 的映射，用神经网络 G_\theta 替代整个重建过程。这种方法简单直接，但完全忽略了MRI物理模型的知识，且训练成本高。

展开方法（Unfolded）将传统迭代优化算法展开为固定深度的网络结构。每一层对应一次迭代，包含梯度下降步骤 I - \eta\nabla f 和学习的正则化模块。这种方法的优点是可以学习先验信息，并且物理模型可微分；缺点是训练代价很高，且在分布外数据上容易失效。

即插即用方法（Plug-and-Play）将预训练的去噪网络作为近端算子 \text{prox}_g \simeq G_\theta 嵌入到迭代算法中。网络与测量矩阵 A 独立训练，因此对不同的采样模式具有通用性。但这种方法面临分布偏移的问题，即训练数据与实际应用数据分布不一致时性能下降。

基于能量的模型（Energy-based models）学习一个能量函数 E_\theta(x): \mathbb{C}^m \to \mathbb{R}，对应的概率分布为：

其中 Z_\theta 是配分函数。这种方法同样具有网络独立训练、对测量矩阵 A 不敏感的优点，但也存在分布偏移问题。

从单域方法到基于能量的模型，方法的设计越来越多地受到优化理论的启发。

基于处理域的分类

根据神经网络作用的数据域，深度学习MRI重建方法可以分为四类。

图像域学习（Image-domain learning）：首先对k空间数据进行逆傅里叶变换（IFT）得到混叠图像，然后用深度神经网络在图像域进行去伪影和增强处理。

跨域学习（Cross-domain learning）：网络同时在k空间和图像空间操作，通过傅里叶变换（FT）和逆变换在两个域之间切换。这种方法可以在k空间进行数据一致性校正，在图像空间进行去噪和结构恢复。

域变换学习（Domain-transform learning）：神经网络直接学习从k空间到图像空间的变换，隐式地包含了傅里叶变换的功能。

传感器域学习（Sensor-domain learning）：首先在k空间用神经网络进行处理（如数据补全或去噪），然后再进行逆傅里叶变换得到图像。

其中跨域学习方法结合了k空间和图像空间处理的优势，是目前研究的重点方向。

图像域单域学习

基本架构

图像域学习的流程是：k空间数据 → 逆傅里叶变换 → 深度神经网络 → 重建图像。最简单的形式直接对IFT结果应用深度网络，但这样网络无法获得关于采样模式等物理信息的知识。

更合理的做法是利用伴随算子 A^H 在图像空间构建初始估计 x_0 = A^H y，然后用U-Net等网络进行增强。U-Net的编码器-解码器结构配合跳跃连接，可以在多个尺度上提取特征并恢复细节。

这种方法的优点是简单且推理速度快。缺点是欠采样导致的混叠伪影在空间上是广泛分布的，需要较深的网络才能有效去除，因为网络需要足够大的感受野来捕捉远距离的伪影关联。

图像域学习的风险

图像域学习方法存在对未见病理的鲁棒性问题。当训练数据中不包含某种病变类型时，网络可能无法正确重建这些结构，甚至可能错误地将病变当作伪影去除。这一问题在ISMRM 2020数据采样与图像重建研讨会上被深入讨论。

相关研究表明，无模型的深度学习重建方法虽然速度快，但在泛化到训练时未见过的病理情况时可能出现严重失败。研究发现，包含MR采集模型及其伴随算子迭代应用的方法在泛化性能上更优，这与迭代重建方法的思想一致。

另一项发表于PNAS 2020的研究系统性地揭示了深度学习图像重建中的不稳定性问题，表明微小的输入扰动可能导致重建结果的显著变化，这对医学应用构成潜在风险。

图像域学习的临床风险案例

一个具体的临床案例展示了图像域学习方法的潜在危险。该案例涉及间变性星形细胞瘤（anaplastic astrocytoma），这是一种罕见的恶性脑肿瘤。对比参考图像、SPARSE-SENSE重建和深度学习（DL）重建的结果，三者均来自相同的4倍加速回顾性欠采样数据。

从全体积平均绝对误差（MAE）来看，深度学习方法（0.0173）优于SPARSE-SENSE方法（0.0263）。然而，仔细观察肿瘤区域可以发现，深度学习重建在肿瘤附近的区域出现了明显的失真，未能正确重建病变结构。这说明全局误差指标可能掩盖局部关键区域的重建失败，而这些区域恰恰是临床诊断最关注的部分。

k空间单域学习

基本思想与架构

k空间学习的流程是：k空间数据 → 深度神经网络 → 处理后的k空间 → 逆傅里叶变换 → 图像。为了构建更有信息量的模型，可以利用伴随算子 \mathcal{A}^H：

神经网络 f_\theta 在k空间中操作，然后通过伴随算子映射到图像空间。这种方法可以理解为非线性版本的GRAPPA并行成像方法，因此有时被称为非线性GRAPPA。

这种方法的优点包括：简单且推理速度快；可以实现无数据库训练，即仅从当前扫描的自校准数据（ACS）中学习网络参数，类似于传统GRAPPA从ACS线估计卷积核的方式。缺点是k空间中的局部操作可能难以表示图像空间中的局部特征，因为k空间中的每个点都包含整个图像的全局信息。

商业应用实例

西门子医疗的Deep Resolve技术是k空间学习在7T MRI上的商业化应用。对比传统重建和Deep Resolve重建：T2 TSE序列，传统方法使用PAT 2加速，分辨率0.2×0.2×2.0 mm³，采集时间4:42分钟；Deep Resolve使用PAT 4加速，相同分辨率，采集时间仅2:14分钟，快了52%。从图像质量来看，Deep Resolve在更高加速倍数下仍能保持清晰的解剖细节。

模型无关学习

最极端的端到端学习方式是完全丢弃物理模型，直接学习从k空间到图像的映射：

网络输入是原始k空间数据，输出是重建图像，不使用任何关于采样模式、傅里叶变换或线圈灵敏度的先验知识。

这种方法的主要缺点是：没有利用已知的物理知识，网络需要从数据中重新学习傅里叶变换等基本操作，效率低下；缺乏可扩展性，当图像尺寸、采样模式或线圈配置改变时，需要重新训练整个网络。

AUTOMAP变换学习

方法框架

AUTOMAP（Automated Transform by Manifold Approximation）是域变换学习的代表性方法，于2018年发表于Nature。该方法通过神经网络学习从传感器数据到图像的变换，目标是最小化KL散度：

其中 P 是真实的联合分布，Q 是模型分布，Y 是测量数据，X 是真实图像，\tilde{X} 是重建图像。

从信号处理角度，传统重建链可以表示为：

其中 \phi_x 是图像域的变换，g 是测量过程，\phi_y 是测量域的变换。AUTOMAP用神经网络直接学习这个复合映射的逆。

网络架构包括：首先是全连接层（FC1, FC2, FC3）处理复数传感器数据，将其映射到适当的特征空间；然后是卷积层（C1, C2）进一步处理并输出最终图像。全连接层的输入维度为 2n^2（实部和虚部），输出图像尺寸为 n \times n。

实验结果

AUTOMAP在多种成像模式上进行了测试，包括：Radon投影（CT成像）、螺旋非笛卡尔傅里叶采样、欠采样傅里叶采样、以及存在误差的傅里叶采样。

对于Radon投影重建，AUTOMAP达到SNR 33.8dB，RMSE 2.6%，显著优于传统方法（SNR 14.2dB，RMSE 5.3%）。螺旋非笛卡尔采样的结果更为突出，AUTOMAP的SNR为42.7dB，RMSE 1.7%，而传统方法仅有SNR 13.8dB，RMSE 5.0%。欠采样傅里叶情况下，AUTOMAP同样表现优异（SNR 59.3dB vs 43.5dB）。即使在存在测量误差的情况下，AUTOMAP也能产生合理的重建结果。

AUTOMAP的不稳定性问题

对抗扰动的脆弱性

尽管AUTOMAP展示了令人印象深刻的重建质量，但后续研究揭示了其严重的稳定性问题。考虑两种类型的扰动：

第一行显示了对k空间数据添加小扰动 e_1, e_2, e_3 后的重建结果 \Psi(Ax + e_i)。可以看到，虽然扰动很小，但重建图像出现了明显的伪影和失真，最后一幅图像几乎完全被噪声淹没。

第二行使用FIRENET方法对相同的扰动进行重建 \Phi(Ax + \tilde{e}_i)。FIRENET是一种设计上考虑了稳定性的网络架构，对相同类型的扰动展现出更好的鲁棒性，重建结果保持了基本的图像结构。

这一对比说明，高精度与高稳定性之间可能存在权衡。相关理论工作（PNAS 2022）进一步分析了计算稳定且精确的神经网络的固有困难，指出某些类型的逆问题可能不存在同时满足精度和稳定性要求的神经网络解。

跨域学习

核心思想

跨域学习结合了k空间和图像空间处理的优势。基本架构包括：k空间数据首先通过解析重建（Analytic Recon）得到初始图像，然后神经网络同时利用k空间数据和图像域信息进行迭代优化，最终输出通过前向求解器（Forward Solver）保证与物理模型一致。

跨域学习的优点包括：基于物理的方法，同时利用k空间数据一致性和图像域先验信息；与优化方法有可解释的联系，网络结构可以理解为展开的迭代算法；在MRI重建挑战赛中取得最佳结果，证明了方法的有效性。

缺点是计算量较大，因为需要多次迭代（对应网络的多层），每次迭代都涉及FFT或NUFFT操作。对于非笛卡尔采样，NUFFT的计算成本更高。

展开重建算法

从优化迭代到网络结构

展开（Unrolling）方法的核心思想是将迭代优化算法的有限次迭代展开为神经网络的层结构。以压缩感知MRI的典型优化问题为例：

其中 \mathcal{F}_K 是欠采样傅里叶算子，y 是k空间测量数据，\Psi 是稀疏变换（如小波变换），\lambda 是正则化参数。

用ISTA求解该问题，每次迭代包含两步。首先是梯度下降步骤（半迭代）：

这一步沿着数据保真项的负梯度方向移动，\eta 是步长，\mathcal{F}_K^* 是 \mathcal{F}_K 的伴随算子。然后是近端步骤：

这一步在变换域应用软阈值，然后逆变换回图像域。

ISTA的图表示

将ISTA迭代写成更紧凑的形式：

第一步是数据一致性步骤（Data Consistency, DC），确保重建结果与测量数据一致；第二步是近端步骤（Proximal, Prx），应用正则化。

这个迭代过程可以表示为计算图：从初始估计 x_0 开始，交替通过DC模块和Prx模块，每次迭代产生新的估计。图中箭头表示数据流向，DC模块利用测量数据 y 和采样模式信息，Prx模块应用稀疏先验。

展开为固定深度网络

将 N 次迭代展开，得到链式结构：

整个链条重复 N 次DC-Prx单元。在传统ISTA中，每次迭代使用相同的操作和参数。展开后可以让不同迭代使用不同的参数，甚至不同的操作。

用神经网络替代近端算子

关键的创新是用可学习的神经网络 f_{\theta_n} 替代固定的近端算子：

DC步骤保持不变，仍然使用物理模型（傅里叶变换和采样模式）确保数据一致性；Prx步骤被神经网络替代，网络参数 \theta_n 通过端到端训练学习。不同迭代的网络可以有不同的参数（\theta_1, \theta_2, \ldots, \theta_N），也可以共享参数以减少模型大小。

展开后的整体结构：

整个链条构成一个完整的神经网络，可以端到端训练。损失函数通常定义为最终输出 x_N 与真实图像之间的距离（如MSE或感知损失）。

通用重建模型框架

双域网络架构

通用重建模型同时在k空间和图像空间进行处理。输入包括采样模式 \Omega 和测量数据 y。处理流程包含两条并行路径。

k空间路径：测量数据 y 首先经过数据一致性模块，然后通过k空间校正网络（k-space net），可以进行数据插值或去噪。

图像路径：通过傅里叶变换 F^T（转置，即逆变换）将k空间数据映射到图像域，图像校正网络（image net）在图像域进行去伪影和增强处理。

两条路径通过傅里叶变换 F 和逆变换 F^T 相互连接，形成迭代结构。最终输出重建图像 \hat{x}。

图像校正网络通常采用全卷积网络，可以是残差结构（学习图像与目标之间的残差）以便于训练。

不同展开网络的变体

展开网络的设计涉及三个主要选择：展开哪种优化算法、神经网络 f_\theta 的架构、以及展开的迭代次数 N。

KIKI-Net在k空间和图像空间都使用CNN，交替在两个域进行处理，充分利用两个域的互补信息。

U-Net方法仅进行图像域校正，只需 N=1 次迭代，结构简单但可能无法充分利用k空间信息。

CascadeNet是多个图像校正网络的级联，每个网络之间插入数据一致性层，是一种简化的展开结构。

PDNet（Primal-Dual Network）基于原始-对偶优化算法展开，同样仅进行图像校正，但增加了来自前次迭代的记忆机制（memory），可以利用迭代历史信息。

重建网络基准测试

实验设置

在大规模公开数据集fastMRI上对不同重建网络进行了系统性比较。评估指标为PSNR，在200个验证体积上计算平均值。比较的模型基于三个维度的不同选择：展开的优化算法、网络架构 f_\theta、以及迭代次数 N。

定量结果

各方法的PSNR性能如下：Zero-filled（零填充基线）为29.61 dB；KIKI-net为31.38 dB；U-net为31.78 dB；Cascade net为31.97 dB；PD-net达到最高的32.15 dB。

从结果可以看出，所有深度学习方法都显著优于零填充基线，提升约2-2.5 dB。基于原始-对偶展开的PD-net表现最好，说明结合优化理论的网络设计确实有效。即使是简单的单次U-net校正也能获得不错的性能，但迭代结构（如CascadeNet和PD-net）可以进一步提升。KIKI-net虽然同时在两个域处理，但性能略低于纯图像域方法，可能是因为k空间网络的设计还有改进空间。

重建网络可视化比较

膝关节MRI重建结果

从膝关节MRI的重建结果可以直观比较各方法的性能。参考图像显示清晰的骨骼和软组织结构。零填充重建存在明显的混叠伪影，表现为图像模糊和条纹状干扰。KIKI-net的重建已经显著改善，但仍有轻微伪影残留。U-net、Cascade-net和PD-net的重建质量逐步提升，其中PD-net最接近参考图像。

误差图（第二行）更清晰地显示了各方法的重建误差分布。零填充的误差图显示广泛分布的高强度误差。各深度学习方法的误差图中，误差逐渐减小且分布更均匀，PD-net的误差图整体最暗，表明重建误差最小。

相关代码和预训练模型权重已开源，分别托管在GitHub（github.com/zaccharieramzi/fastmri-reproducible-benchmark）和HuggingFace（huggingface.co/zaccharieramzi）平台。

XPDNet展开模型

基本架构

XPDNet是一种基于原始-对偶算法展开的深度学习重建网络。其基本结构为：

为了更清晰地描述网络结构，将数据一致性步骤（DC）和神经网络模块（f_{\theta}）组合成一个迭代单元（Iteration Unit, IU）。整个网络由 N 个IU串联组成：

迭代单元的详细结构

每个迭代单元内部包含数据一致性计算和图像校正网络。以第一个IU为例，输入为当前估计 x_0，处理流程如下：

首先计算数据一致性项。将当前图像估计通过前向算子 \mathcal{A} 映射到k空间，得到预测的测量值 \mathcal{A}x_0。然后计算与实际测量 y 的残差 (\cdot - y)，再通过伴随算子 \mathcal{A}^H 映射回图像空间，得到梯度方向的校正项 \mathcal{A}^H(\mathcal{A}x_0 - y)。

然后进行特征融合。将原始图像 x_0 与数据一致性校正项通过拼接（concat）操作合并为多通道输入。这样网络可以同时利用当前图像估计和数据一致性信息。

最后通过图像校正网络 f_{\theta_1} 处理拼接后的特征，输出更新后的图像估计 x_1。网络采用多级小波CNN（MWCNN）架构，该架构通过小波变换进行多尺度特征提取，特别适合图像恢复任务。

灵敏度图精炼

在多线圈并行MRI中，线圈灵敏度图（Sensitivity Maps）的准确性直接影响重建质量。XPDNet引入了灵敏度图精炼模块（Smaps refiner），使用U-Net对初始估计的灵敏度图进行优化。

灵敏度图精炼网络以初始灵敏度估计 S 为输入，输出精炼后的灵敏度图。该U-Net对所有线圈的灵敏度图使用相同的权重，即共享网络参数，这既减少了模型参数量，又保证了不同线圈处理的一致性。精炼后的灵敏度图被用于后续的数据一致性计算中。

这种设计来源于2019年fastMRI挑战赛的获胜方案，证明了灵敏度图精炼对提升重建质量的有效性。

展开优化算法的完整框架

多线圈非笛卡尔MRI重建

展开优化的核心思想是：将若干次已知优化算法的迭代展开为深度神经网络的层结构，从而解决病态逆问题。

对于多线圈并行MRI，完整的重建问题涉及多个组件：非均匀采样模式、NUFFT算子（Non-Uniform Fast Fourier Transform）、多个接收线圈的灵敏度图、k空间测量数据、以及稀疏正则化（如 \ell_1、Group LASSO等配合小波分解）。

优化问题的形式为：

其中 L 是线圈数量，\mathcal{F}_\Omega 是在采样模式 \Omega 上的傅里叶变换（对于非笛卡尔采样是NUFFT），S_\ell 是第 \ell 个线圈的灵敏度对角矩阵，y_\ell 是第 \ell 个线圈的k空间数据，R(x) 是正则项。

近端梯度下降的展开

用近端梯度下降求解上述问题，每次迭代的数据一致性步骤为：

其中 \tau_n 是第 n 次迭代的步长。这个表达式计算所有线圈贡献的加权和，然后沿负梯度方向更新图像估计。

近端步骤用神经网络替代：

展开后的网络结构为：传统ISTA的 N 次DC-Prx迭代被展开为包含 N 个IU的深度网络，每个IU内的近端算子被替换为可学习的神经网络 f_{\theta_n}。

灵敏度图估计模块

处理流程

灵敏度图估计（Sensitivity Map Estimation, SME）是多线圈重建的关键步骤。XPDNet采用基于深度学习的SME方法，处理流程如下：

输入k空间数据首先通过ACS掩模（Auto-Calibration Signal Mask）提取中心区域的全采样数据。ACS区域通常是k空间中心的一小块完全采样区域，包含低频信息，足以估计线圈灵敏度的平滑变化。

对ACS数据进行逆傅里叶变换（IFT），得到每个线圈的低分辨率图像。这些图像包含了各线圈的灵敏度加权信息。

U-Net对低分辨率线圈图像进行处理，去除噪声并增强灵敏度图的平滑性。同一个U-Net以相同的权重应用于所有线圈的粗灵敏度图，实现参数共享。

最后除以RSS（Root Sum of Squares）图像进行归一化。RSS是所有线圈图像幅值平方和的平方根，代表组合图像的强度分布。除以RSS后得到归一化的灵敏度图，满足 \sum_\ell |S_\ell|^2 = 1 的约束。

输出为估计的灵敏度图集合，每个线圈对应一幅复数灵敏度图，显示该线圈在不同空间位置的接收灵敏度。从示例图中可以看到，不同线圈的灵敏度图呈现不同的空间分布模式，反映了线圈阵列中各单元的几何配置。

这种灵敏度图估计方法是2019年fastMRI挑战赛获胜方案的核心组件之一。

2020年fastMRI挑战赛

挑战赛目标与设置

2020年fastMRI挑战赛由Facebook AI和NYU Langone Health联合举办，目标是建立一个国际性基准测试平台来评估深度学习在脑部MR图像重建中的解决方案。挑战赛的采集设置贴近临床实际，采用多线圈采集和多种对比度，包括T1加权、T1增强后（T1POST）、T2加权和FLAIR序列。

训练集规模较大，共包含6970次脑部扫描，原始k空间数据约1.5 TB，其中3001次扫描来自1.5T设备。数据按不同对比度和用途进行划分：训练集包含T1（498）、T1POST（949）、T2（2678）、FLAIR（344）共4469个样本；验证集共1378个样本；测试集分为4倍加速（281个）和8倍加速（277个）两个版本；挑战集同样分为4倍和8倍加速版本。

此外还设置了迁移学习赛道（Transfer Track），使用GE和Philips设备采集的数据测试模型的跨厂商泛化能力。GE数据共211个样本，Philips数据共118个样本（无T1POST序列）。

数据处理与评估

对于多线圈数据，每个线圈 \ell = 1, \ldots, L 的图像通过逆傅里叶变换从k空间恢复：

最终的组合图像采用RSS（Root Sum of Squares）方法计算：

RSS组合将各线圈图像的幅值平方求和后开方，是多线圈MRI中常用的图像组合方法，无需显式估计线圈灵敏度。

fastMRI挑战赛定量结果

放射科医师评估

挑战赛的评估由放射科医师进行，使用质量排名和Likert评分量表。评估维度包括：总体排名（Rank）、伪影程度（Artifacts）、清晰度（Sharpness）和对比噪声比（CNR）。分数越低表示性能越好。

4倍加速赛道结果显示：AIRS Medical团队获得最佳成绩，各项指标均为1.36-1.53；Nspin团队（即NeuroSpin/CEA团队，使用XPDNet）排名第二，各项指标为1.72-1.94；ATB团队排名第三。

8倍加速赛道结果中：AIRS Medical仍然领先（1.28-1.94）；Nspin团队保持第二（2.25-2.72）；ATB团队第三。

从重建图像的可视化比较可以看到，Ground Truth显示T1POST、T2和FLAIR三种对比度的参考图像。AIRS Medical的重建结果SSIM值最高（0.933, 0.946, 0.864）。ATB团队的结果略低（0.907, 0.936, 0.847）。NeuroSpin（Nspin）团队的结果（0.904, 0.935, 0.836）在学术研究机构中排名第一，总体排名第二。

这一成绩表明XPDNet在大规模公开数据集上具有竞争力，特别是在学术研究机构的提交中表现最优。

7T高分辨率数据的鲁棒性测试

跨分辨率泛化

为测试XPDNet的泛化能力，将在fastMRI数据（1.5T/3T，常规分辨率）上训练的模型应用于7T超高场强的高分辨率数据。7T MRI可以获得更高的空间分辨率和信噪比，但数据特性与训练数据存在显著差异。

测试结果显示：左侧为T2加权、GRAPPA采集的7T参考图像，清晰显示了大脑的精细解剖结构，包括小脑的详细纹理。右侧为XPDNet重建结果，该模型在4倍加速（R=4）数据上训练。

从对比可以观察到，XPDNet能够产生合理的重建结果，但存在两个明显问题：分辨率较低，细节不如参考图像清晰；小脑区域缺失或模糊。这些问题反映了训练数据与测试数据之间的域偏移：fastMRI训练数据可能不包含小脑区域，或者小脑在训练数据中的表示不足；分辨率差异导致网络学到的特征无法完全适应高分辨率数据的细节恢复需求。

精度与稳定性的权衡

幻觉问题

展开神经网络虽然比纯端到端方法更具物理约束，但仍然容易产生幻觉（hallucinations），即重建出训练数据中不存在或与真实结构不符的伪结构。

实验展示了这一现象：上排显示不同输入条件下的脑部图像，x_{\text{br}} 表示大脑区域，x_{\text{th}} 表示丘脑区域，x_{\text{mi}} 表示中间结构。当输入完整信号 x_{\text{br}} + x_{\text{th}} + x_{\text{mi}} 时，重建网络 \Psi 产生正确的结果。

下排显示重建结果：当输入缺少某些结构时，网络可能产生幻觉。\Psi(\mathcal{A}(x_{\text{br}} + x_{\text{th}})) 的重建中，红色箭头指示的位置出现了原本不存在的伪结构。\Psi(\mathcal{A}x_{\text{br}}) 的结果同样显示了类似的幻觉现象。这说明网络倾向于根据训练数据的统计特性填补缺失信息，但这种填补可能与实际情况不符。

精度-幻觉屏障

理论分析表明存在精度与幻觉之间的基本权衡。右侧示意图显示：横轴为误差率（Error rate），纵轴为幻觉概率（Chance of hallucinations）。曲线呈现负相关趋势，表明降低误差率会增加幻觉风险，反之亦然。这条曲线代表了精度-幻觉屏障（accuracy-hallucination barrier）。

产生这一现象的根本原因在于：大多数先验（priors）被训练来恢复欠采样数据，或者恢复属于测量矩阵零空间 \mathcal{N}(A) 的数据分量。零空间 \mathcal{N}(A) 包含所有被采样过程丢失的信息，网络需要学习如何从训练数据的统计规律中推断这些缺失信息。当测试数据的统计特性与训练数据不同时（如罕见病变、不同解剖区域），网络的推断可能产生错误的幻觉。

这一发现对临床应用具有警示意义：在追求更高重建精度的同时，必须仔细评估幻觉风险，特别是在涉及诊断决策的场景中。

非笛卡尔MRI的展开神经网络

本节介绍专门针对非笛卡尔采样设计的展开网络NC-PDNet（Non-Cartesian Primal Dual Network）。非笛卡尔采样（如螺旋、径向、SPARKLING等）相比笛卡尔采样具有更好的采样效率和运动鲁棒性，但重建更具挑战性。相关工作包括密度补偿展开网络用于2D和3D非笛卡尔MRI重建，以及开源的MRI-NUFFT Python包使非笛卡尔MR成像更加便捷。

NC-PDNet网络架构

NC-PDNet的整体结构遵循展开优化的范式，由多个迭代单元（IU）串联组成：

每个迭代单元的内部结构针对非笛卡尔采样进行了专门设计。输入图像 x_0 首先通过前向算子 \mathcal{A} 映射到非笛卡尔k空间，然后计算与测量数据 y 的残差 (\cdot - y)。关键的改进是引入了密度补偿算子 \text{DCp}（Density Compensation），用于校正非笛卡尔采样的非均匀密度分布。校正后的残差通过伴随算子 \mathcal{A}^H 映射回图像空间。

密度补偿的作用在于：非笛卡尔轨迹在k空间中心通常采样更密集，外围更稀疏。如果不进行补偿，直接应用伴随算子会导致低频成分被过度加权。密度补偿通过对每个k空间样本乘以与其Voronoi区域面积成正比的权重来校正这种不均匀性。

数据一致性校正项与原始图像估计 x_0 拼接后，送入CNN进行图像域处理，输出更新后的图像 x。此外，网络还包含灵敏度图精炼模块（Smaps refiner），使用U-Net优化初始估计的线圈灵敏度图。

方法比较与定量结果

在膝关节MRI数据上比较不同方法的重建效果。Adj.+DCp是简单的伴随加密度补偿方法，作为基线。DIP（Deep Image Prior）是一种无需训练数据的方法。grid.PDNet是在网格化数据上应用PDNet。PDNet和U-net是笛卡尔采样设计的方法直接应用于非笛卡尔数据。NCPDNet是专门设计的非笛卡尔方法。

从重建图像可以看到，NCPDNet在细节恢复方面表现最好，红框标记的区域显示其他方法存在模糊或伪影，而NCPDNet保持了清晰的结构。定量结果显示，在径向和螺旋采样下，NCPDNet分别达到40.00/0.9191和40.68/0.9258的PSNR/SSIM，显著优于其他方法。箱线图进一步确认了NCPDNet在PSNR和SSIM两个指标上的优势。

多线圈3D数据集训练

Calgary-Campinas数据集

Calgary-Campinas公开脑部MR数据集提供了12通道线圈的原始数据，具有以下特点：原始完全采样的复数k空间数据；T1加权MPRAGE序列；1mm各向同性分辨率；每个通道的采集矩阵大小为 N_x \times N_y \times N_z = 256 \times 218 \times [170, 180]。

数据集划分为：训练集47个3D多线圈体积，验证集20个体积，测试集50个体积。右侧展示了12个通道的图像，每个通道显示不同的灵敏度加权模式。最下方是RSS组合后的最终图像，综合了所有通道的信息。

可扩展且内存高效的3D多线圈NC-PDNet

技术方案

为了实现3D多线圈非笛卡尔MRI的高效重建，采用了三项关键技术：MRI-NUFFT提供快速的非均匀傅里叶变换PyTorch实现；线圈压缩（Coil Compression, CC）减少线圈数量以降低计算和内存需求；线圈无关训练（Coil-agnostic Training）使模型可以处理不同线圈配置的数据。

实验使用6倍回顾性欠采样，测试了四种3D非笛卡尔轨迹：3D Cones（锥形）、3D Radial（径向）、TPI（Twisted Projection Imaging）和GoLF-SPARKLING。

不同轨迹的定量比较

在12通道和32通道测试数据上评估不同轨迹的重建性能，所有模型均采用线圈无关配置训练。

12通道数据结果：TPI达到33.03 dB PSNR和0.932 SSIM；3D Radial为33.65/0.931；3D Cones为33.82/0.934；GoLF-SPARKLING配合线圈压缩达到36.83/0.960；GoLF-SPARKLING不使用线圈压缩达到最高的37.62/0.962。

32通道数据结果：性能普遍提升，GoLF-SPARKLING达到43.33/0.987的最佳性能。

箱线图显示，无论是12通道还是32通道数据，GoLF-SPARKLING轨迹配合NC-PDNet都获得最高的PSNR和SSIM，证明了优化轨迹设计与深度学习重建结合的优势。

重建效率与图像质量

在12通道数据上比较不同轨迹的重建结果。Ground truth显示清晰的矢状位脑部图像。3D Radial、3D Cones和TPI的重建结果存在不同程度的模糊或伪影。GoLF-SPARKLING的重建质量最高，无论是否使用线圈压缩都能获得清晰的图像。

使用线圈压缩时，GoLF-SPARKLING达到PSNR 36.867/SSIM 0.964；不使用线圈压缩时达到PSNR 37.381/SSIM 0.964。两者SSIM相同，PSNR仅相差约0.5 dB，说明线圈压缩对图像质量的影响很小。

32通道体积推理的计算效率比较：不使用线圈压缩时，重建时间为21.1秒，GPU内存使用19.74 GB；使用线圈压缩后，重建时间降至4.95秒，GPU内存使用仅5.49 GB。线圈压缩带来约4倍的速度提升和约3.6倍的内存节省，同时保持相近的图像质量。

从回顾性到前瞻性验证

方法优势

端到端展开深度神经网络在3D多线圈非笛卡尔欠采样数据上的训练取得了显著成果。GoLF-SPARKLING配合NC-PDNet在32通道数据重建中获得+2.43 dB PSNR和+0.01 SSIM的提升。内存和计算效率方面，1mm各向同性体积可以在4.95秒内完成重建，仅使用5.49 GB GPU内存。可扩展方案结合了线圈无关训练、线圈压缩和快速NUFFT PyTorch实现，在性能损失极小的情况下大幅提升了效率。

方法局限

主要局限在于对前瞻性欠采样数据的泛化能力有限。回顾性加速数据的噪声水平低于前瞻性加速数据，因为回顾性方法是从完全采样数据中模拟欠采样，而前瞻性方法直接采集更少的数据，导致信噪比降低。此外，3T NeuroSpin数据与Calgary数据集之间可能存在分布偏移，包括扫描仪差异、序列参数差异等。

从实际测试可以看到，左侧Ground Truth来自笛卡尔完全采样数据的IFFT重建，右侧是NC-PDNet对3T NeuroSpin前瞻性欠采样数据的重建结果。虽然整体结构保持，但细节和对比度存在一定差异，反映了训练数据与实际应用数据之间的域偏移影响。

3T内部原始数据集预处理

数据集描述

NeuroSpin内部数据集包含79例健康被试的3D T1加权MPRAGE脑部扫描，采集参数为1mm各向同性分辨率，使用3T MRI扫描仪和20通道接收线圈。完全采样的笛卡尔采集扫描时间约9分钟，采集矩阵大小为 256 \times 240 \times 176。

预处理步骤

数据预处理包含以下关键步骤。首先使用GoLF-SPARKLING轨迹进行欠采样，加速倍数约为5倍（AF≈5）。然后进行噪声注入，向欠采样k空间数据添加高斯噪声，模拟前瞻性采集的噪声特性。接着进行SVD线圈压缩，将20通道数据压缩至7通道，减少计算和存储需求。最后进行k空间归一化，使用所有线圈中心k空间点的平均幅值对数据进行归一化。

Ground Truth定义为多线圈笛卡尔重建的RSS组合图像，作为训练和评估的参考标准。

归一化的重要性

直方图对比显示了k空间归一化的效果。

左图（无归一化）显示三种数据的幅值分布差异很大：笛卡尔参考（蓝色）、带噪声回顾性数据（橙色）和前瞻性数据（绿色）的分布峰值和范围明显不同。右图（有归一化）显示归一化后，回顾性数据和前瞻性数据的分布更加接近，这有助于减少训练数据与测试数据之间的域偏移。

前瞻性验证

实验设置与结果

前瞻性欠采样验证使用GoLF-SPARKLING轨迹（AF=5），与训练数据集中用于回顾性欠采样的轨迹相同。这种设置使得训练和测试的采样模式一致，减少了由轨迹差异导致的域偏移。

(a) Ground Truth来自完全采样笛卡尔采集，扫描时间9.2分钟。图像显示三个正交视图（轴位、矢状位、冠状位）以及放大的细节区域，展示清晰的脑部解剖结构。

(b) NC-PDNet重建结果来自GoLF-SPARKLING前瞻性欠采样数据，扫描时间仅1.8分钟，重建时间5秒。重建图像与Ground Truth整体一致，细节区域（红框）的放大图显示网络能够恢复大部分解剖细节。

该图展示了GoLF-SPARKLING轨迹的3D可视化，显示其在k空间中的复杂螺旋结构，实现了高效的非笛卡尔采样覆盖。

扫描时间从9.2分钟缩短至1.8分钟，实现了约5倍加速，同时保持可接受的图像质量。这验证了NC-PDNet在前瞻性非笛卡尔MRI重建中的实用价值。

减少训练时间的并行策略

数据分布式并行

为加速NC-PDNet的训练，实现了数据分布式并行（Data Distributed Parallelism, DDP），使用cufinufft和gpunufft后端进行GPU加速的非均匀傅里叶变换。初始版本仅支持单对比度训练。使用4块GPU配合gpunufft后端，每个epoch的训练时间从2小时05分钟减少到28分钟，实现约4.5倍加速。

DDP的工作原理是：数据加载器（Dataloader）将不同的数据批次分发到多块GPU（GPU 0-3），每块GPU独立进行前向传播和反向传播。计算完成后，各GPU的梯度进行同步（Synchronize gradients），然后每块GPU独立更新各自的模型副本。这种方式保证了所有GPU上的模型参数保持一致。

完全分片数据并行

为支持双对比度训练（同时训练多种MRI对比度），引入了完全分片数据并行（Fully Sharded Data Parallelism, FSDP）。FSDP的核心思想是在多GPU之间分片存储模型参数、梯度和优化器状态，从而提高内存效率。双对比度训练的每个epoch时间从2小时11分钟减少到36分钟。

零冗余优化器（Zero Redundancy Optimizer, ZeRO）是FSDP的理论基础，通过将模型状态分布存储到多GPU来减少内存占用。传统方法中每块GPU存储完整的模型参数、梯度和优化器状态，而ZeRO Stage 3将这三者都进行分片，每块GPU只存储其中一部分。

FSDP的具体流程是：每块GPU存储模型参数的一个分片；前向传播时，通过Get weights操作收集其他GPU的参数分片；完成前向传播后，进行反向传播计算梯度；梯度同步后各GPU独立更新其负责的参数分片。

当前方法仍然受限于内存约束，只能使用较小的图像精炼网络架构。

流水线并行

为进一步克服内存限制，正在实现流水线并行（Pipeline Parallelism）。

模型并行（Model Parallelism）将模型的不同层放置在不同GPU上，但存在GPU利用率低的问题：当一层在进行计算时，其他层所在的GPU处于空闲状态。

流水线并行通过将输入分割成多个微批次（microbatch），在多GPU上流水线式执行来提高利用率。图示显示：纵轴为4块GPU分别放置模型的4层（F_0 到 F_3 的前向传播，B_0 到 B_3 的反向传播）；横轴为时间。流水线并行使得当GPU 0处理第一个微批次的 F_{0,1} 时，GPU 1可以同时处理第零个微批次的 F_{1,0}，以此类推。但流水线中仍存在气泡（Bubble），即某些时刻某些GPU处于空闲状态，这是流水线启动和结束阶段不可避免的开销。该功能目前仍在调试中。

堆叠ΔB₀-PDNet架构

离共振校正的信号模型

在非笛卡尔MRI中，B_0 场不均匀性会导致离共振效应，产生图像模糊和几何畸变。信号模型可以表示为：

其中 s(t) 是在时刻 t 采集的信号，f(r) 是位置 r 处的图像值，\omega(r) 是位置 r 处的离共振频率（单位为rad/s），k(t_m) 是时刻 t_m 的k空间位置。离共振项 e^{-i\omega_n t_m} 导致信号相位随时间累积偏移。

为了高效计算，将离共振相位项近似为有限个基函数的线性组合：

其中 L \ll L_{\max}，b_{m,l} 是时间相关的系数，c_{l,n} 是空间相关的系数。代入信号模型得到：

这种分解将原本需要对每个体素单独处理的计算转化为 L 次标准NUFFT操作的加权和，大幅降低了计算复杂度。

网络架构

堆叠ΔB₀-PDNet采用展开结构，进行 N_i 次迭代。每次迭代包含两个主要模块。

第一个模块是数据一致性层。输入图像通过伴随Pseudo-NUFFT（考虑离共振的非均匀傅里叶变换的伴随算子）映射到k空间，经过神经网络处理后，再通过前向Pseudo-NUFFT进行数据一致性校正。

第二个模块是图像精炼网络。数据一致性校正后的图像再次通过伴随Pseudo-NUFFT映射到图像空间，经过神经网络进行去伪影和增强处理。

加速策略

加速该架构的关键在于减少NUFFT调用次数 N_F，具体策略包括：减少迭代次数 N_i，但可能影响重建质量；减少线圈通道数 N_c，通过线圈压缩实现；减少插值器数量 L，即使用更少的基函数近似离共振相位。这些策略需要在计算效率和重建质量之间权衡。

NC-PDNet中的ΔB₀校正效果

插值器数量的影响

实验比较了不同参数配置下的重建质量，关键参数包括迭代次数 N_i、线圈通道数 N_c 和离共振插值器数量 L。NUFFT调用次数 N_F 与这些参数直接相关，影响计算成本。

第一组对比展示了从无校正到完全校正的过渡。参考图像（无校正，N_i=20, N_c=20，N_F=780）显示存在明显的离共振伪影。使用网络重建时，L=1（N_F=45）的结果RMSE为0.033，PSNR为30.02 dB，SSIM为0.911，伪影仍然明显。L=3（N_F=135）的结果改善至RMSE 0.025，PSNR 32.21 dB，SSIM 0.944。L=5（N_F=225）达到最佳平衡，RMSE 0.022，PSNR 33.37 dB，SSIM 0.951，用绿色边框标记。参考图像（带校正，N_i=20, N_c=20, L=20，N_F=15600）显示完全校正的理想结果。

随着插值器数量 L 增加，离共振校正效果改善，但计算成本也相应增加。L=5 在质量和效率之间取得了良好平衡。

与传统方法的对比

第二组对比将网络方法与传统小波方法进行比较。参考图像（无校正）作为基线。小波方法（N_i=5, N_c=5, L=5，N_F=225）达到RMSE 0.027，PSNR 31.78 dB，SSIM 0.932。网络方法（无校正，N_i=5, N_c=5，N_F=45）的结果为RMSE 0.028，PSNR 31.26 dB，SSIM 0.926。网络方法（带校正，N_i=5, N_c=5, L=5，N_F=225）达到最佳结果RMSE 0.022，PSNR 33.37 dB，SSIM 0.951。

结果表明，带离共振校正的网络方法明显优于传统小波方法和无校正的网络方法，证明了将物理模型（离共振校正）集成到深度学习框架中的有效性。

实际扫描仪数据的离共振校正

高分辨率SWI重建

在实际扫描仪数据上验证方法的有效性。采集参数为0.6mm各向同性分辨率的SPARKLING轨迹，扫描时间3分钟，加速倍数AF=15。这是一个具有挑战性的场景，因为高分辨率和高加速倍数都会加剧离共振效应。

无离共振校正的结果：压缩感知（CS）重建耗时15分钟，展开网络重建耗时8分钟。两种方法的重建图像都存在明显的模糊和畸变（红色和蓝色箭头标记的区域），SSIM为0.8947。

带离共振校正的结果：CS重建需要8小时的极长计算时间，展开网络重建仅需10分钟。重建质量显著提升，SSIM达到0.9541。黄色和绿色箭头标记的区域显示细节恢复良好，之前的模糊和畸变得到有效校正。

展开网络方法在保持高重建质量的同时，将计算时间从8小时缩短到10分钟，实现了约48倍的加速。这对于临床应用具有重要意义。

通过数据欠采样的自监督学习

SSDU方法原理

传统监督学习需要完全采样的参考数据作为训练目标，但在许多实际场景中这类数据不可用。自监督学习通过数据欠采样（Self-Supervision via Data Undersampling, SSDU）提供了一种替代方案。

核心思想是将已采集的k空间位置 \Omega 分割成两个不相交的子集：

其中 \Theta 用于网络输入（训练集），\Lambda 用于损失计算（验证集）。网络仅看到 \Theta 中的数据，然后在 \Lambda 位置评估重建结果与实际测量的一致性。

训练流程

网络输入为在 \Theta 位置的k空间数据 y_\Theta^i 和对应的编码矩阵 E_\Theta^i。数据通过伴随算子 E_\Theta^H y_\Theta 映射到图像空间，结合灵敏度图进行多线圈处理。

展开网络包含多个单元（Unit 1, Unit 2, ..., Unit T），每个单元执行数据一致性（DC）和图像精炼操作。网络输出经过编码算子 E_\Lambda 映射到 \Lambda 位置的k空间，与实际测量 y_\Lambda^i 比较计算损失。

端到端最小化的目标函数为：

其中 f(\cdot; \theta) 是参数为 \theta 的展开网络，\mathcal{L} 是损失函数。通过反向传播更新网络参数。

这种方法的理论基础与Noisier2Noise类似，通过变密度子采样实现有效的自监督学习。

监督学习与自监督学习对比

前瞻性欠采样场景

实验场景为前瞻性欠采样数据（基础加速R=2），没有完全采样的参考数据可用，因此无法进行传统的监督深度学习训练。自监督方法在这种情况下提供了可行的训练策略。

测试了不同加速倍数（R=2, 4, 6, 8）下的重建效果。上排显示重建的脑部图像，即使在高达8倍加速的情况下，自监督方法仍能产生合理的重建结果。下排显示对应的k空间采样模式，加速倍数越高，采样越稀疏。

绿色边框标记的区域显示细节保持良好，证明自监督学习在高加速率下也能实现成功的重建。

非笛卡尔展开网络总结

核心要点

密度补偿是非笛卡尔展开神经网络的关键。非笛卡尔采样的非均匀密度分布必须通过适当的补偿来校正，否则会导致重建偏差。

快速且内存高效的训练和推理可以通过MRI-NUFFT库配合线圈压缩实现。MRI-NUFFT提供了GPU加速的非均匀傅里叶变换，线圈压缩减少了需要处理的数据维度。

可扩展且通用的神经网络模型通过线圈无关训练实现。模型可以处理不同线圈配置的数据，无需针对每种配置重新训练。

重建性能显著依赖于k空间轨迹的选择。GoLF-SPARKLING等优化轨迹配合深度学习重建可以获得最佳性能。

训练集和测试集之间采样掩模的不相交分割是自监督学习的最优策略。这确保了网络必须学习真正的图像先验，而不是简单地记忆采样模式。

自监督学习的进一步发展

自监督学习使用随机方向k空间带的方法（K-band）是一个新兴方向。该方法通过在k空间子集上进行随机梯度下降实现自监督MRI重建，进一步扩展了无需完全采样参考数据的训练策略。

即插即用方法用于MRI重建

即插即用（Plug-and-Play, PnP）方法是一类将预训练去噪器集成到迭代优化框架中的图像重建方法。这类方法的核心思想是用学习得到的去噪神经网络替代传统优化算法中的近端算子，从而结合物理模型的数据一致性约束和深度学习的图像先验。相关研究涵盖了理论、算法和应用多个层面，包括可证明的预条件PnP方法和等变PnP图像重建等。

基本框架

PnP方法求解的逆问题形式为 y = Ax + n，对应的优化问题为：

其中 f(x) 是数据保真项，g(x) 是正则项。传统的前向-后向分裂（Forward-Backward Splitting）算法迭代格式为：

第一步沿数据保真项的负梯度方向移动，第二步应用正则项的近端算子。PnP方法的关键创新是用神经网络 G_\theta 替代近端算子：

图像对比显示了PnP-ISTA（40.58 dB）和带谱归一化的PnP-ISTA-SN（34.87 dB）的重建结果。谱归一化用于约束网络的Lipschitz常数以保证稳定性，但会牺牲一定的去噪性能。这体现了PnP方法中稳定性与精度之间的权衡：对去噪器施加Lipschitz约束可以保证算法收敛，但会降低去噪器 G_\theta 的去噪性能。

解决这一权衡的方案包括：半二次分裂与预条件方法、等变学习方法。

预条件即插即用方法

预条件数据一致性

预条件方法通过引入预条件矩阵 P 来重新加权数据一致性项：

其中 \|\cdot\|_{2,P}^2 表示带预条件的加权范数。对应的前向-后向分裂变为：

PnP版本将近端算子替换为神经网络：

半二次分裂方法

半二次分裂（Half-Quadratic Splitting, HQS）是另一种策略，其迭代格式为：

这里 \text{prox}_{\gamma_k f}^P 是带预条件的数据保真项近端算子，G_{\theta_k} 是去噪网络。

预条件矩阵设计

预条件矩阵 P 的选择目标是使数据一致性步骤更稳定。优化问题为：

其中 \rho(\cdot) 表示谱半径。一个有效的选择是F-1预条件：

适当选择参数 \alpha 可以使谱半径最小化，从而稳定PnP方法的迭代。

去噪器的无监督训练

Noise2Noise方法

传统的去噪器训练需要成对的干净-噪声图像，但在许多实际场景中干净图像不可用。Noise2Noise方法的关键观察是：训练网络 f_\theta 最小化两幅独立噪声图像之间的差异：

同样可以得到去噪神经网络。这是因为当 y_1 和 y_2 是同一干净图像加独立噪声时，最小化上述目标等价于学习去除噪声恢复底层干净图像。

Neighbor2Neighbor方法

Neighbor2Neighbor方法更进一步，仅需单幅噪声图像即可训练去噪器。方法通过采样器 G 从原始图像 y 中提取两个子采样版本 g_1y 和 g_2y。如果 g_1 和 g_2 选取的是相邻但不重叠的像素，则这两个子采样图像的噪声近似独立。

训练目标为：

第一项是基本的Noise2Noise损失，第二项是一致性正则项，确保去噪器对不同子采样的输出保持一致。

从示例图像可以看到，输入的噪声人脸图像经过训练后的网络处理，输出清晰的去噪结果。MRI应用中，输入的带噪脑部图像同样可以被有效去噪。

2D解剖MRI的回顾性实验

方法比较

在2D解剖MRI数据上比较不同重建方法。测试了两个加速倍数：AF=4和AF=16。方法包括FISTA-Wavelet（传统压缩感知）、HQS-F1（带F-1预条件的半二次分裂PnP）、PNP-Id（标准PnP）、NCPDNET（展开网络）以及Ground Truth参考。

AF=4结果：FISTA-Wavelet达到33.234 dB/0.923 SSIM；HQS-F1达到最佳的36.000 dB/0.921 SSIM；PNP-Id为34.025 dB/0.891 SSIM；NCPDNET为37.765 dB/0.977 SSIM。

AF=16结果：FISTA-Wavelet降至29.771 dB/0.867 SSIM；HQS-F1保持较好的33.356 dB/0.909 SSIM；PNP-Id为32.113 dB/0.874 SSIM；NCPDNET降至29.912 dB/0.916 SSIM。

收敛曲线显示，HQS配合F1预条件（绿色实线）在200次迭代内达到最高PSNR且收敛稳定。其他预条件方案（Id、Cheb、Dyn等）的收敛性能各有差异。

核心发现

HQS-F1方法在各种条件下表现最佳，特别是在高加速倍数下保持稳定。NCPDNET在分布外数据（如高加速倍数）上出现性能崩溃，这与展开网络对训练数据分布敏感的特性一致。PnP方法对加速倍数变化具有鲁棒性，无需针对不同加速倍数进行微调。

前瞻性3D MPRAGE重建

多种方法对比

在前瞻性8倍加速的3D MPRAGE数据上比较不同采样和重建策略。采样方案包括：GRAPPA 2x1（传统并行成像，扫描5分钟）、Poisson Disk AF=4（随机欠采样，扫描2分钟）、以及SPARKLING GoLF配合GRAPPA 2x2（优化非笛卡尔轨迹，扫描1分钟）。

重建方法包括：CS（压缩感知）、DC Adjoint（即时伴随重建）、Plug-n-Play（即插即用）、NC-PDNet（非笛卡尔展开网络）。

结果显示：传统GRAPPA需要5分钟扫描获得基线质量；Poisson Disk采样配合CS或DC Adjoint在2分钟或1分钟内完成，但图像质量有限；SPARKLING GoLF配合GRAPPA 2x2仅需1分钟扫描，使用Plug-n-Play方法（用绿色边框标记）获得最佳图像质量，该方法将3D DRUNet嵌入重建流程；NC-PDNet作为端到端训练的展开网络同样在1分钟扫描数据上取得良好结果。

Plug-n-Play方法的优势在于：使用预训练的3D DRUNet作为去噪器，无需针对特定采样模式重新训练，对不同加速倍数和轨迹设计具有更好的泛化能力。

SNAKE fMRI模拟器的研究背景

模拟器框架

SNAKE（Simulator from Neuro-Activation to K-space Evaluation）是一个模块化的fMRI数据模拟器，可以从时空域到k空间进行端到端的模拟。模拟器的输入包括三个主要组件：采集设置（Acquisition Setting）定义了各种k空间轨迹，包括笛卡尔、径向、EPI和3D SPARKLING等；大脑模型（Brain, B_0, Activated Regions）提供解剖结构、场图和激活区域定义；实验范式设计和血流动力学响应函数（HRF）模型定义了刺激时序和神经-血管耦合特性。

这些输入通过SNAKE模拟器处理，生成模拟的k空间数据。然后数据经过重建算法处理，最后进行统计分析验证，评估不同采集和重建策略对激活检测的影响。这个框架使得研究者可以在没有真实扫描的情况下系统地评估各种方法学选择。

加速功能MRI的即插即用方法

完整处理流程

fMRI数据处理包含四个主要阶段：采集阶段获取随时间变化的k空间数据；重建阶段从k空间恢复BOLD信号时间序列；统计分析阶段计算z-score激活图；神经科学解释阶段进行功能定位和认知推断。

实验配置

去噪器网络训练使用Calgary数据集上训练的3D DRUNet，输入为T1加权1mm分辨率图像，采用单线圈配置（虚拟线圈组合）。采集场景使用SNAKE模拟器生成，参数为3D 3mm各向同性分辨率，TR=0.75秒，AF=4加速，32通道线圈，模拟了Petrov等人2017年的实验设置。

定量结果

比较了三种方法：CG（共轭梯度）、CS（压缩感知）和HQS-F1（带F-1预条件的半二次分裂PnP）。HQS-F1在各项指标上表现最佳：PSNR达到19.338 dB，SSIM为0.577，tSNR（时间信噪比）为6.286，AUC（曲线下面积）为0.574，BACC（平衡准确率）为0.891。

核心发现是首次将SNAKE与3D PnP方法结合用于fMRI重建。HQS-F1方法在图像质量方面表现优异，统计分析性能略有下降，但去噪器是在分布外数据（T1加权，1mm）上训练的，说明方法具有一定的跨域泛化能力。

通过等变性实现PnP稳定性

等变性定义

如果算子 D 对于群 \mathcal{G} 的作用 T_g 满足以下条件，则称 D 关于 \mathcal{G} 是等变的：

这意味着先对输入应用变换 T_g 再去噪，与先去噪再对输出应用相同变换的结果相同。例如，如果 T_g 是旋转变换，等变去噪器对旋转后的图像去噪的结果等于对原图去噪后再旋转。

等变化构造

任何算子 D 都可以通过以下方式等变化：

这个构造对输入应用群中所有变换，分别去噪，然后将结果逆变换后平均。右侧示意图展示了这个过程：原始蝴蝶图像经过去噪器 D 处理；同时图像经过变换 T_g（如旋转）后去噪，然后逆变换 T_g^{-1}。等变化的结果是这些处理的平均。

稳定性改进

理论结果表明，等变去噪器可以在不施加训练时Lipschitz约束的情况下改善PnP算法的稳定性。收敛曲线显示，标准PnP（蓝色）在迭代过程中PSNR先上升后下降，出现不稳定；等变PnP（橙色）保持稳定收敛，最终PSNR更高。

从重建图像可以看到：反投影结果存在严重伪影；标准去噪器重建有明显残留伪影；等变去噪器重建质量最好，细节清晰且伪影最少。等变性通过对多个变换的平均有效地抑制了去噪器的不稳定行为。

PnP方法总结

核心优势

PnP方法是对前向算子不敏感的重建方法。这意味着方法对采集设置的变化（线圈灵敏度、采样模式、加速倍数）不太敏感，因为去噪器与物理模型是分离训练的。

稳定重建且不降低图像质量可以通过以下技术实现：半二次分裂配合预条件确保数据一致性步骤的稳定收敛；等变学习通过对称性约束改善去噪器的泛化能力。

去噪器可以采用监督或无监督方式训练，适用于2D和3D MRI（线圈组合后）。无监督方法如Noise2Noise和Neighbor2Neighbor在没有干净参考数据时特别有用。

相对于展开网络的优势

PnP方法在以下情况下优于展开网络：训练和测试之间欠采样倍数变化时，NC-PDNet等展开网络性能显著下降，而PnP方法保持稳定；从回顾性加速到前瞻性加速的噪声分布偏移时，PnP方法表现更鲁棒。

未来方向

训练和推理时的可扩展性仍是挑战。R2D2深度神经网络系列为可扩展的非笛卡尔磁共振成像提供了新的解决思路。

基于能量的MRI重建模型

本节介绍基于能量的模型（Energy-Based Models, EBM）在MRI重建中的应用。这类方法通过学习图像的能量函数来隐式定义图像先验，相关工作包括多尺度能量（MuSE）框架等。

能量函数与概率分布

基于能量的模型用神经网络 E_\theta(x): \mathbb{C}^m \to \mathbb{R} 建模图像的负对数先验（能量）。能量函数高的图像概率低，能量函数低的图像概率高。

对应的概率分布通过Boltzmann分布定义：

其中 Z_\theta = \int \exp(-E_\theta(x)) dx 是配分函数，确保概率归一化。左侧图示展示了能量函数 E_\theta(x)（上）与对应概率密度 p_\theta(x)（下）的关系：能量函数的谷对应概率密度的峰。

Fisher散度训练

训练能量模型的一种方法是最小化Fisher散度：

这个目标函数比较模型分布和真实分布的分数函数（对数概率密度的梯度）。分数函数 \nabla_x \log p_\theta(x) 指向概率密度增加最快的方向，与能量函数的负梯度相关。

去噪分数匹配训练能量模型

从Fisher散度到去噪分数匹配

直接计算Fisher散度需要知道真实分布 p(x) 的分数函数，这通常不可用。去噪分数匹配（Denoising Score Matching, DSM）提供了一种实用的替代方案。

考虑加噪版本 \tilde{x} = x + \sigma z，其中 z \sim \mathcal{N}(0, I)。DSM目标函数为：

由于 p_\sigma(\tilde{x}|x) = \mathcal{N}(\tilde{x}; x, \sigma^2 I)，其分数函数为 \nabla_{\tilde{x}} \log p_\sigma(\tilde{x}|x) = -(\tilde{x} - x)/\sigma^2 = -z/\sigma。

定义 H_\theta(x) = \nabla_x E_\theta(x) 为能量函数的梯度，则 \nabla_x \log p_\theta(x) = -H_\theta(x)/\sigma^2。代入后，DSM目标简化为：

训练过程解释

训练过程可以理解为：取干净图像 x，加入噪声得到 \tilde{x} = x + \sigma z；网络学习从噪声图像预测噪声 \sigma z；能量函数的梯度 H_\theta 实际上学习了去噪方向。

网络架构包括：输入 \tilde{x} 通过能量网络 E_\theta(x) 得到标量能量值；通过自动微分计算梯度 H_\theta(x) = \nabla_x E_\theta(x)；输出预测的噪声，与真实噪声比较计算损失。

单尺度能量的初始化敏感性

重建问题形式

基于能量的MRI重建问题为：

第一项是数据保真项，权重由 \eta^2 控制；第二项是能量先验项，权重由 \sigma^2 控制。能量函数定义为：

其中 \psi_\theta(x) 是去噪网络的输出。这个能量形式意味着能量越低，图像越接近去噪器认为的"干净"图像。

初始化影响

实验显示单尺度能量方法对初始化高度敏感。使用不同初始化的重建结果差异显著：Init-Atb（伴随初始化，A^H y）达到39.84 dB；Init-Org（真实图像初始化）达到最佳的41.16 dB；Init-Randn（随机初始化）仅有9.38 dB，完全失败。

误差图（第二行，红色椭圆标记区域）显示：伴随初始化的结果存在局部误差；真实图像初始化的误差最小；随机初始化的误差极大，覆盖整个图像。

这种初始化敏感性源于能量函数的非凸性。当初始点远离最优解时，优化可能陷入局部极小值或发散。这是单尺度能量模型的主要局限，需要通过多尺度方法来克服。

显式多尺度能量改进收敛性

尺度参数化的平滑分布

单尺度能量方法的初始化敏感性源于能量景观的复杂性。多尺度方法通过学习由噪声标准差 \sigma（尺度）参数化的平滑概率分布来解决这个问题。

基本思想是考虑加噪版本 x + \sigma z，其对应的概率分布是原始分布 p(x) 与高斯核 \mathcal{N}(0, \sigma^2) 的卷积。较大的 \sigma 产生更平滑的分布，较小的 \sigma 保留更多细节。

训练目标是学习尺度相关的能量函数：

网络参数 \theta 依赖于尺度 \sigma，不同尺度对应不同的能量函数。

从粗到精的能量演化

图示展示了不同尺度下能量函数等高线的变化。以二维优化问题为例，红点表示局部极小值位置。\sigma_1 = 0.5（最粗尺度）时，能量景观非常平滑，只有少数明显的吸引盆地，便于找到全局结构。\sigma_2 = 0.3 时，开始出现更多细节，但整体仍较平滑。\sigma_3 = 0.2 时，局部极小值变得更加明显，但吸引盆地仍然较宽。\sigma_4 = 0.1（最细尺度）时，能量景观呈现完整的复杂结构，包含所有细节。

重建过程从粗尺度开始，逐步过渡到细尺度。粗尺度帮助找到正确的吸引盆地，细尺度恢复精细细节。这种从粗到精的策略避免了单尺度方法容易陷入局部极小值的问题。

显式多尺度：多个去噪分数匹配

多尺度网络架构

显式多尺度方法为每个尺度训练独立的能量网络。

在尺度1（Scale 1，最粗），输入图像加上大噪声后，通过能量网络 E_1(x) 计算能量，然后求梯度 H_1(x) = \nabla_x E_1(x)，减去噪声估计得到初步去噪结果。

在尺度2（Scale 2），使用上一尺度的输出加上中等噪声，通过另一个能量网络 E_2(x) 处理。梯度 H_2(x) = \nabla_x E_2(x) 提供更精细的校正。

这个过程持续到尺度k（Scale k，最细），每个尺度都有专门的能量网络 E_k(x)。尺度顺序满足 Scale 1 < Scale 2 < ... < Scale k，即从粗到精逐步精化。

从示意图可以看到，随着尺度增加，输入噪声减小，输出图像逐渐清晰。每个尺度的网络专注于恢复特定频率范围的信息。

多尺度优化对初始化不敏感

鲁棒性验证

多尺度能量方法的重建问题形式为：

能量项现在是尺度相关的 E_{\theta(\sigma)}，定义为：

实验验证显示，多尺度方法对初始化高度鲁棒。三种初始化方案（Init-Atb伴随初始化、Init-Org真实图像初始化、Init-Randn随机初始化）都收敛到几乎相同的结果：分别达到41.51 dB、41.50 dB和41.51 dB。

误差图（第三行，红色椭圆标记区域）显示三种初始化的误差分布几乎一致，都呈现低且均匀的误差。这与单尺度方法形成鲜明对比，单尺度方法中随机初始化导致灾难性失败（仅9.38 dB）。

多尺度策略通过从平滑能量景观开始优化，确保所有初始点都能找到正确的吸引盆地，然后逐步精化到最终解。

隐式多尺度：单一能量适用所有尺度

统一网络架构

显式多尺度方法需要为每个尺度训练独立网络，参数量和训练成本较高。隐式多尺度方法用单一网络处理所有尺度。

训练目标修改为在多个尺度上的期望：

外层期望对尺度 \sigma 取平均，使得单一网络 I_\theta 学习在所有尺度上有效的能量表示。

网络架构显示，同一个网络 I_\theta 配合梯度计算 H_\theta 应用于不同尺度的输入。对于小噪声 \sigma_1、中等噪声 \sigma_2 到大噪声 \sigma_K，网络输出对应尺度的去噪结果。这种设计大大减少了参数量，同时保持了多尺度的优势。

隐式能量的σ二次特性

理论分析

隐式能量函数对 \sigma 呈二次关系。

左图显示能量随噪声标准差的变化：隐式方法（蓝色）呈现光滑的二次曲线；显式方法在不同 \sigma 值（0.09和0.03）下是独立的曲线（黄色和红色）。

右图显示分数函数（负对数概率密度梯度）与噪声 z 内积的关系。隐式方法在零点附近呈线性（见插图放大），符合理论预期。

隐式能量的数学形式为：

其中 d(\tilde{x}, \mathcal{M}) 是加噪图像到干净图像流形 \mathcal{M} 的距离。这意味着隐式能量本质上测量的是图像偏离"干净"流形的程度。

关键等式是：

隐式能量的梯度等价于显式尺度相关能量的梯度，这解释了为什么单一隐式网络可以替代多个显式网络。

隐式MuSE接近展开神经网络性能

方法比较

比较隐式MuSE（i-MuSE）、显式MuSE（e-MuSE）、单尺度能量模型（EM）和端到端展开网络（E2E-MoDL）。重建问题统一为：

定量结果显示：i-MuSE达到40.07 dB；e-MuSE为38.99 dB；EM（\sigma=0.01）为38.91 dB；E2E-MoDL达到最高的40.88 dB。

从重建图像（第一行和第二行，红框标记的细节区域）可以看到，i-MuSE的图像质量接近E2E-MoDL，优于e-MuSE和单尺度EM。误差图（第三行）显示i-MuSE的误差分布与E2E-MoDL相近，整体误差低且均匀。

值得注意的是，i-MuSE作为一种基于能量的方法，不需要端到端训练，却能达到接近E2E-MoDL的性能。这表明隐式多尺度能量模型是一种有效的图像先验表示方法，结合了基于能量方法的灵活性和深度学习的重建质量。

隐式MuSE对前向算子失配不敏感但无法去除结构化伪影

前向算子失配实验

实验比较了i-MuSE和E2E-MoDL在前向算子失配情况下的表现。

当测试时的采样模式与训练时不同，i-MuSE达到35.01 dB，而E2E-MoDL仅有30.43 dB。

从重建图像可以看到，i-MuSE在各个视图（轴位、矢状位、冠状位）都能保持较好的图像质量，红框标记的细节区域显示解剖结构清晰。E2E-MoDL则出现明显的伪影和失真，特别是在冠状位图像中。k空间图（最下排）显示i-MuSE的频谱特性更接近参考，而E2E-MoDL出现异常模式。

这说明基于能量的方法对前向算子（采样模式、线圈配置等）的变化更加鲁棒，因为能量先验是独立于物理模型训练的。

结构化扰动的挑战

虽然i-MuSE对前向算子失配鲁棒，但实验揭示了其对结构化伪影的局限性。右侧展示了扰动分析：原始图像 x（Org. Img）是干净的脑部MRI；高斯扰动 z（Gaussian pert.）是随机噪声模式；结构化扰动 s（Structural pert.）具有与图像相关的空间结构。

组合扰动图像为 \dot{x} = x + \alpha_z^* z + \alpha_s^* s，其中 \alpha_z^* 和 \alpha_s^* 是扰动系数。能量函数 E_\theta^{\text{MuSE}}(x + \alpha_z z + \alpha_s s) 的等高线图显示：在 (\alpha_z^*, \alpha_s^*) 位置（虚线交叉点）的能量值约为0.36，表明能量模型难以区分结构化扰动和真实图像结构。

差异图 |x - \dot{x}| 显示扰动主要分布在边缘和高对比度区域，这正是结构化伪影（如混叠、截断伪影）通常出现的位置。这种结构化扰动与图像本身的统计特性相似，使得基于高斯噪声假设训练的去噪器难以有效去除。

基于能量模型的总结

核心要点

基于能量的方法学习先验概率分布（或其分数函数）。去噪分数匹配（DSM）提供了一种无需计算配分函数的实用训练方法。多尺度策略比单尺度更稳定，通过从粗到精的优化避免陷入局部极小值。隐式多尺度比显式多尺度更高效且自动化，用单一网络替代多个尺度特定网络。

能量模型与PnP方法的关系密切，都是通用的、对前向算子不敏感的方法。两者都将图像先验与物理模型分离，使得同一先验可以应用于不同的成像配置。

未来方向

结构化扰动的去噪与采样是一个重要研究方向。DEEPEN（Memory-efficient deep end-to-end posterior network）方法针对逆问题中的结构化扰动提出了解决方案。

收敛性保证也是活跃的研究领域。局部凸多尺度能量模型（LC-MUSE）为MAP图像恢复提供了理论收敛保证。

第六讲结论：方法对比

综合评估表格

四类方法（CS压缩感知、PnP即插即用、E2E MoDL端到端展开网络、能量模型）在七个维度上的表现如下。

基于能量（Energy based）：CS方法本质上基于能量（稀疏性正则化可视为能量项）；PnP方法不是显式基于能量；E2E MoDL不基于能量；能量模型显然基于能量。

收敛保证（Conv. guarantees）：CS方法在凸问题下有保证；PnP方法在CNN满足压缩条件时有保证，但这会损害图像质量；E2E MoDL无收敛保证；能量模型有收敛保证。

通用性（Universal，适用于任意前向模型）：CS方法通用；PnP方法通用（如果CNN是压缩的）；E2E MoDL不通用，需要针对特定配置重新训练；能量模型通用。

图像质量（Image quality）：CS方法质量一般；PnP方法质量好；E2E MoDL质量最好但仅在训练分布内；能量模型质量好。

训练内存需求（Memory demand during training）：CS方法不适用（无需训练）；PnP方法低；E2E MoDL高（需要反向传播通过整个展开网络）；能量模型低。

结果对初始化的依赖（Results dependent on initialization）：CS方法有依赖；PnP方法无依赖；E2E MoDL不适用；能量模型（多尺度）无依赖。

支持采样（Enables sampling）：CS方法不支持；PnP方法不支持；E2E MoDL不支持；能量模型支持（可以从学习的分布中采样）。

能量模型在所有维度上都表现良好，是唯一一类在所有评估标准上都获得肯定评价的方法。

本教程未涵盖的主题

生成模型

变分自编码器（Variational Autoencoder, VAE）是一类重要的生成模型。VAE包含编码器 \text{Encoder}(\phi) 将输入图像 x_0 = A^H y 编码为潜在分布的参数（均值 \mu 和标准差 \sigma），从中采样得到潜在变量 Z；解码器 \text{Decoder}(\theta) 从潜在空间重建图像。VAE可以学习图像的低维表示，用于MRI重建中的先验建模。

生成对抗网络（Generative Adversarial Network, GAN）由生成器 G（参数 \theta_G）和判别器 D（参数 \theta_D）组成。生成器从伴随重建 x_0 = A^H y 生成重建图像 x_{\text{Ref}}；判别器区分生成图像 x_{\text{Gen}} 和真实图像。通过对抗训练 \nabla_{\theta_G} G 和 \nabla_{\theta_D} D，生成器学习产生逼真的重建结果。GAN可以生成视觉上非常逼真的图像，但可能产生不存在于真实数据中的幻觉细节。

扩散模型

扩散模型（Diffusion Model）是近年来非常成功的生成模型。前向过程逐步向图像添加噪声，从干净图像（Step 0）到纯噪声（Step N）。反向过程学习逐步去噪，从 x_0 = A^H y 开始，通过多步去噪恢复清晰图像。每一步都结合测量数据 y 进行数据一致性约束。

扩散模型的优势在于生成质量高、训练稳定，但推理速度较慢（需要多步迭代）。在MRI重建中，扩散模型可以作为强大的图像先验，特别是在高加速倍数下的重建任务中。

归一化流模型

归一化流模型（Normalizing Flow-based Model）通过可逆变换建立复杂分布与简单分布之间的映射。从复杂的图像分布 p(x) 通过逆变换 f^{-1}(x; \theta) 映射到简单的潜在分布 p_Z(z)（通常是高斯分布）；通过正向流 f(z; \theta) 从潜在空间生成图像。

归一化流的优势是可以精确计算似然函数，便于训练和评估。在MRI重建中，可以用于建模图像的精确后验分布，支持不确定性量化。

联合学习采集轨迹与重建网络

MRI成像的两大核心问题

在前面几讲中，我们分别讨论了k空间采样轨迹的优化设计以及基于深度学习的图像重建方法。这一部分的核心思想是将这两个原本独立处理的问题融合在一起，通过端到端的方式同时优化采样轨迹和重建网络，从而实现整体性能的最大化。

从图像质量的演进可以清晰看到这种融合的价值：传统径向采样得到的脑部图像存在明显伪影；使用SPARKLING优化轨迹配合压缩感知重建后图像质量有所提升；进一步使用SPARKLING配合学习的重建网络效果更好；而当我们同时学习轨迹和重建网络时，可以获得最佳的图像质量。这种渐进式的改进说明，采集和重建两个环节的协同优化能够带来超越单独优化的效果。

MRI成像本质上包含两个相互关联的子问题。第一个是采集问题：在硬件约束条件下如何高效地对k空间进行采样。SPARKLING方法通过优化问题

来设计满足硬件约束集合 Q 的轨迹 K，使其尽可能逼近目标采样密度分布 \pi。这里 C(K, \pi) 衡量实际轨迹采样分布与目标分布之间的差异。

第二个是重建问题：如何从欠采样的k空间数据高效地恢复高质量图像。传统的非线性迭代重建方法虽然效果不错，但重建时间通常需要15到20分钟。而基于深度学习的方法可以将重建时间缩短到1分钟以内。深度网络的训练通过最小化重建损失来完成：

其中 \mathcal{R}_\theta 是参数为 \theta 的重建网络。然而，这种分离式的方法存在一个根本问题：轨迹设计时并不知道后续会用什么重建算法，而重建网络训练时使用的是固定的采样模式。这种信息的不对称限制了整体性能的提升。

联合学习的通用框架

为了解决上述问题，需要建立一个能够同时优化轨迹和重建网络的统一框架。

整个流程从训练数据开始，输入是复数MR图像 x \in \mathbb{C}^{N \times N}。首先通过虚拟线圈组合将多线圈数据合并。然后图像经过采集模型 F_{K_{t+1}} 的前向变换，该模型使用当前迭代步的轨迹 K_{t+1} 进行非均匀傅里叶变换，生成模拟的k空间数据 y \in \mathbb{C}^{N_c \times N_s \times \frac{\Delta t}{\delta t}}。这里 N_c 是线圈数，N_s 是采样点数，\frac{\Delta t}{\delta t} 表示时间维度的离散化。

k空间数据随后被送入重建模型 \mathcal{R}_\theta^{K_{t+1}}，该模型的参数 \theta 和所使用的轨迹信息 K_{t+1} 都会影响重建结果。重建得到的图像 \hat{x} 与原始图像 x 进行比较，计算损失函数 \mathcal{L}。

在轨迹优化路径上，每次迭代得到的中间轨迹 K_{t+\frac{1}{2}} 需要经过硬件约束投影算子 \Pi_\Omega 的处理，确保轨迹满足梯度幅值、摆率等物理限制，得到可行的轨迹 K_{t+1}。同时，轨迹的变化会影响密度补偿因子 D_{K_{t+1}}，这些补偿因子需要根据新轨迹重新计算，以校正非笛卡尔采样的非均匀密度分布。

整个框架的核心优化问题可以写成：

这个公式的含义是同时寻找最优轨迹 \hat{K} 和最优网络参数 \hat{\theta}，使得重建损失最小化。这里 F_{S(K)} 表示依赖于轨迹 K 的前向采集算子，S(K) 将轨迹映射为对应的采样模式。重建网络 \mathcal{R}_K^\theta 同时依赖于轨迹和网络参数，因为不同的轨迹可能需要不同的重建策略。损失函数 \mathcal{L}_r 通常采用L1范数、L2范数和多尺度结构相似度MSSIM的组合，以同时优化像素级精度和感知质量。

在重建网络的选择上，可以使用无参数的方法如密度补偿伴随重建，这种方法简单快速但精度有限；也可以使用学习的网络如U-Net或专门为非笛卡尔采样设计的NC-PDNet，这些网络能够获得更高的重建质量。

反向传播与梯度计算

联合优化框架的训练通过反向传播来实现，需要同时更新轨迹和重建网络参数。给定当前的轨迹 K 和网络参数 \theta，首先计算重建结果：

然后通过链式法则计算轨迹的梯度更新：

这个更新公式包含三个梯度分量的乘积。第一项 \nabla \mathcal{L}_r(x, \hat{x}) 是损失函数对重建图像的梯度，衡量当前重建与目标的差距。第二项 \frac{\partial \mathcal{R}_K^\theta}{\partial K} 是重建网络输出对轨迹的偏导数，反映轨迹变化如何通过重建网络影响最终结果。第三项 \frac{\partial (F_{S(K)}x)}{\partial K} 是k空间数据对轨迹的偏导数，这是关键的一环，因为采集过程直接依赖于轨迹的几何形状。

计算 \frac{\partial (F_{S(K)}x)}{\partial K} 需要对非均匀傅里叶变换求导。如果使用直接的非均匀离散傅里叶变换NDFT，计算复杂度为 O(p^2)，其中 p 是采样点数。这个二次复杂度在大规模问题中是不可接受的。实际应用中使用非均匀快速傅里叶变换NUFFT，复杂度降为 O(p \log p)，但NUFFT涉及插值操作，其自动微分得到的梯度可能不够精确。

PILOT方法的局限性

PILOT方法采用NUFFT进行前向计算并通过自动微分获取梯度。从实验结果来看，学习得到的轨迹在k空间中呈现出类似径向的模式，但梯度波形在时间域上显示出高频振荡。这种振荡的根本原因是NUFFT自动微分产生的梯度不够精确。当我们比较通过自动微分NUFFT得到的梯度与直接计算NDFT梯度时，可以发现两者存在明显差异：NUFFT梯度呈现出较大的波动，而NDFT梯度更加平滑。这种梯度误差会累积并影响轨迹优化的收敛性和最终质量。

BJORK方法的改进

BJORK方法针对上述问题提出了三个改进策略。首先是使用B样条参数化来表示轨迹，而不是直接优化每个采样点的位置。B样条提供了一种平滑的曲线表示方式，天然地限制了轨迹的高频变化，使得优化更加稳定。其次是将硬件约束作为惩罚项加入损失函数，而不是通过投影来强制满足。这种软约束的方式使得梯度流更加平滑，避免了投影操作带来的不连续性。第三是采用更精确的梯度计算方法，减少NUFFT近似带来的误差。

从BJORK的实验结果可以看到，学习得到的轨迹同样呈现合理的k空间覆盖模式，但梯度波形更加规整，没有PILOT中观察到的高频振荡。这说明B样条参数化和改进的梯度计算确实提升了优化的稳定性。

硬件约束的两种强制策略

在联合优化轨迹和重建网络时，必须确保学习得到的轨迹满足MRI扫描仪的物理限制。扫描仪只能沿着满足特定约束的曲线进行采样，这些约束主要包括梯度幅值上限和梯度摆率上限。硬件约束集合可以数学表示为：

这里 \dot{k} 是轨迹的一阶导数，对应于梯度波形，其无穷范数受 \alpha 限制，即峰值梯度强度 G_{max}；\ddot{k} 是二阶导数，对应于梯度的变化率即摆率，其无穷范数受 \beta 限制，即最大摆率 S_{max}。梯度波形曲线通常呈现先快速上升到峰值、维持一段时间、然后下降的模式。

惩罚方法

第一种策略是将硬件约束作为惩罚项添加到损失函数中。总损失变为重建损失与约束惩罚的加权和：

约束惩罚项 \mathcal{L}_c(K) 的具体形式为：

这个惩罚对所有 N_c 条轨迹的所有 N_s 个采样点进行累加。每个采样点的梯度幅值和摆率分别通过惩罚函数 \phi_{G_{max}} 和 \phi_{S_{max}} 进行约束，权重系数 \lambda_1 和 \lambda_2 控制两种约束的相对强度。惩罚函数采用铰链形式：

这个函数在 x \leq a 时为零，在 x > a 时线性增长。因此只有当梯度或摆率超过限制时才会产生惩罚。

轨迹的更新规则相应变为：

惩罚方法的优点是实现简单，梯度计算直接。但存在三个主要缺点：需要调节超参数 \lambda_1 和 \lambda_2，不同的值会导致不同的优化轨迹；即使训练收敛，也不能保证最终轨迹严格满足约束，可能存在轻微违反；惩罚项的梯度可能与重建损失的梯度方向冲突，导致优化过程被扭曲。

投影方法

第二种策略采用SPARKLING中使用的投影方法。原始SPARKLING的轨迹更新为：

这里先沿着目标采样密度匹配损失的负梯度方向移动，然后通过投影算子 \Pi_Q 将结果映射到可行集 Q 内。

将这种思想应用到联合学习框架，得到PROJeCTOR方法的更新规则：

PROJeCTOR的全称是PRojection for Jointly lEarning non-Cartesian Trajectories while Optimizing Reconstructor，即在优化重建器的同时联合学习非笛卡尔轨迹的投影方法。梯度下降步骤可能产生不满足约束的中间轨迹，但投影步骤 \Pi_{Q^{N_c}} 确保最终轨迹落在可行域内。投影算子作用于所有 N_c 条轨迹，将k空间中的任意轨迹映射为满足硬件约束的可行轨迹。

与现有方法的实验比较

BJORK和PILOT的对比

为了公平比较，所有方法使用相同的轨迹规格，并且PILOT和BJORK也采用NC-PDNet作为重建网络。实验在20倍加速因子下进行，测试了T1加权和T2加权两种对比度。

在T1加权图像上，PILOT方法达到SSIM 0.845和PSNR 28.874 dB；SPARKLING配合学习的采样密度达到SSIM 0.889和PSNR 29.869 dB；PROJeCTOR达到最高的SSIM 0.915和PSNR 32.614 dB。随着采样次数从16增加到32，对应欠采样因子从5.0降低到2.5，所有方法的性能都有提升，但PROJeCTOR始终保持领先。

在T2加权图像上观察到类似趋势：BJORK达到SSIM 0.874和PSNR 25.044 dB；SPARKLING为SSIM 0.903和PSNR 27.43 dB；PROJeCTOR为SSIM 0.934和PSNR 29.757 dB。箱线图显示PROJeCTOR在SSIM和PSNR两个指标上的分布都优于对比方法，中位数更高且方差更小。

轨迹形态与约束满足分析

通过可视化不同约束策略得到的轨迹，可以直观理解投影方法的优势。在3D k空间中，红色曲线代表不满足硬件约束的部分，蓝色曲线代表满足约束的部分。

情况(i)是既不使用惩罚也不使用投影的结果，轨迹几乎全部为红色，说明学习得到的轨迹严重违反硬件约束，在实际扫描仪上无法执行。情况(ii)使用较小的惩罚权重 \lambda = 0.0001，轨迹仍然大部分不可行。情况(iii)是PROJeCTOR的结果，轨迹完全为蓝色，严格满足所有硬件约束。情况(iv)使用较大的惩罚权重 \lambda = 0.001，轨迹基本可行但仍有少量违反。

从定量曲线来看，随着惩罚权重 \lambda 从 10^{-5} 增加到 10^{-1}，SSIM从约0.88下降到0.80，PSNR从约35 dB下降到30 dB。这说明增强约束惩罚会损害重建质量。与此同时，最大梯度 G_{max} 从约 4 \times 10^{-2} mT/m逐渐接近可行水平（红色虚线），最大摆率 S_{max} 从超过 10^0 T/m/s降低到接近可行水平。PROJeCTOR方法（黑色圆点）在保持最高图像质量的同时严格满足约束，实现了SSIM约0.88和PSNR约35 dB，比惩罚方法提升了约0.07 SSIM和3.69 dB PSNR。

回顾性研究的定量分析

在Calgary Brain数据集的160幅T1加权图像上进行了更大规模的回顾性研究，加速因子为20倍，重建时间控制在1分钟以内。

SPARKLING方法达到PSNR 31.03 dB和SSIM 0.92，从重建图像可以看到整体质量尚可，但细节区域存在一定模糊。PROJeCTOR方法达到PSNR 35.14 dB和SSIM 0.95，图像更加清晰锐利，解剖结构的边界更加分明。

残差图显示了重建图像与参考图像之间的差异。SPARKLING的残差图中可以观察到明显的结构化误差，特别是在脑沟和灰白质边界处。PROJeCTOR的残差图整体更加均匀，结构化误差明显减少，说明联合学习确实提升了重建的准确性。

箱线图的统计分析进一步确认了这一结论。在PSNR方面，SPARKLING的中位数为34.486 dB，PROJeCTOR为35.368 dB，两者差异的p值小于 1 \times 10^{-5}，统计显著。在SSIM方面，SPARKLING的中位数为0.935，PROJeCTOR为0.955，同样具有统计显著性（p \leq 1 \times 10^{-5}）。PROJeCTOR相比SPARKLING在PSNR上平均提升约0.9 dB，在SSIM上提升约0.02，这种提升在高度优化的基线上是相当可观的。

本讲总结

3D-SPARKLING方法从 B_0 场不均匀性校正中获益显著：校正后图像质量得到改善，时间信噪比tSNR增加，对BOLD效应的敏感性提高，BOLD相位图的准确性也得到提升。

头部运动仍然是一个需要解决的问题，相关研究正在探索运动校正的方法。一个可能的扩展方向是使用GoLF-SPARKLING结合动态估计的 \Delta B_0 图，实时补偿场不均匀性的变化。

当前正在进行的工作包括：基于自监督深度学习的快速图像重建方法，例如展开低秩加稀疏模型，利用SnakefMRI模拟器进行开发和验证；滑动窗口重建中诱导时间相关性的统计控制；以及在序贯学习范式中表征快速血流动力学响应的EXPLORE+方法。

迭代软阈值算法的特殊推导

L1正则化最小二乘问题

LASSO问题是压缩感知和稀疏信号恢复中最基础的优化模型。给定矩阵 A \in \mathbb{R}^{m \times n} 和观测向量 b \in \mathbb{R}^m，目标是求解 x \in \mathbb{R}^n 使得：

这个优化问题由两部分组成。第一项 \frac{1}{2}\|Ax - b\|_2^2 是数据拟合项，衡量线性模型 Ax 对观测数据 b 的逼近程度，这一项是可微的，可以直接使用梯度下降方法处理。第二项 \lambda\|x\|_1 是正则化项，其中 \lambda \geq 0 是正则化参数，控制稀疏性惩罚的强度。L1范数 \|x\|_1 在零点处不可微，但它是凸函数，能够诱导解的稀疏性，即倾向于让解向量中的许多分量精确等于零。

求解带L1正则化的最小二乘问题

首先回顾不带正则化的情况。当没有正则化项时，梯度下降的每一步可以理解为在当前点 x_k 处对目标函数 f 建立一个局部二次模型，然后求这个二次模型的最小值点作为下一个迭代点：

这里 f(x) = \frac{1}{2}\|Ax - b\|_2^2 是数据拟合项，t_k > 0 是第 k 步的步长，\nabla f(x_k) 是 f 在 x_k 处的梯度，k = 1, 2, \ldots 是迭代计数器。上述最小化问题的解显然就是 x_{k+1} = x_k - t_k\nabla f(x_k)，这正是标准的梯度下降更新。

当加入L1正则化项后，子问题变为：

目标函数 F 是两个凸函数的和，因此 F 本身也是凸函数。更关键的是，F 关于 x 的各个分量是可分的，即 F(x) = \sum_i F_i(x_i)，每个 F_i 只依赖于标量 x_i，这使得高维优化问题可以分解为多个独立的一维标量优化问题。

可分问题的分解

引入记号 y = x_k - t_k\nabla f(x_k)，即当前点沿负梯度方向移动后的位置，则子问题可以简化为：

利用范数的定义 \|x\|_2^2 = \sum_i x_i^2 和 \|x\|_1 = \sum_i |x_i|，目标函数可以写成各分量的和：

由于各分量之间没有耦合，整体最小化等价于对每个分量分别最小化。因此我们得到逐坐标的标量表达式：

问题转化为：给定 y_i，求标量函数 \frac{1}{2t_k}(x - y_i)^2 + \lambda|x| 的最小值点。

标量问题的求解

考虑标量函数：

这个函数是凸的，因为 |x| 虽然在零点不可微但是凸函数，(x-y)^2 是可微的凸函数，两个凸函数的非负加权和仍然是凸函数。

记 f 的最小值点为 x^*。由于 |x| 在零点不可微，不能直接对整个函数求导并令其为零，而需要使用次梯度的最优性条件。对于凸函数，x^* 是最小值点当且仅当 f 在 x^* 处是次可微的，并且零属于 f 在 x^* 处的次微分集合：

次梯度最优性条件

对于标量函数 f(x) = \frac{1}{2t_k}(x - y)^2 + \lambda|x|，需要分别计算两部分的次微分。绝对值函数 |x| 的次微分是符号函数 \text{sgn}(x)，它是一个集值映射：当 x > 0 时取值为 1，当 x < 0 时取值为 -1，当 x = 0 时取值为整个区间 [-1, 1]。二次项 (x - y)^2 处处可微，其梯度为 2(x - y)。

因此最优性条件 0 \in \partial f(x^*) 变为：

为简化推导，令 t = 1（一般情况只需将 \lambda 替换为 \lambda t 即可），最优性条件变为：

整理后得到：

这个等式建立了给定值 y 与最优解 x^* 之间的关系，是推导软阈值算子的关键。

从最优性条件推导软阈值算子

方程 y = x^* + \lambda\text{sgn}(x^*) 将 y 表示为 x^* 的函数，但实际问题中给定的是 y，需要求出对应的最优解 x^*。因此需要将这个关系反转，把 x^* 表示为 y 的函数。

几何上，这种反转可以通过交换函数图像的 x 和 y 坐标轴来实现。为简化分析，取 \lambda = 1。

首先看符号函数 \text{sgn}(x) 的图像：当 x > 0 时函数值为 1，当 x < 0 时函数值为 -1，在 x = 0 处有一个从 -1 到 1 的跳跃。恒等函数 x 的图像是过原点的45度直线。将两者相加得到 x + \lambda\text{sgn}(x) 在 \lambda = 1 时的图像：当 x > 0 时，函数值为 x + 1，这是一条斜率为1、截距为1的直线；当 x < 0 时，函数值为 x - 1，这是一条斜率为1、截距为-1的直线；在 x = 0 处，由于 \text{sgn}(0) \in [-1, 1]，函数值可以取遍区间 [-1, 1] 中的任何值。

软阈值算子

将上述图像的坐标轴交换，即原来的横轴变为纵轴、纵轴变为横轴，就得到了软阈值算子 \mathcal{T} 的图像。

从图中可以看到：当输入 x 的绝对值小于1时，输出为0；当 x > 1 时，输出为 x - 1；当 x < -1 时，输出为 x + 1。这种操作将绝对值小于阈值的分量置零，将绝对值大于阈值的分量向零收缩一个阈值的距离。

软阈值算子的数学表达式为：

这里 (\cdot)_+ 表示取正部运算，即 (a)_+ = \max(a, 0)。对于一般的阈值参数 \lambda，软阈值算子为：

在MATLAB中可以用向量化的方式实现：sign(x).*(max(abs(x)-lambda,0))。

迭代软阈值算法

综合前面的推导，L1正则化最小二乘问题

可以通过以下更新规则求解：

其中软阈值算子逐分量作用于向量。第 i 个分量的更新公式为：

这个算法称为迭代收缩阈值算法ISTA。每一步迭代首先沿着数据拟合项的负梯度方向移动，然后对结果应用软阈值操作。从近端算子的角度来看，ISTA就是将近端梯度更新应用于L1正则化最小二乘问题。

算法效果演示

通过MATLAB仿真可以观察ISTA的收敛行为。

目标函数值随迭代次数的变化曲线显示，算法在前几百次迭代中快速下降，然后逐渐趋于平稳。真实信号 x_i 是一个稀疏向量，只在少数坐标位置有非零值。初始猜测 x_0 是一个全零向量或随机向量。经过1000次迭代后，估计结果 x_{1000} 几乎完美地恢复了真实信号的稀疏结构，非零分量的位置和幅值都与真实值高度吻合。

近端算子的一般定义

近端算子是处理非光滑优化问题的核心工具。对于函数 f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}，其近端算子定义为：

其中参数 \lambda > 0 控制正则项与距离项之间的权衡。这个定义的直观含义是：在 v 的附近找一个点，使得该点的函数值 f(x) 较小，同时又不能离 v 太远。\lambda 越大，对 f(x) 的重视程度越高；\lambda 越小，越倾向于待在 v 附近。

近端算子的适用范围很广：f 可以是非光滑函数，可以包含隐式约束（通过在约束集外取值 +\infty 实现），也可以是任意凸函数。计算近端算子本身是一个凸优化问题，可以用标准方法如BFGS求解，但对于许多常见的函数，近端算子有解析解或高效的专用算法。例如L1范数的近端算子就是软阈值算子，可以在线性时间内计算。

近端算子与投影

近端算子是凸集投影的推广。考虑凸集 \mathcal{C} 的指示函数 I_\mathcal{C}(x)，它在集合内取值为0，在集合外取值为 +\infty。这个指示函数的近端算子就是到凸集的投影：

投影的含义是在凸集 \mathcal{C} 中找到离给定点 v 最近的点。由于指示函数在集合外为无穷大，最小化问题的解必然落在集合内部或边界上。

一个具体例子是向盒子约束集 \mathcal{C} = \{x \mid l \preceq x \preceq u\} 的投影，其中 l 和 u 分别是下界和上界向量。投影可以逐分量独立计算：

这个操作将每个分量截断到对应的上下界范围内：如果分量已经在界内则保持不变，如果超出上界则截断到上界，如果低于下界则截断到下界。

几种常见函数的近端算子

二次函数的近端算子

对于二次函数 f(x) = \frac{1}{2}x^T P x + q^T x + r，其中 P 是对称正定矩阵，q 是向量，r 是常数，近端算子有解析解：

这个结果来自于对近端算子定义中的优化问题求解。将 f(x) 代入并对 x 求导令其为零，可以得到一个线性方程组，其解就是上述形式。

计算这个近端算子的代价取决于矩阵 P 的性质。如果 P 是稠密矩阵且使用直接方法，首先需要 O(n^3) 次浮点运算来分解矩阵 I + \lambda P，之后每次计算近端算子只需 O(n^2) 次运算来做回代求解。如果 P 是稀疏矩阵或具有特殊结构（如对角、三对角、循环矩阵等），可以利用这些结构大幅降低计算成本。另一种选择是使用共轭梯度CG或LSQR等迭代方法，此时可以用 v 作为初始点进行热启动，因为连续迭代之间的 v 值变化通常不大，这样可以加速收敛。

可分和函数的近端算子

如果函数 f 关于变量的各个分块是可分的，即 f(x) = \sum_{i=1}^{N} f_i(x_i)，那么整体的近端算子可以分解为各分块近端算子的独立计算：

这种可分性是并行和分布式近端算法的基础，因为各分块的近端计算可以同时进行而无需通信。

以L1范数 f = \|\cdot\|_1 为例，由于 \|x\|_1 = \sum_i |x_i|，L1范数是完全可分的。其近端算子可以表示为：

这正是之前推导的软阈值算子的分段表示形式。

更一般地，如果 f = \|\cdot\| 是某个范数，\mathcal{B} 是该范数对偶范数的单位球，则近端算子可以用投影来表示：

这个公式建立了范数的近端算子与对偶单位球投影之间的关系。

近端梯度算法的一般形式

考虑复合优化问题：

其中 f 是光滑函数，g: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\} 是闭的真凸函数，可以是非光滑的。近端梯度算法的迭代格式为：

每一步首先沿着光滑项 f 的负梯度方向移动，然后对结果应用非光滑项 g 的近端算子。当 \nabla f 是Lipschitz连续的，Lipschitz常数为 L，且步长选取为 \lambda^k = \lambda \in (0, 1/L] 时，算法以 O(1/k) 的速率收敛，即目标函数值与最优值的差距以迭代次数的倒数衰减。

当 g 取为凸集 \mathcal{C} 的指示函数 I_\mathcal{C} 时，近端算子退化为投影，近端梯度算法变为投影梯度法。这说明投影梯度法是近端梯度算法的特例。近端梯度算法的理论基础可以追溯到1970年代Bruck、Lions和Mercier的工作。

深度学习基础

深度学习的表达能力

在图像处理和计算机视觉中，许多任务需要建模复杂的视觉先验或映射关系，这些关系难以用简单的数学公式显式描述。深度学习提供了一种通过数据驱动的方式构建复杂函数的能力。例如图像分类任务中，深度神经网络 f_\theta 可以将一张狗的图片映射为标签DOG。网络的本质是一系列基本线性变换和非线性函数的链式组合，通过大量参数 \theta 的调整，可以逼近任意复杂的输入输出映射。

神经网络的基本结构

神经网络优化的基本形式是最小化预测输出与真实标签之间的差异：

其中 x 是输入，y 是目标输出，w 是网络的权重参数，Network(w|x) 是给定参数 w 时网络对输入 x 的输出。

单个人工神经元的结构包含以下组成部分：多个输入 x_1, x_2, \ldots, x_n 分别乘以对应的权重 w_1, w_2, \ldots, w_n，然后求和得到 \sum_{i=1}^{n} x_i w_i；这个加权和再通过一个非线性激活函数 f，产生输出 y_j = f(\sum_{i=1}^{n} x_i w_i)。激活函数的作用是引入非线性，常见的选择包括sigmoid、tanh和ReLU等。

将多个神经元按层组织就构成了深度神经网络。网络包含输入层、若干隐藏层和输出层。在输入层，数据通过权重 w_i 传递到第一隐藏层的各个节点 j；第一隐藏层的输出为 y_j = f(\sum x_i w_i)。然后通过权重 w_j 传递到第二隐藏层的节点 k，输出为 y_k = f(\sum x_j w_j)。最后通过权重 w_k 传递到输出层的节点 l，产生最终输出 y_l = f(\sum x_k w_k)。层数越深，网络能够表达的函数越复杂。

卷积神经网络

卷积神经网络CNN是处理图像数据的标准架构。以手写数字识别为例，输入是 28 \times 28 \times 1 的灰度图像。

网络的第一层Conv_1是卷积层，使用 5 \times 5 的卷积核和valid填充方式，将输入转换为 24 \times 24 \times n_1 的特征图，其中 n_1 是卷积核的数量。接着是 2 \times 2 的最大池化层Max-Pooling，将特征图下采样为 12 \times 12 \times n_1。

第二层Conv_2再次应用 5 \times 5 卷积核，得到 8 \times 8 \times n_2 的特征图。经过第二个 2 \times 2 最大池化后，特征图大小变为 4 \times 4 \times n_2。

然后将特征图展平Flattened为一维向量，包含 n_3 个单元。这个向量经过全连接层fc_3，配合ReLU激活函数和dropout正则化，映射到fc_4层。最终输出层有10个节点，分别对应数字0到9的分类概率。

卷积层的优势在于参数共享和局部连接：同一个卷积核在整个图像上滑动，大幅减少了参数数量；每个输出只依赖于输入的局部区域，符合图像的局部相关性特点。池化层通过下采样减少特征图的尺寸，提供一定程度的平移不变性。

监督学习的形式化框架

监督学习的目标是从带标签的数据中学习一个从输入到输出的映射关系。整个框架涉及以下核心要素。

训练集 \mathcal{D} 由 n 个输入-输出对组成：

其中 \mathcal{X} 是输入空间，\mathcal{Y} 是输出空间。每个样本 (x_i, y_i) 独立同分布地采样自某个未知的联合分布 \mathbb{P}_{(X,Y)}。在图像分类任务中，x_i 是一张图片（如狗的照片），y_i 是对应的标签（如DOG）。

预测函数 f_\theta: \mathcal{X} \to \mathcal{Y} 将输入映射到输出，其行为由参数 \theta 决定。在深度学习中，f_\theta 是一个神经网络，\theta 包含网络中所有可学习的权重和偏置。参数 \theta 属于某个参数空间 \Theta。

监督学习的优化目标是找到最优参数 \theta，使得预测函数在训练数据上的损失最小：

这个公式有几层含义。求和符号 \sum_{(x_i, y_i) \in \mathcal{D}} 遍历训练集中的所有样本。对于每个样本，x_i 是输入（图片），通过神经网络 f_\theta 得到预测输出 f_\theta(x_i)。损失函数 \mathcal{L} 衡量预测 f_\theta(x_i) 与真实标签 y_i 之间的差距，同时可能还依赖于参数 \theta 本身（例如正则化项）。整个优化过程在参数空间 \Theta 中搜索使总损失最小的参数值。

训练集上的求和实际上是对未知数据分布期望值的经验估计。由于无法直接计算真实分布下的期望 \mathbb{E}_{(X,Y) \sim \mathbb{P}}[\mathcal{L}(f_\theta(X), Y, \theta)]，我们用有限样本的平均来近似它。当样本量足够大时，根据大数定律，这个经验平均会收敛到真实期望。

求解监督学习问题的核心工具

求解上述优化问题需要两个主要工具。

第一个是随机梯度下降SGD。由于训练集可能包含数百万个样本，每次迭代计算所有样本的梯度和是不现实的。SGD的思想是每次只用一个或一小批样本来估计梯度，从而大幅降低每次迭代的计算成本。虽然每步的梯度估计有噪声，但在期望意义下仍然指向正确的下降方向，最终可以收敛到最优解的邻域。

第二个是链式法则。神经网络是由多层函数复合而成的复杂映射，要计算损失函数对网络参数的梯度，需要通过链式法则将梯度从输出层逐层传播回输入层。链式法则的数学表达为：

这个公式说明，如果 f 通过中间变量 x 依赖于 y，那么 f 对 y 的偏导数等于 f 对 x 的偏导数乘以 x 对 y 的偏导数。在神经网络中，这种链式关系可以递归应用于任意深度的网络结构，使得即使是包含数十亿参数的模型也能高效计算梯度。这正是反向传播算法的理论基础。

第六讲总结

监督学习的成熟性

深度学习MRI重建在监督学习设置下已经成熟：测试阶段可以在较低计算成本下获得更高的图像质量；对不同成像对比度、信噪比和场强具有鲁棒性；4倍加速和8倍加速赛道学习了不同的网络架构；FIRENET等展开架构具有最好的稳定性。

XPDNet解决方案

XPDNet在2020年Brain fastMRI挑战赛中排名第二，在学术界排名第一。该方法结合了基于物理的知识和深度学习的进展（如MWCNN架构）。XPDNet支持非笛卡尔采样和3D成像，并且可以扩展用于离共振伪影校正。

自监督数据欠采样

当没有完全采样数据可用时，自监督方法与展开神经网络兼容。零样本学习和数据增强是相关的研究方向，但需要注意潜在的数据偏差问题。