扩散光学层析成像（DOT）逆问题病态性严重。传统方法成像精度不高，计算耗时，制约了 DOT 技术的临床应用。因此，本文提出一种基于栈式自编码器（SAE）的 DOT 逆问题求解方法。首先采用传统 SAE 方法代替迭代方法进行逆问题计算，其次改进了 SAE 神经网络的输出结构，使用单输出 SAE 降低单个网络负担，最后将改进 SAE 方法与传统列文伯格-马夸尔特（LM）迭代方法、传统 SAE 方法进行仿真比较。结果表明，本文所提方法逆问题求解平均用时只有 LM 迭代方法的 1.67%，实验模型下均方误差（MSE）值较迭代方法降低了 46.21%，较传统 SAE 方法降低了 61.53%，图像相关系数（ICC）值较传统方法提升了 4.03%，较传统 SAE 方法提升了 18.7%，并且在 3% 噪声条件下具有良好的抗噪性。通过本文的研究结果证明，改进后的 SAE 方法相较于传统 SAE 方法具有更高的图像质量及抗噪性，同时相较传统迭代方法具有较快的计算速度，有利于神经网络在 DOT 逆问题计算中的应用。
引用本文： 田文旭, 杨丹, 魏竹林, 王骄. 基于改进栈式自编码器的扩散光学层析成像逆问题求解方法研究. 生物医学工程学杂志, 2021, 38(4): 774-782. doi: 10.7507/1001-5515.202010041 复制
引言
扩散光学层析成像（diffuse optical tomography, DOT）作为一种新兴的断层成像技术，其原理是利用近红外光在组织体内传播时与组织的吸收散射效应进行成像，并结合算法获取被测对象光学参数（吸收系数和散射系数）。正常和病变组织的光学参数具有明显差异，且 DOT 具有无损、实时、造价低等优点，因此 DOT 技术在动作识别、乳腺肿瘤、脑出血等疾病诊断方面具有潜在的应用价值[1-3]。DOT 技术的研究可分为正问题和逆问题[4]，正问题是指已知光源和被测组织的光学参数分布，根据光在组织中的传播模型，估计发射光的边界辐射强度；逆问题则是已知组织体表面光源分布，通过获取测量边界的辐射强度，重构被测目标成像区域内光学参数的分布。然而，逆问题具有严重的病态性，测量值任意微小的扰动下都会产生很大的波动，快速有效的逆问题求解方法是本领域的研究热点[5]。
传统的 DOT 逆问题求解多以迭代算法为主，其中有代表性的是代数重建算法（algebraic reconstruction technique, ART），它是一种基于矩阵的快速稳定逆问题求解算法，包含稀疏矩阵大规模线性方程组求解，1995 年被 Arridge[6]引入 DOT 逆问题中，取得较好效果，可在投影数据不完备且投影角度不均匀的情况下也能获得髙质量图像，但该方法计算量大，占用空间内存多，导致逆问题计算速度慢。佟珊珊[7]和唐锦萍[8]提出了混合总变差与 L1、Lp范数正则化的迭代策略，在保持边界锐利性及稳定性的同时，提高了计算精度，但传统迭代方法计算速度缓慢及强扰动计算误差大的问题并没有得到较好的解决。Vidal-Rosas 等[9-10]提出了在逆问题中使用降阶模型用于加速逆问题的迭代速度，通过求解光学参数与辐射强度的非线性关系简化逆问题计算过程，但降阶模型仍然存在可解释范围差的问题。近年来，神经网络的方法被应用到 DOT 逆问题求解中，Feng 等[11]提出了基于反向传播神经网络（back propagation neural networks, BPNN）的 DOT 逆问题计算方案，通过直接建立边界光学信息与成像区域内部的神经网络求解成像区域光学参数分布，避免了逆问题求解的不适定和速度慢的问题，但 BPNN 扰动吸收系数特异性要求高，泛化能力不足的问题无法得到解决。Yoo 等[12]和 Yedder 等[13]提出了基于卷积神经网络（convolutional neural networks, CNN）的逆问题求解方案，提高了神经网络的泛化能力，但 CNN 仍然存在网络参数过多，计算量大，对计算机要求高等问题。
综上，为了解决现有的 DOT 逆问题求解泛化能力差、计算精度低等问题，本文提出了一种基于改进栈式自编码器（stacked autoencoder, SAE）的逆问题求解方法。通过 SAE 对输入数据进行降维及特征提取，并使用单输出 SAE 降低传统 SAE 网络的任务量，最终期望本文所提的改进 SAE 方法的逆问题求解方法相较于传统神经网络能够具有更高的图像精度及更好的泛化能力，有利于神经网络在 DOT 领域中的实际应用。
1 DOT 求解
DOT 求解正问题是根据光在组织中的传播模型，给定成像区域内不同位置 r 处的光学参数 μ(r)，得到边界光辐射强度 M，其中 μ（r）与 M 的关系如式（1）～式（2）所示：
 |
 |
其中，μa(r)、μs(r) 分别表示位于成像区域 r 处的吸收系数和散射系数，S、D 表示光源和探测器的数量，
表示光学参数μ(r) 映射到边界辐射强度M的正向算子，使用辐射传输方程（radiative transfer equation, RTE）近似求解，具体形式如式（3）所示：
 |
在体元 X 内，u(r,ω) 表示在 r 处 ω 方向上的单位立体角 dσ(ω) 发射的光子密度，q(r) 表示 r 处入射的光子密度。r 处的扩散系数 κ(r) 如式（4）所示：
 |
其中，g 为已知的各向异性系数。在成像区域边界 ξ 处光子传输如式（5）～式（7）所示：
 |
 |
 |
其中，A 表示边界处内部与外部折射率之差，n表示边界处的法向向量，θ（ξ）表示边界 ξ 处所有方向光子密度的积分。DOT 求解逆问题为，给定边界光辐射强度M，计算成像区域内的光学参数μ（r）分布。选择吉洪诺夫（Tikhonov）正则化，正则化矩阵取单位矩阵I时，称为列文伯格-马夸尔特（levenberg-marquardt, LM）方法，此时 DOT 逆问题求解如式（8）所示：
 |
其中，
表示成像区域内光学参数分布，P′（μ（k））在此简写为J。每次迭代计算 P（μ（k+1））与M的差值，当满足极小值要求时，迭代结束。
2 提出的方法
2.1 SAE 的 DOT 逆问题求解
SAE 由多个自编码器（auto-encoder, AE）堆栈产生，其主要应用在特征提取、数据降维去噪及图像识别等问题中[14-16]，训练过程如图 1 所示。
图1
SAE 神经网络基本结构
Figure1.
Basic structure of SAE neural network
图 1 中由 AE 网络构成传统 SAE 神经网络并完成神经网络训练流程如下：
步骤 1：计算 AE 网络输入层到隐藏层的输出。从输入层到隐藏层的过程称为编码过程，输入层到隐藏层之间的关系如式（9）～式（11）所示：
 |
 |
 |
其中，上式及图 1 中
为输入层一个样本的 m 维向量，
为隐藏层的 n 维向量，W为从输入层到隐藏层的权值矩阵，
为隐藏层的偏置向量。
为隐藏层的激活函数，本文使用 tanh 函数作为激活函数，tanh 函数如式（12）所示：
 |
步骤 2：计算 AE 网络隐藏层到输出层的输出。隐藏层到输出层之间的关系如式（13）～式（15）所示：
 '/> |
 '/> |
 '/> |
其中，
为与输入向量维度相同的向量，W'为隐藏层到输出层的权值矩阵，b'为偏置向量。
步骤 3：优化 AE 网络。通过反向传播算法迭代更新权值矩阵神经网络参数，使输出向量
尽可能等于输入向量x，损失函数 L（·）如式（16）所示：
 |
步骤 4：AE 网络堆栈。将隐藏层h作为输入，重复步骤 1～步骤 3 训练图 1 中新的隐藏层
，得到更低维度的信息空间。训练完成后，舍弃当前输出层，在隐藏层后加入一个新的输出层，此时得到完整的传统 SAE 神经网络结构。
步骤 5：输出层的有监督训练。如图 1 所示，在无监督训练完成后，利用有标签数据集采用有监督学习算法对神经网络参数进行调优，最终输出 t 维数据集
。
2.2 改进 SAE 的 DOT 逆问题求解
使用 SAE 方法完成 DOT 逆问题求解的过程如图 2 所示。其中 SAE 神经网络的数据输入为已知的 S × D 维归一化光辐射强度 M，输出数据为 n 维有限元节点的 μa，传统的 SAE 方法将 DOT 逆问题视为回归问题，将输出数据视为连续值，按照图 1 的步骤完成 SAE 神经网络训练。
图2
基于传统 SAE 方法的 DOT 逆问题求解
Figure2.
The calculation of DOT inverse problem based on traditional SAE method
SAE 方法拟合 DOT 逆问题参数分布时，需要拟合所有节点参数分布，SAE 神经网络负载高，导致输出拟合结果差。因此，本文提出了一种改进 SAE 方法求解 DOT 逆问题。该方法采用多个子神经网络计算多个有限元节点，其中每个 SAE 子网络只计算一个节点，计算过程如图 3 所示。
图3
基于改进 SAE 方法的逆问题计算求解
Figure3.
The calculation process of inverse problem based on improved SAE method
在图 3 中输入光强信号仍为 S×D 维归一化辐射强度 M，输出数据为 1 维的单一有限元节点，此时将输入数据同时送入 n 个 SAE 网络，每个网络都是独立的，同时每个网络只计算某一特定有限元节点。
改进 SAE 的数据处理方式如图 4 所示，在改进的 SAE 方法中，将得到的输出数据分段表示为独热编码，SAE 子网络将 μa(r) 的值 x 按大小分为 z 个区间，同时判断节点 μa（r）所处区间，从而将回归问题改善为分类问题。改进 SAE 中输出数据采用 softmax 函数，如式（17）所示：
图4
μa 分布的独热编码表示
Figure4.
One-hot encoding representation of μa
 |
其中，z 表示将 μa(r) 分为 z 个区间，f(xi) 表示 μa(r) 属于第 i 区间的概率。
3 仿真实验
3.1 参数设置
仿真成像区域如图 5 所示，其具体设计如下：
图5
成像区域光源与探测器放置位置（单光源多探测器工 作模式）
Figure5.
The placement of the light sources and detectors in the imaging area (single light source and multiple detectors working mode)
① 传统 SAE 网络使用四层结构，每层神经元个数设置为{240, 480, 800, 1 785}。初始学习率为 0.01，迭代 1 000 次。
② 改进 SAE 网络使用四层结构，每层神经元个数设置为{240, 130, 80, 10}。初始学习率为 0.01，迭代 1 000 次，共 1 785 个单输出 SAE 子神经网络。
③ 在图中背景 μa 和 μs 分别设置为 0.01 mm–1和 1 mm–1。成像区域半径为 43 mm，在边界位置处按角度均匀分配 16 个光源和 16 个探测器。
④ 光源的辐射强度为 1 W/mm2，有限元网格数量为 1 785，使用有限元分析软件 NIRFAST 9.1（Kitware Inc., 美国）得到的 10 000 组数据作为训练样本[17-18]。实验模型中的扰动半径为 6～18 mm 随机值，μa 为 0.01～0.1 mm–1随机值，扰动位置随机分布在成像区域内，扰动的数量随机 1～2 个。
3.2 性能比较及结果分析
为了验证所提方法对 DOT 逆问题的求解性能，选取四种扰动模型，扰动的吸收系数分别为 0.09、0.08、0.07、0.02 mm–1，四种模型的设置如图 6 所示。其中模型一为高吸收系数单一扰动，扰动中心位于 (10, 0)，模型二为相切的两个扰动球，用于验证不同方法对邻近的扰动体分辨能力。模型三为两个相距 3 mm 的扰动球，用于验证两个扰动存在一定距离的条件下不同方法的分辨能力。模型四为位于中心的低吸收系数扰动球，用于验证在中心深区域低吸收系数扰动的分辨能力。
图6
模型背景及扰动设置
Figure6.
Model background and disturbance settings
在无噪声情况下，改进 SAE 方法、LM 方法及 SAE 方法的结果比较如图 7 所示。对于模型一，由 LM 方法得到的图像边界较为模糊，存在分层现象；这是由于 LM 方法采用 L2范数作为惩罚项，导致出现平滑现象。而 SAE 方法会过大估计实际扰动范围，这是由于网络对正问题数据进行预估判断时产生的错误计算结果。改进 SAE 方法扰动轮廓具有更小形变，但是会存在部分有限元节点判断错误的情况，进而产生暗点。由于改进 SAE 方法由多个子系统构成，某个系统的错误判断并不会影响到其他系统对输出的判断，虽然改进 SAE 方法仍然会有部分节点判断错误产生暗点，但是整体图像的轮廓和特征并不会受到影响，大部分子系统仍能判断正确，保持锐利。对于模型二，LM 方法由于 L2范数惩罚限制导致过于平滑的两个扰动几乎融合到一起，传统 SAE 方法仍然会错误估计扰动的大小和轮廓，也导致了两个扰动的融合，而改进 SAE 方法则能较好获得图像边界，图像较为稳定，两个扰动的边界十分清晰。模型三中两个扰动相距较远，此时 LM 及传统 SAE 方法得到改善，能够分辨两个扰动的位置，但图像分层及错误估计现象仍然严重，改进 SAE 方法仍然保持稳定。对于模型四，扰动吸收系数十分接近背景，此时传统 LM 方法没有受到过多影响，但传统 SAE 方法已经无法判断扰动的大小，改进 SAE 方法也受到一定影响，有更多的节点存在判断错误的情况，但图像边界仍然锐利，大小判断正确，总体的稳定性仍然超过传统 SAE 方法。
图7
未加噪声仿真结果比较
Figure7.
Comparison of simulation results without noise
DOT 的硬件系统的噪声通常在 1% 以下，即使在未添加低频调制的直流信号检测中，噪声也可以限制在 3% 以内[19]。在 3% 随机噪声情况下三种方法的结果比较如图 8 所示，在低信噪比条件下，LM 方法在高吸收扰动条件下几乎不会受到影响，但随着扰动吸收系数降低，LM 方法开始产生大量伪影，传统 SAE 方法则受噪声影响大，图像变形严重，当存在两个扰动时，图像的错误估计直接影响了对目标的识别。改进 SAE 方法表现出良好的抗噪能力，在低信噪比条件下也能有清晰的边界，图像变形很小。
图8
3% 噪声仿真结果比较
Figure8.
Comparison of simulation results without 3% noise
为了比较所提方法逆问题求解性能，采用了图像相关系数（image correlation coefficient, ICC）和均方误差（mean square error, MSE）评价求解的图像，ICC 计算方法如式（18）所示：
 |
MSE 计算方法如式（19）所示：
 |
其中，
表示实际图像和结果图像的 μa 值，
、
表示实际图像和结果图像 μa 的平均值。
如表 1 所示，在未添加噪声时，三种方法在单一扰动条件下都具有较高的 ICC 值及较低的 MSE 值；但在两个扰动距离较近时，LM 方法及传统 SAE 方法的 ICC 值下降，MSE 值变大，图像质量严重下降，四个模型条件下，改进 SAE 方法的 MSE 值较 LM 方法降低了 46.21%，较传统 SAE 方法降低了 61.53%，ICC 值较 LM 方法提升了 4.03%，较传统 SAE 方法提升了 18.79%，这与图 7 结果一致。如表 2 所示，加入噪声后对传统 SAE 方法图像质量影响严重，ICC 值明显下降，MSE 值变大，而改进 SAE 方法的 ICC 值与 MSE 值则无明显变化，表现出更强的抗噪声能力，四个模型条件下，改进 SAE 方法的 MSE 值较 LM 方法降低了 49.98%，较传统 SAE 方法降低了 71.94%，ICC 值较 LM 方法提升了 6.38%，较传统 SAE 方法提升了 57.15%，与图 8 结果相符合。
表1
MSE 及 ICC 评价结果（未添加噪声）
Table1.
MSE and ICC evaluation results (no noise added)
表2
MSE 及 ICC 评价结果（3% 随机噪声）
Table2.
MSE and ICC evaluation results (3% random noise)
为了验证不同方法的逆问题求解速度，本文采用了 200 个包含随机扰动的模型进行实验，此时 LM 方法、SAE 方法及改进 SAE 方法逆问题计算所用时间 t 如图 9 所示，改进 SAE 方法平均用时只有 LM 方法的 1.67%，传统 SAE 方法由于更低的计算量，其时间只有 LM 方法的 0.01%，SAE 及其改进方法的求解速度远高于 LM 方法。
图9
不同扰动模型下逆问题求解时间结果
Figure9.
Results of solving inverse problems under different distributions
4 结论
本文提出了改进 SAE 方法的 DOT 逆问题求解。该方法基于神经网络对抽象特征的学习，得到光辐射强度与单一节点光学参数之间的非线性特征。虽然基于 SAE 神经网络的方法需要长久的训练，但训练结束后网络具有极高的计算速度，且病态性的问题也得以避免。另一方面，SAE 神经网络的抗噪性和图像质量受到了制约，少量的噪声就能对网络产生很大影响，单一网络完成所有节点的求解极大加重了网络负担。为了降低网络负担，本文使用多个 SAE 子网络完成单一节点重建，有效降低了网络计算量，同时保证了各个节点计算的独立性。单一网络的估计错误不会影响到其他网络，从而使图像质量的失真以像素点为单位，而非以区域为单位，提升了 SAE 网络的重建精准度。本文提出的新方法避免了 LM 方法逆问题病态性与高耗时的同时，提高了传统 SAE 方法逆问题求解的 ICC 值，降低了 MSE 值，提升了图像质量及抗噪声能力，促进了 DOT 逆问题计算中神经网络的应用。未来为了实现所提方法的硬件化，需要进一步研究降低单输出 SAE 网络的复杂度，提高计算速度。
利益冲突声明：本文全体作者均声明不存在利益冲突。
引言
扩散光学层析成像（diffuse optical tomography, DOT）作为一种新兴的断层成像技术，其原理是利用近红外光在组织体内传播时与组织的吸收散射效应进行成像，并结合算法获取被测对象光学参数（吸收系数和散射系数）。正常和病变组织的光学参数具有明显差异，且 DOT 具有无损、实时、造价低等优点，因此 DOT 技术在动作识别、乳腺肿瘤、脑出血等疾病诊断方面具有潜在的应用价值[1-3]。DOT 技术的研究可分为正问题和逆问题[4]，正问题是指已知光源和被测组织的光学参数分布，根据光在组织中的传播模型，估计发射光的边界辐射强度；逆问题则是已知组织体表面光源分布，通过获取测量边界的辐射强度，重构被测目标成像区域内光学参数的分布。然而，逆问题具有严重的病态性，测量值任意微小的扰动下都会产生很大的波动，快速有效的逆问题求解方法是本领域的研究热点[5]。
传统的 DOT 逆问题求解多以迭代算法为主，其中有代表性的是代数重建算法（algebraic reconstruction technique, ART），它是一种基于矩阵的快速稳定逆问题求解算法，包含稀疏矩阵大规模线性方程组求解，1995 年被 Arridge[6]引入 DOT 逆问题中，取得较好效果，可在投影数据不完备且投影角度不均匀的情况下也能获得髙质量图像，但该方法计算量大，占用空间内存多，导致逆问题计算速度慢。佟珊珊[7]和唐锦萍[8]提出了混合总变差与 L1、Lp范数正则化的迭代策略，在保持边界锐利性及稳定性的同时，提高了计算精度，但传统迭代方法计算速度缓慢及强扰动计算误差大的问题并没有得到较好的解决。Vidal-Rosas 等[9-10]提出了在逆问题中使用降阶模型用于加速逆问题的迭代速度，通过求解光学参数与辐射强度的非线性关系简化逆问题计算过程，但降阶模型仍然存在可解释范围差的问题。近年来，神经网络的方法被应用到 DOT 逆问题求解中，Feng 等[11]提出了基于反向传播神经网络（back propagation neural networks, BPNN）的 DOT 逆问题计算方案，通过直接建立边界光学信息与成像区域内部的神经网络求解成像区域光学参数分布，避免了逆问题求解的不适定和速度慢的问题，但 BPNN 扰动吸收系数特异性要求高，泛化能力不足的问题无法得到解决。Yoo 等[12]和 Yedder 等[13]提出了基于卷积神经网络（convolutional neural networks, CNN）的逆问题求解方案，提高了神经网络的泛化能力，但 CNN 仍然存在网络参数过多，计算量大，对计算机要求高等问题。
综上，为了解决现有的 DOT 逆问题求解泛化能力差、计算精度低等问题，本文提出了一种基于改进栈式自编码器（stacked autoencoder, SAE）的逆问题求解方法。通过 SAE 对输入数据进行降维及特征提取，并使用单输出 SAE 降低传统 SAE 网络的任务量，最终期望本文所提的改进 SAE 方法的逆问题求解方法相较于传统神经网络能够具有更高的图像精度及更好的泛化能力，有利于神经网络在 DOT 领域中的实际应用。
1 DOT 求解
DOT 求解正问题是根据光在组织中的传播模型，给定成像区域内不同位置 r 处的光学参数 μ(r)，得到边界光辐射强度 M，其中 μ（r）与 M 的关系如式（1）～式（2）所示：
|
|
其中，μa(r)、μs(r) 分别表示位于成像区域 r 处的吸收系数和散射系数，S、D 表示光源和探测器的数量，表示光学参数μ(r) 映射到边界辐射强度M的正向算子，使用辐射传输方程（radiative transfer equation, RTE）近似求解，具体形式如式（3）所示：
|
在体元 X 内，u(r,ω) 表示在 r 处 ω 方向上的单位立体角 dσ(ω) 发射的光子密度，q(r) 表示 r 处入射的光子密度。r 处的扩散系数 κ(r) 如式（4）所示：
|
其中，g 为已知的各向异性系数。在成像区域边界 ξ 处光子传输如式（5）～式（7）所示：
|
|
|
其中，A 表示边界处内部与外部折射率之差，n表示边界处的法向向量，θ（ξ）表示边界 ξ 处所有方向光子密度的积分。DOT 求解逆问题为，给定边界光辐射强度M，计算成像区域内的光学参数μ（r）分布。选择吉洪诺夫（Tikhonov）正则化，正则化矩阵取单位矩阵I时，称为列文伯格-马夸尔特（levenberg-marquardt, LM）方法，此时 DOT 逆问题求解如式（8）所示：
|
其中，表示成像区域内光学参数分布，P′（μ（k））在此简写为J。每次迭代计算 P（μ（k+1））与M的差值，当满足极小值要求时，迭代结束。
2 提出的方法
2.1 SAE 的 DOT 逆问题求解
SAE 由多个自编码器（auto-encoder, AE）堆栈产生，其主要应用在特征提取、数据降维去噪及图像识别等问题中[14-16]，训练过程如图 1 所示。
图1
SAE 神经网络基本结构
Figure1.
Basic structure of SAE neural network
图 1 中由 AE 网络构成传统 SAE 神经网络并完成神经网络训练流程如下：
步骤 1：计算 AE 网络输入层到隐藏层的输出。从输入层到隐藏层的过程称为编码过程，输入层到隐藏层之间的关系如式（9）～式（11）所示：
|
|
|
其中，上式及图 1 中为输入层一个样本的 m 维向量，为隐藏层的 n 维向量，W为从输入层到隐藏层的权值矩阵，为隐藏层的偏置向量。为隐藏层的激活函数，本文使用 tanh 函数作为激活函数，tanh 函数如式（12）所示：
|
步骤 2：计算 AE 网络隐藏层到输出层的输出。隐藏层到输出层之间的关系如式（13）～式（15）所示：
| '/> |
| '/> |
| '/> |
其中，为与输入向量维度相同的向量，W'为隐藏层到输出层的权值矩阵，b'为偏置向量。
步骤 3：优化 AE 网络。通过反向传播算法迭代更新权值矩阵神经网络参数，使输出向量尽可能等于输入向量x，损失函数 L（·）如式（16）所示：
|
步骤 4：AE 网络堆栈。将隐藏层h作为输入，重复步骤 1～步骤 3 训练图 1 中新的隐藏层，得到更低维度的信息空间。训练完成后，舍弃当前输出层，在隐藏层后加入一个新的输出层，此时得到完整的传统 SAE 神经网络结构。
步骤 5：输出层的有监督训练。如图 1 所示，在无监督训练完成后，利用有标签数据集采用有监督学习算法对神经网络参数进行调优，最终输出 t 维数据集。
2.2 改进 SAE 的 DOT 逆问题求解
使用 SAE 方法完成 DOT 逆问题求解的过程如图 2 所示。其中 SAE 神经网络的数据输入为已知的 S × D 维归一化光辐射强度 M，输出数据为 n 维有限元节点的 μa，传统的 SAE 方法将 DOT 逆问题视为回归问题，将输出数据视为连续值，按照图 1 的步骤完成 SAE 神经网络训练。
图2
基于传统 SAE 方法的 DOT 逆问题求解
Figure2.
The calculation of DOT inverse problem based on traditional SAE method
SAE 方法拟合 DOT 逆问题参数分布时，需要拟合所有节点参数分布，SAE 神经网络负载高，导致输出拟合结果差。因此，本文提出了一种改进 SAE 方法求解 DOT 逆问题。该方法采用多个子神经网络计算多个有限元节点，其中每个 SAE 子网络只计算一个节点，计算过程如图 3 所示。
图3
基于改进 SAE 方法的逆问题计算求解
Figure3.
The calculation process of inverse problem based on improved SAE method
在图 3 中输入光强信号仍为 S×D 维归一化辐射强度 M，输出数据为 1 维的单一有限元节点，此时将输入数据同时送入 n 个 SAE 网络，每个网络都是独立的，同时每个网络只计算某一特定有限元节点。
改进 SAE 的数据处理方式如图 4 所示，在改进的 SAE 方法中，将得到的输出数据分段表示为独热编码，SAE 子网络将 μa(r) 的值 x 按大小分为 z 个区间，同时判断节点 μa（r）所处区间，从而将回归问题改善为分类问题。改进 SAE 中输出数据采用 softmax 函数，如式（17）所示：
图4
μa 分布的独热编码表示
Figure4.
One-hot encoding representation of μa
|
其中，z 表示将 μa(r) 分为 z 个区间，f(xi) 表示 μa(r) 属于第 i 区间的概率。
3 仿真实验
3.1 参数设置
仿真成像区域如图 5 所示，其具体设计如下：
图5
成像区域光源与探测器放置位置（单光源多探测器工 作模式）
Figure5.
The placement of the light sources and detectors in the imaging area (single light source and multiple detectors working mode)
① 传统 SAE 网络使用四层结构，每层神经元个数设置为{240, 480, 800, 1 785}。初始学习率为 0.01，迭代 1 000 次。
② 改进 SAE 网络使用四层结构，每层神经元个数设置为{240, 130, 80, 10}。初始学习率为 0.01，迭代 1 000 次，共 1 785 个单输出 SAE 子神经网络。
③ 在图中背景 μa 和 μs 分别设置为 0.01 mm–1和 1 mm–1。成像区域半径为 43 mm，在边界位置处按角度均匀分配 16 个光源和 16 个探测器。
④ 光源的辐射强度为 1 W/mm2，有限元网格数量为 1 785，使用有限元分析软件 NIRFAST 9.1（Kitware Inc., 美国）得到的 10 000 组数据作为训练样本[17-18]。实验模型中的扰动半径为 6～18 mm 随机值，μa 为 0.01～0.1 mm–1随机值，扰动位置随机分布在成像区域内，扰动的数量随机 1～2 个。
3.2 性能比较及结果分析
为了验证所提方法对 DOT 逆问题的求解性能，选取四种扰动模型，扰动的吸收系数分别为 0.09、0.08、0.07、0.02 mm–1，四种模型的设置如图 6 所示。其中模型一为高吸收系数单一扰动，扰动中心位于 (10, 0)，模型二为相切的两个扰动球，用于验证不同方法对邻近的扰动体分辨能力。模型三为两个相距 3 mm 的扰动球，用于验证两个扰动存在一定距离的条件下不同方法的分辨能力。模型四为位于中心的低吸收系数扰动球，用于验证在中心深区域低吸收系数扰动的分辨能力。
图6
模型背景及扰动设置
Figure6.
Model background and disturbance settings
在无噪声情况下，改进 SAE 方法、LM 方法及 SAE 方法的结果比较如图 7 所示。对于模型一，由 LM 方法得到的图像边界较为模糊，存在分层现象；这是由于 LM 方法采用 L2范数作为惩罚项，导致出现平滑现象。而 SAE 方法会过大估计实际扰动范围，这是由于网络对正问题数据进行预估判断时产生的错误计算结果。改进 SAE 方法扰动轮廓具有更小形变，但是会存在部分有限元节点判断错误的情况，进而产生暗点。由于改进 SAE 方法由多个子系统构成，某个系统的错误判断并不会影响到其他系统对输出的判断，虽然改进 SAE 方法仍然会有部分节点判断错误产生暗点，但是整体图像的轮廓和特征并不会受到影响，大部分子系统仍能判断正确，保持锐利。对于模型二，LM 方法由于 L2范数惩罚限制导致过于平滑的两个扰动几乎融合到一起，传统 SAE 方法仍然会错误估计扰动的大小和轮廓，也导致了两个扰动的融合，而改进 SAE 方法则能较好获得图像边界，图像较为稳定，两个扰动的边界十分清晰。模型三中两个扰动相距较远，此时 LM 及传统 SAE 方法得到改善，能够分辨两个扰动的位置，但图像分层及错误估计现象仍然严重，改进 SAE 方法仍然保持稳定。对于模型四，扰动吸收系数十分接近背景，此时传统 LM 方法没有受到过多影响，但传统 SAE 方法已经无法判断扰动的大小，改进 SAE 方法也受到一定影响，有更多的节点存在判断错误的情况，但图像边界仍然锐利，大小判断正确，总体的稳定性仍然超过传统 SAE 方法。
图7
未加噪声仿真结果比较
Figure7.
Comparison of simulation results without noise
DOT 的硬件系统的噪声通常在 1% 以下，即使在未添加低频调制的直流信号检测中，噪声也可以限制在 3% 以内[19]。在 3% 随机噪声情况下三种方法的结果比较如图 8 所示，在低信噪比条件下，LM 方法在高吸收扰动条件下几乎不会受到影响，但随着扰动吸收系数降低，LM 方法开始产生大量伪影，传统 SAE 方法则受噪声影响大，图像变形严重，当存在两个扰动时，图像的错误估计直接影响了对目标的识别。改进 SAE 方法表现出良好的抗噪能力，在低信噪比条件下也能有清晰的边界，图像变形很小。
图8
3% 噪声仿真结果比较
Figure8.
Comparison of simulation results without 3% noise
为了比较所提方法逆问题求解性能，采用了图像相关系数（image correlation coefficient, ICC）和均方误差（mean square error, MSE）评价求解的图像，ICC 计算方法如式（18）所示：
|
MSE 计算方法如式（19）所示：
|
其中，表示实际图像和结果图像的 μa 值，、 表示实际图像和结果图像 μa 的平均值。
如表 1 所示，在未添加噪声时，三种方法在单一扰动条件下都具有较高的 ICC 值及较低的 MSE 值；但在两个扰动距离较近时，LM 方法及传统 SAE 方法的 ICC 值下降，MSE 值变大，图像质量严重下降，四个模型条件下，改进 SAE 方法的 MSE 值较 LM 方法降低了 46.21%，较传统 SAE 方法降低了 61.53%，ICC 值较 LM 方法提升了 4.03%，较传统 SAE 方法提升了 18.79%，这与图 7 结果一致。如表 2 所示，加入噪声后对传统 SAE 方法图像质量影响严重，ICC 值明显下降，MSE 值变大，而改进 SAE 方法的 ICC 值与 MSE 值则无明显变化，表现出更强的抗噪声能力，四个模型条件下，改进 SAE 方法的 MSE 值较 LM 方法降低了 49.98%，较传统 SAE 方法降低了 71.94%，ICC 值较 LM 方法提升了 6.38%，较传统 SAE 方法提升了 57.15%，与图 8 结果相符合。
表1
MSE 及 ICC 评价结果（未添加噪声）
Table1.
MSE and ICC evaluation results (no noise added)
表2
MSE 及 ICC 评价结果（3% 随机噪声）
Table2.
MSE and ICC evaluation results (3% random noise)
为了验证不同方法的逆问题求解速度，本文采用了 200 个包含随机扰动的模型进行实验，此时 LM 方法、SAE 方法及改进 SAE 方法逆问题计算所用时间 t 如图 9 所示，改进 SAE 方法平均用时只有 LM 方法的 1.67%，传统 SAE 方法由于更低的计算量，其时间只有 LM 方法的 0.01%，SAE 及其改进方法的求解速度远高于 LM 方法。
图9
不同扰动模型下逆问题求解时间结果
Figure9.
Results of solving inverse problems under different distributions
4 结论
本文提出了改进 SAE 方法的 DOT 逆问题求解。该方法基于神经网络对抽象特征的学习，得到光辐射强度与单一节点光学参数之间的非线性特征。虽然基于 SAE 神经网络的方法需要长久的训练，但训练结束后网络具有极高的计算速度，且病态性的问题也得以避免。另一方面，SAE 神经网络的抗噪性和图像质量受到了制约，少量的噪声就能对网络产生很大影响，单一网络完成所有节点的求解极大加重了网络负担。为了降低网络负担，本文使用多个 SAE 子网络完成单一节点重建，有效降低了网络计算量，同时保证了各个节点计算的独立性。单一网络的估计错误不会影响到其他网络，从而使图像质量的失真以像素点为单位，而非以区域为单位，提升了 SAE 网络的重建精准度。本文提出的新方法避免了 LM 方法逆问题病态性与高耗时的同时，提高了传统 SAE 方法逆问题求解的 ICC 值，降低了 MSE 值，提升了图像质量及抗噪声能力，促进了 DOT 逆问题计算中神经网络的应用。未来为了实现所提方法的硬件化，需要进一步研究降低单输出 SAE 网络的复杂度，提高计算速度。
利益冲突声明：本文全体作者均声明不存在利益冲突。
