羅 輝， 李 彤， 熊凱文

(1. 華中科技大學(xué) 土木工程與力學(xué)學(xué)院， 湖北 武漢 430074； 2. 深圳市華陽國際建筑產(chǎn)業(yè)化有限公司， 廣東 深圳 518048)

可靠性是對結(jié)構(gòu)安全性、適用性及耐久性的總稱，對可靠性進行度量的指標(biāo)稱之為可靠度[1]。結(jié)構(gòu)可靠度問題的研究始于20世紀(jì)40年代，前蘇聯(lián)學(xué)者首次提出一次二階矩的基本概念后，各國結(jié)構(gòu)可靠性研究在此基礎(chǔ)上蓬勃發(fā)展，取得了豐碩的成果。當(dāng)前，對于橋梁結(jié)構(gòu)，構(gòu)件可靠度的研究已經(jīng)相對成熟，而體系可靠度的研究仍在發(fā)展中[2]。實際工程中，橋梁體系失效路徑繁多，難以直接采用功能函數(shù)表達；采用串并聯(lián)模型法可以實現(xiàn)從結(jié)構(gòu)到構(gòu)件的演繹，但難以衡量失效模式間的相關(guān)系數(shù)[3]，進而影響體系失效概率求解的準(zhǔn)確性。針對串并聯(lián)模型的近似估計，現(xiàn)有經(jīng)典理論為界限估計法[2]；Cornell提出一階區(qū)間理論，僅考慮相關(guān)程度的兩個極值狀況，即完全相關(guān)與相互獨立，對串聯(lián)、并聯(lián)體系分別給出計算公式，稱之為寬界限公式；Ditlevsen[4]提出窄界限法，考慮了兩個失效模式同時失效的聯(lián)合失效概率，但未能有效估計聯(lián)合失效概率，使得該方法應(yīng)用困難。

Copula理論最早由Sklar[5]提出，本質(zhì)上是通過一個函數(shù)溝通多元聯(lián)合分布與邊緣概率分布，可以描述邊緣分布間的相關(guān)關(guān)系，目前主要運用于金融、數(shù)理等領(lǐng)域[6,7]。在結(jié)構(gòu)可靠性方面，Copula理論的應(yīng)用更集中于巖土工程，許曉亮等[8]對均質(zhì)邊坡提出了Copula理論下基于G-line失效域的可靠性分析方法，并對不同Copula函數(shù)下的失效概率進行了比較；黃達等[9]利用Copula理論在建立多變量聯(lián)合分布函數(shù)上的優(yōu)勢，構(gòu)造了擬合粗粒土臨界水力比降、孔隙比、級配不均勻系數(shù)以及曲率系數(shù)共4個參數(shù)間相關(guān)關(guān)系的最優(yōu)Copula函數(shù)，并應(yīng)用于估計粗粒土滲透破壞臨界水力比降的值。

本文引入Copula理論，用于對構(gòu)件失效模式間相關(guān)關(guān)系的度量，并構(gòu)建聯(lián)合分布函數(shù)從而得到聯(lián)合失效概率。采用Copula理論，相對于局限于二元相關(guān)關(guān)系窄界限理論，具有更好的適用性；同時不同于現(xiàn)有大部分理論僅考慮線性相關(guān)關(guān)系，Copula理論既可以衡量線性相關(guān)關(guān)系，也可以衡量非線性相關(guān)關(guān)系。在構(gòu)件可靠度計算的基礎(chǔ)上，得到構(gòu)件功能函數(shù)表達式；結(jié)合串并聯(lián)模型法，將橋梁體系失效概率分解為構(gòu)件失效間的串-并聯(lián)關(guān)系，代入Copula計算可得到考慮了相關(guān)性的橋梁體系失效概率。

1 研究方法

1.1 Copula理論

Copula函數(shù)[10]是把隨機向量X1,X2,…XN的聯(lián)合分布函數(shù)F(x1,x2,…xN)與各自的邊緣分布函數(shù)FX1(x1),FX2(x2),…FXN(xN)相連接的連接函數(shù)。Copula函數(shù)滿足Sklar定理：

F(x1,x2,…xN)=C(FX1(x1),FX2(x2),…，F(xiàn)XN(xN))

(1)

若FX1(x1),FX2(x2),…,FXN(xN)是連續(xù)函數(shù)，則C(u1,u2,…,uN)唯一確定；反之，若FX1(x1),FX2(x2),…,FXN(xN)為一元分布函數(shù)，C(u1,u2,…,uN)是一個Copula函數(shù)，則由式(1)可確定具有邊緣分布FX1(x1),FX2(x2),…,FXN(xN)的N元聯(lián)合分布函數(shù)FX1(x1),FX2(x2),…FXN(xN)。

不同類型的Copula函數(shù)在描述變量間的相關(guān)關(guān)系上具有不同的特點。常見的Copula函數(shù)有高斯Copula，t Copula及阿基米德Copula函數(shù)，將上述常用Copula函數(shù)的二元分布函數(shù)C(u,v)及其參數(shù)范圍匯總于表1[11]。由Copula函數(shù)參數(shù)可計算Copula中所包含的邊緣分布相關(guān)關(guān)系，可采用秩相關(guān)系數(shù)、Spearman相關(guān)系數(shù)及描述尾部相關(guān)性的上尾、下尾相關(guān)系數(shù)進行衡量，計算公式如表2所示[11]。

表1 常用Copula函數(shù)

表2 Copula參數(shù)與相關(guān)系數(shù)解析表達式

1.2 體系可靠度

串并聯(lián)模型法是主要的體系可靠度計算方法之一。實際上，大部分結(jié)構(gòu)都是串、并聯(lián)并存的組合模型。串并聯(lián)模型法較全面地考慮了整體結(jié)構(gòu)失效和單個構(gòu)件失效之間的關(guān)系，能較好地解決構(gòu)件到結(jié)構(gòu)的歸納分析問題[3]。

串并聯(lián)模型常見的近似計算方法主要包括寬界限法和窄界限法[2]。寬界限法僅考慮體系中失效模式相關(guān)性的兩種極端狀況，一是所有失效模式完全相關(guān)，二是所有失效模式完全獨立；記體系失效概率為Pf，第i個模式失效概率為Pfi，可將體系失效概率的最大界限估計表示為：

(1)

Ditlevsen對此提出了改進，給出了窄界限公式：

(3)

式中：Pfij為模式i，j同時失效的概率。

本文將Copula函數(shù)與串并聯(lián)模型法結(jié)合，將Copula函數(shù)在描述多元變量相關(guān)性方面的優(yōu)良性能，及其構(gòu)建聯(lián)合概率分布函數(shù)的便捷性運用到體系失效概率的計算中，推導(dǎo)結(jié)合Copula函數(shù)的串聯(lián)、并聯(lián)及混合體系中失效概率的表達式如下所示：

(1)串聯(lián)體系

對于兩個失效模式串聯(lián)的二元串聯(lián)體系，任一失效模式失效則體系失效，記失效功能函數(shù)為Z1，Z2，其邊緣概率分布函數(shù)分別為F1(z1)，F(xiàn)2(z2)，聯(lián)合概率分布函數(shù)為F(z1,z2)，則由概率論隨機事件計算的加法公式及式(1)可知，二元串聯(lián)體系失效概率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元串聯(lián)體系，記第i個失效功能函數(shù)Zi的邊緣概率分布函數(shù)為Fi(Z)，構(gòu)造n維聯(lián)合概率分布函數(shù)的Copula函數(shù)，記為Cn，則體系失效概率Pf為：

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

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

(5)

(2)并聯(lián)體系

兩個失效模式均失效則體系失效，故二元并聯(lián)體系失效概率Pf為：

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

(6)

同理，多元并聯(lián)體系有：

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

(7)

(3)混聯(lián)體系

工程實際中結(jié)構(gòu)體系不一定能劃分能單純串聯(lián)或是并聯(lián)模型，很可能存在串-并聯(lián)均有的混合模型。將m個失效模式先串聯(lián)為1個子系統(tǒng)而后n個子系統(tǒng)并聯(lián)成體系的模型稱為串-并聯(lián)體系；將m個失效模式先并聯(lián)為1個子系統(tǒng)而后n個子系串聯(lián)成體系的模型稱為并-串聯(lián)模型。將第i個子系統(tǒng)的第j個失效模式記為Zij，不考慮子系統(tǒng)間的相關(guān)性，則對于串-并聯(lián)體系有：

(8)

對于并-串聯(lián)體系有：

(9)

結(jié)合式(5)，(7)即可進行求解。

2 Copula建模

2.1 基于MATLAB的Copula建模

本文在MATLAB平臺編程構(gòu)建Copula模型，由變量的邊緣分布獲得二元及多元聯(lián)合分布。按以下步驟構(gòu)建模型：

(1)獲取樣本數(shù)據(jù)?？梢詫嶋H調(diào)查數(shù)據(jù)，或通過蒙特卡羅法計算，得到需要構(gòu)建聯(lián)合分布的邊緣分布函數(shù)的樣本數(shù)據(jù)。

(2)確定邊緣分布。根據(jù)樣本數(shù)據(jù)確定邊緣分布可采用兩種方法：其一是參數(shù)法，即假設(shè)隨機變量符合常見的分布類型，通過分布擬合工具箱估計參數(shù)，而后對假設(shè)進行檢驗；其二是非參法，該方法無需假定樣本服從某個分布，而是依據(jù)樣本點求總體分布的近似，常用的有經(jīng)驗分布函數(shù)法與核密度估計法。

(3)選取Copula函數(shù)。通常由二元頻率直方圖的形狀，在常見Copula函數(shù)選取概率密度圖近似的Copula函數(shù)。

(4)參數(shù)估計。選取的Copula函數(shù)包含有未知參數(shù)，需對其進行估計。常用的方法有最大似然估計、分步估計及半?yún)?shù)估計[11]。

本文選取半?yún)?shù)估計，以邊緣分布函數(shù)的近似估計替代邊緣分布函數(shù)。選取Copula函數(shù)為C(u,v,α)，其密度函數(shù)為c(u,v,α)，且有c(u,v,α)=C(u,v,α)/(?u?v)，則按式(10)對α進行估計。

(10)

2.2 Copula評價準(zhǔn)則

以上述步驟建立Copula模型后需進行模型評價，以便從一組Copula函數(shù)中選擇擬合效果最好的Copula函數(shù)。常用評價準(zhǔn)則如下：

(1)AIC信息準(zhǔn)則

由日本統(tǒng)計學(xué)家赤池弘次提出并完善，他結(jié)合K-L距離與最大似然法，定義：

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

(11)

式中：k為獨立參數(shù)維度。AIC越小則模型越優(yōu)。

(2)BIC信息準(zhǔn)則

貝葉斯信息準(zhǔn)則，在AIC的基礎(chǔ)上考慮了樣本數(shù)量，防止樣本量過大導(dǎo)致的模型復(fù)雜化。

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

(12)

式中：N為樣本容量。BIC越小則模型越優(yōu)。

(3)平方歐氏距離[12]

平方歐氏距離表示理論累計分布函數(shù)值和經(jīng)驗累計分布函數(shù)值在原始樣本點處偏差的平方和，定義為：

(13)

AIC，BIC信息準(zhǔn)則和平方歐氏距離方法依賴于參數(shù)估計，在小樣本數(shù)據(jù)條件下，參數(shù)估計不準(zhǔn)確可能導(dǎo)致錯誤的識別結(jié)果，在樣本量較多的情況下這三種方法原理簡單、穩(wěn)定性好，都可以采用。

本文算例中引入最小二乘向量機對有限元計算數(shù)據(jù)進行拓展，符合大樣本數(shù)據(jù)條件，AIC，BIC信息準(zhǔn)則和平方歐氏距離方法均適用，本文選擇平方歐氏距離方法。

3 算例分析

3.1 桁架算例

3.1.1 桁架體系失效概率計算

某單層鋼桁架結(jié)構(gòu)如圖1所示，僅考慮桿件的拉、壓破壞。各桿件的橫截面積為確定性參數(shù)，另將隨機變量分別列入表3。

圖1 某單層鋼桁架/m

表3 體系隨機變量

由結(jié)構(gòu)受力分析可知，該結(jié)構(gòu)為靜定結(jié)構(gòu)即任一桿件破壞都將造成體系失效。其中桿受拉，桿，受壓，桿不受力，故該體系為3個失效模式構(gòu)成的串聯(lián)體系，失效功能函數(shù)可寫作：

Zi=Ai[σc]-Gi

(14)

式中：Zi表示編號為i的桿件的失效功能函數(shù)；Ai，Gi分別表示編號為i的桿件的橫截面積和內(nèi)力；[σc]為容許應(yīng)力。

結(jié)合式(5)，基于Copula的體系失效概率可按下式計算：

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

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

(15)

采用蒙特卡羅法計算構(gòu)件失效概率，并得到失效功能函數(shù)的樣本數(shù)據(jù)；據(jù)前文中所述流程構(gòu)建Copula模型，可得到由Copula函數(shù)計算得到的聯(lián)合失效概率，而后由式(15)可得到體系失效概率，如表4所示。

表4 聯(lián)合失效概率

3.1.2 與現(xiàn)有界限估計的比較

現(xiàn)有界限估計也需要求解聯(lián)合失效概率。Thoft-Christensen和Murotsu[13]提出了針對正相關(guān)的近似計算公式；Feng[14]提出了考慮失效模式構(gòu)成的超平面夾角近似估計；董聰?shù)萚15]提出考慮兩個失效模式的失效概率值存在一定差異狀況下的近似計算公式。但仍過于粗糙，不能帶入窄界限估計中進行計算。

董聰和夏人偉[15]提出采用數(shù)值積分的方式求解聯(lián)合概率，該方法要求功能函數(shù)正態(tài)分布Xi，Xj，然后利用二維正態(tài)密度函數(shù)可顯示表達聯(lián)合分布，從而求解聯(lián)合失效概率值，通常認為該方法為近似精確解。相關(guān)系數(shù)采用線性相關(guān)系數(shù)ρij和二維聯(lián)合分布Φ(xi,xj;ρij)，φ(ti,tj;ρij)分別按式(16)(17)計算。

(16)

(17)

該方法計算得到的二元聯(lián)合失效概率值與Copula函數(shù)值結(jié)果基本一致，差值在0.05%～6%之間，差距非常微弱，如圖2所示。

圖2 不同計算方法下二元聯(lián)合失效概率值

將本文計算得到的體系失效概率值與寬界限法及窄界限法進行對比如表5所示。

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

由本例的計算可見，由寬界限至窄界限，雖然區(qū)間范圍變化巨大，但基于Copula模型得到的體系可靠度點估計值始終在區(qū)間范圍內(nèi)，說明本文采用的計算方式具有較高的可信度。同時，Copula建模的方式不限制衡量線性相關(guān)關(guān)系，采用非參法亦不限制功能函數(shù)所服從的分布類型，且建模計算的復(fù)雜程度亦不受邊緣分布維度的影響，故具有較高的適用性。

3.2 連續(xù)梁橋算例

以武漢至孝感高速公路上的某寬幅預(yù)應(yīng)力空心板五跨連續(xù)梁橋為例，主梁由預(yù)制預(yù)應(yīng)力混凝土空心板和現(xiàn)澆混凝土橋面板組合而成，先簡支后連續(xù)，每半橋?qū)捴髁汉?塊空心板。半橋?qū)挋M截面如圖3所示。

圖3 半跨寬橋梁橫幅布置/cm

本文采用ANSYS進行有限元建模對連續(xù)梁橋的受力狀態(tài)進行計算，采用MATLAB計算結(jié)構(gòu)抗力；由蒙特卡羅法計算構(gòu)件的失效概率，僅考慮空心板的抗彎與抗剪失效；并引入最小二乘向量機對有限元計算數(shù)據(jù)進行拓展，以便滿足計算的樣本需求。流程如圖4所示。

圖4 可靠度計算流程

選擇表6[16,17]中參數(shù)作為不確定性參數(shù)進行計算。參考文獻[18]，任意相鄰兩塊空心板失效則體系失效，則以邊跨為例，可將兩塊相鄰空心板視作1個子系統(tǒng)，即1個子系統(tǒng)由2塊空心板并聯(lián)，一片主梁由7個子系統(tǒng)串聯(lián)，形成圖5所示串并聯(lián)模型。

表6 隨機變量分布特征

圖5 串-并聯(lián)模型

本文僅考慮子系統(tǒng)內(nèi)部相關(guān)性，在假設(shè)子系統(tǒng)之間相互獨立的前提下，按考慮相關(guān)性程度不同，分為以下3種計算模式。

(1)完全不相關(guān)：假定各個空心板失效相互獨立，假定單個空心板的兩種失效模式也相互獨立。

(2)僅考慮單塊板的不同失效模式相關(guān)：假定各個空心板的失效相互獨立，考慮單塊板不同失效模式間的相關(guān)性。

(3)考慮子系統(tǒng)內(nèi)相關(guān)性：考慮單塊板不同失效模式間的相關(guān)性，也考慮相鄰空心板的抗彎抗剪失效模式間存在的相關(guān)性，且考慮兩塊板同時失效存在的相關(guān)性；即考慮子系統(tǒng)內(nèi)部任意一塊空心板任意一種失效模式之間均存在相關(guān)性。

以上考慮相關(guān)性的失效模式均需要采用Copula模型衡量相關(guān)關(guān)系，并計算聯(lián)合失效概率，而后按第1部分公式求解體系失效概率。以ZM，ZV分別表示抗彎、抗剪失效，數(shù)字表示空心板編號，則由Copula模型得到子系統(tǒng)1內(nèi)的聯(lián)合失效概率及相關(guān)關(guān)系，如表7所示。

表7 二元Copula聯(lián)合分布結(jié)果

比較各個不同失效模式間的聯(lián)合失效概率值，作出柱狀圖如圖6。

圖6 子系統(tǒng)1不同失效模式間的二元聯(lián)合失效概率值

由圖6可見，子系統(tǒng)內(nèi)不同失效模式的聯(lián)合失效概率存在巨大差異。其中不同空心板抗彎失效的聯(lián)合失效概率(ZM_1&ZM_2)最大，其次是不同空心板抗剪失效的聯(lián)合失效概率(ZV_1&ZV_2)，二者相差約21倍，二者邊緣分布間的Kendall相關(guān)系數(shù)均在0.95左右，而邊緣失效概率值則分別為(3.57×10-5,3.74×10-5)、(8.30×10-6,1.60×10-6)，可見在邊緣分布相關(guān)程度相當(dāng)?shù)那闆r下，邊緣分布失效概率值越大則聯(lián)合失效概率值越大。

另一方面，將ZV_1&ZM_2的聯(lián)合失效概率與ZV_1&ZV_2比較，前者邊緣分布概率為(8.30×10-6, 3.74×10-5)，后者為(8.30×10-6,1.60×10-6)，發(fā)現(xiàn)前者的邊緣失效概率遠大于后者，而后者聯(lián)合失效概率卻遠大于前者，比較二者的Kendall相關(guān)系數(shù)，前者為0.22，后者為0.94，可見邊緣分布的相關(guān)性越強，則聯(lián)合分布概率值越大。

將不同計算模式下各個子系統(tǒng)及邊跨體系失效概率列入表8。

表8 不同計算模式下的體系失效概率

對比計算模式①與計算模式②中各個子系統(tǒng)的失效概率值，發(fā)現(xiàn)二者差值非常微弱，累計差值約0.03%。原因在于同一塊空心板的抗彎失效與抗剪失效的聯(lián)合失效概率值極小，而該值小又是由抗彎失效與抗剪失效的相關(guān)性很弱導(dǎo)致，故由此可知，體系中的弱相關(guān)關(guān)系可以忽略，當(dāng)作相互獨立進行處理對失效概率的計算影響不大。

對比計算模式②③可見，計算模式③得到的子系統(tǒng)失效概率約為計算模式②得到的子系統(tǒng)失效概率的2×104，失效概率大幅增長意味著可靠性大幅下降，結(jié)構(gòu)存在安全隱患的可能性大幅提高。這一顯著差異的來源在于，計算模式②將并聯(lián)的空心板作為相互獨立關(guān)系處理，將單塊板的失效概率直接相乘，而計算模式③則是認為并聯(lián)空心板存在相關(guān)性，通過聯(lián)合概率分布求取兩塊板同時失效的概率；通過表7的計算結(jié)果可知，并聯(lián)空心板之間確實存在強相關(guān)關(guān)系，故而聯(lián)合失效概率值遠大于二者乘積。由此可見，強相關(guān)關(guān)系對體系失效概率值的影響十分顯著，在計算中應(yīng)予以關(guān)注，不可忽略。

另外對比計算模式①②③對體系失效概率的計算公式可知，對于串聯(lián)結(jié)構(gòu)而言，不考慮相關(guān)性，將構(gòu)件及其失效模式當(dāng)作相互獨立處理，得到的體系失效概率偏大，可靠性偏小，計算結(jié)果偏保守；對于并聯(lián)結(jié)構(gòu)而言，不考慮相關(guān)性則使得體系失效概率偏小，可靠度偏大，即計算結(jié)果較危險，若是并聯(lián)結(jié)構(gòu)的強相關(guān)關(guān)系被忽略，則很可能無法對結(jié)構(gòu)存在安全隱患做出正確判斷。

4 結(jié) 論

(1)本文利用Copula函數(shù)對體系可靠度進行計算，在衡量相關(guān)性時可以考慮到非線性相關(guān)，同時可以得到與數(shù)值積分精確度相近的聯(lián)合失效概率，沒有參數(shù)符合正態(tài)分布的限制，更貼合工程實際也更具適用性；

(2)本文利用Copula函數(shù)構(gòu)建多元聯(lián)合分布，結(jié)合串并聯(lián)模型可將體系層面的失效概率降維到構(gòu)件失效模式，避開體系失效模式難以顯示表達的問題；

(3)本文將Copula理論運用于空心板連續(xù)梁橋體系可靠度計算發(fā)現(xiàn)，忽略弱相關(guān)關(guān)系對體系失效概率的影響不顯著，但忽略強相關(guān)關(guān)系體系失效概率可能偏差若干數(shù)量級，使計算結(jié)果不可信賴，故在對體系進行失效概率計算時充分考慮體系內(nèi)相關(guān)性是非常有必要的。