CN105258684A - 基于激光点云为约束的多波束低掠射角波束归位方法 - Google Patents

Abstract

本发明公开了一种基于激光点云为约束的多波束低掠射角波束归位方法，属于海洋测绘技术领域，其首先根据沿岸区域内的若干个声速剖面和沿航迹线附近内的表层声速值，利用经验正交函数分析法构建区域内浅层修正后的三维声速场，完成测深数据初始声速改正；然后以同目标激光点云提取的特征点为最外侧声波束位置约束，获取初始声速改正后特征点的坐标偏差，据此计算低掠射层内初始旋转角和初始放大系数，进而按角度线性加权内插60度以上剩余各波束对应的旋转角和放大系数，最终实现波束几何改正。本发明通过这种办法，实现了浅层海水中声波束的准确归位,解决了水岸线上下点云无缝拼接的问题。

Description

基于激光点云为约束的多波束低掠射角波束归位方法

技术领域

本发明涉及海洋测绘技术领域，具体涉及一种基于激光点云为约束的多波束低掠射角波束归位方法。

背景技术

远海岛礁、海岸、滩涂等浅水区海底地形及近岸水上地形测量一直是国内外测绘领域的一个热点。为此，国内外相继出现了以船载激光技术和机载雷达技术为代表的海陆一体化测量手段。其中机载雷达技术由于存在研发难度大、价格昂贵，且激光在海水中传播能量损失快等问题，因此在实际的测量工作中难以得到应用。而搭载多波束的船载激光技术则利用激光和声波分别测量水上和水下，同时获取近岸水上地形和浅水区海底地形，与传统的靠单波束和RTK测量方式相比，具有易于集成、测量精度高、效率快等优势。但是一方面由于多波束换能器声波发射具有一定的开角局限性，如R2Sonic2024开角最大为160度，因此水下部分不能够得到完整的测量；另一方面，由于收到波束导向作用的影响，边缘波束脚印分辨率相对较低，使得水下部分地形不能得到精细测量。鉴于此，很多研发单位在安装多波束换能器时偏转一定的角度β，使得水下靠近水面处的部分能够得到精确的测量。但在后期数据处理时，声速改正部分仍旧按照以往换能器正常安装时的方式利用一个或几个声速剖面根据时间就近原则或则距离就近原则进行校正，忽视了声速的时空可变性；声速的变化主要受水中温度、盐度、压强和密度等因素的影响，理想状态下，假设水中密度、温度、压强、盐度等均没有变化，多波束换能器发射的声波的路径应为一条直线，但是，实际情况下这些因素都是存在且变化的，由其是海水表层易受光照、海面风等不可抗拒因素的影响，表层中水平与垂直方向上易形成温度不一、密度不同的水团，从而导致海水表层中的声速复杂多变且无规律可循，利用常规的声速剖面进行声线跟踪不能反映声波的实际传播路径，从而导致水上下同目标坐标存在异常差值，即使存在微小的声速误差，因波束入射角大，也对波束垂直位置产生较大影响，在水上下点云拼接时导致存在较大的缝隙,声速不规则变化导致水上下同目标错位示意图(如图3所示)。

由于声速剖面在时间上和空间上的可变性，采用有限的声速剖面数据必然不能满足水深数据进行声速改正的精度要求，从而引入测深误差。利用经验正交函数(empiricalorthogonalfunction，简称EOF)能够利用海洋中声速剖面数据的时空相关性，通过对具有相同特性的声速剖面群(同一水域内多个声速剖面)进行特征向量提取(模态分解)，进而结合采样数据构建测区内三维声速场。LRLeblanc等的研究表明EOF是描述声速剖面最有效的基函数，近年来的相关研究证明了这种方法的可行性和有效性。此外，因为海水表层声速受到多方因素的影响，声速结构复杂多变，所以可根据NewBrunswick大学的海洋制图组提出的一种经验改正法，对多波束测深横剖面上点叠加统计，用几何旋转的方法实现校正。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种基于激光点云为约束的多波束低掠射角波束归位方法，设计合理，克服了现有技术的不足，解决了水岸线上下点云无缝拼接问题。

为了实现上述目的，本发明采用如下技术方案：

一种基于激光点云为约束的多波束低掠射角波束归位方法，采用船载水岸线一体化测量系统，其水上部分包括三维激光扫描仪和组合导航系统，水下部分包括多波束测深仪，多波束测深仪带有具有一定偏转角度的多波速换能器，水上部分与水下部分通过一个平台支架固连；其中的三维激光扫描仪用于测量水面以上地物地貌，多波束测深仪用于测量水下地物地貌，组合导航系统用于为三维激光扫描仪和多波束测深仪提供定位信息、时间信息、姿态信息和航向信息；船载水岸线一体化测量系统还包括表层声速仪和声速剖面仪两个辅助设备；

所述的多波束低掠射角波束归位方法按照如下步骤进行：

步骤1：对船载水岸线一体化测量系统进行标定，确定多波束测深仪、三维激光扫描仪和组合导航系统之间的空间相对位置关系；

步骤2：选取一段具有水上下公共目标的测区，并利用船载水岸线一体化测量系统沿岸匀速进行测量，获取水上下点云同一坐标系下坐标；

步骤3：对存储的多波束数据按时间分区，分区间隔为10秒，分别记为A₁,A₂,A₃...A_t；

步骤4：根据表层声速仪和GPS数据存储的时间标签，获取每一个分区内行船附近的平均表层声速值c(A_ts)和船附近内一点的地理坐标(x_At,y_At)；

步骤5：利用声速剖面仪均匀获取沿岸测量区域内多个声速剖面数据svp_n，并记录每一个采样点的地理位置(x_N,y_N)；

步骤6：对步骤5中得到的声速剖面数据svp_n进行层化处理，获得不同位置点处同一水层上的声速c(z_i)和每一层声速的平均值

步骤7：通过步骤6层化处理后的声速剖面基于公式(1)得到测区内每一个点处的声速值：

$c_n(x,y,z)=\stackrel{\‾}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\Σ}}{\mathrm{\α}}_k(x,y)f_k\left(z_i\right)---\left(1\right)$

其中加权系数α_k(x,y)为地理坐标；f_k(z_i)为正交函数；为每一层声速的平均值；

步骤8：根据各分区内平均表层声速值c(A_ts)，对该区域内(x_At,y_At)处拟合的声速剖面表层部分(多波束换能器吃水深度D)进行修正；

步骤9：根据各区域修正后的声速剖面，进行声速改正，获取多波束换能器坐标系下水上下同目标特征点处激光点与多波束点的平均坐标差，计算低掠射层内初始几何旋转角和初始放大系数进而按角度线性加权内插计算低掠射层内60度以上剩余波束角对应的旋转角和放大系数

步骤10：对低掠射层的固定层部分内的波束进行初始旋转角为和初始放大系数为的几何改正，对60度以上的剩余波束即自由层部分，进行旋转角为放大系数为的几何改正。

优选地，在所述的步骤3中，对存储的多波束数据按时间进行分区的具体方法如下：

步骤3.1：记录船载水岸线一体化测量系统开始工作和结束工作的时间，精确到秒；

步骤3.2：然后根据带有时间标签的每一个发射接收周期多波束信息进行工作区间上的分区，分区间隔为t秒。

优选地，在所述的步骤4中，一个分区内行船附近的平均表层声速值c(A_ts)和船附近一点的地理坐标(x_At,y_At)的具体获取方法如下：

步骤4.1：根据带有时间标签的每一个发射接收周期多波束信息和对应每一个分区时间节点，找到每一分区内记录的声速值c_i，根据公式(2)则可获取每一个分区内行船附近的平均表层声速值c(A_ts)，即：

c(A_ts)＝(c₁+c₂+…+c_m)/m；(2)

步骤4.2：根据GPS记录数据的时间和对应每一个分区时间节点，获取分区时间中间节点处的坐标(b,l,h)，并转到地理坐标系下得到(x_At,y_At)。

优选地，在所述的步骤7中，获得每一点处的声速c_n(x,y,z_i)和加权系数α_k(x,y)的具体方法如下：

步骤7.1：将根据步骤6得到的多个声速剖面数据写成声速值矩阵，形式如式(13)，再根据公式(15)求得协方差矩阵；

步骤7.2：求取声速值协方差矩阵的特征值和对应的特征向量，特征向量即为经验正交函数；

步骤7.3：将特征值按照从大到小排列，求取各个特征值的贡献率，选择较大的特征值所对应的特征向量，利用正交三角分解的方法在最小二乘意义下求取α_k(x,y)，这样就可以利用求得的经验正交函数和X＝[a₀,a₁,a₂,a₃,]^T系数矩阵，表示测深区域任意一点的声速剖面。

优选地，在所述的步骤8中，对分区内(x_At,y_At)处拟合的声速剖面表层部分(多波束换能器吃水深度D)进行修正具体获得方法如下：

步骤8.1：根据公式(3)计算区域A_t内(x_At,y_At)处深度D(多波束换能器吃水深度)处声速修正量

$\begin{array}{c}{c}_n(x_{At},y_{At},D)=\stackrel{\‾}{c}\left(D\right)+\underset{k=0}{\overset{k}{\Σ}}{\mathrm{\α}}_k(x_{At},y_{At})f_k\left(D\right)\\ {\mathrm{\Δc}}_t^0=c\left(A_{ts}\right)-c_n(x_{At},y_{At},D)\end{array}---\left(3\right);$

步骤8.2：根据以及公式(4)计算修正后的测区内每一点的声速值即：

$\begin{array}{c}{c}_n(x,y,z)=\stackrel{\‾}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\Σ}}{\mathrm{\α}}_k(x,y)f_k\left(z_i\right)+{\mathrm{\Δc}}_t^0(z_i\≤D)\\ c_n(x,y,z)=\stackrel{\‾}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\Σ}}{\mathrm{\α}}_k(x,y)f_k\left(z_i\right)(z_i>D)\end{array}---\left(4\right).$

优选地，在所述的步骤9中，低掠射层的定义如下：

低掠射层由多波束换能器偏转方向一侧60度以上波束角组成，其包括固定层和自由层两部分，固定层部分由2倍多波束换能器最边缘波束与水平线之间夹角φ构成，自由层部分则由60度以上剩余部分角度组成。

根据多波束换能器坐标系下水上下同目标特征点处激光点与多波束点的平均坐标差，计算低掠射层固定层部分内波束初始几何旋转角和初始放大系数以及自由层部分各波束对应的旋转角和放大系数获得的具体方法如下：

步骤9.1：在区域A_t内选取水上下形状规则的目标，拾取换能器坐标中心下K组约束激光点和多波束测深点坐标值(x_l,y_l,z_l)和(x_e,y_e,z_e)，根据几何关系，首先得到XOY平面内平均旋转角θ₁；

$\begin{array}{c}{\mathrm{\θ}}_e=\left(\underset{e=1}{\overset{K}{\Σ}}\arctan \frac{y_e}{x_e}\right)/K\\ {\mathrm{\θ}}_l=\left(\underset{l=1}{\overset{K}{\Σ}}\arctan \frac{y_l}{x_l}\right)/K\\ {\mathrm{\θ}}_1={\mathrm{\θ}}_l-{\mathrm{\θ}}_e\end{array}---\left(5\right);$

步骤9.2：根据几何关系，得到与XOY平面垂直方向上的平均旋转角θ₂：

$\begin{array}{c}{d}_l=\sqrt{x_l^2+y_l^2}\\ {\mathrm{\θ}}_2=\left(\underset{l=1}{\overset{K}{\Σ}}\arctan \frac{z_l}{d_l}\right)/K\end{array}---\left(6\right);$

则初始几何旋转角即：

${\mathrm{\θ}}_{At}^0=\left(\begin{array}{c}{\mathrm{\θ}}_1\\ {\mathrm{\θ}}_2\end{array}\right)---\left(7\right);$

假设换能器工作时开角为向右上方偏转角度为β，扫岸一侧某波束角为γ_i，则可以得到低掠射层内自由层部分各波束角对应的旋转角

${\mathrm{\&theta;}}_{At}^i=\left\{\begin{array}{c}{\mathrm{\&theta;}}_{At}^0(180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;}\&le;{\mathrm{\&gamma;}}_i\&le;\frac{\mathrm{\&alpha;}}{2})\\ \frac{{\mathrm{\&gamma;}}_i+\mathrm{\&beta;}-60}{120-\frac{\mathrm{\&alpha;}}{2}-\mathrm{\&beta;}}\&CenterDot;{\mathrm{\&theta;}}_{At}^0(60-\mathrm{\&beta;}<{\mathrm{\&gamma;}}_i<180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;})\end{array}\right.---\left(8\right);$

步骤9.3：根据两种点云的坐标，计算初始放大系数得：

同理可得低掠射层内自由层部分各波束角对应的放大系数

优选地，在所述的步骤10中，进行几何改正的具体方法如下：

步骤10.1：根据步骤9中获得的各波束放大系数，首先对低掠射层内波束点进行距离改正：

步骤10.2：根据步骤9中获得的各波束旋转角，对低掠射层内各波束点进行角度改正，最终得到几何改正后的坐标(x'_e,y'_e,z'_e)：

$\left\{\begin{array}{c}{x}_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\cos (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1\left)\right)\\ y_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\sin (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1\left)\right)\\ z_e^{,}=s_e^i\&CenterDot;\sin \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\end{array}\right.---\left(12\right).$

本发明所带来的有益技术效果：

本发明提出了一种基于激光点云为约束的多波束低掠射角波束归位方法，与现有技术相比，一种基于激光点云为约束的多波束低掠射角波束归位方法，考虑到沿岸区域内海水声速的时空变化变化，并结合表层声速数据对拟合的声速剖面中浅层声速进行了修正；结合三维激光扫描仪与多波束测得的同目标点云对多波束60度以上的声波进行几何改正，纠正了多波束每一发射和接收中掠射角较低的波束的位置，解决了船载水岸线一体化测量系统中无缝拼接问题。

附图说明

图1为本发明一种基于激光点云为约束的多波束低掠射角波束归位方法的流程框图。

图2为本发明中船载水岸线一体化测量系统的工作示意图。

图3为声速不规则变化导致水上下同目标错位的示意图。

图4为本发明中几何改正的示意图。

图5为本发明中低掠射层角度改正角度范围的示意图。

图6为本发明中声速剖面站的分布示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

实施例：

首先介绍一下本发明所需要的理论基础：

1、多波束勘测声速剖面场的经验正交函数分析方法

国内外研究表明在近海和远洋海域，声速变化呈现一定的规律性，在一定程度上是可以预估的。声速剖面结构参数化模型能够实现对声速剖面数据的简化描述与扩展，经验正交函数分析方法(empiricalorthogonalfunction,EOF)也称特征向量分析(eigenvectoranalysis)或主成分分析(principalcomponentanalysis,PCA)，是分析矩阵数据的结构特性、提取主要数据特征量的一种方法,目前在气象学研究中有较多应用。在多波束勘测声速剖面场的数学模型建立中，EOF展开的基本思想就是将声速c_n要素场序列分解成正交的时间函数与正交的空间函数的乘积之和,把声速剖面样本采集点的地理坐标看作声速场的权重系数,其他要素看作典型的声速场。对于浅水的情况下，由于深度有限，水面附近的声速己成为声速剖面的重要组成部分，局部对流己经变成了重要的热交换。阳光，大气温度对于浅水温度的影响很大，使得声速剖面在浅水环境下变化相当剧烈，沈海远等则利用实测数据证实在一定区域内的浅水声速剖面已经验证叫函数(EOF)近似表示，可以达到较好的精度。其基本原理如下：

设经过层化处理到的N个声速剖面等深度点处的深度为ih(i＝0,1,…,K)，将N个声速剖面样本表示为矩阵的形式：

$C=\left|\begin{array}{cccc}{c}_1\left(0\right)& c_2\left(0\right)& \mathrm{...}& c_3\left(0\right)\\ c_1\left(h\right)& c_2\left(h\right)& \mathrm{...}& c_3\left(h\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ c_1\left(kh\right)& c_2\left(kh\right)& \mathrm{...}& c_3\left(kh\right)\end{array}\right|---\left(13\right)$

定义C的协方差矩阵为R_c，R_c可表示为：

$R_c=\frac{1}{N}\underset{N=1}{\overset{N}{\&Sigma;}}{C}_N{C}_N^T---\left(14\right)$

其中的每一个元素r₁₁为：

$r_{ij}=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\&lsqb;c_n\left(z_i\right)-\stackrel{\&OverBar;}{c}\left(z_i\right)\&rsqb;*\&lsqb;c_n\left(z_j\right)-\stackrel{\&OverBar;}{c}\left(z_j\right)\&rsqb;,(i,j=0,L,K)---\left(15\right)$

上式中，为N个等深度声速采样数据的平均值。

求出此协方差矩阵的特征值和特征向量，其特征向量就是所求的经验正交函数,有：

RF＝λF(16)

式中，λ为协方差矩阵的特征值，F为与特征值相对应的纵向量，即为特征向量，则：

F(z_ih)＝(f₀(z_ih),…,f_k(z_ih)),i＝0,…,K(17)

f_k(z_ih)，k,j＝0,…K就是所确定的经验正交函数。于是，声速场中任意一点的声速可以表示为：

$c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{K}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)---\left(18\right)$

式中，加权平均系数α_k(x,y)是地理坐标(x,y)的函数。上式写成矩阵的形式为：

$\left[\begin{array}{c}{c}_i\left(0\right)\\ .\\ .\\ .\\ c_i\left(Kh\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&OverBar;}{c}\left(0\right)\\ .\\ .\\ .\\ \stackrel{\&OverBar;}{c}\left(Kh\right)\end{array}\right]+\&lsqb;f_0\left(z\right)\mathrm{...}{f}_K\left(z\right)\&rsqb;\left[\begin{array}{c}{\mathrm{\&alpha;}}_0(x,y)\\ .\\ .\\ .\\ {\mathrm{\&alpha;}}_K(x,y)\end{array}\right]---\left(19\right)$

假设欲确定的α_k(x,y)函数的形式为：

α(x,y)＝a₀+a₁x+a₂y+a₃xy(20)

其中，a₀、a₁、a₂、a₃是待求的参数。对同一层c_i(k)，有：

$\left[\begin{array}{c}{\mathrm{\&alpha;}}_1(x_1,y_1)\\ .\\ .\\ .\\ {\mathrm{\&alpha;}}_N(x_N,y_N)\end{array}\right]=\left[\begin{array}{cccc}1& x_1& y_1& x_1{y}_2\\ 1& x_2& y_2& x_2{y}_2\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& x_N& y_N& x_N{y}_N\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right]---\left(21\right)$

用矩阵表示为：

α(x,y)＝AX(22)

根据最小二乘法原理，可得：

X＝(A^TPA)^-1(A^TPα)(23)

式中，X＝[a₀,a₁,a₂,a₃]^T即为所求定的参数。

为了衡量选取的正交函数对声速剖面的贡献程度，引入方差贡献计算公式：

${\mathrm{Cont}}_j={\mathrm{\&lambda;}}_j/\underset{i=1}{\overset{K}{\&Sigma;}}{\mathrm{\&lambda;}}_i---\left(24\right)$

式中，λ_j为前几阶空间函数对应的特征值，K为特征值求解的个数。

公式(19)是由K+1个未知数、K+1个线性方程组成的，因此可求出不同的函数α_k(x,y),然后可以依照式(21)进行拟合。若求得的矩阵A中的特征值有少数几个特征值的绝对值远远大于其余的特征值，这时用较少的参数来表征声速剖面是可行的。可以由大到小选定较大绝对值的特征值所确定的经验正交函数。假设由此而选定了p个经验正交函数，一般选定p＜＜K，方程组则变为由p个未知数、K+1个线性方程组成的一个超定的系统，对此利用QR分解的方法来求出在最小二乘意义下的最优结果。于是可用求得的参数拟合不同的声速剖面。

2、多波束声速改正之几何改正原理

声速改正后处理是为进一步提高多波束声速改正精度而引入的室内多波束侧身数据的处理方法，其本质是利用更精确或更接近真实的声速剖面以及在无法确切了解声速结构变化的情况下通过引入一些原则对多波束测深数据进行在处理。多波束声速改正后处理方法可以分为两类。一类是以改变声速剖面为思路的处理方法，它涉及对多波束折射路径解的重新计算；在一类是根据声速对波束折射路径解的控制规律进行声速改正，即几何改正。

在多波束测量过程中有时会面临无法确知声速结构变化的状况。这种状况常发生在因某些原因(如联合调查，船尾有其他设备或者电缆拖拽等)测量过程中无法经常停船进行声速测量或在声速变化复杂的区域作业。针对这种情况，NewBrunswick大学的海洋制图组建立了一种经验改正法，以改正声速不准产生的假地形或地形畸变。这种方法不从波束折射路径解方法入手，而是通过对波束在测深横剖面上的叠加统计，用几何旋转的方法改正地形畸变，从而达到改善多波束测量精度的目的。这种方法基于如下假设：

(1)通过叠加足够数量的波束测深横剖面，所有局部地形特征均可以被正负抵消二平均化；

(2)声速不准确导致的地形畸变可以通过波束在测深横剖面上的单一旋转来改正。

这种方法通常包括三个步骤：

①进行足够数量的波束测深剖面叠加。由于中央波束受声速影响最小，因此先将测深横剖面上的波束深度标准化至中央波束深度。

②统计给定波束与中央的深度差、侧向中心距之比，并将其转化为角度(反正切关系)，统计各波束左右舷角度的均值即为标准化了测深横剖面的旋转改正角。

③用该旋转改正角对每一波束进行一次以中央波束为原点的旋转，并计算各波束新的深度和侧向中心距，完成波束归位。

由于沿岸海区海水表层易受多种因素影响，如太阳光照、海表面风、涌浪等，特别是岛屿附近，海水表层内易形成温度不一、密度不同的水团，使得该层的声速值变化更为复杂且无规律可循，利用常规的声速剖面进行声线跟踪不能反映声波的实际传播路径。因此通过结合三维激光扫描仪获得的水上下同目标特征点，利用几何改正法的思想对多波束测深仪中掠射角较低的波束进行归位，可以达到较好的结果，解决水上下无缝拼接问题。

3、实际实施与计算流程

该实施例采用Sonic2024多波束测深仪、VZ1000三维激光扫描仪和NovatelSpan-Lci组合导航系统集成的船载水岸线一体化测量系统(其中多波束测深仪中的多波束换能器偏转角度为30度)。

一种基于激光点云为约束的多波束低掠射角波束归位方法，采用船载水岸线一体化测量系统(如图2所示)，其水上部分包括三维激光扫描仪和组合导航系统，水下部分包括多波束测深仪，多波束测深仪带有具有一定偏转角度的多波束换能器，水上部分与水下部分通过一个平台支架固连；其中的三维激光扫描仪用于测量水面以上地物地貌，多波束测深仪用于测量水下地物地貌，组合导航系统用于为三维激光扫描仪和多波束测深仪提供定位信息、时间信息、姿态信息和航向信息；船载水岸线一体化测量系统还包括表层声速仪和声速剖面仪两个辅助设备；

所述的多波束低掠射角波束归位方法(如图1所示)按照如下步骤进行：

步骤1：对船载水岸线一体化测量系统进行标定，确定Sonic2024多波束测深仪、VZ1000三维激光扫描仪和组合导航系统之间的空间相对位置关系；

步骤2：选取一段具有水上下公共目标的测区，并利用船载水岸线一体化测量系统沿岸匀速进行测量，根据水质和地形情况实际情况，合理选择多波束测深仪的发射频率、ping率和波束开角，这里选择400KHz、30ping/s，波束开角为160度，获取水上下信息，并记录测量工作开始和结束的时间点；

步骤3：根据数据结果，对测量区域内水下多波束测深数据进行分区，分区间隔为10秒，即正常情况下每一分区内包含300ping数据；对存储的多波束数据按时间进行分区的具体方法如下：

步骤3.1：记录船载水岸线一体化测量系统开始工作和结束工作的时间，精确到秒；

步骤3.2：然后根据带有时间标签的每一个发射接收周期多波束信息进行工作区间上的分区，分区间隔为t秒；

步骤4：根据存储的表层声速仪和GPS数据存储的时间标签和每一个分区时间节点，获取每一个分区内行船附近的平均表层声速值c(A_ts)和船附近内一点的地理坐标(x_At,y_At)；一个分区内行船附近的平均表层声速值c(A_ts)和船附近一点的地理坐标(x_At,y_At)的具体获取方法如下：

步骤4.1：根据带有时间标签的每一个发射接收周期多波束信息和对应每一个分区时间节点，找到每一分区内记录的声速值c_i，根据公式(2)则可获得一个分区内行船附近的平均表层声速值c(A_ts)，即：

c(A_ts)＝(c₁+c₂+…+c_m)/m；(2)

步骤4.2：根据GPS记录数据的时间和对应每一个分区时间节点，获取分区时间中间节点处的坐标(b,l,h)，并转到地理坐标系下得到(x_At,y_At)；

步骤5：利用声速剖面仪以早中晚为三个时段，均匀获取沿岸测量区域内15个声速剖面数据svp_n，并记录每一个采样点的地理位置(x_N,y_N)，声速剖面站分布如图6所示；

步骤6：对所得到的15个声速剖面数据svp_n进行层化处理，获得不同位置点处同一水层上的声速c(z_i)和每一层声速的平均值对15个声速剖面进行层化处理，是由于在实际的采样数据中，受到实验条件的限制，采集的声速剖面数据不一定分布同一个深度上，一次通过对每一个声速剖面的采样点之间采用三次样条插值，从而得到等深点处的声速值，便于后续矩阵运算；

步骤7：通过步骤6层化处理后的声速剖面基于公式(1)得到测区内每一个点处的声速值：

$c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)---\left(1\right)$

其中加权系数α_k(x,y)为地理坐标，f_k(z_i)为正交函数，为每一层声速的平均值；

获得每一点处的声速c_n(x,y,z_i)和加权系数α_k(x,y)的具体方法如下：

步骤7.1：将根据步骤6得到的多个(例如15个)声速剖面数据写成声速值矩阵，形式如式(13)，再根据公式(15)求得协方差矩阵；

步骤7.2：求取声速值协方差矩阵的特征值和对应的特征向量，特征向量即为经验正交函数；

步骤7.3：将特征值按照从大到小排列，求取各个特征值的贡献率，选择较大的特征值所对应的特征向量，假设选定了P个特征值，并且相加后远远超过其他特征值，确定了1到P阶的经验正交函数，则式(19)变成了由P个未知数，K+1个线性方程组成的一个超定系统，利用正交三角分解的方法在最小二乘意义下求取α_k(x,y)，这样就可以利用求得的1到P阶的经验正交函数和X＝[a₀,a₁,a₂,a₃,]^T系数矩阵，表示测深区域任意一点的声速剖面；

步骤8：根据各分区内平均表层声速值c(A_ts)，对该区域内(x_At,y_At)处拟合的声速剖面表层部分进行修正；

对分区内(x_At,y_At)处拟合的声速剖面表层部分(多波束换能器吃水深度D)进行修正具体获得方法如下：

步骤8.1：根据公式(3)计算区域A_t内(x_At,y_At)处深度D(多波束换能器吃水深度)处声速修正量

$\begin{array}{c}{c}_n(x_{At},y_{At},D)=\stackrel{\&OverBar;}{c}\left(D\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x_{At},y_{At})f_k\left(D\right)\\ {\mathrm{\&Delta;c}}_t^0=c\left(A_{ts}\right)-c_n(x_{At},y_{At},D)\end{array}---\left(3\right);$

步骤8.2：根据以及公式(4)计算修正后的测区内每一点的声速值即：

$\begin{array}{c}{c}_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)+{\mathrm{\&Delta;c}}_t^0(z_i\&le;D)\\ c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)(z_i>D)\end{array}---\left(4\right);$

步骤9：选取水上下形状规则的目标区域，分别拾取20组约束处激光点和多波束测深点在多波束换能器坐标系下的坐标值，根据如图4所示的几何关系，利用公式(5)、(6)、(7)、(8)、(9)、(10)计算低掠射层固定层部分内波束初始几何旋转角和初始放大系数以及自由层部分各波束对应的旋转角和放大系数获得的具体方法如下：

步骤9.1：在区域A_t内选取水上下形状规则的目标，拾取多波束换能器坐标中心下K组约束激光点和多波束测深点坐标值(x_l,y_l,z_l)和(x_e,y_e,z_e)，根据几何关系，首先得到XOY平面内平均旋转角θ₁；

$\begin{array}{c}{\mathrm{\&theta;}}_e=\left(\underset{e=1}{\overset{K}{\&Sigma;}}\arctan \frac{y_e}{x_e}\right)/K\\ {\mathrm{\&theta;}}_l=\left(\underset{l=1}{\overset{K}{\&Sigma;}}\arctan \frac{y_l}{x_l}\right)/K\\ {\mathrm{\&theta;}}_1={\mathrm{\&theta;}}_l-{\mathrm{\&theta;}}_e\end{array}---\left(5\right);$

步骤9.2：根据几何关系，得到与XOY平面垂直方向上的平均旋转角θ₂：

$\begin{array}{c}{d}_l=\sqrt{x_l^2+y_l^2}\\ {\mathrm{\&theta;}}_2=\left(\underset{l=1}{\overset{K}{\&Sigma;}}\arctan \frac{z_l}{d_l}\right)/K\end{array}---\left(6\right)$

则初始几何旋转角即：

${\mathrm{\&theta;}}_{At}^0=\left(\begin{array}{c}{\mathrm{\&theta;}}_1\\ {\mathrm{\&theta;}}_2\end{array}\right)---\left(7\right);$

假设多波束换能器工作时开角为向右上方偏转角度为β，扫岸一侧某波束角为γ_i，则可以得到低掠射层内自由层部分各波束角对应的旋转角

${\mathrm{\&theta;}}_{At}^i=\left\{\begin{array}{c}{\mathrm{\&theta;}}_{At}^0(180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;}\&le;{\mathrm{\&gamma;}}_i\&le;\frac{\mathrm{\&alpha;}}{2})\\ \frac{{\mathrm{\&gamma;}}_i+\mathrm{\&beta;}-60}{120-\frac{\mathrm{\&alpha;}}{2}-\mathrm{\&beta;}}\&CenterDot;{\mathrm{\&theta;}}_{At}^0(60-\mathrm{\&beta;}<{\mathrm{\&gamma;}}_i<180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;})\end{array}\right.---\left(8\right);$

步骤9.3：根据两种点云的坐标，计算初始放大系数得：

同理可得低掠射层内自由层部分各波束角对应的放大系数

步骤10：对低掠射层固定层部分内的波束进行旋转角为放大系数为的几何改正，对自由层部分各波束，进行旋转角为放大系数为的几何改正，进行几何改正的具体方法如下：

步骤10.1：根据步骤9中获得的各波束放大系数，首先对低掠射层内波束点进行距离改正：

步骤10.2：根据步骤9中获得的各波束旋转角，对低掠射层内波束点进行角度改正，最终得到几何改正后的坐标(x'_e,y'_e,z'_e)：

$\left\{\begin{array}{c}{x}_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\cos (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1\left)\right)\\ y_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\sin (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1\left)\right)\\ z_e^{,}=s_e^i\&CenterDot;\sin \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\end{array}\right.---\left(12\right);$

利用公式(11)、(12)对该区域内300ping波束进行几何改正，其他区域同理可得。

图4为几何改正示意图。如图4所示，在多波束换能器坐标系中，A和B两点分别是三维激光扫描仪和多波束测深仪获得的水上下同目标特征点坐标，理想情况下A和B两点坐标应该一致，但由于海水表层声速变化复杂致使A和B坐标存在偏差。图中A'点为扫描仪激光点A点在XOY面上的投影点，B'点为多波束测深点在XOY面上的投影，θ_L为A'O与X轴之间的夹角，θ_e为B'O与X轴之间的夹角，θ₁为θ_L和θ_e之间的角度差，θ₂为AO与XOY面之间的夹角。

图5为低掠射层角度改正角度范围示意图。如图5所示，低掠射层由多波束换能器偏转方向一侧60度以上波束角组成，其包括固定层和自由层两部分，固定层部分由2倍多波束换能器最边缘波束与水平线之间夹角φ构成，自由层部分则由60度以上剩余部分角度组成。

图6为声速剖面站分布示意图。如图6所示，航行箭头方向即指测量船行进方向，测量箭头方向即一体化测量系统沿岸测量方向，声速剖面点分布要求均匀，具体数量不必参照图示，应根据测区实际大小合理选择并布设。

本发明一种基于激光点云为约束的多波束低掠射角波束归位方法，考虑到沿岸区域内海水声速的时空变化，并结合表层声速数据对拟合的声速剖面中浅层声速进行了修正；结合三维激光扫描仪与多波束测深仪测得的同目标点云对多波束60度以上的声波进行几何改正，纠正了多波束每一发射和接收中掠射角较低的波束的位置，解决了船载水岸线一体化测量系统中无缝拼接问题。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (7)

1.一种基于激光点云为约束的多波束低掠射角波束归位方法，采用船载水岸线一体化测量系统，其水上部分包括三维激光扫描仪和组合导航系统，水下部分包括多波束测深仪，多波束测深仪带有具有一定偏转角度的多波束换能器，水上部分与水下部分通过一个平台支架固连；其中的三维激光扫描仪用于测量水面以上地物地貌，多波束测深仪用于测量水下地物地貌，组合导航系统用于为三维激光扫描仪和多波束测深仪提供定位信息、时间信息、姿态信息和航向信息；船载水岸线一体化测量系统还包括表层声速仪和声速剖面仪两个辅助设备；

其特征在于：所述的多波束低掠射角波束归位方法按照如下步骤进行：

步骤1：对船载水岸线一体化测量系统进行标定，确定多波束测深仪、三维激光扫描仪和组合导航系统之间的空间相对位置关系；

步骤2：选取一段具有水上下公共目标的测区，并利用船载水岸线一体化测量系统沿岸匀速进行测量，获取水上下点云同一坐标系下坐标；

步骤3：对存储的多波束数据按时间分区；

步骤4：根据表层声速仪和组合导航系统中的GPS数据存储的时间标签，获取每一个分区内行船附近的平均表层声速值c(A_ts)和船附近内一点的地理坐标(x_At,y_At)；

步骤5：利用声速剖面仪均匀获取沿岸测量区域内多个声速剖面数据svp_n，并记录每一个采样点的地理位置(x_N,y_N)；

步骤6：对步骤5中得到的声速剖面数据svp_n进行层化处理，获得不同位置点处同一水层上的声速c(z_i)和每一层声速的平均值

步骤7：通过步骤6层化处理后的声速剖面基于公式(1)得到测区内每一个点处的声速值：

$c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)---\left(1\right)$

其中加权系数α_k(x,y)为地理坐标，f_k(z_i)为正交函数，为每一层声速的平均值；

步骤8：根据各分区内平均表层声速值c(A_ts)，对该区域内(x_At,y_At)处拟合的声速剖面表层部分进行修正；

步骤9：根据各区域修正后的声速剖面，进行声速改正，获取多波束测深仪中的多波束换能器坐标系下水上下同目标特征点处激光点与多波束点的平均坐标差，计算低掠射层内初始几何旋转角和初始放大系数进而按角度线性加权内插计算低掠射层内60度以上剩余波束角对应的旋转角和放大系数

步骤10：对低掠射层的固定层部分内的波束进行初始旋转角为和初始放大系数为的几何改正，对低掠射层的自由层部分各波束进行旋转角为放大系数为的几何改正。

2.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤3中，对存储的多波束数据按时间进行分区的具体方法如下：

步骤3.1：记录船载水岸线一体化测量系统开始工作和结束工作的时间，精确到秒；

步骤3.2：然后根据带有时间标签的每一个发射接收周期多波束信息进行工作区间上的分区，分区间隔为t秒。

3.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤4中，一个分区内行船附近的平均表层声速值c(A_ts)和船附近一点的地理坐标(x_At,y_At)的具体获取方法如下：

步骤4.1：根据带有时间标签的每一个发射接收周期多波束信息和对应的每一个分区的时间节点，找到每一分区内记录的声速值c_i，根据公式(2)则可获取每一个分区内行船附近的平均表层声速值c(A_ts)，即：

c(A_ts)＝(c₁+c₂+…+c_m)/m(2)；

步骤4.2：根据GPS记录数据的时间和对应每一个分区时间节点，获取分区时间中间节点处的坐标(b,l,h)，并转到地理坐标系下得到(x_At,y_At)。

4.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤7中，获得每一点处的声速c_n(x,y,z_i)和加权系数α_k(x,y)的具体方法如下：

步骤7.1：将根据步骤6得到的多个声速剖面数据写成声速值矩阵，然后再求协方差矩阵；

步骤7.2：求取声速值协方差矩阵的特征值和对应的特征向量；

步骤7.3：将特征值按照从大到小排列，求取各个特征值的贡献率，选择较大的特征值所对应的特征向量，利用正交三角分解的方法在最小二乘意义下求取加权系数α_k(x,y)，利用求得的经验正交函数和系数矩阵，表示测深区域任意一点的声速剖面。

5.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤8中，对分区内(x_At,y_At)处拟合的声速剖面表层部分即多波束换能器吃水深度为D处进行修正的具体方法如下：

步骤8.1：根据公式(3)计算区域A_t内(x_At,y_At)处多波束换能器吃水深度D处声速修正量

$\begin{array}{c}{c}_n(x_{At},y_{At},D)=\stackrel{\&OverBar;}{c}\left(D\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x_{At},y_{At})f_k\left(D\right)\\ {\mathrm{\&Delta;c}}_t^0=c\left(A_{ts}\right)-c_n(x_{At},y_{At},D)\end{array}---\left(3\right);$

步骤8.2：根据以及公式(4)计算修正后的测区内每一点的声速值即：

$c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)+{\mathrm{\&Delta;c}}_t^0(z_i\&le;D)$

$c_n(x,y,z)=\stackrel{\&OverBar;}{c}\left(z_i\right)+\underset{k=0}{\overset{k}{\&Sigma;}}{\mathrm{\&alpha;}}_k(x,y)f_k\left(z_i\right)(z_i>D)---\left(4\right).$

6.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤9中，低掠射层的定义如下：

低掠射层由多波束换能器偏转方向一侧60度以上波束角组成，其包括固定层和自由层两部分，固定层部分由2倍多波束换能器最边缘波束与水平线之间夹角φ构成，自由层部分则由60度以上剩余部分角度组成；

根据多波束换能器坐标系下水上下同目标特征点处激光点与多波束点的平均坐标差，计算低掠射层固定层部分内波束初始几何旋转角和初始放大系数以及自由层部分各波束对应的旋转角和放大系数获得的具体方法如下：

步骤9.1：在区域A_t内选取水上下形状规则的目标，拾取多波束换能器坐标中心下K组约束激光点和多波束测深点坐标值(x_l,y_l,z_l)和(x_e,y_e,z_e)，根据几何关系，首先得到XOY平面内平均旋转角θ₁；

${\mathrm{\&theta;}}_e=\left(\underset{e=1}{\overset{K}{\&Sigma;}}arctan\frac{y_e}{x_e}\right)/K$

${\mathrm{\&theta;}}_l=\left(\underset{l=1}{\overset{K}{\&Sigma;}}arctan\frac{y_l}{x_l}\right)/K$

θ₁＝θ_l-θ_e(5)；

步骤9.2：根据几何关系，得到与XOY平面垂直方向上的平均旋转角θ₂：

$d_l=\sqrt{x_l^2+y_l^2}$

${\mathrm{\&theta;}}_2=\left(\underset{l=1}{\overset{K}{\&Sigma;}}arctan\frac{z_l}{d_l}\right)/K---\left(6\right);$

则初始几何旋转角为：

${\mathrm{\&theta;}}_{At}^0=\left(\begin{array}{c}{\mathrm{\&theta;}}_1\\ {\mathrm{\&theta;}}_2\end{array}\right)---\left(7\right);$

假设多波束换能器工作时开角为向右上方偏转角度为β，扫岸一侧某波束角为γ_i，则可以得到低掠射层内自由层部分各波束对应的旋转角

${\mathrm{\&theta;}}_{At}^i=\left\{\begin{array}{c}{\mathrm{\&theta;}}_{At}^0(180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;}\&le;{\mathrm{\&gamma;}}_i\&le;\frac{\mathrm{\&alpha;}}{2})\\ \frac{{\mathrm{\&gamma;}}_i+\mathrm{\&beta;}-60}{120-\frac{\mathrm{\&alpha;}}{2}-\mathrm{\&beta;}}\&CenterDot;{\mathrm{\&theta;}}_{At}^0(60-\mathrm{\&beta;}<{\mathrm{\&gamma;}}_i<180-\frac{\mathrm{\&alpha;}}{2}-2\mathrm{\&beta;})\end{array}\right.---\left(8\right);$

步骤9.3：根据两种点云的坐标，计算初始放大系数得：

$S_l=\left(\underset{l=1}{\overset{20}{\&Sigma;}}\sqrt{x_l^2+y_l^2+z_l^2}\right)/K$

$S_e=\left(\underset{e=1}{\overset{20}{\&Sigma;}}\sqrt{x_e^2+y_e^2+z_e^2}\right)/K$

同理可得低掠射层内自由层部分各波束对应的放大系数

7.根据权利要求1所述的基于激光点云为约束的多波束低掠射角波束归位方法，其特征在于：在所述的步骤10中，进行几何改正的具体方法如下：

步骤10.1：根据步骤9中获得的各波束放大系数，首先对低掠射层内波束点进行距离改正：

步骤10.2：根据步骤9中获得的各波束旋转角，对低掠射层内波束点进行角度改正，最终得到几何改正后的坐标(x'_e,y'_e,z'_e)：

$\left\{\begin{array}{c}{x}_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\cos (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1))\\ y_e^{,}=s_e^i\&CenterDot;\cos \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\&CenterDot;\sin (\arctan \frac{y_e}{x_e}+{\mathrm{\&theta;}}_{At}^i({\mathrm{\&theta;}}_1))\\ z_e^{,}=s_e^i\&CenterDot;\sin \left({\mathrm{\&theta;}}_{At}^i\right({\mathrm{\&theta;}}_2\left)\right)\end{array}\right.---\left(12\right).$