简介这份MATLAB水文预报代码包面向水文水资源、水利工程方向的研究人员与学生用于流域产流量计算与洪水过程线预测。包内共21个文件以17个csv降雨与蒸发时序数据、3个m脚本含主预报程序与函数模块和1个xls数据表为主压缩包约93KB结构紧凑便于直接运行与二次开发。模型分日模型与次洪模型两部分日模型采用泰森多边形法估算流域面雨量并利用蒸发资料率定作物系数Kc以刻画水分平衡次洪模型则基于降雨径流转换关系推求洪水过程线。读者可借此理解从降雨输入到产流预测的完整算法链路查看具体实现细节或替换自有站点数据定制参数。目前已有1313人学习下载适合作为水文预报建模的入门参考与教学案例。1. 拿到这份 MATLAB 水文模型代码我第一件事是跑通一场历史洪水很多人做水文预报卡住的地方不是公式不会推而是从「公式」到「能跑出流量过程线」之间缺一套能直接运行的工程代码。这份data.zip给的就是这么一套东西MATLAB 版的水文预报日模型和次洪模型前者按天滚动算长期径流后者按场次洪水做短历时模拟。它适合两类人——一类是水文水资源方向的学生需要拿真实数据把产汇流流程走一遍另一类是做防汛调度、水库预报的工程人员想找个能改参数、能换流域的代码底座。我拿到手第一件事不是读代码而是先找一份历史洪水过程把模型跑通看输出和实测对得上多少再回头拆结构。2. 先看清这套代码的骨架日模型和次洪模型到底怎么分工2.1 两个模型不是二选一而是上下游关系日模型和次洪模型在这套代码里是配合使用的。日模型负责长系列连续模拟时间步长是 1 天主要用来做年径流、月径流预报以及为次洪模型提供前期土壤含水量。次洪模型负责单场洪水时间步长通常是 1 小时或几小时用来做洪峰流量、峰现时间的预报。常见做法是先用日模型把整个汛期跑一遍输出每天的土壤含水量、蒸散发、径流深然后把某场洪水开始那天的土壤含水量作为次洪模型的初始条件再驱动次洪模型算逐时流量。这样次洪模型不用从零猜初始状态精度会稳很多。代码里对应的两个入口文件一般命名类似daily_model.m和flood_event_model.m或者放在daily/和flood/两个子目录下。你拿到手先确认哪个是主脚本哪个是函数库。2.2 目录结构数据、模型、工具、输出四块这类项目通常不会只有一个.m文件而是按职责分目录。我拆过的同类代码结构大致是这样data.zip ├── data/ │ ├── rainfall_daily.xlsx % 日降雨列日期、站点、雨量 │ ├── flow_daily.xlsx % 日流量列日期、流量 │ ├── rainfall_hourly.xlsx % 次洪降雨列时间、雨量 │ ├── flow_hourly.xlsx % 次洪流量列时间、流量 │ └── evaporation.xlsx % 蒸发皿或潜在蒸散发 ├── model/ │ ├── daily_model.m % 日模型主入口 │ ├── flood_event_model.m % 次洪模型主入口 │ ├── runoff_generation.m % 产流计算 │ ├── confluence.m % 汇流计算 │ └── parameters.m % 参数定义与默认值 ├── utils/ │ ├── read_data.m % 读 Excel/CSV │ ├── nse.m % Nash-Sutcliffe 效率系数 │ ├── rmse.m % 均方根误差 │ └── plot_process.m % 画过程线 └── output/ ├── daily_result.mat └── flood_result.mat这个结构的好处是数据、模型、工具、输出分开你换流域只需要改data/里的文件换算法只需要动model/里的函数不会牵一发动全身。2.3 数据格式时间列必须能解析缺测要提前处理MATLAB 读 Excel 用readtable最省事但有两个坑一是日期列可能被读成字符串或数字二是缺测值可能是NaN、-99或空字符串。我一般会在读数据后加一段清洗% read_data.m 片段读日降雨并清洗 function T read_daily_rain(file) T readtable(file); % 日期列统一转 datetime兼容 Excel 序列号和字符串 if isnumeric(T.date) T.date datetime(T.date, ConvertFrom, excel); else T.date datetime(T.date); end % 雨量缺测统一转 NaN-99 和空值都算缺测 T.rain(T.rain -99 | isnan(T.rain)) NaN; % 按日期排序防止原始表乱序 T sortrows(T, date); end逻辑说明readtable自动识别列名但日期格式因 Excel 版本而异所以先判断是不是数字是数字就按 Excel 序列号转雨量缺测用-99是水文数据常见约定统一转NaN后后续插补或跳过都方便。参数file是文件路径返回的T是排序后的 table。2.4 参数文件所有可调参数集中放别散落在各函数里我见过太多代码把参数写死在函数内部改一个值要翻五个文件。这套代码如果做得规范应该有一个parameters.m把日模型和次洪模型的参数都列出来% parameters.m集中管理参数 function p parameters() % 日模型参数 p.daily.K 0.8; % 蒸散发折算系数 p.daily.WM 120; % 张力水容量mm p.daily.B 0.3; % 蓄水容量曲线指数 p.daily.IM 0.02; % 不透水面积比例 p.daily.SM 30; % 自由水蓄量mm p.daily.KG 0.4; % 地下水出流系数 p.daily.KI 0.3; % 壤中流出流系数 % 次洪模型参数 p.flood.n 2.0; % 纳什汇流参数 n p.flood.K 3.5; % 纳什汇流参数 K小时 p.flood.CS 0.6; % 地表径流消退系数 end逻辑说明p是结构体日模型和次洪模型各一个子结构。WM、SM这类参数有物理意义改的时候要对照流域特征n、K是汇流参数一般用实测洪水率定。把参数集中后主脚本开头一句p parameters();就能拿到全部默认值后面想批量调参也方便。3. 跑通日模型从降雨蒸发到逐日流量过程3.1 主脚本流程读数据、算产流、算汇流、评精度日模型的主脚本一般不长核心就是按时间步循环。我把它拆成四步% daily_model.m日模型主流程 clear; clc; p parameters(); % 加载参数 T read_daily_rain(data/rainfall_daily.xlsx); E read_daily_evap(data/evaporation.xlsx); Qobs read_daily_flow(data/flow_daily.xlsx); % 对齐时间轴取三者共同日期 [t, R, Ep, Qo] align_series(T, E, Qobs); n length(t); Qsim zeros(n, 1); W p.daily.WM * 0.5; % 初始土壤含水量取半容量 S p.daily.SM * 0.5; % 初始自由水蓄量 for i 1:n % 产流算径流深 R 和更新土壤含水量 [Rs, W, S] runoff_generation(R(i), Ep(i), W, S, p.daily); % 汇流把径流深转成流量 Qsim(i) confluence_daily(Rs, p.daily); end % 精度评价 NSE nse(Qo, Qsim); RMSE rmse(Qo, Qsim); fprintf(日模型 NSE %.3f, RMSE %.3f\n, NSE, RMSE); plot_process(t, Qo, Qsim, 日模型);逻辑说明align_series把降雨、蒸发、流量按共同日期对齐避免时间错位。循环里先产流后汇流W和S是状态变量每一步更新后传给下一步。nse和rmse是评价函数输出精度指标。参数p.daily控制产汇流行为。3.2 产流函数蓄满产流 自由水蓄量产流是日模型的核心。常见做法是蓄满产流即土壤含水量达到田间持水量后才产流。代码里一般这样写% runoff_generation.m蓄满产流 自由水蓄量 function [Rs, W, S] runoff_generation(P, Ep, W, S, p) % 蒸散发先算实际蒸发不超过可用水量 Ea min(Ep * p.K, W); W W - Ea; % 降雨补充土壤水 W W P; % 超过张力水容量的部分产流 if W p.WM Rs W - p.WM; W p.WM; else Rs 0; end % 自由水蓄量再分配一部分直接径流一部分壤中流一部分地下径流 S S Rs; if S p.SM excess S - p.SM; S p.SM; Rs excess; % 超过自由水容量的部分直接产流 else Rs 0; end end逻辑说明Ea是实际蒸发受p.K折算W先减蒸发再加降雨超过p.WM的部分就是径流深Rs。自由水蓄量S再分配超过p.SM的部分直接产流。参数p.K、p.WM、p.SM决定产流多少。3.3 汇流函数线性水库叠加日模型汇流常用线性水库把径流深按消退系数逐日释放% confluence_daily.m线性水库汇流 function Q confluence_daily(Rs, p) persistent storage if isempty(storage) storage 0; end storage storage Rs; Q p.KG * storage; % 地下水出流 storage storage - Q; end逻辑说明persistent让storage在函数调用间保持模拟水库蓄量。p.KG是出流系数越大退水越快。注意这个简化版只算地下水完整版还会加壤中流和地表径流。3.4 精度评价NSE 和 RMSE 怎么看NSE 越接近 1 越好大于 0.7 算合格大于 0.8 算好。RMSE 是均方根误差单位是流量单位越小越好。代码里% nse.mNash-Sutcliffe 效率系数 function v nse(obs, sim) obs obs(:); sim sim(:); v 1 - sum((obs - sim).^2) / sum((obs - mean(obs)).^2); end逻辑说明obs和sim是列向量mean(obs)是实测均值。NSE 对洪峰误差敏感如果洪峰差太多NSE 会明显下降。4. 次洪模型实战一场洪水的逐时模拟与参数率定4.1 次洪模型主流程用日模型输出做初始条件次洪模型的主脚本和日模型类似但时间步是小时且初始土壤含水量从日模型结果里取% flood_event_model.m次洪模型主流程 clear; clc; p parameters(); % 读次洪降雨和流量 T read_hourly_rain(data/rainfall_hourly.xlsx); Qobs read_hourly_flow(data/flow_hourly.xlsx); % 从日模型结果加载初始土壤含水量 load(output/daily_result.mat, W_end, S_end); W W_end; S S_end; n length(T.rain); Qsim zeros(n, 1); for i 1:n [Rs, W, S] runoff_generation(T.rain(i), 0, W, S, p.daily); Qsim(i) confluence_flood(Rs, p.flood); end NSE nse(Qobs, Qsim); fprintf(次洪 NSE %.3f\n, NSE); plot_process(T.time, Qobs, Qsim, 次洪模型);逻辑说明W_end和S_end是日模型最后一步的状态作为次洪模型初始条件。次洪期间蒸发忽略所以Ep传 0。confluence_flood用纳什瞬时单位线汇流。4.2 纳什汇流n 和 K 怎么调纳什汇流有两个参数n是线性水库个数K是每个水库的滞时。代码里% confluence_flood.m纳什瞬时单位线汇流 function Q confluence_flood(Rs, p) persistent UH if isempty(UH) % 预计算纳什单位线长度取 5*n*K t 0:1:(5 * p.n * p.K); UH (1 / (p.K * gamma(p.n))) .* (t / p.K).^(p.n - 1) .* exp(-t / p.K); UH UH / sum(UH); % 归一化保证水量平衡 end % 卷积径流深与单位线卷积得流量 Q sum(Rs * UH(1:min(end, numel(UH)))); end逻辑说明gamma(p.n)是伽马函数UH是纳什单位线归一化后保证总水量为 1。n越大过程线越对称K越大洪峰越滞后。调参时先调K对峰现时间再调n对过程线形状。4.3 参数率定手动试算 自动优化手动试算就是改parameters.m里的值跑一遍看 NSE。自动优化可以用fminsearch% calibrate_flood.m用 fminsearch 率定 n 和 K function calibrate_flood() Qobs read_hourly_flow(data/flow_hourly.xlsx); R read_hourly_rain(data/rainfall_hourly.xlsx); obj (x) -nse(Qobs, run_flood(R, x(1), x(2))); x0 [2.0, 3.5]; xopt fminsearch(obj, x0); fprintf(最优 n %.3f, K %.3f\n, xopt(1), xopt(2)); end逻辑说明obj是负 NSE因为fminsearch求最小。x0是初值一般取经验值。run_flood是封装好的次洪模拟函数输入降雨和参数输出模拟流量。4.4 结果对比洪峰、峰现时间、径流深次洪模型评价不只看 NSE还要看洪峰相对误差、峰现时间误差、径流深相对误差。代码里可以加% 洪峰误差 [Qmax_obs, i_obs] max(Qobs); [Qmax_sim, i_sim] max(Qsim); peak_err (Qmax_sim - Qmax_obs) / Qmax_obs * 100; % 峰现时间误差小时 time_err i_sim - i_obs; fprintf(洪峰误差 %.2f%%, 峰现时间误差 %d 小时\n, peak_err, time_err);逻辑说明max返回洪峰值和位置peak_err是相对误差time_err是时间步差。一般要求洪峰误差在 20% 以内峰现时间误差在 1-2 个时间步内。5. 避坑与排查跑水文模型最容易翻车的五个地方5.1 现象NSE 是负数过程线完全对不上原因时间轴没对齐降雨和流量错位一天或几小时。解决用align_series强制按共同时间戳合并打印前 10 行检查。5.2 现象洪峰模拟偏高或偏低很多原因产流参数WM或SM不合适或者汇流参数K太大导致峰值被削平。解决先固定汇流参数调产流参数让径流深接近实测再调汇流参数对洪峰和峰现时间。5.3 现象日模型跑着跑着土壤含水量爆掉原因蒸发或降雨缺测被当成 0 或极大值。解决读数据时把缺测转NaN循环里遇到NaN就跳过或插补别直接参与计算。5.4 现象次洪模型初始土壤含水量不知道取多少原因没跑日模型或者日模型结果没保存。解决先跑日模型把洪水开始前一天的W和S存到daily_result.mat次洪模型直接load。5.5 现象fminsearch率定结果不稳定每次跑都不一样原因初值太远或者目标函数有多个局部极小。解决先用经验值手动试几组选 NSE 最高的那组做初值或者改用patternsearch加边界约束。6. 二次开发换流域、换算法、接实时数据的改造思路6.1 换流域只改数据文件和参数换流域不需要动模型代码只需要把data/里的 Excel 换成新流域的然后重新率定parameters.m里的参数。我一般会建一个params/目录每个流域一个参数文件主脚本开头用p parameters(流域名);加载。这样同一套代码可以管多个流域。6.2 换算法把产流和汇流函数替换掉如果你想试新安江模型以外的算法比如 SCS-CN 或 TOPMODEL只需要重写runoff_generation.m和confluence.m保持输入输出接口不变。主脚本不用改。接口约定产流函数输入(P, Ep, W, S, p)输出(Rs, W, S)汇流函数输入(Rs, p)输出Q。6.3 接实时数据把读文件改成读数据库或 API实时预报场景下数据不是 Excel而是数据库或接口。改造点在read_data.m把readtable换成fetch或webread返回同样格式的 table 即可。注意加异常处理接口超时或返回空时要能降级。6.4 批量率定用循环遍历多场洪水单场洪水率定容易过拟合我一般会选 5-10 场洪水一起率定目标函数改成平均 NSE% 多场洪水平均 NSE function v mean_nse(events, x) v 0; for i 1:numel(events) Qsim run_flood(events(i).R, x(1), x(2)); v v nse(events(i).Q, Qsim); end v v / numel(events); end逻辑说明events是结构体数组每场洪水有R和Q。mean_nse返回平均 NSEfminsearch最大化它。这样率定的参数泛化能力更好。6.5 一个我踩过的坑别在循环里反复读文件早期我写代码在时间循环里每一步都readtable读一次降雨结果跑一次要十几分钟。后来改成循环外一次性读入循环内只索引速度提升几十倍。从那以后我每次写循环都强制检查有没有重复 IO、有没有重复计算。希望帮到你。本文还有配套的精品资源点击获取
