简介面向水文水资源领域研究与应用的三层蒸发蓄满产流模型新安江模型Python计算程序用于基于降水、蒸发能力数据驱动完成产流过程模拟。程序主体为单个.py脚本压缩包内共1个文件体积约1KB轻量简洁适合水文专业学生、科研人员及模型初学者阅读与调试。运行时需自行准备或从作者处获取包含降水、蒸发能力数据的PEdata.xlsx文件以正确驱动模型计算。目前已有2642人学习下载可见其在水文建模入门与教学场景中具有一定参考价值。读者通过该源码可快速掌握三层蒸发与蓄满产流的核心算法结构并能够依据个人数据调整参数、扩展功能是学习新安江模型Python化实现的实用工具。 水文预报圈子里有个老问题新安江模型都用了快五十年现成工具也多为什么还有人要自己写Python程序我的答案是软件包帮你算结果但不会帮你理解流域。三层蒸发蓄满产流模型新安江模型的Python实现看似是写一段水文计算代码实际上是把包气带的蓄水、蒸发、产流、汇流机制从头到尾复述一遍。这篇文章我想分享一份可运行的教学版实现思路和核心代码适合刚接触水文模拟的年轻工程师也适合那些想把自己单位旧平台上的预报方案迁移到Python环境的同学。程序不长但涉及的知识点非常密集我会把关键公式的来龙去脉、代码的组织方式、以及我实际调试时踩过的坑都讲清楚。1. 新安江模型的分块逻辑先分清蓄和流两本账1.1 三套蓄水容器与三股径流来源新安江模型本质上是一套水量账本。整个流域在垂向上被划分为三层张力水蓄水容器上层WU、下层WL和深层WD分别代表植被截留层与表层土、根系层、深层包气带。三者的蓄水容量差异很大典型取值是上层WUM为5~20毫米下层WLM为60~90毫米深层则为总张力水容量WM减去前两者。产流阶段还有第二套容器自由水蓄量S。张力水蓄满之后多余的水进入自由水层再按比例分成三股径流——地表径流RS、壤中流RI和地下径流RG。之所以要分三股是因为它们的汇流速度差别太大地表水几天内就能到断面壤中流持续几周基流则能维持几个月。如果只算总径流退水段永远模拟不像。1.2 模型主流程先蒸散发、再产流、最后汇流每个时段的计算顺序是固定的先根据蒸发能力EM和当前土壤含水量扣减蒸散发更新三层张力水然后计算净雨PE用蓄水容量曲线算产流量R接着把R引入自由水蓄量按自由水蓄水容量曲线分水源最后三股水分别经过线性水库调蓄后叠加得到断面流量。这个顺序不是随意定的它反映的是物理过程的先后蒸发发生在降雨之前还是之后会影响土壤含水量的初始状态进而影响产流量。实际程序中我把蒸发放在时段最前面这样净雨PE P - EM只有降雨大于蒸发能力时才可能产流。2. 三层蒸散发逐层清算代码里最容易出错的环节2.1 蒸发能力折算与三层容量的含义蒸发数据通常来自蒸发皿观测需要乘折算系数K得到流域蒸发能力EM。K在不同流域差异很大湿润地区一般在0.8~1.0干旱半干旱地区可能降到0.5~0.7。这个系数直接决定水量平衡的闭合程度是第一个要率定的参数。三层蒸散发的核心思想是逐层剥夺上层含水量少但蒸发不受限制只要有水就优先蒸发上层蒸干后下层按土壤含水量占其容量的比例蒸发只有当下层也满足不了剩余蒸发能力时才考虑动深层的水而且深层蒸发受系数C控制。2.2 逐层扣减的顺序与WU、WL、WD状态更新用伪代码描述这个逻辑上层蒸散发量EU min(WU, EM)然后WU减去EU。如果EU小于EM说明上层已经蒸干剩余蒸发能力EM_left EM - EU。下层需求为EM_left * WL / WLM实际下层蒸散发量EL min(WL, 需求)然后WL减去EL。如果EL小于需求说明下层也蒸干此时才判断是否启动深层蒸发。启动条件是下层含水量与下层容量之比小于C由于下层已干这个条件通常满足深层蒸散发量ED min(C * EM_left_remaining, WD)然后WD减去ED。我见过很多实现把深层蒸散发的触发条件写错有的不判断下层是否蒸干就直接从深层扣水有的把C乘到了总蒸发能力上。这两种写法都会导致深层水被过度消耗长期运行后基流明显偏小退水过程过于干瘪。2.3 一个可以复制的蒸散发子函数def evap(self, em): 三层蒸散发输入em为当日蒸发能力(mm) eu min(self.WU, em) self.WU - eu remain em - eu el 0.0 ed 0.0 if remain 0: demand remain * self.WL / self.WLM el min(self.WL, demand) self.WL - el remain - el # 只有下层蒸干且仍有剩余蒸发能力时才启用深层 if el demand - 1e-6 and remain 0: ed min(self.C * remain, self.WD) self.WD - ed return eu, el, ed注意最后一行的min(ed, WD)很多人会漏掉深层水容量是有限的不能无限蒸发。另外浮点比较要留容差否则在边界状态下容易进错分支。3. 蓄满产流与三水源划分蓄水容量曲线如何驱动径流3.1 张力水蓄水容量曲线反解初始蓄量蓄满产流的前提假设是流域内某点包气带蓄满之前不产流蓄满之后降雨全部产流。但流域内各处蓄水容量并不均匀赵人俊团队用一条抛物线来描述蓄水容量的空间分布[ \frac{f}{F} 1 - \left(1 - \frac{Wm}{W{MM}}\right)^B ]其中W_m是单点蓄水容量WMM WM * (1 B)是最大点蓄水容量B是曲线指数一般取0.2~0.4。B越大说明流域蓄水容量空间差异越大。编程时不能直接拿这个公式做产流因为不知道当前时刻各点的蓄水状态。标准做法是先根据当前张力水总蓄量W0 WU WL WD反解一个虚拟纵坐标AA WMM * (1 - (1 - W0 / WM) ** (1 / (1 B)))A的物理含义是在当前流域蓄水状态下蓄水容量曲线横轴上对应的分界蓄量。净雨PE到来后所有蓄水容量小于APE的点都已经蓄满这些点面积上产生的径流就是总产流量。3.2 产流量两种情形的统一处理产流量的计算分两种情况。如果A PE WMM说明还有部分面积未蓄满[ R PE - WM \left[ \left(1 - \frac{A}{W_{MM}}\right)^{1B} - \left(1 - \frac{APE}{W_{MM}}\right)^{1B} \right] ]如果A PE ≥ WMM说明全流域蓄满此时R PE - (WM - W0)简单直接。这里要注意A在每次产流后必须更新因为张力水蓄量变了。更新公式就是W0_new W0 PE - R即降雨中扣除蒸发和产流后剩余的全部蓄在张力水里。实际代码里我把PE减掉的这部分补充到WU里这样三层张力水容量的约束自动保持。3.3 自由水蓄水容量曲线与三水源分配产流量R不是直接变成地表径流而是先进入自由水蓄量S。自由水蓄水容量SM同样用抛物线分布描述指数为EX。计算思路和张力水类似先根据当前S反解AUAU SMM * (1 - (1 - S / SM) ** (1 / (1 EX)))SMM SM * (1 EX)。然后计算自由水蓄量增量DS和地面径流RS若AU R SMMDS R - SM * [(1-AU/SMM)^(1EX) - (1-(AUR)/SMM)^(1EX)]若AU R ≥ SMMDS R - (SM - S)地面径流RS R - DS即自由水未蓄住的那部分水量。随后壤中流和地下径流按各自出流系数从S中出流ri KI * self.S rg KG * self.S self.S - (ri rg)KI、KG一般合起来取0.5~0.8两者之比决定壤中流和基流的分配。SM则是影响产流面积和洪峰形态的关键参数取值通常在5~50毫米之间SM偏大会让产流面积变小、洪峰变缓。4. Python实现全流程一个可运行的模型类4.1 参数集与状态变量的组织方式我建议把所有参数收进一个字典或配置类状态变量作为模型实例的属性。这样进行多方案对比时只需要复制实例并替换参数即可。模型类大致骨架如下class XajModel: def __init__(self, params): self.load_params(params) self.reset_state()参数加载时注意做合法性检查WM必须大于WUMWLMSM必须大于0C在0~1之间。这些检查能省掉后面一大堆莫名其妙的结果。4.2 蒸散发、产流、分水源、汇流四个核心方法蒸散发方法在第2节已经给出。产流和分水源方法对应第3节公式。汇流部分最直接的实现是三个线性水库def route(self, rs, ri, rg): self.QS self.CS * self.QS (1 - self.CS) * rs self.QI self.CI * self.QI (1 - self.CI) * ri self.QG self.CG * self.QG (1 - self.CG) * rg return self.QS self.QI self.QGCS、CI、CG是三个消退系数取值在0~1之间越接近1表示调蓄能力越强、退水越慢。地表水库的CS一般取0.1~0.4反映地表径流快速退水壤中流CI取0.5~0.9地下水CG取0.9~0.999反映基流缓慢消退。4.3 逐时段驱动主循环与水量平衡检查主循环很简单for p, e in zip(rain, evap): em self.KC * e eu, el, ed self.evap(em) pe max(p - em, 0.0) r self.gen_runoff(pe) rs, ri, rg self.split_runoff(r) q self.route(rs, ri, rg) qs.append(q)循环末尾强烈建议做水量平衡检查累计降雨减去累计蒸散发、产流、出流和土壤蓄水变化误差应小于0.1毫米量级。我实际调试中发现很多bug比如蒸散发重复扣减、产流后没有更新WU都能被这道检查当场抓出来。新安江模型虽然概念简单但状态变量之间的耦合关系非常容易出错。4.4 初始状态与预热期处理冷启动时三层张力水和自由水蓄量都设为0但实际流域在雨季前包气带往往有前期含水量。直接把初始蓄量设0模拟初期会有一段干土层吸水过程导致前几个月的模拟流量系统性偏小。解决办法有两个一是用实测前期影响雨量或前期径流估算初始蓄量二是模型前面加一年预热期用历史数据跑一遍后再开始统计效率系数。第二个办法更省事我在代码里预留了一个spinup参数默认365天。5. 调试与率定中的实战经验哪些参数先调哪些坑先躲5.1 参数敏感性排序先水量后过程我自己的经验是参数敏感性从高到低大致是WM、SM、K、KG、KI、B、CS、CI、WUM、WLM、C、EX。率定时不要一上来就自动优化先把最不敏感的参数定下来再一步步缩小范围。第一步是调水量平衡K和各层容量决定年总蒸散发和总径流量K偏大一年下来模拟径流偏小K偏小则水量盈余。第二步看过程形态SM和B控制产流面积动态直接影响洪峰大小和涨水段形态。第三步才是汇流参数CS、CI、KG它们控制洪峰滞后和退水段形态。5.2 五个容易翻车的细节第一是量纲。蒸发能力EM的单位必须是毫米/时段跟降雨一致。如果蒸发资料是逐日、降雨是逐小时必须先把蒸发换算成小时尺度否则蒸散发量级完全不对。第二是负值处理。PE P - EM可能为负必须截图成0。同样自由水蓄量在扣除RI和RG后也可能出现微小负值加一个max(0)兜底。第三是深层蒸散发启动条件。很多版本不判断下层是否蒸干就直接算深层蒸发模型在湿润期会把深层水白白蒸发掉基流模拟严重偏低。第四是S的上下界。自由水蓄量S严格来说不会超过SM但在分水源公式的近似表达中可能出现S略大于SM的情况需要钳制。同时蓄水容量曲线公式中S/SM不能等于1否则幂运算分母为0。第五是预热期不够。新安江模型每层蓄量都有记忆初始状态错一点后面几十天都会受影响。至少跑够一个水文年再开始统计模拟序列前面一段直接丢弃。5.3 自动率定的目标函数建议如果用SCE-UA或遗传算法做自动率定目标函数不要只盯NSE。NSE对洪峰峰值敏感、对基流过程不敏感容易出现过拟合洪峰而年径流总量偏小的情况。推荐用复合目标函数[ F (1 - NSE) \alpha \cdot |\text{水量平衡误差}| ]水量平衡误差是模拟总径流与实测总径流的相对偏差通常控制到5%以内。还需要给参数加物理约束比如KG必须大于0.9、CS必须小于0.5否则优化算法为了极小化目标函数可能跑出物理上毫无意义的结果。结尾一点建议我在实际调试这个程序时最大的体会是新安江模型的代码并不难难的是把每个公式的前提条件搞清楚。蓄水容量曲线不是拿来就算的先反解A再用A推产流这个中间环节很多人会跳过去结果算出来的产流量要么偏大要么为负。如果初学者想把这份代码作为基础扩展我建议下一步加入产流面积FR的动态计算即把集总式模型改造为考虑部分产流面积的分布式结构这会让洪峰模拟效果上一个台阶。另外调试时一定把逐时段的WU、WL、WD、S这些状态变量打印出来亲眼看着它们如何随降雨和蒸发变化比只看最终流量过程线有用得多。本文还有配套的精品资源点击获取