简介本资源是面向水文与水资源工程专业本科生及水文模型初学者的课程设计实践代码包聚焦新安江模型日尺度模拟与三水源划分核心任务完整覆盖地面产流、汇流计算及预报精度评估等关键环节可直接用于河海大学《水文预报》课程设计作业或模型原理验证。压缩包仅含2个文件1个Python主程序1个说明文本体积精简至6KBpy文件实现自动读取Excel输入数据、按年分段计算误差指标如确定性系数NSE与三水源比例txt文件提供参数说明与运行指引结构简洁、即装即用。目前已有295人学习下载适合需要快速上手新安江模型Python实现、理解日模型参数率定逻辑与精度评价方法的学习者。1. 新安江模型日模型三水划分为什么用 Python 复现河海大学水文预报课程设计成了水文专业学生绕不开的“硬核入门关”2022年河海大学水文预报课程设计里那个被反复调试、手写公式、Excel 校核、最后在 Python 里跑通的新安江模型日模型三水划分实现——它不是一段炫技代码而是一道真实存在的工程门槛你得亲手把产流计算、水源划分、汇流演算三个模块串成闭环让一组实测雨量、蒸发、流量数据真正输出符合水文物理机制的日径流过程线。很多同学卡在“三水源怎么分”这一步表层流、地下流、壤中流不是按比例切蛋糕而是依赖蓄水容量曲线、自由水蓄水库结构、以及逐日动态水量平衡。更关键的是这个模型对初值敏感、对蒸散发估算鲁棒性差、对暴雨峰形响应存在相位滞后——这些都不是教科书习题能暴露的问题只有在 Python 里一行行调试、画出中间变量时才突然意识到“原来新安江模型的‘经验性’就藏在那几个查算表和迭代初值里。”适合刚学完《水文水利计算》《水资源系统分析》、正要接触实际建模流程的本科生也适合需要快速搭建教学演示原型、验证参数敏感性的青年教师。它不追求工业级精度但必须可追溯、可修改、可解释——这才是课程设计该有的样子。2. 从概念到代码新安江模型日模型三水划分的核心逻辑与 Python 实现路径新安江模型日模型不是黑匣子它的“三水”本质是产流机制的空间解耦表层流坡面汇流、壤中流包气带侧向渗流、地下流深层基流。但实现上它通过一套嵌套的“蓄水—产流—分水源—汇流”链条完成。核心不在算法多复杂而在状态变量如何传递、查算表如何嵌入、迭代如何收敛。我一般会先画一张“日推进”状态图前一天的土壤含水量 W、自由水蓄水库蓄量 SM、地下水库蓄量 EX共同决定当天的产流量、三水源比例、以及下一天的初值。整个过程没有微分方程求解全是代数迭代查表插值——这正是 Python 能胜任且优于 Fortran 的地方NumPy 向量化处理日序列Pandas 管理实测输入输出Matplotlib 直观诊断中间过程。2.1 输入数据准备与格式规范为什么必须严格按“日期-雨量-蒸发-流量”四列组织课程设计给定的数据通常是 Excel 表格但直接读进 Python 容易翻车。常见错误是日期列识别失败、空值未填充、单位不统一如蒸发用 mm/d但模型内部默认 mm。我坚持用以下结构初始化 DataFrameimport pandas as pd import numpy as np # 读取原始数据假设文件名为 input_data.xlsx df pd.read_excel(input_data.xlsx, parse_dates[date], index_coldate) # 强制重采样为日频填充缺失用前向填充后向填充兜底 df df.asfreq(D).fillna(methodffill).fillna(methodbfill) # 检查并标准化单位所有输入统一为 mm/d df[rainfall] df[rainfall].apply(lambda x: max(0, x)) # 雨量不能为负 df[evaporation] df[evaporation].apply(lambda x: max(0, x)) # 流量需转为等效深度mm/d需提供流域面积km² area_km2 1250.0 # 示例值需根据题目给定替换 df[flow_mm] (df[flow_m3s] * 86400 / (area_km2 * 1e6)) * 1000 # m³/s → mm/d提示asfreq(D)是关键。新安江模型严格按日历日推进跳过某天或日期重复会导致状态链断裂。flow_m3s列名需按你实际数据调整若原始数据无流量该列可全设为 NaN仅用于后续率定验证。2.2 核心参数初始化与物理意义说明哪些参数必须手调哪些可默认新安江模型日模型有 15 个左右参数但课程设计聚焦三水划分只需关注以下 9 个核心参数。它们不是随便填的数字每个都对应明确的水文物理含义参数符号物理含义典型取值范围课程设计建议初值调参敏感度WM土壤蓄水容量mm30–12075.0⭐⭐⭐⭐B蓄水容量曲线指数0.1–0.30.2⭐⭐⭐IM初始不透水面积比0.01–0.150.05⭐⭐SM自由水蓄水库容量mm10–5025.0⭐⭐⭐⭐⭐EX地下水蓄水库容量mm5–3015.0⭐⭐⭐KI壤中流出流系数0.05–0.30.15⭐⭐⭐KG地下流出流系数0.001–0.020.005⭐⭐⭐⭐C河网汇流系数日模型中常简化0.5–0.950.8⭐⭐L河网滞时日0.5–3.01.0⭐⭐⭐注意SM和EX是自由水蓄水库和地下水蓄水库的“最大容量”不是初始蓄量。初始蓄量需单独设置如SM00.5*SM,EX00.3*EX这是率定时的重要自由度。C和L在纯日模型中常合并为一个“汇流单位线”参数课程设计常用简化线性汇流Q_out[t] C * Q_in[t] (1-C) * Q_out[t-1]。2.3 三水源划分的逐日计算逻辑从产流到分水的完整推导链三水划分发生在产流之后、汇流之前是新安江模型最体现经验性的环节。其逻辑链如下以第 t 日为例计算当日净雨P_eff[t] max(0, P[t] - E[t])忽略前期影响课程设计常简化计算土壤缺水量Def[t] WM * (1 - (W[t-1]/WM)^B)W[t-1] 为昨日土壤含水量产流计算RO[t] P_eff[t] if P_eff[t] Def[t] else 0超渗产流自由水蓄水库入流IN_SM[t] RO[t] * (1 - IM)扣除不透水区三水源分配表层流地表径流RO_S[t] IN_SM[t] * min(1, SM[t-1]/SM)当 SM 库满时全为表层流壤中流RO_SS[t] IN_SM[t] * KI * (1 - min(1, SM[t-1]/SM))地下流RO_G[t] IN_SM[t] * KG * (1 - min(1, SM[t-1]/SM))更新水库蓄量SM[t] SM[t-1] IN_SM[t] - RO_S[t] - RO_SS[t] - RO_G[t]EX[t] EX[t-1] RO_G[t] - KG * EX[t-1]地下水库出流按线性水库这个链条必须严格按顺序执行且所有变量需初始化如W[0] 0.5*WM,SM[0] 0.5*SM,EX[0] 0.3*EX。下面给出核心计算函数def xaj_daily_step(t, df, params, state): 新安江模型日模型单步计算 :param t: 当前日索引整数 :param df: 输入数据DataFrame含rainfall,evaporation,flow_mm :param params: 参数字典含WM,B,IM,SM,EX,KI,KG,C,L :param state: 状态字典含W,SM,EX均为标量t-1时刻值 :return: tuple (RO_S, RO_SS, RO_G, new_state) P df[rainfall].iloc[t] E df[evaporation].iloc[t] # 1. 净雨 P_eff max(0, P - E) # 2. 土壤缺水量 Def params[WM] * (1 - (state[W]/params[WM])**params[B]) # 3. 产流 RO P_eff if P_eff Def else 0 # 4. 自由水蓄水库入流 IN_SM RO * (1 - params[IM]) # 5. 三水源划分 # 表层流当SM库满时全部为表层流 if state[SM] params[SM]: RO_S IN_SM RO_SS 0.0 RO_G 0.0 else: # 壤中流和地下流按比例分配剩余入流 frac (params[SM] - state[SM]) / params[SM] # 可用空间比例 RO_S IN_SM * (1 - frac) RO_SS IN_SM * params[KI] * frac RO_G IN_SM * params[KG] * frac # 6. 更新状态 new_SM state[SM] IN_SM - RO_S - RO_SS - RO_G new_EX state[EX] RO_G - params[KG] * state[EX] # 土壤含水量更新经典新安江W[t] W[t-1] P_eff - E - RO new_W state[W] P_eff - E - RO # 状态约束不能为负不能超容 new_W max(0, min(new_W, params[WM])) new_SM max(0, min(new_SM, params[SM])) new_EX max(0, min(new_EX, params[EX])) new_state {W: new_W, SM: new_SM, EX: new_EX} return RO_S, RO_SS, RO_G, new_state # 初始化状态 params { WM: 75.0, B: 0.2, IM: 0.05, SM: 25.0, EX: 15.0, KI: 0.15, KG: 0.005, C: 0.8, L: 1.0 } state {W: 0.5 * params[WM], SM: 0.5 * params[SM], EX: 0.3 * params[EX]} # 主循环 n_days len(df) RO_S_arr np.zeros(n_days) RO_SS_arr np.zeros(n_days) RO_G_arr np.zeros(n_days) for t in range(n_days): RO_S, RO_SS, RO_G, state xaj_daily_step(t, df, params, state) RO_S_arr[t] RO_S RO_SS_arr[t] RO_SS RO_G_arr[t] RO_G逻辑说明此函数严格遵循新安江模型日模型原始文献赵人俊1984的三水源划分逻辑。关键点在于RO_S的触发条件——当SM库满时所有入流变为表层流这是模拟“超渗产流主导”的物理基础。RO_SS和RO_G的分配比例KI和KG并非固定而是随frac可用空间动态缩放体现“蓄满产流”下壤中流/地下流生成能力随土壤湿润程度增强的机制。new_W更新采用经典公式确保产流后土壤含水量同步变化。3. 汇流演算与输出整合如何把三水源结果变成可验证的日流量过程线有了三水源日分量RO_S_arr,RO_SS_arr,RO_G_arr下一步是将它们分别汇流再叠加得到出口断面日流量。课程设计中常采用线性水库汇流或单位线法。由于是日模型单位线通常简化为 23 日的 S 曲线但为降低复杂度我推荐用双线性水库串联表层流快汇、壤中流慢汇、地下流最慢汇。这样既能体现三水汇流差异又避免引入过多参数。3.1 三水源独立汇流为什么不能简单相加而要分别演算直接Q_total RO_S RO_SS RO_G是严重错误。三水源形成位置不同、运移路径不同、响应时间不同表层流几小时到1天内到达河道壤中流需13天地下流则可能滞后数周。若不区分汇流模型将完全失去“三水”意义退化为单水源模型。课程设计要求体现这一物理差异因此必须为每类水源配置独立汇流参数。我采用以下简化方案表层流汇流1 阶线性水库滞时L_s 0.5日系数C_s 0.9壤中流汇流1 阶线性水库滞时L_ss 1.5日系数C_ss 0.7地下流汇流1 阶线性水库滞时L_g 5.0日系数C_g 0.3线性水库演算公式为Q_out[t] C * Q_in[t] (1-C) * Q_out[t-1]注意Q_in[t]是当日产生的水源量mm/dQ_out[t]是当日到达出口的流量mm/d。需为每个水源单独维护一个Q_out序列并初始化首日值如Q_out[0] Q_in[0]。def linear_reservoir(Q_in, C, L, dt1.0): 线性水库汇流离散化 :param Q_in: 输入序列mm/d :param C: 汇流系数0C1 :param L: 滞时日用于初始化延迟 :param dt: 时间步长日默认1 :return: 输出序列mm/d n len(Q_in) Q_out np.zeros(n) # 初始化前 L 天用 Q_in[0] 填充简化处理 init_len int(np.ceil(L)) for i in range(min(init_len, n)): Q_out[i] Q_in[0] * (C ** (i1)) # 主循环Q_out[t] C * Q_in[t] (1-C) * Q_out[t-1] for t in range(init_len, n): Q_out[t] C * Q_in[t] (1 - C) * Q_out[t-1] return Q_out # 为三水源分别汇流 Q_s linear_reservoir(RO_S_arr, C0.9, L0.5) # 表层流 Q_ss linear_reservoir(RO_SS_arr, C0.7, L1.5) # 壤中流 Q_g linear_reservoir(RO_G_arr, C0.3, L5.0) # 地下流 # 叠加得总出口流量mm/d Q_total_mm Q_s Q_ss Q_g # 转回 m³/s用于与实测流量对比 Q_total_m3s Q_total_mm * (area_km2 * 1e6) / 86400 / 1000参数说明C越大响应越快表层流C0.9几乎瞬时L越大初始延迟越长地下流L5.0前5天输出极小。linear_reservoir函数中init_len的处理是课程设计常用简化——不引入复杂单位线但保留了滞时物理概念。实际中L可通过试错调整使模拟峰现时间接近实测。3.2 输出整合与可视化一张图看懂三水贡献占比与过程匹配度课程设计报告的核心图之一就是三水源流量过程线叠加图。它不仅是结果展示更是诊断工具若模拟峰远滞后于实测说明表层流汇流太慢调大C_s或减小L_s若枯季基流衰减过快说明地下流C_g过大或KG过小。import matplotlib.pyplot as plt # 创建时间轴 dates df.index plt.figure(figsize(12, 8)) # 绘制三水源 plt.fill_between(dates, Q_s, label表层流 (mm/d), alpha0.6, colorred) plt.fill_between(dates, Q_s Q_ss, Q_s, label壤中流 (mm/d), alpha0.6, colororange) plt.fill_between(dates, Q_s Q_ss Q_g, Q_s Q_ss, label地下流 (mm/d), alpha0.6, colorblue) # 绘制总流量mm/d plt.plot(dates, Q_total_mm, k-, linewidth2, label总流量 (mm/d)) # 绘制实测流量mm/d plt.plot(dates, df[flow_mm], g--, linewidth2, label实测流量 (mm/d)) plt.xlabel(日期) plt.ylabel(流量 (mm/d)) plt.title(新安江模型三水源划分与汇流过程) plt.legend() plt.grid(True, alpha0.3) plt.xticks(rotation30) plt.tight_layout() plt.show() # 计算三水源总量占比整个模拟期 total_RO np.sum(RO_S_arr RO_SS_arr RO_G_arr) pct_s np.sum(RO_S_arr) / total_RO * 100 pct_ss np.sum(RO_SS_arr) / total_RO * 100 pct_g np.sum(RO_G_arr) / total_RO * 100 print(f三水源总量占比表层流 {pct_s:.1f}%壤中流 {pct_ss:.1f}%地下流 {pct_g:.1f}%)技巧用fill_between而非plot绘制三水源能直观显示各水源对总流量的“堆叠贡献”。百分比计算基于产流量RO_*_arr而非汇流后流量Q_*因为汇流是线性过程不改变总量。这个占比是检验模型物理合理性的第一关——湿润地区地下流占比常超 40%干旱区则表层流主导。4. 避坑指南新安江模型 Python 实现中 4 个血泪教训与排查方法新安江模型看似公式清晰但在 Python 中落地时新手常因细节疏忽导致结果完全失真。以下是我在指导多个模拟项目X时学生反复踩中的 4 个典型坑附带现象、原因与可立即执行的排查命令。4.1 现象模拟流量全程为 0或仅在暴雨日突增巨大尖峰原因RO产流计算中P_eff max(0, P - E)未做导致蒸发大于雨量时P_eff为负进而Def计算错误负缺水量最终RO恒为 0或RO_S分配时未加max(0, ...)保护负值参与后续计算。解决在xaj_daily_step函数开头强制校验P_eff max(0, P - E) # 必须 # 所有水源输出前加保护 RO_S max(0, RO_S) RO_SS max(0, RO_SS) RO_G max(0, RO_G)排查打印前 5 日的P_eff,Def,ROprint(fDay {t}: P{P:.1f}, E{E:.1f}, P_eff{P_eff:.1f}, Def{Def:.1f}, RO{RO:.1f})4.2 现象模拟流量持续增长数月后达数千 mm/d明显发散原因水库状态更新未加约束new_SM或new_EX超过容量SM或EX导致后续日IN_SM被无限放大因frac计算失效形成正反馈发散。解决状态更新后必须截断new_SM max(0, min(new_SM, params[SM])) # 必须上下限 new_EX max(0, min(new_EX, params[EX]))排查监控状态变量极值print(fSM: {state[SM]:.1f} → {new_SM:.1f}, EX: {state[EX]:.1f} → {new_EX:.1f})若发现new_SM params[SM]立即检查min是否漏写。4.3 现象三水源占比异常如地下流占 90% 以上且枯季无基流衰减原因KG参数过小如设为 0.0001导致地下水库出流极慢蓄量持续累积或EX0初始值过大如设为EX使模型从第一天就处于“超饱和”状态。解决KG取值需与EX匹配经验公式KG ≈ 0.005 ~ 0.02EX0设为0.2*EX ~ 0.4*EX。排查单独绘制EX序列plt.plot(df.index, EX_arr); plt.title(地下水库蓄量变化)正常应呈“锯齿状波动”而非单调上升。4.4 现象模拟峰现时间比实测晚 23 天且峰高偏低原因表层流汇流参数C_s过小如 0.5或L_s过大如 2.0导致响应过慢或KI过大将本该为表层流的量分给了壤中流。解决优先调C_s目标C_s ≥ 0.85其次微调KI增大KI会减少表层流慎用。排查关闭壤中流和地下流只运行表层流汇流Q_s_only linear_reservoir(RO_S_arr, C0.85, L0.3)观察其峰现时间是否接近实测峰——这是定位问题的第一步。注意所有排查均需在单日调试模式下进行注释掉主循环手动设t10选一个典型暴雨日逐行打印中间变量。这是比盲目调参高效十倍的方法。5. 参数率定实战用 NSE 和相对误差双指标驱动的手动调参技巧课程设计不要求全自动优化但必须体现“率定”思维——即通过有限次手动调整让模拟结果在统计指标和过程形态上逼近实测。我从不推荐一上来就扔进scipy.optimize因为那会掩盖模型行为。真正的率定是用眼睛看过程、用指标验结果、用物理猜方向的三步闭环。5.1 两个必算指标NSE 效率系数与枯季相对误差NSENash-Sutcliffe Efficiency是水文模型黄金标准定义为NSE 1 - sum((Q_sim - Q_obs)^2) / sum((Q_obs - mean(Q_obs))^2)NSE ∈ (-∞, 1]0.5 为可接受0.75 为好。但它对洪峰不敏感。因此必须补充枯季相对误差取实测流量最小的连续 30 天计算mean(|Q_sim - Q_obs| / Q_obs)。若该值 30%说明基流模拟严重失真需重点调KG和EX。def calculate_metrics(Q_sim, Q_obs): 计算 NSE 和枯季相对误差 # NSE numerator np.sum((Q_sim - Q_obs) ** 2) denominator np.sum((Q_obs - np.mean(Q_obs)) ** 2) NSE 1 - numerator / denominator if denominator ! 0 else -999 # 枯季相对误差找实测最小的连续30天 window 30 min_idx np.argmin(np.convolve(Q_obs, np.ones(window), valid)) Q_obs_dry Q_obs[min_idx:min_idxwindow] Q_sim_dry Q_sim[min_idx:min_idxwindow] # 避免除零 valid_mask Q_obs_dry 0.1 if np.any(valid_mask): rel_err_dry np.mean(np.abs(Q_sim_dry[valid_mask] - Q_obs_dry[valid_mask]) / Q_obs_dry[valid_mask]) else: rel_err_dry np.nan return NSE, rel_err_dry NSE_val, dry_err calculate_metrics(Q_total_m3s, df[flow_m3s].values) print(fNSE {NSE_val:.3f}, 枯季相对误差 {dry_err:.2%})技巧np.convolve用于高效计算滑动窗口和比循环快百倍。valid_mask防止实测枯季流量为 0 时除零错误。这两个指标必须同时看NSE 高但枯季误差大说明模型“靠洪峰撑场面”反之则是“基流准但峰不准”。5.2 手动调参的三阶策略从宏观到微观每次只动一个参数我教学生的口诀是“先调峰再调量最后调尾巴”。具体步骤调峰表层流响应固定其他参数只调C_s0.8→0.9→0.95。目标是让模拟峰现时间与实测对齐。若仍滞后再微调L_s0.3→0.5。每次只改一个改后立刻画图看峰。调量总产流规模峰对齐后若峰高偏低优先调WM减小WM增加产流若峰高偏高调B增大B使蓄水曲线更陡减少产流。WM是最敏感参数±5mm 就可能引起 NSE 变化 0.1。调尾巴枯季基流用枯季误差指导。若误差大调KG增大KG加快地下流出流若调KG后峰变矮说明地下流“抢了”表层流需同步微调KI减小KI让更多水进地下库。血泪经验永远不要同时调WM和SM。WM控制产流总量SM控制三水分配比例二者耦合极强。我一般先用WM75固定调出基本形态再放开WM微调。5.3 一个实用技巧用“参数影响热力图”替代盲目试错与其随机改参数不如做一个轻量级敏感性分析。对每个核心参数在 ±10% 范围内取 5 个值固定其他参数跑一次模型记录 NSE。用seaborn.heatmap可视化import seaborn as sns params_base {WM: 75.0, B: 0.2, SM: 25.0, KG: 0.005} param_names [WM, B, SM, KG] variations np.linspace(0.9, 1.1, 5) # ±10% results np.zeros((len(param_names), len(variations))) for i, pname in enumerate(param_names): for j, ratio in enumerate(variations): p_temp params_base.copy() p_temp[pname] p_temp[pname] * ratio # 运行单次模型略去细节返回NSE nse_temp run_single_simulation(df, p_temp, state_init) results[i, j] nse_temp # 绘图 plt.figure(figsize(10, 4)) sns.heatmap(results, xticklabels[f{r:.1f} for r in variations], yticklabelsparam_names, cmapRdYlBu_r, center0.6) plt.title(参数敏感性热力图NSE值) plt.xlabel(参数缩放比例) plt.ylabel(参数名称) plt.show()效果一眼看出哪个参数对 NSE 影响最大颜色越深越敏感以及最优值在哪个区间。比如若WM行显示 0.95 比例处颜色最深就锁定WM71.25作为下一步起点。这比“凭感觉调”靠谱得多。希望帮到你。这些年带过不少模拟项目X最深的体会是新安江模型的代码本身不难难的是每一行背后都站着赵人俊先生当年在淮河流域实测的几百组数据。我们写的不是 Python是把那些雨滴、土壤、岩层的记忆重新编译成可运行的物理逻辑。调试时多看一眼中间变量就少走十天弯路。本文还有配套的精品资源点击获取