刘阔¹, 李晓雷², 王健³
1. 大连理工大学 机械工程学院, 辽宁 大连 116024;
2. 哈尔滨工业大学 机电工程学院, 黑龙江 哈尔滨 150001;
3. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2016-01-20
基金项目：中央高校基本科研业务费专项资金资助项目(DUT16RC(3)122)。
作者简介：刘阔(1983-), 男, 河北石家庄人, 大连理工大学副教授。

摘要：通过分析现有机械结构可靠性抽样方法存在的不足以及影响失效概率估计精度的主要因素, 提出了一种基于Kriging模型及自适应抽样方法的机械结构可靠性分析方法.该抽样方法将随机抽样与聚类算法相结合, 能够在概率上保证新增样本点落在对失效概率贡献较大的区域, 避免对非重要区域的不必要抽样.以大数定律及中心极限定理为基础, 推导了所提出的Kriging模型的收敛条件.通过两个算例说明所提出方法的迭代收敛过程、准确性及稳定性, 结果表明, 该方法能够在较少调用结构功能函数情况下得到失效概率较准确的估计值.
关键词：Kriging模型机械结构可靠性失效概率自适应抽样方法蒙特卡罗方法
An Analysis Method of Mechanical Structural Reliability Based on the Kriging Model
LIU Kuo¹, LI Xiao-lei², WANG Jian³
1. School of Mechanical Engineering, Dalian University of Technology, Dalian 116024, China;
2. School of Mechatronics Engineering, Harbin Institute of Technology, Harbin 150001, China;
3. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: LIU Kuo, E-mail:liukuo@dlut.edu.cn
Abstract: Based on an analysis of the drawbacks of the existing mechanical structure reliability sampling methods and the main factors influencing the estimation accuracy of failure probability, an analysis method of mechanical structural reliability based on the Kriging model and adaptive sampling strategy is proposed. The proposed sampling strategy combines random sampling and clustering algorithm, and ensures in probability that the new sample points locate themselves in the area that makes significant contribution to failure probability and avoids unnecessary sampling in the unimportant areas. The condition of convergence for the proposed Kriging model is deduced mainly based on the law of large numbers and central limit theorem. Two examples are adopted to illustrate the convergence process, accuracy and stability of the proposed method. The results show that the proposed method can estimate failure probability with high accuracy in the condition that the number of calls to structural performance function is small.
Key Words: Kriging modelmechanical structural reliabilityfailure probabilityadaptive sampling strategyMonte Carlo method
机械结构可靠性分析的主要目的是在考虑结构各输入参数随机性情况下计算其失效概率.近年来国内外学者将一阶矩、二阶矩、Monte Carlo(MC)^[1]、代理模型等方法应用于机械结构可靠性分析.一阶矩及二阶矩方法难以适应工程结构非线性隐式的特点, 精度有限.MC方法鲁棒性最好, 然而耗时过长.代理模型方法是目前应用最广泛的机械结构可靠性分析方法.
众多学者通过调整代理模型函数形式^[2-6]及选取样本点策略^[7-10]提高代理模型精度及效率.Blatman^[2], Roussouly^[3]等将稀疏多项式应用于结构可靠性分析, 通过不同方法筛选“重要项”.此外, 为提高模型的非线性程度, Kriging^[4-5]、神经网络^[6]、SVM^[7]也被广泛采用.Kriging模型是一种高效的差值模型.因此, 本文采用Kriging作为基本代理模型.
样本点的选取策略对代理模型精度和效率有很大影响.为提高样本点的质量, 很多自适应抽样方法被提出^{[3, 5, 7-10]}.现有抽样方法中鲜有考虑在迭代抽样过程中可能出现由于个别样本点间距离过近甚至矩阵奇异问题, 尤其是当算法接近收敛时, 设计点附近很小范围内往往集中过多样本点.本文提出一种新的自适应抽样方法，避免对“不重要区域”没必要抽样; 保证抽取的样本点保持与现有样本点的距离.
1 主要理论假设某机械结构的输入变量为X(X为N维随机向量), X的联合概率密度函数为f(x), 该结构的功能函数为G(x), 其安全域和失效域分别为S_s={x|G(x)>0, x∈R^N}和S_f={x|G(x)≤0, x∈R^N}, 则该结构失效概率为

(1)

若已知样本集Ω={(x_i, y_i), i=1, 2, …, M}, 且功能函数G(x)可表示为

(2)

其中：g(x)为多项式基函数；β为g(x)的系数向量, 本文中取g(x)的次数为1；z(x)为零均值同方差高斯过程, z(x_i), z(x_j)的协方差为

(3)

其中：σ_z²是高斯过程的方差；R(x_i, x_j; θ)表示z(x_i)和z(x_j)的相关系数, θ为参数.应用最广泛的是高斯相关函数, 其形式为

(4)

其中, x_{i, n}表示向量x_i的第n个元素.
在给定样本集Ω情况下, 式(2)~式(4) 中的未知参数β, σ_z², θ可通过极大似然法估计得到,

(5)

式中:为广义最小二乘法估计的残差；R为M×M相关系数矩阵；分别表示θ, β, σ_z²的极大似然估计值.点x处结构的功能函数可表示为Y的线性组合形式,

(6)

式中，Y=[y₁, y₂, …, y_M]^T.
要实现的最小方差无偏性^[11], 则c(x)需满足

(7)

求解式(7), 可得的Kriging表达式为

(8)

则该结构失效概率P_f的估计值为

(9)

当X的维数较高时, 直接求式(9) 的积分值很困难.本文采用MC法通过随机抽样近似计算, 根据大数定律,

(10)

其中：N_MC为随机抽样次数；x_{MC, i}(i=1, 2, …, N_MC)是来自f(x)的独立同分布随机序列.由中心极限定理知, 当N_MC足够大时, 的均值依分布收敛于正态分布, 可得
其中:Ф(·)为标准正态分布函数; λ₀为正数.
若要保证在95%的置信水平下与间的相对误差不超过5%, N_MC需满足

(11)

通常情况下工程机械结构的失效概率数量级很小, (1－)应接近1, 因此忽略分母中的(1－), 式(11) 可写成

(12)

N_MC表示当结构的极限状态函数为=0时N_MC次随机抽样中失效样本个数的期望.
2 本文的方法及收敛条件2.1 本文的方法本文提出在=0上选取具有“代表性的”点, 计算这些点的结构功能函数值, 更新=0, 在迭代过程中使=0逐渐接近G(x)=0，直至满足收敛条件.
本文提出的选取样本点方法的主要步骤为
步骤1 ?t=0, 应用拉丁超立方抽样方法^[12]随机产生初始样本点并计算各点结构功能函数值, 计算初始样本点个数为M₀, 令
步骤2 ?t=t+1, 随机抽取K个=0上的点.给定, 生成N维服从f(x)的随机向量, 若该随机向量满足式(13), 则认为其在=0上.当满足式(13) 随机向量数达到K个时, 停止生成随机向量过程, 并记

(13)

步骤3 ?对作聚类分析, 将其分成k类并得到k个类别的中心, 并将k个中心点映射到=0上.令{s_{t－1, 1}, s_{t－1, 2}, …, s_{t－1, k}}表示k个聚类中心.通常情况下=0是非线性曲面, 因此不能保证k个聚类中心点在=0上, 因此要将其映射到=0上.映射的方法是在=0上寻找满足式(14) 的点, 得到, 其中为的设计点.

(14)

其中, i=1, 2, …, k.
步骤4 ?调整集合S_t-1中各点位置.如式(15) 定义D₀, 本文中认为两个样本点间距离小于D₀时是不可接受的, 需要调整个别点位置.

(15)

式中，e₀为给定常数.
步骤5 ?计算S_t-1中各样本点的结构功能函数值.令
步骤6 ?根据Ω_t计算, .若满足式(16) 的收敛条件, 则停止迭代过程, 即为P_f估计值; 否则返回步骤2, 直至满足收敛条件.
2.2 收敛条件当迭代过程收敛时, 可以认为
由中心极限定理, 若能保证抽样次数满足
则
根据式(11) 和式(12),
收敛条件为

(16)

3 实例分析3.1 实例1选用文献[10]中某二维输入变量结构作为分析对象, 该算例的主要目的是进一步用于说明本文所提出的方法.该结构的功能函数为
其中, x₁, x₂相互独立且均服从标准正态分布.

(17)

其中,
通过MATLAB软件, 采用数值积分的方法计算式(17), 得到失效概率为
采用拉丁超立方抽样方法在[-5, 5]×[-5, 5]内随机产生6个初始样本点(M₀=6), 计算其功能函数值, 得到Ω₀, 建立初始Kriging代理模型, 如图 1a所示.取K=200, ε=0.01, 随机产生K个=0上的点.取e₀=0.15, k=2, 依据2.1节步骤3的方法, 对K个随机点进行聚类分析, 并将k个聚类中心映射到=0上.图 1为算法迭代过程中G(x)=0与=0的对比及每步迭代产生的新样本点.算法在第4次迭代时收敛, 计算得到=3.59×10^－3, 图 2为收敛后G(x)=0与=0的对比.
图 1(Fig. 1)

图 1 迭代过程中G(x)=0与=0对比Fig.1 Comparison between G(x)=0 and =0 during iteration (a)—第1次迭代；(b)—第2次迭代; (c)—第3次迭代；(d)—第4次迭代.

图 2(Fig. 2)

图 2 收敛后G(x)=0与=0对比Fig.2 Comparison between G(x)=0 and =0 after convergence

3.2 实例2本节选用文献[2-3]中的桁架结构作为分析对象, 如图 3所示.该结构有10个输入随机变量, 分别为A₁, A₂, E₁, E₂, P₁, …, P₆.其中A₁, A₂分别为水平杆及倾斜杆的横截面积；E₁, E₂为水平杆及倾斜杆的杨氏模量；P₁, …, P₆为随机载荷.涉及到的10个随机变量的分布类型及参数如表 1所示.
图 3(Fig. 3)

图 3 桁架结构示意图Fig.3 Figure of the truss structure

表 1(Table 1)

表 1 实例2随机变量分布类型及参数Table 1 Distribution patterns and parameters for 2nd example 变量 分布类型 均值 标准差 P1~P6/N Gumbel 5×104 7.5×103 A1/m2 对数正态 2×10－3 2×10－4 A2/m2 对数正态 1×10－3 1×10－4 E1/Pa 对数正态 2.1×1011 2.1×1010 E2/Pa 对数正态 2.1×1011 2.1×1010

表 1 实例2随机变量分布类型及参数 Table 1 Distribution patterns and parameters for 2nd example

该结构在P₁, …, P₆随机载荷作用下点E在竖直方向位移最大, 因此以E点位移(s(x))作为结构响应, 为与文献[2]保持一致, 位移阀值设为0.14 m, 因此该结构功能函数为
此桁架结构的失效概率约为3.45×10^-5, 该结果来自500 000次重要抽样^[2].
令K=3 000, ε=0.001, k=8, e=0.15, 首先将10个随机变量变换为标准正态随机变量, 再应用第2节所提出方法计算该结构的失效概率.拉丁超立方抽样在[-5, 5]¹⁰立方体内进行(M₀=15).图 4为应用本文所提方法运行5次所得的迭代次数与失效概率预测值折线图, 由于失效概率预测值数量级变化比较大, 纵坐标采用对数形式.表 2为本文方法与其他方法精度对比.
图 4(Fig. 4)

图 4 Pf收敛折线图Fig.4 Broken lines of Pf during convergence

表 (Table )

表 实例2结果对比Table Comparison of results for 2nd example 来源 ×105 相对误差/% Mcall 参考值[2] 3.45 — — 文献[2]方法 2.35 -31.9 207 文献[3]方法 3.30 -4.3 142 本文方法 3.36~3.64 -2.6~5.5 81~89

表 实例2结果对比 Table Comparison of results for 2nd example

3.3 结果分析实例1主要用来具体说明本文所提出的结构可靠性分析方法.图 1表明, 本文所提出方法能够保证选取的样本点落在“重要区域”, 避免了对非重要区域的不必要抽样.同时, 各样本点间能够保持一定距离, 这也在很大程度上提高了算法的效率.实例2选用的桁架结构在其他文献[2-3]中已有提到, 通过对比可知, 本文所提出方法能够在较少调用结构功能函数情况下得到较高精度失效概率估计值.通过图 4可知, 本文方法收敛速度快且稳定.
4 结论1)?提出了一种自适应的选点方法, 通过实例表明该方法能够保证所选样本点落在对失效概率贡献较大的重要区域, 有效减少了对非重要区域的抽样次数, 提高各样本点的质量.
2)?提出的自适应方法能够保证各样本点间距离不会太小, 结合Kriging模型, 能够更充分利用各点统计信息, 提高实效概率计算精度.
3)?提出的机械结构可靠性分析方法能够在较少调用结构功能函数情况下实现对单设计点机械结构失效概率的较准确估计.
参考文献

[1]	Alvarez D A, Hurtado J E. An efficient method for the estimation of structural reliability intervals with random sets, dependence modeling and uncertain inputs[J].Computers and Structures, 2014, 142: 54–63.DOI:10.1016/j.compstruc.2014.07.006
[2]	Blatman G, Sudret B. An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[J].Probabilistic Engineering Mechanics, 2010, 25: 183–197.DOI:10.1016/j.probengmech.2009.10.003
[3]	Roussouly N, Petitjean F, Salaun M. A new adaptive response surface method for reliability analysis[J].Probabilistic Engineering Mechanics, 2013, 32: 103–115.DOI:10.1016/j.probengmech.2012.10.001
[4]	Echard B, Gayton N, Lemaire M. AK-MCS:an active learning reliability method combining Kriging and monte carlo simulation[J].Structural Safety, 2011, 33: 145–154.DOI:10.1016/j.strusafe.2011.01.002
[5]	佟操, 孙志礼, 杨丽, 等. 一种基于Kriging和Monte Carlo的主动学习可靠度算法[J].航空学报, 2015, 36(9): 2992–3001. ( Tong Cao, Sun Zhi-li, Yang Li, et al. An active learning reliability method based on Kriging and Monte Carlo[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 2992–3001.)
[6]	杨多和, 安伟光, 李铁钧. 基于人工神经网络的结构可靠性分析[J].兵工学报, 2007, 28(4): 495–498. ( Yang Duo-he, An Wei-guang, Li Tie-jun. Structural reliability analysis based on artificial neural network[J].Acta Armamentarii, 2007, 28(4): 495–498.)
[7]	Alibrandi U, Alani A M, Ricciardi G. A new sampling strategy for SVM-based response surface for structural reliability analysis[J].Probabilistic Engineering Mechanics, 2015, 41: 1–12.DOI:10.1016/j.probengmech.2015.04.001
[8]	张琦, 李兴斯. 基于Kriging模型的结构可靠性分析[J].计算力学学报, 2006, 23(2): 175–179. ( Zhang Qi, Li Xing-si. Analysis of structural reliability based on Kriging model[J].Chinese Journal of Computational Mechanics, 2006, 23(2): 175–179.)
[9]	Xiong F F, Liu Y, Xiong Y, et al. A double weighted stochastic response surface method for reliability analysis[J].Journal of Mechanical Science and Technology, 2012, 26(8): 2573–2580.DOI:10.1007/s12206-012-0425-4
[10]	Nguyen X S, Sellier A, Duprat F, et al. Adaptive response surface method based on a double weighted regression technique[J].Probabilistic Engineering Mechanics, 2009, 24: 135–143.DOI:10.1016/j.probengmech.2008.04.001
[11]	Gao Y, Turng L S, Wang X. Adaptive geometry and process optimization for injection molding using the Kriging surrogate model trained by numerical simulation[J].Advances in Polymer Technology, 2008, 27(1): 1–16.DOI:10.1002/adv.v27:1
[12]	Shi X, Teixeira A P, Zhang J, et al. Structural reliability analysis based on probabilistic response modelling using the maximum entropy method[J].Engineering Structures, 2014, 70(9): 106–116.