基于有限元与深度学习的BLT图像重建方法

张暄轩，曹 旭，张久楼，张 琳

(1.西安邮电大学 通信与信息工程学院，陕西 西安 710121;2.西安电子科技大学 生命科学技术学院，陕西 西安 710071；3.南京医科大学 医学影像学院，江苏 南京 210029；4.山东师范大学 信息科学与工程学院，山东 济南 250014)

生物发光断层成像(Bioluminescence Tomography,BLT)是一种预临床成像技术，其通过采集生物体表透射出的生物发光，能够非侵入式地揭示活体动物体内的功能活动[1-2]。生物发光位于可见光或近红外波段，产生自生物体内的某些生物发光反应，如荧光素酶催化反应。生物发光的机理来源于分子层面的化学反应，因此生物发光光源的分布能够反映这些特异性分子的分布，成为一种重要的分子影像技术。作为一种分子影像技术，BLT通过编码生物发光蛋白基因，能够可视化地体现多种多样的细胞层面和亚细胞层面的生物过程，业已在诸多生物医学研究中展现出强大的潜力，如肿瘤放疗[3]、癌细胞检测[4]和治疗评价[5]。

BLT的成像理念类似于诸如扩散光学断层成像[6](Diffuse Optical Tomography,DOT)和荧光分子断层成像[7](Fluorescence Molecular Tomography,FMT)的其他光学断层成像技术，都是通过探测生物体表的光子密度，再利用数学模型反向推算出以某种物理参数为对比度的断层图像。不同的是，BLT的光源位于生物体内并且不需要外部激发。考虑不需要外部激发光源，BLT的测量数据不会像FMT那样被背景光干扰，使得BLT具有更高的灵敏度。在BLT的成像过程中，需要先在生物体表测量生物发光的光子密度，再建立某种数学模型描述生物发光在生物体内的扩散过程。最后，利用测量数据求解该数学模型的逆向问题重建出生物发光光源的分布。扩散方程[8]是最常用的正向模型，其是辐射传输方程的一种近似，方程形式为一个椭圆型偏微分方程，能够被许多数值方法求解，如有限元方法[8](Finite Element Method,FEM)。基于正向模型，能够通过最小化测量数据与正向模型的预测值之间的差距建立逆向问题。类似于其他光学断层成像技术，BLT的逆向问题也具有高度病态性，其由测量数据的维度与待重建参数的维度之间的不平衡造成[9]。在BLT中，生物体表所获取的测量数据具有部分相关性，且不能像DOT和FMT中那样通过增加投影角度扩增，使得BLT的病态性问题要比DOT和FMT更加严重。为了增加BLT测量数据的多样性，通常需要采集多个波长下的生物发光图像[10]。此外，如果直接使用DOT或FMT中的算法，如采用牛顿型算法[8]进行BLT图像重建，逆向问题的病态性会恶化图像重建结果的质量。为了获取可靠的重建结果，BLT重建通常需要添加额外的信息，如利用来自其他成像模态的解剖结构信息正则迭代过程[11-13]，或是利用可行域缩小待重建参数的维度[14-16]。但是，这些方法并不总是可用的，解剖结构信息需要额外的成像模态[9]，而可行域的选择需要一定的先验信息[15]。

近年来，深度学习[17](Deep Learning，DL)受到了各个领域的广泛关注，特别是图像领域，而近些年关于深度学习在医学成像中应用的报道主要集中在图像处理或分析[18]。尝试利用深度学习构建图像重建方法的研究则集中出现于近两三年，其结果显示，基于深度学习的图像重建方法具有巨大潜力[19]。这些研究涵盖了广泛的成像模态，包括X射线计算机断层成像[20-22](X-ray Computed Tomography,X-CT)、磁共振成像[23-24](Magnetic Resonance Imaging,MRI)、正电子发射断层成像[25](Positron Emission Tomography,PET)以及多种光学断层成像技术，如DOT[26]、FMT[27]和BLT[28]。深度学习是一种数据驱动型机器学习方法，能够通过学习从数据中抽取出两个相关数据域之间的映射关系。利用机器学习策略实现图像重建时，深度学习既可以被用在图像域也可以被用在数据域。对于图像域，神经网络可以被用来获取两种不同条件下图像的转换关系，如被噪声污染的图像与无噪声图像之间的转换关系。严格来说，这种策略是一种后处理方法而不是图像重建，但通常却被用来解决某些传统图像重建问题，如X-CT中的低剂量重建[20]、稀疏角度重建[21]、有限角度重建问题[22]以及MRI中的稀疏重建[23]。对于数据域，深度学习可以被用来预处理数据或是构建一个端到端的模型直接将原始数据转换为断层图像，如Würfl等[22]提出了一种含有锥形束反投影层的深度学习方法实现端到端的X-CT有限角度重建；Zhu等[24]针对MRI提出了一种由若干卷积层和全连接层组成的端到端神经网络，称为流形逼近自动变换(Automated Transform by Manifold Approximation,AUTOMAP)，实现数据域到图像域的直接重建； Yedder等[26]尝试使用端到端的卷积神经网络实现从DOT投影数据到断层图像的直接重建；Guo等[27]提出一种三维深度编码-解码网络实现端到端的FMT图像重建；Yao等[28]使用一种4层的端到端全连接神经网络实现BLT图像重建。

在以往的研究中，BLT的图像重建策略被迭代法所主导。迭代法计算耗时，且需要一个较精确的数学模型描述生物发光的传输过程。与迭代法相比，DL作为一种数据驱动型机器学习方法，一旦神经网络被训练好，即具有极快的推断速度，且能够直接从数据中抽取映射关系，不需要人为选择特殊的正向模型。DL的应用难点在于收集大量数据用以训练神经网络。文献[28]的研究显示了通过深度学习实现BLT图像重建的可行性，其所使用的数据集由蒙特卡洛方法[29]生成。蒙特卡洛方法通过大量重复单个光子在生物体内的传输过程并统计结果实现对光场分布的估计，其精度由所重复的光子个数所决定，通常需要重复上亿个光子。虽然大量的光子能够保证精度，但同时也耗费了较长的计算时间。考虑深度学习需要大量的训练数据和测试数据，计算时间问题的考量对于深度学习的应用尤为重要。BLT中被成像物体和被成像目标会随着所研究问题的不同有所变化，需要对不同的研究问题生成不同的数据集，这就增强了对数据生成速度的要求。使用FEM求解扩散方程[8]的策略可以作为另一种数据生成的可行方案，相比蒙特卡洛方法，这种方案具有极快的计算速度。因此，更适合频繁生成大规模数据集。

1 基于FEM的BLT数据生成

1.1 扩散方程

扩散方程配以FEM是一种被广泛使用的光子扩散模拟方法。其中扩散方程是一个椭圆型偏微分方程，源自辐射传输方程，是其近似形式，稳态扩散方程式[8]为

(1)

式中：D和μa分别表示扩散系数和吸收系数；Φ(r)表示空间位置r处的光子密度，r∈Ω；X(r)表示光源项，在BLT中即待重建的未知参数；Ω表示求解的空间域，在BLT中即生物体内的空间。

1.2 有限元方程

作为一个偏微分方程，扩散方程在无限介质中的解析解有其确定形式，但具有复杂边界的介质上的解析解形式却难以获得，因而通常使用数值方法求解扩散方程。FEM是求解偏微分方程的常用数值方法之一，其需要将空间域Ω离散为一个由若干节点和单元组成的网格，并将方程中的函数使用基函数展开，由基函数所构成的线性空间被称为有限元空间。通过离散空间域形成的有限元空间具有有限维度，那么在有限元空间上展开的任意函数都可以被一组展开系数表示，由此能将一个偏微分方程转化为一个线性方程组求解，称为有限元方程。

为了描述复杂边界情况，需要引入以下Robin边界条件[8]，其表达式为

(2)

通过引入Robin边界条件，能够通过FEM将扩散方程离散为有限元方程[8]，其表达式为

KΦ=X

(3)

式中：Φ和X分别表示光子密度向量和源向量，都是维度为N的列向量，N表示网格中节点的数量；K是一个N×N的方阵，在FEM中一般被称为刚度矩阵。矩阵K以及向量X的元素是基函数的函数，表达式分别为[8]

离了林强信，景花厂没有倒。林强信不甘心，他是个精明的生意人，从不做赔了夫人又折兵的买卖，他不甘心白吃了这个哑巴亏。景花厂后来渐渐有了起色，订单多了，员工多了，林强信更不甘心了。他要击败景花厂，要阿花像只无家可归的猫，乖乖回到他的怀抱。

(4)

(5)

式中：υ表示有限元空间的基函数；i和j分别表示行或列的索引。

刚度矩阵K是一个稀疏、对称且满秩的方阵，因此有限元方程是一个具有唯一解的非其次线性方程组，可以通过诸多成熟的数值方法进行快速求解。对于某一被成像物体，只要先构建相应的几何体并且进行网格剖分，再构建有限元空间产生刚度矩阵K。同时，通过人为设计构建特定光源分布的X，最后求解有限元方程即可实现大量BLT数据的快速获取。

2 基于深度学习的BLT图像重建

2.1 深度学习

深度学习[17]是一种数据驱动型机器学习方法，其能够使用由多个层级结构结合激活单元所组成的神经网络近似任意非线性函数。深度学习有两种标准训练模式，分别称为监督学习和无监督学习。在监督学习中，神经网络需使用成对的输入数据和标签数据训练解决分类或回归问题，而无监督学习则不需要标签数据，通常用来寻找数据中的范式。总体来说，无监督学习因其自动挖掘数据范式的能力而更适合人工智能任务，但其研究相对较少，理论体系尚不完善，在许多问题中仍不能使用。与无监督学习相比，监督学习应用的更广泛，并且只要有标签数据可用，通常可以获取可靠的结果。与计算机视觉中的一些问题不同，生物医学成像对模型的稳定有着较高的要求，如X-CT和MRI图像常被用来诊断人的疾病，细微的成像图像差异会造成难以估计的损失。因此，目前大部分有关使用深度学习实现图像重建的研究均使用监督学习。

2.2 基于深度学习的光子密度分布补全

为了使真实空域的复杂形状边界能够和离散后的网格边界更好的贴合，FEM通常使用非结构网格构建有限元空间，二维情形一般使用三角形网格，而三维情形则使用四面体网格。与结构化网格不同，非结构化网格由不规则连接的节点组成，没有规则的拓扑结构，所有节点的拓扑关系需要进行记忆存储。而非结构化网格的不规则性会对图像重建产生影响[7]，其通常会因重建目标的稀疏性增加而加剧。因此，预先测试了直接使用一个全连接网络构建端到端的BLT图像重建，结果显示网格的不规则性会显著地恶化图像重建质量。为了消除网格不规则性的影响，利用神经网络实现光子密度分布的补全而非直接重建出图像。在生物体表获取的测量数据实际上是光子密度在表面的采样值。因此，光子密度分布的补全就是要从体表部分光子密度值推测出整个空间的光子密度分布，在获取完整的光子密度分布后，再通过将刚度矩阵与光子密度向量相乘获得生物发光光源分布。光子密度分布是一个连续变化的函数，基本不会受到网格不规则性的影响。受网格拓扑结构的启发，设计了一种网格拓扑层接入体表测量数据。网格拓扑层由多个子级层级结构组成，子层级之间的连接关系又由网格的拓扑关系决定。所使用的网格拓扑层和神经网络的组成结构如图1所示，使用了一个具有少量节点的粗糙三角形网格解析网络结构。根据拓扑结构，图1(左)所示的网格的节点可以被分为5层，最外层的节点是表面采样节点，测量数据是在这些节点上获取的，根据拓扑连接关系，其余内部节点可以分为4层。网格拓扑层的子层级数量及每个子层级的神经元数量与内部节点的层级结构相对应，相邻两个子层级之间的神经元使用全连接。网格拓扑层的输出是所有子层级输出的合集，该输出被接入一个由若干全连接层组成的子网络，子网络每层的神经元数为内部节点的总数。最终，子网络的输出和输入的测量数据组成完整的光子密度分布，将其与刚度矩阵相乘得到光源分布。

图1 所使用的网格拓扑层和神经网络示意图

3 实验结果与分析

3.1 实验设置

为了验证所提方法，进行了仿真实验。使用一个直径为3 cm的圆形物体作为成像介质，其被进一步离散成一个由2 168个节点和4 206个单元组成的三角形网格。通过FEM生成了3个不同数据集分别用于训练、验证和测试。训练集由2 219个样本组成，每个样本都是通过在介质中嵌入一个直径为0.2 cm的圆形光源获得，光源的位置规则排列在介质内。验证集和测试集分别由500个样本组成，这些样本均使用介质内的随机光源位置获得。每个样本都是由光子密度表面采样值和内部值所组成的成对数据。神经网络使用Python 3.6和Tensorflow搭建，训练的硬件平台为一台配有Nvidia GeForce RTX 2080 Ti显卡和16 GB内存的个人电脑。使用平方损失函数为

(6)

还使用了一个全连接网络(Fully Connected Network，FCN)[28]和基于可行域的迭代法(Optimized Permission Source Region，OPSR)[14]对测试集数据进行重建。FCN由4个全连接层组成，每层的神经元数与网格的节点数相同。与所提出的方法不同，FCN使用端到端的训练方式，即网络的输入为表面测量数据，而输出直接为光源的分布。FCN的训练策略与所提的方法相同。OPSR中的Nf参数设置为1，迭代次数设置为50。

3.2 实验结果分析

实验的重建结果如图2所示，共展示了7组具有代表性的图像。考虑所提方法使用神经网络实现光子密度分布的补全，图2也展示了光子密度的分布，所有图像均使用最大值进行了归一化。为了定量评价重建结果，引入了均方根误差(Root Mean Square Error，RMSE)，其结果如表1所示。进一步对比了处理单幅图像所需要的计算时间，其结果如表2所示。

图2 BLT重建结果

表1 训练集重建结果的RMSE

表2 单幅图像的处理时间

由图2可以看出，所提方法能从表面采样数据较好地预测介质内整个光子密度分布，即图2第3行，这种高质量的光子密度补全保证了与K相乘所得的重建结果也具有较高质量，平均RMSE仅为0.033 01，远小于FCN和OPSR，即图2第4行和表1第1行。相比之下，FCN的重建结果表现较差，即图2第5行和表1第2行，重建出的目标趋向于集中在某个或某几个孤立的点上，这种现象常出现在网格比较细小的光学断层成像重建中，其主要原因是网格的不规则性。此外，FCN在重建某些目标时会失效，如图2(e4)，这些样本的FCN输出会是一个全为零的空向量。图2第6行显示，OPSR能够重建出目标的位置，但却很难重建出其轮廓。考虑OPSR是一种迭代法，决定目标轮廓的可行域会在每次迭代时通过移除低贡献节点的方式进行收缩。因此，可行域的最终形状很大程度上依赖于算法的参数设置，包括最终节点数量及迭代次数。最终节点数量的推荐值为1[14]，也就是说可行域会随着迭代尽可能地缩小，这样的设置不需要任何先验信息。但是，最终结果会趋向于集中在一个很小的区域上，结合一定的先验信息设置最终节点数量才能使可行域最终收敛在一个合适的大小上。

由图2和表1可以看出，所提方法表现优于FCN和OPSR。FCN的结果表明，对于FEM框架下的BLT重建，由于FEM网格的不规则性，直接训练端到端的全连接神经网络很难取得较好的效果。所提方法和FCN的差别能够解释这种表现上的差异，其使用神经网络实现光子密度分布的补全，而非直接端到端的重建出光源分布。光源分布和光子密度分布的区别在于分布的形式不同，前者具有明确的边界，是一个分片连续的函数且包含大量的零值，后者则是一个连续函数且没有零值。相比稀疏的分片连续函数，连续函数受到网格不规则性的影响更弱。此外，结果表明网格拓扑层能够改善网络的性能，其模拟了光子密度分布的补全应该是一个由外向内逐步递进的过程。

与传统迭代法相比，基于深度学习的BLT图像重建方法在计算效率方面表现卓越。如表2所示，基于深度学习的两种方法重建一副图像只需不到1 ms的时间，而作为迭代法的OPSR则需要超过10 s，OPSR的大量耗时主要源于必不可少的迭代过程。相比之下，基于深度学习的方法只需要将测量数据输入至神经网络，神经网络的预测过程则只需要若干次矩阵乘法即可完成，这样不仅计算效率极高，且可以同时馈送多组测量数据进行同时计算。深度学习的这种计算高效性，源于其预测时较低的时间复杂度。与迭代法不同，深度学习的预测过程仅由若干个矩阵乘法组成，而不存在求逆过程，也不需要计算梯度，降低了计算的时间复杂度。但是，深度学习的网络模型通常较大，需要存储大量的权重参数，因此较迭代法具有更高的空间复杂度。

此外，还对初始化、初始学习率以及NFCL等3种因素对网络性能的影响进行了测试，结果如图3所示。初始化的测试通过重复4次训练过程实现，结果如图3(a)所示，权重值的初始值由正态分布确定，因此每次训练都具有不同的初始化。初始学习率的测试中分别测试了10-6、10-5、10-4和10-3等4种不同初始学习率下的网络性能，结果如图3(b)所示。NFCL的测试中分别为网格拓扑层接续了1～5层的全连接层，结果如图3(c)所示。

图3 初始化、初始学习率以及NFCL对网络性能的影响

神经网络的训练是深度学习的关键步骤，一个神经网络的表现很大程度上取决于训练数据集和训练策略。对于BLT，基于FEM的扩散方程求解长久以来被广泛用于模拟仿真和图像重建。FEM的高效计算特性适应深度学习对大规模数据和频繁训练场景的需求。对于训练策略，作者测试了3种可能会影响网络性能的因素。如图3(a)所示，网络参数的初始化对网络性能的影响非常微弱，不同初始化的验证损失基本上变化一致。图3(b)显示，初始学习率会显著影响网络性能，一个较大的初始学习率会导致过早的训练停滞，而一个较小的初始学习率会降低收敛速度，初始学习率为10-4时会带来最好的网络性能。图3(c)表明，3个全连接层能够充分应付当前数据集的训练，少于3个全连接层会导致网络性能的明显降低，而大于3个全连接层基本上不会带来性能的提升但会增加训练时间。

利用有限元方法产生训练样本可以快速产生大量训练样本。空间域固定，则K就是固定的，对于不同的样本仅需计算一次，只要设计不同的X分布即可同时获得大量样本。而利用蒙特卡洛方法产生训练样本，每个样本都需要进行独立的光子输运模拟，计算时间会随着样本数的增加成倍增长。利用有限元方法产生训练样本的缺陷是有限元的计算精度与网格的细分程度直接相关，网格越精细，则计算的精度越高，但网格密度的提升会同时增大神经网络的规模，将会增加网络训练的难度。目前，神经网络的快速训练依赖于图形处理器(Graphics Processing Unit,GPU)加速，而GPU的显存容量有明确的界限，限制了单显示卡下所能训练的网络规模，间接限制了有限元网格的精细化程度。分布式计算能够在一定程度上解决该问题，但其计算效率提升不如单显示卡训练，并且会增加硬件成本。

4 结 语

提出一种基于FEM和深度学习的BLT图像重建方法。该方法使用FEM求解扩散方程生成数据集，采用一种由网格拓扑层和全连接层组成的神经网络，通过该神经网络实现光子密度分布的补全，而非构造直接端到端的图像重建。在获取完整的光子密度分布后，再通过将刚度矩阵与光子密度向量相乘获得生物发光光源分布。仿真实验结果表明，所提方法的表现优于直接端到端的FCN和基于可行域的迭代法OPSR。