高堆石坝瞬变-流变参数三维全过程 联合反演方法及变形预测

马 刚 ，常晓林 ，周 伟 ，花俊杰

（1. 武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2. 武汉大学 水工岩石力学教育部重点实验室，武汉 430072）

1 引 言

堆石体的变形是长期而复杂的过程，其变形包括瞬时变形和随时间变化的流变变形，采用有限元等数值方法分析面板堆石坝的应力与变形已得到广泛运用。在面板堆石坝的应力与变形分析中，选用合理的堆石体本构模型以及准确的模型参数是整个分析的关键，堆石体的参数一般由室内或现场试验获得，然而受试验条件、缩尺效应的限制和堆石材料自身性质的离散性，使得测定的力学特性参数与实际值存在一定的差异，由此计算的堆石坝应力、变形与实测值差别较大，如室内流变试验一般几个月就完成了，而坝体的后期变形往往持续数年。因此，有必要利用坝体实测位移资料对堆石体的参数进行反演分析，并进行堆石坝后期变形预测。

以实测位移为基础的位移反分析方法已成为岩土工程学科的一个重要组成部分。以往，岩土工程中的反演方法常采用直接法，将参数反演问题转换成优化问题，文献[1]根据数座已建土石坝的原型观测资料采用复合形法进行了流变参数的反演分析。但由于该类优化问题的复杂性，传统的优化方法很难收敛到全局最优解，随着智能优化算法的飞速发展，遗传算法[2－4]、蚁群算法[5－6]、粒子群算法[7]、序列二次规划法[8]等全局智能优化算法被用于求解该问题，取得了较好的效果。在参数反演过程中需要反复进行有限元正分析，对于高堆石坝来说计算量非常大，而利用神经网络建立待反演参数与坝体位移之间的映射关系代替有限元正分析将大大提高计算效率[6－9]。径向基函数(radial basis function，RBF)网络是一种新型的神经网络模型[10]，与BP 网络相比，具有生物学基础和数学基础，隐节点具有局部特性，逼近能力更强，具有较好的泛化能力和函数逼近能力，可实现性能优越的非线性预测器。

在以往的研究中，堆石体参数反演分析主要是针对静力本构模型参数[2－7]，也有对流变参数进行反演的[1，8]。反演流变参数时，通常选择一个间歇期长的临时断面进行流变参数反演[8]，或者将满蓄以前的变形看作是瞬时变形，将满蓄以后的变形作为流变变形，以此为前提进行流变参数的反演分析。考虑到堆石坝实际施工过程和变形机制的复杂性，很难将流变变形与瞬时变形分开，因此，有必要对静力本构模型和流变模型的参数进行综合反演。本文在对静力本构模型参数和流变参数进行敏感性分析的基础上，选择对坝体变形敏感的参数作为待反演参数，采用基于粒子迁徙的粒子群算法和径向基函数神经网络构建反演平台进行参数反演，对水布垭面板堆石坝进行了参数反演分析，并预测了坝体变形稳定时的沉降值。

2 堆石体静力本构模型和流变模型

2.1 南水模型

南水模型[11]克服了邓肯-张E-B 模型不能考虑堆石体剪胀(缩)的缺点，能反映堆石体的压硬性和各向异性。南水模型假定应力空间中存在双屈服面，分别采用椭圆和幂函数表示，其屈服函数为

式中：p、q 分别为平均正应力和八面体剪应力；对于堆石体r、s 可取2。

南水模型中切线模量tE 的计算和邓肯-张模型相同。

式中：iE 为初始切线模量；c 为材料凝聚力；φ 为内摩擦角；lS 为应力水平，它反映材料强度发挥程度；fR 为破坏比；k、n 为模型参数；ap 为一个标准大气压；1σ 、3σ 分别为大小主应力。

南水模型采用抛物线描述体积应变与轴向应变的关系，可以统一考虑堆石体的剪胀(缩)性。采用切线体积比代替E-B 模型中的切线体积模量tB ：

式中： Rs= RfS； cd、 nd、 Rd为模型参数。

2.2 长科院幂函数流变模型

目前，在高堆石坝流变分析中常采用的是能模拟高围圧条件的幂函数流变模型[12]：

式中：sfε 、vfε 为某个应力状态下最终轴向流变量和最终体积流变量；τ 为时间；s( )ε τ 、v( )ε τ 分别为0 τ→ 时段内累积的轴向和体积流变量；sλ 、vλ为反映轴向和体积流变速率的参数。

最终轴向流变量sfε 和应力水平lS 、围圧3σ 的关系如下：

最终体积流变量vfε 和应力水平lS 、围圧3σ 可用线性函数拟合：

式(4)～(6)及参数c、d、η 、m、ac 、ad 、cβ、dβ、vλ 完整地表达了堆石体的流变特性。

堆石坝都是采用分层分区填筑，其加载过程复杂，要准确确定某一浇筑层具体的初始流变发生时间以及应力状态变化后的后续流变发生时间十分困难，因此，采用文献[1]的方法，采用相对时间代替绝对时间。

2.3 参数敏感性分析

南水模型有8 个参数，长科院流变模型有9 个参数，若都进行反演，则工作量很大，因此，有必要进行参数敏感性分析。修正的Morris 法采用灵敏度指标S 来反映参数的敏感性[6]：

式中：S 为灵敏度指标；iy 为模型第i 次计算输出值；0y 为初始参数对应的模型输出值；ix 为第i 次计算时参数；0x 为初始参数；n 为计算次数。

文献[6]对南水模型的参数进行了敏感性分析，认为k、φ、fR 、dc 、dn 、dR 对坝体的沉降比较敏感。本文采用一均质坝对长科院幂函数流变模型进行了参数敏感性分析，结果如图1 所示，可以看出，参数c、η 、cα、cβ、vλ 对坝体沉降较为敏感。

图1 长科院幂函数流变模型参数敏感性分析 Fig.1 Parameter sensitivity analysis of creep model

3 基于MPSO-RBF 网络的参数反演

3.1 基于粒子迁徙的粒群算法

在粒群算法[13]中，每个粒子代表待优化问题在多维空间中的一个潜在解。每个粒子具有位置和速度两个特征，粒子位置对应的目标函数值即可作为该粒子的适应度值。每个粒子根据它自身的“经验”和同伴的“经验”在搜索空间中向更好的位置“飞行”，直到在整个搜索空间中找到最优解或达到最大迭代次数为止。

在经典粒群算法(classic particle swarm optimization，CPSO)中，惯性权重ω 和加速因子1c ，2c 在优化过程中保持不变。在PSO 算法中，合理地调整全局搜索和局部开发的关系是提高算法性能的关键之一，较大的ω 值有利于粒子充分的在搜索空间中探索，较小的ω 值有助于粒子在当前位置的附近搜索。根据文献[14]提出的线性递减策略，在线性递减权重粒群算法中，在演化的早期阶段，采用较大的ω 值，在演化后期逐渐减小ω 值，惯性权重在优化过程中采用下式描述：

式中：iter 为当前迭代次数；i termax为允许最大迭代次数； ωmax、 ωmin分别为最大和最小惯性权重。

在群体智能算法中，提高算法性能的另一个关键是在演化的早期鼓励粒子在整个搜索空间“漫游”，不至于聚集在超级粒子周围而导致早熟收敛；在演化的后期，鼓励粒子“飞向”搜索到的最优值，以加速收敛。基于以上认识，文献[15]提出了线性变化的加速因子，在该策略中代表“认知”部分1c 随迭代次数逐渐减小，代表“社会”部分的2c 随迭代次数逐渐增大：

式中：1ic ，1fc ，2ic 和2fc 分别是初始和最终“认知”、“社会”加速因子。

受自然界物种迁徙能保持种群多样性的启示，本文采用一种新的改进的粒群算法(MPSO)。算法初始化为一群随机粒子，然后粒子被随机划分为若干个子粒群。每个子粒群独立演化，演化策略采用考虑了线性递减的惯性权重、线性变化的加速因子和自适应的变异算子。在演化的过程中，每隔若干迭代次数，进行一次粒子迁徙。粒子迁徙时，不仅将粒子目前的位置代入新的粒群中，还将粒子的个体极值pBest 引入新的粒群，以此加强子粒群间的信息交流并提高粒群的多样性。

3.2 RBF 神经网络

径向基函数(RBF)网络是一种两层前向神经网络，包括一个具有径向基函数的隐层和一个具有线性神经元的输出层，它能以任意精度逼近任意函数，具有较强的逼近能力，特别适合解决函数逼近问题，网络结构如图2 所示。

图2 RBF 神经网络结构图 Fig.2 Structure of RBF ANN

选用高斯函数作为径向基函数：

式中：σ 决定径向基函数的形状，σ 越大则基函数越平滑。

3.3 确定待反演参数

根据参数敏感性分析可知，k、φ、fR 、dc 、dn 、dR 、c、η 、cα、cβ、vλ 对坝体的变形比较敏感。由于堆石坝有多种材料分区，如果将上述所有敏感的参数均作为待反演的变量，工作量仍然较大，而且无法保证反演结果收敛到正确值，所以应该减少待反演参数，去掉测试方法较成熟的φ 和工程经验丰富的fR[2－3]。计算中的待反演参数共13个，分别为1k 、d1c 、d1n 、d1R 、2k 、d2c 、d2n 、d2R 、c、η 、cα、cβ、vλ ，下标1、2 分别对应主堆石和次堆石，主、次堆石采用相同的流变参数。

3.4 反演分析目标函数

堆石坝的参数反演就是寻找一组参数使计算位移值与实测位移值最佳逼近，由于堆石坝的测点众多，因此，上面所说的最佳逼近是指总体上和平均意义上的最好近似。目标函数可取为监测点的计算位移值与实测位移值差的二范数式，由于每个监测点的监测数据都是一个时间序列，故目标函数取为

式中： x1， x2， ，x13对应一组待反演的堆石体参数；m 为监测断面的个数；wi为第i 个断面的权重系数；ni为第i 个断面上监测点的个数；为第i 个断面上第j 个监测点在第k 个时间点的沉降计算值；为相应的实测值。

4 水布垭面板堆石坝变形参数的综合反演及变形预测

水布垭面板堆石坝，坝高为233 m，坝体和混 凝土面板的有限元网格分别见图3、4，均采用8 节点等参实体单元离散，面板与垫层之间，设置接触，接触模型为Coulomb 摩擦，摩擦系数取0.8。堆石体南水模型参数见表1，流变模型参数由长江科学院根据室内流变试验得出，见表2，反演参数取值范围见表3。

4.1 训练样本的选取

为了保证训练样本具有足够的代表性，用正交试验设计方法生成27 组样本，同时添加一些随机样本以保持样本的多样性与均匀性，共300 组样本，训练样本占总样本数的80%，测试样本占20%，样本集见表4。

图3 堆石坝有限元网格 Fig.3 Finite element model of CFRD

图4 混凝土面板有限元网格 Fig.4 Finite element model of concrete face

表1 堆石体南水模型材料参数表 Table 1 Parameters for double-yield-surface model of rockfill

表2 堆石料幂函数流变模型试验参数 Table 2 Experimental rheological parameters of rockfill

表3 待反演参数取值范围 Table 3 Ranges of parameters to be inversed

4.2 神经网络的训练

对样本数据进行归一化处理，将网络的输入数据限制在[0， 1]区间内。采用提前终止法和基于正则化的贝叶斯方法提高网络的泛化能力[10]，图5 是神经网络训练过程误差图，可见RBF 神经网络学习速度很快，到第6 次迭代时测试误差和训练误差均达到设计要求。

图5 神经网络训练过程 Fig.5 Neural network training process

4.3 参数反演

本文采用基于粒子迁徙的粒群算法对水布垭面 板堆石坝进行参数反演分析。算法参数如下：粒群规模为50，子粒群的个数 p= 4，粒群的迁徙率r= 0.4，迁徙间隔等于5，最大迭代次数取500，ωmin= 0.4，ωmax= 0.9， c1i= c2f=2.0， c1f= c2i=0.5，变异概率 pm= 0.2，iterMax = 10。目标函数中最大监测断面的权重系数为0.6，其余2 个较小断面的权重系数为0.2。图6 为反演过程中目标函数收敛过程。反演分析得到的最优参数组合见表5。

图6 反演计算收敛过程 Fig.6 Iteration process of parametric inversion

表5 堆石体参数反演结果 Table 5 Inversion parameters of CFRD dam

由表5 可以看出，反演得到的参数与试验参数差别较大，主堆石的k 值较试验值小很多，次堆石的k 值略小于试验值，反映最终流变量的参数c、cα、cβ较试验值增大，表明最终流变变形量增加，反映流变速率的η 、vλ 较试验值减小，表明流变变形速率变慢，坝体达到变形稳定所需要的时间增加。

基于反演得到的参数进行水布垭面板堆石坝的应力、变形分析，由图7 可以看出，监测点的沉降计算值和实测值在数值和发展规律上均吻合得较好。图8 为坝体在基本稳定期的沉降和水平位移等值线图，图9 为面板在基本稳定期的扰度图和轴向变形图，基于反演参数的计算结果表明，水布垭面板坝坝体沉降速率在2012 年底降低到12 mm/a，认为坝体变形达到基本稳定，此时坝体的最大沉降位移为2.63 m，流变变形约占坝体总变形的18%，流变变形显著，面板的最大扰度为0.9 m。

图7 坝体最大断面测点沉降实测值与计算值对比 Fig.7 Comparison between measured and calculated settlements of points in the largest section of dam

图8 变形稳定期最大断面位移图 Fig.8 Predicted displacements at largest section of dam in the deformation stable phase

图9 变形稳定期面板变形图 Fig.9 Predicted displacements of concrete face in the deformation stable phase

5 结 论

（1）本文采用基于粒子迁徙的粒子群算法和径向基函数神经网络构建参数反演平台，克服了粒子群算法易陷入局部最优和早熟收敛的缺点，采用经过训练的神经网络来描述模型参数和位移之间的映射关系，节省了参数反演的计算时间，提高了反演效率。

（2）对长科院幂函数流变模型进行了参数敏感性分析，认为参数c、η 、cα、cβ、vλ 对坝体沉降较为敏感。

（3）由于堆石坝填筑过程和变形机制复杂，很难将瞬时变形与流变变形分开，因此，对静力本构参数和流变参数进行综合反演。采用本文构建的参数反演平台对水布垭混凝土面板堆石坝进行了参数反演分析，基于反演参数的堆石坝应力变形分析结果表明，测点沉降计算值与实测值在数值和发展规律上均吻合地较好；在2012 年底坝体沉降速率降低到12 mm/a，认为坝体变形达到基本稳定，在变形稳定期坝体的最大沉降为2.63 m，约占坝高的1.12%，流变变形约占坝体总变形的18%，流变变形显著，面板的最大扰度为0.9 m，均在合理的范围以内。

[1] 沈珠江， 赵魁芝. 堆石坝流变变形的反馈分析[J]. 水利学报， 1998， 6： 1－6. SHEN Zhu-jiang， ZHAO Kui-zhi. Back analysis of creep deformation of rockfill dams[J]. Journal of Hydraulic Engineering， 1998， 6： 1－6.

[2] 郭雪莽， 田俊明， 秦理曼. 土石坝位移反分析的遗传方法[J]. 华北水利水电学院学报， 2001， 22(3)： 94－98. GUO Xue-mang， TIAN Jun-ming， QIN Li-man. Displacement back analysis of embankment dam using genetic algorithm[J]. Journal of North China Institute of Water Conservancy and Hydroelectric Power， 2001， 22(3)： 94－98.

[3] 张社荣， 何辉. 改进的遗传算法在堆石体参数反演中的应用[J]. 岩土力学， 2005， 26(2)： 182－186. ZHANG She-rong， HE Hui. Application of improved genetic algorithm to back analyzing parameters of rockfill[J]. Rock and Soil Mechanics， 2005， 26(2)： 182－186.

[4] 朱晟， 张美英， 戴会超. 土石坝沥青混凝土心墙力学参数反演分析[J]. 岩土力学， 2009， 30(3)： 635－644. ZHU Sheng， ZHANG Mei-ying， DAI Hui-chao. Back analysis of mechanical parameters for asphalt-concrete core earth-rock dam[J]. Rock and Soil Mechanics， 2009， 30(3)： 635－644.

[5] 田俊明， 周晶. 基于蚁群算法的土石坝土体参数反演[J]. 岩石力学与工程学报， 2005， 24(8)： 1411－1416. TIAN Jun-ming， ZHOU Jing. Inversing soil mechanical parameters of embankment dam using ant colony algorithm[J]. Chinese Journal of Rock Mechanics and Engineering， 2005， 24(8)： 1411－1416.

[6] 康飞， 李俊杰， 许青. 堆石坝参数反演的蚁群聚类RBF网络模型[J]. 岩石力学与工程学报， 2009， 28(增刊2)： 3639－3640. KANG Fei， LI Jun-jie， XU Qing. Ant colony clustering radial basis function network model for inverse analysis of rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering， 2009， 28(Supp.2))： 3639－3640.

[7] 李金凤， 杨启贵， 徐卫亚. 基于改进粒子群算法CHPSO-DS 的面板坝堆石体力学参数反演[J]. 岩石力学与工程学报， 2008， 27(6)： 1229－1235. LI Jin-feng， YANG Qi-gui， XU Wei-ya. Back analyzing mechanical parameters of rockfill based on modified particle swarm optimization CHPSO-DS[J]. Chinese Journal of Rock Mechanics and Engineering， 2008， 27(6)： 1229－1235.

[8] 周伟， 徐干， 常晓林， 等. 堆石体流变本构模型参数的智能反演[J]. 水利学报， 2007， 38(4)： 389－394. ZHOU Wei， XU Gan， CHANG Xiao-lin， et al. Intelligent back analysis on parameters of creep constitutive model[J]. Journal of Hydraulic Engineering， 2007， 38(4)： 389－394.

[9] 张丙印， 袁会娜， 李全明. 基于神经网络和演化算法的土石坝位移反演分析[J]. 岩土力学， 2005， 26(4)： 547－552. ZHANG Bing-yin， YUAN Hui-na， LI Quan-ming. Displacement back analysis of embankment dam based on neural network and evolutionary algorithm[J]. Rock and Soil Mechanics， 2005， 26(4)： 547－552.

[10] 张德丰. MATLAB 神经网络应用设计[M]. 北京： 机械工业出版社， 2005.

[11] 沈珠江. 土体应力应变分析中的一种新模型[C]//第五届土力学及基础工程学术讨论会论文集. 北京： 中国建筑工业出版社， 1990： 101－105.

[12] 程展林， 丁红顺. 堆石料蠕变特性试验研究[J]. 岩土工程学报， 2004， 26(4)： 473－476. CHENG Zhan-lin， DING Hong-shun. Creep test for rockfill[J]. Chinese Journal of Geotechnical Engineering， 2004， 26(4)： 473－476.

[13] KENNEDY J， EBERHART R C. Particle swarm optimization[C]//Proceedings of IEEE International Conference on Neural Networks. Piscataway： [s. n.]， 1995： 1942－1948.

[14] CLERC M， KENNEDY J. The particle swarm-explosion， stability， and convergence in a multidimensional complex space[C]//Proceedings of IEEE Transactions on Evolutionary Computation. [S. l.]： [s. n.]， 2002， 6： 58－73.

[15] YANG X， YUAN J， YUAN J， et al. A modified particle swarm optimizer with dynamic adaptation[J]. Applied Mathematics and Computation， 2007， 189(2)： 1205－1213.