基于深度能量模型的低剂量CT重建.docx

基于深度能量模型的低剂量 CT 重建

摘要：降低计算机断层扫描（CT）的剂量对于降低临床应用中的辐射风险至关重要，深度学 习的快速发展和广泛应用为低剂量 CT 成像算法的发展带来了新的方向。与大多数受益于手动 设计的先验函数或有监督学习方案的现有先验驱动算法不同，本文使用基于深度能量模型来 学习正常剂量 CT 的先验知识，然后在迭代重建阶段，将数据一致性作为条件项集成到低剂 量 CT 的迭代生成模型中，通过郎之万动力学迭代更新训练的先验，实现低剂量 CT 重建。实 验比较，证明所提方法的降噪和细节保留能力优良。
关键词：低剂量CT；深度学习；郎之万动力学；深度能量模型

计算机断层扫描（computed tomography，CT）是一种准确、非侵入式的体内异常检测 技术，可应用于生物学、医学和其他领域。CT 检查的广泛使用引起了人们对 X 射线辐射致 癌和遗传损伤潜在风险的担忧[1] ，于是，减少 CT 扫描次数降低 CT 的辐射剂量成为科研人 员关注的焦点。针对如何在尽可能低的扫描剂量条件下获取与常规剂量 CT 质量相近的 CT 图像，即低剂量 CT（low-dose CT，LDCT）可以对许多临床适应症进行成像，以最大限度地减 少与辐射相关的风险，而不会显著影响筛查或诊断的性能。然而，降低辐射剂量会增加数 据噪声并将伪影引入重建图像，如果不解决这些问题，则会对其诊断性能产生不利影响[2]。
为了解决这个问题，许多研究工作致力于设计算法来提高低剂量 CT 的图像质量。传统 的重建算法一般可以分为解析法和迭代法[3] 。其中解析法具备算法简单、计算速度快等优 点，常用的解析算法有滤波反投影算法（filtered back projection，FBP）；但是，解析 法需要较完整的投影数据，且对噪声敏感，容易产生伪影和噪声污染。迭代算法将重建图 像当作未知变量，以此建立目标函数进行迭代求解，比如全变分模型（total variation， TV）[4]和噪声统计模型[5-6]；虽然这些迭代算法确实提高了图像质量，但是算法复杂度大，运 行时占用内存等问题也限制了其应用。
随着深度学习的快速发展，卷积神经网络（convolutional neural networks，CNN） 在图像去噪领域取得了很好进展[7-9] 。最近，几种基于深度学习的算法已被应用于 CT 重建。

这些方法主要分为两种，一种是使用端到端的有监督深度学习技术作为后处理方法[10-11]， 例如Würfl 等[12]采用神经网络来进行CT 图像重建；Cheng 等[13]使用基于深度学习的跨越策 略模拟了迭代重建过程；Han 等[14]提出了一种基于深度残差学习网络的稀疏投影CT 重建方 法，结果表明，仅由条纹伪影组成的残差图像比原始 CT 图像简单得多；Jin 等[15]结合FBP 和 U-net[16]的深度卷积网络，并充分利用了多分辨率 CNN 和残差学习[17]。
另一种是将深度学习驱动的先验知识集成到迭代方案，该类方法可以将神经网络学 习的先验用来处理各种任务，而无需重新训练或者修改网络。具体来说，Baguer 等[18]使用 深度图像先验进行CT 重建，在低维数据时取得优异的结果。同时，基于生成模型的特征 学习，从网络结构和目标函数的角度提出了许多图像重建方法。最新的进展主要有两种方法： 变分 自编码器（ variational auto-encoder，VAE）[19-20]和生成对抗 网络 [21]（ generative adversarial networks，GAN），前者使用对数似然函数作为训练目标，而后者使用对抗性 训练来最小化模型和数据分布之间的f 散度[22]或积分概率度量[23] 。蔡宁等[24]通过构建和分 析普通CNN 与GAN 两者的性能差异，提出了一种创新的扫描方式，可有效获取低噪声的微 米级计算机断层扫描（Micro-CT）数据，并取得了很好的重建效果；Niu 等[25]提出的基于 相似性的无监督深度去噪方法，能够以非局部和非线性的方式抑制图像中噪声，结果证明 该方式从低剂量 CT 图像中恢复原有的特征效果；Wu 等[26]采用 K-稀疏自动编码器用于无监 督的特征学习任务，并在重建过程中最小化重建图像和流形表示之间的距离以及数据保真 度来提高图像重建质量；Zhang 等[27]提出的REDAEP 算法在去噪自编码先验模型的基础上， 采用变量增强技术和 Lp 正则化约束以挖掘图像更高维的先验信息。遗憾的是，尽管生成模 型在处理自然图像方面取得了成功，但在医学成像领域的研究却不是很多，尤其是在低剂 量 CT 方面。
最近，基于深度能量的模型[29]（ energy-based model，EBM）在生成建模方面取得了可 喜的成果。EBM 通过将标量能量（兼容性的度量）与变量相关联来捕获相关性，做出预测 或决策的推理包括设置观察变量的值并找到使能量最小化的剩余变量的值。EBM 在于找到 一个能量函数，该函数将低能量与剩余变量的正确值相关联，将较高能量与不正确的值相 关联。因此，本文提出一种基于深度能量模型的新型 CT 重建方法，对本方法和现有技术之 间的实验比较，证明本方法的降噪和细节保留能力优良。
1 基于深度能量低剂量 CT 成像
1.1 深度能量模型
2019 年，Du 等[29]的工作已经证明，EBM 作为生成模型应用于自然图像恢复表现出了优 秀的性能，具备显著的优势。首先，它提供明确的非标准化的密度信息，具备组合性和灵 活性；此外，与自回归模型和基于流的模型不同的是，EBM 不需要特殊的模型架构。最近， 在 Nijkamp 等[30]和Du 等[29]的工作中可见，EBM 能够成功地进行最大似然训练，并用于后续 图像重建任务中。
基于能量的模型在机器学习中有着悠久的历史，将无监督生成模型与EBM 相结合是一 种新颖的深度学习方法。它由玻尔兹曼分布显式地定义概率，通常对密度函数建模，旨在 学习一个函数 E，该函数将低能量值分配给数据分布中的输入，并将高能量值分配给其他

输入。它们允许使用样本来生成目标中间过程，其中通过马尔科夫链蒙特卡洛采样方法从 x ∼ exp (−E(x))中采样样本 x 。与使用隐式函数生成样本的 GAN 和求对数似然函数下界的 VAE 等方法相比，显式采样与基于深度能量的模型相结合进行生成建模具有简单、稳定和 灵活等优势。
本文将能量函数用一个神经网络 E(x) 表示， x 是输入数据， θ 为神经网络权重，则概 率分布为 pθ (x) = exp (−Eθ (x))/Z(θ) ，其中 Z(θ) = ∫ exp (−Eθ (x))dx 。从这种分布中对样本采 样是很困难的，以前的工作依赖于随机游走或吉布斯采样等方法。这些方法混合时间长， 特别是对于图像等高维复杂数据就难以处理。
1.2 郎之万动力学采样
郎之万动力学[31]可以通过计算概率密度 p(x) 的对数的梯度 ∇x (logp(x))从而生成样本。
给定一个固定的步长 ε > 0 ，一个初始值 0 ∼ π(x) ，其中 π 是一个先验分布，郎之万动力
学方法递归地通过公式（1）计算：

其中 zt ∼ N(0, I) ，当 ε → 0 和 t → ∞ 时， t 的分布等于 p(x) 。在这种情况下， t 在某些条
件下成为来自 p(x) 的精确样本。当 ε > 0 和 t < ∞ 时，需要一个马尔可夫链蒙特卡罗方法来
更新纠正公式（1），但在实践中这经常被忽略。本文假设在步长 ε 足够小且 t 足够大时，真
实和采样数据间误差可忽略不计。
将以郎之万动力学为基础的公式（1）与深度能量模型概率分布相结合，可得：
（2）

公式（2）中定义概率分布为 qθ , 使得 ˜(x)K ∼ qθ 。正如上述所说，当 ε → 0 和 t → ∞ 时，
qθ → pθ 。此过程通过郎之万动力学公式从能量函数定义的概率分布中生成样本。因此，样
本由能量函数 Eθ (x) 采样生成，而不是仅由前馈网络生成。
1.3 基于深度能量模型的训练
希望 由能量 函数 E定义 的分布对数据分布 pD 进行建模 ，通过最小化负对数似然 LML (θ) = Ex∼pD (−logpθ (x))来实现，其中 −logpθ (x) = Eθ (x)− logZ(θ) 。该目标函数的梯度为： ∇θLML (θ) = Ex+ ∼pD (∇θ Eθ (x+ ))− Ex− ∼pθ (∇θ Eθ (x− )) 。 （3）
公式（3）直观地表明这个梯度减小了正样本 x+ 的能量，同时增加了来自模型 pθ 的能 量，通过公式（2）来生成 qθ 作为 pθ 的近似值。则公式（3）可近似为：
∇θLML (θ) ≈ Ex+ ∼pD (∇θ Eθ (x+ ))− Ex− ∼qθ (∇θ Eθ (x− )) 。 （4）
1.4 低剂量 CT 深度能量成像模型
本文用以下成像模型来表示低剂量 CT 图像的成像过程：


= argm s. t. y = Ax + μ ,

（5）

其中 μ 是随机噪声。低剂量 CT 观测数据 y 不是直接由 p(x) 采样而来，而是由后验数据分 布 p(x|y) 采样而来。基于贝叶斯法则得 p(x|y) = p(y|x)p(x)/p(y) ，结合公式（1）得到公式（6）：

（6）

其中 p (y|t ) 由数据模型给出，并由事先已知的有关真实模型参数的信息的先验模型给出， 超参数 λ 平衡了先验和数据保真度之间的权衡。结合能量函数定义 pθ (x) = exp (−Eθ (x))/Z(θ) 和公式（6），得出基于深度能量模型的低剂量 CT 成像模型公式：



其次，采用轮换更新方法来求公式（7）。将公式（7）分为两个子迭代步骤：


（7）

（8）

（9）

其中第 2 项为解决不适定逆问题而引入的正则化项， β = β(λ) 是正则化参数，与 λ 相关。 值得注意的是，公式（9）的第 1 项 ∥y− Ax∥2 很难直接最小化，采用可分离二次迭代（SQS） 算法[32]从凸问题中获得CT 重建数据：

其中 AT 是 A转置即反投影，1 是一个全 1 向量。式中的除法指的是逐分量运算符，它使得 算法迭代快速。
首先初始化 z1 为标准高斯噪声向量，然后对 t 和 t 分别进行轮换更新，假设其中一项 已知来求解另外一项。此外，为了加速交替更新过程，应用了先前下降方向的 Nesterov 动 量[33] 。动量项可以是 wt+1 = t+1 + Y (t+1 − t ) ，其中 Y 是介于 0 和 1 之间的松弛因子。
综上所述，表 1 给出的算法 1 展示了 EBM-LDCT 算法的重建流程，具体来说，通过对 t 、 t 和 wt 进行迭代更新来完成。首先，将随机产生的标准高斯噪声输入网络以获得能量 标量，之后采用郎之万动力学和SQS 轮换迭代更新以获得重建结果。
结合数据测量拟合项约束下的郎之万动力学采样过程如图 1 所示。
1.5 卷积神经网络结构
He 等[17]在2016 年首次提出了ResNet 结构，该结构的使用能够有效缓解梯度消失问题

从而训练更深的网络结构。本研究采用 ResNet 结构作为标准的特征提取器，提出的低剂 量 CT 深度能量模型的训练流程和网络结构如图 2 所示。
表 1 EBM-LDCT 算法的重建过程
Table 1 The reconstruction process of the EBM-LDCT algorithm

算法 1
初始化：输入数据 ˜(x)0 ∼ N(0; I) ，动量项 w0 ，参数 β , 松弛因子 Y ，迭代次数 T
开始循环： t = 1; 2; ··· ; T
设定 zt ∼ N(0; I)
根据公式 更新： t−1 − ∇xEθ
最后更新： wt = xt + Y(xt − xt−1)
t t
根据公式（8）更新˜ : ˜t = w˜t−1 −
← w
结束

在网络中首先对输入数据进行 3 × 3 的 卷积以达到升维目的，并采用6 个Resblock 对数据逐步提取深度语义信息，其卷积核大 小均为 3 × 3，前 4 个 Resblock 会进行步长 为 2 的下采样运算，第 6 个 Resblock 后的 全连接层能够对提取的特征进行线性求和， 最终得到网络输出的能量标量。
1.6 评价指标
为了评估重建图像的质量，选择平均绝


xt
xT−1
x0
x(^)
xT

x1
x0 初始数据
xt 中间迭代数据 xT 最终采样数据 x(^) 真实数据分布 朗之万采样
图 1 郎之万动力学采样迭代过程
Fig.1 Iteration process of the
Langevin dynamics sampling

对误差（mean absolute error，MAE）、峰值信噪比（peak signal to noise ratio，PSNR）

朗之万采样 优化函数
采样数据 x-
（a）训练流程
(256)
F(x)+x
3×3 3×3 3×3 3×3 Conv ResBlocklResBlock2 ResBlock3
(512)
输出能量标量Eθ (x)
3×3 ResBlock6
（b）网络结构
Conv
F(x) R(L)e(e)LU(ak)y x
Conv
x
ResBlock
ResBlock4
(512)

3×3 ResBlock5
(64) (64) (128)
原始教据 x+
输入~数据
x
(512)
3×3
×T

图 2 模型训练流程与网络结构
Fig.2 Model training process and network structure

和结构相似性指数（structural similarity，SSIM）3 个指标进行定量评估。在统计学中， MAE 是对表达相同现象的成对观察之间的误差的度量，它的定义为：
MAE = （11）
其中 N 是重建图像像素的个数，MAE 的值越趋近于 0 说明重建的图像越接近原图。
PSNR 度量描述了信号的最大功率与噪声功率之间的关系。更高的 PSNR 意味着更好的 图像质量。 x 和 对应重建图像和真实图像，PSNR 表示为：
PSNR = 20lg （12）
SSIM 用于衡量原始图像和重建图像之间的相似性，它从 3 个方面进行评估：亮度、对 比度和结构相关性。SSIM 值在 0 和 1 之间，1 表示两个图像相等的情况，越高的SSIM 值表 示更好的结构相似性。SSIM 的定义为：

其中 μx 和 σx(2)分别是是 x 的均值和方差， σxˆ(x) 是 x 和 协方差。 c1 和 c2 是用于维持平衡的常
数。为了保证本文所提模型和算法的重复可实现性，我们已经将代码上传至Github 网站：
2 实验结果与分析
2.1 CIRS 仿真数据集
本文使用了两个训练数据集来测试所提出的模型，首先采用的是 CIRS 数据集。它是 从 GE Discovery HD750CT 系 统 中 获 得 的 一 套 拟 人 CIRS 体 模 的 高 质 量 CT 体 积 （（512 × 512 × 100）mm3 ）体素信息，体素尺寸（（0.78 × 0.78 × 0.625）mm3），其中管电流值设 置 600mAs 以保证 150mAs 以下低剂量模拟良好的图像质量。在重建阶段中，通过将噪声等 级为 b = 5× 104 的泊松噪声加入到正弦投影数据中，可以模拟相应的低剂量 CT 图像。实验 给出的重建结果及误差图从左至右分别说明如下，（a1）和（a2）是参考图像，（b1）和（b2）是 FBP 方法重建结果，（c1）和（c2）是 TV 方法重建结果，（d1）和（d2）是K-SVD 方法重建结 果，（e1）和（e2）是 RED-CNN 方法重建结果，（f1）和（f2）是 EBM-LDCT 方法重建结果。
表 2 使用CIRS 数据集的多种低剂量 CT 重建方法定量结果比较
Table 2 Comparison of quantitative results of multiple LDCT
reconstruction methods using CIRS dataset

方法
MAE
PSNR/dB
SSIM/10-2
FBP（Ramp-filter）
9.48 ± 0.64
40.67 ± 0.59
94.36 ± 0.87
TV
7.64 ± 0.41
42.12 ± 0.35
96.58 ± 0.51
K-SVD
7.01 ± 0.19
42.86 ± 0.34
97.01 ± 0.36
RED-CNN
6.74 ± 0.07
41.76 ± 0.12
97.47 ± 0.16
EBM-LDCT
6.29 ± 0.23
42.73 ± 0.20
97.81 ± 0.27

从表 2 的实验结果可以看出，在本研究进行的对比实验中，EBM 用于 CT 任务得到的重 建图像质量是最好的，并产生较高的 MAE 和 SSIM。可见与端到端的成像方法相比，无监督 生成方法重建的低剂量 CT 图像具有更好的细节与结构相似性。
为了直观地说明重建性能，图 3 描绘了由红色矩形标记的放大感兴趣区域（region of interest，ROI）。观察图 3 的（b1）和（b2）中的 FBP 重建结果，我们可以看到 FBP 重建 导致低剂量 CT 图像严重退化，噪声和伪影明显放大。作为传统方法，图 3 中（c1）和（c2） 表现了 K-SVD 能够生成具有视觉上更平滑外观的重建图像。但是一些微小结构可能会被平 滑从而导致组织对比度降低。此外，还可以看出深度学习方法有效地降低了噪声，它们提 高了降噪效果并抑制了大多数伪影。然而，它们对细节和纹理信息的保存不完整。

（a1） （b1） （c1） （d1） （e1） （f1）

（a2） （b2） （c2） （d2） （e2） （f2）
图 3 不同方法的 CIRS 挑战数据的重建结果
Fig.3 Reconstruction results of CIRS challenge data by different methods
比较图 3 中（f1）和（f2）的结果，可见所提 EBM-LDCT 方法在噪声伪影抑制和组织特 征保留方面实现了最佳性能。相应的残差图像如图4 所示，EBM-LDCT 在边缘与细节方面可 以实现比其他算法更好的重建精度。
1.0 0.8 0.6 0.4 0.2 0
图 4 原始 CT 图像与不同算法重建的 CT 图像之间的残差图
Fig.4 The residual plot of the reference CT images and the CT
images reconstructed by different algorithms

为了更好地说明EBM-LDCT 去除噪声的有效性，图 5 绘制了通过图 3 中（a1）的红线的 一维线强度分布图，它比较来自不同方法重建的 CT 图像的相同线强度分布。通过观察发现， 所提方法中的线强度分布与正常剂量 CT 图像中的线强度分布最相似，比较证明了所提出的 重建方法在边缘保留方面优于其他重建算法。

图 5 在图 3（a1）中红线一维数值强度分布（对比所有实验结果）
Fig.5 The one-dimensional numerical intensity distribution of the red
line in Fig.3 (a1) (comparison of all experimental results)

2.2 AAPM 挑战数据集
本研究使用的第 2 个数据集是 Mayo Clinic 提供的 AAPM 低剂量 CT 挑战赛的人体 CT 扫 描重建模拟数据[34] 。数据包括来自 10 位患者在高剂量 CT 下的扫描数据，其中 9 位患者数 据用于训练，另外 1 位患者数据用于评估。本研究使用单切片厚度进行重建，每张图像的 大小为 512 × 512 像素。CT 成像数据采集系统带有 1000 个采集角度，放射源到中心轴的距 离为 500mm，探测器到中心轴的距离为 500mm，成像扫描方式为二维扇形光束。
因为AAPM 数据集比 CIRS 数据集更为复杂，相对来说数据的噪声更大。从表 3 可知， EBM-LDCT 取得了最高的 PSNR 值 39.56 dB。观察图 6 可以看出，得益于 RED-CNN 方法优秀 的去噪能力，在噪声较大的时候，RED-CNN 能有效地抑制大噪声，但是却将一些图像细节 平滑化了。观察残差图 7 可知，K-SVD 算法性能优秀，RED-CNN 能较好地去除边缘突出处
表 3 使用AAPM 数据集的多种 LDCT 重建方法定量结果比较
Table 3 Comparison of quantitative results of multiple LDCT
reconstruction methods using AAPM dataset

方法
MAE
PSNR/dB
SSIM/10-2
FBP
（Ramp-filter）
67.69 ± 13.16
28.68 ± 2.03
44.43 ± 7.18
TV
29.65 ± 4.60
34.98 ± 1.32
86.10 ± 3.48
K-SVD
21.00 ± 1.04
35.68 ± 2.26
81.98 ± 2.41
RED-CNN
16.97 ± 1.02
39.37 ± 1.87
94.78 ± 0.56
EBM-LDCT
18.75 ± 1.31
39.56 ± 0.54
94.16 ± 0.68

（a1） （b1） （c1） （d1） （e1） （f1）

（a2） （b2） （c2） （d2） （e2） （f2）
图 6 不同方法的 AAPM 数据集的重建结果
Fig.6 Reconstruction results of AAPM challenge data by different methods

1.0 0.8 0.6 0.4 0.2 0
图 7 使用AAPM 数据集的原始 CT 图像与不同算法重建结果之间的残差图
Fig.7 The residual plot of the original CT images using AAPM
dataset and the CT images reconstructed by different
algorithms

Pixel localtion
图 8 在图 6（a1）中红线一维数值强度分布（对比所有实验结果）
Fig.8 The one-dimensional numerical intensity distribution of the red
line in Fig.6 (a1) (comparison of all experimental results)

的噪声，EBM-LDCT 则更注重区域块的整体去噪。根据图 8 可知，通过观察实验重建结果的 一维强度图，EBM-LDCT 依旧表现出了良好的去噪效果。
3 结束语
本文提出了一种基于深度能量模型的低剂量 CT 重建方法，该方法通过将深度能量模 型、郎之万动力学和低剂量 CT 成像模型结合，通过深度能量模型对低剂量 CT 概率密度函 数建模，并用郎之万动力学采样公式，对复杂的概率密度函数采样，并且有效地训练神经 网络，生成了高质量的重建图像，体现了该方法的优越性。但是，郎之万动力学采样的时 间往往比端对端的模型需要更多的时间，更有效的对数据概率密度采样将是我们以后提升 的方向。
参考文献
[1] BRENNER D J, HALL E J. Computed tomography: An increasing source of radiation exposure[J]. New England Journal of Medicine, 2007, 357(22): 2277-2284.
[2] KALRA M K, MAHER M M, TOTH T L, et al. Strategies for CT radiation dose optimiza- tion[J]. Radiology, 2004, 230(3): 619-628.
[3] 韩泽芳, 上官宏, 张雄, 等. 基于深度学习的低剂量 CT 成像算法研究进展[J]. CT 理论与应用研究, 2022, 31(1): 118-136. DOI:10.15953/j.1004-4140.2022.31.01.14.
HAN Z F, SHANGGUAN H, ZHANG X, et al. Research progress of low-dose CT imaging algorithm based on deep learning[J]. CT Theory and Applications, 2022, 31(1): 118-136. DOI: 10.15953/j.1004-4140.2022.31.01.14. (in Chinese).
[4] TIAN Z, JIA X, YUAN K, et al. Low-dose CT reconstruction via edge-preserving total variation regularization[J]. Physics in Medicine & Biology, 2011, 56(18): 5949.
[5] MA J, LIANG Z, FAN Y, et al. A study on CT sinogram statistical distribution by information divergence theory[C]//2011 IEEE Nuclear Science Symposium Conference Record. IEEE, 2011: 3191-3196.
[6] ELBAKRI I A, FESSLER J A. Statistical image reconstruction for polyenergetic X-ray computed tomography[J]. IEEE Transactions on Medical Imaging, 2002, 21(2): 89-99.
[7] ZHANG K, ZUO W, CHEN Y, et al. Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[J]. IEEE Transactions on Image Processing, 2017, 26(7): 3142-3155.
[8] LEFKIMMIATIS S. Universal denoising networks: A novel CNN architecture for image denoising[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018: 3204-3213.
[9] ZHANG K, ZUO W, ZHANG L. FFDNet: Toward a fast and flexible solution for CNN-based image denoising[J]. IEEE Transactions on Image Processing, 2018, 27(9): 4608-4622.
[10] CHEN H, ZHANG Y, ZHANG W, et al. Low-dose CT via convolutional neural network[J]. Biomedical Optics Express, 2017, 8(2): 679-694.
[11] CHEN H, ZHANG Y, KALRA M K, et al. Low-dose CT with a residual encoder-decoder convolutional neural network[J]. IEEE Transactions on Medical Imaging, 2017, 36(12): 2524-2535.
[12] WÜRFL T, GHESU F C, CHRISTLEIN V, et al. Deep learning computed tomography[C]// International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Cham, 2016: 432-440.
[13] CHENG L, AHN S, ROSS S G, et al. Accelerated iterative image reconstruction using a