《地理信息系统》课程教学资源（Arcview讲义）第十四章 ArcView DEM地形分析精度.doc

第十四章ArcViewDEM地形分析精度能有效地利用DEM数据进行地形定量因子的自动提取，是ArcVieWGIS软件空间分析模块的重要功能。数字高程模型（DEM）是地理信息系统地理数据库中最为重要的空间信息资料和赖以进行地形分析的核心数据系统。目前世界各主要发达国家都纷纷建立了覆盖全国的DEM数据系统，DEM已经在测绘、资源与环境、灾害防治、国防等与地形分析有关的科研及国民经济各领域发挥着越来越巨大的作用。但是，由于DEM原始信息源精度、DEM空间分辨率、以及研究区地形复杂度的差异，DEM所提取的地形因子的精度存在着相当大的差异。本章重点介绍DEM地形分析精度与不确定性方面的部分研究成果。第一节DEM地形描述误差的量化模拟一、DEM地形描述误差的概念：DEM精度是指所建立的DEM对真实地面描述的准确程度。DEM误差的大小被普遍视为衡量DEM精确性的标准。然而，人们在该问题上存在着明显的片面认识。以往的研究普遍重视在DEM采样点上出现的高程采样误差，而相对忽视由于DEM离散采样所造成的地形描述误差。无疑，高程采样误差是影响DEM精度的重要因素，但决不是唯一因素。因为，即使DEM在所有高程采样点上的误差均为零，有限的DEM栅格采样点所DEM模拟地面构成的高程模型也只能是对实际地面的近似模拟。我们将E这种在假定DEM高程采样误差为零条件下，模拟地面与CB实际地面Etr实际地面之差异，定义为DEM地形描述误差（以后简称atDEt）。如图1所示，A、B两点为DEM地面高程采样点，A、B两点的连线为DEM模拟地面，假定在该两点的高程采样误差为零，则Etc、Etp及Ete分别为在C、D、E三点的地形描述误差。无疑，DEM栅格分辨率与地形起伏的栅格分辨率复杂程度是影响Et大小的两个关键因子，建立该2因子图14-1.Et分布示意图与DEM地形描述误差之间的量化关系，是对误差进行定量模拟的关键。二。实验样区与原始DEM数据精度在我国选择有代表性的6个不同地面复杂度的地区作为试验样区，实验区面积均为10km×10km，其主要地形因子及原始信息源精度如表1所示。采用全数字摄影测量方法所建立的DEM作为基本信息源，DEM栅格分辨率均为5m。在每个实验区的地形图上随机选择50个左右高程控制点，并视其高程值为准值，分别对各个实验区所建立DEM的高程采样精度进行测定。表1显示该组DEM具有较高的高程采样精度，便于作为基本信息源进行DEM地形描述误差的研究

第十四章 ArcView DEM 地形分析精度能有效地利用 DEM 数据进行地形定量因子的自动提取，是 ArcView GIS 软件空间分析模块的重要功能。数字高程模型（DEM）是地理信息系统地理数据库中最为重要的空间信息资料和赖以进行地形分析的核心数据系统。目前世界各主要发达国家都纷纷建立了覆盖全国的 DEM 数据系统， DEM 已经在测绘、资源与环境、灾害防治、国防等与地形分析有关的科研及国民经济各领域发挥着越来越巨大的作用。但是，由于 DEM 原始信息源精度、DEM 空间分辨率、以及研究区地形复杂度的差异，DEM 所提取的地形因子的精度存在着相当大的差异。本章重点介绍 DEM 地形分析精度与不确定性方面的部分研究成果。第一节 DEM 地形描述误差的量化模拟一、 DEM 地形描述误差的概念: DEM 精度是指所建立的 DEM 对真实地面描述的准确程度。DEM 误差的大小被普遍视为衡量 DEM 精确性的标准。然而，人们在该问题上存在着明显的片面认识。以往的研究普遍重视在 DEM 采样点上出现的高程采样误差，而相对忽视由于 DEM 离散采样所造成的地形描述误差。无疑，高程采样误差是影响 DEM 精度的重要因素，但决不是唯一因素。因为，即使 DEM 在所有高程采样点上的误差均为零，有限的 DEM 栅格采样点所构成的高程模型也只能是对实际地面的近似模拟。我们将这种在假定 DEM 高程采样误差为零条件下，模拟地面与实际地面之差异，定义为 DEM 地形描述误差（以后简称 Et）。如图 1 所示, A、B 两点为 DEM 地面高程采样点，A、 B 两点的连线为 DEM 模拟地面，假定在该两点的高程采样误差为零，则 EtC 、EtD 及 EtE分别为在 C、D、E 三点的地形描述误差。无疑，DEM 栅格分辨率与地形起伏的复杂程度是影响 Et 大小的两个关键因子，建立该 2 因子与 DEM 地形描述误差之间的量化关系，是对误差进行定量模拟的关键。二．实验样区与原始 DEM 数据精度在我国选择有代表性的6个不同地面复杂度的地区作为试验样区，实验区面积均为10km×10km，其主要地形因子及原始信息源精度如表 1 所示。采用全数字摄影测量方法所建立的 DEM 作为基本信息源，DEM 栅格分辨率均为 5m。在每个实验区的地形图上随机选择 50 个左右高程控制点，并视其高程值为准值，分别对各个实验区所建立 DEM 的高程采样精度进行测定。表 1 显示该组 DEM 具有较高的高程采样精度，便于作为基本信息源进行 DEM 地形描述误差的研究。 A C B EtD D E EtE 栅格分辨率实际地面 DEM模拟地面图 14-1. Et 分布示意图

表14-1试验区主要地形因子及DEM精度平原低丘丘陵中山高山混合类型地理位置关中平原东北曼岗丘陵江南丘陵北京军都山秦岭首阳山陕西骊山实验区中心点109°2704"E126°2119"E114°30°28E116°1916"E108°2509"E109°10'31"E试验区地理坐标34°35'12"N47°12'18"N27°25'46"N40°32'31"N33°58"15"N34°21'51"N平均高程 (m)4252242278242614662地形因子平均坡度(m)2.427.1515.120.727.514.315.4721.24剖面曲率 (度)5.968.8734.8018.48原始DEM均方差（m）0.390.641.151.522.821.350.280.571.041.411.23标准差（m）2.16精度平均误差（m）0.240.410.911.032.091.11三，提取DEM地形描述误差Et的方法栅格中点的高程与该栅格四个角点高程平均高程之差可以被定义为该栅格的地形描述误差。因此，采用栅格窗口I.如图2所示，对于DEM单元栅格分析法实现Et的提取。O为栅格中心点，A、B、C、D、O为对应地面点abcd，位，假定在该栅格四个采样点A、B、C、D上的高程采样误差均为零；Ha,HB.Hc,Hp及Ho分别为在A,B,C,D及O各点的高程，O'为A，B.CD四点的平均高程位置，这样，我们将O点与O点的高差作为该栅格的地形描述误差Et，即：Et =Ho-Ho-=Ho-(Ha+HB +Hc+HD)/4(1)按照公式（1)，已知单元栅格中点与其周围四个相邻点的图2.DEM高程采样栅格元高程是获取该单元Et值的必要条件。由于单元栅格中点的实Fig.14-2ADEM grid unit model际高程是未知的，在实际计算中采用窗口分析法实现E的提取。对于3x3的正方形分析窗口，如果DEM空间分辨率为d，在该分析窗口内，i行i列Et值可以通过下式求得Et(.j) = H(.j) - (H(-1,j-1)+H(-1,j+1)+H(+1,j-1)+H(i+1, j+1))/4(2)显然，3x3窗口的分析分辨率为2d（见图14-3）。当顺序移动该分析窗口对整个DEM进行逐点Et值计算，即能获得相应的误差矩阵。当将分析窗口依次扩大到5x5、7×7...，可以同理提取分析分辨率分别等于4d、6d、....的误差矩阵。对于5m分辨率的DEM，其相应的分析分辨率依次为10m、20m、30m、.....100m。(i+1, j-1)(i+1,j+1)(i+2.-2)(i+2, j+2)+I-1--1-(i.i)(i)1(i-1, j+1)(i-1, j-1)(i-2, j-2)(i-2, j+2)4d图14-3栅格分析窗口示意图四、DEM地形描述误差的量化模拟根据上述方法所得到的不同分辨率条件下的误差矩阵，便可以利用统计与比较分析的方法揭示Et随分辨率及地形复杂度的变化而变化的规律，表14-2及图14-4为对误差矩阵的数据进行统计的结

表 14-1 试验区主要地形因子及 DEM 精度平原低丘丘陵中山高山混合类型试验区地形因子地理位置关中平原东北曼岗丘陵江南丘陵北京军都山秦岭首阳山陕西骊山实验区中心点地理坐标 109°27’04’’E 34°35’12’’N 126°21’19’’E 47°12’18’’N 114°30’28’’E 27°25’46’’N 116°19’16’’E 40°32’31’’N 108°25’09’’E 33°58’15’’N 109°10’31’’E 34°21’51’’N 平均高程 (m) 425 224 227 824 2614 662 平均坡度 (m) 2.42 7.15 15.1 20.7 27.5 14.3 剖面曲率 (度) 5.96 8.87 15.47 21.24 34.80 18.48 原始 DEM 精度均方差（m） 0.39 0.64 1.15 1.52 2.82 1.35 标准差（m） 0.28 0.57 1.04 1.41 2.16 1.23 平均误差（m） 0.24 0.41 0.91 1.03 2.09 1.11 三．提取 DEM 地形描述误差 Et 的方法栅格中点的高程与该栅格四个角点高程平均高程之差, 可以被定义为该栅格的地形描述误差。因此，采用栅格窗口分析法实现 Et 的提取。如图 2 所示，对于 DEM 单元栅格 abc d， O´´为栅格中心点，A、B、C、D、O 为对应地面点位，假定在该栅格四个采样点 A、B、C、D 上的高程采样误差均为零；HA , HB , HC , HD 及 HO 分别为在 A, B, C, D 及 O 各点的高程，O´为 A, B, C, D 四点的平均高程位置，这样，我们将O点与O´点的高差作为该栅格的地形描述误差Et，即： Et = HO - HO’ = HO - (HA + HB + HC + HD) / 4 (1) 按照公式(1)，已知单元栅格中点与其周围四个相邻点的高程是获取该单元 Et 值的必要条件。由于单元栅格中点的实际高程是未知的，在实际计算中采用窗口分析法实现 Et 的提取。对于 33 的正方形分析窗口，如果 DEM 空间分辨率为 d, 在该分析窗口内，i 行 j 列 Et 值可以通过下式求得: Et(i, j) = H(i, j) - (H(i-1, j-1)+H(i-1, j+1)+H(i+1, j-1)+H(i+1, j+1) )/4 (2) 显然，33 窗口的分析分辨率为 2d（见图 14-3）。当顺序移动该分析窗口对整个 DEM 进行逐点 Et 值计算，即能获得相应的误差矩阵。当将分析窗口依次扩大到 55、77.，可以同理提取分析分辨率分别等于 4d、6d、. 的误差矩阵。对于 5m 分辨率的 DEM，其相应的分析分辨率依次为 10m、20m、30m、.100m。四、DEM 地形描述误差的量化模拟根据上述方法所得到的不同分辨率条件下的误差矩阵，便可以利用统计与比较分析的方法揭示 Et 随分辨率及地形复杂度的变化而变化的规律，表 14-2 及图 14-4 为对误差矩阵的数据进行统计的结图 2. DEM 高程采样栅格元 Fig. 14-2 A DEM grid unit model (i, j) (i-1, j-1) (i-1, j+1) (i+1, j-1) (i+1, j+1) (i, j) (i-2, j-2) (i+2, j+2) (i-2, j+2) (i+2, j-2) 2d 4d 图 14-3 栅格分析窗口示意图

果。表 14-2 不同地貌类型区 Et 均方差值（RMSE）统计表分辨率（m）平原低丘丘陵中山高山混合 10 0.599 0.678 0.856 1.012 1.378 0.938 20 0.975 1.237 1.831 2.350 3.571 2.102 30 1.350 1.796 2.805 3.687 5.763 3.266 40 1.726 2.355 3.779 5.025 7.955 4.431 50 2.101 2.914 4.754 6.363 10.148 5.595 60 2.476 3.474 5.728 7.701 12.340 6.759 70 2.852 4.033 6.703 9.039 14.533 7.924 80 3.227 4.592 7.677 10.376 16.725 9.088 90 3.602 5.151 8.651 11.714 18.917 10.252 100 3.978 5.710 9.626 13.052 21.110 11.417 DEM 的均方差值（RMSE）被是描述 DEM 误差的重要统计指标[10, 11]。图 4 显示在不同实验区内， Et 的均方差值（RMSE）随分辨率的降低而升高并呈很好的线性相关关系。整理图 5 中的回归方程，得到式（4）。如果将以上方程视为 Y=aX+b 模式，可以发现以上方程的系数 a，b 也分别同试验样区地面平均剖面曲率呈较好的线形相关（见图 14-5、图 14-6）。综合以上结果，可以进一步得出以下方程： RMS Et =（0.0063 V + 0.0066）R - 0.022 V + 0.2415 (5) 其中，R、V 分别表示 DEM 的空间分辨率与平均剖面曲率。公式（5）显示 Et 同 DEM 分辨率与反映地形复杂度的因子平均剖面曲率成正相关；混合地貌类型实验区同样适合运用以上误差模拟方程。如果将公式改写为： (高山) (中山) (丘陵) (4) (低丘) (平原) (混合地形) RMS Et = Y=0.2139X - 0.4278 Y=0.1334X - 0.3101 Y=0.0999X - 0.2253 Y=0.0551X+0.1094 Y=0.0378X+0.2063 Y=0.1181X - 0.2770 高山：y = 0.2139x - 0.4278 r 2 = 0.9956* 中山：y = 0.1334x - 0.3101 r 2 = 0.9972* 低丘：y = 0.0551x + 0.1094 r 2 = 0.9959 平原：y = 0.0378x + 0.2063 r 2 = 0.9916 混合地形：y = 0.1181x - 0.277 r 2 = 0.9995 丘陵：y = 0.0999x - 0.2253 r 2 = 0.9975* 0 5 10 15 20 0 10 20 30 40 50 60 70 80 90 100 DEM分辨率 RMS Et 平原低丘丘陵中山高山混合地形 * The correlation is significant at 0.01 level 图 14-4 Et 随 DEM 分辨率及地形复杂度变化回归模型 ( r = 相关系数)

R=(RMSEt+0.022V-0.2415)/(0.0063V+0.0066)(6)0. 30. 3a= 0.0061V+0.0052b=-0.022v+0.24150.20. 22=0.9899ra=0.84250. 1 X0.01030-0. 1200. 2 0. 1-0. 30. 0-0. 4 4001020300. 5平均剖面曲率V平均剖面曲率V图14-5方程系数a与面曲率相关性图14-6方程系数b与剖面曲率相关性(r为相关系数）(r为相关系数）可根据DEM误差的限定指标直接推算适宜的DEM分辨率。其中剖面曲率为地面高程变化的2阶导数，其实质为栅格分析窗口内，坡度在水平方向的最大变化率。Tang[10]已经证明，在ARC/INFO、ArcView等GIS软件中，地面剖面曲率的数字矩阵可以直接通过对DEM数据求取slopeofslope（地面坡度的坡度）而获得：亦可经过统计计算获得研究区的平均剖面曲率。选择陕西汉中老君（平原区）、铜川高坪（黄土丘陵区）和潼关太要（山区）为检验样区（4kmx4km），采用高精度航测的方法提取DEM（1m分辨率）并测定其地形描述误差值，通过与模拟方程结果的对比，显示该误差模拟方程对3个检验样区Et估算的准确性分别为91.3%、94.7%和95.3%，显示该方程以上具有良好的误差估算效果。五、误差地图虽然对于某一区域DEM误差的统计值有助于总体上了解不同地貌类型区的误差量值，但在任何试验区内，DEM误差的空间分布是否存在有特定的规律，是否随着地形部位的不同而变化，是普遍关心的问题。Guth（1992)、Li(1993)及Monckton(1994)曾强调误差在空间分布具有一定的结构化特征。误差地图无疑是最能直观反映误差空间分布规律的技术手段。Monckton（1994)曾探索利用离散点位的专题制图法绘制DEM的误差地图，但是，由于离散点位在描述连续现象的局限性，误差地图的效果不甚理想。本研究采用窗口分析法获取了在空间上连续的误差矩阵，因而能够在ARC/VIEW地理信息系统软件的支持下，利用质地法绘制误差地图。图14-7为黄土丘陵沟区的DEM地形描述误差地图（DEM分辨率为25米）。通过与图中等高线的对比分析，反映误差随地形变化的基本规律。DEM地形描述误差在空间分布呈较为明显的自相关性，较大误差主要分布在山谷、山脊以及地面坡度转折处：图中反映出误差值的大小在很大的程度上同地面垂直曲率、水平曲率有密切关系。在今后的研究中应对其相关关系进行定量测算，从而为误差模拟方程的建立提供基本依据。误差地图能直观地提供误差在空间分布的规律与特征，而量化DEM地形描述误差在空间的分布规律，对02m2-6m等高线图例于以后所进行的误差的微观模拟，是必不可少的前提6-10m>10m条件。而自相关分析法是其中十分有效的方法（Wood,图14-7DEM误差地图（25米分辨率）

R= (RMS Et+ 0.022 V - 0.2415) /（0.0063 V + 0.0066） (6) 可根据 DEM 误差的限定指标直接推算适宜的 DEM 分辨率。其中剖面曲率为地面高程变化的 2 阶导数，其实质为栅格分析窗口内，坡度在水平方向的最大变化率。Tang [10] 已经证明，在 ARC/INFO、ArcView 等 GIS 软件中，地面剖面曲率的数字矩阵可以直接通过对 DEM 数据求取 slope of slope（地面坡度的坡度）而获得；亦可经过统计计算获得研究区的平均剖面曲率。选择陕西汉中老君（平原区）、铜川高坪（黄土丘陵区）和潼关太要（山区）为检验样区（4km4km），采用高精度航测的方法提取DEM (1m分辨率)并测定其地形描述误差值，通过与模拟方程结果的对比，显示该误差模拟方程对 3 个检验样区 Et 估算的准确性分别为 91.3%、94.7%和 95.3%, 显示该方程以上具有良好的误差估算效果。五、误差地图虽然对于某一区域 DEM 误差的统计值有助于总体上了解不同地貌类型区的误差量值，但在任何试验区内，DEM 误差的空间分布是否存在有特定的规律，是否随着地形部位的不同而变化，是普遍关心的问题。Guth (1992)、Li (1993)及 Monckton (1994)曾强调误差在空间分布具有一定的结构化特征。误差地图无疑是最能直观反映误差空间分布规律的技术手段。Monckton (1994)曾探索利用离散点位的专题制图法绘制 DEM 的误差地图，但是，由于离散点位在描述连续现象的局限性，误差地图的效果不甚理想。本研究采用窗口分析法获取了在空间上连续的误差矩阵，因而能够在 ARC/VIEW 地理信息系统软件的支持下，利用质地法绘制误差地图。图 14-7 为黄土丘陵沟壑区的 DEM 地形描述误差地图（DEM 分辨率为 25 米）。通过与图中等高线的对比分析，反映误差随地形变化的基本规律。DEM 地形描述误差在空间分布呈较为明显的自相关性，较大误差主要分布在山谷、山脊以及地面坡度转折处；图中反映出误差值的大小在很大的程度上同地面垂直曲率、水平曲率有密切关系。在今后的研究中应对其相关关系进行定量测算，从而为误差模拟方程的建立提供基本依据。误差地图能直观地提供误差在空间分布的规律与特征，而量化 DEM 地形描述误差在空间的分布规律，对于以后所进行的误差的微观模拟，是必不可少的前提条件。而自相关分析法是其中十分有效的方法（Wood, 图例 2-6m 6-10m >10m 等高线图 14-7 DEM 误差地图（25 米分辨率） a = 0.0061V + 0.0052 r 2 = 0.9899 0.0 0.1 0.1 0.2 0.2 0.3 0 10 20 30 40 平均剖面曲率 V 方程系数 a 图 14-5 方程系数 a 与剖面曲率相关性（r 为相关系数） b = -0.022v + 0.2415 r 2 = 0.8425 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0 10 20 30 40 平均剖面曲率 V 方程系数 b 图 14-6 方程系数 b 与剖面曲率相关性（r 为相关系数）

1996）。地理学的一个基本规律既是：空间的事物总在不同程度上相互联系与制约，而相近的事物之间的影响通常大于较远事物的影响。在地形学的研究中，地形的自相关值往往被用以描述地面粗糙程度 (Tobler, W. 1979)。DEM 误差的空间自相关可以定义为某一栅格单元的误差值与其相邻栅格误差值的趋近程度。通过计算 DEM 地形描述误差的空间自相关值，可借以判断 DEM 误差在水平方向的集聚度，从而进一步反映 DEM 误差在空间的分布特征。空间自相关通常采用 Moran 统计算子求算，Moran 自相关算式可以表达为：    = = = = = =  − n i n i n j i ij n i n j ij ij v w w c I 1 1 1 1 1 ( ) .(7) 其中 wij 为给予每一个栅格测量单元的权重，cij 和 vI 的量值可以表达为： c (z z)(z z) ij = i − j − . . (8) n z z z z v i i i ( − )( − ) = . . (9) 其中 zi 、zj 分别为在位置 i、j 的测量值，z 是在所有 i、j 位置点测量值的均值，n 为所有测量点的数目。自相关 I 的值域在+1、-1 之间。数据越接近+1，表示正自相关愈强；越接近-1，表示负自相关愈强，0 表示非自相关随机分布。以上公式在空间属性分布研究中得到广泛的应用。然而，Wood (1996) 的研究工作证明，对于栅格数据自相关的计算，自相关计算公式可以简化为：     = = = = = − − − n i i j i j n j n i n j j i j i j i z z z z z z I 1 2 1 1 1 ( ) ( )( )     .(10) 根据以上原理，我们以所获得的 DEM 误差数字矩阵为基本数据源，计算 DEM 误差在不同地形样区及不同空间分辨率条件下的误差自相关值。表 14-3 及图 14-8 为试验结果。表 14-3. 误差空间自相关量测结果地貌类型 DEM 空间分辨率 m 5 10 15 20 25 30 平原 0.120 0.286 0.415 0.482 0.520 0.527 低丘 0.180 0.360 0.478 0.522 0.548 0.54 丘陵 0.201 0.350 0.534 0.531 0.513 0.37m 中山 0.262 0.386 0.527 0.548 0.510 0.409 高山 0.351 0.492 0.564 0.562 0.426 0.344