目录
前言
新安江模型是我国应用最广泛的概念性流域水文模型。本文根据模型原理，编制各个模块程序，最后组合成完整的模型程序，使程序更加规范简明。
上图为新安江模型计算流程图，线上的变量为参数。
一、蒸散发计算
土壤在垂向分为三层，用三层蒸发模型计算蒸散发量。
% evap 三层蒸发模型
% 输入：降雨量p，水面蒸发e0，初始土壤含水量wu,wl,wd
% 参数：蒸散发折算系数kc,下层张力水容量lm,深层蒸散发折算系数c
% 输出：各层蒸发eu,el,ed，总蒸发e，净雨量pe
function [eu,el,ed,e,pe]=evap(p,e0,wu,wl,wd,kc,lm,c)
ep=kc*e0;
if wu+p>ep
eu=ep;el=0;ed=0;
elseif wl>c*lm
eu=wu+p;el=(ep-eu)*wl/lm;ed=0;
% 只有当(ep-eu)>lm时会出现el>wl的情况
elseif wl>c*(ep-wu-p) % 这里不能是(ep-eu)
eu=wu+p;el=c*(ep-eu);ed=0;
else
eu=wu+p;el=wl;ed=c*(ep-eu)-el;
if ed>wd
ed=wd;
end
end
e=eu+el+ed;
pe=p-e;二、产流量计算
采用蓄满产流模型计算产流量，需注意的是：
（1）只有当PE>0时才会有产流量，否则产流量R为0。
（2）有人会如流程图中在不透水面积比IM上计算产流量RB再加到R中，不需要。
% runoff 蓄满产流模型
% 输入：净雨量pe，初始土壤含水量w
% 参数：流域平均张力水容量wm，张力水蓄水容量曲线方次b，不透水面积比例im
% 输出：产流量r
function r=runoff(pe,w,wm,b,im)
if pe>0 % 只有pe>0时才产流
mm=wm*(1+b)/(1-im);
a=mm*(1-(1-w/wm)^(1/(1+b)));
if pe+a<mm
r=pe-wm+w+wm*(1-(pe+a)/mm)^(1+b);
else
r=pe-wm+w;
end
else
r=0;
end三、土壤含水量更新
降雨量、蒸发量、产流量以及土壤含水量的变化量之间存在水量平衡关系。利用前面计算得到的蒸发量及产流量，进行各层土壤含水量的更新。多余的水量先补充上层，再补充下层及深层。
这里单独用一个模块来计算土壤含水量的变化与更新，使程序更加简明易读。
% update_w 土壤水更新
% 输入：初始土壤含水量，各层蒸发，降雨量，产流量
% 参数：各层含水能力um,lm,dm
% 输出：各层土壤含水量和总量
function [wu,wl,wd,w]=update_w(wu0,wl0,wd0,eu,el,ed,p,r,um,lm,dm)
wu=wu0+p-eu-r;
wl=wl0-el;
wd=wd0-ed;
% 检查土壤水超量
if wu>um
wl=wl+(wu-um);wu=um;
if wl>lm
wd=wd+(wl-lm);wl=lm;
if wd>dm
wd=dm;
end
end
end
w=wu+wl+wd;四、分水源计算
产流量R进入*水蓄水库，根据出流情况划分为地面径流RS、壤中流RI和地下径流RG。三分水源计算是新安江模型中较难理解和编程易出错的部分，我总结了有以下要点：
（1）底宽的变化
*水蓄水库的底宽为产流面积比FR，并且FR是随时段而变化的。可以理解为S*FR为实际蓄量，或者说RS、RI、RG都是发生在产流面积比FR之上。
（2）差分误差的处理
为了减小因向前差分所造成的误差，每个计算时段的入流量R，按5mm为一段划分为N段进入水库进行水量平衡计算，计算时段DT也分为了N段。
（3）参数的时段转换
参数KI，KG，包括后面的CI、CG、CR，都是以日（24h）为时段长定义的，它们需要根据实际计算时段长而改变。
% separate_r 三分水源
% 输入：净雨pe，产流r，时段初*水蓄量s，产流面积比fr
% 参数：不透水面积比im，*水蓄水容量sm，*水蓄水容量曲线方次ex，
% 壤中流出流系数ki，地下水出流系数kg
% 输出：三水源rs,ri,rg,时段末s,fr
function [rs,ri,rg,fr,s]=separate_r(pe,r,fr,s,im,sm,ex,ki,kg)
ms=sm*(1+ex);
if pe>0
x=fr;
fr=(r-im*pe)/pe; %扣除不透水面积比来计算产流面积比
s=s*x/fr; %s有变化底宽fr的影响
ss=s; %保存初始s用于后面计算
q=r/fr; %产流面积比上的产流深
% 处理向前差分误差
n=fix(q/5)+1; %按5mm一段划分n段
q=q/n;
kid=(1-(1-(ki+kg))^(1/n))/(1+kg/ki); %出流系数按时段长改变
kgd=kid*kg/ki;
rs=0;ri=0;rg=0;
for i=1:n
if s>sm;s=sm;end;
au=ms*(1-(1-s/sm)^(1/(1+ex)));
if q+au<ms
rs=fr*(q+s-sm+sm*(1-(q+au)/ms)^(1+ex))+rs;
else
rs=fr*(q+s-sm)+rs;
end
s=s+q*i-rs/fr;
ri=kid*s*fr+ri;
rg=kgd*s*fr+rg;
s=ss+q*i-(rs+ri+rg)/fr;
end
else
rs=0;
ri=s*ki*fr;
rg=s*kg*fr;
s=s*(1-ki-kg);
end说明：
（1）在计算产流面积比FR时，需要考虑不透水面积比IM的影响，所以仍需要引入此参数
（2）输入和输出的S采用同一个变量，FR也是如此，在差分误差的处理计算中，S、RS、RI、RG是不断更新变化的
（3）这里的KI、KG是已经根据计算时段DT转换后的值，在后面调用此模块的程序中会有所体现。在这个程序段中，可以看到KI、KG仍需要根据DT分段后的新时段进行再一次转换。
五、汇流计算
（1）坡地汇流
采用线性水库法，将RS、RI、RG转变为QS、QI、QG。相应的消退系数参数分别为CS、CI、CG，一般CS取0，所以略去参数CS。
（2）单元面积河网汇流
单元面积河网总入流QT=QS+QI+QG，采用滞后演算法计算单元面积河网汇流QTR。
（3）河道汇流
采用马斯京根分段连续演算法，计算QTR由单元面积出口至流域总出口的河道汇流过程。
% conflu 三水源汇流计算
% 输入：三水源产流rs,ri,rg，流域面积F
% 参数：壤中地下消退系数ci,cg,河网水流消退系数cr,滞时L,马斯京根参数k,x,n
% 输出：流量q
function [qt,q]=conflu(rs,ri,rg,F,cs,ci,cg,cr,L,k,x,n,qt0)
dt=k;
U=F/(3.6*dt); %单位换算系数
%坡地汇流
qs(1)=rs(1)*(1-cs)*U; %初始流量计算
qi(1)=ri(1)*(1-ci)*U;
qg(1)=rg(1)*(1-cg)*U;
qt(1)=qs(1)+qi(1)+qg(1);
for i=2:length(rs);
qs(i)=qs(i-1)*cs+rs(i)*(1-cs)*U;
qi(i)=qi(i-1)*ci+ri(i)*(1-ci)*U; %线性水库法
qg(i)=qg(i-1)*cg+rg(i)*(1-cg)*U;
qt(i)=qs(i)+qi(i)+qg(i);
end
%单元面积河网汇流
exn=10; %扩展
qtt=[ones(1,exn)*qt0,qt];
qtr=zeros(size(qtt));
for i=(exn+1):length(qtr)
qtr(i)=cr*qtr(i-1)+(1-cr)*qtt(i-L); %滞后演算法
end
qtr=qtr(exn+1:end);
if n>0
%单元面积以下河道汇流
c0=(0.5*dt-k*x)/(0.5*dt+k-k*x);
c1=(0.5*dt+k*x)/(0.5*dt+k-k*x);
c2=(k-0.5*dt-k*x)/(0.5*dt+k-k*x);
qq=zeros(length(qtr)+1,n+1);
qq(2:end,1)=qtr';
%用分段马斯京根法
for j=2:size(qq,2)
for i=2:size(qq,1)
qq(i,j)=c0*qq(i,j-1)+c1*qq(i-1,j-1)+c2*qq(i-1,j);
end
end
q=qq(2:end,n+1);
else
q=qtr';
end
qt=qtt(exn:end-1);六、单元的综合计算
新安江模型根据雨量站的分布，采用泰森多边形法，将整个流域分为若干单元。每个单元将上面的各个模块串联起来计算，便可以得到该单元对流域出口的流量过程。
这里输入、输出的降雨、蒸发、产流等，都是时间序列的数组。前面模块函数中定义的是单个值，所以在各计算时段调用这些模块得到时间序列值。该程序段的流程包括：参数的读取，参数的时段转换，初值的设置，产流计算，汇流计算。
% 新安江模型主模块
% 输入：降雨，水面蒸发，参数，初值，面积权重，至出口断面河段数
% 输出：流域蒸发，产流，土壤含水量，流量
function [e,r,w,wu,wl,wd,fr,s,qt,q] = xaj_mdl(p,e0,Para,S0,FW,NE)
%disp('Input parameters ...');
KC = Para(1) ; UM = Para(2) ; LM = Para(3) ; C = Para(4);
WM = Para(5) ; B = Para(6) ; IM = Para(7) ;
SM = Para(8) ; EX = Para(9) ; KI = Para(10); KG = Para(11);
CI = Para(12); CG = Para(13); CR = Para(14); LT = Para(15);
XE = Para(16); KE = Para(17); F = Para(18)*FW;
DM=WM-UM-LM; CS=0; DT=KE; % time step (h)
KIt=(1-(1-KI-KG)^(DT/24))/(1+KG/KI);
KGt=KIt*KG/KI;
CIt=CI^(DT/24);
CGt=CG^(DT/24);
CRt=CR^(DT/24);
%disp('Calculate initial state ...');
% initial value
[wu(1),wl(1),wd(1),w(1),fr(1),s(1),qt0]=ini_state(S0,WM,UM,LM,DM,SM);
%disp('Calculate runoff product ...');
for i=1:length(p)
[eu,el,ed,e(i),pe]=evap(p(i),e0(i),wu(i),wl(i),wd(i),KC,LM,C);
r(i)=runoff(pe,w(i),WM,B,IM);
[wu(i+1),wl(i+1),wd(i+1),w(i+1)]=update_w(wu(i),wl(i),wd(i),eu,el,ed,p(i),r(i),UM,LM,DM);
[rs(i),ri(i),rg(i),fr(i+1),s(i+1)]=separate_r(pe,r(i),fr(i),s(i),IM,SM,EX,KIt,KGt);
end
%disp('Calculate flow confluence ...');
[qt,q]=conflu(rs,ri,rg,F,CS,CIt,CGt,CRt,LT,KE,XE,NE,qt0);
end七、总流域的综合计算
前面得到的是计算单元的结果，将各个计算单元结果综合得到总流域的结果。
function [re_s,re_dt,re_val] = xaj_cal(dh,ST,ET,Pdat,Edat,Qdat,Para,Zdat,Sdat,Sdat0)
DT=Para(17); F=Para(18); ZN=size(Zdat,1); % 时段，面积，分块数
%disp('Extract data ...');
[T,TL,p,e0,qob,Tq,TLq,Qobq]=extra_dat(dh,ST,ET,DT,Pdat,Edat,Qdat,ZN);
N=length(e0);
% 分区域计算结果------------------------------------------------------------
e_z=[]; r_z=[]; w_z=[]; q_z=[]; re_s=[];
for zn=1:ZN
if isempty(Sdat)
disp(['zone ',num2str(zn),' has no initial state value, use default']);
S0=Sdat0;
else
idx=find(Sdat(:,1)==str2num(datestr(ST,'yyyymmdd'))&Sdat(:,2)==zn);
if isempty(idx)
disp(['zone ',num2str(zn),' have no state value, use default']);
S0=Sdat0;
else
S0=Sdat(idx,3:end);
end
end
[e,r,w,wu,wl,~,fr,s,qt,q] = xaj_mdl(p(:,zn),e0,Para,S0,Zdat(zn,1),Zdat(zn,2));
e_z=[e_z,e'];r_z=[r_z,r'];w_z=[w_z,w(2:end)'];q_z=[q_z,q]; %保存各分区结果
if (dh=='d')
re_s=[re_s;TL',zn*ones(size(TL')),w(1:end-1)',wu(1:end-1)', ...
wl(1:end-1)',fr(1:end-1)',s(1:end-1)',qt']; %日模型保存初始状态值
end
end
% 计算流域平均结果----------------------------------------------------------
zn_w=Zdat(:,1); % 面积权重矩阵（zn*1）
p_m=p*zn_w; e_m=e_z*zn_w; r_m=r_z*zn_w; w_m=w_z*zn_w; qcal=sum(q_z,2);
% 计算特征值并显示----------------------------------------------------------
[r_cal,r_obs,r_er,qm_cal,qm_obs,qm_er,Dc]=obj_fun(TL,qcal,TLq,Qobq,DT,F);
p_t=sum(p_m); e_t=sum(e_m); r_t=3.6*sum(qcal)*DT/F; qm_t=max(qcal); %计算期总量
disp(strcat('Total-----------------------Rainfall (mm): ',num2str(p_t)));
disp(strcat(' Evaporation (mm): ',num2str(e_t)));
disp(strcat(' Runoff (mm): ',num2str(r_t)));
disp(strcat(' Peak (m3/s): ',num2str(qm_t)));
disp(strcat('Measure Time---------Observed Runoff (mm): ',num2str(r_obs)));
disp(strcat(' Calculated Runoff (mm): ',num2str(r_cal)));
disp(strcat(' Runoff Difference: ',num2str(r_er),'%'));
disp(strcat(' Observed Peak (m3/s): ',num2str(qm_obs)));
disp(strcat(' Calculated Peak (m3/s): ',num2str(qm_cal)));
disp(strcat(' Peak Difference: ',num2str(qm_er),'%'));
disp(strcat('--------------------------N-S Coefficient: ',num2str(Dc)));
disp(' ');
re_val=[TL(1),TL(end),p_t,e_t,r_t,qm_t,r_obs,r_cal,r_er,qm_obs,qm_cal,qm_er,Dc];
re_dt=[TL',p_m,e_m,r_m,w_m,qcal,qob'];
% 画图---------------------------------------------------------------------
figure;
if isempty(Tq)
subplot(3,1,1); bar(1:N,p_m); ylabel('P（mm）');
set(gca,'XTick',[]); set(gca,'YDir','reverse'); %对Y方向反转
title([num2str(TL(1)),' - ',num2str(TL(end)),' hydrograph']);
subplot(3,1,2:3); plot(T,qcal,'b-'); ylabel('Q（m^3/s）');
legend('Calculated');
else
subplot(3,1,1); bar(1:N,p_m); ylabel('P（mm）');
set(gca,'XTick',[]); set(gca,'YDir','reverse'); %对Y方向反转
title([num2str(TL(1)),' - ',num2str(TL(end)),' hydrograph']);
subplot(3,1,2:3); plot(Tq,Qobq,'r.',T,qcal,'b-'); ylabel('Q（m^3/s）');
legend('Observed','Calculated');
end
end说明：
（1）这里输入的降雨是二维数组，即各雨量站（计算单元）、各时段的值。
（2）该程序段的流程可概括为：读取参数，提取所选时间的资料，状态初值的处理，新安江模型主模块的运行，流域结果的综合，特征值计算，画图。
（3）该程序段中还用到其他模块程序，后文说明。
八、其他模块
以上所述的是新安江模型的核心原理及模块程序。对于完整的程序系统，还有其他功能需实现，所以还有其他模块程序编制。
上图为各程序文件，简述如下：
（1）输入和输出数据库
输入：input.xlsx，存放日模型和次洪模型的参数、流域分区面积权重以及到流域出口的河段数、洪水场次、降雨、蒸发、实测流量资料
输出：result.xlsx，存放日模型状态结果（作为初值调用），以及日模型库中没有结果时的默认初始值，日模型和次洪模型的时段结果（包括时间、降雨、蒸发、产流、土壤含水量、计算流量、实测流量），日模型和次洪模型的结果特征值（包括计算起始时间、降雨量、蒸发量、径流量、洪峰值、径流误差、洪峰误差、峰现时差、确定性系数）
（2）核心模块
核心模块程序在上文中已讲述，包括evap（蒸散发计算），runoff（产流量计算），update_w（土壤含水量更新），separate_r（分水源），conflu（汇流计算），xaj_mdl（单元的综合），xaj_cal（流域的综合）
（3）其他模块
extra_dat：根据设置的起止时间从数据库提取降雨、蒸发及实测流量资料，并进行缺测值的识别与处理
ini_state：初值的读取与判定
obj_fun：计算相对误差、确定性系数等特征值
XAJ_DP：主运行程序，包括日模型或次洪模型的选择，数据的读取，场次洪水的选择，新安江模型的运行，数据的保存
九、结果展示
选择几场洪水的运行结果如下：
模型参数表
运行的屏幕输出
出图
结果的数据输出
总结
完整的新安江模型程序系统是一个复杂的系统，体现在：
（1）模型分为蒸散发计算、产流计算、水源划分、汇流计算这四个层次。
（2）为了考虑降雨空间分布不均匀性，将流域分为若干单元计算。
（3）初值的处理以及状态值的保存，特征值的计算。
（4）分为日模型和次洪模型的计算。
（5）洪水场次的选择，数据的提取，缺测数据的处理，数据的保存。
这里采用模块化的Matlab编程，实现了上述功能，并使程序更加规范。