简介三层蒸发蓄满产流模型新安江模型Python计算程序是一份面向水文水资源专业学生、科研人员与水利工程师的代码资源用于实现基于降水与蒸发能力数据的三层蒸发蓄满产流计算可直接运行并输出产流结果。资源仅包含1个Python脚本.py包体仅1KB结构精简代码逻辑聚焦三层蒸发与蓄满产流核心算法便于阅读修改和二次开发。需要配合作者另行上传的PEdata.xlsx数据文件使用该文件提供降水、蒸发能力等输入数据适合教学演示、课程设计或实际项目快速验证。目前已有2642人学习说明该代码在相关领域具有较好的参考价值。通过运行该程序读者可以掌握新安江模型三层蒸发模式的程序化实现思路理解蓄满产流计算流程并在此基础上扩展参数率定或与其他水文模型耦合。1. 暴雨洪水的账可以提前算新安江模型的 Python 落地做水文预报的人对“产流”这两个字都不会陌生。同样一场暴雨落在山区和落在平原水网形成的洪水过程线完全不同同一个流域前期下了十天雨和旱了两个月暴雨的产流率可能差出一倍多。新安江模型出现的意义就是解决这样一类问题在缺乏完整物理方程描述的流域上用一套“蓄满产流 三层蒸散发”的经验性结构把降雨转成径流的过程算出来。这套 Python 计算程序针对的是新安江模型中三层蒸发蓄满产流这一主流程把原来靠 Excel 或手算的逐日、逐时水文账变成可复现、可批量率定的代码。它是国内洪水预报、水文模拟中最常被拿来做二次开发的底子适合做水文预报的工程师、搞模型复现的研究生以及想把传统水文模型移植进现代数据管线的人。2. 把模型拆成四个计算块蒸散发、产流、分水源、汇流的代码实现新安江模型从结构上来讲并不复杂它把流域产流过程抽象成“蒸散发—产流—水源划分—汇流”四个环节。Python 程序也是按这四个环节来组织的。动手改代码之前先把每一段的物理含义和计算顺序弄清楚后面调参才不会像在黑匣子里乱摸。2.1 三层蒸散发先算上层再补下层深层只在极端情况启动三层蒸散发是整套模型的“扣损”环节它回答的问题很简单降下来的水有多少先被蒸发掉了剩下的才能进入产流计算。它的核心逻辑是分三层蓄水结构上层WUM、下层WLM、深层WDM蒸发时先抽上层上层不够再抽下层下层还不够才动深层。深层蒸发有一个额外系数 C它的作用是限制深层水的蒸发强度避免模型把深层水当免费的水库随便抽。程序里对应的蒸散发计算块长这样def calc_evaporation(wu, wl, wd, p, e0, wum, wlm, c, kc): # 单元流域蒸散发能力kc为折算系数 ep e0 * kc pe p - ep # 上层蒸散发优先消耗上层张力水 if pe 0: eu ep else: eu max(0.0, wu pe) # 剩余蒸发能力给下层 rest ep - eu if rest 0: # 下层蒸散发受下层蓄水和上层蓄水比例约束 el rest * wu / wum if el wl: el wl else: el 0.0 # 深层蒸散发只有上下层都干了才启动 ed 0.0 rest2 ep - eu - el if rest2 0: # 深层蒸发强度受系数c限制 ed c * rest2 if ed wd: ed wd return eu, el, ed这一段代码里最值得注意的地方是两个第一个是el rest * wu / wum这个公式它表达的是“下层蒸发能力与上层蓄水程度成正比”——上层土越湿蒸发就越容易穿透到下层第二个是深层的c系数它默认取 0.1 左右作用是让深层蒸发非常缓慢。实际率定时先调kc让总量蒸发对上再动c修正极端干旱期的低水模拟。2.2 张力水蓄泄与产流判断蓄水容量曲线的面积逻辑产流计算是整套模型的心脏它的物理基础是“蓄满产流”流域内每一点都有一个张力水容量降雨先补充张力水补满了才产流。但流域内各点的蓄水容量不一致所以引入了一条蓄水容量曲线用参数B控制曲线形状。WM WUM WLM WDM是流域平均张力水容量W是当前张力水蓄量PE P - 蒸发是扣损后的净雨。一个时段内的产流量通过面积积分来计算def calc_runoff(w, pe, wm, b, im): # 不透水面积直接产流 runoff_imp im * pe # 透水面积按蓄水容量曲线积分 a (1.0 b) * wm if pe 0: if w pe wm: runoff pe - (wm - w) wm * (1.0 - ((1.0 - w / wm) ** (1.0 b))) / (1.0 b) else: runoff (pe - (wm - w)) wm * (((w pe) / wm) ** (1.0 b) - (w / wm) ** (1.0 b)) / (1.0 b) else: runoff 0.0 # 返回总产流与不透水面积产流 return max(0.0, runoff), runoff_imp这一段是很多人第一次接触会看懵的地方两个分支分别对应“PE 把整个流域蓄满”和“只蓄满一部分”的情况。第一个分支里wm * (1.0 - ((1.0 - w / wm) ** (1.0 b))) / (1.0 b)是蓄水容量曲线在最大蓄水处的面积第二个分支是在当前蓄水W处到W PE处的面积增量。IM是流域内不透水面积比城市流域可以给到 0.3 以上天然流域一般在 0.01 到 0.05 之间。2.3 三水源划分自由水蓄水库怎么把总产流切成三份产流算出来之后还要把它分成地表径流、壤中流、地下径流三部分。新安江模型用一个“自由水蓄水库”来完成这个切分总产流先进入自由水蓄水库蓄满后溢出成为地表径流同时水库向壤中流和地下径流两个出口出流出流系数分别由KI和KG控制。分水源后的三部分水流汇入河网的速度不同所以后续汇流计算会分别处理。典型代码如下def calc_split(r, sm, ex, ki, kg): # 自由水蓄水库分水源 # 先按蓄水容量曲线求产流面积比例fr if r 0: fr 1.0 - (1.0 - ((sm * (1.0 ex)) / (sm * (1.0 ex) r))) ** (1.0 / (1.0 ex)) else: fr 0.0 # 地表径流自由水蓄满部分溢出 rs max(0.0, r - sm * fr) # 壤中流和地下径流按出流系数分配 ri ki * sm * fr rg kg * sm * fr return rs, ri, rg这里ex是自由水蓄水容量曲线的方次一般取 1.5 左右sm是自由水容量它的大小直接影响洪峰的形状。sm调大自由水蓄水库能存住更多水壤中流和地下径流占比更高洪峰被拉平sm调小地表径流占比大洪峰更尖更高。这个参数对过程线形状的影响非常敏感是我每次率定最先试的参数之一。2.4 汇流计算单元面积上的坡面汇流与河网演进汇流部分在新安江模型里有两种实现层次单元流域上的坡面汇流和河道汇流。单元流域的汇流通常用分段线性水库或者瞬时单位线来做把三水源分别演算到流域出口河道汇流则常用马斯京根法把上一段的出流过程平滑地推进到下一段。汇流阶段代码示例以线性水库为例def calc_concentration(qs_in, qi_in, qg_in, cs, ci, cg, dt): # 三水源分别做消退演算 qs_out qs_in * cs qi_out qi_in * ci qg_out qg_in * cg # 单元流域出口总流量 total_q qs_out qi_out qg_out return total_q, qs_out, qi_out, qg_outcs、ci、cg分别对应地表、壤中流、地下径流的消退系数取值在 0.9 到 0.998 之间。地表水消退快取低值地下水消退慢取接近 1 的值。这个代码块看起来简单但它是在时段单位比如 1 小时或 6 小时下的状态更新时间步长变了消退系数必须跟着重率定——这是很多人直接把参数从日模型搬到小时模型后模拟结果惨不忍睹的根本原因。3. 参数率定与物理意义三水源模型的 12 个关键参数怎么调三段式产流结构背后是十几个参数。新手最容易犯的错是拿到程序就急着跑优化算法把 12 个参数一起丢给程序自动寻优。结果往往收敛不到合理值就算收敛了参数之间互相补偿物理意义全丢。我一般的做法是分阶段率定先蒸散发再产流最后分水源和汇流。3.1 参数表每个参数的物理量纲、范围与调试方向先列一份完整的参数清单调参的时候对照这张表看不容易迷失方向参数物理意义取值范围调试方向KC蒸散发折算系数0.5 ~ 1.5控制总水量平衡偏大则总量偏干WUM上层张力水容量(mm)10 ~ 30上层蒸发敏感影响雨季初期的蒸散发WLM下层张力水容量(mm)50 ~ 100影响连续无雨期的退水WDM深层张力水容量(mm)10 ~ 30影响极端干旱模拟C深层蒸散发系数0.05 ~ 0.2干旱期低水调节WM平均张力水容量(mm)WUMWLMWDM由三个分项加和得到B蓄水容量曲线方次0.1 ~ 0.5影响产流面积分布越大越不均匀IM不透水面积比例0.01 ~ 0.3城市流域调大SM自由水容量(mm)5 ~ 50最敏感参数之一影响洪峰形态EX自由水曲线方次1.0 ~ 2.0一般取 1.5敏感度低KI/KG自由水出流系数0.1 ~ 0.7KIKG决定壤中流和地下径流出流速度CS/CI/CG消退系数0.9 ~ 0.998时间步长变化时必须重率定3.2 率定顺序把参数分成“总量控制”和“形状控制”两组我的率定顺序分三步走。第一步只调KC、WUM、WLM、WDM蒸散发组目标是让多年累计径流量与实际基本吻合第二步调B、IM、SM、KI、KG产流与分水源组目标是让场次洪水的洪峰量级和峰现时间对上第三步才碰CS、CI、CG汇流参数目标是优化洪水过程线的形状。调参时最好固定其中一个参数组只动另一个。比如调产流参数时蒸散发参数固定不变这样每次改动的影响才能从模拟结果中分辨出来。如果所有参数一起动经常出现“蒸发多了但产流参数也改了最终结果碰巧很好”的假象换个时段验证立刻露馅。3.3 数据输入格式整理降雨、蒸发、实测流量怎么对齐程序能跑起来的前提是输入数据格式能对上。我通常准备三个 CSV 文件降雨序列、蒸发序列用 E601 蒸发皿或卫星反演数据折算以及用于率定的实测流量序列。时间轴必须完全对齐缺测数据不能简单填 0否则程序会把“缺测”当成“无雨”引入虚假的干旱期。读取数据的代码建议用 pandas 统一处理import pandas as pd # 读取降雨、蒸发、实测流量时间列统一为索引 rain pd.read_csv(rain.csv, index_coltime, parse_datesTrue) evap pd.read_csv(evap.csv, index_coltime, parse_datesTrue) qobs pd.read_csv(qobs.csv, index_coltime, parse_datesTrue) # 对齐到相同时间索引缺测位置用前向填充阈值判断 df pd.concat([rain[rain], evap[e0], qobs[q]], axis1) df.columns [p, e0, q] df df.fillna(methodffill).dropna() # 设置初始土壤含水量为0.5倍田间持水量 wu 0.5 * wum wl 0.5 * wlm wd 0.0时间对齐这一步是最容易翻车的CSV 里的时间戳如果一个是字符串、一个是 datetime 对象pd.concat时索引对不上后面所有结果全部错位。我的习惯是先df.index pd.to_datetime(df.index)强制转一次再检查df.isna().sum()输出缺测数量确认没有隐藏的 NaN 再往下走。3.4 输出数据与水量平衡一段代码检查你的模拟是否“漏了水”率定之前先做一个水量平衡检查能过滤掉一大半的数据问题。水量平衡公式是降雨总量 - 蒸发总量 - 径流总量 流域蓄水变化量。在长序列模拟里蓄水变化量很小如果算出来差了一大截说明代码有 bug 或者参数错得离谱。total_p df[p].sum() total_e df[e0].sum() * kc total_q q_sim.sum() storage_change total_p - total_e - total_q print(f总降雨: {total_p:.1f} mm) print(f总蒸发: {total_e:.1f} mm) print(f总径流: {total_q:.1f} mm) print(f蓄水变化: {storage_change:.1f} mm)storage_change一般应该在一个很小的范围内。如果出现几百毫米的偏差优先检查是不是初始土壤含水量设得太离谱。初始含水量设为 0 的话模型需要跑相当长的预热期才能让土壤湿度达到合理水平前几年模拟结果会系统性偏低——这就是有名的预热期spin-up问题。4. 运行排查与避坑笔记五类高频报错和数据异常实录三层蒸发蓄满产流模型在代码层面并不复杂但真正跑起来之后各种“看起来没错但结果不对”的情况非常多。我把自己反复踩过的坑整理成几条都是现象→原因→解决的完整链路。4.1 模拟结果整体偏低洪水过程线被“削平”了现象模拟径流总量比实测少 30% 以上洪水过程线明显比实测平坦但峰现时间基本吻合。原因这是我最常遇到的一种情况。排查后发现多数时候是初始土壤含水量设得过低模型在前几年一直处于“先填土壤水库”的状态降雨被大量吸收根本没有多余的水产流。还有一个可能是IM不透水面积设得太小城市型流域的快速产流没有被表达出来。解决把初始WU、WL设为田间持水量的 50% 左右模型最开始的 1 到 2 年强制作为预热期丢弃不参与率定。检查IM是否贴近流域实际下垫面情况有条件的用卫星遥感不透水面积数据核对一下。做完这两步绝大多数“整体偏低”问题能解决。4.2 模拟流量出现负值现象明明降雨和蒸发都是正数模拟流量序列里却出现了负数程序在某一步直接崩溃。原因自由水蓄水库在分水源时如果处理不当可能出现sm * fr大于当前蓄水量的情况这时r - sm * fr变成负数地表径流被算成了负值。还有一个隐蔽的原因是蒸散发计算里el rest * wu / wum没有做上限约束下层蒸发把下层张力水抽成了负数。解决在分水源函数里强制rs max(0.0, r - sm * fr)同时在蒸发函数里对el和ed都加上限约束。我一般还会在代码入口处加一个统一的非负检查任何负值直接断言报错宁可崩掉也不要带着负数往下算否则错误会被汇流段的消退系数一路放大。4.3 率定算法把参数调到物理范围之外现象用自动优化算法比如 scipy 的 differential_evolution 或遗传算法率定时SM被优化到 90B变成 2.5模拟结果看似很好但换一个时段验证立刻崩溃。原因优化算法只知道让目标函数最小化不知道参数的物理意义。如果参数范围给得太宽算法会用极端的参数组合去“硬拟合”训练期的数据产生严重的过拟合。尤其是SM和B之间的相互作用可以在训练期互相补偿到惊人的程度。解决优化前把所有参数锁死在物理范围里。我的做法是把参数范围写成字典优化算法只在这个字典定义的范围内搜索而且在目标函数里加上“参数越界惩罚项”任何越过物理范围的候选解直接给一个极大的损失值比在约束条件里写死更稳。4.4 时间步长从日模型切到小时模型所有参数失效现象同一个流域日模型率定的参数拿去做 1 小时步长模拟洪峰被严重高估退水过程乱成一团。原因新安江模型的参数与时间步长强相关。CS、CI、CG是时段消退系数它们的数学形式决定了取值必须随时间步长变化SM自由水容量虽然量纲是毫米但在不同步长下其“等效容量”也需要重调。日模型的参数直接搬进小时模型等于用低频标定的体系去算高频过程必然出问题。解决切换时间步长后把CS、CI、CG按公式换算CS_hour CS_day ** (dt_hour / dt_day)然后重新率定SM和KI/KG。蒸散发参数可以保持基本不变但汇流和自由水水库相关参数必须从头来过。4.5 降雨数据缺测被当成 0模型旱得离谱现象模拟径流在汛期突然断流但实测流量明明还有持续的基流检查输入数据发现某雨量站连续 20 天缺测代码把缺测当成了 0。原因数据预处理阶段直接把 NaN 填充成了 0。这在数学上合法但对水文模拟来说等于人为制造了一个虚假的干旱期土壤水库会被模型抽干后续场次洪水的产流条件彻底失真。解决把缺测分成两类处理。短时间缺测比如 6 小时以内用前后时段的线性插值长时间缺测直接标记为无效时段模拟时跳过这段时间或不参与率定评分。我还在代码里加了一个判断任何时段p为 0 且e0大于 0 时检查前后 24 小时的累计降雨如果前后都在下雨而中间为零大概率是缺测打一个警告。5. 模型检验与结果验证确定性系数、水量平衡与敏感性分析怎么做模型跑通了不等于模型能用了。新安江模型的模拟结果要真正用于洪水预报至少过三关过程线形状对不对、总量平不平衡、参数稳不稳定。这一章把三关的检查方法逐一过一遍。5.1 确定性系数 NSE一个数字评价过程线拟合好坏纳什效率系数NSE是水文模拟最常用的评价指标公式是NSE 1 - Σ(Qobs - Qsim)^2 / Σ(Qobs - Qmean)^2。NSE 接近 1 说明模拟完美小于 0 说明模拟比直接取平均还差。import numpy as np def nse(qobs, qsim): # 去掉无效时段后计算NSE mask ~np.isnan(qobs) ~np.isnan(qsim) qo qobs[mask] qs qsim[mask] if len(qo) 10: return float(nan) qo_mean np.mean(qo) r 1.0 - np.sum((qo - qs) ** 2) / np.sum((qo - qo_mean) ** 2) return rNSE 在 0.7 以上算基本合格0.85 以上算良好。但如果只看 NSE 会踩一个坑NSE 对洪峰非常敏感只要洪峰对准了即便枯水期模拟得一团糟NSE 数字还是很好看。所以我通常同时看相对误差RE和洪峰误差RE_peak不要单一指标定生死。5.2 分时段检验汛期、枯水期、场次洪水分开算全局 NSE 会掩盖掉很多局部问题。同一套参数可能在汛期表现优秀但在枯水期完全失真原因是基流的模拟主要受CG和深层蒸发参数控制而汛期暴雨的产流主要受SM控制两组参数被全局 NSE 评价时互相稀释。# 按月份拆分检验 month pd.to_datetime(qobs.index).month flood month.isin([5, 6, 7, 8, 9]) # 汛期月份 dry month.isin([11, 12, 1, 2]) # 枯水期月份 nse_flood nse(qobs[flood], qsim[flood]) nse_dry nse(qobs[dry], qsim[dry]) print(f汛期 NSE: {nse_flood:.3f}) print(f枯水期 NSE: {nse_dry:.3f})如果汛期 NSE 高、枯水期 NSE 低先调CG和C深层蒸发系数反过来枯水期好但汛期差重心放在SM、KI、KG上。场次洪水的检验更严格选 10 场以上的大洪水单独计算洪峰误差和峰现时间差平均绝对误差控制在 10% 以内峰现时间差不超过一个时段才敢拿去用。5.3 敏感性分析找出“动一点就天翻地覆”的参数敏感性分析的目的是搞清楚哪些参数对结果的影响最大。做法很简单固定其他参数不变把某个参数在其范围内上下浮动 10%观察 NSE 的变化幅度。变化幅度大的参数就是敏感参数率定时优先保证它们准确不敏感的参数固定在最合理的物理值即可不要浪费算力去精细调优。def sensitivity_analysis(base_params, vary_keys, pct0.1): result {} base_nse run_model(base_params) for key in vary_keys: params_up base_params.copy() params_down base_params.copy() params_up[key] base_params[key] * (1 pct) params_down[key] base_params[key] * (1 - pct) nse_up run_model(params_up) nse_down run_model(params_down) result[key] max(abs(nse_up - base_nse), abs(nse_down - base_nse)) return result # 对SM、B、CS、CG做敏感性分析 sens sensitivity_analysis(params, [SM, B, CS, CG])实践下来SM、CS、KG是最敏感的B和EX敏感度较低。如果自动率定算法同时调SM和B很容易出现它们互相抵消的现象这时候把B锁死在 0.3只让算法调SM结果更稳定。5.4 参数移植性原则从一个流域到另一个流域的边界条件把某流域率定好的参数搬到邻近流域时有一条重要的引用原则蒸散发组参数KC、WUM、WLM、WDM可以基本沿用因为它们反映的是气候和下垫面的大尺度特征产流参数中的SM、KI、KG必须重新率定因为它们对土壤结构、地形坡度极其敏感。IM则要单独评估城市扩张前后的流域IM可能相差数倍。实操时我一般先把邻流域的实测洪水场次挑 5 场用原参数直接跑一遍看 NSE 和水量平衡误差的初始水平再决定哪些参数需要重新率定而不是把整套参数推倒重来。6. 从单流域到批量模拟脚本化率定与过程线绘制的完整技巧模型在单流域上跑通之后真正的产出效率来自批量化。水文分析很少只做一个流域一上来就是十几个甚至几十个断面。每个断面单独手工率定不现实这时候要把模型封装成可批量调用的函数配合自动率定算法和可视化出图形成一套“输入数据→自动率定→输出图表”的处理链路。我习惯把所有流域的配置写成一个 CSV每一行是一个流域 ID、参数初值、率定时段。读取这个 CSV 后循环调用模型import json, csv from scipy.optimize import differential_evolution def optimize_for_basin(basin_id, params_init, bounds, rain, evap, qobs): def objective(x): params dict(zip(param_keys, x)) qsim run_newanjiang_model(rain, evap, params) return -nse(qobs, qsim) # 取负号最小化即最大化NSE result differential_evolution(objective, bounds, maxiter300, seed42) return dict(zip(param_keys, result.x))批量跑完之后出图是关键。我一般不直接看数字而是把模拟和实测过程线叠在一起画。一张好的对比图比 10 个评价指标更能暴露问题。import matplotlib.pyplot as plt def plot_hydrograph(qobs, qsim, title): plt.figure(figsize(12, 5)) plt.plot(qobs.index, qobs.values, k-, lw1.2, labelObserved) plt.plot(qsim.index, qsim.values, r--, lw1.0, labelSimulated) plt.legend() plt.title(title) plt.ylabel(Discharge (m³/s)) plt.grid(alpha0.3) plt.tight_layout() plt.savefig(f{title}.png, dpi150)绘图时我习惯把 NSE 和水量平衡误差直接写在图标题或图注上这样每次检查结果一眼就能看到关键指标。图中如果出现“模拟过程线始终慢半拍”的视觉特征优先怀疑汇流参数如果洪峰高度一致但宽度偏胖优先检查SM是不是偏大。最后讲一个我自己的教训以前做参数移植时把某流域的完整参数组直接套到邻域NSE 只有 0.3当时以为是数据问题反复检查输入数据耗了两天。最终发现是SM没重率定邻域的土壤偏砂质SM应该调小我却沿用了原流域的大值。从那以后我每次做参数移植都强制走一遍敏感性分析先把最敏感的三个参数识别出来重新率定再逐步放开其余参数。希望这套方法对你的模型调试同样有帮助。本文还有配套的精品资源点击获取