新安江模型Matlab代码

% 两水源新安江模型程序核心源代码

function [objqqrr]=XAJ_fun（xajPara）
%obj为模拟的各目标函数值，qq为时段实际流量与模拟流量，rr为时段地面净雨地下净雨总净雨
% XAJ是新安江的运行程序用于参数优化算法调用也用于新安江模型的预报
% xajPara是调用的优化参数;
global DATA
uh=load（‘.\data\uh.txt‘）;%导入地面径流汇流单位线，为一列或一行数据
runoffPara=load（‘.\data\runoffPara.txt‘）;%读入优化好的产流参数
startn=100;%统计预报结果起算时段序号*****************
% 输入起始值 WWUWLWDQG
WU=10;WL=61;WD=20;  % WU-流域初始上层张力水蓄量;WL-流域初始下层张力水蓄量;WD-流域初始深层张力水蓄量；
AREA=12200;  % 流域面积，km2;
N=6;   % 一个时段长（h）；
U=AREA/3.6/N;   % U-单位换算系数
W=WU+WL+WD;   % 流域初始蓄水量mm;

%输入雨量P蒸散发能力EM实测流量QS
TIME=DATA（:1）;
P=DATA（:2）;   % 输入雨量P
EM=DATA（:3）;  % 流域日蒸散发能力EM
QS=DATA（:4）;  % 实测流量QS

TRG0=0.4.*QS（1）;
%=================================
%产流参数
K=runoffPara（1）;%流域蒸发折算系数，是流域蒸散发能力与蒸发皿蒸发量之比
FC=runoffPara（2）;% 稳渗率
WUM=runoffPara（3）;% 流域平均上层蓄水容量 mm ；
WLM=runoffPara（4）;% 流域平均下层蓄水容量 mm ；
WDM=runoffPara（5）;% 流域平均深层蓄水容量 mm
WM=WUM+WLM+WDM;
IMP=runoffPara（6）;  % % 不透水面积占全流域面积的百分比，％，一般为0.01-0.02:[3664]p151倒数第一行；
B=runoffPara（7）;% 流域蓄水容量分布曲线指数；
C=runoffPara（8）; % 流域蒸发扩散系数 OR 深层蒸散发系数???
%------------------------------
%汇流参数
KKG=xajPara（1）; % 地下水消退系数；
CS=xajPara（2）;% CS－河网蓄水消退系数；
L=xajPara（3）;% L－滞后时间h；
L=round（L）;
KKS=xajPara（4）;% 地面径流消退系数
%=================================
WMM=（1+B）.*WM/（1-IMP）;  % 是流域内点最大蓄水容量 ；
M=size（P1）;   % M - 降雨数据的数量
PE=P-K.*EM;  % PE-扣除雨期蒸发后的降雨量mm;[3664]p145
for T=1:M               %%  T以时段为单位计算
   %===========================================
    %以下为产流计算;输入数据：PE，W； 输入参数：WM，WMM，B; 输出数据：R
    if PE（T）<0
        R=0; % 流域总径流量，mm
    else
        if  W>=WM
            A=WMM;  % 流域起始蓄水量为W0时，所对应的最大点蓄水容量mm或为对应的前期影响雨量[9343]p84;
        else
            A=WMM*（1-（1-W/WM）.^（1/（1+B）））;  % [3664]p145：式（5-9）
        end
        if A+PE（T）>0
            if A+PE（T）                R=PE（T）-WM+W+WM.*（1-（PE（T）+A）./WMM）.^（1+B）;  % [3664]p145：式（5-11）
            else
                R=PE（T）+W-WM;                               % [3664]p145：式（5-13）
            end
        else
            R=0;
        end
    end
    %===========================================
    % 以下为蒸发计算;;输入数据：PE，W，R； 输入参数：WM，WMM，B; 输出数据：EW
    if PE（T）<0   % 当降雨量小于蒸发量时时无产流R
        if WU+PE（T）>0
            EU=K*EM（T）;   % EU等于蒸散发能力，即在充分供水的条件下，流域蒸散发达到最大;[3664]p143：式（1）
            ED=0;
            EL=0;
            WU=WU+PE（T）;
        else              % [3664]p143：式（2）
            EU=WU+P（T）;
            WU=0;
            if WL>C*WLM   % [3664]p143：式（3）;C-深层蒸散发扩散系数;此式的物理意义是如果下层