2.2 海水信道光传输散射特性的蒙特卡洛模拟

蒙特卡洛模拟法在分析和研究海水信道光传输特性中的应用最为广泛[16-18]。本节详细介绍光子在海水介质中随机传输的蒙特卡洛模拟法，并且对蒙特卡洛模拟法的散射相位函数进行分析，最后根据蒙特卡洛模拟过程对影响接收光功率的因素进行仿真，并对仿真结果进行分析。

蒙特卡洛模拟法对研究光在海水中的传输情况，以及对海水信道特性的分析与研究具有重要的意义。由于海洋环境是非常复杂的，所以光在水下进行传输时，所能涉及的不同海水水质参数具有非常大的随机性。光在海水传输的过程中会受到多次散射的影响，几乎不可能完全用解析法来描述光信号在水下的传输过程。在实验室的环境下，海水中所含的不同杂质成分和这些成分的衰减系数等因素的选择会受到实验条件的限制，因此要完全通过实验的方法得到所有数据几乎是不可能的。

光子在水下发生散射的过程，可看作是光子在水中的随机碰撞，这个随机运动过程可用蒙特卡洛模拟法来进行模拟。同时蒙特卡洛模拟法具有以下几个优点[19]。

① 算法设计简单，精度可控。

② 可以通过模拟不同发射机和接收机的特性，对海水不同介质的边界条件和散射条件进行模拟，并且将发散角、光子在碰撞过程中发生的散射等一系列影响因素全部加以考虑。

③ 容易确定光子散射的平均路径。

蒙特卡洛模拟法最大的不足之处就是计算时间太长，且单个光子的运动特性容易对其造成影响，从而在一定时间内发射的光子被接收机接收到的概率很小，所以为了得出比较准确的结果，通常需要模拟大量光子的运动轨迹。

根据光的粒子性，通常将光脉冲当成由很多个光子组成的光子包，所以光脉冲的传输问题就成了光子的运动问题[20-21]。由于海水中颗粒性质的不同，其在碰撞过程中将会被吸收或发生散射。若光子在传输过程中被吸收，则其传输过程就会终止；若光子发生散射，则其运动方向就会发生改变。

光子在海水运动的过程中，经常与粒子发生吸收性或者散射性的随机碰撞[22-24]。根据随机碰撞原理，构建蒙特卡洛模拟法的模型，并利用大量的光子进行模拟试验，则光束在海水介质中的传输规律就是对这些光子在运动过程中状态的统计。

海水中粒子的散射主要分为两类，第一类为瑞利散射，第二类为米氏散射。海水水分子的散射遵循的是瑞利散射规律，而其他颗粒的散射遵循的是米氏散射规律。

2.2.1 海水信道中的光传输散射效应

通常情况下，对海水介质中颗粒散射特性的分析，都是相对于单个粒子来说的。为了简便起见，通常利用米氏散射理论，将粒子与光子间的碰撞转化为球形粒子间的随机碰撞。根据海水介质颗粒密度的不同，一般把散射效应归纳为4类，分别是单次散射、一级多次散射、多次散射和漫射。分析发现，多次散射综合了单次、二次、高次散射以及在散射路径上的衰减，适用于大多数海水介质，所以它是分析和研究光子在海水介质中传输时最典型的一类散射。多次散射的效应分布如图2-10所示。

如图2-10所示，光子在海水介质中传输的过程中，会遇到很多颗粒而发生散射，因此对于非散射部分的直射光会变得很少。图2-10所示的多次散射光指的是传输过程的光发生散射后会偏离光轴，但是经过多次散射后部分光子又再次进入光轴中的光。

光子在海水介质中传输时，会发生多次散射，这会导致光能量发生衰减，如图2-11所示。

如图2-11所示，当光束入射到海水介质中时，会在各个方向上发生散射，由于海水具有强烈的前向散射特性，所以大部分的光会在光传输方向上发生散射。但是光束在海水中的传输距离增加时，直射的光束会变得越来越少，当光的传播路径非常长时，传播方向上未发生碰撞而直接到达接收面的光很少，所以海水中的光束基本上都是多次散射光。另一方面，即使散射光会发生散射，接近光轴中心处的光能量也仍然很大。因此，光束在传输过程中，会在时间和空间上发生扩展，而发生多次散射后，光束的能量在传输过程中也会发生衰减。

2.2.2 海水信道的光传输散射相位函数

散射相位函数能准确模拟出影响辐射传输的特性，同时被用来比较在不同介质中，散射角度的分布状况。

不含有任何杂质的海水为纯水，纯水水分子发生的散射服从的是瑞利散射，其散射相位函数如式（2-14）所示。

其中，θ表示散射角。

由于海水介质中含有大量悬浮颗粒，所以光束在海水中传输时一般都服从米氏散射。米氏散射相位函数P(θ)一般常用亨利－格林斯坦（Henyey-Greentein，H-G）相位函数[25-26]表示。H-G相位函数的优势是表达式相对简便，求解非常方便，模拟近似结果也比较好，并且可以更好地体现米氏散射相位函数前向峰值的特点。H-G相位函数是将散射角θ和散射因子g作为参数的函数表达式，如式（2-15）所示。

其中，表示非对称因子。g=0时，H-G相位函数表示的是散射各向同性；g=0.9时，表示的是强前向散射。Mobley[8]的分析结果表明，当θ非常接近0°和180°时，该函数存在一定的误差，当选取合适的g时，拟合H-G相位函数的匹配度在90%以上。当θ接近π时，P(θ)与介质实际的散射相位函数相比较小。这样光子发生后向散射的概率就会减小。这不但增加计算时间，而且还会降低精度。

因此可以发现，虽然H-G相位函数可以更好地体现出米氏散射前向峰值的主要特点，但是它却不能准确地模拟后向散射。1975年，Kattawar采用双H-G相位函数（TTHG）[27]来改善H-G相位函数的缺陷，如式（2-17）所示。

当g1的取值趋近1，g2的取值为负数时，双H-G相位函数的变化较为尖锐。从式（2-16）中可以发现，虽然双H-G相位函数能够更好地表示出后向散射的峰值，但是由于它是由f、g1、g2这3个不同的参数确定的，所以它的计算过程相对来说比较复杂。

Riewe与Green[28]将H-G相位函数和勒让德（Legendre）多项式结合起来，用来近似表示散射相位函数，如式（2-17）所示。

当设定参数g=0.72、f=0.5时，式（2-17）比较接近米氏散射相位函数。

Cornette和Shanks[29]定义了只含有一个参数的H-G相位函数，如式（2-18）所示。

此相位函数可以准确地表示海水的散射函数，但是该函数表达式相对比较复杂，而且利用该函数不能直接求解出散射角θ的解析式，因此很难应用在蒙特卡洛模拟中。

综上所述，由于本节只考虑前向散射，而且要选择能够应用在蒙特卡洛模拟中的函数，所以在仿真中选择的是式（2-15）所示的H-G相位函数。

2.2.3 海水信道中光传输散射特性的蒙特卡洛仿真

1. H-G相位函数的散射角

由文献[25-26]可知信道散射角θ的解析表达式为

其中，g表示非对称因子，r表示在[0,1]上均匀分布的随机数。

2. H-G相位函数的非对称因子

由于海水中含有的粒子直径大小不一致，所以一般用平均散射相位函数来表示特定的海水散射特性。Petzold通过对海水中平均粒子的散射相位函数中的散射相位余弦cosθ进行测量，测量的结果表明它的平均值为0.924。则非对称因子g=cosθ=0.924。

3. 海水中光信号蒙特卡洛模拟过程

蒙特卡洛模拟法被广泛应用于光束在海水介质的仿真研究中。在蒙特卡洛模拟法中，通常将光看成一系列的光子包，然后对这些光子包的状态进行追踪，最后统计出光子包的分布情况，以此来表示光的分布情况。

光子在海水中由于随机运动发生散射，利用蒙特卡洛模拟法来模拟光子的运动轨迹，通常可归纳为4步：第一步是光子的随机产生，第二步是光子的随机迁移，第三步是光子的吸收或散射，第四步是判定光子是否终止。判断光子终止的条件如下。

① 由光子权重决定。当光子的权重小于阈值时，可认为此光子已经消亡，从而不对其进行追踪。一般在模拟中，阈值主要由实验设备中探测器的灵敏度决定。

② 由介质几何限制决定。在随机碰撞过程中，如果光子被散射到介质外就不再对其进行追踪。

图2-12表示的是本节的蒙特卡洛模拟的仿真流程。在仿真中对光子的运动方向和传输距离进行判断。当光子运动方向Uz<0时，就终止对这个光子的追踪，而开始对另一个光子进行追踪；当光子运动方向Uz≥0时，再判断此光子是否已经到达接收面，即对此次光子的传输距离z进行判断。若光子未到达接收面，则让其继续碰撞后再判断，直到到达接收面。若光子到达接收面，则对光子的位置和权重进行统计。

光子在传输中的状态是由光子的坐标R(X,Y,Z)、运动方向Ω0、运动路径长度L以及权值W等决定的，用式（2-20）表示。

Sn表示从发射端发射出的光子经过n次随机碰撞的状态。

其中，Rn是指光子在第n次随机碰撞后的位置坐标；Ωn表示的是光子在第n次碰撞后的运动方向；Ln表示的是光子在第n次碰撞后已经传输的路径总长度；Wn指的是光子在第n次碰撞后的权值，通常用它来表示由于散射而造成的能量损失。

光子在第n次散射发生前的位置坐标可由(Xn-1,Yn-1,Zn-1)来表示，图2-13所示为光子在海水介质中的传输模型。

对蒙特卡洛模拟法[30]模拟的过程进行分析后，总结的步骤如下。

（1）光子的产生

在蒙特卡洛模拟过程中，假定光束是一个垂直向下的准直极细的入射光束，光束在发散角范围内是呈高斯分布的。此次模拟仿真中光子的位置是用三维坐标(X,Y,Z)来表示，假设光子的起始坐标是(0,0,0)，初始方向Ω0=(0,0,1)，本书中规定向下为正方向。

（2）确定光子坐标

光子在传输过程中发生碰撞时，它会用新的方向角(θn,φn,dn)进行下一次的传输。其中，θn指散射角，表示入射光与散射方向上的光在散射面上的夹角；φn指方向角，表示光子的散射方向在水平面上投影的旋转角度；dn指光子的自由路径，表示在下一次散射发生之前光子的自由路径。

光子从第n次到第n+1次散射之间的路径dn的分布用概率密度函数p(dn)来表示，如式（2-22）所示。

其中，c表示介质的衰减系数，即海水的总衰减系数。

光子的自由路径dn如式（2-23）所示。

通常情况下，蒙特卡洛模拟法模拟的自由路径用dn表示，通常选择如式（2-24）所示的表达式。

其中，r1表示的是[0,1]上均匀分布的随机数。

假设已知光子当前位置为(Xn,Yn,Zn)，它的传输方向为（μxn,μyn,μzn），预测达到下一个散射点的位置坐标为(Xn+1,Yn+1,Zn+1)，其如式（2-25）所示。

（3）确定散射后光子的运动方向

当光子发生散射后，光子的散射方向(φn,θn)可由式（2-26）来表示。

其中，光子发生第n次碰撞后的运动方向相对于碰撞前的运动方向的散射角用θn来表示，方位角用φn来表示，前者通过对H-G相位函数抽样来获得，后者则是(0,2π)上均匀分布函数的一个抽样值，r2、r3各是[0,1]上均匀分布的随机数。

对式（2-26）进行坐标变换，得到光子发生碰撞后新的方向矢量，用式（2-27）来表示。

若光子的运动方向非常靠近传输方向z轴，那么光子新的运动方向可用式（2-28）表示。

（4）权值改变

假设所有光子碰撞前的初值都为1，若光子在海水运动的过程中发生碰撞，那么其权值就会减少。因此，光子发生第n次碰撞之后的权值W可以由式（2-29）表示。

其中，Wn是光子在第n次碰撞发生前的权值，ω0是介质的单次散射率，其表示的是海水中所含杂质的散射系数与总衰减系数的比值。

（5）判决过程

在光子每一次发生散射之前，都必须比较光子当前权值和阈值的关系，若该光子当前权值比所确定的阈值小，则必须终止对该光子运动的模拟过程。若该光子当前的权值不小于所设定的阈值，就返回“确定光子坐标”的过程，此时继续进行蒙特卡洛模拟过程。如果经过判断，某个光子终止以后，则继续对下一个光子的模拟过程进行判断，直至全部光子都完成整个蒙特卡洛模拟过程。