Zero-Shot 低剂量CT图像去噪的扩散概率先验

论文链接：https://arxiv.org/abs/2305.15887

项目链接：https://github.com/DeepXuan/Dn-Dp

Abstract

低剂量CT图像去噪是医学图像计算中的一项关键任务。近年来，基于监督的深度学习方法在这一领域取得了重大进展。然而，这些方法通常需要对低剂量和正常剂量的CT图像进行训练，这在临床环境中很难获得。现有的基于无监督深度学习的方法通常需要使用大量低剂量CT图像进行训练，或者依靠专门设计的数据采集流程来获得训练数据。为了解决这些限制，我们提出了一种新的无监督方法，该方法在训练期间仅使用正常剂量的CT图像，从而实现对低剂量CT图像的Zero-Shot去噪。我们的方法利用了扩散模型，一个强大的生成模型。我们首先训练一个级联无条件扩散模型，该模型能够从低分辨率到高分辨率生成高质量的正常剂量CT图像。级联结构使得高分辨率扩散模型的训练更加可行。随后，我们将低剂量CT图像作为似然引入扩散模型的反向过程中，结合扩散模型提供的先验，迭代求解多个最大后验(MAP)问题实现去噪。此外，我们提出了自适应调整MAP估计中平衡似然和先验系数的方法，允许适应低剂量CT图像中的不同噪声水平。我们在不同地区不同剂量水平的低剂量CT数据集上测试了我们的方法。结果表明，我们的方法优于目前最先进的无监督方法，并超过了几种基于监督的深度学习方法。

I. INTRODUCTION

低剂量计算机断层扫描(CT)因其显著降低辐射剂量而受到越来越多的关注。虽然低剂量CT降低了患者暴露于X射线的风险，但它通常会导致质量较差的重建图像和明显的噪声，这可能会妨碍准确的诊断。为了提高低剂量CT的图像质量，人们开发了多种算法，包括正弦图滤波[1]、迭代重建[2]和图像去噪[3]。在这些方法中，低剂量CT图像去噪应用最广泛，因为它不涉及重建过程，可以应用于不同的CT成像系统。

近十年来，在深度学习技术快速发展的推动下，深度神经网络已成为低剂量CT图像去噪的有力工具，并取得了显著的成功。通常，基于深度学习的方法以监督的方式学习从低剂量到正常剂量CT图像的端到端映射。最初，卷积神经网络(CNN)是低剂量CT图像去噪最常用的方法[3]。为了提高去噪图像的视觉质量，生成对抗网络(GANs)也被用于低剂量CT图像去噪[4]，[5]。此外，基于Transformer的网络进一步增强了低剂量CT图像去噪的性能[6]，[7]。

然而，监督深度学习算法的成功在很大程度上依赖于大量的成对训练数据[8]。不幸的是，由于遵守“尽可能低的合理可行”原则，在临床环境中获取此类数据具有挑战性[9]。因此，开发能够在不需要大量标记数据的情况下利用深度神经网络高性能的方法至关重要。

一些现有的方法试图使用不同的策略来解决这个问题。一类方法依赖于通过特殊获取方案获得的训练数据。Yuan等人提出Half2Half，利用配对的半剂量CT图像训练神经去噪器[10]。Wu等人开发了一种利用奇偶投影重建的图像对训练低剂量CT图像去噪网络的方法[11]。另一类方法进一步放宽了训练数据的条件，但仍然需要大量的低剂量CT图像进行训练。例如，Du等人提出从噪声图像中学习不变表示并重建干净的观察结果[12]。Liu等人采用可逆神经网络模拟配对正常剂量和低剂量CT图像，并进行非配对训练[8]。Niu等人提出了一种基于相似度的无监督深度去噪方法来处理CT图像中的相关噪声[13]。

与上述无监督方法不同，本文提出的方法不需要配对正常剂量和低剂量CT图像，也不需要任何低剂量CT图像作为训练集。我们仅使用正常剂量的CT图像训练去噪扩散概率模型[14]，[15] (下文简称扩散模型)，然后利用预训练扩散模型中嵌入的先验信息实现Zero-Shot低剂量CT图像去噪。扩散模型定义了一个恒定的正向过程，并通过迭代求解反向过程来实现图像生成。扩散模型不仅实现了最先进的生成结果，而且还提供了生成过程中数据分布的分析表示，这是GAN等模型无法实现的。由于扩散模型具有强大的生成和表征能力，已有一些研究探索了其先验在图像重建和恢复等下游任务中的应用[16]-[19]。然而，这些方法大多要么忽略了噪声的存在[19]，要么只考虑简单的噪声场景[16]-[18]。在我们的方法中，我们首先提出将扩散先验纳入最大后验(MAP)框架，以解决低剂量CT图像去噪问题。此外，我们引入了一种动态调整先验和似然权重的策略来处理复杂的低剂量CT噪声，这些噪声在不同的CT切片上表现出不同的水平。

首先，我们使用正常剂量的CT图像训练扩散模型。先前的研究强调了直接生成尺寸为512×512的高分辨率图像的挑战。采用级联生成方法来解决这个问题[20]-[22]。具体来说，我们从随机噪声中生成低分辨率图像，然后通过对低分辨率图像的调节生成大小为512×512的高分辨率图像。然后，我们提出了一种将预训练扩散模型纳入MAP去噪框架的算法。在扩散模型生成过程的每一次迭代中，我们引入低剂量CT图像并解决MAP估计问题，确保生成的图像与输入的低剂量CT图像具有良好的似然性。最终，该扩散模型可以生成输入低剂量CT图像的高质量对应物，实现低剂量CT图像去噪。此外，考虑到不同的低剂量CT切片表现出不同的噪声水平，我们设计了两种自适应策略来改进MAP估计中的系数，以平衡似然和先验。该算法利用改进后的系数在中间时间步恢复去噪过程，可以自适应处理不同噪声水平的低剂量CT图像。

我们的贡献可以概括如下：

1. 我们提出了在大规模CT图像数据集上以级联方式训练扩散模型，该模型可以从随机噪声中逐步生成真实的低分辨率和高分辨率正常剂量CT图像。
2. 通过求解扩散模型反向过程中的多个MAP问题，提出了一种基于扩散先验的完全无监督Zero-Shot低剂量CT图像去噪算法。
3. 引入了两种自适应策略来平衡MAP估计系数的似然性和先验性，以处理不同水平的低剂量CT噪声。
4. 我们的方法在视觉效果和指标方面都表现出色，优于最先进的无监督方法，甚至超过了几种基于监督的深度学习方法。

II. METHOD

A. 基于扩散模型的NDCT图像生成

我们采用两步级联扩散模型来生成正常剂量的CT图像。第一个模型是无条件扩散模型，它从随机噪声中生成低分辨率正常剂量CT图像。第二种模型是条件扩散模型，以低分辨率图像为条件生成尺寸为512×512的高分辨率正常剂量CT图像。整个模型是无监督的，因为它只使用正常剂量的CT图像，没有标签进行训练。

1) 无条件扩散模型(Unconditional Diffusion Model)：给定一组数据$x_0\sim q(x_0)$，扩散模型定义了一个固定的正向过程，并旨在将其反转，以实现从随机噪声中生成图像。前向过程$q(x_{1:T}|x_0)$被定义为一个马尔可夫链，其中每个时间步的概率分布只依赖于前一个时间步。具体来说，我们有：
$$\begin{array}{rl}q\left(x_{1:T}\mid x_0\right)& :=\prod _{t=1}^{T}q\left(x_t\mid x_{t-1}\right),q\left(x_t\mid x_{t-1}\right)\\ & :=\mathcal{N}\left(x_t;\sqrt{1-{\beta }_t}{x}_{t-1},{\beta }_{t}I\right),\end{array}\text{\hspace{1em}(1)}$$
式中，$x_t$为时间步长t的数据，$\mathcal{N}\left(\cdot ;\mu ,\Sigma \right)$表示高斯分布，平均值$\mu$和协方差矩阵$\Sigma$， ${\beta }_1,\dots ,{\beta }_T$是固定方差表。前向过程逐步向$x_0$中加入噪声[15]。根据(1)，我们可以推导出以下两个条件分布：
$$\begin{array}{r}q\left(x_t\mid x_0\right)=\mathcal{N}\left(x_t;\sqrt{{\bar{\alpha }}_t}{x}_0,\left(1-{\bar{\alpha }}_t\right)I\right),\\ q\left(x_{t-1}\mid x_t,x_0\right)=\mathcal{N}\left(x_{t-1};{\tilde{\mu }}_t\left(x_t,x_0\right),{\tilde{\beta }}_{t}I\right),\end{array}\text{\hspace{1em}(2)(3)}$$
其中：
$$\begin{array}{rl}& {\bar{\alpha }}_t:=\prod _{i=1}^t(1-{\beta }_i),\\ & {\tilde{\mu }}_t\left(x_t,x_0\right):=\frac{\sqrt{{\bar{\alpha }}_{t-1}}{\beta }_t}{1-{\bar{\alpha }}_t}{x}_0+\frac{\sqrt{{\alpha }_t}\left(1-{\bar{\alpha }}_{t-1}\right)}{1-{\bar{\alpha }}_t}{x}_t,\\ & {\tilde{\beta }}_t:=\frac{(1-{\bar{\alpha }}_{t-1})}{1-{\bar{\alpha }}_t}{\beta }_t.\end{array}\text{\hspace{1em}(4)(5)(6)}$$
将反向过程定义为具有可学习参数$\theta$的马尔可夫链，
$$p_{\theta }\left(x_{0:T}\right):=p\left(x_T\right)\prod _{t=1}^T{p}_{\theta }\left(x_{t-1}\mid x_t\right),\text{\hspace{1em}(7)}$$
其中：
$$\begin{array}{rl}& p\left(x_T\right)=\mathcal{N}\left(x_T;0,I\right),\\ & p_{\theta }\left(x_{t-1}\mid x_t\right):=\mathcal{N}\left(x_{t-1};{\mu }_{\theta }\left(x_t,t\right),{\sigma }_t^{2}I\right).\end{array}\text{\hspace{1em}(8)(9)}$$
这里，${\mu }_{\theta }\left(x_t,t\right)$一般由U-nets预测，${\sigma }_t^2$定义为$\frac{1-{\bar{\alpha }}_t-1}{1-{\bar{\alpha }}_t}{\beta }_t$。训练损失是通过最小化负对数似然$-\log p_{\theta }\left(x_0\right)$的变分界得到的[14]，

2) 条件扩散模型(Conditional Diffusion Model)：给定图像对$(x_0,c)$的数据集，条件扩散模型与无条件模型保持相同的正向过程，其中正向过程与$c$无关，遵循(1)。而反向过程以$c$为附加条件，可定义为$p_{\theta }\left(x_{0:T}\mid c\right):=p\left(x_T\right){\prod }_{t=1}^T{p}_{\theta }\left(x_{t-1}\mid x_t,c\right)$，其中
$$p_{\theta }\left(x_{t-1}\mid x_t,c\right):=\mathcal{N}\left(x_{t-1};{\mu }_{\theta }\left(x_t,t,c\right),{\sigma }_t^2\right)\text{\hspace{1em}(10)}$$
3) 级联扩散模型用于NDCT图像生成(Cascaded Diffusion Model for NDCT Image Generation)：我们建议使用级联框架生成高分辨率正常剂量CT图像，因为先前的研究表明直接生成高分辨率图像存在挑战[20]-[22]。给定一组正常剂量的CT图像$x_0^{hr}\in {\mathbb{R}}^{512\times 512}$，我们通过对$x_0^{hr}$进行k的降采样得到图像对$(x_0^{hr},x_0^{lr})$，其中$x_0^{lr}\in {\mathbb{R}}^{512\times 512}$。我们首先定义了一个前向过程参数 { α ˉ t l r } t = 1 , … , T \{\bar{\alpha}^{lr}_t\}_{t=1,…,T} {αˉtlr​}t=1,…,T​和反向过程参数${\theta }^{lr}$，用于从随机噪声中生成低分辨率正常剂量CT图像。第二种扩散模型是以低分辨率CT图像为条件，实现正常剂量CT图像超分辨率的条件扩散模型。其正向过程参数为 { α ˉ t s r } t = 1 , … , T ∞ \{\bar{\alpha}_{t}^{sr}\}_{t=1,\ldots,T}^{\infty} {αˉtsr​}t=1,…,T∞​和反向过程参数为${\theta }^{sr}$。通过级联两种扩散模型，首先由随机噪声生成低分辨率正常剂量CT图像，然后将生成的低分辨率CT图像作为超分辨率扩散模型的条件生成高分辨率正常剂量CT图像。总体生成过程如图1所示。

B. 基于扩散先验的低剂量CT图像去噪

在本节中，我们提出了一种称为Dn-Dp(扩散先验去噪)的算法，用于使用扩散先验对低剂量CT图像进行去噪。我们首先在一般情况下描述该算法，然后利用我们的预训练扩散模型将其应用于低剂量CT图像去噪的具体问题。

1）使用扩散先验的迭代去噪：在图像去噪问题中，损坏可以建模如下：
$$y_0=x_0+n,\text{\hspace{1em}(11)}$$
其中$y_0$是观察到的噪声图像，$x_0$是我们要估计的底层干净图像。为了实现这一点，我们利用一个训练有素的扩散模型，该模型使用前向过程参数 { α ˉ t } t = 1 , … , T \{\bar{\alpha}_t\}_{t=1,…,T} {αˉt​}t=1,…,T​准确捕获$x_0$的先验分布和反向过程参数$\theta$。目标是使用扩散先验有效地去噪。具体方法概述如下。

在扩散模型中，图像$x_t$在时间步长$t=1,\dots ,T$根据正向过程，由$x_0$生成：
$$x_t=\sqrt{{\bar{\alpha }}_t}{x}_0+\sqrt{1-{\bar{\alpha }}_t}{ϵ}_t,\text{\hspace{1em}(12)}$$
其中${ϵ}_t\sim \mathcal{N}\left(0,I\right)$取自正态分布。为了在每个时间步$t$上引入输入噪声图像$y_0$的似然，根据之前的工作[18]我们使用与(12)相同的${ϵ}_t$从$y_0$中生成$y_t$(对于$t=1,\dots ,T$)：
$$y_t=\sqrt{{\bar{\alpha }}_t}{y}_0+\sqrt{1-{\bar{\alpha }}_t}{ϵ}_t.\text{\hspace{1em}(13)}$$
将式(11)和式(12)代入式(13)，设${\bar{\alpha }}_0=1$，我们得到$t=0,\dots ,T$时的方程：
$$y_t=x_t+\sqrt{{\bar{\alpha }}_t}n.\text{\hspace{1em}(14)}$$
式(14)表明，在每个时间步长$t$，存在一个新的噪声$\sqrt{{\bar{\alpha }}_t}n$的去噪问题。由于${\bar{\alpha }}_t$随着$t$的增加而减小，并且${\bar{\alpha }}_t\approx 0$[14]，我们可以使用${\hat{x}}_T=y_T$作为初始化。然后，对于$t=T\dots ,1$，我们可以将扩散先验$p_{\theta }\left(x_{t-1}\mid x_t\right)$纳入MAP框架[23]来求解$x_{t-1}$：
$$\begin{array}{rl}{\hat{x}}_{t-1}& =\arg \underset{x_{t-1}}{\min }[\Vert x_{t-1}-y_{t-1}{\Vert }^2\\ & -{\lambda }_{t-1}\log p_{\theta }(x_{t-1}\mid {\hat{x}}_t)],\end{array}\text{\hspace{1em}(15)}$$
其中${\lambda }_t$是在时间步长t处权衡似然和先验的系数。在扩散模型的反向过程中，$p_{\theta }\left(x_{t-1}\mid {\hat{x}}_t\right)$服从高斯分布，其均值为${\mu }_{\theta }\left({\hat{x}}_t,t\right)$，对角协方差矩阵为${\sigma }_t^{2}I$。当$t>1$时，${\sigma }_t>0$[14]，则有：
$$\begin{array}{rl}{\hat{x}}_{t-1}& =\arg \underset{x_{t-1}}{\min }[\Vert x_{t-1}-y_{t-1}{\Vert }^2\\ & +\frac{{\lambda }_{t-1}}{{\sigma }_t}\Vert x_{t-1}-{\mu }_{\theta }({\hat{x}}_t,t){\Vert }^2],t>1.\end{array}\text{\hspace{1em}(16)}$$
当$t=1$时，${\sigma }_t=0$，表明扩散的最终生成步长是确定的：
$${\hat{x}}_{t-1}={\mu }_{\theta }\left({\hat{x}}_t,t\right),t=1.\text{\hspace{1em}(17)}$$
式(16)是一个具有封闭解的凸优化问题。通过求解(16)并与(17)结合，我们可以得到每次迭代中${\hat{x}}_{t-1}$的解如下：
$${\hat{x}}_{t-1}=\{\begin{array}{ll}\frac{y_{t-1}+{\lambda }_{t-1}{\mu }_{\theta }\left({\hat{x}}_t,t\right)/{\sigma }_t}{1+{\lambda }_{t-1}/{\sigma }_t},& t>1\\ {\mu }_{\theta }\left({\hat{x}}_t,t\right),& t=1\end{array}.\text{\hspace{1em}(18)}$$
同理，若预训练扩散模型以$c$为条件，则解为：
$${\hat{x}}_{t-1}^c=\{\begin{array}{ll}\frac{y_{t-1}+{\lambda }_{t-1}{\mu }_{\theta }\left({\hat{x}}_t,t,c\right)/{\sigma }_t}{1+{\lambda }_{t-1}/{\sigma }_t},& t>1\\ {\mu }_{\theta }\left({\hat{x}}_t,t,c\right),& t=1\end{array}.\text{\hspace{1em}(19)}$$
式(18)和(19)说明了使用预训练的无条件或条件扩散模型去噪的过程。在每个时间步，(13)中的$y_{t-1}$有助于校正扩散模型的输出，确保生成图像的准确性。具体算法用算法1表示。

给定低剂量CT图像$y_0^{hr}$和两个参数为 { θ l r , { α ˉ t l r } t = 1 , . . . , T } \{θ^{lr},\{\bar{α}^{lr}_t\}_{t=1,...,T}\} {θlr,{αˉtlr​}t=1,...,T​}的预训练扩散模型和 { θ s r , { α ˉ t s r } t = 1 , … , T } \left\{\theta^{s r},\left\{\bar{\alpha}_t^{s r}\right\}_{t=1, \ldots, T}\right\} {θsr,{αˉtsr​}t=1,…,T​}，我们首先通过$k\times$下采样得到低分辨率低剂量CT图像$y_0^{lr}$。接下来，我们输入$y_0^{lr}$, $c=None$, ${\theta }^{lr},{\left\{{\bar{\alpha }}_t^{lr}\right\}}_{t=1,\dots ,T},{\left\{{\lambda }_t^{lr}\right\}}_{t=1,\dots ,T-1}$代入算法1，得到去噪后的低分辨率CT图像${\hat{x}}_0^{\text{lr }}$，如图2(a)所示。然后，输入$y_0^{hr}$，${\hat{x}}_0^{lr}$，${\theta }^{sr}$， { α ˉ t s r } t = 1 , … , T \left\{\bar{\alpha}_t^{s r}\right\}_{t=1, \ldots, T} {αˉtsr​}t=1,…,T​， { λ t s r } t = 1 , … , T − 1 \left\{\lambda_t^{s r}\right\}_{t=1, \ldots, T-1} {λtsr​}t=1,…,T−1​ (其中${\hat{x}}_0^{lr}$为条件)转化为算法1。得到去噪后的高分辨率CT图像${\hat{x}}_0^{hr}$，如图2(b)所示。图2中的蓝色块表示我们去噪算法的一次迭代。

2) λ的自适应细化：在式1中，超参数 { λ t } t = 1 , … , T − 1 \left\{\lambda_t\right\}_{t=1, \ldots, T-1} {λt​}t=1,…,T−1​用于平衡似然和先验之间的权衡。根据经验，噪音水平较低λ需要更小的值。由式(14)可知，时间步长$t$处的噪声电平由$$\sqrt{{\bar{\alpha }}_t}n$$给出，其中当$t$从0增加到T时，${\bar{\alpha }}_t$从1减小到接近0。因此，我们暂定${\lambda }_t$为${\bar{\alpha }}_t$与常数${\lambda }_0$的乘积，
$${\lambda }_t={\lambda }_0\cdot \sqrt{{\bar{\alpha }}_t}.\text{\hspace{1em}(20)}$$
我们把这种策略称为ConsLam。然而，不同的低剂量CT切片通常表现出不同的噪声水平，使得固定的${\lambda }_0$很难达到每个图像的最佳去噪结果。为了对不同噪声水平的低剂量CT图像进行自适应降噪，首先提出了估计单幅低剂量CT图像的噪声水平的方法。

由于低剂量CT噪声的非解析性，从单幅低剂量CT图像中直接估计噪声水平是困难的。因此，我们首先使用固定的${\lambda }_0$对低剂量CT图像进行去噪，得到初始去噪结果${\hat{x}}_0$。然后将估计的噪声计算为原始图像$y_0$与去噪图像${\hat{x}}_0$之间的绝对差，
$$\hat{n}=|y_0-{\hat{x}}_0|.\text{\hspace{1em}(21)}$$
虽然使用固定的${\lambda }_0$获得的估计$\hat{n}$可能不是最优的，但它仍然提供了输入低剂量CT图像$y_0$中的噪声的可靠表示。我们基于$\hat{n}$调整${\lambda }_0$，然后从中间某一步恢复去噪算法，以获得更精确的去噪。具体来说，我们提出了两种精炼${\lambda }_0$的策略。

第一种策略是使用估计噪声的标准差作为噪声级来动态确定${\lambda }_0$。它被表示为${\lambda }_0^{ada}$，
$${\lambda }_0^{ada}=a\cdot std(\hat{n})+b,\text{\hspace{1em}(22)}$$
其中std(·)为图像中所有像素的标准差，a和b为人工选择的参数。我们将此策略称为AdaLam-I。

第二种策略考虑到低剂量CT图像的不同区域表现出不同程度的噪声。为了解决这个问题，我们在不同位置引入了一个不均匀的${\lambda }_0$，得到了一个矩阵${\Lambda }_0^{ada}$，
$${\Lambda }_0^{ada}=c\cdot \hat{n},\text{\hspace{1em}(23)}$$
其中c是手动选择的系数。为了将矩阵${\Lambda }_0^{ada}$代替标量${\lambda }_0$合并到算法1中，我们执行Hadamard积。我们把第二个策略称为AdaLam-II。在第IV . b节中，我们验证了分别使用这两种策略中的每一种策略以及同时使用它们的有效性。

图3给出了${\lambda }_0$细化的详细过程。值得注意的是，以更大的步长恢复可以从$\hat{n}$中得到更准确的${\lambda }_0^{ada}$或${\Lambda }_0^{ada}$。理论上，在时间步长T处恢复将是最优选择。然而，这种方法会使算法的计算时间增加一倍。因此，为了平衡准确性和计算时间，我们选择只回滚几个步骤。

在我们的实验中，我们将迭代步骤设置为3，从${\hat{x}}_3$恢复，这只会增加推理时间大约10%。值得注意的是，我们只在降噪算法的第二阶段采用了$\lambda$的自适应细化策略，因为常数${\lambda }_0$足以用于低分辨率图像的降噪。

3) 加速：扩散模型依赖于迭代生成过程，其耗时的性质造成了显著的限制。同样的问题也影响着我们基于扩散模型的去噪算法。为了解决这个问题，我们采用了DDIM[15]提出的方法，通过区间抽样加快扩散模型的生成过程。具体地说，我们提取了一个子序列$\tau =[{\tau }_1,\dots ,{\tau }_S]$的长度为S，从[1，…,T]，然后从$x_{{\tau }_s}$到$x_{{\tau }_{s-1}}$进行采样。对于我们的去噪算法，迭代次数按照以下关系递减：
$${\hat{x}}_{{\tau }_{s-1}}=\frac{y_{{\tau }_{s-1}}+{\lambda }_{{\tau }_{s-1}}{\mu }_{\theta }\left({\hat{x}}_{{\tau }_s},{\tau }_s\right)/{\sigma }_{{\tau }_s}}{1+{\lambda }_{{\tau }_{s-1}}/{\sigma }_{{\tau }_s}},s>1.\text{\hspace{1em}(24)}$$
对于加速算法，当t > 1时，式(24)代替(18)。在条件扩散模型的情况下，实现是类似的，将${\mu }_{\theta }\left({\hat{x}}_{{\tau }_s},{\tau }_s\right)$替换为${\mu }_{\theta }\left({\hat{x}}_{{\tau }_s},{\tau }_s,c\right)$。

在典型的生成任务中，加速采样的时间步长通常被均匀地减小。然而，在本文中，我们提出采用一种非均匀策略从2000个时间步长中提取子序列。具体来说，在前500个时间步内，我们每20个时间步提取一个。在随后的1500个时间步中，我们每500个时间步提取一个。最后，我们加入最后的时间步长T，以保证算法初始条件的有效性，即${\hat{x}}_T=y_T$。因此，我们设$\tau =[1,21,41,\dots ,501,1001,1501,2000]$并将总迭代次数从2000减少到29。设计这种非均匀采样方法是因为似然项在较小的时间步长中发挥更突出的作用，而更密集的采样可以更好地利用低剂量CT输入，从而获得更准确的去噪结果。

III. IMPLEMENTATION DETAILS

A. 数据集

我们分别在低剂量 CT 和投影数据 (LDCT-PD) 数据集 [24]、[25] 中的 10 名患者的腹部和胸部正常剂量 CT 图像上预训练两个扩散模型。对于腹部和胸部CT图像，高分辨率图像尺寸均为512×512，低分辨率图像尺寸分别为128×128和256×256。为了评估我们提出的Dn-Dp算法的去噪性能，我们使用了来自2016年AAPM Grand Challenge数据集[3]的4名患者的四分之一剂量腹部CT图像和来自LDCT-PD数据集的4名患者的10%剂量胸部CT图像，以确保这些图像与训练数据不同。

B. 网络结构及参数

对于我们的扩散模型，我们采用来自SR3[26]的U-Net架构，而残差块则来自BigGAN[27]。在条件扩散模型的情况下，使用双三次插值对低分辨率图像进行上采样，然后沿着通道维度与U-Net的先前输出连接。表1提供了初始信道维度(Dim $d$)，信道乘数(Muls { m [ 1 ] , … , m [ k ] } \{m[1],…,m[k]\} {m[1],…,m[k]})的具体数量，残差块的数量(ResBlocks)，使用自注意残差块的分辨率(Res-Attn)，以及dropout rate (dropout)。

我们的扩散模型中使用的条件U-Net架构如图4所示。为了获得无条件版本，只需从模型中删除条件图像c。

对于预训练的两个级联扩散模型，我们将时间步长数设为定义为T的2000。在正演过程中，我们设置${\beta }_1=1\times 10^{-6}$，${\beta }_T=1\times 10^{-2}$，${\beta }_i(1<i<T)$在中间时间步长线性增加。本文提出的去噪算法涉及到(20)、(22)和(23)中的${\lambda }_0$、$a$、$b$和$c$几个超参数，分别通过人工调整来实现腹部和胸部CT图像的去噪。我们在表二中提供了选定的参数。

此外，在正演过程中引入的变化噪声产生的去噪结果显示一致的分布，但像素值不同。为了提高性能，我们计算10个去噪结果的平均值[16]。

IV. EXPERIMENTS

A. 正常剂量CT图像生成

图5展示了由我们预先训练的级联扩散模型生成的高质量正常剂量CT图像。这些生成的图像显示出显著的真实感和多样性，表明我们训练有素的扩散模型对正常剂量CT图像具有很强的先验理解。

B. λ的消融研究

在本小节中，我们将评估(22)和(23)中提出的$\lambda$改进策略的有效性，包括AdaLam-I、AdaLam-II和表示同时使用的AdaLam-I& ii。为了进行消融研究，我们从4例患者的腹部和胸部数据中统一提取了80个CT切片，形成了两个迷你数据集。表III给出了去噪结果的统计度量。度量是在[- 1024,3072]HU窗口内计算的。

从表III的结果可以看出，对于腹部CT图像去噪，AdaLam-I和AdaLam-II的平均PSNR(峰值信噪比)都提高了约0.1dB，两者联合使用效果同样良好。另一方面，对于噪声水平较高的胸部CT图像去噪，与ConsLam相比，AdaLam-I和AdaLam-II在PSNR和SSIM(结构相似性指数测量)方面表现出更显著的增强。此外，与单独使用任何一种相比，使用adalam - i和ii进行胸部CT图像可获得更好的结果。

由于在不同的低剂量CT图像中观察到不同的噪声水平，使用自适应细化策略可以增强去噪结果。通常，低噪声水平的低剂量CT图像更适合较小的${\lambda }_0$，而高噪声水平的图像则适合较大的${\lambda }_0$。图6(a)和图6©分别绘制了两幅腹部低剂量CT图像和两幅胸部低剂量CT图像去噪后图像PSNR与常数${\lambda }_0$的关系。

很明显，在不同的CT图像中使用ConsLam可以获得具有不同${\lambda }_0$值的最佳去噪结果。通过计算自适应${\lambda }_0^{ada}$和${\Lambda }_0^{ada}$，然后恢复去噪算法进行几次迭代，我们可以对不同噪声水平的低剂量CT图像获得更好的去噪效果。如图6(b)和6(d)所示，同时适应AdaLam-I和AdaLam-II产生的去噪性能超过或接近使用ConsLam获得的最佳结果。这说明了为什么采用自适应$\lambda$策略可以增强算法在整个数据集上的性能，如表III所示。

C. 降噪性能比较

基于λ的消融研究，我们已经确定同时使用AdaLam-I和AdaLamII可以为我们的去噪算法产生最佳性能。在接下来的实验中，我们使用Adalam - i&ii进行DnDp。此外，我们评估了两种版本的Dn-Dp: Dn-Dp (Single)，它对输入图像进行一次去噪，以及Dn-Dp (Average)，它对相同的输入图像产生多个去噪结果并对它们进行平均。后一种方法可以提高性能，但需要更多的计算时间。

我们将该方法与各种低剂量CT图像去噪算法进行了比较，包括两种经典的非学习方法:块匹配3D (BM3D)[28]和非局部均值(NLM)[29]。此外，我们还考虑了在低剂量CT图像去噪方面取得优异效果的两个基于监督学习的神经网络:RED-CNN[3]和CTformer[6]，其中RED-CNN基于CNN, CTformer是基于Transformer的模型。由于我们的方法以完全无监督的方式运行，我们还将其与两种基于无监督学习的方法：LIR-IR[12]和Noise2Sim[13]进行了比较，后者是目前最先进的无监督低剂量CT图像去噪方法。

表IV给出了不同方法的定量去噪结果，度量使用窗口[- 1024,3072]HU计算。我们提出的方法在PSNR和SSIM方面显示出有竞争力的定量结果，如表4所示。

具体而言，在腹部CT数据集上，我们的方法优于所有无监督方法。值得注意的是，我们的方法比经典的基于监督的深度学习方法RED-CNN实现了更高的PSNR，并且接近最新的基于Transformer的网络CTformer的性能。这个结果令人惊讶，因为用MSE(均方误差)损失训练的监督网络通常优先考虑更高的PSNR值。然而，对于胸部CT数据集，低剂量CT图像表现出更明显的噪声，基于CNN的监督神经网络(如REDCNN)表现出最好的性能。由于噪声水平较高，无监督算法和有监督算法之间的性能差距变得更加明显。尽管如此，我们的方法仍然是无监督方法中表现最好的，超过了LIR-IR和Noise2Sim。

图7显示了腹部CT图像的去噪结果。

很明显，我们的方法产生的结果与正常剂量CT图像的视觉质量非常相似。通过对比©和(d)可以看出，传统方法无法有效消除低剂量CT中存在的噪声。另一方面，(e)和(f)表明，基于无监督深度学习的方法比传统方法产生更好的结果，但由于其端到端性质，仍然容易丢失正常剂量CT图像中的一些纹理细节。对于有监督的深度学习算法，这个问题在(g)和(h)中变得更加明显，导致明显的过平滑。相比之下，我们的方法利用MAP框架结合先验和似然，在视觉上获得了更真实的结果。

图8显示了胸部CT图像的视觉结果。

在黄色感兴趣区域(ROI)内，比较了不同算法恢复骨结构的能力。由于10%剂量胸部CT的噪声水平较高，可以观察到即使是监督算法也会丢失一些轮廓信息。然而，我们的方法力求在去噪过程中尽可能地保留轮廓。红色ROI内，箭头示实性非钙化肺结节。用BM3D、NLM和LIRIR去噪会使结节模糊，妨碍临床诊断。相比之下，我们的方法和监督方法有效地保持了结节的形状。

D. 讨论

虽然我们的方法在视觉效果和定量指标方面取得了很好的效果，但仍然存在一些值得注意的局限性。我们的方法的主要限制在于模型参数的大小和所需的计算时间。为了训练生成高分辨率正常剂量CT图像的扩散模型，需要级联生成方法和非常大的模型尺寸。尽管使用加速算法将迭代次数从2000次减少到29次，但与端到端网络相比，单个图像去噪的推理时间仍然要高得多。表V给出了所用模型的参数大小以及使用这些模型去噪所需的相应推理时间。其中Time(1)和Times(10)分别表示使用单个NVIDIA GeForce RTX 3090显卡对单个图像和10个图像并行去噪的推理时间。

在对低分辨率图像进行Dn-Dp去噪后，另一种方法是直接训练超分辨率神经网络来获得高分辨率的干净图像。虽然超分辨率网络需要监督训练，但训练对可以通过下采样从正常剂量的CT图像中产生。使用这样的超分辨率网络代替Dn-Dp的第二阶段将不可避免地导致性能损失，因为它没有利用高分辨率低剂量CT图像的信息，但它可以节省一些推理时间并提供合理的去噪结果 。

此外，我们选择使用级联方式生成高分辨率图像的原因是为了便于扩散模型的预训练。然而，这种方法可能会导致性能下降，因为降采样实际上会进一步降低图像质量，导致一些信息丢失。如果我们可以直接训练无条件扩散模型用于高分辨率图像生成，那么我们的算法对于单阶段去噪也是有效的，并且结果可能会更好。然而，训练这样一个高分辨率的无条件扩散模型往往需要更深入的网络、更多的训练数据和更长的训练时间。

V. CONCLUSION

本文提出了一种基于扩散先验的低剂量CT图像去噪算法Dn-Dp。Dn-Dp是一种完全无监督的算法，它只需要用正常剂量的CT图像训练网络。该算法包括训练级联扩散模型来生成正常剂量的CT图像。随后，在扩散模型的反向过程中，我们迭代求解了多个MAP估计问题。这确保了生成的图像与输入的低剂量CT图像之间的高似然性，同时保持了良好的图像质量。为了解决不同低剂量CT图像中不同的噪声水平，我们采用了两种自适应策略来调整系数$\lambda$，以平衡MAP估计中的似然项和先验项。我们还介绍了一种以更精确的λ在中间阶段恢复去噪的方法。我们的方法在定量指标和视觉效果方面表现出色，特别是在四分之一剂量的腹部CT图像去噪方面，在PSNR和SSIM方面甚至超过了监督方法。此外，我们的目标是将基于扩散先验的方法扩展到CT重建算法中，从正弦图域开始。