一维河道水动力模型：圣维南方程离散与MATLAB实现

简介：这份资源是一套面向水文学、水利工程方向研究人员与高年级学生的Matlab一维水动力模型代码，用于模拟矩形河道内流速与水深变化，可服务于河流管理、洪水预警与水利工程设计等场景。模型基于圣维南方程组，涵盖连续性方程与动量方程，并考虑入口、出口、自由表面等边界条件及糙率系数对摩擦损失的影响，通过数值方法求解河道水位与流速分布。压缩包共4个文件，包含2个m脚本与2个xlsx数据表，前者承担初始化、主程序求解与结果输出，后者用于输入参数与边界条件配置，整体约323KB，结构精简便于二次修改。目前已有1603人学习下载，适合具备Matlab编程基础与水动力学知识的读者参考，可据此理解一维流动的离散求解思路、时间步进方案与结果可视化流程，为河道水流特性分析与水环境管理提供可复用的计算工具。

1. 一维河道水动力模型：从圣维南方程到可跑的 MATLAB 代码

做洪水演进、灌区输配水或者河道整治方案的人，迟早会碰到同一个问题：手上只有断面地形和上下游流量过程，怎么算出沿程水位随时间的变化。商业软件能算，但授权、二次开发、批量工况都卡脖子。一维河道水动力模型就是干这件事的——把圣维南方程组在河道断面上离散，用数值方法推水位和流量。MATLAB 在这里不是玩具，它的矩阵运算、稀疏求解器和绘图能力，足够撑起一条几百个断面、几天模拟期的河道计算。这篇文章面向已经懂一点水力学、但还没自己写过求解器的工程师，从方程离散讲到代码结构，再到参数怎么调、结果怎么验。读完你应该能自己搭一个能跑通、能出图、能改边界的一维模型，而不是只会点商业软件的按钮。

2. 圣维南方程怎么离散才能不炸：Preissmann 四点隐式格式拆解

2.1 为什么显式格式在一维河道里容易翻车

一维河道水动力模型的核心是圣维南方程组，连续方程和动量方程各一条。显式格式听起来简单，时间步长受柯朗条件限制，河道水深变浅或者断面间距拉大时，时间步长会被压到很小，算一天洪水可能要跑几个小时。更麻烦的是，显式格式对急变流和干湿边界非常敏感，稍微有点底坡突变就发散。我早期用显式格式做过一个山区河道，上游涨水快，下游断面间距 200 米，时间步长被迫降到 0.5 秒，算 6 小时洪水花了 40 分钟，而且中间还炸了两次。

隐式格式把时间步长从柯朗条件里解放出来，代价是每个时间步要解一个大型线性方程组。Preissmann 四点隐式格式是河道计算里最常用的，因为它对时间层和空间层都做加权平均，权重系数可以调，既能保证稳定性，又能控制数值耗散。MATLAB 的稀疏矩阵求解器\对带状矩阵效率很高，几百个断面的方程组求解时间可以忽略不计。

2.2 Preissmann 格式的离散形式和系数含义

Preissmann 格式对任意变量 f 在空间点 j 和时间层 n 的加权平均写成：

f ≈ θ(f_j^{n+1} + f_{j+1}^{n+1})/2 + (1-θ)(f_j^n + f_{j+1}^n)/2

时间导数写成：

∂f/∂t ≈ (f_j^{n+1} + f_{j+1}^{n+1} - f_j^n - f_{j+1}^n) / (2Δt)

空间导数写成：

∂f/∂x ≈ θ(f_{j+1}^{n+1} - f_j^{n+1})/Δx + (1-θ)(f_{j+1}^n - f_j^n)/Δx

θ 取 0.5 到 1.0 之间，θ=1 是全隐式，稳定性最好但精度一阶；θ=0.5 是二阶精度但可能产生振荡。实际河道计算里 θ 取 0.6 到 0.75 比较稳。Δt 和 Δx 不需要满足柯朗条件，但 Δx 太大时河道地形细节被抹平，一般建议 Δx 在 50 到 500 米之间，视河道坡度和断面变化剧烈程度定。

把连续方程和动量方程按上述格式展开，每个断面产生两个未知数：水位 Z 和流量 Q。相邻断面之间的方程线性化后，形成一组三对角或者五对角的稀疏线性方程组。边界条件通常给上游流量过程 Q(t) 和下游水位过程 Z(t)，或者下游用水位-流量关系曲线。

2.3 用 MATLAB 组装系数矩阵的最小代码框架

下面这段代码展示单个时间步内如何组装系数矩阵。假设有 N 个断面，未知数排列为 [Z1, Q1, Z2, Q2, ..., ZN, QN]，共 2N 个未知数。连续方程和动量方程各提供 N-1 个方程，加上上游和下游各一个边界方程，刚好 2N 个。

% 一维河道 Preissmann 格式系数矩阵组装示例 % 输入：N 断面数，dt 时间步长，dx 断面间距向量，B 水面宽，A 过水面积， % Q_old, Z_old 上一时间层流量和水位，theta 权重系数 % 输出：稀疏矩阵 M 和右端项 rhs N = length(Z_old); M = spalloc(2*N, 2*N, 8*N); % 预分配稀疏矩阵 rhs = zeros(2*N, 1); for j = 1:N-1 % 连续方程：对断面 j 和 j+1 的 Z 和 Q 做加权平均 % 方程编号 2*j-1 row = 2*j - 1; dx_j = dx(j); % 对 Z_j^{n+1} 的系数 M(row, 2*j-1) = -theta * B(j) / (2*dt) - theta / dx_j; % 对 Z_{j+1}^{n+1} 的系数 M(row, 2*j+1) = -theta * B(j+1) / (2*dt) + theta / dx_j; % 对 Q_j^{n+1} 的系数 M(row, 2*j) = theta / (2*dx_j); % 对 Q_{j+1}^{n+1} 的系数 M(row, 2*j+2) = theta / (2*dx_j); % 右端项来自上一时间层 rhs(row) = (1-theta) * (Q_old(j+1) - Q_old(j)) / dx_j ... - (1-theta) * (B(j+1)*Z_old(j+1) - B(j)*Z_old(j)) / (2*dt) ... - (B(j+1)*Z_old(j+1) - B(j)*Z_old(j)) / (2*dt); % 动量方程：方程编号 2*j row = 2*j; % 动量方程系数涉及摩阻项、对流项和压力项 % 这里只给出结构示意，实际系数需要根据断面几何和流速计算 M(row, 2*j-1) = -theta * 9.81 * A(j) / dx_j; M(row, 2*j+1) = theta * 9.81 * A(j+1) / dx_j; M(row, 2*j) = theta / (2*dt) + theta * abs(Q_old(j)) * 2 * 0.025^2 / (A(j)*1.0); M(row, 2*j+2) = theta / (2*dt); rhs(row) = (1-theta) * (Q_old(j) - Q_old(j+1)) / dx_j ... + (1-theta) * 9.81 * (A(j+1) - A(j)) / dx_j ... + (Q_old(j) - Q_old(j+1)) / (2*dt); end % 上游边界：给定流量 Q1 = Q_up M(2*N-1, 2) = 1.0; rhs(2*N-1) = Q_up; % 下游边界：给定水位 ZN = Z_down M(2*N, 2*N-1) = 1.0; rhs(2*N) = Z_down; % 求解 solution = M \ rhs; Z_new = solution(1:2:end); Q_new = solution(2:2:end);

这段代码的关键点在于稀疏矩阵的预分配和索引映射。spalloc提前分配非零元空间，避免循环中动态扩容。未知数排列采用 Z 和 Q 交替，使得每个断面的两个未知数在矩阵中相邻，带宽控制在 3 以内。动量方程里的摩阻项用了曼宁公式的线性化形式，糙率系数 0.025 是示例值，实际要根据河道材质调整。边界方程直接替换最后两行，上游给流量、下游给水位是最常见的组合。

提示：动量方程的对流项线性化容易出错，建议先用小流速工况验证，确认矩阵条件数在 1e8 以下再放大流量。

3. 断面几何和糙率怎么设：从实测大断面到水力半径的映射

3.1 实测断面数据的预处理和插值

实际河道测量给的是离散点坐标 (x, z)，x 是起点距，z 是高程。计算过水面积和水面宽时，需要根据当前水位在断面上插值。MATLAB 的interp1可以快速完成，但要注意断面点排序和重复点处理。我一般先把断面点按起点距排序，去掉重复的 x 值，然后对水位做二分查找确定淹没范围。

% 断面几何计算：给定水位 Z，返回过水面积 A、水面宽 B、湿周 P function [A, B, P] = section_geometry(x_sec, z_sec, Z) % x_sec: 起点距向量，z_sec: 高程向量，Z: 当前水位 % 确保按起点距排序 [x_sec, idx] = sort(x_sec); z_sec = z_sec(idx); % 找到水位以下的区间 below = z_sec <= Z; if ~any(below) A = 0; B = 0; P = 0; return; end % 在断面两端插值出水位与地形的交点 x_left = interp1(z_sec, x_sec, Z, 'linear', 'extrap'); x_right = x_left; % 需要根据实际断面形状调整 % 简化处理：对淹没区域做数值积分 x_fine = linspace(min(x_sec), max(x_sec), 500); z_fine = interp1(x_sec, z_sec, x_fine, 'linear'); depth = max(Z - z_fine, 0); A = trapz(x_fine, depth); B = sum(depth > 0) * (x_fine(2) - x_fine(1)); P = sum(sqrt(diff(x_fine).^2 + diff(z_fine).^2) .* (depth(1:end-1) > 0)); end

这个函数每次调用都做一次 500 点插值，对于 1000 个断面、1000 个时间步的模拟，计算量在可接受范围。如果嫌慢，可以预先把断面离散成固定网格，用查表法替代实时插值。水面宽 B 的计算用了淹没网格数乘以网格间距，精度取决于网格密度，500 点对于一般河道够用。

3.2 糙率参数的空间分布和率定策略

曼宁糙率 n 是一维河道模型里最玄学的参数。同一个断面，枯水期和洪水期的等效糙率能差一倍。常见做法是分段给常数，比如主槽 0.025，滩地 0.045，然后根据实测水位流量关系反推。率定时先固定其他参数，只调糙率，用下游水位过程的纳什效率系数做目标函数。

率定时的血泪经验是：不要一次性调所有断面的糙率。先按河道分段，每段给一个初始值，然后从下游往上游逐段调。下游水位拟合好了再调上游，因为下游边界误差会向上游传播。如果上下游同时调，很容易陷入局部最优，调半天还不如初始值。

3.3 干湿边界的处理技巧

河道计算里最头疼的是干湿边界。水位低于河底时，过水面积为零，方程退化。直接令 A=0 会导致矩阵奇异。常见做法是给一个最小水深，比如 0.01 米，对应的面积和宽按三角形或梯形断面折算。当水深小于这个阈值时，强制流量为零，只解水位方程。

% 干湿边界处理：最小水深阈值 h_min = 0.01; % 最小水深，单位米 for j = 1:N h = Z_new(j) - z_bed(j); if h < h_min Z_new(j) = z_bed(j) + h_min; Q_new(j) = 0; end end

这个处理简单粗暴但有效。缺点是质量不守恒，因为强制抬高水位相当于凭空加水。改进方法是用窄缝法，把干河床当成一个很窄的矩形槽，宽度取 0.001 米，这样面积连续变化，矩阵不会奇异。窄缝法的代价是时间步长不能太大，否则窄缝里的水位变化会引入数值振荡。

4. 时间步长和空间步长怎么配：稳定性与精度的平衡

4.1 隐式格式下的步长选择经验

隐式格式理论上无条件稳定，但实际计算里 Δt 和 Δx 的比值影响精度和收敛速度。经验公式是 Δt ≤ Δx / (2sqrt(gh_max))，这个条件比显式的柯朗条件宽松 2 倍左右。对于水深 5 米、Δx=200 米的河道，Δt 可以取到 20 秒左右。如果取 60 秒，相位误差会明显增大，洪峰到达时间可能偏差半小时以上。

我一般先用 Δt=10 秒跑一个短时段，看水位过程是否光滑。如果出现锯齿状振荡，要么减小 Δt，要么增大 θ。θ 从 0.6 调到 0.8 能压住大部分振荡，代价是洪峰坦化更严重。对于洪水预报，洪峰坦化 5% 以内可以接受；对于水位流量关系率定，坦化超过 2% 就要警惕。

4.2 空间步长的非均匀划分

河道断面间距往往不均匀，弯道和卡口处需要加密。Preissmann 格式对非均匀网格的处理很直接，每个断面对用各自的 Δx。但相邻断面的 Δx 比值不要超过 2，否则截断误差会集中在突变处。如果实测断面间距从 50 米跳到 500 米，中间要插值补断面。

% 非均匀断面间距检查和平滑 dx = diff(x_sec); ratio = dx(2:end) ./ dx(1:end-1); bad = find(ratio > 2 | ratio < 0.5); if ~isempty(bad) warning('断面间距突变位置：%s，建议插值加密', mat2str(bad)); end

插值加密用线性插值就够，因为一维模型本身对地形细节不敏感。加密后的断面几何按距离加权平均，糙率取相邻断面的较大值，这样偏安全。

4.3 收敛判据和迭代控制

每个时间步内，Preissmann 格式的系数矩阵依赖上一时间层的流速和水深，严格来说需要迭代。但实际计算里，如果 Δt 不大，单步线性化（不迭代）的误差可以接受。判断标准是：连续方程残差和动量方程残差都小于 1e-4。如果残差大，要么减小 Δt，要么加一步迭代。

% 单步残差检查 residual_cont = norm(M(1:2:end-2,:) * solution - rhs(1:2:end-2)); residual_mom = norm(M(2:2:end-1,:) * solution - rhs(2:2:end-1)); if residual_cont > 1e-4 || residual_mom > 1e-4 fprintf('时间步 %d 残差偏大：连续 %.2e，动量 %.2e\n', ... t_step, residual_cont, residual_mom); end

残差突然增大通常意味着干湿边界切换或者边界条件跳变。这时候不要硬算，回退到上一时间步，把 Δt 减半重试。我一般设一个最小 Δt，比如 1 秒，如果减到 1 秒还不行，就检查边界输入有没有异常值。

5. 避坑与排查：一维河道模型跑不通的五个典型场景

5.1 矩阵奇异：水位低于河底导致面积为零

现象：求解时 MATLAB 报「矩阵接近奇异或缩放错误」，或者解出来的水位是 NaN。

原因：某个断面水位低于河底，过水面积 A=0，动量方程里的 1/A 项爆炸。或者干湿边界处理时强制 Q=0 但没同步修改矩阵系数。

解决：在组装矩阵前检查每个断面的水深，小于 h_min 时用窄缝法替代。窄缝宽度取 0.001 到 0.01 米，面积按矩形计算。同时把该断面的动量方程替换成 Q=0 的强制方程，避免除以零。

5.2 水位振荡：θ 太小或 Δt 太大

现象：水位过程线出现高频锯齿，相邻时间步水位差超过 0.5 米。

原因：Preissmann 格式在 θ<0.5 时可能不稳定，θ=0.5 时数值耗散最小但容易激发高频振荡。Δt 超过 Δx/(2sqrt(gh)) 时，相位误差累积也会导致振荡。

解决：先把 θ 调到 0.7 到 0.8，如果振荡消失，说明是格式问题。如果还有振荡，把 Δt 减半。如果减半后振荡频率变了但幅度没降，检查边界条件是否平滑，上游流量过程有没有阶跃。

5.3 质量不守恒：干湿边界强制抬高水位

现象：模拟结束后，上游来水总量和下游出流总量差超过 5%。

原因：窄缝法或最小水深法在干河床上强制抬高水位，相当于每个时间步都往河道里加水。干湿交替频繁的河段，累积误差很大。

解决：记录每个时间步的强制加水量，在总水量平衡里扣除。如果加水量超过来水量的 2%，说明最小水深阈值设大了，从 0.01 降到 0.005。或者改用真正的干湿边界算法，比如把干单元从方程组里剔除，只解湿单元。

5.4 洪峰坦化严重：数值耗散过大

现象：模拟的洪峰流量比上游边界小 15% 以上，洪峰到达时间滞后超过 1 小时。

原因：θ 太大（接近 1.0）导致一阶精度，数值耗散强。Δx 太大，河道调蓄作用被高估。糙率偏大，阻力耗散能量。

解决：θ 降到 0.6，Δx 加密到 100 米以内，糙率按实测水位流量关系重新率定。如果还不行，检查动量方程里的对流项有没有被错误忽略。一维河道里对流项虽然小，但洪峰陡峭时不能省。

5.5 边界反射：下游边界处理不当

现象：下游水位过程出现虚假波动，或者上游流量过程被反射回来。

原因：下游给水位边界时，如果水位过程有高频噪声，会向上游传播。或者下游边界离感兴趣区域太近，反射波还没耗散完就到达计算域。

解决：下游边界至少放在感兴趣区域下游 5 到 10 公里。水位边界先做低通滤波，截止频率取 1/3600 Hz。如果下游是水位流量关系，用宽矩形断面近似，避免关系曲线斜率突变。

6. 用实测水位流量关系反推糙率：一个可复现的率定脚本

率定糙率的核心思路是：给定实测水位过程，调整各河段糙率，使模拟水位和实测水位的均方根误差最小。下面这个脚本用 fminsearch 做单参数优化，适合先率定一个河段的平均糙率，再手动微调。

% 糙率率定：单河段平均糙率优化 % 假设已有函数 run_model(n_manning) 返回模拟水位矩阵 Z_sim % Z_obs 是实测水位矩阵，同维度 Z_obs = load('observed_stage.mat'); % 实测水位，单位米 n0 = 0.030; % 初始糙率 % 目标函数：均方根误差 obj_fun = @(n) sqrt(mean((run_model(n) - Z_obs).^2, 'all')); % 优化范围限制在 0.015 到 0.060 lb = 0.015; ub = 0.060; options = optimset('Display', 'iter', 'TolX', 1e-4); [n_opt, fval] = fminsearch(@(n) obj_fun(max(min(n, ub), lb)), n0, options); fprintf('最优糙率：%.4f，均方根误差：%.4f 米\n', n_opt, fval);

这个脚本的关键是run_model函数要封装完整的时间步进循环，输入糙率返回水位矩阵。优化时每次迭代跑一遍完整模型，如果模型跑一次要 10 秒，优化 20 次就是 200 秒，可以接受。如果河段多，把糙率向量作为优化变量，用fmincon加约束，避免糙率出现负值。

率定完的验证方法是：用另一场洪水做检验，如果均方根误差在 0.1 米以内，说明糙率可靠。如果检验误差比率定误差大 50% 以上，说明糙率是「过拟合」了，需要减少率定参数，比如把三个河段合并成两个。

我自己的习惯是：率定前先把实测水位过程画出来，看看有没有明显异常点。有一次率定结果总是不好，后来发现是下游水位计在洪峰期间被冲歪了，数据本身有问题。模型跑不通的时候，先怀疑数据，再怀疑代码。希望帮到你。

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