结构可靠性分析的主要任务在于计算结构失效概率^[1]，迄今为止，研究者们已经提出了许多高效的可靠性计算方法。基于设计点的方法^[2-4]是最为常见的一种可靠性分析方法，其中，一次二阶矩方法^[4]最具有代表性，它通过将极限状态函数在设计点处进行泰勒展开，保留展开式一阶项，从而实现对极限状态函数的线性近似。但该方法需要计算极限状态函数在设计点处的梯度，且此方法的准确度严重依赖于所研究问题的线性特征。对于具有复杂失效域的问题，其失效域往往是数量众多的不连通域，且每个失效域相对于整个取值域占比较小，失效域的个数难以预先得到，设计点也难以预先求解。因此该方法无法解决复杂失效域、多设计点等问题。数字模拟方法中最基本的蒙特卡罗模拟(Monte Carlo Simulation, MCS)方法通过抽样并判断样本所处状态，计算结构失效概率。对于实际工程问题，失效概率通常较小(10^-3量级及更小)，导致MCS方法往往需要大量的样本(一般为10²~10⁴/P_f，P_f为结构真实失效概率)，才能获得收敛解，计算效率低下。因此，研究者们提出了一系列改进数字模拟方法。其中，重要抽样(Importance Sampling, IS)法^[5-8]是一种应用广泛的方差缩减方法，它通过构建重要抽样密度函数，使所抽取的样本尽可能多地落入失效域中，从而大幅缩小抽样样本池规模，提高失效概率的计算效率。

重要抽样方法的关键在于重要抽样密度函数的构建。Melcher^[9]提出了基于设计点的方法，通过将原随机变量概率密度函数的中心平移至设计点处构建重要抽样密度函数。但该方法无法处理多设计点、复杂失效域的问题。吴斌等^[10]提出了扩大方差的方法，这种方法难以确定方差扩大的倍数。Priebe和Marchette^[11]提出了混合重要抽样法，对每个失效模式单独构造重要抽样密度函数后加权混合，有效解决了多失效模式问题。Au和Beek^[12]提出了基于核密度估计的方法，通过马尔可夫链对失效域进行预抽样，获得一定量的失效样本后，通过核密度估计得到重要抽样密度函数；但核密度估计方法计算成本较大，且精度难以保证。Dubourg等^[13]提出了元模型重要抽样(Meta-IS)方法，使用Kriging模型近似理论上最优重要抽样密度函数，同时，将失效概率转化为失效概率估计值和修正系数乘积的形式；但该方法在Kriging构建时收敛标准的物理意义并不明确，且需要大量调用极限状态函数以计算样本的真实响应值。Cadini等^[14]提出了改进的Meta-IS-AK²方法，但该方法通过单高斯模型获得的重要抽样密度函数针对多失效域问题效率较低。Kurtz和Song^[15]提出了交叉熵重要抽样(Cross Entropy Importance Sampling, CE-IS)方法，使用高斯混合模型作为重要抽样函数，通过交叉熵原理，迭代更新高斯混合模型(Gaussian Mixture Model, GMM)的参数，逐步逼近最优重要抽样密度函数。该方法较Meta-IS-AK²获得的重要抽样密度函数有更高的抽样效率，但在GMM的参数计算过程中，需大量调用极限状态函数计算样本点处真实输出响应值，使得最终的计算效率大幅下降。且Kurtz和Song^[15]所提方法以落入失效域的重要样本比率作为收敛准则，没有考虑失效概率估计值的变异系数，对于部分问题，会存在解不稳定，迭代次数过多，甚至无收敛解的问题。

对于实际工程问题，极限状态函数往往是基于有限元模型^[16]的隐式函数，每次调用都会带来大量的计算量，因此，应尽可能减少极限状态函数的调用次数。代理模型^[17-22]可以有效避免极限状态函数多次调用问题，其通过一定的选点准则，筛选出对失效概率计算影响较大的点，添加至模型训练集，构建代理模型预测样本点处的响应值，计算结构失效概率。其中，Kriging模型由于其在可靠性分析领域存在特有的优势受到了最为广泛的关注。Kriging模型假设样本点之间服从高斯随机过程，通过极大似然估计模型参数值。由于高斯过程的特性，Kriging模型不仅可以预测样本点处的响应值，还可以预测样本点处的方差值。因此，Kriging模型有别于其他代理模型，对于样本点的预测具有随机属性，可以提供样本点的统计信息，与可靠性分析通过输入变量统计信息求解输出响应部分统计信息的思想相吻合。基于此，研究者们提出了一系列的学习函数，自适应地筛选出训练样本点，以提高Kriging模型的构建效率。常用的自适应学习函数有期望可行性函数^[23]、U学习函数^[24]、H学习函数^[25]、期望最大学习函数^[26]等。其中，U学习函数综合考虑样本点距极限状态面的距离和样本点的预测误差，筛选出最易错误判断失效状态的点，简单易行，本文选择U函数作为学习函数。通过自适应Kriging(Adaptive Kriging, AK)代理模型方法，来用少量训练样本即可构建高精度的收敛的代理模型替代极限状态函数，显著提高失效概率的计算效率。常用的AK-MCS^[24]依旧需要构建大规模的样本池，对于小失效概率问题，计算效率依旧低下。

本文将CE-IS方法与AK方法相结合，通过少量训练样本，自适应地拟合极限状态面，减少CE-IS的极限状态函数调用次数，同时通过GMM，构建高效的重要抽样函数，大幅缩小AK-MCS方法的样本池规模。相较于现有方法，本文所提方法有以下创新性：

1) 改进CE-IS方法收敛条件，避免冗余迭代，同时增强CE-IS方法的适用性。

2) 通过AK模型拟合极限状态函数，可显著减少CE-IS方法的计算量。

3) 通过GMM近似最优重要抽样密度函数，相较于现有的重要抽样密度函数的构造方法，更为高效。

1 基于交叉熵的重要抽样法 1.1 重要抽样法

结构可靠性分析的关键是失效概率的计算，对于已知结构，其极限状态函数为g(x)，其中，x为随机输入变量，其概率密度函数为fX(x)，结构的失效域为F={x|g(x)≤0}，失效概率P_f即为输入随机变量落入失效域的概率，可由积分求得，即

$ {P_{\rm{f}}} = P\left\{ {\mathit{\boldsymbol{x}}\left| {g\left( \mathit{\boldsymbol{x}} \right) \le 0} \right.} \right\} = \int_{x \in \mathit{\boldsymbol{F}}} {{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}x} $	（1）

式(1)通常无解析解，故常采用MCS方法作为参照解。引入失效域指示函数I_F(x) :

$ {I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right) = \left\{ {\begin{array}{*{20}{l}} 1&{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{F}}}\\ 0&{\mathit{\boldsymbol{x}} \notin \mathit{\boldsymbol{F}}} \end{array}} \right. $	（2）

则通过MCS方法估计结构失效概率的表达式为

$ \begin{array}{l} {{\hat P}_{\rm{f}}} = \int_{x \in \mathit{\boldsymbol{F}}} {{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \int_{{{\boldsymbol{R}}^d}} {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right){f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \\ \;\;\;\;\;\;\;E\left[ {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right)} \right] = \frac{1}{N}\sum\limits_{i = 1}^N {{I_{\rm{F}}}\left( {{\mathit{\boldsymbol{x}}_i}} \right)} \end{array} $

（3）

式中：x_i(i=1, 2, …, N)为N个按fX(x)抽取的MCS样本点; R^d为d维的实数域。

然而，对于小失效概率问题，样本点大量落入安全域内，导致MCS方法抽样效率极低，故需要大量的样本点数才能保证计算结果的收敛。针对这个问题，研究人员提出了IS方法，其核心思想在于构建重要抽样密度函数h_X(x)代替原密度函数fX(x)，提高抽样效率。IS方法得到结构失效概率的估计值${\hat{P}} $_f表达式为

$ \begin{array}{l} {{\hat P}_{\rm{f}}} = \int_{x \in \mathit{\boldsymbol{F}}} {{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \int_{x \in \mathit{\boldsymbol{F}}} {\frac{{{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)}}{{{h_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)}}{h_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \\ \;\;\;\;\;\;{E_{{h_\mathit{\boldsymbol{X}}}}}\left[ {{I_F}\left( \mathit{\boldsymbol{x}} \right) \cdot \frac{{{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)}}{{{h_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)}}} \right] = \\ \;\;\;\;\;\;\frac{1}{N}\sum\limits_{i = 1}^N {{I_{\rm{F}}}\left( {\mathit{\boldsymbol{x}}_i^h} \right)\frac{{{f_\mathit{\boldsymbol{X}}}\left( {\mathit{\boldsymbol{x}}_i^h} \right)}}{{{h_\mathit{\boldsymbol{X}}}\left( {\mathit{\boldsymbol{x}}_i^h} \right)}}} \end{array} $

（4）

式中：h_X(x)为重要抽样密度函数，x_i^h(i=1, 2, …, N)为按h_X(x)抽取的N个重要抽样密度函数的样本。失效概率估计值的方差表达式为

$ \text{Var}\left[ {{{\hat P}_{\rm{f}}}} \right] \approx \frac{1}{{N - 1}}\left[ {\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{f_\mathit{\boldsymbol{X}}^2\left( {\mathit{\boldsymbol{x}}_i^h} \right)}}{{h_\mathit{\boldsymbol{X}}^2\left( {\mathit{\boldsymbol{x}}_i^h} \right)}}{I_{\rm{F}}}\left( {\mathit{\boldsymbol{x}}_i^h} \right)} - \hat P_{\rm{f}}^2} \right] $	（5）

理论上可以推导最优的重要抽样密度函数h_opt(x)为

$ {h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right) = {I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right) \cdot {f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)/{P_{\rm{f}}} $	（6）

以h_opt(x)作为重要抽样密度函数时，失效概率估计值的方差为零。通过对最优重要抽样密度函数进行近似拟合，即可得到抽样效率最高的重要抽样密度函数。

1.2 高斯混合模型

假设随机向量x服从K元高斯混合分布，其概率密度函数p(x)为

$ p\left( \mathit{\boldsymbol{x}} \right) = \sum\limits_{k = 1}^K {{{\rm{ \mathit{ π} }}_k}} \cdot f\left( {\mathit{\boldsymbol{x}}\left| {{\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k}} \right.} \right) $	（7）

式中：K为高斯混合成分个数；f(x|μ_k, Σ_k)为第k个混合成分的混合高斯分布密度函数，μ_k和Σ_k分别为相应均值向量和协方差矩阵; π_k为混为混合成分系数，由密度函数的归一化条件知，π_k应满足:

$ \sum\limits_{k = 1}^K {{\pi _k}} = 1 $	（8）

$ 0 \le {\pi _k} \le 1 $	（9）

式中:混合系数π_k表示随机变量x来自于第k个混合成分的概率。

记二值随机变量z_k(k=1, 2, …, K)表示在给定样本x的情况下，其是否由第k个混合变量生成，z_k的后验概率为

$ P\left( {{z_k} = 1\left| x \right.} \right) = \frac{{{\pi _k} \cdot f\left( {x\left| {{\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k}} \right.} \right)}}{{\sum\limits_{k = 1}^K {{\pi _k}} \cdot f\left( {x\left| {{\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k}} \right.} \right)}} = {\gamma _k}\left( \mathit{\boldsymbol{x}} \right) $

（10）

对于一个未知K元混合高斯分布模型

$ p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right) = \sum\limits_{k = 1}^K {{\pi _k}} \cdot f\left( {x\left| {{\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k}} \right.} \right) $	（11）

其未知参数v共3×K组，每一元高斯分布的参数为π_k, μ_k, Σ_k，则

$ \mathit{\boldsymbol{v}} = \left\{ {{\pi _k}, {\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k};k = 1, 2, \cdots , K} \right\} $	（12）

通过交叉熵准则，可以迭代求出K元混合高斯分布的参数，使其尽可能地逼近理论最优重要抽样密度函数。

1.3 交叉熵准则

在信息论中，Kullback-Leibler交叉熵(K-L Cross-Entropy，KL-CE) D, 简称交叉熵(Cross-Entropy，CE)常用来度量函数的差异，其定义为

$ D\left( {f\left( \mathit{\boldsymbol{x}} \right), g\left( \mathit{\boldsymbol{x}} \right)} \right) = \int {f\left( \mathit{\boldsymbol{x}} \right)} \cdot \ln \frac{{f\left( \mathit{\boldsymbol{x}} \right)}}{{g\left( \mathit{\boldsymbol{x}} \right)}}{\rm{d}}x $	（13）

由式(13)可知，交叉熵越小，函数f(x)与函数g(x)之间的差异就越小。为了使p(x; v)逼近最优重要抽样密度函数h_opt(x)，计算它们之间的交叉熵:

$ \begin{array}{l} D\left( {{h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right), p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) = \int {{h_{{\rm{opt}}}}} \left( \mathit{\boldsymbol{x}} \right) \cdot \ln \frac{{{h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right)}}{{p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)}}{\rm{d}}\mathit{\boldsymbol{x}} = \\ \;\;\;\;\;\;\int {{h_{{\rm{opt}}}}} \left( \mathit{\boldsymbol{x}} \right) \cdot \ln {h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}x - \int {{h_{{\rm{opt}}}}} \left( \mathit{\boldsymbol{x}} \right) \cdot \\ \;\;\;\;\;\;\ln p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right){\rm{d}}\mathit{\boldsymbol{x}} \end{array} $

（14）

将式(6)代入式(14)可得

$ \begin{array}{*{20}{c}} {D\left( {{h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right), p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) = H\left( {{h_{{\rm{opt}}}}\left( \mathit{\boldsymbol{x}} \right)} \right) - P_{\rm{f}}^{ - 1} \cdot }\\ {\int {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right)} \cdot \ln \left( {p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) \cdot {f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} \end{array} $

（15）

式中：H(h_opt(x))=∫h_opt(x)·lnh_opt(x)dx为最优重要抽样密度函数的信息熵; P_f为结构真实失效概率，对于给定结构，二者均为定值。最小化交叉熵等价于寻找满足式(16)的参数v：

$ \begin{array}{l} \mathit{\boldsymbol{v}} = \arg \max \int {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right)} \cdot \ln \left( {p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) \cdot {f_x}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}} = \\ \;\;\;\;\;\arg \max {E_\mathit{\boldsymbol{X}}}\left[ {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right) \cdot \ln \left( {p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right)} \right] \end{array} $

（16）

式(16)为嵌套形式，通常无解析解，一般采用迭代算法求解，逐步更新GMM参数，直至满足收敛要求。

为了提高计算效率，每一步迭代也采用重要抽样法进行计算，即，使用前一步迭代所得GMM p(x; w)作为当前重要抽样密度函数，w为其参数。式(16)可改写为

$ \begin{array}{l} \mathit{\boldsymbol{v}} = \arg \max \int {{I_{\rm{F}}}} \left( \mathit{\boldsymbol{x}} \right) \cdot \ln \left( {p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) \cdot \frac{{{f_\mathit{\boldsymbol{X}}}\left( \mathit{\boldsymbol{x}} \right)}}{{p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{w}}} \right)}} \cdot \\ \;\;\;\;\;p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{w}}} \right){\rm{d}}\mathit{\boldsymbol{x}} = \arg \max {E_{\mathit{\boldsymbol{X}} \sim p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{w}}} \right)}}\left[ {{I_{\rm{F}}}\left( \mathit{\boldsymbol{x}} \right) \cdot } \right.\\ \;\;\;\;\;\left. {\ln \left( {p\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{v}}} \right)} \right) \cdot W\left( \mathit{\boldsymbol{x}} \right)} \right] \end{array} $

（17）

式中：W(x; w)= $\frac{{{f}_{\boldsymbol{X}}}\left( \boldsymbol{x} \right)}{p\left( \boldsymbol{x;w} \right)} $。用样本均值代替总体均值:

$ \mathit{\boldsymbol{v}} = \arg \max \sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot \ln \left( {p\left( {\mathit{\boldsymbol{x}}_i^p;\mathit{\boldsymbol{v}}} \right)} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) $	（18）

式中：x_i^p是根据p(x; w)抽取的样本。通常情况下，式(18)的右端是关于v的可导函数，其极值在梯度为0处取得，即

$ \begin{array}{l} \mathit{\boldsymbol{v}} = \left\{ {\mathit{\boldsymbol{v}}\left| {\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\nabla _\mathit{\boldsymbol{v}}}\ln \left( {p\left( {\mathit{\boldsymbol{x}}_i^p;\mathit{\boldsymbol{v}}} \right)} \right) \cdot } \right.} \right.\\ \;\;\;\;\;\left. {W\left( {\mathit{\boldsymbol{x}}_i^p} \right) = 0} \right\} \end{array} $

（19）

将GMM的定义式(7)代入式(19)，得到关于模型参数v的方程：

$ \begin{array}{l} \sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\nabla _\mathit{\boldsymbol{v}}}\ln \left( {\sum\limits_{k = 1}^N {{\pi _k}} \cdot } \right.\\ \;\;\;\;\left. {f\left( {\mathit{\boldsymbol{x}}_i^p\left| {{\mathit{\boldsymbol{\mu }}_k}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k}} \right.} \right)} \right) = 0 \end{array} $

（20）

分别对3×K组模型参数求导，以均值向量μ_k(k=1, 2, …, K)为例，可建立方程为

$ \begin{array}{*{20}{c}} {\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot }\\ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k^{ - 1}\left( {\mathit{\boldsymbol{x}}_i^p - {\mathit{\boldsymbol{\mu }}_k}} \right) = 0} \end{array} $	（21）

解得

$ {\mathit{\boldsymbol{\mu }}_k} = \frac{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot \mathit{\boldsymbol{x}}_i^p}}{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right)}} $

（22）

同理可得Σ_k和π_k(k=1, 2, …, K)的表达式为

$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k} = \frac{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot \left( {\mathit{\boldsymbol{x}}_i^p - {\mathit{\boldsymbol{\mu }}_k}} \right) \cdot {{\left( {\mathit{\boldsymbol{x}}_i^p - {\mathit{\boldsymbol{\mu }}_k}} \right)}^{\rm{T}}}}}{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right)}} $

（23）

$ {{\rm{ \mathit{ π} }}_k} = \frac{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot {\gamma _k}\left( {\mathit{\boldsymbol{x}}_i^p} \right)}}{{\sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {\mathit{\boldsymbol{x}}_i^p} \right) \cdot W\left( {\mathit{\boldsymbol{x}}_i^p} \right)}} $	（24）

式(22)~式(24)即为基于交叉熵原理得到的GMM参数更新准则，通过多次迭代，GMM即可较为准确地近似最优重要抽样密度函数h_opt(x)。但注意到，式(22)~式(24)中失效域指示函数I_F(x)的计算需要调用结构极限状态函数，而多次迭代会使得计算量显著加大，因此，本文将AK方法与CE-IS方法相结合，从而提高算法的效率。

2 CE-IS-AK算法 2.1 Kriging模型

根据Kriging理论，给定包含n₀个训练点的训练集X_t=[x₁^t   x₂^t … x^t_n0]^T及相应的输出响应值集合Y_t=[y₁^t   y₂^t … y^t_n0]^T，x_i^t为d维随机输入变量，Kriging模型可表示为一个多项式与随机过程之和的形式：

$ Y = \mathit{\boldsymbol{f}}{\left( \mathit{\boldsymbol{x}} \right)^{\rm{T}}}\mathit{\boldsymbol{\beta }} + \mathit{\boldsymbol{z}}\left( \mathit{\boldsymbol{x}} \right) $	（25）

式中：f(x)=[f₁(x)   f₂(x) … f_M(x)]^T为M×1维基函数向量; β=[β₁   β₂ … β_M]^T为相应的系数向量；z(x)为零均值高斯随机过程，协方差函数为

$ \text{Cov}\left[ {z\left( {{\mathit{\boldsymbol{x}}_i}} \right), z\left( {{\mathit{\boldsymbol{x}}_j}} \right)} \right] = \sigma _z^2R\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_j};\mathit{\boldsymbol{\theta }}} \right) $	（26）

式中：σ_z²为过程方差；R(x_i, x_j; θ)为相关函数，θ为参数。通过最小二乘法可以确定参数β和σ_z²的估计值分别为${\boldsymbol{\hat{ }\!\!\beta\!\!\text{ }}} $和$\hat{\sigma }_{z}^{2} $：

$ \mathit{\boldsymbol{\widehat \beta }} = {\left( {{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{F}}} \right)^{ - 1}}{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{Y}}_{\rm{t}}} $	（27）

$ \hat \sigma _z^2 = {\left( {{\mathit{\boldsymbol{Y}}_{\rm{t}}} - \mathit{\boldsymbol{F\hat \beta }}} \right)^{\rm{T}}}{{\boldsymbol{R}}^{ - 1}}\left( {{{\boldsymbol{Y}}_{\rm{t}}} - {\boldsymbol{F}}\widehat {\boldsymbol{ \pmb{\mathit{ β}} }}} \right)/{n_0} $	（28）

式中：F为n₀×M维基函数矩阵，F_ij=f_j(x_i^t)(i=1, 2, …, n₀; j=1, 2, …, M); R为n₀×n₀维相关函数矩阵：

$ \begin{array}{l} \mathit{\boldsymbol{R}} = \\ \;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {R\left( {\mathit{\boldsymbol{x}}_1^{\rm{t}}, \mathit{\boldsymbol{x}}_1^{\rm{t}}} \right)}&{R\left( {\mathit{\boldsymbol{x}}_1^{\rm{t}}, \mathit{\boldsymbol{x}}_2^{\rm{t}}} \right)}& \cdots &{R\left( {\mathit{\boldsymbol{x}}_1^{\rm{t}}, \mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}} \right)}\\ {R\left( {\mathit{\boldsymbol{x}}_2^{\rm{t}}, \mathit{\boldsymbol{x}}_1^{\rm{t}}} \right)}&{R\left( {\mathit{\boldsymbol{x}}_2^{\rm{t}}, \mathit{\boldsymbol{x}}_2^{\rm{t}}} \right)}& \cdots &{R\left( {\mathit{\boldsymbol{x}}_2^{\rm{t}}, \mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}} \right)}\\ \vdots & \vdots &{}& \vdots \\ {R\left( {\mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}, \mathit{\boldsymbol{x}}_1^{\rm{t}}} \right)}&{R\left( {\mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}, \mathit{\boldsymbol{x}}_2^{\rm{t}}} \right)}& \cdots &{R\left( {\mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}, \mathit{\boldsymbol{x}}_{{n_0}}^{\rm{t}}} \right)} \end{array}} \right] \end{array} $

（29）

由上述过程可知，Kriging模型的构建即为相关函数参数θ的确定，采取最大似然估计计算最优参数为

$ \mathit{\boldsymbol{\theta }} = \arg \min \left\{ {\sigma _z^2{{\left[ {\text{det}\left( \mathit{\boldsymbol{R}} \right)} \right]}^{1/{n_0}}}} \right\} $	（30）

确定模型参数后，即可通过Kriging模型给出测试样本点x处的预测值${{\mu }_{{\hat{g}}}} $(x)和预测误差$\sigma _{{\hat{g}}}^{2} $(x)：

$ \begin{array}{l} \sigma _{\hat g}^2\left( \mathit{\boldsymbol{x}} \right):\\ {\mu _{\hat g}}\left( \mathit{\boldsymbol{x}} \right) = \mathit{\boldsymbol{f}}{\left( \mathit{\boldsymbol{x}} \right)^{\rm{T}}}\mathit{\boldsymbol{\hat \beta }} + \mathit{\boldsymbol{r}}{\left( \mathit{\boldsymbol{x}} \right)^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\left( {{\mathit{\boldsymbol{Y}}_{\rm{t}}} - \mathit{\boldsymbol{F\hat \beta }}} \right) \end{array} $

（31）

$ \begin{array}{l} \sigma _{\hat g}^2\left( \mathit{\boldsymbol{x}} \right) = \hat \sigma _z^2\left\{ {1 - \mathit{\boldsymbol{r}}{{\left( \mathit{\boldsymbol{x}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}}\left( \mathit{\boldsymbol{x}} \right) + } \right.\\ \;\;\;\;\;\;\;\;{\left[ {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{x}} \right)} \right]^{\rm{T}}}{\left( {{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{F}}} \right)^{ - 1}} \cdot \\ \;\;\;\;\;\;\;\;\left. {\left[ {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{x}} \right)} \right]} \right\} \end{array} $

（32）

式中：r(x)=[R(x, x₁^t)   R(x, x₂^t) … R(x, x^t_n0)]^T为n₀×1维预测点x与训练点集X_t之间的相关系数向量。

2.2 自适应Kriging模型

为了高效地选择模型训练点，通常需采取自适应的方法，选择恰当学习函数与收敛条件。本文选择的U学习函数表达式为

$ U\left( \mathit{\boldsymbol{x}} \right) = \left| {\frac{{{\mu _{\hat g}}\left( \mathit{\boldsymbol{x}} \right)}}{{{\sigma _{\hat g}}\left( \mathit{\boldsymbol{x}} \right)}}} \right| $	（33）

U函数表征了样本输出响应y=g(x)的正负状态被误分的概率。U函数值越小，样本点被错误分类的概率越大，故而筛选出U函数值较小的样本点作为新的训练点：

$ {\mathit{\boldsymbol{x}}_{{\rm{new}}}} = \mathop {\arg }\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{S}}} \min U\left( \mathit{\boldsymbol{x}} \right) $	（34）

式中：S为候选样本池。

U(x)=2时，样本点被正确分类的概率为1-Φ(2)=97.7%，故可以认为当U(x)≥2时，所构建的Kriging模型对样本点的正负号判断有超过97.7%的概率预测正确。因此，确定Kriging模型构建的收敛条件为

$ \mathop {\min }\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{S}}} \left[ {U\left( \mathit{\boldsymbol{x}} \right)} \right] \ge 2 $	（35）

2.3 本文所提算法

本文所提方法结合了CE-IS方法与AK方法，避免了CE-IS方法计算量过大的缺陷，并提出了新的收敛准则，避免了冗余迭代，增强了CE-IS算法的适用性。同时，较AK-MCS方法缩减了样本池，使小失效概率的计算更为高效。以下给出本文所提方法的详细计算步骤。

1) 初始化混合模型参数π_k⁽⁰⁾、μ_k⁽⁰⁾和Σ_k⁽⁰⁾   k=1, 2, …, K。本文设定为等概率混合，即π₁⁽⁰⁾=π₂⁽⁰⁾=…=π_K⁽⁰⁾=1/K，协方差矩阵Σ_k⁽⁰⁾设定为单位矩阵，均值向量μ_k⁽⁰⁾通过空间均匀抽样(本文采用Sobol序列)取得。初始化迭代步数l=0。

2) 根据初始GMM，抽取n₀个样本点作为Kriging模型初始训练点集X_t=[x₁^t   x₂^t … x^t_n0]^T，计算响应值Y_t=[y₁^t   y₂^t … y^t_n0]^T，构建Kriging模型。

3) 使用当前GMM作为重要抽样密度函数$h\left( \boldsymbol{x} \right)=\sum\limits_{k=1}^{K}{\pi _{k}^{\left( l \right)}}\cdot f\left( \boldsymbol{x }\!\!|\!\!\text{ }\!\!\mu\!\!\text{ }_{k}^{\left( l \right)}, \boldsymbol{ }\!\!\Sigma\!\!\text{ }_{k}^{\left( l \right)} \right) $，抽取当前重要样本池S^(l)=[x₁^(l)   x₂^(l) … x_N^(l)]。

4) 用Kriging模型计算样本池S^(l)中样本点U函数值。

5) 判断$\mathop {\min }\limits_{\boldsymbol{x} \in {\boldsymbol{S}^{\left( l \right)}}} $[U(x)]≥2是否满足，如满足，执行第6)步；否则，根据式(34)筛选新训练样本点x_new，计算其真实响应值g(x_new)，并添加至训练点集X_t。即，X_t=X_t∪x_new，Y_t=Y_t∪g(x_new)，更新Kriging模型，且n₀=n₀+1，返回第4)步。

6) 用Kriging模型预测S^(l)中样本点响应值，根据式(4)计算结构失效概率估计值$ \hat{P}_{\text{f}}^{\left( l \right)}$和失效概率估计值的方差$\text{Var}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right] $，之后求得失效概率估计值的变异系数：

$ Cov\left[ {\hat P_{\rm{f}}^{\left( l \right)}} \right] = \sqrt {\text{Var}\left[ {\hat P_{\rm{f}}^{\left( l \right)}} \right]} /\left[ {\hat P_{\rm{f}}^{\left( l \right)}} \right] $	（36）

判断$ \text{Var}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right]$≤ε^*(ε^*为预先设定的变异系数上限值)是否满足，如满足，则最终失效概率计算结果为$ {{{\hat{P}}}_{\text{f}}}=\hat{P}_{\text{f}}^{\left( l \right)}$，变异系数为$\text{Var}\left[ {{{\hat{P}}}_{\text{f}}} \right]=\text{Var}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right] $；如不满足，则执行第7)步。

7) 计算抽样效率：

$ \eta = \sum\limits_{i = 1}^N {{I_{\rm{F}}}} \left( {{\mathit{\boldsymbol{x}}_i}} \right)/N $	（37）

判断η是否大于指定值η^*，如满足，则执行第8)步；否则，根据式(22)~式(24)更新GMM参数，增加迭代次数l=l+1，转跳至第3)步。

8) 扩大样本池规模，如：N=10N，转跳至第3)步。

方法流程图如图 1所示。

2.4 一些讨论及说明

文献[15]以其所提ρ分位数作为收敛准则，ρ分位数的概念与抽样效率η概念相同，均表征重要样本落入失效域中的比例；以其作为收敛准则，即以抽样效率η作为收敛准则，存在以下3个主要问题：

1) 对于部分复杂失效域问题，由于理论最优重要抽样密度函数的高度非线性以及不连续性，GMM难以十分准确地逼近理论最优重要抽样密度函数，单纯以抽样效率作为收敛准会造成方法迭代次数过多，甚至不收敛。

2) 对于部分失效概率较大的问题，只需较少的样本点即可得到稳定的解，以抽样效率作为收敛指标会使算法进行多次无意义迭代。

3) 对于部分极小失效概率问题，即使使用重要抽样的方法，也需要较大规模的样本池，仅通过抽样效率难以判断是否得到了稳定的收敛解。

导致这3个问题的主要原因是，在使用数字模拟方法计算失效概率过程中，失效概率估计值的变异系数$\text{Cov}\left[ {{{\hat{P}}}_{\text{f}}} \right] $≤ε^*，是衡量方法是否收敛的最终指标。抽样效率与失效概率估计值的变异系数有一定的关系，例如在样本点数一定的情况下，使用抽样效率较高的重要抽样密度函数，通常能得到更小的失效概率估计值的变异系数。重要抽样法的最终目的是使用尽可能少的样本点数，来获得满足失效概率估计值变异系数要求的结果。

因此，本文采用变异系数阈值ε^*作为方法收敛的最终标准；抽样效率阈值η^*作为参考指标。每次迭代更新得到失效概率估计值${\hat{P}} $_f和其$ \text{Cov}\left[ {{{\hat{P}}}_{\text{f}}} \right]$后，判断$ \text{Cov}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right]$≤ε^*是否满足，如满足，无论抽样效率η≥η^*是否满足，则都认为方法得到了收敛解。如$\text{Cov}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right] $≤ε^*没有被满足，则认为此时的解并不是收敛解，此时再判断不等式η≥η^*是否被满足，若η≥η^*被满足，则认为重要抽样密度函数已经较为高效，不收敛是由于样本点数目不足引起，则此时需扩大样本池规模；若η≥η^*没有被满足，则认为是由于当前抽样效率η太低，众多样本落入安全域内而导致的方法不收敛，此时首先应继续迭代，提高抽样效率。

本文以失效概率估计值的变异系数作为方法收敛准则，可以准确衡量解的稳定性，以及样本池规模是否恰当，同时避免方法出现冗余迭代。

本文构建AK模型用于拟合极限状态函数，因此，每次更新重要样本池时，无需重新构建模型，仅需在原有Kriging模型的基础上，通过学习函数添加少量训练样本点，修正现有模型，即可完成收敛。

根据文献[15]的研究，本文建议高斯混合模型中混合成分数目K的取值应满足如下要求：

$ \max \left\{ {d, {N_{{\rm{component}}}}} \right\} \le K \le 5\% \eta {\rm{N}}{N_{{\rm{component}}}} $	（38）

式中：N_component为结构失效模式数目，对于单模式问题N_component=1。建议迭代样本池规模N为10³~10⁵。

对于复杂失效问题，往往难以预先判断其失效域个数，因此建议选择d < K≤100d，以尽可能保证每个失效域内都能抽取重要样本。

3 算例分析

本节给出5个算例证明所提算法的高效与适用性，第1个算例为小失效概率问题；第2个算例为串联系统问题，具有多失效域；第3个算例为复杂多失效域问题，不连续失效域较多，失效边界复杂；第4个算例为一个动态非线性振子问题；第5个算例为工程有限元算例。本文将所提方法与AK-MCS和CE-IS等方法所得结果进行对比。

3.1 算例1

一个非线性单设计点问题的极限状态函数为

$ g\left( {{x_1}, {x_2}} \right) = 0.5{\left( {{x_1} - 2} \right)^2} - 1.5{\left( {{x_2} - 5} \right)^3} - 3 $	（39）

输入随机变量为标准正态变量，图 2给出了算例一通过CE-IS-AK方法最终所得的重要样本点和极限状态面以及真实极限状态面。表 1给出了算例一CE-IS-AK方法具体迭代结果，表 2给出了算例一的各方法计算对比结果，表中，N_call为极限状态函数调用次数，N_condidate为备选样本点总数。

针对此算例，本文选择混合成分的个数K=20，每次迭代的样本池规模为N=10³，收敛条件为ε^*=5%，η^*=80%。从算例1结果可看出，本文所提方法具有较高的效率，这是由于GMM模型本身的特性，可以对不连续的多失效域进行近似。通过3次迭代，方法即收敛。相较于CE-IS方法，本文所提方法在大幅减小极限状态函数调用次数的情况下，依旧保持了与其相当的高抽样效率；与AK-MCS方法相比，本文在极限状态函数调用次数基本相当的情况下，大幅度缩减了样本池规模，在实际计算中，大幅提升了计算效率。

3.2 算例2

串联系统具有4个失效域，其极限状态函数为

$ \begin{array}{l} g\left( {{x_1}, {x_2}} \right) = \\ \;\;\;\;\;\;\;\min \left\{ {3 + {{\left( {{x_1} - {x_2}} \right)}^2}/10 - \left( {{x_1} + {x_2}} \right)/\sqrt 2 } \right., \\ \;\;\;\;\;\;\;3 + {\left( {{x_1} - {x_2}} \right)^2}/10 + \left( {{x_1} + {x_2}} \right)/\sqrt 2 , \\ \;\;\;\;\;\;\;\left. {\left( {{x_1} - {x_2}} \right) + 7/\sqrt 2 , - \left( {{x_1} - {x_2}} \right) + 7/\sqrt 2 } \right\} \end{array} $

（40）

式中：输入随机变量均为标准正态变量。

图 3给出了算例2通过CE-IS-AK方法最终所得的重要样本点、极限状态面以及真实极限状态面。表 3给出了算例2 CE-IS-AK方法具体迭代结果，表 4给出了各方法计算对比结果。

针对此算例，本文选择混合成分的个数K=50，每次迭代的样本池规模为N=10³，ε^*=5%，迭代门限值为η^*=80%。

CE-IS方法只考虑抽样效率，故出现了冗余迭代的情况。而由于本文提出了更加合理的收敛条件，在失效概率估计值的变异系数满足要求后即收敛，虽然最终得到的抽样效率略低于CE-IS方法，但求解精度并无显著区别，且计算量有所减小，若再进行一次迭代即可得到与CE-IS方法相近的求解情况。与AK-MCS方法相比，本文所提方法虽然调用极限状态函数次数有所增加，但样本池规模远远小于AK-MCS方法，避免了对非常大一部分安全域内的样本点进行无意义的Kriging模型的反复多次预测。

3.3 算例3

一经典的复杂高度非线性多失效域问题的极限状态函数为

$ g\left( {{x_1}, {x_2}} \right) = 10 - \sum\limits_{i = 1}^2 {\left( {x_i^2 - 5\cos \left( {2{\rm{ \mathit{ π} }}{x_i}} \right)} \right)} $	（41）

式中：输入随机变量均为标准正态变量。表 5给出了算例3的各种方法的计算结果。图 4给出通过所提方法最终所得的重要样本点、近似极限状态面以及真实极限状态面。

针对算例3，本文选择混合成分的个数K=100，每次迭代的样本池规模为N=10³，ε^*=5%，η^*=80%。此算例为一高度非线性多失效域算例，失效域数量较多，面积均较小且较为分散，理论上难以获得抽样效率较高的重要抽样密度函数，但此算例失效概率较大，因此并不需要特别大量的重要样本点即可得到收敛解。CE-IS方法由于仅以抽样效率作为收敛指标，故没有得到收敛的解。本文所提方法虽然最终抽样效率较之前算例有明显下降，但依旧得到了较为准确的收敛解，同时，抽样效率仍为所有对比方法中最高且样本池规模最小。故而，相对于其他对比方法而言，仍然具有较大优势。

3.4 算例4

一单自由度无阻尼弹簧振子模型，系统结构示意图见图 5，系统的极限状态函数为

$ g = 3r - \left| {{Z_{\max }}} \right| = 3r - \left| {\frac{{2{F_1}}}{{m\omega _0^2}}\sin \left( {\frac{{{\omega _0}{t_1}}}{2}} \right)} \right| $	（42）

式中：$ {{\omega }_{0}}=\sqrt{\frac{{{c}_{1}}+{{c}_{2}}}{m}}$, c₁和c₂为弹簧刚度系数，ω₀为谐振频率，m为质量；r为屈服极限；F₁为外加激振力；Z_max为位移；t₁为对应时刻。表 6中给出了六维独立输入随机变量分布参数，表 7中给出了此算例采用不同方法的计算结果。

针对此算例，选择混合成分的个数K=100，每次迭代的样本池规模为N=10³，ε^*=5%，η^*=80%。此算例的各方法对比结果证明了本文所提算法对较高维数的问题依然可以高效处理。

3.5 算例5

一个二维桁架结构，如图 6所示。该结构具有23根杆，13个节点，几何形状确定，材料属性和载荷X=[A₁   A₂   E₁   E₂   P₁   P₂ … P₆]^T具有随机不确定性，其中，[A₁   E₁^T]为横杆的横截面积和杨氏模量，[A₂   E₂^T为斜杆的横截面积和杨氏模量，[P₁   P₂ … P₆^T]为节点载荷，u为位移，表 8中给出了输入变量具体的分布形式。

以桁架中点位移u(X)为研究对象，结构极限状态函数为：

$ g\left( \mathit{\boldsymbol{X}} \right) = \tau - u\left( \mathit{\boldsymbol{X}} \right) $	（43）

式中: τ=12 cm为位移u(X)的阈值。

相应地，结构失效域为F={X|g(X)≤0}。算例5无法通过结构力学知识解析得到u(X)的具体表达式，因此，采用有限元仿真对其进行计算。

本文采用规模为1×10⁶的MCS样本池，虽然在计算过程中需大量调用有限元模型，但由于该结构较为简单，有限元模型计算迅速，因此，仍可通过MCS方法得到参照解。表 9中给出了该工程算例的计算结果。同时，各方法均由同一计算机进行计算验证，其计算时间T一并列出。

针对此算例，选择混合成分的个数K=100，每次迭代的样本池规模为N=10³，ε^*=5%，η^*=80%。

3.6 一些讨论及说明

文中Kriging模型采用Matlab中的Dace工具包构建，该工具包自适应地优化求解参数θ。Dace工具包需人工设置如下参数：回归项类型、相关函数类型、θ初始值、θ下界和θ上界。本文的统一设置为：0阶回归多项式、高斯型相关函数、θ初值为0.1×[1 1 … 1]_d^T、θ下界为10^{- 6}×[1 1 … 1]_d^T、θ上界为10⁶×[1 1 … 1]_d^T, 其中，[1 1 … 1]_d^T为d维向量，每一维均为1，d为输入变量的个数。

关于抽样阈值η^*的选取，从2.3节的分析可以看出，过小的η^*会使算法提前扩大样本池规模，增大了计算量，与重要抽样法缩减样本池规模的初衷相违背。过大的η^*有可能会使抽样效率指标失去参考意义(如取η^*=100%，当$\text{Cov}\left[ \hat{P}_{\text{f}}^{\left( l \right)} \right] $≤ε^*不满足要求时，算法会自动更新重要抽样密度函数，若初始样本池规模选择过小，算法还会存在无法收敛的风险)。因此，本文结合文献[15]所提收敛指标，选取η^*=80%。

对于算例2，其在第3次迭代(l=2)时，虽然抽样效率只有η=76%，但其失效概率估计值的变异系数$\text{Cov}\left[ {{{\hat{P}}}_{\text{f}}} \right] $=3.66% < 5%=ε^*，已经收敛。作为比较，增加了一次重要抽样密度函数的迭代，从结果可以看出，虽然迭代后抽样效率有所增加，但计算结果并无显著提升，可认为此次迭代为冗余迭代。对于算例3和算例5，其最终计算结果，抽样效率均未达到η=80%，但失效概率估计值的变异系数均满足阈值要求，因此，得到了失效概率的收敛解。与MCS所得参照解的对比也可看出，本文所提方法求得的失效概率是满足精度要求的。

本文所提方法是重要抽样与代理模型方法的结合，所提方法采用高斯混合模型作为重要抽样密度函数来抽取样本，该方法在参数迭代过程中(式(22)~式(24))计算W(x; w)= $ \frac{{{f}_{\boldsymbol{X}}}\left( \boldsymbol{x} \right)}{p\left( \boldsymbol{x;w} \right)}$时，需要用到输入变量的原始密度函数fX(x)，但其对fX(x)的类型并无限制。因此所提方法是适用于非正态变量的情况的。算例5中变量分布类型为对数正态分布和耿贝尔分布，其计算结果的正确性也证明了本文所提方法是适用于非正态分布情况的。

本文所提算法可高效处理复杂失效域与小失效概率耦合问题，即结构系统的失效域较复杂且同时失效概率也较小的问题。

对于小失效概率问题，常规的数字模拟方法，需要规模巨大的样本池来保证计算的精度。所提方法使用基于高斯模型的重要抽样法，大幅度缩减了样本池规模。同时，引入自适应Kriging代理模型辅助计算，减少了功能函数的调用次数，进一步提高了计算效率。算例1所示结果也证明了本文所提方法能在精度范围要求内高效处理小失效概率问题。

对于复杂失效域问题，往往是串并联系统问题、隐式函数问题，甚至是有限元模型问题，针对此类问题，常规代理模型方法虽然能减少计算成本巨大的功能函数的调用次数，但由于样本池规模所限，计算时间成本依旧较大。本文所提方法引入基于高斯混合模型的重要抽样法，极大地缩减了样本池的规模，使计算效率进一步提高。算例2、算例3、算例4和算例5的计算结果也验证了本文所提方法的高效性。

4 结论

1) 本文针对复杂失效域和小失效概率耦合的失效概率计算问题，提出了交叉熵重要抽样结合自适应Kriging模型的CE-IS-AK方法。通过交叉熵准则指导高斯混合模型逐步逼近理论最优重要抽样密度函数，参数更新中使用AK模型近似计算样本响应值，从而大大减少了计算过程中的结构极限状态函数调用次数，在保证方法计算精度的同时提高了CE-IS方法的计算效率。

2) 本文改进了CE-IS方法的收敛准则，以失效概率估计值的变异系数作为最终指标，抽样效率作为参考指标，避免了冗余迭代，也克服了CE-IS方法对于部分复杂失效域问题无法得到收敛解的缺点，扩大了方法的适用范围。

3) 相较于CE-IS方法，由于AK模型的引入，结构极限状态函数的调用次数大幅缩减，计算效率显著提高。相较于AK-MCS方法，由于采用了高斯混合模型代替重要抽样密度函数，抽样效率较MCS方法得到了显著提高，样本池规模大幅缩小，适用于小失效概率问题。同时，由于利用了高斯混合模型，方法对于多失效域和复杂失效域等问题也有良好的适用性。最后，算例分析结果也证明了本文所提方法的高效性。