简介这是一份面向水文水资源领域的三层蒸发蓄满产流模型新安江模型Python计算程序适用于水利工程、水文预报及流域模拟方向的从业者、研究人员与学生帮助解决蓄满产流过程的核心计算与编程实现问题。压缩包内仅含1个py脚本文件整体约1KB程序结构紧凑、逻辑清晰适合直接阅读与二次修改。运行程序需同步使用配套的PEdata.xlsx数据文件其中包含计算蓄满产流所需的降水、蒸发能力等输入数据。通过本程序可系统理解新安江模型三层蒸发模式的计算流程、参数衔接关系及Python代码组织方式既可当作水文模型课程的教学参考也可作为科研中快速搭建产流计算模块的基础脚本。该资源已有2642人浏览学习适合需要上手新安江模型编程或开展流域水文模拟的读者下载使用。1. 三层蒸发蓄满产流模型这份 Python 代码为什么值得你认真跑一遍做水文模拟的人迟早会撞上一个尴尬局面拿着现成的降雨、蒸发、流量序列想算一个流域的产流量却不知道该选哪种模型。用超渗产流算湿润地区洪水过程畸大用统计回归算旬月尺度又没法解释日尺度的涨水退水。而标题里的「三层蒸发蓄满产流模型」也就是大家常说的新安江模型恰恰是湿润半湿润地区日尺度连续模拟最稳的一条路。它把流域看成一块会蓄水的海绵按张力水蓄量判断是否产流再按上、下、深三层分配蒸散发最终输出逐日径流过程。这篇文章不讲模型史直接讲怎么用 Python 把这套逻辑落地成可复现的程序适合正在做水文模拟、洪水预报或水资源评估的工程师和研究生。你会看到完整的参数表、逐行解释的代码以及我最开始跑这个程序时踩过的几个隐蔽的坑。2. 蓄满产流与三层蒸发的物理逻辑先立住新安江模型的四个概念2.1 蓄满产流的本质张力水蓄量与流域蓄水容量曲线新安江模型的核心假设是蓄满产流意思是只有当某点的土壤张力水蓄量达到田间持水量即蓄满后后续降雨才会产生径流。这个假设在湿润地区非常贴近实际因为湿润地区的土壤含水量常年偏高降雨稍大就容易达到蓄满状态地表径流和壤中流的生成主要由蓄满后的多余水量贡献。要描述整个流域并非同时蓄满模型引入了一条蓄水容量曲线。常用的是指数型或抛物线型分布表示流域内各点张力水容量 Wm 的空间不均匀性。曲线形态由参数 B 控制B0 意味着全流域蓄水容量均匀B 越大表示空间差异越大。初始时刻流域实际张力水蓄量 W 对应的最大点蓄水容量为 AU我一般把这个过程叫做「先找阈值再算产流」。实际的产流公式如下。[ WMM WM \times (1 B) ][ AU WMM \times \left(1 - \left(1 - \frac{W}{WM}\right)^{\frac{1}{1B}}\right) ]当降雨扣除蒸发后的净雨 PE 加上 AU 仍小于 WMM 时产流量 R 为[ R PE - (WM - W) WM \times \left(1 - \left(1 - \frac{AU PE}{WMM}\right)^{1B}\right) ]当 PE AU 大于等于 WMM 时全流域蓄满R PE - (WM - W)。这两个公式是整个产流计算的原点。很多人第一次实现时容易把 WM 和 W 搞混WM 是流域平均张力水容量是固定参数W 是当前时刻的实际张力水蓄量逐日更新。两者的差就是流域还能继续吸纳的水量。2.2 三层蒸发的分配逻辑上层、下层、深层各管什么蒸散发部分为什么非要分三层因为单层模型无法刻画土壤水分沿剖面的不均匀分布。晴天时表层土壤先干下层还有水如果只用一层模型会把整个土柱的水分均匀消耗导致表层很快耗尽、蒸发骤停和实测的连续蒸发过程对不上。新安江模型把张力水蓄量分成上、下、深三层各自有容量上限 WUM、WLM、WDM三层之和等于 WM。蒸散发的分配规则是先按蒸散发能力 KC 乘以蒸发皿蒸发 EP 得到流域蒸散发能力 EM优先从上层消耗。上层蓄水 WU 足够时上层蒸发 EU 就等于 EM下层和深层不蒸发。若 WU 不够上层先耗尽剩余蒸发能力转入下层下层蒸发 EL 受下层蓄水比例和深层系数 C 的限制公式为[ EL (EM - EU) \times \frac{WL}{WLM} ]但 EL 不能超过 (EM - EU) × C。若下层也耗尽剩余能力才进入深层深层蒸发 ED 等于剩余能力减去 EL同时受深层蓄水 WD 约束。这个三层规则的实际效果是上层蒸发不受限制下层和深层蒸发按蓄水比例衰减。C 参数一般取 0.1 到 0.2代表深层水向上补给的能力很弱。2.3 输入输出与时间步长日尺度的标准管道新安江模型最常用的时间步长是日也有用于次洪的三小时步长版本。日尺度模型的输入只有三个序列日降雨量 P毫米、日蒸发皿蒸发量 EP毫米、日实测流量 Q立方米每秒或毫米。输出是逐日模拟流量。这里有一个重要的单位换算问题降雨和蒸发都是毫米深度流量却是立方米每秒在计算水量平衡时要乘上流域面积并除以 86400 秒才能对齐。我一般先把流量也换算成毫米每天这样所有中间量都在同一套单位下调试时不容易乱。数据文件建议用 CSV 或 Excel列名固定为 date、P、EP、Q。日期列必须是标准格式否则后续怎么对齐都会出问题。预热期的问题也要在这个阶段想清楚模型初始张力水蓄量并不知道通常的做法是提前跑一年作为预热期把模拟结果的前一年丢弃用这一年的数据来消除初始蓄量的影响。3. 用 Python 实现三层蒸散发与蓄水更新核心函数逐行拆解3.1 参数初始化与状态变量把模型配置写成一个字典我不建议把参数散落在代码各处最稳的写法是在文件开头集中定义一个参数字典所有公式都从字典取值。这样后面率定时只需要改一处避免出现「改了 A 函数里的参数、B 函数还在用旧值」这种低级错误。下面这份字典覆盖了新安江模型三层蒸发和蓄满产流所需的核心参数。# 模型参数配置单位统一为 mm 和 日 params { KC: 1.0, # 蒸散发折算系数蒸发皿蒸发转换为流域蒸散发能力 WM: 120.0, # 流域平均张力水容量mm常取 80~150 WUM: 20.0, # 上层张力水容量mm常取 15~25 WLM: 80.0, # 下层张力水容量mm常取 60~90 WDM: 20.0, # 深层张力水容量mmWM - WUM - WLM B: 0.3, # 蓄水容量曲线指数反映空间不均匀性 C: 0.16, # 深层蒸散发系数常取 0.1~0.2 IM: 0.01, # 不透水面积比例 SM: 15.0, # 自由水蓄水容量mm EX: 1.5, # 自由水蓄水容量曲线指数 KG: 0.35, # 地下径流出流系数 KI: 0.35, # 壤中流出流系数 KKG: 0.98, # 地下径流消退系数 KKI: 0.90, # 壤中流消退系数 KKSS: 0.80, # 表层自由水消退系数汇流简化用 }参数初始化的逻辑说明上层容量 WUM 控制雨后地表附近的蓄水能力它在三层中最小因为表层土壤厚度有限。WM 减去 WUM 和 WLM 后得到 WDM深层水容量决定了干旱期基流的持续能力。注意 IM 是不透水面积比例这部分面积上降雨直接变成径流不参与蓄满产流计算。看到这你可能发现新安江模型的参数不少但实际需要手动率定的并不到一半很多参数可以从流域土壤类型和 land use 估算初值。状态变量我用一个字典或者 NumPy 数组保存每日的 WU、WL、WD、W、自由水蓄量 S逐日循环更新。初始时刻 WU、WL、WD 分别取各自容量的 0.5 倍左右即可反正预热期会把这部分误差消化掉。如果数据长度不够一年建议把初始蓄量设为 WUM、WLM 的一半WDM 按零算这样冷启动偏差会小一些。3.2 三层蒸散发计算函数严格按规则顺序执行蒸散发的计算顺序十分关键写错一个分支就会导致下层蒸发永远算不出来。下面这个函数接收当日蒸发能力 EM 和当前三层张力水蓄量返回三层蒸发量和更新后的蓄量。def evap_3layer(EM, WU, WL, WD, WUM, WLM, WDM, C): 三层蒸散发计算严格按上层→下层→深层的顺序消耗 # 上层蒸发优先满足但不超过上层蓄水量 EU min(EM, WU) WU WU - EU remain EM - EU # 剩余蒸发能力 # 下层蒸发按蓄水比例分配同时受 C 限制 EL 0.0 if remain 0 and WL 0: EL remain * (WL / WLM) EL min(EL, remain * C, WL) WL WL - EL remain remain - EL # 深层蒸发剩余能力还在时受 C 和 WD 双重限制 ED 0.0 if remain 0 and WD 0: ED min(remain, remain * C, WD) WD WD - ED remain remain - ED # 完全消耗不掉的部分直接舍弃流域蒸发能力未被充分利用 return EU, EL, ED, WU, WL, WD逻辑说明上层蒸发用 min(EM, WU)保证上层蓄水不够时不会强行蒸发到负值。下层蒸发的公式里remain * (WL / WLM) 表示下层蓄水越满蒸发越多符合土壤含水量越高、蒸发越接近潜在值的物理规律。乘 C 是限制深层补给能力防止下层长期干涸时还硬算蒸发。这三个 min 嵌套的顺序不能调换否则可能出现 EL 算完、remain 还没扣减的问题。参数说明C 在这个函数里的作用有两处一处限制 EL一处限制 ED。很多教材只写「深层蒸发量为 C 乘以剩余蒸发能力」但实际落地时你会发现若下层蓄水接近 WLMEL 会先吃掉大部分剩余能力ED 自然变小。这是合理的下层有水时优先下层蒸发深层只作最后补充。3.3 蓄满产流与张力水更新逐日主循环的实现三层蒸发更新完张力水蓄量后紧接着做蓄满产流计算。这里有一个容易忽略的细节降雨到达地面后要先扣除雨期蒸发剩余部分 PE 才参与产流。实际程序中我把当日降雨 P 与三层蒸发总量做先后处理蒸发先消耗蓄水降水再补给蓄水然后判断是否产流。def runoff_generate(PE, W, WM, B, IM): 蓄满产流计算返回产流量 R 和更新后的张力水蓄量 W # 不透水面积直接产流 R PE * IM # 透水面积按蓄满产流公式计算 WMM WM * (1 B) AU WMM * (1 - (1 - W / WM) ** (1 / (1 B))) if PE AU WMM: R R PE - (WM - W) WM * ((1 - (1 - AU / WMM) ** (1 B))) else: R R PE - (WM - W) # 更新透水面积张力水蓄量 W W PE - (R - PE * IM) if W WM: # 数值保护理论上 W 不会超过 WM超过说明上游计算有误差 R R (W - WM) W WM return R, W逻辑说明不透水面积 IM 的产流是直接加在总产流量上的这些面积没有蓄水能力降雨即产流。透水面积部分先根据当前 W 和蓄水容量曲线求 AU再判断 PE 是否超过剩余蓄水能力。当 PE 很大时公式自动落到全流域蓄满的情形R PE - (WM - W)物理含义是降雨扣除补齐缺水量后的剩余全部变成径流。最后的数值保护非常必要浮点误差累积可能导致 W 略微超过 WM此时把多余水量当作径流确保蓄量状态始终在物理可行的范围内。参数说明B 越大AU 的计算结果越大意味着同一蓄水量下达到蓄满所需降雨越小。这从物理上也说得通B 大说明流域内蓄水容量小的区域占比高这些地方先蓄满先产流。率定时一般先定 WM再调 BB 对洪峰影响明显对基流影响较小。主循环的完整骨架如下它把前面两个函数串成逐日迭代import pandas as pd def run_xinanjiang(df, params): # 初始化状态变量 WU params[WUM] * 0.5 WL params[WLM] * 0.5 WD 0.0 W WU WL WD S 0.0 # 自由水蓄量初始为 0 qsim_list [] for i in range(len(df)): P df.loc[i, P] EP df.loc[i, EP] # 蒸散发能力 EM params[KC] * EP EU, EL, ED, WU, WL, WD evap_3layer( EM, WU, WL, WD, params[WUM], params[WLM], params[WDM], params[C] ) W WU WL WD # 净雨 PE P - (EU EL ED) if PE 0: PE 0.0 # 蓄满产流 R, W runoff_generate(PE, W, params[WM], params[B], params[IM]) # 自由水蓄水库与汇流计算下一章展开 q routing(R, S, params) # 重新分配三层蓄量按各自的容量上限限制 WU min(WU, params[WUM]) WL min(WL, params[WLM]) W WU WL WD qsim_list.append(q) df[Qsim] qsim_list return df逻辑说明每个时间步内先算蒸发再算净雨再算产流最后做汇流。PE 小于零时强制归零因为蒸发已经扣除降雨为负没有物理意义。三层蓄量更新后要重新按容量上限约束一次防止浮点误差把 WU 推到 WUM 之上。这种防御式写法对长序列尤其重要连续模拟十年以上时微小误差会不断累积。4. 产流与汇流的完整链路从自由水蓄水库到逐日流量过程4.1 自由水蓄水库三水源划分的关键中间层蓄满产流算出来的 R 是总径流量还要把它划分成地表径流、壤中流和地下径流否则后续汇流没法体现退水过程的几个不同衰减段。新安江模型在这里引入了一个自由水蓄水库蓄量 S容量 SM。当 R 补入自由水蓄水库后超过 SM 的部分直接成为地表径流 RS蓄水库内的水按出流系数 KI 和 KG 分别生成壤中流 RI 和地下径流 RG。划分的规则是先按自由水蓄水容量曲线指数 EX计算产流比例 FR再把 R 分配到 S 中若 S 超过 SM溢出部分为 RS。随后[ RI KI \times S ][ RG KG \times S ][ S S - RI - RG ]注意 KI 和 KG 之和一般小于 1多出来的部分留在蓄水库中供下一步继续消退。这两个系数决定了壤中流和地下径流的出流速度直接影响退水曲线的陡缓。代码实现时我把自由水蓄水库和汇流放在同一个函数里因为它们在时间步进上是耦合的。def routing(R, S, params): 自由水蓄水库三水源划分 线性水库汇流 返回当日模拟流量 q_mm SM params[SM] EX params[EX] KI params[KI] KG params[KG] KKSS params[KKSS] KKI params[KKI] KKG params[KKG] # 1. 计算产流面积比例自由水蓄量越满产流面积越大 if S SM: FR 1.0 else: FR (S / SM) ** EX # 2. 蓄水库补水 S S R / max(FR, 0.01) # 3. 超过容量部分为地表径流 RS 0.0 if S SM: RS S - SM S SM # 4. 划分壤中流和地下径流 RI KI * S RG KG * S S S - RI - RG # 5. 线性水库消退简化汇流 # 实际程序中这些变量应作为持久状态跨日保存 return RS RI RG逻辑说明这版 routing 函数做了汇流简化处理把当天 RS、RI、RG 直接加起来作为当日径流输出。如果你想看到更平滑的退水过程需要对 RI 和 RG 分别做消退演算即# 线性水库消退的写法示例 q_gi q_gi_prev * KKG RG * (1 - KKG) q_ii q_ii_prev * KKI RI * (1 - KKI)地表径流 RS 一般用单位线或滞后演算法推演在日尺度模型中可以直接用当日 RS 代表当日直接径流误差可以接受。KKG 取 0.98 左右时地下径流的消退非常缓慢这正是湿润地区基流长退水段的特征KKI 取 0.9壤中流退水稍快。这几个消退系数在率定阶段是调节洪峰后退水形态的主要旋钮。4.2 汇流演算的几种常见处理从滞后演算到单位线日尺度新安江模型的汇流部分不需要做得很重常见的做法是地表径流走滞后演算法壤中流和地下径流走线性水库消退。滞后演算法很简单把当天的 RS 延迟 L 天输出即可L 是流域汇流时间参数一般取 1 到 3 天。这样做的好处是能模拟出洪峰滞后于降雨的现象而线性水库消退则刻画退水段的指数衰减。如果你要模拟场次洪水而不仅仅是日径流序列建议把时间步长缩到 1 小时或 3 小时汇流部分改用瞬时单位线或地貌单位线。日尺度模型对汇流参数不敏感率定出来的 L 换到次洪模型往往要重率定这是正常的不必追求一套参数通吃所有时间尺度。4.3 核心参数表量纲、推荐范围与调参优先级下面这张参数表是我实际调参时最常参照的量纲统一为毫米和日。优先级列代表率定顺序的建议先定优先级高的参数再调低的。参数含义常用范围优先级调参影响KC蒸散发折算系数0.8 ~ 1.2高影响总水量平衡WM平均张力水容量80 ~ 150 mm高影响产流阈值控制径流总量B蓄水容量曲线指数0.1 ~ 0.5中影响洪峰大小与产流时机C深层蒸散发系数0.1 ~ 0.2低影响枯水期蒸发和基流WUM上层张力水容量15 ~ 25 mm低影响雨后短时蒸发WLM下层张力水容量60 ~ 90 mm低影响中层蒸发过程SM自由水蓄水容量5 ~ 30 mm高影响地表/地下径流比例EX自由水蓄水曲线指数1.0 ~ 2.0中影响洪峰削减KI / KG壤中流/地下流出流系数0.3 ~ 0.5中影响退水速度KKG / KKI消退系数0.8 ~ 0.99中影响基流退水段参数说明KC 是最先要定的它直接控制流域蒸散发总量率定目标通常是让模拟年径流量和实测年径流量接近。如果模拟流量整体偏大先降 KC 或升 WM两者效果相似但物理含义不同。SM 的大小决定了降雨产流中有多少变成地表径流SM 越大蓄水库越不容易漫溢地表径流占比越低洪峰越平缓。这套参数之间的耦合关系是调参中最容易翻车的地方后面第 5 章会具体展开。5. 跑通代码后的四个常见问题数据、单位与参数冷启动的避坑记录5.1 时间序列错位模拟流量比实测流量整体滞后一天现象程序跑完NSE 看起来还行但逐日对比发现模拟洪峰总比实测晚一天退水段也整体偏移。原因数据文件里的降雨日期和流量日期不是同一套时间口径。雨量站记录的是前一天的降雨量而蒸发和流量是当天 8 点的观测值这在很多水文资料里是常态。第二个可能原因是代码里用了pd.shift()但对齐方向搞反了。解决先把三条序列画在同一张图上人工检查。然后在代码里统一日期口径比如明确「当日降雨对应当日 8 时至次日 8 时」把降雨序列整体前移或后移一天再对齐。我通常的做法是写一个小的对齐函数比较降雨与流量的互相关自动找出最优滞后天数但更稳的做法还是直接用肉眼核对原始数据说明。5.2 蒸散发输入错误直接用了蒸发皿蒸发却没乘 KC现象模拟径流量系统性偏小尤其在非汛期模拟流量长期贴近零。原因数据提供的蒸发值通常是大口蒸发皿蒸发量而新安江模型需要的是流域蒸散发能力两者之间差一个折算系数 KC。直接拿蒸发皿蒸发当 EM 用相当于把蒸散发能力放大了土壤蓄水被过度消耗产流自然偏少。解决在数据预处理阶段把 EP 乘以 KC 再进模型或者把 KC 参数调小到 0.6 左右来补偿。我建议前者因为 KC 的取值应该反映蒸发皿到流域的折算关系而不是作为误差吸收器。最容易踩坑的地方是同时用了自动率定工具率定算法会把 KC 和 WM 一起调导致参数失衡——这时候要固定 WM 或 KC 中的一个先单独率定另一个。5.3 冷启动模拟 NSE 为负初始蓄量误差吞噬了前段模拟结果现象模拟序列前三个月完全对不上实测NSE 为负但把前半年丢掉后评价指标明显变好。原因初始张力水蓄量 WU、WL、WD 设置不合理。流域前期偏干但初始蓄量设大了模型一上来就大量产流和实测的枯水段对不上。解决预留预热期。把模拟起点往前推一年跑完后丢弃前 365 天的结果用后面的数据计算评价指标。更讲究一点的做法是让模型先单独跑两遍第一遍用默认初始蓄量第二遍用第一遍结束时的蓄量作为初始条件这样预热效果更好代价只是多跑一遍循环而已。5.4 张力水蓄量溢出与负值浮点误差累积导致状态量越界现象程序跑长序列时报错或者某个中间状态变量 W 变成负值、WU 超过 WUM产流量出现异常大值。原因日尺度循环里反复做乘方和减法浮点误差逐日累积加上 PE 特别大时产流公式的中间结果可能接近机器精度边界。解决在每次状态更新后做钳制W 不超过 WM、WU 不超过 WUM、WD 不为负。这类防御代码虽然看起来多余但跑十年日序列时完全不是多余的。另外所有涉及到(1 - W/WM)的表达式要先判断 W/WM 是否在 0 到 1 之间越界时强制赋 0 或 1避免对负数开根号。6. 率定参数的顺序与一个硬指标怎么判断这份代码值不值得继续调参数率定我强烈建议手动优先自动率定放第二位。原因是新安江模型的参数之间耦合很强自动算法容易把参数推到物理上不合理但数值上最优的角落。我一般按水量平衡、洪峰过程、退水过程三步走。第一步固定 B、C、WUM 等次要参数只调 KC 和 WM 让模拟年径流总量和实测误差小于 5%。第二步调 SM 和 KI、KG让洪峰大小和涨水速率对上。第三步调 KKG、KKI 和汇流滞后 L让退水段的形态贴近实测。每一步都画三张图实测对模拟的全年过程线、汛期放大过程线、退水段半对数图半对数图最能暴露基流消退系数的问题。验证指标只看一个还不够我会同时输出 NSE、水量平衡误差和洪峰相对误差。NSE 反映整体拟合度但它在低流量段不敏感可能全年过程线看起来不错、洪峰却明显偏低。洪峰相对误差单独看控制在正负 20% 以内才算合格。水量平衡误差是硬约束超过 5% 说明蒸散发或产流总量的分配有问题后面再怎么调峰形都是白费力气。我自己第一次做这个模型时就是先自动率定了一轮参数确实好看NSE 到了 0.8 以上但换到另一年的验证期直接掉到 0.2连洪峰方向都对不上。后来老老实实按顺序手动调虽然花了两天时间但验证期稳定在 0.75 以上。这段血泪经验让我形成了一个固定习惯每次改参数后必须重新检查水量平衡否则模型的模拟结果就只是个黑匣子永远不知道它为什么好、为什么坏。这套 Python 计算程序本身并不复杂复杂的是对每个参数物理含义的尊重。希望这份笔记能帮你少走一段弯路跑通它然后真正理解它。本文还有配套的精品资源点击获取