标签：disp end Para 程序实现 wl wu Matlab 计算 新安江
目录
前言
新安江模型是我国应用最广泛的概念性流域水文模型。本文根据模型原理，编制各个模块程序，最后组合成完整的模型程序，使程序更加规范简明。
上图为新安江模型计算流程图，线上的变量为参数。
一、蒸散发计算
土壤在垂向分为三层，用三层蒸发模型计算蒸散发量。
% 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编程，实现了上述功能，并使程序更加规范。
标签：disp,end,Para,程序实现,wl,wu,Matlab,计算,新安江 来源： https://blog.csdn.net/weixin_53441521/article/details/117383531
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。