光学学报, 2023, 43 (18): 1829001, 网络出版: 2023-09-14   

基于偏振蒙特卡罗法的散射环境前向传输优化算法【增强内容出版】

Forward Transmission into Scattering Media by an Improved Polarized Monte Carlo Program
作者单位
1 鲁东大学交通学院,山东 烟台 264025
2 山东德衡人力资源股份有限公司,山东 烟台 264025
摘要
偏振蒙特卡罗法在模拟大量光子传输时不存在由有限离散引起的计算误差,因而通常作为验证其他方法计算精度的标准。通常情况下,偏振保持特性良好的光子的信息损耗少且传输距离远,分析这些光子的传输特性对于提取优良信号具有潜在的应用价值。为此,基于偏振子午面蒙特卡罗法,提出一种散射环境中前向传输光子偏振态的统计方法。模拟计算偏振光在聚苯乙烯悬浊液中的传输特性,结果表明,优化算法不仅能够反映光子偏振态的变化情况,而且能够统计偏振保持特性良好的光子数占比。
Abstract
Objective

Vector radiation transport in scattering media is a research hotspot. The exact analytical solution to the vector radiative transfer equation cannot be obtained without the simplified treatment, and it needs to be calculated by numerical calculation methods. The polarized Monte Carlo program simulates the transmission of mass photons. As it does not lead to computational errors due to finite dispersion, the program is usually employed as a standard to verify the computational accuracy of other methods. At present, this method can obtain the outgoing polarization state of each light wave after the solution is obtained, but it is difficult to know the photon transmission situation. This limits the analysis of the polarization state retention of scattering light. However, photons with good retention characteristics of the polarization state usually have a small information loss and a long transmission distance. The study and analysis of these photons can be potentially applied to extract good transmission signals.

Methods

This paper proposes a method to count the photon polarization states during forward transmission into scattering media on the basis of the polarization meridian Monte Carlo program. The original algorithm and improved algorithm are shown in Fig. 1. The white part is required for calculating both the original algorithm and the improved algorithm, and the blue part is the additional flow for calculating the improved algorithm. One hundred thousand parallelly polarized photons or right-handed circularly polarized photons are sent into a slab represented by one particular particle distribution for each environment, and photons are transmitted at a given distance. Then, the aggregated polarization is calculated from the photons that arrive in a given area, and photons on the front face of the slab are considered the transmitted photons. This process continues for all the launched photons, and the result is calculated. For the original algorithm, a photon completes forward transmission when it passes through a specific forward distance. The improved algorithm adds a link to output the polarization state of each photon. Moreover, the improved algorithm can count the total number of received photons and the number of photons similar to the polarization state of the initial photon.

Results and Discussions

This paper uses an example to compare the original algorithm and the improved algorithm. The following simulations are performed in the polystyrene suspension with a mass concentration of 2.08 μg/μm3. The wavelength of incident light is 532 nm, and polystyrene's refractive index is 1.597. The particle diameter of the polystyrene suspension is 1 μm in simulation, and the transmission distance is 10 cm. One hundred thousand parallelly polarized photons or right-handed circularly polarized photons are launched for simulations of the original algorithm and the improved algorithm separately. There are 66358 received photons after the transmission of parallelly polarized photons and 66367 after the transmission of right-handed circularly polarized photons. After that, the first 100 received photons are selected as samples (Figs. 3 and 4). Calculations show that for the sample data of parallelly polarized photons, the NRoPS is 0.13, and for the sample data of right-handed circularly polarized photons, the NRoPS is 0.53. The calculation results are both similar to the overall situation. The simulated calculation of polarized light transmission in the polystyrene suspension demonstrates that the optimized method can not only reflect the change in the photon polarization state but also count the percentage of photons with good retention characteristics of the polarization state. The comparison between the original results and the optimization results shows that the results of the optimized algorithm are less than the polarization degree value. This is because the optimized algorithm not only avoids the error introduced in the calculation of the intensity difference in the orthogonal component but also excludes photons that are not similar to the polarization state of the initial photon.

Conclusions

Compared with the polarized Monte Carlo program, the improved method can not only reflect the change in the photon polarization state but also count the percentage of photons with good retention characteristics of the polarization state. It reflects the change in the polarization state of the transmitted photons in multiple dimensions. This study can provide technical support for research on the extraction of excellent transmission signals.

1 引 言

诸如云、雾、霾、烟、沙尘、浑浊水域等散射环境普遍存在1-3,光在散射环境中传输时受到介质的吸收、散射等作用,表现出光强度、频率、传播方向、偏振等性质改变4。构造辐射传输方程能够描述光在散射介质中的传输特性,传统求解方法为规避求解复杂的矢量辐射传输方程,通常将光当作标量处理,但近似简化不仅造成了计算精度下降,而且丢失了偏振维度信息。例如,我国第一颗海洋水色卫星(HY-1A)采用标量辐射理论近似算法处理大气校正,导致其信息提取精度和应用范围受到限制5。由于标量辐射传输理论所构建的仪器已趋于技术极限,偏振技术逐渐发展起来。众多生物、理论和应用案例证明,散射偏振传输特性是丰富且敏感的6-9,相关研究有助于提高信号质量10-11、增加传输距离12乃至反演环境属性13

偏振光具有矢量性,而矢量计算比标量计算更为复杂,因此在散射偏振传输计算问题中矢量辐射传输方程在没有简化处理时无法得到精确解析解,需要借助数值计算方法进行计算。目前,发展较为成熟的数值计算方法包括离散坐标法14、球谐函数法15、倍加累加法16、蒙特卡罗法、二流近似法17-18等。其中,蒙特卡罗法19-20模拟大量光子统计传输规律,通过增加模拟光子数量提高计算精度,不存在由有限离散引起的计算误差,它通常作为验证其他方法计算精度的标准。

偏振蒙特卡罗法通过统计大量光子获取出射光的偏振态,继而统计不同类型入射光的出射偏振态,即可求解散射环境的穆勒矩阵,但该算法在求解过程中难以知晓光子的传输情况,这就限制了其对散射偏振保持情况的分析。为此,本文提出一种基于偏振子午面蒙特卡罗法的优化算法,并在优化算法中增加输出每个光子偏振态的环节,能够统计与初始光子偏振态相近的光子数和总接收光子数。所提优化算法不仅能够反映散射环境中光子前向传输时偏振态的变化情况,而且能够统计偏振保持特性良好的光子数占比。

2 所提优化算法构建

所提优化算法基于偏振子午面蒙特卡罗法19-20,流程如图1所示。所提优化算法在偏振子午面蒙特卡罗法的基础上增加了输出每个光子偏振态的环节,并统计总接收光子数和与初始光子偏振态相近的光子数。

图 1. 基于偏振子午面蒙特卡罗法的优化算法

Fig. 1. Improved algorithm of polarized meridian plane Monte Carlo program

下载图片 查看所有图片

所提优化算法的计算步骤为:①发射光子,此时需要定义光子的初始参考平面,确保微粒与散射矩阵的相互作用在同一参考系中。②光子通过散射介质。由于吸收作用,光子的权重因子会有所下降;由于散射作用,在单次散射发生后,参考平面需要重新定义,以确保偏振光进入正确的子午平面。散射后光子的电场旋转到新的子午平面。③判定光子是否达到死亡条件。④判定光子是否穿过介质,若穿过介质,调整子午线平面,输出每个光子的偏振态,否则循环步骤②和③。⑤判定光子是否与初始光子偏振态相近,若相近,则另置输出,统计计数。⑥利用“轮盘赌”算法,所有光子完成传输后终止算法运行。算法中的光子并非现实中的光子,实为能量可约的小光束,这些小光束在传输过程中,能量会有所衰减,而探测器有阈值限定,如果小光束能量低于探测器的阈值,则判定光子消亡。算法采用权重因子判断能量的衰减,权重因子被分配给小光束,初始值为1,小光束每次散射后进行更新,并在整个蒙特卡罗模拟中进行跟踪,当权重因子降低到一定阈值,设定小光束被完全吸收。

子午面转换如图2所示。从图2(a)可以看到,光子的电场E最初是相对于子午面COA定义的。在散射事件发生之前,电场E需通过旋转矩阵Rψ)旋转到平面BOA图2(b)],以确保微粒与散射矩阵的相互作用在同一参考系中,其中ψ为面COA到面BOA的旋转角,旋转方向为逆时针。从图2(c)可以看到,光子移动并发生散射,光子位于散射参考平面BOA,散射后光子从A点移动到B点,微粒的散射作用通过散射矩阵Ma)转换,其中a为散射角。从图2(d)可以看到,当单次散射发生后,参考平面需重新定义以确保偏振光进入正确的子午平面,散射后光子的电场E旋转到新的子午平面COBΘ为面BOA到面COB的旋转角,此时旋转方向为顺时针,旋转矩阵为R(-Θ)。单次散射后光子的Stokes矢量为

Sout=R(-Θ)M(a)R(ψ)Sin

式中:Sin为初始入射光子的Stokes矢量。

图 2. 子午面转换图。(a)光子电场E最初是相对于子午面COA定义的;(b)在散射发生之前,电场E旋转到参考散射平面BOA;(c)光子在散射平面BOA发生散射;(d)散射后,电场E旋转到参考散射平面COB

Fig. 2. Diagrams of meridian plane conversion. (a) E is initially defined with respect to the meridian plane COA; (b) before scattering, E must be referenced to the scattering plane BOA; (c) scattering at the scattering plane BOA; (d) after scattering, the reference plane is aligned to a new meridian plane COB

下载图片 查看所有图片

所提优化算法统计与初始光子偏振态相近的光子数和总接收光子数,并不统计与初始光子偏振状态一致的光子数,这是因为散射偏振前向传输过程中会出现弹道光子、蛇形光子和散射光子21,仅有少部分弹道光子能够与初始光子偏振态一致,统计时占比很小。秉持偏振保持特性良好的光子的信息损耗通常较少的原则,权衡增加统计光子数和尽量规避偏振保持特性劣质的光子数的实现难度,将与初始光子偏振态相近的接收光子作为统计对象。偏振态相近定义为接收的光子类型与初始发射的光子类型相近,相近程度依据接收光子的Stokes参量变化程度判定,判定标准视具体情形调整。例如,初始光子为完全水平线偏振光,Stokes参量S1a,归一化后s1=S1/S0=1,将接收光子的s1是否在[0.97,1.00]区间作为判定条件,当接收光子的s1在[0.97,1.00]区间时,判定为相近,否则判定为不相近。

3 模拟计算

选用聚苯乙烯悬浊液为散射环境,它由粒径为1 μm的聚苯乙烯微粒与水混合而成,质量浓度为2.08 μg/cm3。模拟时采用高纯单分散的颗粒分布。聚苯乙烯的折射率与波长的关系见文献[22],模拟时选用的激光波长为532 nm,聚苯乙烯的折射率为1.597。模拟水平线偏振光和右旋圆偏振光在聚苯乙烯悬浊液中前向传输,传输距离设定为10 cm。初始入射光为同一偏振态光子的集合,分别模拟传输100000个水平线偏振光子和右旋圆偏振光子,最终接收的水平线偏振光子为66358个,接收的右旋圆偏振光子为66367个。

由于总体数据容量较大,本实验进行取样分析。为规避接收光子能量不均的干扰,采样之前,对接收光子偏振态进行归一化处理。采用能量比值的方式规避能量的干扰,计算过程则无需考虑光子能量变化以及光子能量不可约的因素。为此,对Stokes矢量进行归一化处理。

任意光的偏振态都可用4个Stokes参量构成的Stokes矢量表示,即

表 1. 右旋圆偏振光子的归一化Stokes参量

Table 1. Normalized Stokes parameters of right-handed circularly polarized photons

Numbers0s1s2s3
110.09780.06030.9934
210.07970.07340.9941
310.0493-0.12960.9940
41-0.51370.45820.7253
510.1399-0.14230.9799
61-0.43920.79970.4093
71-0.10390.86710.4870
810.01580.00970.9998
910.19430.14820.9697
101-0.1985-0.09570.9754
111-0.06230.30110.9516
1210.2486-0.35470.9013
131-0.0314-0.28580.9577
141-0.2427-0.32310.9147
1510.51210.04080.8579

查看所有表

S=S0S1S2S3

式中:S0为光的总强度;S1为水平线偏振光与垂直线偏振光的光强差值;S2为45°线偏振光与-45°线偏振光的光强差值;S3为右旋圆偏振光与左旋圆偏振光的光强差值。

式(2)归一化处理,得到

S=S0S0    S1S0    S2S0   S3S0T=1   s1   s2   s3T

式中:s1为水平线偏振光和垂直线偏振光的光强差值占总光强的比例;s2为45°线偏振光与-45°线偏振光的光强差值占总光强的比例;s3为右旋圆偏振光与左旋圆偏振光的光强差值占总光强的比例。

随机截取右旋圆偏振光传输结果作为样本,其归一化处理结果如表1所示。初始光子为右旋圆偏振光,接收光子中未有与初始光子偏振状态一致的接收光子。归一化Stokes参量s3越接近1,接收光子的类型与初始发射光子的类型越相近。对于1、2、3、5、8、9、10号接收光子,其归一化Stokes参量s3在[0.965,1.000]区间,s1s2均在[-0.2,0.2]区间,转化为偏振度(PDoP)和线偏振度(PDoLP),分别得到

PDoP=S12+S22+S32S0=s12+s22+s33PDoLP=S12+S22S0=s12+s22

PDoLP/PDoP结果分布区间为[0,0.08],可见线偏振光对偏振度的干扰较小。因此,设定归一化Stokes参量s3是否在[0.965,1.000]区间为判定接收光子与初始光子相近的条件,当接收光子的s3在[0.965,1.000]区间时,判定为相近,否则判定为不相近。该区间可视情况调整。

水平线偏振光子和右旋圆偏振光子传输后,分别选取前100个接收光子作为样本。水平线偏振光的传输样本如图3所示,其中Nphoton表示第N号光子。设定接收光子与初始光子相近的指标为归一化Stokes参量s1在[0.965,1.000]区间。统计每个接收光子的s1分布情况,采用柱状图描述整体样本态势,采用雷达图描述每个接收光子的态势。由图3可知,水平线偏振光传输后,光子归一化Stokes参量s1发生明显变化,s1值在[0.965,1.000]区间的水平线偏振光子为13个。

图 3. 水平线偏振光子传输样本分布。(a)柱状图;(b)雷达图

Fig. 3. Sample distribution of horizontal polarization photons. (a) Histogram; (b) radar chart

下载图片 查看所有图片

右旋圆偏振光传输样本如图4所示。设定接收光子与初始光子相近的指标为归一化Stokes参量s3在[0.965,1.0000]区间,统计每个接收光子的s3分布情况。由图4可知,右旋圆偏振光传输后,光子归一化Stokes参量s3变化缓和,s3值在[0.965,1.000]区间的水平圆偏振光子为53个。

图 4. 右旋圆偏振光子传输样本分布。(a)柱状图;(b)雷达图

Fig. 4. Sample distribution of right-handed circularly polarization photons. (a) Histogram; (b) radar chart

下载图片 查看所有图片

为描述前向散射光的偏振态变化,先前研究23-24定义了新参量——RoPS,其定义为偏振状态保持率。当以光子形式传输时,RoPS的表达式为

NRoPS=±NζNtotal

式中:Ntotal表示发生前向散射的全部光子数;Nζ表示ζ型光子的数量。ζ型光子与初始入射光子具有相近的偏振态,±表示ζ型光子的类型,其判定标准与Stokes参量判定标准一致。

求解得到水平线偏振光传输样本的NRoPS=0.13,将样本容量扩展至总体,求得NRoPS=0.14;求解得到右旋圆偏振光传输样本的NRoPS=0.53,采用相同的计算方式,将样本容量扩展至总体,求得NRoPS=0.55。样本计算结果与总体计算结果相近。

4 结果分析

图5所示为所提优化算法与原算法的比较。入射光子经由散射介质后,原算法的输出结果为Stokes矢量;优化算法统计总接收光子数和与初始光子偏振态相近的光子数。

图 5. 原算法与改进算法的比较

Fig. 5. Comparison between original algorithm and improved algorithm

下载图片 查看所有图片

采用偏振子午面蒙特卡罗法计算出射光的偏振态。初始入射光为同一偏振态光子的集合,分别模拟传输100000个水平线偏振光子和右旋圆偏振光子。表2所示为出射光的Stokes矢量。

表 2. 出射光的Stokes矢量

Table 2. Stokes vectors of output light

LightHorizontal linearly polarized lightRight circularly polaried light
Stokes vector0.65000.4148-0.0038-0.00240.6502-0.00180.00210.5448

查看所有表

求解得到水平线偏振光传输后的线偏振度PDoLP=0.4073,右旋圆偏振光传输后的圆偏振度PDoCP=0.8379。比较优化结果可知,优化算法的偏振度数值均小于原算法的结果。这是因为优化算法不仅规避了正交分量光强度差值计算引入的误差,而且规避了将与初始光子偏振状态不相近的光子纳入计算。

此外,原算法输出结果仅能解算整体出射光的偏振态,结果由所有接收光子会聚产生,而每个接收光子的偏振状态却无法观察,这限制了对散射光偏振保持情况的分析。优化算法增加了输出每个光子偏振态的环节,并且统计与初始光子偏振态相近的光子数和接收的总光子数。优化算法不仅能够反映光子偏振态的变化情况,而且能够分析偏振保持特性良好的光子数占比,这从多维度反映了传输光子偏振态的变化情况。

5 结 论

提出一种基于偏振子午面蒙特卡罗法的优化算法,该算法能够反映单个输出光子的偏振态,并且能够统计与初始光子偏振态相近的光子数和总接收光子数。优化算法不仅能够反映光子偏振态的变化情况,而且能够统计偏振保持特性良好的光子数占比。这弥补了求解过程中难以知晓光子传输情况的缺陷。模拟完全水平线偏振光和完全右旋圆偏振光在聚苯乙烯悬浊液中的前向传输特性,通过算法对比,验证了优化算法的可行性,该算法可为优良传输信号的提取研究提供新的思路。

参考文献

[1] 厉祥, 张然, 路淑芳, 等. 基于米散射蒙特卡罗法的紫外偏振透云研究[J]. 激光与光电子学进展, 2021, 58(17): 1701001.

    Li X, Zhang R, Lu S F, et al. Ultraviolet polarization employing Mie scattering Monte-Carlo method for cloud-based navigation[J]. Laser & Optoelectronics Progress, 2021, 58(17): 1701001.

[2] 徐俊杰, 卜令兵, 刘继桥, 等. 机载高光谱分辨率激光雷达探测大气气溶胶的研究[J]. 中国激光, 2020, 47(7): 0710003.

    Xu J J, Bu L B, Liu J Q, et al. Airborne high-spectral-resolution lidar for atmospheric aerosol detection[J]. Chinese Journal of Lasers, 2020, 47(7): 0710003.

[3] 唐军武, 朱培志, 刘秉义, 等. 海洋剖面激光雷达探测中颗粒物偏振散射问题[J]. 光学学报, 2022, 42(12): 1200001.

    Tang J W, Zhu P Z, Liu B Y, et al. Polarized light scattering of particulate matter in detection of oceanographic profiling LiDAR[J]. Acta Optica Sinica, 2022, 42(12): 1200001.

[4] 刘红林. 关于散射成像研究现状的一些思考[J]. 红外与激光工程, 2022, 51(8): 20220261.

    Liu H L. Some thoughts about the current research situation of imaging through scattering media[J]. Infrared and Laser Engineering, 2022, 51(8): 20220261.

[5] 何贤强, 潘德炉. 海洋-大气耦合矢量辐射传输模型及其遥感应用[M]. 北京: 海洋出版社, 2010.

    HeX Q, PanD L. Ocean-atmosphere coupling vector radiation transmission model and its remote sensing application[M]. Beijing: Ocean Press, 2010.

[6] LernerA, ShasharN. Polarization vision in cephalopods[M]//Horváth G. Polarized light and polarization vision in animal sciences. Springer series in vision research. Heidelberg: Springer, 2014, 2: 217-224.

[7] ZengX W, ChenX Y, LiY H, et al. Polarization enhancement of linearly polarized light through foggy environments at UV-NIR wavelengths[J]. Applied Optics, 2021, 60(26): 8103-8108.

[8] 郭忠义, 汪信洋, 李德奎, 等. 偏振信息传输理论及应用进展[J]. 红外与激光工程, 2020, 49(6): 20201013.

    Guo Z Y, Wang X Y, Li D K, et al. Advances on theory and application of polarization information propagation[J]. Infrared and Laser Engineering, 2020, 49(6): 20201013.

[9] HouW Z, LiZ Q, WangJ, et al. Improving remote sensing of aerosol microphysical properties by near-infrared polarimetric measurements over vegetated land: information content analysis[J]. Journal of Geophysical Research: Atmospheres, 2018, 123(4): 2215-2243.

[10] SankaranV, WalshJ T, Jr, MaitlandD J. Polarized light propagation through tissue phantoms containing densely packed scatterers[J]. Optics Letters, 2000, 25(4): 239-241.

[11] van der LaanJ D. Evolution and persistence of circular and linear polarization in scatting environments[D]. Tucson: University of Arizona, 2015.

[12] FadeJ, PanigrahiS, CarréA, et al. Long-range polarimetric imaging through fog[J]. Applied Optics, 2014, 53(18): 3854-3865.

[13] LiuJ, LiuJ H, HeX Q, et al. Retrieval of marine inorganic particle concentrations in turbid waters using polarization signals[J]. International Journal of Remote Sensing, 2020, 41(13): 4901-4922.

[14] CollinC, PattanaikS, LiKamWaP, et al. Discrete ordinate method for polarized light transport solution and subsurface BRDF computation[J]. Computers & Graphics, 2014, 45: 17-27.

[15] TapimoR, KamdemH T T, YemeleD. A discrete spherical harmonics method for radiative transfer analysis in inhomogeneous polarized planar atmosphere[J]. Astrophysics and Space Science, 2018, 363(3): 52.

[16] 张颖, 张熠, 赵慧洁. 多种天气条件下的天空光偏振模型[J]. 红外与毫米波学报, 2017, 36(4): 453-459.

    Zhang Y, Zhang Y, Zhao H J. A skylight polarization model of various weather conditions[J]. Journal of Infrared and Millimeter Waves, 2017, 36(4): 453-459.

[17] MarkelV A. Two-stream theory of light propagation in amplifying media[J]. Journal of the Optical Society of America B, 2018, 35(3): 533-544.

[18] YinQ, SongC. Fundamental definition of two-stream approximation for radiative transfer in scattering atmosphere[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 2003614.

[19] Ramella-RomanJ C, PrahlS A, JacquesS L. Three Monte Carlo programs of polarized light transport into scattering media: part I[J]. Optics Express, 2005, 13(12): 4420-4438.

[20] Ramella-RomanJ C, PrahlS A, JacquesS L. Three Monte Carlo programs of polarized light transport into scattering media: part Ⅱ[J]. Optics Express, 2005, 13(25): 10392-10405.

[21] AlfanoR R, WangW B, WangL M, et al. Light propagation in highly scattering turbid media: concepts, techniques, and biomedical applications[J]. Photonics: Biomedical Photonics, Spectroscopy, and Microscopy, 2015: 367-412.

[22] MaX Y, LuJ Q, BrockR S, et al. Determination of complex refractive index of polystyrene microspheres from 370 to 1610 nm[J]. Physics in Medicine and Biology, 2003, 48(24): 4165-4172.

[23] 曾祥伟. Mie散射环境中光前向传输偏振特性的关键技术研究[D]. 大连: 大连理工大学, 2020.

    ZengX W. Key technology for polarization characteristics of light forward transmission in Mie scattering environments[D]. Dalian: Dalian University of Technology, 2020.

[24] ChuJ K, WuQ M, ZengX W, et al. Forward transmission characteristics in polystyrene solution with different concentrations by use of circularly and linearly polarized light[J]. Applied Optics, 2019, 58(25): 6750-6754.

曾祥伟, 张燕, 杨钧秀. 基于偏振蒙特卡罗法的散射环境前向传输优化算法[J]. 光学学报, 2023, 43(18): 1829001. Xiangwei Zeng, Yan Zhang, Junxiu Yang. Forward Transmission into Scattering Media by an Improved Polarized Monte Carlo Program[J]. Acta Optica Sinica, 2023, 43(18): 1829001.

本文已被 2 篇论文引用
被引统计数据来源于中国光学期刊网
引用该论文: TXT   |   EndNote

相关论文

加载中...

关于本站 Cookie 的使用提示

中国光学期刊网使用基于 cookie 的技术来更好地为您提供各项服务,点击此处了解我们的隐私策略。 如您需继续使用本网站,请您授权我们使用本地 cookie 来保存部分信息。
全站搜索
您最值得信赖的光电行业旗舰网络服务平台!