0 引言

激光在介质中的散射特性是许多光学测量技术的基础，例如大气气溶胶中悬浮粒子成分的测量^［1］，生物医学中组织细胞特征的检测和识别^［2］，工业加工中微粒尺寸分布的质量监控等^［3］。通过对浑浊介质中散射光强的测量，可以获得散射颗粒的大小、形状和浓度等信息^［4-5］。当前用于介质中散射颗粒测量的理论存在一定的局限性，主要是因为介质中不可避免的多次散射效应难以量化描述。如何描述激光通过介质的传递情况，量化多次散射效应是当前的重要研究课题。考虑到描述光子迁移的辐射传递方程没有合适的解析解，基于统计的蒙特卡罗方法具有很强的灵活性，能考虑各种复杂的光源、介质散射特性以及检测器特性，可通过跟踪大量光子轨迹来得到与实际情况相符的结果。如徐兰青等^［6］通过蒙特卡罗方法研究了散射过程中逃逸出介质的后向散射光子数与散射次数的关系，得到了接收到的后向散射光子数随着光子散射次数的增加有先增大后减小的趋势。张硕等^［7］通过蒙特卡罗方法研究了在相同消光系数和光学厚度下，不同粒径的水雾对红外辐射的衰减效应，获得了当雾滴半径远大于辐射波长时，减小雾滴半径会加剧衰减效应的结果。上述研究讨论了各种颗粒物参数对散射光强大小的影响，没有得到二维散射光的分布信息，也没有具体实验组作为对比。管风等^［8］基于蒙特卡罗方法建立了水下激光在随机粗糙表面和水体悬浮粒子条件下的散射模型，但是由于计算效率不高仿真得到的二维图像相对离散与实际结果差别较大。此外，由于精确的Mie散射相位函数形式复杂，多数基于蒙特卡罗方法仿真光子散射的模型采用H-G相位函数来代替Mie散射相位函数，虽然可以提升计算效率，但是会对计算精度有一定影响。

本文通过实验记录了激光在不同光学厚度，不同粒径的均匀悬浮液中的侧向散射光强，根据实际参数建立了激光多次散射蒙特卡罗仿真模型，通过对Mie散射相位函数的累积概率密度函数插值的方法求取散射角。首先介绍了实验装置和原理，描述了激光在浑浊介质中传输的蒙特卡罗仿真模型，然后给出了实验记录和仿真中的激光侧向散射光强二维分布和衰减对比，分析比较了不同参数浑浊介质对结果的影响，最后根据仿真记录的光子散射次数讨论了不同散射阶次对侧向散射总光强的贡献。

1 实验装置及原理

激光束在均匀单分散散射介质中的侧向散射光强分布实验装置如图1，入射激光为532 nm激光器，功率为80 mW，出光口束腰半径为1 mm。激光通过6 cm×4 cm×5 cm的样品槽传输，槽中的散射介质为聚苯乙烯颗粒在蒸馏水中配置的均匀悬浮液，聚苯乙烯颗粒粒径为10 µm和5 µm。

图 1.

Fig. 1. Experimental setup

下载图片 查看所有图片

激光通过恒定消光系数的均匀介质，光强呈指数衰减，可表示为

I=I0⋅e-μel

式中，是通过散射介质前的光强，是光子运动的距离，为介质的消光系数，等于散射系数与吸收系数之和。光子在介质中与颗粒前后两次相互作用的平均距离可用平均自由程表示，与消光系数成反比。介质的光学厚度OD定义为光子在介质中与散射颗粒相互作用的平均次数^［9］，可表示为

OD=l/lfp=l⋅μe

可以看出，光学厚度为无量纲。当激光穿过固定距离的散射介质，介质的光学厚度与消光系数成正比。介质的消光系数与介质中单位体积所包含的颗粒数N成正比，可通过改变样品槽中聚苯乙烯颗粒的数量浓度来改变介质的光学厚度。

对于10 µm与5 µm两种粒径的聚苯乙烯颗粒，分别配置了光学厚度为2、5和8的样品进行实验。选用型号为Watec WAT-902H2的CCD相机，像素为768×576，镜头焦距f =12 mm。调整相机位置使得包含样品槽的64 mm×48 mm成像面清晰地呈现在视野正中^［10］，此时CCD像素的亮度值对应照射区颗粒在一个立体角范围内的散射光能。此时镜头距离激光束为9 cm，镜头的检测接收角为18.4°。在实验前对CCD响应进行标定，在无增益时散射光强与CCD像素的灰度值响应近似为线性。

2 浑浊介质中光子传输的蒙特卡罗模型

采用蒙特卡罗方法模拟激光在浑浊介质中的传输，通过将激光抽象成大量独立光子，追踪每个光子的路径^［11］。模拟步骤为：以激光与介质初次接触的表面中心为原点建立坐标系^［12］，假设光子从坐标原点沿Z轴正方向进入散射介质，光子初始坐标为（0，0，0），方向余弦为（0，0，1），如图2所示。

图 2.

Fig. 2. Monte Carlo simulation diagram

下载图片 查看所有图片

根据Beer-Lambert定律，光子每次与介质相互作用前运动的随机步长为

L=-ln(ξ)μe

式中，为0和1之间均匀分布的随机数。根据MA Xiaoyan等^［13］的实验结果，聚苯乙烯球在532 nm激光的照射下，折射率m近似为1.588+0.0i，因此认为在聚苯乙烯球中激光不发生吸收，只进行散射。由于介质的散射系数等于消光系数，介质的反照率为1，即光子的权重不随散射次数的增加而下降。

为了确定光子散射后的坐标，除了随机步长外，还需要知道光子的运动方向。光子的运动方向由散射角和方位角确定。散射微粒在方位角上表现为各向同性，取0到间均匀分布的随机数。散射角是光子入射方向与散射方向之间的夹角，可以由散射相位函数进行描述。散射相位函数的物理意义是入射能量经粒子散射后的能量角分布^［14］。根据Mie散射理论计算^［15］，直径10 µm和5 µm的聚苯乙烯颗粒球对于波长为532 nm激光的散射相位函数如图3所示。

图 3.

Fig. 3. Scattering phase function of a single polystyrene sphere (logarithmic scale)

下载图片 查看所有图片

由于Mie散射得到的散射相位函数无法解析表示，所以无法用求反函数的方法获得散射角^［16］。因此，在仿真中采用Mie散射相位函数的累积概率密度函数（Cumulative Probability Density Function ， CPDF）来描述。如图4所示，由于CPDF是单调递增的，虽然无法得到解析式，但其函数值和散射角有一一对应的关系。通过建立散射角与CPDF之间的查找表，就能够用0和1之间均匀分布的随机数对CPDF进行插值反求散射角，在仿真中控制散射角的分辨率为0.1°。

图 4.

Fig. 4. CPDF of Mie scattering phase function

下载图片 查看所有图片

综上所述，光子发生散射的位置可表示为

xn+1=xn+L×μxyn+1=yn+L×μyzn+1=zn+L×μz

式中，为光子当次散射前的坐标，为光子当次散射的方向余弦。通过再次计算得到散射角和方位角后，建立过渡坐标系，可得新的方向余弦为^［17］

μx1μy1μz1=sinθ(μxμzcosφ-μysinφ)1-(μz)2+μxcosθsinθ(μyμzcosφ+μxsinφ)1-(μz)2+μycosθ-sinθcosφ1-(μz)2+μzcosθ

若散射后的光子仍在介质中，则继续计算随机步长散射角和方位角得到下一个位置，重复此过程直到光子离开散射介质。对于每一个光子，记录其在离开介质前的散射次数，坐标以及最后的运动方向。最后根据物理光学理论设置合理的光子检测区域，若光子在接受视场内，则记录光子在检测区域的坐标。本实验对于每种情况的介质参数，每次试验5×10⁷个光子数量。

3 实验与仿真结果对比分析

对于用直径分别为10 µm和5 µm的聚苯乙烯颗粒配置的均匀介质，得到侧向散射光强二维分布如图5所示，左边两列表示10 µm数据，右边两列表示5 µm数据，结果图像保留靠近入射端16.78.8 mm包含散射光强信息的部分，对应X轴为237到342像素，Z轴为20到220像素，共200106像素。其中第一、第三列表示实验光强，第二、第四列表示仿真光强，对应的刻度表示每个像素检测到的光强与入射光强比值。从上往下各行的介质浓度分别与光学厚度为2、5和8时相对应，介质浓度依次增大。由于图像仅呈现16.7 mm介质厚度的侧向散射光情况，对应实际光学厚度值为0.56，1.39和2.22。激光沿着横坐标Z从左向右穿透介质，为了避免实验中入射端处样品槽壁反射对实验结果的影响，选择Z轴的起始刻度为距离实际入射端20像素处。对于10 µm粒径，实验与仿真得到的最大散射光强与入射光强比值在OD=0.56时为4.91×10^-5和3.44×10^-5，OD=1.39时为7.81×10^-5和5.24×10^-5，OD=2.22时为1.14×10^-4和8.74×10^-5。对于5 µm粒径，实验与仿真得到的最大散射光强与入射光强比值在OD=0.56时为7.36×10^-5和5.95×10^-5，OD=1.39时为1.08×10^-4和1.11×10^-4，OD=2.22时为1.29×10^-4和1.51×10^-4。可以看出对于相同粒径的介质，最大散射光强随着光学厚度的增大而增大。在相同光学厚度下，5 µm聚苯乙烯颗粒的散射光强要大于10 µm颗粒，这是因为在90^o附近的散射相位函数值5 µm颗粒要大于10 µm颗粒，与图3所示的散射相位函数相符合。由此可以看出，在消光系数相同的情况下，介质颗粒的粒径对散射光强分布有直接影响，导致这种情况的直接原因是不同粒径的介质颗粒在不同的角度上散射强度分布不同。对于所探究的不同介质参数，仿真与实验得到的最大散射光强与初始光强比均在同一数量级上，证明该蒙特卡罗仿真能够较好地符合实际情况。

图 5.

Fig. 5. Comparison of two-dimensional side scattered light intensity distribution in experiment and simulation

下载图片 查看所有图片

为了进一步比较散射光强在Z轴方向的分布情况，统计了图5中各光强分布图中间3行沿Z轴的平均光强大小，结果如图6。从图6中可以看出，介质的光学厚度越大，散射光强也越大，散射光强沿着Z轴方向不断减小。由于实验样品槽中颗粒物的不均匀和环境反射杂光的干扰，实验得到曲线有明显的振荡现象。仿真模拟则是在理想环境下，采用了足够多光子数，所以显得平滑。为了衡量散射光强减小的程度，统计了图6各条光强曲线前10个像素与末尾10个像素的平均光强并计算它们的差值。对于10 µm粒径，实验与仿真沿Z轴减小的光强与入射光强的比值在OD=0.56时为1.35×10^-5和1.28×10^-5，OD=1.39时为3.33×10^-5和2.98×10^-5，OD=2.22时为4.37×10^-5和4.24×10^-5。对于5 µm粒径，实验与仿真沿Z轴减小的光强与入射光强的比值在OD=0.56时为2.26×10^-5和2.64×10^-5，OD=1.39时为5.13×10^-5和5.78×10^-5，OD=2.22时为6.06×10^-5和6.31×10^-5。可以看到散射光强的衰减也随着光学厚度的增大而增大。光学厚度越大，光子与介质中散射颗粒作用的平均次数就越多，初始散射光强越大，同时散射光强衰减的程度也越大。定义光强衰减的百分比差异为实验与仿真衰减光强之差的绝对值与衰减光强平均值之比。对于10 µm粒径，实验与仿真光强衰减的百分比差异在OD为0.56，1.39和2.22时分别为5.2%，10.8%和3.0%。对于5 µm粒径，实验与仿真光强衰减的百分比差异在OD为0.56，1.39和2.22时分别为15.8%，12.0%和4.1%。在相同介质条件下，实验与仿真得到的散射光强衰减百分比差异均小于16%，且散射光强曲线走势基本一致，表明该仿真能够较好地反映实际实验。

图 6.

Fig. 6. Comparison of experimental and simulated one-dimensional light intensity curves

下载图片 查看所有图片

在仿真中通过记录每个光子到达检测区域前经历的散射次数，进一步得到各个散射阶次检测的光强大小，即可量化各个散射阶次对最后总光强的贡献，结果如图7所示。对于10 µm粒径的聚苯乙烯颗粒，在光学厚度分别为2、5和8的介质中，侧向散射的最大阶次分别为2次、4次和7次，对应概率分别为33.1%，24.5%和17.3%。对于5 µm粒径的聚苯乙烯颗粒，在相应光学厚度的介质中，侧向散射的最大阶次也分别为2次、4次和7次，对应概率分别为32.8%，22.4%和16.4%。可以看到对于不同粒径的介质颗粒，侧向散射的最大阶次均与光学厚度相等或相近。

图 7.

Fig. 7. Scattering orders of particles with different sizes under different optical thicknesses

下载图片 查看所有图片

在相同的光学厚度下，不同粒径聚苯乙烯颗粒的各个散射阶次的分布几乎相同。这说明激光在介质中的散射次数分布由介质的光学厚度决定，与颗粒粒径大小无关。同时，还能看到随着光学厚度的增大，多次散射明显增强，主要散射阶次的贡献不断下降，散射次数分布曲线趋于平滑。结合散射次数分布与光强分布，可以看到多次散射对光强增大影响显著，必须加以考虑。

4 结论

通过实验探究了激光在均匀单分散聚苯乙烯悬浮液中的侧向散射光强分布情况，研究了不同介质参数对实验结果的影响。通过改变浑浊介质的浓度使其光学厚度分别等于2、5和8，使用5 µm和10 µm两种粒径来改变介质颗粒的散射相位函数，建立了激光在浑浊介质中传输的蒙特卡罗仿真模型。在所有研究的不同介质参数条件下，蒙特卡罗仿真与实验得到的侧向散射光强二维分布一致性良好。仿真与实验得到的最大散射光强与入射光强的比值和散射光强的衰减均随光学厚度的增大而增大，相同介质条件下散射光强的衰减百分比差异均小于16%。仿真还能提供多次散射情况下不同散射阶次的分布情况，表明该模型能够准确分析多次散射对结果的影响。