基于实测资料和非侵入式随机有限元法的结构可靠度分析及其应用

季 昀，孟亚运

(1. 国家能源局大坝安全监察中心，杭州 310000；2. 中国电建集团华东勘测设计研究院有限公司，杭州 310000； 3. 中国电建集团中南勘测设计研究院有限公司，长沙 410000)

0 引 言

结构可靠度理论是研究定量描述系统不确定性因素的一种理论。一般认为，1946年A.M. Freudenthal 发表的《The safety of structures》标志着可靠度理论在结构设计领域研究的开始。之后，出现了考虑基本变量概率分布类型的一次二阶矩法(FORM)、二阶可靠度方法(SORM)、Monte-Carlo法等实用化方法。然而，尽管以FORM、SORM等方法为代表的结构可靠度算法得到了广泛的应用，但它们都仅适用于功能函数简单且为显式形式的情形[1]。事实上，现实中的大量结构系统属于复杂非线性系统[2]，其功能函数通常是隐式的，此时FORM、SORM等方法将不再适用。Monte-Carlo方法虽然同时适用于求解功能函数为隐式与显式时的结构可靠指标与失效概率，但由于它在分析小失效概率的可靠度问题时要求的抽样次数太多、计算效率太低而难以适用。

正是由于使用上述计算方法进行复杂结构系统可靠指标与失效概率求解有诸多不便，Faravelli[3]、Bucher与Bourgund[4]、Gavin与Yau[5]和Nguyen[6]等人将Box与Wilson[7]提出的响应面法(Response Surface Method, RSM)应用于结构可靠度分析领域，并取得了良好的效果。经典的响应面法是将原有结构系统的功能函数用序列二次多项式进行替代，随后以中心复合抽样或者拉丁超立方抽样等抽样方法获得样本点，应用回归方法确定该序列多项式的待定系数，最后针对获得的显式替代功能函数求解结构可靠指标与失效概率。从本质上来说，经典响应面方法的基本思路为采用序列多项式替代原有结构系统模型并拟合系统输入输出，这与近年来发展起来的Metamodel(或称surrogate model，model of model，代理模型)方法在基本思路上是一致的，只是相对而言，Metamodel的概念内涵更加宽泛一些[8]。除了经典响应面法外，Metamodel还包括神经网络模型[9,10]、Kriging模型[11]、支持向量机[12,13]以及随机响应面等多种近似模型。

随机响应面法凭借其明确的物理意义和数学意义上的收敛特性(在功能函数平方可积的条件下)，近年来在结构可靠度分析领域得到了一定的推广和应用[14,15]，但目前在水工结构工程中应用尚不多见。将随机响应面法与工程结构分析的有限元法采用“非侵入式”方法相结合，可得非侵入式随机有限元法，这有别于20世纪末提出的随机有限元法[16]。此外，目前结构可靠度分析中，结构的物理力学参数大多直接选取设计值进行分析，然而，结构在运行期的物理力学参数往往较设计值有一定差别。鉴于此，本文提出了一种基于实测资料反演力学参数与随机响应面法相结合的结构可靠度计算模式，并将其应用于重力坝结构可靠度分析。

1 结构可靠度随机响应面法简介

随机响应面法与经典响应面法最大的区别在于其响应面函数形式为服从标准正态分布的随机变量构成的Hermite正交多项式，而非直接由输入随机变量构成的序列多项式。1938年，Wiener首次采用随机多项式(Wiener's polynomial chaos)来研究系统随机过程的基本理论，这可以看作是有关随机响应面研究的最早案例。Cameron和Martin[17]从数学理论上证明，任意一个随机过程的二阶矩可以分解为服从标准高斯分布的随机变量组成的有限项多项式，并且该多项式是收敛的。Hermite正交多项式正是最早被证明具有该特性的多项式，因此，一般情况下，Hermite多项式又直接被称作随机多项式。

一维Hermite正交多项式的具体表达式为：

(1)

其递推公式为Hi(x)=xHi-1(x)-(i-1)Hi-2(x),i≥3，即H0(x)=1，H1(x)=x，H3(x)=x3-3x，H4(x)=x4-6x2+3。

根据Cameron-Martin定理[17]，对于平方可积空间L2(R,μ)内的任意函数f(x)，都存在一个Fourier-Hermite展开，满足：

(2)

式中：fi为对应Hermite多项式Hi(x)的系数。

当式(2)中的变量x用标准正态随机变量ξ替代时，可将其改写为：

(3)

式中：fi为对应Hermite多项式Hi(ξ)的系数，可以通过适当的抽样方法获得一定的样本点，然后由回归方法求得。

由此可见，任意平方可积空间内的函数都可以用Hermite多项式进行近似表达，且选取的Hermite多项式阶数越高，则逼近误差越小，因为Hermite多项式构成了平方可积空间的正交基，这一性质是经典响应面采用的序列多项式所不具备的。此外，由于一般工程中所遇到的问题中的函数基本都是满足平方可积这一特性，这也使得将其用于结构系统响应或者是功能函数的拟合成为可能。

总体说来，随机响应面法就是通过采用随机多项式展开来近似拟合结构的随机输入与随机响应之间关系的方法。采用随机多项式拟合功能函数后，可采用FORM、SORM、几何法等传统的结构可靠度计算方法求解相应的失效概率和可靠指标。具体计算流程如下。

(1)收集有关信息和数据，分析结构可靠度问题类型，确定输入随机变量向量x。将输入随机变量转换为独立的标准正态变量随机向量ξ。

(2)选用合适形式的Hermite展开多项式构成随机响应面中的响应面函数，以此近似拟合原结构系统的复杂的功能函数。

(3)选用合适的配点方法获得随机响应面函数的概率配点。将各配点作为原模型的输入，计算原模型相应的输出值。

(4)用回归方法计算随机响应面函数中的待定系数，获得随机响应面函数的最终表达式。

(5)针对上一步确定的随机响应面函数表达式，采用几何法或Monte-Carlo模拟法等方法求解结构可靠指标与失效概率。

(6)计算结束。

2 基于实测资料和非侵入式随机有限元法的结构可靠度分析模式

本文提出基于实测资料和非侵入式随机有限元法的结构可靠度分析方法，主要分为力学参数反演、构建功能函数、结构可靠度求解3个方面。一般地，重力坝结构可靠度分析包括抗滑稳定可靠度、坝踵抗拉可靠度以及坝趾抗压可靠度等，以下以重力坝的坝踵抗拉可靠度(以下简称重力坝抗拉可靠度)为例进行分析。具体操作流程如下：

(1)建立位移监测点实测值的统计模型，分离出其中的水压分量。大坝位移统计模型一般如下：

δ(δx、δy或δz)=δH+δT+δθ=

(4)

(2)根据现场试验得到的坝体与地基弹性模量取值范围，确定若干对坝体与地基弹性模量的组合。

(3)建立坝体与地基结构分析的有限元模型，分别计算第(2)步中对应的不同坝体与地基弹性模量组合情况下的大坝在监测点的位移值。

(4)将第(3)步得到的位移值作为输入，对应的坝体及地基弹性模量作为输出，应用非线性模型(如径向基函数神经网络[10]，RBFNN)建立3者的映射关系。

(5)将第(1)步中得到的监测点位移统计模型中的水压分量作为第(4)步建立的非线性模型的输入，即可得到一组坝体与地基弹性模量。可以认为，该组坝体与地基弹性模量即为真实的坝体与地基弹性模量。

(6)确定输入随机变量向量x，一般为上游水位、坝体弹模和坝基弹模等，其中，坝体和坝基弹模均值取为第(5)步的反演值。将各输入随机变量转换为独立的标准正态变量随机向量ξ。

(7)确定重力坝的失效模式，在此基础上建立重力坝的结构功能函数。

表1示，在合并组和单纯CRC组中，Duke分期分布差异有统计学意义，χ2=8.93，P=0.03；经两两比较，A期与B期差异有统计学意义，P=0.002，校准后的检验水准α=0.007。

(8)选用合适形式的Hermite展开多项式构成随机响应面中的响应面函数，以此近似拟合重力坝结构系统的隐式功能函数。

(9)选用合适的配点方法获得随机响应面函数的概率配点。将各配点作为原模型的输入，调用有限元软件计算重力坝有限元模型相应的输出值。

(10)用回归方法计算随机响应面函数中的待定系数，获得随机响应面函数的最终表达式。

(11)针对上一步确定的随机响应面函数表达式，采用几何法或Monte-Carlo模拟法等方法求解重力坝的结构可靠指标与失效概率。

(12)计算结束。

3 工程应用

3.1 工程概况

某水电站挡水建筑物为混凝土重力坝，大坝为2级水工建筑物等级，正常蓄水位为150.00 m。大坝自左岸向右岸分为10个坝段，最大坝高为112.00 m，位于5号坝段。5号坝段最大坝高为112.00 m，坝顶宽度6.0 m，坝体上游面高程84.00 m以上为竖直面，以下坡度为1∶0.3，下游面折坡点高程为145.00 m，折坡点以下坡度为1∶0.75。地基岩性主要为石英砂岩、粉砂质泥质岩石、细砂岩等，上覆第四系松散堆积物，主要有残坡积块及河床冲积漂块石等。

3.2 力学参数反演

建立该坝段的二维有限元分析网格模型，坝基在上下游方向各延伸200 m，坝基深度取为200 m，整体有限元网格及坝体有限元网格分别见图1和图2。模型中一共包括1 955个节点，1 840个单元，且均为四节点四边形单元。x向为大坝上下游方向；y向为大坝高度方向。在地基的上下游侧各节点施加法向约束，在地基的底部各节点施加固定约束。分析中仅考虑上游水压力以及坝体、地基自重荷载的作用。

图1 坝体及坝基有限元网格模型

图2 坝体有限元网格模型 图3 位移测点EX2-4埋设位置

首先，应用式(4)对5号坝段水平位移测点EX2-4的测值进行统计回归分析，得到统计回归方程的相关系数为0.929 3，表明回归模型正确，同时也表明大坝目前基本处于弹性工作状态。当上游水位为141.33 m时，从该统计回归模型中分离出该测点顺河向位移的水压分量为8.435 mm。

然后，按照坝体混凝土型号及坝基岩体试验分别得到的坝体与坝基弹性模量取值范围，将坝体与坝基弹性模量组合为36种工况，通过有限元分析得到各种组合情况下的测点EX2-4处的顺河向位移值，见表1。

表1 坝体与地基弹性模量组合方案

将表1中的位移值向量作为输入，对应的坝体与地基弹模组合作为输出，然后应用径向基神经网络进行非线性拟合，得到3者之间的映射关系。将统计回归模型中分离得到的测点EX2-4顺河向位移的水压分量8.435 mm代入训练后的神经网络，得到坝体弹性模量为21.88 GPa，坝基弹性模量为4.64 GPa。

利用上述反演得到的坝体与坝基弹性模量进行有限元正分析，得到测点EX2-4位置的顺河向位移值为8.367 mm，与统计回归模型中的水压分量的相对误差为0.8%。

综上所述，通过对该水电站5号坝段挡水段的力学参数反演分析，得到该坝段的坝体弹性模量(综合弹模)均值为21.88 GPa，坝基弹性模量为4.64 GPa。

3.3 构建功能函数

在构建功能函数前，首先需确定基本随机变量。结合该水电站重力坝的受力特性，取大坝上游水位Hu、坝体弹性模量Ec以及坝基弹性模量Er3个变量作为结构可靠度分析的基本输入随机变量。大坝上游水位的均值取为正常蓄水位150.00 m，坝体与地基弹性模量的均值取为3.2节中根据监测资料反演分析的结果。各输入随机变量的统计特性见表2。

表2 各输入随机变量的统计特性

以重力坝抗拉可靠度分析为例，重力坝抗拉可靠度的真实功能函数为：

Z=g(x)=min {0.07wb,wc}-wt

(5)

式中：x为重力坝结构可靠度分析的输入随机变量向量；wt为踵的拉应力区宽度；wb为坝底宽度；wc为坝踵至帷幕中心线的距离。

根据式(5)，建立本实例的重力坝结构抗拉可靠度功能函数为：

Z=g(x)=g(Hu,Ec,Er)=6.79-wt

(6)

式中：x为输入随机变量向量；wt为重力坝坝踵的拉应力区宽度，m。

与式(6)意义类似，如果坝踵的拉应力区是连续的，还可以以坝基面上距离坝踵6.79 m处的竖向应力σy的正负来构建重力坝的抗拉功能函数：

Z=g(x)=σy

(7)

式中：σy为坝基面上距离坝踵0.07倍坝底宽度处的竖向应力，应力符号以拉为正，以压为负。

在应用随机响应面法计算结构可靠度之前，首先须将输入的非标准正态分布随机变量向量x=[x1,x2,x3](x1代表大坝上游水位Hu，x2代表坝体弹性模量Ec，x3代表坝基弹性模量Er)转换为标准正态随机向量ξ=[ξ1,ξ2,ξ3]。

采用自编结构可靠度随机响应面法程序，得到二阶随机响应面函数如下：

(8)

3.4 结构可靠度求解

以式(8)作为重力坝结构抗拉可靠度的功能函数，应用几何法求得重力坝抗拉可靠指标为1.629，对应的失效概率为5.17%。根据文献[14]的研究，对于失效概率大于等于10-3的情况，二阶随机响应面法可满足结构可靠度分析的要求，因此，本文采用二阶随机响应面函数作为功能函数是合适的。

标准正态空间的设计验算点为[-0.765, -0.387, -1.385]，将其转换至原始变量空间为[143.115 m, 20.609 GPa, 3.355 GPa]。采用Monte-Carlo抽样并结合MATLAB内置的概率密度拟合函数求得重力坝结构功能函数响应概率密度曲线，见图4。

图4 某重力坝结构抗拉可靠度概率密度曲线

4 结 语

本文提出将实测监测资料反演和非侵入式随机有限元法相结合的结构可靠度分析模式，并给出了具体的算法流程。将该方法应用于工程实例，得到了良好的结果。研究表明：采用实测资料反演力学参数，可以使得输入随机变量的取值更为接近结构真实运行状态；本文提出的算法适用于重力坝等复杂结构的结构可靠度分析。

□