引　　言

工程结构往往存在大量不确定因素，可靠性方法是处理上述不确定性因素的有效手段之一。传统的概率可靠性方法已得到广泛应用，然而由于条件限制，许多情况下难以获得用于确定概率分布的足够样本信息。近年来，能够处理有限样本信息的非概率可靠性方法受到广泛关注。

Ben-Haim^[1]首次提出了非概率可靠性的概念。Elishakoff等^[2]针对此概念提出了一种基于非概率安全系数的可靠性度量，该度量不是一特定数值而是一区间值。与概率可靠性指标类似，郭书祥等^[3]提出了一种基于区间模型的非概率可靠性指标，并阐述了该指标的物理和几何意义。Jiang等^[4]发展了针对该可靠性指标的半解析法。曹鸿钧和段宝岩^[5]提出了一种基于椭球模型的非概率可靠性指标，该指标采用标准空间中原点到失效面的最短欧式距离度量可靠性。Jiang等^[6]针对上述可靠性指标提出中心点法和设计点法两种求解方法。Meng等^[7]提出了一种基于p范数的超参数凸模型可靠性指标。Meng等^[8]进一步提出了一种基于指数凸模型的可靠性指标。借鉴概率可靠度的概念，Wang等^[9]、Jiang等^[10]进一步提出了区间模型和椭球模型的非概率可靠度概念，得到了不确定变量域落入可靠域的体积与整个不确定变量域的体积的比值度量可靠性。Jiang等^[10]根据此概念进一步提出了椭球模型非概率可靠度的一阶近似法和二阶近似法。

上述可靠性分析方法均是针对显式结构功能函数，在复杂工程结构中通常面临的是隐式结构功能函数，对其进行可靠性分析通常需要采用计算成本较高、样本需求量大的有限元分析或Monte-Carlo模拟。鉴于此，采用计算量小、样本需求少的代理模型来代替原分析模型对结构进行可靠性分析成为了一种有效手段。在基于代理模型非概率可靠性分析方面，陈江义等^[11]利用响应面法对椭球模型与区间模型进行非概率可靠性分析。Bai等^[12]提出了一种基于不含交叉项二次响应面法的非概率可靠性分析方法。马超等^[13]用支持向量机法拟合结构功能函数，并将该方法应用于实际的飞机机翼可靠性分析中。潘林锋等^[14]用Kriging模型拟合回归极限状态方程，并用HL-RF修正法计算了非概率可靠性指标。郑严^[15]利用支持向量机和三种改进粒子群算法进行了结构可靠性分析及优化设计，主要解决了隐式结构概率及非概率可靠性分析问题。但Kriging模型的搜索方式对初始值具有依赖性，初始值的设置误差会导致后续模型建立时预测结果陷入局部最优的情况，且Kriging模型采用模式搜索法进行单点搜索，搜索路径单一。

针对Kriging模型的上述不足，本文提出了一种粒子群优化(particle swarm optimization，PSO)与Kriging模型相结合的非概率可靠性分析方法。本文的结构安排如下：第1节对基于椭球模型的非概率可靠性指标进行概述；第2节给出了本文所提方法及其实施流程；第3节通过三个算例验证了本文方法的有效性；第4节对本文进行总结。

1.   基于椭球模型的非概率可靠性指标

考虑某工程结构功能函数：

式中，${\boldsymbol{X}} = {\left( {{X_1},{X_2}, \cdots ,{X_n}} \right)^{\text{T}}}$表示结构的$n$维不确定向量，失效面$g\left( {\boldsymbol{X}} \right) = 0$为结构状态，分成安全$g\left( {\boldsymbol{X}} \right) \geqslant 0$与失效($g\left( {\boldsymbol{X}} \right) \lt 0$)两种。

本文将不确定参数X的不确定域用一个多维椭球来描述：

式中，E为$n$维椭球，X ⁰为椭球的中心，Ω为椭球的特征矩阵，其表达式为

式中，$x_n^{\text{r}}$表示变量的区间半径，${\rho _{ij}}$表示变量${x_i}$与变量${x_j}$的相关系数。更多有关多维椭球模型的信息，可以参考文献[16]。

为方便对结构可靠性进行定义，需首先将不确定参数${\boldsymbol{ X}}$线性变换为其标准化形式${\boldsymbol{\delta}}$，相应的结构功能函数和不确定域分别变换为$G\left( {\boldsymbol{\delta }} \right)$和一个中心在原点的单位超球${{\boldsymbol{\delta}} ^{\text{T}}}{\boldsymbol{\delta }} \leqslant 1$。图1给出了二维空间的线性变换。通过上述处理，可有效避免由于工程实际问题中不确定参数的数量级不同导致病态矩阵的产生，同时变量的无量纲化也便于结构可靠性的定义与比较。

图  1  线性变换与可靠性指标

Figure  1.  The linear transform and the reliability index

下载: 全尺寸图片 幻灯片

此时，结构的非概率可靠性指标$\eta$可定义为标准化空间内原点到失效面$G({\boldsymbol{\delta}} ) = 0$的最短距离，可通过如下优化问题求解：

上述可靠性指标与传统的概率可靠性指标几何意义相同，可采用HL-RF法进行求解。

2.   PSO-Kriging方法的非概率可靠性分析

考虑到许多工程实际问题中结构功能函数均为隐式的，以下发展一种将上述非概率可靠性指标与PSO-Kriging方法相结合的非概率可靠性分析方法。该方法主要由以下三个部分组成：① 构建Kriging模型；② 运用PSO方法对模型相关参数进行寻优；③ 基于代理模型求解可靠性指标。

对于代理模型的建立，首先需要选取合理分布的样本点。样本点的选取对后续建模的精度有一定影响，因此对于样本采样的最佳策略应在符合采样精度的前提下，采样时间最短、采样点应均匀分布在整个待测区域。Sobol序列抽样方法^[17]无论样本点$N$为多少，均可在面域内取得均匀分布的点，且不会出现局部团簇现象，故本文采用Sobol序列抽样方法进行抽样。

2.1   Kriging模型的构建

Kriging模型将需代理求解的未知函数看作是一个随机过程，其表达的隐函数为

式中，${\boldsymbol{f}}( {\boldsymbol{x}} ) = [ {{f_1}({\boldsymbol{ x}} ),{f_2}( {\boldsymbol{x}} ), \cdots } ,{{f_p}( {\boldsymbol{x}} )} ]^{\rm T}$为回归基函数，$\tilde {\boldsymbol{\beta}} = {[ {{\beta _1},{\beta _2}, \cdots ,{\beta _p}} ]^{\rm T}}$为基函数回归系数，二者乘积为确定部分，${\boldsymbol{z}}( {\boldsymbol{x}} )$为一个均值为0、方差为${\sigma ^2}$、协方差为${\sigma ^2}R( {\theta ,{{\boldsymbol{x}}_i},{{\boldsymbol{x}}_j}} )$的随机过程，${{\boldsymbol{x}}_i}$和${{\boldsymbol{x}}_j}$为两任意输入向量，R为不同样本点组成的相关模型，表示各样本点之间的相关性，可表示为如下的Gauss函数：

式中，${\theta _k}$表示函数输出对第$k$维输入变量的敏感程度，${\theta _k}$增加时，R随着样本距离$| {{\boldsymbol{x}}_i^{( k )} - {\boldsymbol{x}}_j^{( k )}} |$的增加而急剧减小。可见，相关参数影响着Kriging模型的输入和输出特性，对其进行优化从而保证预测精度是必要的。

假设$m$个样本点的响应值为${\boldsymbol{Y}} = {[ {{\boldsymbol{y}}( {{{\boldsymbol{s}}_1}} ),{\boldsymbol{y}}( {{{\boldsymbol{s}}_2}} ), \cdots ,{\boldsymbol{y}}( {{{\boldsymbol{s}}_m}} )} ]^{\rm T}}$。根据模型无偏估计的需要，使得预测值${\boldsymbol{y}}( {\boldsymbol{x}} )$的方差最小化，得到$\beta$和${\sigma ^2}$的无偏估计值为

其中，F为回归函数的基函数矩阵，R为相关矩阵。在极大似然估计的条件下相关参数的最优解为

上述最优相关参数的确定为构建最优Kriging模型奠定了基础，进而可求解相应的可靠性指标。

2.2   PSO方法求解模型相关参数

PSO方法采用多点并行搜索，利用群体中的各个粒子对信息进行共享，每个粒子根据个体极值和群体极值的变化方向判断自身的搜索方向，从而使得整个群体的运动在求解空间中产生从无序到有序的演变过程，无论初始值的位置如何，都可大概率寻找到最优解。本文采用PSO对相关参数$\theta$进行寻优。

在寻优过程中，以式(9)作为适应度函数，以Kriging模型的$n$个样本点作为粒子，其中第$i$个粒子位置为${x_i}$，根据适应度函数，将第$i$个粒子的个体极值${p_{\rm{best}}}$记为${p_i}$，群体最优位置极值${g_i}$为所有${p_{\rm{best}}}$的最优解。每个粒子的移动速度记为${v_i}$，则粒子根据每次位置判断下次更新的速度和位置为

式中，$k$为当前迭代数，${\rm{rand}}$为$\left( {0,1} \right)$之间的随机数，$\omega$为惯性权重因子，更新的学习因子分别为${c_1}$，${c_2}$。当粒子群搜索满足式(9)时，最优$\theta$值即为该粒子的全局最优位置${g_{\rm{best}}}$。此时构建的Kriging模型为最优模型，可据此模型拟合出功能函数，利用式(4)计算设计点和可靠性指标$\eta$。

由于代理模型是近似拟合，因此需要多次拟合结构功能函数并求解可靠性指标，直至满足收敛条件$\left| {{\eta _i} - {\eta _{i - 1}}} \right| \leqslant \varepsilon$，$\varepsilon$为给定的允许值。

2.3   所提方法实施过程及流程

本文所提方法流程如图2所示。

图  2  算法流程图

Figure  2.  The flowchart for the proposed method

下载: 全尺寸图片 幻灯片

步骤1　利用Sobol序列抽样选取$n$个样本点，根据隐式函数确定对应的函数值。

步骤2　PSO搜索最优$\theta$值的过程：

1) 设置初始参数。设置PSO初始种群个数为$N$，设定最大迭代次数为$I_{\rm{ger }}$，并将惯性权重设为$\omega$，学习因子为${c_1}$，${c_2}$，初始种群的位置以及粒子的位置参数及速度限制范围。

2) 通过式(9)计算每个粒子的个体适应度值，比较当前粒子的适应度与历史最佳适应度${p_{\rm{best}}}$，如果当前位置的适应值更高，更新粒子群的粒子的个体最优适应度${p_{\rm{best}}}$。

3) 将粒子当前位置的适应值与其全局最佳位置${g_{\rm{best}}}$对应的适应值进行比对，如果当前位置的适应值更高，则用当前位置更新全局最佳位置${g_{\rm{best}}}$。

4) 更新每个粒子的速度与位置。

5) 判断是否满足式(9)，若满足，则得出最优的$\theta$；若不满足，则继续迭代。

步骤3　得出最优的$\theta$值，构建相应的Kriging模型。

步骤4　由预测出的Kriging模型求解非概率可靠性指标及设计点坐标。

步骤5　若可靠性指标满足收敛条件$\left| {{\eta _i} - {\eta _{i - 1}}} \right| \leqslant \varepsilon$，则停止；若不满足，则在当前设计点附近重新确定样本范围，返回步骤1。

3.   算　　例

为验证本文方法的精度和效率，对三个算例进行分析，其中包括一个数值算例和两个工程算例，计算结果与Kriging方法和Monte-Carlo (MCS)方法的结果进行对比。算例中任意两个变量相关系数均设定为$\rho$，相对误差的计算公式为$\left| {\dfrac{{\eta - {\eta _{{\text{exact}}}}}}{{{\eta _{{\text{exact}}}}}}} \right| \times 100{\text{% }}$，$\eta$为采用本文方法或者Kriging方法给出的分析结果；${\eta _{{\text{exact}}}}$表示HL-RF法或者MCS方法(样本为10⁸)给出的参考解，为增加MCS给出结果的可信度，针对不同相关参数，进行100次MCS计算，根据所得结果求解MCS法的均值、离差以及95%置信区间的上下界。本文所有程序均在Intel core i3-10100 3.60 GHz 四核处理器上运行。

算例1　数值算例

式中，${X_1} \in \left[ {0.0,1.0} \right]$，${X_2} \in \left[ {0.0,1.0} \right]$。在本算例中，样本$M = 50$，参数${c_1}$，${c_2}$均设置为0.3，惯性权重$\omega$均设为0.75，最大迭代次数为10次，速度限制为$V \in \left[ { - 1,1} \right]$；迭代收敛标准为$\left| {{\eta _i} - {\eta _{i - 1}}} \right| \leqslant 0.001$。

表1给出了相关系数为0.2时本文方法的迭代过程；表2列出了不同相关系数下，本文方法、Kriging方法和HL-RF法(参考值)的对比结果；图3给出了不同相关系数下两种方法的相对误差。

表  1  相关系数为0.2的迭代过程

Table  1.  The iterative process for a correlation coefficient of 0.2

iteration step	design point	reliability index
0	(−2.4442, 4.2332)	4.8882
1	(−1.8607, 4.1064)	4.5083
2	(−1.8641, 4.1039)	4.5074
reference value	(−1.8692, 4.1017)	4.5075

下载: 导出CSV 

| 显示表格

表  2  不同相关系数的可靠性指标（算例1）

Table  2.  Reliability indexes for different correlation coefficients (example 1)

correlation coefficient	PSO-Kriging	iterations	CPU time	Kriging	iterations	CPU time	reference value
0	4.1704	2	0.45703 s	4.1653	5	1.00112 s	4.1701
0.2	4.5073	2	0.38524 s	4.5002	6	1.20133 s	4.5075
0.5	5.2134	3	0.68554 s	5.2039	5	1.01367 s	5.2138
0.7	5.9237	3	0.88563 s	5.9355	7	1.40156 s	5.9235
0.9	7.0337	2	0.57889 s	7.0302	6	1.16778 s	7.0339

下载: 导出CSV 

| 显示表格

算例2　考虑图4所示的刮板输送机GB/T 12718链轮，在其链窝处受到一个120 kN的切向力。在本算例中，齿弧半径${R_1}$，齿根圆弧半径${R_2}$，链窝部的平面圆弧半径${R_3}$，短齿根部的圆弧半径${R_4}$以及平环中面至链轮中心距离H均为未知但有界的不确定变量，其区间分别为：$R_1^{\text{I}} = \left[ {30,40} \right]{\text{ mm}}$，$R_2^{\text{I}} = \left[ {8,10} \right]{\text{ mm}}$，$R_3^{\text{I}} = \left[ {24,36} \right]{\text{ mm}}$，$R_4^{\text{I}} = \left[ {7,11} \right]{\text{ mm}}$和${H^{\text{I}}} = \left[ {125,151} \right]{\text{ mm}}$，齿顶位移最大许用值为${d_{{\text{max}}}}$。

图  3  相对误差对比图（算例1）

Figure  3.  Comparison of relative errors (example 1)

下载: 全尺寸图片 幻灯片

图  4  链轮

Figure  4.  A chain wheel

下载: 全尺寸图片 幻灯片

链轮的功能函数为

本算例中样本$M = 30$，参数${c_1}$，${c_2}$均设置为0.2，惯性权重$\omega$均设为0.75，最大迭代次数为10次，速度限制为$V \in \left[ { - 1,1} \right]$，迭代收敛标准为$\left| {{\eta _i} - {\eta _{i - 1}}} \right| \leqslant 0.001$。

表3给出了相关系数为0.9时本文方法的迭代过程；表4列出了不同相关系数下，本文方法、Kriging方法和参考值(MCS法)的对比结果；表5给出了MCS求解的95%置信区间，表明参考值是有效的；图5给出了不同相关系数下两种方法的相对误差。

表  3  相关系数为0.9的迭代过程

Table  3.  The iterative process for a correlation coefficient of 0.9

iteration step	design point	reliability index
0	(0.0662, 0.0265, 0.0088, 0.00650.4059)	0.4123
1	(−0.0031, −0.0004, 0.0146, 0.0141, 0.3458)	0.3464
2	(−0.0040, −0.0005, 0.0124, 0.0176, 0.3415)	0.3422
3	(−0.0044, −0.0007, 0.00132, 0.0179, 0.3441)	0.3448
4	(−0.0042, −0.0007, 0.0134, 0.0181, 0.3439)	0.3447
reference value	(−0.0058, −0.0015, 0.0146, 0.0155, 0.3437)	0.3444

下载: 导出CSV 

| 显示表格

算例3　如图6所示的25杆桁架结构^[10]，桁架弹性模量$E = 199\;949.2{\text{ MPa}}$，Poisson比$\,\mu = 0.3$，横向、纵向杆件的长度为$L$，杆件① ~ ④、⑯ ~ 、⑪ ~ ⑮、⑤ ~ ⑩的截面积分别表示为${A_1}$，${A_2}$，${A_3}$和${A_4}$。作用在节点7、9和11的竖直载荷分别为${F_3} = 1\;779.2{\text{ kN}}$，${F_2} = 2\;224{\text{ kN}}$和${F_1} = 1\;779.2{\text{ kN}}$，节点1的水平载荷${F_4} = 1\;334.4{\text{ kN}}$，节点6水平位移的最大许用值为${d_{{\text{max}}}}$。由于制造与测量过程中出现的误差，杆件的截面积${A_i}\left( {i = 1,2,3,4} \right)$以及杆件的长度L为不确定变量，其不确定区间如表6所示。

表  4  不同相关系数的可靠性指标（算例2）

Table  4.  Reliability indexes for different correlation coefficients (example 2)

correlation coefficient	PSO-Kriging	iterations	CPU time	Kriging	iterations	CPU time	reference value
0	0.3469	4	16.15572 s	0.3406	7	27.16953 s	0.3465
0.3	0.3504	5	20.19465 s	0.3466	9	34.93226 s	0.3509
0.5	0.3517	4	17.23898 s	0.3612	7	26.35562 s	0.3513
0.7	0.3487	5	21.56773 s	0.3454	8	31.05089 s	0.3490
0.9	0.3447	4	18.06577 s	0.3404	7	28.09947 s	0.3444

下载: 导出CSV 

| 显示表格

表  5  MCS求解的95%置信区间（算例2）

Table  5.  The 95% confidence intervals of MCS solutions (example 2)

MCS value	mean value	deviation	lower bound	upper bound
0.3465	0.3469	0.0029	0.3463	0.3475
0.3509	0.3506	0.0029	0.3500	0.3511
0.3513	0.3511	0.0030	0.3505	0.3517
0.3490	0.3491	0.0030	0.3485	0.3497
0.3444	0.3448	0.0031	0.3442	0.3454

下载: 导出CSV 

| 显示表格

图  5  相对误差对比图（算例2）

Figure  5.  Comparison of relative errors (example 2)

下载: 全尺寸图片 幻灯片

图  6  25杆桁架

Figure  6.  A 25-bar truss structure

下载: 全尺寸图片 幻灯片

表  6  不确定变量的分布参数

Table  6.  The distribution parameters of uncertain variables

uncertain variable	mean value	radius
${A_1}$/mm2	700	70
${A_2}$/mm2	6200	620
${A_3}$/mm2	5300	530
${A_4}$/mm2	8800	880
$L$/mm	15000	1500

下载: 导出CSV 

| 显示表格

该桁架极限状态函数为

本算例中样本$M = 30$，参数${c_1}$，${c_2}$均设置为0.2，惯性权重$\omega$均设为0.75，最大迭代次数为10次，速度限制为$V \in \left[ { - 1,1} \right]$，迭代收敛标准为$\left| {{\eta _i} - {\eta _{i - 1}}} \right| \leqslant 0.01$。

表7给出了相关系数为0时本文方法的迭代过程；表8列出了不同相关系数下，本文方法、Kriging方法和参考值(MCS法)的对比结果；表9给出了MCS求解的95%置信区间，表明参考值是有效的；图7给出了不同相关系数下两种方法的相对误差。

表  7  相关系数为0的迭代过程

Table  7.  The iterative process for a correlation coefficient of 0

iteration step	design point	reliability index
0	(−0.0083, −0.4266, −0.3646, −0.0519, 0.8256)	0.9997
1	(−0.0109, −0.5652, −0.4831, −0.0587, 0.9338)	1.1952
2	(−0.0100, −0.5201, −0.4446, −0.0540, 0.8593)	1.0997
3	(−0.0089, −0.4827, −0.4126, −0.0478, 0.7875)	1.0870
4	(−0.0084, −0.5107, −0.4222, −0.0509, 0.8504)	1.0793
5	(−0.0105, −0.5065, −0.4289, −0.0465, 0.8372)	1.0694
6	(−0.0094, 0.5032, −0.4301, −0.0431, 0.8356)	1.0669
reference value	(−0.0092, −0.5043, −0.4280, −0.0484, 0.8340)	1.0656

下载: 导出CSV 

| 显示表格

表  8  不同相关系数的可靠性指标（算例3）

Table  8.  Reliability indexes for different correlation coefficients (example 3)

correlation coefficient	PSO-Kriging	iterations	CPU time	Kriging	iterations	CPU time	reference value
0	1.0669	6	27.02193 s	0.9944	9	37.73403 s	1.0656
0.3	1.2683	7	31.52558 s	1.1793	10	41.92578 s	1.2679
0.5	1.4906	6	28.02193 s	1.3997	10	40.44356 s	1.4917
0.7	1.9043	6	26.02193 s	1.8253	9	38.37882 s	1.9023
0.9	3.1439	7	31.52558 s	2.9934	11	46.11937 s	3.1431

下载: 导出CSV 

| 显示表格

表  9  MCS求解的95%置信区间（算例3）

Table  9.  The 95% confidence intervals of MCS solutions (example 3)

MCS value	mean value	deviation	lower bound	upper bound
1.0656	1.0659	0.0029	1.0653	1.0665
1.2679	1.2672	0.0028	1.2667	1.2679
1.4917	1.4919	0.0025	1.4914	1.4924
1.9023	1.9023	0.0028	1.9017	1.9029
3.1431	3.1430	0.0029	3.1424	3.1435

下载: 导出CSV 

| 显示表格

由表1、3、7可以看出，本文所提的PSO-Kriging方法均能逐渐收敛于参考值附近，且与参考值误差较小。

由表2、4、8可知，不同相关系数下，相比Kriging方法，PSO-Kriging方法的计算时间较短，且给出的可靠性指标更接近参考解，表明本文方法具备更高的计算效率与精度。

由图3、5、7可以看出，PSO-Kriging方法的相对误差对相关系数的变化不敏感，具有较强的鲁棒性，而Kriging方法则对相对系数的变化较敏感。

图  7  相对误差对比图（算例3）

Figure  7.  Comparison of relative errors (example 3)

下载: 全尺寸图片 幻灯片

4.   结　　论

本文提出了一种PSO与Kriging模型相结合的非概率可靠性分析方法，该方法采用PSO方法得到最优相关参数，避免了Kriging模型的单向搜索及对初始值敏感的缺点。三个算例分析结果表明，本文方法有效可行，且精度和效率均高于Kriging模型非概率可靠性分析方法。