基于DEM的水系提取及集水阈值确定方法研究

罗大游, 温兴平, 沈 攀, 张皓楠, 张丽娟, 李进波, 郁 智

(1.昆明理工大学 国土资源工程学院, 云南 昆明 650093;2.贺州学院, 广西 贺州 542899; 3.云南省矿产资源预测评价工程实验室, 云南 昆明 650093)

基于DEM的水系提取及集水阈值确定方法研究

罗大游1,3, 温兴平1,3, 沈 攀1,3, 张皓楠1,3, 张丽娟2, 李进波1,3, 郁 智1,3

(1.昆明理工大学国土资源工程学院,云南昆明650093;2.贺州学院,广西贺州542899; 3.云南省矿产资源预测评价工程实验室,云南昆明650093)

[目的] 确定适宜的集水阈值，使自动提取的水系河网与实际河道相符。[方法] 以DEM数据为基础，利用ArcGIS水文分析模块对研究区流域河网水系进行自动提取，基于河网密度与集水阈值的相关性，在已提出的拟合函数一阶导数求转折点、二阶导数求拐点法确定集水阈值的基础上，以龙川江流域为例，提出求解河网密度变化率等于集水阈值变化率的数值方法，得到适宜的集水阈值。[结果] 通过设置不同集水阈值生成河网，发现不同集水阈值对主河道长度及地理空间位置影响较小，但对提取的河网特征影响较大，最终确定龙川江流域集水阈值设置为0.12 km2为宜。[结论] 集水阈值的确定影响着河网提取的精度，通过变化率确定集水阈值的方法主观因素较少，避免了人为干扰，相对客观，可对区域水土流失监测提供一定依据。

DEM; ArcGIS; 河网密度; 集水阈值

文献参数： 罗大游, 温兴平, 沈攀, 等.基于DEM的水系提取及集水阈值确定方法研究[J].水土保持通报，2017,37(4)：189-193.DOI:10.13961/j.cnki.stbctb.2017.04.032； Luo Dayou, Wen Xingping, Shen Pan, et al. Information extraction of river networks and determination of drainage area threshold using DEM data[J]. Bulletin of Soil and Water Conservation, 2017,37(4)：189-193.DOI:10.13961/j.cnki.stbctb.2017.04.032

DEM(digital elevation model)是对地形表面形态数字化表达的一种数字高程模型，它在水文领域中有着广泛的应用，是地形划分、水文分析的重要基础数据。通过利用ArcGIS水文分析模块，结合DEM数据自动获取包含河网在内的水文参数，对数据加以分析，最终实现地形模型可视化。目前水文特征提取与分析在水土保持中起到了重要作用，在暴雨洪水灾害中，通过DEM生成的水文模型提取的河网作为研究区危险性评价指标[1-2]；在泥石流、滑坡等地质灾害频发，水土流失严重地区，为灾后研究区的水土治理提供基础依据[3]；在河流、湖水流域作为生态环境治理实践数据支撑[4-5]。

以DEM数据为基础，利用GIS空间分析方法提取水系以提取研究区水文特征。国内外在这方面已开展了大量的工作[6-8]。国内在集水阈值的确定方面也取得了长足的进步。传统方法中孔凡哲等[9-10]观察集水阈值与河网密度的相关性，取河网密度变化趋于平缓时对应的值作为集水阈值，此方法简单，易于操作但人为主观性较高；数学方法中叶章蕊等[11]通过拟合函数一阶导数的求解得到转折点。关颖慧等[12]利用拟合函数二阶倒数求拐点以获取集水阈值。该方法将阈值确定在较小范围但可能出现拟合函数如幂函数无法求得转折点、拐点的情况；其他方法中王林等[13]基于网格法、Horton定理的估算方法计算流域水系分维数，结合水系分维与集水阈值的相关关系，把分维值趋于平缓的点作为集水阈值，此方法同样存在人为主观影响。在前人的基础上，提出基于河网密度与集水阈值的相关性，求解河网密度变化率等于集水阈值变化率的数值方法，将河网密度与集水阈值的幂函数同切线方程联立，求解唯一值作为河网提取的阈值。此方法避免了水系特征提取的人为主观性，具有速度快，成本低的优点，真实反映水系的空间分布规律的优点。

本文拟以龙川江流域为研究对象，DEM数据为基础，利用ArcGIS水文分析模块对流域进行水文特征提取分析，对了解地区地形发育特征、水土流失监测、地质灾害提取等有着重要意义。

1 研究区概况

龙川江为金沙江南岸一级支流，水源丰富，发源于云南省楚雄彝族自治州西部，由西向东流经楚雄市，横穿元谋坝区，最终在元谋北部的江边乡汇入金沙江。流域面积约6 109 km2，境内最高高程与入江口相对落差约2 000 m，流域平均海拔约1 850 m。该流域以山区和丘陵为主，属低纬度高原季风气候，受大气环流和季风影响，干湿季分明，降水量年内分配不均匀，雨季6—10月降水量占年降水量的80%以上，多年平均降水量847.2 mm，年最大降水量为1 075.2 mm，最小降水量为634.6 mm。2016年9月17日，因持续强降雨，研究区发生近年来最大山洪泥石流灾害，在龙川江内形成朱布、海洛2个堰塞湖，严重威胁着上下游区域人民群众的生命财产安全，因此龙川江是云南省防汛抗灾的重点区域。

2 数据来源及研究方法

2.1 数据来源

本文采用ASTER GDEM作为基础数据。ASTER GDEM(星载热发射和反射辐射仪全球数字高程模型)是由美国航天局(NASA)与日本经济产业省(METI)共同推出的地球电子地形数据,根据NASA的新一代对地观测卫星Terra的详尽观测结果制作完成的。其数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，达到了地球陆地表面的99%。目前，中国科学院计算机网络信息中心科学数据中心已经加工产生了中国及周边区域范围内ASTER GDEM 30 m分辨率系列数据产品。从中国科学院国际科学数据服务平台找到龙川江下游所在图幅，选择数字高程数据下载即获得所需的DEM数据。

ASTER GDEM数据可免费获取，为检验数据的可行性，对比分析了ASTER GDEM(裸露地表高程+植被和人工建筑物)与地形图制作的DEM(裸露地表高程)。由于研究区地面多为低矮植被，结果是二者的差距不明显，高程信息基本一致，提取的水系格局也基本一致。因此可用ASTER GDEM替代地形图制作的DEM提取水系特征。

2.2 研究方法

以DEM数据为基础，利用ArcGIS水文分析模块提取地表水流径流模型的水流方向、汇流累积量、河网提取、研究区流域的划分。通过对基本水文因子的提取分析，可以再现地表水的流动过程，最终实现研究区流域水文特征的数字化分析。图1简要描述了水文要素提取的基本流程。

图1 水文信息提取流程

2.2.1 水流方向提取与无洼地DEM生成 水流方向是指水流经过每一个栅格单元时的指向，由计算中心栅格与领域栅格的最大距离权落差来确定。ArcGIS中水流方向是利用D8算法[14](最大坡降法)，用邻域栅格的坡降来确定水流流向，坡降最大的单元格即为水流出的方向。

一般认为DEM是较光滑的地形表面模拟，但由于真实存在的凹陷部分(如喀斯特地貌)，可导致在进行水流流向计算时在研究区内得到不合理或错误的水流方向。因此在进行水流方向的计算前，应先对DEM数据进行填洼处理，得到无洼地的DEM[15]。洼地填充的基本过程是先计算出DEM数据中的洼地区域，然后计算出洼地区域的洼地深度，最后以最大洼地深度为参考而设定填充阈值进行研究区洼地填充，反复填洼直至无洼地DEM生成。

2.2.2 汇流累积量及河网提取 汇流累积量的基本思想是以规格网表示的DEM每点处有一个单位水量，按照自然水流由高处流向低处的自然规律，根据区域流向栅格计算每点流过的水量值，从而得到该区域的汇流累积量。当汇流量累积到一定程度的时候，地表就有水流产生，所有汇流量超过临界值得栅格就是潜在的水流路径，这些流水路径构成河网。目前常用的河网提取方法是地表径流漫流模型。计算步骤是以水流方向及汇流累积量为基础前提，设置合理的集水阈值提取河网。目前较为常用“试误法”确定集水阈值。阈值设定的大小影响河网的密度见图2。阈值越大，河道数目越少，河网密度越小；阈值越小，河道数目越多，河网密度越大。

图2 研究区河网随集水阈值变化

2.2.3 流域的分割 流域是由分水岭分割而成的汇水区域，通过对水流方向数据的分析确定出所有相互连接并处于同一流域盆地的栅格。流域的分割可以减少大尺度空间上不同因素差异导致的误差，其步骤是确定小级别流域的出水口位置，然后分析找出所有上游出水点流过出水口的栅格，结合流向数据，确定集水区栅格的位置，从而实现子流域的划分。

3 结果与分析

设置不同的集水阈值，获取的河网也不同(图2)。选取0.01～1.0 km2范围内的13个集水阈值，分别计算出对应的河网密度和子流域数量，进而得出集水阈值与河网密度、子流域数量关系图(图3—4)。

从图3—4中看出随着集水阈值的变化，河网密度、子流域数量的变化趋势大致相同。

图3 研究区集水阈值与子流域数量的关系

图4 研究区集水阈值与河网密度的关系

3.1 集水阈值对河网提取的影响

从图2可看出集水阈值取值较小时，生成的河网密集，在平地区域容易生成平行状河道(伪河道)，且生成大量破碎的子流域。随着集水阈值取值逐渐增大，河网密度和子流域数量呈下降趋势，但是降低的趋势越来越缓慢。当集水阈值由0.01 km2变化至0.1 km2，关系曲线下降明显。河网密度由0.834 km/km2减少到0.261 km/km2，子流域数量由4 686下降到464。当集水阈值从0.1 km2增至1.0 km2时，河网密度减少到0.087 km/km2，子流域数量下降到47。

3.2 集水阈值的确定

分别采用幂函数、指数、对数、线性、多项式、移动平均6种函数进行趋势拟合，发现采用乘函数进行趋势拟合时两者的拟合度最好，得到的拟合方程如下：

y=0.084 6x-0.495(R2=0.999 8)

(1)

幂函数的一阶导数、二阶导数曲线关系见图5。河网密度随集水阈值的增大无限逼近于0，即不能通过幂函数的一阶、二阶求导的方法求出转折点以提取集水阈值。因此本文提出用幂函数与直线相切的方法提取集水阈值。

图5 研究区河网密度随集水阈值变化幂函数的一阶导数、二阶导数曲线

研究中采用公式(2)为切线方程：

y=-x+a

(2)

式中：y——河网密度(km/km2)；x——集水阈值(km2)；a——截距。

幂函数与直线相切的点为拐点，拐点在数学上也称为稳定点，出现拐点之前，集水阈值变化率Δx小于河网密度变化率Δy，反之Δy<Δx；尝试求拐点即河网密度随集水阈值变化达到稳定时，得到合理集水阈值。当x有唯一解时，Δy=Δx，所取的值即为最佳集水阈值。经计算，当a=0.359时，x有唯一解0.12；当a<0.359时，x无解；当a>0.359时，x有双解。a取值0.359后，x<0.12时，Δy>Δx；x=0.12时，Δy=Δx；x>0.12时，Δy<Δx。通过上述确定集水阈值的方法与传统方法相比，避免了集水阈值取值的主观性，运用数学求值的方法简化了反复试验阈值的繁琐过程。

3.3 提取河网与真实地形对比

随着集水阈值的不断增大，河道逐渐向地势平坦处缩回，平行状河道逐渐消失。当集水阈值为0.12 km2时，地势平坦地区未出现平行线段，即伪河道全部被消除，生成的河网曲线较为光滑符合流域的河道走向。采用部分参考文献[9-10]中河网密度变化趋于稳定提取阈值的方法，可将研究区的集水阈值定为0.2 km2；在本文运用的数值方法中求得的集水阈值为0.12 km2。通过与传统方法比较，发现集水阈值为0.2 km2的河网中部分常流河道缺失，而本文提出的方法保持了常流河道的完整性，生成的河网水系与元谋县中部SPOT6遥感影像吻合。

4 结果与讨论

本文以龙川江流域为例，应用DEM和ArcGIS的水文分析模块对研究区进行水文特征提取与分析。研究区地形起伏较大，在没有高精度的DEM数据的情况下，提取的河网也能较好的与自然水系吻合。集水阈值对主干河道的空间地理位置影响小，主河流长度基本不发生变化，对研究区流域提取的河网特征影响极大。采用求解河网密度变化率等于集水阈值变化率的数值方法，确定研究区适宜的集水阈值为0.12 km2，进而得出自动提取的河网密度为0.239 km/km2，生成的河网与河网密度见图6。从河网随集水阈值变化趋势可以看出，当集水阈值较小时河网密集，存在众多小级别河道，集水阈值超过0.3 km2时小级别河道消失，主河道长度基本不发生改变，说明龙川江存在季节性河流。当集水阈值取值偏小时，提取的水文特征可用于检测季节性河流，进而对地质灾害点的提取起到积极作用。河网密度最密集的地区(图6b)分别是元谋县姜驿乡、江边乡、黄瓜园镇，这与地质灾害统计数据元谋县泥石流频发区域相吻合；当集水阈值取值偏大时，可用以识别非季节性河流的长度与地理位置。

图6 研究区河网与河网密度

不同DEM分辨率下提取的流域水系形态有差异，分辨率越高，河网越复杂精细，但河网结构大体上是一致的,提取的流域面积、河网密度相差不大[16]。随着分辨率的降低，在地势起伏较平坦地区河网特征变化大[17-18],研究区地形起伏较大，用低分辨率DEM替代高分辨率DEM，仍可取得良好的水系特征[19]。生成的河网存在部分缺失的情况与DEM数据来源有一定关系，因研究区海拔差异较大，且区内多为低矮植被，最终取得的水系特征良好。

综上所述，利用DEM和ArcGIS进行的水文特征提取具有速度快，成本低，提取效果良好的优点。通过此方法可为研究区的地形研究、水土流失防治等工作提供研究基础和参考依据。

致谢：本文得到中国科学院计算机网络信息中心提供的DEM数据，昆明理工大学地质过程与矿产资源省创新团队和昆明理工大学遥感地球化学学科方向团队联合资助，在此表示感谢！

[1] 张游,王绍强,葛全胜,等.基于GIS的江西省洪涝灾害风险评估[J].长江流域资源与环境,2011(S1):166-172.

[2] 李慧敏,张建军,黄明,等.基于DEM的黄土高原典型流域特征参数分析[J].北京林业大学学报,2012,34(2):90-95.

[3] 曾超.GIS支持下岷江上游水文特征定量分析[D].成都:四川师范大学,2011.

[4] 李晶,张征,朱建刚,等.基于DEM的太湖流域水文特征提取[J].环境科学与管理,2009,34(5):138-142.

[5] 宋向阳,吴发启,赵龙山,等.基于DEM的延河流域水文特征提取与分析[J].干旱地区农业研究,2012,30(4):200-206.

[6] Murphy P N C, Ogilvie J, Meng Fanrui, et al. Stream network modelling using lidar and photogrammetric digital elevation models: A comparison and field verification[J]. Hydrological Processes, 2008,22(12):1747-1754.

[7] Ariza-Villaverde A B, Jiménez-Hornero F J, Gutiérrez de Ravé E. Influence of DEM resolution on drainage network extraction: A multifractal analysis[J]. Geomorphology, 2015,241:243-254.

[8] 任立良,刘新仁.基于数字流域的水文过程模拟研究[J].自然灾害学报,2000,9(4):45-52.

[9] 孔凡哲,李莉莉.利用DEM提取河网时集水面积阈值的确定[J].水电能源科学,2005(4):65-67.

[10] 朱海玲,杨晓晖,张学培,等.基于DEM的密云水库上游流域特征提取与分析[J].中国水土保持科学,2013,11(3):66-72.

[11] 叶章蕊,卢毅敏,张永田.基于曲线割线斜率法的水文特征提取[J].人民黄河,2016,38(2):28-31.

[12] 关颖慧,郑粉莉,王彬,等.基于DEM的黑龙江宾州河流域水系提取试验研究[J].水土保持通报,2012,32(1):127-131.

[13] 王林,陈兴伟.基于DEM的晋江流域水系分维估算方法探讨[J].地球科学信息,2007,9(4):133-137.

[14] O’Callaghan J F, Mark D M. The extraction of drainage networks from digital elevation data[J]. Computer Vision Graphics & Image Processing, 1984,28(3):323-344.

[15] Martz L W, Garbrecht J. An outlet breaching algorithm for the treatment of closed depressions in a raster DEM[J]. Computers & Geosciences, 1999,25(7):835-844.

[16] 王培法.栅格DEM的尺度与水平分辨率对流域特征提取的分析:以黄土岭流域为例[J].江西师范大学学报:自然版,2004,28(6):549-554.

[17] 吴险峰,刘昌明,王中根.栅格DEM的水平分辨率对流域特征的影响分析[J].自然资源学报,2003,18(2):148-154.

[18] 林声盼,荆长伟, MOORE Nathan,等.数字高程模型分辨率对流域地形特征参数的影响[J].水科学进展,2012,23(4):457-463.

[19] 康婷婷,王茜,赵苗琦,等.不同比例尺水系修正DEM的流域河网提取对比[J].地理与地理信息科学,2011,27(1):111-112.

Information Extraction of River Networks and Determination of Drainage Area Threshold Using DEM Data

LUO Dayou1,3, WEN Xingping1,3, SHEN Pan1,3, ZHANG Haonan1,3, ZHANG Lijuan2, LI Jinbo1,3, YU Zhi1,3

(1.Faculty of Land Resource Engineering, Kunming University of Science and Technology, Kunming, Yunnan 650093, China; 2.Hezhou University, Hezhou, Guangxi 5428999, China； 3.Mineral Resources Prediction and Evaluation Engineering Laboratory of Yunnan Province, Kunming, Yunnan 650093, China)

[Objective] Method of determining appropriate drain area threshold was demonstrated to reflect the actual river drainage network. [Methods] Based on the DEM data of Longchuan River basin, the automatic extraction module of watershed drainage network in ArcGIS software was used to extract relevant drainage information. After that, a function between river network density and catchment threshold was fitted, which was used to determine the appropriate drainage threshold by finding the turning point of its first derivative function and /or the inflection point of the second order derivative. This was simultaneously exemplified by the Longchuan River basin. In this case, the threshold was determined by finding where the ratio change of river network density equal to the change rate of the water drainage threshold. [Results] The river networks under different catchment drainage threshold were compared and found that no obvious differences existed with regard to the catchment threshold length of the main channel and the geographic location, but other features of river network were apparently different. In the case of the Longchuan River basin, its appropriate threshold was thought to be 0.12 km2. [Conclusion] The river network extraction precision varied with the drain area threshold. The proposed method in the determination of the catchment threshold is less subjective, and is more objective, which can be used to provide some bases for the monitoring of regional soil erosion.

DEM;ArcGIS;riverdensity;drainageareathreshold

A

： 1000-288X(2017)04-0189-05

： P627

2016-12-05

：2017-02-28

国家自然科学基金项目“天空开阔度对城市热岛效应的影响研究”(41101343)

罗大游(1993—),女(汉族),四川省成都市人,硕士研究生,主要研究方向为遥感地质。E-mail:746206934@qq.com。

温兴平(1970—),男(汉族),山西省兴县人,博士,教授,博士生导师,主要从事遥感地质及应用研究。E-mail:wfxyp@qq.com。