规则网格DEM数据提取山脊线的应用

栗 业

(河北省地质矿产勘查开发局第二地质大队/河北省矿山环境修复治理技术中心，河北 唐山 063000)

1 引言

规则网格DEM实质是一种栅格形式的数据，网格大小即为像素大小(分辨率)，高程值为灰度值。数字高程模型是通过有限的高程值对地形曲面的数字化模拟或者是地表形态的数字化表达[1]，英文是Digital Elevation Model，简称DEM。从狭义角度说DEM是区域地表海拔的数字化表达；广义上是地理空间中的地理对象表面海拔高度的数字化表达。

数学意义上的DEM是表示区域A上的三维向量有限序列，函数表达为：

Ai=(Xi,Yi,Zi) (i=1,2,3…，n)

(1)

式(1)中，(Xi,Yi)是平面坐标，Zi是(Xi,Yi)对应的高程。

数字高程模型(DEM)包含了大量的地性线信息。山脊线是构成地势起伏变化的重要特征线，对地形地貌研究有重要意义。基于DEM提取山脊线的研究在区域水文分析、土木工程建设、地质勘探选定资源靶区、遥感影像自动配准等诸多领域具有重要的实践应用价值和科学意义[2]。

目前国内外关于山脊山提取的研究众多并取得了良好的结果。研究热点主要集中在提取算法和DEM不同尺度两方面。从算法设计分类来说，大致有基于图像处理的算法、基于地表几何形态的算法、基于地形表面流水物理模型分析算法、基于地形表面几何形态分析和地表流水物理模型分析结合的算法、平面曲率与坡形组合法等5种方法[3,4,13]。其中，平面曲率与坡形组合法、基于地形表面流水物理模型分析算法是目前最常用的2种地性线提取方法。DEM是地性线提取的基础数据，DEM不同尺度单元网格的大小不同，代表了不同的地形模拟精度。但对山脊线的提取，并不是对地表模拟越逼真提取的越精确。因为地表模拟越逼真就越容易出现过度提取数据冗余的现象，因此选择适当尺度的DEM是当前研究的另一热点问题。本文分别基于网格大小为90 m、30 m和内插生成网格大小为20 m、10 m的规则格网DEM数据利用平面曲率和坡向变率(SOA)方法及水文分析方法提取山脊线，并对各方法各分辨率提取下的结果进行研究和分析。

2 提取应用实验

2.1 实验区概述

本文以西南横断山区某地为典型实验区域，此区域雪山林立、江河纵横、草原遍布、湖泊众多。地质处于“三江”褶皱系与扬子准台地交接地带，全县10%的面积为扬子准台地，90%为“三江”褶皱系。地质结构复杂，构造运动活跃，地势自西北向东南倾斜,地貌形态分为山地、高原、盆地、河谷地貌。对该地区山脊线信息的提取研究为今后西南地区的地质地貌分析、水文和水利工程建设有重要意义。

2.2 数据及软件

数字高程模型(DEM)。SRTM 90 m数据，空间分辨率为900 m，SRTM(Shuttle Radar Topography Mission)，由美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量。GDEM 30 m数据，本数据利用ASTER GDEM第一版本(V1)的数据进行加工得来，数字高程数据产品的空间分辨率为30 m。基于网格大小为90 m、30 m内插生成网格大小为20 m、10 m的规则格网DEM数据。

对比分析数据：分辨率为5 m的法国spot5卫星影像数据。

软件平台为地理信息系统软件ArcGIS10.3的空间分析工具。

2.3 提取方法

2.3.1 基于 DEM 平面曲率与坡形组合提取山脊线

2.3.1.1 基于平面曲率提取山脊线的方法

求取DEM的平面曲率和地表的正负地形，正地形上曲率的大值是山脊，负地形上曲率的大值是山谷。这种方法提取的宽度可由平面曲率的大小来调节，方法简便效果好，但提取的图形是栅格形式，在转为矢量线性方面还需进一步研究改进[5,6]。

利用DEM计算平面曲率并求出正地形，平面曲率图正地形图叠加，计算取得正地形上平面曲率上的大值即为山脊线。提取流程见图1所示。通过改变平面曲率的大小可以调节提取山脊线的宽度，此方法简便，效果好[12,14,15]。阈值需要根据实验区的晕渲图辅助确定，但主观干扰因素较多，需要丰富经验。规则格网DEM的曲率可以利用曲率分析模块直接计算得来，包括总曲率、剖面曲率及平面曲率。求交工具可以使用栅格计算器(Raster Calculator)。

图1 基于平面曲率提取流程

2.3.1.2 基于无误差坡向变率提取地形特征线的方法

坡向是一个环形变量数据，经过内插和代数运算后再应用原始的坡向变率会存在明显的误差。例如，坡向为330°的坡面和坡向为30°的坡面，两坡面的坡向差显然是60°，而不是330°-30°=300°。所以在实际提取中，利用无误差的坡向变率(SOA)代替平面曲率来提取山脊线的效果会更好。如图2所示，基于无误差的坡向变率提取流程有别与平面曲率的提取。利用该方法提取山脊线时，如何获取无误差坡向变率是问题的关键。无误差坡向变率SOA 可以利用栅格计算器通过式(2)计算得到。合理的阈值仍需通过晕渲图和等值线图不断比对判定得到。

SOA=(([SOA1]+[SOA2])-Abs([SOA1]

(2)

式(2)中，SOA1表示提取的坡度数据，SOA2表示利用SOA方法求出的反地形的坡向变率。

2.3.2 基于DEM 水文分析提取山脊线

根据水往低处流，最终顺谷而流，汇集成水系网的自然规律，Chakrevavanieh率先提出模拟地表流水物理模型。D8法是该模型中的最常用算法：假设只有1个单一的网格的水流进8个相邻的网格单元。它用最陡坡度法来确定水流的方向，即在3×3 的DEM栅格上，计算中心栅格与各相邻栅格间的距离

图2 基于坡向变率提取流程

权落差取距离权落差最大的栅格为中心栅格的流出栅格。以数值表示每个单元的流向。数字变化范围是1～255。其中1:东；2:东南；4:南；8:西南；16:西；32:西北；64:北；128:东北。然后在流向栅格图的基础上生成汇流栅格图。汇流栅格上每个单元的值代表上游汇流区内流入该单元的栅格点的总数，既汇入该单元的流入路径数较大者，可视为河谷，值为0，则是较高的地方，可能为流域分水岭。有了流域汇流栅格图就可以很容易地提取流域的各种特征参数。如模拟流域，选取一个阈值，大于该值的格点为沟谷线上的点，把各个点相连接就形成了沟谷线，即河道水系。但在应用时需要解决以下3个问题：①水道起始点位置的确定；②伪负地形的识别及填充；③凹陷与平坦处水流方向的设定。而且基于地表水流模拟方法存在2个缺陷：①洼地的填充和平坦地区水流方向的设定十分麻烦；②该方法所计算的汇水量与高程有关,计算的结果必然是高程值大的汇水量小，高程值小的汇水量大，由于自然特性的水从高处连续性的向下流动，使得在一个相对高位置山谷线上的某一点因其汇水量较小而被遗漏，而处于低位置非山谷处的某点因其汇水量较大而被误判为山谷线上的点。在后来的研究中，许多学者提出利用平滑DEM，垫高填平凹陷区，标识洼地集水区的方法解决上述问题并取得了良好效果。改进这种算法的基本思路是首先按照地表流水模型在稀疏的DEM中概略提取地性线；然后对周围地域进行几何地理形态分析，在精确定位地性线。该方法的关键是求出概略地性线与DEM网格的交点，在该点周围的小范围内进行几何分析，找出正交方向地形断面上的高程极值点，正是这一点为精确点。这种综合算法在特征点提取阈值的选取、DEM网格大小的确定方面还有待进一步研究[7～13,16,17]。

山脊线提取实质上是分水线提取，因此可以利用水文分析的方法提取。水文分析的主要内容是利用水文分析工具提取地表水、径流模型的水流方向、汇流累积量、水流长度、河网以及对研究区的流域进行分割等，还可以结合几何分析提取研究区域的山脊[16]。通过模拟地表流水的实际情况，以山脊线和山谷线的物理特性为依据，地表水文分析提取方法的基本原理，利用ArcGIS软件Spatial analyst tools的水文分析工具，操作步骤见图3，由于山脊线没有一个精确的提取对比标准，所以阈值可根据研究区域晕渲图、等值线图或卫星影像数据辅助判断进行选取。通过实验表明阈值直接影响汇流累积量，所以选取合理的阈值至关重要。

图3 水文分析提取方法流程

3 实验提取结果

经过不同方法，不同网格大小DEM的提取实验，得出以下结果。

基于无误差坡向变率提取山脊线结果(图4)：这4幅图的DEM的网格大小分别为10×10、20×20、30×30、90×90。通过晕渲图对阈值的赋值判定，经过反复试验，分别取值40、50、30、20。经对比分析，90×90的DEM提取效果最优，提取的山脊线连贯性好，几乎没有断裂，破碎度最低，噪音最少。但是提取的山脊线对应实际的地形位置有所偏差，偏差范围在2、3个网格之间，即实际180～270 m的偏离，侧山脊线提取不明显。30×30的DEM提取位置相对精确，主山脊线提取效果较好，侧山脊线出现了明显的断裂，噪声点增多。10×10、20×20的DEM提取的主山脊线符合实际地貌形态，精度高。但侧山脊线全是破碎的断裂，只有主山脊线提取较完好。10×10的DEM数据中主山脊线也出现了断裂。由于10×10、20×20的DEM数据是重采样数据，在提取坡度时产生了明显的误差和数据丢失，结果产生了呈格网状的噪声。

图4 基于无误差坡向变率提取山脊线结果

基于DEM 水文分析提取山脊线结果(图5)：这4幅图的DEM的网格大小分别为10×10、20×20、30×30、90×90。同样通过晕渲图对阈值的赋值判定，分别

图5 基于DEM 水文分析提取山脊线结果

取值0.5037、0.3379、0.6612、0.5558。经对比分析，90×90的DEM提取出的山脊线连贯性好，断裂少。但有些地区的山脊线未被提取出来，还有些地区提取过度。提取结果也与实际山脊线位置有偏差，在地貌复杂和相对海波较高的区域偏差小，在相对平坦地区偏差较大。20×20、30×30的DEM数据提取效果较好，虽然出现了断裂和破碎噪声，但山脊线基本被提取出来而且，位置相对准确。10×10的DEM数据提取位置基本和实际地形相符，但是断裂太多，破碎噪声太多，后期处理繁琐。

对于破碎噪声，本次实验利用众数滤波进行处理，即根据相邻像元数据值的众数替换栅格中的像元。众数滤波工具需要满足2个条件才能发生替换：一是，具有近似值的相邻像元数必须足够多，并且这些像元在滤波器内核周围必须是连续的；二是，与像元的空间连通性有关，目的是将像元的空间模式的破坏程度降到最低。在具体操作中，相邻像元数使用FOUR会保留矩形区域的拐角。使用EIGHT将使矩形区域的拐角变得平滑。如果将替换阈值指定为HALF，并且2个值的出现次数相等，则当处理的像元值与其中某一半的值相同时将不会发生替换。HALF比MAJORITY的过滤范围广泛。当边和角栅格像元的相邻条件相同时，它们会遵循不同的MAJORITY和HALF规则。使用FOUR内核时，边或角像元始终要求存在2个匹配的相邻像元才能发生替换。使用EIGHT内核时，角像元在所有相邻像元均具有相同值时才能发生更改，而边像元需要3个相邻像元(包括边上的像元)具有相同值才发生更改。运行几次众数滤波后，输出栅格将会稳定下来不再发生变化。图6为滤波处理后的效果图。

图6 众数滤波处理效果

4 结论与展望

提取的山脊线重复使用众数滤波工具可有效减少提取山脊线的破碎噪声。通过2种提取方法和不同网格大小的比较发现：并不是网格越小分辨率越高，地表模拟越精确提取的效果越好。阈值选择仍靠遥感影像或晕渲图辅助判断，没有定量算法，个人主观因素多，需要丰富的经验。不同方法应根据各自的特点选取不同网格大小的DEM数据和阈值。内插DEM产生的网格状噪音是由于网格间坡度信息不断丢失，地面表达能力逐步下降造成的。网格细化内插的过程可以近似看作对DEM数据不断平滑填洼的一个过程，在这个过程中虽然DEM平滑性得到了保障，但坡度信息也在平滑过程中不断丢失。由此看出，一个合理的DEM平滑度与阈值结合对地性线的提取至关重要。

在未来的研究中，提取过程阈值的定量确定的问题，提取成果线性断裂和破碎噪音较多的问题，批量自动化提取的问题需进一步完善和改进。山脊线、山谷线等地性线的提取，在分析区域水文分析、探矿找矿、水利工程建设等方面有着广泛的应用前景，今后在研究和应用方面还要努力做大量的研究工作。