Administrator
发布于 2026-09-09 / 0 阅读
0
0

医学成像(四):发射断层成像(PET/SPECT)

本节内容聚焦于核医学中的发射成像模态。与MRI重建不同,发射断层成像面对的是泊松断层成像病态定量逆问题,这类问题的统计特性和物理建模方式与MRI有本质区别,因此需要专门的求解方法。

从方法论角度,课程涵盖医学成像逆问题的完整处理流程:首先是逆问题的表征,这需要建立采集过程的物理模型,其次是解空间的设计,通过先验约束(可以是隐式的或显式的,也可以通过学习获得)来限定可行解的范围,最后是针对具体医学应用场景的优化算法设计。

非侵入性体内成像模态

医学影像技术可分为三大类。

image-20260114145217204

第一类是放射学,包括X射线成像、磁共振成像和超声成像,这些技术主要提供结构信息。

image-20260114145245101

第二类是核医学,包括伽马射线成像和正电子成像,这类技术的特点是能够提供功能和分子层面的信息。

image-20260114145316449

第三类是混合成像,将上述两类技术结合,例如PET/CT(正电子与X射线结合)、PET/MR(正电子与磁共振结合)以及PET/US(正电子与超声结合),混合成像能够同时获取解剖信息和功能信息。

透射成像与发射成像的区别

X射线、伽马射线和正电子成像都基于同一物理原理:检测沿直线传播的高能光子。然而,根据辐射源的位置不同,可分为两种成像模式。

透射成像
image-20260114145428561

在透射成像中,X射线源位于人体外部,射线穿透人体后被对侧的探测器接收。探测器测量的是射线穿过组织后的衰减程度,不同组织对X射线的吸收系数不同,因此这种成像方式能够区分骨骼、软组织、空气等不同解剖结构,属于解剖成像。

发射成像
image-20260114145447522

在发射成像中,辐射源位于人体内部。具体做法是将放射性示踪剂注入体内,示踪剂会随血液循环分布到特定组织或参与特定代谢过程。示踪剂发射的伽马射线或正电子湮灭产生的光子从体内向外传播,被环绕人体的探测器捕获。由于示踪剂的分布反映的是生理功能或代谢活动,因此发射成像属于功能成像或分子成像。

平面成像与断层成像

根据成像维度的不同,可分为平面成像和断层成像两种方式。

平面成像

平面成像将三维的人体信息投影到二维平面上,典型例子包括普通X光片(放射摄影)和闪烁显像。

image-20260114145602199

1895年Wilhelm Röntgen拍摄了第一张X光照片,1953年G.L. Brownell和W.H. Sweet进行了早期的正电子成像研究。平面成像的根本局限在于无法提供深度信息,所有沿射线方向的结构都叠加在一起,无法分辨前后位置关系。

断层成像
image-20260114145641757

断层成像通过从多个角度采集数据,利用计算方法重建出三维体积图像,能够提供完整的空间定位信息。CT是X射线计算机断层成像的缩写,利用透射原理获取解剖信息。SPECT是单光子发射计算机断层成像的缩写,利用伽马射线发射原理。PET是正电子发射断层成像的缩写,利用正电子湮灭产生的成对光子进行成像。断层成像的核心技术挑战在于断层重建,即如何从投影数据中恢复出三维分布,这正是本课程后续章节的主要内容。

发射断层成像的分子成像特性

发射断层成像能够在活体内可视化并测量分子和细胞层面的生物过程,同时不干扰机体的正常功能。这种成像方式具有极高的灵敏度,能够检测皮摩尔(picomolar)到纳摩尔(nanomolar)范围内的物质浓度。这一灵敏度水平远超其他成像模态,使得研究者能够追踪极微量的生物分子。发射断层成像的核心应用在于靶向那些在疾病发生或进展过程中起关键作用的生化过程。

临床应用中有两个典型示例。

image-20260114145833830

第一个是 [^{18}\text{F}]-L-DOPA,这种示踪剂用于观察多巴胺的合成过程,主要应用于帕金森病的诊断,因为帕金森病患者的多巴胺能神经元会发生退化,通过观察多巴胺合成能力的下降可以评估病情。

image-20260114145843964

第二个是 [^{18}\text{F}]-FDG,用于观察葡萄糖代谢,主要应用于肿瘤学,因为肿瘤细胞通常具有异常旺盛的糖代谢(Warburg效应),在FDG图像上会呈现高摄取。

核医学成像的技术要素

分子成像原理

分子成像的核心思想是使用分子探针来靶向特定的生化过程。所谓分子探针,是指能够与特定生物分子结合或参与特定代谢途径的化学物质。在实际操作中,注射的示踪剂量极小,属于示踪量级,不会改变患者的生物状态。这一点与药物治疗有本质区别:药物需要达到一定浓度才能产生治疗效果,而示踪剂只需要足够被探测到即可,其浓度远低于能够产生任何生理效应的水平。

放射化学标记

放射化学的任务是将放射性核素标记到分子探针上,形成放射性示踪剂。标记过程必须满足一个关键条件:不能改变分子探针的生物分布特性。换言之,标记后的放射性示踪剂在体内的行为应当与未标记的分子探针完全一致,这样探测到的信号才能真实反映目标生化过程。

image-20260114145852632

FDG的全称是氟代脱氧葡萄糖(fluorodeoxyglucose),化学名为2-脱氧-2-[^{18}\text{F}]氟-D-葡萄糖,用于追踪葡萄糖代谢。

image-20260114145901660

另一个例子是FPCIT(商品名Ioflupane),化学名为 [^{123}\text{I}]N-\omega-氟丙基-2\beta-甲氧羰基-3\beta-(4-碘苯基)降托烷,用于显像多巴胺转运体(DAT),主要应用于帕金森病和相关运动障碍的诊断。

信号检测与成像

放射性示踪剂注射到体内后,会随血液循环分布到目标组织。示踪剂中的放射性核素发生衰变时会发射伽马射线,这些射线从体内向外传播,被放置在体外的探测器捕获。探测器接收到的伽马射线能量范围通常在100到500 keV之间。

image-20260114150612810

通过分析探测器收集到的数据,可以重建体内感兴趣参数的空间分布,最典型的就是放射性示踪剂的浓度分布。

image-20260114150631146

PET完整处理流程

PET成像涉及一个完整的处理链条。首先是放射性核素的生产,通常使用回旋加速器完成。常用的核素包括 ^{18}\text{F}(半衰期110分钟)和 ^{11}\text{C}(半衰期20分钟)。半衰期是指放射性活度衰减到初始值一半所需的时间,较短的半衰期意味着核素需要在生产后尽快使用。

image-20260114150846726

接下来是放射性示踪剂的合成,将生产出的放射性核素标记到分子探针上,例如合成 [^{18}\text{F}]FDG。合成完成后将示踪剂注射到患者体内,等待一定时间让示踪剂分布到目标组织后,使用PET扫描仪进行数据采集。探测器环绕患者排列,记录从体内发出的伽马射线。采集到的原始数据需要经过图像重建算法处理,生成三维的放射性浓度分布图像。最后是图像的读取、处理和定量分析,提取临床或科研所需的信息。

核衰变的物理机制

发射断层成像依赖两种不同的核衰变机制,由此产生了两种主要的成像模态。

SPECT与伽马衰变

SPECT(单光子发射计算断层成像)利用的是直接发射伽马射线的核素。在伽马衰变(\gamma decay)过程中,原子核从激发态跃迁到较低能态时直接发射一个伽马光子,能量通常在100到300 keV范围内。由于每次衰变只产生一个光子,因此称为单光子成像。

PET与正电子衰变

PET(正电子发射计算断层成像)利用的是发射正电子的核素。在 \beta^+ 衰变过程中,原子核内的一个质子转变为中子,同时发射一个正电子(\beta^+)和一个中微子。正电子是电子的反粒子,当它在组织中运动一小段距离后会与周围的电子发生湮灭反应。湮灭时,正电子和电子的质量完全转化为能量,产生两个沿相反方向飞出的伽马光子,每个光子的能量均为511 keV。这个能量值来源于质能方程 E = mc^2,电子的静止质量对应的能量恰好是511 keV。PET探测器通过同时检测这对方向相反的光子(符合探测)来确定湮灭事件发生的位置,这种机制赋予了PET比SPECT更高的空间分辨率和灵敏度。

定量成像的实现

发射断层成像的最终目标是重建注射到体内的放射性核素的放射性浓度空间分布。重建结果的单位是贝克勒尔每单位体积,记作 [\text{Bq/cm}^3],其中贝克勒尔(Bq)是放射性活度的国际单位,1 Bq表示每秒发生一次核衰变。这个浓度分布直接反映了标记探针在体内的分布情况,进而反映了目标生化过程的空间分布和强度。

实现真正的定量成像需要建立准确的正向模型,即理解和建模采集过程中活度浓度与测量数据之间关系的物理现象。这些物理现象包括光子在组织中的衰减、散射、探测器的响应特性等,只有准确地对这些因素进行校正,才能从探测器测量到的计数数据中恢复出真实的放射性浓度值。

发射断层成像的物理学基础

理解发射成像需要掌握放射性衰变的基本物理过程,这是整个成像技术的根基。

核素与核素图

原子核由质子和中子(统称核子)组成,用符号 ^A_Z X 表示一个核素,其中 X 是元素符号,Z 是原子序数(即质子数),A 是质量数(即质子数与中子数之和)。原子序数 Z 决定了元素的化学性质,因为它等于核外电子数,质量数 A 则决定了原子核的质量。

image-20260114151100880

核素图以质子数为纵轴、中子数为横轴绘制。图中黄色方格表示稳定核素,白色方格表示不稳定(放射性)核素。稳定核素在图中形成一条狭窄的稳定带,轻核的稳定带接近 N = Z 的对角线,即中子数与质子数大致相等,随着原子序数增大,稳定核素需要更多的中子来平衡质子间的库仑排斥力,稳定带逐渐偏离对角线。位于稳定带两侧的核素是不稳定的,会通过放射性衰变向稳定带靠近。

β衰变

β衰变是一种放射性衰变过程,衰变后生成的子核与母核具有相同的质量数 A(即核子总数不变),但原子序数 Z 发生改变,因此属于不同的化学元素。具有相同质量数但不同原子序数的核素称为同量异位素(isobar)。

β⁻衰变

\beta^- 衰变发生在中子过剩的核素中。在这个过程中,一个中子转变为一个质子,同时发射一个电子(e^-)和一个反电子中微子(\bar{\nu}_e)。衰变方程的一个例子是:

^{18}_7\text{N} \rightarrow {}^{18}_8\text{O} + e^- + \bar{\nu}_e
image-20260114151228164

在核素图上,\beta^- 衰变表现为向左上方移动:质子数增加1,中子数减少1,质量数保持不变。这使得核素从稳定带的右侧(中子过剩侧)向稳定带移动。

β⁺衰变

\beta^+ 衰变发生在质子过剩的核素中。在这个过程中,一个质子转变为一个中子,同时发射一个正电子(e^+)和一个电子中微子(\nu_e)。衰变方程的一个例子是:

^{18}_9\text{F} \rightarrow {}^{18}_8\text{O} + e^+ + \nu_e
image-20260114151239860

在核素图上,\beta^+ 衰变表现为向右下方移动:质子数减少1,中子数增加1,质量数保持不变。这使得核素从稳定带的左侧(质子过剩侧)向稳定带移动。\beta^+ 衰变是PET成像的物理基础,因为衰变产生的正电子是成像信号的来源。

正电子的性质

正电子 e^+ 是电子的反粒子。反粒子与对应的粒子遵循相同的物理规律,具有相同的质量,但携带相反的电荷。电子带负电荷,正电子则带正电荷。

正电子与电子相遇时会发生湮灭反应,两个粒子的全部质量转化为能量,以光子的形式释放。根据质能方程,电子的静止质量对应的能量为:

E = m_e c^2 = 511 \text{ keV}

其中 m_e = 9.1 \times 10^{-31} kg 是电子的静止质量,c = 299\,792\,458 m/s 是真空中的光速,能量单位的换算关系为 1 \text{ eV} = 1.602 \times 10^{-19} J。

当正电子在组织中运动并逐渐减速至接近静止时,与电子发生湮灭,产生一对能量均为511 keV的伽马光子。由于动量守恒,这两个光子沿几乎完全相反的方向飞出(理想情况下成180°角)。PET探测器正是利用这种成对光子的符合探测来定位湮灭事件的发生位置。

β⁺衰变与电子俘获的竞争

将质子转变为中子的过程并非只有 \beta^+ 衰变一种途径,还存在电子俘获这一竞争机制。对于同一个母核,这两种衰变模式可能同时存在,各自占有一定的分支比。

\beta^+ 衰变发射正电子,衰变方程为:

^{18}_9\text{F} \rightarrow {}^{18}_8\text{O} + e^+ + \nu_e + 0.63 \text{ MeV} \quad (97\%)

电子俘获过程中,原子核俘获一个内层轨道电子(通常是K层电子),与一个质子结合生成中子,同时发射一个中微子。这个过程不产生正电子,衰变方程为:

^{18}_9\text{F} + e^- \rightarrow {}^{18}_8\text{O} + \nu_e + 1.65 \text{ MeV} \quad (3\%)

对于 ^{18}\text{F}\beta^+ 衰变的分支比约为97%,电子俘获的分支比约为3%。由于只有 \beta^+ 衰变产生的正电子才能在湮灭后产生可供探测的511 keV光子对,因此只有 \beta^+ 衰变对PET成像有贡献。电子俘获虽然也消耗放射性核素,但不产生可用于成像的信号,这在进行定量分析时需要考虑。

正电子发射核素的选择条件

只有 \beta^+ 衰变对PET成像有用,但并非所有质子过剩的放射性核素都能发生 \beta^+ 衰变。\beta^+ 衰变需要满足一个能量阈值条件:母核与子核之间的质量差必须大于两个电子质量(即1.022 MeV),因为衰变过程需要产生一个正电子,而正电子的静止质量与电子相同。如果能量不足,核素只能通过电子俘获方式衰变,无法用于PET成像。

image-20260114151330215

在核素图中,能够发生 \beta^+ 衰变的核素用青色标记,同时还需要考虑半衰期 T_{1/2} 大于1分钟的实用性条件(用斜线标记)。半衰期过短的核素在生产和运输过程中会大量衰变,难以实际应用。

常用PET核素的衰变特性

选择PET成像核素时,需要综合考虑 \beta^+ 衰变分支比 I_{\beta^+} 和电子俘获分支比 I_{EC},以及半衰期。理想的PET核素应具有高 \beta^+ 衰变分支比。

核素 I_{\beta^+} [%] I_{EC} [%] 半衰期 [min]
^{11}_6\text{C} 99.76 0.24 20
^{13}_7\text{N} 100 0 10
^{15}_8\text{O} 100 0 2
^{18}_9\text{F} 96.73 3.27 110
^{82}_{37}\text{Rb} 95.45 4.55 1
^{124}_{53}\text{I} 22.8 77.2 6013

^{13}\text{N}^{15}\text{O}\beta^+ 衰变分支比为100%,但半衰期很短,需要就近的回旋加速器才能使用。^{18}\text{F} 具有较长的半衰期(110分钟),允许在生产后运输到较远的成像中心,是目前最常用的PET核素。^{124}\text{I} 虽然半衰期很长(约4天),但 \beta^+ 衰变分支比仅为22.8%,大部分衰变通过电子俘获进行,成像效率较低。

正电子的能谱特性

\beta^+ 衰变释放的能量在正电子和中微子之间分配。以 ^{18}\text{F} 为例,衰变方程为:

^{18}_9\text{F} \rightarrow {}^{18}_8\text{O} + e^+ + \nu_e + 0.63 \text{ MeV}

由于能量在正电子和中微子之间随机分配,正电子的动能呈连续谱分布,而非单一能量值。能谱从零开始,到最大动能(等于衰变释放的总能量减去正电子静止质量能量)截止。不同核素由于衰变释放的总能量不同,其正电子能谱的形状和最大能量也不同。

image-20260114151349692

从能谱图可以看出,^{18}\text{F} 的正电子能量较低(最大约600 keV),而 ^{15}\text{O} 的正电子能量较高(最大约1700 keV)。

正电子的热化过程

正电子发射后携带一定的动能,在组织中运动时会通过一系列电离碰撞逐渐损失能量,这个过程称为热化。正电子在每次碰撞中使周围原子电离,自身动能减少,运动轨迹呈曲折的路径。当动能降至接近热能水平时,正电子与周围电子发生湮灭。

image-20260114151436603

正电子从发射点到湮灭点之间存在一定距离,称为正电子射程。射程的大小取决于正电子的初始动能:动能越高,射程越长。不同核素的平均正电子动能和在水中的平均射程如下:

核素 平均正电子动能 [keV] 水中平均射程 [mm]
^{18}_9\text{F} 250 0.6
^{11}_6\text{C} 386 1.1
^{13}_7\text{N} 492 1.5
^{15}_8\text{O} 735 2.5
^{82}_{37}\text{Rb} 1418 4.7

^{18}\text{F} 的正电子射程最短(约0.6 mm),^{82}\text{Rb} 的射程最长(约4.7 mm)。正电子射程对PET图像的空间分辨率有直接影响,因为探测到的湮灭位置与放射性核素的实际位置之间存在这个距离的偏差。

正电子湮灭与光子产生

正电子在热化后与电子发生湮灭反应。湮灭时两个粒子的质量完全转化为能量,静止状态下的总能量为:

2 \times m_e c^2 = 2 \times 511 \text{ keV}

根据能量守恒和动量守恒定律,湮灭产生两个光子,每个光子能量为511 keV,两光子沿相反方向飞出,夹角为180°。实际上由于湮灭时正电子和电子的残余动量,两光子的夹角会略微偏离180°,这种偏差称为非共线性,约为±0.25°,也会影响PET的空间分辨率。

正电子的探测原理

正电子在体内湮灭,无法直接被探测。真正能够从体内逃逸并被外部探测器捕获的是湮灭产生的511 keV光子对。探测器环绕人体放置,当两个探测器在极短时间窗口内(通常为纳秒量级)同时检测到511 keV光子时,系统判定这是一次符合事件,并推断湮灭发生在连接这两个探测器的直线(称为响应线,Line of Response,LOR)上的某处。

image-20260114151452775

正电子湮灭分布与空间分辨率

PET重建的是湮灭分布,而非放射性核素的空间分布。由于正电子射程的存在,湮灭位置与核素位置之间有一定偏移,因此重建图像相对于真实核素分布存在空间模糊。正电子平均动能越高的核素,射程越长,导致的空间分辨率退化越严重。

为了获得较高的空间分辨率,应优先选择平均正电子能量低于500 keV的核素。^{18}\text{F}(250 keV)、^{11}\text{C}(386 keV)和 ^{13}\text{N}(492 keV)都满足这一条件,而 ^{82}\text{Rb}(1418 keV)的正电子射程达4.7 mm,会显著降低图像分辨率。

SPECT中的伽马发射

SPECT成像探测的是核素直接发射的伽马光子,而非湮灭光子。在SPECT使用的核素衰变过程中,可能伴随 \beta 粒子(e^-e^+)的发射,但这些带电粒子只会在体内沉积能量、增加辐射剂量,对成像没有贡献。理想的SPECT核素应该是纯伽马发射体,即不伴随粒子发射,以最小化患者的辐射剂量。

亚稳态同质异能跃迁

纯伽马发射可以通过亚稳态同质异能跃迁实现。^{99m}\text{Tc}(锝-99m)是最常用的SPECT核素,其中m表示亚稳态(metastable)。^{99m}\text{Tc}^{99}\text{Mo}(钼-99)通过 \beta^- 衰变产生,^{99}\text{Mo} 的半衰期为2.7479天,\beta^- 衰变分支比为100%。

^{99m}\text{Tc} 处于激发态,半衰期为6.0067小时。它通过同质异能跃迁(Isomeric Transition,I.T.)释放能量,跃迁到基态 ^{99}\text{Tc},发射能量为142.683 keV的伽马光子。这个过程不涉及粒子发射,100%为同质异能跃迁。基态 ^{99}\text{Tc} 的半衰期极长(211.5×10³年),可视为准稳定核素。衰变能量参数为 Q^- = 1357.2 keV(^{99}\text{Mo}^{99m}\text{Tc}\beta^- 衰变)和 Q_{IT} = 142.683 keV(^{99m}\text{Tc} 的同质异能跃迁)。

纯伽马发射的另一种途径

除了亚稳态同质异能跃迁外,当 \beta^+ 衰变因能量不足而被禁止时,核素只能通过电子俘获衰变,这也能实现近似的纯伽马发射。^{111}_{49}\text{In}(铟-111)是这类核素的典型代表。

^{111}\text{In} 的半衰期为2.8049天,100%通过电子俘获衰变为 ^{111}_{48}\text{Cd}(镉-111),衰变能量 Q^+ = 861.8 keV。电子俘获后,子核 ^{111}\text{Cd} 处于激发态,通过级联伽马跃迁释放能量。主要的跃迁路径为:从416.6 keV能级跃迁到245.4 keV能级(释放171 keV伽马光子),再从245.4 keV能级跃迁到基态(释放245 keV伽马光子)。99.995%的衰变遵循这一主要路径,仅有0.005%通过其他次要路径。由于电子俘获过程本身不产生带电粒子发射,只伴随特征X射线和俄歇电子,^{111}\text{In} 的辐射剂量主要来自伽马光子,适合SPECT成像。

SPECT中的伽马探测

在SPECT成像中,探测器检测的是核素直接发射的伽马光子或特征X射线。与PET不同,SPECT探测的光子能量取决于所使用的放射性核素,而非固定的511 keV。常用SPECT核素的光子能量和半衰期如下:

核素 光子能量 [keV] 半衰期 [h]
^{99m}_{43}\text{Tc} 141 (\gamma) 6
^{111}_{49}\text{In} 245 & 171 (\gamma) 67
^{123}_{53}\text{I} 159 (\gamma) 13
^{201}_{81}\text{Tl} 71 & 69 & 80 (X) 73

^{99m}\text{Tc} 发射141 keV的伽马光子,^{111}\text{In} 发射245 keV和171 keV两种能量的伽马光子,^{123}\text{I} 发射159 keV的伽马光子,而 ^{201}\text{Tl} 主要发射的是特征X射线(71、69、80 keV)而非伽马射线。探测器需要根据所用核素调整能量窗口,选择性地接收特定能量范围内的光子。

准直技术

要从探测到的光子重建放射性核素的空间分布图像,必须知道每个探测到的光子来自哪个方向,这就需要准直技术。准直有两种实现方式。

机械准直

机械准直使用放置在探测器前方的铅准直器,准直器上有许多平行的小孔。只有沿孔轴方向传播的光子才能穿过准直器到达探测器,其他方向的光子被铅壁吸收。这种方式的优点是结构简单,适用于单光子成像,缺点是灵敏度损失严重,因为绝大多数光子被准直器阻挡,只有极小比例能够通过小孔被探测到。SPECT使用机械准直,探测效率通常低于0.01%。

电子准直

电子准直是PET特有的技术,利用符合探测电路实现。当两个探测器在极短时间窗口 \tau 内(|t_1 - t_2| < \tau,通常为几纳秒)同时检测到511 keV光子时,系统判定这两个光子来自同一次湮灭事件。湮灭点必定位于连接这两个探测器的直线上,这条直线称为响应线(Line of Response,LOR)。电子准直不需要物理准直器,因此不会阻挡光子,灵敏度远高于机械准直。

扫描仪几何结构

与平面成像不同,断层成像需要从 2\pi 范围内的所有方向采集数据,才能重建横断面图像。这意味着需要检测从被检体各个方向发出的光子。

SPECT扫描仪
image-20260114151836489

SPECT通常使用旋转式双头伽马相机。两个探测头相对放置(通常成180°或90°),整个机架绕患者旋转一周,从不同角度采集投影数据。旋转过程中,探测器分别在横断面和矢状面方向移动,覆盖完整的角度范围。双头设计相比单头可以将采集时间减半,或在相同时间内获得更多计数。

PET扫描仪
image-20260114151931619

PET使用完整的探测器环,探测器单元排列成一个或多个环形,完全包围患者。这种设计在横断面方向覆盖完整的 2\pi 角度,无需旋转即可同时采集所有方向的符合事件。探测器环的设计使PET具有比SPECT更高的灵敏度和时间分辨率。

伽马光子探测器的工作原理

标准的伽马光子探测器由致密闪烁体(无机晶体)与光电倍增管耦合组成,探测过程分为两个步骤。

image-20260114152016887

第一步是闪烁转换:高能伽马光子(约100 keV量级)进入闪烁体晶体后,与晶体原子相互作用,将能量沉积在晶体中。晶体被激发后发出大量低能可见光光子(约eV量级),这个过程称为闪烁。产生的闪烁光子数量与入射伽马光子的能量成正比。

第二步是光电转换与信号放大:闪烁光子到达光电倍增管的光阴极,通过光电效应释放出电子。这些电子在倍增管内经过多级倍增极(打拿极)的级联放大,每经过一级倍增极,电子数量增加数倍。最终在阳极收集到的电子形成电脉冲信号,脉冲幅度与入射伽马光子的能量成正比。通过测量脉冲幅度,可以确定光子能量,这对于能量甄别和散射校正至关重要。

SPECT伽马相机结构

SPECT使用的标准相机称为Anger相机,以发明者Hal Anger命名。Anger相机的核心部件是一块大面积的单块NaI(Tl)闪烁晶体(碘化钠掺铊),典型尺寸为40 cm × 40 cm。

image-20260114152531009

相机的结构从下到上依次为:平行孔准直器位于最前端,只允许垂直入射的光子通过,紧接着是NaI(Tl)闪烁晶体,将伽马光子转换为闪烁光,晶体后方是光电倍增管(PMT)阵列,将闪烁光转换为电信号并放大。

光子在晶体中产生闪烁后,闪烁光会被多个相邻的PMT同时探测到。Anger逻辑电路(脉冲算术电路)根据各PMT信号的加权计算,确定闪烁事件在晶体中的二维位置坐标 (s_1, s_2),同时计算信号总和得到能量信息 E。这三个参数(s_1 坐标、s_2 坐标、能量 E)经过模数转换(ADC)后送入计算机进行数据存储和图像重建。

SPECT伽马相机的实际安装

双头E.CAM相机是典型的SPECT扫描设备。

image-20260114152430059

从设备照片可以看到两个大型探测头相对安装在旋转机架上,患者躺在检查床上,床可以移动进入两个探测头之间的成像区域。探测头可以绕患者旋转,从不同角度采集数据。

SPECT数据采集过程

在SPECT采集中,探测器头绕患者旋转,在每个角度位置 \psi 处采集投影数据。对于每个角度 \psi,探测到的光子计数存储在二维直方图 m_\psi(s_1, s_2) 中,其中 s_1s_2 是探测器平面上的二维坐标。

image-20260114152612561

直方图 m_\psi(s_1, s_2) 中每个位置的数值代表在采样探测器位置 (s_1, s_2) 处、能量 E 落在预设能量窗口内的探测光子数。能量窗口的设置是为了选择特定能量的光子(如 ^{99m}\text{Tc} 的141 keV),同时排除散射光子。m_\psi 实际上对应于角度 \psi 处放射性示踪剂分布的投影图像。

随着探测器旋转,从 \psi = 0°\psi = 180°(或更大范围)采集一系列投影图像。每个角度的投影图像反映了从该方向观察时放射性分布的累积效果。将所有角度的投影数据综合起来,就可以通过断层重建算法恢复三维的放射性分布。

PET块状探测器结构

PET扫描仪的标准探测器采用块状探测器设计,由多个探测器环组成。每个块状探测器包含一个体素化的LSO(硅酸镥)晶体矩阵,典型的单个晶体尺寸为 4 mm × 4 mm。晶体矩阵由 2 × 2 排列的光电倍增管(PMT)读出。

与SPECT的Anger相机类似,PET块状探测器也使用Anger逻辑来确定闪烁事件在晶体矩阵中的位置。设四个PMT的信号分别为 P1P2P3P4,则位置坐标的计算公式为:

s_1 = \frac{(P2 + P3) - (P1 + P2)}{P1 + P2 + P3 + P4}
s_2 = \frac{(P1 + P2) - (P3 + P4)}{P1 + P2 + P3 + P4}

这两个公式通过比较相对两侧PMT信号的差异来确定事件发生的位置。s_1 坐标由水平方向两侧信号的差值决定,s_2 坐标由垂直方向两侧信号的差值决定,分母为总信号,起归一化作用。通过这种方式,可以将闪烁事件定位到具体的晶体单元。

image-20260114152717119

数字化PET探测器

最新一代的PET系统采用半导体光电探测器——硅光电倍增管(Silicon Photomultiplier,SiPM)取代传统的PMT。SiPM具有多项优势:响应速度快,可实现更精确的时间测量,体积紧凑,允许更高的探测器密度,支持晶体与光电探测器的一对一耦合,提高空间分辨率,对磁场不敏感,与MRI兼容,可用于PET/MR混合成像系统。

image-20260114152759587

采用SiPM技术的商用PET/CT系统包括GE的Discovery MI、Philips的Vereos和Siemens的Biograph Vision等。这些系统相比传统PMT系统具有更高的时间分辨率和空间分辨率。

PET数据采集模式

列表模式数据

PET采集中,每个符合事件被单独记录,形成列表模式数据。每个事件的记录包含:探测到两个光子的探测器标识 d_1r_1(第一个探测器的环号和晶体号)和 d_2r_2(第二个探测器的环号和晶体号),以及两个光子到达的时间差 t_1 - t_2

image-20260114152906255
有效事件的判定条件

一个符合事件被判定为有效需要满足两个条件。第一个条件是时间符合:两个光子的探测时间差必须在符合时间窗口内,即 |t_1 - t_2| \leq \tau。时间窗口 \tau 的设置取决于视野(Field of View)的大小,计算公式为:

\tau = \frac{R_{\text{FoV}}}{c/2}

其中 R_{\text{FoV}} 是视野半径,c 是光速。对于 R_{\text{FoV}} = 30 cm的典型视野,\tau \approx 2 ns。这个时间窗口确保只有来自视野内的湮灭事件才被接受。

第二个条件是能量甄别:每个光子(称为单事件,single)的能量必须落在以511 keV为中心的能量窗口内。这用于排除散射光子和其他背景事件。

每个有效事件由两个探测器的标识 dr,以及飞行时间值 t_2 - t_1 来表征。飞行时间信息可用于进一步缩小湮灭点的定位范围。

正弦图数据格式

在图像重建时,列表模式数据通常被整理成投影矩阵,称为正弦图(sinogram)。一组平行的响应线对应发射分布的一个投影。例如,角度 \psi = 45° 的所有平行LOR构成该角度的投影数据。将各角度的投影按行排列,就形成正弦图。正弦图中,点源在图像中呈现正弦曲线形状,这也是正弦图名称的由来。

image-20260114152928685

坐标 s_1 表示LOR在投影方向上的位置(即投影坐标),s_2 表示轴向位置。不同角度 \psi 的投影数据组合起来,包含了重建三维图像所需的完整信息。

发射断层成像的正向模型

发射断层成像的核心问题是建立活度浓度与探测事件之间的关系,即正向模型。这个模型对图像重建至关重要,因为重建算法需要知道给定的活度分布会产生什么样的测量数据,才能通过逆向推断从测量数据恢复活度分布。

采集模型的基本原理

当检测到一个有效事件时,可以确定某处沿投影线(SPECT)或响应线(PET)发生了一次衰变。然而,如果暂时不考虑PET的飞行时间信息,仅凭符合探测无法确定湮灭(或衰变)发生在这条线上的具体深度位置。这是断层成像需要从多角度采集数据的根本原因:单一方向的投影无法区分沿射线方向的深度信息。

image-20260114153002788

采集模型的基本表述是:对于给定的一条线(投影线或响应线),探测到的事件数与沿该线的发射分布成正比。这意味着,如果沿某条线的放射性活度越高,在这条线上探测到的计数就越多。这个线性关系是断层重建数学框架的基础。

线积分采集模型

采集模型的数学表述是:对于给定的一条线 L,探测到的事件数 y_L 与发射分布 f(\mathbf{r}) 沿该线的线积分成正比:

y_L \approx \int_L f(\mathbf{r}) \mathrm{d}\mathbf{r}

这里 f(\mathbf{r}) 表示空间位置 \mathbf{r} 处的放射性活度浓度,积分沿着响应线(或投影线)L 进行。这个模型的物理意义是:沿一条线方向的所有放射性衰变事件都会贡献到这条线对应的测量计数中,贡献的权重与该位置的活度成正比。

image-20260114153026598

线积分模型是断层成像的数学基础。它表明测量数据是活度分布的线性变换(具体是Radon变换),从多个角度采集足够多的线积分数据后,理论上可以通过逆变换唯一地重建出原始的活度分布 f(\mathbf{r})

PET飞行时间技术原理

在传统PET中,虽然符合探测确定了湮灭发生在响应线上,但无法确定沿线的具体位置。飞行时间(Time of Flight,ToF)技术利用两个光子到达探测器的时间差来进一步定位湮灭点。

设两个光子分别在时刻 t_1t_2 被探测到,则湮灭点相对于响应线中点的偏移量 \Delta l 为:

\Delta l = \frac{1}{2} \times (t_2 - t_1) \times c

其中 c 是光速。这个公式的推导基于简单的运动学:如果湮灭点偏离中点距离 \Delta l,则两个光子到达各自探测器的路程差为 2\Delta l,对应的时间差为 2\Delta l / c,因此 \Delta l = (t_2 - t_1) \times c / 2

image-20260114153053243
飞行时间定位的不确定性

飞行时间定位的精度受限于符合时间分辨率 \delta t,由此导致的空间定位不确定性 \delta l 为:

\delta l = \frac{1}{2} \times c \times \delta t

要实现4 mm的空间分辨率,需要约26.8 ps的符合时间分辨率。目前的技术水平尚未达到这一要求,最先进的系统约为200 ps,对应约3 cm的空间不确定性。因此,飞行时间信息虽然能缩小湮灭点的定位范围,但仍不足以直接定位,图像重建方法仍然是必需的。

飞行时间采集模型

考虑飞行时间信息后,采集模型需要修改。对于给定的响应线 L 和时间差 t_2 - t_1,探测到的事件数 y_{L, t_2 - t_1} 为:

y_{L, t_2 - t_1} \approx \int_L h(r - (t_2 - t_1)c/2) \, f(\mathbf{r}) \mathrm{d}\mathbf{r}

这里 h(x) 是飞行时间核函数,通常假设为高斯分布,其宽度反映时间分辨率。核函数 h 以位置 (t_2 - t_1)c/2 为中心,对发射分布 f(\mathbf{r}) 进行加权。靠近飞行时间估计位置的区域权重较大,远离的区域权重较小。

image-20260114153111990

与传统的线积分模型相比,飞行时间模型不再是简单的均匀线积分,而是加权线积分。权重由飞行时间核函数决定,使得测量数据更多地反映湮灭点附近的活度信息,而非整条响应线上的总活度。

飞行时间技术的优势

飞行时间技术的主要优势是降低重建图像中的统计噪声,从而提高信噪比(Signal-to-Noise Ratio,SNR)。飞行时间PET与传统PET的信噪比之比近似为:

\frac{SNR_{\text{ToF}}}{SNR_{\text{non-ToF}}} \simeq \sqrt{\frac{D_{\text{obj}}}{\delta l}}

其中 D_{\text{obj}} 是被扫描对象的直径,\delta l 是飞行时间定位的空间不确定性。这个公式表明,飞行时间带来的信噪比增益取决于对象尺寸与定位精度的比值。对象越大(如腹部扫描)或时间分辨率越好(\delta l 越小),增益越显著。

image-20260114153120593

对于直径40 cm的大型对象,当时间分辨率从800 ps改善到100 ps时,信噪比增益可从约2倍提高到约6倍。这意味着相同的图像质量可以用更少的扫描时间或更低的注射剂量获得,或者相同条件下可以获得更好的图像质量。

飞行时间技术的发展历程

飞行时间PET的性能持续改善。2006年商用系统的时间分辨率约为600 ps,2014年提升到400 ps,2023年最先进的系统(如Siemens Biograph Vision)达到约178 ps。图中标注的商用系统包括GE Discovery MI PET/CT和Siemens Biograph Vision。

image-20260114153137189

未来的挑战目标是达到10 ps的时间分辨率,这将使 \delta l 降至约1.5 mm,接近晶体本身的空间分辨率,届时飞行时间信息将能够直接定位湮灭点,从根本上改变重建算法的需求。

线积分模型的局限性

线积分数学模型过于简化,不足以准确描述采集过程中实际发生的物理现象,尤其在需要定量估计时问题更为突出。该模型未考虑三个关键因素:第一,光子不沿直线传播的事件,第二,未被探测到的事件,第三,底层过程的统计特性。

背景事件的影响

理想情况下,有用事件应该是沿直线传播的光子。然而实际采集中,并非所有探测到的事件都是有用事件,存在两类主要的背景事件。

散射事件

光子在穿过组织时可能与原子发生康普顿散射,改变传播方向。散射后的光子如果仍被探测到,其响应线不再通过原始的衰变位置,导致错误的空间定位。散射事件在SPECT和PET中都会发生。

随机符合

随机符合是PET特有的问题。当两个独立的湮灭事件各产生一个光子,这两个无关的光子恰好在符合时间窗口内被探测到时,系统会错误地将它们判定为来自同一次湮灭的光子对。随机符合产生的响应线与真实的湮灭位置完全无关。

image-20260114153216270

考虑背景事件后,采集模型修正为:

y_L \approx \int_L f(\mathbf{r}) \mathrm{d}\mathbf{r} + b_L

其中 b_L 表示响应线 L 上的背景事件贡献,在实际处理中通常假设 b_L 可以通过独立的方法估计并视为已知量。

光子衰减效应

光子穿过物质时会被吸收或散射,导致强度衰减。根据Beer-Lambert定律,初始光子数为 N_0 的光束穿过衰减系数为 \mu、厚度为 d 的均匀物质后,剩余光子数为:

N = N_0 e^{-\mu \cdot d}

在人体内部区域,由于光子需要穿过较长的组织路径,衰减非常严重,超过90%的发射光子可能被衰减而无法到达探测器。

SPECT中的衰减校正
image-20260114153244618

在SPECT中,从位置 \mathbf{r} 发射的光子需要穿过从 \mathbf{r} 到探测器的路径才能被探测,衰减因子取决于这段路径上的组织衰减系数。对于投影线 L(s, \psi),采集模型为:

y_L^{\text{SPECT}} \approx \int_{L(s,\psi)} f(\mathbf{r}) e^{-\int_{L(s,\psi,r)} \mu(\mathbf{r}') \mathrm{d}\mathbf{r}'} \mathrm{d}\mathbf{r} + b_L

这里积分号内的指数因子表示从发射点 \mathbf{r} 到探测器路径上的累积衰减。\mu(\mathbf{r}') 是位置 \mathbf{r}' 处对应140 keV光子的衰减系数 \mu_{140\text{keV}}(\mathbf{r})

PET中的衰减校正

在PET中,湮灭产生的两个光子必须都到达探测器才能形成符合事件。因此衰减因子是整条响应线 L(a,b) 上的总衰减,与湮灭发生的具体位置无关。采集模型为:

y_L^{\text{PET}} \approx e^{-\int_{L(a,b)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}} \int_{L(a,b)} f(\mathbf{r}) \mathrm{d}\mathbf{r} + b_L

这里 \mu(\mathbf{r}) 是511 keV光子的衰减系数 \mu_{511\text{keV}}(\mathbf{r})。PET衰减校正的一个优势是衰减因子可以提到积分外面,因为它只依赖于响应线而不依赖于线上的具体位置。

探测器效率因素

伽马探测器的效率不是100%,存在多种效率损失机制:闪烁晶体的有限阻止能力使部分光子穿透而不被探测,探测器单元之间存在间隙,电子学系统存在死时间,在处理一个事件期间无法记录新事件。这些因素用探测效率乘性因子 \epsilon_L 来建模。

加入探测效率后,完整的采集模型为:

SPECT:

y_L^{\text{SPECT}} \approx \epsilon_L \int_{L(s,\psi)} f(\mathbf{r}) e^{-\int_{L(s,\psi,r)} \mu(\mathbf{r}') \mathrm{d}\mathbf{r}'} \mathrm{d}\mathbf{r} + b_L

PET:

y_L^{\text{PET}} \approx \epsilon_L e^{-\int_{L(a,b)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}} \int_{L(a,b)} f(\mathbf{r}) \mathrm{d}\mathbf{r} + b_L

探测效率 \epsilon_L 对于不同的响应线可能不同,取决于探测器单元的位置、入射角度等因素,通常通过归一化扫描或计算方法确定。

数据的统计特性

采集模型表明探测事件数与沿线的发射分布成正比,但需要进一步明确发射分布的具体含义。在成像检查过程中,存在多个随机过程:注射的放射性核素总数按体积分布、核素的衰变、光子的探测,这些都是随机过程。建立完整的采集模型需要对这些随机过程进行数学建模。

放射性核素分布的泊松过程模型

基本假设

对注射的放射性示踪剂作两个假设。第一,如果重复进行相同的实验,注射的放射性核素总数服从泊松分布。第二,在任意时刻,每个注射核素的空间位置是独立同分布的随机变量,其分布由示踪剂的生理学和化学特性决定,与注射顺序无关。

空间泊松点过程

基于上述假设,可以证明:在任意时刻,按体积索引的注射放射性核素数量构成一个空间泊松点过程。具体而言,对于不相交的体积单元(如体素),各体积内的核素数量是相互独立的泊松随机变量。位置 \mathbf{r} 处体素内核素数量的均值为 \tilde{f}(\mathbf{r}, t),它是时间 t 和空间位置 \mathbf{r} 的函数,反映了示踪剂在体内的动态分布。

探测过程的非齐次空间泊松过程模型

衰变与探测的假设

在放射性核素分布的泊松过程模型基础上,还需要对衰变和探测过程作进一步假设。第一,每个核素的衰变时间是独立的随机变量,不依赖于空间位置,服从均值为 \lambda 的指数分布。这里 \lambda 是平均寿命,与半衰期 T_{1/2} 的关系为 \lambda = T_{1/2} / \ln 2。指数分布的无记忆性意味着核素在任意时刻衰变的概率只取决于从当前时刻起经过的时间,而与核素已经存在了多久无关。第二,假设理想探测器:每次衰变被归属到单一的响应线(LOR),且衰变事件的探测相互独立,不存在堆积效应或死时间影响。

探测率的表达式

基于上述假设,可以证明:按发射位置和时间索引的探测事件数构成一个非齐次空间泊松点过程。该过程的强度(探测率)为:

\lambda(\mathbf{r}, t) \propto \tilde{f}(\mathbf{r}, t) \frac{\exp(-t/\lambda)}{\lambda}

这个表达式反映了两个因素的乘积:\tilde{f}(\mathbf{r}, t) 是位置 \mathbf{r} 处核素数量的期望,\exp(-t/\lambda)/\lambda 是指数分布的概率密度函数,表示核素在时刻 t 衰变的概率密度。

探测数据的泊松噪声模型

实际采集中,数据在注射后从时刻 t_ft_f + \Delta t 的时间段内测量。定义空间泊松过程的均值函数为:

f(\mathbf{r}) = \int_{t_f}^{t_f + \Delta t} \lambda(\mathbf{r}, t) \mathrm{d}t

这个积分将时变的探测率在采集时间段内累积,得到与空间位置相关的期望探测数。基于此,探测事件数可建模为泊松分布,这一模型在实验上足够精确:

SPECT:

y_L^{\text{SPECT}} \sim \text{Po}\left( \epsilon_L \int_{L(s,\psi)} f(\mathbf{r}) e^{-\int_{L(s,\psi,r)} \mu(\mathbf{r}') \mathrm{d}\mathbf{r}'} \mathrm{d}\mathbf{r} + b_L \right)

PET:

y_L^{\text{PET}} \sim \text{Po}\left( \epsilon_L e^{-\int_{L(a,b)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}} \int_{L(a,b)} f(\mathbf{r}) \mathrm{d}\mathbf{r} + b_L \right)

这里 \text{Po}(\cdot) 表示以括号内表达式为均值的泊松分布。实际应用中,还需要额外的乘性标定因子(如衰变常数、分支比等)来建立探测计数与示踪剂浓度分布之间的定量关系。

泊松噪声的统计特性

对于泊松随机变量 y_L,其概率质量函数为:

p(y_L = k) = \frac{\bar{y}_L^k e^{-\bar{y}_L}}{k!}, \quad k \in \mathbb{N}_0

其中 \bar{y}_L 是均值参数。泊松分布有一个关键特性:期望等于方差:

\text{E}(y_L) = \text{Var}(y_L) = \bar{y}_L

由此可得信噪比(SNR)与均值的平方根成正比:

\text{SNR} \propto \frac{\bar{y}_L}{\sqrt{\bar{y}_L}} = \sqrt{\bar{y}_L}

这意味着要使信噪比提高一倍,需要使计数增加四倍。在发射断层成像的典型条件下,每个投影单元平均只探测到 \mathcal{O}(1) 个事件,即个位数量级,这导致数据的统计噪声极为严重。

image-20260114153303475

[^{18}\text{F}]FDG脑部PET研究的图像可以直观看到:3600秒采集时间获得的图像质量良好,但随着采集时间缩短到120秒、60秒、10秒,图像噪声急剧增加,最终几乎无法辨识结构。这与MRI形成鲜明对比:MRI图像基本不受采集时间影响而保持清晰。

噪声等效计数率

并非所有探测到的事件都是有用的真符合事件,背景事件(散射和随机符合)同样会贡献到重建图像的噪声中。为了更好地表征数据中的有效信息量,引入噪声等效计数率(Noise Equivalent Count Rate,NECR)的概念。NECR定义为:有用符合事件(不含随机和散射)的信噪比与单位时间NECR的平方根成正比。

NECR = \frac{T}{1 + kR/T + S/T}

其中 T 是真符合计数率,R 是随机符合计数率,S 是散射计数率。参数 k 取决于随机符合估计方法:使用单计数率估计时 k = 1,使用延迟符合窗口估计时 k = 2

NECR可以理解为等效的无背景真符合计数率,即假想情况下如果没有散射和随机符合,为达到相同信噪比所需的真符合计数率。研究表明,在使用经典重建方法时,全局或局部图像信噪比与NECR或NEC衍生量具有良好的相关性。

NECR曲线的特性

以PET/MR系统为例,计数率曲线随活度浓度变化的典型特性如下。随着活度浓度增加,总计数(Prompts)、真符合(Trues)、随机符合(Randoms)和散射(Scatter)都会增加。然而,随机符合的增长速度比真符合更快(随机符合与活度的平方成正比,而真符合与活度成正比),因此在高活度时随机符合会成为主要噪声源。

image-20260114153328088

NECR曲线(NEC)呈现先上升后下降或趋于平稳的特征。存在一个最优活度浓度使NECR最大化,超过这个浓度后,继续增加活度反而会降低有效信息量。这为临床实践中选择合适的注射剂量提供了指导。

空间分辨率的限制

由于采集数据是离散的,重建得到的是原始活度分布 f(r) 与某个低通滤波器 W_c 的卷积:

W_c * f(r)

实际重建的是离散化图像。MRI的空间分辨率可以达到亚毫米级(<1 mm),而 [^{18}\text{F}]FDG PET脑部研究的典型分辨率约为4 mm。

image-20260114153410319

这个分辨率差异的根本原因在于探测灵敏度:如果灵敏度不提高,单纯提高空间分辨率会使每个体素的计数更少,噪声更大,图像质量反而下降。因此,空间分辨率与灵敏度之间存在权衡关系。

提高信噪比的技术途径

提高探测器灵敏度

探测器的阻止能力决定了多少光子能被有效探测。LSO闪烁晶体密度为7.4 g/cm³,晶体厚度通常为20至30 mm,以确保足够的光子阻止效率。

扩大轴向覆盖范围是另一个提高灵敏度的途径。传统PET的轴向视野约25 cm,而新一代全身PET(total-body PET)的轴向覆盖可达194 cm,能够同时采集全身数据,大幅提高总体灵敏度。

飞行时间分辨率的改善也能有效提高信噪比。技术发展路线为:178 ps(对应2.7 cm空间不确定性)→ 100 ps(对应1.5 cm)→ 目标10 ps(对应1.5 mm)。

先进数据处理技术

除硬件改进外,软件算法的进步同样能提升图像质量。先进的图像重建算法能够更好地利用统计模型,抑制噪声同时保持分辨率。机器学习方法(如深度学习去噪)也越来越多地应用于PET图像后处理。

PET图像质量的历史演进

[^{18}\text{F}]FDG全身和全身PET的图像质量从1996年到2019年有了显著提升。1996年的系统需要60分钟采集时间才能获得可接受的图像质量,2016年的系统将采集时间缩短到10分钟,同时图像质量有所改善,2019年的全身PET系统仅需14分钟即可完成全身扫描,图像质量进一步提高。这种进步来自探测器技术、数据采集方法和重建算法的综合改进。

image-20260114153441601

发射断层成像的临床应用

本节进入原理部分的第二个主题:发射断层成像的应用。以下主要介绍PET在临床各领域的具体应用场景。

FDG示踪剂的代谢原理

[^{18}\text{F}]FDG是PET中使用最广泛的放射性示踪剂。FDG是葡萄糖的类似物,用 ^{18}_9\text{F} 标记。理解FDG成像的关键在于其代谢机制与正常葡萄糖的差异。

正常葡萄糖的代谢路径为:血浆中的葡萄糖通过葡萄糖转运蛋白(GLUT)进入组织细胞,在细胞内被磷酸化生成葡萄糖-6-磷酸,继而转化为果糖-6-磷酸,最终通过糖酵解途径分解为丙酮酸。

image-20260114153509736

FDG的代谢路径在前两步与葡萄糖相同:通过GLUT进入细胞,被磷酸化生成FDG-6-磷酸。然而,FDG-6-磷酸无法继续进行后续的代谢步骤,因为FDG分子结构中2位碳上的羟基被氟原子取代,使其不能被进一步代谢。同时,FDG-6-磷酸带有负电荷,无法穿过细胞膜逃逸。结果是FDG被"困"在细胞内,持续积累。

这种代谢陷阱机制使得FDG的摄取量直接反映组织的葡萄糖代谢水平:代谢越旺盛的组织,摄取并积累的FDG越多,在PET图像上显示越亮。

肿瘤学应用

肿瘤的代谢特征

恶性肿瘤的一个显著特征是葡萄糖代谢增强(Warburg效应)。肿瘤细胞即使在有氧条件下也倾向于通过糖酵解而非氧化磷酸化来获取能量,导致对葡萄糖的需求大幅增加。这使得FDG PET成为检测和评估肿瘤的理想工具。

临床应用场景

PET FDG主要用于癌症分期和治疗效果监测。对于脑肿瘤,由于正常脑组织本身葡萄糖代谢很高,有时需要使用其他示踪剂。成像通常采用静态方式,在注射后约30分钟(示踪剂达到平衡时)进行扫描。

image-20260114153529405

从临床图像可以看到:基线扫描(治疗前)显示多处高摄取灶,提示肿瘤病变,治疗后随访扫描显示这些高摄取灶明显减少或消失,表明治疗有效。图像的定量信息(单位为kBq/cm³)允许对治疗反应进行客观评估。

标准化摄取值

为了比较不同患者、不同时间点的FDG摄取,引入标准化摄取值(Standardized Uptake Value,SUV)作为半定量指标:

SUV = \frac{\text{放射性浓度 [kBq/mL]}}{\text{注射活度 [kBq]} / \text{体重 [g]}} \times 1 \text{[g/mL]}

SUV实际上是组织中放射性浓度与假设活度均匀分布于全身时的平均浓度之比。SUV = 1表示该组织的摄取与全身平均水平相当,SUV > 1表示摄取高于平均,提示代谢活跃。

image-20260114153543129

临床实践中常用SUV阈值来区分良恶性病变,例如SUV阈值 = 5意味着SUV超过5的区域高度怀疑为恶性。图像显示肝脏病变区域的SUV最大值为14.84,远超阈值,强烈提示恶性。

神经学应用

神经退行性疾病的代谢成像

[^{18}\text{F}]FDG PET可以显示脑区的代谢活动。在神经退行性疾病中,受累脑区的神经元功能下降,表现为局部葡萄糖代谢降低。这种代谢降低模式是神经退行性变的生物标志物。

image-20260114153609330

从阿尔茨海默病的图像可以看到疾病进展过程:正常人脑代谢均匀,轻度认知障碍患者的顶叶和颞叶开始出现代谢降低,阿尔茨海默病患者这些区域的代谢明显减低,范围扩大。这种特征性的代谢降低模式有助于早期诊断和鉴别诊断。

多巴胺能系统成像

[^{18}\text{F}]F-DOPA是多巴胺神经递质的前体,用于显像纹状体的多巴胺能神经末梢。在帕金森病中,黑质多巴胺能神经元进行性退化,导致纹状体多巴胺能神经末梢减少。

image-20260114153617473

从F-DOPA图像可以看到:正常人的纹状体(包括尾状核和壳核)摄取均匀,接近100%,早期帕金森病患者壳核后部摄取开始下降,晚期帕金森病患者整个纹状体摄取严重减低。这种成像可用于帕金森病的早期诊断和病情监测。

癫痫灶定位

对于药物难治性癫痫患者,[^{18}\text{F}]FDG PET用于定位引起癫痫活动的脑区,为手术治疗提供依据。在发作间期,癫痫灶通常表现为局部代谢降低(低摄取区)。

image-20260114153650720

图像比较显示:健康志愿者脑代谢对称均匀,癫痫患者可见局部代谢降低区域,与MRI(T1加权和T2-FLAIR)上的结构异常相对应。PET与MRI的融合图像([^{18}\text{F}]FDG/T1w)有助于精确定位癫痫灶。

image-20260114153705121

心脏病学应用

心肌活力评估

PET在心脏病学中的主要应用是灌注/代谢成像,用于区分存活但功能受损的心肌与不可逆损伤(坏死/瘢痕)的心肌组织。这对于心肌梗死后评估受损心肌的活力、决定是否进行血运重建治疗至关重要。

代谢成像

[^{18}\text{F}]FDG用于评估心肌代谢。存活心肌仍有葡萄糖代谢能力,会摄取FDG,坏死心肌丧失代谢功能,不摄取FDG。

image-20260114153721476

图像显示:健康志愿者心肌呈均匀的环形摄取,心肌梗死患者可见摄取缺损区,对应坏死心肌。

灌注成像

心肌灌注成像使用 [^{15}\text{O}]\text{H}_2\text{O}(水)或 [^{13}\text{N}]\text{NH}_3(氨)等示踪剂,评估心肌血流。灌注-代谢匹配(两者均降低)提示不可逆损伤,灌注-代谢不匹配(灌注降低但代谢保留)提示冬眠心肌,血运重建后可能恢复功能。

动态成像与定量分析

静态成像的局限

临床成像通常采用静态方式,在示踪剂分布达到平衡后进行单时间点扫描。这种方法简单快捷,但只能提供示踪剂累积量的信息,无法获得动力学参数。

动态成像的优势

在动态成像中,从示踪剂注射开始连续采集多个时间帧的数据,记录示踪剂在体内分布的时间变化过程。通过分析示踪剂的动力学特性,可以估计定量的生理参数。

image-20260114153743350

[^{11}\text{C}]PE2I(多巴胺转运体配体)为例,图像展示了注射后0-1分钟、4-5分钟、13-15分钟、30-35分钟、55-60分钟的脑部分布变化。早期主要反映血流灌注,后期反映特异性结合。

通过绘制不同脑区的时间-活度曲线(TAC),可以看到纹状体(富含多巴胺转运体)的摄取持续上升并保持高水平,而小脑(参考区,多巴胺转运体密度低)的摄取先升后降。利用动力学模型分析这些曲线,可以提取结合势(binding potential)等定量生理参数,比单纯的SUV更能反映受体/转运体密度。

药代动力学房室模型

动态PET成像引入了一个新的(时间维度上的)逆问题:利用描述测量动态数据与生理参数之间关系的数学模型,从数据中估计感兴趣的生理参数。

一般房室模型结构

房室模型将示踪剂在体内的分布划分为若干个动力学上均匀的"房室"。一个典型的完整模型包含:动脉血浆中的示踪剂 X,通过血流进入组织,组织中的游离示踪剂 X(free compartment),与受体或转运体结合的示踪剂 X(bound compartment),代谢产物 X^m(metabolized compartment),以及流出到静脉血的示踪剂。

image-20260114153856368

各房室之间的转运由速率常数描述:K_1 是从动脉血浆到组织的转运速率(单位为mL/min/g),k_2 是从组织返回血浆的速率,k_3 是从游离态到结合态或代谢态的速率,k_4 是从结合态返回游离态的速率,k_5 是代谢清除速率。

组织中游离示踪剂浓度 C_f 的变化由以下微分方程描述:

\frac{\mathrm{d}C_f}{\mathrm{d}t} = K_1 C_a - k_2 C_f - k_3 C_f + k_4 C_b - k_5 C_f

其中 C_a 是动脉血浆浓度(输入函数),C_b 是结合态浓度。逆问题的目标是从测量数据中估计这些微观参数 k

葡萄糖代谢的房室模型

FDG的动力学可用不可逆两房室模型描述,用于估计局部葡萄糖代谢率。

image-20260114153925321

模型结构为:血浆中的FDG以速率 K_1 进入组织,以速率 k_2 返回血浆,组织中的游离FDG以速率 k_3 被磷酸化为FDG-6-磷酸。由于FDG-6-磷酸不能继续代谢也不能逆转,磷酸化步骤是不可逆的(k_4 = 0)。

PET测量的是整个组织体素内的总放射性浓度,包括血管内和血管外的成分。完整的模型方程为:

C_{\text{PET}}(t) = (1 - V_b) C_a(t) * \left( \frac{K_1 k_2}{k_2 + k_3} e^{-(k_2 + k_3)t} + \frac{K_1 k_3}{k_2 + k_3} \right) + V_b \cdot C_a(t)

其中 V_b 是血管容积分数,* 表示卷积运算,C_a(t) 是动脉输入函数。

宏观参数与代谢率

在实际应用中,通常关注宏观参数而非所有微观参数。FDG的净摄取速率(influx rate)定义为:

K_i = \frac{K_1 k_3}{k_2 + k_3}

K_i 反映了FDG从血浆不可逆进入组织(被磷酸化)的总速率。由于FDG与葡萄糖竞争相同的转运和磷酸化通路,局部脑葡萄糖代谢率(regional Cerebral Metabolic Rate of Glucose consumption,rCMRGLc)与 K_i 成正比:

\text{rCMRGLc} \propto K_i

精确计算rCMRGLc还需要考虑FDG与葡萄糖在转运和磷酸化效率上的差异(集总常数,lumped constant)以及血浆葡萄糖浓度。

参数成像

参数成像实现了从半定量成像(SUV)到全定量成像的跨越。

image-20260114153945758

以头颈部肿瘤的 [^{18}\text{F}]FDG PET扫描为例,SUV图像显示的是放射性浓度的归一化值,而 K_i 图像(单位为min⁻¹)显示的是葡萄糖代谢速率的定量估计值。两种图像虽然都能显示肿瘤,但 K_i 提供的是真正的生理参数,具有更明确的生物学意义。

动态成像的挑战

基于动态成像的参数估计对数据噪声更为敏感。每个时间帧的计数较少,噪声较大,而参数估计需要拟合多帧数据的时间曲线,噪声会在拟合过程中传播和放大。

为了应对这一挑战,可以将动力学参数估计嵌入到图像重建过程中,形成动态参数重建(dynamic parametric reconstruction)方法。这种方法直接从原始投影数据估计参数图像,而非先重建各帧图像再逐体素拟合,能够更好地利用数据的统计特性,获得噪声更低的参数图像。

断层重建概述

课程现在进入第三部分:断层重建。本节首先介绍解析重建方法。

发射断层数据的特殊性

发射断层成像数据具有几个独特的特征,使其重建问题比其他成像模态更具挑战性。

第一,采集数据并非直接得到感兴趣的量。原始数据是投影形式,受到衰减、背景事件、PET中的正电子射程等物理因素的影响,且服从泊松统计。第二,目标是获得定量的功能参数估计,不仅仅是对比度图像。第三,这是一个大规模问题:正弦图数据的规模在 \mathcal{O}(100M, 10B)(即一亿到百亿)量级的数据单元,需要估计 \mathcal{O}(10M, 100M)(即千万到亿)量级的参数。

因此,重建的本质是:求解一个具有泊松数据的大规模定量断层病态逆问题。

正问题与逆问题的关系

重建必须基于采集过程的数学模型。

image-20260114154053272
正问题

正问题(Direct Problem)描述的是:给定未知的放射性分布,通过采集模型预测会得到什么样的测量数据。采集模型包括了之前讨论的所有物理因素:投影几何、衰减、散射、探测器响应、统计噪声等。

逆问题

逆问题(Inverse Problem)是正问题的反向:给定采集到的测量数据,通过重建算法估计未知的放射性分布。重建算法的设计必须与采集模型相匹配,对采集物理的建模越准确,重建结果的定量精度越高。

这种正问题-逆问题的框架是所有断层成像重建的基础。不同的重建算法本质上是对这个逆问题的不同求解策略:解析方法基于连续数学模型的精确逆变换,迭代方法基于离散模型的优化求解。

重建模型的设计原则

一个好的重建模型需要平衡两个相互矛盾的要求。

第一个要求是模型足够简单,以支持快速重建。由于断层重建是大规模计算问题,过于复杂的模型会导致计算代价过高,无法在临床可接受的时间内完成重建。

第二个要求是模型足够精确,以实现准确的定量。精确的采集模型需要描述两方面的物理过程:一是采集过程中的物理现象,包括光子与患者组织的相互作用(衰减、散射)以及探测器的几何和响应特性,二是可能还需要包含示踪剂的药代动力学(用于参数成像)。

改进模型可以显著提升重建后的图像质量和定量精度。这是一个持续的研究方向:在计算效率允许的范围内,不断改进采集模型的精确度。

解析重建方法

图像重建算法可分为知识驱动方法和数据驱动方法两大类。知识驱动方法的特点是采集模型 A 是显式已知并构建的。

解析方法的基本思想

解析方法基于简化的采集模型,寻找其数学上的解析逆变换。如果采集模型为 A,测量数据为 g,则重建图像 f 通过以下公式直接计算:

f = A^{-1}(g)

这里 A^{-1} 表示采集模型的逆算子。解析方法的优点是计算速度快、算法确定性强,缺点是只能处理简化的采集模型(如纯线积分模型),难以融入复杂的物理校正。

PET中的解析重建
image-20260114154207925

在PET解析重建中,正问题采用线积分模型:假设测量数据是活度分布沿响应线的线积分。逆问题则是对线积分模型进行解析求逆。经典的滤波反投影(Filtered Back-Projection,FBP)算法就是基于Radon变换逆公式的解析重建方法。

迭代重建方法

当采集模型过于复杂而无法直接求逆,或者解析逆变换不存在时,需要采用迭代方法。

迭代方法的基本框架

迭代重建的核心思想是:从一个初始估计出发,通过反复迭代逐步改进估计值,直到估计值产生的模拟数据与实际测量数据足够接近。

具体流程如下:首先给定初始估计 f^{[0]},然后在第 p 次迭代中,利用采集模型(此时作为模拟器)计算当前估计 f^{[p]} 对应的模拟采集数据 A(f^{[p]}),将模拟数据与实际采集数据 g 进行比较,计算差异,根据差异信息更新估计值,得到 f^{[p+1]},重复上述过程直到收敛。

image-20260114154231159

更新公式的一般形式为:

f^{[p+1]} = f^{[p]} \times \Delta f^{[p]}

其中 \Delta f^{[p]} 是基于模拟数据与实际数据比较结果计算的校正因子。不同的迭代算法(如MLEM、OSEM等)对应不同的校正因子计算方式。

迭代方法的优势

迭代方法的优势在于可以融入任意复杂的采集模型,包括精确的几何建模、衰减校正、散射校正、探测器响应建模、统计噪声模型等。这使得迭代重建能够实现更高的定量精度,是目前临床PET重建的主流方法。

数据驱动方法

数据驱动方法的特点是采集模型 A 不是显式已知的,而是通过数据学习得到。

直接深度学习

直接深度学习方法使用深度神经网络从训练数据库中学习逆模型。网络输入是测量数据 g,输出是重建图像 f

f = \text{neural network}(g, \text{training database})

网络可以从零开始设计,也可以采用特定的架构(如U-Net等)。训练过程需要大量的配对数据(投影数据-真实图像对)。这种方法的优点是推理速度快,缺点是可解释性差,且对训练数据分布之外的情况泛化能力有限。

混合方法

混合方法结合知识驱动模型和数据驱动模型的优势。

物理信息深度学习

物理信息深度学习(Physics-informed Deep Learning)是一种典型的混合方法。它保留可信的物理模型和噪声分布模型作为采集过程的描述,同时使用数据驱动模型(如深度神经网络)来进行正则化或先验建模。

这种方法的优点是:物理模型部分保证了基本的物理正确性和可解释性,数据驱动部分可以学习复杂的图像先验,提供更好的降噪和细节保持能力。物理信息深度学习是当前医学图像重建领域的研究热点。

解析重建方法

本节详细介绍基于线积分模型的解析重建方法。解析重建是透射断层成像(CT)和发射断层成像(SPECT、PET)共同的数学基础。

线积分模型的数学框架

解析重建的数学基础是线积分模型。设重建的视野(Field of View)为 \mathbb{R}^n 中的紧支撑集 \mathfrak{C},其中 n = 2 对应二维重建,n = 3 对应三维重建。

待重建的函数 f(\mathbf{r}) 是定义在 \mathfrak{C} 上的连续函数,即 f(\mathbf{r}): \mathfrak{C} \subset \mathbb{R}^n \rightarrow \mathbb{R}。在发射断层成像中,f(\mathbf{r}) 表示放射性活度浓度分布。

测量数据 g(L) 是定义在投影线流形上的连续函数,即 g(L): \mathbb{T}^{(2n-2)} \rightarrow \mathbb{R},其中 \mathbb{T}^{(2n-2)} 是所有投影线构成的流形。对于二维情况(n=2),投影线流形是二维的(用角度和偏移量参数化),对于三维情况(n=3),投影线流形是四维的。

线积分模型将数据与图像联系起来:

g(L) = \int_L f(\mathbf{r}) \mathrm{d}\mathbf{r}

即测量数据等于待重建函数沿投影线 L 的线积分。这个变换称为Radon变换(或X射线变换),解析重建的目标就是求解Radon变换的逆。

透射断层成像的线积分模型

在CT中,待重建的物理量是光子线性衰减系数 \mu(\mathbf{r}),对应于 f(\mathbf{r})

采集几何与物理过程

CT扫描仪由旋转的X射线管和旋转的X射线探测器阵列组成。X射线管发出初始强度为 I_0 的射线束,穿过人体后被探测器接收,测量到的强度为 I。射线用两个参数描述:\beta 是X射线管的旋转角度,\alpha 是探测器单元相对于射线束中心的角度。投影线记为 L(\beta, \alpha)

image-20260114154519910

根据Beer-Lambert定律,射线穿过物质后的强度衰减满足:

I = I_0 e^{-\int_{L(\beta,\alpha)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}}
CT线积分模型

对上式取对数并取负,得到CT的线积分模型:

-\ln\left(\frac{I}{I_0}\right)_{\text{CT}} = g(L) = \int_{L(\beta,\alpha)} f(\mathbf{r}) \mathrm{d}\mathbf{r}

这里 f(\mathbf{r}) = \mu(\mathbf{r}) 是衰减系数分布。CT数据经过对数变换后,满足标准的线积分模型,可以直接应用Radon变换逆公式进行重建。

SPECT的线积分模型

在SPECT中,待重建的物理量是单光子发射浓度 f(\mathbf{r})

采集几何

SPECT使用旋转的伽马探测器头,探测器前配有平行孔准直器。探测器在角度 \psi 位置采集数据,投影线用偏移量 s 和角度 \psi 参数化,记为 L(s, \psi)

变量定义

f(\mathbf{r}) 是单光子发射浓度(待重建量),\mu(\mathbf{r}) 是伽马射线线性衰减系数(通常已知或可从CT获取),y 是探测到的伽马光子数。

SPECT衰减线积分模型

由于光子从发射点到探测器的路径上会被组织衰减,SPECT的采集模型不是纯线积分,而是衰减加权线积分:

y_{\text{SPECT}} = g(L) = \int_L f(\mathbf{r}) e^{-\int_{L(s,\psi,r)} \mu(\mathbf{r}') \mathrm{d}\mathbf{r}'} \mathrm{d}\mathbf{r}

指数因子中的积分 \int_{L(s,\psi,r)} \mu(\mathbf{r}') \mathrm{d}\mathbf{r}' 是从发射点 \mathbf{r} 沿投影线到探测器路径上的累积衰减。这个模型比CT复杂,因为衰减因子依赖于发射点在投影线上的位置。

image-20260114154630945

SPECT衰减线积分模型不满足标准Radon变换的形式,其解析逆变换更为复杂,通常称为衰减Radon变换或指数Radon变换。

PET的线积分模型

在PET中,待重建的物理量是正电子发射浓度 f(\mathbf{r})

采集几何

PET使用环形探测器阵列,响应线 L(a,b) 由两个符合探测器 ab 确定。

变量定义

f(\mathbf{r}) 是正电子发射浓度(待重建量),\mu(\mathbf{r}) 是511 keV光子线性衰减系数,y 是探测到的符合事件数。

PET原始采集模型

湮灭产生的两个光子都需要到达探测器才能形成符合事件,因此总衰减因子是整条响应线上的累积衰减:

y_{\text{PET}} = e^{-\int_{L(a,b)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}} \int_{L(a,b)} f(\mathbf{r}) \mathrm{d}\mathbf{r}
PET线积分模型

PET相比SPECT有一个关键优势:衰减因子可以提到积分外面,因为它只依赖于响应线本身,而不依赖于湮灭发生在线上的具体位置。这意味着如果衰减因子已知,可以通过乘以衰减因子的倒数进行预校正:

y_{\text{PET}} \cdot e^{+\int_{L(a,b)} \mu(\mathbf{r}) \mathrm{d}\mathbf{r}} = g(L) = \int_{L(a,b)} f(\mathbf{r}) \mathrm{d}\mathbf{r}

经过衰减校正后,PET数据满足标准的线积分模型,可以直接应用Radon变换逆公式。这是PET相对于SPECT在解析重建方面的一个显著优势。

image-20260114154813528

断层成像的几何结构

根据射线束的几何形态,断层成像可分为三种主要几何结构。

平行束几何

在平行束几何中,所有投影线相互平行。这种几何对应于使用平行孔准直器的PET和SPECT系统。平行束几何的数学处理最为简单,是经典Radon变换理论的直接应用场景。

image-20260114154840887
扇形束几何

在扇形束几何中,投影线从一个点源发出,呈扇形展开。这种几何对应于单层CT扫描仪,以及使用扇形束准直器的SPECT系统。扇形束数据可以通过重排算法转换为等效的平行束数据,也可以直接使用扇形束重建算法。

image-20260114154850494
锥形束几何

在锥形束几何中,投影线从一个点源发出,呈三维锥形展开。这种几何对应于多层CT扫描仪,以及使用针孔准直器的SPECT系统。锥形束重建是真正的三维问题,数学处理比平行束和扇形束复杂得多。

image-20260114154902386

平行束几何的坐标系统

以三维平行束几何为例,建立完整的坐标系统用于描述投影线和投影数据。

方向向量与投影平面

投影线的方向用单位向量 \mathbf{e}_\theta \in \mathbb{S}^2 表示,其中 \mathbb{S}^2 是三维空间中的单位球面。方向向量可以用两个角度参数化:方位角 \psi(azimuthal angle)和极角 \phi(co-polar angle)。

与方向向量 \mathbf{e}_\theta 正交的平面称为投影平面,定义为:

\theta^\perp = \{\mathbf{r} \in \mathbb{R}^3 | \mathbf{r} \cdot \mathbf{e}_\theta = 0\}

投影平面是所有与 \mathbf{e}_\theta 垂直的向量构成的二维子空间。

image-20260114155722682
投影线的参数化

给定方向 \mathbf{e}_\theta 和投影平面上的位置向量 \mathbf{s} \in \theta^\perp,投影线 L(\mathbf{e}_\theta, \mathbf{s}) 定义为:

L(\mathbf{e}_\theta, \mathbf{s}) = \{\mathbf{s} + t\mathbf{e}_\theta | t \in \mathbb{R}^1\}

即从点 \mathbf{s} 出发沿方向 \mathbf{e}_\theta 延伸的直线。

投影线流形

三维空间中所有投影线构成的集合(投影线流形)为:

\mathbb{T}^4 = \{(\mathbf{e}_\theta, \mathbf{s}) | \mathbf{e}_\theta \in \mathbb{S}^2, \mathbf{s} \in \theta^\perp\}

这是一个四维流形:方向向量在二维球面 \mathbb{S}^2 上变化(2个自由度),位置向量在二维投影平面 \theta^\perp 上变化(2个自由度)。

二维与三维成像的区别

二维成像

在二维成像模式下,所有测量的投影线都垂直于断层轴 \mathbf{e}_z。这意味着数据只包含横断面内的信息,不同横断面之间相互独立。

在PET中,二维模式通过在探测器环之间放置隔板(septa)实现,只接受同一环或相邻环之间的符合事件。这样,三维的重建问题被分解为一系列独立的二维问题,每个横断面可以单独重建。

image-20260114155804089

二维成像的数学特点是:待重建函数 f(\mathbf{r}) 被分解为一系列二维平面进行独立重建,所有方向 \mathbf{e}_\theta \in \mathbb{S}^1(圆周)都被测量,方向仅由角度 \psi 定义,对于每个方向 \mathbf{e}_\theta,所有投影位置(由径向位置 s 定义)都被测量。

三维成像

在三维成像模式下,测量的投影线不限于垂直于断层轴,可以有任意倾斜角度。这样可以利用更多的符合事件,大幅提高灵敏度。

image-20260114155812496

三维成像的数学特点是:待重建函数 f(\mathbf{r}) 需要整体重建,不能分解为独立的二维问题,并非所有方向 \mathbf{e}_\theta \in \mathbb{S}^2(球面)都被测量,存在采样不完整的问题,对于给定方向 \mathbf{e}_\theta,并非所有投影位置 \mathbf{s} 都必然被测量。

算法差异

二维和三维解析重建算法共享相同的核心概念(Radon变换及其逆),但具体实现本质上不同。二维重建算法成熟且计算效率高,三维重建需要处理数据不完整性,算法更为复杂,计算量也更大。

n维平行射线变换

定义

n 维平行射线变换(也称为X射线变换或Radon变换)将图像函数 f: \mathbb{R}^n \rightarrow \mathbb{R} 变换为定义在投影线流形上的函数 \mathbf{X}f。投影线流形为:

\mathbb{T}^{(2n-2)} = \{(\mathbf{e}_\theta, \mathbf{s}) | \mathbf{e}_\theta \in \mathbb{S}^{n-1}, \mathbf{s} \in \theta^\perp\}

平行射线变换的定义为:

\mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) = \int_{\mathbb{R}^1} f(\mathbf{s} + t\mathbf{e}_\theta) \mathrm{d}t

即沿方向 \mathbf{e}_\theta、经过点 \mathbf{s} 的直线上 f 的线积分。

数据采集模型

采集数据 g 与待重建图像 f 的关系通过平行射线变换建立:

g(\mathbf{e}_\theta, \mathbf{s}) = \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s})
重建问题

图像重建的数学本质是求解平行射线变换的逆:给定 g = \mathbf{X}f,求 f。这就是Radon变换逆问题,其解析解即为著名的滤波反投影公式。

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

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

\mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) = \int_{\mathbb{R}^1} f(\mathbf{s} + t\mathbf{e}_\theta) \mathrm{d}t
二维情况:一维投影

在二维情况下,图像函数 f(x, y) 定义在 \mathbb{R}^2 上,投影数据 g(\psi, s) 定义在二维流形 \mathbb{T}^2 上。投影方向由单一角度 \psi 确定,投影位置由标量 s 确定。

二维平行射线变换的显式表达式为:

\mathbf{X}f(\psi, s) = \int_{\mathbb{R}^1} f(s\cos\psi - t\sin\psi, s\sin\psi + t\cos\psi) \mathrm{d}t

这里 (s\cos\psi, s\sin\psi) 是投影平面上距原点为 s 的点,(-\sin\psi, \cos\psi) 是投影方向,积分沿该方向进行。

三维情况:二维投影

在三维情况下,图像函数 f(x, y, z) 定义在 \mathbb{R}^3 上,投影数据 g(\psi, \phi, s_1, s_2) 定义在四维流形 \mathbb{T}^4 上。投影方向由两个角度(方位角 \psi 和极角 \phi)确定,投影位置由二维向量 (s_1, s_2) 确定。

投影平面上的位置向量表示为:

\mathbf{s} = s_1 \mathbf{e}_{\theta_1}^\perp + s_2 \mathbf{e}_{\theta_2}^\perp

其中 \mathbf{e}_{\theta_1}^\perp\mathbf{e}_{\theta_2}^\perp 是投影平面内的两个正交基向量。

投影方向向量的显式表达式为:

\mathbf{e}_\theta = -\cos\phi\sin\psi \, \mathbf{e}_x + \cos\phi\cos\psi \, \mathbf{e}_y + \sin\phi \, \mathbf{e}_z

一维投影与正弦图

投影的几何意义

对于二维图像 f(x, y),在角度 \psi 处的投影 \mathbf{X}f(\psi, s) 是一维函数,表示沿该角度方向的所有线积分。从图示可以看到:\psi = 0° 时投影沿水平方向,\psi = 45° 时投影沿对角方向,\psi = 90° 时投影沿垂直方向,\psi = 135°\psi = 180° 分别对应其他方向。

image-20260114155831789
正弦图的构成

将所有角度的一维投影堆叠起来,形成二维数据阵列,称为正弦图(sinogram)。正弦图的定义为:

m(s, \psi) = g(\psi, s) = \mathbf{X}f(\psi, s)

正弦图以 s(投影位置)为横轴,\psi(投影角度)为纵轴。在正弦图中,图像中的一个点源会映射为一条正弦曲线,这是因为点源在不同角度投影中的位置随角度呈正弦变化。正弦图这一名称正是来源于此特性。

image-20260114155853979

从图示可以看到,左边的简单二维图像(包含一个小圆和一个大的椭圆形结构)经过Radon变换后,在右边形成了具有正弦曲线特征的正弦图。正弦图完整地包含了重建所需的全部信息。

n维中心切片定理

中心切片定理(Central Slice Theorem,也称为傅里叶切片定理或投影切片定理)是解析重建理论的核心。它建立了投影数据的傅里叶变换与图像傅里叶变换之间的关系。

傅里叶变换的定义

首先定义两种傅里叶变换。投影数据 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) 对投影平面变量 \mathbf{s}(n-1) 维傅里叶变换为:

\mathcal{F}_{n-1}\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \int_{\theta^\perp} e^{-2\pi i \mathbf{s} \cdot \boldsymbol{\nu}_\perp} \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) \mathrm{d}\mathbf{s}

其中 \boldsymbol{\nu}_\perp \in \theta^\perp 是投影平面内的频率向量。

图像函数 f(\mathbf{r})n 维傅里叶变换为:

\mathcal{F}_n\{f\}(\boldsymbol{\nu}) = \int_{\mathbb{R}^n} e^{-2\pi i \mathbf{r} \cdot \boldsymbol{\nu}} f(\mathbf{r}) \mathrm{d}\mathbf{r}

其中 \boldsymbol{\nu} \in \mathbb{R}^nn 维频率向量。

定理陈述

对于 \mathbb{R}^n 中的函数 f(\mathbf{r}),有:

\mathcal{F}_{n-1}\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \mathcal{F}_n\{f\}(\boldsymbol{\nu} = \boldsymbol{\nu}_\perp)

这个定理说明:对投影数据在投影平面内做傅里叶变换,得到的结果等于图像傅里叶变换在通过原点、垂直于投影方向的超平面上的取值。

定理证明

证明通过直接计算展开。从投影数据傅里叶变换的定义出发:

\mathcal{F}_{n-1}\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \int_{\theta^\perp} e^{-2\pi i \mathbf{s} \cdot \boldsymbol{\nu}_\perp} \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) \mathrm{d}\mathbf{s}

代入平行射线变换的定义:

= \int_{\theta^\perp} e^{-2\pi i \mathbf{s} \cdot \boldsymbol{\nu}_\perp} \int_{\mathbb{R}^1} f(\mathbf{s} + t\mathbf{e}_\theta) \mathrm{d}t \mathrm{d}\mathbf{s}

进行变量替换 \mathbf{r} = \mathbf{s} + t\mathbf{e}_\theta。由于 \boldsymbol{\nu}_\perp\mathbf{e}_\theta 正交,有 \mathbf{r} \cdot \boldsymbol{\nu}_\perp = \mathbf{s} \cdot \boldsymbol{\nu}_\perp(因为 t\mathbf{e}_\theta \cdot \boldsymbol{\nu}_\perp = 0)。变量替换后:

\mathcal{F}_{n-1}\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \int_{\mathbb{R}^n} e^{-2\pi i \mathbf{r} \cdot \boldsymbol{\nu}_\perp} f(\mathbf{r}) \mathrm{d}\mathbf{r} = \mathcal{F}_n\{f\}(\boldsymbol{\nu} = \boldsymbol{\nu}_\perp)

证毕。

二维中心切片定理的几何解释

在二维情况下,中心切片定理有直观的几何解释。对于图像 f(x, y),其二维傅里叶变换 \mathcal{F}_2\{f\}(\nu_x, \nu_y) 定义在二维频率平面上。

对于角度 \psi = 45° 处的一维投影 \mathbf{X}f(\psi = 45°, s),其一维傅里叶变换 \mathcal{F}_1\{\mathbf{X}f\}(\psi = 45°, \nu_s) 是一维函数。

image-20260114155918971

中心切片定理说明:这个一维傅里叶变换恰好等于二维傅里叶变换沿通过原点、与投影方向垂直的直线上的取值。即 \psi = 45° 投影的一维傅里叶变换对应于二维频率平面上倾斜 45° 的径向切片。

从所有角度采集投影,就可以获得二维傅里叶变换在所有通过原点的径向切片上的取值。理论上,如果采集了足够密集的角度,就可以完整地覆盖整个二维频率平面,从而通过逆傅里叶变换重建出原图像。这就是傅里叶重建方法的理论基础。

三维中心切片定理

三维中心切片定理是二维情况的自然推广。对于三维图像 f(x, y, z),其三维傅里叶变换 \mathcal{F}_3\{f\}(\nu_x, \nu_y, \nu_z) 定义在三维频率空间中。

对于给定方向 \mathbf{e}_\theta 的二维投影 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}),其二维傅里叶变换 \mathcal{F}_2\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_s) 等于三维傅里叶变换在通过原点、垂直于 \mathbf{e}_\theta 的平面上的取值。

image-20260114155934805

从图示可以看到:三维物体 f(x, y, z) 沿某方向投影得到二维投影图像 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}),对三维物体做三维傅里叶变换得到三维频谱 \mathcal{F}_3\{f\}(\nu_x, \nu_y, \nu_z),对二维投影做二维傅里叶变换得到 \mathcal{F}_2\{\mathbf{X}f\}(\mathbf{e}_\theta, \boldsymbol{\nu}_s),中心切片定理表明,后者恰好是前者在某个通过原点的平面上的切片。

直接傅里叶变换重建

中心切片定理的一个直接应用是:在二维情况下,如果投影 g(\psi, s) = \mathbf{X}f(\psi, s) 在角度区间 \psi \in [0, \pi[ 上被测量,则可以重建图像 f。这是因为 [0, \pi[ 范围的投影对应的频率切片已经完整覆盖了整个二维频率平面([\pi, 2\pi[ 范围的投影提供的是冗余信息)。

直接傅里叶变换重建流程

直接傅里叶变换重建(Direct Inversion Fourier Transform,DIFT)的步骤为:首先对每个角度的一维投影 g(\psi, s) 做一维FFT,得到 \mathcal{F}_1\{g\}(\psi, \nu_s),根据中心切片定理,这些数据位于二维频率平面的极坐标网格上,将极坐标网格上的数据重采样到笛卡尔网格上,最后做二维逆FFT得到重建图像 f(x, y)

image-20260114155957600
主要困难

DIFT方法的主要困难在于从极坐标网格到笛卡尔网格的重采样。极坐标采样在低频区域密集、在高频区域稀疏,重采样会引入插值误差,尤其是高频部分。这个问题限制了DIFT方法的实际应用,滤波反投影方法提供了更优雅的解决方案。

三维成像中的数据冗余

根据三维中心切片定理,仅使用 \phi = 0(即所有投影方向都在横断面内)的投影数据 g(\psi, \phi = 0, \mathbf{s})\psi \in [0, \pi[)就足以重建三维图像 f(x, y, z)。这些数据对应的频率切片是一系列通过 \nu_z 轴的平面,它们的并集完整覆盖了三维频率空间。

然而,采集 \phi \neq 0 的倾斜投影数据可以提高信噪比。这是因为更多的投影数据意味着更多的计数统计量,虽然频率空间覆盖存在冗余,但额外的数据可以用于降低重建噪声。这就是三维PET相比二维PET灵敏度更高的原因。

image-20260114160015067

需要注意的是,在三维成像中,由于探测器的有限轴向范围,可能无法测量完整的 4\pi 立体角投影,导致对象的截断问题。后续讨论假设投影数据是完整的。

n维反投影算子

反投影的定义

n 维反投影算子将定义在投影线流形 \mathbb{T}^{(2n-2)} 上的函数 g 变换为定义在 \mathbb{R}^n 上的函数 \mathbf{X}^* g,定义为:

\mathbf{X}^* g(\mathbf{r}) = \int_{\mathbb{S}^{n-1}} g(\mathbf{e}_\theta, \mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta) \mathrm{d}\mathbf{e}_\theta

这个定义的几何意义是:对于图像空间中的点 \mathbf{r},遍历所有可能的投影方向 \mathbf{e}_\theta,对于每个方向,找到通过点 \mathbf{r} 的投影线,该线的投影平面坐标为 \mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta(即 \mathbf{r} 在投影平面上的投影),读取投影数据 g 在该位置的值,对所有方向的贡献进行积分。

反投影是X射线变换的伴随算子

反投影算子 \mathbf{X}^* 是X射线变换 \mathbf{X} 的伴随算子(adjoint),满足:

\int_{\mathbb{S}^{n-1}} \int_{\theta^\perp} g(\mathbf{e}_\theta, \mathbf{s}) \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) \mathrm{d}\mathbf{e}_\theta \mathrm{d}\mathbf{s} = \int_{\mathbb{R}^n} \mathbf{X}^* g(\mathbf{r}) f(\mathbf{r}) \mathrm{d}\mathbf{r}

这个伴随关系在泛函分析中具有重要意义,它说明反投影是投影操作在某种内积意义下的"转置"。

反投影不是逆变换

反投影算子不是X射线变换的逆算子:

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) \neq f(\mathbf{r})

直接对投影数据做反投影会得到模糊的图像,而非原始图像。这是因为反投影只是"涂抹"操作,将每条投影线上的值均匀地分配回线上所有点,这会导致低频成分被过度放大。正确的重建需要在反投影之前或之后进行滤波校正。

二维与三维反投影的具体形式

反投影的一般定义为:

\mathbf{X}^* g(\mathbf{r}) = \int_{\mathbb{S}^{n-1}} g(\mathbf{e}_\theta, \mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta) \mathrm{d}\mathbf{e}_\theta
二维反投影

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

\mathbf{X}^* g(x, y) = \int_{\mathbb{S}^1} g(\psi, x\cos\psi + y\sin\psi) \mathrm{d}\psi

对于点 (x, y),表达式 s = x\cos\psi + y\sin\psi 给出了该点在角度 \psi 投影中的位置。反投影将所有通过点 (x, y) 的投影线上的投影值累加起来。

三维反投影

在三维情况下,图像 f(x, y, z) 定义在 \mathbb{R}^3 上,投影数据 g(\psi, \phi, s_1, s_2) 定义在 \mathbb{T}^4 上。反投影为:

\mathbf{X}^* g(x, y, z) = \int_{\mathbb{S}^2} g(\psi, \phi, x\cos\psi + y\sin\psi, x\sin\psi\sin\phi - y\cos\psi\sin\phi + z\cos\phi) \mathrm{d}\psi \mathrm{d}\phi

积分遍历球面 \mathbb{S}^2 上的所有方向,对于每个方向读取投影数据在相应位置的值并累加。

正向投影与反投影的几何关系

从图示可以直观理解正向投影和反投影的几何操作。

正向投影

正向投影 \mathbf{X} 将二维图像 f(x, y) 变换为投影数据。对于每个角度 \psi(如 45°90°135°180°),沿该方向对图像做线积分,得到一维投影曲线 \mathbf{X}f(\psi, s)。所有角度的投影组合成正弦图。

反投影

反投影 \mathbf{X}^* 将投影数据变换回图像空间。对于每个角度 \psi,将该角度的一维投影 g(\psi, s) 沿投影方向"涂抹"回二维平面,即将投影值分配给该投影线上的所有点。所有角度的涂抹结果叠加,形成反投影图像 \mathbf{X}^* g(x, y)

image-20260114160100434

从图示可以看到,反投影结果呈现出从中心向外辐射的条纹状伪影,图像整体模糊。这正是因为反投影不是逆变换,需要配合滤波才能得到正确的重建结果。

二维反投影的扩散效应

image-20260114160117976

图示展示了反投影过程中角度逐渐增加时的累积效果。最上排显示了单个角度(0°)的投影及其反投影,呈现为水平条纹。随着更多角度的投影被加入(2个方向、4个方向),条纹开始交叉形成网格状图案。

下排显示了随着投影角度数量继续增加,反投影图像逐渐接近原始图像的形态,但始终呈现出模糊的特征。即使使用了完整 [0, 2\pi[ 范围内的所有角度,反投影结果仍然是原始图像的模糊版本,边缘不清晰,存在从中心向外辐射的"星形"伪影。这种模糊效应正是因为反投影不是X射线变换的逆变换。

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

卷积的定义

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

(f * h)(\mathbf{r}) = \int_{\mathbb{R}^n} f(\mathbf{r} - \mathbf{u}) h(\mathbf{u}) \mathrm{d}\mathbf{u}

卷积在傅里叶域对应于乘积:

\mathcal{F}_n\{(f * h)\}(\boldsymbol{\nu}) = \mathcal{F}_n\{f\}(\boldsymbol{\nu}) \mathcal{F}_n\{h\}(\boldsymbol{\nu})
核心定理:反投影-投影的卷积表示

对于 \mathbb{R}^n 中的函数 f(\mathbf{r}),有:

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) = 2 \|\mathbf{r}\|^{1-n} * f(\mathbf{r})

这个定理说明:先投影再反投影等价于原图像与核函数 2\|\mathbf{r}\|^{1-n} 的卷积。这解释了为什么反投影会产生模糊:卷积核 \|\mathbf{r}\|^{1-n} 是一个低通滤波器,它压制了高频成分。

定理证明

从反投影的定义出发:

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) = \int_{\mathbb{S}^{n-1}} \mathbf{X}f(\mathbf{e}_\theta, \mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta) \mathrm{d}\mathbf{e}_\theta

代入X射线变换的定义:

= \int_{\mathbb{S}^{n-1}} \int_{\mathbb{R}^1} f(\mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta + t\mathbf{e}_\theta) \mathrm{d}t \mathrm{d}\mathbf{e}_\theta

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

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) = \int_{\mathbb{S}^{n-1}} \int_{\mathbb{R}^1} f(\mathbf{r} + u\mathbf{e}_\theta) \mathrm{d}u \mathrm{d}\mathbf{e}_\theta = 2 \int_{\mathbb{S}^{n-1}} \int_0^\infty f(\mathbf{r} + u\mathbf{e}_\theta) \mathrm{d}u \mathrm{d}\mathbf{e}_\theta

再进行变量替换 \mathbf{v} = -u\mathbf{e}_\theta

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) = 2 \int_{\mathbb{R}^n} f(\mathbf{r} - \mathbf{v}) \|\mathbf{v}\|^{1-n} \mathrm{d}\mathbf{v}

这正是卷积的形式,证毕。

二维反投影的滤波校正

对于二维函数 f: \mathbb{R}^2 \rightarrow \mathbb{R}(即 n = 2),卷积关系简化为:

\mathbf{X}^* \mathbf{X} f(\mathbf{r}) = 2 \frac{1}{\|\mathbf{r}\|} * f(\mathbf{r})

图示展示了原始图像 f(x,y) 经过 \mathbf{X}^* \mathbf{X} 操作后变得模糊的效果。

image-20260114160251991
频域分析

利用傅里叶变换的性质:

\mathcal{F}_2\left\{\frac{1}{\|\mathbf{r}\|}\right\}(\boldsymbol{\nu}) = \frac{1}{\|\boldsymbol{\nu}\|}

1/\|\mathbf{r}\| 的二维傅里叶变换是 1/\|\boldsymbol{\nu}\|。结合卷积定理,在频域中有:

\mathcal{F}_2\{f\}(\boldsymbol{\nu}) = \frac{1}{2} \mathcal{F}_2\{\mathbf{X}^* \mathbf{X} f\}(\boldsymbol{\nu}) \|\boldsymbol{\nu}\|

这表明:要从反投影结果恢复原图像,需要在频域乘以 \|\boldsymbol{\nu}\|,即应用斜坡滤波器(ramp filter)。\|\boldsymbol{\nu}\| 滤波器增强高频、压制低频,正好补偿反投影造成的模糊。

二维滤波反投影定理

滤波反投影方法交换了滤波和反投影的顺序,使得计算更加高效。

定理陈述

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

f(x, y) = \frac{1}{2} \mathbf{X}^* g^f(x, y)

其中 g^f(\psi, s) 是滤波后的投影,定义为:

g^f(\psi, s) = \mathcal{F}_1^{-1}\{|\nu_\perp| \mathcal{F}_1\{g\}(\psi, \nu_\perp)\}(\psi, s)

这里 \mathcal{F}_1 表示对投影变量 s 的一维傅里叶变换,|\nu_\perp| 是斜坡滤波器。

利用投影对称性简化

由于投影满足对称性 \mathbf{X}f(\psi, s) = \mathbf{X}f(\psi + \pi, -s),只需测量 \psi \in [0, \pi[ 范围的投影即可完成重建:

f(x, y) = \int_0^\pi g^f(\psi, x\cos\psi + y\sin\psi) \mathrm{d}\psi
定理证明

从傅里叶逆变换出发:

f(x, y) = \mathcal{F}_2^{-1}\{\mathcal{F}_2\{f\}(\nu_x, \nu_y)\}(x, y) = \int_{\mathbb{R}^2} \mathcal{F}_2\{f\}(\nu_x, \nu_y) e^{+2\pi i(x\nu_x + y\nu_y)} \mathrm{d}\nu_x \mathrm{d}\nu_y

将笛卡尔坐标转换为极坐标:\nu_x = \nu_\perp \cos\psi\nu_y = \nu_\perp \sin\psi,雅可比行列式为 \nu_\perp

= \int_0^{2\pi} \int_0^\infty \mathcal{F}_2\{f\}(\nu_\perp \cos\psi, \nu_\perp \sin\psi) e^{+2\pi i \nu_\perp (x\cos\psi + y\sin\psi)} \nu_\perp \mathrm{d}\nu_\perp \mathrm{d}\psi

应用二维中心切片定理,将图像傅里叶变换替换为投影傅里叶变换:

= \int_0^{2\pi} \int_0^\infty \mathcal{F}_1\{\mathbf{X}f\}(\psi, \nu_\perp) e^{+2\pi i \nu_\perp (x\cos\psi + y\sin\psi)} \nu_\perp \mathrm{d}\nu_\perp \mathrm{d}\psi

这个表达式可以分解为:先对每个角度的投影应用斜坡滤波器 |\nu_\perp| 并做逆傅里叶变换(得到滤波投影 g^f),然后对所有角度做反投影。

二维滤波反投影定理的证明

前面已经从傅里叶逆变换出发,将笛卡尔坐标转换为极坐标,并应用中心切片定理将图像傅里叶变换替换为投影傅里叶变换。现在继续完成证明的最后步骤。

从极坐标积分到滤波反投影

证明的关键在于将双重积分(对频率 \nu_\perp 和角度 \psi 的积分)重新组织为先滤波后反投影的形式。从上一步得到的表达式出发:

f(x,y) = \int_0^{2\pi} \int_0^\infty \mathcal{F}_1\{\mathbf{X}f\}(\psi, \nu_\perp) e^{+2\pi i \nu_\perp (x\cos\psi + y\sin\psi)} \nu_\perp \mathrm{d}\nu_\perp \mathrm{d}\psi

这里积分限是 \nu_\perp \in [0, \infty)\psi \in [0, 2\pi)。利用投影的对称性 \mathbf{X}f(\psi, s) = \mathbf{X}f(\psi + \pi, -s),可以将角度范围减半为 \psi \in [0, \pi),同时将频率范围扩展为 \nu_\perp \in (-\infty, \infty)。这是因为 \psi 增加 \pi 等价于 \nu_\perp 取负值:

f(x,y) = \int_0^{\pi} \int_{-\infty}^{\infty} \mathcal{F}_1\{\mathbf{X}f\}(\psi, \nu_\perp) |\nu_\perp| e^{+2\pi i \nu_\perp (x\cos\psi + y\sin\psi)} \mathrm{d}\nu_\perp \mathrm{d}\psi

这里 \nu_\perp 前面的因子从 \nu_\perp 变成了 |\nu_\perp|,这正是斜坡滤波器的频域表达式。最内层对 \nu_\perp 的积分可以识别为一维逆傅里叶变换,因此:

f(x,y) = \int_0^{\pi} \mathcal{F}_1^{-1}\{|\nu_\perp| \mathcal{F}_1\{\mathbf{X}f\}\}(\psi, x\cos\psi + y\sin\psi) \mathrm{d}\psi

这个表达式的结构清晰地展示了滤波反投影的两个步骤:首先对每个角度 \psi 的投影进行滤波操作(频域乘以 |\nu_\perp| 后做逆傅里叶变换),然后对滤波后的投影在 s = x\cos\psi + y\sin\psi 处取值并对所有角度积分(即反投影操作)。

滤波反投影算法的实现步骤

滤波反投影(Filtered Back-Projection,FBP)算法将上述理论公式转化为可执行的计算流程。设采集到的投影数据为 g(\psi, s) = \mathbf{X}f(\psi, s),重建图像 f(x,y) 的步骤如下。

第一步:投影的傅里叶变换

对每个角度 \psi 处的一维投影 g(\psi, s) 做一维傅里叶变换,得到频域表示:

\hat{g}(\psi, \nu_\perp) = \mathcal{F}_1\{g\}(\psi, \nu_\perp)

这一步将投影从空间域转换到频率域,为后续的滤波操作做准备。

第二步:斜坡滤波

在频域对投影乘以斜坡滤波器 |\nu_\perp|

\hat{g}^f(\psi, \nu_\perp) = |\nu_\perp| \hat{g}(\psi, \nu_\perp)

斜坡滤波器的作用是补偿反投影造成的低频放大效应。从物理上理解,反投影相当于将投影值沿射线方向涂抹回图像空间,这个涂抹过程会使低频成分累积过多,而斜坡滤波器通过增强高频、压制低频来抵消这一效应。

第三步:逆傅里叶变换

对滤波后的频域数据做一维逆傅里叶变换,得到空间域的滤波投影:

g^f(\psi, s) = \mathcal{F}_1^{-1}\{\hat{g}^f\}(\psi, s)

滤波投影 g^f(\psi, s) 与原始投影 g(\psi, s) 的主要区别在于边缘处出现了负值(下冲),这些负值在反投影时会抵消相邻区域的贡献,从而消除模糊。

第四步:反投影

将角度 \psi 处的滤波投影反投影到图像空间。对于图像点 (x, y),计算其在该角度投影中的位置 s = x\cos\psi + y\sin\psi,读取滤波投影在该位置的值,并累加到图像中:

f(x,y) = f(x,y) + g^f(\psi, x\cos\psi + y\sin\psi) \Delta\psi

这里 \Delta\psi 是角度采样间隔,对应于积分的离散化。

第五步:角度循环

对所有采集角度 \psi \in [0, \pi) 重复执行第一步到第四步,完成后的累积结果即为重建图像 f(x,y)

斜坡滤波的效果

均匀圆盘的投影与滤波

以均匀圆盘为例可以直观理解滤波的作用。均匀圆盘在任意角度的投影 g(\psi, s) 具有相同的形状:一条平滑的弧形曲线,中心最高、两端逐渐降低到零。这是因为穿过圆盘中心的射线路径最长,线积分值最大,靠近边缘的射线路径较短,线积分值较小。

经过斜坡滤波后,滤波投影 g^f(\psi, s) 的形状发生显著变化:中心区域仍为正值但幅度降低,而边缘处出现明显的负值下冲。这些负值的物理意义是:在反投影时,边缘处的负贡献会抵消相邻区域的正贡献,从而在最终图像中产生清晰的边界,而非模糊的渐变过渡。

滤波对反投影过程的影响

从反投影的逐步累积过程可以清楚地看到滤波的作用。

image-20260114160501481

在不带滤波的情况下,随着参与反投影的角度数量增加,图像逐渐从单一方向的条纹演变为交叉网格,最终收敛到一个模糊的结果,边缘不清晰,存在向外辐射的星形伪影。这正是前面讨论的反投影不等于逆变换的直观体现。

image-20260114160532988

在带滤波的情况下,每个角度的滤波投影在反投影时会产生正负相间的条纹。随着角度数量增加,这些条纹相互叠加时,正值在目标区域累积,负值在目标区域外相互抵消。最终结果是清晰重建出原始的图像结构,边缘锐利,背景干净。滤波投影中的负值起到了关键的抵消作用,确保反投影过程能够正确恢复原始图像。

离散采样与频率限制

投影数据的离散采样

实际采集系统无法测量连续的线积分,而是按照离散的采样网格进行测量。投影数据在两个维度上被采样:径向方向的采样间隔为 \Delta s,角度方向的采样间隔为 \Delta \psi。离散采样得到的数据记为 g_{j,k} = g(\psi_j, s_k),其中 \psi_j = j \cdot \Delta\psis_k = k \cdot \Delta s。对于直径为 2R 的视野,径向采样范围覆盖 s \in [-R, R]

image-20260114160624708
奈奎斯特-香农采样定理的约束

根据奈奎斯特-香农采样定理,采样间隔 \Delta s 决定了能够无混叠重建的最大空间频率。具体而言,可重建的最大频率为:

\nu_{\max} = \frac{1}{2\Delta s}

超过这个频率的信息会发生频谱混叠,无法正确恢复。这意味着离散采样从根本上限制了重建图像的空间分辨率。

带限重建与低通滤波

由于采样定理的限制,实际重建得到的不是原始图像 f(x,y),而是经过低通滤波后的版本 W_c * f(x,y),其中 W_c 是截止频率为 \nu_c 的低通滤波器,且 \nu_c \leq \frac{1}{2\Delta s}。低通滤波器 W_c 的二维傅里叶变换具有以下形式:

\mathcal{F}_2\{W_c\}(\boldsymbol{\nu}) = \Phi\left(\frac{|\boldsymbol{\nu}|}{\nu_c}\right)

其中 \Phi 是一个窗函数,满足 0 \leq \Phi(|\boldsymbol{\nu}|/\nu_c) \leq 1,且当 |\boldsymbol{\nu}|/\nu_c \geq 1\Phi(|\boldsymbol{\nu}|/\nu_c) = 0。窗函数 \Phi 的作用是将高于截止频率的分量完全滤除。

离散采样下的滤波反投影定理

定理陈述

对于图像空间中的低通滤波器 W_c \in \mathbb{R}^2,存在投影空间中的对应滤波器 w_c \in \mathbb{T}^2,使得:

W_c * f = \mathbf{X}^*(w_c * \mathbf{X}f)

这里 \mathbf{X}^* w_c = W_c,即投影空间滤波器的反投影等于图像空间滤波器。这个定理建立了图像域低通滤波与投影域滤波之间的对应关系。

投影域滤波器的频域表达式

如果图像域低通滤波器的傅里叶变换为 \mathcal{F}_2\{W_c\}(\boldsymbol{\nu}) = \Phi(|\boldsymbol{\nu}|/\nu_c),则投影域滤波器的一维傅里叶变换为:

\mathcal{F}_1\{w_c\}(\psi, \nu_\perp) = \frac{|\nu_\perp|}{2} \Phi\left(\frac{|\nu_\perp|}{\nu_c}\right)

这个表达式表明投影域滤波器是斜坡滤波器 |\nu_\perp| 与低通窗函数 \Phi 的乘积(除以2是归一化因子)。

离散数据的滤波反投影算法

基于上述定理,可以推导出针对离散投影数据 g = \mathbf{X}f 的滤波反投影算法。重建公式为:

f_{\text{FBP}}(x,y) = \frac{1}{2} \mathbf{X}^* g^f(x,y)

其中滤波投影 g^f 的傅里叶变换定义为:

\mathcal{F}_1\{g^f\}(\psi, \nu_\perp) = |\nu_\perp| \Phi\left(\frac{|\nu_\perp|}{\nu_c}\right) \mathcal{F}_1\{g\}(\psi, \nu_\perp)

滤波器由两部分组成:斜坡滤波器 |\nu_\perp| 补偿反投影的模糊效应,低通窗函数 \Phi(|\nu_\perp|/\nu_c) 限制频率范围以避免混叠。

连续数据与离散数据的重建对比

连续数据的理想情况

当投影数据 g = \mathbf{X}f 是连续可得的(理论情况),滤波反投影可以精确重建原始图像:

f(x,y) = \frac{1}{2} \mathbf{X}^* g^f(x,y)

滤波投影的傅里叶变换为:

\mathcal{F}_1\{g^f\}(\psi, \nu_\perp) = |\nu_\perp| \mathcal{F}_1\{g\}(\psi, \nu_\perp)

此时使用的是纯斜坡滤波器,没有频率截止,理论上可以恢复所有频率成分。

离散采样数据的实际情况

当投影数据以间隔 \Delta s 离散采样时,记采样数据为 g_{\Delta s} = \mathbf{X}f,重建结果为:

f_{\text{FBP}}(x,y) = \frac{1}{2} \mathbf{X}^* g^f_{\Delta s}(x,y) = W_c * f(x,y) \neq f(x,y)

滤波投影的傅里叶变换变为:

\mathcal{F}_1\{g^f_{\Delta s}\}(\psi, \nu_\perp) = |\nu_\perp| \Phi\left(\frac{|\nu_\perp|}{\nu_c}\right) \mathcal{F}_1\{g_{\Delta s}\}(\psi, \nu_\perp)

重建结果是原始图像与低通滤波器的卷积,高频细节不可避免地丢失。采样间隔 \Delta s 越大,截止频率 \nu_c 越低,图像越模糊。

矩形低通滤波器

滤波器的频域定义

最简单的低通滤波器是矩形滤波器(也称为理想低通滤波器),其窗函数 \Phi 定义为:

\Phi\left(\frac{|\nu|}{\nu_c}\right) = \begin{cases} 1 & \frac{|\nu|}{\nu_c} \leq 1 \\ 0 & \frac{|\nu|}{\nu_c} > 1 \end{cases}

矩形滤波器在截止频率 \nu_c 以内完全通过,超过截止频率则完全阻断,没有过渡带。

image-20260114160750669
投影域滤波器的形式

结合斜坡滤波器,投影域的总滤波器在频域为:

\mathcal{F}_1\{w_c\}(\psi, \nu_\perp) = \frac{1}{2} |\nu_\perp| \Phi\left(\frac{|\nu_\perp|}{\nu_c}\right)

这个滤波器在 |\nu_\perp| < \nu_c 范围内呈线性增长(斜坡特性),在 |\nu_\perp| = \nu_c 处截断为零。从频域图形上看,它呈现为从原点出发的两条斜线段,在 \pm\nu_c 处突然降为零。

空间域滤波核

矩形滤波器的空间域表达式 w_c(\psi, s) 可以通过逆傅里叶变换得到:

w_c(\psi, s) = \begin{cases} \displaystyle\frac{\nu_c^2}{4\pi^2} \left( \frac{\cos s\nu_c - 1}{(s\nu_c)^2} + \frac{\sin s\nu_c}{s\nu_c} \right) & s \neq 0 \\ \displaystyle\frac{\nu_c^2}{8\pi^2} & s = 0 \end{cases}

空间域滤波核呈现振荡衰减的特性:中心处为正峰值,两侧出现交替的正负振荡,振幅随距离增大而衰减。这种振荡特性是矩形滤波器频域截断的直接后果,会在重建图像的边缘处产生振铃伪影(Gibbs现象)。滤波反投影的完整公式可以写为:

f_{\text{FBP}} = \mathbf{X}^*(w_c * g)

即先对投影数据做卷积滤波,再进行反投影。

病态逆问题

平行射线变换求逆的数学特性

平行射线变换的求逆是一个病态问题(ill-posed problem)。所谓病态,是指解对数据的微小变化不具有连续依赖性。滤波反投影解 f_{\text{FBP}} = \frac{1}{2}\mathbf{X}^* g^f 不连续地依赖于数据 g,这意味着数据中的微小扰动可能导致重建图像中任意大的误差。

发射断层成像中的典型扰动

在发射断层成像中,数据扰动的主要来源是泊松随机噪声。由于每个投影单元的期望计数通常很低(个位数量级),泊松噪声的相对幅度很大。病态性质使得这种噪声在重建过程中被严重放大,导致重建图像质量对计数统计量极为敏感。

噪声与计数统计量的关系

模拟实验可以直观展示噪声对FBP重建质量的影响。以双椭圆幻影为例,原始图像 f(\mathbf{x}) 由两个相邻的椭圆组成。在不同总计数水平下进行泊松噪声模拟并用FBP重建,结果显示:当总计数为 10^2 时,重建图像完全被噪声淹没,无法辨识任何结构,10^3 计数时仍以噪声为主,10^4 计数时开始隐约可见结构轮廓,10^5 计数时噪声明显降低但仍可见,10^6 计数时图像质量较好,10^7 计数时接近理想重建。这一系列图像说明FBP重建的图像质量与计数统计量直接相关,低计数情况下噪声问题极为严重。

image-20260114160925637

斜坡滤波器的噪声放大效应

功率谱分析

理解噪声放大的关键在于分析斜坡滤波器在频域的作用。考虑投影数据的功率谱,信号分量的功率谱通常随频率增加而衰减(因为图像的主要能量集中在低频),而泊松噪声的功率谱近似为白噪声(各频率分量能量相当)。在低频区域,信号功率远大于噪声功率,在高频区域,信号功率衰减而噪声功率保持不变,两者趋于相当甚至噪声占优。

当投影数据通过斜坡滤波器 |\nu_\perp| 时,滤波器对高频的增益大于对低频的增益。这导致滤波后的功率谱发生显著变化:信号分量的高频部分被适度增强,但噪声分量的高频部分被大幅放大。在滤波后的功率谱中,高频区域噪声功率可能远超信号功率,导致重建图像中出现大量高频噪声纹理。

159,000计数的模拟示例

以159,000总计数的模拟为例,原始正弦图在添加泊松噪声后呈现明显的颗粒感。经过斜坡滤波后,正弦图中的高频噪声被显著放大,条纹状噪声结构更加明显。最终的FBP重建图像虽然保留了基本结构,但叠加了大量高频噪声,图像质量明显下降。

image-20260114160955128

高频截止策略

截止频率的选择

抑制高频噪声放大的直接方法是不重建超过某个截止频率 \nu_c 的频率分量。截止频率应满足:

\nu_c < \nu_{\text{Nyquist}} = \frac{1}{2\Delta s}

即截止频率低于奈奎斯特频率。通过选择适当的 \nu_c,可以在信号功率仍显著高于噪声功率的频率范围内截断,避免噪声主导的高频区域对重建的贡献。

image-20260114161039020
截止频率对图像质量的影响

截止频率的选择涉及分辨率与噪声之间的权衡。设截止频率为奈奎斯特频率的某个比例 \nu_c = k \cdot \nu_Nk \in (0,1]),不同 k 值对应不同的重建效果:\nu_c = 0.2\nu_N 时,图像非常平滑但严重模糊,细节丢失,\nu_c = 0.4\nu_N 时,模糊程度有所改善但仍明显,\nu_c = 0.6\nu_N 时,开始能够分辨主要结构,\nu_c = 0.8\nu_N 时,结构更清晰但噪声开始可见,\nu_c = 1.0\nu_N 时,分辨率最高但噪声也最明显。选择最优截止频率需要根据具体应用场景在分辨率和噪声之间取得平衡。

image-20260114161047936

切趾窗函数

矩形窗的问题

矩形窗(即理想低通滤波器)在截止频率处有陡峭的跳变,频域的这种不连续性在空间域会产生振荡(Gibbs现象),表现为图像边缘处的振铃伪影。此外,矩形窗的频率截断是硬截断,在接近截止频率处信号和噪声同时被保留,噪声抑制效果有限。

切趾窗的作用

切趾窗(apodization window)通过在截止频率附近引入平滑过渡来改善这些问题。切趾窗函数在低频区域接近1(完全通过),在接近截止频率时逐渐衰减到0,没有硬边界。这种平滑过渡有两个好处:一是消除频域不连续性,减少空间域振铃伪影,二是在信噪比较低的高频区域给予较小的权重,进一步抑制噪声。

image-20260114161201082
矩形窗与Hann窗的对比

以矩形窗和Hann窗为例进行对比。信号功率谱随频率衰减,噪声功率谱近似平坦。使用矩形窗时,窗函数在 \nu_c 处突然截断,滤波后的噪声功率谱在高频处仍有显著幅度,重建图像边缘锐利但伴有振铃。使用Hann窗时,窗函数从 \nu_c 前就开始平滑衰减,滤波后的噪声功率谱在高频处被有效压制,重建图像边缘略微柔和但更加干净自然。

常用低通滤波器

广义Hamming窗

广义Hamming窗是一类参数化的窗函数,定义为:

\Phi_H(\nu) = \begin{cases} \alpha + (1-\alpha)\cos\left(\frac{\pi|\nu|}{\nu_c}\right) & |\nu| < \nu_c \\ 0 & |\nu| \geq \nu_c \end{cases}

参数 \alpha 控制窗函数的形状。当 \alpha = 0.5 时得到Hann窗(也称为升余弦窗),窗函数在边界处精确到达零且导数为零,过渡最为平滑。当 \alpha = 0.54 时得到Hamming窗,边界处有小的非零值,主瓣稍窄但旁瓣稍高。\alpha 越接近1,窗函数越接近矩形窗。

image-20260114161232635
Shepp-Logan滤波器

Shepp-Logan滤波器是专门为断层重建设计的窗函数,定义为:

\Phi_S(\nu) = \begin{cases} \displaystyle\frac{\nu_c}{|\nu|\pi}\sin\left(\frac{|\nu|\pi}{\nu_c}\right) & |\nu| < \nu_c \\ 0 & |\nu| \geq \nu_c \end{cases}

这个滤波器的特点是在低频处接近1,随频率增加缓慢下降,在截止频率处到达零。Shepp-Logan滤波器与斜坡滤波器的乘积在空间域具有较好的局部化特性,能有效减少条纹伪影。

n阶Butterworth滤波器

Butterworth滤波器是一类最大平坦幅度响应的滤波器,n阶Butterworth滤波器定义为:

\Phi_B(\nu) = \begin{cases} \displaystyle\frac{1}{\sqrt{1 + \left(\frac{|\nu|}{\nu_c}\right)^{2n}}} & |\nu| < \nu_c \\ 0 & |\nu| \geq \nu_c \end{cases}

阶数 n 控制过渡带的陡峭程度。n = 1 时过渡非常平缓,低通效果明显但会损失较多高频信息,n = 4 时过渡较陡但仍平滑,n = 10 时接近矩形窗但仍无硬截断。Butterworth滤波器的优点是通带内响应非常平坦,不会引入幅度失真。

各滤波器的频率响应比较
image-20260114161215462

从频率响应曲线可以看到各滤波器的特性差异。Shepp-Logan滤波器在整个频率范围内相对较高,对高频衰减较少,保留分辨率但噪声抑制较弱。Hann窗(\alpha = 0.5)和Hamming窗(\alpha = 0.54)特性相近,在高频有显著衰减。低阶Butterworth滤波器(1阶)衰减最为平缓,高阶Butterworth滤波器(10阶)则接近矩形特性。实际应用中,滤波器的选择需要根据噪声水平和分辨率要求进行权衡。

广义Hamming窗的参数效应

广义Hamming窗的定义为:

\Phi_H(\nu) = \begin{cases} \alpha + (1-\alpha)\cos\left(\frac{\pi|\nu|}{\nu_c}\right) & |\nu| < \nu_c \\ 0 & |\nu| \geq \nu_c \end{cases}

参数 \alpha 的取值范围为 [0.5, 1.0],控制窗函数从Hann窗(\alpha = 0.5)到矩形窗(\alpha = 1.0)之间的过渡。

切趾参数与截止频率的联合效应

以含噪正弦图的FBP重建为例,展示不同参数组合的效果。

image-20260114161548108

\alpha = 1(矩形窗)且 \nu_c = \nu_N 时,图像噪声最大,但边缘最锐利。将 \alpha 降低到0.54(Hamming窗)同时保持 \nu_c = \nu_N,噪声明显减少,边缘略微柔化。进一步降低截止频率到 \nu_c = 0.5\nu_N,噪声继续减少但图像开始模糊。当 \nu_c = 0.2\nu_N 时,噪声几乎消失但图像严重模糊,细节完全丢失。这一系列对比清楚地展示了切趾参数和截止频率对重建质量的影响。

截止频率与空间分辨率

分辨率的定量度量

空间分辨率通常用点扩散函数(Point Spread Function,PSF)的半高全宽(Full Width at Half Maximum,FWHM)来度量。FWHM越小表示分辨率越高,图像越清晰。降低截止频率会增大PSF的FWHM,导致空间分辨率下降。

不同滤波器的分辨率特性
image-20260114161601351

以截止频率为横轴(单位为 2\nu_N,即奈奎斯特频率的两倍)、FWHM为纵轴,可以绘制不同滤波器的分辨率曲线。三种滤波器(矩形窗、Shepp-Logan、Hann窗)的FWHM都随截止频率降低而增大,但增大的速率不同。矩形窗在相同截止频率下FWHM最小,分辨率最高,Hann窗FWHM最大,分辨率最低,Shepp-Logan居中。

具体数值示例:当FWHM = 5.0时,使用矩形窗重建的点源图像紧凑清晰,使用Hann窗重建的点源更加弥散。当FWHM增大到5.8、6.2、6.9时,两种滤波器重建的点源都逐渐变大变模糊,但Hann窗始终比矩形窗更平滑。

信号与噪声的权衡

权衡曲线

信号(分辨率)与噪声之间存在内在的权衡关系,可以用二维曲线来表示。以信号质量(与分辨率正相关)为横轴、噪声水平为纵轴,不同滤波器对应不同的权衡曲线。理想情况是右下角(高信号、低噪声),但实际上无法同时达到两者的最优。

image-20260114161648847

矩形窗的权衡曲线位于最右侧,在相同噪声水平下能达到最高的信号质量,但在相同信号质量下噪声也最高。Hann窗的曲线位于左侧,噪声抑制效果最好但信号质量受损最大。Shepp-Logan介于两者之间。

均值图像与方差图像

通过多次噪声实现的蒙特卡洛模拟,可以分别计算重建图像的均值(反映信号)和方差(反映噪声)。使用矩形窗时,均值图像保留了较好的细节,但方差图像显示噪声分布不均匀且整体较高。使用Shepp-Logan滤波器时,均值图像细节略有损失,方差降低。使用Hann窗时,均值图像最为平滑,方差最低但细节丢失明显。

image-20260114161657624

这里需要注意,最优的权衡点取决于图像本身的频率内容。如果图像主要包含低频成分(如大的均匀区域),可以使用较强的平滑而不损失太多信息,如果图像包含精细结构,则需要保留较高的截止频率。

滤波器选择的实践考量

经验选择方法

在实际应用中,滤波器的选择通常是经验性的,基于对图像质量的主观评估。一般原则是:适度平滑后图像质量趋于稳定,过度平滑会导致空间分辨率和对比度恢复的严重退化。

image-20260114161722744

以三组参数设置为例:矩形窗(\nu_c = \nu_N)噪声明显但细节保留,Hann窗(\nu_c = 0.8\nu_N)噪声减少且细节基本保留,是常用的折中选择,Hann窗(\nu_c = 0.5\nu_N)噪声很低但图像明显模糊。

定量评价指标

除主观评估外,还可以使用定量指标来指导滤波器选择。常用的定量指标包括:信噪比(Signal-to-Noise Ratio,SNR),衡量信号强度与噪声水平的比值,均方根误差(Root Mean Square Error,RMSE),衡量重建图像与参考图像之间的偏差,ROC曲线分析,通过数值观察者或人类观察者研究来评估特定检测任务的性能。这些定量指标能够为滤波器参数的优化提供客观依据。

三维滤波反投影

与二维FBP的相似性

三维FBP的基本原理与二维情况类似,都基于中心切片定理。根据三维中心切片定理,方向为 \mathbf{e}_\theta 的二维投影 \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) 的二维傅里叶变换 \mathcal{F}_2\{\mathbf{X}f\}(\psi, \phi = 0, \nu_\perp) 等于三维图像傅里叶变换 \mathcal{F}_3\{f\}(\nu_x, \nu_y, \nu_z) 在通过原点、垂直于 \mathbf{e}_\theta 的平面上的切片。

image-20260114161733272

当投影方向限制在横断面内(\phi = 0)时,所有切片都是包含 \nu_z 轴的垂直平面。这些平面的并集覆盖整个三维频率空间,因此仅使用 \phi = 0 的完整投影 g(\psi, \phi = 0, \mathbf{s})\psi \in [0, \pi))就足以重建三维图像 f(x, y, z)

与二维FBP的差异

三维情况与二维情况有两个主要差异。第一是重建的充分条件不同:三维重建对投影数据的完整性有特定要求,需要满足一定的几何条件才能保证精确重建。第二是FBP滤波器不再唯一:在三维情况下,存在多种不同的滤波器选择都能实现精确重建,这为算法设计提供了更多灵活性。

有限孔径导致的截断问题

实际系统的几何限制

在实际的PET或SPECT系统中,无法采集覆盖完整 4\pi 立体角的投影数据。探测器的有限轴向范围导致只能测量部分方向的投影,这种限制称为有限孔径(limited aperture)。

系统孔径可以用方向集合 \Omega 来描述:

\Omega = \{\mathbf{e}_\theta | g(\mathbf{e}_\theta, \mathbf{s}) \text{ 被测量}\}

对于PET系统,探测器通常呈截断圆柱形,孔径可以表示为:

\Omega_{\text{PET}} = \{\mathbf{e}_\theta | |\phi| < \phi_{\max}\}

其中 \phi_{\max} 是最大接收角,由探测器环的轴向长度与直径之比决定。

平移不变性假设

为简化分析,通常假设孔径具有平移不变性,即对于孔径内的任意方向 \mathbf{e}_\theta \in \Omega,该方向的投影数据 g(\mathbf{e}_\theta, \mathbf{s}) 对所有投影平面位置 \mathbf{s} \in \theta^\perp 都被测量,或者缺失数据被设为零。这个假设在实际中并不严格成立(因为视野边缘处存在截断),但它简化了数学处理,使得可以使用类似二维FBP的方法进行三维重建。

三维重建的Orlov条件

频率空间覆盖的几何条件

根据三维中心切片定理,方向为 \mathbf{e}_\theta 的投影的二维傅里叶变换对应于三维频率空间中通过原点、垂直于 \mathbf{e}_\theta 的平面。要完整重建三维图像,必须能够恢复三维频率空间中的每一个频率值 \mathcal{F}_3\{f\}(\boldsymbol{\nu})

image-20260114161905099

Orlov条件给出了精确重建的几何判据:对于三维频率空间中的任意频率向量 \boldsymbol{\nu},在单位球面 \mathbb{S}^2 上以 \boldsymbol{\nu} 方向为法向量的赤道圆必须与孔径区域 \Omega^2 相交。换言之,孔径 \Omega 必须足够大,使得对于任意频率方向,至少存在一个测量方向能够提供该频率的信息。

不同探测器几何的孔径特性

以三种典型的探测器几何为例说明Orlov条件的应用。

image-20260114161929272

二维PET的孔径 \Omega_0 是单位球面赤道上的一个环带,仅包含 \phi = 0(横断面内)的方向。这个环带与所有包含 z 轴的赤道圆都相交,因此满足Orlov条件,可以进行精确重建。

三维圆柱形PET的孔径 \Omega_{\phi_{\max}} 是以赤道为中心、张角为 \phi_{\max} 的球带。当 \phi_{\max} 足够大时,该球带与单位球面上的所有赤道圆都相交,满足Orlov条件。

双平面探测器的孔径 \Omega_{\text{planar}} 由两个分离的小区域组成,对应于两个固定探测器平面的法向量附近的方向。这种孔径通常无法与所有赤道圆相交,不满足Orlov条件,只能进行有限角度重建,存在缺失频率。

三维滤波反投影算法

算法的一般形式

对于三维平行射线投影的最一般线性平移不变逆变换是三维滤波反投影算法。设投影数据 g(\mathbf{e}_\theta, \mathbf{s}) = \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) 在平移不变孔径 \Omega 上被测量,则重建公式为:

f(\mathbf{r}) = \mathbf{X}^*_\Omega g^f(\mathbf{r})

其中 \mathbf{X}^*_\Omega 表示在孔径 \Omega 范围内的反投影算子,g^f 是滤波后的投影。

滤波投影的定义

滤波投影 g^f(\mathbf{e}_\theta, \mathbf{s}) 通过二维傅里叶变换在投影平面上定义:

g^f(\mathbf{e}_\theta, \mathbf{s}) = \mathcal{F}_2^{-1}\{\mathcal{F}_2\{g\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp)\}(\mathbf{e}_\theta, \mathbf{s})

其中 \boldsymbol{\nu}_\perp 是投影平面内的二维频率向量,满足 \mathbf{e}_\theta \cdot \boldsymbol{\nu}_\perp = 0。滤波器 H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) 依赖于孔径 \Omega,由于数据的冗余性,滤波器的选择不是唯一的。

三维FBP滤波器的非唯一性

滤波器有效性条件

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

\iint_\Omega \delta(\boldsymbol{\nu} \cdot \mathbf{e}_\theta) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) \mathrm{d}^2\mathbf{e}_\theta = 1

对于三维频率空间中的任意 \boldsymbol{\nu} \in \mathbb{R}^3 都成立,其中 \boldsymbol{\nu}_\perp\boldsymbol{\nu} 在投影平面 \theta^\perp 上的投影。这个条件确保每个频率分量在重建过程中被正确归一化。

滤波器选择对噪声传播的影响

满足有效性条件的所有滤波器对于一致性(无噪声)数据产生相同的重建结果。然而,不同滤波器在噪声传播方面的表现不同:某些滤波器会放大特定方向或频率的噪声,而另一些滤波器可能具有更好的噪声抑制特性。因此,滤波器的选择会影响重建图像的噪声特性和方差分布。

滤波器有效性条件的证明

证明思路

对三维FBP公式两边取三维傅里叶变换:

\mathcal{F}_3\{f\}(\boldsymbol{\nu}) = \iint_\Omega \mathrm{d}\mathbf{e}_\theta^2 \iiint_{\mathbb{R}^3} \mathrm{d}\mathbf{r} \, g^f(\mathbf{e}_\theta, \mathbf{r} - (\mathbf{r} \cdot \mathbf{e}_\theta)\mathbf{e}_\theta) e^{-2\pi i \mathbf{r} \cdot \boldsymbol{\nu}}

引入变量替换 r_\parallel = \mathbf{r} \cdot \mathbf{e}_\theta\mathbf{r}_\perp = \mathbf{r} - r_\parallel \mathbf{e}_\theta,将积分分解为沿投影方向和垂直于投影方向的两部分:

\mathcal{F}_3\{f\}(\boldsymbol{\nu}) = \iint_\Omega \mathrm{d}\mathbf{e}_\theta^2 \int_{\mathbb{R}} \mathrm{d}r_\parallel e^{-2\pi i r_\parallel \mathbf{e}_\theta \cdot \boldsymbol{\nu}} \mathcal{F}_2\{g\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp)
应用中心切片定理

利用三维中心切片定理 \mathcal{F}_2\{g\}(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp),上式变为:

\mathcal{F}_3\{f\}(\boldsymbol{\nu}) = \iint_\Omega \mathrm{d}\mathbf{e}_\theta^2 \int_{\mathbb{R}} \mathrm{d}r_\parallel e^{-2\pi i r_\parallel \mathbf{e}_\theta \cdot \boldsymbol{\nu}} \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp)
利用δ函数的积分表示

注意到 \delta(x) = \lim_{b \to \infty} \int_{-b}^{b} \mathrm{d}y \exp(-2\pi i xy),对 r_\parallel 的积分给出:

\mathcal{F}_3\{f\}(\boldsymbol{\nu}) = \iint_\Omega \mathrm{d}\mathbf{e}_\theta^2 \, \delta(\mathbf{e}_\theta \cdot \boldsymbol{\nu}) \mathcal{F}_3\{f\}(\boldsymbol{\nu}_\perp) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp)

由于这个等式必须对任意 f(\mathbf{r}) 成立,滤波器必须满足:

\iint_\Omega \mathrm{d}^2\mathbf{e}_\theta \, \delta(\boldsymbol{\nu} \cdot \mathbf{e}_\theta) H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = 1 \quad \forall \boldsymbol{\nu} \in \mathbb{R}^3

这就是滤波器有效性条件。

Colsher滤波器

可分解滤波器族

一类常用的滤波器形式是可分解滤波器,即滤波器可以写成频率相关部分和方向相关部分的乘积:

H_\Omega(\mathbf{e}_\theta, \boldsymbol{\nu}_\perp) = H'_\Omega(\boldsymbol{\nu} = \boldsymbol{\nu}_\perp) w(\mathbf{e}_\theta)

其中 w(\mathbf{e}_\theta) 是方向权重函数。将这种形式代入滤波器有效性条件,可以求解出频率相关部分:

H'_\Omega(\boldsymbol{\nu}) = \frac{1}{\iint_\Omega \mathrm{d}^2\mathbf{e}_\theta \, \delta(\boldsymbol{\nu} \cdot \mathbf{e}_\theta) w(\mathbf{e}_\theta)} = \frac{|\boldsymbol{\nu}|}{\iint_\Omega \mathrm{d}^2\mathbf{e}_\theta \, \delta\left(\frac{\boldsymbol{\nu}}{|\boldsymbol{\nu}|} \cdot \mathbf{e}_\theta\right) w(\mathbf{e}_\theta)}
等权重选择与方差最小化

当选择等权重 w(\mathbf{e}_\theta) = 1 时,滤波器对所有测量方向赋予相同的权重。可以证明,这种选择能够最小化重建图像的方差,是在噪声传播意义下的最优选择。

Colsher滤波器的表达式

等权重条件下得到的滤波器称为Colsher滤波器:

H_{\text{Colsher}}(\boldsymbol{\nu}) = \frac{|\boldsymbol{\nu}|}{\iint_\Omega \mathrm{d}^2\mathbf{e}_\theta \, \delta\left(\frac{\boldsymbol{\nu}}{|\boldsymbol{\nu}|} \cdot \mathbf{e}_\theta\right)}

分母是在孔径 \Omega 内、与频率向量 \boldsymbol{\nu} 垂直的所有投影方向的测度。对于圆柱形PET孔径,这个测度取决于频率向量相对于 z 轴的倾斜角:当频率向量接近水平(\nu_z 小)时,有更多的投影方向贡献,分母较大,滤波器增益较小,当频率向量接近垂直(\nu_z 大)时,贡献的投影方向较少,分母较小,滤波器增益较大。这种自适应增益特性能够补偿不同频率方向数据冗余度的差异。

image-20260114162007541

Colsher滤波器的显式表达式

对于圆柱形PET探测器的孔径 \Omega_{\phi_{\max}},Colsher滤波器有解析的显式表达式。设 \gamma 为频率向量 \boldsymbol{\nu}z 轴的夹角(即 \sin|\gamma| = |\nu_z|/|\boldsymbol{\nu}|),则:

H_{\text{Colsher}}(\boldsymbol{\nu}) = \frac{|\boldsymbol{\nu}|}{\iint_\Omega \mathrm{d}^2\mathbf{e}_\theta \, \delta\left(\frac{\boldsymbol{\nu}}{|\boldsymbol{\nu}|} \cdot \mathbf{e}_\theta\right)} = \begin{cases} \displaystyle\frac{|\boldsymbol{\nu}|}{2\pi} & \text{if } \sin|\gamma| < \sin\phi_{\max} \\ \displaystyle\frac{|\boldsymbol{\nu}|}{4\arcsin(\sin\phi_{\max}/\sin|\gamma|)} & \text{otherwise} \end{cases}

当频率向量与 z 轴的夹角 |\gamma| 较小(即频率接近横断面内)时,分母为 2\pi,此时有完整一圈的投影方向贡献,滤波器简化为 |\boldsymbol{\nu}|/(2\pi),与二维FBP的斜坡滤波器形式一致。当 |\gamma| 较大(频率向量倾斜角超过孔径张角)时,只有部分投影方向能够贡献,分母减小,滤波器增益相应增大以补偿数据冗余度的降低。

Colsher滤波器的特性

Colsher滤波器与二维FBP的斜坡滤波器具有相似性,因此噪声传播特性也相近。该滤波器是三维PET解析重建的标准滤波器。在实际应用中,通常还需要在投影平面的两个方向 \mathbf{e}_{\theta_1}^\perp\mathbf{e}_{\theta_2}^\perp 上施加切趾窗函数,以抑制高频噪声,这与二维情况类似。

数据截断问题

圆柱形探测器的非平移不变性

对于实际的圆柱形PET探测器,孔径 \Omega 并不是严格平移不变的。投影数据 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 \neq 0 的倾斜投影。

image-20260114162023828

响应线 L(\mathbf{e}_\theta, \mathbf{s}) 由方向向量 \mathbf{e}_\theta(由方位角 \psi 和极角 \phi 确定)以及投影平面上的位置向量 \mathbf{s} 定义。对于横断面投影(\phi = 0),所有轴向位置的数据都被完整采集。但对于倾斜投影(\phi \neq 0),由于探测器轴向长度有限,视野边缘处的倾斜响应线无法被完整测量,导致数据在轴向方向截断。

环差与数据覆盖

在PET扫描仪中,响应线可以用两个探测器环的编号差(环差,ring difference)来表征。环差为零对应横断面投影,环差越大对应越倾斜的投影。由于探测器环数量有限,大环差的响应线在轴向边缘处会被截断。数据截断会导致重建伪影,需要通过特殊的算法来处理。

三维重投影算法

非倾斜投影的特殊性

非倾斜投影 g(\psi, \phi = 0, \mathbf{s}) 具有两个特殊性质:它们是平移不变的(不存在轴向截断问题),并且满足Orlov条件(足以进行完整重建)。这意味着仅使用非倾斜投影就可以进行三维重建,倾斜投影提供的是冗余信息。

3DRP算法的基本思想

三维重投影算法(3D Reprojection Algorithm,3DRP)的核心思想是:首先利用非倾斜投影进行初步的二维FBP重建,然后利用这个初步重建来估计缺失的倾斜投影数据。

3DRP算法的步骤

第一步是对非倾斜投影进行二维FBP重建,得到每个横断面的初步图像。第二步是对初步重建图像进行正向投影,计算缺失的倾斜投影数据的估计值。第三步是将估计的投影数据与实际测量的投影数据融合,得到完整的三维投影数据集。第四步是对完整数据集进行三维FBP重建,得到最终的三维图像。

3DRP算法是三维PET解析重建的标准算法,能够有效处理数据截断问题,同时充分利用倾斜投影数据提供的额外统计信息来提高图像质量。

三维到二维重分箱算法

三维采集与二维重建的结合

三维采集模式相比二维模式具有更高的灵敏度,能够提供更好的信噪比。然而,二维重建算法计算速度快且足够成熟。重分箱(rebinning)算法将三维采集的优势与二维重建的效率结合起来。

重分箱的基本流程

重分箱算法分为两步。第一步是将三维模式采集的数据重新组织为一组二维正弦图。三维数据包含不同环差的响应线,重分箱过程将倾斜响应线的数据合理地分配到对应的横断面正弦图中。常用的重分箱方法包括单层重分箱(Single-Slice Rebinning,SSRB)和傅里叶重分箱(Fourier Rebinning,FORE)。第二步是对重分箱后的二维正弦图堆栈应用快速的二维重建算法(如二维FBP或二维迭代重建)。

image-20260114162122915

从示意图可以看到,原始的三维正弦图数据(包含多个环差的信息)经过重分箱后被压缩为单一的二维正弦图堆栈,每个横断面对应一个二维正弦图。这种数据压缩大大减少了重建的计算量,同时保留了三维采集的灵敏度优势。

解析重建的预校正要求

线积分模型的前提条件

3DRP和FORE等解析重建算法都基于线积分模型,即假设测量数据是活度分布沿响应线的线积分。然而,原始的PET采集数据受到多种物理因素的影响,并不直接满足线积分模型。

必要的数据预校正

在进行解析重建之前,原始PET数据必须经过以下校正:随机符合校正,去除由独立湮灭事件偶然形成的虚假符合,散射符合校正,去除由散射光子形成的错误定位事件,衰减校正,补偿光子穿过组织时的吸收损失,探测效率归一化,补偿不同探测器对的灵敏度差异。

这些校正通常在正弦图域进行,校正后的数据近似满足线积分模型,可以用于解析重建。校正的顺序和方法对最终图像质量有显著影响,是PET数据处理流程中的关键环节。

三维PET解析重建方法

3DRP的地位

3DRP算法是三维PET解析重建的金标准方法,能够精确处理完整的三维投影数据,重建质量高。其主要缺点是计算量大,尤其是正向投影步骤需要遍历整个三维图像空间。

替代性三维算法

FAVOR算法是一种基于FBP滤波器的替代方案,其特点是不要求完整的倾斜投影数据,能够处理平移变化的孔径(shift variant),适用于某些特殊的探测器几何。

FORE算法的影响

FORE(傅里叶重分箱)算法通过将三维数据高效地转换为二维数据,使得可以在三维采集模式下使用快速的二维迭代重建算法。这一创新对三维PET的临床应用产生了重大影响,显著降低了重建时间。在许多系统上,3DRP已经被FORE加二维FBP(FORE+2D-FBP)或FORE加二维直接傅里叶逆变换(FORE+2D-DIFT)所取代。

精确重分箱与飞行时间数据

在傅里叶域存在精确的重分箱方程,可以在不损失信息的情况下进行三维到二维的数据转换。对于飞行时间PET数据,也存在近似和精确的傅里叶重分箱方法,能够在保留飞行时间信息优势的同时实现高效的数据处理和重建。

附录:n维Radon变换

Radon变换的定义

n 维Radon变换将函数 f: \mathbb{R}^n \rightarrow \mathbb{R} 变换为定义在 \mathbb{R}^n 中超平面流形上的函数 \mathbf{R}f。超平面 \Pi 由其法向量 \mathbf{e}_\Pi \in \mathbb{S}^{n-1} 和沿法向量方向到原点的有符号距离 l 唯一确定。Radon变换的定义为:

\mathbf{R}f(\mathbf{e}_\Pi, l) = \int_{\mathbf{r} \cdot \mathbf{e}_\Pi = l} f(\mathbf{r}) \mathrm{d}\mathbf{r}

即函数 f 在超平面 \{\mathbf{r} | \mathbf{r} \cdot \mathbf{e}_\Pi = l\} 上的积分。在二维情况下,超平面退化为直线,在三维情况下,超平面是通常意义上的平面。

附录:平行射线变换与Radon变换的关系

两种变换的联系

平行射线变换(X射线变换)\mathbf{X} 和Radon变换 \mathbf{R} 是密切相关的。对于所有与 \mathbf{e}_\Pi \in \mathbb{S}^{n-1} 正交的方向 \mathbf{e}_\theta \in \mathbb{S}^{n-1},有:

\mathbf{R}f(\mathbf{e}_\Pi, l) = \int_{\mathbf{s} \in \theta^\perp, \mathbf{s} \cdot \mathbf{e}_\Pi = l} \mathbf{X}f(\mathbf{e}_\theta, \mathbf{s}) \mathrm{d}\mathbf{s}

即Radon变换可以通过对平行射线变换在满足特定约束的投影平面位置上积分得到。

二维情况的等价性

在二维情况下(n = 2),Radon变换与平行射线变换是等价的,只是参数化方式不同。对于平行射线变换,参数是射线方向 \mathbf{e}_\theta 和投影平面上的位置 s,对于Radon变换,参数是直线的法向量 \mathbf{e}_\Pi 和到原点的距离 l

两种参数化的关系是:\mathbf{e}_\Pi\mathbf{e}_\theta 正交,即 \mathbf{e}_\Pi = \mathbf{e}_{\theta^\perp}。因此在二维情况下:

\mathbf{R}^{2D}f(\mathbf{e}_{\theta^\perp}, l) = \mathbf{X}^{2D}f(\mathbf{e}_\theta, s)

其中 l = s。这种等价性使得二维断层重建的文献中经常交替使用Radon变换和X射线变换的术语。

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

响应线的索引方式

在三维PET中,响应线可以用多种方式参数化。一种实用的参数化方式是基于轴向位置 z 和环差 \Delta。设响应线连接探测器环 A(轴向位置 z_A)和探测器环 B(轴向位置 z_B),则定义:

z = z_A + z_B, \quad \Delta = z_A - z_B

z 表示响应线的平均轴向位置(或响应线中点的轴向坐标的两倍),\Delta 表示两个探测器环的轴向距离差,反映了响应线的倾斜程度。

参数转换公式

(s, \psi, z, \Delta) 参数化到标准的投影参数 g(\psi, \phi, s, s_2) 的转换公式为:

m(s, \psi, z, \Delta) = g\left(\psi, \arctan\left(\frac{\Delta}{2\sqrt{R^2 - s^2}}\right), s, \frac{z}{2\sqrt{1 + \frac{\Delta^2}{4(R^2 - s^2)}}}\right)

其中 R 是探测器环的半径。极角 \phi 由环差 \Delta 和径向位置 s 共同决定,轴向投影坐标 s_2 则是 z 的缩放版本。这种参数化方式便于按环差组织数据,对于重分箱算法的实现非常方便。

附录:单层重分箱算法

SSRB的基本思想

单层重分箱(Single-Slice Rebinning,SSRB)是一种近似的重分箱算法,其核心思想是将所有穿过某一轴向位置的倾斜响应线的数据简单平均,分配给该轴向位置的横断面正弦图。

SSRB公式

SSRB的数学表达式为:

m_{\text{SSRB}}(s, \psi, z) = \frac{1}{2\Delta_{\max}(z)} \int_{-\Delta_{\max}(z)}^{+\Delta_{\max}(z)} \mathrm{d}\Delta \, m(s, \psi, z, \Delta)

其中 \Delta_{\max}(z) 是轴向位置 z 处允许的最大环差,取决于探测器的几何结构和视野范围。这个公式将所有环差的数据加权平均,得到等效的二维正弦图。

从示意图可以看到,多条不同倾斜角度的响应线(对应不同的环差)穿过同一轴向层面,SSRB将这些响应线的计数简单相加并归一化,作为该层面横断面正弦图的估计值。这种方法计算简单快速,但在大环差或视野边缘处会引入模糊。

附录:傅里叶重分箱算法

FORE的理论基础

傅里叶重分箱(Fourier Rebinning,FORE)是在傅里叶域进行的近似重分箱算法,比SSRB更精确。FORE利用了频率-距离关系:倾斜正弦图的二维傅里叶变换与对应横断面正弦图的傅里叶变换之间存在近似的位移关系。

image-20260114162158402
正弦图的傅里叶变换

对正弦图 m(s, \psi, z, \Delta)(s, \psi) 平面上做二维傅里叶变换:

\mathcal{F}_2\{m\}(\nu_s, k, z, \Delta) = \int_0^{2\pi} \int_{-R}^{+R} m(s, \psi, z, \Delta) e^{-2\pi i(k\psi + \nu_s s)} \mathrm{d}s \mathrm{d}\psi

其中 \nu_s 是径向频率,k 是角度频率(整数)。

频率-距离关系

FORE的核心是频率-距离关系:傅里叶变换 \mathcal{F}_2\{m\} 在频率 (\nu_s, k) 处的值主要来自轴向位置 t = -k/(2\pi\nu_s) 处的源。基于这一关系,可以建立近似等式:

\mathcal{F}_2\{m\}(\nu_s, k, z, \Delta) \approx \mathcal{F}_2\{m\}\left(\nu_s, k, z - \frac{k}{2\pi\nu_s} \frac{\Delta}{2R}, 0\right)

这个关系表明,环差为 \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 2R},将傅里叶分量累加到 M(\nu_s, k, z')。完成所有倾斜正弦图的处理后,对 M(\nu_s, k, z) 进行归一化。最后对每个轴向位置 z 做二维逆傅里叶变换,得到重分箱后的二维正弦图 m(s, \psi, z) = \mathcal{F}_2^{-1}\{M(\nu_s, k, z)\}。对这些二维正弦图进行标准的二维重建即可得到三维图像。

FORE相比SSRB能够更准确地保留空间分辨率,尤其是在大环差和高频区域,是三维PET数据处理的重要工具。

image-20260114162213671

评论