ArcGIS中基于DEM提取沟道特征

严建钢，金复鑫，2，周小程，马海洋

（1.海军航空工程学院指挥系，山东烟台264001；2.陆军航空兵学院，北京101123）

从数字高程模型（Digital Elevation Model，DEM）生成的积水流域和河道沟系网络数据是大多数地表水文分析模型的输入数据[1]。因为DEM 中包含有各种分辨率的地形高程信息，可以提取出大量地表形态信息，如流域的坡度、坡向和流域单元格的流向等。由DEM 提取地形特征的研究在国外从20 世纪60年代就已经开始，80年代之前研究比较少，范围也仅限于分水线和山谷的识别和提取；其高峰出现在80～90年代，相继出现了各种提取河网、流域边界以及划分子流域的方法[1]。

随着流域特征提取研究方法的不断成熟，各种专业软件和计算机模型相继推出，所有这些软件和模型中，尤以ESRI 的ArcGIS 最为成熟，应用也最为普及。因而本文以ArcGIS9.3为软件平台对研究区进行流域沟道特征的提取。

1 资料与方法

1.1 研究区概况

大南沟流域地处陕北黄土丘陵沟壑区，位于安塞县城西北约7.5 km 处，属延河一级支流，流域面积约3.6 km2，海拔1 100～1 327 m。流域属暖温带半干旱气候区，年均温9 ℃，年均降雨量549 mm。流域内地形破碎，沟壑密度达6.9 km/km2；地势起伏率和坡度大，沟峁频繁相间，土壤侵蚀强烈，塑造成典型的梁状黄土丘陵地貌和复杂多样的土地类型；流域内部地势相对开阔，多平缓塌地，有大面积的梯田。流域内绝大部分土壤是黄绵土，还有少量二色土和红胶土[2]。

1.2 数据源

本文用于提取流域特征的是大南沟流域5 m分辨率的DEM数据[3]。在这里，DEM数据是通过数字化研究区1∶10 000 地形图，再插值生成TIN 数据，然后由TIN数据生成ESRI的grid格式栅格数据而成。在X方向上有580个栅格，在Y方向上有519个栅格，Z方向上有1个栅格。

1.3 提取过程

ArcGIS9.3 中没有可视化的水文分析模块，其相关的功能要通过在Grid 模块中调用一系列函数来实现。图1的流程图简要地说明了分析的过程。

图1 ArcGIS进行水文分析流程图

1.3.1 数据预处理

预处理主要是对原始DEM中存在的洼地进行填平处理。洼地是指高程低于周围栅格的一个栅格或空间上相联系的栅格的集合，洼地的出现可能是由DEM在数据采集和插值时的误差造成的，也可能是真实地形特征的反映（如喀斯特地区、采石场等）[4]。由于洼地被较高的地形包围，不能确定其水流方向，因而是进行水文分析的一大障碍，必须将其填平[5]。

在ArcGIS 中对洼地的填平处理包括以下步骤[6]：首先，基于原始DEM 数据生成流域的流向数据；然后，结合流向数据调用SINK 函数对洼地进行标定。SINK函数的功能就是标定出流域内高程低于周围的封闭区域。语句格式为SINK（＜flowdir＞），flowdir 为流向数据。

标定出洼地之后要确定洼地的影响区域，即有多少栅格的汇水到达这些洼地中。本步骤的实现调用Watershed函数，该函数需要流向数据和上一步中的洼地作为输入。语句格式为watershed（＜flowdir＞，＜sinks＞），sinks为标定的洼地栅格。

下一步要结合高程数据和影响面积数据分别确定洼地的最低高程（洼地底部）和最高高程（洼地边缘的最低高程），需要分别调用zonalmin 和zonalfill 函数。函数的语句格式分别为zonalmin（＜sink_areas＞，＜elevation＞）和zonalfill（＜sink_areas＞，＜elevation＞）。

完成上一步之后要确定需要填充洼地的最深值，即超过这一深度的洼地被认为是合理的地形特征，而不再进行填充，在此深度范围内的洼地才填充，该值用洼地最高高程栅格减去最低高程栅格即得。

最后一步就是根据以上步骤的结果对原始数据进行填洼，调用的函数是Fill。Fill 的工作原理是：扫描各个单元格的时候，比较该单元格与相邻的8 个单元格的高程，如果是洼地，那么该单元格的高程值将被赋予相邻8个单元格中高程最低的那个；Fill函数的语句格式为FILL＜in_grid＞＜out_grid＞{SINK|PEAK}{z_limit}{out_dir_grid}，in_grid 为待处理DEM 数据，out_grid 为处理后DEM 数据的路径和文件名，SINK PEAK是选择填洼或是削峰，在此处选择填洼，z_limit即为填充洼地的最深值，超过此深度不再填充，out_dir_grid为生成的新流向数据。

因为在填洼的同时改变了原有的栅格高程，可能在洼地的边缘处产生新的洼地，所以以上填洼过程要反复进行，直到用SINK 函数标定不出洼地时才可以认为填洼完成，生成了无洼地DEM数据。

一般认为，公司的盈余质量、信息披露质量、会计稳健性等越好，财务报告质量越好。目前，以审计委员会中独立董事的比例度量审计委员会独立性的研究，发现审计委员会独立性越高，内部控制质量越好（谢海娟等2016、陈文娟等2016），盈余质量越好（王振秀2017、李建红2016），会计稳健性越好（谢香兵2011），信息披露质量越好（蔡卫星等2009、柯明等2011、刘彬2014），审计意见越好（何卫红2016），财务报告质量越好（周国华2011、周国平等2013）。另外，潘珺等（2017）发现审计委员会独立董事和非执行董事比例之和越高，公司的财务报告质量越好。

1.3.2 流向数据的生成

完成原始DEM 数据的填洼预处理之后，则可基于无洼地DEM数据生成流向数据。

在ArcGIS 中生成流向数据用的是D8 算法[7]，该方法假设单个网格中的水流只有8 种可能的流向，用最陡坡度法来确定水流的方向，每个栅格单元水流的流向为其邻近8 个栅格单元中坡度最大的那个单元格，坡度的计算公式为θ=arctan[(hi-hj)/D]，其中hi和hj为2 个单元格的高程值，D为2 个单元格中心的距离，当相邻栅格为共边时，D为栅格像元的尺寸，当相邻栅格为对角时，D为 2倍的栅格尺寸。

在ArcGIS 中调用Flowdirection 函数生成流域8种方向值的流向数据，函数的语句格式为flowdirection（＜elevation＞）。图2 所示为流向栅格数据每个数值所代表的方向。

图2 栅格流向示意图

1.3.3 累积汇水面积数据生成

每个栅格单元的累积汇水面积表现了该栅格汇集上游来水的能力，汇水能力越大的栅格就越可能是河道。在ArcGIS中，这一过程被称为求栅格的累计流量[8]，通过调用Flow accumulation 函数来实现。此功能的原理是假想在集水区的每一网格上降下一单位的水量，而后按网格的流向来向下移动，其移动经过的网格则使其累积流量值提升一个单位，因而每一网格都能计算出其所累积的上游流量值。由于投入每一网格的水量皆为一单位，故流量累积值亦代表各网格的上游集流网格数量，将之乘上网格面积便可得到每一网点的上游集水面积[9]。函数的语句格式为flow accumulation（flowdir）。

1.3.4 临界汇水面积的确定和沟道提取

临界汇水面积，即CSA（Critical Source Area），是区分河道栅格和流域内非河道栅格的临界值，大于该数值的栅格被赋值为1（沟道），小于该数值的赋值为NULL。这一值在不同的气候带和土壤覆盖条件下是不同的，如果取值太小，则会在提取中出现大量伪沟道，如果取值太小则会忽略真实水系。国内外很多学者对CSA的取值进行了研究，有些学者认为在地形复杂区应根据地貌特征选择不同的值，国内一些学者采取与光照晕渲图拟合的方法确定CSA取值[10]，本文即采用此种方法来确定临界值。

确定出汇水临界面积之后，要构造条件语句进行流域汇水能力数据的二值化，即流域沟道的提取。在Grid 中调用con 函数，语句格式为con（＜condition＞，＜true_expression＞，{＜condition＞，＜true_expression＞}，…{＜condition＞， ＜true_expression＞}， {false_expression}）。其中{false_expression}为可选项。可选项为空时，如果＜condition＞语句返回的值为假，那么单元格将被默认赋值为NULL。例如：con（flow_accumulation＞Const，1）表示将栅格数据flow_accumulation 中的汇流的值大于某一常数的全部赋值为1，其余赋值为NULL。

1.3.5 沟道的分级

对沟道的分级有2 种方案[10]：Strahler 水系系统分级方案和Shreve 水系系统分级方案。前者的原理是以源头沟道开始为第1 级，相同级别的沟道相交后级别增1级，不同级别沟道相交后级别不增，保持原沟道中较高的级别，依次向下游追踪；后者考虑河网中所有的结点，最外层沟道赋级别为1，向下游追索到的结点级别为上游级别的累加。分级要调用streamorder函 数，语句格式为streamorder（streamnet，flowdir，strahler/shreve），streamnet 为河道数据，flowdir 为流向数据。strahler/shreve为提供备选的分级方案，本文选择Strahler方案。

2 结果与分析

在本次研究中自动统计生成的关键的数据有3个：填洼之后生成的无洼地DEM 数据，基于无洼地DEM数据生成的流向数据和流域累计汇水面积数据。

2.1 无洼地DEM

在预处理中通过扫描原始DEM 数据标定了5 处洼地区域，最低的洼地底部高程为1 145 m，最高的洼地底部高程为1 195 m。计算得出的填充洼地的最深值为1，以此值为限制用Fill 函数进行填充，对填充后的DEM高程数据用SINK函数查找洼地得到空值，说明流域中影响流向确定的洼地已经被移除，见图3。

图3 流域无洼地DEM图

2.2 累积汇水面积分析

流域累积汇水面积数据（见图4）是进行沟道提取的基础，本文结合生成的地形光照晕渲图来设定合适的CSA 值，经过多次选值与晕渲图叠加分析，认为当CSA 取300 时，所提取的沟道与晕渲图的地形特征吻合较好（图5 白色线状部分为提取的沟道），可以把二值化得出的数据作为最终的沟道提取结果。

图4 流域累积汇水面积图

图5 沟道与地形晕渲图叠加示意图

2.3 沟道分级结果

通过对提取出的流域沟道数据进行Strahler 方案水系分级，可将流域内的沟道系统分为4级，其中从源头开始的1 级水系有58 条，一级水系交汇而成的2 级水系有15 条，3 级水系有3 条，4 级水系有1 条。水系的分级基本符合流域沟道的实际分布特征，见图6。

图6 沟道分级图

3 结论

通过与大南沟小流域1∶10 000地形图所反映的实际流域特征对比分析，可以认为提取出的沟道特征基本符合实际，即用ArcGIS的水文分析功能可以进行流域沟道特征的提取，这在水文分析和沟道系统研究领域有着重要的意义。

但在分析中也存在一些问题有待于做进一步的研究和探讨，主要有以下几方面：

1）填洼预处理过程有可能会导致实际地形的改变。在ArcGIS 中把洼地的成因归为是由于对地形的过低估计而致，故而采用增加高程的填平方案，而实际上洼地还可能是由于周围高程过高估计而成，或是（有坡度地区）被高地阻挡而成（阻挡型洼地），单纯的填平会导致这类洼地实际地形被改变。

2）对平地的处理效果不明显，在地形起伏小的地区，沟道的模拟会出现误差，尤其是在填平处理后产生的平坦地区，流向的确定仍不理想。

3）CSA值的选取，本文采用地形光照晕渲图做参照的方法来确定临界汇水面积值，但这种方法的主观性太强，缺乏定量的数学基础，且光照晕渲图只能模糊地反映地形起伏特征，以其为参照会产生误差，故需要进一步讨论更好的方法，以使提取的沟道能正确地反映实际状况。

[1] 李丽，郝振纯.基于DEM的流域特征提取综述[J].地球科学进展，2003，4（18）：251-252.

LI LI，HAO ZHENCHUN. The automated extraction of catchment properties from digital elevation models[J].Advance in Earth Sciences，2003，4（18）：251-252.（in Chinese）

[2] 许明祥，刘国彬，温仲明，等.黄土丘陵区小流域土壤特性时空动态变化研究[J].水土保持通报，2000，20（1）：21-23.

XU MINGXIANG，LIU GUOBIN，WEN ZHONGMING，et al. Temporal and spatial variation of soil characters in small catchment of loess hilly areas[J]. Bulletin of Soil and Water Conservation，2000，20（1）：21-23.（in Chinese）

[3] 赵健，贾忠华，罗纨.ARCGIS 环境下基于DEM 的流域特征提取[J].水资源与水工程学报，2006，17（1）：74-76.

ZHAO JIAN，JIA ZHONGHUA，LUO WAN. Extracted catchment properties from DEM using ARCGIS[J]. Journal of Water Resources & Water Engineering，2006，17（1）：74-76.（in Chinese）

[4] 肖飞，张百平，凌峰，等.基于DEM的地貌实体单元自动提取方法[J].地理研究，2008，27（2）：459-466.

XIAO FEI，ZHANG BAIPING，LING FENG，et al. DEM based auto-extraction of geomorphic units[J]. Geographical Research，2008，27（2）：459-466.（in Chinese）

[5] JENSON S K，DOMINGUE J O. Extracting topographic structure from digital elevation data for geographical information system analysis[J]. Photogram Metric Engineering and Remote Sensing，1988，54（11）：1593-1600.

[6] 张超，郑钧，张尚弘.ArcGis9.0中基于DEM的水文信息提取方法[J].水利水电技术，2005，36（1）：1-4.

ZHANG CHAO，ZHENG JUN，ZHANG SHANGHONG.Extraction of hydrological information from digital elevation model with ArcGis 9.0[J]. Water Resources and Hydropower Engineering，2005，36（1）：1-4.（in Chinese）

[7] 孙庆艳，余新晓，胡淑萍，等.基于ArcGIS 环境下DEM流域特征提取及应用[J].北京林业大学学报，2008，30（2）：144-147.

SUN QINGYAN，YU XINXIAO，HU SHUPING，et al.Extraction and application of hydrological information based on DEM in ArcGIS environment[J].Journal of Beijing Forestry University，2008，30（2）：144-147.（in Chinese）

[8] 李翀，杨大文.基于栅格数字高程模型DEM的河网提取及实现[J].中国水利水电科学研究院学报，2004，2（3）：208-214.

LI CHONG，YANG DAWEN. Deriving drainage networks and catchment boundaries from Grid Digital Elevation Model[J]. Journal of China Institute of Water Resources and Hydropower Research，2004，2（3）：208-214.（in Chinese）

[9] 易红伟，汤国安，刘咏梅，等. 河网径流节点及其基于DEM 的自动提取[J].水土保持学报，2003，19（3）：113-118.

YI HONGWEI，TANG GUOAN，LIU YONGMEI，et al.Stream runoff nodes and their derivation based on DEM[J].Journal of Soil and Water Conservation，2003，19（3）：113-118.（in Chinese）

[10]李金朝，国庆喜，葛剑平.基于DEM 的黄土高原沟壑区沟道系统的自动提取[J]. 西北林学院学报，2009，24（6）：220-223.

LI JINZHAO，GUO QINGXI，GE JIANPING，et al.DEM based automated extraction system of gully in the loess upland gully are[J].Journal of Northwest Forestry University，2009，24（6）：220-223.（in Chinese）