郭亞洲， 張淑華

(河海大學(xué) 港口海岸與近海工程學(xué)院，南京 210098)

系統(tǒng)可靠度計(jì)算一直存在瓶頸，就是如何找出復(fù)雜系統(tǒng)的失效模式，并在考慮元器件間、各失效模式間的相關(guān)性下計(jì)算系統(tǒng)的可靠度。在可靠度意義上，系統(tǒng)可以看作是多個(gè)可能的失效模式組成的串聯(lián)系統(tǒng)，各失效模式間所包含的隨機(jī)變量至少是部分重合的，因此各失效模式間具有相關(guān)性。失效模式間的相關(guān)性介于完全相關(guān)和相互獨(dú)立兩種極端假定之間，它們之間相關(guān)性的度量是可靠度分析時(shí)的一大難點(diǎn)，以往一般采取的獨(dú)立假定雖然大大簡化了計(jì)算，但結(jié)果偏于危險(xiǎn)。目前考慮變量間、失效模式間相關(guān)性的方法主要有二階界限法[1]，相關(guān)系數(shù)大于0.6時(shí)，界限范圍較大，精度較低；Nataf變換法[2]將非正態(tài)變量轉(zhuǎn)換為正態(tài)變量，通過求解多維正態(tài)積分計(jì)算可靠度，方法簡便但計(jì)算困難。Rosenblatt變換法[3]，將相關(guān)變量轉(zhuǎn)換為獨(dú)立變量，但需獲知精確的聯(lián)合概率分布函數(shù)，很難應(yīng)用于實(shí)際工程中；簡化多維積分計(jì)算的降維分解法[4-6]：結(jié)果精確度足夠，適用于極限狀態(tài)方程易于分解的；近似計(jì)算法：PNET法(概率網(wǎng)絡(luò)估算法)[7]，可以識別出代表失效模式，并用其計(jì)算系統(tǒng)可靠度，減少計(jì)算量且精度較高。Monte-Carlo方法[8]及其改進(jìn)是精確方法，能計(jì)算任何系統(tǒng)的可靠度，對于大型復(fù)雜系統(tǒng)工作量過大，常作為檢驗(yàn)新方法的標(biāo)準(zhǔn)解。這些理論大多建立在線性相關(guān)系數(shù)ρ基礎(chǔ)上，存在以下缺陷：只能描述線性相關(guān)關(guān)系，無法描述非線性的相關(guān)關(guān)系，例如x和x2之間的關(guān)系；ρ的計(jì)算受聯(lián)合概率密度函數(shù)及相關(guān)元個(gè)數(shù)影響較大，對于概率密度函數(shù)較為復(fù)雜或存在多個(gè)相關(guān)元的結(jié)構(gòu)，ρ的求解較為困難；ρ通常用于描述正態(tài)變量之間的相關(guān)性，非正態(tài)變量需要先轉(zhuǎn)換為正態(tài)變量，局限性較大。因此本文提出了一個(gè)基于copula的考慮失效模式間相關(guān)性的系統(tǒng)可靠度分析方法，可以更精確地解決工程可靠度問題。copula方法是近年發(fā)展起來的方法，目前主要應(yīng)用于水文、金融和結(jié)構(gòu)工程等領(lǐng)域[9-12]，由于其在處理變量之間相關(guān)性方面的獨(dú)特優(yōu)勢，逐步被引入到可靠性領(lǐng)域，唐小松等[13]基于二維 Copula 函數(shù)構(gòu)建了變量的聯(lián)合分布函數(shù)，采用直接積分對結(jié)構(gòu)進(jìn)行了可靠度分析；Tang 等[14]研究了不同二維 Copula函數(shù)構(gòu)建的二維變量相關(guān)模型的差異，分析了其對可靠性分析結(jié)果的影響；Jiang等[15]基于Vine Copula函數(shù)構(gòu)建了多維可靠性分析模型，為多維相關(guān)性的結(jié)構(gòu)可靠性問題提出了一個(gè)新的解決途徑。以上是copula在可靠性方面應(yīng)用的初步探索，有待進(jìn)一步發(fā)展。Copula函數(shù)在處理相關(guān)性方面有如下優(yōu)勢：copula函數(shù)可以將隨機(jī)變量聯(lián)合分布的構(gòu)建拆分為邊緣分布形式的估計(jì)和邊緣分布之間的相關(guān)結(jié)構(gòu)兩個(gè)部分，使聯(lián)合分布的構(gòu)建更為方便；copula函數(shù)能夠描述邊緣分布為非正態(tài)分布的變量間的相關(guān)性；copula能描述變量之間的非線性相關(guān)甚至尾部相關(guān)。本文的方法為非正態(tài)變量之間的線性相關(guān)關(guān)系的描述提供了一種新的途徑。

1 copula計(jì)算模型

1.1 copula理論

Copula函數(shù)是聯(lián)合分布與邊緣分布之間的連接函數(shù)，本質(zhì)上它也是一種聯(lián)合分布函數(shù)。已知n維隨機(jī)變量X=(X1，X2，…，Xn);由Sklar[16]定理,存在函數(shù)C使得任意分布的邊緣分布如：F1(X1)，…，F(xiàn)n(Xn)的聯(lián)合分布函數(shù)F表示為

F(X1，X2，…，Xn)=C(F1(X1)，F(xiàn)2(X2)，…，F(xiàn)n(Xn))

(1)

那么隨機(jī)變量X的聯(lián)合概率密度函數(shù)為

(2)

式中：c為copula函數(shù)C的密度函數(shù)，所以c的表達(dá)式為

(3)

Copula函數(shù)通常使用Kendall秩相關(guān)系數(shù)τ或Spearman秩相關(guān)系數(shù)ρ作為相關(guān)性測度，它們與copula的相關(guān)系數(shù)θ(用以表示相關(guān)性)關(guān)系如下

(4)

1.2 常規(guī)的copula函數(shù)

現(xiàn)有的copula應(yīng)用大多針對二維問題，表1中列舉了幾種常見的copula函數(shù)。

表1 常見的copula函數(shù)Tab.1 Common copula function

1.3 系統(tǒng)可靠度的copula計(jì)算模型

本文copula計(jì)算模型的構(gòu)建主要分為兩步：一是copula的選擇，二是模型的構(gòu)建。

1.3.1 copula選擇

本文采用基于非參數(shù)核密度估計(jì)的copula選擇原理[17]，其核心思想是采用非參數(shù)核密度估計(jì)樣本邊緣分布，從待選copula函數(shù)中初步篩選出能描述失效模式間相關(guān)性的copula函數(shù)，同時(shí)采用copulafit函數(shù)估計(jì)copula中的相關(guān)參數(shù)，最終通過比較幾種copula與經(jīng)驗(yàn)copula[18]的歐氏平方距離選出最優(yōu)copula。

1.3.2 計(jì)算模型構(gòu)建

(5)

式中：ui=FXi(i=1,2,…,n)；u=(u1,u2,…,un)。

根據(jù)式(6)在邊緣分布已知情況下，可以選擇合適的copula函數(shù)構(gòu)造聯(lián)合分布函數(shù)。設(shè)系統(tǒng)有k個(gè)不同的失效模式，其相應(yīng)的功能函數(shù)為

Yj=zj(X1，X2，…，Xn)j=1,2,…,k

(6)

則Ej=[zj(X1，X2，…，Xn)≤0]表示第j個(gè)失效模式出現(xiàn)，系統(tǒng)的失效概率表示為

(7)

結(jié)構(gòu)系統(tǒng)主要分為兩類：串聯(lián)系統(tǒng)(只要任意一個(gè)極限狀態(tài)失效，則認(rèn)為系統(tǒng)失效)和并聯(lián)系統(tǒng)(所有極限狀態(tài)全部失效，則認(rèn)為系統(tǒng)失效)，結(jié)合copula理論推導(dǎo)出相應(yīng)的系統(tǒng)可靠度計(jì)算模型。

(1)串聯(lián)系統(tǒng)。

系統(tǒng)失效概率為

(8)

(2)并聯(lián)系統(tǒng)。

系統(tǒng)失效概率為

(9)

注：基準(zhǔn)面采用黃海平均海平面。圖1 秦皇島某直立堤斷面示意Fig.1 Cross section of the vertical breakwater at Qinhuangdao port

2 防波堤可靠度分析

2.1 防波堤概況

本文的實(shí)例計(jì)算采用的是秦皇島某防波堤[19]。此防波堤形式為沉箱式直立堤，沉箱以及上部結(jié)構(gòu)的材料全部采用鋼筋混凝土，沉箱內(nèi)填充塊石。防波堤總長為250 m,每個(gè)沉箱縱向長度為12.5 m,共分為6個(gè)艙格，防波堤各部分高程詳見圖1。

2.2 失效概率計(jì)算

表2 荷載統(tǒng)計(jì)參數(shù)Tab.2 Statistic parameters of the loads

表2中P和Pu分別代表水平波浪力和浮托力，Mp、Mpu為其相應(yīng)的力矩，此外重力G和摩擦系數(shù)f為定值，參考文獻(xiàn)[22]中的數(shù)值，G取945.46 kN·m-1，f取0.6，重力的力矩MG計(jì)算得5 923.59 kN·m·m-1。

本文選取的是較常見的滑動(dòng)(Z1)、傾覆破壞狀態(tài)(Z2)，另一種滑移破壞狀態(tài)(Z3)失效概率過小，對系統(tǒng)可靠度影響較小，可以忽略[21]。所以文中列出的兩種代表性的失效模式完全可以代表全部失效模式，并用其求解系統(tǒng)失效概率，其極限狀態(tài)方程如下

Z1=g(G,P,PU,f)=(G-PU)·f-P=0

(10)

Z2=MG-MP-MPU=0

(11)

Z3=(G+G1)·f2-P(地基為砂)Z3=Su·l1-P(地基為飽和粘性土)

(12)

式中：重力G1為擦破裂面所圍的拋石的重力；系數(shù)f2為拋石基床與天然地基的摩擦系數(shù)；Su為地基土的粘聚力；l1為基床滑裂段的長度。

本文采用蒙特卡羅數(shù)值模擬(Monte Carlo Simulations, MCS)方法[22]計(jì)算每個(gè)失效模式的失效概率。MCS 方法計(jì)算結(jié)構(gòu)失效概率時(shí)，主要分兩步：一根據(jù)各隨機(jī)變量分布形式產(chǎn)生隨機(jī)數(shù)，二將其代入極限狀態(tài)方程得到相應(yīng)失效模式的失效概率。G和f看作定值，水平波浪力和波浪浮托力均采用Gumbel 分布，其分布形式為

F(x)=exp{-exp[-α(x-β)}

(13)

根據(jù)隨機(jī)變量的分布形式(13)采用MATLAB編程產(chǎn)生隨機(jī)數(shù)，帶入極限狀態(tài)方程，得出滑移失效模式以及傾覆失效模式的失效概率，并與以往文獻(xiàn)[23]的計(jì)算結(jié)果進(jìn)行對比，驗(yàn)證本文所得結(jié)果的正確性，詳見表3。

表3 失效概率計(jì)算結(jié)果Tab.3 Calculation results of failure probability

由表中數(shù)據(jù)看，本文失效概率計(jì)算結(jié)果誤差在允許范圍內(nèi)，滿足工程精度要求，可以應(yīng)用于求解系統(tǒng)失效概率。

2.3 失效相關(guān)的系統(tǒng)可靠度計(jì)算

2.3.1 確定邊緣分布

令X、Y分別代表滑移失效和傾覆失效的功能函數(shù)，X、Y的樣本是由上面蒙特卡洛模擬產(chǎn)生的隨機(jī)數(shù)帶入功能函數(shù)所得。用樣條差值法求得X、Y的經(jīng)驗(yàn)分布，將其與采用核密度估計(jì)法得到的核分布估計(jì)形式進(jìn)行對比，見圖2、圖3。

圖2 滑移失效功能函數(shù)分布Fig.2 Distribution of slip failure function圖3 傾覆失效功能函數(shù)分布Fig.3 Distribution of overturn failure function

圖4 二元頻率直方圖Fig.4 Bivariate frequency histogram

由圖2和圖3可知：樣本的經(jīng)驗(yàn)分布和核分布估計(jì)圖像幾乎重合，差別較小，所以經(jīng)驗(yàn)分布和核分布估計(jì)形式都可以作為樣本的邊緣分布，令U=F(X)和V=F(Y)分別代表X和Y的分布形式。

2.3.2 選取適當(dāng)?shù)腸opula函數(shù)

根據(jù)(Ui,Vi)(i=1,2,…n)繪制二元頻率直方圖，觀察其形狀和特征，初步從表1中列舉的copula函數(shù)中選擇符合的copula函數(shù)，詳見圖4。

由圖中可得，尾部基本上對稱，由表1中所列各copula函數(shù)的尾部特性可知，二元正態(tài)copula和二元Frank-copula適合于尾部對稱且漸進(jìn)獨(dú)立；二元t-copula 適合于尾部對稱且相關(guān)，所以根據(jù)圖像初步選擇二元正態(tài)copula或二元t-copula 來擬合之間的相關(guān)結(jié)構(gòu)。

2.3.3 參數(shù)估計(jì)

基于X和Y的樣本數(shù)據(jù)，本文采用MATLAB自帶函數(shù)copulafit分別確定二元正態(tài)copula和二元t-copula 的相關(guān)參數(shù)。

二元正態(tài)copula中的相關(guān)參數(shù)ρ的估計(jì)值為

二元t-copula中的相關(guān)參數(shù)ρ和自由度k的估計(jì)值為

二元正態(tài)copula的表達(dá)式為

(14)

二元t-copula的表達(dá)式為

(15)

將二元正態(tài)copula和二元t-copula的相關(guān)參數(shù)分別帶入式(14)、(15)，計(jì)算得出二元正態(tài)copula和t-copula 的密度函數(shù)和分布函數(shù)值，繪制密度函數(shù)和分布函數(shù)圖，見圖5～圖8。

圖5 二元正態(tài)copula密度函數(shù)圖Fig.5 Density function graph of bivariate normal-copula圖6 二元正態(tài)copula分布函數(shù)圖Fig.6 Distribution function graph of bivariate normal-copula

圖7 二元t-copula密度函數(shù)圖Fig.7 Density function graph of bivariate t-copula圖8 二元t-copula分布函數(shù)圖Fig.8 Distribution function graph of bivariate t-copula

由圖可得：二元正態(tài)copula和二元t-copula對X與Y之間相關(guān)結(jié)構(gòu)的擬合程度都比較理想，二元t-copula尖尾特性明顯，但由于二者之間差距較小，所以需要進(jìn)一步通過與經(jīng)驗(yàn)copula之間的歐氏平方距離篩選出最優(yōu)copula。

2.3.4 模型評價(jià)

計(jì)算二元正態(tài)copula、t-copula 與經(jīng)驗(yàn)copula之間的歐氏平方距離，距離最小的為最優(yōu)copula。記X和Y的經(jīng)驗(yàn)分布函數(shù)分別為Fn(x)和Gn(y)，所以定義樣本的經(jīng)驗(yàn)copula為

(16)

式中：I[ ]為示性函數(shù)(當(dāng)Fn(xi)≤u時(shí)，I[Fn(xi)≤u]=1；否則I[Fn(xi)≤u]=0)。

則二元正態(tài)copula和二元t-copula和經(jīng)驗(yàn)copula之間的歐氏平方距離分別為

(17)

由表4可知，二元正態(tài)copula的歐氏平方距離更小，對于X,Y之間相關(guān)結(jié)構(gòu)的擬合更為理想，所以最終選取二元正態(tài)copula擬合X,Y之間的相關(guān)結(jié)構(gòu)。

2.3.5 系統(tǒng)失效概率

將表3中求解出的X和Y的失效概率和正態(tài)copula中的相關(guān)參數(shù)帶入式(9)和(15)中得系統(tǒng)失效概率Pf，并與采用蒙特卡洛模擬[23]求得的系統(tǒng)失效概率進(jìn)行對比，結(jié)果見表5。

表4 歐氏平方距離計(jì)算結(jié)果Tab.4 Calculation results of Squared Euclidean distance

表5 系統(tǒng)可靠度計(jì)算結(jié)果對比Tab.5 Comparison of system reliability calculation results

由表5可知：本文構(gòu)建的基于copula的系統(tǒng)可靠度計(jì)算模型所得結(jié)果在誤差允許范圍內(nèi)，滿足工程精度要求。

3 結(jié)論

本文基于copula函數(shù)，提出了一種結(jié)構(gòu)系統(tǒng)內(nèi)失效模式相關(guān)的結(jié)構(gòu)體系可靠度的計(jì)算方法，克服了傳統(tǒng)求解結(jié)構(gòu)體系可靠度將失效模式相互獨(dú)立假設(shè)的局限性，為結(jié)構(gòu)體系可靠度的求解提供了一種新思路。以秦皇島某防波堤為例計(jì)算了結(jié)構(gòu)可靠度并與Monte Carlo方法計(jì)算結(jié)果進(jìn)行對了對比，驗(yàn)證了本文理論的實(shí)用性與結(jié)果的正確性。