企业官网建设流程全解析

每次做水文预报，最绕不开的就是新安江模型。我最早接触这个概念是在做中小流域洪水预报的时候，那时手里只有日降雨量和蒸发资料，要估计一场洪水的产流量，翻了很多教材，最后发现真正能落地、能写进程序、参数含义又清晰的，还是三水源新安江模型。这个模型用matlab来实现，代码量不大，但每一步的水量平衡关系都清清楚楚，特别适合入门水文模型的人当作第一个练手项目。

这篇文章就是我在matlab里实现三水源新安江模型的全过程记录。我会把模型怎么拆、每个模块的公式怎么转成代码、参数怎么定、结果怎么调，一条线讲清楚。不管你是做洪水预报、水资源评价，还是单纯想搞懂蓄满产流和三种水源是怎么被分开的，这篇都能给你一个可以直接跑起来的起点。

1. 模型核心逻辑：为什么湿润地区产流要先算蓄满

三水源新安江模型最本质的一个假设是蓄满产流。这和超渗产流是两条完全不同的路子。超渗产流看的是雨强和地面下渗能力谁大，而蓄满产流看的是包气带土壤含水量有没有到田间持水量。在湿润地区，比如南方大部分流域，降雨频繁且总量大，包气带容易蓄满，所以蓄满产流的假设更贴合实际。

那“蓄满”意味着什么？意味着土壤不再能继续存住水分，后续降雨直接在地表形成径流。新安江模型的巧妙之处在于，它把整个产流过程拆成了两步：先算整个流域有多少面积蓄满了，蓄满的面积上产生径流；再把这部分径流用一个自由水蓄水库去分水源——分成地面径流、壤中流和地下径流。这就是“三水源”的来历。

用matlab实现的时候，整场降雨要按小时或天来一步步推。每一时段里，降雨先进入张力水蓄水层，蓄满之后的水进自由水蓄水层，自由水层也满了之后就产流。这个逻辑顺序如果理不清楚，后面代码很容易写得混成一团。我的建议是先把下面这张水量平衡表在草稿纸上画出来，再动手写代码：

这里有个新手容易忽略的点：张力水和自由水是两种不同“性质”的水。张力水被土壤颗粒吸附，很难自由流动，只有当超过田间持水量后才会变成自由水。所以在程序里，这两个蓄水层要分开设置容量、分开计算，不能图省事合并成一个变量。

整场降雨过程中，每个时段的蒸散发损失优先从张力水中扣除。这个顺序也重要，后面对照实测流量过程线时，如果模拟的洪峰偏大，往往就是蒸散发扣得太少，土壤一直处在蓄满状态。

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计算公式涉及一个中间变量A：

A = 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 - RG

RS是地面径流，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.4，EX常见取值在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，水库的调蓄能力越强，出流越平缓。这个函数虽然只有一行，但它是汇流计算里的基本积木。地面径流用CS=0.8到0.9，壤中流用CS=0.9到0.98，地下径流用CS=0.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左右，干旱地区可以到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按同一个比例更新三层张力水蓄量。这是不严格的。严格做法应该把三层分别做水量平衡，上层优先蓄满再进入下层。但工程上很多简化版新安江模型代码就是这么写的，对总产流的影响不显著，因为产流总量取决于WU+WL+WD的总和，而比例更新法保持了总和正确。

运行这段程序之后，你会得到一组流量过程线。如果数据来自实际流域，就把模型输出和实测流量画在一起看拟合效果。我常用的评价指标是两个：确定性系数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应等于张力水蓄满的面积比例，也就是产流面积比。这一步如果省略，模型在干旱时期的模拟误差会明显偏大。我的建议是，刚开始练手时可以用FR=1的简化版，但真正做研究时，要把产流面积比例的计算加回去。

第四个坑是matlab里for循环的效率和数值稳定性。如果你模拟的时间序列很长，比如几十年日尺度数据，for循环里有大量重复计算。可以在循环之前把所有常数参数预先算好，比如WMM = WM * (1+B)、SMM = SM * (1+EX)，避免每次迭代都重复计算。另外，代码里所有除法都要加一个小量保护，防止除零。比如计算ratio的时候，我用max(W0, 1e-6)来避免分母为零。

就个人调试经验而言，三水源新安江模型是水文模型里“性价比”最高的一个。它的物理概念清晰、代码实现简洁、参数不是特别多，但模拟效果在湿润地区普遍不错。用matlab把它写一遍，你对蓄满产流、水源划分、汇流过程的理解会从“看书懂了”变成“真正会用了”。而且这套代码框架，后面改造成分布式新安江模型、或者耦合数据同化模块，都是现成的底子。