第9章DEM与数字地形

第9章DEM与数字地形分析朱莹主要内容基本概念DEM建立数字地形分析基本概念数字高程模型数字高程模型（DigitalElevationModel，简称DEM）是通过有限的地形高程数据实现对地形曲面的数字化模拟（即地形表面形态的数字化表示），它是对二维地理空间上具有连续变化特征地理现象的模型化表达和过程模拟DEM也常常称为DTM（DigitalTerrainModel）从狭义角度定义，DEM是区域表面海拔高程的数字化表达从广义角度定义，DEM是地理空间中地理对象表面海拔高度的数字化表达，该定义将描述对象不再限定在“地表面”，因而具有更大的包容性，有海底DEM、下伏岩层DEM、大气等压面DEM等基本概念数学意义上的数字高程模型是定义在二维空间上的连续函数DEM的数学定义为区域D的采样点或内插点Pj按某种规则连接成的面片M的集合DEM按照结构，可分为规则格网DEM、TIN、基于点的DEM和基于等高线的DEM等由于规则格网结构简单，算法设计明了，在实际运用中被广泛采用数字地形分析数字地形分析DigitalTerrainAnalysis,DTA，是指在数字高程模型上进行地形属性计算和特征提取的数字信息处理技术DTA技术是各种与地形因素相关空间模拟技术的基础数字地形分析根据地形要素的关系特征和计算特征，地形属性可分为地形曲面参数（parameters）具有明确的数学表达式和物理定义，并可在DEM上直接量算，如坡度、坡向、曲率等地形形态特征（features）是地表形态和特征的定性表达，可以在DEM上直接提取，其特点是定义明确，但边界条件有一定的模糊性，难以用数学表达式表达，如在实际的流域单元的划分中，往往难于确定流域的边界数字地形分析地形统计特征（statistics）指给定地表区域的统计学上的特征复合地形属性（compoundattributes）在地形曲面参数和地形形态特征的基础上，利用应用学科（如水文学、地貌学和土壤学）的应用模型而建立的环境变量，通常以指数形式表达数字地形分析数字地形分析的方法提取坡面地形因子地形定量因子是为有效地研究与表达地貌形态特征所设定的具有一定意义的参数或指标常用的坡面地形因子有坡度、坡向、平面曲率、坡面曲率、地形起伏度、粗糙度、切割深度等数字地形分析数字地形分析的方法提取特征地形要素流域分析主要是根据地表物质运动的特性，特别是水流运动的特点，利用水流模拟的方法来提取水系、山脊线、谷底线等地形特征线，并通过线状信息分析其面域特征可视性分析包括两方面内容，一个是两点之间的通视性(Intervisibility)，另一个是可视域(ViewShed)，即对于给定的观察点所覆盖的区域数字地形分析数字地形分析的方法地形统计特征分析地形统计分析是应用统计方法对描述地形特征的各种可量化的因子或参数进行相关、回归、趋势面、聚类等统计分析，找出各因子或参数的变化规律和内在联系，并选择合适的因子或参数建立地学模型数字地形分析方法坡面地形因子提取特征地形要素提取地形统计特征分析流域分析可视域分析提取坡度………提取坡面曲率数字地形分析常用方法DEM建立数字高程模型的建立过程是一个模型建立过程，建模的目的是对复杂的客体进行简化和抽象，并把对客体（源域，DEM中为地形起伏）的研究转移到对模型的研究上来首先要为模型构造一个合适的空间结构spatialframework空间结构是为把特定区域内的空间目标镶嵌在一起而对区域进行的划分，划分出的各个空间范围称为位置区域或空间域。空间结构一般是规则的（如格网），或不规则的（如不规则三角网TIN）DEM建立建立在空间结构基础上的模型是由n个空间域的有限集合组成由于空间数据包含位置特征和属性特征，而属性特征是定义在位置特征上的，因此每一个空间域就是由空间结构到属性域的计算函数或域函数模型构建的内容和过程：采用合适的空间模型构造空间结构；采用合适的属性域函数；在空间结构中进行采样，构造空间域函数；利用空间域函数进行分析DEM建立对于数字高程模型而言，空间结构的构造过程即为DEM的格网化过程（形成格网），属性值为高程，构造空间域函数即为内插函数的确定，利用空间域函数进行分析就是求取格网点的函数值规则格网DEM建立DEM是在二维空间上对三维地形表面的描述构建DEM的整体思路首先在二维平面上对研究区域进行格网划分（格网大小取决于DEM的应用目的），形成覆盖整个区域的格网空间结构然后利用分布在格网点周围的地形采样点内插计算格网点的高程值最后按一定的格式输出，形成该地区的格网DEMx不规则分布点规则分布等高线分布Y对每一格网点求取格网点高程格网DEM建立流程DEM内插方法DEM构建的关键环节是根据采样点的值内插计算格网点上的高程值内插是指根据分布在内插点周围的已知参考点的高程值求出未知点的高程值，它是DEM的核心问题，贯穿于DEM的生产、质量控制、精度评定、分析应用的各个环节

DEM内插分类方法DEM内插数据分布规则分布内插方法不规则分布内插方法等高线数据内插方法内插范围整体内插方法局部内插方法逐点内插方法内插曲面与参考点关系纯二维内插曲面拟合内插内插函数性质多项式内插线性插值双线性插值高次多项式插值样条内插有限元内插最小二乘配置内插地形特征理解克立金内插多层曲面叠加内插加权平均值内插分形内插傅立叶级数内插整体内插在整个区域用一个数学函数来表达地形曲面整体内插函数通常是高次多项式，要求地形采样点的个数大于或等于多项式的系数数目优点：整个区域上函数的唯一性、能得到全局光滑连续的DEM、充分反映宏观地形特征等缺点：保凸性较差、不容易得到稳定的数值解、多项式系数的物理意义不明显、解算速度慢且对计算机容量要求较高、不能提供内插区域的局部地形特征等一般与局部内插方法配合使用，例如在使用局部内插方法前，利用整体内插去掉不符合总体趋势的宏观地物特征。另外也可用来进行地形采样数据中的粗差检测局部分块内插将地形区域按一定的方法进行分块，对每一分块，根据其地形曲面特征单独进行曲面拟合和高程内插一般按地形结构线或规则区域进行分块，分块的大小取决于地形的复杂程度、地形采样点的密度和分布为保证相邻分块之间的曲面平滑连接，相邻分块之间要有一定宽度的重叠，或者对内插曲面补充一定的连续性条件简化了地形的曲面形态，使得每一分块可用不同的曲面表达，同时得到光滑连续的空间曲面不同的分块单元可以使用不同的内插函数常用的内插函数有线性内插、双线性内插、多项式内插、样条函数、多层曲面叠加法等逐点内插以内插点为中心，确定一个邻域范围，用落在邻域范围内的采样点计算内插点的高程值逐点内插本质上是局部内插，但与局部分块内插不同的是，局部内插中的分块范围一经确定，在整个内插过程中其大小、形状和位置是不变的，而逐点内插法的邻域范围大小、形状、位置乃至采样点个数随内插点的位置而变动逐点内插逐点内插法要注意的问题：一是选择合适的内插函数，内插函数决定着DEM精度、DEM连续性、内插点邻域的最小采样点个数和内插计算效率二是确定内插点邻域，内插点的邻域大小和形状、邻域内参加内插计算的数据点的个数、采样点的权重、采样点的分布、附加信息等不仅会影响到DEM的内插精度，也影响到内插速度数字地形分析基本因子分析DEM是地形的一个数学模型，可以看成是一个或多个函数的集合。实际上许多地形因子就是从这些函数进行一阶或二阶推导出来的，也有的通过某种组合或复合运算得到基本地形因子包括斜坡因子（坡度、坡向、坡度变化率、坡向变化率等）、面积因子（表面积、投影面积、剖面积）、体积因子（山体体积、挖填体积）和面元因子（相对高差、粗糙度、凹凸系数、高程变异等）基本因子分析坡度地表面任一点的坡度是指过该点的切平面与水平地面的夹角x

(N)y(E)zPSlopeAspectn地表单元坡度示意图o坡度坡度表示了地表面在该点的倾斜程度，在数值上等于过该点的地表微分单元的法矢量与z轴的夹角slope=当具体进行坡度提取时，常采用简化的差分公式，完整的数学表示为fx是X方向高程变化率，fy是Y方向高程变化率坡度实际应用中坡度的两种表示方式坡度(degreeofslope)：即水平面与地形面之间夹角坡度百分比（percentslope）：即高程增量（rise）与水平增量（run）之比的百分数坡度求坡度的常用方法拟合曲面法，一般采用二次曲面，即在3×3的DEM栅格分析窗口中进行，每个栅格中心为一个高程值，分析窗口在DEM数据矩阵中连续移动完成整个区域的计算工作坡度常用的计算fx、fy的方法是三阶反距离平方权，该算法也用于ArcView和ARC/INFOg为格网间距坡向坡向定义为地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角，其数学公式为对于地面任何一点来说，坡向表征了该点高程值改变量的最大变化方向在输出的坡向数据中，坡向值规定正北方向为0°，顺时针方向计算，取值范围为0°～360°坡向坡向可在DEM数据中用公式直接求取。但由公式求出坡向有与x轴正向和x轴负向夹角之分，此时就要根据fx和fy的符号来进一步确定坡向值采用这种方法求取的坡向分级比较详细，但实际应用中往往需要归并，在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°）注：上述情况假定所建立的DEM数据从南向北获取的，且x轴与正北方向重合，否则上述公式求得的坡向值，还应加上x轴偏离正北方向的夹角值。180+α0～90<01800=0180+α-90～0>0<0360+α-90～0<000=0α0～90>0>0270∕<0-1∕=090∕>0=0坡向示意Aspectfyfx坡向值的判断y(E)360+α180+α180+ααααααx(N)α=原始DEM数据及实验区等高线图ARCVIEW软件下提取的坡度图ARCVIEW软件下提取的坡度图由DEM提取的坡向图曲率曲率是对地形表面一点扭曲变化程度的定量化度量因子，地面曲率在垂直和水平两个方向上分量分别称为平面曲率和剖面曲率地形表面曲率反映了地形结构和形态，同时也影响着土壤有机物含量的分布，在地表过程模拟、水文、土壤等领域有着重要的应用价值和意义剖面曲率是对地面坡度的沿最大坡降方向地面高程变化率的度量，数学表达式为：曲率平面曲率指在地形表面上，具体到任何一点P，过该点的水平面沿水平方向切地形表面所得的曲线在该点的曲率值平面曲率描述的是地表曲面沿水平方向的弯曲、变化情况，也就是该点所在的地面等高线的弯曲程度平面曲率示意图x（N）y（E）zPV水平方向切地表面所得曲线曲率平面曲率公式p是x方向高程变化率q是y方向高程变化率r为x方向高程变化率的变化率s为x方向高程变化率在y方向的变化率t为y方向高程变化率的变化率曲率曲率因子的提取算法的基本原理在DEM数据的基础上，根据其离散的高程数值，把地表模拟成一个连续的曲面，从微分几何的思想出发，模拟曲面上每一点所处的垂直于和平行于水平面的曲线，利用曲线曲率的求算方法的推导得出各个曲率因子的计算公式在DEM中求算高程的微分分量有一套独特的算法，最常用是三阶反距离平方权差分IhgfedcbaI'xh'xg'xf'xe'xd'xc'xb'xa'xI'yh'yg'yf'ye'yd'yc'yb'ya'yq值矩阵p值矩阵地面曲率提取步骤流程图KV剖面曲率Kh平面曲率ArcView提取的剖面曲率ArcView提取的平面曲率宏观地形因子地形起伏度、地形表面粗糙度与地表切割深度等地形因子是描述和反映地形表面较大区域内地形的宏观特征，在较小的区域内并不具备任何地理和应用意义宏观地形因子对于在宏观尺度上的水土保持、土壤侵蚀特征、地表发育、地貌分类等研究中具有重要的理论意义基于栅格DEM计算宏观地形因子时，关键在于确定分析半径的大小，确定一个合适的分析窗口半径或分析区域，使得求取的宏观因子能够准确反映地面的起伏状况与水土流失特征，是提取算法的核心步骤和决定信息提取效果与有效性的关键宏观地形因子地形起伏度地形起伏度是指，在所指定的分析区域内所有栅格中最大高程与最小高程的差在区域性研究中，利用DEM数据提取地形起伏度能够直观的反映地形起伏特征在水土流失研究中，地形起伏度指标能够反映水土流失类型区的土壤侵蚀特征，比较适合区域水土流失评价的地形指标宏观地形因子地形粗糙度地表粗糙度，一般定义为地表单元的曲面面积S曲面与其在水平面上的投影面积S水平之比地表粗糙度能够反映地形的起伏变化和侵蚀程度基于DEM的地表粗糙度的提取主要分为以下两个步骤：根据DEM提取坡度因子S；根据公式R=1/cos(S)计算地表粗糙度宏观地形因子地表切割深度地表切割深度是指地面某点的邻域范围的平均高程与该邻域范围内的最小高程的差值地表切割深度直观的反映了地表被侵蚀切割的情况，并对这一地学现象进行了量化，是研究水土流失及地表侵蚀发育状况时的重要参考指标地形特征分析地形特征提取主要是指地形特征点、线、面的提取，并进而通过基本要素的组合进行地表形态分析地形特征点提取地形特征点主要包括山顶点（peak）、凹陷点（pit）、脊点（ridge）、谷点（channel）、鞍点（pass），平地点（plane）等利用DEM提取地形特征点，可通过一个3×3或更大的栅格窗口，通过中心格网点与8个邻域格网点的高程关系来进行判断该方法假设DEM表面为z=f(x,y),但由于真实地表与数学表面的差别，在利用该方法在DEM上提取特征点，结果常产生伪特征点山顶点是在局部区域内各方向上都没有凹凸性变化的点。平地点（plane）是指在两个相互正交的方向上，一个方向凸起，而另一个方向凹陷的点。鞍点（pass）是指在两个相互正交的方向上，一个方向凹陷，而另一个方向没有凹凸性变化的点。谷点（channel）是指在两个相互正交的方向上，一个方向凸起，而另一个方向没有凹凸性变化的点。脊点（ridge）是指在局部区域内海拔高程的极小值点，表现为在各方向上都为凹陷。凹陷点（pit）是指在局部区域内海拔高程的极大值点，表现为在各方向上都为凸起。山顶点（peak）邻域高程关系定义名称，，，，，或，或，，，或在局部区域内利用x，y方向的凹凸性判断地形特征分析在一个3×3的栅格窗口中，也可以直接利用中心格网点与8个邻域格网点的高程关系来进行判断地形特征点(i,j)(i-1,j-1)(i-1,j)(i-1,j+1)(i,j+1)(i+1,j+1)(i+1,j)(i+1,j-1)(i,j-1)差分算法示意图地形特征分析如果（Zi，j-1-Zi，j）（Zi，j-1-Zi，j）>0当Zi，j+1>Zi，j则VR（i，j）=-1 当Zi，j+1<Zi，j则VR（i，j）=1如果（Zi-1，j-Zi，j）（Zi+1，j-Zi，j）>0当Zi+1>Zi，j则VR（i，j）=-1当Zi+1<Zi，j则VR（i，j）=-1如果上述条件同时成立，则VR（i，j）=2如果以上条件都不成立，则VR（i，j）=0其中等高线山顶点鞍部图例利用ArcView软件及DEM数据提取的山顶、鞍部地形特征分析山脊线和山谷线提取山脊线和山谷线构成了地形起伏变化的分界线（骨架线）对于水文物理过程研究而言，由于山脊、山谷分别表示分水性与汇水性，山脊线和山谷线的提取实质上也是分水线与汇水线的提取地形特征分析对山脊线、山谷线的提取方法基于图像处理技术的原理基于地形表面几何形态分析原理基于地形表面流水物理模拟分析原理的算法基于地形表面几何形态分析和流水物理模拟分析相结合平面曲率与坡位组合法地形特征分析基于图像处理技术原理的主要思路：设计一个2×2窗口以对DEM格网阵列进行扫描；第一次扫描中，将窗口中的具有最低高程值的点进行标记，自始至终未被标记的点即为山脊线上的点；第二次扫描中，将窗口中的具有最高高程值的点进行标记，自始至终未被标记的点即为山谷线上的点以上方法存在两个主要缺陷：取特征点时必须排除DEM中噪声的影响；特征点连接成线时的算法设计较为困难地形特征分析基于地形表面几何形态分析原理的典型算法就是断面极值法基本思想是地形断面曲线上高程的极大值点就是分水点，而高程的极小值点就是汇水点基本过程为：找出DEM的纵向与横向的两个断面上的极大、极小值点，作为地形特征线上的备选点；根据一定的条件或准则将这些备选点划归各自所属的地形特征线地形特征分析缺点由于这种方法对地形特征线上的点的判定与其所属的地形特征线的判定是分开进行的，在确定地形特征线时，全区域采用一个相同的曲率阈值作为判定地形特征线上点的条件。当阈值选择较大时，会丢失许多地形特征线上的点，导致后续跟踪的地形特征线间断且较短；如果选择过小，会产生地形特征线上点的误判，给后续地形特征线的跟踪带来困难由于该方法只选择纵、横两个断面来去确定高程变化的极值点，因此它所确定的地形特征线具有一定的近似性，还会出现遗漏地形特征分析基于地形表面流水物理模拟分析原理的算法基本思想按照流水从高至低的自然规律，顺序计算每一栅格点上的汇水量，然后按汇水量单调增加的顺序，由高到低找出区域中的每一条汇水线根据得到的汇水线，通过计算找出各自汇水区域的边界线，就得到了分水线地形特征分析基于地形表面流水物理模拟分析原理方法的缺点由于该算法所计算的汇水量与高程有关，计算的结果必然是高程值大的地形特征线上的点的汇水量小，高程值小的地形特征线上的点的汇水量大。因此，可能导致低处非地形特征线上的点的汇水量也较大而被误认为地形特征线上的点；而位于高处的地形特征线上的点会因为汇水量小而被排除；造成用该算法所确定的地形特征线（汇水线）的两端效果很差由于该算法降格汇水区域的公共边界视为分水线，因此它所确定的分水线均为闭合曲线，这与实际的地形特征线（山脊线）不符地形特征分析由于基于地形表面几何形态分析原理和基于地形表面流水物理模拟的算法均存在一定的缺陷，因此将两者结合起来实现地形特征线的提取基本思路首先采取较稀疏的DEM格网数据，按流水物理模拟算法去提取区域内概略的地形特征线然后用其引导，在其周围邻近区域对地形进行几何分析，来精确的确定区域的地形特征线关键求出已提取的概略的地形特征线与DEM格网线的交点，在该交点附近的一个小区域内，对DEM数据进行几何分析，即找出该区域内与概略的地形特征线正交方向地形断面上高程变化的极值点，该点即为地形特征线的精确位置地形特征分析基本过程：概略DEM的建立；地形流水物理模拟；概略地形特征线提取；地形几何分析；地形特征线精确确定地形特征分析平面曲率与坡位组合法首先利用DEM数据提取地面的平面曲率及地面的正负地形，取正地形上平面曲率的大值即为山脊，负地形上平面曲率的大值为山谷。该种方法提取的山脊、山谷的宽度可由选取平面曲率的大小来调节，方法简便效果好DEM提取的山脊线图9.18DEM提取的山谷线流域分析两个相邻集水区之间的最高点连接成的不规则曲线，就是两条河流或水系的分水线，因此，流域也可以说是河流分水线以内的地表范围格网DEM实现流域地形分析的步骤DEM洼地填充水流方向确定（flowdirection）水流累积矩阵生成（flowaccumulation）流域网络提取（streamnetworks）流域分析洼地填充