用MATLAB把三水源新安江模型跑起来从湿润过程模拟到流量预报但凡在湿润地区做过洪水预报或径流模拟基本绕不开新安江模型。这个诞生于上世纪七十年代的降雨径流模型到今天依然是国内水文预报的主流工具之一尤其在三水源划分和蓄满产流机制上它把湿润地区的产流过程刻画得相当通透。最近不少同行问我能不能用MATLAB实现一套完整的三水源新安江模型顺带把土壤湿润过程也模拟出来——也就是降雨怎么一步步把干土浸润、蓄满、产流、分水源的全过程。这篇文章我就以三水源新安江模型为核心手把手把模型原理拆开再把MATLAB实现的关键步骤、参数设置、水量平衡检查、常见坑全部过一遍。无论你是正在做课程设计的学生还是刚接到洪水预报任务的工程师参考这套代码和思路都能快速搭起一个可用的小模型。1. 三水源新安江模型的原理与模拟目标1.1 为什么湿润地区绕不开蓄满产流说起产流机制旱地和非旱地的思路完全不同。北方干旱半干旱地区雨强大、下渗能力有限超渗产流是主导但到了湿润地区降雨频繁、土壤含水量长期偏高绝大多数径流来自“蓄满”之后的产流——也就是土壤先喝饱水多余的水再往外出。三水源新安江模型的核心假设就是蓄满产流每个点都存在一个张力水蓄水容量当累计降雨扣除蒸发达到这个容量后该点开始产流。这个概念很关键它意味着我们模拟的并不是一场雨的直接下渗过程而是整个流域土壤水库的蓄泄过程。用MATLAB做湿润模拟时核心就是跟踪每块网格或每个单元的土壤含水量变化看它什么时候从“欠饱和”变成“蓄满”。这也就解释了为什么标题里会同时出现“液体湿润模拟”和“新安江模型”——模型模拟的湿润过程本质上是土壤水分动态平衡而三水源新安江模型恰好把湿润地区的蓄满产流、水量分配和汇流过程串成了一条完整链路。1.2 三水源是怎么划分出来的蓄满产流之后产流的水并不全走地面。模型把总径流分成三部分地面径流RS降雨强度大、土壤表层蓄满后直接形成汇流最快。壤中流RSS下渗到土壤中层的自由水沿水平方向侧向流动速度中等。地下径流RG继续下渗到深层补给地下水后再缓慢汇入河网速度最慢。三水源的划分不是简单按比例切而是通过引入“自由水蓄水库”来实现。你可以把自由水想象成土壤中那些不吸附在颗粒表面、可以在孔隙间自由流动的水。当张力水蓄满后其余水分进入自由水蓄水库再按出流系数分为壤中流和地下径流只有当自由水蓄量超过其容量时多余部分才形成地面径流。这个结构的好处是物理概念清晰、参数不多而且把产流面积变化对径流的影响体现出来了——产流面积自由水蓄满的面积比例会随土壤湿润程度动态变化这是新安江模型最有价值的地方之一。1.3 模型整体结构蒸散发、产流、分水源、汇流三水源新安江模型的标准结构可以概括为四个计算环节蒸散发计算采用三层蒸散发模型按上层、下层、深层张力水蓄量逐层扣减。产流计算用蓄水容量曲线描述流域内张力水蓄水容量的空间不均匀性计算产流量。分水源计算用自由水蓄水库把总产流划分为地面、壤中、地下三部分。汇流计算地面径流用单位线汇流壤中流和地下径流用线性水库消退。整个流程在MATLAB里实现其实就是把一个物理过程拆成若干个状态变量和递推公式逐时段更新。只要把数据结构设计好整个程序可以控制在几百行以内而且很容易扩展成分布式版本。2. MATLAB实现前的准备数据、参数与代码架构2.1 输入数据怎么准备新安江模型需要两类最基本的时间序列逐时段面降雨量和蒸散发能力或实测蒸发皿蒸发量。在实际项目中降雨资料来自雨量站或雷达测雨蒸发则常用E601蒸发皿数据转换。为了模型稳定建议把时间步长统一日模型用一个值次洪模型一般用1小时或更短。数据格式我建议直接使用Excel或CSV原因是MATLAB读取方便而且同行之间交换数据也通用。表格至少包含三列时间、降雨量mm、蒸发量mm。如果做次洪模拟还需要起止时间对应的前期影响雨量或初始土壤含水量这些情况可以在数据文件中加一列初始状态说明。注意所有水量单位必须统一成毫米mm时间步长保持一致。这是新手最容易栽跟头的地方——降雨是毫米、流量是立方米每秒最后对不上账水量平衡一算全是漏洞。2.2 核心参数表与初值设置三水源新安江模型的参数不算多但每个参数对模拟结果的影响都不小。我把核心参数整理成下面这张表方便对照设置参数符号含义单位常见取值范围说明WM流域平均张力水蓄水容量mm100~200土壤能保持的最大水分量WUM / WLM / WDM上、下、深层张力水容量mm20/60/80左右三层分配总和等于WMKC蒸散发折算系数无量纲0.8~1.2蒸发皿换算到流域蒸散发C深层蒸散发系数无量纲0.1~0.2深层水分参与蒸散发的比例SM自由水蓄水容量mm10~50决定三水源划分比例EX自由水蓄水容量分布指数无量纲1.0~1.5越大说明空间分布越不均匀KI / KG壤中流/地下径流出流系数1/时段0.2~0.5 / 0.1~0.3两者之和一般小于1CS / CI / CG地面/壤中/地下汇流消退系数无量纲0.5~0.9 / 0.7~0.9 / 0.9~0.999越接近1消退越慢UH地面汇流单位线无量纲由流域特征确定用实测洪水反推或经验公式初始状态包括三层张力水蓄量WU、WL、WD以及自由水蓄量S。次洪模拟时初始值可以根据前期降雨估算或者用模型预热期自动稳定日连续模拟时从土壤偏干状态开始跑一两年系统自然会收敛到合理水平。2.3 代码模块怎么划分写MATLAB程序最忌讳一堆代码塞在同一个脚本里。我习惯把模型拆成两大部分一个是模型函数负责输入参数和气象序列、输出模拟流量另一个是主脚本负责读数据、调参数、画图、做率定分析。这样后续做参数自动优选、不确定性分析就非常方便。模型函数内部再按环节拆成子函数蒸散发模块、蓄满产流模块、三水源划分模块、汇流模块。每个子函数输入输出尽量只有必要变量方便单独测试。我见过不少人把所有计算写在同一个循环里最后参数率定的时候改一处就要全局排查效率特别低。3. 蓄满产流与湿润过程模拟核心计算环节详解3.1 三层蒸散发先消耗上层再挖下层最后动用深层蒸散发计算是湿润过程模拟的第一环很多初学者忽略它导致土壤含水量一直偏高产流偏大。三水源新安江模型用三层蒸散发模式上层张力水WU首先满足蒸散发如果WU不够再从下层WL扣除当下层也不够且WU已经耗尽时才按比例从深层WD扣除同时用深层蒸散发系数C限制最大消耗量。蒸散发能力EM的表达式是EM KC × E0其中E0是实测蒸发皿蒸发量。实际蒸散发量EU、EL、ED分三种情况计算当WU EM时EU EMEL 0ED 0当WU EM但WL C × (EM - WU)时EU WUEL EM - WUED 0当WU EM且WL C × (EM - WU)时EU WUEL WLED C × (EM - WU) - WL。这套规则背后是对干旱期土壤水分垂直分布的经验总结看起来简单但在MATLAB里用if判断实现时要注意顺序先判断上层够不够再判断下层能不能补最后才动用深层。我见过有人把条件写反导致湿润期蒸发被低估。3.2 蓄水容量曲线把流域的空间不均匀性塞进一个公式新安江模型引人入胜的地方在于它用一条蓄水容量曲线来描述流域内不同位置蓄水容量的差异。曲线公式是f / F 1 - (1 - WMM / WMM) ^ B实际工作中常用的是抛物线形式当蓄水容量小于某个阈值时对应的面积比例呈指数分布。简化处理时产流计算可以用下面的公式直接算若 PE 0不产流 若 PE W0 WMMR PE - WM × (1 - (PE W0) / WMM)^(1 B) W0 - ...完整推导见下这里我不展开全部推导直接给出编程时最常用的递推公式。假设流域平均张力水蓄量初始值为W0产流面积比FR计算公式为如果 W0 WMM 且 PE W0 WMM产流面积比 1 - (1 - (PE W0 - WMM) / (WM PE))^B 之类的修正式。说实话不同教材的公式写法略有差异关键是把握物理意义PE降雨扣除蒸发后的剩余水量先补充张力水蓄量蓄满的部分才产流蓄水容量曲线描述的就是“流域里有多少面积蓄满了”。在MATLAB中实现时建议用土壤含水量W的动态更新来隐含实现不必每次显式计算面积比这样代码更稳定。我实际使用的简化蓄满产流循环如下先计算净雨量 PE P - EM如果PE 0则只更新蒸发不产流。用当前张力水蓄量W0与最大蓄水容量WM比较计算产流量R。为了程序简化很多工程版本直接采用“如果W PE WM则R W PE - WM否则R 0”的“单点蓄满”模式。严格来说这只适用于面积均匀的极小区块流域尺度上要引入蓄水容量曲线。但如果你的资料不足以率定B参数用单点模式做初步模拟也能看趋势。3.3 湿润过程模拟土壤含水量如何逐时段更新回到“液体湿润模拟”这个主题。所谓湿润模拟在模型里就是逐时段更新张力水蓄量W和自由水蓄量S的过程。每次计算完产流后张力水蓄量更新为W_new W_old P - EM - R直观理解降雨和蒸发先影响土壤水库多余的水才变成径流被“挤”出去。自由水蓄量的更新则取决于分水源计算这部分在下一个小节展开。我建议在MATLAB中用数组记录每个时段的W和S哪怕最后输出不需要也要留着做过程分析和水量平衡验证。拿到湿润过程曲线后你能很直观地看到土壤在雨季反复蓄满、旱季逐步消耗这种过程分析对校核模型行为非常有用。3.4 自由水蓄水库与三水源划分代码实现自由水蓄水库的作用模型如下产流R先进入自由水蓄水库同时自由水还有前期蓄量S0。水库出流包括壤中流和地下径流出流量与当前蓄量成正比当蓄量超过SM时超出部分直接成为地面径流。分水源计算编程时用以下递推公式RS max(0, S_new - SM)RSS KI × (S_new - RS)RG KG × (S_new - RS)S_final S_new - RS - RSS - RG其中S_new表示产流进入后的自由水蓄量考虑前期蓄量。需要注意的是KI和KG的量纲是1/时段实际使用时要结合时段长做调整如果时段长从1小时变成1天系数必须重新率定。MATLAB代码实现这一段的典型写法是function [RS, RSS, RG, S_new] splitWater(R, S_old, SM, KI, KG) S_tmp S_old R; % 产流进入自由水蓄水库 RS max(0, S_tmp - SM); % 超蓄部分形成地面径流 S_tmp S_tmp - RS; RSS KI * S_tmp; % 壤中流出流 RG KG * S_tmp; % 地下径流出流 S_new S_tmp - RSS - RG; % 时段末自由水蓄量 end这套代码虽然看起来简单但它是整个三水源划分的发动机。实测中我经常把RS、RSS、RG分别存成数组后面做水源组成分析时直接画堆叠面积图一目了然。3.5 湿润过程与产流面积变化的联动效应三水源新安江模型里还有一个隐藏的联动效应产流面积不是固定比例而是随土壤湿润程度变化的。在湿润期土壤含水量高产流面积迅速扩大地面径流占比上升在干旱期产流面积缩小降雨多用来补充土壤水分地面径流很少。这个现象在MATLAB模拟中会自然地涌现出来因为W在雨季偏高产流R也随之偏高自由水蓄水库更容易蓄满RS占比自然增大。如果你做的湿润模拟结果里RS时间序列和W的变化趋势高度一致说明模型行为是对的如果出现W很低但RS很大的情况就要检查是不是蒸散发模块或初值出了问题。4. 汇流计算与模拟流量输出4.1 地面径流汇流单位线法三水源划分完三股水还不能直接相加成出口流量因为它们汇入河网的时间不同。地面径流最快通常用单位线法做汇流。单位线本质上是一个流域的“脉冲响应”表示一个单位净雨在出口断面形成的过程线。在MATLAB中单位线可以预先存储为一个向量UH然后对地面径流序列RS做卷积。标准实现方式如下QRS conv(RS, UH) * catchmentArea / dt / 1000;这里有个单位换算容易搞错RS单位是mm要转换成流量m³/s才能参与后面求和。换算系数是流域面积Fkm²除以时段长dt小时再乘以1000/3600之类建议单独写一个单位换算函数别在代码里散落魔法数字。单位线本身的推求方法有经验公式和实测资料反演两种。没有实测资料时可以用瞬时单位线或Clark单位线近似但参数需要根据流域面积、河长、坡降估算。4.2 壤中流与地下径流汇流线性水库消退壤中流和地下径流速度慢模型用线性水库来模拟。线性水库的思想是出流量等于当前蓄量乘以消退系数蓄量越大出流越快。递推公式很简洁Q(t) Q(t-1) × CI RSS(t) × (1 - CI) × 面积转换系数其中CI就是壤中流消退系数。地下径流同理用CG代替CI即可。这两个系数越接近1说明水库调节能力越强径流过程越平缓反之则越尖瘦。实现时注意一个细节汇流模块的输入是时段产流量输出是时段流量两者的时间对齐方式会直接影响模拟洪峰的相位。一般推荐产流序列按“时段末”输出汇流结果对应“时段末”流量这样后续与实测流量对比时不会出现系统性错位。4.3 总流量合成与过程线绘制三股汇流结果相加就得到流域出口断面的模拟流量Qsim QRS QRSS QRG然后把它和实测流量放在同一张图上对比再用Nash-Sutcliffe效率系数NSE和相对误差等指标评价模拟效果。MATLAB画流过程线的方式很简单plot(t, Qsim, r-, t, Qobs, k--); legend(模拟流量, 实测流量); xlabel(时间); ylabel(流量 (m^3/s));除了流量过程线我还建议大家把降雨、张力水蓄量、自由水蓄量、三水源拆分结果一起画出来。这样能看到整个流域的湿润状态和径流响应的对应关系对判断模型是否合理、参数是否灵敏非常有帮助。5. 实操过程中的高频问题与参数率定心得5.1 常见问题速查表症状可能原因排查方向模拟流量整体偏大蒸散发被低估或WM偏大检查KC取值、EM序列是否合理降低WM试算洪峰过高、过程尖瘦地面径流占比过大或汇流系数偏大调小SM或降低KI/KG让更多水走地下退水太慢CS或CG偏大降低消退系数加快退水基流几乎为零初始地下水库蓄量偏低或KG太小预热模型或抬高初始RG蓄量水量平衡严重不符单位换算错误或时间步长不一致逐项核对输入数据单位与循环步长NSE很低但过程趋势对初值设置不当导致前期偏差大增加预热期或直接采用实测前期土壤含水量这套表是我在实际项目中迭代出来的排查顺序一般都是先查数据单位再查边界条件最后才动参数。很多人一上来就调参数结果越调越乱。5.2 参数率定的实操建议三水源新安江模型的参数率定有两种路径手动试错和自动优化。手动试错适合刚上手时理解模型行为自动优化适合精度要求高、参数多的场景。MATLAB自带的优化工具箱可以用但我更常用的是简单的SCE-UA算法或者遗传算法因为它们在处理水文模型参数优选时比较成熟。无论用哪种方法都建议分步率定先率定产流参数WM、WUM、WLM、KC、C再率定分水源参数SM、EX、KI、KG最后率定汇流参数单位线、CS、CI、CG。分步率定能大幅度减少参数之间的相互干扰收敛速度也快得多。另外参数率定不能只看NSE一个指标。我习惯同时看洪峰相对误差、洪量相对误差和过程线形态。有些参数组合能让NSE很高但洪峰明显偏低这种结果在洪水预报实战中并不好用。5.3 水量平衡检查模型是否合理的底线水量平衡是检验模型实现是否正确的最可靠手段。连续模拟时整个模拟期的总降雨应当等于总蒸散发加总径流再加土壤蓄量变化量。即使简化模型这个等式也应该基本成立误差一般控制在2%以内。在MATLAB里我通常会在主脚本最后加一段自动校验totalP sum(P); totalE sum(EM); totalR sum(Qsim) * dt * 3600 / (catchmentArea * 1e6) * 1000; dW W_end - W_start; balanceError totalP - totalE - totalR - dW; fprintf(水量平衡误差: %.4f mm\n, balanceError);如果误差超过2%先别急着调参数回头查代码逻辑和单位换算。5.4 MATLAB运行的几个容易踩的坑最后说几个MATLAB环境下的实际坑。第一循环次数多的时候不要用动态增长的数组提前用zeros或NaN初始化可以省非常多时间第二反复调用模型做率定时尽量把模型封装成函数而不是脚本函数可以避免工作区变量污染而且方便用parfor并行率定第三读Excel数据时注意日期格式建议统一用datetime类型否则时间轴对不上画图全是乱的第四MATLAB 2020及以上版本对中文注释支持没问题但函数文件名不要用中文否则偶尔会有奇怪的编码问题。我在实际项目中还养成了一个习惯把参数、输入数据、输出结果全部写成结构化组织并用一个配置文件统一管理。这样做的好处是换流域、换时段时只需要改配置文件模型函数基本不动可复用性很强。根据我个人经验三水源新安江模型在MATLAB里的实现难度其实不高真正考验人的是对物理过程的理解和数据处理的基本功。把蓄满产流、三层蒸散发、自由水蓄水库这几个核心逻辑吃透再配合规范的数据管理和参数率定流程哪怕一个小型流域的洪水预报系统也完全能用这套代码搭建起来。后续如果大家有兴趣我可以再把自动率定、不确定性分析比如GLUE方法和分布式扩展的版本整理出来分享。