蒙特卡罗,光在烟雾中传输.pdf

蒙特卡罗方法模拟光在烟雾中的传输光在烟雾中的传输是一个复杂的过程烟雾粒子对光的吸收作用使光的能量发生衰减， 对光的散射作用使光的传输方向发生变化并且由于烟雾粒子密度和尺度的分布不均匀，光的能量在烟雾中不同区域的衰减是不同的， 散射方向也有很大的差别 通过求解传输方程来解释传输过程比较复杂 而蒙特卡罗方法作为一种随机统计方法，比较灵活，通用性强，可以很好的模拟烟雾对光传输的影响一蒙特卡罗方法蒙特卡罗方法作为一种统计模拟的求解方法，具有清晰的物理意义，适用于模拟辐射传输过程蒙特卡罗方法按照其计算过程，可以分为正向方法和反向方法 正向蒙特卡罗方法是直接模拟光子由源到探测器的物理过程， 这种方法产生的光子位置比较分散；反向蒙特卡罗方法是将光子由接收器发出，反向追踪光子， 这种方法产生的光子位置集中，使得对接收器接收的散射光贡献最大的区域最容易被抽取并可以计算一个特定的方向上或立体角内的光强本文用来模拟光在烟雾中的传输使用的是正向蒙特卡罗方法蒙特卡罗方法求解辐射传输问题是把光在介质中的传输看作是相继的光子与介质中其他分子及气溶胶粒子随机碰撞事件的马尔科夫链求解过程可以简单描述为：首先从光源释放样本，在相应的模拟介质（由适当数量单元组成的模拟区域，每个单元赋予一定的粒子单次散射光学性质）中追踪他们的行走轨迹，在每一个碰撞位置，根据吸收概率确定样本的生存状态， 根据散射相函数确定散射后新的传输方向。

这种方法可以较好地处理包括多次散射、不规则散射相函数、非均匀介质及复杂边界的介质中的辐射传输问题，这是其它方法难以解决的问题蒙特卡罗方法的计算效率比较高，且不受计算机内存限制 但蒙特卡罗方法依然面临着如何在保证计算精度的同时提高数值模拟效率的问题所以在使用蒙特卡罗方法来模拟辐射传输过程时，除了需要获取准确的模拟结果，还需要了考虑蒙特卡罗方法的收敛速度蒙特卡罗方法的收敛速度主要取决于两个方面：模拟时间和计算结果精度其中精度与模拟样本数目的均方根成正比所以提高蒙特卡罗模拟效率的两种途径为： 减少固定样本数目的模拟时间或在相同时间内降低模拟结果的方差 目前使用的比较多的途径是缩减模拟结果的方差本次模拟过程中， 为了解决贡献过小的样本由于频繁抽样而消耗大量计算时间的问题， 采用的俄罗斯轮盘赌这种方差缩减技术二光在烟雾中传输的模拟模拟过程中， 由光源位置和接收器位置确定烟雾的大小，即认为在光源和接收器之间都充满了烟雾整个区域内的光学性质是不变的，包含了多种不同半径的烟雾粒子，给出每种粒子的衰减系数、 吸收系数、散射系数、相函数图1 给出了蒙特卡罗方法模拟光在烟雾中传输的流程图，主要包括四个部分：首先初始化烟雾的光学性质；第二步初始化光源（光子的状态） ；第三步模拟光子在散射介质中的随机行走过程，即自由程抽样过程；第四步模拟光子的碰撞事件，具体包括光子权重更新模拟、碰撞类型抽样、碰撞方向抽样。

重复模拟最后两个步骤，直至光子被散射介质吸收、 被接收器接收或逃逸出散射介质确定下一个碰撞点确定散射方向光子的新坐标和传输方向吸收光子权重 >门限值？碰撞碰撞 ?结束是否开始否是初始化烟雾的光学性质及 光子的能量、位置俄罗斯轮盘赌选择光子权重为0？是结束否图 1 蒙特卡罗方法模拟光在烟雾中传输的流程图（1）初始化过程1）初始化烟雾的光学性质确定烟雾的大小，输入不同半径烟雾粒子的衰减系数、吸收系数、散射系数以及相函数，并进行编号2）初始化光子的散射性质初始化光源的位置、 能量以及光源的传输方向， 也就是光子的发射方向，该方向用天顶角和方位角来表示设置模拟的光子数目， 模拟的光子数越多， 模拟结果越接近真实的传输过程， 但考虑计算机的内存和计算效率，应该设置恰当的光子数2）光子自由程模拟过程1）确定碰撞位置首先要确定相邻两个碰撞点之间的自由路径L在相邻两次碰撞之间，光子的输运问题服从指数分布，所以L为：rrl1ln(1) sar(2) 其中：1r是[0,1]之间的随机数；r是衰减系数；a是吸收系数；s是散射系数2）确定碰撞类型目前程序中采用的是抽样方法是：将不同尺度的烟雾粒子编号， 给出每种烟雾粒子的衰减系数、吸收系数、散射系数；抽取相应的编号，就确定了碰撞粒子的光学性质。

3）确定散射方向光子碰撞后的方向由散射相函数来决定假设烟雾粒子是球型粒子，那么光子与烟雾粒子的散射具有关于Z轴对称的特点， 也就是光子在垂直传输方向的平面内具有各向同性，方位角和散射角分别由以下公式确定：22r(3)222 2 311cos12(12)g arcggggr (4) 公式(4)是由 Henyey-Greenstein相函数得到的，2r和3r是[0,1]之间均匀分布的随机数其中不对称因子 g 的计算方法：不对称因子 g 是与波长和粒子半径有关的物理量，确定光源的波长后， 建立不对称因子查找表 （该表是按烟雾粒子半径的分布建立的） ，抽样得到烟雾粒子的半径后， 选取与之对应的不对称因子， 确定散射角度3）光子的新坐标与传输方向光子的新坐标：111111nnnnnnnnnlZZlYYlXX(5) 上式中：111(,,)nnnXYZ是光子散射前的坐标， 当光子发生碰撞时，u、v、由下式给出：coscossinsinsin(6) （4）光子权重更新定义散射介质的单次散射反照率0W：0/stW(7)按照统计加权的方法，在模拟过程中赋予每个光子一个初始权重01W， 权重与碰撞次数 m 的关系是（单次散射反照率相同的情况下） ：2 1100/()mmstmWWWWW(8)如果光子被吸收，则权重变为0，那么该关注的传输过程立即停止。

5）判断光子的生存状态1）返回步骤 1 继续模拟设置一个阈值，判断光子的权重是否小于该值，如果是大于该值的，则返回步骤 1 继续模拟2）重建概率模型经过多次碰撞后， 光子权重可能会变得非常小 如果直接终止权重小于设定阈值的光子的模拟过程，会对最后的计算结果产生误差，但如果不对这些权重特别小的光子进行处理，会耗费大量的计算时间为了解决这个问题， 引进了一种方差缩减技术， 当光子权重小于一定值时，重建概率模型，并随机产生新的光子权重，使期望值等于光子权重12min121min1/0PPWWPPPWw(9) 其中：w是新的光子权重模拟此概率模型随机产生新的光子权重：wWrwWrWwmin4min4min 0(10) 其中：4r是[0,1]之间的随机数由公式 (9)、(10)不难看出，两者的期望值保持不变， 但新的权重将权重值限制在以上，避免了计算大量的小贡献值，加快了收敛这种方法叫做俄罗斯轮盘赌重新构建光子权重后，如果权重为0，则结束此次模拟过程；如果不为 0，判断光子是否被接收或是否逃逸出烟雾范围，如果没有，则返回步骤 1 继续模拟6) 光子的接收当光子的状态满足一下条件时，该光子被探测器接收：1）光子的天顶角小于接收器的方位角与视场角一半的和，以及大于接收器的方位角与视场角一半的差时，光子进入接收器视场角；2）光子坐标 (px,py,pz) 与接收器坐标(rx,ry,rz)之间满足一下条件（模拟的是光子在y 轴方向的传播）：当 py=ry 时，如果222)()(rrzpzrxpx,则该光子被接收， 其中 r 是接收器的接收半径；当当 py>ry 时，将光子的坐标转换到垂直于ry 的坐标平面上，则：rzpyrypzrxpyrypx)()(其中：),,(是光子的方向向量。

当222)()(rrzpzrxpx时，该光子被接收三模拟结果及分析（1）多次散射对模拟结果的影响当光源的位置、光子的初始方向以及接收器的位置一定时，分析多次散射对模拟结果的影响模拟的光子数为20000，单次散射反照率分别为0.4、0.5、0.6、0.7、0.8 时，统计散射次数为1~6 次时，探测器接收到的信概率，结果如图1 所示12345-0.10.00.10.20.30.40.5接收到的信号的概率散射次数w=0.4w=0.5w=0.6w=0.7w=0.8图 1 多次散射接收的信号概率图 1 中，横坐标是散射次数，纵坐标是接收到的信号概率，折线分别表示单次散射反照率为0.4、0.5、0.6、0.7 和 0.8 时多次散射接收的信号概率由图可知，在其他模拟条件相同时， 单次散射反照率越大，接收器可以接收的光子数目越多；多次散射中，起主要作用的是一、二、三次散射，特别是一次散射，四次以后的散射几乎不起作用2）接收器位置对模拟结果的影响为了分析接收器位置对模拟结果的影响，光源的位置坐标为)5,0, 5(，方位角为 0 度，天顶角为 45 度；接收器的坐标为), 0,10(x，模拟过程中x的值取 5~10m，其他模拟条件相同；模拟的光子数为20000。

模拟结果如图 2 所示56789100.000.050.100.150.200.250.300.350.400.45接收的信号的概率接收器位置（m）图 2 接收器位置对模拟结果的影响图 2 中横坐标表示接收器的垂直高度， 纵坐标表示的是在不同高度上接收到的信号的概率由图可知，接收器的位置不同，对模拟结果的影响很大确定光源位置和发射角度后， 选取合适的探测器位置十分重要。