发射断层扫描图像重建

本课程属于医学成像（Medical Imaging）课程单元中的发射断层扫描（Emission Tomography，简称ET）部分。发射断层扫描是核医学成像的核心技术，主要包括正电子发射断层扫描（PET）和单光子发射计算机断层扫描（SPECT）。整个课程的核心内容围绕断层重建（Tomographic reconstruction）展开，而重建方法主要分为两大类：解析重建（Analytical reconstruction）和迭代重建（Iterative reconstruction）。

发射断层数据的特殊性

发射断层扫描所采集到的数据具有几个区别于其他成像模态的独特性质。首先，采集到的原始数据并不是我们真正想要的物理量本身，而是经过投影变换后的间接测量值。在这个过程中，数据受到多种物理效应的影响，包括光子在组织中的衰减、来自散射和随机符合的背景事件、以及PET成像中正电子在湮灭前的飞行距离。此外，由于核衰变本质上是随机过程，采集到的光子计数服从泊松分布，这意味着数据本身带有统计噪声。

其次，发射断层扫描的目标不仅仅是获得定性的图像对比度，而是要进行定量的功能估计。这意味着重建出的图像中每个像素的数值应当反映真实的生理参数，比如放射性示踪剂的浓度分布。

第三，这是一个大规模计算问题。采集端的数据量（即正弦图，Sinograms）通常在 \mathcal{O}(100M, 10B) 个数据点的量级，其中 M 表示百万，B 表示十亿。重建的目标是估计图像空间中 \mathcal{O}(10M, 100M) 个参数，即每个体素的放射性浓度值。如此庞大的数据量和参数空间对重建算法的计算效率提出了很高的要求。

综合以上特点，可以将发射断层重建问题定义为：求解一个大规模的、定量的、不适定的断层逆问题，且数据服从泊松统计特性。

正问题与逆问题的框架

理解断层重建必须从正问题（Direct Problem）和逆问题（Inverse Problem）的关系入手。正问题描述的是采集过程的物理模型：给定体内放射性示踪剂的真实空间分布（这是我们想知道但未知的量），通过采集系统的物理模型，可以预测探测器将会记录到什么样的数据。这个从未知分布到采集数据的映射过程就是正问题，也称为采集模型（Model of acquisition）。

逆问题则是正问题的反向过程：我们已经拥有了探测器记录的采集数据，目标是反推出体内的放射性分布。这就是重建（Reconstruction）要解决的问题。然而，逆问题通常是不适定的，意味着解可能不唯一、对噪声敏感、或者不稳定，因此需要采用适当的数学方法来获得稳定且有意义的解。

好的重建模型应具备的性质

在设计重建算法时，需要在模型的简单性和准确性之间取得平衡。一方面，模型需要足够简单，以便在合理的计算时间内完成大规模断层重建任务。另一方面，模型需要足够准确，才能实现精确的定量分析。

准确的采集模型需要考虑两个层面的物理过程。第一个层面是采集过程中的物理效应，包括光子与患者组织的相互作用（如衰减和散射）、探测器的响应特性、以及成像系统的几何结构。第二个层面涉及示踪剂在体内的药代动力学行为，这可以用数学模型来描述示踪剂随时间的分布变化规律。

改进采集模型的准确性可以显著提升重建后图像的质量和定量精度。这也是发射断层成像研究中的一个持续发展方向。

图像重建算法分类

图像重建算法可以根据对采集模型的使用方式分为两大类：知识驱动方法和数据驱动方法。

知识驱动方法

知识驱动方法的核心前提是采集模型 A 是显式已知并且可以被数学构建的。这里的 A 代表从图像空间到数据空间的正向映射，描述了给定一个放射性分布 f 时，探测器会采集到什么样的数据 g。在这个框架下，重建问题可以表述为：已知 g = A(f)，求解 f。

知识驱动方法又可细分为解析方法和迭代方法。

解析方法适用于采集模型相对简单、存在闭式逆变换的情况。其基本思想是直接对采集模型进行解析求逆，重建公式可以写成：

其中 f 是待重建的图像（即放射性示踪剂的空间分布），g 是采集到的测量数据，A^{-1} 表示采集算子的逆算子。以PET成像为例，正问题的采集模型采用线积分模型，该模型假设探测器测量到的信号等于沿探测线路径对放射性分布的线积分。相应地，逆问题就是对线积分模型进行解析求逆，这正是滤波反投影算法的数学基础。

然而，当采集模型过于复杂而无法直接求逆，或者根本不存在闭式逆变换时，就需要采用迭代方法。迭代重建不试图一步到位求出解，而是从一个初始估计出发，逐步迭代逼近真实解。迭代过程遵循以下流程：首先给定初始估计 f^{[0]}，然后在第 p 次迭代中，利用当前估计 f^{[p]} 通过采集模型计算出预期的测量数据 A(f^{[p]})，将其与实际采集数据 g 进行比较，根据比较结果计算修正量 \Delta f^{[p]}，最后按照更新公式：

得到新的估计。这个过程不断重复，直到满足某个收敛条件为止。这里需要注意更新公式中使用的是乘法运算，这与某些迭代算法（如最大似然期望最大化算法）的特定形式有关。采集模型在迭代框架中充当模拟器的角色，用于预测当前估计图像会产生什么样的投影数据。

数据驱动方法

与知识驱动方法不同，数据驱动方法不要求显式知道采集模型 A 的数学形式。

直接深度学习方法完全从数据中学习逆映射关系。通过构建一个深度神经网络，并使用包含大量配对样本（测量数据和对应真实图像）的训练数据库进行训练，网络学会直接从测量数据预测重建图像：

这种方法的优势在于无需对物理过程建模，但需要大量高质量的训练数据，且模型的可解释性较差。

混合方法试图结合知识驱动和数据驱动两种范式的优点。物理信息深度学习是混合方法的一个代表，它利用已知且可信的物理模型和噪声分布模型来约束问题的结构，同时引入数据驱动的模型来提供正则化。这种方式既保留了物理模型的可解释性和泛化能力，又能通过学习来补偿模型的不完善之处。

基于线积分模型的解析重建

解析重建技术无论是在透射断层扫描还是发射断层扫描中，都建立在线积分模型的基础之上。理解这个模型是掌握所有解析重建算法的前提。

线积分模型的数学表述

重建的视场（Field of View）定义为 \mathbb{R}^n 空间中的一个紧支撑区域 \mathfrak{C}，其中 n=2 对应二维重建，n=3 对应三维重建。待重建的对象是一个定义在该区域上的连续函数 f(\mathbf{r}): \mathfrak{C} \subset \mathbb{R}^n \to \mathbb{R}，这个函数描述了空间中每个位置 \mathbf{r} 处某个物理量的分布，具体物理含义取决于成像模态。

测量数据也是一个连续函数 g(L): \mathbb{T}^{(2n-2)} \to \mathbb{R}，其中 \mathbb{T}^{(2n-2)} 表示投影线的流形。在二维情况下，每条投影线需要两个参数来确定，因此流形维度为 2 \times 2 - 2 = 2；在三维情况下，流形维度为 2 \times 3 - 2 = 4。线积分模型将图像函数与测量数据联系起来：

这个公式表明，沿着任意一条投影线 L，测量值 g(L) 等于函数 f 沿该线的积分。这就是断层成像的核心数学模型，也称为Radon变换。

透射断层扫描中的线积分模型

在X射线CT中，待重建的函数 f(\mathbf{r}) 代表光子的线性衰减系数 \mu(\mathbf{r})，它描述了组织对X射线的吸收能力。成像系统由旋转的X射线管和探测器阵列组成，X射线管发射初始强度为 I_0 的射线束，穿过人体后被探测器接收，测得强度为 I。

X射线穿过物质时的衰减遵循比尔-朗伯定律（Beer-Lambert Law）：

其中 L(\beta, \alpha) 表示由角度参数 \beta（X射线管的旋转角度）和 \alpha（射线在扇形束中的位置角）确定的一条射线路径。这个公式说明射线强度按指数规律衰减，衰减程度取决于沿路径的衰减系数积分。

对上式两边取对数并取负号，就得到CT的线积分模型：

这样，通过对测量数据进行对数变换，就把非线性的指数衰减关系转化为线性的线积分模型，从而可以应用解析重建算法。

单光子发射断层扫描中的线积分模型

SPECT成像中，f(\mathbf{r}) 表示单光子发射浓度，即注入体内的放射性示踪剂在各位置的浓度分布。探测器是围绕患者旋转的伽马相机，测量从体内发射出的伽马光子数 y。投影线 L(s, \psi) 由位置参数 s 和角度参数 \psi 确定。

与CT不同，SPECT的光子是从体内向外发射的，在到达探测器之前会经历组织的衰减。因此，位于 \mathbf{r} 处发射的光子，需要乘以一个衰减因子才能得到实际被探测到的贡献。这个衰减因子等于 \mathrm{e}^{-\int_{L(s,\psi,\mathbf{r})} \mu(\mathbf{r}')\mathrm{d}\mathbf{r}'}，表示从发射点 \mathbf{r} 到探测器这段路径上的累积衰减。SPECT的衰减线积分模型为：

其中 \mu(\mathbf{r}) 是伽马射线的线性衰减系数。这个模型比CT复杂，因为衰减因子依赖于发射点的位置，使得整个积分不再是简单的线积分形式，这也是SPECT定量重建更具挑战性的原因之一。

正电子发射断层扫描中的线积分模型

PET成像的物理原理与SPECT有本质区别。f(\mathbf{r}) 表示正电子发射浓度，放射性示踪剂（如FDG）在衰变时发射正电子，正电子与周围电子湮灭后产生一对方向相反的511 keV伽马光子。当探测器环上的两个探测器几乎同时（在一个很短的时间窗内）检测到光子时，就记录一次符合事件，测量值 y 是检测到的符合计数。连接这两个探测器的线称为符合响应线（Line of Response，LOR），记为 L(a,b)，其中 a 和 b 是两个探测器的位置。

由于两个光子沿相反方向传播，每个光子各自经历从湮灭点到各自探测器的衰减。但无论湮灭发生在响应线上的哪个位置，两个光子的总衰减路径始终覆盖整条响应线。因此，总的衰减因子为：

这里衰减因子 \mathrm{e}^{-\int_{L(a,b)} \mu(\mathbf{r})\mathrm{d}\mathbf{r}} 可以提到积分号外面，因为它只取决于整条响应线，而与发射点位置无关。如果已知衰减分布 \mu(\mathbf{r})（通常通过CT扫描获得），就可以对测量数据进行衰减校正，得到PET的线积分模型：

经过衰减校正后，PET数据也服从标准的线积分模型，可以使用与CT相同的解析重建算法。这也解释了为什么PET/CT联合扫描设备如此普及：CT提供衰减图用于PET的衰减校正，同时提供解剖结构信息。

扫描几何结构

断层扫描系统根据射线束的几何形状可以分为三种主要类型。

平行束几何是最简单的情况，所有射线彼此平行。这种几何结构出现在使用平行孔准直器的PET和SPECT系统中。平行孔准直器由大量平行排列的小孔组成，只允许垂直于探测器表面的光子通过，从而形成平行的投影线。平行束几何的数学处理最为简单，是理解断层重建理论的基础。

扇形束几何中，射线从一个点源发出呈扇形展开。单层CT扫描仪和配备扇形束准直器的SPECT系统采用这种几何结构。在CT中，X射线管作为点源，发出的射线束在一个平面内呈扇形分布，对面是弧形或直线排列的探测器阵列。

锥形束几何是扇形束的三维扩展，射线从点源出发在空间中呈锥形展开。多层CT（即螺旋CT）和使用针孔准直器的SPECT系统属于这种几何。锥形束能够在单次旋转中覆盖较大的轴向范围，提高了扫描效率，但重建算法也更为复杂。

平行束坐标系统

为了数学上精确描述平行束几何中的投影线，需要建立合适的坐标系统。在三维空间中，设 \mathbf{e}_x、\mathbf{e}_y、\mathbf{e}_z 为笛卡尔坐标系的基向量，其中 \mathbf{e}_z 通常取为断层扫描的轴向方向。

投影线的方向用单位向量 \mathbf{e}_\theta \in \mathbb{S}^2 表示，其中 \mathbb{S}^2 是三维空间中的单位球面。方向向量可以用两个角度参数化：方位角 \psi 和极角 \phi。方位角 \psi 描述方向在 xy 平面上的投影与 x 轴的夹角，极角 \phi 描述方向与 z 轴的夹角。

给定方向 \mathbf{e}_\theta，与之垂直的平面定义为：

这个平面称为投影平面，所有方向为 \mathbf{e}_\theta 的平行射线都垂直穿过这个平面。投影平面内的点用向量 \mathbf{s} 表示，\mathbf{s} 唯一确定了投影平面内的一个位置。

一条具体的投影线 L(\mathbf{e}_\theta, \mathbf{s}) 可以参数化为：

这表示从投影平面上的点 \mathbf{s} 出发，沿方向 \mathbf{e}_\theta 延伸的直线，参数 t 取遍所有实数。

三维空间中所有可能的投影线构成一个四维流形：

流形维度为4是因为：方向 \mathbf{e}_\theta 在二维球面上需要2个参数，投影平面内的位置 \mathbf{s} 需要2个参数，共计4个自由度。在PET成像中，探测器环确定了符合响应线（LOR）的集合，每条LOR对应流形中的一个点。

二维与三维成像的区别

二维成像和三维成像在数据采集和重建策略上存在本质差异。

二维成像

在二维成像模式下，所有测量的投影线都垂直于断层轴 \mathbf{e}_z。在PET中，这通过在探测器环之间放置铅或钨制成的隔板（septa）来实现，隔板阻挡斜向入射的符合事件，只接受同一环内或相邻环之间的符合。这种限制使得三维体数据 f(\mathbf{r}) 可以分解为一系列相互独立的二维平面，每个横断面可以单独重建。

在二维情况下，方向向量 \mathbf{e}_\theta 位于单位圆 \mathbb{S}^1 上，只需要一个角度参数 \psi 即可确定。对于任意给定的方向 \mathbf{e}_\theta，由径向位置 s 确定的所有投影都被完整测量。这意味着二维数据是完备的，每个横断面的重建是一个标准的二维逆问题。

三维成像

三维成像模式下，投影线不再局限于垂直于 \mathbf{e}_z 的方向。在PET中，移除隔板后可以接受任意两个探测器之间的符合事件，大大增加了灵敏度。然而，三维成像带来了新的挑战：整个体数据 f(\mathbf{r}) 必须作为一个整体同时重建，无法像二维那样逐层处理。

三维情况下，方向向量 \mathbf{e}_\theta 位于二维球面 \mathbb{S}^2 上，需要两个角度参数。但实际的断层扫描仪并非球形，探测器的几何结构通常是圆柱形的，这意味着并非所有方向 \mathbf{e}_\theta \in \mathbb{S}^2 都能被测量到。此外，对于任意给定的方向 \mathbf{e}_\theta，也不是投影平面内所有位置 \mathbf{s} 的投影都能被测量，这会产生截断效应。

尽管二维和三维解析重建算法共享相同的基本概念（如傅里叶切片定理、滤波反投影等），但由于数据完备性和几何约束的差异，它们的具体实现方式有本质不同。

n维平行射线变换

为了统一处理二维和三维重建问题，可以定义一般的 n 维平行射线变换。设 f: \mathbb{R}^n \to \mathbb{R} 是待重建的图像函数，其 n 维平行射线变换 \mathbf{X}f 是定义在投影线流形上的函数：

变换的具体定义为：

这个积分沿着方向为 \mathbf{e}_\theta、经过点 \mathbf{s} 的直线对函数 f 进行积分。参数 t 是沿射线方向的位移量，积分范围覆盖整条射线。

数据采集模型可以简洁地表示为：

即测量数据 g 就是图像函数 f 的平行射线变换。断层重建的任务就是对平行射线变换求逆，从测量数据 g 恢复原始图像 f。这个逆问题的解法构成了解析重建理论的核心内容。

二维与三维平行射线变换的具体形式

平行射线变换的一般定义为：

这个抽象定义在不同维度下有具体的参数化形式。

二维情形：一维投影

在二维情况下，图像是定义在 \mathbb{R}^2 上的函数 f(x, y)，而投影数据 g(\psi, s) 定义在二维流形 \mathbb{T}^2 上。方向向量 \mathbf{e}_\theta 由单一角度 \psi 确定，可以写成 \mathbf{e}_\theta = (\sin\psi, \cos\psi)。与之垂直的方向为 \mathbf{e}_\theta^\perp = (\cos\psi, -\sin\psi)，投影平面上的点可以写成 \mathbf{s} = s \cdot \mathbf{e}_\theta^\perp = (s\cos\psi, -s\sin\psi)。

将这些代入一般定义，沿射线的点可以参数化为 \mathbf{s} + t\mathbf{e}_\theta = (s\cos\psi + t\sin\psi, -s\sin\psi + t\cos\psi)。因此二维平行射线变换的显式形式为：

这个积分对固定的角度 \psi 和径向位置 s，沿着一条与 \mathbf{e}_\theta 方向平行的直线对图像函数求积分。对于每个角度 \psi，当 s 取遍所有值时，得到一个一维函数，这就是该角度下的投影。

三维情形：二维投影

在三维情况下，图像是定义在 \mathbb{R}^3 上的函数 f(x, y, z)，投影数据 g(\psi, \phi, s_1, s_2) 定义在四维流形 \mathbb{T}^4 上。方向向量需要两个角度参数 \psi（方位角）和 \phi（极角）：

投影平面是垂直于 \mathbf{e}_\theta 的二维平面，平面内的点用两个坐标 (s_1, s_2) 描述。需要在投影平面内选取两个正交的基向量 \mathbf{e}_{\theta_1}^\perp 和 \mathbf{e}_{\theta_2}^\perp，则投影平面内的位置向量为：

对于每个固定的方向 (\psi, \phi)，当 (s_1, s_2) 取遍投影平面时，得到一个二维函数，这就是该方向下的二维投影。

正弦图

正弦图（Sinogram）是二维断层扫描中存储投影数据的标准格式。将所有角度的一维投影按角度顺序排列，横轴为径向位置 s，纵轴为角度 \psi，就形成了正弦图：

正弦图得名于其特殊的视觉特征：图像中的一个点源在正弦图中会呈现为一条正弦曲线。这是因为当投影角度 \psi 变化时，点源的投影位置 s 随 \psi 呈正弦规律变化。具体来说，位于 (x_0, y_0) 处的点源，在角度 \psi 的投影中出现在位置 s = x_0 \cos\psi + y_0 \sin\psi，这正是关于 \psi 的正弦函数形式。

从图示可以看到，当分别在 \psi = 0°、45°、90°、135°、180° 等角度采集投影时，每个角度对应正弦图中的一行。图像中物体的轮廓在不同角度的投影呈现不同形态，所有这些投影组合起来就构成了完整的正弦图。正弦图包含了重建图像所需的全部信息，解析重建的任务就是从正弦图恢复原始图像。

n维中心切片定理

中心切片定理（又称傅里叶切片定理或投影切片定理）是解析重建理论的核心，它揭示了投影数据与图像在频域之间的关系。

傅里叶变换的定义

首先定义两种傅里叶变换。对投影数据 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) 关于投影平面内的坐标 \mathbf{s} 做 (n-1) 维傅里叶变换：

其中 \boldsymbol{\nu}_\perp \in \theta^\perp 是投影平面内的频率变量。对于二维情况，这是一维傅里叶变换；对于三维情况，这是二维傅里叶变换。

对图像函数 f(\mathbf{r}) 做 n 维傅里叶变换：

其中 \boldsymbol{\nu} \in \mathbb{R}^n 是 n 维频率空间中的变量。

定理陈述

对于 \mathbb{R}^n 中的函数 f(\mathbf{r})，其平行射线变换的 (n-1) 维傅里叶变换等于图像的 n 维傅里叶变换在相应超平面上的取值：

这个等式的含义是：对某一方向 \mathbf{e}_\theta 的投影做傅里叶变换，得到的结果恰好等于图像的 n 维傅里叶变换在通过原点、垂直于 \mathbf{e}_\theta 的超平面上的切片。

定理证明

证明从投影的傅里叶变换定义出发：

将平行射线变换的定义代入：

现在做变量替换 \mathbf{r} = \mathbf{s} + t\mathbf{e}_\theta。由于 \mathbf{s} \in \theta^\perp 且 t \in \mathbb{R}^1，当 \mathbf{s} 遍历整个投影平面、t 遍历整条实数轴时，\mathbf{r} 恰好遍历整个 \mathbb{R}^n 空间。另外，由于 \boldsymbol{\nu}_\perp \in \theta^\perp 而 t\mathbf{e}_\theta 垂直于 \theta^\perp，所以 \mathbf{r} \cdot \boldsymbol{\nu}_\perp = (\mathbf{s} + t\mathbf{e}_\theta) \cdot \boldsymbol{\nu}_\perp = \mathbf{s} \cdot \boldsymbol{\nu}_\perp。因此：

证毕。

二维情形的直观理解

在二维情况下，中心切片定理可以这样理解：对角度为 \psi = 45° 的一维投影 \mathbf{X}f(\psi=45°, s) 做一维傅里叶变换，得到的一维频谱 |\mathcal{F}_1\{\mathbf{X}f\}(\psi=45°, \nu_s)| 恰好等于图像的二维傅里叶变换 |\mathcal{F}_2\{f\}(\nu_x, \nu_y)| 沿着通过原点、与投影方向垂直的那条直线上的取值。换言之，每个角度的投影提供了图像频谱的一条径向切片。当采集足够多角度的投影时，就能覆盖整个二维频率平面，从而可以通过逆傅里叶变换重建图像。

三维中心切片定理

三维情况下的中心切片定理与二维情况类似，但涉及更高维度的频率空间。三维图像 f(x,y,z) 在某一方向 \mathbf{e}_\theta 上的二维投影 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s})，对其做二维傅里叶变换 \mathcal{F}_2\{\mathbf{X}f\}(\mathbf{e}_\theta, \nu_s)，得到的结果等于三维图像的傅里叶变换 \mathcal{F}_3\{f\}(\nu_x, \nu_y, \nu_z) 在通过原点、垂直于投影方向 \mathbf{e}_\theta 的平面上的切片。

从图示可以看到，三维频率空间是一个立方体，每个投影方向对应其中一个通过原点的平面。当采集足够多方向的投影时，这些平面切片就能覆盖整个三维频率空间，从而可以通过三维逆傅里叶变换重建图像。

直接傅里叶变换重建

中心切片定理直接导出了一种重建方法。在二维情况下，如果投影 g(\psi, s) = \mathbf{X}f(\psi, s) 在角度区间 \psi \in [0, \pi[ 上被完整测量，那么就能够重建图像 f。角度只需覆盖半圆周 [0, \pi[ 而非完整的 [0, 2\pi[ 是因为相差 \pi 的两个角度提供的是同一条投影线的重复测量（方向相反但积分值相同）。

直接傅里叶变换重建（Direct Inversion Fourier Transform，DIFT）的步骤如下：首先对每个角度 \psi 的一维投影 g(\psi, s) 做一维FFT，得到 \mathcal{F}_1\{g\}(\psi, \nu_s)；根据中心切片定理，这些一维频谱构成了图像二维频谱在极坐标网格上的采样；然后将极坐标网格上的数据重采样到笛卡尔网格上；最后对笛卡尔网格上的二维频谱做二维逆FFT，得到重建图像 f(x,y)。

这种方法的主要困难在于从极坐标网格到笛卡尔网格的重采样。投影数据在频域中自然形成以原点为中心的径向分布（极坐标采样），而标准的FFT算法要求数据在笛卡尔网格上均匀分布。重采样过程需要插值，会引入误差，特别是在高频区域（远离原点处）采样密度较低，插值误差更大。这个缺点限制了DIFT方法在实际中的应用。

三维数据冗余性

中心切片定理还揭示了三维成像中的数据冗余性。根据定理，仅使用极角 \phi = 0（即投影线垂直于 z 轴）、方位角 \psi \in [0, \pi[ 的数据 g(\psi, \phi=0, \mathbf{s})，就足以重建三维图像 f(x,y,z)。这是因为当 \phi = 0 时，投影的二维傅里叶变换对应三维频率空间中 \nu_z = 0 的平面。当 \psi 从 0 变化到 \pi 时，这些平面绕 \nu_z 轴旋转，覆盖了整个三维频率空间（除了 \nu_z 轴本身，但这是零测度集）。

这意味着二维数据采集模式（所有投影线垂直于断层轴）在理论上已经包含了重建三维图像所需的全部信息。那么采集 \phi \neq 0 的斜向投影有什么意义呢？答案是提高信噪比。斜向投影虽然不提供新的独立信息，但提供了数据的冗余测量，可以通过适当的组合来降低噪声。

这里需要注意，在实际的三维成像中，由于断层扫描仪不是球形的，无法覆盖完整的 4\pi 立体角，可能存在物体的截断效应。上述分析假设投影是完整的，没有截断。

n维反投影算子

定义

反投影（Backprojection）是与投影（正投影或X射线变换）对偶的操作。给定投影数据 g: \mathbb{T}^{(2n-2)} \to \mathbb{R}，其 n 维反投影 \mathbf{X}^* g: \mathbb{R}^n \to \mathbb{R} 定义为：

这个公式的含义是：对于图像空间中的任意一点 \mathbf{r}，找出所有经过该点的投影线，将这些投影线上的测量值 g 累加起来。具体来说，\mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta 是点 \mathbf{r} 在垂直于方向 \mathbf{e}_\theta 的投影平面上的投影点，它确定了经过 \mathbf{r} 且方向为 \mathbf{e}_\theta 的那条投影线。积分遍历所有可能的方向 \mathbf{e}_\theta \in \mathbb{S}^{n-1}。

算子 \mathbf{X}^* 是X射线变换 \mathbf{X} 的伴随算子，这可以通过内积关系来验证：

左边是数据空间中 g 和 \mathbf{X}f 的内积，右边是图像空间中 \mathbf{X}^* g 和 f 的内积。这个等式表明 \mathbf{X}^* 确实是 \mathbf{X} 的伴随。

这里需要特别注意：反投影算子不是X射线变换的逆算子，即 \mathbf{X}^* \mathbf{X} f(\mathbf{r}) \neq f(\mathbf{r})。直接对投影数据做反投影得到的图像是模糊的，这将在后面详细分析。

二维反投影

在二维情况下，图像 f(x,y) 定义在 \mathbb{R}^2 上，数据 g(\psi, s) 定义在 \mathbb{T}^2 上。二维反投影的具体形式为：

其中 x\cos\psi + y\sin\psi 正是点 (x,y) 在角度 \psi 对应的投影线上的径向坐标 s。这个积分遍历所有角度 \psi \in [0, 2\pi[，将经过点 (x,y) 的所有投影值累加。

三维反投影

在三维情况下，图像 f(x,y,z) 定义在 \mathbb{R}^3 上，数据 g(\psi, \phi, s_1, s_2) 定义在 \mathbb{T}^4 上。三维反投影需要在整个方向球面 \mathbb{S}^2 上积分：

积分中的 \cos\phi 因子来自球面坐标系的测度。

正投影与反投影的对比

从图示可以直观理解正投影和反投影的区别。正投影是将图像 f(x,y) 沿不同方向积分，得到各角度的一维投影曲线 \mathbf{X}f(\psi, s)。反投影则是相反的过程：将每个角度的投影曲线沿其对应的方向"涂抹"回图像空间，然后将所有角度的涂抹结果叠加，得到反投影图像 \mathbf{X}^* g(x,y)。

反投影图像与原图像相比是模糊的，呈现出低频增强的特点，这是因为每个投影值被均匀分配到整条投影线上，导致了信息的弥散。

反投影的弥散效应

图示直观展示了二维反投影如何逐步累积形成重建图像。第一行从左到右分别是：只有一个角度的反投影（呈现为沿该方向的条纹），两个角度叠加，四个角度叠加，以及更多角度叠加后的结果。可以看到，随着参与累加的角度增多，图像逐渐显现出原始物体的轮廓。

然而，即使使用了所有角度的投影，反投影得到的图像仍然是模糊的。原始图像中的两个椭圆在反投影后变成了弥散的光晕，边缘不再清晰。这种弥散效应是反投影操作的固有特性，必须通过额外的滤波步骤来校正。

正投影与反投影的卷积关系

卷积的定义与性质

在 \mathbb{R}^n 中，两个函数 f(\mathbf{r}) 和 h(\mathbf{r}) 的卷积定义为：

卷积与傅里叶变换之间存在卷积定理：空间域的卷积对应频域的乘积：

核心定理：反投影-正投影的卷积形式

对于 \mathbb{R}^n 中的函数 f(\mathbf{r})，先做正投影再做反投影的结果等于原函数与特定核函数的卷积：

这个定理揭示了反投影图像模糊的数学本质：它是原图像与 \|\mathbf{r}\|^{1-n} 的卷积结果。在二维情况下，核函数为 \|\mathbf{r}\|^{-1} = 1/\|\mathbf{r}\|，这是一个随距离衰减的函数，造成了图像的低通滤波效应。

定理证明

从反投影的定义出发，将正投影 \mathbf{X}f 代入：

再将正投影的线积分定义展开：

做变量替换 u = t - \mathbf{r} \cdot \mathbf{e}_\theta：

这里利用了被积函数关于 u 的对称性，将积分区间从 \mathbb{R}^1 改为 [0, \infty) 并乘以2。

再做变量替换 \mathbf{v} = -u\mathbf{e}_\theta，当 u 从 0 到 \infty、\mathbf{e}_\theta 遍历整个 (n-1) 维球面时，\mathbf{v} 恰好遍历整个 \mathbb{R}^n 空间。雅可比行列式为 \mathrm{d}u \mathrm{d}\mathbf{e}_\theta = \|\mathbf{v}\|^{n-1} \mathrm{d}\mathbf{v}，因此：

这正是卷积的定义形式，证毕。

二维反投影滤波分析

对于二维情况 f: \mathbb{R}^2 \to \mathbb{R}，上述定理变为：

图示展示了原始图像（左）经过 \mathbf{X}^* \mathbf{X} 操作后变成模糊图像（右）的效果。

为了从反投影图像恢复原图像，需要知道核函数 1/\|\mathbf{r}\| 的傅里叶变换。在二维情况下：

利用卷积定理，可以得到原图像频谱与反投影图像频谱的关系：

这个公式表明，要从反投影结果恢复原图像，需要在频域乘以 \|\boldsymbol{\nu}\|，即进行高通滤波来补偿反投影造成的低通效应。

二维滤波反投影算法

滤波反投影定理

滤波反投影（Filtered Backprojection，FBP）是解析重建中最经典的算法。其核心思想是交换滤波与反投影的顺序：先对投影数据进行滤波，再做反投影。

如果投影 g(\psi, s) = \mathbf{X}f(\psi, s) 在角度区间 \psi \in [0, 2\pi[ 上被测量，则图像可以通过以下公式重建：

其中 g^f(\psi, s) 是滤波后的投影，由下式定义：

这个滤波操作的物理含义是：对每个角度 \psi 的一维投影 g(\psi, s) 做一维傅里叶变换，乘以斜坡滤波器 |\nu_\perp|（其幅度与频率的绝对值成正比），再做一维逆傅里叶变换得到滤波后的投影 g^f(\psi, s)。

由于投影的对称性 \mathbf{X}f(\psi, s) = \mathbf{X}f(\psi + \pi, -s)（相差 \pi 的两个角度给出相同的投影线，只是方向相反），只需测量 \psi \in [0, \pi[ 区间的数据就足够了。此时重建公式简化为：

定理证明

证明从逆傅里叶变换公式开始：

将笛卡尔坐标 (\nu_x, \nu_y) 转换为极坐标 (\nu_\perp, \psi)，其中 \nu_x = \nu_\perp \cos\psi，\nu_y = \nu_\perp \sin\psi：

这里 \nu_\perp \mathrm{d}\nu_\perp \mathrm{d}\psi 是极坐标的面积元素。

根据二维中心切片定理，\mathcal{F}_2\{f\}(\nu_\perp\cos\psi, \nu_\perp\sin\psi) = \mathcal{F}_1\{\mathbf{X}f\}(\psi, \nu_\perp)，代入得：

被积函数中的 \nu_\perp 因子正是斜坡滤波器 |\nu_\perp| 在 \nu_\perp > 0 时的形式。通过适当处理负频率部分，可以将这个积分写成完整的滤波反投影形式。

滤波反投影证明的继续

从上一部分的推导继续，将频域积分从笛卡尔坐标转换到极坐标后：

现在需要处理积分限。原式中 \nu_\perp 只取正值，角度 \psi 从 0 到 2\pi。可以利用对称性将积分改写为角度从 0 到 \pi、频率从 -\infty 到 +\infty 的形式。当 \psi 变为 \psi + \pi 且 \nu_\perp 变为 -\nu_\perp 时，指数项中的 \nu_\perp(x\cos\psi + y\sin\psi) 保持不变，因此：

这里 \nu_\perp 因子变成了 |\nu_\perp|，因为当 \nu_\perp 取负值时需要保持被积函数的正确符号。内层积分恰好是一维逆傅里叶变换的形式，因此：

证毕。这个结果表明图像可以通过对滤波后的投影做反投影来重建。

滤波反投影算法步骤

滤波反投影算法的具体实现分为以下步骤。

第一步，对每个角度 \psi 的投影做一维傅里叶变换：

第二步，将频谱乘以斜坡滤波器：

第三步，对滤波后的频谱做一维逆傅里叶变换，得到滤波后的投影：

第四步，对每个角度 \psi，将滤波后的投影反投影到图像空间并累加：

第五步，对所有角度重复步骤一到四，完成图像重建。

这个算法的计算效率很高，因为一维FFT的复杂度为 O(N \log N)，而反投影可以高效地并行实现。

滤波效果的直观理解

以均匀圆盘的投影为例，可以直观理解斜坡滤波器的作用。均匀圆盘的投影 g(\psi, s) 是一个半椭圆形的曲线（实际上对所有角度都相同，因为圆盘具有旋转对称性），在圆盘边缘处平滑过渡到零。

经过斜坡滤波后，g^f(\psi, s) 的形态发生显著变化：在圆盘内部区域变成近似平坦，而在边缘处出现尖锐的负值峰。这些负值的产生是因为斜坡滤波器增强了高频成分，而投影在边缘处的急剧变化包含丰富的高频信息。滤波后投影的正负值在反投影累加时会相互抵消，消除原本会出现在圆盘外部的模糊光晕。

图示对比了不带滤波和带滤波的反投影结果。上半部分展示不带滤波的情况：随着角度增加，反投影逐渐累积，最终图像虽然能看出两个椭圆的大致形状，但整体非常模糊，背景区域有明显的灰色光晕。下半部分展示带滤波的情况：每个角度的滤波投影在反投影时包含正负值（图中可见明暗交替的条纹），这些正负值在累加过程中相互抵消背景区域的虚假信号，最终得到清晰的图像，椭圆边缘锐利，背景干净。

采样问题

实际的断层扫描系统中，投影数据不是连续测量的，而是按照离散网格采样的。径向方向的采样间隔为 \Delta s，角度方向的采样间隔为 \Delta \psi。离散化的投影数据可以表示为 g_{j,k} = g(\psi_j, s_k)。

根据奈奎斯特-香农采样定理，能够无混叠重建的最大频率为：

这意味着采样间隔 \Delta s 决定了重建图像能够分辨的最小细节尺度。如果原图像包含高于 \nu_{\max} 的频率成分，这些成分会混叠到低频区域，造成伪影。

在实际应用中，重建的是经过低通滤波的图像 W_c * f，其中 W_c 是截止频率为 \nu_c \leq \frac{1}{2\Delta s} 的低通滤波器。W_c 的傅里叶变换具有以下性质：

其中窗函数 \Phi 满足 0 \leq \Phi(|\boldsymbol{\nu}|/\nu_c) \leq 1，且当 |\boldsymbol{\nu}|/\nu_c \geq 1 时 \Phi = 0。这保证了只有低于截止频率的成分被保留。

离散滤波反投影与采样的关系

对于图像空间中的低通滤波器 W_c \in \mathbb{R}^2，存在一个定理将其与投影空间中的滤波器联系起来：

其中 w_c \in \mathbb{T}^2 是投影空间中的滤波器，满足 \mathbf{X}^* w_c = W_c。这个关系表明，在图像空间做低通滤波等价于在投影空间做相应的滤波后再反投影。

如果图像空间滤波器的傅里叶变换为 \mathcal{F}_2\{W_c\}(\boldsymbol{\nu}) = \Phi\left(\frac{|\boldsymbol{\nu}|}{\nu_c}\right)，则投影空间滤波器的傅里叶变换为：

这个结果直接给出了离散情况下滤波反投影的实现方式。对于离散投影 g = \mathbf{X}f，重建公式为：

其中滤波后投影的频谱为：

滤波器 |\nu_\perp| \Phi(|\nu_\perp|/\nu_c) 是理想斜坡滤波器 |\nu_\perp| 与窗函数的乘积。窗函数的作用是在截止频率处平滑地将滤波器截断为零，避免直接截断带来的振铃伪影。

连续数据与采样数据的滤波反投影对比

对于连续数据 g = \mathbf{X}f，滤波反投影给出精确重建：

其中滤波后投影的频谱为：

然而对于采样数据 g_{\Delta s} = \mathbf{X}f，由于采样带来的频率限制，重建结果不再是原图像，而是低通滤波后的版本：

滤波后采样数据的频谱为：

这表明采样数据的重建本质上是带限的，无法恢复高于截止频率的细节。

矩形滤波器

矩形滤波器是最简单的窗函数选择，重建公式为：

图示展示了矩形滤波器的三个相关函数：低通滤波器的频率响应（左）、斜坡滤波器（中）、以及两者乘积后的空间域核函数（右）。

低通滤波器 \Phi(|\nu|/\nu_c) 的定义为：

投影空间滤波器的频率响应为：

对应的空间域滤波核函数 w_c(\psi, s) 可以通过逆傅里叶变换得到：

这个核函数在 s = 0 处取最大值，在两侧呈振荡衰减的形态。振荡的存在是矩形窗函数在频域硬截断的直接后果，会在重建图像中产生振铃伪影。

不适定问题

平行射线变换的求逆是一个不适定问题（ill-posed problem）。通过奇异值分解（SVD）分析X射线变换系统，二维和 n 维情况下的研究表明，奇异值 \sigma_{ml} 的衰减规律为：

奇异值以 m^{-1/2} 的速度趋向零，这意味着对应高频成分的奇异值非常小。

滤波反投影解 f_{\text{FBP}} = \frac{1}{2}\mathbf{X}^* g^f 对数据 g 不具有连续依赖性。数据中的微小扰动可能导致重建图像中任意大的误差。在发射断层扫描中，数据扰动的典型来源是泊松随机噪声。由于光子计数服从泊松分布，计数越少噪声越显著，这对低剂量成像构成严峻挑战。

图示展示了不同计数水平下的FBP重建结果。从左到右、从上到下，总计数分别为 10^2、10^3、10^4、10^5、10^6、10^7。当计数仅为 10^2 时，正弦图几乎完全被噪声淹没，重建图像无法辨认；随着计数增加，图像质量逐步改善；当计数达到 10^7 时，重建图像接近真实图像 f(\mathbf{x})。这说明FBP对噪声极其敏感，需要足够高的统计计数才能获得可接受的图像质量。

滤波反投影与噪声放大

图示详细分析了噪声如何在FBP中被放大。左上角的原始正弦图和添加了泊松噪声（159,000个事件）的正弦图对比显示，噪声使得正弦图出现明显的颗粒感。原始图像在FBP重建后变得模糊且布满噪声斑点。

中间的频谱分析揭示了噪声放大的机制。绿色曲线是信号的功率谱，蓝色曲线是测量数据的功率谱，红色曲线是噪声的功率谱。可以看到，信号主要集中在低频区域，而噪声在整个频率范围内相对平坦。当乘以斜坡滤波器（中间图，红色曲线，随频率线性增长）后，高频区域的噪声被大幅放大，而信号在高频区域本来就很弱。结果是滤波后的频谱中，高频部分几乎完全由放大后的噪声主导。

功率谱的角度清楚地表明：斜坡滤波器本质上是在放大高频噪声。这是FBP算法噪声敏感性的根本原因。

高频截止策略

为了控制噪声放大，一个直接的策略是不重建高于某个截止频率 \nu_c 的成分，其中 \nu_c < \nu_{\text{Nyquist}} = \frac{1}{2\Delta s}。通过选择较低的截止频率，可以在滤波后的频谱中抑制高频噪声，代价是牺牲图像的空间分辨率。

图示展示了这种噪声-分辨率权衡。当 \nu_c = 0.2\nu_N 时，图像非常平滑但极度模糊，细节完全丢失；随着 \nu_c 增加到 0.4\nu_N、0.6\nu_N、0.8\nu_N，图像逐渐变得清晰但噪声也逐渐增加；当 \nu_c = 1.0\nu_N（即等于奈奎斯特频率）时，图像包含最多的细节，但噪声也最严重。

这种权衡是断层重建中的基本困境：高截止频率保留更多细节但放大噪声，低截止频率抑制噪声但模糊图像。在实际应用中，需要根据成像任务的具体需求选择合适的截止频率，或者采用更复杂的正则化方法来更好地平衡噪声和分辨率。

切趾窗函数

矩形窗函数在截止频率处的硬截断会引入振铃伪影。切趾（Apodization）技术通过使用平滑过渡的窗函数来缓解这个问题。

图示对比了矩形窗和汉宁窗（Hann window）的效果。上半部分使用矩形窗：信号和噪声的频谱乘以矩形窗（在 \nu_c 处突然截断），再乘以斜坡滤波器，得到滤波后的结果。由于矩形窗在频域的突变，重建图像边缘会出现尖锐的过渡，可能伴随振铃现象。

下半部分使用汉宁窗：窗函数从 \nu_c 开始平滑地过渡到零，而不是突然截断。乘以斜坡滤波器后，滤波后的信号和噪声都更加平滑。重建图像的边缘呈现平滑过渡，避免了振铃伪影，但代价是稍微降低了空间分辨率。

常用低通滤波器

实际应用中有多种窗函数可供选择，它们在噪声抑制和分辨率保持之间提供不同的权衡。

广义汉明窗

参数 \alpha 控制窗函数的形状。当 \alpha = 0.5 时得到汉宁窗（Hann），当 \alpha = 0.54 时得到汉明窗（Hamming），当 \alpha = 1 时退化为矩形窗。\alpha 越小，高频衰减越强，噪声抑制效果越好，但分辨率损失也越大。

谢普-洛根滤波器

谢普-洛根滤波器（Shepp-Logan filter）使用正弦函数调制，在中低频区域保持较高的增益，在接近截止频率时平滑衰减。它在噪声抑制和分辨率保持之间提供了较好的平衡，是CT重建中常用的选择。

n阶巴特沃斯滤波器

巴特沃斯滤波器（Butterworth filter）的阶数 n 控制过渡带的陡峭程度。n = 1 时过渡最平缓，随着 n 增大，滤波器逐渐接近矩形窗的硬截断特性。n = 4 或 n = 10 是常用的选择。

图示展示了这些滤波器的频率响应曲线对比。谢普-洛根滤波器在低频区域保持最高增益，汉宁窗和汉明窗次之，巴特沃斯滤波器的响应取决于阶数。所有滤波器都在截止频率 \nu_c 处或之前开始衰减，最终在 \nu_{\text{Nyquist}} 处达到零。

广义汉明窗的参数效果

广义汉明窗的完整定义为：

参数 \alpha 在 0.5（汉宁窗）和 1.0（矩形窗）之间变化。

图示展示了不同参数组合下的重建效果。最左边是原始图像（无噪声）。当 \alpha = 1、\nu_c = \nu_N 时（矩形窗，截止频率等于奈奎斯特频率），图像噪声最大但分辨率最高。当 \alpha = 0.54、\nu_c = \nu_N 时（汉明窗），噪声有所减少，边缘稍微模糊。当 \alpha = 0.54、\nu_c = 0.5\nu_N 时，噪声进一步减少，但图像明显变模糊。当 \alpha = 0.54、\nu_c = 0.2\nu_N 时，噪声很低，但图像严重模糊，细节几乎完全丢失。

截止频率与空间分辨率的关系

降低截止频率可以减少随机噪声，但同时会降低空间分辨率。这是断层重建中的经典权衡关系。

空间分辨率通常用点扩散函数的半高全宽（Full Width at Half Maximum，FWHM）来度量。FWHM越大，表示图像越模糊，分辨率越低。

图示左侧的曲线展示了不同滤波器的FWHM与截止频率的关系。横轴是截止频率（以 2\nu_N 为单位），纵轴是FWHM（以像素为单位）。可以看到，当截止频率从 0.5（即 \nu_N）降低到 0.2 时，FWHM从约 4 增加到约 12，即分辨率下降了约 3 倍。

三种滤波器中，矩形窗在相同截止频率下具有最小的FWHM（最高分辨率），谢普-洛根次之，汉宁窗的FWHM最大。这与它们对高频成分的保留程度一致：矩形窗保留最多，汉宁窗衰减最强。

图示右侧展示了实际重建图像的对比。随着FWHM从 5.0 增加到 6.9，图像中的小圆点逐渐变大变模糊。矩形窗重建的图像噪声较多但点更清晰，汉宁窗重建的图像更平滑但点更模糊。

信号与噪声的权衡

不同的滤波器和截止频率选择对应于信号-噪声权衡曲线上的不同位置。图示上方展示了这种权衡关系：横轴是信号（与分辨率相关），纵轴是噪声。矩形窗、谢普-洛根和汉宁窗对应曲线上的不同点。

沿着曲线移动，可以选择不同的工作点。曲线右上方意味着高信号但高噪声（高分辨率但噪点多），左下方意味着低噪声但低信号（图像平滑但模糊）。

图示下方展示了三个不同工作点的重建结果，分别显示平均图像和方差图像。平均图像反映偏差（bias），方差图像反映噪声水平。可以看到：右侧工作点（高信号区域）的平均图像最清晰，但方差图像显示噪声分布很强；左侧工作点（低噪声区域）的方差图像很平滑，但平均图像模糊。

这种权衡依赖于图像本身的频率内容。如果图像主要包含低频成分（如大的均匀区域），可以使用较低的截止频率而不损失太多信息；如果图像包含精细结构（高频成分），则需要较高的截止频率来保留这些细节。

滤波器选择策略

选择合适的噪声-分辨率权衡点通常是一个经验性的过程，依赖于对图像质量的主观评估。一般的指导原则是：在一定程度的平滑后，图像质量趋于稳定；但过度平滑会导致空间分辨率严重下降和对比度恢复退化。

图示展示了三种滤波设置的对比：矩形窗配合 \nu_c = \nu_N（最高分辨率但最多噪声）、汉宁窗配合 \nu_c = 0.8\nu_N（适度平滑）、汉宁窗配合 \nu_c = 0.5\nu_N（强平滑但模糊）。

除了主观评估，还可以使用定量的性能指标来指导滤波器选择：信噪比（Signal-to-Noise Ratio，SNR）衡量信号强度与噪声水平的比值；均方根误差（Root Mean Square Error，RMSE）衡量重建图像与参考图像之间的整体偏差；数值观察者或人类观察者研究通过受试者工作特征曲线（ROC curves）评估特定诊断任务的性能。

在临床应用中，最终的滤波器选择取决于具体的诊断任务。例如，检测小病灶可能需要较高的分辨率（较高截止频率），而评估大区域的摄取均匀性可能更适合使用较强的平滑。

三维滤波反投影

三维重建与二维重建共享相同的理论基础，但存在两个关键差异。根据三维中心切片定理，仅使用 \phi = 0（投影线垂直于 z 轴）的完整投影 g(\psi, \phi=0, \mathbf{s})，其中 \psi \in [0, \pi[，就足以重建三维图像 f(x,y,z)。这些投影的二维傅里叶变换 \mathcal{F}_2\{\mathbf{X}f\}(\psi, \phi=0, \nu_\perp) 对应三维频率空间 \mathcal{F}_3\{f\}(\nu_x, \nu_y, \nu_z) 中通过原点的一系列平面，当 \psi 从 0 到 \pi 变化时，这些平面能够覆盖整个三维频率空间。

然而，三维情况与二维存在两个主要区别：第一，三维重建需要满足特定的几何条件才能保证解的存在性和唯一性；第二，满足重建条件的滤波器不再唯一，存在无穷多个有效滤波器。

截断问题与有限孔径

实际的断层扫描系统无法在完整的 4\pi 立体角上采集数据。系统的有限孔径定义为所有能够被测量的投影方向的集合：

对于PET系统，探测器通常排列成有限长度的圆柱形，只能接受极角 \phi 在一定范围内的符合事件：

其中 \phi_{\max} 取决于探测器圆柱的长度与直径之比。

为简化分析，通常假设孔径 \Omega 是平移不变的，即对于 \Omega 中的任意方向 \mathbf{e}_\theta，投影 g(\mathbf{e}_\theta, \mathbf{s}) 要么对所有 \mathbf{s} \in \theta^\perp 都被测量，要么都为零。这个假设在实际中不完全成立（存在截断效应），但对于理论分析是必要的简化。

三维重建条件：奥尔洛夫条件

根据中心切片定理，为了恢复三维频谱中的任意一个频率值 \mathcal{F}_3\{f\}(\boldsymbol{\nu})，需要至少有一个投影方向 \mathbf{e}_\theta 满足 \boldsymbol{\nu} 位于该方向的投影平面内（即 \boldsymbol{\nu} \cdot \mathbf{e}_\theta = 0）。

奥尔洛夫条件（Orlov's condition）给出了完整重建的充分必要条件：对于三维频率空间中的任意频率向量 \boldsymbol{\nu}，单位球面上与 \boldsymbol{\nu} 垂直的大圆必须与孔径区域 \Omega 相交。换言之，单位球面上的每一个赤道圆都必须与 \Omega 有交集。

图示展示了三种不同孔径配置对奥尔洛夫条件的满足情况。二维PET模式（\Omega_0）对应孔径仅为 xy 平面上的一个圆（\phi = 0），虽然理论上可以重建，但数据冗余度为零。三维圆柱形PET（\Omega_{\phi_{\max}}）的孔径是一个球带，覆盖 |\phi| < \phi_{\max} 的区域，只要 \phi_{\max} 足够大就能满足奥尔洛夫条件。两个固定平面探测器（\Omega_{\text{planar}}）的孔径是两个相对的球冠，这种配置也可以满足重建条件。

三维滤波反投影算法

对于三维平行射线投影，最一般的线性平移不变反演方法是三维滤波反投影算法。如果投影 g(\mathbf{e}_\theta, \mathbf{s}) = \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) 在平移不变孔径 \Omega 上被测量，则图像可以通过以下公式重建：

其中 \mathbf{X}^*_\Omega 表示仅在孔径 \Omega 上的反投影，滤波后的投影 g^f 定义为：

其中 \boldsymbol{\nu}_\perp 满足 \mathbf{e}_\theta \cdot \boldsymbol{\nu}_\perp = 0，即频率向量位于投影平面内。与二维情况不同，三维滤波器 H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) 依赖于孔径 \Omega 的几何形状，而且由于数据冗余性的存在，满足重建条件的滤波器不是唯一的。

滤波器的非唯一性

任何有效的三维FBP滤波器 H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) 必须满足以下条件：

对于 \mathbb{R}^3 中的任意频率向量 \boldsymbol{\nu} 成立，其中 \boldsymbol{\nu}_\perp 是 \boldsymbol{\nu} 在投影平面 \theta^\perp 上的投影。

这个条件的物理含义是：对于频率空间中的每一点 \boldsymbol{\nu}，所有能够测量到该频率的投影方向（由 \delta(\boldsymbol{\nu} \cdot \mathbf{e}_\theta) = 0 确定，即 \boldsymbol{\nu} 垂直于 \mathbf{e}_\theta）上的滤波器响应的加权和必须等于 1。

满足此条件的滤波器有无穷多个。对于一致性数据（无噪声），所有有效滤波器给出相同的重建结果。但它们在噪声传播特性上存在差异：不同滤波器对数据中噪声的放大方式不同，导致重建图像的噪声特性不同。因此，在实际应用中需要根据噪声特性选择合适的滤波器。

滤波器有效性的证明

证明从三维FBP公式两边取三维傅里叶变换开始：

做变量替换，将 \mathbf{r} 分解为平行于 \mathbf{e}_\theta 的分量 r_\parallel = \mathbf{r} \cdot \mathbf{e}_\theta 和垂直于 \mathbf{e}_\theta 的分量 \mathbf{r}_\perp = \mathbf{r} - r_\parallel \mathbf{e}_\theta：

对 r_\parallel 的积分给出关于 \mathbf{e}_\theta \cdot \boldsymbol{\nu} 的狄拉克函数。再利用中心切片定理 \mathcal{F}_2\{g\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp)：

当 \mathbf{e}_\theta \cdot \boldsymbol{\nu} = 0 时，\boldsymbol{\nu}_\perp = \boldsymbol{\nu}，因此右边的 \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp) 可以提到积分号外。为使等式成立，剩余的积分必须等于 1，这正是滤波器有效性条件。

滤波器有效性证明的继续

从上一部分推导继续，利用中心切片定理后得到：

关于 r_\parallel 的积分可以利用狄拉克函数的积分表示来处理。狄拉克函数可以写成：

因此对 r_\parallel 积分后得到 \delta(\mathbf{e}_\theta \cdot \boldsymbol{\nu})，上式变为：

狄拉克函数 \delta(\mathbf{e}_\theta \cdot \boldsymbol{\nu}) 的作用是选出所有满足 \mathbf{e}_\theta \cdot \boldsymbol{\nu} = 0 的方向，即 \boldsymbol{\nu} 垂直于 \mathbf{e}_\theta 的方向。此时 \boldsymbol{\nu}_\perp = \boldsymbol{\nu}，因此 \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp) = \mathcal{F}_3\{f\}(\boldsymbol{\nu}) 可以提到积分号外。

由于这个等式必须对任意函数 f(\mathbf{r}) 成立，剩余的积分必须恒等于 1，即：

这就是三维FBP滤波器必须满足的有效性条件。

科尔舍滤波器

在满足有效性条件的无穷多个滤波器中，需要选择具有最优噪声特性的滤波器。一类重要的滤波器是可分解滤波器族，其形式为：

这里滤波器被分解为两部分：H'_\Omega(\boldsymbol{\nu}) 仅依赖于频率向量，w(\mathbf{e}_\theta) 是方向依赖的权重函数。

将可分解形式代入有效性条件，可以解出：

第二个等式利用了狄拉克函数的齐次性质。

当选择均匀权重 w(\mathbf{e}_\theta) = 1 时，得到的滤波器能够最小化重建图像的方差，这就是科尔舍滤波器（Colsher's filter）：

分母是孔径 \Omega 在单位球面上与 \boldsymbol{\nu} 方向垂直的大圆的交线长度。图示展示了对于圆柱形PET孔径 \Omega_{\phi_{\max}}，给定频率方向 \boldsymbol{\nu} 时，分母对应的几何量是球带边界与垂直于 \boldsymbol{\nu} 的大圆的交集弧长。

科尔舍滤波器的显式形式

对于圆柱形PET孔径 \Omega_{\phi_{\max}} = \{\mathbf{e}_\theta | |\phi| < \phi_{\max}\}，科尔舍滤波器可以写成显式形式。设 \gamma 是频率向量 \boldsymbol{\nu} 与 xy 平面（即 \nu_z = 0 平面）的夹角，则：

当频率向量的倾斜角 |\gamma| 较小（即频率接近 xy 平面）时，\sin|\gamma| < \sin\phi_{\max}，此时垂直于 \boldsymbol{\nu} 的大圆完全落在孔径内，分母为 2\pi（完整大圆的周长），滤波器简化为 |\boldsymbol{\nu}|/(2\pi)，这与二维FBP的斜坡滤波器形式一致。

当频率向量倾斜角较大时，垂直于 \boldsymbol{\nu} 的大圆只有部分落在孔径内，分母变为孔径内的弧长，滤波器响应相应增大以补偿数据覆盖的减少。

科尔舍滤波器与二维FBP的斜坡滤波器具有相似的形式和噪声特性，是三维PET重建的标准滤波器。在实际应用中，还需要在投影平面的两个正交方向 \mathbf{e}_{\theta_1}^\perp 和 \mathbf{e}_{\theta_2}^\perp 上分别进行切趾处理，以控制噪声放大。

数据截断问题

实际的圆柱形PET探测器存在数据截断问题，这使得理论分析中的平移不变性假设不再成立。

圆柱形探测器的孔径 \Omega 并非平移不变的。对于倾斜的投影线（\phi \neq 0），投影数据 g(\psi, \phi, \mathbf{s} \cdot \mathbf{e}_{\theta_1}^\perp, \mathbf{s} \cdot \mathbf{e}_{\theta_2}^\perp) 沿着 \mathbf{e}_{\theta_2}^\perp 方向（即轴向方向在投影平面上的分量）会被截断。这是因为圆柱形探测器的轴向覆盖有限，倾斜的符合响应线在到达探测器边界时就无法被测量。

图示展示了这种截断效应。圆柱形探测器由多个环组成，相邻环之间可以形成符合事件。但对于较大的极角 \phi，只有轴向位置靠近探测器中心的符合线才能被完整测量，边缘区域的数据会被截断。这种截断会在重建图像中引入伪影，需要通过特殊的处理方法来缓解，如数据外推或迭代重建方法。

三维重投影算法

三维重投影算法（3D Reprojection Algorithm，3DRP）是解决数据截断问题的标准方法。该算法的核心思想是利用非倾斜投影（\phi = 0）的完备性来估计缺失的倾斜投影数据。

非倾斜投影 g(\psi, \phi=0, \mathbf{s}) 是平移不变的，并且满足奥尔洛夫条件，因此仅用这部分数据就可以完成重建。3DRP算法利用这一性质，先用非倾斜数据进行初步重建，再通过正投影来估计缺失的倾斜投影。

3DRP算法的具体步骤如下。第一步，使用非倾斜投影数据进行二维滤波反投影重建，得到初始的三维图像估计。第二步，对初始重建图像进行正投影，生成缺失的倾斜投影数据的估计值。第三步，将估计的投影数据与实际测量的投影数据进行融合，在有测量数据的区域使用测量值，在截断区域使用估计值。第四步，使用完整的（融合后的）数据进行三维滤波反投影重建，得到最终图像。

3DRP是三维PET解析重建的标准算法，在许多商业PET系统中得到应用。

三维到二维重排算法

三维数据采集相比二维模式具有更高的信噪比，因为可以接受更多的符合事件。然而二维重建算法计算速度快且对于许多应用已经足够。重排算法（Rebinning Algorithm）试图结合两者的优点：以三维模式采集数据以获得高灵敏度，然后将数据重排为一组二维正弦图，再使用快速的二维重建技术。

重排过程将三维采集的数据（参数化为 \psi、\phi、s_1、s_2）转换为一叠二维正弦图（参数化为 \psi、s），每个正弦图对应一个横断面切片。常用的重排方法包括单切片重排（Single-Slice Rebinning，SSRB）和傅里叶重排（Fourier Rebinning，FORE）。

图示展示了这个转换过程：上方是三维数据的表示，由多个倾斜的正弦图组成，形成一个四维数据结构；下方是重排后的二维正弦图堆栈，每个正弦图对应一个轴向位置 z。重排后可以对每个二维正弦图独立进行快速的二维FBP或迭代重建。

解析重建的数据校正

解析重建算法（如3DRP和FORE）都基于线积分模型，该模型假设测量数据是放射性分布沿投影线的线积分。然而原始PET数据受到多种物理效应的影响，在进行解析重建之前必须进行校正。

需要校正的效应包括：随机符合（两个不相关的光子偶然同时到达两个探测器形成的假符合）、散射符合（光子在体内散射后被检测到的符合，其响应线不再通过湮灭点）、衰减（光子在到达探测器前被组织吸收）、以及探测效率（不同探测器对的灵敏度差异）。只有在完成这些校正后，数据才近似满足线积分模型，才能应用解析重建算法。

三维PET解析重建方法总结

3DRP是三维PET解析重建的标准方法，提供了高质量的重建结果。但由于需要处理完整的三维数据和进行重投影步骤，计算量较大。

FAVOR算法是另一种三维解析算法，它使用的滤波反投影滤波器不要求完整的倾斜投影数据，因此是平移变化的（shift variant），但可以直接处理截断数据而无需重投影步骤。

FORE算法对三维PET产生了重大影响，因为它允许在三维采集数据上使用快速的二维迭代重建算法。在许多系统中，3DRP已被FORE+2D-FBP或FORE+2D-DIFT的组合所取代，这种方法在保持图像质量的同时大大减少了计算时间。

此外，还存在基于傅里叶域的精确重排方程，以及针对飞行时间（Time-of-Flight，TOF）PET数据的近似和精确傅里叶重排方法，这些都是活跃的研究方向。

附录：n维Radon变换

Radon变换是另一种描述断层成像的数学工具，它与平行射线变换（X射线变换）密切相关但有所不同。

对于函数 f: \mathbb{R}^n \to \mathbb{R}，其 n 维Radon变换 \mathbf{R}f 是定义在 \mathbb{R}^n 的超平面流形上的函数。设超平面 \Pi 由其法向量 \mathbf{e}_\Pi \in \mathbb{S}^{n-1} 和沿法向方向的有符号距离 l 确定，则Radon变换定义为：

这个积分是在满足 \mathbf{r} \cdot \mathbf{e}_\Pi = l 的超平面上进行的面积分。在二维情况下，超平面退化为直线；在三维情况下，超平面是通常意义上的平面。

平行射线变换与Radon变换的关系

平行射线变换（X射线变换）沿直线积分，而Radon变换沿超平面积分。两者之间存在明确的数学关系。

对于所有与超平面法向量 \mathbf{e}_\Pi \in \mathbb{S}^{n-1} 正交的射线方向 \mathbf{e}_\theta \in \mathbb{S}^{n-1}，Radon变换可以表示为平行射线变换在相应超平面上的积分：

在二维情况（n = 2）下，Radon变换与平行射线变换是等价的，只是参数化方式不同。超平面（直线）的法向量 \mathbf{e}_\Pi 与射线方向 \mathbf{e}_\theta 正交，\mathbf{e}_\Pi 对应于 \mathbf{e}_{\theta}^\perp。因此有：

其中 l 和 s 表示同一个径向坐标，只是参考方向不同。图示展示了二维情况下两种变换的几何关系：符合响应线 \text{LOR}(\mathbf{e}_\theta, s) 可以等价地用其法向方向 \mathbf{e}_\Pi 和到原点的距离 l 来描述。

在三维及更高维情况下，两种变换不再等价，因为Radon变换在平面上积分而X射线变换在直线上积分，但它们仍然通过上述积分关系相联系。

附录：三维PET中的正弦图参数化

在三维PET系统中，符合响应线（LOR）可以用多种方式进行参数化。一种实用的方式是根据两个探测器的轴向位置来索引。设两个探测器的轴向坐标分别为 z_A 和 z_B，则定义轴向位置 z = z_A + z_B（反映LOR的轴向中心位置）和环差 \Delta = z_A - z_B（反映LOR的倾斜程度）。

使用这种参数化，三维正弦图数据可以表示为 m(s, \psi, z, \Delta)，其中 s 是径向位置，\psi 是方位角。这种表示与之前使用的 g(\psi, \phi, s, s_2) 表示之间存在确定的几何关系：

其中 R 是探测器环的半径。这个公式将环差 \Delta 和径向位置 s 转换为极角 \phi，将轴向位置 z 转换为投影平面内的轴向坐标 s_2。

图示展示了三维正弦图的数据结构。对于固定的轴向位置 z 和环差 \Delta，数据形成一个二维正弦图（s vs \psi）。不同的 z 和 \Delta 组合构成一系列这样的正弦图，形成四维数据结构。当 \Delta = 0 时对应非倾斜投影（\phi = 0），\Delta \neq 0 对应倾斜投影。

附录：单切片重排算法

单切片重排（Single Slice Rebinning，SSRB）是最简单的三维到二维重排方法，它是一种近似算法。

SSRB的基本思想是将所有经过某个轴向位置 z 的倾斜LOR都近似地分配到该位置的二维正弦图中。具体来说，对于给定的 (s, \psi, z)，SSRB正弦图由以下公式定义：

这个公式对所有可能的环差 \Delta（从 -\Delta_{\max}(z) 到 +\Delta_{\max}(z)）进行平均，其中 \Delta_{\max}(z) 取决于探测器的几何约束。

图示展示了SSRB的几何含义：所有以某个轴向位置为中心、但具有不同倾斜角的LOR被平均后分配到单一的二维切片。这种近似在LOR倾斜角较小时是合理的，但对于较大的倾斜角会引入误差，导致轴向分辨率下降。

SSRB的优点是计算简单快速，缺点是重建图像的轴向分辨率会随着环差的增大而降低。

附录：傅里叶重排算法

傅里叶重排（Fourier Rebinning，FORE）是一种在傅里叶域进行的近似重排算法，比SSRB具有更高的精度。

FORE的数学基础是正弦图的二维傅里叶变换。对于三维正弦图 m(s, \psi, z, \Delta)，关于 (s, \psi) 的二维傅里叶变换定义为：

其中 \nu_s 是对应于径向变量 s 的频率，k 是对应于角度变量 \psi 的整数频率（角谐波阶数）。

FORE算法的核心是频率-距离关系（frequency-distance relation）。这个关系指出，傅里叶变换 \mathcal{F}_2\{m\} 在频率点 (\nu_s, k) 处的值主要来自于空间中某个特定位置 S 的贡献，该位置由 t = -k/(2\pi\nu_s) 确定。利用这个性质，可以建立倾斜正弦图与非倾斜正弦图之间的近似关系：

这个近似式表明，倾斜正弦图（\Delta \neq 0）在频率点 (\nu_s, k) 处的傅里叶系数，近似等于非倾斜正弦图（\Delta = 0）在经过轴向位移后的相应频率点的傅里叶系数。轴向位移量为 \frac{k\Delta}{2\pi\nu_s \cdot 2R}，它依赖于频率 (\nu_s, k) 和环差 \Delta。

FORE算法步骤

FORE算法的实现分为以下步骤。

第一步，初始化累加器 M(\nu_s, k, z) = 0，用于存储重排后的傅里叶系数。

第二步，对每个倾斜正弦图（由 z 和 \Delta 确定）进行处理：首先计算其二维傅里叶变换 \mathcal{F}_2\{m\}(\nu_s, k, z, \Delta)；然后对于每个频率点 (\nu_s, k)，计算校正后的轴向位置 z' = z - \frac{k\Delta}{2\pi\nu_s \cdot 2R}；最后将傅里叶系数 \mathcal{F}_2\{m\}(\nu_s, k, z, \Delta) 累加到 M(\nu_s, k, z') 中。

第三步，对累加器 M(\nu_s, k, z) 进行归一化处理，以补偿每个频率点可能接收到不同数量的贡献。

第四步，对归一化后的 M(\nu_s, k, z) 进行二维逆傅里叶变换，得到重排后的二维正弦图 m(s, \psi, z) = \mathcal{F}_2^{-1}\{M(\nu_s, k, z)\}。

第五步，对每个轴向位置 z 的二维正弦图 m(s, \psi, z) 进行独立的二维重建。

FORE算法在傅里叶域进行操作，能够根据频率内容精确地分配倾斜投影的贡献，因此比简单的空间域平均（SSRB）具有更好的轴向分辨率保持能力。FORE已成为三维PET数据处理的标准方法之一，特别是与快速的二维迭代重建算法结合使用时，能够在保持图像质量的同时显著减少计算时间。

迭代重建方法

从解析重建到迭代重建的过渡

前面详细介绍的解析重建方法以线积分模型为基础，通过滤波反投影（FBP）算子实现对采集模型的线性求逆。这种方法的优势在于其数学框架清晰，能够帮助我们理解断层重建问题的本质特性，包括 n 维X射线变换的不适定性、FBP重建中的伪影来源、斜坡滤波器对高频噪声的放大效应、数据截断带来的问题、以及三维情况下的数据冗余性等。

然而，解析重建的核心假设是数据采集过程可以用理想化的线积分模型来描述，且不考虑数据的统计特性。这种理想化假设与实际采集条件之间存在差距：真实的测量数据受到探测器响应、散射、衰减、随机符合等多种物理效应的影响，同时光子计数本身服从泊松分布，具有随机噪声。当实际数据不满足线积分模型的假设时，解析重建的图像质量会明显退化。

解析重建的局限性

解析重建方法基于以下理想化假设：采集过程严格遵循线积分模型、数据是完整的（没有缺失的投影角度或截断）、数据是连续的（无离散化误差）、数据是无噪声的。这些假设在实际成像中均无法完全满足。

实际的发射断层扫描数据具有以下特点：探测器具有有限的空间分辨率，会对真实信号产生模糊效应；由于探测器几何限制或患者运动，部分投影数据可能缺失；数据以离散的探测器单元和时间帧采集，而非连续分布；光子计数服从泊松统计，低计数区域的相对噪声尤其显著。

当重建数据不满足模型假设时，解析方法会产生系统性偏差和噪声放大，导致图像质量下降。解决这一问题的途径是采用更真实的采集模型，这正是迭代重建方法的出发点。

迭代重建的基本思想

迭代重建采用更真实的采集模型，该模型可以纳入探测器分辨率、数据缺失、离散采样、以及随机噪声等实际因素。由于模型的复杂性，通常无法像解析方法那样找到闭式的逆变换公式，因此需要通过迭代的方式逐步逼近解。

迭代重建的核心思想是：从一个初始图像估计出发，通过反复应用采集模型（正问题）来预测该估计图像会产生什么样的投影数据，将预测数据与实际测量数据进行比较，根据比较结果调整图像估计，如此循环往复，直到图像估计序列收敛到满意的解。

迭代重建的工作流程

迭代重建的完整流程可以描述如下。首先，体内真实的放射性示踪剂分布（这是我们想要重建的未知量）通过实际的采集系统产生测量数据，这些数据以正弦图或其他形式存储。

在重建过程中，我们维护一个对放射性分布的当前估计。将这个估计图像输入采集模型（也称为正问题模型或模拟器），该模型模拟采集过程并预测出对应的投影数据。这里的采集模型可以包含各种物理效应的建模，比理想的线积分模型更加精细。

接下来，将模型预测的投影数据与实际测量的投影数据进行比较。比较的方式取决于所采用的统计模型，例如可以计算两者之间的某种距离或似然函数。基于这个比较结果，定义一个目标函数（或称代价函数），该函数量化了当前估计与真实解之间的差距。

优化算法根据目标函数的值及其梯度信息，计算出对当前估计的修正量，更新图像估计。这个更新后的估计再次进入正问题模型，开始下一轮迭代。整个过程不断重复，直到目标函数收敛或达到预设的迭代次数。

迭代重建的关键组成部分包括：采集模型（描述从图像到数据的正向映射）、统计模型（描述数据的噪声特性）、目标函数（量化估计与数据的拟合程度，可能还包含正则化项）、以及优化算法（决定如何更新图像估计）。后续内容将逐一深入讨论这些组成部分。

迭代重建的核心特征

与解析重建通过离散化实现连续解析公式（如FBP）不同，迭代重建从一开始就在离散框架下工作。这种方法依赖于对测量数据和重建图像的离散表示，将原本的连续逆问题转化为有限维的优化问题。

迭代重建的三个核心优势在于：首先，它天然地处理离散数据和离散图像，无需担心连续到离散转换带来的近似误差；其次，它允许在重建框架中纳入测量数据的统计模型，从而能够正确处理泊松噪声等随机性；第三，它能够引入关于图像的先验知识，例如图像的平滑性、稀疏性或解剖结构约束，这些先验信息可以有效地正则化不适定问题，提高重建质量。

迭代重建方法的分类

根据所包含的建模要素，迭代重建方法有不同的名称。代数重建方法是最基本的形式，它仅依赖于测量数据和重建图像的离散表示，将重建问题表述为求解线性方程组，不涉及数据的统计特性。

统计重建方法在代数方法的基础上增加了测量数据的统计模型。由于光子计数服从泊松分布，统计方法通过最大化似然函数来估计图像，能够更好地处理低计数情况下的噪声。

贝叶斯重建方法是最完整的框架，它不仅包含离散表示和统计模型，还引入了图像的先验概率分布。先验分布编码了我们对图像特性的预期，例如相邻像素值应该相似（平滑先验）或图像在某个变换域下应该稀疏（稀疏先验）。贝叶斯方法通过最大化后验概率来估计图像，实现了数据拟合与先验约束之间的平衡。

离散模型的建立

在解析重建中，图像被视为连续函数 f(\mathbf{r})，正弦图也是连续函数 g(\mathbf{e}_\theta, \mathbf{s})。迭代重建则将它们转换为有限维向量：图像表示为体素值组成的向量 \boldsymbol{\theta}，正弦图表示为探测器测量值组成的向量 \mathbf{m}。

以一个 16 \times 16 的二维图像为例，将其按行或列展开后得到一个包含 256 个元素的一维向量 \boldsymbol{\theta}。向量的每个分量 \theta_j 对应图像中一个体素的值。这种向量化处理使得图像重建问题可以用标准的线性代数和优化理论来分析。

迭代重建的五个构成要素

完整的迭代重建方法由以下五个要素组成。

第一个要素是图像的有限参数化。图像被表示为参数向量 \boldsymbol{\theta} = \{\theta_j | j = 1, \cdots, J\}，其中 J 是参数的总数（例如体素的数目）。

第二个要素是离散测量数据均值的模型。测量数据向量 \mathbf{m} = \{m_i \in \mathbb{N} | i = 1, \cdots, I\} 的期望值 q_i 由以下线性模型给出：

这个公式描述了第 i 个探测器单元的期望计数如何依赖于所有体素的放射性浓度。系数 a_{ij} > 0 表示位于第 j 个体素处发生的衰变事件被第 i 条响应线检测到的概率（或比例），这个系数包含了探测器的几何灵敏度和物理效应的建模。常数项 b_i 表示与放射性分布无关的背景计数，例如散射事件和随机符合。

第三个要素是测量数据的概率分布函数，描述了给定期望值 q_i 时实际测量值 m_i 的统计分布。

第四个要素是目标函数，它量化了图像估计与测量数据之间的一致性，可能还包含正则化项。

第五个要素是迭代算法，用于优化目标函数并求解图像估计。

图像参数化

图像参数化是迭代重建的第一步，它决定了如何用有限个参数来表示连续的放射性分布。一般形式是将图像表示为空间基函数的线性组合：

其中 f_j(x, y, z) 是第 j 个基函数，\theta_j 是对应的系数。不同的基函数选择对应不同的参数化方式。

体素基函数

最常用的参数化方式是体素表示。每个基函数 f_j 是一个以点 (x_j, y_j, z_j) 为中心的立方体指示函数：

其中 \Delta x、\Delta y、\Delta z 是体素在三个方向上的半宽度。这种表示假设放射性浓度在每个体素内部是均匀的，体素之间没有重叠。体素表示的优点是简单直观，与图像的存储格式直接对应；缺点是基函数不光滑，可能在体素边界处产生阶梯状伪影。

Blob基函数

另一种选择是使用具有球对称性的平滑重叠基函数，称为Blob。每个Blob基函数 f_j(\mathbf{r}) 是一个仅依赖于到中心点 \mathbf{r}_j 距离的函数：

其中 b(\cdot) 是一个径向函数，典型的选择是广义Kaiser-Bessel窗函数。Blob基函数彼此重叠，具有平滑的过渡，能够产生更光滑的重建图像。空间基函数的选择可以视为对图像的一种隐式约束：体素假设图像是分段常数的，而Blob假设图像具有一定的平滑性。

测量数据均值的模型

迭代重建的第二个核心要素是建立测量数据期望值的数学模型。这个模型描述了图像空间中的放射性分布如何映射到探测器测量的期望计数。

定义 \alpha_i(\mathbf{r}) 为位于空间位置 \mathbf{r} 处发生的衰变事件被第 i 条响应线（LoR）检测到的概率（或检测效率）。这个函数综合了探测器的几何灵敏度、立体角覆盖、以及可能的物理效应（如衰减）。另外定义 b_i 为第 i 条响应线的期望背景事件数，包括散射事件和随机符合事件，这些事件与真实的放射性分布无关。

第 i 条响应线的期望计数 q_i 可以写成对整个视场 \mathfrak{C} 的积分：

将图像的基函数展开式 f(\mathbf{r}) = \sum_{j=1}^{J} \theta_j f_j(\mathbf{r}) 代入上式：

括号内的积分定义为系统矩阵的元素：

这个系数 a_{ij} 的物理含义是：在第 j 个体素中发生的衰变事件被第 i 条响应线检测到的比例。它综合了体素的空间位置、探测器的几何关系、以及基函数 f_j(\mathbf{r}) 的形状。

由此，期望计数可以简洁地写成：

定义系统矩阵 \mathbf{A} = \{a_{ij} | i = 1, \cdots, I; j = 1, \cdots, J\}，其中 I 是响应线的总数，J 是体素的总数。将所有响应线的期望计数组成向量 \mathbf{q}，图像参数组成向量 \boldsymbol{\theta}，背景计数组成向量 \mathbf{b}，则整个模型可以用矩阵形式表示：

这里假设对于每个体素 j，至少存在一条响应线能够检测到它，即 \sum_{i=1}^{I} a_{ij} > 0 对所有 j 成立。这个条件保证了每个体素都对测量数据有贡献，从而原则上可以从数据中估计出该体素的值。

飞行时间PET模型

飞行时间（Time-of-Flight，ToF）PET是一种先进的成像技术，它利用两个湮灭光子到达探测器的时间差来估计湮灭事件在响应线上的位置。在传统PET中，我们只知道湮灭发生在响应线的某处，但不知道具体位置；在ToF-PET中，时间差信息提供了关于湮灭位置的额外约束。

设 h(\cdot) 为飞行时间响应函数，它描述了测量的时间差与真实湮灭位置之间的关系。ToF-PET的系统矩阵元素需要包含这个时间响应模型：

这个公式中，(\mathbf{e}_{\theta_i}, \mathbf{s}_i) 确定了第 i 条响应线的方向和位置，\delta_t 是测量到的时间差，c 是光速。表达式 (\mathbf{r} - \mathbf{s}_i) \cdot \mathbf{e}_{\theta_i} 计算了空间点 \mathbf{r} 沿响应线方向相对于参考点的位移，而 \delta_t c/2 是根据时间差推断的位移量（除以2是因为时间差对应两个光子路径差的一半）。

系数 a_{ij\delta_t} 的物理含义是：在第 j 个体素中发生的、且被第 i 条响应线检测到的、同时具有特定ToF值 \delta_t 的衰变事件的概率。ToF信息将原本均匀分布在整条响应线上的不确定性集中到一个较短的区间内，从而提高了重建的信噪比。

系统矩阵的几何解释

系统矩阵 \mathbf{A} 是一个 I \times J 的矩阵，其中 I 是响应线数目，J 是体素数目。理解系统矩阵的列和行的物理含义有助于深入把握迭代重建的本质。

列的含义

系统矩阵的第 j 列包含元素 \{a_{1j}, a_{2j}, \cdots, a_{Ij}\}，它描述了如果在第 j 个体素中放置一个单位强度的点源，各条响应线会检测到多少计数。换言之，第 j 列就是第 j 个体素的点扩散函数（Point Spread Function）在正弦图空间中的表示，也称为成像系统对该体素中点源的响应。

行的含义

系统矩阵的第 i 行包含元素 \{a_{i1}, a_{i2}, \cdots, a_{iJ}\}，它描述了各个体素对第 i 条响应线的贡献权重。将第 i 行重新排列成图像的空间形状，就得到了第 i 条响应线的灵敏度图像，显示了哪些体素位置的放射性会被该响应线检测到。

正投影与反投影的离散对应

解析重建和迭代重建中的正投影、反投影操作存在直接的对应关系。

在解析框架中，正投影（X射线变换）是一个连续算子，将图像函数 f(\mathbf{r}) 映射到投影数据：

在迭代框架中，正投影对应于矩阵-向量乘法：

这个乘法计算每条响应线的期望计数，本质上是将图像向量 \boldsymbol{\theta} 通过系统矩阵 \mathbf{A} 投影到数据空间。

解析框架中的反投影是正投影的伴随算子：

其中 \mathbf{X}^* 是 \mathbf{X} 的伴随。在迭代框架中，反投影对应于转置矩阵与向量的乘法：

这里 \mathbf{A}^\mathsf{T} 是系统矩阵 \mathbf{A} 的转置。正如连续情况下 \mathbf{X}^* 是 \mathbf{X} 的伴随，离散情况下 \mathbf{A}^\mathsf{T} 是 \mathbf{A} 的转置，两者在数学结构上完全对应。

理想化线积分系统矩阵

当采集模型简化为理想的线积分模型时，系统矩阵元素 a_{ij} 正比于第 i 条响应线 (\mathbf{e}_{\theta_i}, \mathbf{s}_i) 与第 j 个体素的交集长度：

这个公式中的狄拉克函数 \delta(\cdot) 选取出位于响应线上的点：(\mathbf{r} - \mathbf{s}_i) - ((\mathbf{r} - \mathbf{s}_i) \cdot \mathbf{e}_{\theta_i})\mathbf{e}_{\theta_i} 是点 \mathbf{r} 到响应线的垂直距离向量，当且仅当 \mathbf{r} 位于响应线上时该距离为零。

对于体素基函数，a_{ij} 就是响应线穿过体素的路径长度。如果响应线完全不经过某个体素，则对应的 a_{ij} = 0。这种稀疏性是系统矩阵的重要特征：虽然 \mathbf{A} 的维度可能非常大（I 和 J 都可以达到数百万），但绝大多数元素为零，因为每条响应线只穿过少数体素。

图示展示了两种几何配置下 a_{ij} 的计算：环形PET探测器（左）和平板探测器（右）。在两种情况下，a_{ij} 都等于响应线 i 穿过体素 j 的长度。Siddon算法是一种高效计算射线与体素交集长度的方法，广泛应用于系统矩阵的计算。

立体角系统矩阵

除了基于线积分的系统矩阵外，另一种更精确的建模方式是立体角系统矩阵。在这种模型中，系数 a_{ij} 正比于从第 j 个体素看向第 i 对探测器所张的立体角：

这种建模方式更准确地反映了探测器的几何灵敏度：距离探测器较近的体素具有更大的立体角覆盖，因此具有更高的检测效率。与简单的线积分模型相比，立体角模型考虑了探测器的有限尺寸效应。

然而，立体角系统矩阵在三维情况下没有解析的闭式表达式，需要通过数值积分或蒙特卡罗模拟来计算。这增加了计算复杂度，但能够提供更精确的物理建模。

探测器响应函数

实际的探测器模块并非理想的点检测器，而是具有有限的空间响应。每个探测器都有其固有的点扩散函数（Intrinsic detector PSF），描述了探测器对入射光子位置的响应不确定性。

当一个伽马光子入射到探测器阵列时，即使光子实际击中某个特定位置，由于闪烁晶体中的光传播和光电转换过程，探测到的位置信号会在空间上有一定的弥散。这种弥散效应使得真实的响应线（LoR）与记录的响应线之间存在偏差。

图示展示了PET探测器对中两个探测器阵列的固有PSF：响应不是一条理想的直线，而是在两端都有一定的空间展宽。这种展宽效应需要在系统矩阵中加以建模，才能获得准确的重建结果。

更真实的系统矩阵建模

为了获得高质量的重建图像，系统矩阵需要准确建模发射和检测过程中的各种物理效应。主要的物理效应包括：

相互作用深度（Depth of Interaction，DOI）效应源于伽马光子可能在探测器晶体的不同深度发生相互作用，而不是都在表面被吸收。对于斜入射的光子，不同的相互作用深度会导致不同的记录位置，从而产生视差误差。

探测器内散射指的是光子在一个探测器晶体中发生康普顿散射后，可能在相邻晶体中被吸收，导致记录的位置与真实入射位置不同。

非共线性（Acolinearity）是PET特有的效应。理论上，正电子湮灭产生的两个伽马光子应该严格反向传播（夹角180°），但由于湮灭时正电子-电子系统可能具有残余动量，实际夹角会有约0.5°的偏差。这个偏差导致记录的响应线不精确地通过湮灭点。

正电子射程（Positron Range）效应是另一个PET特有的限制因素。正电子在湮灭前会在组织中行进一段距离（取决于示踪剂核素的能量），因此湮灭点并不完全等于示踪剂的位置。

通过在系统矩阵中准确建模这些物理效应，可以在重建过程中补偿它们对图像质量的影响。

空间分辨率退化因素

PET成像的空间分辨率受到多种因素的综合限制。Moses和Derenzo提出了一个半经验公式来估算视场中心的空间分辨率（以半高全宽 \Gamma 表示）：

公式中各项的含义如下：d 是探测器晶体的宽度，d/2 项反映了探测器像素化造成的位置不确定性；D 是探测器环的直径，0.0022D 项来自非共线性效应（约0.5°的角度偏差在直径为 D 的探测器环上造成的位置误差）；r 是正电子射程的有效值；b 是探测器内散射和块效应等其他因素的贡献。

以临床全身PET为例，其空间分辨率约为4毫米，而MRI可以达到约1毫米。这种分辨率差异部分来自上述固有的物理限制。通过在系统矩阵 \mathbf{A} 中准确建模探测器响应函数，可以在重建过程中部分补偿这些退化效应，从而改善图像质量。

探测器响应的解析建模

系统矩阵 \mathbf{A}（维度为 I \times J）可以分解为多个组成部分，分别建模不同层次的物理效应：

矩阵 \mathbf{H}（维度为 J \times J）建模图像空间中的分辨率效应，这通常是一个各向异性、空间非平稳的三维卷积操作。它描述了由于正电子射程、非共线性等因素导致的图像空间模糊。

矩阵 \mathbf{X}（维度为 I \times J）是几何投影矩阵，如基于Siddon算法的线积分模型或立体角模型。如果忽略散射事件（散射通常建模在背景项 \mathbf{b} 中），\mathbf{X} 是一个稀疏矩阵。

矩阵 \mathbf{D}（维度为 I \times I）建模数据空间中的分辨率效应，这是一个各向异性、空间非平稳的四维卷积操作。它描述了探测器固有PSF、相互作用深度等因素导致的数据空间模糊。

完整的系统模型可以写成：

这个分解形式的一个关键点是：矩阵 \mathbf{H} 和 \mathbf{D} 的具体形式依赖于几何投影矩阵 \mathbf{X} 的选择。不同的 \mathbf{X} 定义方式会改变图像空间模糊和数据空间模糊之间的分配。

分辨率建模的效益

在系统矩阵中引入分辨率建模（特别是图像空间分辨率矩阵 \mathbf{H}）能够显著改善重建图像的质量。

图示展示了36次独立采集的重建结果对比。第一列使用简单的几何系统矩阵 \mathbf{A} = \mathbf{X}，重建图像 \hat{\boldsymbol{\theta}}_{A=X} 显示热点区域（高摄取区）的均值图像中热点形状较为清晰，但方差图像显示噪声水平较高。

第二列是对第一列结果进行后平滑处理，噪声明显降低但图像变得模糊，热点边缘不再清晰。

第三列使用包含分辨率建模的系统矩阵 \mathbf{A} = \mathbf{X}\mathbf{H}，重建图像 \hat{\boldsymbol{\theta}}_{A=XH} 在保持热点清晰度的同时，方差图像显示噪声水平与后平滑结果相当甚至更低。

实验表明，分辨率建模带来两个主要优势：首先，它改善了对比度恢复，减少了部分容积效应（即热点区域的摄取值向周围组织"溢出"以及周围组织的信号"溢入"热点区域的现象）；其次，它降低了体素方差，因为重建算法能够"拒绝"与物理模型不一致的噪声成分，只接受符合系统响应特性的信号。

分辨率建模的临床效益

分辨率建模在实际临床应用中的效益可以通过脑部PET成像来说明。图示展示了使用简单几何系统矩阵（\mathbf{A} = \mathbf{X}）和包含分辨率建模的系统矩阵（\mathbf{A} = \mathbf{X}\mathbf{H}）重建的脑部图像对比。

使用分辨率建模后，重建图像 \hat{\boldsymbol{\theta}}_{A=XH} 相比 \hat{\boldsymbol{\theta}}_{A=X} 具有更清晰的功能结构轮廓，脑区之间的边界更加分明。右侧的时间-活度曲线显示了不同脑区（尾状核Caudate、小脑Cerebellum、壳核Putamen、腹侧纹状体Ventral striatum）的放射性浓度随时间的变化。使用分辨率建模（实心方块）得到的曲线与不使用分辨率建模（空心方块）的结果存在系统性差异，这种差异直接影响后续的动力学参数估计。

分辨率建模带来的实际效益包括：更好的功能结构描绘，使得不同脑区的视觉区分更加清晰；以及对定量分析的影响，由于减少了部分容积效应，从感兴趣区域提取的生理参数更加准确。

分辨率建模的潜在问题

尽管分辨率建模能够改善图像质量，但它也可能引入新的伪影。

图示展示了一个仿真研究：左侧是真实的数字体模 \mathbf{f}，包含一个大圆盘和几个小圆点；中间是不使用分辨率建模的重建结果 \hat{\boldsymbol{\theta}}_{A=X}，图像整体较为模糊但没有明显伪影；右侧是使用精确分辨率建模的重建结果 \hat{\boldsymbol{\theta}}_{A=XH}，虽然整体更清晰，但在大圆盘的边缘出现了明显的边缘伪影（edge artifact）——边缘处出现了不应该存在的亮环。

这种边缘伪影的产生原因不仅仅是模型失配或迭代未收敛，更根本的原因在于问题的不适定性。分辨率建模本质上是在进行去卷积（deconvolution）操作，而去卷积是一个严重不适定的逆问题。复合算子 \mathbf{X}\mathbf{H} 的零空间（null space）是非平凡的，这意味着存在多个不同的图像 \boldsymbol{\theta} 产生完全相同的投影数据。边缘伪影与去卷积中的振铃效应（ringing）以及X射线变换的右奇异向量结构有关。

边缘伪影的缓解策略

缓解边缘伪影的一种方法是比较重建估计 \hat{\boldsymbol{\theta}}_{A=XH} 与经过平滑处理的体模 \mathbf{H}\mathbf{f}。

图示展示了不同平滑程度下的剖面线比较：蓝色实线是平滑后的体模 \mathbf{H}\mathbf{f}，红色虚线是重建估计 \hat{\boldsymbol{\theta}}_{A=XH}，黑色虚线是平滑后的估计 \mathbf{H}\hat{\boldsymbol{\theta}}_{A=XH}，紫色线是两者之差。

当 \mathbf{H} 的平滑核半高全宽（FWHM）为0毫米时（即没有分辨率模糊），重建估计与真实体模直接比较，边缘处存在显著差异（紫色线）。随着FWHM增加到4毫米、8毫米、12毫米，平滑后的估计 \mathbf{H}\hat{\boldsymbol{\theta}}_{A=XH} 与平滑后的体模 \mathbf{H}\mathbf{f} 之间的差异逐渐减小。

这个现象说明 \hat{\boldsymbol{\theta}}_{A=XH} 和真实体模 \mathbf{f} 都是方程 \mathbf{X}\mathbf{H}\boldsymbol{\theta} = \mathbf{g} 的解，因为 \mathbf{X}\mathbf{H} 具有非平凡的零空间。虽然两者在高分辨率下看起来不同，但经过 \mathbf{H} 平滑后它们变得一致，因为零空间中的高频成分被平滑操作消除了。

图示还展示了三种重建结果的直接比较：不使用分辨率建模的 \hat{\boldsymbol{\theta}}_{A=X}（左）、使用分辨率建模的 \hat{\boldsymbol{\theta}}_{A=XH}（中）、以及使用分辨率建模并后平滑的 \mathbf{H}\hat{\boldsymbol{\theta}}_{A=XH}（右）。\hat{\boldsymbol{\theta}}_{A=X} 和 \mathbf{H}\hat{\boldsymbol{\theta}}_{A=XH} 具有相同的空间分辨率，但后者的噪声方差更低。这等价于在新的（平滑的）空间基函数上进行重建，获得了更好的噪声特性而不损失有效分辨率。

分辨率建模总结

在图像空间中使用各向同性模型进行分辨率建模是一种近似方法，其核心要点如下：

首先，这是一个改进的模型，更接近真实的采集物理过程。即使是近似的分辨率模型，也比完全忽略分辨率效应要好。

其次，这种建模能够产生改进的结果，包括更准确的定量估计和更低的噪声水平。分辨率建模使得重建算法能够区分真实信号和系统响应不一致的噪声成分。

第三，引入分辨率建模意味着求解一个新的逆问题，这个新问题具有不同的稳定性条件。并非所有信息都能够稳定地恢复——高频成分由于去卷积的不适定性而难以可靠重建。在实际应用中，需要在分辨率恢复和稳定性之间取得平衡。

真实系统矩阵的估计方法

获取真实系统矩阵 \mathbf{A} 有三种主要途径。

第一种是实验测量方法：使用伪点源在扫描仪视场内的不同位置进行实际测量。每个点源位置对应系统矩阵的一列，记录该位置的点扩散函数。这种方法能够获得最真实的系统响应，但实验工作量巨大，且难以覆盖整个视场。

第二种是蒙特卡罗模拟方法：在计算机中模拟点源在视场内不同位置的发射和检测过程。通过详细模拟光子传输、探测器响应等物理过程，可以获得精确的系统矩阵估计。这种方法避免了实际实验的困难，但计算量非常大。

第三种是解析建模方法：基于对扫描仪响应函数的数学描述来计算系统矩阵。这种方法计算效率较高，但依赖于对物理过程的简化假设。

系统矩阵的规模是一个严峻的挑战：典型的临床PET系统中，\mathbf{A} 的维度约为 \mathcal{O}(100M \times 10M)，即约 1000 万亿个元素。即使 \mathbf{A} 是稀疏的，显式存储仍然不现实，因此实际中通常采用在线计算（on-the-fly）的方式：正投影器和反投影器在每次迭代时实时计算所需的矩阵-向量乘积，而不预先存储整个矩阵。

完整的系统矩阵模型

为了避免构建单一的巨型系统矩阵，可以将扫描仪响应分解为多个子矩阵，每个子矩阵建模不同的物理现象。除了之前介绍的分辨率矩阵 \mathbf{H} 和 \mathbf{D} 外，还需要考虑：

探测器效率矩阵 \mathbf{E}（维度 I \times I，对角矩阵）：不同的探测器对具有不同的检测效率，\mathbf{E} 对每条响应线的计数进行效率归一化校正。

衰减矩阵 \mathbf{L}（维度 I \times I，对角矩阵）：建模光子在体内的衰减。在SPECT中，衰减效应与发射位置相关，需要嵌入到几何投影矩阵 \mathbf{X} 中；在PET中，由于符合检测的特性，衰减因子只依赖于响应线而与发射位置无关，因此可以分离出来作为独立的对角矩阵。

运动矩阵 \mathbf{M}（维度 J \times J）：如果采集过程中患者存在运动（如呼吸运动、心脏跳动），\mathbf{M} 建模图像空间中的形变场，将参考帧的图像变换到各个运动状态。

综合所有这些因素，完整的系统模型为：

几何投影矩阵 \mathbf{X} 仍然是最大的组成部分，典型维度为 100M \times 10M = 1000T（万亿）个潜在元素。在实际实现中，\mathbf{X} 及其转置通过投影器（projector）和反投影器（backprojector）算法在线计算，而不是预先存储。这种分解方式的优势在于各个子矩阵可以独立优化和更新，且许多子矩阵（如 \mathbf{E}、\mathbf{L}）是对角矩阵，计算效率很高。

数据的概率分布函数

迭代重建的第三个核心要素是测量数据的概率分布函数。这个分布描述了在给定期望值的条件下，实际测量值的随机波动特性。正确的噪声模型是统计重建方法的基础。

噪声假设与泊松点过程

发射断层扫描中伽马光子的检测过程可以用泊松点过程来描述。这一结论基于以下物理假设：

第一，注入体内的放射性核素数目服从泊松分布。这是因为放射性药物的制备和注射过程涉及大量原子的随机行为。

第二，在任意时刻，放射性核素在体内的空间分布位置是相互独立的随机变量。这是示踪剂假设的核心内容：示踪剂的浓度足够低，使得示踪剂分子之间没有相互作用，各自独立地分布在组织中。

第三，放射性核素衰变发生的时刻是相互独立的随机变量，且服从参数为 \mu = \frac{T_{1/2}}{\ln 2} 的指数分布。这是放射性衰变的基本物理规律，T_{1/2} 是核素的半衰期。

第四，每次衰变的检测过程是相互独立的随机变量。这假设探测器是理想的，没有死时间效应（即探测器在检测到一个事件后立即准备好检测下一个事件）。

在这些假设下，可以证明各条响应线（LoR）上测量到的离散计数是相互独立的泊松随机变量。

泊松噪声模型

测量数据向量 \mathbf{m} 的每个分量 m_i 都是独立的随机变量。第 i 条响应线检测到的事件数 m_i 服从参数为 q_i 的泊松分布，记为 m_i \sim \text{Po}(q_i)。泊松分布的概率质量函数为：

其中 m_i \in \mathbb{N}_0（非负整数），q_i > 0 是期望计数。

泊松分布有一个独特的性质：其方差等于期望值，即 \text{Var}(m_i) = \mathbb{E}[m_i] = q_i。这意味着信噪比（Signal-to-Noise Ratio）为：

信噪比与期望计数的平方根成正比，因此计数越高，相对噪声越小。

图示展示了两种不同期望值下的泊松分布：当 q_i = 1.2 时（左图），分布高度偏斜，集中在小整数值；当 q_i = 12 时（右图），分布更加对称，接近高斯分布的形状。

在发射断层扫描中，典型的期望计数 q_i \ll 10，这意味着泊松分布与高斯分布有显著差异，不能用高斯近似 m_i \not\sim \mathcal{N}(q_i, q_i) 来代替。这是统计重建方法必须显式使用泊松模型的原因。

泊松噪声模型的应用条件

为了使泊松噪声模型成立，必须使用未经预处理的原始数据进行重建。这意味着所有的物理效应校正（如衰减校正、散射校正、随机符合校正等）都应该包含在系统模型 \mathbf{q} 中，而不是预先对数据 \mathbf{m} 进行校正。

如果在重建之前对数据进行预处理（例如减去散射估计、除以衰减因子等），处理后的数据将不再服从泊松分布，这会导致统计模型失配。模型失配会影响重建算法的最优性，可能导致偏差增加或噪声特性恶化。

因此，统计重建的标准做法是：保持原始数据不变，将所有校正因子纳入正向模型 \mathbf{q} = \mathbf{A}\boldsymbol{\theta} + \mathbf{b}，让重建算法在正确的统计框架下工作。

信噪比加权

解析重建和统计重建在处理数据时的一个本质区别在于对不同测量值的加权方式。

解析重建（如滤波反投影）对所有检测事件赋予相同的权重。无论某条响应线的计数是高还是低，它对最终重建图像的贡献方式是相同的。这种均匀加权没有利用数据的统计特性信息。

统计重建则根据信噪比对事件进行差异化加权：信噪比较高的测量值（即期望计数较大的响应线）获得更大的权重，因为这些测量值更可靠；信噪比较低的测量值获得较小的权重，因为它们受噪声影响更大，可靠性较低。

这种加权策略的数学基础来自最大似然估计或贝叶斯估计理论。通过正确利用泊松统计特性，统计重建能够在相同的数据条件下获得比解析重建更好的噪声-分辨率权衡，这也是统计方法在低计数成像中优势明显的原因。

目标函数

迭代重建的第四个核心要素是目标函数（也称代价函数或损失函数）。目标函数量化了图像估计与测量数据之间的一致性程度，优化算法的任务就是找到使目标函数取得最优值的图像。

重建问题可以形式化为以下优化问题：

这个优化问题的各个组成部分含义如下：

\boldsymbol{\theta} \in \mathbb{R}^J 是待估计的图像参数向量，包含 J 个体素的值。

\mathbf{m} \in \mathbb{N}^I 是测量数据向量，包含 I 条响应线的计数值；\mathbf{b} \in \mathbb{R}_+^I 是背景事件的期望值向量。

\Phi(\mathbf{m}, \cdot) 是数据拟合项（或数据保真项），它度量测量数据 \mathbf{m} 与模型预测 \mathbf{A}\boldsymbol{\theta} + \mathbf{b} 之间的差异。这个函数从 \mathbb{R}^I 映射到 \mathbb{R}，其具体形式通常由数据的统计模型推导得出。

U(\cdot) 是图像惩罚项（或正则化项），它编码了关于图像的先验信息或约束条件。这个函数从 \mathbb{R}^J 映射到 \mathbb{R}，用于抑制不期望的图像特性（如过度的噪声或不合理的结构）。

\beta \in \mathbb{R}_+ 是超参数，控制数据拟合项与惩罚项之间的相对权重。\beta 越大，正则化效果越强，图像越平滑但可能偏离数据；\beta 越小，越强调数据拟合，图像细节保留更多但噪声也更大。

\Theta 是解空间或搜索空间，定义了 \boldsymbol{\theta} 的允许取值范围。通过设计 \Theta，可以隐式地施加约束条件，例如非负性约束 \Theta = \mathbb{R}_+^J 要求所有体素值非负。

整个优化问题的目标是：在满足约束条件 \Theta 和先验假设 \beta U(\cdot) 的前提下，找到一个图像估计 \hat{\boldsymbol{\theta}}，使其对应的正向模型预测 \mathbf{A}\hat{\boldsymbol{\theta}} + \mathbf{b} 最好地拟合测量数据 \mathbf{m}。

目标函数的凸性与解的存在唯一性

定义总目标函数为 J(\boldsymbol{\theta}) = \Phi(\mathbf{m}, \mathbf{A}\boldsymbol{\theta} + \mathbf{b}) + \beta U(\boldsymbol{\theta})，则重建问题变为：

在一般情况下，最小值的存在性和唯一性并不能保证。为了确保优化问题有良好的性质，需要对 J 和 \Theta 施加一定的假设。基本假设包括：J 的定义域非空（\text{dom}J \neq \varnothing）、J 不恒为负无穷（-\infty \notin J(\Theta)）、以及定义域与约束集有交集（\text{dom}J \cap \Theta \neq \varnothing）。

关于解的存在性和唯一性有以下结论：如果 J 连续且 \Theta 紧致，则至少存在一个最小值点；如果 J 连续、强制（coercive，即当 \|\boldsymbol{\theta}\| \to \infty 时 J(\boldsymbol{\theta}) \to \infty）且 \Theta 闭合，则至少存在一个最小值点；如果 J 凸且 \Theta 闭合且凸，则任何局部最小值都是全局最小值，且最小值点的集合是凸集；进一步地，如果 J 下半连续且满足强制性或 \Theta 有界，则存在全局最小值；如果 J 严格凸，则全局最小值唯一。

选择惩罚项 U 和约束集 \Theta 的一个重要考虑因素就是确保目标函数的凸性，从而保证全局最小值的存在性。

数据拟合项：距离度量

数据拟合项 \Phi(\mathbf{m}, \mathbf{q}) 可以理解为测量数据 \mathbf{m} 与模型预测 \mathbf{q} = \mathbf{A}\boldsymbol{\theta} + \mathbf{b} 之间的某种距离或差异度量。不同的距离度量对应不同的噪声假设。

最小二乘距离

最小二乘（Least Squares，LS）距离适用于白高斯噪声假设，即各测量值的噪声独立同分布且方差相同：

这是最简单的距离形式，直接计算预测值与测量值之差的平方和。

加权最小二乘距离

加权最小二乘（Weighted Least Squares，WLS）距离，也称为马氏距离的平方，适用于一般的高斯噪声假设，即各测量值的噪声可以有不同的方差且可能相关：

其中 \boldsymbol{\Sigma} 是测量噪声的协方差矩阵。当噪声独立时，\boldsymbol{\Sigma} 是对角矩阵，\boldsymbol{\Sigma}^{-1} 的对角元素 1/\sigma_i^2 作为各测量值的权重。

广义Kullback-Leibler散度

广义Kullback-Leibler（KL）散度适用于泊松噪声假设，是发射断层扫描统计重建中最常用的数据拟合项：

这个表达式也称为与负熵相关的Bregman散度。当 m_i = 0 时，按照连续性约定 0 \log 0 = 0。KL散度总是非负的，且当且仅当 q_i = m_i 对所有 i 成立时取得最小值零。

数据拟合项：似然函数

数据拟合项的另一个重要来源是统计模型的似然函数。给定图像参数 \boldsymbol{\theta}，测量数据 \mathbf{m} 的概率分布为 p(\mathbf{m}|\boldsymbol{\theta})，这个分布被视为 \boldsymbol{\theta} 的函数时称为似然函数。

最大似然估计（Maximum Likelihood Estimation，MLE）的思想是：在所有可能的图像中，选择那个使观测到当前数据 \mathbf{m} 的概率最大的图像。形式上：

其中对数似然函数 L(\boldsymbol{\theta}, \mathbf{m}) = \log p(\mathbf{m}|\boldsymbol{\theta})。由于对数函数是单调递增的，最大化对数似然等价于最大化似然本身，但对数形式在数值计算上更稳定（避免了概率连乘导致的下溢）。

最大化对数似然等价于最小化负对数似然 -L(\boldsymbol{\theta}, \mathbf{m})，因此可以将负对数似然作为数据拟合项：\Phi(\mathbf{m}, \mathbf{q}) = -L(\boldsymbol{\theta}, \mathbf{m})。当不使用正则化项（U = 0）时，最小化目标函数就等价于最大似然估计。

泊松对数似然

对于泊松噪声模型，各响应线的测量值独立且 m_i \sim \text{Po}(q_i)，其中 q_i = \sum_{j=1}^{J} a_{ij}\theta_j + b_i。联合概率为各边缘概率的乘积：

取对数后得到对数似然函数：

最后一项 \log m_i! 不依赖于 \boldsymbol{\theta}，在优化时可以忽略。因此，最大似然估计可以写成：

这个目标函数关于 \boldsymbol{\theta} 是凹函数（因为 \log 是凹函数，线性组合保持凹性），因此最大化问题是凸优化问题，等价于最小化其负值（凸函数）。泊松对数似然的负值与广义KL散度只相差一个不依赖于 \boldsymbol{\theta} 的常数，因此两者在优化意义上等价。

对数似然函数的可微性与凹性

泊松对数似然函数为：

这个函数关于 \boldsymbol{\theta} 是可微的，可以计算其Hessian矩阵（二阶偏导数矩阵）。Hessian矩阵的元素为：

为了判断函数的凹凸性，需要检验Hessian矩阵是否半负定。对于任意向量 \boldsymbol{\nu}，计算二次型：

由于 m_i \geq 0、期望值 q_i = \sum_j a_{ij}\theta_j + b_i > 0，以及平方项 (\sum_k a_{ik}\nu_k)^2 \geq 0，整个表达式总是非正的。这证明了Hessian矩阵是半负定的，因此 L(\boldsymbol{\theta}, \mathbf{m}) 是凹函数。

相应地，负对数似然 -L(\boldsymbol{\theta}, \mathbf{m}) 是可微的、凸的，并且在非负锥（\theta_j \geq 0）上是强制的（coercive）。由于非负锥是闭集，根据凸优化理论，L(\boldsymbol{\theta}, \mathbf{m}) 在非负约束下存在全局最大值。

最大似然与KL散度的等价性

最大似然估计和最小化KL散度实际上给出相同的解。最大似然估计为：

最小化KL散度的估计为：

将KL散度展开：m_i \log(m_i/q_i) + q_i - m_i = m_i \log m_i - m_i \log q_i + q_i - m_i。其中 m_i \log m_i 和 -m_i 都是常数（不依赖于 \boldsymbol{\theta}），因此最小化KL散度等价于最小化 -m_i \log q_i + q_i，即等价于最大化 m_i \log q_i - q_i，这正是对数似然的主要部分。因此：

两种方法得到相同的解集。

最大后验估计

最大似然估计将图像 \boldsymbol{\theta} 视为确定性的未知参数，不引入关于图像的先验信息：

在贝叶斯框架中，图像本身也被视为随机变量，具有先验概率分布 p(\boldsymbol{\theta})。先验分布编码了我们在观测数据之前对图像的信念或期望，例如图像应该是平滑的或具有某种结构特性。

根据贝叶斯定理，给定观测数据 \mathbf{m} 后，图像的后验概率分布为：

其中 p(\mathbf{m}|\boldsymbol{\theta}) 是似然函数，p(\boldsymbol{\theta}) 是先验分布，p(\mathbf{m}) 是证据（边缘似然），起归一化常数的作用。

取对数后：

由于 \log p(\mathbf{m}) 不依赖于 \boldsymbol{\theta}，最大后验（Maximum A Posteriori，MAP）估计为：

MAP估计在最大似然的基础上增加了先验项 \log p(\boldsymbol{\theta})，当先验取均匀分布时（p(\boldsymbol{\theta}) = \text{常数}），MAP退化为ML。

贝叶斯重建的不同层次

MAP估计虽然利用了先验信息，但其结果仍然是一个点估计 \hat{\boldsymbol{\theta}}_{\text{MAP}}，只给出后验分布的众数（最大值点），不提供关于估计不确定性的信息。

完全贝叶斯框架的目标更为宏大：不仅求后验分布的最大值点，而是估计整个后验分布 p(\boldsymbol{\theta}|\mathbf{m})。有了完整的后验分布，可以计算各种统计量，例如后验均值作为点估计、可信区间量化不确定性、后验方差或协方差描述估计的精度等。

完全贝叶斯推断在计算上比MAP困难得多，因为需要对高维空间进行积分或采样。常用的方法包括马尔可夫链蒙特卡罗（MCMC）采样和变分推断。这些方法在医学成像中的应用是活跃的研究方向。

简单先验模型：体素独立假设

最简单的先验模型假设各体素的值相互独立。

独立高斯模型

假设每个体素值独立地服从高斯分布，均值为 \bar{\theta}_j，方差为 \sigma_j^2：

其中 Z 是归一化常数。这个先验假设图像在期望图像 \bar{\boldsymbol{\theta}} 附近波动，波动幅度由方差控制。

最大熵模型

另一种独立先验是最大熵模型：

这个先验来源于信息论，当仅知道图像的期望值时，最大熵原理给出这种形式的先验。

独立先验模型的优点是优化问题可以分解，计算相对简单。但它们的主要缺点是需要预先指定期望图像 \bar{\boldsymbol{\theta}}，而这在实际应用中往往是未知的。更根本的问题是，这类模型假设体素值围绕固定均值独立波动，且交换任意两个体素的残差不改变概率——这与真实医学图像的结构特性严重不符。解剖和功能图像具有复杂的空间相关性和结构信息，独立先验无法捕捉这些特性，因此在发射断层重建中很少使用。

具有空间相关性的先验：马尔可夫随机场

为了建模体素之间的空间相关性，可以使用马尔可夫随机场（Markov Random Fields，MRF）。在MRF框架中，图像被表示为一个无向图：每个体素对应图中的一个顶点，相邻体素之间通过边连接，表示它们的值具有统计相关性。

图中的团（clique）c 是一个全连通子图，即子图中任意两个顶点之间都有边相连。在图像处理中，常用的团包括：单个体素（一阶团）、相邻的两个体素（二阶团，即边）、以及更大的邻域结构。\mathbb{C} 表示所有团的集合。

MRF的概率分布由Gibbs分布给出：

其中 U(\boldsymbol{\theta}) 称为Gibbs能量函数，它是所有团上势函数 \psi(c) 的总和。Z 是配分函数，保证概率归一化。\beta 是逆温度参数，控制分布的集中程度。

Gibbs分布的形式意味着：能量越低的图像配置概率越高。通过设计势函数 \psi(c)，可以编码对图像的各种期望，例如相邻体素值应该相似（平滑先验）或图像应该具有分段常数的特性（边缘保持先验）。

取负对数后，-\log p(\boldsymbol{\theta}) = \beta U(\boldsymbol{\theta}) + \log Z，由于 \log Z 不依赖于 \boldsymbol{\theta}，MAP估计中的先验项 \log p(\boldsymbol{\theta}) 正比于 -U(\boldsymbol{\theta})。因此，目标函数中的惩罚项 \beta U(\boldsymbol{\theta}) 可以直接解释为Gibbs能量，这建立了优化观点（正则化）和统计观点（贝叶斯推断）之间的联系。

马尔可夫随机场的MAP估计

当先验分布采用马尔可夫随机场模型时，Gibbs分布为：

MAP估计变为：

对于基于成对团（二阶团）的MRF模型，能量函数可以写成所有相邻体素对上势函数的加权和：

其中 N_j 表示体素 j 的三维邻域（通常取6邻域、18邻域或26邻域），条件 k > j 避免重复计算同一对体素，\omega_{kj} 是邻域权重（通常与体素间距离相关，距离近的权重大），\psi(\theta_k, \theta_j) 是定义在相邻体素对上的势函数。

如果目标函数 -L(\boldsymbol{\theta}, \mathbf{m}) + \beta U(\boldsymbol{\theta}) 是凸函数、下半连续、强制的，且约束集 \Theta 是闭凸集，则存在全局最大值，可以使用标准的非线性优化方法求解。

势函数的典型形式

势函数 \psi(\theta_k, \theta_j) 的设计决定了先验模型的性质。其基本作用是惩罚相邻体素之间的差异，从而促使图像具有某种平滑性或规则性。

二次势函数

二次势函数对相邻体素的差异施加平方惩罚，\sigma 控制惩罚的强度。这是最简单的形式，对应于高斯马尔可夫随机场。其优点是数学上易于处理（目标函数保持严格凸性），但缺点是对大差异的惩罚过强，会过度平滑图像边缘。

全变分势函数

全变分（Total Variation，TV）势函数使用绝对值惩罚。与二次惩罚相比，TV对大差异的惩罚增长较慢（线性而非二次），因此允许相邻体素之间存在较大的跳变，从而能够保持图像边缘。TV正则化倾向于产生分段常数的图像，这在某些应用中是期望的特性。TV势函数在原点不可微，需要使用次梯度方法或近端算法进行优化。

Huber势函数

Huber势函数是二次函数和绝对值函数的混合：当差异较小时（|\theta_k - \theta_j| \leq \delta）采用二次惩罚，当差异较大时（|\theta_k - \theta_j| > \delta）采用线性惩罚。参数 \delta 是过渡阈值。Huber势函数结合了二次势函数在小差异处的平滑性（保证可微）和TV在大差异处的边缘保持能力。

相对差异势函数

相对差异（Relative Difference，RD）势函数的惩罚强度不仅取决于差异的绝对大小，还取决于体素值本身的大小。当 \theta_k 和 \theta_j 都较大时，分母较大，惩罚相对较轻；当体素值较小时，惩罚相对较重。这种自适应特性在发射断层成像中有一定意义，因为高摄取区域的绝对噪声通常较大。\gamma 是调节参数。这个势函数是非对称的（交换 \theta_k 和 \theta_j 后形式改变）。

图示比较了这四种势函数的形状：横轴是 \theta_k - \theta_j，纵轴是势函数值。二次势函数（蓝色）增长最快，形成尖锐的抛物线；TV（橙色）是V形，增长为线性；Huber（绿色）在原点附近是抛物线，远离原点后变为线性；RD（红色/紫色/棕色）的形状依赖于参数设置和体素值本身。

这些势函数都是数学上方便的选择，在保证优化问题有良好性质的同时实现平滑和边缘保持的平衡。但它们主要是启发式设计，并非直接来源于对真实医学图像统计特性的研究。

迭代算法

至此，迭代重建的前四个要素已经确定：图像的参数化表示（体素或blob）、采集数据均值的模型（系统矩阵 \mathbf{A} 和背景 \mathbf{b}）、数据的统计模型（泊松分布）、以及目标函数（ML或MAP）。

剩下的最后一个要素是设计优化算法来最小化（凸）或最大化（凹）目标函数。这是迭代重建中计算上最具挑战性的部分，因为问题的规模巨大（数百万个未知量）且需要处理各种约束条件。

不完全数据与完全数据

泊松对数似然函数的形式为：

这个似然函数的一个关键特点是：每个测量值 m_i 的概率 p(m_i|\boldsymbol{\theta}) 依赖于多个图像参数 \theta_j（通过求和 \sum_j a_{ij}\theta_j）。由于对数函数内部有求和，无法将对数似然分解为各参数的独立贡献，这使得直接解析最大化变得困难。在统计学术语中，这种数据称为不完全数据（incomplete data）。

为了推导高效的优化算法，引入完全数据（complete data）的概念。假设我们不仅知道每条响应线的总计数 m_i，还知道这些计数的来源分解：设 n_{ij} 为响应线 i 检测到的、来源于体素 j 的事件数，n_i 为来源于背景的事件数。则有：

这是一个多对一的映射：完全数据 \mathbf{n} = \{n_{ij}, n_i\} 可以唯一确定不完全数据 \mathbf{m}，但反过来不成立（给定 m_i 无法唯一确定各 n_{ij}）。

完全数据的对数似然函数为：

这个表达式的关键优势在于：对数函数内部不再有求和，各 \theta_j 的贡献是分离的。这使得完全数据似然函数可以分解为各参数的独立项，极大地简化了最大化过程。期望最大化（EM）算法正是利用这一性质，通过在不完全数据和完全数据之间迭代来求解最大似然估计。

期望最大化算法

完全数据的对数似然函数为：

这个函数相对于 \boldsymbol{\theta} 很容易最大化，因为各 \theta_j 的贡献是分离的。然而问题在于完全数据 \mathbf{n} = \{n_{ij}, n_i\} 是未知的——我们只能观测到不完全数据 \mathbf{m}。

期望最大化（Expectation Maximization，EM）算法提供了一种优雅的解决方案：不直接最大化完全数据似然，而是最大化其在当前估计下的条件期望。算法由两个交替执行的步骤组成。

E步（期望步）

给定当前的参数估计 \boldsymbol{\theta}^{(p)}，计算完全数据对数似然关于完全数据的条件期望：

这个期望是关于完全数据 \mathbf{n} 的，条件是观测到的不完全数据 \mathbf{m} 和当前参数估计 \boldsymbol{\theta}^{(p)}。

M步（最大化步）

最大化 Q 函数得到新的参数估计：

将 \boldsymbol{\theta}^{(p+1)} 作为新的当前估计，重复E步和M步，直到收敛。

EM算法的核心性质是单调性：每次迭代后，不完全数据的对数似然不会减少，即 L(\boldsymbol{\theta}^{(p+1)}, \mathbf{m}) \geq L(\boldsymbol{\theta}^{(p)}, \mathbf{m})。

MLEM单调性证明

为了证明EM算法的单调性，首先建立 Q 函数与不完全数据对数似然之间的关系。

Q 函数可以写成：

利用条件概率的关系 p(\mathbf{n}|\boldsymbol{\theta}) = p(\mathbf{n}|\mathbf{m}, \boldsymbol{\theta}) p(\mathbf{m}|\boldsymbol{\theta})，上式变为：

展开对数并分离积分：

因此：Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)}) = H(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)}) + L(\boldsymbol{\theta})。

根据M步的定义，\boldsymbol{\theta}^{(p+1)} 最大化 Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)})，因此对于任意 \boldsymbol{\theta}：

特别地，取 \boldsymbol{\theta} = \boldsymbol{\theta}^{(p)}：

由此得到：

根据Jensen不等式，H(\boldsymbol{\theta}^{(p)}|\boldsymbol{\theta}^{(p)}) \geq H(\boldsymbol{\theta}^{(p+1)}|\boldsymbol{\theta}^{(p)})，因此右边非负，从而 L(\boldsymbol{\theta}^{(p+1)}) \geq L(\boldsymbol{\theta}^{(p)})。

ML-EM算法的E步

在发射断层成像的泊松模型下，给定总计数 m_i，各来源的分解 \{n_{ij}, n_i\} 服从多项分布。利用当前估计 \boldsymbol{\theta}^{(p)}，可以计算各分量的条件期望：

这些公式的物理含义是直观的：总计数 m_i 按照各来源对期望计数的贡献比例进行分配。体素 j 对响应线 i 的贡献为 a_{ij}\theta_j^{(p)}，占总期望计数 \sum_l a_{il}\theta_l^{(p)} + b_i 的比例就是 n_{ij} 的期望占 m_i 的比例。

将期望值代入完全数据对数似然，得到 Q 函数：

ML-EM算法的M步

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)}) 关于 \boldsymbol{\theta} 是严格凹函数，因此存在唯一的最大值点。在非负约束 \theta_j \geq 0 下，最大化的必要条件是Karush-Kuhn-Tucker（KKT）条件：对于所有 j，\theta_j^* \geq 0，\frac{\partial Q(\boldsymbol{\theta}^*|\boldsymbol{\theta}^{(p)})}{\partial \theta_j} \geq 0，且满足互补松弛条件。

计算偏导数并令其为零（对于内点解）：

求解得到ML-EM更新公式：

这是一个乘性更新公式：新估计等于旧估计乘以一个校正因子。由于所有项都是非负的，如果初始估计 \theta_j^{(0)} > 0，则所有后续估计都保持非负，即 \theta_j^{(p)} \geq 0 对所有 p 成立。

ML-EM算法的矩阵形式

E步和M步可以合并为一个简洁的更新公式。用矩阵-向量记号，ML-EM算法可以写成：

其中除法和乘法都是逐元素运算（Hadamard乘积/除法），\mathbf{1} 是全1向量。

算法的具体步骤如下：

第一步，从一个严格正的初始图像 \boldsymbol{\theta}^{(0)} 开始（例如均匀图像）。

第二步，正投影计算当前估计的期望数据：\mathbf{q}^{(p)} = \mathbf{A}\boldsymbol{\theta}^{(p)} + \mathbf{b}。

第三步，计算当前误差比值：\mathbf{e}^{(p)} = \frac{\mathbf{m}}{\mathbf{q}^{(p)}}（逐元素除法）。

第四步，反投影误差比值：\mathbf{c}^{(p)} = \mathbf{A}^\mathsf{T}\mathbf{e}^{(p)}。

第五步，乘性更新当前估计：\boldsymbol{\theta}^{(p+1)} = \frac{\boldsymbol{\theta}^{(p)}}{\mathbf{A}^\mathsf{T}\mathbf{1}} \mathbf{c}^{(p)}。

第六步，重复步骤二至五，直到满足收敛条件。

分母 \mathbf{A}^\mathsf{T}\mathbf{1} 是灵敏度图像，其第 j 个分量 \sum_i a_{ij} 表示体素 j 对所有响应线的总贡献能力。这个归一化因子确保了在没有噪声的理想情况下，当 \mathbf{m} = \mathbf{A}\boldsymbol{\theta}^{\text{true}} + \mathbf{b} 时，真实图像 \boldsymbol{\theta}^{\text{true}} 是算法的不动点。

ML-EM算法的收敛特性

图示展示了ML-EM算法从均匀初始图像开始的迭代过程。从一个常数图像（所有体素值相同）出发，随着迭代次数增加，图像逐渐显现出物体的结构。

在前几次迭代（1、2、3、4次）中，图像呈现出非常模糊的轮廓，只能看到大致的椭圆形状。随着迭代继续（8、16、32、64次），图像细节逐渐丰富，边缘变得更加清晰，内部结构（两个椭圆形区域）开始显现。

从频率角度观察这一现象：ML-EM算法首先恢复图像的低频成分（整体形状和亮度分布），然后逐步恢复越来越高的频率成分（精细结构和边缘）。这种从低频到高频的渐进恢复特性是ML-EM算法的一个重要特征，也是理解其收敛行为和噪声特性的关键。

ML-EM作为预条件梯度上升算法

ML-EM算法可以从另一个角度理解：它是一种对角预条件的梯度上升算法。将ML-EM更新公式改写为加法形式：

对数似然的梯度为：

用矩阵记号，ML-EM可以写成：

这个形式表明ML-EM本质上是沿梯度方向更新，但使用了一个依赖于当前估计的对角预条件矩阵 \text{diag}(\theta_j^{(p)}/\sum_i a_{ij})。预条件矩阵根据当前体素值和灵敏度自适应地调整步长：值较大的体素获得较大的更新步长，灵敏度较高的体素获得较小的步长（因为它们对数据有更强的约束）。

如果不使用这种预条件而直接进行梯度上升，收敛速度会慢得多，因为标准梯度方法在条件数较差的问题上表现不佳。

ML-EM作为函数替代方法

ML-EM还可以理解为一种函数替代算法（functional substitution algorithm），也称为优化转移（optimization transfer）或极小化-极大化（Minorization-Maximization，MM）算法。

MM算法的基本思想是：在每次迭代中，不直接最大化复杂的目标函数 L(\boldsymbol{\theta})，而是用一个更容易最大化的替代函数 S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)}) 来近似它。替代函数必须满足以下三个条件：

第一，在当前点处相切：L(\boldsymbol{\theta}^{(p)}) = S(\boldsymbol{\theta}^{(p)}, \boldsymbol{\theta}^{(p)})。

第二，在当前点处梯度相同：\left.\nabla L(\boldsymbol{\theta})\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}^{(p)}} = \left.\nabla S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)})\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}^{(p)}}。

第三，全局下界性质：L(\boldsymbol{\theta}) \geq S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)}) 对所有 \boldsymbol{\theta} 成立。

在当前点处的相切性和全局下界性保证了：最大化替代函数 S 会单调地增加原目标函数 L。迭代公式为：

图示展示了一维情况下MM算法的工作原理。黑色曲线是目标函数 L(\theta)，彩色曲线是不同迭代点处的替代函数 S(\theta, \theta^{(p)})。每个替代函数在当前点与 L 相切，且始终位于 L 的下方。通过最大化替代函数（找到彩色曲线的峰值），得到新的迭代点，该点处的 L 值必然不低于前一点。随着迭代进行，估计值 \theta^{(p)} 逐渐接近最大似然估计 \hat{\theta}_{\text{ML}} = 3.00。

设计MM算法的关键是找到易于最大化的替代函数，通常这意味着替代函数应该关于各参数可分离。

基于Jensen不等式的替代函数构造

对于泊松对数似然函数：

难以直接最大化的原因是对数函数内部有求和。为了构造替代函数，首先将期望计数 q_i = \sum_j a_{ij}\theta_j + b_i 改写为加权和的形式。假设 \theta_j^{(p)} > 0 且 q_i^{(p)} = \sum_j a_{ij}\theta_j^{(p)} + b_i > 0，定义权重：

由于权重满足 \frac{a_{ij}\theta_j^{(p)}}{q_i^{(p)}} \in [0, 1]、\frac{b_i}{q_i^{(p)}} \in [0, 1]，且 \sum_j \frac{a_{ij}\theta_j^{(p)}}{q_i^{(p)}} + \frac{b_i}{q_i^{(p)}} = 1，可以应用Jensen不等式。

对于凹函数 \log(\cdot)，Jensen不等式给出：

其中 w_j \geq 0 且 \sum_j w_j = 1。应用到对数似然的第一项：

右边定义为替代函数 S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)})。这个替代函数的关键优势在于：各 \theta_j 的贡献完全分离，可以独立最大化。最大化 S 得到的更新公式正是ML-EM算法。

替代函数的显式形式与ML-EM推导

通过Jensen不等式构造的替代函数可以写成：

可以验证这个替代函数满足MM算法的两个关键条件：在当前点处与原函数相切 L(\boldsymbol{\theta}^{(p)}) = S(\boldsymbol{\theta}^{(p)}, \boldsymbol{\theta}^{(p)})，以及在当前点处梯度相同 \left.\nabla L(\boldsymbol{\theta})\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}^{(p)}} = \left.\nabla S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)})\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}^{(p)}}。

在非负约束下最大化 S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)})，KKT条件为：

求解得到ML-EM更新公式：

值得注意的是，替代函数 S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)}) 与EM算法中的 Q 函数只相差一个不依赖于 \boldsymbol{\theta} 的常数 C，即 S(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)}) = Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)}) + C。这证明了EM算法和MM算法在泊松似然情况下得到相同的更新公式。

列表模式重建

传统的ML-EM算法基于直方图模式（histogram mode）数据，即将检测事件按响应线分组统计，形成计数向量 \mathbf{m} = \{m_i\}。列表模式（list-mode）重建则直接使用事件列表，不进行预先分组。

ML-EM可以从直方图模式转换为列表模式，这基于解关于 \mathbf{m} 的线性性。将更新公式中的求和展开：

由于 m_i 是响应线 i 上的计数，可以将其展开为对单个事件的求和：

进一步简化为对所有检测事件 e 的求和：

其中 i_e 表示事件 e 所在的响应线。

列表模式重建的优势在于不需要在重建前将数据进行直方图化，可以节省存储空间（当事件数远小于可能的响应线数时）和计算时间。这在高分辨率PET和具有连续响应线参数化的系统中特别有用。

ML-EM算法的优缺点

优点

ML-EM算法在满足一定条件时单调收敛到对数似然函数 L(\boldsymbol{\theta}) 的全局最大值。这些条件包括：对所有 i 有 m_i, b_i > 0，对所有 j 有 \sum_i a_{ij} > 0（每个体素至少被一条响应线覆盖）。收敛性基于二阶KKT条件的分析。

当使用正值初始图像时，ML-EM自动产生非负图像。这是乘性更新公式的直接结果，无需显式施加非负约束。

算法实现简单，每次迭代只需要一次正投影和一次反投影。此外，如前所述，可以直接应用于列表模式数据。

从MM算法的角度理解ML-EM，可以自然地扩展到包含先验项的MAP重建，得到广义EM（Generalized EM，GEM）算法。

缺点

ML-EM的收敛速度较慢，通常需要数十到数百次迭代才能接近收敛，这在临床应用中可能不可接受。

更重要的问题是ML-EM没有内置的正则化机制。最大似然估计本身是一个不适定问题，在噪声数据下会产生严重的噪声放大，这将在下一节详细讨论。

最大似然估计的噪声问题

图示比较了滤波反投影（FBP）和ML-EM在不同采集时间下的重建结果。左侧是T1加权MRI作为解剖参考。从左到右，PET采集时间从3600秒减少到10秒，数据中的计数统计量相应减少。

对于3600秒的长时间采集（高计数），FBP和MLEM都能产生较好的图像质量，脑部结构清晰可辨。但随着采集时间缩短，两种方法的图像质量都下降。特别是在10秒采集的极低计数情况下，两种方法都产生了严重的噪声伪影。

MLEM图像呈现出一种特征性的"夜空伪影"（night sky artefacts）：图像中出现许多孤立的高亮点，如同夜空中的星星。这种伪影的原因在于：如果ML估计 \hat{\boldsymbol{\theta}}_{\text{ML}} 是严格正的（即位于非负象限的内部），则KKT条件要求 \mathbf{1}_M - \frac{\mathbf{m}}{\mathbf{A}\hat{\boldsymbol{\theta}}_{\text{ML}} + \mathbf{b}} = \mathbf{0}，即 \mathbf{A}\hat{\boldsymbol{\theta}}_{\text{ML}} + \mathbf{b} = \mathbf{m}。但由于噪声和系统矩阵 \mathbf{A} 的病态性，这个方程通常不能精确满足。因此 \hat{\boldsymbol{\theta}}_{\text{ML}} 往往不是内点，而是位于非负象限的边界上，导致许多体素值为零或接近零，而少数体素具有很高的值。

最大似然估计不一定是最理想的解——需要通过正则化来获得临床上有用的图像。

早停作为隐式正则化

一种简单的正则化策略是早停（early stopping）：不让ML-EM算法完全收敛，而是在适当的迭代次数后终止。

图示上半部分展示了重建误差（RMSE，相对于真实图像）随迭代次数的变化。误差曲线呈现U形：初期误差快速下降（图像从模糊的初始估计向真实图像逼近），在某个迭代次数（约12次）达到最小值，之后误差反而增加（因为噪声被逐渐放大）。

图像序列显示了这一过程：迭代1次时图像非常模糊；迭代4、8次时图像逐渐清晰；迭代12次时图像质量最好；迭代48次时虽然细节更丰富但噪声也明显增加。

这种行为可以理解为：从低频常数初始图像出发，ML-EM首先恢复低频成分，然后逐渐恢复高频成分。真实的图像信息主要在中低频，而噪声则均匀分布在所有频率。因此，早期迭代恢复的主要是信号，后期迭代恢复的主要是噪声。早停相当于对高频成分进行抑制，实现了隐式的平滑正则化。

然而，早停有严重的缺陷：最佳停止点是非平稳的（依赖于图像中的具体位置）和对象依赖的（不同的成像对象有不同的最佳迭代次数）。图示下半部分的脑部图像显示，在 k = 4 次迭代时图像过于模糊，k = 16 时开始显现细节，k = 64 时噪声明显，k = 1600 时严重的夜空伪影。没有一个统一的停止准则能够在所有情况下给出最佳结果。

后平滑正则化

另一种简单的正则化策略是后平滑（post-smoothing）：对ML-EM重建结果应用低通滤波器（如高斯滤波）来抑制噪声。

图示展示了不同平滑程度的效果。上排从左到右是不同迭代次数（k = 4, 16, 64, 1600）的ML-EM结果，下排是对最终ML估计 \hat{\boldsymbol{\theta}}_{\text{ML}} 应用不同半高全宽（FWHM = 2 mm, 4 mm, 6 mm）高斯滤波器后的结果。

选择FWHM为4毫米的滤波器可以有效抑制噪声，同时保留主要的脑部结构。但后平滑的缺点是：它是在重建完成后进行的，平滑核是均匀的，不能区分图像中需要保留的边缘和需要抑制的噪声。

在实际临床应用中，常采用早停和后平滑相结合的策略：运行有限次数的迭代，然后对结果进行适度平滑。值得注意的是，与分辨率建模的情况类似，平滑滤波器也可以整合到重建过程中（称为sieves方法），而不是作为后处理步骤。

有序子集期望最大化（OSEM）

ML-EM算法的一个主要缺点是收敛速度慢。有序子集期望最大化（Ordered Subsets Expectation Maximization，OSEM）是一种加速技术，显著提高了收敛速度。

OSEM属于块迭代方法（block-iterative methods）：每个更新步骤只使用数据的一个子集 S，而不是全部数据 M。这比同时使用所有数据的方法（如标准ML-EM）快得多。

有序子集方法利用了泊松数据的独立性，使得对数似然可以分解为各子集贡献的和：

其中 N_s 是子集数目，L_q(\boldsymbol{\theta}) 是仅基于子集 S_q 中数据的对数似然。子集构成数据的一个划分：\cup_{q=1}^{N_s} S_q = M，且要求每个体素至少被每个子集中的某些响应线覆盖，即 \sum_{i \in S_q} a_{ij} > 0 对所有 j 和 q 成立。

OSEM算法依次对每个子集应用EM更新：

遍历所有子集完成一次完整迭代：\boldsymbol{\theta}^{(p,0)} = \boldsymbol{\theta}^{(p-1,N_s)}。

子集平衡

为了使OSEM有效工作，子集应该具有平衡性：每个子集应该具有相等的检测概率，能够均匀地"看到"整个图像空间。

子集划分的常用策略包括：随机选择响应线或事件分配到各子集；或者按方位角划分，使每个子集内的响应线在角度上有最大的分散。图示展示了按方位角将响应线分成四个子集的例子：原始的全部响应线被分成四组，每组覆盖不同的角度方向，组内的响应线尽可能均匀分布。

更精确地说，关键条件是子集梯度平衡：

当这个条件成立时，将EM视为预条件梯度上升，每个子集的更新大致给出总梯度方向的 1/N_s，从而实现近似线性的加速。

OSEM的加速效果

图示比较了ML-EM和OSEM（使用16个子集）的收敛速度。ML-EM在16、32、48、64次迭代后的图像分别与OSEM在1、2、3、4次迭代后的图像质量相当。

这说明在早期迭代阶段，N 次ML-EM迭代大致等价于 N/16 次OSEM迭代（当使用16个子集时）。换言之，OSEM在早期提供了约16倍（即子集数目倍）的加速。

这种加速的直观解释是：每次子集更新都提供了关于图像的部分信息，虽然每次更新的信息量较小，但更新频率更高。在早期迭代阶段，图像与真实解差距较大，任何方向的更新都是有益的，因此频繁的子集更新能够快速接近解。

OSEM的优缺点

优点

OSEM在早期迭代阶段提供了显著的加速，这对于临床应用至关重要，使得迭代重建在实际扫描流程中变得可行。

OSEM实现简单，只需对ML-EM进行最小的修改：主要是将数据和系统矩阵按子集组织，循环遍历子集进行更新。

与ML-EM一样，OSEM可以应用于列表模式数据，也可以扩展到MAP重建（称为OSMAP或OS-GEM）。

缺点

OSEM的主要理论缺陷是不保证收敛到 \hat{\boldsymbol{\theta}}_{\text{ML}}。在实践中经常观察到极限环（limit cycles）现象：算法在一组解之间振荡，而不是收敛到单一的点。

有序子集方法的收敛性分析

OSEM通常不收敛的原因可以从KKT条件的角度理解。假设算法收敛到一组固定点 \{\boldsymbol{\theta}^{*,q}\}（每个子集对应一个），则对于每个子集 q：

且要求 \theta_j^{(*,N_s)} = \theta_j^{(*,0)}（完成一个周期后回到起点）。

然而，没有保证对于任意 \boldsymbol{\theta}^{(*,q)}，总梯度为零：

如果子集梯度不能相互抵消，算法将在不同子集的"拉扯"下振荡，形成极限环。

为了获得有保证收敛的有序子集方法，需要引入额外的机制来逐渐减小步长。这类方法包括：行动最大似然算法（Row-Action Maximum Likelihood Algorithm，RAMLA）通过逐渐减小的松弛因子来保证收敛；块序贯正则化EM（Block Sequential Regularized EM，BSREM）通过正则化项来稳定收敛。这些方法在额外条件下可以证明全局收敛。

尽管存在收敛性问题，在实际发射断层成像中，早停的OSEM仍然是最广泛使用的算法，因为它在实践中表现良好且实现简单。

全局收敛的松弛有序子集方法：BSREM-II

BSREM-II（Block Sequential Regularized EM-II）是一种保证全局收敛的有序子集算法。假设目标函数 f = L + R 是严格凹的且可微，BSREM-II可以视为投影预条件梯度下降算法。

算法设计需要解决三个关键问题：

第一，确保迭代序列保持正值且有界，即位于有界凸集 \mathcal{B} = \{\boldsymbol{\theta} : 0 \leq \theta_j \leq U, \forall j\} 内。这通过修改预条件矩阵实现：

当体素值接近上界 U 时，预条件因子减小，防止超出边界。

第二，确保迭代序列保持在 \mathcal{B} 的内部（避免卡在边界上）。这通过将迭代结果投影到一个内部闭凸集 \mathcal{T} \in \text{int } \mathcal{B} 来实现。

第三，确保收敛性。这通过选择递减的预条件梯度步长 \alpha_n 来实现，要求 \sum_n \alpha_n = +\infty（保证能够到达任意远的点）且 \sum_n \alpha_n^2 < +\infty（保证步长最终趋于零）。

BSREM-II的更新公式为：

其中 \mathcal{P}_{\mathcal{T}} 是到集合 \mathcal{T} 的投影算子。

其他并行/分布式/随机方法

基于对数似然的可分性 L(\boldsymbol{\theta}) = \sum_{q=1}^{N_s} F_q(\mathbf{A}\boldsymbol{\theta})，现代凸优化技术可以用于处理块/并行最小化问题。

分布式或块策略包括：ADMM（交替方向乘子法）算法将问题分解为多个子问题并行求解；原始-对偶混合梯度方法同时更新原始变量和对偶变量；对偶块坐标前向-后向算法在对偶空间进行块更新。

随机方法是近年来在PET重建中提出的新方向，其核心思想是避免在每次迭代中评估所有子集。具体做法是在每次迭代开始时随机抽取一个子集，仅基于该子集进行更新。为了保证收敛，需要使用方差缩减技术来控制随机梯度的方差。

最大似然估计的噪声与偏差

图示比较了OP-OSEM（有序子集后平滑OSEM）和FBP在不同采集时间下的偏差特性。横轴是采集时间（对数刻度，从1秒到10000秒），纵轴是估计的偏差百分比。测试对象是一个均匀的椭圆体模。

FBP（绿色方块）在所有采集时间下都保持接近零的偏差，这是线性重建方法的固有特性。而OP-OSEM（蓝色圆点）在长采集时间（高计数）下偏差接近零，但随着采集时间缩短（低计数），偏差显著增加，在极短采集时间下可达近10%。

这说明ML-EM在极低计数情况下存在偏差，而FBP是无偏的。

非线性与偏差的进一步分析

图示更详细地展示了FBP和EM-ML的行为差异。对于不同的采集时间（3600秒到10秒），分别显示了单次实现（1 replica）和多次实现的平均图像（Sum image）。

对于FBP：单次实现随着采集时间缩短噪声增加，但平均图像在所有采集时间下都保持一致的结构，没有明显偏差。这反映了FBP的线性性质：噪声的平均值为零，多次平均后消除。

对于EM-ML：单次实现在低计数下出现严重的夜空伪影，且平均图像随着采集时间缩短而变化——不仅噪声增加，图像结构本身也发生改变。这说明ML-EM是非线性的：\mathbb{E}[\hat{\boldsymbol{\theta}}(\mathbf{m})] \neq \hat{\boldsymbol{\theta}}(\mathbb{E}[\mathbf{m}])。

非线性导致偏差的机制是：在低计数情况下，非负约束变得活跃（许多体素被迫取零值），而非负约束的效果是非线性的，取决于噪声的具体实现。

偏差补偿方法

ML-EM在渐近意义下（无噪声极限）是无偏的。研究假设极低计数数据下的偏差与泊松似然的非负性约束有关：当期望计数很小时，非负约束强制许多体素取零值，导致系统性的低估。

一种补偿方法是NEG-ML（Negative-allowing ML）：修改泊松分布在零附近的行为，允许期望值为负。具体做法是在某个阈值 \psi 以下将标准泊松对数似然（蓝色曲线）替换为一个延伸到负值区域的修改函数（红色曲线）。

图示显示了这种修改：原始泊松对数似然在 \bar{n} = 0 处趋向负无穷，强制期望值非负；修改后的函数在小正值处平滑过渡，允许期望值取小的负值。这种修改放松了非负约束，可以减少低计数情况下的偏差。

引入正则化：MAP算法

为了克服ML估计的噪声问题，需要引入正则化。这可以通过两种等价的框架来实现。

最大后验（贝叶斯框架）

在贝叶斯框架中，引入图像的先验分布 p(\boldsymbol{\theta})，MAP估计最大化后验概率：

当先验采用Gibbs分布 p(\boldsymbol{\theta}) = \frac{1}{Z}\mathrm{e}^{-\beta U(\boldsymbol{\theta})} 时：

惩罚最大似然（优化框架）

从优化角度，直接在目标函数中添加惩罚项 U(\boldsymbol{\theta})：

两种框架在数学形式上完全等价，区别在于解释：贝叶斯框架将惩罚项解释为先验信息的负对数，优化框架将其视为对不期望图像特性的惩罚。超参数 \beta 控制数据拟合与正则化之间的权衡。

MAP优化：广义EM算法

广义EM（Generalized EM，GEM）算法将标准EM算法扩展到包含正则化项的MAP估计。

标准ML的EM算法

回顾标准ML-EM，E步计算完全数据对数似然的条件期望：

M步最大化 Q 函数：\boldsymbol{\theta}^{(p+1)} = \arg\max_{\boldsymbol{\theta} \in \Theta} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(p)})。

展开后：

MAP的广义EM算法

对于MAP估计，E步计算包含惩罚项的期望：

由于惩罚项 U(\boldsymbol{\theta}) 不依赖于完全数据 \mathbf{n}，期望直接穿过：

M步变为：

与标准EM相比，唯一的区别是M步需要同时考虑数据拟合项和惩罚项。M步的求解难度取决于惩罚函数 U(\boldsymbol{\theta}) 的形式。

GEM：独立体素先验

当惩罚函数假设体素之间相互独立时，M步可以直接求解。以独立高斯先验为例：

M步的目标函数关于每个 \theta_j 可以分离，对 \theta_j 求导并令其为零，得到一元二次方程。求解得到闭式更新公式：

其中 \alpha_j = \sum_{i=1}^{J} a_{ij} 是灵敏度，\theta_j^{\text{EM}(p)} = \theta_j^{(p)} \sum_{i=1}^{I} a_{ij} \frac{m_i}{\sum_{l=1}^{J} a_{il}\theta_l^{(p)} + b_i} 是标准EM更新的中间结果。

这个闭式解的结构可以理解为：在EM更新结果和先验期望值之间进行加权平均，权重由正则化参数 \beta 和方差 \sigma_j^2 决定。

值得注意的是，这种求解方法也可以用于计算泊松负对数似然的近端算子（proximal operator），在近端算法中有重要应用。

GEM：空间相关先验

当惩罚函数包含体素之间的空间相关性时（如MRF先验），M步的目标函数不再关于各 \theta_j 可分离：

由于 U(\boldsymbol{\theta}) 包含相邻体素之间的耦合项，不能像独立先验那样直接求解闭式解。

单步延迟方法（One-Step-Late，OSL）

一种简单的近似方法是单步延迟（OSL）：在计算惩罚项的梯度时，使用上一次迭代的估计值 \boldsymbol{\theta}^{(p)} 而不是当前值：

OSL的优点是实现简单，几乎与标准ML-EM一样容易。但其主要缺点是不保证收敛到 \hat{\boldsymbol{\theta}}_{\text{MAP}}，且当 \beta 较大时可能出现不稳定。

MAP优化：预条件梯度上升

另一类方法是基于梯度的优化算法。回顾ML-EM可以写成预条件梯度上升的形式：

惩罚预条件共轭梯度（PPGML）

惩罚预条件共轭梯度最大似然（Penalized Preconditioned Gradient Maximum Likelihood，PPGML）算法在预条件中同时考虑似然项和惩罚项的曲率：

其中 \mathcal{P}_{\geq 0} 是到非负象限的投影（将负值截断为零）。预条件矩阵的分母中增加了惩罚项Hessian的对角元素 \theta_j^{(p)} \beta \frac{\partial^2 U}{\partial \theta_j^2}，以适应正则化带来的额外曲率。

PPGML实现简单，但与OSL类似，不保证收敛到 \hat{\boldsymbol{\theta}}_{\text{MAP}}。

MAP优化：二次惩罚的替代函数

对于二次惩罚函数，可以构造使M步可直接求解的替代函数。二次惩罚的形式为：

这个惩罚函数由于涉及相邻体素的耦合，直接优化需要求解联立方程组。De Pierro提出了一种替代函数：

可以验证 V(\boldsymbol{\theta}, \boldsymbol{\theta}^{(p)}) \geq U_{\text{quad}}(\boldsymbol{\theta}) 且在 \boldsymbol{\theta} = \boldsymbol{\theta}^{(p)} 处相等（相切），因此 -V 是 -U_{\text{quad}} 的下界。

使用这个替代函数，M步变为：

由于 V 关于每个 \theta_j 只包含 \theta_j 自身和固定的邻域值 \theta_k^{(p)}，目标函数变得关于各 \theta_j 可分离，可以逐体素求解。

二次惩罚的闭式更新

对于二次替代函数，M步可以得到闭式解：

其中：W_j = \sum_{k \in N_j} \omega_{kj} 是体素 j 的邻域权重总和，\gamma_j^{(p)} = \theta_j^{(p)} W_j + \sum_{k \in N_j} \theta_k^{(p)} \omega_{kj} 是当前估计的邻域加权和，\alpha_j = \sum_{i=1}^{J} a_{ij} 是灵敏度，\theta_j^{\text{EM}(p)} = \theta_j^{(p)} \sum_{i=1}^{I} a_{ij} \frac{m_i}{\sum_{l=1}^{J} a_{il}\theta_l^{(p)} + b_i} 是EM中间结果。

可以验证当 \beta \to 0 时：

即退化为标准ML-EM更新，符合预期。

这种基于替代函数的算法保证收敛，但收敛速度较慢。对于非二次但满足一定条件的势函数（凸、对称、可微、曲率非增），也存在类似的闭式更新公式。

MAP优化：其他算法

当惩罚函数 U(\boldsymbol{\theta}) 是凸的、可微的且具有 \beta-Lipschitz连续梯度时，可以使用广泛的MM类算法来计算MAP估计。对于凸正则化问题，基于分裂的原始-对偶近端算法特别有效。这类算法的核心是能够高效计算近端算子：

近端算子可以理解为在当前点附近寻找一个既接近当前点又使函数值较小的点。

ADMM算法

交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）是一种强大的分裂算法，特别适合处理可分离结构的优化问题。对于简单惩罚函数 U 的凸问题，ADMM的迭代步骤如下：

第一步（重建步骤）：

第二步（正则化步骤）：

第三步（对偶变量更新）：

ADMM将原问题分解为两个子问题：第一步主要处理数据拟合（似然最大化），第二步主要处理正则化。两个子问题通过对偶变量 \mathbf{u} 进行协调。

惩罚函数的替代方案：核EM

除了显式设计惩罚函数 U(\boldsymbol{\theta}) 外，还有一种不同的正则化思路：将图像表示为某个线性算子的输出。

假设图像可以表示为：\boldsymbol{\theta} = \mathbf{K}\mathbf{x}，其中 \mathbf{x} 是潜在表示（latent representation），\mathbf{K} 是一个字典矩阵，其列可以视为原子（atoms）或基向量。

核EM算法

核EM（Kernel EM）使用ML-EM来估计潜在表示 \mathbf{x}：

最终的图像估计为 \hat{\boldsymbol{\theta}}_{\text{KEM}} = \mathbf{K}\hat{\mathbf{x}}_{\text{EM}}。

核矩阵 \mathbf{K} 通常从辅助图像（如MRI或CT）导出，利用核方法构造。典型的构造方式是：

其中 \mathbf{f}_j 是辅助图像中体素 j 的特征向量（如强度值或局部纹理特征），\mathcal{N}_j 是体素 j 的邻域，\alpha 控制核的带宽，Z(j) = \sum_{l \in \mathcal{N}_j} \exp(-\alpha\|\mathbf{f}_j - \mathbf{f}_l\|^2) 是归一化因子。

这种构造的效果是：在辅助图像中特征相似的相邻体素，在核矩阵中具有较大的权重，从而在重建时被联合处理。核EM通过限制解空间（\boldsymbol{\theta} 必须在 \mathbf{K} 的列空间中）实现了隐式正则化。

解剖信息与发射断层成像

发射断层图像的空间形态与底层解剖结构密切相关。由于PET和SPECT的空间分辨率有限，重建图像缺乏高频细节。这些缺失的高频信息可以从其他成像模态（如MRI或CT）的解剖图像中借用。

图示展示了这种关联：左侧的MRI和PET脑部图像显示了明显的结构对应关系；右侧的猴脑自放射影像图更清晰地展示了示踪剂分布与解剖结构之间的联系——神经受体示踪剂的分布精确地跟随脑部的解剖分区。

利用这种结构相关性，可以将解剖图像的信息融入PET/SPECT重建，以提高图像质量。

MRF与解剖信息

一种利用解剖信息的简单方法是修改MRF先验中的邻域权重。基本思想是：如果两个相邻体素属于不同的解剖结构（根据MRI或CT图像判断），则停用它们之间的平滑惩罚（令 \omega_{kj} = 0）。

解剖引导的MAP估计为：

其中邻域权重 \omega_{kj} 由解剖图像确定：当体素 k 和 j 属于同一解剖区域时 \omega_{kj} > 0，否则 \omega_{kj} = 0。

图示展示了这种方法的效果：左侧是MRI图像，中间是PET图像和MRI分割结果的叠加，右侧是使用解剖先验重建的PET图像。通过停用跨越解剖边界的平滑惩罚，重建图像能够更好地保持与解剖结构一致的边缘。

解剖引导重建的偏差来源

使用解剖信息引导PET重建需要谨慎，因为可能引入偏差。偏差的主要来源包括两个方面。

首先是模型有效性问题：解剖先验假设PET信号的边界与解剖边界一致。但功能信息（如代谢活动、受体密度）并不总是严格遵循解剖边界。例如，肿瘤可能跨越正常的解剖分区，或者神经退行性疾病可能在解剖结构完整的区域导致功能异常。在这些情况下，强制PET图像遵循解剖边界会引入偏差。

其次是技术实现问题：解剖引导需要将MRI/CT图像与PET图像进行配准，还需要对解剖图像进行分割以提取区域边界。配准误差和分割误差都会传播到重建结果中。

图示展示了处理流程：MRI图像经过分割得到解剖区域划分，PET图像经过配准与MRI对齐，然后将分割信息应用到PET重建中。每个步骤都可能引入误差。

解剖MAP：Bowsher先验

解剖引导重建的一个实际困难是如何提取解剖边界。MRI图像的分割是一个复杂的任务，尤其对于病理情况，自动分割可能不可靠。

Bowsher先验提供了一种无需显式分割的解剖引导方法。其核心思想是：不进行MRI分割，而是在体素 j 的邻域 N_j 中，选择MRI强度值最接近体素 j 的 M 个体素作为"有效邻域" M_j。

其中 \delta_{jk} 是体素 j 和 k 之间的距离权重。

图示说明了这个过程：体素 j 的完整邻域 N_j 包含27个体素（3×3×3立方体），但Bowsher先验只选择其中MRI强度最相似的 M = 6 个体素组成有效邻域 M_j。只有这6个体素参与平滑惩罚的计算。

Bowsher先验的优点是自动适应局部解剖结构，无需预先分割。它隐式地假设：MRI强度相似的相邻体素很可能属于同一组织类型，因此应该在PET重建中被一起平滑。

先验对重建图像的影响

图示展示了不同先验对重建结果的影响。从左到右分别是：均匀先验（相当于ML估计）、二次先验、边缘保持先验、以及使用MRI信息的解剖先验。

上排显示MAP估计的均值（通过多次噪声实现取平均得到，反映偏差）：均匀先验的均值图像噪声最大但偏差最小；二次先验大幅降低噪声但也明显模糊了边缘；边缘保持先验在降噪的同时较好地保持了边缘；MRI引导先验产生了最清晰的边缘，但如果解剖与功能不完全匹配，可能引入偏差。

下排显示MAP估计的方差（反映估计的不确定性）：均匀先验方差最大；二次先验显著降低方差但空间分布较均匀；边缘保持先验在边缘处方差较大（因为边缘处的正则化较弱）；MRI引导先验在组织内部方差很低，但在解剖边界处方差较高。

这说明先验的选择同时影响对比度恢复（均值）和估计的不确定性（方差），需要根据具体应用场景权衡。

MAP重建的现存问题

尽管MAP方法提供了一个原则性的正则化框架，但在实际应用中仍面临诸多挑战。

首先，先验设计依赖手工制作。目前使用的各种势函数（二次、TV、Huber等）主要是数学上方便的选择，而非基于对真实医学图像统计特性的深入研究。如何从数据中学习更合适的先验是一个活跃的研究方向。

其次，优化算法的速度和收敛性各异。不同算法（OSL、PPGML、替代函数方法、ADMM等）在收敛速度和收敛保证方面差异很大，且并非所有常用算法都有严格的收敛证明。

第三，解剖引导方法对配准、分割以及功能/解剖不匹配的敏感性问题。当PET信号分布与解剖结构不完全一致时，强制解剖约束可能引入临床相关的偏差。

第四，空间分辨率的非均匀性问题。正则化通常导致空间变化的分辨率——高对比度区域的分辨率高于低对比度区域，这可能使定量分析复杂化。

第五，除了解剖信息外，MRI还提供其他有价值的信息（如功能MRI、扩散张量成像），如何有效利用这些信息是开放问题。

最后，协同PET/MRI重建是一个新兴方向：不是将MRI作为固定的先验信息，而是同时重建PET和MRI图像，共享正则化项，以实现真正的多模态融合。

机器学习图像重建

这部分内容属于断层成像重建的第三个分支。在前面的学习中，我们已经接触了解析重建方法和迭代重建方法，现在进入机器学习重建这一现代方法。这一领域的发展可以追溯到近年来深度学习在医学成像中的广泛应用，两篇代表性的综述文章分别发表于2018年的IEEE Transactions on Medical Imaging和2021年的IEEE Transactions on Radiation and Plasma Medical Sciences，标志着这一研究方向的兴起与成熟。

传统方法与深度学习方法的对比

经典的基于模型的方法

传统的图像重建方法建立在物理采集模型的基础上，即我们明确知道探测器如何采集数据、射线如何穿过人体组织、信号如何衰减等物理过程，并将这些过程用数学方程精确描述。在此基础上，通过正则化手段对解施加约束，比如要求图像具有平滑性、允许存在边缘不连续性、保证像素值非负等简单的数学性质。迭代算法在这类框架下运行，并且具备可证明的收敛性质。

这类方法的优势在于可靠性：我们对估计结果有明确的控制能力，包括算法是否收敛、解的稳定性如何、不确定性有多大，这些都可以从理论上分析。然而其劣势在于通用性带来的代价——由于模型是通用设计的，在面对特定的成像任务、特定的扫描协议或特定的患者群体时，往往无法达到最优性能。为了适应具体场景，需要进行大量的超参数调节工作。

基于深度学习的数据驱动方法

深度学习方法采用完全不同的思路：在重建过程中引入深度神经网络作为核心组件或辅助模块。这类方法的核心优势在于能够针对特定任务、特定协议甚至特定患者群体进行优化，从而获得更好的图像对比度、更高的信噪比、更准确的定量结果。此外，深度学习方法还能支持更快的数据采集、更低的辐射剂量、更精细的动态数据采样，并且在重建速度上也可能更快。

然而深度学习方法面临两个核心挑战。第一是稳定性与可靠性问题：发射断层成像本质上是一个服从泊松统计的病态逆问题，深度学习方法如何在这种条件下保持稳定仍是开放性问题。第二是泛化能力问题：训练深度网络需要大量标注数据，但医学成像领域的数据集通常规模有限，每个患者的检查次数也很少，这对监督学习构成了根本性限制。

深度学习的四种应用策略

前处理与后处理

第一种策略是将深度神经网络用于数据流的前端或后端。在前处理阶段，网络作用于原始测量数据，目的是降低噪声或提升空间分辨率；在后处理阶段，网络作用于已经重建好的图像，同样用于去噪或增强。这种方式的特点是网络与传统重建算法相对独立，只是在输入或输出端做改进。

直接重建

第二种策略是用深度网络完全替代传统重建算法，直接学习从原始测量数据到重建图像的映射关系。设 m 为原始测量数据，F 为深度网络所表示的映射函数，\vartheta 为网络参数，则重建图像的估计为：

这种端到端的方式完全依赖网络从数据中学习重建规则，不显式利用物理模型。

带一致性约束的展开网络

第三种策略是将传统的迭代重建算法展开成一系列网络模块。迭代算法的每一次迭代对应网络中的一个块，算法中的某些运算（如正则化项的近端算子）被替换为可学习的网络层。这种方式保留了迭代算法的结构，同时通过学习来提升性能，并且可以施加数据一致性约束，确保重建结果与测量数据在物理上相容。

用于正则化的深度学习

第四种策略是在传统的最大后验概率重建框架内，用深度神经网络来构建先验模型。具体而言，可以用网络来学习正则化项的超参数，或者在某个正则化函数族中学习一个自适应的先验分布。这种方式保持了MAP框架的整体结构，仅将其中的先验建模部分替换为数据驱动的网络。

深度学习图像重建的具体方法

后处理方法

后处理方法是将深度神经网络应用于已经通过传统算法重建好的图像上，目的是进一步降低噪声或提升图像质量。这类方法的网络架构和学习策略多种多样，可以根据是否需要标签数据分为监督方法和无监督方法两大类。

监督学习方法

监督方法需要配对的训练数据，典型应用场景是低剂量成像：用低剂量采集的图像作为输入，用高剂量或长时间采集的高质量图像作为标签进行训练。常用的网络架构包括残差网络 ResNet、编码器-解码器结构的 U-Net、生成对抗网络 GAN 及其各种变体，以及近年来引入自注意力机制的网络结构如 Transformer。残差网络通过跳跃连接让网络学习输入与目标之间的残差，更容易训练深层网络；U-Net 的编码器逐步提取多尺度特征，解码器逐步恢复空间分辨率，跳跃连接保留细节信息；GAN 通过生成器和判别器的对抗训练，能够生成视觉上更逼真的图像；Transformer 的自注意力机制可以捕获图像中的长程依赖关系。

无监督学习方法

无监督方法不需要配对的干净图像作为标签。Deep Image Prior（DIP）利用网络结构本身作为正则化，将随机初始化的网络拟合到单张噪声图像上，由于网络结构对自然图像具有归纳偏置，早期迭代会先拟合图像的主要结构，后期才拟合噪声，因此在适当时机停止训练即可得到去噪结果。Noise2Noise 证明了即使训练数据的输入和目标都是噪声图像，只要噪声是零均值且相互独立的，网络仍然可以学会去噪。Noise2Void 更进一步，只需要噪声图像本身，通过盲点网络的设计让网络无法直接看到要预测的像素，从而避免学习恒等映射。

后处理方法面临的核心问题

这类方法需要回答几个关键问题：重建图像经过网络处理后，分辨率和噪声特性会如何变化？如果使用迭代算法重建，网络输入的图像是否已经收敛？评价指标应该采用定量指标（如均方误差）还是感知指标（如结构相似性或视觉质量）？网络处理后的图像是否仍然与原始测量数据保持一致？这些都是实际应用中需要仔细考量的问题。

直接重建方法

直接重建是指用深度网络直接学习从原始测量数据到重建图像的端到端映射，完全绕过传统重建算法。这种方法的核心挑战在于：投影空间和图像空间的维度、结构完全不同，网络需要学习一个非常复杂的非线性变换。

AUTOMAP

AUTOMAP 是直接重建方法的早期代表性工作，最初为 MRI 重建设计，后来也被应用于 2D PET 重建。其网络架构分为两个阶段：首先使用全连接层学习从传感器数据到图像的映射，这一步相当于学习逆 X 射线变换或逆 Radon 变换；然后使用卷积层对初步结果进行去噪和细化。全连接层的作用是建立投影空间到图像空间的全局映射关系，因为投影数据中每个元素都可能与图像中所有像素相关；卷积层则利用图像的局部相关性进行后处理。这种设计的代价是参数量巨大，达到约 \mathcal{O}(800M) 即 8 亿参数量级，这是因为全连接层的参数数量与输入输出维度的乘积成正比。

DeepPET

DeepPET 采用全卷积的编码器-解码器架构，将 PET 正弦图作为输入，逐步变换到一个学习得到的潜在空间表示，然后通过上采样恢复到 PET 图像域。与 AUTOMAP 相比，DeepPET 完全使用卷积操作，参数量降低到约 \mathcal{O}(60M) 即 6000 万量级，但仍然只能处理 2D 重建。从网络结构图可以看到，编码器部分使用多层卷积逐步降低空间维度、增加特征通道数，解码器部分则逐步恢复空间分辨率。右侧的对比图显示，DeepPET 的重建结果在视觉上与传统 OSEM 重建相当。

DirectPET

DirectPET 的设计目标是减少参数量以支持多层三维重建。其核心思想是：全连接层之所以参数量大，是因为它试图学习一个完全自由的映射，而实际上投影数据到图像的变换具有明确的物理结构——Radon 变换的逆。DirectPET 提出了物理启发的 Radon 反演层，首先对预校正的 TOF 数据进行 2D TOF-FORE 重分箱，然后通过一个带有可学习掩码的 Radon 反演层完成从投影空间到图像空间的变换。这个反演层内嵌了 Radon 变换的数学结构，只需要学习一些调制参数，而不是完全自由的映射。这种设计使参数量减少了约 15 倍。

FastPET

FastPET 是一种混合策略，利用飞行时间信息来简化重建问题。当 TOF 分辨率足够高时，可以直接根据时间差信息在图像空间中近似定位湮灭事件的位置，得到所谓的 histo-image，这张图像沿 TOF 方向存在模糊。重建问题因此从完整的逆 Radon 变换简化为图像空间中的去模糊问题，只需要用 U-Net 进行后处理即可。由于省去了复杂的投影-反投影计算，FastPET 的速度比 OSEM 快 67 倍，比 FBP 快 50 倍。从流程图可以看到：PET/CT 扫描仪输出 list-mode 符合事件流，经过直方图化和最可能湮灭位置计算得到 histo-image，结合衰减图输入神经网络，输出最终图像体积。

直接重建方法的评价与局限

从上述几种方法的发展可以看出，当网络设计中融入更多物理知识时（如 X 射线变换结构、TOF 信息），所需的参数量显著减少。这些方法确实实现了可行性验证，并且重建速度很快。然而，直接重建方法存在几个根本性的局限。

首先，这些方法通常需要对原始数据进行预校正，包括随机符合校正、散射校正、衰减校正等，这意味着仍然依赖传统的数据处理流程。其次，目前没有充分的证据表明直接重建方法在图像质量上优于基于模型的迭代重建方法。第三，这些方法没有显式地强制数据一致性约束，即网络输出的图像前向投影后不一定与原始测量数据吻合。第四，由于缺乏物理约束，网络可能产生不稳定的输出甚至幻觉，即生成实际不存在的结构。综合这些因素，直接重建目前并不是首选的实用策略，但它为理解深度学习在图像重建中的能力和局限提供了有价值的探索。

直接重建的幻觉问题与混合策略

直接重建方法的不稳定性

直接重建方法的一个严重缺陷是容易产生幻觉，即生成输入数据中并不存在的结构。图中展示了一个典型实验：第一行是原始MRI图像 A(x)，后续几列分别是在原始数据上添加了不同扰动 r_1, r_2, r_3, r_4 后的结果 A(x+r_i)，这些扰动在视觉上几乎不可察觉。第二行是AUTOMAP对这些扰动输入的重建结果 f(A(x+r_i))，可以看到即使输入的扰动很小，输出图像却出现了明显的伪影和失真。第三行是传统最先进方法的重建结果，对同样的扰动保持了稳定。这个实验揭示了端到端深度学习方法的脆弱性：网络可能对输入空间的某些方向极其敏感，微小的扰动会被放大成严重的伪影。

PET成像的特殊挑战

PET成像面临着更加严峻的条件。首先，PET重建是一个大规模的病态逆问题，测量数据服从泊松统计，低计数情况下噪声非常严重。图中展示了模拟的FDG脑PET数据切片：D=460M 表示全剂量采集（4.6亿计数），D/5 和 D/20 分别表示将计数降低到五分之一和二十分之一，可以看到低计数数据的噪声急剧增加。其次，医学成像领域的训练数据通常很有限，每种检查类型可能只有 \mathcal{O}(20\text{-}100) 例样本，这样的数据量是否足够覆盖医学场景中的多样性是一个问题。在低计数条件下，幻觉现象预期会更加严重，因为网络需要从噪声更大的输入中推断图像内容，更容易填充实际不存在的结构。

混合重建策略

鉴于直接深度学习重建的不可靠性，研究者提出了混合策略，将传统基于模型的方法与深度学习相结合。这类方法统称为 unrolling（展开）、synthesis（合成）或 Plug-and-Play（即插即用）方法。

混合方法的核心思想是将重建算法分解为两类模块：一类是来自凸优化理论的重建模块，负责执行数据一致性步骤，这部分严格基于物理模型，确保重建结果与测量数据相符；另一类是深度学习架构，充当学习得到的自适应非凸正则化，负责利用图像的先验知识进行去噪或特征增强。这种设计的优势在于：我们对重建过程有更强的数学刻画能力，学习到的参数具有可解释性，同时仍然保留了深度学习针对特定任务优化的能力。一个待研究的问题是这类方法的数值效率如何。

混合方法与传统方法的对比

将经典的基于模型方法与混合机器学习方法并列对比：经典方法基于物理采集模型，使用正则化强制平滑性、不连续性、正定性等数学性质，迭代算法具有收敛性保证，优点是可靠（对收敛性、稳定性、不确定性有控制），缺点是通用性导致在特定场景下性能欠佳，需要调参。

混合机器学习方法则有三种主要形式：第一种是在展开算法中加入数据一致性层；第二种是在合成方法中使用隐式约束；第三种是在 Plug-and-Play 方法中利用不动点收敛性。这些方法的优势包括：针对特定任务可以获得更好的对比度和信噪比；学习效率更高，不需要海量数据；可靠性更强，能保证数据一致性，Plug-and-Play 方法有收敛性结果。仍然存在的问题是：对训练数据库之外的变化的鲁棒性如何？是否仍会产生幻觉？

直接深度学习重建方法（用单一大型网络从高质量数据库学习数据到图像的直接映射）虽然反演速度快、针对特定任务效果好，但存在严重缺陷：需要大量学习数据、不可靠且不稳定、缺乏数据一致性保证、对训练分布外的输入鲁棒性差、容易产生幻觉。因此这种策略被否定，混合策略成为首选方向。

展开迭代重建与深度学习

展开方法的目标

展开方法的目标是保留发射断层成像的物理和统计模型 L(\theta|m)，即似然函数描述的测量数据 m 与未知图像 \theta 之间的概率关系，同时用深度学习来建模我们不确定的部分——整个势函数 U(\theta)，也就是先验项。更准确地说，是用网络来学习算子 \mathcal{D}^{(p)}，这个算子的作用相当于：

其中 \text{Id} 是恒等算子，\nabla U(\theta) 是势函数的梯度，\lambda 是步长。在梯度下降框架中，\theta - \lambda \nabla U(\theta) 正是朝着最小化 U(\theta) 方向迈进一步的结果。用网络 \mathcal{D}^{(p)} 替代这个更新步骤，相当于让网络学习一个隐式的正则化效果。

展开方法的实现

具体实现方式是将传统的迭代重建算法展开成 P 个块，每个块 p 中包含一个神经网络 \mathcal{D}^{(p)}(\theta^{(p)}; \vartheta^{(p)})，其中 \theta^{(p)} 是第 p 块的输入图像，\vartheta^{(p)} 是该块网络的可学习参数。从流程图可以看到：初始图像进入 Block 1，经过一系列块的处理，最终到达 Block N 输出重建结果。每个块内部的结构包括：去噪网络 F^n(x^n) 对当前图像进行处理，同时有基于测量数据 m 和先验的 MAP-EM 重建步骤（执行 k=1,...,K 次 EM 更新），然后将去噪结果和 EM 更新结果进行组合，得到下一个块的输入 x^{n+1}。这种设计确保了每个块都既利用了网络的去噪能力，又通过 EM 更新步骤保持了与测量数据的一致性。

展开算法实例：FBSEM-Net

算法目标与框架

FBSEM-Net 的目标是用深度网络学习势函数 \beta U(\theta)，使得最终的 MAP 估计：

能够匹配高质量图像。这里 L(\theta|m) 是对数似然函数，描述测量数据 m 在给定图像 \theta 下的概率；U(\theta) 是势函数，起正则化作用；\beta 是控制正则化强度的超参数。整个优化问题的含义是：在所有可能的图像中，找到一个既能很好解释测量数据（似然项大），又符合先验假设（势函数小）的图像。

实现这一目标的方法是将前向-后向分裂（Forward-Backward Splitting, FBS）算法展开成网络结构，每个块中用卷积神经网络来建模 \beta \nabla U(\theta)，即势函数的梯度乘以正则化系数。

前向-后向分裂算法

FBS 算法将优化问题分解为两个交替执行的步骤。第一步是正则化更新（前向步骤）：

这一步从上一次迭代的结果 \theta^{(p-1)} 出发，沿着势函数梯度的负方向走一步，步长为 \gamma \beta。直观理解是：让当前图像朝着更符合先验假设的方向移动。\gamma 是算法步长，\beta 是正则化权重，\nabla U(\theta^{(p-1)}) 是势函数在当前点的梯度。

第二步是近端映射（后向步骤）：

这一步在最大化似然函数的同时，要求解不能离第一步得到的 \theta_{\text{Reg}}^{(p)} 太远。第二项 \frac{1}{2\gamma}\|\theta - \theta_{\text{Reg}}^{(p)}\|^2 是近端惩罚项，\gamma 越小惩罚越强，解越接近 \theta_{\text{Reg}}^{(p)}。这个优化问题的解析解或高效数值解通常可以通过 EM 算法的若干次迭代获得。

两步交替的逻辑是：第一步利用先验知识对图像进行正则化处理，第二步将正则化后的图像拉回到与测量数据一致的方向。两者协同作用，最终收敛到 MAP 估计。

FBSEM-Net 的网络化实现

FBSEM-Net 将 FBS 算法中的正则化更新步骤替换为可学习的神经网络。具体而言，第一步变为：

其中 \mathcal{D} 是一个卷积神经网络，\vartheta 是其可学习参数。网络的作用相当于隐式地执行 \theta - \gamma\beta\nabla U(\theta) 这个操作，但 U 的具体形式由数据驱动学习得到，而非人工设计。第二步的近端映射保持不变，仍然通过 EM 更新来求解。

从网络结构图可以看到：初始图像 x^{(0)} 进入第一个块，每个块包含两条并行路径——上方是 EM 更新路径，利用测量数据 m 执行数据一致性步骤；下方是残差学习单元（Residual Learning Unit），这就是替代正则化梯度的 CNN，同时可以利用 MR 图像 x^{MR} 作为解剖先验。两条路径的输出通过融合（Fusion）模块合并，得到下一个块的输入。经过 N 个块后输出最终重建结果 x^{(N)}。

训练时，网络参数 \vartheta 和步长参数 \gamma 在整个展开网络中端到端学习。训练数据是 2 分钟采集的低计数数据，训练目标是使网络输出 \hat{\theta}_{\text{FBSEM-Net}} = \theta^{(P)} 尽可能接近 30 分钟采集的高质量参考图像 \hat{\theta}_{30\text{min-ML}}。

实验结果

图中展示了真实 FDG 数据的重建结果对比。第一列是 MRI 解剖图像作为参考，第二列是 30 分钟采集数据的参考重建。后续几列展示了不同方法对 2 分钟采集数据的重建：OSEM 重建噪声明显；单纯的 U-Net 后处理虽然降噪但可能丢失细节；FBSEM-Net 结合了 MRI 信息和数据一致性约束，在降噪的同时更好地保持了结构细节，视觉效果更接近 30 分钟参考图像。

合成方法：深度学习模型的隐式正则化

合成方法的核心思想

合成方法（Synthesis）采用与展开方法不同的策略。它同样保留发射断层成像的物理和统计模型 L(\theta|m)，但使用一个预训练的深度生成器 \mathcal{D}(x; \vartheta_{\text{FIX}}) 来实现正则化，其中 \vartheta_{\text{FIX}} 表示生成器的参数在重建过程中是固定的，不再更新。

正则化的约束是隐式的：要求解 \theta 必须属于生成器的像空间 \text{Im}\,\mathcal{D}。换句话说，重建得到的图像必须能够由生成器产生。如果生成器是在高质量医学图像上训练的，那么它的像空间就代表了"合理的医学图像"这一流形，噪声和伪影则不在这个流形上。通过限制解在这个流形上，自然就实现了正则化效果。

从网络结构示意图可以看到：左侧是某种编码表示，经过 U-Net 结构的生成器后得到右侧的重建图像。生成器的输入可以是潜在空间中的编码向量，输出是图像域中的重建结果。

生成器的训练

生成器采用 U-Net 架构的卷积神经网络，用 \mathcal{D}(x; \vartheta) 表示对 \theta 的一种表示。训练数据是配对的低计数输入图像和高计数目标图像：网络学习从噪声图像映射到干净图像。训练完成后，生成器的参数 \vartheta 被固定，其像空间就定义了重建问题的可行解集合。

与展开方法的区别在于：展开方法在重建过程中仍然训练或微调网络参数，网络与优化算法深度耦合；合成方法则将网络训练和重建过程分离，预训练好的生成器只提供一个约束流形，重建时只需要在这个流形上搜索最优解。

合成方法的数学形式化

优化问题的构建

合成方法的核心思想是用卷积神经网络来表示 PET 图像，即 \theta = \mathcal{D}(x; \vartheta)，其中 x 是网络的输入（潜在表示），\vartheta 是网络参数，\mathcal{D} 是网络映射。重建问题可以分两步理解：首先在潜在空间中找到最优的输入 \hat{x}，使得生成器输出的图像最大化似然函数：

然后通过生成器得到最终的重建图像：

这两步等价于求解一个带约束的优化问题：

约束条件 \theta = \mathcal{D}(x; \vartheta_{\text{FIX}}) 的含义是：重建图像 \theta 必须能够由生成器从某个潜在表示 x 产生出来，即 \theta 必须位于生成器的像空间 \text{Im}\,\mathcal{D} 中。

增广拉格朗日方法

带等式约束的优化问题可以通过增广拉格朗日方法转化为无约束问题。增广拉格朗日函数为：

这里 u 是对偶变量（拉格朗日乘子的缩放形式），\rho > 0 是惩罚参数。第二项 -\frac{\rho}{2}\|\theta - (\mathcal{D}(x; \vartheta_{\text{FIX}}) - u)\|_2^2 惩罚 \theta 与生成器输出之间的偏差，\rho 越大惩罚越强，约束满足程度越好。第三项 \frac{\rho}{2}\|u\|_2^2 是对偶变量的正则化项。当约束被精确满足时，\theta = \mathcal{D}(x; \vartheta_{\text{FIX}})，对偶变量 u 趋向于零。

ADMM 算法求解

增广拉格朗日问题通过交替方向乘子法（ADMM）求解。算法初始化潜在表示 x^{(0)}、对偶变量 u^{(0)}、惩罚参数 \rho 和迭代次数 K，然后在每次迭代 k = 0, \cdots, K-1 中执行三个步骤。

第一步更新图像 \theta：

这一步在最大化似然函数的同时，要求 \theta 接近当前生成器输出减去对偶变量的值。这是一个近端映射问题，可以通过 EM 算法的若干次迭代求解。

第二步更新潜在表示 x：

这一步寻找一个潜在表示，使得生成器输出尽可能接近 \theta^{(k+1)} - u^{(k)}。由于生成器 \mathcal{D} 是神经网络，这个优化问题可以通过 TensorFlow 或 PyTorch 的自动微分和梯度下降求解。

第三步更新对偶变量 u：

对偶变量累积了约束违反的历史：如果 \theta^{(k+1)} > \mathcal{D}(x^{(k+1)}; \vartheta_{\text{FIX}})，则 u 增大，下一次迭代会更强地惩罚这种偏差。经过 K 次迭代后返回 \theta^{(K)} 作为最终重建结果。

图中展示了不同方法的重建结果对比：高计数参考图像、EM 加滤波、MAP 公平惩罚、字典学习、CNN 去噪、以及迭代 CNN（即合成方法）。可以看到迭代 CNN 方法在保持结构细节的同时有效降低了噪声。

深度图像先验方法

DIP 方法的核心思想

深度图像先验（Deep Image Prior, DIP）方法同样保留发射断层成像的物理统计模型 L(\theta|m)，但使用无监督的深度图像先验框架 \mathcal{D}(x_{\text{FIX}}; \vartheta) 作为正则化手段。与前面的合成方法相比，DIP 方法有一个关键区别：合成方法中网络参数 \vartheta_{\text{FIX}} 是固定的（预训练好的），优化的是网络输入 x；而 DIP 方法中网络输入 x_{\text{FIX}} 是固定的（通常是随机噪声或 MRI/CT 图像），优化的是网络参数 \vartheta。

约束同样是隐式的，但形式不同：解 \theta 应该属于 \text{Im}\,\mathcal{D}(x_{\text{FIX}}; \vartheta)，这是由给定网络架构和固定输入所定义的图像族。DIP 方法的优势是不需要预训练：没有配对的训练数据也可以使用，只需要选择网络架构和固定输入（随机噪声或更好的是 MRI/CT 解剖图像）作为 CNN 输入来生成去噪图像。

DIP 训练与重建的统一

DIP 方法可以先单独训练网络参数。给定一个初始图像 \theta_0（比如传统方法的初步重建结果），DIP 训练寻找网络参数使输出尽可能接近这个初始图像：

然后用训练好的参数进行重建：

这等价于带约束的优化问题：

增广拉格朗日形式为：

ADMM-DIP 算法

DIP 方法的 ADMM 求解过程与合成方法类似，但优化变量不同。算法初始化网络参数 \vartheta^{(0)}（而非潜在表示）、对偶变量 u^{(0)}、惩罚参数 \rho 和迭代次数 K。

第一步更新图像 \theta：

这一步与合成方法相同，通过 EM 算法求解。

第二步更新网络参数 \vartheta：

这一步调整网络参数，使得在固定输入 x_{\text{FIX}} 下，网络输出尽可能接近目标 \theta^{(k+1)} - u^{(k)}。这是一个标准的神经网络训练问题，可以用 TensorFlow 或 PyTorch 求解。

第三步更新对偶变量 u：

图中展示了 DIP 方法的重建结果对比：先验图像 x_{\text{FIX}}（MRI）、高计数参考、EM 加滤波、KMRI、迭代 CNN、以及 DIPRecon。可以看到 DIPRecon 利用 MRI 作为结构先验，在保持解剖结构的同时实现了有效去噪，视觉效果接近高计数参考图像。

Plug-and-Play 方法

通用 ADMM 框架下的 MAP 重建

Plug-and-Play 方法的出发点是标准的 MAP 重建问题。设正则化函数为 R_{\mu}(\theta)，MAP 估计定义为：

其中 L(\theta|m) 是对数似然函数，R_{\mu}(\theta) 是正则化项，\mu 是正则化参数。这个问题可以通过 ADMM 算法求解，在每次迭代 k = 0, \ldots, K 中执行以下步骤。

第一步是重建更新：

这一步在最大化似然函数的同时，要求解 \theta 接近 z^{(k)} - u^{(k)}，其中 z 是辅助变量，u 是对偶变量，\rho^{(k)} 是惩罚参数。约束 \theta \in \mathbb{R}^N_+ 保证图像非负。

第二步是正则化更新：

这里 \text{Prox}_{R_{\mu}} 是正则化函数 R_{\mu} 的近端算子，定义为：

近端算子的作用是：在保持接近输入 \theta^k 的同时，最小化正则化函数。当 R_{\mu} 是某种平滑性惩罚时，近端算子的效果相当于对输入进行去噪处理，去除高斯噪声的同时保持图像的主要结构。

第三步是对偶变量更新：

用去噪器替代近端算子

Plug-and-Play 方法的核心思想是：既然近端算子 \text{Prox}_{R_{\mu}} 的作用相当于高斯去噪器，那么可以直接用现成的去噪算子 \mathcal{D}_{\vartheta} 来替代它，而无需显式定义正则化函数 R_{\mu}。这样，算法变为：

第一步重建更新保持不变：

第二步正则化更新替换为去噪：

第三步对偶变量更新保持不变：

去噪器 \mathcal{D}_{\vartheta} 可以是任何现成的深度学习去噪网络，通过离线的监督学习预先训练好，然后直接插入到 ADMM 迭代中使用。这种模块化的设计使得去噪器的训练和重建算法的设计可以独立进行。

Plug-and-Play 的具体算法

完整的 ADMM Plug-and-Play 算法如下：初始化辅助变量 z^{(0)}、对偶变量 u^{(0)}、惩罚参数 \rho 和迭代次数 K。对于 k = 0, \cdots, K-1，依次执行重建步骤、去噪步骤和对偶更新步骤，最终返回 \theta^{(K)}。其中重建步骤可以通过 EM 算法的若干次迭代求解。

收敛性分析

Plug-and-Play 方法的收敛性是一个核心理论问题。可以将 ADMM 的一次完整迭代表示为一个算子 \mathcal{T}_{\vartheta,\rho}，即 v^{(k+1)} = \mathcal{T}_{\vartheta,\rho}(v^{(k)})，其中 v 包含了所有迭代变量。收敛性分析依赖于以下数学结构：

如果算子 \mathcal{T}_{\vartheta,\rho} 是 1/2-均值的（1/2-averaged），那么它存在不动点，并且根据 Krasnoselskii-Mann 定理，迭代序列会收敛到某个不动点，即 \mathcal{T}_{\vartheta,\rho}(v) = v。

为了保证 \mathcal{T}_{\vartheta,\rho} 是 1/2-均值的，需要定义算子 \mathcal{Q}_{\vartheta} = 2\mathcal{D}_{\vartheta} - \text{Id}。如果 \mathcal{Q}_{\vartheta} 是非扩张的（nonexpansive），即 \|\mathcal{Q}_{\vartheta}(x) - \mathcal{Q}_{\vartheta}(y)\| \leq \|x - y\| 对所有 x, y 成立，那么 \mathcal{T}_{\vartheta,\rho} 就是 1/2-均值的，进而 \mathcal{D}_{\vartheta} 是某个极大单调（Maximal Monotone, MM）算子的预解式。在这种条件下，根据 Douglas-Rachford 分裂算法的收敛性理论，迭代变量 z, \theta 会收敛到某个不动点 \theta^*_{\vartheta,\rho}。

需要注意的是：收敛到的解满足不动点条件，但不一定对应于某个显式目标函数的最小化。这与传统的 MAP 重建不同，后者的解是明确定义的目标函数的最优点。

确保收敛的学习策略

强制非扩张性

为了保证 Plug-and-Play 算法的收敛性，需要确保算子 \mathcal{Q}_{\vartheta} = 2\mathcal{D}_{\vartheta} - \text{Id} 是非扩张的。实现这一目标有两种主要技术手段。

第一种是对网络所有层进行谱归一化（Spectral Normalization）。谱归一化将每一层权重矩阵的谱范数（最大奇异值）限制为不超过 1，从而控制整个网络的 Lipschitz 常数。如果网络的 Lipschitz 常数不超过 1，则网络是非扩张的。

第二种是局部雅可比正则化（Local Jacobian Regularization）。在训练去噪网络时，除了标准的均方误差损失外，还加入一个正则化项来约束网络的局部行为。训练目标变为：

其中 \theta_b 是训练样本，\bar{\theta}_b 是对应的干净目标图像，\tilde{\theta}_b 是在 \theta_b 邻域内随机采样的点。第一项是标准的均方误差损失，衡量去噪效果；第二项是雅可比正则化，\nabla \mathcal{Q}_{\vartheta}(\tilde{\theta}_b) 是算子 \mathcal{Q}_{\vartheta} 在点 \tilde{\theta}_b 处的雅可比矩阵，\|\cdot\| 是其谱范数。这一项惩罚雅可比矩阵谱范数超过 1-\epsilon 的情况，从而在局部强制非扩张性。参数 \beta 控制正则化强度，\epsilon 是一个小的正数，提供一定的松弛。

实验结果对比

图中展示了在真实数据上植入模拟肿瘤后的重建结果对比。左上角 EM(LD) 是低剂量数据的 EM 重建，噪声严重；右上角 EM(HD) 是高剂量参考重建，显示了三个植入的肿瘤位置（Tumor 1、Tumor 2、Tumor 3）。

中间一列展示了三种方法：后处理方法（Post-processing）直接对低剂量 EM 重建结果进行网络去噪；FBSEM 是前面介绍的展开算法；MBIR(Fair) 是使用 Fair 惩罚函数的传统基于模型的迭代重建。

右侧的 PnP-Unet 是 Plug-and-Play 方法，使用 U-Net 作为去噪器。从对比可以看出，PnP-Unet 在降低噪声的同时较好地保持了肿瘤的可见性，尤其是对于较小或对比度较低的肿瘤，表现优于单纯的后处理方法，同时在某些细节上也优于传统的 MBIR 方法。

三种混合方法的对比

学习方式

展开方法（Unrolling）采用端到端学习，整个展开网络作为一个整体进行训练，网络参数和重建过程紧密耦合。合成方法（Synthesis）和 Plug-and-Play 方法都采用离线学习，即先独立训练好生成器或去噪器，然后在重建时直接使用，训练过程和重建过程相互独立。

优化过程中网络的角色

展开方法中，网络参数在整个展开结构中端到端优化，每个块的网络都参与反向传播。合成方法中，虽然网络参数是固定的，但优化过程需要对网络进行前向和反向传播来更新潜在表示或网络参数（DIP 情况），因此网络仍然参与优化计算。Plug-and-Play 方法中，网络仅作为一个黑盒去噪器被调用，优化过程不需要对网络求梯度，网络完全不参与优化的数学结构。

内存占用

展开方法的内存占用与展开块数 N_{unroll} 和网络参数量 N_{params} 的乘积成正比，即 \propto N_{unroll} N_{params}，因为端到端训练需要存储所有块的中间激活值用于反向传播。合成方法和 Plug-and-Play 方法的内存占用仅与网络参数量成正比，即 \propto N_{params}，因为只需要存储单个网络的参数和激活值。

收敛性保证

展开方法在实践中的收敛性依赖于展开块数 N_{unroll}，理论上没有严格的收敛性保证，块数的选择往往是经验性的。合成方法的收敛性目前仍是开放问题，缺乏完善的理论分析。Plug-and-Play 方法在满足前述非扩张性条件时具有不动点收敛性的理论保证，这是其相对于其他两种方法的理论优势。

发射断层成像重建方法总结

解析重建方法

解析重建技术基于线性模型，需要对原始数据进行预校正（包括随机符合校正、散射校正、衰减校正等）。这类方法不考虑数据的统计特性，但有助于理解逆问题的本质特征，包括方向性伪影的产生机制、高频分量的放大效应、切趾和截断的必要性，以及三维重建中的数据冗余性。这些理论认识为设计更好的重建算法提供了基础。

经典迭代重建方法

经典的迭代重建技术采用更接近真实情况的复杂模型，包括完整的物理采集模型和数据的统计特性（泊松分布）。这些模型的引入对重建结果的定性和定量质量都有显著影响。正则化通过手工设计的先验来引导解具有分片平滑的特性，可以结合额外的解剖信息（如 MRI 或 CT）来恢复结构之间丢失的高频细节。

深度学习重建方法

深度学习重建方案倾向于使用自适应模型，这些模型通常从学习数据库中训练得到。其代价是在数据稀缺或低剂量条件下可能出现不稳定性和幻觉。混合方法将基于模型的似然函数与学习得到的正则化算子相结合，通过保留数据一致性约束来确保更好的鲁棒性。对于 Plug-and-Play 方法，在适当条件下甚至可以证明不动点收敛性。