[学习分享]基于matlab的新安江模型_01_模型介绍与蓄满产流

写在前面的 

  最近笔者刚完成水文预报这门课的课程设计，课程设计要求根据课本自行实现新安江模型，完成径流模拟。现在课程设计已经基本全部做完，自己感觉做的也还不错，同时也因为蛮喜欢水文预报这门课的，所以想再对课程设计的整个过程做个整理分享出来，也希望能够帮助到一些困惑于课程设计的同学。

  需要注意的是，本文构建的新安江模型仅是一个简易版本的新安江模型。主要体现为：在分水源计算中只划分了地面径流和地下径流；在坡面汇流中地上部分径流未直接转换法（即没有考虑其汇流过程）；在参数方面仅进行了简易的人为调参，并未进行更进一步的参数调整。

  同时笔者代码能力和语言表达能力有限，可能存在部分错误，如果有发现，也麻烦大家多多指正，谢谢！

内容大纲

  本文完成了从半分布式和集总式两个方面进行新安江模型的构建，模型各个部分具体使用的模块如下：

流程	集总式	半分布式
数据预处理	计算泰森多边形	计算泰森多边形
	计算面雨量
蓄满产流计算	根据E0推求E，并进行三层蒸散发计算	根据E0推求EP，并进行三层蒸散发计算
	根据蒸散发量E及面雨量推求净雨量PE	根据蒸散发量E及各子流域降雨量推求各子流域净雨量PE
	根据净雨量及EU,EL,ED计算得出本时刻径流深及下一时刻初的土壤含水量W	根据净雨量及EU,EL,ED计算得出本时刻径流深及下一时刻初的土壤含水量W
	循环上述三步骤计算得出整个流域每个时刻的径流深R	循环上述三步骤计算得出各个子流域每个时刻的径流深Ri
分水源计算	根据蒸发能力FC及净雨量PE将径流深分为地面径流Rs及地下径流Rg	根据蒸发能力FC及净雨量PE将径流深分为地面径流Rs及地下径流Rg
坡面汇流	地上径流使用单位线转换法直接将Rs转化为Qs 地下径流使用线性水库法将Rg转化为Qg	地上径流使用单位线转换法直接将Rs转化为Qs 地下径流使用线性水库法将Rg转化为Qg
河网汇流	由Qs和Qg线性相加计算Qt，再根据滞后演算法计算模拟流量Q	对每个子流域根据马斯京跟法分段演算法计算其在出口断面处流量并线性相加得到模拟流量Q
精度评定	绘制逐年实测径流和模拟径流及降水图并计算确定性系数DC及均方根误差RSME	绘制逐年实测径流和模拟径流及降水图并计算确定性系数DC及均方根误差RSME

  本文的公式参考主要来源于包为民老师的《水文预报》第五版，后续的原理、公式等若不做额外说明则均摘自于该书。

新安江模型介绍

  新安江模型是华东水利学院(现河海大学)赵人俊教授团队提出的一个水文模型，是中国少有的一个具有世界影响力的水文模型。新安江模型是集总式水文模型（划分子流域时成为分散式水文模型），可用于湿润地区与半湿润地区的湿润季节。当流域面积较小时，新安江模型采用集总式模型，当面积较大时，采用半分布式模型。它把全流域分为许多块单元流域，对每个单元流域作产汇流计算，得出单元流域的出口流量过程。再进行出口以下的河道洪水演算，求得流域出口的流量过程。把每个单元流域的出流过程相加，就求得了流域的总出流过程。

蓄满产流

  蓄满产流的基本原理是:土壤所能够容纳的水量是一定的，当土壤蓄水量达到一定程度时，降落到土壤上的水一部分将不会再进行下渗而是产生径流，即“蓄满产流”。

   三层蒸散发蓄满产流的基本流程为（以逐日为例）：

   (1)首先根据土壤上层含水量WU、降雨P和流域蒸发能力Ep，进行三层蒸散发推求土壤上中下层的土壤蒸发:

    由于本课程设计蒸发数据为蒸发皿蒸发量E0，因此想要得到流域蒸发能力Ep需要对E0乘上蒸发折算系数Kc：

  Ep=Kc·E0

    折算后进行三层蒸散发计算，原理如下：

    计算公式如下：

   (2)根据(1)得到的总蒸发量计算净雨量PE。

PE=P-E

   (3)根据降雨径流关系计算产流量：

   (4)通过上述计算可以得到某一日的三层蒸发量Eu,El,Ed、净雨量PE、产流量R。如果这一日的PE>0，则这一日土壤含水量的变化量为：当日降雨量经过蒸发再经产流所剩余的水量。因此，计算产流后剩余的水量，三层土壤的含水量从上到下三层依次得到补充。但如果这一日的PE<=0，则这一日的土壤含水量变化量为：当日蒸发蒸发完当日降水量后所剩的蒸发量。同样的，剩余土壤的三层土壤的含水量从上到下三层依次进行蒸发。最后得到这一日结束时的土壤含水量作为下一日开始的土壤含水量，再从(1)开始进行计算。

   后续的代码将按照上述四个步骤循环计算，直到完成所有计算，最终对于每一次蓄满产流计算均可以得到该蓄满产流模型产流量计算表：

  对于上述步骤的matlab代码实现：

function [jieguo]=xuman(p,e0,w_chushi,wm,wum,wlm,wdm,b,k,c)
%p,e0请以行向量的形式输入
%w_chushi，请按照wm wu wl的形式以一个1x3的行向量的形式输入
%eg:w_chushi=[0 0 50];
%输如初始值
%wm=120;wum=15;wlm=85;wdm=20;b=0.3;k=0.95;c=0.14;
%读取数据并赋值
ep=e0*k;
w=w_chushi;%初始含水量
%%土壤含水量记录器
w_recoder=w;
%计算最大值a
wmm=wm*(1+b);
%设置各个变量的初始值
[x]=length(p);
eu=zeros(x,1);el=zeros(x,1);ed=zeros(x,1);e=zeros(x,1);r=zeros(x,1);
%计算三层蒸散发
for i=1:length(p)
    if w(1)+p(i)>=ep(i)
        eu(i)=ep(i);
        el(i)=0;
        ed(i)=0;
    elseif w(1)+p(i)<ep(i) 
        if w(2)>=c*wlm
            eu(i)=w(1)+p(i);
            el(i)=(ep(i)-eu(i))*w(2)/wlm;
            ed(i)=0;
        elseif  c*(ep(i)-w(1)+p(i))<=w(2) && w(2)<c*wlm
            eu(i)=w(1)+p(i);
            el(i)=c*(ep(i)-eu(i));
            ed(i)=0;
        elseif w(2)<c*(ep(i)-w(1)+p(i))
            eu(i)=w(1)+p(i);
            el(i)=w(2);
            ed(i)=c*(ep(i)-eu(i))-el(i);
        end
    end
    %计算总蒸散发及pe
    e(i)=eu(i)+el(i)+ed(i);
    pe(i)=p(i)-e(i);
    %径流
    a(i)=wmm*(1-(1-(sum(w)/wm))^(1/(1+b)));
    if pe(i)<=0
        r(i)=0;
    else
        if a(i)+pe(i)<=wmm
            r(i)=pe(i)+sum(w)-wm+wm*(1-(pe(i)+a(i))/wmm)^(b+1);
        else
            r(i)=pe(i)-(wm-sum(w));
        end
    end
    %考虑降水与蒸发对土壤含水量变化
    if w(1)+pe(i)-r(i)>=0
        w(1)=w(1)+pe(i)-r(i);
    else 
        if w(2)+pe(i)-r(i)>=0
            w(2)=w(1)+pe(i)+w(2)-r(i);
            w(1)=0;
        else 
            w(3)=w(1)+w(2)+w(3)+pe(i)-r(i);
            w(1:2)=0;
        end
    end
    %因为每一层土壤含水量都有最大值，保证含水量不会超过最大值
    if w(1)>wum
        w(2)=w(2)+w(1)-wum;
        w(1)=wum;
    end
    if w(2)>wlm
        w(3)=w(3)+w(2)-wlm;
        w(2)=wlm;
    end
    if w(3)>wdm
        w(3)=wdm;
    end
    w_recoder=[w_recoder;w];%记录
end
jieguo(:,1)=p;
jieguo(:,2)=e0;
jieguo(:,3)=ep;
jieguo(:,4)=eu;
jieguo(:,5)=el;
jieguo(:,6)=ed;
jieguo(:,7)=e;
jieguo(:,8)=pe;
jieguo(:,9)=w_recoder(1:length(w_recoder)-1,1);
jieguo(:,10)=w_recoder(1:length(w_recoder)-1,2);
jieguo(:,11)=w_recoder(1:length(w_recoder)-1,3);
jieguo(:,12)=sum(w_recoder(1:(length(w_recoder)-1),:),2);
jieguo(:,13)=r;
end

    

  该函数最终输出的jieguo矩阵以步骤⑤中的表形式保存排列，大家可以自行代入该表的数据进行计算以验证代码计算的正确性。