一种混浊介质中光子传输方法及装置

摘要

本发明涉及一种混浊介质中光子传输方法及装置，在混浊介质中布置多个离散的场节点，求解混浊介质中每个场节点上的移动最小二乘形函数；利用完全变换法修正移动最小二乘形函数，使其具有狄拉克函数特性；利用修正的移动最小二乘形函数进行插值得到混浊介质的近似光通量密度，并利用迦辽金方法获得扩散方程和第三类边界条件的弱解形式，得到矩阵阵方程并进行求解，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。本发明仅布置场节点，而不需要任何场节点连接信息和单元信息，避免复杂的网格剖分过程，对移动最小二乘近似形函数进行修正，使得形函数具有狄拉克函数性质，可直接添加任何边界条件。

权利要求

1.一种混浊介质中光子传输的方法，其特征在于：

步骤1：在混浊介质中布置多个离散的场节点和用于数值积分的背景网格，用于确定混浊介质中每个组成部分的光学特性参数后，在混浊介质中布置一系列的场节点，用移动最小二乘近似求解得到混浊介质中每个场节点上的移动最小二乘形函数；

步骤2：确定背景网格内的高斯积分点：在每个高斯积分点支持域内搜索用于插值的节点，并计算插值节点处的移动最小二乘近似形函数；利用完全变换法对计算得到的移动最小二乘近似形函数进行修正，使其满足狄拉克函数特性，在背景网格内选取高斯积分点，对混浊介质边界和整个区域进行高斯数值积分，用于直接准确的施加任何边界条件；

步骤3：利用修正的移动最小二乘形函数对混浊介质的光通量密度Φ(x)进行插值得h到混浊介质的近似光通量密度Φ(x)，并利用迦辽金方法获得用于描述光子传输的扩散方程和第三类边界条件的弱解形式；所述弱解形式为：

1 1

其中，Ψ(x)∈H(Ω)，H(Ω)为索伯列夫空间；D(x)是与位置相关的光子扩散系数，Φ(x)是在点x处的光通量密度，Ψ(x)是测试函数，μa(x)表示吸收系数，An(x)是一个表示混浊介质和周边介质折射系数不匹配的常数，S(x)表示混浊介质内的光源密度，Ω和分别是混浊介质区域及其边界；

h

步骤4：将混浊介质的近似光通量密度Φ(x)代入扩散方程和第三类边界条件的弱解形式，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。

2.根据权利要求1所述的混浊介质中光子传输方法，其特征在于：将所述混浊介质中每个组成部分的光学特性参数作为先验知识，表现为混浊介质中各个组成部分光学特性参数的个体差异，用于解决复杂混浊介质的非匀质特性。

说明书

一种混浊介质中光子传输方法及装置

技术领域

[0001] 本发明属于分子影像领域，涉及混浊介质中光子传输方法，尤其是涉及一种基于修正无网格迦辽金的光子传输方法。

背景技术

[0002] 光学分子影像是近年来新兴的一种分子影像模态，光学分子影像主要包括自发荧光和激发荧光两种方式。自发荧光是利用荧光素酶标记混浊介质中的目标，在ATP以及氧气存在的情况下，若荧光素酶遇到底物荧光素，荧光素酶将会催化荧光素的氧化反应并产生光子，产生光的强度与标记目标的数量成正比。激发荧光技术则采用荧光报告基团标记目标，荧光报告基团吸收外界激发光源射入混浊介质中的光子，产生能量跃迁至受激虚态，从受激虚态向下跃迁至基态过程中，产生荧光光子。

[0003] 由于混浊介质具有强散射特性，光子在其内部的传输不再沿直线传播，而是经过大量的散射过程，从而远远偏离了原来的运行轨道。因此，光子在混浊介质中的传输是一个非常复杂的过程。在光学分子影像领域，有限元算法是一种经典的解决光子在复杂混浊介质中传播的数值计算方法。然而，在使用有限元算法之前，必须对混浊介质进行网格剖分，并对网格数据进行预处理。虽然目前已研发了多种网格剖分软件，但是对具有复杂内部结构的非规则混浊介质来说，进行网格剖分还是比较困难和耗时的。另外，由于混浊介质的边界并非是一个固定的边界，而是一个随时间移动的边界，因此在光子传输问题求解过程中，也要随时间对混浊介质进行有限元网格剖分。此外，在有限元自适应技术实现过程中，在两次迭代之间需要对混浊介质进行网格重构、调整网格节点间连接信息，因此实现起来也是比较复杂耗时的。

发明内容

[0004] 为了解决光子传输有限元分析存在的网格剖分困难，本发明的目的是基于修正无网格迦辽金法求解光子在非匀质复杂混浊介质中的传输，此外，为了解决现有无网格方法不能直接准确施加第三类边界条件的问题，本发明对移动最小二乘近似形函数进行了修正，使得形函数具有狄拉克函数性质，可以直接准确地施加任何边界条件。

[0005] 为了实现所述的目的，本发明第一方面，提供一种基于修正无网格迦辽金的光子传输方法的技术方案包括如下步骤：

[0006] 步骤1：在混浊介质中布置多个离散的场节点和用于数值积分的背景网格，用于确定混浊介质中每个组成部分的光学特性参数后，来求解混浊介质中每个场节点上的移动最小二乘形函数；

[0007] 步骤2：利用完全变换法修正移动最小二乘形函数，使其满足狄拉克函数特性具有狄拉克(Delta)函数特性，用于直接准确的施加任何边界条件；

[0008] 步骤3：利用修正的移动最小二乘形函数对混浊介质的光通量密度Φ(x)进行插h值得到混浊介质的近似光通量密度Φ(x)，并利用迦辽金方法获得用于描述光子传输的扩散方程和第三类(Robin)边界条件的弱解形式；

[0009] 步骤4：将混浊介质的近似光通量密度Φh(x)代入扩散方程和第三类边界条件的弱解形式，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。

[0010] 为了实现所述的目的，本发明第二方面，提供一种基于修正无网格迦辽金的光子传输装置，包括：一混浊介质模块，为光子传输的载体；一光强探测模块输入端与混浊介质模块的输出端连接，用于测量混浊介质模块的混浊介质表面的光通量密度；一三维外形扫描模块输入端与混浊介质模块的输出端连接，用于重建混浊介质表面的外形；一计算机模块输入端与光强探测模块的输出端和三维外形扫描模块的输出端连接，将光强探测模块和三维外形扫描模块输出的表面光通量密度和三维外形数据输入计算机模块，用于运行混浊介质中光子传输方法，得到混浊介质表面光通量密度；一图形显示模块计算机模块输出端连接，用于显示混浊介质中光通量密度。

[0011] 本发明的有益效果是：利用修正的迦辽金无网格方法求解光子在非匀质复杂混浊介质中的传输，以移动最小二乘近似为基础，仅仅需要在混浊介质中布置一系统的场节点，而不需要任何场节点连接信息和单元信息，这不仅避免了复杂的网格剖分过程，同时也能更加方便准确的描述非匀质复杂混浊介质；此外，为了解决现有无网格方法不能直接准确施加第一类和第三类边界条件的问题，对移动最小二乘近似形函数进行了修正，使得形函数具有狄拉克函数性质，可以直接准确地施加任何边界条件，解决了一般无网格方法不能直接准确施加第一类和第三类边界条件的问题；在求解过程中，结合了光学特性参数的先验知识，有效地解决了荧光光子在非匀质复杂混沌介质中的传输问题。

附图说明

[0012] 图1为混浊介质光子传输装置方框图。

[0013] 图2为用于仿真实验的二维非匀质混浊介质仿体，其中包括图2a和图2b。

[0014] 图3为基于二维非匀质混浊介质仿体，本发明仿真结果与有限元方法计算结果的比较，其中包括图3a和图3b。

[0015] 图4为用于仿真实验的三维非匀质混浊介质仿体，其中包括图4a和图4b。

[0016] 图5为基于三维非匀质混浊介质仿体，本发明仿真结果与有限元方法计算结果的比较，其中包括图5a和图5b。

具体实施方式

[0017] 下面将结合附图对本发明加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

[0018] 本发明的具体实施例如下：

[0019] 1.在混浊介质Ω中布置一系列场节点，并确定混浊介质中每个组成部分的光学特性参数，用作为先验知识，以体现混浊介质中各个组成部分光学特性参数的个体差异，并解决复杂混浊介质的非匀质特性。同时，生成用于数值积分的背景网格，求解混浊介质中每个场节点上的移动最小二乘形函数。移动最小二乘形函数求解在混浊介质中布置一系列的场节点，用移动最小二乘近似求解得到混浊是介质中每个场节点上的移动最小二乘形函数。

[0020] 2.确定背景网格内的高斯积分点：针对高斯积分点，在每个高斯积分点支持域内搜索用于插值的节点，并计算插值节点处的移动最小二乘近似形函数。利用现有的完全变换法对计算得到的移动最小二乘近似形函数进行修正，使其满足狄拉克Delta函数性质。在背景网格内选取高斯积分点，对混浊介质边界和整个区域进行高斯数值积分。

[0021] 3.利用修正的移动最小二乘形函数对混浊介质中的光通量密度Φ(x)进行插值h得到混浊介质中的近似光通量密度Φ(x)，并利用迦辽金方法获得用于描述光子传输的扩散方程和第三类Robin边界条件的弱解形式；

[0022] 4.将 混 浊 介 质 中 的 近 似 光 通 量 密 度 Φh(x) 代 入 扩 散 方 程和第三类边界条件的弱解形式结合迦辽金方法和修正移动最
小二乘近似，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求解混浊介质中任意一点的光通量密度。

[0023] 其中，Ω和 分别是混浊介质区域及其边界；Φ(x)是在点x处的光通量密度；D(x)是与位置相关的光子扩散系数，其表示如下：

[0024]

[0025] 其中，μa(x)为吸收系数，μx(x)为散射系数，g是各向异性参数；S(x)是荧光光源密度；v(x)是混浊介质边界的外法线单位向量；A(x；n，n′)是一个表示混浊介质和周边介质折射系数不匹配的常数，它可以近似的表示为：

[0026]

[0027] 其中，n和n′分别表示为混浊介质和周边介质的折射系数，如若周边介质为空气，那么n′≈1.0，则R(x)可以近似表示为：

[0028] R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n。

[0029] 下面对上述的步骤2加以详细说明，具体形式如下面所述：

[0030] 2.1根据移动最小二乘近似，光通量密度Φ(x)的近似表达式如下：

[0031]

[0032] 其中，Ni(x)为移动最小二乘近似形函数，为光通量密度在第i节点上的参数。Ni(x)可以表示为：

[0033] Ni(x)＝pT(x)A-1(x)Bi(x)

[0034] 其中pT(x)为多项式基函数，Bi(x)为矩阵B(x)的第i列，矩阵A(x)和B(x)的表达式如下：

[0035]

[0036] 式中的w(x-xi)，i＝1，…，Nn为权函数。

[0037] 2.2利用下面修正的第l个场节点上的移动最小二乘形函数Ml(x)的公式对移动最小二乘近似形函数Nl(x)进行修正，使其满足狄拉克函数特性，并可以直接准确的施加边界条件：

[0038]

[0039] 其中，Ml(x)为修正的移动最小二乘近似形函数，Nn为混浊介质中所有场节点的个数，Nl(x)和Ni(x)分别为第l和i个场节点上的移动最小二乘形函数，Ml(x)为修正的第l个场节点上的移动最小二乘形函数，l和i表示场节点的序号。

[0040] 下面对上述的步骤3进行详细说明，具体形式如下面所述：

[0041] 3.1混浊介质对光子的作用主要表现为散射和吸收，辐射传输方程是能够精确描述光子在混浊介质中传输的数学模型，但是由于辐射传输方程是一个积分-微分方程，对该方程进行求解是非常耗时的。然而混浊介质对于光子是高散射低吸收的媒质，所以扩散方程能够对光子传输进行近似模拟。当光学成像实验工作于连续波(Continuous Wave)方式时，光源可以认为是稳定的，于是可以得到稳态的扩散方程：

[0042]

[0043] 当光学成像实验在一个全黑的环境下进行时，进一步考虑到混浊介质折射系数n和外部媒介折射系数n′的非匹配性，上述扩散方程的边界条件可以被表达为：

[0044]

[0045] 3.2基于扩散方程和第三类边界条件，利用迦辽金方法和高斯理论可以得到Φ(x)的弱解形式：

[0046]

[0047]1 1

[0048] 其中，Ψ(x)∈H(Ω)，H(Ω)为索伯列夫空间；D(x)是与位置相关的光子扩散系数，Φ(x)是在点x处的光通量密度，Ψ(x)是测试函数，μa(x)表示吸收系数，An(x)是一个表示混浊介质和周边介质折射系数不匹配的常数，S(x)表示混浊介质内的光源密度。

[0049] 3.3利用修正的移动最小二乘近似形函数，光通量密度Φ(x)的近似表达式如下：

[0050]

[0051] 其中，Mi(x)为第i节点上的修正移动最小二乘形函数；φi为第i个节点上的光通量密度；Nn为混浊介质内所有场节点的个数。h

[0052] 3.4把光通量密度Φ(x)的近似表达式Φ(x)代入到扩散方程和第三类边界条件的弱解形式中，可以得到如下的矩阵方程：

[0053] (K+C+F)Φ＝GΦ＝S

[0054] 其中，矩阵K、C、F和向量S的分量可以表示为：

[0055]

[0056] 3.5由于矩阵G是一个对称正定矩阵，所以可以利用下面的公式来求解混浊介质中场节点上的光流量密度：

[0057] Φ＝G-1S

[0058] 3.6利用混浊介质中场节点上的光通量密度和插值方法，可求得混浊介质中任意一点的光通量密度。

[0059] 运行结果

[0060] 为了验证本发明的方法，我们利用三维非匀质混浊介质仿体进行了实验：

[0061] 利用图1所示的混浊介质光子传输装置可以对三维非匀质混浊介质仿体中光子传输进行实验验证，采用一台具有2.8G赫兹中央处理器和1G字节内存的奔腾4计算机并用Fortran语言编制了光子传输算法程序，实现了本发明的混浊介质光子传输装置，所述的混浊介质光子传输装置包括：混浊介质模块11，为光子传输的载体，光强探测模块12、三维外形扫描模块13、计算机模块14和图形显示模块15，混浊介质模块11，为光子传输的载体；光强探测模块12输入端与混浊介质模块11的输出端连接，利用光强探测模块12测量混浊介质模块11的混浊介质表面的光通量密度；三维外形扫描模块13输入端与混浊介质模块11的输出端连接，利用三维外形扫描模块13对混浊介质模块11的混浊介质外形进行重建，计算机模块14输入端与光强探测模块12的输出端和三维外形扫描模块13的输出端连接，将光强探测模块12和三维外形扫描模块13输出的表面光通量密度和三维外形数据输入计算机模块14，运行二维到三维光强映射程序，得到三维非匀质混浊介质仿体表面光通量密度。同时，利用计算机模块14运行混浊介质中光子传输算法，也得到三维非匀质混浊介质仿体表面光通量密度。最后，图形显示模块计算机模块14输出端连接，利用图形显示模块15将三维非匀质混浊介质仿体表面光通量密度进行显示。

[0062] 第一个实验利用二维非匀质混浊介质仿体，请参见图2a表示的是带有两个光源3的非匀质仿体，它由区域1和区域2两部分组成，并在该非匀质仿体中布置了场节点，对基于修正无网格迦辽金的光子传输方法进行了验证。该非匀质混浊仿体是一个边长为30mm-1 -1的正方形，其中，区域1的光学特性参数为：μa＝0.035mm ，μs＝6.0mm ，g＝0.9，n＝-1 -1
1.37；区域2的光学特性参数为：μa＝0.35mm ，μs＝23.0mm ，g＝0.94，n＝1.37。根据本发明所提出的这种基于基于修正无网格迦辽金的光子传输方法，验证主要包含下列步骤：

[0063] 1.在区域2中嵌入两个正方形荧光光源3，两个光源的中心分别位于(13.5，13.5)2
和(18.5，18.5)，其边长均为1.0mm，密度均为1.0nano-Watts/mm。在混浊介质内布置
31×31均匀分布的场节点，如图2a所示。

[0064] 2.利用本发明所提出的这种基于修正无网格迦辽金的光子传输方法，计算二维非匀质混浊介质中的光通量密度，请参见图2b示出的基于二维非匀质混浊介质仿体，本发明求解的光通量密度分布结果。

[0065] 3.利用有限元算法计算二维非匀质混浊介质中的光通量密度，请参见图3a示的光通量密度分布计算结果。

[0066] 4.比较本发明所提出方法和有限元算法在二维非匀质混浊介质边界处的光通量密度，如图3b示出本发明与有限元方法的计算结果的比较，图中曲线用点“·”表示本发明方法，用方格加直线条 表示为有限元算法。

[0067] 第二个实验利用三维非匀质混浊介质仿体，如图4a表示的是带有两个光源3的非匀质仿体，它由区域1和区域2两部分组成，并在该非匀质仿体中布置了场节点，对基于修正无网格迦辽金的光子传输方法进行了验证。该非匀质混浊仿体是一个边长为20mm的正-1 -1方体。其中，区域1的光学特性参数为：μa＝0.035mm ，μs＝6.0mm ，g＝0.9，n＝1.37；
-1 -1
区域2的光学特性参数为：μa＝0.01mm ，μs＝4.0mm ，g＝0.9，n＝1.37。根据本发明所提出的这种基于基于修正无网格迦辽金的光子传输方法，验证主要包含下列步骤：

[0068] 1.在区域2中嵌入两个正方体荧光光源3，两个光源的中心分别位于(6.25，6.25，2
11.25)和(13.75，13.75，11.25)，其边长均为2.5mm，密度均为1.0nano-Watts/mm。在混浊介质内布置9×9×9均匀分布的场节点，如图4a所示。

[0069] 2.利用本发明所提出的这种基于修正无网格迦辽金的光子传输方法，计算三维非匀质混浊介质中的光通量密度，请参见图4b示出求解的光通量密度分布结果。

[0070] 3.利用有限元算法计算三维非匀质混浊介质中的光通量密度，请参见图5a示出的光通量密度分布计算结果。

[0071] 4.比较本发明所提出方法和有限元算法在三维非匀质混浊介质边界处的光通量密度，请参见图5b示出的计算结果比较，图中曲线用点“×”表示本发明方法，用方格加直线条 表示为有限元算法。

[0072] 前面已经具体描述了本发明的实施方案，应当理解，对于一个具有本技术领域的普通技能的人，在不背离本发明的范围的情况下，在上述的和在附加的权利要求中特别提出的本发明的范围内进行变化和调整能同样达到本发明的目的。