快速荧光扩散断层成像的格子玻尔兹曼前向模型

岑星星，严壮志,

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

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

0 引言

荧光扩散断层成像（Fluorescence Diffuse Optical Tomography,FDOT）是一种新型的医学成像模态，通过体表采集的荧光信号，完成体内分子探针的三维成像功能，能在细胞、分子水平实现相应病灶的定位及监测工作[1]。FDOT所具有的高灵敏度、安全无创、低成本等优势，使得该成像技术在近几年得到迅速发展，并在临床、预临床等医学领域取得广泛应用，如：癌症检查、药物研究等[2]。对于FDOT的发展，关于前向模型和重建算法的研究必不可少，前向模型描述了光在特定介质中的传播过程；重建算法则通过前向模型和边界采集的光学信号，对体内标记物的光场分布进行重建计算[3]。由于生物组织的高散射性和探测器的测量不充分，相应的重建计算属于不适定问题，使得所用前向模型的微小偏离会造成重建结果的巨大误差。因此，建立一个准确且计算可行的前向模型对于FDOT的实现非常重要。

辐射传输方程（Radiative Transfer Equation,RTE）和扩散方程（Diffusion Equation，DE）是FDOT中的两个经典前向模型。RTE是一种被广泛接受的高精度模型，但受制于昂贵的计算成本，使其仅出现于简单的FDOT运用场景中[4]。相较于RTE，RTE的简化模型有着更好的适用性，其中的DE作为RTE的一阶球谐近似，具有更高的计算效率，并在FDOT的相关研究中有着非常广泛的应用。然而，DE也存在着一系列的问题，无法对靠近光源和边界区域的光传播过程进行准确描述[5]。此外，DE经有限元等数值方法求解的计算效率依然无法满足我们的需求，尤其是在增加计算网格、探测器数量以保证FDOT成像质量的情况下。

格子玻尔兹曼（Lattice Boltzmann,LB）方法具有物理含义清晰、计算程序简单、计算结构并行性高等优点[6]，开发基于LB的前向模型，对于FDOT的性能提升有着巨大潜能。LB方法源于格子气自动机（LGA）方法，是一种特殊的离散分子动力学模型，已成功应用于各种物理学研究，包括流体力学、热力学等[7]。近年来，LB方法也被发展用于辐射传输描述，涉及源项、介质和边界上的辐射传输描述[8-10]。研究表明，用于辐射描述的LB能提供精确的计算结果，且有着较高的计算效率。然而，当前所涉及的大多数LB辐射传输研究都是处于封闭长方体容器中的热辐射问题，极少有完整关于FDOT前向模型的研究报道。

基于我们最初的想法[11]，本文提出了基于LB的FDOT前向模型，通过RTE离散向原有LB方程引入光传播描述[12-13]。进一步，本文给还出了基于LB前向模型的FDOT算法流程，选择代数重建技术（Algebra Reconstruction Technique，ART）作为重建算法，寻找成像标记物的相应光源分布，将其代入LB前向模型，使计算得到的模拟测量能与实际的测量达到最佳拟合。

1 基于格子玻尔兹曼方法的前向模型

1.1 光在生物组织内的传播模型

LB前向模型的光传播描述，通过RTE引入。类比于传统LB模型的离散化推导，得到经RTE离散后的LB前向模型。

传统的LB基本模型可由玻尔兹曼方程经时间、速度、相空间的离散化得到[13]。玻尔兹曼方程描述了粒子在介观尺度上的统计行为，相应的Bhatnagar-Gross-Krook（BGK）近似如下所示[14]：

其中，f=f(r,ξ,t)为分子密度函数，代表t时刻，位置r处，速度为ξ的分子数量；λ为松弛时间，代表f趋近于平衡态feq的速度快慢。假设离散的时间间隔为△t，离散后的空间为Ωk，离散速度模型为DmQn（m代表空间维数，n代表离散速度的数量），经玻尔兹曼BGK离散得到的LB基本模型如下所示[13]：

其中，f(r,ξ,t)为离散后的分子密度函数，代表t时刻，位置r∈Ωk处，以离散速度ξi进行运动的分子数量（下标i为离散后速度方向的索引，介于1到n之间）；f(r+ξ△t,ξ,t+△t)为下一时刻t+△t上的分子密度函数；τ≡λ/△t为无量纲的松弛因子；feq(r,ξi,t)为离散后的平衡态分布函数。

同时，RTE可以改写为玻尔兹曼方程的类似形式，公式及物理描述的相似性提供了条件支持。对于物理空间Ω中的单色光传播，RTE可给出精确描述[4]。进一步，RTR可重写为如下形式：

其中，w(r)=μs(r)/β(r)代表反照率；为时刻沿方向运动的光源项；为散射项函数，描述了光子经散射后由方向运动到方向的可能性。对比式（1）和式（3）可以看出，RTE和玻尔兹曼BGK具有相似的数学形式，并且前者的对应后者的平衡态分布函数feq，cβ(r)对应松弛时间λ的倒数。此外，对两者描述的物理量进行观察，发现RTE的辐射率φ和玻尔兹曼BGK的分子密度函数f，都是关于时间、物理空间和速度的函数。虽然，辐射率φ从能量角度进行描述，不同于分子密度函数f中的粒子数量，但光传播过程中的辐射率φ和相应的光子密度分布函数有着直接联系[12]：

其中，h为普朗克常量；v为光的频率大小。

通过上述过程建立了RTE和玻尔兹曼BGK的联系。接下来本文类比于LB基本模型的离散化推导，通过RTE离散得到相应的LB前向模型。对应式（3）、式（4），RTE经时间、速度和相空间离散后的LB形式如下所示：

其中，Θ(si,sj)为离散后的散射相函数；wi为权重系数。下标i、j指向不同的离散速度方向，当取D3Q6离散速度模型（图1）时，为对应箭头方向，且介于1到6之间。

图1 LB方法中D3Q6离散模型Fig.1 The D3Q6 for discrete model in the LB

1.2 格子玻尔兹曼前向模型的建立

1.2.1 平衡态权重系数权重系数决定了LB前向模型中平衡分布函数S(r,t)的具体表达，本文通过局部辐射场的物理守衡进行权重系数求解。假设辐射率在达到平衡态时保持各向同性，且大小为(r,t)，同时假设光在均匀介质中以恒定速度大小c进行传播，k阶辐射强度的矩方程经过高斯求积后可近似为：

对式（8）～（9）进行整理，可以得到如下矩阵方程：

选择DmQn并确定离散方向后，相应的权重系数通过式（10）求解得到，不同离散模型对应的权重系数如下所示：

1.2.2 平衡态的散射相函数

1.2.3 边界条件

为了能让LB前向模型的边界计算变得简便，同时又能维持较好的精度，在边界节点r∈∂Ωk（k为离散空间节点的索引）进行如下光子密度φ(r,t)的计算：

并与DE中的Robin边界条件进行关联，得到如下边界计算公式[18]：

其中，D(r)1/[3(μa(r)+(1-g)μs(r))]为光扩散系数；是边界处的外法线向量；当组织周围填充空气，组织折射率为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。

2 荧光分子断层成像实现

2.1 LB前向模型计算流程

本文使用迭代算法完成LB前向模型的求解。假设模拟光源为稳定光源，通过迭代计算直到输出结果稳定为止。LB前向模型的具体计算流程如下：

（1）选择离散速度模型DmQn，并用适当数量的格子对连续空间Ω进行网格划分；

（2）设置初始参数μa(r)、μs(r)、g和q(r,t)，并初始化φ(r,t)；

（3）根据式（10）和（12），分别计算权重系数和离散后的散射相函数；

根据新经济地理学和空间集聚等理论，在产业集聚过程中，由于生产要素和知识技术的集聚使得资源得到合理利用，基础设施得到共享，进而使边际成本降低；随着大量产业集聚于同一区域，由于能源、劳动力等是有限的，过度集聚会造成边际成本上升。因此，假设产业集聚和边际成本可能具有非线性关系：

（4）使用式（7），计算平衡态分布函数S(r,t)；

（6）根据式（13）和（14），在边界节点使用罗宾（Robin）边界条件进行计算；

（7）通过式（15）计算相对差异，如最大差异值低于阈值（本文取10-2），则终止迭代过程，否则返回步骤4；

（8）计算格子节点r∈Ωk处的光子密度φ(r,t)：

2.2 基于LB前向模型的FDOT重建

荧光分子断层成像方程由下式给出：

（1）初始化待重建光场分布S0；

（2）采集边界测量值φmeas；

（3）通过LB前向模型计算生成灵敏度矩阵W；

（4）计算当前估计值WS'和真实测量值φmeas的差值φmeas-WS'（t代表计算时刻）；

（5）用当前差值矫正下一时刻的光场分布：

（6）重复步骤（4）和（5），使得相应差值φmeas-WSt足够小，从而实现待重建光场分布的求解。

3 实验与结果分析

本文分别进行了数值仿体和物理仿体实验，通过对比不同重建结果，评估了LB前向模型用于FDOT成像时的性能表现。

数值仿体实验以直径3 cm、高5 cm的圆柱作为成像物体，两个高、宽都为0.2 cm的圆柱标记物被放置于成像物体中，使用不同位置的标记物，如表1所示，完成不同实验模型的设置，具体如图2(a)所示。为简化问题，圆柱体被填充均匀介质，相应的光学参数μa(r)、μs(r)、g分别为0.3 cm-1、10 cm-1和0.75。

在实验计算中，LB模型被用于光传播模拟，通过D3Q6进行速度离散，并将物理空间离散为15×15×25的节点网格，在离散节点进行LB迭代计算，并在最大差异值小于10-2时终止，具体计算流程见2.1。其中，实验以DE作为对比，进行相似程度的离散，对应5 831个节点，并在离散节点使用限元方法进行求解，相应计算在Toast++上完成[19]。

表1 数值仿体实验及标记物位置Tab.1 Numerical simulation experiment and targets positions

图2 数值仿体模型和不同前向模型的重建结果Fig.2 Numerical phantom model and reconstruction results from different forward models

对应不同的实验组，图2给出了基于LB和基于DE前向模型的FDOT重建结果，左右两边分别对应重建结果的三维、二维显示，二维结果对应三维空间光源所在层（对应a-c的黑色圆圈）上的重建光场分布。通过对比可以看出，基于LB前向模型的重建结果相似于基于DE模型的重建结果，和正确标记物位置都有着良好重叠；相对而言，基于LB前向模型的FDOT重建结果在正确区域上显示得更加饱满，而基于DE的重建结果有着更少的伪影表现。

更多的细节对比由图3给出，对应图2(e)和2(f)白色虚线上的重建结果。通过对比可以看出，基于LB前向模型的重建结果能和标记物的正确位置有着更好的重叠，非重叠区域的面积更小。另外，在中间非标记物的区域，基于LB前向模型的重建结果值为0，说明相应的重建能将两个标记物的位置进行完全分离。其次，实验还对比了不同前向模型重建结果的对比信噪比（Contrast-to-noise ratio，CNR），将真实标记存在和不存在空间分别定义为感兴趣区域（ROI）和背景区域（BCK），相应的计算公式如下：

图3 不同前向模型对应重建结果的一维对比Fig.3 The 1D comparison of reconstruction results from different forward models

其中，μROI和μBCK分别为ROI和BCK区域的平均重建强度；σ2ROI和σ2BCK分别为相应区域重建强度的方差；wROI和wBCK分别代表相应区域的容积大小。另外，实验还记录了不同前向模型的FDOT计算时间。上述实验的CNR值和计算时间由表2给出。通过观察可以发现，LB和DE模型所对应重建结果的CNR值非常接近；基于LB的FDOT重建比基于DE的FDOT重建快上将近5倍。

表2 不同前向模型重建结果的CNR值和时间Tab.2 The CNR and calculation time from different models

本文除数值仿体实验外，还进行了物理仿体实验。使用FDOT/XCT混合的成像系统，完成成像物体表面光学信号的采集。实验以直径为3 cm、高为4.4cm的圆柱为成像物体，圆柱体内填充均匀介质，相应的光学参数μ a(r)、μ s(r)、g分别为0.3 cm-1、10 cm-1和0.75，一个填充了荧光标记物（ICG）、直径为0.4 cm的管子被放置于仿体模型中，如图4（a）所示。

在实验计算过程中，使用了类似于数值仿体实验的LB、DE计算。相应的重建结果由图4给出，左右两边分别对应重建结果的三维、二维显示，其中的二维显示对应三维圆柱黑色圆圈所在层的内容。通过三维显示对比可以看出，基于LB的重建结果呈现出了类似于真实标记物形状的圆柱形；而基于DE的重建结果呈现出了和真实标记物形状差异较大的球型；稍显不足的是，基于LB的FDOT重建结果有着较高的伪影。通过重建结果的二维对比，可以看出基于LB的重建结果在正确区域显得更亮，能和真实的重建区域进行更好拟合。

图4 物理仿体模型和不同前向模型的重建结果Fig.4 Physical phantom model and reconstruction results from different forward models

更多的细节对比由图5给出，分别对应图4（e）和2（f）白色虚线上的重建结果。通过对比可以看出，基于LB和DE的FDOT重建结果都能和标记物所在位置的轮廓进行很好重合。

图5 不同前向模型的一维重建结果对比Fig.5 The 1D comparison of reconstruction results from different forward models

其次，实验还记录了上述物理仿体实验中不同前向模型FDOT重建结果的CNR值和相应的计算时间，由表3给出。经对比表发现，LB前向模型的使用可以使得FDOT的重建速度提升5倍左右，CNR提升2倍左右。

表3 不同前向模型重建结果的CNR值和时间Tab.3 The CNR and calculation time from different models

4 结论

FDOT作为一种新型成像技术，具有很高的应用价值。其中的前向模型对于FDOT的成像性能有着直接影响。为了进一步提升FDOT的性能，有必要建立一个计算快速且准确的前向模型。本文，我们提出了一种基于LB方法的前向模型，并将其用于FDOT成像。本文通过数字、物理仿体实验，验证了LB前向模型对于FDOT的性能提升。实验证明，LB前向模型的使用可以让FDOT的计算时间得到大量缩短，且具有较好的精度表现。