MATLAB实现基于遗传算法/引力搜索算法优化新安江水文模型

1 新安江模型

1.1 新安江模型结构

新安江模型是分散性模型，可用于湿润地区与半湿润地区的湿润季节。当流域面积较小时，新安江模型采用集总模型，当面积较大时，采用分块模型。它把全流域分为许多块单元流域，对每个单元流域作产汇流计算，得出单元流域的出口流量过程。再进行出口以下的河道洪水演算，求得流域出口的流量过程。把每个单元流域的出流过程相加，就求得了流域的总出流过程。
该模型按照三层蒸散发模式计算流域蒸散发，按蓄满产流概念计算降雨产生的总径流量，采用流域蓄水曲线考虑下垫面不均匀对产流面积变化的影响。在径流成分划分方面，对三水源情况，按“山坡水文学”产流理论用一个具有有限容积和测孔、底孔的自由水蓄水库把总径流划分成饱和地面径流、壤中水径流和地下水径流。在汇流计算方面，单元面积的地面径流汇流一般采用单位线法，壤中水径流和地下水径流的汇流则采用线性水库法。
概念性模型的结构应该反映客观水文规律，参数应该代表流域的水文特征，把模型设计成为分散性的，主要是为了考虑降雨分布不均的影响，其次也便于考虑下垫面条件的不同及其变化。降雨分布不均，不但对汇流产生明显的影响，而且对产流也产生明显的影响。如果采用集总性模型，应用面平均雨量来进行计算，误差可能很大， 而且是系统性的。本次计算采用的降雨即为平均雨量。
新安江模型的流程图如图1所示。图中输入为实测降雨P和实测蒸散发能力，输出为流域出口断面流量Q和流域蒸散发量 E。方框内是状态变量，方框外是常数常量。模型主要由四部分组成，即蒸散发计算、产流量计算、水源划分和汇流计算。

新安江模型可以用于湿润和半湿润地区，当流域面积较小时，新安江模型采用集总模型，面积较大时，采用分单元模型。分单元模型把流域分成若干个单元面积，对每个单元用马斯京根汇流模型计算到达流域出口断面的流量过程。

1.2 模型参数种类及意义

2 新安江模型优化参数

新安江模型把全流域分成若干单元流域, 对每个单元流域分别作产汇流计算, 得出各个单元流域的出口流量过程, 再分别将出口以下的河道洪水演算至流域出口断面, 把同时刻的流量相加即求得流域出口的流量过程。
若忽略河网汇流过程，则需对三水源新安江水文模型的12个参数进行优化率定，可分为四个层次：

• 蒸散发参数: KC、WUM、WLM、C
• 产流量参数： WM、B
• 水源划分参数： SM、EX、 KSS、KG
• 坡地汇流参数：KKSS、KKG

2.1 蒸散发参数: KC、WUM、WLM、C

参数	解释
蒸散发能力折算系数KC	流域蒸散发能力与实测水面蒸发值之比。此参数控制着总水量平衡，因此，对水量计算是十分重要的。
上层蓄水容量WUM	包括植物截留量。在植被与土壤比较发育的流域，约为20mm； 在植被与土壤颇差的流域，约为5~6mm。
下层蓄水容量WLM	可取60~90mm。注意：因WDM=WM-WUM-WLM，故WDM不再是一个独立的参数。
深层蒸散发系数C	决定于深根植物占流域面积的比数，同时也与WUM+WLM值有关，此值越大，深层蒸散发越困难。一般经验，在江南湿润地区C值约为0.15-0.20左右，而在华北半湿润地区则在0.09-0.12左右。

2.2 产流量参数: WM、B

参数	解释
流域蓄水容量WM	流域干湿程度的指标。一般分为上层WUM、下层WLM和深层WDM，约为120~ 180mm。
蓄水容量曲线指数B	反映流域上蓄水容量分布的不均匀性。一般经验，流域越大，各种地质地形配置越多样，B值也越大。在山丘区域，很小面积(几平方公里)的B值为0.1左右，中等面积（300平方公里以内）的B值为0.2-0.3左右，较大面积（数千平方公里）的B值为0.3-0.4左右。

2.3 水源划分参数: SM、EX、KSS、KG

参数	解释
流域平均自由水蓄水容量SM	参数受降雨资料时段均化的影响，当以日为计算时段长时，一般流域的SM值约为10~50mm，当所选取的计算时段长较小时，SM要增大，这个参数对地面径流的多少起着决定性作用，因此十分重要。
自由水蓄水容量曲线指数EX	表示自由水容量分布不均匀性。通常EX取值在1~1.5之间。
自由水蓄水库对壤中流的出流系数KSS	一般取值为0.01~0.7。
自由水蓄水库对地下径流出流系数KG	一般取值为0.01~0.7。KSS和KG两个出流系数是并联的，其和代表着自由水出流的快慢。一般来说，KSS+KG=0.7，相当于从雨止到壤中流止的时间为3天。

2.4 坡地汇流参数：KKSS、KKG

参数	解释
壤中流的消退系数KKSS	如无深层壤中流时，此值趋于0；当壤中流很丰富时，此值趋于0.9，相当于汇流时间为10天。
地下水的日消退系数KKG	如以日为计算时段长，此值一般为0.98-0.998，相当于汇流时间为50-500天。

3 模型目标函数

目标函数用来评价实测流量与模拟流量过程的吻合程度, 不同的目标函数用来评价水文过程的不同的特征, 目标函数的选择对优选结果至关重要。一般可采用水量平衡、全流量过程均方误差、高水流量过程均方误差和枯水流量过程均方误差作为目标函数。

4 代码

主函数代码如下：

%% 清空环境变量
clc;
close all;
clear;
%% 程序与数据说明
%本程序利用新安江模型进行洪水预报演算，数据由源代码文件下目录提供
%E87-90.DAT.dat提供蒸发皿日蒸发量数据
%P87-90.DAT.dat提供雨量计日降雨量数据
%Q87-90.DAT.dat提供水文站日均流量数据
%F.DAT提供洪号、实测流量起涨时间、实测流量退水时间 
%% 将原始的流域统计资料进行赋值存储
Pfile = fopen('P87-90.DAT','rt');  %降雨
Pcell = textscan(Pfile,'%s %f');
fclose(Pfile);
Qfile = fopen('Q87-90.DAT','rt');  %流量
Qcell = textscan(Qfile,'%s %f');
fclose(Qfile);
Efile = fopen('E87-90.DAT','rt');  %蒸发
Ecell = textscan(Efile,'%s %f');
fclose(Efile);
DATE = Pcell{1}';
P = cell2mat(Pcell(:,2));
Q = cell2mat(Qcell(:,2));
E0 = cell2mat(Ecell(:,2));        %蒸发皿测得多年平均蒸发量
FDATE_FILE = fopen('F.DAT','rt');
FDATE_cell = textscan(FDATE_FILE,'%d %s %s');
fclose(FDATE_FILE);

FNUM=FDATE_cell{1};             %FNUM    洪号
FDATES=FDATE_cell{2};           %FDATES  实测流量起涨时间
FDATEE=FDATE_cell{3};           %FDATEE  实测流量退水时间
%% 将给定的参数进行赋值存储
WM = 126;            %流域平均蓄水容量
WUM = 63;            %流域上层土壤含水容量
WLM = 13;            %流域下层土壤含水容量
WDM = 50;            %流域深层土壤含水容量
KC = 0.71;              %流域蒸发折算系数
C = 0.17;               %蒸发扩散系数
B = 3;                     %流域蓄水容量曲线指数
IMP = 0.001;          %流域不透水面积比重
FE = 0.8;                %初始土壤含水容量折算系数 
SM = 36;              %流域平均自由水蓄水容量
EX = 0.46;             %自由水蓄水容量曲线的指数
KG = 0.05;             %自由水对地下水的日出流系数
K = 48;                  %蓄泄系数，单位h(小时)
KKG = 0.995;         %地下水消退系数
KSS = 0.06;           %自由水对壤中流的日出流系数
KKSS = 0.83;         %壤中流消退系数
F = 6448;              %流域面积单位为km^2     
 
%将给定的单位线进行赋值存储
load('q.mat')      %24h10min单位线
UH=q;              %单位线
%初始数据处理
EP=KC.*E0;          %水面蒸发乘系数KC转换成流域蒸发
S_0=SM*FE;          %初始自由水蓄量
WU_0=FE*WUM;
WL_0=FE*WLM;
WD_0=FE*WDM;        %初始土壤含水量
%% 函数调用及计算

%场次洪水统计
[PO,RCE,DATE_CE,QCE,RCE1] = Datastatistics(DATE,P,Q,F,K,FNUM,FDATES,FDATEE);

% 模型运算
[R,E,PE,W,WMM,WU,WL,WD,a] = runoffgenerate(P,EP,WM,WUM,WLM,WDM,WU_0,WL_0,WD_0,C,B,IMP);
[RS,RI,RG,FR,AU,S] = three_water_sources_division(P,EP,PE,W,WM,B,EX,KG,KSS,S_0,SM,R,IMP);
[QRS,QRI,QRG,Q_all,R_all] = confluence(UH,F,RS,RI,RG,KKG,KKSS);

% 预报评价
[Rmo,ERabs,ERrel,acceptrateR ]= EVALUATE_R(DATE,FDATES,FDATEE,R_all,RCE);
[Qmo,EQabs,EQrel,acceptrateQ,Date_savestart,Date_saveend,Date_save,Q_save,Qmo_save,P_save] = EVALUATE_Q(DATE,FDATES,FDATEE,Q_all,QCE,P,Q);
[DATEmo,EDabs,acceptrateD,LIM]=EVALUATE_D(DATE,FDATES,FDATEE,Q_all,DATE_CE);
[DC,DC_all ]=EVALUATE_DC(DATE,FDATES,FDATEE,Q_all,Q);

数据及完整代码请关注作者并私聊获取

5 案例

5.1 场次洪水资料

向家坪水文站控制集水面积6448km2，流域内有太平站、沙沟站、东江口、黄金站等28 个雨量站，各站均有较长期的实测雨量资料（1980-1990）。在洪水预报方案的研制中，已经将1980-1986 年资料作为模型率定期，选出了25 场暴雨与旬河同期实测流量资料作为优选产流参数和汇流过程的依据。
现需要将预留的1987-1990 年资料进行整理，为模型产汇流方案检验提供场次洪水资料。资料的整理和数据文件的组织都已通过编制的软件由计算机完成，1987-1990 年旬河流域面平均时段雨量数据，由数据文件P87-90.DAT 给出（设面平均雨量可采用算术平均法计算）；旬河1987-1990 年流量过程由数据文件Q87-90.DAT 提供。

5.2 模型结果

模拟流量与实测流量过程对比图：

参考