采用ArcGIS平台的地势起伏度自动提取技术研究

王 康，何俊仕，于德浩，龙 凡，李 霞

(1.沈阳农业大学 水利学院，辽宁 沈阳 110866;2.沈阳军区司令部工程科研设计院，辽宁 沈阳 110162)

随着3S技术的快速发展，利用DEM(digital elevation model)数字高程模型和遥感影像等数据源进行地貌特征信息的提取已成为地貌研究的主要手段。其核心思想是根据相对起伏度或海拔高度来快速地确定基本地貌形态类型［1－4］。从理论上分析，只要DEM数据足够精确，就可以很准确地确定基本地貌形态。因此对于地势起伏度的提取，关键在于一定统计窗口的确定。即确定面积有多大，才能恰到好处地反映山体的完整性，并具有一定范围内较强的代表性，亦即普适性［5］。刘振东［6］提出按照地貌发育的基本理论，存在一个使最大高差达到相对稳定的最佳统计窗口，这个最佳统计窗口即为地势起伏度的最佳统计单元。一般采取“拐点”理论［5］来确定窗口大小，这里所定义的“拐点”与数学上的拐点不同，这里“拐点”的定义为，高差由大变小的点为“拐点”。做统计时会发现这样一条规律:那就是随着统计窗口的扩大，窗口中最高点与最低点的高差不断增加，开始增加的很快，当高差增加到一定值时，随着统计窗口的再增加高程差增加不大，而这样的点，则为要寻找的“拐点”。这个“拐点”所对应的面积，即为最终的统计单元，而“拐点”所对应的高差值即为地势起伏高度。

利用ArcGIS软件，加载Spatial Analyst模块，借助Neighborhood Statistic功能，对工作区以3×3起始窗口，到90×90终止窗口，移动步距多为1，计算窗口内最大高程值MAXIMUM和最小高程值MINIMUM［7］。然后利用 Raster Calculator工具计算专题层的MAXIMUM与MINIMUM的差值即为地势起伏度层，其每个栅格的值即为以该栅格为中心的确定领域的地面起伏度值。最后，利用窗口大小与最大高程差值生成所需图表文件，并依此寻求拐点位置。

上述方法的确可行，但是此方法需依次重复从3×3窗口到90×90窗口的领域计算等操作，该过程工作繁琐，效率低，中间生成大量缓存文件，而且大量的繁琐操作容易出现人为操作失误，导致最终数据精度不可靠。

于德浩等［8］学者曾利用遥感特征指数对地表水体实现了自动化提取。基于此成功案例，本文针对常规方法提取地势起伏度操作繁琐、效率低，以及提取过程生成大量缓存文件等问题，提出一种自动提取的方法。该方法对传统“拐点”理论手工确定窗口大小的过程加以改进，利用计算机编程实现不同统计单元起伏度值的自动计算，并依此寻求最佳统计单元，实现了地势起伏度的自动提取。

1 地势起伏度及其作用

地势起伏度(Relief Amplitude)，也称为地形起伏度或地势能量、相对高度。在我国国土大地貌研究和1:400万地貌图工作中，将地势起伏度定义为:一个山体的山脊最高点与顺水流方向最近河流的第一个汇流处为基准面的最大高差［9］。在实际操作中，涂汉明等［10］又将某点的地势起伏度RAi定义为某一确定面积(A)中的最高点(Zimax)与最低点(Zimin)之高差。它是描述一个区域地形宏观性指标之一，在宏观的区域内反映地面的起伏特征。用公式表示如下:

式中:RAi指某区域的地势起伏度;Zimax指该区域内的最大高程值;Zimin则指该区域内的最小高程值;i为自然数，指某一统计区域。

坡度和起伏度是地形描述中最常用的参数，它们能快速而直观地反映地势起伏特征。其中，坡度是划分平原和非平原的重要依据之一;而地势起伏度又可进一步划分为平原、台地、丘陵、小起伏山地、中起伏山地、大起伏山地和极大起伏山地等7级［11］。地形起伏是构造作用与地表剥蚀过程相互作用的结果，是地貌研究的必要手段［1］。地势起伏度的提取为获取地表信息，进行地形定量化分析提供有效手段。

2 地势起伏度自动提取的实现

2.1 自动计算3×3窗口到90×90窗口的最大值与最小值的差值

应用ArcGIS的Python模块进行批处理，输入代码:

在工作空间目录下新建一个文件夹命名为Output，用来保存输出结果。

运行结束之后会在output文件夹下面自动生成r_3至r_90的栅格数据集。操作界面如图1所示。

图1 Python操作界面

2.2 读取所生成数据集的起伏度值(最大值属性)

在事先准备好的Excel中A1输入“r_3”;B1输入“MAXIMUM”。通过下拉生成88行如图2所示。

将表格内容复制，然后应用ArcGIS10在Arc-Toolbox中选择“获取栅格属性”工具，选择“批处理”。添加行，添加至88行。将88行全部选中然后粘贴可简化对设置每行的数据输入和所获取属性等的复杂操作，之后在“环境…”中把当前工作空间和临时工作空间都设置为:E:/exe/output。

图2 Excel所生成内容

确定执行后返回的MAXIMUM属性将显示在地理处理移动窗口中，见图3。

图3 消息窗口

右击“消息”，将其复制，粘贴在Word中。

将第一行“消息”、第三行“开始时间……”和最后一行“成功在……”等内容删掉。将“执行:GetRasterProperties E:/exe/output/”、 “ MAXIMUM”和“最大值 =”均替换为空，之后转换成2列的表格。粘贴到Excel中，便可根据要求生成图表来寻求“拐点”。

提取结果见图4、图5。从图中均可看到在一定的范围内随着移动窗口的扩大，起伏度也在不断扩大，但当窗口达到一定程度后起伏度增长已十分缓慢，之后随着移动窗口的增大起伏度已趋于稳定。这里的“拐点”即为所要寻求的最佳点——地势起伏度。

通过对测试区提取结果进行分类，划出地势起伏度分级图，见图6a，较符合实际情况;同时，对局部区域放大并结合等高线分析，见图6b，地势起伏度较大地区集中在等高线密集的山区地带，而地势起伏度较小的地区主要集中在等高线稀疏的河谷等平缓地带，与实际情况基本吻合。

图4 试验区I最佳统计单元提取结果

图5 试验区II最佳统计单元提取结果

图6 地势起伏度提取结果

3 结论

本文提出了一种地势起伏度自动提取的方法，对传统的利用“拐点”理论确定最佳统计单元的过程加以改进，通过编程实现了不同统计单元起伏度值的自动计算，并依此确定最佳统计单元，实现了地势起伏度的自动提取。与常规方法比较，此方法具有如下优点:

(1)操作简便，效率高:常规方法需要重复从3×3窗口到90×90窗口的邻域计算等操作，操作复杂，工作量大，本方法通过Python模块编程能够实现此过程的自动化运行，减少了繁琐而复杂的操作，效率得到了明显提高。

(2)减少大量不必要缓存文件:常规方法计算中，首先求算分析窗口内DEM数据的最大与最小高程值，分别记作专题层Max与Min，之后利用Raster Calculater工具计算专题层Max与Min的差值，得到地势起伏度专题层。这意味着将要生成三倍于原始文件大小的缓存文件，如果重复计算3×3窗口到90×90窗口的邻域计算则需要生成264倍于原始文件大小的缓存文件。本研究一步到位，能够直接计算出起伏度值专题层，可使缓存文件缩小到常规方法的1/3。

该研究有利于地势起伏度研究工作的快速开展，可为该研究提供有效的技术支撑，同时也必将促进地势起伏度研究工作的快速发展。此外，由于地势起伏度值计算过程均是计算机自动计算，因此难免出现许多面积相对较小、很难在图上得以显示的破碎图斑，这也是以后工作中有待提高的方面。

［1］A Kühni，O.A Pfiffener.The relief of the Swiss Alps and adjacent areas and structure:topographic analysis from a250-m DEM［J］.Geomorphology，2001，41:285－307.

［2］G Onorati，R Ventura，M Poscolieri，et al.The Digital Elevation Model of Italy for geomorphology and structural geology［J］.CATENA，1992，19(2):147 －178.

［3］闾国年，钱亚东，陈钟明.基于栅格数字高程模型提取特征地貌技术研究［J］.地理学报，1998，53(6):562－569.

［4］Balázs，Székely，Dávid Karátson.DEM-based morphometry as a tool for reconstructing primary volcanic landforms:examples from the Borzsony Mountains，Hungary［J］.Geomorphology，2004，63:25 －37.

［5］涂汉明，刘振东.中国地势起伏度最佳统计单元的求证［J］.湖北大学学报(自然科学版)，1990，12(3):266－271.

［6］刘振东.利用DTM编制小比例尺地势起伏度图的初步研究［J］.测绘学报，1990，19(1):57－62.

［7］中国科学院地理科学与资源研究所资源与环境信息系统国家重点实验室.中国1∶100万数字地貌遥感解译与集成流程［S］.(试行本).2005.

［8］于德浩，龙凡，邓正栋，等.基于遥感特征指数的地表水体自动提取技术研究［J］.地球物理学进展，2009(2):707－713.

［9］刘振东，涂汉明.中国地势起伏度统计单元的初步研究［J］.热带地理，1989，9(1):31 －38.

［10］涂汉明，刘振明.中国地势起伏度研究［J］.测绘学报，1991，20(4):311 －319.

［11］中科院地理科学与资源研究所(资源与环境信息系统国家重点实验室).中华人民共和国1:100万数字地貌图制图规范［Z］.(征求意见稿).2005.