编程语言
首页 > 编程语言> > 新安江模型的Matlab程序实现

新安江模型的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