LISFLOOD-FP避坑指南:从DEM单位到结果解析的完整排查链路
2026/9/21 1:23:49 网站建设 项目流程

坐标系里的“度”混进模型后,水深直接冲到几百米。这不是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³/s,L/s要除以1000。问题在于,当你面对一个几千个时段的边界文件时,很容易只在文件开头改了一个单位,后面时段全是原始值,模型的入流过程整体偏大或偏小,下渗、淹没范围全都会跟着错。

1.3 我的单位自检五步法

踩了几次坑之后,我给自己定了一套单位自检流程,每次模拟前都走一遍:

  1. 把所有外部输入文件列一张表,包括DEM、降雨、入流、初始水深,每项标注原始单位。
  2. 统一转换成SI单位,并把转换后的文件单独放在一个units_converted目录里,绝不改动原始数据。
  3. 在.par参数文件的注释部分(如果版本支持)写清每个输入文件对应的单位。
  4. 写一个简单的Python脚本,读入所有转换后文件的位置和全局统计值,检查最小值、最大值是否在合理范围,比如降雨强度在1e-8到1e-4 m/s之间,入流在0到几百m³/s之间。
  5. 跑一个理想化案例,比如在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_timeinitial_water_depthboundary_file等字段,具体以你编译的版本手册为准。我的经验是,拿到一个新版本先跑通这个最小案例,再逐步加复杂功能,别直接上大流域,否则一旦报错很难判断是参数文件格式还是物理条件的问题。

2.2 参数字段与单位的对应关系

我用一张表整理一下常见字段和单位,这张表我贴在自己工位对面:

字段名含义常用单位备注
DEMfile高程数据文件路径必须为投影坐标系
sim_time模拟总时长若写成小时会缩小3600倍
initial_tstep初始时间步长受CFL约束,宁小勿大
resstep结果输出步长太大会漏掉洪水过程细节
massint质量统计间隔守恒诊断用
fpfric洪泛区曼宁n-无量纲,常取0.02-0.15
fpfricfile空间曼宁n栅格-栅格数值对应n值
rainfallfile降雨过程文件m/s不是mm/hr
qfile边界流量文件m³/s可以是过程线

很多人会在sim_time那里直接写7200,觉得单位无所谓,结果预期是两小时,模型却只跑了两秒。也有人在resstep里写了300秒,但sim_time2小时,导致只有一个结果帧,输出文件小得可怜。这些数字看起来不起眼,却直接决定结果文件有没有意义。

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_*.ascxxx_Qy_*.asc:x方向和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写的,需要指定dtype='>f4'才能正确读。如果你发现数据整体乱序,先检查np.fromfile(..., dtype='f4').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, dtype='>f4'): arr = np.fromfile(path, dtype=dtype) if arr.size != rows * cols: arr = arr.byteswap() return arr.reshape(rows, cols)

注意这里的rowscols通常可以从模型运行的日志或者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_Qmxxx_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, axis=0) max_depth = stack.max(axis=0) 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-fpfric=0.03
结果目录output/res-二进制输出,6小时间隔

这样一个表格,等三个月后你回来看项目,或者同事接手你的模拟,都不需要再从头猜单位。这比我踩过的那些坑更值钱:让错误不再重复。

5.3 最后补充:永远保留中间版本

模拟过程中会反复修改参数、修正单位。我强烈建议每改一次参数文件就整个复制一份,命名为v01v02v03,而不是在原文件上覆盖。因为经常出现同一个参数,尝试了十种组合,最后你发现最合理的反而是最早的v01,但如果没有版本管理,你只能重新试错。

我还会把每次运行后results文件夹里的*.res日志文件按版本号归档。日志里记录了实际使用的时间步长、水量误差、文件输出情况,这些信息比结果图本身更能帮助回溯问题。

这是我跑LISFLOOD-FP几年下来最深的体会:模型本身不算难学,难的是把数据和参数管理成可追溯、可复现的状态。单位检查、参数模板、结果解析脚本这套组合拳打好了,洪水模拟的坑基本能避开八成。剩下的两成,就交给实测数据和运气吧。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询