CN102446367B - 基于多波束声纳海底测量数据构建三维地形矢量模型的方法 - Google Patents

Abstract

基于多波束声纳海底测量数据构建三维地形矢量模型的方法，它涉及一种用多波束声纳的测量数据进行三维矢量地形模型构建的方法。它解决了多波束测距传感器只能提供离散地形高程且仅能查询测量点地形参数的缺点。用于非结构化海底地形建模和矢量数据检索。该方法首先采集海底区域地形的采样点水平坐标及所对应的高程数据，根据多波束声纳地形测量数据的特征定义地形三角网模型的数据结构，最后根据多波束声纳地形测量的各点高程数据和Delaunay三角网的数据结构，建立Delaunay三角网失量模型。

Description

基于多波束声纳海底测量数据构建三维地形矢量模型的方法

技术领域

本发明涉及是一种基于地形测量数据构建三维地形模型的方法，特别是涉及一种多波束声纳的测量数据进行三维矢量地形模型构建的方法。

背景技术

由于海洋环境地形复杂、海况恶劣且存在难以预测的作业风险，无法采用人工方式直接完成海底地形测量与作业，近年来基于无人自主水下航行器AUV(AutonomousUnderwater Vehicle)及其载荷的声学(如多波束声纳)等传感器实施海底地形、地貌测量作业，已成为海洋勘测与海洋学调查领域内安全有效的手段之一和技术发展趋势。然而，声学传感器的专用软件仅能为用户提供被测量地形的显示界面，其软件平台的依赖性使得测量数据无法直接生成海底地形模型(特别是生成矢量化的地形模型)，无法基于离散测量数据完成对非结构化地形的连续处理与分析，难以达到海洋环境勘测与地形绘制的作业目标。因此，寻求基于声纳测量数据构建非结构化海底的矢量地形模型，用于可视化地形显示与数据分析，在AUV及其声学传感器工程实际应用与海洋勘测技术研究领域，就显得尤为迫切。

目前，国内外试图寻求一种有效的手段用于多波束声纳地形测量数据的建模。如高金耀等人在文献《多波束数据的海底数字地形模型构建》(海洋通报，2003年，第23卷第1期)中，在权重配赋网格插值基础上，提出分形fBM和张力样条配合使用的方程式插值方法，构建了格网数据模型，该方法适合于结构化地形环境构建。马建林等人在文献《多波束测深海底数字地形模型的建立》(海洋测绘，2005年，第25卷第5期)中针对多波束数据的特点，采用改进的IDW算法构建规则格网地形数据模型，该方法适合于结构化地形环境构建。海军海洋测绘研究所的贾俊涛等人在文献《基于海量多波束数据的海底地形模型的构建与可视化》(测绘科学技术学报，2008年，第25卷第4期)中亦采用格网数据模型构建海底地形，该方法适合结构化环境。上述研究中所提出的地形模型构建方法均仅适合于结构化环境建模，而海底地形均呈现非结构化特性，因此，这些结构化地形建模方法难以满足地形测绘的工程实际应用需求。

上海大学的顾晨等人在发明专利《基于多波束声纳数据的水下三维地形重建方法》中公开了一种基于多波束声纳数据的水下三维地形重建方法，其特征在于：首先对Seabeam2000多波束声纳系统采集的测深数据进行内插，计算得到各点的波束与水平线之间的夹角

波束射程R_s和波束入射角α；然后根据多波束声纳系统采集的回波强度数据I_s与波束入射角α拟合Lambertian反射模型参数ρ，利用有限差分法校正波束与水平线之间的夹角

迭代N次，使估算出的回波强度数据与采集的回波强度数据平均相对误差绝对值最小，计算此时的测深数据，重建水下三维地形，否则就迭代拟合模型参数来估算回波强度数据。该研究方法提高了测深数据的空间分辨率，侧重于数据处理，而不是地形建模算法。

发明内容

本发明的目的是提供一种基于多波束声纳海底测量数据构建三维地形矢量模型的方法，以解决多波束测距传感器只能提供离散地形高程且仅能查询测量点地形参数的缺点。本发明包括下述步骤：

一、利用无人自主水下航行器采用多波束声纳，采集海底区域地形的采样点水平坐标及所对应的高程数据，作为多波束声纳地形测量数据点；

二、构建Delaunay三角网；具体过程如下：对所有的多波束声纳地形测量数据点，按照水平地理坐标的横坐标X和纵坐标Y，进行排列；连接横坐标X值最小的两点p₁和p₂，并将此两点中纵坐标y值大的点指向纵坐标y值小的点的方向作为初始基线的方向；

对前述的初始基线的两端点p₁和p₂，按照其地理坐标p₁(x₁，y₁)，p₂(x₂，y₂)，确定初始基线p₁p₂的方程为：F(x，y)＝y-(ax-b)                                (7)

式中a＝(y₂-y₁)/(x₂-x₁)，b＝(x₁y₂-x₂y₁)/(x₂-x₁)；

将任意一个多波束声纳测量数据点p_j中地形的地理坐标x_j、y_j带入初始基线p₁p₂的方程F(x，y)中，判断F(x，y)的正负，从而决定点p_j是否被扩展；

若F(x，y)＞0，p_j点位于初始基线p₁p₂的正区，则p_j点为待扩展的点；

若F(x，y)＜0，则p_j点位于初始基线p₁p₂的负区，则p_j点不扩展；

若F(x，y)＝0，则p_j点位于初始基线p₁p₂上，则p_j点不扩展；

在所有待扩展的点中找到某一点p₃，使点p₃和点p₂的连线p₃p₂与点p₃和点p₁的连线p₃p₁间张角最大，将该点p₃作为第三点；

连接点p₂与点p₃，并以p₂p₃作为基线并按照指向纵坐标y值小的方向确定基线的方向，利用上述基线正负区域判别方法查找第四点的待选集合，利用张角最大准则确定第四点；重复上述操作直到所有的多波束声纳地形测量数据点都被编入Delaunay三角网；

三、根据多波束声纳地形测量数据的结构特征及Delaunay三角网算法数据结构的需求，建立Delaunay三角网的点数据结构、边数据结构和三角形面片数据结构。

本发明的有益效果在于针对多波束声纳的地形测量数据，通过定义矢量模型的拓扑数据结构并提出基于分治点排序思想的基线生成与基线快速扩展方法，所构建的三维海底地形矢量模型，适合于非结构化海底地形建模和矢量数据检索。由于建立Delaunay三角网矢量模型，使本发明解决多波束测距传感器只能提供离散地形高程且仅能查询测量点地形参数的缺点。

附图说明

图1为多波束声纳地形测量数据点的位置计算原理图，图中的附图标记为H的是换能器；

图2为多波束船参考坐标系与当地水平坐标系示意图；

图3为点数据结构示意图；

图4为边数据结构示意图；

图5为三角形面片数据结构示意图；

图6-1和图6-2是确定基线过程示意图；

图7为优化的Delaunay三角网生长算法中初始基线生成过程流程图；

图8为在基线外快速搜索第三点的流程图；

图9为线模式显示的海试区域海底地形矢量模型图；

图10为面模式显示的海试区域海底地形矢量模型图。

具体实施方式

具体实施方式一：

下面结合图1至图8具体说明本实施方式。本实施方式包括下述步骤：

一、利用无人自主水下航行器采用多波束声纳，采集海底区域地形的采样点水平坐标及所对应的高程数据，作为多波束声纳地形测量数据点；

二、构建Delaunay三角网；具体过程如下：对所有的多波束声纳地形测量数据点，按照水平地理坐标的横坐标X和纵坐标Y，进行排列；连接横坐标X值最小的两点p₁和p₂，并将此两点中纵坐标y值大的点指向纵坐标y值小的点的方向作为初始基线的方向；

对前述的初始基线的两端点p₁和p₂，按照其地理坐标p₁(x₁，y₁)，p₂(x₂，y₂)，确定初始基线p₁p₂的方程为：F(x，y)＝y-(ax-b)                     (7)

式中a＝(y₂-y₁)/(x₂-x₁)，b＝(x₁y₂-x₂y₁)/(x₂-x₁)；

将任意一个多波束声纳测量数据点p_j中地形的地理坐标x_j、y_j代入初始基线p₁p₂的方程F(x，y)中，判断F(x，y)的正负，从而决定点p_j是否被扩展；

若F(x，y)＞0，p_j点位于初始基线p₁p₂的正区，则p_j点为待扩展的点；

若F(x，y)＜0，则p_j点位于初始基线p₁p₂的负区，则p_j点不扩展；

若F(x，y)＝0，则p_j点位于初始基线p₁p₂上，则p_j点不扩展；

在所有待扩展的点中找到某一点p₃，使点p₃和点p₂的连线p₃p₂与点p₃和点p₁的连线p₃p₁间张角最大，将该点p₃作为第三点；

连接点p₂与点p₃，并以p₂p₃作为基线并按照指向纵坐标y值小的方向确定基线的方向，利用上述基线正负区域判别方法查找第四点的待选集合，利用张角最大准则确定第四点；重复上述操作直到所有的多波束声纳地形测量数据点都被编入Delaunay三角网；

三、根据多波束声纳地形测量数据的结构特征及Delaunay三角网算法数据结构的需求，建立Delaunay三角网的点数据结构、边数据结构和三角形面片数据结构。

下面对本发明进行详细介绍：

(一)地形三角网模型及其数据结构定义

1)地形三角网模型描述

根据数字高程模型(Digital Elevation Model，DEM)的数学定义，本发明所构建的海底地形三角网模型可以描述为：海底地形区域D中任意一点地形数据P_j由多波束声纳测量得到，包括三个数据元素(x_j，y_j，H_j)，前述的地形测量数据P_j按不规则三角网生成规则ζ所连接成的三角形面片M_i的集合便构成了地形三角网模型MAP，其表达式如下：

MAP＝{M_i＝ζ(P_j)|P_j(x_j，y_j，H_j)∈D，j＝1，…，n；i＝1，…，m}(1)

式中x_j、y_j——多波束声纳测量数据中地形的地理坐标；

H_j——多波束声纳测量数据中地形的高程数据。

前述的多波束声纳地形测量数据P_j按下面步骤计算得到：

①建立多波束船参考坐标系和当地水平坐标系，如表1所示。

表1多波束船参考坐标系与当地水平坐标系

LLS的原点E与VFS的原点o重合，三坐标轴不同，两者之间的转换关系为：

$\left\{\begin{array}{c}{X}_{\mathrm{LLS}}=x_{\mathrm{VFS}}\left(\cosh \cos p\right)+y_{\mathrm{VFS}}(\cosh \sin p\sin r-\sinh \cos r)+z_{\mathrm{VFS}}(\cosh \sin p\cos r+\sinh \sin r)\\ Y_{\mathrm{LLS}}=x_{\mathrm{VFS}}\left(\sinh \cos p\right)+y_{\mathrm{VFS}}(\sinh \sin p\sin r+\cosh \cos r)+z_{\mathrm{VFS}}(\sinh \sin p\cos r-\cosh \sin r)\\ Z_{\mathrm{LLS}}=x_{\mathrm{VFS}}(-\sin p)+y_{\mathrm{VFS}}\left(\cos p\sin r\right)+z_{\mathrm{VFS}}\left(\cos p\cos r\right)\end{array}\right.---\left(2\right)$

式中(X_LLS，Y_LLS，Z_LLS)、(x_VFS，y_VFS，z_VFS)——LLS和VFS下任意一点的坐标；

h——艏摇角，ox轴在水平面XEY上的投影与EX轴间的夹角，见图2；

P——纵摇角，ox轴与水平面XEY之间的夹角；

R——横摇角，xoz平面与通过ox轴的垂直平面xoZ之间的夹角。

②处理姿态传感器测量数据

首先剔除姿态数据中的异常值，然后对删除(或缺失)的点，根据前后两次采集数据进行插值。

③处理多波束声纳传感器测量数据

针对多波束声纳的二进制(BINARY)编码方式，参照多波束采集系统配有的数据输出格式说明指导书，根据多波束声纳原始地形测量数据的特点定义2个字节整型、4个字节整型、4个字节浮点型三个联合进行多波束地形测量数据(波束入射角θ₀、波束的射程R)进行提取和格式转换。

④进行测量数据的空间配准

VFS的原点o位于换能器中心时，波束脚印在VFS下的坐标是：

[xyz]^T _VFS＝[0Rsinθ₀Rcosθ₀]^T                       (3)

由式(2)和(3)，根据三角和差公式化简得到合并后的波束脚印计算式：

$\left\{\begin{array}{c}{x}_j=R\sin p\cosh \cos ({\mathrm{\θ}}_0-r)-R\sinh \sin ({\mathrm{\θ}}_0-r)\\ y_j=R\sin p\sinh \cos ({\mathrm{\θ}}_0-r)+R\cosh \sin ({\mathrm{\θ}}_0-r)\\ z_j=R\cos p\cos ({\mathrm{\θ}}_0-r)\end{array}\right.---\left(4\right)$

式中θ₀、R——分别为波束入射角、波束的射程，见图1；

h、r、p——分别为艏摇角，横摇角，纵摇角。

多波束声纳测量数据中地形的高程数据计算式为：

H_j＝h_tide-(z_LLS+h_ss+h_ds)＝h_tide-Rcospcos(θ₀-r)-h_ss-h_ds                 (5)

式中H_j——多波束声纳测量数据中地形的高程数据；

h_ss、h_ds、h_tide——静态吃水(常量)、动态吃水(测得的Heaven值)、潮位。

⑤进行测量数据的时间配准

若姿态传感器在t₀时刻未采集到数据，只在t₀的邻近时刻t_i、t_i+1采集到数据时，基于线性插值的时间配准公式，进行时间配准。

$y_0=\frac{t_0-t_i}{t_{i+1}-t_i}{y}_{i+1}+\frac{t_{i+1}-t_0}{t_{i+1}-t_i}{y}_i---\left(6\right)$

式中y₀——t₀时刻姿态传感器采集到的数据；

y_i、y_i+1——分别为t_i、t_i+1时刻的姿态传感器采集到的数据。

2)数据结构建立

根据多波束声纳地形测量数据的特征以及Delaunay三角网算法的数据结构需求，建立点数据结构、边数据结构、三角形面片数据结构来描述式(1)的地形三角网模型。

①点数据结构：存储多波束声纳地形测量数据点的结构体。该结构体记录多波束声纳地形测量数据的地理坐标和高程值。

点数据结构如下，见图3。

②边数据结构：存储边信息的结构体。该结构体记录边的索引、边的两个端点、边的长度以及边的使用次数。构网过程中以该结构体中的边作为基线进行扩展，形成新边。边数据结构如下，见图4。

③三角形面片数据结构：存储三角形信息的结构体。该结构体记录三角形的索引、三角形三条边的索引、三角形三个顶点的指针、三角形所在面的法向量。三角形数据结构存储所生成的Delaunay三角网模型中三角形之间的拓扑关系。

三角形面片数据结构如下，见图5。

(二)基于分治点排序思想生成初始基线

对多波束声纳测量数据中地形的地理坐标x_j、y_j，按x_j坐标从小到大进行排序：

vertex[0]，vertex[1]，...，vertex[i]，...，vertex[MAX-1]。

其中i为排序后点的索引值。

比较vertex[0].y和vertex[1].y的大小，若vertex[0].y＞vertex[1].y，vertex[0]为第一点，记为p₁；vertex[1]为第二点，记为p₂；若vertex[0].y＜vertex[1].y，vertex[1]为第一点，记为p₁；vertex[0]为第二点，记为p₂。

连接p₁和p₂，生成初始基线p₁p₂，并规定从p₁到p₂的方向是初始基线p₁p₂的前进方向。

图6是确定基线的示意图，箭头方向是初始基线的前进方向。

(三)基于正负区域判别原理及张角最大准则快速扩展基线

①基于正负区域判别原理查找第三点

对前述的初始基线两端点p₁和p₂，设其地理坐标为p₁(x₁，y₁)，p₂(x₂，y₂)，初始基线p₁p₂的方程为：

F(x，y)＝y-(ax-b)                      (7)

式中  a＝(y₂-y₁)/(x₂-x₁)，b＝(x₁y₂-x₂y₁)/(x₂-x₁)。

将任意一个多波束声纳测量数据p_j中地形的地理坐标x_j、y_j带入初始基线p₁p₂的方程F(x，y)中，判断F(x，y)的正负，决定点p_j是否被扩展。

若F(x，y)＞0，说明p_j点位于初始基线p₁p₂的正区，则p_j点为待扩展的点；

若F(x，y)＜0，则p_j点位于初始基线p₁p₂的负区，则p_j点不扩展。

若F(x，y)＝0，则p_j点位于初始基线p₁p₂上，则p_j点不扩展。

②张角最大准则确定第三点

在初始基线p₁p₂正区中找到与初始基线p₁p₂张角最大，即余弦最小的点，作为第三点p₃。

针对初始基线p₁p₂与基线正区中某点p_j构成的三角形ΔP₁P₂P_j，利用余弦定理计算式(8)，计算待扩展点p_j与初始基线p₁p₂的夹角余弦，找到位于初始基线p₁p₂正区中，且与初始基线p₁p₂夹角余弦值最小的点，选取该点为第三点，记为p₃(x₃，y₃)。

$P_1{P}_2^2=P_1{P}_j^2+P_2{P}_j^2-2P_1{P}_j\·P_2{P}_j\·\cos \∠P_1{P}_j{P}_2,\cos \∠P_1{P}_j{P}_2=\frac{P_1{P}_j^2+P_2{P}_j^2-P_1{P}_2^2}{2P_1{P}_j\·P_2{P}_j}---\left(8\right)$

式中  P₁、P₂、P_j——分别为ΔP₁P₂P_j的三个顶点；

∠P₁P_jP₂——初始基线p₁p₂所对的角。

③定义边矢量vector<Edge>edges，存储已经参与构网的边，定义三角形面片矢量vector<Triangle>triangles，用于存储已经构成的三角形信息。

连接初始基线p₁p₂两个端点p₁、p₂与第三点p₃，构成第一个三角形ΔP₁P₂P₃，将p₁p₂存入edges中，将Δp₁p₂p₃的三个顶点存入triangles。将形成的两条新边p₁p₃和p₂p₃作为新的基线，进行扩展。

扩展新基线时，先判断该基线的两端点是否与edges中某条边的两端点相同，如果不同，则该边未标记，可以继续参与构网；如果相同，则该边已标记，停止。

④重复所述步骤(①～③)，直到所有多波束声纳地形测量数据都被遍历，都被编入三角网，且生成的三角网符合Delaunay三角网为止。

对多波束声纳所有地形测量数据进行排序，可以保证扩展基线总是沿着基线前进方向的左边进行，限定了搜索方向，减少了搜索范围，缩短了搜索时间。图7是生成初始基线流程图；图8为在基线外快速搜索并确定第三点的流程图。

具体实施方式二：

下面结合图9和10具体说明本实施方式。为了验证本发明方法的合理性、可行性，2010年10月30日在大连进行DVL(Doppler Velocity Log，多普勒计程仪)辅助导航海试试验时，获取了一组原始多波束声纳地形测量数据并存储于.837文件中。该实施例包含以下步骤：

先使用多波束声纳自带软件DeltaT_TCP中的Data_Output项将获得的.837文件转换为.83P文件。.83P文件中用于表达每个Ping数据的字节是0～1215，其中256～424字节及1122～1215字节内容为0，剔除这部分冗余信息。针对425～1121字节数据，每隔24个字节提取一个原始多波束声纳地形测量数据。然后依据前述多波束声纳地形测量数据P_j的计算步骤，计算得到多波束声纳地形测量数据P_j的地理坐标x_j、y_j和高程值H_j。最后针对多波束声纳地形测量数据P_j，基于本发明中优化的Delaunay三角网生长算法构建海试区域的三维海底地形矢量模型。

由于数值仿真计算量较大，本发明采用了Visual C++2005平台作为仿真环境。图9是线模式显示的海试区域海底地形矢量模型，图10是面模式显示的海试区域海底地形矢量模型。

Claims (2)

1.基于多波束声纳海底测量数据构建三维地形矢量模型的方法，其特征在于它包括下述步骤：

一、利用无人自主水下航行器采用多波束声纳，采集海底区域地形的采样点水平坐标及所对应的高程数据，作为多波束声纳地形测量数据点；

二、构建Delaunay三角网；具体过程如下：对所有的多波束声纳地形测量数据点，按照水平地理坐标的横坐标X和纵坐标Y，进行排列；连接横坐标X值最小的两点p₁和p₂，并将此两点中纵坐标y值大的点指向纵坐标y值小的点的方向作为初始基线的方向；

对前述的初始基线的两端点p₁和p₂，按照其地理坐标p₁(x₁，y₁)，p₂(x₂，y₂)，确定初始基线p₁p₂的方程为：F(x，y)＝y-(ax-b)                            (7)

式中a＝(y₂-y₁)/(x₂-x₁)，b＝(x₁y₂-x₂y₁)/(x₂-x₁)；

将任意一个多波束声纳测量数据点p_j中地形的地理坐标x_j、y_j带入初始基线p₁p₂的方程F(x，y)中，判断F(x，y)的正负，从而决定点p_j是否被扩展；

若F(x，y)＞0，p_j点位于初始基线p₁p₂的正区，则p_j点为待扩展的点；

若F(x，y)＜0，则p_j点位于初始基线p₁p₂的负区，则p_j点不扩展；

若F(x，y)＝0，则p_j点位于初始基线p₁p₂上，则p_j点不扩展；

在所有待扩展的点中找到某一点p₃，使点p₃和点p₂的连线p₃p₂与点p₃和点p₁的连线p₃p₁间张角最大，将该点p₃作为第三点；

连接点p₂与点p₃，并以p₂p₃作为基线并按照指向纵坐标y值小的方向确定基线的方向，利用上述基线正负区域判别方法查找第四点的待选集合，利用张角最大准则确定第四点；重复上述操作直到所有的多波束声纳地形测量数据点都被编入Delaunay三角网；

三、根据多波束声纳地形测量数据的结构特征及Delaunay三角网算法数据结构的需求，建立Delaunay三角网的点数据结构、边数据结构和三角形面片数据结构。

2.根据权利要求1所述的基于多波束声纳海底测量数据构建三维地形矢量模型的方法，其特征在于多波束声纳地形测量数据点通过下述步骤获得：①建立多波束船参考坐标系和当地水平坐标系；②处理姿态传感器测量数据：首先剔除数据中的异常值，然后对删除或缺失的点，根据前后两次采集数据进行插值；③处理多波束声纳传感器测量数据；④进行测量数据的空间配准；⑤进行测量数据的时间配准。

Images