CN106649466B - 数字地图中典型地形几何参数获取方法 - Google Patents

数字地图中典型地形几何参数获取方法 Download PDF

Info

Publication number
CN106649466B
CN106649466B CN201610853996.XA CN201610853996A CN106649466B CN 106649466 B CN106649466 B CN 106649466B CN 201610853996 A CN201610853996 A CN 201610853996A CN 106649466 B CN106649466 B CN 106649466B
Authority
CN
China
Prior art keywords
landform
point
boundary
going
mountain
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
CN201610853996.XA
Other languages
English (en)
Other versions
CN106649466A (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.)
Xidian University
Original Assignee
Xidian University
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 Xidian University filed Critical Xidian University
Priority to CN201610853996.XA priority Critical patent/CN106649466B/zh
Publication of CN106649466A publication Critical patent/CN106649466A/zh
Application granted granted Critical
Publication of CN106649466B publication Critical patent/CN106649466B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases

Landscapes

  • Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Data Mining & Analysis (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明涉及一种数字地图中典型地形几何参数获取方法，包括地形分类及识别、地形简化、地形几何参数提取等过程。针对ASTER GDEM类型的GeoTIFF格式的数字地图，借助ArcGIS工具，通过计算归一化的高程、坡度、坡度变率、曲率、地形起伏度，然后以其为特性指标采用ISO(Iterative SelfOrganization)非监督聚类得到地形的初步分类，再按照各类地形平均起伏度的大小，将地形分为平原、丘陵、低山、中山、高山等几种地形；最终得到这几种典型地形的位置、底面及高度尺寸等几何参数。本发明具有可行、实用的优点。

Description

数字地图中典型地形几何参数获取方法
技术领域
本发明属于地理信息系统(GIS)和计算机图学领域，涉及地形的分类方法和地形参数提取及简化模型的建立方法，具体地说就是数字地图中典型地形几何参数获取方法，可用于野外大型地形场景的电波传播近似计算和快速实时显示。
背景技术
随着近年来计算机、卫星遥感等科技的迅猛发展，数字化浪潮正在席卷全球。地球作为人类活动和生活的首要载体，人们在认识自然过程中以及对自然进行改造过程中，对周围地形地貌的环境信息一直不断试着利用不同的方法进行描述及表达。在众多的方法中，数字地形图这种表达方法能够较好对地形表面形态状况进行描述和表达。通过地理测绘等相关手段得到数字化坐标数据，地貌仿真后产生三维可视化效果，地形地貌结构可以很好被表达，这种方法在二十世纪以来应用较为广泛。上个世纪四十年代计算机技术诞生及其不断蓬勃发展，以及伴随着计算机图学(Computer Graphics)、现代数学理论、数据挖掘、模式识别等相关技术理论完善和应用，都极大地促进了数字地图的发展及应用。除了在互联网和手机等行业数字地图有着广泛的应用之外，数字地图还有许多其他的重要应用前景，例如我们如果能把各种物体、个人、建筑、交通、景观等动态的和静态的数字信息和三维数字地图有机结合成一个整体，就会形成一个数字化的城市应用服务场景，将对我们的生活产生无可比拟的便捷和享受。对于数字地图，常常是通过读取各种不同精度的高程原始坐标数据构成多边形(常用三角形)面片，然后运用多边形面片无限逼近真实地貌。数字高程是利用有序的数据阵列来存储地面高程，其模型(DigitalElevation Model，简称DEM)是归属于数字地形模型(Digital TerrainModel，简称DTM)的一个分支。由于DEM数据是一个DTM分析研究所需的基本数据源，所以我们希望从中提取各种各样有用的信息，可为国民经济发展、生态农业、资源利用和作战指挥等相关工作提供参考。DEM就是用数字建模的形式来表达现实地球表面环境的方法，自从二十世纪五十年代被采用，在各种应用领域得到了极为广泛的应用。随着计算机软硬件水平以及计算机图学等技术爆发式的发展，DEM从获取方式、存储结构以及处理速度等方面都获得了质的飞跃。如今，伴随大量可靠高精度DEM数据的产生，亟需根据各种现实需求场景开展DEM数据应用研究。
如今的数字地图的发展促使我们通过各种途径获得越来越精确的大量地图信息，伴随着虚拟数字地球概念的提出，数字地图的用途不断从军用逐渐延伸到民用，日常的生活场景更是逐渐离不开其相关的应用。三维地图数据挖掘就是要以实际需求为出发点，从三维地图数据中挖掘出可行的、规律的、有价值的信息，使之服务于实际研究和现实生产。三维地图数据的挖掘可以通过采用诸多挖掘方法，例如通过统计、专家系统(经验体系)、机器学习、归纳法和模式识别等。地形地貌分类的传统获取方式主要利用人工手段辨识地形信息，其工作量大、周期长、难以完成大范围地形信息的快速低成本获取，针对遥感地图自动识别的问题，目前也尚无完全自动化的针对性的实用自动分类方法。现代社会生活中通信电子设备越来越多，不断面临着电子设备之间的耦合问题，也就是电磁干扰与兼容的问题。电磁兼容(ElectromagneticCompatibility，EMC)作为现代热门的重要前言研究，主要包括了两部分：EMI(电磁干扰)和EMS(电磁耐受性)。随着研究不断深入，各种电磁分析软件不断专业化、自动化，其中关于电波传播计算的场景建模主要有两种方式：一种是依靠人机交互的方式手工输入复杂地形环境数据实现建模；另一种是基于数字地图的自动建模，但其场景限于城市及其周围范围，场景中待计算的物及地形较为单一、尺寸不是很大。
数字地图的地形信息提取研究内容包括对山顶点的提取，对山底面信息的提取，对地形的分类，以及对所提取的地形进行简化。文献“地理信息系统原理及应用”(电子工业出版社，2011，作者：胡祥培等)从理论的角度，详细阐述了空间地理数据的获取以及处理，并论述了空间数据的管理和地理空间分析；硕士论文“基于DEM的地形特征提取算法研究及应用”(西安建筑科技大学，2012，作者：易炜)研究了目前现有的各类山顶点、鞍部点和山脊线的提取算法，对各种地形进行了定性和定量的分析，并通过对DEM格网数据进行等高距分层以模拟其在等高线地图中呈现的特性和拓扑结构，提出了相应的地形特征提取算法；硕士论文“利用图像处理技术提取地形特征线的方法研究”(西安建筑科技大学，2011，作者：方莉)在认真分析地形特征线的物理形态和结构特征的基础上，提出了一种利用灰度形态学算子提取山脊线和山谷线的新方法，另外针对山脚线的提取利用其处于山体和平地交界处这一典型物理特征，考虑将灰度图像的统计原理与地形坡度相结合设计一种将平地区域分割出来的方法，最后通过提取平地的边缘得到山脚线；论文“地貌认知及空间剖分的山顶点提取”(测绘科学，2010年9月，126-127，作者：罗明良等)通过对山顶局部形态特征及空间定量特征分析，结合传统地貌学认知思路，给出了基于空间剖分的山顶点快速提取方法，空间剖分标准对应于地貌学山地分类，通过高差限制剖分实现DEM分块，并在此基础上拟合函数，给出符合地貌认知的山顶点；论文“山顶点类型及其形态特征数字表达”(南京师范大学学报，2010年3月，136-130，作者：苍学智等)在分析了山顶点所具有的空间特征、成因特征、尺度特征及点群特征的基础上，依照科学性、系统性、实用性、可实现性的原则，对山顶点进行了系统的分类，并阐述了其定量描述方法，并且还深入探讨了山顶点群的基本类型、定量描述指标与地学意义。可以看出，已有研究成果存在以下问题：(1)研究仅涉及到地形特殊点(如山顶点、鞍部点)、特殊线(山脊线、山谷线、山脚线)等局部几何特征的识别，缺乏对地形整体几何形状识别、简化和几何参数提取的研究；(2)研究方法较为单一，主要基于几何图形学，未充分利用图像学方法进行研究；(3)未充分利用地形学中的坡度、剖面曲率、正切曲率、地形起伏度等地形因子开展研究，存在一定的局限性。
发明内容
本发明的目的在于克服上述现有研究中存在的问题，提供一种数字地图中典型地形几何参数获取方法，以满足野外大型地形场景的电波传播近似计算和快速实时建模及显示的需要。
本发明的目的是这样实现的：数字地图中典型地形几何参数获取方法，其特征是：至少包括如下步骤：
步骤101：使用ArcGIS工具，打开一个ASTER GDEM类型的GeoTIFF地图文件；
步骤102：将所有栅格点的高程值作归一化处理，结果保存在DEM_1.tif文件中；归一化处理的方法为：设x为任一栅格点处的量值，该量值是高程、坡度、坡度变率、曲率、地形起伏度中的任一种，xmax是所有栅格点量值的最大值，xmin是所有栅格点量值的最小值，x′是x的归一化值，则x′＝(x-xmin)×255/(xmax-xmin)；
步骤103：计算所有栅格点的坡度值，再作归一化处理，归一化的坡度值保存在Slope_1.tif文件中；
步骤104：计算所有栅格点的坡度变率，再作归一化处理，归一化的坡度变率值保存在PoDuBianLv_1.tif文件中；
步骤105：计算所有栅格点的曲率，并对曲率求绝对值，再作归一化处理，归一化的曲率值保存在Curvature_1.tif文件中；
步骤106：计算所有栅格点的地形起伏度，结果保存在QiFuDu.tif文件中，再对其作归一化处理，归一化的地形起伏度值保存在QiFuDu_1.tif文件中；
步骤107：将归一化的高程、坡度、坡度变率、曲率和地形起伏度做为特性指标，采用ISO(Iterative SelfOrganization)非监督聚类得到地形的初步分类结果，每类地形有一个类别编号，初步分类结果保存在JuLei.tif文件中；
步骤108：使用GDAL库函数打开步骤101的地图文件，打开
QiFuDu.tif文件和JuLei.tif文件；
步骤109：依据地形平均起伏度对步骤107初步分类结果进行地形判定分类；
步骤110：用四位正整数标志地形类别及区域，前一位表示地形类别，平原、丘陵、低山、中山、高山分别用1、2、3、4、5表示；后三位表示此种类别地形的区域个数编号，从001开始累积编号，地形判定分类的结果用此种方式表示成数值之后，保存在地形分类结果文件FenLei.tif之中；
步骤111：读取分类结果文件FenLei.tif，根据地形类别及区域标志数字选择待研究地形区域；
步骤112：判断地形类别标志是否为1，若是，转至步骤113；若不是，转至步骤115；
步骤113：将平原地形简化为平面，求取地形类别标志为平原的所有栅格点高程的平均值，记其为平均高程Ea；
步骤114：求取地形类别标志为平原的所有栅格点地形起伏度的平均值，记其为平均起伏度Ha，转至步骤127；
步骤115：对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中；
步骤116：求地形区域的面积S、地形区域的中心C(x0,y0)和地形的体积V：由链表Point_List存储的边界点序列连线形成一个边界多边形，用图形学中的“交点计数法”判断栅格点是否在边界多边形内，统计处在边界多边形内和位于边界多边形边上的栅格点数目num，设单个栅格点所占面积为a，则地形区域的面积S＝num×a；地形区域的中心点C的坐标xi、yi是边界多边形内及边界多边形边上的任一栅格点的坐标；地形的体积hi是边界多边形内及边界多边形边上的任一栅格点的高程；
步骤117：地形类别标志是否为2，若是，转至步骤118；若不是，转至步骤120；
步骤118：将丘陵地形简化为球缺，求取球缺的底面圆半径
步骤119：求取球缺的高度h1：解方程(h1)3+3(r1)2(h1)-6V/π＝0即得球缺的高度转至步骤127；
步骤120：求地形区域的周长L：由链表Point_List存储的边界点序列连线形成一个边界多边形，边界点即为边界多边形的顶点，设多边形顶点个数为n，令第n+1顶点等于第1个顶点，边界多边形任一边的边长|PiPi+1|由边的两个顶点Pi、Pi+1通过两点距离公式求得，则边界多边形的周长即为地形区域的周长L，
步骤121：求地形区域的形状因子f：形状因子f＝4πS/L2；
步骤122：判断形状因子f是否大于或等于0.7，若是，转至步骤123；若不是，转至步骤125；
步骤123：将地形简化为锥形山，求锥形山底面圆的半径
步骤124：求取锥形山的高度转至步骤127；
步骤125：将地形简化为楔形山，求取楔形山底面区域的最小外接矩形；
步骤126：求取楔形山的高度其中a、b为步骤125求取的外接矩形的两个边长，转至步骤127；
步骤127：输出地形几何参数：若地形为平原，输出平均高程Ea、平均起伏度Ha；若地形为丘陵，输出球缺底面圆中心C(x0,y0)、半径r1和球缺高度h1；若地形为锥形山，输出圆锥底面圆中心C(x0,y0)、半径r1和圆锥高度h2；若地形为楔形山，输出楔形体底面矩形中心C(x0,y0)，边长a、b，矩形方向角an和楔形体高度h3。
所述的步骤109中依据地形平均起伏度对步骤107初步分类结果进行地形判定分类，包括以下步骤：
步骤201：读取JuLei.tif文件，遍历地形类别编号，记该编号为n；
步骤202：遍历类别n的地形范围，从QiFuDu.tif文件中读取相对应位置的地形起伏度数据；
步骤203：求取该类别地形区域的地形起伏度平均值U；
步骤204：判断U是否小于或等于20米，若是，转至步骤205；若不是，转至步骤206；
步骤205：该类别地形判定为平原；
步骤206：判断U是否小于或等于200米，若是，转至步骤207；若不是，转至步骤208；
步骤207：该类别地形判定为丘陵；
步骤208：判断U是否小于或等于500米，若是，转至步骤209；若不是，转至步骤210；
步骤209：该类别地形判定为低山；
步骤210：判断U是否小于或等于1500米，若是，转至步骤211；若不是，转至步骤212；
步骤211：该类别地形判定为中山；
步骤212：该类别地形判定为高山。
所述的步骤115中对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中，包括以下步骤：
步骤301：基于GDAL读取地形分类结果文件FeiLei.tif；
步骤302：根据待研究的地形区域的四位正整数标志，将其地形类别标志数字记作type；
步骤303：按照从上至下、由左到右的顺序，从文件FeiLei.tif中图像左上角依次读取栅格点的地形类别标志数字，当第一次读到的栅格点的地形类别标志数字等于type时，则将该点存入边界点序列链表Point_List中，并记录该点为边界起点；
步骤304：令当前边界点等于边界起点，令下一个边界点的搜索方向码Directcode＝0；边界点搜索采用八邻域方法，每个栅格点有八个其他的栅格点与其相邻，该栅格点正上方栅格点的搜索方向码记为0，该栅格点右上方、正右方、右下方、正下方、左下方、正左方、左上方栅格点的搜索方向码依次记为1、2、3、4、5、6、7；
步骤305：读取当前边界点沿Directcode方向的搜索点的地形类别标志数字，将其记作value；
步骤306：判断value是否等于type，若是，转至步骤308；若不是，转至步骤307；
步骤307：令Directcode＝Directcode+1，转至步骤305；
步骤308：令当前边界点等于此搜索点；
步骤309：判断当前边界点是否等于边界起点，若是，则边界跟踪结束；若不是，转至步骤310；
步骤310：将当前边界点存入链表Point_List中；
步骤311：令Directcode＝Directcode-2；
步骤312：判断Directcode是否大于或等于0，若是，转至步骤305；若不是，转至步骤313；
步骤313：令Directcode＝Directcode+8，转至步骤305。
所述的步骤125中将地形简化为楔形山，求取楔形山底面区域的最小外接矩形，包括以下步骤：
步骤401：读取地形区域边界点序列链表Point_List，令链表Rot_List等于链表Point_List，ang＝3°,定义二维数组Data[30][4]；
步骤402：对链表Rot_List中的多边形各顶点做以地形区域中心C(x0,y0)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Rot_List中；
步骤403：求链表Rot_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；
步骤404：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；
步骤405：将步骤402～404循环29次；
步骤406：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m；
步骤407：读取数组Data中第m行的前3列值，分别用an、a、b记录；令an＝90°-an，定义an为矩形方向角，它是长度为a的矩形边与正东方向的夹角；
步骤408：输出地形底面区域的最小外接矩形边长a、b及矩形方向角an。
本发明有如下优点：
(1)提出了图像处理与图形处理相结合的典型地形整体几何形状识别的方法；
(2)通过地形简化，解决了基于数字地图的典型地形快速几何建模问题；
(3)所提方法可行、实用。
附图说明
图1是本发明的总流程图；
图2是地形判定分类的流程图；
图3是边界跟踪提取的流程图；
图4是楔形山底面区域最小外接矩形求取的流程图；
图5是一座山的实际形状；
图6是图5所示的山进行地形分类和边界跟踪提取后图示化的结果；
图7是图5所示山简化并提取几何参数后再显示的结果。
具体实施方式
本发明所读取的数据源，是一种数字高程模型(DEM)，数据格式是ASTER GDEM(Advanced Spaceborn Thermal Emission and Reflection Radiometer Global DigitalElevation Model)类型的GeoTIFF遥感影像数据(.tif文件)，该图像栅格点对应的空间分辨率约为30m×30m。借助ArcGIS工具，计算归一化的高程、坡度、坡度变率、曲率、地形起伏度，以其为特性指标采用ISO(Iterative SelfOrganization)非监督聚类得到地形的初步分类结果；通过GDAL(Geospatial Data Abstraction Library)这种API读取初步分类结果文件并处理，按照各类地形平均起伏度的大小，将地形分为平原、丘陵、低山、中山、高山等几种地形；之后，对每一类地形，通过八邻域边界跟踪和地形区域底面的面积、中心位置、周长、形状因子、最小外接矩形的计算以及区域的体积计算，将地形简化为丘陵、锥形山、楔形山等典型地形并得到底面的中心位置和几何尺寸，再根据体积和底面尺寸求出地形的高度尺寸，最终得到这几种典型地形的位置、底面及高度尺寸等几何参数。
参照图1，本发明的数字地图中典型地形几何参数获取方法包括如下步骤：
步骤101：使用ArcGIS工具，打开一个ASTER GDEM类型的GeoTIFF地图文件，如打开ASTGTM2_N32E115_dem.tif文件；
步骤102：将所有栅格点的高程值作归一化处理，结果保存在DEM_1.tif文件中；归一化处理的方法为：设x为任一栅格点处的量值，该量值是高程、坡度、坡度变率、曲率、地形起伏度中的任一种，xmax是所有栅格点量值的最大值，xmin是所有栅格点量值的最小值，x′是x的归一化值，则x′＝(x-xmin)×255/(xmax-xmin)；归一化处理的目的是将有量纲的数据转化为无量纲的数据并使转化后的数据值范围为0～255,为下面地形的模糊聚类分类做准备；
步骤103：计算所有栅格点的坡度值，再作归一化处理，归一化的坡度值保存在Slope_1.tif文件中；
步骤104：计算所有栅格点的坡度变率，再作归一化处理，归一化的坡度变率值保存在PoDuBianLv_1.tif文件中；
步骤105：计算所有栅格点的曲率，并对曲率求绝对值，再作归一化处理，归一化的曲率值保存在Curvature_1.tif文件中；
步骤106：计算所有栅格点的地形起伏度，结果保存在QiFuDu.tif文件中，再对其作归一化处理，归一化的地形起伏度值保存在QiFuDu_1.tif文件中；地形起伏度统计范围取60×60栅格点，即统计范围为1.8km×1.8km；
步骤107：将归一化的高程DEM_1.tif、坡度Slope_1.tif、坡度变率PoDuBianLv_1.tif、曲率Curvature_1.tif和地形起伏度QiFuDu_1.tif做为特性指标，设定类数目如25，采用ISO(Iterative SelfOrganization)非监督聚类得到地形的初步分类结果，共有17类，每类地形有一个类别编号，初步分类结果保存在JuLei.tif文件中；
步骤108：使用GDAL库函数打开步骤101的地图文件ASTGTM2_N32E115_dem.tif，打开QiFuDu.tif文件和JuLei.tif文件；
步骤109：依据地形平均起伏度对步骤107初步分类结果进行地形判定分类；
步骤109中依据地形平均起伏度对步骤107初步分类结果进行地形判定分类，参照图2，包括以下步骤：
步骤201：读取JuLei.tif文件，遍历地形类别编号，记该编号为n；
步骤202：遍历类别n的地形范围，从QiFuDu.tif文件中读取相对应位置的地形起伏度数据；
步骤203：求取该类别地形区域的地形起伏度平均值U；
步骤204：判断U是否小于或等于20米，若是，转至步骤205；若不是，转至步骤206；
步骤205：该类别地形判定为平原；
步骤206：判断U是否小于或等于200米，若是，转至步骤207；若不是，转至步骤208；
步骤207：该类别地形判定为丘陵；
步骤208：判断U是否小于或等于500米，若是，转至步骤209；若不是，转至步骤210；
步骤209：该类别地形判定为低山；
步骤210：判断U是否小于或等于1500米，若是，转至步骤211；若不是，转至步骤212；
步骤211：该类别地形判定为中山；
步骤212：该类别地形判定为高山。
步骤110：用四位正整数标志地形类别及区域，前一位表示地形类别，平原、丘陵、低山、中山、高山分别用1、2、3、4、5表示；后三位表示此种类别地形的区域个数编号，从001开始累积编号，地形判定分类的结果用此种方式表示成数值之后，保存在地形分类结果文件FenLei.tif之中；
步骤111：读取分类结果文件FenLei.tif，根据地形类别及区域标志数字选择待研究地形区域；
步骤112：判断地形类别标志是否为1，若是，转至步骤113；若不是，转至步骤115；
步骤113：将平原地形简化为平面，求取地形类别标志为平原的所有栅格点高程的平均值，记其为平均高程Ea；
步骤114：求取地形类别标志为平原的所有栅格点地形起伏度的平均值，记其为平均起伏度Ha，转至步骤127；
步骤115：对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中；
步骤115中对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中，参照图3，包括以下步骤：
步骤301：基于GDAL读取地形分类结果文件FeiLei.tif；
步骤302：根据待研究的地形区域的四位正整数标志，将其地形类别标志数字记作type；
步骤303：按照从上至下、由左到右的顺序，从文件FeiLei.tif中图像左上角依次读取栅格点的地形类别标志数字，当第一次读到的栅格点的地形类别标志数字等于type时，则将该点存入边界点序列链表Point_List中，并记录该点为边界起点；
步骤304：令当前边界点等于边界起点，令下一个边界点的搜索方向码Directcode＝0；边界点搜索采用八邻域方法，每个栅格点有八个其他的栅格点与其相邻，该栅格点正上方栅格点的搜索方向码记为0，该栅格点右上方、正右方、右下方、正下方、左下方、正左方、左上方栅格点的搜索方向码依次记为1、2、3、4、5、6、7；
步骤305：读取当前边界点沿Directcode方向的搜索点的地形类别标志数字，将其记作value；
步骤306：判断value是否等于type，若是，转至步骤308；若不是，转至步骤307；
步骤307：令Directcode＝Directcode+1，转至步骤305；
步骤308：令当前边界点等于此搜索点；
步骤309：判断当前边界点是否等于边界起点，若是，则边界跟踪结束；若不是，转至步骤310；
步骤310：将当前边界点存入链表Point_List中；
步骤311：令Directcode＝Directcode-2；
步骤312：判断Directcode是否大于或等于0，若是，转至步骤305；若不是，转至步骤313；
步骤313：令Directcode＝Directcode+8，转至步骤305。
步骤116：求地形区域的面积S、地形区域的中心C(x0,y0)和地形的体积V：由链表Point_List存储的边界点序列连线形成一个边界多边形，用图形学中的“交点计数法”判断栅格点是否在边界多边形内，统计处在边界多边形内和位于边界多边形边上的栅格点数目num，设单个栅格点所占面积为a，则地形区域的面积S＝num×a；地形区域的中心点C的坐标xi、yi是边界多边形内及边界多边形边上的任一栅格点的坐标；地形的体积hi是边界多边形内及边界多边形边上的任一栅格点的高程；
步骤117：地形类别标志是否为2，若是，转至步骤118；若不是，转至步骤120；
步骤118：将丘陵地形简化为球缺，求取球缺的底面圆半径
步骤119：求取球缺的高度h1：解方程(h1)3+3(r1)2(h1)-6V/π＝0即得球缺的高度转至步骤127；
步骤120：求地形区域的周长L：由链表Point_List存储的边界点序列连线形成一个边界多边形，边界点即为边界多边形的顶点，设多边形顶点个数为n，令第n+1顶点等于第1个顶点，边界多边形任一边的边长|PiPi+1|由边的两个顶点Pi、Pi+1通过两点距离公式求得，则边界多边形的周长即为地形区域的周长L，
步骤121：求地形区域的形状因子f：形状因子f＝4πS/L2；
步骤122：判断形状因子f是否大于或等于0.7，若是，转至步骤123；若不是，转至步骤125；
步骤123：将地形简化为锥形山，求锥形山底面圆的半径
步骤124：求取锥形山的高度转至步骤127；
步骤125：将地形简化为楔形山，求取楔形山底面区域的最小外接矩形；
步骤125中将地形简化为楔形山，求取楔形山底面区域的最小外接矩形，参照图4，包括以下步骤：
步骤401：读取地形区域边界点序列链表Point_List，令链表Rot_List等于链表Point_List，ang＝3°,定义二维数组Data[30][4]；
步骤402：对链表Rot_List中的多边形各顶点做以地形区域中心C(x0,y0)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Rot_List中；
步骤403：求链表Rot_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；
步骤404：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；
步骤405：将步骤402～404循环29次；
步骤406：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m；
步骤407：读取数组Data中第m行的前3列值，分别用an、a、b记录；令an＝90°-an，定义an为矩形方向角，它是长度为a的矩形边与正东方向的夹角；
步骤408：输出地形底面区域的最小外接矩形边长a、b及矩形方向角an。
步骤126：求取楔形山的高度其中a、b为步骤125求取的外接矩形的两个边长，转至步骤127；
步骤127：输出地形几何参数：若地形为平原，输出平均高程Ea、平均起伏度Ha；若地形为丘陵，输出球缺底面圆中心C(x0,y0)、半径r1和球缺高度h1；若地形为锥形山，输出圆锥底面圆中心C(x0,y0)、半径r1和圆锥高度h2；若地形为楔形山，输出楔形体底面矩形中心C(x0,y0)，边长a、b，矩形方向角an和楔形体高度h3。
仿真实例
利用本发明对国内某座名为孤峰山的地形进行几何参数提取。图5是读取处理孤峰山的DEM数据并以网格形式显示的孤峰山的实际形状。通过对孤峰山所在的ASTER GDEM类型的GeoTIFF格式的数字地图文件用ArcGIS工具计算归一化的高程、坡度、坡度变率、曲率、地形起伏度并以其为特性指标采用ISO非监督聚类得到地形初步分类；再对初步分类结果按照各类地形平均起伏度的大小进行地形判定分类，选取孤峰山作为待研究地形区域，图6是该地形区域分类和边界跟踪提取后图示化的结果，有平原和中山两种地形类别，品红色表示中山，绿色表示平原；通过对提取的边界求面积、中心坐标、地形体积、边界周长和计算形状因子，将此中山简化为锥形山，再求出底面圆半径和锥形山高度，最终得到锥形山底面圆中心、半径和圆锥高度几何参数，图7是此孤峰山简化并按提取的几何参数以网格形式绘制的图形。

Claims (3)

1.数字地图中典型地形几何参数获取方法，其特征是：包括如下步骤：
步骤101：使用ArcGIS工具，打开一个ASTER GDEM类型的GeoTIFF地图文件；
步骤102：将所有栅格点的高程值作归一化处理，结果保存在DEM_1.tif文件中；归一化处理的方法为：设x为任一栅格点处的量值，该量值是高程、坡度、坡度变率、曲率、地形起伏度中的任一种，xmax是所有栅格点量值的最大值，xmin是所有栅格点量值的最小值，x′是x的归一化值，则x′＝(x-xmin)×255/(xmax-xmin)；
步骤103：计算所有栅格点的坡度值，再作归一化处理，归一化的坡度值保存在Slope_1.tif文件中；
步骤104：计算所有栅格点的坡度变率，再作归一化处理，归一化的坡度变率值保存在PoDuBianLv_1.tif文件中；
步骤105：计算所有栅格点的曲率，并对曲率求绝对值，再作归一化处理，归一化的曲率值保存在Curvature_1.tif文件中；
步骤106：计算所有栅格点的地形起伏度，结果保存在QiFuDu.tif文件中，再对其作归一化处理，归一化的地形起伏度值保存在QiFuDu_1.tif文件中；
步骤107：将归一化的高程、坡度、坡度变率、曲率和地形起伏度做为特性指标，采用ISO(Iterative SelfOrganization)非监督聚类得到地形的初步分类结果，每类地形有一个类别编号，初步分类结果保存在JuLei.tif文件中；
步骤108：使用GDAL库函数打开步骤101的地图文件，打开QiFuDu.tif文件和JuLei.tif文件；
步骤109：依据地形平均起伏度对步骤107初步分类结果进行地形判定分类；
步骤110：用四位正整数标志地形类别及区域，前一位表示地形类别，平原、丘陵、低山、中山、高山分别用1、2、3、4、5表示；后三位表示此种类别地形的区域个数编号，从001开始累积编号，地形判定分类的结果用此种方式表示成数值之后，保存在地形分类结果文件FenLei.tif之中；
步骤111：读取分类结果文件FenLei.tif，根据地形类别及区域标志数字选择待研究地形区域；
步骤112：判断地形类别标志是否为1，若是，转至步骤113；若不是，转至步骤115；
步骤113：将平原地形简化为平面，求取地形类别标志为平原的所有栅格点高程的平均值，记其为平均高程Ea；
步骤114：求取地形类别标志为平原的所有栅格点地形起伏度的平均值，记其为平均起伏度Ha，转至步骤127；
步骤115：对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中；
步骤116：求地形区域的面积S、地形区域的中心C(x0,y0)和地形的体积V：由链表Point_List存储的边界点序列连线形成一个边界多边形，用图形学中的“交点计数法”判断栅格点是否在边界多边形内，统计处在边界多边形内和位于边界多边形边上的栅格点数目num，设单个栅格点所占面积为a，则地形区域的面积S＝num×a；地形区域的中心点C的坐标xi、yi是边界多边形内及边界多边形边上的任一栅格点的坐标；地形的体积hi是边界多边形内及边界多边形边上的任一栅格点的高程；
步骤117：地形类别标志是否为2，若是，转至步骤118；若不是，转至步骤120；
步骤118：将丘陵地形简化为球缺，求取球缺的底面圆半径
步骤119：求取球缺的高度h1：解方程(h1)3+3(r1)2(h1)-6V/π＝0即得球缺的高度转至步骤127；
步骤120：求地形区域的周长L：由链表Point_List存储的边界点序列连线形成一个边界多边形，边界点即为边界多边形的顶点，设多边形顶点个数为n，令第n+1顶点等于第1个顶点，边界多边形任一边的边长|PiPi+1|由边的两个顶点Pi、Pi+1通过两点距离公式求得，则边界多边形的周长即为地形区域的周长L，
步骤121：求地形区域的形状因子f：形状因子f＝4πS/L2；
步骤122：判断形状因子f是否大于或等于0.7，若是，转至步骤123；若不是，转至步骤125；
步骤123：将地形简化为锥形山，求锥形山底面圆的半径
步骤124：求取锥形山的高度转至步骤127；
步骤125：将地形简化为楔形山，求取楔形山底面区域的最小外接矩形；
所示的步骤125中将地形简化为楔形山，求取楔形山底面区域的最小外接矩形，包括以下步骤：
步骤401：读取地形区域边界点序列链表Point_List，令链表Rot_List等于链表Point_List，ang＝3°,定义二维数组Data[30][4]；
步骤402：对链表Rot_List中的多边形各顶点做以地形区域中心C(x0,y0)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Rot_List中；
步骤403：求链表Rot_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；
步骤404：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；
步骤405：将步骤402～404循环29次；
步骤406：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m；
步骤407：读取数组Data中第m行的前3列值，分别用an、a、b记录；令an＝90°-an，定义an为矩形方向角，它是长度为a的矩形边与正东方向的夹角；
步骤408：输出地形底面区域的最小外接矩形边长a、b及矩形方向角an；
步骤126：求取楔形山的高度其中a、b为步骤125求取的外接矩形的两个边长，转至步骤127；
步骤127：输出地形几何参数：若地形为平原，输出平均高程Ea、平均起伏度Ha；若地形为丘陵，输出球缺底面圆中心C(x0,y0)、半径r1和球缺高度h1；若地形为锥形山，输出圆锥底面圆中心C(x0,y0)、半径r2和圆锥高度h2；若地形为楔形山，输出楔形体底面矩形中心C(x0,y0)，边长a、b，矩形方向角an和楔形体高度h3。
2.根据权利要求1所述的数字地图中典型地形几何参数获取方法，其特征是：所示的步骤109中依据地形平均起伏度对步骤107初步分类结果进行地形判定分类，包括以下步骤：
步骤201：读取JuLei.tif文件，遍历地形类别编号，记该编号为n；
步骤202：遍历类别n的地形范围，从QiFuDu.tif文件中读取相对应位置的地形起伏度数据；
步骤203：求取该类别地形区域的地形起伏度平均值U；
步骤204：判断U是否小于或等于20米，若是，转至步骤205；若不是，转至步骤206；
步骤205：该类别地形判定为平原；
步骤206：判断U是否小于或等于200米，若是，转至步骤207；若不是，转至步骤208；
步骤207：该类别地形判定为丘陵；
步骤208：判断U是否小于或等于500米，若是，转至步骤209；若不是，转至步骤210；
步骤209：该类别地形判定为低山；
步骤210：判断U是否小于或等于1500米，若是，转至步骤211；若不是，转至步骤212；
步骤211：该类别地形判定为中山；
步骤212：该类别地形判定为高山。
3.根据权利要求1所述的数字地图中典型地形几何参数获取方法，其特征是：所示的步骤115中对待研究的地形区域进行八邻域边界跟踪，提取出边界点序列存入链表Point_List中，包括以下步骤：
步骤301：基于GDAL读取地形分类结果文件FeiLei.tif；
步骤302：根据待研究的地形区域的四位正整数标志，将其地形类别标志数字记作type；
步骤303：按照从上至下、由左到右的顺序，从文件FeiLei.tif中图像左上角依次读取栅格点的地形类别标志数字，当第一次读到的栅格点的地形类别标志数字等于type时，则将该点存入边界点序列链表Point_List中，并记录该点为边界起点；
步骤304：令当前边界点等于边界起点，令下一个边界点的搜索方向码Directcode＝0；边界点搜索采用八邻域方法，每个栅格点有八个其他的栅格点与其相邻，该栅格点正上方栅格点的搜索方向码记为0，该栅格点右上方、正右方、右下方、正下方、左下方、正左方、左上方栅格点的搜索方向码依次记为1、2、3、4、5、6、7；
步骤305：读取当前边界点沿Directcode方向的搜索点的地形类别标志数字，将其记作value；
步骤306：判断value是否等于type，若是，转至步骤308；若不是，转至步骤307；
步骤307：令Directcode＝Directcode+1，转至步骤305；
步骤308：令当前边界点等于此搜索点；
步骤309：判断当前边界点是否等于边界起点，若是，则边界跟踪结束；若不是，转至步骤310；
步骤310：将当前边界点存入链表Point_List中；
步骤311：令Directcode＝Directcode-2；
步骤312：判断Directcode是否大于或等于0，若是，转至步骤305；若不是，转至步骤313；
步骤313：令Directcode＝Directcode+8，转至步骤305。
CN201610853996.XA 2016-09-27 2016-09-27 数字地图中典型地形几何参数获取方法 Active CN106649466B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610853996.XA CN106649466B (zh) 2016-09-27 2016-09-27 数字地图中典型地形几何参数获取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610853996.XA CN106649466B (zh) 2016-09-27 2016-09-27 数字地图中典型地形几何参数获取方法

Publications (2)

Publication Number Publication Date
CN106649466A CN106649466A (zh) 2017-05-10
CN106649466B true CN106649466B (zh) 2019-09-13

Family

ID=58853401

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610853996.XA Active CN106649466B (zh) 2016-09-27 2016-09-27 数字地图中典型地形几何参数获取方法

Country Status (1)

Country Link
CN (1) CN106649466B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107368839B (zh) * 2017-06-22 2020-05-12 南京师范大学 一种基于dem的断层面的自动提取方法
CN109118939B (zh) * 2017-06-23 2022-02-18 薛富盛 三维地形地图及其制作方法
CN107680154B (zh) * 2017-09-28 2020-01-10 西安电子科技大学 基于视图的体素几何参数提取方法
CN108389255B (zh) * 2018-02-12 2021-05-11 西安电子科技大学 基于分层高程云图的地形几何参数提取方法
CN109633542B (zh) * 2018-11-30 2021-10-15 中交遥感载荷(江苏)科技有限公司 一种用于无人机的局域性地图网络的布置方法
CN109657713A (zh) * 2018-12-11 2019-04-19 武汉大学 一种基于众源路网数据的多因子路网匹配方法及系统
CN109448539B (zh) * 2018-12-14 2021-03-19 深圳市置辰海信科技有限公司 一种基于qt框架的多波束色阶图绘制方法
CN110635856B (zh) * 2019-09-29 2022-04-19 北京电子工程总体研究所 一种计算机模拟在林区地形中进行无线电波通信的方法
CN112102435B (zh) * 2020-09-24 2023-08-01 安徽文香科技股份有限公司 一种几何图形绘制的方法、装置、设备及存储介质
CN113298927A (zh) * 2020-12-14 2021-08-24 阿里巴巴（中国）有限公司 数据处理、应用界面显示及辅助操作方法和装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105469061A (zh) * 2015-08-04 2016-04-06 电子科技大学中山学院 地形特征线提取方法及装置
CN105787457A (zh) * 2016-03-08 2016-07-20 浙江工商大学 一种modis卫星集成dem提高植被分类遥感精度的估算方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105469061A (zh) * 2015-08-04 2016-04-06 电子科技大学中山学院 地形特征线提取方法及装置
CN105787457A (zh) * 2016-03-08 2016-07-20 浙江工商大学 一种modis卫星集成dem提高植被分类遥感精度的估算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Modelling land cover change using Geographic Information Systems, remote sensing, an urbanization-suitability layer and SLEUTH-3r for Groningen, the Netherlands;Venema;《TU/e》;20160831;全文 *
基于数字地图的地形要点选取方法;张德;《火力与指挥控制》;20131031;第38卷(第10期);全文 *

Also Published As

Publication number Publication date
CN106649466A (zh) 2017-05-10

Similar Documents

Publication Publication Date Title
CN106649466B (zh) 数字地图中典型地形几何参数获取方法
Shi et al. Road detection from remote sensing images by generative adversarial networks
CN108389255B (zh) 基于分层高程云图的地形几何参数提取方法
CN112070769B (zh) 一种基于dbscan的分层点云分割方法
CN102289991B (zh) 一种基于视觉变量的地图注记自动分类配置方法
CN105513127A (zh) 基于密度峰值聚类的杆状物规则化三维建模方法及系统
CN110021072B (zh) 面向全息测绘的多平台点云智能处理方法
Shen et al. A polygon aggregation method with global feature preservation using superpixel segmentation
CN108921943A (zh) 一种基于车道级高精度地图的道路三维模型建模方法
Li et al. Intelligent map reader: A framework for topographic map understanding with deep learning and gazetteer
CN103065361A (zh) 三维海岛沙盘实现方法
Yan et al. Description approaches and automated generalization algorithms for groups of map objects
CN106780586A (zh) 一种基于地面激光点云的太阳能潜力评估方法
Walter et al. Automatic interpretation of digital maps
Guo et al. Large-scale and refined green space identification-based sustainable urban renewal mode assessment
Chen et al. A Voronoi interior adjacency-based approach for generating a contour tree
STROBL Segmentation-based terrain classification
CN102142155A (zh) 面向网络交互可视化的地面三维模型数据组织方法
CN113838199A (zh) 一种三维地形生成方法
Zhou GeoAI-Enhanced Techniques to Support Geographical Knowledge Discovery from Big Geospatial Data
Douglas Experiments to locate ridges and channels to create a new type of digital elevation model
Spark et al. Digital Terrain Models and the visualization of structural geology
Wang et al. Learning visual features from figure-ground maps for urban morphology discovery
Minchilli et al. The geographical distribution of nuraghi in north-western Sardinia: Analysis and evaluation of the influence of anthropic and natural factors
Hassan Comparing Geomorphometric Pattern Recognition Methods for Semi-Automated Landform Mapping

Legal Events

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