马君明，李 惠,2，兰成明，刘彩平

(1. 北京科技大学国家材料服役安全科学中心，北京 100083；2. 哈尔滨工业大学土木工程学院，黑龙江，哈尔滨 150090；3. 北京科技大学土木与资源工程学院，北京 100083)

随着社会经济的发展，工程结构的安全性日益为人们所关注，土木工程结构可靠度研究的重点是量化结构失效事件中的不确定性[1-2]。不确定性本质上是对各类信息或认识缺失的表现，如表征材料性能、几何参数、抗力计算模式及作用统计模型的随机变量都包含不确定性[3]，随机变量的不确定性可由样本数据的统计分析进行量化，通过检测、监测及验证荷载试验等手段获取的可减小随机变量不确定性的额外信息一般称为观测信息，结合观测信息采用Bayesian 理论对随机变量已有信息进行更新，可以提升对随机变量不确定性的认知[4-9]，降低不确定性的影响，实现基于观测信息的结构可靠度更新。

基于观测信息的可靠度更新求解算法方面，Schall 等[10]提出使用二阶矩方法求解条件失效概率分母中的积分常数，但在处理多个观测信息时，确定联合事件域内的“设计点”存在困难，且该近似方法的求解误差难以定量评估。对于特定荷载作用下结构未发生失效的条件失效概率，即验证荷载下的条件失效概率，Der Kiureghian 等[11]采用随机模拟选择安全域内的随机样本实现可靠度更新，研究未获取随机变量后验的分布特征。Straub 和Der Kiureghian[12-13]首 先 将Bayesian 网络引入结构可靠度更新，通过对等式观测信息离散化建立条件概率表进而构建Bayesian 网络，提出节点消除算法实现基于Bayesian 网络的结构可靠度更新。孙鸿宾等[14]在中国较早提出采用Bayesian网络处理已知检测信息条件下的可靠度更新，研究中未考虑随机变量相关性对更新结果的影响。由于Bayesian 网络在处理相关随机变量时需要将随机变量的联合分布转化为随机变量间的条件分布，转化过程较为复杂，因此Bayesian 网络在处理含有相关随机变量观测信息的可靠度更新时计算量偏大且计算效率不高。随机变量后验概率密度函数分母中的积分常数是结构可靠度Bayesian更新计算的难点，而通过抽样得到服从后验分布的随机样本可以回避求解积分常数[15-17]。MCMC方法作为一般随机抽样方法可用于解决该问题，但是MCMC 算法存在较长的初始波动段，对于多随机变量概率密度更新问题，其计算量随着随机变量数目的增加而急剧增大[18]，且进行统计量估计时，样本的自相关性会影响估计值的准确性，设置采样间隔可以降低样本的自相关性，但计算量显著增加。理论上拒绝抽样算法可以很好地解决该问题，拒绝抽样可产生服从不具有完备形式的概率密度函数(实数域上积分不为1)的随机样本，根据目标分布构造相应的样本接受域，选择符合目标分布的随机样本，实现随机变量后验概率密度函数中分子核密度的随机抽样[19-22]。

本文首先基于Bayesian 理论建立考虑观测信息的结构体系失效概率更新模型，根据观测信息事件类型分别建立不等式和等式观测信息条件下随机变量的似然函数并推导其后验概率密度函数，建立不等式和等式观测信息条件下的观测信息域，基于观测信息域确定随机变量后验样本的拒绝抽样策略，在理论层面将不等式观测信息与等式观测信息统一在结构体系可靠度更新模型中；其次研究拒绝抽样算法的抽样效率，推导更新后结构体系失效概率估计值及其标准差的计算公式，形成完整的结构体系可靠度Bayesian 更新模型及其拒绝抽样策略。以刚架结构发生塑性破坏体系可靠度更新为例，深入研究抗力检测值、检测误差及验证荷载对随机变量后验概率密度与体系失效概率的影响，并与MCMC 计算结果进行对比，验证基于拒绝抽样的体系可靠度更新结果的准确性。

1 结构体系可靠度更新模型

弹塑性结构体系通常以塑性铰增加，结构转变为机构定义结构体系失效[23-24]。由于结构体系失效模式数量庞大，一般通过对失效概率较大的显著失效模式进行研究确定结构体系的失效概率，某一显著失效模式是由失效单元(构件)组成的并联模型，结构体系失效可视为由显著失效模式组成的串联模型[25]。定义第i个显著失效模式为事件Fi，则结构体系失效事件F为其显著失效事件Fi的并集，可表示为：

式中，k为结构形成显著失效模式的个数。则结构体系失效概率Pf(F)可表示为：

式中：f(x) 为随机变量X={X1,X2,···,Xn}的联合概率密度函数； ΩF为失效域[26]。对于串联体系，Ditlevsen 推导体系失效概率的窄界估计如下[27]：

式中，Pf(Fi) 为第i个显著失效模式的失效概率，按降序排列。

结构服役过程中可以通过各种观测手段获得影响可靠度随机变量的观测信息，根据随机变量联合概率密度函数(先验概率密度函数)f(x)及随机变量观测信息的似然函数L(x)，采用Bayesian理论更新随机变量的后验概率密度函数f′(x)[28]：

后文将根据观测信息似然函数的表达形式将其分为不等式观测信息和等式观测信息逐一讨论。

1.1 不等式观测信息的似然函数

不等式观测信息一般是指观测到某一确定荷载工况条件下(部分随机变量取确定值)结构某一特定功能未发生失效，即所谓“验证荷载”。这类观测信息的似然函数可通过功能函数g(x)进行表达，不等式观测信息的似然函数L1(x)记为：

式中：x0表示某一确定荷载工况条件下部分随机变量取确定值；I[·]为示性函数：

式中，g(x0)>0和g(x0)≤0分别表示在该确定荷载条件下结构特定功能未失效和失效两种状态。该类观测信息的似然函数以不等式形式表示，因此称为不等式观测信息。

1.2 等式观测信息的似然函数

等式观测信息是指通过特定观测手段获得结构服役过程中随机变量X或其函数h(X)的观测值hm(x) ，获取的信息必然存在观测误差 ε，这类观测信息的似然函数采用观测误差 ε的概率密度函数fε(·)的等式形式表达，因此称为等式观测信息。等式观测信息的似然函数L2(x)记为：

由式(6)和式(8)可以看出，观测信息的似然函数表征观测信息的范围。其中，不等式观测信息的似然函数为示性函数，提供随机变量取值范围的信息；等式观测信息的似然函数为观测误差的概率密度函数，提供随机变量取值的概率信息。不等式观测信息一般用于表征荷载随机变量的观测信息，而等式观测信息一般用于表征抗力随机变量的观测信息。

2 结构体系可靠度更新的拒绝抽样算法

理论上，式(4)中随机变量后验概率密度函数分母中的积分常数难以通过解析计算获得，为避免复杂的积分运算，可对式(4)中分子核密度函数进行随机模拟获得服从后验概率密度函数的随机样本。通过拒绝抽样算法获取服从式(4)的随机样本的抽样策略如下[19]：随机样本x¯服从先验概率密度函数f(x) ，存在满足A≥max[L(x)]的任意常数A，若满足Au≤L(x¯) ，则随机样本x¯服从后验概率密度函数f′(x) ，其中，u为(0, 1)均匀分布的随机数，证明如下：

显然，常数A与后验概率密度函数的分母共同决定拒绝抽样算法的抽样效率，因为，后者是常数，所以，A在满足A≥max[L(x)]前提下，应尽可能取接近似然函数上确界的数值以提高拒绝抽样算法获取后验样本的效率。

定义式(9)中样本选择域为观测信息事件域ΩZ={Au≤L(x)|x,u}，则可靠度更新可表达为观测信息事件Z条件下的结构失效概率：

体系失效是由其显著失效模式控制，根据式(11)可以看出，体系可靠度更新过程中假定显著失效模式未发生变化。实际观测信息会导致体系不同失效模式的失效概率发生变化，因此，更新过程中需要保证可靠度更新前、后所选取的显著失效模式未发生变化。

2.1 不等式观测信息条件下后验样本的拒绝抽样

根据式(4)和式(6)，不等式观测信息条件下随机变量X的后验概率密度函数可表示为：

从不等式观测信息的后验概率密度函数可以看出，先验概率密度函数与似然函数的乘积表示以示性函数重新定义先验概率密度函数的定义域。采用式(9)的抽样策略时，由于似然函数为示性函数，其最大值为1，所以A取为1 时，式(9)可表示为：

2.2 等式观测信息条件下后验样本的拒绝抽样

根据式(4)和式(8)，等式观测信息条件下随机变量X的后验概率密度函数可表示为：

为了获取服从f′(x)的随机样本，可采用式(9)抽样策略。其中，先验样本根据随机变量的先验概率密度函数产生，即使随机变量间存在相关性，也无需确定随机变量之间的条件概率分布。由于，拒绝抽样算法需要满足A≥max[L2(x)]，因此，需要确定似然函数上确界，若A不满足此要求，则获取的样本无效。

观测信息域 ΩZ是在随机变量X基础上加入均匀分布随机变量u构成的事件域：

由式(16)可知，拒绝抽样算法将等式观测信息转化为不等式观测信息。

2.3 体系可靠度更新与误差分析

可以看出，结构体系条件失效概率估计值的标准差是 Pr(Z)的函数，观测信息事件数量越多，则似然函数对应的观测信息域交集越小， Pr(Z)减小，对应的条件失效概率估计值的标准差将增大。因此，当观测信息事件数量较大时，为降低条件失效概率估计值的标准差，需要增加服从先验概率密度函数随机样本的数量。根据式(19)，可得出满足条件失效概率估计值变异系数VPˆf(F|Z)所需的先验样本数目ns：

显然，当 Pr(F∩Z)处于(10-4, 10-3)时，为满足条件失效概率估计值的变异系数不超过0.1，至少需要先验样本的数量达到(105, 106)个，后续拒绝抽样随机模拟过程中先验样本数量取为107个。后续计算分析中，为验证拒绝抽样算法模拟的准确性，进一步采用MCMC 算法中的改进Metropolis-Hastings (MH)算法对结构体系的失效概率进行更新[29]。改进MH 算法通过构造稳态分布为目标分布的马尔可夫链生成随机样本，其中目标分布为随机变量的后验分布，如式(12)和式(15)所示。为便于比较拒绝抽样算法与改进MH 算法的计算效率，且保证更新后失效概率具有相近的精度，计算时两种算法均获取相同的后验样本数量，即调用功能函数的次数相同，根据所需抽样样本的数量对比拒绝抽样算法与改进MH 算法的计算效率。

3 单层刚架体系可靠度更新

如图1 所示某单层刚架受到水平荷载H和竖向荷载V作用，其中5 个关键控制截面的塑性抵抗矩Ri(i=1,2,···,5)互为相关的随机变量，各随机变量分布类型及参数取值在表1 中给出，具体详见文献[11]。

图1 单层刚架模型及其3 种失效模式Fig. 1 Single-story planar frame and its failure mechanisms

表1 各随机变量分布类型及分布参数Table 1 Distribution types and parameters for random variables

刚架发生塑性破坏时3 种显著失效模式如图1所示，体系失效是由失效单元和显著失效模式组成的串-并联模型，相应的功能函数及失效域为：

采用直接Monte Carlo 方法计算得到3 种显著失效模式对应的失效概率分别为：1.016×10-2、1.755×10-3和2.342×10-2，相应结构体系发生塑性破坏的失效概率为2.644×10-2，根据式(3)估计得到失效概率的界限为(2.639×10-2, 2.647×10-2)。

3.1 不等式观测信息的可靠度更新

通过验证荷载对刚架承载能力进行检测，已知竖向荷载V取v0=100 kN时刚架未失效，则由式(14)与式(21)可确定其不等式观测信息事件域为ΩZ={∩gi>0|v0=100 kN},(i=1,2,3)。首 先 根据表1 给出的随机变量分布类型及分布参数生成水平荷载H和截面塑性抵抗矩Ri(i=1,2,···,5)的先验样本，根据观测信息事件域 ΩZ和式(12)计算各随机变量的后验概率密度函数，即满足观测信息事件域随机变量的先验样本作为该观测信息条件下随机变量的后验样本。根据拒绝抽样算法获得的各截面塑性抵抗矩的先验、后验样本及其概率密度函数(PDF)在图2 中给出。

从图2(a)～图2(e)可以看出，各截面塑性抵抗矩R1、R2、R3、R4和R5的后验概率密度均向右偏移，刚架在验证荷载v0=100 kN作用下未失效，各截面塑性抵抗矩取大值的概率增加。根据观测信息域可知，截面塑性抵抗矩R3对功能函数g2和g3的影响最大，截面塑性抵抗矩R3对验证荷载v0最敏感，因此，截面塑性抵抗矩R3的后验概率密度变化最大。从图2(f)可以看出，水平荷载H的后验概率密度略向左偏移，即在验证荷载v0=100 kN 作用下结构未失效，水平荷载H取大值的概率降低。

图2 不等式信息观测条件下各随机变量的先验和后验样本及其概率密度函数Fig. 2 Prior and posterior samples and their PDFs of random variables under inequality observation information condition

由式(17)计算更新后结构体系的失效概率为2.279×10-4，此时根据后验样本估计显著失效模式g1的失效概率等于0，因此采用式(3)得到结构体系失效概率的上、下界限与更新后结构体系的失效概率相同。可以看出，由于增加观测信息，即在验证荷载v0=100 kN作用下结构未失效，更新后结构体系的失效概率与原失效概率相比显著减小。采用改进MH 算法对结构体系的失效概率进行更新，计算得到失效概率为2.320×10-4，两种算法得到更新后体系失效概率较接近，但为降低MH 算法后验样本的自相关性并获得与拒绝抽样相同的后验样本数量，MH 算法需要抽样样本数量远大于拒绝抽样。

3.2 等式观测信息的可靠度更新

采用相关检测方法对截面塑性抵抗矩R4和R5进行检测，得到R4和R5的检测值分别为r4,m=150 kN·m和r5,m=200 kN·m ，已 知 检 测 误 差εRi服从正态分布 εRi～N(0,15)，i=4,5，即检测误差的均值为0，标准差为 15 kN·m ，且R4和R5的检测误差互相独立。根据式(8)，观测信息的似然函数表示为：

为提高后验样本的接受概率，常数A选取似然函数的最大值A=max[Lr4,r5(r4,r5)]=7.074×10-4。图3 给出1000 个先验样本的拒绝抽样过程，图中全部数据点为服从先验概率密度函数的R4、R5的随机样本，曲面下方为样本选择域：Au¯≤Lr4,r5(r¯4,r¯5)，则曲面下方的数据点即为选择得到的服从后验概率密度函数的随机样本。各截面塑性抵抗矩的先验、后验样本及其概率密度函数在图4 中给出。

图3 拒绝抽样随机样本选择示意图Fig. 3 Schematic diagram of rejection sampling for random samples

从图4(d)和图4(e)中可看出，由于增加截面塑性抵抗矩R4和R5的观测信息，因此两者的后验概率密度集中于相应的检测值附近。虽然增加的观测信息与R4和R5直接相关，但由于各截面塑性抵抗矩的先验分布是正相关的，使得截面塑性抵抗矩R1、R2和R3的后验概率密度也发生变化，相对于先验概率密度向右偏移，均值略有增大，见图4(a)～图4(c)。由于荷载随机变量V和H与各截面塑性抵抗矩互相独立，因此，这两个随机变量的后验概率密度不受R4和R5观测信息的影响。

图4 等式观测信息条件下截面塑性抵抗矩的先验和后验样本及其概率密度函数Fig. 4 Prior and posterior samples and their PDFs of resistance moments under equality observation information condition

由式(17)计算得到更新后结构体系的失效概率为6.971×10-3，估计结构体系失效概率的界限为(6.967×10-3, 6.972×10-3)，采用改进的MH 算法获得更新后结构体系的失效概率为6.900×10-3，两种算法得到更新后体系失效概率较接近，但为降低MH 算法后验样本的自相关性，需要增加其采样间隔，因此抽样样本数量远大于拒绝抽样。可以看出，由于增加观测信息，截面塑性抵抗矩R5的后验均值显著增大，且R4和R5后验变异性降低，使得更新后结构体系的失效概率显著降低。

显然，检测误差将对随机向量后验概率密度及失效概率产生较大影响，其他条件不变的前提下，以检测误差 εRi的标准差 σεR为变量，研究检测误差对随机变量后验均值、标准差及失效概率的影响，相应结果在图5 和图6 中给出。从图5 中可以看出，随着截面塑性抵抗矩检测误差标准差的增大，R1、R2、R3和R5的后验均值逐渐减小，由于截面塑性抵抗矩检测误差的不确定性增大导致R4和R5的后验标准差逐渐增大，使得更新后结构体系的失效概率增大(见图6)。从图5 可以看出R4的后验均值随检测误差标准差的增加略有增大，其原因是随着观测信息的不确定性增大，随机变量相关性对其后验分布影响增大，导致其后验均值略有增大。从上述分析可以看出，为获取更有价值的观测信息更新随机变量的后验分布，进而更新结构体系的失效概率，必须严格控制相关检测误差的标准差，否则观测信息不确定性增大对结构体系失效概率的更新没有实际意义。

图5 检测误差标准差对随机变量后验参数的影响Fig. 5 Influences of standard deviation for detection error on parameters of posterior distribution for random variables

图6 检测误差标准差对更新后结构体系失效概率的影响Fig. 6 Influences of standard deviation for detection error on updated failure probabilities

4 双层刚架体系可靠度更新

某双层刚架计算简图如图7 所示，各杆件的塑性抵抗矩M1、M2、···、M5及荷载S1、S2、···、S4均为相互独立的正态随机变量，相应分布参数见表2。根据文献[30]，该双层刚架的8 种显著塑性失效模式的功能函数及失效域如下：

表2 各随机变量分布参数Table 2 Parameter values of random variables

图7 双层刚架计算简图 /mFig. 7 Schematic diagram of two-story planar frame

采用直接Monte Carlo 方法计算结构体系发生塑性破坏的失效概率为3.971×10-2，相应界限估计为(3.823×10-2, 4.058×10-2)。

4.1 不等式观测信息的可靠度更新

已知验证荷载S1、S2和S3分别为s1,0=200 kN、s2,0=120 kN和s3,0=150 kN时刚架未失效，则其观测信息事件域 ΩZ根据式(14)和式(23)确定，拒绝抽样过程详见3.1 节，拒绝抽样算法获得的各杆件塑性抵抗矩的先验、后验样本及其PDF 在图8中给出。从图8 可以看出，M1、M2、M3、M4和M5的概率密度略向右偏移，而S4的概率密度向左偏移，由于在该验证荷载条件下结构未发生失效，说明表征抗力的随机变量取大值的概率增加，荷载S4取大值的概率降低。获取后验样本后，采用拒绝抽样算法更新后结构体系失效概率为1.761×10-3，估计结构体系失效概率的界限为(1.756×10-3, 1.768×10-3)，更新后结构体系的失效概率显著减小。采用改进的MH 算法获得更新后结构体系的失效概率为1.764×10-3，但为降低MH算法后验样本的自相关性，需要增加其采样间隔，因此抽样样本数量远大于拒绝抽样。

图8 不等式信息观测条件下各随机变量的先验和后验样本及其概率密度函数Fig. 8 Prior and posterior samples and their PDFs of random variables under inequality observation information condition

4.2 等式观测信息的可靠度更新

采用相关检测方法对杆件塑性抵抗矩M1和M3进行检测，得到M1和M3的检测值分别为m1,m=120 kN·m和m3,m=110 kN·m，已知检测方法的检测误差 εMi服从正态分布 εMi～N(0,10)，且检测误差互相独立，则观测信息的似然函数与式(22)类似。根据拒绝抽样算法获得M1和M3的先验、后验样本及其概率密度函数在图9 中给出，由于各随机变量相互独立，因此观测信息仅影响M1和M3的后验概率密度，同时M1和M3的先验概率密度相同，而M1的检测值比M3的检测值大，在相同检测误差条件下，M1的后验概率密度相比于M3更趋向取大值，从图9 中可以清晰看到这点。采用拒绝抽样算法更新后结构体系的失效概率为3.862×10-3，相应上、下界限为(3.662×10-3, 3.956×10-3)。采用改进MH 算法获得更新后结构体系的失效概率为3.899×10-3，两种算法得到的更新后体系失效概率较接近。仅改变抗力检测误差的标准差 σεM，随机变量后验均值和标准差随检测误差标准差的变化曲线在图10 中给出，相应更新后结构体系失效概率及其界限估计在图11中给出。

图9 等式观测信息条件下随机变量的先验和后验样本及其概率密度函数Fig. 9 Prior and posterior samples and their PDFs of random variables under equality observation information condition

从图10 可以看出，随着检测误差的标准差增大，M1和M3的后验均值逐渐减小，后验标准差逐渐增大，这是由于检测误差的标准差增大，观测信息的不确定性增大，M1和M3先验信息的影响增大，其后验分布逐渐趋向先验分布。显然，M1和M3的后验均值减小且标准差增大，导致更新后的失效概率逐渐增大(如图11 所示)。由于随机变量M1和M3的先验分布及似然函数均为正态分布，且相互独立，根据共轭性质可知其后验分布依然服从正态分布，其后验标准差 σ′M是 σM和 σεM的函数[31]：

图10 检测误差的标准差对随机变量后验参数的影响Fig. 10 Influences of standard deviation for detection error on parameters of posterior distribution for random variables

图11 抗力检测误差标准差对结构体系失效概率的影响Fig. 11 Influences of standard deviation for detection error on updated failure probability

已知M1和M3的先验标准差 σM相同，检测误差的标准差 σεM亦相同，由式(24)得出M1和M3后验标准差也应该相同。图10 采用拒绝抽样得到的M1和M3后验标准差几乎完全相同，进一步验证了拒绝抽样算法的准确性。

5 结论

本文着重研究基于观测信息的结构体系可靠度更新模型及其拒绝抽样算法。基于Bayesian 理论建立结构体系可靠度更新模型，进一步根据观测信息事件类型建立不等式和等式观测信息的随机变量似然函数及其后验概率密度函数，确定观测信息域及后验样本的拒绝抽样策略，估计拒绝抽样算法的抽样效率，推导更新后结构体系失效概率估计值及其标准差的表达式。以刚架结构塑性失效为例研究不同观测信息对随机变量后验概率密度及体系失效概率的影响。得到如下主要结论：

(1) 不等式与等式观测信息的区别在于观测信息似然函数的形式不同，不等式观测信息对应的是示性函数，提供随机变量取值范围信息，而等式观测信息对应的似然函数是概率密度函数，提供随机变量取值的概率信息，拒绝抽样算法最终将等式观测信息转化为不等式观测信息进行处理，根据先验样本满足相应观测事件域实现后验样本的抽样。

(2) 为提高等式观测信息后验样本的接受率，需要控制拒绝抽样算法中常数A尽可能取等大于观测信息似然函数上确界的某一小值。

(3) 拒绝抽样算法与改进MH 算法计算得到的更新后结构体系失效概率较为接近，但拒绝抽样无需考虑后验样本的自相关性，可以显著提高抽样效率。

(4) 与抗力相关的随机变量检测值提高可以降低更新后结构体系的失效概率，但需要降低检测误差的标准差，以减小随机变量后验分布的变异性，否则观测信息不确定性增大导致失效概率的更新没有实际意义。

(5) 验证荷载值增加可以显著降低更新后结构体系的失效概率，此时与相应验证荷载较敏感的表征抗力的随机变量后验概率密度变化较显著。