 Geomatics Science and Technology 测绘科学技术, 2013, 1, 11-17 http://dx.doi.org/10.12677/GST.2013.13003 Published Online October 2013 (http://www.hanspub.org/journal/gst.html) The Study on Terrain Index Extraction Method and Practice Based on ArcGIS Model Builder* Yan Xiao1, Peng Xie2 1School of Environment and Planning, Liaocheng University, Li aocheng 2College of Urban and Environment, Liaoning Normal University, Dalian Email: xy790511@163.com Received: Sep. 2nd, 2013; revised: Sep. 10th, 2013; accepted: Sep. 15th, 2013 Copyright © 2013 Yan Xiao, Peng Xie. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use , distribution, and re production in any medium, provided the original work is properly cit ed. Abstract: River terrain index TOPMODEL hydrological model has been applied to Land surface model (Land Surface Models, LSMs), to improve Land surface of hydrological process simulation model. The core of TOPMODEL is terrain index (ln(α/tanβ)). And it is an important index on the hydrologic analysis. Therefore, effective extracting terrain index information is very important. In this paper, the huaihe river basin as a case area, the terrain index extraction process model is build, using the ArcGIS9 Model Builder tools. This model is based on a single flow algorithm for calculating process. Model has higher versatility, because of fully considering the effectiveness of the parameters such as the slope and the flow accumulation in the model design process. Using this model can be very simple and quick to get terrain exponential di st ri but i o n dat a of terrain index. Keywords: Terrain Index; ArcGIS9; Model Builder; Single Flow Algorithm 基于 ArcGIS Model Builder的地形指数提取方法及实践研究* 肖 燕1，解 鹏2 1聊城大学环境与规划学院，聊城 2辽宁师范大学城市与环境学院，大连 Email: xy790511@163.com 收稿日期：2013 年9月2日；修回日期：2013 年9月10 日；录用日期：2013 年9月15 日 摘 要：流域地形指数水文模型 TOPMODEL 已被应用于陆面模式(Land Surface Models, LSMs)以改进陆面模式 对水文过程的模拟。地形指数(ln(α/tanβ))是TOPMODEL 的核心，是进行水文分析的重要指标。因此，高效的提 取地形指数信息非常重要。本文将淮河流域作为案例区，以地形指数单流向算法为理论基础，利用 ArcGIS9 的 Model Builder建模工具，建立地形指数提取模型。在模型设计过程中充分考虑到汇流面积、坡度等参数的有效 性，使模型具有较高的通用性。使用者启动模型加载数据即可简单又快速地获得地形指数分布数据。 关键词：地形指数；ArcGIS9；Model Builder；单流向算法 1. 引言 1979 年Beven 和Kirky 提出以地形为基础的半分 布式流域水文模型 TOPMODEL (Topography-based hydrological model)，充分考虑了地形对产流区形成和 变化的影响，以地形指数 ln(α/tanβ)的空间格局来反映 流域饱和缺水量的空间分布[1]。其中 α为流经坡面任 一点单位等高线长度的汇流面积，β为该点处的坡度。 此模型及其理论在水文学中一直被广泛应用。地形指 *本研究得到国家自然科学基金面上项目(41075060)和山东省自然 科学基金(ZR2011DQ015)的支持。 Copyright © 2013 Hanspub 11  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 数作为该模型的重要参数，其计算和提取非常关键。 因此，正确而又高效的提取地形指数信息对于水文模 型的发展和应用尤其重要。国内外学者对于地形指数 计算的研究也有不少。雍斌、张万昌、陈艳华[2]，基 于多流向算法对流动累积分配中的有效等高线长度 精确计算提出了几何锥面内切圆算法，改进了传统的 地形指数中单位等高线汇流面积α的计算方法，增强 了多流向算法对 DEM 中异常栅格的处理能力。孙凡 哲、芮孝芳[3]对地形指数计算方法做出改进，计算下 坡单元累积汇流面积时包括了欲计算 ln(α/tanβ)的单 元网格的面积。但这些研究大都建立在专门的程序算 法基础上。对于应用人员的专业素质要求比较高，且 可视化效果也不好。本研究以地形指数的生成和步骤 为基础，探索使用 ArcGIS 软件地理处理 (Geoprocessing)的可视化建模平台(Model Builder)进 行地形指数提取。并通过本模型的建立，以期有利于 模型推广和应用。 2. 理论基础 TOPMODEL 模型中的地形指数 ln(α/tanβ)计算方 法大体有两种：多流向算法和单流向算法。一般使用 基础数据都是数字高程模型。ArcGIS9 中流域水文分 析工具内嵌算法是基于单流向的D8 算法。因此本研 究地形指数的提取算法基础为单流向算法。 2.1. 单流向算法 目前已广泛使用 DEM 数据作为计算地形指数空 间分布的基础数据。计算地形指数 ln(α/tanβ)，要 分 别 计算流经坡面任一点单位等高线长度的汇流面积 α， 及该点处的坡度 β[4]。 地形指数单流向算法是由 O Callaghan和Mark 于 1984 年提出的，此方法计算简单，通常认为单位等高 线长度等于单元格网长度，坡度为顺坡方向最大坡 度。其中 α定义为 1N A (1) 式中：N——流入一个网格单元的上坡汇流栅格数； A——单个栅格面积[5]。坡度定义为 tan h l (2) 式中： ——相邻栅格间的高程差； ——相邻栅 格中心的距离，主方向 hl l 等于单元网格长度，对角线 方向 l 等于 2×单元格网长度[6]。 2.2. ArcGIS建模介绍 2.2.1. 地理处理 地理处理的基本目的在于使用户能够自动执行 GIS 任务以及执行空间分析和建模任务。几乎 GIS 的 所有使用情况都会涉及重复的工作，因此需要创建可 自动执行、记录及共享多步骤过程(即工作流)的方法。 地理处理通过提供一组丰富的工具和机制来实现工 作流的自动化操作，这些工具和机制能够使用模型和 脚本将一系列的工具和一组按顺序的操作结合在一 起。 2.2.2. ModelB ui l der建模平台 模型主要有三部分组成，输入数据、空间处理工 具和输出数据。输入数据和输出数据类型多样，可以 为矢量数据也可为栅格数据。空间处理工具是模型的 核心，参与建模的空间处理工具可以为系统默认的工 具也可为用户自定义的工具模型。将上述三个主要元 素通过连接工具连接起来，就构成了模型。 在ArcGIS 软件的 ArcToolbox 窗口内提供许多系 统工具箱，包括数据管理、数据转换、空间分析、地 理编码等等许多复杂的空间处理工具。这些功能都是 固定用户无法修改。创建地形指数提取模型之前应当 首先创建地形指数工具集。建立模型后可以进行编 辑。在拖拽 ArcToolbox 中的模型至当前模型编辑窗 口，就可以使用相应的功能。 3. 应用案例 3.1. 案例区介绍 淮河流域地处我国东部，介于长江和黄河两流域 之间，位于 111˚55'~121˚25'E，33˚55'~36˚36'N，面积 为27 万km2。淮河流域西部、西南部及东北部为山区、 丘陵区，其余为广阔的平原。山丘区面积约占总面积 的三分之一，平原面积约占总面积的三分之二。流域 西部的伏牛山、桐柏山区，一般高程 200~500 m (85 黄海高程，下同)；南部大别山区高程在300~1774 m； 东北部沂蒙山区高程在200~1155 m。丘陵区主要分布 在山区的延伸部分，西部高程一般为 100~200 m，南 部高程为 50~100 m，东北部高程一般在 100 m左右。 Copyright © 2013 Hanspub 12  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 淮河干流以北为广大冲、洪积平原，地面自西北向东 南倾斜，高程一般 15~50 m；淮河下游苏北平原高程 为2~10 m；南四湖湖西为黄泛平原，高程为 30~50 m。 流域内除山区、丘陵和平原外，还有为数众多、星罗 棋布的湖泊、洼地。 3.2. 数据源介绍 SRTM (Shuttle Radar Topography Mission)，由美 国太空总署(NASA)和国防部国家测绘局(NIMA)联合 测量。2000 年2月11 日，美国发射的“奋进”号航 天飞机上搭载 SRTM 系统，共计进行了 222 小时 23 分钟的数据采集工作，获取北纬 60 度至南纬 60 度之 间总面积超过 1.19亿平方公里的雷达影像数据，覆盖 地球 80%以上的陆地表面。SRTM 系统获取的雷达影 像的数据量约 9.8 万亿字节，经过两年多的数据处理， 制作成数字地形高程模型(DEM)。本研究使用 SRTMDEM V4.1版本。此数据该版本是由 CIAT (国际 热带农业中心)利用新的插值算法得到的 SRTM地形 数据，更好的填补了 SRTM 90的数据空洞。 3.3. 建模过程 本建模过程主要采用ArcGIS9的水文分析工具集 和空间分析工具集。基于单流向算法进行案例区地形 指数的计算。首先为了保证数据的完整性，在空间分 析模块中对数据进行空缺内插单元的内插填补处理， 然后分别进行坡度分析和水文分析。其中为了保证地 形指数 ln(α/tanβ)的有效性，将坡度小于 1˚的栅格赋值 1˚以此去除掉坡度为零的区域；将流水累积量为零的 区域赋值为 1。通过设置河道起始临界值 CIT (Channel Initiation Threshold, CIT)生成河道，集水面积大于 CIT 的网格作为河道网格，不参与地形指数分布函数的统 计。最后使用 ArcGIS9 空间分析工具中的栅格计算工 具进行地形指数的提取。 3.3.1 DEM填洼 在填洼处理过程中，闭合洼地和平坦区域均被认 为是伪地形，在提取地形指数之前必须进行移除。 Jensonhe 和Domingue 算法[7]及Martz 和Jong 算法[8] 都有自己的假设和相应的处理方法。他们的共同思路 都是采用填平，然后逐步垫高最终使流域内每个网格 水流都能够流向流域出口。通过参考案例区各种水文 相关资料，确定填充最低阈值，通过 ArcGIS9 软件中 ArcToolbox 中的 Fill 工具进行填充。 3.3.2 计算汇流累积量 汇流累积量数值矩阵表示的是区域地形每点的 流水累积量的空间分布，其汇流累积量的数值等于其 上游最终汇流经过该单元格的栅格总数。地表径流漫 流模型中计算汇流累积量的基本思想是：以规则格网 表示的数字地面高程模型每点处有一个单位的水量 基于水流方向数字矩阵计算每点处所流过的水量数 值，便可以得到该区域的汇流累积量数值矩阵[9]。此 过程用到 Flow Direction与Flow Accumulation两个工 具，将其加入到计算汇流累计量模型中。见图 1。 3.3.3. 计算有效汇流面积 对于含河道的网格由于其汇流面积包含了河岸 两边坡面的汇流面积以及该网格以上的汇流面积，如 不将河道网格进行一定的处理将会造成流域的平均 地形指数以及高地形指数的网格数的增加，因此使地 形指数分布函数发生变化。若消除河道网格对地形指 数分布函数的影响，可设置汇流累积量 a的上限值。 对于汇流面积大于该上限值的网格，在计算地形指数 时，汇流面积取该上限值。该上限值就是生成河道的 起始临界值(Channel Initiation Threshold, CIT)，汇流面 积大于该上限值的网格被定义为河道[10]。 根据以上理论建立上限控制操作模型，用于控制 栅格河网数据的上限值，如图 2。采用 DEM 数据为 90 m × 90 m数据，单个栅格单元面积为 8,100 m2，CIT 按5 km2计算得到上坡汇流累计栅格数上限值为 617。 以汇流累积量数值矩阵为基础数据，使用 Minus 工具将此汇流累积量数值矩阵为与上坡汇流累计栅 格数上限值作差，得到差额数据层，其中超额区域为 正值，未超额区域为负值。使用Greater Than工具， 将此汇流累积量数值矩阵与上坡汇流累计栅格数上 限值作逻辑大于运算，得到超额判定图层，其中超额 区域为 1，未超额的区域为 0。在使用 Times 工具将 上述两步结果作乘运算，得到超额数据图层，其中超 额区域为正值，且此值即为上坡汇流累计栅格数大于 上限值的差值，未超额区域为0。最后利用 Minus 工 具，将汇流累积量数值矩阵减去超额数据图层，即得 到上限值处理后的有效汇流累积量数值矩阵。 Copyright © 2013 Hanspub 13  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 Copyright © 2013 Hanspub 14 Figure 1. Confluence accumulation model 图1. 汇流累计量模型 Figure 2. River limit control model 图2. 河网上限控制模型 波，张万昌[12]的论文中提到对坡度为 0的点用该点毗 邻的坡度非零栅格中的坡度值最小的那个来替代当 前栅格的坡度值。本模型构造中根据实际情况，将小 于1˚的坡度值替换为1˚。提取有效坡度的模型如图 4 所示。 为了得到有效汇流面积还需要将栅格河网数据 转化为汇流面积数据，其处理过程如图 3。以汇流累 积量数值矩阵为基础数据，使用 Plus 工具将此数据做 加1处理，得到加1后的汇流累积量数值矩阵，再利 用上述上限控制操作模型得到有效汇流累积量数值 矩阵，最后利用 Times 工具将此有效汇流累积量数值 矩阵乘以单个栅格面积，即得到有效汇流面积。 以坡度图层为基础数据，使用 Greater Than Equal 工具、Less Than Equal工具获得分别获得大于 0与小 于1的部分，然后使用 Boolean And工具得到坡度介 于0~1 的二值数据图层。对该数据做取反运算，再与 原坡度图层相乘获得抠图数据图层。同时该二值数据 图层使用 Times 工具计算获得新值区域。最后将新值 区域与抠图区域做 Plus 运算，获得有效坡度图层。 3.3.4. 坡度提取 地表面任一点的坡度是指过该点的切平面与水 平平面的夹角，坡度表示了该表面在该点的倾斜程度 [11]。在 ArcGIS9 软件中可直接使用Arc Toolbox中的 Slope 工具提取坡度。 3.3.6. 计算地形指数 3.3.5. 计算有效坡度 得到有效汇流面积与有效坡度后，便可以根据公 式ln(α/tanβ)进行地形指数的计算。首先使用 Times 工 具将有效坡度与 在计算ln(α/tanβ)时不允许出现坡度值为0的情 况。通常的处理方法是把这个点去掉不算，认为此点 对产流贡献不大。然而事实上平坦的地区也是有产流 的，而且这样的处理方法也增加了模型的误差。在吴 180做乘法运算，使坡度由原来 的角度制转化为弧度制，以方便三角函数的使用。将  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 Figure 3. Calculating effective confluence area model 图3. 计算有效汇流面积模型 Figure 4. Calculating effective slope model 图4. 计算有效坡度模型 此弧度制的坡度图层做Tan 运算，在使用 Divide 工具 将有效汇流面积与 Tan 三角函数工具计算结果做除法 运算，最后使用 Ln 函数工具进行运算，得到最终的 地形指数。得到如图 5所示运算模型。 上述为计算地形指数各个子过程模型，分步使用 可得到处理过程的中间数据，方便查错修改。同时为 Copyright © 2013 Hanspub 15  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 简化中间环节操作，方便用户使用，笔者将以上个子 运算模块整合建立完整的地形指数提取工具模型如 图6。 3.3.7. 模型运算结果 加载地形指数提取工具模型，得到地形指数分布 栅格数据如图 7所示将数据做取整处理得到地形指数 分布数据，利用属性表中的导出数据功能，将地形指 数分布数据导出为 dbf格式，再使用Excel加载数据， 绘制频率分布图如图 8所示。 从图 7可以看出西南区域地形指数较小，东北区 域地形指数较大。图中出现局部连续的地形指数较高 值，经过与水系图对比发现此区域为湖区，符合地形 指数的物理意义。从图 8可以看出地形指数频率分布 整体呈正态分布，分布范围为 9到20，频率最大值在 14 位置处，最大频率为 0.35。 4. 结语 使用 ArcGIS9 软件中的Model Builder工具建模， 通过构建合理的流程，以SRTMDEM 数据为源数据进 行淮河流域地形指数的提取，主要得到以下几点结 论。 1) 在ArcGIS9 软件中建立地形指数提取模型完 全可行，模型采用 ArcGIS9 自带建模工具制作，无需 使用开发语言进行编码，简单易用，且可以直接得到 Figure 5. ln(α/tanβ) calculation model 图5. ln(α/tanβ)计算模型 Figure 6. Terrain index extraction model 图6. 地形指数提取模型 Copyright © 2013 Hanspub 16  基于 ArcGIS Model Builder的地形指数提取方法及实践研究 Figure 7. Terrain exponential distribution 图7. 地形指数分布 Figure 8. Terrain index frequency distribution 图8. 地形指数频率分布图 可视化地形指数分布图，在很大程度上减少了科研人 员的工作强度。 2) 本研究的模型通用性强，模型使用时，需要根 据实际情况输入 DEM 数据的分辨率及汇流累积量栅 格上限数量，可以适用于不同地形类型不同分辨率的 DEM 数据。 参考文献 (References) [1] K. Beven, M. J. Kirkby. A physically based, variable contribut- ing area model of basin hydrology. Hydrological Science Bulle- tin, 1979, 24(1): 43-69. [2] 雍斌, 张万昌, 陈艳华. TOPMODEL中地形指数 ln(α/tanβ)的 新算法[J]. 地理研究, 2007, 26(1): 37-46. [3] 孙凡哲, 芮孝芳. TOPMODEL 中地形指数计算方法的探讨 [J]. 水科学进展, 2003, 14(1): 41-45. [4] J. F. O. Callaghan, D. M. Mark. The extraction of drainage net- works from digital elevation data. Computer Vision, Graphics and Image Processing, 19 84, 28: 323-344. [5] 张渡光, 王克林, 陈洪松, 李新, 王书功. 基于 DEM 的地形 指数提取方法及应用[J].长江流域资源与环境, 2005, 14(6): 715-719. [6] 解河海, 黄国如. 地形 指数若干计算方法探讨[J]. 河海大学 学报(自然科学版), 2006, 34(1): 46-50. [7] S. K. Jenson, J. O. Domingue. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing, 1988, 54(11), 1593-1600. [8] L. W. Martz, E. de Jong. CATCH: A fortran program for meas- uring catchment area from digital elevation models. Computer & Geosciences, 1988, 14(5): 627-640. [9] 马兰艳, 周春平, 胡卓玮, 王志恒. 马国斌基于SRTM DEM 和ASTER GDEM的辽河流域河网提取研究[J]. 安徽农业科 学, 2011, 39(5): 2692-2695. [10] 邓慧萍, 李秀彬. 地形指数 的物理意义分析[J]. 地理科学进 展, 2002, 21(2): 103-110. [11] 汤国安, 杨昕. ArcGIS地理信息系统空间分析实验教程[M]. 北京: 科学出版社, 2006: 267-26 8. [12] 吴波, 张万昌. TOPMODE地形指数的计算[J]. 水土保持通报, 2006, 26(5): 49-52. Copyright © 2013 Hanspub 17 |