基于Copula理论的桥梁可靠度计算

罗 辉， 李 彤， 熊凯文

(1. 华中科技大学 土木工程与力学学院， 湖北 武汉 430074； 2. 深圳市华阳国际建筑产业化有限公司， 广东 深圳 518048)

可靠性是对结构安全性、适用性及耐久性的总称，对可靠性进行度量的指标称之为可靠度[1]。结构可靠度问题的研究始于20世纪40年代，前苏联学者首次提出一次二阶矩的基本概念后，各国结构可靠性研究在此基础上蓬勃发展，取得了丰硕的成果。当前，对于桥梁结构，构件可靠度的研究已经相对成熟，而体系可靠度的研究仍在发展中[2]。实际工程中，桥梁体系失效路径繁多，难以直接采用功能函数表达；采用串并联模型法可以实现从结构到构件的演绎，但难以衡量失效模式间的相关系数[3]，进而影响体系失效概率求解的准确性。针对串并联模型的近似估计，现有经典理论为界限估计法[2]；Cornell提出一阶区间理论，仅考虑相关程度的两个极值状况，即完全相关与相互独立，对串联、并联体系分别给出计算公式，称之为宽界限公式；Ditlevsen[4]提出窄界限法，考虑了两个失效模式同时失效的联合失效概率，但未能有效估计联合失效概率，使得该方法应用困难。

Copula理论最早由Sklar[5]提出，本质上是通过一个函数沟通多元联合分布与边缘概率分布，可以描述边缘分布间的相关关系，目前主要运用于金融、数理等领域[6,7]。在结构可靠性方面，Copula理论的应用更集中于岩土工程，许晓亮等[8]对均质边坡提出了Copula理论下基于G-line失效域的可靠性分析方法，并对不同Copula函数下的失效概率进行了比较；黄达等[9]利用Copula理论在建立多变量联合分布函数上的优势，构造了拟合粗粒土临界水力比降、孔隙比、级配不均匀系数以及曲率系数共4个参数间相关关系的最优Copula函数，并应用于估计粗粒土渗透破坏临界水力比降的值。

本文引入Copula理论，用于对构件失效模式间相关关系的度量，并构建联合分布函数从而得到联合失效概率。采用Copula理论，相对于局限于二元相关关系窄界限理论，具有更好的适用性；同时不同于现有大部分理论仅考虑线性相关关系，Copula理论既可以衡量线性相关关系，也可以衡量非线性相关关系。在构件可靠度计算的基础上，得到构件功能函数表达式；结合串并联模型法，将桥梁体系失效概率分解为构件失效间的串-并联关系，代入Copula计算可得到考虑了相关性的桥梁体系失效概率。

1 研究方法

1.1 Copula理论

Copula函数[10]是把随机向量X1,X2,…XN的联合分布函数F(x1,x2,…xN)与各自的边缘分布函数FX1(x1),FX2(x2),…FXN(xN)相连接的连接函数。Copula函数满足Sklar定理：

F(x1,x2,…xN)=C(FX1(x1),FX2(x2),…，FXN(xN))

(1)

若FX1(x1),FX2(x2),…,FXN(xN)是连续函数，则C(u1,u2,…,uN)唯一确定；反之，若FX1(x1),FX2(x2),…,FXN(xN)为一元分布函数，C(u1,u2,…,uN)是一个Copula函数，则由式(1)可确定具有边缘分布FX1(x1),FX2(x2),…,FXN(xN)的N元联合分布函数FX1(x1),FX2(x2),…FXN(xN)。

不同类型的Copula函数在描述变量间的相关关系上具有不同的特点。常见的Copula函数有高斯Copula，t Copula及阿基米德Copula函数，将上述常用Copula函数的二元分布函数C(u,v)及其参数范围汇总于表1[11]。由Copula函数参数可计算Copula中所包含的边缘分布相关关系，可采用秩相关系数、Spearman相关系数及描述尾部相关性的上尾、下尾相关系数进行衡量，计算公式如表2所示[11]。

表1 常用Copula函数

表2 Copula参数与相关系数解析表达式

1.2 体系可靠度

串并联模型法是主要的体系可靠度计算方法之一。实际上，大部分结构都是串、并联并存的组合模型。串并联模型法较全面地考虑了整体结构失效和单个构件失效之间的关系，能较好地解决构件到结构的归纳分析问题[3]。

串并联模型常见的近似计算方法主要包括宽界限法和窄界限法[2]。宽界限法仅考虑体系中失效模式相关性的两种极端状况，一是所有失效模式完全相关，二是所有失效模式完全独立；记体系失效概率为Pf，第i个模式失效概率为Pfi，可将体系失效概率的最大界限估计表示为：

(1)

Ditlevsen对此提出了改进，给出了窄界限公式：

(3)

式中：Pfij为模式i，j同时失效的概率。

本文将Copula函数与串并联模型法结合，将Copula函数在描述多元变量相关性方面的优良性能，及其构建联合概率分布函数的便捷性运用到体系失效概率的计算中，推导结合Copula函数的串联、并联及混合体系中失效概率的表达式如下所示：

(1)串联体系

对于两个失效模式串联的二元串联体系，任一失效模式失效则体系失效，记失效功能函数为Z1，Z2，其边缘概率分布函数分别为F1(z1)，F2(z2)，联合概率分布函数为F(z1,z2)，则由概率论随机事件计算的加法公式及式(1)可知，二元串联体系失效概率Pf为：

Pf=P(Z1≤0∪Z2≤0)

=P(Z1≤0)+P(Z2≤0)-P(Z1≤0,Z2≤0)

=F1(0)+F2(0)-F(0,0)

=F1(0)+F2(0)-C(F1(0),F2(0))

(4)

同理，对于N元串联体系，记第i个失效功能函数Zi的边缘概率分布函数为Fi(Z)，构造n维联合概率分布函数的Copula函数，记为Cn，则体系失效概率Pf为：

(-1)N-1P(Z1≤0,Z2≤0,…,ZN≤0)

(-1)N-1CN(F1(0),F2(0),…,FN(0))

(5)

(2)并联体系

两个失效模式均失效则体系失效，故二元并联体系失效概率Pf为：

Pf=P(Z1≤0,Z2≤0)=F(0,0)=C(F1(0),F2(0))

(6)

同理，多元并联体系有：

=F(0,0,…,0)=C(F1(0),F2(0),…,FN(0))

(7)

(3)混联体系

工程实际中结构体系不一定能划分能单纯串联或是并联模型，很可能存在串-并联均有的混合模型。将m个失效模式先串联为1个子系统而后n个子系统并联成体系的模型称为串-并联体系；将m个失效模式先并联为1个子系统而后n个子系串联成体系的模型称为并-串联模型。将第i个子系统的第j个失效模式记为Zij，不考虑子系统间的相关性，则对于串-并联体系有：

(8)

对于并-串联体系有：

(9)

结合式(5)，(7)即可进行求解。

2 Copula建模

2.1 基于MATLAB的Copula建模

本文在MATLAB平台编程构建Copula模型，由变量的边缘分布获得二元及多元联合分布。按以下步骤构建模型：

(1)获取样本数据。可以实际调查数据，或通过蒙特卡罗法计算，得到需要构建联合分布的边缘分布函数的样本数据。

(2)确定边缘分布。根据样本数据确定边缘分布可采用两种方法：其一是参数法，即假设随机变量符合常见的分布类型，通过分布拟合工具箱估计参数，而后对假设进行检验；其二是非参法，该方法无需假定样本服从某个分布，而是依据样本点求总体分布的近似，常用的有经验分布函数法与核密度估计法。

(3)选取Copula函数。通常由二元频率直方图的形状，在常见Copula函数选取概率密度图近似的Copula函数。

(4)参数估计。选取的Copula函数包含有未知参数，需对其进行估计。常用的方法有最大似然估计、分步估计及半参数估计[11]。

本文选取半参数估计，以边缘分布函数的近似估计替代边缘分布函数。选取Copula函数为C(u,v,α)，其密度函数为c(u,v,α)，且有c(u,v,α)=C(u,v,α)/(∂u∂v)，则按式(10)对α进行估计。

(10)

2.2 Copula评价准则

以上述步骤建立Copula模型后需进行模型评价，以便从一组Copula函数中选择拟合效果最好的Copula函数。常用评价准则如下：

(1)AIC信息准则

由日本统计学家赤池弘次提出并完善，他结合K-L距离与最大似然法，定义：

AIC=-2ln(最大似然函)+2k

(11)

式中：k为独立参数维度。AIC越小则模型越优。

(2)BIC信息准则

贝叶斯信息准则，在AIC的基础上考虑了样本数量，防止样本量过大导致的模型复杂化。

BIC=-2ln(最大似然函)+kln(N)

(12)

式中：N为样本容量。BIC越小则模型越优。

(3)平方欧氏距离[12]

平方欧氏距离表示理论累计分布函数值和经验累计分布函数值在原始样本点处偏差的平方和，定义为：

(13)

AIC，BIC信息准则和平方欧氏距离方法依赖于参数估计，在小样本数据条件下，参数估计不准确可能导致错误的识别结果，在样本量较多的情况下这三种方法原理简单、稳定性好，都可以采用。

本文算例中引入最小二乘向量机对有限元计算数据进行拓展，符合大样本数据条件，AIC，BIC信息准则和平方欧氏距离方法均适用，本文选择平方欧氏距离方法。

3 算例分析

3.1 桁架算例

3.1.1 桁架体系失效概率计算

某单层钢桁架结构如图1所示，仅考虑杆件的拉、压破坏。各杆件的横截面积为确定性参数，另将随机变量分别列入表3。

图1 某单层钢桁架/m

表3 体系随机变量

由结构受力分析可知，该结构为静定结构即任一杆件破坏都将造成体系失效。其中杆受拉，杆，受压，杆不受力，故该体系为3个失效模式构成的串联体系，失效功能函数可写作：

Zi=Ai[σc]-Gi

(14)

式中：Zi表示编号为i的杆件的失效功能函数；Ai，Gi分别表示编号为i的杆件的横截面积和内力；[σc]为容许应力。

结合式(5)，基于Copula的体系失效概率可按下式计算：

C(F1(0),F3(0))-C(F2(0),F3(0))+

C(F1(0),F2(0),F3(0))

(15)

采用蒙特卡罗法计算构件失效概率，并得到失效功能函数的样本数据；据前文中所述流程构建Copula模型，可得到由Copula函数计算得到的联合失效概率，而后由式(15)可得到体系失效概率，如表4所示。

表4 联合失效概率

3.1.2 与现有界限估计的比较

现有界限估计也需要求解联合失效概率。Thoft-Christensen和Murotsu[13]提出了针对正相关的近似计算公式；Feng[14]提出了考虑失效模式构成的超平面夹角近似估计；董聪等[15]提出考虑两个失效模式的失效概率值存在一定差异状况下的近似计算公式。但仍过于粗糙，不能带入窄界限估计中进行计算。

董聪和夏人伟[15]提出采用数值积分的方式求解联合概率，该方法要求功能函数正态分布Xi，Xj，然后利用二维正态密度函数可显示表达联合分布，从而求解联合失效概率值，通常认为该方法为近似精确解。相关系数采用线性相关系数ρij和二维联合分布Φ(xi,xj;ρij)，φ(ti,tj;ρij)分别按式(16)(17)计算。

(16)

(17)

该方法计算得到的二元联合失效概率值与Copula函数值结果基本一致，差值在0.05%～6%之间，差距非常微弱，如图2所示。

图2 不同计算方法下二元联合失效概率值

将本文计算得到的体系失效概率值与宽界限法及窄界限法进行对比如表5所示。

表5 不同估计方法下体系失效概率值 ×10-2

由本例的计算可见，由宽界限至窄界限，虽然区间范围变化巨大，但基于Copula模型得到的体系可靠度点估计值始终在区间范围内，说明本文采用的计算方式具有较高的可信度。同时，Copula建模的方式不限制衡量线性相关关系，采用非参法亦不限制功能函数所服从的分布类型，且建模计算的复杂程度亦不受边缘分布维度的影响，故具有较高的适用性。

3.2 连续梁桥算例

以武汉至孝感高速公路上的某宽幅预应力空心板五跨连续梁桥为例，主梁由预制预应力混凝土空心板和现浇混凝土桥面板组合而成，先简支后连续，每半桥宽主梁含8块空心板。半桥宽横截面如图3所示。

图3 半跨宽桥梁横幅布置/cm

本文采用ANSYS进行有限元建模对连续梁桥的受力状态进行计算，采用MATLAB计算结构抗力；由蒙特卡罗法计算构件的失效概率，仅考虑空心板的抗弯与抗剪失效；并引入最小二乘向量机对有限元计算数据进行拓展，以便满足计算的样本需求。流程如图4所示。

图4 可靠度计算流程

选择表6[16,17]中参数作为不确定性参数进行计算。参考文献[18]，任意相邻两块空心板失效则体系失效，则以边跨为例，可将两块相邻空心板视作1个子系统，即1个子系统由2块空心板并联，一片主梁由7个子系统串联，形成图5所示串并联模型。

表6 随机变量分布特征

图5 串-并联模型

本文仅考虑子系统内部相关性，在假设子系统之间相互独立的前提下，按考虑相关性程度不同，分为以下3种计算模式。

(1)完全不相关：假定各个空心板失效相互独立，假定单个空心板的两种失效模式也相互独立。

(2)仅考虑单块板的不同失效模式相关：假定各个空心板的失效相互独立，考虑单块板不同失效模式间的相关性。

(3)考虑子系统内相关性：考虑单块板不同失效模式间的相关性，也考虑相邻空心板的抗弯抗剪失效模式间存在的相关性，且考虑两块板同时失效存在的相关性；即考虑子系统内部任意一块空心板任意一种失效模式之间均存在相关性。

以上考虑相关性的失效模式均需要采用Copula模型衡量相关关系，并计算联合失效概率，而后按第1部分公式求解体系失效概率。以ZM，ZV分别表示抗弯、抗剪失效，数字表示空心板编号，则由Copula模型得到子系统1内的联合失效概率及相关关系，如表7所示。

表7 二元Copula联合分布结果

比较各个不同失效模式间的联合失效概率值，作出柱状图如图6。

图6 子系统1不同失效模式间的二元联合失效概率值

由图6可见，子系统内不同失效模式的联合失效概率存在巨大差异。其中不同空心板抗弯失效的联合失效概率(ZM_1&ZM_2)最大，其次是不同空心板抗剪失效的联合失效概率(ZV_1&ZV_2)，二者相差约21倍，二者边缘分布间的Kendall相关系数均在0.95左右，而边缘失效概率值则分别为(3.57×10-5,3.74×10-5)、(8.30×10-6,1.60×10-6)，可见在边缘分布相关程度相当的情况下，边缘分布失效概率值越大则联合失效概率值越大。

另一方面，将ZV_1&ZM_2的联合失效概率与ZV_1&ZV_2比较，前者边缘分布概率为(8.30×10-6, 3.74×10-5)，后者为(8.30×10-6,1.60×10-6)，发现前者的边缘失效概率远大于后者，而后者联合失效概率却远大于前者，比较二者的Kendall相关系数，前者为0.22，后者为0.94，可见边缘分布的相关性越强，则联合分布概率值越大。

将不同计算模式下各个子系统及边跨体系失效概率列入表8。

表8 不同计算模式下的体系失效概率

对比计算模式①与计算模式②中各个子系统的失效概率值，发现二者差值非常微弱，累计差值约0.03%。原因在于同一块空心板的抗弯失效与抗剪失效的联合失效概率值极小，而该值小又是由抗弯失效与抗剪失效的相关性很弱导致，故由此可知，体系中的弱相关关系可以忽略，当作相互独立进行处理对失效概率的计算影响不大。

对比计算模式②③可见，计算模式③得到的子系统失效概率约为计算模式②得到的子系统失效概率的2×104，失效概率大幅增长意味着可靠性大幅下降，结构存在安全隐患的可能性大幅提高。这一显著差异的来源在于，计算模式②将并联的空心板作为相互独立关系处理，将单块板的失效概率直接相乘，而计算模式③则是认为并联空心板存在相关性，通过联合概率分布求取两块板同时失效的概率；通过表7的计算结果可知，并联空心板之间确实存在强相关关系，故而联合失效概率值远大于二者乘积。由此可见，强相关关系对体系失效概率值的影响十分显著，在计算中应予以关注，不可忽略。

另外对比计算模式①②③对体系失效概率的计算公式可知，对于串联结构而言，不考虑相关性，将构件及其失效模式当作相互独立处理，得到的体系失效概率偏大，可靠性偏小，计算结果偏保守；对于并联结构而言，不考虑相关性则使得体系失效概率偏小，可靠度偏大，即计算结果较危险，若是并联结构的强相关关系被忽略，则很可能无法对结构存在安全隐患做出正确判断。

4 结 论

(1)本文利用Copula函数对体系可靠度进行计算，在衡量相关性时可以考虑到非线性相关，同时可以得到与数值积分精确度相近的联合失效概率，没有参数符合正态分布的限制，更贴合工程实际也更具适用性；

(2)本文利用Copula函数构建多元联合分布，结合串并联模型可将体系层面的失效概率降维到构件失效模式，避开体系失效模式难以显示表达的问题；

(3)本文将Copula理论运用于空心板连续梁桥体系可靠度计算发现，忽略弱相关关系对体系失效概率的影响不显著，但忽略强相关关系体系失效概率可能偏差若干数量级，使计算结果不可信赖，故在对体系进行失效概率计算时充分考虑体系内相关性是非常有必要的。