CN113361080A - 基于gpu的多层水体光子传输半解析蒙特卡洛仿真方法 - Google Patents

Abstract

本发明提供一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，包括：在CPU的内存中分配固定的输入输出数据储存空间；在GPU的显存中分配固定的输入输出数据储存空间；将在CPU端赋初值的参数数组输入GPU；在GPU端进行并行光子计算；在GPU中对接收到的能量进行记录，并进行累加；所有光子都湮灭或逸出边界后，将记录数组从GPU输出到CPU端，使用CPU对记录数组进行统计绘图，输出不同散射次数的回波信号强度与深度的变化情况，获得最终仿真结果。本发明通过对多层水体光子传输半解析蒙特卡洛模型使用GPU并行运算，提高了多层水体光子传输半解析蒙特卡洛仿真的计算速度。

Description

基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法

技术领域

本发明涉及激光雷达技术领域，尤其涉及一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法。

背景技术

蒙特卡洛方法是一种常用的数学方法，可以用于仿真多种物理过程。用于弹性激光雷达多次散射模拟的蒙特卡洛模型假设光子只有粒子性而忽略光的波动性，认为激光在水体中的传输路径由许多完全随机的光子轨迹构成。对于半解析蒙特卡洛模型而言，其与标准蒙特卡洛模型主要的不同在于在半解析蒙特卡洛模型中，能量的接收发生在所有满足激光雷达接收器的孔径和视场角条件的散射事件中，而不是最终能够返回接收器的非常少的光子，因此，半解析蒙特卡洛模型显著降低了数据的统计不确定性，使用较少的光子就能获得较高的仿真可靠性。

目前多层水体光子传输半解析蒙特卡洛仿真主要采用CPU进行运算，虽然基于CPU的仿真算法成熟，可选的编程环境也很多，但是由于CPU更注重指令运算，并行运算能力不足，所以采用CPU对多层水体光子传输半解析蒙特卡洛模型进行仿真往往需要耗费大量的时间。而GPU以图形类数值计算为核心，虽然GPU执行每个数值计算的速度并没有比CPU快，但是GPU可以容纳上千个数值计算线程，这使GPU相较于CPU更适合处理数据密集型的计算任务。多层水体光子传输半解析蒙特卡洛模型中包含大量数据密集型的计算，因此，使用GPU对多层水体光子传输半解析蒙特卡洛模型进行仿真能够提高模型仿真速度。目前民用级别的GTX1080显卡，可实现9T次每秒的浮点运算能力，比intel core i7运算能力高10倍以上，因此，基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法具有很大优势。

发明内容

本发明在于克服现有技术的不足，提供一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，通过对多层水体光子传输半解析蒙特卡洛模型使用GPU并行运算，提高了多层水体光子传输半解析蒙特卡洛仿真的计算速度，同时，使用CPU对GPU输出的记录数组进行统计，可以获得不同散射次数的回波信号强度与深度的变化情况。

本发明采用如下技术方案：

一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，包括：

S101，在CPU的内存中分配固定的输入输出数据储存空间；

S102，在GPU的显存中分配固定的输入输出数据储存空间；

S103，在CPU端对输入参数赋初值，并将初值传入GPU端；所述输入参数包括水体状态参数、光子初始状态参数和其它配置参数；

S104，在GPU端进行并行光子计算；

S105，记录接收能量及接收深度；

S106，判断所有光子是否都完成模拟，如果否，重复对所述输入参数赋初值，并将初值传入GPU端，如果是，累计记录结果并传回CPU端；

S107，在CPU端对记录结果进行统计绘图，输出不同散射次数的回波信号强度与深度的变化情况，获得最终仿真结果。

优选的，在CPU的内存中分配输入输出数据储存空间的函数为Malloc；在GPU的显存中分配输入输出数据储存空间的函数为cudaMalloc；将CPU端所赋初值传入GPU端所使用的函数为cudaMemcpy。

优选的，所述水体状态参数包括水体层厚、层数、各层的吸收系数和散射系数；所述光子初始状态参数包括光子位置、光子方向和光子能量；所述其它配置参数包括望远镜视场角大小、望远镜孔径、激光雷达离水高度和发射的光子总数。

优选的，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的能量大小并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

优选的，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

产生一个随机波长作为激发的荧光的波长，根据波长、吸收能量以及水体状态参数，计算荧光能量大小，根据荧光的散射相函数计算荧光的方向，将荧光光子的波长、方向、位置及能量入栈；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的波长和能量大小并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

判断栈是否为空，如果不为空，荧光光子的波长、方向、位置及能量出栈，继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

优选的，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

根据光子位置、接收器位置、入射光Stokes矢量以及水体状态参数，计算后向散射的Stokes矢量；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的能量大小和后向散射的Stokes矢量并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

优选的，所述随机步长和光子能量的计算方法包括：

S701，设置一个初始值为1的自然数变量i，判断

是否小于z，如果小于，对i进行累加操作，即i＝i+1；如果大于等于，输出此时的i值，即光子当前位置所处的层数；其中，z表示z方向的光子位置；

S702，将得到变量i值赋给变量j，判断μ_z>0，如果为真，执行S703；如果为假，执行S704；其中，μ_z表示当前位置z方向的方向余弦；

S703，判断

如果为假，令j＝j+1，如果为真，输出此时的j值，则光子随机步长的计算方式如下：

其中，ξ₁表示在(0,1)上均匀分布的随机数；c_j表示第j层水体上的光衰减系数；D_j表示第j层水体层厚；D_k表示第k层水体层厚；c_k表示第k层光衰减系数；

所述光子能量的计算方式如下：

其中，W表示初始光子权重；S表示随机步长；a_k表示第k层水体上的吸收系数；c_k表示第k层水体上的光衰减系数；a_j表示第j层水体上的吸收系数；

S703，判断

如果为假，j＝j-1，如果为真，输出此时的j值

则光子随机步长的计算方式如下：

所述光子能量的计算方式如下：

优选的，光子位置的更新方法如下：

其中，(x,y,z)表示光子的初始位置，(x′,y′,z′)表示光子更新后的位置，(μ_x,μ_y,μ_z)表示光子初始方向余弦。

优选的，光子方向的更新方法包括：

S901，计算方位角：

ψ＝2πξ₂

其中，ξ₂表示在(0,1)上均匀分布的随机数；

计算散射角：

其中，g表示水体各向异性常数；ξ₃表示在(0,1)上均匀分布的随机数；

S902，判断|μ_z|>0.99，如果为真，散射更新的方向余弦公式如下：

其中，(μ_x′,μ_y′,μ_z′)表示更新后的光子方向余弦；

如果为假，散射更新的方向余弦公式如下：

优选的，激光雷达接收器能接收到的能量大小的计算方法包括：

计算从光子当前位置前往接收器方向的方向余弦μ_bz，如下：

其中，H表示激光雷达离水高度；

计算光子前进方向与从光子当前位置前往接收器方向的夹角cosω，如下：

计算散射相函数p(ω)，如下：

计算激光雷达接收器能接收到的能量大小，如下：

其中，D_k表示第k层水体层厚；T_s表示水面透射率。

由上述本发明提供的技术方案可以看出，使用GPU代替CPU对多层水体光子传输半解析蒙特卡洛模型进行仿真，提高了该模型仿真的计算速度。同时，使用经过优化的程序架构，提高了GPU并行运算的计算效率，此外，使用CPU对GPU输出的记录数组进行统计，可以获得不同散射次数的回波信号强度与深度的变化情况，便于在此仿真方法的基础上开展其他研究工作。

根据下文结合附图对本发明具体实施例的详细描述，本领域技术人员将会更加明了本发明的上述及其他目的、优点和特征。

附图说明

图1为本发明实施例的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法的流程图；

图2为本发明实施例一的多层水体光子传输半解析蒙特卡洛模型流程图；

图3为本发明实施例一的多层水体光子传输半解析蒙特卡洛仿真方法的输出结果与基于CPU的仿真输出结果的比对；其中(a)为基于CPU的仿真输出结果，(b)为基于GPU的仿真输出结果；

图4为本发明实施例二的多层水体光子传输荧光半解析蒙特卡洛模型流程图；

图5为本发明实施例三的多层水体光子传输偏振半解析蒙特卡洛模型流程图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

2007年英伟达NVIDIA公司推出了统一计算设备架构(CUDA)，将显卡(GPU)变成可以用作数据并行计算的设备。经过十多年的发展，CUDA技术已经大规模商用，完全兼容C/C++等常见编程语言，具有高级而完备的软件开发工具和环境，且GPU广泛用于娱乐、高性能计算领域，容易采购且性价比高。目前民用级别的GTX1080显卡，单块不足5000元，可实现9T次每秒的浮点运算能力，比intel core i7运算能力高10倍以上，因此，基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法具有很大优势。本实施例中，使用CUDA编写光子单步长状态变化函数以及接收函数，并使用GPU多线程对这些函数进行并行计算。具体的，所述的单步长状态变化函数包括随机步长函数、散射函数、吸收函数以及位置变化函数，所述随机补充函数用于实现光子随机步长的计算，所述散射函数用于实现光子方向的更新，所述吸收函数用于实现光子能量的更新，所述位置变化函数用于实现光子位置的更新。所示接收函数用于实现接收器能接收到的能量大小。

实施例一

参见图1所示，本发明一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，包括：

S1，使用Malloc在CPU的内存中分配固定的输入输出数据储存空间；

S2，使用cudaMalloc在GPU的显存中分配固定的输入输出数据储存空间；

S3，使用cudaMemcpy将在CPU端赋初值的参数数组输入GPU；

S4，使用CUDA编写光子单步长状态变化函数以及接收函数，使用GPU多线程对这些函数进行并行计算；

S5，在GPU中对接收函数接收到的能量进行记录，并进行累加；

S6，所有光子都湮灭或逸出边界后，使用cudaMemcpy将记录数组从GPU输出到CPU端，使用CPU对记录数组进行统计绘图，输出不同散射次数的回波信号强度与深度的变化情况，获得最终仿真结果。

本实施例中，显存中分配固定的输入输出数据储存空间至少为两块。

参见图2所示，为本发明实施例一提供的多层水体光子传输半解析蒙特卡洛仿真模型的流程图，具体如下：

S201，输入水体状态参数、光子初始状态参数以及其它配置参数；如下表1。

表1

S202，计算随机步长，依据输入的水体参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

S203，判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算接收器能接收到的能量大小并记录，如果不满足，不进行记录；

S204，判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，光子湮灭或逸出边界，反之光子继续循环，进行下一个步长计算，直至湮灭或逸出边界；

S205，光子湮灭或逸出边界后判断所有光子是否都完成模拟，如果所有光子都完成，则模型运行结束，反之则继续发射下一个光子进入水中。

本实施例中，所述随机步长和光子能量的计算方法包括：

S701，设置一个初始值为1的自然数变量i，判断

是否小于z，如果小于，对i进行累加操作，即i＝i+1；如果大于等于，输出此时的i值，即光子当前位置所处的层数；其中，z表示z方向的光子位置；

S702，将得到变量i值赋给变量j，判断μ_z>0，如果为真，执行S703；如果为假，执行S704；其中，μ_z表示当前位置z方向的方向余弦；

S703，判断

如果为假，令j＝j+1，如果为真，输出此时的j值，则光子随机步长的计算方式如下：

其中，ξ₁表示在(0,1)上均匀分布的随机数；c_j表示第j层水体上的光衰减系数；D_j表示第j层水体层厚；D_k表示第k层水体层厚；c_k表示第k层光衰减系数；

所述光子能量的计算方式如下：

其中，W表示初始光子权重；S表示随机步长；a_k表示第k层水体上的吸收系数；c_k表示第k层水体上的光衰减系数；a_j表示第j层水体上的吸收系数；

S703，判断

如果为假，j＝j-1，如果为真，输出此时的j值

则光子随机步长的计算方式如下：

所述光子能量的计算方式如下：

进一步的，光子位置的更新方法如下(光子移动)：

其中，(x,y,z)表示光子的初始位置，(x′,y′,z′)表示光子更新后的位置，(μ_x,μ_y,μ_z)表示光子初始方向余弦。

进一步的，光子方向的更新方法包括(光子散射)：

S901，计算方位角：

ψ＝2πξ₂

其中，ξ₂表示在(0,1)上均匀分布的随机数；

计算散射角：

其中，g表示水体各向异性常数；ξ₃表示在(0,1)上均匀分布的随机数；

S902，判断|μ_z|>0.99，如果为真，散射更新的方向余弦公式如下：

其中，(μ_x′,μ_y′,μ_z′)表示更新后的光子方向余弦；

如果为假，散射更新的方向余弦公式如下：

更进一步的，光子到接收机(如接收望远镜)的直线与接收机(如接收望远镜)轴线之间的夹角α的计算方式如下：

判断α≤FOV(视场角)，如果为真，接收机接收光子部分能量W_E，记录W_E以及光子从发射到接收的总散射次数和总路程，此时光子剩余权重W_l需要将接收部分扣除。

W_E的计算方法包括：

计算μ_bz，如下：

其中，H表示激光雷达接收器距离水面的高度；

计算cosω，如下：

计算p(ω)，如下：

计算激光雷达接收器能接收到的能量大小，如下：

其中，D_k表示第k层水体层厚；T_s表示水面透射率。

则：

如果α>FOV，判断光子能量是否小于权重阈值，如果为真，光子湮灭，进行下一个光子的发射；如果为假，光子继续运动。

参见图3所示，为本发明的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法的输出结果与基于CPU的仿真输出结果的比对；其中(a)为基于CPU的仿真输出结果，(b)为基于GPU的仿真输出结果。

从图3可以看出，本实施例基于英伟达公司的统一计算架构(CUDA)算法，使用显卡(GPU)并行运算代替之前多层水体光子传输半解析蒙特卡洛仿真的CPU运算，使用CPU对GPU输出的记录数组进行统计和处理。相对于基于CPU的多层水体光子传输半解析蒙特卡洛仿真，使用GPU运算速度更快。使用CPU对GPU输出的记录数组进行统计，可以获得不同散射次数的回波信号强度与深度的变化情况，便于在此仿真方法的基础上开展其他研究工作。

实施例二

本实施例一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法参见图1所示。本实施例与实施例一的区别在于，所提供的多层水体光子传输半解析蒙特卡洛仿真模型有一定的区别，具体的，参见图4所示，模型包括：

S401，输入水体状态参数、光子初始状态参数以及其它配置参数；输入的水体状态参数、光子初始状态参数以及其它配置参数参见表1。

S402，计算随机步长，依据输入的水体参数计算光子因为吸收散射导致的状态参数的变化；

S403，产生一个随机波长作为激发的荧光的波长，根据波长、吸收能量以及水体状态参数，计算荧光能量大小，根据荧光的散射相函数计算荧光的方向，将荧光光子的波长、方向、位置及能量入栈；

S404，判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算接收器能接收到的波长和能量大小并记录，如果不满足，不进行记录；

S405，判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，光子湮灭，反之光子继续循环，进行下一个步长计算，直至湮灭或逸出边界；

S406，判断栈是否为空，如果不为空，荧光光子的波长、方向、位置及能量出栈，继续循环，反之荧光光子湮灭或逸出边界；

S407，光子湮灭或逸出边界后判断所有光子是否都完成模拟，如果所有光子都完成，则模型运行结束，反之则继续发射下一个光子进入水中。

具体的，所述随机步长和光子能量的计算方法、光子方向的更新方法、方向余弦的更新方法以及接收能量的计算方法参见实施例一，本实施例不再重复说明。

相对于基于CPU的多层水体光子传输半解析蒙特卡洛仿真，本实施例使用GPU运算速度更快。同时使用CPU对GPU输出的记录数组进行统计，可以获得不同散射次数的回波信号强度与深度的变化情况，便于在此仿真方法的基础上开展其他研究工作。

实施例三

本实施例一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法参见图1所示。本实施例与实施例一的区别在于，所提供的多层水体光子传输半解析蒙特卡洛仿真模型有一定的区别，具体的，参见图5所示，模型包括：

S501，输入水体状态参数、光子初始状态参数以及其它配置参数；

S502，计算随机步长，依据输入的水体参数计算光子因为吸收导致的光子包权重的变化；

S503，根据光子位置、接收机位置、入射光Stokes矢量以及水体状态参数，计算后向散射的Stokes矢量；

S504，判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算接收器能接收到光子包能量大小并记录，如果不满足，不进行记录；

S505，根据偏振的散射相函数和入射光Stokes矢量，产生散射的极角和方位角，计算散射后的Stokes矢量、光子包传输的方向。

S506，判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，光子湮灭或逸出边界，反之光子继续循环，进行下一个步长计算，直至湮灭或逸出边界；

S507，光子湮灭或逸出边界后判断所有光子是否都完成模拟，如果所有光子都完成，则模型运行结束，反之则继续发射下一个光子进入水中。

具体的，所述随机步长和光子能量的计算方法、光子方向的更新方法、方向余弦的更新方法以及接收能量的计算方法参见实施例一，本实施例不再重复说明。

相对于基于CPU的多层水体光子传输半解析蒙特卡洛仿真，本实施例使用GPU运算速度更快。同时使用CPU对GPU输出的记录数组进行统计，可以获得不同散射次数的回波信号强度与深度的变化情况，便于在此仿真方法的基础上开展其他研究工作。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (10)

1.一种基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，包括：

S101，在CPU的内存中分配固定的输入输出数据储存空间；

S102，在GPU的显存中分配固定的输入输出数据储存空间；

S103，在CPU端对输入参数赋初值，并将初值传入GPU端；所述输入参数包括水体状态参数、光子初始状态参数和其它配置参数；

S104，在GPU端进行并行光子计算；

S105，记录接收能量及接收深度；

S106，判断所有光子是否都完成模拟，如果否，重复对所述输入参数赋初值，并将初值传入GPU端，如果是，累计记录结果并传回CPU端；

S107，在CPU端对记录结果进行统计绘图，输出不同散射次数的回波信号强度与深度的变化情况，获得最终仿真结果。

2.根据权利要求1所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，在CPU的内存中分配输入输出数据储存空间；在GPU的显存中分配输入输出数据储存空间；将CPU端所赋初值传入GPU端。

3.根据权利要求1所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，所述水体状态参数包括水体层厚、层数、各层的吸收系数和散射系数；所述光子初始状态参数包括光子位置、光子方向和光子能量；所述其它配置参数包括望远镜视场角大小、望远镜孔径、激光雷达离水高度和发射的光子总数。

4.根据权利要求1所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的能量大小并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

5.根据权利要求1所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

产生一个随机波长作为激发的荧光的波长，根据波长、吸收能量以及水体状态参数，计算荧光能量大小，根据荧光的散射相函数计算荧光的方向，将荧光光子的波长、方向、位置及能量入栈；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的波长和能量大小并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

判断栈是否为空，如果不为空，荧光光子的波长、方向、位置及能量出栈，继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

6.根据权利要求1所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，所述S104～S108，具体包括：

发射光子到水中；

计算随机步长，依据输入的水体状态参数计算光子因为吸收散射导致的状态参数的变化；变化的状态参数包括光子位置、光子方向和光子能量；

根据光子位置、接收器位置、入射光Stokes矢量以及水体状态参数，计算后向散射的Stokes矢量；

判断此时的光子是否满足激光雷达接收器的孔径和视场角条件，如果满足，计算激光雷达接收器能接收到的能量大小和后向散射的Stokes矢量并记录；如果不满足，不进行记录；

判断此时的光子能量是否小于设定阈值，如果光子能量比阈值小，说明光子湮灭或逸出边界，反之光子继续循环，进行下一个随机步长计算，直至湮灭或逸出边界；

光子湮灭或逸出边界后判断是否所有光子都完成模拟，如果所有光子都完成，则累计记录结果并传回CPU端，反之则继续发射下一个光子进入水中。

7.根据权利要求4、5或6所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，所述随机步长和光子能量的计算方法包括：

S701，设置一个初始值为1的自然数变量i，判断

是否小于z，如果小于，对i进行累加操作，即i＝i+1；如果大于等于，输出此时的i值，即光子当前位置所处的层数；其中，z表示z方向的光子位置；

S702，将得到变量i值赋给变量j，判断μ_z>0，如果为真，执行S703；如果为假，执行S704；其中，μ_z表示当前位置z方向的方向余弦；

S703，判断

如果为假，令j＝j+1，如果为真，输出此时的j值，则光子随机步长的计算方式如下：

其中，ξ₁表示在(0,1)上均匀分布的随机数；c_j表示第j层水体上的光衰减系数；D_j表示第j层水体层厚；D_k表示第k层水体层厚；c_k表示第k层光衰减系数；

所述光子能量的计算方式如下：

其中，W表示初始光子权重；S表示随机步长；a_k表示第k层水体上的吸收系数；c_k表示第k层水体上的光衰减系数；a_j表示第j层水体上的吸收系数；

S703，判断

如果为假，j＝j-1，如果为真，输出此时的j值

则光子随机步长的计算方式如下：

所述光子能量的计算方式如下：

8.根据权利要求7所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，光子位置的更新方法如下：

其中，(x,y,z)表示光子的初始位置，(x′,y′,z′)表示光子更新后的位置，(μ_x,μ_y,μ_z)表示光子初始方向余弦。

9.根据权利要求8所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，光子方向的更新方法包括：

S901，计算方位角：

ψ＝2πξ₂

其中，ξ₂表示在(0,1)上均匀分布的随机数；

计算散射角：

其中，g表示水体各向异性常数；ξ₃表示在(0,1)上均匀分布的随机数；

S902，判断|μ_z|>0.99，如果为真，散射更新的方向余弦公式如下：

其中，(μ_x′,μ_y′,μ_z′)表示更新后的光子方向余弦；

如果为假，散射更新的方向余弦公式如下：

10.根据权利要求9所述的基于GPU的多层水体光子传输半解析蒙特卡洛仿真方法，其特征在于，激光雷达接收器能接收到的能量大小的计算方法包括：

计算从光子当前位置前往接收器方向的方向余弦μ_bz，如下：

其中，H表示激光雷达离水高度；

计算光子前进方向与从光子当前位置前往接收器方向的夹角cosω，如下：

计算散射相函数p(ω)，如下：

计算激光雷达接收器能接收到的能量大小，如下：

其中，D_k表示第k层水体层厚；T_s表示水面透射率。

Images