简介:面向洪水灾害评估、城市防洪规划及地理信息系统应用人员,这份研究报告系统梳理了基于ArcGIS的洪水淹没分析与三维模拟方法。内容先对比基于水位与基于水量两种分析思路,选用无源淹没模型,并说明其适用于洪水源易确定、地势相对平坦的区域。在ArcGIS实现层面,文档详细介绍了从TIN数据预处理生成数字高程模型,到利用Spatial Analyst中的栅格计算器提取低于设定水位的淹没区,再到统计淹没面积并建立洪水水位与淹没面积关系公式的完整流程;结合某建成区案例,展示了按五十年、二百年一遇设防水位划分区域天然防洪能力的做法,使分析结果更贴合实际应用。随后演示了使用ArcScene对水位渐进抬升过程进行三维模拟,直观呈现不同水位下的淹没范围变化。资源为单个doc格式电子文档,压缩包仅11KB,内容紧凑、可操作性强,适合GIS初、中级学习者对照练习。目前已有330人浏览学习,有助于快速掌握洪水淹没分析的核心思路与ArcGIS关键工具操作。
1. 洪水淹没分析不只是画一条等水位线
遇到过不少用ArcGIS做洪水淹没分析的人,第一反应是在DEM上按高度值做符号化分级,认为设到某个高程以上就是淹没区。这个做法应付汇报图勉强够用,但一旦要算淹没面积、与土地利用叠加统计、导出风险图,误差会大得没法交代。原因在于符号化只是把连续地形切成了几段颜色,既没处理洼地、也没有考虑水体在平面上的连通关系,更谈不上有多少水量、从哪里来、往哪里去。基于ArcGIS的洪水淹没分析与三维模拟,核心是把“地形淹没”这件事拆成两步:先在二维栅格上用确定性的方法算出淹没范围,再把范围叠回DEM生成三维体块。二维部分常见有两条路线——无源淹没按水位抬升,有源淹没按水流路径生长;三维部分则用TIN或栅格拉伸来表达水位面与地形的交切关系。这篇文章就按“数据准备—无源淹没—有源淹没—三维模拟—验证与排错”的顺序,把每一步的原理和ArcGIS里的具体操作讲透,适合正在做防洪规划、灾损评估或应急预案的技术人员参考。
2. 数据准备与无源淹没分析的基本范式
2.1 DEM 预处理必须先于一切分析
ArcGIS里的任何淹没分析都以数字高程模型DEM为底座,栅格质量直接决定后续计算的可信度。我一般要求DEM至少满足三个条件:一是投影坐标系,单位最好是米,避免用经纬度直接算面积导致结果不可用;二是浮点型存储,整数型高程在处理缓坡时会丢失精度;三是已完成填洼或至少知道洼地的数量和分布。具体做法是在ArcToolbox里走Spatial Analyst Tools → Hydrology → Fill,填洼阈值默认是无穷大,即把所有汇水洼地都抬平。实际项目里不建议一上来就用默认值,而要看DEM的噪声水平。如果原始DEM来自5米或10米格网,微小的人工沟渠、路基形成的假洼地会被无差别填平,导致淹没范围外扩。这时可以先做Focal Statistics取3×3邻域的中值,把孤立噪点压掉,再用Fill处理真正的凹陷。
import arcpy from arcpy.sa import * arcpy.CheckOutExtension("Spatial") dem = r"D:\flood\dem_5m.tif" dem_filled = Fill(dem) dem_filled.save(r"D:\flood\dem_filled.tif")代码里Fill函数返回的是填洼后的栅格,参数ZLimit可选,表示填洼的最大深度限制,设为None表示全部填平。实际使用时我常用ZLimit来控制填洼深度,比如设定为2米,只填掉2米以内的洼地,保留真实的地形凹陷。这样后续河道和堰塞体周围的淹没分析更接近自然状态。填洼后的栅格需要检查属性表中的最小值和最大值,确认没有负值异常——负高程在沿海地区合法,在内陆山区多半是原始数据的问题。
2.2 无源淹没:水平面法下的面积与边界提取
无源淹没指假设水位瞬间抬高到某个高度H,凡高程低于H的像元都被淹没,不关心水源位置和连通路径。这在蓄滞洪区启用、水库溃坝的极端情景估算里很常用,优点是计算简单、物理意义明确。实现上就是一条栅格计算器表达式:Con("dem_filled" <= H, 1, 0),得到0和1构成的二值栅格,1代表淹没。为了后续统计方便,我通常把Con的结果再转成整型,并做一次RegionGroup把连通的淹没斑块赋予唯一编号,这样每个斑块的面积可以用ZonalGeometry或RasterToPolygon转面后按属性表计算。
in_dem = r"D:\flood\dem_filled.tif" water_level = 35.5 inundated = Con(in_dem <= water_level, 1, 0) inundated.save(r"D:\flood\inundated_355.tif") polygons = RasterToPolygon(inundated, simplify=True, raster_field="VALUE") polygons.save(r"D:\flood\inundated_355.shp")这段代码里Con函数走的是Spatial Analyst的条件判断,in_dem <= 35.5生成布尔栅格,满足条件赋1,否则0。RasterToPolygon默认会把值为1的区域转成面,simplify=True表示对边界做抽稀,减少面要素的顶点数量。这里要特别说明:无源淹没并没有区分“连不连通”,如果一个低洼山谷和一个盆地的高程都低于水位,它们都会显示为淹没区;如果两个区域被山脊隔开,也想合并成同一片洪水体,那就得换有源淹没的思路。另外,水位H的取值通常来自水文频率计算的洪峰水位或某重现期下的控制断面水位,没做水动力模拟的项目里,H一般用设计洪水位代替。
2.3 无源淹没结果的两个常见误读与修正
第一个误读是把淹没栅格的像元数量直接当面积。栅格分辨率如果是10米,一个像元代表100平方米,统计时用GetRasterProperties取COUNT再乘像元面积才对;如果转了矢量面,属性表里的Shape_Area字段会直接给出平方米数,但转面时空白区也就是NoData区域会被跳过,边界处的面积略小于栅格真实面积,这在极高精度要求下需要注意。第二个误读是忽略高程基准面不一致。同一个研究区里,DEM的高程基准若是1985国家高程基准,而水位数据来自吴淞基准,两者相差约1.8米左右的常数偏移,必须先做统一,否则1米级别的淹没范围判断会失真。
修正做法:在转出淹没面后,用Tabulate Intersection叠加土地利用或人口分布数据,把淹没范围落到承灾体上。
landslide = r"D:\flood\landuse_2020.shp" inundated_poly = r"D:\flood\inundated_355.shp" stats = TabulateIntersection(inundated_poly, "OBJECTID", landslide, "LANDUSE", r"D:\flood\inund_landuse.dbf")TabulateIntersection参数第一个是区域要素,第二个是区域标识字段,第三个是被统计图层,第四个是统计字段,输出是DBF表。表里每一行对应区域要素的每个被统计类别,AREA字段给出相交面积。这样可以量化不同地类里被淹的面积,为灾损粗估提供直接输入。
3. 有源淹没分析的实现路径与关键参数
3.1 有源淹没与无源淹没的本质差异
有源淹没考虑的是洪水从某个入口持续流入,像元被淹不仅因为高程低,还因为有一条从源点到达它的路径。这个路径由地形坡度决定,水往低处流,但不能爬过高山,所以有源淹没天然排除了被山脊隔开的低洼盆地。工程上常见做法以种子蔓延为核心:从指定的进水口像元出发,向四周八邻域扩展,当相邻像元高程低于“起始水位加上沿程水头损失”时纳入淹没区,否则停止。严格来说这属于简化后的静水淹没模型,没有解圣维南方程,但它比无源淹没多了一个“水流可达性”约束,用来评估溃口周边范围的快速淹没非常合适。
ArcGIS里没有现成的“有源淹没”按钮,一般两种实现思路。一种是用Cost Distance或Path Distance加Con组合出累积代价面,把“到达某个像元需要翻越的最大地形阻力”算出来,再与水位比较;另一种是直接用Python配合Raster Calculator做区域生长,或者借助Watershed工具给定出水口反推汇流区。前一种思路稳定,适合批量水位方案比较;后一种思路更接近物理直觉,但需要写循环。多数情况下我推荐前者,下面给出具体参数和表达式。
3.2 用 Cost Distance 实现最小阻力可达性分析
先基于填洼后的DEM生成每个像元到源点的最小成本距离,这里的成本不是路程,而是“途经像元时相对源点需要克服的高程差”。具体做法是用Path Distance工具,把Cost raster设置为全1栅格,把Surface raster设置为原始DEM,Vertical factor选择Binary,这样高程每升高一点,成本就大幅增加;更简单的办法是先用Minus("dem_filled", source_level)生成相对高差栅格,再把正的高差部分作为成本。我常写这样一个步骤:
先创建源点栅格,进水口位置赋0,其他像元为NoData。然后用CostDistance得到累计最小成本表面。这个成本表面的含义是:从源点到达每个像元,沿途需要翻越的高程代价之和。注意这里不是简单的高程差值,因为路径可以选择从山脊两侧绕过去。
source_point = r"D:\flood\source.shp" # 点要素表示进水口 source_raster = RasterToPoint(source_point) # 需配合提取赋值 # 更可靠的做法是先转成栅格 arcpy.env.extent = dem_filled.extent arcpy.env.snapRaster = dem_filled source_r = arcpy.conversion.PointToRaster(source_point, "OBJECTID", r"D:\flood\source_raster", "MOST_FREQUENT", 10) cost_surface = CostDistance(source_r, 1, dem_filled, "Binary")CostDistance的第一个参数是源栅格,第二个是权重栅格即成本栅格,第三个是表面栅格用于考虑实际地表距离,第四个Vertical factor用Binary表,含义是只要上坡,成本即乘以一个极大系数,而下坡成本为1。这样生成的成本面数值越大,说明从源点到达这里越需要爬升。然后与水位做比较:假设源点的起始水位为H0,沿程水面坡降近似忽略,那么凡cost_surface小于某个阈值的区域就是可淹没区。阈值不是水位值本身,而是“允许的最大爬升高度”,在平原区往往取0.5-1米,代表洪水漫过微小田埂的能力;山区取2米以上,代表快速上涨时水头对局部鞍部的漫越。
{{INSERT}}### 3.3 区域生长法的ArcGIS Python实现
如果研究区小、源点数量有限,更精细的是区域生长。用Python循环处理,每次取当前边界像元的八邻域,判断相邻像元高程与当前水位的差值来决定是否纳入。这里要注意,区域生长如果写成逐像元循环,10米分辨率、10000×10000的栅格会非常慢;所以我一般把生长过程向量化:维护一个候选队列,优先处理最低高程的像元,这样避免了大面积重复扫描。
import numpy as np from scipy import ndimage dem_arr = arcpy.RasterToNumPyArray(dem_filled, nodata_to_value=-9999) inund_arr = np.zeros_like(dem_arr, dtype=np.int8) queue = [(r0, c0)] # 源点行列号 water_level = 40.8 while queue: r, c = queue.pop(0) for dr in [-1, 0, 1]: for dc in [-1, 0, 1]: nr, nc = r + dr, c + dc if 0 <= nr < dem_arr.shape[0] and 0 <= nc < dem_arr.shape[1]: if inund_arr[nr, nc] == 0 and dem_arr[nr, nc] <= water_level: inund_arr[nr, nc] = 1 queue.append((nr, nc))这段逻辑是广度优先搜索,每纳入一个新像元就把它加入队列,继续向外探。循环条件是水位恒定,如果模拟随时间上涨的水位,就要把water_level改成随迭代次数增加的变量,并且每次增长时重新扫描边界像元。不过要注意,Python队列处理大范围区域会慢,实际工作流建议先用CostDistance粗筛,再用区域生长做局部细化。
3.4 参数敏感性:水位、源点位置和DEM分辨率的影响
有源淹没对源点位置非常敏感。同一个DEM,源点设在下游河道和设在溃口处,结果差异极大;源点所在的像元如果是洼地底部,洪水会先填满洼地再外溢,生长范围受限;源点如果选在山脊上,则几乎不生长。我一般在源点生成前用Snap Pour Point工具把点对齐到流向累积量最大的像元,保证起始水位出现在河流的谷底线。水位参数的敏感性在不同坡度下差别很大——陡峭山区水位抬高5米,水平扩展范围可能只有几十米;平原地区水位抬高0.5米,淹没范围就成片扩张。做方案比选时,建议做成参数栅格图,水位H按0.5米间隔生成10组淹没面,叠加成一张风险图,比单一一套结果更有说服力。
4. ArcGIS 三维模拟:从淹没面到体块表达
4.1 三维场景的两种基础表达:栅格拉伸与 TIN
二维淹没分析得到的是一个面范围,三维模拟要表达的是“水位面以下、地形表面以上的空间体”。ArcGIS里常用ArcScene或ArcGIS Pro的Local Scene实现。第一步是构建地形表面,小范围用原始DEM直接作为高程源,大范围先转TIN再做Terrain优化。第二步是把淹没区做垂直拉伸,拉伸方式有两种:基于图层属性拉伸,和基于栅格值拉伸。如果淹没范围是从Con表达式里来的二值栅格,给0值设透明、1值设置基础高度为水位值,符号系统里选“按属性拉伸”,拉伸字段设为水位常数;更好的做法是先用RasterCalculator生成一个常量栅格:高程值等于水位的平面,再对这个平面做ExtrudeBetween,把地形面和平面之间的空间生成实体。
我用Python走一遍完整流程,核心是把DEM和水位面都转成TIN或带Z值的多面体。这里可以借助GP工具RasterSurface和ExtrudeBetween,但ArcGIS原生没有直接生成水体的工具,常用做法是把淹没区的矢量面要素增强为三维要素。
arcpy.env.workspace = r"D:\flood\scene.gdb" dem_tin = arcpy.ddd.RasterTin(r"D:\flood\dem_filled.tif", "TIN_DEM") water_surface = arcpy.ddd.CreateConstantRaster(water_level, "FLOAT", 10, extent) # 把水位面转成多面体 water_tin = arcpy.ddd.RasterTin(water_surface, "TIN_WATER") # 用ExtrudeBetween生成水位面与地形之间的体块 arcpy.ddd.ExtrudeBetween(dem_tin, water_tin, r"D:\flood\flood_body", "BETWEEN")RasterTin把栅格转成不规则三角网,ExtrudeBetween接受两个TIN,输出两者之间的闭合体块。BETWEEN模式生成地形与水位面之间的体积,ABOVE和BELOW则分别只保留某一侧的体块。生成后的多面体可以直接在ArcScene中设置透明度,做成半透明水体叠在影像上。不过ExtrudeBetween在数据量大时容易产生产状物,因为两个TIN的三角形边长不一致,交切处会出现细长三角形。常见解法是在转TIN时把z_tolerance设为0.5到1米,减少三角形数量,交切更干净。
4.2 场景符号化与飞行漫游导出
三维场景的视觉效果取决于两个设置:光照和水体透明度。ArcScene里水体图层的符号系统选择Simple Fill,把颜色改为浅蓝色,透明度调到40%左右;地形表面用带山体阴影的立体符号,叠加影像底图。淹没体块的底面和顶面如果是分开的,需要做一次Union确保体块闭合,否则在场景里旋转视角时会发现内部是空的,不利于评审汇报。导出动画用ArcGIS Pro的View → Animation功能,设定相机的路径和朝向,顺时针绕研究区一圈,帧数设为300帧左右,导出MP4。
4.3 体积计算:淹没水量和灾损关联
三维体块除了看,还能量算体积。ArcGIS的Surface Volume工具可以计算水体体积和表面积。输入是表示水位的栅格表面,基准面设为ABOVE,结果表给出参考高度以下的体积。这个体积就是淹没水量,可以与水文站的洪量数据对照,反过来验证淹没范围的合理性。也可以用PolygonVolume直接对三维体块要素计算体积。体积计算结果的单位取决于DEM的投影单位,如果DEM是米,体积就是立方米;如果是英尺,就得乘0.028317转换。
vol_table = r"D:\flood\volume.dbf" arcpy.ddd.SurfaceVolume(r"D:\flood\dem_filled.tif", "ABOVE", vol_table, 10, 35.5)SurfaceVolume的四个参数依次是输入表面、方向(ABOVE表示计算表面以上到参考平面的体积)、输出表、参考高程间隔和参考基准面高程。参考高程35.5米即水位,输出表里的Volume字段就是淹没水量。与前面无源淹没水面法对比,这个体积值更准确,因为它在每个像元上动态计算面积和高差,而不只是平面面积乘平均水深。
5. 结果验证与参数调试的接地气技巧
5.1 与历史洪水痕迹对照的三步验证法
计算完的淹没范围不能直接采信,要和历史数据对照。常见做法是找水文站记录的“历史最高水位”的痕迹线,比如墙上的水痕、桥墩上的标尺记录,在ArcGIS里把这些痕迹点或线叠加到淹没范围上,看吻合度。验证分三步:第一步,把淹没范围的边界提取出来,转成线要素,FeatureToLine,叠加痕迹线看偏移距离;第二步,统计痕迹点中有多少落在淹没范围内,用Select By Location的INTERSECT,淹中率低于70%说明参数偏保守;第三步,反查错误:把漏淹的点导出来,看它们的流向和DEM坡向,通常原因是水位偏低或者源点位置不对。
5.2 水位步长与栅格分辨率的最佳配合
做多水位方案时,水位步长和DEM分辨率之间有个经验配比。DEM水平分辨率为10米时,水位步长选0.5米产生的面积增量往往在平原区突变明显,而山区变化平缓;若水位步长小于DEM垂直精度,比如DEM垂直误差有0.3米,那0.1米步长就是虚假精度,算出来的面积差异全是噪声。我通常先查DEM属性表里的标准差,取标准差的一半作为最小水位步长。研究区是陡峭山区,DEM标准差可能有30米,那水位步长取1米甚至2米就够了;平原地区标准差只有2米,水位步长取0.2米才能反映微地形的淹没差异。分辨率方面,如果原始DEM是30米,不建议直接重采样成10米来追求精度,因为插值出来的假细节会误导淹没边界;正确做法是保持30米分析,但在出图前做一步Aggregate的平滑,只改善视觉效果,不改变计算值。
5.3 三条常见的失败路径和对应的排除方法
图表:常见问题与检查顺序
| 现象 | 直接原因排查 | 处理手段 |
|---|---|---|
| 淹没范围完全空白 | 坐标系不一致或源点落在NoData区 | 用Project Raster统一到DEM坐标系,IsNull检查源点像元 |
| 无源淹没面积远超预期 | DEM没有填洼,内部盆地全部积水 | 做Fill后重跑Con,并对比填洼前后DEM的差异 |
| 三维体块有破洞 | 水位面与地形面三角形尺寸差异过大 | 调整z_tolerance重建TIN,或把水位面重采样为DEM两倍像元大小 |
三个问题都出在“数据级别不匹配”。Coordinated Universal Time的错位、NoData填充值不同、像元对齐不一致,是最容易忽略的三个点。展开说,Con判断时如果水位值落入NoData区,结果是不参与计算的;而在ArcScene里显示时,NoData默认是透明,看起来像被掏空。最实用的检查方式是打开栅格的属性表,看COUNT和唯一值分布,如果唯一值只有0和1,且COUNT远小于整个范围的像元数,说明有NoData混入;用IsNull做掩膜,把NoData填成0或极大值,再做分析。
5.4 一套可复用的批处理模板
最后给出一套可以改参数直接跑的批处理骨架,把DEM路径、水位列表、输出目录三者分离,适合做多情景淹没分析。水位列表可以是多个断面水位,也可以按时间增长的演进序列。这个模板我在多个防洪评价项目里用过,替换数据路径后基本不用改逻辑。
import arcpy from arcpy.sa import * arcpy.CheckOutExtension("Spatial") dem = r"D:\flood\dem_filled.tif" out_dir = r"D:\flood\scenarios" levels = [32.5, 33.0, 33.5, 34.0, 34.5, 35.0] for h in levels: out_name = f"inund_{str(h).replace('.', '_')}" inund = Con(dem <= h, 1, 0) inund.save(fr"{out_dir}\{out_name}.tif") poly = RasterToPolygon(inund, simplify=True) poly.save(fr"{out_dir}\{out_name}.shp") print(f"水位{h}完成,范围面积={poly.getArea()}")打印的面积实际需要从结果中读取,这里展示的是循环写法;多情景的关键是水位列表的生成方式:可以用range(320, 361, 5)除以10生成等差序列,也可以直接从CSV文件读取水文站的预报水位。批量输出的栅格可以用MosaicToNewRaster合并成多层栅格,提供后续制图使用。三维模拟方面,把每个水位对应的淹没范围批量拉伸并存入同一个场景文件,逐层显示,就可以制作出随时间上涨的动画效果,比单张静态图清晰得多。
本文还有配套的精品资源,点击获取