干GIS的,基本都躲不开洪水淹没分析这道坎。不管是做防灾减灾、洪水风险评估,还是水库泄洪影响分析,只要手上有一份DEM数据,大概率会被问到两个问题:“洪水会淹到哪里?”、“淹了多少水?”第一个问题对应提取淹没区,第二个问题对应计算淹没区洪水量。用ArcGIS基于DEM做洪水淹没分析,本质上就是把这两个问题落到栅格计算上,中间涉及填洼、流向、流量、连通性判断等一系列操作。
这篇教程不是泛泛讲工具,我尽量把每一步背后的道理讲清楚,包括为什么必须填洼、为什么有源淹没比无源淹没更贴近实际、为什么算水量的时候不能用填洼后的DEM。内容适合正在学ArcGIS水文分析的学生,也适合实际工作中需要快速出淹没范围图、算淹没水量的规划、水文、应急相关从业者。软件操作我按ArcGIS 10.x的Spatial Analyst扩展模块来写,ArcGIS Pro的操作路径也大同小异,可以直接对应。
1. 洪水淹没分析的整体思路与工具选型
1.1 无源淹没与有源淹没:如果连这个都没分清,后面全白做
网上很多教程一上来就教你用栅格计算器写一句Con("DEM" <= 水位, 1),这叫无源淹没分析。它的逻辑很简单:把所有高程低于水位的像元全部当成淹没区。这种做法在什么场景下勉强能用?大范围、快速初判,比如整个平原区在某个水位线以下的土地面积有多大。但它的致命缺陷也很明显:它会淹没掉所有低于水位的独立洼地,哪怕这些洼地和河流之间隔着一道山脊,水根本流不过去,也被算成淹水。
有源淹没才是更贴近实际的思路。它要求淹没区必须从某个“源点”(比如河道出水口、决口位置、水文站断面)开始,通过连续的地形路径向外扩展,扩展过程中不能跨越高于水位的阻隔。简单说,水是从河里漫出来的,而不是从天上均匀落下来的。
这两种方法在ArcGIS里都能实现,但工具组合完全不同。无源淹没一条栅格计算器语句就能解决,有源淹没则需要借助填洼、流向、区域分组等工具做连通性分析。我这篇教程会把两条路都讲清楚,重点说透有源淹没的做法,因为实际项目中真正用得上的大都是有源淹没。
1.2 技术路线总览与工具选择
做一次完整的洪水淹没分析,常规流程分为四步:数据准备与预处理、河网提取与源点确定、淹没区提取、洪水量计算。
数据准备阶段需要一份质量可靠的DEM,然后依次做填洼、计算流向、计算流量。河网提取阶段需要根据流量阈值提取河道栅格,再矢量化成河网线,并在合适的位置定出淹没分析的源点。淹没区提取阶段是核心,我会采用“潜在淹没区+区域分组+源点连通”这一套逻辑,比单纯的阈值判断靠谱很多。洪水量计算阶段则是把水位面与原始DEM做差,得到水深栅格,再做体积统计。
工具方面,基础需求只需要ArcGIS的Spatial Analyst扩展模块,里面包含了填洼、流向、流量、区域分组、栅格计算器等所有关键工具。如果后续想用三维方式核验体积,可以再加一个3D Analyst模块。
整个流程中我特别想强调两个容易被忽略的细节,第一,流向和河网分析要基于填洼后的DEM,但洪水量计算必须基于原始DEM,否则填洼把洼地填平了,算出来的体积会明显偏小。第二,所有面积和体积计算都应该在投影坐标系下进行,直接拿经纬度的DEM算面积会出错,这点后面会专门展开。
2. 数据准备与预处理:从DEM下载到能算的水域
2.1 DEM数据源与坐标系校验:投影不对,后面全是灾难
做洪水淹没分析,DEM数据的分辨率和坐标系是两个决定性因素。常用的公开DEM包括30米分辨率的ASTER GDEM和SRTM,部分地区能拿到12.5米的ALOS数据,如果项目精度要求高,还可以用无人机航测或LiDAR点云生成的高分辨率DEM。分辨率越高,地形细节越丰富,淹没边界越准确,但数据量和计算时间也会成倍增加,实际项目里30米DEM做区域级淹没分析已经够用。
拿到DEM后的第一件事,不是急着填洼,而是检查坐标系。这里要特别强调“定义投影”和“投影”这两个操作的区别。定义投影只是给数据打上一个坐标系标签,并不改变数据的实际位置和形状;投影才是真正把数据从地理坐标系转换到投影坐标系。如果一份DEM没有正确的投影信息,或者坐标系标错了,流向计算、填洼、面积统计这些后续操作都会跟着出错,而且这种错很难直观发现,往往是最后算体积时才发现结果离谱。
我的习惯是先把DEM统一投影到等面积投影坐标系,比如UTM或阿尔伯斯等积圆锥投影。原因很简单,洪水淹没分析最终要算淹没面积和洪水量,这两个指标都强依赖面积精度,只有等面积投影才能保证栅格每个像元的实际面积基本一致。投影完成后,再用“按掩膜提取”把DEM裁剪到研究区范围,顺带把边缘的Nodata问题处理掉。Nodata区域在流向计算时会变成黑洞,直接影响河网提取结果。
2.2 填洼(Fill):为什么这个步骤死活不能省
DEM数据里几乎必然存在洼地,可能是原始地形的真实凹陷,也可能是数据采集和处理过程中产生的伪洼地。如果不填洼,水流在流向计算时会直接停在洼地里,形成断头路,后续的流量累加和河网提取都会断成一截一截,根本没法用。所以填洼是水文分析中最基础、也最不可跳过的一步。
ArcGIS里的填洼工具路径是Spatial Analyst Tools → Hydrology → Fill,操作上只需要输入DEM,工具会自动把所有洼地填充到与周围最低出水口平齐。这里有一个常被忽略的参数Z limit,它限制最大填充深度。如果不设这个参数,工具会把一些真实地貌上的深坑也一并填平,比如喀斯特地貌的漏斗、人工采石坑,填完之后地形变得过于平坦,与实际不符。我在实际项目中通常会给Z limit设一个经验值,例如10到20米,这样既能消除数据噪声造成的浅洼地,又不会抹掉真实的地形凹陷。
填洼完成后,建议用栅格计算器看一眼填洼前后的差值,"Fill_DEM" - "DEM",如果差值栅格里出现大面积的高值区,就要警惕是不是Z limit设置得太激进,把不该填的地形也填掉了。
2.3 流向与流量:连通性分析的底层依赖
填洼之后紧接着算流向。ArcGIS的Flow Direction工具默认采用D8单流向算法,也就是假设水流只从中心像元流向周围8个邻居中坡度最陡的那个方向。D8算法简单高效,非常适合洪水淹没分析这种需要明确“水往哪走”的场景。
流向栅格的每个像元值都是一个特定的方向编码,1表示东,4表示南,16表示西,以此类推。这个栅格本身并不直观,但后续所有水文分析几乎都要依赖它,包括流量累加、水流长度、河网提取。
流量累加工具Flow Accumulation是在流向栅格的基础上,对每个像元统计有多少上游像元的水汇入该像元。输出栅格里数值越大的像元,代表汇水能力越强,也就是河道的主干位置。默认情况下工具会使用权重栅格为1,也就是说每个像元贡献一份水量,这样得到的流量栅格值实际是上游汇水像元的数量,把它乘以每个像元的面积,就能得到理论汇水面积。
这里提醒一句,流量累加的输出值经常出现很大的数字,几万甚至几十万都很正常,不用觉得数据出错了。后面提取河网时,可以把这个值当作河网密度的调节参数来用。
3. 提取淹没区的完整实操:从河网到连通淹没区
3.1 提取河网并确定淹没源点
有了流量栅格,提取河网就很简单了。栅格计算器里写一句Con("FlowAcc" > 阈值, 1),把流量大于某个阈值的像元标记为河道,其余的设为NoData。这个阈值的选取没有固定标准,阈值越小,河网越密集;阈值越大,留下的只有主河道。对于30米分辨率的DEM,我一般先从流量像元数的1%开始试,比如总像元数是100万,阈值就取10000,然后根据生成结果和实际地形不断调整。
河道栅格可以用Stream to Feature工具转成矢量线,工具箱路径是Spatial Analyst Tools → Hydrology → Stream to Feature。输入河道栅格时建议先执行一次栅格清理,把一些小碎块去掉,避免矢量化之后出现大量细碎的短线,影响成图效果。
接下来是确定淹没源点。源点就是你希望模拟洪水从哪里开始漫溢的位置,可以是河道出水口、堤防决口点,也可以直接选在河道中段某个位置。如果手头有水文站的实测断面位置,那是最好的源点选择。如果没有,可以在河网矢量线上手动新建一个点要素,落在河道上即可。
这里推荐使用Snap Pour Point工具来做源点整理,它能把手动画的点自动捕捉到流量最大、汇水能力最强的河道像元上,避免因为点在河道边缘导致后续连通性分析落空。
3.2 水位怎么定:这一步决定了分析结果是否可信
水位是整个分析中输入最敏感、也最需要专业判断的参数。定水位的方法主要有三种。有实测数据的情况下,直接用水文站实测的最高洪水位或设计洪水位,这是最可靠的方案。没有实测数据时,可以用源点位置的地面高程加上预估淹水深度,比如源点高程是32米,预估淹水3米,水位就取35米。极端情况下如果连预估深度都没有,可以做多水位情景模拟,分别计算35米、36米、37米的水位下淹没范围和水量,形成不同重现期下的淹没等级图。
我在实操中经常发现一个误区,很多新手把水位取成某个固定海拔高度,却忽略了水面在河道上下游之间是有坡降的。对于小范围、短河段的淹没模拟,固定水位精度尚可接受;如果研究区跨度很大,沿河水位变化明显,就需要把水位面做成随距离变化的栅格,或者分段设置水位。这个进阶做法可以先记下来,我后面在洪水量计算部分还会提到。
3.3 核心实操:区域分组法提取有源淹没区
无源淹没只需要一句栅格计算器语句,这里不浪费篇幅。有源淹没的提取,我强烈推荐用“区域分组法”,它不需要写复杂的循环代码,整个逻辑清晰,结果也可靠。
第一步,生成潜在淹没区。假设水位为35米,在栅格计算器里执行Con("Fill_DEM" <= 35, 1),得到所有高程低于水位的像元。之所以用填洼后的DEM,是因为在有源淹没分析中,我们希望水流可以穿过这些原本是洼地的区域继续向下游流动,如果直接用原始DEM,洼地边缘的地形会让连通性关系变得特别复杂。
第二步,对潜在淹没区做区域分组。使用Region Group工具,路径是Spatial Analyst Tools → Generalization → Region Group。这里有一个关键设置,把Number of neighbors to use即邻域类型改成EIGHT,也就是八邻域连通。为什么不用默认的FOUR?因为水流漫延是可以斜着走的,八邻域连通更符合洪水横向漫延的真实情况。
第三步,找出源点所在的区域组号。这一步有两种做法,一种是用Extract Multi Values To Points工具直接把源点所在栅格的组号提取到点属性表里,另一种是用Zonal Statistics as Table统计每个组的分区属性,再通过源点位置去对应。前者更直接,我个人习惯用前者。
第四步,把源点所在的分组提取出来。假设提取到的组号是7,在栅格计算器里执行Con("RegionGroup" == 7, 1),得到的就是从源点出发能够连续连通到、且高程低于水位线的区域,这就是有源淹没区。
看明白这个方法的精妙之处了吗?它用一次区域分组,就把“高程低于水位”和“与源点连通”这两个条件同时满足了。那些虽然高程低于水位、但被山脊隔开的孤立洼地,因为不在源点所在的连通组里,会被自动排除掉,这正是有源淹没和无源淹没的差别所在。
如果不想用工具界面,这段流程用ArcPy脚本实现也很快,核心代码大概是这样的:
import arcpy from arcpy.sa import * arcpy.CheckOutExtension("Spatial") arcpy.env.workspace = r"D:\flood" arcpy.env.overwriteOutput = True # 输入数据 fill_dem = Raster("fill_dem") source_point = "source_point.shp" # 生成潜在淹没区 potential = Con(fill_dem <= 35, 1) potential.save("potential_inundation") # 区域分组 groups = RegionGroup(potential, "EIGHT") groups.save("region_group") # 提取源点所在组号 group_value = ExtractMultiValuesToPoints(source_point, [groups], "NONE") group_id = 0 with arcpy.da.SearchCursor(source_point, [groups.name]) as cursor: for row in cursor: group_id = row[0] # 提取连通淹没区 inundation = Con(groups == group_id, 1) inundation.save("inundation_area")3.4 淹没区栅格转面与后处理
淹没区栅格生成后,可以转成面要素方便制图和后续统计。用Raster to Polygon工具,输入淹没区栅格,字段选择VALUE,输出即得到面要素。需要注意的是,栅格转面后边界是锯齿状的,直接出图不够美观,可以用Cartography工具条里的Smooth Polygon或者ArcToolbox里的Smooth Polygon工具做一次平滑处理,平滑容差可以根据制图比例尺来定,一般10到30米比较合适。
另外还有一个常见问题,淹没区面和河网线、原始等高线叠加时,可能会存在因为栅格分辨率导致的位置偏移,这不是操作错误,是栅格离散化的正常现象。如果项目要求高精度边界,建议用更高分辨率的DEM重新跑一遍流程,而不是盲目后处理。
4. 淹没区洪水量计算与成果表达
4.1 洪水量计算的原理与单位陷阱
洪水量计算的本质是:把淹没区内每一个像元对应的水深加总,再乘以单个像元的面积。水深就是水位减去该像元的地面高程,换算成公式就是洪水量等于Σ(水位−高程) × 像元面积。
这里有一个极其容易踩的坑,就是我前面反复强调的,算水深和体积必须用原始DEM,不能用填洼后的DEM。填洼把洼地填平了,相当于抹掉了一部分库容,直接用填洼DEM算出来的体积会小于真实洪水量。很多新手在这里栽跟头,算出来的水量怎么都对不上实测值,排查到最后才发现是用了fill_dem。
单位换算也特别重要。如果投影坐标系是米制,像元大小是30米×30米,那么单个像元面积就是900平方米。高程单位一般也是米,所以水深乘以像元面积得到的体积单位是立方米。但有些DEM的高程单位是厘米,比如某些精细化测绘成果,这时候如果不统一单位,计算结果会差100倍。拿到任何DEM,第一件事就确认单位和坐标系,这个习惯能帮你避开很多低级错误。
另一个常见的单位陷阱是经纬度坐标。如果直接在WGS84经纬度坐标下做分析,像元大小是0.000277度这种,无法直接换算成米,算出来的面积和体积就毫无意义。所以在分析前一定要投影到米制等面积坐标系。
4.2 实际操作:水深栅格与体积统计
计算洪水量最直观的做法分三步。
第一步,生成水深栅格。在栅格计算器里执行"水位面" - "原始DEM",水位面可以直接填数字,比如35米,也可以是一个表示水面的栅格数据。得到的栅格里正值代表该像元会被水淹没,负值代表地面高于水位。如果研究区内水面有坡降,这里就体现出水位面栅格的优势了,我们可以先构造一个随距离变化的水位栅格,再做整个计算。
第二步,把水深栅格裁剪到淹没区。用Extract by Mask工具,输入水深栅格,用淹没区作为掩膜,得到只包含淹没区内水深像元的栅格。
第三步,统计水深总和。打开这个裁剪后水深栅格的属性表,直接看Sum字段,再把Sum乘以像元面积,就是总洪水量。如果属性表里的Sum不好找,也可以用Zonal Statistics as Table,分区栅格用淹没区,统计类型选SUM,输出表里的SUM字段就是水深总和,乘上像元面积就是洪水量。
举个例子,某研究区DEM是30米分辨率,投影坐标系为UTM,水位取35米,计算出来的水深栅格Sum值是1250000,那洪水量就是1250000 × 900 = 11.25亿立方米,同时淹没面积可以通过水深大于0的像元数来算,假设是50000个像元,那么淹没面积为50000 × 900 = 4500万平方米,也就是45平方公里。
如果安装了3D Analyst扩展,还可以用Surface Volume工具做交叉验证。输入原始DEM,选择参考平面ABOVE,平面高度填水位值,Z因子填1,工具会输出参考平面与DEM之间的容量体积。需要注意的是Surface Volume算的是整个DEM范围内低于水位面的空间体积,如果DEM范围大于淹没区,结果会比实际洪水量偏大,所以最好先把DEM裁剪到淹没区范围再用这个工具验证。
4.3 成果整理:从地图到统计表
淹没分析和洪水量算完之后,成果输出通常要包含三个部分,淹没范围图、洪水量统计表、淹没区属性清单。
淹没范围图的核心要素是淹没边界和背景地形。把淹没区面要素叠在山体阴影图上,透明度设置为40%到50%,用蓝色系填充,看起来既专业又直观。在水文分析里,山体阴影和淹没范围的叠加能让决策者一眼看出淹到了哪条等高线。
洪水量统计表建议至少包含淹没面积、平均水深、最大水深、总洪水量这四项指标。这些指标都可以从水深栅格的属性表或Zonal Statistics输出表中提取。平均水深等于水深Sum除以像元数,最大水深可以用水深栅格属性表里的Max值。
淹没区属性清单可以按行政村、地块或流域分区统计淹没面积和水量,这一步用Zonal Statistics as Table就能做到,分区栅格导入行政区划面转成的栅格,值栅格用淹没区或水深栅格,输出表格里每一行对应一个分区,面积和水量一目了然。这种分区统计结果在应急管理、保险定损、灾后评估里非常实用。
5. 高频报错与避坑指南
5.1 DEM预处理阶段的典型问题:投影、Nodata和填洼异常
很多人在DEM预处理阶段就卡住了。最常见的是填洼工具或者流向工具直接报错或者输出结果全黑。排查顺序一般是:先确认Spatial Analyst扩展模块是否启用,再确认输入DEM的坐标系是否有效,最后确认DEM范围里有没有大面积的Nodata区域。Nodata边缘在填洼时会被当作洼地边界参与计算,导致边缘出现异常抬升。
投影坐标系问题方面,如果同事发来的数据看起来没问题,但填洼结果特别奇怪,不妨检查一下是不是数据被错误地定义了投影。有些数据下载下来是经纬度坐标,但原始文件里没有写坐标系信息,ArcGIS默认把它当作未知坐标系,后续所有操作都建立在错误基础上。这时不要急着Project,先用Define Projection正确申明坐标系,再投影到目标坐标系。
使用Arc Hydro工具做DEM reconditioning报错也是高频问题,很多做淹没分析的人习惯先把DEM和河网套合处理再进入水文分析。如果这一步报错,常见原因有三个:一是DEM范围没有完全覆盖河网线,二是河网线是MultiPart要素且没有做拓扑检查,三是DEM像元和河网线的坐标系不一致。我的建议是,如果只是做淹没范围分析,不一定要用Arc Hydro的reconditioning,用前述的填洼+流向+流量流程就足够了,工具越少,出错的概率越低。
5.2 淹没分析和体积计算阶段的坑
提取淹没区阶段,我最常被问到的问题是,区域分组后源点所在组号提取不到,或者组号是NoData。出现这种情况,九成原因是源点位置的高程高于设定水位,也就是源点本身就不在潜在淹没区里。遇到这种问题别急着怀疑工具,先检查源点的Z值、水位设定和源点是否落在河道像元上。
还有一个隐蔽问题,Region Group默认的邻域类型是FOUR,如果忘记改成EIGHT,生成的连通组会比较破碎,淹没区边缘会往外扩散不足,导致面积偏小。这种错误从结果图上看往往不明显,但和真实淹没范围一对比就会露馅。
洪水量计算阶段,如果发现体积明显偏大或偏小,优先检查三件事。第一,是否误用了填洼后的DEM作为高程基准。第二,投影坐标系是不是米制,像元面积是否计算正确,如果研究区在高纬度用了Web Mercator,面积会被严重放大。第三,水深栅格的属性表里Sum字段是否被Nodata污染,如果淹没区掩膜没有覆盖完整,部分像元会变成NoData,统计结果就会偏小。用Extract by Mask时,可以顺手用Con("水深栅格" > 0, 水深栅格)把小于0的像元过滤掉,这样的统计结果更干净。
5.3 我个人的几个习惯,分享给大家
做这套分析做多了,我有几个固定习惯。第一,所有中间栅格数据都命名为fill_dem、flow_dir、flow_acc这种带下划线的英文名,绝对不用中文名和空格,ArcGIS对中文路径和特殊字符的兼容性时好时坏,没必要给自己挖坑。第二,每一步重要输出都备份一份原始版本,比如填洼前后的DEM、投影前后的DEM,后面排查问题时会省很多时间。第三,在正式分析前,先小范围跑一遍流程,确认每个工具的输入输出都正确,再全区域跑大范围,避免等了一个小时的运算最后发现参数设错了。
另外,如果项目允许,尽量用ArcGIS Pro来做,Pro在栅格计算、区域分组、水文分析这些工具的稳定性比10.x好不少,大范围数据运算也更流畅。而且Pro 3.0以上的水文工具集中还自带Flood工具,可以一步完成有源淹没模拟,思路和区域分组法一脉相承,适合用来交叉验证手动流程的结果。
最后分享一个我自己常用的验证方法:把算出来的淹没区边界和现场照片或者已知历史淹没范围叠在一起看看,哪怕只是目视对比,也能提前发现很多数据层面的问题。技术分析的最终目的不是为了输出一张漂亮的图,而是这个结果能不能经得起实际检验。