基于格子波尔兹曼模型的光在介质中传播的描述方法与流程

本发明涉及一种光在介质中传播的描述方法及其应用，具体地，涉及一种基于格子波尔兹曼模型的光在介质中传播的描述方法及其在荧光扩散断层成像中的应用。

背景技术：

随着光学成像技术的发展，荧光扩散断层成像(FDOT)技术已经受到越来越多的关注。FDOT利用近红外光激发组织体内的荧光团，并根据在组织体边界检测到的光强值重建出荧光团浓度。该技术能够准确定位和量化在体荧光团，可用于肿瘤的早期诊断和分子、细胞活性的监测，具有无创/微创、相对深度深、高灵敏性和特异性等优点。

FDOT技术是对扩散断层成像(DOT)技术的直接扩展，其性能主要依赖于两部分：准确描述光在组织中传播的前向模型和重建算法。重建对于前向模型引入的噪声敏感，且整个成像过程所需时间大部分由前向模型消耗。目前主要的前向模型求解方法有蒙特卡罗(Monte Carlo，MC)法、解析法和数值法。

MC法是20世纪40年代发展起来的一门计算科学，该方法是一种基于统计理论的数值方法。1983年Wilson等人将其用于组织光学领域。MC法将光在组织体中的运动看作一个随机过程，在这个过程中光子或光子包在组织中随机行走从而进行传播，并利用一系列的随机数来模拟光子在介质中的运动和行径，追踪每个光子或光子包的历史，然后对获得的信息加以分析。

MC法是对实际物理过程的真实描述，能对任何复杂几何形状域进行计算，并能获得精确的结果。但蒙特卡罗法耗时非常严重。

解析法又称作格林函数法，各个位置的解可看作是格林函数和光源分布的卷积，从而易于计算。格林函数虽然能够基本求解任意形状下的扩散方程(Diffuse Equation,DE)，但耗时严重。因此，Jorge等人进一步提出了基尔霍夫近似(Kirchhoff approximation，KA)。KA的本质是对格林函数的规范化，它是一种能模拟任意几何体光强分布的线性方法，在求解正问题时，KA不考虑矩阵求逆问题，所以能用来生成系统的敏感函数或者加权值，因此，诸如代数重建技术(algebraic reconstruction techniques,ART)之类的逆问题都可用KA求解。由于KA求解过程中不含矩阵求逆问题，所以计算时间在所述的这些方法中是最少的，但是准确性问题是一个弊端。

数值法主要包括有限元法、有限差分法、有限体积法和边界元法等。这些方法一般都是对近似的DE求解，其中有限元法(Finite Element Method,FEM)最为常用。该方法利用变分原理，建立与DE初值问题等价的积分方程，然后对整个高散射介质进行整体离散以获得等价的有限元方程组，最后求解该方程组获得相应解。还有一些基于有限元的仿真软件也常用于多种光学应用中，如COMSOL等。

光在生物组织中的传播可以近似看作是各向异性扩散，而格子波尔兹曼(Lattice Boltzmann，LB)模型是一种常用的描述各向异性扩散的模型。加之其具有清晰的物理思想、固有的稳定性和并行性，使之非常适合于FDOT技术中前向模型的建立。通过设计LB模型模拟光在组织体中的运动，进而利用重建算法重建出光学图像。因此LB模型方法为实现FDOT的快速高效和准确性提供了实现途径。

技术实现要素：

为了解决以上诸多问题，本发明提出了一种基于格子波尔兹曼模型的光在介质中传播的描述方法，其由于LB模型具有并行性结构，易于实现并行运算，相比于MC法，大大提高了荧光团重建效率；相对于FEM，在具有相当的重建精度基础上，大幅提高了计算效率。

本发明的技术解决方案如下：一种基于格子波尔兹曼模型的光在介质中传播的描述方法，其包括以下步骤：

步骤一，初始化组织体的尺寸、获取光源和检测器的位置，设置组织体的组织参数，包括吸收系数μ_am、散射系数μ_sm和各向异性因子g；

步骤二，荧光团受激发出荧光；

步骤三，设计扩散系数，与距离光源点的距离有关，当距离光源点远时，扩散系数小；反之，当距离光源点近时，扩散系数大，这样能实现光的强度随距离光源点的距离增大而衰减；

步骤四，建立格子波尔兹曼(Lattice Boltzmann，LB)模型的D3Q7模型(表示3维空间7个离散速度方向)，初始化平衡态分布函数；

步骤五，计算LB模型(指“格子波尔兹曼模型”)中的碰撞过程；

步骤六，判断光颗粒是否运动到边界，若是则对边界处的光颗粒进行边界处理，否则直接进入步骤七；

步骤七，计算LB模型中的迁移过程；

步骤八，更新各节点的光强值，即对各节点处强度值求和；

步骤九，判断是否到达最大的迭代次数，如果是则显示输出，否则返回步骤三；

步骤十，更新输出结果。

优选地，所述扩散系数为利用朗伯-比尔定律设计的扩散系数，表示为下式：

其中，L为荧光射入组织体的厚度，该扩散系数能添加与初始荧光光强有关的系数；表示为未加入组织体厚度因素的扩散系数。

对于散射作用远大于吸收作用的生物组织，用约化散射系数(reduced scattering coefficient)μ_sm′＝μ_sm(1-g)代替散射系数μ_sm，则上式能改写为下式：

优选地，所述步骤四中，提出对LB模型演化方程的吸收项的设计，辐射传输方程(RTE)、扩散方程DE、LB方法的宏观方程分别写成以下三式：

其中，式中为组织体吸收和散射造成的辐射率的减小；μ_a(r)为组织体位于空间节点r上的吸收系数；μ_s(r)为组织体位于空间节点r上的散射系数；为位于空间节点r处方向上的源项；式中的μ_aΦ(r，t)+Q(r，t)表示组织体吸收造成的辐射率的减小和源项；因此，对于LB方法，设计吸收项F表示吸收造成的节点处辐射值的减少和源项，表示为：

F＝-μ_af_α(r，t)+Q(r，t)；

其中，μ_a表示为吸收系数，f_α(r，t)表示为在t时刻、节点r处的光颗粒数，Q(r，t)表示为在t时刻、节点r处的源项光子密度。

优选地，所述吸收项保证了能量守恒，由于光在组织体扩散时不断被组织体吸收，该部分光子不参与后续的扩散，本发明设计吸收项以确保能量守恒。

优选地，所述步骤六中，提出采用罗宾(Robin)边界条件，能准确求解LB模型，该边界条件考虑了组织体的折射率n_i大于周围介质的折射率n_t时边界出射光的反射和折射效应，反射概率能够通过反射系数R表示，反射效应和折射效应的边界条件上的辐射率表示为：

其中，r表示位于组织边界上的空间点，时，表示由方向反射到方向上的辐射率；时，表示由方向折射到方向上的辐射率，组织体与周围介质的折射率相差不大时，折射角相对于入射的偏移量不大，并且出射的辐射率不会再对组织体内产生影响。

对于由反射回到组织体内的光颗粒仍然遵守RTE描述，在组织体内继续传播，而折射出组织体的光颗粒能被探测器测量得到，表示为：

在LB方法中，光颗粒始终沿着格线在节点间运动，即与e_α相隔θ_i，即此时，表示为：

光颗粒并非位于节点时，利用LB方法的插值格式处理。

优选地，所述LB模型采用D3Q7模型模拟光子在组织体中的扩散，所述模型包括所有DmQn类模型(表示m维空间n个离散速度方向)，常用的LB模型包括D2Q5、D3Q7、D3Q15、D3Q19等，LB模型的步骤主要包括：

步骤十一，初始化组织体尺寸，设置组织体的参数，包括吸收系数、散射系数和各向异性因子，并获取光源和探测器位置；

步骤十二，使用近红外光激发荧光团发出荧光；

步骤十三，分别使用提出的LB模型模拟激发光和荧光的扩散过程；

步骤十四，根据步骤十三得到的结果，采用重建算法重建组织体内的荧光团浓度。

优选地，所述步骤十一中的组织体能够模拟所有组织的光学参数，所述的荧光团能够设置为任意数量、尺寸大小及位置。

优选地，所述组织体包括第一组织体、第二组织体、第三组织体，组织体是填充有脂肪乳的参数确定的仿体。

优选地，所述扩散系数控制LB模型的扩散速度，本发明对常用扩散系数进行优化，利用朗伯-比尔定律设计出与组织体厚度和荧光光强相关的扩散系数。

优选地，所述组织体是基于FDOT(荧光扩散断层成像)这一应用使用的具体介质。该发明广泛适用于多种介质。所述的模型为一种新型的描述光在介质中传播的模型，并基于该模型构建了新型FDOT系统。所述的方法可广泛应用于诸如光学成像、光声成像、电磁源成像、微波成像、电阻抗成像等图像重建中前向模型的建立和求解。

本发明的积极进步效果在于：本发明由于LB模型具有并行性结构，易于实现并行运算，相比于MC法，大幅提高了荧光团重建效率；相对于FEM，在具有相当的重建精度基础上，大幅提高了计算效率。本发明所述的模型为一种新型的描述光在介质中传播的模型，并基于该模型构建了新型FDOT系统。所述的模型可广泛应用于图像重建中前向模型的建立和求解。

附图说明

图1是本发明基于格子波尔兹曼模型的光在介质中传播的描述方法中基于D3Q7的LB模型的各矢量方向图。

图2是本发明基于格子波尔兹曼模型的光在介质中传播的描述方法的流程示意图。

图3是本发明中基于LB方法的FDOT重建算法流程图。

图4是本发明中Robin边界条件下光颗粒在格子边界的反射与折射图。

图5是本发明仿真模型示意图。

图6a至6h是第一组织体基于本发明和COMSOL方法得到的三维扩散结果图。

图7a至7h是第二组织体基于本发明和COMSOL方法得到的三维扩散结果图。

图8a至8h是第三组织体基于本发明和COMSOL方法得到的三维扩散结果图。

图9a至9d是本发明中基于LB方法得到的三维重建结果图。

具体实施方式

下面结合附图给出本发明实施例，以具体说明本发明的技术方案。

如图2、图4所示，本发明公开了一种基于格子波尔兹曼模型的光在介质中传播的描述方法，其主要包括以下步骤：

步骤一，初始化组织体的尺寸、获取光源和检测器的位置，设置组织体的组织参数，包括吸收系数μ_am、散射系数μ_sm和各向异性因子g；

步骤二，荧光团受激发出荧光；

步骤三，设计扩散系数，与距离光源点的距离有关，当距离光源点远时，扩散系数小；反之，当距离光源点近时，扩散系数大，这样能实现光的强度随距离光源点的距离增大而衰减；

步骤四，建立LB模型的D3Q7模型(表示3维空间7个离散速度方向)，初始化平衡态分布函数；

步骤五，计算LB模型中的碰撞过程；

步骤六，判断光颗粒是否运动到边界，若是则对边界处的光颗粒进行边界处理，否则直接进入步骤七；

步骤七，计算LB模型中的迁移过程；

步骤八，更新各节点的光强值，即对各节点处强度值求和；

步骤九，判断是否到达最大的迭代次数100，如果是则显示输出，否则返回步骤三；

步骤十，更新输出结果。

所述步骤三中的扩散系数为利用朗伯-比尔定律设计的扩散系数，表示为式(1)：

μ_sm表示为荧光波长λ_m下的散射系数，μ_am表示为荧光波长λ_m下的吸收系数，L表示为荧光射入组织体的厚度，g(κ，L)表示为利用朗伯-比尔定律设计的扩散系数。该扩散系数能添加与初始荧光光强有关的系数。

对于散射作用远大于吸收作用的生物组织，用约化散射系数(reduced scattering coefficient)μ_sm′＝μ_sm(1-g)代替散射系数μ_sm，则式(1)能改写为式(2)：

μ_sm表示为荧光波长λ_m下的散射系数，μ_am表示为荧光波长λ_m下的吸收系数，L表示为荧光射入组织体的厚度，g表示为各向异性因子，g(κ，L)表示为利用朗伯-比尔定律设计的扩散系数。

所述步骤四中，提出对LB模型演化方程的吸收项的设计，RTE、DE、LB方法的宏观方程分别写成式(3)、式(4)、式(5)：

式(3)中，c表示为光在组织体中的传播速度，表示为组织体中r点，某个确定方向在某时间点t的辐射率，μ_α(r)表示为该组织体位于空间节点r上的吸收系数；μ_s(r)表示为该组织体位于空间节点r上的散射系数，表示为散射相位函数，表示为位于空间节点r处，某个确定方向上的源项辐射率。

式(3)中表示为组织体吸收和散射造成的辐射率的减小；表示为源项带来的辐射率增加。

式(4)中，c表示为光在组织体中的传播速度，Φ(r，t)表示为节点r处，在某时间点t的光子密度，κ表示为扩散系数，μ_a表示为吸收系数，Q(r，t)表示为节点r处，在某时间点t的光源项的光子密度。

式(5)中，τ表示为松弛时间，e_α表示为光颗粒速度矢量，表示为局部平衡态分布函数。

因此，对于LB方法，设计吸收项F表示吸收造成的节点处辐射值的减少和源项，表示为式(6)：

F＝-μ_af_α(r，t)+Q(r，t) (6)；

μ_a表示为吸收系数，f_α(r，t)表示为在t时刻、节点r处的光颗粒数，Q(r，t)表示为在t时刻、节点r处的源项光子密度。

所述吸收项保证了能量守恒，由于光在组织体扩散时不断被组织体吸收，该部分光子不参与后续的扩散，本发明设计吸收项以确保能量守恒。

所述步骤六中，提出采用Robin边界条件，该边界条件考虑了组织体的折射率n_i大于周围介质的折射率n_t时边界出射光的反射和折射效应，反射概率能够通过反射系数R表示，反射效应和折射效应的边界条件上的辐射率表示为式(7)：

r表示为边界上的空间点，表示为出射光方向，表示为该边界点上的外向单位法向向量。时，表示由方向反射到方向上的辐射率；时，表示由方向折射到方向上的辐射率，组织体与周围介质的折射率相差不大时，折射角相对于入射的偏移量不大，并且出射的辐射率不会再对组织体内产生影响。

对于由反射回到组织体内的光颗粒仍然遵守RTE描述，在组织体内继续传播，而折射出组织体的光颗粒能被探测器测量得到，表示为式(8)：

r表示为边界上的空间点，表示为该边界点上的外向单位法向向量，表示为在某空间点r上，某时刻t时在方向上的辐射率。

在LB方法中，光颗粒始终沿着格线在节点间运动，即与e_α相隔θ_i，即此时，表示为式(9)：

e_α表示为光颗粒速度矢量。

光颗粒并非位于节点时，利用LB方法的插值格式处理。

所述模型通过利用LB模型模拟光子在组织体中的扩散，所述模型包括所有DmQn类模型(表示m维空间n个离散速度方向)，常用的LB模型包括D2Q5、D3Q7、D3Q15、D3Q19等，LB模型的步骤主要包括：

步骤十一，初始化组织体尺寸，设置组织体的参数，包括吸收系数、散射系数和各向异性因子，并获取光源和探测器位置；

步骤十二，使用近红外光激发荧光团发出荧光；

步骤十三，分别使用提出的LB模型模拟激发光和荧光的扩散过程；

步骤十四，根据步骤十三得到的结果，采用重建算法重建组织体内的荧光团浓度。

所述步骤十一中的组织体能够模拟所有组织的光学参数，所述的荧光团能够设置为任意数量、尺寸大小及位置。

所述组织体包括第一组织体、第二组织体、第三组织体，是填充有脂肪乳的参数确定的仿体，本发明中用到3种参数不同的第一组织体、第二组织体、第三组织体。如下表1所示；所述光源为产生发射波长为λ_x激发光的激光器，所述探测器为CCD相机；所述的LB模型：通过利用LB模型模拟光子在组织体中的扩散，常用的LB模型包括D2Q5—二维5个方向扩散，D3Q7—三维7个方向扩散。所述的基于LB模型的FDOT系统在包括PC或服务器在内的硬件计算平台中运行。

表1：圆和圆柱体仿真实验的光学参数设置

所述扩散系数(即松弛因子，与松弛时间τ成反比)控制LB模型的扩散速度，本发明对常用扩散系数进行优化，利用朗伯-比尔定律设计出与组织体厚度和荧光光强相关的扩散系数。

如图5所示，本发明在第一组织体、第二组织体、第三组织体中设计4个case(实验例)，利用LB模型的D3Q7模型模拟光子扩散过程，

如图1所示，LB模型的D3Q7模型是由均匀网格组成，其中网格节点对应于组织体离散后的节点，即光颗粒沿着网格线在相邻的网格点间运动。4个实验例的仿真模型如图5，所述D3Q7模型如图1，包含有7个方向的离散速度，分别表示为式(10)：

α表示为光颗粒运动方向下标，e_αD3Q7表示为光颗粒速度矢量。

D3Q7模型拥有7个离散速度，模型中微观粒子扩散演化方程为式(11)：

f_α(r，t)表示为在t时刻、节点r处的光颗粒数，表示为依赖于局部强度和速度的局部平衡分布函数，ρ为光颗粒密度，e_α表示为光颗粒速度矢量，τ表示为松弛时间，δt表示为时间步长，F_α表示为吸收项。

其中，τ为松弛时间，与前面所述扩散系数的关系为式(12)：

g(κ，L)表示为利用朗伯-比尔定律设计的扩散系数，τ表示为松弛时间。

在D3Q7模型中，A_α表示为式(13)：

α表示为光颗粒运动方向下标，A_α表示为平衡态加权系数。

FDOT的完整过程包含两次DOT的过程，第一次模拟光在组织体中的传播，第二次模拟荧光团发出的荧光在组织体中的传播，两次过程极为相似。因此本发明假设荧光团已经被激发，第一组织体中的三维仿真方法直接采用如图2所示的流程，待扩散完成后，输出结果。如图3所示，FDOT重建算法包括以下步骤:

步骤二十，初始化；

步骤二十一，实验得测量值；

步骤二十二，求解正问题，得值；

步骤二十三，ART求解方程；

步骤二十四，判断条件是否成立，是则转步骤二十五，否则转步骤二十一；

步骤二十五，结束。

第一组织体扩散结果如图6a至6h所示，其重建结果如图9a至9d所示，上述实施的硬件条件为处理器Xeon W3550(英特尔至强W3550处理器)，主频3.07GHz，内存8G；软件配置为Window7系统和matlab2012b。对于3中组织体的四个实验例，仿真所需时间与COMSOL(多物理场仿真软件)结果对比如下表2所示。

表2：LB方法和COMSOL进行三维仿真实验所需时间

第二组织体实施步骤与上述一致，利用上述4种实验例进行基于LB模型的光子扩散，扩散结果如图7a至7h、9a至9d所示，仿真所需时间与COMSOL结果对比如表2所示。

第三组织体实施步骤与上述一致，利用上述4种实验例进行基于LB模型的光子扩散，扩散结果如图8a至8h、9a至9d所示，仿真所需时间与COMSOL结果对比如表2所示。

实验显示，针对FDOT，采用基于LB模型前向模型能够得到较好的重建效果，并且大幅提高了计算效率，该方法有望用到实际应用中。

以上所述的具体实施例，对本发明的解决的技术问题、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。