结构失效概率分析的球层抽样方法

白冰

doi:10.3969/j.issn.0258-2724.20210848

摘要:

传统Monte Carlo抽样方法应用于小失效概率等复杂可靠度问题时，存在效率低下、精度有限等不足，针对这一问题提出一种球层抽样分析方法. 首先，通过分离距离与方向参量，对标准正态随机向量进行重构，并验证其标准正态性与相互独立性；然后，采用分层抽样策略，将半径大于一阶可靠度指标外的正态空间划分为多个球层，进而利用所重构的向量逐层进行抽样，并结合全概率公式，形成估计结构失效概率的球层抽样算法；最后，通过对3个典型算例进行对比分析，验证算法性能. 分析表明：所提出算法具有较高的抽样效率与收敛性能，算例计算结果误差在3%以内；与其他算法相比，其估计方差更小，且可有效解决多设计验算点等复杂可靠度问题；算法在抽样效率、适用范围以及稳定性等方面具有优势，更适用于实际复杂结构可靠度的求解分析.

关键词:

Abstract:

When the traditional Monte Carlo sampling method is applied to complex reliability problems such as small failure probability, there are some shortcomings such as low efficiency and limited accuracy. To solve this problem, a spherical layer sampling analysis method is developed. Firstly, by dividing the distance and direction parameters, the standard normal random vector is reconstructed, and its standard normality and mutual independence are verified. Thereafter, based on a layered sampling strategy, the standard normal space with the radius beyond first order reliability index is divided into multiple spherical layers, which are then sampled by the reconstructed vector layer by layer. Combined with the full probability formula, a spherical layer sampling algorithm is developed to estimate the structural failure probability. Finally, three typical examples are taken as objects of interest, and the performance of the algorithm is verified through comparative analysis. The results show that, the proposed algorithm has high sampling efficiency and convergence performance, and the error of calculation results is within 3%. Compared with other algorithms, its estimation variance is smaller, and it can effectively solve complex reliability problems such as multiple design check points. The algorithm has advantages in sampling efficiency, scope of application, and stability, and is more suitable for solving and analyzing the reliability of actual complex structures.

Key words:

现代大型桥梁和结构构造组成复杂，对结构安全性要求高，采用基于概率的结构可靠度理论进行安全评估已渐成共识. 然而，由于结构本身及随机影响因素的复杂性，传统的一阶（FORM）和二阶（SORM）可靠度计算方法^[1]在处理复杂非线性功能函数时精度较差，适用性有限. 尤其是对于一些重要结构，如核电站保护壳和生命线工程等，其失效概率计算精度要求较高，传统方法已较难满足相应要求. 为解决这一问题，基于抽样模拟技术的Monte Carlo法由于其独特优势而获得了长足发展.

Monte Carlo法以大数定律为理论基础，通过数值抽样来估计结构失效概率，并以抽样模拟频率代替失效概率，具有对随机变量维数不敏感的优势. 然而，面对小失效概率问题时，为确保计算精度，需要增加抽样次数，这是当前Monte Carlo法主要面临的挑战. 为解决这个问题，许多学者提出了各种抽样方差缩减技术.

Melchers^[2]基于Kahn的重要抽样思想，以设计验算点为重心，使用n维独立正态分布密度函数作为重要抽样密度进行抽样估计，以获取结构的失效概率，但对重要抽样密度函数参数的确定未进行深入研究. Geyer等^[3]基于Rubinstein的交叉熵方法对重要抽样密度函数进行优化探讨，取得了良好进展. 线抽样方法是另一类有效的高效抽样方法，由Koutsourelakis等^[4-5]提出，通过适当变换结构失效概率估计函数并沿重要方向进行抽样，结合逐步计算策略，进一步提升算法的方差缩减程度，具有较好的估计精度和适用性，但是算法同样也存在一些问题. 吕震宙等^[1,6]提出解决方案并取得了良好效果. 不同于上述方法，Au 等^[7-8]提出一种新的子集抽样法，通过定义一系列中间失效事件并结合Markov chain Monte Carlo （MCMC）方法，将一个小失效概率问题转化为多个较大条件失效概率的乘积来缩减抽样方差. 这一思想被Katafygiotis等^[9]借鉴，研发了一种球层子集抽样模拟方法. 该方法将标准正态空间划分为若干互斥的球层子区域，在引入各球层内锥形失效域空间由外向内逐渐递减的假定基础上，将各子区域失效概率结果转化为多个参数（每一参数相当于一个中间失效事件）的乘积，并通过MCMC条件抽样方法对各层锥形失效域进行抽样搜索以估计相应参数，实现结构失效概率的高效模拟. Katafygiotis等^[10]进一步优化球层子区域的选取，提高了算法效率. 研究表明，其提出的方法具有出色的方差缩减能力和估计精度，适用于小失效概率等问题. 然而，该方法主要针对一般失效域问题，并对失效域形态提出了一定要求，对于多失效域等复杂可靠度问题，其所采用的Markov链搜索方法可能存在困难，并且MCMC方法具有样本相关性高等缺陷. 因此，需要进一步研究改进，以解决复杂情形的可靠度问题.

总体而言，上述代表性抽样方法均具有较高的计算效率和估计精度，尤其在处理复杂高维、小失效概率问题时展现出独特优势和较好适用性，推动了Monte Carlo方法的发展. 但是对于多设计验算点（多重要失效域）问题，由于这些方法均具有一定的方向侧重性，可能较难解决这一问题. 因此，本文提出一种新的结构失效概率估计的球层抽样方法. 首先，对多维随机向量进行重新构造变换，保证其标准正态性与相互独立性不变；然后，将标准正态空间外围区域划分为多个同心球层，并通过重点对这些球层进行均匀无方向抽样，来缩减抽样方差并提升精度效率. 这种方法具有良好的估计性能，且由于抽样的均匀无方向性，可有效解决前述多设计验算点问题，更适用于复杂小失效概率和体系可靠度等问题的求解.

1. 抽样向量的构造及其正态性与独立性研究

1.1 抽样向量构造

针对结构随机变量，通常可以通过正态化和线性变换将其转化至标准正态空间中进行分析^[11]. 经处理后，可以得到标准正态空间的n维基本向量X=[X₁X₂… X_n]. 为构建后续球层抽样方法，在这个空间中，借鉴Katafygiotis等^[9]的做法，对抽样向量进行改写构造.

在标准正态空间中，随机向量X的参考坐标系为直角坐标系，因此，可以使用一般形式描述随机抽样点x，也可以使用其与原点的距离R及单位方向向量α进行描述^[1]，如式（1）所示. 这种向量描述形式对于后续抽样方法的建立具有较多方便之处. 二维向量示意如图1所示.

图 1 标准正态空间中样本点

Figure 1. Sampling point in standard normal space

下载: 全尺寸图片幻灯片

${\boldsymbol{x}}{\text{ = }}R{\boldsymbol{\alpha}} {\text{ = }}\left\| {\boldsymbol{x}} \right\|\left(\frac{{{x_1}}}{{\left\| {\boldsymbol{x}} \right\|}},\frac{{{x_2}}}{{\left\| {\boldsymbol{x}} \right\|}}, \cdot \cdot \cdot ,\frac{{{x_n}}}{{\left\| {\boldsymbol{x}} \right\|}}\right),$

(1)

$R{\text{ = }}\left\| {\boldsymbol{x}} \right\| = \sqrt {x_1^2 + x_2^2 + \cdot \cdot \cdot + x_n^2},$

(2)

${\boldsymbol{\alpha}} {\text{ = (}}{{{x_1}} / {\left\| {\boldsymbol{x}} \right\|}}{\text{, }}{{{x_2}} / {\left\| {\boldsymbol{x }}\right\|}}, \cdot \cdot \cdot {\text{, }}{{{x_n}} / {\left\| {\boldsymbol{x}} \right\|}}{\text{)}} ,$

(3)

式中：x= (x₁, x₂,…, x_n).

根据上述描述，结构随机向量可仿照式（1）～（3）进行构造.

对于距离R，定义随机变量L为

$L{\text{ = }}Q_1^2 + Q_2^2 + \cdots + Q_i^2 \cdots + Q_n^2{\text{ }},$

(4)

式中：Q_i～N(0, 1)，L服从χ²(n)分布^[12].

对于单位方向向量α，定义随机向量ξ为

${\boldsymbol{ \xi }}{\text{ = (}}{{{X_1}} / {\left\| {\boldsymbol{X}} \right\|}}{\text{, }}{{{X_2}} / {\left\| {\boldsymbol{X}} \right\|}}, \cdot \cdot \cdot {\text{, }}{{{X_n}} / {\left\| {\boldsymbol{X}} \right\|}}{\text{) ,}}$

(5)

式中：X_i～N(0, 1).

由于引入了归一化约束，因而ξ自由度降为n−1维，ξ服从n−1维球面上的均匀分布.

将随机变量L开方并与 $\boldsymbol{\xi}$ 相乘，以引入表征距离影响的随机因素，可得新构造向量Z为

$\begin{split} & {\boldsymbol{Z}} = \sqrt L {\boldsymbol{\xi}} {\text{ = (}}{{\sqrt L {X_1}} / {\left\| {\boldsymbol{X}} \right\|}}{\text{, }}{{\sqrt L {X_2}} / {\left\| {\boldsymbol{X}} \right\|}}, \cdots {\text{, }}\\ &\quad {{\sqrt L {X_n}} / {\left\| {\boldsymbol{X}} \right\|}}{\text{) }}{\text{.}} \end{split}$

(6)

由式（6）可以看出，该向量除ξ外，新增加了一维随机变量L. 因而，其维度与n维标准正态空间维度相等，且Z中各分量Z_i= ${{\sqrt L {X_i}} \mathord{\left/ {\vphantom {{\sqrt L {X_i}} {\left\| X \right\|}}} \right. } {\left\| {\boldsymbol{X}} \right\|}}$ 保持同分布. 其抽样频率分布与标准正态分布对比如图2所示，图中，n=100，抽样次数为10⁵次. 可以发现，Z_i分布特性与标准正态分布概率密度高度一致.

图 2 Z_i抽样频率密度与标准正态分布概率密度对比

Figure 2. Fig. 2 Comparison of probability density between Z_i sampling and standard normal distribution

下载: 全尺寸图片幻灯片

Katafygiotis等^[9]指出向量Z服从n维独立的标准正态分布，但并未给出严格证明，而这是后续算法研究的重要基础. 因此，有必要对其特性进行验证. 下面将对随机向量Z的标准正态性与相互独立性进行验证.

1.2 向量标准正态性验证

由式（6）定义可知，随机向量Z为 $\sqrt L$ 与 ${\boldsymbol{\xi}}$ 的乘积，且二者相互独立. 由于 ${\boldsymbol{\xi}}$ 服从n−1维自由度球面上的均匀分布，因此其分布域实质上是一个以原点为中心、半径为1的多维球面 $S_1^{(n-1)}$ ，如图3所示. 球面上各点概率密度h(ξ)（式（7））相同.

图 3 n维标准正态空间中n−1维球面均匀分布

Figure 3. Spherical uniform distribution of n−1 dimensions in standard normal space

下载: 全尺寸图片幻灯片

$h({\boldsymbol{\xi}} ){\text{ = }}\frac{1}{{\displaystyle\int_{S_1^{(n - 1)}} {{\text{ }}{\mathrm{d}}{A_1}} }}.$

(7)

对式（7）进行变形，考虑标准正态随机向量旋转不变性，则

$\begin{split} & h({\boldsymbol{\xi}} )\text=\frac{\phi ({\boldsymbol{x}}|{\boldsymbol{x}}\in {S}_{r}^{(n-1)})}{{\displaystyle {\int }_{{\boldsymbol{\xi}} \in {S}_{1}^{(n-1)}}\phi ({\boldsymbol{x}}|{\boldsymbol{x}}\in {S}_{r}^{(n-1)}){\mathrm{d}}{A}_{1}}}=\\ &\quad \left\{\begin{array}{l}\dfrac{\phi ({\boldsymbol{x}})}{{\displaystyle {\int }_{{\boldsymbol{\xi}} \in {S}_{1}^{(n-1)}}\phi ({\boldsymbol{x}}){\mathrm{d}}{A}_{1}}}\text{, }\quad \Vert {\boldsymbol{x}}\Vert =r,\\ \text{ }0\text{ , }\quad 其他\text{,}\end{array}\right.\\ \end{split}$

(8)

式中：ϕ(•)为标准正态分布概率密度函数， S_r⁽ⁿ⁻¹⁾为n维标准正态空间中半径为r的多维球面（r可取大于0的任意值）.

由于h(ξ)仅在||x||=r上取值非0，因此，仅对式（8）中||x||=r进行分析，通过对半径r遍历即可覆盖整个标准正态空间. 结合式（5）可得

$\begin{split} & h({\boldsymbol{\xi}} ) = \frac{{\phi ({\boldsymbol{x}})}}{{\displaystyle\int_{r{\boldsymbol{\xi}} \in S_r^{(n - 1)}} {\phi ({\boldsymbol{x}})/{r^{n - 1}}{\mathrm{d}}{A_r}} }} = \frac{{\phi ({\boldsymbol{x}})}}{{\displaystyle\int_{{\boldsymbol{x}} \in S_r^{(n - 1)}} {\phi ({\boldsymbol{x}})/{r^{n - 1}}{\mathrm{d}}{A_r}} }} = \\ &\quad \frac{{\phi ({\boldsymbol{x}})}}{{P({r^2} \leqslant {{\boldsymbol{x}}^2} \leqslant {r^2} + {\mathrm{d}}{r^2})/({r^{n - 1}}{\mathrm{d}}r)}}{\text{ = }} \\ &\quad \frac{{\phi ({\boldsymbol{x}})}}{{{\mathrm{d}}F({r^2})/({r^{n - 1}}{\mathrm{d}}r)}} = \frac{{\phi ({\boldsymbol{x}})}}{{2rf({r^2}) /{r^{n - 1}}}} , \end{split}$

(9)

式中：F(•)和f(•)分别为χ²(n)分布的累积分布函数和概率密度函数.

整理式（9）可得

$\phi ({\boldsymbol{x}})= 2rh(\xi )f({r^2})/{r^{n - 1}}{\text{ .}}$

(10)

对式（10）两边进行积分，则

$\begin{split} & \int_V {\phi ({\boldsymbol{x}}){\mathrm{d}}{\boldsymbol{x}}} {\text{ = }}\int_{S_r^{(n - 1)}} {h({\boldsymbol{\xi}} )/{r^{n - 1}}} {\mathrm{d}}{A_r}\int_r {2rf({r^2})} {\mathrm{d}}r{\text{ = }} \\ &\quad \int_{S_1^{(n - 1)}} {h({\boldsymbol{\xi}} )} {\mathrm{d}}{A_1}\int_l {f(l){\mathrm{d}}l} , \end{split}$

(11)

式中：l为L～χ²(n)的抽样样本点，V为标准正态空间；右侧两项分别对应于随机向量ξ以及变量L.

对比式（6）、（11）可以看出，向量Z= $\sqrt L$ ξ概率密度分布与n维标准正态随机向量相同，其标准正态性由此得证.

事实上，向量Z的标准正态性亦可通过其矩信息来反映^[13]. 由于Z中各变量Z_i与Z_j （i≠j）始终保持同分布，因此，只需对其中任一变量Z_i进行分析即可. 考虑到理论计算的难度与复杂性，采用数值抽样方法分别计算不同维度n下Z_i的前四阶矩（均值、方差、偏度、峰度），其结果如图4所示.

图 4 不同维度下Z_i前四阶矩抽样结果

Figure 4. Sampling results of first four-order moments of Z_i in different dimensions

下载: 全尺寸图片幻灯片

由图4可以看出，各阶矩抽样结果与标准正态分布理论结果符合良好，其抽样值均收敛至精确值，这从另一方面也说明了向量Z的标准正态性.

1.3 变量独立性验证

在构造的向量Z中，为确保抽样的正确性，除需验证向量的标准正态性之外，还需要确认向量中各个变量的独立性.

由于Z服从标准正态分布，根据正态分布特性可知，检验其独立性可转化为检验变量间相关性^[13]. 若变量线性不相关，即等价于其相互独立，由此，计算变量间相关系数为

$\begin{split} & \rho ({Z_i},{\text{ }}{Z_j}) = \frac{{{\text{cov}}({Z_i},{\text{ }}{Z_j})}}{{\sqrt {D({Z_i})} \sqrt {D({Z_j})} }} = {\text{cov}}({Z_i},{\text{ }}{Z_j})=\\ &\quad E({Z_i}{Z_j}) - E({Z_i})E({Z_j}){\text{ = }}E({Z_i}{Z_j}) = \\ &\quad E(L{{{X_i}{X_j}} / {{{\left\| {\boldsymbol{X}} \right\|}^2}}}) = E(L)E({{{X_i}{X_j}} / {{{\left\| {\boldsymbol{X}} \right\|}^2}}}), \end{split}$

(12)

式中：D（•）为方差.

X_i和X_j独立同分布，且均服从N(0, 1)，则由对称性可知

$E({{{X_i}{X_j}} / {{{\left\| {\boldsymbol{X}} \right\|}^2}}}){\text{ = }}0{\text{ }}.$

(13)

将式（13）代入式（12），则

$\rho ({Z_i},{\text{ }}{Z_j}){\text{ = cov}}({Z_i},{\text{ }}{Z_j}) = 0{\text{ }}.$

(14)

根据这一结果，对于向量Z中任意不同变量，可证其相关系数ρ=0，即各变量均相互独立.

结合前述正态性验证结果可知，向量Z= $\sqrt L$ ξ为标准正态向量，且变量间相互独立. 由这一特性，Z与基本向量X在标准正态空间中实质为等价向量. 但与向量X相比，对Z进行抽样操作更加灵活，其组成变量L控制了样本点相对原点之间的距离，而ξ控制了样本点相对原点的角度方位，二者因素相互剥离. 这为后续抽样方法的构造提供了基础.

2. 球层抽样算法

2.1 抽样方法原理及策略

在标准正态空间中，结构的一阶可靠度指标β_form定义了结构失效域相对原点最近的距离^[11]. 因此，对于各种极限状态方程形式，以β_form为半径的多维球形空间均包含于结构的可靠域内，且为可靠域的空间下限^[14]，如图5所示. 图中，g(X)为功能函数.

图 5 可靠域中多维球形空间

Figure 5. Mlti-dimensional spherical space in reliable domain

下载: 全尺寸图片幻灯片

对于半径大于β_form的正态空间，其空间特点使其成为重点关注区域. 针对这部分空间，借鉴Gong等^[15]方法的思想，对其进行分层处理. 通过将半径大于β_form外正态空间按球形逐层划分，可以得到若干同心球层区域，如图6所示. 此时，每个球层中的部分区域都会落入失效域，且相对各球层子空间，这些区域的条件失效概率会显著增加. 这一特点为缩减抽样方差、克服小失效概率抽样效率低的问题提供了解决思路. 根据这一思路，结构失效概率可表达为

图 6 标准正态空间按球层划分

Figure 6. Spherical layer division in standard normal space

下载: 全尺寸图片幻灯片

$\begin{split} & {P_{\text{f}}}{\text{ = }}P(g({\boldsymbol{X}}) < 0){\text{ = }} \\ &\quad P(g({\boldsymbol{X}}) < 0|{{\boldsymbol{X}}^2} \leqslant \beta _{{\text{form}}}^2) P({{\boldsymbol{X}}^2} \leqslant \beta _{{\text{form}}}^2) + \\ &\quad \sum_k {P(g({\boldsymbol{X}}) < 0|\eta _k^2 \leqslant {{\boldsymbol{X}}^2} \leqslant \eta _{k + 1}^2) P(\eta _k^2 \leqslant {{\boldsymbol{X}}^2} \leqslant \eta _{k + 1}^2)}= \\ &\quad \sum_k {P(g({\boldsymbol{X}}) < 0|\eta _k^2 \leqslant {{\boldsymbol{X}}^2} \leqslant \eta _{k + 1}^2) P(\eta _k^2 \leqslant {{\boldsymbol{X}}^2} \leqslant \eta _{k + 1}^2)} , \end{split}$

(15)

式中：P(•)为概率；η_k为球层k与球层k−1的交界半径，k=1, 2, …, K；X²服从χ²(n)分布.

故式（15）可化为

${P_{\text{f}}}{\text{ = }}\sum_k {P(g({\boldsymbol{X}}) < 0|\eta _k^2 \leqslant {{\boldsymbol{X}}^2} \leqslant \eta _{k + 1}^2) [F(\eta _{k + 1}^2) - F(\eta _k^2)]}.$

(16)

式（16）中g(X)<0| $\eta_k^2$ ≤X²≤ $\eta_{k+1}^2$ 区域并非一定是一个连续锥形区域，这一点对于多失效域可靠度问题尤为重要.

根据第1节向量构造与验证结果可知，所构造向量Z与基本向量X在标准正态空间中分布等价. 因此，将Z= $\sqrt L$ ξ = $\sqrt L$ X/||X||代替X代入式（16）中，结果如式（17）所示.

$\begin{split} & P_{\mathrm{f}}= \sum_k P\left(g(\boldsymbol{Z})<0 \mid \eta_k^2 \leqslant L \leqslant \eta_{k + 1}^2\right) \left(F\left(\eta_{k + 1}^2\right)- F\left(\eta_k^2\right)\right)= \sum_k \iint I(g(\sqrt{l} \boldsymbol{x} /\|\boldsymbol{x}\|)) \phi(\boldsymbol{x}) f\left(l \mid \eta_k^2 \leqslant l \leqslant \eta_{k + 1}^2\right) {\mathrm{d}} \boldsymbol{x} {\mathrm{d}} l \;\times \\ &\quad \left(F\left(\eta_{k + 1}^2\right)-F\left(\eta_k^2\right)\right)= \sum_k E\left(I(g(\sqrt{l} \boldsymbol{x} /\|\boldsymbol{x}\|)) \mid \eta_k^2 \leqslant l \leqslant \eta_{k + 1}^2\right) \left(F\left(\eta_{k + 1}^2\right)-F\left(\eta_k^2\right)\right), \end{split}$

(17)

式中：I(g(z))为示性函数，如式（18）所示， ${\boldsymbol{z}} = \sqrt{l}{\boldsymbol{x}}/||{\boldsymbol{x}}||$ .

$I(g({\boldsymbol{z}})) = \left\{ \begin{gathered} 1，{\text{ }}g({\boldsymbol{z}}) < 0 , \\ 0，{\text{ }}g({\boldsymbol{z}}) \geqslant 0. \\ \end{gathered} \right.$

(18)

由此，失效概率P_f可通过抽样估计为^[10]

$\begin{split} & P_{\mathrm{f}} \approx P_{\text {estimate }}= \sum_k\left[F\left(\eta_{k + 1}^2\right)-F\left(\eta_k^2\right)\right]\times \\ &\quad \sum_m \frac{I\left(\sqrt{l^{(m)}} \boldsymbol{x}^{(m)} \Big/\left\|\boldsymbol{x}^{(m)}\right\| \in \boldsymbol{F}^* \Big| \eta_k^2 \leqslant l \leqslant \eta_{k + 1}^2\right)}{N_k}, \end{split}$

(19)

式中：P_estimate为失效概率P_f的抽样估计结果，F^*为功能函数g(z)对应的失效域，l^(m)和x^(m)分别为l和x的第m个抽样样本，m=1, 2, …, N_k.

式（19）即为结构失效概率的抽样理论公式. 可以看出，所建立抽样方法的一项前提在于如何将抽样点按所服从分布均匀控制在所要求半径η_k≤R≤η_k+1的球层中，而这正是前述构造向量Z= $\sqrt L$ X/||X||的目的所在. 通过对基本向量Z的重构，可将表征距离的变量L和表征方向的向量X/||X||分离开来，从而实现基于特定球层的分析处理.

对球层[η_k, η_{k + 1}]空间内样本获取采用如下方法：对L在概率区间[F( $\eta_k^2$ ), F( $\eta_{k+1}^2$ )]内抽样，生成随机样本l，在此基础上随机抽取标准正态分布样本x，构造方向向量x/||x||^[10]，得到球层内随机抽样点z= $\sqrt{l}$ x/||x||. 该样本即为式（19）估计结构失效概率P_f的随机抽样点，如图7所示.

图 7 各球层内随机样本点生成

Figure 7. Generation of random sample points in each spherical layer

下载: 全尺寸图片幻灯片

对于球层的划分，由于半径R = β_form的多维球形空间是结构可靠域的空间下限，因此，仅需对R> β_form的空间进行划分. 这种方式相比其他方法^[15]既可以确保不遗漏失效域，又能有效降低抽样工作量. 此时可按照等概率原则合理确定各球层，该原则相对于部分研究中所建议的预抽样反复迭代半径^[9]等算法，效率已能够满足要求. 考虑到变量 ${X^2} = \displaystyle\sum\limits_{i = 1}^n {X_i^2}$ ～χ²(n)，各球层（等概率划分为K层）的半径下限η_k为

${\eta _k} = \sqrt {{F^{ - 1}}\left(F(\beta _{{\text{form}}}^2) + \frac{{1 - F(\beta _{{\text{form}}}^2)}}{k} (k - 1)\right)},$

(20)

式中：F⁻¹(•)为χ²(n)累积分布的逆函数.

根据上述分析结果，结构失效概率的无偏估计可以通过球层抽样来实现. 在确定R > β_form的外部正态空间球层划分的基础上，通过控制L的抽样概率区间，在各球层空间内分别进行结构随机抽样并计算相应球层的条件失效概率. 在此基础上，利用式（19）将各层抽样结果进行叠加，最终可得到结构失效概率的估计值. 与其他算法相比，所提出方法避开了结构可靠概率最高的R ≤ β_form的球形区域（这些区域的概率可知且无效），并将抽样重点分层转移到结构失效域邻近的区域. 这有效改善了小失效概率事件的抽样效率低下的问题. 同时，由于抽样的无方向侧重性和均匀性，算法极大地提升了对于多设计验算点和多重要失效域问题的适用性.

2.2 算法步骤

在获取结构一阶可靠度指标β_form的基础上，对抽样向量Z分析构造，并结合标准正态空间球层划分以及逐层抽样策略，最终可得到结构失效概率的估计值. 基于上述思路，所构造抽样方法的主要步骤总结如下：

步骤1　根据所分析问题，确定结构基本向量X及其所对应标准正态空间维度n.

步骤2　对R > β_form外的标准正态空间，以原点为中心进行球层划分，可以按照等概率原则或其他合理原则划分为k层；当采用等概率原则时，各球层半径下限η_k按式（20）确定.

步骤3　在第k层空间内（半径为[η_k，η_k+1]）对变量L和X分别进行抽样，并生成样本点l和x；其中，L服从截尾的χ²(n)分布，截尾上、下限分别为 $\eta_k^2$ 和 $\eta_{k+1}^2$ ，X服从n维独立的标准正态分布；对样本l和x进行组合，获取抽样点z = $\sqrt l$ x/||x||，并判断此时g(z)是否小于0.

步骤4　重复抽样N_k次，检验每次样本点z是否落入失效域F^*（g(z)<0），并据此估计空间第k层球层内结构的条件失效概率.

步骤5　k=k + 1，重复步骤3～4，直至对所有球层均完成抽样，按式（19）将各层条件失效概率进行叠加组合，计算结构整体失效概率的估计值，抽样计算完毕.

在获取结构失效概率后，即可对其失效行为特点进行分析，评估结构的可靠性.

3. 算法验证及算例分析

为验证所提出算法的计算精度、效率以及适用性，分别采用所提出算法对3个典型算例进行失效概率抽样估计. 在此基础上，将计算结果与其他典型算法进行分析对比，以检验算法针对复杂结构可靠度问题的解决能力.

3.1 算例1

某一结构的功能函数^[16-17]为g(Z) = 5−Z₂−0.5 × (Z₁−0.1)²，计算结构的失效概率，式中，Z₁和Z₂均服从标准正态分布N(0, 1)，且二者之间相互独立.

首先，利用Hasofer Lind- Rackwitz Fiessler （HL-RF）算法^[18]对结构进行分析，经试算发现，结构共有2个设计验算点，分别为 ${\boldsymbol{z}}_1^*$ （−2.7409, 0.9648）和 ${\boldsymbol{z}}_2^*$ （2.9159, 1.0355），算例为多设计验算点问题，其极限状态曲面如所示. 图中， ${{\beta }_{\text{form}1}}$ 、 ${{\beta }_{\text{form}2}}$ 分别为 ${\boldsymbol{z}}_1^*$ 、 ${\boldsymbol{z}}_2^*$ 距原点的距离.

图 8 算例1极限状态曲面及球层划分

Figure 8. Limit state surface and spherical layer division in standard normal space for case 1

下载: 全尺寸图片幻灯片

由图8可以看出，β_form1= 2.9057，β_form2=3.0943. 根据结构一阶可靠度指标定义可得β_form= 2.9057，其所对应设计验算点为（−2.7409, 0.9648）. 利用上述信息，分别采用不同抽样方法对结构失效概率进行计算，其结果如表1所示. 表中，MCS为Monte Carlo仿真.

表 1 算例1不同方法计算结果对比

Table 1. Comparison of calculation results with different sampling methods for case 1

计算方法	抽样数/次	P_f/×10⁻³	β_form	相对误差/%	变异系数
HL-RF		1.832	2.9057	38.3
重要抽样	1000	2.090	2.8642	29.6	0.360
线抽样	1000	1.957	2.8850	34.1	0.004
球层抽样	1000	2.998	2.7480	0.9	0.062
MCS	1×10⁶	2.970	2.7511
注：各算法变异系数结果为重复100次计算得到，后续算例同.

下载: 导出CSV

| 显示表格

对比表1计算结果可以发现，由于算例具有2个设计验算点且可靠域为双侧有界区域，因此，其极限状态曲面极为复杂，采用传统线抽样及重要抽样方法均会产生过大误差，从而导致计算失败，即使增加抽样次数，其估计结果也难以有效改善. 这一现象主要是由于二者算法的抽样方向或区域的局部侧重性所导致的. 而所提出的球层抽样算法由于抽样原理的无方向性，因而可有效避免上述问题的产生，并具有较高的计算效率，在重复1000次后即可得到失效概率P_f的可靠估计值，其结果相对于Monte Carlo解，误差仅为0.9%，表明算法具有良好的适用性和计算精度.

进一步分析算法的抽样效率，可得P_f估计值变异系数随抽样次数的变化规律如图9所示. 由图可以看出，算法具有良好的收敛性，在抽样超过2000次后，其变异系数已不足4%，结构失效概率基本达到收敛. 这一结果表明，算法具有良好的抽样估计效率，在解决复杂可靠度问题时，其计算工作量将大为降低.

图 9 算例1球层抽样失效概率变异系数随抽样次数变化

Figure 9. Coefficient of variation in failure probability of spherical layer sampling with sampling times for case 1

下载: 全尺寸图片幻灯片

3.2 算例2

某一屋架结构^[19]如所示，屋架上弦杆和其他压杆采用钢筋混凝土构件，下弦杆和其他拉杆采用钢构件. 设屋架承受均布荷载q（kN/m）作用，经简化并由结构受力分析可得，屋顶点C竖直方向的挠度 ${\varDelta _{{C}}}{\text{ = }}\dfrac{{q{l_{\mathrm{s}}^2}}}{2}\left( {\dfrac{{3.81}}{{{E_{\text{c}}}{A_{\text{c}}}}} + \dfrac{{1.13}}{{{E_{\text{s}}}{A_{\text{s}}}}}} \right)$ ，其中，A_c、E_c、A_s、E_s和l_s分别为混凝土和钢杆件的截面积、弹性模量和长度. 以点C挠度不超过3 cm为约束条件，构造结构的功能函数g=0.03−Δ_C. 设结构中所有变量均服从正态分布且相互独立，其变量分布参数如表2所示，计算结构的失效概率.

图 10 算例2屋架结构示意

Figure 10. Illustration of case 2 truss structure

下载: 全尺寸图片幻灯片

表 2 算例2随机变量的分布参数

Table 2. Distribution parameters of random variables for case 2

项目	q/ （kN•m⁻¹）	l_s/m	A_s/m²	A_c/m²	E_s/ GPa	E_c/ GPa
均值	20	12	9.82 × 10^-4	0.04	100	20
变异系数	0.07	0.01	0.06	0.12	0.06	0.06

下载: 导出CSV

| 显示表格

算例2为一单榀桁架结构，针对这一结构的功能函数，首先采用HL-RF算法进行分析，可得结构一阶可靠度指标及相应设计验算点分别为β_form=2.4287及z^*（1.5063, 0.4735, −0.9533, −0.4334, −1.0744, −1.0744）. 基于上述结果，分别采用不同典型抽样方法对其失效概率进行估计分析，并对比相应结果，其结果如表3所示.

表 3 算例2不同方法计算结果对比

Table 3. Comparison of calculation results with different sampling methods for case 2

计算方法	抽样数/次	P_f/×10⁻³	β_form	相对误差/%	变异系数
HL-RF		7.577	2.4287	19.0
子集抽样^[19]	1500	9.317	2.3528	0.4	0.242
球层抽样	1500	9.396	2.3496	0.5	0.148
MCS^[19]	5 × 10⁷	9.352	2.3514

下载: 导出CSV

| 显示表格

对比表3计算结果发现，采用子集抽样算法与所提出球层抽样算法均可实现对算例2结构失效概率的准确估计，二者误差较小且十分接近. 但分析其变异系数可以看出，当抽样次数为1500次时，球层抽样算法P_f变异系数仅为0.148，与子集抽样算法相比减小39%. 这一结果表明，所提出球层抽样算法具有更加优良的抽样效率与收敛性，当抽样次数超过4500次后，其变异系数已不足0.090，结构失效概率收敛迅速，而采用子集抽样算法则需要更多的抽样次数方能收敛，其速度略慢. 由上述对比结果可知，所提出算法计算精度及抽样效率均较为优良.

3.3 算例3

某施工中钢斜拉桥结构如图11所示，其施工期最大单悬臂状态是一个潜在的控制工况. 此时结构主要承受恒荷载和施工临时荷载作用. 结构功能函数为g=R(k_m, k_f, k_q)−S_D−S_Q. 其中：R(k_m, k_f, k_q)为结构的抗力函数，k_m、k_f和k_q分别为其材料强度不定性变量、几何尺寸不定性变量以及计算模式不定性变量，S_D和S_Q分别为结构的恒载效应和施工临时荷载效应. 结构抗力和荷载变量服从正态分布，结构抗力分布参数如表4所示，表中，t为钢材厚度. 施工临时荷载变量分布参数如表5所示（结构自重比例系数均值为1.0237、变异系数为0.0400^[17]），计算结构的失效概率.

图 11 算例3钢斜拉桥最大单悬臂施工工况示意（单位：m）

Figure 11. Illustration of case 3 steel cable-stayed bridge in longest single cantilever construction phase （unit: m）

下载: 全尺寸图片幻灯片

表 4 算例3抗力不确定性变量分布参数

Table 4. Distribution parameters of resistance random variables for case 3

项目	k_m				k_f	k_p
项目	t≤16 mm	16 mm<t≤25 mm	25 mm<t≤38 mm	38 mm<t≤50 mm	k_f	k_p
偏差系数	1.1511	1.1228	1.0991	1.2074	1.0000	1.0890
变异系数	0.0656	0.0670	0.0843	0.0917	0.0200	0.1080

下载: 导出CSV

| 显示表格

表 5 算例3荷载变量分布参数

Table 5. Distribution parameters of temporary load random variables during construction for case 3

项目	桥面吊机和焊接设备	焊接平台	张拉平台
均值/kN	2402	150	399
变异系数	0.0500	0.0500	0.0500

下载: 导出CSV

| 显示表格

根据结构描述，算例3为一悬臂施工钢斜拉桥结构，其构造组成十分复杂，因而可有效检验所提出算法在实际结构中的适用性. 针对这一结构，首先计算其一阶可靠度指标，可得β_form=6.9533. 利用上述信息，采用所提出方法对结构失效概率进行求解计算，可得其结果如表6所示.

表 6 算例3不同方法计算结果对比

Table 6. Comparison of calculation results with different sampling methods for case 3

计算方法	抽样数/ 次	P_f/× 10⁻¹²	β_form	相对误差%
HL-RF		1.784	6.9533	43.1
球层抽样	1×10⁴	3.048	6.8774	2.8
改进 MCS^[17]	1×10⁵	3.137	6.8733

下载: 导出CSV

| 显示表格

由表6计算结果可以看出，算例3结构失效概率相对前2个算例更小，应用传统Monte Carlo法将极大增加抽样工作量. 而采用所提出的球层抽样方法可有效降低抽样计算量，其结果吻合较好，误差满足工程精度要求，算法具有良好的抽样效率与估计性能. 由这一结果，在针对复杂结构失效问题分析时，算法同样具有良好的适用性.

综合上述分析结果可以看出，算法在所有算例中均展现出良好的计算精度与性能，其正确性得到有效验证. 进一步地总结分析各算法特点，可以发现：

1）球层抽样算法针对结构多设计验算点等复杂体系可靠度问题，具有更好的求解性能与适用性. 除β_form数值外，算法无需其他额外信息. 而线抽样等估计方法除需对算法进行改进与修正^[17]外，尚需较多额外设计验算点信息，这在实际很多情况下（判断并识别出所有设计验算点信息）较难实现，因而其求解面临较多困难. 所提出球层抽样算法具有更广泛的适用性与优良求解性能.

2）与重要抽样、子集抽样等目前流行算法相比，所提出球层抽样算法具有更小的P_f估计变异系数与波动性，其抽样方差缩减能力相对更加优良，从而可有效实现结构失效概率的快速稳健收敛. 在解决复杂结构、小失效概率问题时，算法具有优势.

由上述对比分析，所提出球层抽样算法具有良好的抽样效率与问题适用性，其应用前景良好广阔.

4. 结　语

1）针对标准正态空间随机向量，基于距离与角度参量分离思想，对其进行了重构说明，并验证所构造向量的标准正态性与相互独立性，为结构随机样本的灵活抽取与方法建立提供了理论基础.

2）基于所构造向量并结合半径R>β_form外正态空间球形分层抽样策略，提出一种计算结构失效概率的球层抽样算法，该算法具有良好的无方向侧重性，可有效改善小失效概率事件抽样效率低下的问题，其性能较为优良.

3）相对其他传统抽样方法，所提出球层抽样算法具有更好的抽样效率与适用性. 算例对比结果表明：相同抽样次数条件下，球层抽样算法P_f估计变异系数较小，且可有效解决结构多设计验算点等复杂可靠度问题，更适用于小失效概率以及体系可靠度等问题的求解.

4）由算法解决复杂结构可靠度问题的独特优势所决定，结合实际工程结构，开展大跨度钢斜拉桥等复杂结构体系可靠度规律的探讨，是下一步研究的重点.

图 1 标准正态空间中样本点

Figure 1. Sampling point in standard normal space

下载: 全尺寸图片幻灯片

图 2 Z_i抽样频率密度与标准正态分布概率密度对比

Figure 2. Fig. 2 Comparison of probability density between Z_i sampling and standard normal distribution

下载: 全尺寸图片幻灯片

图 3 n维标准正态空间中n−1维球面均匀分布

Figure 3. Spherical uniform distribution of n−1 dimensions in standard normal space

下载: 全尺寸图片幻灯片

图 4 不同维度下Z_i前四阶矩抽样结果

Figure 4. Sampling results of first four-order moments of Z_i in different dimensions

下载: 全尺寸图片幻灯片

图 5 可靠域中多维球形空间

Figure 5. Mlti-dimensional spherical space in reliable domain

下载: 全尺寸图片幻灯片

图 6 标准正态空间按球层划分

Figure 6. Spherical layer division in standard normal space

下载: 全尺寸图片幻灯片

图 7 各球层内随机样本点生成

Figure 7. Generation of random sample points in each spherical layer

下载: 全尺寸图片幻灯片

图 8 算例1极限状态曲面及球层划分

Figure 8. Limit state surface and spherical layer division in standard normal space for case 1

下载: 全尺寸图片幻灯片

图 9 算例1球层抽样失效概率变异系数随抽样次数变化

Figure 9. Coefficient of variation in failure probability of spherical layer sampling with sampling times for case 1

下载: 全尺寸图片幻灯片

图 10 算例2屋架结构示意

Figure 10. Illustration of case 2 truss structure

下载: 全尺寸图片幻灯片

图 11 算例3钢斜拉桥最大单悬臂施工工况示意（单位：m）

Figure 11. Illustration of case 3 steel cable-stayed bridge in longest single cantilever construction phase （unit: m）

下载: 全尺寸图片幻灯片

表 1 算例1不同方法计算结果对比

Table 1. Comparison of calculation results with different sampling methods for case 1

计算方法	抽样数/次	P_f/×10⁻³	β_form	相对误差/%	变异系数
HL-RF		1.832	2.9057	38.3
重要抽样	1000	2.090	2.8642	29.6	0.360
线抽样	1000	1.957	2.8850	34.1	0.004
球层抽样	1000	2.998	2.7480	0.9	0.062
MCS	1×10⁶	2.970	2.7511
注：各算法变异系数结果为重复100次计算得到，后续算例同.

下载: 导出CSV

表 2 算例2随机变量的分布参数

Table 2. Distribution parameters of random variables for case 2

项目	q/ （kN•m⁻¹）	l_s/m	A_s/m²	A_c/m²	E_s/ GPa	E_c/ GPa
均值	20	12	9.82 × 10^-4	0.04	100	20
变异系数	0.07	0.01	0.06	0.12	0.06	0.06

下载: 导出CSV

表 3 算例2不同方法计算结果对比

Table 3. Comparison of calculation results with different sampling methods for case 2

计算方法	抽样数/次	P_f/×10⁻³	β_form	相对误差/%	变异系数
HL-RF		7.577	2.4287	19.0
子集抽样^[19]	1500	9.317	2.3528	0.4	0.242
球层抽样	1500	9.396	2.3496	0.5	0.148
MCS^[19]	5 × 10⁷	9.352	2.3514

下载: 导出CSV

表 4 算例3抗力不确定性变量分布参数

Table 4. Distribution parameters of resistance random variables for case 3

项目	k_m				k_f	k_p
项目	t≤16 mm	16 mm<t≤25 mm	25 mm<t≤38 mm	38 mm<t≤50 mm	k_f	k_p
偏差系数	1.1511	1.1228	1.0991	1.2074	1.0000	1.0890
变异系数	0.0656	0.0670	0.0843	0.0917	0.0200	0.1080

下载: 导出CSV

表 5 算例3荷载变量分布参数

Table 5. Distribution parameters of temporary load random variables during construction for case 3

项目	桥面吊机和焊接设备	焊接平台	张拉平台
均值/kN	2402	150	399
变异系数	0.0500	0.0500	0.0500

下载: 导出CSV

表 6 算例3不同方法计算结果对比

Table 6. Comparison of calculation results with different sampling methods for case 3

计算方法	抽样数/ 次	P_f/× 10⁻¹²	β_form	相对误差%
HL-RF		1.784	6.9533	43.1
球层抽样	1×10⁴	3.048	6.8774	2.8
改进 MCS^[17]	1×10⁵	3.137	6.8733

下载: 导出CSV

[1]	吕震宙,宋述芳,李洪双,等. 结构机构可靠性及可靠性灵敏度分析[M]. 北京: 科学出版社,2009: 13-85,194-242.
[2]	MELCHERS R E. Importance sampling in structural systems[J]. Structural Safety, 1989, 6(1): 3-10. doi: 10.1016/0167-4730(89)90003-9
[3]	GEYER S, PAPAIOANNOU I, STRAUB D. Cross entropy-based importance sampling using Gaussian densities revisited[J]. Structural Safety, 2019, 76: 15-27. doi: 10.1016/j.strusafe.2018.07.001
[4]	KOUTSOURELAKIS P S, PRADLWARTER H J, SCHUËLLER G I. Reliability of structures in high dimensions, part Ⅰ: algorithms and applications[J]. Probabilistic Engineering Mechanics, 2004, 19(4): 409-417. doi: 10.1016/j.probengmech.2004.05.001
[5]	SCHUËLLER G I, PRADLWARTER H J, KOUTSOURELAKIS P S. A critical appraisal of reliability estimation procedures for high dimensions[J]. Probabilistic Engineering Mechanics, 2004, 19(4): 463-474. doi: 10.1016/j.probengmech.2004.05.004
[6]	吕召燕,吕震宙,张磊刚,等. 基于条件期望的改进线抽样方法及其应用[J]. 工程力学,2014,31(4): 34-39. LÜ Zhaoyan, LÜ Zhenzhou, ZHANG Leigang, et al. An improved line sampling method and its application based on conditional expectation[J]. Engineering Mechanics, 2014, 31(4): 34-39.
[7]	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. doi: 10.1016/S0266-8920(01)00019-4
[8]	AU S K, WANG Y. Engineering risk assessment with subset simulation[M]. Singapore: Wiley, 2014: 157-204.
[9]	KATAFYGIOTIS L S, CHEUNG S H. Application of spherical subset simulation method and auxiliary domain method on a benchmark reliability study[J]. Structural Safety, 2007, 29(3): 194-207. doi: 10.1016/j.strusafe.2006.07.003
[10]	KATAFYGIOTIS L, CHEUNG S H, YUEN K V. Spherical subset simulation (S³) for solving non-linear dynamical reliability problems[J]. International Journal of Reliability and Safety, 2010, 4(2/3): 122-138. doi: 10.1504/IJRS.2010.032442
[11]	NOWAK A S, COLLINS K R. Reliability of structures [M]. 2nd edition. New York: CRC Press, 2013: 97-100,126-161.
[12]	何平. 数理统计与多元统计[M]. 成都: 西南交通大学出版社,2004: 18-27.
[13]	盛骤,谢式千,潘承毅. 概率论与数理统计[M]. 4版. 北京: 高等教育出版社,2008: 106-110,135-208.
[14]	巩祥瑞,吕震宙,孙天宇,等. 一种新的矩独立重要性测度分析方法及高效算法[J]. 北京航空航天大学学报,2019,45(2): 283-290. GONG Xiangrui, LÜ Zhenzhou, SUN Tianyu, et al. A new moment-independent importance measure analysis method and its efficient algorithm[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(2): 283-290.
[15]	GONG W P, JUANG C H, MARTIN J R II, et al. New sampling method and procedures for estimating failure probability[J]. Journal of Engineering Mechanics, 2016, 142(4): 04015107.1-04015107.10.
[16]	DER KIUREGHIAN A, DAKESSIAN T. Multiple design points in first and second-order reliability[J]. Structural Safety, 1998, 20(1): 37-49. doi: 10.1016/S0167-4730(97)00026-X
[17]	白冰. 大跨度钢斜拉桥施工及运营过程系统可靠度研究[D]. 成都: 西南交通大学,2015.
[18]	CREMONA C. Structural performance: probability-based assessement[M]. London: ISTE,2011.
[19]	薛国峰. 结构可靠性和概率失效分析数值模拟方法[D]. 哈尔滨: 哈尔滨工业大学,2010.