简介：本资源是面向水文与水资源工程专业本科生及水文模型初学者的课程设计实践代码包，聚焦新安江模型日尺度模拟与三水源划分核心任务，完整覆盖地面产流、汇流计算及预报精度评估等关键环节，可直接用于河海大学《水文预报》课程设计作业或模型原理验证。压缩包仅含2个文件（1个Python主程序+1个说明文本），体积精简至6KB，py文件实现自动读取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）。我坚持用以下结构初始化 DataFrame：

import pandas as pd import numpy as np # 读取原始数据（假设文件名为 'input_data.xlsx'） df = pd.read_excel('input_data.xlsx', parse_dates=['date'], index_col='date') # 强制重采样为日频，填充缺失（用前向填充+后向填充兜底） df = df.asfreq('D').fillna(method='ffill').fillna(method='bfill') # 检查并标准化单位：所有输入统一为 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 个核心参数。它们不是随便填的数字，每个都对应明确的水文物理含义：

注意：SM和EX是自由水蓄水库和地下水蓄水库的“最大容量”，不是初始蓄量。初始蓄量需单独设置（如SM0=0.5*SM,EX0=0.3*EX），这是率定时的重要自由度。C和L在纯日模型中常合并为一个“汇流单位线”参数，课程设计常用简化线性汇流：Q_out[t] = C * Q_in[t] + (1-C) * Q_out[t-1]。

2.3 三水源划分的逐日计算逻辑：从产流到分水的完整推导链

三水划分发生在产流之后、汇流之前，是新安江模型最体现经验性的环节。其逻辑链如下（以第 t 日为例）：

1. 计算当日净雨：P_eff[t] = max(0, P[t] - E[t])（忽略前期影响，课程设计常简化）
2. 计算土壤缺水量：Def[t] = WM * (1 - (W[t-1]/WM)^B)（W[t-1] 为昨日土壤含水量）
3. 产流计算：RO[t] = P_eff[t] if P_eff[t] > Def[t] else 0（超渗产流）
4. 自由水蓄水库入流：IN_SM[t] = RO[t] * (1 - IM)（扣除不透水区）
5. 三水源分配：
   • 表层流（地表径流）：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))
6. 更新水库蓄量：
   • 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，下一步是将它们分别汇流，再叠加得到出口断面日流量。课程设计中，常采用线性水库汇流或单位线法。由于是日模型，单位线通常简化为 2～3 日的 S 曲线，但为降低复杂度，我推荐用双线性水库串联：表层流快汇、壤中流慢汇、地下流最慢汇。这样既能体现三水汇流差异，又避免引入过多参数。

3.1 三水源独立汇流：为什么不能简单相加，而要分别演算

直接Q_total = RO_S + RO_SS + RO_G是严重错误。三水源形成位置不同、运移路径不同、响应时间不同：表层流几小时到1天内到达河道，壤中流需1～3天，地下流则可能滞后数周。若不区分汇流，模型将完全失去“三水”意义，退化为单水源模型。课程设计要求体现这一物理差异，因此必须为每类水源配置独立汇流参数。

我采用以下简化方案：

• 表层流汇流：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/d），Q_out[t]是当日到达出口的流量（mm/d）。需为每个水源单独维护一个Q_out序列，并初始化首日值（如Q_out[0] = Q_in[0]）。

def linear_reservoir(Q_in, C, L, dt=1.0): """ 线性水库汇流（离散化） :param Q_in: 输入序列（mm/d） :param C: 汇流系数（0<C<1） :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 ** (i+1)) # 主循环：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, C=0.9, L=0.5) # 表层流 Q_ss = linear_reservoir(RO_SS_arr, C=0.7, L=1.5) # 壤中流 Q_g = linear_reservoir(RO_G_arr, C=0.3, L=5.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越大，响应越快（表层流C=0.9，几乎瞬时）；L越大，初始延迟越长（地下流L=5.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)', alpha=0.6, color='red') plt.fill_between(dates, Q_s + Q_ss, Q_s, label='壤中流 (mm/d)', alpha=0.6, color='orange') plt.fill_between(dates, Q_s + Q_ss + Q_g, Q_s + Q_ss, label='地下流 (mm/d)', alpha=0.6, color='blue') # 绘制总流量（mm/d） plt.plot(dates, Q_total_mm, 'k-', linewidth=2, label='总流量 (mm/d)') # 绘制实测流量（mm/d） plt.plot(dates, df['flow_mm'], 'g--', linewidth=2, label='实测流量 (mm/d)') plt.xlabel('日期') plt.ylabel('流量 (mm/d)') plt.title('新安江模型三水源划分与汇流过程') plt.legend() plt.grid(True, alpha=0.3) plt.xticks(rotation=30) 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,RO：print(f"Day {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(f"SM: {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.02；EX0设为0.2*EX ~ 0.4*EX。
排查：单独绘制EX序列：plt.plot(df.index, EX_arr); plt.title('地下水库蓄量变化')，正常应呈“锯齿状波动”，而非单调上升。

4.4 现象：模拟峰现时间比实测晚 2～3 天，且峰高偏低

原因：表层流汇流参数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, C=0.85, L=0.3)，观察其峰现时间是否接近实测峰——这是定位问题的第一步。

注意：所有排查均需在单日调试模式下进行：注释掉主循环，手动设t=10（选一个典型暴雨日），逐行打印中间变量。这是比盲目调参高效十倍的方法。

5. 参数率定实战：用 NSE 和相对误差双指标驱动的手动调参技巧

课程设计不要求全自动优化，但必须体现“率定”思维——即通过有限次手动调整，让模拟结果在统计指标和过程形态上逼近实测。我从不推荐一上来就扔进scipy.optimize，因为那会掩盖模型行为。真正的率定，是用眼睛看过程、用指标验结果、用物理猜方向的三步闭环。

5.1 两个必算指标：NSE 效率系数与枯季相对误差

NSE（Nash-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_idx+window] Q_sim_dry = Q_sim[min_idx:min_idx+window] # 避免除零 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(f"NSE = {NSE_val:.3f}, 枯季相对误差 = {dry_err:.2%}")

技巧：np.convolve用于高效计算滑动窗口和，比循环快百倍。valid_mask防止实测枯季流量为 0 时除零错误。这两个指标必须同时看：NSE 高但枯季误差大，说明模型“靠洪峰撑场面”；反之，则是“基流准但峰不准”。

5.2 手动调参的三阶策略：从宏观到微观，每次只动一个参数

我教学生的口诀是：“先调峰，再调量，最后调尾巴”。具体步骤：

1. 调峰（表层流响应）：固定其他参数，只调C_s（0.8→0.9→0.95）。目标是让模拟峰现时间与实测对齐。若仍滞后，再微调L_s（0.3→0.5）。每次只改一个，改后立刻画图看峰。
2. 调量（总产流规模）：峰对齐后，若峰高偏低，优先调WM（减小WM增加产流）；若峰高偏高，调B（增大B使蓄水曲线更陡，减少产流）。WM是最敏感参数，±5mm 就可能引起 NSE 变化 0.1。
3. 调尾巴（枯季基流）：用枯季误差指导。若误差大，调KG（增大KG加快地下流出流）；若调KG后峰变矮，说明地下流“抢了”表层流，需同步微调KI（减小KI，让更多水进地下库）。

血泪经验：永远不要同时调WM和SM。WM控制产流总量，SM控制三水分配比例，二者耦合极强。我一般先用WM=75固定，调出基本形态，再放开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], yticklabels=param_names, cmap='RdYlBu_r', center=0.6) plt.title('参数敏感性热力图（NSE值）') plt.xlabel('参数缩放比例') plt.ylabel('参数名称') plt.show()

效果：一眼看出哪个参数对 NSE 影响最大（颜色越深越敏感），以及最优值在哪个区间。比如若WM行显示 0.95 比例处颜色最深，就锁定WM=71.25作为下一步起点。这比“凭感觉调”靠谱得多。

希望帮到你。这些年带过不少模拟项目X，最深的体会是：新安江模型的代码本身不难，难的是每一行背后都站着赵人俊先生当年在淮河流域实测的几百组数据。我们写的不是 Python，是把那些雨滴、土壤、岩层的记忆，重新编译成可运行的物理逻辑。调试时多看一眼中间变量，就少走十天弯路。

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