考虑抗力非平稳劣化和荷载相关性的桥梁可靠性分析

金聪鹤，钱永久，徐望喜，黎 璟

(1.西南交通大学土木工程学院，四川，成都 610031；2.蜀道投资集团有限责任公司，四川，成都 610093)

我国公路桥梁受到交通量和负载持续增长的车辆荷载作用[1]，长期处于“带病工作”的状态，安全问题日益突出[2-3]。据统计，2008 年-2010 年，我国公路危桥数量由6000 座猛增至9.35 万座[4-5]；2005 年-2018 年期间，全国至少发生28 起大中型桥梁垮塌事故，其中35%由车辆超载引发，80%的桥梁剩余使用寿命不足30 年[6]。车辆超重、交通量急剧增长等因素使服役桥梁损失的承载力远超预期，增大了车辆和行人通勤的风险，如不采取措施，可能造成难以挽回的损失[7]。然而，对所有公路桥梁都进行彻底的检测和加固，既不经济，也在技术上难以实现[8]。可取的做法是：先对桥梁结构进行基于可靠性理论的安全性能评估，以此作为其后续养护维修的科学依据[3-4,8-9]。

对于在役桥梁，由于设计时对材料强度、车辆通勤频率和荷重的增长缺乏预估，以致不少桥梁在设计使用年限内已出现了混凝土开裂、钢筋锈蚀等病害[10]。考虑到我国公路桥梁服役状况十分严峻，若沿用《公路工程结构可靠性设计统一标准》[11]中基于极限状态设计原则的瞬时可靠度指标计算方法，不能考虑结构抗力随时间劣化的效应，不利于对在役桥梁进行安全状况评估。

为了考虑桥梁承载力的时变效应，1993 年MORI 和ELLINGWOOD[12]基于考察时段内结构的首次穿越问题，假设荷载服从平稳Poisson 过程且彼此独立，提出了时变可靠度理论，并用MCS法对计算结果进行了验证。时变可靠度理论是基于参数随机性评估结构剩余寿命的最佳工具之一[13]，具备广泛的工程前景与重大经济社会效益，成为了国内外十分关注的课题[14-15]。王草等[3-4]和LI等[16]提出了独立荷载基于非平稳Poisson 过程的时变可靠度计算方法；叶新一等[9]将可靠度积分转化为近似度较高的代数运算；李全旺和王草[17]、金聪鹤等[8]基于MCS 实验讨论了荷载相关性对结构可靠度的影响；YUAN 等[18-19]基于Markov 过程和Gamma 过程建立了桥梁时变抗力的非平稳劣化模型；WANG[20]提出了基于复合Poisson 过程的时变抗力冲击衰减模型；金聪鹤等[15]将Poisson参数λ 表示为频率函数λ(t)，提出了荷载发生频率随时间增大的时变可靠度分析方法；徐望喜等[21]基于时变可靠度理论中车载分布的非平稳极值效应，运用改进广义帕累托分布模型对随机车流作用下的荷载效应进行极值外推分析；另有学者[14,22-24]考虑历史车辆荷载对服役桥梁时变抗力的验证作用，提出了基于历史荷载信息的结构时变可靠度分析方法。

除了上文提到的穿越概率法和数值模拟法，近些年国内外学者还提出了许多新颖的时变可靠度分析模型，包括代理模型法和准静态法等[25-26]。HAWCHAR、DAI 和ZHENG 等[27-29]提出基于多项式混沌理论(polynomial chaos theory)和Karhunen-Loeve 展开式的时变可靠度代理模型；ZHAO 等[26]基于包络函数算法提出了考虑认知不确定性(epistemic uncertainty)的结构使用寿命预测方法；周建方等[30]对基于模糊可靠度理论的几种准静态算法进行了对比；宋帅等[31]基于Copula 函数建立了桥梁地震易损性分析模型；AMINI 等[32]研究了受益于多项式混沌Kringing(PCK)模型的自适应可靠性算法，并利用C-vine 和D-vine Copula 理论研究了随机变量之间的非线性依赖关系。

多荷载分布的概率计算涉及复杂的积分运算，且时变抗力和荷载通常均不服从正态分布[7,21,33]，因此，提出合理的可靠度计算方法，兼顾运算效率和解决实际工程问题是当下的研究重点，亦是本文的核心。当斯—托马斯悖论(Downs-Thomson paradox)指出，过路费的增加和路面运力的改善可能造成交通拥堵[34]。若服役桥梁处于拥堵路径之中，则其负载与交通主管的定期策略相关；且作用于工程结构的荷载通常认为具备时间相关性[8,17,35-37]。因此，桥梁安全性能评价不能忽视车辆荷载时间相关性的影响。另外，构件抗力随时间的劣化是单调不可逆的[38]，更符合一个“经历‘λT’次衰减总共需要消耗多少抗力”的Gamma 过程[18,39]。综上，本文基于时变可靠度理论，将劣化抗力表示为单调递减的Gamma 过程；考虑结构服役期间荷载均值增大的效应以及它们之间的时间相关性，提出了基于时变抗力和车辆荷载同为非平稳过程的桥梁时变可靠度分析方法，并给出了计算公式。

1 桥梁时变可靠度分析理论

在时变可靠度分析中，穿越概率指考察时段内荷载首次超越抗力的概率，如图1 所示。R(t)和S(t)分别表示时变抗力和连续荷载过程，功能函数写为Z(t)=R(t)-S(t)，考察时段为(0,T)，记T0=0。实际上，可能使结构失效的荷载通常为强度较高、频率较低的瞬时荷载，因此通常采用Poisson 随机过程予以描述。假设这段时间发生了n次荷载作用S1,S2,···,Sn，任意t∈(0,T)时刻的时变抗力R(t)可以表示为初始抗力R0与抗力衰减函数g(t)的乘积：R(t)=R0·g(t)。荷载Si对应的发生时刻Ti服从(0,T)的均匀分布[9](i=1, 2, ···,n；T1

图1 时变抗力与荷载Poisson 过程Fig.1 Time-dependent resistance and Poisson process of load

若考察时段任意荷载均不大于对应时刻的时变抗力值，则结构在考察时段内可靠，它的概率Pl(T)=P[R(T1)>S1∩R(T2)>S2∩···∩R(Tn)>Sn] =P[Z(t)>0,∀t∈(0,T)]表示结构在(0,T)的时变可靠度。从而穿越概率即结构在(0,T)的时变失效概率，记为Pf(T)=1-Pl(T)。设n个荷载的联合累积分布函数为FS1,S2,···,Sn(·)，Pl(T)可以表示为：

对于任意荷载Si，其发生时刻Ti的联合密度函数fTi(t)为：

从而式(2)可写作：

联立式(2)和式(5)，由全概率公式得[4,40]：

式(6)成立的必要条件为P[N(T)=0]；结构在考察时段发生首次超越事件即告失效，因此抗力大 于 荷 载 的 时 长Ts可 以 表 示 为将Ts代入式(1)同样可以得到上述结果。考虑初始抗力R0的随机性，如图1 所示，设其概率密度函数为fR0(r)，则式(6)可写为：

式(7)即假定荷载为平稳Poisson 过程的时变可靠度计算公式。将FS(·)的均值参数表示为一个与时间t相关的随机变量，可以扩展到荷载非平稳Poisson 过程的可靠度计算中[3]，不过这种算法仅考虑了荷载均值增大的情况，并未考虑抗力劣化的非平稳性。

对式(7)的计算结果进行MCS 实验，由于荷载的频率已知，而发生时间未知，因此将目标考察时段划分为n个相等时长的子区间，(0,T1], (T1,T2], … , (Tn-1,T]，如图2 所示。其中第i个子区间的最大荷载为Si(0

图2 时变可靠度MCS 模型Fig.2 MCS model of time-dependent reliability

若荷载之间具备相关性，由于FS1,S2,···,Sn(R1,此时式(7)不再成立，需要寻求其他结构时变可靠度评估方法。Copula 函数理论能够准确地描述变量的相关性，是解决高维随机变量联合概率分布问题的常用方法[41]。根据Sklar 定理的n维情形，对实数域内的任意随机变量si∈，总是存在一个n元Copula 函数满足：

式(8)表明：n维荷载变量的联合分布总是客观存在的。设ri=r·g(ti)表示时变抗力ri的随机变量形式，则任意子时段荷载发生时刻的联合密度函数为：联立式(2)与式(8)得：

式(9)即考虑n维相关荷载的时变可靠度计算公式。当n>2 时，若关联的多元荷载服从正态分布或t分布，可以采用金聪鹤等[8]提出的Copula随机数方法给出式(9)基于MCS 的模拟解，并同样可以扩展到荷载非平稳过程的模拟，除此以外很难通过直接计算Copula 函数积分获得数值解。基于动态称重系统(WIM)的研究表明，作用于桥梁的车辆荷载效应服从极值I 型分布[7,21,42]，上述方法便不再适用。Gamma 过程可以很好地应用于模拟结构性能的退化，例如腐蚀、开裂和服役性能降低[43]。YUAN 等[18]和YANG 等[44]提出了利用Gamma 过程更新结构抗力的方法；WANG[45]将子区间(Ti-1,Ti]抗力的劣化值 Δi表示为服从Gamma分布的随机变量，将每个子时段内结构的可靠度作为独立事件，依据条件概率分布，提出了将抗力表示为可叠加的Gamma 分布的时变可靠度计算方法：

式(10)是一个贝叶斯公式，表示P(Ri-Δi+1>Si+1)是P(Ri>Si)的条件概率，并且与其他抗力大于荷载事件相独立。根据Gamma 分布的可加性，若Δi服从形状参数和统一尺度参数分别为 ηi和ξ 的Gamma 分布，那么服从均值和方差分别为的Gamma 分布。不过，该模型并未涉及荷载之间的相关性，以及尺度参数变化对可靠度计算结果的影响。为了避免三重以上积分的求解，本文假设仅相邻荷载之间具备相关性，并与其他相邻荷载的结构可靠性互为独立事件，对服从极值I 型分布车辆荷载的二元Copula 联合分布函数(JCDF)进行积分，基于贝叶斯原理，计算并讨论了荷载相关性和抗力劣化的非平稳性对时变可靠度结果的影响；并采用MCS 方法对荷载完全独立时的计算结果进行了验证，具体步骤可参考文献[3]。

2 考虑抗力Gamma 劣化和关联荷载的时变可靠度模型

服从极值I 型分布的二维随机变量Si和Si+1的Copula JCDF 可写为[46]：

式中，ui和αi为边缘分布Si的参数，按下式计算：

Si的边缘分布写为：

式中：μSi和 σSi分别表示荷载Si的均值和标准差；γ ≈0.5772为Euler 常数。为了考虑荷载随结构服役增大的特性，设年均增长率为ε，荷载均值按照线性增长进行估算：

式中：μS0为初始荷载均值；系数χi,i+1按下式计算：

式中，ρi,i+1表示相邻荷载的相关系数。假设相邻荷载具备时间相关性，可以按照指数衰减规律进行描述[17,36]：

式中：m为参数；ti为服从(0,T)均匀分布的随机变量。当相邻荷载发生时刻间隔越长，相关性越低，符合期望。根据贝叶斯原理，若P(Ri+1=Ri-Δi+1>Si+1)可视为P(Ri>Si)的条件概率， 则P(Ri>Si,Ri-Δi+1>Si+1)也可作为P(Ri>Si)的条件概率。从而将可靠度表示为：

式中：n=λT；fRi(r)为第i个子时段服从对数正态分布抗力Ri的边缘分布：

式中，θi和ζ分别为抗力Ri的分布参数，按下式计算：

式中：μR0为初始抗力R0的均值；VR为与混凝土材料相关的抗力变异系数；fΔi(δ)为服从Gamma分布的抗力劣化值 Δi的密度函数：

即以每个时间间隔基于衰减函数的劣化值作为形状参数。式中，抗力衰减函数g(t)按照下式进行计算[16]：

式中，a和b为参数。如图3 所示，为了使抗力劣化的Gamma 过程与衰减函数g(t)的1 倍相对应，根据Gamma 分布的可加性，fRi,Gam(r)作为服从Gamma 分布时变抗力的均值和方差分别为表示为：

图3 抗力Gamma 劣化模型Fig.3 Gamma degradation model of resistance

依据Gamma 过程的性质，当尺度参数ξ 足够小时，Ri的Gamma 分布集中在均值附近，使得劣化过程非常接近衰减函数；反之，ξ 越大，则Ri的分布越松散，抗力在短时间骤降的可能性增大，突变性增强。注意到对x＜0.3 有e-x≈1-x，从而式(11)可简化为：

联立式(17)、式(18)、式(20)、式(24)，得到关联荷载基于抗力非平稳劣化的桥梁时变可靠度算式：

在具体应用时，先计算n的值及每个子时段Ri>Si的概率fRi(r)FSi(r)dr，将结果保存起来；然后计算分母的二重积分结果，再将所得子区间的条件概率值连乘，并与fR1(r)FS1(r)dr的结果相乘，最终得到结构的可靠度。当二重积分FSi,Si+1(r,r-δ)fRi(r)fΔi(δ)drdδ的结果不小于0.7 时，可以依据式(25)计算结构基于二维关联荷载和抗力非平稳劣化时变可靠度。式(25)可从三个参数角度来分析结构的可靠性变化，分别是：抗力劣化的尺度参数ξ、荷载年均增长率ε 以及荷载相关性参数m。在适用性和一般性放面已达到足够的评估需求；荷载频率λ只影响发生次数n的期望，并不直接参与式(25)的积分运算，因此不再考虑其为时间函数λ(t)的情形。

3 实例分析

四川省某II 类环境条件钢筋混凝土简支梁桥，上部结构为T 型梁，计算跨径20 m，横截面尺寸如图4 所示。混凝土强度等级C40，纵向钢筋和箍筋等级分别为HRB335 和 HPB300。依据《公路桥涵设计通用规范》[47]和《公路桥梁技术状况评定标准》[48]对该桥主梁进行自然电位法检测。结果表明：下部受拉钢筋表面有锈斑，但锈蚀深度较浅，初步判定为结构抗力开始衰减。现基于时变可靠度理论，对该桥梁未来30 年内的安全性能进行基于本文可靠性模型的评估。

图4 主梁横截面图 /mmFig.4 Cross section of main girder

依据《公路钢筋混凝土及预应力混凝土桥涵设计规范》[49]对主梁进行正截面承载力复核，得到除去自重的初始抗力为R0=6000 kN·m，以跨中截面弯矩进行表示，服从变异系数为0.15 的对数正态分布[33]。由于主梁为适筋梁，仅考虑下部受拉钢筋完全拉断的正截面承载力失效，忽略剪切变形和偏心受压等复杂的失效情形。根据当地气象条件，得出g(t)的参数a=0.0005，b=1.8。据当地其他桥梁历史数据，得到2 年一遇(λ=0.5 a-1)的最大初始荷载S0为R0的33%，且年均增长率ε 介于1%～3%。

图5 和图6 给出了结构在服役第10 年(i=5)时抗力的密度函数fR5(r)以及采用尺度参数 ξ由1 增大至1000 时对应的Gamma 分布fR5_Gam(r)的比较结果。由结果可知， ξ越小，R5在其均值附近取值的概率越大，分布越密集；反之，则R5的分布越离散，取值范围增大。

图5 R5 的分布对比图(ξ=1, 10)Fig.5 Distribution comparison of R5 (ξ=1, 10)

图6 R5 的分布对比图(ξ=100, 1000)Fig.6 Distribution comparison of R5 (ξ=100, 1000)

图7 所示为不同尺度参数下对应的抗力可能发生的Gamma 劣化过程，但未示出ξ=1的结果，因为此时Gamma 过程与R(t)几乎完全重合；当ξ由10 增大到1000，每个子区间抗力劣化的突变性增大，这与图5 和图6 的结果相符。 ξ增大意味着抗力突变的概率增加，劣化的非平稳特性更显著，可能导致结构失效概率增大。在尺度参数取值方面， ξ取1000 的时变抗力在20 年～22 年并未发生衰减，然而从22 年～24 年由5131 kN·m 骤降至3716 kN·m，这非常不符正常运营的桥梁承载力状态。试算表明， ξ过大可能使得式(20)中δ/ξ在n和m较小时无解。因此本文分别取 ξ取1、10、50 进行时变可靠度分析。

图7 多尺度参数的抗力可能劣化路径Fig.7 Possible resistance deterioration of multiple scale parameters

3.1 计算结果与分析

表1 展示了当ε=3%时，结构在4 年～30 年的部分计算结果，以时变失效概率Pf(T)表示。

表1 部分时变失效概率计算结果Table 1 Partial calculation results of time-dependent failure reliability

T=0 时结构完全可靠；T=(0, 2)的失效概率可通过计算P(R1>S1)直接获得。其中m=0.001 和m=10 分别表示荷载几乎完全相关以及荷载趋近于完全独立两种情形，而ξ=1 和ξ=50 则分别表示抗力衰减由平稳过程向非平稳过程的过渡。由结果可知，较强的荷载相关性降低了服役期内荷载出现极值的可能，因此结构失效概率降低，可靠性增加，这与相关研究[8,17]得出的结论一致。当ξ 增大时，结构的失效概率增加，这表明抗力劣化的非平稳性会增大结构的失效概率；ξ=1，m=0.001时，结构在(0,Ti) (i=2, 3, 4, ···, 15)的Pf(Ti)与ξ=50，m=0.001 的Pf(Ti)的比率介于1.019～1.205；而在m=10 时，ξ=50 的Pf(Ti)与ξ=1 的Pf(Ti)的比率介于1.009～1.021，说明荷载相关性是影响结构时变可靠度的主要原因，但是当相关性较强时，结构抗力劣化的非平稳特性对时变可靠度的影响不能忽略。

对式(25)进行计算时，本文采用了图2 所示的均分时间段方法，在任意子区间荷载的二维分布FSi,Si+1(r,r-δ)、抗力的分布fRi(r)和抗力劣化值的分布fΔi+1(δ)是唯一确定的(i=1, 2, 3, ···,n-1)。为了讨论时间分布的随机性对计算结果产生影响，在计算相关参数ρi,i+1时，本文分别采取第i个荷载Si发生时刻TSi为(Ti-1,Ti)中间时刻的方式及(0,T)均匀分布的方式(先生成随机数，然后对其排序)。前者将固定荷载相关参数，而后者则将荷载相关参数随机化。对比结果如图8 所示，其中ε=0.03， ξ=50，m=0.1。三个对照组均为任意荷载发生时刻服从(0,T)均匀分布的结果。

图8 荷载发生时刻固定与随机分布的结构时变失效概率Fig.8 Structural time-dependent failure probability with fixed and stochastic time of load occurrence

由图8 可知，将荷载相关参数随机化后，计算所得的时变失效概率为一个范围。考虑到荷载非平稳性参数ξ的取值受限，对照组方法可能使抗力非平稳劣化对结构可靠度的影响无法体现。因此后续计算均采用荷载发生时刻为该子区间中间时刻的方式。不过，在确定性的概率分布下，通过随机化相关性参数，便可得到一个可靠度范围，对于解决荷载中度相关的问题十分有意义，例如经常拥堵桥梁的安全性能评估问题。

图9 和图10 展示了ε=3%的所有时变失效概率计算结果。可知，当ξ 由1 分别增大到10 和50，不论荷载的相关性强弱如何，结构的失效概率都有略微的增大，但相较于荷载相关性的影响偏弱。m越大，荷载独立性越强，ξ 变化对结构失效概率的影响越不显著。这表明：不论荷载是否具备相关性，忽略桥梁抗力的非平稳劣化特性都会高估桥梁的可靠性；ξ 越大，表明抗力发生骤降的程度越高，因此结构失效的概率也越大。由表1可知，考虑抗力非平稳劣化的结构时变失效概率与不考虑的比值大于1.01，因此桥梁实际服役期间的可靠性要偏低一些。

图9 结构在16 年以内的时变失效概率Fig.9 Structural time-dependent failure probability within 16 years

图10 结构在18 年～30 年的时变失效概率Fig.10 Structural time-dependent failure probability within 18 to 30 years

图11 给出了当ξ 和m分别为50 和10 时，结构在ε 分别为1%、2%和3%的时变失效概率结果。可知：荷载均值的增长率对结构失效概率影响较大。即使不考虑荷载强度几何增长的情形，结构在T=30 a 的时变失效概率Pf(30)分别为0.023、0.128 和0.466，如果考虑荷载为几何级数的增长，或者λ 为时间相关的增函数，那么失效概率还会进一步增大。不过，当ε=3%时，荷载在T=26、28 和30 的时变失效概率增长过快，这表明该子区间的FSi,Si+1(r,r-δ)fRi(r)fΔi(δ)drdδ的积分可能小于0.7，使得计算结果不准。另外，2%的年均荷载增长率可以视为桥梁服役安全与否的阈值。若超过该值，应当限制重车超重车的负载，并定期对桥梁结构主要承重构件进行检测、维修或更换。否则结构服役30 年内的失效风险较大。

图12 给出了式(25)与MCS 结果的比较，其中m=10，表示荷载几乎彼此独立；ξ=50；MCS 的模拟次数为200 万。ε=0.03 时桥梁在T=26 a～30 a 的时变失效概率较MCS 结果偏离较大，原因是T>24 a 时FSi,Si+1(r,r-δ)fRi(r)fΔi(δ)drdδ的积分结果小于0.7，因此e-x>1-x，被积函数偏小使得计算结果较MCS 模拟结果存在较大误差；而ε=0.01 时则不存在上述问题，因此得出的计算结果与MCS 结果十分接近，证明了本文所提公式的正确性。

图12 公式解与MCS 对比Fig.12 Comparison between equation solution and MCS

4 结论

本文基于时变可靠度理论，将抗力的劣化表示为Gamma 非平稳随机过程；考虑相邻荷载之间的时间相关性，将二维荷载的可靠度作为条件概率，基于贝叶斯原理，提出了考虑抗力非平稳劣化和荷载相关性的桥梁时变可靠度分析模型，给出了计算公式，并采用MCS 方法进行了验证。结果表明，当时变失效概率计算结果小于0.2 时，本文所提公式精度较高。得出以下结论：

(1) 不考虑结构抗力劣化的非平稳性会一定程度上高估结构的可靠性，荷载相关性越强，高估的程度越大。荷载相关性对结构可靠性的影响大于抗力劣化非平稳性。

(2) 其他条件相同的情况下，劣化抗力的尺度参数越大，结构失效概率越高。因此，经受了较高强度荷载作用的桥梁结构，例如地震、大风、交通事故、渡轮撞击等，要确保对关键承重构件的及时检测、养护和更换，以防止承载力骤降威胁桥梁安全。

(3) 荷载相关性越高，结构的失效概率越低，安全性能越有保障。对此合理的解释类似本文抗力Gamma 劣化过程的分析，即关联的荷载过程降低了其在桥梁后继服役期的变异性，同时降低了其大于骤降抗力的概率。因此道路主管部门可以对超重车进行限制，增大其相关性，以提升公路桥梁安全性，例如在桥梁服役的某段时期允许重车超重车通勤。

(4) 对于在役桥梁，若已服役期的年均交通增长率超过2%，那么结构在30 年内失效风险较高，应当限制后继服役期车辆的通勤数量和轴重。