GIS课件第9章 DEM与数字地形分析.doc

第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等。数学意义上的数字高程模型是定义在二维空间上的连续函数。由于连续函数的无限性，DEM通常是将有限的采样点用某种规则连接成一系列的曲面或平面片来逼近原始曲面，因此DEM的数学定义为区域D的采样点或内插点Pj按某种规则连接成的面片M的集合： （9.1）DEM按照其结构，可分为规则格网DEM、TIN、基于点的DEM和基于等高线的DEM等。由于规则格网结构简单，算法设计明了，在实际运用中被广泛采用。本书中的DEM仅指规则格网DEM。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）。x不规则分布点规则分布等高线分布 Y对每一格网点求取格网点高程图9.2 格网DEM建立流程9.2.3 DEM内插方法DEM建立过程中的关键环节是根据采样点的值内插计算格网点上的高程值。内插是指根据分布在内插点周围的已知参考点的高程值求出未知点的高程值，它是DEM的核心问题，贯穿于DEM的生产、质量控制、精度评定、分析应用的各个环节。随着DEM的发展和完善，已经提出了多种高程内插方法。根据不同的分类标准，有不同的内插方法分类，例如按数据分布规律分类，有基于规则分布数据的内插方法、基于不规则分布的内插方法和适合于等高线数据的内插方法等；按内插点的分布范围，内插方法分为整体内插、局部内插和逐点内插法；从内插函数与参考点的关系方面，又分为曲面通过所有采样点的纯二维插值方法和曲面不通过参考点的曲面拟合插值方法；从内插曲面的数学性质来讲，有多项式内插、样条内插、最小二乘配置内插等内插函数；从对地形曲面理解的角度，内插方法有克立金法、 多层曲面叠加法、加权平均法、分形内插等。表9.1对各种DEM内插分类方法进行了简要的总结和归纳。本小节仅从内插点的分布范围来看，简要介绍整体内插法、局部内插法和逐点内插法。详细介绍参见第十章。表9.1 DEM内插分类方法DEM内插数据分布规则分布内插方法不规则分布内插方法等高线数据内插方法内插范围整体内插方法局部内插方法逐点内插方法内插曲面与参考点关系纯二维内插曲面拟合内插内插函数性质多项式内插线性插值双线性插值高次多项式插值样条内插有限元内插最小二乘配置内插地形特征理解克立金内插多层曲面叠加内插加权平均值内插分形内插傅立叶级数内插整体内插是指在整个区域用一个数学函数来表达地形曲面。整体内插函数通常是高次多项式，要求地形采样点的个数大于或等于多项式的系数数目。整体内插方法有整个区域上函数的唯一性、能得到全局光滑连续的DEM、充分反映宏观地形特征等优点。但由于整体内插函数往往是高次多项式，它也有保凸性较差、不容易得到稳定的数值解、多项式系数的物理意义不明显、解算速度慢且对计算机容量要求较高、不能提供内插区域的局部地形特征等缺点。在DEM内插中，一般是与局部内插方法配合使用，例如在使用局部内插方法前，利用整体内插去掉不符合总体趋势的宏观地物特征。另外也可用来进行地形采样数据中的粗差检测。局部分块内插是将地形区域按一定的方法进行分块，对每一分块，根据其地形曲面特征单独进行曲面拟合和高程内插。一般按地形结构线或规则区域进行分块，分块的大小取决于地形的复杂程度、地形采样点的密度和分布。为保证相邻分块之间的曲面平滑连接，相邻分块之间要有一定宽度的重叠，或者对内插曲面补充一定的连续性条件。这种方法简化了地形的曲面形态，使得每一分块可用不同的曲面表达，同时得到光滑连续的空间曲面。不同的分块单元可以使用不同的内插函数。常用的内插函数有线性内插、双线性内插、多项式内插、样条函数、多层曲面叠加法等。逐点内插是以内插点为中心，确定一个邻域范围，用落在邻域范围内的采样点计算内插点的高程值。逐点内插本质上是局部内插，但与局部分块内插不同的是，局部内插中的分块范围一经确定，在整个内插过程中其大小、形状和位置是不变的，凡是落在该块中的内插点，都用该块中的内插函数进行计算，而逐点内插法的邻域范围大小、形状、位置乃至采样点个数随内插点的位置而变动，一套数据只用来进行一个内插点的计算。逐点内插法要注意两个问题，一是选择合适的内插函数，内插函数决定着DEM精度、DEM连续性、内插点邻域的最小采样点个数和内插计算效率。二是确定内插点邻域，内插点的邻域大小和形状、邻域内参加内插计算的数据点的个数、采样点的权重、采样点的分布、附加信息等不仅会影响到DEM的内插精度，也影响到内插速度。逐点内插方法计算简单，内插效率较高，应用比较灵活，是目前较为常用的一类DEM内插方法。在建立DEM时，要根据情况选择合适的、运算效率高的方法。而众多内插方法并不是独立的，而往往是相互结合使用，这在后续的章节里会讲到。9.3 数字地形分析地形分析是地形环境认知的一种重要手段，传统的地形分析是基于二维平面地图进行的。从基于纸质地图的地形分析发展到到基于数字地图的地形分析，计算机取代了大量的人工计算和绘制，地形分析的手段、功能发生了一次飞跃；可视化技术和虚拟现实技术的发展，使得建立三维实时、交互的仿真地形环境成为可能，同时也需要实现三维地形环境中的地形分析。特别是DEM的出现和大量应用，使得从地形属性中提取各类地形参数和特征因子更加的简便和准确。用来描述地形特征和空间分布的地形参数很多，不同的应用目的，不同的学科和领域对此的理解和分类也不同。本章将综合相关知识，着重介绍基本因子分析、地形特征提取、水文分析和可视域分析。9.3.1 基本因子分析本质上讲，DEM是地形的一个数学模型，可以看成是一个或多个函数的集合。实际上许多地形因子就是从这些函数进行一阶或二阶推导出来的，也有的通过某种组合或复合运算得到。基本地形因子包括斜坡因子（坡度、坡向、坡度变化率、坡向变化率等）、面积因子（表面积、投影面积、剖面积）、体积因子（山体体积、挖填体积）和面元因子（相对高差、粗糙度、凹凸系数、高程变异等）。本节将阐述一些常用的基本地形因子，为了方便起见，并从实际应用角度考虑，本节这些地形因子的计算都是基于格网DEM。1. 坡度x (N)y (E)zPSlopeAspectn图9.5 地表单元坡度示意图o严格地讲，地表面任一点的坡度是指过该点的切平面与水平地面的夹角。坡度表示了地表面在该点的倾斜程度，在数值上等于过该点的地表微分单元的法矢量与z轴的夹角（如图9.5所示），即：Slope = （9.1） 当具体进行坡度提取时，常采用简化的差分公式，完整的数学表示为： （9.2）式中，fx是X方向高程变化率，fy是Y方向高程变化率。地面坡度实质是一个微分的概念，地面上每一点都有坡度，它是一个微分点上的概念，是地表曲面函数z = f(x,y)在东西、南北方向上的高程变化率的函数。实际应用中，坡度有两种表示方式（如图9.6）：l 坡度(degree of slope)：即水平面与地形面之间夹角。l 坡度百分比（percent slope）：即高程增量（rise）与水平增量（run）之比的百分数。图9.6 坡度的两种表示方法拟合曲面法是解求坡度的最常用的方法。拟合曲面法，一般采用二次曲面，即在33的DEM栅格分析窗口中（如图9.7）进行，每个栅格中心为一个高程值，分析窗口在DEM数据矩阵中连续移动完成整个区域的计算工作。常用的计算fx、 fy的方法是三阶反距离平方权，该算法也用于ArcView和ARC/INFO。其计算方法为： （9.3）式（9.3），g为格网间距。2. 坡向坡向定义为：地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角（如表7.1中的坡向示意图所示，x轴为正北方向）。其数学表达公式为： （9.4）对于地面任何一点来说，坡向表征了该点高程值改变量的最大变化方向。在输出的坡向数据中，坡向值有如下规定：正北方向为0，顺时针方向计算，取值范围为0360。坡向可在DEM数据中用式9.4直接提取。但应注意，由于式9.4求出坡向有与x轴正向和x轴负向夹角之分，此时就要根据fx和fy的符号来进一步确定坡向值（如表9.2所示）图9.7 33分析窗口表9.2 坡向值的判断fxfy = Aspect坡向示意=0090y (E)360+180+180+x (N)=0-100090=0000-900360+0-900180+=001800（1）当Zi，j+1 Zi，j 则VR（i，j）= -1（2）当Zi，j+10等高线山顶点鞍部图例图9.16 利用ArcView软件及DEM数据提取的山顶、鞍部（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其中，2. 山脊线和山谷线提取山脊线和山谷线构成了地形起伏变化的分界线（骨架线），因此它对于地形地貌研究具有重要的意义。另一方面，对于水文物理过程研究而言，由于山脊、山谷分别表示分水性与汇水性，山脊线和山谷线的提取实质上也是分水线与汇水线的提取。这一特性又使得山脊线和山谷线在许多工程应用方面有着特殊的意义。在对山脊线、山谷线的提取方法中，基于规则格网DEM的方法是主要。从原理上来分，主要分为以下四种：（1） 基于图像处理技术的原理因为规则格网DEM数据事实上是一种栅格形式的数据，可以利用数字图像处理中的技术来设计算法。利用数字图像处理技术设计的算法大都采用各种滤波算子进行边缘提取。基于该原理有一种简单移动窗口的算法，其主要思路是： 计一个22窗口以对DEM格网阵列进行扫描； 第一次扫描中，将窗口中的具有最低高程值的点进行标记，自始至终未被标记的点即为山脊线上的点； 第二次扫描中，将窗口中的具有最高高程值的点进行标记，自始至终未被标记的点即为山谷线上的点。以上方法存在两个主要缺陷： 取特征点时必须排除DEM中噪声的影响； 特征点连接成线时的算法设计较为困难。（2） 基于地形表面几何形态分析原理基于地形表面几何形态分析原理的典型算法就是断面极值法。其基本思想就是地形断面曲线上高程的极大值点就是分水点，而高程的极小值点就是汇水点。该方法的基本过程为：找出DEM的纵向与横向的两个断面上的极大、极小值点，作为地形特征线上的备选点；根据一定的条件或准则将这些备选点划归各自所属的地形特征线。这种算法存在两个主要缺陷： 由于这种方法对地形特征线上的点的判定与其所属的地形特征线的判定是分开进行的，在确定地形特征线时，全区域采用一个相同的曲率阈值作为判定地形特征线上点的条件。因此它忽略了每条地形特征线必然存在的曲率变化现象。当阈值选择较大时，会丢失许多地形特征线上的点，导致后续跟踪的地形特征线间断且较短；如果选择过小，会产生地形特征线上点的误判，给后续地形特征线的跟踪带来困难。由于该方法只选择纵、横两个断面来去确定高程变化的极值点，因此它所确定的地形特征线具有一定的近似性，与实际的地形特征线有一定的差异，有时候还会出现遗漏。（3） 基于地形表面流水物理模拟分析原理的算法这种算法的基本思想是：按照流水从高至低的自然规律，顺序计算每一栅格点上的汇水量，然后按汇水量单调增加的顺序，由高到低找出区域中的每一条汇水线。根据得到的汇水线，通过计算找出各自汇水区域的边界线，就得到了分水线。算法采用了DEM的整体追踪分析的思路与方法，分析结果具有系统性好，还便于进行相应的径流成因分析，但是，该方法也存在以下两个明显的缺陷：由于该算法所计算的汇水量与高程有关，计算的结果必然是高程值大的地形特征线上的点的汇水量小，高程值小的地形特征线上的点的汇水量大。因此，可能导致低处非地形特征线上的点的汇水量也较大而被误认为地形特征线上的点；而位于高处的地形特征线上的点会因为汇水量小而被排除；这就造成用该算法所确定的地形特征线（汇水线）的两端效果很差。由于该算法降格汇水区域的公共边界视为分水线，因此它所确定的分水线均为闭合曲线，这与实际的地形特征线（山脊线）不符。（4） 基于地形表面几何形态分析和流水物理模拟分析相结合由于基于地形表面几何形态分析原理和基于地形表面流水物理模拟的算法均存在一定的缺陷，因此可将两者结合起来以实现地形特征线的提取。这种算法的基本思路是：首先采取较稀疏的DEM格网数据，按流水物理模拟算法去提取区域内概略的地形特征线；然后用其引导，在其周围邻近区域对地形进行几何分析，来精确的确定区域的地形特征线。这一算法的关键在于：求出已提取的概略的地形特征线与DEM格网线的交点，在该交点附近的一个小区域内，对DEM数据进行几何分析，即找出该区域内与概略的地形特征线正交方向地形断面上高程变化的极值点，该点即为地形特征线的精确位置。这一算法的基本过程可归纳为：概略DEM的建立；地形流水物理模拟；概略地形特征线提取；地形几何分析；地形特征线精确确定。（5） 平面曲率与坡位组合法即首先利用DEM数据提取地面的平面曲率及地面的正负地形，取正地形上平面曲率的大值即为山脊，负地形上平面曲率的大值为山谷。该种方法提取的山脊、山谷的宽度可由选取平面曲率的大小来调节，方法简便效果好。图9.17 DEM提取的山脊线图9.18 DEM提取的山谷线9.3.3 流域分析1. 流域定义降水汇集在地面低洼处，在重力作用下经常或周期性地沿流水本身所造成的槽形谷地流动，形成所谓的河流。河流沿途接纳很多支流，水量不断增加。干流和支流共同组成水系。每一个河流或每一个水系都从一部分陆地面积上获得补给，这部分陆地面积就是河流或水系的流域，也就是河流或水系在地面的集水区。把两个相邻集水区之间的最高点连接成的不规则曲线，就是两条河流或水系的分水线，因此，流域也可以说是河流分水线以内的地表范围。高程格网和栅格数据运算用于流域分析，以获取流域和河网等在水文过程中非常重要的地形要素。2. 流域提取在格网DEM实现流域地形分析，需要顺序执行如下的步骤：第一步：DEM洼地填充。由于数据噪音、内插方法的影响，DEM数据中常常包含一些“洼地”，“洼地”将导致流域水流不畅，不能形成完整的流域网络，因此在利用模拟法进行流域地形分析时，要首先对DEM数据中的洼地进行处理。填充洼地最常用的方法是之一是把 其单元值加高至周围的最低单元值。第二步：水流方向确定（flow direction）。水流方向是指水流离开格网时的流向。流向确定目前有单流向和多流向两种，但在流域分析中，常是在33局部窗口中找出八个周边单元中一个最陡的坡度在流域分析中（如图9.19），水流方向矩阵是一个基本量，这个中间结果要保存起来，后续的几个环节都要用到水流方向矩阵。原始DEM与中心格网高差中心单元流向图9.19 33窗口中心单元流向确定第三步：水流累积矩阵生成（flow accumulation）。水流累积矩阵是指流向该格网的所有的上游格网单元的水流累计量（将格网单元看作是等权的，以格网单元的数量或面积计），它是基于水流方向确定的，是流域划分的基础。流水累计矩阵的值可以是面积，也可以是单元数，取决于具体的软件，如ArcView中采用的是格网单元数。两者之间的关系是面积 = 格网单元数目 单位格网面积。无洼地DEM、水流方向矩阵、流水累计矩阵是DEM流域分析的三个基础矩阵。第四步：流域网络提取（stream networks）。流域网络是在水流累计矩阵基础上形成的，它是通过所设定的阈值，即沿水流方向将高于此阈值的格网连接起来，从而形成流域网络。9.3.4 可视性分析可视性分析也称通视分析，它实质属于对地形进行最优化处理的