顾敬桥,李高杰,胡鹏伟,钱建强 *
(1.北京航空航天大学 物理学院,北京 102200;2.北京航空航天大学 自动化科学与电气工程学院,北京 100191)
现有导航方式多依赖卫星等有源信号,在水下导航环境中缺少自主性强的导航模块。水下偏振光存在规律性分布模式,可以根据太阳位置进行预测[1]。许多海洋生物对偏振敏感,并能够利用偏振信息实现导航[2]。水下仿生偏振导航技术由此被提出。这种导航方式具有无源自主、抗电磁干扰、误差不积累、不依赖先验数据等优势。但相较于大气环境,水下环境复杂、测量困难,且对偏振信息存在影响。因此,研究水下偏振光分布对水下导航的发展具有重要意义。
水下偏振光的分布主要受到太阳位置、水的光学特性、观测深度和水面波浪情况的影响[3]。太阳光经过大气散射形成天空光,再经气-水界面折射和水分子散射等光学过程,最终在水下形成具有一定规律的偏振分布模式[4]。从水下观察大气时,天球被压缩到一个锥形区域内,这个锥形区域被称为Snell 窗[5]。水下环境的复杂性给水下偏振模式测量带来了难度[6]。Waterman 于1954 年最早研究了水下偏振模式分布特征[7]。Horváth 和Varjú于1995 年利用单次Rayleigh 散射和天空偏振模型,计算了折射对平静水面下Snell 窗内偏振模式的影响[8]。对于偏振光的空气传输过程,已有较多基于大气气溶胶的光学参数确定方法[9]。对于水下光场,何大华等提出了一种求解光场的迭代方法[10]。在水下偏振分布模式方面,国内的褚金奎课题组在近年来进行了系统性的研究,包括基于水分子单次Rayleigh 散射的静态水面下偏振模式[11]和基于Mie 散射的浑浊水下偏振模式[12]。提出了一种基于水下偏振光导航的实时引导方法[13]。分析了太阳位置、波长变化与浊度变化等多种因素对水下偏振模式的影响,特别针对水面波动问题,通过对水面波时间平均模拟得到了波浪水面下偏振的平均趋势,但未能模拟波浪存在时水下Snell 窗的形变现象[14-15]。现有研究中,仿真结果与波浪对水下偏振光分布模式的实际作用效果差距较大。因此,需要建立更为有效的模型,来描述不同形状的波浪对水下偏振模式的影响,使水下偏振模式更趋近于实测结果。对波面下透射光的偏振分布进行理论模拟,可揭示波面下偏振光分布模式的一般规律,为水下导航增加理论依据。
本工作的目的是建立一个简单有效的模型来描述波动水下偏振光分布模式,通过蒙特卡罗法仿真获得多次Rayleigh 散射下的大气偏振分布模式与波浪特征,并利用折射定律对波浪水下偏振光分布模式进行计算,得到不同太阳高度下的水下偏振光分布模式,最后,将仿真结果分别与在北京航空航天大学的实验数据及文献报道的在Santa Barbara 实测实验结果[16]进行对比,验证该模型的可行性。
太阳光进入大气层后,被大气中的分子和气溶胶散射,产生偏振光,在大气中形成特定的偏振分布模式。大气中的偏振光在水和空气的界面发生折射,形成了水下的偏振分布模式。在偏振光的散射和折射过程中,可利用Stokes 矢量和Mueller 矩阵来描述偏振光特性与光学过程。
在大气-水体的介质组成中,光束经过大气分子和气溶胶的散射后,被空气-水界面折射,最终形成水下的偏振模式。如图1(彩图见期刊电子版)所示。
图1 太阳光的水下偏振传输模型[13]Fig.1 Underwater polarized transmission model of sunlight [13]
在晴朗的天气条件下,太阳光主要被大气中的氮、氧等气体分子散射,该散射过程被称为Rayleigh 散射。光的偏振态可用Stokes 矢量描述:
其中,S为光的Stokes 矢量,I为光强,Q和U分别代表两个垂直方向的线偏振光,V代表圆偏振光。大气分子散射过程中V=0。大气中发生的Rayleigh 散射对光的偏振状态的影响通常用散射矩阵MS(Mueller 矩阵)来描述,散射光的Stokes矢量可以表示为:
其中,θs表示散射角。S0用来描述入射天空光的偏振态。
实际大气环境中,光子发生多次散射的情况广泛存在,因此为准确模拟实际光学过程,采用蒙特卡洛方法仿真太阳光在大气中的多次散射。蒙特卡洛方法的仿真原理是将光在大气中的传输过程转化为光子随机运动结果的统计平均,如图2(彩图见期刊电子版)所示[17]。
图2 Monte Carlo 方法模拟的实际大气环境中的光学过程原理图[17]Fig.2 Schematic diagram of the light transmission of actual atmosphere simulated by Monte Carlo method[17]
入射光经过分子散射后得到新的散射方向,新散射方向的Stokes 矢量S′′计算方式如下:
式中,i1和i2分别表示入射子午面和散射子午面与散射面的夹角,L(-i1)和L(-i2)均为旋转矩阵,分别表示为:
在此基础上,可以计算出大气各个方向的偏振 度(Degree of Polarization,DOP)P''和偏振角(Angle of Polarization,AOP)A'',得到大气偏振模式:
大气偏振光入射水中时,会在水体表面发生折射,平静水面下,大气天球的视场会被压缩到视角为97.2°的锥形区域内,如图3 所示,这个区域被称为Snell 窗。然而,当水面存在波浪时,Snell窗的边缘会变得参差不齐。
当大气偏振光在水面发生折射时,偏振态的 变化可以用水面折射的Mueller 矩阵来表示:
其中,θi表示入射角,θr表示折射角。根据Snell 折射定律,n1sinθi=n2sinθr,n1和n2分别表示空气和水的折射率。天空偏振光经水面折射后的Stokes 矢量可以表示为:
水下偏振分布模式会受到水面波浪的影响。采用经典P-M 海浪谱,可对水面波浪进行仿真[18]。基于仿真得到的海面波浪,建立反向追踪模型,模型如图4 所示。
图4 反向追踪仿真方法原理图Fig.4 Schematic diagram of the reverse tracing simulation method
具体过程如下:(1)首先建立坐标系,设定水下观测点,并设置水下成像面;(2)选取成像面上的成像点,以观测点为起点,成像点为终点,计算入射光线矢量;(3)计算入射矢量与波浪水面的交点,计算该交点位置的法向量;(4)根据折射定律计算出折射光线矢量,求出大气偏振分布上折射光线矢量对应点的Stokes矢量S′′;(5)以海面折射点为中心,反求入射光与出射光,计算对应折射矩阵MR;(6)将S′′与折射矩阵MR相乘,得到水下Stokes 矢量S′′′,记录该点;(7)在靶面上遍历100×100个点,重复上述操作,获得水下偏振分布。
图5(彩图见期刊电子版)为不同太阳位置的太阳光在大气中发生单次Rayleigh 散射后的偏振度和偏振角分布。可见,DOP 分布中,最小偏振度位置随太阳高度的降低而降低。30°,60°,90°的AOP 分布中仅有太阳位置处为中性点。
图5 太阳光在大气中发生单次Rayleigh 散射的DOP 和AOP 分布,太阳天顶角分别为30°,60°,90°Fig.5 The DOP and AOP distributions after single Rayleigh scattering of sunlight in the atmosphere.The zenith angles of the sun are 30°,60°,and 90°,respectively
图6(彩图见期刊电子版)为不同太阳位置的太阳光在大气中发生多次Rayleigh 散射后的DOP 和AOP 分布。和图5 相比,大气分子的多次散射导致太阳位置的中性点发生了分裂,DOP 全局范围略有降低,在太阳位置,DOP 接近零的区域变大,DOP 和AOP分布仍然呈现明显的对称性。
图6 太阳光在大气中发生多次Rayleigh 散射的DOP 和AOP 分布,太阳天顶角分别为30°,60°,90°Fig.6 The DOP and AOP distributions after multiple Rayleigh scattering of sunlight in the atmosphere.The zenith angles of the sun are 30°,60°,and 90°,respectively
图7(彩图见期刊电子版)为太阳光经大气多次散射后,在平静的气-水界面折射下的水下DOP 和AOP 分布。图像为Snell 窗区域内的仿真图像。与图6 相比,偏振光经过水面的折射,太阳位置的中性点分裂范围变大,最大DOP 带向中间位置略有凹陷,最大DOP 值变小。在AOP 分布中,偏振角接近0°的区域变大。
图7 大气多次散射光被平静水面折射到水中的DOP 和AOP 分布,太阳天顶角分别为30°,60°,90°Fig.7 The distribution of DOP and AOP after multiple scattering of sunlight and refraction by calm water.The solar zenith angles are 30°,60°,and 90°,respectively
图8(彩图见期刊电子版)为大气偏振光经过正弦形状的水面波进入水下后,水下偏振模式的DOP 和AOP 分布。可以看到,图像周围明显出现了波浪形状的边界,这是由于波浪导致Snell 窗出现不规则的边缘,与水下实拍图像中的Snell 窗相符合。与图7 相比,DOP 和AOP 分布的主要特点没有明显的改变,但中心位置存在额外凹陷,太阳位置上仍存在中性点分裂现象,DOP 和AOP 分布呈现明显的对称性。另外,水面的波动影响了DOP 和AOP 的分布,由于水面波折射的影响,图像出现了和水面波一致的起伏,这种变化在图像边界上尤为明显。
图8 (a)正弦波形水面波;(b)大气多次散射光被正弦形状的水面折射到水中的DOP 和AOP 分布,太阳天顶角分别为30°,60°,90°Fig.8 (a) Sinusoidal surface waves;(b) the underwater distribution of DOP and AOP after multiple scattering of sunlight and refraction by sinusoidal surface waves.The solar zenith angles are 30°,60°,and 90°,respectively
随机波模拟采用了P-M 海浪谱,该模型由Pierson 与Moskowitz 提出。随机波计算设定风速为10 m/s。图9(a)(彩图见期刊电子版)为仿真得到的结果,图9(b)(彩图见期刊电子版)为大气偏振光经过随机水面波进入到水下的DOP 和AOP 分布。DOP 和AOP 分布图像的特点与图7和图8 仍然具有一致性,即图像仍然出现了和水面波一致的起伏,在边界尤为明显。
图9 (a)随机波形水面波;(b)大气多次散射光被随机波形的水面折射到水中的DOP 和AOP 分布,太阳天顶角分别为30°,60°,90°Fig.9 (a) Random surface wave;(b) the underwater distribution of DOP and AOP after multiple scattering of sunlight and refraction by random waveforms surface waves.The solar zenith angles are 30°,60°,and 90°,respectively
实测实验采用的偏振成像系统由鱼眼镜头(FE185C057HA-1)、偏振相机(PHX050S-P)以及带有密封腔的水平平台组成。其中,偏振相机配备了索尼IMX250MZR 芯片,鱼眼镜头视场大小为185°,观测地点位于北京航空航天大学学院路校区(E116.353°,N39.986°),实验时间为2019 年7 月1 日,当日空气质量指数为36,风力方向为东北风1 级,实验于17 时,18 时无风状态下共测得6 组数据,其对应太阳高度角分别为:65.054°,68.054 °,71.029°,75.431°,79.747°,85.185°(0°为垂直方向)。实际海面波浪情况下观测难度较大,囿于实验条件限制,将随机波波浪情况下的仿真结果与Bhandari 等人在Santa Barbara 测量的数据[16]进行对比,如图10(彩图见期刊电子版)所示。
图10 平静水面下的实测与模拟结果对比。(a)平静水面下的实测图片;(b)实测数据的DOP 分布;(c)实测数据的AOP 分布;(d) DOP 模拟结果;(e) AOP 模拟结果Fig.10 Comparison of measured and simulated results under calm water.(a) Measured images under calm water;(b) DOP distribution of measured data;(c) AOP distribution of measured data;(d) results of DOP simulation;(e) results of AOP simulation
图10(a)为无风状态下实测得到的水下图像,10(b),10(c)为对应的DOP 分布和AOP 分布,对应角度的模拟结果如图10(d),10(e)所示。其中,实测图像圆圈内为Snell 窗范围,模拟图像仅模拟Snell 窗内部图像,Snell 窗外图像偏振模式成因较为复杂,为全反射所成水下图像,其偏振特征取决于光照情况、水面波纹情况与水下反射物的光学性质与分布,因而Snell 窗外无明显特征。在Snell 窗范围内,DOP 和AOP 实测图像均具有良好的对称性。其中,DOP 最低区域位于太阳附近,AOP 对称轴方向指向太阳方向,并在太阳附近处出现明显的中性点分裂情况,其特征与仿真结果预测一致。在实际情况中,大气存在灰尘、水汽等气溶胶,这些气溶胶会造成更明显的中性点分裂。
图11(彩图见期刊电子版)从左到右分别为实测图像,理想大气下模拟图像,与加入气溶胶后的模拟图像。其中气溶胶垂直光学厚度为0.2,模拟波长范围为400~750 nm,气溶胶颗粒是平均有效半径为6 μm 的球状水滴。结果表明,气溶胶导致中性点出现明显分裂,降低了DOP 图片的偏振度,并使DOP 图像偏振度最高区域向远离太阳一侧移动。实际大气中气溶胶成分更加复杂,因缺少当日实时气溶胶光学深度(aerosol optical depth,AOD)观测数据,未针对当日气象情况进行模拟。
图11 高度角为65.054°时平静水面下的实测结果、理想大气模拟结果、加入气溶胶后的模拟结果对比Fig.11 Comparison of measured results under calm water,simulation results under ideal atmosphere,and simulation results with adding aerosols when zenith angle is 65.054°
图12(彩图见期刊电子版)和图13(彩图见期刊电子版)分别为仿真与实际测量得到的水下DOP、AOP 分布,其中,实测数据由Bhandari和Voss 等人于Santa Barbara 海洋下方拍摄得到[16]。将仿真得到的偏振度分布中最大偏振度调整到和实测结果一致,分别为28%,45%,55%,65%。太阳天顶角分别为34°,58°,77°,88°。拍摄时水下深度分别为30 m,11 m,6 m,1.5 m;风速分别为2.5 m/s,6 m/s,6 m/s,6 m/s.
图12 随机波水面下DOP 分布及Santa Barbara 实测得到的DOP 分布的对比图[16]。太阳天顶角分别为34°,58°,77°,88°Fig.12 The DOP distribution under the random wave surface compared with the DOP distribution obtained from the actual measurement in Santa Barbara[16].The solar zenith angles are 34°,58°,77°,88°,respectively
图13 随机波水面下AOP 分布及Santa Barbara 实测得到的AOP 分布的对比图[16]。太阳天顶角分别为34°,58°,77°,88°Fig.13 The AOP distribution under the random wave surface compared with the AOP distribution obtained from the actual measurement in Santa Barbara[16].The solar zenith angle is 34°,58°,77°,88°,respectively
结果表明,天顶角较小时,受到水中其他粒子和波动的影响,实际测量得到的DOP 和AOP 分布特点不明显。太阳天顶角较大时,实际测量结果中,DOP 和AOP 的分布特征明显。从DOP 分布可以明显看到随着太阳天顶角的变大,中间的最大偏振度带变得更加明显,并逐渐向中间凹陷。AOP 分布中,天顶位置出现了偏振角的反转。仿真结果中,随着太阳天顶角的增加,最大偏振度带同样变得更加明显,并逐渐向中间凹陷。偏振角分布中,中性点位置出现分裂特征,天顶位置出现偏振角的反转。实际测量中,海浪波面情况随时间不断变化,故模拟可对整体特征进行描述,但无法对相位等实时特征造成的细节影响进行模拟,除上述时变特征外,AOP 和DOP 仿真结果与实际测量结果一致。
在太阳光经过大气多次Rayleigh 散射的基础上,仿真了大气偏振光经过平静水面和波动水面折射后的DOP 和AOP 分布,并得到了与波浪条件相符合的Snell 窗现象。仿真结果表明:大气分子对太阳光的多次Rayleigh 散射形成了大气偏振模式,大气偏振光经过折射进入水中,形成了水下偏振模式。与平静水面下的偏振模式相比,波动水面下的偏振模式分布出现了与水面波形状一致的变化。不同形状水面下的DOP 和AOP 在主要特征上呈现一致性:DOP 和AOP 的整体分布始终保持明显的对称性,水面的折射增加了AOP 中性点的分裂程度,并且使DOP 分布中的最大偏振带出现了凹陷,偏振分布模式随着太阳高度的变化发生了明显的变化。因此,水下偏振模式分布不仅与太阳位置有关,还受到水面波浪较大影 响。将仿真结果与实测结果进行比较,验证了该 仿真方法的正确性,为提高水下偏振导航在 水面波动条件下的环境适应性提供了理论模型 基础。