CN113029342A - 基于双向反射理论的水下目标反射偏振光的仿真方法 - Google Patents

Abstract

本发明提供一种基于双向反射理论的水下目标反射偏振光的仿真方法，解决目前多为基于特定环境下通过实验的方式获得先验知识并对图像进行复原，存在局限性，对不同环境先验知识不足的问题。包括：1)设置初始光源、散射环境以及反射目标的参数；2)根据初始光源参数，采样初始光子坐标、运动矢量、斯托克斯矢量；根据散射环境参数，计算衰减系数、散射系数、单次反照率；3)模拟光子移动，采样光子移动步长、能量权值；4)根据光子移动步长和能量权值，判断光子行为；5)根据接收到的光子的运动矢量，计算最终接收光子的斯托克斯矢量；6)对最终接收到的光子进行统计分析，获得设定环境参数下的光子传输特性以及偏振特性的分析结果。

Description

基于双向反射理论的水下目标反射偏振光的仿真方法

技术领域

本发明属于水下光学成像技术领域，具体涉及一种基于双向反射理论的水下目标反射偏振光的仿真方法。

背景技术

水下光学成像技术利用对海水穿透能力很强的蓝、绿激光代替声波进行通信或探测活动。水下成像系统模型如下：发射光子—光子传输—目标反射—光子传输—接收光子。水下光学成像具有探测目标直观、成像分辨率高、信息含量高等优点，其作为一种新兴的技术与声学成像手段相互补充，在水下成像探测中发挥着重要的作用。

但是由于水介质及水下散射粒子如泥沙、藻类等悬浮颗粒的尺寸与蓝、绿激光波长近似，造成光子发生强烈散射，最终导致水下成像质量降低。目前的水下光学成像质量优化技术包括图像增强和图像复原；前者作用范围有限，不能很好的解决因强烈散射造成的图像雾化现象；后者为基于物理模型的方法，需要足够的假设条件或先验知识作为基础。

Ishimaru等人在《Polarized pulse waves in random discrete scatterers》论文中利用矢量辐射传输方程，分析了聚苯乙烯球形散射介质中偏振度与光学厚度的关系。Laan等人在《Detection range enhancement using circularly polarized light inscattering environments for infrared wavelengths》论文中采用线偏振光和圆偏振光作为入射光源，研究了红外波段偏振光在霾粒子和撒哈拉沙尘粒子中的传输特性。Zhang等人在《Effect of underwater suspended particles on the transmissioncharacteristics of polarized lasers》论文中基于蒙特卡洛模型计算了不同复折射率下，透射距离变化对线偏振光和圆偏振光的前向散射偏振度的影响。

可见，目前对水下光学成像物理模型的建立还停留在对光源在传输介质中的前向散射以及后向散射的特性分析步骤，并未考虑目标反射光的散射和衰减这一因素，即现有的水下光子传输模型无法为基于物理模型的图像复原方法提供足够的假设条件或先验知识。因此目前多为基于特定环境下通过实验的方式获得先验知识并对图像进行复原，存在局限性。

因此，有必要提供一种用于水下光子散射及反射特性模拟分析的仿真方法，解决对不同环境先验知识不足的问题。

发明内容

本发明的目的在于解决目前多为基于特定环境下通过实验的方式获得先验知识并对图像进行复原，存在局限性，对不同环境先验知识不足的问题，而提供了一种基于双向反射理论的水下目标反射偏振光的仿真方法。

为实现上述目的，本发明所提供的技术解决方案是：

基于双向反射理论的水下目标反射偏振光的仿真方法，其特殊之处在于，包括以下步骤：

1)设置初始光源、散射环境以及反射目标的参数；

2)根据初始光源参数，采样初始光子坐标、运动矢量和斯托克斯矢量；根据散射环境参数，计算衰减系数、散射系数和单次反照率；

3)模拟光子移动，根据衰减系数，采样光子移动步长和能量权值；

4)根据光子移动步长和能量权值，判断光子行为：

若光子根据移动步长移动后到达反射面，则为反射，继续判断反射后的光子行为；

若光子根据移动步长移动后到达接收面，则为接收，进入步骤5)；

若光子根据移动步长移动后未到达接受面或反射面，则为散射，返回步骤3)；

若能量权值小于阈值，则消亡；

5)根据接收到的光子的运动矢量，计算最终接收光子的斯托克斯矢量；

6)对最终接收到的光子的坐标、运动方向、能量权值以及斯托克斯矢量进行统计分析，获得设定环境参数下的光子传输特性以及偏振特性的分析结果。

进一步地，步骤1)中，参数包括：发射总光子数N、发射光源频率λ、激光束发散半角φ_div和波束宽度w_width、发射光子斯托克斯矢量S₀＝[I₀ Q₀ U₀ V₀]、水下环境的散射粒子单位体积浓度N₀、粒子半径r、复折射率m＝m1-m2·i、介质的折射率n₁、传输距离d、反射面粗糙度σ、反射材料的折射率n₂、反射材料的吸收系数k_re、反射材料的反射系数r_re以及反射材料的尺寸。

进一步地，步骤2)具体是：

定义光子位于直角坐标系的原点(0，0，0)，z轴向右的正方向为光子传输的方向；

由步骤1)设置的激光束发散半角φ_div和波束宽度w_width，得到波束的焦距为f_l＝-w_width/φ_div；

光子的初始散射角θ₀根据激光束的初始条件而定，θ₀＝-r₀/f_l，其中

为波束半径的采样值，光子的初始方位角为φ₀＝2π·rnd2，其中rnd1，rnd2为0～1均匀分布的随机数，因此在发射平面上光子的初始横坐标和纵坐标分别为：

则光子的初始传播矢量余弦表达式为：

将全局参考系设置为XOZ平面，已知初始斯托克斯矢量为S₀＝[I₀ Q₀ U₀ V₀]，根据光子初始传播矢量，将光子初始斯托克斯矢量S₀旋转至相应子午面参考系，即L(ψ)·S₀，其中L(ψ)为旋转矩阵：

旋转角ψ通过下式解得：

光子在水中移动，碰到水中悬浮粒子，发生散射，基于步骤1)中设置的水下环境的散射粒子单位体积浓度N₀、粒子半径r、复折射率m＝m1-m2·i，根据米氏散射理论，计算出相应的衰减系数c，散射系数b，以及单次反照率ω＝b/c，并根据单次反照率设置能量权值阈值w_T为：

进一步地，步骤3)具体是：

光子在各向异性的介质中移动的步长由光子自由程的概率分布决定，即：P_l(l)＝1-e^-l；

式中：P_l(l)表示光子移动步长为l的概率；

根据光子步长的定义可得l＝cr，c为介质中的衰减系数，r_d为散射前后光子位置的几何距离，则：r_d＝-ln[1-P_l(l)]/c；其中，P_l(l)取0～1之间均匀分布的随机数；

新的光子坐标为(x′,y′,z′)，其中：

假设所有光子的初始权值为w₀＝1，则光子每次发生散射后的能量权值为w_n＝w_n-1·ω。

进一步地，步骤4)中，判断光子行为是散射时，根据散射环境参数采样光子散射角、旋转角、散射后运动方向以及散射后斯托克斯矢量；具体是：

光子到达散射节点与水下悬浮粒子发生碰撞后发生散射，需对散射角θ以及旋转角β进行采样；

散射角θ根据散射相位函数进行采样；使用散射相函数p(θ)描述光散射时的能量分布；由步骤1)可知水下散射环境为单一色散系统，即均匀粒径的粒子系统，散射相函数表示为：

其中，式中S₁和S₂均为散射角θ的散射函数，分别表示散射光复振幅函数的垂直分量和平行分量；

α＝2πr/λ表示悬浮粒子的尺寸参数，λ为发射光源频率，r为粒子半径；

σ_sca为散射截面，由米氏散射理论计算获得；

散射相函数p(θ)满足归一化条件

记散射角为θ_s，则

R_θ为0～1之间均匀分布的随机数，通过数值方法得到θ_s，作为下一次散射的散射角；

下一次散射的旋转角β通过拒绝法进行采样；

根据米氏散射理论可知，散射后的斯托克斯矢量满足以下公式：

I＝m₁₁(θ)I₀+m₁₂(θ)(Q₀cos(2β)+U₀sin(2β)) (6)

根据三角函数公式可得：

记旋转角β_s＝2πR_β，R_β为0～1之间均匀分布的随机数；

将上述采样散射角θ_s和旋转角β_s代入式(6)、(7)，比较I和ξ·I_max值的大小，其中ξ为0～1之间均匀分布的随机数；

若I＞ξ·I_max，则旋转角β_s被采用作为下一次散射的旋转角；

若I≤ξ·I_max，则旋转角β_s被拒绝，重新生成新的旋转角β_s并重复上述步骤；

在对散射角θ_s和旋转角β_s采样完成后，计算出光子新的运动方向：

当|μ_z|≈1时，公式(8)可简化为

光子散射后，更新其斯托克斯矢量：

S_new＝L(-γ)M(θ_S)L(β_S)S_old (10)

其中，L为步骤2)中所述旋转矩阵，M为基于米氏散射理论所得的散射矩阵，γ为将斯托克斯矢量旋转至基准子午面的旋转角：

当π＜β_s＜2π时，取正号；0＜β_s＜π时，取负号；

返回步骤3)。

进一步地，步骤4)中，判断光子行为是反射时，根据反射目标参数，采样光子反射类型、反射后运动矢量、能量权值以及斯托克斯矢量；具体是：

当光子传播至反射平面，即z′＝z+lμ_z≥d时，计算光子沿原运动方向传播至反射平面后的剩余步长，并计算光子反射后新的运动方向；

当光子到达反射目标表面，首先转换光子相对于反射面的入射天顶角θ_i和方位角

满足以下几何关系：

其次是随机抽样微面元法线向量：

①对法向量n＝(0,0,-1)进行扰动，定义微元面法向量为n′，生成两个相互独立的服从均匀随机分布的随机变量u₁和u₂，令：

若

则u₁和u₂被拒绝，重复上述步骤；反之则接受；

②将被接受的U₁和U₂转换为高斯变量：

若

则拒绝，从①开始重新生成随机变量u₁和u₂；反之则接受，找到扰动法线的第三个坐标：

从而获得微元面法线n′＝(n′_x,n′_y,n′_z)；

最后，确定反射类型以及入射天顶角θ_i和方位角

出射天顶角θ_r和方位角

令μ＝(μ_x,μ_y,μ_z)，计算向量μ和n′的点积；

若μ·n′＜0，为镜面反射，反之则为漫反射；

当反射类型为镜面反射时，(μ′_x,μ′_y,μ′_z)＝μ-2n′(μ·n′)，从而出射天顶角θ_r和方位角

当反射类型为漫反射时，生成一个相对于全局法线的随机方向的漫反射；

完成对反射类型及反射后的运动方向的采样后，计算反射后的斯托克斯矢量，并对反射后光子的能量权值进行计算；

若反射类型为镜面反射，根据微面元理论，得到PBRDF的镜面反射部分为：

式中，θ_N为微面元法向与目标表面法向之间的夹角，

为遮蔽因子，表示不同粗糙度表面的遮蔽效应，典型的遮蔽模型认为相邻面元之间的结构呈V形，公式如下：

为4×4的穆勒矩阵；

由于PBRDF模型所描述的几何关系建立于宏观坐标系，而整体的粗糙表面不满足菲涅尔反射公式，因此需要在宏观坐标系和微面元坐标系之间建立其联系；

θ_ni为入射光线在微面元上的入射角，对于微面元模型，其各角度满足

而θ_N由β计算得到：

cosθ_N＝cosθ_i+cosθ_r/(2cosβ) (19)

η_i表示入射方向分别与粗糙表面法线和微面元法线组成的平面的夹角，η_r表示出射方向分别与粗糙表面法线和微面元法线组成的平面的夹角，满足关系：

通过琼斯矩阵表示入射光和散射光s，p振动方向上的电场矢量传输关系：

式中，r_s，r_p为菲涅尔反射系数，其数值大小与介质的折射率n₁、粗糙表面的折射率n₂和θ_ni有关；公式如下：

首先求得微面元上折射角

从而得到

从而得到：

综上，在微面元坐标系中利用菲涅尔反射公式计算微面元的反射后，通过琼斯矩阵与穆勒矩阵的转换关系，可建立微面元坐标系和宏观坐标系之间的联系，从而得到偏振PBRDF的镜面反射部分；

最后，将光子斯托克斯矢量转至全局参考系，与镜面反射穆勒矩阵相乘后，转至反射后的子午面参考系即可；

光子能量权值为w_n＝w_n-1·G；

若为漫反射，PBRDF的漫反射部分为：

式中，ρ_d为漫反射和表面完全Lambertian表面散射的光强比，表示为

式中，R_i代表光线从水中入射到物体表面的菲涅尔反射率，可以由式R_i＝(r_s+r_p)/2计算得到；

K₂是材料表面漫反射光的菲涅耳反射，是表面粗糙度、折射率和散射角的函数；K₂等于R_i或被视为常数；

r_re为当物体表面厚度看为无穷时的反射率系数，对于固定入射波长下给定的材料，为一个常数，由步骤1)可知；

G_d表示物体表面朗伯体漫散射强度的角度分布，G_d＝cosθ_r/π；

穆勒矩阵

描述的是漫散射的退偏过程，表示为

同样将光子斯托克斯矢量转至全局参考系，与镜面反射穆勒矩阵相乘后，转至反射后的子午面参考系即可；

光子能量权值为w_n＝w_n-1·G_d·ρ_d；

返回步骤4)。

进一步地，步骤5)具体是：

将所有成功抵达接收平面的光子的斯托克斯矢量转至全局参考系，与旋转矩阵L(Ω)相乘，Ω＝-arctan(μ_y/μ_x)。

进一步地，所述激光束为高斯激光束。

本发明的优点是：

1.本发明提供的基于双向反射理论的水下目标反射偏振光的仿真方法，模拟光子在水下环境的传输、散射及反射的过程，进而得到接收到的光子的空间坐标、能量、偏振状态，为水下成像、识别等领域的分析研究提供重要且便捷的仿真工具。

2.本发明在设定环境参数以及反射目标参数的情况下可通过模拟仿真的形式对光子的传输特性以及偏振特性进行分析，在减少实验成本的同时保证了光学物理规律的严谨性。

附图说明

图1为实施例1《烟雾浓度对偏振光传输特性的影响》论文中的实验原理图；

图1中的附图标号如下：1-激光器，2-衰减片，3-偏振片，4-1/4波片，5-烟雾模拟环境，6-偏振态测量仪；

图2为实施例1《烟雾浓度对偏振光传输特性的影响》论文中的实验结果图；

图3为本发明仿真方法的流程示意图；

图4为本发明仿真方法的仿真模型图；

图5为实施例1针对《烟雾浓度对偏振光传输特性的影响》论文中的实验环境，采用本发明仿真方法进行模拟实验的仿真结果图；

图6为本发明实施例2的仿真结果图。

具体实施方式

以下结合附图和具体实施例对本发明的内容作进一步的详细描述：

实施例1有效性试验

为验证本发明仿真模型的可靠性，本发明对张肃等人的《烟雾浓度对偏振光传输特性的影响》论文中的实验环境进行模拟，并对比仿真与实验结果。论文中，实验原理如图1所示，激光器发出的光通过衰减片对光束进行衰减，之后经过偏振片进行起偏，通过调节偏振片角度可产生不同的线偏振光，1/4波片用来产生圆偏振光，可根据需要随时旋入或旋出光路，起偏后的偏振光通过烟雾模拟环境，在烟雾模拟环境中燃烧烟煤粒子，并通过底部风扇系统的搅拌使烟雾环境中粒子浓度保持均匀，在出射端由偏振态测量仪接收出射光，并探测其偏振态。

实验通过控制每次燃烧烟煤粒子的量来控制烟雾环境中介质的浓度，从而实现控制散射环境的光学厚度，实验结果如图2所示。

为验证本发明仿真模型的可靠性，针对上述论文的实验环境，应用本发明方法的仿真模型进行模拟对比实验(流程示意图以及仿真模型图参见图3和图4)。仿真中，选择可见光波长532nm，光子个数10⁶个，烟煤气溶胶粒子半径0.37μm，复折射率为1.75+0.44i。光子从原点出发，经散射后到达接收平面，调节与实验对应的光学厚度以便于对比。对不同光学厚度下的不同偏振光进行仿真，对比线偏振光和圆偏振光的前向散射的偏振特性的影响，仿真结果如图5所示。

对比图2和图5可以看出，随着光学厚度的增加，DOP均减小。同时，圆偏振光的偏振性能优于线偏振光，因为圆偏振光的光学厚度值较高。图5中的仿真结果与图2中的实验结果基本一致，因此可以证明，本发明的模型的仿真结果与DOP变化趋势中实际测量值的变化趋势是一致的，从而证明采用本发明仿真方法仿真模型的有效性。

实施例2

为分析相同散射环境下，使用线偏振光对目标进行偏振差分成像，随着光学厚度的变化，后向散射光和目标反射光的偏振度区分情况的变化，使用本发明的仿真模型进行试验。发射总光子数10⁶个、发射光源频率532nm、高斯激光束发射半角0.0015/2和波束宽度0.001、发射光子斯托克斯矢量为线偏振光S₀＝[1 1 0 0]、水下环境的散射粒子单位体积个数为10⁹、粒子半径2μm、复折射率m＝1.52-0.001i、光学厚度为1、反射面粗糙度0.2、反射材料折射率1.46、反射材料的吸收系数1.32反射材料的反射系数0.052。仿真结果如图6所示。

仿真结果可以看出，随着光学厚度的增加，经水下悬浮粒子散射的后向散射光的偏振度逐渐降低，同时，目标反射光的偏振度也因到达其反射表面的光子偏振度整体降低的原因，出现小幅降低的现象。同时，基于水下悬浮粒子散射的后向散射光的偏振度较高，目标反射光较低，这一差异为基础进行的偏振差分成像方法，在随着光学厚度的增加，效果将大大缩减。通过本发明的仿真模型进行试验获得的仿真结果可以作为分析偏振差分算法在普通光学厚度下的可实施性的参考依据。

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

Claims (8)

1.基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于，包括以下步骤：

1)设置初始光源、散射环境以及反射目标的参数；

2)根据初始光源参数，采样初始光子坐标、运动矢量和斯托克斯矢量；根据散射环境参数，计算衰减系数、散射系数和单次反照率；

3)模拟光子移动，根据衰减系数，采样光子移动步长和能量权值；

4)根据光子移动步长和能量权值，判断光子行为：

若光子根据移动步长移动后到达反射面，则为反射，继续判断反射后的光子行为；

若光子根据移动步长移动后到达接收面，则为接收，进入步骤5)；

若光子根据移动步长移动后未到达接受面或反射面，则为散射，返回步骤3)；

若能量权值小于阈值，则消亡；

5)根据接收到的光子的运动矢量，计算最终接收光子的斯托克斯矢量；

6)对最终接收到的光子的坐标、运动方向、能量权值以及斯托克斯矢量进行统计分析，获得设定环境参数下的光子传输特性以及偏振特性的分析结果。

2.根据权利要求1所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于：

步骤1)中，参数包括：发射总光子数N、发射光源频率λ、激光束发散半角φ_div和波束宽度w_width、发射光子斯托克斯矢量S₀＝[I₀ Q₀ U₀ V₀]、水下环境的散射粒子单位体积浓度N₀、粒子半径r、复折射率m＝m1-m2·i、介质的折射率n₁、传输距离d、反射面粗糙度σ、反射材料的折射率n₂、反射材料的吸收系数k_re、反射材料的反射系数r_re以及反射材料的尺寸。

3.根据权利要求2所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于，步骤2)具体是：

定义光子位于直角坐标系的原点(0，0，0)，z轴向右的正方向为光子传输的方向；

由步骤1)设置的激光束发散半角φ_div和波束宽度w_width，得到波束的焦距为f_l＝-w_width/φ_div；

光子的初始散射角θ₀根据激光束的初始条件而定，θ₀＝-r₀/f_l，其中

为波束半径的采样值，光子的初始方位角为φ₀＝2π·rnd2，其中rnd1，rnd2为0～1均匀分布的随机数，因此在发射平面上光子的初始横坐标和纵坐标分别为：

则光子的初始传播矢量余弦表达式为：

将全局参考系设置为XOZ平面，已知初始斯托克斯矢量为S₀＝[I₀ Q₀ U₀ V₀]，根据光子初始传播矢量，将光子初始斯托克斯矢量S₀旋转至相应子午面参考系，即L(ψ)·S₀，其中L(ψ)为旋转矩阵：

旋转角ψ通过下式解得：

光子在水中移动，碰到水中悬浮粒子，发生散射，基于步骤1)中设置的水下环境的散射粒子单位体积浓度N₀、粒子半径r、复折射率m＝m1-m2·i，根据米氏散射理论，计算出相应的衰减系数c，散射系数b，以及单次反照率ω＝b/c，并根据单次反照率设置能量权值阈值w_T为：

4.根据权利要求3所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于，步骤3)具体是：

光子在各向异性的介质中移动的步长由光子自由程的概率分布决定，即：P_l(l)＝1-e^-l；

式中：P_l(l)表示光子移动步长为l的概率；

根据光子步长的定义可得l＝cr，c为介质中的衰减系数，r_d为散射前后光子位置的几何距离，则：r_d＝-ln[1-P_l(l)]/c；其中，P_l(l)取0～1之间均匀分布的随机数；

新的光子坐标为(x′,y′,z′)，其中：

假设所有光子的初始权值为w₀＝1，则光子每次发生散射后的能量权值为w_n＝w_n-1·ω。

5.根据权利要求4所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于：

步骤4)中，判断光子行为是散射时，根据散射环境参数采样光子散射角、旋转角、散射后运动方向以及散射后斯托克斯矢量；具体是：

光子到达散射节点与水下悬浮粒子发生碰撞后发生散射，需对散射角θ以及旋转角β进行采样；

散射角θ根据散射相位函数进行采样；使用散射相函数p(θ)描述光散射时的能量分布；由步骤1)可知水下散射环境为单一色散系统，即均匀粒径的粒子系统，散射相函数表示为：

其中，式中S₁和S₂均为散射角θ的散射函数，分别表示散射光复振幅函数的垂直分量和平行分量；

α＝2πr/λ表示悬浮粒子的尺寸参数，λ为发射光源频率，r为粒子半径；

σ_sca为散射截面，由米氏散射理论计算获得；

散射相函数p(θ)满足归一化条件

记散射角为θ_s，则

R_θ为0～1之间均匀分布的随机数，通过数值方法得到θ_s，作为下一次散射的散射角；

下一次散射的旋转角β通过拒绝法进行采样；

根据米氏散射理论可知，散射后的斯托克斯矢量满足以下公式：

I＝m₁₁(θ)I₀+m₁₂(θ)(Q₀cos(2β)+U₀sin(2β)) (6)

根据三角函数公式可得：

记旋转角β_s＝2πR_β，R_β为0～1之间均匀分布的随机数；

将上述采样散射角θ_s和旋转角β_s代入式(6)、(7)，比较I和ξ·I_max值的大小，其中ξ为0～1之间均匀分布的随机数；

若I＞ξ·I_max，则旋转角β_s被采用作为下一次散射的旋转角；

若I≤ξ·I_max，则旋转角β_s被拒绝，重新生成新的旋转角β_s并重复上述步骤；

在对散射角θ_s和旋转角β_s采样完成后，计算出光子新的运动方向：

当|μ_z|≈1时，公式(8)可简化为

光子散射后，更新其斯托克斯矢量：

S_new＝L(-γ)M(θ_S)L(β_S)S_old (10)

其中，L为步骤2)中所述旋转矩阵，M为基于米氏散射理论所得的散射矩阵，γ为将斯托克斯矢量旋转至基准子午面的旋转角：

当π＜β_s＜2π时，取正号；0＜β_s＜π时，取负号；

返回步骤3)。

6.根据权利要求4所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于：

步骤4)中，判断光子行为是反射时，根据反射目标参数，采样光子反射类型、反射后运动矢量、能量权值以及斯托克斯矢量；具体是：

当光子传播至反射平面，即z′＝z+lμ_z≥d时，计算光子沿原运动方向传播至反射平面后的剩余步长，并计算光子反射后新的运动方向；

当光子到达反射目标表面，首先转换光子相对于反射面的入射天顶角θ_i和方位角

满足以下几何关系：

其次是随机抽样微面元法线向量：

①对法向量n＝(0,0,-1)进行扰动，定义微元面法向量为n′，生成两个相互独立的服从均匀随机分布的随机变量u₁和u₂，令：

若

则u₁和u₂被拒绝，重复上述步骤；反之则接受；

②将被接受的U₁和U₂转换为高斯变量：

若

则拒绝，从①开始重新生成随机变量u₁和u₂；反之则接受，找到扰动法线的第三个坐标：

从而获得微元面法线n′＝(n′_x,n′_y,n′_z)；

最后，确定反射类型以及入射天顶角θ_i和方位角

出射天顶角θ_r和方位角

令μ＝(μ_x,μ_y,μ_z)，计算向量μ和n′的点积；

若μ·n′＜0，为镜面反射，反之则为漫反射；

当反射类型为镜面反射时，(μ′_x,μ′_y,μ′_z)＝μ-2n′(μ·n′)，从而出射天顶角θ_r和方位角

当反射类型为漫反射时，生成一个相对于全局法线的随机方向的漫反射；

完成对反射类型及反射后的运动方向的采样后，计算反射后的斯托克斯矢量，并对反射后光子的能量权值进行计算；

若反射类型为镜面反射，根据微面元理论，得到PBRDF的镜面反射部分为：

式中，θ_N为微面元法向与目标表面法向之间的夹角，

为遮蔽因子，表示不同粗糙度表面的遮蔽效应，典型的遮蔽模型认为相邻面元之间的结构呈V形，公式如下：

为4×4的穆勒矩阵；

由于PBRDF模型所描述的几何关系建立于宏观坐标系，而整体的粗糙表面不满足菲涅尔反射公式，因此需要在宏观坐标系和微面元坐标系之间建立其联系；

θ_ni为入射光线在微面元上的入射角，对于微面元模型，其各角度满足

而θ_N由β计算得到：

cosθ_N＝cosθ_i+cosθ_r/(2cosβ) (19)

η_i表示入射方向分别与粗糙表面法线和微面元法线组成的平面的夹角，η_r表示出射方向分别与粗糙表面法线和微面元法线组成的平面的夹角，满足关系：

通过琼斯矩阵表示入射光和散射光s，p振动方向上的电场矢量传输关系：

式中，r_s，r_p为菲涅尔反射系数，其数值大小与介质的折射率n₁、粗糙表面的折射率n₂和θ_ni有关；公式如下：

首先求得微面元上折射角

从而得到

从而得到：

综上，在微面元坐标系中利用菲涅尔反射公式计算微面元的反射后，通过琼斯矩阵与穆勒矩阵的转换关系，可建立微面元坐标系和宏观坐标系之间的联系，从而得到偏振PBRDF的镜面反射部分；

最后，将光子斯托克斯矢量转至全局参考系，与镜面反射穆勒矩阵相乘后，转至反射后的子午面参考系即可；

光子能量权值为w_n＝w_n-1·G；

若为漫反射，PBRDF的漫反射部分为：

式中，ρ_d为漫反射和表面完全Lambertian表面散射的光强比，表示为

式中，R_i代表光线从水中入射到物体表面的菲涅尔反射率，可以由式R_i＝(r_s+r_p)/2计算得到；

K₂是材料表面漫反射光的菲涅耳反射，是表面粗糙度、折射率和散射角的函数；K₂等于R_i或被视为常数；

r_re为当物体表面厚度看为无穷时的反射率系数，对于固定入射波长下给定的材料，为一个常数，由步骤1)可知；

G_d表示物体表面朗伯体漫散射强度的角度分布，G_d＝cosθ_r/π；

穆勒矩阵

描述的是漫散射的退偏过程，表示为

同样将光子斯托克斯矢量转至全局参考系，与镜面反射穆勒矩阵相乘后，转至反射后的子午面参考系即可；

光子能量权值为w_n＝w_n-1·G_d·ρ_d；

返回步骤4)。

7.根据权利要求5或6所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于，步骤5)具体是：

将所有成功抵达接收平面的光子的斯托克斯矢量转至全局参考系，与旋转矩阵L(Ω)相乘，Ω＝-arctan(μ_y/μ_x)。

8.根据权利要求7所述的基于双向反射理论的水下目标反射偏振光的仿真方法，其特征在于：所述激光束为高斯激光束。

Images