基于DEM地形要素提取方法的实验

孙嘉骏,黄天勇,陈 科

（1.中国水电顾问集团 昆明勘测设计研究院有限公司 测绘地理信息分院,云南 昆明 650041）

DEM是通过有限的地形高程数据实现对地形曲面的数字化模拟或是地形表面形态的数字化表示的量化模型[1]。近年来，随着空间数据基础设施的建设和“数字地球”战略的实施，加快了DEM与地理信息系统、遥感等的一体化进程，为DEM的应用开辟了更广阔的天地。

1 DEM概述

1.1 DEM的含义

DEM表示三维向量有限序列，是对地球表面地形地貌的一种离散的数字表达，也是测绘学中最常用的地形表示模型，用函数形式描述为：iiii

式中，（Xi,Yi）是平面坐标；Zi是（Xi,Yi）对应的高程。当该序列中各平面向量的平面位置呈规则格网排列时，其平面坐标可省略，此时DEM就简化为一维向量序列{Zi, i=1,2,3,…,n}[1]。

1.2 DEM的表示模型及特点

随着对DEM需求的不断增加，对DEM数据的多样性、实用性及精度提出了更高的要求。同时由于其数据来源不同，导致数据存储结构的差别。DEM数据结构表示模型主要包括：

1）不规则离散点型。由不规则离散点表示的DEM是将连续地球表面形态离散成在某一个区域D上的以Xi、Yi、Zi三维坐标形式存储的高程点Zi((Xi,Yi)ID)的集合[2]。离散点DEM往往是通过测量直接获取的地球表面的原始特征点数据。这些特征点之间相互独立，彼此间没有任何联系，它是DEM中最简单的数据组织形式。

2）规则格网型。规则格网DEM是指在水平方向和垂直方向按等间隔方式记录格网点的坐标，格网点的高程Z表示地形，格网点的平面坐标隐含在行列号中，故适宜用矩阵形式进行存储，即按行(或列)逐一记录每一个格网单元的高程值。通常，为了压缩栅格DEM的冗余数据，可采取用游程编码或四叉树编码方法对数据进行处理；为准确表示地形细部，可采用附加地形特征数据描述地形结构[3]。

3）等高线型。等高线通常被存成一个有序的坐标点对序列，是地形表达最为常用的形式，能较为直观科学地反映地面高程、山体、坡度、坡形、山脉走向等基本的地貌形态及其变化。利用等高线重新构建地表模型时，只能近似表示真实的地球表面，在水平切割过程中丢失的详细地表信息是不可能从等高线中恢复的。同时，单条等高线无法直接反映地貌形态，必须通过等高线组间接地表示地貌形态。

4）不规则三角网（TIN）型。为克服规则格网的缺点，采用附加地形特征数据的方法，按一定的规则将离散点连接成覆盖整个区域且互不重叠、结构最佳的三角形，从而构成完整的DEM，这就是不规则三角网数字高程模型[4]。TIN克服了高程矩阵中冗余数据的问题，其最主要的优点就是具有可变的分辨率，可根据不同地形，选取合适的采样点数。另外，TIN还具有考虑重要表面数据点的能力，能利用地貌的特征点线，较好地表示复杂地形。但其缺点是数据结构复杂，构成TIN算法复杂，计算时间长，大区域的TIN的管理十分困难。

DEM具有表达样式多、精度恒定、自动化、实时性强、尺度综合性好等特点，在实际研究与应用中，应根据区域的研究目的和大小，对数据存储的结构格式进行转换，以达到最优化的存储和快速查询与编辑。

2 基于DEM的地形要素提取实验

2.1 实验区概况

武夷山地处中国福建省的西北部，江西省东部，位于福建与江西的交界处。地理坐标为：北纬27°32'36"～27°55'15"，东经 117°24'12"～118°02'50"，总面积999.75 km2，核心面积63 575 hm2，同时划定了外围保护缓冲区地带，面积27 888 hm2。武夷山风景名胜区主要景区方圆70 km2，平均海拔350 m，属典型的丹霞地貌，素有“碧水丹山”、“奇秀甲东南”之美誉。

2.2 DEM应用分析

本实验使用软件为美国ESRI公司的ArcGIS9.X。首先对实验区元数据进行浏览，得知该数据覆盖范围，投影为横轴墨卡托投影，单位为m，其中央经线为东经117°；再查看图层属性信息，在Source选项卡中，可看到该栅格图层的相关信息，如该图层有行列值，栅格大小约为26.17 m。栅格图层的行列数和栅格大小是2个很重要的信息，对后面的参数设定有最直接的影响。

2.2.1 地形指标的提取

地形指标是描述和评价一个地区地形情况的有效工具，是进行多要素叠加分析中一个可选条件，也是最短路径分析中获取成本图层的一个可行方法。地形指标一般包括坡度、坡度变率、坡向、坡向变率、起伏度和粗糙度等。

通过ArcMap的“Aspect”工具可先提取出原DEM的坡向信息（见图1）。而提取坡向变率（SOA）的过程较为复杂，首先要进行反地形的提取，然后对反地形（见图2）提取坡向信息，再分别对正、反地形的坡向进行坡度提取（所得结果为“SOA1”和“SOA2”），最后利用Raster Calculator输入模型表达式：SOA=(([SOA1]+ [SOA2]) - Abs([SOA1]- [SOA2]))/ 2，从而得到误差纠正后的“SOA”。这里需要注意的是，在提取反地形坡向变率时，需用一个常数值减去原DEM，从而得到反地形数据。由于求算坡向数据与反地形的绝对高程无关，仅与相对高程有关，因此这里的常数可任意设置（一般取原DEM的最大高程值）。

图1 实验区坡向图层

图2 实验区反地形图层

2.2.2 等高线的提取

一般来说，在提取等高线时等高距设置为大于原DEM栅格大小减30%的区间为好，这样可有效避免出现孤立的等高线小尖角、等高线小圆圈等情况。例如，在本实验中原DEM的栅格大小约为26.17 m，将等高距设置为大于20 m的数值，得到的结果比较符合实际情况。这也从另一个角度说明了查看原数据信息的重要性（见图3、图4）。

图3 提取后查看100 m等高距

2.2.3 山脊线、山谷线的提取

山脊线和山谷线是地形的重要特征。从地貌学的角度，它们是坡向变率较大的点的连线；从水文学的角度，它们则分别是分水线和汇水线[5]。相应地，提取这2种地形特征的方法就有空间分析和水文分析2种。

1）空间分析法。先用邻域统计分析工具，求算出邻域平均值图层；然后用原DEM减去邻域平均值图层，得到差值图层[c]；再通过栅格计算器，分别求出“sj1=[cha]> 0 & [SOA]> 70，sj2=[cha]> 0 & [SOA]>60”，以及sj3=[cha]> 20 & [SOA]> 60的逻辑结果。

比较“sj1”、“sj2”和“sj3”3个图层，可知设定不同的阈值，所得到的结果也是不同的。一般来说，阈值范围越大，山脊线就越粗越密集，伪山脊线存在的可能性也越大；反之，则山脊线就越细越稀疏，伪山脊线存在的可能性也就越小。根据需要和具体情况对阈值进行设定，从而提取适当密度的山脊线和山谷线（见图5、图6）。

图5 提取sj1计算结果

图6 提取sj2计算结果

2）水文分析法。先对原DEM进行流向的提取，并利用ArcToolbox中水文分析Sink工具，对该流向数据进行洼地的检查。如存在洼地，则用Fill工具进行洼地填充，以消除错误的流向。对洼地填充后的结果进行流向提取，并利用流向数据进行汇流累积量提取，运用Raster Calculator提取汇流累积量为0的区域，并对其进行邻域平均值平滑。然后，对平滑后的结果进行重分类。经反复实验，0.555 4为最佳分界阈值，接近1的部分赋值为1，接近0的部分赋值为0。最后，计算正地形和前面重分类结果的乘积，将所得结果为0的部分重分类为Nodata,1的部分保留为1，从而完成了水文分析法的山脊线提取。山谷线的提取方法与山脊线基本类似，只是用反地形数据当作原数据，重复山脊线提取步骤，则可得到反地形的山脊线，即原地形的山谷线数据（见图7、图8）。

图7 重分类后提取山脊线图层

图8 重分类后提取山谷线图层

3 结 语

1）对于相似的要素提取过程或应用模型，应当掌握其共同的算法思想[6,7]，同时注意区别其不同之处。在具体应用中，需清楚其原理，了解更改每一步的参数设置之后会得到哪些相应的变化，对所得到的结果应结合已有的数据进行正误的判断和分析。

2）对于数据量较大的工程，在计算机设备有限的条件下，应当灵活设置邻域窗口、输出栅格尺寸等参数，在运算时间和结果精度之间寻找平衡点。

3）对于水文分析，应先对原DEM的流向数据作“Sink”检验，发现有洼地存在时，应及时填充，以消除不合理的河流方向。

4）阈值的设置是实验结果好坏的关键，有些可借鉴其他研究者的经验数值，有些则需实验比较，从而得出最优的阈值。

[1]李志林,朱庆.数字高程模型[M]. 武汉：武汉大学出版社,2008

[2]汤国安,刘学军,闾国年.数字高程模型及地学分析的原理与方法[M].北京:科学出版社,2005

[3]汤国安,杨昕.ArcGIS地理信息系统空间分析实验教程[M].北京:科学出版社,2008

[4]艾廷华,祝国瑞,张根寿.基于Delaunay三角网模型的等高线地形特征提取及谷地树结构化组织[J].遥感学报,2004,31(2):43-47

[5]靳海亮,康建荣,高井祥.利用等高线数据提取山脊(谷)线算法研究[J].武汉大学学报：信息科学版,2005(9):810-812

[6]朱庆,赵杰,钟正.基于规则格网DEM的地形特征线的提取算法[J].测绘学报.2004,33(1):77-82

[7]刘学军.基于规则网数字高程模型解译算法误差分析与评价[D].武汉：武汉大学，2002