CN106682355A - 一种基于pso‑ga混合算法的水文模型参数率定方法 - Google Patents

Abstract

本发明提出的一种基于PSO‑GA混合算法的水文模型参数率定方法，属于水文模型参数率定领域。该方法首先在水文模型的初始化阶段输入降雨量，水文模型所有参数以及每个参数的最大值和最小值；然后执行水文模型算法程序，得到输出的预报流量值；根据计算出的预报流量值与实际值进行校验比对，得到确定性系数并进行判定；当确定性系数大于等于0.2时，说明预报流量值与实际值之间的误差大于等于20％，则通过GA‑PSO混合算法进行参数率定。该方法结合了GA和PSO两种算法各自寻优特点，在避免早熟收敛，提高搜索精确度。

Description

一种基于PSO-GA混合算法的水文模型参数率定方法

技术领域

本发明属于水文模型中参数率定领域，具体来说是一种结合了PSO粒子群算法及GA遗传算法的水文模型参数率定方法。

背景技术

自然界的水文现象是一种非常复杂的现象，它受到降雨特性、流域下垫面、人类活动等诸多因素的影响。在难以弄清水文现象的规律之前，通过建立模型对水文过程进行模拟(试验)是一种行之有效的途径，该模型即水文模型。水文模型参数可分为两类：一类参数具有明确的物理含义，可以根据实际情况进行确定，比如不透水面积占比，另一类是没有物理含义或者物理含义不明确的参数，地下水消退系数，壤中流日出流系数，这些参数需要根据以往的观测数据进行率定。第二类模型参数往往表现出高维，多峰值，非线性，不连续，非凸性及带噪声等复杂特征。

然而，受当前计算能力和模型设计水平的限制，传统水文模型参数率定仍然以人工试错法为主。人工试错法是通过比较模拟值与实测值拟合度，人工调整参数的一种参数率定方法，这种方法过多依赖个人经验主观性较强，率定效果欠佳，不利于模型的推广应用。因此，近年发展起来的自动参数率定方法逐渐受到人们重视，并且开始成为一种发展方向。自动参数率定是由计算机按照一定的规则自动优选参数。在传统的自动参数率定方法中有最速下降法，步长加速法，模式法等，但这些算法可能只对某一类方法奏效，且对模型要求较为苛刻，比如要求模型具有连续，可导，单峰等特性；而一些新兴优化算法就在此显示出优势来，在这之中，GA遗传算法是一种模拟遗传学激励的生物进化的计算模型，通过模拟自然进化过程搜索最优解的方法，具有较强的全局搜索能力,但局部搜索能力较弱,且变异具有随机性并在选择操作时会损失一部分适应度较小的粒子，而这些粒子身上可能也有好的经验可以学习，由于缺乏有效的局部区域搜索机制，GA算法在接近最优解时收敛缓慢甚至出现停止现象，另一方面，PSO粒子群算法是受飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型，在算法早期，也存在精度低，易发散的缺点，而在收敛的情况下，由于所有粒子都朝着最优的方向飞去，所以粒子趋向于同一化，使得后期收敛速度明显变慢，造成了容易陷入局部最优的缺陷。

发明内容

本发明的目的是针对现有技术的不足之处，提出一种基于PSO-GA混合算法的水文模型参数率定方法。该方法能够充分发挥两种算法各自寻优特点，避免早熟收敛，提高搜索精确度，并有效解决了参数自动率定中的模型要求苛刻、限制多的问题，具有较强的理论依据，可以在较短的时间内快速稳定找到全局最优或近似最优解。

本发明提出的一种基于PSO-GA混合算法的水文模型参数率定方法，其特征在于，包括以下步骤：

1)选择水文模型，在水文模型的初始化阶段输入降雨量，水文模型所有参数以及所述水文模型每个参数的最大值和最小值；

2)执行水文模型算法程序，得到输出的预报流量值；

3)根据步骤2)中计算出的预报流量值与实际值进行校验比对，遵循水利部的水文情报预报规范的标准，采取确定性系数计算方法如式(1)所示：

式中，DC为确定性系数,y为实际值，yc为预报流量值，n为序列长度即降雨量的次数，i为当前序列值；

4)对步骤3)得到的确定性系数DC的值进行判定：根据步骤3)中的标准，若DC<0.2，说明预报流量值与实际值之间的误差小于20％，误差在标准范围内，则无需参数率定，本次率定流程终止；否则，若DC≧0.2，说明预报流量值与实际值之间的误差大于等于20％，则转入步骤5)，进行参数率定；

5)通过GA-PSO混合算法对水文模型参数进行率定，之后重新返回步骤2)。

本发明的特点及有益效果：

本发明提出的一种基于PSO-GA混合算法的水文模型参数率定方法，所采用的混合算法适用性强，可靠性高，其主要步骤均由计算机完成，利用了GA算法的随机性以增加搜索范围，然后利用PSO算法在找到的最优微粒附近进行更细致的搜索，并结合水文模型算法及其参数，初始值及范围从而得到最优解；该混合算法能够充分发挥两种算法各自寻优特点，避免早熟收敛，提高搜索精确度。本发明有效解决了传统人工参数率定的费时费力和早前方法的精度低，效果差等问题。达到了快速高效参数率定的效果。

附图说明

图1为本发明方法的流程图。

图2为本发明实施例的新安江水文模型算法流程图。

图3为本发明方法中的GA-PSO混合算法流程图。

具体实施方式

本发明提出的一种基于PSO-GA混合算法的水文模型参数率定方法，下面结合附图和具体实施例进一步详细说明如下。

本发明实施例以三水源新安江流域水文模型为例对本发明提出的一种基于PSO-GA混合算法的参数率定方法进一步详细说明。

本发明提出的一种基于PSO-GA混合算法的参数率定方法，总体流程图如图1所示，包括以下步骤：

1)选择水文模型，本实施例选择的是新安江水文模型；在新安江水文模型的初始化阶段输入降雨量，水文模型所有参数包括：蒸散发系数K,地下水出流系数KG,地下水消退系数CG等，及所述水文模型每个参数的最大值和最小值。

2)执行水文模型算法程序，包括蒸散发，水源划分，产汇流等步骤，最终得到输出的预报流量值。本实施例执行的是新安江水文模型算法程序，所涉及的新安江水文模型根据河海大学建立的国内第一个完整的流域水文模型思想，将其结合JAVA编程语言进行编写得到新安江水文模型算法程序；在实际操作时，可由本领域内编程人员通过计算机编程实现。

3)根据步骤2)中计算出的预报流量值与实际值进行校验比对，遵循水利部的水文情报预报规范的标准，采取确定性系数计算方法如式(1)所示：

式中，DC为确定性系数,y为实际值，yc为预报流量值，n为序列长度即降雨量的次数，i为当前序列值。

4)对步骤3)得到的确定性系数DC的值进行判定：根据步骤3)中的标准，若DC<0.2，说明预报流量值与实际值之间的误差小于20％，误差在标准范围内，可以接受，则无需参数率定，本次率定流程终止；否则，若DC≧0.2，说明预报流量值与实际值之间的误差大于等于20％，则转入步骤5)，进行参数率定。

5)通过GA-PSO混合算法对水文模型参数进行率定，之后重新返回步骤2)。

上述步骤2)中执行新安江水文模型算法程序，具体实现流程如图2所示，首先进行一次降雨，按照降雨所在的区域不同分为透水面积和不透水面积，透水面积进入土壤并进行蒸散发及相应的产流，最终将产流汇流为总入流，具体包括以下步骤：

2-1)在水文模型中输入降雨量及蒸散发能力等值(即上层张力水量WU,下层张力水量WL,深层张力水量W-(WU+WL),蒸散发能力折算系数K,张力水蓄水容量曲线指数B,不透水面积占比IMP,自由水蓄水量S,自由水蓄水容量曲线指数EX,自由水蓄水水库地面径流RS,地下水径流RG及日出流系数KG,壤中流径流RI及日出流系数KI,地面径流消退系数CS，地下水消退系数CG,壤中流消退系数CI，时段步长DT,单位线UH,降雨时间序列P,蒸散发能力EP,单位出流量Q,初始产流面积占比FR)，其中待率定的为：蒸散发能力折算系数K、自由水蓄水量S、地下水日出流系数KG、壤中流日出流系数KI、地面径流消退系数CS、地下水消退系数CG、壤中流消退系数CI 6个参数；

2-2)将一场降雨的降雨量下分为透水面积与不透水面积，透水面积水量下渗入土壤等，不透水面积则直接产流；透水面积形成产流流域，产流的流量为R(1-IMP)，该产流流量进入土壤，式中R为单元面积，IMP为不透水占比面积；

2-3)张力水进行蒸散发，自由水进行水源划分及产流。将步骤2-2)得到的产流流量进入到上层张力水中，张力水更新得到新的上层张力水量WU＝WU+R*(1-IMP),然后土壤中负责蒸散发部分的张力水按照上下深层的系数进行蒸散发：

若P+WU>＝EP,则EU＝EP,EL＝0，ED＝0；若P+WU<EP,则EU＝P+WU；

若WL>C*LM,则EL＝(EP-EU)*WL/LM,ED＝0；若WL<C*LM且WL>＝C*(EP-EU),则

EL＝C*(EP-EU),ED＝0；若WL<C*LM且WL<C*(EP-EU),则EL＝WL,ED＝C*(EP-EU)-WL；

C为深层蒸散发扩散系数，LM为下层张力水容量，EL为下层蒸散发量，EU为上层蒸散发量，ED为深层蒸散发量，PE为净雨量；

2-4)土壤中自由水进行水源划分，若净雨量PE小于0，则自由水S按S＝S-(RI_t+RG_t)/FR进行更新,按照步骤2-1)的待率定参数，即地面径流消退系数CS，壤中流消退系数CI和地下水消退系数CG分别按公式

QS_t＝CS*QS_t-1+(1-CS)*RS_t*U,QI_t＝CS*QS_t-1+(1-CI)*RI_t*U,QG_t＝CS*QS_t-1+(1-CG)*RG_t*U,产生本次地面径流量QS，壤中流流量QI与地下径流量QG，其中QS_t-1，QI_t-1，QG_t-1分别表示上次的地面径流量，上次的壤中流流量与上次的地下径流量，RS_t,RI_t,RG_t分别表示本次的地面径流，本次的壤中流，本次的地下径流,U为单位转换系数；

2-5)将步骤2-4)得到的地面径流量QS，壤中流流量QI和地下径流流量QG进行汇流得到总流量，即预报流量值。之后通过计算由步骤21)传入的待率定参数而得到的预报流量值与实测值进行校验比对。

上述步骤5)中通过GA-PSO混合算法对水文模型参数进行率定，该混合算法可由计算机完成，流程如图3所示，具体包括以下步骤：

5-1)初始化两个种群P1，P2，其中P1种群用于PSO粒子群算法，P2种群用于GA遗传算法。其中，种群中每个“个体”都是一个一维向量，向量中的每个元素称为“基因”，本发明中“个体”即为所有待辨识的水文模型参数，“基因”则为每个待辨识的水文模型参数。

设置最大迭代次数(最大迭代次数可取数千至数万，并无严格限制，具体取值需结合具体情况，本实施例中取5000)以作为判断终止的条件之一，初始化代数T为0。

5-2)计算P2种群的适应度。针对于水文模型，构造适应度函数计算公式如式(2)所示：

式中，f表示适应度，Q_obs，i为实际值，Q_sim，i为计算出的预报流量值。i表示初始代数，n表示当前代数。根据构造的适应度函数，分别计算P2种群中每个粒子的适应度值。适应度值反映的是该粒子实际经过路径与目标路径的距离即预报流量值距离实际值的差值。

5-3)根据步骤5-2)得到得P2中每个粒子的适应度值判断终止条件，即是否已经找到了最优解(最优解即最符合期望的值，在本实施例适应度函数式(2)中，最优解可设置为f＝10000)或已经达到了迭代的最大次数。若判断为是，则转步骤5-9)，否则执行步骤54)。

5-4)对P2中每个粒子按照适应度值从大到小排序进行划分种群；筛选出种群P2中粒子基因较好的个体，即f越大越好(如可取999或9999)的粒子，将排序后筛选出的粒子以φ比例(φ取0-1之间的实数，本实施例中取0.2)选取P2中粒子加入到种群P1，得到新一代种群P1’。

5-5)根据步骤5-4)得到新一代种群P1’，采用PSO粒子群算法，根据如式(3)和(4)所示的粒子群寻优方法得到自身和全局最优粒子，更新种群P1，得到新的种群P1”。根据式(2)计算更新后的种群P1”所有粒子的适应度函数，得到P1”的每个粒子适应度值。上述所采用的粒子群寻优方法具体如式(3)(4)所示：

V_id＝wV_id+c₁r₁(P_id-X_id)+c₂r₂(P_gd-X_id) (3)

X_id＝X_id+v_id (4)

式中，c1和c2分别是非负数学习因子，r1和r2是介于0-1的两个随机数，v_id和x_id分别代表当前的速度和位置。从上式可知，速度的更新分为三个部分：历史速度v_id，粒子自身在运动中的最优组合，粒子在全局中运动的最优组合；w为惯性权重因子，该惯性权重用于平衡全局和局部搜索能力，较大的惯性权重更倾向于全局搜索，而较小的惯性权重适于局部搜索。比较各个粒子的当前适应值与自身历史中的最优适应值,如果当前适应值优于自身历史最优适应值,则置当前适应值为自身历史最优适应值，该粒子为自身最优粒子。比较各个粒子的当前适应值与该种群的全局中的最优适应值,如果某粒子的当前适应值优于该种群的全局最优适应值,则置该粒子当前适应值为该种群的全局最优适应值，该粒子为全局最优粒子。

5-6)根据步骤5-5)得到的P1”中每个粒子的适应度值判断终止条件(最优解即最符合期望的值，在本实施例适应度函数式(2)中，最优解为f＝10000)，即是否已经找到了最优解或已经达到了迭代的最大次数。若判断为是，转步骤5-9)，否则执行步骤5-7)。

5-7)按照φ比例将P1”产生的随机粒子补充到种群P2，并除去种群P2中在步骤5-4)迁移到P1中参加粒子群算法的粒子，形成种群P2’；对种群P2’采用GA遗传算法寻优，按照选择，交叉，变异操作，得到新一代的种群P2”。

其中选择方法采用轮盘赌(Roulette Wheel Selection)选择法，则种群P2’中粒子被选中的概率如式(5)所示：

式中，i表示初始代数，n表示当前代数，f_i表示各染色体的适应度值(同之前的每个粒子的适应度值)，可见适应度越大的个体被选中的概率越大；而交叉即在把个体进行二进制编码后随机把其中几个位于同一位置的编码进行交换，产生新的个体；变异则按照一定突变的概率将之前得到的二进制编码的部分进行取反，基因串上的“0”或“1”有一定几率变成与之相反的“1”或“0”。

5-8)迭代次数加1，T＝T+1；之后重新返回步骤5-2)。

5-9)输出最终的最优解，即为率定得到的参数。

Claims (3)

1.一种基于PSO-GA混合算法的水文模型参数率定方法，其特征在于，包括以下步骤：

1)选择水文模型，在水文模型的初始化阶段输入降雨量，水文模型所有参数以及所述水文模型每个参数的最大值和最小值；

2)执行水文模型算法程序，得到输出的预报流量值；

3)根据步骤2)中计算出的预报流量值与实际值进行校验比对，遵循水利部的水文情报预报规范的标准，采取确定性系数计算方法如式(1)所示：

$DC=1-\frac{\underset{i-1}{\overset{n}{\&Sigma;}}{\&lsqb;yc\left(i\right)-y\left(i\right)\&rsqb;}^2}{\underset{i-1}{\overset{n}{\&Sigma;}}{\&lsqb;yc\left(i\right)-\stackrel{\&OverBar;}{y}\left(i\right)\&rsqb;}^2}---\left(1\right)$

式中，DC为确定性系数,y为实际值，yc为预报流量值，n为序列长度即降雨量的次数，i为当前序列值；

4)对步骤3)得到的确定性系数DC的值进行判定：根据步骤3)中的标准，若DC<0.2，说明预报流量值与实际值之间的误差小于20％，误差在标准范围内，则无需参数率定，本次率定流程终止；否则，若DC≧0.2，说明预报流量值与实际值之间的误差大于等于20％，则转入步骤5)，进行参数率定；

5)通过GA-PSO混合算法对水文模型参数进行率定，之后重新返回步骤2)。

2.如权利要求1所述的方法，其特征在于，所述步骤2)中执行水文模型算法程序，具体包括以下步骤：

2-1)在水文模型中输入降雨量及蒸散发能力的参数值，具体包括：上层张力水量WU,下层张力水量WL,深层张力水量W-(WU+WL),蒸散发能力折算系数K,张力水蓄水容量曲线指数B,不透水面积占比IMP,自由水蓄水量S,自由水蓄水容量曲线指数EX,自由水蓄水水库地面径流RS,地下水径流RG及日出流系数KG,壤中流径流RI及日出流系数KI,地面径流消退系数CS，地下水消退系数CG,壤中流消退系数CI，时段步长DT,单位线UH,降雨时间序列P,蒸散发能力EP,单位出流量Q,初始产流面积占比FR；其中待率定的参数为：蒸散发能力折算系数K、自由水蓄水量S、地下水日出流系数KG、壤中流日出流系数KI、地面径流消退系数CS、地下水消退系数CG、壤中流消退系数CI 6个参数；

2-2)将一场降雨的降雨量下分为透水面积与不透水面积，透水面积水量下渗入土壤，不透水面积则直接产流；透水面积形成产流流域，产流的流量为R(1-IMP)，该产流流量进入土壤，式中R为单元面积，IMP为不透水占比面积；

2-3)张力水进行蒸散发，自由水进行水源划分及产流；将步骤2-2)得到的产流流量进入到上层张力水中，张力水更新得到新的上层张力水量WU＝WU+R*(1-IMP),然后土壤中负责蒸散发部分的张力水按照上下深层的系数进行蒸散发：

若P+WU>＝EP,则EU＝EP,EL＝0，ED＝0；若P+WU<EP,则EU＝P+WU；

若WL>C*LM,则EL＝(EP-EU)*WL/LM,ED＝0；若WL<C*LM且WL>＝C*(EP-EU),则EL＝C*(EP-EU),ED＝0；若WL<C*LM且WL<C*(EP-EU),则EL＝WL,ED＝C*(EP-EU)-WL；

C为深层蒸散发扩散系数，LM为下层张力水容量，EL为下层蒸散发量，EU为上层蒸散发量，ED为深层蒸散发量，PE为净雨量；

2-4)土壤中自由水进行水源划分，若净雨量PE小于0，则自由水S按S＝S-(RI_t+RG_t)/FR进行更新,按照步骤2-1)的待率定参数，即地面径流消退系数CS，壤中流消退系数CI和地下水消退系数CG分别按公式：

QS_t＝CS*QS_t-1+(1-CS)*RS_t*U,QI_t＝CS*QS_t-1+(1-CI)*RI_t*U,

QG_t＝CS*QS_t-1+(1-CG)*RG_t*U,

产生本次地面径流量QS，壤中流流量QI与地下径流量QG；其中，QS_t-1，QI_t-1，QG_t-1分别表示上次的地面径流量，上次的壤中流流量与上次的地下径流量，RS_t,RI_t,RG_t分别表示本次的地面径流，本次的壤中流，本次的地下径流,U为单位转换系数；

2-5)将步骤2-4)得到的地面径流量QS，壤中流流量QI和地下径流流量QG进行汇流得到总流量，即预报流量值。

3.如权利要求1所述的方法，其特征在于，所述步骤5)中通过GA-PSO混合算法对水文模型参数进行率定，具体包括以下步骤：

5-1)初始化两个种群P1，P2，其中P1种群用于PSO粒子群算法，P2种群用于GA遗传算法；设置最大迭代次数以作为判断终止的条件之一，初始化代数T为0；

5-2)计算P2种群的适应度；针对于水文模型，构造适应度函数计算公式如式(2)所示：

$f=\underset{i=1}{\overset{n}{\&Sigma;}}{Q}_{obs,i}/\underset{i=1}{\overset{n}{\&Sigma;}}{Q}_{obs,i}-Q_{sim,i}---\left(2\right)$

式中，f表示适应度，Q_obs，i为实际值，Q_sim，i为计算出的预报流量值，i表示初始代数，n表示当前代数；根据构造的适应度函数，分别计算P2种群中每个粒子的适应度值，适应度值反映的是该粒子实际经过路径与目标路径的距离即预报流量值距离实际值的差值；

5-3)根据步骤5-2)得到得P2中每个粒子的适应度值判断终止条件，即是否已经找到了最优解或已经达到了迭代的最大次数；若判断为是，则转步骤5-9)，否则执行步骤5-4)；

5-4)对P2中每个粒子按照适应度值从大到小排序进行划分种群；筛选出种群P2中粒子基因较好即f值较大的个体；将排序后筛选出的的粒子以φ比例选取P2中粒子加入到种群P1，得到新一代种群P1’，φ为0-1之间的实数；

5-5)根据步骤5-4)得到新一代种群P1’，采用PSO粒子群算法，根据如式(3)和(4)所示的粒子群寻优方法得到自身和全局最优粒子，更新种群P1，得到新的种群P1”；根据式(2)计算更新后的种群P1”所有粒子的适应度函数，得到P1”的每个粒子适应度值；上述所采用的粒子群寻优方法具体如式(3)(4)所示：

V_id＝wV_id+c₁r₁(P_id-X_id)+c₂r₂(P_gd-X_id) (3)

X_id＝X_id+v_id (4)

式中，c1和c2分别是非负数学习因子，r1和r2是介于0-1的两个随机数，v_id和x_id分别代表当前的速度和位置；w为惯性权重因子；

比较各个粒子的当前适应值与自身历史中的最优适应值,如果当前适应值优于自身历史最优适应值,则置当前适应值为自身历史最优适应值，该粒子为自身最优粒子；比较各个粒子的当前适应值与该种群的全局中的最优适应值,如果某粒子的当前适应值优于该种群的全局最优适应值,则置该粒子当前适应值为该种群的全局最优适应值，该粒子为全局最优粒子；

5-6)根据步骤5-5)得到的P1”种群中每个粒子的适应度值判断终止条件，即是否已经找到了最优解或已经达到了迭代的最大次数；若判断为是，则转步骤5-9)，否则执行步骤5-7)；

5-7)按照φ比例将种群P1”产生的随机粒子补充到种群P2，并除去步骤5-4)中迁移到P1中参加粒子群算法的粒子，形成种群P2’；对种群P2’采用GA遗传算法寻优，按照选择，交叉，变异操作，得到新一代的种群P2”；

5-8)迭代次数加1，T＝T+1，之后重新返回步骤5-2)；

5-9)输出最终的最优解，即为率定得到的参数。