辐射能传播蒙特卡洛算法中的局部网格及光子加密方法

摘要

本发明提供一种辐射能传播蒙特卡洛算法中的局部网格及光子加密方法：首先，对传播介质内不同组成部分依据对应组成部分的光学特性及形状结构选用尺寸不同的网格；其次，根据最稀疏网格的网格尺寸选取所述传播介质内光传播计算所需的总光子数，总光子数的大小应满足按所述最稀疏网格的网格尺寸计算时网格内光子能量分布的统计学要求，当光子由网格尺寸较大的稀疏网格进入网格尺寸较小的细密网格时，光子进行分裂，本发明在现有“基于网格的蒙特卡洛方法”的基础上，依据传播介质的光学特性及形状结构，提出局部网格及光子加密方法，极大地提高了蒙特卡洛方法的计算效率。

说明书

技术领域

本发明属于热辐射计算领域，主要涉及针对不同传播介质或具 有多种组成的传播介质内辐射能传播蒙特卡洛计算方法的改进算法。

背景技术

辐射传播是热量传递的主要三种形势之一，在工业生产及生活领 域具有广泛的应用。例如利用复色光辐射能——太阳能进行发电，或 利用单色光辐射能——激光进行工业材料处理或临床激光热疗。由于 辐射能在介质内部(如太阳能集热器、工业材料或生物组织内)的传 播并不可见，因而数值计算成为了一种有效地“可视化”方法，以增强 对辐射能与传播介质间相互作用的理解并以此提高其工程实际应用 的效果。在众多的热辐射计算方法中蒙特卡洛方法是目前应用最广的 辐射能传播及分布求解方法之一。

蒙特卡洛方法是利用统计学模型求解物理过程的数学方法。由于 辐射能传播介质的几何形状往往具有高度复杂性，光子在介质内不同 组分界面(以下简称界面)上的折射及散射计算是蒙特卡洛方法实施 中的关键问题。Wang等的多层蒙特卡洛模拟方法(Multi-layered  Monte Carlo method)是目前应用较广的光传播蒙特卡洛方法(Lihong  W,Jacques SL,Liqiong Z.MCML-Monte Carlo modeling of light  transport in multi-layered tissues.Computer Methods and Programs in  Biomedicine,1995,47(2):131-146.)，然而该方法只能处理与坐标平面 垂直的界面。Pfefer等人开发的“基于网格的蒙特卡洛方法” (Voxel-based Monte Carlo method)是目前公认的可以较好处理这类 问题的蒙特卡洛方法(Pfefer TJ,Barton JK,Chan EK,et al.A  three-dimensional modular adaptable grid numerical model for light  propagation during laser irradiation of skin tissue.Ieee Journal of  Selected Topics in Quantum Electronics,1996,2(4):934-942.)。与传统 蒙特卡洛方法中网格单纯的能量统计功能不同，基于网格的蒙特卡洛 方法最大的特点就是光子的运行与网格设置紧密相连，连续的复杂曲 面被离散的微小的网格界面所替代，利用光子与网格界面的相互作用 替代与曲面间的相互作用。通过该方法极大地简化了复杂界面上折 射、反射的求解难度。

Binzoni等人研究结果表明，“基于网格的蒙特卡洛计算方法”中， 由于网格尺寸选取不当造成的计算误差可达20％以上(Binzoni T, Leung TS,Giust R,et al.Light transport in tissue by 3D Monte Carlo: Influence of boundary voxelization.Computer Methods and Programs in  Biomedicine,2008,89(1):14-23.)。另一方面，蒙特卡洛方法具有统计 学特性，其结果依赖于光子数的选择，光子数越多其计算结果越接近 真实结果，网格数越多所需要的光子数也越多。因而，蒙特卡洛方法 的计算效率与网格及光子数的选取有着密切的关系。但是，截至目前 为止，针对蒙特卡洛方法中网格尺寸及光子数选择判据的研究依然未 见报道。蒙特卡洛计算中网格尺寸及光子数的选择多根据经验而无科 学依据。

发明内容

本发明的目的在于提供一种辐射能传播蒙特卡洛算法中的局部 网格及光子加密方法，可以显著提高基于网络的蒙特卡洛算法计算辐 射能传播及分布时的效率。

为达到上述目的，本发明采用了以下技术方案：

首先，对辐射能传播介质(以下简称介质)不同组成部分依据该 介质对应组成部分的光学特性及形状结构选用尺寸不同的网格；其 次，根据最稀疏网格的网格尺寸选取所述介质内辐射能传播蒙特卡洛 计算所需的总光子数，总光子数的大小应满足按所述最稀疏网格的网 格尺寸计算时网格内光子能量分布的统计学要求，当光子由网格尺寸 较大的稀疏网格进入网格尺寸较小的细密网格时，光子进行分裂，分 裂后光子的数目应满足所述细密网格内的光子能量分布的统计学要 求。

对所述介质中光吸收较弱的组成部分采用网格尺寸较大的稀疏 网格，而对所述介质中光吸收较强的组成部分采用网格尺寸较小的细 密网格；对所述介质中形状结构较简单的组成部分采用网格尺寸较大 的稀疏网格，而对所述介质中形状结构较复杂的组成部分采用网格尺 寸较小的细密网格。

所述介质不同组成部分的网格尺寸的最大允许值依据如下关系 式进行选取：

$\left(\begin{array}{cc}\Delta w^{*}=0.2r_{\max }^{*}& r_{\max }^{*}\le 1\\ \Delta w^{*}=0.2& r_{\max }^{*}\ge 1\end{array}\right)---\left(1\right)$

其中，Δw^*是无量纲网格尺寸，Δw^*＝Δwμ_t，Δw为实际网格尺寸， μ_t为辐射能在介质中的衰减系数，是介质中不同组成部分边界上 无量纲曲率半径最大值，r_max是介质中不同组成部分边 界上曲率半径最大值。

为满足光子能量分布的统计学要求，单个网格内平均光子数应大 于等于5。

本发明的有益效果体现在：

本发明在现有基于网格的蒙特卡洛方法的基础上，依据介质的光 学特性及形状结构，提出局部网格及光子加密方法，极大地提高了蒙 特卡洛方法的计算效率。

附图说明

图1是含有单根圆柱辐射能吸收体的传播介质示意图；

图2是网格加密区示意图；

图3是非加密区网格划分示意图；

图4是加密区网格划分示意图；

图5是光子数加密示意图；

图6是网格局部加密对蒙特卡洛计算时间(Computation time) 的影响；

图7是光子增值倍数对蒙特卡洛计算时间的影响；

图8是无光子分裂(a)及光子分裂10倍(b)及100倍(c)后，加密区 光密度分布图。

具体实施方式

下面结合附图和实施例对本发明做详细描述。

利用蒙特卡洛方法计算介质辐射能传播及吸收时网格及光子加 密的基本思路是：依据介质各组成部分(简称成分)的形状结构及光 学参数，针对各个组成部分划分不同尺寸的网格。光学吸收性强，且 形状复杂的组成部分要选取较小的网格尺寸。随后，依据计算区域内 最大网格选取蒙特卡洛计算所需的光子数，当光子运行过程中进入较 细密的网格区域时，光子分裂为若干个相同权值的子光子，以满足细 密网格内蒙特卡洛方法对于光子数的需求。

本发明首先从吸收能量的准确统计及形状不规则界面的准确刻 画这两方面要求出发，推导出网格尺寸选取的判据：网格大小的选取 与介质的光学性质(μ_t)及介质内界面的曲率半径(r)有关，光在介质 内衰减的越快或介质内界面曲率越大，则所需网格尺寸越小，反之则 越大。无量纲网格尺寸可依据如下关系式进行选取：

$\left(\begin{array}{cc}\Delta w^{*}=0.2r_{\max }^{*}& r_{\max }^{*}\le 1\\ \Delta w^{*}=0.2& r_{\max }^{*}\ge 1\end{array}\right)---\left(1\right)$

其中，Δw^*是无量纲网格尺寸，Δw^*＝Δwμ_t，Δw为实际网格尺寸， r_max是介质中不同组成部分边界上曲率半径最大值，即 介质不同组成部分边界上曲率最大值，是介质不同组成部分边界 上无量纲曲率半径最大值，μ_t为辐射能在所述介质内的衰减系数。

另外，利用单相高吸收介质内能量密度分布一维分析解与不同光 子数蒙特卡洛方法计算结果进行对比，引出光子数的选择判据：计算 中光子数的选择与网格尺寸及数量关系紧密，网格越小，网格数量越 大，则计算所需光子数越多。单个网格内平均光子数应大于等于5， 即光子使用总数应大于等于5倍网格数。

随后，本发明根据所得判据发展出“网格及光子加密技术”，可显 著提高辐射能传输蒙特卡洛算法的计算效率。通过以上判据可以看 出：介质各成分间网格尺寸选取值并不相同，辐射能吸收较强及形状 较复杂的介质成分所需的网格尺寸小于辐射能吸收较弱、结构简单的 介质成分；并且光子数的选取与网格大小相关，为了满足统计学特性， 网格越小所需的光子数就越多。在复杂介质内辐射能传播及分布蒙特 卡洛计算中，对不同介质成分采用均匀网格时，网格尺寸选取必须以 强吸收、形状复杂的成分为依据，但这对于弱吸收、形状简单的成分 而言显然是一种计算资源的浪费。因而，本发明在介质内建立非均匀 网格，依据介质的辐射能吸收特性及形状对网格进行局部加密，可以 提高“基于网格的蒙特卡洛方法”的计算效率。另一方面，依据判据使 用非均匀网格时，若光子数的选取以小网格为依据，则对于大网格而 言也是一种计算资源的浪费。因而在使用非均匀网格模拟时，应对网 格局部加密区域辅以光子分裂技术，以进一步提高“基于网格的蒙特 卡洛方法”的计算效率。

以下，本发明以典型单色光辐射能应用——血管性皮肤病激光热 疗中辐射能传输及分布模拟为例，说明网格及光子局部加密技术实施 的基本思路及其对计算效率的影响。针对图1所示的单根血管皮肤模 型，选择对治疗激光吸收性较强的血管区域为网格加密区。在网格划 分中，介质各成分间网格尺寸的选择应遵循网格选取准则。含有单根 血管皮肤组织的非均匀网格划分的具体过程如下：

首先，确定加密区的范围。加密区应该包含整个血管，可以略大 于血管面积，如图2所示。

随后，对吸收特性较弱的表皮及真皮区域采用稀疏网格进行网格 划分。网格划分结果如图3所示。

最后，对包含血管的区域采用细密网格进行网格局部加密。整体 网格划分结果二维图如图4所示。本发明所采用的算例中，血管轴平 行于y轴，且与xz平面垂直，因而在本发明的网格划分时对其他xz 平面采用相同操作即可。

光子局部加密：计算中光子数的选择与网格尺寸及数量关系紧 密，网格越小，所需的光子数越多，单个网格内平均光子数应大于等 于5。光子使用总数应大于等于5倍网格数。在局部网格及光子加密 技术中，非加密区光子数的选择，以及加密区光子数的加密倍数亦应 遵循各选取准则。

首先，根据网格尺寸较大区域(如图4表皮及真皮区域)计算所 需的总光子数，其大小为以当地网格尺寸计算所得的整个计算区域网 格个数的五倍。

另一方面，为了满足细密网格区域(图4中血管附近区域)每个 网格内的统计值具有统计学意义，光子进入该区域时应进行光子分 裂，光子的分裂倍数应使分裂后光子数满足细密网格内网格数需求。

光子加密的具体实施方法是：首先，规定一个光子加密区域，该 区域可以大于也可以等于网格加密区，如图5所示。假设，模拟中第 k个光子碰到该区域边界，则该光子分裂成n个，同时为了保证能量 守恒，每个新光子的能量为原来的1/n。新光子一个接一个进入光子 加密区进行传播，直到n个光子走完为止再重新生成第k+1个光子。 由皮肤表面进入组织。

本发明首先采用无网格加密及光子加密的原始程序对图1所示 组织模型的光传播进行计算，网格尺寸为了满足血液内的吸收统计特 性取为10μm，光子数为10⁷个。在HP xw6600Workstation上所需的 计算时间为17.5个小时。以此结果作为标准结果，来说明网格及光 子加密对于计算效率的提高及计算精度的影响。

采用非均匀网格技术，不同真皮网格尺寸下蒙特卡洛计算所需的 时间如图6所示。对直径为120μm的血管而言，衰减系数为191cm^-1， 对于圆柱形血管来说，r_max＝60μm，根据公式(1)计算可得网格尺寸Δw＝0.2/μ_t＝10.47μm。同 时对于表皮而言，衰减系数为79.38cm^-1，其与真皮组织的边界为平 面，因而r_max趋近于无穷大，根据公式(1)，Δw＝0.2/μ_t＝25.19μm。 对于真皮而言，与表皮类似，Δw＝0.2/μ_t＝0.2/13.0＝143.9μm，即真 皮网格尺寸小于等于143.9μm即可。因而，计算中选取表皮层内网格 尺寸为20μm，加密区域网格尺寸为10μm，计算采用的光子数为10⁷个。从图中可以看出，当采用非均匀网格时计算时间明显短于均匀网 格。并且，随着真皮内网格尺寸的增加，程序运行的时间不断缩短。 当真皮中网格选取为80μm时，程序运行时间仅需要5.5个小时。说 明了网格加密技术对计算效率的提高程度。

将光子加密技术与网格加密技术相结合，光子加密倍数n其对计 算时间的影响如图7所示。从图中可以看出，仅仅采用网格加密可将 计算时间由17.5小时缩短至5.5个小时。若采用加密网格的同时，在 加密区内将光子增值10倍，可将计算时间缩短为3.5个小时。然而 将n提高到100时，计算时间为2.3个小时，比原始计算时间17.5小 时明显降低。

当计算中选取表皮层内网格尺寸为20μm，真皮区域网格尺寸为 80μm，加密区域网格尺寸为10μm时，从图8中可以看出，当光子加 密数n分别为10倍(图8b)及100倍(图8c)时，加密区光密度分 布与图8a无光子加密时基本相同。因而可以说明光子加密可以保证 蒙特卡洛方法的计算精度。

本发明的特点如下：针对形状简单、光吸收特性较弱的介质区域 采用较粗的计算网格，而在形状复杂、光吸收特性较强的组成部分实 施网格局部加密。蒙特卡洛计算时光子数仅需满足粗网格计算需求即 可，仅对局部网格加密区域计算时实施光子分裂，使局部加密后的光 子数满足加密网格区域的计算需求。

总之，针对复杂结构介质内辐射能传播蒙特卡洛计算时，依据本 发明提出的蒙特卡洛方法网格及光子数选取判据，对介质不同组成部 分采取局部网格及光子加密，满足计算精度的同时降低计算量，提高 了计算效率。