基于ASTER GDEM的县域尺度地形起伏度分析——以华蓥市为例

陈珂，黄小羽，罗明良

(西华师范大学 国土资源学院，四川 南充 637009)

1 引 言

起伏度是一定面积内最高和最低点之高差[1]，其研究源于前苏联学者提出的切割深度，现已成为划分地貌类型的重要指标，也是地貌图编制的基本划分依据[2]。利用数字高程模型作为基本信息源对区域地势起伏度的研究逐渐兴盛起来[3]，业已完成不同比例尺起伏度最佳统计单元求解[4-7]。起伏度研究在诸多方面应用广泛，如区域人口分布和经济发展[8-9]、区域土地时空变化和土地利用/土地覆盖空间结构[10]，以及用于水土流失定量评价及区域滑坡灾害评价等[5，11]。

基于DEM提取地形起伏度的关键在于确定一个合适的分析窗口，使得在宏观范围内求取的起伏度能够准确反映地面的起伏状况[1]。用地形起伏度分析使用窗口递增方法，得到平均起伏度趋于平稳的曲线，通过定性判别或均值变点等方法得到最佳分析单元[3]。现有研究已经表明窗口递增法适用于大区域、较低分辨率的DEM[12]；面临县域尺度、较高分辨率的DEM，平均起伏度趋向不断增加，如何确定窗口递增法的最大分析窗口？ASTER GDEM作为全球公开共享数据，具有约30m水平分辨率[13]，在地貌研究中具有越来越广泛的应用空间[13-14]，探索基于30m×30m栅格分辨率DEM提取县域尺度地形起伏度的方法，并结合实际地形对其结果进行地学解释和分析显得尤为必要。

本文以ASTER GDEM数据，以川东平行岭谷区的华蓥市为研究区域，使用窗口递增方法，考察最大起伏度随分析窗口变化趋势，确定最大分析窗口，得到最大起伏度随分析窗口趋于平稳的曲线；基于均值变点法选取适宜分析窗口，得到华蓥市地形起伏度。研究对于深入认识华蓥市地貌形态有借鉴意义。

2 研究样区与数据

2.1 研究样区概况

研究样区为华蓥市，部分位于国家级地质公园内，地处四川盆地川东平行岭谷区的华蓥山中段西缘、渠江东岸，海拔介于206m～1704m；高登山为最高点。襄渝铁路沿华蓥山脚纵贯南北，将全市自然分为大致相等的东西两大部分，东西两部分的地形地貌有着明显差异。西部低丘广布，地势偏低，以丘陵和低丘为主，深丘较少；平坝、台地散布其间；多紫色土，土壤肥沃，灌溉便利，是主要的农耕区。东部山地为主，地势较高。东南部是华蓥山中段的西翼部分，低山为主，海拔超过1000m的华蓥山山脊一带出现中山地貌。东部山区石灰岩广布，喀斯特地貌突出，有成片的石林和众多溶洞，如瓦店的石林、溪口的仙鹤洞等。土壤为黄壤，土层薄，肥力低，农耕条件差，自然资源丰富。

2.2 实验数据

本文所采用的数据是ASTER GDEM(先进星载热发射和反射辐射仪全球数字高程模型)数据，该数据是由美国国家航空航天局(NASA)与日本经济产业省(METI)于2009年6月共同推出的，是目前较精确的地球陆地海拔地形数据。空间分辨率为1×1弧秒(约30m×30m)，每个分片包含3601行×3601列。全球范围内，在置信度为95%时，ASTER GDEM的垂直精度约为20m，水平精度达30 m[13]。

3 研究方法与步骤

3.1 起伏度求解

基于ArcMap空间分析(Spatial Analyst)提取地形起伏度，使用邻域分析(Neighborhood Statistics)实现。分析中使用n×n(n=3，5，7，…，83)分析窗口，计算窗口内最大高程与最小高程之差；并求取正方形窗口下的平均地形起伏度及最大起伏度，并依此数据对其进行拟合，得到网格单元面积与平均地形起伏度的关系。

3.2 最佳统计单元分析

利用窗口递增法可以得到地形起伏度随分析窗口变化曲线，曲线形态变化有助于推断起伏度最佳统计单元。为了避免目视解译带来的主观性，沿用了均值变点分析法，以合理计算出曲线上由陡变缓的点，即最佳统计单元。均值变点分析法是一种数理统计学上的对非线性数据进行处理的方法，其计算步骤如下[3]：

(1)令i=2，…，N，对每个i将样本分为两段：X1，X2，Xi-1和Xi，Xi+1，…，XN。计算每段样本的算术平均值Xi1和Xi2及统计量Si：

(1)

(2)计算统计量：

(2)

(3)计算期望值：E(S-Si)，i=2，3，…，N；其中S为原始样本的统计量，Si为样本分段后的统计量。

(3)

基于上述方法，作出不同窗口下S与Si差值的变化曲线，找到曲线上S与Si差值的最大值点，通过对比分析，得到适合本区域的最佳统计单元面积。

4 实验结果

4.1 不同窗口与平均、最大地形起伏度的对应关系

最大分析窗口的确定是窗口递增法的重要步骤之一。华蓥市面积仅470km2，区域较小；研究使用的DEM分辨率较高；属于典型的小区域、较高分辨率条件下的地形起伏度研究。研究表1及图1(a)可以发现，在分析窗口递增条件下，地形起伏度均值由陡变缓趋势不明显。

按照地貌发育的基本理论，该区域存在使最大高差达到相对稳定的最佳分析面积[4]。基于此，研究摈弃了根据平均起伏度判定终止分析窗口的一贯做法，转而考察区域最大起伏度何时趋于稳定。从表1和图1(b)可以看出，分析窗口扩大至79～83时，最大地形起伏度保持1182m不变。据此，分析使用的最大分析窗口为83×83，面积约为6.2km2。图1(a)为分析窗口与平均地形起伏度对应关系曲线；图1(b)为分析窗口与最大地形起伏度对应关系曲线。

图1 分析窗口与平均、最大起伏度对应关系图

表1分析窗口与平均起伏度(m)、最大起伏度(m)及窗口面积(104m2)对应关系

窗口大小平均起伏度最大起伏度窗口面积窗口大小平均起伏度最大起伏度窗口面积3×316.08200.000.8145×45199.66936.00182.255×531.10289.002.2547×47205.41958.00198.817×744.88353.004.4149×49211.05974.00216.099×957.52396.007.2951×51216.58988.00234.0911×1169.18431.0010.8953×53222.01995.00252.8113×1380.03466.0015.2155×55227.351005.00272.2515×1590.18518.0020.2557×57232.631046.00292.4117×1799.75552.0026.0159×59237.811061.00313.2919×19108.81589.0032.4961×61242.911075.00334.8921×21117.43628.0039.6963×63247.941087.00357.2123×23125.67703.0047.6165×65252.881096.00380.2525×25133.57756.0056.2567×67257.741110.00404.0127×27141.16768.0065.6169×69262.521137.00428.4929×29148.47789.0075.6971×71267.221151.00453.6931×31155.53814.0086.4973×73271.851163.00479.6133×33162.35833.0098.0175×75276.411175.00506.2535×35168.97847.00110.2577×77280.901181.00533.6137×37175.40860.00123.2179×79285.331182.00561.6939×39181.68872.00136.8981×81289.711182.00590.4941×41187.80894.00151.2983×83294.021182.00620.0143×43193.79916.00166.41

4.2 变点法寻求最佳统计单元

从图1(b)可见，拟合曲线存在由陡变缓的拐点，该点所对应的区域面积即所要求取的最佳分析区域面积。为进一步精确判断拐点的大小，研究采用均值变点法对其进行分析(图2)。从图2可以看出，系列分析窗口3×3、5×5，…，83×83对应的S与Si差值最大点为第13个点，其窗口大小为27×27，正方形窗口边长约为810m，对应的区域面积约为0.65km2。这一结果略小于朱红春等基于1∶1万DEM得到的黄土高原绥德样区的1km2[12]，大于程维明等基于1∶10万DEM得到的0.45km2[7]。

4.3 华蓥山地形起伏度分析

根据数字地貌图制图规范，地形起伏度可划分为7级：<30m、30m～70m、70m～200m、200m～500m、500m～1000m、1000m～2500m、>2500m。本文按此标准对华蓥山地形起伏度进行分级，用上述分析结果，以27×27分析窗口作为最佳统计单元，得到华蓥山地形起伏度分级图(图3)。

从华蓥山整体范围来看东部为山区，西部为阡-陵平坝。中低山地集中在东部地区，中东部地区存在岩溶中山、低山及岩溶槽谷组成；西部多丘陵、平坝等。其中东北部低起伏区为华蓥山天池湖所在地；中部及南部起伏度超过500m的地方为华蓥山山脊部分。华蓥山处于川东帚状构造带，由一系列近于平行的狭长不对称的箱状高背斜组成，轴部由于断裂破碎带存在，强烈差异侵蚀形成陡峻山体和深窄谷底，在图3中体现出较大的起伏度，最大起伏度768m，其中中低山部分约为全区面积的40.85%；这与华蓥市全市470km2，中低山面积198km2较为相符[15]。从图3中还可以看出，华蓥市起伏度在70m以下的部分约为总面积的17%，这与实际情况有一定差异。考察这部分区域，大多为植被覆盖的山麓地带或城市建成区。由于ASTER GDEM利用遥感影像立体像对建立DEM，在一定程度上不能完全消除建筑物、水系、云层及植被覆盖等的高度影响[16-17]，导致部分地区高程垂直精度异常。

图2 S与Si差值曲线

图3 华蓥山地形起伏度分级图

5 结束语

地形起伏度的量测是重要的地形因子之一，利用ASTER GDEM计算研究区的起伏度，并将结果与常规测量结果进行比较也是衡量DEM可用性和可靠性的重要方法。本文以川东平行岭谷华蓥市为例，基于ARSTER GDEM数据，通过窗口递增分析法、结合均值变点，对县域尺度华蓥市地形起伏度进行分析，结果表明，县域尺度、较高分辨率条件下，平均地形起伏度曲线趋向增大；但最大起伏度可寻找到趋于稳定的区间；均值变点探测的华蓥市的最佳栅格分析窗口边长约为27×27，对应最佳统计单元为0.65km2，以这一分析单元得到的华蓥市中低山面积约为全区的40.85%；其结果与建立在常规测量基础上的地理统计结果基本吻合，证明了分辨率为30m×30m的国际共享数据ASTER GDEM可以完成类似川东华蓥山脉及谷地的地貌起伏度计算。但研究区内起伏度小于70m的区域面积偏大，除受ASTER GDEM自身精度的影响，分析使用的窗口分析法也存在较大的边界效应，在边界区域得到的地形起伏度存在一定误差。后续研究应当侧重从地貌基础认知出发，考虑基于DEM的山地地貌实体表达方法[12，18]，以求更加客观地反映区域地形起伏度。

参考文献：

[1] 汤国安，李发源，刘学军.数字高程模型教程[M].北京：科学出版社，2010.

[2] 周成虎，程维明，钱金凯.数字地貌遥感解析与制图[M].北京：科学出版社，2009.

[3] 韩海辉，高婷，易欢，等.基于变点分析法提取地势起伏度——以青藏高原为例[J].地理科学，2012，32(1)：101-104.

[4] 涂汉明，刘振东.中国地势起伏度研究[J].测绘学报，1991，20(4)：311-319.

[5] 刘新华，杨勤科，汤国安.中国地形起伏度的提取及在水土流失定量评价中的应用[J].水土保持通报，2001，21(1)：57-59+62.

[6] 刘爱利，汤国安.中国地貌基本形态DEM的自动划分研究[J].地球信息科学，2006，8(4)：8-14.

[7] 郎玲玲，程维明，朱启疆，等.多尺度DEM提取地势起伏度的对比分析——以福建低山丘陵区为例[J].地球信息科学，2007，9(6)：1-6+135-136.

[8] 封志明，唐焰，杨艳昭，等.中国地形起伏度及其与人口分布的相关性[J].地理学报，2007，62(10)：1073-1082.

[9] 周自翔，李晶，任志远.基于GIS的关中-天水经济区地形起伏度与人口分布研究[J].地理科学，2012，32(8)：951-957.

[10] 谷建立，张海涛，陈家赢，等.基于DEM的县域土地利用空间自相关格局分析[J].农业工程学报，2012，28(23)：216-224.

[11] 郭芳芳，杨农，孟晖，等.地形起伏度和坡度分析在区域滑坡灾害评价中的应用[J].中国地质，2008，35(1)：131-143.

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

[13] 康晓伟，冯钟葵.ASTER GDEM数据介绍与程序读取[J].遥感信息，2011，36(6)：69-72.

[14] 张学儒，官冬杰，牟凤云，等.基于ASTER GDEM数据的青藏高原东部山区地形起伏度分析[J].地理与地理信息科学，2012，28(3)：11-14.

[15] 周申立.华蓥山滑坡研究[M].成都：西南交通大学出版社，2008.

[16] 郭笑怡，张洪岩，张正祥，等.ASTER-GDEM与SRTM3数据质量精度对比分析[J].遥感技术与应用，2011，26(3)：334-339.

[17] 朱圣军.ASTER GDEM数据精度分析及其在石油勘探中的应用[J].物探装备，2011，21(1)：54-59.

[18] 罗明良，汤国安.地貌认知及空间剖分的山顶点提取[J].测绘科学，2010，35(5)：126-127+253.