简介一份基于MATLAB实现的水文预报工具包面向水文专业学生、研究人员与工程技术人员用于完成流域日模型与次洪模型的产流预测和洪水过程线估算。资源共21个文件包含17个csv雨量/蒸发等数据文件、3个m脚本文件和1个xls表格压缩包大小仅93KB代码与数据配套完整便于快速运行和二次开发。模型融合泰森多边形法计算面雨量结合Kc作物系数率定修正蒸发输入再通过次洪模型实现降雨到径流的转换预测借助MATLAB脚本可直观查看核心算法实现理解参数调整对预报结果的影响也可替换自有数据集开展定制化模拟。资源已有1300余人学习下载适合具备基础水文知识和MATLAB操作经验的读者用于课程设计、科研实验或工程洪水预报方案的快速搭建与验证。通过研读这套代码能够掌握从降雨数据处理到洪水过程线输出的完整技术链条为后续开展更复杂的分布式水文建模打下基础。 做水文预报的同学应该都有同感论文里产流、汇流公式写得清清楚楚但真要把这套东西用 MATLAB 跑起来出一个能看的日径流过程、场次洪水过程中间隔着无数个“坑”。有人卡在数据结构上有人卡在单位换算更多人是不知道模型代码该怎么组织。这次分享的是一套 matlab 水文预报代码包data.zip包含日模型、次洪模型两部分输入降雨、蒸发、流量资料就能直接率定参数、模拟出流非常适合水文专业研究生做论文、基层防汛人员做预报方案时参考。这套代码的核心思路是把“概念性水文模型”用最朴素的方式落地日模型走连续模拟逐日更新土壤蓄水量、计算蒸发与产流次洪模型聚焦单场暴雨用前期影响雨量刻画初始湿润状态再配合单位线完成汇流。下面我从模型原理、代码结构、实操调参到常见报错把这套东西完整拆一遍。1. 整体设计思路先搞清楚模型要解决什么问题1.1 日模型与次洪模型的本质区别很多初学者会把日模型和次洪模型混为一谈其实两者解决的是完全不同的问题。日模型关心的是“长时段”的水量平衡模拟逐日径流过程常用于年径流模拟、枯季径流预报、水资源评价。它的时间步长是“天”输入是日降雨量、日蒸发能力输出是日径流深或日流量。次洪模型关心的是“单场暴雨”的洪水过程模拟从降雨开始到洪水退尽这一段的流量过程线主要服务于汛期洪水预报。它的时间步长通常是“小时”甚至更短输入是时段降雨输出是时段流量。这个区别直接决定了模型结构的差异。日模型只需要一个能逐日更新状态的“蓄水库”结构土壤含水量、地表径流、地下径流都通过退水系数和蓄泄关系来刻画次洪模型则必须显式处理“产流”和“汇流”两个环节——先由降雨扣除损失得到净雨再把净雨通过单位线或等流时线推到出口断面。代码包里把这两套逻辑分开写正是为了让使用者能按需取用而不是在一个庞大模型里绕来绕去。1.2 为什么用 MATLAB 而不是其他语言水文领域的老传统是 FortranPython 这几年也起来了但 MATLAB 在水文预报里依然有不可替代的位置。最大的优势是“开箱即用”矩阵运算天然适合逐时段循环的模型计算绘图功能极其方便率定时调用 fminsearch 等优化函数也顺手。更重要的是很多高校和设计院的水文计算环境就是 MATLAB研究生拿到的历史代码、师兄师姐留下的模型脚本绝大多数都是 .m 文件。这套代码选 MATLAB 实现意味着你大概率能在实验室的旧电脑上直接跑起来不需要额外搭环境。2. 核心原理拆解日模型和次洪模型分别怎么算2.1 日模型三层蒸发与蓄满产流日模型采用的是概念性水文模型中常见的“蓄满产流”机制也就是著名的张力水蓄水库结构。模型把流域的土壤蓄水量分成三层上层WUM、下层WLM和深层WDM总和为流域平均张力水容量 WM。每一层都有自己的蓄水容量和蒸发折算系数降雨先补充上层上层蓄满后依次向下层、深层渗透。每日的计算顺序是先算蒸发再算产流。蒸发时优先消耗上层张力水上层不够就动用下层下层不够再动用深层这模拟了土壤水的垂直分布规律。产流则遵循“蓄满才产流”的原则——当土壤蓄水量达到田间持水量即 WM时后续降雨全部变成径流未蓄满时则不产流或产流很少。这种机制特别适合湿润半湿润地区比如我国南方的大部分流域。以典型的参数初值为例WUM 取 20mmWLM 取 80mmWDM 取 40mmWM 合计 140mm。蒸发折算系数 K 在 0.8~1.2 之间按季节调整深层蒸发系数 C 通常取 0.1~0.2。这些初值不是拍脑袋定的而是来自长期水文资料分析的统计规律。代码里把这些值都放在了参数数组里你可以按自己流域的实际情况去率定。2.2 次洪模型前期影响雨量与单位线汇流次洪模型的产流环节通常采用初损后损法或降雨径流相关图法但不管哪种方法都绕不开一个关键变量——前期影响雨量 Pa。Pa 反映的是降雨开始前土壤的湿润程度它直接影响产流损失的大小。Pa 的递推公式很简单Pa(t1) K × (Pa(t) P(t))其中 P 是日降雨量K 是消退系数一般在 0.85~0.95 之间干旱地区取小值湿润地区取大值K 由流域蒸发能力和土壤含水特性共同决定。Pa 有一个上限即流域蓄水容量 WM超过后不再增加。汇流环节这里用的是“单位线法”。单位线本质上是一个线性时不变系统的脉冲响应函数——单位时段内、单位净雨深通常取 10mm在流域出口形成的流量过程线。有了单位线任意净雨过程都能通过卷积叠加得到出流过程这就是经典的单位线叠加原理。单位线有三要素洪峰流量、峰现时间、总历时。实际操作中时段的选择很关键——时段太短单位线抖动剧烈模拟结果不稳定时段太长又无法刻画洪峰的陡涨陡落。经验上时段长度取流域汇流时间的 1/5 到 1/3 比较稳妥。代码里提供了单位线输入接口你可以根据流域的实测洪水资料推求也可以先用经验公式如纳希瞬时单位线概化一条。3. 代码包结构解读data.zip 打开后你会看到什么3.1 文件清单与各自职责拿到 data.zip解压后建议先别急着跑花十分钟把目录结构看一遍。一个规范的 MATLAB 水文预报代码包通常长这样文件/文件夹作用main_daily.m日模型主程序负责数据读取、调用模型、绘制结果main_storm.m次洪模型主程序按场次洪水逐场模拟model_daily.m日模型核心函数输入气象数据与参数输出模拟径流model_storm.m次洪模型核心函数输入降雨、Pa、单位线与参数输出流量过程unit_hydrograph.m单位线处理函数负责净雨与单位线的卷积计算calib_daily.m日模型参数率定脚本可调用 fminsearch 自动寻优calib_storm.m次洪模型参数率定脚本read_data.m通用数据读取函数支持 txt、Excel 格式plot_result.m绘图函数输出实测与模拟对比图data/原始资料目录存放降雨、蒸发、流量、面积等数据文件这种“主程序 函数 数据目录”的结构是工程上比较推荐的写法。主程序只负责流程控制模型逻辑全部封装在函数里这样后续调整参数、替换数据都不需要动核心代码。3.2 数据文件的组织格式数据文件一般放在 data/ 子目录下命名要有明确的时间或场次标识。日模型常用的数据文件包括日降雨序列P_daily.txt、日蒸发序列E_daily.txt、实测日径流Q_daily.txt次洪模型则准备场次降雨、场次洪水流量外加一个记录场次起止时间的索引文件。数据格式建议统一为三列或两列第一列是日期或时段序号第二列是降雨量第三列如果有是蒸发量或流量。单位一定要统一模型代码里默认降雨为 mm、蒸发为 mm、流量为 m³/s、时间步长为日或小时。如果不统一比如降雨你给的是 m蒸发给的是 cm模型跑出来就是一团乱麻。这里有个小经验拿到任何水文数据第一步先做单位一致性检查再开始建模。4. 关键代码实现模型函数到底怎么写的4.1 日模型主循环的核心逻辑日模型的函数体其实非常简洁核心就是一个逐日更新的循环。下面是一段精简后的示意代码function Qsim model_daily(P, E, params, W0) % P: 日降雨序列 (mm) % E: 日蒸发能力序列 (mm) % params: [WUM, WLM, WDM, KC, C] % W0: 初始土壤蓄水量 (mm) WUM params(1); WLM params(2); WDM params(3); WM WUM WLM WDM; KC params(4); C params(5); WU 0; WL 0; WD 0; % 初始各层蓄水量可优化 n length(P); Qsim zeros(n, 1); for t 1:n % 蒸发计算上层先蒸发 EM KC * E(t); EU min(WU, EM); WU WU - EU; EL 0; ED 0; if EU EM % 上层不够下层补 EL min(WL, (EM - EU) * WL / WLM); WL WL - EL; end if EU EL EM % 还不够深层按系数 C 补 ED C * WLM * (EM - EU - EL) / WM; WD WD - ED; end % 产流计算蓄满产流 W WU WL WD; if W P(t) WM R W P(t) - WM; % 蓄满部分产流 else R 0; end % 土壤水量分配降雨先补上层再中下层简化示意 WU WU P(t) - R; if WU WUM WL WL (WU - WUM); WU WUM; if WL WLM WD WD (WL - WLM); WL WLM; end end % 汇流简化出流按线性水库退水实际可使用单位线或汇流系数 Qsim(t) R K_route * (W - W0); % 示意实际按模型结构展开 end end代码的逻辑顺序是“先蒸发、后产流、再分水”这正是蓄满产流模型的顺序。注意初始土壤蓄水量 W0 的取值会影响模型前一段的预热过程通常的做法是取 0.5×WM 作为初值让模型先跑一年作为“预热期”预热期的模拟结果不参与率定。4.2 次洪模型的净雨与汇流实现次洪模型的代码有两段核心一段是扣损求净雨另一段是单位线卷积。扣损的简化版本如下% 初损后损法计算净雨 function Rn net_rain(P, Pa0, I0, f, dt) % P: 时段降雨 (mm) % Pa0: 前期影响雨量 (mm) % I0: 初损量 (mm) % f: 后损率 (mm/h) % dt: 时段长 (h) Pa Pa0; Rn zeros(size(P)); loss I0; % 先扣初损 for t 1:length(P) if loss 0 consume min(loss, P(t)); loss loss - consume; P_remain P(t) - consume; else P_remain P(t); end % 后损按 f*dt 扣除 infil min(f * dt, P_remain); Rn(t) P_remain - infil; % 更新前期影响雨量次洪过程中简化处理 Pa min(WM, Pa P(t)); end end单位线卷积计算是次洪模型另一个关键点。单位线法实际上是一个离散卷积假设单位线为 u(j)净雨为 h(i)则出口断面的流量为Q(i) (F / (3.6 × Δt)) × Σ u(j) × h(i-j1)其中 F 是流域面积km²Δt 是时段长h系数 3.6 实现了深度mm到流量m³/s的单位换算。这段卷积在 MATLAB 里可以用 conv 函数直接搞定Qsim conv(unit_h, net_rain) * F / (3.6 * dt); Qsim Qsim(1:length(net_rain)); % 截取有效长度我刚开始用单位线的时候常常忘记乘 F/(3.6×Δt) 这个换算系数结果模拟流量和实测流量差了好几个量级当时还以为是参数问题折腾了很久才发现是单位换算漏了。这个坑必须单独提醒。4.3 参数率定与目标函数手工调参太痛苦代码包里也提供了自动率定接口。常用的目标函数是纳什效率系数 NSENSE 1 - Σ(Qobs - Qsim)² / Σ(Qobs - Qmean)²NSE 越接近 1模拟效果越好小于 0 说明模拟还不如直接用平均值基本属于失败。配合 MATLAB 的 fminsearch 函数可以自动搜索参数。但要注意fminsearch 是局部优化算法对初值敏感所以正确做法是先用人工试错或网格搜索确定一组“比较靠谱”的初值再交给 fminsearch 精调。另外一个容易忽略的率定指标是水量平衡误差 RE也就是模拟总径流量和实测总径流量的相对误差。有时候 NSE 很高但 RE 偏大说明过程形状拟合得好但水量对不上这往往是因为蒸发参数或者初损值设置有问题。建议率定时两个指标一起看。5. 实操全流程从数据读入到出图对比5.1 数据读入与预处理拿到代码包后第一步是把你的流域资料整理成代码能读的格式。推荐用 readtable 函数读取 Excel 或 txt 文件下面是一段读数示范% 读入日降雨、蒸发、流量 data readtable(data/P_daily.txt); P data{:, 2}; % 第二列为降雨 E readtable(data/E_daily.txt); E E{:, 2}; Qobs readtable(data/Q_daily.txt); Qobs Qobs{:, 2}; % 读入流域面积 (km²) area readtable(data/area.txt); F area{1, 1};数据读进来之后最关键的预处理是“对齐时间序列”。不同数据文件的行数必须一致日期必须一一对应。我遇到过好几次数据错位的问题降雨序列从 1961 年开始流量序列从 1961 年 1 月 2 号才开始两者相差一天模型模拟结果怎么看都不对后来逐行检查数据才发现。建议先用 length 核对行数再用 head 查看前几行确认起始日期。5.2 模型运行与绘图运行主程序比较简单核心是把参数传给模型函数然后画图对比。绘图函数的示意如下% 绘制实测与模拟流量对比 t 1:length(Qobs); figure(Color, w); plot(t, Qobs, b-, LineWidth, 1.5); hold on; plot(t, Qsim, r--, LineWidth, 1.5); xlabel(时间序列日/时段); ylabel(流量 (m^3/s)); legend(实测, 模拟, Location, NorthWest); title(水文预报模型模拟效果对比); grid on;画图时我习惯把汛期局部放大再看一段。整体效果看着不错不代表洪峰附近也对得上——NSE 容易被大流量主导小流量段拟合差可能被掩盖。所以建议除了全序列对比图再单独画几个典型场次的洪水对比。5.3 调参的几个方向性经验调参这事新手容易陷入“乱试一通”的误区。我的经验是先看水量平衡——如果模拟总径流量偏小优先调大产流参数如减小 WM 或降低初损如果总径流量偏大优先调大蒸发参数增加 KC。水量对上了再看洪峰形状——洪峰偏低或偏缓调整汇流参数或重新推求单位线退水段拖得太长调大退水系数。参数敏感性也是有规律可循的。在湿润地区产流参数对结果的影响最显著在半干旱地区蒸发参数和初损参数的影响明显上升。所以调参顺序建议是先产流参数再蒸发参数最后汇流参数。不要在汇流参数上花太多时间往往调半天还不如把单位线重新推一次有效。6. 常见问题与排查技巧实录6.1 运行报错的四类典型情况现象可能原因排查方向矩阵维度不匹配降雨、蒸发、流量序列长度不一致检查各数据文件行数模拟流量为 NaN 或 Inf数据中有缺失值或除零用 isnan 检查数据检查分母是否为零NSE 为负或极度偏低参数初值不合理、数据单位错误先检查单位换算再调整参数初值模拟过程整体滞后/超前时间序列错位、单位线时段不匹配核对起始日期检查单位线时段与降雨时段是否一致这里多说一句 NaN 的问题。很多实测水文资料里缺测值用 -99 或 9999 表示readtable 读进来之后如果不处理模型计算时这些异常值会传播导致整条模拟序列报废。最稳妥的做法是在读数后统一做一次清洗P(P 0 | P 500) NaN; % 超过合理范围的降雨置为缺失 P fillmissing(P, linear); % 线性插值填补 Qobs(Qobs 0) NaN; Qobs fillmissing(Qobs, linear);6.2 次洪模型的特殊坑Pa 初值、场次分割、单位线时段次洪模型有几个日模型没有的“专属坑”。第一个是 Pa 初值。如果场次洪水是连续发生的上一场结束时的 Pa 就是下一场开始的 Pa不能每场都从零开始。如果只有单场资料可以用前期的降雨天数来估算连续晴天 15 天以上 Pa 取 0最近几天有降雨则按递推公式计算。第二个是场次分割。洪水过程线的起涨点和退水终止点如何确定直接影响计算成果。起涨点一般取流量开始明显抬升的时刻退水终止点用退水曲线延长法或直接取流量回落到基流的时刻。如果分割不干净把上一场洪水的退水段算进来了模拟效果肯定差。第三个是单位线时段。单位线的时段必须与降雨资料的时段一致如果降雨是 1 小时时段单位线也必须是 1 小时单位线。实际操作中经常遇到降雨资料是 3 小时时段但手头单位线是 6 小时的这时候需要用 S 曲线法进行单位线时段转换不能直接拿过来用。6.3 提高模型稳定性的三个实用习惯最后分享三个我长期用下来觉得特别实用的习惯。第一个是在模型函数开头加参数合法性检查比如 WM 必须是正数、K 必须在 0 到 1 之间一有问题直接报错避免带着错误参数跑出一堆垃圾结果。第二个是每次率定前固定随机种子保证结果可复现。第三个是保留率定过程记录把每次调整的参数、对应的 NSE 和 RE 记下来这能帮你快速判断参数调整方向是否有效而不是凭感觉瞎试。说回这套代码包它的价值不在于“多么先进”而在于“结构清晰、能改能用”。日模型帮你把连续模拟的水量平衡捋顺次洪模型帮你把单场洪水的产汇流算明白两者的数据接口、函数组织方式都是可以沿用到后续研究里的。我自己在多个流域试过类似的代码框架只要数据质量过关、单位统一、参数初值给得合理模拟出来的过程线都能达到实用精度。希望这份拆解能让你少走点弯路跑通代码之后再根据自己的流域特性去改参数、换模块慢慢就能搭出属于自己的预报模型了。本文还有配套的精品资源点击获取
