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

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

0 引 言

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

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

圖1 拉丁超立方抽樣原理

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

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

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

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

4)對于已知秩相關(guān)矩陣的多維相關(guān)隨機(jī)變量,可將其轉(zhuǎn)換為具有相同秩相關(guān)矩陣及相同維度的標(biāo)準(zhǔn)正態(tài)分布的相關(guān)隨機(jī)變量,再根據(jù)相關(guān)隨機(jī)變量的邊緣分布生成所需要的樣本,使得每個維度上相關(guān)隨機(jī)變量樣本從大到小的排序與對應(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語言編程實(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ī)變量均可以通過標(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是一個對稱的正定矩陣,可以通過Choleshy分解成2個矩陣的乘積,因此,只要得到V的Choleshy分解矩陣就可以求得矩陣A。

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

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

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

(2)

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

令yi=lnxi

(3)

則,式(2)可以寫成:

y=A·ξ+μ

(4)

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

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

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

(5)

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

(6)

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

得到了正態(tài)分布y的A與μ,后續(xù)進(jìn)行拉丁超立方抽樣的計算步驟與正態(tài)分布的拉丁超立方抽樣一致。得到相關(guān)隨機(jī)變量yi的樣本后,通過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)的空間,這可通過正態(tài)分布相關(guān)隨機(jī)變量的拉丁超立方抽樣來解決。

原始空間向相關(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)已知,通過等概率轉(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個相關(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ù)值計算中,式(8)的雙重積分上下限選擇[-4,4]就足夠精確了[13]。在ρxixj、μxi、μxj已知的情況下,幾乎不可能直接計算ρyiyj,因此,需采用數(shù)值計算法求解ρyiyj。當(dāng)ρxixj>0時,具體步驟如下:

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

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

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

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

以巖土工程中土體常見的黏聚力c、內(nèi)摩擦角φ為例進(jìn)行拉丁超立方抽樣,抽取10 000個樣本。均值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)分布)、對數(shù)正態(tài)分布相關(guān)隨機(jī)變量抽樣(c、φ均滿足對數(shù)正態(tài)分布)、非正態(tài)分布相關(guān)隨機(jī)變量抽樣(c滿足對數(shù)正態(tài)分布,φ滿足正態(tài)分布),抽樣結(jié)果分布如圖2,抽樣參數(shù)c、φ概率密度如圖3。

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

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

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

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

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

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

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

Spearman秩相關(guān)系數(shù)只對變量的秩次大小進(jìn)行線性相關(guān)分析,所以沒有要求變量的分布,但秩相關(guān)系數(shù)與樣本的排序有關(guān),而與樣本具體的值關(guān)聯(lián)不大,因此,秩相關(guān)隨機(jī)變量的拉丁超立方抽樣可通過秩相關(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ī)變量具體的分布沒有關(guān)系,那么可以設(shè)ξ為標(biāo)準(zhǔn)正態(tài)分布,即Z為正態(tài)分布。

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

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

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

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

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

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

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

2.1 實(shí)例背景

目前,隧道工程仍以確定性分析為主,但其圍巖的材料參數(shù)往往具有不確定性,因此,有必要對某些情況進(jìn)行不確定性分析或可靠度設(shè)計[17]。不確定性分析方法可以是蒙特卡洛抽樣和數(shù)值計算方法的結(jié)合[18]。為了展示拉丁超立方抽樣在工程應(yīng)用中的優(yōu)越性,筆者用一個直徑D=6 m,埋深H=6 m的圓柱形隧道掌子面進(jìn)行失效概率分析。隧道周圍土體采用摩爾-庫倫模型,襯砌采用彈性模型。將土體的黏聚力c和內(nèi)摩擦角φ考慮為服從正態(tài)分布的隨機(jī)變量,兩個參數(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。隧道采用全斷面開挖,不考慮超前支護(hù)措施。

2.2 三維數(shù)值模型

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

2.3 模擬過程

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

N≥100/Pf

(10)

式中:Pf為失效概率。

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

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

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

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

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

2.4 結(jié)果分析

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

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

由圖7可見,當(dāng)樣本數(shù)量N=97 000個時,掌子面失效樣本的數(shù)量為100個,此時的失效概率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ù)量僅為蒙特卡洛抽樣樣本數(shù)的16.5%。由此可見,拉丁超立方抽樣更具代表性,可以顯著降低概率分析所用的樣本數(shù)量。

3 結(jié) 論

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

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

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