动力学系统的结构可靠性问题是海上平台、桥梁工程等领域的研究热点,受到诸如地震、风及海浪等随机激励的动力学系统的首穿失效问题是其研究的焦点之一．首穿失效概率作为首穿失效问题的重要指标, 自Rice[1-2]1945年的开拓性研究以来,迄今为止,即使是线性结构动力学系统,也无可利用的精确解,因此，数值方法、模拟技巧成为解决此问题的关键．研究者们先后提出路径积分法[3]、胞映射法[4]、随机平均法[5-6]等数值方法，然而这些方法的复杂度会随着状态空间维数的增加呈指数增长趋势．近年来,研究者们又提出了基于方差缩减技术的Monte-Carlo模拟方法(MCS),此后，该法迅速成为计算首穿失效概率的高效方法[7-9]．在诸多方差缩减技术中,重要抽样法[10]被广泛使用,针对线性系统，2001年Au等深入分析线性动力学系统失效区域,利用设计点构建重要抽样密度函数,提出基于重要抽样的估计方法[11]．2005年Ka-Veng等提出概率简单加方法[12]；次年Lambros等提出区域分解法[13]；2009年何军提出了非平稳Gauss白噪声激励下线性系统条件首次穿越概率的近似解析解[14]；2011年Zuev等利用复数值的Ornstein-Uhlenbeck过程，构造符合已给噪声功率谱的稳态噪声，基于Gisranov定理，提出horseracing模拟算法[15]．2014年Valdebenito等提出不确定线性系统的首穿失效概率的估计方法[16]．2016年Wang等[17]基于重要抽样法，提出风激励下线性系统的首穿失效概率的估计方法．对非线性动力学系统,2004年Pradlwarter等利用线性抽样提出平均失效概率流方法[18]；2016年He(何军)等利用极值分布的尾近似，提出非线性结构系统的首穿失效概率的估计方法[19]．而借助等效线性化方法和线性系统失效概率的估计方法，也是解决非线性系统失效问题的方法之一，例如2007年Olsen等提出两步迭代算法[20]；2012年任丽梅等基于测度变换，提出首穿失效概率估计方法[21]．
本文基于重要样本法，提出了线性与非线性动力学系统首穿失效概率的估计方法．其关键是将系统失效域重构为一系列互斥的基本失效区域之和,并利用基本失效域概率构造重要抽样密度函数，使得首穿失效概率的估计转换为估计一系列落在不同基本失效区域的样本数之比．而对非线性系统,重点是寻求与之等效的线性系统．最后给出线性与非线性数值算例，并将本文方法与MCS及区域分解法比较,结果显示本文方法是正确有效的．
考虑Gauss白噪声激励下线性结构动力学系统：
(1)
其中β是阻尼系数,ω0是自然频率，W(t)是单边谱密度为G0的Gauss白噪声,Y(t)为系统响应，激励持时为T．令b是安全域边界,则
是失效域．系统的首穿失效概率为
PF=P(F)=P{∃t∈(0,T),|Y(t)|>b}．
(2)
根据Duhamel积分,系统(1)的激励与响应关系可表示为Y(t)=
h(t,τ)W(τ)dτ，其中h(t,τ)是单位脉冲响应函数,表示τ时刻输入的单位脉冲在t时刻的输出．令Δt=T/n，tk=kΔt,k=1,2,…,n，且当Δt→0时，g(k,s)→h(tk,ts)．令Zk(1),Zk(2),…,Zk(s)表示相互独立的正态随机变量，则离散化后系统响应Y(t)为
(3)
相应的失效域也分解为
此时称Fk(k=1,2，…,n)为基本失效域．则tk时刻的极限状态曲面方程为Y(k)=b，tk时刻的可靠性指标βk为
(4)
解tk时刻的优化问题
可得tk时刻的局部设计点：
(5)
其中
且
当k=n时，称
为局部设计点；当k≠n时，称为全局设计点．基本失效域可被其设计点完全刻画[22-24]．则tk时刻基本失效域的概率为
(6)
其中Φ(·)为标准正态分布函数．
由于基本失效域Fk相互重叠，令
(7)
其中“-”表示集合减法．易知Gi是Fi的子集，Gi由属于Fi但不属于Fj(j=1,2,…,i-1)的样本点构成．则失效域F可重构为
即首穿失效概率可为
即
(8)
假如可以有效地估计ηk，则首穿失效概率的估计为
(9)
实际上
其中Nk表示落在基本失效域Fk中的样本数量,
表示落在Gk中的样本数量．则关键为怎样生成基本失效域Fk中的样本(具体见第2节中的step 4)及怎样缩减样本方差．现解决缩减样本方差的问题．由于
(10)
构建重要抽样密度函数：
(11)
则首穿失效概率的估计可为
(12)
其中M表示模拟样本的数量，
显然
是均值为pkj的(0-1)随机变量．估计的变异系数为
(13)
计算时pkj用估计值
代替，
表示同时落在Gkj和Fkj中的样本数目．
Step 1 得到形如式(3)的动力学系统离散表示．
Step 2 根据式(5)计算设计点
将
代入式(4)计算可靠性指标βk(k=1,2，…,n),根据式(6)计算基本失效域的概率P(Fk)．
Step 3 根据式(11)选择一个失效域Fk：产生[0,1]上均匀分布随机数U1，并计算区间
视U1落入哪个区间，则选k=i．
Step 4 产生Fk中的样本Zs：产生n维标准正态随机向量Z及[0,1]上的均匀分布随机数U2，计算
Step 5 根据式(3)及Zs计算Y(k)(k=1,2，…，n)．若Y(k)满足Y(k)>b，Y(k-1)<b，…，Y(1)<b，则Y(k)落在Gk内,此时
否则
Step 6 根据式(12)计算首穿失效概率,式(13)计算变异系数．
Gauss白噪声激励的线性动力学系统:
(14)
其中
安全域边界分别取b=2.24 cm和b=2.8 cm，激励持时T=20 s．脉冲响应函数为
其中
取Δt=0.02，则n=1 000．对每个安全域边界进行10次试验，每次100个样本．模拟结果见表1、2，表中第1、3列是每次试验同时落入Fk与Gk的样本数量；第2、4列是每次试验的首穿失效概率估计值．将10次试验结果的算术平均值作为首穿失效概率的估计值．当b=2.24 cm时,估计值为3.025×10-3，变异系数为25.11%；当b=2.8 cm时,估计值为3.678×10-5，变异系数为23.90%．
表1 b=2.24 cm，M=100时的试验结果
Table 1 Results withb=2.24 cm andM=100
mPFmPF113.422 855 4E-361.867 012 0E-3103.111 686 8E-3123.734 024 2E-392.800 518 1E-3113.422 855 4E-3103.111 686 8E-382.489 349 3E-3134.045 192 8E-3113.422 855 4E-3
图1是表1中第一次试验的首穿失效概率随样本的变化情况，此时100个样本中有11个样本同时属于Fk与Gk，每增加一个这样的样本, 首穿失效概率估计值就上一台阶，估计值会相差
随着这样样本数的增多,台阶的高度则越来越矮．图2给出了本文方法所得估计值随样本数(M=500)的变化情况，直线是MCS方法(106个样本)得到两个安全域边界下首穿失效概率的估计值：2.961×10-3和3.701×10-5．图2显示随样本数的增多,利用本文方法(虚线)所得概率估计值快速逼近2.961×10-3和3.701×10-5，大概200个样本之后,首穿概率估计值与MCS方法模拟值几乎相差无几．
表2 b=2.8 cm，M=100时的试验结果
Table 2 Results withb=2.8 cm andM=100
mPFmPF143.456 189 4E-5143.456 189 4E-5194.690 542 8E-5184.443 672 2E-5122.962 448 2E-5163.949 930 8E-5184.443 672 2E-5122.962 448 2E-5102.468 706 8E-5163.949 930 8E-5
虽然线性动力学系统的首穿失效概率无解析解，但目前研究者已提出几个高效的模拟算法，则对非线性动力学系统，寻求它的等效线性系统，不失为解决问题的方法．本文采用的线性化原理是非线性与等效线性系统对指定的安全域边界具有相同的平均上穿率[25-27]．
考虑非线性动力学系统：
(15)
其中β=0.05和ω0=1.0表示当ε=0时的线性振子的阻尼与自然频率，ε(ε≥0)是非线性参数，γ=0.3，激励持时T=15 s．系统(15)的失效域
和
的联合概率密度函数为
(16)
其中
分别表示当ε=0时X(t)与的稳态方差,A是归一化常数，
K1/4是修正的Bessel函数．
系统(15)的平均上穿率为
(17)
令等效线性系统为
(18)
ωe是等效线性化参数．则系统(18)的平均上穿率为
(19)
其中
则当
(20)
时可得非线性系统的等效线性化参数ωe．根据式(20)绘制平均上穿率随非线性参数的变化图，图3显示随安全域边界的增大，平均上穿率逐渐减小，对同一个安全域边界，非线性参数越大平均上穿率越小．本文考虑3组不同安全域边界、不同非线性参数情形：
1) 安全域边界b=3σ0，非线性参数ε=0.1，由式(20)可得等效线性参数ωe=1.310 634；
2) 安全域边界b=2.5σ0，非线性参数ε=0.5，等效线性参数ωe=1.858 361；
3) 安全域边界b=1.5σ0，非线性参数ε=2，等效线性参数ωe=2.164 058．
则非线性系统(15)的首穿失效概率近似为
PF(T)=P{∃t∈(0,T),|X(t)|>b}≈P{∃t∈(0,T),|Y(t)|>b}，
(21)
等效线性化系统的脉冲响应函数为
其中
取Δt=0.05，n=300．
图3不同非线性化参数及不同安全域边界下的平均上穿率
Fig. 3 Mean high level crossing rates for different threshold levelsb and non-linearity parametersε
利用本文方法每组参数进行10次试验，考虑到3.1小节的试验效果，每次试验的样本数取M=200．3种参数下的试验结果见表3～5，表中第1、3列仍是每次试验200个样本中同时落入Fk与Gk的样本数量；第2、4列分别是每次试验的首穿失效概率估计值．10次试验的首穿失效概率算术平均值及变异系数见表6．为了比较，利用MCS(106(或108)个样本)方法，对非线性系统进行Euler(欧拉)离散，可得3组参数下首穿失效概率估计值，见表6；对等效线性系统，本文还利用区域分解法对3组参数情况给出了首穿失效概率的估计值，每组参数试验10次，每次100个样本，10次试验的首穿失效概率算术平均值及变异系数见表6．图4是本文方法(虚线)与区域分解法(实线)给出的概率估计值随样本数的变化情况．
表3 第一组参数下的试验结果
Table 3 Results of case 1
mPFmPF203.887 448 5E-5265.053 683 0E-5193.693 076 0E-5234.470 565 5E-5214.081 821 0E-5183.498 703 5E-5316.025 545 0E-5244.664 938 0E-5173.304 331 0E-5193.693 076 0E-5
表4 第二组参数下的试验结果
Table 4 Results of case 2
mPFmPF307.772 795 2E-7235.959 142 9E-7318.031 888 1E-7287.254 608 9E-7359.068 260 7E-7318.031 888 1E-7297.513 701 8E-7266.736 422 6E-7391.010 463 4E-6379.586 447 0E-7
表5 第三组参数下的试验结果
Table 5 Results of case 3
mPFmPF302.124 923 7E-3251.770 769 8E-3322.266 585 4E-3211.487 446 6E-3251.770 769 8E-3221.558 277 4E-3241.699 939 0E-3261.841 600 6E-3261.841 600 6E-3191.345 785 1E-3
表6 本文方法与MCS方法、区域分解法的结果比较
Table 6 Comparison of the proposed method, the MCS and the domain decomposition method
method and parametercase123the proposed methodestimated value4.431 7×10-58.005 9×10-71.770 8×10-3coefficient variation14.02%15.75%16.25%domain decomposition methodestimated value4.487 1×10-59.702 7×10-71.826 1×10-3coefficient variation10.67%12.25%11.73%MCSestimated value4.209 4×10-58.437 3×10-72.011 9×10-3
由表6及图4可见，与MCS相比，本文方法是正确有效的；与区域分解法相比，表6显示本文方法正确，但区域分解法的变异系数更小，图4也显示了区域分解法给出的估计值在大概50个样本后便趋于稳定，显示出区域分解法的稳健性，而本文方法至少需200个样本，但本文方法思想简单、算法简单、重要抽样密度函数构造简单，也是有效的估计方法．
图43组参数下，本文方法(虚线)与区域分解法(实线)的首穿失效概率随样本数的变化图
Fig. 4 Variation of the first passage failure probability with the sample sizeM
(the proposed method(dot line), the domain decomposition method(solid line))
系统失效域是基本失效区域之和，而由局部设计点描述的基本失效区域之间存在着大量重叠．本文重构失效域，并基于重要抽样法、等效线性化及平均上穿率，提出了线性与非线性系统首穿失效概率的估计方法．关键点是失效域的重构和等效线性化原理．本文给出线性和非线性数值算例，并与MCS及区域分解方法相比较，显示出本文方法的正确性及有效性，虽不及区域分解法稳健高效，但由于思想简单、重要抽样密度函数构造简单，也不失为首穿失效概率估计的有效方法，并且算例结果也显示了本文等效线性化方法的有效性；另一方面，考虑到变异系数及工程应用的要求，利用本方法至少需200个样本，估计效果更好．
[1] RICE O C. Mathematical analysis of random noise[J].Bell System Technology, 1944,23(3): 282-332.
[2] RICE O C. Mathematical analysis of random noise[J].Bell System Technology, 1945,24(3): 46-156.
[3] ROBERTS J B. First passage probability for nonlinear oscillators[J].Journal of Engineering Mechanics, 1976,102: 851-866.
[4] SPENCER B F, BERGMAN L A. On the numerical solution of the Fokker-Planck equation for nonlinear stochastic systems[J].Nonlinear Dynamics, 1993,4(4): 357-372.
[5] WEI L, WEI X. Stochastic optimal control of first-passage failure for coupled Duffing-van der Pol system under Gaussian white noise excitations[J].Chaos,Solitons &Fractal, 2005,25(5): 1221-1228.
[6] WEI L, WEI X. First passage problem for strong nonlinear stochastic dynamical system[J].Chaos,Solitons &Fractal, 2006,28(2): 414-421.
[7] PRADLWARTERH J, SCHUELLER G I. Assessment of low probability events of dynamical systems by controlled Monte Carlo simulation[J].Probability Engineering Mechanics, 2009,14(3): 213-227.
[8] PROPPE C, PRADLWARTER H J, SCHUELLER G. Equivalent linearization and Monte Carlo simulation in stochastic dynamics[J].Probability Engineering Mechanics, 2003,18(1): 1-15.
[9] JUU O, HIROAKI T. Importance sampling for stochastic systems under stationary noise having a specified power spectrum[J].Probability Engineering Mechanics, 2009,24(4): 537-544.
[10] NEWTON N J. Variance reduction for simulate diffusions[J].SIAM Journal on Applied Mathematics, 1994,54(6): 1780-1805.
[11] AU K L, BECK J L. First excursion probability for linear systems by very efficient importancesampling[J].Probability Engineering Mechanics, 2001,16(3): 193-207.
[12] KA-VENG Y, LAMBROS S K. An efficient simulation method for reliability analysis of linear dynamical systems using simple additive rules of probability[J].Probability Engineering Mechanics, 2005,20(1): 109-114.
[13] LAMBROS K, SAIHUNG H. Domain decomposition method for calculating the failure probability of linear dynamic systems subjected to Gaussian stochastic loads[J].Journal of Engineering Mechanics, 2006,20: 475-486.
[14] 何军. 非平稳随机激励下系统首次穿越概率的近似解法[J]. 应用数学和力学, 2009,30(2): 245-252.(HE Jun. Approximation for the first passage probability of systems under non-stationary random excitation[J].Applied Mathematics and Mechanics, 2009,30(2): 245-252.(in Chinese))
[15] ZUEV K M, KATAFYGIOTIS L S. The horseracing simulation algorithm for evaluation of small failure probabilities[J].Probability Engineering Mechanics, 2011,26(2): 157-164.
[16] VALDEBENITO M A, JENSEN H A, LABARCA A A. Estimation of first excursion probabilities for uncertain stochastic linear systems subject to Gaussian load[J].Computers and Structures, 2014,138: 36-48.
[17] WANG J, KATAFYGIOTIS L S, FENG Z Q. An efficient simulation method for the first excursion problem of linear structures subjected to stochastic wind loads[J].Computers &Structures, 2016,166: 75-84.
[18] PRADLWARTER H J, SCHUELLER G J. Excursion probabilities of non-linear systems[J].International Journal of Non-Linear Mechanics, 2004,39(9): 1447-1452.
[19] HE J, GONG J H. Estimate of small first passage probabilities of nonlinear random vibration systems by using tail approximation of extreme distributions[J].Structural Safety, 2016,60: 28-36.
[20] OLSEN I A, NAESS A. An importance sampling procedure for estimating failure probabilities of non-linear dynamic systems subjected to random noise[J].International Journal of Non-Linear Mechanics, 2007,42(6): 848-853.
[21] 任丽梅, 徐伟, 肖玉柱, 等. 基于重要样本法的结构力学系统的首次穿越[J]. 力学学报, 2012,44(3): 648-652.(REN Limei, XU Wei, XIAO Yuzhu, et al. First excursion probabilities of dynamical systems by importance sampling[J].Chinese Journal of Theoretical and Applied Mechanics, 2012,44(3): 648-652.(in Chinese))
[22] KOO H, KIUREGHIAN A D, FUJIMURA K. Design-point excitation for non-linear random vibrations[J].Probabilistic Engineering Mechanics, 2005,20(2): 136-147.
[23] AU S K. Critical excitation of sdof elasto-plastic systems[J].Journal of Sound and Vibration, 2006,296(4/5): 714-733.
[24] VALDEBENITO M A, PRADLWARTER H J. The role of the design point for calculating failure probabilities in view of dimensionality and structural nonlinearities[J].Structural Safety, 2010,32(2): 101-111.
[25] TAKEWAKI I. Critical excitation for elastic-plastic structures via statistical equivalent linearization[J].Probability Engineering Mechanics, 2002,17(1): 73-84.
[26] SOONG T T, GRIGORIU M.Random Vibration of Mechancial and Structural Systems[M]. Englewood Cliffs, NJ: Prentice-Hall Inc,1997.
[27] OKSENDAL B.Stochastic Differential Equations:an Introduction With Application[M]. Berlin: Springer, 1998.
应用数学和力学编委会，ISSN 1000-0887
任丽梅. 基于失效域重构和重要抽样法的结构动力学系统首穿失效概率[J]. 应用数学和力学, 2019,40(4): 463-472.
REN Limei. The first passage failure probabilities of dynamical systems based on the failure domain reconstruction and important sampling method[J].Applied Mathematics and Mechanics, 2019,40(4): 463-472.