基于响应面法的堆石体流变模型参数反演

2023-11-10解枫赞包腾飞李子阳张玉霞

人民黄河 2023年11期

解枫赞，包腾飞，李子阳，张玉霞，任杰

（1.河南省河口村水库运行中心，河南济源 454650；2.河海大学水利水电学院，江苏南京 210098；3.南京水利科学研究院水灾害防御全国重点实验室，江苏南京 210029）

堆石体变形除受材料应力状态的影响外，还受时间的影响。大量工程监测资料分析表明，混凝土面板堆石坝在蓄水运行后一定时期内会发生较大变形，说明堆石体具有流变特性［1］。堆石体的流变变形对面板及趾板等防渗体结构安全造成很大的影响［2］。例如：罗马尼亚里苏混凝土面板堆石坝（坝高为60 m）的流变变形导致面板产生大量裂缝［3］；天生桥混凝土面板堆石坝的流变变形造成面板脱空、产生裂缝［4］。因此，有必要研究堆石体的流变特性并预测其时空变化规律。堆石体流变变形的预测精度取决于流变参数的准确性，流变参数主要通过工程类比或室内试验获得［5］，但工程类比法精度低，室内试验只能采用缩尺三轴试验，无法准确获得原级配堆石体的流变参数。

采用实测变形资料反演堆石体流变参数是获取真实材料力学参数的主要方法之一。反演参数方法有单纯形法［6］、复合形法［7］、遗传算法［8］、蚁群算法［9］等。若待反演参数较多，采用上述反演参数方法时需要进行大量的有限元分析计算，费时费力。响应面法采用响应面模型代替有限元模型，可以大大提高计算效率。响应面法是一种数学统计方法，能够处理复杂系统的输入（基本变量）和输出（系统响应）的转换关系。当系统输入与输出之间为隐式函数关系时，可以通过响应面法建立两者间的显式函数关系［10－11］。本文利用响应面法建立流变参数与面板堆石坝运行期流变位移之间的响应关系，利用有限元正分析计算结果求解响应面函数的系数，建立流变模型参数与流变位移间的显式函数关系，基于面板堆石坝变形监测资料，建立目标寻优函数并求解，从而反演流变参数。

1 堆石体流变模型及响应面函数构建

当前，指数模型、九参数模型、Kelvin 模型、Maxwell 模型、Burgers 模型等流变模型常被用来描述堆石体材料的流变本构关系［12］。本文采用Burgers 模型，Burgers 模型由Kelvin 模型和Maxwell 模型组成，其形式为

响应面法处理复杂系统的输入与输出的转换关系，实际上是构造一个包含未知参数的多项式函数，以此替代隐式函数，借助统计方法求解最优的响应面函数［13－14］。因此，应用响应面法最关键的步骤是构造高精度、高效率的响应面函数。

综合考虑计算精度和计算量等，对于n个随机变量x1，x2，…，xn，通常采用不含交叉项的二次多项式来描述系统输入与输出间的响应关系，公式为

式中：s（X）为系统输出；X为系统输入，X＝{x1，x2，…，xn} ；a、bi、ci为待定系数。

对于堆石体流变参数反演分析而言，主要是用响应面函数建立流变参数与大坝流变变形间的非线性函数关系。参考堆石体Burgers 流变模型和响应面法的基本原理，构造含有指数项的堆石体流变变形响应面函数，表达式为

式中：sk（X）为测点k竖向位移的流变分量，xi为待反演参数。

然而，上述各反演参数的数量级相差较大，主堆石体材料与次堆石体材料的参数差异也很显著，如果直接将上述参数值代入响应面函数，必定会造成函数计算误差较大。因此，有必要对参数进行无量纲化处理，即将每个模型参数都与其初始值相除。需要指出的是，用于有限元计算的参数是未经无量纲化处理的，参数无量纲化只在求解响应面函数时进行，对最优反演结果进行无量纲化逆操作，即可得到真正的反演值［15］。

参数无量纲化处理公式为

综上，堆石体流变变形响应面函数形式如下：

2 基于响应面法的Burgers 模型参数反演分析

2.1 响应面函数求解

以测点P的流变变形为例，详细阐述响应方程待定系数的求解步骤。

1）对模型待反演参数进行组合，基于参数组合对各参数加、减扰动量，以此增加组合数量，共构造13 个参数组合。参数组合公式为

2）利用有限元法计算各组合中坝体的流变变形，将流变变形结果作为式（9）等号的左边项，形成测点P的响应面方程组。

3）若待定系数的个数为Q，那么其与待反演参数个数n的关系为Q＝2n＋1。采用上述参数组合方法恰好可形成2n＋1 个方程，则待定系数个数和线性方程个数相等，基于线性方程组唯一解的原则，求解得到测点P的唯一响应面函数。同理，可得到其他测点的响应面函数。

2.2 堆石体流变参数反演分析

确定各测点的响应面函数后，采用最优化方法确定响应面函数系数，目标寻优函数为流变分量计算值与实测值的均方根误差最小：

采用Levenberg－Marquart（L－M）算法求解目标函数的最小值，进而反演得出堆石体流变模型参数。L－M算法是介于梯度下降法与高斯－牛顿法之间的非线性优化算法，既具有梯度下降法的全局寻优特性又具有高斯－牛顿法的局部收敛特性，能够处理冗余参数问题［16］，对过参数化问题敏感性低，迭代收敛速度快，不易陷入局部最优。计算公式如下：

将反演得到的参数作为有限元模型参数，再次计算流变变形，比较流变变形计算值与实测值，若两者相差较大，则以反演参数为基准值重新进行模型参数组合、响应面方程求解和参数反演等，直到目标函数满足精度要求。

总结反演分析基本流程如下：1）采用有限元法计算面板堆石坝流变变形［17］，即求解流变位移的正分析过程；2）分析面板堆石坝运行期的流变位移数据，选择用于反分析的测点和时间段；3）基于Burgers 流变模型的数学形式，构造响应面函数；4）设计模型参数组合并进行有限元正分析，采用组合参数和计算的流变位移求解响应面函数系数；5）构造测试样本参数组，检验响应面函数的计算精度；6）定义参数反演的目标函数，采用L－M 算法搜索最优解，确定反演参数；7）由反演参数进行有限元分析得到特征点流变位移，与实测值比较，若两者相差较大，则以反演参数为基准值重复进行步骤3）～步骤6），直到目标函数满足精度要求。

3 工程实例

3.1 工程概况

某混凝土面板堆石坝由面板、垫层、过渡层、主堆石区、次堆石区、下游块石护面及坝前辅助防渗料区等组成。最大坝高120.0 m，坝顶高程760.0 m，坝顶长259.8 m、宽9.0 m，上下游坝坡分别为1 ∶1.40、1 ∶1.35。面板顶部厚0.3 m、底部厚0.7 m。水库校核洪水位759.1 m，设计洪水位756.2 m，正常蓄水位755.0 m，总库容2.65 亿m3，为多年调节水库。堆石坝于1998 年1 月13 日开始填筑，至2000 年9 月填筑至760.0 m，浇筑基本完成。水库于1999 年10 月11 日下闸蓄水，至2000 年9 月24 日库水位为745.08 m。

3.2 有限元模型建立

依据工程资料建立堆石坝三维有限元模型（见图1），模型包括面板、垫层、过渡层、主堆石区和次堆石区，共有17 854 个单元、19 545 个节点。

图1 某混凝土面板堆石坝三维有限元模型

3.3 特征点选取

在该混凝土面板堆石坝布置水平和垂直位移监测仪器，利用垂直位移进行反演分析。坝体垂直位移监测仪器为安装在最大断面的NSC－1 型水管式沉降仪，测点布置见图2。选取W2、W3、W6、W7、W10、W11 测点作为反演分析的特征点。这6 个测点位于坝高1／3～2／3 处，其沉降位移较大，能很好地反映堆石坝的沉降变形特征。

图2 坝体监测仪器布置

3.4 流变模型参数组合

流变模型参数组合见表1。组合1 为堆石体参数初始值，采用式（9）参数组合方法生成另外12 个参数组合。

表1 流变模型参数组合

3.5 有限元正演计算

该混凝土面板堆石坝于2000 年9 月24 日蓄水至正常水位745.08 m，在此之后水位维持在正常水位附近。本文选取2000 年9 月24 日至2004 年9 月17 日作为流变参数反演时段，该时段内起止2 d 的水位相同，2 d 的坝体竖向位移增量可认为是面板堆石坝4 a间竖向位移的流变变形。W2、W3、W6、W7、W10、W11 测点的实测流变分量分别为5.89、6.13、5.48、5.73、3.77、3.95 cm。利用表1 的13 个参数组合进行有限元分析，得到6 个测点在13 个参数组合的流变位移。

3.6 各测点响应面函数求解

按照式（6）和式（7）对模型参数进行无量纲化处理，将无量纲化后的模型参数和其对应的流变位移代入已定义好的响应面方程，构建方程组并求解待定系数。以测点W2 为例，W2 测点的响应面函数系数a、bi、ci组成待求矩阵A1×13，无量纲化后的模型参数构成矩阵M13×13，不同参数组合下的W2 测点流变位移计算值构成矩阵S1×13，则有：

由式（12）得到W2 测点的响应面函数系数，其响应面函数表达式见式（13）。同理，可得其他测点的响应面函数。