专利名称：一种混浊介质中光子传输方法及装置的制作方法
技术领域：
本发明属于分子影像领域，涉及混浊介质中光子传输方法，尤其是涉及一种基于
修正无网格迦辽金的光子传输方法。
背景技术：
光学分子影像是近年来新兴的一种分子影像模态，光学分子影像主要包括自发荧光和激发荧光两种方式。自发荧光是利用荧光素酶标记混浊介质中的目标，在ATP以及氧气存在的情况下，若荧光素酶遇到底物荧光素，荧光素酶将会催化荧光素的氧化反应并产生光子，产生光的强度与标记目标的数量成正比。激发荧光技术则采用荧光报告基团标记目标，荧光报告基团吸收外界激发光源射入混浊介质中的光子，产生能量跃迁至受激虚态，从受激虚态向下跃迁至基态过程中，产生荧光光子。 由于混浊介质具有强散射特性，光子在其内部的传输不再沿直线传播，而是经过大量的散射过程，从而远远偏离了原来的运行轨道。因此，光子在混浊介质中的传输是一个非常复杂的过程。在光学分子影像领域，有限元算法是一种经典的解决光子在复杂混浊介质中传播的数值计算方法。然而，在使用有限元算法之前，必须对混浊介质进行网格剖分，并对网格数据进行预处理。虽然目前已研发了多种网格剖分软件，但是对具有复杂内部结构的非规则混浊介质来说，进行网格剖分还是比较困难和耗时的。另外，由于混浊介质的边界并非是一个固定的边界，而是一个随时间移动的边界，因此在光子传输问题求解过程中，也要随时间对混浊介质进行有限元网格剖分。此外，在有限元自适应技术实现过程中，在两次迭代之间需要对混浊介质进行网格重构、调整网格节点间连接信息，因此实现起来也是比较复杂耗时的。
发明内容
为了解决光子传输有限元分析存在的网格剖分困难，本发明的目的是基于修正
无网格迦辽金法求解光子在非匀质复杂混浊介质中的传输，此外，为了解决现有无网格方
法不能直接准确施加第三类边界条件的问题，本发明对移动最小二乘近似形函数进行了修
正，使得形函数具有狄拉克函数性质，可以直接准确地施加任何边界条件。 为了实现所述的目的，本发明第一方面，提供一种基于修正无网格迦辽金的光子
传输方法的技术方案包括如下步骤 步骤1 :在混浊介质中布置多个离散的场节点和用于数值积分的背景网格，用于确定混浊介质中每个组成部分的光学特性参数后，来求解混浊介质中每个场节点上的移动最小二乘形函数； 步骤2 :利用完全变换法修正移动最小二乘形函数，使其满足狄拉克函数特性具有狄拉克(Delta)函数特性，用于直接准确的施加任何边界条件； 步骤3 :利用修正的移动最小二乘形函数对混浊介质的光通量密度①(X)进行插值得到混浊介质的近似光通量密度Oh(X)，并利用迦辽金方法获得用于描述光子传输的扩
4散方程和第三类(Robin)边界条件的弱解形式； 步骤4 :将混浊介质的近似光通量密度Oh(X)代入扩散方程和第三类边界条件的弱解形式，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。 为了实现所述的目的，本发明第二方面，提供一种基于修正无网格迦辽金的光子传输装置，包括一混浊介质模块，为光子传输的载体；一光强探测模块输入端与混浊介质模块的输出端连接，用于测量混浊介质模块的混浊介质表面的光通量密度；一三维外形扫描模块输入端与混浊介质模块的输出端连接，用于重建混浊介质表面的外形；一计算机模块输入端与光强探测模块的输出端和三维外形扫描模块的输出端连接，将光强探测模块和三维外形扫描模块输出的表面光通量密度和三维外形数据输入计算机模块，用于运行混浊介质中光子传输方法，得到混浊介质表面光通量密度；一图形显示模块计算机模块输出端连接，用于显示混浊介质中光通量密度。 本发明的有益效果是利用修正的迦辽金无网格方法求解光子在非匀质复杂混浊介质中的传输，以移动最小二乘近似为基础，仅仅需要在混浊介质中布置一系统的场节点，而不需要任何场节点连接信息和单元信息，这不仅避免了复杂的网格剖分过程，同时也能更加方便准确的描述非匀质复杂混浊介质；此外，为了解决现有无网格方法不能直接准确施加第一类和第三类边界条件的问题，对移动最小二乘近似形函数进行了修正，使得形函数具有狄拉克函数性质，可以直接准确地施加任何边界条件，解决了一般无网格方法不能直接准确施加第一类和第三类边界条件的问题；在求解过程中，结合了光学特性参数的先验知识，有效地解决了荧光光子在非匀质复杂混沌介质中的传输问题。
图1为混浊介质光子传输装置方框图。 图2为用于仿真实验的二维非匀质混浊介质仿体，其中包括图2a和图2b。
图3为基于二维非匀质混浊介质仿体，本发明仿真结果与有限元方法计算结果的比较，其中包括图3a和图3b。 图4为用于仿真实验的三维非匀质混浊介质仿体，其中包括图4a和图4b。
图5为基于三维非匀质混浊介质仿体，本发明仿真结果与有限元方法计算结果的比较，其中包括图5a和图5b。
具体实施例方式
下面将结合附图对本发明加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。
本发明的具体实施例如下 1.在混浊介质Q中布置一系列场节点，并确定混浊介质中每个组成部分的光学特性参数，用作为先验知识，以体现混浊介质中各个组成部分光学特性参数的个体差异，并解决复杂混浊介质的非匀质特性。同时，生成用于数值积分的背景网格，求解混浊介质中每个场节点上的移动最小二乘形函数。移动最小二乘形函数求解在混浊介质中布置一系列的场节点，用移动最小二乘近似求解得到混浊是介质中每个场节点上的移动最小二乘形函
5
2.确定背景网格内的高斯积分点针对高斯积分点，在每个高斯积分点支持域内
搜索用于插值的节点，并计算插值节点处的移动最小二乘近似形函数。利用现有的完全变
换法对计算得到的移动最小二乘近似形函数进行修正，使其满足狄拉克Delta函数性质。
在背景网格内选取高斯积分点，对混浊介质边界和整个区域进行高斯数值积分。 3.利用修正的移动最小二乘形函数对混浊介质中的光通量密度①(x)进行插值
得到混浊介质中的近似光通量密度①h(x)，并利用迦辽金方法获得用于描述光子传输的扩
散方程和第三类Robin边界条件的弱解形式； 4.将混浊介质中的近似光通量密度Oh(X)代入扩散方程
(—▽ po)v①(匀]+ = S(X)， x e Q )和第三类边界条件的弱解形式
(d>0) + 2^0;",W)DO)(vO:). V(DO)) = 0， x e )，结合迦辽金方法和修正移动最小二乘近似，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求解混浊介质中任意一点的光通量密度。 其中，Q和3 Q分别是混浊介质区域及其边界；①(X)是在点X处的光通量密度；D(X)是与位置相关的光子扩散系数，其表示如下= :~~、 , 、n 其中，i!a(x)为吸收系数，i!s(X)为散射系数，g是各向异性参数；S(X)是荧光光源密度；v(x)是混浊介质边界的外法线单位向量；A(x ;n， r/)是一个表示混浊介质和周边气關 ,.^閱，
1 + 7 00
丌质折射系数不匹配的常数，它可以近似的表示为
那么r/
l一雖)
其中，n和r/分别表示为混浊介质和周边介质的折射系数，如若周边介质为空气，- l.O，则R(x)可以近似表示为:R(x) " -1. 4399n—2+0. 7099n—、0. 6681+0. 0636n。下面对上述的步骤2加以详细说明，具体形式如下面所述2. 1根据移动最小二乘近似，光通量密度①(x)的近似表达式如下
i.=i 其中，Ni(x)为移动最小二乘近似形函数，《为光通量密度在第i节点上的参数。
K(X)可以表示为: & (x) = pT (x) A—1 (x) Bi (x) 其中pT(x)为多项式基函数，Bi(x)为矩阵B(x)的第i列，矩阵A(x)和B(x)的表达式如下<formula>formula see original document page 6</formula><formula>formula see original document page 7</formula>Nn为权函数。 2.2利用下面修正的第l个场节点上的移动最小二乘形函数Mjx)的公式对移动最小二乘近似形函数& (x)进行修正，使其满足狄拉克函数特性，并可以直接准确的施加边界条件
w M, O) = J] M (JC)乂 O,. )_1
,■=1 其中，Mjx)为修正的移动最小二乘近似形函数，Nn为混浊介质中所有场节点的个数，Njx)和Ni(x)分别为第l和i个场节点上的移动最小二乘形函数，Mjx)为修正的第l个场节点上的移动最小二乘形函数，1和i表示场节点的序号。
下面对上述的步骤3进行详细说明，具体形式如下面所述 3. 1混浊介质对光子的作用主要表现为散射和吸收，辐射传输方程是能够精确描述光子在混浊介质中传输的数学模型，但是由于辐射传输方程是一个积分-微分方程，对该方程进行求解是非常耗时的。然而混浊介质对于光子是高散射低吸收的媒质，所以扩散方程能够对光子传输进行近似模拟。当光学成像实验工作于连续波(Conti皿ous Wave)方式时，光源可以认为是稳定的，于是可以得到稳态的扩散方程
-▽ . ["(X)VO)(X)] + //。 O)O(义)=6"(X) (x e Q ) 当光学成像实验在一个全黑的环境下进行时，进一步考虑到混浊介质折射系数n和外部媒介折射系数r/的非匹配性，上述扩散方程的边界条件可以被表达为(D(x) + 2^0; w， )D(x)(v(jc). V①(x)) = 0 O e 3.2基于扩散方程和第三类边界条件，利用迦辽金方法和高斯理论可以得到①(X)的弱解形式
+ 其中，W(x) G tf(Q)，tf(Q)为索伯列夫空间；D(x)是与位置相关的光子扩散系数，O (x)是在点x处的光通量密度，W (x)是测试函数，y a(x)表示吸收系数，An(x)是一个表示混浊介质和周边介质折射系数不匹配的常数，S(x)表示混浊介质内的光源密度。
3. 3利用修正的移动最小二乘近似形函数，光通量密度①(x)的近似表达式如下 0^(X)-ZM,.0)《 其中，Mi(x)为第i节点上的修正移动最小二乘形函数；^为第i个节点上的光通量密度；Nn为混浊介质内所有场节点的个数。 3. 4把光通量密度①(x)的近似表达式Oh(X)代入到扩散方程和第三类边界条件的弱解形式中，可以得到如下的矩阵方程
(K+C+F)①=GO = S 其中，矩阵K、C、F和向量S的分量可以表示为
7
;=!^o)(vm,o)).(vm)cq
& =工qM,(x)M)(jc)/(2」0;"，"'))^q
3.5由于矩阵G是一个对称正定矩阵，所以可以利用下面的公式来求解混浊介质
中场节点上的光流量密度
①=6—5 3. 6利用混浊介质中场节点上的光通量密度和插值方法，可求得混浊介质中任意
一点的光通量密度。 运行结果 为了验证本发明的方法，我们利用三维非匀质混浊介质仿体进行了实验
利用图1所示的混浊介质光子传输装置可以对三维非匀质混浊介质仿体中光子传输进行实验验证，采用一台具有2. 8G赫兹中央处理器和1G字节内存的奔腾4计算机并用Fortran语言编制了光子传输算法程序，实现了本发明的混浊介质光子传输装置，所述的混浊介质光子传输装置包括混浊介质模块ll，为光子传输的载体，光强探测模块12、三维外形扫描模块13、计算机模块14和图形显示模块15，混浊介质模块ll，为光子传输的载体；光强探测模块12输入端与混浊介质模块11的输出端连接，利用光强探测模块12测量混浊介质模块11的混浊介质表面的光通量密度；三维外形扫描模块13输入端与混浊介质模块11的输出端连接，利用三维外形扫描模块13对混浊介质模块11的混浊介质外形进行重建，计算机模块14输入端与光强探测模块12的输出端和三维外形扫描模块13的输出端连接，将光强探测模块12和三维外形扫描模块13输出的表面光通量密度和三维外形数据输入计算机模块14，运行二维到三维光强映射程序，得到三维非匀质混浊介质仿体表面光通量密度。同时，利用计算机模块14运行混浊介质中光子传输算法，也得到三维非匀质混浊介质仿体表面光通量密度。最后，图形显示模块计算机模块14输出端连接，利用图形显示模块15将三维非匀质混浊介质仿体表面光通量密度进行显示。 第一个实验利用二维非匀质混浊介质仿体，请参见图2a表示的是带有两个光源3的非匀质仿体，它由区域1和区域2两部分组成，并在该非匀质仿体中布置了场节点，对基于修正无网格迦辽金的光子传输方法进行了验证。该非匀质混浊仿体是一个边长为30mm的正方形，其中，区域1的光学特性参数为ii a = 0. 035mm 、 P s = 6. 0mm 、 g = 0. 9， n =1. 37 ;区域2的光学特性参数为y a = 0. 35mm—、 P s = 23. 0mm—、 g = 0. 94， n = 1. 37。根据本发明所提出的这种基于基于修正无网格迦辽金的光子传输方法，验证主要包含下列步
1.在区域2中嵌入两个正方形荧光光源3，两个光源的中心分别位于(13.5，13. 5)
骤
和(18. 5，18. 5)，其边长均为1.0mm，密度均为1. 0nano-Watts/mm2。在混浊介质内布置31 X 31均匀分布的场节点，如图2a所示。 2.利用本发明所提出的这种基于修正无网格迦辽金的光子传输方法，计算二维非匀质混浊介质中的光通量密度，请参见图2b示出的基于二维非匀质混浊介质仿体，本发明
8求解的光通量密度分布结果。 3.利用有限元算法计算二维非匀质混浊介质中的光通量密度，请参见图3a示的光通量密度分布计算结果。 4.比较本发明所提出方法和有限元算法在二维非匀质混浊介质边界处的光通量密度，如图3b示出本发明与有限元方法的计算结果的比较，图中曲线用点" "表示本发明方法，用方格加直线条"~B~"表示为有限元算法。 第二个实验利用三维非匀质混浊介质仿体，如图4a表示的是带有两个光源3的非匀质仿体，它由区域1和区域2两部分组成，并在该非匀质仿体中布置了场节点，对基于修正无网格迦辽金的光子传输方法进行了验证。该非匀质混浊仿体是一个边长为20mm的正方体。其中，区域l的光学特性参数为iia = 0. 035mm 、 P s = 6. Omm 、 g = 0. 9， n = 1.37 ;区域2的光学特性参数为:y a = 0. Olmm 、 P s = 4. Omm 、g = 0. 9，n = 1. 37。根据本发明所提出的这种基于基于修正无网格迦辽金的光子传输方法，验证主要包含下列步骤
1.在区域2中嵌入两个正方体荧光光源3，两个光源的中心分别位于(6. 25，6. 25，11. 25)和(13. 75， 13. 75， 11. 25)，其边长均为2. 5mm，密度均为1. 0nano-Watts/mm2。在混浊介质内布置9X9X9均匀分布的场节点，如图4a所示。 2.利用本发明所提出的这种基于修正无网格迦辽金的光子传输方法，计算三维非
匀质混浊介质中的光通量密度，请参见图4b示出求解的光通量密度分布结果。 3.利用有限元算法计算三维非匀质混浊介质中的光通量密度，请参见图5a示出
的光通量密度分布计算结果。 4.比较本发明所提出方法和有限元算法在三维非匀质混浊介质边界处的光通量密度，请参见图5b示出的计算结果比较，图中曲线用点"X "表示本发明方法，用方格加直线条""B""表示为有限元算法。 前面已经具体描述了本发明的实施方案，应当理解，对于一个具有本技术领域的普通技能的人，在不背离本发明的范围的情况下，在上述的和在附加的权利要求中特别提出的本发明的范围内进行变化和调整能同样达到本发明的目的。
9
权利要求
一种混浊介质中光子传输的方法，其特征在于步骤1在混浊介质中布置多个离散的场节点和用于数值积分的背景网格，用于确定混浊介质中每个组成部分的光学特性参数后，来求解混浊介质中每个场节点上的移动最小二乘形函数；步骤2利用完全变换法修正移动最小二乘形函数，使其满足狄拉克函数特性具有狄拉克(Delta)函数特性，用于直接准确的施加任何边界条件；步骤3利用修正的移动最小二乘形函数对混浊介质的光通量密度Φ(x)进行插值得到混浊介质的近似光通量密度Φh(x)，并利用迦辽金方法获得用于描述光子传输的扩散方程和第三类(Robin)边界条件的弱解形式；步骤4将混浊介质的近似光通量密度Φh(x)代入扩散方程和第三类边界条件的弱解形式，得到矩阵方程，并求解矩阵方程，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。
2. 根据权利要求1所述的混浊介质中光子传输方法，其特征在于所述移动最小二乘 形函数求解是在混浊介质中布置一系列的场节点，用移动最小二乘近似求解得到混浊介质 中每个场节点上的移动最小二乘形函数。
3. 根据权利要求1所述的混浊介质中光子传输方法，其特征在于所述对移动最小二乘形函数Njx)进行修正是利用修正第l个场节点上的移动最小二乘形函数Mjx):其中，Nn为混浊介质中所有场节点的个数，Njx)为第1个场节点上的移动最小二乘形 函数，Mjx)为修正的第i个场节点上的移动最小二乘形函数，1和i表示场节点的序号。
4. 根据权利要求1所述的混浊介质中光子传输方法，其特征在于针对扩散方程和第三类边界条件，使用迦辽金方法和高斯理论，获得弱解形式为<formula>formula see original document page 2</formula>其中，W(X) G tf(Q)， tf(Q)为索伯列夫空间；D(X)是与位置相关的光子扩散系数， ①(x)是在点x处的光通量密度，W (x)是测试函数，y a(x)表示吸收系数，An(x)是一个表 示混浊介质和周边介质折射系数不匹配的常数，S(x)表示混浊介质内的光源密度。
5. 根据权利要求1所述的混浊介质中光子传输方法，其特征在于在背景网格内选取高斯积分点，对混浊介质边界和整个区域进行高斯数值积分。
6. 根据权利要求1所述的混浊介质中光子传输方法，其特征在于将所述混浊介质中每个组成部分的光学特性参数作为先验知识，表现为混浊介质中各个组成部分光学特性参 数的个体差异，用于解决复杂混浊介质的非匀质特性。
7. —种混浊介质中光子传输装置，其特征在于，包括一混浊介质模块，为光子传输的载体；一光强探测模块输入端与混浊介质模块的输出端连接，用于测量混浊介质模块的混浊介质表面的光通量密度；一三维外形扫描模块输入端与混浊介质模块的输出端连接，用于重建混浊介质表面的 外形；一计算机模块输入端与光强探测模块的输出端和三维外形扫描模块的输出端连接，将 光强探测模块和三维外形扫描模块输出的表面光通量密度和三维外形数据输入计算机模 块，用于运行混浊介质中光子传输方法，得到混浊介质表面光通量密度；一图形显示模块计算机模块输出端连接，用于显示混浊介质中光通量密度。
8.根据权利要求7所述混浊介质中光子传输装置，其特征在于，在所述混浊介质模块 的混浊介质中布置多个离散的场节点和背景网格，用于确定混浊介质中每个组成部分的光 学特性参数，求解混浊介质中每个场节点上的移动最小二乘形函数。
全文摘要
本发明涉及一种混浊介质中光子传输方法及装置，在混浊介质中布置多个离散的场节点，求解混浊介质中每个场节点上的移动最小二乘形函数；利用完全变换法修正移动最小二乘形函数，使其具有狄拉克函数特性；利用修正的移动最小二乘形函数进行插值得到混浊介质的近似光通量密度，并利用迦辽金方法获得扩散方程和第三类边界条件的弱解形式，得到矩阵阵方程并进行求解，得到混浊介质中场节点上的光通量密度，进而利用插值方法，求得混浊介质中任意一点的光通量密度。本发明仅布置场节点，而不需要任何场节点连接信息和单元信息，避免复杂的网格剖分过程，对移动最小二乘近似形函数进行修正，使得形函数具有狄拉克函数性质，可直接添加任何边界条件。
文档编号G01J1/00GK101738253SQ200810225780
公开日2010年6月16日 申请日期2008年11月12日 优先权日2008年11月12日
发明者刘凯, 杨鑫, 田捷, 秦承虎 申请人:中国科学院自动化研究所