随着机械结构的日趋复杂，基于代理模型(如响应面RSM^[1]、神经网络NNs^[2]、支持向量机SVM^[3]、Kriging^[4-7]等)的可靠性分析方法由于其高效性和准确性而越来越受到重视.Kriging模型作为一种精确的内插方法，在提供预测响应的同时还给出了预测响应的方差，这使得Kriging模型在可靠性分析中得到了广泛应用并逐渐成为研究热点.

近年来，为改善Kriging模型的精度和可靠性分析效率，众多学者提出了各种自适应试验设计策略和可靠性分析方法，其中，基于学习函数的选点策略是研究最广泛的.例如，Bichon等^[8]提出了学习函数EFF，用以测量样本接近极限状态的程度；Echard等^[9]定义了学习函数U，用以选择最容易错分类的样本；Sun等^[4]通过提出的学习函数LIF来选择对失效概率预测精度改善最大的样本.此外，学者们还从采样区域^[10]及交叉验证^[11]等方面发展了各种选点策略.然而，现有自适应试验设计的停止准则主要聚焦于单个样本满足学习函数的要求，而不是从预测的失效概率满足可靠性分析精度需求的角度出发，这往往会增加不必要的功能函数调用次数，降低可靠性分析效率.为此，Hu等^[12]和Wang等^[13]基于失效概率的相对误差估计提出了两种新的停止准则；然而，由于保守或错误地估计真实失效样本数的取值范围，新的停止准则仍难以应用到工程实际中.

针对上述存在的问题，本文基于失效概率的相对误差估计提出了一种更精确的真实失效样本数的区间估计方法，从而有效地避免了不必要的功能函数调用和错误的失效概率评估结果，并采用两个实例来说明所提方法的通用性、准确性及高效性.

1 可靠性分析基本理论

在可靠性领域中，失效概率被定义为

(1)

式中：x为N维随机向量(即N个影响结构输出的随机输入因素)；f(x)为x的联合概率密度函数；G(x)为结构的功能函数且G(x)将x的空间分为两部分，即失效域{x|G(x)≤0}和安全域{x|G(x)>0}.

在Kriging模型中，G(x)被表示为

(2)

式中：g(x)为多项式回归基函数；β为相应的回归系数，本文令g(x)的阶次为0；z(x)为高斯过程.z(x)具有如下特征：

式中：σ是高斯过程z(x)的标准差；R(x, w; θ)表示z(x)和z(w)间含未知参数θ的相关系数.本文采用应用最广泛的高斯相关函数：

(3)

式中：x^k，w^k，θ^k分别表示x，w和θ的第k个元素.

在给定含M个训练样本的初始DoE(design of experiments)下，β和σ的估计分别为

(5)

式中：g=[g(x₁), g(x₂), …, g(x_M)]^T, R=[R(x, w; θ)]_N×N.

参数θ由极大似然估计得到：

(6)

为确保预测响应最小方差无偏性，样本x处的Kriging预测响应及其方差分别为

(7)

(8)

式中：

因此，结构失效概率的估计值为

(9)

对复杂结构而言，直接求解积分(9)很困难，甚至是不可能的.因此，本文使用Monte Carlo(MC)法近似估计，即

(10)

式中：N_MC为MC样本数；x_i是第i个MC样本；(x)为失效指示函数，当(x)≤0时，(x)=1，当(x)>0时，(x)=0.

的变异系数为

(11)

的相对误差为

(12)

式中，为通过Monte Carlo Simulation(MCS)和真实功能函数得到的结构失效概率.

2 失效概率的相对误差估计 2.1 学习函数U

Echard等^[9]从样本符号预测精度的角度出发定义了一个学习函数U来量化符号预测的不确定性，其表达式为

(13)

则该样本x处符号预测错误的概率为

(14)

式中, Φ(·)为标准正态分布的累积分布函数.

显然，μ_G(x)越小(样本x越接近极限状态G(x)=0)或σ_G(x)越大(样本x的不确定性越高)或两者兼之，则U(x)值越小，样本x处符号越容易预测错误.因此，定义U(x)最小的点为新训练样本，即

(15)

基于函数U的自适应可靠性分析流程如下：

① 在[-5, 5]^N空间内采用拉丁超立方抽样(LHS)定义初始训练样本集Ω，并应用MC法产生候补样本集Ω_CS；

② 基于样本集Ω，构建Kriging代理模型，并根据式(7)、式(8)和式(13)找到样本集Ω_CS中对应的最小U值min(U(x))，如果min(U(x))>2，停止迭代，由式(10)~式(12)计算可靠性结果，否则继续步骤③；

③ 计算新训练样本x^*处的功能函数值G^*，将其加入到初始训练样本集Ω中，并返回步骤②.

2.2 相对误差估计方法

根据学习函数U，U(x)>2的样本的符号可以认为是预测正确的.因此，式(12)中的相对误差可以改写为

(16)

式中：分别表示域和Ω_f2内的样本总数，且域分别代表如下所示的样本点.

显然，N_f2等于域内符号预测正确的样本数和域内符号预测错误的样本数之和.类似于失效指示函数，本文定义一个指示函数I_w(x)，用于识别样本符号是否预测错误，即

因此，域内符号预测正确的样本数S_fc和域内符号预测错误的样本数S_sw可分别表示为

(17)

式中, 为域内的样本总数.

忽略Kriging预测之间的相关性，则S_sw和S_fc均服从泊松二项分布，其分布参数为

(18)

(19)

式中：U_i=U(x_i); U_j=U(x_j).

随着Kriging模型精度的逐渐改善，U(x)≤2的样本数逐渐减少，文献[13]采用的林德伯格-费勒中心极限定理不再适用.此外，根据Le cam定理定义的误差界限可知，只有在样本符号预测正确或错误的概率较小时，泊松分布才能较好地近似泊松二项分布.因此，本文采用适用性更广泛且近似效果优于泊松分布的正态分布来模拟S_sw和S_fc的分布，即

(20)

因此，N_f2的分布为

(21)

式中：μ_Nf2=μ_sw+μ_fc-1；σ²_Nf2=σ_sw²+σ_fc².

在置信水平为0.95时，N_f2的置信区间为

(22)

文献[12]将N_f2的取值区间定义为

(23)

显然，式(23)是极端情况下(域内的样本符号全部预测错误且域内样本符号全部预测正确，即N_f2取值为0；域内的符号全部预测正确且域内符号全部预测错误，即N_f2取值为+ )得到的N_f2的取值区间，因此过于保守.

文献[13]将N_f2的取值区间定义为

(24)

显然，式(24)中的取值下限忽略了域内符号预测错误的样本数，取值上限则认为域内样本符号全部预测正确，因此式(24)中N_f2的取值区间不精确.

此外，在Kriging模型精度较差时，由于域，内样本(即U(x)≤2的样本)总数较多，导致式(22)和式(24)中N_f2的下限过高，即U(x)≤2区域内的真实失效样本数没有落在N_f2的取值区间内，从而错误地估计相对误差和失效概率.

随着Kriging模型精度逐渐改善，域内的样本数基本呈增长趋势且逐渐收敛于整个样本空间内真实的失效样本数N_f，即N_f2与间的比值逐渐减小.因此，为使U(x)≤2内的真实失效样本数落在N_f2的取值区间内，本文提出引入一个比例因子α，通过α 重新定义N_f2的取值下限，即

(25)

本文令α=0.1，则失效概率估计值的相对误差为

(26)

因此，只要式(26)中失效概率的最大相对误差ε_max满足式(27)中的条件，就说明满足可靠性分析的精度要求.

(27)

式中, ε_thr为最大相对误差的阈值，工程实际中ε_thr一般为0.01~0.05^[12-13]，本文取ε_thr=0.03.

2.3 基于相对误差估计的可靠性分析

基于相对误差估计的自适应可靠性分析方法的基本步骤如下：

步骤1   采用LHS法产生初始训练样本集S_DoE，并计算相应的真实功能函数值G_DoE.初始训练集S_DoE一般只需少量的样本，在后续迭代选点过程中不断更新训练集.

步骤2   基于现有DoE和MATLAB中的DACE工具箱构建Kriging模型，并应用MC法产生候补样本集Ω_CS.

步骤3   计算候补样本集Ω_CS相应的μ_G(x)，σ_G(x)和U(x)；根据式(26)获取失效概率估计值的最大相对误差ε_max，如果不满足式(27)，将样本集Ω_CS中最小U值min(U(x))对应的样本作为新训练样本加入到DoE中并返回步骤2；否则，继续下一步.

步骤4   应用MC法估计失效概率及其变异系数，如果变异系数满足式(28)，则停止迭代并输出可靠性分析结果；否则，扩大MC样本数并重复步骤2~4，直至满足式(28).

(28)

3 算例验证 3.1 算例1

选取文献[8]中的二元非线性功能函数作为研究对象，其表达式为

(29)

式中，x₁~N(0, 1)，x₂~N(0, 1)，且相互独立.

初始LHS训练样本点设置为6个，并采用10⁶个MC样本点计算该二元函数的失效概率.由不同真实失效样本数N_f取值范围的估计方法及相应停止条件得到的可靠性分析结果如表 1所示，其中N_call为功能函数调用次数.显然，本文的可靠性分析方法在满足精度要求的同时，比其他方法需要更少的功能函数调用次数.

图 1、图 2分别显示了由不同估计方法得到的N_f和失效概率相对误差ε_f随迭代次数N_it的变化.可以发现，本文所估计的N_f的取值范围和失效概率的最大相对误差比文献[12]和文献[13]的更精确.

3.2 算例2

以图 3所示的具有9个随机变量的悬臂式管结构^[14]为分析对象，其功能函数为结构屈服强度S和最大应力σ_max之差，即

(30)

式中各随机变量的分布特征如表 2所示，σ_max为原点处悬臂管上表面所受的最大等效应力：

(31)

式中

初始LHS训练样本点设置为12个，并采用2.0×10⁵个MC样本点计算该悬臂管的失效概率.表 3和图 4列举了采用不同方法得到的结果.相较于其他方法，在满足可靠性分析精度要求时，本文方法在效率上更具优势.

4 结论

1) 通过两个算例表明，本文确定的取值区间不仅能够准确地包含真实失效样本数，而且还能得到更精确的失效概率相对误差估计结果.

2) 本文提出的基于失效概率相对误差估计的自适应可靠性分析方法在满足分析精度要求的同时极大地减少了功能函数的调用次数，为解决实际工程中复杂结构的可靠性问题提供了一种切实可行的分析方法.