CN102594440A - 一种光子传输性能仿真方法 - Google Patents

Abstract

本发明公开了一种光子传输性能仿真方法，包括：步骤1：设定待仿真光子数目、各光子的初始状态和介质光学特性参数；步骤2：根据初始状态发射光子，对每一光子传输路径进行追踪，追踪过程中，对经过介质界面或接收平面的光子采用修正的蒙特卡罗方法实现追踪；步骤3：基于追踪结果进行统计计算，生成光斑分布图、多径时延概率图、最大多径时延、63％光子平均多径时延和最大时延、光子所剩数目并显示。本发明是针对光子运动路径上的介质界面及接收平面，提出了修正的蒙特卡罗方法，使得仿真结果更加准确，为进一步开发激光对水下接收平台通信系统提供了更精确的仿真数据。

Description

一种光子传输性能仿真方法

技术领域

本发明属于激光通信领域，特别是一种对激光对水下接收平台通信过程中，光子传输性能进行仿真的方法。

背景技术

在激光通信领域中，由于大气、海水的不稳定性，光传输过程中涉及的各种环境参数具有很大的随机性，而且光在大气和海水中传输由于多次散射的影响，几乎不可能用解析法描述光信号的传输过程。此外，星载激光对水下接收机通信试验工具和试验费用高，并且在实验室环境下，一些参数选择受实验条件限制(如海水和大气各层衰减系数等)，完全通过实验方法获取所有的数据也是不现实的。因此，目前是利用蒙特卡罗(Monte Carlo)方法对激光在大气和海水中的传输情况进行仿真。

蒙特卡罗方法，又称随机抽样方法，是一种基于“随机数”的计算方法，它利用随机数进行统计试验，以求得的统计特征值作为待解问题的数值解。具体地，光子在大气或海水中发生散射的过程，可视为光子的随机运动过程，其过程完全可以利用蒙特卡罗方法正确地模拟。蒙特卡罗方法用于光子运动轨迹的计算机模拟时，光子以一个特定方向进入随机介质，通过散射或吸收系数确定光子下一步的运动距离。如果散射现象发生，可以由适当的散射相位函数选取散射后光子新的运动方向。重复进行以上过程，直到光子被完全吸收或从边界逃逸为止。

目前，国内的许多高校如华中科技大学、电子科技大学、长春理工大学、桂林电子科技大学等开展了蓝绿激光的应用研究。但是，现有研究资料采用蒙特卡罗方法对星载或者机载激光对水下接收平台的通信仿真中，均无提到遇到边界条件时所产生的误差。

发明内容

本发明针对上述问题，提出了一种光子传输性能仿真方法。本发明采用的技术手段如下：

一种光子传输性能仿真方法，其特征在于包括：

步骤1：设定待仿真光子数目、各光子的初始状态和介质光学特性参数。初始状态包括光子的坐标r(X，Y，Z)，光子的运动方向ψ，光子运动路径长度l，光子能量权值W；介质光学特性参数包括介质吸收系数μ_a、介质散射系数μ_s；

步骤2：根据初始状态发射光子，对每一光子传输路径进行追踪，追踪过程中，对经过介质界面或接收平面的光子采用修正的蒙特卡罗方法实现追踪：

步骤3：基于追踪结果进行统计计算，生成光斑分布图、多径时延概率图、最大多径时延、63％光子平均多径时延和最大时延、光子所剩数目并显示。

本发明通过建立激光对水下接收平台通信传输的蒙特卡罗模型，统计接收面的光子信息，输出光子到达接收平面的多径时间延迟和光斑分布信息，模拟接收面光斑接收图以及多径时延图，其优点是针对光子运动路径上的介质界面及接收平面，提出了修正的蒙特卡罗方法，使得仿真结果更加准确，为进一步开发激光对水下接收平台通信系统提供了更精确的仿真数据。

附图说明

图1为光子传输性能仿真方法流程图。

图2为532nm绿光水下30米采用修正方法时接收面光斑分布情况。

图3为532nm绿光水下30米不采用修正方法时接收面光斑分布情况。

图4为532nm绿光水下30米采用修正方法时多径时延概率分布。

图5为532nm绿光水下30米不采用修正方法时多径时延概率分布。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

如图1所示，本发明方法包括：

步骤1：设定待仿真光子数目、各光子的初始状态和介质光学特性参数。初始状态包括光子的坐标r(X，Y，Z)，光子的运动方向ψ，光子运动路径长度l，光子能量权值W；介质光学特性参数包括介质吸收系数μ_a、介质散射系数μ_s。

步骤2：根据初始状态发射光子，对每一光子传输路径进行追踪，追踪过程中，对经过介质界面或接收平面的光子采用修正的蒙特卡罗方法实现追踪。

由于光子运动的步长很有可能使光子逃逸大气-海水界面或者接收平面，即遇到边界的情况下，会对仿真结果造成巨大的误差，因此采用修正的蒙特卡罗方法对经过介质界面或接收平面的光子进行追踪。修正的蒙特卡罗方法包括修正的光子步长公式：d′＝d_v/cosθ，修正的多径时延公式：

修正的光子行走被吸收的概率：

d_v是光子前坐标为与大气-水界面的垂直距离，散射角θ为光子当前坐标方向与接收平面法向方向的夹角，d′是修正后的光子步长，d是遇到边界条件时的光子步长，Δt′是修正后的多径时延，Δt是遇到边界条件时的多径时延，p′是修正后的光子行走被吸收的概率，p是遇到边界条件时光子行走被吸收的概率。

步骤3：基于追踪结果进行统计计算，生成光斑分布图、多径时延概率图、最大多径时延、63％光子平均多径时延和最大时延、光子所剩数目并显示。

以光子的传输路径是由大气进入水面下为例，步骤2又包括：

步骤21：发射光子，计算光子在大气运动中每一步移动的步长，并计算光子运动产生的多径时延和能量损失。

在发射光子时，光子初始位置为仿真坐标原点，并服从高斯分布：设l＝cotγ，γ为激光发散半角，求得光子包的初始方向为：

$\left\{\begin{array}{c}{X}_0=\frac{x}{\sqrt{x^2+y^2+l^2}}\\ Y_0=\frac{y}{\sqrt{x^2+y^2+l^2}}\\ Z_0=\frac{l}{\sqrt{x^2+y^2+l^2}}\end{array}\right.$

在计算光子在大气运动中每一步移动的步长时，首先计算光子衰减系数μ_t满足μ_t＝μ_a+μ_s，则根据Beer定理确定随机步长的表达式为：d＝-lnξ/μ_t，ξ是(0，1)之间均匀分布的随机数。

在计算光子在每次行走中被吸收的概率p时，满足：

$p=\frac{{\mathrm{\μ}}_a}{{\mathrm{\μ}}_a+{\mathrm{\μ}}_s}=1-\frac{{\mathrm{\μ}}_s}{{\mathrm{\μ}}_a+{\mathrm{\μ}}_s}=1-a.$

在计算光子能量损失时，假设光子行走n个步长后的剩余能量是W_n，满足：W_n＝W×aⁿ。当光子的能量W_n小于给定的阈值时，表示光子被完全吸收，光子行走结束。

假设光子在第n-1次散射发生前的坐标方向是(X_n-1，Y_n-1，Z_n-1)，当光子发生散射后，光子以新的方向角传输，假设散射后光子坐标方向是(X_n，Y_n，Z_n)，方位角

表示光子散射方向在水平面投影的旋转角度，则有：

步骤22：在光子到达大气-水界面时，采用修正的蒙特卡罗方法计算光子步长、多径时延和光子行走被吸收的概率，并由光子行走被吸收的概率得到光子的能量损失。

步骤23：计算光子在大气-水界面发生折射后的传输方向。

当光子发生水面折射时，水面受风吹动形成海浪，改变了光线的传播方向。沿用前面的直角坐标系，用矢量表示海浪表面的法线方向n＝(n_x，n_y，n_z)，其中， $\left\{\begin{array}{c}{n}_x=\cos \mathrm{\α}\sin \mathrm{\β}\\ n_y=\sin \mathrm{\α}\sin \mathrm{\β}\\ n_z=\cos \mathrm{\β}\end{array}\right.,$ α是折射方位角，它在(0，2π)内服从均匀分布；β为俯仰角，是法线方向与z轴之间的夹角。现有技术根据观测得出海浪在一定风速下俯仰角的经验归一化概率密度p(β)表示为：

式中，σ＝(0.003+0.00512V)^0.5，V为风速(m/s)。

对概率密度p(β)进行采样可得到相应俯仰角β，表示为：

ξ在(0，1)内服从均匀分布。

已知海浪的法线方向n＝(n_x，n_y，n_z)，假设入射光线方向为R_i＝(X_n，Y_n，Z_n)，经海水折射后的方向为

设入射角为θ_i，折射角为θ_o，已知大气的折射率n₁＝1，海水的折射率n₂＝1.33，根据几何关系及斯涅尔定理求入射角θ_i及折射角θ_o：

由

|n₁|＝1，|r_i|＝1，求得θ_i的标量形式：θ_i＝arccos(n_xX_n+n_yY_n+n_zZ_n)；由n_i sinθ_i＝n_o sinθ_o求得：

根据光学知识可知，入射光线、法线、折射光线在同一平面上，折射角即法线与折射光线的夹角，用公式表示为： $\left(\begin{array}{ccc}{X}_n& Y_n& Z_n\\ n_x& n_y& n_z\\ X_n^{\′}& Y_n^{\′}& Z_n^{\′}\end{array}\right)=0,$ $\cos {\mathrm{\θ}}_o=n_x{X}_n^{\′}+n_y{Y}_n^{\′}+n_z{Z}_n^{\′},$ $X_n^{\′2}+Y_n^{\′2}+Z_n^{\′2}=1,$ 由该三个公式可求得折射后的坐标方向令：

$w_1=\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}$

$w_2=\frac{X_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}$

$m=\frac{\cos {\mathrm{\θ}}_o}{n_x+n_z\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}}$

$p=-\frac{n_y-n_z\frac{X_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}}{n_x+n_z\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}}$

$l=\frac{\cos {\mathrm{\θ}}_o}{n_x+n_z\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}}\·\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}$

$q=\frac{Y_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}p-\frac{X_n{n}_z-n_y{Z}_n}{Y_n{n}_x-n_y{X}_n}$

则折射方向满足：

$Y_{n\mathrm{1,2}}^{\′}=-\frac{\mathrm{lq}+\mathrm{mp}}{1+p^2+q^2}\±\frac{\sqrt{{(\mathrm{mp}+\mathrm{lq})}^2-(1+p^2+q^2)(l^2+m^2-1)}}{1+p^2+q^2}$

$X_{n\mathrm{1,2}}^{\′}=m+p{Y}_{n\mathrm{1,2}}^{\′}$

$Z_{n\mathrm{1,2}}^{\′}=w_1{X}_{n\mathrm{1,2}}^{\′}-w_2{Y}_{n\mathrm{1,2}}^{\′}$

根据折射定律，海水的折射率大于大气的折射率，所以海水中的折射角小于大气中的入射角。分别用

与r_i的夹角相比较，夹角绝对值小的为正确的折射方向。

步骤24：计算光子在水中运动时的衰减系数，该衰减系数随海水深度变化而变化，采用叶绿素浓度对衰减系数影响的模型：

$K\left(\mathrm{\λ}\right)=A_w\left(\mathrm{\λ}\right)+A_c^0\left(\mathrm{\λ}\right){\left(\frac{C_c}{C_c^0}\right)}^{0.602}+B_w\left(\mathrm{\λ}\right)+\left(\frac{550}{\mathrm{\λ}}\right)0.3C_c^{0.602}$

式中，K(λ)是波长为λ的光的衰减系数，A_w(λ)是光吸收系数，B_w(λ)是光散射系数，为叶绿素的光谱吸收系数，C_c是叶绿素总浓度(单位为mg/m³)，常量

等于1mg/m³。

步骤25：在光子到达接收平面时，采用修正的蒙特卡罗方法计算光子步长、多径时延和光子行走被吸收的概率，重复上述步骤，直到追踪完全部光子。

步骤3又包括：

步骤31：将到达接收面的光子数进行统计，得出接收面光子数；

步骤32：将到达接收面的光子的光斑分布进行统计，得到光斑分布图；

步骤33：将所有光子多径时延进行统计，得到最大平均时延；

步骤34：将所有光子多径时延的总和除以所有光子数，得到平均多径时延；

步骤35：统计大约63％的光子落在的区域定义为63％半径，得到63％光子最大多径时延以及63％光子平均时延；

步骤36：将到达接收面的光子的能量权值进行相加，得到接收面光子权值总和；

步骤37：将63％光子的多径时延进行统计并按概率排列，得到多径时延概率分布图。

为了更好说明上述方法的效果，如下表一是532nm绿光采用修正公式以及不采用修正公式时相关仿真数据：

表一

如下表二是美国“魔灯”机载激光探雷系统主要性能：

表二

工作物质	Nd:YAG(倍频)
		工作波长	532nm(绿光)
重复率	大于或等于5kHz
		脉冲宽度	小于30ns

脉冲能量	大于2mJ
		帧速	大于5kHz
分辨率	0.0076-0.0203m
		工作水深	24.38m
探测概率	0.9-0.95
		误警	小于1次/h
工作环境	适应2-4级海况

如下表三是采用本方法推算美国“魔灯”机载激光探雷系统的水下光接收性能数据：

表三

通过表一、表三以及图2至图5的比较说明，在模拟光传输的环境已考虑美国机载激光探雷系统最坏环境的情况下，采用修正方法的仿真结果与现有美国机载激光探雷系统的工作数据更加接近，没有采用修正方法的模拟数据在一定程度上忽略了边界条件产生的影响。故修正方法在一定程度上改进了原有蒙特卡罗方法。

本发明通过建立激光对水下接收平台通信传输的蒙特卡罗模型，统计接收面的光子信息，输出光子到达接收平面的多径时间延迟和光斑分布信息，模拟接收面光斑接收图以及多径时延图，其优点是针对光子运动路径上的介质界面及接收平面，提出了修正的蒙特卡罗方法，使得仿真结果更加准确，为进一步开发激光对水下接收平台通信系统提供了更精确的仿真数据。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (4)

1.一种光子传输性能仿真方法，其特征在于包括：

步骤1：设定待仿真光子数目、各光子的初始状态和介质光学特性参数，初始状态包括光子的坐标r(X，Y，Z)，光子的运动方向ψ，光子运动路径长度l，光子能量权值W，介质光学特性参数包括介质吸收系数μ_a、介质散射系数μ_s；

步骤2：根据初始状态发射光子，对每一光子传输路径进行追踪，追踪过程中，对经过介质界面或接收平面的光子采用修正的蒙特卡罗方法实现追踪：

步骤3：基于追踪结果进行统计计算，生成光斑分布图、多径时延概率图、最大多径时延、63％光子平均多径时延和最大时延、光子所剩数目并显示。

2.根据权利要求1所述的方法，其特征在于步骤2又包括：

步骤21：发射光子，计算光子在大气运动中每一步移动的步长，并计算光子运动产生的多径时延和能量损失；

步骤22：在光子到达大气-水界面时，采用修正的蒙特卡罗方法计算光子步长、多径时延和光子行走被吸收的概率，并由光子行走被吸收的概率得到光子的能量损失；

步骤23：计算光子在大气-水界面发生折射后的传输方向；

步骤24：计算光子在水中运动时的衰减系数；

步骤25：在光子到达接收平面时，采用修正的蒙特卡罗方法计算光子步长、多径时延和光子行走被吸收的概率，重复上述步骤，直到追踪完全部光子。

3.根据权利要求2所述的方法，其特征在于修正的蒙特卡罗方法包括：修正的光子步长公式：d′＝d_v/cosθ，修正的多径时延公式：

修正的光子行走被吸收的概率：

d_v是光子前坐标为与大气-水界面的垂直距离，d′是修正后的光子步长，θ是光子当前坐标方向与接收平面法向方向的夹角，Δt′是修正后的多径时延，Δt是遇到边界条件时的多径时延，d是遇到边界条件时的光子步长，p′是修正后的光子行走被吸收的概率，p是遇到边界条件时光子行走被吸收的概率。

4.根据权利要求2或3所述的方法，其特征在于步骤3又包括：

步骤31：将到达接收面的光子数进行统计，得出接收面光子数；

步骤32：将到达接收面的光子的光斑分布进行统计，得到光斑分布图；

步骤33：将所有光子多径时延进行统计，得到最大平均时延；

步骤34：将所有光子多径时延的总和除以所有光子数，得到平均多径时延；

步骤35：统计大约63％的光子落在的区域定义为63％半径，得到63％光子最大多径时延以及63％光子平均时延；

步骤36：将到达接收面的光子的能量权值进行相加，得到接收面光子权值总和；

步骤37：将63％光子的多径时延进行统计并按概率排列，得到多径时延概率分布图。

Images