Python辅助下的坡度提取不确定性分析

吴 馁，李朝奎，方 军，褚 楠，吴柏燕

(1. 湖南科技大学地理空间信息技术国家地方联合工程实验室，湖南 湘潭 411201； 2. 湖南科技大学地理空间信息湖南省工程实验室，湖南 湘潭 411201)

滑坡和泥石流是主要的地质灾害形式，其基本形成条件离不开地形地貌[1-2]，而坡度是地貌学中描述地貌形态的两个重要因子之一[3]，对地表过程有强烈影响，具有重要的地理学意义[4]。

在区域尺度的研究中，坡度通常基于DEM提取[5]，DEM只是真实地表形态的近似表达[6]，而不同分辨率DEM数据提取的坡度精度存在差异，地形起伏等因素导致差异增大。因此，了解坡度提取误差的成因、大小和空间分布规律有利于改进坡度提取方法，提高坡度提取的精度。

国内许多学者在基于DEM提取坡度的不确定性分析中作了大量的研究。赵牡丹等选取陕北黄土丘陵沟壑区、黄土梁峁丘陵区、黄土低丘区及黄土破碎塬区4种地貌类型区作为试验样区，采用高精度1∶1万比例尺5 m分辨率的DEM为基准数据，应用比较分析、回归分析和相关分析等方法研究黄土高原不同分辨率DEM提取地面平均坡度精度的量化估算方法[7]；汤国安等选择陕北黄土高原6个典型地貌类型区为试验样区，采用野外实测及高精度的1∶1万比例尺DEM为基准数据，研究栅格分辨率及地形粗糙度对DEM所提取地面平均坡度精度的影响[8]；刘敏等以陕西省671幅1∶5万比例尺，25 m分辨率的DEM数据为信息源，分析DEM分辨率对提取坡度信息的不确定性影响[9]。

上述研究多从理论层面研究坡度提取的不确定性，但很少考虑整个计算过程的方法误差，以及可能存在的人为误差。本文针对现有坡度提取方法，运用以Python为开发语言的地理建模技术，从实践层面研究DEM分辨率对坡度提取的影响，整合分析过程，提高分析效率，减少人为误差。研究结论对完善现有坡度提取理论具有显著的理论意义。

1 坡度计算及不确定性分析方法

1.1 基于DEM的坡度计算

地表面上任一点的坡度是指过该点的切平面与水平地面的夹角[10]，表示地表面在该点的倾斜程度。因此，坡度是一个既有大小又有方向的量，即矢量。在数值上等于过该点的地表微分单元的法矢量n与地表曲面函数Z轴的夹角,如图1所示。用公式表示为

图1 地表单元坡度示意图

(1)

式中，n为地面一点的微分单元的法矢量；z为Z轴分量；|n|与|z|分别为各自方向的模。地面上每一点都有坡度，它是一个点上的微分概念，是地表曲面函数z=f(x,y)在东西、南北方向上的高程变化率的函数。因此，进行坡度提取常采用简化差分公式计算

(2)

式中，fx为X方向的高程变化率；fy为Y方向的高程变化率。由式(2)可知，fx与fy的求解是获取坡度值的关键。由于DEM数据是离散的，曲面函数f(x,y)未知，故坡度的提取通常是在3×3的栅格分析窗口中(如图2所示)，采用几何平面来拟合或差分计算的方法进行。

图2 三阶差分

基于DEM格网计算fx、fy可以有以下3种不同的计算方法

(3)

(4)

(5)

式中，d表示DEM的分辨率，即栅格格网的距离大小；fx与fy的含义与式(2)同。从式(3)—式(5)中可以看出，fx与fy的值主要受格网大小d的影响，因此坡度的提取精度也将明显受到格网分辨率的限制。为了提高坡度计算精度，式(5)常被用于求解fx与fy。

1.2 坡度计算的不确定性分析方法

Python作为一种通用的开源编程语言，易于学习和使用，在ArcGIS地理处理中用作脚本语言。而ArcGIS中的Arcpy站点包给Python处理地理数据提供了基础[11-12]。本文基于Arcpy站点包编写Python脚本程序对DEM数据进行批量重采样、批量坡度提取及批量掩膜提取。利用numpy包对数据进行整理，scipy包对数据进行科学计算并分析，matplotlib包将数据分析结果制成图表。

在Model Builder中构建模型工具。Model Builder也叫模型构建器，是ArcGIS提供的一种通过构造工作流与处理脚本的图形化建模工具。利用模型构建器将需要用到的工具连接，前一个工具的输出作为后一个工具的输入，形成工作流，从而不用频繁调用工具。本文使用Model Builder构建栅格河网模型计算沟壑密度。为了便于操作，将脚本工具和模型工具都集成到ArcGIS的工具箱中统一管理与调用。

2 试验及结果分析

2.1 试验数据及试验流程

采用平均坡度来表征区域坡度的一般水平[13]。以天津蓟县为试验区，地貌类型如图3所示。根据地貌类型将其划分为侵蚀中低山、侵蚀丘陵、侵蚀堆积山间盆地、平原洼地、洪积冲积倾斜平原、冲积平原、冲积洪积倾斜平原7个地貌类别区域。

以30 m空间分辨率的DEM数据为基础数据，对上述7类地貌类型区坡度提取的不确定性进行分析，具体试验流程如下：

(1) 在批量重采样Python脚本代码中，调用arcpy.Resample_management()函数，使用循环将30 m分辨率DEM重采样成40、50、60、70、80、90、100、110、120 m共10组不同分辨率的DEM数据。采样结果如图4所示。

注：1 miles=1.609 344 km图3 天津蓟县地貌类型

图4 不同分辨率DEM

(2) 在批量坡度提取Python脚本中，调用arcpy.Slope()函数，使用循环功能对上述重采样的DEM数据分别提取坡度值，如图5所示。调用掩膜提取函数arcpy.ExtractByMask()提取出每种地貌类型区域的平均坡度，结果见表1。

图5 不同分辨率DEM提取的坡度

(°)

(3) 运用分辨率与平均坡度分析模型进行分析计算，得到以DEM分辨率为横坐标、7种地貌类型区域的平均坡度值为纵坐标的散点图，并对散点图进行拟合，通过回归算法求解得到回归方程的系数及常数项。

(4) 建立不同地貌类型区域的4种地形因子，即地形曲率、地表起伏度、地面粗糙度及沟壑密度，以此作为坡度提取的地形因子并构建对应的模型工具。分别拟合上述回归方程中求得的常数项，建立平均坡度与DEM分辨率、区域地貌特征的关系。以沟壑密度模型提取为例，其提取模型及过程[14-15]如图6所示，其他因子提取与此相同。

图6 沟壑密度提取模型

2.2 试验结果分析

以DEM的分辨率为自变量，以不同地貌区域内的平均坡度值为因变量，定量分析平均坡度的变异特征。采用分辨率与平均坡度模型计算功能，其计算结果如图7所示。

图7 不同地貌单元的分辨率与平均坡度回归关系

由图7可知，平均坡度值与分辨率呈现很强的线性关系，且是正相关关系。因此以线性方程来描述这种关系，方程式如下

y=ax+b

(6)

式中，x为DEM分辨率；y为平均坡度；a为回归方程系数；b为常数项。则7种地貌类型区域内的平均坡度与分辨率的回归关系计算结果见表2。

表2 不同地貌单元内平均坡度与分辨率的回归关系

从表2可见，除了侵蚀中低山地貌类型，其他地貌类型的曲率基本一致，故取其他6种地貌类型曲率的平均值作为研究区的平均曲率，即a=-0.012 4，分别对4种地形因子模型进行计算，其计算结果如图8所示。

由图8可见，不同地貌单元内的地形因子沟壑密度值与相应的回归方程的常数项b存在显著的二次函数关系，其二次回归方程如下

b=0.004 7x2-0.095 5x+1.247 1

(7)

式中，x为不同地貌类型所对应的沟壑密度值；b为不同地貌类型区域内平均坡度与分辨率线性回归方程的常数项。由此得到试验区域内平均坡度的变异特征如下

F(g，x)=-0.012 4g+(0.004 7x2-0.095 5x+1.247 1)

(8)

式中，g为DEM空间分辨率；F(g，x)为对应的平均坡度值；x为沟壑密度值。

图8 不同地形因子与常数项b的关系

假定以分辨率g=30 m的DEM提取的平均坡度值为参考M，则M可以表示为

M(x)=-0.372+0.004 7x2-0.095 5x+1.247 1

(9)

式中，x的含义与式(7)、式(8)同。则以其他分辨率提取的平均坡度的误差E估算模型如下

E=M(x)-F(g,x)= -0.372+0.012 4g

(10)

式中，g为DEM空间分辨率；E为试验区域内平均坡度提取误差。可见坡度的提取平均误差只与DEM的空间分辨率有关。

3 结 论

针对现有基于DEM的坡度提取方法，通过地理建模技术，从实践层面研究DEM分辨率对坡度提取的影响，整合分析过程，提高分析效率，减少人为误差。以天津蓟县为研究区，采用ArcGIS平台和Python脚本语言功能分析了基于DEM的坡度提取方法的不确定性。研究结果表明：

(1) 研究区在不同地貌类型下提取的平均坡度与DEM空间分辨率之间呈现很强的线性关系，平均坡度随着DEM分辨率的减小而减小，但衰减速率基本不变。

(2) 不同地貌单元平均坡度与DEM分辨率之间回归方程的常数项与沟壑密度呈显著的二次函数特征。

(3) 坡度提取的精度与DEM的分辨率呈正相关，即坡度的提取精度随DEM分辨率的升高而升高。

(4) 基于Python的地理建模技术能有效地整合坡度提取精度分析内容，极大地提高坡度提取的不确定性分析效率。

研究结论进一步验证了当前基于DEM的坡度提取理论方法的实用性，是对现有理论方法的有效补充，具有一定的理论价值。

参考文献：

[1] 刘成，徐刚，杨娟，等.激发坡面泥石流发生的坡度因素探讨——以重庆北碚地区为例[J].中国地质灾害与防治学报，2006，17(4)：165-168.

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

[3] 杨勤科，贾大韦，李锐，等.基于DEM的坡度研究——现状与展望[J].水土保持通报 ,2007，27(1)：146-150.

[4] 李天文，刘学军，汤国安.地形复杂度对坡度坡向的影响[J].山地学报，2004，22(3)：272-277.

[5] 马思煜，王春梅，土祥，等.坡度随DEM分辨率变化的地统计学解释——以黄土丘陵沟壑区为例[J].地下水，2015，37(5)：278-280.

[6] 陈楠.基于信息损失量的坡度精度分析[J].地球信息科学，2014，16(6)：852-858.

[7] 赵牡丹，陈正江，晋锐.DEM提取地面平均坡度误差的量化模拟[J].水土保持通报，2002，22(3)：55-57.

[8] 汤国安，赵牡丹，李天文，等.DEM提取黄土高原地面坡度的不确定性[J].地理学报，2001，58(6)：824-830.

[9] 刘敏，汤国安，王春，等.DEM提取坡度信息的不确定性分析[J].地球信息科学，2007，19(2)：65-69.

[10] 汤国安，李发源，刘学军.数字高程模型教程[M].3版.北京:科学出版社，2016.

[11] 刘雄，李宇，陈刚.Arcpy支持下的线要素几何自动接边与属性检查[J].测绘通报,2016(8)：114-117.

[12] 余咏胜,彭艳丽,尹言军,等.基于Arcpy的影像地图自动处理技术研究[J].测绘通报,2015(3)：82-85.

[13] 胡卓玮，李洋，王志恒.基于DEM的四川省低山丘陵区坡度提取不确定性分析[J].山地学报，2012，30(5)：636-640.

[14] 师长兴，周园园，范小黎，等.利用DEM进行黄河中游河网提取及河网密度空间差异分析[J].测绘通报,2012(10)：24-27.

[15] 刘学军，卢华兴,卞璐，等.基于DEM的河网提取算法的比较[J].水利学报，2006，37(9):1134-1141.

欢迎订阅《测绘通报》

《测绘通报》创刊于1955年，是由国家测绘地理信息局主管、中国地图出版社(测绘出版社)主办的反映我国测绘地理信息科技发展现状和指导全国测绘地理信息生产业务的综合性、技术性刊物。现为《中文核心期刊要目总览》收录期刊、中国科技核心期刊(中国科技论文统计源期刊)、中国科学引文数据库(CSCD)来源期刊、中国学术期刊网络出版总库(CAJD)全文收录期刊、中国学术期刊综合评价数据库(CAJCED)核心统计源期刊、《中国学术期刊文摘》收录期刊。曾被评为第二、三届国家期刊奖百种重点期刊，多次荣获全国优秀测绘地理信息期刊奖。

《测绘通报》的主要内容包括大地测量、全球导航卫星系统(GNSS)、摄影测量、遥感(RS)、地图制图、地理信息系统(GIS)、工程测量、矿山测量、地籍测绘、海洋测绘、测绘仪器、信息传输、图形图像处理等方面的新成果和新技术；行业管理、科研、教学、生产的先进经验；计算机、通信等相关理论技术在测绘地理信息领域里的应用及测绘地理信息科技在国家经济建设各个方面的应用；国内外测绘地理信息学术动态及有关测绘地理信息科技信息等。

《测绘通报》的主要栏目有学术研究、技术交流、测绘地理信息论坛、经验介绍、测量员之窗、国外测绘地理信息、行业管理、行业研究、行业调查、知识窗、新书介绍、企业之窗、测绘地理信息教学等。

《测绘通报》为月刊，2018年定价：35.00元，邮发代号：2-223。

编辑部地址：北京市西城区三里河路50号，邮编：100045，订阅电话：010-68531192(金老师)，68531317(传真)。

网址：http:∥tb.sinomaps.com