基于数字高程模型的河网和流域边界提取

贾岸斌，佘高杰 ，卢 舟，蒯星龙

(1.常德市气象局，湖南常德 415000；2.桃源县气象局，湖南桃源 415700)

河网和流域边界的描述是水文气象相关业务与研究的前提，是做好流域面雨量预报预测及预警服务等各项工作的前提条件[1]。传统上(甚至现在)，这项工作一般是通过地形图或等高线图手工绘制完成的[2]。随着遥感、GIS等行业及数字高程模型(Digital Elevation Modell，简称DEM)的发展，为数字化提取水文信息提供了可能。周成虎、杜云艳利用AVHRR影像，提出了基于水体光谱的水体自动提取识别模型[3]。DEM本身不包含河流、湖泊水库、堤坝等信息，通过DEM直观判读流域水系及边界较遥感影像困难，需要借助算法实现。目前最常用的是基于地形表面流水物理模拟分析算法中的D8算法。利用遥感影像识别水体和DEM提取流域边界无论在精度上还是效率上都高于传统地形图手工勾绘方法。遥感影像具有获取成本低、覆盖率高、更新快、现势性强等优点，但也存在云、积雪覆盖、山体阴影等不利影响[2]。DEM在山地丘陵区和平均地形坡度不小于3°的区域所生成河网具有很高的可靠性，而在平均地形坡度小于3°的平坦区域河网生成中产生的虚拟河网与自然水系偏差较大[4]，但DEM还可以很方便地提取其它流域参数，仍然被广为使用。本文以沅水流域为例，尝试利用数字高程模型制作沅水主要河流的河网及流域边界信息化文件，并将其应用于流域面雨量预报服务工作，希望能对气象部门开展以流域为对象的相关预报服务及研究工作，比如针对河(湖)长的预报服务[5]等，起到一定的帮助。

1 研究区域与数据资料

沅水，又称沅江，位于26°N～30°N、107°E～112°E之间，属于长江流域洞庭湖水系，干流全长1 033 km，流域面积89 163 km2，跨湘黔渝鄂四省(市)。沅水自河源至黔城为上游，多深山幽谷，黔城至沅陵为中游，为丘陵地区，沅陵以下称下游，桃源以下为冲积平原[6]，如图1所示。

覆盖研究区域的DEM数据，用于水文特征的提取。可供选择的DEM数据源较多，如GTOPO30 (https://lta.cr.usgs.gov/GTOPO30)、SRTM DEM(http://srtm.csi.cgiar.org)、ASTER GDEM (http://gdem. ersdac.jspacesystems.or.jp/)、World DEM (http://www.astrium-geo.com) 等[7]。可从网格分区中选取目标区域的Geo TIFF文件。若目标区域涵盖几个分区，则需做合并处理。本文选取美国奋进号航天飞机的雷达地形测绘SRTM(Shuttle Radar Topography Mission)数据，目前最新为V4.1版本[8]，分辨率3″(约90 m)，高程垂直误差小于16 m，水平分辨率高于现行全国智能网格预报(1～5 km)，提取数据精度能达到日常业务应用要求。

2 主要算法原理及关键技术

在GIS环境下，基于DEM的河网及流域边界等水文信息的数字化提取流程如图2所示，主要有DEM预处理，流向、汇流累积量计算等水文分析。最后将栅格数据矢量化，获得GIS交换文件。

图1 沅水流域地理位置及地形示意图

图2 基于DEM的数字化河网及流域边界提取流程

2.1 DEM数据预处理

DEM数据预处理主要包括投影、拼接、裁剪与填洼。

一般DEM数据，只定义了地理坐标，无投影坐标，需要从Project工具，选择WGS_1984_World_Mercator将原地理坐标GCS_WGS_1984投影到平面上。否则，虽不影响提取，但长度、面积属性将无法使用。拼接(Mosaic To New Raster)之前，先导入一个DEM数据文件，查看相关属性，方便拼接时填写。主要有像素类型(16_BIT_SIGNED)及色带数量(1)。DEM裁剪工具位于Spatial Analyst Tools-Extraction-Extract by Mask处。

一般用原始DEM数据计算流向(Flow Direction)，可以发现结果并不是期望的8个预定方位，而是被赋值1～255。因为有汇的存在，流向被赋予了所有可能流向的和。汇查找工具(Sink)可以计算得到所有汇。可以认为所有汇，均不是自然存在[9]，而是数据原因导致，全部需要填充。填洼Fill，输入原始DEM作为参数。可选参数z limit不填，即默认全部填充。完成填洼运算后，得到无凹陷DEM数据。

2.2 水流方向计算

最先出现的流域特征提取方法是1975年Peuker等提出的识别谷点法[10]，但是该方法即便经过(Band和Donglas)改进仍然不能适用于大多数地貌类型[11]。目前应用最广泛的是由奥克拉芳(O’Challaghan)和马克(Mark)两人于在1984年提出的无限坡面流累计方法—D8算法[12]。D8算法是最早出现的经典的单流向算法。该方法对自然状态的水流方向进行了极大的概括，计算简单，可操作性强。算法假设单个网格中的水流只有 8 种可能流向，即流入与之相邻的8个网格中。它用最陡坡度法来确定水流的方向，在3×3的 DEM网格中，计算中心网格与各相邻网格间的距离权落差(即网格中心点落差除以网格中心点之间的距离)，取距离权落差最大的网格为中心网格的流出网格，该方向即为中心网格的流向。流向分析如图3所示。可以发现若定义网格距为1单位，则与其正交相邻4个单元格中心距离也为1单位，与其对角相邻4个单元格中心距离则为21/2，约1.414 21。

图3 方向编码及流向

根据流向分析结果，可以计算每个单元格的上游汇流能力，如图4所示。依据需要(如集水面积)设定汇流能力阈值，不低于该阈值的单元格标记为河谷。不难发现，最大值处应为该流域出水口。该方法简单，可直接产生连续的河网。

流向计算采用Flow Direction 工具。以填洼后的无凹陷DEM为输入参数，可以选择同时输出drop raster，即D8算法中的距离权落差。这一步可以检验填洼是否完全。正常情况，应得到预定的8个方位编码值，结果只能是1，2，4，……128，如若结果范围为1～255，则证明填洼不完全，需要重复填洼步骤，直到结果为预定8方位编码为止。

图4 基于流向的汇流累积量分析

2.3 汇流累积量计算

采用Flow Accumulation工具，将流向结果作为参数输入，权重参数默认相同。因为河网分级与流域划分，一般按照一定的控制面积来确定，而像元大小均一。若通过径流量的大小来确定，则需要设定每个网格的权重参数。因为每个网格的降水量和地质地貌等条件都将影响该网格最后的径流量，全区域并不均一。另外，此部分计算相对较慢，尤其当DEM精度较高且研究区域较大时。

2.4 河网提取分级与矢量化

提取河网采用基于汇流累积量设定阈值的方法。阈值依据实际需要设定。例如网格分辨率为90 m×90 m，提取流域面积100 km2以上的河网数据，则汇流累计量需要不小于100×106/(90 m×90 m)个网格。通过Raster Calculator工具，SetNull("FlowAcc"<12345,1)，将河谷栅格设置为1，非河谷栅格设置为空(也可使用Con 工具)，得到河网栅格。河网连接Stream Link和河网分级Stream Order分别为每段河流分配唯一编号与标定相应级别。

河网栅格矢量化，Stream to Feature可以将设定汇流累积量后的河网，编号后的河网连接以及分级河网等栅格数据转化为矢量数据，以shapefile格式存储，方便交换使用。

另外栅格矢量化Raster to Polyline工具，也可以矢量化河网，区别是前者主要针对河流网络矢量化或者说任何一类方向明确的代表线条网络的栅格数据的矢量化，经过优化。表现在基于方向数据的辅助下，该算法将相同值的两个相邻线条特征矢量化为两条平行线，而后者会将两条线条重叠到一起，建议选择Stream to Feature工具。

2.5 流域划分与矢量化

流域提取工具有两个，分别是流域盆地(Basin)和集水区(Watershed)划分，主要区别是Basin从全局数据范围划定流域，而Watershed可以通过设定出水口位置的方法，提取任意(子)流域。若无出水点位置栅格或者矢量数据，可以利用已经生成的Stream Link作为汇水区的出水点，此方法因出水点较多，将得到当前条件下最细小的子流域划分。最后流域面可以通过矢量化工具Raster to Polygon转为shapefile文件存储。

2.6 提取检验

为了检验提取效果，将沅水流域面积[13]与其对应提取面积对比，如表1，可以发现沅水流域面积误差率在0.71%，其主要支流提取误差除武水、酉水偏差4%左右略大外，其他均在±1%内。

表1 沅水流域面积提取误差

3 分级河网及流域文件在气象中的应用

不同于单一图形文件，利用数字化的河网和流域边界信息，可制作GrADS掩膜mask文件。设定子流域编码为1～n，在全域范围内按照气象数值模式网格，将所有格点赋值其所在子流域编码值，流域外则赋值为0。基于此特定mask文件，可方便利用GrADS函数(aave(maskout(YS,0.5-abs(mask.2-YsSubCode)),g))计算所有子流域面雨量，利用Draw String和Draw shp绘制填值、填色图。图5是欧洲中心细网格模式预报沅水流域未来12～84 h降水量的原始格点值与流域面雨量统计图。

图5 欧洲中心2017-05-19T08模式预报沅水流域未来12～84 h降水量(a 模式格点降水填色，b 子流域面雨量统计)

可以发现预报降水呈北多南少的态势，沅水流域大部分地区有中雨，南源龙头江与南侧支流渠水、巫水为小雨，各子流域统计结果与模式原始格点数据相符。子流域统计后的预报结论，汇流关系清晰，更利于防汛抗旱与水库调度工作的开展。此方法可直接使用模式原始数据，充分利用了GrADS内部函数，无需另行统计各子流域面雨量，且简便可行，符合气象行业标准[14]。

4 结论与展望

本文以沅水流域为例，通过ArcGIS软件，运用D8算法，详尽介绍了其中原理，展示了相关流程与实现步骤，并结合实践中可能遇到的问题，重点提出了几点细节。此方法，相对纯手工绘制简便快捷。提取的河网及流域边界文件，可直接应用于MICAPS、SWAN等气象业务系统。结合特定mask文件，还可利用气象常用绘图软件GrADS的内部函数，简便的应用于流域面雨量预报的图形展示。由于DEM精度等各方面原因，提取结果与实际情况存在一些细微出入，尤其在平原地区，因高程差小，相对误差大，容易产生平行河流。充分利用高分辨率DEM数据结合卫星影像的比对订正，是提高水系流域提取精度的有效方法。