CN119652425B - 基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法及系统 - Google Patents

Abstract

本发明涉及一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法及系统，该方法包括步骤：S10，加载光子并进行初始化设置；S20，计算光子在信道中随机游走过程中的位置变化；S30，根据光子的位置变化判断光子是否达到探测面，如果是则返回步骤S10，否则进入步骤S40；S40，计算光子的吸收率、局部正交坐标转换；S50，通过生存轮盘选择决定光子是否存活，如果是则返回步骤S20，如果否则进入步骤S60；S60，判断所有光子是否计算完毕，如果是则进入步骤S70，如果否则返回步骤S10；S70，计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存。本发明可以更准确模拟涡旋光跨海面传输过程，且得到累加电场信息。

Description

基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法及系统

技术领域

本发明涉及光通信技术领域，尤其涉及一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法及系统。

背景技术

光束在海洋环境中的传输特性对海洋光通信、海上激光雷达等应用具有重要意义。蒙特卡罗方法(Monte Carlo Method)是一种基于随机采样的数值计算方法，广泛应用于模拟各种复杂的物理过程，包括光的传播和散射。在光学领域，尤其是光在复杂介质中的传播(如空气、海水等)，蒙特卡洛法通过跟踪大量光子轨迹来模拟光子在介质中的吸收、散射等行为。通过这种方式，蒙特卡罗法能够在统计上准确地描述光的传输过程，特别适用于处理多次散射、随机介质中的光传输等问题。

然而传统的蒙特卡罗方法通常只关心光子的宏观状态，如位置、方向、能量等，例如公开号为CN 114884599A的发明《蓝绿涡旋光气海跨介质准直入射单程传输模型的构建方法》，因此，在模拟过程中，光子经过多次散射后的相干性和干涉效应常常被忽略，导致涡旋光的传输过程模拟不够精确。

发明内容

本发明的目的在于提供一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法及系统，以改善现有技术中所存在的涡旋光跨海面传输模拟不够精确的问题。

为了实现上述目的，本发明提供了以下技术方案：

第一方面，本发明提供了一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，包括以下步骤：

S10，加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦、初始的复电场平行分量和复电场垂直分量，初始的复电场平行分量和复电场垂直分量分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数，为方位角；

S20，计算光子在信道中随机游走过程中的位置变化；

S30，根据光子的位置变化判断光子是否达到探测面，如果是则返回步骤S10，模拟下一个光子的散射过程；如果否则进入步骤S40；

S40，计算光子的吸收率、局部正交坐标转换，然后进入步骤S50；

S50，通过生存轮盘选择决定光子是否存活，如果是，则返回步骤S20，继续模拟光子的后续行为，如果否则光子不再存活，进入步骤S60；

S60，判断所有光子是否计算完毕，如果是则进入步骤S70，如果否则返回步骤S10模拟下一个光子的散射过程；

S70，计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存。

上述方案主要是针对涡旋光在海洋水下信道中传输时受到的吸收、散射的影响，使用光子作为光信息的载体，并记录光子在每次Mie散射后的电场信息,以及位置和方向，从而更准确地计算光在水下的传播和散射过程，弥补了现有仿真方法的不足,并结合起伏海面对光束的影响，得到了涡旋光束在跨海面传输时对涡旋光的电场强度E_ins，为光束跨介质传输的提供了新的思路和方法。

第二方面，本发明提供了一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟系统，包括：

初始化模块，用于加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦、初始的复电场平行分量和复电场垂直分量，初始的复电场平行分量和复电场垂直分量分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数，为方位角；

位置变化探测模块，用于计算光子在信道中随机游走过程中的位置变化；

第一判断模块，用于根据光子的位置变化判断光子是否达到探测面，如果是则返回初始化模块，模拟下一个光子的散射过程；如果否则进入第一计算模块；

第一计算模块，用于计算光子的吸收率、局部正交坐标转换；

第二判断模块，用于通过生存轮盘选择决定光子是否存活，如果是，则返回位置变化探测模块，继续模拟光子的后续行为，如果否则光子不再存活，进入第三判断模块；

第三判断模块，用于判断所有光子是否计算完毕，如果是则进入第二计算模块，如果否则返回初始化模块模拟下一个光子的散射过程；

第二计算模块，用于计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存。

第三方面，本发明提供了一种计算机程序产品，包括计算机可读指令，其特征在于，所述计算机可读指令在被处理器执行时实现本发明基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。

第四方面，本发明提供了一种包括计算机可读指令的计算机可读存储介质，其特征在于，所述计算机可读指令在被处理器执行时实现本发明基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。

第五方面，本发明提供了一种电子设备，包括：存储器，存储程序指令；处理器，与所述存储器相连接，执行存储器中的程序指令，实现本发明基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。

与现有技术相比，本发明具有以下显著的技术优势：

尤其在结合海洋与空气界面形成一个适用的系统方面，通过详细记录光子的电场信息和位置，本发明增强了对光传播过程的计算精度。这种记录方式使得每次Mie散射后的电场信息得以保留，从而能够准确反映光束在传播过程中的变化，尤其是在光从水下到空气界面和从空气界面到水下的传输过程中，均能实现有效的模拟。

此外，本发明综合考虑了海面起伏对光束传播路径和强度的影响，特别是在实际海洋环境中，海面波动会显著影响光的传播。通过对这些因素的全面分析，本发明能够准确预测涡旋光在跨海面传输中的电场强度，使得系统能够适用于海洋和空气界面之间的双向传输场景。

为了更有效地处理复杂的散射和吸收问题，本发明采用了电场蒙特卡罗方法，能够反映多重散射效应对电场强度的影响。在电场的初始化和更新过程中，本发明所采用的公式发挥了重要作用。电场初始化公式用于设定光子进入散射介质前的电场状态，其特殊效果在于为后续的传播和散射过程提供了清晰的起点，使得光场的模拟更加准确。初始化时，电场分布受到入射光束类型(如平面波或涡旋光束)的影响，确保模型能够真实反映实际情况，这为有效模拟提供了基础，帮助研究者更好地理解光在介质中的传播特性。

其次，光子经过散射后的电场更新公式在每次散射后更新光子电场的变化，实时反映光子与散射介质相互作用后的电场强度变化。该公式考虑了散射角和散射振幅等因素，确保电场信息的动态更新。这一动态更新机制不仅提高了模型的精度，还能有效模拟多次散射对光传播的影响，更真实地反映光束的衰减特性和传播路径，从而实现了从海洋到空气界面和从空气界面到海洋界面的双向适用性。

总而言之，本发明通过结合海洋和空气界面，形成了一个适用于两种传输场景的方法及系统，具备了高效、准确和可靠的技术优势。

本发明还具有的其他优势请见实施例部分的相关描述。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为实施例中举例的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法的流程图。

图2a、图2b、图2c分别是无风、风速6m/s、风速8m/s时光束强度的变化光场分布图。

图3为实施例中举例的基于电场蒙特卡罗法的涡旋光跨海面传输模拟系统的组成框图。

图4为实施例中电子设备的组成框图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图1，本实施例中提供了一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，包括以下步骤：

S10，加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦。

使用EMC仿真方法模拟OAM(涡旋光)光束在水下散射信道传输时，其中最重要的一个步骤就是对光子初始坐标的设置。设初始发射的轴向坐标为z＝0，在极坐标系下，需要对光子初始发射面上发射坐标进行选取，ρ与OAM光束光强的径向分布特性有关，光强随ρ分布的大小决定了抽取ρ的概率。

OAM光束在发射面上的光强分布与径向基函数有关。对于贝塞尔高斯光束,其径向基函数可表示为：

其中，E₀为常数，ω₀为束腰半径，k_ρ为径向波矢，l为贝塞尔高斯光束的阶数，也即拓扑荷数,可取任意整数，ρ为初始径向坐标，其中E₀、ω₀、k_ρ、l可以直接设定，为了方便观察，这里将光束设定成一个圆环的光束分布样式，J_l(·)表示的是l阶的第一类Bessel函数，函数中的“·”表示参数。

根据此径向基函数,可以得到归一化的光强密度函数：

这个归一化的光强密度函数描述了OAM光束在发射面上的径向光强分布情况。在仿真中,可以通过拒绝算法采样光强密度函数来确定每个光子的初始径向坐标ρ。

EMC方法发射光子ρ的随机抽样服从归一化光强密度函数对应的概率密度分布；为了方便采样，这里设置入射OAM光束的截止半径为R_m，令ρ为[0,R_m]范围内的均匀分布。归一化光强密度函数I(ρ)在该ρ位置上发射的光子数量越多，发射的概率越大。

由此，采用拒绝算法对输入OAM光束的ρ₀采样的流程如下：

(1)产生一个随机数χ，满足χ～U(0,1)，计算得到一个径向坐标ρ₀＝χR_m。

(2)产生一个随机数满足如果满足则接受此ρ₀。

对于方位角的采样,由于OAM光束在极坐标系中具有径向对称性,因此可以直接采用0到2π之间的均匀分布。

此外,为了保证光子的初始电场分量,此处根据OAM光束的特性，对于初始的复电场平行分量和复电场垂直分量，分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数。

使用笛卡尔坐标系，光子初始坐标可以表示为：

初始方向余弦可以表示为：

其中θ₀和是光子发射方向上的仰角和方位角,θ₀分别由下式得到：

其中r₃、r₄表示0～1之间的随机数，在0～1内均匀分布，θ₁为光束的发散角。

S20，光子在信道中随机游走过程中，计算光子的位置变化。

定义光子散射前位置坐标为(x,y,z)，散射后位置坐标为(x′,y′,z′)，光子散射前的方向余弦为(u,v,w)，假设光子发生1次散射的传播距离为l_sct，根据Beer's Law，光子散射步长表示为l_sct＝-ln(ζ)/Q_sca，其中Q_sca为散射系数，ζ为一个符合0～1之间均匀分布的随机数。由此，光子由散射前位置坐标(x,y,z)传输了距离l_scc后,位置坐标(x′,y′,z′)可表示为：

x′＝x+ul_sct

y′＝y+vl_sct

z′＝z+wl_sct；

方向余弦在散射后可由下式更新：

θ、分别为散射角和方位角。

光子在不断散射过程中，散射角和方位角会发生变化，每散射一次更新一次散射角和方位角。

在不同的散射角度范围内，选择最合适的相位函数。常用的HG相位函数的参量较少，所以简单易用。HG相位函数的表达式为：

通过相位函数分布规律可以设置散射角θ：

θ为散射角；g为不对称因子，其大小影响前向散射与后向散射的比重。ζ为一个符合0～1之间均匀分布的随机数。给定ζ和g即可计算得到散射角θ。

方位角在0-2π之间均匀分布，可以表示为在散射循环中使用。

S30，根据光子的位置变化判断光子是否达到探测面，如果是则返回步骤S10，模拟下一个光子的散射过程；如果否则进入步骤S40。

例如，假设一个光子从原点(0.1,0.5,0)出发，经过几次传播后，其位置更新为(1,1,6)。在这种情况下，因为z坐标(6)大于探测面的z轴坐标z_d(z_d＝5)，所以判断为光子已到达探测面。

如果散射角不为零，则判断为已经发生散射。

S40，计算光子的吸收率、局部坐标等参数，然后进入步骤S50。

当光子与球形介质粒子相互作用时，一部分能量会被介质粒子吸收和散射。粒子散射系数、吸收系数、消光系数分别为Q_sca、Q_abs、Q_ext，其表达式如下：

Q_ass＝Q_ext-Q_sca。

a_N与b_N分别是Mie系数，符号表示取实部，x_r表示散射球形粒子的相对粒径大小。N表示不同阶数的项，用于表示不同阶的球谐函数的贡献，通过累加这些不同阶数的项，可以更准确地描述散射、吸收和消光等参数。

为了表示传播过程中光束的吸收,可以通过调整每一次散射后光子的权重W_a来实现，权重变化可表示为W_a＝W(1-Q_sca/Q_ext)，Q_sca/Q_ext表示单次散射反照率，表示散射在消光中所占比例，W为光子初始权重,吸收率即为

将散射前光子的局部正交坐标系记为(m,n,s)，与介质粒子相互作用后的坐标系记为(m',n',s')，m和n分别为散射前电场的平行分量、垂直分量，且为相对于上一次散射的散射平面方向上的单位矢量。电场分量是复电场分量的实部，电场分量描述了电场在位置上的实际振幅，物理测量仪器可以直接测量的量，复电场分量提供了完整的振幅和相位信息。s为光子在散射前的运动方向(可称为入射方向)，s'为与介质粒子散射后的光子运动方向(可称为散射方向)。散射方向s'与入射方向s决定了散射平面。单位矢量m'和n'分别为散射后散射光平行于于散射面的平行分量、垂直于散射面的垂直分量，并且满足m'×n'＝s'。E₁、E₁′、E₂、E₂′分别为散射前电场平行方向上的复电场分量、散射后电场平行方向上的复电场分量、散射前电场垂直方向上的复电场分量、散射后电场垂直方向上的复电场分量。散射前后的光子局部坐标变换规则可以通过向量运算抽象出来，并可用矩阵形式表示。每次光子被散射时，局部坐标系的变化满足以下关系：

其中θ表示散射角,表示方位角。

同理,光子经过散射后复电场的更新可表示为：

其中，S₁和S₂是关于散射角θ的散射振幅函数，其与散射系数、散射角函数有关。根据Mie散射理论，其分别可以表示为：

其中a_N与b_N分别是Mie系数，是与散射介质粒子实际尺寸、入射光波长以及散射介质粒子复折射率有关的函数，但与散射角θ无关。π_N与τ_N是散射角函数，且只与散射角θ有关，它们分别可以表示为：

其中，x_r＝πd/λ，表示散射球形粒子的相对粒径大小，其中d为粒子的实际直径尺寸，单位为米，λ为波长。散射球形粒子的复折射率为：n_sct＝n_real+i*n_imag，其中实部n_real与介质粒子散射作用相关，虚部n_imag与介质粒子吸收作用相关，n_imag不为0时，表示光子会受到介质粒子的吸收作用而造成能量衰减。ψ_N(x)与ξ_N(x)分别为半整数阶Bessel函数和第二类汉克尔(Hankel)函数。P_N(cosθ)与P_N′(cosθ)分别为关于cosθ的勒让德函数和一阶缔合勒让德函数。表示沿方向传播的散射光的强度归一化因子，表示为：

其中，符号表示取实部，*表示取共轭，为方位角，同时满足|E₁|²+|E₂|²＝1。

S50，通过生存轮盘选择决定光子是否存活，如果是，则返回步骤S20，继续模拟光子的后续行为，如果否则光子不再存活，进入步骤S60。

当光子在传输过程中未到达接收面时，光子权重值会随着光子游走不断减小，当其小于设置的湮灭值时，即代表光子在运动过程中湮灭，不予在接收面记录，因为即使当光子到达接收面时，对光束强度的贡献也几乎可以忽略不计，将湮灭的光子忽略掉可以加快程序的执行效率，即加快模拟效率。而当光子到达接收面，且权重大于湮灭值时，光子终止运动并记录累加的电场信息。

S60，判断所有光子是否计算完毕，如果是则进入步骤S70，如果否则返回步骤S10模拟下一个光子的散射过程。

S70，计算透射系数、海面光束影响因子并作用到累加的电场上，得到电场强度，并记录保存。

光子经过n次介质粒子的散射作用后，其最终局部正交坐标系为(m⁽ⁿ⁾,n⁽ⁿ⁾,s⁽ⁿ⁾)，为了方便统计,这里限定所有光子的最后运动方向s⁽ⁿ⁾均为(u_d,v_d,w_d)＝(0,0,1)即向z轴正方向传播达到探测面。此时m⁽ⁿ⁾和n⁽ⁿ⁾方向上的平行电场分量和垂直电场分量分别为 对应的x、y、z三个方向上的电场分量可由根据矢量变换得到,同时考虑到光子与散射介质的相互作用以及在散射介质传播过程中而产生的相位延迟,经过光子叠加以后x、y、z三个方向的电场分量可表示为：

其中G为x、y、z三个方向，k为波数，l⁽ⁿ⁾为光子散射n次后移动的总路径长度。

本实施例中探测面设置在z＝z_d的位置上，而实际光子经过n次散射后到达最终散射位置z轴方向坐标为z⁽ⁿ⁾，一般来说，光子最终很大概率不会直接落在探测面(z轴方向坐标为z_d)上被接收，即z⁽ⁿ⁾＞z_d或者z⁽ⁿ⁾≤z_d。这种情况下，需要考虑光子最终散射位置到探测面的相延，可表示为另外，光子从最终散射位置到探测位置存在未统计到的散射过程造成光子能量的衰减，因此可设置一个权重衰减因子为来等效该衰减作用，其中散射截面l_sct为散射随机步长。

综合考虑上述所有因素，探测面接收到的光子在x,y,z方向上的电场分量可表示为：

w_d表示光子到达接收平面处的方向余弦。

由菲涅耳公式可知，光束从水下向海面上传输的总透射率T_up为：

T_up＝T_wa1T_wa2；

其中T_wa1表示气海界面本身的透过率，满足α表示光束从海洋-空气界面传输的入射角，β表示折射角。根据折射定理可知，折射角的表达式n₁、n₂分别表示海水折射率、空气折射率。T_wa2表示海面上泡沫以及其他因素造成的透过率，当风速p≤9m/s时，满足T_wa2＝1-1.2×10^-5×p^3.3，当风速为p＞9m/s时，满足T_wa2＝1-1.2×10^-5×p^3.3×(0.225p-0.99)。

贝塞尔高斯光束穿过海面后的x,y,z方向上的电场分量为：

其中表示从海面下方的光束经过海水多次散射传输之后的贝塞尔高斯光束各分量电场,是贝塞尔高斯光束的波数，λ是光束波长，z_h表示海面高度起伏函数，n₁为海水折射率，T为总透射率。总的电场强度E_ins的大小可以使用以下公式表示：

其中E_x、E_y、E_z分别是电场在x、y、z方向上的电场分量。

前面已经探讨了光在水下传播并穿过海洋界面的上行传输过程。现在，当光线从空气中传输，透过海洋界面下行传输时的传输过程为：

S100，加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦。

使用EMC仿真方法模拟OAM(涡旋光)光束在空气散射信道传输时，其中最重要的一个步骤就是对光子初始坐标的设置。设初始发射的轴向坐标为z＝0，在极坐标系下，需要对光子初始发射面上发射坐标进行选取，ρ与OAM光束光强的径向分布特性有关，光强随ρ分布的大小决定了抽取ρ的概率。

OAM光束在发射面上的光强分布与径向基函数有关。对于贝塞尔高斯光束,其径向基函数可表示为：

其中，E₀为常数，ω₀为束腰半径，k_ρ为径向波矢，l为贝塞尔高斯光束的阶数，也即拓扑荷数,可取任意整数，ρ为初始径向坐标，其中E₀、ω₀、k_ρ、l可以直接设定，为了方便观察，这里将光束设定成一个圆环的光束分布样式，J_l(·)表示的是l阶的第一类Bessel函数，函数中的“·”表示参数。

根据此径向基函数,可以得到归一化的光强密度函数：

这个归一化的光强密度函数描述了OAM光束在发射面上的径向光强分布情况。在仿真中,可以通过拒绝算法采样光强密度函数来确定每个光子的初始径向坐标ρ。

EMC方法发射光子ρ的随机抽样服从归一化光强密度函数对应的概率密度分布；为了方便采样，这里设置入射OAM光束的截止半径为R_m，令ρ为[0,R_m]范围内的均匀分布。归一化光强密度函数I(ρ)在该ρ位置上发射的光子数量越多，发射的概率越大。

由此，采用拒绝算法对输入OAM光束的ρ₀采样的流程如下：

(1)产生一个随机数χ，满足χ～U(0,1)，计算得到一个径向坐标ρ₀＝χR_m。

(2)产生一个随机数满足如果满足则接受此ρ₀。

对于方位角的采样,由于OAM光束在极坐标系中具有径向对称性,因此可以直接采用0到2π之间的均匀分布。

此外,为了保证光子的初始电场分量,此处根据OAM光束的特性，对于初始的复电场平行分量和复电场垂直分量，分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数。

使用笛卡尔坐标系，光子初始坐标可以表示为：

初始方向余弦可以表示为：

其中θ₀和是光子发射方向上的仰角和方位角,θ₀分别由下式得到：

其中r₃、r₄表示0～1之间的随机数，在0～1内均匀分布，θ₁为光束的发散角。

S200，光子在信道中随机游走过程中，计算光子的位置变化。

定义光子散射前位置坐标为(x,y,z)，散射后位置坐标为(x′,y′,z′)，光子散射前的方向余弦为(u,v,w)，假设光子发生1次散射的传播距离为l_sct，根据Beer's Law，光子散射步长表示为l_sct＝-ln(ζ)/Q_sca，其中Q_sca为散射系数，ζ为一个符合0～1之间均匀分布的随机数。由此，光子由散射前位置坐标(x,y,z)传输了距离l_scc后,位置坐标(x′,y′,z′)可表示为：

x′＝x+ul_sct

y′＝y+vl_sct

z′＝z+wl_sct；

方向余弦在散射后可由下式更新：

θ、分别为散射角和方位角。

光子在不断散射过程中，散射角和方位角会发生变化，每散射一次更新一次散射角和方位角。

在不同的散射角度范围内，选择最合适的相位函数。常用的HG相位函数的参量较少，所以简单易用。HG相位函数的表达式为：

通过相位函数分布规律可以设置散射角θ：

θ为散射角；g为不对称因子，其大小影响前向散射与后向散射的比重。ζ为一个符合0～1之间均匀分布的随机数。给定ζ和g即可计算得到散射角θ。

方位角在0-2π之间均匀分布，可以表示为在散射循环中使用。

S300，根据光子的位置变化判断光子是否达到探测面，如果是则返回步骤S100，模拟下一个光子的散射过程；如果否则进入步骤S400。

例如，假设一个光子从原点(0.1,0.5,0)出发，经过几次传播后，其位置更新为(1,1,6)。在这种情况下，因为z坐标(6)大于探测面的z轴坐标z_d(z_d＝5)，所以判断为光子已到达探测面。

如果散射角不为零，则判断为已经发生散射。

S400，计算光子的吸收率、局部坐标等参数，然后进入步骤S500。

当光子与球形介质粒子相互作用时，一部分能量会被介质粒子吸收和散射。粒子散射系数、吸收系数、消光系数分别为Q_sca、Q_abs、Q_ext，其表达式如下：

Q_ass＝Q_ext-Q_sca。

a_N与b_N分别是Mie系数，符号表示取实部，x_r表示散射球形粒子的相对粒径大小。N表示不同阶数的项，用于表示不同阶的球谐函数的贡献，通过累加这些不同阶数的项，可以更准确地描述散射、吸收和消光等参数。

为了表示传播过程中光束的吸收,可以通过调整每一次散射后光子的权重W_a来实现，权重变化可表示为W_a＝W(1-Q_sca/Q_ext)，Q_sca/Q_ext表示单次散射反照率，表示散射在消光中所占比例，W为光子初始权重,吸收率即为

将散射前光子的局部正交坐标系记为(m,n,s)，与介质粒子相互作用后的坐标系记为(m',n',s')，m和n分别为散射前电场的平行分量、垂直分量，且为相对于上一次散射的散射平面方向上的单位矢量。电场分量是复电场分量的实部，电场分量描述了电场在位置上的实际振幅，物理测量仪器可以直接测量的量，复电场分量提供了完整的振幅和相位信息。s为光子在散射前的运动方向(可称为入射方向)，s'为与介质粒子散射后的光子运动方向(可称为散射方向)。散射方向s'与入射方向s决定了散射平面。单位矢量m'和n'分别为散射后散射光平行于于散射面的平行分量、垂直于散射面的垂直分量，并且满足m'×n'＝s'。E₁、E₁′、E₂、E₂′分别为散射前电场平行方向上的复电场分量、散射后电场平行方向上的复电场分量、散射前电场垂直方向上的复电场分量、散射后电场垂直方向上的复电场分量。散射前后的光子局部坐标变换规则可以通过向量运算抽象出来，并可用矩阵形式表示。每次光子被散射时，局部坐标系的变化满足以下关系：

其中θ表示散射角,表示方位角。

同理,光子经过散射后复电场的更新可表示为：

其中，S₁和S₂是关于散射角θ的散射振幅函数，其与散射系数、散射角函数有关。根据Mie散射理论，其分别可以表示为：

其中a_N与b_N分别是Mie系数，是与散射介质粒子实际尺寸、入射光波长以及散射介质粒子复折射率有关的函数，但与散射角θ无关。π_N与τ_N是散射角函数，且只与散射角θ有关，分别可以表示为：

其中，x_r＝πd/λ，表示散射球形粒子的相对粒径大小，其中d为粒子的实际直径尺寸，单位为米，λ为波长。散射球形粒子的复折射率为：n_sct＝n_real+i*n_imag，其中实部n_real与介质粒子散射作用相关，虚部n_imag与介质粒子吸收作用相关，n_imag不为0时，表示光子会受到介质粒子的吸收作用而造成能量衰减。ψ_N(x)与ξ_N(x)分别为半整数阶Bessel函数和第二类汉克尔(Hankel)函数。P_N(cosθ)与P_N′(cosθ)分别为关于cosθ的勒让德函数和一阶缔合勒让德函数。表示沿方向传播的散射光的强度归一化因子，表示为：

其中，符号表示取实部，*表示取共轭，为方位角，同时满足|E₁|²+|E₂|²＝1。

S500，通过生存轮盘选择决定光子是否存活，如果是，则返回步骤S200，继续模拟光子的后续行为，如果否则光子不再存活，进入步骤S600。

当光子在传输过程中未到达接收面时，光子权重值会随着光子游走不断减小，当其小于设置的湮灭值时，即代表光子在运动过程中湮灭，不予在接收面记录，因为即使当光子到达接收面时，对光束强度的贡献也几乎可以忽略不计，将湮灭的光子忽略掉可以加快程序的执行效率，即加快模拟效率。而当光子到达接收面，且权重大于湮灭值时，光子终止运动并记录累加的电场信息。

S600，判断所有光子是否计算完毕，如果是则进入步骤S700，如果否则返回步骤S100模拟下一个光子的散射过程。

S700，计算透射系数、海面光束影响因子并作用到累加的电场上，得到电场强度，并记录保存。

光子经过n次介质粒子的散射作用后，其最终局部正交坐标系为(m⁽ⁿ⁾,n⁽ⁿ⁾,s⁽ⁿ⁾)，为了方便统计,这里限定所有光子的最后运动方向s⁽ⁿ⁾均为(u_d,v_d,w_d)＝(0,0,1)即向z轴正方向传播达到探测面。此时m⁽ⁿ⁾和n⁽ⁿ⁾方向上的平行电场分量和垂直电场分量分别为 对应的x、y、z三个方向上的电场分量可由根据矢量变换得到,同时考虑到光子与散射介质的相互作用以及在散射介质传播过程中而产生的相位延迟,经过光子叠加以后x、y、z三个方向的电场分量可表示为：

其中G为x、y、z三个方向，k为波数，l⁽ⁿ⁾为光子散射n次后移动的总路径长度。

本实施例中探测面设置在z＝z_d的位置上，而实际光子经过n次散射后到达最终散射位置z轴方向坐标为z⁽ⁿ⁾，一般来说，光子最终很大概率不会直接落在探测面(z轴方向坐标为z_d)上被接收，即z⁽ⁿ⁾＞z_d或者z⁽ⁿ⁾≤z_d。这种情况下，需要考虑光子最终散射位置到探测面的相延，可表示为另外，光子从最终散射位置到探测位置存在未统计到的散射过程造成光子能量的衰减，因此可设置一个权重衰减因子为来等效该衰减作用，其中散射截面l_sct为散射随机步长。

综合考虑上述所有因素，探测面接收到的光子在x,y,z方向上的电场分量可表示为：

w_d表示光子到达接收平面处的方向余弦。

由菲涅耳公式可知，通过海面的光束的总透射率由下式给出：

T_down＝T_wa1T_wa2；

其中T_wa1表示气海界面本身的透过率，满足α表示光束从空气-海洋界面传输的入射角。根据折射定理可知，折射角的表达式n₁、n₂分别表示海水折射率、空气折射率。T_wa2表示海面物质(例如海面上的泡沫以及其他因素)造成的透过率，当风速p≤9m/s时，满足T_wa2＝1-1.2×10^-5×p^3.3，当风速为p＞9m/s时，满足T_wa2＝1-1.2×10^-5×p^3.3×(0.225p-0.99)。T_down表示为光束从空气向下传输的总透射率；

贝塞尔高斯光束穿过海面后的x,y,z方向上的电场分量为：

其中表示从空气上方的光束经过空气多次散射传输之后的贝塞尔高斯光束各分量电场,是贝塞尔高斯光束的波数，λ是光束波长，z_h表示海面高度起伏函数，n₁为海水折射率，T为总透射率。总的电场强度E_ins的大小可以使用以下公式表示：

其中E_x、E_y、E_z分别是电场在x、y、z方向上的电场分量。

如图2a、图2b、图2c所示，图2a、图2b、图2c分别是无风、风速6m/s、风速8m/s时光束强度的变化光场分布，随着海面风速的增加，贝塞尔高斯光束的强度分布会逐渐发生畸变。当海面风速达到8m/s时，光强分布失真严重。这是由于海面风速越大，海面的透射率越低，透射光线也受到海面自身遮挡的影响。随着风速的增大，这种遮挡作用变得更加显著。从透射光强分布图来看，海面风速越高，激光透射光强的分布越不均匀。特别是在海面相对高度差较大的情况下，光线的透射受到了明显的遮挡，导致透射出的光线极少。

在一个试验例中，发射光束使用实验室常见且容易产生的贝塞尔高斯光束，为了平衡代码的执行效率和仿真精度，使用了下面的参数：

在光束从水下传输到海面或从空气向水下传输的过程中，会经历多次散射，现有的类似技术未能有效体现这一现象，本发明采用电场蒙特卡罗方法，将光束划分为多个光子，模拟这些光子在传播过程中与海水中散射介质的多次碰撞，从而真实反映光束从水下到海面传输时，由于多次散射导致的光束强度衰减。

请参阅图3，基于相同的发明构思，本实施例中同时提供了一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟系统，包括：

初始化模块，用于加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦、初始的复电场平行分量和复电场垂直分量，初始的复电场平行分量和复电场垂直分量分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数，为方位角；

位置变化探测模块，用于计算光子在信道中随机游走过程中的位置变化；

第一判断模块，用于根据光子的位置变化判断光子是否达到探测面，如果是则返回初始化模块，模拟下一个光子的散射过程；如果否则进入第一计算模块；

第一计算模块，用于计算光子的吸收率、局部正交坐标转换；

第二判断模块，用于通过生存轮盘选择决定光子是否存活，如果是，则返回位置变化探测模块，继续模拟光子的后续行为，如果否则光子不再存活，进入第三判断模块；

第三判断模块，用于判断所有光子是否计算完毕，如果是则进入第二计算模块，如果否则返回初始化模块模拟下一个光子的散射过程；

第二计算模块，用于计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存。

上述系统中各个组成模块的具体执行过程可以参见前述方法的相关描述，此处不再赘述。

如图4所示，本实施例同时提供了一种电子设备，该电子设备可以包括处理器41和存储器42，其中存储器42耦合至处理器41。值得注意的是，该图是示例性的，还可以使用其他类型的结构来补充或替代该结构，实现数据提取、报告生成、通信或其他功能。

如图4所示，该电子设备还可以包括：输入单元43、显示单元44和电源45。值得注意的是，该电子设备也并不是必须要包括图4中显示的所有部件。此外，电子设备还可以包括图4中没有示出的部件，可以参考现有技术。

处理器41有时也称控制器或操作控件，可以包括微处理器或其他处理器装置和/或逻辑装置，该处理器41接收输入并控制电子设备的各个部件的操作。

其中，存储器42例如可以是缓存器、闪存、硬驱、可移动介质、易失性存储器、非易失性存储器或其他合适装置中的一种或多种，可存储上述处理器41的配置信息、处理器41执行的指令等信息。处理器41可以执行存储器42存储的程序，以实现信息存储或处理等。在一个实施例中，存储器42中还包括缓冲存储器，即缓冲器，以存储中间信息。

本发明实施例还提供一种计算机程序产品，包括计算机可读指令，当在电子设备中执行所述计算机可读指令时，所述程序产品使得电子设备执行本发明方法所包含的操作步骤。

本发明实施例还提供一种存储有计算机可读指令的存储介质，所述计算机可读指令使得电子设备执行本发明方法所包含的操作步骤。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述实施例，仅为本发明具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改、替换和改进等等，这些修改、替换和改进都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (8)

1.一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，其特征在于，包括以下步骤：

S10，加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦、初始的复电场平行分量和复电场垂直分量，初始的复电场平行分量和复电场垂直分量分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数，为方位角；

S20，计算光子在信道中随机游走过程中的位置变化；

S30，根据光子的位置变化判断光子是否达到探测面，如果是则返回步骤S10，模拟下一个光子的散射过程；如果否则进入步骤S40；

S40，计算光子的吸收率、局部正交坐标转换，然后进入步骤S50；

S50，通过生存轮盘选择决定光子是否存活，如果是，则返回步骤S20，继续模拟光子的后续行为，如果否则光子不再存活，进入步骤S60；

S60，判断所有光子是否计算完毕，如果是则进入步骤S70，如果否则返回步骤S10模拟下一个光子的散射过程；

S70，计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存；

S10中，先采样获得初始径向坐标ρ₀，然后使用笛卡尔坐标系，光子初始坐标表示为：

初始方向余弦表示为：

其中θ₀和是光子发射方向上的仰角和方位角；θ₀分别由下式得到：

其中r₃、r₄表示0～1之间的随机数，在0～1内均匀分布，θ₁为光束的发散角；

采用拒绝算法对输入OAM光束的初始径向坐标采样的流程如下：

(1)产生一个随机数χ，满足χ～U(0,1)，计算得到一个径向坐标ρ₀＝χR_m；R_m为入射OAM光束的截止半径；

(2)产生一个随机数满足如果满足则接受此ρ₀，其中，E₀为常数，ω₀为束腰半径，k_ρ为径向波矢，J_l(·)表示的是l阶的第一类Bessel函数，函数中的“·”表示参数，E₀、ω₀、k_ρ、l可以直接设定。

2.根据权利要求1所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，其特征在于，S20中，定义光子散射前位置坐标为(x,y,z)，散射后位置坐标为(x′,y′,z′)，光子散射前的方向余弦为(u,v,w)，由此有：

x′＝x+ul_sct

y′＝y+vl_sct

z′＝z+wl_sct；

方向余弦在散射后由下式更新：

l_sct为光子发生1次散射的传播距离，l_sct＝-ln(ζ)/Q_sca，Q_sca为散射系数，ζ为一个符合0～1之间均匀分布的随机数；θ、分别为散射角和方位角。

3.根据权利要求2所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，其特征在于，S40中，将散射前光子的局部正交坐标系记为(m,n,s)，与介质粒子相互作用后的坐标系记为(m',n',s')，m'×n'＝s'；每次光子被散射时，局部坐标系的变化满足以下关系：

其中θ表示散射角,表示方位角；

光子经过散射后复电场的更新可表示为：

其中，E₁、E₁′、E₂、E₂′分别为散射前电场平行方向上的复电场分量、散射后电场平行方向上的复电场分量、散射前电场垂直方向上的复电场分量、散射后电场垂直方向上的复电场分量；S₁和S₂是关于散射角θ的散射振幅函数，

其中a_N与b_N分别是Mie系数，π_N与τ_N是散射角函数，

其中，x_r＝πd/λ，表示散射球形粒子的相对粒径大小，其中d为粒子的实际直径尺寸，单位为米，λ为波长，散射球形粒子的复折射率为：n_sct＝n_real+i*n_imag，其中n_real表示实部，n_imag表示虚部；ψ_N(x)与ξ_N(x)分别为半整数阶Bessel函数和Hankel函数，P_N(cosθ)与P_N′(cosθ)分别为关于cosθ的勒让德函数和一阶缔合勒让德函数，表示沿方向传播的散射光的强度归一化因子，表示为：

其中，符号表示取实部，*表示取共轭，|E₁|²+|E₂|²＝1。

4.根据权利要求3所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法，其特征在于，S70中，探测面接收到的光子在x,y,z方向上的电场分量可表示为：

其中，分别为m⁽ⁿ⁾和n⁽ⁿ⁾方向上的平行电场分量和垂直电场分量，i为虚数单位，是贝塞尔高斯光束的波数，λ是光束波长，l⁽ⁿ⁾为光子散射n次后移动的总路径长度，表示光子最终散射位置到探测面的相延，表示光子能量权重衰减因子；z_d为z轴方向坐标，散射截面w_d表示光子到达接收平面处的方向余弦，z⁽ⁿ⁾为光子经过n次散射后到达最终散射位置z轴方向坐标；G为x、y、z三个方向，光子经过n次介质粒子的散射作用后最终局部正交坐标系为(m⁽ⁿ⁾,n⁽ⁿ⁾,s⁽ⁿ⁾)；

当光束从水下向海面上传输时，贝塞尔高斯光束穿过海面后的x,y,z方向上的电场分量为：

其中T_wa1表示气海界面本身的透过率，满足α表示光束从海洋-空气界面传输的入射角，β表示折射角，T_wa2表示海面物质造成的透过率；z_h表示海面高度起伏函数，n₁、n₂分别表示海水折射率、空气折射率，T_up表示光束从水下向海面上传输的总透射率；

当光束从空气向海面下传输时，贝塞尔高斯光束穿过海面后的x,y,z方向上的电场分量为：

T_down＝T_wa1T_wa2；

总的电场强度E_ins的大小使用以下公式表示：

其中E_x、E_y、E_z分别是电场在x、y、z方向上的电场分量。

5.一种基于电场蒙特卡罗法的涡旋光跨海面传输模拟系统，其特征在于，包括：

初始化模块，用于加载光子并进行初始化设置，包括设置光子的初始坐标、初始方向余弦、初始的复电场平行分量和复电场垂直分量，初始的复电场平行分量和复电场垂直分量分别设置为和E₂＝0，其中i表示虚数单位，l表示涡旋光束的拓扑荷数，为方位角；

位置变化探测模块，用于计算光子在信道中随机游走过程中的位置变化；

第一判断模块，用于根据光子的位置变化判断光子是否达到探测面，如果是则返回初始化模块，模拟下一个光子的散射过程；如果否则进入第一计算模块；

第一计算模块，用于计算光子的吸收率、局部正交坐标转换；

第二判断模块，用于通过生存轮盘选择决定光子是否存活，如果是，则返回位置变化探测模块，继续模拟光子的后续行为，如果否则光子不再存活，进入第三判断模块；

第三判断模块，用于判断所有光子是否计算完毕，如果是则进入第二计算模块，如果否则返回初始化模块模拟下一个光子的散射过程；

第二计算模块，用于计算海面光束影响因子作用下的累加电场，得到电场强度，并记录保存；

初始化模块在设置光子的初始坐标、初始方向余弦时执行以下操作：先采样获得初始径向坐标ρ₀，然后使用笛卡尔坐标系，光子初始坐标表示为：初始方向余弦表示为：其中θ₀和是光子发射方向上的仰角和方位角；θ₀分别由下式得到：其中r₃、r₄表示0～1之间的随机数，在0～1内均匀分布，θ₁为光束的发散角；

采用拒绝算法对输入OAM光束的初始径向坐标采样的流程如下：(1)产生一个随机数χ，满足χ～U(0,1)，计算得到一个径向坐标ρ₀＝χR_m；(2)产生一个随机数满足如果满足则接受此ρ₀，其中， E₀为常数，ω₀为束腰半径，k_ρ为径向波矢，J_l(·)表示的是l阶的第一类Bessel函数，函数中的“·”表示参数，E₀、ω₀、k_ρ、l可以直接设定。

6.一种计算机程序产品，包括计算机可读指令，其特征在于，所述计算机可读指令在被处理器执行时实现权利要求1-4任一项所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。

7.一种包括计算机可读指令的计算机可读存储介质，其特征在于，所述计算机可读指令在被处理器执行时实现权利要求1-4任一项所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。

8.一种电子设备，其特征在于，包括：存储器，存储程序指令；处理器，与所述存储器相连接，执行存储器中的程序指令，实现权利要求1-4任一项所述的基于电场蒙特卡罗法的涡旋光跨海面传输模拟方法中的步骤。