企业官网建设流程全解析

简介：这是一份面向水利工程、水文水资源专业研究人员与防洪排涝规划从业者的研究报告，聚焦江苏省兴化市平旺湖退圩还湖工程，借助MIKE21FM二维水动力模型，对比分析20年一遇洪水条件下工程实施前后湖区水动力条件、河道槽蓄能力与洪水特征值的变化，可为同类退圩还湖方案论证与湖泊水环境、水生态修复研究提供参考依据。压缩包共1个文件，为21KB的docx文档，正文系统梳理了研究区概况、退圩还湖工程概况、模型构建与数据方法、结果分析等完整章节，并配有网格设置、地形重塑、排泥场布置等关键内容说明。文中涉及退圩还湖对调蓄库容、淹没水深、蓄洪量以及水体南北往复流动的影响分析，也给出了MIKE21FM模型在流速流场与洪水淹没分布模拟中的具体应用思路，便于读者借鉴建模流程与评价指标。目前已有96人学习，适合需要掌握二维水动力建模方法或从事湖泊治理、防洪评价工作的读者参考。

1. 退圩还湖洪水特征分析为什么绕不开 MIKE21FM

一个退圩还湖方案做完，湖面面积增加十几平方公里，很多人会默认洪峰水位必然明显下降。真跑一遍模型才发现，削掉的洪峰可能只有几厘米，局部水位甚至不降反升。问题通常不在水动力方程，而在网格尺度、圩区开挖高程和糙率分区这些前处理环节。退圩还湖的洪水效应本质上是“水面变大、调蓄容积变大、阻力变小”三个作用叠加，但它们对最高水位、淹没面积、流速场的贡献方向并不一致，靠经验判断很容易翻车。MIKE21FM（MIKE 21 Flow Model FM）用非结构三角网格做二维浅水方程求解，天然适配“岸线曲折、圩堤纵横、闸站密布”的湖泊场景，干湿动边界能直接处理退圩前后水陆交替。这篇讲的是从地形改造到情景对比的完整链条：网格怎么控尺度，退圩区高程怎么改，边界和干湿参数怎么给，最高水位、淹没历时、最大流速这些指标怎么从结果文件里提出来，以及结果明显不合理时先查哪几处。面向做湖泊治理、防洪评价和水动力数值模拟的工程师，也适合刚接触 MIKE 体系、想跑通第一个退圩还湖算例的人。

2. MIKE21FM 非结构网格与退圩区地形概化

2.1 非结构网格生成与尺度控制

MIKE21FM 的网格走的是 MIKE Zero 里的 Mesh Generator，导入岸线 shp 和地形散点 xyz 之后生成三角网格，网格文件是 .mesh。湖区这种大水面加复杂岸线的场景，用结构网格很难贴合岸线和圩堤，非结构网格可以做到“该密的地方密、该疏的地方疏”。

尺度控制的经验值：开阔湖面 150~300 m，岸线附近和主要行洪通道 30~50 m，圩堤、闸口、涵洞附近 10~20 m，入湖河道 5~15 m。判断密不密的办法不是看数字，而是看网格能不能把关键地形表达出来——一条宽度 8 m 的圩堤，如果网格尺度是 30 m，它几乎不存在于模型里。

三角形质量比尺度更容易被忽略。MIKE21FM 对网格内角敏感，建议内角控制在 30°~130°，避免狭长三角形和钝角单元。狭长单元会让对流项计算产生伪振荡，表现为局部水位出现不合理的尖峰。跑之前用 Mesh Generator 的质量检查工具过一遍，比跑完之后对着异常水位发愁省事得多。

还有一个容易踩的坑：现状方案和退圩还湖方案分别生成两套网格。两套网格的节点位置、单元数量都不同，算出来的水位差异里混进了离散误差，最后根本说不清是工程效应还是网格效应。常见做法是只做一套网格，计算域覆盖包含圩区陆地在内的整个区域，靠干湿处理区分水陆；两个情景共用一个 .mesh，只换地形文件。这样对比出来的差值才是干净的。

2.2 退圩还湖前后地形高程的差异处理

地形是退圩还湖分析里最核心的输入。现状情景里，圩区是农田或村庄，田面高程通常比湖底高 2~4 m，圩堤堤顶再高出 1~3 m。退圩还湖情景里，这部分区域开挖到设计湖底高程，同时拆除或部分拆除圩堤。

最省事的做法是在 GIS 里用退圩区多边形裁切地形散点，把多边形内的节点高程统一替换成设计湖底高程。这个做法能跑，但会在开挖边界处形成垂直台阶，MIKE21FM 在这种地形上容易出现反复干湿震荡和负水深报错。合理做法是在圩堤内侧留一条过渡带，带内高程从原地面线性过渡到设计湖底。

import numpy as np import geopandas as gpd import shapely pts = np.loadtxt("terrain_current.xyz") # N x 3，列顺序 x y z，单位 m zone = gpd.read_file("tuwei_zone.shp") # 退圩区范围，坐标系需与地形一致 poly = shapely.union_all(zone.geometry.values) x, y = pts[:, 0], pts[:, 1] inside = shapely.contains_xy(poly, x, y) # 判断点是否落在圩区内 d_edge = shapely.distance(shapely.points(x, y), poly.exterior) # 点到圩堤内边线的距离 band = 80.0 # 过渡带宽度 m w = np.clip(d_edge / band, 0.0, 1.0) # 0 在堤脚，1 在核心区 bed = 4.5 # 设计湖底高程 m z_new = pts[:, 2].copy() z_new[inside] = (1 - w[inside]) * pts[inside, 2] + w[inside] * bed np.savetxt("terrain_return_lake.xyz", np.column_stack([x, y, z_new]), fmt="%.3f")

这段代码的逻辑是：圩区外的节点高程原样保留，圩区内从堤脚的原地面高程线性过渡到核心区的设计湖底高程。band取 80 m 意味着过渡带跨越 80 m 水平距离，比一个网格单元略大，数值上不会形成突变；如果局部网格只有 10 m，可以把这个值降到 30~50 m。bed必须用工程设计给出的湖底高程，不能凭经验“降个一两米”——退圩还湖的调蓄容积增加量几乎完全由这个值决定，差 0.5 m 就能让洪峰水位变化翻倍。

圩堤的拆除也要落到地形上。保留堤段把高程抬到堤顶，拆除堤段把高程压到湖底，并且保证拆除段在网格里有至少 3~4 个单元跨越，否则水流会从“数字上的堤”旁边绕过去。

2.3 圩堤、闸站与涵洞的概化方式

圩堤在 MIKE21FM 里有三种表达路径。一是抬高地形，用网格单元的高程表达堤身，实现简单、鲁棒性好，缺点是堤顶宽度小于网格尺度时会被“抹平”。二是在 mesh 里加 Dike 结构，相当于把一维堤线嵌进二维网格，阻水效果准确，但需要额外定义堤线几何和堤顶高程。三是用 MIKE 21 的 Structure 模块做堰流、闸孔，适合需要模拟闸门调度的情况。

闸站通常按功能分开处理。排涝闸、节制闸这类有调度规则的，用 Structure 里的闸孔公式，给定闸底高程、闸宽和开度过程。涵洞用 Culvert，给定断面尺寸和进出口高程。如果分析重点只是湖区的洪水特征，闸站细节影响有限，可以把闸址处的过流能力折算成一条水位—流量关系曲线，加在边界或内部连接上，这也是防洪评价里常见的简化。

不管用哪种方式，都要保证概化后的过流能力与设计资料一致。一个检查办法是单独跑一个恒定流算例，固定上下游水位，看计算流量与设计流量差多少，偏差超过 10% 就回去改参数。

2.4 糙率分区与干湿判定

退圩还湖会显著改变下垫面，糙率必须跟着改。MIKE21FM 里糙率可以用曼宁数 M（m^(1/3)/s）或曼宁 n，两者互为倒数，PFS 里选 Manning number 时要看清单位。常用的分区取值：

分区可以按 mesh 的 region 赋值，也可以叠加一张 dfs2 糙率场。现状方案里圩区按农田取值，退圩方案里改成湖面或浅水沼泽，这一项对最高水位的影响经常和地形开挖量级相当，不要只改地形不改糙率。

干湿参数决定水陆边界怎么移动。Drying depth 建议 0.005~0.02 m，取值过小会让节点在干湿之间反复切换，日志里出现大量警告，计算也慢；取值过大则淹没范围被系统性低估。Flooding depth 和 Wetting depth 一般取 0.05~0.10 m，且两者相等，避免出现“淹没但不流动”的异常单元。

3. 退圩还湖算例的边界条件与 MIKE21FM 参数配置

3.1 开边界水位过程与流量过程怎么选

湖区洪水特征分析的边界通常由三部分组成：入湖河道给流量过程，出湖口给水位过程，区间来水用降雨径流模型折算后作为旁侧入流。入湖流量用 dfs0 时间序列，出湖口如果受外江水位顶托控制，也给 dfs0 水位过程；如果外江水位影响小，可以给水位—流量关系。

用 mikeio 生成一个开边界水位文件：

import numpy as np import pandas as pd import mikeio t = pd.date_range("2020-07-01 00:00", "2020-08-01 00:00", freq="1h") i = np.arange(len(t)) wl = 8.5 + 3.2 * np.exp(-((i - 72) / 60.0) ** 2) # 示意洪水位过程，单位 m ds = mikeio.Dataset( data=[pd.DataFrame({"Water level": wl}, index=t)], items=[mikeio.ItemInfo("Water level", mikeio.EUMType.Water_Level)], ) ds.to_dfs0("open_boundary_wl.dfs0")

逻辑很简单：构造一个带洪峰的逐小时水位序列，写成 MIKE 的 dfs0 格式。freq="1h"要和 mdf 里的时间步输出设置一致，否则插值会引入误差。8.5是起涨水位，3.2是洪峰涨幅，这两个数必须来自实测或设计洪水过程，不能随便写。mikeio 的 Dataset 构造签名在不同大版本间有差异，如果报参数错误，先print(mikeio.__version__)再对照官方示例调整，读到 Dataset 之后取变量名这一步是通用的。

入湖流量边界同理，把 item 换成mikeio.EUMType.Discharge，单位 m³/s。注意流量过程的总水量要和设计洪水的洪量对得上，跑完之后拿 log 里的累计入流核对一遍，差太多说明序列本身有问题。

3.2 时间步长与库朗数约束

MIKE21FM 用自适应时间步长，由库朗数控制。库朗数定义是流速乘时间步长除以网格尺度，建议取值范围 0.7~0.9。MaxTimeStep 给一个上限，防止在流速极小的区域步长无限放大。

第一步跑完先看 log 里的实际时间步长分布。如果步长被压到 0.1 s 还在报库朗数超限，说明某处网格过小或者出现了异常流速，先查加密区周围的地形是不是有孤立的深坑或尖峰。

3.3 干湿水深、涡黏与风场设置

干湿参数按上一节的推荐值给。除此之外还有两个开关值得注意：一是 Flood and Dry 必须打开，退圩还湖分析里水陆边界是核心；二是涡黏模型，湖区推荐用 Smagorinsky，系数 0.28，比恒定涡黏更适应流速梯度变化大的区域。如果计算域内有明显的回流区、弯道，恒定涡黏容易把涡抹平，算出来的流速场偏平滑。

风场在大湖面上不能完全忽略。洪水期如果恰好遇到大风，风生流会把水位在迎风岸堆高几十厘米，这和退圩还湖的效应量级接近。有实测风速风向资料时，按 dfs2 给空间分布的风场；没有的话至少给一个常数风速做敏感性测试，看看结果对风是否敏感，敏感就在报告里说明。

3.4 用脚本批量生成和提交算例

退圩还湖分析至少要跑现状、退圩后两个情景，往往还要加一个“只退不挖”的敏感性情景。手工在界面里改参数容易出错，也难复现。可行做法是先用 MIKE Zero 界面生成一个能跑通的 mdf 作为模板，再用脚本做文本替换。

import re from pathlib import Path tpl = Path("template.m21fm").read_text(encoding="utf-8") cases = { "S0_baseline": {"MANNING": "25", "MAXDT": "20"}, # 现状：农田糙率 "S1_return_lake": {"MANNING": "40", "MAXDT": "20"}, # 退圩：湖面糙率 "S2_no_dredge": {"MANNING": "30", "MAXDT": "20"}, # 只退不挖 } for name, cfg in cases.items(): txt = tpl txt = re.sub(r"(ManningNumber\s*=\s*)\S+", r"\g<1>" + cfg["MANNING"], txt) txt = re.sub(r"(MaxTimeStep\s*=\s*)\S+", r"\g<1>" + cfg["MAXDT"], txt) Path(f"{name}.m21fm").write_text(txt, encoding="utf-8") print(name, "written")

这里用正则直接替换 PFS 文本里的字段值，好处是不依赖任何特定版本的 Python 接口，坏处是字段名和缩进必须和模板完全一致。实际用的时候先grep -n "ManningNumber" template.m21fm确认字段存在，再写正则。地形文件路径、结果文件路径、边界文件路径也可以一起替换，把三个情景的输出目录分开。

提交计算用命令行引擎，路径按本机安装目录调整：

set EXE="C:\Program Files\DHI\MIKE Zero\2023\bin\x64\Mike21FM.exe" for %C in (S0_baseline S1_return_lake S2_no_dredge) do %EXE% -x %C.m21fm

-x表示无界面运行，适合放在后台批量跑。三个情景如果共用一个 mesh，只有地形、糙率和边界不同，跑完的差异就可以直接归因到工程本身。

4. 退圩还湖情景对比：MIKE21FM 洪水特征指标提取

4.1 情景设置与对比逻辑

情景设置要围绕“想回答什么问题”来定，而不是越全越好。退圩还湖分析通常关心三类问题：退圩后最高水位降多少、淹没范围怎么变、流速场有没有出现新的高风险区。对应的情景至少要有现状和退圩后两个，再加一个只退不挖的敏感性情景用来分离“水面扩大”和“容积增加”的贡献。

S2 是关键。很多项目只做 S0 和 S1，得出“最高水位下降 0.15 m”就收工，但说不清这 0.15 m 里有多少来自开挖、多少来自拆堤。加上 S2 之后，S1 与 S2 的差就是开挖的净贡献，这个数在方案论证时比总量更有说服力。

4.2 dfsu 结果的读取与统计

MIKE21FM 的结果是 dfsu 格式，节点或单元上的时间序列。用 mikeio 读取并统计：

import numpy as np import mikeio ds = mikeio.read("S1_return_lake.dfsu") # 结果文件，含 H、U、V、Total water depth d = ds["Total water depth"].to_numpy() # 形状 (nt, n_elem)，水深 m u = ds["U-velocity"].to_numpy() # 单元中心 x 向流速 m/s v = ds["V-velocity"].to_numpy() area = ds.geometry.get_element_area() # 每个单元的面积 m^2 dt_h = (ds.time[1] - ds.time[0]).total_seconds() / 3600.0 # 输出间隔，小时 thr = 0.30 # 淹没判定水深阈值 m inundated = d > thr # 布尔矩阵 (nt, n_elem) area_max = area[inundated.any(axis=0)].sum() # 计算期内最大淹没面积 m^2 duration = inundated.sum(axis=0) * dt_h # 每个单元的淹没历时 h spd = np.sqrt(u ** 2 + v ** 2) spd_max = np.nanmax(spd, axis=0) # 每个单元的最大流速 h_max = np.nanmax(d, axis=0) # 每个单元的最大水深 print(f"最大淹没面积 {area_max/1e6:.2f} km²") print(f"最大流速 {spd_max.max():.2f} m/s")

关键在thr这个阈值。防洪评价里常用 0.3 m，人员风险分析用 0.5 m，湿地和生态淹没用 0.1 m。阈值一变，淹没面积能差出百分之十几，报告里必须写清楚用的是什么标准。duration用输出间隔乘以出现次数来近似，输出间隔给到 1 h 时，历时统计的误差在 ±1 h 以内，够用；如果要分析短历时内涝，间隔要缩到 10 min。

spd是单元中心流速，MIKE21FM 也可以输出节点上的流速，两者在浅水边缘差异较大。做流速场图时建议用单元值，因为单元中心不会落在干单元上。

4.3 洪水特征指标的对比与解读

把三个情景的指标并排放在一张表里，对比才有意义。下面这张表是格式示意，数值需要按自己的算例填：

解读时要抓住两点。第一，退圩还湖后最高水位下降、淹没面积增大是正常现象，因为工程本身把洪水从“高水位窄水面”变成“低水位宽水面”，淹没面积增加不等于防洪能力下降，关键看新增淹没区是不是规划中的退圩区。第二，S2 与 S0 的水位差往往很小，说明拆堤本身对调蓄的贡献有限，真正起作用的是开挖容积。如果 S1 的水位降幅也小，先回头看开挖高程有没有按设计值给、糙率有没有从农田改成湖面，这两处是最常见的“改了个寂寞”。

流速指标容易被忽略。退圩后过水断面变大，主流区流速一般下降，但圩堤拆除口门附近可能出现局部流速增大，形成新的冲刷风险点。把spd_max画成等值面图，对照口门位置检查一遍。

4.4 结果合理性检查的四个入口

第一看水量平衡。log 里的 mass balance error 应该小于 1%，超过 5% 说明边界流量、降雨或干湿处理有问题，后面的水位结论都不用信。

第二看库朗数和时间步长。如果实际步长一直在下限附近徘徊，说明某处网格或地形异常，查加密区的孤立深坑和尖峰。

第三看边界拟合。把计算得到的出湖口流量过程与边界给定的水位过程对照，检查是否自洽；入湖总水量与边界给定序列的积分对比，偏差应小于 2%。

第四看网格收敛性。把湖面网格从 200 m 加密到 100 m 重跑，最高水位变化应小于 5 cm。如果变化超过 10 cm，说明当前网格还没收敛，对比结论不可靠。这个检查成本不高，但能挡掉大部分“数字好看、结论站不住”的情况。

5. MIKE21FM 糙率率定与退圩情景结果验证的几个技巧

参数率定要按影响量级排序，不要一上来就调一堆参数。糙率对水位的影响最大，先率定曼宁 M；涡黏主要影响流速场，对水位影响次之，放在第二步；干湿水深影响淹没边界，放最后。率定目标是实测站水位过程的纳什效率系数、均方根误差和洪峰水位误差，工程报告的常见门槛是 NSE 大于 0.85、洪峰水位误差小于 0.10 m。

单站率定很容易过拟合。至少选 2~3 个站，且分布要覆盖湖心、入湖口和出湖口，避免所有站都在同一片水域。如果湖心站拟合很好而出湖口站偏差大，问题多半出在出湖口边界条件或口门附近的网格概化，而不是糙率。

退圩还湖情景本身没有实测资料可验证，只能靠现状年率定出的参数外推。这里有一条底线：除糙率分区随下垫面变化外，其他参数在情景之间必须保持一致。为了让退圩后的水位降得更“好看”而偷偷调小涡黏或干水深，属于典型的结果导向操作，评审时一问就露馅。

把计算结果插值到实测站位置时，直接取最近的网格单元会带来几十米的错位误差，在岸线附近尤其明显。用单元中心坐标做线性插值更稳：

import numpy as np from scipy.interpolate import griddata coords = ds.geometry.element_coordinates[:, :2] # 单元中心 x, y h_last = ds["Surface elevation"].to_numpy()[-1] # 末时刻水位 sx, sy = 112.4832, 30.6715 # 实测站坐标，需与网格同坐标系 h_station = griddata(coords, h_last, (sx, sy), method="linear") print(f"站点插值水位 {h_station:.3f} m")

method="linear"在站点落到凸包外时会返回 nan，这时改用"nearest"并检查站点是否真的在计算域内。插值前一定要确认站点坐标和网格用的是同一套坐标系，投影不一致会得到一个完全离谱的水位值，而且不会报错。

最后是评价指标的计算，NSE 的写法固定，注意观测和模拟的序列要先按时间对齐、剔除缺测：

def nse(obs, sim): obs, sim = np.asarray(obs, float), np.asarray(sim, float) m = ~(np.isnan(obs) | np.isnan(sim)) obs, sim = obs[m], sim[m] return 1 - np.sum((obs - sim) ** 2) / np.sum((obs - obs.mean()) ** 2)

率定期 NSE 能到 0.9 而验证期掉到 0.7，先别急着改糙率，去看验证年的边界流量是不是用了另一套产汇流参数——退圩还湖洪水分析里，边界条件的年际不一致比模型参数问题更常见。

本文还有配套的精品资源，点击获取