水文预报课设：用新安江模型与单位线法实现暴雨洪水全过程模拟

简介：本资源是一套面向水文与水资源工程专业本科生的《水文预报》课程设计实践材料，聚焦降雨—径流过程模拟与洪水预报建模核心能力训练，适用于课程设计、毕业设计及水文模型入门实践。压缩包共11个文件，含7个历史年份（1987–1992）降雨数据xlsx文件，用于统计分析降雨时空特征；2个MATLAB脚本（part1.m、part2.m），实现绿-阿姆斯特朗等理论模型的编程构建与参数率定；1份完整课设报告（docx），系统呈现数据处理、模型推导、代码实现与结果验证全过程；另附1张次洪流量过程线jpg图，支撑模型精度评估与洪水规律分析。资源包仅675KB，轻量实用，结构紧凑、模块清晰，便于对照学习与复现。已有645人学习下载，可直接获取从原始数据整理、理论模型应用到可视化验证的全链条实践范例，显著降低水文模型上手门槛。

1. 水文预报课设不是抄模型，而是用真实水文逻辑跑通一场暴雨洪水全过程

很多学生拿到“水文预报课设”作业第一反应是找现成代码改个参数交差，结果在老师追问“为什么选这个产流模型？”“你校验过汇流参数的物理合理性吗？”时哑口无言。实际上，一份合格的水文预报课程设计，核心不是调包跑出一条流量过程线，而是完整复现从降雨输入→流域产流→坡面汇流→河道演进→出口断面预报的全链条水文物理过程。它要求你理解新安江模型的三层蒸散发结构、掌握单位线法中峰现时间与流域坡度的定量关系、能用实测洪峰流量反推马斯京根法中的K和x参数——这些才是水利类专业课设区别于普通编程作业的关键。适合本科高年级或研究生初学水文模型的学生，尤其需要动手调试参数、比对实测与模拟误差、解释偏差成因，而不是仅输出一张漂亮图表。

2. 用新安江模型+单位线法搭建水文预报最小可行系统

水文预报课设最常被选用的技术路径，是将新安江模型（XAJ）作为产流模块，叠加瞬时单位线（IUH）或马斯京根法（Muskingum）作为汇流模块。这种组合兼顾了物理机制可解释性与计算实现可行性：新安江模型能反映湿润地区非线性产流特征，单位线法则避免了复杂河道地形建模。我们不依赖商业软件或黑箱平台，全程使用Python生态实现，确保每一步计算逻辑透明、参数可调、中间结果可查。

2.1 新安江模型产流计算：从降雨到净雨的关键四步

新安江模型的核心在于将降雨分配为蒸发、张力水、自由水三部分，再经自由水蓄满后产流。其关键参数包括流域平均张力水容量Wm、自由水蓄水容量Sm、蒸散发折算系数Kc等。以下是最简产流计算流程（以日尺度为例）：

import numpy as np def xaj_produce_runoff(precip, et, Wm=100.0, Sm=20.0, Kc=1.0, Im=0.1): """ 新安江模型日尺度产流计算（简化版） precip: 日降雨量 (mm) et: 日潜在蒸散发 (mm) Wm: 张力水容量 (mm), 典型值80-120 Sm: 自由水容量 (mm), 典型值15-30 Kc: 蒸散发折算系数, 实测值常取0.8-1.0 Im: 初始张力水相对蓄水量 (0~1), 默认0.1表示初始较干 返回: 净雨量 (mm) """ # 1. 计算实际蒸散发（二层模型简化） W0 = Im * Wm # 初始张力水储量 if et <= W0: Ea = et W1 = W0 - et else: Ea = W0 + (et - W0) * Kc W1 = 0.0 # 2. 降雨入渗与张力水补给 P_eff = max(0, precip - Ea) # 有效降雨 W2 = min(Wm, W1 + P_eff * 0.9) # 张力水层补给（90%入渗率假设） # 3. 自由水生成（蓄满产流） free_water = P_eff * 0.1 + (W1 + P_eff * 0.9 - W2) # 超渗+超蓄部分 R_free = max(0, free_water - Sm * 0.3) # 自由水层产流阈值（取Sm的30%） # 4. 净雨输出（地表径流） runoff = max(0, R_free) return runoff # 示例：某日降雨80mm，ET=3mm，参数取典型值 net_rain = xaj_produce_runoff(precip=80.0, et=3.0, Wm=100.0, Sm=20.0, Kc=0.95) print(f"该日净雨量：{net_rain:.2f} mm")

提示：此代码省略了深层蒸散发与地下径流分摊，但保留了新安江模型最关键的“张力水—自由水”双蓄水结构。Wm决定流域持水能力，Sm控制自由水产流强度，Kc影响干旱期蒸散发抑制程度——这三者是课设中必须手动调试并说明取值依据的参数。

2.2 单位线法汇流：用S曲线反推瞬时单位线

产流得到净雨后，需将其转化为出口断面流量过程。单位线法因其物理意义清晰、参数少、易校准，成为课设首选。我们采用经典的S曲线法由时段单位线（如6小时单位线）推导瞬时单位线（IUH），再与净雨序列卷积得流量过程。

def derive_iuh_from_scurve(duration_hours=6, n_steps=20): """ 由6小时单位线推导瞬时单位线（IUH） duration_hours: 原单位线时段长（小时） n_steps: 时间步数（每步1小时） 返回: IUH序列（长度n_steps），归一化为面积=1 """ # 假设已知6小时单位线（示例数据，实际需根据流域特征拟定） # U6h = [0, 0.1, 0.3, 0.4, 0.15, 0.05] # 6个时段，总和=1.0 U6h = np.array([0, 0.1, 0.3, 0.4, 0.15, 0.05]) # 构造S曲线：累加平移的U6h S_curve = np.zeros(n_steps) for i in range(len(U6h)): shift = i * duration_hours if shift < n_steps: S_curve[shift:] += U6h[i] # IUH = S(t) - S(t-Δt)，Δt=1小时 iuh = np.zeros(n_steps) iuh[1:] = S_curve[1:] - S_curve[:-1] iuh = iuh / iuh.sum() # 归一化 return iuh def convolve_runoff_to_discharge(net_rain_series, iuh, dt_hours=1): """ 净雨序列与IUH卷积得流量过程（m³/s） net_rain_series: 净雨序列（mm），长度N iuh: 瞬时单位线（无量纲），长度M dt_hours: 时间步长（小时） 假设流域面积A=100 km²，产流效率η=0.8 """ A_km2 = 100.0 # 流域面积（km²） eta = 0.8 # 产流效率（考虑下渗损失） # 将净雨(mm)转为径流深(m³/km²) → 再乘面积得总径流体积(m³) # 1mm = 1m³/km² → 总体积 = net_rain * A_km2 * 1e6 (m³) runoff_volume = net_rain_series * A_km2 * 1e6 * eta # m³ # 卷积：Q(t) = Σ runoff_volume[i] * iuh[t-i] Q = np.convolve(runoff_volume, iuh, mode='full')[:len(net_rain_series)] # 转换为流量（m³/s）：体积/时间步长（秒） dt_seconds = dt_hours * 3600 discharge = Q / dt_seconds # m³/s return discharge # 示例：3天净雨序列（每6小时一个值，共12个时段） net_rain_6h = np.array([0,0,5,12,25,30,20,10,5,0,0,0]) # mm iuh = derive_iuh_from_scurve(duration_hours=6, n_steps=24) Q_sim = convolve_runoff_to_discharge(net_rain_6h, iuh, dt_hours=1) print(f"模拟出口流量峰值：{Q_sim.max():.2f} m³/s，出现在第{np.argmax(Q_sim)+1}小时")

注意：单位线参数本质是流域汇流时间的量化表达。duration_hours对应流域响应速度——山区小流域常用3小时单位线，平原大流域用12小时以上。课设中若缺乏实测资料，可依据《水利水电工程设计洪水计算规范》SL44-2006中表3.2.3估算：峰现时间Tp ≈ 0.278 × L / (J^0.5 × v)，其中L为河长（km）、J为河道比降（‰）、v为流速（m/s）。这个公式必须写入课设报告的参数确定依据章节。

3. 用实测水文数据完成参数率定与误差分析

课程设计的价值不在于“跑通”，而在于“证伪与修正”。必须引入真实水文站实测数据（如中国水文信息网公开的长江支流站点日径流数据），通过对比模拟与实测过程线，定量评估模型性能，并针对性调整参数。这是区分优秀课设与应付作业的核心环节。

3.1 获取与预处理实测数据：以长江上游某水文站为例

中国水文信息网（http://www.hydroinfo.gov.cn）提供历史逐日径流数据。以“岷江高场站”为例，下载2020年汛期（6–9月）日流量数据，需进行如下清洗：

import pandas as pd import matplotlib.pyplot as plt # 模拟读取实测数据（实际需从CSV或API获取） # 格式：date, Q_obs (m³/s) data = { 'date': pd.date_range('2020-06-01', '2020-09-30', freq='D'), 'Q_obs': np.random.normal(800, 200, 122) # 占位数据，实际替换为真实值 } df_obs = pd.DataFrame(data) # 添加人工制造的典型暴雨事件（用于重点分析） # 2020-07-15至2020-07-18发生强降雨，实测洪峰1250 m³/s df_obs.loc[(df_obs['date'] >= '2020-07-15') & (df_obs['date'] <= '2020-07-18'), 'Q_obs'] = \ [950, 1120, 1250, 1180] # 生成对应时段的模拟降雨输入（需匹配实测降雨，此处简化为三角形降雨过程） rain_input = np.array([0, 0, 15, 45, 60, 30, 10, 0]) # mm，8个时段（6小时/时段） et_input = np.full_like(rain_input, 3.0) # 日均ET=3mm，按6小时折算为1.5mm/时段 # 运行模型获取模拟流量 net_rain_series = [] for p, e in zip(rain_input, et_input): net_rain_series.append(xaj_produce_runoff(p, e, Wm=95.0, Sm=18.0, Kc=0.92)) iuh_test = derive_iuh_from_scurve(duration_hours=6, n_steps=24) Q_sim_event = convolve_runoff_to_discharge(np.array(net_rain_series), iuh_test, dt_hours=1) # 对齐时间：模拟从2020-07-15 00:00开始，每小时输出 hours_since_start = np.arange(len(Q_sim_event)) df_sim = pd.DataFrame({ 'datetime': pd.date_range('2020-07-15', periods=len(Q_sim_event), freq='H'), 'Q_sim': Q_sim_event })

3.2 五维误差指标体系：不止看RMSE

课设报告中仅列RMSE（均方根误差）是严重不足的。必须构建包含时效性、量级、过程形态的综合评价体系：

def evaluate_forecast(Q_obs, Q_sim): """计算五维误差指标""" # 提取洪峰与峰现时间（需确保序列对齐） Q_obs_peak = Q_obs.max() Q_sim_peak = Q_sim.max() Tp_obs = Q_obs.idxmax() if hasattr(Q_obs, 'idxmax') else np.argmax(Q_obs) Tp_sim = Q_sim.idxmax() if hasattr(Q_sim, 'idxmax') else np.argmax(Q_sim) PE = abs(Q_sim_peak - Q_obs_peak) / Q_obs_peak * 100 TTE = abs(Tp_sim - Tp_obs) # NSE计算（要求长度一致） Q_obs_arr = np.array(Q_obs) Q_sim_arr = np.array(Q_sim)[:len(Q_obs_arr)] # 截断对齐 nse = 1 - np.sum((Q_obs_arr - Q_sim_arr)**2) / np.sum((Q_obs_arr - Q_obs_arr.mean())**2) # 水量平衡误差（假设Δt=1小时） wbe = (Q_sim_arr.sum() - Q_obs_arr.sum()) / Q_obs_arr.sum() * 100 # 相关系数 r = np.corrcoef(Q_obs_arr, Q_sim_arr)[0,1] return {'PE': PE, 'TTE': TTE, 'NSE': nse, 'WBE': wbe, 'R': r} # 执行评估 eval_result = evaluate_forecast(df_obs['Q_obs'].iloc[30:38], Q_sim_event) print("暴雨事件模拟评估结果：") for k, v in eval_result.items(): print(f" {k}: {v:.2f}")

提示：若PE > 20% 或 TTE > 12小时，说明参数严重失配。此时应优先调整Sm（影响洪峰量级）和单位线duration_hours（影响峰现时间），而非盲目修改Wm。课设报告中必须呈现“参数调整→指标变化”对照表，例如：

Sm (mm)	duration (h)	PE (%)	TTE (h)	NSE
15	6	28.3	14	0.42
22	6	12.1	8	0.71
22	4	9.5	3	0.79

4. 马斯京根法替代单位线：当河道几何数据缺失时的稳健选择

并非所有课设都能获取足够精度的流域汇流参数。当缺乏实测单位线或S曲线推导条件时，马斯京根法（Muskingum）是更鲁棒的替代方案——它仅需出口断面实测流量与上游断面流量，即可率定出描述河道演进的两个核心参数K（传播时间）和x（权重系数）。这对仅有上下游水文站数据的课程设计场景尤为实用。

4.1 用线性回归法率定K与x参数

马斯京根方程离散形式为：
Q₂(t) = C₀·Q₁(t) + C₁·Q₁(t−Δt) + C₂·Q₂(t−Δt)
其中C₀、C₁、C₂由K、x、Δt决定。实际操作中，我们直接对历史上下游流量数据做线性回归求解系数，再反推K与x：

def muskingum_calibrate(Q_up, Q_down, dt_hours=6): """ 用线性回归率定马斯京根参数 Q_up: 上游断面流量序列 (m³/s) Q_down: 下游断面流量序列 (m³/s) dt_hours: 计算时段长（小时） 返回: K (小时), x, 以及各系数 """ # 构造回归矩阵：Q_down(t) ~ Q_up(t) + Q_up(t-1) + Q_down(t-1) n = len(Q_up) X = np.column_stack([ Q_up[1:n], # Q_up(t) Q_up[0:n-1], # Q_up(t-1) Q_down[0:n-1] # Q_down(t-1) ]) y = Q_down[1:n] # Q_down(t) # 最小二乘求解 coeffs, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None) C0, C1, C2 = coeffs # 验证系数合理性：C0+C1+C2≈1，且均>0 if not (0.95 < C0 + C1 + C2 < 1.05 and all(c > 0 for c in [C0,C1,C2])): raise ValueError("马斯京根系数不合理，请检查数据质量或时段长") # 反推K与x（公式推导见《水文学原理》P187） K = dt_hours / (C0 + C1) x = (C0 - C1) / (2 * K * (C0 + C1)) if K > 0 else 0.2 return K, x, (C0, C1, C2) # 示例：用上游高场站与下游宜宾站2020年汛期日流量数据 # Q_up = [...] # 高场站日流量 # Q_down = [...] # 宜宾站日流量 # K_est, x_est, coeffs = muskingum_calibrate(Q_up, Q_down, dt_hours=24) # print(f"率定结果：K={K_est:.1f}小时，x={x_est:.2f}")

4.2 马斯京根法在课设中的实操边界

马斯京根法虽易用，但有明确适用前提：河道顺直、断面变化平缓、无显著支流汇入。若课设流域存在大型水库、分汊河道或频繁溃堤，则必须注明模型局限性。常见误用包括：

• 用日尺度数据率定却用于小时尺度预报（K值需按比例缩放）；
• 忽略x参数物理意义（x∈[0,0.5]，x=0为运动波，x=0.5为扩散波）；
• 未验证C₂<0.5（否则数值不稳定，需减小Δt或改用隐式格式）。

课设报告中应明确写出：“本流域河道比降1.2‰，主槽宽深比稳定，无大型支流，符合马斯京根法应用条件；率定K=18.3h，x=0.28，表明以运动波为主导，与实地勘查一致。”

5. 课设成果交付：三张图+一张表构成技术闭环

一份能体现专业深度的水文预报课设，最终交付物绝非冗长文字报告，而是用三张核心图表+一张参数表自洽闭环：它们分别回答“发生了什么”、“为什么这样发生”、“如何证明可信”三个本质问题。这是评审老师快速判断工作量与思考深度的锚点。

5.1 必须包含的三张技术图

1. 降雨-径流过程叠置图：横轴时间，双纵轴（左：降雨mm，右：流量m³/s），实测与模拟流量线重叠，降雨柱状图置于下方。关键标注：洪峰时刻、洪量、起涨点。
2. 参数敏感性热力图：以Sm为横轴（10–30mm）、K为纵轴（10–30h），每个格子填入对应组合下的NSE值，用颜色深浅直观显示最优参数区间。
3. 误差分布直方图：横轴为相对误差（(Qsim−Qobs)/Qobs×100%），纵轴频次，叠加正态分布拟合线。若峰值偏左（负误差多），说明系统性低估，需检查产流效率η。

5.2 不可省略的参数配置表

所有模型参数必须集中呈现在一张表中，并注明来源依据，禁止只写“参考文献[3]”或“经验值”。例如：

注意：表中“依据说明”栏必须具体到数据来源页码、公式编号或实测值出处。若引用《SL44-2006》，需写明“第5.3.2条：x值宜取0.2–0.3，本流域取0.28”。这是课程设计与课程论文的本质分野——前者强调可追溯的工程决策，后者允许理论推演。

课设最后一步，是把模拟流量过程线导入Excel，用“数据验证”功能设置洪峰误差自动标红（PE>15%单元格填充红色），让老师一眼看到你是否真正完成了误差控制闭环。

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