一种新的矩独立重要性测度分析方法及高效算法

巩祥瑞， 吕震宙， 孙天宇， 张雷雷， 封雷

(1. 西安现代控制技术研究所, 西安 710065； 2. 西北工业大学 航空学院, 西安 710072)

灵敏度分析方法主要研究结构系统的输入变量对输出响应的影响程度[1]，已广泛应用于可靠性工程[2]、危险分析[3-4]、环境科学[5]、计算物理化学[6-7]等领域。灵敏度分析一般分为局部灵敏度分析方法和全局灵敏度分析方法。局部灵敏度分析方法研究模型输入变量的分布参数对结构系统输出响应的影响程度[8]，只能反映输入变量在某一名义值处对输出响应的影响，具有很大的局限性[9]。

全局灵敏度分析，也称为重要性测度分析，研究一个或多个输入变量对模型输出响应的影响程度，能够反映输入变量对结构系统输出响应的平均影响程度，因此，该分析方法在实际工程中得到了更广泛的应用。重要性测度分析方法的核心思想是对影响结构系统输出响应的输入变量进行重要性排序，从而得到对结构系统输出响应影响较大的输入变量，为结构系统的设计改进提供理论依据。近年来，研究人员已经提出了很多重要性测度分析方法。Saltelli 和 Helton提出了无参的方法，但缺乏模型的独立性[9-10]。Sobol[11]和 Iman[12]等提出了基于方差的重要性测度分析方法和相应的求解算法，该方法隐含地假设了方差可以充分描述结构系统输出响应的不确定性。然而，相关学者指出随机变量的任何一阶矩只能反映其部分信息，必然会丢失一些信息[13]。为了分析输入变量在其整个分布区域上对结构系统输出响应的影响程度，Borgonovo[14]、Liu[15]和 Li[16]等又分别提出了矩独立重要性测度分析(Moment-independent Importance Measures analysis，M-IM)方法。

目前，Borgonovo[14]所提出的矩独立重要性测度分析方法应用最为广泛。但在可靠性分析中，结构系统的失效概率通常是工程技术人员所关注的焦点。Li等[16]又提出了基于失效概率的矩独立重要性测度分析方法，该方法能够计算输入变量在其分布区域的任一固定点时，结构系统的无条件失效概率和条件失效概率的平均偏差。然而，在实际工程中将输入变量固定于某一点，这是很难实现的，只能使输入变量在其分布区域的某一区间上缩减变化。

为了解决这一问题，Wei等[17]提出了一种新的输入变量缩减区间的获取方法，该方法的基本思想是：首先产生2个随机数q1和q2，q1和q2服从均匀分布，q1～U(0,1)，q2～U(q1,1)，然后根据输入变量累积分布函数的逆函数F-1(·)，得到输入变量的缩减区间[F-1(q1),F-1(q2)]。该方法的优点是很容易获取输入变量的缩减区间，然而它并没有考虑输入变量缩减区间获取的任意性和等可能性。

本文提出了一种新的基于失效概率的矩独立重要性测度分析方法，并给出了一种更加合理的输入变量缩减区间的获取方法。同时，引入自适应超球重要抽样(Adaptive Radial-Based Importance Sampling，ARBIS)方法[18-19]，构建所提新指标的高效求解算法。应用数值算例和工程算例，验证本文所提新指标的合理性和求解算法的高效性。

1 新的矩独立重要性测度指标

1.1 新指标的定义

EQ(E(IF)-E(IF|Xi∈Ui(Q)))2

(1)

式中：Pf为非条件失效概率；Q为分位数矩阵，用来生成输入变量Xi的缩减区间；当Xi的分布区间缩减到其所有可能的子区间Ui(Q)时，就能够得到条件失效概率Pf|Xi∈Ui(Q)；EQ(·)为对Q求期望；IF为指示函数，当X位于失效域(即g(X)≤0)时，IF(x)=1，否则IF(x)=0。

1.2 输入变量缩减区间获取方法

为了等概率地获取输入变量所有可能的缩减区间，首先定义一个分位数矩阵Hk(3×2)：

(2)

式中：hk(1)和hk(2)为[0,1]区间上的2个随机数，满足0≤hk(1)

Hk(k=1,2,…,N)组成一个(3N×2)的分位数矩阵Q：

(3)

式中：Hk(k=1,2,…,N) 是被等概率获取的，例如P{Hk}=1/N。

为了更好地理解式(3)，可表示为

(4)

式中：q3k-2(1)=0，q3k-2(2)=q3k-1(1)=hk(1)，q3k-1(2)=q3k(1)=hk(2)，q3k(2)=1 (k=1,2,…,N)。

因此，根据分位数矩阵Q就可以定义缩减区间矩阵U为

(5)

对于输入变量Xi，其缩减区间可表示为

(6)

上述是一种新的区间划分技术来获取输入变量在其分布区域上所有可能的缩减区间。本文所提出的这种新的区间获取方法具有2个优点：①输入变量所有可能的缩减区间能够被等可能地获取，从而就能够更加合理地计算当输入变量在其分布区域的缩减区间上变化时对结构系统失效概率的平均影响程度；②全期望公式在这种区间划分技术中成立，基于全期望公式，第2节中就能够将本文提出的新的基于失效概率的矩独立重要性测度指标转化为更加便于计算的基于方差的重要性测度指标。

2 新指标与基于方差的指标的关系

由于1.2节中提出的输入变量缩减区间获取方法中全期望公式成立，就能够推导出新的矩独立重要性测度指标与基于方差的重要性测度指标间的关系，同时得到一种新的基于方差的矩独立重要性测度指标表示方式。

输入变量在缩减区间上的全期望公式表示为

E(IF)=EQ(E(IF|Xi∈U(Q)))

(7)

证明过程详见附录A。

VQ(E(IF|Xi∈Ui(Q)))

(8)

从式(8)中可以看出，新的矩独立重要性测度指标能够很容易地转换为基于方差的重要性测度指标，从而新的矩独立重要性测度指标就能够采用基于方差的重要性测度分析方法来进行求解。

3 新指标的高效求解算法

新的矩独立重要性测度指标传统的计算方法是双层重复抽样蒙特卡罗(Double-Loop-Repeat-Set Monte Carlo，DLRS MC)[17]方法，该方法先用N个样本点计算非条件失效概率，然后在每一个子区间上，再产生N个新的样本点来计算条件失效概率，最后计算非条件失效概率和条件失效概率差异的平均值。

DLRS MC方法可以获得高精度的结果值，但其计算量很大，计算效率较低。本节提出了一种新的求解算法——ARBIS方法，旨在对新的矩独立重要性测度指标进行高效求解。

3.1 自适应超球重要抽样方法

ARBIS方法最初是用来计算失效概率的，而计算新的基于失效概率的矩独立重要性测度指标的核心就是计算非条件失效概率和条件失效概率。因此，本节基于ARBIS方法构建新指标的求解算法。

ARBIS方法[18-19]的基本思想是自适应地寻找到一个落在功能函数安全域内的最大的超球。当抽样点落入超球内部时，就将该抽样点归类为安全点，不必再代入结构系统功能函数进行求解。因此，该方法在计算结构系统输出响应的失效概率时，计算量将大幅度降低。同时， ARBIS方法采用自适应的方法来获取满足条件的最大超球半径，相比于传统的梯度搜索算法，效率更高，稳健性也更好。

任何一个随机变量通过合适的变换方法，都可以转换到标准正态空间中[20-21]。因此，以下步骤都假定在标准正态空间中进行。为了更好地理解ARBIS方法，给出该方法的一些关键步骤。

步骤1设定初始化超球半径β0。 为了确保超球与失效域相交，初始半径β0可通过式(9)来确定：

(9)

步骤2根据随机变量的分布类型，通过转化得到标准正态空间的样本点uk(k=1,2,…,N)。

(10)

图1 自适应策略获取最优化半径Fig.1 Adaptive strategy for obtaining optimal radius

步骤6重复步骤4和步骤5，直到收敛条件被满足[22]，即可得到最优的超球半径βopt。

3.2 采用自适应超球重要抽样方法求解新指标

当输入变量Xi缩减到区间Ui(Q)时，Xi就是一个截断变量，则所有输入变量的联合概率密度函数可表示为

(11)

式中：K为截断系数，可表示为

(12)

对于式(12)的估计，并不需要计算结构系统的功能函数值，其计算量可以忽略。因此，条件失效概率Pf|Xi∈Ui(Q)可重新表示为

(13)

从式(13)中可以看出，条件失效概率Pf|Xi∈Ui(Q)的计算已转化为一个不含截断变量的多模态并联系统问题。经过这种转化，构建ARBIS方法计算新的矩独立重要性测度指标的算法如下：

步骤1自适应搜索获取超球半径βopt，并计算无条件失效概率Pf；同时得到半径为βopt的超球外失效点组成的样本点矩阵B。

步骤2当输入变量Xi缩减到区间Ui(Q)时，在样本点矩阵B中搜索输入变量Xi落在缩减区间Ui(Q)的样本点，记样本点个数为m。

步骤3计算条件失效概率：

(14)

式中：N为采用ARBIS方法的初始样本点个数。

4 算例分析

4.1 数值算例

某一结构系统的功能函数可表示为

(15)

式中：X1、X2和X3为相互独立的输入随机变量，均服从标准正态分布，即Xi～N(0,1),i=1,2,3。

表1 DLRS MC方法和ARBIS方法计算的新的矩独立重要性测度指标值Table 1 New moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

图2 DLRS MC方法和ARBIS方法计算的新的矩独立重要性测度指标的收敛曲线Fig.2 Convergence curves of new moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

4.2 工程算例

在汽车结构中，车桥承载着大部分汽车重量，其通过悬臂与车架相接，将来自车轮的牵引力和制动力，还有侧向力经过悬架传递给车架，起主要承载作用的就是汽车前轴[23]。由于工字型截断梁能够提高抗弯强度，因此前轴通常采用工字结构梁。

图3为前轴的结构示意图。危险截面处的最大正应力为σ=M/Wx，最大切应力为τ=T/Wρ，M为弯矩，T为扭矩，Wx为截面系数，Wρ为极截面系数，且有

(16)

Wρ=0.8bt2+0.4[a3(h-2t)/t]

(17)

考虑前轴结构静强度失效，有以下极限状态方程：

(18)

式中：σS为静强度屈服极限，根据前轴材料特性有σS=460 MPa;g为裕度值。

前轴结构尺寸参数和承受外载看作独立正态随机变量，分布参数如表2所示。

采用DLRS MC方法和ARBIS方法计算的新的矩独立重要性测度指标结果如表3所示，同时也给出了计算功能函数的总次数。可以看出，2种算法计算结果相近。图4给出了DLRS MC方法和ARBIS方法计算结果标准差的收敛曲线。可见，估计值的标准差都比较小，因此计算结果可靠稳健。DLRS MC方法的总次数为446 000，而ARBIS方法计算总次数仅为4 475，因此ARBIS方法的计算效率很高。

2种算法计算的新的矩独立重要性测度指标的排序均相同：t>T>b>a>h>M。这表明当分别固定6个输入随机变量在其各自的分布区域内缩减变化时，t对汽车前轴结构的失效概率影响程度最大，而M对其失效概率的影响几乎可以忽略不计。

图3 汽车前轴结构示意图Fig.3 Schematic of automobile front axle structure

表2 汽车前轴结构各输入变量分布参数Table 2 Distribution parameters of input variables of automobile front axle structure

表3 DLRS MC方法和ARBIS方法计算汽车前轴指标结果Table 3 Results of indices of automobile front axle computed by DLRS MC and ARBIS methods

图4 DLRS MC方法和ARBIS方法计算汽车前轴的收敛曲线Fig.4 Convergence curves of automobile front axle computed by DLRS MC and ARBIS methods

5 结 论

本文提出了一种新的基于失效概率的矩独立重要性测度分析方法来计算输入随机变量对结构系统输出响应失效概率的平均影响程度。

1) 建立了一种新的区间划分技术来等可能地获得输入随机变量所有可能的缩减区间，并给出了相应证明。

2) 引入了自适应超球重要抽样方法来进行新指标的求解，提高了计算效率。

3) 给出了一个数值算例和一个工程算例，说明了本文所提新指标的意义，同时也验证了新算法的高效性。