随着结构可靠性评估和稳健性设计概念在现代工程可靠性分析中的日益深入，可靠性分析方法也得到了广泛发展和进一步研究.由于在实际工程问题中其功能函数往往是隐式的(强非线性或耗时的)，采用经典分析技术如一阶可靠度法(FORM)、二阶可靠度法(SORM)^[1]、蒙特卡罗模拟法(MCS)在精度和时间成本上往往难以接受.因此，采用代理模型法(如响应面法^[2]、人工神经网络^[3]、支持向量机^[4]、Kriging^[5]等)来近似功能函数计算可靠性得到了迅速发展.

Kriging代理模型^[6]是一种精确的插值方法且具有随机性，不仅能提供未采样点的预测值，还能对预测方差进行估计.因此，该方法在结构可靠性分析中得到了广泛应用.为了进一步提高模型的计算精度和效率，一些学者对传统的Kriging模型进行了改进研究，如在建立Kriging模型的过程中加入样本点函数值的同时也考虑了其函数值的梯度信息值，这种模型被称为Co-Kriging模型^[7]或梯度增强Kriging(gradient-enhanced Kriging，GEK)^[8]，在此基础上，Han等^[9]为进一步提高GEK的计算效率，提出了改进的梯度增强Kriging方法(weighted gradient-enhanced Kriging, WGEK).Schoebi等^[10]提出了一种基于多项式混沌展开(polynomial-chaos-expansion，PCE)和Kriging模型相结合的可靠性方法(polynomial-chaos-based Kriging，PC-Kriging).在分析比较复杂工程结构问题时，往往采用试验设计的思想.因此，若干种自适应试验设计(DoE)策略已被构建，如Bichon等^[11]通过所提的EFF函数选择距极限状态最近的样本点作为新增训练点.Echard等^[12]提出U函数来衡量未测点符号预测错误的概率，并将U值最小的点定义为最佳样本点.

上述这些基于自适应DoE策略的学习函数能够有效地提高Kriging模型的准确性和效率.然而，由文献[13]可知具有不同基函数的Kriging模型的准确性是不同的，使得Kriging模型在高阶的可靠性分析中难以计算.为了解决这类问题，提出了一种基于PC-Kriging和自适应k-means相结合的结构可靠性分析方法(APC-Kriging).首先，PC-Kriging是一种改进的Kriging算法，其回归基函数采用稀疏多项式最优截断集合来近似数值模型全局行为，而用Kriging来处理模型输出的局部变化，在保证精度的同时提高了计算效率.其次，常见的可靠性方法的采集样本点为逐个采集，而本文的自适应k-means聚类分析将空间分成若干个区域，并从每个区域选取一个最佳样本点，从而使多个区域同时达到提高PC-Kriging模型精度的目的，从而再次提升模型的计算效率.

1 PC-Kriging模型

在构造结构功能函数G(x)的近似代理模型时，Kriging将其假定为由确定性g(x)^Tβ和随机性z(x)两部分组成，即

(1)

其中：向量x为包含M个对结构状态(即功能函数)产生影响的基本随机变量的集合；g(x)为模型基函数，即在PC-Kriging模型中改进的部分；β为基函数的系数矢量；z(x)为均值为0的高斯随机过程，其协方差为

(2)

式中：σ为z(x)的标准差；R(x_i, x_j; θ)为z(x_i)和z(x_j)间的带有参数θ的相关系数，通常为应用最广泛的高斯相关函数.

(3)

式中，θ_m，x_i^m，x_j^m分别为矢量θ，x_i，x_j的第m个元素.

给定N个训练样本S_DoE=[x₁, x₂, …, x_N]及其对应的结构状态值Y=[y₁, y₂, …, y_N]，G(x)预测值的无偏估计及其方差为

(4)

(5)

式中：

其中，参数向量θ需采用极大似然估计获取.

PC-Kriging采用多项式混沌展开替代传统Kriging模型的回归基函数来增强预测模型的全局近似精度，并利用Kriging模型来捕捉预测模型局部特征的能力.采用最小角回归(LAR)构建回归基函数的最优多项式数量集，同时用Akaike信息准则(AIC)来确定最优的截断集合.

令π_j^(m) (j=1, 2, …)表示希尔伯特空间L²(R, f_Xi)的完全正交基，这意味着π_j^(m) (j=1, 2, …)是f_Xi上的正交函数系列.

(6)

式中，f_Xi(x)为x的第i个边缘概率密度函数.

一个完备的Hilbert空间L²(R, f)的标准正交基(f表示x的联合概率密度函数)是

其中，α=[α₁, α₂, …, α_M]是一个自然数的M维向量.这里，满足总项数|α|=α₁+α₂+…+α_M不超过给定阈值T₀的项将被保留作为Kriging模型基函数的候选项.基函数的候选项Α^{M, T₀}可表示为A^{M, T₀}={α∈N^M, |α|≤T₀}，A^{M, T₀}中项数的数量可由P表示:

(7)

然后，考虑基函数的所有候选项的“完备”设计矩阵：

式中：α_i∈A(i=0, 1, …, P-1)；x_n∈S_DoE (n=1, 2, …, N).

根据式(7)，如果将Α^{M, T₀}中的所有候选函数用作Kriging模型的基函数，则函数调用的次数将随着T₀的增加而急剧增加.为了避免这类问题，在构造稀疏多项式Kriging模型的基函数时，采用LAR^[14]理论定义了函数的基函数集个数，并用AIC^[15]准则来确定哪一个是最佳的.最后，保留中间的候选项以确保Kriging模型的准确性，同时大大减少函数调用的数量(即保留对模型贡献更多的候选项).选择Kriging模型基函数的主要步骤：

步骤1  设置相关参数值，即保留多项式T₀的最大阶次和基函数p_max的最多项数，本文设T₀=3，p_max=0.5card(S_DoE).H=min{P, p_max}.

步骤2  初始化所有候选项中的系数a_αi=0 (i=0, 1, …, P-1).根据LAR理论，初始化后的剩余项等于结构响应Y.

步骤3  找出ψ_αi(i=0, 1, …, P-1)与当前残差之间相关系数最大的矢量ψ_α1^′, G₁=ψ_α1^′.

步骤4  h=2.调整a_α1^′在G₁上当前残差的最小二乘系数的方向，直到另一个向量ψ_α2^′与G₁具有相同的相关系数.G_h=[G₁, ψ_α2^′]，a=[a_α1^′, a_α2^′].

步骤5  h=h+1.共同移动a朝向当前残差为G_h-1的联合最小二乘系数，直到向量ψ_αh与当前残差G_k-1具有相同的相关系数.

步骤6  重复步骤5，直到满足h=H.

步骤7  计算G_h (h=1, …, H)的AIC值：

其中, SSE_h=[G_h(G_h^TG_h)^-1G_hY-Y]^T×[G_h(G_h^TG_h)^-1G_hY-Y].

步骤8  找到最小的AIC_h (h=1, …, H), p=arg min_h{AIC_h; h=1, …, H}, 对于S_DoE和Y，PC-Kriging模型的最优基函数为

(8)

2 自适应PC-Kriging方法

构造完近似代理模型后，结构的失效概率可通过式(9)近似计算：

(9)

式中：N_MC为蒙特卡罗样本数；为失效指示函数，当时，；否则，=0.

然而，失效概率估计值的精度取决于与G(x)=0间的“距离”，在的区域，只要与G(x)的符号相同，就不会对的精度产生影响；当与G(x)=0完全重合时，.因此，本文提出自适应策略的思想是在上选取对失效概率贡献大的若干个点，并计算出所选若干个点的结构响应值，在迭代过程中不断更新，进而使逐渐接近G(x)=0直至满足收敛条件.

2.1 自适应PC-Kriging方法

所提出的k-means聚类分析自适应策略选取样本点方法步骤：

步骤1  t=0，最初的试验设计样本点是由拉丁超立方随机采样生成的，并精确地计算出相应的功能函数响应值，即计算，.设初始样本点个数为M₀，则有

步骤2  t=t+1，在上通过马尔科夫链蒙特卡罗模拟法(MCMC)产生K个点.给定，采用MCMC法生成M维服从f(x) (f(x)为x的联合概率密度函数)且满足式(10)的随机向量，则认为随机抽取的点在上.当抽取的点的数量达到K个时随机抽取过程停止，则生成的随机向量为，本文令K为2 000，[ε]为0.01.

(10)

步骤3  采用k-means聚类分析方法将分成k个类别，并将这k类别的中心点映射到上.令{s_{t-1, 1}, s_{t-1, 2}, …, s_{t-1, k}表示k个聚类中心.当为非线性曲面时，则不能保证这k类别的中心点都在上，这时就需要将未在上的中心点映射到上.映射方法为找到满足式(11)的点，并得到，其中为的设计点.

(11)

其中，i=1, 2, …, k.

步骤4  调整集合S_t-1中各点位置.定义距离D₀，如式(12)所示.假设在集合S_t-1中如果任意两个样本点之间的距离小于D₀时可视为是不能接受的，此时需要将集合S_t-1中个别点位置进行调整.

(12)

式中e为给定常数.

集合S_t-1中各点可能出现两种情况：① S_t-1内部某些点间的距离可能过小. 与, …, 中某些点距离很小的可能性较大；② S_t-1中的某个点与X_t-1中某些点间的距离过小.若出现情况①，比如与的距离小于D₀，则要改变概率密度函数较小点的位置，概率密度函数较大点的位置不变；若出现情况②，则改变S_t-1中对应点.样本点位置调整的方法为，假设先要改变的位置，将中点按照与距离升序排列，依次将变换至新序列各点位置，直至满足与所有S_t-1及X_t-1中所有点距离都大于D₀.

步骤5  计算出集合S_t-1中各样本点相对应的功能函数值.Ω_t-1⁰={ , i=0, 1, …, k}，Ω_t=Ω_t-1∪Ω_t-1⁰.

步骤6  根据Ω_t并结合式(4)、式(8)和式(9)计算若满足式(13)的收敛条件，则迭代过程停止，即为P_f估计值；否则返回步骤2，直至满足收敛条件.

2.2 收敛条件

收敛条件采用文献[16]中的学习停止条件，其基本思想为随着迭代过程进行，符号预测错误的样本点数占总失效样本点数的比重很小时，失效概率估计值满足精度要求，学习过程停止，其表达式为

(13)

式中：N_un为符号预测错误的样本总数；N_fail则代表总失效样本数，e′为的许用误差，其中

(14)

式中，用N_{U < P}表明符号预测错误概率高的样本点，这样的点可以看作一定失效的点.用N_P≤U≤Q表明符号预测错误概率较高的样本点，这样的样本点可用N_P≤U≤Q·Φ(-U)表明其失效预测错误的总个数.本研究中，P=1，Q=2，e′=0.03.

3 算例 3.1 算例1

选取文献[12]中具有多个设计点的状态函数，其表达式为

其中，x₁, x₂为服从标准正态N(0, 1)的独立同分布随机变量.

适当抽取构建初始PC-Kriging模型时所需样本点数，设N₀=6.为取得分布较为均匀的初始样本点，采用拉丁超立方抽样(Latin hypercube sampling, LHS)在[-5, 5]区域内抽取初始样本点.通过Matlab中的工具箱根据所提算法建立PC-Kriging预测模型，并通过主动学习，更新样本空间DoE，并不断提高模型精度，重复该过程，直到满足迭代停止条件，结果见图 1、表 1.其中，N_it表示迭代次数，N_call为调用样本点数，ε为失效概率估计值与标准值的相对误差，失效概率标准值由MCS方法计算得到.

通过图 1和表 1的结果对比，可以发现在N_call项所提方法需要最少的样本点就可以达到足够高精度，而通过对N_it项，因为应用到聚类分析，使得在满足精度的要求下迭代次数得到大幅度的降低，从而提高计算效率.

3.2 算例2

选用文献[17]悬臂式圆柱筒结构(图 2)，它是一个拥有9个随机变量的高维非线性的工程结构实例，该结构的9个输入随机变量分别为t，d，L₁，L₂，F₁，F₂，P，T，S.其分布特征如表 2所示.

图 2中悬臂式圆筒结构受到外力F₁，F₂，P和扭矩T的作用，其功能函数表示为屈服强度S和最大应力σ_max的差：

其中：σ_max表示在原点处筒上表面所受的最大等效应力；σ_x为正应力；τ_zx为扭应力.

其中:θ₁=5°; θ₂=10°，M为弯矩，M=F₁L₁cosθ₁+F₂L₂cosθ₂；；c=d/2；

首先，采用Nataf变换将上述随机变量映射到标准正态空间，并在[-5, 5]⁹立方体内抽取N₀=11个拉丁超立方随机样本点，再应用所提算法评估悬臂式圆柱筒的可靠性.所提算法与其他现有算法所得失效概率估计值随迭代次数的变化趋势如图 3所示，其中横坐标为迭代次数N_it，纵坐标为失效概率估计值的对数形式.此外，表 3列举了本文算法与其他方法所得结果.通过比较，不难发现本文算法相较于其他方法在效率和精度方面更具优势.

4 结论

1) 本文提出了一种新的改进Kriging基函数与k-means聚类分析相结合的可靠性分析方法.使得PC-Kriging方法在基函数回归方面具有更大的灵活性.结合k-means聚类分析方法的自适应策略在每次迭代时添加多个新增样本点来提高计算效率.

2) 通过与两个算例中MCS抽样计算的结果对比可知，所提出的算法与MCS模拟较为接近，从而验证了本文所提算法计算精度的正确性.

3) 通过对比的可靠性分析结果可以发现，所提方法相对于AK-MCS+U和AK-MCS+EFF将迭代次数降低同时提高了估计精度，说明所提方法具有更快的效率以及更高的精度.

4) 所提方法亦可适用于解决功能函数为隐式多维非线性问题，这为解决实际工程中的可靠性评估问题提供了重要的参考价值.