简介面向水利、水文及计算流体方向学习者这份资源演示了基于 MATLAB 的一维河道水动力学建模流程核心围绕圣维南方程组与 Preissmann 离散法展开可用于理解流速、水位和流量变化模拟以及上、下游边界条件的设定方法。压缩包共5个文件以 .m 主程序、.mat 数据文件和 .txt 说明与地形资料为主整体约401KB结构精简。已有59人学习下载。内容包含完整建模主程序、河道地形数据、输入参数与边界条件示例并配有使用说明与地形资料文本对照代码与说明可掌握一维水动力模型从数据准备、算法离散到结果解读的完整思路对洪水预警、水资源管理和水利工程设计具有直接参考价值。1. 一维河道水动力学建模别一上来就追二维先把圣维南方程组跑稳很多做防洪预报、水环境影响评估的工程师拿到一个河道项目的第一反应是上二维模型。但实际干过几个项目就知道几公里甚至几十公里的天然河道二维网格动辄十万级率定参数调到怀疑人生最后算出来的水位可能还不如一维模型加一套靠谱的糙率。一维河道水动力学建模核心就是沿着河道中心线求解圣维南方程组把流量和水位沿程算出来。它的价值不在于精度碾压二维而是用最小的数据成本和计算成本解决“洪水从哪来、水位到多高、流量怎么演进”这类最刚需的问题。适合刚接触水动力模拟的从业者也适合经费有限、实测断面数据不全的中小型项目。2. 圣维南方程组怎么简化才合理连续方程、动量方程与三种波的取舍2.1 一维河道模型到底在算什么两个方程、两个未知数的闭式系统一维河道水动力模型的控制方程是法国科学家圣维南在1871年提出的浅水方程组由连续方程和动量方程组成。连续方程表达质量守恒河道里某一段水体的水量变化等于上游流进来减去下游流出去还要考虑旁侧入流。动量方程表达牛顿第二定律在水体上的体现惯性力、压力梯度力、重力沿坡降的分量以及河床和岸壁的阻力必须达到平衡。写成常见的守恒形式# 连续方程∂A/∂t ∂Q/∂x q_lateral # 动量方程∂Q/∂t ∂(Q²/A)/∂x g·A·∂Z/∂x g·A·(S0 - Sf) # A: 过水断面面积 (m²) # Q: 流量 (m³/s) # Z: 水位 (m) # q_lateral: 单位长度旁侧入流 (m²/s) # S0: 河床底坡 (m/m) # Sf: 摩阻坡降用曼宁公式计算 Sf n²·Q·|Q| / (A²·R^(4/3))这里有两个未知量——流量 Q 和水位 Z或者断面面积 A两个方程构成闭式系统所以给定边界条件和初始条件理论上就能求出沿程各断面的水位和流量随时间的演变。实际编码时常用水位 Z 和流量 Q 做未知量而不是用断面面积 A因为天然河道断面不规则A 和 Z 的关系是非线性的函数表用 Z 做变量更方便在断面数据表里插值。2.2 运动波、扩散波与动力波不是越完整越好要看你的问题尺度完整圣维南方程组包含所有物理项数值求解难度最大对数据要求也最高。实际工程里会根据问题的物理特性做简化。如果河道底坡陡、洪水波拉长、回水影响小比如山区性河流或城市排水明渠可以把惯性项和压力项都丢掉只保留重力项与摩阻项平衡这就是运动波近似。运动波模型的核心假设是流量与水位之间存在单值关系也就是流量只取决于当地水深和坡降下游的水位根本影响不到上游所以计算简单、速度快适合做长时段洪水演进。如果下游有闸坝、支流汇入、潮汐顶托这类回水影响运动波就失效了因为水位流量关系被迫变成绳套曲线同一水位对应两个流量。这时要保留压力梯度项这就成了扩散波近似。扩散波能够刻画回水影响但依然忽略惯性项适用于坡度适中、回水问题不严重的河道。如果要做溃坝洪水、闸门骤启骤停这类惯性效应显著的场景必须用完整动力波也就是求解完整圣维南方程组但这时的参数率定难度和数值稳定性问题都会陡增。2.3 离散化的数学本质从偏微分方程到差分方程的那一步决定了成败偏微分方程没法直接让计算机求解必须把连续空间和时间离散成网格节点用差分或有限体积逼近导数项。常见做法是空间上沿河道按断面间距划分网格时间上按固定步长推进。这里有三个关键选择空间差分格式、时间积分格式和网格间距。# 圣维南方程组的通用离散框架 # 空间导数三点中心差分或迎风差分 dQ_dx (Q[i1] - Q[i-1]) / (2 * dx) # 中心差分二阶精度 dQ_dx_upwind (Q[i] - Q[i-1]) / dx # 迎风差分适合单向传播 # 时间导数显式或隐式推进 # 显式dQ_dt (Q_new[i] - Q_old[i]) / dt # 隐式dQ_dt (Q_new[i] - Q_old[i]) / dt右端项用 Q_new 取值显式格式实现简单但稳定性受限于 CFL 条件水流在一个时间步内不能跨过超过一个空间网格。也就是说 dt 必须小于 dx 除以波速否则数值解发散。天然河道里波速可能到每秒几米如果断面间距是 100 米时间步长最大也就 20 秒左右。隐式格式没有这个限制时间步长可以放宽到几十倍代价是每步要解一个大型方程组。工程模型普遍选隐式格式原因很现实天然河道断面数据粗糙没必要在时间分辨率上较劲把计算步长放到分钟级一天的水流演进只要几百步就能算完。3. 用 Python 手写一维求解器Preissmann 四点隐式格式与稀疏矩阵组装3.1 断面数据怎么组织河道不是一根线是一串断面动手写代码前先把河道数据整理成稳定、可索引的结构。一维模型中的“河道”在经济意义上是一条中心线但实际水力计算需要的是每个断面桩号处的水位-面积-湿周关系表。野外测量通常给的是断面上若干测点的起点距和高程需要预处理成水深-面积和水深-水力半径的关系函数。# 断面数据处理从散点测量数据到水力要素表 import numpy as np from scipy.interpolate import interp1d # 原始测量点起点距 x 和河底高程 z x_pts np.array([0.0, 5.0, 10.0, 15.0, 20.0]) z_pts np.array([102.5, 101.2, 100.1, 101.0, 102.3]) def cross_section_properties(z_level, x_pts, z_pts): # 计算给定水位下的过水面积 A 和湿周 P # 思路水位以下的测点参与计算水位以上的点用线性插值找水面交点 if z_level np.min(z_pts): return 0.0, 0.0 # 分段线性积分的简化实现 area 0.0 wetted_perimeter 0.0 for i in range(len(x_pts) - 1): z_left z_pts[i] z_right z_pts[i1] if z_left z_level or z_right z_level: # 三角形或梯形面积 h_left max(z_level - z_left, 0.0) h_right max(z_level - z_right, 0.0) width x_pts[i1] - x_pts[i] area 0.5 * (h_left h_right) * width # 湿周只计算淹没部分的底边长 wet_frac max(h_left / (h_left h_right 1e-9), 0.0) wetted_perimeter wet_frac * np.hypot(width, z_right - z_left) return area, wetted_perimeter # 生成水位-面积插值表对每个断面预计算一组离散水位对应的面积 water_levels np.linspace(np.min(z_pts), np.max(z_pts) 2.0, 50) areas [cross_section_properties(zl, x_pts, z_pts)[0] for zl in water_levels] A_interp interp1d(water_levels, areas, fill_valueextrapolate)这段代码的核心思想是先对每个野外实测断面构建一张水位-面积关系表之后求解器在迭代中通过插值函数快速获取任一水位对应的过水面积而不用每次重新积分。参数说明水位数组的上下限决定了模型的适用水位范围上限至少要覆盖设计洪水位以上 1 米否则高水位时插值外延误差不可控。实测断面中如果有明显的滩地或主槽分界最好在 z_pts 里保留地形突变点否则湿周计算会出现跳变。3.2 Preissmann 四点隐式格式的离散模板半个网格处的空间差商Preissmann 格式是一维河道模型中使用最广泛的隐式格式1961 年由 Preissmann 提出特点是空间差商取在两个相邻网格点之间的中间位置时间导数项也取同一位置。这种“半网格”处理让格式在空间上有二阶精度同时保持了隐式格式的无条件稳定性——至少理论上是这样实际中时间步长太大一样会振荡。# Preissmann 四点格式的离散系数模板 # 设变量 f 代表水位 Z 或流量 Q在网格点 i 和 i1 之间取加权平均 # 权重系数 θ控制时间加权工程上通常取 θ 0.6~0.7 theta 0.65 # 空间导数在 n1/2 时刻的离散形式 dZ_dx theta * (Z[i1] - Z[i]) / dx (1 - theta) * (Z_old[i1] - Z_old[i]) / dx # 时间导数在 i1/2 位置的离散形式 dZ_dt 0.5 * ((Z[i1] - Z_old[i1]) / dt (Z[i] - Z_old[i]) / dt) # 连续性方程的最终离散形式 # 0.5 * ((A[i1] A[i]) - (A_old[i1] A_old[i])) / dt # theta * (Q[i1] - Q[i]) / dx # (1 - theta) * (Q_old[i1] - Q_old[i]) / dx # q_lateral (旁侧入流项)θ 参数的物理含义是当前时刻和上一时刻的权重分配。θ 1 时是全隐式虽然最稳定但数值扩散严重会把洪水波拉平θ 0 时退化为显式稳定性条件又回来了。工程经验一般取 0.6 到 0.7折中稳定性和精度。这段代码展示的是单个网格对的离散实际求解时会把所有网格点的离散方程联立成一个矩阵。3.3 组装稀疏矩阵并推进求解从牛顿迭代到稀疏LU分解Preissmann 格式每个网格点会生成两个方程一个来自离散后的连续方程一个来自动量方程。如果有 N 个断面就有 2N 个未知量每个断面的 Z 和 Q对应 2(N-1) 个方程剩余两个方程由上下游边界条件补齐。这个方程组是高度非线性的因为断面面积 A 和水力半径 R 是水位 Z 的强非线性函数所以每个时间步内都要做牛顿迭代。# 完整求解器骨架牛顿迭代 稀疏矩阵求解 from scipy.sparse import diags from scipy.sparse.linalg import spsolve import numpy as np def solve_one_step(Z_old, Q_old, boundary_up, boundary_down, sections, dt, dx, theta0.65): n len(Z_old) # 未知量排列先所有断面的 Z再所有断面的 Q共 2n 个 x np.concatenate([Z_old, Q_old]) for iter in range(20): # 牛顿迭代上限 # 组装残差向量 F 和雅可比矩阵 J F np.zeros(2 * n) J np.zeros((2 * n, 2 * n)) for i in range(n - 1): # 连续方程残差第一个方程 A_i sections[i].area(Z_old[i]) A_ip1 sections[i1].area(Z_old[i1]) A_new_i sections[i].area(x[i]) A_new_ip1 sections[i1].area(x[i1]) F[i] (0.5 * ((A_new_ip1 A_new_i) - (A_ip1 A_i)) / dt theta * (x[ni1] - x[ni]) / dx (1 - theta) * (Q_old[i1] - Q_old[i]) / dx) # 对 Z[i]、Z[i1]、Q[i]、Q[i1] 的偏导数填入雅可比 dA_dZ_i sections[i].dA_dZ(Z_old[i]) # 面积对水位的导数由插值表数值求导 J[i, i] 0.5 * dA_dZ_i / dt J[i, i1] 0.5 * dA_dZ_ip1 / dt J[i, ni] - theta / dx J[i, ni1] theta / dx # 动量方程残差第二个方程公式较长此处省略展开 # 关键点对每个网格对补充一组动量残差方程 # 加入边界条件上游流量边界下游水位边界 F[n] x[n] - boundary_up # 上游断面流量 给定流量过程 F[2*n-1] x[n-1] - boundary_down # 下游断面水位 给定水位 # 稀疏矩阵求解增量 dx_new delta spsolve(diags(J), -F) x x delta # 收敛判断残差小于容差 if np.max(np.abs(F)) 1e-4: break Z_new x[:n] Q_new x[n:] return Z_new, Q_new代码逻辑每一步先拿上一时刻的结果做初值组装残差向量和雅可比矩阵然后解线性方程组得到修正量反复迭代到残差收敛。代码里的 dA_dZ 是面积对水位的导数这个值如果直接解析求导很难因为天然断面是不规则的工程上的做法是对插值表做数值差分。注意这里 scipy.sparse 的 diags 只是构造一个占位矩阵真实项目中雅可比是五对角加边界行需要手动构造稀疏矩阵的 row/col/data 数组。雅可比矩阵组装是这段代码的复杂度瓶颈实测 500 个断面时牛顿迭代在 scipy 稀疏求解下大概需要 0.1 到 0.3 秒如果单步超过 1 秒就要检查是不是时间步长设得过大导致迭代不收敛。4. 边界条件、初始场与参数设定糙率 n 值怎么标定才不翻车4.1 上游来流过程、下游水位过程和内部闸堰三类边界条件的处理一维河道模型至少需要上下游各一个边界条件。上游边界通常是流量过程线把水文预报或实测的入流过程按时间序列给定到首断面下游边界通常给水位过程线比如湖泊水位、潮位或者控制站的水位实测。这两类边界条件直接决定了模型在边界处怎么推进代码中就是替换雅可比矩阵中的对应行。如果研究河段内有闸坝、堰、桥梁等水工建筑物就不能简单把河道断面连续处理。闸坝的位置上上下游断面的流量是连续的但水位会出现跃变必须在模型里加内部边界条件。常见做法是在闸位断面布置“流量-水位差”关系用宽顶堰公式或闸孔出流公式计算过闸流量。# 闸坝内部边界条件的处理流量连续性 堰流公式 def gate_boundary(Z_up, Z_down, gate_width, gate_open, coeff0.62): # 自由出流或淹没出流判断 delta_z Z_up - Z_down g 9.81 if delta_z 0: # 自由出流流量取决于上游水头 H max(Z_up - gate_crest_elev, 0.0) Q coeff * gate_width * gate_open * np.sqrt(2 * g * H) else: # 淹没出流流量取决于上下游水位差 Q coeff * gate_width * gate_open * np.sqrt(2 * g * abs(delta_z)) return Q这段代码是一个工程化的简化处理内部边界条件的真实做法是在雅可比矩阵中把闸位前后的两对网格方程替换成流量连续方程和堰流公式并将水位差引入残差。coeff 的取值不是死的宽顶堰通常在 0.85 左右锐缘堰在 0.62 左右闸孔出流在 0.60 到 0.65如果现场有条件建议用实测流量反推。4.2 糙率 n 值的经验取值和率定路径从“查表”到“调参”曼宁糙率 n 是整个一维模型里最“玄学”的参数。它不是一个纯物理量而是综合反映了河道阻力、断面不规则性、植被、蜿蜒度的经验系数。查 Chow 的《Open-Channel Hydraulics》的糙率表是起步做法顺直天然河道 n 在 0.025 到 0.035 之间弯曲且有杂草的河道在 0.035 到 0.05卵石河床可以到 0.04 以上。但这些表值只是起点真实项目里必须用实测水文资料率定。率定的标准流程是先按断面特征初定一个 n 值比如主槽 0.03、滩地 0.05跑一场历史洪水把计算水位和实测水位对比如果计算水位偏高说明糙率偏大降低 n偏低则增加 n。单个断面调好后再看沿程各站点的误差分布。这里有个血泪经验千万别为了拟合某一个站的水位把 n 值调到离谱程度比如城市河道弄出 0.08 的糙率——虽然可能拟合上了但模型就失去外推预测能力了。4.3 初始条件怎么给冷启动和预热策略模型启动时需要各断面初始水位和流量。最简单的做法是给定恒定流初始场假设全河段为恒定均匀流用曼宁公式推算出各断面初始水位。但天然河道特别是下游有回水影响的河段恒定均匀流的假设不成立冷启动会导致前几个时间步内水位剧烈震荡。稳妥的做法是“预热”法用边界条件中起涨前的小流量持续跑一段虚拟时间让模型自动调整到一个自洽的初始状态。预热时长一般取模型时间步长的几百倍比如时间步长 60 秒预热跑 6 小时到 12 小时。判断预热是否完成的标准看全河段流量和水位不再随时间漂移稳定在初始边界条件下的平衡态。跑洪水过程前先做预热这步能省掉后面大量排查“莫名其妙振荡”的时间。5. 一维河道建模的常见坑数值振荡、干河床和回水顶托的排查路径5.1 干湿交替导致负水深现象计算到河段滩地退水阶段水位降到河底以下程序报出负水深或者 NaN。原因天然河道断面在低水位时可能只有主槽过流滩地是干的。但 Preissmann 格式假设断面始终有水离散方程中的面积函数在水位低于河底时返回 0这会导致连续性方程和动量方程之间出现不匹配进而出现负水深。解决在面积插值表中把枯水位以下的面积曲线做平滑延伸不让它直接归零同时在水位低于某个阈值时把该断面的动量方程退化处理忽略压力梯度项。另一种做法是引入“临界深度”判定当某断面水深小于阈值比如 1 厘米就把这个断面的动量方程替换成流量和水位的线性关系保证不出现负水深。5.2 下游回水顶托导致动量方程不收敛现象下游有闸坝或湖泊顶托模型计算到回水末端附近时牛顿迭代不收敛残差在某个断面附近来回震荡。原因回水影响区的流态是缓流物理上信息既可以向上游传播也可以向下游传播而 Preissmann 格式在这种流态下如果 θ 取小了可能出现解的非唯一性。解决把 θ 调高到 0.7 或 0.8 试试再不行就缩减该河段的空间步长。排查方法画出断面流量-水位关系线如果出现明显振荡优先怀疑 θ 和时间步长。回水段的初始场也要格外小心冷启动时给错水位场回水区的调整会非常慢甚至影响整个模拟期的结果。5.3 糙率不分区导致流量分配失真现象漫滩洪水时计算的主槽和滩地流量比例与实测偏差很大水位对不上。原因天然河道主槽和滩地的阻力特性差异悬殊主槽糙率可能只有 0.025滩地因为有植被和农田。从 0.05 到 0.08 不等。如果建模时整个断面用一个糙率值高水位的计算结果必然失真。解决每个断面把过水面积拆成主槽、左滩、右滩三块分别计算面积和湿周各自乘以对应的糙率。代码实现上在断面对象中维护三个糙率值而不是一个。5.4 时间步长过大导致的长周期振荡现象河道较长、步长拉大后模拟曲线出现周期两三个小时的波浪形振荡且不随时间衰减。原因虽然隐式格式“无条件稳定”但物理上可分辨的时间尺度受限于波速和网格长度。时间步长过大模型无法分辨洪水波在相邻断面间的传播过程数值上表现为典型的伪振荡。解决检查 CFL 数即便用隐式格式也要保证 dt 不超过 dx 除以波速的三到五倍。实测中 100 米断面间距、河道波速 3 到 4 米每秒时dt 超过 20 秒就可能有轻微振荡但很多项目用 60 秒也能跑。如果必须用大步长考虑改用工字型断面概化让波速更均匀或者加密局部断面。5.5 断面数据缺失时硬插值的坑现象两条实测断面间距超过 1 公里模型算出来的水面线在中间出现不符合实际的锯齿状跳变。原因河道地形在长距离内变化很大中间可能有卡口、急滩或蜿蜒用线性插值补齐断面间的河床高程和断面形状和实际情况对不上。解决宁可选择水文意义上的控制断面做有限数量断面建模也不要中间硬插。常见做法是把河段按“顺直段-卡口段-弯道段”分区每个区段选一个代表性断面通过控制断面间的能量损失局部水头损失系数来修正。如果你依赖的天然河道数据只有 DEM可以用 GIS 工具沿中心线每隔 200 到 500 米提取剖面线做断面但务必和实际地形图对照DEM 提取的断面在水下地形部分是严重的黑匣子。6. 一个值得养成的收尾习惯用实测水面线做“反向糙率校核”别只看水文站水位模型算完一场洪水水位过程线和流量过程线都对上了是不是就能出报告了还不够。我吃过一次亏某个项目模型率定阶段只看下游水文站的洪峰水位拟合得非常好结果甲方要求对沿程三个临时水尺点位的洪痕做对比算出最大偏差将近 0.8 米。后来查明是上游一段卡口断面概化太平滑导致卡口壅水没体现出来。现在我在每个项目里都会养成一个收尾习惯把实测洪痕哪怕只有三五个点点绘到计算水面线上逐个断面标注偏差。如果某个断面的偏差系统性偏大先别急着调糙率去查这个断面上下游是不是有地形突变——比如桥梁、卡口、弯道——这些位置用圣维南方程组的一维形式根本刻画不了合理的做法是在模型中增加局部水头损失系数来表达而不是用调糙率去掩盖。局部水头损失系数的量级渐变段取 0.1 到 0.2桥墩束窄处取 0.3 到 0.6急弯段取 0.2 左右具体最好用实测压差反推。这个“反向校核”只用一场洪水的数据就能暴露模型中绝大部分的结构性问题比盲目调 n 值更值得投入时间。对预算有限的河段哪怕只有一场小洪水的洪痕也能把模型的可信度从“能跑通”提升到“能用于决策”。一维河道建模归根结底是一个数据工程和参数工程结合的活儿模型结构选对、边界条件给准、糙率分区合理这个方向就值得投入反之如果你连一场完整的水位流量实测都没有建议先补数据再做模型否则算出来的结果只是看起来精致的数字游戏。以上就是我在一维河道水动力学建模这条路上最大的几个教训希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/3 8:31:34 · 阅读 1.2k 一维河道水动力建模：从圣维南方程组到Preissmann格式的Python实践
简介面向水利、水文及计算流体方向学习者这份资源演示了基于 MATLAB 的一维河道水动力学建模流程核心围绕圣维南方程组与 Preissmann 离散法展开可用于理解流速、水位和流量变化模拟以及上、下游边界条件的设定方法。压缩包共5个文件以 .m 主程序、.mat 数据文件和 .txt 说明与地形资料为主整体约401KB结构精简。已有59人学习下载。内容包含完整建模主程序、河道地形数据、输入参数与边界条件示例并配有使用说明与地形资料文本对照代码与说明可掌握一维水动力模型从数据准备、算法离散到结果解读的完整思路对洪水预警、水资源管理和水利工程设计具有直接参考价值。1. 一维河道水动力学建模别一上来就追二维先把圣维南方程组跑稳很多做防洪预报、水环境影响评估的工程师拿到一个河道项目的第一反应是上二维模型。但实际干过几个项目就知道几公里甚至几十公里的天然河道二维网格动辄十万级率定参数调到怀疑人生最后算出来的水位可能还不如一维模型加一套靠谱的糙率。一维河道水动力学建模核心就是沿着河道中心线求解圣维南方程组把流量和水位沿程算出来。它的价值不在于精度碾压二维而是用最小的数据成本和计算成本解决“洪水从哪来、水位到多高、流量怎么演进”这类最刚需的问题。适合刚接触水动力模拟的从业者也适合经费有限、实测断面数据不全的中小型项目。2. 圣维南方程组怎么简化才合理连续方程、动量方程与三种波的取舍2.1 一维河道模型到底在算什么两个方程、两个未知数的闭式系统一维河道水动力模型的控制方程是法国科学家圣维南在1871年提出的浅水方程组由连续方程和动量方程组成。连续方程表达质量守恒河道里某一段水体的水量变化等于上游流进来减去下游流出去还要考虑旁侧入流。动量方程表达牛顿第二定律在水体上的体现惯性力、压力梯度力、重力沿坡降的分量以及河床和岸壁的阻力必须达到平衡。写成常见的守恒形式# 连续方程∂A/∂t ∂Q/∂x q_lateral # 动量方程∂Q/∂t ∂(Q²/A)/∂x g·A·∂Z/∂x g·A·(S0 - Sf) # A: 过水断面面积 (m²) # Q: 流量 (m³/s) # Z: 水位 (m) # q_lateral: 单位长度旁侧入流 (m²/s) # S0: 河床底坡 (m/m) # Sf: 摩阻坡降用曼宁公式计算 Sf n²·Q·|Q| / (A²·R^(4/3))这里有两个未知量——流量 Q 和水位 Z或者断面面积 A两个方程构成闭式系统所以给定边界条件和初始条件理论上就能求出沿程各断面的水位和流量随时间的演变。实际编码时常用水位 Z 和流量 Q 做未知量而不是用断面面积 A因为天然河道断面不规则A 和 Z 的关系是非线性的函数表用 Z 做变量更方便在断面数据表里插值。2.2 运动波、扩散波与动力波不是越完整越好要看你的问题尺度完整圣维南方程组包含所有物理项数值求解难度最大对数据要求也最高。实际工程里会根据问题的物理特性做简化。如果河道底坡陡、洪水波拉长、回水影响小比如山区性河流或城市排水明渠可以把惯性项和压力项都丢掉只保留重力项与摩阻项平衡这就是运动波近似。运动波模型的核心假设是流量与水位之间存在单值关系也就是流量只取决于当地水深和坡降下游的水位根本影响不到上游所以计算简单、速度快适合做长时段洪水演进。如果下游有闸坝、支流汇入、潮汐顶托这类回水影响运动波就失效了因为水位流量关系被迫变成绳套曲线同一水位对应两个流量。这时要保留压力梯度项这就成了扩散波近似。扩散波能够刻画回水影响但依然忽略惯性项适用于坡度适中、回水问题不严重的河道。如果要做溃坝洪水、闸门骤启骤停这类惯性效应显著的场景必须用完整动力波也就是求解完整圣维南方程组但这时的参数率定难度和数值稳定性问题都会陡增。2.3 离散化的数学本质从偏微分方程到差分方程的那一步决定了成败偏微分方程没法直接让计算机求解必须把连续空间和时间离散成网格节点用差分或有限体积逼近导数项。常见做法是空间上沿河道按断面间距划分网格时间上按固定步长推进。这里有三个关键选择空间差分格式、时间积分格式和网格间距。# 圣维南方程组的通用离散框架 # 空间导数三点中心差分或迎风差分 dQ_dx (Q[i1] - Q[i-1]) / (2 * dx) # 中心差分二阶精度 dQ_dx_upwind (Q[i] - Q[i-1]) / dx # 迎风差分适合单向传播 # 时间导数显式或隐式推进 # 显式dQ_dt (Q_new[i] - Q_old[i]) / dt # 隐式dQ_dt (Q_new[i] - Q_old[i]) / dt右端项用 Q_new 取值显式格式实现简单但稳定性受限于 CFL 条件水流在一个时间步内不能跨过超过一个空间网格。也就是说 dt 必须小于 dx 除以波速否则数值解发散。天然河道里波速可能到每秒几米如果断面间距是 100 米时间步长最大也就 20 秒左右。隐式格式没有这个限制时间步长可以放宽到几十倍代价是每步要解一个大型方程组。工程模型普遍选隐式格式原因很现实天然河道断面数据粗糙没必要在时间分辨率上较劲把计算步长放到分钟级一天的水流演进只要几百步就能算完。3. 用 Python 手写一维求解器Preissmann 四点隐式格式与稀疏矩阵组装3.1 断面数据怎么组织河道不是一根线是一串断面动手写代码前先把河道数据整理成稳定、可索引的结构。一维模型中的“河道”在经济意义上是一条中心线但实际水力计算需要的是每个断面桩号处的水位-面积-湿周关系表。野外测量通常给的是断面上若干测点的起点距和高程需要预处理成水深-面积和水深-水力半径的关系函数。# 断面数据处理从散点测量数据到水力要素表 import numpy as np from scipy.interpolate import interp1d # 原始测量点起点距 x 和河底高程 z x_pts np.array([0.0, 5.0, 10.0, 15.0, 20.0]) z_pts np.array([102.5, 101.2, 100.1, 101.0, 102.3]) def cross_section_properties(z_level, x_pts, z_pts): # 计算给定水位下的过水面积 A 和湿周 P # 思路水位以下的测点参与计算水位以上的点用线性插值找水面交点 if z_level np.min(z_pts): return 0.0, 0.0 # 分段线性积分的简化实现 area 0.0 wetted_perimeter 0.0 for i in range(len(x_pts) - 1): z_left z_pts[i] z_right z_pts[i1] if z_left z_level or z_right z_level: # 三角形或梯形面积 h_left max(z_level - z_left, 0.0) h_right max(z_level - z_right, 0.0) width x_pts[i1] - x_pts[i] area 0.5 * (h_left h_right) * width # 湿周只计算淹没部分的底边长 wet_frac max(h_left / (h_left h_right 1e-9), 0.0) wetted_perimeter wet_frac * np.hypot(width, z_right - z_left) return area, wetted_perimeter # 生成水位-面积插值表对每个断面预计算一组离散水位对应的面积 water_levels np.linspace(np.min(z_pts), np.max(z_pts) 2.0, 50) areas [cross_section_properties(zl, x_pts, z_pts)[0] for zl in water_levels] A_interp interp1d(water_levels, areas, fill_valueextrapolate)这段代码的核心思想是先对每个野外实测断面构建一张水位-面积关系表之后求解器在迭代中通过插值函数快速获取任一水位对应的过水面积而不用每次重新积分。参数说明水位数组的上下限决定了模型的适用水位范围上限至少要覆盖设计洪水位以上 1 米否则高水位时插值外延误差不可控。实测断面中如果有明显的滩地或主槽分界最好在 z_pts 里保留地形突变点否则湿周计算会出现跳变。3.2 Preissmann 四点隐式格式的离散模板半个网格处的空间差商Preissmann 格式是一维河道模型中使用最广泛的隐式格式1961 年由 Preissmann 提出特点是空间差商取在两个相邻网格点之间的中间位置时间导数项也取同一位置。这种“半网格”处理让格式在空间上有二阶精度同时保持了隐式格式的无条件稳定性——至少理论上是这样实际中时间步长太大一样会振荡。# Preissmann 四点格式的离散系数模板 # 设变量 f 代表水位 Z 或流量 Q在网格点 i 和 i1 之间取加权平均 # 权重系数 θ控制时间加权工程上通常取 θ 0.6~0.7 theta 0.65 # 空间导数在 n1/2 时刻的离散形式 dZ_dx theta * (Z[i1] - Z[i]) / dx (1 - theta) * (Z_old[i1] - Z_old[i]) / dx # 时间导数在 i1/2 位置的离散形式 dZ_dt 0.5 * ((Z[i1] - Z_old[i1]) / dt (Z[i] - Z_old[i]) / dt) # 连续性方程的最终离散形式 # 0.5 * ((A[i1] A[i]) - (A_old[i1] A_old[i])) / dt # theta * (Q[i1] - Q[i]) / dx # (1 - theta) * (Q_old[i1] - Q_old[i]) / dx # q_lateral (旁侧入流项)θ 参数的物理含义是当前时刻和上一时刻的权重分配。θ 1 时是全隐式虽然最稳定但数值扩散严重会把洪水波拉平θ 0 时退化为显式稳定性条件又回来了。工程经验一般取 0.6 到 0.7折中稳定性和精度。这段代码展示的是单个网格对的离散实际求解时会把所有网格点的离散方程联立成一个矩阵。3.3 组装稀疏矩阵并推进求解从牛顿迭代到稀疏LU分解Preissmann 格式每个网格点会生成两个方程一个来自离散后的连续方程一个来自动量方程。如果有 N 个断面就有 2N 个未知量每个断面的 Z 和 Q对应 2(N-1) 个方程剩余两个方程由上下游边界条件补齐。这个方程组是高度非线性的因为断面面积 A 和水力半径 R 是水位 Z 的强非线性函数所以每个时间步内都要做牛顿迭代。# 完整求解器骨架牛顿迭代 稀疏矩阵求解 from scipy.sparse import diags from scipy.sparse.linalg import spsolve import numpy as np def solve_one_step(Z_old, Q_old, boundary_up, boundary_down, sections, dt, dx, theta0.65): n len(Z_old) # 未知量排列先所有断面的 Z再所有断面的 Q共 2n 个 x np.concatenate([Z_old, Q_old]) for iter in range(20): # 牛顿迭代上限 # 组装残差向量 F 和雅可比矩阵 J F np.zeros(2 * n) J np.zeros((2 * n, 2 * n)) for i in range(n - 1): # 连续方程残差第一个方程 A_i sections[i].area(Z_old[i]) A_ip1 sections[i1].area(Z_old[i1]) A_new_i sections[i].area(x[i]) A_new_ip1 sections[i1].area(x[i1]) F[i] (0.5 * ((A_new_ip1 A_new_i) - (A_ip1 A_i)) / dt theta * (x[ni1] - x[ni]) / dx (1 - theta) * (Q_old[i1] - Q_old[i]) / dx) # 对 Z[i]、Z[i1]、Q[i]、Q[i1] 的偏导数填入雅可比 dA_dZ_i sections[i].dA_dZ(Z_old[i]) # 面积对水位的导数由插值表数值求导 J[i, i] 0.5 * dA_dZ_i / dt J[i, i1] 0.5 * dA_dZ_ip1 / dt J[i, ni] - theta / dx J[i, ni1] theta / dx # 动量方程残差第二个方程公式较长此处省略展开 # 关键点对每个网格对补充一组动量残差方程 # 加入边界条件上游流量边界下游水位边界 F[n] x[n] - boundary_up # 上游断面流量 给定流量过程 F[2*n-1] x[n-1] - boundary_down # 下游断面水位 给定水位 # 稀疏矩阵求解增量 dx_new delta spsolve(diags(J), -F) x x delta # 收敛判断残差小于容差 if np.max(np.abs(F)) 1e-4: break Z_new x[:n] Q_new x[n:] return Z_new, Q_new代码逻辑每一步先拿上一时刻的结果做初值组装残差向量和雅可比矩阵然后解线性方程组得到修正量反复迭代到残差收敛。代码里的 dA_dZ 是面积对水位的导数这个值如果直接解析求导很难因为天然断面是不规则的工程上的做法是对插值表做数值差分。注意这里 scipy.sparse 的 diags 只是构造一个占位矩阵真实项目中雅可比是五对角加边界行需要手动构造稀疏矩阵的 row/col/data 数组。雅可比矩阵组装是这段代码的复杂度瓶颈实测 500 个断面时牛顿迭代在 scipy 稀疏求解下大概需要 0.1 到 0.3 秒如果单步超过 1 秒就要检查是不是时间步长设得过大导致迭代不收敛。4. 边界条件、初始场与参数设定糙率 n 值怎么标定才不翻车4.1 上游来流过程、下游水位过程和内部闸堰三类边界条件的处理一维河道模型至少需要上下游各一个边界条件。上游边界通常是流量过程线把水文预报或实测的入流过程按时间序列给定到首断面下游边界通常给水位过程线比如湖泊水位、潮位或者控制站的水位实测。这两类边界条件直接决定了模型在边界处怎么推进代码中就是替换雅可比矩阵中的对应行。如果研究河段内有闸坝、堰、桥梁等水工建筑物就不能简单把河道断面连续处理。闸坝的位置上上下游断面的流量是连续的但水位会出现跃变必须在模型里加内部边界条件。常见做法是在闸位断面布置“流量-水位差”关系用宽顶堰公式或闸孔出流公式计算过闸流量。# 闸坝内部边界条件的处理流量连续性 堰流公式 def gate_boundary(Z_up, Z_down, gate_width, gate_open, coeff0.62): # 自由出流或淹没出流判断 delta_z Z_up - Z_down g 9.81 if delta_z 0: # 自由出流流量取决于上游水头 H max(Z_up - gate_crest_elev, 0.0) Q coeff * gate_width * gate_open * np.sqrt(2 * g * H) else: # 淹没出流流量取决于上下游水位差 Q coeff * gate_width * gate_open * np.sqrt(2 * g * abs(delta_z)) return Q这段代码是一个工程化的简化处理内部边界条件的真实做法是在雅可比矩阵中把闸位前后的两对网格方程替换成流量连续方程和堰流公式并将水位差引入残差。coeff 的取值不是死的宽顶堰通常在 0.85 左右锐缘堰在 0.62 左右闸孔出流在 0.60 到 0.65如果现场有条件建议用实测流量反推。4.2 糙率 n 值的经验取值和率定路径从“查表”到“调参”曼宁糙率 n 是整个一维模型里最“玄学”的参数。它不是一个纯物理量而是综合反映了河道阻力、断面不规则性、植被、蜿蜒度的经验系数。查 Chow 的《Open-Channel Hydraulics》的糙率表是起步做法顺直天然河道 n 在 0.025 到 0.035 之间弯曲且有杂草的河道在 0.035 到 0.05卵石河床可以到 0.04 以上。但这些表值只是起点真实项目里必须用实测水文资料率定。率定的标准流程是先按断面特征初定一个 n 值比如主槽 0.03、滩地 0.05跑一场历史洪水把计算水位和实测水位对比如果计算水位偏高说明糙率偏大降低 n偏低则增加 n。单个断面调好后再看沿程各站点的误差分布。这里有个血泪经验千万别为了拟合某一个站的水位把 n 值调到离谱程度比如城市河道弄出 0.08 的糙率——虽然可能拟合上了但模型就失去外推预测能力了。4.3 初始条件怎么给冷启动和预热策略模型启动时需要各断面初始水位和流量。最简单的做法是给定恒定流初始场假设全河段为恒定均匀流用曼宁公式推算出各断面初始水位。但天然河道特别是下游有回水影响的河段恒定均匀流的假设不成立冷启动会导致前几个时间步内水位剧烈震荡。稳妥的做法是“预热”法用边界条件中起涨前的小流量持续跑一段虚拟时间让模型自动调整到一个自洽的初始状态。预热时长一般取模型时间步长的几百倍比如时间步长 60 秒预热跑 6 小时到 12 小时。判断预热是否完成的标准看全河段流量和水位不再随时间漂移稳定在初始边界条件下的平衡态。跑洪水过程前先做预热这步能省掉后面大量排查“莫名其妙振荡”的时间。5. 一维河道建模的常见坑数值振荡、干河床和回水顶托的排查路径5.1 干湿交替导致负水深现象计算到河段滩地退水阶段水位降到河底以下程序报出负水深或者 NaN。原因天然河道断面在低水位时可能只有主槽过流滩地是干的。但 Preissmann 格式假设断面始终有水离散方程中的面积函数在水位低于河底时返回 0这会导致连续性方程和动量方程之间出现不匹配进而出现负水深。解决在面积插值表中把枯水位以下的面积曲线做平滑延伸不让它直接归零同时在水位低于某个阈值时把该断面的动量方程退化处理忽略压力梯度项。另一种做法是引入“临界深度”判定当某断面水深小于阈值比如 1 厘米就把这个断面的动量方程替换成流量和水位的线性关系保证不出现负水深。5.2 下游回水顶托导致动量方程不收敛现象下游有闸坝或湖泊顶托模型计算到回水末端附近时牛顿迭代不收敛残差在某个断面附近来回震荡。原因回水影响区的流态是缓流物理上信息既可以向上游传播也可以向下游传播而 Preissmann 格式在这种流态下如果 θ 取小了可能出现解的非唯一性。解决把 θ 调高到 0.7 或 0.8 试试再不行就缩减该河段的空间步长。排查方法画出断面流量-水位关系线如果出现明显振荡优先怀疑 θ 和时间步长。回水段的初始场也要格外小心冷启动时给错水位场回水区的调整会非常慢甚至影响整个模拟期的结果。5.3 糙率不分区导致流量分配失真现象漫滩洪水时计算的主槽和滩地流量比例与实测偏差很大水位对不上。原因天然河道主槽和滩地的阻力特性差异悬殊主槽糙率可能只有 0.025滩地因为有植被和农田。从 0.05 到 0.08 不等。如果建模时整个断面用一个糙率值高水位的计算结果必然失真。解决每个断面把过水面积拆成主槽、左滩、右滩三块分别计算面积和湿周各自乘以对应的糙率。代码实现上在断面对象中维护三个糙率值而不是一个。5.4 时间步长过大导致的长周期振荡现象河道较长、步长拉大后模拟曲线出现周期两三个小时的波浪形振荡且不随时间衰减。原因虽然隐式格式“无条件稳定”但物理上可分辨的时间尺度受限于波速和网格长度。时间步长过大模型无法分辨洪水波在相邻断面间的传播过程数值上表现为典型的伪振荡。解决检查 CFL 数即便用隐式格式也要保证 dt 不超过 dx 除以波速的三到五倍。实测中 100 米断面间距、河道波速 3 到 4 米每秒时dt 超过 20 秒就可能有轻微振荡但很多项目用 60 秒也能跑。如果必须用大步长考虑改用工字型断面概化让波速更均匀或者加密局部断面。5.5 断面数据缺失时硬插值的坑现象两条实测断面间距超过 1 公里模型算出来的水面线在中间出现不符合实际的锯齿状跳变。原因河道地形在长距离内变化很大中间可能有卡口、急滩或蜿蜒用线性插值补齐断面间的河床高程和断面形状和实际情况对不上。解决宁可选择水文意义上的控制断面做有限数量断面建模也不要中间硬插。常见做法是把河段按“顺直段-卡口段-弯道段”分区每个区段选一个代表性断面通过控制断面间的能量损失局部水头损失系数来修正。如果你依赖的天然河道数据只有 DEM可以用 GIS 工具沿中心线每隔 200 到 500 米提取剖面线做断面但务必和实际地形图对照DEM 提取的断面在水下地形部分是严重的黑匣子。6. 一个值得养成的收尾习惯用实测水面线做“反向糙率校核”别只看水文站水位模型算完一场洪水水位过程线和流量过程线都对上了是不是就能出报告了还不够。我吃过一次亏某个项目模型率定阶段只看下游水文站的洪峰水位拟合得非常好结果甲方要求对沿程三个临时水尺点位的洪痕做对比算出最大偏差将近 0.8 米。后来查明是上游一段卡口断面概化太平滑导致卡口壅水没体现出来。现在我在每个项目里都会养成一个收尾习惯把实测洪痕哪怕只有三五个点点绘到计算水面线上逐个断面标注偏差。如果某个断面的偏差系统性偏大先别急着调糙率去查这个断面上下游是不是有地形突变——比如桥梁、卡口、弯道——这些位置用圣维南方程组的一维形式根本刻画不了合理的做法是在模型中增加局部水头损失系数来表达而不是用调糙率去掩盖。局部水头损失系数的量级渐变段取 0.1 到 0.2桥墩束窄处取 0.3 到 0.6急弯段取 0.2 左右具体最好用实测压差反推。这个“反向校核”只用一场洪水的数据就能暴露模型中绝大部分的结构性问题比盲目调 n 值更值得投入时间。对预算有限的河段哪怕只有一场小洪水的洪痕也能把模型的可信度从“能跑通”提升到“能用于决策”。一维河道建模归根结底是一个数据工程和参数工程结合的活儿模型结构选对、边界条件给准、糙率分区合理这个方向就值得投入反之如果你连一场完整的水位流量实测都没有建议先补数据再做模型否则算出来的结果只是看起来精致的数字游戏。以上就是我在一维河道水动力学建模这条路上最大的几个教训希望帮到你。本文还有配套的精品资源点击获取