CN108389255B - 基于分层高程云图的地形几何参数提取方法 - Google Patents

Abstract

本发明涉及一种基于分层高程云图的地形几何参数提取方法，包括地形分类及识别、地形简化、地形几何参数提取过程。针对ASTER GDEM类型的GeoTIFF格式的数字地图中指定的一个401×401DEM格网地形，通过编程处理，首先计算相对高度将地形分为平原、丘陵、低山、中山、高山五种类型；然后对高程分层和颜色映射，得到指定地形的分层高程云图；之后利用图像处理中的八邻域边界跟踪算法，求得地形底面边界、底面面积和底面的中心，再识别代表最大高程区间的颜色以求取地形最高区域的平均高程；最后结合地形分类、底面形状因子计算、坡度计算，求取地形简化几何模型平原、球冠、圆锥体、楔形体的几何参数。本发明具有方法简单、计算量小、实时性和实用好的优点。

Description

基于分层高程云图的地形几何参数提取方法

技术领域

本发明属于地理信息系统(GIS)和计算机图学、图像领域，涉及地形地貌的分类方法和地形几何参数提取及简化模型的建立方法，具体地说就是基于分层高程云图的地形几何参数提取方法，可用于地形类型自动判别、电波传播近似计算和地形模型快速实时显示。

背景技术

基于地理信息系统(GIS)的数字地图以地理信息系统和数字制图系统为技术支撑，其功能和作用随着GIS技术和数字制图系统的发展而不断进步。数字地图具体含义是指制图元素和离散数据在可读取存储介质上的广义和有序集合，同时这些集合与数据在某个系统中具有确切的位置和属性。与传统的地图比较有很多优点：数字地图可存放大量的地形地貌特征信息，并且这些信息传播快，同时信息更替方便，使用与之匹配的遥感(RS)，可实现任何时候都能改变并更替时空信息，建立数字地图的一个很关键目的就是查阅，数字地图信息查阅更加多样化，通过PC端或其他电子设备就能实现放大、缩小，同时动态地比例尺显示，使信息检阅更方便快捷，它的一些特性已经远远超越了过去的纸质地图。

数字高程模型(Digital Elevation Model，简称DEM)在获取地理信息过程中，提供了全新的数据源，地表的各种特征在DEM中也能充足并且真实的表现。数字高程模型采用数字表达的方式对地球上的地形进行表示，其具体是由一些地面上点的平面位置以及点的高程信息，采用一定的结构进行组合，完成表达地形特征的空间布局模型。DEM在日常生活、科学研究和工程中具有极大的使用价值，例如：(1)地图产品的生成，如对数字地图处理可以产生印刷地图、电子地图、导航地图、断面图、天气预报图、国家防汛图、灾害分析图等地图产品；(2)工程设计、道路规划和城市规划等，利用DEM可以获取坐标、角度、长度、面积、体积及土方、高度、密度、梯度等设计和规划所需的数据；(3)航空工业中，飞机的飞行和地面导航的建立，都可以用DEM作为基本的数据源，它们都需要模拟真实地面；(4)可应用于无线电波传播计算，移动通信、野战通信等的信号覆盖范围和通信距离需要电波传播计算来确定，在电波传播过程中，山地、丘陵、平原等传播环境对信号损耗有着重要的影响，而利用DEM可以实现较为准确的电波传播计算；(5)在研究机器人移动技术的时候，DEM可以提供大量的室外地形数据，包括地形起伏度、坡度等等地形属性信息。

数字地图的地形信息提取研究内容包括对坡度、坡向、起伏度、山顶点、山脊线、山谷线、山底线的提取，对地形的分类，以及对所提取的地形进行简化。论文“基于DEM山脉线的提取及地貌类型方法判别研究”(西安建筑科技大学硕士学位论文，2015，作者：雷红涛)将多尺度分析与图像融合理论相结合，设计了一种山脉线自动提取方法，该方法以不同参数的高斯核函数对规则格网DEM数据做多尺度分解，先抽取各尺度下的山脊线，再构造合理的山脊线融合规则，形成能体现宏观地形特征的山脉线；同时借鉴图像分割及连通域分析等图像处理方法，采用地形因子建立盆地和高原的形态描述模型，并结合地势起伏度和海拔高度这两个全国普适性指标，设计了平原、丘陵、山地、盆地和高原地貌类型的自动判别方法；论文“基于DEM的防界线提取与路径规划方法研究”(西安建筑科技大学硕士学位论文，2014，作者：高黄玮))通过对坡度、坡向等微观地形因子进行组合运用，设计了一种防界线提取算法，该方法着眼于防界线局部的微观形态特征，利用滑动窗口分析山体曲面的坡度变化情况，筛选出坡度改变较大的防界点并结合其它约束条件形成防界线，同时针对现有路径规划方法存在环境代价模型不够准确、算法未考虑路面宽度因素等问题，从影响车辆通行的地形要素出发，对路径规划中的环境代价模型进行了重新研究，给出了改进的宽路径搜索方案；论文“一种新的山脊线和山谷线自动提取方法”(中国图像图形学报,2001年6月，1230-1234，作者：陈永良,刘大有)采纳DEM每个点各个方向的剖面线特征并识别山脊点和山谷点，对这些点进行一番筛选并连接筛选后的点，完成山脊线和山谷线的提取，该方法并不能实现完全自动化，存在一定的局限性；论文“提取山脊线与山谷线的一种新方法”(武汉大学学报·信息科学版,2001年3月，247-251，作者：黄培之)针对现有的提取山脊线和山谷线含有其几何特性和物理特性，将地形表面流水模拟分析法和等高线几何分析法进行有机结合，克服各自所包含的缺点，设计一种提取山脊线和山谷线的方法，该方法中提取特征会用到阈值的选取、网格的确定，如果控制不好，则会引起较大误差；论文“基于数字地图的地形信息提取技术研究”(西安电子科技大学硕士学位论文，2015，作者：李杨)从数学的角度利用模糊聚类的思想对地形特征提取进行了研究，但是此种方法需要借助第三方软件来实现，不能自动识别地形分类、地形特征信息提取，在聚类过程中需要输入多个初始化参数，如期望得到的聚类数、初始聚类中心、一个类别中最小样本数目、类的分散程度的参数(如标准差、方差)、类间距离的参数(如最小距离)、每次迭代允许合并的最大聚类对数、允许迭代的次数等，人工干预性比较大，如果参数选取不当，很难到达预期的效果。可以看出，已有研究成果存在以下问题：(1)研究仅涉及到地形特殊点(如山顶点、鞍部点)、特殊线(山脊线、山谷线、山脉线、山脚线、防界线)等局部几何特征的识别，缺乏对地形整体几何形状识别、分割、简化和几何参数提取的研究；(2)研究方法因有人工设置的一些参数，自动化程度偏低；(3)对于地形实时应用的场合，未有地形距离范围限定的研究，方法的实用性难于保证。

发明内容

本发明的目的在于克服上述现有研究中存在的问题，提供一种基于分层高程云图的地形几何参数提取方法，以满足地形场景计算和地形模型快速实时显示的需要。

本发明的目的是这样实现的：基于分层高程云图的地形几何参数提取方法，其特征是：至少包括如下步骤：

步骤101：编程打开一个地区的ASTER GDEM类型的GeoTIFF地图文件；

步骤102：设置参照点的位置(m₀,n₀)和观察方向，(m₀,n₀)为图像坐标，观察方向包括东经度观察方向和北纬度观察方向，东经度观察方向分为经度值增加方向和经度值减小方向，北纬度观察方向分为纬度值增加方向和纬度值减小方向；

步骤103：确定待读取地形的起点坐标(m₁,n₁)，起点为图像左上角位置，若观察方向为东经度值减小方向和北纬度值增加方向，则m₁＝m₀-400，n₁＝n₀-400；若观察方向为东经度值增加方向和北纬度值增加方向，则m₁＝m₀，n₁＝n₀-400；若观察方向为东经度值减小方向和北纬度值减小方向，则m₁＝m₀-400，n₁＝n₀；若观察方向为东经度值增加方向和北纬度值减小方向，则m₁＝m₀，n₁＝n₀；

步骤104：从点(m₁,n₁)开始，按照先水平向右方向、后竖直向下方向顺序读取401×401个栅格点，栅格点依次存入数组A(i,j)，i,j＝0,1,2,…,400，栅格点存放在数组中的行位置即对应下标i、列位置即对应下标j与其在401×401图像中的行、列位置对应，数组元素为Point类，Point类记录了每个栅格点的行数、列数、经度值、纬度值、高程值、坡度值、高程分层颜色值和是否为地形底面边界点的标志值，栅格点坡度值的初值为0，栅格点高程分层颜色值的初值为蓝色，栅格点是否为地形底面边界点的标志值初值为0；

步骤105：遍历数组A(i,j)，找出所有栅格点高程的最小值h_min和最大值h_max，计算地形的相对高度H＝h_max-h_min；

步骤106：根据H值确定地形类型参数k_d和高程分层高度Δh；

步骤107：判断k_d是否为1，若是，转至步骤108；若不是，转至步骤110；

步骤108：遍历数组A(i,j)，求所有栅格点高程之和H_sum，则平原地形的平均高程H_a＝H_sum/160801；

步骤109：遍历数组A(i,j)并访问打开的GeoTIFF地图文件，对于A(i,j)中的每个栅格点，以其为中心取周围60×60栅格点为统计范围，求取该范围内栅格点的最小高程和最大高程，则最大高程减去最小高程即为该栅格点的地形起伏度U_d，数组A(i,j)中所有栅格点的U_d之和除以160801即为平原地形的平均起伏度U_da，转至步骤123；

步骤110：根据h_min、h_max和Δh计算并形成分层高程云图；

步骤111：对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记；

步骤112：按照先行后列顺序遍历数组A(i,j)，对于有地形底面边界点标志值为-1的一行栅格点，记录该行中边界点个数m′、边界点及边界点之间的栅格点数目n′、高程分层颜色值为红色的栅格点数目k′及高程值之和h′，将所有具有边界点的行的m′相加得到地形底面边界点总数B_n，将所有具有边界点的行的n′相加得到地形底面边界包围的栅格点总数G_n，将所有具有边界点的行的k′相加、h′相加分别得到地形最高区域含有的栅格点总数k和地形最高区域的高程值之和H_t，则该地形最高区域的平均高程H_ta＝H_t/k；

步骤113：设单个栅格点所占面积为a，则地形底面的面积S＝G_n×a；地形底面的中心点C的坐标

x_i、y_i是地形底面边界及边界内的任一栅格点的经度、纬度坐标，z_j是地形底面边界上栅格点的高程值；

步骤114：判断k_d是否为2，若是，转至步骤115；若不是，转至步骤116；

步骤115：将丘陵地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤116：求地形底面的形状因子

步骤117：判断f是否大于或等于0.65，若是，转至步骤118；若不是，转至步骤122；

步骤118：求地形底面边界围成区域的平均坡度sl_a；

步骤119：判断sl_a是否大于15°，若是，转至步骤120；若不是，转至步骤121；

步骤120：将地形简化为圆锥形山，圆锥底面圆半径为

圆锥高度为h₂＝H_ta-h_min，转至步骤123；

步骤121：将山体地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤122：将地形简化为楔形山，楔形山的几何模型为楔形体，该楔形体是一个底面为矩形，四个侧面中两个相对侧面为等腰三角形、另两个相对侧面为矩形的几何体；求楔形体底面矩形边长l和w，矩形方向角an和楔形体高度h₃；

步骤123：输出地形几何参数：若地形为平原，输出平均高程H_a、平均起伏度U_da；若地形为丘陵，输出球冠底面圆中心C(x₀,y₀,z₀)、底面圆半径r₁和球冠高度h₁；若地形为锥形山，输出圆锥底面圆中心C(x₀,y₀,z₀)、底面圆半径r₂和圆锥高度h₂；若地形为楔形山，输出楔形体底面矩形中心C(x₀,y₀,z₀)，底面矩形边长l和w，矩形方向角an和楔形体高度h₃。

所述的步骤106中，根据H值确定地形类型参数k_d和高程分层高度Δh，包括以下步骤：

步骤201：判断H是否小于或等于20，若是，转至步骤202；若不是，转至步骤203；

步骤202：该地形判定为平原，令地形类型参数k_d＝1，高程分层高度Δh＝20；

步骤203：判断H是否小于或等于200，若是，转至步骤204；若不是，转至步骤205；

步骤204：该地形判定为丘陵，令地形类型参数k_d＝2，高程分层高度Δh＝30；

步骤205：判断H是否小于或等于500，若是，转至步骤206；若不是，转至步骤207；

步骤206：该地形判定为低山，令地形类型参数k_d＝3，高程分层高度Δh＝40；

步骤207：判断H是否小于或等于1000，若是，转至步骤208；若不是，转至步骤209；

步骤208：该地形判定为中山，令地形类型参数k_d＝4，高程分层高度Δh＝50；

步骤209：该地形判定为高山，令地形类型参数k_d＝5，高程分层高度Δh＝50。

所述的步骤110中，根据h_min、h_max和Δh计算并形成分层高程云图，包括以下步骤：

步骤301：根据h_min、h_max和Δh确定地形高程分层数L_n，L_n＝[(h_max-h_min)/Δh]+1，符号[]表示取整数；

步骤302：将区间(h_min，h_max]等分为L_n个小区间并按值由小到大排列，若L_n＝2，则值较小的小区间分配的颜色为青色，值较大的小区间分配的颜色为红色；若L_n≥3，则值较小的小区间分配的颜色为青色，值最大的小区间分配的颜色为红色，中间位置的小区间按照绿、黄、橙色调依次与小区间对应分配颜色，每个小区间分配一个颜色；

步骤303：遍历数组A(i,j)，判断每个栅格点高程值属于哪个小区间，将该小区间分配的颜色值赋值给Point类的栅格点高程分层颜色变量，数组A(i,j)记录了所关注地形的分层高程云图。

所述的步骤111中，对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记，包括以下步骤：

步骤401：根据步骤302确定地形底面边界就是L_n个小区间中值最小的区间所对应的栅格点围成的封闭区域的外轮廓，地形底面边界栅格点的颜色为青色；

步骤402：按照从下往上、从左往右的顺序遍历数组A(i,j)，遇到的高程分层颜色为青色的第一个栅格点即为地形底面边界的起点，记该点为P_s，其Point类存储的栅格点行数、列数分别记为i_s、j_s，令该点Point类中的地形底面边界点的标志值为-1；

步骤403：设定八邻域搜索的方向编码，记搜索到的地形底面边界上的当前点为P_c，其栅格点行数、列数分别为i_c、j_c，将指向行数、列数分别为i_c-1、j_c-1的栅格点的搜索方向编码规定为0，将指向行数、列数分别为i_c-1、j_c的栅格点的搜索方向编码规定为1，将指向行数、列数分别为i_c-1、j_c+1的栅格点的搜索方向编码规定为2，将指向行数、列数分别为i_c、j_c+1的栅格点的搜索方向编码规定为3，将指向行数、列数分别为i_c+1、j_c+1的栅格点的搜索方向编码规定为4，将指向行数、列数分别为i_c+1、j_c的栅格点的搜索方向编码规定为5，将指向行数、列数分别为i_c+1、j_c-1的栅格点的搜索方向编码规定为6，将指向行数、列数分别为i_c、j_c-1的栅格点的搜索方向编码规定为7；

步骤404：令底面边界上搜索到的当前栅格点P_c＝P_s，搜索方向码变量S_d＝0；

步骤405：查找P_c沿S_d方向上的栅格点P_n的高程分层颜色值C_n；

步骤406：判断C_n是否为青色，若是，转至步骤407；若不是，转至步骤409；

步骤407：将数组A(i,j)中P_n位置的栅格点的地形底面边界点的标志值置为-1，令P_c＝P_n；

步骤408：判断P_n是否等于P_s，若不是，转至步骤412；若是，则地形底面边界跟踪结束；

步骤409：令S_d＝S_d+1；

步骤410：判断S_d是否大于或等于8，若是，转至步骤411；若不是，转至步骤405；

步骤411：令S_d＝S_d-8，转至步骤405；

步骤412：令S_d＝S_d-2；

步骤413：判断S_d是否小于0，若是，转至步骤414；若不是，转至步骤405；

步骤414：令S_d＝S_d+8，转至步骤405。

所述的步骤118中，求地形底面边界围成区域的平均坡度sl_a，包括以下步骤：

步骤501：按照先行后列顺序遍历数组A(i,j)，对于地形底面边界点标志值为-1的栅格点以及一行内处在两个边界点之间的栅格点，逐个计算这些栅格点的坡度值sl：设待计算的栅格点高程为h_e，该栅格点左上方的栅格点高程为h_a，该栅格点正上方的栅格点高程为h_b，该栅格点右上方的栅格点高程为h_c，该栅格点正左方的栅格点高程为h_d，该栅格点正右方的栅格点高程为h_f，该栅格点左下方的栅格点高程为h_g，该栅格点正下方的栅格点高程为h_h，该栅格点右下方的栅格点高程为h_i，令f_x＝(h_i-h_g+2(h_f-h_d)+h_c-h_a)/(8×cellsize)，f_y＝(h_a-h_g+2(h_b-h_h)+h_c-h_i)/(8×cellsize)，cellsize表示DEM格网的分辨率，则该栅格点的坡度

步骤502：对边界区域内的每个栅格点，将其Point类中坡度变量值置为各自的计算值sl；

步骤503：遍历数组A(i,j)，计算地形底面边界区域内所有栅格点坡度值之和sl_sum，则地形底面边界围成区域的平均坡度sl_a＝sl_sum/G_n。

所述的步骤122中，求楔形体底面矩形边长l和w，矩形方向角an和楔形体高度h₃，包括以下步骤：

步骤601：建立以(m0,n0)为原点，水平向右为X轴正向、竖直向下为Y轴正向的局部坐标系；

步骤602：形成以直角坐标表示点的地形边界点序列链表Point_List：先按行从上往下、列从左往右遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，然后按行从下往上、列从右往左遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，该链表存储的边界点连线形成一个多边形；

步骤603：令ang＝3°,定义二维数组Data[30][4]；

步骤604：对链表Point_List中的多边形各顶点做以点(x0,y0)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Point_List中；

步骤605：求链表Point_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；

步骤606：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；

步骤607：将步骤604～606循环29次；

步骤608：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m′；

步骤609：读取数组Data中第m′行的前3列值，分别用an、l、w记录；令an＝90°-an，定义an为楔形体底面矩形方向角，它是长度为l的矩形边与正东方向的夹角，w为底面矩形的另一边长；

步骤610：楔形体高度为h₃＝H_ta-h_min。

本发明有如下优点：

(1)提出了基于分层高程云图的典型地形整体几何形状自动识别的方法，此方法简单、实用；

(2)以参考点一侧设定的任务范围内的地形为研究对象，方法计算量小、实时性好；

(3)通过地形简化，解决了基于数字地图的典型地形快速几何建模问题。

附图说明

图1是本发明的总流程图；

图2是地形分类判定的流程图；

图3是地形分层高程云图形成的流程图；

图4是地形八邻域边界跟踪流程图；

图5是地形平均坡度计算流程图；

图6是楔形山底面区域最小外接矩形求取的流程图；

图7是某座山的实际形状；

图8是图7所示山对应的分层高程云图；

图9是图7所示山简化并提取几何参数后再显示的结果。

具体实施方式

本发明所读取的数据源，是一种数字高程模型(DEM)，数据格式是ASTER GDEM(Advanced Spaceborn Thermal Emission and Reflection Radiometer Global DigitalElevation Model)类型的GeoTIFF遥感影像数据(.tif文件)，该图像栅格点对应的空间分辨率约为30m×30m。在指定图像上参照点位置(比如说通信车、飞行器位置)和观察方向后，从参照点开始沿观察方向读取一个401×401DEM格网作为地形几何参数提取的区域。本发明地形几何参数提取的方法包括地形分类及识别、地形简化、地形几何参数提取几个过程。通过编程处理，首先计算相对高度将地形分为平原、丘陵、低山、中山、高山五种类型；然后对高程分层和颜色映射，得到指定地形的分层高程云图；之后利用图像处理中的八邻域边界跟踪算法，求得地形底面边界、底面面积和底面的中心，再识别代表最大高程区间的颜色以求取地形最高区域的平均高程；最后结合地形分类、底面形状因子计算、坡度计算，求取地形简化几何模型平原、球冠、圆锥体、楔形体的几何参数。

参照图1，本发明的基于分层高程云图的地形几何参数提取方法包括如下步骤：

步骤101：编程打开一个地区的ASTER GDEM类型的GeoTIFF地图文件；

步骤102：设置参照点的位置(m₀,n₀)和观察方向，(m₀,n₀)为图像坐标，观察方向包括东经度观察方向和北纬度观察方向，东经度观察方向分为经度值增加方向和经度值减小方向，北纬度观察方向分为纬度值增加方向和纬度值减小方向；

步骤103：确定待读取地形的起点坐标(m₁,n₁)，起点为图像左上角位置，若观察方向为东经度值减小方向和北纬度值增加方向，则m₁＝m₀-400，n₁＝n₀-400；若观察方向为东经度值增加方向和北纬度值增加方向，则m₁＝m₀，n₁＝n₀-400；若观察方向为东经度值减小方向和北纬度值减小方向，则m₁＝m₀-400，n₁＝n₀；若观察方向为东经度值增加方向和北纬度值减小方向，则m₁＝m₀，n₁＝n₀；

步骤104：从点(m₁,n₁)开始，按照先水平向右方向、后竖直向下方向顺序读取401×401个栅格点，栅格点依次存入数组A(i,j)，i,j＝0,1,2,…,400，栅格点存放在数组中的行位置即对应下标i、列位置即对应下标j与其在401×401图像中的行、列位置对应，数组元素为Point类，Point类记录了每个栅格点的行数、列数、经度值、纬度值、高程值、坡度值、高程分层颜色值和是否为地形底面边界点的标志值，栅格点坡度值的初值为0，栅格点高程分层颜色值的初值为蓝色，栅格点是否为地形底面边界点的标志值初值为0；

步骤105：遍历数组A(i,j)，找出所有栅格点高程的最小值h_min和最大值h_max，计算地形的相对高度H＝h_max-h_min；

步骤106：根据H值确定地形类型参数k_d和高程分层高度Δh；

步骤106中根据H值确定地形类型参数k_d和高程分层高度Δh，参照图2，包括以下步骤：

步骤201：判断H是否小于或等于20，若是，转至步骤202；若不是，转至步骤203；

步骤202：该地形判定为平原，令地形类型参数k_d＝1，高程分层高度Δh＝20；

步骤203：判断H是否小于或等于200，若是，转至步骤204；若不是，转至步骤205；

步骤204：该地形判定为丘陵，令地形类型参数k_d＝2，高程分层高度Δh＝30；

步骤205：判断H是否小于或等于500，若是，转至步骤206；若不是，转至步骤207；

步骤206：该地形判定为低山，令地形类型参数k_d＝3，高程分层高度Δh＝40；

步骤207：判断H是否小于或等于1000，若是，转至步骤208；若不是，转至步骤209；

步骤208：该地形判定为中山，令地形类型参数k_d＝4，高程分层高度Δh＝50；

步骤209：该地形判定为高山，令地形类型参数k_d＝5，高程分层高度Δh＝50。

步骤107：判断k_d是否为1，若是，转至步骤108；若不是，转至步骤110；

步骤108：遍历数组A(i,j)，求所有栅格点高程之和H_sum，则平原地形的平均高程H_a＝H_sum/160801；

步骤109：遍历数组A(i,j)并访问打开的GeoTIFF地图文件，对于A(i,j)中的每个栅格点，以其为中心取周围60×60栅格点为统计范围，求取该范围内栅格点的最小高程和最大高程，则最大高程减去最小高程即为该栅格点的地形起伏度U_d，数组A(i,j)中所有栅格点的U_d之和除以160801即为平原地形的平均起伏度U_da，转至步骤123；

步骤110：根据h_min、h_max和Δh计算并形成分层高程云图；

步骤110中根据h_min、h_max和Δh计算并形成分层高程云图，参照图3，包括以下步骤：

步骤301：根据h_min、h_max和Δh确定地形高程分层数L_n，L_n＝[(h_max-h_min)/Δh]+1，符号[]表示取整数；

步骤302：将区间(h_min，h_max]等分为L_n个小区间并按值由小到大排列，若L_n＝2，则值较小的小区间分配的颜色为青色，值较大的小区间分配的颜色为红色；若L_n≥3，则值较小的小区间分配的颜色为青色，值最大的小区间分配的颜色为红色，中间位置的小区间按照绿、黄、橙色调依次与小区间对应分配颜色，每个小区间分配一个颜色；

步骤303：遍历数组A(i,j)，判断每个栅格点高程值属于哪个小区间，将该小区间分配的颜色值赋值给Point类的栅格点高程分层颜色变量，数组A(i,j)记录了所关注地形的分层高程云图。

步骤111：对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记；

步骤111中对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记，参照图4，包括以下步骤：

步骤401：根据步骤302确定地形底面边界就是L_n个小区间中值最小的区间所对应的栅格点围成的封闭区域的外轮廓，地形底面边界栅格点的颜色为青色；

步骤402：按照从下往上、从左往右的顺序遍历数组A(i,j)，遇到的高程分层颜色为青色的第一个栅格点即为地形底面边界的起点，记该点为P_s，其Point类存储的栅格点行数、列数分别记为i_s、j_s，令该点Point类中的地形底面边界点的标志值为-1；

步骤403：设定八邻域搜索的方向编码，记搜索到的地形底面边界上的当前点为P_c，其栅格点行数、列数分别为i_c、j_c，将指向行数、列数分别为i_c-1、j_c-1的栅格点的搜索方向编码规定为0，将指向行数、列数分别为i_c-1、j_c的栅格点的搜索方向编码规定为1，将指向行数、列数分别为i_c-1、j_c+1的栅格点的搜索方向编码规定为2，将指向行数、列数分别为i_c、j_c+1的栅格点的搜索方向编码规定为3，将指向行数、列数分别为i_c+1、j_c+1的栅格点的搜索方向编码规定为4，将指向行数、列数分别为i_c+1、j_c的栅格点的搜索方向编码规定为5，将指向行数、列数分别为i_c+1、j_c-1的栅格点的搜索方向编码规定为6，将指向行数、列数分别为i_c、j_c-1的栅格点的搜索方向编码规定为7；

步骤404：令底面边界上搜索到的当前栅格点P_c＝P_s，搜索方向码变量S_d＝0；

步骤405：查找P_c沿S_d方向上的栅格点P_n的高程分层颜色值C_n；

步骤406：判断C_n是否为青色，若是，转至步骤407；若不是，转至步骤409；

步骤407：将数组A(i,j)中P_n位置的栅格点的地形底面边界点的标志值置为-1，令P_c＝P_n；

步骤408：判断P_n是否等于P_s，若不是，转至步骤412；若是，则地形底面边界跟踪结束；

步骤409：令S_d＝S_d+1；

步骤410：判断S_d是否大于或等于8，若是，转至步骤411；若不是，转至步骤405；

步骤411：令S_d＝S_d-8，转至步骤405；

步骤412：令S_d＝S_d-2；

步骤413：判断S_d是否小于0，若是，转至步骤414；若不是，转至步骤405；

步骤414：令S_d＝S_d+8，转至步骤405。

步骤112：按照先行后列顺序遍历数组A(i,j)，对于有地形底面边界点标志值为-1的一行栅格点，记录该行中边界点个数m′、边界点及边界点之间的栅格点数目n′、高程分层颜色值为红色的栅格点数目k′及高程值之和h′，将所有具有边界点的行的m′相加得到地形底面边界点总数B_n，将所有具有边界点的行的n′相加得到地形底面边界包围的栅格点总数G_n，将所有具有边界点的行的k′相加、h′相加分别得到地形最高区域含有的栅格点总数k和地形最高区域的高程值之和H_t，则该地形最高区域的平均高程H_ta＝H_t/k；

步骤113：设单个栅格点所占面积为a，则地形底面的面积S＝G_n×a；地形底面的中心点C的坐标

xi、yi是地形底面边界及边界内的任一栅格点的经度、纬度坐标，z_j是地形底面边界上栅格点的高程值；

步骤114：判断k_d是否为2，若是，转至步骤115；若不是，转至步骤116；

步骤115：将丘陵地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤116：求地形底面的形状因子

步骤117：判断f是否大于或等于0.65，若是，转至步骤118；若不是，转至步骤122；

步骤118：求地形底面边界围成区域的平均坡度sl_a；

步骤118中求地形底面边界围成区域的平均坡度sl_a，参照图5，包括以下步骤：

步骤501：按照先行后列顺序遍历数组A(i,j)，对于地形底面边界点标志值为-1的栅格点以及一行内处在两个边界点之间的栅格点，逐个计算这些栅格点的坡度值sl：设待计算的栅格点高程为h_e，该栅格点左上方的栅格点高程为h_a，该栅格点正上方的栅格点高程为h_b，该栅格点右上方的栅格点高程为h_c，该栅格点正左方的栅格点高程为h_d，该栅格点正右方的栅格点高程为h_f，该栅格点左下方的栅格点高程为h_g，该栅格点正下方的栅格点高程为h_h，该栅格点右下方的栅格点高程为h_i，令f_x＝(h_i-h_g+2(h_f-h_d)+h_c-h_a)/(8×cellsize)，f_y＝(h_a-h_g+2(h_b-h_h)+h_c-h_i)/(8×cellsize)，cellsize表示DEM格网的分辨率，则该栅格点的坡度

步骤502：对边界区域内的每个栅格点，将其Point类中坡度变量值置为各自的计算值sl；

步骤503：遍历数组A(i,j)，计算地形底面边界区域内所有栅格点坡度值之和sl_sum，则地形底面边界围成区域的平均坡度sl_a＝sl_sum/G_n。

步骤119：判断sl_a是否大于15°，若是，转至步骤120；若不是，转至步骤121；

步骤120：将地形简化为圆锥形山，圆锥底面圆半径为

圆锥高度为h₂＝H_ta-h_min，转至步骤123；

步骤121：将山体地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤122：将地形简化为楔形山，楔形山的几何模型为楔形体，该楔形体是一个底面为矩形，四个侧面中两个相对侧面为等腰三角形、另两个相对侧面为矩形的几何体；求楔形体底面矩形边长l和w，矩形方向角an和楔形体高度h₃；

步骤122中求楔形体底面矩形边长l和w，矩形方向角an和楔形体高度h₃，参照图6，包括以下步骤：

步骤601：建立以(m₀,n₀)为原点，水平向右为X轴正向、竖直向下为Y轴正向的局部坐标系；

步骤602：形成以直角坐标表示点的地形边界点序列链表Point_List：先按行从上往下、列从左往右遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，然后按行从下往上、列从右往左遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，该链表存储的边界点连线形成一个多边形；

步骤603：令ang＝3°,定义二维数组Data[30][4]；

步骤604：对链表Point_List中的多边形各顶点做以点(x0,y0)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Point_List中；

步骤605：求链表Point_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；

步骤606：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；

步骤607：将步骤604～606循环29次；

步骤608：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m′；

步骤609：读取数组Data中第m′行的前3列值，分别用an、l、w记录；令an＝90°-an，定义an为楔形体底面矩形方向角，它是长度为l的矩形边与正东方向的夹角，w为底面矩形的另一边长；

步骤610：楔形体高度为h₃＝H_ta-h_min。

步骤123：输出地形几何参数：若地形为平原，输出平均高程H_a、平均起伏度U_da；若地形为丘陵，输出球冠底面圆中心C(x₀,y₀,z₀)、底面圆半径r₁和球冠高度h₁；若地形为锥形山，输出圆锥底面圆中心C(x₀,y₀,z₀)、底面圆半径r₂和圆锥高度h₂；若地形为楔形山，输出楔形体底面矩形中心C(x₀,y₀,z₀)，底面矩形边长l和w，矩形方向角an和楔形体高度h₃。

仿真实例

利用本发明对国内某座山的地形进行几何参数提取。图7是读取该山的DEM数据并以网格形式显示的山的实际形状。通过编程读取该山所在的ASTER GDEM类型的GeoTIFF格式的数字地图文件中401×401格网，首先计算相对高度确定该地形为高山类型；然后对其高程分层和颜色映射，得到图8所示的分层高程云图；之后利用图像处理中的八邻域边界跟踪算法，求得地形底面边界、底面面积和底面的中心，再识别代表最大高程区间的颜色以求取地形最高区域的平均高程；最后经过底面形状因子计算确定该山为楔形山，求取楔形体底面矩形边长和方向角、楔形体高度，从而得到该山简化模型—楔形体的几何参数。图9是该山简化后并按提取的几何参数绘制的图形。

Claims (5)

1.基于分层高程云图的地形几何参数提取方法，其特征是：包括如下步骤：

步骤101：编程打开一个地区的ASTER GDEM类型的GeoTIFF地图文件；

步骤102：设置参照点的位置(m₀,n₀)和观察方向，(m₀,n₀)为图像坐标，观察方向包括东经度观察方向和北纬度观察方向，东经度观察方向分为经度值增加方向和经度值减小方向，北纬度观察方向分为纬度值增加方向和纬度值减小方向；

步骤103：确定待读取地形的起点坐标(m₁,n₁)，起点为图像左上角位置，若观察方向为东经度值减小方向和北纬度值增加方向，则m₁＝m₀-400，n₁＝n₀-400；若观察方向为东经度值增加方向和北纬度值增加方向，则m₁＝m₀，n₁＝n₀-400；若观察方向为东经度值减小方向和北纬度值减小方向，则m₁＝m₀-400，n₁＝n₀；若观察方向为东经度值增加方向和北纬度值减小方向，则m₁＝m₀，n₁＝n₀；

步骤104：从点(m₁,n₁)开始，按照先水平向右方向、后竖直向下方向顺序读取401×401个栅格点，栅格点依次存入数组A(i,j)，i,j＝0,1,2,…,400，栅格点存放在数组中的行位置即对应下标i、列位置即对应下标j与其在401×401图像中的行、列位置对应，数组元素为Point类，Point类记录了每个栅格点的行数、列数、经度值、纬度值、高程值、坡度值、高程分层颜色值和是否为地形底面边界点的标志值，栅格点坡度值的初值为0，栅格点高程分层颜色值的初值为蓝色，栅格点是否为地形底面边界点的标志值初值为0；

步骤105：遍历数组A(i,j)，找出所有栅格点高程的最小值h_min和最大值h_max，计算地形的相对高度H＝h_max-h_min；

步骤106：根据H值确定地形类型参数k_d和高程分层高度Δh；

步骤107：判断k_d是否为1，若是，转至步骤108；若不是，转至步骤110；

步骤108：遍历数组A(i,j)，求所有栅格点高程之和H_sum，则平原地形的平均高程H_a＝H_sum/160801；

步骤109：遍历数组A(i,j)并访问打开的GeoTIFF地图文件，对于A(i,j)中的每个栅格点，以其为中心取周围60×60栅格点为统计范围，求取该范围内栅格点的最小高程和最大高程，则最大高程减去最小高程即为该栅格点的地形起伏度U_d，数组A(i,j)中所有栅格点的U_d之和除以160801即为平原地形的平均起伏度U_da，转至步骤123；

步骤110：根据h_min、h_max和Δh计算并形成分层高程云图；

步骤111：对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记；

步骤112：按照先行后列顺序遍历数组A(i,j)，对于有地形底面边界点标志值为-1的一行栅格点，记录该行中边界点个数m′、边界点及边界点之间的栅格点数目n′、高程分层颜色值为红色的栅格点数目k′及高程值之和h′，将所有具有边界点的行的m′相加得到地形底面边界点总数B_n，将所有具有边界点的行的n′相加得到地形底面边界包围的栅格点总数G_n，将所有具有边界点的行的k′相加、h′相加分别得到地形最高区域含有的栅格点总数k和地形最高区域的高程值之和H_t，则该地形最高区域的平均高程H_ta＝H_t/k；

步骤113：设单个栅格点所占面积为a，则地形底面的面积S＝G_n×a；地形底面的中心点C的坐标

x_i、y_i是地形底面边界及边界内的任一栅格点的经度、纬度坐标，z_j是地形底面边界上栅格点的高程值；

步骤114：判断k_d是否为2，若是，转至步骤115；若不是，转至步骤116；

步骤115：将丘陵地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤116：求地形底面的形状因子

步骤117：判断f是否大于或等于0.65，若是，转至步骤118；若不是，转至步骤122；

步骤118：求地形底面边界围成区域的平均坡度sl_a；

步骤119：判断sl_a是否大于15°，若是，转至步骤120；若不是，转至步骤121；

步骤120：将地形简化为圆锥形山，圆锥底面圆半径为

圆锥高度为h₂＝H_ta-h_min，转至步骤123；

步骤121：将山体地形简化为球冠，球冠底面圆半径为

球冠高度为h₁＝H_ta-h_min，转至步骤123；

步骤122：将地形简化为楔形山，楔形山的几何模型为楔形体，该楔形体是一个底面为矩形，四个侧面中两个相对侧面为等腰三角形、另两个相对侧面为矩形的几何体；求楔形体底面矩形边长l和w，楔形体底面矩形方向角an和楔形体高度h₃；

步骤123：输出地形几何参数：若地形为平原，输出平均高程H_a、平均起伏度U_da；若地形为丘陵，输出球冠底面圆中心C(x₀,y₀,z₀)、底面圆半径r₁和球冠高度h₁；若地形为锥形山，输出圆锥底面圆中心C(x₀,y₀,z₀)、底面圆半径r₂和圆锥高度h₂；若地形为楔形山，输出楔形体底面矩形中心C(x₀,y₀,z₀)，底面矩形边长l和w，楔形体底面矩形方向角an和楔形体高度h₃；

所示的步骤122中，求楔形体底面矩形边长l和w，楔形体底面矩形方向角an和楔形体高度h₃，包括以下步骤：

步骤601：建立以(m₀,n₀)为原点，水平向右为X轴正向、竖直向下为Y轴正向的局部坐标系；

步骤602：形成以直角坐标表示点的地形边界点序列链表Point_List：先按行从上往下、列从左往右遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，然后按行从下往上、列从右往左遍历数组A(i,j)，每行遇到第1个地形底面边界点标志值为-1的栅格点时，将该栅格点所属数组元素的下标i、j分别乘以DEM格网的分辨率cellsize即得该栅格点的x、y直角坐标，按顺序将边界点直角坐标依次存入边界点序列链表Point_List中，该链表存储的边界点连线形成一个多边形；

步骤603：令ang＝3°,定义二维数组Data[30][4]；

步骤604：对链表Point_List中的多边形各顶点做以点(x₀,y₀)为旋转中心、3°为转角的旋转变换，旋转后的多边形各顶点保存在链表Point_List中；

步骤605：求链表Point_List中的多边形的轴向包围盒的坐标x_max、x_min、y_max、y_min，计算此包围盒边长L1＝|x_max-x_min|、L2＝|y_max-y_min|，计算此包围盒的面积s0＝L1×L2；

步骤606：将ang、L1、L2、s0四个值作为一行依次存入数组Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；

步骤607：将步骤604～606循环29次；

步骤608：求数组Data中第4列包围盒面积的最小值，记录该面积最小值所在的数组行号为m′；

步骤609：读取数组Data中第m′行的前3列值，分别用an、l、w记录；令an＝90°-an，定义an为楔形体底面矩形方向角，它是长度为l的矩形边与正东方向的夹角，w为底面矩形的另一边长；

步骤610：楔形体高度为h₃＝H_ta-h_min。

2.根据权利要求1所述的基于分层高程云图的地形几何参数提取方法，其特征是：所示的步骤106中，根据H值确定地形类型参数k_d和高程分层高度Δh，包括以下步骤：

步骤201：判断H是否小于或等于20，若是，转至步骤202；若不是，转至步骤203；

步骤202：该地形判定为平原，令地形类型参数k_d＝1，高程分层高度Δh＝20；

步骤203：判断H是否小于或等于200，若是，转至步骤204；若不是，转至步骤205；

步骤204：该地形判定为丘陵，令地形类型参数k_d＝2，高程分层高度Δh＝30；

步骤205：判断H是否小于或等于500，若是，转至步骤206；若不是，转至步骤207；

步骤206：该地形判定为低山，令地形类型参数k_d＝3，高程分层高度Δh＝40；

步骤207：判断H是否小于或等于1000，若是，转至步骤208；若不是，转至步骤209；

步骤208：该地形判定为中山，令地形类型参数k_d＝4，高程分层高度Δh＝50；

步骤209：该地形判定为高山，令地形类型参数k_d＝5，高程分层高度Δh＝50。

3.根据权利要求1所述的基于分层高程云图的地形几何参数提取方法，其特征是：所示的步骤110中，根据h_min、h_max和Δh计算并形成分层高程云图，包括以下步骤：

步骤301：根据h_min、h_max和Δh确定地形高程分层数L_n，L_n＝[(h_max-h_min)/Δh]+1，符号[]表示取整数；

步骤302：将区间(h_min，h_max]等分为L_n个小区间并按值由小到大排列，若L_n＝2，则值较小的小区间分配的颜色为青色，值较大的小区间分配的颜色为红色；若L_n≥3，则值较小的小区间分配的颜色为青色，值最大的小区间分配的颜色为红色，中间位置的小区间按照绿、黄、橙色调依次与小区间对应分配颜色，每个小区间分配一个颜色；

步骤303：遍历数组A(i,j)，判断每个栅格点高程值属于哪个小区间，将该小区间分配的颜色值赋值给Point类的栅格点高程分层颜色变量，数组A(i,j)记录了所关注地形的分层高程云图。

4.根据权利要求1所述的基于分层高程云图的地形几何参数提取方法，其特征是：所示的步骤111中，对待研究的地形区域进行八邻域边界跟踪，并在数组A(i,j)中对地形底面边界点做标记，包括以下步骤：

步骤401：根据步骤302确定地形底面边界就是L_n个小区间中值最小的区间所对应的栅格点围成的封闭区域的外轮廓，地形底面边界栅格点的颜色为青色；

步骤402：按照从下往上、从左往右的顺序遍历数组A(i,j)，遇到的高程分层颜色为青色的第一个栅格点即为地形底面边界的起点，记该点为P_s，其Point类存储的栅格点行数、列数分别记为i_s、j_s，令该点Point类中的地形底面边界点的标志值为-1；

步骤403：设定八邻域搜索的方向编码，记搜索到的地形底面边界上的当前点为P_c，其栅格点行数、列数分别为i_c、j_c，将指向行数、列数分别为i_c-1、j_c-1的栅格点的搜索方向编码规定为0，将指向行数、列数分别为i_c-1、j_c的栅格点的搜索方向编码规定为1，将指向行数、列数分别为i_c-1、j_c+1的栅格点的搜索方向编码规定为2，将指向行数、列数分别为i_c、j_c+1的栅格点的搜索方向编码规定为3，将指向行数、列数分别为i_c+1、j_c+1的栅格点的搜索方向编码规定为4，将指向行数、列数分别为i_c+1、j_c的栅格点的搜索方向编码规定为5，将指向行数、列数分别为i_c+1、j_c-1的栅格点的搜索方向编码规定为6，将指向行数、列数分别为i_c、j_c-1的栅格点的搜索方向编码规定为7；

步骤404：令底面边界上搜索到的当前栅格点P_c＝P_s，搜索方向码变量S_d＝0；

步骤405：查找P_c沿S_d方向上的栅格点P_n的高程分层颜色值C_n；

步骤406：判断C_n是否为青色，若是，转至步骤407；若不是，转至步骤409；

步骤407：将数组A(i,j)中P_n位置的栅格点的地形底面边界点的标志值置为-1，令P_c＝P_n；

步骤408：判断P_n是否等于P_s，若不是，转至步骤412；若是，则地形底面边界跟踪结束；

步骤409：令S_d＝S_d+1；

步骤410：判断S_d是否大于或等于8，若是，转至步骤411；若不是，转至步骤405；

步骤411：令S_d＝S_d-8，转至步骤405；

步骤412：令S_d＝S_d-2；

步骤413：判断S_d是否小于0，若是，转至步骤414；若不是，转至步骤405；

步骤414：令S_d＝S_d+8，转至步骤405。

5.根据权利要求1所述的基于分层高程云图的地形几何参数提取方法，其特征是：所示的步骤118中，求地形底面边界围成区域的平均坡度sl_a，包括以下步骤：

步骤501：按照先行后列顺序遍历数组A(i,j)，对于地形底面边界点标志值为-1的栅格点以及一行内处在两个边界点之间的栅格点，逐个计算这些栅格点的坡度值sl：设待计算的栅格点高程为h_e，该栅格点左上方的栅格点高程为h_a，该栅格点正上方的栅格点高程为h_b，该栅格点右上方的栅格点高程为h_c，该栅格点正左方的栅格点高程为h_d，该栅格点正右方的栅格点高程为h_f，该栅格点左下方的栅格点高程为h_g，该栅格点正下方的栅格点高程为h_h，该栅格点右下方的栅格点高程为h_i，令f_x＝(h_i-h_g+2(h_f-h_d)+h_c-h_a)/(8×cellsize)，f_y＝(h_a-h_g+2(h_b-h_h)+h_c-h_i)/(8×cellsize)，cellsize表示DEM格网的分辨率，则该栅格点的坡度

步骤502：对边界区域内的每个栅格点，将其Point类中坡度变量值置为各自的计算值sl；

步骤503：遍历数组A(i,j)，计算地形底面边界区域内所有栅格点坡度值之和sl_sum，则地形底面边界围成区域的平均坡度sl_a＝sl_sum/G_n。

Images