基于GRU-WPA法的重力坝地基力学参数反演分析方法

徐喆成，刘晓青，季 威，王雪红

(1.河海大学水利水电工程学院，江苏 南京 210024；2.南京市长江河道管理处，江苏 南京 210011)

0 引 言

大坝在建设期和运行期的监测数据是评价大坝安全性能的重要指标。利用监测数据对大坝及其地基的力学参数进行反演对大坝安全运行具有重要意义[1]。反演分析的目标函数往往是一个非线性多峰函数，监测数据与力学参数之间的显式关系是难以识别的，常采用正演优化反分析方式实现参数识别过程，而传统优化算法易产生局部极值问题[2]。对于许多需要在庞大的搜索空间中搜寻最优解的多极值实际工程问题，合理设计的智能优化算法可以实现搜寻全局最优。同时，地基力学参数反演分析过程中需要调用的正演分析模型往往是有限元模型，计算耗时较长，效率较低，在实际工程中的应用较为局限，而基于机器学习算法的代理模型来代替有限元模型，可以极大地提高计算效率，近年来在反演分析领域得到了越来越多的应用。代理模型的基本原理是通过寻找一组输入变量和响应变量之间的响应关系得出近似模型，从而替代真实系统快速给出求解，常用的代理模型有神经网络模型、SVR模型、MARS模型等[3]。

本文针对某实际工程大坝建立有限元模型，随机在一定范围内均匀抽取部分地基弹性模量和水压状态组合，利用自研程序进行有限元计算，得出各个状态下模型节点的位移，并将结果作为训练和验证样本对GRU神经网络模型进行训练，构建代理模型以替代有限元模型作为调用正分析模型，并结合狼群算法提高弹性模量的空间搜索能力，实现对大坝地基弹性模量的反演分析。

1 研究理论

1.1 参数反演理论

综合考虑多种外部荷载的作用，坝体顺河向位移disp可看作水压分量δH、时效分量δH以及温度分量δT三部分累加[4]，即

(1)

(2)

(3)

(4)

(5)

式中，βi、β1i、β2i、C1、C2、C3、C4均为回归系数；Hi为坝前水深的i次方；t为测日期数值；t0为起始日期；τ为随机误差项。

随着上游水位的升高，坝体承受的顺河向水体压力同步增长，相应的顺河向位移加大。因此水压分量δH的大小与水位的变动、结构以及地基材料弹性模量和目标位置的坐标存在着强关联性。混凝土坝本构模型数学表达为

uc=F(E，H，x，y)

(6)

E=[E1，E2，…，En]

(7)

式中，F(*)为有限元模型节点在不同的材料参数和环境荷载作用下与节点位移之间的映射关系；E为有限元模型各个区域不同弹性模量组成的向量；H为坝前水位高度(不考虑下游水位高度)；(x，y)为有限元模型目标节点的坐标；uc为目标节点的顺河向位移值。

反演的最终目标是找出合适的力学参数，使得通过有限元计算的目标节点位移值序列同实际测点测量所得位移值所经过分离得来的荷载分量的误差fe达到最小。误差fe为

(8)

式中，ucal为计算位移值；utrue为实测位移值。

1.2 GRU神经网络模型理论

图1 GRU神经网络重复模块

GRU神经网络模型的核心数学表达[5]为

rt=σ(Wr·[ht-1，xt])

(9)

zt=σ(Wz·[ht-1，xt])

(10)

(11)

(12)

式中，Wr、Wz、W为网络参数；σ为igmoid激活函数。

1.3 狼群算法理论

狼群算法结构流程如图2所示。

图2 狼群算法结构流程示意

狼群算法通过5个步骤来实现求解最优化问题[6]，分别是初始化狼群、游走行为、头狼召唤、围攻猎物以及狼群更新，其算法步骤为：

(1)初始化狼群。以D作为目标变量维数，N作为狼的数量，构建一个D×N的解空间。在解空间中初始化第i头人工狼位置Xi=(xi1，xi2，xi3，…，xiD)，i=(1，2，3，…，N)，若初始目标函数值为Yi=f(X)，则设使目标函数值最大的人工狼为头狼Ylead。

(2)游走行为。将一部分目标函数值较大的人工狼设为探狼，每个探狼将分别以活动步长stepa在解空间中的h个方向进行游走，每到达一个新位置都将重新计算该位置的目标函数值并返回起点重新出发，当发现某个位置的Yi>Ylead时或游走次数达到上限时，更新头狼位置。该过程的数学表达为

(13)

(3)头狼召唤。头狼会对猛狼进行召唤，猛狼会以步长stepb向头狼逼近，若在逼近过程中猛狼的目标函数值大于头狼的目标函数值，当头狼与猛狼的距离小于设定阈值时，进入算法下一步。猛狼的逼近过程计算式为

(14)

(4)围攻猎物。将头狼的位置标记为猎物位置，猛狼联合探狼对猎物进行步长为stepc攻，围攻过程为

(15)

式中，λ为[-1，1]之间的随机数。stepa、stepb、stepc的关系为

(16)

式中，S为步长因子；[mind，maxd]为第d个变量的取值范围。

通过以上5个步骤不断地往复循环，算法在迭代次数达到预设的最大值或算法的寻优个体达到精度要求时结束。

1.4 代理模型构建与弹性模量反演分析

利用系统抽样的方法抽取部分1.1节所述的计算分量，确定最佳的GRU网络参数，从而构建GRU代理模型，并结合狼群算法对地基弹性模量进行反演。

1.4.1 GRU代理模型构建

选取目标工程的典型断面，通过Hypermesh软件，建立该断面的有限元模型F。在合理的材料力学参数范围内，选取m组不同的参数组合，同时再从断面的坝底到坝顶的范围内随机选取n个水位高度H。

在有限元计算步骤中，通过依次选取不同的混凝土结构材料的弹性模量和不同高度的坝前水位，相应修改有限元模型F中的单元材料参数和水压荷载参数。经过有限元方法进行力学运算之后，根据目标节点的坐标(x，y)选取有限元模型对应位置的节点的计算位移值uc，依次选取所有目标节点的位移值之后，将计算结果和输入变量以[E，H，x，y⋮uc]的格式储存为样本，即

uc=F(E，H，x，y)

(17)

根据上述样本构造流程可知，样本[E，H，x，y⋮uc]中，[E，H，x，y]为自变量，确定了网络的输入节点的数量，[uc]为因变量，确定了网络的输出节点的数量。因为不同的网络结构训练结果影响较大，且该结构的确定不具备明显的推理性，需要人为选取不同的的网络结构等超参数，即中间各层的节点数量，以及各层的激活函数，根据损失函数的种类来试算不同的迭代步长α，最大迭代次数N，以及迭代阈值，经比较后选择最佳的网络结构状态。因为一批训练中，较大的样本数量可以促使网络快速收敛至较优参数状态，但是对于服务器的计算内存要求较高，而同一批次中较小的样本数量会使得模型训练时间加长，且网络收敛状态有偏差，所以根据实际GPU的内存大小来确定每批合适的训练数量m。在0～1范围内对网络的权值和阈值参数矩阵W进行随机初始化，为了防止网络失效，每个参数的初始化均需独立计算。将样本中自变量组成的向量[E，H，x，y]导入输入节点，进入后向传播阶段，通过每层节点之间的网络参数矩阵W对自变量进行组合计算，并从网络输入端传至输出端，得到每个样本对应的网络计算值O，即

O=f(E，H，x，y)

(18)

将O结合因变量uc计算本批次样本的损失值loss，即

(19)

最后，进入网络前向传播步骤，对各个网络参数W进行偏导计算，结合步长α对网络权重及阈值参数进行调整。进而继续重复后向传播和前向传播，直至N次迭代终止。当网络参数W达到当前最佳状态，保存网络结构及参数，即GRU代理模型构造完毕。

1.4.2 地基弹性模量反演流程

构建GRU代理模型后，根据狼群算法理论，先对弹性模量进行初始化设置，将初始弹性模量通过GRU代理模型计算得到位移值并与实际位移进行对比，若误差不满足要求且迭代次数未达到设定的上限，则更新作为头狼的弹性模量，在多次迭代寻优后将输出表现最优的头狼的位置定为最佳弹性模量，该计算流程如图3所示，该计算流程简称为GRU-WPA法。

图3 GRU-WPA法反演分析计算流程示意

表1 地基材料岩体力学参数

2 GRU-WPA法反演实例分析

选取某碾压混凝土坝为模型验证研究对象，该坝坝底高程41.00 m，坝顶高程153.00 m，最大坝高112 m，坝顶宽度6.0 m，上游折坡顶点高程为84.0 m，坡度为1∶0.3，下游面坡度1∶0.75，折坡点高程为145.00 m。该工程共计10个坝段。首次蓄水日期为2011年4月1日。地基材料力学参数如表1所示。在坝体压力以及地下水长期作用下，该工程的顺河向位移在监测期内呈现缓慢上升状态，因此对该处坝段地基材料参数应当重点关注，本文模型验证目标为地基弹性模量。

2.1 研究思路

本文计算选取该建筑物的某一个坝段，该位置地基高程为45.5 m，坝高107.5 m，坝基长88 m，选取地基长488 m，地基深300 m。针对该坝段的倒垂线进行分析，该仪器于2014年布置于靠近坝体上游面位置。选取该倒垂线从2014年7月25日到2019年10月31日监测的共计221个顺河向位移测值。选定地基主要由石英砂岩构成，材料密度为2 500 kg/m3，泊松比为0.167；坝体材料密度为2 400 kg/m3，泊松比为0.167。

本文抽取部分监测的位移值通过处理输入有限元软件进行状态组合计算，将计算结果结合假定的弹性模量作为神经网络的训练和验证样本，构造初始参数的GRU代理模型，如图3的步骤，利用狼群算法继续对弹性模量进行搜索寻优，从而得到最终的反演弹性模量。

为验证模型反演结果的合理性，本文利用工程常用的粒子群算法(PSO法)对该地基弹性模量进行反演分析，对比反演结果的差异。

2.2 GRU代理模型样本构造

因为地基的实际弹性模量未知，训练样本应包含可能的目标数值，所以根据表1所示地基可能的总体弹性模量，从3～10 GPa范围内随机均匀抽取200个不同数值。针对共计44 200种材料状态与水压状态组合，采用有限元软件GeHoMadrid对上述模型进行所有的状态组合计算，得出各个状态组合下模型节点的位移值disp，储存格式为[E，H，x，y，disp]，作为网络模型的训练和验证样本。

2.3 弹性模量反演

模型采用的多元线性回归模型见式(5)。根据该式提取水压分量构建目标样本格式为[E，H，x，y，δH]，其中，E为随机设定初始值；H为实测水位高度；(x，y)为点D有限元模型中的坐标；δH为上述模型提取的水压分量。

将取得的样本排列顺序进行随机打乱以后，所有数据按特征分别归一化到[0， 1]范围，抽取样本对网络进行训练验证，经过训练后保留网络参数以便后续调用。

经过如1.4节所示的过程得出训练完成后的GRU代理模型后，代替有限元模型进行基于狼群算法的弹性模量反演分析计算。利用系统抽样的方法选取图4中有限元模型测点A的221个不同水位的样本序列作为计算样本，其水位变化范围为125.33～145.96 m，迭代弹性模量预测位移与实测位移的平均绝对误差变化过程见图5。

图4 有限元模型

图5 WPA-GRU代理模型反演计算适应度迭代过程

由图5可知，迭代弹性模量预测位移值与实测值的平均绝对误差值从0.151 mm下降到0.137 mm，平均绝对误差在计算全程都不大，说明计算精度较高，且初定50次迭代的计算在迭代30次时数据已趋于稳定，主要是因为狼群算法具有较强的全局搜索能力。

为了进行对比，将同样的计算样本利用PSO法进行反演计算，其中设定粒子数为50个，个体学习因子为0.8，群体学习因子为0.2，惯性权重为0.5。由于每需要调用有限元模型50次迭代计算用时约4 min，而GRU-WPA反演分析的方法在30 s内即可完成50次迭代计算。

GRU-WPA代理模型反演计算得到的地基弹性模量为5.27 GPa，PSO法反演计算得到的弹性模量为5.16 GPa，结合不同的水位与两种不同反演弹性模量进行有限元计算，得到的计算位移结果与实测位移进行对比，位移值和相对误差结果如图6所示。由图6a可知，实测位移值与GRU-WPA反演弹性模量有限元计算位移值整体贴合度较好，仅在少数样本点处出现了较为明显的分离情况，说明这两组数据序列在相同水位状态下的数值较为接近；而PSO算法的反演弹性模量有限元计算位移值与实测数据整体变化趋势相近但有一定误差。由图6b可知，GRU-WPA反演弹性模量的计算位移值误差大部分都处于较低水平，大部分都远低于10%，部分峰值也不会超过40%，而PSO法反演弹性模量计算位移相对误差大部分略超10%，最大相对误差达到了68.9%，显然GRU-WPA反演分析计算结果具有更高的精度。

图6 实测位移与不同反演方法计算位移值的结果对比

3 结 论

(1)与PSO法相对比，为提高重力坝地基弹性模量反演计算的效率而提出的GRU-WPA代理模型的参数空间搜索能力会更高，由于不需要调用有限元模型，计算用时缩短，提高了计算效率。

(2)通过GRU-WPA反演分析计算得出的反演弹性模量与PSO法反演分析计算得出的反演弹性模量相比相差不大，说明其结果是合理的，且反演弹性模量的有限元计算位移值的相对误差整体更小，因此，GRU-WPA代理模型的计算精度是可以保证的。

(3)GRU-WPA代理模型在单参数反演中进行了尝试，后续可继续开展对于多参数高维度的反演计算研究，探究如何在高维度的情况下保证计算精度和求解效率。