CN106446910B - 一种复杂地质曲面特征提取与重构方法 - Google Patents

一种复杂地质曲面特征提取与重构方法 Download PDF

Info

Publication number
CN106446910B
CN106446910B CN201610814861.2A CN201610814861A CN106446910B CN 106446910 B CN106446910 B CN 106446910B CN 201610814861 A CN201610814861 A CN 201610814861A CN 106446910 B CN106446910 B CN 106446910B
Authority
CN
China
Prior art keywords
edge
entering
point
catchment
terrain
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.)
Active
Application number
CN201610814861.2A
Other languages
English (en)
Other versions
CN106446910A (zh
Inventor
姚兴苗
王运璇
胡光岷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201610814861.2A priority Critical patent/CN106446910B/zh
Publication of CN106446910A publication Critical patent/CN106446910A/zh
Application granted granted Critical
Publication of CN106446910B publication Critical patent/CN106446910B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/30—Polynomial surface description
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2215/00—Indexing scheme for image rendering
    • G06T2215/06—Curved planar reformation of 3D line structures

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明公开了一种复杂地质曲面特征提取与重构方法，通过借鉴水文分析中河网追踪方面的相关知识并引入面追踪，解决了特征线提取时的碎片化问题，同时运用哈希表极大地减少了地形特征提取部分的计算量，极大地提高了提取效率。此外，本发明针对现有地形重构方法计算量大、地形骨架点无法确定会被保留等缺陷，采用了基于四阶偏微分方程的地形曲面平滑重构思想，在平滑过程中为保留原始曲面的地形特征而引入了特征点约束，在曲面平滑和曲面特征保留之间取得了较好的平衡。

Description

一种复杂地质曲面特征提取与重构方法
技术领域
本发明属于地形特征提取与重构技术领域，具体涉及一种复杂地质曲面特征提取与重构方法的设计。
背景技术
数字高程模型(DEM，Digital Elevation Model)的概念是在1958年由美国麻省理工学院的Miller教授提出来的，主要是针对于物理地形的一种数学建模，以便能运用计算机对现实世界中真实存在的地形进行分析和数字化处理。依据当前计算机的强大处理能力，相比于传统纸质地形模型(如地图)及图像，基于数字高程模型的数字地形至少具有以下优点：
(1)便于地形数据的存储与管理；
(2)便于对地形进行分析和计算；
(3)可以无损的缩放，便于分析地形细节；
(4)便于设计算法对地形进行智能化处理；
(5)便于实时更新。
从狭义角度而言，数字高程模型定义为区域地形表面海拔高度的数字化表达；从广义角度，数字高程模型是地理空间中地理对象表面海拔高度的数字化表达。数学上，数字高程模型为定义在二维平面上的函数，其参数为地形点在XOY平面上的点坐标，其值为该地形点的高程值。数字高程模型根据网格的不同分为基于格网的DEM和基于三角网的DEM。
在数字地形分析领域，基于数字高程模型的地形分析已经成为最主要的方向之一，其主要内容包括DEM的地形特征提取、地形重建等，同时这两者也是数字地形分析领域的重点与难点，相关文章已大量存在。数字地形分析在经济建设、社会发展、精细农业、智能化交通、现代化战争及人民生活等方面发挥着重要作用，这不仅给现代科学技术的发展带来了新的契机，在深层次上实现计算机科学、信息科学、地球科学、系统科学的有机融合，而且将对人们的生存方式乃至国家的发展战略产生深远的影响。
根据所选取的数字高程模型的类型不同，DEM地形特征的提取分为：基于规则格网DEM的地形特征提取、基于等高线DEM的地形特征提取和基于三角网DEM的地形特征提取。
(1)规则格网数字高程模型(格网DEM)。
规则格网DEM是地形分析中最典型的数据模型之一，许多国家的DEM数据都是以规则格网的数据矩阵形式提供的(如美国USGS提供分辨率为30m、90m的DEM；澳大利亚、日本、英国均提供50m分辨率的DEM；而加拿大为93m的DEM，法国为100m的DEM)。规则格网DEM以格网方式放置原始点，每个子单元为一个标准矩形，各地形点根据其X，Y坐标值与格网坐标系的映射关系置于各子格网顶点中，地形点在格网中的Z值为该点的高程值。规则格网DEM因其数据结构的规范性与简便性，使其成为了目前数字地形研究中主要采取的数据源格式。对于规则格网DEM中的地形点，与其拓扑相关的其他地形点主要为其上、下、左、右、左下、左上、右下、右上共八个领域点，点结构拓扑结构相对三角网DEM而言更容易构建。针对规则格网DEM的研究已经较为成熟。
(2)三角网(TIN)数字高程模型(三角网DEM)。
三角网数字高程模型也称为三角网DEM，是另一种在地形研究中被普通采用的模型结构，主要以三角网形式组织原始地形数据，即将原始地形点根据其X，Y坐标放置于平面三角网中某个三角形顶点处，高程对应于该点的Z值。三角网DEM通过互相邻接的三角形集来表示原始地形，组成三角网的三角形之间互不相交、互不包含，相比格网DEM对原始地形点分布位置的严格要求，三角网DEM在表示无规律、分布不均匀的原始地形点方面具有明显的优势。但鉴于TIN在数据结构、构建TIN、算法设计、管理等方面较为复杂，使得针对于三角网DEM的研究相对较少，继续研究的空间较大。
值得一提的是，规则格网DEM可非常简便地转为三角网DEM，即只需将规则格网的每个矩形中某条对角线连接起来即可。当然，以上三角化方法有较多缺陷，如未考虑连接的最佳性问题等，其他连接算法如DELAUNAY(德劳内)三角化算法能较好地解决三角化的最优问题。
根据提取算法的不同，地形特征的提取算法又分为：断面极值法、邻域比较法和流水模拟法。
(1)断面极值法提取地形特征点。
地形几何分析中提取地形特征点的一个典型算法为断面极值法，该方法认为地形断面曲线上的极大值点为地形脊点，而极小值点为地形谷点。算法的实现过程可总结为：先找出DEM的系列纵横断面上的极大值、极小值点，作为地形特征上的候选点，然后通过设计合适的特征辨别算法，滤除不符合条件的特征点以及将特征点归类到对应的特征线中。
算法缺点：不能顾及区域地形变化的特殊规律，对地形噪声的容错能力不高；易导致特征点(或线)提取的碎片化；可供求导的方向有限，易导致地形点属性的判断失误；无法对边界点进行判断。
(2)邻域比较法提取地形特征点。
邻域比较法的思想为对目标地形点与局部窗口的其他各点进行高程比较，比较方法和窗口大小的选取在各文献中稍有不同，其核心为：如果该点为窗口内的最高点，则该点为候选山脊点，反之则为候选山谷点；其特征属性的最终确认还需配合其他方法。
算法缺点：对特殊地形地貌的提取效果较差，对噪声的容错能力较差；生成的特征线相当零碎，在细微地方存在特征线交叉问题；阈值的大小设定对提取结果影响较大，需人工反复测试阈值的合适大小。
(3)流水模拟法提取地形特征点.
流水模拟法是三种地形特征提取算法中最复杂的一个，同时在大多数情况下也是提取效果最好的一个，算法经过简单变形就能应用于规则格网DEM和三角网DEM，因此算法的适用范围较广。流水模拟法也是目前研究相对不足的算法，可供继续研究的空间较大。
流水模拟法根据水流总是由高处流向低处这一物理规则来提取地形的山脊线与山谷线(两者也是最难准确提取的地形特征)，流水模拟法认为山脊线也即分水线，山谷线也即汇水线，即水流流入山脊点或山脊线时将向邻域出流，反之，水流将从邻域流向该区域内的山谷点或山谷线。考虑到山脊线与山谷线的定义和物理特性，该判断规则与实际情况基本相符。算法认为每点自身的水流量为1，经过水流的流入与流出，如果某点最终的汇水量超过某一设定的阈值，则该点为山谷点，山脊点的判断与提取和山谷点类似。
地形重构技术也称为地形生成技术，指的是根据已有地形点(或线)通过适当的重建算法生成原地形的三维模拟地形的技术。地形重构技术在地质勘探、军事活动、农业生产及水利方程等方面发挥着重要作用，同时也是逆向工程领域研究的重点内容之一。
地形重构可分为直接重构与间接重构，直接重构即通过将已有地形点直接映射到对应的空间位置，并依据点间拓扑关系将空间点连接成地形表面来重构原始地形。因为直接重构技术对地形点的严格要求，使得原始已有地形点极少能符合条件，这导致直接重构技术无法单独用于地形重构。
近年来，偏微分方程(Partial Differential Equation，PDE)方法作为一种高效的曲面设计工具被引进了计算机辅助设计领域。从理论上分析，偏微分方程曲面造型技术具有一大优势，即定义一张空间曲面的大部分信息来源于曲面的边界条件，这就可以仅通过少量的参数来重构和控制曲面。通过改变边界条件和偏微分方程中的形状控制参数，可以生成各种曲面形状。
哈希关系的本质是一种映射关系，将索引和索引对应的值以哈希规则连接起来，之后的查询只要知道了索引(称为哈希函数的键)值便可通过已知的哈希映射规则找到对应的值(称为哈希值)。这里的映射规则即为哈希函数(Hash Function)，索引值通过哈希函数得到的映射值称为哈希地址(Hash Address)，所有索引值与对应的哈希函数值组成一张对应表，称为哈希表(Hash Table)。哈希表冲突是指在给定的相同哈希规则下，不同的键值对应相同的哈希值。哈希冲突是设计哈希函数时必须考虑的重点和难点，理论上而言，哈希冲突不可能完全避免，所以设计的方向应该是考虑如何让哈希冲突发生的机率尽可能地小。
发明内容
本发明的目的是为了解决现有技术中基于三角网DEM的地形特征提取算法计算量大、效率低，特征点或特征线的提取易出现碎片化，且地形重构方法具有不平滑性的问题，提出了一种复杂地质曲面特征提取与重构方法。
本发明的技术方案为：一种复杂地质曲面特征提取与重构方法，包括以下步骤：
S1、提取满足唯一性约束的三角网候选山谷线；
S2、对候选山谷线进行追踪与结构建立；
S3、通过对原地形进行翻转后提取翻转地形的山谷线提取三角网山脊线；
S4、采用四阶偏微分方程进行曲面平滑重构。
进一步地，步骤S1包括以下分步骤：
S1-01、遍历边哈希表lineMap中的每条边；
S1-02、判断边是否遍历结束，若是则进入步骤S1-11，否则进入步骤S1-03；
S1-03、获取当前边及其左右三角形ID；
S1-04、判断左右三角形ID是否均不为空，若是则进入步骤S1-06，否则进入步骤S1-05；
S1-05、判定当前边为边界边，不考虑汇水性，返回步骤S1-01；
S1-06、根据左右三角形ID从三角网哈希表tinMap中获取当前边的左右邻接三角形；
S1-07、计算当前边左右邻接三角形的水流方向；
S1-08、计算当前边左右邻接三角形水流方向与该边的流向关系值；
S1-09、判断当前边的汇水特性并更新边哈希表lineMap；
S1-10、将汇水边加入候选汇水边集合，返回步骤S1-01；
S1-11、遍历三角网哈希表tinMap中的每个三角形；
S1-12、判断三角形是否遍历结束，若是则进入步骤S2，否则进入步骤S1-13；
S1-13、获取当前三角形的汇水边条数；
S1-14、判断汇水边是否多于一条，若是则进入步骤S1-15，否则返回步骤S1-11；
S1-15、对每条汇水边，计算其左右邻接三角形流向该边的汇水量之和；
S1-16、将原汇水边中汇水量最大的边汇水属性保持不变，其他边汇水属性改为未判断；
S1-17、更新边哈希表lineMap及候选汇水边集合，返回步骤S1-11。
进一步地，步骤S2包括以下分步骤：
S2-01、遍历边哈希表lineMap，获取边界点集；
S2-02、更新以给定点为起点的边哈希表linesWithStartpointMap及以给定点为终点的边哈希表linesWithEndpointMap；
S2-03、遍历点哈希表pointMap获取洼地点集；
S2-04、遍历候选汇水边，将终止于边界点或洼地点的边加入待追踪汇水边；
S2-05、遍历待追踪汇水边；
S2-06、判断待追踪汇水边是否遍历结束，若是则进入步骤S3，否则进入步骤S2-07；
S2-07、设当前遍历边为curLine；
S2-08、新建结果汇水边集；
S2-09、判断curLine是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-10；
S2-10、将curLine加入结果汇水边集；
S2-11、判断curLine是否为边界边，若是则进入步骤S2-24，否则进入步骤S2-12；
S2-12、从linesWithEndpointMap中获取流入curLine起点的汇水边集；
S2-13、判断获取的汇水边集是否为空，若是则进入步骤S2-15，否则进入步骤S2-14；
S2-14、遍历每一条获取到的汇水边，置curLine为当前遍历汇水边，返回步骤S2-09；
S2-15、边汇流追踪中断，考虑面汇流；
S2-16、判断连续追踪汇流面个数是否达到设定值，若是则进入步骤S2-24，否则进入步骤S2-17；
S2-17、判断出流对象是边还是面，若是边则进入步骤S2-18，若是面则进入步骤S2-19；
S2-18、获取流入curLine的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-19、获取流入三角形三边的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-20、判断curTIN是否存在，若是则进入步骤S2-21，否则进入步骤S2-24；
S2-21、判断curTIN是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-22；
S2-22、将curTIN的三边加入结果汇水边集；
S2-23、获取流入curTIN三个顶点的汇水边集，返回步骤S2-13；
S2-24、判定当前边的所有流入分支追踪完毕，返回步骤S2-05。
进一步地，步骤S4中若通过三角网数据提取出的山谷点或山脊点不在任何一个格网点上，则以该山谷点或山脊点为中心生成单个格网区域，该格网区域内的格网点的高程值被赋值为该特征点的高程值，并且不参与平滑，即该山谷点或山脊点单格网影响区域内的格网点均被认为是该山谷点或山脊点本身。
本发明的有益效果是：
(1)通过借鉴水文分析中河网追踪方面的相关知识并引入面追踪，解决了特征线提取时的碎片化问题。
(2)运用哈希表极大地减少了地形特征提取部分的计算量，极大地提高了提取效率。
(3)基于三角网数据直接重构地形的不平滑性，采用了基于四阶偏微分方程的地形曲面平滑重构思想，在平滑过程中为保留原始曲面的地形特征而引入了特征点约束，在曲面平滑和曲面特征保留之间取得了较好的平衡。
附图说明
图1为本发明提供的一种复杂地质曲面特征提取与重构方法流程图。
图2为本发明实施例的三角形特征边唯一性约束示意图。
图3为本发明步骤S1的分步骤流程图。
图4为本发明实施例的特征线面追踪示意图。
图5为本发明步骤S2的分步骤流程图。
图6为本发明实施例的格网点添加特征约束示意图。
图7为本发明实施例的无特征约束的四阶偏微分曲面重构效果图。
图8为本发明实施例的特征约束下的四阶偏微分曲面平滑重构效果图。
具体实施方式
下面结合附图对本发明的实施例作进一步的说明。
本发明实施例中用到的主要数据结构如表1所示：
表1
Figure BDA0001112731070000061
在表1中，点数据结构的点ID为字符串类型，由1自增长；三角网常规链表结构中的三角形ID为整型，由1自增长，三角形的三点ID为对应点结构的ID；三角网边结构中的边ID为字符串类型，由边两端点对应的点ID拼接而成，其中高程值较大的点ID在前，高程值较小的点ID在后，边起点和终点ID为对应点结构ID，边左、右邻接三角形属性值为相应的三角形ID值(约定边方向为由高程值较大的端点指向高程值较小的端点)，边汇水性为枚举类型，枚举值中的枚举元素包括汇水边、分水边、过水边、未判断。
本发明实施例中用到的主要哈希函数(均为全局变量)如下表2所示：
表2
Figure BDA0001112731070000062
Figure BDA0001112731070000071
本发明对哈希函数的运用广泛存在于三角网边结构的构建、特征线的判断与追踪、特征线结构的构建及特征线编码等各个方面，哈希思想的运用极大地减少了本发明地形特征提取部分的计算量，极大地提高了提取效率。
本发明提供的一种复杂地质曲面特征提取与重构方法，如图1所示，包括以下步骤：
S1、提取满足唯一性约束的三角网候选山谷线。
通过对每一条边进行遍历，判断其汇水特性，提取所有的汇水段作为候选山谷线，分水段作为候选山脊线。然而，通过上述方法有可能出现一个三角形上出现两条甚至三条山谷线段(或山脊线段)的情况，依据水文分析理论，三角网中任一三角形最多只能有一条山谷线段(或山脊线段)。
如图2所示，假设求得三角形AB与BC边均为汇水边，此时可通过比较边AB与边BC的过水面积，过水面积较大的边保持边汇水属性判断，而过水面积较小的边其属性改为未判断。比较AB边与BC边的过水面积等同于比较图2中边AE与CF的大小，判定方法如下：以点B为坐标原点，此时矢量S的方程为ax+by＝0，点A与点C的坐标分别变换为(xA-xB,yA-yB)和(xC-xB,yC-yB)，根据点到直线的距离公式
Figure BDA0001112731070000072
可得，当a(xA-xB)+b(yA-yB)+c|＞|a(xC-xB)+b(yC-yB)+c|时，边AB的过水面积大于边BC的过水面积，反之亦然。考虑到本发明采用的是三角网的边拓扑结构，则边的汇水量计算调整为该边左右两个邻接三角形的水流流向该边的汇水量之和。
因此，如图3所示，步骤S1具体包括以下分步骤：
S1-01、遍历边哈希表lineMap中的每条边；
S1-02、判断边是否遍历结束，若是则进入步骤S1-11，否则进入步骤S1-03；
S1-03、获取当前边及其左右三角形ID；
S1-04、判断左右三角形ID是否均不为空，若是则进入步骤S1-06，否则进入步骤S1-05；
S1-05、判定当前边为边界边，不考虑汇水性，返回步骤S1-01；
S1-06、根据左右三角形ID从三角网哈希表tinMap中获取当前边的左右邻接三角形；
S1-07、计算当前边左右邻接三角形的水流方向；
S1-08、计算当前边左右邻接三角形水流方向与该边的流向关系值；
S1-09、判断当前边的汇水特性并更新边哈希表lineMap；
S1-10、将汇水边加入候选汇水边集合，返回步骤S1-01；
S1-11、遍历三角网哈希表tinMap中的每个三角形；
S1-12、判断三角形是否遍历结束，若是则进入步骤S2，否则进入步骤S1-13；
S1-13、获取当前三角形的汇水边条数；
S1-14、判断汇水边是否多于一条，若是则进入步骤S1-15，否则返回步骤S1-11；
S1-15、对每条汇水边，计算其左右邻接三角形流向该边的汇水量之和；
S1-16、将原汇水边中汇水量最大的边汇水属性保持不变，其他边汇水属性改为未判断；
S1-17、更新边哈希表lineMap及候选汇水边集合，返回步骤S1-11。
S2、对候选山谷线进行追踪与结构建立。
候选山谷线集与地形的实际山谷线集并非两个完全相等的集合。地形山谷线提取算法也有零碎性的问题和洼地的处理问题，这也是大部分基于流水法的地形特征提取算法中的普遍性问题。综合上述情况考虑，需要在已提取出的候选山谷线中追踪出地形的实际山谷线，这里的追踪包含以下三层含义：(1)从候选山谷线中甄选出符合预期的实际山谷线；(2)连接甄选出来的山谷线，即解决山谷线的零碎性问题；(3)洼地情况的处理。
本发明实施例中，山谷线提取算法是基于山谷线段必须是三角网中某个三角形边的假设求得的，即河网分析中的边汇流处理，这也是导致山谷线(即汇水性)零碎化的主要原因。因此，在山谷线追踪阶段，将面汇流考虑进来，即考虑某个三角形的三边均为汇水边，水流通过该三角形面流出而非通过该三角形的某条边流出。如图4所示汇流边AO，但AO既不是边界点，也不是洼地底点，实际上，AO将通过与O点相连的三角形面与下游河道汇合，倘若不考虑三角形的面汇流，则局部流域将无法正常出流，与事实不符。
处理办法是选择与点O相连的三角形中向下坡度最大的三角形面作为AO的出流面，即将该三角形的三边作为山谷线与AO相连。这里的三角形面要符全以下三点要求：1、该三角形其余两点高程均应低于O点的高程以保证水流是从O点流出；2、该三角形是符合高程相求的三角形中坡度最大的；3、该三角形不能包含AO边(否则违背了山谷线唯一性约束)。但是，坡度汇流只是解决山谷线零碎性的辅助手段，因此，单条山谷线中考虑的连续坡面应有个限度，即如果三角形面三点均无汇水出流边与其相连时应继续考虑面汇流的三角形个数的限制。
上述考虑面汇流的关键点是最大坡度的三角形获取，三角形的坡度可用其相对于XOY平面的倾斜角来代替，而根据正弦函数在
Figure BDA0001112731070000091
角度范围的单调性，其倾斜角可用该角度的正弦值代替，因此三角形面的坡度可用其相对于XOY平面倾斜角的正弦值代替。该倾斜角的正弦值可通过以下方式求得：假设空间直角坐标系统中的某个三角形的三个顶点坐标分别为：A(x1,y1,z1)，B(x2,y2,z2)，C(x3,y3,z3)，则该三角形面的法向量计算如公式(1)所示：
Figure BDA0001112731070000092
其中M＝(y1-y2)(z1-z3)-(y1-y3)(z1-z2)，N＝(x1-x3)(z1-z2)-(x1-x2)(z1-z3)，L＝(x1-x2)(y1-y3)-(x1-x3)(y1-y2)。该法向量与天顶方向的夹角就是ABC三角面相对于水平面的倾斜角度，设该角度为θ，则由正弦定理及空间向量知识可知该角的正弦值可由下式(2)得到：
Figure BDA0001112731070000093
本发明中山谷线追踪算法是同时考虑了边汇流与面汇流而得到的，同时以假设山谷线(汇流线)终止于边界点或洼地点为追踪前提，这是甄别候选山谷线是否是实际地形的山谷线的重要依据。而边界点集和洼地点集可通过以下方法获取：遍历每一条边，如果当前遍历边的左邻接三角形或右邻接三角的为空，则该边为边界边，将该边的两个端点加入边界点集；遍历点，对于当前遍历点，从以给定点为起点的边中获取以当前遍历点为起点的边集，以及从以给定点为终点的边中获取以当前遍历点为终点的边集，如以当前遍历点为起点的边集的大小为0而以当前遍历点为终点的边集的大小大于0，则当前遍历点为洼地点，将其加入洼地点集。
综上所述，如图5所示，步骤S2包括以下分步骤：
S2-01、遍历边哈希表lineMap，获取边界点集；
S2-02、更新以给定点为起点的边哈希表linesWithStartpointMap及以给定点为终点的边哈希表linesWithEndpointMap；
S2-03、遍历点哈希表pointMap获取洼地点集；
S2-04、遍历候选汇水边，将终止于边界点或洼地点的边加入待追踪汇水边；
S2-05、遍历待追踪汇水边；
S2-06、判断待追踪汇水边是否遍历结束，若是则进入步骤S3，否则进入步骤S2-07；
S2-07、设当前遍历边为curLine；
S2-08、新建结果汇水边集；
S2-09、判断curLine是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-10；
S2-10、将curLine加入结果汇水边集；
S2-11、判断curLine是否为边界边，若是则进入步骤S2-24，否则进入步骤S2-12；
S2-12、从linesWithEndpointMap中获取流入curLine起点的汇水边集；
S2-13、判断获取的汇水边集是否为空，若是则进入步骤S2-15，否则进入步骤S2-14；
S2-14、遍历每一条获取到的汇水边，置curLine为当前遍历汇水边，返回步骤S2-09；
S2-15、边汇流追踪中断，考虑面汇流；
S2-16、判断连续追踪汇流面个数是否达到设定值，若是则进入步骤S2-24，否则进入步骤S2-17；
S2-17、判断出流对象是边还是面，若是边则进入步骤S2-18，若是面则进入步骤S2-19；
S2-18、获取流入curLine的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-19、获取流入三角形三边的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-20、判断curTIN是否存在，若是则进入步骤S2-21，否则进入步骤S2-24；
S2-21、判断curTIN是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-22；
S2-22、将curTIN的三边加入结果汇水边集；
S2-23、获取流入curTIN三个顶点的汇水边集，返回步骤S2-13；
S2-24、判定当前边的所有流入分支追踪完毕，返回步骤S2-05。
S3、通过对原地形进行翻转后提取翻转地形的山谷线提取三角网山脊线。
山脊线的提取与处理与山谷线基本一致，其依据是：如果将地形翻转，即地形最高点翻转为地形最低点，而地形最低点翻转为地形最高点，则山脊线便成了山谷线。这可通过将每个地形点高程值减去目标地形区的最大高程值实现，因此，山脊线其实可通过对原地形进行翻转后提取翻转地形的山谷线而得到。
S4、采用四阶偏微分方程进行曲面平滑重构。
本发明实施例中，采用的是四阶偏微分方程进行曲面平滑重构，主要原因是四阶偏微分方程能很好地解决重构曲面中因地形噪声而产生的毛刺现象，保证了重构曲面的光滑性。为了在曲面平滑重构时保留地形特征，如果目标点为地形特征点(包括山谷点、山脊点、山顶点等)，则其高程值保持不变，不对其进行平滑调整，而非特征点的处理不受影响。
通过三角网数据提取出的特征点(在本发明中为山谷点与山脊点)可能不在任何一个格网点上，因此为保留该特征点的约束，可以以该特征点为中心生成单个格网区域，该格网区域内的格网点的高程值被赋值为该特征点的高程值并且不参与平滑，即该特征点单格网影响区域内的格网点均被认为是该特征点本身。如图6所示，格网点B在特征点A的影响区域内，则认为格网点B即为特征点A，点B的高程值设置为点A的高程值，且点B不参与四阶偏微分平滑。
直接利用四阶偏微分方程重构的曲面俯视效果图如图7所示，在地形约束下利用四阶偏微分方程重构的曲面俯视效果图如下图8所示，图8中深色点为高程值较大的点，浅色点为高程值较小的点。由图7及图8可知特征约束下的四阶偏微分方程平滑方法生成的重构地形曲面很好地同时解决了地形平滑及保留地形特征的问题。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种复杂地质曲面特征提取与重构方法，其特征在于，包括以下步骤：
S1、提取满足唯一性约束的三角网候选山谷线；
S2、对候选山谷线进行追踪与结构建立；所述步骤S2包括以下分步骤：
S2-01、遍历边哈希表lineMap，获取边界点集；
S2-02、更新以给定点为起点的边哈希表linesWithStartpointMap及以给定点为终点的边哈希表linesWithEndpointMap；
S2-03、遍历点哈希表pointMap获取洼地点集；
S2-04、遍历候选汇水边，将终止于边界点或洼地点的边加入待追踪汇水边；
S2-05、遍历待追踪汇水边；
S2-06、判断待追踪汇水边是否遍历结束，若是则进入步骤S3，否则进入步骤S2-07；
S2-07、设当前遍历边为curLine；
S2-08、新建结果汇水边集；
S2-09、判断curLine是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-10；
S2-10、将curLine加入结果汇水边集；
S2-11、判断curLine是否为边界边，若是则进入步骤S2-24，否则进入步骤S2-12；
S2-12、从linesWithEndpointMap中获取流入curLine起点的汇水边集；
S2-13、判断获取的汇水边集是否为空，若是则进入步骤S2-15，否则进入步骤S2-14；
S2-14、遍历每一条获取到的汇水边，置curLine为当前遍历汇水边，返回步骤S2-09；
S2-15、边汇流追踪中断，考虑面汇流；
S2-16、判断连续追踪汇流面个数是否达到设定值，若是则进入步骤S2-24，否则进入步骤S2-17；
S2-17、判断出流对象是边还是面，若是边则进入步骤S2-18，若是面则进入步骤S2-19；
S2-18、获取流入curLine的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-19、获取流入三角形三边的坡度最大的汇水面curTIN，进入步骤S2-20；
S2-20、判断curTIN是否存在，若是则进入步骤S2-21，否则进入步骤S2-24；
S2-21、判断curTIN是否已处理过，若是则进入步骤S2-24，否则进入步骤S2-22；
S2-22、将curTIN的三边加入结果汇水边集；
S2-23、获取流入curTIN三个顶点的汇水边集，返回步骤S2-13；
S2-24、判定当前边的所有流入分支追踪完毕，返回步骤S2-05；
S3、通过对原地形进行翻转后提取翻转地形的山谷线提取三角网山脊线；
S4、采用四阶偏微分方程进行曲面平滑重构。
2.根据权利要求1所述的复杂地质曲面特征提取与重构方法，其特征在于，所述步骤S1包括以下分步骤：
S1-01、遍历边哈希表lineMap中的每条边；
S1-02、判断边是否遍历结束，若是则进入步骤S1-11，否则进入步骤S1-03；
S1-03、获取当前边及其左右三角形ID；
S1-04、判断左右三角形ID是否均不为空，若是则进入步骤S1-06，否则进入步骤S1-05；
S1-05、判定当前边为边界边，不考虑汇水性，返回步骤S1-01；
S1-06、根据左右三角形ID从三角网哈希表tinMap中获取当前边的左右邻接三角形；
S1-07、计算当前边左右邻接三角形的水流方向；
S1-08、计算当前边左右邻接三角形水流方向与该边的流向关系值；
S1-09、判断当前边的汇水特性并更新边哈希表lineMap；
S1-10、将汇水边加入候选汇水边集合，返回步骤S1-01；
S1-11、遍历三角网哈希表tinMap中的每个三角形；
S1-12、判断三角形是否遍历结束，若是则进入步骤S2，否则进入步骤S1-13；
S1-13、获取当前三角形的汇水边条数；
S1-14、判断汇水边是否多于一条，若是则进入步骤S1-15，否则返回步骤S1-11；
S1-15、对每条汇水边，计算其左右邻接三角形流向该边的汇水量之和；
S1-16、将原汇水边中汇水量最大的边汇水属性保持不变，其他边汇水属性改为未判断；
S1-17、更新边哈希表lineMap及候选汇水边集合，返回步骤S1-11。
3.根据权利要求2所述的复杂地质曲面特征提取与重构方法，其特征在于，所述步骤S4中若通过三角网数据提取出的山谷点或山脊点不在任何一个格网点上，则以该山谷点或山脊点为中心生成单个格网区域，该格网区域内的格网点的高程值被赋值为该特征点的高程值，并且不参与平滑，即该山谷点或山脊点单格网影响区域内的格网点均被认为是该山谷点或山脊点本身。
CN201610814861.2A 2016-09-12 2016-09-12 一种复杂地质曲面特征提取与重构方法 Active CN106446910B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610814861.2A CN106446910B (zh) 2016-09-12 2016-09-12 一种复杂地质曲面特征提取与重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610814861.2A CN106446910B (zh) 2016-09-12 2016-09-12 一种复杂地质曲面特征提取与重构方法

Publications (2)

Publication Number Publication Date
CN106446910A CN106446910A (zh) 2017-02-22
CN106446910B true CN106446910B (zh) 2020-10-27

Family

ID=58168929

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610814861.2A Active CN106446910B (zh) 2016-09-12 2016-09-12 一种复杂地质曲面特征提取与重构方法

Country Status (1)

Country Link
CN (1) CN106446910B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107341493B (zh) * 2017-07-14 2020-03-06 电子科技大学中山学院 基于分支相似度的地形特征线提取方法、装置及电子设备
CN107464287B (zh) * 2017-08-14 2021-04-27 电子科技大学 基于多目标优化的曲面重构方法
CN109191566B (zh) * 2018-08-22 2022-11-15 攀枝花学院 基于tin的三维梯田建模方法
CN110046397A (zh) * 2019-03-22 2019-07-23 广东辰誉电力科技有限公司 一种地表水汇流量定量计算方法
CN111735430B (zh) * 2020-06-04 2021-11-26 长江水利委员会长江科学院 一种河道断面地形重构方法
CN112529811A (zh) * 2020-12-17 2021-03-19 中国地质大学（武汉） 一种保留地形表面结构特征的dem数据去噪方法
CN113421341B (zh) * 2021-06-25 2023-04-21 中国电建集团昆明勘测设计研究院有限公司 一种水库冲刷淤积计算方法
CN113360594B (zh) * 2021-07-05 2023-09-26 中煤航测遥感集团有限公司 基于数字高程模型的汇水区提取方法、装置、设备及介质
CN114926599A (zh) * 2022-05-27 2022-08-19 甘肃省水利水电勘测设计研究院有限责任公司 一种地形曲面轻量化方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103177258A (zh) * 2013-03-29 2013-06-26 河南理工大学 一种根据矢量等高线数据自动提取地性线的方法
CN103345589A (zh) * 2013-07-19 2013-10-09 吴立新 一种顾及约束特征的城区汇水单元划分方法
CN104715507A (zh) * 2015-04-09 2015-06-17 武汉大学 一种基于曲面片的三维地理实体自动构建方法
CN105303615A (zh) * 2015-11-06 2016-02-03 中国民航大学 一种图像二维拼接与三维表面重建的组合方法
CN105760588A (zh) * 2016-02-04 2016-07-13 国家海洋局第海洋研究所 一种基于二层规则网格的sph流体表面重建方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103177258A (zh) * 2013-03-29 2013-06-26 河南理工大学 一种根据矢量等高线数据自动提取地性线的方法
CN103345589A (zh) * 2013-07-19 2013-10-09 吴立新 一种顾及约束特征的城区汇水单元划分方法
CN104715507A (zh) * 2015-04-09 2015-06-17 武汉大学 一种基于曲面片的三维地理实体自动构建方法
CN105303615A (zh) * 2015-11-06 2016-02-03 中国民航大学 一种图像二维拼接与三维表面重建的组合方法
CN105760588A (zh) * 2016-02-04 2016-07-13 国家海洋局第海洋研究所 一种基于二层规则网格的sph流体表面重建方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Hash Function for Triangular Mesh Reconstruction;Václav Skala等;《Proceedings of Spie the International Society for Optical Engineering 》;20021231;第1-6页 *
Hash functions and triangular mesh reconstruction;Jan Hrádek等;《Computers & Geosciences》;20031231;第741–751页 *
基于 Hash 函数的 TIN 拓扑关系重建;潘胜玲等;《地理与地理信息科学》;20060331;第22卷(第2期);第21-24,29页 *
基于三⻆网DEM的地形特征提取算法;苏丹阳等;《应用基础与工程科学学报》;20091130;第38-42页 *
基于不规则三角网TIN的流域特征自动提取算法与原型系统设计研究;任政;《中国优秀硕士学位论文全文数据库（电子期刊）基础科学辑》;20090115;第三章 *
基于偏微分方程的复杂地质曲面光滑重构;邓世武;《中国优秀硕士学位论文全文数据库（电子期刊）基础科学辑》;20160415;第三章、第四章 *

Also Published As

Publication number Publication date
CN106446910A (zh) 2017-02-22

Similar Documents

Publication Publication Date Title
CN106446910B (zh) 一种复杂地质曲面特征提取与重构方法
CN103236086B (zh) 一种顾及地表水文上下文的多尺度dem建模方法
US20120166160A1 (en) Block model constructing method for complex geological structures
CN106097450B (zh) 一种等高线生成方法及装置
CN110276732B (zh) 一种顾及地形特征线要素的山区点云空洞修复方法
CN101838958B (zh) 道路坡度的检测方法
CN111858810B (zh) 一种面向道路dem构建的建模高程点筛选方法
CN102521884A (zh) 一种基于LiDAR数据与正射影像的3维屋顶重建方法
CN104574512A (zh) 一种顾及地形语义信息的多尺度dem构建方法
CN105160658A (zh) 一种基于子流域边界和流路特征的山脊线提取方法
CN110189409B (zh) 一种基于plaxis的快速真三维地质建模方法及系统
CN114926602B (zh) 基于三维点云的建筑物单体化方法及系统
CN109671149A (zh) 基于dem的地形素描图自动绘制方法
Tseng et al. Extraction of building boundary lines from airborne lidar point clouds
CN115861247A (zh) 一种高分辨率遥感影像轮廓多级正则化方法、系统及应用
KR20140024590A (ko) 연속지적도의 다축척 모델 생성 방법
CN110660027B (zh) 一种针对复杂地形的激光点云连续剖面地面滤波方法
CN105894553B (zh) 一种基于格栅选择的街巷空间形态布局方法
CN114119902A (zh) 一种基于无人机倾斜三维模型的建筑物提取方法
CN114463494B (zh) 一种地形特征线自动提取方法
KR101063827B1 (ko) 한국토지정보시스템 연속지적도와 수치지형도의 기하학적 지도 변환을 위한 반자동화된 공액점 쌍 추출방법
Xu et al. Methods for the construction of DEMs of artificial slopes considering morphological features and semantic information
CN115359149A (zh) 一种泛化拓扑结构约束下的连续图斑化简方法
Yao et al. Research on three-dimensional model reconstruction of slope erosion based on sequence images
CN111260758B (zh) 一种面状行政区域层次关系的构建方法及系统

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant