做过地形分析的人都知道山脊线和山谷线是地貌形态的两个骨架一个管分水、一个管汇水。以前我在项目里拿到一块DEM第一反应就是开ArcGIS一顿操作提取坡度、坡向、曲率最后却卡在山脊山谷线上——这东西看着简单真想自动化提取出来传统等高线目视勾绘效率太低人工误差也大。后来我把水文分析和表面分析的思路打通配合ArcGIS里的所谓“图解建模工具”——也就是ModelBuilder——把整套流程固化成模型参数。输入一个DEM点一下运行山脊线、山谷线、河网全部出来还能批量处理整个研究区。这篇文章就把这套方法从原理到实操完整拆开适合正在做水文分析、地形地貌分析、或者需要批量出图写报告的GIS从业者。哪怕你只熟悉ArcMap的界面按照步骤走也能把模型搭起来以后换DEM只改输入文件就能复跑。1. 先把原理讲明白为什么水文分析能提取山脊线和山谷线很多初学者搞不明白一件事水文分析明明是研究水怎么流的怎么会跟提取山脊线、山谷线扯上关系这里面的核心逻辑其实就一句话——山脊线是水流的“发源地”山谷线是水流的“集结地”。放大地形来看雨水落在山坡上会沿着坡度最陡的方向往低处流。而山脊线附近的地形恰好是两侧坡面都往外倾斜水在山脊上汇不起来只会向两边分流。所以山脊线本质上就是分水线也就是水流路径的起点。反过来山谷线是两侧坡面都往里收的位置水从四面八方流过来聚在这里形成合水线也就是天然的沟谷、河道所在。既然水流的起点和路径分别对应了山脊和山谷那用水文分析来提取就很自然了先算每个像元的水流向哪里再统计每个像元有多少上游水流量汇入。流量大的地方就是山谷线流量为零或者极其小的地方就是山脊线的候选区域。这个思路比单纯靠曲率或者坡度阈值硬切割要稳定得多因为它自带上下游拓扑关系出来的线是连续的、有方向感的不会碎成一片一片的散点。不过有个问题需要解决水文分析天然提取的是汇水线也就是山谷线。山脊线是分水线和汇水线正好相反怎么用同一套算法提出来这就引出了整个方法里最巧妙的一步——反向DEM。1.1 山脊线与山谷线的本质先说山脊线。从地貌学角度山脊线是流域的边界一侧的水流向这个流域另一侧的水流向旁边的流域。在栅格DEM里山脊线通常表现为局部的高程脊线横切面向上凸出两侧坡度方向相反。山谷线则相反它是流域内部水流的汇集通道横切面向下凹陷两侧坡度方向指向谷底。用等高线看更直观山脊线上等高线向低处凸出山谷线上等高线向高处凸出。这个规律在地形图判读里叫“凸高为谷凸低为脊”。在栅格计算里我们通常用曲率来判断剖面曲率大于0是凸坡接近山脊小于0是凹坡接近山谷。但曲率计算对DEM噪声极其敏感哪怕是很小的数据抖动都会产生大量伪山脊、伪山谷。水文分析的优势在于它用了“汇水累积量”这个积分量天然做了平滑和滤波小的地形噪声对最终结果影响小得多。1.2 D8流向算法与“反向DEM”这个关键操作ArcGIS水文分析里最基础的流向算法是D8也就是单流向法。原理很简单对中心像元周围的8个邻居像元计算每个方向的高程差除以距离代价距离然后选择落差最大的方向作为水流方向。对角方向距离要乘以根号2这一点ArcGIS内部已经处理了不用我们操心。D8是所有后续分析的基础。有了流向栅格再统计每个像元被多少个上游像元指向就得到流量累积栅格。流量累积值越大代表这个位置越可能是河道、沟谷。这就是山谷线的雏形。但山脊线怎么用同样的逻辑提答案是把DEM整个翻转过来——用栅格计算器给DEM乘以-1原来的山脊变成“倒置的山谷”原来的山谷变成“倒置的山脊”。这样一来水流在反向DEM上会自然汇聚到原来的山脊位置。然后再走一遍水文分析的完整流程得到的汇水线就是原DEM的山脊线。这个思路是我第一次用时觉得最妙的地方等于把“分水岭”问题转换成了“汇水区”问题一劳永逸。1.3 整条技术路线的地图在动手之前先把整条流程链路列清楚后面每一步都是在链路上加节点。山谷线河网链路原始DEM → 填洼 → 流向 → 流量累积 → 栅格计算器阈值提取 → 河网栅格 → 矢量化 → 简化平滑。山脊线链路原始DEM → 栅格计算器乘以-1 → 反向填洼 → 流向 → 流量累积 → 栅格计算器阈值提取 → 山脊栅格 → 矢量化 → 简化平滑。对比一下就能发现山脊线就是比山谷线多了“乘-1”这一步。理解了这条链路ModelBuilder其实就是在画这个流程的图形化连接图。2. 数据准备与预处理DEM选料和填洼的门道这套方法的输入只有一个——DEM。但DEM的选型和预处理直接决定输出质量这一步偷懒后面全白干。我把这里容易踩的坑按顺序说一遍。2.1 DEM数据来源、投影与裁剪DEM的来源一般是公开的全球数据比如ASTER GDEM、SRTM或者是ALOS这类精度从30米、12.5米到1米都有。项目里如果只做宏观分析30米够用如果做小流域的精细沟谷提取最好用12.5米或者更高的分辨率不然细小的山谷线会被抹掉。数据拿到手第一件事是检查坐标系。水文分析要求DEM必须是投影坐标系不能是经纬度坐标系。原因很简单D8算法要计算像元之间的距离和对角线长度经纬度下每个像元的实际距离随纬度变化算出来的流向和流量在南北方向上会失真。所以要么在数据下载时选择投影坐标系产品要么用“投影栅格”工具把DEM转到合适的高斯投影或UTM投影下。还有一个细节如果研究区比较大建议先用矢量边界把DEM裁剪出来。裁剪有两个好处一是减少计算量流量累积这种工具在几千万像元的大范围上跑起来非常慢二是避免边界外的地形对边界附近流向造成干扰。裁剪推荐用“按掩膜提取”工具注意勾选“使用输入要素的几何特征”选项避免裁剪结果边缘出现锯齿状的NoData像元。2.2 填洼Fill为什么必做、Z limit参数怎么给填洼是水文分析里默认必须做的一步。DEM数据不管来自哪里原始值都带有大量局部洼地——就是比周围一圈都低的像元。这些洼地可能是真实的地形比如喀斯特漏斗、冰川湖但更多是数据插值产生的噪声点和错误。如果不填洼水流到洼地里就停了流向分析到此中断流量累积在这里断掉后面的河网提取就会碎得一塌糊涂。ArcGIS的填洼工具在Spatial Analyst工具箱下名字就叫“填洼”Fill有一个Z limit参数值得说。Z limit的含义是“允许被填的最大深度差”默认是无穷大也就是所有洼地全部填平。这个默认值在大多数地形区域没问题但如果你研究区有真实存在的洼地地形比如岩溶区的天坑、人工水库建议把Z limit设成一个合适的高程值比如3到5米。这样只填掉小于这个深度的小噪声洼地保留真实的大洼地。我实际用下来常规山地地形直接用默认值不会出错平原微丘陵区域反而要更谨慎因为那里的DEM噪声引起的伪装洼地太多填完之后地形被抹成一块平板流向结果会变得很不稳定。2.3 浮点型DEM的两个隐藏坑填洼和流向工具对输入的DEM数据类型没那么讲究但有两个隐藏坑还是需要提醒。第一个坑是DEM为浮点型时个别版本的ArcGIS填洼工具会出现结果没变化的现象。严格说这不是工具bug而是浮点精度导致的——浮点型DEM里几乎每个像元都有微小的高程差洼地的定义在这种精度下形同虚设。处理办法是先使用“整数”工具或者栅格计算器里的“Int()”函数把DEM转成整型再继续。整型DEM的流向结果也更干净下游矢量化时线段更规则。第二个坑是DEM存在NoData空洞。有的DEM在湖泊、云覆盖区域是空值这些NoData区域在流向计算里不参与运算导致周边像元的流向全部指向NoData结果栅格会出现一个黑洞似的不连通的区域。处理办法是在填洼之前先用“焦点统计”或者“按区域填补”把NoData区域插值填上。项目里如果空洞面积不大也可以用“栅格转点反距离权重插值”的方式补洞这一步做完地形平滑度反而更好。3. 正式提取山谷线与河网手把手流程预处理做完就可以正式进入水文分析链路了。我用ArcGIS Pro 3.x的界面来讲ArcMap 10.x的操作基本一致只是工具栏位置稍有不同。按顺序一共六步每步都是独立工具可以单跑也可以拼进模型。3.1 流向计算和流量累积第一步是“填洼”Fill输入预处理后的DEMZ limit按前面说的设置。输出命名为fill_dem。第二步是“流向”Flow Direction输入fill_dem输出flow_dir。这个工具默认使用D8算法我们保持默认即可。得到一个整型栅格值从1到255分别代表东、东南、南等8个方向。第三步是“流量”Flow Accumulation输入flow_dir输出flow_acc。这一步统计的是每个像元的上游汇水像元数量也可以理解成汇水面积。如果输入的是整型DEM这里的输出值会大得吓人——百万、千万级别别慌这是正常的。做完这三步打开flow_acc的属性表看最大最小值。最小值通常为0最大值取决于研究区大小和DEM分辨率。比如一个100平方公里、30米分辨率的区域最大汇水单元格数大约在11万左右。这个最大值在后面定阈值时用于比例估算先记下来。3.2 栅格计算器提取水系第四步是整个流程里最有技术含量的一步——从连续变化的流量累积栅格里提取出真正代表河网的那些像元。说白了就是定一个阈值流量大于等于阈值的像元算河道小于阈值的算坡面。实际操作在栅格计算器里写表达式SetNull(flow_acc 1000, 1)这个表达式的意思是flow_acc中小于1000的像元设为NoData其余的赋值为1。得到的栅格就是河网二值图1代表河道NoData代表非河道。阈值怎么定我个人的经验是“先小后大、逐步目视校准”。第一次先取一个较小的值比如最大流量累积值的0.1%跑出河网后叠加在卫星影像或山体阴影上看河道是否完整连接、有没有过多分支。如果细枝末节太多调大阈值如果主河道都断了调小阈值。一般取最大值的0.1%到1%之间个别平坦区域可能要放宽到5%以上。这个阈值不建议一次定死做成模型参数反而更方便。3.3 矢量化与平滑处理第五步是把河网栅格转成矢量线。工具是“栅格河网矢量化”Stream to Feature输入河网栅格和流向栅格flow_dir输出一个线要素类。这个工具会沿流向方向把栅格像元连成矢量线保证河网线的方向是从上游到下游这个方向信息后面做水文拓扑分析还能用到。第六步是简化和平滑。直接转出来的矢量线是锯齿状的因为栅格像元是方形的。在制图出图前建议先用“平滑线”Smooth Line工具做一次PAEK平滑容差给两个像元大小左右例如30米DEM就设60米。注意不要一轮做完就直接出图要看平滑后的线有没有穿过山脊或者偏离实际沟谷如果偏离明显说明原始河网定位就有问题需要回溯检查阈值。到这里山谷线河网提取就完成了。我把参数整理成一个速查表环节工具关键参数输出填洼FillZ limit可选默认全填fill_dem流向Flow DirectionD8算法flow_dir流量Flow Accumulation输入flow_dirflow_acc阈值提取栅格计算器SetNull(flow_acc 阈值, 1)stream_grid矢量化Stream to Feature输入stream_grid flow_dirstream_line平滑Smooth LinePAEK容差约2倍像元尺寸stream_smooth4. 提取山脊线的完整操作把DEM翻过来再做一遍山脊线提取的核心思想前文说过——把DEM反向。但实际操作中反向之后的处理和普通河网提取有一个重要区别这个区别是我踩坑多次才总结出来的先卖个关子。4.1 反向DEM制作第一步用栅格计算器把DEM乘以-1公式是Float(dem) * -1这里我习惯加Float转换防止部分版本直接把int乘-1输出成整型四舍五入丢精度。得到reverse_dem。第二步对reverse_dem做填洼。这是关键——很多人直接忽略了这一步想着正向DEM填过了反向就不用填了。大错特错。反向之后原来DEM上的山脊变成了“谷底”而真实地形中的山谷在反向DEM里变成“脊线”但更重要的是原本DEM上的山顶平台、鞍部这些区域在反向后会形成大量虚假的洼地。如果不填流向会在这些位置打转出来的山脊线断断续续甚至会出现一圈一圈的闭合环。所以一定要对reverse_dem单独做一次填洼。填洼参数和正向一样Z limit可以适当调大一点因为在反向地形里真实山脊两侧的高差往往比较大小阈值填不干净。第三步剩下的步骤就和河网提取一模一样了流向、流量累积、阈值提取、矢量化。最后得到的线就是原始DEM上的山脊线。4.2 山脊线提取后的修饰用这个方法直接跑出来的山脊线和普通河网比起来通常更碎、分支更多。原因是山脊在DEM上常常比较“钝”是一个宽缓的脊面而不是一条明显的脊线D8算法会把脊面上所有近似分水的位置都算作起点结果就产生一堆平行的短线。针对这个问题我常用的处理方法是做“密度聚类”式的筛选。先把山脊线按长度分类保留长度大于某个阈值的线把短线过滤掉。这个阈值没法给固定值得看研究区实际地形一般取平均线长的三分之一到二分之一。在ArcGIS里操作是用“按属性选择”加Length字段筛选然后导出。也可以用“消除”或者“简化线”来合并短线段。另一个修饰技巧是结合坡度真正的山脊线通常伴随明显的坡度转折把坡度栅格和山脊矢量叠在一起凡是穿过高坡度区域的线段保留完全位于平坦区域的短线大概率是伪山脊。虽然没有一个一键工具但做一次目视筛选加属性筛选能很大程度提升成果质量。5. 用ModelBuilder把整套流程固化下来前面写的流程如果每次都手动点工具一个区域至少要操作十几步。要是手里有20块DEM要处理光点鼠标就能点到怀疑人生。这个场景就是ModelBuilder发挥价值的地方。5.1 建模思路哪些参数要暴露哪些要锁死ModelBuilder说白了就是一个可视化流程编排窗口把工具和数据用箭头连起来让ArcGIS按顺序自动执行。它最核心的优势有三个可重复、可批处理、可共享给别人用。建模时要先想清楚哪些参数要作为模型参数即每次运行都可以改哪些参数要锁死。我的建议是作为模型参数输入的DEM、提取阈值、输出山脊线/山谷线的保存位置。锁死参数填洼的Z limit、流向算法D8、平滑容差。这些一般不需要频繁调整锁死可以让模型更稳定。这样做的好处是别人拿到模型后只需要填四个参数就能运行不容易改乱内部设置。5.2 模型搭建的分步操作打开ArcGIS Pro在分析选项卡里点击ModelBuilder新建一个模型。在目录窗格里展开Spatial Analyst工具箱里的水文工具集把填洼、流向、流量、栅格河网矢量化拖进画布再把“栅格计算器”或“提取分析”里的相应工具拖进来。然后用箭头连接工具和数据的顺序具体连接方式可以直接从工具的右键菜单选择输入输出。以山谷线链路的模型为例流程是DEM → 填洼 → 流向 → 流量 → 栅格计算器阈值 → 栅格河网矢量化 → 平滑线。山脊线链路则在最前面插入一个栅格计算器做乘-1。这里出现了一个会卡死很多新手的问题如果模型中同时存在两条链路怎么让它们各自引用正确的输入输出解决办法是右键中间数据比如fill_dem把它的中间数据属性打开或者在工具连接上改成用“预条件”Precondition控制执行顺序。更粗暴的办法是建两个模型分开跑山谷线一个模型、山脊线一个模型虽然多一个文件但运行结果更清晰不易出错。5.3 用模型批量跑多个DEM把DEM这个变量设置为模型参数后运行模型时可以直接输入多个DEM吗答案是不行模型参数界面一次只能填一个栅格。要真正批处理有两个办法。办法一右键模型里的DEM变量打开“迭代栅格器”Iterate Rasters设置一个存放多个DEM的文件夹模型会自动遍历文件夹里的每个栅格依次执行。这个功能看着很酷但有一个坑如果中间数据没有前缀区分上一次循环产生的中间文件会被下一次覆盖所以最好在中间数据的输出路径里加上“%名称%”这样的动态变量。办法二生成脚本。右键模型选择“导出为Python脚本”ArcGIS会把整个模型转成arcpy脚本。然后自己在脚本外层套一个for循环遍历文件夹里的所有DEM每次更新输入路径。这个方法的优点是运行速度快而且能加入异常处理某一个DEM跑挂了不影响后面的。我自己现在都用脚本方式ModelBuilder更多是用来快速搭原型验证流程。6. 常见翻车现场与排查技巧无论按我前面的步骤走得多顺实际项目里总会冒出各种奇怪问题。这一节专门收录我这些年遇到的典型问题和排查思路。6.1 填洼前后无变化怎么排查症状填洼工具运行结束输入输出两个栅格的属性表一对比像元值统计几乎一模一样。排查顺序第一步检查输入DEM是不是浮点型前面说了浮点型DEM可能存在精度问题先转整型再试。第二步双击填洼工具看Z limit是不是设成了0Z limit为0等于不填。第三步放大查看局部地形看DEM是不是已经被填得很平了如果是说明工具本身运行正常只是之前的DEM在这个区域确实没有洼地这是正常情况不是问题。6.2 河网断开是阈值问题还是数据问题症状提取出来的河网主线断断续续中游出现较长的空白段。先在断口处叠加山体阴影查看DEM如果断口处地形明显是沟谷但河网没提出来说明这个位置的流量累积值没达到阈值多半是上游汇水面积太小直接调低阈值即可。但要注意调低阈值会带来大量支毛沟所以更好的办法是先检查填洼是否彻底有时候是洼地切断了上游汇水路径导致下游流量骤降。这时候把Z limit调大重新填洼效果比盲目调阈值更本质。6.3 提取结果锯齿严重怎么处理症状矢量化出来的山脊线和山谷线呈明显的“阶梯状”制图效果很差。这类问题多发生在低分辨率DEM上比如30米或者更大像元尺寸。处理办法是按顺序来做第一步先把DEM重采样到更高分辨率用“重采样”工具双线性插值比如重采样到10米第二步再走一遍水文流程。但重采样不增加真实地形信息只让线条更光滑本质是用算力换视觉效果。第二步用“简化线”工具算法选“消除伪节点”或“弯曲简化”容差根据线层比例尺来设通常1:5万制图用50米左右容差。如果还嫌不够可以用“平滑线”的PAEK算法二次处理。还有一个小技巧在矢量化之前先把河网栅格用“栅格清理”工具的加厚或细化功能处理一下栅格的骨骼化结果会比原始像元连线更规则。这套流程跑通之后整个山脊线山谷线提取就不再是每次项目里的“玄学步骤”而是一个可复用的标准操作。我后来在多个项目里都用同一套模型只换DEM输入和阈值参数十分钟内出完整套地形骨架线。让我比较意外的是把这个方法教给团队里刚接触GIS的实习生他们也能顺利跑通可见这套方法的稳定性还是不错的。最后说一个我自己的习惯建模时不要把所有中间数据都设成临时文件。我倾向于把填洼结果和流向结果保留下来因为后面调试阈值、排查问题都要反复用每次重跑整个模型的时间成本远大于那点磁盘空间。等整个项目结束再统一删掉中间数据只留最终的山脊线、山谷线和必要的中间成果这样既方便过程复盘又不至于让数据目录乱成一锅粥。
