简介这是一套基于MATLAB的一维水质模型模拟程序包面向环境工程、水利及科研人员用于解决河流、渠道等一维流动水体中污染物输移与水质变化预测问题。包内共2个文件均为.m源文件包含shuizhi.m及Untitled.m代码涵盖河流参数定义、离散网格构建、有限差分格式选择、边界条件设定及时间步长迭代求解等关键环节可供相关专业学生与工程师参考运行或二次开发。压缩包仅1KB体型轻量便于快速部署验证。当前已有285人学习适合具备水文学与水力学基础、需要初步掌握一维水质差分建模的MATLAB使用者。通过阅读和调试这两段程序使用者可直观理解一维水质模型从方程离散到数值求解的完整流程并据此开展扩展研究或教学演示。1. 一维水质模型与一维差分法河道模拟先从一维算起拿到shuizhi.zip这类资源包里面往往塞满了脚本、数据和一个又一个调试记录。拆开看真正核心的东西大多不是某个现成软件而是“用一维差分法求解一维水质模型”这行字背后的数学处理。河道水质预测并不总需要立刻上 EFDC、MIKE 这类大型三维模型在一条断面混合均匀的中小河流上研究污染物浓度沿程变化、估算纳污能力、比选排放口位置一维对流扩散方程就是最常用的理论框架而有限差分法是把方程变成可执行代码的最短路径。本文从方程离散开始写出一套能运行、能调参、能验证的 Python 实现并给出稳定性控制和参数率定的实操经验。无论是刚接触水质模型的工程师还是准备二次开发自己工具的算法工程师都能直接照着改。2. 一维水质模型的控制方程与有限差分离散格式2.1 一维对流扩散方程中的三个物理项一维水质模型的控制方程由对流、扩散和衰减三个物理过程组成。对于一条沿 x 方向流动的河道若过水断面面积 A 恒定污染物浓度 C(x,t) 满足$$\frac{\partial C}{\partial t} -u\frac{\partial C}{\partial x} D\frac{\partial^2 C}{\partial x^2} - kC$$其中 u 是断面平均流速m/sD 是纵向离散系数m²/sk 是一级衰减系数1/d。对流项描述污染物随水流平移扩散项描述浓度梯度引起的弥散衰减项代表降解、挥发等去除过程。实际应用中还会增加源汇项 S例如点源排放或支流汇入此时方程右侧再加 S。这个方程在数学上属于抛物线型偏微分方程解析解只在恒定流、恒定源强且边界无限等理想条件下存在。真实河道断面变化、流速波动、源强随时间变化必须用数值方法求近似解。一维差分法的基本思路是把连续的 x 和 t 离散成网格节点用节点浓度的线性组合替代偏导数得到一个可迭代推进的代数方程组。2.2 显式、隐式与 Crank-Nicolson 格式的取舍对时间导数常见做法是用一阶前向差分$$\frac{\partial C}{\partial t} \approx \frac{C_i^{n1} - C_i^n}{\Delta t}$$对空间导数显式格式在 n 时刻取值隐式格式在 n1 时刻取值。三种常用格式的离散方程如下显式FTCS$$C_i^{n1} C_i^n - \frac{u\Delta t}{2\Delta x}(C_{i1}^n - C_{i-1}^n) \frac{D\Delta t}{\Delta x^2}(C_{i1}^n - 2C_i^n C_{i-1}^n) - k\Delta t C_i^n$$隐式BTCS需要解方程组把空间导数全部写为 n1 时刻$$C_i^{n1} \frac{u\Delta t}{2\Delta x}(C_{i1}^{n1} - C_{i-1}^{n1}) - \frac{D\Delta t}{\Delta x^2}(C_{i1}^{n1} - 2C_i^{n1} C_{i-1}^{n1}) k\Delta t C_i^{n1} C_i^n$$Crank-Nicolson 是显式和隐式的算术平均时间精度为二阶空间精度也为二阶但同样需要解方程组。选型取决于对稳定性和计算成本的权衡。显式格式每步只做简单赋值实现最快但时间步长受 Courant 数和扩散数约束隐式格式无条件稳定允许放大步长代价是每一步求解三对角方程组Crank-Nicolson 精度高但在浓度骤变的边界可能出现振荡。实际工程中如果网格较粗且源强变化平稳显式格式配合合适的步长完全够用要做长周期模拟或推进大时间步我一般会用隐式格式配合稀疏矩阵求解。2.3 差分格式的最小 Python 示意与参数表下面是最小化的显式离散核心代码只展示推进部分import numpy as np # 参数设置 dx 100.0 # 空间步长m dt 60.0 # 时间步长s u 0.3 # 流速m/s D 5.0 # 纵向离散系数m2/s k 1e-5 # 衰减系数1/s注意单位换算 # 稳定性检查 courant u * dt / dx diffus D * dt / (dx * dx) assert courant 1.0, Courant 数必须小于等于1 assert diffus 0.5, 扩散数必须小于等于0.5 # C 为浓度数组C_old 保存上一时刻 C_new C_old.copy() # 内部节点推进 C_new[1:-1] (C_old[1:-1] - 0.5 * courant * (C_old[2:] - C_old[:-2]) diffus * (C_old[2:] - 2 * C_old[1:-1] C_old[:-2]) - k * dt * C_old[1:-1])这段代码里C_old[2:]、C_old[:-2]分别对应 i1 和 i-1 节点的浓度。用数组切片代替显式循环计算速度接近编译语言。courant和diffus是两个无量纲数分别控制对流项和扩散项的稳定性。三种格式的特性对比如下表格式时间精度稳定性限制每步计算量适用场景显式 FTCS一阶Courant ≤ 1扩散数 ≤ 0.5最低教学演示、快速试算隐式 BTCS一阶无条件稳定低需解三对角阵长周期、大时间步模拟Crank-Nicolson二阶无条件稳定中对浓度峰值位置精度要求高的场景3. 用 Python 实现一维水质模型最小可运行代码3.1 网格划分与初始/边界条件设定构建一维网格时先把河道抽象成一条线段。设定河长 L 10 km从 x0 到 x10000 m空间步长 dx100 m网格节点数为 101。初始浓度取背景值 0.5 mg/L上游边界x0恒定浓度为 5 mg/L代表一个持续排放的污染源下游边界采用自由出口条件浓度梯度为零。完整求解之前需要把单位统一。水力停留时间、衰减系数常用天表示但计算中流速是 m/s时间步长是 s。若 k0.5 1/d换算成 1/s 需除以 86400。这类错误是 shuizhi.zip 里最常见的调试失败原因浓度曲线形状完全对但量级差几个数量级先查单位。3.2 显式有限差分求解器完整代码以下代码可以直接保存为water_quality_1d.py运行求解一个恒定源排放 12 小时后的浓度沿程分布import numpy as np import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt # 物理参数 L 10000.0 # 河长m dx 100.0 # 空间步长m u 0.3 # 流速m/s D 5.0 # 纵向离散系数m2/s k_per_day 0.5 # 衰减系数1/d k k_per_day / 86400.0 # 换算为 1/s # 数值参数 dt 60.0 # 时间步长s sim_time 12 * 3600.0 # 模拟时长12小时 n_steps int(sim_time / dt) # 网格 nx int(L / dx) 1 x np.linspace(0, L, nx) # 初始条件背景浓度 0.5 mg/L C np.full(nx, 0.5) # 上游恒定浓度边界 C[0] 5.0 # 稳定性检查 courant u * dt / dx diffus D * dt / (dx * dx) print(fCourant数{courant:.3f}, 扩散数{diffus:.3f}) assert courant 1.0 and diffus 0.5 # 时间推进 for n in range(n_steps): C_old C.copy() # 内部节点 C[1:-1] (C_old[1:-1] - 0.5 * courant * (C_old[2:] - C_old[:-2]) diffus * (C_old[2:] - 2 * C_old[1:-1] C_old[:-2]) - k * dt * C_old[1:-1]) # 上游边界保持恒定 C[0] 5.0 # 下游自由出口复制前一个节点浓度 C[-1] C[-2] # 绘制浓度沿程曲线 plt.plot(x / 1000.0, C, b-, labelt12h) plt.xlabel(距上游距离 (km)) plt.ylabel(浓度 (mg/L)) plt.legend() plt.savefig(concentration_profile.png, dpi150)运行后生成concentration_profile.png能看到浓度从 5 mg/L 沿程衰减扩散作用让峰值变平缓。这个求解器的核心是第 25 行的数组运算显式格式每步只需一次切片和加减乘除没有矩阵组装过程因此 12 小时模拟 720 步在普通笔记本上耗时可忽略。3.3 参数修改对浓度曲线的影响把u从 0.3 提高到 0.8峰面会更快向下游迁移衰减距离变长把D从 5 提高到 20浓度梯度变小上游边界影响范围增大峰形更宽平把k_per_day从 0.5 提高到 2.0沿程浓度快速下降下游可能接近背景值。建议每次只改一个参数观察曲线形态变化记录峰浓度和峰位置这样能快速建立参数敏感度概念。参数增大后的现象常见取值范围u浓度峰向下游移动速度加快0.1 ~ 1.5 m/sD浓度曲线更平滑峰值降低1 ~ 20 m²/sk整体浓度水平下降0.1 ~ 3.0 1/ddx空间分辨率变粗容易振荡10 ~ 500 m4. 一维差分法在水质模拟中的稳定性控制与标定4.1 Peclet 数与 Courant 数先算再跑网格距离和流速决定了数值解的“保真度”。Peclet 数定义为$$Pe \frac{u \Delta x}{D}$$它衡量对流项与扩散项的网格尺度之比。显式 FTCS 格式要求 Pe 不能过大工程上通常控制在 2 以内否则空间离散会产生数值振荡。Courant 数定义为$$Cr \frac{u \Delta t}{\Delta x}$$它表示一个时间步内物质迁移的网格距离比例。显式格式要求 Cr ≤ 1实际取值 0.5 左右较好。很多初学者的代码发散不是物理模型错误而是这两个数不满足约束。调试时第一时间打印这两个数比盯浓度曲线效率高。若 Pe 超过 2优先减小 dx而不是增大 D。dx 减半时扩散数会变为原来的 1/4因扩散数限制dt 也需要相应减小计算量按网格数乘步数上升。这是显式格式的天然瓶颈也是隐式格式在实际工程中占优的原因。4.2 边界条件与源项的处理坑上游边界是浓度边界条件代码里每次赋值即可。真正的坑在下游。用零梯度条件C[-1] C[-2]模拟无限远出口只适合出口处浓度已经接近均匀分布的情况。如果污染物运动到边界时浓度梯度仍然很大这个边界条件会反射虚假波导致下游区域浓度虚高或虚低。解决方法是把下游边界放足够远让模拟范围超出污染物实际迁移距离或改用对流边界条件。点源源项的处理也容易出错。若在 x2000 m 处有一个恒定流量 Q0.5 m³/s、浓度为 20 mg/L 的污水排入河流主流流量 Q_r10 m³/s排放口下游浓度应为C_mix (10.0 * 0.5 0.5 * 20.0) / (10.0 0.5) # 完全混合浓度然后再把C_mix赋给 x2000 m 对应的网格点。这里的常见错误是直接用C[20] C_mix覆盖原浓度值正确做法是每一时间步都让源项在背景浓度基础上叠加否则源强会被错误地随时间稀释。4.3 实测标定三步走参数率定、误差分析、结果校验实际项目里没有哪个参数能直接查表。我一般按下面三步做标定第一步手动率定。先固定 u 为断面实测流速D 根据经验公式估算初值用过去 5 天的实测浓度序列对比模拟值调整 k 让下游断面的衰减趋势匹配。第二步计算均方根误差RMSE和纳什效率系数NSE量化模拟效果。第三步敏感性分析每个参数上下浮动 20%重新模拟看哪个参数对峰值浓度影响最大优先率定它。# 误差评估代码 obs np.array([4.8, 3.2, 2.1, 1.4, 0.9]) sim np.array([4.5, 3.0, 2.3, 1.5, 0.8]) rmse np.sqrt(np.mean((obs - sim) ** 2)) nse 1 - np.sum((obs - sim) ** 2) / np.sum((obs - np.mean(obs)) ** 2) print(fRMSE{rmse:.2f} mg/L, NSE{nse:.3f})NSE 接近 1 表示模拟效果优秀低于 0 则模型还不如直接用均值预测需要重新考虑边界条件和参数初值。标定时不要只盯着 NSE还要检查浓度峰值的出现时间是否偏差过大时间偏差往往代表 u 设置不准确。5. 验证与封装的实用技巧用解析解给一维水质模型兜底5.1 稳态解析解验证数值解编写差分求解器时先用稳态解析解验证代码是最稳妥的做法。恒定流、恒定衰减、忽略下游反射的场景下一维稳态水质模型有解析解$$C(x) C_0 \exp\left[\frac{u x}{2D}\left(1 - \sqrt{1 \frac{4kD}{u^2}}\right)\right]$$运行数值模型到足够长时间达到稳态把x位置处的浓度与该解析解比较。差值小于 5% 说明离散化和边界条件正确。一个典型错误是把衰减系数 k 的代数值符号搞反解析解会立刻暴露这一点因为两条曲线方向相反。5.2 把求解器封装成命令行函数重复使用时把求解逻辑封装成函数输入参数用字典传递输出结果保存为 CSVdef run_water_quality(u, D, k_per_day, dt60.0, sim_hours24.0): dx 100.0 L 10000.0 # ... 内部实现与第三章相同... return x, C这样在参数寻优时可以批量循环不同k_per_day值自动计算 NSE 并选用最优参数。注意函数内部必须做 Courant 数和扩散数检查否则批量跑的时候隐蔽发散会浪费调试时间。5.3 最后一招用三对角矩阵求解隐式格式显式格式步长受限隐式格式才能放开步长做长周期模拟。隐式格式离散后得到一个三对角线性方程组用scipy.linalg.solve_banded或scipy.sparse求解都可以。这里给出基于稀疏矩阵的核心构造from scipy.sparse import diags # 组装三对角矩阵 A对应隐式格式 alpha D * dt / dx**2 beta u * dt / (4 * dx) # 中心差分对流项 gamma k * dt / 2 main_diag 1 2 * alpha gamma off_diag -alpha beta A diags([off_diag, main_diag, off_diag], offsets[-1, 0, 1], shape(nx, nx)).tocsr() # 时间推进内部节点 b C.copy() b[0] C[0] # 边界条件 C solver(A, b)这里的off_diag在对流项上取了中心差分Courant 数影响矩阵对角占优程度步长太大会产生振荡但不会发散。隐式格式允许时间步长放大到显式格式的 5 到 10 倍对长时间模拟特别有用。对比显式和隐式结果两者差异缩小到可接受范围后再交给下游水质分析使用。本文还有配套的精品资源点击获取