基于蒙特卡罗的光散射传输三维动态仿真*

张阳阳，高 隽，范之国，吴良海，高丽娟

(合肥工业大学 计算机与信息学院，安徽 合肥 230009)

天空中存在许多不同偏振度、不同偏振化方向的散射光，便形成了特定的包含大量偏振信息的偏振分布模式[1]，这种偏振分布模式蕴含丰富的方位信息和大气参数信息。这种偏振特性而主要是由于光在传输过程中与介质中粒子散射碰撞而产生的，研究光在散射介质中的传输和散射过程对偏振导航、探测都有重要的意义[2-3]。

光在散射介质中的传输是一个非常复杂的过程。典型的传输理论为辐射传输理论，该理论一般难以得到解析解。对散射介质中的光传输问题也可以用统计方法进行描述。目前解决光线在大气中的辐射传输问题比较有效的方法是蒙特卡罗方法，其思路清晰，能方便地应用到实际计算中。WILSON和ADAM首次将蒙特卡罗思想引入激光同组织相互作用方面[4]，随后，蒙特卡罗模拟方法就被广泛地应用到光在组织中的传输模拟。SEBASTIAN B应用斯托克斯-穆勒(Stokes-Mueller)形式模拟了偏振光在半无限介质中的传输，将蒙特卡罗模拟方法的应用扩展到了偏振领域[5]。JESSICA总结并给出了半无限介质中的偏振光传输的三种蒙特卡罗模拟方法[6]，蒙特卡罗方法得到的结果依赖于对大量光子的统计，由于模拟得到的结果与实验结果比较符合，因此蒙特卡罗方法被越来越多的科研工作者所采用。

本文采用蒙特卡罗统计方法模拟光束在散射介质中的传输和散射过程。同时，利用开源图形函数库OpenGL工具在VC++环境下设计实现了三维动态仿真，提供了三维空间中光散射传输的可视化系统，并对光散射传输后的偏振特性进行了分析。

1 偏振蒙特卡罗模型

英国物理学家Stokes于1852年提出利用一种四维向量(即Stokes矢量)，来描述光的偏振状态，由于其具有可加性，因此非常适用于对光线的分析工作。用Stokes矢量来描述光的偏振状态：S=[I Q U V]T,其中，I表示光术的总强度，±Q反映线偏振光水平或垂直的程度；±U反映线偏振光与Q成 ±45°的程度；±V反映光束中圆偏振光的左旋或右旋的程度。其偏振度可以表示为：

入射光Sin在介质层中传输时，输出光为Sout，可以用一个4×4的穆勒矩阵来描述其中的变化过程，有：

其中,πn和 τn是散射角度 α 的勒让德多项式，系数 an和bn是复折射率m=n+n′i(球形粒子相对于外界环境的折射系数)以及尺寸参数的函数。

偏振蒙特卡罗模型中粒子运动示意图如图1所示。其描述的是光子在半无限单一介质层中的散射传输过程，并用一定数量的光子依次循环来模拟现实光束散射的效果。

散射介质层可以看作是一个无限大但厚度有限的面板，在本文中厚度取4/us(us为散射系数)，介质层中为同一介质。如图1所示，沿着Z轴正方向从Z=0(即X-Y面)到Z=4/us面为介质层，可以认为在发射点 O有一个光子发射器沿着Z轴正方向发射光子，在Z=4/us面上放置一个正方形接收器(图中灰色区域所示)，本文中接收器的边长取14/us。为了较真实地反应一束光穿过介质层之后偏振态的变化情况，在本文中，如果没有特殊说明，取1e6个光子连续发射。接收器取100×100个方格点采样接收。

蒙特卡罗方法的算法流程图如图2所示，其中灰色步骤只有在偏振蒙特卡罗模型中需要。

图2 蒙特卡罗算法流程图

下面对蒙特卡罗方法的关键步骤做简要解释（可参考JESSICA总结的偏振蒙特卡罗子午面模型[6]）：

(1)初始参考平面为X-Z平面;

(2)启动时刻光的E矢量是在X-Z平面内定义的;

(3)根据已知的散射物质和拒绝方法的相位函数，随机确定角度α(散射角度)和角度β(入射平面和散射平面的夹角)。

(4)根据散射参数确定穆勒矩阵计算出射光;

(5)在光子散射出介质层后，若打在接收器上，进行累计至所有光子散射完。

另外，当光子其权重低于阈值时，光子运动终止。

2 光散射传输的三维动态仿真系统设计实现

按照上述的偏振蒙特卡罗模型，利用OpenGL图形库和VC++工具设计实现光在介质中散射传输的三维动态仿真系统，图3所示为仿真系统的功能框图。

仿真系统采用MFC框架的单文档结构，使用双线程来处理散射传输计算与可视化交互，界面线程负责光子散射传输过程的三维动态显示，结合OpenGL图形库，将设备描述符DC和渲染描述符RC联系起来，以完成显示功能；后台工作线程负责光子散射传输过程的路经计算，根据所述的偏振蒙特卡罗模型实现光子能量、散射方向、偏振态等信息的计算，两个线程通过消息队列相互通信，将后台计算的光子路径和接收器偏振信息显示在三维视图中。图4所示为仿真系统演示图。

利用此三维动态仿真系统可以直观地在三维空间中查看光子在介质中的散射传输过程和接收器的累积过程，且操作方便，可扩展性强。

3 系统仿真结果分析

利用所述仿真系统来模拟光束散射传输的过程，查看分析接收器上接收累积的光子偏振特性，即前向散射的偏振特性。参数为：入射光是自然光[1,0,0,0]，粒子半径 1.0 μm，光子波长 0.6 μm，粒子复折射系数(1.59+0i)。图5所示分别为接收器上 Stokes矢量（I,Q,U,V）的各个参数、偏振度DOP和Q分量的退偏振度DeP-Q。

从图5中可以发现，当入射光为自然光[1,0,0,0]时，经过介质散射传输后，强度I主要集中在中心区域且从中心向外扩散，这是符合一般常理的；对于Q分量和U分量，在相应区域均有一个线偏振较大的区域，其他区域基本趋于零；圆偏振强度V分量的分布比较散乱，无规律，且强度较低；对于偏振度OPP和Q分量的退偏振度DeP-Q分布基本上是以中心点对称的，且退偏振度较大的区域退偏现象较明显，故相应区域偏振度较低。

U分量和Q分量的分布结果呈一定的规律性，可以将其利用于偏振应用的相关方面。对于V分量，含量较少，且获取有一定的难度，故一般在偏振信息获取时可以将其忽略。结果中多数图像都是对称的，这些对称性包含了基于蒙特卡罗的光散射传输的一些偏振特性，值得后期更详细地研究。

文中采用斯托克斯-穆勒(Stokes-Mueller)形式描述偏振光散射传输过程，介绍了偏振蒙特卡罗模型，并以

此为基础，利用VC和OpenGL工具设计实现了三维

动态仿真系统，可以直观地在三维空间中查看光子在介质中的散射传输过程和接收器的累积过程，为光散射传输过程中的偏振特性分析提供了可视化工具。同时，对系统仿真的存储结果进行了分析，数据结果表现了偏振信息的对称性，具有一定的后续研究价值。

[1]COULSON K.Polarization and intensity of light in the atmosphere[M].Hampton VA:A.Deepak Pub，1988.

[2]黄旭锋，步扬，王向朝，等.基于米氏散射理论的太阳光散射偏振特性[J].中国激光，2010，37(12):3002-3006.

[3]吴良海，高隽，范之国，等.大气粒子散射特性及其对空间偏振分布的影响[J].光学学报，2011，31(7):27-33.

[4]WILSON B C，ADAM G.A Monte Carlo model for the absorption and flux distributions of light in tissue[J].Med Phys，1983，10(6):824-830.

[5]SEBASTIAN B，HIELSCHER A H.Monte Carlo simulations of the diffuse backscattering mueller matrix for highly scattering media[J].Appl Opt，2000，39(10):1580-1588.

[6]JESSICA R R，SCOTT P，STEVE J.Three Monte Carlo programs of polarized light transport into scattering media:part I[J].Opt.Express，2005，13(12):4420-4438.