海流对海水养殖、航行运输、海洋工程和海洋渔业等活动有显著影响, 实时快速、高精度、大面积的海流观测技术是海洋强国战略实施的重要保障。近年来, 利用海流计、声学多普勒流速剖面仪等仪器进行现场观测的传统方法已经不能满足要求, 遥感观测作为一种新的观测技术逐渐受到关注^[1-2]。

常用的遥感观测海流的方式主要有以下几种。星载雷达高度计和合成孔径雷达(Synthetic Aperture Radar, SAR)可以获得大范围的流场, 但是雷达高度计只能获得地转流^[3-4]、SAR只能观测径向流速^[5-6], 不能得到真实的流场, 并且卫星轨道的重访周期长, 不能满足实时观测海洋的需求。岸基X波段雷达^[7]和高频地波雷达^[8]等能够长期观测波浪和流场, 但是岸基观测的覆盖范围较小, 并且雷达站的选址困难、安装使用不方便。

自从Cox和Munk^[9]利用机载光学相机观测海面粗糙度以来, 光学遥感也被用于海浪和海流的观测。除了机载观测方式, 刘海江等^[10]基于视频图像提取了波浪的相速度; 姜文正等^[11]提出了一种无海面控制点的立体摄影海浪测量方法, 并将测量结果与波浪骑士的结果进行对比, 验证了该系统的可靠性, 但是这两种观测方式需要固定位置安装。光学遥感测流常用的方法主要是三维傅里叶变换方法(3-Dimensional Fast Fourier Transformation, 3-D FFT)^[12]、最大相关系数法(Maximum Cross-Correlation, MCC)、粒子图像测速技术(Particle Image Velocimetry, PIV)和激光雷达多普勒测速等。Dugan等^[13]将3-D FFT算法应用于处理机载光学相机拍摄海面图像序列; 将MCC应用于机载红外成像系统上^[14], 获得了海流和水深。Anderson等^[15]开发了机载远程海洋洋流成像系统(Remote Ocean Current Imaging System, ROCIS)拍摄海浪图像, 并利用3-D FFT算法获得海洋表面流。

MCC通过两幅以上遥感图像中特定示踪物的移动来获取海流, 所用的示踪物可以是海面温度(sea surface temperature, SST)^[16]、叶绿素浓度^[17]等海洋水色参数。但是, 影响海表温度变化的因素有很多, 除了平流之外海气交换、垂向对流等也会造成海表面温度的变化, 并且由于海洋流体运动会造成严重的几何变形, 这样利用MCC法求得的位移矢量不一定是海流的流动方向。Wu等^[18]用相关松弛法(Correlation Relaxation, CR)改进了MCC法, 将单个位置处的海流与周围海流联系起来, 保证了局部流场的一致性。毛志华等^[19]比较了MCC和CR法, 认为CR法计算海表流场具有一定的可靠性。总之, MCC法能初步获取海流, 但获取的数据较为粗糙, 示踪物数据的空间分辨率对MCC估计的影响较大, 很难得到海洋锋、中尺度涡附近的流场, 也会受到云层覆盖、不良的卫星观测条件的影响。

PIV技术是通过匹配随流体运动的示踪粒子来测量流速, 可用于实验室环境下各种复杂流动的观测^[20]。Fujita等^[21]在PIV技术的基础上提出了大尺度粒子图像测速技术(Large Scale Particle Image Velocimetry, LSPIV), 利用归一化互相关方法对河流表面瞬时流速进行估计。Scarano等^[22]在PIV方法的基础上引入了迭代多重网格方法, 用于提高全局测量精度。Bechle等^[23]在对运动矢量进行估算时, 利用最小二次差分方法来适应非均匀分布的水面模式, 但计算量较大。张振等^[24]认为基于时空图像的测速方法对输入噪声的鲁棒性更强, 可以降低算法的时间复杂度, 但是这种方法是面向流向的一维测速, 不能测量不同方向的循环流动。受实验现场复杂性的影响, 例如水流示踪物密度低、时空分布不均、低空间分辨率、水面光学噪声大、流场定标困难等, PIV技术在天然河流流速、流量监测中的应用还有较大的不确定性, 得到的瞬时流场矢量正确率较低。

激光雷达因其特殊的工作波段(一般为蓝绿波段), 对海水有较强的穿透能力^[25]。海水中的悬浮粒子可以散射激光雷达发射的入射光, 假设悬浮粒子只受海流的影响, 从而可以根据多普勒原理估算海水的流速, 其优点是动态响应快、可非接触测量和实时性佳等, 但其测量范围和精度受环境和仪器的影响较大, 并且不同海域、不同海水的光学特性有很大差异, 因此激光雷达测速多用于实验室流场分析^[26]。

近年来, 无人机遥感技术不断发展, 通过搭载可见光相机、红外相机、高光谱成像仪等实时提供高分辨率遥感数据。基于无人机平台可以对海洋灾害、海洋污染物等进行监测^[27-28]。在海洋动力要素观测方面, Dugan等^[29-30]多年来在利用小型飞机采集海面图像提取海面浪、流参数方面做了许多卓越的工作, 设计研发了机载远程光学聚焦系统(Airborne Remote Optical Spotlight System, AROSS), 可以将其应用到无人机平台上。Streßer等^[31]利用无人机搭载光学相机垂直下视拍摄河流表面视频图像, 并基于3-D FFT算法从河流表面视频图像中提取出了河流表面流速, 测量结果与ADCP实测数据相比具有较好的一致性。Lund等^[32]利用无人机光学海面视频数据和3-D FFT算法对海表面流场进行测量, 并针对无人机的航向和位置偏移提出了一种基于波浪信号的图像校正方法。Tauro等^[33]利用四旋翼无人机搭载相机拍摄河流表面图像, 基于LSPIV法来获得了较高精度的时间平均水面流速, 但是LSPIV法受水环境的影响较大。Fulton等^[34]将多普勒测速雷达集成在无人机平台上, 根据水面上小尺度表面波的多普勒频移来测量河流流速, 多普勒雷达可全天候测流, 受水质影响较小, 但是当流速较小(海面较为平静)时信噪比急剧降低, 因此无法测量较低的流速^[35]。

与传统方法相比较, 小型无人机具有起降环境要求低、性能可靠、机动灵活等特点, 更加适合于海岛、海岸带遥感调查的需要^[36]。但是, 对小型无人机观测海流的最佳飞行指标、使用环境要素的研究还不成熟。本文研究利用小型无人机搭载光学相机观测海面流场的方法, 并分析各种环境因素对海流估算的影响, 包括无人机和相机参数(时间分辨率、空间分辨率、航速、飞行高度等)、反演算法参数的选取等, 为无人机观测海流技术的应用提供基础。

1 实验数据

本文的实验地点选在山东威海的褚岛附近海域, 如图 1(a)所示为实验区域的水深图, 红框区域为实验地点, 使用的数据来源于全球海洋测深数据集(General Bathymetric Chart of the Oceans, GEBCO), 并利用克里金插值对数据进行处理。褚岛的北面海域开阔, 该处海区的年平均有效波高为0.6~0.8 m, 平均有效波周期为3~6 s, 波高向外海逐渐增大, 且受地形和季风的影响^[37]。

实验中, 无人机搭载的光学相机用于拍摄海面图像, 并利用GPS等设备记录无人机的飞行姿态数据(无人机俯仰角、滚转角、航向角、飞行高度)和相机拍摄参数(相机俯仰角、航向角、拍摄时间、经纬度)。每幅海面图像覆盖的海域面积与无人机的飞行高度、搭载的光学相机的视场角以及相机俯仰角有关。实验过程中使用同一光学相机进行拍摄, 实验中用到的相机固有参数见表 1, 另外, 相邻两幅图像的平均采样时间间隔为3 s, 相机在x、y方向上的视场角分别为32°、23°, 拍摄海面时相机主光轴方向与垂直方向的角度为30°左右。

为了验证无人机观测海流的精度, 实验期间无人机围绕浮标飞行, 该浮标在实验期间观测了海流的流速和流向。图 1(b)分别显示了18日和19日5个飞行架次的航线图, 红色星号为浮标所在位置, 红色虚线为浮标周围半径500 m的海域, 具体的飞行参数及浮标实测海浪和海流数据如表 2所示。

2021年11月18和19日无人机分别从船上和岸上起飞, 利用无人机搭载高分辨率光学相机拍摄光学海面图像序列。原始的光学海面图像成像过程中会受到许多因素的影响, 主要有以下几点: 第一, 相机以倾斜姿态拍摄的图像会发生几何畸变。第二, 太阳与相机视线方向在一定的角度下, 图像上部分区域会出现太阳耀斑, 造成图像亮度分布不均, 如图 2所示是经过几何校正以后的海面灰度值图像, 图中由于太阳耀斑的影响, 不同位置处的图像亮度分布有很大差异。与太阳活动中太阳表面局部区域突然变亮的“耀斑”现象不同, 本文中的太阳耀斑是指当太阳光在海面上发生镜面反射, 并在一定观测角度下被光学传感器捕捉到, 在成像图中形成了非常亮的光斑, 则将图像中的光斑称为太阳耀斑。第三, 当拍摄海面上有浮标、浮球、海鸟、船只、浪花等物体时, 会在图像上显示出不同的光学特征, 而海浪作为背景物与这些异物相比灰度值反差很大, 图像整体会显得较暗。除此以外, 为了获取相同海面区域处的子图像序列, 需要通过飞机的飞行参数和图像拍摄时的经纬度信息等数据来对不同图像进行图像配准。在第3节我们会给出图像预处理的相关过程。

2 方法

为了利用无人机获取的海面光学图像序列反演海流, 首先对光学图像做预处理获得海浪信息, 然后利用三维傅里叶变换和频散关系估算海流的流速和流向。

2.1 图像预处理

为了从无人机拍摄的光学图像中准确提取海浪条纹, 需要根据无人机的飞行姿态数据和相机拍摄参数等对无人机光学图像数据做几何校正、坐标转换、图像插值、图像增强等预处理。另外, 为了减小太阳耀斑对海浪成像的影响, 本节还提出了太阳耀斑的校正方法。

2.1.1 几何校正

相机倾斜拍摄海面时, 海面上各点在像片上的投影会根据其在像片上的位置不同而产生不同程度的像点位移, 因此需要对倾斜拍摄的图像进行几何校正。为了减小海面定标的困难, 本文根据摄像机内部的物理模型, 结合相机拍摄海面目标系统的几何关系对图像进行校正^[38], 具体方法如下。

利用光学相机对海面进行拍摄, 海面各点反射或散射来自太阳和天空的光线, 经由相机物镜投射到相机底片的感光层上, 形成物像。各光线汇聚的焦点S称为摄影中心, 摄影中心点S距离平均海平面的高度H为飞行高度, 摄影机焦距为f ^[39-40]。在无人机相机的镜头与被测海面目标构成的系统(图 3)中, A、B两点分别为在某一方向上(x或y方向)相机拍摄海面时单幅取景框内能够拍摄到的最大范围所在点, a、b两点分别为A、B两点在相机底片上的物像点, o点为a、b两点的中点。当像片水平, 拍摄的海面水平时, 由相似三角形理论可知, 像片上a、b两像点间的距离l与海面上A、B两物点间的距离L有以下关系:

$\frac{l}{L} = \frac{f}{H} .$

(1)

$\overline {ab}$的中点o和$\overline {AB}$的中点O, 两点的连线为主光轴, 直线oO分别平分∠aSb和∠ASB, ∠OSB为视场角α, 根据几何关系, 视场角$\alpha = \arctan \frac{L}{{2H}}$。有x方向视场角和y方向视场角, 若探测器x方向分辨率为N_x, 探测器y方向分辨率为N_y, 像元大小为a, 则x和y方向的视场角分别为:

$\left\{ \begin{gathered} {\alpha _x} = \arctan \frac{{a{N_x}}}{{2f}} \hfill \\ {\alpha _y} = \arctan \frac{{a{N_y}}}{{2f}} \hfill \\ \end{gathered} \right. .$

(2)

相机的主光轴与平均海平面的法线即铅垂光线之间的夹角, 称为像片倾角β。在相机倾斜拍摄海面目标的系统(图 4)中, 摄影中心S在海平面上的投影为S^′, SS^′垂直于平均海平面。M为平均海平面上某一被测目标点, 物点M在像平面上的像点为m。假设m、o两像点间距离为$\overline {om}$, 则其在平均海平面上相应的物点间距离为$\overline {OM}$, $\overline {MP}$为$\overline {OM}$垂直于主光轴方向的正射投影, ∠α为直线Mm与主光轴oO之间的夹角, 由相似三角形理论有$\frac{{\overline {om} }}{{\overline {MP} }} = \frac{f}{{SP}}$, 根据几何关系, $\overline {MP}$可以由以下公式得到:

$\overline {MP} = \frac{H}{{\cos \left( {\beta - \Delta \alpha } \right)}}\cos \Delta \alpha \cdot \frac{{\overline {om} }}{f} .$

(3)

倾斜拍摄时, 相机拍摄到的海面某一方向上的最大距离为L, 由几何关系计算L:

$L = \left( {\frac{H}{{\cos \left| {\alpha - \beta } \right|}} + \frac{H}{{\cos \left| {\alpha + \beta } \right|}}} \right) \cdot \frac{{\sin \alpha }}{{\cos \beta }} .$

(4)

若将不同方向上的总视场角α_x(或α_y)按像元个数划分为N_x(或N_y)个单位视场角, 则在平均海平面上最大拍摄距离L内, 海面上任意一点M_n和摄影中心S的连线SM_n与主光轴oO之间的夹角∠α_n都能依次被表示。那么可以由公式(3)和公式(4)依次得到平均海平面上任意相邻物点间距离$\overline{O M_{n}}$在主光轴方向上的投影大小$\overline{M_{n} P_{n}}$, 并将其作为正射投影下的海面像素大小来对倾斜拍摄的海面图像进行校正, 这样就可以得到几何校正以后的海面图像。经过几何校正以后的海面图像, 不同像点间的距离大小表示的是实际平均海平面上两物点间的距离, 即可以得到海面图像上各点的相对地理坐标。

2.1.2 坐标转换

无人机搭载光学相机对海面波浪进行拍摄, 在连续的多幅海面图像序列中, 相邻的N幅图像上会有一定的区域被重复拍摄, 即存在重叠区域, 从重叠区域选取某个子区, 即可以得到N幅子图序列, 再进行海表面流的反演。为了将连续的海面图像序列联系起来, 本文以拍摄第一幅海面图像时相机在平均海平面上的投影中心点为坐标原点$o$, 相机航向为y轴正方向, 垂直相机航向向右为x轴正方向, 建立图 5中的第一坐标系即$x'oy'$坐标系。图 5中, 短虚线表示的是无人机航向, 点虚线为相机航向, 在实际的无人机飞行实验中, 相机的拍摄方向与无人机的飞行方向并不是完全一致的, 可能会存在一定角度。假设无人机沿相同方向飞行, 若拍摄第一幅海面图像时图像上各点在平均海平面上的相对地理坐标为(x₀, y₀), 那么第i张图像的相对地理坐标(x_i, y_i)可以用以下公式来计算:

$\left\{ \begin{gathered} {x_i} = {x_0} + \left[ {\left( {{d_i} - {d_0}} \right) \pm v \cdot \text{d}t \cdot \cos \left( {{\theta _u}} \right) \cdot i} \right] \cdot \cos \left( {\frac{\pi }{2} - {\theta _u}} \right) \hfill \\ {y_i} = {y_0} + \left[ {\left( {{d_i} - {d_0}} \right) \mp v \cdot \text{d}t \cdot \sin \left( {{\theta _u}} \right) \cdot i} \right] \cdot \sin \left( {\frac{\pi }{2} - {\theta _u}} \right) \hfill \\ \end{gathered} \right. ,$

(5)

其中, (d_i–y_i)为第i张图像的相对地理坐标中心点与第一张图像的相对地理坐标中心点之间的距离, 用每张图像拍摄时的经纬度数据转换到平面直角坐标系下作差来获得。$\pm v \cdot \text{d}t \cdot \cos \left( {{\theta _u}} \right) \cdot i$或$\mp v \cdot \text{d}t \cdot \sin \left( {{\theta _u}} \right) \cdot i$表示在不同方向上对无人机定位的校正项, 与无人机的飞行速度v和相邻图像的采集时间间隔dt有关。在正北方向下, 顺时针旋转为正, θ_u是无人机航向。θ_c是相机航向与正北方向之间的夹角, 当$\left| {{\theta _c}} \right| < \frac{\pi }{2}$时, x和y两个方向上的校正项分别为–$v \cdot \text{d}t \cdot \cos \left( {{\theta _u}} \right) \cdot i$和$v \cdot \text{d}t \cdot \sin \left( {{\theta _u}} \right) \cdot i$; 当$\left| {{\theta _c}} \right| \geqslant \frac{\pi }{2}$时, 两个方向上的校正项分别为$+ v \cdot \text{d}t \cdot \cos \left( {{\theta _u}} \right) \cdot i$和$- v \cdot \text{d}t \cdot \sin \left( {{\theta _u}} \right) \cdot i$。

如果直接对第一坐标系下的海面图像进行海表面流的反演, 那么得到的海流参数是相对于第一坐标系下的, 而浮标测量的海流参数是以浮标为坐标原点, 正北方向为0°, 顺时针旋转为正进行记录的, 本文同样以拍摄第一张图像时相机在平均海平面上的投影点为坐标原点, 正北方向为y轴正方向, 正东方向为x轴正方向建立第二坐标系, 即图 5中的xoy坐标系。

为了能够将利用光学图像数据反演的海流值与浮标实测值进行对比, 需要将反演的海表面流矢量从第一坐标系旋转到第二坐标系下, 利用旋转矩阵R来实现这个过程:

$\left[\begin{array}{l} U_x \\ U_y \end{array}\right]=\boldsymbol{R}\left[\begin{array}{l} U_x^{\prime} \\ U_y^{\prime} \end{array}\right]=\left[\begin{array}{c} \cos \varphi-\sin \varphi \\ \sin \varphi-\cos \varphi \end{array}\right]\left[\begin{array}{l} U_x^{\prime} \\ U_y^{\prime} \end{array}\right],$

(6)

其中, (U_x, U_y)为旋转后的海表面流矢量坐标; $({U'_x}, \;{U'_y})$为旋转前的海表面流矢量坐标; φ = –θ_c为旋转角度, φ > 0时海表面流矢量逆时针旋转, φ < 0时海表面流矢量顺时针旋转。

2.1.3 图像插值

在对海面图像进行几何校正以后, 获得的海面图像数据相邻点之间的距离分布不均匀。例如, 图 6是当飞行高度H=320 m, 像片倾角β=30°时, 经过几何校正以后的海面像素点坐标网格, 从x方向上看, 相邻海面像素点之间的距离在x方向相等; 从y方向上看, 相邻像素点之间的距离沿着y轴正方向不断增加。因此, 本文采用双线性插值将海面图像数据插值到规则网格。

2.1.4 图像增强

为了减少海浪以外的物体在海面图像中的影响, 需要先对原始海面图像进行灰度化处理, 即将真彩色图像RGB转换为灰度强度图像I。当灰度图像中除海浪外有船只、浮标等异物时, 海浪作为背景物与这些异物相比灰度值反差很大, 整个图像会显得很暗; 当存在太阳耀斑时, 图像的亮度分布不均匀。所以需要对灰度图像进行直方图均衡化处理, 提高图像的对比度, 从而抑制噪声, 更加突出了海浪信号。

2.2 太阳耀斑的影响

为了减小太阳耀斑对图像的影响, 本文提出基于灰度值分布的校正方法, 具体如下: 首先, 在较大的海面区域内, 选取一个和反演子区大小相同的“窗”, 每间隔几米分别计算“窗”内各子图像序列的灰度值均值再取平均, 最终得到整个海面区域内的灰度值均值分布; 其次, 对相同纵坐标下的子图像灰度值均值进行二次函数拟合; 最后, 用式(7)对子图像的灰度值进行处理, 再对海表面流进行反演:

$I = {I_o} - {I_{ft}} ,$

(7)

其中, I_o为原子图像的灰度值; I_ft为拟合后的灰度值均值, 用处理过的子图像灰度值$I$计算海表面流。

2.3 海流的反演算法

对无人机拍摄的海面图像预处理后, 可以获得位置校准后的清晰海面图像。根据镜面反射理论, 在拍摄的光学图像上, 海浪在波峰和波谷位置处对应的灰度值大小各不相同, 光学图像的灰度值分布情况与海面粗糙度和海浪的斜率分布有关^[9]。因此, 可以利用3-D FFT方法从海面图像序列中反演出海表面流信息^[12]。

对海面图像序列的灰度值I(x, y, t)做3-D FFT, 将时间-空间域转换为波数-频率域, 得到三维图像谱$F_I^{(3)}({k_x}, \;{k_y}, \;\omega )$:

$F_I^{\left( 3 \right)}\left( {{k_x}, \;{k_y}, \;\omega } \right)\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {I\left( {x, \;y, \;t} \right){e^{ - j\left( {k, \;x + {k_y}y - \omega t} \right)}}\text{d}x\;\text{d}y\;\text{d}t} } } .$

(8)

对应的谱分辨率分别为:

$\Delta \omega = \frac{{2\pi }}{T}, \;\Delta {k_x} = \frac{{2\pi }}{{{L_x}}}, \;\Delta {k_y} = \frac{{2\pi }}{{{L_y}}} ,$

(9)

其中, Δω为角频率分辨率; Δk_x, Δk_y为波数分辨率; L_x, L_y分别为光学海面图像中选取的反演区域在水平方向或垂直方向上的长度, 是反演子图在某方向上的像素个数N和在该方向上的图像分辨率Δx(或Δy)的乘积, 即满足L_x=N·Δx或L_y=N·Δy; T为时间序列的总长度, 是图像幅数n和时间分辨率Δt的乘积, 满足T=n·Δt。图像功率谱E(k_x, k_y, ω)由下面的公式来定义:

$E({k_x}, \;{k_y}, \;\omega ) = \frac{{{{\left| {F_I^{(3)}({k_x}, \;{k_y}, \;\omega )} \right|}^2}}}{{{L_x}{L_y}T}} .$

(10)

结合线性波频散关系, 并利用最小二乘法计算海表面流的公式如下:

$\left[ \begin{gathered} {U_x} \hfill \\ {U_y} \hfill \\ \end{gathered} \right] = {\left[ \begin{gathered} \sum {Ek_x^2} \;\;\;\;\;\sum {E{k_x}{k_y}} \hfill \\ \sum {E{k_x}{k_y}\;\;\sum {Ek_y^2} } \hfill \\ \end{gathered} \right]^{ - 1}}\left[ \begin{gathered} \sum {E{k_x}} (\omega - \sqrt {gk\tanh (kd)} ) \hfill \\ \sum {E{k_y}\;(\omega - \sqrt {gk\tanh (kd)} )} \hfill \\ \end{gathered} \right] .$

(11)

综上所述, 本文提出的利用无人机拍摄的光学海面图像序列反演海表面流场的方法如图 7所示。

3 实验结果与分析

为了验证本文提出的海流观测方法, 首先结合实例分析了图像预处理的影响, 然后利用预处理后的图像反演了流速, 并与浮标实测的流速比较, 分析了不同因素对海流观测的影响。

3.1 无人机光学海面图像反演海表面流 3.1.1 海面图像数据的预处理

以实验F4的数据为例来对海面图像数据的预处理过程进行介绍。实验F4中的图像数据受太阳耀斑影响明显, 图像亮度分布不均, 根据上文3.1节的预处理流程依次对其进行几何校正、坐标转换、图像插值、图像增强等处理。图 8(a)为经几何校正以后的海面灰度图像, 此时图像上各点的坐标为该单幅图像拍摄海面区域内各点的相对地理坐标; 如果以此航次无人机的起点位置为坐标原点$o$, 那么该图像相对于起点位置的相对地理坐标如图 8(b)所示; 选取图 8(b)中的红框区域为一个反演子区, 将其插值到512×512像素的规则网格上, 得到图 8(c), 图像分辨率为0.3 m, 图像尺寸为153.6 m×153.6 m; 由于海面图像整体亮度分布不均, 海浪特征不明显, 我们将子图像做直方图均衡化处理以后得到的图像如图 8(d)所示, 海浪信号被加强, 更有利于下一步从海面图像序列中提取海表面流信息。

3.1.2 海表面流的反演

根据上文2.3节的海表面流反演流程(图 7)依次对连续4幅海面图像进行预处理后, 接着选取4幅海面图像重叠区域内的某个反演区域, 得到子图像的时间序列, 对它们做3-D FFT并计算图像功率谱E(k_x, k_y, ω), 再由公式(8)—(11)计算子图像区域内的海表面流。图 9为某角频率下的图像功率谱, 图中叠加了在该角频率下满足的无海表面流时的频散关系曲线(红色虚线圆圈), 从图中可以很清楚地看出, 海表面流的存在使得频谱能量在波数平面上产生多普勒频移。根据图中频谱位置的改变, 采用最小二乘法, 计算出来的海表面流流速大小U = 0.56 m/s, 与浮标实测流速大小U_t = 0.58 m/s较为接近, 相对误差为–2.94%; 计算出来的海表面流流向相对于正北方向为136°, 而在图像拍摄时间1 h范围内, 浮标实测流向在94°~140°变化。因此, 可以从无人机光学海面图像中提取出合理的海表面流信息。

3.2 不同因素对海表面流反演的影响

为了研究不同因素对无人机观测海流的影响, 本节分析无人机的飞行参数和算法中的参数选取等对海流反演结果的影响。

3.2.1 图像分辨率

从实验F2中选取距离浮标位置接近的无人机光学海面图像, 用相同尺寸、不同分辨率的子图像序列计算海表面流。为了方便计算, 在进行海流反演时, 不同方向上(x或y方向)所取的图像尺寸L_x或L_y相等, 用L来表示; 不同方向上的图像分辨率Δx或Δy也相等, 用dD来表示。图 10(a)为相同图像尺寸L=38.4 m×38.4 m, 不同图像分辨率dD时的反演流速与实测流速之间的均方根误差RMSE(U), 从图中可以看出, 当图像分辨率dD < 0.15 m或dD > 0.15 m时, 反演流速精度较差; 在图像分辨率dD = 0.15 m时均方根误差达到最小, RMSE(U) > 0.15 m。由公式(10)可知, 不同方向上的波数分辨率$\Delta {k_x}$或$\Delta {k_y}$与反演子区的长度L_x或L_y有关, 而L = L_x = L_y = N·dD。那么当图像分辨率dD增大, 像素个数N增加, 则子图上单个波长内的采样点数量增加, 反演效果比较好, 但如果继续提高图像分辨率, 没有体现出更好的反演海流效果。因此, 本文后续研究使用的图像分辨率dD为0.15 m。

3.2.2 图像尺寸

图像尺寸大小的选取与反演海区波浪的主波波长λ有关, 本文研究发现图像尺寸大小L在波浪主波波长λ的2~3倍反演效果较好。例如, 实验F2中浮标实测波浪的主波波长λ在15~20 m, 图 10(b)为相同分辨率、不同图像尺寸下的反演流速与实测流速之间的均方根误差RMSE(U), 可以看出, 当图像尺寸L < 20 m×20 m时, 反演流速精度较差, 而图像尺寸L = 38.4 m×38.4 m在主波波长的2~3倍, 反演流速的均方根误差较小。当图像尺寸L小于或接近λ时无法得到有效的海流反演结果, L过大可能会引入由于非均匀波浪场造成的误差。

3.2.3 图像幅数

在图像分辨率和图像尺寸合适的前提下, 图像幅数$n$也会影响海表面流的反演。从实验F1数据中挑选出连续的4、8、16幅海面图像, 分别计算距离浮标位置较近的子区域海表面流, 得到不同图像幅数反演的流速大小如表 3所示。浮标在x和y方向上的实测流速分别为: U_tx = 0.55 m/s和U_ty = –0.01 m/s。从表 3可以看出, 随着图像幅数的增加, 反演结果与浮标实测结果之间的相对误差逐渐降低, 这是因为增加图像幅数提高了角频率分辨率, 更有利于获得更高精度的流速值。

3.2.4 时间分辨率

为了研究时间分辨率对海流反演的影响, 本文从实验F1数据中选取不同时间分辨率dt下(即Δt)的海面图像序列, dt范围为2.3~8 s, 分别从这些图像序列中计算海表面流流速大小, 并与浮标实测流速进行对比。图 11为反演流速与浮标实测流速之间的相对误差大小, 从结果来看, 相邻两张图片的平均采样时间间隔在2.3~3.4 s反演结果较好, 平均采样时间间隔大于3.4 s以后, 反演结果较差。分析其原因, 根据采样定理, 如果图像的采样时间间隔大于主波周期, 会降低海浪谱反演的精度。从表 2可知, 实验F1中波浪的主波周期为3.8 s, 而实验结果中平均采样时间间隔小于3.8 s的图像数据反演海流效果较好, 因此, 图像的采样时间间隔应满足小于波浪的主波周期。但是, 在观测前一般无法获得波浪的主波周期、主波波长等参数, 所以可以在设置相机拍摄参数时适当提高时间分辨率, 以保证能够获得较好的海流反演结果。

3.2.5 飞行速度和飞行高度

无人机的飞行速度决定了航程, 飞行高度决定了单幅光学海面图像所覆盖的区域大小, 而不同图像间的重叠度与无人机的飞行速度和高度以及相机的采样时间间隔有关。图像重叠度越高, 获取海面上相同区域的图像序列数量越多。

本研究分别采用了不同的飞行速度和飞行高度来进行实验, 表 4给出了不同飞行高度和飞行速度下反演流速与浮标实测流速之间的相对误差。从表中可以看出: 第一, 飞行高度越高, 单幅海面图像所覆盖的区域范围越大。第二, 与反演图像尺寸的选取相类似, 飞行高度的选取也与反演海区波浪的主波波长有关, 实验F1、F2中波浪主波波长λ=15~20 m, 无人机的飞行高度较低时反演的流速精度较高; 实验F3、F4和F5中波浪主波波长λ=40~50 m, 无人机的飞行高度增高时流速反演情况较好, 即在本实验的条件下, 无人机飞行高度大致在反演图像尺寸的3~4倍左右时比较合适。第三, 本研究分别用4、8、16幅海面图像序列均获得了有效的海表面流流速, 我们认为反演海表面流至少需要4幅图像序列, 因此在无人机拍摄图像的过程中, 应满足至少海面上同一反演区域在连续的4幅图像序列中都被拍摄到。另外, 随着图像幅数和时间分辨率的增加, 反演流速精度有所提高, 可以通过降低无人机飞行速度或提高图像采样频率的方式来获取更多的海面图像序列, 但计算量也会随之增大。

3.2.6 太阳耀斑

若无人机搭载的光学相机与太阳之间的相对位置发生改变, 拍摄的光学海面图像会受到不同程度太阳光的影响。在研究过程中, 我们发现太阳耀斑也会对海流的反演造成影响。为研究太阳耀斑对海流反演的影响, 本文选择不同相机航向角、不同高度下拍摄的海面图像, 按照受太阳耀斑影响程度的强弱对其依次进行流场的反演。图 12(a)、(c)、(e)依次是从受太阳耀斑影响越来越大的海面图像中挑选的反演区域, 对其进行图像预处理和图像增强处理后计算海表面流, 图 12(b)、(d)、(f)分别是它们反演的流场图。从图中可以看出, 随着海面图像受太阳耀斑影响的程度增加, 反演出来的流场整体流速越来越小。

4 讨论 4.1 图像定位误差

由于相机搭载在飞行器平台上, 会受到无人机的飞行动作(爬升、下降、悬停和飞行)、飞行姿态(俯仰、翻滚、偏航)和风速的影响(风速越大, 无人机平衡保持能力会有所下降), 从而造成拍摄的海面图像产生位置偏移、几何畸变等问题, 因此需要对拍摄的图像进行几何校正和图像配准。

本文采用基于三维傅里叶变换和频散关系理论的海流反演算法(第3.3节), 不同于粒子图像测速(PIV)等技术, 算法不需要选取海面的标志物。但是, 如果存在海面标志物, 可以帮助提高图像正射校正和无人机定位的精度。由于本文实验中无人机从船上起飞, 距离岸边有一定距离, 无法选取陆地作为图像校正的标志物。但是相机拍摄海面时, 部分图像会覆盖海面上的浮标, 为了分析图像校正和无人机定位的精度, 本文选取观测区域中锚定浮标的位置作为参考。

利用实验F2的数据为例, 选取一组连续的无人机和相机参数较为稳定的海面图像序列, 无人机航向角θ_u = 58°, 相机航向角θ_c = 54°, 相机倾角β大约为30°, 其余无人机飞行参数(飞行速度、飞行高度)如表 2所示。通过公式(2)—(5)对原海面图像进行几何校正和坐标转换的计算, 图 13是经过预处理之后挑选的覆盖海上浮标的子区域图像序列, 从图中可以看出, 浮标在图像序列中的位置基本不变, 图像漂移产生的定位误差在1 m以内, 具有较高的定位精度。不过, 由于海浪、海流的运动, 浮标实际上也有较小的位置偏移, 这也会带来一定的误差, 如果能进一步提高定位精度, 将更有利于海流的反演。

4.2 太阳耀斑的校正

第3.2.6节表明, 海面图像受太阳耀斑影响会造成反演的海流流速偏小。为了探究太阳耀斑对海流反演造成影响的机制, 本文从受不同程度太阳耀斑影响的海面图像的波数能量谱和灰度值分布特征进行分析, 并尝试基于灰度值分布的校正方法来减小太阳耀斑的影响。

如图 14(a)和(b)分别是选取距离太阳耀斑中心区域较远和较近位置处的子图像, 对其进行图像对比度调整以后, 做二维傅里叶变换得到波数谱F(k_x, k_y), 进一步由$\boldsymbol{\psi} ({k_x}, \;{k_y}) = {\left| {F({k_x}, \;{k_y})} \right|^2}$得到波数能量谱, 再从波数能量谱中提取谱峰最高值对应的波数即峰值波数k_px和k_py, 通过下面的公式可以计算主波波长:

$L=\frac{2 \pi}{\sqrt{k_{\mathrm{p}x}^2+k_{\mathrm{p} y}^2}}.$

(12)

图 14(c)和(d)分别是受太阳耀斑影响较小和较大时子图像的波数能量谱, 图 14(c)中峰值波数的模k_p= 0.13 rad/m($k_{\mathrm{p}}=\sqrt{k_{\mathrm{p} x}^2+k_{\mathrm{p} y}^2}$), 图 14(d)中k_p=0.17 rad/m。即受太阳耀斑影响较大时, 计算出来的峰值波数会变大, 那么对应的主波波长L变小。而从灰度值分布图来看, 受太阳耀斑影响较小的子图的灰度值分布形状较集中[图 14(e)], 而受太阳耀斑影响较大的子图的灰度值分布形状较分散[图 14(f)], 这些可能是导致海流反演较小的原因之一。

基于灰度值分布的校正方法来减小太阳耀斑影响的过程主要如下。首先, 根据2.2节中的步骤得到整个海面区域内的灰度值均值图[图 15(a)], 图 15(b)是相同y坐标下灰度值均值在x方向上的变化, 可以看出, 当子图像受太阳耀斑影响越大时, 灰度值均值也越大。用类似的方法计算灰度值方差的变化[图 15(c)和(d)], 得到的结果与前者类似。

其次, 先对实验F4的数据不拟合处理计算海表面流, 图 16(b)是得到的流场图。图中左下角黄色太阳标志表示原海面图像上受太阳耀斑影响最大的位置。将灰度值方差图[图 16(a)]与反演得到的流场[图 16(b)]旋转到同一方向下比较, 发现用于反演的子图受太阳耀斑影响越大时灰度值方差越大, 并且计算出来的海表面流精度也越差(浮标实测流速U_t=0.37~ 0.60 m/s, 流向为88°~109°)。最后, 用公式(7)对实验F4图像数据的灰度值进行拟合处理计算海表面流, 图 16(c)是拟合处理后得到的流场图。图 16(d)为拟合处理后的流速U_fit与原方法反演的流速U_o之间的差值, 从反演结果来看, 在太阳耀斑影响较小的图像区域, 经过拟合处理后的海表面流流速U_fit更加接近浮标实测流速, 而在太阳耀斑影响较大的图像区域, 反演的海表面流速变化不明显。

5 总结

小型无人机操作方便、使用灵活, 可以搭载光学相机获取高分辨率的海面图像, 观测海洋参数。本文首先对无人机光学海面图像作预处理, 获取几何校正和图像增强的海面图像, 然后利用3-D FFT方法来反演海表面流场, 反演结果与浮标实测结果一致。为了进一步研究方法的适用性, 本文分别分析了不同因素对于海表面流反演的影响, 包括图像分辨率、图像尺寸、图像幅数、时间分辨率、无人机飞行速度和飞行高度等, 主要结论有以下几个方面: (1)提高图像的空间分辨率、时间分辨率有利于海表面流的反演, 在本文的实验条件下(有效波高0.2~1.5 m, 波浪主波周期3.8~5.7 s), 图像空间分辨率为0.15 m、时间分辨率选取2.3~3.4 s, 可以获得较好的流速反演结果。(2)图像尺寸和飞行高度的选取与反演海区波浪的主波波长有关, 在本实验中, 图像尺寸大小满足波浪的主波波长的2~3倍时反演海流的精度较高、飞行高度为反演图像尺寸的3~4倍时比较合适。(3)前人在太阳耀斑对海流反演影响方面尚未开展广泛研究, 本文利用海面图像灰度值的特征统计变量如方差、均值对其做了初步研究, 结果发现, 海面图像灰度值方差越大, 反演出来的海表面流越小。最后, 本文利用二次函数拟合对海面图像灰度值进行修正, 用拟合后海面灰度图像进行海流反演, 发现拟合处理后的海面灰度图像中受太阳耀斑影响较小的区域提高了海流反演精度。

由于实验条件等限制, 目前利用无人机观测的光学海面图像反演海流还存在一些问题。第一, 本文采集无人机海面图像数据时, 无人机按照预定的航线执行飞行命令, 选取飞行状态较为稳定的图像序列, 忽略无人机俯仰角和翻滚角对几何校正的影响, 仅根据飞行高度、相机的视场角和像片倾角来对拍摄的海面图像进行几何校正。第二, 光学相机获取海面图像受太阳耀斑、云和海雾等环境影响较大。第三, 当海况较高, 海面上有大量白浪时, 可能也会影响海表面流的反演。下一步将开展更多实验, 进一步验证本文提出的无人机光学观测系统和反演算法的适用性, 以备实际观测应用。