基于ASTER GDEM数据的青藏高原东部山区地形起伏度分析

张 学 儒 ，官 冬 杰，牟 凤 云，陈 春

基于ASTER GDEM数据的青藏高原东部山区地形起伏度分析

张 学 儒 ，官 冬 杰，牟 凤 云，陈 春

（重庆交通大学河海学院，重庆 400074）

以青藏高原东部山区为研究区，基于高空间分辨率的ASTER GDEM数据，通过AML语言编程调用ArcGIS中用于邻域分析的focal函数，计算不同邻域尺度单元下地形起伏度。研究表明：邻域尺度单元大小对地形起伏度计算至关重要，起伏度值先随邻域尺度单元面积增大而快速增大，当邻域尺度单元面积达一定阈值后，其增大速度开始减缓并趋于平稳，且在增速减慢过程中存在一明显拐点，即最佳邻域尺度单元。通过高差显著性变化检验法，确定最佳邻域尺度单元为5.0625 km2，据此制作地形起伏度分级图，发现研究区自西北向东南地形起伏度逐步增加，地势以中度起伏（200～500 m）为主。

地形起伏度;ASTER GDEM;邻域分析;青藏高原东部山区

0 引言

地形起伏度作为反映地表起伏程度的一个指标，在地貌制图、资源环境评价等领域有广泛应用，是地貌类型提取和划分的重要指标，也是区域地貌对比研究的客观依据［1］。近年来随着遥感、地理信息系统以及全球定位系统等技术的快速发展，基于数字高程模型（Digital Elevation Model，DEM）的地形起伏度计算研究工作得以开展。国际地理学会地貌调查制图委员会基于DEM数据采用16 km2的统计面积，计算地形起伏度并编制了欧洲1∶250万地貌图［2］。涂汉明等利用国家数字地形数据库，结合中国地貌类型的基本特征，提出地形起伏度最佳统计单元应满足山体完整性与区域普适性两条原则，论证了中国地形起伏度最佳统计单元为21 km2［3］;汤国安等基于DEM数据，使用窗口分析法计算多尺度地形起伏度［1］;刘新华等基于我国1∶100万DEM数据，在全国采集6个样区，运用窗口分析法得出我国地形起伏度最佳统计窗口为25 km2［4］;唐飞等通过对准噶尔盆地及其西北山区地形起伏度的研究，得出新疆克拉玛依地区地形起伏度的最佳统计单元面积为4 km2［5］;郎玲玲等比较了福建低山丘陵区不同尺度DEM提取的地形起伏度，认为基于1∶10万、1∶25万DEM计算地形起伏度的最佳统计单 元 分 别 为 0.41 km2、4.41 km2［6］;龙 恩 基 于SRTM-DEM数据对长白山区地形起伏度的最佳断点进行了判断［7］;张磊基于1∶25万DEM数据，采用邻域分析和Focal函数计算了多尺度地形起伏度，实现了京津冀地区地貌类型的自动划分［8］;王玲等采用变点分析法处理1∶25万DEM数据，实现了新疆地区地貌类型的自动划分［9］。传统方法由于地形数据量大、计算繁琐，使定量的地表形态、地貌研究发展缓慢，DEM的快速精确获取和GIS技术的发展为基于DEM进行地形起伏度提取与分析提供了有效手段。本文基于ASTER GDEM数据，探讨一种高效率计算不同邻域尺度单元下地形起伏度的方法。

1 研究区概况与数据

研究区地处青藏高原东南边缘与四川盆地西部边缘山地交接地带（北纬31°04′～35°64′、东经96°92′～104°40′），横跨青海、四川两省，涉及玛多、达日、班玛、色达、壤塘、马尔康、金川7县，面积74 865 km2，呈西北东南走向的条带状，属于高原与高山峡谷过渡地带［10］，也是我国第一个阶梯与第二个阶梯的过渡地区，按地貌分为高原区、山原区和高山峡谷区，地势自西北向东南逐渐降低。西北部高原区多为各水系上游切割地区，岭平谷宽，山体较大，山坡较缓，地形起伏度较小;东南部高山峡谷区多为各水系切割剧烈地区，山脊多呈锯齿状，地势险要，地形垂直差异明显［11，12］。研究区平均海拔 4 245 m，其中4 200～4 400 m的高程所占比重最高，最高峰是位于玛多县南部的巴颜喀拉山，海拔5 286 m［10］。

本研究采用ASTER GDEM数据，该数据是美国航空航天管理局（NASA）和日本经贸及工业部（METI）2009年7月2日共同发布的全球空间分辨率为30 m的最新栅格型数字地形数据，是根据NASA新一代对地观测卫星TERRA的详尽观测结果制作完成的，数据覆盖北纬83°至南纬83°之间的所有陆地区域。DEM的空间分辨率对地形起伏度的计算精度至关重要，在ASTER GDEM推出以前，可以获取到的大范围DEM数据是空间分辨率为90 m的SRTM DEM数据，该数据在小尺度地形提取时一直受到空间分辨率过低的影响。ASTER GDEM数据最大的优势在于其拥有较高的空间分辨率，在小尺度精细地形特征提取时明显优于SRTM DEM数据;另外ASTER GDEM数据空间分辨率与TM遥感影像一致，为两者的结合使用提供了便利条件。

2 地形起伏度计算

2.1 地形起伏度计算原理

地形起伏度（Relief Amplitude）是表征地表单元地势起伏复杂程度的一个指标，也是定量描述地貌形态、划分地貌类型的重要指标，其表达式如下［13］：

式中：RA表示地形起伏度，Hmax和Hmin分别是单位面积内最大和最小高程值。可以利用ArcGIS用于邻域分析的focal函数实现上式的计算，采用Focalmax和Focalmin分别计算单位面积内的最大高程值和最小高程值，二者相减即可求得地形起伏度。

2.2 邻域分析

邻域分析又称窗口分析，是一种局部运算的空间分析，涉及一个中心栅格单元和一组环绕栅格单元，根据中心单元和环绕单元的栅格数值为中心单元位置生成一个新的函数值或邻域统计量，即输出栅格图中每个格子的值均是输入栅格图中对应格子邻域单元的函数。在ArcGIS中基于邻域函数的统计方法有10种，在地形起伏度计算中仅涉及Maximum和Minimum两种邻域分析方法。

focal函数的两个重要参数（即邻域形状和邻域单元大小）对地形起伏度的影响至关重要，因此需要对邻域分析问题展开深入探讨。在ArcGIS软件中基于邻域分析函数包括空间分析模块下Neighborhood Statistic工具和GRID模块下的focal函数。由于focal函数可以通过AML语言编程反复调用，更适用于运算量较大的计算，因此本文采用focal函数计算地形起伏度。运用focal函数展开邻域计算时，先审查当前单元格邻域内的所有单元，此时邻域相当于一个滤波窗口。根据所选统计类型开展邻域计算，并将统计结果赋给输出单元，然后滤波窗口平移到下一个格子，重复上述过程，直到所有栅格计算完毕，生成新栅格图层，邻域计算完成。Arc GIS软件提供了矩形、圆形、环形、扇形及自定义的不规则多边形5种邻域，在地形起伏度计算中多采用矩形或圆形邻域［2］，本研究采用矩形邻域。

2.3 多尺度地形起伏度计算

邻域范围直接影响地形起伏的计算精度，一般情况下，地形起伏度开始会随着邻域统计范围的增大而显著增大，当增加到一定阈值后，其随邻域统计范围增大不会发生显著变化。因此基于邻域分析方法计算地形起伏度时，应将邻域统计范围逐步扩大，计算不同尺度下的地形起伏度。目前最被认可的方法是规则窗口递增法［14］。根据式（1），以n×n（n＝3，7，11，…，167）增加步长为4的栅格矩形作为模版算子，3×3为起始窗口，对整个研究区进行遍历计算，到167×167终止。由于运算量很大，因此采用AML语言编程调用focal函数来实现，程序设计流程如图1所示。

图1 地形起伏度计算程序设计流程Fig.1 Program design flow of calculating relief amplitude

3 结果分析

3.1 地形起伏度与邻域单元变化分析

邻域统计单元不同所计算出的地形起伏度差异显著，而且多数邻域单元的计算效果并不理想。因此定量分析地形起伏度与邻域统计范围关系是确定最佳邻域尺度单元的基础。将地形起伏度与邻域统计范围关系综合成表1，在EXCEL中对表1中数据绘制散点图并进行Logarithmic对数方程拟合（图2）。

表1 地形起伏度与邻域面积关系Table 1 The relationship between relief amplitude and neighborhood area

图2 地形起伏度与邻域尺度拟合曲线Fig.2 The fitting curve of relief amplitude and neighborhood scales

图2的拟合统计结果经拟合优度检验，拟合度较好。拟合曲线总体变化规律是：起伏度值先随统计格网单元面积增大而快速增大，当格网单元面积达一定阈值后，其增大速度开始减缓，最后趋于平稳，且在增速开始减慢过程中存有一明显拐点，即为曲线由陡变缓的阈值，该值所对应的邻域统计单元即为最佳尺度单元。

3.2 最佳邻域尺度单元

常规的确定最佳尺度单元方法是人工作图判定法和最大高差法。人工作图判定法虽然可以快速地判定曲线上的拐点，但其判定受人为主观因素的影响较大，往往精度不高;最大高差法容易产生多个拐点，多作为人工作图判定法的辅助方法。针对常规方法的不足，本文采用基于计算相邻尺度最大起伏度差与邻域面积差的比值的高差显著性变化检验法来确定最佳尺度单元。

基于正方形邻域获得的多尺度地形起伏度高差显著性变化检验公式为∶

式中：Ij表示地形起伏度随面积变化的强度，Rj表示j×j尺度地形起伏度的最大值，Sj表示j×j尺度地形起伏度计算时的邻域面积，（Sj＋4－Sj）/Sj＋4是将面积变化幅度标准化。Ij越大，表示单位邻域增长面积对应的地形起伏度变化越大;当Ij最大时，即出现对应的地形起伏度拐点，该点是地形起伏度所对应的最佳面积单元。

依据上述公式，得到高差显著性值随尺度变化图（图3）。可以看出，正方形邻域的最佳尺度单元位于第19个点，对应邻域单元大小为75×75个栅格单元，面积为5.0625 km2。采用同样的计算与检验方法，基于90 m空间分辨率的SRTM-DEM数据再次计算研究区地形起伏度，得到最佳邻域尺度单元为14.29 km2，由此可见，DEM的精度对地形起伏度的计算至关重要。

图3 变化强度Ij随邻域尺度变化Fig.3 The change of Ij with neighborhood scales

3.3 地形起伏度分级

参考中国1∶100万数字地貌制图规范中的分级标准［15］，对最佳尺度下的地形起伏度进行分级，生成地形起伏度分级图（图4）。研究区地形起伏度的总体趋势是：从西北向东南方向起伏度逐步增加。玛多、达日、色达3县地势较平坦，地形起伏度小于500 m，壤塘、马尔康、金川3县地势起伏剧烈，绝大部分地区地形起伏度在500 m以上。从各级地形起伏度所占比例结构看，200～500 m起伏度所占比例高达40.12%;其次是500～1 000 m，比例为26.05%;起伏度小于30 m的区域所占比例最小，仅有1.53%;起伏度大于1 000 m的区域有3 924.12 km2，占5.01%。由此可见，研究区地势以中度起伏为主。

图4 地形起伏度分级Fig.4 The class of relief amplitude

4 结论

基于ASTER GDEM数据和ArcGIS平台AML编程，调用GRID模块中用于邻域分析的focal函数，可以高效计算不同邻域尺度统计单元下地形起伏度，是一种有效的窗口递增分析方法。邻域统计单元大小对地形起伏度计算至关重要。两者之间的拟合关系表明：起伏度值先随邻域统计单元面积增大而快速增大，当邻域统计单元面积达一定阈值后，其增大速度开始减缓，最后趋于平稳，且在增速减慢过程中存在一明显拐点。经高差显著性变化检验法确定最佳邻域统计单元面积为5.06 km2。研究区自西北向东南地形起伏度逐步增加，地势以中度起伏为主。由于研究区面积较大，且ASTER GDEM数据空间分辨率高，导致地形起伏度计算过程运算量很大，所以窗口分析的步长设置为4，如果将步长缩小，可能会得到更精确的结果。

［1］汤国安，杨玮莹，杨昕，等.对DEM地形定量因子挖掘中若干问题的探讨［J］.测绘科学，2003，28（1）：28－32.

［2］汤国安，刘学军，闾国年.数字高程模型及地学分析的原理与方法［M］.北京：科学出版社，2005.

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

［4］刘新华，张晓萍，杨勤科，等.不同尺度下影响水土流失地形因子指标的分析与选取［J］.西北农林科技大学学报（自然科学版），2004，32（6）：107－110.

［5］唐飞，陈曦，程维明，等.基于DEM的准噶尔盆地及其西北山区地势起伏度研究［J］.干旱区地理，2006，29（3）：388－392.

［6］郎玲玲，程维明，朱启疆，等.多尺度DEM提取地势起伏度的对比分析——以福建低山丘陵区为例［J］.地球信息科学，2007，9（6）：1－6.

［7］龙恩.基于Srtm-DEM与遥感的长白山基本地貌类型提取方法［J］.山地学报，2007，25（5）：557－565.

［8］张磊.基于地形起伏度的地貌形态划分研究——以京津冀地区为例［D］.河北师范大学，2009.

［9］王玲，同小娟.基于变点分析的地形起伏度研究［J］.地理与地理信息科学，2007，23（6）：65－67.

［10］张学儒.基于地貌和植被的青藏高原东部样带土地类型研究［D］.中国科学院地理科学与资源研究所，2011.

［11］刘林山.黄河源地区高寒草甸退化研究——以达日县为例［D］.中国科学院地理科学与资源研究所，2006.

［12］阎建忠，张镱锂，摆万奇，等.大渡河上游生计方式的时空格局与土地利用/覆被变化［J］.农业工程学报，2005，21（3）：83－89.

［13］朱红春，陈楠，刘海英，等.自1∶10000比例尺DEM提取地形起伏度——以陕北黄土高原的实验为例［J］.测绘科学，2005，30（4）：86－88.

［14］周侗，龙毅，汤国安，等.面向DEM地形复杂度分析的分形方法研究［J］.地理与地理信息科学，2006，22（1）：26－30.

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

Analysis on the Relief Amplitude Based on ASTER GDEM Data in Mountain Area of Eastern Tibetan Plateau

ZHANG Xue－ru，GUAN Dong－jie，MOU Feng－yun，CHEN Chun
（CollegeofRiverSea，ChongqingJiaotongUniversity，Chongqing400074，China）

Relief amplitude is an indicator that can reflect change of rolling surface topography.It is widely used in terrain mapping，evaluation of resources and the environment and so on.This paper takes mountain area of the eastern Tibetan Plateau as study area.The relief amplitude in different neighborhood scale were calculated by AML programming invoking focal functions in ArcGIS base on ASTER GDEM data with high spatial resolution.The researches show that：Neighborhood scales have a critical effect on relief amplitude.In beginning，the values of relief amplitude increase rapidly with growing in area of neighborhood statistics units.But when the area of neighborhood units reaches a certain threshold value，the growth rate begins to slow down and stabilize.And there is a significant inflection point that is best neighborhood-scale unit when the growth rate begins to slow down.The best neighborhood-scale unit for measuring relief amplitude is 5.0625 km2by significance test.Class of relief amplitude is got，and it is found that the relief amplitude gradually increases from northwest to southeast in study area.The relief amplitude between 200 m to 500 m has the highest proportion，about 40.12%，and level of 500～1 000 m accounts for 26.05%.Therefore fluctuation of terrain is mainly moderate in study area.

relief amplitude;ASTER GDEM;neighborhood analysis;mountain area of eastern Tibetan Plateau

P 931;P208

A

1672－0504（2012）03－0011－04

2011－11- 30;

2012－02－13

国家自然科学基金项目（40901057）;国家科技支撑计划项目（2006BAB15B02-04）;河北省重点基础研究项目（08966712D）

张学儒（1982－），男，博士，讲师，主要从事土地利用与GIS应用研究。E－mail：zhangxueru5＠126.com