CN103389077B - 一种基于mbes的海底沙波地貌运动探测方法 - Google Patents

Abstract

本发明公开海洋测绘、海底地形地貌调查研究的方法，具体是指一种基于MBES的海底沙波地貌运动探测方法。本发明是采用分米级测量精度的多波束测深技术和亚米级定位精度的导航定位系统，通过基于MBES的海底沙波探测前准备、基于MBES的海底沙波地貌的探测、基于MBES实测数据的DDM构建、经二次探测与处理、以及基于DDM的沙波运动速率计算得到海底沙波地貌的运动速率和方向。该发明给出了海底沙波运动的精细探测方案和海底网格的快速构建方法，进而给出了海底沙波运动速率的准确计算方法，在海底高活动区进行海洋测绘、海洋调查与海底科学研究均具有重要的实际应用价值。

Description

一种基于MBES的海底沙波地貌运动探测方法

技术领域

本发明涉及到海洋测绘、海底地形地貌调查研究和海底科学等技术领域，具体是指一种基于MBES的海底沙波地貌运动探测方法。

背景技术

多波束测深(MBES,Multi-Beam Echo Sounding)技术代表当代测深领域最新发展的高新技术成果，是多传感器复杂集成与多源数据高度融合的综合处理系统。多波束测深系统于20世纪60年代萌芽于美国海军需求，80年代开始商业化，90年代引进中国并大规模应用。多波束测深与传统的单波束测深相比具有高分辨率、高精度、高效率和全覆盖的显著特点。该技术目前广泛应用于海底水深测量与地形地貌调查研究，如海底沙波、沙脊地貌的探测等。

海底沙波是一种在潮流环境下广泛分布的海底形态，是一种活动性极强的海底地貌类型，一般而言，海底沙波的出现表征了海底底形的强烈活动性。沙波发育在海底床面上的特征形态，被认为是海底床面沉积物为适应水流剪切作用而产生变形、从而使床面发生规律变化的形态。我国邻近海域砂质底形十分发育，从渤海东岸浅滩、北黄海、苏北岸外、东海扬子滩，到台湾海峡和南海北部海区都有大面积的水下沙波分布。研究沙波的运动变化规律，不仅对沙质海岸地貌演变过程和沙质陆架的沉积演变具有重要意义，同时对于海底管线与石油平台等海洋工程建设的安全评估具有重要的应用价值。

国内外对于海底沙波的研究历史悠久，现有的沙波运动研究方法主要有：数值模型计算方法、遥感探测方法和海底定位观测方法等。数值计算方法精度较低，且不同的数值计算方法结果差异很大，遥感方法难以精确的定位且精度较低，在海底沙波年运动仅几米或十余米的海区，这两种方法显得无能为力。海底定位观测方法精度高，但缺点是难以实施且成本较高，仅适用在小范围的海区进行海底沙波运动探测，对于沙波在垂向变化观测比较适合，对于沙波水平运动观测较为困难。多波束测深是当前最先进的海底地形地貌探测方法，如何使用该技术进行海底沙波短周期运动的精确探测与计算显得尤为重要。从检索的公开资料分析，目前尚无与本发明完全一致的方法适用于海底沙波地貌的探测与计算问题。

发明内容

本发明针对海底沙波地貌运动难以准确探测与定量计算问题，提出基于MBES的海底沙波地貌运动探测与计算方法，从而实现海底沙波地貌运动的高效探测与运动速率的准确计算。

本发明是通过下述技术方案得以实现的：

基于MBES的海底沙波地貌运动探测与计算方法，包括下列步骤：

步骤1：基于MBES的海底沙波探测前准备

（1）设备检测与标定：多波束测深系统是一种复杂的海底地形地貌探测系统，在测量之前所有仪器设备均需按照国家标准GB/T12763.10-2007进行自校或送法定计量机构进行标定，以保证所有传感器处于正常工作状态，保证仪器精度符合国标要求。多波束系统的测深精度达到分米级甚至厘米级，GPS定位精度达到亚米级或分米级。

（2）潮位站布设：围绕沙波测量区布设2～4个临时验潮站，获取临时潮位资料Tide＝{tide_i}_i＝2,4，以控制测区潮位。

（3）声速剖面测量：开始勘测前，必须在测区附近进行声速剖面测量,至少获取测区一个声速剖面

n1为声速剖面层数。

步骤2：基于MBES的海底沙波地貌的探测

可采取十字正交测线或矩形区域两种布线方式进行海底沙波的多波束探测，具体步骤如下：

(1)十字测线探测法：在沙波区布设十字正交测线，其中一条平行沙波走向，一条垂直沙波走向。测线的宽度D应大于沙波的波长d，测线的长度L应大于沙波的横向延伸距离l。勘测时要求直线航行，船速6节左右，偏航小于5m,多波束全开角勘测，全程GPS有差分信号。

或(2)矩形区域探测法：在沙波区布设矩形区域，矩形的长边应平行沙波走向，短边垂直沙波走向，矩形的长边A应大于沙波的横向延伸距离l，矩形的短边B应大于沙波波长d的两倍。要求直线航行，船速6节左右，偏航小于5m,多波束全开角勘测，全程GPS有差分信号。

形成首次探测原始数据集合Raw_t1＝{raw1_i}，t1为海底地形测量时间。

步骤3：基于MBES实测数据的DDM(Digital Depth Model，简称网格)构建

(1)数据处理：对获取的原始数据集合Raw_t1＝{raw1_i}进行潮位改正、吃水改正、声速改正和噪声点编辑等处理后，形成处理后的离散水深数据集合Proc_t1＝{(x_i,y_i,z_i)}_i＝1,n,处理后的数据要求尽量保留所有可用离散波束点。

(2)数据分辨率评估：将离散水深数据集合正投影到二维平面上，按抽样法量算离散数据点的空间距离，要求95％的离散水深数据与其周围点的距离小于沙波运动估计速率v，否则降低船速，返回步骤2重新进行探测。

(3)构建DDM:采取改进的距离反比加权(FIDW,Fast Inverse DistanceWeighted)方法对离散水深数据集合Proc_t1进行处理，形成DDM_t1＝{dep_(i,j)}_{i＝1,n,j＝1,m}。

FIDW方法计算公式为：

${\mathrm{dep}}_{(i,j)}=\left[\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k{z}_k\right]/\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k;$

w_k＝1/d_k ²;

$d_k=\sqrt{{(x_{(i,j)}-x_k)}^2+{\left(y_{(i,j)}\right)-y_k}^2}$

式中，x_k，y_k，z_k为离散水深点的横坐标、纵坐标和水深值，来自集合Proc_t1。w_k为离散水深数据点的权重值。dep_(i,j)为网格值，x_(i,j)和y_(i,j)为网格横坐标和纵坐标值。

设计网格点结构体DEP＝{dep_(i,j),A_(i,j),B_(i,j)},

在读入离散数据的同时，根据每个离散水深数据的坐标（x_k，y_k），根据预先设定的作用距离dis（dis由用户设定，一般为网格间距的2～3倍），判断该点可作用的网格点，也就是说当d_k≤d时，该点参与网格化运算，否则直接跳过。当所有离散点遍历完成后，遍历一次网格即可获取每个网格点的水深差值dep_(i,j)＝A_(i,j)/B_(i,j)，当B_(i,j)＝0时该点为无效点。当参与网格dep_(i,j)运算的点数少于n2时，该点也为无效点，n2为自然数，由系统缺省给定，也可由用户后期交互修改；

步骤4：二次探测与处理

(1)二次探测：按照步骤2中（1）或（2）的测量方法，采取一定的时间间隔t进行二次重复测量，前后重复测量的航迹保持一致。测量时间间隔t要求为：t≥Δd÷v，Δd为Gps定位中误差，v为沙波运动估计速度。或者要求在一次风暴潮发生前后分别测量两次海底地形，用于研究风暴引起的海底沙波运动。重复探测形成原始数据集合Rwa_t2＝{raw2_i}，t2为二次海底地形测量时间。

(2)二次构建DDM：按照步骤3中（1）和（2）之方法，采取FIDW方法构建二次探测后的DDM_t2。

DDM_t1和DDM_t2为沙波两次测量形成的网格，两个网格的范围、行数、列数、行间距和列间距要完全一致。

(3)构建差值ΔDDM：DDM_t1和DDM_t2按网格点逐点进行相减运算，生成新

的差值网格ΔDDM＝DDM_t2-DDM_t1。

步骤5：基于DDM的沙波运动速率计算

采用剖面对比法或差值法两种方法来判断海底沙波是否运动，并计算海底沙波运动速率。

（1）剖面对比法：垂直沙波走向，设计地形剖面线L₀，分别在网格DDM_t1、DDM_t2和ΔDDM中进行相交运算，获取3条地形剖面线L_t1、L_t2和L_Δt，在同一平面坐标系下按照相同参数用不同颜色绘制三条地形剖面线，对比沙波峰或谷的迁移距离Δl₁,可以获取一处沙波的运动速率v1_i＝Δl₁÷t；按照上述方法，垂直一条沙波的不同位置，可获取该沙波n处位置的运动速率，通过平均的方法可以获取该沙波的平均运动速率

或（2）差值法：基于期沙波差值网格ΔDDM，绘制海底三维地形图，如果沙波出现运动，海底地形表现为谷峰相间条纹，通过人机交互模式追踪条纹的脊线和谷底线，测量距离最近的脊线和谷底线空间距离Δl₂，可以获取一处沙波的运动速率v2_i＝Δl₂÷t；按照上述方法，可获取n条沙波的运动速率，通过平均的方法可以获取一片沙波的平均运动速率

或(3子窗相关法：通过建立滑动子窗口的方式来分割两期海底网格的互相关来判断海底沙波的运动。在两期海底网格DDM_t1和DDM_t2中，建立相同大小的移动的矩形子窗口ΔD_t1和ΔD_t2，其行数和列数分别为m和n（由系统给定初始值，但可由用户自动调整大小），通过子窗口的移动来计算两个时期沙波网格的相关系数，从而判断海底沙波的运动方向和速率。

具体步骤如下：

（a）两个子窗口的互相关计算公式为：

$R=\frac{R_{\mathrm{ab}}}{R_a{R}_b}$

其中，R_a、R_b分别为子窗口ΔD_t1和ΔD_t2的方差，R_ab为子窗口ΔD_t1和ΔD_t2的协方差，计算公式分别为：

$R_a=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(a_{\mathrm{ij}}-\stackrel{\‾}{a})}^2}$ $R_b=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(b_{\mathrm{ij}}-\stackrel{\‾}{b})}^2}$

$R_{\mathrm{ab}}=\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(a_{\mathrm{ij}}-\stackrel{\‾}{a})(b_{\mathrm{ij}}-\stackrel{\‾}{b})$

a_ij、b_ij为两个子窗口ΔD_t1和ΔD_t2相应位置的水深值，为子窗口ΔD_t1和ΔD_t2的平均水深值。

(b)按照空间顺序在DDM_t1和DDM_t2中分别设置初始的子窗口ΔD_t1和ΔD_t2。

(c)固定ΔD_t1，在DDM_t2中按先行后列的空间顺序滑动子窗口ΔD_t2，按照步骤（a）中之公式分别计算子窗口ΔD_t1与移动子窗口ΔD_t2的互相关系数R，保存最大相关系数所在的子窗口ΔD_t2。通过两个子窗口ΔD_t1和ΔD_t2的中心坐标值可获取两者的距离dis₁₂和方向关系ang₁₂。子窗口ΔD_t1所在的沙波在t1至t2时间内的移动速率：v3＝dis₁₂÷t。

(d)在DDM_t1中按先行后列的空间顺序滑动子窗口ΔD_t1，然后按照步骤（c）可获取每个滑动子窗口ΔD_t1的最大移动速率v3_(i,j)和移动方向ang_(i,j)。

有益效果：

本发明的显著特点是基于多波束测深技术，采用亚米级精度的导航定位系统，实现了海底沙波地形地貌的精确探测，为揭示沙波的运动规律，提出了十字交叉线和矩形区域的方式进行重复探测，在此基础上，采用FIDW方法进行多波束测深数据的网格化，并基于两期DDM及差值，采用剖面法、差值法和滑动子窗相关法计算海底沙波地貌的运动速率。该发明不仅解决了海底沙波运动的精细探测方案，还给出了海底沙波运动速率的准确计算方法。该发明在海洋测绘、海洋调查与海底科学研究中具有重要的实际应用价值。

附图说明

图1本发明的工作流程图

图2本发明实施的沙波探测方案

图3本发明实施的剖面法计算沙波运动速率

图4本发明实施的差值法计算沙波运动速率

图5本发明实施的滑动子窗相关法计算沙波运动速率与方向

具体实施方式

下面结合本发明的实施作具体说明：

实施例1

基于MBES的海底沙波地貌运动探测与计算方法是按照下述步骤实现的。

本发明详细步骤与流程见图1。

步骤1：基于MBES的海底沙波探测前准备

（1）设备检测与标定：多波束测深系统是一种复杂的海底地形地貌探测系统，在测量之前所有仪器设备均需按照国家标准GB/T12763.10-2007进行自校或送法定计量机构进行标定，以保证所有传感器处于正常工作状态，保证仪器精度符合国标要求。多波束系统的测深精度达到分米级甚至厘米级，GPS定位精度达到亚米级或分米级。

（2）潮位站布设：围绕沙波测量区布设2个临时验潮站，获取临时潮位资料Tide＝{tide_i}_i＝1,4，以控制测区潮位。

（3）声速剖面测量：开始勘测前，必须在测区附近进行声速剖面测量,至少获取测区一个声速剖面

n1为声速剖面层数。

步骤2：基于MBES的海底沙波地貌的探测

可采取十字正交测线或矩形区域两种布线方式进行海底沙波的多波束探测，具体步骤如下：

(1)十字测线探测法：在沙波区布设十字正交测线，其中一条平行沙波走向，一条垂直沙波走向。测线的宽度D大于沙波的波长d，测线的长度L大于沙波的横向延伸距离l。勘测时要求直线航行，船速6节左右，偏航小于5m,多波束全开角勘测，全程GPS有差分信号。

或(2)矩形区域探测法：在沙波区布设矩形区域，矩形的长边应平行沙波走向，短边垂直沙波走向，矩形的长边A应大于沙波的横向延伸距离l，矩形的短边B应大于沙波波长d的两倍。要求直线航行，船速6节左右，偏航小于5m,多波束全开角勘测，全程GPS有差分信号。

形成首次探测原始数据集合Raw_t1＝{raw1_i}，t1为海底地形测量时间。

沙波探测方案见图2。

步骤3：基于MBES实测数据的DDM(Digital Depth Model，简称网格)构建

(4)数据处理：对获取的原始数据集合Raw_t1＝{raw1_i}进行潮位改正、吃水改正、声速改正和噪声点编辑等处理后，形成处理后的离散水深数据集合Proc_t1＝{(x_i,y_i,z_i)}_i＝1,n,处理后的数据要求尽量保留所有可用离散波束点。

(5)数据分辨率评估：将离散水深数据集合正投影到二维平面上，按抽样法量算离散数据点的空间距离，要求95％的离散水深数据与其周围点的距离小于沙波运动估计速率v，否则降低船速，返回步骤2重新进行探测。

(6)构建DDM:采取改进的距离反比加权(FIDW,Fast Inverse DistanceWeighted)方法对离散水深数据集合Proc_t1进行处理，形成DDM_t1＝{dep_(i,j)}_{i＝1,n,j＝1,m}。

FIDW方法计算公式为：

${\mathrm{dep}}_{(i,j)}=\left[\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k{z}_k\right]/\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k;$

w_k＝1/d_k ²;

$d_k=\sqrt{{(x_{(i,j)}-x_k)}^2+{(y_{(i,j)}-y_k)}^2}$

式中，x_k，y_k，z_k为离散水深点的横坐标、纵坐标和水深值，来自集合Proc_t1。w_k为离散水深数据点的权重值。dep_(i,j)为网格值，x_(i,j)和y_(i,j)为网格横坐标和纵坐标值。

设计网格点结构体DEP＝{dep_(i,j),A_(i,j),B_(i,j)},

在读入离散数据的同时，根据每个离散水深数据的坐标（x_k，y_k），根据预先设定的作用距离dis（dis由用户设定，一般为网格间距的2～3倍），判断该点可作用的网格点，也就是说当d_k≤d时，该点参与网格化运算，否则直接跳过。当所有离散点遍历完成后，遍历一次网格即可获取每个网格点的水深差值dep_(i,j)＝A_(i,j)/B_(i,j)，当B_(i,j)＝0时该点为无效点；当参与网格dep_(i,j)运算的点数少于n2时，该点也为无效点，n2为自然数，由系统缺省给定，也可由用户后期交互修改。

步骤4：二次探测与处理

(1)二次探测：按照步骤2中（1）或（2）的测量方法，采取一定的时间间隔t进行二次重复测量，前后重复测量的航迹保持一致。测量时间间隔t要求为：t≥Δd÷v，Δd为Gps定位中误差，v为沙波运动估计速度。或者要求在一次风暴潮发生前后分别测量两次海底地形，用于研究风暴引起的海底沙波运动。重复探测形成原始数据集合Rwa_t2＝{raw2_i}，t2为二次海底地形测量时间。

(2)二次构建DDM：按照步骤3中（1）和（2）之方法，采取FIDW方法构建二次探测后的DDM_t2。

DDM_t1和DDM_t2为沙波两次测量形成的网格，两个网格的范围、行数、列数、行间距和列间距要完全一致。

(3)构建差值ΔDDM：DDM_t1和DDM_t2按网格点逐点进行相减运算，生成新

的差值网格ΔDDM＝DDM_t2-DDM_t1。

步骤5：基于DDM的沙波运动速率计算

采用剖面对比法或差值法两种方法来判断海底沙波是否运动，并计算海底沙波运动速率。

（1）剖面对比法：垂直沙波走向，设计地形剖面线L₀，分别在网格DDM_t1、DDM_t2和ΔDDM中进行相交运算，获取3条地形剖面线L_t1、L_t2和L_Δt，在同一平面坐标系下按照相同参数用不同颜色绘制三条地形剖面线，对比沙波峰或谷的迁移距离Δl₁,可以获取一处沙波的运动速率v1_i＝Δl₁÷t；按照上述方法，垂直一条沙波的不同位置，可获取该沙波n处位置的运动速率，通过平均的方法可以获取该沙波的平均运动速率

剖面对比法见图3。

或（2）差值法：基于期沙波差值网格ΔDDM，绘制海底三维地形图，如果沙波出现运动，海底地形表现为谷峰相间条纹，通过人机交互模式追踪条纹的脊线和谷底线，测量距离最近的脊线和谷底线空间距离Δl₂，可以获取一处沙波的运动速率v2_i＝Δl₂÷t；按照上述方法，可获取n条沙波的运动速率，通过平均的方法可以获取一片沙波的平均运动速率

网格差值法见图4。

或(3子窗相关法：通过建立滑动子窗口的方式来分割两期海底网格的互相关来判断海底沙波的运动。在两期海底网格DDM_t1和DDM_t2中，建立相同大小的移动的矩形子窗口ΔD_t1和ΔD_t2，其行数和列数分别为m和n（由系统给定初始值，但可由用户自动调整大小），通过子窗口的移动来计算两个时期沙波网格的相关系数，从而判断海底沙波的运动方向和速率。

具体步骤如下：

（a）两个子窗口的互相关计算公式为：

$R=\frac{R_{\mathrm{ab}}}{R_a{R}_b}$

其中，R_a、R_b分别为子窗口ΔD_t1和ΔD_t2的方差，R_ab为子窗口ΔD_t1和ΔD_t2的协方差，计算公式分别为：

$R_a=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(a_{\mathrm{ij}}-\stackrel{\‾}{a})}^2}$

$R_b=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(b_{\mathrm{ij}}-\stackrel{\‾}{b})}^2}$

$R_{\mathrm{ab}}=\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(a_{\mathrm{ij}}-\stackrel{\‾}{a})(b_{\mathrm{ij}}-\stackrel{\‾}{b})$

a_ij、b_ij为两个子窗口ΔD_t1和ΔD_t2相应位置的水深值，为子窗口ΔD_t1和ΔD_t2的平均水深值。

(b)按照空间顺序在DDM_t1和DDM_t2中分别设置初始的子窗口ΔD_t1和ΔD_t2。

(c)固定ΔD_t1，在DDM_t2中按先行后列的空间顺序滑动子窗口ΔD_t2，按照步骤（a）中之公式分别计算子窗口ΔD_t1与移动子窗口ΔD_t2的互相关系数R，保存最大相关系数所在的子窗口ΔD_t2。通过两个子窗口ΔD_t1和ΔD_t2的中心坐标值可获取两者的距离dis₁₂和方向关系ang₁₂。子窗口ΔD_t1所在的沙波在t1至t2时间内的移动速率：v3＝dis₁₂÷t。

(d)在DDM_t1中按先行后列的空间顺序滑动子窗口ΔD_t1，然后按照步骤（c）可获取每个滑动子窗口ΔD_t1的最大移动速率v3_(i,j)和移动方向ang_(i,j)。

子窗口相关法见图5。

Claims (1)

1.一种基于MBES的海底沙波地貌运动探测方法，其特征在于包括下列步骤：

步骤1：基于MBES的海底沙波探测前准备

（1）设备检测与标定：选用多波束系统的测深精度达到分米级的仪器，GPS定位精度达到亚米级；

（2）潮位站布设：在海底沙波测量区布设2～4个临时验潮站，获取临时潮位资料Tide＝{tide_i}_i＝2,4，以控制测区潮位；

（3）声速剖面测量：开始勘测前，必须在测区附近进行声速剖面测量,至少获取测区一个声速剖面

n1为声速剖面层数；

步骤2：基于MBES的海底沙波地貌的探测

采取十字正交测线探测法或矩形区域探测法两种布线方式进行海底沙波的多波束探测，

其中，十字正交测线探测法是指：在沙波区布设十字正交测线，其中一条平行沙波走向，一条垂直沙波走向；测线的宽度D应大于沙波的波长d，测线的长度L应大于沙波的横向延伸距离l；勘测时沿直线移动，多波束全开角勘测，全程GPS有差分信号；

矩形区域探测法是：在沙波区布设矩形区域，矩形的长边应平行沙波走向，短边垂直沙波走向，矩形的长边A应大于沙波的横向延伸距离l，矩形的短边B应大于沙波波长d的两倍；勘测时沿直线移动，多波束全开角勘测，全程GPS有差分信号；

形成首次探测原始数据集合Raw_t1＝{raw1_i}，t1为海底地形测量时间；

步骤3：基于MBES实测数据的DDM构建

(1)数据处理：对获取的原始数据集合Raw_t1＝{raw1_i}进行潮位改正、吃水改正、声速改正和噪声点编辑处理后，形成处理后的离散水深数据集合Proc_t1＝{(x_i,y_i,z_i)}_i＝1,n；

(2)数据分辨率评估：将离散水深数据集合正投影到二维平面上，按抽样法量算离散数据点的空间距离，要求95％的离散水深数据与其周围点的距离小于沙波运动估计速率v，否则降低直线移动速度，返回步骤2重新进行探测；

(3)构建DDM:采取改进的距离反比加权(FIDW,Fast Inverse DistanceWeighted)方法对离散水深数据集合Proc_t1进行处理，形成DDM_t1＝{dep_(i,j)}_{i＝1,n,j＝1,m}；

FIDW方法计算公式为：

${\mathrm{dep}}_{(i,j)}=\left[\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k{z}_k\right]/\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{w}_k;$

w_k＝1/d_k ²;

$d_k=\sqrt{{(x_{(i,j)}-x_k)}^2+{(y_{(i,j)}-y_k)}^2}$

式中，x_k，y_k，z_k为离散水深点的横坐标、纵坐标和水深值，来自集合Proc_t1；w_k为离散水深数据点的权重值；dep_(i,j)为网格值，x_(i,j)和y_(i,j)为网格横坐标和纵坐标值；

设计网格点结构体DEP＝{dep_(i,j),A_(i,j),B_(i,j)},

在读入离散数据的同时，根据每个离散水深数据的坐标（x_k，y_k），根据预先设定的作用距离dis，判断该点可作用的网格点；若d_k≤dis时，该点参与网格化运算，否则直接跳过；当所有离散点遍历完成后，遍历一次网格即可获取每个网格点的水深差值dep_(i,j)＝A_(i,j)/B_(i,j)，当B_(i,j)＝0时该点为无效点；当参与网格dep_(i,j)运算的点数少于n2时，该点也为无效点，n2为自然数，由系统缺省给定，也可由用户后期交互修改；

步骤4：二次探测与处理

(1)二次探测：按照步骤2中的测量方法，采取一定的时间间隔t进行二次重复测量，前后重复测量的航迹保持一致；测量时间间隔t要求为：t≥Δd÷v，Δd为Gps定位中误差，v为沙波运动估计速度；或者要求在一次风暴潮发生前后分别测量两次海底地形，用于研究风暴引起的海底沙波运动；重复探测形成原始数据集合Rwa_t2＝{raw2_i}，t2为二次海底地形测量时间；

(2)二次构建DDM：按照步骤3中（1）和（2）之方法，采取FIDW方法构建二次探测后的DDM_t2；

DDM_t1和DDM_t2为沙波两次测量形成的网格，两个网格的范围、行数、列数、行间距和列间距要完全一致；

(3)构建差值ΔDDM：DDM_t1和DDM_t2按网格点逐点进行相减运算，生成新的差值网格ΔDDM＝DDM_t2-DDM_t1；

步骤5：基于DDM的沙波运动速率计算

采用剖面对比法、差值法、或子窗相关法的方法来判断海底沙波是否运动，并计算海底沙波运动速率；

其中，剖面对比法是指：垂直沙波走向，设计地形剖面线L₀，分别在网格DDM_t1、DDM_t2和ΔDDM中进行相交运算，获取3条地形剖面线L_t1、L_t2和L_Δt，在同一平面坐标系下按照相同参数用不同颜色绘制三条地形剖面线，对比沙波峰或谷的迁移距离Δl₁,可以获取一处沙波的运动速率v1_i＝Δl₁÷t；按照上述方法，垂直一条沙波的不同位置，可获取该沙波n处位置的运动速率，通过平均的方法可以获取该沙波的平均运动速率

差值法是指：基于期沙波差值网格ΔDDM，绘制海底三维地形图，如果沙波出现运动，海底地形表现为谷峰相间条纹，通过人机交互模式追踪条纹的脊线和谷底线，测量距离最近的脊线和谷底线空间距离Δl₂，可以获取一处沙波的运动速率v2_i＝Δl₂÷t；按照上述方法，可获取n条沙波的运动速率，通过平均的方法可以获取一片沙波的平均运动速率

子窗相关法是指：通过建立滑动子窗口的方式来分割两期海底网格的互相关来判断海底沙波的运动，在两期海底网格DDM_t1和DDM_t2中，建立相同大小的移动的矩形子窗口ΔD_t1和ΔD_t2，其行数和列数分别为m和n，通过子窗口的移动来计算两个时期沙波网格的相关系数，从而判断海底沙波的运动方向和速率；

具体步骤如下：

（a）两个子窗口的互相关计算公式为：

$R=\frac{R_{\mathrm{ab}}}{R_a{R}_b}$

其中，R_a、R_b分别为子窗口ΔD_t1和ΔD_t2的方差，R_ab为子窗口ΔD_t1和ΔD_t2的协方差，计算公式分别为：

$R_a=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(a_{\mathrm{ij}}-\stackrel{\‾}{a})}^2}$

$R_b=\sqrt{\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(b_{\mathrm{ij}}-\stackrel{\‾}{b})}^2}$

$R_{\mathrm{ab}}=\frac{1}{m-1}\frac{1}{n-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(a_{\mathrm{ij}}-\stackrel{\‾}{a})(b_{\mathrm{ij}}-\stackrel{\‾}{b})$

a_ij、b_ij为两个子窗口ΔD_t1和ΔD_t2相应位置的水深值，

为子窗口ΔD_t1和ΔD_t2的平均水深值；

(b)按照空间顺序在DDM_t1和DDM_t2中分别设置初始的子窗口ΔD_t1和ΔD_t2；

(c)固定ΔD_t1，在DDM_t2中按先行后列的空间顺序滑动子窗口ΔD_t2，按照步骤（a）中之公式分别计算子窗口ΔD_t1与移动子窗口ΔD_t2的互相关系数R，保存最大相关系数所在的子窗口ΔD_t2；通过两个子窗口ΔD_t1和ΔD_t2的中心坐标值可获取两者的距离dis₁₂和方向关系ang₁₂；子窗口ΔD_t1所在的沙波在t1至t2时间内的移动速率：v3＝dis₁₂÷t；

(d)在DDM_t1中按先行后列的空间顺序滑动子窗口ΔD_t1，然后按照步骤（c）可获取每个滑动子窗口ΔD_t1的最大移动速率v3(_i,j)和移动方向ang(_i,j)。

Images