相關(guān)隨機(jī)變量拉丁超立方抽樣及隧道掌子面 穩(wěn)定性概率分析

章從旭,藍(lán)元盛,李 斌

(1. 中交路橋建設(shè)有限公司,北京 100027; 2. 武漢理工大學(xué) 交通與物流工程學(xué)院,湖北 武漢 430063)

0 引 言

隨著可靠度理論的不斷發(fā)展,可靠度設(shè)計(jì)在巖土工程的應(yīng)用越來(lái)越廣。可靠度設(shè)計(jì)方法包括一次二階矩法、一階可靠度法、響應(yīng)面法、數(shù)據(jù)表法和蒙特卡洛法等。其中,蒙特卡洛法是最為直接和準(zhǔn)確的方法[1-2]。蒙特卡洛法采集的樣本可通過(guò)有限元計(jì)算得到失效概率,分析結(jié)果可作為其它方法計(jì)算結(jié)果的唯一驗(yàn)證[3-9]。然而,蒙特卡洛法的計(jì)算成本較高,對(duì)于失效概率較小的問(wèn)題,有時(shí)幾乎不可能實(shí)現(xiàn)。

M. D. MCKAY等[10]在1979年指出拉丁超立方抽樣技術(shù)可降低計(jì)算消耗,該技術(shù)具有分層均勻、抽樣少的優(yōu)點(diǎn)。具體計(jì)算步驟:① 確定抽樣樣本數(shù)N;② 將變量的累積分布分成相同的N個(gè)小區(qū)間,即將變量累積分布(0,1)區(qū)間均分成N段,使得每個(gè)區(qū)間有相同的概率;③ 從N個(gè)區(qū)間段中分別隨機(jī)抽取一個(gè)值,將抽取的N個(gè)值分別通過(guò)變量的分布反函數(shù)映射為變量的樣本,共計(jì)得到N個(gè)樣本。隨機(jī)變量的拉丁超立方抽樣原理如圖1。圖1(a)表示將樣本的概率分布范圍劃分為5個(gè)概率相等的區(qū)間,每個(gè)區(qū)間采集一個(gè)樣本點(diǎn);圖1(b)表示2個(gè)獨(dú)立樣本經(jīng)分層抽樣后的組合情況,2個(gè)變量在任何一個(gè)分割的概率區(qū)間只出現(xiàn)一次,以保證所采集的樣本最具代表性。拉丁超立方抽樣保證了尾部層中的樣本也會(huì)被抽取,而這些樣本對(duì)于失效概率計(jì)算來(lái)說(shuō)是非常重要的;拉丁超立方抽取樣本數(shù)量越多,分層就越多,極端的尾部值就更容易被抽取。

圖1 拉丁超立方抽樣原理

拉丁超立方抽樣假定所有隨機(jī)變量是相互獨(dú)立的,而在巖土工程中很多參數(shù)往往具有相關(guān)性,因此,無(wú)法直接進(jìn)行拉丁超立方抽樣,需要進(jìn)行相應(yīng)的變換。具體變換方式可參考相關(guān)隨機(jī)變量蒙特卡洛抽樣。常見(jiàn)4種不同類型相關(guān)隨機(jī)變量蒙特卡洛抽樣的實(shí)現(xiàn)步驟如下:

1)對(duì)于均為正態(tài)分布的多維相關(guān)隨機(jī)變量,可基于協(xié)方差矩陣的Choleshy因子分解將其轉(zhuǎn)換為獨(dú)立標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量進(jìn)行蒙特卡洛抽樣[11]。

2)對(duì)于均為對(duì)數(shù)正態(tài)分布的多維相關(guān)隨機(jī)變量,可將其轉(zhuǎn)換為正態(tài)分布的多維相關(guān)隨機(jī)變量進(jìn)行蒙特卡洛抽樣,轉(zhuǎn)換后協(xié)方差矩陣跟隨改變[12]。

3)對(duì)于非正態(tài)分布的多維相關(guān)隨機(jī)變量,可采用Nataf變換分2步將其轉(zhuǎn)換為獨(dú)立標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量進(jìn)行蒙特卡洛抽樣[13]。

4)對(duì)于已知秩相關(guān)矩陣的多維相關(guān)隨機(jī)變量,可將其轉(zhuǎn)換為具有相同秩相關(guān)矩陣及相同維度的標(biāo)準(zhǔn)正態(tài)分布的相關(guān)隨機(jī)變量,再根據(jù)相關(guān)隨機(jī)變量的邊緣分布生成所需要的樣本,使得每個(gè)維度上相關(guān)隨機(jī)變量樣本從大到小的排序與對(duì)應(yīng)維度上標(biāo)準(zhǔn)正態(tài)分布相關(guān)隨機(jī)變量蒙特卡洛抽樣樣本從大到小的排序一致,從而實(shí)現(xiàn)秩相關(guān)隨機(jī)變量的蒙特卡洛抽樣[14]。

筆者參考相關(guān)隨機(jī)變量蒙特卡洛抽樣的實(shí)現(xiàn),結(jié)合拉丁超立方抽樣原理,采用Python語(yǔ)言編程實(shí)現(xiàn)了上述4種類型的相關(guān)隨機(jī)變量拉丁超立方抽樣。

1 相關(guān)隨機(jī)變量拉丁超立方抽樣原理及實(shí)現(xiàn)

1.1 相關(guān)隨機(jī)變量Pearson相關(guān)

1.1.1 相關(guān)隨機(jī)變量的正態(tài)分布

正態(tài)分布進(jìn)行線性變化后仍然是正態(tài)分布,且相關(guān)性不變,因此,任意滿足正態(tài)分布的相關(guān)隨機(jī)變量均可以通過(guò)標(biāo)準(zhǔn)正態(tài)分布的線性組合得到,表達(dá)式如式(1):

x=A·ξ+μ

(1)

式中:x為一組正態(tài)分布的相關(guān)隨機(jī)變量;μ為各相關(guān)隨機(jī)變量的均值矩陣;ξ為一組相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布變量;A為轉(zhuǎn)換矩陣。

G. ZEROVNIK等[12]證明了V=AAT,其中V為相關(guān)隨機(jī)變量x的協(xié)方差系數(shù)矩陣。又知V是一個(gè)對(duì)稱的正定矩陣,可以通過(guò)Choleshy分解成2個(gè)矩陣的乘積,因此,只要得到V的Choleshy分解矩陣就可以求得矩陣A。

在Python語(yǔ)言中,使用numpy庫(kù)實(shí)現(xiàn)Choleshy分解,得到A=np.linalg.cholesky(V)。如果是給出的相關(guān)系數(shù)矩陣R,則V=σRσT,σ為對(duì)角陣,對(duì)角線上的元素為各相關(guān)隨機(jī)變量的標(biāo)準(zhǔn)差。所以,對(duì)x的拉丁超立方抽樣可以轉(zhuǎn)化為對(duì)獨(dú)立標(biāo)準(zhǔn)正態(tài)分布ξ的拉丁超立方抽樣。實(shí)現(xiàn)拉丁超立方抽樣用到Python的scipy.stats中的qmc和norm。

1.1.2 相關(guān)隨機(jī)變量的對(duì)數(shù)正態(tài)分布

對(duì)數(shù)正態(tài)分布拉丁超立方抽樣與正態(tài)分布拉丁超立方抽樣具有相似性,若xi服從對(duì)數(shù)正態(tài)分布,則lnxi服從正態(tài)分布,即

(2)

式中:xi為第i個(gè)相關(guān)隨機(jī)變量;Aij為轉(zhuǎn)換系數(shù);ξj為第j個(gè)獨(dú)立標(biāo)準(zhǔn)正態(tài)分布變量;μi為lnxi的均值。

令yi=lnxi

(3)

則,式(2)可以寫(xiě)成:

y=A·ξ+μ

(4)

可見(jiàn),式(4)的形式同式(1)。

因此,一組滿足對(duì)數(shù)正態(tài)分布的相關(guān)隨機(jī)變量的拉丁超立方抽樣,可以轉(zhuǎn)換為一組滿足正態(tài)分布相關(guān)隨機(jī)變量的拉丁超立方抽樣。

根據(jù)G. ZEROVNIK等[12]研究,原對(duì)數(shù)正態(tài)分布x與轉(zhuǎn)換后正態(tài)分布y之間的協(xié)方差系數(shù)的關(guān)系如式(5):

(5)

根據(jù)式(6)可以求得轉(zhuǎn)換后yi的均值:

(6)

式中:E(yi)為第i個(gè)y變量的均值;E(xi)為第i個(gè)x變量的均值;D(xi)為第i個(gè)x變量的方差。

得到了正態(tài)分布y的A與μ,后續(xù)進(jìn)行拉丁超立方抽樣的計(jì)算步驟與正態(tài)分布的拉丁超立方抽樣一致。得到相關(guān)隨機(jī)變量yi的樣本后,通過(guò)xi=eyi轉(zhuǎn)換,即得到相關(guān)隨機(jī)變量x的樣本。

1.1.3 相關(guān)隨機(jī)變量的非正態(tài)分布

LIU Peiling等[15-16]提出的Nataf轉(zhuǎn)換法,將非正態(tài)分布轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)分布空間下,只需要知道變量的邊緣分布以及變量間的相關(guān)系數(shù),即可以把原始分布空間轉(zhuǎn)換到獨(dú)立標(biāo)準(zhǔn)正態(tài)空間。Nataf轉(zhuǎn)換分2步:①將原始空間轉(zhuǎn)換成相關(guān)的標(biāo)準(zhǔn)正態(tài)空間;②將相關(guān)的標(biāo)準(zhǔn)正態(tài)空間轉(zhuǎn)換成獨(dú)立的標(biāo)準(zhǔn)正態(tài)的空間,這可通過(guò)正態(tài)分布相關(guān)隨機(jī)變量的拉丁超立方抽樣來(lái)解決。

原始空間向相關(guān)的標(biāo)準(zhǔn)正態(tài)空間轉(zhuǎn)換的操作為:輸入一組相關(guān)隨機(jī)變量X=(x1,x2,…,xn),若相關(guān)隨機(jī)變量xi(i=1,2,…,n)的概率密度函數(shù)fi(xi)和累積密度函數(shù)Fi(xi)已知,通過(guò)等概率轉(zhuǎn)換原則[15-16]可得:

(7)

式中:yi為一組標(biāo)準(zhǔn)正態(tài)相關(guān)隨機(jī)變量Y=(y1,y2,…,yn)中的變量;Φ(·)、Φ-1(·)分別為標(biāo)準(zhǔn)正態(tài)累積密度分布函數(shù)和逆累積密度分布函數(shù)。

設(shè)X中任意2個(gè)相關(guān)隨機(jī)變量xi、xj的相關(guān)系數(shù)為ρxixj,則

(8)

式中:μxi、μxj分別為xi、xj的均值;Δxi、Δxj分別為xi、xj的標(biāo)準(zhǔn)差;ρyiyj為yi、yj的相關(guān)系數(shù)。

在二元標(biāo)準(zhǔn)正態(tài)分布的聯(lián)合概率密度函數(shù)數(shù)值計(jì)算中,式(8)的雙重積分上下限選擇[-4,4]就足夠精確了[13]。在ρxixj、μxi、μxj已知的情況下,幾乎不可能直接計(jì)算ρyiyj,因此,需采用數(shù)值計(jì)算法求解ρyiyj。當(dāng)ρxixj>0時(shí),具體步驟如下:

Step 1在[0,1]區(qū)間假設(shè)一個(gè)ρyiyj值。

Step 2將假設(shè)的ρyiyj值代入式(8),計(jì)算得到ρxixj。

當(dāng)ρxixj<0時(shí),在[-1,0]區(qū)間進(jìn)行類似操作。

1.1.4 相關(guān)隨機(jī)變量拉丁超立方抽樣舉例

以巖土工程中土體常見(jiàn)的黏聚力c、內(nèi)摩擦角φ為例進(jìn)行拉丁超立方抽樣,抽取10 000個(gè)樣本。均值c=10 kPa,φ=30°,標(biāo)準(zhǔn)差Δc=2 kPa,Δφ=2°,Pearson相關(guān)系數(shù)為-0.5。

分別進(jìn)行正態(tài)分布相關(guān)隨機(jī)變量抽樣(c、φ均滿足正態(tài)分布)、對(duì)數(shù)正態(tài)分布相關(guān)隨機(jī)變量抽樣(c、φ均滿足對(duì)數(shù)正態(tài)分布)、非正態(tài)分布相關(guān)隨機(jī)變量抽樣(c滿足對(duì)數(shù)正態(tài)分布,φ滿足正態(tài)分布),抽樣結(jié)果分布如圖2,抽樣參數(shù)c、φ概率密度如圖3。

(a) 正態(tài)分布抽樣

(b) 對(duì)數(shù)正態(tài)分布抽樣

(c) 非正態(tài)分布抽樣

圖3 正態(tài)分布、對(duì)數(shù)正態(tài)分布、非正態(tài)分布抽樣概率密度

由圖3可見(jiàn),當(dāng)Pearson相關(guān)系數(shù)不等于0時(shí),c與φ的抽樣概率密度分布曲線與理論概率密度曲線存在些許差異。分析原因,一方面可能是拉丁超立方抽樣樣本數(shù)量相對(duì)較少造成的;另一方面,可能是每個(gè)區(qū)間抽樣的隨機(jī)性造成的,尤其是在均值附近。

綜上,隨機(jī)變量Pearson相關(guān)的3種不同類型的拉丁超立方抽樣保持了原變量的分布,在保持抽樣精度的同時(shí),大大降低了樣本數(shù)量。

1.2 隨機(jī)變量Spearman秩相關(guān)

Spearman秩相關(guān)系數(shù)只對(duì)變量的秩次大小進(jìn)行線性相關(guān)分析,所以沒(méi)有要求變量的分布,但秩相關(guān)系數(shù)與樣本的排序有關(guān),而與樣本具體的值關(guān)聯(lián)不大,因此,秩相關(guān)隨機(jī)變量的拉丁超立方抽樣可通過(guò)秩相關(guān)矩陣得到。

相關(guān)隨機(jī)變量Z如果滿足Z=ξPT,而VS=PPT,VS為變量Z秩相關(guān)矩陣,則Z與ξPT具有相同的秩相關(guān)系數(shù),其中ξ為一組相互獨(dú)立的隨機(jī)變量[14]。

既然秩相關(guān)系數(shù)與相關(guān)隨機(jī)變量具體的分布沒(méi)有關(guān)系,那么可以設(shè)ξ為標(biāo)準(zhǔn)正態(tài)分布,即Z為正態(tài)分布。

假設(shè)一組包含m個(gè)秩相關(guān)的相關(guān)隨機(jī)變量X=(x1,x2,…,xm)需抽取N個(gè)樣本,秩相關(guān)矩陣為VS,通過(guò)Python語(yǔ)言實(shí)現(xiàn)拉丁超立方抽樣過(guò)程如下:

Step 1通過(guò)拉丁超立方抽樣對(duì)隨機(jī)變量ξ進(jìn)行抽樣,ξ=(ξ1,ξ2,…,ξm),其中:ξ1、ξ2、…、ξm為相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量,每個(gè)變量共抽取N個(gè)樣本。

Step 2對(duì)VS進(jìn)行Choleshy分解得到PT,通過(guò)Z=ξPT轉(zhuǎn)換得到正態(tài)分布Z的N個(gè)樣本,分別對(duì)Z=(z1,z2,…,zm)中每個(gè)變量的N個(gè)樣本進(jìn)行排序,得到排序序號(hào)。

Step 3通過(guò)拉丁超立方抽樣分別對(duì)xi(i=1,2,…,m)抽取N個(gè)樣本,并對(duì)這N個(gè)樣本進(jìn)行排序,排序序號(hào)與第i個(gè)變量zi的N個(gè)樣本排序號(hào)相同,從而得到秩相關(guān)隨機(jī)變量X=(x1,x2,…,xm)的N個(gè)拉丁超立方抽樣樣本。

以巖土工程中土體常見(jiàn)的參數(shù)c、φ為例進(jìn)行拉丁超立方抽樣,抽取10 000個(gè)樣本。兩個(gè)參數(shù)的均值u分別為uc=10 kPa、uφ=30°,標(biāo)準(zhǔn)差Δc=2 kPa、Δφ=2°,其中c滿足對(duì)數(shù)正態(tài)分布、φ滿足正態(tài)分布,Spearman秩相關(guān)系數(shù)為-0.5。抽樣結(jié)果如圖4。

圖4 Spearman秩相關(guān)抽樣結(jié)果

2 隧道掌子面失效概率分析

2.1 實(shí)例背景

目前,隧道工程仍以確定性分析為主,但其圍巖的材料參數(shù)往往具有不確定性,因此,有必要對(duì)某些情況進(jìn)行不確定性分析或可靠度設(shè)計(jì)[17]。不確定性分析方法可以是蒙特卡洛抽樣和數(shù)值計(jì)算方法的結(jié)合[18]。為了展示拉丁超立方抽樣在工程應(yīng)用中的優(yōu)越性,筆者用一個(gè)直徑D=6 m,埋深H=6 m的圓柱形隧道掌子面進(jìn)行失效概率分析。隧道周圍土體采用摩爾-庫(kù)倫模型,襯砌采用彈性模型。將土體的黏聚力c和內(nèi)摩擦角φ考慮為服從正態(tài)分布的隨機(jī)變量,兩個(gè)參數(shù)的均值u分別為uc=10 kPa、uφ=30°,標(biāo)準(zhǔn)差Δc=2 kPa、Δφ=2°,Pearson相關(guān)系數(shù)為-0.5。土體的重度γ=18 kN/m3,彈性模量E=240 MPa,泊松比ν=0.3。隧道采用全斷面開(kāi)挖,不考慮超前支護(hù)措施。

2.2 三維數(shù)值模型

圖5為概率分析所用的三維數(shù)值模型,考慮隧道結(jié)構(gòu)和荷載的對(duì)稱性,僅建立了1/2模型進(jìn)行分析。為避免邊界效應(yīng)的影響,根據(jù)圣維南原理,隧道開(kāi)挖邊界距模型邊界應(yīng)為3～5倍洞徑D,因此分析時(shí),隧道開(kāi)挖邊界至模型的水平邊界以及模型底部邊界的距離均取3D,掌子面距模型縱向邊界的距離也為3D。在模型的縱向上,掌子面前方一倍洞徑的土體網(wǎng)格單元為0.5 m,其它部分的網(wǎng)格縱向尺寸為2 m;模型單元為57 400個(gè),節(jié)點(diǎn)為62 118個(gè);邊界條件上,模型的左右和前后邊界為法向約束,模型的底部為全約束。

2.3 模擬過(guò)程

要計(jì)算圖5數(shù)值模型的隧道掌子面的失效概率,需抽取一定數(shù)量的樣本進(jìn)行確定性分析。為保證一定的精度,樣本的數(shù)量N應(yīng)滿足

N≥100/Pf

(10)

式中:Pf為失效概率。

筆者采用批量采樣的方式,進(jìn)行確定性分析,統(tǒng)計(jì)失效樣本的數(shù)量,每批次采集1 000個(gè)樣本;重復(fù)這一步驟,直到滿足NPf≥100,終止采樣,計(jì)算最終的失效概率。

對(duì)采集的每個(gè)樣本(ci,φi),傳統(tǒng)方法是通過(guò)FLAC3D自帶的強(qiáng)度折減法來(lái)計(jì)算安全系數(shù)f的。f≥1,表示隧道掌子面穩(wěn)定;f<1,表示隧道掌子面不穩(wěn)定。

為提高計(jì)算效率,筆者不計(jì)算安全系數(shù)具體值,而是通過(guò)運(yùn)行命令“solve fos bracket 1.0 1.0”得到計(jì)算的最大不平衡力比值r。當(dāng)r<1×10-5時(shí),表明迭代計(jì)算是收斂的,即在不進(jìn)行折減的情況下,隧道掌子面是穩(wěn)定的;當(dāng)r≥1×10-5時(shí),掌子面不穩(wěn)定。

圖6為數(shù)值模型在某一黏聚力c與內(nèi)摩擦角φ組合下的計(jì)算結(jié)果?？梢?jiàn),模型基本沒(méi)有塑性區(qū),r=9.96 × 10-6,因此,掌子面是穩(wěn)定的。

圖6 三維數(shù)值模型計(jì)算結(jié)果

2.4 結(jié)果分析

圖7為用蒙特卡洛法計(jì)算得到的隧道掌子面失效概率分析結(jié)果。

圖7 樣本數(shù)量與失效概率

由圖7可見(jiàn),當(dāng)樣本數(shù)量N=97 000個(gè)時(shí),掌子面失效樣本的數(shù)量為100個(gè),此時(shí)的失效概率Pf=100/97 000≈0.001。在蒙特卡洛模擬中,95%置信區(qū)間允許的誤差為[19]:

(11)

即,滿足95%置信區(qū)間的失效概率范圍為:[0.8Pf, 1.2Pf]=[0.000 82, 0.001 24]。

根據(jù)拉丁超立方抽樣采樣結(jié)果,當(dāng)樣本數(shù)量達(dá)到16 000次時(shí),已能滿足所需的計(jì)算精度。所需樣本數(shù)量?jī)H為蒙特卡洛抽樣樣本數(shù)的16.5%。由此可見(jiàn),拉丁超立方抽樣更具代表性,可以顯著降低概率分析所用的樣本數(shù)量。

3 結(jié) 論

拉丁超立方抽樣從多元分布中獨(dú)立抽取最具代表性的樣本,在保證計(jì)算精度的同時(shí),可大大減少抽樣樣本數(shù),但拉丁超立方抽樣不適用于相關(guān)隨機(jī)變量的直接抽樣,需進(jìn)行一定的轉(zhuǎn)換。筆者通過(guò)Python編程將4種不同類型的相關(guān)隨機(jī)變量轉(zhuǎn)換為獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量進(jìn)行拉丁超立方抽樣。具體為:正態(tài)分布與對(duì)數(shù)正態(tài)分布用協(xié)方差矩陣通過(guò)Choleshy分解實(shí)現(xiàn)相關(guān)隨機(jī)變量的轉(zhuǎn)換;非正態(tài)分布通過(guò)Nataf變換實(shí)現(xiàn)相關(guān)隨機(jī)變量的轉(zhuǎn)換;秩相關(guān)通過(guò)秩相關(guān)矩陣分解實(shí)現(xiàn)相關(guān)隨機(jī)變量的轉(zhuǎn)換。同時(shí),基于獨(dú)立標(biāo)準(zhǔn)正態(tài)分布相關(guān)隨機(jī)變量拉丁超立方抽樣,開(kāi)展了隧道掌子面穩(wěn)定性概率分析,并與蒙特卡洛模擬結(jié)果進(jìn)行了對(duì)比。研究得到以下主要結(jié)論:

1)相關(guān)隨機(jī)變量拉丁超立方抽樣的抽樣結(jié)果與其理論分布結(jié)果基本吻合,驗(yàn)證了拉丁超立方抽樣方法的準(zhǔn)確性。

2)隧道掌子面穩(wěn)定性概率分析表明,拉丁超立方抽樣所需樣本數(shù)為蒙特卡洛抽樣樣本數(shù)的16.5%即可滿足計(jì)算精度要求,極大地提高了計(jì)算效率。