每次做水文预报最绕不开的就是新安江模型。我最早接触这个概念是在做中小流域洪水预报的时候那时手里只有日降雨量和蒸发资料要估计一场洪水的产流量翻了很多教材最后发现真正能落地、能写进程序、参数含义又清晰的还是三水源新安江模型。这个模型用matlab来实现代码量不大但每一步的水量平衡关系都清清楚楚特别适合入门水文模型的人当作第一个练手项目。这篇文章就是我在matlab里实现三水源新安江模型的全过程记录。我会把模型怎么拆、每个模块的公式怎么转成代码、参数怎么定、结果怎么调一条线讲清楚。不管你是做洪水预报、水资源评价还是单纯想搞懂蓄满产流和三种水源是怎么被分开的这篇都能给你一个可以直接跑起来的起点。1. 模型核心逻辑为什么湿润地区产流要先算蓄满三水源新安江模型最本质的一个假设是蓄满产流。这和超渗产流是两条完全不同的路子。超渗产流看的是雨强和地面下渗能力谁大而蓄满产流看的是包气带土壤含水量有没有到田间持水量。在湿润地区比如南方大部分流域降雨频繁且总量大包气带容易蓄满所以蓄满产流的假设更贴合实际。那“蓄满”意味着什么意味着土壤不再能继续存住水分后续降雨直接在地表形成径流。新安江模型的巧妙之处在于它把整个产流过程拆成了两步先算整个流域有多少面积蓄满了蓄满的面积上产生径流再把这部分径流用一个自由水蓄水库去分水源——分成地面径流、壤中流和地下径流。这就是“三水源”的来历。用matlab实现的时候整场降雨要按小时或天来一步步推。每一时段里降雨先进入张力水蓄水层蓄满之后的水进自由水蓄水层自由水层也满了之后就产流。这个逻辑顺序如果理不清楚后面代码很容易写得混成一团。我的建议是先把下面这张水量平衡表在草稿纸上画出来再动手写代码蓄水层输入输出平衡条件张力水层降雨扣除雨期蒸发下渗补给自由水层不得超蓄容量WM自由水层张力水下渗地面径流、壤中流、地下径流不得超蓄容量SM河网三种水源汇集出口断面流量汇流滞后这里有个新手容易忽略的点张力水和自由水是两种不同“性质”的水。张力水被土壤颗粒吸附很难自由流动只有当超过田间持水量后才会变成自由水。所以在程序里这两个蓄水层要分开设置容量、分开计算不能图省事合并成一个变量。整场降雨过程中每个时段的蒸散发损失优先从张力水中扣除。这个顺序也重要后面对照实测流量过程线时如果模拟的洪峰偏大往往就是蒸散发扣得太少土壤一直处在蓄满状态。2. 蒸散发计算模块三层蒸散发模型的代码实现蒸散发是产流计算里最容易“差不多”但又影响很大的部分。新安江模型的蒸散发按土层分为上层、下层和深层。三层之间不是平均分摊而是有一个“从上往下、不够再下一层”的转移逻辑。用代码写这个逻辑其实就是判断嵌套。先定义土壤含水量变量上层张力水蓄量WU下层张力水蓄量WL深层张力水蓄量WD对应容量为UM、LM、DM总的张力水容量WM UM LM DM。蒸散发的目标是三层的蓄水都要被消耗但只有上层还有水时优先蒸发上层上层水不够了才轮到下层和深层。这一段是matlab函数我通常单独写成evap_transfer.m方便在产流循环里反复调用function [WU, WL, WD, E] evap_transfer(P, EP, WU, WL, WD, UM, LM, DM) % 输入: P 时段降雨, EP 时段蒸发能力, WU/WL/WD 三层张力水蓄量 % 输出: 更新后的WU/WL/WD以及实际蒸发量E E 0.0; % 先蒸上层 if WU P EP E EP; WU WU P - EP; % 上层水量足够下层和深层不蒸发 return; else E WU P; % 上层全部蒸掉 WU 0.0; rest EP - E; % 剩余蒸发能力 end % 下层蒸发 if WL rest WL WL - rest; E E rest; rest 0.0; else E E WL; rest rest - WL; WL 0.0; end % 深层蒸发 if rest 0 if WD rest WD WD - rest; E E rest; else E E WD; WD 0.0; end end % 此时如果还有剩余蒸散发能力直接放弃因为已经没有水可蒸 end这里有一个细节降雨P先补进上层再被蒸发。也就是说这个时段下的雨优先留在上层供蒸发消耗而不是先渗到下层。这也是实际物理过程的近似降雨先湿润表层。我记得第一次调试这个函数的时候忘记在“上层水量足够”的分支里直接return结果下层和深层也被扣了水模拟出来的蒸散发总量偏大土壤蓄量一直偏低洪峰模拟得又矮又胖。原因是把“返回”和“继续”的逻辑写错了。所以在写水位转移类代码时每个分支的出口都要检查一遍别让水流到了不该去的地方。深层蒸散发的存在还有另一层物理意义在久旱季节上层和下层都干了但深层土壤中仍然存着水分植被根系可以直接从深层吸水蒸腾。如果模型里没有这一层旱季的基流会模拟得偏大因为水分没有被植被消耗掉全部变成了枯季径流。3. 产流计算模块蓄水容量曲线和自由水蓄水库蒸散发算完接下来是核心的产流计算。三水源新安江模型之所以能模拟出地面径流、壤中流和地下径流三种水源关键在于两个蓄水容量曲线一个是张力水的蓄水容量曲线用它来算产流面积另一个是自由水蓄水库用它来分配三种水源的比例。3.1 张力水蓄水容量曲线全流域各点的张力水蓄水容量并不一致有的地方薄、有的地方厚。新安江模型用一条抛物线来概化这个不均匀性。蓄水容量曲线的表达式是f/F 1 - (1 - Wm/WMM)^B其中WMM是流域最大点蓄水容量B是蓄水容量曲线指数。B越大表示流域蓄水容量越不均匀。有了这条曲线再结合当前时段的流域平均张力水蓄量W就能算出当前降雨下有多少面积产生了径流。在matlab里计算蓄满产流我用的是赵人俊先生教材里的经典算法。设时段初张力水蓄量为W0降雨P蒸散发已经在前面扣除了那么产流量R计算公式涉及一个中间变量AA WMM * (1 - (1 - W0 / WM)^(1 / (1 B)))如果 P A WMM R P - (WM - W0) WM * (1 - (A P) / WMM)^(1 B)这个公式看着复杂推导起来其实就是对蓄水容量曲线做空间积分。我不建议死记公式要在代码里把物理含义注释清楚。以下是产流模块函数其中W0表示时段初张力水总蓄量function [R, W1] runoff_generate(P, W0, WM, B, IMP) % 不透水面积占比IMP直接产流 % 张力水蓄水容量曲线参数 WMM WM * (1 B); % 最大点蓄水容量 A WMM * (1 - (1 - W0 / WM)^(1 / (1 B))); if A P WMM R P - (WM - W0) WM * (1 - (A P) / WMM)^(1 B); else R P - (WM - W0); end % 扣除不透水面积的影响 R R * (1 - IMP) P * IMP; W1 W0 P - R; % 时段末张力水蓄量 end这个函数跑通之后你会得到一个很直观的现象当土壤含水量W0越高同样一场雨的产流量R就越大。这就解释了同一场暴雨为什么前期湿润的条件下容易发洪水前期干旱的时候降雨都被土壤吸收了。3.2 自由水蓄水库分水源产流量R不会全部成为地面径流。模型假设产流面积上存在一个自由水蓄水库这个水库的容量是SM。它接受来自张力水的下渗补给同时侧向排出壤中流向下渗漏补给地下径流溢流部分才是地面径流。自由水蓄水库的平衡方程S2 S1 R - RS - RI - RGRS是地面径流RI是壤中流RG是地下径流。地面径流的出流系数KG和壤中流出流系数KI是模型参数。这三个水源的分配比例与自由水蓄量的关系写成代码function [RS, RI, RG, S2] free_water_reservoir(R, S1, SM, KG, KI, EX) % R: 进入自由水水库的产流量 % S1: 时段初自由水蓄量 % KG: 地下径流出流系数 % KI: 壤中流出流系数 % EX: 自由水蓄水容量曲线指数 if SM 0 RS R; RI 0; RG 0; S2 0; return; end SMM SM * (1 EX); FR 1.0; % 产流面积比例单点模型简化为1 % 根据自由水蓄水容量曲线积分 AU SMM * (1 - (1 - S1 / SM)^(1 / (1 EX))); if R AU SMM RS (R - SM S1 SM * (1 - (R AU) / SMM)^(1 EX)) * FR; else RS (R - SM S1) * FR; end S2 S1 R - RS; RI S2 * KI * FR; RG S2 * KG * FR; S2 S2 - RI - RG; end这段代码里有个经验参数EX自由水蓄水容量曲线指数。它和张力水的B参数类似都用来表示蓄水容量的空间不均匀性。但两者的取值逻辑不同B一般取0.2到0.4EX常见取值在1.0到1.5之间。如果你刚开始试算EX取1.0就行后面再根据模拟效果微调。分完三种水源后各自还要经过汇流才能到达流域出口。这就涉及下一个模块。4. 汇流计算模块从坡面到河网的三个滞后环节三水源各自汇流的路径不一样。地面径流走得最快壤中流是横向的浅层流动地下径流则是漫长的基流过程。新安江模型用“线性水库”来模拟这三种汇流每一种水源都对应一个蓄泄关系。先说单位线或纳什瞬时单位线。在实际编程中许多教材版本用滞后演算法lag-and-route或线性水库串联来模拟地表汇流。我在这里用的是最常用的纳什单位线法把地面径流过程先转换为流域出口断面的流量过程再叠加壤中流和地下径流。线性水库的通用函数写成这样function Q_out linear_reservoir(Q_in, CS, Q_prev) % 线性水库Q_out CS * Q_out_prev (1 - CS) * Q_in Q_out CS * Q_prev (1 - CS) * Q_in; end其中CS是消退系数反映水库对洪水的调蓄削峰作用。CS越接近1水库的调蓄能力越强出流越平缓。这个函数虽然只有一行但它是汇流计算里的基本积木。地面径流用CS0.8到0.9壤中流用CS0.9到0.98地下径流用CS0.98到0.998。从数值上就能看出来地下基流被调蓄得非常慢一场暴雨后很久地下径流还在缓慢释放。如果把整个流域当作一个整体不考虑空间分布可以用一个“三水源汇流”的完整函数把三种水源分别过一遍线性水库最终叠加得到出口断面总流量function [Q_total, QS_prev, QI_prev, QG_prev] route_three_sources(RS, RI, RG, ... CS, CI, CG, QS_prev, QI_prev, QG_prev) % 地面径流汇流 QS CS * QS_prev (1 - CS) * RS; % 壤中流汇流 QI CI * QI_prev (1 - CI) * RI; % 地下径流汇流 QG CG * QG_prev (1 - CG) * RG; % 出口总流量 Q_total QS QI QG; end这里有一个汇流计算中特别常见的坑你模拟出来的洪峰形状总是“太尖”或者“太胖”。太尖说明CS取值太大洪水被调蓄得不够太胖说明CS太小把洪水过程拉得太长。调参的时候先看洪峰出现的时间是否吻合再看退水段是否贴合实测。如果洪峰提前说明地面径流的汇流时间被低估了需要增大CS如果洪峰滞后就需要减小CS。在实际操作中我一般把单位时段取为1小时。如果你的资料是日尺度那么地面径流的CS基本可以取0.98以上因为一天的时段内地面径流基本已经流完了。用日资料来率定CS时参数敏感性会变差这一点在做日模型时要心里有数别把日尺度的率定参数拿去用于小时尺度的预报。5. 参数表与敏感性分析哪些参数要用实测率定哪些可以直接取经验值三水源新安江模型的参数可以分成四类每一类的获取难度和敏感性都不一样。写代码之前先把参数表定下来后面调参会省很多事。我常用的参数整理成下表参数物理含义常见取值范围获取方式KC蒸散发折算系数0.8 ~ 1.2率定UM上层张力水容量10 ~ 30 mm经验LM下层张力水容量60 ~ 100 mm经验DM深层张力水容量10 ~ 40 mm经验WM总张力水容量80 ~ 200 mm经验率定B张力水蓄水容量曲线指数0.1 ~ 0.5率定IMP不透水面积比0 ~ 0.05地理信息SM自由水蓄水容量5 ~ 50 mm率定EX自由水蓄水容量曲线指数1.0 ~ 1.5经验KG地下径流出流系数0.1 ~ 0.5率定KI壤中流出流系数0.1 ~ 0.5率定CS地面径流消退系数0.5 ~ 0.95率定CI壤中流消退系数0.5 ~ 0.95率定CG地下径流消退系数0.95 ~ 0.999率定L河网滞后时间0 ~ 3 h地理信息这里单独说一下KC。它是蒸发皿实测蒸发与流域实际蒸散发能力的比值。一般取0.8左右干旱地区可以到0.6湿润地区可以到1.0以上。这是影响水量平衡的第一敏感参数它偏大流域蓄量就偏低它偏小土壤一直偏湿后一场小雨也能产生很大的径流。所以率定参数时我习惯从KC开始调先把多年水量平衡对上再调洪峰形状相关参数。SM也是一个很关键的自由水参数。SM越大意味着土壤能存住更多的自由水地面径流占比就小洪水过程就越平缓。SM越小地面径流占比大洪峰尖瘦。实际率定时你会发现SM和KG、KI之间有明显相关性存在异参同效的问题——即不同组合的参数都能得到相似的拟合结果。对待这个问题我的原则是不要盲目相信自动优化得到的最优参数要结合下垫面情况约束参数范围让参数在物理上说得通。6. 主程序框架如何把各个模块串起来跑一场洪水模块都准备好之后主程序就是按时间步长逐时段推进。整个过程可以理解成一个“加工流水线”降雨和蒸发进来先过蒸散发层再过张力水产流层再过自由水水库分水源最后三路汇流叠加得到出口流量。这段主循环代码是整个模型的主干% 主程序示例newanjiang_simulate.m % 读取降雨序列P、蒸发序列EP以及初始蓄量状态 n length(P); % 时段数 Q_total zeros(n, 1); % 出口总流量序列 QS zeros(n, 1); QI zeros(n, 1); QG zeros(n, 1); % 初始状态正常情况下湿润地区模型预热期至少3个月 WU0 20; WL0 50; WD0 10; S_free 5; % 自由水蓄量初值 QS_prev 0; QI_prev 0; QG_prev 0; % 模型参数 KC 0.9; UM 20; LM 80; DM 20; WM UM LM DM; B 0.3; IMP 0.01; SM 15; EX 1.0; KG 0.3; KI 0.3; CS 0.85; CI 0.92; CG 0.99; WU WU0; WL WL0; WD WD0; for t 1:n % 蒸散发计算 [WU, WL, WD, E_actual] evap_transfer(... P(t), EP(t) * KC, WU, WL, WD, UM, LM, DM); % 张力水总蓄量 W0 WU WL WD; % 产流计算注意这里产流前的降雨P(t)已经扣除蒸散发消耗部分 P_net P(t) - E_actual; % 进入土壤的水量包含土壤蓄水增量与产流量 [R, W1] runoff_generate(P_net, W0, WM, B, IMP); % 更新张力水蓄量这里按三层各自恢复并不唯一简化处理为按比例分配 ratio W1 / max(W0, 1e-6); WU WU * ratio; WL WL * ratio; WD WD * ratio; % 注意三层按比例分配是一种简化严格的做法是分层计算但工程上比例法可用 % 自由水蓄水库分水源 [RS, RI, RG, S_free] free_water_reservoir(R, S_free, SM, KG, KI, EX); % 三水源汇流 [Q_total(t), QS_prev, QI_prev, QG_prev] ... route_three_sources(RS, RI, RG, CS, CI, CG, QS_prev, QI_prev, QG_prev); QS(t) QS_prev; QI(t) QI_prev; QG(t) QG_prev; end % 输出结果 plot(Q_total, b-, LineWidth, 1.5); hold on; plot(QS, r--); plot(QI, g-.); plot(QG, k:); legend(总流量, 地面径流, 壤中流, 地下径流);上面的主程序里有一个简化每时段计算完产流后新算出的W1按同一个比例更新三层张力水蓄量。这是不严格的。严格做法应该把三层分别做水量平衡上层优先蓄满再进入下层。但工程上很多简化版新安江模型代码就是这么写的对总产流的影响不显著因为产流总量取决于WUWLWD的总和而比例更新法保持了总和正确。运行这段程序之后你会得到一组流量过程线。如果数据来自实际流域就把模型输出和实测流量画在一起看拟合效果。我常用的评价指标是两个确定性系数NSE和水量平衡误差RE。NSE大于0.7算基本可用大于0.85就是很好的模拟效果了。RE控制在5%以内说明水量平衡没跑偏。NSE 1 - sum((Q_sim - Q_obs).^2) / sum((Q_obs - mean(Q_obs)).^2); RE sum(Q_sim) / sum(Q_obs) - 1;7. 参数率定的实际操作手动率定还是自动率定参数率定是水文模型里最考验经验的一环。三水源新安江模型有十几个参数理论上每个都要率定。但实际操作中我不建议一上来就自动化率定。先用“手动经验”把模型调到大致合理再考虑用优化算法微调这个顺序能让你的思路始终保持清醒。手动率定有个口诀叫“先水量、后过程先大后小、逐步逼近”。具体来说第一步调KC和WM。让模拟的总径流量和实测的总径流量大致相等。水量都不平衡后面调什么都是白搭。第二步调SM和B。这两个参数主要影响地面径流和地下径流的比例分配进而决定洪峰的陡缓和基流大小。调的时候看两个特征洪峰的高度和退水段的形状。第三步调CS、CI、CG。这三个消退系数控制汇流速度直接影响洪峰出现的时间和退水段的坡度。第四步微调KI和KG。这两个参数对总径流量的影响比较小但对基流的比例影响明显。KI大一些洪峰后的一段“腰部”流量就会更突出KG大一些基流更充沛。手动率定到NSE大于0.6之后可以接上自动率定。matlab的优化工具箱里fmincon、ga都可以用。目标函数一般取NSE最大化或NSE和水量平衡误差的加权组合。不过自动率定有个陷阱参数容易跑到物理上不合理的范围。比如SM调到80mm这在湿润地区几乎不可能。所以自动率定时一定要给参数加边界约束。我用过粒子群算法来率定新安江模型初始种群设置为500迭代次数200效果比fmincon好不容易陷入局部最优。但粒子群也有随机性每次运行的结果会略有不同。稳妥的做法是多次运行取NSE最高且参数最合理的一组。8. matlab实现中的典型坑湿润初始化和时间尺度一致性问题最后说几个我在matlab实现三水源新安江模型时踩过的坑。这些坑在教科书上不会写但对调试效率影响非常大。第一个坑是模型预热期不足。新安江模型的初始土壤含水量对前几场洪水的模拟影响极大。如果土壤初始含水量设置得偏高第一场小雨也会算出很大的洪峰如果设置偏低第一场大洪水的洪峰又会被严重低估。解决办法是用模拟年份前至少半年或一年的日资料做预热spin-up让模型的土壤含水量自动调整到一个合理的平衡状态。正式统计模拟效果的时段从预热期之后开始。我第一次做的时候嫌麻烦直接用估算的初始蓄量开跑结果前两场洪水的NSE惨不忍睹白白浪费了调参的时间。第二个坑是时间单位不一致。降雨如果是mm/h蒸发如果是mm/d那产流计算就全错了。我曾经因为蒸发资料是日值、降雨是小时值直接把两者放在同一个循环里计算水量平衡误差高达30%一开始还以为是参数问题后来才发现是单位没统一。解决方法是进入模型之前统一转换为mm/时段并且保证降雨量级和蒸散发量级匹配。第三个坑是free_water_reservoir函数里的产流面积比例FR。在主程序示例中我把FR简化为1这是基于全流域蓄满的假设。但在半湿润地区或者大流域产流面积并不等于全部面积。严格的三水源新安江模型里FR应等于张力水蓄满的面积比例也就是产流面积比。这一步如果省略模型在干旱时期的模拟误差会明显偏大。我的建议是刚开始练手时可以用FR1的简化版但真正做研究时要把产流面积比例的计算加回去。第四个坑是matlab里for循环的效率和数值稳定性。如果你模拟的时间序列很长比如几十年日尺度数据for循环里有大量重复计算。可以在循环之前把所有常数参数预先算好比如WMM WM * (1B)、SMM SM * (1EX)避免每次迭代都重复计算。另外代码里所有除法都要加一个小量保护防止除零。比如计算ratio的时候我用max(W0, 1e-6)来避免分母为零。就个人调试经验而言三水源新安江模型是水文模型里“性价比”最高的一个。它的物理概念清晰、代码实现简洁、参数不是特别多但模拟效果在湿润地区普遍不错。用matlab把它写一遍你对蓄满产流、水源划分、汇流过程的理解会从“看书懂了”变成“真正会用了”。而且这套代码框架后面改造成分布式新安江模型、或者耦合数据同化模块都是现成的底子。