随着现代结构产品功能和复杂性的增加，结构的可靠性越来越受到人们的重视. 实际工程结构的极限状态函数(limit state function，LSF)通常无法显式表达，为了精确评估失效概率，所需计算量太大，从而限制了大多数可靠性分析方法在实际工程结构中的应用. 代理模型也称为元模型或近似模型，是提高实际工程结构优化与分析效率的重要方法^[1]. 常用的代理模型方法有响应面法、径向基函数方法、Kriging模型、神经网络法和支持向量机方法等. 其中支持向量机因其具有较强的小样本学习能力和泛化能力，已广泛应用于工程结构可靠性分析中^[2].

结构可靠性分析主要通过判断结构系统的极限状态函数是否超过许用值来评估结构的失效概率. 小于许用值则认为结构处于安全域，反之则认为结构处于失效域^[3]. 代理模型方法通常利用少量样本点来建立精度较高的显式或者近似极限状态函数. 不同的试验设计方法对代理模型精度有显著影响. Jin等^[4]指出对于确定性的试验设计方法，样本点应分布在整个参数空间内，同时表明在相同的样本量与代理模型方法下，样本点分布越均匀，模型精度越高. 拉丁超立方抽样方法(Latin hypercube sampling，LHS)是一种广泛应用的试验设计方法. 虽然能够同时考虑可行域与不可行域，但难以保证训练样本的空间均匀性. 近年发展起来的质心Voronoi网格试验设计方法^[5-8](centroidal Voronoi tessellation，CVT)和拉丁质心Voronoi网格试验设计方法^[9-10](Latinized centroidal Voronoi tessellation，LCVT)较其他试验设计方法有更好的均匀性.

本文结合支持向量机回归方法与拉丁质心Voronoi网格试验设计方法，并将该方法应用于结构可靠性分析领域. 与其他抽样方法进行对比，在相同数量的训练样本集下，获取更多的设计变量空间信息，以建立更高精度的结构极限状态函数代理模型. 利用1个数值算例及2个工程算例，验证所提方法的计算性能和可行性.

1 质心Voronoi网格

给定一个开区域为 $\varOmega \subseteq {{\bf R}^N}$ ,定义一组初始点集为 ${{Z}} = \left\{ {{z_i}} \right\}_{i = 1}^n$ ， ${{Z}} \in \varOmega $ 根据下式对区域 $\varOmega $ 进行划分：

$\left.\begin{aligned} {v_i} = \left\{ {y \in \varOmega \left| {{{\left\| {y - {z_i}} \right\|}_2} < {{\left\| {y - {z_j}} \right\|}_2}} \right.} \right\}{\text{ }}{\text{;}} \\ j = 1,2,\,\cdots\,,n,\;j \ne i. \end{aligned}\right\}$

(1)

式中： $\left\| \bullet \right\|$ 为 ${{\bf R}^N}$ 中的欧式范数， ${z_i}$ 为Voronoi网格的发生元. 每一个点 ${z_i}$ 对应于一个Voronoi区域 ${{v}}$ ，且 ${v_i} \cap {v_j} = \phi $ ， $ \cup _{i = 1}^n{v_i} = \varOmega $ ，

$z_i^ * = \frac{{\int\limits_v {y\rho \left( y \right){\rm d}y} }}{{\rho \left( y \right){\rm d}y}}{\text{ }}.$

(2)

当 ${z_i} = z_i^ * $ 时， $i = 1,2,\,\cdots\,,n_i^ * $ ，即当Voronoi网格的发生元等于质心时，称为质心Voronoi网格(centroidal Voronoi tessellation，CVT). 能量函数(3)具有最小值：

$F\left( {{z_i},{v_i}} \right) = \sum\limits_{i = 1}^n {\int\limits_{y \in {v_i}} {\rho \left( y \right)} } {\left\| {y - {z_i}} \right\|^2}{\rm d}y{\text{ }}{\text{.}}$

(3)

如图1所示为样本点数为20的随机 Voronoi 图和 CVT 图.

2 拉丁质心Voronoi网格试验设计

CVT抽样方法是在质心Voronoi网格(CVT)的基础上发展出来的一种试验设计方法^[5-6]. Voronoi网格是对区域按空间距离进行的一个划分，CVT抽样方法是寻找一种最优的Voronoi划分，使得CVT抽样样本点为质心Voronoi网格中的质心点. CVT抽样方法具有良好的稳定性，在确定采样点之后，对应的质心Voronoi网格也确定了. Saka等^[9]结合LHS方法和CVT抽样方法，提出新的抽样方法，称为拉丁质心Voronoi网格抽样方法. 对质心Voronoi网格进行特殊处理，使其具有拉丁超立方体的性质.

若LHS的样本量为 $n$ ，维数为 $d$ ，则可以定义 $H = {[0,1]^d}$ 表示 $d$ 维[0,1]区间的拉丁超立方体. ${z_i}\left( j \right)$ 为在第 $j$ 维空间内第 $i$ 点的坐标. 定义 ${P_j}$ 表示由数字1~ $n$ 构成的集合 $\left\{ {1,2, \cdots ,n} \right\}$ 中的一种排列组合方式，通常为顺序排列， ${P_j}\left( i \right)$ 为在第 $j$ 种排列组合中第 $i$ 个点的坐标，

$\left.\begin{aligned} {z_i}\left( j \right) = \frac{{{P_j}\left( i \right) - {U_{ji}}}}{n}{\text{ }}; \hfill \\ {\text{ }}i = 1,2,\,\cdots\,,n{\text{ }},\;j = 1,2,\,\cdots\,,d{\text{ }}{\text{.}} \hfill \\ \end{aligned} \right\}$

(4)

式中： ${U_{ij}}$ 为在第 $j$ 维空间内均匀分布的随机数，通常采用中值拉丁超立方抽样，即 ${U_{ij}} = 0.5$ . 可得一个拉丁超立方体的中心坐标为

$\begin{gathered} {z_i}\left( j \right) = \frac{{{P_j}\left( i \right) - 0.5}}{n}; \; i = 1,2,\,\cdots\,,d,\;j = 1,2,\,\cdots\,,n{\text{ }}{\text{.}} \hfill \\ \end{gathered} $

(5)

经过转换后，将非拉丁超立方体的抽样点调整为拉丁超立方体的抽样点，该调整称为拉丁化. 具体步骤如下：首先定义一个抽样点集 ${{Z}} = \left\{ {{z_i}} \right\}_{i = 1}^n$ ，正整数 $k$ ，满足 $1 \leqslant k \leqslant d$ .

1)第 $k$ 次对抽样点集 $Z$ 中的点依次进行转换，记为 $\left\{ {{S_k}{z_i}} \right\}_{i = 1}^n$ ，转换方式为

$\begin{aligned} {S_k}{z_i}\left( j \right) = \left\{ {\begin{aligned} &{{z_i}\left( j \right),{\text{ }}{\text{ }}j \ne k}; \\ &{\displaystyle\frac{{i - {U_{ji}}}}{n},{\text{ }}{\text{ }}j = k}. \end{aligned}} \right.\end{aligned}$

(6)

2)第 $k$ 次对抽样点集 $Z$ 中的点进行重新排列，记为 $\left\{ {{R_k}{z_i}} \right\}_{i = 1}^n$ ，根据第 $i$ 个 ${z_i}$ 在第 $k$ 维空间内的坐标进行排列.

3)对于任意一个抽样点集 $Z = \left\{ {{z_i}} \right\}_{i = 1}^n$ ，对应的拉丁化抽样点集 $\left\{ {{{\rm L{z}}_i}} \right\}_{i = 1}^n$ 可由下式计算得到：

${{\rm L{z}}_i} = \left( {\displaystyle\Pi _{k = 1}^d\left( {{S_k}{R_k}} \right)} \right){z_i}{\text{ }}{\text{.}}$

(7)

对质心Voronoi网格进行拉丁化，获得的新方法称为拉丁质心Voronoi网格. 如图2所示为CVT抽样点与LCVT抽样点.

3 基于LCVT-SVR的可靠性分析方法 3.1 支持向量机

支持向量机方法(support vector machine，SVM)是建立在结构风险最小化原则和统计学习理论的VC维理论基础上的一种机器学习方法，可以用于分类问题与回归估计. 当SVM应用于回归估计的时候称为支持向量机回归(support vector regression，SVR). 利用SVR进行回归估计.

给定训练样本集：

$\begin{gathered} S = \left\{ {{{\left( {{{{x}}_i},{{{y}}_i}} \right)}^l}} \right\};\; {{{x}}_i} \in {{\bf R}^n}\,,\;{{{y}}_i} \in {{\bf R}^m}\,,\; i = 1,2, \,\cdots\, ,m{\text{ }}{\text{.}} \hfill \\ \end{gathered} $

(8)

式中： $n$ 为输入参数维数， $m$ 为输出响应维数， ${{{x}}_i}$ 为输入向量， ${{{y}}_i}$ 为输出响应. 当训练样本为线性时，SVR的回归函数可由下式表示：

$f({{x}}) = ({{w}} \cdot {{x}}) + {w_0}{\text{ }}{\text{.}}$

(9)

根据结构风险最小化原则，对式(9)的求解转化为优化问题：

$\left.\begin{gathered} \left\| {{{ w}^2}} \right\| = {{\left( {{{w}} \cdot {{w}}} \right)} / 2}; \hfill \\ {\text{s}}{\text{.t}}{\text{. }}\left| {{y_i} - \left( {{{w}} \cdot {{{x}}_i}} \right) + {w_0}} \right| \leqslant \varepsilon ,i = 1,2,\cdots\,,m{\text{ }}{\text{.}} \hfill \\ \end{gathered} \right\}$

(10)

式中： ${{w}}$ 和 ${w_0}$ 为算法控制参数， $\varepsilon $ 为不敏感损失函数.

当约束条件不满足时，引入松弛变量 $\xi $ 和惩罚因子 $c$ ， $c$ 表示对误差超过 $\varepsilon $ 的惩罚程度. 采用Lagrange函数，将该优化问题转换为二次规划问题：

$\begin{split} \min &{\text{ }}\frac{1}{2}\sum\limits_{i,j}^m {\left( {{\alpha _i} - \alpha _i^*} \right)} \left( {{\alpha _i} - \alpha _j^*} \right)\left( {{{{x}}_i} \cdot {{{x}}_j}} \right) +\\ &\varepsilon \sum\limits_{i = 1}^m {\left( {{\alpha _i} + \alpha _i^*} \right)} - \sum\limits_{i = 1}^m {{y_i}\left( {{\alpha _i} - \alpha _i^*} \right)} {\text{ }}{\text{.}} \end{split} $

(11)

式中： ${\alpha _i}$ 和 $\alpha _i^*$ 为Lagrange乘子.

约束条件为

$\left.\begin{aligned} {\text{s}}{\text{.t}}{\text{. }}\sum\limits_{i = 1}^m {\left( {{\alpha _i} - \alpha _i^*} \right)} = 0;{\text{ }} \hfill \\ {\text{ 0}} \leqslant {\alpha _i},\;{\text{ }}\alpha _i^* \leqslant c{\text{ }}\left( {i = 1,2,\,\cdots\,,m} \right){\text{ }}{\text{.}} \hfill \\ \end{aligned} \right\}$

(12)

求得最优回归函数：

$f\left( x \right) = \sum\limits_{{x_i} \in { S}_{\rm V}} {\left( {{\alpha _i} - \alpha _i^*} \right)} \left( {{{{x}}_i} \cdot {{{x}}_j}} \right) + {w_0}{\text{ }}{\text{.}}$

(13)

式中： ${ S}_{\rm V}$ 为支持向量.

对于非线性训练样本，通过非线性变换 $x \to \varPhi \left( x \right)$ ，将输入空间映射到高维特征空间 $Z$ ，在特征空间内构造线性回归估计. 内积的计算可以用满足Mercer条件的核函数 $K\left( {{{{x}}_i},{{{x}}_j}} \right)$ 得到. 将线性SVM中的 $\left( {{{{x}}_i} \cdot {{{x}}_j}} \right)$ 用核函数 $K\left( {{{{x}}_i},{{{x}}_j}} \right)$ 替换，此时求得最优回归函数：

$f\left( {{x}} \right) = \sum\limits_{{x_i} \in SV} {\left( {{\alpha _i} - \alpha _i^*} \right)} K\left( {{{{x}}_i},{{{x}}_j}} \right) + {w_0}{\text{ }}{\text{.}}$

(14)

3.2 抽样策略

在结构可靠性分析中，各随机变量通常服从某一种概率分布. 逆变换法是目前应用最广泛的随机变量产生方法. 令随机变量 ${{ x}_i}$ 的累积分布函数为 ${F_X}\left( {{{{x}}_i}} \right)$ . 根据定义 ${F_X}\left( {{{{x}}_i}} \right) \in \left[ {0,1} \right]$ ，令 ${u_i}$ 为区间[0,1]上均匀分布的随机数，则由逆分布函数 ${F^{ - 1}}\left( {{u_i}} \right)$ 求得随机变量 ${{{x}}_i}$ 为

${{{x}}_i} = {F^{ - 1}}\left( {{u_i}} \right){\text{ }}{\text{.}}$

(15)

对于SVR代理模型，训练样本的均匀性和覆盖性是至关重要的. 对于小失效概率的结构极限状态函数，上述简单随机抽样方法可能经过多次抽样不出现失效点，不宜作为SVR的训练样本. 参考文献[11]，采用LCVT抽样方法在 $d$ 维[0,1]区间得到样本点 ${{Z}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ 后，不采用逆变化方法得到 ${{{x}}_i}$ ，而是将 ${{Z}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ 映射到区间 $[{{{\mu}} _i} - m{{{\sigma}} _i},{{{\mu}} _i} + m{{{\sigma}} _i}]$ 中，以获得更多的失效样本点. 构建更高精度的代理模型. 其中， ${{{\mu}} _i}$ 和 ${{{\sigma}} _i}$ 分别为随机变量 ${{{x}}_i}$ 的均值和标准差， $m$ 的大小要保证样本中包含一定数量的失效样本点.

对于 $m$ 的选取，采用“试错法”. 先取 $m = 3.5$ ，采样后根据正负响应值进行调整(正的响应值位于安全域). 若负的响应值占抽样总数的10%~40%，则认为本次抽样符合要求，否则将对 $m$ 进行调整，保证正、负响应值在上述范围内.

3.3 实现步骤

结构可靠性分析的LCVT-SVR方法的实现步骤如下.

1)确定结构输入变量数目 $n$ 及其概率分布，确定构建SVR模型所需要的样本数 $k$ ，每个样本都是n维向量.

2)随机生成样本 ${{Z}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ 作为Voronoi区域的发生元，根据式(2)计算各区域质心，通过式(3)不断调整直至生成满足质心Voronoi网格的样本 ${{{Z}}_{\rm {cvt}}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ .

3)根据式(7)不断调整 ${{{Z}}_{\rm {cvt}}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ ，得到LCVT抽样下的样本点 ${{{Z}}_{\rm {lcvt}}} = \left\{ {{z_i}} \right\}_{i = 1}^k$ .

4)根据各随机变量的概率分布，由抽样策略将LCVT抽样的样本点转换为SVR的训练样本 $X_{\rm{tr}}$ ，同时得到LSF的响应值 $Y_{\rm{tr}}$ ，SRS抽样生成测试样本 $X_{\rm{t}}$ ，得到LSF的响应值 $Y_{\rm{t}}$ .

5)对 $\left\{ {X_{\rm{tr}},Y_{\rm{tr}}} \right\}$ 和 $\left\{ {X_{\rm t},Y_{\rm t}} \right\}$ 进行归一化，得到 $\left\{ {X^{'}_{\rm{tr}},Y^{'}_{\rm{tr}}} \right\}$ 和 $\left\{ {X^{'}_{\rm{t}},Y^{'}_{\rm{t}}} \right\}$ ，以避免不同量纲的输入变量对SVM训练结果的影响. 归一化公式为

${{x}}_i^* = \frac{{2\left( {{{{x}}_i} - {{{x}}_{i{\min }}}} \right)}}{{{{{x}}_{i{\max }}} - {{{x}}_{i{\min }}}}} - {\bf 1}{\text{ }}{\text{.}}$

(16)

式中：x_i为归一化前的数据，x_i^*为归一化后的数据，x_imax和x_imin分别为归一化前数据的最大值和最小值.

6)通过 $\left\{ {X^{'}_{\rm{tr}},Y^{'}_{\rm{tr}}} \right\}$ 训练SVR，通过 $\left\{ {X^{'}_{\rm t},Y^{'}_{\rm t}} \right\}$ 测试满足精度后得到SVR代理模型，同时进行反归一化得到真实的输出值. 反归一化与归一化互逆.

7)根据得到的SVR代理模型，运用蒙特卡罗方法计算失效概率：

${P_{\rm f}} = P\left( {g\left( {{x}} \right) \leqslant 0} \right) \approx P\left( {f\left( {{x}} \right) \leqslant 0} \right) = \frac{{{N_{\rm f}}}}{N}{\text{ }}{\text{.}}$

(17)

式中： $g\left( { x} \right)$ 为真实的极限状态函数， $f\left( { x} \right)$ 为SVR代理模型近似的极限状态函数， $N$ 为按随机变量概率密度函数抽取的样本数， ${N_{\rm f}}$ 为落入失效域 $f\left( {{x}} \right) \leqslant 0$ 的样本数.

4 算例分析

分别采用1个数值算例和2个工程算例，研究所提方法在结构可靠性中的应用，考察计算精度. SVR代理模型基于LIBSVM^[12]工具箱建立.

4.1 算例1

考察4次极限状态函数的可靠性分析，通过蒙特卡罗计算，得到失效概率为0.001 84^[13].

$g\left( {{x}} \right) = 2 + \exp \left( { - \frac{{x_1^2}}{{10}}} \right) + {\left( {\frac{{x_1^2}}{5}} \right)^4} - {x_2}{\text{ }}{\text{.}}$

(18)

式中： ${x_1}$ 和 ${x_2}$ 为标准正态随机变量.

对于常用的抽样方法，如简单随机抽样(simple random sampling，SRS)、LHS抽样、CVT抽样、LCVT抽样、均匀抽样(uniform sampling，US)以及Sobol、Halton、Hammersley等序列抽样方法，分别生成20个样本点，并构建SVR代理模型. 随机采样10 000个点作为测试样本. 其中 $m$ 取6. 为了便于比较，SVR参数固定，分别取 $C = 1\,000,g = 0.01$ . 在得到SVR代理模型之后，采用相同的MCS样本计算失效概率，样本量为10⁶. 极限状态方程的近似情况如图1所示.

图3给出不同抽样方法下LSF的近似情况. 可以看出，SRS抽样与US抽样方法LSF的近似情况较差，其他抽样方法的近似情况较好. 表1给出不同抽样方法的计算结果. 表中， ${\text{MSE}}$ 为均方误差， ${{{R}}^2}$ 为回归相关系数， $\delta $ 为利用不同抽样方法计算得到的 ${P_{\rm f}}$ 与理论失效概率的相对误差. 由表1可知，SRS和US抽样的 ${\text{MSE}}$ 较大， ${{{R}}^2}$ 较小， $\delta $ 较大. SRS和US抽样的SVR方法精度较低，效果较差，与图1的结果对应. 表1中，LCVT抽样方法的各评价指标值均较好. 此外，LHS和Halton抽样方法表现较好. 为了避免一次计算的偶然性，表2给出不同抽样方法独立运行30次的统计结果. 表中，SD为标准差，COV为变异系数.

由表2可知，SRS和US抽样的SD和COV较大，表明SRS-SVR和US-SVR计算结果不稳定，计算精度不高，鲁棒性较差，这与图1和表1的结果相互验证. LCVT抽样均值更加接近理论值，同时SD较小，表明LCVT-SVR计算结果稳定，同时精度较高，鲁棒性较好，可以认为计算结果是可靠的.

4.2 算例2

如图4所示为各向异性的悬臂梁结构^[14]，可以建立如下的极限状态函数：

$g = g\left( {{{w}},{{L}},{{E}},{{I}}} \right) = {\varDelta _{\max }} - \frac{{w{L^4}}}{{8EI}}.$

(19)

式中： $w$ 为横向分布载荷， $L$ 为悬臂梁的长度， $E$ 为悬臂梁的弹性模量， $I$ 为悬臂梁的惯性矩. 基本随机变量服从正态分布且相互独立，分布参数见表3.

采用本文方法，取均匀映射参数 $m = 3.5$ ，抽取40个样本点作为训练样本来训练建立SVR代理模型，通过SRS方法抽取100个样本点作为测试样本. 为了便于比较，SVR参数固定，分别取 $C = 10,g = 0.01$ . 在得到SVR代理模型后，采用相同的MCS样本计算失效概率，样本量为10⁶. 计算结果如表4所示.

从表4可以看出，LCVT方法的计算结果与文献[15]的理论值接近，尤其是在小失效概率区域能够更好地拟合结构极限状态函数.

4.3 算例3

算例来自文献[15]，如图5所示为九盒段机翼翼盒，该结构由64个杆单元和42个板单元组成. 由设计准则可以建立如下的极限状态函数：

$g = g\left( {{{L}},{{A}},{{E}},{{P}},{{H}}} \right) = {\varDelta _{\max }} - {\varDelta _y}{\text{ }}{\text{.}}$

(20)

式中： ${\varDelta _{\max }}$ 为机翼翼盒结构的最大许用位移， ${\varDelta _y}$ 为结构在 $y$ 方向上的实际位移， $L$ 为结构的基本长度单位， $A$ 为杆单元的横截面积， $E$ 为杆单元和板单元的弹性模量， $P$ 为机翼翼盒结构的外部载荷， $H$ 为板单元的厚度. 基本随机变量服从正态分布且相互独立，分布参数如表5所示.

采用本文所提方法，取均匀映射参数 $m = 3.5$ ，抽取400个样本点作为训练样本来训练建立SVR代理模型，利用SRS方法抽取100个样本点作为测试样本. 为了便于比较，SVR参数固定，分别取 $C = 15,g = 0.1$ . 在得到SVR代理模型后，采用相同的MCS样本计算失效概率，样本量为10⁶. 计算结果如表6所示.

从表4、6可以看出，采用LCVT抽样方法得到的结构可靠性计算结果相较于其他抽样方法更接近MCS计算结果. 在CVT试验设计方法与LHS方法的基础上发展起来的LCVT试验设计方法兼具CVT方法的均匀性与LHS方法一定的随机性，能够更好地反映出参数设计空间的分布情况，从而避免了主观因素的干扰.

5 结　论

(1)采用支持向量机回归建立代理模型，以较小的训练样本构建代理模型，可以用于结构可靠性分析中隐式极限状态函数计算量较大的问题.

(2)利用质心Voronoi网格产生分布均匀的CVT样本，对CVT样本进行拉丁化操作，得到兼具均匀性与随机性的LCVT样本. 在得到LCVT样本后，没有利用常用的逆变换法得到特定概率分布下的训练样本，而是将LCVT样本均匀映射至均值周围 $m$ 倍方差以内. 在同样的抽样数目下，可以得到更多的失效样本.

(3)对结构可靠性分析中常用的抽样方法进行比较研究，包括SRS抽样、LHS抽样、CVT抽样、US抽样以及Sobol、Halton和Hammersly 3种序列抽样方法. 计算结果表明：在相同的抽样样本数及SVR参数下，LCVT抽样方法有着更高的均匀性和鲁棒性，基于LCVT构建的支持向量回归代理模型逼近精度高，可以用于结构可靠性分析.