基于DEM的数字地形分析有哪些应用

第9章 DEM与数字地形分析

数字地面模型于1958年提出，特别是基于DEM的GIS空间分析方法的出现，使传统的地形分析方法产生了革命性的变化，数字地形分析方法逐步形成和完善。目前，基于DEM的数字地形分析已经成为GIS空间分析中最具特色的部分，在测绘、遥感及资源调查、环境保护、城市规划、灾害防治及地学研究各方面发挥越来越重要的作用。本章首先介绍了数字高程模型的基本概念和建立步骤，然后从基本坡面因子、特征地形因子、水文因子和可视域等方面简述数字地形分析的主要内容和研究方法。

9.1 基本概念

9.1.1 数字高程模型

数字高程模型（Digital Elevation

Model，简称DEM）是通过有限的地形高程数据实现对地形曲面的数字化模拟（即地形表面形态的数字化表示），它是对二维地理空间上具有连续变化特征地理现象的模型化表达和过程模拟。由于高程数据常常采用绝对高程（即从大地水准面起算的高度），DEM也常常称为DTM（Digital

Terrain

Model）。“Terrain”一词的含义比较广泛，不同专业背景对“Terrain”的理解也不一样，因此DTM趋向于表达比DEM更为广泛的内容。

从研究对象与应用范畴角度出发，DEM可以归纳为狭义和广义两种定义。从狭义角度定义，DEM是区域表面海拔高程的数字化表达。这种定义将描述的范畴集中地限制在“地表”、“海拔高程”及“数字化表达”内，观念较为明确。从广义角度定义，DEM是地理空间中地理对象表面海拔高度的数字化表达。这是随着DEM的应用不断向海底、地下岩层以及某些不可见的地理现象（如空中的等气压面等）延伸，而提出的更广义的概念。该定义将描述对象不再限定在“地表面”，因而具有更大的包容性，有海底DEM、下伏岩层DEM、大气等压面DEM等。

数学意义上的数字高程模型是定义在二维空间上的连续函数H?f(x,y)。由于连续函数的无限性，DEM通常是将有限的采样点用某种规则连接成一系列的曲面或平面片来逼近原始曲面，因此DEM的数学定义为区域D的采样点或内插点Pj按某种规则?连接成的面片M的集合：

（9.1）

DEM按照其结构，可分为规则格网DEM、TIN、基于点的DEM和基于等高线的DEM等。由于规则格网结构简单，算法设计明了，在实际运用中被广泛采用。本书中的DEM仅指规则格网DEM。

DEM?{Mi??(Pj)Pj(xj,yj,Hj)?D,j?1,?n,i?1,?,m}

9.1.2 数字地形分析

数字地形分析（Digital Terrain Analysis,

DTA），是指在数字高程模型上进行地形属性计算和特征提取的数字信息处理技术。DTA技术是各种与地形因素相关空间模拟技术的基础。

地形属性根据地形要素的关系特征和计算特征，可以归纳为地形曲面参数（parameters）、地形形态特征（features）、地形统计特征（statistics）和复合地形属性（compound

attributes）。

地形曲面参数具有明确的数学表达式和物理定义，并可在DEM上直接量算，如坡度、坡向、曲率等。地形形态特征是地表形态和特征的定性表达，可以在DEM上直接提取，其特点是定义明确，但边界条件有一定的模糊性，难以用数学表达式表达，如在实际的流域单元的划分中，往往难于确定流域的边界。地形统计特征是指给定地表区域的统计学上的特征。复合地形属性是在地形曲面参数和地形形态特征的基础上，利用应用学科（如水文学、地貌学和土壤学）的应用模型而建立的环境变量，通常以指数形式表达。

数字地形分析的主要内容有两方面，一是在复杂的现实世界地理过程中各影响因子和简单、高效、精确、易于理解的抽象与计算机实现中找到平衡。简单地说，就是提取描述地形属性和特征的因子，并利用各种相关技术分析解释地貌形态、划分地貌形态等。二是DTM的可视化分析。数字地形分析中可视化分析的重点在于地形特征的可视化表达和信息增强，以帮助传达地形曲面参数、地表形态特征和复合地形属性的信息。

根据分析内容，常用的数字地形分析的方法有以下几种（图9.1）：

1.提取坡面地形因子

地形定量因子是为有效地研究与表达地貌形态特征所设定的具有一定意义的参数或指标。从地形地貌的角度考虑，地表是由不同的坡面组成的，而地貌的变化，完全源于坡面的变化。常用的坡面地形因子有坡度、坡向、平面曲率、坡面曲率、地形起伏度、粗糙度、切割深度等。

2.提取特征地形要素

（1）流域分析

流域分析主要是根据地表物质运动的特性，特别是水流运动的特点，利用水流模拟的方法来提取水系、山脊线、谷底线等地形特征线，并通过线状信息分析其面域特征。

（2）可视域分析

可视性分析包括两方面内容，一个是两点之间的通视性(Intervisibility)，另一个是可视域(ViewShed)，即对于给定的观察点所覆盖的区域。

3.地形统计特征分析

地形统计分析是应用统计方法对描述地形特征的各种可量化的因子或参数进行相关、回归、趋势面、聚类等统计分析，找出各因子或参数的变化规律和内在联系，并选择合适的因

子或参数建立地学模型，从更深层次探讨地形演化及其空间变异规律。

提取坡度

提取坡面曲率

流域分析

地形统计特征分析 可视域分析

图

9.1

数字地形分析常用方法

9.2 DEM建立

9.2.1 DEM建立的一般步骤

数字高程模型的建立过程是一个模型建立过程。从模型论角度讲，就是将源域（地形）表现在另一个域（目标域或DEM）中的一种结构，建模的目的是对复杂的客体进行简化和抽象，并把对客体（源域，DEM中为地形起伏）的研究转移到对模型的研究上来。

模型建立之初，首先要为模型构造一个合适的空间结构（spatial

framework）。空间结构是为把特定区域内的空间目标镶嵌在一起而对区域进行的划分，划分出的各个空间范围称为位置区域或空间域。空间结构一般是规则的（如格网），或不规则的（如不规则三角网TIN）。

建立在空间结构基础上的模型是由n个空间域的有限集合组成。由于空间数据包含位置特征和属性特征，而属性特征是定义在位置特征上的，因此每一个空间域就是由空间结构到属性域的计算函数或域函数。模型的可计算性要求有两点，一是空间域的数量、属性域和空间结构是有限的，二是域函数是可计算的。构筑模型的一般内容和过程为：

①采用合适的空间模型构造空间结构；

②采用合适的属性域函数；

③在空间结构中进行采样，构造空间域函数；

④利用空间域函数进行分析。

当空间结构为欧几里德平面，属性域是实数集合时，模型为一自然表面。将欧几里德平面充当水平的XY平面，属性域给出Z坐标（或高程），模型即为数字高程模型。

对于数字高程模型而言，空间结构的构造过程即为DEM的格网化过程（形成格网），属性值为高程，构造空间域函数即为内插函数的确定，利用空间域函数进行分析就是求取格网点的函数值。

9.2.2 规则格网DEM的建立

DEM是在二维空间上对三维地形表面的描述。构建DEM的整体思路是首先在二维平面上对研究区域进行格网划分（格网大小取决于DEM的应用目的），形成覆盖整个区域的格网空间结构，然后利用分布在格网点周围的地形采样点内插计算格网点的高程值，最后按一定的格式输出，形成该地区的格网DEM（图9.1）。

图9.2 格网DEM建立流程

9.2.3 DEM内插方法

DEM建立过程中的关键环节是根据采样点的值内插计算格网点上的高程值。内插是指根据分布在内插点周围的已知参考点的高程值求出未知点的高程值，它是DEM的核心问题，贯穿于DEM的生产、质量控制、精度评定、分析应用的各个环节。

随着DEM的发展和完善，已经提出了多种高程内插方法。根据不同的分类标准，有不同的内插方法分类，例如按数据分布规律分类，有基于规则分布数据的内插方法、基于不规则分布的内插方法和适合于等高线数据的内插方法等；按内插点的分布范围，内插方法分为整体内插、局部内插和逐点内插法；从内插函数与参考点的关系方面，又分为曲面通过所有采样点的纯二维插值方法和曲面不通过参考点的曲面拟合插值方法；从内插曲面的数学性质来讲，有多项式内插、样条内插、最小二乘配置内插等内插函数；从对地形曲面理解的角度，内插方法有克立金法、

多层曲面叠加法、加权平均法、分形内插等。表9.1对各种DEM内插分类方法进行了简要的总结和归纳。

本小节仅从内插点的分布范围来看，简要介绍整体内插法、局部内插法和逐点内插法。详细介绍参见第十章。

表9.1 DEM内插分类方法

整体内插是指在整个区域用一个数学函数来表达地形曲面。整体内插函数通常是高次多项式，要求地形采样点的个数大于或等于多项式的系数数目。整体内插方法有整个区域上函数的唯一性、能得到全局光滑连续的DEM、充分反映宏观地形特征等优点。但由于整体内插函数往往是高次多项式，它也有保凸性较差、不容易得到稳定的数值解、多项式系数的物理意义不明显、解算速度慢且对计算机容量要求较高、不能提供内插区域的局部地形特征等缺点。在DEM内插中，一般是与局部内插方法配合使用，例如在使用局部内插方法前，利用整体内插去掉不符合总体趋势的宏观地物特征。另外也可用来进行地形采样数据中的粗差检测。

局部分块内插是将地形区域按一定的方法进行分块，对每一分块，根据其地形曲面特征单独进行曲面拟合和高程内插。一般按地形结构线或规则区域进行分块，分块的大小取决于地形的复杂程度、地形采样点的密度和分布。为保证相邻分块之间的曲面平滑连接，相邻分块之间要有一定宽度的重叠，或者对内插曲面补充一定的连续性条件。这种方法简化了地形的曲面形态，使得每一分块可用不同的曲面表达，同时得到光滑连续的空间曲面。不同的分块单元可以使用不同的内插函数。常用的内插函数有线性内插、双线性内插、多项式内插、样条函数、多层曲面叠加法等。

逐点内插是以内插点为中心，确定一个邻域范围，用落在邻域范围内的采样点计算内插点的高程值。逐点内插本质上是局部内插，但与局部分块内插不同的是，局部内插中的分块范围一经确定，在整个内插过程中其大小、形状和位置是不变的，凡是落在该块中的内插点，都用该块中的内插函数进行计算，而逐点内插法的邻域范围大小、形状、位置乃至采样点个数随内插点的位置而变动，一套数据只用来进行一个内插点的计算。

逐点内插法要注意两个问题，一是选择合适的内插函数，内插函数决定着DEM精度、DEM连续性、内插点邻域的最小采样点个数和内插计算效率。二是确定内插点邻域，内插点的邻域大小和形状、邻域内参加内插计算的数据点的个数、采样点的权重、采样点的分布、附加信息等不仅会影响到DEM的内插精度，也影响到内插速度。逐点内插方法计算简单，内插效率较高，应用比较灵活，是目前较为常用的一类DEM内插方法。

在建立DEM时，要根据情况选择合适的、运算效率高的方法。而众多内插方法并不是独立的，而往往是相互结合使用，这在后续的章节里会讲到。

9.3 数字地形分析

地形分析是地形环境认知的一种重要手段，传统的地形分析是基于二维平面地图进行的。从基于纸质地图的地形分析发展到到基于数字地图的地形分析，计算机取代了大量的人工计算和绘制，地形分析的手段、功能发生了一次飞跃；可视化技术和虚拟现实技术的发展，使得建立三维实时、交互的仿真地形环境成为可能，同时也需要实现三维地形环境中的地形分析。特别是DEM的出现和大量应用，使得从地形属性中提取各类地形参数和特征因子更加的简便和准确。

用来描述地形特征和空间分布的地形参数很多，不同的应用目的，不同的学科和领域对此的理解和分类也不同。本章将综合相关知识，着重介绍基本因子分析、地形特征提取、水文分析和可视域分析。

9.3.1 基本因子分析

本质上讲，DEM是地形的一个数学模型，可以看成是一个或多个函数的集合。实际上许多地形因子就是从这些函数进行一阶或二阶推导出来的，也有的通过某种组合或复合运算得到。基本地形因子包括斜坡因子（坡度、坡向、坡度变化率、坡向变化率等）、面积因子（表面积、投影面积、剖面积）、体积因子（山体体积、挖填体积）和面元因子（相对高差、粗糙度、凹凸系数、高程变异等）。

本节将阐述一些常用的基本地形因子，为了方便起见，并从实际应用角度考虑，本节这些地形因子的计算都是基于格网DEM。 1. 坡度

严格地讲，地表面任一点的坡度是指过该点的切平面与水平地面的夹角。坡度表示了地表面在该点的倾斜程度，在数值上等于过该点的地表示），即：

(N)

?

n微分单元的法矢量与z轴的夹角（如图9.5所

??

z?n

ArcCos(z?n

Slope =

图9.5 地表单元坡度示意图

（9.1）

当具体进行坡度提取时，常采用简化的差分公式，完整的数学表示为：

Slope?arctan

（9.2）

式中，fx是X方向高程变化率，fy是Y方向高程变化率。

地面坡度实质是一个微分的概念，地面上每一点都有坡度，它是一个微分点上的概念，是地表曲面函数z =

f(x,y)在东西、南北方向上的高程变化率的函数。实际应用中，坡度有两种表示方式（如图9.6）：

? 坡度(degree of slope)：即水平面与地形面之间夹角。 ? 坡度百分比（percent

slope）：即高程增量（rise）与水平增量（run）之比的百分数。

图9.6 坡度的两种表示方法

拟合曲面法是解求坡度的最常用的方法。拟合曲面法，一般采用二次曲面，即在3×3的DEM

栅格分析窗口中（如图9.7）进行，每个栅格中心为一个高程值，分析窗口在DEM数据矩阵中连续移动完成整个区域的计算工作。常用的计算fx、

fy的方法是三阶反距离平方权，该算法也用于ArcView和ARC/INFO。其计算方法为：

zi?1,j?1?2zi,j?1?zi?1,j?1?zi?1,j?1?2zi,j?1?zi?1,j?1?fx??8g?

?

zi?1,j?1?2zi?1,j?zi?1,j?1?zi?1,j?1?2zi?1,j?zi?1,j?1?fy?

?8g? （9.3）

式（9.3），g为格网间距。 2. 坡向

??n

坡向定义为：地表面上一点的切平面的法线矢量n在水平面的投影xoy与过该点的正北

Aspect?arctg(

)

fx （9.4）

y

方向的夹角（如表7.1中的坡向示意图所示，x轴为正北方向）。其数学表达公式为： f

对于地面任何一点来说，坡向表征了该点高程值改变量的最大变化方向。在输出的坡向数据中，坡向值有如下规定：正北方向为0°，顺时针方向计算，取值范围为0°～360°。

坡向可在DEM数据中用式9.4直接提取。但应注意，由于式9.4求出坡向有与x轴正向和x轴负向夹角之分，此时就要根据fx和fy的符号来进一步确定坡向值（如表9.2所示）

注：上述情况假定所建立的DEM数据从南向北获取的，且x轴与正北方向重合，否则上述公式求得的坡向值，还应加上x轴偏离正北方向的夹角值。

采用这种方法求取的坡向分级比较详细，但实际应用中往往需要给予归并，在ArcView和ArcGIS软件中，通常把坡向综合成九种坡向：平缓坡（－1）、北坡（0°

- 22.5°, 337.5° - 360°）、东北坡（22.5° - 67.5°）、东坡（67.5° - 112.5°）、东南坡（112.5° -

157.5°）、南坡（157.5° - 202.5°）、西南坡（202.5° - 247.5°）、西坡（247.5° - 292.5°）、西北坡（292.5°

- 337.5°）。

图9.8 原始DEM数据及实验区等高线图

图9.9 ARCVIEW软件下提取的坡度图

图9.10 由DEM提取的坡向图

3. 曲率

曲率是对地形表面一点扭曲变化程度的定量化度量因子，地面曲率在垂直和水平两个方向上分量分别称为平面曲率和剖面曲率。地形表面曲率反映了地形结构和形态，同时也影响着土壤有机物含量的分布，在地表过程模拟、水文、土壤等领域有着重要的应用价值和意义。

剖面曲率是对地面坡度的沿最大坡降方向地面高程变化率的度量。数学表达式为：

）

y（E）

图9.11平面曲率示意

Kv??

p2r?2pqs?q2t

(p2?q2)?p2?q2 （9.5）

平面曲率指在地形表面上，具体到任何一点P，指用过该点的水平面沿水平方向切地形表面所得的曲线在该点的曲率值（图9.11所示）。平面曲率描述的是地表曲面沿水平方向的弯曲、变化情况，也就是该点所在的地面等高线的弯曲程度。从另一个角度讲，地形表面上一点的平面曲率也是对该点微小范围内坡向变化程度的度量。数学表达式为：

Kh??

q2r?2pqs?p2t

(p2?q2)?p2?q2 （9.6）

曲率数学表达式中，利用离散的DEM数据把地表曲面数学模拟为一个连续的曲面H(x,y)，x和y地面点的平面坐标值，H(x,y)为地面点高程值，式中其它符号所表示的意义为：

?H

?x，是x方向高程变化率； ?Hq?

?y，是y方向高程变化率； p?

r?

?2H

?x2，对高程值在x方向上的变化率进行同方向求算变化率，即x方向高程变化率

的变化率；

?2Hs?

?x?y，对高程值在x方向上的变化率进行y方向上求算变化率，即x方向高程变化

率在y方向的变化率；

t??2H

?y2，对高程值在y方向上的变化率同方向上求算变化率，即y方向高程变化率的

变化率。

曲率因子的提取算法的基本原理为：在DEM数据的基础上，根据其离散的高程数值，把地表模拟成一个连续的曲面，从微分几何的思想出发，模拟曲面上每一点所处的垂直于和平行于水平面的曲线，利用曲线曲率的求算方法的推导得出各个曲率因子的计算公式。利用公式求算出每一点的曲率值的关键在于确定得出式中各个参量的值，在DEM中求算高程的微分分量有一套独特的算法，最常用是三阶反距离平方权差分。对每一个栅格点都确定一个3×3的分析窗口，其过程如图9.12所示。利用ArcView所提取的剖面曲率与平面曲率图如图9.13和9.14所示。

(ay?2by?cy)?(gy?2hy??

8*Cellsize

"""""

图9.12 地面曲率提取步骤流程图

4. 宏观地形因子

地形起伏度、地形表面粗糙度与地表切割深度等地形因子是描述和反映地形表面较大区

域内地形的宏观特征，在较小的区域内并不具备任何地理和应用意义。这些参数对于在宏观尺度上的水土保持、土壤侵蚀特征、地表发育、地貌分类等研究中具有重要的理论意义。基于栅格DEM计算宏观地形因子时，关键在于确定分析半径的大小。不同地貌类型、不同分辨率的数据，计算宏观地形因子所取的分析半径大小是不一。因此，确定一个合适的分析窗口半径或分析区域，使得求取的宏观因子能够准确反映地面的起伏状况与水土流失特征，是提取算法的核心步骤和决定信息提取效果与有效性的关键。

⑴ 地线起伏度 地形起伏度是指，在所指定的分析区域内所有栅格中最大高程与最小高程的差。可表示为如下公式：

(9.7)

式中，RFi指分析区域内的地面起伏度，Hmax指分析窗口内的最大高程值，Hmin指分析窗口内的最小高程值。

地形的起伏是反映地形起伏的宏观地形因子，在区域性研究中，利用DEM数据提取地形起伏度能够直观的反映地形起伏特征。在水土流失研究中，地形起伏度指标能够反映水土流失类型区的土壤侵蚀特征，比较适合区域水土流失评价的地形指标。

⑵ 地形粗糙度

地表粗糙度，一般定义为地表单元的曲面面积S曲面与其在水平面上的投影面积S水平之比。用数学公式表达为：

R = S曲面 / S水平 (9.8)

(9.10)

式中，Di指地面每一点的地表切割深度，Hmean指一个固定分析窗口内的平均高程，Hmin

指一个固定分析窗口内的最低高程。

地表切割深度直观的反映了地表被侵蚀切割的情况，并对这一地学现象进行了量化，是研究水土流失及地表侵蚀发育状况时的重要参考指标，其提取算法可参照地表起伏度的提取。

9.3.2 地形特征分析

虽然地表形态各式各样，但地形点、地形线、地形面等地形结构的基本特征构成了地形的骨架，因此一般的地形特征提取主要是指地形特征点、线、面的提取，并进而通过基本要素的组合进行地表形态分析。特征地形要素的提取更多地应用较为复杂的技术方法，其中山谷线、山脊线的提取采用了全域分析法，成为数字高程模型地学分析中很具特色的数据处理内容。

1. 地形特征点提取

地形特征点主要包括山顶点（peak）、凹陷点（pit）、脊点（ridge）、谷点（channel）、鞍点（pass），平地点（plane）等。利用DEM提取地形特征点，可通过一个3×3或更大的栅格窗口，通过中心格网点与8个邻域格网点的高程关系来进行判断会获取。即在一个局部区域内，用x方向和y方向上关于高程z的二阶导数的正负组合关系来判断（见表9.3）。该方法假设DEM表面为z

= f(x,y),但由于真实地表与数学表面的差别，在利用该方法在DEM上提取特征点，结果常产生伪特征点。

表9.3

中的关于地形特征点的判断是在局部区域内利用x，y方向的凹凸性判断的，该

?2z

2

判断法十分适合利用在DEM上判断地形特征点。在DEM中可以利用差分的方法得到?x和

?2z

?x2的值。

除上述算法外，在一个3×3的栅格窗口中，也可以直接利用中心格网点与8个邻域格网点的高程关系来进行判断地形特征点。具体方法为：

假设有一个如图9.10所示的3×3

窗口。则：如果（Zi，j-1 - Zi

，j

）（Zi，j-1 - Zi，j）>0

（1）当Zi，j+1> Zi，j 则VR（i，j）= -1 （2）当Zi，j+1< Zi，j 则VR（i，j）= 1

(i,j+1)

图9.15 差分算法示意图

如果（Zi-1，j - Zi，j）（Zi+1，j - Zi，j）（3）当Zi+1> Zi，j 则VR（i，j）= -1 （4）当Zi+1<

Zi，j 则VR（i，j）= -1

如果（1）和（4）或（2）和（3）同时成立，则VR（i，j）= 2 如果以上条件都不成立，则VR（i，j）= 0

??1，表示谷点?

?1，表示脊点

VR?i，j???

?2，表示鞍点?0，表示其他点?其中，

2. 山脊线和山谷线提取

山脊线和山谷线构成了地形起伏变化的分界线（骨架线），因此它对于地形地貌研究具

有重要的意义。另一方面，对于水文物理过程研究而言，由于山脊、山谷分别表示分水性与汇水性，山脊线和山谷线的提取实质上也是分水线与汇水线的提取。这一特性又使得山脊线和山谷线在许多工程应用方面有着特殊的意义。

在对山脊线、山谷线的提取方法中，基于规则格网DEM的方法是主要。从原理上来分，主要分为以下四种： （1） 基于图像处理技术的原理

因为规则格网DEM数据事实上是一种栅格形式的数据，可以利用数字图像处理中的技术来设计算法。利用数字图像处理技术设计的算法大都采用各种滤波算子进行边缘提取。基于该原理有一种简单移动窗口的算法，其主要思路是：

① 计一个2×2窗口以对DEM格网阵列进行扫描；

② 第一次扫描中，将窗口中的具有最低高程值的点进行标记，自始至终未被标记的点即为山脊线上的点；

③ 第二次扫描中，将窗口中的具有最高高程值的点进行标记，自始至终未被标记的点

行，交给我搞定。