岑星星，严壮志,，吴焕迪

1 上海大学 通信与信息工程学院，上海市，200444

2 上海大学 生物医学工程研究所，上海市，200444

0 引言

荧光扩散断层成像（Fluorescence Diffuse Optical Tomography,FDOT）等光学分子影像技术，被认为是现阶段生物医学取得进一步突破的最大可能[1]。FDOT通过外部光源来激发体内的荧光分子探针使其发出荧光，并利用体表采集到的光学信号完成体内分子探针分布的重建工作。相较于其它传统成像方式，FDOT具有高敏感性、高通量、低成本、安全无创等优势[2]。近年来，随着成像系统、成像算法、光学分子探针等领域的快速发展，FDOT逐步走向成熟，并有望在预临床和临床领域取得广泛应用。

FDOT相关系统的开发一直是当前生物医学工程领域的研究热点。独立开源系统，例如SHENGHAN等[3]基于蒙特卡罗（Monte Carlo,MC）前向模型提出的用于组织体光传播模拟的MOSE平台；SCHWEIGER等[4]基于扩散方程（Diffusion Equation,DE）前向模型提出的用于组织体光传播模拟及相应光学重建的Toast++平台，都取得了广泛应用。另一方面是基于商用系统，如通过COMSOL、ANSYS等有限元分析工具进行辐射传输方程（Radiative Transfer Equation,RTE）、DE等前向模型求解，实现光成像系统的搭建，具有更好的操作体验。光学成像的相关系统平台已有很多，但大多数基于传统前向模型进行搭建，包括DE、RTE等前向模型。

传统的前向模型有着各自的优势但也存在着明显的缺陷。RTE作为麦克斯韦方程的近似，能对生物组织体中的光传播过程进行准确描述[5]。然而，RTE同时也是结构复杂的微积分方程，需要很高的计算成本，并且在复杂的生物组织应用环境中，往往无法直接通过解析法进行求解[6]。另一方面，DE作为RTE的一阶近似，具有结构简单、易于求解、计算效率高等优势，是当前FDOT应用中使用最多的前向模型[7-8]。然而，作为最低阶数的近似方程，高度的简化处理，极大地损耗了该模型的精度，并带给了DE局限性：①不能用于非散射介质中的光传播描述；②不能对光源和边界处的光传播过程进行准确描述[9]。RTE、DE等传统模型存在的明显缺陷，限制了相应成像系统的更好发展。

格子玻尔兹曼（Lattice Boltzmann,LB）方法基于分子运动理论，利用规则网格中的离散粒子运动来对复杂的物理现象进行模拟，实现相应物理量的求解[10]。相比于传统方法，LB方法的对应数学模型具有计算结构简单、程序易于实现和并行化、计算稳定性好等优势。因此，构建基于LB方法的前向模型，并将其应用于FDOT系统的搭建，是进一步推动FDOT发展的极佳选择。基于上述目的，笔者提出了基于LB前向模型的FDOT系统。

1 方法

1.1 LB前向模型

FDOT中的前向模型描述了光在生物组织中的传播过程。基于LB方法的前向模型如下所示：

其中，ω(r)=μs(r)+β(r)代表反照率，wi为权重系数，可通过局部辐射场的守恒进行求解，如下所示：

其中，ψk(s)为s的多项式。不同离散速度模型所对应的权重系数如表1所示。

表1 LB方法中不同离散速度模型对应的权重系数Tab.1 The weights of different discrete models in LB method

LB前向模型的边界节点计算由罗宾（Robin）边界条件给出。求出边界节点r∈eΩk（k为离散空间节点的索引）上的光子密度Φ(r,t)：

并和罗宾边界条件进行关联，得到如下边界计算公式[18]：

其中，D(r)=1/[3(μa(r)+(1-g)μs(r))]为光扩散系数，g为各向异性因子；是边界处的外法线向量；当组织周围填充空气，组织折射率为n时，A(r;n,n')≈(1+R(r)/(1-R(r))，其中的R(r)≈-1.439 9n-2+0.709 9n-1+0.668 1+0.063 6n。

设置不同位置的光源并进行LB光传播模拟求解，构建出灵敏度矩阵W，并得到FDOT重建计算的相关线性方程，具体方程如下：

其中，Φmeas为边界测量值；S为待重建的光场分布。

1.2 ART重建算法

ART可用于线性方程的求解，对式（6）的线性方程WS=Φmeas进行ART求解，具体公式如下：

其中，n是当前光源的索引值；t是迭代索引值；wn代表灵敏度矩阵W的第n行元素；φm是边界测量向量Φmeas的第m个分量；λ'n代表松弛参数。其基本思想是通过当前边界估计值wn[st]T和真实边界测量值Φm的差值Φm-wn[st]T，去校正下一时刻的未知光源分布估计值st+1。通过多次迭代计算，使得估计值st+1能满足上述线性方程，从而实现相应方程的求解。

2 系统设计

FDOT系统主要用于实现光传播模拟和FDOT重建，系统的具体框架如图1所示。对应不同功能，将用户界面划分成两个部分，包括光传播模拟和FDOT重建，并单独进行操作。虽然两个界面所要实现的功能不同，但相应的程序操作流程有着极大的重合性，包括光学参数设置、LB前向模型计算等。因此，对两个不同的界面使用一样的核心处理模块，包括：参数模块、算法模型、结果显示模块和数据交互模块，并进行了详细介绍。

图1 FDOT系统框架图Fig.1 FDOT system framework

2.1 参数模块

参数模块涉及系统中的各种输入参数设置，具体的结构如图2所示。

图2 FDOT系统参数结构图Fig.2 Parameter structure diagram in FDOT system

由图2可见，参数模块的主要工作内容分为四块：完成实验模型参数、光源参数、组织体的光学参数和算法参数。实验模型参数由ExperimentModel类进行管理，包含模型形状及相应的长、宽、高。光源参数由Source类进行管理，包含模拟光源的空间位置。组织体的光学参数由Optical类进行管理，包含光吸收系数、散射系数和g。算法参数由Algorithm类文件进行管理，包括LB前向模型计算时的迭代次数、DmQn模型、空间离散尺度、边界条件和ART重建计算时的松弛参数等。

2.2 算法模块

算法模块与系统所用算法的计算实现相关，包括LB前向模型和ART重建算法的计算流程，算法相关结构如图3所示。

图3 LB和ART算法结构图Fig.3 Algorithm structure diagram of LB and ART

由图3可见，算法模块的主要工作内容分为两个部分：LB前向模型的计算及相应参数的读取和结果存储；ART的重建计算及相应参数的读取及结果存储。LB前向模型的计算对应方法RhoCalculation( )，ART重建算法的计算对应方法ImDataCalculation( )，具体的计算流程如图4所示。

图4 两种算法流程图Fig.4 Two kinds of algorithm flow chart

经LB计算后的结果存储在变量Rho中，经ART重建计算后的结果存储于ImData中。

2.3 结果显示模块

结果显示模块，用于显示计算后的结果，包括光传播模拟时，经LB前向模型计算后的光场分布；FDOT重建时，经LB、ART计算后的重建结果。

观察图5可以知道，结果显示模块的主要功能分为两个部分：对LB计算后的光场分布进行三维、两维和一维的显示（二维、一维显示内容以光源位置为中心），分别对应Show3D( )、Show2D( )和Show1D( )方法；对LB、ART重建计算后的重建结果进行多个维度的显示。

图5 LB和ART重建结果显示结构图Fig.5 Reconstruction result display diagram of LB and ART

2.4 数据交互模块

数据交互模型主要用于完成系统计算结果的数据导出及FDOT重建时的数据导入工作。导出内容包括：光传播模拟时，经LB前向模型计算后的光场分布数据；FDOT重建时，经LB、ART计算后的重建结果数据。导入内容为FDOT的待重建数据，包括FDOT重建时的激发光源位置、探测器位置及探测器上的测量值。

2.5 用户图形界面

用户图形界面分为光传播模拟界面和FDOT重建界面。光传播模拟界面分为参数设置区、操作区、模型显示区和模拟结果显示区，具体如图6所示。在参数设置区设置相应参数，得到如模型显示区所显示的实验模型；在操作区进行LB光传播模拟计算后显示在结果显示区内。

图6 光传播模拟界面Fig.6 Optical propagation simulation interface

FDOT重建的界面设计风格和光传播模拟界面保持一致，也分为参数设置区、操作区、模型显示区和模拟结果显示区，具体如图7所示。操作流程类似于光传播模拟界面的使用流程。

图7 FDOT重建界面Fig.7 FDOT reconstruction interface

通过图形界面的使用，可以使得FDOT系统的操作变得更加直观、简洁，并提升了该系统的可操作性和交互性。

3 系统验证和结果

3.1 光传播模拟的可靠性验证

本研究通过对比MCxyz的MC光传播模拟结果及Toast++的DE光传播计算结果，来验证本系统是否能可靠用于光传播模拟计算。实验模型采用直径为3 cm、高为5 cm的圆柱，在圆柱中心放置初始能量密度为1 nw/sr且持续发光的点光源，如图8所示。以MC的计算结果作为标准。

图8 FDO可靠性模拟实验模型Fig.8 The simulation model of FDOT to test reliability

在实验过程中，实验模型被依次填充不同光学参数的介质，模拟常规和低吸收组织环境，分别对应Case1和Case2，如表2所示。

表2 实验介质的光学参数Tab.2 Optical parameters of mediums in experiment

本系统中的LB前向模型经D3Q26离散化，离散的网格尺寸为0.1 cm，具体的计算流程由章节2.2的图4给出。

3.2 FDOT重建的可靠性验证

通过与待重建标记物的真实位置进行对比，来验证本系统是否能可靠用于FDOT重建。实验以直径为3 cm、高为5 cm的圆柱模型作为成像物体，在成像物体内部放置两个高、直径都为0.2 cm的圆柱形标记物，标记物的中心位置分别为(－0.4 cm,0.4 cm,2.7 cm)、(0.4cm,－0.4 cm,2.7 cm)，并在成像物体周围放置探测器数据采集点。成像物体被填充均匀介质，介质的光学参数μa(r)、μs(r)、g分别为0.3 cm-1、10 cm-1和0.75。

3.3 验证结果分析

不同系统的光传播模拟计算结果由图9给出，显示内容为圆柱实验模型中间虚线圆圈所在层上的光场分布对数值。通过对比可以发现，经本系统LB模型计算的光场分布和经MCxyz系统MC计算的光场分布更加接近，特别是在低吸收组织环境中，两者的结果几乎一致，如图9（a2）～（c2）所示；经Toast++系统DE计算的光场分布要明显亮于其它系统的模型计算结果，图9（b1）～（b2）所示。通过上述实验，验证了本系统用于光传播模拟时的可靠性。

实验模型如图10（a）所示。经本系统计算得到的FDOT重建结果，如图10（b）（c）所示，分别对应三维、二维显示。为了进行更直观的对比，成像标记物的真实位置被放置于重建结果中，如图10（c）中的白色圆圈所示。通过对比可以看出，经本系统计算得到的重建光源位置能和真实标记物的正确位置进行良好重合，无论是在三维显示中，图10（a）（b）所示，还是在二维显示中，图10（c）所示。通过上述实验，验证了本系统用于FDOT重建的可靠性。

4 总结

图9 不同前向模型的二维计算结果对比Fig.9 The comparison of 2D results from different forward models

图10 实验模型和重建结果的三维显示Fig.10 Experimental model and 3D display of reconstruction results

本研究以格子玻尔兹曼前向模型作为核心，搭建了新型的FDOT系统，经光传播模拟实验和FDOT重建实验，评估了所搭建系统的可靠性。主要内容包括系统的核心算法、框架、核心模块及用户可视化界面。对该系统的实用性进行分析，和其它商用系统相比，该系统的开源性和简易性，能极大地降低FDOT成像的使用成本；和其它开源系统相比，由于新型前向模型的使用，使得该系统具有更佳的并行化及速度提升优势。因此，提出的FDOT系统具有很高的推广价值。