CN102831645A - 一种应用于海底地形的数字高程模型的建立方法 - Google Patents
一种应用于海底地形的数字高程模型的建立方法 Download PDFInfo
- Publication number
- CN102831645A CN102831645A CN2012102485347A CN201210248534A CN102831645A CN 102831645 A CN102831645 A CN 102831645A CN 2012102485347 A CN2012102485347 A CN 2012102485347A CN 201210248534 A CN201210248534 A CN 201210248534A CN 102831645 A CN102831645 A CN 102831645A
- Authority
- CN
- China
- Prior art keywords
- point
- line segment
- interpolation
- end points
- node
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Processing Or Creating Images (AREA)
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单元确定。自然邻点插值方法具有插值域稳定、平滑性高的特点,除了在自然邻点上是C0连续(曲线连续但一阶导数不连续),在Delaunay三角形外接圆边界上是C1连续(曲线连续且一阶导数连续),在定义域其它地方都是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 PNj PNj+1的外心得到V(X),构建二阶常规Voronoi图。
步骤6、构建二阶约束Voronoi图。
步骤7、确定插值点的属性值,建立海底地形的数字高程模型。
本发明的优点和积极效果在于:
(1)本发明采用基于Voronoi图的Sibson坐标,通过协调Delauanay三角化,保证了约束线段的对偶边为Voronoi边,使得由插值点二阶常规Voronoi单元构建二阶约束Voronoi单元过程产生影响的约束线段的端点只在插值点的准自然邻点范围内出现,再经过简单的几何运算就可得到插值点的二阶约束Voronoi单元。
(2)本方法可用于建立海底规则格网数字高程模型(Grid-DEM),海底地形对船舶航行安全具有重要意义,矢量电子海图是常用的航海设备,但它只能以水深点和等深线等二维信息表达海底地形,存在表达不直观的缺陷。利用矢量电子海图中包含的水深点和岛屿海岸线等二维信息,构建海底2.5维数字高程模型,可使可观测维数增加,有助于更直观地进行航线规划。
(3)Grid-DEM建立过程的核心环节是格网点高程的内插计算。地形插值有整体拟合法、局部拟合法和逐点内插法,其中逐点内插法由于内插效率较高已成为目前DEM生产中常采用的方法。采用本方法可使海图中的岛屿海岸线等特征线作为约束线在插值过程中保持,且插值曲面可获得全局C0、大部分区域C∞的光滑性,最大程度的保持了曲面的原有特性。
附图说明
图1为本发明的应用于海底地形的数字高程模型的建立方法的步骤流程图;
图2A为本发明实施例中插值点的初始影响域的示意图;
图2B为图2A嵌入约束线段后建立的协调Delaunay三角化后的三角形网格示意图;
图2C为在图2B基础上构建的一阶Voronoi图和二阶常规Voronoi图;
图2D为在图2C基础上嵌入约束线段P2P4后得到的插值点的二阶约束Voronoi单元;
图2E为在图2D基础上嵌入约束线段P1P6后得到的插值点的二阶约束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,Г)。其中Г={Г1,…,ГM}为约束线段集,P={P1,P2,…,PN}表示内点集和约束线段端点集组成的离散节点集。X为待插点。
确定以待插点X为中心,r倍节点平均边长为边长的正方形区域,把正方形区域内的节点和与正方形交集不为空的约束线段作为待插点X的初始影响域Ω(PI,ГI),PI为正方形区域内的节点和约束线段端点的合集,ГI为约束线段集。r设置太大会影响搜索效率,太小可能会漏掉准自然邻点,r的大小应根据离散点的平均分布密度而定,最好使落在正方形区域内的节点为20~30个。
步骤2、判断待插点X是否落在初始影响域的约束线段集ГI上,如果是则转步骤7执行,如果不是则进行步骤3。
步骤3、构建嵌入特征约束后符合Delaunay三角化特性的三角网。
对初始影响域中的PI进行标准DT剖分,首先构建不考虑约束线段的无约束Delaunay三角网,然后利用端点三角形外接圆法构建协调Delaunay三角化后的三角形网格。
根据文献“地磁匹配用于水下载体导航的初步分析”(郝燕玲等,2008)提供的协调Delaunay三角化算法对没有被DT网格包含的约束线段进行细分,构建Conforming-DT网格。细分点作为Steiner点PS,PX=PI∪PS,ГX为对ГI细分后的约束线段集,此时数据域(使插值点X满足delaunay条件的点和约束边)变为Ω(PX,ГX),满足DT(PX)=CDT(PX,Г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[Ph Pe]是一条没有被DT网格包含的特征线段,△h是Ph的端点三角形,△e是Pe的端点三角形,Cir(△h)和Cir(△e)分别为△h和△e的外接圆。
端点三角形外接圆法具体是:设s[PhPe]是一条没有被DT网格包含的特征线段,Ph、Pe分别为线段的两个端点,△h是Ph的端点三角形,△e是Pe的端点三角形,Cir(△h)和Cir(△e)分别为△h和△e的外接圆。设Jh=Cir(△h)∩s[PhPe],即线段PhPe包含在△h的外接圆内的部分,Je=Cir(△e)∩s[PhPe],即线段PhPe包含在△e的外接圆内的部分,∩表示交集,如果则比较线段Jh和Je长度的大小。若线段Jh的长度小于线段Je的长度,则取线段Je作为可嵌入的部分,线段Jh所对应的交点(外接圆Cir(△h)与特征线段s[PhPe]的交点)作为Steiner点(反之取线段Jh作为可嵌入的部分,线段Je所对应的交点(外接圆Cir(△e)与特征线段s[PhPe]的交点)作为Steiner点),特征线段去掉可嵌入部分后剩余的部分继续使用上面的方法,直到外接圆与特征线段的交点形成的线段Jh′和Je′满足此时取Jh′∩Je′的中点作为Steiner点;若当前处理的特征线段影响到已嵌入的特征线段,则把受影响的特征线段按上面端点三角形外接圆法重新嵌入。
如图2A,设插值点X的初始影响域为Ω(PI,ГI),其中PI={P1,P2,P3,P4,P5,P6,P7}为初始影响域内的节点和约束线段端点的集合,约束线段ГI={P2 P4,P1 P6}以黑色实线表示。首先构建不考虑约束线段的无约束Delaunay三角网,然后利用端点三角形外接圆法构建协调Delaunay三角化后的三角形网格。如图2B为图2A嵌入约束线段P2P4,P1P6后建立的协调Delaunay三角化后的三角形网格,节点P8为Steiner点。
端点三角形外接圆法与现有方法相比的优点是能够保证剖分结果的稳定性,可以有效提高***附加点的速度。保证了不论以线段的哪个端点为起点,附加点的位置是唯一的;并且在***Steiner点的过程中,迅速找到包含附加点的三角形或四边形是减少计算量的关键,鉴于附加点在特征线段的影响域内,ETCM一次性找到特征线段的影响域,然后在影响域中寻找,不必每次都要在全局进行寻找,减少了查找范围,可以快速完成附加点的***。相对Sapidis算法中需要寻找特征线段剩余部分的端点三角形,一次性找到特征线段的影响域并没有增加计算量。
步骤4、寻找插值点的准自然邻点。
加入附加点P8后的节点集合记为PX=PI∪P8,ГX为对ГI细分后的约束线段集,ГX={P2 P4,P1 P8,P8 P6}。此时数据域变为Ω(PX,ГX)。
将在PX中距点X最近的节点P1作为X的第一个准自然邻点PN1。若加入附加点后的节点集合PX中存在节点位于向量X PN1的左侧,则设定寻找方向为左,若加入附加点后的节点集合PX中不存在节点位于向量X PN1的左侧,则设定寻找方向为右。如图2B所示,由于PX中存在节点(P2,P3)位于向量X PN1的左侧,因此寻找方向为左。以X与刚找到的准自然邻点PNi构成的向量XPNi为基边,把PX中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点PNi+1。按照上述方法不断寻找下一个准自然邻点PNi+1,直到新找到的下一个准自然邻点PNi+1为PN1时停止,N=i,{PN1,PN2,…,PNn}构成X的准自然邻点集PN,PN中的节点顺次相连构成对X的闭包。以X与第一个准自然邻点PN1构成的向量X PN1为基边,把PX中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点PN2,如图2B中的P2点。所述的空外接圆准则为:在任意一个三角形的外接圆范围内(除圆周外)不包含点集PX中的任何点。然后以向量X PN2为基边,按照同样的方法,寻找准自然邻点PN3。以此类推,直到新找到的下一个准自然邻点为PN1时停止寻找。图2B中,记PN={P1,P2,P3,P4,P8}为X的准自然邻点集。节点P1,P2,P3,P4,P8顺次相连构成对X的闭包Cell(PN)。
步骤5、根据对偶法则,顺次连接△X PNj PNj+1的外心得到V(X),构建二阶常规Voronoi图。
根据***X前三角网格的拓扑关系,按对偶法则构建PI的一阶Voronoi图V(PNi)(i=1~n),PNi表示准自然邻点集中第i个节点,n为准自然邻点集的节点个数,如图2C中的细虚线所示。顺次连接X与Cell(PN)上两个相邻节点构成的三角形△X PNj PNj+1的外心得到X的一阶常规Voronoi单元V(X),如图2C中的粗虚线所示。V(X)与V(PNi)重合得到插值点X的二阶常规Voronoi图V(X,PNi)。图2C中,V(X,P1)=deh,V(X,P2)=eafgh,V(X,P3)=abf,V(X,P4)=bcgf,V(X,P8)=cdhg。
步骤6、构建二阶约束Voronoi图。
本发明中提出的由二阶常规Voronoi图构建二阶约束Voronoi图的方法描述如下:
给定待插值点的二阶常规Voronoi图和一条约束线段ΓXk,设Dual(PiPj)表示在Conforming-DT网格的对偶Voronoi图中,节点连线PiPj的对偶Voronoi边;Cell(PN)是PN中的节点顺次相连形成对X的闭包。对于任意指定的约束线段ΓXk,如果ΓXk的两个端点是PN中相邻节点,则删去V(X,Pej)(j=1,2)顶点中与X处于ΓXk异侧的,然后ΓXk与Dual(XPej)的交点和ΓXk的平分点,这两个新产生的点与V(X,Pej)中剩余的顶点构成VC(X,Pej);否则,如果ΓXk的两个端点不是PN中相邻节点,则对于PN中的任意节点,若PNi不可与X互视并且V(X,PNi)与ΓXk不相交,令VC(X,PNi)=0,且Dual(XPej)靠近ΓXk的一个端点由Dual(XPej)所在直线与EXk的交点代替,这个新点与V(X,PNi)中剩余点构成VC(X,PNi);若PNi可与X互视并且V(X,PNi)与ΓXk不相交,令VC(X,PNi)=V(X,PNi);若PNi不可与X互视且V(X,PNi)与ΓXk相交,令VC(X,PNi)=0,Dual(XPej)靠近ΓXk的一个端点由Dual(XPej)所在直线与ΓXk的交点代替,Dual(ΓXk)靠近ΓXk的一个端点由ΓXk的平分点代替,这两个新点与V(X,Pej)中剩余顶点构成VC(X,Pej)。重复以上过程,最终得到插值点嵌入ΓXk后的二阶约束Voronoi单元。
逐条加入ГX中的约束线段,构建插值点的二阶约束Voronoi单元。
由二阶常规Voronoi单元构建二阶约束Voronoi单元的算法实现流程如下:
图2D为嵌入约束线段P2P4后得到的插值点的二阶约束Voronoi单元,ΓX1=P2P4,VC(X,P1)=deh,VC(X,P2)=epogh,VC(X,P3)=0,VC(X,P4)=oqcg,VC(X,P8)=cdhg。
图2E为嵌入约束线段P1P6后得到的插值点的二阶约束Voronoi单元。ΓX2=P1P8,VC(X,P1)=ehwv,VC(X,P2)=epogh,VC(X,P3)=0,VC(X,P4)=oqcg,VC(X,P8)=cywhg。ΓX3=P8P6的两个端点不是PN中的点所以不会对X的二阶约束Voronoi单元产生影响。
图2F为插值点X完整的二阶约束Voronoi单元。VC(X,P1)=ehwv,VC(X,P2)=epogh,VC(X,P3)=0,VC(X,P4)=oqcg,VC(X,P8)=cywhg。
步骤7、确定插值点的属性值。
前面已经建立好了插值点的二阶约束Voronoi单元VC(X,Pi),根据VC(X,Pi)和Pi的属性值就可以确定插值点的属性值了。
根据式(1),二维约束域中协调自然邻点插值点X的形函数φC(X,PNi)为:
根据式(2),二维约束域中协调自然邻点插值点X的试探函数uh(X)为:
如果插值点落在某条约束线段上,由于约束线段不透明,相当于插值点落在数据域边界上,此时插值点的一阶约束Voronoi单元没有有形边界,约束线段端点的二阶约束Voronoi单元也没有有形边界,它们的二维测度都为无穷大,根据式(3)可推出此时插值点的其它可视自然邻点所对应的形函数趋近于零,插值点属性值是该约束线段两个端点属性值的线性比例之和。
uh(X)=(1-ξ)·uh(P1)+ξ·uh(P2) (4)
ξ为插值点X与P1点的距离。φ(X,P1)、φ(X,P2)、φ(X,P3)分别表示点P1、P2、P3的属性值对插值点X的属性值的影响权值,
见图5,ξ为X到P1的距离。δ1为多边形cdea的面积、δ2为多边形aefb的面积、δ3为多边形cab的面积。L表示线段P1P2到点集区域边界的距离。
uh(P1)、uh(P2)分别表示节点P1、P2的试探函数。因为协调自然邻点插值计算基于约束Voronoi图的Sibson坐标,所以与自然邻点插值拥有相似的连续性,但由于插值点穿越约束线段前后可视自然邻点发生变化,导致插值结果的一阶导数在约束线段处不连续,所以协调自然邻点插值在节点处和约束线段上是C0连续,其它区域是C1连续(其中大多数地方是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的初始影响域Ω(PI,ГI),PI为正方形区域内的节点和约束线段端点的合集,ГI为约束线段集;r的大小应根据离散点的平均分布密度而定;
步骤2、判断待插点X是否落在初始影响域的约束线段集ГI上,如果是则转步骤7执行,如果不是则进行步骤3;
步骤3、构建嵌入特征约束后符合Delaunay三角化特性的三角网;
步骤4、寻找插值点的准自然邻点;
步骤5、根据对偶法则,顺次连接△X PNj PNj+1的外心得到V(X),构建二阶常规Voronoi图;
步骤6、构建二阶约束Voronoi图;
步骤7、确定插值点的属性值,建立海底地形的数字高程模型。
2.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法,其特征在于:所述的r的设置使落在正方形区域内的节点为20~30个。
3.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法,其特征在于:所述的步骤3中采用端点三角形外接圆法,特征线段的一个端点三角形是仅以特征线段的一个端点为顶点,且一条边与特征线段相交的三角形;
设s[Ph Pe]是一条没有被DT网格包含的特征线段,△h是Ph的端点三角形,△e是Pe的端点三角形,Cir(△h)和Cir(△e)分别为△h和△e的外接圆;
端点三角形外接圆法具体是:设s[Ph Pe]是一条没有被DT网格包含的特征线段,Ph、Pe分别为线段的两个端点,△h是Ph的端点三角形,△e是Pe的端点三角形,Cir(△h)和Cir(△e)分别为△h和△e的外接圆;设Jh=Cir(△h)∩s[PhPe],即线段PhPe包含在△h的外接圆内的部分,Je=Cir(△e)∩s[PhPe],即线段PhPe包含在△e的外接圆内的部分,∩表示交集,如果则比较线段Jh和Je长度的大小:若线段Jh的长度小于线段Je的长度,则取线段Je作为可嵌入的部分,线段Jh所对应的外接圆Cir(△h)与特征线段s[PhPe]的交点作为Steiner点;反之取线段Jh作为可嵌入的部分,线段Je所对应的外接圆Cir(△e)与特征线段s[PhPe]的交点作为Steiner点,特征线段去掉可嵌入部分后剩余的部分继续使用上面的方法,直到此时取Jh′∩Je′的中点作为Steiner点;若当前处理的特征线段影响到已嵌入的特征线段,则把受影响的特征线段按上面端点三角形外接圆法重新嵌入。
4.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法,其特征在于:所述的步骤4具体为:
将在节点集合PX中距点X最近的节点P1作为X的第一个准自然邻点PN1;若加入附加点后的节点集合PX中存在节点位于向量X PN1的左侧,则设定寻找方向为左,若加入附加点后的节点集合PX中不存在节点位于向量X PN1的左侧,则设定寻找方向为右;以X与刚找到的准自然邻点PNi构成的向量XPNi为基边,把PX中位于基边寻找方向一侧且满足空外接圆准则的一个点作为X的下一个准自然邻点PNi+1;
按照上述方法不断寻找下一个准自然邻点PNi+1直到新找到的PNi+1为PN1时停止。
5.根据权利要求4所述的一种应用于海底地形的数字高程模型的建立方法,其特征在于:所述的空外接圆准则为:在任意一个三角形的外接圆范围内不包含点集PX中的任何点。
6.根据权利要求1所述的一种应用于海底地形的数字高程模型的建立方法,其特征在于:所述的步骤6具体为:
给定待插值点的二阶常规Voronoi单元和一条约束线段ΓXk,设Dual(PiPj)表示在Conforming-DT网格的对偶Voronoi图中,节点连线PiPj的对偶Voronoi边;Cell(PN)是PN中的节点顺次相连形成对X的闭包;对于任意指定的约束线段ΓXk,如果ΓXk的两个端点是PN中相邻节点,则删去V(X,Pej)顶点中与X处于ΓXk异侧的,j=1,2;然后ΓXk与Dual(XPej)的交点和ΓXk的平分点,这两个新产生的点与V(X,Pej)中剩余的顶点构成VC(X,Pej);否则,如果ΓXk的两个端点不是PN中相邻节点,则对于PN中的任意节点,若PNi不可与X互视并且V(X,PNi)与ΓXk不相交,令VC(X,PNi)=0,且Dual(XPej)靠近ΓXk的一个端点由Dual(XPej)所在直线与EXk的交点代替,这个新点与V(X,PNi)中剩余点构成VC(X,PNi);若PNi可与X互视并且V(X,PNi)与ΓXk不相交,令VC(X,PNi)=V(X,PNi);若PNi不可与X互视且V(X,PNi)与ΓXk相交,令VC(X,PNi)=0,Dual(XPej)靠近ΓXk的一个端点由Dual(XPej)所在直线与ΓXk的交点代替,Dual(ΓXk)靠近ΓXk的一个端点由ΓXk的平分点代替,这两个新点与V(X,Pej)中剩余顶点构成VC(X,Pej);重复以上过程,最终得到插值点嵌入ΓXk后的二阶约束Voronoi单元。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2012102485347A CN102831645A (zh) | 2012-07-18 | 2012-07-18 | 一种应用于海底地形的数字高程模型的建立方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2012102485347A CN102831645A (zh) | 2012-07-18 | 2012-07-18 | 一种应用于海底地形的数字高程模型的建立方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN102831645A true CN102831645A (zh) | 2012-12-19 |
Family
ID=47334758
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN2012102485347A Pending CN102831645A (zh) | 2012-07-18 | 2012-07-18 | 一种应用于海底地形的数字高程模型的建立方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102831645A (zh) |
Cited By (18)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103093410A (zh) * | 2013-01-06 | 2013-05-08 | 国家***第二海洋研究所 | 一种海底地形六维网格的测绘方法 |
CN103256914A (zh) * | 2012-12-27 | 2013-08-21 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淹没面积的方法及*** |
CN103278115A (zh) * | 2012-12-27 | 2013-09-04 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淤积量的方法及*** |
CN103824510A (zh) * | 2013-12-31 | 2014-05-28 | 厦门雅迅网络股份有限公司 | 一种基于Voronoi图的筛选电子地图点状要素的方法 |
CN103854302A (zh) * | 2013-12-23 | 2014-06-11 | 哈尔滨工程大学 | 一种多约束条件下的auv航行环境空间构建方法 |
CN103955932A (zh) * | 2014-05-05 | 2014-07-30 | 张立华 | 一种不漏水深浅点的海图水深自动选取方法 |
CN105973246A (zh) * | 2016-04-29 | 2016-09-28 | 海尔优家智能科技(北京)有限公司 | 一种地磁地图的绘制方法、装置及机器人 |
CN106485244A (zh) * | 2016-10-12 | 2017-03-08 | 上海联影医疗科技有限公司 | 采样方法及装置 |
CN106875486A (zh) * | 2017-02-22 | 2017-06-20 | 哈尔滨工程大学 | 一种基于节点信息量统计的多波束地形分块方法 |
CN107038308A (zh) * | 2017-04-18 | 2017-08-11 | 南京工程学院 | 一种基于线性内插的规则格网地形建模方法 |
CN107798339A (zh) * | 2017-09-29 | 2018-03-13 | 天津大学 | 一种基于泰森多边形图的多波束异常值数据处理算法 |
CN108681561A (zh) * | 2018-04-20 | 2018-10-19 | 贾帅东 | 一种航海图等深线平滑性的评估方法 |
CN108874932A (zh) * | 2018-05-31 | 2018-11-23 | 哈尔滨工程大学 | 一种基于改进的光线投射算法的海洋水声场三维可视化方法 |
CN109524990A (zh) * | 2018-12-03 | 2019-03-26 | 三峡大学 | 一种基于Voronoi图重心内插法的虚拟惯量配置方法 |
CN110118526A (zh) * | 2019-03-08 | 2019-08-13 | 浙江中海达空间信息技术有限公司 | 一种支持实时监测的船载砂石体积自动计算方法 |
CN110927792A (zh) * | 2019-10-24 | 2020-03-27 | 广州海洋地质调查局 | 一种地震构造图的建立方法及处理终端 |
CN113806951A (zh) * | 2021-09-23 | 2021-12-17 | 华侨大学 | 一种基于半边数据结构的自然邻近点搜索的弹性仿真方法 |
CN116127792A (zh) * | 2023-04-17 | 2023-05-16 | 北京世冠金洋科技发展有限公司 | 一种散乱数据的插值方法及装置 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20020035553A1 (en) * | 2000-07-20 | 2002-03-21 | Kim Seung Bum | Intelligent interpolation methods for automatic generation of an accurate digital elevation model |
CN102521882A (zh) * | 2011-12-05 | 2012-06-27 | 西北工业大学 | 基于离散高程和自适应混合加权得到海床地形数据的方法 |
-
2012
- 2012-07-18 CN CN2012102485347A patent/CN102831645A/zh active Pending
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20020035553A1 (en) * | 2000-07-20 | 2002-03-21 | Kim Seung Bum | Intelligent interpolation methods for automatic generation of an accurate digital elevation model |
CN102521882A (zh) * | 2011-12-05 | 2012-06-27 | 西北工业大学 | 基于离散高程和自适应混合加权得到海床地形数据的方法 |
Non-Patent Citations (3)
Title |
---|
GEORGE Y. LU ET AL.: "An adaptive inverse-distance weighting spatial interpolation technique", 《COMPUTERS & GEOSCIENCES》 * |
田峰敏: "水下地形导航模型求解与导航区初选策略研究", 《中国博士学位论文全文数据库 信息科技辑》 * |
田峰敏等: "一种建立海底格网数字高程模型的插值方法", 《中国航海》 * |
Cited By (28)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103256914A (zh) * | 2012-12-27 | 2013-08-21 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淹没面积的方法及*** |
CN103278115A (zh) * | 2012-12-27 | 2013-09-04 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淤积量的方法及*** |
CN103256914B (zh) * | 2012-12-27 | 2015-11-18 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淹没面积的方法及*** |
CN103278115B (zh) * | 2012-12-27 | 2016-06-08 | 北京地拓科技发展有限公司 | 一种基于dem计算淤地坝淤积量的方法及*** |
CN103093410A (zh) * | 2013-01-06 | 2013-05-08 | 国家***第二海洋研究所 | 一种海底地形六维网格的测绘方法 |
CN103854302A (zh) * | 2013-12-23 | 2014-06-11 | 哈尔滨工程大学 | 一种多约束条件下的auv航行环境空间构建方法 |
CN103854302B (zh) * | 2013-12-23 | 2016-09-28 | 哈尔滨工程大学 | 一种多约束条件下的auv航行环境空间构建方法 |
CN103824510A (zh) * | 2013-12-31 | 2014-05-28 | 厦门雅迅网络股份有限公司 | 一种基于Voronoi图的筛选电子地图点状要素的方法 |
CN103955932A (zh) * | 2014-05-05 | 2014-07-30 | 张立华 | 一种不漏水深浅点的海图水深自动选取方法 |
CN103955932B (zh) * | 2014-05-05 | 2016-08-24 | 张立华 | 一种不漏水深浅点的海图水深自动选取方法 |
CN105973246A (zh) * | 2016-04-29 | 2016-09-28 | 海尔优家智能科技(北京)有限公司 | 一种地磁地图的绘制方法、装置及机器人 |
CN106485244A (zh) * | 2016-10-12 | 2017-03-08 | 上海联影医疗科技有限公司 | 采样方法及装置 |
CN106875486B (zh) * | 2017-02-22 | 2020-04-07 | 哈尔滨工程大学 | 一种基于节点信息量统计的多波束地形分块方法 |
CN106875486A (zh) * | 2017-02-22 | 2017-06-20 | 哈尔滨工程大学 | 一种基于节点信息量统计的多波束地形分块方法 |
CN107038308B (zh) * | 2017-04-18 | 2019-09-13 | 南京工程学院 | 一种基于线性内插的规则格网地形建模方法 |
CN107038308A (zh) * | 2017-04-18 | 2017-08-11 | 南京工程学院 | 一种基于线性内插的规则格网地形建模方法 |
CN107798339B (zh) * | 2017-09-29 | 2021-04-27 | 天津大学 | 一种基于泰森多边形图的多波束异常值数据处理算法 |
CN107798339A (zh) * | 2017-09-29 | 2018-03-13 | 天津大学 | 一种基于泰森多边形图的多波束异常值数据处理算法 |
CN108681561A (zh) * | 2018-04-20 | 2018-10-19 | 贾帅东 | 一种航海图等深线平滑性的评估方法 |
CN108874932A (zh) * | 2018-05-31 | 2018-11-23 | 哈尔滨工程大学 | 一种基于改进的光线投射算法的海洋水声场三维可视化方法 |
CN108874932B (zh) * | 2018-05-31 | 2022-03-25 | 哈尔滨工程大学 | 一种基于改进的光线投射算法的海洋水声场三维可视化方法 |
CN109524990A (zh) * | 2018-12-03 | 2019-03-26 | 三峡大学 | 一种基于Voronoi图重心内插法的虚拟惯量配置方法 |
CN109524990B (zh) * | 2018-12-03 | 2022-06-17 | 三峡大学 | 一种基于Voronoi图重心内插法的虚拟惯量配置方法 |
CN110118526A (zh) * | 2019-03-08 | 2019-08-13 | 浙江中海达空间信息技术有限公司 | 一种支持实时监测的船载砂石体积自动计算方法 |
CN110927792A (zh) * | 2019-10-24 | 2020-03-27 | 广州海洋地质调查局 | 一种地震构造图的建立方法及处理终端 |
CN113806951A (zh) * | 2021-09-23 | 2021-12-17 | 华侨大学 | 一种基于半边数据结构的自然邻近点搜索的弹性仿真方法 |
CN113806951B (zh) * | 2021-09-23 | 2023-05-26 | 华侨大学 | 一种基于半边数据结构的自然邻近点搜索的弹性仿真方法 |
CN116127792A (zh) * | 2023-04-17 | 2023-05-16 | 北京世冠金洋科技发展有限公司 | 一种散乱数据的插值方法及装置 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102831645A (zh) | 一种应用于海底地形的数字高程模型的建立方法 | |
Zou et al. | Iso-parametric tool-path planning for point clouds | |
Han et al. | Point cloud simplification with preserved edge based on normal vector | |
CN102194253B (zh) | 一种面向三维地质层面结构的四面体网格生成方法 | |
CN111581776B (zh) | 一种基于几何重建模型的等几何分析方法 | |
Marchandise et al. | Cardiovascular and lung mesh generation based on centerlines | |
CN102867332A (zh) | 基于复杂边界约束的多级细分网格曲面拟合方法 | |
CN103927783A (zh) | 一种三维三角网构建填挖空间的图割方法 | |
CN101866015A (zh) | 石油地震构造图的建立装置 | |
Zhang et al. | Adaptive generation of hexahedral element mesh using an improved grid-based method | |
CN105869210A (zh) | 三维地质表面模型中的插值数据处理方法 | |
CN105931297A (zh) | 三维地质表面模型中的数据处理方法 | |
Held | VRONI and ArcVRONI: Software for and applications of Voronoi diagrams in science and engineering | |
CN105487116A (zh) | 层面模型建立方法 | |
Xiao-Ping et al. | An algorithm for generation of DEMs from contour lines considering geomorphic features | |
Lee | Automatic metric 3D surface mesh generation using subdivision surface geometrical model. Part 1: Construction of underlying geometrical model | |
Wang et al. | Point cloud hole filling based on feature lines extraction | |
Lee | Automatic metric 3D surface mesh generation using subdivision surface geometrical model. Part 2: Mesh generation algorithm and examples | |
CN114170388A (zh) | 一种基于八叉树的局部异向性搜索椭球体动态建模方法 | |
Dakowicz et al. | A unified spatial model for GIS | |
Liu et al. | 3D seabed terrain establishment based on moving fractal interpolation | |
CN105893492A (zh) | 三维地质表面模型中的曲面求交数据处理方法 | |
CN105894580A (zh) | 三维地质表面模型中的曲面延伸数据处理方法 | |
Gold | The dual is the context: Spatial structures for GIS | |
Zhu et al. | Computing the minimum distance between a point and a curve on mesh |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C02 | Deemed withdrawal of patent application after publication (patent law 2001) | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20121219 |