河道洪水演算是根据河道上游水情推求下游水情的一种洪水预报方法, 在防洪规划和水库调洪计算中有着重要的作用. 马斯京根模型是自然河道和河流中最常用的洪水演算方法之一, 它是美国著名科学家McCarthy GT在1938年研究了马斯京根河流域后首次提出的^[1]. 采用马斯京根模型进行洪水演算时, 由于参数值不能从历史入流和出流水文图得到, 所以模型的参数估计问题成为最重要的任务. 求解马斯京根模型参数的传统方法有试错法^[2]、最小二乘法^[3]和非线性规划法^[4]等, 这些方法在实际应用中计算过程复杂, 不具有足够的稳定性和高效性. 因此, 近年来许多学者将智能优化算法应用于马斯京根模型参数估计, 并取得了较好的优化效果, 如遗传算法^[5]、差分进化算法(DE)^[6]、入侵杂草优化算法^[7]、飞蛾火焰算法^[8]和免疫克隆选择算法^[9]等. 然而, 由于马斯京根模型本身的近似性及这些算法存在早熟收敛和易陷入局部极值等问题, 使得洪水演算的精度并不十分理想.

受自然界中灰狼的社会等级制度和捕食行为的启发, 学者Mirjalili S等人^[10]在2014年提出了一种新型元启发式算法, 称为灰狼优化算法 (Grey Wolf Optimizer, GWO). 算法的优点主要是结构简单、需要设置的参数少和在实验编码中容易实现等, 且研究结果表明其优化性能优于DE算法、PSO算法和引力搜索算法. 然而, 基本GWO算法存在依赖初始种群、易陷入局部最优和难以协调其探索和开发能力等缺点, 针对这些缺点, 本文对基本GWO算法进行了改进, 利用混沌Chebyshev映射产生初始灰狼种群, 增强全局搜索过程中的种群多样性; 给出一种随机分布调整收敛因子策略, 以平衡算法的全局搜索和局部搜索能力. 采用7个基准测试函数对比仿真验证了改进GWO算法(IGWO)的有效性. 最后用IGWO算法来估计马斯京根模型中的参数, 对1961年称沟湾至临清段洪水进行演算, 与现有的其它参数估计方法进行了对比分析.

1 马斯京根模型理论

马斯京根模型是一种常用的河道洪水演算方法,该方法利用河段水量槽蓄方程代替复杂的水动力方程, 极大地简化了计算. 其基本方程为^[9]:

$\left\{ {\begin{aligned} & {\displaystyle\frac{{{\rm{d}}W}}{{{\rm{d}}t}} = I - Q} \\ & {W = K\left[ {xI + \left( {1 - x} \right)Q} \right]} \end{aligned}} \right.$

(1)

式(1)中, W为河段的槽蓄量, t为时间, I和Q分别为河段的入流量和出流量, x为流量比重因子, K为槽蓄系数.

式(1)的差分解为

$\left\{ {\begin{array}{*{20}{l}} {\tilde Q\left( 1 \right) = Q\left( 1 \right)} \\ {\tilde Q\left( i \right) = {C_0}I\left( i \right) + {C_1}I\left( {i - 1} \right) + {C_2}\tilde Q\left( {i - 1} \right),i = 2,3,\cdots,n} \end{array}} \right.$

(2)

式(2)中, $\tilde Q\left( i \right)$ 和 $Q\left( i \right)$ 分别为i时刻下断面的演算出流量和实际出流量, $I\left( i \right)$ 为第i个演算时段的入流量, n为演算时段个数, ${C_0}$ 、 ${C_1}$ 、 ${C_2}$ 为流量演算系数, 满足

${C_0} + {C_1} + {C_2} = 1$

(3)

根据式(2)、(3)建立马斯京根模型的参数优化模型, 优化模型如下:

$\begin{split}& \min F = {\displaystyle\sum\limits_{i = 2}^n {\left( {\left| {{C_0}I\left( i \right) + {C_1}I\left( {i - 1} \right) + {C_2}\tilde Q\left( {i - 1} \right) - Q\left( i \right)} \right|} \right)} ^q}\\& {\rm{s}}.{\rm{t}}.\;\;\;\;\;\;{g_1}:{C_0} \in \left[ {0,1} \right]\\& \;\;\;\;\;\;\;\;\;\;\;{g_2}:{C_1} \in \left[ {0,1} \right]\\& \;\;\;\;\;\;\;\;\;\;\;{g_3}:1 - {C_0} - {C_1} \in \left[ {0,1} \right]\end{split}$

(4)

文中取常数 $q = 1$ , 此时为最小一乘准则, 对约束条件 ${g_3}$ 采用惩罚函数进行处理, 优化准则函数取为:

$\min f = F + p\left( {{g_3}\left( {{C_0},{C_1}} \right)} \right)$

(5)

式(5)中 $p\left( {{g_3}\left( {{C_0},{C_1}} \right)} \right)$ 为惩罚项, 当 ${g_3}$ 满足时取0, 否则取 ${10^8}$ .

2 灰狼优化算法 2.1 基本灰狼优化算法

GWO算法^[10]是通过模拟自然界中灰狼群体的社会等级机制和捕食行为而提出的一种新型群体智能优化算法. 灰狼种群具有严格的等级制度, 如图1所示.

图1中, α是最高级别的灰狼, β是α的下属灰狼, δ服从α及β的命令, 等级最低的灰狼称为ω. 灰狼群体严格的社会等级机制在实现群体高效捕杀猎物的过程中发挥着至关重要的作用. 在捕食猎物时, 群体中其它灰狼在头狼α 的带领下有组织地对猎物进行围攻. 首先, 狼群通过气味等信息追踪、靠近猎物; 然后, 狼群包围猎物, 直到猎物停止移动; 最后, 攻击猎物.

假设灰狼种群数目为N, 灰狼的位置为X. 群体最优解为α, 次优解为β, 第三最优解为δ, 其它个体为ω. 则灰狼捕食行为的数学模型描述如下:

${{D}} = \left| {{{C}} \cdot {{{X}}_p}\left( t \right) - {{X}}\left( t \right)} \right|$

(6)

${{X}}\left( {t + 1} \right) = {{{X}}_p}\left( t \right) - {{A}} \cdot {{D}}$

(7)

其中, t表示当前迭代次数, A和C是系数向量, ${{{X}}_p}$ 是猎物的位置向量, X是灰狼的位置向量.

${{A}} = 2a \cdot {{{r}}_1} - a$

(8)

${{C}} = 2 \cdot {{{r}}_2}$

(9)

其中, ${r_1}$ 和 ${r_2}$ 均是[0, 1]之间的随机数, $a$ 是收敛因子, 随着迭代次数从2线性递减到0.

种群中其它灰狼的位置由α、β和δ的位置共同决定:

$\left\{ {\begin{array}{*{20}{c}}{{{{D}}_\alpha } = \left| {{{{C}}_1} \cdot {{{X}}_\alpha } - {{X}}} \right|,\;\;\;{{{X}}_1} = {{{X}}_\alpha } - {{{A}}_1} \cdot {{{D}}_\alpha }}\\{{{{D}}_\beta } = \left| {{{{C}}_2} \cdot {{{X}}_\beta } - {{X}}} \right|,\;\;\;{{{X}}_2} = {{{X}}_\beta } - {{{A}}_2} \cdot {{{D}}_\beta }}\\{{{{D}}_\delta } = \left| {{{{C}}_3} \cdot {{{X}}_\delta } - {{X}}} \right|,\;\;\;\;{{{X}}_3} = {{{X}}_\delta } - {{{A}}_3} \cdot {{{D}}_\delta }}\end{array}} \right.$

(10)

${{X}}\left( {t + 1} \right) = \frac{{{{{X}}_1} + {{{X}}_2} + {{{X}}_3}}}{3}$

(11)

2.2 改进的灰狼优化算法

GWO算法自提出以来已在许多领域得到了应用, 但其仍存在依赖初始种群和易陷入局部最优的缺点. 为了解决这些问题, 对基本GWO算法作如下改进:

(1) 混沌种群初始化

基本GWO算法在搜索空间中随机初始化灰狼种群, 容易造成灰狼位置分布不均匀, 导致种群多样性差. 混沌映射具有较好的遍历性和不重复性, 因此可以用来代替种群随机初始化. 本文采用Chebyshev混沌映射产生灰狼种群的初始位置, 映射方程^[11]为:

${x_{n + 1}} = \cos \left( {k{\rm{cos}^{ - 1}}\left( {{x_n}} \right)} \right)$

(12)

式(12)中k为阶次, 当 $k \ge 2$ 时无论初始值选择如何接近, 迭代出来的序列在此范围内都是混沌和遍历的.

(2) 随机分布调整收敛因子

由文献[10]可知, 当|A|>1时, 灰狼群体将扩大搜索范围寻找猎物, 即全局搜索; 当|A|<1时, 灰狼群体将收缩搜索范围对猎物进行攻击, 即局部搜索. 因此,A的大小与GWO算法的全局搜索和局部搜索能力有很大关系. 由公式(8)可以看出, A随着收敛因子 $a$ 的变化而变化.

若将收敛因子a设定为服从某种分布的随机数, 利用随机变量的特性调整控制参数, 既能在进化初期有机会取得较大或较小的值, 又能在进化后期取得较小或较大的值, 从而协调GWO算法的探索和开发能力. 基于上述分析, 受文献[12]的启发, 本文提出一种随机分布调整收敛因子策略:

$a = \left( {{a_{\max }} - {a_{\min }}} \right)rand() + \sigma randn()$

(13)

式(13)中, a_max、a_min是收敛因子的最大值和最小值, rand()是[0, 1]均匀分布的随机数, randn()是正态分布的随机数, σ是标准差.

2.3 仿真实验

为了验证IGWO算法的性能, 利用其对国际上通用的7个经典测试函数进行求解, 并与基本GWO算法、PSO算法和细菌觅食算法(BFA)进行对比, 比较结果如表1所示. 其中PSO算法是通过模拟自然界鸟群觅食行为来寻找最优解的一种演化算法. BFA算法是基于大肠杆菌觅食行为提出的一种仿生类群体智能算法. 同时给出了Sphere函数和Ackley函数的收敛曲线图, 如图2和图3所示.

算法的参数设置: 各算法的种群规模均为30, 最大迭代次数为500, PSO算法中粒子速度V_max=5, c₁=c₂=0.5; BFA算法中最大趋化代数N_c=20, 最大前进步数N_s=4, 复制操作执行次数Nre=10, 迁徙操作执行次数N_ed=2, 迁徙概率P_ed=0.25.

从表1的数据结果可以看出, 在种群规模和最大迭代次数设置相同的情况下, 与基本GWO算法、PSO算法和BFA算法相比, IGWO算法在7个测试函数上均表现出了较好的寻优性能, 收敛精度有较大的提高, 特别是Rastrigin 函数收敛到了理论最优值. 从图2和图3可以看出, IGWO算法具有更高的求解精度和更快的收敛速度.

下面验证Chebyshev混沌初始化策略对算法性能的影响. 图4和图5给出了Schwefel 1.2函数和Quartic函数在不同初始化策略下的收敛曲线, GWO表示随机初始化种群, Chebyshev-GWO表示Chebyshev混沌初始化. 从图4和图5可以看出, 本文所采用的Chebyshev映射初始化种群策略, 有效地提高了算法的收敛精度和收敛速度.

此外, 讨论文中给出的随机调整收敛因子策略对算法性能的影响. 图6和图7是Schwefel 2.22函数和Ackley函数采用不同收敛因子时的收敛曲线, GWO表示收敛因子线性递减的GWO算法, a-GWO表示随机调整收敛因子的GWO算法. 从图6和图7可以看出, 本文提出的随机调整收敛因子策略能有效地提高算法的收敛精度.

3 马斯京根模型的参数估计步骤

利用改进的灰狼优化算法(IGWO)估计马斯京根模型参数时, 每头灰狼的位置对应一组解, 以实际流量和演算流量的最小绝对误差和为目标函数, 运用IGWO对目标函数进行求解, 即以式(5)作为目标函数, 寻找f的最小值, 对应于最优准则函数值. 灰狼的个体位置为一个2维矢量 ${{X}} = \left( {{x_{i1}},{x_{i2}}} \right)$ , 其中 $i = 1,2, \cdots ,N$ 代表种群数目, 2代表所要估计的两个参数, 即C₀和C₁. 算法的具体执行步骤如下:

Step1. 设置种群规模N=30, 维数dim=2, 最大迭代次数t_max=500, 搜索区间为[0, 1], 随机生成 $a$ 、A、C等参数.

Step2. 利用Chebyshev混沌映射策略初始化灰狼种群;

Step3. 将式(5)设置为适应度函数, 计算种群中所有灰狼个体的适应度值, 选择前三个最好的狼, 记录其位置X_α、X_β和X_δ;

Step4. 利用公式(10)和(11)更新种群中其它灰狼个体的位置;

Step5. 利用公式(13)计算 $a$ , 然后利用公式(8)和(9)更新A, C的值;

Step6. 判断算法是否满足结束条件, 若达到预定的最大迭代次数t_max, 则停止计算, 输出最优位置 ${{{X}}_\alpha }$ 和最优函数值; 否则, 重复执行Step3–Step5.

4 应用实例

以南运河称沟湾至临清河段1961年8月的一次洪水过程的流量演算为例(该河段长83.8 km, 中间支流影响不明显), 数据来源于文献[9], 演算时段 $\Delta t$ 取12 h. 利用IGWO算法估计马斯京根模型的参数, 并将演算结果的统计指标与现有参数估计方法估计的结果进行比较. 洪水演算结果如表2所示. 图8更直观、更清晰的展现了表2中的演算结果. 各方法比较结果如表3所示. 图9给出了GWO算法和IGWO算法估计马斯京根模型参数时的迭代收敛图.

从表2和表3的计算结果和各项评价指标评价的对比结果可以看出, 与现有的参数估计方法相比, IGWO算法所求得的最优化准则函数值, 平均绝对误差和平均相对误差更小, 表明利用该算法估计马斯京根模型参数具有更高的求解精度. 从图9可以看出, IGWO算法估计参数时收敛速度比基本GWO更快.

5 结论

本文将改进的灰狼优化算法引入马斯京根模型, 对模型中的相关参数进行估计. 以南运河称沟湾至临清段1961年的洪水资料为例, 验证了IGWO算法求解马斯京根模型参数的有效性. 相比较于现有的其它参数估计方法, IGWO算法能很好地解决马斯京根模型的参数估计问题, 具有更高的求解精度和更好的实际应用价值.