ArcGIS基于DEM提取微流域单元完整实操指南

前阵子做项目需要把一片山区划分为一个个相对独立的汇水单元本来想直接用现成的子流域边界数据结果要么精度太粗、要么覆盖范围对不上折腾半天还是决定自己从DEM算一遍。用ArcGIS基于DEM提取微流域单元这套流程跑通之后回头整理了一下发现核心逻辑其实就一句话让水流在栅格DEM上“走”一遍然后把每条水流路径上游的集水区切出来。这篇文章就把完整流程拆开讲清楚从填洼、流向、汇流累积、河网提取到最终的微流域矢量面每一步我都会说明“为什么这么做”“参数怎么定”“踩过什么坑”并给出练习数据的准备方式。适合GIS相关专业的学生、做水文分析或面源污染估算的从业者以及所有想把一片复杂地形快速切成一个个“小流域管理单元”的人。1. 先搞清楚微流域单元到底是什么为什么要自己做1.1 三个概念怎么分流域、子流域、微流域流域这个概念大家都不陌生简单说就是降水落到地面后顺着地形汇流最终都流向同一个出水口的整个区域。子流域则是在大流域内部沿着次级水系、支沟继续划分出来的相对独立的汇水区。而微流域是更往下一级的空间单元往往以某条小支沟或者某个小汇水洼地为边界面积可以从几公顷到几十平方公里不等。在实际项目中微流域单元之所以好用是因为它足够小、足够完整。比如做面源污染估算时如果拿整个大流域来统计上游和下游的土地利用差异太大污染物负荷被平均掉了结果根本没参考意义。但切成一个个微流域单元后每个单元内部的土地利用、坡度、土壤类型相对均一统计出来的数据才真正能支撑决策。水土保持规划、生态水文模拟、山地灾害评价等领域也都有类似需求。1.2 从DEM到微流域的原理一条水流走出来的逻辑链DEM也就是数字高程模型本质上是一张每个栅格都带有高程值的“数字地形沙盘”。ArcGIS提取微流域的整套逻辑就是在这张沙盘上模拟水流运动。首先需要确定每个栅格的水往哪个方向流这就是流向计算。ArcGIS默认使用D8算法也就是把每个栅格周围8个相邻栅格中坡度最陡的方向作为水流方向。水流方向确定后再统计每个栅格到底“接收”了多少上游栅格汇入的水这就是汇流累积。汇流累积值大的栅格往往就是沟道、河谷的位置。基于这些沟道位置可以提取出河网再以每条独立河段作为出水口把上游所有汇入该河段的区域划出来就得到了一个个微流域单元。可以这样理解把整个流域看成一座教学楼微流域单元就是一间间独立教室。雨水落在每间教室的屋顶通过各自的排水沟汇入走廊最后从大门排出。水流路径就是“屋顶—排水沟—走廊—大门”的层级关系ArcGIS要做的就是把每一间教室的范围精确圈出来。1.3 哪些场景非用微流域单元不可我自己接触过的项目里微流域单元主要用在三类场景。第一类是面源污染估算。上游是农田还是林地、施肥量多少、坡度多大这些因素直接影响污染物入河量。把研究区切成微流域单元后可以按单元聚合土地利用和污染源数据再结合降雨径流系数算各单元的污染负荷。第二类是水文模型的汇水区准备。SWAT、HEC-HMS这类模型需要先把流域划分成子汇水区。虽然模型自带划分工具但直接用ArcGIS提取好的微流域单元导入往往更可控尤其是当你对单元大小和边界有特定要求时。第三类是水土保持和生态评价中的“小单元治理”。每条小支沟就是一个治理单元治理前后做对比监测也需要以微流域为统计口径。这时候如果手头没有现成边界基于DEM自己提取是成本最低、周期最短的方案。2. 开工前准备DEM数据源、坐标系与练习数据2.1 常见公开DEM数据源怎么选提取微流域单元之前得先有一份可靠的DEM。不同数据源的空间分辨率直接决定了微流域划分的精细程度选数据时不能只看“能用”还得看“够不够用”。下面是我用过的几个公开DEM数据源对比数据源空间分辨率覆盖范围获取途径适用场景SRTM 1弧秒30米全球60°N~56°S地理空间数据云、USGS EarthExplorer中小流域、区域级分析ASTER GDEM V330米全球83°N~83°S地理空间数据云高纬度补充、无SRTM区域ALOS PALSAR12.5米全球ASF DAAC小流域精细分析Copernicus GLO-3030米全球Copernicus官网、AWS欧洲及全球通用无人机LiDAR点云生成的DEM0.5~2米局部自有航测或测绘单位精细微流域、小地块做微流域单元我个人建议优先选12.5米或更高精度的DEM。分辨率越低地形细节丢失越多原本清晰的小支沟可能直接“消失”在栅格里划分出来的微流域边界也会比较粗糙。如果只是练手或者做大区域的初步分析30米SRTM也够用后面实操步骤完全一致。这里顺便提一句DEM和DSM不是一个东西。DSM是地物表面模型包含了建筑、树木的高度做水文分析时会造成虚假的阻挡和水流中断所以一定不要用DSM代替DEM。从开源数据里下载时注意看清产品类型。2.2 坐标系这个坑别等算面积时才后悔这是整个流程里最容易被忽略、但后果最严重的一步。绝大多数公开DEM原始数据都是WGS84经纬度坐标系单位是度直接用这种数据做水文分析流程能跑通但到了最后计算面积、长度时结果全是错的因为经度1度在不同纬度代表的实际距离不同栅格像元根本不是正方形。正确的做法是先做投影。用ArcGIS的Project Raster工具把DEM从WGS84投影到适合研究区的投影坐标系。中国地区我常用UTM分区坐标系比如UTM 50N覆盖东经117°至123°区域做省级以上大范围分析则建议用Albers等积投影。判断标准只有一个投影后栅格像元的单位必须是米而且研究区内面积变形要尽量小。这个坑我在刚接触水文分析时踩过当时提取完微流域后顺手算面积结果一个明显几千亩的小流域属性表里面积却只有几百检查了半天才发现是坐标系惹的祸。所以现在每次拿到DEM第一件事就是看属性里的坐标系绝不嫌麻烦。2.3 练习数据准备与软件版本说明这篇实操使用的练习数据我建议你从“地理空间数据云”免费下载网站提供SRTM和ASTER GDEM两种常用DEM注册后按区域框选下载即可。练习区选一块地形起伏明显的山地面积大约20km×20km内部包含山脊、沟谷、坡面这样提取出来的微流域单元层次分明方便对比效果。软件部分ArcMap 10.2到10.8都可以ArcGIS Pro同样适用水文分析工具组的位置和参数基本一致。我用的是ArcMap 10.8版本整个流程在ArcGIS Pro 3.x里也跑通过区别只是界面布局略有不同工具名称和参数含义没变化。下载好DEM后先在ArcMap里加载打开图层属性确认坐标系。如果是WGS84经纬度坐标先执行Project Raster投影到UTM。如果已经是投影坐标系且单位是米就可以直接进入下一步了。3. ArcGIS基于DEM提取微流域单元完整实操3.1 填洼最容易被忽略却影响全局的第一步在开始水文分析之前几乎所有教程都会强调先做填洼处理但很多人只是机械执行并不理解为什么要填。DEM虽然是数字地形模型但受数据采集精度和地表真实形态影响里面往往存在一些“伪洼地”——也就是周围栅格都比它高、水流进去就出不来的低洼格。如果不把这些洼地填平水流方向计算到这里就会中断提取出的河网会支离破碎后面的微流域划分自然一塌糊涂。操作路径是ArcToolbox → Spatial Analyst Tools → Hydrology → Fill。输入DEM输出结果命名Fill_DEM关键参数Z limit可以留空也可以自己设定。留空代表把所有洼地全部填平是最省事的做法设定Z limit则代表只填深度小于该值的洼地保留真实的自然洼地。我个人的经验是在喀斯特地貌或者有明显人工水库、天然湖泊的区域不能无脑全填。有一次我在喀斯特区域做分析溶蚀洼地本身就是真实地形的一部分全填之后地形几乎被“磨平”了微流域边界完全偏离实际。后来设置Z limit为50米只填掉小于50米的伪洼地结果才合理。普通山地练习数据一般不会遇到这个问题直接用默认全填即可不影响学习流程。3.2 流向计算与汇流累积搞懂数据在“算”什么填洼完成后接下来就是流向计算。运行Flow Direction工具输入Fill_DEM输出方向栅格。这个栅格每个像元的值是一个数字编码对应D8算法确定的8个方向之一比如1代表向东流128代表向西流等等。ArcGIS默认使用D8算法也是水文分析最经典的算法它假设每个栅格的水只流向周围8个邻域中坡度最陡的那一个。方向栅格出来后接着做汇流累积。运行Flow Accumulation工具输入方向栅格输出结果里每个像元的值代表有多少个上游像元的水汇流经过这个位置。可以这么理解在分水岭附近的山脊线上汇流累积值往往是0或者很小的数字因为那里是水流的起点而在沟道和河谷里汇流累积值会非常大。实操中我习惯在得到汇流累积结果后先右键查看符号系统做一次拉伸显示或者分位数分类大致了解数值分布区间。如果最大汇流累积值只有几百说明研究区不大如果有几万甚至几十万说明这是一个比较大的流域。这个数值能直接帮你判断下一步河网提取的阈值范围不用靠猜。3.3 河网提取阈值怎么定直接影响微流域大小河网提取是整条流程中最需要“手感”的一步。ArcGIS并没有单独一个“抽取河网”的按钮而是用栅格计算器对汇流累积结果做条件筛选当像元值大于某个阈值时就认为是河道赋值为1否则为NoData。具体操作是打开Spatial Analyst工具里的栅格计算器输入表达式以我常用的图层名为例Con(FlowAcc 500, 1)这里的500就是阈值。阈值越小被判定为河道的像元越多提取出的河网越密集最终划分出的微流域单元也越多、越小阈值越大河网越稀疏微流域单元越少、越大。到底选多少取决于你想要的微流域尺度。我一般会做2~3组对比测试。比如在约50平方公里的山地小流域里阈值500提取出来的是主沟道和几条大支沟微流域单元面积相对较大阈值100时细小的坡面沟道也进来了微流域数量陡增单元会非常碎。练习阶段建议从较小的阈值开始试比如200或300然后叠加卫星影像或者地形图对比看提取的河网是否与真实沟道吻合。这个调整过程没有标准答案完全是基于研究目的和对地形判断的经验活。这里有一个小技巧栅格计算器里的图层名如果包含空格或特殊字符必须用双引号引起来。很多新手在这里报错多半就是图层名没加引号或者引号用成了中文全角。3.4 用河网链接当倾泻点切出微流域单元河网提取完成后下面这一步才是微流域单元划分的核心。这里要用到两个工具Stream Link和Watershed。先运行Stream Link工具输入河网栅格和流向栅格输出“河网链接”栅格。这个栅格给每条独立的河段分配了一个唯一编号每条河段代表一段连续的河道。接着运行Watershed工具输入流向栅格倾泻点数据选择上一步得到的河网链接栅格输出结果就是微流域栅格。为什么用河网链接而不是直接用河网栅格作为倾泻点我最早做的时候图省事直接拿河网栅格去划分流域结果每个河道像元都被当成了出水口划出来的区域碎得像马赛克完全没法用。而河网链接把每条河段作为整体单元划分出的每个区域正好是这条河段的直接汇水区这才是“微流域单元”的意义所在。这一步得到的流域栅格值就是河段编号。同一编号的区域就是对应河段上游的集水区。栅格里每个独立区域代表一个微流域单元单元数量等于河段数量相互之间不会重叠边界沿着分水岭走效果非常直观。3.5 栅格转矢量与后处理让结果能直接用微流域栅格虽然已经能看出划分效果但实际项目中通常需要矢量面用于制图和空间分析所以还要把栅格转为面要素。运行Raster to Polygon工具输入微流域栅格不要勾选简化面选项输出得到一个矢量面图层属性表里带一个Gridcode字段代表每个微流域单元的编号。栅格转矢量只是第一步后处理才是让结果能真正交付的关键。我一般会做以下几步第一用Dissolve工具以Gridcode字段为融合字段把相邻的同编号图斑合并成一个完整面。有时候栅格转面会把一个单元拆成多个碎片融合能解决这个问题。第二在属性表里新增Area字段右键Calculate Geometry计算每个微流域单元的面积。注意这一步之前必须确保图层是投影坐标系单位是米否则面积就是错的。这也是我在2.2节反复强调坐标系的原因。第三筛选删除边缘碎片。研究区边界处会出现一些面积特别小、形状不完整的图斑用Select By Attributes选出面积小于某个阈值的图斑根据实际情况删除或者合并到相邻单元中。第四如果对边界平滑度有要求可以用Simplify Polygon或Smooth工具做轻度平滑。但不要平滑太狠否则边界偏移会导致面积失真。对于绝大多数分析场景不平滑也完全可以接受。4. 实操避坑常见问题与排查方法4.1 填洼填过头地形全“平”了怎么办症状填洼后的DEM在三维显示或山体阴影下看着特别“平”真实谷地、盆地特征消失了提取的微流域边界也跟实际地形对不上。原因Z limit设置过大或者把真实的自然洼地也一并填掉了。解决办法查看原始DEM中洼地的深度分布。可以在填洼前先用Zonal Fill或焦点统计对比原DEM与填洼结果的差值栅格找出哪些区域被填得最深。如果这些区域恰好是真实的喀斯特洼地、采石坑或者水库库区就需要重新执行Fill工具设置一个合理的Z limit只填掉浅层噪声造成的伪洼地。4.2 提取的河网断断续续、不连贯症状河网栅格里有大量断点河道不连续划分出来的微流域也是七零八落。原因最常见的原因是原始DEM存在NoData区域或者填洼步骤没有执行。NoData会把流向计算“截断”水流到无效值边界就停了后面的河网自然不连续。解决办法填洼前先检查DEM是否有NoData区域有的话先用栅格计算器把NoData替换成周围有效高程值或者用按掩膜提取裁掉无效边界。然后确保填洼步骤执行成功再重新计算流向和汇流累积。另外一个容易被忽略的原因是DEM范围太小流路还没汇成河就出了边界这种情况需要扩大DEM范围。4.3 微流域碎片太多、面积忽大忽小症状划分出来的微流域单元数量巨大很多单元面积只有几百平方米有些又有几十平方公里两级分化严重。原因河网提取阈值太小把坡面细沟都当成了河道或者直接用河网栅格而不是河网链接作为倾泻点又或者栅格转矢量后没有做融合和碎片筛选。解决办法调大河网提取阈值重新提取河网确认Watershed工具用的是Stream Link的结果转矢量后执行Dissolve融合并按面积筛选删除边缘碎片。如果希望单元更均匀可以考虑用Stream Order工具对河网分级只保留特定级别的河段作为划分依据这样可以人为控制微流域的尺度。4.4 倾泻点跑偏出水口不在河道上症状自定义倾泻点做流域提取时得到的流域面积特别小甚至只有一个栅格。原因倾泻点位置没有落在河道像元上而是落在坡面上。水流方向从坡面汇入河道需要一段距离如果倾泻点本身不在河道内Watershed工具只能找到该点正上方极小的汇水区域。解决办法使用Snap Pour Point工具把倾泻点捕捉到指定范围内汇流累积值最大的像元上。设置合适的捕捉距离比如100米或200米确保点被吸到河道中心。实操中我习惯在捕捉前先加载汇流累积栅格目视确认点的位置是否与高值沟道吻合。4.5 面积计算结果不对先查坐标系症状计算出的面积数值明显偏小或者偏大和常识对不上。原因图层坐标系是WGS84经纬度或者栅格转矢量后坐标系信息丢失。解决办法查看图层属性里的源坐标系。如果是地理坐标系需要投影变换后再计算面积如果坐标系信息显示为Unknown需要用Define Projection工具重新定义坐标系再做面积计算。这类问题排查顺序永远是坐标系优先不要先怀疑工具出了问题。5. 微流域单元做完之后还能干什么5.1 水文模型汇水区准备微流域单元最常见的后续应用是水文模型的汇水区输入。SWAT模型导入子流域边界时可以直接使用ArcGIS提取的微流域矢量面每个面作为一个水文响应单元的基础。配合坡度、土地利用、土壤类型数据就可以构建完整的产汇流模拟框架。我实际做过的一个项目里把研究区划分出127个微流域单元导入HEC-HMS后按单元设置降雨参数和糙率系数模拟精度明显好于直接用单一流域出口做总量计算。原因很简单微流域单元能体现空间异质性把降雨和地表条件差异保留在模型里而不是平均掉。5.2 叠加统计分析把“地形单元”变成“管理单元”提取完微流域单元下一步通常是叠加其他专题数据做统计。比如把土地利用类型和微流域边界做相交分析计算每个单元内部的耕地占比、林地占比把坡度数据按微流域做分区统计算出每个单元的平均坡度和最大坡度或者把降雨量插值结果按微流域聚合得到每个单元的面雨量。这些统计结果会自动成为微流域属性表里的字段后续做生态评价、污染负荷估算、治理优先级排序时直接按字段筛选排序就可以了。这时候你会真正感觉到微流域单元已经从“地形产物”变成了“管理工具”。5.3 进阶结合ArcGIS Pro与自动化批处理如果你用的是ArcGIS Pro水文分析工具组里还有不少增强功能比如可以更方便地做参数设置和结果符号化。更进阶的用法是使用ArcGIS Pro的模型构建器或者Python地理处理脚本把填洼、流向、汇流累积、河网提取、流域划分这条链路封装成一个自动化流程。这样不同区域、不同DEM数据只要更换输入文件就能一键跑完批量处理效率会高出很多。我后来在处理多个县区的微流域划分时就是写好一个Python脚本循环处理每个区域的DEM输出矢量微流域边界。整个过程从手工一个个点工具变成了自动化批量运行不仅省时间还避免手动操作导致的参数不一致问题。最后说一点个人体会。微流域提取这套流程本身并不复杂难的是每一步都要知道自己在做什么。填洼时要知道填的是伪洼地还是真实洼地调阈值时要知道阈值控制的是河网密度和单元尺度选坐标系时要知道它直接影响面积计算的正确性。我见过不少朋友把默认参数从头跑到尾做出来的结果支离破碎回头怪数据不好其实多半是没理解流向和汇流累积这两步到底在算什么。把这些原理吃透换任何软件、任何数据源你都能快速上手。