CN102999914B - 一种基于地形网格的大陆坡脚点自动识别方法 - Google Patents

Abstract

本发明涉及一种海底地形的识别方法，具体是指一种基于地形网格的大陆坡脚点自动识别方法。本发明是通过用直线或折线对地形网格模型进行分割生成二维地形剖面线，再对二维地形剖面首次求导生成坡度剖面线和二阶导数剖面线，再获取二阶导数剖面线极值点，再使用D-P算法获取二次简化后的D-P地形剖面，然后对D-P地形剖面进行二次求导，以及使用地形与坡度判断法识别并消除D-P剖面中凹包地形，最后在D-P地形、坡度和二阶导数剖面基础上，综合使用坡度、水深、二阶导数、凹凸性、连续性和分段性等判定方法形成自动识别方法。本发明的优点在我国沿海的测绘过程中，具有测绘数据准确、操作方便、以及受自然天气影响较上等。

Description

一种基于地形网格的大陆坡脚点自动识别方法

技术领域

本发明涉及一种海底地形的识别方法，具体是指一种基于海底地形网格生成地形剖面线，并能对海底地形剖面进行自动分析，从而自动识别出大陆坡脚点的技术方法，涉及到计算机图形学、海底科学和海洋划界领域。

背景技术

大陆坡脚点是200海里以外大陆架划界中最重要的界限点，将影响大陆坡脚点外推60海里线和1％沉积物厚度线及最终外部界限的精度。对二维地形剖面进行智能分析，从而自动识别大陆坡脚点，在200海里以外大陆架划界、二维地形剖面特征自动识别与二维地形剖面绘制方面均具有重要的应用价值。

1982年签署生效的《海洋法公约》第七十六条规定：“沿海国的大陆架包括其领海以外依其陆地领土的全部自然延伸，扩展到大陆边外缘的海底区域的海床和底土，如果从测算领海宽度的基线量起到大陆边的外缘的距离不到二百海里，则扩展到二百海里的距离”；沿海国如果主张从测算领海宽度的基线量起超过200海里，则应按照《公约》第七十六条和附件二第四条的相关要求划定200海里以外大陆架外部界限，并向大陆架界限委员会提交划界案。

200海里以外大陆架划界的关键证据是系列划界界限，包括：大陆坡脚点(简称FOS)、公式线(FOS+60M线和1％沉积物厚度线)、限制线(350M线及2500m+100M线)和外部界限等。其中，FOS是最重要的界限点，因为FOS是确定200海里以外大陆架系列界限的起点，将直接影响FOS+60M线和1％沉积物厚度线的准确性，并最终影响外部界限的坐标和所圈定的面积，如果FOS错误或者误差过大将对最终外部界限产生重大影响。

从检索的公开资料分析，目前国内无成熟的技术方法用于大陆坡脚点的自动识别，尽管有少数文献对200海里以外大陆架划界进行了介绍，但无翔实的大陆坡脚点自动生成方法之内容，现有技术方法不足以支持大陆坡脚点的识别和生成。

发明内容

本发明针对现有技术中的不足，提出一种利用海底地形网格自动生成地形剖面线，并通过对地形剖面进行二次简化和二次求导，进而对简化剖面进行综合分析，从而实现大陆坡脚点的自动分析与提取的技术方法。该发明可应用于网格模型自动生成任意二维剖面线、地形剖面特征自动识别与综合剖面绘制、海洋划界中关键界限点的生成与海底地貌类型的自动识别等技术环节。

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

一种基于地形网格的大陆坡脚点自动识别方法，其特征在于，包括下列步骤：

(1)网格切割：用直线f(x,y)对地形网格模型z_i,j=Grid(i,j)进行分割，生成二维地形剖面线；直线f(x,y)的起点和终点坐标分别为O_(x1,y1)和E_(x2,y2)。

x_i,j和y_i,j为网格模型第i行第j列网格点之横坐标和纵坐标值；z_i,j为网格模型第i行第j列网格点之水深值；

直线f(x,y)斜率k=(x1-x2)/(y1-y2),当y1=y2时，斜率表达为k=(y1-y2)/(x1-x2)；

(2)首次求导：对二维地形剖面首次求导，在海底地形网格模型Grid(i,j)基础上，通过直线f(x,y)与之相交运算，获得f(x,y)与网格模型Grid(i,j)交点坐标及水深值点集，生成坡度剖面线和二阶导数剖面线；

由步骤(1)和(2)生成初始数据点集G₀={g_i}，每个数据点包括距离、水深、坡度和二阶导数值；

f(x,y)与网格Grid(i,j)的交点之坐标和水深值dep(x,y)分下述情形计算：

(a)交点位于网格点Grid(i,j)：直接返回网格点坐标(x_i,j,y_i,j)与水深z_i,j，与起点O(x₀,y₀)的距离: $\mathrm{dis}(x,y)=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2};$

(b)交点位于行上：交点位于第i行、第j至j+1列，距离计算同步骤(a)；

交点纵坐标:y=y_i,j，当y1=y2时，y=y1；

交点横坐标:x=x1+(y-y1)×k,当y1=y2时，x=x_i,j；

交点水深值:dep(x,y)=z_i,j+(x-x_i,j)×(z_i,j+1-z_i,j)÷(x_i,j+1-x_i,j)；

(c)交点位于列上：交点位于第j列、第i至i+1行，距离计算同步骤(a)；

交点横坐标:x=x_i,j；

交点纵坐标:y=y1+(x-x1)×k，当y1=y2时,y=y1；

交点水深值:dep(x,y)=z_i,j+(y-y_i,j)×(z_i,j+1-z_i,j)÷(y_i,j+1-y_i,j)；

(d)交点位于网格中：交点位于第i行至i+1行、第j至j+1列之间，距离计算同步骤(a)；交点水深值dep(x,y)可按照距离反比加权方法求取；交点被4个相邻网格点所包围，则该点水深值为：

$\mathrm{dep}(x,y)=\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i{z}_i/\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i,$ 其中， $w_i=1/d_i^2$

z_i、w_i和d_i分别为4个网格点的水深值、计算权重值及与交点的距离值；

上述交点的横坐标和纵坐标值为O_(x1,y1)或E_(x2,y2)；

(3)首次简化：获取二阶导数剖面线极值点，以极值点为拐点生成新的首次简化地形剖面线；由步骤(3)生成首次简化后的数据点集G₁={g_i}；

(4)二次简化：使用D-P算法对极值剖面进行运算，保留符合D-P算法的数据点，从而获取二次简化后的D-P地形剖面；

(5)二次求导：使用步骤(2)之方法，对D-P地形剖面进行二次求导，获得新的基于D-P地形剖面的坡度剖面线和二阶导数剖面线；由步骤(4)和(5)生成二次简化并求导后的数据点集G₂={g_i}；

(6)消除凹包：使用地形与坡度判断法，识别并消除D-P剖面中凹包地形；由步骤(6)生成消除凹包地形后的数据点集G₃={g_i}；

(7)综合判断：在D-P地形、坡度和二阶导数剖面基础上，综合使用坡度、水深、二阶导数、凹凸性、连续性和分段性判定方法，自动识别地形剖面中的大陆坡脚点，并在二维剖面中自动标识；

分项判断步骤：

(a)坡度法：对点集G₃={g_i}中的平均坡度值ms_i进行分区统计，分别获取陆架与海盆平均坡度以及陆坡区平均坡度遍历点集G₃={g_i}，平均坡度为陆架或海盆，平均坡度为陆坡区；

(b)水深法：遍历点集G₃={g_i}，对平均坡度的数据点进行深度分类，并获得陆架平均水深值和海盆平均水深值再次遍历点集G₃={g_i}，水深值为陆架，水深值 ${\mathrm{dep}}_i\∈|\stackrel{\‾}{d_2}-\mathrm{\Δd},\stackrel{\‾}{d_2}+\mathrm{\Δd}|$ 为海盆；

(c)二阶导数法：大陆坡脚点FOS是陆坡至海盆海底坡度最大变化之点，也就是二阶导数极值点；

(d)凸包特性：大陆坡脚点FOS位于陆坡至海盆转折处，因此在地形上具有凸包特征，具体表现为二阶导数值为正值的数据点；

(e)分段性法：遍历点集G₃={g_i}，上坡度符合且下坡度符合的数据点初步判断为大陆坡脚点FOS；

(f)连续性法：遍历点集G₃={g_i}，记录每个点向剖面首点生长的距离bd_i及剖面尾点生长的距离pd_i，再次遍历点集G₃={g_i}，比较每个数据点的bd_i和pd_i值，生长距离最远者为大陆坡脚点FOS；

综合运用步骤(a)～(f)，同时满足步骤(c)～(f)条件的数据点为大陆坡脚点FOS。

根据联合国《海洋法公约》第七十六条及大陆架划界委员会技术准则之规定，大陆坡脚点位于大陆坡坡底坡度变化最大点之位置，也就是大陆坡底部二阶导数极值点位置。本发明在地形剖面中识别符合该要求的地形剖面数据点。本发明提出基于二阶导数极值点和D-P算法的二次拟合方法，以及原始剖面和D-P剖面的二次求导方法，并给出坡度、水深、二阶导数、凹凸性、连续性和分段性等判定方法等综合判定方法，最终实现了大陆坡脚点的快速自动识别。

有益效果：本发明在我国沿海的测绘过程中，具有测绘数据准确、操作方便、以及受自然天气影响较上等。

附图说明

图1本发明的算法总体流程示意图

图2本发明实施例1中典型地形网格模型及地形示意图，其中OE为实例地形剖面线f(x,y)所在位置；虚线网为网格模型Grid(i,j)，为M行、N列的海底地形网格

图3本发明实施例1中的原始地形剖面示意图

图4本发明实施例1中的首次简化后的极值点地形剖面示意图

图5本发明实施例1中二次简化后的D-P地形剖面示意图

图6本发明实施例1中消除凹包后的地形剖面及识别的大陆坡脚点（TOS:大陆坡折点；FOS:大陆坡脚点）的示意图

图7本发明实施例1中的标准地形剖面及识别的大陆坡脚点示意图，其中TOS:大陆坡折点；FOS:大陆坡脚点；图中波浪线为二阶导数曲线；

图8本发明实施例1中的复杂陆坡的地形剖面及识别的大陆坡脚点示意图，其中TOS:大陆坡折点；FOS:大陆坡脚点；图中波浪线为二阶导数曲线；

图9本发明实施例1中的存在海山的地形剖面及识别的大陆坡脚点示意图，其中TOS:大陆坡折点；FOS:大陆坡脚点；图中波浪线为二阶导数曲线；

图10本发明实施例1中的陆坡有隆起的地形剖面及识别的大陆坡脚点示意图，其中TOS:大陆坡折点；FOS:大陆坡脚点；图中波浪线为二阶导数曲线；

具体实施方式

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

实施例1

如图1所示的流程：

(1)切割地形网格生成地形剖面线

首先获取一个大陆边缘海底地形网格模型z_i,j=Grid(i,j)，该网格共M行、N列，由M×N个规则排列的水深点组成，如：z_i,j为网格模型第i行第j列之水深值。该网格模型需包括陆架、陆坡和海盆等地形单元。陆架、陆坡和海盆为地貌学术语。陆架，一般指靠近大陆且水深小于200m的平坦海底地形区域。陆坡指紧邻陆架，且海底水深急剧加深的海底地形区，水深约百余米至数千米。海盆指紧邻陆坡且海底较为平坦的海底地形区，水深可达千余米或数千米。陆坡区及海盆区水深与大陆边缘性质相关，在世界各大陆边缘中，陆坡区及海盆区的具体水深是变化的。可使用多种现有工具或者方法生成海底地形网格模型。

本实施例中，首先网格切割：用直线f(x,y)对地形网格模型z_i,j=Grid(i,j)进行分割，生成二维地形剖面线；直线f(x,y)的起点和终点坐标分别为O_(x1,y1)和E_(x2,y2)。

x_i,j和y_i,j为网格模型第i行第j列网格点之横坐标和纵坐标值；z_i,j为网格模型第i行第j列网格点之水深值；

直线f(x,y)斜率k=(x1-x2)/(y1-y2),当y1=y2时，斜率表达为k=(y1-y2)/(x1-x2)；

然后进行首次求导：对二维地形剖面首次求导，在海底地形网格模型Grid(i,j)基础上，通过直线f(x,y)与之相交运算，获得f(x,y)与网格模型Grid(i,j)交点坐标及水深值点集，生成坡度剖面线和二阶导数剖面线；

由步骤(1)和(2)生成初始数据点集G₀={g_i}，每个数据点包括距离、水深、坡度和二阶导数值；

f(x,y)与网格Grid(i,j)的交点之坐标和水深值dep(x,y)分下述情形计算：

(d)交点位于网格点Grid(i,j)：直接返回网格点坐标(x_i,j,y_i,j)与水深z_i,j，与起点O(x₀,y₀)的距离: $\mathrm{dis}(x,y)=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2};$

(e)交点位于行上：交点位于第i行、第j至j+1列，距离计算同步骤(a)；

交点纵坐标:y=y_i,j，当y1=y2时，y=y1；

交点横坐标:x=x1+(y-y1)×k,当y1=y2时，x=x_i,j；

交点水深值:dep(x,y)=z_i,j+(x-x_i,j)×(z_i,j+1-z_i,j)÷(x_i,j+1-x_i,j)；

(f)交点位于列上：交点位于第j列、第i至i+1行，距离计算同步骤(a)；

交点横坐标:x=x_i,j；

交点纵坐标:y=y1+(x-x1)×k，当y1=y2时,y=y1；

交点水深值:dep(x,y)=z_i,j+(y-y_i,j)×(z_i,j+1-z_i,j)÷(y_i,j+1-y_i,j)；

在海底地形网格模型Grid(i,j)基础上，通过直线或折线f(x,y)与之相交运算，获得f(x,y)与网格模型Grid(i,j)交点坐标及水深值点集G₀={g_i}，按点集G₀与f(x,y)首点O_(x1,y1)距离的顺序进行保存，并形成横坐标为距离、纵坐标为水深值的初始地形剖面线g₀(x,y)。为方便回溯地形剖面点的平面坐标值，需设计合理结构同时保存点集G₀的平面坐标、距离和水深值。

f(x,y)与网格Grid(i,j)的交点之水深值dep(x,y)分下述情形计算：

(a)交点位于网格点。直接返回网格点坐标与水深。与起点O(x₀,y₀)的距离: $\mathrm{dis}(x,y)=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}.$

(b)交点位于行上。交点位于第i行、第j至j+1列，距离计算同步骤(a)。交点水深值:dep(x,y)=z_i,j+(x-x_i,j)×(z_i,j+1-z_i,j)÷(x_i,j+1-x_i,j)。x(i,j)和z(i,j)为网格第i行第j列的横坐标及水深值。

(c)交点位于列上。交点位于第j列、第i至i+1行，距离计算同步骤(a)。交点水深值:dep(x,y)=z_i,j+(y-y_i,j)×(z_i,j+1-z_i,j)÷(y_i,j+1-y_i,j)。y(i,j)为网格第i行第j列的纵坐标。

(d)交点位于网格中。交点位于第i行至i+1行、第j至j+1列之间，距离计算同步骤(a)。交点水深值dep(x,y)可按照距离反比加权方法求取。交点被4个相邻网格点所包围，则该点水深值为：

$\mathrm{dep}(x,y)=\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i{z}_i/\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i,$ 其中， $w_i=1/d_i^2$

z_i、w_i和d_i分别为4个网格点的水深值、计算权重值及与交点的距离值；上述交点的横坐标和纵坐标值为O_(x1,y1)或E_(x2,y2)。

形成的原始地形剖面要符合“陆架－陆坡－海盆”特征，也就是说地形剖面线必须穿过大陆边缘的陆架、陆坡和海盆区，才会存在大陆坡脚点，该条件也是后续工作的基础。

(2)首次求导生成坡度和二阶导数剖面

对地形剖面线g₀(x,y)进行首次求导，获得坡度剖面g₀’(x,y)和二阶导数剖面g₀”(x,y)。

(a)形成坡度剖面。对地形剖面中第i个数据点g_i而言，前后点分别为g_i-1和g_i+1。已知该点水深值(z_i)和距离值(d_i)，对g_i点进行求导，可获得该点的上坡度(us_i)、下坡度(ds_i)和平均坡度(ms_i)。

us_i=(z_i-1-z_i)÷|d_i-1-d_i|

ds_i=(z_i+1-z_i)÷|d_i+1-d_i|

ms_i=(|us_i|+|ds_i|)÷2

地形剖面的首点上坡度无法直接计算，假定与剖面第2点相同。尾点下坡度也无法计算，假定与倒数第2点下坡度值相同。

坡度剖面g₀’(x,y)中每个点g’_i包括坐标、距离、上坡度、下坡度、水深和平均坡度等数值，其坐标、距离、水深与原地形剖面相同。

(b)形成二阶导数剖面。对坡度剖面g’(x,y)进行求导，可以获取该剖面的二阶导数剖面g₀”(x,y)，第i点二阶导数sec_i为：

sec_i=(z_i+1-z_i)÷|d_i+1-d_i|²

二阶导数剖面g₀”(x,y)中每个点g”_i包括坐标、距离、水深和二阶导数数值，其坐标、距离、水深与原地形剖面相同。

(3)首次简化形成极值地形剖面

根据上述步骤(2)中(a)获取了二阶导数剖面g”(x,y)，该剖面中存在很多极值点，由极值点水深组成新的简化地形剖面g₁(x,y)。

二阶导数g”_i前后点分别为g”_i-1和g”_i+1，如果g”_i-1和g”_i+1数值符号相同，则该点为极值点，并形成新的点集G₁={g_i}，每个点包括距离、水深、坡度和二阶导数值，点集G₁={g_i}形成新的极值点地形剖面g₁(x,y)、坡度剖面g₁’(x,y)和二阶导数剖面g₁”(x,y)。

极值点地形剖面g₁(x,y)与原地形剖面相比，进行了首次简化，仅保留了符合二阶导数极值点特征的数据点。

(4)二次简化形成D-P地形剖面

Douglas-Peucker算法(简称D-P算法)是由DavidDouglas和ThomasPeucker于1973年提出的一种曲线抽稀算法，可对曲线大量冗余点进行精简并保留曲线最基本特征。简单的说，根据曲线离散度给定一个初始偏差值，将曲线首尾点相连形成直线，并查询曲线所有拐点中距离该直线最远的点。如果该点与直线的距离小于初始偏差值，则删除曲线所有拐点并返回。如果该点与直线的距离大于初始偏差值，则保留该点作为特征点，并以曲线首点和该点形成新的直线并重新距离最远的点。依此类推，直至搜索完曲线所有点并保留下特征点。整体算法是D-P算法的突出优点，可保留曲线最大弯曲形态的点。D-P算法已是计算机图形学曲线简化常规方法。

由步骤(3)获取的极值点地形剖面g₁(x,y)，采用D-P算法，在地形剖面g₁(x,y)中计算，产生符合D-P算法和偏差值的数据点集G₂={g_i}，并形成新的地形剖面g₂(x,y)，该剖面是二次简化后的地形剖面，仅保留了极少数数据点，每个点的水深值和距离值与原始剖面g₀(x,y)中同位置点相同。

(5)二次求导形成新的坡度和二阶导数剖面

由步骤(4)形成二次简化后的D-P地形剖面g₂(x,y)，采用步骤(2)之方法，对g₂(x,y)所表达的地形进行再次求导，形成新的坡度剖面g₂’(x,y)和二阶导数剖面g₂”(x,y)。

(6)识别并消除凹包地形

凹包地形影响大陆坡脚点的识别和判断，要在最终确定前予以消除。所谓凹包地形是指在地形剖面中隆起地形，也就是水深值比相邻点水深浅的地形。按照笛卡尔坐标系表达，横轴为距离、纵轴为水深的地形剖面，隆起地形在笛卡尔坐标系中表达为下凹形态，因此，将之定义为凹包。在地形剖面中下凹的地形,在笛卡尔坐标系中实际为凸包。在该步骤要消除二次简化后的D-P地形剖面中的凹包地形。

由步骤(4)形成二次简化后的D-P地形剖面g₂(x,y)，由步骤(5)形成二次求导后的坡度剖面g₂’(x,y)，如果g₂’(x,y)中第i点上坡度(us_i)和下坡度(ds_i)符号相同，且都为正值，则该点为凹包点，应该消除。采用二次循环，遍历地形剖面g₂(x,y)中所有点，消除符合凹包特征的点后形成新的点集G₃={g_i}，每个数据点包括地形、坡度和二阶导数值，从而形成新的地形剖面g₃(x,y)、坡度剖面g₃’(x,y)和二阶导数剖面g₃”(x,y)。

(7)综合识别并判断大陆坡脚点

通过步骤(1)至步骤(6)，原始地形剖面已大幅简化但保留最基本特征，消除了局部小地形和凹包地形对大陆坡脚点识别的影响，海底地形剖面已经简化为易于识别的“陆架-陆坡-海盆”三段式地形剖面，陆架和海盆地形平坦，海底坡度较小，陆坡地形呈现斜坡形状，坡度稳定，在陆坡与海盆交接处呈现凸包特征。

基于上述简化后的地形特征，可采用坡度、水深、二阶导数、凸包特征、连续性和分段性对步骤(6)形成的地形剖面g₃(x,y)进行查询，遍历点集G₃={g_i}识别并判断出大陆坡脚点，分项判断依据如下：

(a)坡度法。经2次简化后的地形剖面，其平均坡度值ms_i分布两个坡度区间，分别对应陆架与海盆区，以及陆坡区。对点集G₃={g_i}中的平均坡度值ms_i进行分区统计，分别获取陆架与海盆平均坡度以及陆坡区平均坡度

遍历点集G₃={g_i}，平均坡度为陆架或海盆，平均坡度为陆坡区，其中Δg为坡度限差值。通过该步骤可识别出属于陆坡地形区的数据点。

(b)水深法。根据步骤(a)，获取了平坦区的平均坡度但仅根据坡度无法区分陆架和海盆区。陆架水深和海盆水深差别很大，陆架水深一般小于250m，海盆水深一般大于1000m，甚至达到4000m以深。遍历点集G₃={g_i}，对平均坡度的数据点进行深度分类，并获得陆架平均水深值和海盆平均水深值

遍历点集G₃={g_i}，水深值定义为陆架，水深值定义为海盆。通过该步骤，识别出属于陆架和海盆地形区的数据点。

(c)二阶导数。大陆坡脚点FOS是陆坡至海盆中海底坡度最大变化之点，也就是二阶导数极值点。

(d)凸包法。大陆坡脚点FOS位于陆坡至海盆转折处，因此在地形上具有凸包特征。表现为二阶导数值为正值的数据点。

(e)分段性法。大陆坡脚点FOS位于陆坡至海盆转折处，在地形上具有“陆坡-海盆”的分段性特征。大陆坡脚点FOS点向陆方向(向剖面首点)，为陆坡，地形陡峭，坡度符合特征。大陆坡脚点FOS点向海方向(向剖面尾点)，坡度符合 ${\mathrm{ms}}_i\∈|\stackrel{\‾}{g_1}-\mathrm{\Δg},\stackrel{\‾}{g_1}+\mathrm{\Δg}|$ 特征。

遍历点集G₃={g_i}，符合上坡度且下坡度的数据点初步判断为大陆坡脚点FOS。

(f)连续性法。经步骤(a)至(e)，可能有多个点符合条件。可根据曲线连续性进一步判断并排除异常点。经简化后的地形剖面单一，因此根据平均坡度值ms_i，从大陆坡脚点FOS向曲线首点追溯，并按照步骤(e)所述分段性特征，可以追溯出完整的陆坡。向曲线尾点追溯，按照步骤(e)所述分段性特征，可以追溯出完整的海盆。

遍历点集G₃={g_i}，记录每个点向首点追溯的生长距离bd_i和向尾点追溯的生长距离pd_i。再次遍历点集G₃={g_i}，比较每个数据点的bd_i和pd_i值，距离最远者为大陆坡脚点FOS。经该步骤确定的FOS点具有唯一性。

综合运用步骤(a)至(f)，同时满足(c)～(f)条件的数据点为大陆坡脚点FOS。

针对一个具体的地形剖面，如果无法正确识别出符合条件的大陆坡脚点FOS，可能因两种原因造成。情形之一，可能所选的地形剖面不符合“陆架-陆坡-海盆”假设，不存在“陆坡-海盆”转折之地形。在该种情形下，可返回步骤(1)重新生成符合条件的地形剖面，然后重复(2)至(7)之步骤。情形之二，是采用D-P算法筛选了过多的数据点，导致剩余的数据点不足以描述“陆架-陆坡-海盆”的基本地形特征。在该种情形下，返回步骤(4)，重新设置较小的距离偏差值，重复(4)至(7)之步骤，直至识别出符合步骤(7)要求的大陆坡脚点FOS。

通过步骤(1)可以基于海底地形网格生成任意的海底二维地形剖面，因此，本发明可适用于基于网格的二维地形剖面的自动生成。基于相同方法，也可适用其他任意类型网格的任意二维剖面线生成。

通过步骤(2)可以计算地形剖面的坡度和二阶导数，因此，本发明可适用基于地形剖面的坡度和二阶导数剖面的自动生成，可应用于包括地形、坡度和二阶导数的综合剖面绘制。

通过步骤(3)和(4)可以自动生成保留原剖面骨架特征的简化地形剖面，因此，本发明适用复杂地形剖面的简化，可应用于二维剖面的骨架点生成。

通过地形、坡度和二阶导数可以判别曲线的凹凸特性，通过步骤(6)可以识别凹包地形，通过步骤(a)和(b)可以识别出陆架、陆坡和海盆地形。因此，本发明适用于二维地形剖面特征的自动分析和识别。

通过步骤(1)～(7)，可自动识别大陆坡脚点，因此，本发明适用于基于地形网格的大陆坡脚点自动判断和识别，可应用于海洋划界中FOS+60M线与1％沉积物厚度线的生成。考虑到陆架至陆坡转折处的大陆坡折点与大陆坡脚点特征的类似性，因此，本发明也适用于大陆坡折点的自动识别与生成。

实施例2

对于给定的海底地形网格模型，我们通过切割网格、首次求导、首次简化地形、二次简化地形、二次求导、消除凹包地形、综合判断等7个步骤实现大陆坡脚点的自动识别(总体技术流程见图1)，具体实施方式如下：

(1)切割网格

使用直线f(x,y)切割海底地形网格模型Grid(i,j)，进行相交运算，获得原始地形数据点集G₀={g_i}，每个数据点包括距离、水深和平面坐标值。点集G₀={g_i}按与f(x,y)首点O_(x1,y1)距离的顺序进行保存，并形成横坐标为距离、纵坐标为水深值的原始地形剖面线g₀(x,y)。

地形剖面g₀(x,y)要符合“陆架－陆坡－海盆”特征，也就是说地形剖面线必须穿过大陆边缘的陆架、陆坡和海盆区，否则重新采用人机交互方式生成新的地形剖面。

典型海底地形网格模型见图2，剖面线位置见图2中直线生成的地形剖面g₀(x,y)见图3。

(2)首次求导

对地形剖面线g₀(x,y)进行首次求导，获得坡度剖面g₀’(x,y)和二阶导数剖面g₀”(x,y)。

坡度剖面g₀’(x,y)中每个数据点g’_i包括坐标、距离、上坡度、下坡度、水深和平均坡度等数值，其坐标、距离、水深与原地形剖面相同。

二阶导数剖面g₀”(x,y)中每个点g”_i包括坐标、距离、水深和二阶导数数值，其坐标、距离、水深与原地形剖面相同。

由g₀(x,y)、g₀’(x,y)和g₀”(x,y)中距离、地形、坡度和二阶导数值共同组成点集G₀={g_i}。

(3)首次简化

根据上述步骤(a)获取了二阶导数剖面g₀”(x,y)，该剖面存在很多极值点。仅保留原剖面极值点，并由极值点组成新的简化地形剖面g₁(x,y)。极值点判断依据为：二阶导数g”_i前后点分别为g”_i-1和g”_i+1，如果g”_i-1和g”_i+1数值符号相同，则该点为极值点。形成新的点集G₁={g_i}，每个点包括距离和水深值，点集G₁={g_i}形成新的极值点地形剖面g₁(x,y)。极值点地形剖面g₁(x,y)与原地形剖面相比，进行了首次简化，仅保留了部分符合二阶导数极值点特性的水深点。

首次简化后的地形剖面见图4。

(4)二次简化

由步骤(3)获取的极值点地形剖面g₁(x,y)，采用D-P算法，在极值点地形剖面g₁(x,y)中计算符合要求的数据点集G₂={g_i}，并形成新的地形剖面g₂(x,y)，该剖面是二次简化后的地形剖面，仅保留了极少数数据点，每个点的水深值和距离值与原始剖面g₀(x,y)相同。

二次简化后的D-P剖面见图5。

(5)二次求导

由步骤(4)形成二次简化后的D-P地形剖面g₂(x,y)，采用步骤(2)之方法，对g₂(x,y)所表达的地形进行再次求导，形成新的坡度剖面g₂’(x,y)和二阶导数剖面g₂”(x,y)。

由g₂(x,y)、g₂’(x,y)和g₂”(x,y)中距离、地形、坡度和二阶导数值共同组成数据点集G₂={g_i}。

(6)消除凹包

由步骤(4)形成二次简化后的D-P地形剖面g₂(x,y)，由步骤(5)形成二次求导后的坡度剖面g₂’(x,y)，如果g₂’(x,y)中第i点上坡度(us_i)和下坡度(ds_i)符号相同，且都为正值，则该点为凹包点，应该消除。

采用二次循环，遍历地形剖面g₂(x,y)中所有点，消除符合凹包特征的点后形成新的点集G₃={g_i}，每个数据点包括地形、坡度和二阶导数值，从而形成新的地形剖面g₃(x,y)、坡度剖面g₃’(x,y)和二阶导数剖面g₃”(x,y)。

消除凹包后的剖面见图6。

(7)综合判断

通过步骤(1)至步骤(6)，获取了简化后的综合剖面，我们采用水深、坡度、二阶导数、凸包特征、连续性和分段性对步骤(6)形成的地形剖面g₃(x,y)进行查询，遍历点集G₃={g_i}识别并判断出大陆坡脚点。分项判断依据如下：

(a)坡度法。对点集G₃={g_i}中的平均坡度值ms_i进行分区统计，分别获取陆架与海盆平均坡度以及陆坡区平均坡度遍历点集G₃={g_i}，平均坡度为陆架或海盆，平均坡度为陆坡区。

(b)水深法。遍历点集G₃={g_i}，对平均坡度的数据点进行深度分类，并获得陆架平均水深值和海盆平均水深值再次遍历点集G₃={g_i}，水深值为陆架，水深值 ${\mathrm{dep}}_i\∈|\stackrel{\‾}{d_2}-\mathrm{\Δd},\stackrel{\‾}{d_2}+\mathrm{\Δd}|$ 为海盆。

(c)二阶导数。大陆坡脚点FOS是陆坡至海盆海底坡度最大变化之点，也就是二阶导数极值点。

(d)凸包特性。大陆坡脚点FOS位于陆坡至海盆转折处，因此在地形上具有凸包特征。具体表现为二阶导数值为正值的数据点。

(e)分段性法。大陆坡脚点相邻的前后点分别为陆坡和海盆，陆坡和海盆通过步骤(a)和(b)已经识别。遍历点集G₃={g_i}，上坡度符合且下坡度符合的数据点初步判断为大陆坡脚点FOS。

(f)连续性法。遍历点集G₃={g_i}，记录每个点向剖面首点生长的距离bd_i及剖面尾点生长的距离pd_i。再次遍历点集G₃={g_i}，比较每个数据点的bd_i和pd_i值，生长距离最远者为大陆坡脚点FOS。

综合运用步骤(a)至(f)，同时满足(c)～(f)条件的数据点为大陆坡脚点FOS。最终识别的大陆坡脚点见图6。

我们使用C++编程语言在2.40GHz的Intel(R)Core2,2GB内存的PC机上实现了本发明之算法，并使用典型大陆边缘之海底地形网格模型(见图2)，自陆架至海盆切割典型海底地形剖面，采用本发明之技术方法准确识别出大陆坡脚点FOS。本发明总体技术流程见图1，典型剖面识别出大陆坡脚点的具体步骤见图2～图6。还对不同类型地形剖面进行了自动识别大陆坡脚点实验，结果如图7所示。采用本发明相同之方法，也实现了大陆坡折点的自动识别，见图6和图7中大陆坡折点TOS。

实施例3

与实施例1相同方法，对复杂陆坡的地形剖面进行自动识别，可得图8所示结果。

实施例4

与实施例1相同方法，对存在海山的地形剖面进行自动识别，可得图9所示结果。

实施例5

与实施例1相同方法，对陆坡有隆起的地形剖面进行自动识别，可得图10所示结果。

Claims (1)

1.一种基于地形网格的大陆坡脚点自动识别方法，其特征在于，包括下列步骤：

(1)网格切割：用直线f(x,y)对地形网格模型z_i,j＝Grid(i,j)进行分割，生成二维地形剖面线；直线f(x,y)的起点和终点坐标分别为O_(x1,y1)和E_(x2,y2)；

x_i,j和y_i,j为网格模型第i行第j列网格点之横坐标和纵坐标值；z_i,j为网格模型第i行第j列网格点之水深值；

直线f(x,y)斜率k＝(x1-x2)/(y1-y2),当y1=y2时，斜率表达为k＝(y1-y2)/(x1-x2)；

(2)首次求导：对二维地形剖面线首次求导，在地形网格模型Grid(i,j)基础上，通过直线f(x,y)与之相交运算，获得f(x,y)与网格模型Grid(i,j)交点坐标及水深值点集，生成坡度剖面线和二阶导数剖面线；其中，对地形剖面线g₀(x,y)进行首次求导，获得坡度剖面g₀’(x,y)和二阶导数剖面g₀”(x,y)；

由步骤(1)和(2)生成初始数据点集G₀＝{g_i}，每个数据点包括距离、水深、坡度和二阶导数值；

f(x,y)与网格Grid(i,j)的交点之坐标和水深值dep(x,y)分下述情形计算：

(a)交点位于网格点Grid(i,j)：直接返回网格点坐标(x_i,j,y_i,j)与水深z_i,j，与起点O(x₁,y₁)的距离: $\mathrm{dis}(x,y)=\sqrt{{(x-x_1)}^2+{(y-y_1)}^2};$

(b)交点位于行上：交点位于第i行之上、第j至j+1列之间，距离计算同步骤(a)；

交点纵坐标:y＝y_i,j，当y1=y2时，y=y1；

交点横坐标:x＝x1+(y-y1)×k,当y1=y2时，x＝x_i,j；

交点水深值:dep(x,y)＝z_i,j+(x-x_i,j)×(z_i,j+1-z_i,j)÷(x_i,j+1-x_i,j)；

(c)交点位于列上：交点位于第j列之上、第i至i+1行之间，距离计算同步骤(a)；

交点横坐标:x＝x_i,j；

交点纵坐标:y＝y1+(x-x1)×k，当y1=y2时,y=y1；

交点水深值:dep(x,y)＝z_i,j+(y-y_i,j)×(z_i,j+1-z_i,j)÷(y_i,j+1-y_i,j)；

(d)交点位于网格中：交点位于第i行至i+1行、第j至j+1列之间，距离计算同步骤(a)；交点水深值dep(x,y)可按照距离反比加权方法求取；交点被4个相邻网格点所包围，则该点水深值为：

$\mathrm{dep}(x,y)=\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i{z}_i/\underset{i=1}{\overset{i=4}{\mathrm{\Σ}}}{w}_i,$ 其中， $w_i=1/d_i^2$

z_i、w_i和d_i分别为4个网格点的水深值、计算权重值及与交点的距离值；

上述步骤（d)中所述交点的横坐标和纵坐标值为O(_x1,y1)或E(_x2,y2)；

(3)首次简化：获取二阶导数剖面线极值点，以极值点为拐点生成新的首次简化地形剖面线；由步骤(3)生成首次简化后的数据点集G₁＝{g_i}；

(4)二次简化：使用D-P算法对极值剖面进行运算，保留符合D-P算法的数据点，从而获取二次简化后的D-P地形剖面；

(5)二次求导：使用步骤(2)之方法，对D-P地形剖面进行二次求导，获得新的基于D-P地形剖面的坡度剖面线和二阶导数剖面线；由步骤(4)和(5)生成二次简化并求导后的数据点集G₂＝{g_i}；

(6)消除凹包：使用地形与坡度判断法，识别并消除D-P剖面中凹包地形；由步骤(6)生成消除凹包地形后的数据点集G₃＝{g_i}；

(7)综合判断：在D-P地形、坡度和二阶导数剖面基础上，综合使用坡度、水深、二阶导数、凹凸性、连续性和分段性判定方法，自动识别地形剖面中的大陆坡脚点，并在二维剖面中自动标识；

分项判断步骤：

(a)坡度法：对点集G₃＝{g_i}中的平均坡度值ms_i进行分区统计，分别获取陆架与海盆平均坡度以及陆坡区平均坡度遍历点集G₃＝{g_i}，平均坡度 ${\mathrm{ms}}_i\∈|\stackrel{\‾}{g_1}-\mathrm{\Δg},\stackrel{\‾}{g_1}+\mathrm{\Δg}|$ 为陆架或海盆，平均坡度 ${\mathrm{ms}}_i\∈|\stackrel{\‾}{g_2}-\mathrm{\Δg},\stackrel{\‾}{g_2}+\mathrm{\Δg}|$ 为陆坡区；其中Δg为坡度限差值；

(b)水深法：遍历点集G₃＝{g_i}，对平均坡度的数据点进行深度分类，并获得陆架平均水深值和海盆平均水深值再次遍历点集G₃＝{g_i}，水深值为陆架，水深值为海盆；其中Δd为水深限差值；

(c)二阶导数法：大陆坡脚点FOS是陆坡至海盆海底坡度最大变化之点，也就是二阶导数极值点；

(d)凸包特性：大陆坡脚点FOS位于陆坡至海盆转折处，因此在地形上具有凸包特征，具体表现为二阶导数值为正值的数据点；

(e)分段性法：遍历点集G₃＝{g_i}，上坡度符合且下坡度符合的数据点初步判断为大陆坡脚点FOS；

(f)连续性法：遍历点集G₃＝{g_i}，记录每个点向剖面首点生长的距离bd_i及剖面尾点生长的距离pd_i，再次遍历点集G₃＝{g_i}，比较每个数据点的bd_i和pd_i值，生长距离最远者为大陆坡脚点FOS；

综合运用步骤(a)～(f)，同时满足步骤(c)～(f)条件的数据点为大陆坡脚点FOS。