企业官网建设流程全解析

简介：这是一份面向水文水资源领域的三层蒸发蓄满产流模型（新安江模型）Python计算程序，适用于水利工程、水文预报及流域模拟方向的从业者、研究人员与学生，帮助解决蓄满产流过程的核心计算与编程实现问题。压缩包内仅含1个py脚本文件，整体约1KB，程序结构紧凑、逻辑清晰，适合直接阅读与二次修改。运行程序需同步使用配套的PEdata.xlsx数据文件，其中包含计算蓄满产流所需的降水、蒸发能力等输入数据。通过本程序可系统理解新安江模型三层蒸发模式的计算流程、参数衔接关系及Python代码组织方式，既可当作水文模型课程的教学参考，也可作为科研中快速搭建产流计算模块的基础脚本。该资源已有2642人浏览学习，适合需要上手新安江模型编程或开展流域水文模拟的读者下载使用。

1. 三层蒸发蓄满产流模型：这份 Python 代码为什么值得你认真跑一遍

做水文模拟的人迟早会撞上一个尴尬局面：拿着现成的降雨、蒸发、流量序列，想算一个流域的产流量，却不知道该选哪种模型。用超渗产流算湿润地区，洪水过程畸大；用统计回归算旬月尺度，又没法解释日尺度的涨水退水。而标题里的「三层蒸发蓄满产流模型」，也就是大家常说的新安江模型，恰恰是湿润半湿润地区日尺度连续模拟最稳的一条路。它把流域看成一块会蓄水的海绵，按张力水蓄量判断是否产流，再按上、下、深三层分配蒸散发，最终输出逐日径流过程。这篇文章不讲模型史，直接讲怎么用 Python 把这套逻辑落地成可复现的程序，适合正在做水文模拟、洪水预报或水资源评估的工程师和研究生。你会看到完整的参数表、逐行解释的代码，以及我最开始跑这个程序时踩过的几个隐蔽的坑。

2. 蓄满产流与三层蒸发的物理逻辑：先立住新安江模型的四个概念

2.1 蓄满产流的本质：张力水蓄量与流域蓄水容量曲线

新安江模型的核心假设是蓄满产流，意思是只有当某点的土壤张力水蓄量达到田间持水量（即蓄满）后，后续降雨才会产生径流。这个假设在湿润地区非常贴近实际，因为湿润地区的土壤含水量常年偏高，降雨稍大就容易达到蓄满状态，地表径流和壤中流的生成主要由蓄满后的多余水量贡献。

要描述整个流域并非同时蓄满，模型引入了一条蓄水容量曲线。常用的是指数型或抛物线型分布，表示流域内各点张力水容量 Wm' 的空间不均匀性。曲线形态由参数 B 控制，B=0 意味着全流域蓄水容量均匀，B 越大表示空间差异越大。初始时刻流域实际张力水蓄量 W 对应的最大点蓄水容量为 AU，我一般把这个过程叫做「先找阈值，再算产流」。实际的产流公式如下。

[ WMM = WM \times (1 + B) ]

[ AU = WMM \times \left(1 - \left(1 - \frac{W}{WM}\right)^{\frac{1}{1+B}}\right) ]

当降雨扣除蒸发后的净雨 PE 加上 AU 仍小于 WMM 时，产流量 R 为：

[ R = PE - (WM - W) + WM \times \left(1 - \left(1 - \frac{AU + PE}{WMM}\right)^{1+B}\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, # 深层张力水容量，mm，WM - 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 乘以剩余蒸发能力」，但实际落地时你会发现若下层蓄水接近 WLM，EL 会先吃掉大部分剩余能力，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，再调 B，B 对洪峰影响明显，对基流影响较小。

主循环的完整骨架如下，它把前面两个函数串成逐日迭代：

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 是最先要定的，它直接控制流域蒸散发总量，率定目标通常是让模拟年径流量和实测年径流量接近。如果模拟流量整体偏大，先降 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 计算程序本身并不复杂，复杂的是对每个参数物理含义的尊重。希望这份笔记能帮你少走一段弯路，跑通它，然后真正理解它。

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