p13第十三章 数字高程模型(简版).ppt

第十三章 数字高程模型DEM 表面分析 DEM的表示和生成方法DEM基本地形因子的计算方法DEM复杂地形分析 第一节 数字高程模型DEM的表示和生成方法 概念区分 数字地面模型 DigitalTerrainModel DTM 是地表形态等多种信息的一个数字表示 其中的属性包括地形 x y z 地貌 地物 自然资源 环境 社会经济等等信息的定量或定性描述 因而DTM是GIS中一种非常重要的产品 DTM也称为数字地形模型 数字高程模型 DigitalElevationModel DEM 是表示区域上地形三维向量的有限序列 Vi Xi Yi Zi Xi Yi D是平面坐标 Zi是 Xi Yi 对应的高程 对于规则格网 Vi Zi DHM DigitalHeightModel 是一个与DEM等价的概念 4D产品 数字高程模型 DEM 数字正射影像图 DigitalOrthoimageMap DOM 数字线划图 DigitalLineGraphic DLG 和数字栅格地图 DigitalRasterGraphic DRG 等 前3D为国家空间数据基础设施 NSDI 的框架数据 一 DEM的表示方法 数字地形的表达可分为 矢量点等高线三角网栅格 连续栅格 二 DEM GRID的数据来源和生成方法 DEM的来源DEM的生成方法人工网格法立体像对分析等值线插值三角网转换法高程点插值方法 DEM的数据来源 直接定位测量 GPS接收仪 普通测量设备 如经纬仪 用于小范围内各种大比例尺高精度的地形建模 这种数据获取方法的工作量很大 效率不高 费用高昂 地形图 数字化或由等高线自动生成 对于经济发达地区 由于土地开发利用使得地形地貌变化剧烈而且迅速 既有地图往往也不宜作为DEM的数据源 但对于其他经济落后地区如山区 因地形变化小 既有地图无疑是物美价廉的数据源 航空影像 航空影像测量方法可直接得到DEM 航空摄影测量一直是地形图测绘和更新最有效也是最主要的手段 其获取的影像是高精度大范围DEM生产最有价值的数据源 利用该数据源 可以快速获取或更新大面积的DEM数据 从而满足对数据现势性的要求 DEM的数据来源 遥感卫星 卫星扫描系统 如SPOT卫星上的立体扫描仪 获取的图像也能提供DEM SAR 雷达干涉波可直接得到DEM 机载激光扫描仪也可直接获取DEM 1 人工网格法 将地形图蒙上格网 逐格读取中心或交叉点的高程值 构成数字高程模型 2 摄影测量法 利用遥感立体像对 根据视差模型 通过选配左右影像的同名点 可建立数字高程模型 3 等值线插值 根据各局部等值线上的高程点 通过插值公式计算各点的高程 得到DEM 等值线插值法是比较常用的方法 输入等值线后 可在矢量格式的等值线数据基础上进行 插值效果较好 等高线插值算法生成DEM的方法 采用移动拟合加权平均插值方法 设A点为待内插的点 从A点按45 的方位间隔引出八条搜索射线 八条射线与A点相邻等高线的交点为C1 C2 Ci 其高程分别为Z1 Z2 Zi 它们到P点的距离设为d1 d2 di 则P点的插值高程Zp为 4 高程点插值方法 以不规则点图元组织的Z变量的数据 并不适合于图形显示 也不适于进行分析 高级曲面分析要求将Z值转换成一个规则间距空间格网 或者转换成不规则三角形网 栅格法可用来将不规则的空间数据转化为规则格网的空间模式 栅格法以一系列不规则分布的点图元为基础 它涉及某个有效范围内任意点Z值的内插 这个生成的点通常在在格网垂直线与水平线的交叉处 反距离加权法样条函数法趋势面法克利金法 5 三角网转换法 对有限个离散点 每三个邻近点联接成三角形 构成三角网 每个三角形代表一个局部平面 再根据每个平面的拟合方程 可计算各格网点高程 生成DEM 1 TIN的生成方法 首先取其中任一点P 在其余各点中寻找与此点距离最近的点P1 连接P和P1构成三角形的第一边 然后在其余所有点中寻找与这条边最近的点P2 找到后即构成第一个三角形 再以这个三角形新生成的两边为底边分别寻找距它们最近的点构成第二个 第三个三角形 依此类推 直到把所有的点全部连入三角网中 2 三角网的插值 双线性方法 求落在各个三角形内的网络点高程值 包括落在三角形边上的点 待求点落在三角形ABC内 先用线性插值的方法 求D E两点的值 设A B C D E P处的值分别为VA VB VC VD VE 其中VA VB VC为已知 在DEM中实质上为高程值 则D E两点处的插值为 则P点的插值为 ArcGIS的命令 TINLATTICE LINEAR QUINTIC z factor FLOAT INT 待转换的三角网 要生成的DEM LINEAR 插值方法选择线性插值 QUINTIC 插值方法采用双变量五次插值方程 z factor Z坐标的调正因子 FLOAT INT GRID的输出数据类型 生成方法小结 DEM数据采集的几点结论摄影测量是DEM的重要数据源 是进行数据库更新的最有效方法之一 现有地形图是DEM的另一重要数据源 从等高线生产DEM的方法已经完全成熟 被广泛地用于生产 使用GPS 激光扫描 干涉雷达等新型技术进行DEM数据采集是很有发展前景的方法 利用基于不规则三角网TIN的方法进行数据建模和栅格转换 是快速可靠地生产高精度DEM的切实可行方案 无论哪种方法 最大限度地采集重要的地形特征点是保证DEM质量和提高作业效率的基本前提 第二节 DEM基本地形因子计算 地形分析分类 基本地形因子计算 复杂地形分析 地形分析 坡度计算 地形水文分析 地形特征提取 可视性分析 坡向计算 粗糙度计算 可从DEM数据中通过地形分析方法得到的一些主要地形属性 据I D MOORE等 1 高程分级 高程影响地表物质和能量的分布 等间距或不等间距划分为若干高程等级 如用来区分丘陵 低山 中山 高山等 划分地貌类型 ArcGIS命令 LATTICEPOLY lookup table z factor 转换栅格中的坡度 坡向 高程范围和边界成多边形SLOPE 根据查找表中的坡度等级生成多边形 ASPECT 根据查找表中的坡向等级生成多边形 RANGE 根据查找表中的高程等级生成多边形 NODATA 根据栅格单元是否有值或空值来生成多边形 BOX 生成最小矩形区域包含输入格网中的所有栅格EXTENT 按栅格的extent生成单个矩形框 查找表 lookup table 包含分类范围和相关输出多边形代码lookupitem PERCENT SLOPE DEGREE SLOPE ASPECTorRANGE作为lookup table的索引字段 类型可以为I B F N例如PERCENT SLOPESLOPE CODEDataranges310 0 slope 3523 slope 51035 slope 1025410 slope 2540525 slope 4060640 slope 6080760 slope 80 2 平均高程 式中n的计算单元内栅格个数 h Pk 为第k点的高程 ArcGIS中使用Describe查找或者使用ZonalMean 获取 3 极值高程和高差 4 区域相对高程 设参考高程为hm 则各栅格点上相对高程为 k 1 NArcGIS下使用newdem dem heght命令 其中Height可以为常量 5 局部相对高差 格点面元是在格网DEM的水平投影面上 以四个相邻格点 i j i j 1 i 1 j 1 i 1 j 为顶点的面积范围 格点面元的相对高差指在格点面元的四个格点中 最高点与最低点之差 按下式求算 h MAX h00 h01 h02 h03 MIN h00 h01 h02 h03 式中hij i j 0 1 为四个格点的高程 ArcGIS Grid的实现方法 Focalmax Focalmin 6 剖面 Profile 用于检查和量测沿某条线的高度变化SURFACEPROFILE profile info table sample distance ArcView下的图示 剖面积计算 根据某穿越地形表面的线路 求该线路垂直剖面积 具体计算时先求得该线路与DEM格网边的所有交点Pi Xi Yi Zi 再按下面公式计算得到 n为交点数 Di i 1为Pi与Pi 1间距离 计算曲面沿线长度 SURFACELENGTH z factor sample distance surface length item 计算沿弧段的曲面长度 7 坡度和坡向 坡度和坡向是表示地表面在地面某一点处的倾斜程度和倾斜方位的一个量 它是一个矢量 既有大小又有方向 某个点的坡度是其数值等于地表曲面在该点的切平面与水平面夹角 某个点的坡向是指该切平面上沿最大倾斜方向矢量在水平面上的投影的方位角 坡度的表示 输出的坡度网格单位为度或百分数 注意 坡度表示为0 90 0为代表水平面的情况 当坡度为90 百分数为无贫大 坡度坡向的计算 拟合曲面法坡度的计算方法有很多 经证明 拟合曲面法是解求坡度的最佳方法 一般采用二次曲面 使用一个3 3窗口 每个窗口单元的中心有一高程点 坡度坡向的计算公式 坡度坡向的计算公式 例子 dz dx a 2d g c 2f i 8 x mesh spacing dz dy a 2b c g 2h i 8 y mesh spacing rise run SQRT SQR dz dx SQR dz dy degree slope ATAN rise run 57 29578 ESRI ArcGIS中的SLOPE函数 SLOPE DEGREE PERCENTRISE ArcGIS中的Aspect函数 ESRI ASPECT 注意坡向的规定 8 Curvature地表曲率 profilecurvature 坡度方向的曲率 区别凸坡和凹坡 剖面曲率影响水流的加速和减速 从而影响侵蚀和沉积过程 在凸坡区域 水流加速 水流可聚集更多的能量 输送量也增加 发生侵蚀过程 相反 在凹坡区域 水流减速 水流失去能量 发生沉积过程 Plancurvature 表示等高线方向的曲率 影响水流的聚集和分散 Converging divergingflow 以及土壤水分的含量 soilwatercontent 用于根据坡形的水平凹凸变化情况 确定沟谷线 山脊和鞍部的位置 划分流域范围 Tangentialcurvature切线曲率 水平曲率和坡度的乘积 Plancurvaturemultipliedbyslope Providesalternativemeasureoflocalflowconvergenceanddivergence 影响局部水流的聚集和分散 剖面曲率的算法原理 待算格点面元的四个格点中 最高点与其对角点的连线称为格点主轴 主轴两端点高程的平均值与格点面元平均高程的比 称为格点面元凹凸系数 CD 式中hmax为最高格点高程 h max为最高格点的对角格点高程 h为格点面元高程的平均值 当CD为正时 格点面元的实际表面为凸形坡 为负时为凹形坡 剖面曲率 表示坡面是凸坡或凹坡 可导致水流的加速或减速 分别导致侵蚀和沉积 ArcGIS计算方法 CURVATURE out profile curve out plan curve out slope out aspect 正的曲率值表示凸坡 负值表示凹坡 0表示平地 9 表面积计算 即求算格网表面的面积对于格网分解为两个三角形 然后求三角形的表面积 整个DEM的表面积为所有格网单元面积之和 对于TIN则直接求三角形面积 求三角形面积使用每点的 x y z 值 表面积计算公式 具体的计算公式为 式中 Di表示第i对三角形两顶点间的曲面距离 S表示三角形的表面积 P表面三角形的周长的一半 10 投影面积计算 投影面积指的是任意多边形在水平面上投影的面积 可直接采用海伦公式进行计算 一种简单的方法是根据梯形法则 如一个多边形由顺序排列的N个点 xi yi i 1 N 组成 并且第N点与第1点相同 则水平投影面积计算公式如下 如果多边形顶点按顺时针方向排列 则计算的面积值为负 否则为正 11 地表粗糙度 反映某一面积单元内地势伏变化的复杂程度格点面元的粗糙度 roughness 指格点面元所对应的DEM上表面积与其水平投影面积之比 记为CZ CZ S表面积 S投影面积当CZ 1时 粗糙度最小 格点面元的实际表面即为水平面 12 体积计算 DEM的体积可由四棱柱和三棱柱的体积进行累加得到 四棱柱上表面可用抛物双曲面拟合 三棱柱上表面可用斜平面拟合 下表面均为水平面或参考平面 计算公式分别为 其中S3与S4分别是三棱柱与四棱柱的底面积 体积计算的ARCGIS命令 VOLUME base value out info file z factor 输出结果 NAMEWIDTHOUTPUTTYPEN DECTIN3232C ZMIN412F3ZMAX412F3DATUM412F3AREA818F5planimetric平面面积VOLUME818F5UNITS3232C ZUNITS3232C ZFACTOR412F3 14 挖方填方计算 通过对两个输入曲面进行对比 计算出挖方区域和填方区域CUTFILL z factor 计算流程 15 山体阴影 HILLSHADE azimuth altitude ALL SHADE SHADOW z factor Azimuth 太阳方位角Altitude 太阳高度角SHADE 仅考虑localilluminationangles 不考虑遮盖的影响 输出范围在0and255 0表示黑 255表示亮 SHADOW 输出为0表示遮盖区 1表示非遮盖区 Shade选项 Shadow选项 地表辐照度 计算辐照度需考虑日照条件 太阳赤纬 高度角 时角及大气状况 与坡面几何条件的相互关系 辐照度由下式决定 式中 大气透过率 与太阳高度和大气状况有关 Sc为太阳常数 Sa为太阳高度角 可由球面三角公式求出 t是时角 a b为坡面方程系数 为坡度 16 生成等高线 TINCONTOURLATTICECONTOUR 二 地表水再分配分析 利用数字地形模型作为输入 通过确定栅格任意点上坡区域的贡献和下坡区的水流方向 勾画出水系并且对其相关特征定量化 利用DEM生成的流域和水系 是大多数水文分析的数据源 可用于决定洪水淹没范围分析 进行流域渗透水文的模拟研究 1 地表水再分配的基本过程 水系和流域的概念水系连接的类型径流过程洼地和峰的问题 1 水系和流域的概念 水系 drainagesystem 为水在地表流动并且流向出口的网络 流域 drainagebasin watershed basin catchment 是地表水和其他物质流向统一出口点 outlet pourpoint 的区域 出口点为流域边界的最低点 两个流域的边界被定义为分水岭 drainage 2 水系连接的类型 水系可以理解为树 出口点为树根 树干为河沟 河沟的交点被称为节点或交点 nodeorjunction 连接两个内流的节点 或连接节点和出口点的河沟叫做内流连接 interiorlinks 外连接 Exteriorlinks 定义为树最外面的分支 3 径流过程 水流的方向取决于每个点的坡向 水流的能量取决于地表坡度 陡的坡度产生大的能量 可以输送更多和更大的固体物质 并具有较大的侵蚀潜力 水平曲率用于区别斜坡是山脊或山沟 山脊和山沟分别导致水流的聚集和分散 4 洼地 sinks 和峰 peaks DEMS中的错误通常划分为洼地 sinks 和峰 peaks 洼地 sinks 是比较高的邻域所包围的区域 洼地可以是自然形成的 峰 peaks 是被比较低的邻域所包围的地区 通常是自然地貌 洼地的数量在低分辨率时较多 另外当珊格的数据格式采用整形时 在低高差的地区洼地数量也会增多 当计算流向时 洼地会产生不合适的结果 在提取流向时剔除 DEM可能包含明显的水平条带 这来源于生成dem时的系统取样误差 这在平坦地表的整形珊格中非常明显 2 流域河网 流域边界提取算法 基于栅格的地形分析技术在水文学中的应用 一般采用了O Callaghan和Mark的坡面流模拟方法 Jense和Domingue 1988 Martz和DeJong 1988 Garbrecht 1997 在此基础上做了改进 这种算法一般称之为D8 Deterministiceight neighbours 算法 D8算法示例 D8模型 高程格网示例 格网流向模型 水流聚集模型 算法流程图 1 洼地的填充 洼地的出现可导致流向的计算错误 在某些情况下 可能有合理的洼地存在 这就需要对当地的地形情况有较多了解 以辨别真的洼地和错误 洼地和流向不确定的情况 1 周围单元格的值高于处理的单元格的值 2 两个单元格间的水流互相流动 3 单元格在多个方向上有同样的高差 要精确的表示流向 应该对洼地进行填充 洼地可以用Sink函数定位 并进行填充 填充洼地的过程 计算洼地深度 1 利用SINK函数识别洼地Grid sinks sink flowdir 2 利用WATERSHED函数生成每个洼地代表的流域区域Grid sink areas watershed flowdir sinks 3 生成每个洼地流域范围内的最低高程值栅格Grid sink min zonalmin sink areas elevation 4 将洼地填充到沿洼地边界的最低高程值 使水流能流出洼地 表示为Sink maxGrid sink max zonalfill sink areas elevation 5 最低高程值网格减去最低高程值网格得到洼地深度栅格Grid sink depth sink max sink min洼地深度的计算便于下述fill指令确定z limit的取值 填充洼地的过程 填充 FILL SINK PEAK z limit out dir grid 表示连续曲面的网格 剔除洼地或山顶的输出网格 SINK PEAK 表示洼地将被填充或山顶将被削去 z limit 洼地和其出口点的最大高差或山顶与其相邻最高单元格之间的高差 洼地深度大于这个值 则该洼地不被填充 如果山顶与其相邻最高单元格之间的高差大于这个值 则山顶不被移去 默认情况下时填充所有洼地或移去所有山顶 而不考虑深度 out dir grid 可选输出流向栅格 2 FlowDirection 按流向选最大坡度的原理 计算每个单元格流向其邻域的方向 流出方向 如果到相邻栅格的坡度相同 则放大邻域以选择最陡的方向 流向的规定 A 判断指标 距离权落差 距离权落差按下式计算 对角方向距离D为sqrt 2 格网尺寸 其它方向D为格网尺寸 B 确定水流方向步骤 1 对所有DEM边缘的格网 赋以指向边缘的方向值 可增加判断流向内部格网 2 计算中心格网对8个邻域网的距离权落差值 3 确定具有最大落差值的格网 分以下几种情况分别处理 如果最大落差 0 则赋负值 如果最大落差 0 赋最大值方向编码 如果最大落差 0 有一个以上最大值 按逻辑判断作为流向编码 如果最大落差 0 则0值对应的方向值相加 作为其方向代码 4 对没有赋以负值或水流方向值为1 2 4 128 的格网 检查最大落差邻域格网 如果此格网正确编码 且方向没有指向中心格网 则以此格网的方向值作为中心格网的方向值 5 重复第4步 直至没有任何格网能被赋以方向值 对方向值不为1 2 4 128的格网赋以负值 ArcGIS流向命令 Grid flow dir flowdirection elevation 3 FlowAccumulation 累积流表示流入该单元格的栅格数目 高的累积流可用于生成水系网 低的累积流表示局部的地形高区 可用于识别山脊 ArcGIS中累积流的计算 FLOWACCUMULATION函数计算流入每一个上坡方向的栅格累积值作为这个栅格的累积流 也可以考虑权重因子 比如降雨量 Grid flow acc flowaccumulation flow dir rainfall 4 流域河网的提取 确定上游集水区面积阈值 标注集水量累计值大于阈值的格网 将得到的图形进行矢量化处理 通过阈值调整控制生成河网密度 水系网提取的ArcGIS 水系网可从FLOWACCUMULATION函数中得到的结果中勾画出来 通过使用地图代数表达式对累积流网格设置阈值 可以勾画出水系网 例如 要生成一个在无值背景上值1代表的水系网 可以使用 Grid stream net con flowacc 100 1 Grid stream net setnull flowacc 100 1 上述单元格取值大于100的所有单元格被赋值为1 其它单元格被赋值为无值NODATA 5 水系网级别的划分 水系定级是一种给水系网中的连接赋级别代码 目的是按水系的支流数识别和分类河谷类别 水系的某些特征可以通过其级别推断出来 例如 一级支流仅受坡面漫流的控制 它没有其它支流的水流贡献 它们极易受非点源污染的影响 Usage STREAMORDER STRAHLER SHREVE STRAHLER SHREVE STRAHLER外连接总被赋值为1 但当同样级别的水系交叉时 水系级别增加 所以 两个第二级的水系交叉时 会生成第三级的水系 但当两个不同级别的水系交叉时 则不会产生水系别的增加 例如 第一和第二级的水系交叉不会产生第三级的水系 而保持较高级别的水系等级 Shreve法计算水系网中的所有连接数 外连接被赋值为1 但对内连接其级别数是增加的 例如 两个一级支流交叉时会产生二级支流 而一级支流和二级支流交叉时则生成三级支流 6 水系矢量化 STREAMLINE out item weed 7 ArcGIS流域的划分 watersheds watershed flow dir pour points BASIN ArcGIS中的流域命令 将流向网格作为输入 可以利用WATERSHED或BASIN函数生成流域 Grid watersheds watershed flow dir pour points Grid wsheds watershed flowdirection elev selectpoint elev 对任意单元格可以勾画流域 在WATERSHED函数中可以利用SELECTPOINT函数标记流域出口点 pourpoints 在交互选择流域出口点时 显示累积流网格 有助于选择主要的水系分段 否则 可能选出接近主水系的网格 生成的流域要比期望的小 小结 ArcGIS中的处理流程 3 其它基