[发明专利]一种基于MBES的海底沙波地貌运动探测方法

权利要求书

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

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

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

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

（3）声速剖面测量：开始勘测前，必须在测区附近进行声速剖面测量,至少获取测区一个声速剖面Svp={v_i}_i=1,n1，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,FastInverseDistance Weighted)方法对离散水深数据集合Proc_t1进行处理，形成DDM_t1={dep_(i,j)}_i=1,n,j=1,m；

FIDW方法计算公式为：

dep(i,j)=[Σk=1nwkzk]/Σk=1nwk;]]>

w_k=1/d_k²

dk=(x(i,j)-xk)2+(y(i,j)-yk)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=RabRaRb]]>

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

Ra=1m-11n-1Σi=1mΣj-1n(aij-a‾)2]]>

Rb=1m-11n-1Σi=1mΣj=1n(bij-b‾)2]]>

Rαb=1m-11n-1Σi=1mΣj=1n(aij-a‾)(bij-b‾)]]>

a_ij、b_ij为两个子窗口△D_t1和△D_t2相应位置的水深值，a、b为子窗口△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)。