一种基于证据理论和随机变量的主动学习克里金方法，用于量化边际值与不确定性

《Applied Mathematical Modelling》：An active learning Kriging method for quantification of margins and uncertainties using evidence theory and aleatory variables

　　凌兆|张辉|张志林|齐阳|索斌|文金鹏•主动学习克里金法可量化混合不确定性下的安全裕度与不确定性。•通过基于基本概率赋值的权重函数来指导样本更新。•高精度评估的次数可减少4.4%至93.6%。•该方法支持多输出故障模式的安全裕度与不确定性评估。•楔环连接案例证明了其在工程问题中的适用性。1. 引言高重要性工程结构通常在复杂的载荷条件及极端环境中运行，因此其可靠性评估不仅需要估算故障概率，还需判断安全裕度是否充足[1,2]。由于几何尺寸、材料参数、载荷条件、接触状态以及模型误差等因素的影响，结构中往往同时存在随机不确定性与认知不确定性，这使得结构的响应与失效边界呈现出区间特性[[3], [4], [5]]。因此，在混合不确定性条件下获得可靠的故障概率界限，同时降低高精度模型计算的成本，是此类结构可靠性评估中的关键问题。安全裕度与不确定性量化（QMU）已成为高重要性系统安全裕度评估、认证及风险决策的重要方法[[6], [7], [8], [9]]。然而在混合不确定性条件下，系统响应通常不再对应唯一的概率分布，而是呈现区间或分布包络特征。因此，将QMU与混合不确定性传播相结合，是高重要性工程结构可靠性评估的重要方向。针对混合不确定性下的可靠性分析，现有研究主要采用证据理论[[10], [11], [12]]、区间模型[[13], [14], [15]]、概率箱模型[[16], [17], [18]]以及模糊集[[19], [20], [21]]来描述认知不确定性。其中，Dempster-Shafer（DS）证据理论不依赖精确的概率分布，能够通过焦点元素区间与基本概率赋值来描述有限信息条件下的不确定性，因此在工程可靠性分析中被广泛应用。杜等人[22]提出了统一不确定性分析（UUA）框架，将随机不确定性与认知不确定性纳入统一的可靠性分析流程中。张等人[23]将该框架应用于工业机器人定位精度的可靠性分析。江等人[24]提出了最可能焦点元素的概念，并通过均匀性方法将证据变量转化为随机变量。肖等人[25]基于等面积原理对证据变量进行了概率等价变换。段等人[26]进一步提出了全局主动学习克里金可靠性分析方法，该方法通过单一的全局替代模型同时获取所有联合焦点元素的响应，进而估算信念值（Bel）与合理性值（Pl）。不过现有研究主要集中在估算Bel与Pl可靠性界限上，而对于高重要性结构的QMU评估尚未深入探讨。在此基础之上，混合不确定性条件下的QMU评估更侧重于故障概率界限、响应分布包络以及安全裕度判定。Sentz等人[27]将不完全信息条件下的不确定性表示方法拓展为概率界限分析，为QMU中的认知不确定性传播提供了依据。Shah等人[28]将QMU与DS证据理论相结合，并采用非侵入式多项式混沌方法降低计算成本。谢等人[29]通过将证据理论与克里金替代模型结合，建立了混合不确定性条件下的QMU传播流程。尽管替代模型能有效提升计算效率，但当需同时考虑多个认知变量时，现有方法仍多依赖于单个焦点元素的传播或局部替代建模，重复建模与传播的成本依然较高。此外，大多数现有研究仅关注单输出响应或单一故障模式，而实际工程结构系统往往受多种故障模式共同影响，这些故障模式之间还可能存在竞争与切换现象[30,31]。因此，有必要建立一种针对多输出响应、多重边界及多种故障模式的系统级QMU评估方法，并结合敏感性分析识别关键参数，从而为参数控制、结构优化及不确定性降低提供依据。主动学习克里金法为减少高精度性能函数评估的次数提供了一种有效途径。主动学习克里金蒙特卡洛模拟（AK-MCS）[32,33]将克里金替代模型与蒙特卡洛候选样本库相结合，利用预测均值与方差引导新样本向极限状态边界靠拢，从而以更少的真实函数评估次数估算故障概率。此后，通过引入重要性采样、子集模拟以及概率密度函数信息，又相继发展出了基于重要性采样的主动学习克里金方法（AK-IS）[[34], [35], [36]）、基于子集模拟的主动学习克里金方法（AK-SS）[37,38]以及基于概率密度函数信息的主动学习克里金方法（AK-PDF）[39,40]，用以处理低故障概率问题或复杂故障域问题。但这些方法主要聚焦于随机不确定性下的故障概率点估计，其学习准则通常侧重于边界分类的不确定性，难以直接适配证据理论中的焦点元素条件传播与BPA加权融合机制。特别是在混合不确定性条件下，Bel与Pl的估算不仅取决于极限状态边界附近的分类情况，还取决于每个联合焦点元素内响应极值的可靠近似值。因此，基于AK-MCS的主动学习理念，有必要构建一种同时考虑分类不确定性、替代模型预测不确定性以及焦点元素证据权重的学习准则，以便在混合不确定性条件下以较低的计算成本满足QMU评估的需求。为解决上述问题，本文提出了一种用于混合不确定性条件下QMU评估的主动学习克里金方法。首先构建由随机变量与认知变量组成的混合输入空间，并生成初始训练样本。随后建立全局替代模型，通过加权主动学习准则依次选取样本以更新模型。模型收敛后，在每个联合焦点元素内进行条件传播，从而得到全局故障概率界限，再构建响应分布曲线以计算QMU指标。最后通过敏感性分析识别影响评估结果的主要不确定性来源。本文的主要贡献如下：（1）仅需一个全局替代模型，所有联合焦点元素的条件传播与故障概率界限估算均在替代模型层面完成，从而减少了为每个焦点元素单独建模以及高精度性能函数评估所带来的计算成本。（2）提出了基于BPA加权的主动学习函数，使主动学习过程直接指向Bel与Pl估算的精度要求，进而提升了混合不确定性传播的效率。（3）所提方法可直接应用于涉及多输出响应、多重边界及多种故障模式的系统级QMU评估，能够描述不同故障模式之间的竞争关系，并通过敏感性分析识别关键参数，为设计改进提供依据。本文的其余部分结构如下：第2节介绍混合不确定性条件下QMU评估的基本定义；第3节阐述所提出的方法论及其实施流程；第4节通过两个基准示例以及一个楔环连接结构的工程案例对该方法进行验证；第5节对全文进行总结，并探讨未来的研究方向。2. 混合不确定性条件下的QMU问题定义2.1 安全裕度与不确定性指标QMU评估是一种基于传统安全裕度概念的复杂系统风险评估决策支持方法，它充分考虑了系统运行与服务过程中的各种不确定性，以性能安全裕度及其相关不确定性作为核心量[1,7]。图1展示了QMU的概念框架与关键指标。下载：下载高分辨率图像（133KB）下载：下载完整尺寸图像图1 QMU指标示意图。在QMU框架中，系统可靠性通常用置信因子（CF）[8]来表征，其定义为(1)CF=M/U其中M表示性能响应与性能阈值之间的安全裕度，U则表示评估过程中量化的不确定性。如果系统评估中的不确定性能够被充分量化，且CF>1，则认为系统是可靠的；否则则存在潜在的故障风险。对于复杂系统，QMU分析过程通常分为三个阶段，如图2所示。首先，建立可靠性观测清单，以识别潜在的故障模式、故障机理以及关键特征参数。其次，通过定义性能指标、性能阈值及其允许范围来构建性能通道，然后从理论模型、仿真模型或实验数据中获取系统响应。最后，对材料、几何形状、载荷、模型以及测试等方面的不确定性进行表征与传播，进而计算出M、U和CF的值，为系统可靠性评估与工程决策提供定量支持。下载：下载高分辨率图像（556KB）下载：下载完整尺寸图像图2 QMU评估的分析流程图。2.2 混合不确定性表示与故障事件定义在混合不确定性评估中，随机变量向量记为X=(X1,X2,…,XnX)，其中nX表示随机变量的数量，X的一个具体取值记为x=(x1,x2,?,xnX)。认知不确定性通常由信息不足、实验数据匮乏或专家判断不完整引起，很难用唯一且精确的概率密度函数来描述。因此，本研究采用DS证据理论[12]来表征认知不确定性。设Θ为判别框架，任意命题A则是Θ的子集。基本概率赋值被定义为一种映射m:2Θ→[0,1]，且满足m(?)=0；m(A)≥0，?A∈2Θ；∑A∈2Θm(A)=1。其中m(A)表示赋予命题A的BPA值，而那些满足m(A)>0的子集被称为焦点元素。对于任意命题B?Θ，其信念函数Bel(B)与合理性函数Pl(B)分别定义为(2)Bel(B)=∑A?Bm(A)(3)Pl(B)=∑A∩B≠?m(A)因此，命题B的可信度可以用区间[Bel(B),Pl(B)]来表示。信念度量、合理性度量与不确定性区间之间的关系如图3所示[41]。下载：下载高分辨率图像（44KB）下载：下载完整尺寸图像图3 Bel与Pl函数的示意图。根据证据理论，认知变量向量记为Y=(Y1,Y2,?,YnY)，其中nY表示认知变量的数量。第i个认知变量Yi由其证据空间(Уi,Ωi,mi)来表示，其中Уi是Yi的可能取值集合，Ωi是在Уi上定义的焦点元素集合，mi则是对应的BPA值。如果Yi包含ni个焦点元素，那么Ωi={Usi1,Ui2,?,Uini}，其中Uis表示Yi的第s个焦点元素，s=1,2,?,ni，mi(Uis)则是对应的BPA值。如果各认知变量相互独立，那么可以通过单个变量焦点元素的笛卡尔积来构建联合焦点元素，所有联合焦点元素再重新编号为Uj，j=1,2,…,Nf，其中Nf=∏i=1nYni。因此，第j个联合焦点元素可以表示为(4)Uj=U1κ1(j)×U2κ2(j)×?×UnYκnY(j)，其中κi(j)∈{1,2,?,ni}表示在第j个联合焦点元素中选自认知变量Yi的单个变量焦点元素的索引。相应地，第j个联合焦点元素Uj的联合BPA值则是通过将各个对应单个变量焦点元素的BPA值相乘得到的，即(5)mY(Uj)=∏i=1nYmi(Uiκi(j))表1给出了两个认知变量Y1与Y2构建联合焦点元素及联合BPA值的示例。表1 联合焦点元素与联合BPA值的构建示例。Y2Y1U11= [3.0,3.4]U12= [3.4,3.8]U13= [3.8,4.2]空单元格m1(U11)=0.3m1(U12)=0.3m1(U13)=0.4U21=[0.1,0.4]m2(U21)=0.2mY(U1)=0.06mY(U2)=0.06mY(U3)=0.08U22=[0.4,0.7]m2(U22)=0.5mY(U4)=0.15mY(U5)=0.15mY(U6)=0.20U23=[0.7,1.0]m2(U23)=0.3mY(U7)=0.09mY(U8)=0.09mY(U9)=0.12对于涉及多种故障模式的复杂系统，工程实践中通常采用串联系统的故障判定标准，即一旦出现任何一种故障模式，整个系统就被视为发生故障。在传统的QMU处理方法中，系统的CF通常简化为所有单个故障模式中最小的CF值[30,31]，即(6)CFsys=min(CF1,CF2,?,CFnm)=min(M1U1,M2U2,?,MnmUnm,)其中nm为故障模式的数量，Mp与Up分别表示第p种故障模式的安全裕度与不确定性。这种方法实现简单，但实际上只是在结果层面进行标量聚合，因此难以反映不同故障模式之间的竞争关系，以及相同输入状态下主导故障模式的切换过程。为了在响应传播层面统一表示多种故障模式，将故障模式p的原始性能响应定义为Rp(X,Y)，相应的上界性能阈值记为ξp。极限状态函数定义为(7)Gp(X,Y)=ξp?Rp(X,Y)，其中Gp(X,Y)>0表示安全状态，Gp(X,Y)<0则表示故障状态。对于每一种独立的失效模式，都会建立一种克里金模型，用以获取极限状态函数的预测值，记为G^p(X,Y)。为了将不同的失效模式转化为由失效驱动的统一响应，首先从替代预测的极限状态函数中得到失效模式p的响应估计值，如下式所示：R^p(X,Y)=ξp?G^p(X,Y)，其中ξp是对应于失效模式p的性能边界的一个代表值。该值是根据QMU分析中性能阈值的最佳估计YT确定的，其作用在于消除不同失效模式之间的维度与尺度差异。失效模式p的标准化响应进一步定义为：Г^p(X,Y)=R^p(X,Y)ξp=1?G^p(X,Y)ξp。系统标准化响应则是所有独立失效模式中最大标准化响应的值，即：Г^sys(X,Y)=maxp=1,2,?,nmГ^p(X,Y)。因此，系统失效事件可表示为Fsys={Г^sys(X,Y)>1}，相应的系统失效阈值则为ГT=1。在此基础上，对系统标准化响应Г^sys(X,Y)进行混合不确定性传播，从而得到系统裕度Msys、系统不确定性Usys以及系统置信因子，公式为：CFsys=Msys/Usys。

3. 提出的方法论
3.1 在焦点元素约束下的Bel与Pl传播
基于2.2节所建立的混合不确定性表示方法，随机变量X由概率分布描述，而认知变量Y则由联合焦点元素Uj及对应的联合BPA mY(Uj)来描述。对于给定的阈值ξ，阈值事件定义为：F={G(X,Y)< />< />< />< />< />< />< /><ξ}=∑j=1nfmy(uj)ppl,j。这样，混合不确定性下的失效概率便可以区间形式表示为pf∈[bel(f),pl(f)]。通过改变阈值ξ，可构建系统响应的累积信念函数(cbf)与累积合理性函数(cpf)，以此描述混合不确定性下的响应分布范围。不过，若直接在每个联合焦点元素内进行极值搜索并开展嵌套随机传播，就需要在每个焦点元素中反复调用高精度模型，这将带来较高的计算成本。因此，随后会在统一的输入空间中构建一个基于焦点元素条件的全局主动学习克里金替代模型，从而在替代模型层面近似进行极值搜索与随机传播，进而高效地估算出bel(f)和pl(f)。 3.2 用于焦点元素条件传播的全局主动学习克里金替代模型 设组合输入向量为z=(x,y)。全局替代模型是在统一的输入空间ΩX×ΩYbox中构建的，其中ΩX是随机变量的建模覆盖区域，而ΩYbox则由认知变量所有焦点元素的最小值与最大值端点确定，可视为覆盖所有联合焦点元素的最小边界框。在初始的替代模型构建阶段，采用极大极小拉丁超立方抽样(LHS) [42]在标准化空间[0,1]nx+ny中生成n0个初始实验设计(doe)样本，这些样本随后会被映射到物理输入空间ωz中。初始训练集记为d0=>< /><ξ]，其中i[·]为指示函数。若括号内的条件成立，则i[·]=1，否则i[·]=0。最后，利用联合bpa对所有焦点元素上的条件失效概率进行加权聚合，即可得到全局的失效概率区间估算值，即式(31)bel^(f)=∑j=1nfmy(uj)p^bel,j，式(32)pl^(f)=∑j=1nfmy(uj)p^pl,j。 3.3 qmu指标的计算 为计算cf，需要从混合不确定性下系统响应与性能阈值的分布信息中提取裕度与不确定性。3.2节中得到的bel^(f)与pl^(f)估算值，对应的是固定阈值ξ时的失效概率区间估算值。为了进一步开展qmu评估，需将固定阈值ξ替换为响应扫描阈值c，此时相应事件定义为式(33)：ec=>
3.4 基于变异系数的收敛准则
由于最终的qmu评估依赖于通过对所有联合焦点元素进行bpa加权聚合得到的全局失效概率区间，而非单个焦点元素对应的局部失效概率，因此需要对全局估算值bel^(f)与pl^(f)实施统计精度控制。为评估全局失效概率区间估算值的统计稳定性，定义了全局采样变异系数，即式(44)covbel^=1?bel^(f)nrbel^(f)，式(45)covpl^=1?pl^(f)nrpl^(f)。此处采用如下的全球统计精度准则：式(46)max(covbel^,covpl^)≤εcov，其中εcov为预设的采样精度阈值。在本研究中，收敛阈值取值为εcov=0.05。

3.5 sobol敏感性分析
为识别影响系统响应与qmu评估结果的主要不确定性来源，基于已收敛的全局克里金替代模型开展了sobol全局敏感性分析。sobol方法基于方差分解原理，能够定量评估各输入变量对输出方差的独立贡献与总贡献，常被用于非线性模型的全局敏感性分析[44,45]。统一的输入向量记为z=(x,y)，最终全局克里金替代模型的输出值则作为sobol分析的响应值。设系统响应为h=g^(z)。对于随机变量，需根据其概率分布生成样本；而对于认知变量，则先根据联合bpa对联合焦点元素进行抽样，再在对应的焦点元素区间内抽取具体数值。随后采用saltelli抽样方法来估算一阶sobol指数sh与总效应sth。sh用于衡量变量zh对系统响应方差的独立贡献，而sth则表示变量zh的总贡献，包括其与其它变量的交互作用。通过比较sh与sth，便可识别出主要的不确定性来源以及变量间的交互效应。

3.6 所提算法的实现
综上所述，本文提出了一种基于混合不确定性下主动学习克里金的qmu评估方法，其实现流程如下。图5展示了所提计算框架的流程图。下载：下载高分辨率图像（429kb）下载：下载完整尺寸图像图5. 所提方法的流程图。步骤1：已知认知变量y的焦点元素区间与bpa，构建联合焦点元素uj，计算它们的联合bpa(my(uj))，并生成随机样本池sx。步骤2：在统一的输入空间ωx×ωybox中，采用极大极小lhs方法生成初始的实验设计样本，随后通过高精度模型计算得到响应值，从而构成初始训练集d0。步骤3：基于训练集d构建全局克里金模型，进而得到预测均值g^(z)与预测标准差σ^g(z)。步骤4：在每个uj内部生成认知候选样本池sy,j，再将该池与随机样本池sx结合起来，得到全局主动学习候选池э。步骤5 上边界裕度与不确定性的计算。在获得系统响应与性能阈值的最佳估计值及不确定性之后，上边界问题的性能裕度、不确定性与置信因子可定义为式(41)m=YT?Yfunction，式(42)U=UT+Ufunction，式(43)CF=MU=YT?YfunctionUT+Ufunction。如果性能阈值本身存在随机不确定性、认知不确定性，或是二者结合的混合不确定性，可采用与系统响应相同的分位数提取策略。此时，性能阈值的最佳估计值YT与不确定性UT可从性能阈值的CBF/CPF分布曲线范围内确定。对于下性能边界问题，分位数选取的方向则相反，即最佳估计值与不确定性应从系统响应与性能阈值均不利的方向中提取。相关的QMU分位数提取方法可参考文献[28,29]。 3.4 基于变异系数的收敛准则 由于最终的qmu评估依赖于通过对所有联合焦点元素进行bpa加权聚合得到的全局失效概率区间，而非单个焦点元素对应的局部失效概率，因此需要对全局估算值bel^(f)与pl^(f)实施统计精度控制。为评估全局失效概率区间估算值的统计稳定性，定义了全局采样变异系数，即式(44)covbel^=1?Bel^(F)NRBel^(F)，式(45)COVPl^=1?Pl^(F)NRPl^(F)。此处采用如下的全球统计精度准则：式(46)max(COVBel^,COVPl^)≤εCOV，其中εCOV为预设的采样精度阈值。在本研究中，收敛阈值取值为εCOV=0.05。 3.5 sobol敏感性分析 为识别影响系统响应与qmu评估结果的主要不确定性来源，基于已收敛的全局克里金替代模型开展了sobol全局敏感性分析。sobol方法基于方差分解原理，能够定量评估各输入变量对输出方差的独立贡献与总贡献，常被用于非线性模型的全局敏感性分析[44,45]。统一的输入向量记为z=(x,y)，最终全局克里金替代模型的输出值则作为Sobol分析的响应值。设系统响应为H=G^(z)。对于随机变量，需根据其概率分布生成样本；而对于认知变量，则先根据联合BPA对联合焦点元素进行抽样，再在对应的焦点元素区间内抽取具体数值。随后采用Saltelli抽样方法来估算一阶Sobol指数Sh与总效应STh。Sh用于衡量变量zh对系统响应方差的独立贡献，而Sth则表示变量zh的总贡献，包括其与其它变量的交互作用。通过比较Sh与Sth，便可识别出主要的不确定性来源以及变量间的交互效应。 3.6 所提算法的实现 综上所述，本文提出了一种基于混合不确定性下主动学习克里金的qmu评估方法，其实现流程如下。图5展示了所提计算框架的流程图。下载：下载高分辨率图像（429kb）下载：下载完整尺寸图像图5.>
3.4 基于变异系数的收敛准则
由于最终的qmu评估依赖于通过对所有联合焦点元素进行bpa加权聚合得到的全局失效概率区间，而非单个焦点元素对应的局部失效概率，因此需要对全局估算值bel^(f)与pl^(f)实施统计精度控制。为评估全局失效概率区间估算值的统计稳定性，定义了全局采样变异系数，即式(44)covbel^=1?bel^(f)nrbel^(f)，式(45)covpl^=1?pl^(f)nrpl^(f)。此处采用如下的全球统计精度准则：式(46)max(covbel^,covpl^)≤εcov，其中εcov为预设的采样精度阈值。在本研究中，收敛阈值取值为εcov=0.05。

3.5 sobol敏感性分析
为识别影响系统响应与qmu评估结果的主要不确定性来源，基于已收敛的全局克里金替代模型开展了sobol全局敏感性分析。sobol方法基于方差分解原理，能够定量评估各输入变量对输出方差的独立贡献与总贡献，常被用于非线性模型的全局敏感性分析[44,45]。统一的输入向量记为z=(x,y)，最终全局克里金替代模型的输出值则作为sobol分析的响应值。设系统响应为h=g^(z)。对于随机变量，需根据其概率分布生成样本；而对于认知变量，则先根据联合bpa对联合焦点元素进行抽样，再在对应的焦点元素区间内抽取具体数值。随后采用saltelli抽样方法来估算一阶sobol指数sh与总效应sth。sh用于衡量变量zh对系统响应方差的独立贡献，而sth则表示变量zh的总贡献，包括其与其它变量的交互作用。通过比较sh与sth，便可识别出主要的不确定性来源以及变量间的交互效应。

3.6 所提算法的实现
综上所述，本文提出了一种基于混合不确定性下主动学习克里金的qmu评估方法，其实现流程如下。图5展示了所提计算框架的流程图。下载：下载高分辨率图像（429kb）下载：下载完整尺寸图像图5. 所提方法的流程图。步骤1：已知认知变量y的焦点元素区间与bpa，构建联合焦点元素uj，计算它们的联合bpa(my(uj))，并生成随机样本池sx。步骤2：在统一的输入空间ωx×ωybox中，采用极大极小lhs方法生成初始的实验设计样本，随后通过高精度模型计算得到响应值，从而构成初始训练集d0。步骤3：基于训练集d构建全局克里金模型，进而得到预测均值g^(z)与预测标准差σ^g(z)。步骤4：在每个uj内部生成认知候选样本池sy,j，再将该池与随机样本池sx结合起来，得到全局主动学习候选池э。步骤5>对于每个候选点zmjl∈Э，利用公式(24)计算BPA加权的主动学习函数LAL(zmjl)，再根据公式(25)选择新的样本z*。步骤6：若满足公式(26)，则进入下一步；否则，使用高精度模型评估z*的真实响应，更新训练集，然后返回步骤3。步骤7：利用克里金模型，根据公式(27)和(28)估算每个Uj内的最大响应G^j,max(xm)和最小响应G^j,min(xm)，再根据公式(29)和(30)计算失效概率的下限估计值P^Bel,j和上限估计值P^Pl,j。步骤8：判断全局采样系数变异度是否满足公式(46)。若不满足，则扩大随机样本池SX，返回步骤7。步骤9：根据公式(31)和(32)，用mY(Uj)对所有联合焦点元素的结果进行加权聚合，从而得到全局失效概率区间估计值Bel^(F)和Pl^(F)。步骤10：扫描响应阈值c，构建不确定性分布曲线，并根据公式(38)–(43)计算M、U和CF。步骤11：进行索伯尔敏感性分析，识别不确定性的主要来源，为设计改进提供依据。

4. 案例研究
为验证所提方法的有效性、计算效率及工程适用性，选取了两个基准示例和一个楔环连接结构的工程案例。首先通过多输出数值示例，检验该方法在涉及多个边界和多种失效模式的系统级QMU评估中的适用性；其次将曲柄滑块机构案例与已有研究成果对比，评估其在混合不确定性传播、失效概率区间估算以及QMU指标计算方面的优势；最后将该方法应用于实际的高后果工程结构——楔环连接结构，并通过与现有方法的对比及索伯尔敏感性分析，进一步验证其在高后果工程结构安全裕度评估中的适用性与工程价值。

4.1 示例1：数学示例
本小节采用参考文献[28]中的数学示例。该示例包含两个输出响应F1和F2，其中F1由四维罗森布鲁克函数表示，F2则由麦考密克函数加上耦合项F1/100构成：
$$F1(X,Y)=(1?X1)2+100(Y1?X12)2+(1?Y1)2+100(X2?Y12)2+(1?X2)2+100(Y2?X22)2$$
$$F2(X,Y)=sin(X1+Y1)+(X1?Y1)2?1.5X1+2.5Y1+1+F1/100$$
其中X1和X2为随机变量，Y1和Y2为证据变量。此外，性能边界函数还涉及随机变量X3和X4以及证据变量Y3。为保证公平比较，随机变量的概率分布、认知变量的焦点元素区间以及对应的BPA值均取自参考文献[28]中的基准设定。这种设定可理解为在信息受限条件下基于证据理论的离散化处理，其中焦点元素区间描述了认知变量可能的取值范围，BPA值则表示分配给各区间的相对证据支持程度。相关的不确定性信息列于表2和表3中。

表2 示例1的随机变量
| 变量 | 均值 | 标准差 | 分布 |
|------|------|--------|------|
| X1 | 0.15 | 0.38 | 均匀分布[-0.5,0.8] |
| X2 | 0.25 | 0.03 | 正态分布 |
| X3 | 3 | 0.00 | 正态分布 |
| X4 | -100 | 0.00 | 正态分布 |

表3 示例1的证据变量与BPA值
| 变量 | 区间 | BPA值 |
|------|------|-------|
| Y1 | [-0.50,-0.10] | 0.167 |
| | [0.10,0.40] | 0.166 |
| | [0.00,0.50] | 0.111 |
| | [-0.40,0.20] | 0.111 |
| | [-0.10,0.10] | 0.111 |
| | [0.25,0.35] | 0.117 |
| | [-0.45,-0.29] | 0.217 |
| Y2 | [0.20,1.00] | 0.150 |
| | [-1.00,0.40] | 0.350 |
| | [-0.20,0.30] | 0.167 |
| | [-0.50,-0.15] | 0.167 |
| | [0.15,0.90] | 0.166 |
| Y3 | [6.00,6.50] | 0.500 |
| | [6.30,6.75] | 0.300 |
| | [5.90,6.20] | 0.200 |

注：Y1和Y2的BPA值是根据参考文献[28]中的平均规则，融合多种证据信息得到的；而Y3则是来自单一来源的证据变量。该示例同时包含下限和上限性能边界。F1的下限边界定义为：
$$FL,1=X4$$
F1的上限边界定义为：
$$FU,1=(X3+2Y1?7)2+(2X3+Y1?5)2$$
F2的上限边界定义为：
$$FU,2=(X1?1)?+(Y3?1)?$$
因此，三种失效模式对应的极限状态函数可定义为：
$$G1(X,Y)=F1(X,Y)?FL,1$$
$$G2(X,Y)=FU,1?F1(X,Y)$$
$$G3(X,Y)=FU,2?F2(X,Y)$$
运用所提方法，可得到响应、性能边界以及系统归一化响应的CBF/CPF分布，如图6所示。

为降低单次主动学习过程中随机采样的影响，分析时采用了10次独立运行的平均值。在置信水平γ=0.95时，QMU评估结果列于表4中。此处Ncall表示高精度性能函数评估的次数。对于系统级结果，Ncall指的是系统评估过程中实际需要的唯一调用次数，而非各失效模式调用次数的简单相加。

表4 示例1的QMU评估结果汇总
| 评估对象 | 方法 | UCF | Ncall（1次运行） | Ncall（10次运行） |
|----------|------|-----|----------------|------------------|
| FL,1 | DSTE-NIPC [28] | 109.859 | -0.832 | -140 |
| | 所提方法 | 115.416 | 115.288 | 0.963 | 0.963 | 998 | 8.7 |
| FU,1 | DSTE-NIPC [28] | 1292.169 | -1.644 | -152 |
| | 所提方法 | 1441.165 | 1435.142 | 2.360 | 2.376 | 423 | 7.4 |
| FU,2 | DSTE-NIPC [28] | 339.712 | -1.684 | -170 |
| | 所提方法 | 342.967 | 342.957 | 1.818 | 1.818 | 413 | 7.8 |
| 系统整体 | DSTE-NIPC [28] | min(CFp) | -0.832 | -322 |
| | 所提方法 | min(CFp) | -0.963 | 0.963 | 142 | 123.9 |
| | 所提方法计算的系统G^sys值 | 0.375 | 0.397 | 2.516 | 2.343 | 142 | 123.9 |

注：参考文献[28]中，FL,1和FU,1的DSTE-NIPC QMU指标来自表8，FU,2的指标来自表9，Ncall值则来自表7，其中152=140+12，170=140+30，322=152+170。系统的CF值为min(0.832,1.644,1.684)=0.832。如表4所示，所提方法与DSTE-NIPC对三种单独的失效模式得出的QMU安全判定结果一致，但U和C的具体数值存在差异。对于FL,1，其CF值从DSTE-NIPC得到的0.832上升至所提方法得到的0.963。虽然该值仍小于1，但更接近安全边界，说明FL,1仍是系统的主要风险控制模式。对于FU,1和FU,2，其CF值大于1，表明安全裕度相对充足。这些差异主要源于两种方法所采用的混合不确定性传播策略、替代模型更新流程以及焦点元素极值估算策略的不同。在系统级评估中，若采用传统的min(CFp)准则，所提方法得到的系统CF值为0.963，高于DSTE-NIPC的0.832，不过主导失效模式并未发生变化。此外，还基于系统归一化响应Γsys对多种失效模式进行了响应级竞争性聚合，得到Usys=0.397，CFsys=2.343。这一结果并非通过取各失效模式的最低CF值得到，而是通过在统一的归一化响应空间中传递各失效模式之间的竞争关系得以得出。因此，Γsys可作为min(CFp)准则的补充，揭示相同输入状态下不同失效模式之间的竞争关系，以及它们对系统安全裕度和不确定性的综合影响。在计算效率方面，所提方法大幅减少了高精度性能函数评估的次数。与DSTE-NIPC相比，三种单独失效模式的评估次数分别减少了约36.6%、75.4%和77.8%，系统级评估的次数则减少了约61.5%。这表明所提方法在保持一致QMU判定结果的同时，有效降低了计算成本，在单个失效模式评估和系统级评估中都具有较高的计算效率。

4.2 示例2：曲柄滑块机构
第二个示例是工程机械中常用的曲柄滑块机构[29]，其结构示意图如图7所示。

图7 曲柄滑块机构的结构示意图

曲柄长度a、联轴器长度b、外部作用力p、联轴器材料的杨氏模量E以及联轴器的屈服强度S均被建模为随机变量，具体参数列于表5中。

表5 示例2的随机变量
| 变量 | 参数符号 | 均值 | 标准差 |
|------|----------|------|--------|
| X1 | a（毫米） | 1000 | 0.01 |
| X2 | b（毫米） | 4000 | 0.01 |
| X3 | p（千牛） | 280 | 28 |
| X4 | E（吉帕） | 200 | 10 |
| X5 | S（兆帕） | 290 | 29 |
| X6 | d1（毫米） | 603 | |
| X7 | d2（毫米） | 252.5 | |

偏移位移ζ和摩擦系数μ被视为认知变量，它们的焦点元素区间及BPA值列于表6中。为保证比较的一致性，随机变量的概率分布、认知变量的焦点元素区间以及对应的BPA值均采用参考文献[29]中的设定。

表6 示例2的证据变量与BPA值
| 变量 | 参数符号 | 区间 | BPA值 |
|------|----------|------|-------|
| Y1 | ζ | [100,120] | 0.2 |
| | [120,140] | 0.4 |
| | [140,150] | 0.4 |
| Y2 | μ | [0.15,0.18] | 0.3 |
| | [0.18,0.23] | 0.3 |
| | [0.23,0.25] | 0.4 |

曲柄滑块机构的极限状态函数定义为屈服强度与联轴器最大组合应力之差：
$$g1(X,Y)=S?4p(b?a)π((b?a)2?θ2?μθ)(d22?d12)$$

在QMU评估框架中，联轴器的最大组合应力g2(X,Y)被视为系统响应Y函数，而屈服强度S则作为性能阈值YT，即：
$$Y函数=g2(X,Y)=4P(b?a)π((b?a)2?e2?μe)(d22?d12)$$
$$YT=S$$

运用所提方法，可得到混合不确定性条件下系统响应和性能阈值的不确定性分布，如图8所示。

图8 系统响应与性能阈值的不确定性分布

为降低单次主动学习过程中随机采样的影响，分析时采用了10次独立运行的平均值。不同置信水平γ下的QMU评估结果列于表7中。

表7 示例2的QMU评估结果汇总
| γ值 | UCF | CF |
|----|-----|----|
| 1次运行 | 117.06 | 128.46 |
| 10次运行 | 115.288 | 1.1956 |
| 1次运行 | 0.963 | 0.963 |
| 10次运行 | 138.00 | 0.9674 |
| 0.95 | 40.909 | 8.6010 |
| 1.00 | 17.50 | 1.1360 |
| 0.80 | 77.71 | 86.38 |
| 1.93 | 17.779 | 293.50 |
| 1.42 | 80.00 | |

注：参考值U和CF取自参考文献[29]中的表6。如表7所示，所提方法在不同置信水平下都能得到稳定的QMU评估结果。随着γ值减小，U值下降而CF值上升，这一趋势与QMU评估规律相符。以γ=0.95为例，参考文献[29]中报告的U值为138.0，CF值为0.9674，而所提方法得到的平均值为128.46和1.1956。这说明所提方法能够降低不确定性，提高置信度，同时使CF值从小于1上升至大于1。这些结果表明，BPA加权的主动学习策略能更有效地引导样本向临界边界区域靠近，进而提升QMU安全评估能力。此外，还计算了失效概率区间以及高精度性能函数评估的次数，结果列于表8中。

表8 示例2的可靠性区间与计算成本
| 方法 | Bel（×10??） | Pl（×10??） | (Pl - Bel)（×10??） | Ncall |
|------|-------------|-------------|-------------------|--------|
| 参考文献[29] | 4.44 | 1.44 | 9.969 | 900 |
| 所提方法 | 3.88 | 1.04 | 6.525 | 57.6 |

注：参考区间和Ncall值取自参考文献[29]第6.1节的内容。如表8所示，通过10次重复运行得到的失效概率区间与参考文献[29]报告的区间量级一致，区间宽度为6.52×10??，小于参考文献中给出的9.96×10??，说明不确定性传播的区间更窄。同时，所提方法所需的高精度性能函数平均评估次数为57.6，相比参考文献中的900次减少了约93.6%，这充分体现了所提方法在保持合理可靠性区间估计的同时，具备极高的计算效率。为验证焦点元素极值估算的准确性，还进行了额外的数值验证。在相同条件下，分别使用最终的克里金替代模型和解析性能函数计算焦点元素内的极值，进而得到相应的Bel、Pl值以及失效概率区间宽度，结果列于表9中。此外，为考察极值本身的近似精度，表10给出了Gmin和Gmax的均值相对误差以及失效分类误差。失效分类误差是指替代模型得到的极值与解析极值在安全状态或失效状态分类上出现差异的案例比例。

表9 基于焦点元素极值估算Bel和Pl的验证结果
| 项目 | 克里金替代模型 | 解析函数 | 绝对误差 | 相对误差 |
|------|--------------|----------|----------|----------|
| Bel（×10??） | 3.88 | 4.01 | 0.13 | 3.24% |
| Pl（×10??） | 1.04 | 1.10 | 0.06 | 5.45% |
| (Pl - Bel)（×10??） | 6.52 | 6.99 | 0.47 | 6.72% |

表10 使用克里金替代模型估算的焦点元素极值的误差统计结果
| 项目 | 均值相对误差 | 失效分类误差 |
|------|------------|--------------|
| Gmin | 1.26% | 2.17×10?? |
| Gmax | 3.46% | 1.01×10?? |

如表9所示，克里金替代模型得到的Bel、Pl值以及区间宽度与解析函数的结果非常接近，相对误差分别为3.24%、5.45%和6.72%。表10进一步显示，Gmin和Gmax的均值相对误差分别为1.26%和3.46%，而失效分类误差仅分别为2.17×10??和1.01×10??。这些结果表明，克里金替代模型能够有效捕捉用于估算Bel和Pl的焦点元素极值，且极值估算误差对最终失效概率区间的影响较小。为研究BPA赋值过程对QMU结果的影响，还开展了多层扰动敏感性分析。具体而言，分别将每个焦点元素的BPA值在其原始值基础上扰动±5%、±10%和±20%，之后对扰动后的BPA值进行归一化处理，确保每个认知变量的BPA值之和仍为1。同时，为考察认知变量区间边界变化的影响，还分析了两种额外情况：焦点元素区间扩大10%，以及同时对BPA值和区间边界进行扰动的情况。在置信水平γ=0.95时，每种扰动情况重复10次，结果如表11所示。表11. 示例2中BPA赋值与认知区间的敏感性分析。案例Bel(×10?5)Pl(×10?4)(Pl- Bel)(×10?5)UCFBaseline3.881.046.52128.461.1956BPA±5%3.891.046.51128.481.1958BPA±10%3.901.046.50128.491.1959BPA±20%3.921.046.48128.521.1962区间+10%3.440.976.26128.631.1945BPA±10% +区间+10%3.460.976.24128.641.1939如表11所示，在BPA±5%、BPA±10%和BPA±20%的扰动下，Bel、Pl、U和CF与基准情况相比变化较小，表明BPA赋值的变化对QMU结果的影响有限。当焦点元素区间宽度扩大10%时，U略有上升，CF略有下降，这说明区间扩大使得评估结果更为保守。在综合扰动情况下，CF仍保持在1.19左右且大于1。总体而言，在所有扰动情况下CF均大于1，说明QMU的安全判断并未发生反转，所提出的方法对认知输入设置具有较好的鲁棒性。

4.3. 示例3：楔环连接结构的安全裕度评估楔环连接是一种新型的圆柱壳连接方式。与传统的螺旋连接和管螺纹连接相比，它具有诸多优势，比如表面更光滑、径向尺寸更小、附加质量更低、结构更为紧凑，且安装与拆卸更为方便。因此，它被应用于鱼雷、水雷以及一些特殊结构中。不过，这种结构通常包含多个接触界面，在载荷传递过程中会出现显著的接触非线性、局部应力集中以及多参数耦合效应，使得结构响应对几何尺寸、材料特性和接触参数的不确定性极为敏感。目前，关于楔环连接的可靠性、安全裕度及不确定性量化研究还比较有限。因此，有必要在混合不确定性条件下对楔环连接进行可靠性评估和QMU分析，这有助于识别关键的不确定性来源和潜在的失效风险，为高重要性工程结构的安全设计及可靠运行提供理论和工程支持。由于多个接触界面带来的强非线性，楔环结构的极限状态响应高度依赖于高精度的有限元仿真。基于FEM的直接混合不确定性传播与QMU评估会导致较高的计算成本。为构建后续不确定性建模的输入集，首先从几何尺寸、材料参数和接触参数中筛选出潜在的敏感参数，然后使用这些参数的名义值进行基准有限元仿真。整个有限元分析流程如图9所示。下载：下载高分辨率图像（783KB）下载：下载全尺寸图像图9. 拉伸载荷作用下楔环结构的有限元分析过程。在此基础上，分别在四种工况下进行对比响应分析，即基准工况、楔环厚度变化、预加载变化以及泊松比变化。以名义参数为基准，当拉伸载荷逐渐增加到200 kN时，不同工况下最大间隙和最大应变的变动情况如图10所示。下载：下载高分辨率图像（386KB）下载：下载全尺寸图像图10. 拉伸载荷作用下不同参数条件下楔环结构的响应曲线。如图10所示，在预加载条件下，最大间隙和最大应变均显著增加，这说明预加载改变了结构的初始承载状态，会对变形和应变响应产生持续的放大效应。而泊松比变化下的响应曲线与名义条件下的几乎一致，仅在局部受力区域存在轻微偏差，这表明泊松比对结构的宏观响应影响有限，因此在后续建模中不将泊松比视为不确定输入。楔环结构所选的随机变量和认知变量分别列于表12和表13中。随机变量的均值和标准差是根据结构设计参数、制造公差以及基准有限元输入确定的。由于杨氏模量E和摩擦系数μ缺乏足够的统计样本，因此将其视为认知变量。它们的焦点元素区间和BPA是通过参考材料手册、现有测试信息、初步有限元分析结果以及工程经验，借助专家意见确定的。表12. 示例3的随机变量变量参数符号均值标准差X1步高：h(mm)20.30.233X2半径：R(mm)109.90.033X3内环形槽宽度：w1(mm)8.40.033X4外环形槽宽度：w2(mm)8.40.033X5预加载（过盈配合）：P(mm)0.10.033表13. 示例3的认知变量与BPA变量参数符号区间BPAY1杨氏模量：E(GPa)[199,203]0.5[198,204]0.3[196,206]0.2Y2摩擦系数：μ[0.13,0.17]0.5[0.12,0.18]0.3[0.1,0.2]0.2基于这些不确定输入，最大间隙和最大应变被选为两种独立的失效模式。根据初步有限元仿真、物理试验观测以及工程判断，最大间隙的失效阈值设定为0.28毫米，最大应变的失效阈值设定为8.5×10?3毫米/毫米。随后，构建了最大间隙、最大应变以及系统响应的不确定性分布曲线，如图11所示。下载：下载高分辨率图像（433KB）下载：下载全尺寸图像图11. 多输出响应的不确定性分布曲线。为降低单次主动学习运算中随机采样的影响，分析时采用了10次独立运算的平均结果。在γ=0.95时，QMU评估结果列于表14中。系统级的Ncall并非各失效模式评估次数的简单相加，而是指系统级评估过程中实际需要的高精度评估次数。表14. 示例3的QMU评估结果汇总评估对象方法BelPlUCFNcall系统ET-ALK [26]0.01280.02220.12701.0465105.1所提方法0.01220.02210.12661.0520100.5最大间隙ET-ALK [26]4.3375×10?50.02220.04521.050145.6所提方法4.6728×10?50.02200.04521.050041.9最大应变ET-ALK [26]0.01280.01780.00141.0726105.4所提方法0.01220.01780.00141.069284.7注：ET-ALK的结果是通过将参考文献[26]中的方法应用于当前示例得到的。如表14所示，ET-ALK方法和所提方法在楔环结构的QMU评估结果上总体一致，说明这两种方法在混合不确定性传播方面，无论是对于系统整体还是各个单独的失效模式，都具有良好的一致性。在系统层面，所提方法的Bel略低，Pl几乎相同，U略低，而CF略高于ET-ALK，这说明在保持QMU判断一致性的同时，获得了稳定的系统级评估结果。对于最大间隙，两种方法得到的Pl、U和CF值几乎一致，表明所提方法能够准确识别间隙响应的潜在失效风险。对于最大应变，Pl和U值基本一致，而CF存在轻微差异，说明在应变响应评估方面也具有较好的一致性。在计算效率方面，所提方法所需的有限元评估次数低于ET-ALK方法，分别减少了约4.4%、8.1%和19.6%，这表明所提方法在保持QMU评估结论一致性的同时，进一步降低了有限元性能函数评估的次数，展现了良好的计算效率与工程适用性。同样，为了研究BPA赋值过程对系统级QMU评估结果的影响，对系统输出进行了多层扰动敏感性分析。在置信水平γ=0.95时，每种扰动情况重复10次，结果如表15所示。表15. 示例3系统级QMU评估中BPA赋值与认知区间的敏感性分析。案例Bel(×10?2)Pl(×10?2)(Pl- Bel)(×10?2)UCFBaseline1.22142.20950.98810.12661.0520BPA±5%1.22132.20900.98770.12661.0520BPA±10%1.22122.20920.98800.12661.0519BPA±20%1.22092.21180.99090.12671.0518区间+10%1.20302.54871.34570.12931.0240BPA±10% +区间+10%1.20152.54171.34020.12931.0245如表15所示，在BPA±5%、BPA±10%和BPA±20%的扰动下，Bel、Pl、U和CF与基准情况相比变化较小，CF始终保持在1.05左右且大于1，这说明BPA赋值的变化不会改变QMU的安全判断。当区间宽度扩大10%时，(Pl?Bel)和U上升，而CF下降至1.0240，这说明区间扩大使得评估结果更为保守。在综合扰动情况下，CF仍大于1。总体而言，在所有扰动情况下CF均大于1，说明系统级QMU的安全判断并未发生反转，所提出的方法在系统级评估结果方面具有较好的鲁棒性。在完成QMU评估后，基于收敛后的全局克里金替代模型进行了索伯尔敏感性分析，以识别影响CF的关键输入参数。结果如图12所示。下载：下载高分辨率图像（360KB）下载：下载全尺寸图像图12. 索伯尔敏感性指数排名。如图12所示，系统输出及各个单独失效模式输出的索伯尔敏感性排序总体一致，这说明系统级评估并未改变对主要敏感参数的识别结果。在所有变量中，参数P的一阶索伯尔指数和总效应索伯尔指数最高，因此是影响CF的主要因素。据此，后续的设计改进可以优先考虑仅调整参数P，而保持其他参数不变，以此提升CF，避免因多参数耦合调整而导致的不确定性放大和可追溯性丧失。

5. 结论针对在高随机不确定性和认知不确定性共同作用下的高重要性工程结构的可靠性评估及QMU安全裕度评估，本文提出了一种用于混合不确定性条件下QMU分析的主动学习克里金方法。通过两个基准示例和一个楔环连接结构的工程案例，验证了该方法的有效性、计算效率及工程适用性。主要结论如下：（1）所提方法仅需一个全局主动学习克里金替代模型，即可实现联合焦点元素上的条件传播及失效概率区间估算，从而避免了因逐个焦点元素建模及反复传播而产生的冗余计算。（2）在保证QMU评估结果可信的前提下，与参考方法相比，所提方法可将高精度性能函数评估的次数减少约4.4%–93.6%，大幅降低了混合不确定性条件下的可靠性评估计算成本。（3）所提方法能够处理涉及多输出响应、多重边界条件及多种失效模式的系统级QMU评估。通过系统归一化响应，可以进一步分析不同失效模式之间的竞争关系，避免仅依赖单一失效模式或传统评判标准所带来的信息丢失。（4）楔环连接结构的工程案例证明了所提方法在实际高重要性工程结构可靠性评估中的应用价值。结合索伯尔敏感性分析，所提方法能够识别影响系统级评估结果的关键参数，从而为工程设计优化和不确定性控制提供依据，而不仅仅是像传统评估方法那样获得置信因子。总体而言，所提方法能够在保证可信的失效概率区间和QMU结论的同时，有效降低高精度模型评估的成本。它还将可靠性评估结果转化为关键参数识别、参数控制及结构优化的依据。需注意的是，所提方法更适用于随机变量和认知变量数量适中、联合焦点元素数量可控，且结构响应能够通过克里金替代模型有效近似的问题。而对于涉及高维混合不确定性、随机与认知变量数量急剧增加，或极限状态函数极度不光滑的问题，替代模型构建和主动学习的成本可能依然会上升。未来的工作可以通过变量筛选、焦点元素缩减、自适应焦点元素划分、局部替代模型构建以及稀疏高维替代模型构建等方法，进一步提升所提方法在复杂高维工程系统中的适用性和鲁棒性。CRediT作者贡献说明赵凌：写作——审阅与编辑，写作——初稿撰写，可视化，验证，软件，项目管理，形式分析，数据整理。张辉：资源，方法论。张志林：资源，调研。齐阳：数据整理。索斌：指导，方法论，资金获取。文金鹏：指导，概念构思。