水文圈里有一句老话：模型的骨架是产流，灵魂是分水源。最近我把三水源新安江模型用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，地下水库调蓄能力强，流量衰减非常慢，所以前一天的水量几乎都留到了今天。

为了方便查参数，我整理了下面这个速查表，也是我调试时的初始值参考：

汇流计算里还有一个容易被忽略的点：单位换算。模型算出来的径流深单位是毫米每天，要变成流量单位立方米每秒，必须乘以流域面积和换算系数。这个系数等于面积乘以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)^(B+1)) = 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值。如果你懒得解这个方程，也可以牺牲一点精度，把蓄水容量曲线近似成均匀分布，令B=0，这样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。但自动率定有个问题：容易出现过拟合，模拟期效果好，验证期一塌糊涂。我的做法是把序列分成率定期和验证期，率定结束后必须用另一段时间检查，两段效果都好才算通过。

下面把常见问题整理成速查表：

4.5 一个常被忽略的细节：预热期与初始蓄量

模型初始蓄量如果直接设为零，前几天的模拟误差会非常大，尤其是流域初始土壤较湿润时，模型需要一段时间“填平”蓄水缺口。解决方法是把模拟期往前扩展半年到一年，等模型进入稳定状态后再截取结果。

我踩过一次很深的坑：有一个流域前期有连续降雨，我却把WU、WL、WD初始值设成0，导致模拟初期产量明显偏小，无论怎么调参数都改善不了。后来在时间序列前面加了三个月的预热期，问题立刻解决。做水文模拟的都懂一句话：宁可把预热期加长，也不要让初始状态干扰参数率定的判断。

跑完这版模型，我最大的体会是：三水源新安江模型看起来只有十几个参数，但每个参数背后都有明确的物理含义，调试过程其实是在不断加深对流域水文过程的理解。MATLAB实现本身不难，难点在于怎样用代码准确表达“上层先湿润、下层再补给、蓄满才产流、产流再分家”这个逻辑链条。如果你也是刚开始做水文模型，别急着追求代码高级，先把水量平衡算清楚，把湿润过程画出来，你会发现很多模型问题都能从图上直接看出来。