简介：一个基于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:t+m-1) = Qs(t:t+m-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) |
|---|---|---|---|---|
| 1 | 10 | 10 | 0 | 0 |
| 2 | 10 | 20 | 9 | 1.8 |
| 3 | 10 | 30 | 9 | 6.3 |
| 4 | 10 | 40 | 9 | 7.2 |
| 5 | 10 | 50 | 9 | 6.3 |
| 6 | 10 | 60 | 9 | 5.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--');最后一个习惯帮了我很多次：每次率定只动两个参数，改完记录下改了什么、结果怎么变。自动率定工具看起来省事，但参数之间相互补偿，调完经常得到一组无法解释的数值。手动记录能让你清楚知道每个参数在控制模型的哪一部分，下次换流域时效率明显更高。希望帮到你。
本文还有配套的精品资源，点击获取
