简介这份资源是一份面向水文水资源专业学生与从业者的三层蒸发蓄满产流模型新安江模型Python计算程序用于完成基于降水与蒸发能力数据的产流模拟计算。程序核心实现新安江模型中的三层蒸发与蓄满产流机制适合学习流域水文模型、开展课程设计或进行科研验证的初级至中级Python使用者。资源包为1个py文件体积约1KB内容紧凑主程序需配合作者提供的PEdata.xlsx降水蒸发数据运行可直接复制使用或在此基础上修改参数。截至目前已有2642人浏览学习说明该程序在同类资源中具有较高实用参考价值。通过运行这份代码读者可以直观理解新安江模型的计算流程快速上手实现产流模拟为后续开展流域水文分析提供可复用的脚本基础。1. 三层蒸发蓄满产流模型新安江模型python计算程序.py从降雨到径流的确定性推演做水文预报的人迟早会碰到这个场景手里有流域的逐日降雨和蒸发数据想算出对应的径流过程但商业化软件黑匣子太多参数一改就不知道内部发生了什么。三层蒸发蓄满产流模型常说的新安江模型就是干这个事的——它把流域降水转化成蒸散发、土壤蓄水、地表径流和壤中流是一个完整闭合的水量平衡计算框架。标题里这个python计算程序.py本质是一份可以直接运行的脚本解决的是没有现成水文模型库时如何用几百行代码实现经典蓄满产流计算的现实问题。适合水文专业学生做课程设计、科研人员做对比实验以及工程师做中小流域的洪水预报方案初筛。只要你有逐日降雨、蒸发和流量数据这套代码就能从零开始推演出模拟径流过程还能人工调参逼近实测值。核心并不神秘蓄满产流的核心就是先算土壤缺水量再决定产不产流。2. 蓄满产流与三层蒸发模型运行的底层逻辑2.1 张力水和自由水蓄满产流的核心区分蓄满产流模型能成立是因为它区分了两种性质完全不同的土壤水。张力水是被土粒吸附、能被植物根系吸收的水分它不参与重力排水自由水则是填充在土壤孔隙中、受重力支配可以横向或纵向流动的水分。降雨先补充张力水张力水蓄满后才开始产生自由水自由水超过一定容量就形成径流。这个先蓄满再产流的机制是模型最关键的判断条件。实际计算中每个计算单元都有一个张力水蓄水容量WM它等于上层张力水容量WUM、下层张力水容量WLM和深层张力水容量WDM之和。降雨先补上层上层满了溢到下层下层满了溢到深层。只有当整个土柱的张力水都蓄满降雨才转为自由水蓄量进而划分产流。这个顺序保证了水量平衡的闭合性也让参数物理意义明确——WUM和WLM可以根据土壤类型和经验值初定再通过率定微调。2.2 三层蒸发为什么不能只算一层很多初学模型的人会问一个流域的总蒸发量乘以一个系数不就完了为什么要分三层答案在于干旱期和湿润期的蒸发规律截然不同。湿润期土壤含水充足蒸发接近蒸散发能力干旱期表层先干蒸发衰减很快但深层仍有水分缓慢供给。三层蒸发结构正是为了模拟这种表层快、深层慢的衰减过程。模型里每层有自己的蒸发能力分配系数。上层蒸发EU取蒸散发能力EM乘以系数K如果上层有足够水量就直接消耗不够时下层蒸发EL取剩余蒸散发能力乘以一个折减系数但最多不能超过下层张力水蓄量的某个比例再不够才动用深层ED深层蒸发速率接近常数。这一套规则写进代码就三个条件分支却决定了连续干旱期的模拟精度。2.3 从降雨到径流水源划分和汇流张力水全部蓄满后超出的降雨进入自由水蓄量SM。自由水蓄量有一个最大的容量超过了就产地面径流RS同时自由水还按一个固定的消退系数KG向下补给深层地下水形成地下径流RG。壤中流RI则和自由水蓄量成正比计算为KI乘以自由水蓄量其中KI和KG是消退系数。降雨产流过程结束后各水源按不同路径进入河道地面径流几乎即时汇入壤中流衰减稍慢地下径流衰减最慢。模型通常用一个线性水库或单位线把三种水源分别演算到流域出口。这个先分水源、再分路径的做法让模型既能模拟洪水过程的陡涨陡落又能模拟退水过程的缓降比单一水库模型灵活得多。3. 编写新安江模型Python计算程序蓄满产流演算模块的完整实现3.1 数据输入和初始化直接读取Excel更方便我一般用pandas读取Excel里的三列数据日期、降雨量、实测流量。蒸散发能力可以单独一列也可以由程序根据月份自动调整。初始化时先定义土壤含水量初值、蓄水容量参数、蒸发系数以及各水源的消退系数。代码定义得越清晰后面调参越容易。import pandas as pd import numpy as np # 读取数据 df pd.read_excel(流域日数据.xlsx, sheet_nameSheet1) dates df[日期].values P df[降雨].values.astype(float) # 日降雨单位mm EM_input df[蒸发].values.astype(float) # 日蒸散发能力单位mm Q_obs df[流量].values.astype(float) # 实测流量单位m3/s用于对比 # 模型参数初始值 WUM, WLM, WDM 20, 90, 40 # 上、下、深层张力水容量(mm) WM WUM WLM WDM # 总张力水容量(mm) K 1.0 # 蒸散发能力折算系数 SM 15.0 # 自由水蓄水容量(mm) KG, KI 0.3, 0.35 # 地下水和壤中流出流系数 IM 0.1 # 不透水面积比例 FC 10.0 # 稳定下渗率(mm/day) # 注原文此处为帮助记忆的占位符 # 初值 WU, WL, WD 10.0, 50.0, 10.0 # 三层初始张力水蓄量(mm) S_free 5.0 # 初始自由水蓄量(mm) Q_sim []这段代码解决的问题是模型跑起来之前所有状态变量从哪来。参数初值不要求精确后续会率定。比较关键的是前三层初始蓄水量的设置——如果给得太小模拟初期会有一段填洼过程模拟流量偏低如果给得太大初期径流会虚高。3.2 三层蒸发计算模块核心的日水量平衡每天的演算遵循固定顺序先算蒸发再算产流最后更新各层蓄量。蒸发计算的先后不能颠倒因为蒸发消耗的是当天的土壤蓄水而产流用的是当天的降雨增量。# 三层蒸发计算逻辑 def calc_evap(P_day, EM_day, WU, WL, WD): # 1. 上层蒸发优先消耗上层水量 EU EM_day * K if WU EU: WU WU - EU else: EU WU WU 0 # 2. 下层蒸发剩余蒸散发能力按折减系数分配 EL (EM_day * K - EU) * 0.6 # 0.6是下层蒸发折减系数 if WL EL: WL WL - EL else: EL WL WL 0 # 3. 深层蒸发只占剩余蒸发的小部分 ED (EM_day * K - EU - EL) * 0.1 if WD ED: WD WD - ED else: ED WD WD 0 return EU, EL, ED, WU, WL, WD函数内部的三个if判断体现的是能供多少就蒸发多少的规则。比如上层含水量只有5mm但当天蒸发能力是8mm那就只能蒸发5mm余下的蒸发需求依次往下层、深层传递。这也解释了为什么连续干旱天模拟蒸发衰减快——上层很快耗尽下层折减系数0.6进一步限制了蒸发速率。3.3 蓄满产流和自由水蓄量降雨转化过程蒸发算完后再算产流。核心操作是判断加雨后的总张力水蓄量是否超过WM。如果超过超出部分进入自由水蓄量如果没超过降雨全部用于补充土壤缺水量不产流。这部分的顺序、条件、更新逻辑是模型的重中之重。# 产流计算模块 def calc_runoff(P_day, WU, WL, WD, S_free): # 净雨扣除不透水面积直接产流 Pe P_day * (1 - IM) P_imp P_day * IM # 不透水面积上的直接径流 # 1. 降雨先补充张力水 remaining P_day * (1 - IM) # 上层补充 if WU WUM: dWU min(WUM - WU, remaining) WU dWU remaining - dWU # 下层补充 if WL WLM and remaining 0: dWL min(WLM - WL, remaining) WL dWL remaining - dWL # 深层补充 if WD WDM and remaining 0: dWD min(WDM - WD, remaining) WD dWD remaining - dWD # 2. 张力水蓄满后剩余水量进入自由水蓄量 dS remaining # 剩余水量全部进入自由水蓄量 # 3. 自由水蓄量超过容量时产生地面径流 if S_free dS SM: RS (S_free dS - SM) S_free SM else: RS 0 S_free S_free dS # 4. 自由水的壤中流和地下径流分配 RI KI * S_free RG KG * S_free # 更新自由水蓄量扣除流出部分 S_free S_free - RI - RG # 总径流 地面径流 壤中流 地下径流 不透水面积直接径流 R_total RS RI RG P_imp return R_total, WU, WL, WD, S_free, RS, RI, RG产流划分的关键在于自由水蓄量的库容有限。降雨水超过SM的部分直接成为地面径流而壤中流和地下径流的量由KI、KG和当前S_free共同决定。如果把SM调大地面径流会减少壤中流增多洪水过程会变得坦化如果把SM调小地面径流占比大洪峰更尖更陡。这部分参数直接控制模拟的洪峰形态。3.4 汇流计算三种水源分别演算产流计算结束后地面径流RS直接进入出口断面壤中流RI经过一个线性水库衰减地下径流RG经过另一个衰减更慢的线性水库。线性水库的退水公式是今天出流量等于昨天的剩余蓄量乘以消退系数。# 汇流模块线性水库演算 def route_flow(RS_day, RI_day, RG_day, S_RI_prev, S_RG_prev, CS, CG): # 地面径流直接汇入但做单位线或线性水库处理 # 这里使用一个简单线性水库模拟地面汇流 Q_surface CS * S_RI_prev (1 - CS) * RS_day S_surface S_RI_prev * (1 - CS) (1 - CS) * RS_day # 简化示意 # 壤中流线性水库 Q_inter CS * S_RI_prev (1 - CS) * RI_day S_RI_prev_new S_RI_prev * (1 - CS) RI_day * (1 - CS) # 地下径流线性水库 Q_ground CG * S_RG_prev (1 - CG) * RG_day S_RG_prev_new S_RG_prev * (1 - CG) RG_day * (1 - CG) # 总出流量转换为m3/s Q_total (Q_surface Q_inter Q_ground) * 流域面积_km2 * 1000 / 86400 return Q_total, S_RI_prev_new, S_RG_prev_new汇流部分在实际工程中经常被简化。如果只有一个日尺度的流量序列用线性水库就够了不需要复杂的单位线。参数CS反映地面汇流速度取值0.3~0.8越大表示汇流越快CG一般取0.9以上表示地下水消退非常缓慢。这里单位换算要注意每天的径流深(mm)乘以流域面积再除以86400秒才能转换为m3/s。4. 参数率定和敏感性分析让模拟结果靠近实测值4.1 人工试错法和纳什效率系数模型参数有七八个直接上来就全自动寻优容易陷入局部最优而且你不知道哪个参数起了主要作用。常见做法是先手动调一遍让模拟过程线大致贴合实测再用自动优化微调。评价模拟好坏的指标最常用纳什效率系数NSE它衡量模拟和实测的吻合程度计算公式是1减去误差方差除以实测方差。def calc_nse(Q_obs, Q_sim): obs_mean np.mean(Q_obs) numerator np.sum((Q_obs - Q_sim) ** 2) denominator np.sum((Q_obs - obs_mean) ** 2) return 1 - numerator / denominatorNSE等于1表示完美拟合0表示模拟效果等同用实测均值预测负数表示模型比均值预测还差。实际工程中NSE达到0.7以上就算可用0.8以上效果良好。调参顺序也有讲究先调WM和K让年水量平衡闭合再调SM控制洪峰形态最后调KI、KG控制退水过程。4.2 参数敏感性哪个参数最影响结果我做过一组实验把各个参数单独拉偏20%观察NSE变化幅度。结果最敏感的是蒸发折算系数K它对年径流总量影响极大——K从1.0调到1.2年径流可能会减少20%以上因为多算的蒸散发直接减少了产流。其次是WM它控制着流域的总体蓄水能力WM太小模拟洪峰偏高WM太大则洪峰偏低。SM和KI、KG的敏感性相对中等它们主要影响过程的形态而不是总量。参数率定建议先跑一次基准模拟输出年径流量和实测年径流量对比。如果模拟偏大大概率是K偏小或者WM偏小优先调这两个。如果年总量对上了但洪峰偏高再调SM。这种分步调整比一次性盲目搜索效率高得多也更容易攒出调参经验。4.3 率定脚本批量搜索加人工判断# 简单的网格搜索率定 best_nse -999 best_params {} for test_WM in [100, 120, 150, 180]: for test_SM in [10, 15, 20, 25]: for test_K in [0.8, 1.0, 1.2]: # 运行模型得到Q_sim nse calc_nse(Q_obs, Q_sim) if nse best_nse: best_nse nse best_params {WM: test_WM, SM: test_SM, K: test_K}网格搜索看起来粗暴但胜在直观可靠。三层嵌套的循环把所有组合跑一遍每组参数都调用一次完整的模拟计算最后选出NSE最高的一组。注意每组参数跑模拟之前要让模型有一段预热期比如前30天的模拟结果不参与NSE计算用来消除初始条件的影响。5. 模型调试避坑运行报错和结果失真排查5.1 数值溢出Python报numpy.inf或nan现象运行一段时间后模拟流量变成inf或nan过程线断掉。原因通常是某个状态变量被算成了负值或无穷大。比如蒸发计算时上层含水量已经是0但代码里仍然尝试减去蒸发量导致负值产生。累计下去自由水蓄量或张力水蓄量变成无意义的大数最终溢出。解决在每个计算函数入口加一个保护判断任何蓄量变量如果出现负值就强制归零另一个常见原因是蒸发输入数据有缺测或零值导致EM_day * K为0但后续计算中分母出现了0。检查原始Excel数据里的空值用pandas的dropna或者在读取后直接用fillna(0)填充。WU max(WU, 0) WL max(WL, 0) WD max(WD, 0) S_free max(S_free, 0)5.2 模拟径流一直为0检查蓄满条件是否永远不满足现象模型跑完模拟流量全是0一条直线。原因大概率是WM过大或降雨太小所有降雨都被土壤吸收从未蓄满。另一种情况是蒸散发输入值异常大导致土壤始终处于亏水状态。解决第一步在产流计算模块里加一行print把每天的WUWLWD和WM对比输出观察蓄水量是否在逼近WM。如果连续几十天蓄量都在WM以下说明降雨和蒸发不匹配。第二步把SM调小让它更容易产流。第三步检查蒸发数据是否使用了总的潜在蒸散发而不是实际蒸散发能力。如果确实用了潜在的考虑把K调小到0.5~0.7来折减。5.3 洪水过程模拟不错退水过程偏高或偏低现象洪峰对上了但退水尾巴太慢或者太快。原因在于KI和KG的比例不对。KG太大地下水补充太快退水过程降不下来KG太小退水过程又太平缓。另一个因素时线性水库的消退系数CS和CG设置不合理CS取值太小造成地面汇流延迟过大导致洪峰后还有一个不小的拖尾流量。解决单独拿出退水段的数据画出半对数图。如果实测退水段在半对数坐标下是一条直线说明退水主要受单一线性水库控制这时候就只需要调CG。如果实测退水是折线前段陡后段缓说明地面-壤中流-地下水三分水源的消退系数都需要调整。具体到代码先固定CS0.5、CG0.95再单独微调KG从0.2到0.4看NSE变化。5.4 年径流总量偏差大引入水量平衡误差检查现象模拟年径流总量比实测偏大30%以上。原因可能是不透水面积比例IM设置过高或者蒸散发能力数据偏小。IM直接决定有多少降雨不参与蓄满产流而直接成为径流IM每增加0.1年径流就会显著增加。解决在代码末尾加一个年度水量平衡的统计输出打印降雨总量、蒸散发总量、径流总量三项。如果降雨减去蒸发减去径流后剩余的水量不为零说明存在计算bug或蓄量变化不闭合。逐项检查后优先把IM调回0.05以下同时把K往1.1~1.2方向调整让蒸散发消耗掉更多水分。5.5 初值影响前几天的模拟结果预热期不够现象模拟前10天结果和实测差距很大后面逐渐正常。原因是初始土壤含水量设置不合理。如果初始WU、WL、WD设得太小模型需要一段时间填蓄这段时间模拟径流几乎为0。解决常见的处理方式有两种。第一种是把模拟起点提前一年最后使用后半段的结果。第二种是在代码中增加一个预热循环使用实测降雨的前90天数据连续计算让状态变量收敛到合理范围然后再开始记录模拟结果。第二种更简单不用额外拼接数据只需在计算循环开始前多跑几十次同样的模拟逻辑。6. 参数边界和率定顺序的实操心得用好模型的关键技巧调参顺序比参数值本身更重要。这条建议反复说因为很多同学上来就把全部参数扔给自动寻优算法结果几千次迭代后NSE还没到0.6回头问为什么。原因是参数之间存在明显的相关性比如WM和K都会影响年径流总量如果不固定其中一个直接一起搜算法会陷入斜向的山脊里怎么跳都找不到最优。我现在的习惯是先固定各层张力水容量比例为2:9:4K取1.0SM取15KG取0.3也就是用一套普通参数把模型跑通看过程线的基本形态再进入调参阶段。第一步只调K紧盯年径流总量让水量平衡误差在5%以内。第二步调WM确保湿润期产流不偏早不偏晚。第三步调SM修正洪峰的坦化和锐化。第四步调KG和KI修正退水过程。第五步才允许微调各层蒸发系数。中间还有几个细节容易被忽略自由水蓄量的初始值S_free如果设置过大前几天的壤中流会异常偏大这个问题在北方干旱流域尤其严重另一个常见问题是蒸发数据的量纲日蒸散发能力一般小于8mm如果某人给的数据单位是月蒸发量而你直接当成日值用了模型会一直蓄不满水还有流域面积的单位换算调用汇流模块前确认面积单位是km2而不是m2或者亩这个错了输出流量会差好几个数量级。我的做法是在代码开头统一加一个assert断言检查关键参数的范围比如0KG1、WM50、SM0超出范围直接报错避免带着错误参数跑几百行代码最后才发现问题。跑完模型之后还有一件值得做的事把模拟和实测的双累积曲线画出来对比。双累积曲线的斜率如果在中后段明显改变大概率是参数K需要分段取值或者流域内有水利工程调蓄导致下垫面条件变化。我自己做过一个案例某河流域前半段模拟效果很好后半段模拟流量整体偏低查了资料发现后段有多个小型水库开始蓄水蒸发和调度都变了再调K也救不回来只能分段率定。水文模型不是万能钥匙参数只能反映下垫面的平均状态子流域的调蓄改变属于结构性误差这时候往参数上找补只会让模型过拟合。经验上一篇完整的新安江模型Python程序应该包含参数配置文件建议用yaml或json分开写、核心计算模块蒸发、产流、汇流三个函数、率定模块网格搜索或遗传算法、结果输出模块生成过程线图和NSE、水量平衡误差报告。这些代码加起来大概500行不算长但结构清晰后续无论是换流域还是改时间尺度都能快速适配。希望这篇笔记能把你的新安江模型程序从能跑带到跑得明白、调得有方向祝你好运。本文还有配套的精品资源点击获取 新安江模型Python实现：三层蒸发蓄满产流从降雨到径流模拟
简介这份资源是一份面向水文水资源专业学生与从业者的三层蒸发蓄满产流模型新安江模型Python计算程序用于完成基于降水与蒸发能力数据的产流模拟计算。程序核心实现新安江模型中的三层蒸发与蓄满产流机制适合学习流域水文模型、开展课程设计或进行科研验证的初级至中级Python使用者。资源包为1个py文件体积约1KB内容紧凑主程序需配合作者提供的PEdata.xlsx降水蒸发数据运行可直接复制使用或在此基础上修改参数。截至目前已有2642人浏览学习说明该程序在同类资源中具有较高实用参考价值。通过运行这份代码读者可以直观理解新安江模型的计算流程快速上手实现产流模拟为后续开展流域水文分析提供可复用的脚本基础。1. 三层蒸发蓄满产流模型新安江模型python计算程序.py从降雨到径流的确定性推演做水文预报的人迟早会碰到这个场景手里有流域的逐日降雨和蒸发数据想算出对应的径流过程但商业化软件黑匣子太多参数一改就不知道内部发生了什么。三层蒸发蓄满产流模型常说的新安江模型就是干这个事的——它把流域降水转化成蒸散发、土壤蓄水、地表径流和壤中流是一个完整闭合的水量平衡计算框架。标题里这个python计算程序.py本质是一份可以直接运行的脚本解决的是没有现成水文模型库时如何用几百行代码实现经典蓄满产流计算的现实问题。适合水文专业学生做课程设计、科研人员做对比实验以及工程师做中小流域的洪水预报方案初筛。只要你有逐日降雨、蒸发和流量数据这套代码就能从零开始推演出模拟径流过程还能人工调参逼近实测值。核心并不神秘蓄满产流的核心就是先算土壤缺水量再决定产不产流。2. 蓄满产流与三层蒸发模型运行的底层逻辑2.1 张力水和自由水蓄满产流的核心区分蓄满产流模型能成立是因为它区分了两种性质完全不同的土壤水。张力水是被土粒吸附、能被植物根系吸收的水分它不参与重力排水自由水则是填充在土壤孔隙中、受重力支配可以横向或纵向流动的水分。降雨先补充张力水张力水蓄满后才开始产生自由水自由水超过一定容量就形成径流。这个先蓄满再产流的机制是模型最关键的判断条件。实际计算中每个计算单元都有一个张力水蓄水容量WM它等于上层张力水容量WUM、下层张力水容量WLM和深层张力水容量WDM之和。降雨先补上层上层满了溢到下层下层满了溢到深层。只有当整个土柱的张力水都蓄满降雨才转为自由水蓄量进而划分产流。这个顺序保证了水量平衡的闭合性也让参数物理意义明确——WUM和WLM可以根据土壤类型和经验值初定再通过率定微调。2.2 三层蒸发为什么不能只算一层很多初学模型的人会问一个流域的总蒸发量乘以一个系数不就完了为什么要分三层答案在于干旱期和湿润期的蒸发规律截然不同。湿润期土壤含水充足蒸发接近蒸散发能力干旱期表层先干蒸发衰减很快但深层仍有水分缓慢供给。三层蒸发结构正是为了模拟这种表层快、深层慢的衰减过程。模型里每层有自己的蒸发能力分配系数。上层蒸发EU取蒸散发能力EM乘以系数K如果上层有足够水量就直接消耗不够时下层蒸发EL取剩余蒸散发能力乘以一个折减系数但最多不能超过下层张力水蓄量的某个比例再不够才动用深层ED深层蒸发速率接近常数。这一套规则写进代码就三个条件分支却决定了连续干旱期的模拟精度。2.3 从降雨到径流水源划分和汇流张力水全部蓄满后超出的降雨进入自由水蓄量SM。自由水蓄量有一个最大的容量超过了就产地面径流RS同时自由水还按一个固定的消退系数KG向下补给深层地下水形成地下径流RG。壤中流RI则和自由水蓄量成正比计算为KI乘以自由水蓄量其中KI和KG是消退系数。降雨产流过程结束后各水源按不同路径进入河道地面径流几乎即时汇入壤中流衰减稍慢地下径流衰减最慢。模型通常用一个线性水库或单位线把三种水源分别演算到流域出口。这个先分水源、再分路径的做法让模型既能模拟洪水过程的陡涨陡落又能模拟退水过程的缓降比单一水库模型灵活得多。3. 编写新安江模型Python计算程序蓄满产流演算模块的完整实现3.1 数据输入和初始化直接读取Excel更方便我一般用pandas读取Excel里的三列数据日期、降雨量、实测流量。蒸散发能力可以单独一列也可以由程序根据月份自动调整。初始化时先定义土壤含水量初值、蓄水容量参数、蒸发系数以及各水源的消退系数。代码定义得越清晰后面调参越容易。import pandas as pd import numpy as np # 读取数据 df pd.read_excel(流域日数据.xlsx, sheet_nameSheet1) dates df[日期].values P df[降雨].values.astype(float) # 日降雨单位mm EM_input df[蒸发].values.astype(float) # 日蒸散发能力单位mm Q_obs df[流量].values.astype(float) # 实测流量单位m3/s用于对比 # 模型参数初始值 WUM, WLM, WDM 20, 90, 40 # 上、下、深层张力水容量(mm) WM WUM WLM WDM # 总张力水容量(mm) K 1.0 # 蒸散发能力折算系数 SM 15.0 # 自由水蓄水容量(mm) KG, KI 0.3, 0.35 # 地下水和壤中流出流系数 IM 0.1 # 不透水面积比例 FC 10.0 # 稳定下渗率(mm/day) # 注原文此处为帮助记忆的占位符 # 初值 WU, WL, WD 10.0, 50.0, 10.0 # 三层初始张力水蓄量(mm) S_free 5.0 # 初始自由水蓄量(mm) Q_sim []这段代码解决的问题是模型跑起来之前所有状态变量从哪来。参数初值不要求精确后续会率定。比较关键的是前三层初始蓄水量的设置——如果给得太小模拟初期会有一段填洼过程模拟流量偏低如果给得太大初期径流会虚高。3.2 三层蒸发计算模块核心的日水量平衡每天的演算遵循固定顺序先算蒸发再算产流最后更新各层蓄量。蒸发计算的先后不能颠倒因为蒸发消耗的是当天的土壤蓄水而产流用的是当天的降雨增量。# 三层蒸发计算逻辑 def calc_evap(P_day, EM_day, WU, WL, WD): # 1. 上层蒸发优先消耗上层水量 EU EM_day * K if WU EU: WU WU - EU else: EU WU WU 0 # 2. 下层蒸发剩余蒸散发能力按折减系数分配 EL (EM_day * K - EU) * 0.6 # 0.6是下层蒸发折减系数 if WL EL: WL WL - EL else: EL WL WL 0 # 3. 深层蒸发只占剩余蒸发的小部分 ED (EM_day * K - EU - EL) * 0.1 if WD ED: WD WD - ED else: ED WD WD 0 return EU, EL, ED, WU, WL, WD函数内部的三个if判断体现的是能供多少就蒸发多少的规则。比如上层含水量只有5mm但当天蒸发能力是8mm那就只能蒸发5mm余下的蒸发需求依次往下层、深层传递。这也解释了为什么连续干旱天模拟蒸发衰减快——上层很快耗尽下层折减系数0.6进一步限制了蒸发速率。3.3 蓄满产流和自由水蓄量降雨转化过程蒸发算完后再算产流。核心操作是判断加雨后的总张力水蓄量是否超过WM。如果超过超出部分进入自由水蓄量如果没超过降雨全部用于补充土壤缺水量不产流。这部分的顺序、条件、更新逻辑是模型的重中之重。# 产流计算模块 def calc_runoff(P_day, WU, WL, WD, S_free): # 净雨扣除不透水面积直接产流 Pe P_day * (1 - IM) P_imp P_day * IM # 不透水面积上的直接径流 # 1. 降雨先补充张力水 remaining P_day * (1 - IM) # 上层补充 if WU WUM: dWU min(WUM - WU, remaining) WU dWU remaining - dWU # 下层补充 if WL WLM and remaining 0: dWL min(WLM - WL, remaining) WL dWL remaining - dWL # 深层补充 if WD WDM and remaining 0: dWD min(WDM - WD, remaining) WD dWD remaining - dWD # 2. 张力水蓄满后剩余水量进入自由水蓄量 dS remaining # 剩余水量全部进入自由水蓄量 # 3. 自由水蓄量超过容量时产生地面径流 if S_free dS SM: RS (S_free dS - SM) S_free SM else: RS 0 S_free S_free dS # 4. 自由水的壤中流和地下径流分配 RI KI * S_free RG KG * S_free # 更新自由水蓄量扣除流出部分 S_free S_free - RI - RG # 总径流 地面径流 壤中流 地下径流 不透水面积直接径流 R_total RS RI RG P_imp return R_total, WU, WL, WD, S_free, RS, RI, RG产流划分的关键在于自由水蓄量的库容有限。降雨水超过SM的部分直接成为地面径流而壤中流和地下径流的量由KI、KG和当前S_free共同决定。如果把SM调大地面径流会减少壤中流增多洪水过程会变得坦化如果把SM调小地面径流占比大洪峰更尖更陡。这部分参数直接控制模拟的洪峰形态。3.4 汇流计算三种水源分别演算产流计算结束后地面径流RS直接进入出口断面壤中流RI经过一个线性水库衰减地下径流RG经过另一个衰减更慢的线性水库。线性水库的退水公式是今天出流量等于昨天的剩余蓄量乘以消退系数。# 汇流模块线性水库演算 def route_flow(RS_day, RI_day, RG_day, S_RI_prev, S_RG_prev, CS, CG): # 地面径流直接汇入但做单位线或线性水库处理 # 这里使用一个简单线性水库模拟地面汇流 Q_surface CS * S_RI_prev (1 - CS) * RS_day S_surface S_RI_prev * (1 - CS) (1 - CS) * RS_day # 简化示意 # 壤中流线性水库 Q_inter CS * S_RI_prev (1 - CS) * RI_day S_RI_prev_new S_RI_prev * (1 - CS) RI_day * (1 - CS) # 地下径流线性水库 Q_ground CG * S_RG_prev (1 - CG) * RG_day S_RG_prev_new S_RG_prev * (1 - CG) RG_day * (1 - CG) # 总出流量转换为m3/s Q_total (Q_surface Q_inter Q_ground) * 流域面积_km2 * 1000 / 86400 return Q_total, S_RI_prev_new, S_RG_prev_new汇流部分在实际工程中经常被简化。如果只有一个日尺度的流量序列用线性水库就够了不需要复杂的单位线。参数CS反映地面汇流速度取值0.3~0.8越大表示汇流越快CG一般取0.9以上表示地下水消退非常缓慢。这里单位换算要注意每天的径流深(mm)乘以流域面积再除以86400秒才能转换为m3/s。4. 参数率定和敏感性分析让模拟结果靠近实测值4.1 人工试错法和纳什效率系数模型参数有七八个直接上来就全自动寻优容易陷入局部最优而且你不知道哪个参数起了主要作用。常见做法是先手动调一遍让模拟过程线大致贴合实测再用自动优化微调。评价模拟好坏的指标最常用纳什效率系数NSE它衡量模拟和实测的吻合程度计算公式是1减去误差方差除以实测方差。def calc_nse(Q_obs, Q_sim): obs_mean np.mean(Q_obs) numerator np.sum((Q_obs - Q_sim) ** 2) denominator np.sum((Q_obs - obs_mean) ** 2) return 1 - numerator / denominatorNSE等于1表示完美拟合0表示模拟效果等同用实测均值预测负数表示模型比均值预测还差。实际工程中NSE达到0.7以上就算可用0.8以上效果良好。调参顺序也有讲究先调WM和K让年水量平衡闭合再调SM控制洪峰形态最后调KI、KG控制退水过程。4.2 参数敏感性哪个参数最影响结果我做过一组实验把各个参数单独拉偏20%观察NSE变化幅度。结果最敏感的是蒸发折算系数K它对年径流总量影响极大——K从1.0调到1.2年径流可能会减少20%以上因为多算的蒸散发直接减少了产流。其次是WM它控制着流域的总体蓄水能力WM太小模拟洪峰偏高WM太大则洪峰偏低。SM和KI、KG的敏感性相对中等它们主要影响过程的形态而不是总量。参数率定建议先跑一次基准模拟输出年径流量和实测年径流量对比。如果模拟偏大大概率是K偏小或者WM偏小优先调这两个。如果年总量对上了但洪峰偏高再调SM。这种分步调整比一次性盲目搜索效率高得多也更容易攒出调参经验。4.3 率定脚本批量搜索加人工判断# 简单的网格搜索率定 best_nse -999 best_params {} for test_WM in [100, 120, 150, 180]: for test_SM in [10, 15, 20, 25]: for test_K in [0.8, 1.0, 1.2]: # 运行模型得到Q_sim nse calc_nse(Q_obs, Q_sim) if nse best_nse: best_nse nse best_params {WM: test_WM, SM: test_SM, K: test_K}网格搜索看起来粗暴但胜在直观可靠。三层嵌套的循环把所有组合跑一遍每组参数都调用一次完整的模拟计算最后选出NSE最高的一组。注意每组参数跑模拟之前要让模型有一段预热期比如前30天的模拟结果不参与NSE计算用来消除初始条件的影响。5. 模型调试避坑运行报错和结果失真排查5.1 数值溢出Python报numpy.inf或nan现象运行一段时间后模拟流量变成inf或nan过程线断掉。原因通常是某个状态变量被算成了负值或无穷大。比如蒸发计算时上层含水量已经是0但代码里仍然尝试减去蒸发量导致负值产生。累计下去自由水蓄量或张力水蓄量变成无意义的大数最终溢出。解决在每个计算函数入口加一个保护判断任何蓄量变量如果出现负值就强制归零另一个常见原因是蒸发输入数据有缺测或零值导致EM_day * K为0但后续计算中分母出现了0。检查原始Excel数据里的空值用pandas的dropna或者在读取后直接用fillna(0)填充。WU max(WU, 0) WL max(WL, 0) WD max(WD, 0) S_free max(S_free, 0)5.2 模拟径流一直为0检查蓄满条件是否永远不满足现象模型跑完模拟流量全是0一条直线。原因大概率是WM过大或降雨太小所有降雨都被土壤吸收从未蓄满。另一种情况是蒸散发输入值异常大导致土壤始终处于亏水状态。解决第一步在产流计算模块里加一行print把每天的WUWLWD和WM对比输出观察蓄水量是否在逼近WM。如果连续几十天蓄量都在WM以下说明降雨和蒸发不匹配。第二步把SM调小让它更容易产流。第三步检查蒸发数据是否使用了总的潜在蒸散发而不是实际蒸散发能力。如果确实用了潜在的考虑把K调小到0.5~0.7来折减。5.3 洪水过程模拟不错退水过程偏高或偏低现象洪峰对上了但退水尾巴太慢或者太快。原因在于KI和KG的比例不对。KG太大地下水补充太快退水过程降不下来KG太小退水过程又太平缓。另一个因素时线性水库的消退系数CS和CG设置不合理CS取值太小造成地面汇流延迟过大导致洪峰后还有一个不小的拖尾流量。解决单独拿出退水段的数据画出半对数图。如果实测退水段在半对数坐标下是一条直线说明退水主要受单一线性水库控制这时候就只需要调CG。如果实测退水是折线前段陡后段缓说明地面-壤中流-地下水三分水源的消退系数都需要调整。具体到代码先固定CS0.5、CG0.95再单独微调KG从0.2到0.4看NSE变化。5.4 年径流总量偏差大引入水量平衡误差检查现象模拟年径流总量比实测偏大30%以上。原因可能是不透水面积比例IM设置过高或者蒸散发能力数据偏小。IM直接决定有多少降雨不参与蓄满产流而直接成为径流IM每增加0.1年径流就会显著增加。解决在代码末尾加一个年度水量平衡的统计输出打印降雨总量、蒸散发总量、径流总量三项。如果降雨减去蒸发减去径流后剩余的水量不为零说明存在计算bug或蓄量变化不闭合。逐项检查后优先把IM调回0.05以下同时把K往1.1~1.2方向调整让蒸散发消耗掉更多水分。5.5 初值影响前几天的模拟结果预热期不够现象模拟前10天结果和实测差距很大后面逐渐正常。原因是初始土壤含水量设置不合理。如果初始WU、WL、WD设得太小模型需要一段时间填蓄这段时间模拟径流几乎为0。解决常见的处理方式有两种。第一种是把模拟起点提前一年最后使用后半段的结果。第二种是在代码中增加一个预热循环使用实测降雨的前90天数据连续计算让状态变量收敛到合理范围然后再开始记录模拟结果。第二种更简单不用额外拼接数据只需在计算循环开始前多跑几十次同样的模拟逻辑。6. 参数边界和率定顺序的实操心得用好模型的关键技巧调参顺序比参数值本身更重要。这条建议反复说因为很多同学上来就把全部参数扔给自动寻优算法结果几千次迭代后NSE还没到0.6回头问为什么。原因是参数之间存在明显的相关性比如WM和K都会影响年径流总量如果不固定其中一个直接一起搜算法会陷入斜向的山脊里怎么跳都找不到最优。我现在的习惯是先固定各层张力水容量比例为2:9:4K取1.0SM取15KG取0.3也就是用一套普通参数把模型跑通看过程线的基本形态再进入调参阶段。第一步只调K紧盯年径流总量让水量平衡误差在5%以内。第二步调WM确保湿润期产流不偏早不偏晚。第三步调SM修正洪峰的坦化和锐化。第四步调KG和KI修正退水过程。第五步才允许微调各层蒸发系数。中间还有几个细节容易被忽略自由水蓄量的初始值S_free如果设置过大前几天的壤中流会异常偏大这个问题在北方干旱流域尤其严重另一个常见问题是蒸发数据的量纲日蒸散发能力一般小于8mm如果某人给的数据单位是月蒸发量而你直接当成日值用了模型会一直蓄不满水还有流域面积的单位换算调用汇流模块前确认面积单位是km2而不是m2或者亩这个错了输出流量会差好几个数量级。我的做法是在代码开头统一加一个assert断言检查关键参数的范围比如0KG1、WM50、SM0超出范围直接报错避免带着错误参数跑几百行代码最后才发现问题。跑完模型之后还有一件值得做的事把模拟和实测的双累积曲线画出来对比。双累积曲线的斜率如果在中后段明显改变大概率是参数K需要分段取值或者流域内有水利工程调蓄导致下垫面条件变化。我自己做过一个案例某河流域前半段模拟效果很好后半段模拟流量整体偏低查了资料发现后段有多个小型水库开始蓄水蒸发和调度都变了再调K也救不回来只能分段率定。水文模型不是万能钥匙参数只能反映下垫面的平均状态子流域的调蓄改变属于结构性误差这时候往参数上找补只会让模型过拟合。经验上一篇完整的新安江模型Python程序应该包含参数配置文件建议用yaml或json分开写、核心计算模块蒸发、产流、汇流三个函数、率定模块网格搜索或遗传算法、结果输出模块生成过程线图和NSE、水量平衡误差报告。这些代码加起来大概500行不算长但结构清晰后续无论是换流域还是改时间尺度都能快速适配。希望这篇笔记能把你的新安江模型程序从能跑带到跑得明白、调得有方向祝你好运。本文还有配套的精品资源点击获取