基于参数化概率盒模型的系统可靠性分析方法

System Reliability Analysis Method Based on Parametric Probability Box Model

摘要

针对机械系统存在的一类随机和认知不确定性混合的问题，提出一种基于参数化概率盒（probability box， p-box）模型的系统可靠性分析方法. 基于序列迭代解耦方法获得单失效模式的最小可靠度指标；针对随机与认知不确定性混合的多失效模式问题，建立基于参数化p-box不确定性的系统可靠性分析模型；考虑各失效模式之间的相关性，通过线性相关度分析方法计算得到各失效模式间的相关系数矩阵；提出串联系统和并联系统可靠性求解方法，与传统双层蒙特卡罗采样方法相比，本文方法具有较高的系统可靠性计算效率，能够满足实际工程需求.

Abstract

An efficient system reliability analysis method with a parametric probability box （p-box） model is proposed to solve the problem that the mixture of aleatory and epistemic uncertainty exists in mechanical systems. Firstly， the minimum reliability index of the single failure mode is obtained based on the sequential iteration decoupling method. Then， aiming at the multiple failure modes problem under the mixture of aleatory and epistemic uncertainty， the system reliability analysis model is established based on the parametric p-box uncertainty. Considering the correlation among the failure modes， the correlation coefficient matrix between each failure mode is calculated by the linear correlation degree analysis method. Finally， a reliability calculation method for both series and parallel systems is proposed， which has high reliability calculation efficiency and can meet the practical engineering requirements compared with the traditional double-loop Monte Carlo sampling method.

Graphical abstract

关键词

参数化概率盒 / 可靠性 / 最大失效概率 / 多失效模式

实际工程中广泛存在着与环境载荷、物理特性、几何特性和边界条件等有关的各种不确定性^［1］，针对机械结构不确定性的研究，人们多采用概率模型^［2］、模糊模型^［3］、非概率凸集模型^［4］等来描述系统的不确定性，并根据各自的理论建立相应的可靠性分析方法. 然而，随机可靠性和模糊可靠性方法都是以概率为基础的，在样本数据不足或样本无法获得的情况下，则难以构建变量的概率分布函数. 非概率凸集模型虽然只需要不确定性变量的上下界限，并能在一定程度上解决概率方法无法解决的不确定性问题，但会导致过于保守的分析结果. 因此，鉴于概率模型和非概率凸集模型各自的优缺点，近年来发展出了一种新的可靠性分析模型—概率盒（p-box）模型^［5-6］. 相比于概率模型和非概率凸集模型，p-box模型既不需要太多的样本点，也不会浪费相关的不确定性信息，即使在概率信息不充分的条件下也能构建由上下边界函数所构成的累积分布函数. 因此，p-box模型在结构可靠性分析领域具有重要的工程应用价值，并受到国内外研究人员的关注和重视.

目前，基于p-box模型的不确定性量化与可靠性分析研究已取得了一定的成果. Zhang等^［7］提出了一种直接区间蒙特卡罗采样法，并将该方法应用于含p-box不确定变量的结构有限元分析. Xiao等^［8］针对含参数化p-box不确定性的结构可靠性问题，开发了基于顶点分析的采样方法. Chen等^［9］针对参数化p-box不确定均质材料结构-声场耦合系统的可靠性分析问题，提出了基于一阶矩阵分解摄动有限元的区间蒙特卡罗方法. Wei等^［10］提出了扩展蒙特卡罗采样方法进行p-box可靠性分析及全局敏感性分析. 上述采样类可靠性分析方法往往需要大量重复的采样计算以获得可靠的计算结果. 为避免重复耗时的采样分析，提高p-box可靠性分析的计算效率，Crespo等^［11］提出一种基于Bernstein多项式展开和区间算法的p-box可靠性分析方法，用于计算结构失效概率的变化范围. 刘鑫等^［12］针对不确定因素对载人空降气囊座椅防护性能的影响，提出一种基于概率-概率盒混合模型的气囊座椅防护特性可靠性分析方法. 黄洪钟等^［13］基于多变量函数降维技术和通用生成函数法则，提出p-box不确定性条件下机床主轴可靠性建模与分析方法. Wang等^［14］基于线性规划理论，提出一种结构p-box可靠性窄界限求解方法. 丁家满等^［15］基于概率盒模型和笛卡儿运算法则，提出一种新的不确定性信息条件下系统可靠性分析方法. 项涌涌等^［16］针对参数p-box不确定性情况的模型确认问题，提出面向响应准确度的模型确认方法. Liu等^［17］基于最大熵准则，提出了一种p-box不确定性建模及混合可靠性分析方法. Zhang等^［18］基于序列优化可靠性评估策略，提出了一种自适应Kriging参数化p-box可靠性优化设计方法. Xie等^［19］针对非参数化概率盒的不确定性传播问题，提出了一种基于降维积分和最大熵的响应函数概率边界计算方法. Xiao等^［20］提出了一种计算模型概率边界的协同区间拟蒙特卡罗方法. 但是实际工程结构往往复杂多样，通常涉及多种失效模式或者多个元件的失效状态，如航空发动机涡轮叶片的失效，可能同时存在蠕变、疲劳和磨损等失效模式，因此结构可靠性实际上是系统可靠性.

相比单一失效模式的可靠性问题，结构系统的可靠性问题通常十分复杂，对于含参数化p-box不确定性的结构可靠性问题目前已有少量研究^［21-22］，主要是针对简单的线性失效模式. 然而在实际工程问题中，失效模式往往是多样相关且非线性的，直接假设线性失效模式则会导致较大的计算误差，故发展能解决非线性失效模式的含p-box不确定性的系统可靠性分析方法具有重要的工程意义.

本文用参数化p-box变量来表示系统中的随机与认知混合不确定参数，并提出了一种新的基于参数化p-box模型的系统可靠性分析方法，可以获得多个非线性失效模式下系统可靠度的区间. 首先，基于序列迭代解耦方法计算单个失效模式下的结构参数化p-box可靠度；再针对多失效模式问题，建立基于参数化p-box不确定性的系统可靠性分析模型；考虑各失效模式之间的相关性，通过线性相关度方法计算各失效模式间的相关系数；最后分别提出了串、并联系统可靠性计算方法. 3个数值算例表明，本文方法可以实现含参数化p-box不确定性的多个非线性失效模式下系统可靠度的计算，具有较高的计算效率和精度，能够满足工程需求.

1 单失效模式下系统可靠性分析

结构的可靠概率是一个与随机变量数目有关的高维积分，且极限状态方程通常是非线性函数，因此给实际工程的应用与分析带来了较多困难. 为了计算方便，研究者提出了多种计算方法，其中一阶可靠性方法^［23］（the first order reliability method， FORM）因其方便有效性而得到广泛应用. 下面对其主要思想进行简单介绍.

若系统的输入为随机变量，且变量之间相互独立，则单失效模式下结构的失效概率为：

\begin{array}{l}{p}_{\mathrm{f}}=\mathrm{P}\mathrm{r}\left\{g\left(\mathit{X}\right)\le \mathrm{0}\right\}={\int }_{g\left(\mathit{X}\right)\le \mathrm{0}}{f}_{\mathit{X}}\left(\mathit{x}\right)\mathrm{d}x=\\ {\int }_{g\left(\mathit{X}\right)\le \mathrm{0}}{f}_{X_{\mathrm{1}}}(x_{\mathrm{1}})f_{X_{\mathrm{2}}}(x_{\mathrm{2}})\cdots f_{X_n}(x_n)\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}}\cdot \cdot \cdot \mathrm{d}{x}_n\end{array}（1）

式中：\mathrm{P}\mathrm{r}\{\cdot \}表示概率；g(\mathit{X})为功能函数；\mathit{X}=(X_{\mathrm{1}},X_{\mathrm{2}},\cdots ,X_n)表示n维随机变量；f_{X_i}(x_i)是随机变量X_i的概率密度函数，i=\mathrm{1,2},\cdots ,n. FORM方法可以简化式（1），近似求解结构的可靠度. 首先要对随机变量开展全概率变换，将原始空间中的随机向量\mathit{X}映射到标准正态空间中，得到新的标准正态向量\mathit{U}=(U_{\mathrm{1}},U_{\mathrm{2}},\cdots ,U_n{)}^{\mathrm{T}}，其映射关系为：

U_i={\Phi }^{-\mathrm{1}}\left[F_{X_i}\left(x_i\right)\right]（2）

式中：F_{X_i}(\cdot )为随机变量X_i的累积分布函数；{\Phi }^{-\mathrm{1}}为标准正态变量累积分布函数的逆函数.可进一步获得如下功能函数映射表达式：

g\left(\mathit{X}\right)=g\left[T\left(\mathit{X}\right)\right]=G\left(\mathit{U}\right)（3）

式中：T(\mathit{X})为转换函数. 如果随机变量X_i累计分布函数F_{X_i}(x_i)中存在区间值的分布参数向量\mathit{\theta }={\left({\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},\cdots ,{\theta }_m\right)}^{\mathrm{T}}， m为分布参数的数量，则变量X_i的参数化p-box形式可以表示为：

{\mathit{F}}_{X_i}^{\mathrm{P}}=\left\{F_{X_i}\left(x_i,\mathit{\theta }\right):\mathit{\theta }\in \left({\mathit{\theta }}^{\mathrm{L}},{\mathit{\theta }}^{\mathrm{U}}\right)\right\}（4）

式中：{\mathit{F}}_{X_i}^{\mathrm{P}}的上标代表参数化的p-box；F_{X_i}(x_i,\mathit{\theta })为其累计分布函数；上标 L 和 U 分别表示区间的下界和上界，在实际工程中通常采用区间估计方法获得分布参数的区间. 在式（2）全概率变换过程中，不可避免地将参数化p-box的参数\mathit{\theta }引入标准正态空间中的极限状态方程中，此时式（3）中的转换关系为：

g\left(\mathit{X}\right)=g\left[T\left(\mathit{X},\mathit{\theta }\right)\right]=G\left(\mathit{U},\mathit{\theta }\right)（5）

若没有\mathit{\theta }的存在，极限状态曲线映射到标准正态空间中能得到唯一的一条曲线及其唯一对应的可靠度指标\beta. 然而在p-box模型中，由于分布参数区间化，原空间中的功能函数映射到标准正态空间后构成的极限状态曲面不再是一条曲线，而是由两个边界面\underset{\mathit{\theta }}{\mathrm{m}\mathrm{a}\mathrm{x}}G\left(\mathit{U},\mathit{\theta }\right)=\mathrm{0}和\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}G\left(\mathit{U},\mathit{\theta }\right)=\mathrm{0}构成的极限状态带^［24］，如图1所示.

极限状态带的两条边界上最大概率点（most probable point， MPP）对应的β不再是确定值，而是一个变化区间：

\beta =\left[{\beta }^{\mathrm{L}},{\beta }^{\mathrm{U}}\right]（6）

式中：{\beta }^{\mathrm{L}}为可靠性指标最小值；{\beta }^{\mathrm{U}}为可靠性指标最大值.基于FORM方法，可构造如下两个优化问题^［25］，获得上述极限状态带可靠度指标的上、下边界：

\left\{\begin{array}{l}{\beta }^{\mathrm{L}}=\underset{\mathit{U}}{\mathrm{m}\mathrm{i}\mathrm{n}}‖\mathit{U}‖\\ \mathrm{s}.\mathrm{t}.\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}G\left(\mathit{U},\mathit{\theta }\right)=\mathrm{0}\end{array}\right.（7）

\left\{\begin{array}{l}{\beta }^{\mathrm{U}}=\underset{\mathit{U}}{\mathrm{m}\mathrm{i}\mathrm{n}}‖\mathit{U}‖\\ \mathrm{s}.\mathrm{t}.\underset{\mathit{\theta }}{\mathrm{m}\mathrm{a}\mathrm{x}}G\left(\mathit{U},\mathit{\theta }\right)=\mathrm{0}\end{array}\right.（8）

式（7）和式（8）为双层嵌套优化问题，可采用高效的序列迭代解耦方法^{［24，26］}进行求解，每次迭代过程中依次进行概率分析和区间分析，经过多次迭代最终达到稳定解.

下面以求解最小可靠度指标{\beta }^{\mathrm{L}}为例，给出详细的求解过程. 令第i步迭代计算得到{\mathit{U}}^i和{\mathit{\theta }}^i，在下一步迭代计算中，固定区间分布参数向量{\mathit{\theta }}^i，再利用改进的HL-RF迭代法，即iHL-RF^［27］求得{\mathit{U}}^{i+\mathrm{1}}：

\left\{\begin{array}{l}{\beta }^i=\frac{G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)-{\left(\nabla G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)\right)}^{\mathrm{T}}{\mathit{U}}^i}{‖\nabla G\left({\mathit{U}}^k,{\mathit{\theta }}^k\right)‖}\\ {\mathit{U}}^{i+\mathrm{1}}={\mathit{U}}^i+\alpha {\mathit{d}}^i\end{array}\right.（9）

式中：{\mathit{d}}^i为搜索方向.

{\mathit{d}}^i=\frac{\nabla G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right){\mathit{U}}_i^{\mathrm{T}}-G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)}{{‖\nabla G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)‖}^{\mathrm{2}}}\nabla G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)-{\mathit{U}}^i（10）

式中：

\nabla G\left({\mathit{U}}^i,{\mathit{\theta }}^i\right)={\left(\frac{\partial G}{\partial U_{\mathrm{1}}},\frac{\partial G}{\partial U_{\mathrm{2}}},\cdots ,\frac{\partial G}{\partial U_n}\right)}_{{\mathit{U}}^i,{\mathit{\theta }}^i}（11）

\alpha为迭代步长，由最小化价值函数m确定^［27］.

\left\{\begin{array}{l}m\left(\mathit{U},\mathit{\theta }\right)=\frac{\mathrm{1}}{\mathrm{2}}‖\mathit{U}‖+h\left|G\left(\mathit{U},\mathit{\theta }\right)\right|\\ h>\frac{‖\mathit{U}‖}{‖\nabla G\left(\mathit{U},\mathit{\theta }\right)‖}\end{array}\right.（12）

式中：h为常数. 获得{\mathit{U}}^{i+\mathrm{1}}后，再通过内层优化求解得到{\mathit{\theta }}^{i+\mathrm{1}}.

\left\{\begin{array}{c}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}G\left({\mathit{U}}^{i+\mathrm{1}},\mathit{\theta }\right)\\ \mathrm{s}.\mathrm{t}.\mathit{\theta }\in [{\mathit{\theta }}^{\mathrm{L}},{\mathit{\theta }}^{\mathrm{U}}]\end{array}\right.（13）

通过上述步骤多次迭代，假如满足‖{\mathit{U}}^{i+\mathrm{1}}-{\mathit{U}}^i‖/‖{\mathit{U}}^i‖\le {\epsilon }_{\mathrm{1}}且\left|G({\mathit{U}}^{i+\mathrm{1}},{\mathit{\theta }}^{i+\mathrm{1}})\right|\le {\epsilon }_{\mathrm{2}}，其中{\epsilon }_i\left(i=\right.\left.\mathrm{1,2}\right)为给定收敛误差，则迭代停止，获得最小可靠性指标{\beta }^{\mathrm{L}}=‖{\mathit{U}}^{i+\mathrm{1}}‖；否则，i=i+\mathrm{1}进入下一次迭代；根据p_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\Phi (-{\beta }^{\mathrm{L}})求得最大失效概率. 同理，类似上述过程求解式（8），并由p_{\mathrm{f}}^{\mathrm{m}\mathrm{i}\mathrm{n}}=\Phi (-{\beta }^{\mathrm{U}})计算出最小失效概率. 在实际工程中，最大失效概率通常是人们最为关心和重视的指标，因此，本文后续将以最大失效概率来分析含参数化p-box不确定性的系统可靠性.

2 多失效模式下系统可靠性分析

第1节中可靠度计算方法解决了单一失效模式问题，但在实际工程中，结构的失效模式通常不止一个，如航空发动机涡轮叶片的失效，可能同时存在蠕变、疲劳和磨损等失效模式；各失效模式因受相同的输入载荷、几何和物理参数等因素影响而具有一定的相关性^［28］. 实际上，系统可靠性模型一般可分为以下3类^［29］：串联、并联和混联系统模型，如图2所示.

图2中，B_i （i=1，2，…，a₁）和B_j （j=1，2，…，b₁）分别为串联系统和并联系统中的失效模式.为方便后续理论推导，串联、并联和混联系统的失效概率可统一记为：

（14）

式中：a_{\mathrm{1}}和b_{\mathrm{1}}分别表示串联和并联单元数目，当a_{\mathrm{1}}=b_{\mathrm{1}}=\mathrm{1}时，式（14）表示串联系统失效概率；当a_{\mathrm{1}}>\mathrm{1},b_{\mathrm{1}}=\mathrm{1}时，式（14）表示并联系统失效概率；当a_{\mathrm{1}}=\mathrm{1},b_{\mathrm{1}}>\mathrm{1}时，式（14）表示混联系统失效概率. 实际上，串联系统和并联系统是用来建立任何系统的两个基本系统，而在混联系统中，一般可将一个并联系统作为一个失效模式看待，并可以建立与其对应的功能函数，从而可将并联-串联系统可靠性简化为串联系统可靠性计算^［30-31］，这里不再赘述.

计算含参数化p-box不确定性的系统可靠度的主要困难在于需要考虑存在区间分布参数向量的各失效模式之间的相关性.

2.1 多失效模式间的相关性分析

设含有参数化p-box不确定性的系统具有K_a个失效模式，各失效模式对应的功能函数为：

g_i\left(\mathit{X}\right)=g\left(X_{\mathrm{1}},X_{\mathrm{2}},\cdots ,X_n\right),i=\mathrm{1,2},\cdots ,K_a（15）

基于式（5），将式（15）进行全概率变换，得：

Z_i=g_i\left(\mathit{X}\right)=g_i\left[T\left(\mathit{U},\mathit{\theta }\right)\right]=G_i\left(\mathit{U},\mathit{\theta }\right)（16）

为计算各功能函数之间的相关性，可在Z_i的最大概率点({\mathit{U}}^{\mathit{*}},{\mathit{\theta }}^{\mathit{*}})处固定区间分布参数向量{\mathit{\theta }}^{\mathrm{*}}，将 式（16）对随机向量 U 进行一阶Taylor级数展开：

Z_{Li}=‖\nabla G_i({\mathit{U}}^{\mathrm{*}},{\mathit{\theta }}^{\mathrm{*}})‖({\beta }_i^{\mathrm{L}}-{\mathit{\alpha }}_{U,i}^{\mathrm{T}}\mathit{U})（17）

式中：{\beta }_i^{\mathrm{L}}为对应于功能函数Z_{Li}的最小可靠度指标，可通过式（9）~式（13）求得；{\mathit{\alpha }}_{U,i}^{\mathrm{T}}为第i个失效模式的单位梯度向量，表示为：

{\mathit{\alpha }}_{U,i}^{\mathrm{T}}=\left(\frac{\partial G({\mathit{U}}^{\mathrm{*}},{\mathit{\theta }}^{\mathrm{*}})}{\partial U_{\mathrm{1}}}\right.,\frac{\partial G({\mathit{U}}^{\mathrm{*}},{\mathit{\theta }}^{\mathrm{*}})}{\partial U_{\mathrm{2}}},\cdots ,\left.\frac{\partial G({\mathit{U}}^{\mathrm{*}},{\mathit{\theta }}^{\mathrm{*}})}{\partial U_n}\right)（18）

由式（18）可进一步计算第i个和第j个功能函数间的相关系数^［30］.

{\rho }_{ij}=\frac{\mathrm{C}\mathrm{o}\mathrm{v}(Z_{Li},Z_{Lj})}{{\sigma }_{Z_{Li}}{\sigma }_{Z_{Lj}}}={\mathit{\alpha }}_{U,i}^{\mathrm{T}}{\mathit{\alpha }}_{U,j},i\ne j（19）

式中：\mathrm{C}\mathrm{o}\mathrm{v}(Z_{Li},Z_{Lj})表示Z_{Li}和Z_{Lj}的协方差；{\sigma }_{Z_{Li}}和{\sigma }_{Z_{Lj}}分别表示Z_{Li}和Z_{Lj}的标准差.

2.2 串联系统可靠性求解

系统极限状态带边界如图3所示.分析含有参数化p-box不确定性的串联系统（包含a₁个失效模式），由于存在区间分布参数，原空间中的极限状态方程映射到标准正态空间后将生成极限状态曲面“带”，如图3（a）所示.

串联系统曲面带的下边界和上边界可分别表示为：

S_{\mathrm{s}}^{\mathrm{L}}:{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cup }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}{G}_i(\mathit{U},\mathit{\theta })=\mathrm{0},S_{\mathrm{s}}^{\mathrm{U}}:{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cup }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{a}\mathrm{x}}{G}_i(\mathit{U},\mathit{\theta })=\mathrm{0}（20）

因此系统失效概率P_f为一区间，其边界可通过求解如下两个优化获得：

（21）

（22）

为计算串联系统最大失效概率，可将式（22）转换为：

\mathrm{1}-\left\{{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cap }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}{G}_i\left(\mathit{U},\mathit{\theta }\right)>\mathrm{0}\right\}=

\mathrm{1}-\left\{\left.{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cap }}}\left[‖\nabla G_i({\mathit{U}}^{\mathit{*}},{\mathit{\theta }}^{\mathit{*}})‖({\beta }_i^{\mathrm{L}}-{\mathit{\alpha }}_{U,i}^{\mathrm{T}}\mathit{U})>\mathrm{0}\right]\right\}\right.（23）

结合式（16）和式（17），记L_i={\beta }_i^{\mathrm{L}}-{\mathit{\alpha }}_{U,i}^{\mathrm{T}}\mathit{U}，则式（23）可变为：

\begin{array}{l}{P}_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1}-\mathrm{P}\mathrm{r}\left\{\left.{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cap }}}\left[({\beta }_i^{\mathrm{L}}-{\mathit{\alpha }}_{U,i}^{\mathrm{T}}\mathit{U})>\mathrm{0}\right]\right\}\right.=\\ \begin{array}{cc}& \end{array}\mathrm{1}-\mathrm{P}\mathrm{r}\left\{\left.{\displaystyle \underset{i=\mathrm{1}}{\overset{a_{\mathrm{1}}}{\cap }}}\left[L_i>\mathrm{0}\right]\right\}\right.\end{array}（24）

根据串联系统失效概率计算方法，式（24）可变为：

P_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1}-{\Phi }_{a_{\mathrm{1}}}\left({\mathit{\beta }}^{\mathrm{L}},\mathit{\rho }\right)（25）

式中：{\Phi }_{a_{\mathrm{1}}}(\cdot )为a_{\mathrm{1}}维标准高斯分布函数；{\mathit{\beta }}^{\mathrm{L}}=\left({\beta }_{{}_{\mathrm{1}}}^{\mathrm{L}},{\beta }_{{}_{\mathrm{2}}}^{\mathrm{L}},\cdots ,{\beta }_{{}_{a_{\mathrm{1}}}}^{\mathrm{L}}\right)为极限状态方程G_i\left(\mathit{U},\mathit{\theta }\right),i=\mathrm{1,2},\cdots ,a_{\mathrm{1}}对应的最小可靠度指标{\beta }_{{}_i}^{\mathrm{L}},i=\mathrm{1,2},\cdots ,a_{\mathrm{1}}组成的向量；\mathit{\rho }为G_i\left(\mathit{U},\mathit{\theta }\right),i=\mathrm{1,2},\cdots ,a_{\mathrm{1}}的相关系数矩阵，通过计算式（19）可获得式（25）中\mathit{\rho }的每一个元素. 通过计算式（25）中的多维高斯分布可获得串联系统最大失效概率，具体可参考文献［32-33］，也可通过MATLAB软件中相应的库函数进行计算.

2.3 并联系统可靠性求解

本小节分析含有参数化p-box不确定性的并联系统（包含b₁个失效模式）. 由于存在区间分布参数，原空间中的极限状态方程映射到标准正态空间后将生成极限状态曲面带，如图3（b）所示，则并联系统曲面带的下边界和上边界可表示为：

S_{\mathrm{s}}^{\mathrm{L}}:{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cap }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}{G}_k(\mathit{U},\mathit{\theta })=\mathrm{0},S_{\mathrm{s}}^{\mathrm{U}}:{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cap }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{a}\mathrm{x}}{G}_k(\mathit{U},\mathit{\theta })=\mathrm{0}（26）

因此，并联系统失效概率P_{\mathrm{f}}为一个区间，其边界可通过求解如下两个优化获得：

（27）

（28）

为计算并联系统的最大失效概率，可将式（28）转变为：

\begin{array}{l}{P}_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1}-\mathrm{P}\mathrm{r}\left\{{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cup }}}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}{G}_k\left(\mathit{U},\mathit{\theta }\right)>\mathrm{0}\right\}=\\ \mathrm{1}-\mathrm{P}\mathrm{r}\left\{\left.{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cup }}}\left[({\beta }_k^{\mathrm{L}}-{\mathit{\alpha }}_{U,k}^{\mathrm{T}}\mathit{U})>\mathrm{0}\right]\right\}\right.=\\ \mathrm{1}-\mathrm{P}\mathrm{r}\left\{\left.{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cup }}}\left[L_k>\mathrm{0}\right]\right\}\right.\end{array}（29）

记式（29）中的L_k>\mathrm{0}为事件D_k，则式（29）可变为：

P_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1}-P\left\{{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cup }}}{D}_k\right\}（30）

根据并联事件概率计算方法，可将式（30）变为：

\left.{\left(-\mathrm{1}\right)}^{b_{\mathrm{1}}-\mathrm{1}}P\left\{{\displaystyle \underset{k=\mathrm{1}}{\overset{b_{\mathrm{1}}}{\cap }}}{D}_k\right\}\right\}（31）

其中，D_{k_{\mathrm{1}}}\bigcap D_{k_{\mathrm{2}}}\bigcap \cdots \bigcap D_{k_l}发生的概率为：

P\left\{D_{k_{\mathrm{1}}}\bigcap D_{k_{\mathrm{2}}}\bigcap \cdots \bigcap D_{k_l}\right\}=

P\left\{\left\{L_{k_{\mathrm{1}}}>\mathrm{0}\right\}\bigcap \left\{L_{k_{\mathrm{2}}}>\mathrm{0}\right\}\bigcap \dots \bigcap \left\{L_{k_l}>\mathrm{0}\right\}\right\}（32）

因为L_k~N\left({\beta }_k^{\mathrm{L}},\mathrm{1}\right)，其中{\beta }_k^{\mathrm{L}}(k=\mathrm{1,2},\cdots ,l)为G_k\left(\mathit{U},\mathit{\theta }\right)(k=\mathrm{1,2},\cdots ,l)对应的最小可靠度指标，则式（32）可转变为：

P\left\{D_{k_{\mathrm{1}}}\bigcap D_{k_{\mathrm{2}}}\bigcap \cdots \bigcap D_{k_l}\right\}={\Phi }_l({\beta }^{\mathrm{L}},\mathit{C})（33）

式中：{\Phi }_l为l维高斯累积分布函数；{\beta }^{\mathrm{L}}=\left({\beta }_{\mathrm{1}}^{\mathrm{L}},{\beta }_{\mathrm{2}}^{\mathrm{L}},\cdots ,{\beta }_k^{\mathrm{L}}\right)； C 为相关系数矩阵，可通过式（19）计算得到. 根据式（33）可获得式（31）中的所有概率项，最终求得最大失效概率P_{\mathrm{f}}^{\mathrm{m}\mathrm{a}\mathrm{x}}.

3 数值算例

本文分析了3个数值算例，算例中以系统的最大失效概率来分析基于参数化p-box模型的系统可靠度，并通过对比传统双层蒙特卡罗采样方法来阐明本文方法的精度和效率.

3.1 两单元Daniels系统算例

如图4所示为一个并联的两单元Daniels系统，该算例在文献［34］的基础上改进而来. 单元1的截面长宽分别为{\tilde{c}}_{\mathrm{1}}和{\tilde{d}}_{\mathrm{1}}，单元2的截面长宽分别为{\tilde{c}}_{\mathrm{2}}和{\tilde{d}}_{\mathrm{2}}，F为系统所受载荷，当两单元均屈服时，结构失效. 两单元对应的功能函数分别表示为：

g_{\mathrm{1}}({\tilde{c}}_{\mathrm{1}},{\tilde{d}}_{\mathrm{1}},{\tilde{\sigma }}_{\mathrm{1}},F)={\tilde{c}}_{\mathrm{1}}{\tilde{d}}_{\mathrm{1}}{\tilde{\sigma }}_{\mathrm{1}}-F/\mathrm{2}（34）

g_{\mathrm{2}}({\tilde{c}}_{\mathrm{2}},{\tilde{d}}_{\mathrm{2}},{\tilde{\sigma }}_{\mathrm{2}},F)={\tilde{c}}_{\mathrm{2}}{\tilde{d}}_{\mathrm{2}}{\tilde{\sigma }}_{\mathrm{2}}-F/\mathrm{2}（35）

式中：{\tilde{\sigma }}_{\mathrm{1}}和{\tilde{\sigma }}_{\mathrm{2}}分别为单元1和单元2的屈服强度. 考虑存在随机和认知的混合不确定性，式（34）和 式（35）中的不确定参数如表1所示.

本算例考虑了单元1与单元2失效模式相关的并联系统结构可靠度指标. 悬臂梁可靠度指标如 表2所示. 采用双层采样方法（double loop sampling， DLS）的计算结果作为参考解，表2中括号内的数为本文方法与DLS计算结果的相对误差. DLS内外层蒙特卡罗采样次数分别为10⁵次和10³次，由表2可知，本文获得的最小可靠度指标与DLS计算结果的偏差为1.17%.从计算效率来看，当对两单元的并联系统进行可靠度分析时，DLS对两功能函数的调用次数为2×10⁸次，而本文方法的调用次数为280次. 综上所述，所提方法具有很高的系统可靠性计算精度和效率.

3.2 10杆桁架

考虑10杆平面桁架结构^［35］，如图5所示. 竖直杆和水平杆的长度均为L=\mathrm{9.144}\mathrm{m}，密度\rho =\mathrm{2}\mathrm{768}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，弹性模量E=\mathrm{68}\mathrm{948}\mathrm{M}\mathrm{P}\mathrm{a}. 该桁架左端固定，受到两个垂直力F_{\mathrm{1}}、F_{\mathrm{2}}以及水平力F_{\mathrm{3}}的作用. 考虑13个不确定量，相应的分布类型及分布参数取值情况如表3所示，其中杆件的横截面积A_i(i=\mathrm{1,2},\cdots ,\mathrm{10})以及节点处所承受的外载荷F_{\mathrm{1}}、F_{\mathrm{2}}和F_{\mathrm{3}}均为参数化p-box变量.

节点2处的垂直位移d_y不能超过允许值d_{y\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{50}\mathrm{m}\mathrm{m}，最大应力发生在7杆，构件7的许用拉应力为{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{250}\mathrm{M}\mathrm{P}\mathrm{a}，所以结构的功能函数可表示为：

g_{\mathrm{1}}\left(\mathit{X}\right)=d_{y\mathrm{m}\mathrm{a}\mathrm{x}}-d_y\left(\mathit{X}\right)（36）

g_{\mathrm{2}}\left(\mathit{X}\right)={\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}-\sigma \left(\mathit{X}\right)（37）

根据平衡方程及兼容性方程，各杆的轴向力N_i,i=\mathrm{1,2},\cdots ,\mathrm{10}分别为：

\left\{\begin{array}{l}{N}_{\mathrm{1}}=F_{\mathrm{2}}-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{8}}\\ N_{\mathrm{2}}=-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{10}}\\ N_{\mathrm{3}}=-F_{\mathrm{1}}-\mathrm{2}{F}_{\mathrm{2}}+F_{\mathrm{3}}-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{8}}\\ N_{\mathrm{4}}=-F_{\mathrm{2}}+F_{\mathrm{3}}-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{10}}\\ N_{\mathrm{5}}=-F_{\mathrm{2}}-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{8}}-\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{10}}\\ N_{\mathrm{6}}=\sqrt[]{\mathrm{2}}/\mathrm{2}{N}_{\mathrm{10}}\\ N_{\mathrm{7}}=\sqrt[]{\mathrm{2}}\left(F_{\mathrm{1}}+F_{\mathrm{2}}\right)+N_{\mathrm{8}}\\ N_{\mathrm{8}}=\left(a_{\mathrm{22}}{b}_{\mathrm{1}}-a_{\mathrm{12}}{b}_{\mathrm{2}}\right)/\left(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}\right)\\ N_{\mathrm{9}}=\sqrt[]{\mathrm{2}}{F}_{\mathrm{2}}+N_{\mathrm{10}}\\ N_{\mathrm{10}}=\left(a_{\mathrm{11}}{b}_{\mathrm{2}}-a_{\mathrm{21}}{b}_{\mathrm{1}}\right)/\left(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}\right)\end{array}\right.（38）

式中：

\left\{\begin{array}{l}{a}_{\mathrm{11}}=\left(\mathrm{1}/A_{\mathrm{1}}+\mathrm{1}/A_{\mathrm{3}}+\mathrm{1}/A_{\mathrm{5}}+\mathrm{2}\sqrt[]{\mathrm{2}}/A_{\mathrm{7}}+\mathrm{2}\sqrt[]{\mathrm{2}}/A_{\mathrm{8}}\right)L/\mathrm{2}E\\ a_{\mathrm{12}}=a_{\mathrm{21}}=L/\mathrm{2}{A}_{\mathrm{5}}E\\ a_{\mathrm{22}}=\left(\mathrm{1}/A_{\mathrm{2}}+\mathrm{1}/A_{\mathrm{4}}+\mathrm{1}/A_{\mathrm{6}}+\mathrm{2}\sqrt[]{\mathrm{2}}/A_{\mathrm{10}}\right)L/\mathrm{2}E\\ b_{\mathrm{1}}=\left[F_{\mathrm{2}}/A_{\mathrm{1}}-\left(F_{\mathrm{1}}+\mathrm{2}{F}_{\mathrm{2}}-F_{\mathrm{3}}\right)/A_{\mathrm{3}}\right.-F_{\mathrm{2}}/A_{\mathrm{5}}-\\ \left.\mathrm{2}\sqrt[]{\mathrm{2}}\left(F_{\mathrm{1}}+F_{\mathrm{2}}\right)/A_{\mathrm{7}}\right]\sqrt[]{\mathrm{2}}L/\mathrm{2}E\\ b_{\mathrm{2}}=\left[\sqrt[]{\mathrm{2}}\left(F_{\mathrm{3}}-F_{\mathrm{2}}\right)/A_{\mathrm{4}}-\sqrt[]{\mathrm{2}}{F}_{\mathrm{2}}/A_{\mathrm{5}}-\mathrm{4}{F}_{\mathrm{2}}/A_{\mathrm{9}}\right]L/\mathrm{2}E\end{array}\right.（39）

最大位移d_y和构件7的许用应力{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}分别如下：

d_y=\left[{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}\frac{N_i^{\mathrm{0}}{N}_i}{A_i}+\sqrt[]{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{7}}^{\mathrm{10}}}\frac{N_i^{\mathrm{0}}{N}_i}{A_i}\right]\frac{L}{E}（40）

{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{\sqrt[]{\mathrm{2}}\left(F_{\mathrm{1}}+F_{\mathrm{2}}\right)\left(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}\right)+a_{\mathrm{22}}{b}_{\mathrm{1}}-a_{\mathrm{12}}{b}_{\mathrm{2}}}{A_{\mathrm{7}}\left(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}\right)}（41）

式中：N_i^{\mathrm{0}}由式（38）令F_{\mathrm{1}}=F_{\mathrm{3}}=\mathrm{0}，F_{\mathrm{2}}=\mathrm{1}求得.

考虑双失效模式串联系统的10杆桁架结构可靠度指标，如表4所示，表4中括号内的数为本文方法与DLS计算结果的相对误差. 由表4可知，本文获得的考虑位移失效模式、杆7应力失效模式串联系统的最小可靠度指标与DLS计算结果的偏差为1.54%. 此外，从计算效率来看，当对双失效模式系统进行可靠性分析时，DLS对两功能函数的调用次数为2×10⁸次，而本文方法的调用次数为198次，体现了本文方法具有很高的计算精度和效率.

3.3 平板电脑

本算例考虑一个平板电脑的系统可靠性问题，如图6所示，该平板电脑主要包括如下构件：后壳体、主板、前壳体、触摸屏、显示屏、支架、电池. 涉及如下3种极限工况：环境温度变化过程、自由跌落以及高温环境.

工况1：在环境温度变化过程（0~40 ℃）中，因各种材料热膨胀系数不一致而引起的热应力可能造成电子元器件失效，考虑电池的最大热应力{\Gamma }_{}^{\mathrm{B}\mathrm{A}}不应超过额定值{\Gamma }_{\mathrm{0}}^{\mathrm{B}\mathrm{A}}=24 MPa，以保证电池在使用过程中的安全性.

工况2：平板电脑在自由跌落的撞击过程中，触摸屏上的最大应力{\Gamma }_{}^{\mathrm{T}\mathrm{L}}不应超过其材料的破裂强度{\Gamma }_{\mathrm{0}}^{\mathrm{T}\mathrm{L}}=\mathrm{900}\mathrm{M}\mathrm{P}\mathrm{a}；自由跌落的高度设定为0.5 m，平板电脑正面向下进行撞击.

工况3：在高温环境下，平板电脑的工作温度设置45 ℃；主板上的芯片温度T_{}^{\mathrm{C}\mathrm{H}}不能超过其许用的工作温度T_{\mathrm{0}}^{\mathrm{C}\mathrm{H}}=\mathrm{70}℃.

平板电脑的不确定参数分布如表5所示，3种工况的有限元模型如图7所示，平板电脑的有限元模型信息如表6所示. 为提升效率，对各仿真模型分别采样60次用以建立二次响应面，平板电脑的功能函数响应面如表7所示.

表8为平板电脑可靠度指标.由表8可知，本文获得的考虑3种工况下失效模式串联的结构系统最小可靠度指标为2.432，与其DLS计算结果对比，相对误差为0.61%.此外，从计算效率来看，当对失效模式串联系统进行可靠度分析时，DLS方法对功能函数的调用次数为3×10⁸次，而本文方法需要调用功能函数的次数为241次，体现了本文方法不仅能获得较精确的计算结果而且具有很高的效率.

4 结 论

本文针对含参数化p-box不确定性的结构系统问题，提出了一种有效的系统可靠性分析方法. 首先，基于序列迭代解耦方法获得单失效模式下结构的最小可靠度指标；然后建立了含参数化p-box不确定性的系统可靠性模型；考虑系统各失效模式之间的相关性，通过在最大可能失效点处固定区间分布参数并结合线性相关度方法计算各失效模式间的相关系数；最后提出了考虑参数化p-box不确定性的串联系统和并联系统可靠性计算方法. 本文方法能够处理一类存在随机和认知混合不确定性的复杂机械系统可靠性问题. 3个数值算例分析表明，通过对比传统双层蒙特卡罗采样方法的计算结果，本文方法具有较好的计算精度和效率.

参考文献

[1] TANG J C， LI X， FU C M，et al ．A possibility-based solution framework for interval uncertainty-based design optimization［J］．Applied Mathematical Modelling，2024，125： 649-667．

[2] LIU J F， JIANG C， LIU H B，et al ．A time variant uncertainty propagation method for high-dimensional dynamic structural system via K–L expansion and Bayesian deep neural network［J］．Philosophical Transactions of the Royal Society A：Mathematical，Physical and Engineering Sciences，2023，381（2260）： 20220388．

[3] 吕震宙，岳珠峰 ．模糊随机可靠性分析的统一模型［J］．力学学报，2004，36（5）：533-539．

[4] LÜ Z Z， YUE Z F. Unified reliability model for fuzziness and randomness of the basic variables and state variables in structure［J］. Chinese Journal of Theoretical and Applied Mechanics， 2004， 36（5）： 533-539．（in Chinese）

[5] TANG J C， FU C M， MI C J，et al ．An interval sequential linear programming for nonlinear robust optimization problems［J］．Applied Mathematical Modelling，2022，107：256-274．

[6] LIU H B， JIANG C， XIAO Z ．Efficient uncertainty propagation for parameterized p-box using sparse-decomposition-based polynomial chaos expansion［J］．Mechanical Systems and Signal Processing，2020，138：106589．

[7] LIU H B， CHEN M， DU C，et al ．A copula-based uncertainty propagation method for structures with correlated parametric p-boxes［J］．International Journal of Approximate Reasoning，2021，138：89-104．

[8] ZHANG H， MULLEN R L， MUHANNA R L ．Interval Monte Carlo methods for structural reliability［J］．Structural Safety，2010，32（3）：183-190．

[9] XIAO Z， HAN X， JIANG C，et al ．An efficient uncertainty propagation method for parameterized probability boxes［J］．Acta Mechanica，2016，227（3）：633-649．

[10] CHEN N， YU D J， XIA B Z，et al ．Uncertainty analysis of a structural-acoustic problem using imprecise probabilities based on p-box representations［J］．Mechanical Systems and Signal Processing，2016，80：45-57．

[11] WEI P F， SONG J W， BI S F，et al ．Non-intrusive stochastic analysis with parameterized imprecise probability models：Ⅱ．reliability and rare events analysis［J］．Mechanical Systems and Signal Processing，2019，126：227-247．

[12] CRESPO L G， KENNY S P， GIESY D P ．Reliability analysis of polynomial systems subject to p-box uncertainties［J］．Mechanical Systems and Signal Processing，2013，37（1/2）：121-136．

[13] 刘鑫，何泽波，周振华，等 ．基于概率-概率盒混合模型的气囊座椅防护特性可靠性分析方法［J］．机械工程学报， 2022， 58（24）： 324-333．

[14] LIU X， HE Z B， ZHOU Z H，et al ．Hybrid reliability approach for airbag seat protection performance based on probability and probability box models［J］．Journal of Mechanical Engineering，2022，58（24）： 324-333．（in Chinese）

[15] 黄洪钟，刘征，米金华，等 ．混合不确定性下机床主轴可靠性建模与分析［J］．中国科学（物理学 力学 天文学），2018，48（1）：42-53．

[16] HUANG H Z， LIU Z， MI J H，et al ．Reliability modeling and analysis of heavy-duty CNC machine tool spindle under hybrid uncertainty［J］. Scientia Sinica （Physica，Mechanica & Astronomica），2018，48（1）：42-53．（in Chinese）

[17] WANG C， ZHANG H， BEER M. Computing tight bounds of structural reliability under imprecise probabilistic information［J］．Computers & Structures，2018，208：92-104．

[18] 丁家满， 原琦， 任东磊， 等. 不确定性信息条件下系统可靠性分析［J］．仪器仪表学报，2019，40（4）：153-162．

[19] DING J M， YUAN Q， REN D L，et al ．System reliability analysis under uncertain information［J］．Chinese Journal of Scientific Instrument，2019，40（4）：153-162．（in Chinese）

[20] 项涌涌， 潘柏松， 罗路平， 等. 面向响应准确度的参数不确定性模型确认方法［J］. 中国机械工程， 2019， 30（7）： 811-817．

[21] XIANG Y Y， PAN B S， LUO L P，et al ．Validation method of parameter uncertainty models for response accuracy［J］．China Mechanical Engineering， 2019， 30（7）： 811-817．（in Chinese）

[22] LIU X， WANG X Y， XIE J，et al ．Construction of probability box model based on maximum entropy principle and corresponding hybrid reliability analysis approach［J］. Structural and Multidisciplinary Optimization， 2020， 61（2）： 599-617．

[23] ZHANG J H， GAO L， XIAO M，et al ．An active learning Kriging-assisted method for reliability-based design optimization under distributional probability-box model［J］. Structural and Multidisciplinary Optimization，2020，62（5）： 2341-2356．

[24] XIE H C， LI J W， LIU H B，et al. An efficient uncertainty propagation analysis method of non-parameterized P-boxes based on dimension-reduction integral and maximum entropy estimation［J］. Structural and Multidisciplinary Optimization，2023， 66（12）： 258．

[25] XIAO Z， ZHANG Q W， ZHANG Z，et al ．A collaborative quasi-Monte Carlo uncertainty propagation analysis method for multiple types of epistemic uncertainty quantified by probability boxes［J］．Structural and Multidisciplinary Optimization，2023，66（5）：109．

[26] SALOMON J， WINNEWISSER N， WEI P F，et al ．Efficient reliability analysis of complex systems in consideration of imprecision［J］．Reliability Engineering & System Safety，2021，216： 107972．

[27] FAES M G R， VALDEBENITO M A， MOENS D，et al ．Operator norm theory as an efficient tool to propagate hybrid uncertainties and calculate imprecise probabilities［J］．Mechanical Systems and Signal Processing，2021，152：107482．

[28] KARL B. Asymptotic approximations for multinormal integrals［J］．Journal of Engineering Mechanics， 1984， 110（3）： 357-366．

[29] 刘海波，姜潮，郑静，等.含概率与区间混合不确定性的系统可靠性分析方法［J］. 力学学报， 2017， 49（2）： 456-466.

[30] LIU H B， JIANG C， ZHENG J， et al. A system reliability method for structures with probability and interval mixed uncertainty［J］.Chinese Journal of Theoretical and Applied Mechanics， 2017， 49（2）： 456-466. （in Chinese）

[31] JIANG C， LI W X， HAN X，et al ．Structural reliability analysis based on random distributions with interval parameters［J］．Computers & Structures， 2011， 89（23/24）： 2292-2302．

[32] DU X P， SUDJIANTO A， HUANG B Q ．Reliability-based design with the mixture of random and interval variables［J］．Journal of Mechanical Design， 2005， 127（6）： 1068-1076．

[33] ZHANG Y， KIUREGHIAN A ．Two improved algorithms for reliability analysis［M］//Reliability and Optimization of Structural Systems．Boston，MA：Springer US，1995：297-304．

[34] 韩文钦，周金宇，孙奎洲 ．失效模式相关的机械结构可靠性的Copula分析方法［J］．中国机械工程， 2011， 22（3）： 278-282．

[35] HAN W Q， ZHOU J Y， SUN K Z ．Copula analysis of structural systems reliability with correlated failure mode［J］．China Mechanical Engineering， 2011， 22（3）： 278-282．（in Chinese）

[36] ADDURI P R， PENMETSA R C ．Bounds on structural system reliability in the presence of interval variables［J］．Computers & Structures，2007，85（5/6）：320-329．

[37] 张明. 结构可靠度分析：方法与程序［M］. 北京：科学出版社， 2009： 130-136.

[38] ZHANG M. Structural reliability analysis：methods and procedures［M］. Beijing： Science Press， 2009：130-136. （in Chinese）

[39] 贡金鑫 ．工程结构可靠度计算方法［M］．大连：大连理工大学出版社，2003： 234-268．

[40] GONG J X ．Computational methods for reliability of engineering structures［M］．Dalian：Dalian University of Technology Press，2003： 234-268．（in Chinese）

[41] GENZ A， TRINH G. Numerical computation of multivariate normal probabilities using bivariate conditioning［J］. Journal of Computational & Graphical Statistics， 1992， 1（2）： 141-149.

[42] TANG L K， MELCHERS R E ．Improved approximation for multinormal integral［J］．Structural Safety，1986，4（2）：81-93．

[43] HU Z， DU X P ．First order reliability method for time-variant problems using series expansions［J］. Structural and Multidisciplinary Optimization， 2015， 51（1）： 1-21．

[44] GANZERLI S， PANTELIDES C P ．Optimum structural design via convex model superposition［J］．Computers & Structures，2000，74（6）：639-647．

基金资助

国家自然科学基金资助项目（52205263,51975577）; 中国博士后科学基金资助项目（2021M690988）