简介：这份资源面向水文学、水利工程方向的研究人员与工程技术人员，提供一套基于Matlab实现的一维河道水动力计算代码，用于模拟矩形河道内流速与水深变化，可服务于河流管理、洪水预警与水利工程设计等场景。模型以圣维南方程组为核心，涵盖连续性方程与动量方程，并考虑边界条件、糙率系数及非恒定流入流出等因素，通过数值方法求解水位与流速分布。压缩包共4个文件，包含2个m脚本与2个xlsx数据表，前者承担初始化、主程序求解与结果输出等计算逻辑，后者用于存放输入参数与算例数据，整体约323KB，结构精简便于二次修改。目前已有1603人学习下载，适合具备Matlab编程基础与水动力学知识的读者参考，可据此理解一维水动力建模思路、数值离散方案与结果可视化流程，并在此基础上开展河道水流特性分析与工程应用验证。
1. 一维河道水动力模型：这套 MATLAB 代码到底能算什么
如果你手头有一条几十公里的天然河道，上游有实测流量过程，下游有水位控制，领导要你三天内给出沿程各断面的水位变化——商业软件授权还没批下来，这时候一套能跑的一维河道水动力模型 MATLAB 代码就是救命的。它干的事很明确：把 Saint-Venant 方程组在河道断面上离散，用数值方法推演水位和流量随时间的演化，输出每个断面、每个时刻的 Z 和 Q。适合水文水资源专业的研究生做论文复现，也适合设计院工程师做洪水演进快速估算。但别指望它替代 MIKE 11 或 HEC-RAS 做正式报告，它的定位是「可控、可改、可拆」——你能看到每一行离散格式，能改边界条件，能换求解器。这套代码的核心价值在于透明，不是精度碾压。
2. 从 Saint-Venant 方程到离散格式：代码里到底写了什么
2.1 控制方程与 Preissmann 四点隐格式的选型理由
一维河道水动力模型的数学基础是明渠非恒定流的 Saint-Venant 方程组，包括连续方程和动量方程。连续方程描述水量守恒，动量方程描述能量和动量平衡。写成守恒形式后，对任意断面之间的控制体做积分，就得到有限差分或有限体积的离散基础。
这套代码大概率用的是 Preissmann 四点隐格式，这是一维河道计算里最经典的方案。为什么选隐式而不是显式？显式格式受 CFL 条件限制，时间步长必须很小，对于长河道、慢变洪水波，算到天荒地老。隐式格式理论上无条件稳定，时间步长可以取到 10 分钟甚至 1 小时，计算效率高出一个量级。代价是每个时间步要解一个大型稀疏线性方程组，但 MATLAB 的稀疏矩阵求解器足够应付几百到几千个断面的规模。
Preissmann 格式的核心是在时空网格的四个点上对变量和导数做加权平均。权重系数 θ 一般取 0.6 到 1.0，取 1.0 就是全隐式，稳定性最好但时间精度一阶；取 0.5 是 Crank-Nicolson 格式，时间二阶精度但可能产生数值振荡。代码里通常默认 θ=0.6 或 0.7，这是精度和稳定性的折中。
空间离散上，河道被切成 N 个断面，相邻断面之间是一个计算河段。每个河段上对连续方程和动量方程分别离散，得到两个代数方程。N 个断面产生 2N 个未知量（每个断面的 Z 和 Q），N-1 个河段提供 2(N-1) 个方程，加上上下游各一个边界条件，刚好 2N 个方程闭合。
2.2 核心离散代码拆解与参数含义
下面这段代码展示的是单个河段上连续方程的 Preissmann 离散实现，这是整个模型最核心的一步：
% Preissmann四点隐格式 - 连续方程离散 % 输入: 河段索引 j, 时间层 n, 权重系数 theta % 输出: 系数矩阵 A 和右端项 b 的贡献 function [A_cont, b_cont] = continuity_eq(j, n, theta, dt, dx, B, Q, Z, ... Q_prev, Z_prev, Q_next, Z_next) % 四点加权系数 % theta 控制时间离散: 0.5=Crank-Nicolson, 1.0=全隐式 % 空间权重统一取 0.5 (中心差分) % 连续方程: dZ/dt + (1/B)*dQ/dx = 0 % 离散后: (Z(j,n+1)-Z(j,n))/dt + theta*(...) % 时间导数项系数 coef_Zj = 1/dt; % 当前断面 j 的水位系数 coef_Zj1 = 1/dt; % 下游断面 j+1 的水位系数 % 空间导数项系数 (对 Q 做中心差分) coef_Qj = theta/(2*dx*B(j)); % 断面 j 的流量系数 coef_Qj1 = theta/(2*dx*B(j+1)); % 断面 j+1 的流量系数 % 组装到全局矩阵 (此处为局部贡献示意) A_cont = [coef_Zj, coef_Qj, coef_Zj1, coef_Qj1]; % 右端项: 包含上一时间层的已知量 b_cont = (Z_prev(j) + Z_prev(j+1))/dt ... - (1-theta)*(Q_next(j) - Q_next(j+1))/(2*dx*B(j)); end这段代码的逻辑是：把连续方程中的时间导数用前向差分近似，空间导数用中心差分近似，然后对时间层做 θ 加权。theta是关键参数，取 0.5 时时间精度二阶但可能振荡，取 1.0 时稳定但耗散大。dt是时间步长，单位秒，一般取 60 到 3600。dx是河段长度，单位米，取决于断面间距。B是断面水面宽，单位米，随水位变化，所以是非线性项，需要迭代求解。
动量方程的离散类似，但多了摩阻项和压力项。摩阻坡度通常用 Manning 公式计算，糙率 n 是率定参数，天然河道一般在 0.025 到 0.045 之间。压力项涉及水位对空间的偏导，离散后耦合到相邻断面的水位差。
组装完全部河段的方程后，得到一个块三对角形式的稀疏矩阵。MATLAB 里用sparse存储，用\求解。对于 1000 个断面，矩阵规模 2000×2000，稀疏求解耗时在毫秒级，整个洪水过程的计算可以在几分钟内完成。
2.3 边界条件的代码实现与选择
边界条件是水动力模型的命门。上游通常给流量过程线，下游给水位-流量关系或水位过程线。代码里一般用两种方式实现：一是直接把已知值代入方程，二是用虚拟河段法。
% 上游边界: 给定流量过程线 Q(t) % 方式: 将 Q(1,n+1) 设为已知值, 从方程组中消去该未知量 Q_upstream = interp1(t_series, Q_series, (n+1)*dt, 'linear', 'extrap'); % 修改系数矩阵: 第1个方程替换为 Q(1) = Q_upstream A(1,:) = 0; A(1,1) = 1; b(1) = Q_upstream; % 下游边界: 给定水位-流量关系 Z = f(Q) % 常见做法: 用曼宁公式反算, 或直接给水位过程线 Z_downstream = interp1(t_series, Z_series, (n+1)*dt, 'linear', 'extrap'); % 替换最后一个方程 A(end,:) = 0; A(end,end) = 1; b(end) = Z_downstream;上游给流量是最常见的，因为水文站实测流量比水位可靠。下游给水位关系要注意：如果河道坡度很缓，下游水位控制会向上游传播很远，边界误差会影响整个计算域。我一般建议下游边界至少放在实测水位站位置，不要凭空外推。
提示：边界条件的时间序列要和模型时间步长对齐，用
interp1插值时注意外推方式，'extrap'可能在序列末尾产生异常值，建议手动检查首尾。
3. 从零跑通：断面数据准备、参数率定与结果验证
3.1 断面地形数据的组织与导入
一维河道模型需要两类几何数据：断面间距和每个断面的高程-宽度关系。断面间距决定了空间步长 dx，一般从河道地形图上量取，或者用实测断面桩号。高程-宽度关系描述了断面在不同水位下的过水面积和水面宽，这是模型计算的基础。
常见的数据格式是每个断面一个文本文件，第一列是高程，第二列是宽度，第三列是面积。代码里通常写一个读取函数，把全部断面读进来存成结构体数组：
% 读取断面地形数据 % 文件格式: 每行 "高程(m) 宽度(m) 面积(m2)" % 每个断面一个文件, 命名如 "section_001.txt" num_sections = 50; % 断面数量 sections = struct('z', [], 'B', [], 'A', []); for i = 1:num_sections filename = sprintf('section_%03d.txt', i); data = load(filename); sections(i).z = data(:,1); % 高程向量 sections(i).B = data(:,2); % 宽度向量 sections(i).A = data(:,3); % 面积向量 end % 断面间距 (单位: m) dx = 200 * ones(1, num_sections-1); % 假设均匀, 实际应从桩号计算这段代码把每个断面的地形数据读入结构体。z是高程数组，从河底到堤顶；B是对应高程的水面宽；A是过水面积。计算时根据当前水位插值得到实际的水面宽和面积。注意高程数据要统一基准面，不同断面用不同基准面是新手常犯的错误，会导致水面线出现台阶。
断面间距dx不要求均匀，但相邻断面的间距变化不宜超过 20%，否则离散误差会增大。如果断面数据来自 CAD 或 GIS，导出时注意坐标系和单位，高程用米、距离用米，别混用厘米。
3.2 糙率率定：模型能不能用就看这一步
糙率 n 是水动力模型唯一的率定参数，也是最玄学的部分。理论上糙率反映河床阻力，和床沙粒径、植被、河道弯曲度有关。实际上，你永远无法精确知道每个断面的糙率，只能通过实测水位流量过程反推。
率定的基本流程是：先给一个初始糙率分布（比如全河段 0.03），跑一遍模型，把计算水位和实测水位对比，然后调整糙率。如果计算水位偏高，说明糙率给大了，阻力太大；如果偏低，说明糙率给小了。调整时不要全河段统一改，要分段调，因为不同河段的阻力特性不同。
% 糙率率定: 简单的手动调整循环 n_rough = 0.03 * ones(1, num_sections-1); % 初始糙率 n_obs = [0.028, 0.032, 0.035, ...]; % 根据河段特性分段的经验值 for iter = 1:max_iter % 运行模型 [Z_sim, Q_sim] = run_model(n_rough, boundary, sections, dx, dt); % 计算与实测水位的误差 Z_obs = load('observed_stage.txt'); rmse = sqrt(mean((Z_sim - Z_obs).^2)); % 如果误差满足要求则退出 if rmse < tolerance fprintf('率定完成, RMSE = %.4f m\n', rmse); break; end % 否则调整糙率 (简化示意, 实际需根据误差分布分段调整) n_rough = n_rough + 0.001 * sign(Z_obs - Z_sim); end这段代码展示的是最朴素的手动率定思路。实际工作中，糙率调整需要结合洪峰传播时间、水位过程线形状综合判断。如果洪峰到达时间偏早，说明上游糙率偏小；如果洪峰坦化严重，说明糙率偏大。我一般会先调上游段，再调下游段，因为上游误差会向下游传播。
注意：糙率不是万能调节器。如果模型结构有误（比如断面数据错误、边界条件不合理），再怎么调糙率也拟合不上。先检查数据，再调参数。
3.3 结果验证：不是跑通就行
模型跑通不等于结果可信。验证至少要做三件事：水量平衡检查、水位过程线对比、洪峰传播时间校核。
水量平衡检查是最基本的。计算时段内，上游流入总量减去下游流出总量，应该等于河道槽蓄量的变化。如果误差超过 5%，说明离散格式或边界处理有问题。
% 水量平衡检查 V_in = trapz(t, Q_upstream); % 上游总入流 V_out = trapz(t, Q_downstream); % 下游总出流 V_storage_change = sum(A_total_end - A_total_start); % 槽蓄变化 balance_error = abs(V_in - V_out - V_storage_change) / V_in * 100; fprintf('水量平衡误差: %.2f%%\n', balance_error);水位过程线对比是把计算水位和实测水位画在同一张图上，看相位和幅值是否一致。相位差反映糙率或河道长度误差，幅值差反映断面形状或边界条件误差。
洪峰传播时间校核是看洪峰从上游到下游的到达时间是否合理。天然河道洪峰传播速度一般在 1 到 3 m/s，如果计算出来的传播速度是 10 m/s，那肯定是哪里搞错了。
4. 避坑与排查：那些让模型崩溃的细节
4.1 常见问题一：矩阵奇异或求解失败
现象：运行时报错「Matrix is singular to working precision」或者求解结果全是 NaN。
原因：边界条件没有正确闭合，或者某个断面的水面宽为零导致除零。常见于干河床初始条件，水位低于河底高程时水面宽为零。
解决：检查边界条件是否替换了对应的方程，确保方程组满秩。初始条件给一个最小水深（比如 0.1 m），避免干河床。如果河道有漫滩，确保断面数据覆盖到最高水位以上。
4.2 常见问题二：水位出现非物理振荡
现象：计算结果中水位随时间剧烈振荡，出现负水深或锯齿状波动。
原因：Preissmann 格式在 θ<0.5 时可能产生振荡，或者时间步长过大导致数值色散。另外，断面间距突变也会引起反射振荡。
解决：把 θ 调到 0.6 以上，最好 0.7。减小时间步长试试，如果振荡消失说明是时间离散问题。检查断面间距，相邻断面间距比不要超过 1.5 倍。
4.3 常见问题三：洪峰坦化过度
现象：计算出的洪峰比实测洪峰低很多，过程线变得很平缓。
原因：数值耗散过大，通常是因为 θ 取 1.0 且时间步长太大。另外，糙率给大了也会加速洪峰坦化。
解决：把 θ 降到 0.6 到 0.7，减小时间步长到 60 秒试试。如果还是坦化，检查糙率是否合理，用实测洪峰传播时间反推糙率。
4.4 常见问题四：下游边界反射
现象：下游边界附近水位出现异常波动，向上游传播。
原因：下游边界条件与内部计算格式不兼容，或者边界位置太靠近感兴趣区域。
解决：下游边界尽量放在实测站位置，不要凭空外推。如果必须外推，用水位-流量关系而不是固定水位。把下游边界向下游延伸几个断面，让反射波衰减掉。
4.5 常见问题五：初始条件导致长时间预热
现象：模型开始运行后，前几个小时结果不可用，需要很长时间才能稳定。
原因：初始条件给的是静水，与实际流态差异大，模型需要时间调整。
解决：用稳态计算的结果作为初始条件。先跑一个恒定流计算，把水位和流量作为非恒定流的初始场。这样预热时间可以缩短到几分钟。
5. 进阶用法：从单河道到河网，以及我踩过的那些坑
单河道模型跑通之后，下一步通常是河网。河网的核心难点在汊点连接：三条或更多河道交汇，水位相同但流量分配未知。常见做法是把汊点作为边界，用迭代法求解流量分配，或者把汊点水位作为额外未知量加入全局矩阵。
代码里实现河网，需要把每个河段的系数矩阵按连接关系组装成更大的稀疏矩阵。汊点处增加连续性方程（流入等于流出）和能量方程（水位相等）。MATLAB 的稀疏矩阵组装用sparse的 triplet 形式最方便：
% 河网汊点连接: 三条河段交汇于一个节点 % 河段1和河段2流入, 河段3流出 % 汊点方程: Q1 + Q2 = Q3, Z1 = Z2 = Z3 % 在全局矩阵中, 汊点对应的行替换为连接方程 % 假设汊点水位未知量索引为 iz_node, 流量索引为 iq1, iq2, iq3 A(iz_node, iq1) = 1; A(iz_node, iq2) = 1; A(iz_node, iq3) = -1; b(iz_node) = 0; % 水位相等方程 A(iz_node+1, iz1) = 1; A(iz_node+1, iz2) = -1; b(iz_node+1) = 0;河网计算最容易翻车的地方是汊点初始流量分配。如果初始分配太离谱，迭代可能不收敛。我一般用河道断面面积比作为初始分配比例，效果比较稳。
另一个进阶方向是耦合水质或泥沙。水动力是基础，水质模型在水动力结果上做对流扩散。如果只是做洪水演进，单河道模型足够；如果要做污染物输移，就需要在水动力基础上加输移方程。
最后说一个血泪经验：每次改完代码，先跑一个最简单的测试案例——比如矩形棱柱河道，有解析解的那种。确认代码没改坏，再上实际河道。我见过太多人直接在实际数据上调试，出了问题根本分不清是代码错还是数据错。从那以后我每次改完核心求解器，都强制走一遍矩形河道测试，对比解析解，误差在 1% 以内才继续。希望帮到你。
本文还有配套的精品资源，点击获取
