基于改进鲸鱼优化算法的重力坝变形参数反演

高林钢，李同春,2，林潮宁，朱致远，郑 斌

(1.河海大学 水利水电学院，江苏 南京 210098; 2.水安全与水科学协同创新中心，江苏 南京 210098)

1 研究背景

随着水利资源的持续开发，越来越多的大坝投入建设与运行，建坝的选址条件也愈发复杂，而大坝的规模也越来越高。大坝的建设投资十分巨大，但创造的效益亦非常可观，在实现其工程价值的同时，安全管理工作也很重要。在混凝土坝的安全评价中，常采用数值模拟方法研究大坝的强度和稳定性[1]，而坝体和基岩力学参数选取的合理程度[2]，是影响数值模拟准确性的关键因素[3]。由于环境、施工、老化、徐变等复杂因素的影响，坝体及基岩材料的力学参数设计值往往与实际值存在一定差异[4]。如果仅利用原位试验或实验室测得的参数值，由于受到采样范围和试验条件等因素的限制[5-6]，有失代表性，难以满足数值分析的精度需求[7]。此外，现场取样试验成本高昂，工作复杂，并会对原有结构造成破坏[8]。因此，基于位移实测资料的参数反演分析方法成为识别混凝土坝坝体及基岩力学参数的主要方法之一[9]。在早期的研究中，Bonaldi等[10]使用确定性模型反演混凝土弹性模量、温度线膨胀系数；吴中如等[11]利用临界荷载法和小概率试件法，反演混凝土的断裂韧度；苏怀智等[12]采用遗传模拟退火算法反演大坝及坝基的物理力学参数。目前，反演分析方法主要有传统比值法和优化反演法。

由于难以建立待求解力学参数与观测数据之间的显性关系[13]，因此在反演过程中需要重复进行数值模拟过程[14]，其计算量大，耗时长，收敛速度慢，且容易陷入局部极值[15]。近年来，随着计算机技术的快速发展，遗传算法(GA)、粒子群算法(PSO)、BP神经网络等智能算法被应用于各类参数反演分析中[16]。龚晓雯等[13]利用改进的粒子群算法对混凝土重力坝参数进行反演；魏博文等[14]将ANSYS有限元程序作为子模块嵌入到自适应遗传粒子群算法(GA-PSO)中对混凝土坝参数进行反演；Zhuang等[17]提出了基于支持向量机回归的位移反分析模型，并将多群策人工鱼群算法(MAFSA)应用于隧洞参数识别；刘健[18]将BP神经网络进行改进之后，应用于坝体位移测值的处理，从而反演相应参数。在上述这些算法中，遗传算法经常遇到计算效率低、耗时长的情况，而粒子群算法在求解复杂函数的优化问题时，也往往出现收敛过早、精度低、鲁棒性差等问题[9]。

鲸鱼优化算法(WOA)是一种新型的群智能优化算法，具有参数容忍度强、鲁棒性好、收敛速度快等优点[19]。但是智能算法在求解寻优过程中都会表现出一定的趋同性，易使搜索过程陷入局部最优。一般情况下，增加种群数量，可以降低搜索过程的趋同性，提高求解质量。然而在串行计算中，增加种群会增加计算时间，使计算效率低下，尤其是在反演基岩和坝体力学参数问题中，优化算法需调用有限元程序进行正分析，当种群增大时，有限元计算的次数也随之增加，计算时间增加更加明显。利用计算机多核处理器进行并行计算，线程控制简单，计算环境稳定，正逐渐被应用于科学计算中。

本文将多核并行计算的方法引入标准的鲸鱼优化算法中，并加入非线性收敛因子，提出了基于多核处理器的并行式改进鲸鱼优化算法。在此基础上建立了重力坝坝体、基岩力学参数与监测点位移之间的非线性映射模型。最后，将所建立的模型应用于某混凝土重力坝力学参数反演中，并与用粒子群算法反演的结果进行比较，以验证该模型的可行性与合理性。

2 重力坝参数反演基本原理

在水压力、泥沙压力、扬压力、温度荷载等外部荷载的共同作用下，混凝土重力坝任一点的变形δ按照形成原因可分为水压分量δH、温度分量δT和时效分量δθ3部分[20]，即：

δ=δH+δT+δθ

(1)

将位移监测点实测变形数据进行回归分析，分离出水压分量、温度分量和时效分量，可建立重力坝位移的HST统计模型。在总位移δ中扣除温度分量δT和时效分量δθ之后，得到的水压分量δH为：

δH=δ-δT-δθ

(2)

(3)

坝体变形是最能直接反映大坝安全状态的指标，一般采用倒垂线作为混凝土坝变形监测的基准线，倒垂线从坝体廊道一直贯穿到基岩，不易受外界干扰，较为可靠，其测值可以反映坝基的变形。在外部荷载相似的情况下，材料的弹性模量是影响大坝变形的主要因素，当数值计算中采用的弹性模量与真实的弹性模量较为接近时，计算得到的大坝位移与实测值也比较接近；当数值计算中采用的弹性模量与真实弹性模量相差较大时，计算得到的测点位移值也会与实测值相差较大。

(4)

X=(Ec,Er)T

(5)

式中：P为荷载组合，有限元计算结果中的水压分量利用每年高水位和低水位情况下的测点位移差值计算得到。

对有限元计算得到相应测点的位移与从实测值中分离得到的水压分量求残差加权平方和，作为优化算法的目标函数。

(6)

式中：n为目标函数中包含的位移监测点的数量；wi为第i观测点测量值的权重。

坝体材料及基岩的弹性模量反演的步骤如下：

(1)对待反演坝段的倒垂线测值进行回归分析，建立重力坝位移的HST统计模型，分离得到δH、δT和δθ。

3 改进鲸鱼优化算法

鲸鱼优化算法(whale optimization algorithm, WOA)是Mirjalili等[19]在2016年提出的一种新型元启发式算法，该算法模拟了自然界中座头鲸捕食近水面鱼虾的行为，采用了29个测试函数对该算法进行测试，结果表明，无论是收敛速度还是寻优能力，与PSO、GA、GSA等算法相比，WOA算法均具有一定优势。

鲸是世界上最大的哺乳动物，座头鲸是鲸类中极具代表性的一种。座头鲸捕食时具有一种特殊的的狩猎方法，座头鲸在海洋中水深12 m左右的位置游动并搜寻食物[21]，当它发现目标后，会开始沿着螺旋轨迹上升，并不停吐出气泡，形成一张螺旋状的气泡网，将目标食物驱赶至近水面处的气泡网中心，最终捕获食物。图1为鲸鱼螺旋上升进行捕猎的示意图。WOA算法模拟了座头鲸捕猎的行为，通过迭代更新位置来搜寻猎物，最终获得最优解。

图1 座头鲸螺旋上升捕食方法示意图

在WOA算法中，每一头座头鲸都作为种群的搜索代理，可视为一个候选解。座头鲸捕猎过程分为3个步骤，分别为：猎物包围、气泡网狩猎策略、猎物搜索。

3.1 猎物包围

座头鲸可以识别猎物的位置，并对其进行包围。本算法以当前迭代循环最优解作为目标猎物模拟位置，其余搜索代理根据最优解进行更新迭代。该搜索行为的数学表达如下[22]：

D=|C·X*(t)-X(t)|

(7)

X(t+1)=X*(t)-A·D

(8)

式中：t为当前迭代次数;D为搜索代理与当前最优解的差值;X*(t)为当前最优解位置向量;X(t)为搜索代理当前位置向量;A、C为系数向量。

当某次迭代过程中，最优解发生更新时，X*(t)更新为新的最优解。

系数向量A由公式(9)、(10)得到，当|A|≥1时，算法会扩大搜索范围，进行全局搜索，寻找更优的猎物位置，由随机产生的搜索代理充当关键点，当|A|<1时，算法会缩小搜索范围，进行局部搜索，以当前最优解位置作为枢纽点。系数向量C由公式(11)得到。

A=2a·r1-a

(9)

(10)

C=2·r2

(11)

式中：r1,r2为[0,1]上的随机向量;Tmax为总迭代次数，在迭代过程中，a从2到0呈线性递减。

3.2 气泡网狩猎策略

座头鲸在搜索猎物过程中盘旋上升，在上升时吐出大量气泡形成陷阱，这个阶段被称为气泡网狩猎策略。公式(9)中a的减小即为收缩环绕机制的模拟，而A的波动范围也随着a的减小而缩小。环绕收缩过程可表达为：

X(t+1)=D′·ebl·cos(2πl)+X*(t)

(12)

式中：D′=|X*(t)-X(t)|，b为定义对数螺旋线形状的常数;l为[-1,1]之间的随机值。

座头鲸搜索猎物环绕收缩过程可建立如下更新模型：

(13)

式中：变量p为[0,1]之间的随机数，p随机在0和1直接切换。

3.3 猎物搜索

在初始阶段，座头鲸处于随机搜索猎物阶段，在数学表达中，利用随机数原理进行初始化。

D=|C·Xrand(t)-X(t)|

(14)

X(t+1)=Xrand(t)-A·D

(15)

式中：Xrand(t)为从种群中随机选择的搜索代理位置向量。

3.4 权重因子策略

在群智能优化问题设计中，如何平衡全局寻优能力和局部寻优能力一直是研究的热点。当算法的全局寻优能力较强时，搜索种群能尽可能广地分布在每个区域，从而能快速找出全局最优解所在的位置范围，而局部寻优能力较强时，能够从最优解的范围内对最优解位置进行快速定位，提高算法的收敛速度和结果精确程度。为了同时保证鲸鱼优化算法的全局寻优能力和局部寻优能力，引入权重因子w。当权重因子较大时，算法的全局寻优能力较强，有利于搜索代理进行全局搜索；当权重因子较小时，算法的局部寻优能力较强，能快速收敛，并保证结果精度。

(16)

利用权重因子进行对鲸鱼优化算法进行改进之后，座头鲸向猎物螺旋环绕收缩过程由公式(11)更新为公式(15)。

X(t+1)=

(17)

权重因子w随着迭代的进行不断减小，使搜寻过程从前期的有利于全局寻优过渡到后期的有利于局部寻优，鲸鱼优化算法改进前后权重因子变化过程对比如图2所示，由图2可见，改进后的权重因子w在搜寻初期下降较快，这样的趋势更有利于算法在确定最优解范围之后进行局部寻优，提高算法的收敛速度和结果精确程度。

图2 鲸鱼优化算法改进前后权重因子变化过程对比

3.5 算法的并行改进

随着计算机软硬件技术的发展，近年来，并行计算技术逐渐被应用于科学计算中。并行计算与串行计算的区别在于并行计算通过提高空间复杂度的方式来降低时间复杂度，提高算法效率，而串行计算保持算法顺序逻辑的简单，但需耗费更多的时间。并行计算中，计算线程数与CPU内核数相同，通过将任务分解为小部分，以并发方式，分配给CPU的各个内核同时进行计算，并最终将结果合并。标准鲸鱼优化算法反演力学参数的过程是串行计算，在一个迭代步中，一个搜索代理搜寻结束，下一个搜索代理才能启动搜寻，因此同一时间计算机只能进行一次有限元正分析，反演过程耗时极长。为了提高该优化算法反演混凝土坝力学参数的效率，利用并行计算的原理，提出基于并行式改进鲸鱼算法反演混凝土坝力学参数的方法，将搜索代理种群分为若干个子种群，将子种群分配到CPU不同内核上执行有限元正分析计算，待各子种群计算结束，将之合并形成新的种群，然后进入下一迭代步。

并行式改进鲸鱼优化算法流程如图3所示。

图3 并行式改进鲸鱼优化算法流程图

4 实例分析

某水利枢纽主要由碾压混凝土重力坝、泄水闸、引水建筑物、厂房和开关站等建筑物组成。水库校核洪水位151.88 m，设计洪水位150.00 m，正常蓄水位130.00 m，死水位130.00 m，总库容7.173×108m3。碾压混凝土重力坝坝底高程41.00 m，坝顶高程153.00 m，最大坝高112.00 m，坝顶宽度6.00 m，坝顶长568.00 m。选用该重力坝6#非溢流坝段作为实例进行验证，并选用该坝段60.00 m高程处倒垂线测点IP8的观测资料进行计算分析。

该实例计算所建立的有限元模型如图4所示，地基分别向上、下游延伸100 m，地基深度取为100 m。上游最高水位为150.00 m，最低水位为132.08 m。有限元模型采用四节点四边形单元，节点1 220个，单元1 135个，x向为顺河向，y向为竖直向。

图4 实例计算所建立的有限元模型

坝体材料为C20碾压混凝土，密度为2 400 kg/m3，泊松比为0.167，基岩材料参数如表1所示。荷载考虑上游坝面水荷载和库水重。地基侧面边界为法向约束，底面边界为固定约束。由于倒垂线IP8深入建基面以下40 m，故数值计算中以倒垂测点所在60 m高程上的点A相对于地基深度40 m处点B的位移作为IP8测点处的计算位移值。

表1 基岩材料参数

先以数值算例验证本算法的可行性。假设坝体混凝土材料弹性模量为30.0 GPa，基岩岩体弹性模量为15.0 GPa，其余参数均使用默认值。设定水位变化范围为130.00～150.00 m，水位变化间隔为1.00 m，使用有限元程序计算各水位工况下测点位移，由于外荷载只存在水压作用，因此有限元计算值全部为水压分量，可直接作为基准值。使用改进鲸鱼优化算法进行参数反演，图5为反演迭代过程。

图5 改进鲸鱼优化算法反演迭代过程

由图5可以看出，在迭代次数达到30时，适应度已收敛于0。反演结果为坝体混凝土弹性模量30.0 GPa，基岩岩体弹性模量为15.0 GPa，与预设模量相同，表明算法可应用于混凝土坝参数反演。

将算法应用于本工程实例中，利用实测数据进行反演。倒垂线IP8于2014年7月开始有测值，故分别选取2015、2016、2017、2018 共4个年份所对应的最低水位和最高水位值进行计算。表2为2015-2018年最低和最高水位对应测点的测值，以及最低和最高水位工况的水位差和测值变化。

对60 m高程的6#坝段倒垂线IP8测点的实测数据进行回归分析，分离得到δH、δT和δθ。在有限元计算中考虑上游坝面水荷载和库底水压力，分别计算每年对应最高水位和最低水位情况下，IP8测点对应位置的位移，计算其差值。

将有限元计算得到的水压分量差值与实测水压分量差值的离差平方和作为目标函数，分别用鲸鱼优化算法和粒子群算法求解最低适应度时的坝体和基岩弹性模量。为提高计算效率，对粒子群算法同样进行并行化处理，不影响其计算精度。图6为两种算法的反演过程对比，表3为2015-2018年6#坝段弹模反演计算结果。

由图6可以看出，粒子群算法较早进入局部极值，一直在局部范围内寻找最优解，而改进鲸鱼优化算法寻找的适应度比粒子群算法寻找的适应度小。从表3可得，由改进鲸鱼优化算法对2015-2018年倒垂线IP8监测数据进行反演得到结果的适应度分别为0.013，0.012，0.018，0.015，由粒子群算法反演得到结果的适应度分别为0.036，0.086，0.031，0.075，表明改进鲸鱼优化算法收敛速度比粒子群算法快，其反演结果的适应度也明显优于粒子群算法，且改进鲸鱼算法得到得反演结果更接近于该大坝的实际参数。

图6 改进鲸鱼优化算法和粒子群算法的反演过程对比

表3 改进鲸鱼优化算法和粒子群算法的2015-2018年6#坝段弹性模量反演计算结果

5 结 论

(1)通过对鲸鱼优化算法进行改进，结合有限元正分析，提出了基于改进鲸鱼优化算法反演混凝土坝变形参数的新方法，通过对某碾压混凝土重力坝实例进行验证，计算得到该重力坝坝体和基岩弹性模量，其反演结果优于粒子群算法，表明基于改进鲸鱼优化算法反演大坝力学参数的方法具有搜索精度高，反演速度快等特点，是合理可行的。

(2)基于多核CPU实现算法的并行计算，以子种群的形式对搜索代理种群进行细分，使复杂的计算能在CPU的多个内核上同时进行，极大提高了利用该算法进行复杂问题计算的求解速率，当程序在计算集群上进行计算时，可在极短时间内得到力学参数反演结果，可应用于大坝监测系统中以实时反馈各部位力学参数。

(3)多个力学参数同时反演时，解具有高维度的特点，如何继续改进该优化反演方法，使该算法在同时反演多个力学参数时，仍具有较高的搜索精度和收敛速度，是进一步研究的重点。