土木工程结构在实际使用中承受的外部作用、构件截面尺寸和材料性能等一般都有随机性，因此在结构的设计中，有必要考虑结构本身特性及所受外部作用的随机性。复杂结构在地震作用下的可靠性分析比较困难。对于复杂或新型结构，常不能用简化模型来计算响应，需要使用有限元法进行分析，功能函数不能显式表达且有时是高非线性的，导致可靠度难以由矩方法计算。通过Monte Carlo法(MC)计算结构可靠度不受复杂功能函数的限制，且结果准确，但需要对结构响应进行很多次计算。大型复杂结构的有限元分析一般比较耗时，例如巨-子型有控结构(mega-sub controlled structural system，MSCSS)，由于其振动控制原理与常规结构不同，响应难以用反应谱法或静力弹塑性分析法计算，而需要由耗时的动力时程分析得到。成千上万次分析会耗费大量的计算时间，这限制了MC法在实际工程中的应用。代理模型是目前受到关注的一种有效提高计算效率的方法，该方法使用运算量小且误差在可接受范围内的近似函数代替运算量大的函数来进行计算，从而使可靠度分析效率大幅提高。

常见的代理模型有响应面法、径向基函数、Kriging模型、支持向量机等。Kriging模型是由南非地质学者Krige提出并由法国数学家Matheron完善而形成的一种插值方法^[1]，它具有较强的全局近似和非线性拟合能力。近年来，人们对Kriging模型在可靠性分析中的应用进行了大量研究。Kaymaz对比了Kriging模型和响应面在可靠度计算中的使用效果，认为Kriging模型存在参数难以选取等问题^[2]。Bichon等参考EGO方法提出高效全局可靠性分析(EGRA)方法，该方法将EF函数作为目标函数^[3]，使用全局优化算法寻找加点位置^[4]。Echard等提出AK-MCS方法，其结合U函数和MC法来选取Kriging模型的训练样本，该方法由于计算精度高且不需要借助优化算法^[5]，引起了人们的关注，很多改进的AK-MCS方法被提出^[7]。Jia等将Kriging模型用于框架的抗震可靠性分析，验证了基于Kriging的可靠度计算法在涉及时程分析的工程问题中的有效性^[8]。Zhou等结合活跃子空间法提高了AK-MCS在高维问题中的适用性，并将该方法用于非线性框架的抗震可靠性分析^[9]。

目前，已有多种基于Kriging模型的可靠度计算方法，其中，EGRA和AK-MCS是2种经典的基于Kriging模型和数字模拟的可靠性分析方法，计算精度较高。二者均采用在初次采样基础上增加训练样本点的方式来保证代理模型的精度，由于加点停止条件较保守，用于土木工程结构时，可能会产生不必要的计算成本。根据土木工程结构的特点，本文提出一种基于Kriging模型的可靠性分析方法，其训练样本加点停止条件的工程含义直观，便于设计人员根据需要进行调整。通过2个算例验证了该方法的准确性和高效性，并结合时程分析对巨-子型有控结构体系的抗震可靠度进行了计算，分析结果表明本文方法适用于大型复杂结构体系的可靠度计算分析。

1 Kriging模型

Kriging是一种插值近似方法, 由Kriging模型不仅可以预测未知函数y(x)在某点x的值, 还可以给出该点的预测均方误差(mean squared error, MSE)σ_y²(x), 且函数值y(x)满足正态分布

(1)

式中, μ_y(x)为预测值的均值。

已知m个样本点X=[x⁽¹⁾, x⁽²⁾, …, x^(m)]及相应函数响应值Y=(y₁, y₂, …, y_m)^T, 点x的Kriging预测值为样本点函数响应值的线性加权

(2)

式中：x=(x₁, x₂, …, x_n)为n维向量；c_i(x)(i=1, 2, …, m)为加权系数。

Kriging模型假设未知函数可表示为

(3)

(3) 式由两部分组成, 第一部分为回归多项式, p为多项式的项数, f_j(x)为回归多项式的基函数, β_j为回归系数; 第二部分Z(x)是均值为0、方差为σ_Z²的随机过程, 协方差为

(4)

式中: R(x, w; θ)为任意两点x和w的相关函数；θ=(θ₁, θ₂, …, θ_n)^T为其参数, 相关函数通常采用高斯型, 可表示为

(5)

通过在满足Kriging预测值无偏性条件下最小化预测均方误差, 可以得到未知点x处的预测均值和均方误差为

(6)

(7)

式中

求解优化问题

使lnL(θ)最大化, 可得到θ的最优取值。

2 基于Kriging模型的可靠性分析方法 2.1 结合Kriging和Monte Carlo模拟的可靠性分析

根据时程分析中结构最大响应是否超过限值来构造功能函数的方法^[9]便于工程设计人员使用, 而且在时程分析中不仅可以考虑结构的随机性, 还可以通过调整地面运动时程的参数来反映地震作用的随机性。本文采用这种方法来构造功能函数。不过, 对于复杂结构, 会得到一个计算耗时的隐式功能函数。为能简便、高效地计算结构可靠度且保证结果的准确性, 本文提出一种基于Kriging模型的可靠度计算方法, 该方法的基本思路是根据功能函数在少量训练点的值建立其Kriging代理模型, 并在极限状态曲面附近不断改善代理模型的精度直到满足工程需要, 用代理模型计算功能函数值从而避免直接计算复杂的真实功能函数, 然后结合MC法计算失效概率, 以下对该方法进行详细阐述。

根据功能函数G(x)是否大于0可将变量空间分为安全域和失效域, G(x)=0表示极限状态曲面。MC法以大数定律为基础, 随机抽取N个样本点, 以落入失效域的样本数N_F与N的比值作为真实失效概率的估计值。MC法鲁棒性好, 且易于编程实现, 但是需要对G(x)值进行多次计算。如果根据少量训练样本点建立G(x)的Kriging代理模型, 得到近似功能函数, 用替代G(x)来计算随机样本点的功能函数值, 则可得到失效概率近似值

(8)

式中, 为N_F的近似值, 当时, 当时。这样, 计算可靠度需要调用真实功能函数的次数仅取决于训练样本点的数量。

近似失效概率的准确程度取决于近似函数的精度, 该精度与训练样本的选取有很大关系。在选取训练点时, 通过实验设计方法(如超拉丁方抽样)可使各点均匀地分布于变量空间中, 训练点从而更具有代表性, 能够建立精度更高的代理模型。但即使如此, 一次抽样选取的训练点也往往不能保证Kriging模型满足可靠性分析的需要, 还需要加点采样。由(8)式可知, 提高可靠度计算精度的关键是提高失效域中样本点数量近似值与真实值N_F的接近程度。由于极限状态曲面是安全域与失效域的分界面, 近似极限状态曲面与真实极限状态曲面G(x)=0越接近, 则与N_F越接近。因此, 代理模型在极限状态曲面处的精度对可靠度计算的准确性有十分重要的影响。

本文通过在极限状态曲面附近不断增加训练点来提高代理模型在曲面处的精度。使用全局优化算法求解(9)式即可在近似极限状态曲面上搜索预测均方误差最大点x^*

(9)

x_lb, x_ub为采样范围的界限, 若在x^*处的模型精度未达到加点停止条件, 则将该点作为新增样本点加入到训练样本中以提高Kriging模型的精度。通过这种方式每次增加一个训练点, 直到代理模型在极限状态曲面处的精度达到加点停止条件。(9)式的求解可使用Matlab遗传算法工具箱实现。

由于失效概率计算结果的精度取决于近似极限状态曲面与真实曲面的接近程度, 当近似曲面上各点的预测值与实际值的误差小于一定限值ε时, 可认为代理模型满足工程计算的要求。根据(1)式, 可知, 若某点预测标准差σ_G(x)≤ε/2, 则该点功能函数预测值与真实值间误差小于ε的置信度可达95.4%以上。因此, 为使近似曲面上各点预测值的误差小于ε, 需要使曲面上预测误差最大点x^*的MSE小于ε²/4。另外, x^*处的实际误差也应小于ε, 即|≤ε。因而, 以在x^*处同时满足σ_G²(x^*)≤ε²/4和作为加点停止条件。误差限ε可取小于初次采样点功能函数值标准差1/50的数, 也可根据工程问题的精度需要进行调整。

综上, 本文所提可靠性分析方法的步骤为:

步骤1  用实验设计方法取少量训练样本点, 并计算各训练点的功能函数值。采样范围应能让大部分的Monte Carlo抽样点落入其中, 如[μ_x-5σ_x, μ_x+5σ_x]

步骤2  利用训练样本及其函数值建立Kriging模型

步骤3  求解(9)式中的优化问题得到最优点x^*

步骤4  求x^*处的功能函数值, 若满足且, 则停止加点; 否则, 将x^*作为新增样本点加入训练点集合, 返回步骤2

步骤5  通过式(8)计算失效概率近似值。

2.2 算例验证

算例1

一串联系统的功能函数为

x，y相互独立, 服从标准正态分布。

该算例选择文献[5], 极限状态曲面较复杂, 具有多个设计点。分别使用MC法、一阶可靠度法(FORM)、EGRA法、AK-MCS法和本文方法计算了该算例的失效概率, 结果列于表 1中。其中, N_call表示调用真实功能函数的次数, 以MC法计算结果为失效概率的参考值, ε_f表示失效概率计算结果相对于参考值的误差, 。

使用本文方法进行了2次可靠度计算, 2次分析的初始训练点分别为11个和17个, 调用真实功能函数的总次数分别为64次和46次, 均小于表中其他各种方法的调用次数, 计算的结果却与参考值十分接近, 而由FORM得到的结果与参考值偏差很大。

图 1给出了加点过程中的变化(初始训练点为11个), 由于功能函数是二维的, 极限状态曲面实际上是曲线。可见, 未加点时, 近似曲线与实际曲线相差较大, 随着训练点的增加, 近似曲线逐渐接近实际曲线。达到加点停止条件后, MC样本点的分布如图 2所示, 近似曲线上随机样本分布较多的部位, 已与真实极限状态曲线接近重合。

算例2

图 3所示非线性系统^[7]的功能函数为

式中, 功能函数含6个相互独立的随机变量, 均为正态分布, 其分布参数见表 2, 这是一个中等维度的可靠度计算问题。

由各方法得到的可靠度计算结果列于表 3中。其中, 使用本文方法进行了2次可靠度计算, 二者初始样本点分别为11和17个, 计算结果与参考值相差很小, 而调用真实功能函数的次数仍是表中各方法中最少的。加点过程中(初始点为17个)失效概率计算值随训练点增加的变化如图 4所示, 可见, 未加点时计算值与参考值偏差较大, 随着样本点数的增加, 计算结果的精度会得到改善。

通过2个算例可以看出, 使用本文提出的方法可得到较准确的可靠度分析结果, 且调用真实功能函数的次数少于其他方法, 该方法的准确性和高效性得到了验证。另外, 尽管算例中随机变量是相互独立的, 但本文方法并不要求变量相互独立, 因为在可靠性分析中代理模型仅被用于求功能函数的近似值且建模过程不涉及概率密度函数或其相关性的信息。

3 MSCSS的抗震可靠性分析

MSCSS是结合巨型框架结构(mega frame structure, MFS)和调频子结构控制原理形成的一种新结构体系。该结构在MFS基础上放松主结构与子结构之间的部分连接, 并在主、子结构之间设置阻尼装置, 通过主、子结构的相对运动来减小结构在地震作用下的响应, 从而提高了抗震性能^[10]。根据前文比较分析的结果, 将所提方法用于MSCSS的抗震优化中以验证其在工程结构中的适用性。

图 5为MSCSS的有限元模型, 共32层, 层高4 m, 含有3个巨层。巨型梁采用截面为4 m×4 m×0.04 m×0.04 m的工字钢, 巨型柱采用壁厚0.05 m的方钢, 巨型柱的中心间距为26 m, 子结构柱截面为0.6 m×0.6 m×0.02 m的方钢, 子结构梁截面为0.5 m×0.2 m×0.02 m×0.02 m的工字钢。2, 3巨层为调频子结构, 其5、10层设置与巨型柱相连的黏滞阻尼器, 阻尼系数c_d为1.5×10⁶ N/(m·s^-1), 阻尼指数e为0.4, 阻尼器的非线性力f_d=c_d·v^e, v为阻尼器变形速率。抗震设防烈度为9°。

本文以结构在设计基准期内由于地震作用层间位移角ϕ超过变形限值ϕ_b为失效破坏准则。参考《建筑抗震设计规范》(GB50011-2010), 变形限值ϕ_b取1/250。结构的变形值通过SAP2000中非线性动力时程分析得到。可靠性分析中考虑地面最大加速度、钢材模量和巨型柱截面宽度的随机性。随机变量的分布参数如表 4所示, 其中, 地面运动最大加速度a服从极值Ⅱ型分布^[12], 其余服从正态分布。结构的功能函数可表示为为

式中: 变量x=(a, E, b₁, b₂, b₃); ϕ_m(x)为地震过程中结构的最大层间位移角。在Matlab中通过SAP2000提供的接口可调用该软件进行时程分析并获得其结果数据以计算ϕ_m(x), 从而得到G(x)的值。

使用本文方法计算结构的失效概率, 初始训练样本点数为11个, 相应功能函数值的方差为3.61×10^-3, 误差限ε取相对其较小的数, 为3×10^-5。增加了41个训练点后, 代理模型达到精度要求, 可靠度计算结果列于表 5中。为验证该方法在工程问题分析中的准确性, 使用AK-MCS对结构的可靠度进行了计算以作对比。使用本文方法得到的可靠度结果与AK-MCS法的结果十分接近, 而调用真实功能函数的次数却少很多, 可见本文方法适合用于复杂建筑结构的可靠性分析中, 可以高效地得到较准确的计算结果。

4 结论

提出一种基于Kriging模型的可靠性分析方法，2个数值算例的计算结果表明，该方法能够保证可靠度计算的精度且计算效率较高。将该方法应用于MSCSS的抗震可靠度计算中，使用较少的有限元分析次数得到了准确的分析结果，从而验证了该方法在工程问题中的适用性和高效性。该方法的加点停止条件工程含义直观，便于工程应用，可避免过多的仿真分析，是一种适用于大型复杂结构可靠性分析的有效方法。