简介一个基于MATLAB的水文预报代码包涵盖日模型与次洪模型两大模块面向水文学研究者、水利工程师及高校相关专业学生用于流域产流量预测与洪水过程线模拟。压缩包共21个文件包含17个CSV数据文件1989—2003年多年度实测降雨、蒸发及单位线数据、3个MATLAB脚本主程序、功能函数与辅助计算及1个Excel数据文件整体仅93KB轻量易部署。目前已有1313人浏览/学习。代码完整实现了从泰森多边形法估算面雨量、Kc作物系数率定到次洪过程模拟的关键算法用户可直接运行查看输出也可替换自定义数据集进行二次开发借助自带多年实测数据与模块化程序还能快速理解水文预报建模逻辑、掌握参数调整方法整个模型覆盖了从降雨输入、产流计算到洪水过程线输出的完整预报链路。1. 收到一份Matlab版水文预报日模型、次洪模型代码包先别急着解压压缩包名叫data.zip听起来像普通数据文件实际里面装的是两套时间尺度完全不同的水文预报模型代码日模型按天演算产汇流给出逐日径流过程次洪模型按小时甚至分钟盯一场暴雨到一场洪水的完整响应。很多人拿到这类包的第一反应是解压后直接跑主函数结果要么报错要么模拟过程和实测对不上。问题往往不在代码而在没分清两套模型的运行逻辑有什么区别。这篇文章按一线处理这类代码包的顺序往下拆先看目录和数据格式再分别跑通日模型和次洪模型最后给出一份参数设置和排错记录。适合刚拿到类似压缩包、想快速复现结果的从业者也适合想从代码层面完整理解两类模型差异的同学。2. data.zip里有什么目录结构、数据格式与检查清单2.1 先看目录和主函数模型入口在哪解压之后先别急着点开m文件。常见做法是压缩包按功能拆成几个子目录日模型、次洪模型、数据、参数、结果。日模型目录里通常有一个主函数负责读数据、调产流和汇流子函数、输出模拟流量旁边挂着两三个子函数比如产流计算、蒸散发计算、分水源计算。次洪模型目录里一般按场次组织数据文件加上净雨计算和单位线汇流两个核心函数。data子目录放样例数据有的包会把日模型和次洪模型的样例数据分开放。拿到代码包后我一般先做三件事顺序固定先找README或说明文件读三行就能知道作者预期的操作路径再找主函数入口看名字带main或run的m文件确认它的输入输出最后确认数据文件格式是mat、txt还是xlsx。这三件事决定了后面所有操作。很多代码包不是不能跑而是入口文件和数据格式对不上导致第一步就报错后面无从下手。如果你打开压缩包发现没有README也别慌。直接搜索所有m文件中体积最大、名称含main或run的那个通常就是主函数。用Matlab打开后在编辑器里看函数声明行的输入参数对照data目录下的文件名基本能推断出数据怎么喂进去。目录/文件常见作用重点看什么daily/日模型主函数与子函数主函数输入输出、蒸散发计算方式storm/次洪模型主函数与子函数单位线如何保存、净雨计算在哪一步data/样例降雨、蒸发、流量数据单位、时间步长、缺测标记params/模型参数文件参数是否被主函数引用、量级是否合理result/模拟输出与对比图输出格式、是否含率定统计量2.2 日模型与次洪模型的数据格式差异这是整个代码包里最容易让人翻车的地方。日模型和次洪模型的数据格式几乎不共享。日模型需要的是逐日序列日雨量、日蒸散发能力、日实测流量用于率定再加上一个代表初始状态的土壤含水量初值。次洪模型需要的是单场洪水的时段雨量序列和时段长、降雨起止时间、前期土壤含水量初值以及这场洪水的实测流量过程用于对比。单位不一致是第一个真正的坑。日雨量常用mm/d日流量用m3/s蒸散发能力有时给mm/d有时给的是折算系数次洪模型的时段雨量可能是mm/时段也可能已经换算成mm/h。代码里如果没有显式的单位换算模拟结果和观测值经常差一个数量级。碰到陌生数据包先确认单位再谈跑通这句话值得贴在屏幕边上。除了单位还有个容易被忽略的差异时间轴的对齐方式。日模型要求降雨、蒸发、流量严格按同一天对齐缺一天都可能导致全序列错位次洪模型则要求时段雨量和实测流量用同一个起始时刻起涨点差一两个小时峰现时间误差直接不合格。2.3 跑代码前先过一遍数据检查清单动手之前把样例数据逐项过一遍能省下后面大量排错时间。我每次处理这类代码包都按固定清单查缺测值标记是什么。NaN、-999还是空行直接影响代码里读数据的逻辑。时间排序是正序还是倒序。有些老数据文件从后往前排读进来之后产流计算顺序完全颠倒。降雨、蒸发、流量是否在同一时间轴上。日模型常见的问题是雨量和蒸发差一天因为一个是北京时间8时起算一个是0时起算。次洪模型数据里有没有包含前期土壤湿度信息比如前期降雨指数API值。数据的空间尺度。单位是mm还是万m3有些数据包给出的已经是面平均雨量有些给的是单站雨量两者产流计算结果差异很大。这套检查并不需要运行代码用Excel打开数据文件扫一遍就够。关键是别跳步。数据格式没问题后面跑代码才是真正的模型问题数据格式有问题后面所有结果都是无效的白调参数。3. 跑通日模型代码产流机制、参数设置与最小示例3.1 日模型在算什么蓄满产流与逐日水量平衡日模型本质上是一个逐日水量平衡账号。每天的雨量进入流域一部分补充土壤含水量一部分变成径流同时蒸散发从土壤含水层里扣水。土壤含水量达到蓄水容量后后续降雨基本全部产流这就是蓄满产流的思路。更完整的实现会加自由水层、分水源、汇流演算但核心还是这个每日水量平衡循环。为什么用日步长而不是小时步长因为日模型服务的是中长期趋势旬月径流预报、水库调度、气候情景评估。它不追求单场洪水的峰现时间精度而是追求连续数旬、数月的总量和过程趋势一致。所以日模型的率定目标通常选纳什效率系数NSE但NSE之外还要看累计径流误差两个指标一起评才能说明模型真的抓住了流域的水量特征。蒸散发计算是个容易被忽略但影响很大的模块。有的代码直接用实测蒸发皿数据乘折算系数有的用气温公式估算。不管是哪种都要确认蒸散发的输入序列和降雨序列时间对齐。很多日模型模拟结果夏半年偏低、冬半年偏高问题基本出在蒸散发输入上而不是产流参数。3.2 日模型最小运行流程与代码下面这段代码是日模型产流循环的最小骨架代表了这类代码包里最常见的核心逻辑先算土壤含水量补给再算产流最后扣蒸散发。实际代码包里的实现会更复杂但看懂这一段主函数就能看懂一半。% 日模型主循环逐日计算产流量 % 输入: P(日雨量), Ep(日蒸散发能力), W0(初始张力水), WM(张力水容量) % 输出: R(逐日产流量), W(逐日张力水) function [R, W] daily_runoff(P, Ep, W0, WM) n length(P); R zeros(n, 1); W zeros(n, 1); W(1) W0; for t 1:n % 先补给张力水降雨先补充土壤缺水超出部分才产流 dW WM - W(t); % 当前缺水量 if P(t) dW R(t) P(t) - dW; % 满足缺水量后剩余降雨成为产流 W(t) WM; else W(t) W(t) P(t); % 降雨全用于补充土壤水 end % 蒸散发按张力水比例折算实际项目中会用蒸发皿或气温估算 if t n W(t) max(0, W(t) - Ep(t) * W(t) / WM); end end end这段代码的逻辑很简单每天先看土壤还有多少缺水容量雨大就把缺口填满剩余部分直接产流雨小就全存进土壤里。蒸散发按当前土壤含水量占容量比例扣减含水量越高蒸发扣得越多这个近似在多数概念性模型里都成立。WM在这里是最关键参数它控制着一个流域能存多少水直接影响产流总量和退水快慢。跑通整个日模型一般需要四行代码从加载数据到出图对比% 日模型最小运行流程 load(data/daily_basin.mat); % 样例数据包含P, E, Q par load(params/daily_default.mat); % 参数文件包含WM等 % 调用主函数输入参数结构体、降雨、蒸发、初始土壤含水量 [R, W, Qsim] daily_model(par, P, E, Q(1)); % 画图对比模拟与实测 figure; plot(Q, k-, LineWidth, 1); hold on; plot(Qsim, r--, LineWidth, 1); legend(实测流量, 模拟流量); xlabel(日序号); ylabel(流量);如果数据是txt或xlsx把load换成readtable或readmatrix即可但要注意列名和单位转换。主函数返回的Qsim单位要和田实测Q一致常见做法是代码内部把径流深算完后再乘流域面积换成流量如果Qsim和Q差了好几个数量级优先查这一步。3.3 必调参数表与率定策略日模型的核心参数并不算多真正需要反复调的是下面几个。拿到一个流域先按量级给初值再逐个调不要一上来全部一起动。参数含义常见量级调参影响WM张力水容量100~250 mm控制产流总量和退水速度偏大则产流偏少K蒸散发折算系数0.8~1.2影响无雨期退水速率偏大则枯季流量偏低C深层蒸散发系数0.05~0.2影响枯季基流和连续干旱期的土壤水消耗X蓄水容量曲线方次0.2~0.5控制产流面积随土壤含水量的非线性变化率定顺序我一般这样来先固定X和C调WM和K让过程线整体形状接近再调X改善洪峰的陡峭程度最后调C修正枯季基流。每一步只动一个参数改完看NSE和累计径流误差两个指标。% 纳什效率系数计算日模型和次洪模型通用 NSE 1 - sum((Qsim - Qobs).^2) / sum((Qobs - mean(Qobs)).^2);NSE接近1说明模拟和实测高度一致接近0说明模型只比用均值预测好一点负值说明模型完全不可信。但NSE对洪峰敏感对枯季低流量不敏感所以还要同时看累计径流误差误差在正负10%以内算合格。自动率定不是不能用但直接用遗传算法整包调参很容易陷入过拟合调出来的参数换一场洪水就失效。人工先调一轮再用自动算法小范围精修是最稳妥的路径。4. 次洪模型代码实战从时段雨量到洪水过程线4.1 次洪模型与日模型的本质差异次洪模型和日模型看起来都在做产汇流但目标完全不同。日模型关心年度汛期和非汛期趋势时间步长一天模拟结果对单场暴雨不敏感次洪模型关心一场洪水的洪峰流量、峰现时间和洪水过程线形状时间步长最小可以到10分钟。同样的降雨输入日模型算出来的总产流量可能合理但过程线形状和峰现时间完全不能用这就是为什么两个模型的代码要分开看。最本质的差异在于初始状态的敏感度。次洪模型模拟的那几个小时里流域初始土壤含水量决定了这场洪水能产多少流。同样一场50mm的降雨如果前期土壤已经饱和洪峰可能是400m3/s如果前期干旱洪峰可能只有80m3/s。所以次洪模型的代码里一定有初始含水量参数而且这个参数对结果的敏感度排在所有参数最前面。另一个差异是产流机制。次洪模型更常采用初损后损的思路降雨开始后先满足初损之后超过下渗能力的部分才成为净雨。日模型里用张力水容量平滑处理的事次洪模型要按小时逐个时段判断代码写法也因此不同。4.2 净雨计算与单位线汇流的Matlab实现下面这段代码是次洪模型最常见的处理流程先扣初损后损把总降雨变成净雨过程再用单位线把净雨推成流量过程。单位线本质是流域对单位净雨的响应函数线性叠加后得到整个洪水过程。% 次洪模型净雨计算与单位线汇流 % 输入: Pdt(各时段雨量), dt(时段长小时), Ia(初损), f(稳渗率) % UH(单位线长度m), dtUH(单位线时段长小时) % 输出: Qs(模拟流量过程线) function Qs storm_flow(Pdt, dt, Ia, f, UH, dtUH) n length(Pdt); % 第一步扣初损后损得到时段净雨 In zeros(n, 1); cumP 0; for t 1:n cumP cumP Pdt(t); % 累计降雨量 if cumP Ia In(t) 0; % 初损未扣完不产流 else In(t) max(0, Pdt(t) - f * dt); % 超过稳渗率的部分成为净雨 end end % 第二步净雨经单位线汇流线性叠加 m length(UH); Qs zeros(n m - 1, 1); for t 1:n if In(t) 0 % 每个时段净雨单独乘单位线叠加到对应时段上 Qs(t:tm-1) Qs(t:tm-1) In(t) * UH * dt / dtUH; end end end这段代码有两个细节决定成败。第一初损是按累计降雨量判断的不是按每个时段独立判断否则同一场雨的初损会被重复扣除。第二单位线的时段长可能和降雨时段长不一致代码里用dt/dtUH做换算很多实现在这一步偷懒结果峰现时间全偏。净雨为负数的时段要直接置零否则后续叠加会把流量越算越小。单位线的获取方式在代码包里通常有两种一种是从实测单峰洪水反推另一种是代码内建的经验公式。反推单位线需要用孤立单峰洪水资料过程比较长至少要拿到一场净雨时段少、洪峰独立的实测洪水才能推得稳。如果代码包里直接给了UH向量先确认它的时段长再和雨量数据对齐。下面用一个简化例子看整个过程。某流域一场降雨分6个时段每时段10mm初损15mm稳渗率1mm/h时段长1小时单位线为[0.2, 0.5, 0.3]单位m3/s per mm时段雨量(mm)累计(mm)净雨(mm)单位线叠加后流量(m3/s)11010002102091.83103096.34104097.25105096.36106095.4这个例子说明初损只扣一次后续时段只要雨强超过稳渗率就产生净雨。单位线长度决定洪水过程拉多长如果实测洪峰偏晚看单位线峰值位置是否合理如果洪峰偏大优先查稳渗率和初损。4.3 次洪模型的参数与率定目标参数含义常见量级敏感度调整方向Ia初损5~30mm高影响起涨时间和洪峰相位f稳渗率0.5~5mm/h中影响洪峰大小偏大则洪峰偏小UH单位线由实测推求高影响峰现时间和过程形状W0初始土壤含水量0~Wmax极高影响总产流量误差被直接放大次洪模型的率定目标不能只用NSE。NSE对整体过程敏感但对洪峰和峰现时间的约束不够。我一般会同时看三个数峰现时间误差以小时计、洪峰相对误差、NSE。目标设定为峰现时间误差不超过1个时段洪峰相对误差控制在20%以内NSE大于0.7。三段目标全部满足才认为这场洪水模拟合格。自动率定在次洪模型上要更谨慎。单位线本身就带有很强的不确定性再把初损和稳渗率一起交给优化算法很容易得到一组数值很漂亮但物理意义完全说不通的参数。先固定单位线调Ia和f再微调单位线这个顺序能避开大部分过拟合问题。5. 水文预报模型常见问题与避坑记录5.1 流量单位与产流面积换算差一个数量级的翻车现场现象日模型算出来的模拟流量和实测流量过程线形状完全一致但数值整体差了100倍看起来像是模型产流太少。 原因代码里降雨单位是mm产流结果也是mm但最后转流量时忘记乘流域面积。或者反过来实测流量是m3/s模拟流量是mm/d两者没有统一量纲就画在一张图上。 解决先把流域面积确认清楚单位是km2还是m2再检查主函数里有没有产流深转流量的换算行一般是用面积乘以产流深再除以86400秒。没有换算行就手动补上换算后两条曲线在同一量级再谈率定。这个坑在日模型和次洪模型各自的主函数里都出现过排错时要先看输出变量注释。5.2 次洪模型初始状态第一场洪水为什么经常偏低现象多场洪水的率定结果都很好唯独把模拟范围扩展到第一场新洪水时模拟流量系统性偏低过程线起涨迟缓。 原因初始土壤含水量被设成了一个固定经验值没有考虑这场洪水前的实际干旱程度。前期土壤越干这场雨被土壤吸收的越多产流越少。固定初值的模型等于默认每一场洪水前期条件相同这在水文上站不住脚。 解决在运行次洪模型前单独用一个预热期计算初始土壤含水量取这场洪水前15到30天的降雨和蒸发序列用日模型或一个简单的水平衡公式推算出W0。如果没有预热数据至少用前期降雨指数API来估算代码里把这个值作为参数传入而不是写死一个常数。5.3 蒸散发数据缺失的替代方案与季节偏差现象日模型在无雨期模拟的退水过程比实测陡很多夏季尤其明显连续无雨20天后模拟流量几乎归零实测还在缓慢退水。 原因蒸散发输入用的是多年平均固定值或者干脆用了一个粗糙的常数没有反映夏季高蒸发、冬季低蒸发的季节变化。退水阶段水量平衡完全由蒸散发控制输入偏大退水就偏快。 解决优先用实测蒸发皿数据乘折算系数没有实测资料时用气温法估算逐日蒸散发能力Hamori公式这类简单方法在温带流域的精度够用。如果代码不支持外部蒸散发序列就找代码里蒸散发计算函数入口把固定值替换成读入序列这个改动量不大但能把退水过程从毛糙变平滑。5.4 单位线时段长与雨量时段长不一致过程线变形现象次洪模型模拟的洪水过程线形状看着合理但整体向左平移峰现时间比实测早了2个小时换一场洪水后可能又偏晚。 原因单位线是从特定时段长比如1小时的雨洪资料推出来的而当前数据文件的时段长是30分钟或2小时代码直接拿来用没有做单位线时段转换。单位线代表的是流域响应换时段长必须重新采样或做S曲线转换。 解决看代码里是否有dtUH和dt的换算逻辑没有就补上。最稳妥的做法是确认样例数据的时段长和单位线时段长一致如果数据文件里明确规定时段长就按实际值重新生成单位线不硬套默认。5.5 Matlab版本与文件读取老代码跑不起来的兼容坑现象代码逻辑完全没问题但一运行就报错提示xlsread未定义或中文表头读出来乱码更常见的是load一个mat文件后变量名对不上。 原因老代码是多年以前的Matlab版本写的那时xlsread还能用后来版本推荐readtable接口另外数据文件里保存的变量名和主函数里引用的变量名大小写不一致Matlab对变量名大小写敏感这种问题藏得很深。 解决先看报错行把xlsread整体替换成readmatrix或readtable表头处理统一用readtable的VariableNamingRule参数。mat文件里的变量名不匹配就在命令行输入whos -file 文件名查看实际存储的变量名再决定改数据还是改代码。中文路径也会导致读取失败把所有路径改成纯英文目录这一步能消掉一批莫名其妙的报错。排查这类问题有个通用顺序先确认数据能读进来再确认变量名能对上最后才看数值合理性。很多时间耗在调参上结果问题出在最基础的数据读取环节属于排查顺序错了。还有一点如果代码是用老版本脚本写的里面用了大量全局变量新版本运行前先查全局变量是否在每次运行时被正确重置这个坑会让模拟结果变得和上一次运行相关属于典型的玄学问题。6. 结果可信度验证三个土办法和一个率定习惯模型跑通不等于结果可信率定参数调完也不代表模型能用。我每次交付结果前必做三组检查都不需要额外写复杂代码。第一组是总量闭合检验。把汛期逐日累计实测径流和累计模拟径流做散点图看斜率在不在0.9到1.1之间。斜率偏大说明模拟产流过多偏小说明产流参数整体偏低。这个检验比NSE直观一眼就能看出水量守恒有没有被模型破坏。% 总量闭合检验累计径流对比 figure; plot(cumsum(Qobs), cumsum(Qsim), o); hold on; plot([0 max(cumsum(Qobs))], [0 max(cumsum(Qobs))], k--); xlabel(累计实测径流); ylabel(累计模拟径流);第二组是逐场洪水的峰现时间表。把每场洪水的峰现时间误差和洪峰相对误差列成表格而不是只算所有场次平均性能。平均NSE好看但掩盖单场失败一张逐场次的表能直接看出模型在哪种降雨型下失效。第三组是退水段对数化对比。把实测和模拟的退水段取对数画在同一个坐标里看斜率是否一致。斜率一致代表蓄泄关系对得上斜率差太多说明汇流参数有问题。% 退水段对比选取洪峰后的退水段取对数比较 idx find(Qobs max(Qobs)):find(Qobs 10, 1, last); plot(log(Qobs(idx)), k-); hold on; plot(log(Qsim(idx)), r--);最后一个习惯帮了我很多次每次率定只动两个参数改完记录下改了什么、结果怎么变。自动率定工具看起来省事但参数之间相互补偿调完经常得到一组无法解释的数值。手动记录能让你清楚知道每个参数在控制模型的哪一部分下次换流域时效率明显更高。希望帮到你。本文还有配套的精品资源点击获取
