CN114595420B - 基于多波束测深系统的沙波特征参数确定方法及装置 - Google Patents

Abstract

本发明提供基于多波束测深系统的沙波特征参数确定方法及装置，可在获取后处理数据时实时自动进行计算，直接提取沙波形态参数及确定运动方向，处理效率较高，结果准确可靠。方法包括：步骤1.利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据；步骤2.基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型；步骤3.求解每个中心单元格e的坡度和坡向；步骤4.确定背流面和迎流面、波峰和波谷所对应的单元格；步骤5.计算沙波的波高和平均背流角；步骤6.计算沙垄的波长、波高、背流角。

Description

基于多波束测深系统的沙波特征参数确定方法及装置

技术领域

本发明属于沙波参数获取技术领域，具体涉及基于多波束测深系统的沙波特征参数确定方法及装置。

技术背景

沙波运动是推移质泥沙输移的主要外在表现形式，也是河床变形的最常见形式之一，影响着床面阻力和水流泥沙运动，其几何形态参数的计算在行洪论证、取水工程、航道整治等生产实践活动中有着很大的影响，因此，观测沙波形态特征对于工程实践具有十分重要的意义。

天然河道沙波观测具有难度，对于天然河道沙波形态的研究本身较少。目前河道地形地貌测量方法以多波束测深(Multi-Beam Echo Sounding,MBES)系统为主，然而多波束测深系统经过后处理生成的数据和图像无法直接反映沙波的形态和运动方向，因此沙波形态的提取方法尤为重要。目前天然沙波形态的提取方法较为局限，常用的方法主要有：利用连续的河床地形影像或软件中可视化的测深数据，人工统计和整理沙波信息；利用遥感影像提取沙波信息；利用地质统计学和信号处理技术自动检测。但是这些方法大都存在精度较低，处理效率不高，不能真正反映床面沙波形态，只限于单一剖面而没有分析整个区域床面沙波的问题。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种基于多波束测深系统的沙波特征参数确定方法及装置，可在获取后处理数据时实时自动进行计算，直接提取沙波形态参数及确定运动方向，处理效率较高，结果准确可靠。

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

<方法>

本发明提供了一种基于多波束测深系统的沙波特征参数确定方法，其特征在于，包括以下步骤：

步骤1.利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据；

步骤2.基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型，包括如下子步骤：

步骤2-1.设置网格步长，找出原始测量数据中横坐标最小值、横坐标最大值、纵坐标最小值、纵坐标最大值的4个角点，依次记为A(x_a,y_a)、B(x_b,y_b)、C(x_c,y_c)、D(x_d,y_d)；

步骤2-2.根据4个角点求出沿流向的网格总长度、垂直流向的网格总宽度，将计算范围划分成n×m个网格，其中n、m向上取整数；

求解沿流向的AD段长度并划分n列的过程如下：

式中：n为河段沿流向被划分的列数；AD为河段沿流向的长度；dx为网格沿流向的步长；x_a、y_a为横坐标最小的角点A的横、纵坐标；x_d、y_d为纵坐标最大的角点D的横、纵坐标；

求解垂直流向的DB’段长度并划分m行的过程如下：

式中：m为河垂直流向被划分的行数；DB’为河段垂直流向的长度；dy为网格垂直流向的步长；x_b、y_b为横坐标最大的角点B的横、纵坐标；x、y为网格角点B’的横、纵坐标；α为河段沿流向与大地水平坐标系所成的角度；tanα＝(y_d-y_a)/(x_d-x_a)；

步骤2-3.设置网格方向和沙波运动方向一致；初始沙波运动方向假定为与流向一致，后续循环计算时沙波运动方向为与波峰点传播方向一致；求解各子单元格中心坐标{x(i,j),y(i,j)}：

式中：x(i,j)、y(i,j)为各子单元格中心坐标；x(1,1)、y(1,1)为第1行第1列的单元格中心点坐标；i、j为列数和行数；

步骤2-4.求解每个子单元格4条边界线的表达式，统计每个子单元格覆盖的数据点数p_ij，并将每个单元格内的数据高程光滑化，计算这一单元格内所有数据点高程的均值作为这一单元格中心点的高程z_ij；

步骤3.求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)；

步骤4.确定背流面和迎流面、波峰和波谷所对应的单元格；

步骤5.计算沙波的波高和平均背流角；

步骤6.计算沙垄的波长、波高、背流角。

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤1中，对待测区域采用多波束测深系统及其配套装置测量三维测深数据、水文数据、泥沙数据，并获取测船位置信息和姿态数据，然后对测量后的数据通过潮位误差消除、声速剖面误差消除、噪声滤波、坐标转换后处理，形成基于地理坐标系的三维离散地形数据点集{X_i,Y_i,Z_i}。

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤2-4中，单元格的4条边界线表达式分别为：

y＝tanα·x+a (9)

y＝tanα·x+b (10)

式中：x、y为边界线上点的坐标；截距a、b、c、d为未知数，代入单元格中心点坐标x(i,j)、y(i,j)可以得到单元格4条边界线中各截距a、b、c、d的值为：

统计每个子单元格覆盖的数据点数p_ij，若数据点(x₁,y₁)满足以下4个条件则在此单元格内，否则不在此单元格内：

y₁≥y＝tanα·x₁+a (17)

y₁<y＝tanα·x₁+b (18)

x₁≥x＝tanα·(c-y₁) (19)

x₁<x＝tanα·(d-y₁) (20)

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤3中，设定3×3的窗口，其中各单元格分别为a,b,c；d,e,f；g,h,i，单元格e为中心单元格，计算各中心单元格e的坡度S(i,j)和坡向C(i,j)：

中心单元格e的坡度S(i,j)为：

式中：S(i,j)为坡度；dz为高程变化量；dz/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

中心单元格e的坡向角度C(i,j)为：

式中：dz为高程变化量；dz/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤4中，根据坡向C(i,j)是否在0°～180°范围内，判断每个中心单元格e在背流面还是迎流面上；判断结束后，验证背流面上的数据点数与迎流面上的数据点数相加是否等于网格中各子单元格共覆盖的数据点数确定每一行的波峰、波谷数量和三维空间位置，其中波峰是从迎流面向背流面变化的这一单元格，波谷是从背流面向迎流面变化的这一单元格。

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤5中，计算各沙波波高、背流角：

h_s＝h_crest-h_trough (28)

式中：h_crest为某波峰高程；h_trough为下一个波谷高程；h_s为沙波波高；cresti为某波峰所在的单元格位置；troughi为下一个波谷所在的单元格位置；为对波峰到波谷的背流面所有单元格坡度求和；β为平均背流角。

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以具有以下特征：在步骤6中，分别计算大尺度沙垄和小尺度沙垄的波长、波高、背流角；根据沙波高度与沙波高度阈值的大小关系，判别沙波是大尺度沙垄还是小尺度沙垄；

小尺度沙垄满足：

其余为大尺度沙垄；

式中：h_s为某沙波高度，为沙波高度平均值，σ_s为沙波高度标准差；

判别后，重新计算大尺度沙垄和小尺度沙垄的波高和背流角；

根据上下两个波谷间的距离，分别计算大尺度沙垄和小尺度沙垄的波长：

式中：λ_l为大尺度沙垄波长，x_l2、y_l2为后一个大尺度沙垄波谷的横纵坐标，x_l1、y_l1为前一个大尺度沙垄波谷的横纵坐标，λ_s为小尺度沙垄波长，x_s2、y_s2为后一个小尺度沙垄波谷的横纵坐标，x_s1、y_s1为前一个小尺度沙垄波谷的横纵坐标。

优选地，本发明提供的基于多波束测深系统的沙波特征参数确定方法，还可以包括：

步骤8.数据验证与循环计算：

将计算得到的沙波形态参数值，与实测的沙波形态参数值作部分对比，用相对误差σ和相关系数ρ来衡量计算值和实测值两者的误差；若误差满足精度要求，则输出计算值；若不满足，则根据沙波波峰点的传播方向，判断沙波运动方向，然后进行第二轮计算，缩小网格尺度，但网格尺度应不小于地形数据平面精度，网格方向设置为和第一轮计算中的沙波运动方向一致，重复进行后续计算；如此循环直至计算值与实测值满足精度要求。这样可以更进一步提高精度，特别适合于对精度要求非常高的情况。

<装置>

进一步，本发明还提供了一种自动实现上述<方法>的基于多波束测深系统的沙波特征参数确定装置，其特征在于，包括：

地形数据点集形成部，利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据；

模型构件部，基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型，包括：

步骤2-1.设置网格步长，找出原始测量数据中横坐标最小值、横坐标最大值、纵坐标最小值、纵坐标最大值的4个角点，依次记为A(x_a,y_a)、B(x_b,y_b)、C(x_c,y_c)、D(x_d,y_d)；

步骤2-2.根据4个角点求出沿流向的网格总长度、垂直流向的网格总宽度，将计算范围划分成n×m个网格，其中n、m向上取整数；

求解沿流向的AD段长度并划分n列的过程如下：

式中：n为河段沿流向被划分的列数；AD为河段沿流向的长度；dx为网格沿流向的步长；x_a、y_a为横坐标最小的角点A的横、纵坐标；x_d、y_d为纵坐标最大的角点D的横、纵坐标；

求解垂直流向的DB’段长度并划分m行的过程如下：

式中：m为河垂直流向被划分的行数；DB’为河段垂直流向的长度；dy为网格垂直流向的步长；x_b、y_b为横坐标最大的角点B的横、纵坐标；x、y为网格角点B’的横、纵坐标；α为河段沿流向与大地水平坐标系所成的角度；tanα＝(y_d-y_a)/(x_d-x_a)；

步骤2-3.设置网格方向和沙波运动方向一致；初始沙波运动方向假定为与流向一致，后续循环计算时沙波运动方向为与波峰点传播方向一致；求解各子单元格中心坐标{x(i,j),y(i,j)}：

式中：x(i,j)、y(i,j)为各子单元格中心坐标；x(1,1)、y(1,1)为第1行第1列的单元格中心点坐标；i、j为列数和行数；

步骤2-4.求解每个子单元格4条边界线的表达式，统计每个子单元格覆盖的数据点数p_ij，并将每个单元格内的数据高程光滑化，计算这一单元格内所有数据点高程的均值作为这一单元格中心点的高程z_ij；

求解部，求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)；

单元格确定部，确定背流面和迎流面、波峰和波谷所对应的单元格；

计算部，计算沙波的波高和平均背流角以及沙垄的波长、波高、背流角；

控制部，与地形数据点集形成部、模型构件部、求解部、单元格确定部、计算部均通信相连，控制它们的运行。

发明的作用与效果

本发明提出的基于多波束测深系统的沙波特征参数确定方法及装置，在获取后处理数据后，实时自动进行计算，直接提取沙波形态参数及确定运动方向，实现沙波几何形态特征参数和运动方向的实时高效确定，从而实现对整个测量河道床面沙波形态的有效监测，处理过程简单且高效，得到的结果客观且可靠，对水利水电生产建设、大型河流沙波研究、航道整治等工程或生产实践具有重要意义。

附图说明

图1为本发明实施例涉及的基于多波束测深系统的沙波特征参数确定方法的流程图；

图2为本发明实施例涉及的三维离散地形数据点集可视化地形图；

图3为本发明实施例涉及的网格范围确定过程示意图；

图4为本发明实施例涉及的网格划分示意图；

图5为本发明实施例涉及的各单元格中心坐标值求解示意图；

图6为本发明实施例涉及的中心单元格e边界线示意图；

图7为本发明实施例涉及的双尺度沙垄示意图；

图8为本发明实施例涉及的沙波形态参数不同范围的占比分布直方图和累计分布曲线图，其中，(a)对应小尺度沙垄波长；(b)对应大尺度沙垄波长；(c)对应小尺度沙垄波高；(d)对应大尺度沙垄波高；(e)对应小尺度沙垄背流角；(f)对应大尺度沙垄背流角；

图9为本发明实施例涉及的深泓线典型纵剖面的计算高程与实测高程误差对比图；

图10为本发明实施例涉及的波峰点平面图。

具体实施方式

以下结合附图对本发明涉及的基于多波束测深系统的沙波特征参数确定方法及装置进行详细地说明。

<实施例>

如图1所示，本实施例所提供的基于多波束测深系统的沙波特征参数确定方法包括以下步骤：

步骤1：确定航行或试航行范围(待测量区域)。

测量前，可在测量范围内选取试航行区域，进行小范围多波束测深。若试航行区域不明显存在沙波，即区域内沙波尺度较小甚至几乎不存在，则可以避免耗费人力物力精力在此测量范围内进行测量。若明显存在沙波，再进而研究整个测量范围内天然河道沙波形态的特征。

试航行时，根据国内河道沙波波长的经验值，确定试航行的长度为3倍波长经验值，以保证试航行过程中，沿流向能够计算至少2个完整沙波形态、垂直流向能够连接至少2个波峰线。

航行方向为从靠左岸处顺流向开始不断往返测量。为避免测量时远波束测量不够准确，每次掉头往返的航迹与前一次航行航迹的距离为多波束测量宽度的1/4，从而保证前后两次航迹有一部分地形数据重叠。航行宽度为从左岸至右岸的河道主槽宽度。测船的航行速度不超过10km/h。

本实施例中河段总长约4km，总宽约1.5km，测船施测范围的长约4km，宽约500m，测量范围基本为河道主槽；每次测量至少完成一次往返。本实施例中所选工况对应的流量为34233m³/s。

步骤2：利用多波束测深系统及其配套装置对待测区域进行测量，经后处理后读入原始测量数据离散点集。

用多波束测深系统及其配套装置测量待测区域的三维测深数据、水文数据、泥沙数据，并获取测船位置信息和姿态数据(包括横摇参数、纵摇参数、艏摇参数、涌浪参数等)，测量后的数据通过潮位误差消除、声速剖面误差消除、噪声滤波、坐标转换等后处理，形成基于地理坐标系的三维离散地形数据点集{X_i,Y_i,Z_i}，各测次的后处理数据可视化结果如图2所示。将该三维离散地形数据点集输入程序进行后续的实时自动计算。

步骤3：构建基于三维离散地形数据点集的数字水深模型。

步骤3-1：设置网格步长，找出原始数据点集中横坐标最小值、横坐标最大值、纵坐标最小值、纵坐标最大值的角点，如图3中A(x_a,y_a)、B(x_b,y_b)、C(x_c,y_c)、D(x_d,y_d)所示。

步骤3-2：根据4个角点求出沿流向的网格总长度、垂直流向的网格总宽度，见图3中ADB'C'计算范围。如图4所示，将计算范围划分成n×m网格，其中n、m向上取整数。

求解沿流向的AD段长度并划分n列的过程如下：

式中：n为河段沿流向被划分的段数(即列数)；AD为河段沿流向的长度；dx为网格沿流向的步长(m)；x_a、y_a为横坐标最小的角点A的横、纵坐标；x_d、y_d为纵坐标最大的角点D的横、纵坐标。

求解垂直流向的DB’段长度并划分m行的过程如下：

式中：m为河垂直流向被划分的段数(即行数)；DB’为河段垂直流向的长度；dy为网格垂直流向的步长(m)；x_b、y_b为横坐标最大的角点B的横、纵坐标；x、y为网格角点B’的横、纵坐标；α为河段沿流向与大地水平坐标系所成的角度。其中，x-x_d、y-y_d是根据图3中虚线所示的3个三角形互为相似三角形得到的，tanα＝(y_d-y_a)/(x_d-x_a)。

步骤3-3：设置网格方向和沙波运动方向一致；初始沙波运动方向假定为与流向一致，后续循环计算时沙波运动方向为与波峰点传播方向一致；如图5所示，求解各子单元格中心坐标{x(i,j),y(i,j)}：

式中：x(i,j)、y(i,j)为各子单元格中心坐标；x(1,1)、y(1,1)为第1行第1列的单元格中心点坐标；i、j为列数和行数；其余各参数含义同前。

步骤304：各单元格内高程数据光滑化。求解每个子单元格4条边界线的表达式，统计每个子单元格覆盖的数据点数p_ij，并将每个单元格内的数据高程光滑化，计算这一单元格内所有数据点高程的均值作为这一单元格中心点的高程z_ij。

单元格的4条边界线见图6，表达式分别为：

y＝tanα·x+a (9)

y＝tanα·x+b (10)

式中：x、y为边界线上点的坐标，可以由单元格中心点坐标x(i,j)、y(i,j)转化；截距a、b、c、d为未知数，代入单元格中心点坐标x(i,j)、y(i,j)可以得到单元格4条边界线中各截距a、b、c、d的值为：

故可以统计各单元格内的数据点数，若数据点(x₁,y₁)满足以下4个条件则在此单元格内，否则不在此单元格内：

y₁≥y＝tanα·x₁+a (17)

y₁<y＝tanα·x₁+b (18)

x₁≥x＝tanα·(c-y₁) (19)

x₁<x＝tanα·(d-y₁) (20)

需要注意的是：

1、数据点数太大时，可以先把数据在i方向分段，统计各段内的数据点数；再把各段细分为单元格，统计各单元格内的数据点数。这样可以大量减少程序代码循环的次数，进一步提高处理的速度和效率。

2、若网格内部没有点，此时Z坐标值为NaN。当网格单元较大时，内部没有数据点的网格较少，可用周围的点进行线性插值；但是当网格单元较小、与数据测量精度相近时，用线性插值精确度会大大降低，插值不够平滑，此时可计算网格周围8个网格内的数据平均值作为此网格的数据值。

3、在计算时，最后一行(第m行)和最后一列(第n列)的单元格把压在下线、左线、上线、右线的数据计算在内；除此以外的每个单元格把压在下线、左线的数据计算在内，压在上线、右线的数据不计算在内，这样可以避免各单元格数据统计的重复。

步骤305：网格构建精度检验

网格划分应使网格区域尽量覆盖实测范围，从而使网格内实测数据点数更多，程序计算更准确；同时应尽量避免覆盖实测范围以外的空白范围，以避免空白范围高程插值的不准确。为检验网格构建的精度，在构建了基于三维离散地形数据点集的DDM后，输出网格中各子单元格共覆盖的数据点数计算其占原始数据点数的百分比/>占比在70％以上视为网格划分结果较合理；若占比在70％以下，则调整4个角点的值，使网格尽可能覆盖原始数据点。

步骤4：求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)。

设定3×3的窗口，其中各单元格分别为a,b,c；d,e,f；g,h,i，单元格e为中心单元格，计算各中心单元格e的坡度S(i,j)和坡向C(i,j)。

坡度S(i,j)是指坡面倾角的正切值，表示地形单元的陡缓程度；中心单元格e的坡度S(i,j)为：

式中：S(i,j)为坡度；dz为高程变化量；d_z/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

坡向C(i,j)是指坡面法线投影于水平面上的方位角，即坡度S(i,j)所面对的方向。中心单元格e的坡向值为：

坡向值根据下式转换为坡向角度C(i,j)(0°到360°)：

式中：C(i,j)为坡向角度；dz为高程变化量；d_z/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

步骤5：确定背流面和迎流面、波峰和波谷。

根据坡向C(i,j)是否在0°～180°范围内，判断每个e单元格在背流面还是迎流面上。判断结束后，验证背流面上的数据点数与迎流面上的数据点数相加是否等于网格中各子单元格共覆盖的数据点数

确定每一行的波峰、波谷数量和三维空间位置。其中波峰是从迎流面向背流面变化的这一单元格，波谷是从背流面向迎流面变化的这一单元格。

步骤6：计算各沙波波高、背流角。

h_s＝h_crest-h_trough (28)

式中：h_crest为某波峰高程；h_trough为下一个波谷高程；h_s为沙波波高，即波峰与下一个波谷的高度差；cresti为某波峰所在的单元格位置；troughi为下一个波谷所在的单元格位置；为对波峰到波谷的背流面所有单元格坡度求和；β为平均背流角，即沙波背流面上所有单元格的坡度平均值。

步骤7：区分大尺度沙垄和小尺度沙垄，分别计算大尺度沙垄和小尺度沙垄的波长、波高、背流角。

天然河流床面有沙纹和沙垄存在，将波高小于5cm的沙波划分为沙纹，波高大于等于5cm的沙波划分为沙垄。在河床表面十分发育的大尺度沙垄的迎、背水面上，常常叠加着一些尺度介于沙纹和大尺度沙垄之间的小尺度沙垄，形成双尺度沙垄，如图7所示。这种小尺度的沙波在河流中是很常见的，但大尺度沙垄对航道整治等工程实践影响更大。若不区分开，则小尺度沙垄的形态参数将影响大尺度沙垄的形态参数统计结果，导致形态参数计算偏小。故统计时应区分大尺度沙垄和小尺度沙垄。

根据沙波高度与沙波高度阈值的大小关系，判别沙波是大尺度沙垄还是小尺度沙垄。其中，小尺度沙垄满足：

其余为大尺度沙垄。式中：h_s为某沙波高度；为沙波高度平均值；σ_s为沙波高度标准差。

判别后，重新计算大尺度沙垄和小尺度沙垄的波高和背流角。

根据上下两个波谷间的距离，分别计算大尺度沙垄和小尺度沙垄的波长：

式中：λ_l为大尺度沙垄波长；x_l2、y_l2为后一个大尺度沙垄波谷的横纵坐标；x_l1、y_l1为前一个大尺度沙垄波谷的横纵坐标；λ_s为小尺度沙垄波长；x_s2、y_s2为后一个小尺度沙垄波谷的横纵坐标；x_s1、y_s1为前一个小尺度沙垄波谷的横纵坐标。

本实施例的工况下，统计的小尺度沙垄和大尺度沙垄的波长、波高、平均背流角等形态参数在不同范围的占比分布直方图和累计分布曲线见图8。根据直方图可以直观看出，小尺度沙垄的波长大多在6～21m范围内，占所有小尺度沙垄的70％，大尺度沙垄的波长大多在20～40m范围内，占所有大尺度沙垄的77％；小尺度沙垄的波高大多在0～0.45m范围内，占所有小尺度沙垄的43％，大尺度沙垄的波高大多在1.9～2.45m范围内，占所有大尺度沙垄的65％；小尺度沙垄的背流角在5°～6°范围内最多，占所有小尺度沙垄的15％，大尺度沙垄的背流角在6°～7°范围内最多，占所有大尺度沙垄的25％。

本实施例的工况下，最终计算得到小尺度沙垄和大尺度沙垄的平均波长分别为18.34m和34.62m，平均波高分别为0.74m和2.41m，与实测地形数据的抽样统计值和长江沙波尺度的经验值相符。

步骤8：数据验证与循环计算。

将编程统计得到的沙波形态参数值，与人工抽样统计的沙波形态参数值作部分对比。如针对某一典型纵剖面，对比编程输出和原始数据这一典型纵剖面的波峰和波谷。

本实施例中，输出深泓线所在典型纵剖面的计算高程与实测高程对比，见图9。用相对误差σ和相关系数ρ来衡量计算值和实测值两者的误差：

/>

式中：Z_ic为高程计算值；Z_it为高程真值；为高程计算值的均值；/>为高程真值的均值；σ_c为高程计算值标准差；σ_t为高程真值标准差；N为数据个数。

计算高程与实测高程的平均相对误差为7.09％，相关系数为0.80，相对误差达到精度要求，并且相关程度较高。

若要求更高的精度，则根据沙波波峰点的传播方向，判断沙波运动方向，然后进行第二轮计算。缩小网格尺度(网格尺度应不小于地形数据平面精度)，网格方向设置为和第一轮计算中的沙波运动方向一致，重复进行第一轮的后续计算。如此循环直至编程统计值与人工抽样统计的实测沙波形态参数值吻合度达到精度要求。

本实施例中，第一轮计算后，为了判断沙波的运动方向，把不同纵剖面的第1个波峰点依次相连、第2个波峰点依次相连……以此类推。每9个纵剖面输出1个纵剖面的波峰点连线，得到波峰点平面图，见图10。波峰点的传播方向如图10中虚线所示，由此可以判断本实施例中第一轮计算后，沙波的运动方向基本和水流流向平行，故无需修改初始沙波的运动方向、进行循环计算，可以直接输出程序的计算值。

本实施例是针对一次多波束测深系统测量的沙波形态参数及运动方向确定方法，应当指出，在较短时间间隔内若有多次多波束测深系统测量，还可在上述步骤进行完毕后，继续进行沙波运动速度的计算，这在研究沙波运动对于水流运动的响应、推移质输沙率的估算等方面具有重要作用。具体计算方法为：在较短时间间隔Δt内相邻两次测量后，输出根据本发明计算得到的同一典型纵剖面形态；计算2个纵剖面高程图中各沙波的波峰前移距离L_i，并进行加权平均，得到沙波前移的距离：

式中：n为纵剖面沙波个数；L_i为各沙波的波峰前移距离；ρ_i为加权系数，考虑沙波波高大小的差异对各前移距离，加权系数按如下公式计算：

式中：h_i为纵剖面各沙波的波高；为纵剖面各沙波波高之和。则沙波运动速度为沙波前移距离与相邻两次测量时间间隔之比，即v＝L/Δt。

进一步，本实施例还提供能够自动实现上述方法的基于多波束测深系统的沙波特征参数确定装置，该装置包括待测区域确定部、地形数据点集形成部、模型构件部、求解部、单元格确定部、计算部、验证与循环部、输入显示部以及控制部。

待测区域确定部执行上文步骤1所描述的内容，确定待测区域。

地形数据点集形成部执行上文步骤2所描述的内容，利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据。

模型构件部执行上文步骤3所描述的内容，基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型。

求解部执行上文步骤4所描述的内容，求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)。

单元格确定部执行上文步骤5所描述的内容，确定背流面和迎流面、波峰和波谷所对应的单元格。

计算部执行上文步骤6和7所描述的内容，计算沙波的波高和平均背流角以及沙垄的波长、波高、背流角。

验证与循环部执行上文步骤8所描述的内容，进行数据验证与循环计算。

输入显示部用于让用户输入操作指令，并根据操作指令对相应部的数据和文件以文字、表格、或图形方式进行显示。

控制部与待测区域确定部、地形数据点集形成部、模型构件部、求解部、单元格确定部、计算部、验证与循环部、输入显示部均通信相连，控制它们的运行。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的基于多波束测深系统的沙波特征参数确定方法及装置并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

Claims (10)

1.基于多波束测深系统的沙波特征参数确定方法，其特征在于，包括以下步骤：

步骤1.利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据；

步骤2.基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型，包括如下子步骤：

步骤2-1.设置网格步长，找出原始测量数据中横坐标最小值、横坐标最大值、纵坐标最小值、纵坐标最大值的4个角点，依次记为A(x_a,y_a)、B(x_b,y_b)、C(x_c,y_c)、D(x_d,y_d)；

步骤2-2.根据4个角点求出沿流向的网格总长度、垂直流向的网格总宽度，将计算范围划分成n×m个网格，其中n、m向上取整数；

求解沿流向的AD段长度并划分n列的过程如下：

式中：n为河段沿流向被划分的列数；AD为河段沿流向的长度；dx为网格沿流向的步长；x_a、y_a为横坐标最小的角点A的横、纵坐标；x_d、y_d为纵坐标最大的角点D的横、纵坐标；

求解垂直流向的DB’段长度并划分m行的过程如下：

式中：m为河垂直流向被划分的行数；DB’为河段垂直流向的长度；dy为网格垂直流向的步长；x_b、y_b为横坐标最大的角点B的横、纵坐标；x、y为网格角点B’的横、纵坐标；α为河段沿流向与大地水平坐标系所成的角度；tanα＝(y_d-y_a)/(x_d-x_a)；

步骤2-3.设置网格方向和沙波运动方向一致；初始沙波运动方向假定为与流向一致，后续循环计算时沙波运动方向为与波峰点传播方向一致；求解各子单元格中心坐标{x(i,j),y(i,j)}：

式中：x(i,j)、y(i,j)为各子单元格中心坐标；x(1,1)、y(1,1)为第1行第1列的单元格中心点坐标；i、j为列数和行数；

步骤2-4.求解每个子单元格4条边界线的表达式，统计每个子单元格覆盖的数据点数p_ij，并将每个单元格内的数据高程光滑化，计算这一单元格内所有数据点高程的均值作为这一单元格中心点的高程z_ij；

步骤3.求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)；

步骤4.确定背流面和迎流面、波峰和波谷所对应的单元格；

步骤5.计算沙波的波高和平均背流角；

步骤6.计算沙垄的波长、波高、背流角。

2.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：

其中，在步骤1中，对待测区域采用多波束测深系统及其配套装置测量三维测深数据、水文数据、泥沙数据，并获取测船位置信息和姿态数据，然后对测量后的数据通过潮位误差消除、声速剖面误差消除、噪声滤波、坐标转换后处理，形成基于地理坐标系的三维离散地形数据点集{X_i,Y_i,Z_i}。

3.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：其中，在步骤2-4中，单元格的4条边界线表达式分别为：

y＝tanα·x+a (9)

y＝tanα·x+b (10)

式中：x、y为边界线上点的坐标；截距a、b、c、d为未知数，代入单元格中心点坐标x(i,j)、y(i,j)可以得到单元格4条边界线中各截距a、b、c、d的值为：

统计每个子单元格覆盖的数据点数p_ij，若数据点(x₁,y₁)满足以下4个条件则在此单元格内，否则不在此单元格内：

y₁≥y＝tanα·x₁+a (17)

y₁<y＝tanα·x₁+b (18)

x₁≥x＝tanα·(c-y₁) (19)

x₁<x＝tanα·(d-y₁) (20)。

4.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：

其中，步骤2还包括：步骤2-5.网格构建合理性检验，在构建了基于三维离散地形数据点集的DDM后，输出网格中各子单元格共覆盖的数据点数计算其占原始数据点数的百分比/>占比在70％以上视为网格划分结果视为合理；若占比在70％以下，则调整4个角点的值，使网格尽可能覆盖原始数据点。

5.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：

其中，在步骤3中，设定3×3的窗口，其中各单元格分别为a,b,c；d,e,f；g,h,i，单元格e为中心单元格，计算各中心单元格e的坡度S(i,j)和坡向C(i,j)：

中心单元格e的坡度S(i,j)为：

式中：S(i,j)为坡度；dz为高程变化量；dz/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

中心单元格e的坡向角度C(i,j)为：

式中：dz为高程变化量；dz/dx为中心单元格e在x方向上的变化率；dz/dy为中心单元格e在y方向上的变化率：

6.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：

其中，在步骤4中，根据坡向C(i,j)是否在0°～180°范围内，判断每个中心单元格e在背流面还是迎流面上；判断结束后，验证背流面上的数据点数与迎流面上的数据点数相加是否等于网格中各子单元格共覆盖的数据点数

确定每一行的波峰、波谷数量和三维空间位置，其中波峰是从迎流面向背流面变化的这一单元格，波谷是从背流面向迎流面变化的这一单元格。

7.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：其中，在步骤5中，计算各沙波波高、背流角：

h_s＝h_crest-h_trough (28)

式中：h_crest为某波峰高程；h_trough为下一个波谷高程；h_s为沙波波高；cresti为某波峰所在的单元格位置；troughi为下一个波谷所在的单元格位置；为对波峰到波谷的背流面所有单元格坡度求和；β为平均背流角。

8.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于：

其中，在步骤6中，分别计算大尺度沙垄和小尺度沙垄的波长、波高、背流角；

根据沙波高度与沙波高度阈值的大小关系，判别沙波是大尺度沙垄还是小尺度沙垄；

小尺度沙垄满足：

其余为大尺度沙垄；

式中：h_s为某沙波高度，为沙波高度平均值，σ_s为沙波高度标准差；

判别后，重新计算大尺度沙垄和小尺度沙垄的波高和背流角；

根据上下两个波谷间的距离，分别计算大尺度沙垄和小尺度沙垄的波长：

式中：λ_l为大尺度沙垄波长，x_l2、y_l2为后一个大尺度沙垄波谷的横纵坐标，x_l1、y_l1为前一个大尺度沙垄波谷的横纵坐标，λ_s为小尺度沙垄波长，x_s2、y_s2为后一个小尺度沙垄波谷的横纵坐标，x_s1、y_s1为前一个小尺度沙垄波谷的横纵坐标。

9.根据权利要求1所述的基于多波束测深系统的沙波特征参数确定方法，其特征在于，还包括：

步骤8.数据验证与循环计算：

将计算得到的沙波形态参数值，与实测的沙波形态参数值作部分对比，用相对误差σ和相关系数ρ来衡量计算值和实测值两者的误差；若误差满足精度要求，则输出计算值；若不满足，则根据沙波波峰点的传播方向，判断沙波运动方向，然后进行第二轮计算，缩小网格尺度，但网格尺度应不小于地形数据平面精度，网格方向设置为和第一轮计算中的沙波运动方向一致，重复进行后续计算；如此循环直至计算值与实测值满足精度要求。

10.基于多波束测深系统的沙波特征参数确定装置，其特征在于，包括：

地形数据点集形成部，利用多波束测深系统对待测区域进行测量，然后对测量数据进行后处理，形成基于地理坐标系的三维离散地形数据点集作为输入模型的原始测量数据；

模型构件部，基于三维离散地形数据点集构建利于提取沙波特征参数的数字水深模型，包括：

步骤2-1.设置网格步长，找出原始测量数据中横坐标最小值、横坐标最大值、纵坐标最小值、纵坐标最大值的4个角点，依次记为A(x_a,y_a)、B(x_b,y_b)、C(x_c,y_c)、D(x_d,y_d)；

步骤2-2.根据4个角点求出沿流向的网格总长度、垂直流向的网格总宽度，将计算范围划分成n×m个网格，其中n、m向上取整数；

求解沿流向的AD段长度并划分n列的过程如下：

式中：n为河段沿流向被划分的列数；AD为河段沿流向的长度；dx为网格沿流向的步长；x_a、y_a为横坐标最小的角点A的横、纵坐标；x_d、y_d为纵坐标最大的角点D的横、纵坐标；

求解垂直流向的DB’段长度并划分m行的过程如下：

式中：m为河垂直流向被划分的行数；DB’为河段垂直流向的长度；dy为网格垂直流向的步长；x_b、y_b为横坐标最大的角点B的横、纵坐标；x、y为网格角点B’的横、纵坐标；α为河段沿流向与大地水平坐标系所成的角度；tanα＝(y_d-y_a)/(x_d-x_a)；

步骤2-3.设置网格方向和沙波运动方向一致；初始沙波运动方向假定为与流向一致，后续循环计算时沙波运动方向为与波峰点传播方向一致；求解各子单元格中心坐标{x(i,j),y(i,j)}：

式中：x(i,j)、y(i,j)为各子单元格中心坐标；x(1,1)、y(1,1)为第1行第1列的单元格中心点坐标；i、j为列数和行数；

步骤2-4.求解每个子单元格4条边界线的表达式，统计每个子单元格覆盖的数据点数p_ij，并将每个单元格内的数据高程光滑化，计算这一单元格内所有数据点高程的均值作为这一单元格中心点的高程z_ij；

求解部，求解每个中心单元格e的坡度S(i,j)和坡向C(i,j)；

单元格确定部，确定背流面和迎流面、波峰和波谷所对应的单元格；

计算部，计算沙波的波高和平均背流角以及沙垄的波长、波高、背流角；

控制部，与所述地形数据点集形成部、所述模型构件部、所述求解部、所述单元格确定部、所述计算部均通信相连，控制它们的运行。