基于光传输蒙特卡洛方法的机器学习增强光学断层扫描图像重建

《Journal of Biomedical Optics》：Machine-learning-enhanced image reconstruction in optical tomography using the Monte Carlo method for light transport

　　

图1：一个卷积神经网络（记为 $G_{\theta_i}$）的示意图，表示一次深度随机高斯-牛顿迭代。带有5×5卷积核的卷积层，后接LReLU激活函数，用红色箭头标注。带有5×5卷积核的卷积层，后接标量乘法 $\lambda_j$，用蓝色箭头标注。产生的通道数标注在方框内。跳跃连接（灰色箭头）将残差更新量通过最后一个ReLU（绿色箭头标注）投影到正数域。

## 1. 引言

光断层成像（Optical Tomography, OT）是一种生物医学成像技术，通过边界测量的可见红光或近红外光来重建成像目标的光学参数。$^{1,2}$ 这种无创方法能够提供生物目标的独特功能性和结构性信息，其应用例如包括功能性脑成像，$^{3,4}$ 乳腺癌成像，$^{5}$ 甲状腺癌成像，$^{6}$ 类风湿关节炎诊断，$^{7,8}$ 治疗监测，$^{9}$ 以及小动物研究。$^{10}$

光断层成像中的光学参数重建是一个高度病态的逆问题，即测量或建模中的微小误差都可能导致重建图像出现显著误差。$^{1,11}$ OT中的逆问题通常以差分成像问题或绝对成像问题的形式提出。$^{1}$ 在差分成像中，使用两组测量数据来重建两组之间光学参数变化的图像。在绝对成像中，使用一组测量数据来重建光学参数的（绝对）值图像。在本文中，我们考虑光断层成像中吸收系数和散射系数的绝对重建。

OT中的计算图像重建需要在图像重建算法中对成像域中的光传播进行建模。在生物组织中，可以使用辐射传输方程（RTE）来建模光传播。通常，在OT中，RTE使用其扩散近似进行近似。然而，当对尺寸小于几个散射长度的小目标或含有低散射区域的目标进行成像时，扩散光近似不再有效，需要用辐射传输来建模光传播。$^{1,9}$ 在本文中，RTE的解使用光传输的蒙特卡罗方法进行近似。在蒙特卡罗方法中，光传播通过模拟光子（或光包）并追踪其在介质中的轨迹来近似。$^{12}$ 该方法通常被认为是一种模拟生物组织中光传播的精确方法。有关光传输蒙特卡罗方法的更多信息，参见例如文献$^{13-17}$及其中的参考文献。

在本文中，OT的图像重建问题被公式化为一个最小化问题，其中蒙特卡罗方法作为前向模型。以往在OT中，蒙特卡罗方法已被用于差分成像，例如文献$^{18-24}$。现在考虑绝对成像，图像重建问题按照我们先前提出的随机高斯-牛顿方法来解决。$^{25,26}$ 在该方法中，吸收的雅可比矩阵通过对光包权重的公式进行微分得到。$^{18,25,27}$ 此外，在构建散射的雅可比矩阵时，使用了蒙特卡罗方法的扰动近似。$^{25,26,28,29}$ 正如已注意到的，$^{25,26}$ 光包的数量需要足够大，才能使随机高斯-牛顿算法收敛到正确的最小值。随机噪声的量可以通过调整模拟中使用的光包数量来控制，但仍可能导致计算量庞大的模拟。$^{25,26,30}$

在本文中，我们提出利用机器学习来修正随机高斯-牛顿算法的搜索方向。以往在OT中，机器学习已被用于完全学习算法，$^{31-35}$ 其中从边界数据学习到空间分布光学参数的非线性映射；在图像后处理中，$^{36}$ 其中机器学习用于在初始重建后增强图像质量；以及将机器学习与基于模型的重建方法相结合。$^{37-41}$ 在这些先前的研究中，差分成像被用于重建吸收系数，$^{33,36,37,39,40}$ 或者同时重建吸收和散射系数。$^{34,38}$ 在文献$^{32}$中，估计了吸收和散射的体绝对值。此外，文献$^{31}$、$^{35}$和$^{41}$研究了吸收的绝对成像。在所有这些研究中，RTE的扩散近似被用作光传播模型。在本文中，使用了蒙特卡罗方法。基本上，所提出的方法可以被解释为一种去噪吸收和散射估计的方法；类似地，正如文献$^{42}$中所述，机器学习被用于从随机蒙特卡罗噪声中去除光子密度图像的噪声。

在本文中，使用一种迭代式基于模型的深度学习技术，$^{43-46}$ 有时称为算法展开（algorithm unrolling），来增强以蒙特卡罗方法作为前向模型时的随机高斯-牛顿方法。以往，在OT中，类似的方法论已被用于高斯-牛顿算法中，利用扩散近似和模型不确定性来重建三维绝对吸收和散射分布。$^{43}$ 类似的方法也在文献$^{41}$中被使用，其中使用了一种带有深度学习组件的迭代Levenberg-Marquardt算法来重建吸收系数。我们指出，本文的关键贡献是将深度学习整合到随机高斯-牛顿重建中，以减轻蒙特卡罗噪声，而非开发新的神经网络架构。

本文的其余部分组织如下。在第2节中，介绍了OT图像重建问题和光传输的蒙特卡罗方法，以及随机高斯-牛顿方法。此外，还介绍了所提出的机器学习增强的随机高斯-牛顿方法。在第3节中，展示了数值模拟，其结果随后在第4节中进行讨论和总结。

## 2. 方法

在本文中，研究频域光断层成像，其中强度调制光被引导至成像目标，并测量传输信号的幅度和相位延迟。$^{2}$ 图像重建可以公式化为如下形式的最小化问题，见公式(1)

$$\arg\min_x \left\{ \frac{1}{2} \|L_e (Y_{meas} - Y(x) - \eta_e)\|^2 + \frac{1}{2} \|L_x (x - \eta_x)\|^2 \right\}$$

其中 $Y_{meas} = (Y_{meas,1}, \ldots, Y_{meas,M}) \in \mathbb{R}^M$ 是测量向量，$Y: \mathbb{R}^{2K} \to \mathbb{R}^M$ 是离散化的前向模型，$x = (\mu_a, \mu_s)^T$ 表示未知的光学参数，即吸收 $\mu_a = (\mu_{a,1}, \ldots, \mu_{a,K}) \in \mathbb{R}^K$ 和散射 $\mu_s = (\mu_{s,1}, \ldots, \mu_{s,K}) \in \mathbb{R}^K$。此外，$\eta_e$ 和 $L_e$ 分别是噪声协方差矩阵逆的期望值和Cholesky分解，其中 $\Gamma_e^{-1} = L_e^T L_e$，$\eta_x$ 和 $L_x$ 分别是先验协方差矩阵逆的期望值和Cholesky分解，其中 $\Gamma_x^{-1} = L_x^T L_x$。

在本文中，前向模型 $Y(x)$ 的数值计算基于光传输的蒙特卡罗方法。蒙特卡罗方法模拟辐射传输方程（RTE）的解，其在频域中的形式为，$^{1,47,48}$ 见公式(2)

$$\begin{cases} (i\omega/c + \hat{s} \cdot \nabla + \mu_s(r) + \mu_a(r)) \phi(r, \hat{s}) = \mu_s(r) \int_{S^{d-1}} \Theta(\hat{s} \cdot \hat{s}') \phi(r, \hat{s}') \, d\hat{s}', & r \in \Omega \\ \phi(r, \hat{s}) = \begin{cases} \phi_0(r, \hat{s}), & r \in \bigcup_j \epsilon_j, \, \hat{s} \cdot \hat{n} < 0 \\ 0, & r \in \partial\Omega \setminus \bigcup_j \epsilon_j, \, \hat{s} \cdot \hat{n} < 0 \end{cases} \end{cases}$$

其中 $\Omega \subset \mathbb{R}^d$ 是一个边界为 $\partial\Omega$ 的域，维度 $d=2$ 或 $3$，$i$ 是虚数单位，$\omega = 2\pi f$ 是输入信号的角调制频率，其中 $f$ 是调制频率，$c = c_0/n$ 是介质中的光速，其中 $n$ 是折射率，$c_0$ 是真空中的光速。此外，$\hat{s} \in S^{d-1}$ 是感兴趣方向的单位矢量，$\mu_a(r)$ 是吸收系数，$\mu_s(r)$ 是散射系数，$\phi(r, \hat{s})$ 是辐射度，$\phi_0(r, \hat{s})$ 是位置 $\epsilon_j$ 处的边界光源，$\hat{n}$ 是边界上的向外单位法向量，$\Theta(\hat{s} \cdot \hat{s}')$ 是散射相位函数。在OT中，通常使用Henyey-Greenstein相位函数，$^{49}$ 见公式(3)

$$\Theta(\hat{s} \cdot \hat{s}') = \begin{cases} \frac{1}{2\pi} \frac{1 - g^2}{1 + g^2 - 2g\hat{s} \cdot \hat{s}'}, & d = 2 \\ \frac{1}{4\pi} \frac{1 - g^2}{(1 + g^2 - 2g\hat{s} \cdot \hat{s}')^{3/2}}, & d = 3 \end{cases}$$

其中 $-1 < g < 1$ 是散射各向异性参数。所给出的边界条件假设在边界上，光子只能在光源位置 $\epsilon_j \subset \partial\Omega$ 处沿向内方向传播。在OT中，数据由域边界上的辐射出射度 $Y$ 组成，定义为见公式(4)

$$Y(r) = \int_{S^{d-1}} \phi(r, \hat{s}) (\hat{s} \cdot \hat{n}) \, d\hat{s}$$

### 2.1 光传输的蒙特卡罗方法

在本文中，使用光包蒙特卡罗方法。$^{12}$ 在光包方法中，模拟具有初始权重 $w_0$ 的光包，并追踪它们在介质中经历吸收和散射事件时的路径。$^{12,14,27}$ 假设在传播方向上的小长度 $ds$ 内，光子吸收的概率为 $\mu_a ds$，光子散射的概率为 $\mu_s ds$。光子吸收的概率通过沿轨迹 $s$ 减小光包的权重来考虑。当考虑频域模拟时，权重 $w(s)$ 为见公式(5)

$$w(s) = w_0 \exp\left(-\int_0^s (\mu_a(s') + i\omega/c) \, ds'\right)$$

光子在散射事件前传播的路径长度由散射长度 $l$ 描述，它服从指数概率密度函数，见公式(6)

$$f(l) = \mu_s(l) \exp\left(-\int_0^l \mu_s(l') \, dl'\right)$$

散射角的概率密度函数（在本文中为Henyey-Greenstein相位函数[公式(3)]）定义了散射事件后光包的方向。在蒙特卡罗模拟期间，追踪光包的路径，直到它离开模拟域或被终止。$^{12}$ 然后，使用新的光包重复此模拟，直到获得模拟光传播的足够统计量。

边界元素 $b$ 上的辐射出射度[公式(4)]可以计算为见公式(7)

$$Y_b = \frac{W_b}{PA_b} = \frac{1}{PA_b} \sum_{\text{phot}} w_{\text{phot}}(s) = \frac{1}{PA_b} \sum_{\text{phot}} w_0 \exp\left(-\sum_h (\mu_{a,h} + i\omega/c) l_{h,\text{phot}}\right)$$

其中 $A_b$ 是边界元素 $b$ 的宽度（$d=2$）或面积（$d=3$），$P$ 是模拟的光包数量，$W_b$ 是通过边界元素 $b$ 离开模拟域的光包的权重之和，"phot"表示通过边界元素 $b$ 逃逸的光包，$l_{h,\text{phot}}$ 是在分段常数离散化光学参数的离散化元素 $h$ 中离散化的光子路径长度。$^{18,25,27}$

### 2.2 随机高斯-牛顿方法

由于蒙特卡罗方法的随机性，最小化问题[公式(1)]也是随机的，可以使用随机优化方法$^{30,50}$来处理，例如随机高斯-牛顿（SGN）方法。$^{25,26,30}$ 在SGN方法中，光学参数 $x_i$ 通过公式(8)更新：

$$x_{i+1} = x_i + \Delta x_i$$

其中 $\Delta x_i = \alpha_i \delta(x_i)$ 是高斯-牛顿步，$\delta(x_i)$ 是用 $P$ 个光包计算的高斯-牛顿搜索方向，$\alpha_i$ 是步长参数，在本文中设为 $\alpha = 1$。SGN搜索方向的形式为见公式(9)

$$\delta(x_i) = (J^T(x_i)\Gamma_e^{-1}J(x_i) + \Gamma_x^{-1})^{-1} (J^T(x_i)\Gamma_e^{-1}(Y_{meas} - Y(x_i) - \eta_e) - \Gamma_x^{-1}(x_i - \eta_x))$$

其中 $Y(x_i)$ 是随机前向模型，$J(x_i) = (J_{\mu_a}(x_i), J_{\mu_s}(x_i))$ 是在更新 $x_i$ 处用 $P$ 个光包计算的前向模型的雅可比矩阵。

要构建雅可比矩阵，需要辐射出射度的吸收和散射导数。辐射出射度的吸收导数可以直接从公式(4)和(7)计算得到，$^{18,25}$ 即，在离散化元素 $k$ 中的边界元素 $b$ 上，辐射出射度 $Y_b$ 对吸收 $\mu_{a,k}$ 的导数为见公式(10)

$$\frac{\partial Y_b}{\partial \mu_{a,k}} = \frac{1}{PA_b} \sum_{\text{phot}} (-l_{k,\text{phot}}) w_{\text{phot}}(s) = \frac{1}{PA_b} \sum_{\text{phot}} (-l_{k,\text{phot}}) w_0 \exp\left(-\sum_h (\mu_{a,h} + i\omega/c) l_{h,\text{phot}}\right)$$

辐射出射度的散射导数可以利用扰动蒙特卡罗方法计算，$^{27,28,51}$ 该方法通过重用未扰动模拟的轨迹来评估光学参数扰动的影响。$^{19,29}$ 扰动散射系数 $\tilde{\mu}_s$ 的扰动权重 $\tilde{w}$ 可以通过公式(11)获得：

$$\tilde{w} = w \left(\frac{\tilde{\mu}_s}{\mu_s}\right)^{n_s} \exp\left(-(\tilde{\mu}_s - \mu_s) L_{\text{tot}}\right)$$

其中 $w$ 是未扰动的权重，$n_s$ 是扰动区域内的散射事件数，$L_{\text{tot}}$ 是光包在扰动区域内传播的总距离。$^{19,29}$ 使用扰动近似，在离散化元素 $k$ 中的边界元素 $b$ 上，辐射出射度 $Y_b$ 对散射系数 $\mu_{s,k}$ 的导数为$^{25}$ 见公式(12)

$$\frac{\partial Y_b}{\partial \mu_{s,k}} = \frac{1}{PA_b} \sum_{\text{phot}} \left(\frac{n_s}{\mu_{s,k}} - l_{k,\text{phot}}\right) w_{\text{phot}}(s) = \frac{1}{PA_b} \sum_{\text{phot}} \left(\frac{n_s}{\mu_{s,k}} - l_{k,\text{phot}}\right) w_0 \exp\left(-\sum_h (\mu_{a,h} + i\omega/c) l_{h,\text{phot}}\right)$$

### 2.3 深度随机高斯-牛顿方法

在本文中，我们提出了一种利用机器学习方法$^{43-46}$ 改进随机高斯-牛顿更新公式(8)的方法。在所提出的方法中，随机高斯-牛顿更新[公式(8)]被替换为一个学习到的更新函数，见公式(13)

$$x_{i+1} = G_{\theta_i}(x_i, \Delta x_i)$$

其中函数 $G_{\theta_i}$ 对应于具有选定架构和每次迭代 $i$ 上不同学习参数的卷积神经网络（CNN）。因此，我们可以写出一个深度随机高斯-牛顿（DSGN）算法，其中光学参数的（含噪声的）更新 $\Delta x_i$ 在每次迭代中使用CNN进行修正。$^{43,44}$

DSGN的训练算法如算法1所示。为了考虑OT图像重建问题的非线性，采用迭代展开方法。$^{44}$ 这意味着在每个迭代中，CNN逐层应用于光学参数 $x_i = (\mu_{a,i}, \mu_{s,i})^T$ 和随机高斯-牛顿更新 $\Delta x_i = (\Delta\mu_{a,i}, \Delta\mu_{s,i})^T$。在训练阶段，在每个迭代中，提供真实的吸收和散射分布、当前的高斯-牛顿估计 $x_i$ 以及用 $P$ 个光包计算的当前高斯-牛顿更新 $\Delta(x_i)$ 作为流水线的输入。

为了提高数值稳定性，吸收和散射都通过变量变换重新参数化，使其值处于相同的范围内。52 训练数据由真实吸收和散射图像 \(\{\mu_{a,\text{true}}, \mu_{s,\text{true}}\}_j, j=1,\ldots,N_{\text{samp}}\) 组成，其中 \(N_{\text{samp}}\) 是样本数量。这些样本用于通过蒙特卡洛方法模拟OT数据和SGN更新，然后用加性噪声进行污染。在本工作中，网络参数 \(\theta_i\)（\(i=1,\ldots,I\)，其中 \(I\) 是迭代次数）通过最小化所有样本 \(j\) 的 \(L2\) 损失函数依次训练，由公式 (14) 给出：
\[
\min_{\theta_i} \sum_j \|\mu_{j,a,i+1} - \mu_{j,a,\text{true}}\|^2 + \|\mu_{j,s,i+1} - \mu_{j,s,\text{true}}\|^2,
\]
其中更新值 \(\mu_{j,a,i+1}\) 和 \(\mu_{j,s,i+1}\) 由公式 (13) 给出。本工作使用的CNN架构基于之前的研究，44 如图1所示。该架构的变体也曾在其他断层扫描应用中使用过，例如参考文献43、53和54。该架构包括特征提取部分和重组部分。在网络中，吸收和散射参数在不同的通道中处理。首先，参数图通过一个5×5卷积核的卷积层扩展到20个通道，然后通过维度 \(d=2\) 的卷积层扩展到40个通道。使用“泄漏”整流线性单元（LReLU）作为激活函数，以允许输入参数为负值。在重组部分，这两个流水线的输出相加，首先使用LReLU作为激活函数减少到20个通道，然后使用标量乘法减少到1个通道。最后，结果与当前迭代中的参数值相加，并使用整流线性单元（ReLU）作为激活函数投影到正数。由于本工作的主要贡献不是提出特定的CNN架构，因此架构保持简单。

算法1 训练深度随机高斯-牛顿

1: 从先验模型中抽取集合 \(\{x_{\text{true}}\}_j, j=1,\ldots,N_{\text{samp}}\)。

2: 使用蒙特卡洛方法模拟数据（第2.1节），并用加性噪声进行污染。

3: 为所有 \(j\) 将 \(x_{j,1}\) 设置为选定的初始值，例如先验的均值。

4: 设置 \(i \leftarrow 1\)。

5: 当 \(i < I\) 时：

6: 使用 \(P\) 个光子包和公式 (9) 为所有 \(j\) 计算 \((\delta x_i)_j\)，并设置 \(\alpha_i\)。

7: 函数 TRAIN(\((x_{j,i}, (\Delta x_i)_j, x_{j,\text{true}})\), 对所有 \(j\))

8: 通过最小化公式 (14) 训练参数 \(\theta_i\)。

9: 结束函数 返回 \(\theta_i\)

10: 对所有 \(j\) 设置 \(x_{j,i+1} \leftarrow G_{\theta_i}(x_{j,i}, (\Delta x_i)_j)\)。

11: 设置 \(i \leftarrow i+1\)

12: 结束当

图1 下载完整尺寸图像

一个表示一次深度随机高斯-牛顿迭代的卷积神经网络 \(G_{\theta_i}\) 的示意图。带有5×5卷积核的卷积层，后面跟着LReLU，用红色箭头标记。带有5×5卷积核的卷积层，后面跟着标量乘法 \(\lambda_j\)，用蓝色箭头标记。得到的通道数标记在方块内。通过跳跃连接（灰色箭头），残差更新通过最后一个ReLU（用绿色箭头标记）投影到正数。

DSGN算法可以通过将学习到的网络 \(G_{\theta_i}\) 和训练好的参数集 \(\theta_i\) 应用于每次迭代 \(i\) 的光学参数和高斯-牛顿更新来进行评估。评估对应于算法1，从设置 \(x_{j,1}\) 开始，并跳过“TRAIN”函数。

3. 仿真

二维数值仿真用于评估所提出的方法。最小化算法和仿真在MATLAB（R2024a, MathWorks Inc., Natick, Massachusetts, United States）中实现。在深度随机高斯-牛顿的实现和训练中，使用了Python库TensorFlow（版本2.18）。蒙特卡洛仿真在C++中实现。仿真在AMD Ryzen Threadripper PRO 7965WX 24核CPU计算服务器上进行，具有48个逻辑线程和250 GiB RAM。DSGN的训练在NVIDIA RTX 6000 Ada Generation GPU上进行，具有49,140 MiB内存。

3.1. 仿真几何和离散化

仿真域是一个5 mm × 5 mm的正方形。在仿真设置中，考虑了20个光源（每边5个）和100个探测器（每边25个），宽度为0.2 mm。光源放置在每边上，两个光源之间间隔0.5 mm，居中放置，使得最外侧光源到正方形角落的距离为1 mm。探测器并排放置，覆盖仿真域的边长。一次使用一个光源，并在正方形相邻和相对的边上记录数据。

对于数据仿真，吸收和散射以及光子密度在具有5000个元素的三角形网格中分段常数表示，用于训练和验证。对于测试数据，对于分布内目标（图2第一列所示）使用具有5000个元素的三角形网格。对于分布外目标（图4中稍后展示），圆形目标使用9006个元素，血管目标使用20,000个元素。在图像重建和训练中，吸收和散射在50×50像素网格的分段常数离散化中表示，即2500个像素，光子密度在由5000个元素组成的三角形网格中仿真。

3.2. 数据仿真

目标的光学参数被选择以模拟生物组织的光学特性。55,56 在所有仿真中，折射率为 \(n=1\)，散射各向异性参数为 \(g=0.9\)，强度调制光的频率为 \(f=100\ \text{MHz}\)。此外，在蒙特卡洛实现中，光源形状被建模为角余弦和空间均匀。

训练数据、验证目标和分布内（训练）测试目标通过分割从平方指数先验57中绘制的图像进行K-means聚类来仿真。然后，背景和包含物的吸收和散射值从均匀分布中抽取。这些吸收和散射图像的示例如图2所示。总共仿真了800个训练目标和200个验证目标。此外，考虑了5个分布内测试目标和2个分布外目标（图4），即圆形目标和血管目标。血管目标是从数值乳房体模58修改而来。

光子密度和出口度使用蒙特卡洛方法仿真，如第2.1节所述。数据使用每个光源 \(10^8\) 个光子包进行仿真。作为数据类型，考虑了复出口度的振幅和相位[公式 (7)]。模拟的无噪声数据被污染了零均值、标准差对应于无噪声数据相对幅度0.5%的加性高斯分布随机噪声，以模拟现实的高噪声水平。

3.3. 训练

网络按照第2.3节描述的方法使用训练集的吸收和散射图像以及相应的蒙特卡洛仿真数据进行训练。在训练和验证期间，DSGN的更新方向[公式 (9)]使用 \(10^6\) 个光子包计算。作为初始猜测 \(\mu_{j,1}\)，使用了先验的缩放均值。参数 \(\theta_i\) 通过使用TensorFlow实现的Adam优化器59最小化公式 (14)进行训练，批大小为8，25个epochs，学习率为 \(5 \times 10^{-4}\)。迭代次数 \(I\) 设置为20，根据我们的仿真，这提供了该仿真设置下随机高斯-牛顿迭代的收敛。本工作中所有DSGN重建研究都使用相同的训练网络。

3.4. 重建

使用DSGN近似最小化问题[公式 (1)]的解，由损失函数[公式 (14)]定义（第2.3节）。结果与使用SGN方法（第2.2节）求解的解进行比较。在图像重建中，噪声被建模为零均值的高斯噪声，标准差 \(\sigma_e\) 指定为模拟噪声数据相对幅度的0.5%。此外，Ornstein–Uhlenbeck模型被用作先验模型。57 Ornstein–Uhlenbeck先验属于Matérn协方差函数类，根据作者的先前研究，它是不同成像目标的有效和灵活先验模型。Ornstein–Uhlenbeck先验的协方差矩阵定义为：
\[
\Gamma_x(i,j) = \sigma_x^2 \exp(-\|r_i - r_j\|/\tau),
\]
其中 \(\sigma_x\) 是标准差，\(\tau\) 表示控制空间平滑度的特征长度尺度参数，\(r_i\) 和 \(r_j\) 是离散化点。在所有仿真中，吸收的先验均值为 \(\eta_{\mu_a} = 0.02\ \text{mm}^{-1}\)，散射为 \(\eta_{\mu_s} = 5\ \text{mm}^{-1}\)。吸收的标准差为 \(\sigma_{\mu_a} = 0.0267\ \text{mm}^{-1}\)，散射为标准差 \(\sigma_{\mu_s} = 1.67\ \text{mm}^{-1}\)。标准差值被设定使得吸收和散射的最大可能对比度对应于距背景的三个标准差。特征长度尺度为 \(\tau = 0.5\ \text{mm}\)。

为了提高吸收和散射同时重建的数值稳定性，光学参数通过除以先验均值重新参数化到相同范围内。52 缩放后的先验均值用作所有仿真中高斯-牛顿的初始值。在DSGN方法中，每次迭代中，每个光源仿真 \(10^6\) 个光子包。在SGN方法中，重建评估使用每个光源每次迭代 \(10^6\) 和 \(10^9\) 个光子包。在所有重建中，迭代次数为20，这是为了使所有重建结果具有可比性而选择的。需要注意的是，在实际应用该方法时，会选择收敛准则，在最小化问题收敛后立即结束迭代，例如参考文献60。

3.5. 结果

五个分布内测试目标的绝对吸收和散射重建图像如图2所示。可以看出，DSGN方法提供的重建与使用 \(10^9\) 个光子包计算的SGN参考重建高度相似。使用DSGN和SGN都可以区分包含物的位置。此外，散射图像的对比度与真实目标相比是相似的。在吸收重建中，与使用大量光子包的SGN重建相比，DSGN方法的对比度略好。将DSGN重建与使用相同光子包数量（\(10^6\)）的SGN相比，重建结果噪声更小，使用DSGN时包含物区分得更好，尤其是在比较散射重建时。在所有重建图像中，DSGN重建的背景比SGN重建更平滑。

图2 下载完整尺寸图像

与训练集对应的测试目标的吸收系数 \(\mu_a\ \text{(mm}^{-1})\)（第1至4列）和散射系数 \(\mu_s\ \text{(mm}^{-1})\)（第5至8列）重建图像。真实目标在第1和第5列，深度随机高斯-牛顿重建（DSGN）在第2和第6列，使用 \(10^6\) 个光子包的随机高斯-牛顿重建（SGN）在第3和第7列，使用 \(10^9\) 个光子包的随机高斯-牛顿重建（SGNref）在第4和第8列。

除了视觉比较外，通过计算吸收和散射的相对误差来比较结果：
\[
E_{\mu_a} = \frac{\|\mu_a - \mu_{a,\text{true}}\|}{\|\mu_{a,\text{true}}\|} \cdot 100\%, \quad E_{\mu_s} = \frac{\|\mu_s - \mu_{s,\text{true}}\|}{\|\mu_{s,\text{true}}\|} \cdot 100\%,
\]
其中范数是欧几里得范数，\(\mu_a\) 和 \(\mu_s\) 是重建值，\(\mu_{a,\text{true}}\) 和 \(\mu_{s,\text{true}}\) 是插值到重建离散化的真实吸收和散射系数。图2第一行目标的吸收和散射重建相对于迭代的相对误差如图3第一列所示。根据相对误差，DSGN收敛到与使用大量光子包的SGN相同的水平。在实践中，使用大量光子包的SGN在六次迭代内收敛，而DSGN在16次迭代内收敛，当应用连续三次迭代估计之间的相对差异小于5%的收敛准则时。60 使用 \(10^6\) 个光子包的SGN受到更高的随机噪声影响，并且未能使用与其他迭代相同的准则收敛。还请注意，在使用 \(10^6\) 个光子包的SGN重建中，第一次迭代的吸收相对误差高于初始猜测。我们认为这是由于吸收和散射估计之间的串扰效应造成的，而DSGN可以纠正这种效应。43

图3 下载完整尺寸图像

吸收 \(E_{\mu_a}(\%)\)（a）和散射 \(E_{\mu_s}(\%)\)（b）相对于迭代的相对误差。

从左到右各列：图2第一行中目标的相对误差（第一列）、圆形目标的相对误差（第二列）和血管目标的相对误差（第三列），分别使用深度随机高斯-牛顿方法（DSGN，红色实线）以及随机高斯-牛顿方法（SGN，使用10^6个光子包，蓝色虚线）和10^9个光子包（SGNref，蓝色点划线）进行重构。图2(a)中目标的前向（数据）模拟、训练和重构的计算时间分别使用DSGN和SGN方法进行了记录。使用10^8个光子包的数据模拟耗时321秒。前向和雅可比矩阵评估、从MATLAB运行Python训练/评估脚本、神经网络推理、一次高斯-牛顿迭代、训练20次迭代的总计算时间以及重构算法直至收敛的总时间列于表1中。$t_{\text{python}}$对应于使用MATLAB内置函数system执行Python训练（或评估）脚本所需的总时间。这包括启动Python、导入所需的Python库（如TensorFlow）以及网络计算时间$t_{\text{network}}$所需的时间。如图所示，当训练DSGN时，最耗时的部分是针对800个训练目标和200个验证目标的前向和雅可比矩阵评估，一次DSGN迭代的总耗时约为9600秒。此外，运行Python训练脚本大约需要75秒，而训练本身每次迭代大约需要15秒。DSGN一次迭代的计算时间比使用10^9个光子包的SGN迭代快约160倍。由于DSGN收敛较慢（16次迭代），而SGN参考重构收敛较快（6次迭代），DSGN的总计算时间比参考SGN快约60倍。另一方面，当将DSGN的迭代时间与使用相同光子包数量（10^6）的SGN进行比较时，DSGN大约比SGN慢三倍。

表1 使用P个光子包，针对图2第一行中目标的训练、深度随机高斯-牛顿（DSGN）和随机高斯-牛顿（SGN）算法的计算时间。前向和雅可比矩阵评估时间（$t_{Y,J}$）、从MATLAB运行Python训练/评估脚本的时间（$t_{\text{python}}$）、网络计算时间（$t_{\text{network}}$）、一次高斯-牛顿迭代时间（$t_{\text{iter}}$，20次运行的平均值），以及训练20次迭代和重构算法直至收敛的总计算时间$t_{\text{tot}}$。使用10^6个光子包的SGN算法未收敛到期望标准。

| | $P$ | $t_{Y,J}$(s) | $t_{\text{python}}$(s) | $t_{\text{network}}$(s) | $t_{\text{iter}}$(s) | $t_{\text{tot}}$(s) |
|---|---|---|---|---|---|---|
| 训练 | | | | | 10^6 | 975 | 15 | 9,613 | 192,260 |
| DSGN | | | | | 10^6 | 9 | 18 | 128 | 448 |
| SGN | | | | | 10^6 | 9 | — | — | 10 | — |
| $SGN_{ref}$ | | | | | 10^9 | 4,546 | — | — | 4,547 | 27,282 |

由于SGN方法的随机性，当使用不同的光子包实现集进行图像重构时，结果可能会有所不同。为了提供方法的性能统计信息，使用10^6个光子包计算的DSGN和SGN重构重复了100次。然后计算了吸收$E_{\mu_a}$和散射$E_{\mu_s}$的相对误差均值和标准差，以及吸收$SSIM_{\mu_a}$和散射$SSIM_{\mu_s}$的结构相似性指数的均值和标准差。使用10^9个光子包获得的SGN重构的相对误差被视为参考值。参考值以及相对误差和结构相似性指数的均值和标准差列于表2中。此外，为了在一组目标上研究该方法的性能，针对来自不同分布内目标的100个吸收和散射重构数据，计算了相对误差和结构相似性指数的均值和标准差，同样列于表2中。如图所示，与使用10^6个光子包评估的SGN方法相比，DSGN方法提供了更低的相对误差和更高的SSIM值；与使用10^9个光子包的SGN相比，则具有相似的相对误差和SSIM值。此外，与使用10^6个光子包的SGN相比，DSGN的相对误差标准差更小。在比较结构相似性指数的标准差时，可以看到DSGN的吸收重构标准差更小，而DSGN的散射重构标准差与使用10^6个光子包的SGN相同。需要注意的是，SSIM已被证明在有噪图像中表现不佳。此外，SSIM对对比度高度敏感，这在边缘错位的情况下尤为明显。

表2 使用深度随机高斯-牛顿（DSGN）和随机高斯-牛顿（SGN）方法（使用10^6个光子包）计算的吸收和散射估计值的相对误差$E_{\mu_a}$（%）和$E_{\mu_s}$（%）以及结构相似性指数$SSIM_{\mu_a}$和$SSIM_{\mu_s}$的均值和标准差，以及其参考值（$SGN_{ref}$，使用10^9个光子包计算的随机高斯-牛顿方法）。第1至5行：从分布内测试目标1至5（图2）的数据中100次重复重构计算得到的统计结果。第6行：一组100个分布内测试目标的重构统计结果。

| | DSGN | SGN | $SGN_{ref}$ | DSGN | SGN | $SGN_{ref}$ |
|---|---|---|---|---|---|---|
| $E_{\mu_a}$（%） | | | | $SSIM_{\mu_a}$ | | |
| 1 | 37.6±2.3 | 54.2±3.7 | 39.2 | 0.90±0.01 | 0.75±0.03 | 0.86 |
| 2 | 34.8±1.9 | 48.6±3.1 | 33.8 | 0.88±0.01 | 0.78±0.02 | 0.88 |
| 3 | 35.8±2.3 | 54.1±3.6 | 36.2 | 0.86±0.01 | 0.74±0.03 | 0.86 |
| 4 | 39.7±2.7 | 51.2±5.3 | 32.7 | 0.89±0.01 | 0.80±0.03 | 0.92 |
| 5 | 35.4±1.6 | 51.1±3.7 | 35.9 | 0.90±0.01 | 0.74±0.02 | 0.85 |
| 集合 | 36.9±7.2 | 52.8±8.5 | | 0.89±0.04 | 0.76±0.04 | |
| $E_{\mu_s}$（%） | | | | $SSIM_{\mu_s}$ | | |
| 1 | 8.3±0.2 | 19.4±0.6 | 8.2 | 0.42±0.01 | 0.20±0.01 | 0.40 |
| 2 | 4.9±0.2 | 17.5±0.5 | 5.8 | 0.30±0.01 | 0.08±0.01 | 0.24 |
| 3 | 12.7±0.4 | 20.9±0.7 | 11.7 | 0.55±0.01 | 0.38±0.01 | 0.55 |
| 4 | 5.4±0.2 | 17.8±0.5 | 5.7 | 0.49±0.01 | 0.14±0.01 | 0.48 |
| 5 | 8.9±0.4 | 19.0±0.6 | 9.9 | 0.42±0.01 | 0.22±0.01 | 0.33 |
| 集合 | 8.9±1.9 | 18.9±1.6 | | 0.36±0.06 | 0.19±0.05 | |

### 3.5.1. 对分布外目标的泛化能力

使用两个分布外目标——圆形目标和血管目标——测试了DSGN方法的泛化能力。圆形目标以及使用10^6和10^9个光子包的DSGN和SGN方法的重构结果如图4第一行所示，血管目标及其重构结果如图4第二行所示。可以看到，DSGN对圆形目标的重构结果与使用10^9个光子包的SGN具有相似的质量。另一方面，与使用相同光子包数量的SGN相比，DSGN重构中包层的形状更圆，定位更好，尤其是在吸收重构中。此外，DSGN改善了吸收包层的对比度。DSGN的散射重构明显比使用相同光子包数量的SGN噪声更少。在DSGN对血管目标的重构中，吸收包层的血管被分割成圆形形状。此外，与使用10^9个光子包的SGN估计值相比，散射图像中的血管更宽且对比度更低。当将DSGN重构与使用相同光子包数量的SGN进行比较时，可以看到在吸收和散射重构中包层的形状都能更好地被区分。在应用收敛标准（连续三次迭代的估计值之间相对差异小于5%）时，DSGN在圆形数据上经过13次迭代收敛，在血管数据上经过16次迭代收敛。使用相同光子包数量的SGN重构未能收敛。此外，使用10^9个光子包的SGN在圆形数据上经过9次迭代收敛，在血管数据上经过8次迭代收敛。

**图4** 圆形和血管目标（分别显示在第1行和第2行）的重构吸收系数$\mu_a$（mm$^{-1}$）（第1至4列）和散射系数$\mu_s$（mm$^{-1}$）（第5至8列）。第1列和第5列显示真实目标，第2列和第6列显示深度随机高斯-牛顿（DSGN）重构结果，第3列和第7列显示使用10^6个光子包的随机高斯-牛顿（SGN）重构结果，第4列和第8列显示使用10^9个光子包的随机高斯-牛顿（SGNref）重构结果。

吸收和散射估计值相对于迭代次数的相对误差在图3的第二列（圆形目标）和图3的第三列（血管目标）中给出。此外，从100次重复重构中计算的相对误差和SSIM的均值及标准差，以及使用10^9个光子包的SGN获得的参考值，列于表3中。DSGN在圆形目标上的相对误差和结构相似性指数与使用10^9个光子包的SGN处于同一水平。在血管目标的重构中，与使用10^9个光子包的SGN相比，DSGN的相对误差略高，散射的结构相似性指数略低。DSGN吸收重构的结构相似性指数与参考SGN重构值处于同一水平。当将DSGN与使用相同光子包数量的SGN进行比较时，吸收和散射重构中的相对误差均值都更低，结构相似性指数均值都更高。此外，DSGN估计值的相对误差和结构相似性指数的标准差更小，唯一例外是血管目标的散射SSIM，其标准差与使用10^6个光子包的SGN相同。

表3 使用深度随机高斯-牛顿（DSGN）和随机高斯-牛顿（SGN）方法（使用10^6个光子包）从100次重复重构中计算的吸收和散射估计值的相对误差$E_{\mu_a}$（%）和$E_{\mu_s}$（%）、结构相似性指数$SSIM_{\mu_a}$和$SSIM_{\mu_s}$、位置误差$E_{\mu_{a},pos}$和$E_{\mu_{s},pos}$以及空间重叠度$Ov_{\mu_a}$和$Ov_{\mu_s}$的均值和标准差，以及使用10^9个光子包的随机高斯-牛顿方法（$SGN_{ref}$）计算的圆形和血管目标的参考值。

| | 圆形 | 血管 | 圆形 | 血管 | 圆形 | 圆形 |
|---|---|---|---|---|---|---|
| $E_{\mu_a}$（%） | $SSIM_{\mu_a}$ | $E_{\mu_{a},pos}$（mm） | $Ov_{\mu_a}$ | | | |
| DSGN | 38.5±2.0 | 50.0±1.7 | 0.89±0.01 | 0.85±0.01 | 0.031±0.011 | 0.55±0.07 |
| SGN | 53.3±4.7 | 65.0±4.4 | 0.74±0.03 | 0.67±0.04 | 0.054±0.027 | 0.27±0.17 |
| $SGN_{ref}$ | 38.1 | 47.8 | 0.87 | 0.82 | 0.032 | 0.46 |
| $E_{\mu_s}$（%） | $SSIM_{\mu_s}$ | $E_{\mu_{s},pos}$（mm） | $Ov_{\mu_s}$ | | | |
| DSGN | 9.5±0.4 | 15.7±0.3 | 0.38±0.01 | 0.38±0.01 | 0.005±0.002 | 0.81±0.05 |
| SGN | 20.1±0.7 | 21.3±0.7 | 0.17±0.01 | 0.28±0.01 | 0.007±0.004 | 0.26±0.09 |
| $SGN_{ref}$ | 10.6 | 13.1 | 0.33 | 0.45 | 0.005 | 0.49 |

对于圆形目标，计算了另外两个误差度量——位置误差和空间重叠度——用于评估定位误差。吸收$E_{\mu_{a},pos}$和散射$E_{\mu_{s},pos}$的位置误差通过计算重构图像的中心质量与模拟目标图像的中心质量之间的欧几里得距离得到。吸收$Ov_{\mu_a}$和散射$Ov_{\mu_s}$的空间重叠度使用Dice系数量化，定义为重构图像和目标图像中最突出包层之间重叠体素的比率。此外，这些度量值也是从100次重复重构中计算得出的。这些度量值列于表3中，可以观察到DSGN和参考SGN重构显示出相当的位置误差，而DSGN相比SGN重构提供了改进的空间重叠度（Dice系数）。

### 3.5.2. 对3D数据的泛化能力

此外，研究了2D训练模型对3D数据的泛化能力。为此，在5×5×5 mm3域中模拟了蒙特卡罗数据。光源和探测器放置在立方体边界中间层，遵循2D成像几何的光源和探测器设置，即使用20个光源（立方体每边5个）和100个探测器（立方体每边25个）。该域使用35,639个四面体元素离散化，如图5所示，同时也展示了立方体一侧上光源和探测器的位置。数据使用一个由两个吸收圆柱和四个散射圆柱组成的目标进行模拟，圆柱高度为5 mm。真实目标参数在中间层处的横截面也在图5中进行了可视化。数据模拟中每个光源使用了10^9个光子包。然后，数据经过预校准以匹配2D和3D模型之间的光源强度，使用恒定光学属性的2D和3D模拟之间的数据差异进行校准。

使用DSGN方法计算的重构吸收和散射系数以及与SGN方法的比较结果如图5所示。此外，在将真实3D吸收和散射参数插值到2D重构网格后，计算了吸收和散射估计值的相对误差（公式(16)）和SSIM，列于表4中。如图所示，当将DSGN方法（2D训练）应用于3D数据时，可以获得合理的重构结果。

**图5** (a) 3D模拟域，展示立方体一侧上的光源（红色）和探测器（蓝色）位置。(b) 从3D模拟数据中重构的吸收$\mu_a$（mm$^{-1}$）（第一行）和散射$\mu_s$（mm$^{-1}$）系数（第二行）的2D重构结果。从左到右：目标中间层处的真实目标参数横截面（第一列），使用深度随机高斯-牛顿（DSGN，第二列）和随机高斯-牛顿（SGN，第三列）方法计算的重构结果。

表4 使用深度随机高斯-牛顿（DSGN）和随机高斯-牛顿（SGN）方法计算的吸收和散射估计值的相对误差$E_{\mu_a}$（%）和$E_{\mu_s}$（%）以及结构相似性指数$SSIM_{\mu_a}$和$SSIM_{\mu_s}$。

泛化研究结果：在3D中模拟的数据，使用不同的加性噪声水平（1%和5%）进行数据模拟，以及使用不同数量的光子包（10?和10?）进行数据模拟。

| 数据类型 | Eμa (%) | SSIMμa | Eμs (%) | SSIMμs |
|---|---|---|---|---|
| | DSGN | SGN | DSGN | SGN | DSGN | SGN | DSGN | SGN |
| 3D | 79.6 | 164.8 | 0.65 | 0.23 | 21.4 | 41.9 | 0.30 | 0.13 |
| 1%噪声 | 40.7 | 51.4 | 0.89 | 0.80 | 8.4 | 15.1 | 0.43 | 0.25 |
| 5%噪声 | 53.1 | 56.2 | 0.83 | 0.81 | 10.9 | 14.1 | 0.33 | 0.25 |
| 10?光子包 | 66.0 | 209.3 | 0.77 | ?0.01 | 16.0 | 57.1 | 0.15 | 0.05 |
| 10?光子包 | 48.7 | 81.9 | 0.86 | 0.58 | 9.2 | 27.4 | 0.37 | 0.15 |

3.5.3. 对不同噪声水平的泛化

DSGN方法的泛化能力也在不同加性噪声水平和不同光子包数量的模拟数据上进行了测试。使用10?光子包模拟并添加1%和5%加性噪声的数据，以及使用10?和10?光子包模拟并添加0.5%加性噪声的数据均被研究。注意，在重建算法的噪声模型中[噪声协方差矩阵Γ?，公式(1)和(8)]，噪声的标准差是根据模拟的加性噪声水平进行建模的。

使用DSGN方法计算得到的重建吸收系数和散射系数与SGN方法的对比结果如图6所示。此外，吸收和散射估计的相对误差（公式16）和结构相似性指标（SSIM）见表4。可以看出，DSGN方法在较高加性噪声和较高蒙特卡洛噪声的情况下均能改善SGN的重建结果。然而，当噪声统计特性与训练模型和/或图像重建算法的噪声模型的统计特性存在显著差异时，其泛化能力较差，例如使用10?光子包模拟数据的情况。

**图6**

**下载 全尺寸图像**

图2第一行中目标的吸收系数μa (mm?1)（a）和散射系数μs (mm?1)（b）的重建结果。第1列显示真实目标。第2至3列显示1%噪声模拟数据下的深度随机高斯-牛顿（DSGN）和随机高斯-牛顿（SGN）重建结果，第4至5列显示5%噪声模拟数据下的重建结果。第6至7列显示使用10?光子包模拟数据下的DSGN和SGN重建结果，第8至9列显示使用10?光子包模拟数据下的重建结果。

3.5.4. 与深度学习的后处理方法比较

我们还将该方法与文献63中提出的深度学习后处理方法进行了比较。对于后处理方法，采用了与文献63中相同的网络架构，对使用10?光子包获得的SGN吸收和散射估计进行后处理。用于吸收和散射的后处理网络使用与训练DSGN相同的800个目标样本（吸收和散射图像对）及其对应的含噪声SGN估计进行训练。三个目标（一个分布内和两个分布外）的后处理吸收和散射重建示例如图7所示。这些估计的相对误差分别为Eμa=78.9, 65.4, 78.8%和Eμs=89.8, 88.5, 86.6%，结构相似性指标分别为SSIMμa=0.69, 0.73, 0.70和SSIMμs=0.22, 0.18, 0.26。将图7中的重建结果与图2和图4中的DSGN重建结果以及表2和表3中的误差进行比较，可以看出后处理估计的误差更高，尤其是在分布外情况下，散射估计显示出非现实的高对比度。

**图7**

**下载 全尺寸图像**

图2第一行中目标的吸收系数μa (mm?1)（第1至3列）和散射系数μs (mm?1)（第4至6列）的重建结果，以及两个分布外案例（圆形和管状）。第1行显示真实目标，第2行显示后处理（POST）重建结果。

4. 讨论与结论

在本工作中，提出了一种基于模型的迭代深度学习方法来重建光声层析成像（OT）中吸收和散射的图像，以应对由于使用蒙特卡洛方法作为正向模型进行光传输计算而产生的随机噪声。在该方法中，高斯-牛顿算法的搜索方向在每次迭代中被一个学习到的函数所替代。该学习到的函数基于卷积神经网络（CNN），使用一组吸收和散射图像以及使用较少光子包计算得到的高斯-牛顿搜索方向进行训练。该方法的性能通过模拟进行了研究，并与传统的随机高斯-牛顿（SGN）算法和深度学习的后处理方法进行了比较。

结果表明，所提出的DSGN算法可用于补偿由于使用较少光子包而产生的随机噪声。如图2所示，与使用相同少光子包（10?）的SGN重建结果相比，DSGN提供了更平滑的图像，具有更好的对比度和包埋物定位能力。此外，DSGN提供了质量更好的重建结果、更低的平均相对误差值以及更高的平均结构相似性指标值（表2）。此外，DSGN的重建结果与使用大量光子包（10?）的SGN重建结果相当，具有相似的相对误差和结构相似性指标水平。然而，需要注意的是，所有吸收和散射估计的相对误差都相对较高。因此，虽然OT可以用于重建较低对比度的目标，且DSGN明显改善了重建质量，但仍需进一步研究以提高该方法的定量准确性。

DSGN方法的泛化能力通过使用两个训练集外的测试目标（图3和图4）、3D模拟数据（图5）以及不同噪声水平的模拟数据（图6）进行了评估。结果表明，该方法能够良好地泛化到分布外目标，可用于重建各种吸收和散射分布。使用经过校准的3D数据的2D DSGN重建显示了学习模型在源和探测器位于2D平面上且包埋物穿过该平面时向3D的良好泛化水平。尽管如此，这种2D近似模型并不适用于大多数OT的实际应用，因此，未来的工作包括将该方法扩展到完整的3D。这在计算上更加昂贵，因为已经证明，例如，蒙特卡洛模拟的准确性取决于离散化。此外，3D方法可以通过机器学习方法用于模型降阶来增强，例如文献43。该方法对不同噪声水平数据的泛化能力也得到了研究。发现该方法在不同加性噪声和随机噪声水平下泛化良好，但当噪声统计特性与训练模型和/或图像重建算法的噪声模型的统计特性存在显著差异时，泛化能力较差。总体而言，通过仔细选择训练数据的规模和风格，以及进一步优化训练参数（如批大小、迭代次数和学习率），可以提高泛化能力。

在所有重建中，散射重建中的噪声明显高于吸收重建中的噪声。我们认为这种较高的噪声可能与在散射的雅可比矩阵评估中使用蒙特卡洛方法的扰动近似有关。另一种替代方法是使用正向和对偶光子密度场构建对偶雅可比矩阵。然而，这需要进一步研究来比较这些策略以及其他可能的雅可比矩阵评估方法，并开发补偿随机噪声的方法。

DSGN算法单次迭代的计算时间为28秒（使用10?光子包），而使用10?光子包的SGN需要4,547秒。另一方面，DSGN的收敛速度略慢，即平均15次迭代，而参考SGN方法平均需要8次迭代。因此，当迭代在收敛时停止，平均总计算时间从36,376秒缩短到420秒，这意味着DSGN比使用10?光子包的参考SGN快约85倍。我们还注意到，DSGN方法的训练是耗时操作，但可以在离线状态下完成，之后即可高效应用训练好的模型。

综上所述，已证明机器学习可用于补偿光声层析成像中高斯-牛顿搜索方向因蒙特卡洛随机噪声而产生的不确定性。此外，该方法能够使用比传统随机高斯-牛顿算法显著少得多的光子包以良好质量重建吸收和散射系数。