RCC坝力学参数反演不唯一性概率统计分析方法探讨

2019-07-29

长江科学院院报 2019年7期

(三峡大学水利与环境学院，湖北宜昌 443002)

1 研究背景

高寒地区修建的碾压混凝土(RCC)坝所处环境条件恶劣、年内和年际温度变化明显，导致大坝运行期工作性态复杂。而在室内开展混凝土力学试验时，试件尺寸效应、湿筛、标准养护和“点参数”等因素，导致大坝混凝土实际力学参数与室内试验参数存在较大差异。为了获得大坝当前真实的物理力学参数，目前基于原型观测资料对坝体及坝基物理力学参数进行反演分析已在混凝土坝工程上得到广泛应用。向衍等[1]根据监测资料，提出基于MSC.Marc平台的遗传算法，反演得到大坝坝体弹性模量和坝基变形模量；康飞等[2]提出了一种用于材料参数反演分析的混合单纯形人工蜂群算法；Gu等[3]将均匀设计法、偏最小二乘回归及最小二乘支持向量机引入碾压混凝土坝多组参数反演中；雷鹏等[4]针对混凝土大坝坝体和岩基参数的区间不确定性，构造具有区间分析功能的RNN(粗糙神经网络)模型，并运用该模型反演坝体和岩基区间参数；牛景太[5]将混沌遗传算法应用于碾压混凝土坝横观各向同性参数和渐变参数反演；Wang等[6]提出了一种基于大坝地震响应的参数反演分析方法；Su等[7]提出基于参数灵敏度的最优选择方法，并以某混凝土坝的坝体及坝基物理力学参数反演为例，验证了该方法的准确性和高效性；曹明杰等[8]将量子遗传算法引入碾压混凝土坝综合力学参数反演中；冯帆等[9]提出了以基于施工期大坝和基础监测的垂向压缩变形来反演施工期真实力学参数。

研究表明，当采用传统优化算法进行多参数反演时，反演结果依赖于初值的选取，优化结果容易陷入局部极值等。虽然人工神经网络、遗传算法等仿生算法引入反分析领域后，对推动反演分析的发展起到积极的作用，但是基于仿生算法反演的结果仍然存在容易早熟、计算量大和收敛慢等缺点。而且无论是传统优化算法还是仿生算法，均没有从理论上解决多参数反演不唯一性的问题。

《水工混凝土试验规程》(SL352—2006，DL/T5150—2001)在测定混凝土轴向抗拉强度时，同时测定抗拉弹性模量，以及在测定混凝土圆柱体抗压强度时测定抗压弹性模量。而在水电工程中，大体积混凝土的强度标准值定义为：按照标准方法制作、养护的边长为15 cm立方体试件，在设计龄期用标准试验方法测得的具有80%保证率的抗压强度；对于常态混凝土坝，设计龄期一般采用90 d；对于碾压混凝土坝，设计龄期一般采用180 d。考虑到大体积混凝土强度参数和弹性模量密切相关，而采用80%保证率的抗压强度作为标准值，其中一个原因在于不同试件测得的抗压强度存在一定的离散性，这与反分析参数不唯一有一定相似性。为了较好地解决反分析不唯一性问题，本文根据类比大体积混凝土强度标准值的定义，采用概率统计的方法对反演不唯一性数值解进行分析，并以高寒地区某碾压混凝土重力坝为例，将碾压混凝土作为横观各向同性材料考虑，采用工程类比法，引入包含周期项与滞后项的冻胀分量[10-11]，建立考虑冻胀变形的统计模型，分离出大坝真实的水压分量相对值；然后基于均匀设计[12]、神经网络及有限元法对坝体及坝基物理力学参数进行多次反演；最后对反演值进行概率统计分析，依据80%保证率确定坝体和坝基反演值，获得相对稳定的反演值。

2 高寒地区碾压混凝土坝正分析基本原理

高寒地区碾压混凝土坝正分析主要包括2方面：①碾压混凝土本构模型，即考虑碾压混凝土坝的施工特点，假设碾压混凝土为横观同性材料，按照横观各向同性本构进行有限元计算；②基于大坝实测变形分离出真实水压分量，即针对高寒地区大坝实测变形受冻胀因子的影响，首先建立考虑冻胀因子的变形统计模型，然后分离获得大坝真实水压分量。以下对这2方面主要内容的基本原理进行介绍。

2.1 碾压混凝土横观各向同性本构

工程实践表明，碾压混凝土坝结构复杂，其计算参数和计算模型难以确定，以致坝工设计难以与工程实际相吻合。据已有研究成果，碾压混凝土本构关系包括等效横观各向同性本构、内时损伤本构、正交各向异性本构、基于神经网络的剪切本构关系等[13]。对于碾压混凝土坝，典型的成层结构体系特点导致大坝结构性态分析不能简单套用常态混凝土坝的计算理论与方法。因此，本文为便于计算，将碾压混凝土等效为横观各向同性材料，建立空间直角坐标系，使坐标面xOy与横观各向同性体的各向同性面重合，z轴垂直于层面。由文献[14]推导，碾压混凝土满足的应力-应变关系如式(1)所示，设σ=[σxσyσzσxyσyzσzx]T，ε=[εxεyεzγxyγyzγzx]T，其中：σx，σy，σz分别为x，y，z方向的正应力；σxy，σyz，σzx分别为xOy面、yOz面、zOx面的剪应力；εx，εy，εz，γxy，γyz，γzx分别为正应力和剪应力对应的应变。

σ=Dε。

(1)

式中：E1为横观各向同性面(xOy平面)内的弹性模量；E2为垂直横观各向同性面方向(Oz轴方向)的弹性模量；μ1为横观各向同性面内的泊松比；μ2为垂直横观各向同性面方向的泊松比；G1为横观各向同性面内的剪切模量，不是独立的弹性常数；G2为垂直横观各向同性面方向的剪切模量。

2.2 考虑冻胀的变形统计模型

在长效服役期间，普通环境条件下的大坝变形主要受到水压、温度、时效的影响[15]；而当处于高寒恶劣环境条件下，且坝顶位置不设保温措施时，坝顶测点附近易受温度场周期性冻胀影响。从数学角度出发，借鉴丰满大坝(老坝)工程实例，在传统统计模型基础上引入包含周期项与滞后项的冻胀分量表达式[10-11]；同时为体现冻胀变形的影响，建立适应冻胀变形的统计模型，将冻胀因子考虑在模型内，对冻胀变形数据系列采用基于Heaviside阶跃函数的自适应模型，即得到高寒地区混凝土坝变形统计模型统一表达式为

δ=δH+δT+δθ+H(λ0-λ)δD。

(2)

式中：δH为水压变形分量；δT，δθ分别为温度变形分量、时效变形分量；δD为高寒地区需考虑的冻胀分量，包括周期项δD1和滞后项δD2，即δD=δD1+δD2。

Heaviside函数表达式为

式中：λ0为混凝土内基准温度，λ0=0 ℃；λ为温度。该函数“微分”是狄拉克δ函数，考虑在λ=λ0处有跳跃趋势，所以H(λ0-λ)称为阶跃函数。

周期项描述了大坝冻胀变形的复杂过程，非冻胀阶段周期项按最小值计算；滞后项反映了较长时段平均气温变化的影响，可以描述气温滞后引起的冻胀变形。各分量表达式为：

(4)

(5)

δθ=c1(lnθ-lnθ0)+c2(e-0.1θ-e-0.1θ0)+

c3(e-0.01θ-e-0.01θ0) ;

(6)

di2Icos[2πi(x-x0)/365]} ;

(7)

δD2=d1I20-10+d2I20-20+d3I20-30+d4I30-10+

d5I30-20+d6I30-30。

(8)

式中：ai为水压分量回归系数；重力坝m=3，拱坝m=4；H为监测日水深；H0为起测日水深；b1i,b2i为温度分量回归系数；m1为年周期、半年周期，m1=1～2；t为起测日至监测日的累积天数；t0为起测日至所计算时段起始日的天数；c1,c2,c3为时效分量回归系数；θ为从起测日开始的天数折算1%；θ0为从初始开始到起测日的累积天数折算1%；di1,di2为冻胀分量中的周期项回归系数；I为冻胀因子，以区别同一模型中的周期函数；x为时间长度；x0为同年监测序列分析的起始日到开始有负温的时日天数；d1,…,d6为冻胀分量中的滞后项回归系数；Ii-j为气温滞后作用因子，其中i表示气温滞后时段，j表示滞后时段前按平均气温计算时的天数，一般以周、旬、月为单位，选取滞后天数要根据坝内实测温度与滞后气温进行试算，该因子计算结果为正值时按零等效。

3 参数优化反演步骤

针对优化反分析不唯一性问题，本文基于大体积混凝土强度标准值的定义方法，考虑到大坝混凝土弹性模量和强度之间密切相关，基于不同工况下的实测变形进行多次物理力学参数反演分析，然后对反演结果进行数理统计K-S检验，确定反演参数的概率分布函数，进而依据80%保证率确定反演参数。

在进行每次物理力学参数反演分析时，本文将均匀设计、BP神经网络与有限元法应用于RCC重力坝及坝基物理力学参数反演，步骤如下。

Step1：结合已有类似工程及该工程试验资料，确定待反演物理力学参数变化区间，然后采用均匀设计法在待反演物理力学参数x={x1,x2,…,xn}的可能空间中构造参数取值组合，形成待反演参数若干取值集合。

Step2：建立碾压混凝土坝三维准有限元模型，选定相应的计算工况，将设计的待反演参数样本值输入到有限元模型中计算，得到相应的水压分量计算值。

Step3：将已知参数样本值与相应的水压分量计算相对值作为神经网络训练样本组，即将有限元计算的水压分量相对值作为输入，而将均匀设计得到的物理力学参数作为输出。

Step4：将大坝关键监测点统计模型分离的实测水压分量相对值输入到已训练好的神经网络模型，即自动得到相应的物理力学参数值。

Step5：将步骤Step4中的物理力学参数输入到有限元模型进行正分析计算，若计算值与实测值误差较小，则认为该力学参数为所求；若不满足，继续重复Step2—Step4，直至满足要求为止。

据此，可以得到基于均匀设计、BP神经网络及有限元法的大坝及坝基物理力学参数反演技术路线，见图1。

图1 物理力学参数反演技术路线Fig.1 Technical route for inversion of physical-mechanical parameters

4 实例分析

4.1 工程概况及监测布置

4.1.1 工程概况

某水利枢纽工程位于西北高寒地区，是一项大型水利枢纽工程。枢纽水库总库容24.19亿m3，水库正常蓄水位739.00 m，死水位680.00 m。枢纽主坝采用全断面碾压混凝土重力坝，坝长1 489 m，主坝最大坝高121.50 m，坝顶高程745.50 m。水利枢纽规模为Ⅰ等大(1)型工程，该枢纽坝址处多年平均气温2.7 ℃，历史上曾观测到极端最低气温达到-49.8 ℃，极端最高气温达到40.1 ℃，环境条件极为恶劣。

4.1.2 监测点布置

选取典型挡水坝段作为分析坝段，该坝段共布置3个正垂线监测点，所在地基深处布置一个倒垂线锚固点。即测点PL5-1布置在高程675.10 m处，测点PL5-2布置在高程706.50 m处，坝顶测点PL5-3布置在742.70 m处，坝基测点IP5布置在高程586.30 m处。典型坝段垂线布置见图2。其中监测点编号下方的数字表示正倒垂线的长度，单位为m。

图2 典型坝段垂线布置Fig.2 Vertical line of a typical dam section

4.2 考虑冻胀因子的变形统计模型建立

高寒地区建设的混凝土坝所处环境较为恶劣。如丰满大坝(老坝)的监测资料分析表明[14]，大坝所处环境温度年均值为6.38 ℃，最高温度为31.5 ℃，最低温度为-29.5 ℃。对比丰满大坝(老坝)坝址处环境条件可知，本工程所在地区环境条件较丰满大坝坝址环境条件恶劣。依据大坝监测资料和温度场仿真计算，坝顶附近出现负温，最低温度达到-17.8 ℃，则坝顶位移会受到周期性冻胀的影响；而PL5-1,PL5-2,IP5测点附近年周期温度为正温。为此，建立统计模型时，PL5-3测点考虑冻胀变形的影响，其他测点不考虑冻胀变形的影响。各测点统计模型计算系数见表1，由各测点建立的统计模型分离出各分量(相对于起测日)，见图3。其中，PL5-3测点的冻胀分量回归系数为：d21=-6.23×10-1，d22=0.00，d41=-9.12×10-2，d42=-1.26×10-1，d81，d82，…=0.00；d1=d2=d3=d5=0.00，d4=-4.32×10-2，d6=-2.32×10-2。

由图3可见：低温季节坝顶测点存在冻胀位移，非低温季节基本上呈现为一条直线段，无冻胀位移。非坝顶测点各位移分量呈现一般的变化规律。由表1可知，典型坝段各测点建立的统计模型复相关系数达到0.95以上，拟合效果很好，说明所建立的位移变形统计模型精度总体较高，基本上能拟合出大坝实际测值的变化规律，满足分析要求。

表1 各监测点模型计算系数Table 1 Values of model coefficients of monitoring points

图3 典型坝段各测点位移分量Fig.3 Curves of displacement components ofmonitoring points

4.3 碾压混凝土横观各向同性参数反演

4.3.1 有限元模型建立

模型计算范围：向上下游及地基深度方向取2倍坝高，坐标系x向为顺河向，上游指向下游为正向；y向为垂直向，铅直向上为正；z向为横河向，由右岸指向左岸为正。计算域上下游施加顺河向连杆约束，左右岸施加横河向连杆约束，底部施加完全位移约束。采用六面体八节点等参单元及少量的退化四面体单元进行网格剖分，网格剖分时尽量将测点布置在节点上。有限元模型单元数共13 620个，节点数17 274，其中大坝10 540个单元，典型挡水坝段有限元模型如图4。

图4 典型挡水坝段有限元模型Fig.4 Finite element model of typical dam section

4.3.2 参数上下限确定

假设碾压混凝土为横观各向同性材料，坝基为各向同性材料。考虑到碾压混凝土坝横观各向同性体5个独立参数中泊松比对坝体位移影响较小，故坝体与坝基泊松比不参与反演。实践表明，坝体横向弹性模量与竖向弹性模量比值为1∶2左右，结合地质情况及碾压混凝土试验资料，选定坝体横向综合弹性模量取值范围为20～65 GPa，横向综合弹性模量与竖向综合弹性模量比值为1∶2，竖向综合剪切模量取值范围为10～28 GPa，坝基综合变形模量Er为20～50 GPa，坝基泊松比取0.30，横向泊松比取0.168，竖向泊松比取0.18。

4.3.3 计算工况选取

针对反演分析存在不唯一性的问题，结合典型坝段实测资料，选取不同时段相应的上下游水位，共8种不同计算工况，如表2所示，每种工况下选取前后2个不同时刻上下游水位，进行多次参数反演得到多组反演参数值。每一种工况下的实测值均为前后2个不同时刻的相对位移值，即δ1=Δu′1-Δu1，δ2=Δu′2-Δu2，δ3=Δu′3-Δu3。其中，Δu1,Δu2,Δu3为由统计模型分离出的前一时刻正垂线测点PL5-1，PL5-2，PL5-3的水压位移分别相对倒垂线测点IP5的相对值；Δu1′ ,Δu′2,Δu′3为后一时刻的相对值。计算结果见表2。

4.3.4 反演模型建立

表2 各测点位移水压分量实测相对值Table 2 Relative values of measured displacement ofwater pressure component at each monitoring point

练样本，分别计算不同工况下的水压分量相对值，计算结果见表3。表中δ′1，δ′2，δ′3分别表示由有限元计算的测点PL5-1,PL5-2,PL5-3相对坝基测点IP5的位移。

选择3层神经网络模型进行反分析，将表3中测点位移水压分量相对值作为输入值，相应的坝体、坝基物理力学参数作为输出值，建立神经网络模型。为较好地防止计算过程中出现“过拟合”等问题，在网络训练前，对数据进行“归一化”处理。经多次试算，确定隐含层为10单元，传递函数采用S型正切函数，输出为purelin函数。将表3中的样本值输入到BP神经网络进行训练，经过4 133次学习训练，建立水压分量相对值与力学参数之间3-10-4的高度非线性映射关系。神经网络模型结构示意见图5。

图5 神经网络结构模型示意图Fig.5 Structure of neural network model

4.3.5 参数反演结果

基于上述建立的参数优化反演模型及选取的计算工况，按照第2节参数优化反演步骤进行计算，若不满足，则添加样本训练，直至得到合理的待反演物理力学参数。8种不同工况反演计算结果见表4。

表4 8种不同工况力学参数反演值
Table 4 Values of mechanical parameters obtained frominversion under eight different working conditions

工况E1/GPaE2/GPaEr/GPaG2/GPa148.9334.1745.2325.72236.2130.3343.2020.02333.2024.3034.2018.20447.3226.9045.8317.25547.8228.8442.4315.41648.4025.5037.5013.20744.2034.3032.7019.40840.1230.1037.5017.50

由表4可见，上述8种工况的反演结果不唯一，甚至不同工况下反演值存在较大的差异。考虑到大坝碾压混凝土抗压强度标准值采用80%保证率，而反演的弹性模量与大坝混凝土强度密切相关。为此，从数理统计角度考虑，依据80%保证率对表4获得的样本进行分析。依据数理统计K-S检验法可知，碾压混凝土横向综合弹性模量E1、竖向综合弹性模量E2、坝基变形模量Er、竖向综合剪切模量G2均服从正态分布，即分别满足E1～N(43.275,5.073 8)，E2～N(29.305,2.920)，Er～N(39.824,4.349)，G2～N(18.338,2.532)，然后分别按80%保证率取值，得到E1=39 GPa，E2=26.85 GPa，Er=36.16 GPa，G2=16.21 GPa。横向综合弹性模量与竖向综合弹性模量之比约为1.45，对比已有文献，由于本工程规模较大及环境条件较为复杂，故反演的坝体力学参数值偏大。

4.3.6 反演结果合理性判断

为验证反演计算的正确性，将得到的坝体及坝基物理力学参数输入有限元中进行正分析计算，得到计算相对值与实测相对值，相应的各测点位移值对比见表5。从计算结果看，有限元法正分析结果与水压分量实测值较为接近。由于实测值相对较小，一定程度影响反演精度。但总体来看，基于反演值的计算位移与各测点实测值误差基本上<5%，这说明本文反演的坝体及坝基物理力学参数是可行的，而从概率统计角度处理优化反分析不唯一问题，可以得到相对稳定的反演值。

表5 位移水压分量实测值与计算值对比
Table 5 Comparison of displacement of water pressurecomponent between measured value and calculated value

工况水压分量实测相对值/mm有限元计算相对值/mm相对误差/%δ1δ2δ3δ′1δ′2δ′3Δ1Δ2Δ310.721.542.160.751.562.054.171.305.0921.342.662.121.262.632.145.971.130.9430.921.890.750.951.950.83.263.176.6740.420.961.990.441.012.064.765.213.5250.571.212.090.561.162.181.754.134.3160.631.361.510.631.291.530.005.151.3271.002.061.210.941.961.196.004.851.6581.042.191.941.102.281.855.774.114.64