水文圈里有一句老话模型的骨架是产流灵魂是分水源。最近我把三水源新安江模型用MATLAB完整实现了一遍顺带把降雨后水分在土壤里的湿润过程、张力水蓄量变化、自由水蓄水库的分配逻辑全部做了可视化。很多人看到“matlab液体湿润模拟”这个词会以为是个流体仿真其实放到水文模型语境下它指的就是降雨落地后水分如何下渗、蓄存、横向运动最终变成不同水源的径流。三水源新安江模型恰好是讲清楚这件事最经典的框架之一。这篇文章我会从模型设计思路讲起拆解三水源划分的原理再给出完整的MATLAB实现代码和可视化方案最后把我调试过程中踩过的坑、参数率定的经验一并整理出来。无论你是水文专业的学生、做洪水预报的工程师还是刚接触水文模拟的程序员照着这篇文章的思路都能把模型跑起来。1. 整体设计思路从湿润过程到三水源划分1.1 三水源新安江模型在算哪笔账新安江模型是蓄满产流模型的代表主要用在湿润和半湿润地区的降雨径流模拟。它把流域看成一个“大水箱”降雨下来以后先补充土壤缺水量等土壤蓄满后再产生径流。这里的“蓄满”不是指土壤物理意义上的完全饱和而是指包气带能够继续蓄存的水量达到上限。三水源新安江模型和原始版最大的区别在于它把产流后的水量进一步划分为地表径流、壤中流和地下径流三部分。为什么要这么分因为这三部分水在流域里的运动路径完全不同地表径流汇流速度最快几小时到几天就能到达出口断面壤中流在土壤里横向流动速度稍慢地下径流经过深层含水层的调蓄可以在降雨结束后持续补给河道很多天。如果只算一个总径流量就无法模拟出洪水的退水过程尤其是雨停之后流量缓慢衰减的那一段必须靠地下水和壤中流来支撑。所以在程序实现上三水源新安江模型的输入是逐日或逐小时降雨序列和蒸发能力序列输出是模拟流量序列。核心状态变量包括上层张力水蓄量、下层张力水蓄量、深层张力水蓄量和自由水蓄量整个模拟过程就是在不断更新这几个蓄量。1.2 “湿润模拟”到底模拟了什么把“matlab液体湿润模拟”和新安江模型放在一起理解其实很贴切。降雨之后水分先湿润表层土壤然后逐步下渗到下层这个过程中包气带的张力水蓄量不断增加这就是“湿润过程”。当土壤蓄水能力耗尽多余的水分进入自由水蓄水库再按照不同的出流路径被分成三股水这就是“再分配过程”。我用MATLAB做这件事时最直观的感受是湿润过程不是一条平滑曲线而是和降雨脉冲强相关的锯齿状上升。每次降雨后张力水蓄量跳升然后被蒸散发慢慢消耗下降这个动态过程在代码里就是一次循环里对WU、WL、WD三个变量的顺序操作。模型的“湿润”不仅仅是土壤变湿还包括了湿润锋的推进——上层蓄满后才开始补给下层下层蓄满后才动用深层水。所以这篇文章里的“湿润模拟”本质上是模拟流域包气带的蓄水状态随时间的变化它决定了产流的时机和产流量的大小。理解了这一点后面看代码就顺了。1.3 为什么选MATLAB而不是其他语言做水文模拟可选的语言很多Python有丰富的第三方库R也有水文包但我还是习惯用MATLAB原因有三。第一MATLAB的数组操作天然适合处理时间序列。降雨、蒸发、流量都是等间隔的一维序列在MATLAB里可以直接做向量化运算不需要显式写循环的地方就能省很多事。第二调试方便。水文模型最怕算到一半出现NaN或者负值蓄量MATLAB的编辑器可以随时打断点查看每个变量的实时值比改完一整段再运行要高效得多。第三画图顺手。模拟结果出来以后需要把实测流量和模拟流量画在一起对比把张力水蓄量和降雨画在一起看过程这些MATLAB都有现成的绘图函数几步就能出图。当然现在热词里经常看到“matlab下载”“matlab安装教程”这类搜索说明初学者入门的门槛主要在环境配置上。其实只要装好MATLAB本体水文模型用不到额外的工具箱核心代码几十行就能写完。2. 模型原理与核心参数拆解2.1 蓄满产流机制与张力水蓄水容量曲线蓄满产流的核心假设是流域内任意一点当包气带蓄水量达到田间持水量之后后续降雨全部产流。但流域不是一块均匀的板子土壤厚度、下垫面条件、地形坡度都不一样所以各点的蓄水容量不同。新安江模型用一条抛物线型的蓄水容量曲线来描述这种不均匀性。曲线横坐标是流域面积比例纵坐标是各点的蓄水容量。用一个参数B来控制曲线的凹凸程度B越大代表流域蓄水容量越不均匀B等于0时就是均匀流域所有点同时蓄满。WM是流域平均张力水容量通常取80到200毫米湿润地区取大值干旱地区取小值。产流计算的思路是已知当前流域平均蓄水量通过蓄水容量曲线反推出流域内已经蓄满的面积比例然后叠加降雨量凡是蓄水容量小于累计水量的区域都会产流。为了精确表达这个叠加过程标准做法是先求一个中间变量A再分两种情况计算产流量R。这个计算过程不复杂但公式推导容易绕晕代码实现里我建议封装成独立函数一次性测试通过后再放进主循环。2.2 三层蒸散发先上层、再下层、最后深层蒸散发是湿润过程里消耗水量的主要途径直接决定两次降雨之间土壤蓄量的下降速度。三水源新安江模型采用三层蒸散发模式把张力水分为上层、下层、深层三部分蒸散发按顺序消耗。上层张力水容量UM通常是15到20毫米分布在地表附近。降雨后上层先蓄满蒸散发也优先消耗上层的水。上层蓄量不足时蒸散发需求转向下层但下层的蒸发能力要打一个折扣这个折扣用系数C来表示C越大代表下层水分越容易被蒸发一般取0.1到0.2。如果下层也不够最后才动用深层蓄水。深层水一旦被动用就很难再补充回来所以深层蒸散发量通常很小。这个“先上层、后下层、再深层”的顺序在程序里就是几个if判断。我最初实现的时候把顺序写反了导致夏季模拟的土壤蓄量一直偏高后来检查蒸散发部分才发现问题。三层蒸散发的意义在于它让模型既能模拟出雨季表层土壤频繁干湿交替的规律又能让旱季深层水缓慢释放维持基流两件事同时兼顾。2.3 自由水蓄水库与三水源划分三水源划分是整个模型最精彩的部分。产流量先进入一个“自由水蓄水库”这个水库的特点是蓄满之前不出流蓄满之后按比例同时向壤中流和地下径流供水超过水库容量的部分就变成地表径流直接汇入河道。自由水蓄水库的容量SM是关键参数通常取10到30毫米。SM小的时候降雨稍微大一点就产生地表径流洪水过程线尖瘦SM大的时候更多水分被蓄在土壤里慢慢释放洪水过程线肥胖。壤中流出流系数KI和地下径流出流系数KG决定了自由水向两条路径的分配比例且KI加上KG必须小于1因为还有一部分水要留在水库里继续调节。具体计算时产流面积上自由水蓄量如果超过SM超出的部分直接作为地表径流RS然后当前蓄量分别乘以KI和KG得到壤中流RI和地下径流RG。这里特别要注意产流面积比例FR的换算自由水蓄水库只在产流面积上存在非产流面积上不发生水分交换所以计算时要先把流域平均的S换算成产流面积上的蓄量算完再折算回全流域尺度。2.4 汇流计算与关键参数速查表三股水源产生之后各自经过一个线性水库的调蓄再汇合到出口断面。线性水库的消退系数CS、CI、CG分别控制地表径流、壤中流和地下径流的衰减速度。地表径流的消退系数CS一般取0.7到0.9因为地表汇流快前一天的蓄量留存比例低壤中流CI取0.6到0.9地下径流CG则要取0.98到0.998地下水库调蓄能力强流量衰减非常慢所以前一天的水量几乎都留到了今天。为了方便查参数我整理了下面这个速查表也是我调试时的初始值参考参数含义常见取值范围初始参考值KC蒸散发折算系数0.8—1.21.0UM上层张力水容量(mm)15—2520LM下层张力水容量(mm)60—9080DM深层张力水容量(mm)40—8060B蓄水容量曲线指数0.1—0.40.3C深层蒸散发系数0.1—0.20.16SM自由水蓄水库容量(mm)10—3028EX自由水蓄水容量曲线指数1.0—1.51.0KI壤中流出流系数0.2—0.50.35KG地下径流出流系数0.2—0.50.40CS地表径流消退系数0.7—0.90.85CI壤中流消退系数0.6—0.90.85CG地下径流消退系数0.98—0.9980.99汇流计算里还有一个容易被忽略的点单位换算。模型算出来的径流深单位是毫米每天要变成流量单位立方米每秒必须乘以流域面积和换算系数。这个系数等于面积乘以1000再除以86400漏掉这一步的结果就是模拟流量偏差几个数量级。3. MATLAB代码从零实现3.1 数据准备与参数初始化写代码之前先把数据整理好。模型需要三个时间序列降雨量P单位mm、蒸发能力EM单位mm、实测流量Qobs用于对比单位m³/s。时间步长这里按天处理实际项目里如果是小时尺度参数取值要做相应调整特别是消退系数时间步长越小系数越接近1。参数我建议用结构体保存方便后续修改和调用。初始化代码如下% 参数设置 para.KC 1.0; para.UM 20; para.LM 80; para.DM 60; para.WM para.UM para.LM para.DM; % 总张力水容量 para.B 0.3; para.C 0.16; para.SM 28; para.EX 1.0; para.KI 0.35; para.KG 0.40; para.CS 0.85; para.CI 0.85; para.CG 0.99; % 状态变量初始化 WU 0; WL 0; WD 0; S 0; % 汇流蓄量初始化 QS 0; QI 0; QG 0; % 流域面积单位km2 Area 1000; % 单位换算因子mm/day - m3/s conv Area * 1000 / 86400; N length(P); % 模拟总时长 Qsim zeros(N, 1); WU_rec zeros(N, 1); WL_rec zeros(N, 1); WD_rec zeros(N, 1); S_rec zeros(N, 1);初始蓄量设置成0是一种简化处理。实际应用中如果模拟期前面有明显的退水段最好用前几天的反推法确定初始蓄量或者干脆把模拟期前面加一段预热期让模型自动调整到合理状态。3.2 主循环蒸散发与产流计算主循环是整个模型的心脏每天做四件事算蒸散发、算产流、分水源、汇流。蒸散发的计算顺序是优先消耗上层再下层最后深层for t 1:N % 蒸散发计算 EP para.KC * EM(t); if WU EP EU EP; EL 0; ED 0; else EU WU; D EP - EU; if WL para.C * D EL para.C * D; ED 0; else EL WL; ED para.C * D - EL; end end WU WU - EU; WL WL - EL; WD max(0, WD - ED);这里deep部分的处理我用了一个max(0, ...)来防止深层蓄量出现负值这在干旱条件下可能出现。接着做产流计算把当前的张力水总蓄量带进蓄水容量曲线公式% 产流计算调用子函数 W0 WU WL WD; R calc_runoff(P(t), W0, para);calc_runoff函数内部要解蓄水容量曲线的A值完整代码如下function R calc_runoff(PE, W0, para) WM para.WM; B para.B; WMM WM * (1 B); if PE 0 R 0; return; end if W0 0 A 0; elseif W0 WM A WMM; else % 数值求解 A - WM*(1-(1-A/WMM)^(B1)) W0 fun (A) A - WM * (1 - (1 - A/WMM)^(B 1)) - W0; A fzero(fun, [0, WMM]); end if A PE WMM R PE - (WM - W0); else R WM * ((1 - A/WMM)^(B 1) - (1 - (A PE)/WMM)^(B 1)); end R max(R, 0); end这段代码最关键的地方是用fzero数值求解A值。如果你懒得解这个方程也可以牺牲一点精度把蓄水容量曲线近似成均匀分布令B0这样A就等于W0产流公式退化成“蓄满前不产流、蓄满后全部产流”的简单形式。但对于水文模拟来说B参数对洪峰形状的影响很明显还是保留为好。3.3 分水源与汇流实现产流量R产生之后进入自由水蓄水库进行三水源划分。这里要计算产流面积比例FR可以用前面求出的A值间接得到——蓄满面积比例等于(1 - A/WMM)的B次方。由于calc_runoff函数内部才有A值我这里采用重新计算FR的简化方式% 分水源计算 WM para.WM; WMM WM * (1 para.B); if W0 0 FR 1; elseif W0 WM FR 0; else fun (A) A - WM * (1 - (1 - A/WMM)^(para.B 1)) - W0; A fzero(fun, [0, WMM]); FR (1 - A / WMM)^para.B; end % 产流面积上的自由水蓄量 S_area S / max(FR, 0.01) R / max(FR, 0.01); if S_area para.SM RS (S_area - para.SM) * FR; S_area para.SM; else RS 0; end RI para.KI * S_area * FR; RG para.KG * S_area * FR; S_new S_area * FR - RI - RG; S max(S_new, 0);这段代码比教科书上的写法更直观但要注意RS、RI、RG三者的单位都是毫米每天代表全流域平均的径流深。实测中合流后要验证一下水量平衡R应该等于RS加RI加RG再加上S的变化量如果不满足说明中间有蓄量计算的逻辑错误。汇流部分就比较简单了三个线性水库各自消退% 汇流计算 QS para.CS * QS (1 - para.CS) * RS * conv; QI para.CI * QI (1 - para.CI) * RI * conv; QG para.CG * QG (1 - para.CG) * RG * conv; Qsim(t) QS QI QG;汇流的物理含义是当天产生的水量不会全部当天流到出口而是部分留在河道或水库中第二天继续流出。消退系数就是描述这个“留存比例”的。3.4 结果可视化和精度评估模拟完以后最重要的就是看图。一张标准的对比图上面画实测和模拟流量过程线下面画降雨柱状图倒置这是水文模型验证的通用画法。MATLAB里直接用tiledlayout可以实现figure; tiledlayout(2, 1); % 上半部分流量对比 nexttile; plot(t, Qobs, k-, LineWidth, 1.2); hold on; plot(t, Qsim, r--, LineWidth, 1.2); legend(实测流量, 模拟流量); ylabel(流量 (m3/s)); title(三水源新安江模型模拟结果); % 下半部分降雨倒置 nexttile; bar(t, P, FaceColor, [0.3 0.6 0.9]); set(gca, YDir, reverse); ylabel(降雨 (mm)); xlabel(天数);精度评估常用纳什效率系数NSE和相对误差BIAS。NSE计算公式是NSE 1 - sum((Qobs - Qsim).^2) / sum((Qobs - mean(Qobs)).^2);NSE大于0.7就说明模型模拟效果可以接受大于0.85算优秀。但要注意NSE对大流量敏感如果洪峰对得很准但退水段偏大NSE也会虚高。我一般会同时看模拟和实测的总水量相对误差控制在正负10%以内才算合格。4. 常见问题与调参避坑实录4.1 水量平衡对不上先查单位换算我做这套模型时第一次跑出来的模拟流量比实测大了100多倍排查了半天才发现是汇流时单位换算错了。径流深是毫米每天单位面积换算因子是1000/86400再乘以面积这个因子算出来大约是每秒、每平方毫米多少立方米基础单位没理清就会出问题。建议每跑完一个时段就手动核对一次水量平衡累计降雨量减去累计蒸散发量应该等于累计产流量加土壤蓄量变化量。把这一步写成代码自动检查可以省掉大量无谓的调试时间。还有一个坑是面积单位和时间步长不匹配。如果用小时步长换算因子分母就要从86400变成3600消退系数也要整体调大。很多初学者拿日模型的参数直接跑小时模型流量过程线震荡得像锯齿。4.2 洪峰流量偏高或偏低调整哪些参数模拟洪峰偏高的常见原因是SM设得太小。自由水蓄水库容量小降雨很快蓄满并产生地表径流洪峰自然就高。反过来SM调大更多水分被土壤蓄住洪峰就变矮变胖。这里的SM相当于一个缓冲器容量越大对洪峰的削减作用越强。如果洪峰时间对不上优先检查CS。地表径流消退系数CS控制着洪峰汇流速度CS越大峰值出现越晚、过程越平缓。还有一种情况是洪峰形状对但退水段掉得太快这多半是CG设得太小地下径流维持不住后期的基流。我调试时习惯从大到小逐步消减峰值一次只调一个参数看过程线变化趋势不追求一步到位。4.3 湿润过程模拟失真问题出在蒸散发张力水蓄量的变化轨迹反映了流域湿润状态。我发现模拟的土壤蓄量在雨后恢复得太快或者太慢问题通常出在KC参数上。KC是蒸发能力的折算系数把蒸发皿测得的EM换算成实际蒸散发能力。KC调大会让湿润后的土壤快速变干蓄量下降斜率变陡KC调小则相反。还有一种情况是下层和深层蓄量常年不增加说明LM和DM设置偏大水分始终蓄在上层。我记得有一次模拟半湿润区流域基流一直偏低后来把LM从90调到70地下径流立刻涨了不少。原因是下层蓄水容量太大水分被截留在中层到不了深层也就无法形成稳定的地下径流。这提醒我张力水容量的分配直接影响水分垂直运动的路径绝对不能只看WM总量。4.4 程序调试与参数率定的个人经验写MATLAB代码时我强烈建议把每个中间状态变量都保存下来。WU_rec、WL_rec、WD_rec、S_rec这些数组不仅用于画图也是排查逻辑错误的关键。比如产流计算出错时先看R的序列是否和降雨对应再看S是否有异常累积。参数率定方面人为手工调参效率太低。我在实测中通常先用遗传算法自动率定一轮再利用水文经验手动微调。MATLAB Optimization Toolbox里提供了ga函数目标函数可以用NSE或KGE。但自动率定有个问题容易出现过拟合模拟期效果好验证期一塌糊涂。我的做法是把序列分成率定期和验证期率定结束后必须用另一段时间检查两段效果都好才算通过。下面把常见问题整理成速查表现象优先排查参数调整方向洪峰整体偏高SM、CS增大SM或CS洪峰整体偏低SM、KG减小SM检查KG是否过大退水段掉得太快CG增大CG接近0.99基流整体偏小LM、DM、KG减小LM/DM增大KG土壤蓄量恢复过慢KC增大KC提高蒸发消耗雨后湿润期过长C、KC增大C或KC流量过程线锯齿单位换算、步长检查时间步长与换算因子4.5 一个常被忽略的细节预热期与初始蓄量模型初始蓄量如果直接设为零前几天的模拟误差会非常大尤其是流域初始土壤较湿润时模型需要一段时间“填平”蓄水缺口。解决方法是把模拟期往前扩展半年到一年等模型进入稳定状态后再截取结果。我踩过一次很深的坑有一个流域前期有连续降雨我却把WU、WL、WD初始值设成0导致模拟初期产量明显偏小无论怎么调参数都改善不了。后来在时间序列前面加了三个月的预热期问题立刻解决。做水文模拟的都懂一句话宁可把预热期加长也不要让初始状态干扰参数率定的判断。跑完这版模型我最大的体会是三水源新安江模型看起来只有十几个参数但每个参数背后都有明确的物理含义调试过程其实是在不断加深对流域水文过程的理解。MATLAB实现本身不难难点在于怎样用代码准确表达“上层先湿润、下层再补给、蓄满才产流、产流再分家”这个逻辑链条。如果你也是刚开始做水文模型别急着追求代码高级先把水量平衡算清楚把湿润过程画出来你会发现很多模型问题都能从图上直接看出来。