CN116385554A - 一种基于双无人机视频拼接的近岸海域水深测绘方法 - Google Patents

Abstract

本发明公开了一种基于双无人机视频拼接的近岸海域水深测绘方法，步骤如下：完成两架无人机的相机内参标定；两架无人机携带的相机面朝大海，沿着海岸线分布，保持间距拍摄视频，保证相机视野有相交部分，通过光学相机拍摄视频记录海浪的运动特征；再通过视频拼接的技术，将两架无人机拍摄的视频进行拼接；通过无人机相机位置和姿态信息，获取每一帧的正射图片；沿垂直海岸线方向固定出一条直线作为研究区域，并得到所拍摄视频的图像关键帧的正射校正图片在这条固定直线所在位置上的时间堆栈图；最后通过cBathy水深估计方法得到近岸海域水深。

Description

一种基于双无人机视频拼接的近岸海域水深测绘方法

技术领域

本发明涉及海岸带测绘技术领域，具体涉及一种基于双无人机视频拼接的近岸海域水深测绘方法。

背景技术

近岸海域水深测绘的传统方法主要是由物理测量仪器完成测量，如声纳系统、激光雷达、合成孔径雷达和卫星图像等，但是由于这些测量仪器需要消耗大量的人力和财力，使得测绘的成本非常高。因此需要更多的低成本高效率的测绘方法来达到测量的目的。随着摄影测量技术的发展，出现了在岸边固定摄像头或激光雷达等传感器来记录海浪波的运动特征，但是这种方式依然需要在现场部署摄像头等测量工具以及需要合适的选址，受到安装地的限制同时其安装和拆除也需要较大的成本。

如今，随着自主系统无人机的研究日趋成熟，利用这种搭载有摄像头的机载系统的优势日益突出，目前亟待提出利用无人机和视频处理技术的低成本、高效率测绘方法。

发明内容

本发明的目的是为了解决现有技术中的上述缺陷，提供一种使用双无人机机载相机图像完成近岸海域水深估计的测量方法。

本发明的目的可以通过采取如下技术方案达到：

一种基于双无人机视频拼接的近岸海域水深测绘方法，所述测绘方法包括如下步骤：

S1、完成两架无人机的相机摄像头内参标定，两架无人机均搭载GPS-RTK测量模块和IMU惯性测量单元，分别记录相机位置信息和相机姿态信息；

S2、操控两架无人机沿着待测量海域的海岸线飞行，保持悬停或以相同方向保持间距相同速度匀速移动拍摄视频，保证两相机之间拍摄视频的视野有相交部分；

S3、将两架无人机的相机所拍摄的视频进行视频拼接，得到全景视频，再对全景视频进行图像预处理，其中图像预处理包括图像灰度化和图像滤波；

S4、从经过图像预处理的全景视频中选取待测绘的海域，利用所记录的相机位置信息和相机姿态信息，对所述待测绘的海域生成正射校正后的图像；

S5、利用所生成的正射校正图像，在垂直海岸线方向选择一条直线，并对其进行图像处理，得到所有帧的时间堆栈图像；

S6、通过cBathy水深估计方法，对时间堆栈图像估计出相对应像素坐标点的水深信息。

进一步地，所述步骤S1中，通过数学工具完成相机内参的标定，得到内参矩阵，内参矩阵

其中f_x和f_y是描述相机传感器在图像坐标系x轴水平方向和y轴方向竖直方向上的像素密度，x_o和y_o表示相机光轴在图像坐标系中的像素偏移量。相机的内部参数(简称内参)描述了相机本身的固有属性，决定了相机从三维场景中获取的二维图像的形状和大小，因此是进行图像处理和计算几何变换的重要输入。

进一步地，所述步骤S3过程如下：定义第一架无人机拍摄的视频为A，第二架无人机拍摄的视频为B。对于每架无人机拍摄的视频，将视频每一帧按照分辨率的宽度和高度各均分为m份，视频的每一帧分为m²个网格，并用i表示第i个网格。设F_i(t)表示视频的第i个网格的第t帧和第t+1帧之间的单应性变换矩阵，而视频的路径定义为第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，并用C_i(t)表示：

C_i(t)＝F_i(t)·F_i(t-1)···F_i(1),1≤t≤T-1,3≤T,1≤i≤m²

其中T为单个视频的总帧数，用C_i(t)表示的视频路径通过视频帧之间的变换反映出相机的运动轨迹，也记录了相机运动过程中产生抖动，能够有效描述视频的运动状态，相对应的，视频A中每个网格i的视频路径为

视频B中每个网格i的视频路径为

设定视频路径优化公式Θ(P_i)如下：

其中，Ω_t表示第t帧的邻近帧r的范围，P_i(t)表示视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，是经过C_i(t)迭代优化而来的，P_i(r)也是视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第r个单应性变换矩阵的连乘，是经过C_i(r)迭代优化而来的，将视频A每个网格i的平滑视频路径为定义为P_i ^A(t)，将视频B每个网格i的平滑视频路径为定义为P_i ^B(t)，λ表示整体权重，用于平衡||P_i(t)-C_i(t)||和||P_i(t)-P_i(r)||这两项，w_t,r用于保持处于快速位移或场景变换下的运动不连续性，通过高斯函数G计算：w_t,r＝G(||r-t||)·G(||C_i(r)-C_i(t)||)

其中，高斯函数G的定义为

x0是高斯函数的输入变量，μ为高斯函数的均值，σ则是高斯函数的标准方差，用高斯函数计算w_t,r能够反应出第t帧和第r帧的关系，w_t,r越大说明第t帧和第r帧的运动关系较密切，w_t,r越小说明第t帧和第r帧的运动关系不密切，如果场景变化过快，C_i(r)和C_i(t)的差异会变大，w_t,r会变小，视频路径优化公式Θ(P_i)使平滑视频路径专注于保持与原视频路径一致，而不会让平滑视频路径专注于快速场景切换，因此引入w_t,r的计算利于解决快速位移、旋转和场景切换的情况；视频路径优化公式的迭代是通过计算C_i(t)来更新P_i(t)：

其中，ξ为迭代指数，N(i)表示第i个网格的邻近网格的集合，j是集合N(i)中第j个网格，令

定义E^stable(P)为视频路径稳定公式：

其中，P＝{P_i(t)}，是所有网格i对应平滑视频路径P_i(t)的集合，E^stable(P)对平滑视频路径P迭代优化，让平滑视频路径接近原视频路径，减少裁剪和失真，也让每个网格的平滑视频路径更加平滑。对于视频A和视频B的平滑视频路径P^A和P^B，视频拼接公式如下：

E(P^A,P^B,H)＝E^stable(P^A)+E^stable(P^B)+β·E^stitch(P^A,P^B,H)

其中，β表示拼接系数，E^stitch表示视频的拼接方式，采用特征匹配的方法，

为视频A中第t帧的第L个特征点，

为视频B中第t帧的第L个特征点，特征点表现形式为矩阵

其中u_fp和v_fp的含义为特征点在像素坐标系上的水平坐标和垂直坐标，

和

分别表示特征点

和

所在的网格，H表示网格

和

的经由特征点

和

计算得到的单应性矩阵；

假设对第一架无人机拍摄的视频A的视频路径P^A不进行视频路径稳定优化，将视频拼接公式简化为：

E(P^A,P^B,H)＝E^stable(P^B)+β·E^stitch(P^A,P^B,H)

简化后的视频拼接公式，在迭代的过程中完成视频路径的平滑，同时使视频B能够与视频A拼接，由于只需要对平滑视频路径P_i ^B(t)迭代，使得整体计算量降低。将视频B的平滑视频路径P_i ^B(t)迭代事先指定次数后，再通过P_i ^B(t)·C_i ^B(t)^-1对视频B中的每帧每个网格进行扭曲变换，得到的结果与视频A进行逐帧融合得到全景视频。再对全景视频的每一帧进行图像灰度化以及图像滤波去噪处理。

进一步地，所述步骤S4过程如下：

设世界坐标系的坐标点为(x_w,y_w,z_w)，采用的导航坐标系为北东地坐标系，单位为米，机体坐标系的坐标(x_b,y_b,z_b)，采用的载体坐标系为前右下坐标系，无人机相机坐标系的坐标点为(x_c,y_c,z_c)，采用的摄像机坐标系为右下前坐标系，单位为米，世界空间中的点投影到相机坐标中的过程为：

其中

为相机的外部参数，描述了相机在三维场景中的位置、方向和观察角度，决定了相机从哪个角度观察场景，因此是进行三维重建和姿态估计等任务的重要输入，R_3×3为旋转矩阵，由相机姿态信息决定，设相机姿态信息的横滚角为roll,俯仰角为pitch,偏航角为yaw。那么旋转矩阵的计算公式为：

为相机坐标系相对于机体坐标系的姿态所组成的三维旋转矩阵，将相机坐标系的原点与机体坐标系的原点重合，设机体坐标系可以通过绕其z轴旋转角度θ，再绕其y轴旋转角度

最后绕其x轴旋转角度ψ与相机坐标系三轴重合，那么

t_3×1为平移矩阵，表示相机在世界坐标系下的坐标，可通过计算相机位置与世界坐标系的原点得到，而相机坐标系的坐标点(x_c,y_c,z_c)到像素平面坐标(u,v,1)的转换为：

其中K为步骤S1获得的相机内参矩阵，用于确认相机的投影性质。

相机成像模型的变换公式为：

通过以上变换，世界坐标系的点投影到二维的像素平面上，然后对图像完成正射校正。

进一步地，所述步骤S4中利用相机位置信息和相机姿态信息对选定海域生成正射校正图像的过程如下：

将北东地坐标转换到当地的沿着海岸与垂直海岸方向的坐标系，设当地的沿着海岸的方向为X轴，沿着跨海岸的方向为Y轴。将原来的北东地坐标沿着XY轴所在的平面旋转；

沿着海岸与垂直海岸两个方向定义一个矩形区域，在矩形区域中等间距的在X轴与Y轴采集世界坐标系中的三维空间点，即沿着海岸线与垂直海岸线方向采集世界坐标系中的三维空间点；

通过相机成像模型得到每个点对应的像素坐标，在三维空间中均匀地采集一系列点后，得到对应的像素坐标，然后对像素坐标进行采样插值，重新排列构成新的图像，该图像即为正射校正图片。

进一步地，所述步骤S5中，对于经过正射校正的像素图片，设像素图片总数为m0、列数为cols、行数为rows，设col的取值范围为1～cols,以第一列像素值开始直到第cols列像素值结束，依次从正射图片中第一张图片到第m0张图片的第col列像素值取出并排，组成时间堆栈图，该时间堆栈图的宽度为时间长度，列数为m0，高度为跨海岸长度，行数为rows，时间堆栈图的总张数为cols。常用的水深测绘法多为基于频域的方法，上面步骤处理得到的时间堆栈图有利于从中分析出海浪波的传播方向、传播速度和传播频率等物理量，使用此数据结构也是适用于基于频域的方法来进行水深估计。

进一步地，所述步骤S6过程如下：根据色散关系式中在线性波浪中波浪角频率ω、波浪波数k、水深h和重力加速度g之间存在以下关系：

ω²＝gktanh(kh)

推导得到水深h的计算公式如下：

另外，波浪速度c、波浪波数k和波浪角频率ω、波浪频率f之间存在恒等式：

ω＝2πf，那么水深公式亦可写成：

根据上述公式，仅需知道波浪速度c、波浪波数k、波浪角频率ω和波浪频率f这四个物理量中的其中两个，即可推导出水深的估计值。然后基于步骤S5获得的时间堆栈图，使用cBathy水深估计方法，完成对波浪频率f和波浪波数k的估计从而计算出水深h。cBathy水深估计法是一种二维的频域估计法，主要包括信号处理和容错处理，由于存在大量非线性拟合优化步骤与反馈环节，在估计参数时设置多个因素加权来进行考量，并且在最终得到的参数上采用容错校验操作，所以该算法的鲁棒性较强。

本发明相对于现有技术具有如下的优点及效果：

1.本发明通过双无人机移动录制视频数据进行近岸海域水深测绘，测绘的海域范围可通过无人机沿着海岸线移动不断更新测绘区域，同时记录已测绘区域，以此扩大测绘区域。而已有的方法是通过一架无人机定点悬停在空中或者使用固定点相机进行测绘。相比之下，本发明的测绘范围要大，测绘效率要高。

2.由于单架无人机相机录制视频时，左右远端处的跨海岸海域距离无人机相机较远且拍摄角度不佳，存在信息丢失影响测绘结果的问题。引入两架无人机的相机，令相机拍摄到距前一架无人机一端较远的跨海岸海域，使得整体获取的视频信息尽可能完整。

3.本发明在对无人机录制的视频进行拼接时选择使用捆绑视频路径和联合路径优化的方法，并不是将每个视频同一个时间戳的对应帧进行拼接的传统图像拼接算法，而是结合视频的路径进行联合优化，消除图像拼接后产生的抖动和重影。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明实施例中公开的一种基于双无人机视频拼接的近岸海域水深测绘方法的流程图；

图2是本发明实施例中公开的一种基于双无人机视频拼接的近岸海域水深测绘方法中视频拼接后截取的一帧效果图；

图3是本发明实施例中公开的一种基于双无人机视频拼接的近岸海域水深测绘方法中选定的海域范围图(200米×200米)；

图4是本发明实施例中公开的一种基于双无人机视频拼接的近岸海域水深测绘方法中的正射校正图；

图5本发明实施例中公开的一种基于双无人机视频拼接的近岸海域水深测绘方法中的时间堆栈图；

图6本发明实施例一种基于双无人机视频拼接的近岸海域水深测绘方法的水深估计结果图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

如图1所示，本实施例公开的一种基于双无人机视频拼接的近岸海域水深测绘方法，包括如下步骤：完成两架无人机的相机摄像头内参标定、操作两架无人机保持间距拍摄视频、对两架无人机拍摄的视频进行全景视频拼接、通过每一帧的相机信息得到正射校正图片、划定好研究海域并获取研究海域的时间堆栈图、使用cBathy水深估计法估计出近岸海域水深。下面结合图1具体介绍实施过程：

步骤S1、通过OpenCV库等数学工具计算出两架无人机机载相机的内参矩阵，内参矩阵

其中f_x和f_y是描述相机传感器在图像坐标系x轴水平方向和y轴方向竖直方向上的像素密度，x_o和y_o表示相机光轴在图像坐标系中的像素偏移量。常见的方法有张正友标定法，为求便捷，亦可使用Matlab自带的工具箱完成内参的标定。

步骤S2、由于单架无人机距离远处的跨海岸海域较远，海浪波会显示地比较密集，如果此时海上有雾气或有短风波影响，使得无人机相机拍摄的信噪比过低。操控第二架无人机的另一个目的是令其拍摄到第一架无人机拍摄不佳的海域，提升拍摄到该海域的信噪比。在需要测绘水深的海滩区域内，操作两架无人机，使其分布满足一定要求，如两架无人机沿着海岸线飞行、两架无人机的相机面向大海、两架无人机的相机拍摄范围需存在公共视野、两架无人机在拍摄视频的过程中同时保持悬停或者沿着同一方向匀速移动。其中留有公共视野的目的是为了方便后续视频拼接的操作，以扩大测绘视角，两架无人机拍摄的视频的公共视野至少占单个无人机相机视野的30％，以及无人机移动的速度不超过1米每秒。

步骤S3、通过捆绑视频路径和联合优化路径的方法对两架无人机拍摄的视频进行全景视频拼接，具体操作如下：

定义第一架无人机拍摄的视频为A，第二架无人机拍摄的视频为B。对于每架无人机拍摄的视频，将视频每一帧按照分辨率的宽度和高度各均分为m份，视频的每一帧分为m²个网格，并用i表示第i个网格。设F_i(t)表示视频的第i个网格的第t帧和第t+1帧之间的单应性变换矩阵，而视频的路径定义为第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，并用C_i(t)表示：

C_i(t)＝F_i(t)·F_i(t-1)···F_i(1),1≤t≤T-1,3≤T,1≤i≤m²

其中T为单个视频的总帧数。相对应的，视频A中每个网格i的视频路径为

视频B中每个网格i的视频路径为

设定视频路径优化公式Θ(P_i)如下：

其中，Ω_t表示第t帧的邻近帧r的范围，这个范围取±30帧，P_i(t)表示视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，是经过C_i(t)迭代优化而来的，P_i(r)也是视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第r个单应性变换矩阵的连乘，是经过C_i(r)迭代优化而来的。将视频A每个网格i的平滑视频路径定义为P_i ^A(t)，将视频B每个网格i的平滑视频路径定义为P_i ^B(t)，λ表示整体权重，用于平衡||P_i(t)-C_i(t)||和||P_i(t)-P_i(r)||这两项，初始值设为5，如果希望优化后的视频整体跟原视频相差不大，可将λ调小，如果希望优化后的视频减少失真，可将λ调大。而w_t,r用于保持处于快速位移或场景变换下的运动不连续性，通过高斯函数G计算：w_t,r＝G(||r-t||)·G(||C_i(r)-C_i(t)||)

其中，高斯函数G的定义为

x0是高斯函数的变量，μ为高斯函数的均值，一般取0，σ则是高斯函数的标准方差，初始值可设为10。视频路径优化公式的迭代是通过计算C_i(t)来更新P_i(t)：

其中，ξ为迭代指数，N(i)表示第i个网格的邻近网格的集合，j是集合N(i)中第j个网格，令

定义E^stable(P)为视频路径稳定公式：

其中，P＝{P_i(t)}，是所有网格i对应平滑视频路径P_i(t)的集合。通过迭代P_i(t)，视频的抖动将得到有效去除，这是实现视频拼接稳定的一个前提。对于视频A和视频B的平滑视频路径P^A和P^B，视频拼接公式如下：

E(P^A,P^B,H)＝E^stable(P^A)+E^stable(P^B)+β·E^stitch(P^A,P^B,H)

其中，β表示拼接系数，E^stitch表示视频的拼接方式，使用的是特征匹配的方法来拼接，为视频A中第t帧的第L个特征点，

为视频B中第t帧的第L个特征点，特征点表现形式为矩阵

其中u_fp和v_fp的含义为特征点在像素坐标系上的水平坐标和垂直坐标，

和

分别表示特征点

和

所在的网格，H表示网格

和

的经由特征点

和

计算得到的单应性矩阵；

其中，选择其中第一架无人机拍摄的视频不进行视频路径优化，第二架无人机拍摄的视频需要路径优化，视频拼接公式简化为：

E(P^A,P^B,H)＝E^stable(P^B)+β·E^stitch(P^A,P^B,H)

通过简化后的视频拼接公式，将视频B的平滑视频路径P_i ^B(t)迭代事先指定次数时，再通过P_i ^B(t)·C_i ^B(t)^-1对视频B中的每帧每个网格进行扭曲变换，得到的结果与视频A逐帧图像融合得到全景视频后进行图像灰度化处理、高斯滤波去噪操作。

注意这里对视频的路径优化实际上某种程度上讲是对图像的一种单应性变换，这种变换会使得S1步骤标定的相机内参不可用，从而无法使三维世界坐标映射到二维像素坐标，或者映射过程中出现很大的误差。因此才需要对其中一个无人机拍摄的视频不进行路径优化，以保证其不受单应性变换的影响，确保后续计算可以用到该无人机的相机内参。

当使用联合路径优化的方法完成视频拼接后，整个全景视频共享一个内参，而这个内参正是没有经过路径优化过的第一架无人机拍摄的视频A对应的内参，后续在三维世界坐标划定研究范围映射到二维平面时所需要的内参就满足了条件。然后，全景视频的正射校正将会根据第一架无人机的相机位置信息和相机姿态信息进行后续操作。如果视频拼接的过程中选择的是第二架无人机拍摄的视频B不进行视频路径优化，第一架无人机拍摄的视频进行视频路径优化，那么相对应地，视频拼接公式就变成：

E(P^A,P^B,H)＝E^stable(P^A)+β·E^stitch(P^A,P^B,H)

后续全景视频的处理将会根据第二架无人机的数据进行操作。

步骤S4、通过三维空间到二维平面的投影关系，对于世界坐标转换到图像坐标系的变换矩阵：

其中K为步骤S1获得的相机内参矩阵，其中

为相机的外部参数，R_3×3为旋转矩阵，由相机姿态信息决定，设相机姿态信息的横滚角为roll,俯仰角为pitch,偏航角为yaw。那么旋转矩阵的计算公式为：

为相机坐标系相对于机体坐标系的姿态所组成的三维旋转矩阵，将相机坐标系的原点与机体坐标系的原点重合，设机体坐标系可以通过绕其z轴旋转角度θ，再绕其y轴旋转角度

最后绕其x轴旋转角度ψ与相机坐标系三轴重合，那么

其中，t_3×1为平移矩阵，表示相机在世界坐标系下的坐标，可通过计算相机位置与世界坐标系的原点得到。通过以上变换，最终可以将三维空间中的点投影到二维图像平面上，然后对图像完成正射校正。根据当前帧的相机位置信息和相机姿态信息，可以划定需要进行正射校正的区域。

由于需要通过图像空间中的信息确定海域范围信息，因此需要完成对图像的正射校正。为表示方便起见，将北东地坐标转换到当地的沿着海岸与跨海岸方向的坐标系。设当地的沿着海岸的方向为X轴方向，沿着跨海岸方向为Y轴方向。只需要将原来的北东地坐标系沿着XY轴所在的平面旋转一定角度即可，这个角度可以通过当地的指南针方位角计算出，具体角度可依照当地实际情况取舍。

当坐标系的转换得到后，此时只需要沿着海岸线与跨海岸方向两个方向定义一个矩形区域即可。

当矩形区域选择完毕后，由于需要生成最终的正射校正图片，采取这样一种简单的方式：通过在前面定义的矩形区域中等间距的在X轴与Y轴，即沿着海岸线与跨海岸方向，采集三维空间点。

通过上述方法，此时所要研究的区域的三维空间点坐标此时便可以完全得到，通过相机成像模型便可以得到每个点对应的像素坐标。当在三维空间中均匀地采集一系列点后，其对应的像素坐标也已得到，然后对像素坐标进行插值采样，插值算法可以选择双三次插值或双线性插值，重排使其构成新的图像，该图像即为正射校正图片，此外需要注意的是，三维空间中距离采样之间的间隔一般为0.5米左右，一般不会超过1米，否则正射校正图片会丢失很多信息影响水深估计的精度。采样间隔可根据无人机相机的焦距、分辨率和相机距海面的高度决定的。因此对像素坐标重采样后的正射图片，像素与像素之间存在着距离信息。

步骤S5、全景视频经过正射校正后，在垂直海岸线方向即跨海岸方向选择一条直线，这条直线在三维空间某一个维度上具有固定的坐标，这样通过图像空间想三维空间逆变换时可以得到唯一解，从而在图像空间上可以度量三维空间上的位置信息。然后利用所选区域生成全景视频的校正图片并进行处理得到所有帧在该区域的一个时间堆栈图像。

由于已经得到所有的正射校正图片，为了方便转换为时间序列相关性分析需要的数据结构，在所有经过正射校正后的图片上选择一列像素值(三维空间上表现为垂直海岸线的跨海岸方向)组成时间堆栈图像。具体操作为：对于经过正射变换的像素图片，设像素图片的总数为m0、列数为cols、行数为rows，设col的取值范围为1～cols,以第一列像素值开始直到第cols列像素值结束，依次从正射图片中第一张图片到第m0张图片的第col列像素值取出并排，组成时间堆栈图。该时间堆栈图像的宽度为时间长度，高度正比于跨海岸的长度。

步骤S6、根据上述得到的时间堆栈图，可以使用cBathy水深估计法估计水深：

根据色散关系式中在线性波浪中波浪角频率ω、波浪波数k、水深h和重力加速度g之间存在以下关系：ω²＝gktanh(kh)

推导得到水深h的计算公式如下：

另外，波浪速度c、波浪波数k和波浪角频率ω、波浪频率f之间存在恒等式：

ω＝2πf

那么水深公式进一步写成：

基于上一个步骤获得的时间堆栈图，使用cBathy水深估计方法，完成对波浪频率f和波浪波数k的计算估计出水深h，水深估计结果的示意图如图6所示。其中cBathy水深估计法于2013年发表自期刊Journal of geophysical research:Oceans，第118卷，2595-2609页，文章标题：cBathy:A robust algorithm for estimating nearshore bathymetry。

实施例2

本实施例继续公开一种基于双无人机视频拼接的近岸海域水深测绘方法，包括如下步骤：

S1、准备两架无人机，为了分别获取相机的位置信息和姿态信息，两架无人机需要携带高精度GPS-RTK装置和IMU惯性测量单元。通过Matlab数学工具或OpenCV库计算无人机机载相机的内参矩阵，内参矩阵

其中f_x和f_y是描述相机传感器在图像坐标系x轴水平方向和y轴方向竖直方向上的像素密度，x_o和y_o表示相机光轴在图像坐标系中的像素偏移量。

S2、在需要测绘水深的海滩区域内，操作无人机，飞行高度50～150米，使其分布满足一定要求，如两架无人机沿着海岸线分布、两架无人机的相机面向大海、两架无人机的相机拍摄范围存在公共视野、两架无人机在拍摄视频的过程中同时保持悬停或者沿着同一方向保持间距匀速移动，视频的拍摄对象是往海岸或沙滩方向传播的周期性海浪波。

S3、将多架无人机拍摄的视频数据通过捆绑视频路径和联合优化路径的方法进行视频拼接：

定义第一架无人机拍摄的视频为A，第二架无人机拍摄的视频为B，对于每架无人机拍摄的视频，将视频每一帧按照分辨率的宽度和高度各均分为m份，视频的每一帧分为m²个网格，并用i表示第i个网格。在计算过程中如果视频分辨率过大，可将m设为10，视频分辨率小，则设m为4。设F_i(t)表示视频的第i个网格的第t帧和第t+1帧之间的单应性变换矩阵，而视频的路径定义为第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，并用C_i(t)表示：

C_i(t)＝F_i(t)·F_i(t-1)···F_i(1),1≤t≤T-1,3≤T,1≤i≤m²

其中T为单个视频的总帧数，相对应的，视频A中每个网格i的视频路径为

视频B中每个网格i的视频路径为

设定视频路径优化公式Θ(P_i)，其定义为：

其中，Ω_t表示第t帧的邻近帧r的范围，P_i(t)表示视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，是经过C_i(t)迭代优化而来的，P_i(t)的初始值可以等于C_i(t)，P_i(r)也是视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第r个单应性变换矩阵的连乘，是经过C_i(r)迭代优化而来的，相对应的，将视频A每个网格i的平滑视频路径定义为P_i ^A(t)，将视频B每个网格i的平滑视频路径定义为P_i ^B(t)，λ表示整体权重，初始值可设为5，用于平衡||P_i(t)-C_i(t)||和||P_i(t)-P_i(r)||这两项，而w_t,r用于保持处于快速位移或场景变换下的运动不连续性，通过高斯函数G计算：

w_t,r＝G(||r-t||)·G(||C_i(r)-C_i(t)||)

其中，高斯函数G的定义为

x0是高斯函数的变量，μ为高斯函数的均值，可为0，σ则是高斯函数的标准方差，视频路径优化公式的迭代是通过计算C_i(t)来更新P_i(t)：

其中，ξ为迭代指数，N(i)表示第i个网格的邻近网格的集合，j是集合N(i)中第j个网格，令

定义E^stable(P)为视频路径稳定公式：

其中，P＝{P_i(t)}，是所有网格i对应平滑视频路径P_i(t)的集合，对于视频A和视频B的平滑视频路径P^A和P^B，视频拼接公式如下：

E(P^A,P^B,H)＝E^stable(P^A)+E^stable(P^B)+β·E^stitch(P^A,P^B,H)

其中，β表示拼接系数，用于平衡E^stable和E^stitch，如果希望拼接的视频稳定、去除抖动，则可以将β设小一点，让E^stable的权重占比大一点，如果希望视频拼接的效果好、重影少，则可以将β设大一点，让E^stitch的权重占比大一点，如果希望综合效果好，可取β＝0.01，E^stitch表示视频的拼接方式，使用的是特征匹配的方法来拼接，

为视频A中第t帧的第L个特征点，

为视频B中第t帧的第L个特征点，特征点表现形式为矩阵

其中u_fp和v_fp的含义为特征点在像素坐标系上的水平坐标和垂直坐标。

和

分别表示的是特征点

和

所在的网格。H表示网格

和

的经由特征点

和

计算得到的单应性矩阵；

假设对第一架无人机拍摄的视频A的视频路径P^A不进行视频路径稳定优化，将视频拼接公式简化为：

E(P^A,P^B,H)＝E^stable(P^B)+β·E^stitch(P^A,P^B,H)

通过简化后的视频拼接公式，将视频B的平滑视频路径P_i ^B(t)迭代事先指定次数时，再通过P_i ^B(t)·C_i ^B(t)^-1对视频B中的每个网格进行扭曲变换，得到的结果与视频A逐帧图像融合得到全景视频。然后对全景视频进行预处理，主要包括：图像灰度化、高斯滤波去噪，其具体操作可以使用Matlab图像处理函数、调用OpenCV开源视觉库或者通过手写源代码实现。

其中选择某一个视频不进行视频路径稳定优化是因为视频路径优化本质上是一种扭曲变换或者单应性变换，会改变视频图像的扭曲程度和透射视角，导致后续计算正射校正区域时所需要的相机内参矩阵K不可用，选择第一架无人机拍摄的视频A的视频路径P^A不进行视频路径稳定优化是为了后续正射校正能够使用到第一架无人机的飞行数据，如相机的位置信息和相机的姿态信息。

此外如果单纯使用传统图像拼接的方法，使用特征点匹配的方法将两个视频逐帧缝合，那么拼接的效果会出现极大的抖动以及重影，也无法使用该视频拼接结果生成正射校正图像，因此引入的捆绑视频路径和联合路径优化的方法能很好解决这个问题。

S4、通过三维空间到二维平面的投影关系，对于世界坐标转换到图像坐标系的变换矩阵：

其中

为相机的外部参数，R_3×3为旋转矩阵，由相机姿态信息决定，设相机姿态信息的横滚角为roll,俯仰角为pitch,偏航角为yaw。那么旋转矩阵的计算公式为：

为相机坐标系相对于机体坐标系的姿态所组成的三维旋转矩阵，将相机坐标系的原点与机体坐标系的原点重合，设机体坐标系可以通过绕其z轴旋转角度θ，再绕其y轴旋转角度

最后绕其x轴旋转角度ψ与相机坐标系三轴重合，那么

t_3×1为平移矩阵，表示相机在世界坐标系下的坐标，可通过计算相机位置与世界坐标系的原点得到。通过如上的一系列的变换矩阵操作，最终可以将三维世界坐标系中的点投影到二维平面像素坐标系上，根据当前帧的相机位置和姿态信息，选取所要研究的水深估计区域，利用前述的操作对此区域生成正射校正后的图像：

由于需要通过图像空间中的信息确定海域范围信息，因此需要完成对图像的正射校正。如图3全景图像所示沿着海岸线方向划定200米长度，沿着跨海岸方向划定200米长度，面积为40000平方米，为表示方便起见，将北东地坐标转换到当地的沿着海岸与跨海岸方向的坐标系。设当地的沿着海岸的方向为X轴，沿着跨海岸的方向为Y轴，只需要将原来的北东地坐标系沿着XY轴所在的平面旋转一定角度即可，这个角度可以通过当地的指南针方位角计算出，具体角度可依照当地实际情况取舍。

当坐标系的转换得到后，此时只需要沿着海岸线与跨海岸方向两个方向定义一个矩形区域即可。

当矩形区域选择完毕后，由于需要生成最终的正射校正图片，采取这样一种简单的方式：通过在前面定义的矩形区域中等间距的在X轴与Y轴即沿着海岸线与跨海岸方向，采集三维空间点。

通过上述方法，此时所要研究的区域的三维空间点坐标此时便可以完全得到，通过相机成像模型便可以得到每个点对应的像素坐标。当在三维空间中均匀地采集一系列点后，其对应的像素坐标也已得到，然后对像素坐标进行插值采样，重排使其重新构成新的图像，该图像即为正射校正图片，此外需要注意的是，三维空间中距离采样之间的间隔一般为0.5米左右，一般根据无人机机载相机的分辨率和飞机的高度决定的。因此对像素坐标重采样后的正射图片，像素与像素之间存在着距离信息。

S5、全景视频经过正射校正后，在垂直海岸线方向(即跨海岸方向)选择一条直线，这条直线在三维空间某一个维度上具有固定的坐标，这样通过图像空间想三维空间逆变换时可以得到唯一解，从而在图像空间上可以度量三维空间上的位置信息。然后利用所选区域生成全景视频的校正图片并进行处理得到所有帧在该区域的一个时间堆栈图像。

由于已经得到所有的正射校正图片，为了方便转换为时间序列相关性分析需要的数据结构，对于经过正射校正的像素图片，设像素图片的总数为m0、列数为cols、行数为rows，设col的取值范围为1～cols，以第一列像素值开始直到第cols列像素值结束，依次从正射图片中第一张图片到第m0张图片的第col列像素值取出并排，即在所有经过正射校正后的图片上选择一列像素值(三维空间上表现为垂直海岸线的跨海岸方向)组成时间堆栈图像。该时间堆栈图像的宽度为时间长度(单位为秒)，高度正比于跨海岸的长度。

S6、cBathy水深估计法是一种二维的频域估计法，主要包括信号处理和容错处理，由于存在大量非线性拟合优化步骤与反馈环节，在估计参数时设置多个因素加权来进行考量，并且在最终得到的参数上采用容错校验操作，所以该算法鲁棒性较强。根据上个步骤得到的时间堆栈图，可以使用cBathy水深估计法来估计出水深h。cBathy水深估计法于2013年发表自期刊Journal of geophysical research:Oceans，第118卷，2595-2609页，文章标题：cBathy:A robust algorithm for estimating nearshore bathymetry。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (7)

1.一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述测绘方法包括如下步骤：

S1、完成两架无人机的相机摄像头内参标定，两架无人机均搭载GPS-RTK测量模块和IMU惯性测量单元，分别记录相机位置信息和相机姿态信息；

S2、操控两架无人机沿着待测量海域的海岸线飞行，保持悬停或以相同方向保持间距相同速度匀速移动拍摄视频，保证两相机之间拍摄视频的视野有相交部分；

S3、将两架无人机的相机所拍摄的视频进行视频拼接，得到全景视频，再对全景视频进行图像预处理，其中图像预处理包括图像灰度化和图像滤波；

S4、从经过图像预处理的全景视频中选取待测绘的海域，利用所记录的相机位置信息和相机姿态信息，对此海域生成正射校正后的图像；

S5、利用所生成的正射校正图像，在垂直海岸线方向选择一条直线后进行图像处理，得到所有帧的时间堆栈图像；

S6、使用cBathy水深估计方法，通过时间堆栈图像估计出相对应像素坐标点的水深信息。

2.根据权利要求1所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述步骤S1中，使用图像处理工具得到无人机相机的内参矩阵，而内参矩阵

其中f_x和f_y是描述相机传感器在图像坐标系x轴水平方向和y轴方向竖直方向上的像素密度，x_o和y_o表示相机光轴在图像坐标系中的像素偏移量。

3.根据权利要求1所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述步骤S3过程如下：

定义第一架无人机拍摄的视频为A，第二架无人机拍摄的视频为B。对于每架无人机拍摄的视频，将视频每一帧按照分辨率的宽度和高度各均分为m份，视频的每一帧分为m²个网格，并用i表示第i个网格，设F_i(t)表示视频的第i个网格的第t帧和第t+1帧之间的单应性变换矩阵，而视频的路径定义为第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，并用C_i(t)表示：

C_i(t)＝F_i(t)·F_i(t-1)···F_i(1),1≤t≤T-1,3≤T,1≤i≤m²

其中T为单个视频的总帧数，设定视频A中每个网格i的视频路径为C_i ^A(t)，视频B中每个网格i的视频路径为C_i ^B(t)，视频路径优化公式Θ(P_i)如下：

其中，Ω_t表示第t帧的邻近帧r的范围，P_i(t)表示视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第t个单应性变换矩阵的连乘，是经过C_i(t)迭代优化得到，P_i(r)也是视频第i个网格的平滑视频路径，是视频的第i个网格从第1个单应性变换矩阵到第r个单应性变换矩阵的连乘，是经过C_i(r)迭代优化得到，将视频A每个网格i的平滑视频路径为定义为P_i ^A(t)，将视频B每个网格i的平滑视频路径为定义为P_i ^B(t)，λ表示整体权重，用于平衡||P_i(t)-C_i(t)||和||P_i(t)-P_i(r)||这两项，w_t,r用于保持处于快速位移或场景变换下的运动不连续性，通过高斯函数G计算：

w_t,r＝G(||r-t||)·G(||C_i(r)-C_i(t)||)

其中，高斯函数G的定义为

x0是高斯函数的自变量，μ为高斯函数的均值，σ则是高斯函数的标准方差；视频路径优化公式的迭代是通过计算C_i(t)来更新P_i(t)：

其中，ξ为迭代指数，N(i)表示第i个网格的邻近网格的集合，j是集合N(i)中第j个网格，令

定义E^stable(P)为视频路径稳定公式：

其中，P＝{P_i(t)}，是所有网格i对应平滑视频路径P_i(t)的集合，对于视频A和视频B的平滑视频路径P^A和P^B，视频拼接公式如下：

E(P^A,P^B,H)＝E^stable(P^A)+E^stable(P^B)+β·E^stitch(P^A,P^B,H)

其中，β表示拼接系数，E^stitch表示视频的拼接方式，采用特征匹配的方法，

为视频A中第t帧的第L个特征点，

为视频B中第t帧的第L个特征点，特征点表现形式为矩阵

其中u_fp和v_fp的含义为特征点在像素坐标系上的水平坐标和垂直坐标，

和

分别表示特征点

和

所在的网格，H表示网格

和

的经由特征点

和

计算得到的单应性矩阵；

假设对第一架无人机拍摄的视频A的视频路径P^A不进行视频路径稳定优化，将视频拼接公式简化为：

E(P^A,P^B,H)＝E^stable(P^B)+β·E^stitch(P^A,P^B,H)

通过简化后的视频拼接公式，将视频B的平滑视频路径P_i ^B(t)迭代事先指定次数后，再通过P_i ^B(t)·C_i ^B(t)^-1对视频B中的每帧每个网格进行扭曲变换，得到的结果与视频A进行逐帧融合得到全景视频，再对全景视频的每一帧进行图像灰度化以及图像滤波去噪处理。

4.根据权利要求1所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述步骤S4过程如下：

设世界坐标系的坐标点为(x_w,y_w,z_w)，采用的导航坐标系为北东地坐标系，机体坐标系的坐标(x_b,y_b,z_b)，采用的载体坐标系为前右下坐标系，无人机相机坐标系的坐标点为(x_c,y_c,z_c)，采用的摄像机坐标系为右下前坐标系，世界空间中的点投影到相机坐标中的过程为：

其中

为相机的外部参数，R_3×3为旋转矩阵，由相机姿态信息决定，设相机姿态信息的横滚角为roll，俯仰角为pitch，偏航角为yaw，那么旋转矩阵的计算公式为：

为相机坐标系相对于机体坐标系的姿态所组成的三维旋转矩阵，将相机坐标系的原点与机体坐标系的原点重合，设机体坐标系可以通过绕其z轴旋转角度θ，再绕其y轴旋转角度

最后绕其x轴旋转角度ψ与相机坐标系三轴重合，那么

t_3×1为平移矩阵，表示相机在世界坐标系下的坐标，可通过计算相机位置与世界坐标系的原点得到，而相机坐标系的坐标点(x_c,y_c,z_c)到像素平面坐标(u,v,1)的转换为：

其中K为相机内参矩阵，用于确认相机的投影性质；相机成像模型的变换公式为：

通过以上变换，世界坐标系的点投影到二维的像素平面上，然后对图像完成正射校正。

5.根据权力要求4所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述步骤S4中利用相机位置信息和相机姿态信息对选定海域生成正射校正图像的过程如下：

将北东地坐标转换到当地的沿着海岸与垂直海岸方向的坐标系，设当地的沿着海岸的方向为X轴，沿着跨海岸的方向为Y轴，将原来的北东地坐标沿着XY轴所在的平面旋转；

沿着海岸与垂直海岸两个方向定义一个矩形区域，在矩形区域中等间距的在X轴与Y轴采集世界坐标系中的三维空间点，即沿着海岸线与垂直海岸线方向采集世界坐标系中的三维空间点；

通过相机成像模型得到每个点对应的像素坐标，在三维空间中均匀地采集一系列点后，得到对应的像素坐标，然后对像素坐标进行采样插值，重新排列构成新的图像，该图像即为正射校正图片。

6.根据权利要求1所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，所述步骤S5中，对于经过正射校正的像素图片，设像素图片的总数为m0、列数为cols、行数为rows，设col的取值范围为1～cols，以第一列像素值开始直到第cols列像素值结束，依次从正射图片中第一张图片到第m0张图片的第col列像素值取出并排，组成时间堆栈图，该时间堆栈图的宽度为时间长度，列数为m0，高度为跨海岸长度，行数为rows，时间堆栈图的总张数为cols。

7.根据权利要求1所述的一种基于双无人机视频拼接的近岸海域水深测绘方法，其特征在于，步骤S6过程如下：

根据色散关系式中在线性波浪中波浪角频率ω、波浪波数k、水深h和重力加速度g之间存在以下关系：ω²＝gktanh(kh)

推导得到水深h的计算公式如下：

另外，波浪速度c、波浪波数k和波浪角频率ω、波浪频率f之间存在恒等式：

ω＝2πf，水深公式进一步写成：

根据上述公式，仅需知道波浪速度c、波浪波数k、波浪角频率ω和波浪频率f这四个物理量中的其中两个，即可推导出水深的估计值，然后基于步骤S5获得的时间堆栈图，使用cBathy水深估计算法，完成对波浪频率f和波浪波数k的估计从而计算出水深h。

Images