CN102831645A - 一种应用于海底地形的数字高程模型的建立方法 - Google Patents

Abstract

本发明公开了一种海底数字高程模型的建立方法，采用内插方法，该方法把协调Delaunay三角化的思想融入自然邻点插值，提供一种新的含非凸本质边界条件的插值思想。本方法主要包括五部分内容：构建嵌入特征约束后符合Delaunay三角化特性的三角网，寻找插值点的准自然邻点，构建二阶常规Voronoi单元，由二阶常规Voronoi单元构建二阶约束Voronoi单元和计算插值点的属性值。本发明提出的方法以光滑的插值函数对地形函数进行逼近，具有计算区域固定、光滑度高、满足非凸边界条件、易于实施等优点。

Description

一种应用于海底地形的数字高程模型的建立方法

技术领域

本发明属于地理信息系统空间技术领域，涉及海底地形的拟合，特别涉及一种应用于海底地形的数字高程模型的建立方法。

背景技术

数字高程模型是描述地表起伏形态特征的空间数据模型，是由地面规则格网点的高程值构成的矩阵，形成栅格结构数据集。数字高程模型在建立过程的关键环节是格网点上高程的计算，在数学上属于数值分析中的插值问题。任意一种插值方法都是基于原始地形起伏变化的连续性即邻近的数据点间的相关性来求解待定点的高程。

基于地形的插值方法从数据分布规律来讲，分为基于规则分布的插值方法和基于不规则分布的插值方法；按插值点分布范围，分为整体内插、局部内插和逐点内插；按插值函数与采样点的关系，可分为曲面通过所有采样点的二维插值方法和曲面不通过所有采样点的拟合方法；从内插曲面的数学性质来讲，有多项式、样条、最小二乘配置等插值函数；从对地形函数理解角度，有克里金法、多层曲面叠加法、加权平均法等。

整体内插是在整个区域用一个数学函数来表达地形曲面，要求采样点的个数大于或等于多项式的待定系数，缺点是不易得到稳定的数值解，且求解速度较慢。

局部内插是将区域按一定方法进行分块，对每块区域根据地形曲面特征单独进行曲面拟合，分块插值简化了地形曲面的形态，使得每块都可用不同曲面进行表达，缺点是如何进行分块和不能保证各块曲面之间的连续性。局部内插方法有最小二乘配置法、克里金法、样条函数法和多层曲面叠加法。最小二乘配置法认为一个采样数据由趋势、信号和误差三部分构成，分块地形曲面通过多项式来确定变化趋势，信号反映局部数据点之间的相关性，由数据点之间的协方差函数表达，最后通过误差平方和最小原则求解各个参数，最小二乘配置法虽然数学理论基础严密，但实验结果表明它不一定能够得到好的拟合效果，原因在于实施最小二乘配置法的前提是处理对象必须属于遍历性随机平稳随机过程，实际地形曲面不一定满足这一条件，而且地形之间的自相关性不仅与距离有关，也与方向有关。克里金法又称空间自协方差最佳内插法，采用半方差来构造推估用的方差协方差阵。样条函数法是采用三次多项式对采样曲线进行分段修匀，每次分段拟合利用少数采样点的观测值，为保证各个分块之间的平滑过度，要按照弹性力学条件设立分块之间的连续性条件，然而地面并不是狭义的刚体，也不具备弹性力学光滑性的条件，因此虽然样条函数具有严密的理论基础，但不是数字高程模型内插的理想数学模型。多层曲面叠加法的理论依据是“任何一个光滑的数学表面总可以用一系列简单的数学表面来叠加逼近”，简单数学表面也称为核函数，多层曲面叠加法的优点是核函数设计的灵活性和可控性，虽然核函数选择比较灵活，但难以通过一个确定的函数表示各种地形变化，同时多层曲面叠加函数的处理过程繁琐，所以在数字高程模型建立中并不常用。

逐点插值是以插值点为中心，确定一个邻域范围，用落在邻域范围内的采样点计算内插点的高程。逐点内插与局部内插的区别是，局部内插的分块范围一经确定，凡是落在该范围内的点都要进行曲面拟合，而逐点内插法邻域的范围大小、形状、采样点的个数都随内插点变动，一套数据只能用来进行一个内插点的计算。逐点内插法由于内插效率较高而成为目前数字高程模型生产中常采用的方法。逐点内插法有加权平均法和Voronoi图法。加权平均法通常采用与距离有关的权函数，而事实上距离难以很好地描述空间相邻性。对于离散数据点之间的空间相邻性描述，Voronoi图是一种很好的数学工具。Voronoi图又称Dirichlet镶嵌（Dirichlet tessellation），在二维空间中也称为Thiessen多边形。其概念由Dirichlet于1850年提出（见参考文献[1]：Gustav Lejeune Dirichlet.die Reduktion derpositiven quadratischen Formen mit drei unbestimmten ganzen Zahlen[J].Journal für die Reineund Angewandte Mathematik,1850,40:209–227）。1975年，Voronoi图作为一种数据结构引入到计算机领域，成为计算几何领域研究热点之一。基于Voronoi图的插值方法称为自然邻点插值(Natrual Neighbor Interpolation,NNI)由Sibson于1980年提出（见参考文献[2]：Sibson R.A vector identity for the Dirichlet tessellation[A].Mathematical Proceedigns of the CambridgePhilosophical Scociety 1980[C],87(1):151-155，Bowyer A.Computing Dirichlet Tessellations[J].Computing Journal,1981,24(2):162-166），他发现插值点的位置向量可由其自然邻点的位置向量的凸组合来表示，系数由二阶Voronoi单元确定。自然邻点插值方法具有插值域稳定、平滑性高的特点，除了在自然邻点上是C⁰连续（曲线连续但一阶导数不连续），在Delaunay三角形外接圆边界上是C¹连续（曲线连续且一阶导数连续），在定义域其它地方都是C^∞连续（曲线连续且任意阶导数连续）。在含有约束条件的数据域中，插值结果应满足约束条件，这时应进行约束域插值。

海底地形数据除了水深点之外还有海岸线、岛屿等约束条件，这些约束条件在插值计算过程中需要保持，然而常规Voronoi图定义在凸区域上，对于含有非凸边界条件的区域自然邻点插值不适用。J.Yvonnet等提出约束自然单元法(constrained natural element method,C-NEM，参见参考文献[3]：Yvonnet J,Ryckelynck D,Lorong P，Chinesta F.A new extension of thenatural element method non-convex and discontinuous problem:the constrained natural elementmethod(C-NEM)[J].International Journal for Numerical Methods in Engineering,2004,60(8):1451-1474)，计算基于约束Voronoi图的Sibson坐标，该方法默认约束Voronoi图已知，但约束Voronoi图实现过程复杂。

发明内容

本发明为了减轻约束自然邻点插值中计算约束Voronoi图的压力，把协调Delaunay三角化(Conforming Delaunay Triangulation)的方法融入自然邻点插值，提出一种新的应用于海底地形的数字高程模型的建立方法，所述的建立方法采用插值方法，即协调自然邻点插值(Conforming Natural Neighbor Interpolation,Conforming-NNI)，该插值方法以光滑的插值函数对地形函数进行逼近，具有计算区域固定、光滑度高、满足非凸边界条件、易于实施等优点。

本发明的一种应用于海底地形的数字高程模型的建立方法，主要采用插值方法，具体包括如下步骤:

步骤1、根据实际情况获得确定待插点X的初始影响域。

步骤2、判断待插点X是否落在初始影响域的约束线段集Г^I上，如果是则转步骤7执行，如果不是则进行步骤3。

步骤3、构建嵌入特征约束后符合Delaunay三角化特性的三角网。

步骤4、寻找插值点的准自然邻点。

步骤5、根据对偶法则，顺次连接△X P_Nj P_Nj+1的外心得到V(X)，构建二阶常规Voronoi图。

步骤6、构建二阶约束Voronoi图。

步骤7、确定插值点的属性值，建立海底地形的数字高程模型。

本发明的优点和积极效果在于：

(1)本发明采用基于Voronoi图的Sibson坐标，通过协调Delauanay三角化，保证了约束线段的对偶边为Voronoi边，使得由插值点二阶常规Voronoi单元构建二阶约束Voronoi单元过程产生影响的约束线段的端点只在插值点的准自然邻点范围内出现，再经过简单的几何运算就可得到插值点的二阶约束Voronoi单元。

(2)本方法可用于建立海底规则格网数字高程模型（Grid-DEM），海底地形对船舶航行安全具有重要意义，矢量电子海图是常用的航海设备，但它只能以水深点和等深线等二维信息表达海底地形，存在表达不直观的缺陷。利用矢量电子海图中包含的水深点和岛屿海岸线等二维信息，构建海底2.5维数字高程模型，可使可观测维数增加，有助于更直观地进行航线规划。

(3)Grid-DEM建立过程的核心环节是格网点高程的内插计算。地形插值有整体拟合法、局部拟合法和逐点内插法，其中逐点内插法由于内插效率较高已成为目前DEM生产中常采用的方法。采用本方法可使海图中的岛屿海岸线等特征线作为约束线在插值过程中保持，且插值曲面可获得全局C⁰、大部分区域C^∞的光滑性，最大程度的保持了曲面的原有特性。

附图说明

图1为本发明的应用于海底地形的数字高程模型的建立方法的步骤流程图；

图2A为本发明实施例中插值点的初始影响域的示意图；

图2B为图2A嵌入约束线段后建立的协调Delaunay三角化后的三角形网格示意图；

图2C为在图2B基础上构建的一阶Voronoi图和二阶常规Voronoi图；

图2D为在图2C基础上嵌入约束线段P₂P₄后得到的插值点的二阶约束Voronoi单元；

图2E为在图2D基础上嵌入约束线段P₁P₆后得到的插值点的二阶约束Voronoi单元;

图2F为本发明实施例中插值点完整的二阶约束Voronoi单元；

图3A为本发明测试实例二中采用不规则三角网插值的平面显示结果；

图3B为本发明测试实例二中采用本发明方法的平面显示结果；

图4A为本发明测试实例二中采用不规则三角网插值的2.5维显示结果；

图4B为本发明测试实例二中采用本发明方法的2.5维显示结果；

图5为实施例插值点属性值计算公式中参数表达含义的示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步的详细说明。

如图1所示，本发明的一种应用于海底地形的数字高程模型的建立方法，采用插值方法进行建模，具体步骤包括：构建嵌入特征约束后符合DT（Delaunay Triangulation）特性的三角网，寻找插值点的准自然邻点,构建二阶常规Voronoi单元，由二阶常规Voronoi单元构建二阶约束Voronoi单元，计算插值点的属性值。下面结合实施例对本发明的各个步骤进行详细的说明。

步骤1、根据实际情况获得确定待插点X的初始影响域。

给定二维实数空间中平面直线图（Planar Straight Line Graph,简称PSLG）Ω(P,Г)。其中Г={Г₁,…,Г_M}为约束线段集，P={P₁,P₂,…,P_N}表示内点集和约束线段端点集组成的离散节点集。X为待插点。

确定以待插点X为中心，r倍节点平均边长为边长的正方形区域，把正方形区域内的节点和与正方形交集不为空的约束线段作为待插点X的初始影响域Ω(P^I,Г^I)，P^I为正方形区域内的节点和约束线段端点的合集，Г^I为约束线段集。r设置太大会影响搜索效率，太小可能会漏掉准自然邻点，r的大小应根据离散点的平均分布密度而定，最好使落在正方形区域内的节点为20～30个。

步骤2、判断待插点X是否落在初始影响域的约束线段集Г^I上，如果是则转步骤7执行，如果不是则进行步骤3。

步骤3、构建嵌入特征约束后符合Delaunay三角化特性的三角网。

对初始影响域中的P^I进行标准DT剖分，首先构建不考虑约束线段的无约束Delaunay三角网，然后利用端点三角形外接圆法构建协调Delaunay三角化后的三角形网格。

根据文献“地磁匹配用于水下载体导航的初步分析”（郝燕玲等，2008）提供的协调Delaunay三角化算法对没有被DT网格包含的约束线段进行细分，构建Conforming-DT网格。细分点作为Steiner点P_S，P_X=P^I∪P_S，Г_X为对Г^I细分后的约束线段集，此时数据域（使插值点X满足delaunay条件的点和约束边）变为Ω(P_X,Г_X)，满足DT(P_X)=CDT(P_X,Г_X)。所述的Steiner点是指协调Delaunay三角化过程产生的附加点。

目前，现有的协调Delaunay三角化方法中，都是最优化问题所对应的决策解法，还没有针对附加点最优分布的理论求法。现有的协调Delaunay三角化方法可分为两类：一类是先对特征线段按一定规则细分，然后进行DT（Delaunay Triangulation）剖分，代表是Boissonnat（参考文献[4]：Boissonnat J D.Shape Reconstruction from Planar Cross Sections[J].ComputerVision,Graphics,and Image Processing,1988,44(1):1-29）、Faugeras（参考文献[5]：Faugeras O D,Le Bras-Mehlman E,Bossionnat J D.Representing Stereo Data with the Delaunay Triangulation[J].Artifical Intelligence,1990,44(1-2):41-87）和Edelsbrunner（参考文献[6]：Edelsbrunner H,TiowSeng Tan.An Upper Bound for Conforming Delaunay Triangulations[J].Discrete & ComputationalGeometry,1993,10(2):197-213）算法；另一类是先进行无约束DT剖分，然后按照一定规则在未被包含的特征线段上插入附加点，代表是Boissonnat（参考文献[4]）、卢朝阳（参考文献[7]：卢朝阳,吴成柯,周幸妮.满足全局Delaunay特性的带特征约束的散乱数据最优三角剖分[J].计算机学报,1997,20(2):118-124）、Sapidis（参考文献[8]：Sapidis N,Perucchio R.DelaunayTriangulation of Arbitrarily Shaped Planar Domains[J].Computer Aided Geometric Design,1991,8(6):421-437）、Tsai（参考文献[9]：Tsai V J D.Delaunay Triangulations in TIN Creation:anOverview and Linear Time Algorithm[J].International Journal of Geographical InformationSystems,1993,7(6):501-524）、易法令（参考文献[10]：易法令,韩德志.带特征线约束的Delaunay三角剖分最优算法的研究及实现[J].计算机工程,2001,27(6):32-34）和梅承力（参考文献[11]：梅承力,肖高逾,周源华.一种新的带特征约束的Delaunay三角剖分算法[J].电子学报,2001,29(7):895-898）算法。受Sapidis、易法令和梅承力算法的启发，本发明提出一种改进的协调Delaunay三角化方法：端点三角形外接圆法(Endpoint Triangle′s Circumcircle Method,简称ETCM)，该方法能够保证剖分结果的稳定性，可以有效提高插入附加点的速度。

定义特征线段的一个端点三角形是仅以特征线段的一个端点为顶点，且一条边与特征线段相交的三角形。

设s[P_h P_e]是一条没有被DT网格包含的特征线段，△_h是P_h的端点三角形，△_e是P_e的端点三角形，Cir(△_h)和Cir(△_e)分别为△_h和△_e的外接圆。

端点三角形外接圆法具体是：设s[P_hP_e]是一条没有被DT网格包含的特征线段，P_h、P_e分别为线段的两个端点，△_h是P_h的端点三角形，△_e是P_e的端点三角形，Cir(△_h)和Cir(△_e)分别为△_h和△_e的外接圆。设J_h=Cir(△_h)∩s[P_hP_e]，即线段P_hP_e包含在△_h的外接圆内的部分，J_e=Cir(△_e)∩s[P_hP_e]，即线段P_hP_e包含在△_e的外接圆内的部分，∩表示交集，如果

则比较线段J_h和J_e长度的大小。若线段J_h的长度小于线段J_e的长度，则取线段J_e作为可嵌入的部分,线段J_h所对应的交点（外接圆Cir(△_h)与特征线段s[P_hP_e]的交点）作为Steiner点（反之取线段J_h作为可嵌入的部分,线段J_e所对应的交点（外接圆Cir(△_e)与特征线段s[P_hP_e]的交点）作为Steiner点），特征线段去掉可嵌入部分后剩余的部分继续使用上面的方法，直到外接圆与特征线段的交点形成的线段J_h′和J_e′满足此时取J_h′∩J_e′的中点作为Steiner点；若当前处理的特征线段影响到已嵌入的特征线段，则把受影响的特征线段按上面端点三角形外接圆法重新嵌入。

如图2A，设插值点X的初始影响域为Ω(P^I,Г^I)，其中P^I={P₁,P₂,P₃,P₄,P₅,P₆,P₇}为初始影响域内的节点和约束线段端点的集合，约束线段Г^I={P₂ P₄,P₁ P₆}以黑色实线表示。首先构建不考虑约束线段的无约束Delaunay三角网，然后利用端点三角形外接圆法构建协调Delaunay三角化后的三角形网格。如图2B为图2A嵌入约束线段P₂P₄,P₁P₆后建立的协调Delaunay三角化后的三角形网格，节点P₈为Steiner点。

端点三角形外接圆法与现有方法相比的优点是能够保证剖分结果的稳定性，可以有效提高插入附加点的速度。保证了不论以线段的哪个端点为起点，附加点的位置是唯一的；并且在插入Steiner点的过程中，迅速找到包含附加点的三角形或四边形是减少计算量的关键，鉴于附加点在特征线段的影响域内，ETCM一次性找到特征线段的影响域，然后在影响域中寻找，不必每次都要在全局进行寻找，减少了查找范围，可以快速完成附加点的插入。相对Sapidis算法中需要寻找特征线段剩余部分的端点三角形，一次性找到特征线段的影响域并没有增加计算量。

步骤4、寻找插值点的准自然邻点。

加入附加点P₈后的节点集合记为P_X=P^I∪P₈，Г_X为对Г^I细分后的约束线段集，Г_X={P₂ P₄,P₁ P₈,P₈ P₆}。此时数据域变为Ω(P_X,Г_X)。

将在P_X中距点X最近的节点P₁作为X的第一个准自然邻点P_N1。若加入附加点后的节点集合P_X中存在节点位于向量X P_N1的左侧，则设定寻找方向为左，若加入附加点后的节点集合P_X中不存在节点位于向量X P_N1的左侧，则设定寻找方向为右。如图2B所示，由于P_X中存在节点（P₂,P₃）位于向量X P_N1的左侧，因此寻找方向为左。以X与刚找到的准自然邻点P_Ni构成的向量XP_Ni为基边，把P_X中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点P_Ni+1。按照上述方法不断寻找下一个准自然邻点P_Ni+1，直到新找到的下一个准自然邻点P_Ni+1为P_N1时停止，N=i，{P_N1,P_N2,…,P_Nn}构成X的准自然邻点集P_N，P_N中的节点顺次相连构成对X的闭包。以X与第一个准自然邻点P_N1构成的向量X P_N1为基边，把P_X中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点P_N2，如图2B中的P₂点。所述的空外接圆准则为：在任意一个三角形的外接圆范围内（除圆周外）不包含点集P_X中的任何点。然后以向量X P_N2为基边，按照同样的方法，寻找准自然邻点P_N3。以此类推，直到新找到的下一个准自然邻点为P_N1时停止寻找。图2B中，记P_N={P₁,P₂,P₃,P₄,P₈}为X的准自然邻点集。节点P₁,P₂,P₃,P₄,P₈顺次相连构成对X的闭包Cell(P_N)。

步骤5、根据对偶法则，顺次连接△X P_Nj P_Nj+1的外心得到V(X)，构建二阶常规Voronoi图。

根据插入X前三角网格的拓扑关系，按对偶法则构建P^I的一阶Voronoi图V(P_Ni)(i=1～n)，P_Ni表示准自然邻点集中第i个节点，n为准自然邻点集的节点个数，如图2C中的细虚线所示。顺次连接X与Cell(P_N)上两个相邻节点构成的三角形△X P_Nj P_Nj+1的外心得到X的一阶常规Voronoi单元V(X),如图2C中的粗虚线所示。V(X)与V(P_Ni)重合得到插值点X的二阶常规Voronoi图V(X，P_Ni)。图2C中，V(X,P₁)=deh，V(X,P₂)=eafgh，V(X,P₃)=abf，V(X,P₄)=bcgf，V(X,P₈)=cdhg。

步骤6、构建二阶约束Voronoi图。

本发明中提出的由二阶常规Voronoi图构建二阶约束Voronoi图的方法描述如下：

给定待插值点的二阶常规Voronoi图和一条约束线段Γ_Xk，设Dual(P_iP_j)表示在Conforming-DT网格的对偶Voronoi图中，节点连线P_iP_j的对偶Voronoi边；Cell(P_N)是P_N中的节点顺次相连形成对X的闭包。对于任意指定的约束线段Γ_Xk，如果Γ_Xk的两个端点是P_N中相邻节点，则删去V(X,P_ej)(j=1,2）顶点中与X处于Γ_Xk异侧的，然后Γ_Xk与Dual(XP_ej)的交点和Γ_Xk的平分点，这两个新产生的点与V(X,P_ej)中剩余的顶点构成V^C(X,P_ej)；否则，如果Γ_Xk的两个端点不是P_N中相邻节点，则对于P_N中的任意节点，若P_Ni不可与X互视并且V(X,P_Ni)与Γ_Xk不相交，令V^C(X,P_Ni)=0，且Dual(XP_ej)靠近Γ_Xk的一个端点由Dual(XP_ej)所在直线与E_Xk的交点代替，这个新点与V(X,P_Ni)中剩余点构成V^C(X,P_Ni)；若P_Ni可与X互视并且V(X,P_Ni)与Γ_Xk不相交，令V^C(X,P_Ni)=V(X,P_Ni)；若P_Ni不可与X互视且V(X,P_Ni)与Γ_Xk相交，令V^C(X,P_Ni)=0，Dual(XP_ej)靠近Γ_Xk的一个端点由Dual(XP_ej)所在直线与Γ_Xk的交点代替，Dual(Γ_Xk)靠近Γ_Xk的一个端点由Γ_Xk的平分点代替，这两个新点与V(X,P_ej)中剩余顶点构成V^C(X,P_ej)。重复以上过程，最终得到插值点嵌入Γ_Xk后的二阶约束Voronoi单元。

逐条加入Г_X中的约束线段，构建插值点的二阶约束Voronoi单元。

由二阶常规Voronoi单元构建二阶约束Voronoi单元的算法实现流程如下：

图2D为嵌入约束线段P₂P₄后得到的插值点的二阶约束Voronoi单元，Γ_X1=P₂P₄，V^C(X,P₁)=deh，V^C(X,P₂)=epogh，V^C(X,P₃)=0，V^C(X,P₄)=oqcg，V^C(X,P₈)=cdhg。

图2E为嵌入约束线段P₁P₆后得到的插值点的二阶约束Voronoi单元。Γ_X2=P₁P₈，V^C(X,P₁)=ehwv，V^C(X,P₂)=epogh，V^C(X,P₃)=0，V^C(X,P₄)=oqcg，V^C(X,P₈)=cywhg。Γ_X3=P₈P₆的两个端点不是P_N中的点所以不会对X的二阶约束Voronoi单元产生影响。

图2F为插值点X完整的二阶约束Voronoi单元。V^C(X,P₁)=ehwv，V^C(X,P₂)=epogh，V^C(X,P₃)=0，V^C(X,P₄)=oqcg，V^C(X,P₈)=cywhg。

步骤7、确定插值点的属性值。

前面已经建立好了插值点的二阶约束Voronoi单元V^C(X,P_i)，根据V^C(X,P_i)和P_i的属性值就可以确定插值点的属性值了。

根据式（1），二维约束域中协调自然邻点插值点X的形函数φ^C(X,P_Ni)为：

式中：

表示

空间中元素*的n维Lebesgue测度。

表示实数。

根据式（2），二维约束域中协调自然邻点插值点X的试探函数u^h(X)为：

$u^h\left(X\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\φ}}_i^C(X,P_{\mathrm{Ni}})\·u\left(P_{\mathrm{Ni}}\right)---\left(2\right)$

式中，P_Ni表示X的第i个准自然邻点，

表示第i个准自然邻点的形函数，u(P_Ni)表示点P_Ni的属性值。n为X准自然邻点的个数。

如果插值点落在某条约束线段上，由于约束线段不透明，相当于插值点落在数据域边界上，此时插值点的一阶约束Voronoi单元没有有形边界，约束线段端点的二阶约束Voronoi单元也没有有形边界，它们的二维测度都为无穷大，根据式（3）可推出此时插值点的其它可视自然邻点所对应的形函数趋近于零，插值点属性值是该约束线段两个端点属性值的线性比例之和。

$\left\{\begin{array}{c}\mathrm{\φ}(X,P_1)=\underset{L\→\∞}{\lim }\frac{L\·(1-\mathrm{\ξ})+2{\mathrm{\δ}}_1}{L+2{\mathrm{\δ}}_1+2{\mathrm{\δ}}_2+2{\mathrm{\δ}}_3}=1-\mathrm{\ξ}\\ \mathrm{\φ}(X,P_2)=\underset{L\→\∞}{\lim }\frac{L\·\left(\mathrm{\ξ}\right)+2{\mathrm{\δ}}_2}{L+2{\mathrm{\δ}}_1+2{\mathrm{\δ}}_2+2{\mathrm{\δ}}_3}=\mathrm{\ξ}\\ \mathrm{\φ}(X,P_3)=\underset{L\→\∞}{\lim }\frac{{2\mathrm{\δ}}_3}{L+{2\mathrm{\δ}}_1+{2\mathrm{\δ}}_2+{2\mathrm{\δ}}_3}=0\end{array}\right.---\left(3\right)$

u^h(X)=(1-ξ)·u^h(P₁)+ξ·u^h(P₂)                 （4）

ξ为插值点X与P₁点的距离。φ(X,P₁)、φ(X,P₂)、φ(X,P₃)分别表示点P1、P2、P3的属性值对插值点X的属性值的影响权值，

见图5，ξ为X到P₁的距离。δ₁为多边形cdea的面积、δ₂为多边形aefb的面积、δ₃为多边形cab的面积。L表示线段P1P2到点集区域边界的距离。

u^h(P₁)、u^h(P₂)分别表示节点P₁、P₂的试探函数。因为协调自然邻点插值计算基于约束Voronoi图的Sibson坐标，所以与自然邻点插值拥有相似的连续性，但由于插值点穿越约束线段前后可视自然邻点发生变化，导致插值结果的一阶导数在约束线段处不连续，所以协调自然邻点插值在节点处和约束线段上是C⁰连续，其它区域是C¹连续（其中大多数地方是C^∝连续）。

为进一步说明本发明提出的应用于海底地形的数字高程模型采用的插值方法的可行性和有效性，选取图号为C11010的实际电子海图数据作为测试实例二，以不规则三角网（triangulated irregular network,TIN）插值作为对比，阐明本发明方法的有益效果。

选取的测试用海图参数由1255个水深点、10094条约束线段和69*98个插值点组成。实验中Grid DEM（网格数字高程模型）的网格宽度设为100个图上单位距离。

实验结果：不规则三角网插值和协调自然邻点插值的平面显示与2.5维显示见图3A、3B和图4A、4B。采用不规则三角网插值对海底地形插值逼近的粗糙度为2.5439，坡度标准差为0.04°，水深标准差为17.6465m。而应用本发明方法对应得到的三项数据分别为：2.5926，0.0423°，17.6564m。

在协调自然邻点插值生成的Grid DEM中，岛屿和海岸线边界得到保持。由于已有的海底地形数据是散乱排列的离散点，所以无法从插值均方根误差的大小来分析插值算法的性能，分别从视觉和反应地形信息丰富程度的地形参数上对TIN插值和协调自然邻点插值进行比较。从视觉上看，协调自然邻点插值反应出TIN插值没有反映出的细节内容。从粗糙度、坡度标准差和水深标准差等反映地形信息丰富程度的地形参数上看，协调自然邻点插值得到的地形信息也比TIN插值得到的地形信息要丰富些。

Claims (6)

1.一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的建立方法包括如下步骤：

步骤1、根据实际情况获得确定待插点X的初始影响域：确定以待插点X为中心，r倍节点平均边长为边长的正方形区域，把正方形区域内的节点和与正方形交集不为空的约束线段作为待插点X的初始影响域Ω(P^I,Г^I)，P^I为正方形区域内的节点和约束线段端点的合集，Г^I为约束线段集；r的大小应根据离散点的平均分布密度而定；

步骤2、判断待插点X是否落在初始影响域的约束线段集Г^I上，如果是则转步骤7执行，如果不是则进行步骤3；

步骤3、构建嵌入特征约束后符合Delaunay三角化特性的三角网；

步骤4、寻找插值点的准自然邻点；

步骤5、根据对偶法则，顺次连接△X P_Nj P_Nj+1的外心得到V(X)，构建二阶常规Voronoi图；

步骤6、构建二阶约束Voronoi图；

步骤7、确定插值点的属性值，建立海底地形的数字高程模型。

2.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的r的设置使落在正方形区域内的节点为20～30个。

3.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的步骤3中采用端点三角形外接圆法，特征线段的一个端点三角形是仅以特征线段的一个端点为顶点，且一条边与特征线段相交的三角形；

设s[P_h P_e]是一条没有被DT网格包含的特征线段，△_h是P_h的端点三角形，△_e是P_e的端点三角形，Cir(△_h)和Cir(△_e)分别为△_h和△_e的外接圆；

端点三角形外接圆法具体是：设s[P_h P_e]是一条没有被DT网格包含的特征线段，P_h、P_e分别为线段的两个端点，△_h是P_h的端点三角形，△_e是P_e的端点三角形，Cir(△_h)和Cir(△_e)分别为△_h和△_e的外接圆；设J_h=Cir(△_h)∩s[P_hP_e]，即线段P_hP_e包含在△_h的外接圆内的部分，J_e=Cir(△_e)∩s[P_hP_e]，即线段P_hP_e包含在△_e的外接圆内的部分，∩表示交集，如果

则比较线段J_h和J_e长度的大小：若线段J_h的长度小于线段J_e的长度，则取线段J_e作为可嵌入的部分,线段J_h所对应的外接圆Cir(△_h)与特征线段s[P_hP_e]的交点作为Steiner点；反之取线段J_h作为可嵌入的部分,线段J_e所对应的外接圆Cir(△_e)与特征线段s[P_hP_e]的交点作为Steiner点，特征线段去掉可嵌入部分后剩余的部分继续使用上面的方法，直到

此时取J_h′∩J_e′的中点作为Steiner点；若当前处理的特征线段影响到已嵌入的特征线段，则把受影响的特征线段按上面端点三角形外接圆法重新嵌入。

4.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的步骤4具体为：

将在节点集合P_X中距点X最近的节点P₁作为X的第一个准自然邻点P_N1；若加入附加点后的节点集合P_X中存在节点位于向量X P_N1的左侧，则设定寻找方向为左，若加入附加点后的节点集合P_X中不存在节点位于向量X P_N1的左侧，则设定寻找方向为右；以X与刚找到的准自然邻点P_Ni构成的向量XP_Ni为基边，把P_X中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点P_Ni+1；

按照上述方法不断寻找下一个准自然邻点P_Ni+1直到新找到的P_Ni+1为P_N1时停止。

5.根据权利要求4所述的一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的空外接圆准则为：在任意一个三角形的外接圆范围内不包含点集P_X中的任何点。

6.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法，其特征在于：所述的步骤6具体为：

给定待插值点的二阶常规Voronoi单元和一条约束线段Γ_Xk，设Dual(P_iP_j)表示在Conforming-DT网格的对偶Voronoi图中，节点连线P_iP_j的对偶Voronoi边；Cell(P_N)是P_N中的节点顺次相连形成对X的闭包；对于任意指定的约束线段Γ_Xk，如果Γ_Xk的两个端点是P_N中相邻节点，则删去V(X,P_ej)顶点中与X处于Γ_Xk异侧的，j=1,2；然后Γ_Xk与Dual(XP_ej)的交点和Γ_Xk的平分点，这两个新产生的点与V(X,P_ej)中剩余的顶点构成V^C(X,P_ej)；否则，如果Γ_Xk的两个端点不是P_N中相邻节点，则对于P_N中的任意节点，若P_Ni不可与X互视并且V(X,P_Ni)与Γ_Xk不相交，令V^C(X,P_Ni)=0，且Dual(XP_ej)靠近Γ_Xk的一个端点由Dual(XP_ej)所在直线与E_Xk的交点代替，这个新点与V(X,P_Ni)中剩余点构成V^C(X,P_Ni)；若P_Ni可与X互视并且V(X,P_Ni)与Γ_Xk不相交，令V^C(X,P_Ni)=V(X,P_Ni)；若P_Ni不可与X互视且V(X,P_Ni)与Γ_Xk相交，令V^C(X,P_Ni)=0，Dual(XP_ej)靠近Γ_Xk的一个端点由Dual(XP_ej)所在直线与Γ_Xk的交点代替，Dual(Γ_Xk)靠近Γ_Xk的一个端点由Γ_Xk的平分点代替，这两个新点与V(X,P_ej)中剩余顶点构成V^C(X,P_ej)；重复以上过程，最终得到插值点嵌入Γ_Xk后的二阶约束Voronoi单元。

Images