1.   引　言

现有导航方式多依赖卫星等有源信号，在水下导航环境中缺少自主性强的导航模块。水下偏振光存在规律性分布模式，可以根据太阳位置进行预测^[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]进行对比，验证该模型的可行性。

2.   水下偏振模式理论与仿真模型

2.1   水下偏振传输模型

太阳光进入大气层后，被大气中的分子和气溶胶散射，产生偏振光，在大气中形成特定的偏振分布模式。大气中的偏振光在水和空气的界面发生折射，形成了水下的偏振分布模式。在偏振光的散射和折射过程中，可利用Stokes矢量和Mueller矩阵来描述偏振光特性与光学过程。

在大气-水体的介质组成中，光束经过大气分子和气溶胶的散射后，被空气-水界面折射，最终形成水下的偏振模式。如图1（彩图见期刊电子版）所示。

图  1  太阳光的水下偏振传输模型^[13]

Figure  1.  Underwater polarized transmission model of sunlight^[13]

下载: 全尺寸图片 幻灯片

2.2   天空光偏振分布

在晴朗的天气条件下，太阳光主要被大气中的氮、氧等气体分子散射，该散射过程被称为Rayleigh散射。光的偏振态可用Stokes矢量描述：

其中，S为光的Stokes矢量，I为光强，Q和U分别代表两个垂直方向的线偏振光，V代表圆偏振光。大气分子散射过程中V = 0。大气中发生的Rayleigh散射对光的偏振状态的影响通常用散射矩阵${{\boldsymbol{M}}_S}$（Mueller矩阵）来描述,散射光的Stokes矢量可以表示为：

其中，${\theta _{\rm{s}}}$表示散射角。${{\boldsymbol{S}}_0}$用来描述入射天空光的偏振态。

实际大气环境中，光子发生多次散射的情况广泛存在，因此为准确模拟实际光学过程，采用蒙特卡洛方法仿真太阳光在大气中的多次散射。蒙特卡洛方法的仿真原理是将光在大气中的传输过程转化为光子随机运动结果的统计平均，如图2（彩图见期刊电子版）所示^[17]。

图  2  Monte Carlo方法模拟的实际大气环境中的光学过程原理图^[17]

Figure  2.  Schematic diagram of the light transmission of actual atmosphere simulated by Monte Carlo method^[17]

下载: 全尺寸图片 幻灯片

入射光经过分子散射后得到新的散射方向，新散射方向的Stokes矢量${\boldsymbol{S}}''$计算方式如下：

式中，${i_1}$和${i_2}$分别表示入射子午面和散射子午面与散射面的夹角，${\boldsymbol{L}}(-{i_1})$和${\boldsymbol{L}}(-{i_2})$均为旋转矩阵，分别表示为：

在此基础上，可以计算出大气各个方向的偏振度（Degree of Polarization, DOP）P''和偏振角（Angleof Polarization, AOP）A''，得到大气偏振模式：

2.3   气-水界面折射

大气偏振光入射水中时，会在水体表面发生折射，平静水面下，大气天球的视场会被压缩到视角为97.2°的锥形区域内，如图3所示，这个区域被称为Snell窗。然而，当水面存在波浪时，Snell窗的边缘会变得参差不齐。

图  3  Snell窗^[5]

Figure  3.  Snell’s window^[5]

下载: 全尺寸图片 幻灯片

当大气偏振光在水面发生折射时，偏振态的变化可以用水面折射的Mueller矩阵来表示：

其中，${\theta _{\rm{i}}}$表示入射角，${\theta _{\rm{r}}}$表示折射角。根据Snell折射定律，${n_1}\sin {\theta _{\rm{i}}} = {n_2}\sin {\theta _{\rm{r}}}$ ，${n_1}$和${n_2}$分别表示空气和水的折射率。天空偏振光经水面折射后的Stokes矢量可以表示为：

2.4   水下偏振模式仿真

水下偏振分布模式会受到水面波浪的影响。采用经典P-M海浪谱，可对水面波浪进行仿真^[18]。基于仿真得到的海面波浪，建立反向追踪模型，模型如图4所示。

图  4  反向追踪仿真方法原理图

Figure  4.  Schematic diagram of the reverse tracing simulation method

下载: 全尺寸图片 幻灯片

具体过程如下：（1）首先建立坐标系，设定水下观测点，并设置水下成像面；（2）选取成像面上的成像点，以观测点为起点，成像点为终点，计算入射光线矢量；（3）计算入射矢量与波浪水面的交点，计算该交点位置的法向量；（4）根据折射定律计算出折射光线矢量，求出大气偏振分布上折射光线矢量对应点的Stokes矢量${\boldsymbol{S}}''$；（5）以海面折射点为中心，反求入射光与出射光，计算对应折射矩阵${{\boldsymbol{M}}_R}$；（6）将${\boldsymbol{S}}''$与折射矩阵${{\boldsymbol{M}}_R}$相乘，得到水下Stokes矢量${\boldsymbol{S}}'''$，记录该点；（7）在靶面上遍历$100 \times 100$个点，重复上述操作，获得水下偏振分布。

3.   水下偏振模式特征分析

3.1   大气偏振模式

图5（彩图见期刊电子版）为不同太阳位置的太阳光在大气中发生单次Rayleigh散射后的偏振度和偏振角分布。可见，DOP分布中，最小偏振度位置随太阳高度的降低而降低。30°，60°，90°的AOP分布中仅有太阳位置处为中性点。

图  5  太阳光在大气中发生单次Rayleigh散射的DOP和AOP分布，太阳天顶角分别为30°,60°,90°

Figure  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°

Figure  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

下载: 全尺寸图片 幻灯片

3.2   平静水面下的偏振模式

图7（彩图见期刊电子版）为太阳光经大气多次散射后，在平静的气-水界面折射下的水下DOP和AOP分布。图像为Snell窗区域内的仿真图像。与图6相比，偏振光经过水面的折射，太阳位置的中性点分裂范围变大，最大DOP带向中间位置略有凹陷，最大DOP值变小。在AOP分布中，偏振角接近0°的区域变大。

图  7  大气多次散射光被平静水面折射到水中的DOP和AOP分布，太阳天顶角分别为30°,60°,90°

Figure  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

下载: 全尺寸图片 幻灯片

3.3   正弦水面波下的偏振模式

图8（彩图见期刊电子版）为大气偏振光经过正弦形状的水面波进入水下后，水下偏振模式的DOP和AOP分布。可以看到，图像周围明显出现了波浪形状的边界，这是由于波浪导致Snell窗出现不规则的边缘，与水下实拍图像中的Snell窗相符合。与图7相比，DOP和AOP分布的主要特点没有明显的改变，但中心位置存在额外凹陷，太阳位置上仍存在中性点分裂现象，DOP和AOP分布呈现明显的对称性。另外，水面的波动影响了DOP和AOP的分布，由于水面波折射的影响，图像出现了和水面波一致的起伏，这种变化在图像边界上尤为明显。

图  8  （a）正弦波形水面波；（b）大气多次散射光被正弦形状的水面折射到水中的DOP和AOP分布，太阳天顶角分别为30°，60°，90°

Figure  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

下载: 全尺寸图片 幻灯片

3.4   随机水面波下的偏振模式

随机波模拟采用了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°

Figure  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

下载: 全尺寸图片 幻灯片

3.5   仿真和实际测量结果

实测实验采用的偏振成像系统由鱼眼镜头(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模拟结果

Figure  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°时平静水面下的实测结果、理想大气模拟结果、加入气溶胶后的模拟结果对比

Figure  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°

Figure  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°

Figure  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仿真结果与实际测量结果一致。

4.   结　论

在太阳光经过大气多次Rayleigh散射的基础上，仿真了大气偏振光经过平静水面和波动水面折射后的DOP和AOP分布，并得到了与波浪条件相符合的Snell窗现象。仿真结果表明：大气分子对太阳光的多次Rayleigh散射形成了大气偏振模式，大气偏振光经过折射进入水中，形成了水下偏振模式。与平静水面下的偏振模式相比，波动水面下的偏振模式分布出现了与水面波形状一致的变化。不同形状水面下的DOP和AOP在主要特征上呈现一致性：DOP和AOP的整体分布始终保持明显的对称性，水面的折射增加了AOP中性点的分裂程度，并且使DOP分布中的最大偏振带出现了凹陷，偏振分布模式随着太阳高度的变化发生了明显的变化。因此，水下偏振模式分布不仅与太阳位置有关，还受到水面波浪较大影响。将仿真结果与实测结果进行比较，验证了该仿真方法的正确性，为提高水下偏振导航在水面波动条件下的环境适应性提供了理论模型基础。