企业官网建设流程全解析

简介：面向水文水资源领域研究与应用的三层蒸发蓄满产流模型（新安江模型）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状态更新

用伪代码描述这个逻辑：

1. 上层蒸散发量EU = min(WU, EM)，然后WU减去EU。
2. 如果EU小于EM，说明上层已经蒸干，剩余蒸发能力EM_left = EM - EU。
3. 下层需求为EM_left * WL / WLM，实际下层蒸散发量EL = min(WL, 需求)，然后WL减去EL。
4. 如果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{W'm}{W{MM}}\right)^B ]

其中W'_m是单点蓄水容量，WMM = WM * (1 + B)是最大点蓄水容量，B是曲线指数，一般取0.2~0.4。B越大，说明流域蓄水容量空间差异越大。

编程时不能直接拿这个公式做产流，因为不知道当前时刻各点的蓄水状态。标准做法是：先根据当前张力水总蓄量W0 = WU + WL + WD，反解一个虚拟纵坐标A：

A = WMM * (1 - (1 - W0 / WM) ** (1 / (1 + B)))

A的物理含义是：在当前流域蓄水状态下，蓄水容量曲线横轴上对应的分界蓄量。净雨PE到来后，所有蓄水容量小于A+PE的点都已经蓄满，这些点面积上产生的径流就是总产流量。

3.2 产流量两种情形的统一处理

产流量的计算分两种情况。如果A + PE < WMM，说明还有部分面积未蓄满：

[ R = PE - WM \left[ \left(1 - \frac{A}{W_{MM}}\right)^{1+B} - \left(1 - \frac{A+PE}{W_{MM}}\right)^{1+B} \right] ]

如果A + PE ≥ WMM，说明全流域蓄满，此时R = PE - (WM - W0)，简单直接。

这里要注意，A在每次产流后必须更新，因为张力水蓄量变了。更新公式就是W0_new = W0 + PE - R，即降雨中扣除蒸发和产流后剩余的，全部蓄在张力水里。实际代码里，我把PE减掉的这部分补充到WU里，这样三层张力水容量的约束自动保持。

3.3 自由水蓄水容量曲线与三水源分配

产流量R不是直接变成地表径流，而是先进入自由水蓄量S。自由水蓄水容量SM同样用抛物线分布描述，指数为EX。计算思路和张力水类似：先根据当前S反解AU：

AU = SMM * (1 - (1 - S / SM) ** (1 / (1 + EX)))

SMM = SM * (1 + EX)。然后计算自由水蓄量增量DS和地面径流RS：

• 若AU + R < SMM：DS = R - SM * [(1-AU/SMM)^(1+EX) - (1-(AU+R)/SMM)^(1+EX)]
• 若AU + R ≥ SMM：DS = 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必须大于WUM+WLM，SM必须大于0，C在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.QG

CS、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这些状态变量打印出来，亲眼看着它们如何随降雨和蒸发变化，比只看最终流量过程线有用得多。

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