，， ， ，

(1. 长江勘测规划设计研究院，武汉 430010；2. 武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072)

1 研究背景

高堆石坝变形控制是大坝安全运行的核心，其变形包括瞬时变形和流变变形。作为堆石坝的主要筑坝材料，合适的堆石体本构模型以及准确的模型参数是堆石坝应力变形分析及其工作性态评价的关键。堆石体模型参数的确定方法，目前常用的有试验法、工程类比法和经验法等。受试验试样尺寸限制、试样扰动、取样随机性以及缩尺效应等因素的影响，室内试验很难准确地测出原级配堆石体的材料参数，其应力变形结果难以真实反映坝体的工作性状。工程类比法和经验法虽然简单，但受到实际工程具体施工工艺、施工方法和施工质量的影响，预测精度较低。因此，基于现场监测资料的模型参数反演方法成为确定堆石体材料参数的一条有效途径。

岩土工程中的反演常采用直接法，将岩土体参数反演问题转化成优化问题，将计算位移值与监测位移值的最小误差函数作为目标函数，通过反复迭代逼近待定参数的最优值。在以往的研究中，堆石体本构模型参数的反演主要针对静力模型参数，也有对流变参数进行反演分析。例如，张社荣等[1]采用改进遗传算法对公伯峡面板堆石坝静力模型参数进行了反演分析；朱晟等[2]采用免疫遗传算法对三峡茅坪溪沥青混凝土心墙堆石坝的邓肯E-B模型参数进行了反演分析；董威信等[3]采用基于人工神经网络与演化算法的反演方法对糯扎渡邓肯E-B模型进行动态反演分析；周伟等[4]基于演化神经网络和序列二次规划算法的智能反演方法对堆石坝流变本构模型参数进行反演分析；常晓林等[5]采用改进粒子群算法对堆石体幂函数流变模型参数进行了反演分析；张雷等[6]根据公伯峡面板堆石坝坝体填筑竣工到蓄水前的沉降监测资料，进行反演分析得到相应的流变参数；朱晟等[7]在利用建设阶段变形监测资料反演堆石体静力模型参数的基础上, 对“5·12”汶川地震紫坪铺面板堆石坝动力模型参数进行了初步反演；马刚等[8]对静力本构模型参数和流变模型参数进行综合反演分析。在以上堆石体模型参数反演中，或仅对静力模型参数进行反演，或仅针对流变参数进行反演，或对瞬变和流变参数进行联合反演，但都将整个坝体的堆石体视为均一的流变材料，均没有严格按照坝体实际材料分区进行反演。其根本原因在于，堆石坝存在多个材料分区，不同分区的材料参数均不相同，若同时对不同材料分区的静力模型和流变模型参数进行反演，反演难度会急剧增加甚至导致反演失败。

本文通过考虑瞬时变形和流变变形的相互关系，提出了一种瞬变-流变解耦的参数反演方法，在对静力模型和流变模型参数进行敏感性分析的基础上，选择敏感性较强的参数作为待反演参数，基于改进遗传算法和RBF神经网络对瀑布沟心墙堆石坝进行了瞬变-流变参数解耦反演分析。

2 堆石体静力瞬变模型和流变模型

2.1 邓肯-张E-B模型

邓肯-张E-B模型[9]作为非线性弹性模型之一，因概念明确、形象直观、应用方便，在参数的确定和应用方面积累了丰富的试验资料和工程经验，在国内应用最为广泛。本文静力模型采用邓肯-张E-B模型，切线模量Et和体积模量Bt的计算公式分别为：

(1)

(2)

其中:

(3)

φ=φ0-Δφlg(σ3/pa) 。

(4)

式中：K为切线模量系数；pa为大气压；n为切线模量指数；SL为应力水平；Rf为破坏比；Kb和m为材料常数；c为材料凝聚力；φ为内摩擦角；φ0为当围压为一个标准大气压时的内摩擦角；Δφ为围压相对于标准大气压增大10倍时的内摩擦角递减量。

2.2 南京水利科学研究院Merchant流变模型

南京水利科学研究院通过对流变试验的研究，提出了基于指数衰减Merchant模型的七参数流变模型[10]。根据这一模型，流变变形、最终体积流变和最终剪切流变分别表示为：

εt=εf(1-e-αt) ；

(5)

(6)

(7)

式中：εt为流变变形；εf为最终流变量；εvf为最终体积流变；εsf为最终剪切流变；α为流变随时间的衰减指数；t为时间；e为自然常数；q为偏应力；b，β，mc，nc，d，lc为模型参数。

3 堆石坝瞬变-流变参数解耦反演方法

3.1 基本方法

在堆石坝的参数反演中，静力模型参数一般有10个甚至更多，流变模型参数有5～9个不等，不同材料分区材料参数均不相同，同时获得各个材料分区的静力模型和流变模型参数还很难实现。为了获得主要堆石分区的静力模型和流变模型参数，本文利用施工期、第一次满蓄水(蓄水位850 m)和第二次满蓄水(蓄水位850 m)过程，考虑材料之间以及各种参数之间的解耦关系，根据第一次满蓄水至第二次满蓄水期内的变形增量反演流变参数；得到流变参数后，根据全过程变形增量反演邓肯-张E-B模型参数。瞬变-流变参数反演方法见图1。

综合考虑灯桩建设周期（通常1年）、总投资、勘察设计时间及成本等因素，结合灯桩建设选址位置，根据不同的地质条件，选择几个代表性的地质类型进行勘察，同类型灯桩勘察设计互相借鉴，以此类推，编制灯桩标准化勘察报告。

图1 瞬变-流变参数解耦反演方法Fig.1 Decoupling inversion method of instantaneous and rheological parameters

3.2 可行性分析

为了分析瞬变-流变参数解耦反演方法的可行性，以瀑布沟心墙堆石坝为例，进行了考虑流变和不考虑流变2种情况下的应力变形分析。瀑布沟心墙堆石坝最大坝高176 m，上游坝坡1∶2～1∶2.25，下游坝坡1∶1.8，坝顶宽度14 m。0+240 m监测断面监测点布置及沉降变化分别见图2和图3。

图2 0+240 m断面监测点示意图Fig.2 Schematic diagram of observation points in 0+240 m section

图3 第一次满蓄水至第二次满蓄水期间考虑流变和不考虑流变沉降增量对比Fig.3 Subsidence increments in the presence and in the absence of creep from the first impoundment to the second impoundment

图3给出了第一次满蓄水至第二次满蓄水期间考虑流变和不考虑流变沉降增量对比。从图3中可以看出，是否考虑流变对第一次满蓄水至第二次满蓄水期间的变形影响很大。在不考虑流变时，各监测点变形增量趋近于0，可认为该时段的变形增量主要是由流变引起的，因此可通过该时段变形单独反演流变参数。

3.3 参数反演基本流程

在瞬变-流变参数解耦反演过程中，首先单独反演流变参数，在获得流变参数的基础上，利用施工-蓄水全过程变形增量对静力模型参数进行反演，2次反演流程基本一样，仅待反演参数不同。基于改进遗传算法[11](MGA)和径向基函数神经网络(RBF)的参数反演方法基本流程见图4。

图4 MGA-RBF参数反演流程Fig.4 Flow diagram of parameter inversion by MGA-RBF

由图4可知，参数反演的基本流程为:

(1)构建参数样本。为保证样本具有足够的代表性，利用正交设计助手生成正交设计样本n1组；同时，为保持样本的多样性，用随机函数生成随机样本n2组。

(2)有限元软件正演计算。对(n1+n2)组参数样本运用有限元软件进行变形计算分析，提取监测点沉降计算值，结合实测数据，计算适应度值。

(4)遗传算法演化寻优。利用训练好的RBF神经网络计算适应度值，通过MGA的选择、交叉、变异等遗传操作找到最佳个体，得到反演参数。

4 瀑布沟心墙堆石坝参数反演分析

4.1 参数敏感性分析

在堆石坝应力变形分析中，堆石体的内凝聚力c一般取为0.0。由卸载试验确定回弹模量系数Kur和nur，在数值计算中分别取n的1.0倍和K的1.2～3.0倍，因此对这3个参数可不进行敏感性分析。根据瀑布沟工程资料，采用修正的Morris法[12-13]对邓肯E-B模型中的7个主要参数和流变模型参数分别进行敏感性分析，选取敏感性较强的参数，作为待反演参数提供依据。修正的Morris法公式为

(8)

式中：S为灵敏度因子；q为模型运行次数；Yi为模型第i次运行的输出值；Y0为由初始参数计算的输出值；Pi为第i次模型运算的参数值；P0为初始参数。当S为正时，表示计算值与对应参数呈正相关，反之呈负相关，灵敏度因子绝对值越大则表示变形对该参数的变化越敏感。

E-B模型参数和流变参数敏感性分析结果见图5，可以看出，邓肯-张E-B模型参数K，n，Kb，m的敏感性较强，其他参数敏感性相对较弱。流变模型中，参数mc，nc，d对坝体变形比较敏感。其中mc和nc主要影响最终流变量εf；d影响最终体积流变量εvf。

图5 E-B模型和流变模型参数敏感性分析Fig.5 Sensitivity of parameters in E-B model and creep model

4.2 待反演参数确定及试验样本生成

根据流变模型参数敏感性分析，模型参数mc，nc，d对坝体的变形比较敏感，因此将mc，nc，d作为流变模型待反演参数。此外，虽然α对坝体变形敏感性较弱，但由于室内试验的流变收敛速率与实际工程差异较大，室内几个小时流变即稳定，但在实际工程中，流变长达几年，甚至几十年以上，因此将反映流变速率的参数α也作为流变模型待反演参数。本文使用在室内试验参数的基础上缩放λ倍的方法，对心墙、次堆石区以及主堆石区流变参数α，mc，nc，d进行反演。在归纳总结以往不同工程反演参数与试验参数比例系数的基础上，确定本次待反演参数的缩放系数，见表1。表中下标1，2，3依次代表心墙、次堆石、主堆石，各分区流变参数使用相同的缩放范围。在给定的缩放系数范围内，使用正交试验设计方法生成27组样本，同时使用随机法生成173组样本以保证样本的多样性和均匀性，共200组，见表2。

表1 流变模型待反演参数缩放系数取值范围Table 1 Ranges of the scale factors of parameters to be inversed in creep model

表2 流变参数反演样本集Table 2 Sample set of creep parameters in inversion

表3 E-B模型待反演参数缩放系数取值范围Table 3 Ranges of the scale factors of parameters to be inversed in E-B model

表4 堆石体E-B模型试验参数Table 4 Parameters for E-B model of rockfill

表5 堆石体流变模型试验参数Table 5 Parameters for creep model of rockfill

在获得流变反演参数的基础上，采用上述方法确定E-B模型待反演参数K，n，Kb，m的缩放系数范围以及样本集，限于篇幅，仅给出了E-B模型参数缩放范围，见表3。堆石体E-B模型试验参数见表4，流变模型参数由南京水利科学研究院和长江科学院根据室内流变试验得出，见表5。

4.3 目标函数

堆石坝力学参数反演的实质就是寻找一组参数使计算位移值最佳逼近实测值，对于堆石坝来说，由于监测点众多，这种最佳逼近为总体上和平均意义上的最好近似。目标函数采用监测点位移计算值与位移实测值的二范数，即

f(λ1,λ2,…,λn)=

(9)

4.4 神经网络训练及参数反演

瀑布沟心墙堆石坝三维有限元模型见图6。

图6 大坝有限元网格Fig.6 Finite element model of the dam

对流变参数样本集在高性能计算系统上运用ABAQUS有限元软件进行坝体变形分析，运用Python程序提取监测点沉降并计算适应度值，用于训练RBF神经网络。利用改进的遗传算法进行全局寻优，并根据最优缩放系数组换算得到流变模型反演参数。在获得流变参数的基础上，采用类似方法对邓肯-张E-B模型参数进行反演分析，得到邓肯-张E-B模型反演参数。瀑布沟心墙堆石坝流变模型和邓肯-张E-B模型反演参数见表6。

表6 瀑布沟心墙堆石坝瞬变-流变参数解耦反演结果Table 6 Decoupling inversion parameters of Pubugou high core-wall rockfill dam

图7 0+240 m断面部分测点实测值与反演参数有限元计算值对比Fig.7 Comparison between measured and calculated settlements of some monitoring points in the 0+240 m section

基于反演参数对瀑布沟心墙堆石坝进行应力变形分析，0+240 m断面部分测点实测值与反演参数有限元计算值对比如图7所示。由图7以看出，不同高程的监测点反演参数计算值和沉降实测值在数值和发展规律上均吻合较好。说明本文得到的反演参数是合理的，能够反映坝料的真实力学变形特性,因此本文提出的瞬变-流变解耦的参数反演方法用于堆石体参数反演是有效的。

5 结 论

(1)提出了堆石体瞬变-流变参数解耦反分析方法，在对沉降增量进行分析的基础上，论证了解耦反分析方法是可行的。

(2)对邓肯-张E-B模型参数及南京水利科学研究院七参数流变模型参数进行了参数敏感性分析，认为静力模型参数中K，n，Kb，m对坝体沉降较敏感，流变模型参数中mc，nc，d对坝体沉降较敏感。选择敏感性较强的参数，利用第一次满蓄水至第二次满蓄水期间的变形增量和全过程变形增量分别对流变模型和静力模型参数进行了解耦反演，得到了该方法下E-B模型和流变模型的反演参数。

(3)利用反演参数对瀑布沟心墙堆石坝进行了应力变形仿真分析，监测点沉降计算值与实测值在数值和变形规律上基本一致，说明该方法用于堆石体参数反演是有效的。