简介：三层蒸发蓄满产流模型（新安江模型）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_out

cs、ci、cg分别对应地表、壤中流、地下径流的消退系数，取值在 0.9 到 0.998 之间。地表水消退快，取低值；地下水消退慢，取接近 1 的值。这个代码块看起来简单，但它是在时段单位（比如 1 小时或 6 小时）下的状态更新，时间步长变了，消退系数必须跟着重率定——这是很多人直接把参数从日模型搬到小时模型后模拟结果惨不忍睹的根本原因。

3. 参数率定与物理意义：三水源模型的 12 个关键参数怎么调

三段式产流结构背后是十几个参数。新手最容易犯的错，是拿到程序就急着跑优化算法，把 12 个参数一起丢给程序自动寻优。结果往往收敛不到合理值，就算收敛了，参数之间互相补偿，物理意义全丢。我一般的做法是分阶段率定，先蒸散发，再产流，最后分水源和汇流。

3.1 参数表：每个参数的物理量纲、范围与调试方向

先列一份完整的参数清单，调参的时候对照这张表看，不容易迷失方向：

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_col='time', parse_dates=True) evap = pd.read_csv('evap.csv', index_col='time', parse_dates=True) qobs = pd.read_csv('qobs.csv', index_col='time', parse_dates=True) # 对齐到相同时间索引，缺测位置用前向填充+阈值判断 df = pd.concat([rain['rain'], evap['e0'], qobs['q']], axis=1) df.columns = ['p', 'e0', 'q'] df = df.fillna(method='ffill').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被优化到 90，B变成 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 r

NSE 在 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, pct=0.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, maxiter=300, seed=42) 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-', lw=1.2, label='Observed') plt.plot(qsim.index, qsim.values, 'r--', lw=1.0, label='Simulated') plt.legend() plt.title(title) plt.ylabel('Discharge (m³/s)') plt.grid(alpha=0.3) plt.tight_layout() plt.savefig(f'{title}.png', dpi=150)

绘图时我习惯把 NSE 和水量平衡误差直接写在图标题或图注上，这样每次检查结果一眼就能看到关键指标。图中如果出现“模拟过程线始终慢半拍”的视觉特征，优先怀疑汇流参数；如果洪峰高度一致但宽度偏胖，优先检查SM是不是偏大。

最后讲一个我自己的教训：以前做参数移植时，把某流域的完整参数组直接套到邻域，NSE 只有 0.3，当时以为是数据问题，反复检查输入数据耗了两天。最终发现是SM没重率定，邻域的土壤偏砂质，SM应该调小，我却沿用了原流域的大值。从那以后我每次做参数移植都强制走一遍敏感性分析，先把最敏感的三个参数识别出来重新率定，再逐步放开其余参数。希望这套方法对你的模型调试同样有帮助。

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