CN109856092A - 一种海洋荧光探测半解析蒙特卡罗模拟方法 - Google Patents

Abstract

本发明提供了一种海洋荧光探测半解析蒙特卡罗模拟方法，该方法基于荧光物质垂向密度分布，利用蒙特卡罗方法生成荧光物质光源，模拟荧光光子在海水随机介质中的传输，模型考虑海水和水中悬浮物质的散射和吸收特性；结合半解析方法，快速计算每次散射过程中水上探测器感应的能量，统计每个光子的光程并形成时域探测信号。本发明可用于分析不同水体的水色信息，也有助于开展海洋光学散射特性的研究，节省了大量成本，从而为遥感探测和抗干扰算法设计、参数选择、性能评估提供参考。

Description

一种海洋荧光探测半解析蒙特卡罗模拟方法

技术领域

本发明涉及海洋环境光散射特性研究领域，特别涉及一种运用半解析蒙特卡罗方法的海洋荧光传输特性仿真技术。

背景技术

在海洋光学领域，水下的荧光效应与海洋溶解物质的特性息息相关，利用荧光光谱特性可以反演出丰富的水质信息，对于不同来源的水团示踪，生物活动示踪，水质污染分等应用可起到至关重要的作用。海洋荧光探测模拟技术可以为海洋水色遥感和水下探测的数据分析、特征提取、目标识别等研究提供理论依据。

对于光子在媒介中的传输问题，一些国际上科研技机构已经进行了有效的研究和开发，如University of Texas开发的MCML、University of California 开发的MCNP、National Research Council of Canada联合开发的EGS4等。虽然这些软件可以模拟粒子输运问题，但是均存在一些不足，MCML主要用于模拟外部准直光源发射的光子在多层平板型3D媒介中的传输，采用的光源模型和媒介模型都是简单的经典模型；MCNP可用于模拟多种粒子，如中子、光子、电子在3D媒介中的传输，但MCNP没有直接模拟随机光源参数的粒子在2D媒介中传输的功能；EGS4是模拟电子和光子在不同媒介中传输的软件开发包，但也缺乏直接模拟粒子在2D媒介中传输的功能。虽然这些光子传输特性仿真技术在生物光学等领域被广泛应用，但不能快速、有效地模拟大范围海洋荧光的传输过程。

在检索到的国内外公开及有限范围发表的文献中，传统的蒙特卡罗方法中探测器接收到的光子概率低，结果难以收敛，而对于静态扩散方程的解析方法，复杂形状媒介的解析解通常难以求得。

发明内容

本发明的目的在于提出一种海洋荧光探测半解析蒙特卡罗模拟方法，既可用于分析不同水体的水色信息，也有助于开展海洋光学散射特性的研究，节省了大量成本，从而为遥感探测和抗干扰算法设计、参数选择、性能评估提供参考。

为了达到上述目的，本发明的技术方案是提供一种海洋荧光探测半解析蒙特卡罗模拟方法，用于实现荧光物质垂向分布不均匀的水体受光激发的荧光能量传输特性仿真。

所述方法，包含以下过程：

步骤1：将海水划分为深度等间距的光子多层传输媒介，归一化海水垂向深度的荧光物质分布谱，并将连续分布谱转化为等间距深度采样的离散概率分布；

步骤2：统计实例海域中荧光物质的平均浓度，依据荧光物质浓度与激发光波长的经验公式关系计算受光激发的荧光峰波长；基于该荧光峰波长，利用Mie理论球形粒子近似，计算海水各层媒介的固有光学特性，包括散射系数、吸收系数、不对称因子及折射率；

步骤3：基于荧光物质离散概率分布，利用蒙特卡罗方法随机生成荧光光源；考虑水下悬浮粒子及水体的散射、吸收效应，不同折射率的媒介层直接的反射、透射效应，模拟荧光光子的传输过程；结合解析法，计算每次散射后，荧光光子被水上探测器感应的能量；

步骤4：基于光子输运过程，进行统计运算，转化为其他所需的传输信息。输出探测器收集的荧光光子能量，生成时域探测信号。

本发明带来以下有益效果：

本发明提出半解析蒙特卡罗方法的海洋荧光探测仿真方法，相比传统的蒙特卡罗光子输运模拟，实现了在任意不均匀垂向分布的荧光物质辐射传输特性的快速模拟。解决了海洋遥感和探测对水下荧光辐射特性的需要，和海上实验相比节省了大量的人力、物力，能够有效地为海洋水色遥感、水质分析和探测抗干扰算法设计、参数选择、性能评估提供参考。

附图说明

图1是本发明中的海洋荧光探测半解析蒙特卡罗方法仿真流程图。

具体实施方式

本发明提供了一种运用荧光物质垂向分布谱来模拟荧光光源的半解析蒙特卡罗仿真方法，为了使本发明更急明显易懂，以下结合附图和具体实施方式对本发明做进一步说明。

仿真的具体步骤如下：

步骤1：以海平面为基准，将基于海水垂向深度的荧光物质分布谱归一化，并转化为等间距的离散概率分布谱，其中每单位间距水深的海水视为一层光子传输媒介，在层内的固有光学特性如荧光物质浓度C、吸收系数σ_a、散射系数σ_b不随着深度变化。

步骤2：计算荧光物质受太阳光激发的荧光峰，根据Gitelson的经验回归公式得到荧光峰

λ_peak＝a+b·C_chl (13)

其中C_chl为该海域荧光物质平均浓度，a和b为经验参数，依据荧光物质不同而变化。在本实施例中，荧光物质来源于海洋中的叶绿素，a在661.7- 666.4nm之间变化，b在0.001-0.03之间变化。

根据荧光峰波长，基于Mie理论球形近似，计算光子在各层传输媒介中的吸收系数σ_a、散射系数σ_b、折射率n、不对称因子g。

步骤3：基于荧光物质的垂向离散分布谱，生成不同位置的荧光光源。通过蒙特卡罗方法模拟光子在海水中的输运过程，结合半解析方法计算光子被水上探测器感应到的能量。

步骤3具体包括以下步骤：

步骤3.1：荧光光源位置的抽样。通过蒙特卡罗直接抽样法，生成[0,1] 区间均匀分布随机数ξ₁，判断随机数落在的概率区间

Z_F＝Z_I,当

其中P_i为第i层海水传输媒介的荧光物质分布概率密度，Z_I为第I层海水的上边界深度。进一步地，假设光源在层内服从均匀分布，直接抽样得到光子的出射位置(x,y,z)如下

其中ξ_i是[0,1]区间均匀分布的随机数。

步骤3.2：荧光光子传输方向的抽样。光子的传输方向由方位角和偏向角θ表示，在各向同性介质中通过蒙特卡罗采样得到

转化为三维向量形式

步骤3.3：解析法计算光子直接被探测的能量。探测器位于水平面上方 M处(x_m,y_m,z_m)，忽然多次折射效应，荧光光子不经过散射直接被探测的能量为

其中E_j为第j次散射时光子剩余能量，初始能量被设定为1，f(α)为 Henyey-Greenstein相函数，为光子坐标指向探测器坐标的归一化向量，μ_ti为第i层的消光系数，L_i为第i层的传输路径长度，ΔS为探测器感应面积，d_j为探测器到光子的垂直距离，T_i,i-1(θ)为基于Fresnel方程的从i到第i-1层的光子透过率。统计完直接探测的能量后进入步骤3.4。

步骤3.4：生成光子传输步长。首先判断光子是否有剩余步长S_left，如果有则将剩余步长赋值给光子传输步长，否则根据光子在介质中的自由程概率密度函数计算，可以得到

进一步，判断在现有的光子传输方向下，以该传输步长移动是否会碰撞媒介层边界。如果光子不会碰到边界进入步骤步骤3.5，否则进入步骤3.7。

步骤3.5：如果光子不会碰到边界，则以该步长和方向更新光子的坐标

步骤3.6：光子吸收和散射。光子经过传输步长到一个新位置后，光子的能量会被吸收，剩余光子能量为

同时，光子的传输方向受到散射点影响也会改变。其中方位角在[0,2π]中均匀分布，偏转角服从Henyey-Greenstein相函数概率分布，重复步骤3.2 更新光子的传输方向，并进入步骤3.10。

步骤3.7：如果光子会碰到边界，则根据Fresnel方程判断光子在边界会发生反射或是透射过程：

其中，θ_i和θ_t分别为光子的入射角和透射角；透射角可以由Snell公式计算。

取[0,1]区间均匀分布的随机数ξ，如果ξ≤R(θ)，光子发生反射，否则发生透射。然后，根据Snell公式更新光子运动方向和光子所在媒介层。

如果发生透射进入步骤3.8，发生反射则进入步骤3.9。

步骤3.8：判断光子所在位置是否位于多层媒介的最外层，如果是则进入步骤3.11，否则进入步骤3.9

步骤3.9：将光子依据传输方向移动步长s'到该层媒介的边界，并计算剩余步长

s_left＝s-s' (23)

步骤3.10：判断光子是否死亡。如果光子死亡则进入步骤3.11，否则进一步判断光子剩余能量是否大于阈值。如果光子能量大于阈值则进入步骤 3.3，否则进行Russianroulette轮盘判定光子是否继续输运过程，如果继续光子输运则进入步骤3.3，否则进入步骤3.11。

步骤3.11：判断现在处理的光子是否为程序设定的最后一个光子，如果还有光子剩余，则重复步骤3.1-1.11，否则程序结束。

步骤4：基于光子输运过程，进行统计运算，输出探测器收集的荧光光子能量，将每个光子的输运总光程转化成光子回波时间，生成时域探测信号。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (7)

1.一种海洋荧光探测半解析蒙特卡罗模拟方法，用于实现荧光物质垂向分布不均匀的水体受光激发的荧光能量传输特性仿真，其特征在于，所述方法，包含以下过程：

步骤1：将海水划分为深度等间距的光子多层传输媒介，归一化海水垂向深度的荧光物质分布谱，并将连续分布谱转化为等间距深度采样的离散概率分布；

步骤2：统计实例海域中荧光物质的平均浓度，依据荧光物质浓度与激发光波长的经验公式关系计算受光激发的荧光峰波长；基于该荧光峰波长，利用Mie理论球形粒子近似，计算海水各层媒介的固有光学特性，包括散射系数、吸收系数、不对称因子及折射率；

步骤3：基于荧光物质离散概率分布，利用蒙特卡罗方法随机生成荧光光源；考虑水下悬浮粒子及水体的散射、吸收效应，不同折射率的媒介层直接的反射、透射效应，模拟荧光光子的传输过程；结合解析法，计算每次散射后，荧光光子被水上探测器感应的能量；

步骤4：基于光子输运过程，进行统计运算，转化为其他所需的传输信息；输出探测器收集的荧光光子能量，生成时域探测信号。

2.如权利要求1所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

步骤1中，荧光物质分布谱为等深度间距采样的离散概率分布谱；其中每单位间距水深的海水视为一层光子传输媒介，在层内的固有光学特性包含荧光物质浓度C、吸收系数σ_a、散射系数σ_b不随着深度变化。

3.如权利要求2所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

步骤2中，海水被处理为多层光子传输媒介，其固有光学特性依据荧光峰波长计算。

4.如权利要求3所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

步骤2：计算荧光物质受太阳光激发的荧光峰，根据Gitelson的经验回归公式得到荧光峰

λ_peak＝a+b·C_chl (1)

其中C_chl为该海域荧光物质平均浓度，a和b为经验参数，依据荧光物质不同而变化。

5.如权利要求4所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

荧光物质来源于海洋中的叶绿素时，a在661.7-666.4nm之间变化，b在0.001-0.03之间变化。

6.如权利要求3所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

步骤3中，蒙特卡罗方法模拟每次光子散射后，水上探测器直接感应的能量的半解析计算公式为

其中E_j为第j次散射时光子剩余能量，初始能量被设定为1，f(α)为Henyey-Greenstein相函数，为光子坐标指向探测器坐标的归一化向量，μ_ti为第i层媒介的消光系数，L_i为第i层的传输路径长度，ΔS为探测器感应面积，d_j为探测器到光子的垂直距离，T_i,i-1(θ)为基于Fresnel方程的从i到第i-1层的光子透过率。

7.如权利要求3所述的海洋荧光探测半解析蒙特卡罗模拟方法，其特征在于：

步骤3进一步包括以下步骤：

步骤3.1：荧光光源位置的抽样；通过蒙特卡罗直接抽样法，生成[0,1]区间均匀分布随机数ξ₁，判断随机数落在的概率区间

其中P_i为第i层海水传输媒介的荧光物质分布概率密度，Z_I为第I层海水的上边界深度；假设光源在层内服从均匀分布，直接抽样得到光子的出射位置(x,y,z)如下

其中ξ_i是[0,1]区间均匀分布的随机数；

步骤3.2：荧光光子传输方向的抽样；光子的传输方向由方位角和偏向角θ表示，在各向同性介质中通过蒙特卡罗采样得到

转化为三维向量形式

步骤3.3：解析法计算光子直接被探测的能量；探测器位于水平面上方M处(x_m,y_m,z_m)，忽然多次折射效应，荧光光子不经过散射直接被探测的能量为

其中E_j为第j次散射时光子剩余能量，初始能量被设定为1，f(α)为Henyey-Greenstein相函数，为光子坐标指向探测器坐标的归一化向量，μ_ti为第i层的消光系数，L_i为第i层的传输路径长度，ΔS为探测器感应面积，d_j为探测器到光子的垂直距离，T_i,i-1(θ)为基于Fresnel方程的从i到第i-1层的光子透过率；统计完直接探测的能量后进入步骤3.4；

步骤3.4：生成光子传输步长；首先判断光子是否有剩余步长S_left，如果有则将剩余步长赋值给光子传输步长，否则根据光子在介质中的自由程概率密度函数计算，可以得到

进一步，判断在现有的光子传输方向下，以该传输步长移动是否会碰撞媒介层边界；如果光子不会碰到边界进入步骤步骤3.5，否则进入步骤3.7；

步骤3.5：如果光子不会碰到边界，则以该步长和方向更新光子的坐标

步骤3.6：光子吸收和散射；光子经过传输步长到一个新位置后，光子的能量会被吸收，剩余光子能量为

同时，光子的传输方向受到散射点影响也会改变；其中方位角在[0,2π]中均匀分布，偏转角服从Henyey-Greenstein相函数概率分布，重复步骤3.2更新光子的传输方向，并进入步骤3.10；

步骤3.7：如果光子会碰到边界，则根据Fresnel方程判断光子在边界会发生反射或是透射过程：

其中，θ_i和θ_t分别为光子的入射角和透射角；透射角可以由Snell公式计算；

取[0,1]区间均匀分布的随机数ξ，如果ξ≤R(θ)，光子发生反射，否则发生透射；然后，根据Snell公式更新光子运动方向和光子所在媒介层；

如果发生透射进入步骤3.8，发生反射则进入步骤3.9；

步骤3.8：判断光子所在位置是否位于多层媒介的最外层，如果是则进入步骤3.11，否则进入步骤3.9

步骤3.9：将光子依据传输方向移动步长s'到该层媒介的边界，并计算剩余步长

s_left＝s-s' (12)

步骤3.10：判断光子是否死亡；如果光子死亡则进入步骤3.11，否则进一步判断光子剩余能量是否大于阈值；如果光子能量大于阈值则进入步骤3.3，否则进行Russian roulette轮盘判定光子是否继续输运过程，如果继续光子输运则进入步骤3.3，否则进入步骤3.11；

步骤3.11：判断现在处理的光子是否为程序设定的最后一个光子，如果还有光子剩余，则重复步骤3.1-1.11，否则程序结束。