重要抽样方法

重要抽样方法
直接蒙特卡罗方法用于估计结构的失效概率，通用性强，简单直观，但计算效率很低，为了达到一定的模拟精度，需要抽取大量的样本点。

为提高计算效率，一系列的方差缩减技术得到应用，从而产生了一系列的数值模拟算法。

重要抽样方法是其中最基本也是最重要的一种，能够大幅度的缩减模拟所需抽取的样本点数量，广泛应用于各种结构的可靠性分析之中。

重要抽样方法的基本概念
造成直接蒙特卡罗方法估计失效概率效率低的主要原因，是按原概率密度函数生成的样本点，落入失效域的比例过低，也就是说，样本点中绝大部分对失效概率值的贡献是零。

如果能够提高样本点落入失效域的比例，从而提高样本点对失效概率的贡献，则能够提高模拟效率，降低样本点数量。

重要抽样方法的基本原理是在保持原有样本期望值不变的情况下，改变随机变量的抽样重心，改变了现有样本空间的概率分布，使其方差减小，这样，使对最后结果贡献大的抽样出现的概率增加，抽取的样本点有更多的机会落在感兴趣的区域，使抽样点更有效，以达到减小运算时间的目的。

我们假设一个关于随机变量x 的概率密度函数为()h x ，那么，()1x X
h x ∈=∑可
以得到以下的式子：
()()()()()()[()/]()x X x X
E g x g x f x g x f x h x h x ∈∈==∑∑
从上式中我们可以这样假设出一个新的函数()*g x ，并且，令
()()()()*g x g x f x h x =，
这个新的函数的期望值和原先的()g x 的期望值完全相同，但是它们的方差却是完全不相同的：
()()()()()()*2[()/]()x X
Var g x g x f x h x E g x h x ∈=-∑
可以很明显的看出，如果方差越小，那么计算的工作量就会越小，我们把()h x 例称作为重要抽样分布函数。

重要抽样的基本思路，那就是寻找一个新的分布函数，通过这个函数，我们能够达到保持期望值不变而方差减小的目的，从而提高运算的效率。

当()h x 取最优值的时候，会出现()()*0Var g x =的情况，而这种情况是最理
想的一种情况，在这样的情况下，只需要进行一次抽样就可以完成估计。

但是，这种情况毕竟只是一种理想的状态，因为()()Var g x 的真实值我们是无法知道的，所以在实际的计算中我们只能采用近似的最优分布来逼近()()Var g x ，从而来实现减小方差的目的。

而重要抽样分布函数算法的选择就成了制约计算时间的关键问题。

根据重要抽样方法的原理，重要抽样方法的具体描述如下：
将结构失效概率积分公式改写为
()()()[()]()[()]()[()][()]()()()f f x f x P I g x f x dx I g x h x dx E I g x E I g x R x h x h x +∞+∞-∞-∞⎛⎫==== ⎪⎝
⎭⎰⎰
()R x 称为重要抽样权函数：
()()()
f x R x h x =
()h x 必须满足如下条件： ()0,,()1h x x F h x dx +∞
-∞>∀∈=⎰ {}|()0x h x Ω=>称为抽样域，
易知F ∈Ω，按照重要抽样概率密度函数()h x 抽取的样本点12,,...,N x x x ，全部落在Ω内。

抽样域Ω和失效域F 越是接近，则按照概率密度函数()h x 抽取的样本点落入失效域的比例越大。

结构失效概率通过下式进行估计：
1111()
[()]()[()]()N N f f k k f x P P I g x R x I g x N N h x ==≈==∑∑
容易验证估计值f P 是无偏的，即
[]111[()]()[()]()N f f k E P E I g x R x N E I g x R x P N N
=⎡⎤⎡⎤==⋅⋅=⎢⎥⎣⎦⎣⎦∑ 估计的方差为
()()11[()]()f Var P Var I g x R x N N
==∆ 其中△为重要抽样估计的单位方差，与样本量N 无关，
()
()()()2222
22
[()]()[()]()[()]()[()]()[()]()()f f Var I g x R x E I g x R x E I g x R x E I g x R x P I g x R x h x dx P +∞
-∞∆=⎡⎤=-⎡⎤⎣
⎦⎣⎦=-=-⎰ 使()f Var P 取得最小值的ISD 函数被称为最优ISD 函数，易知最优ISD 函数为
下式所示的概率密度函数：
()()()()()
|opt f I g x f x h x f x F P ==
()opt h x 对应的估计方差为0。

尽管()opt h x 存在，但用()opt h x 作为ISD 函数的重要抽样方法无法应用，其原因有：（1）()opt h x 中包含失效概率f P ，其值为未知；（2）生成服从概率密度
函数为()opt h x 的分布的样本点，通常是一件很困难的事。

重要抽样方法的核心在于ISD 函数的构造，如果构造的ISD 函数不合适，可能导致对失效概率进行估计的方差非常大。

构造ISD 函数应当兼顾以下两个方面：
（1）增加有效抽样的比例；（2）增加有效抽样中对f P 值贡献大的抽样出现的
概率。

增加有效抽样的比例，即增加抽取的样本点落在失效域内的比例，可以通过使抽样域。

尽量接近失效域F 来达到；对失效概率f P 值贡献大的抽样点，即使
()|f x F 值较大的样本点。

此外，构造的ISD 函数还需要满足能够较容易地生成样本点的要求。

重要抽样法在构造抽样函数的时候不需要考虑极限状态曲面的形状。

由于计算简便，并不需要在计算之前做过多的准备工作，实际中的应用比较广泛。

常规的抽样方法关注了抽样点之间的相对位置，而没有考虑样本点的整体位置，这样在随机变量整个定义域内无中心点的抽样方式会使大部分样本点落入安全域内，失效样本过少，效率太低。

重要抽样法通过改变抽样重心，可使抽样区域整体移向失效域，使样本点落入失效域的概率大大增加，大大提高了抽样效率和结果的可信度。

简单抽样，单纯以参数的统计分布进行抽样，使得样本点大部分落于安全域内，抽样效率低，这种缺点在可靠性较高的产品上更为明显；重要抽样方法将抽样中心移至验算点P 点，使样本点落入失效域的概率大大增加，提
高抽样效率。

利用中心正态重要抽样方法，每次抽样中心选为上次迭代计算的验算点处，抽样密度函数为正态分布函数。

这样随着循环计算的进行，验算点越来越精确，抽样中心也会逐步移向真实验算点，使得抽样点以真实验算点为中心正态分布，构造相应的Kriging 模型在验算点附近的拟合精度也会很高，从而提高可靠度计算精度。

中心正态重要抽样法（MCCNIS ）
获得重要抽样密度函数的较简单方法是确定重要抽样域并把密度函数的中心放在这个区域内。

可以选择均匀分布的重要抽样密度函数并将中心放在设计点处，也可以选择将原随机向量的密度函数中心平移在极限状态面上，或是将具有与原分布相同的相关矩阵的n 维正态分布密度函数中心平移到设计点处，且相关系数不小于原随机变量之间相关系数。

ISD 函数的类型多种多样，在实际的研究工作中，为便于理论分析和简化计算，在一些情况下采用n 维正态分布概率密度函数作为ISD 函数是一个合理的选择。

正态分布密度函数的参数为采样中心μ和协方差矩阵S ，未知参数为n(n+3)/2个，当选择二维独立正态分布密度函数作为ISD 函数时，协方差矩阵S 为对角矩阵，未知参数减少为2n 个，且能够较容易的生成样本点。

采用近似算法中的H-L 算法求得设计点*x ，然后选取ISD 函数为以*x 为采样中心的n 维独立正态分布概率密度函数进行重要抽样模拟是合适的。

中心正态重要抽样法的密度函数为：
()*()h x x x φ=-
如果极限状态面非线性不是很高、接近平面，则这样抽样使得约有一半的抽样点落在失效域，抽样效率很好。

但是如果极限状态面曲率较大特别是如果极限状态面为一个以原点为中心的以β为半径的球面时，使用这种方法的效率就非常低了。

β-球面外的结尾正态重要抽样法(MCROIS)
重要抽样密度函数的定义域不一定是整个随机向量空间，失效概率是在失效域上的积分，因此重要抽样密度函数的定义域只需要包含失效域就可以了，最理想的重要抽样密度函数定义域等于失效域且重要抽样密度函数取值与随机向量
的概率密度函数取值成比例，当然通常这种重要抽样密度函数是难以构造的。

如果已经获得了设计点，由设计点的定义知道设计点到原点的距离是所有极限状态面上的点到原点距离中最小的，令这个最小距离为β，则
()()
()()()()()()
()()(
)()()()22222222
00000f P P g x P g x x P x P g x x P x P g x x P x P g x x P x ββββββββ=≤=≤≤≤+≤>>=≤≤≤+≤>> 又 ()()222221n i n i P x P x βββ=⎛⎫≤=≤=Γ ⎪⎝⎭
∑ ()()
2200P g x x β≤≤=
所以 ()()()()
22201f n P P g x x ββ=≤>-Γ 式中()n Γ为自由度为n 的2χ分布函数，x 为向量的模。

图1 β-球外重要抽样法
如图1示，圆心在原点以β为半径的球面(简称β一球面)内，全部处在安全域内，而在独立标准正态随机空间中随机向量取值在这部分区域的概率很高。

如果使用独立标准正态密度函数抽样，就会使得大量抽样点落在β一球内，为了避免这种情况出现，构造如下截尾正态重要抽样密度函数
()()0K x x h x x φββ
⎧>⎪=⎨>⎪⎩
式中，()()
211n K β=-Γ，重要抽样密度函数的抽样域如图1示，β一球面
外的截尾正态重要抽样可以使更多的抽样点落在失效域，从而使得抽样效率增加。

比较两种方法，MCCNIS 的抽样域是整个随机相向量空间，它是将重要抽样密度函数（即正态密度函数）的中心移动到设计点位置，从而抽样点以较高概率落在以设计点为中心的附近区域内，当极限状态面接近平面时有500%左右的抽样点落在失效域内；而MCROIS 的抽样域是中心在原点、半径为β的球面外的随机向量空间，抽样点以较高概率落在球外表面的附近区域内。

当极限状态面接近平面时MCROIS 的抽样点落到失效域的概率要低于MCCNIS ，但是当极限状态面接近球面时MCROIS 的抽样点落到失效域的概率要明显高于MCCNIS 。

基于Kriging 模型的重要抽样方法
作为最重要的数值模拟方法之一，重要抽样方法被广泛应用于结构的可靠性分析中。

相比直接蒙特卡罗方法，重要抽样方法在估计失效概率时所需要的样本点数量大幅度减小。

尽管如此，其所需样本点仍然较多，在确定采样中心之后，要使估计的变异系数低于0.07，所需的样本点数量大于2000。

在结构功能函数没有显式表达的情况下，尤其是进行一次结构分析需要较大计算代价的情况下，采用替代分析模型来代替实际的结构分析，可以进一步减少重要抽样所需结构分析的次数，从而降低模拟的计算代价。

基本思想
已有的基于Kriging 模型的重要抽样方法(Importance Sampling based onKriging model, ISK)，通过实验设计方法选取训练样本，建立Kriging 模型后，对随后所抽取的样本点，直接采用Kriging 模型的预测值代替实际的功能函数值来进行失效概率的估计。

其缺点在于，建立的Kriging 模型对随后抽取的样本点的功能函数进行预测，没有对其预测的准确度进行检验，如果选取训练样本不当，可能导致失效域边界附近的样本点的功能函数值出现较大的偏差，从而使得对失效概率的估计出现较大偏差。

采用学习函数，自动寻找并添加最适合的训练样本，不断修正Kriging 模型的参数值。

即使最初的训练样本选取不恰当，通过不断的修正，最终也能够得到一个可以对样本点进行精确分类的Kriging 模型。

算法描述
重要抽样方法的核心是构造重要抽样密度(ISD)函数，而构造ISD 函数需要获得失效域的信息。

在结构功能函数没有显式表达的情况下，获取失效域的信息只能通过对一些离散的样本点进行相应的结构分析来实现，相对来说比较困难。

Hurtado 提出了一种利用改进的Metropolis 算法，从而在结构功能函数没有显式表达的情况下获得近似设计点的方法：既然设计点是失效域内概率密度函数值最大的点，从失效域内的一个点出发，利用MMA 生成候选状态点，如果候选状态点的概率密度函数值比当前点的概率密度函数值大，则执行接受/拒绝阶段，否则从当前点出发重新生成候选状态点。

在获得近似设计点后，自然的做法是将ISD 函数选取为以近似设计点为采样中心的无关正态概率密度函数。

Hurtado 所提方法对只有单设计点的情况适用，当存在多个设计点且事先并不知道这一点时，该方法可能只能求解其中某一个设计点，从而由一个设计点构造的ISD 函数用于重要抽样，对失效概率的估计可能偏差较大。

在获得近似设计点过程中，那些执行接受/拒绝阶段的样本点，计算了其真实功能函数值，可以作为Kriging 模型的初始训练样本。

当MMA 中建议分布函数的方差选取合适时，寻找近似设计点过程中计算了真实功能函数值的点，一般比较接近失效域边界，将其作为训练样本是较为合适的。

选择式学习函数，以及停止学习准则。

和直接蒙特卡罗方法不同，重要抽样方法所估计的失效概率不仅和样本点是否落入失效域有关，还和样本点的重要抽样权函数()()()R x f x h x =有关。

尽管学习函数仅仅涉及到样本点的功能函数预测值和预测准确度，而与()R x 无关，但采用重要抽样方法时，已经确保了生成的样本点大部分落在对失效概率值贡献大的区域内，即()R x 取较大值的区域内。

因此，采用学习函数寻找到的训练样本也都落在对失效概率值贡献大的区域内，采用该学习函数是合理的。

算法流程
基于Kriging 模型的重要抽样方法的算法流程如下：
(1)利用Hurtado 所提方法求得近似设计点，并用计算了真实功能函数值的样本点作为初始训练样本建立Kriging 模型；
(2) ISD 函数()h x 选取为以近似设计点为采样中心的独立正态概率密度函数，生成N 个样本点12,,,N x x x ⋅⋅⋅，用Kriging 模型预测其对应的功能函数值()k g x 和
Kriging 方差()2
k k x σ，k=1,2,…,N ；
(3)利用学习函数，主动寻找并添加训练样本，不断修正Kriging 模型的参数值，直到满足停止学习准则为止，建立一个能够对样本点精确分类的Kriging 模型；
学习函数
()()()g g x U x x σ=
停止学习准则 ()min 2U x >⎡⎤⎣⎦
(4)用下式估计结构失效概率：
()()()()()()
()().1.1100k k f f k k k T k N i T i N k k f x f x P P I g x I g x N h x h x ∈≤≤∈≤≤⎡⎤≈=≤+≤⎢⎥⎣⎦∑∑ 其中T 为训练样本指标集，
()()0k I g x ≤为示性函数:如果、()0k g x ≤则()()01k I g x ≤=，否则()()00k I g x ≤=。