参数可靠性全局灵敏度高效分析的代理模型法

员婉莹，吕震宙

1.西北工业大学 重庆科创中心，重庆 400000

2.西北工业大学 航空学院，西安 710072

3.西北工业大学 深圳研究院，深圳 518063

灵敏度分析主要研究的是模型的输出不确定性是如何分配到输入不确定性的［1］。常见的灵敏度分析根据所关注的输出统计量的不同，将其分为基于方差的灵敏度分析［2-3］、基于概率密度函数的灵敏度分析［4-5］以及基于可靠性的灵敏度分析［6-7］。可靠性是在充分考虑不确定性因素基础上，对结构安全程度进行量化的指标。精确的不确定性概率分布是进行精细化可靠性分析的基础。结构系统中的不确定性主要分为2 类，即主观不确定性和客观不确定性。主观不确定性主要是认识不足和信息匮乏导致的，客观不确定性是事物内部的固有随机性。如疲劳寿命，大量试验结果表明其服从对数正态分布，但对于缺乏试验数据的情况，其分布参数将无法准确确定，即具有随机性，这个随机性会随着试验样本量的增加而减小。因此，分布参数可靠性全局灵敏度分析［8］也广受关注，其衡量的是分布参数的不确定性对可靠性的影响，通过分布参数可靠性全局灵敏度分析，识别影响可靠性的重要分布参数，可以通过收集更多对应的基本输入变量的样本信息，消除重要分布参数的不确定性，减小分布参数不确定性对可靠性结果的影响。

分布参数可靠性全局灵敏度分析的直接数值模拟法涉及3 层嵌套分析过程，且计算量与不确定性分布参数的维数相关，3 层嵌套分析过程庞大的计算量无法满足工程中快速灵敏度分析的要求。因此，文献［8］提出了代理模型的方法，但由于建立的是失效概率函数的代理模型，其计算过程亦是1 个双层循环过程，且失效概率函数的代理精度会影响分布参数可靠性全局灵敏度分析的准确性。文献［9］提出了Bayes 公式结合代理模型的方法，该方法仅涉及1 个单层分析过程，将分布参数可靠性全局灵敏度分析等价转化为样本状态（失效或安全）识别问题，且仅需1 组无条件样本矩阵，但该方法中的抽样是简单随机抽样，对于小失效概率问题，抽样样本规模庞大，庞大的候选样本池规模会导致代理模型训练过程极其费时。在文献［9］的基础上，本文提出一种分步自适应代理模型结合重要抽样的单层分析法，该方法通过等概率转换引入标准正态辅助变量，实现分布参数嵌入功能函数的目的，进而实现分布参数重要抽样概率密度函数的建立。通过建立低维随机变量空间中功能函数的Kriging 代理模型，实现高维辅助变量及不确定性分布参数变量近似最优抽样概率密度函数的建立，及相应重要抽样样本点的产生；进一步对构造近似最优抽样概率密度函数的低维随机变量空间中功能函数的Kriging 代理模型进行学习，以准确识别重要抽样样本点的状态（失效或安全），高效完成分布参数可靠性全局灵敏度分析。

本文所提方法极大地提高了样本的利用率，减缩了Kriging 代理模型的训练时间。数值算例和工程算例的计算结果表明，在保证计算精度的同时，本文所提方法相对文献［9］方法极大地减缩了Kriging 代理模型训练过程中的候选样本池规模，减少了计算时间。

1 参数可靠性全局灵敏度定义

定义结构功能函数为Y=g(X)，其中：Y 为结构或系统的输出；X=[ X1，X2，…，Xn]为结构的基本随机输入变量；n 表示基本随机输入变量的个数；X 的联合概率密度函数为fX(x|θ)；θ 为基本随机输入变量的分布参数。当基本随机输入变量的样本较少或对其认识有限时，其分布参数亦会存在不确定性，该不确定性可以用区间不确定性描述、模糊不确定性描述或随机不确定性描述。本文主要研究当分布参数为随机不确定性变量时，其不确定性对失效概率分散性的影响程度。当分布参数θ 具有不确定性时，结构的失效概率将不再是一个定值，其是分布参数θ的函数，即

式中：PF(θ)为失效概率函数；IFg(x)为失效域指示函数，定义为。

为了衡量不确定性分布参数对失效概率分散性的影响，文献［9］提出了分布参数可靠性全局灵敏度，其定义式为

式中：E(·)表示期望 算子；θi表示第i 个不 确定性分布参数；θ-i表示除了θi变量之外的其他不确定性分布参数变量的向量。

式（2）表明，直接模拟法求解第i 个不确定性分布参数θi的可靠性全局灵敏度是一个3 层耦合分析过程，最外层是中间层计算结果|Eθ(PF(θ))-Eθ-i(PF(θ)|θi)|关于θi的期望，中间层的求解关键是Eθ-i(PF(θ)|θi)，其是最内层失效概率函数在θi固定时关于θ-i的期望，最内层是求解不同分布参数取值下的失效概率。3 层法求解分布参数可靠性全局灵敏度指标的真实功能函数调用次数为：m×Nθi×Nθ-i×NX+Nθ×NX，其中：m 表示不确定性分布参数的个数；Nθi表示计算外层期望时θi的样本个数；Nθ-i表示计算中间层期望时θ-i的样本个数；NX表示计算内层失效概率抽取的输入变量X 的样本个数；Nθ为计算失效概率函数无条件均值Eθ(PF(θ))时抽取的不确定性分布参数θ 的样本。

为了克服分布参数可靠性全局灵敏度计算的3 层嵌套过程以及计算量与不确定性分布参数维数的相关性，文献［9］提出了基于Bayes 公式的单层分析法，该方法主要采用了Monte Carlo 简单随机抽样，对于小失效概率问题，需要大量的计算样本。因此，本文将在单层计算公式的基础上，研究基于重要抽样的求解算法，并通过2 步Kriging 代理模型实现近似最优重要抽样概率密度函数的构造以及重要抽样样本状态（失效或安全）的识别。

2 参数可靠性全局灵敏度求解的单层Monte Carlo 法

首先，根据等概率转换［10-11］建立标准正态变量与任意随机变量及其分布参数的关系，即

式中：ui表示标准正态变量；θXi表示Xi变量的分布参数；Φ(·)为标准正态变量的累积分布函数；FXi(·)为Xi变量的累积分布函数。

根据式（3）可以得到随机变量Xi与其分布参数θXi及ui之间的关系，即

将式（4）代入功能函数Y=g(X)中，得到式（5）所示的等效功能函数为

式中：等效功能函数中随机变量的维数为n+m；n 为标准正态辅助变量的个数；m 为不确定性分布参数的个数。

基于等效功能函数gG(u，θ)可得Eθ(PF(θ))和Eθ-i(PF(θ)|θi)的等价计算式如式（6）及式（7）所示：

式 中：失 效 域F 的定义为F={u，θ：gG(u，θ)≤0}； IFG(u，θ) 为 失 效 域 指 示 函 数，定 义 为为n 个相互独立的标准正态变量的联合概率密度函数；fθ(θ)为不确定性分布参数的联合概率密度函数；P{·}为概率算子。

根据Bayes 公式可建立P{F|θi}与P{F}的关系，即

求解式（9）的关键是计算P{F}和fθi(θi|F)，根据Monte Carlo模拟法（Monte Carlo Simulation， MCS）抽取变量u和θ的样本得到Su，θ样本矩阵，即

式中：NMCS表示MCS 样本规模，将样本矩阵Su，θ中样本代入等效功能函数中得到P{F}的估计式为

fθi(θi|F)可以根据Su，θ中的失效样本结合概率密度函数估计方法［12-14］求解，无需额外功能函数的计算。因此，基于Bayes 公式的单层样本法的计算量为NMCS，其与不确定性分布参数的维数无关，由于P{F}的求解也即识别样本矩阵Su，θ中样本的状态（失效或安全），因此文献［9］亦将自适应Kriging 代理模型嵌入基于Bayes 公式的单层样本法中进一步降低功能函数的调用次数。但对于小失效概率问题，大规模的样本矩阵Su，θ会使得Kriging 代理模型自适应学习过程较为耗时，因此，本文将提出基于重要抽样的单层计算法。

3 所提方法

3. 1 参数可靠性全局灵敏度求解的单层重要抽样法

式（9）的 计算关键是估计P{F}及fθi(θi|F)，本节将介绍如何利用重要抽样法计算P{F}及fθi(θi|F)。首 先，引 入 重 要 抽 样 概 率 密 度 函 数hu，θ(u，θ)，根 据 重 要 抽 样 法 可 得P{F}的 计 算式为

式中：Eh(·)表示在重要抽样概率密度函数抽样下的期望算子。式（12）的估计通过hu，θ(u，θ)抽取变量u 和θ 的样本，得到对应的样本矩阵为

式中：NIS表示重要抽样样本规模。

通过功能函数判断重要样本的失效域指示函数值进而得到P{F}的估计值，即

3.2 节将详细介绍如何采用分步自适应更新的Kriging 代理模型构造重要抽样概率密度函数以及如何准确预测各个重要样本的失效域指示函数值。

通过引入Metropolis-Hastings 准则［15-16］，可以 将样 本 矩 阵 中 服 从hu，θ(u，θ|F)的 样 本 转换 为 服 从fu，θ(u，θ|F)的 样 本，进 而 通 过Edgeworth 级数［14］得到fθi(θi|F)，实 现重复利用中样本信息得到各个不确定性分布参数的可靠性全局灵敏度，具体执行过程如下。

第2步 定义初始化样本z(k)=，其中k=1。

第3步 计算接受概率并确定z(k+1)。根据Metropolis-Hastings 准 则［15-16］，接 受 概 率 的 定义为

式中：π(·)是目标分布；p(·)是建议分布。在本文中 目 标 分 布 是 fu，θ(u，θ|F)，建 议 分 布 是hu，θ(u，θ|F)。fu，θ(u，θ|F)和hu，θ(u，θ|F)又 可 分别表示为

式中：P { F }=∫ fu，θ(u，θ)IFG(u，θ)dudθ；Ph ｛F｝=∫hu，θ(u，θ)IFG(u，θ)dudθ。

根 据 式（16）构 建 服 从hu，θ(u，θ|F)的 样 本d(k+1)={uhF(k+1)，θhF(k+1)}向z(k+1)转化的接受概率为

根据式（19）的判据，可以得到z(k+1)，即

式中：random [0，1]表示[0，1]区间均匀分布的随机数。

根据以上流程即可得到服从fu，θ(u，θ|F)的样 本，将 样 本 放 入 矩 阵中，即

式 中：H3(·)、H4(·) 和H6(·) 是3 阶、4 阶、6 阶Chebychev-Hermit 多项式。

综上，可以通过重复利用重要抽样样本点得到各个不确定性分布参数的可靠性全局灵敏度指标。重要抽样法中的关键是如何高效地构造最优抽样概率密度函数和准确识别重要抽样样本的失效或安全状态。3.2 节将建立通过低维自适应Kriging 代理模型构建高维重要抽样概率密度函数的方法以及如何在构建高维重要抽样概率密度函数的低维Kriging 代理模型的基础上通过继续自学习实现重要抽样样本安全或失效状态的识别。

3. 2 自适应代理模型结合单层重要抽样法

3.2.1 基于代理模型的近似最优抽样概率密度函数

重要抽样概率密度函数hu，θ(u，θ)的最优形式为

由于P { F }是个待求量，因此，文献［17］提出了基于Kriging 代理模型的近似最优重要抽样概率密度函数构造方法，即

式中：概率分类函数π(u，θ)表示根据当前Kriging 代理模型gG(K)(u，θ)得到样本点{u，θ }处预测值gG(K)(u，θ)≤0 的概率，定义为π(u，θ)=P { gG(K)(u，θ)≤0}=

式 中：μgG(K)(u，θ) 表 示 Kriging 代 理 模 型gG(K)(u，θ)的均 值［18-19］；σgG(K)(u，θ)表 示Kriging代理模型gG(K)(u，θ)的标准差［18-19］。

从式（27）中可以看出，为了构造概率分类函数π(u，θ)，需要建立n+m 维的Kriging 代理模型，相对建立真实功能函数的n 维代理模型g(K)(X)，gG(K)(u，θ)的建立增加了代理模型的维度。因此，基于等概率转换，本文提出利用n 维Kriging 代理模型g(K)(X)来构造n+m 维概率分类函数的方法，进而产生u 和θ 的重要样本，具体推导如下。

式中：X=[ X1，X2，…，Xn]；；μg(K)(X)表示Kriging 代理模型g(K)(X)的均值；σg(K)(X)表示Kriging代理模型g(K)(X)的标准差。

通过式（28）的等式关系，可以通过构造g(K)(X)结合等概率转换得到变量{u，θ }的概率分类函数。具体流程如下：

第1 步 构造功能函数的初始Kriging 代理模型，将功能函数g(X)在该阶段的Kriging 代理模型记为。首先由u 和θ 的概率密度函数fU(u)及fθ(θ)产生变量{u，θ }的少量样本点，并将其放入矩阵中，即

根据式（4）的转换关系可得到对应的X 空间的样本点，即

第2 步 基于π(u，θ)，利用简单拒绝抽样法［21］抽取变量{u，θ }的重要抽样样本。具体过程如下：

第2. 1 步 令q=1；并通过u 和θ 的概率密度函数fU(u)和fθ(θ)产生变量{u，θ }的N 个样本，计算这N 个样本的概率分类函数值，从而确定简单拒绝抽样中的常数 c，即 c=，其中π(u(i)，θ(i))通过式（28）确定。

第2. 2 步 通过变量u 和θ 的概率密度函数 fU(u) 和 fθ(θ) 产 生 变 量{u，θ } 的 样 本{u(q)，θ(q)}。

第2. 3 步 产生标准均匀分布变量的样本p(q)。

第2. 4 步 当p(q)-cπ(u(q)，θ(q))≤0 接受 样本{u(q)，θ(q)}，将 样 本{u(q)，θ(q)}放入重要抽 样 样本矩阵的第q 行，并执行第2.5 步；否则，返回第2.2 步。

第2. 5 步 当q=NIS时结束简单拒接抽样过程，输出重要抽样样本矩阵；否则，令q=q+1，并返回第2.2 步。

第3步 更新训练样本集。对重要抽样样本进行K-means 聚类分析以覆盖所有可能的失效区域，得到K 个形心

第4 步 判断代理模型的收敛性。由交叉验证法［22］计算修正因子αcorr的留一法估计值，的计 算式为

式中：NT表示训练样本集T 中训练样本个数；；IFg(x(i))=表示除去训练样本集T 中第i 个训练样本后构建的Kriging 代理模型所建立的概率分类函数。当概率分类函数精确近似失效域指示函数时，留一法交叉验证误差应为1，此时需要较多的更新训练样本点来更新Kriging代理模型，文献［17］提出了以下折中的停止准则，即当收敛条件满足，输出当前重要抽样样本矩阵，若收敛条件0.1 ＜不 满足，返 回第2 步。

3.2.2 重要抽样样本失效域指示函数值预测的自适应代理模型法

式中：xh(i)表示样本矩阵中的第i 个样本；表示Kriging 代理模型的均值；表示Kriging 代理模型的标准差；Φ(-U(x))表示样本点x 状态（失效或安全）被误判的概率。误判概率越大的点越需将其加入训练样本集中，因此更新训练样本点可根据式（33）选出，即

4 案例分析

4. 1 数值算例

数值算例的功能函数为

式中：输入变量X1，X2和X3均服从正态分布，即Xi～N(μi，1)(i=1，2，3)；μi为Xi的均值，其具有不确定性，均服从均值为4，标准差为1 的正态分布。

表1 为数值算例的分布参数可靠性全局灵敏度分析结果。从表1 的计算结果可以看出，本文所提方法和文献［9］所提方法的计算精度和稳健性随着候选样本池的增加而增加，且随着样本池的增加，2 个方法的计算精度逐渐收敛于双层MCS 采用大样本计算的结果。在基本相同的计算精度和稳健性下，重要抽样样本池的规模是1 000，而文献［9］方法的候选样本池规模是524 288，本文所提方法相对文献［9］中方法降低了99.81% 的样本池规模，降低了78.31%的计算时间。

表1 数值算例的分布参数可靠性全局灵敏度分析结果Table 1 Results of parameter global reliability sensitivity indices of numerical example

4. 2 导弹弹翼

弹翼和舵面是产生气动升力和法向力（垂直导弹纵轴方向）的重要部件，用以克服重力和实现导弹的机动飞行，保证导弹具有良好的操纵性和稳定性。舵面是指在气流中利用偏转而产生导弹平衡力和控制力来操纵导弹飞行的气动翼面，又称操纵面。按功能可分为升降舵、方向舵和副翼；按作用介质可分为空气舵，燃气舵；按照结构和部位安排分为全动舵、位于翼面后缘的舵、翼尖舵、陀螺舵和扰流片。图1 为某型导弹舵面结构的示意图，其由5 条桁条、6 条翼肋、蒙皮以及作用于和弹体相连接的固定端所组成。各个部件的材料均为某钛合金材料，材料的弹性模量为E，泊松比为ν。翼肋与桁条的厚度均为trs，蒙皮的厚度为tsk。在导弹的飞行过程中，舵面的作用主要用于调节和稳定飞行姿态，其所承受的载荷主要为气动载荷，大小为P，垂直作用于上蒙皮。以舵面结构在承受气动载荷情况下的变形不超过δc=8 mm 为门限值，建立该舵面结构的功能函数为

图1 某型导弹舵面结构示意图Fig. 1 Missile wing structure

式中：max(δ)为结构在承受气动载荷情况下的最大变形，其通过有限元分析得到，确定性有限元分析结果如图2 所示。舵面结构材料的弹性模量、泊松比、桁条、翼肋、蒙皮的厚度以及所承受的气动载荷是相互独立的正态分布随机变量，输入随机变量的均值具有不确定性，分别服从正态分布且相互独立。输入随机变量及其均值的分布参数如表2 及表3 所示。

图2 某型导弹舵面结构有限元模型Fig. 2 Finite element model of missile wing structure

表2 随机变量的分布类型及分布参数Table 2 Distribution types and distribution parameters of variables

表3 随机变量均值的分布类型及分布参数Table 3 Distribution types and distribution parameters of means of variables

图3为不同样本池规模下各个变量分布参数可靠性全局灵敏度分析结果的收敛图，从图3 中可以看出，随着样本池规模的增加本文所提方法的计算稳健性不断提高，且计算结果逐步收敛于大样本容量下文献［9］方法的计算结果，且本文所提方法在样本池规模为5 000 时的计算结果与文献［9］方法在2 097 152 样本池规模下的计算结果的精度和稳健性基本一致。表4 为不同样本池规模下，2 种方法所需的功能函数调用次数和选择更新训练样本点的时间，从表4 的结果中可以看出，在计算精度和稳健性基本一致的情况下，本文所提方法与文献［9］所提方法所需的功能函数调用次数相同，但本文所提方法由于采用重要抽样，提高了样本的抽样效率，降低了候选样本池规模，从而降低了在候选样本池迭代选择更新训练样本点的时间。在计算精度和稳健性基本一致的情况下，本文所提方法相对文献［9］方法降低了91.05%的Kriging 代理训练过程中选择更新训练样本点的时间。2 种方法得到的分布参数可靠性全局灵敏度排序一致，如图4 所示，分布参数可靠性全局灵敏度排序为：μE＞μP＞μtsk＞μtrs＞μν。因此，可以通过获取更多弹性模量E 的数据，提高对弹性模量E 统计规律的认识，进而最大程度地降低失效概率的分散性。

图3 分布参数可靠性全局灵敏度计算结果随候选样本池规模变化的收敛图Fig. 3 Convergence plot of parameter reliability sensitivity indices with changes in the size of the candidate sample pool

图4 分布参数可靠性全局灵敏度排序图Fig. 4 Importance rank of parameter global reliability sensitivity indices

表4 所提方法与文献［9］方法不同样本池规模的计算量Table 4 Computational cost of proposed method and method in Ref.［9］ with different sizes of candidate sampling pool

5 结 论

1）本文通过等概率转换，构建了包含不确定性分布参数的等效功能函数，为建立分布参数可靠性全局灵敏度分析的重要抽样法建立基础。

2）本文提出了分布参数可靠性全局灵敏度分析的分步自适应代理模型的单层重要抽样法，实现了利用原始空间功能函数代理模型构建辅助变量与不确定性分布参数的重要抽样概率密度函数，并通过进一步更新构造最优重要抽样概率密度函数的原始空间功能函数，实现辅助变量与分布参数重要样本状态的准确识别，完成参数可靠性全局灵敏度分析。

3）该方法仅需1 组无条件重要抽样样本即可实现各个不确定性分布参数可靠性全局灵敏度分析。

4）算例分析结果表明，本文所提方法相对文献［9］中方法节省了＞70%的代理模型训练时间，说明了本文所提方法的高效性。