一种解耦合荧光蒙特卡罗模型的建立方法

摘要

本发明涉及一种解耦合荧光蒙特卡罗模型的建立方法，包括：构建发射密度积分方程用于描述激发光在生物组织中的传输和荧光激发并传输的过程；通过荧光激发和荧光传播过程中的核函数变换，解耦合出荧光激发和传播过程的光学参数与路径概率密度函数的关联，使抽样的荧光光子路径与激发光光子路径一致，同时使荧光光子的权重函数与荧光激发和传播过程的光学参数相关联；从而沿着激发光光子路径，根据荧光权重函数直接计算荧光光子的权重。本发明建立的模型具有精度高、可适用于任意浑浊介质任意荧光分布的特点，在荧光成像、荧光光谱测量与分析方面有及其广泛的应用前景。

说明书

技术领域

本发明涉及数学仿真和生物医学工程领域，具体涉及一种解耦合荧光 蒙特卡罗模型的建立方法。

背景技术

荧光成像技术和荧光光谱测量在生物医学领域中占据重要地位，其广 泛应用于癌症诊断，药物研发以及基因表达可视化的研究。随着荧光探针 技术长足的发展，荧光成像技术开始应用于小动物模型内部特异生物大分 子活动规律的在体跟踪和测量。在外部激发光源的作用下，生物组织内部 荧光团吸收一定波长的光而发出荧光，所探测到的荧光光谱可以反映分子 的功能信息。近年来，荧光分子断层成像技术得到了广泛的关注。

为了提高荧光分子断层成像的计算效率，许多研究小组均聚焦于荧光 蒙特卡罗法加速。Liebert提出了运用于层状浑浊介质中，高效的时间分辨 荧光蒙特卡罗方法[1]。在这种方法中，通过保存大量光子路径信息，建立 了透射光子权重与组织光学参数的解析关系。Jin Chen提出了一种基于计 算时间门控荧光雅克比矩阵的高效荧光蒙特卡罗方法[2]。从激发源到探测 器的背景权重矩阵可由光子路径和荧光吸收系数计算得出。Anand T.N. Kumar提出了一种基于保存光子路径信息的吸收微扰蒙特卡罗方法，适用 于时间分辨多参数的荧光模型[3]。以上三种方法在荧光蒙特卡罗模拟过程 中，均做了一些假定如荧光激发时的各向同性散射由各向异性散射代替 激发光和荧光在组织中随机行走的路径相同。为了荧光分子断层成像更好 地定量，建立一个高精度和高效率的计算模型去描述光子在生物组织中的 传输是非常关键的。以上三种方法在模拟过程中做了较多假定，方法的适 用情况因此受到了一定限制，故无法完全满足荧光分子断层成像重建的需 求。因此发明一种通过荧光激发和荧光传播过程中的核函数变换，解耦合 出荧光激发和传播过程的光学参数与路径概率密度函数的关联，使抽样的 荧光光子路径与激发光光子路径一致，同时使荧光光子的权重函数与荧光 激发和传播过程的光学参数相关联；从而沿着激发光光子路径，根据荧光 权重函数直接计算荧光光子的权重的荧光蒙特卡罗模型，可适用于任意浑 浊介质任意荧光分布的情况，并满足此需求。

[1]Liebert A，Wabnitz H，Obrig H，et al.Non-invasive detection of  fluorescence from exogenous chromophores in the adult human brain[J]. Neuroimage，2006，31(2)：600-608.

[2]Chen J，Venugopal V，Intes X.Monte Carlo based method for  fluorescence tomographic imaging with lifetime multiplexing using time  gates[J].Biomedical optics express，2011，2(4)：871-886.

[3]Kumar A T N.Direct Monte Carlo computation of time-resolved  fluorescence in heterogeneous turbid media[J].Optics letters，2012，37(22)： 4783-4785.

发明内容

本发明目的在于通过荧光激发和荧光传播过程中的核函数变换，解耦 合出荧光激发和传播过程的光学参数与路径概率密度函数的关联，使抽 的荧光光子路径与激发光光子路径一致，同时使荧光光子的权重函数与荧 光激发和传播过程的光学参数相关联；从而沿着激发光光子路径，根据荧 光权重函数直接计算荧光光子的权重。本发明建立的模型计算精度高、且 可适用于任意浑浊介质任意荧光分布的情况。

一种解耦合荧光蒙特卡罗模型的建立方法，其特征在于包括以下步骤：

步骤1：构建粒子输运方程描述稳态下光子在生物组织中的传输过程；

步骤2：构建输运方程等价的发射密度积分方程描述激发光光子在生 物组织中传输的过程和激发光光子被吸收并激发出荧光光子，荧光光子在 生物组织中继续传输的过程；

步骤3：将积分方程的诺伊曼级数解中荧光激发过程和荧光传播过程 中的核函数均用激发光传播过程的核函数表示，从而解耦合出荧光激发和 传播过程的光学参数与路径概率密度函数的关联；

步骤4：投放激发光光子，追踪光子在生物组织中的传输，利用激发 光散射系数抽样，沿着激发光路径，计算探测器上的荧光强度。

步骤1中稳态下光子在生物组织中的传输过程描述为：

$\left(\begin{array}{c}\hat{s}·▿\phi (\overrightarrow{r},\hat{s},v,t)+{\mu }_t(\overrightarrow{r},v,t)\phi (\overrightarrow{r},\hat{s},v,t)=S(\overrightarrow{r},\hat{s},v,t)+\\ {\int }_{4\pi }{\mu }_t(\overrightarrow{r},v,t)\phi (\overrightarrow{r},{\hat{s}}^{\prime },v,t)·C({\hat{s}}^{\prime }\to \hat{s}){\mathrm{d\Omega }}^{\prime }\end{array}\right)$

为光子当前状态的方向；

为光子当前状态的下一个状态的方向；

r为光子当前状态的位置；

v为光子当前状态的频率；

dΩ′为光子当前状态的下一个状态的立体角元；

为在t时刻内，在r点附近单位频率间隔和单位立体角元内 方向在附近的光子数；

为在t时刻内为光子，由外源发出在r点附近单位频率间隔和 单位立体角元内，方向在附近的光子数；

为方向为光子，在r点处进入碰撞，碰撞后方向在方向的 光子概率数；

μ_t为消光系数。

步骤2中激发光传输和荧光激发并传输的过程分别描述为：

$x\left(p\right)=\int x\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})d{p}^{\prime }+S\left(p\right)$

$y\left(p\right)=\int y\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}}){\mathrm{dp}}^{\prime }+\int x\left(p^{\prime }\right)K\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_{\mathrm{af}}){\mathrm{dp}}^{\prime }$

p为光子传输状态，p′为p之后的一个光子传输状态。

x(p)为激发光的发射密度、y(p)为荧光的发射密度。

S(p)为源发射密度。

K为光传播过程的核函数。

为激发光的散射系数、为激发光的吸收系数。

为荧光的散射系数、为荧光的吸收系数。

μ_af为荧光团的吸收系数。

步骤3中解耦合出荧光光学参数(激发和传播过程)与路径概率密度 函数的关联，荧光激发过程和荧光传播过程中的核函数均可用激发光传播 过程的核函数表示：

(1)荧光激发过程的核函数用激发光传播过程的核函数表示为：

$\left(\begin{array}{c}K(p^{\prime }\to p,{\mu }_{\mathrm{af}})=T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})C({\hat{s}}^{\prime }\to \hat{s}|\hat{r},{\mu }_{\mathrm{af}})\\ =K({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\times \\ \frac{{\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)+{\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right)}{{\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\frac{C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_{\mathrm{af}})}{C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})}\\ =K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\times \\ \frac{\eta {\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\frac{P_I({\hat{s}}^{\prime }·\hat{s})}{P_A({\hat{s}}^{\prime }·\hat{s})}\end{array}\right)$

r′为光子当前状态的下一个状态的位置；

为从r′处发出的光子，在r处进入碰撞的密度。

为激发光的消光系数。

l为光子每步的步长。

η为量子效率。

P_I为各向同性相位函数、P_A为各向异性相位函数。

(2)荧光传播过程的核函数用激发光传播过程的核函数表示为：

$\left(\begin{array}{c}K(p^{\prime }\to p,{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})=T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})\\ =T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})\times \\ C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\frac{{\mu }_s^{\mathrm{em}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\\ =\frac{{\mu }_s^{\mathrm{em}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\end{array}\right)$

为荧光的消光系数。

步骤4中抽样的荧光光子路径与激发光光子路径一致，荧光光子的权 重函数与荧光光学参数(激发和传播过程)相关联；沿着激发光光子路径， 根据荧光权重函数直接计算荧光光子的权重。

步骤4中探测器上的荧光强度表示为：

$\left(\begin{array}{c}D={\int y}_m\left(p\right)g\left(p\right)\mathrm{dp}\\ =\underset{m=0}{\overset{\infty }{\Sigma }}\int ...\int {\tau }_m\left(p\right)w_m\left(p\right)\underset{k=0}{\overset{m}{\Pi }}{\mathrm{dp}}_k\mathrm{dp}\end{array}\right)$

D为探测上接收到的荧光光子强度。

m为到达探测器前光子所碰撞的次数。

g(p)为探测器上探测到的荧光统计量的函数。

τ_m(p)为光子的概率密度。

w_m(p)为光子的权重函数。

(1)光子的概率密度τ_m(p)表示为：

${\tau }_m\left(p\right)=S\left(p_0\right)K(p_0\to p_1,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})...K(p_{m-1}\to p_m,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})$

其中S(p)＝S(p₀)，S(p₀)为源发射密度；

(2)光子的权重函数w_m(p)表示为：

$\left(\begin{array}{c}{w}_m\left(p\right)=g\left(p\right)\underset{i=0}{\overset{m}{\Sigma }}\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\frac{{\mathrm{\eta \mu }}_{\mathrm{af}}\left({\overrightarrow{r}}_i\right)}{{\mu }_s^{\mathrm{ex}}\left({\overrightarrow{r}}_i\right)}\\ \times \frac{P_I({\hat{s}}_i·{\hat{s}}_{i+1})}{P_A({\hat{s}}_i·{\hat{s}}_{i+1})}\underset{j=i+1}{\overset{m}{\Pi }}\left(\frac{{\mu }_s^{\mathrm{em}}\left({\overrightarrow{r}}_j\right)}{{\mu }_s^{\mathrm{ex}}\left({\overrightarrow{r}}_j\right)}\right)\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})\end{array}\right)$

附图说明

图1为本发明的基本流程图。

图2为本发明的具体构建的圆柱模型图。

图3a为在标准荧光蒙特卡罗(sfMC)法模拟中得到的探测器上归一化荧 光强度分布图。

图3b为解耦合荧光蒙特卡罗(dfMC)法模拟中得到的探测器上归一化 荧光强度分布图。

图4为标准和解耦合荧光蒙特卡罗法模拟中得到的探测器上归一化 光强度分布图像的等高线比较图。

具体实施方式

结合附图对本发明作进一步的描述。

如图1所示，本发明的实施步骤如下：

(1)构建粒子输运方程描述稳态下光子在生物组织中的传输过程；

$\left(\begin{array}{c}\hat{s}·▿\phi (\overrightarrow{r},\hat{s},v,t)+{\mu }_t(\overrightarrow{r},v,t)\phi (\overrightarrow{r},\hat{s},v,t)=S(\overrightarrow{r},\hat{s},v,t)+\\ {\int }_{4\pi }{\mu }_t(\overrightarrow{r},v,t)\phi (\overrightarrow{r},{\hat{s}}^{\prime },v,t)·C({\hat{s}}^{\prime }\to \hat{s}){\mathrm{d\Omega }}^{\prime }\end{array}\right)$

(2)构建输运方程等价的发射密度积分方程描述激发光光子在生物组织中 传输，激发光光子被吸收并激发出荧光光子，荧光光子继续在生物组织中 传输最终被探测器接收到的过程；

激发光光子在生物组织中传输等价的发射密度积分方程描述为：

$x\left(p\right)=\int x\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})d{p}^{\prime }+S\left(p\right)$

激发光光子被吸收并激发出荧光光子，荧光光子在生物组织中继续传 输的过程描述为：

$y\left(p\right)=\int y\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}}){\mathrm{dp}}^{\prime }+\int x\left(p^{\prime }\right)K\left(p^{\prime }\right)K(p^{\prime }\to p,{\mu }_{\mathrm{af}}){\mathrm{dp}}^{\prime }$

(3)将积分方程的诺伊曼级数解中荧光激发过程和荧光传播过程中的核函 数均用激发光传播过程的核函数表示，建立荧光激发过程和荧光传播过程 与激发光传播过程的联系，从而解耦合出荧光激发和传播过程的光学参数 与路径概率密度函数的关联；

核函数K表示从状态p到p′的转换，可将其表示为碰撞核和迁移核的 乘积：

$K(p^{\prime }\to p,{\mu }_s,{\mu }_a)=T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s,{\mu }_a)C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s,{\mu }_a)$

T为迁移核，表示从r′处发出的光子，在r处进入碰撞的密度。

$T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s,{\mu }_a)={\mu }_t\left(\overrightarrow{r}\right)\exp (-{\int \mu }_t({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\delta (\hat{s}-(\overrightarrow{r}-+{\overrightarrow{r}}^{\prime })/|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }\left|\right)/{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}^2$

C为碰撞核，表示方向为光子，在r点处进入碰撞，碰撞后方向在方 向的dΩ立体角元内的光子概率数。

激发光传播过程的碰撞核表示为：

$C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})={\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)P_A({\hat{s}}^{\prime }·\hat{s})/({\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)+{\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right))$

荧光激发过程的碰撞核表示为：

$C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_{\mathrm{af}})={\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right)P_A({\hat{s}}^{\prime }·\hat{s})/({\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)+{\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right))$

荧光传播过程的碰撞核表示为：

$C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})={\mu }_s^{\mathrm{em}}\left(\overrightarrow{r}\right)P_A({\hat{s}}^{\prime }·\hat{s})/({\mu }_t^{\mathrm{em}}\left(\overrightarrow{r}\right)+{\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right))$

在非荧光区域式中为0。

荧光激发过程的核函数用激发光传播过程的核函数表示为：

$\left(\begin{array}{c}K(p^{\prime }\to p,{\mu }_{\mathrm{af}})=T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})C({\hat{s}}^{\prime }\to \hat{s}|\hat{r},{\mu }_{\mathrm{af}})\\ =K({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\times \\ \frac{{\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)+{\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right)}{{\mu }_t^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\frac{C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_{\mathrm{af}})}{C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})}\\ =K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\times \\ \frac{\eta {\mu }_{\mathrm{af}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\frac{P_I({\hat{s}}^{\prime }·\hat{s})}{P_A({\hat{s}}^{\prime }·\hat{s})}\end{array}\right)$

荧光传播过程的核函数用激发光传播过程的核函数表示为：

$\left(\begin{array}{c}K(p^{\prime }\to p,{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})=T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})\\ =T({\overrightarrow{r}}^{\prime }\to \overrightarrow{r}|{\hat{s}}^{\prime },{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})\times \\ C({\hat{s}}^{\prime }\to \hat{s}|\overrightarrow{r},{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\frac{{\mu }_s^{\mathrm{em}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\\ =\frac{{\mu }_s^{\mathrm{em}}\left(\overrightarrow{r}\right)}{{\mu }_s^{\mathrm{ex}}\left(\overrightarrow{r}\right)}\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})K(p^{\prime }\to p,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})\end{array}\right)$

(4)投放激发光光子，追踪光子在生物组织中的传输，利用散射系数抽样 沿着激发光路径，计算探测器上的荧光强度；

追踪光子在生物组织中的传输过程具体步骤为：

(4.1)投放激发光光子，激发光光子初始权重设置为1。

(4.2)始终利用激发光的散射系数抽样，光子每步的步长为 ，散射方向由Henyey-Greenstein函数确定。

(4.3)光子包沿着激发光的散射路径连续衰减。

(4.4)激发光光子在非荧光区域，光子包权重的衰减比例为激发光光子在荧光区域，光子包权重的衰减比例为

(4.5)依据荧光激发过程的核函数与激发光传播过程的核函数间关系， 激发光光子散射的概率与激发的荧光光子散射的概率比例为荧光 光子激发后，依旧沿着原激发光散射方向移动。

(4.6)如果荧光的光学参数(吸收系数和散射系数)与激发光的光学 参数不一致，依据荧光激发过程的核函数与激发光传播过程的核函数间关 系，荧光光子在非荧光区域内光子包权重的衰减与激发光子包权重衰减的 比例为荧光光子在荧光区域光子包权重的衰减 与激发光子包权重衰减的比例为$\exp \left((-{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})l\right)\exp \left(({\mu }_t^{\mathrm{em}}-{\mu }_t^{\mathrm{ex}})l\right)\frac{{\mu }_s^{\mathrm{em}}}{{\mu }_s^{\mathrm{ex}}}.$

(4.7)沿着激发光路径，计算探测器上的荧光强度。经过m次碰撞后， 在p点的发射密度函数为y_m(p)：

$\left(\begin{array}{c}{y}_m\left(p\right)=\underset{m=0}{\overset{\infty }{\Sigma }}\underset{i=0}{\overset{m}{\Sigma }}\int ...\int S\left(p_0\right)K(p_0\to p_1,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})...\\ \times K(p_{i-1}\to p_i,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}}+{\mu }_{\mathrm{af}})\times K^{\mathrm{xm}}(p_i\to p_{i+1},{\mu }_{\mathrm{af}})\times K(p_{i+1}\to p_{i+2},{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})...\\ K(p_{m-1}\to p_m,{\mu }_s^{\mathrm{em}},{\mu }_a^{\mathrm{em}})\underset{k=0}{\overset{m}{\Pi }}{\mathrm{dp}}_k\end{array}\right)$

与此同时，可得到光子的概率密度τ_m(p)表示为：

${\tau }_m\left(p\right)=S\left(p_0\right)K(p_0\to p_1,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})...K(p_{m-1}\to p_m,{\mu }_s^{\mathrm{ex}},{\mu }_a^{\mathrm{ex}})$

光子的权重函数w_m(p)表示为：

$\left(\begin{array}{c}{w}_m\left(p\right)=g\left(p\right)\underset{i=0}{\overset{m}{\Sigma }}\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\mu }_{\mathrm{af}}({\overrightarrow{r}}^{\prime }+l\hat{s})\mathrm{dl})\frac{{\mathrm{\eta \mu }}_{\mathrm{af}}\left({\overrightarrow{r}}_i\right)}{{\mu }_s^{\mathrm{ex}}\left({\overrightarrow{r}}_i\right)}\\ \times \frac{P_I({\hat{s}}_i·{\hat{s}}_{i+1})}{P_A({\hat{s}}_i·{\hat{s}}_{i+1})}\underset{j=i+1}{\overset{m}{\Pi }}\left(\frac{{\mu }_s^{\mathrm{em}}\left({\overrightarrow{r}}_j\right)}{{\mu }_s^{\mathrm{ex}}\left({\overrightarrow{r}}_j\right)}\right)\exp (-{\int }_0^{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}({\mu }_t^{\mathrm{em}}({\overrightarrow{r}}^{\prime }+l\hat{s})-{\mu }_t^{\mathrm{ex}}({\overrightarrow{r}}^{\prime }+l\hat{s}))\mathrm{dl})\end{array}\right)$

g(p)为探测器上探测到的荧光统计量的函数。探测器上的荧光强度可 以计算得出：

$\left(\begin{array}{c}D={\int y}_m\left(p\right)g\left(p\right)\mathrm{dp}\\ =\underset{m=0}{\overset{\infty }{\Sigma }}\int ...\int {\tau }_m\left(p\right)w_m\left(p\right)\underset{k=0}{\overset{m}{\Pi }}{\mathrm{dp}}_k\mathrm{dp}\end{array}\right)$

下面通过实例进一步阐述本发明。

实施例：

模拟中所构建的圆柱模型如图2所示，该圆柱模型充满1％的脂肪乳溶 液，充满荧光染料的玻璃棒插在脂肪乳溶液中。对圆柱模型进行体素剖分， 剖分后得到的体素模型中共含有301401个体素，每个体素的尺寸为 0.05cm，模型的尺寸为3.05cm×3.05cm×4.05cm。光源的位置取在图示位 置，探测器取在源对面180度的区域内，其均匀分布在80层，每层180个 探测器。在表1中我们列出了激发光和荧光波段下的组织光学参数值，其 余参数如g设为0.9左右，折射率设为1.37，均为近红外光谱下生物组 的典型值。模拟的光子数为10⁹。

表1 激发光和荧光波段下的组织光学参数值

图3a为在标准荧光蒙特卡罗(sfMC)法模拟中得到的探测器上归一化荧光强 度分布图，图3b为解耦合荧光蒙特卡罗(dfMC)法模拟中得到的探测器上 归一化荧光强度分布图，图4为标准和解耦合荧光蒙特卡罗法模拟中得到 的探测器上归一化荧光强度分布图像的等高线比较。从探测器上归一化荧 光强度分布图和等高线图可以看出两种方法结果符合得很好。