基于证据理论和代理模型的QMU分析

谢朝阳，李贵杰，彭忠明，黄洪钟

(1. 中国工程物理研究院总体工程研究所 四川 绵阳 621000；2. 电子科技大学机械电子工程学院 成都 611731)

在结构可靠性以及系统性能分析中，需要通过不确定性传播研究输入参数的不确定性对系统响应的影响，获得系统响应的不确定性。工程实际中随机和认知不确定性同时存在，设结构系统的响应函数抽象为：

不确定性普遍存在于工程结构的分析和设计中。根据不确定性的不同来源和人们的认知程度，在工程实际研究中将不确定性分为随机不确定性(aleatory uncertainty)和认知不确定性(epistemic uncertainty)[1-2]。传统的结构可靠性分析主要考虑随机不确定性，一般基于概率方法，需要大量实验样本构造参数的精确概率分布[3-4]。对于信息难以获取或掌握信息较少的情况，国内外学者提出了基于模糊数学、区间分析、证据理论等非概率模型来描述处理结构可靠性分析中的认知不确定性，并发展了相应的非概率可靠性指标和分析方法[5-8]。近年来，美国圣地亚等国家实验室提出了裕量与不确定度量化(quantification of margins and uncertainty, QMU)的概念用于具有随机和认知不确定性的禁试武库可靠性、安全性评估认证以及风险决策[9-10]。国内外学者围绕随机和认知不确定性分别研究了QMU理论及其应用于系统性能评估和分析时的具体实施流程。文献[9]描述了如何采用QMU处理复杂系统的可靠性以及代替传统可靠性估计的置信度，提出了QMU方法的关键要素即性能通道、裕量、不确定性。文献[11]提出的QMU度量指标考虑了系统操作区域和性能要求的不确定性。文献[12]应用QMU方法研究了环状结构的可靠性。针对混合不确定性，文献[13]集成概率区间分析和证据理论研究了可应用于QMU的计算工具。文献[14]基于Bayes网络充分考虑复杂分层系统各级组件的随机和认知不确定性，建立了QMU方法的Bayes网络框架。

从上述文献看，不确定性量化与传播是QMU方法应用于工程实际的关键环节，而实际工程中系统性能响应常常由隐式模型表示，需进行复杂的数值模拟计算，因此需要建立代理模型进行不确定性传播分析从而降低计算成本。本文在QMU概念的基础上，考虑随机和认知混合不确定性，提出了基于证据理论和Kriging代理模型的结构性能裕量与不确定性量化计算方法，以提高不确定性分析精度与效率。

1 QMU的基本概念与分析流程

QMU是一套考虑不确定性并以性能阈值和裕量为关键点的复杂系统风险决策支持方法，如图1所示，其裕量(M)的定义为系统名义响应值和性能通道边界阈值之间的距离[9-10]。系统设计或服役结构的性能是否满足要求通过置信因子CF(confidence factor)进行定义度量。置信因子CF定义为性能裕量M与该性能的不确定度U之比：

如果对系统评估过程中的不确定性认识量化足够充分，当置信因子CF＞1时，就判定为可靠；当置信因子CF≤1时，就判定系统存在失效风险。

对于复杂系统或大型结构，运用QMU方法评估其性能状态的分析流程涉及3个环节(如图1所示)：1)建立观测清单；2)建立性能通道；3)性能裕量预测及其不确定性量化。对于一般的结构可靠性分析，如果已知结构性能通道特征方程，即性能响应的功能函数，则采用QMU进行评估的关键是通过不确定性传播计算结构响应的不确定性分布。

2 基于证据理论的混合不确定性分析

2.1 证据理论的基本原理

D-S证据理论源于Dempster和Shafer提出的一套基于“证据”和“组合”来处理不确定性推理问题的数学方法。由于证据理论能够对各种不完全信息、不确定信息、不可靠信息甚至冲突信息进行合理的描述和处理，因而被广泛应用于各种领域的不确定性研究。其基本原理主要包括基本可信度分配函数、信任与似然函数[15]。

图1 QMU方法的关键要素流程示意图

设Θ为识别框架，即对于某一问题或事物所能认识到的所有可能结果的集合。函数m为从Θ的幂集2Θ到[0,1]的映射且满足以下3条性质：则函数m为识别框架Θ上的基本概率分配(BPA)。m(A)称为集合A的基本概率分配值，表示对A的信任程度，其中m(A)＞0的集合称为该分配函数的焦元。若B表示识别框架的任一子集，则B的真实程度由信任函数和似然函数表示：

式中，Bel(B)表示完全支持命题B的证据的基本概率函数值之和；Pl(B)是完全或者部分支持命题B的基本概率函数之和。Bel(B)和Pl(B)分别构成对B信任度的下限和上限，记为[Bel(B),Pl(B)]。

2.2 基于证据理论的不确定性传播模型

根据证据理论，设不确定性变量Y的描述用证据空间(Y,Y,m)表示。其中Y表示不确定性变量Y的所有可能值的集合，构成Y的样本空间。令U表示集合Y的子集，Y表示子集U的可数集合，即Y的焦元的集合。m为子集U的基本概率分配函数，满足当U∈Y时，m(U)≥0；当U∉Y时，m(U)=0且对所有U∈Y，

在结构可靠性以及系统性能分析中，需要通过不确定性传播研究输入参数的不确定性对系统响应的影响，获得系统响应的不确定性。工程实际中随机和认知不确定性同时存在，设结构系统的响应函数抽象为：

式中，X为客观不确定性向量，采用概率分布来描述其不确定性；Y为n个证据空间(Y1,Y1,m1),(Y2,Y2,m2),…,(Yn,Yn,mn)描述的不确定性变量，其联合证据空间可表示为(Y,Y,mY)。其中Y=Y1×Y2×…×YnY，U=U1×U2×…×UnY，U∈Y，Uk∈Yk，(k=1,2,…,nY)，mY(U)为联合BPA，可以通过下式计算：

根据全概率理论，G小于某个阈值c的概率P可写为：

式中，乘积空间Y中的焦元Ui的概率等于Ui的联合BPA，即Pr{Yi∈Ui}=mY(Ui)。根据认知不确定性变量的性质，响应G小于某个阈值c的概率将会是一个区间，其边界用信任函数和似然函数度量：

改变不同的阈值c，即可得到系统响应G的不确定性描述。

3 混合不确定性传播求解方法

3.1 Kriging代理模型

式中，fT(x)=[f1(x),f2(x),…,fb(x)]为回归模型的基函数，是一个确定性部分提供对模拟全局的近似，一般用x的多项式表示；α为回归模型对应基函数的系数；S(x)是一个高斯随机过程，提供对模拟局部偏差的近似，其均值为0，方差为σ2。通常情况下多项式f(x)可以取为平均响应值μ，而不影响模拟的精度，式(10)可改写为：

对于空间某点x′，系统性能函数的响应可通过Kriging代理模型的预测值和预测均方差分别表示为：

式中，r为点x′和样本点之间的相关向量；A是n×1的单位向量。

3.2 基于最大置信水平期望提高加点方法

如何确定Kriging代理模型训练样本的位置和数量是保证代理模型计算精度和效率的关键。文献[18]计算结构可靠性时采用自适应抽样方法，采用最大置信水平期望提高准则建立Kriging代理模型训练样本的加点方法。

基于已构建的Kriging代理模型，可通过MCS计算结构的响应不确定性分布。采用MCS方法通过代理模型M计算时，假设通过随机抽样得到N个输入向量Xm=[xm,1,xm,2,…,xm,N]，结构响应函数在点xm,i的预测值G(xm,i,)为服从分布的随机变量。对于样本空间中的第i个样本向量xm,,的结构响应预测值可按其是否小于或大于阈值c进行分类。针对该点的预测置信水平CL(xm,i)可根据正确分类的概率表示如下：

式中，Φ(⋅)表示标准正态分布函数；|⋅|为求绝对值符号。式(14)表示在响应预测值的条件下，实际真实值G(xm,i)＜c的概率。对于由MCS产生的N个样本点，根据阈值c正确分类概率获得的置信水平期望为：

ECL表示代理模型替代原响应函数的近似精度。如果Kriging代理模型M计算响应预测值的近似精度低于设定的目标值，则需添加新的训练样本点构建新的代理模型M*。设响应函数在点xm,i处的真实响应值为G(xm,i)，增加(xm,i，G(xm,i))为训练数据后构建的代理模型M*，其近似精度的置信水平期望提高可表示为：

由于尚未获得样本点xm,i处的真实响应值G(xm,i)，式(16)无法直接计算获得，因此通过置信水平定义EI计算表达式为[18]：

式中，CL(xm,i)是样本xm,i的响应值正确分类的概率；fX(xm,i)是样本xm,i的概率密度函数值；为响应预测值的均方差。式中第一项表示xm,i成为训练样本后代理模型针对该点预测置信水平提高的潜在值，前两项相乘表示增加样本点xm,i后获得置信水平提高潜在值的可能性，第三项表示响应预测值的均方差对置信水平期望提高的影响。通过式(17)计算MCS样本xm中每个点xm,i的置信水平的期望提高，其中使得该值最大的点作为下一个加入Kriging代理模型的训练样本点，即：

3.3 混合不确定性传播计算流程

基于最大置信水平期望提高准则的自动加点Kriging代理模型和证据理论，以求解式(9)的似然函数为例，计算混合不确定性传播的主要步骤如下：

1)根据证据理论，计算响应函数中主观不确定性变量Y的联合BPA，并确定每个焦元Ui的边界。

2)根据随机变量X的分布类型，采用拉丁超立方抽样(Latin Hypercube sample)方法取其中n为随机变量X的维数。

3)采用优化方法计算所抽取样本点在焦元Ui内对应的响应函数的极值Gmin(xt)。

4)根据xt和Gmin(xt)建立Kriging代理模型M。

5)根据MCS抽样和代理模型，采用MCEI准则搜寻点x*并计算Gmin(x*)，将该数据加入到原样本中返回步骤4)，更新代理模型。

6)重复步骤4)、步骤5)，直至满足ECL大于设定的限制。

7)根据MCS样本计算响应函数小于阈值c的概率边界。

8)重复步骤3)～步骤7)步计算出每个焦元Ui内响应函数的输出不确定性小于某一阈值c的概率边界。

9)根据联合BPA和步骤8)步计算结果求解在随机和认知混合不确定性输入下响应函数的不确定性分布。

式(8)表示的信任函数则需将步骤3)～步骤5)中的求最小值Gmin(xt)替换为求最大值Gmax(xt)，即可获得系统响应小于某一特定阈值c的概率下限。

4 基于证据理论框架的置信因子计算

从QMU的概念和流程可知，性能裕量M由系统响应和阈值的估计值描述，不确定性U也由系统响应和阈值的不确定性进行确定。在证据理论的框架下，结构性能响应的不确定性分布用似然函数和信任函数表示。引入安全系数γ，从支撑风险决策的角度保守估计，针对性能通道上边界阈值要求定义系统响应估计值和不确定性取值见表1。下边界阈值要求对应的系统响应估计值和不确定性取值见表2。

基于表1和表2的取值，上下边界的性能裕量M和不确定性U均可由下式计算：

如果性能通道同时有上下边界阈值要求，表征系统QMU度量的置信因子则为上下边界置信因子的最小值。

表1 系统响应估计值和不确定性取值(上边界阈值)

5 算例研究

图2 曲柄滑块机构示意图

图2为曲柄滑块机构示意图[19]。曲柄的长度为a，连杆的长度为b，其内径为d1，外径为d2，受到外力载荷为F，连杆的杨氏模量为E，连杆的屈服强度为S。它们均为随机变量，其分布参数列于表3中。此外，滑块与地面NN的摩擦力系数μ和偏移量e的精确分布不可获知。基于专家经验和有限的历史数据，变量μ和e的区间和BPA可以获得，它们的BPA及其联合BPA分别列于表4和表5中。针对连杆的强度，传统的结构可靠性分析定义极限状态函数为材料的强度与最大应力之间的差，即：

文献[19]基于证据理论采用UAA方法求解了上述极限状态函数的非精确概率可靠性。本文根据QMU方法的框架，以连杆的最大应力作为系统的性能响应函数Y，S作为性能响应函数的上边界阈值，采用置信因子评估结构的可靠性。

表3 曲柄滑块机构中随机变量X的分布参数

表4 曲柄滑块机构算例认知不确定性变量

表5 曲柄滑块机构认知变量的联合BPA

由上述自动加点Kriging代理模型方法计算得到的系统响应函数的不确定性分布由信任函数和似然函数描述如图3所示，其中虚线表示系统响应Y的信任函数和似然函数，实线表示阈值YT的累积概率分布函数。

根据表1和式(10)～式(20)提供的计算方法，取不同的安全系数γ得到的系统QMU分析结果如表6所示。根据QMU的概念，当置信因子CF＞1时，判定为可靠；当置信因子CF≤1时，判定系统存在失效风险。表6中当风险决策时的安全系数γ=0.95时，置信因子CF＜1，说明系统存在失效风险，当可接受的决策风险变大，即安全系数降低时，置信因子CF＞1，说明系统可靠为可接受状态。算例表明安全系数的大小与QMU分析得到的置信因子密切相关，从一定程度上表示了风险决策者对系统决策风险的可接受程度。

图3 系统响应函数和阈值的不确定性分布

表6 不同安全系数取值计算得到的CF

从结构可靠性分析的视角出发，根据蒙特卡洛方法式(21)定义的极限状态方程，可得到响应失效概率的真实值介于信任函数与似然函数值之间，即4.68×10-4≤Pf≤1.164 8×10-3。说明滑块存在较低的失效概率。γ=0.95时的置信因子CF能够反映出系统的失效状态，这与传统可靠性分析的结果相一致。

由于本例中系统响应函数是认知不确定性变量e和μ的单调函数，函数的极值出现在两个变量的端点处。因此不考虑优化求极值过程，采用基于Kriging方法调用系统响应函数的次数为2×9×50=900(每个焦元建立代理模型的训练样本数为50)，而采用MCS方法9个焦元计算信任与似然函数分布共需极限状态函数次数为9×106。这说明本文建立的基于Kriging方法用于随机和认知混合不确定性下的QMU分析具有很好的可行性。

6 结 束 语

本文建立的基于证据理论和代理模型能够有效处理混合不确定性的传播计算，从而支撑系统QMU分析过程中的不确定性量化这一关键要素。算例研究表明，基于最大置信水平期望提高加点的Kriging模型在处理混合不确定性分析方面具有较高的精度和效率。QMU是一套考虑不确定性并以性能阈值和裕量为关键点的复杂系统风险决策支持方法，算例分析表明QMU方法可作为传统可靠性分析的补充，用于度量系统可接受风险时的可靠性分析。

[1]HELTON J C, JOHNSON J D, OBERKAMPF W L, et al.Representation of analysis results involving aleatory and epistemic uncertainty[J]. International Journal of General Systems, 2010, 39(6)： 605-646.

[2]HELTON J C. Alternative representations of epistemic uncertainty[J]. Reliability Engineering and System Safety,2004, 85(1-3)： 1-10.

[3]RACKWITZ R. Reliability analysis-a review and some perspectives[J]. Structural Safety, 2001, 23(4)： 365-395.

[4]AU S K, BECK J L. Estimation of small failure probabilities in high dimensions by subset simulation[J]. Probabilistic Engineering Mechanics, 2001, 16(4)： 263-277.

[5]DU X, VENIGELLA P K, LIU D. Robust mechanism synthesis with random and interval variables[J]. Mechanism& Machine Theory, 2009, 7(7)： 1321-1337.

[6]BEN-HAIM Y. A non-probabilistic concept of reliability[J].Structural Safety, 1994, 14(4)： 227-245.

[7]ELISHAKOFF I, ELISSEEFF P, GLEGG S A L.Nonprobabilistic, convex-theoretic modeling of scatter in material properties[J]. AIAA Journal, 1994, 32(4)： 843-849.

[8]HUANG H Z, WANG Z L, LI Y F, et al. A nonprobabilistic set model of structural reliability based on satisfaction degree of interval[J]. Mechanika, 2011, 8(1)： 85-92.

[9]EARDLEY D. Quantification of margins and uncertainties(QMU)[R/OL].[2016-03-05]. https：//fas.org/irp/agency/dod/jason/margins.pdf.

[10]NAS, NRC. Evaluation of quantification of margins and uncertainties for assessing and certifying the reliability of the nuclear stockpile[R]. Washington： National Academy Press, 2008.

[11]PEPIN J E, RYTHERFORD A C, HEMEZ F M. Define a practical QMU metric[C]//49th AIAA Structures, Structural Dynamics and Materials Conference. [S.l.]： [s.n.], 2008.

[12]LUCAS L J, OWHADI H, ORTIZ M. Rigorous verification, validation, uncertainty quantification and certification through concentration-of-measure inequalities[J].Computer Methods in Applied Mechanics and Engineering,2008, 197(51-52)： 4591-4609.

[13]SENTZ K, FERSON S. Probabilistic bounding analysis in the quantification of margins and uncertainties[J]. Reliability Engineering and System Safety, 2011, 96： 1126-1136.

[14]URBINA A, MAHADEVAN S, PAEZ T L. Quantification of margins and uncertainties of complex systems in the presence of aleatory and epistemic uncertainty[J].Reliability Engineering and System Safety, 2011, 96：1114-1125.

[15]BAE H R, GRANDHI R V, CANFIELD R A. An approximation approach for uncertainty quantification using evidence theory[J]. Reliability Engineering and System Safety, 2004, 86(3)： 215-225.

[16]MARTIN J D, SIMPSON T W. Use of kriging models to approximate deterministic computer models[J]. AIAA Journal, 2005, 43(4)： 853-863.

[17]BANG I K, HAN D S, HAN G J, et al. Structural optimization for a jaw using iterative Kriging metamodels[J].Journal of Mechanical Science and Technology, 2008,22(9)： 1651-1659.

[18]WANG Z, WANG P. A maximum confidence enhancement based sequential sampling scheme for simulation-based design[J]. Journal of Mechanical Design, 2014, 136(2)：1-10.

[19]DU X P. Unified uncertainty analysis by the first order reliability method[J]. Journal of Mechanical Design, 2008,130(9)： 1-10.