发明名称 一种基于MBES的海底沙波地貌运动探测方法
摘要 本发明公开海洋测绘、海底地形地貌调查研究的方法,具体是指一种基于MBES的海底沙波地貌运动探测方法。本发明是采用分米级测量精度的多波束测深技术和亚米级定位精度的导航定位系统,通过基于MBES的海底沙波探测前准备、基于MBES的海底沙波地貌的探测、基于MBES实测数据的DDM构建、经二次探测与处理、以及基于DDM的沙波运动速率计算得到海底沙波地貌的运动速率和方向。该发明给出了海底沙波运动的精细探测方案和海底网格的快速构建方法,进而给出了海底沙波运动速率的准确计算方法,在海底高活动区进行海洋测绘、海洋调查与海底科学研究均具有重要的实际应用价值。
申请公布号 CN103389077B 申请公布日期 2014.05.07
申请号 CN201310317429.9 申请日期 2013.07.24
申请人 国家海洋局第二海洋研究所 发明人 吴自银;余威;李守军;尚继宏;赵荻能;周洁琼;金肖兵
分类号 G01C13/00(2006.01)I 主分类号 G01C13/00(2006.01)I
代理机构 浙江英普律师事务所 33238 代理人 陈小良
主权项 1.一种基于MBES的海底沙波地貌运动探测方法,其特征在于包括下列步骤:步骤1:基于MBES的海底沙波探测前准备(1)设备检测与标定:选用多波束系统的测深精度达到分米级的仪器,GPS定位精度达到亚米级;(2)潮位站布设:在海底沙波测量区布设2~4个临时验潮站,获取临时潮位资料Tide={tide<sub>i</sub>}<sub>i=2,4</sub>,以控制测区潮位;(3)声速剖面测量:开始勘测前,必须在测区附近进行声速剖面测量,至少获取测区一个声速剖面<img file="FDA0000467679030000011.GIF" wi="315" he="70" />n1为声速剖面层数;步骤2:基于MBES的海底沙波地貌的探测采取十字正交测线探测法或矩形区域探测法两种布线方式进行海底沙波的多波束探测,其中,十字正交测线探测法是指:在沙波区布设十字正交测线,其中一条平行沙波走向,一条垂直沙波走向;测线的宽度D应大于沙波的波长d,测线的长度L应大于沙波的横向延伸距离l;勘测时沿直线移动,多波束全开角勘测,全程GPS有差分信号;矩形区域探测法是:在沙波区布设矩形区域,矩形的长边应平行沙波走向,短边垂直沙波走向,矩形的长边A应大于沙波的横向延伸距离l,矩形的短边B应大于沙波波长d的两倍;勘测时沿直线移动,多波束全开角勘测,全程GPS有差分信号;形成首次探测原始数据集合Raw<sub>t1</sub>={raw1<sub>i</sub>},t1为海底地形测量时间;步骤3:基于MBES实测数据的DDM构建(1)数据处理:对获取的原始数据集合Raw<sub>t1</sub>={raw1<sub>i</sub>}进行潮位改正、吃水改正、声速改正和噪声点编辑处理后,形成处理后的离散水深数据集合Proc<sub>t1</sub>={(x<sub>i</sub>,y<sub>i</sub>,z<sub>i</sub>)}<sub>i=1,n</sub>;(2)数据分辨率评估:将离散水深数据集合正投影到二维平面上,按抽样法量算离散数据点的空间距离,要求95%的离散水深数据与其周围点的距离小于沙波运动估计速率v,否则降低直线移动速度,返回步骤2重新进行探测;(3)构建DDM:采取改进的距离反比加权(FIDW,Fast Inverse Distance Weighted)方法对离散水深数据集合Proc<sub>t1</sub>进行处理,形成DDM<sub>t1</sub>={dep<sub>(i,j)</sub>}<sub>i=1,n,j=1,m</sub>;FIDW方法计算公式为:<maths num="0001"><![CDATA[<math><mrow><msub><mi>dep</mi><mrow><mo>(</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>)</mo></mrow></msub><mo>=</mo><mo>[</mo><munderover><mi>&Sigma;</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mi>w</mi><mi>k</mi></msub><msub><mi>z</mi><mi>k</mi></msub><mo>]</mo><mo>/</mo><munderover><mi>&Sigma;</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mi>w</mi><mi>k</mi></msub><mo>;</mo></mrow></math>]]></maths>w<sub>k</sub>=1/d<sub>k</sub><sup>2</sup>;<maths num="0002"><![CDATA[<math><mrow><msub><mi>d</mi><mi>k</mi></msub><mo>=</mo><msqrt><msup><mrow><mo>(</mo><msub><mi>x</mi><mrow><mo>(</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>)</mo></mrow></msub><mo>-</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mo>(</mo><msub><mi>y</mi><mrow><mo>(</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>)</mo></mrow></msub><mo>-</mo><msub><mi>y</mi><mi>k</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mrow></math>]]></maths>式中,x<sub>k</sub>,y<sub>k</sub>,z<sub>k</sub>为离散水深点的横坐标、纵坐标和水深值,来自集合Proc<sub>t1</sub>;w<sub>k</sub>为离散水深数据点的权重值;dep<sub>(i,j)</sub>为网格值,x<sub>(i,j)</sub>和y<sub>(i,j)</sub>为网格横坐标和纵坐标值;设计网格点结构体DEP={dep<sub>(i,j)</sub>,A<sub>(i,j)</sub>,B<sub>(i,j)</sub>},<img file="FDA0000467679030000023.GIF" wi="641" he="135" />在读入离散数据的同时,根据每个离散水深数据的坐标(x<sub>k</sub>,y<sub>k</sub>),根据预先设定的作用距离dis,判断该点可作用的网格点;若d<sub>k</sub>≤dis时,该点参与网格化运算,否则直接跳过;当所有离散点遍历完成后,遍历一次网格即可获取每个网格点的水深差值dep<sub>(i,j)</sub>=A<sub>(i,j</sub>)/B<sub>(i,j)</sub>,当B<sub>(i,j)</sub>=0时该点为无效点;当参与网格dep<sub>(i,j)</sub>运算的点数少于n2时,该点也为无效点,n2为自然数,由系统缺省给定,也可由用户后期交互修改;步骤4:二次探测与处理(1)二次探测:按照步骤2中的测量方法,采取一定的时间间隔t进行二次重复测量,前后重复测量的航迹保持一致;测量时间间隔t要求为:t≥Δd÷v,Δd为Gps定位中误差,v为沙波运动估计速度;或者要求在一次风暴潮发生前后分别测量两次海底地形,用于研究风暴引起的海底沙波运动;重复探测形成原始数据集合Rwa<sub>t2</sub>={raw2<sub>i</sub>},t2为二次海底地形测量时间;(2)二次构建DDM:按照步骤3中(1)和(2)之方法,采取FIDW方法构建二次探测后的DDM<sub>t2</sub>;DDM<sub>t1</sub>和DDM<sub>t2</sub>为沙波两次测量形成的网格,两个网格的范围、行数、列数、行间距和列间距要完全一致;(3)构建差值ΔDDM:DDM<sub>t1</sub>和DDM<sub>t2</sub>按网格点逐点进行相减运算,生成新的差值网格ΔDDM=DDM<sub>t2</sub>-DDM<sub>t1</sub>;步骤5:基于DDM的沙波运动速率计算采用剖面对比法、差值法、或子窗相关法的方法来判断海底沙波是否运动,并计算海底沙波运动速率;其中,剖面对比法是指:垂直沙波走向,设计地形剖面线L<sub>0</sub>,分别在网格DDM<sub>t1</sub>、DDM<sub>t2</sub>和ΔDDM中进行相交运算,获取3条地形剖面线L<sub>t1</sub>、L<sub>t2</sub>和L<sub>Δt</sub>,在同一平面坐标系下按照相同参数用不同颜色绘制三条地形剖面线,对比沙波峰或谷的迁移距离Δl<sub>1</sub>,可以获取一处沙波的运动速率v1<sub>i</sub>=Δl<sub>1</sub>÷t;按照上述方法,垂直一条沙波的不同位置,可获取该沙波n处位置的运动速率,通过平均的方法可以获取该沙波的平均运动速率<img file="FDA0000467679030000031.GIF" wi="332" he="128" />差值法是指:基于期沙波差值网格ΔDDM,绘制海底三维地形图,如果沙波出现运动,海底地形表现为谷峰相间条纹,通过人机交互模式追踪条纹的脊线和谷底线,测量距离最近的脊线和谷底线空间距离Δl<sub>2</sub>,可以获取一处沙波的运动速率v2<sub>i</sub>=Δl<sub>2</sub>÷t;按照上述方法,可获取n条沙波的运动速率,通过平均的方法可以获取一片沙波的平均运动速率<img file="FDA0000467679030000032.GIF" wi="348" he="128" />子窗相关法是指:通过建立滑动子窗口的方式来分割两期海底网格的互相关来判断海底沙波的运动,在两期海底网格DDM<sub>t1</sub>和DDM<sub>t2</sub>中,建立相同大小的移动的矩形子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>,其行数和列数分别为m和n,通过子窗口的移动来计算两个时期沙波网格的相关系数,从而判断海底沙波的运动方向和速率;具体步骤如下:(a)两个子窗口的互相关计算公式为:<maths num="0003"><![CDATA[<math><mrow><mi>R</mi><mo>=</mo><mfrac><msub><mi>R</mi><mi>ab</mi></msub><mrow><msub><mi>R</mi><mi>a</mi></msub><msub><mi>R</mi><mi>b</mi></msub></mrow></mfrac></mrow></math>]]></maths>其中,R<sub>a</sub>、R<sub>b</sub>分别为子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>的方差,R<sub>ab</sub>为子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>的协方差,计算公式分别为:<maths num="0004"><![CDATA[<math><mrow><msub><mi>R</mi><mi>a</mi></msub><mo>=</mo><msqrt><mfrac><mn>1</mn><mrow><mi>m</mi><mo>-</mo><mn>1</mn></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><munderover><mi>&Sigma;</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>a</mi><mi>ij</mi></msub><mo>-</mo><mover><mi>a</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mrow></math>]]></maths><maths num="0005"><![CDATA[<math><mrow><msub><mi>R</mi><mi>b</mi></msub><mo>=</mo><msqrt><mfrac><mn>1</mn><mrow><mi>m</mi><mo>-</mo><mn>1</mn></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><munderover><mi>&Sigma;</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>b</mi><mi>ij</mi></msub><mo>-</mo><mover><mi>b</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mrow></math>]]></maths><maths num="0006"><![CDATA[<math><mrow><msub><mi>R</mi><mi>ab</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>m</mi><mo>-</mo><mn>1</mn></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><munderover><mi>&Sigma;</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mo>(</mo><msub><mi>a</mi><mi>ij</mi></msub><mo>-</mo><mover><mi>a</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow><mrow><mo>(</mo><msub><mi>b</mi><mi>ij</mi></msub><mo>-</mo><mover><mi>b</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow></mrow></math>]]></maths>a<sub>ij</sub>、b<sub>ij</sub>为两个子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>相应位置的水深值,<img file="FDA0000467679030000045.GIF" wi="132" he="80" />为子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>的平均水深值;(b)按照空间顺序在DDM<sub>t1</sub>和DDM<sub>t2</sub>中分别设置初始的子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>;(c)固定ΔD<sub>t1</sub>,在DDM<sub>t2</sub>中按先行后列的空间顺序滑动子窗口ΔD<sub>t2</sub>,按照步骤(a)中之公式分别计算子窗口ΔD<sub>t1</sub>与移动子窗口ΔD<sub>t2</sub>的互相关系数R,保存最大相关系数所在的子窗口ΔD<sub>t2</sub>;通过两个子窗口ΔD<sub>t1</sub>和ΔD<sub>t2</sub>的中心坐标值可获取两者的距离dis<sub>12</sub>和方向关系ang<sub>12</sub>;子窗口ΔD<sub>t1</sub>所在的沙波在t1至t2时间内的移动速率:v3=dis<sub>12</sub>÷t;(d)在DDM<sub>t1</sub>中按先行后列的空间顺序滑动子窗口ΔD<sub>t1</sub>,然后按照步骤(c)可获取每个滑动子窗口ΔD<sub>t1</sub>的最大移动速率v3(<sub>i,j</sub>)和移动方向ang(<sub>i,j</sub>)。
地址 310012 浙江省杭州市西湖区保俶北路36号