土石坝地震永久变形参数反演方法研究

汪 旭，康 飞，李俊杰

（大连理工大学 水利工程学院，辽宁 大连 116024）

1 引 言

堆石体地震残余变形模型参数选取的准确与否，对土石坝的抗震性能计算以及坝体地震安全评价具有重要意义。但堆石体永久变形参数一般由试验获得，由于受试验条件、施工工艺和施工质量等因素的影响，试验测得的参数往往与实际值存在较大的差异，使得数值计算的坝体地震残余变形与实测值偏离较大。因此，有必要利用坝体观测值进行堆石体地震永久变形参数反演，从而为评价坝体抗震安全与稳定提供依据。

常用的土体参数反演方法有遗传算法、蚁群算法、粒子群算法等直接法和神经网络法[1]。直接法将土体参数反演问题转化为优化问题，目前主要用于求解堆石体静力参数，取得了较好的效果，但存在循环迭代次数多、计算时间长的问题，尤其在有限单元数较多时，这些方法更显吃力。人工神经网络由于能够很好地映射岩土工程中复杂的非线性关系，也较多地应用于多参数反演问题中[2]。常用的神经网络模型是反向传播（BP）网络和径向基（RBF）网络。其中，BP网络是一种全局逼近网络，存在易陷入局部最优、训练速度慢、网络结构难以确定等问题。而RBF网络结构简单，学习速度快，具有局部响应特性，逼近能力更强。

本文提出了一种基于 RBF网络的土石坝地震永久变形参数反演分析模型，该模型在对永久变形参数进行灵敏度分析的基础上，建立了考虑参数灵敏度的网络训练目标函数，进一步提高了反演精度。以紫坪铺面板堆石坝为例，验证了所提出永久变形参数反演分析方法的有效性和实用性。

2 土石坝永久变形计算方法

2.1 静力计算

堆石体的静力计算采用邓肯E-B本构模型。切线模量Et和体积模量B按下式计算：

式中：Ei为初始切线模量；φ为内摩擦角；Rf为破坏比；K为切线模量系数；n为切线模量指数；Pa为大气压。

其中：

式中：Kb为体积模量系数；m为体积模量指数；φ0为当围压为一个标准大气压时的摩擦角；Δφ为围压相对于标准大气压增大10倍时的摩擦递减量。

2.2 动力反应计算

堆石体的动力计算采用等价黏弹性模型，用等效剪切模量G和等效阻尼比λ两个参数，反映土体动应力-应变关系的非线性和滞后性。对堆石料进行动三轴试验，测得动剪切模量比 G/Gmax和动阻尼比λ与动剪应变γ的关系曲线，并采用参考剪应变归一化，得到较为单一的 G/Gmax-γ/γr关系曲线和 λ-γ/γr关系曲线。根据试验结果，最大动剪切模量Gmax可表示为

式中：σ0′为平均有效主应力；Gmax、σ0′、Pa采用同一量纲；C为系数；n1为指数，其值由试验确定。

关于动水压力，本文采用Westergaard附加集中质量法，把动水压力对坝体地震反应的影响用一等效的附加质量来考虑，将其与坝体质量相叠加来进行动力分析。作用于节点i的附加质量为

式中：yi为从水面到节点i的水深；Ai为节点i的控制面积；H0为节点i所在断面从水面到库底的水深；ψ为上游坝坡与水平面的夹角；ρ为水的密度。

2.3 永久变形计算

坝体永久变形计算采用等价节点力法[3]，其基本思想是，由某一经验公式确定坝体各单元地震过程中的残余应变，然后把残余应变转换为单元的等效节点力，以此节点力进行静力计算所得的位移，代替单元残余应变对坝体永久变形的贡献。就残余变形计算模式而言，大多采用水科院模型[4]、谷口模型、沈珠江提出的永久变形模型[5]（以下简称沈珠江模型）以及这些模型的改进模型[6－8]。本文采用改进的沈珠江模型。

沈珠江模型[5]同时考虑了土体剪切变形和体积变形，且只需要一套参数就可以求得不同振次、不同动剪应变、不同应力水平下的残余变形，概念清楚，使用方便，应用较为广泛。残余体应变和残余剪应变表示为振次、动剪应变和应力水平的函数：

式中：Sl为应力水平；N为振动次数；c1、c2、c3、c4、c5均为试验参数；γd为动剪应变幅值；cvr、cdr为初始应力状态如应力水平Sl、动剪应变γd等的函数。该模型假定Sl对cvr无影响，即式中c3= 0。由此，本文对剩余4个参数进行反演。

改进的沈珠江模型[6]进一步探讨了应力水平对残余剪应变的影响，克服了沈珠江模型计算结果有时偏大的问题。在用有限元法计算地震永久变形时，将式（7）、（8）表示成增量形式为

式中：Δεvr、Δγr为残余体应变、残余剪应变增量。

3 永久变形参数反演径向基网络模型

3.1 径向基函数网络模型

RBF网是一种前馈神经网络，一般由输入层、隐含层和输出层3部分组成。图1所示为m-k-p结构的RBF网络，即网络具有m个输入，k个隐节点，p个输出。其中，为网络输入向量，W为权值矩阵,为输出单元的偏移量，y =(y1,… ,yp)T∈Rp,为网络输出向量，di为第i个径向基函数中心，为第i个隐节点的激活函数，常用的是高斯函数：

式中：δ为该基函数的扩展常数，其作用是调节高斯函数神经元的敏感程度。

图1 RBF网络结构图Fig.1 Sketch of structure of RBF neutral networks

RBF网络的第i个输出可表达为

式中：wij为输出层权值；δj为基函数的扩展常数。

利用 Matlab神经网络工具箱编写程序，RBF神经网络训练的关键问题是确定基函数中心、隐含层神经元数量、权值矩阵和扩展常数。前两者的确定采用的方法是每次增加一个神经元，使得网络误差最小的输入向量成为一个新的神经元（基函数中心），直到建立的RBF网络模型输出误差达到预定的要求。权值矩阵根据隐含层输出和网络输出采用最小二乘法计算[1]。另外一个重要的参数是扩展常数，扩展常数的调整有利于提高网络泛化能力，其值越大，函数拟合得就越光滑，但也并非越大越好，太大的扩展常数使输入空间的输出趋于一致，径向基函数趋于完全重叠。其取值一般取决于输入向量之间的距离。扩展常数和网络训练误差的确定采用试算法，兼顾网络对已有样本的拟合能力和对未知样本的预测能力，选取合适的扩展常数值和训练误差。

3.2 永久变形参数灵敏度分析

灵敏度分析可以检验参数的变化对模型输出结果的影响程度，从而对所考察的参数变量有比较直观的认识。Morris法[9]是目前应用较多的一种灵敏度分析方法，它通过对某一参数进行扰动，计算自变量对因变量的影响值。修正的Morris法将自变量以固定步长变化，取多个Morris法影响值的平均值作为灵敏度判别因子：

式中：S为灵敏度判别因子；Yi为模型第i次运行输出值；Y0为由初始参数计算的输出值； Pi为第 i次模型运算的参数值相对于初始参数变化的百分率；q为模型运行次数。

3.3 样本集生成

样本的质量在一定程度上决定着神经网络预测的精度。所以，选择的样本要准确、均衡，具有足够的代表性，能够正确反映模型系统的内在规律。常用的样本生成方法有均匀设计、正交试验设计以及随机生成。本文采用中心组合试验设计[10]和正交试验设计共同生成样本集。

中心组合设计由Box和Wilson提出[10]，是一种国际上较为常用的响应面法。以三因子设计为例，它的试验点分布如图2所示，由立方体角点、中心点和星点三部分组成，可以看作是五水平部分因子试验。星点位于立方体轴上，α表示星点到中心点的距离，通常 α= 2k/4（k为试验因子数）。由此可见，中心组合设计与正交全因子试验设计在角点的样本点是重合的，由这两种方法共同生成的样本点不仅可以考虑立方体内的情况，还可以考虑外部一定范围的情况，具有较好的样本分布密度和较强的代表性。

图2 三因子中心组合设计布点示意图Fig.2 Central composite design scheme for three factors

同时，为了提高网络的泛化能力，本文采用样本分集训练评价相结合，逐步调整中心集宽度（扩展常数）寻优的训练方法[11]。将样本分成训练样本集、验证样本集和测试样本集三部分，当训练样本集的误差性能指标满足后，根据验证样本集的误差性能指标确定最优网络。

3.4 永久变形参数反演分析步骤

永久变形参数反演分析步骤可概括如下：

（1）根据已有资料进行坝体静力和动力反应分析，为地震永久变形计算做好保障工作，同时为后续计算提供必要的数据：各单元应力水平Sl和动剪应变γd。

（2）确定参数反演范围，并做参数敏感度分析。

（3）生成样本集。用中心组合试验设计和正交试验设计生成训练集样本，另随机生成验证集和测试集样本。

（4）用训练集、验证集样本进行网络训练，调整扩展常数δ的值，验证集样本达到目标精度后，用测试集样本进行检测，直到测试集样本的预测精度满足工程要求。

（5）以坝体实际观测位移为输入向量，进行坝体永久变形参数反演。

（6）用反演获得的坝体参数进行正演计算，并对比分析计算结果与实际永久变形。

步骤（4）中，以位移值为输入向量，永久变形参数值为输出向量，样本均方根误差作为评价网络性能的指标。然而，利用神经网络进行永久变形参数反演，是一个多输入、多输出问题，如果单纯以均方根误差来评价网络性能，即

则不能体现不同参数的误差对沉陷量的影响。也就是说，由式（16）计算的总误差较小时，可能因为高灵敏度的参数预测误差相对大，而使正演计算结果误差很大。所以，本文考虑根据步骤（2）灵敏度分析结果，将灵敏度因子作为验证集样本误差权值，确定网络性能评价函数，使网络偏重对高灵敏度参数有更好的拟合精度：

式中：RMSE为均方根误差；n2为验证集样本数量；p为输出向量的维数；yij为第i个样本点的第j个参数输出值；tij为第i个样本点的第j个参数实际值；Sj为第j个参数的灵敏度因子。

4 实例分析

4.1 紫坪埔面板坝概况

紫坪铺大坝最大坝高为 156 m，坝顶全长为663.77 m，是我国新世纪建成的超过150 m高混凝土面板堆石坝之一，经历了远高于其设计水平的浅源近震——汶川大地震，坝体出现了明显损伤，为面板堆石坝的抗震研究提供了宝贵资料[12－13]。

目前，静力参数的获得方法相对成熟。刘小生等[14]已通过试验获得紫坪铺大坝相关参数，朱晟等[15]也通过免疫遗传算法反演了堆石体部分参数。综合采用前人的研究成果，得到大坝静力计算主要参数如表1所示。编号1～5分别代表垫层料、过渡料、主堆石料、次堆石料、覆盖层。混凝土面板的计算采用线弹性模型，其密度为2.4 g/cm3，强度等级为 C25，取弹性模量 E = 28 GPa，泊松比μ=0.167。

表1 静力计算主要参数Table 1 Principal parameters for static calculation

紫坪铺大坝没有坝址基岩实测加速度记录，据多位学者分析[16－17]，估计坝基实际地震动输入在0.5 g以上。朱晟参考于海英等[18]给出的地震衰减关系，考虑上、下盘效应，推求坝址基岩EW（坝轴向）、NS(顺河向）、UD（竖向）三向峰值加速度分别为0.52g、0.46g、0.43g。据文献[18]有关资料，地震测站分布与坝址的相对位置如图3所示。由此本文地震波选择距离坝址较近的卧龙测站（051WCW）记录，并沿用朱晟推求的加速度峰值，采用等比例法缩放获得紫坪铺大坝基岩输入加速度曲线见图4。

图3 “5·12” 汶川地震震中及测站分布Fig.3 Distribution of “5·12” Wenchuan earthquake epicenter and observation stations

图4 坝址基岩输入地震加速度-时程Fig.4 Acceleration-time histories of dam bedrock

4.2 静力与动力分析结果

堆石体三维有限元计算采用中点增量法，分20级加载以模拟施工、蓄水过程，计算地震发生前的大坝应力变形状态。有限元网格剖分见图 5，总单元数为6 772个，总节点数为6 846个。图6、7是大坝在震前运行水位下堆石区典型剖面主要计算结果等值线分布。

图5 大坝计算网格剖分Fig.5 3D mesh of dam for calculation

图6 蓄水期坝体垂直沉降等值线（单位：m）Fig.6 Contours of vertical displacement of the dam in storage period (unit: m)

图7 蓄水期坝体大主应力等值线（单位：MPa）Fig.7 Contours of major stress of the dam in storage period (unit: MPa)

根据动力计算结果，紫坪铺大坝震动基频约为1.65 Hz，坝顶峰值加速度分别为0.86g（顺河向）、0.74g（竖向）、1.36g（坝轴向），这与孔宪京等[17]对坝顶实测加速度反应分析结果基本一致，竖向加速度值偏小。三向地震作用下坝体顺河向加速度最大值位于下游面坝顶处，当加速度反应超过屈服加速度时，便可能发生堆石料滑动。这从定性上解释了紫坪铺坝顶下游面颗粒发生松动、滚落的现象。

4.3 永久变形参数灵敏度分析与反演分析结果

坝体不同分区的材料，料源相同，施工控制指标相近，动力性质不会有大的差异。为减少参数反演的数目，坝体堆石料采用同一组永久变形参数。文献[19]给出了国内外几座高面板堆石坝的永久变形参数的试验值，但不同工程实例的参数差异很大。因此，参考类似工程确定一组初始参数，并初步确定参数反演的范围如表2所示。

表2 待反演参数取值范围和初始值Table 2 Ranges of parameters to be inversed and initial values

以10%为固定步长对某一参数扰动，分别取其初始值±10%、±20%、±30%、±40%的变化率，而固定其他参数不变，对坝体最大沉陷量进行灵敏度分析，结果如图8、9所示。

图8 竖向沉陷量变化趋势Fig.8 Trend of vertical settlement

图9 参数灵敏度分析结果Fig.9 Analysis results of parameter sensitivity

用中心组合试验和正交试验设计共生成 90组样本作为训练集样本，另随机生成9组样本，其中6组作为验证集样本，剩余3组作为测试集样本。图 10给出了验证集均方根误差随扩展常数的变化情况，当扩展常数为1.6时，均方根误差最小，之后无明显降低。

图10 均方根误差与扩展常数曲线Fig.10 Curve of root-mean-square error and expand constant

为检验δ=1.6时网络预测精度，用测试集样本进行参数预测，结果见表3。由表可见，3组测试样本的预测误差均在 5%左右，且参数c5预测误差普遍较小，说明以式（17）作为误差评价函数对于提高高灵敏度参数预测精度效果显著。

为进一步验证预测参数对坝体沉陷量影响的误差精度，将预测参数代入有限元模型进行计算。限于篇幅，仅给出第2组样本典型剖面坝轴线上的正演计算结果与实际值对比情况，见图11。结果表明，由预测参数计算的坝体沉陷量与实际值接近，网络预测效果基本能够满足工程精度要求。

表3 测试样本参数预测结果与实际值比较Table 3 Comparison of predicted results and actual values of test samples’ parameters

图11 坝体沉陷量与坝高关系曲线Fig.11 Relationships of dam subsidence and dam height

以紫坪铺大坝实测地震变形为输入，用训练好的神经网络反演坝体堆石料永久变形参数，结果见表4。

表4 堆石料变形参数反演结果Table 4 Back analysis results of rockfill deformation parameters

用反演的参数计算坝体永久位移，结果如图12（图中所标数值为实测位移）。从图中可以看出，随着坝体高程增加，竖向沉陷量逐渐增大，与实测沉陷规律一致，二者数值大小也很接近。图13显示了有限元计算的大坝震后网格变形（已放大显示），坝体断面变小，上下游边坡均向内部收缩，体现了堆石体“震缩”特性。受上游水荷载作用，坝体最大水平向永久变形指向下游，但地震永久变形表现为以沉降为主。由此可见，计算结果定性上是合理的，说明了基于改进沈珠江模型的等价节点力法计算地震永久变形的可行性。

图12 坝体沉陷量对比（单位：mm）Fig.12 Comparison of dam subsidence (unit: mm)

图13 坝体网格变形图Fig.13 Deformation of dam grid

5 结 论

将径向基网络引入土石坝动力参数反演分析领域，提出一种地震永久变形参数反演的径向基网络模型，并用于紫平铺面板堆石坝地震残余变形参数的反演分析，结果表明：

（1）径向基网络模型用于土石坝地震永久变形参数反演简便、有效。此方法根据类似工程确定参数反演的范围，通过试验设计构造若干学习样本，网络训练完毕后，将实测地震残余变形输入网络即可获得坝体材料残余变形参数。

（2）现有土石坝动力计算理论和方法基本可以反映其抗震特性，但土石坝抗震研究是一个极为复杂的课题，地震动输入、计算模型等众多因素都可能导致计算值与实际值有较大出入，进一步做好动力分析方法研究，会使永久变形参数反演更有意义。

[1]康飞, 李俊杰, 许青. 堆石坝参数反演的蚁群聚类RBF网络模型[J]. 岩石力学与工程学报, 2009, 28(增刊 2):3639－3644.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－3644.

[2]YU Y Z, ZHANG B Y, YUAN H N. An intelligent displacement back-analysis method for earth-rockfill dams[J]. Computers and Geotechnics, 2007, 34(6): 423－434.

[3]刘汉龙, 陆兆溱, 钱家欢. 土石坝地震永久变形分析[J].河海大学学报, 1996, 24(1): 91－96.LIU Han-long, LU Zhao-zhen, QIAN Jia-huan.Earthquake-induced permanent deformation of earth-rock dams[J]. Journal of Hohai University, 1996, 24(1): 91－96.

[4]王昆耀, 常亚屏, 陈宁. 往返荷载下粗粒土的残余变形特性[J]. 土木工程学报, 2000, 33(3): 48－53.WANG Kun-yao, CHANG Ya-ping, CHEN Ning.Residual deformation characteristics of coarse-grained soils under cyclic loading[J]. China Civil Engineering Journal, 2000, 33(3): 48－53.

[5]沈珠江, 徐刚. 堆石料的动力变形特性[J]. 水利水运科学研究, 1996, 6(2): 143－150.SHEN Zhu-jiang, XU Gang. Deformation behavior of rock material under cyclic loading[J]. Hydro-Science and Engineering, 1996, 6(2): 143－150.

[6]邹德高, 孟凡伟, 孔宪京, 等. 堆石料残余变形特性研究[J]. 岩土工程学报, 2008, 30(6): 807－812.ZOU De-gao, MENG Fan-wei, KONG Xian-jing, et al.Residual deformation behavior of rock-fill materials[J].Chinese Journal of Geotechnical Engineering, 2008,30(6): 807－812.

[7]朱晟, 周建波. 粗粒土筑坝材料的动力变形特性[J]. 岩土力学, 2010, 31(5): 1375－1380.ZHU Sheng, ZHOU Jian-bo. Deformation behavior of coarse grained materials under cyclic loading[J]. Rock and Soil Mechanics, 2010, 31(5): 1375－1380.

[8]王玉赞, 迟世春, 邵磊, 等. 堆石料残余变形特性与参数敏感性分析[J]. 岩土力学, 2013, 34(3): 856－862.WANG Yu-zan, CHI Shi-chun, SHAO Lei, et al. Residual deformation behavior of rockfill materials and sensitivity analysis of parameters [J]. Rock and Soil Mechanics,2013, 34(3): 856－862.

[9]黄金良, 杜鹏飞, 何万谦, 等. 城市降雨径流模型的参数局部灵敏度分析[J]. 中国环境科学, 2007, 27(4): 549－553.HUANG Jin-liang, DU Peng-fei, HE Wan-qian, et al.Local sensitivity analysis for urban rainfall runoff modeling[J]. China Environmental Science, 2007, 27(4):549－553.

[10]BEZERRA M A, SANTELLI R E, OLIVEIRA E P, et al.Response surface methodology (RSM)as a tool for optimization in analytical chemistry[J]. Talanta, 2008,76(5): 965－977.

[11]DUAN Shao-jie, HE Chao, XU Li-xin, et al. A training method for improving the generalization performance of radical basis function networks[C]//Proceedings of the 3rd World Congress on Intelligent Control and Automation. Piscataway: IEEE Press, 2000.

[12]陈生水, 霍家平, 章为民. “5·12”汶川地震对紫坪铺混凝土面板坝的影响及原因分析[J]. 岩土工程学报, 2008,30(6): 795－801.CHEN Sheng-shui, HUO Jia-ping, ZHANG Wei-min.Analysis of effects of “5·12” Whenchuan earthquake on Zipingpu concrete face rock-fill dam[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(6): 795－801.

[13]杨泽艳, 张建民, 高希章. 汶川地震中紫坪铺面板坝抗震特性初步分析[J]. 水力发电, 2009, 35(7): 30－33.YANG Ze-yan, ZHANG Jian-min, GAO Xi-zhang. A primary analysis of seismic behavior and damage for Zipingpu CFRD during Wenchuan earthquake[J]. Water Power, 2009, 35(7): 30－33.

[14]刘小生, 王钟宁, 汪小刚, 等. 面板坝大型振动台模型试验与动力分析[M]. 北京: 中国水利水电出版社, 2005.

[15]朱晟, 杨鸽, 周建平, 等. “5·12”汶川地震紫坪铺面板堆石坝静动力初步反演研究[J]. 四川大学学报（工程科学版）, 2010, 42(5): 113－119.ZHU Cheng, YANG Ge, ZHOU Jian-ping, et al. Analysis on static and dynamic characteristics of Zipingpu CFRD under “5·12” Wenchuan earthquake[J]. Journal of Sichuan University (Engineering Science Edition),2010, 42(5): 113－119.

[16]陈厚群, 徐泽平, 李敏. 汶川大地震和大坝抗震安全[J].水利学报, 2008, 39(10): 1158－1167.CHEN Hou-qun, XU Ze-ping, LEE Min. Wenchuan earthquake and seismic safety of large dams[J]. Journal of Hydraulic Engineering, 2008, 39(10): 1158－1167.

[17]孔宪京, 邹德高, 周扬, 等. 汶川地震中紫坪铺混凝土面板堆石坝震害分析[J]. 大连理工大学学报, 2009,49(5): 667－674.KONG Xian-jing, ZOU De-gao, ZHOU Yang, et al.Earthquake damage analysis of Zipingpu concrete face rock-fill dam during Wenchuan earthquake[J]. Journal of Dalian University of Technology, 2009, 49(5): 667－674.

[18]于海英, 王栋, 杨永强, 等. 汶川 8.0级地震强震动加速度记录的初步分析[J]. 地震工程与工程震动, 2009,29(1): 1－13.YU Hai-ying, WANG Dong, YANG Yong-qiang, et al.The preliminary analysis of strong ground motion records from the Ms 8.0 Wenchuan earthquake[J]. Journal of Earthquake Engineering and Engineering Vibration,2009, 29(1): 1－13.

[19]顾淦臣, 沈长松, 岑威钧. 土石坝地震工程学[M]. 北京: 中国水利水电出版社, 2009.