LISFLOOD-FP避坑指南：从DEM单位到结果解析的完整排查链路

坐标系里的“度”混进模型后水深直接冲到几百米。这不是LISFLOOD-FP模型本身离谱而是数据单位在背后捣鬼。很多人第一次跑这个模型栽得最惨的往往不是参数不会调而是DEM单位、降雨强度单位、参数文件里那些“看起来不用管”的字段以及results文件夹里躺着的二进制文件根本读不出来。这篇东西我想用自己踩过的坑把LISFLOOD-FP的数据单位、参数文件、结果解析一次讲清楚给正在被洪水模拟折磨的人一个能直接抄作业的避坑路径。LISFLOOD-FP作为经典二维洪水淹没模型在水利、城市内涝、溃坝分析领域用得极广但它的文档写得直白到近乎简陋很多细节都得靠试错。无论你是刚接触模型的新手还是已经从ASCII结果里画过几张图的老手只要碰过这几类问题这篇应该都能对你有用。1. 数据单位是第一个坑别让DEM和入流单位悄悄毁掉你的模拟我先说结论LISFLOOD-FP内部默认使用国际单位制长度是米时间是秒。但模型不会主动检查你的输入数据单位它只会按照数字去算。你给它的DEM如果是“度”它就当成“米”来算你给的降雨强度如果是mm/hr它也当成默认单位直接参与运算。结果不是模型崩了而是它非常礼貌地给你算出一个毫无物理意义的水深。1.1 DEM的水平和垂直单位必须都是米但现实中它往往是“度”我第一次用SRTM数据跑一个小流域时下载的是WGS84坐标系下的经纬度格网DEM范围大概是东经120.0到120.1度北纬30.0到30.05度cellsize写的是0.0009。当时心里想0.0009这个精度很高啊就直接丢进模型。跑完出来的水深最大值是423米我当时还以为是模型版本有问题后来才发现LISFLOOD-FP会拿这个cellsize当作网格长度去计算流速、时间步长和淹没面积。你想想0.0009度在赤道附近大概对应100米但在北纬30度经度方向的实际距离只有约96米而模型不会管这个它直接把0.0009当成0.0009米来计算。这就等于把真实地球上相距100米的两个栅格压缩成了1毫米不到水平网格尺度严重失真流速和水深自然全乱套。所以无论你的DEM来自SRTM、ASTER还是ALOS只要它是经纬度坐标第一步永远是投影到你所在区域的米制坐标系最常用的是UTM带或者你所在国家的平面坐标系统。投影之后务必检查.asc文件头部里的cellsize它应该是一个合理数值比如30、90、100而不是0.000几。还有个隐蔽问题如果DEM的垂直单位是厘米或英尺而水平单位是米也会让坡度和汇流关系出错。我见过有人从某省测绘局拿到等高线生成的DEM垂直单位标的是cm但文件名写的是DEM_Meter结果模拟出来的淹没范围比实际小一大圈。检查方法很简单打开DEM的元数据或者在一个小范围内计算最大最小高程差如果山地区域动辄几万米那基本是单位错了。1.2 降雨和入流边界单位是容易忽略的“双单位”陷阱DEM问题通常一次就能发现真正反复出错的是降雨强度和入流边界。LISFLOOD-FP的许多版本中降雨强度单位是m/s但气象部门给你的降雨资料几乎都是mm/hr。这两者差着数量级1 mm/hr等于2.77778e-7 m/s。如果你把50 mm/hr直接写成50模型会觉得天上每分钟掉下来50米高的水结果当然恐怖。如果是流域内均匀降雨通常用一个降雨强度值即可如果是时空变化的降雨要用rainfall文件。无论哪种我都建议你在参数文件旁边写一个单位换算备注把原始观测值、转换因子、最终值都列清楚。我以前吃过一次亏脚本里把转换因子写成2.78e-4多了三个数量级结果模型第一小时降雨就积了3米深的水排查了整整两天才发现是乘错了指数。入流边界同样有鬼。模型里河道入流或边界流量通常要求m³/s但你在实际工程中拿到的可能是水文站日均流量单位是m³/s没错可有人会拿到m³/day或者从流量过程线工具里导出的L/s。这些单位换算不复杂m³/day要除以86400才变成m³/sL/s要除以1000。问题在于当你面对一个几千个时段的边界文件时很容易只在文件开头改了一个单位后面时段全是原始值模型的入流过程整体偏大或偏小下渗、淹没范围全都会跟着错。1.3 我的单位自检五步法踩了几次坑之后我给自己定了一套单位自检流程每次模拟前都走一遍把所有外部输入文件列一张表包括DEM、降雨、入流、初始水深每项标注原始单位。统一转换成SI单位并把转换后的文件单独放在一个units_converted目录里绝不改动原始数据。在.par参数文件的注释部分如果版本支持写清每个输入文件对应的单位。写一个简单的Python脚本读入所有转换后文件的位置和全局统计值检查最小值、最大值是否在合理范围比如降雨强度在1e-8到1e-4 m/s之间入流在0到几百m³/s之间。跑一个理想化案例比如在10km长、1km宽的平底矩形河道里给固定入流看稳定后的水深是否和曼宁公式手算接近。如果连理想案例都对不上说明单位或公式理解还有问题。这套办法看似繁琐但能帮你省掉后面数不清的返工。模型跑一次可能几小时单位错了跑完才发现才是最痛的。2. 参数文件看似自由格式每个关键字都在左右结果LISFLOOD-FP的.par文件不是标准配置文件它更像一串“关键词: 值”的行。不同版本能接受的关键字不完全一样你多写了不认识的字段它会直接报错你少写了必需字段它就用默认值有些默认值大到足以毁掉模拟。2.1 一份最小可运行的.par长什么样下面是我经常用的一套模板尽量保持简单DEMfile: input/dem_meter.asc resroot: output/flood_res sim_time: 7200 initial_tstep: 2.0 resstep: 300 massint: 300 fpfric: 0.03每个字段的含义DEMfile输入DEM的ASCII栅格文件路径单位必须是米。resroot结果文件前缀模型会在这个前缀后面加上各种后缀所以它的路径就决定了results文件夹内容的位置。sim_time总模拟时间单位秒。7200就是两小时。initial_tstep初始时间步长单位秒。这个值如果太大模型可能一开始就因CFL条件不稳定而崩溃。resstep结果输出时间间隔单位秒。这里设置300意味着每5分钟输出一个结果帧。massint质量守恒统计的输出间隔单位秒。模型会每隔这么多秒记录一次总水量判断模拟是否守恒。fpfric洪泛区曼宁糙率。对于单一糙率场景这个值就够了如果要空间变化有另外的栅格文件字段。有些版本还会支持end_time、initial_water_depth、boundary_file等字段具体以你编译的版本手册为准。我的经验是拿到一个新版本先跑通这个最小案例再逐步加复杂功能别直接上大流域否则一旦报错很难判断是参数文件格式还是物理条件的问题。2.2 参数字段与单位的对应关系我用一张表整理一下常见字段和单位这张表我贴在自己工位对面字段名含义常用单位备注DEMfile高程数据文件路径米必须为投影坐标系sim_time模拟总时长秒若写成小时会缩小3600倍initial_tstep初始时间步长秒受CFL约束宁小勿大resstep结果输出步长秒太大会漏掉洪水过程细节massint质量统计间隔秒守恒诊断用fpfric洪泛区曼宁n-无量纲常取0.02-0.15fpfricfile空间曼宁n栅格-栅格数值对应n值rainfallfile降雨过程文件m/s不是mm/hrqfile边界流量文件m³/s可以是过程线很多人会在sim_time那里直接写7200觉得单位无所谓结果预期是两小时模型却只跑了两秒。也有人在resstep里写了300秒但sim_time写2小时导致只有一个结果帧输出文件小得可怜。这些数字看起来不起眼却直接决定结果文件有没有意义。2.3 最容易踩的三个参数坑第一个坑是路径分隔符。Windows环境下如果你写DEMfile: D:\models\dem.asc\d会被解析成转义字符模型十有八九找不到文件。我建议在参数文件里统一用正斜杠D:/models/dem.asc或者使用Linux环境跑模型否则你会在文件读取错误上浪费大量时间。第二个坑是initial_tstep的设置。LISFLOOD-FP会自动调整时间步长但初始步长如果和网格尺度、水深不匹配很容易让模型在前几步就发散。比如DEM格子是100米初始水深也是100米量级步长给了10秒弗劳德数可能直接破表结果文件里出现Inf或NaN。我的经验是先给一个保守的初始步长比如0.5 * cellsize / sqrt(g * max_depth)估算一下然后取一半等模型跑稳再让它自适应。第三个坑是resroot不能带上“results”这种目录前缀后忘记创建目录。模型通常不会自动创建文件夹如果你指定的路径里目录不存在它会直接报错或者静默写不进去。所以每次新建案例我会先手动把output/和results/目录建好并且在参数文件里确保resroot只是前缀而不是一个完整文件名。3. results文件夹不是打开即用先读懂输出文件族LISFLOOD-FP跑完之后results目录里的文件长得五花八门。新手最容易犯的错就是直接拿文本编辑器打开一个二进制文件看到一堆乱码后当场懵掉然后跑来群里问“模型是不是坏了”。其实这些文件都有固定结构只是需要按对应格式去读。3.1 results文件夹里到底会有什么按文件类型拆解以我用得最多的版本为例输出大致分为这几类xxx_WD_*.asc水深栅格ASCII格式时会生成一系列文件每个文件对应一个输出时间步。文件内容是网格水深值单位米。xxx_Qx_*.asc、xxx_Qy_*.ascx方向和y方向的单宽流量单位通常是m²/s。如果你关注洪水演进路径可以看这两个量的模。xxx_MF_*.asc动量通量或水流方向相关结果有些版本不输出。xxx.Qm这种小数点后缀的文件往往不是地形图而是边界出口或极值统计的时间序列具体要看文件头注释。xxx.res全局模拟日志里面记录的时间、水量、质量守恒误差都在这里这是第一步要看的文件而不是看水深图。注意不同版本的后缀名差异很大。比如旧版本可能用Qm代表轴向流量新版本用Qx/Qy。所以我每次拿到一个新编译版本都会先跑一个极小案例然后ls -la results/看生成哪些后缀再对照用户手册确认含义。这样能避免把Qx当初WD使用。3.2 二进制结果文件的结构与读取难点很多编译版本默认输出二进制格式可能是因为文件小、读写快。但二进制意味着不能直接用文本编辑器看。它的典型结构是开头若干字节是文件头包含行列数、时间戳、数据类型等信息之后才是纯数据体。不同版本的文件头长度和编码方式不一样有的甚至不写头只有裸浮点数组。我踩过的一个坑是从某台集群上拷贝回来的结果文件用Python的numpy.fromfile读出来后发现数据完全不对后来检查才知道是字节序问题。那台集群是little-endian但结果文件是按big-endian写的需要指定dtypef4才能正确读。如果你发现数据整体乱序先检查np.fromfile(..., dtypef4).byteswap()能不能恢复正常如果还是乱再看文件头里是否藏着网格信息。下面是一段我在读取ASCII水深文件时常用的脚本它可以直接读取LISFLOOD-FP输出的ASCII水深栅格import numpy as np def read_lisflood_ascii(path): with open(path, r) as f: header {} for _ in range(6): line f.readline().strip() key, value line.split() header[key] float(value) data np.loadtxt(f) return header, data如果你拿到的是二进制结果可以尝试import numpy as np def read_bin_flood(path, rows, cols, dtypef4): arr np.fromfile(path, dtypedtype) if arr.size ! rows * cols: arr arr.byteswap() return arr.reshape(rows, cols)注意这里的rows和cols通常可以从模型运行的日志或者DEM的头部拿到。如果不知道自己该填多少行可以试几个值直到arr.size能整除且reshape后画面合理。3.3 结果中的空间参考问题为什么你的水深图“漂浮”在错误位置LISFLOOD-FP输出的ASCII水深文件本身不包含投影坐标信息它只有一行一行的数值。如果你直接把水深文件拖进GIS软件软件会把它当成一个没有地理参考的普通格网叠加到底图上时永远对不上位置。我通常的做法是从原始DEM的.asc文件里复制头部信息再和水深数据合并写成一个新的ASCII或者直接用rasterio写GeoTIFF。下面是一段组合代码import numpy as np from osgeo import gdal def asc_to_geotiff(dem_asc, water_asc, out_tif): with open(dem_asc) as f: header {} for _ in range(6): k, v f.readline().split() header[k] float(v) cols, rows int(header[ncols]), int(header[nrows]) xll, yll header[xllcorner], header[yllcorner] cell header[cellsize] nodata header.get(NODATA_value, -9999) _, water read_lisflood_ascii(water_asc) driver gdal.GetDriverByName(GTiff) ds driver.Create(out_tif, cols, rows, 1, gdal.GDT_Float32) ds.SetGeoTransform((xll, cell, 0, yll rows * cell, 0, -cell)) band ds.GetRasterBand(1) band.WriteArray(water) band.SetNoDataValue(nodata) ds.FlushCache()这样生成的GeoTIFF才能和DEM、底图严格对齐后续做淹没面积统计、与实测水深对比都不会错位。4. 问题导向的完整排查链路当结果水深大到离谱我做了什么前面讲了很多预防措施但模型跑完后结果不对的情况依然常见。我想用一次真实事故的排查过程把整个链路完整走一遍你对这个概念会理解得更透。4.1 一次“水深300米”事故的排查过程去年帮朋友调试一个南方小流域的城市内涝模拟他用的DEM是直接从公开数据源下载的没做投影。第一次跑完一看结果最大水深307米连楼层最高的楼顶都淹掉了。他第一反应是“LISFLOOD-FP不适合城市洪水”我劝他别急先把单位链走一遍。排查第一步查看参数文件。sim_time设的是86400秒对应24小时没问题降雨强度写的是0.000036换算一下就是约130mm/hr虽然偏大但也不是完全不可能fpfric取0.03正常范围。于是先把参数排除。第二步查看DEM头部。打开dem_meter.asccellsize那一行写的是0.0009。这就是破绽。如果DEM单位是米cellsize不可能只有0.0009米这明显是度。我再一看投影信息果然是地理坐标系WGS84没有经过投影。第三步把DEM重投影到UTM 50N。重新查看cellsize变成大约91.7米。这个值就合理多了。此时我重新检查降雨强度它按m/s为单位0.000036换算成mm/hr是129.6虽然很大但对短时强降雨来说不是不可能。然后重跑模型最大水深变成1.8米和一个实测积水点深度吻合。这个case就这么解决了。整个排查过程其实只用了一个小时但核心点在于不要急着改模型参数先确认所有输入数据的空间单位和物理单位是否已经统一。很多时候问题就藏在DEM的cellsize和降雨强度的e-7量级上。4.2 另一个坑结果在边界处出现异常负水深单位排查完之后我后来又遇到一个更隐蔽的问题模拟结果整体合理但在下游边界附近出现一大片负水深。负水深在物理上是没有意义的这通常是初始条件或边界处理出了bug。那次是因为我在初始水深文件里给下游河道赋了一个高于地表的初始水位而河道边界只有一个流量过程线没有对应的水位—流量关系。模型在下游边界倒推水深时产生了震荡个别网格出现负值。解决办法是把初始水深置为0让水流自然演进如果必须要初始水位就在边界额外加一个水位过程线保证边界条件闭合。这类跟边界有关的坑靠肉眼很难看出来最好的办法是检查massint生成的质量守恒记录看模拟结束时总水量和入流量是否平衡。如果质量误差超过0.1%优先怀疑边界条件。4.3 如何验证你的结果可信除了质量守恒我还会从结果文件中提取关键断面的最大水深、淹没时间、流量峰值和实测水文数据对比。如果没有实测数据就做一个水量平衡估算总降水量乘以流域面积减去下游出口累计径流量和入渗量应该等于模拟结束时的地表积水增量。只要这个误差在5%以内我就认为结果可信。你可以用xxx_Qm或xxx_Qy出口文件读取出口流量过程线然后手动积分得到总径流量。这个过程虽然是粗验但对发现单位和参数错误非常有效。我见过有人在河道入流边界写错了单位结果出口流量是正常值的10倍一眼就能识别。5. 多年跑模型后我自己的后处理脚本和避坑清单后面这部分我分享几个每次跑模型都离不开的小工具和习惯。它们不算高深但能帮你省下大量重复劳动。5.1 Python可视化后处理三件套我常用的第一个脚本是把所有ASCII水深文件读取成一个三维数组然后做逐帧最大水深合成。第二个脚本是把结果转成GeoTIFF方便叠加到地图上。第三个脚本是批量绘制某几个网格单元的逐时水深曲线。核心代码如下import glob import numpy as np def max_water_depth(root, output_npy): files sorted(glob.glob(f{root}_WD_*.asc)) arrays [] for f in files: header, data read_lisflood_ascii(f) arrays.append(data) stack np.stack(arrays, axis0) max_depth stack.max(axis0) np.save(output_npy, max_depth)这段代码的精髓是直接利用read_lisflood_ascii把所有结果合成一个max_depth一份报告里放一张最大水深图就够了。5.2 放进项目文件夹的README模板我习惯在每个模拟项目里放一个README.md开头就是一张参数、数据、结果的对照表。表格大概长这样项目文件名/路径单位备注DEMinput/dem_meter.asc米已投影到UTM50N降雨input/rain_hourly.csvmm/hr已转换为m/s边界入流input/q_in.csvm³/s每15分钟一个值参数文件params.par-fpfric0.03结果目录output/res-二进制输出6小时间隔这样一个表格等三个月后你回来看项目或者同事接手你的模拟都不需要再从头猜单位。这比我踩过的那些坑更值钱让错误不再重复。5.3 最后补充永远保留中间版本模拟过程中会反复修改参数、修正单位。我强烈建议每改一次参数文件就整个复制一份命名为v01、v02、v03而不是在原文件上覆盖。因为经常出现同一个参数尝试了十种组合最后你发现最合理的反而是最早的v01但如果没有版本管理你只能重新试错。我还会把每次运行后results文件夹里的*.res日志文件按版本号归档。日志里记录了实际使用的时间步长、水量误差、文件输出情况这些信息比结果图本身更能帮助回溯问题。这是我跑LISFLOOD-FP几年下来最深的体会模型本身不算难学难的是把数据和参数管理成可追溯、可复现的状态。单位检查、参数模板、结果解析脚本这套组合拳打好了洪水模拟的坑基本能避开八成。剩下的两成就交给实测数据和运气吧。