0.   引言

随着高土石坝建设的快速发展,越来越多的高坝已建或在建,比如糯扎渡、双江口、如美等,这些高坝主要分布在中国西南多地震区域,大坝的抗震安全计算尤为重要,其中抗震安全动力分析主要通过有限元方法,然而土石坝的筑坝材料动力模型参数又影响有限元分析的结果,因此动力模型参数的选择是否合理对于抗震安全十分重要。目前土石坝的动力计算模型参数一般是通过室内动三轴试验或者借鉴相关工程的参数来确定,然而由于室内试验的不确定性、缩尺效应和不同工程场地条件的差异性,上述方法确定的动力计算模型参数和实际值往往会存在差异。因此可以结合实测的地震波数据通过反演来确定动力模型计算参数。目前对于高土石坝强震记录监测到的较少,中国台湾集集地震中鲤鱼潭大坝^[1]、日本Tohoku地震中Fujinuma大坝^[2]都有记录,汶川地震中紫坪铺大坝只监测到了坝顶的数据^[3],相对于强震,弱震和余震资料比较丰富^[4-5],因此可以利用坝体余震或者弱震监测资料来反演土石坝的动力模型计算参数,为土石坝抗震安全分析提供依据。

对于大坝的材料动力参数的反演,目前对于拱坝,混凝土重力坝反演较多^[6-9],而对于土石坝的参数反演,其中主要是静力参数反演^[10-14],动力参数反演较少。朱晟等^[15]、汪旭等^[16]对紫坪铺大坝永久变形参数进行了反演,刘振平等^[17]对鲤鱼潭大坝动力参数进行了反演,他们反演计算时都是采用刚性边界模型,忽略了地基的相互作用和行波效应以及地震动非一致性输入的影响,因此其反演计算仍有待改善。并且这些学者反演时都是利用坝体及附近基岩场地监测到的强震记录进行反演,然而对于绝大部分的大坝很难获取到强震记录,这使坝体动参数反演受到一定的限制。

土石坝的地震反应分析一般采用刚性边界一致性输入的方法,但是随着坝高的增加,坝体尺寸变大,行波效应很明显。Zhang等^[18]通过对汶川地震中紫坪铺面板堆石坝震害调查,发现了地震的非一致性会引起坝体变形的差异及周边缝的不连续变位,从而威胁坝体的防渗,因此土石坝地震反应分析有必要采用非一致地震动输入方法。对于高土石坝考虑河谷地基以及坝体模型分析时,地基和坝体之间的相互作用影响,外行的散射波会向无限域地基辐射形成辐射阻尼效应,已有学者研究表明辐射阻尼效应对坝体反应有一定的影响^[19],因此反演分析时也有必要考虑坝基相互作用和辐射阻尼的影响。另外在河谷地形的自由场分析中,当已知的地震动信息为地表处时,截取地基底部的地震动输入一般是采取折减一半输入或者一维波动理论^[20]反演得到。但是由于地震波在河谷地形中反射和折射等影响,河谷自由场中和地表测点峰值、频谱特性及持时等特征均会发生一定变化^[21-23],因此当已知的是河谷中的地震动信息,再采用上述方法就不太合理,盛谦等^[24]、崔臻等^[25]将传递函数引入地下工程地震响应研究,并对响应的地震动输入频谱进行修正,蒋良潍等^[26]利用传递函数对不同地震波作用下的边坡动力响应分析做了研究。本文通过传递函数方法进行三维河谷自由场的地震动反演得到模型底部边界的输入地震动,其精度满足计算要求。

对于岩土工程的材料参数反演,大量的全局智能优化算法被用于解决该问题^{[10, 12-14, 17]},并取得了较好的反演结果。然而在参数搜索优化的过程中,需要多次进行有限元计算,对于包括河谷、地基以及坝体的高土石坝有限元模型来说计算量特别大。为了提高计算效率,径向基神经网络^[11]、响应面模型等^[27]传统机器学习方法被用来代替有限元计算,然而神经网络和响应面模型等方法都是基于经验风险最小化准则,存在收敛速度慢,局部极小点等缺陷,而支持向量机（support vector machine,SVM）基于结构风险最小准则,可以较好地解决小样本、非线性等问题^[28]。传统的SVM是单输出模型,在参数反演分析时,为了更准确地得到材料的整体性能,使用到的监测数据都会来自多个观测点,为了处理这种多输出问题,通常是独立地为每个输出建立单输出SVM,但是这样会增加计算时间,而且没有考虑各测点数据间的相关性从而带来误差。多输出支持向量机（multi-outputsupport vector machine,MSVM）^[29]可以有效解决以上不足。

本文提出了一种弱震情况下,三维河谷地形中考虑相互作用影响的高土石坝动力参数反演的方法。首先通过传递函数方法结合黏弹性人工边界模型解决了三维河谷地形的地震动输入,并考虑了坝基相互作用和辐射阻尼的影响,然后通过MSVM多输出支持向量机代替有限元计算提高了反演计算效率,并采用PSO粒子群算法进行参数优化选取。以糯扎渡心墙堆石坝为例,验证了所提出的动力参数反演方法的可靠性和有效性,可以为高土石抗震安全评价提供有效的动力计算依据。

1.   堆石料的静动力本构模型

静力本构模型采用工程上常用的邓肯-张EB模型,这个模型属于非线性弹性模型,是以弹性增量理论为基础,满足增量的广义胡克定律,即$\left\{\delta \sigma \right\}=$ $[D]\left\{\delta \epsilon \right\}$,弹性模量矩阵[D]中的切线弹性模量E_t和切线体积模量B_t是随应力水平而变化,其表达式如下：

式中 E_t,B_t分别为弹性模量和体积模量；R_f为破坏比；K为切向模量系数；n切向模量指数；$p_{\text{a}}$为标准大气压；c为凝聚力；$\phi$为内摩擦角；K_b,m分别为体积模量和指数；${\sigma }_1$,${\sigma }_3$分别为大、小主应力。

堆石料的动力本构模型采用基于等价黏弹性的等效线性分析方法^[30],这个模型主要考虑了动剪模量比G/G_max以及阻尼比$\lambda$随动剪应变幅$\gamma$的关系,反映了在动荷载作用下土的非线性与黏滞性两个应力应变关系特征。在具体计算时通过迭代方法,使得计算得到动剪模量G与阻尼比$\lambda$与等效剪应变幅${\gamma }_{\text{e}}$相适应,可以在一定程度上反映材料的非线性。

根据大量的试验结果,堆石料的最大动剪切模量G_max可表示为

式中,$p_{\text{a}}$为标准大气压,${{\sigma }^{\prime }}_{\text{m}}$为平均有效应力,${{\sigma }^{\prime }}_m=$$\frac{{\sigma }_1+{\sigma }_2+{\sigma }_3}{3}$；C,n分别为材料的最大动剪切模量系数和指数。其中,最大动剪切模量系数C和指数n在常规计算时可以根据动三轴试验得到,在本文中作为反演分析的研究对象。

2.   基于黏弹性人工边界结合传递函数的地震动输入方法

2.1   黏弹性人工边界和传递函数

黏弹性人工边界因为其易于与有限元程序结合,可以较好地模拟地基相互作用和辐射阻尼的影响,并且较好地模拟波在河谷地形中的传播,而被广泛应用^[31]。本文通过在截取的河谷模型外侧拓展一层无厚度单元来实现黏弹性人工边界,如图1,边界单元上的弹簧刚度系数和阻尼系数表达式采用文献[32]所述,地震动输入采用等效节点力方法,将输入地震动转化为直接作用于人工边界上的等效节点荷载^[33],表达式采用文献[33]所述。

图  1  三维一致黏弹性人工边界及边界单元

Figure  1.  3D consistent viscous-spring artificial boundaries and boundary unit

下载: 全尺寸图片 幻灯片

传递函数表达了被研究系统的频域形式的动力特性,从频域内描述了系统（三维河谷地形）对输入信号（地震波）的传递特性^[24],可以较好地分析地下工程和边坡的地震动力响应。本文通过传递函数方法进行三维河谷自由场的地震动反演得到模型底部边界的输入地震动,其精度满足计算要求。具体步骤如下：

（1）首先采用各频带频谱密度近似为常数的有限带宽白噪声^[24]（如图2）作为模型边界底部的地震动输入进行一次预计算,得到河谷中控制点的地震动信息。其中预计算是采用黏弹性边界模型进行计算分析。

图  2  有限带宽随机白噪声加速度时程及傅氏幅值谱

Figure  2.  Acceleration-time curves and Fourier amplitude spectra of random white noise with limited band width

下载: 全尺寸图片 幻灯片

（2）然后利用傅里叶变换,得到模型底部输入和控制点的加速度时程的傅里叶谱值$F(\omega )$,$F(\omega ,A)$,根据式（5）可得到控制点的传递函数：

式中,$F(\omega )$为输入加速度的傅里叶谱值,$F(\omega ,A)$为模型中控制点加速度的傅里叶谱值,$H(\omega ,\text{A})$为控制点的传递函数,A为模型中的控制点,$\omega$为圆频率。

（3）根据第二步得到的控制点传递函数,将监测的已知河谷地震动信息进行傅里叶变换,仍按照式（5）可以得到模型边界相应输入的地震动的傅里叶谱,将该傅立叶谱进行IFFT变换,就可以得到模型底部要求输入的地震波。

2.2   算例验证

图1（a）是糯扎渡大坝简化的三维河谷地形有限元模型,用八面体单元以及少部分五面体和四面体单元离散,其底面和4个侧面拓展一层无厚度单元作为一致黏弹性边界单元。计算采用线弹性模型,材料弹性模量为15 GPa,泊松比为0.25,密度为2.5 g/cm³。图中A点为河谷控制点,已知测点的地震动信息是景谷地震中距离坝轴线下游侧500 m河谷自由场处的测站的实测值,限于文章篇幅这里只用测站顺河向地震波进行计算验证。

运用2.1节所描述的方法可以得到模型底部输入的顺河向加速度时程,如图3所示,并利用图3加速度时程,进行正演计算,得到控制点A的计算加速度时程,图4是计算的和已知加速度信息对比,可以看出无论是加速度时程还是反应谱,吻合的比较好,精度基本上可以满足计算要求,因此当已知河谷中的地震动信息时,可以利用传递函数结合黏弹性人工边界的地震动输入方法来进行河谷地形自由场动力分析。

图  3  模型底部顺河向输入加速度时程

Figure  3.  Time histories of input stream direction acceleration at bottom of model

下载: 全尺寸图片 幻灯片

图  4  河谷控制点顺河向实测和计算加速度时程和反应谱

Figure  4.  Calculated and recorded time histories of stream direction acceleration and their response spectra at valley control point

下载: 全尺寸图片 幻灯片

3.   基于MSVM-PSO的反演方法

3.1   多输出支持向量机模型

一般的支持向量机往往是单一的输出模型^[34],比如一组训练样本{x_l,y_l},(l=1,2,…,k),x_l∈R^m是m个输入值,y_l∈R是一个输出值,即建立1个支持向量机模型,当求解岩土工程中多输出问题,在本文中输出值为坝体上各个测点3个方向的加速度反应谱值,即y_l∈Rⁿ是有n个输出值,则需要建立n个单一支持向量机模型,这样会增加计算量,并且容易忽略数据间的相关性。为了较好地解决这一问题,Fernando等^[29]提出了一种多输出向量机模型,其非线性回归映射关系和优化的目标函数为

式中$W=\left[w^1,w^2,\mathrm{...}\text{ , }{w}^n\right]$；$b=[b^1,b^2,b^3,\cdots ,$ $b^l{]}^T$；$u_i=\Vert e_i\Vert =\sqrt{e_i^T{e}_i},e_i^T=y_i^T-\phi {(x_i)}^{T}W-b^T$；C为惩罚因子。

为了求解上述优化问题的最优解,采用迭代加权最小二乘算法（IRWLS）进行计算^[29],可以得到

式中$\tau$为与W或b无关的常数项；$a_i=$ $\{\begin{array}{c}0\text{ ,}\\ 2C(u_i-\epsilon )/u_i\end{array}\begin{array}{c}(u_i\le \epsilon )\\ (u_i＞\epsilon )\end{array}$,$\epsilon$是一个很小的误差因子。

根据${L^{″}}_p(W,b)$对于$w^j$与$b^j$的驻点条件：$\frac{\partial {L^{″}}_p}{\partial w^j}=0$,$\frac{\partial {L^{″}}_p}{\partial b^j}=0$可以得到

式中,$D_{\alpha }=diag({\alpha }_1,{\alpha }_2,{\alpha }_3,\cdots ,{\alpha }_l),$$\Phi ={[\phi (x_1),\cdots ,\phi (x_l)]}^T,$ $\alpha ={\left[{\alpha }_1,{\alpha }_2,{\alpha }_3,\cdots ,{\alpha }_l\right]}^T,$$y^j=$ $\left[y_{j1},\cdots ,y_{jl}\right]$,I为单位矩阵,1为单位列向量。

根据表示定理可以得到

式中,${\beta }^j={\left[{\beta }_{j1},\cdots ,{\beta }_{jl}\right]}^T$。

将式（9）代入式（8）,可以得到

式中,$K_{ij}=k(x_i,x_j)=\phi {(x_i)}^T\phi (x_j)$为核函数矩阵。本文采用的核函数是经典Gaussian函数。

通过求解式（11）的线性方程组就可以求得w^j和b^j,进而可以求得多输出支持向量机的回归函数。具体求解步骤可以详见文献[29]。其中在训练MSVM时,核函数参数$\sigma$和惩罚因子C在较大程度上影响了MSVM建模的精度和泛化学习能力。本文采用网格搜索法并结合交叉验证对这两个参数的进行优化选择。

3.2   粒子群优化算法

粒子群优化算法是除了蚁群算法、鱼群算法之外的另一种群体智能的优化算法^[35],该算法最早由Kennedy和Eberhart提出的,算法首先在可解空间中初始化一群粒子,每个粒子都代表参数优化问题的一个潜在解,在PSO中,位置代表着当前粒子的优化数值,速度决定了粒子在迭代过程中的次数,适应度值是根据所定义的目标函数计算得到的。

假设在一个D维的搜索空间中,有n个粒子组成的种群$X=(X_1,X_2,\cdots ,X_n)$,其中第i个粒子表示为一个D维的向量$X_i={[X_{i1},X_{i2},\cdots ,X_{iD}]}^T$,代表第i个粒子在D维搜索空间中的位置,在本文中代表大坝坝料动力参数C和n的参数优化值。根据目标函数即可计算出每个粒子位置X_i对应的适应度值。第i个粒子的速度为$V={[V_{i1},V_{i2},\cdots ,V_{iD}]}^T$,其个体极值为$P_i=[P_{i1},$ $P_{i2},\cdots ,$$P_{iD}{]}^T$,种群全局极值为$P_g={[P_{g1},P_{g2},\cdots ,P_{gD}]}^T$。在每次迭代过程中,粒子通过个体极值和全局极值更新自身的速度和位置,更新公式：

式中$\omega$为惯性权重；d=1,2,…,D；i=1,2,…,n；k为当前的迭代次数,V_id为粒子速度；c₁和c₂是加速度因子,取值范围通常在（0,2）；r₁和r₂是分布于（0,1）之间的随机数。

在PSO方法对参数选取优化时,目标函数的选取是找到一组最优参数使得坝体上各个测点实测和计算的加速度反应谱数据点值相对误差较小,表达式如下：

式中 F(X)为目标函数；$X=\left\{x_1,x_2,\cdots ,x_d,\cdots ,x_D\right\}$为待反演参数,D为待反演参数的个数；N,L,T分别为测点数、加速度方向维数、反应谱周期的点数,$s_{ijk}$,$s_{ijk}^{*}$分别为测点计算和实测的反应谱值；M为质量矩阵,C为阻尼矩阵,K为刚度矩阵；$u$,$\dot{u}$,$\ddot{u}$分别为位移、速度和加速度向量；F为地震动输入的等效荷载向量。

4.   糯扎渡大坝坝料动力参数反演

4.1   工程介绍

糯扎渡水电站位于澜沧江下游普洱市思茅区和澜沧县交界处,装机容量5850 MW,该工程主体挡水结构为261.5 m高的心墙堆石坝,该坝是目前国内已建成的最高的土石坝,并且居同类坝型世界第三。坝顶高程824.1 m,坝顶长630 m,宽18 m,上游坡度比为1.9∶1,下游坡度比为1.8∶1,大坝水库正常蓄水位812 m。为了监测大坝受地震情况下的影响,在糯扎渡大坝上安置了若干个强震监测仪,其中河谷最大剖面上坝顶是T2测点,下游坝坡780高程处是T6测点,下游坝坡高程701高程处是T7测点,在距离坝轴线下游侧约500 m河谷自由场处测点T10,在2014年10月7日发生的6.6级景谷地震中均测到并记录了完整的加速度时程。

简化后的大坝典型剖面有限元模型如图5,其中堆石料Ⅰ和堆石料Ⅱ,由于考虑上游蓄水的影响,又分为上下游,因此具体分为心墙料,上游堆石料Ⅰ,上游堆石料Ⅱ,下游堆石料Ⅰ,下游堆石料Ⅱ。其中为了考虑坝基相互作用以及地基辐射阻尼的影响,向上下游、竖向、左右岸截取约1倍坝高距离计算分析,在截取地基4个侧面及底边界施加黏弹性人工边界。

图  5  糯扎渡大坝简化的典型剖面图

Figure  5.  Simplified typical section of Nuozhadu dam

下载: 全尺寸图片 幻灯片

三维有限元模型如图6,节点数37954,单元数量为36513。

图  6  糯扎渡大坝三维有限元模型

Figure  6.  3D FEM of Nuozhadu dam

下载: 全尺寸图片 幻灯片

4.2   动力模型参数反演

有限元的静力计算坝体采用邓肯-张EB模型,按照施工顺序分级加载,并逐级蓄水至景谷地震时坝体的蓄水位,坝基基岩采用线弹性模型材料参数同2.2节。坝体动力计算利用等价线性模型,地震动的输入采用第2节介绍的方法。在反演分析中计算目标函数采用的坝顶、780 m及701 m高程处的测点加速度计算值和实测值,河谷自由场中测点用来正演计算对比分析。

利用本文提出的反演方法进行糯扎渡心墙坝的动力参数优化,包括心墙料,堆石料,共需反演10个动力参数。坝体材料动力模型参数C和n的数值范围根据室内试验以及类似的工程进行取值,心墙料C值取1000～4000,堆石料C值取2000～5000,n值范围为0.4～0.7。根据正交试验设计生成81组数据,为了保证样本的多样性,又添加了69组随机数据,共150组。首先用样本数据对MSVM进行训练,模型输入值为150组参数组合,输出值对应为每个参数组合有限元计算的坝体3个测点3个方向的加速度反应谱值共9个,由于反应谱采用4 s间隔0.02 s共200个数据点,所以一共200个MSVM模型,模型初始的参数范围分别为${\sigma }^2\in (2^{-2},2^5),C\in (2^{-2},2^{13})$,根据网格搜索法和10折交叉验证,并以每个模型反应谱数值点预测值和计算值的RMSE最小为目标函数对其进行训练,图7是每个训练好的模型的RMSE值,可以看出除了个别训练的MSVM模型的RMSE值偏大,最大值为0.085,整体来看,训练的MSVM模型可以较好地满足计算要求。然后采用PSO粒子群优化算法进行糯扎渡大坝的动力参数优化,参数设置如下：粒子群规模20,维数即待优化参数个数10,加速度因子c₁=2,c₂=2,惯性权重$\omega$=0.7,终止条件是连续6代前后两次适应度值误差小于10^-5。

图  7  每个训练好的MSVM的RMSE值

Figure  7.  RMSE values of each trained MSVM

下载: 全尺寸图片 幻灯片

4.3   反演结果分析

反演得到的最优参数组合见表1。由表1可以看出,反演得到的最大动剪切模量系数C基本上都比试验值偏大,其中心墙料是室内试验值的1.52倍,堆石料Ⅰ平均值是室内试验值的1.50倍,堆石料Ⅱ平均值是室内试验值的1.32倍,其中下游堆石料的最大动剪切模量系数反演值明显比室内试验值明显偏大,分别为1.89倍和1.60倍,并且上游的堆石料反演值相对于下游值偏小,减小范围在35.3%～42.1%,分析原因主要是由于上游蓄水后导致上游堆石料的初始应力（大、中、小主应力）减小^[36],使得反演得到C值偏小。对于反演得到最大动剪模量系数n相对于室内试验值变化不是很大。

表  1  动力参数反演结果

Table  1.  Back-analysis results of dynamic parameters

动力参数		心墙料	坝料分区
			上游堆石料Ⅰ	上游堆石料Ⅱ	下游堆石料Ⅰ	下游堆石料Ⅱ
系数C	反演值	2661	2694	2300	4650	3555
	室内试验值	1754	2455	2216	2455	2216
指数n	反演值	0.412	0.565	0.440	0.544	0.651
	室内试验值	0.451	0.600	0.609	0.600	0.609

下载: 导出CSV 

| 显示表格

为了验证反演得到材料参数的合理性,利用最优参数组合进行糯扎渡大坝的正演动力计算分析。对比分析了坝体上3个测点实测和计算的加速度时程曲线和反应谱曲线,3个测点变化规律一致,限于文章篇幅,只给出坝顶T2测点,如图8所示,可以看出加速度时程曲线和反应谱曲线拟合良好,说明了反演得到的参数合理性。表2为坝体3个测点实测和计算得到的加速度放大倍数,可以看出,坝体顺河向和竖向放大倍数基本一致,随着坝高增加放大倍数变大,坝顶测点放大3倍左右,其中横河向加速度放大倍数计算值比实测值偏大,可能与计算时没有考虑行波效应等因素有关。文献[4]分析了紫坪铺大坝在汶川地震后监测到的余震和弱震信息,坝体的放大规律和糯扎渡大坝在景谷地震中的放大规律基本一致,符合土石坝的地震反应一般特征。

图  8  坝顶T2测站实测和计算的加速度时程和反应谱

Figure  8.  Calculated and recorded time histories of acceleration and their response spectra at T2 station of dam

下载: 全尺寸图片 幻灯片

表  2  坝体加速度放大倍数计算值与实测值对比

Table  2.  Measured and calculated acceleration amplifications

方位	T2		T6		T7
	计算	实测	计算	实测	计算	实测
顺河向	3.48	3.36	1.53	2.06	1.78	2.03
竖向	2.93	3.18	2.45	2.57	1.59	1.53
横河向	2.12	1.66	1.83	1.17	2.47	1.74

下载: 导出CSV 

| 显示表格

为了进一步验证黏弹性人工边界结合传递函数的地震动输入方法的合理性,对比分析了河谷自由场中的T10测点的实测和计算加速度时程曲线和反应谱曲线,如图9,从图9中可以看出拟合良好,说明采用2.2节无坝体的有限元模型得到的地震动信息基本上可以满足计算精度的要求,这也说明反演得到参数的合理性。

图  9  自由场T10测站实测和计算的加速度时程和反应谱

Figure  9.  Calculated and recorded time histories of acceleration and their response spectra at T10 station of dam

下载: 全尺寸图片 幻灯片

综上所述采用本文提出的地震动输入方法和反演方法可以较好地得到坝料的动力参数,并且精度可以满足计算要求。

5.   结论

本文提出了一种弱震情况下,三维河谷地形中考虑相互作用影响的高土石坝坝料动力参数反演方法。利用传递函数结合黏弹性人工边界方法较好地解决了河谷地形的地震动输入问题,并将多输出支持向量机模型和粒子群优化方法相结合,将其应用于糯扎渡大坝坝料动力参数反演分析,结果表明：

（1）取三维河谷自由场中测点为地震响应控制点,基于传递函数结合黏弹性人工边界得到模型输入的地震动,其计算精度满足要求,可以保证河谷地形动力响应分析中地震动输入机制的合理性。并且多输出支持向量机可以在保证计算精度的情况下,只需要相对较少的样本和训练次数,有效的节约有限元计算的时间。

（2）对于弱震情况下,考虑坝基相互作用影响的堆石料动力参数反演方法是可行的,计算精度基本满足实际工程要求,可以用于缺乏强震或震害资料的坝体参数反演,为现行土石坝动力参数确定方法提供一定的参考。

（3）通过对比本文反演结果、室内动三轴试验的最大动剪模量系数C值,室内试验值偏小。对于同种坝料,由于上游受到蓄水等原因的影响,反演得到上游堆石料的最大动剪切模量系数C值相对于下游减小,工程应用时应值得注意并予以适当的修正。