近年来我国现代化路网交通建设取得了长足进步，桥梁作为必要的组成部分，其服役状态会随着时间的推移而逐渐劣化。为保证在役桥梁的通行能力和运营安全，对其进行科学管理是非常必要的，这既需要对桥梁当前的技术状况进行准确评估，也需要建立合理的退化模型预测其未来的劣化发展情况^{[1 − 2]}。桥梁退化模型的研究主要集中在两个方向：基于机理的退化模型和基于统计规律的退化模型^{[3 − 6]}。前者从影响桥梁退化的各种因素出发，通过研究其作用机理，再利用现场数据加以修正建立退化模型^{[3 − 4]}。后者利用反映结构退化的实测桥检数据，通过统计历史数据归纳出桥梁退化规律^{[5 − 6]}。由于影响桥梁退化的因素较多，且退化过程的随机性和时变性都较为复杂，因此基于机理的模型适用于对具体某座桥梁的退化预测，若要反映某一类桥梁的总体退化规律，则基于统计规律的退化模型往往更加适用^[7]。

马尔科夫链模型^[8]是应用最普遍的基于统计规律的桥梁退化模型。该模型利用历史桥检数据归纳出某类桥梁的技术状况转移概率矩阵，用于预测该类桥梁未来的劣化发展^{[8 − 10]}。传统的马尔科夫链模型假设转移概率不随时间变化，即桥梁未来的状态仅取决于桥梁当前状态，而与停留在当前状态的时间长短无关，这显然与常识不符^[11]。为了修正这个缺陷，学者在马尔科夫链模型的框架下引入了状态停留时间，TC (Time-in-Condition)，的概念，二者相结合能够较好地反映桥梁退化的规律^[12]。然而目前估计 TC的方法对桥梁检测数据的完备性有较高的要求，需要完整记录桥梁退化到各个劣化状态的时间^[12]。如果桥梁的历史检测记录不全，或者检测数据的时间间隔过大，已有方法对桥梁退化情况的预测往往带有较大误差^[13]。

目前我国很大一部分既有桥梁的检测数据保存得并不完善，往往只有最近一次的检测记录，无法利用已有方法将状态停留时间模型与马尔科夫链模型相结合，影响了对大量既有桥梁退化规律的研究。针对这个现状，本文提出了面向技术状况停留时间、基于检测结果发生概率的贝叶斯更新方法，该方法降低了对数据完备性的要求，可以利用不完备的历史桥检数据对既有桥梁各个技术状况的停留时间进行合理估计；然后结合马尔科夫链模型，给出了基于当前状况、停留时间和后续服役时间的状态转移概率表达式。通过数值算例验证了方法的准确性，并讨论了桥检数据量、数据完备性和数据误差对模型准确性的影响。最后结合某地区既有钢筋混凝土梁式桥上部结构的实测数据，对本文方法的应用进行了展示。

1   基于历史数据的技术状况退化模型

1.1   技术状况停留时间模型与贝叶斯更新

桥梁的退化状态通常用技术状况等级CR (Condition Rating)表示，例如在中国，桥梁技术状况等级由1~5(JTG/T H21−2011)^[14]依次变差，美国的桥梁状态等级为0~9越来越好(FHWA, 1995)^[15]。考虑到桥梁的设计、施工、维护等的差异，以及所处荷载条件及自然环境的不同，TC可描述成随机变量。已有研究表明技术状况的停留时间基本符合威布尔分布^[16]，其概率密度函数表示为：

式中：θ 为尺度参数；α 为形状参数。但为了研究方便，本文首先假设TC服从正态分布，其概率密度函数和分布函数分别为：

式中：$\mu$为正态分布的均值；$\sigma$为正态分布的标准差。这里假定TC服从正态分布，会在数据分析方面带来便利，但同时也会带来一定的分析误差，后面会对此展开讨论。

对一座既有桥梁而言，历史桥检数据主要包括桥梁的第i次检测的检测时间t_i和基于检测结果评定的技术状况等级I_i。采用贝叶斯更新的方法估计TC分布的参数时，可令D表示检测数据集、c表示模型的未知参数向量(这里c表示不同技术状况的停留时间TC的均值和标准差)、$P\left( c \right)$表示模型参数向量$\boldsymbol{c}$的先验分布，根据贝叶斯更新的原理，模型参数向量$\boldsymbol{c}$的后验分布由下式给出^[17]：

式中：$P\left( {\boldsymbol{D}} \right)$为桥梁检测数据集D出现的概率；$P\left( {{\boldsymbol{D}}\left| {\boldsymbol{c}} \right.} \right)$为给定模型参数条件下数据集D出现的概率，也称似然函数。显然对于给定的D，可以认为$P\left( {\boldsymbol{D}} \right)$是一个常数，因此有：

1.2   停留时间模型的更新方法

设${\mu _i}$和${\sigma _i}$为技术状况$i$的停留时间(TC_i)的均值和标准差。由于在实际分析时，TC_i不同，所用到的数据也来自于与之对应的退化到不同阶段的桥梁(例如：估计TC₁的概率分布用到的是处在技术状况1和2的桥梁的检测数据，估计TC₂用到的主要是处在技术状况2和3的桥梁的检测数据)，这些桥梁大部分是不同的，因此假设TC_i相互独立。

若某座桥的桥检数据只包括最近的一次检测记录，包括：检测时间t和桥梁所处的技术状态等级I。ZHANG等^[18]基于TC服从威布尔分布的假设，推导出了似然函数表达式如下：

式中，${f_k}\left( {{x_k}} \right)$为技术状况k的停留时间(TC_k)的概率密度函数。式(6)遵循了TC威布尔分布的假设，但它包含I重积分运算，计算上存在一定困难，不利于工程推广。

为了简化似然函数的计算，本文假设TC服从正态分布，则桥梁从刚刚建成到桥梁进入技术状况等级$i\left( {i > 1} \right)$所需的时间${T_i}$仍服从正态分布，其均值${\mu _{{T_i}}}$与标准差${\sigma _{{T_i}}}$分别为：

基于该假设，似然函数可表示为：

式中，Φ为标准正态分布的分布函数。与式(6)相比，正态分布的假设使似然函数的计算大为简化。

若某座桥保留了多次历史桥检数据，则需首先对这些数据进行处理：如果几个检测时间对应同一个技术状况等级I，则只保留最小和最大的检测时间，t_I,min和t_I,max(其余的检测记录不包含有用信息)；如果技术状况等级I只对应一个检测时间，则该检测时间记为t_I。多次桥检数据的似然函数为：

式中：A为有多个时间检测到了同一个技术状况的情况；B为某技术状况只对应一个检测时间。

以上是对应单个桥梁的似然函数。面对多个桥梁的检测数据时，可将所有桥梁的似然函数连乘获得总体似然函数，进行贝叶斯更新，即：根据桥梁检测数据完备性的差异，将式(8)或式(9)表示的似然函数代入式(4)，即可获得状态停留时间TC_i的后验分布参数值。在具体计算中，可以采用简单拒绝抽样的方法获得符合检测数据的样本^[17]。为了提高计算效率，本文采用基于Metropolis-Hasting算法的马尔科夫链蒙特卡罗模拟方法^[11]，从先验分布P(c)中抽取样本，通过不断地衍生新的样本并比较衍生前后似然函数的变化情况，逐渐得到可接收的样本，然后通过对接收样本的统计，得到各参数值的后验分布。

1.3   桥梁技术状况退化模型

利用1.2节的贝叶斯更新方法获得桥梁技术状况停留时间TC_i的后验分布参数值(${\mu _i}$和${\sigma _i}$)后，考虑到TC_i实际更符合威布尔分布，根据其后验分布的均值和标准差反算威布尔分布的参数值：尺度参数θ _i和形状参数α _i，反算所参照的表达式为：

式中，$\varGamma$为伽马函数，$\varGamma \left( x \right) = \int_0^\infty {{{\mathrm{e}}^{ - u}}{u^{x - 1}}{\mathrm{d}}u}$。将反算得到的θ _i和α _i代入式(1)即可得到TC_i的概率密度函数。

然后即可根据桥梁的当前状况和状态停留时间预测未来的退化情况。例如：桥梁目前的技术状况等级为I，且在该等级已经停留了a年，则下1年退化至 I+1 级的概率可表示为：

它是状态停留时间a的函数，随停留时间的增加，劣化到下一个状态的概率也随之增加。需要指出的是，式(13)假设桥梁的技术状况不会在1年内退化两个等级，因而忽略了一年中技术状况从I级退化到I+2级的可能性。

考虑到状态转移概率随停留时间的变化，可进一步推导出若干年后桥梁进入不同劣化等级的概率。依然以这座桥为例，它在技术状况等级I已经停留了a年，则b年之后，桥梁仍处在等级I的概率为：

该桥退化到技术状况等级I+1的概率分两种情况考虑：① 若I+1不是最差的技术状况等级，则：

② 若I+1是最差的技术状况等级，即${T_{I + 1}} = \infty$，则$a + b < {T_I} + {T_{I + 1}}$的条件一定满足，退化概率则为：

该桥退化到技术状况等级I+2的概率也分为两种情况：① 若I+2不是最差的技术状况等级，则：

② 若I+2是最差的技术状况等级，${T_{I + 2}} = \infty$，则$a + b < {T_I} + {T_{I + 1}} + {T_{I + 2}}$必然满足，退化概率则为：

计算该桥退化到技术状况等级I+3的概率。我国在对既有桥梁的技术状况进行评定时，当桥梁劣化到等级4(4类)，即认为其使用性能或安全性已不满足要求，需要进行维修加固，因此本研究把技术状况4作为最差的等级。即使当前的等级I = 1，I+3也是最后一个技术状况等级，退化到该等级的概率为：

利用上述公式，即可根据当前的技术状况等级和停留时间，对未来的劣化情况进行预测。图1给出了利用桥检数据更新技术状况停留时间模型，并对未来劣化情况进行预测的流程图。方法的具体应用将在后面的算例中展示。

图  1  本文方法流程图

Figure  1.  Flowchart of the proposed method

下载: 全尺寸图片 幻灯片

2   方法验证

2.1   模拟数据生成

为验证本文方法的准确性，首先假设一个已知的桥梁技术状况退化模型，用于模拟桥梁的性能退化。这样做一方面为本文的更新方法提供输入数据，另一方面通过“已知模型”与分析结果的对比，进行方法准确性的评价。

假设状态停留时间符合威布尔分布，其概率密度函数的尺度参数θ_i和形状参数α_i在表1中给出。

表  1  “已知的”技术状况停留时间模型的概率分布参数

Table  1.  Probability distribution parameters of assumed TC model

分布参数	技术状况1	技术状况2	技术状况3
形状参数αi	2.0	1.8	1.6
尺度参数θi	20	22	24

下载: 导出CSV 

| 显示表格

由于劣化到技术状况4的桥梁往往需要进行维修加固，导致实际工程中技术状况5的桥梁(5类桥)很少，对技术状况4的停留时间的评估也往往缺乏足够数据。因此在本文研究中只关注前3个技术状况。基于“已知的”的退化模型，随机生成500座桥梁的检测数据。首先生成这500座桥的桥龄，为1至80年的均匀分布；再根据表1的模型参数随机模拟出不同桥龄对应的技术状况等级，部分结果如图2所示。

图  2  生成的部分检测结果

Figure  2.  Simulated inspection results

下载: 全尺寸图片 幻灯片

利用本文的方法对这500座桥梁的数据进行分析，以验证方法的准确性。根据工程经验，状态1、2、3的停留时间大体在5年~30年，且离散性较大，因此假设TC_1、TC_2、TC₃先验分布的均值(μ₁、μ₂、μ₃)都是5年~30年间的均匀分布，标准差(σ₁、σ₂、σ₃)都是5年~15年间的均匀分布。按照第1.2节的方法得到接收样本，部分结果如图3所示。根据这些样本统计出参数(μ₁、σ₁、μ₂、σ₂、μ₃、σ₃)的平均值，分别为：μ₁ = 17.68、σ₁ = 9.94；μ₂ = 19.72、σ₂ = 10.47；μ₃ = 20.74、σ₃ = 11.17。

图  3  贝叶斯更新的部分接收样本

Figure  3.  Accepted samples in Bayesian updating

下载: 全尺寸图片 幻灯片

将上述三个退化状态的统计值代入式(11)，反算得到威布尔分布参数值：θ₁=19.70、α₁=1.82；θ₂=22.39、α₂=1.98；θ₃=23.24、α₃=1.88。根据上述参数值画出三个退化状态停留时间的概率分布曲线，如图4中的“正态假设”所示。为了验证方法的准确性，表1的“精确结果”也画在了图4中，可以看出对于这三个技术状况，基于正态假设的贝叶斯更新结果与精确结果都比较接近。相比而言，技术状况1的更新结果稍好于技术状况2，技术状况2的更新结果又好于技术状况3。这是由于桥梁结构总是由前一状态逐渐退化到后一状态，后一状态的数据里面也隐含着前一个状态的信息，并且前一状态停留时间的估计误差会影响到对后面状态停留时间的估计，因此状态越靠前，所获得的有用的信息量越多，数据的误差越小，更新效果也越好。以这500座桥为例，处在技术状况等级1~4的桥梁分别为128座、124座、152座、96座，即用来估计状况3的数据只有96条(进入状况4代表了状况3的结束，是估计状况3的停留时间所必不可少的)，而有助于评估状况1的数据有372条，数据量的差异是这两个技术状况评估准确性存在差异的主要原因。

图  4  状态停留时间的概率分布与精确解的对比

Figure  4.  Comparison of probability distribution of TC between updated results and exact results

下载: 全尺寸图片 幻灯片

为了讨论TC正态分布的假设对评估结果的影响，以式(6)替代式(8)作为似然函数，重新对上述数据进行贝叶斯更新计算，得到状态停留时间的概率分布，表示为图4中的“威布尔假设”。可以看出，由于精确结果是通过TC威布尔分布生成的，基于TC威布尔假设的更新结果几乎与精确结果重合，二者的微小差异是由数据量有所不足导致的。因此，相比于TC威布尔分布的假设，正态分布的假设会带来一定的误差，如图4所示，该误差随着数据量的增加而减小。当用于评估某技术状况的桥梁数量在100座以上时，利用本文方法可对该技术状况的停留时间进行较为准确的估计。

3   更新效果讨论

为给工程应用提供参考，本节讨论数据量、数据完备性和数据准确性对本文方法计算精度的影响。

3.1   数据量的影响

前面的例子用了500座桥的检测数据。这里将桥梁数量分别增加到了1000座和减少到了200座，利用本文方法对状态停留时间进行估计，其概率分布参数的更新结果见表2，可以看出随检测数据量的增加，得到的状态停留时间模型与精确解也愈发接近。

表  2  不同数据量下的参数更新结果

Table  2.  Updated parameters under different amount of data

数据量	技术状况1	技术状况2	技术状况3
200	α =1.80, θ =19.6	α =1.91, θ =21.1	α =2.25, θ =22.6
500	α =1.82, θ =19.7	α =1.98, θ =22.4	α =2.04, θ =22.9
1000	α =1.88, θ =19.8	α =1.96, θ =22.3	α =1.88, θ =23.5
精确解	α =2.00, θ =20.0	α =1.80, θ =22.0	α =1.60, θ =24.0

下载: 导出CSV 

| 显示表格

图5、图6分别给出了利用1000座和200座桥检数据更新得到的状态停留时间的概率分布情况。这里只给出了状态1和状态3，状态2结果的准确性介乎其间。由图5可以看出，当桥梁数量由500座增加到1000座后(处在技术状况1~4的桥梁分别为241座、264座、288座、207座)，对技术状况3停留时间的估计有了明显改善，但对状况1的估计仅略有改善，表明当桥检数据达到一定数量后，本文方法中状态停留时间的估计误差主要来自于似然函数表达式中用正态分布代替威布尔分布的假设，数据量的影响不大。可是当桥梁数量由500座减小到200座后(处在技术状况1~4的桥梁分别为42座、55座、62座、41座)，由图6可以看出，数据量的影响变得显著，状态1停留时间的估计精度相比500座和1000座的情况都有所降低，而状态3的估计则误差显著。综合来看，要获得对某个技术状况的停留时间的准确估计，用于评估该状况的桥梁数量应在100座以上。

图  5  1000座桥检数据的更新结果

Figure  5.  Updated results using inspection data of 1000 bridges

下载: 全尺寸图片 幻灯片

图  6  200座桥检数据的更新结果

Figure  6.  Updated results using inspection data of 200 bridges

下载: 全尺寸图片 幻灯片

3.2   数据完备性的影响

数据完备性表征桥梁检测历史是否完全、数据保留是否完整。假设这500座桥梁每两年检测一次，历史记录保存完整。这样，每一座桥梁的数据除了最后一次检测的时间和所处的技术状况等级之外，还包括进入和离开前面所有技术状况的大致时间。利用完备的桥检数据，代入似然函数表达式(式(9))，通过贝叶斯更新，得到状态1和3的停留时间的概率分布情况，如图7所示。

图  7  完备数据的更新结果

Figure  7.  Updated results using complete inspection data

下载: 全尺寸图片 幻灯片

相比图4中对技术状况1和3的估计，图7中的更新结果更加准确，这是由于完备数据较为准确地记录了进入和离开(从开始到现在)所有技术状况的时间，每一座桥的检测记录包含多个数据用于停留时间的更新，所以对各个技术状况，特别是技术状况3，更新效果也更加明显。

3.3   数据准确性的影响

工人基于检测结果对桥梁技术状况的评定往往带有一定的主观性，在面对同样的检测数据时，不同人员可能会给出不同的评定。下面讨论评定误差对停留时间模型的准确性的影响。假设三种评定准确率：90%、80%和70%。90%的准确率表示只有90%的桥梁被正确评定了技术状况等级，错误地评高或评低了一个等级的概率均为5%(如果技术状况的实际等级是1，则只能错误评为等级2，该概率是10%)；另外两个评定准确率也类似定义。

依然利用前面随机生成的500座桥梁的劣化数据(图2)，这里再根据评定准确率进一步随机生成桥检数据，然后进行状况停留时间参数的贝叶斯更新，结果如图8所示。可以看出，随评定准确率的降低，状况停留时间模型与实际情况的偏差也越来越明显。特别是对技术状况1的停留时间的估计，错误地把状况等级1评定为状况等级2，导致桥检数据低估了技术状况1的实际停留时间，带来系统性的偏差。当评定准确率为90%时，状况等级1的停留时间比实际情况少了3年左右；当评定准确率进一步降低到80%和70%，数据的系统性偏差导致的估计误差也愈发明显。对于技术状况3，随评定准确率的降低，结果误差也越来越显著，但当评定准确率在80%以上时，评估结果依然较为准确，且由于桥检数据不存在系统性偏差，停留时间整体偏大或偏小的现象并不明显，因此可以通过增大桥检数据量进一步改善评估结果的准确性。

图  8  桥检数据准确性对更新结果的影响

Figure  8.  Impact of inspection accuracy on updating accuracy

下载: 全尺寸图片 幻灯片

4   方法应用实例

本节从北方某城市2010年−2019年间约600份混凝土城市桥梁检测报告中选取185座钢筋混凝土梁式桥，分析其上部结构的退化规律。所选桥梁的桥龄及上部结构BCI (Bridge Condition Index)评分如图9所示。基于《公路桥梁技术状况评定标准(JTG/T H21−2011)》^[14]，其技术状况等级也标在了图中。

图  9  钢筋混凝土梁式桥的桥龄和上部结构BCI

Figure  9.  Age and BCI of the superstructure of concrete bridges

下载: 全尺寸图片 幻灯片

由于处在等级4的桥梁数据有限，不足以对技术状况等级3的停留时间做出准确评估，因此在本示例中把等级3作为最后一个技术状况，研究等级1和等级2的停留时间。

利用本文方法对图9中的数据进行分析，得到：θ₁=8.47、α₁=3.72；θ₂=9.36、α₂=3.21，可知技术状况2的停留时间(平均8.4年)稍大于技术状况1(平均7.6年)。代入式(13)和式(16)可得到停留在当前状况的时间与1年后或5年后退化到下一个状况的概率之间的关系，如图10所示。可以看出，在当前技术状况的停留时间越久，劣化到更低等级的概率越大。

进一步地，可对桥梁未来的退化状况进行预测。例如：桥梁A和B的上部结构都处在技术状况1，停留时间分别为0年(新桥)和5年，桥梁C和D 的上部结构都处在技术状况2，停留时间也分别为0年和5年。随后续服役时间的延长，桥梁处在不同技术状况的概率分布如图11所示。

图  10  劣化到下一个等级的概率与当前状态停留时间的关系

Figure  10.  Relationship between degradation probability and TC

下载: 全尺寸图片 幻灯片

与新桥A相比，桥梁B已有5年桥龄，虽然目前仍然处在技术状况1级，但未来劣化情况明显比新桥A更加严重。根据预测结果(图11(a)和图11(b))，9年之后，桥梁A处在技术状况1~3的概率分别为0.286、0.689、0.025；桥梁B停留在技术状况1的概率降低到了接近于0(0.002)，劣化到技术状况3的概率则增加到了0.232。桥梁C由于刚刚劣化到技术状况2，根据预测结果(图11(c)和图11(d))，9年之后它停留在该状态的概率为0.414；而桥梁D由于5年前已劣化到了技术状况2，9年之后依然保持在该状态的概率仅为0.030。

图  11  未来桥梁处在不同状况的概率与继续服役时间的关系

Figure  11.  Changing of the probabilities of bridges in different conditions with the service time

下载: 全尺寸图片 幻灯片

5   结论

本文提出了利用桥梁检测数据对未来技术状况劣化情况进行预测的方法。该方法基于对状态停留时间的估计，既可以使用记录完整的，也可以使用单次记录的桥检数据，弥补了已有方法无法使用不完整数据的缺陷。通过数值算例验证了方法的准确性，讨论了数据量、数据完备性和数据准确性的影响，并将该方法应用于某城市混凝土桥梁上部结构的技术状况劣化预测。研究表明：

(1)当用于评估某技术状况的桥梁数量在100座以上时，本文方法可对该技术状况的停留时间进行较为准确的估计。相比于单次桥检数据，完备桥检数据的更新准确性更好。

(2)技术状况的评定误差对退化模型的准确性有重要影响，特别是对技术状况1，评定误差会造成桥检数据的系统性偏差，导致对停留时间的明显低估；而对其它技术状况，当评定准确率在80%以上时，评估结果较为准确，且可通过增大桥检数据量进一步提高准确性。

(3)对实际混凝土桥上部结构的退化预测，发现技术状况1、2的停留时间平均为7.6和8.4年；与新桥相比，服役5年的桥梁未来保持在技术状况1的概率明显减小，劣化到技术状况3的概率明显增加。由于数据量有限(185座)，该预测存在一定误差，本文方法需要增加更多桥检数据以提高退化模型的准确性。