緱變彩, 張 天, 王 帆

(1. 武漢科技大學 城市學院， 湖北 武漢 430083； 2. 武漢工程大學 資源與土木工程學院， 湖北 武漢 430073)

土體抗剪強度參數(shù)具有較強不確定性，需要采取可靠性分析手段估計隧道、邊坡等巖土結(jié)構(gòu)的失效概率，從而為設(shè)計及施工提供決策依據(jù)。常用的可靠性分析方法包括FORM、SORM和蒙特卡羅模擬等[1]。

然而在實踐中采用上述方法進行可靠性分析時常常會遇到兩個難點。第一個難點在于隱性極限狀態(tài)方程的顯性化表達，通常采用響應(yīng)面法進行極限狀態(tài)曲面的刻畫，如采用二次多項式[2]、神經(jīng)網(wǎng)絡(luò)[3]、Kriging模型[4]等。然而傳統(tǒng)的響應(yīng)面法是一種侵入式算法，需要與原有數(shù)值模型進行聯(lián)合迭代運算。為提高計算效率，目前常采用隨機響應(yīng)面法構(gòu)建響應(yīng)面方程。相比傳統(tǒng)響應(yīng)面法，隨機響應(yīng)面[5]屬于非侵入式算法，其采用Hermite多項式在隨機空間內(nèi)對系統(tǒng)響應(yīng)進行擬合，因此基于隨機響應(yīng)面法構(gòu)建的替代模型可直接用于相關(guān)可靠性分析。

第二個難點在于概率信息的獲取。可靠性分析需要知道隨機變量的聯(lián)合概率分布函數(shù)。而實上，受限于試驗條件等難以準確獲取該聯(lián)合概率分布函數(shù)。相比之下，工程實踐中較易獲取的是隨機變量的邊緣分布及隨機變量之間的相關(guān)系數(shù)。當隨機變量相互獨立時，隨機變量的聯(lián)合概率分布函數(shù)為其邊緣分布函數(shù)的簡單乘積；但當隨機變量之間存在相關(guān)關(guān)系時，基于邊緣分布及相關(guān)系數(shù)不能唯一確定相應(yīng)的聯(lián)合概率分布函數(shù)[6]，因此該概率信息也被稱之為不完備概率信息[7]?；诓煌陚涓怕市畔⒌目煽啃苑治龀２捎肗ataf變換獲取標準隨機空間與原隨機空間的映射關(guān)系，如Li等[8]采用Nataf變換獲取隨機響應(yīng)面配點在原隨機空間內(nèi)的位置以構(gòu)建替代模型。但目前已證明Nataf變換實際上隱含假設(shè)隨機變量之間的相關(guān)結(jié)構(gòu)為高斯相關(guān)結(jié)構(gòu)[9]。換言之，其結(jié)果僅僅只是眾多可能結(jié)果中的一種，因此有必要研究當隨機變量之間的相關(guān)結(jié)構(gòu)為非高斯相關(guān)結(jié)構(gòu)時的巖土結(jié)構(gòu)可靠性。

本文提出了不完備概率信息下可靠性分析的隨機響應(yīng)面法，采用Copula理論建立隨機變量的聯(lián)合概率分布函數(shù)，基于Rosenblatt變換實現(xiàn)非高斯相關(guān)結(jié)構(gòu)下隨機響應(yīng)面在原隨機空間的配點選擇，最后通過兩個案例研究了不同相關(guān)結(jié)構(gòu)對可靠性的影響。

1 基于配點的隨機響應(yīng)面法

1.1 隨機響應(yīng)面法基本原理

隨機響應(yīng)面法一般采用Hermite多項式進行擬合。設(shè)Y是輸出隨機變量，X=[x1,x2, …,xn]是n維輸入隨機變量，U=[U1,U2, …,Un]是相互獨立的標準正態(tài)隨機向量，且X=T(U)表示U到X的等概率變換，則極限狀態(tài)Y=G(X)可表示為：

Y=G(X)=G(T(U))=H(U)

(1)

其中Y=H(U)為Hermite多項式：

(2)

式中：ai1i2…in為待定系數(shù)，Γp(·)為p階Hermite多項式：

Γp(Ui1,Ui2,...,Uip)

(3)

理論上，隨機響應(yīng)面的擬合精度隨著階數(shù)p的增加而增加；然而高階Hermite多項式將會導(dǎo)致展開項的快速增長。因此，通常需要對式(2)進行截斷。表1列舉了從2階到4階的Hermite多項式的展開形式。

表1 2～4階Hermite多項式的展開形式

Li等[8]總結(jié)了建立隨機響應(yīng)面模型的四個主要步驟：(1)通過獨立標準正態(tài)隨機變量表示輸入隨機變量；(2)使用Hermite多項式表示輸出隨機變量；(3)確定Hermite多項式的系數(shù)；(4)用蒙特卡羅模擬等方法估算失效概率。其中建立隨機響應(yīng)面模型的關(guān)鍵在于第三步即確定Hermite多項式的系數(shù)，其求解方法包括Galerkin法和配點法。相較于Galerkin法，配點法操作簡便，其有效性不受原響應(yīng)面的非線性和復(fù)雜性影響，且配點相互獨立，可同時計算每個配點對應(yīng)的系統(tǒng)響應(yīng)值。因此，本文采用配點法建立隨機響應(yīng)面模型。

1.2 配點的選擇

構(gòu)建隨機響應(yīng)面的關(guān)鍵在于求解待定系數(shù)ai1i2…in，可通過配點法確定，其基本原理是，對于p階Hermite多項式，采用p+1階Hermite多項式的根來確定標準隨機空間內(nèi)的輸入隨機變量，然后，通過等概率變換將這些配點從標準隨機空間轉(zhuǎn)換到原隨機空間中，再采用傳統(tǒng)結(jié)構(gòu)分析方法確定這些配點所對應(yīng)的系統(tǒng)響應(yīng)值，最后通過下式求解待定系數(shù)：

a=(HTH)-1HTF

(4)

式中：H和F分別為基于配點給出的Hermite多項式矩陣和配點對應(yīng)的系統(tǒng)響應(yīng)值所構(gòu)成的向量。

p階Hermite多項式的根的數(shù)量也為p個，而p階Hermite多項式的可選配點是p+1階Hermite多項式根的組合，因此對于n維問題，p階隨機響應(yīng)面的可選配點數(shù)量應(yīng)為(p+1)n。Mollon等[10]注意到對于奇數(shù)階隨機響應(yīng)面，原點不在可選配點之中，但考慮到原點在標準隨機空間中占據(jù)高概率密度區(qū)域，因此需要將原點也作為可選配點之一，于是可選配點的數(shù)量Ncp為：

(5)

如果選取所有的候選配點來確定Hermite多項式系數(shù)，那么當n和p較大時就必須進行大量的結(jié)構(gòu)分析來確定配點所對應(yīng)的系統(tǒng)響應(yīng)。實際上，確定系數(shù)所需的最小配點數(shù)目是：

(6)

一般地，Na遠小于Ncp。因此，合理選擇配點可以大大減少結(jié)構(gòu)分析的工作量。對于低維問題，可以通過確保H為滿秩矩陣的方法來選取配點[8]。

在確定Hermite多項式系數(shù)的過程中需要將配點從標準隨機空間映射到原隨機空間中，以便確定配點對應(yīng)的系統(tǒng)響應(yīng)。對于輸入隨機變量相互獨立的情況，可直接通過解析的方式進行表達[5]；對于輸入隨機變量之間存在相關(guān)關(guān)系且相關(guān)結(jié)構(gòu)為高斯相關(guān)結(jié)構(gòu)的情況，Li等[8]則給出了基于Nataf變換的映射方法；而對于不完備概率信息下非高斯相關(guān)結(jié)構(gòu)的相關(guān)隨機變量的配點映射將在下一節(jié)進行討論。

2 不完備概率信息下非高斯相關(guān)結(jié)構(gòu)隨機變量的配點映射

不完備概率信息條件下的配點映射其難點在于聯(lián)合概率分布模型的構(gòu)建。根據(jù)Sklar理論，一個二維聯(lián)合概率分布可以表示為該分布的邊緣分布和一個Copula函數(shù)：

F(x1,x2)=C(u1,u2,θ)

(7)

式中：u1=F(x1)和u2=F(x2)分別為隨機變量x1和x2的邊緣分布函數(shù)；C為Copula函數(shù)；θ是Copula函數(shù)的參數(shù)。

相應(yīng)的聯(lián)合概率分布函數(shù)如下所示：

=c(u1,u2,θ)f1(x1)f2(x2)

(8)

式中：f1(x1)和f2(x2)分別為x1和x2的邊緣概率密度函數(shù)，c(u1,u2,θ)是Copula密度函數(shù)。

當已知隨機變量間的Pearson相關(guān)系數(shù)ρ時，Copula參數(shù)θ可以用下式來確定：

(9)

可見，給定某一特定copula函數(shù)，則二維聯(lián)合概率分布函數(shù)可以基于其邊緣分布和相關(guān)系數(shù)構(gòu)造。特別的，當Copula函數(shù)取Normal時，所對應(yīng)的二維聯(lián)合概率分布函數(shù)為Nataf分布，而Nataf變換實際上就是假設(shè)聯(lián)合概率分布為Nataf分布下的等概率變換。

由于Nataf變換僅適用于高斯相關(guān)結(jié)構(gòu)的隨機變量，因此本文在建立聯(lián)合概率分布函數(shù)之后，采用Rosenblatt變換作為配點映射的方法。對于二維問題，Rosenblatt變換為：

(10)

式中：Φ-1(·)為標準正態(tài)分布累積分布函數(shù)Φ(·)的反函數(shù)；F(x2|x1)為邊緣條件分布，可由下式表示：

(11)

式中：h(u2,u1,θ)為條件Copula函數(shù)。表1列舉了一些常用Copula函數(shù)及相應(yīng)的條件Copula函數(shù)。因此，采用Rosenblatt變換的逆變換即可實現(xiàn)配點從標準隨機空間到原隨機空間的映射X=T(U)：

(12)

注：Φθ為相關(guān)系數(shù)為θ的二維標準高斯分布；Tθ,λ為相關(guān)系數(shù)為θ且自由度為λ的二維標準t分布；Tλ為自由度為λ的一維標準t分布

3 算例1

算例1研究如圖1所示的某一靜水壓力場下圓形隧道的開挖可靠性。假設(shè)巖體是同質(zhì)、均勻和連續(xù)的，靜水壓力場為p0，若開挖一半徑為R的隧道，則塑性區(qū)半徑Rpl為：

(13)

(14)

(15)

(16)

式中：c和φ分別為粘聚力和內(nèi)摩擦角；pi為支護壓力。

假設(shè)p0=2.5 Mpa，pi=0.1 Mpa，c和φ為隨機變量，其統(tǒng)計分布如表3所示，并假設(shè)c和φ之間的相關(guān)系數(shù)ρ=-0.5。

表3 算例1 隨機變量的分布

圖1 隧道開挖示意

假設(shè)允許的塑性區(qū)半徑為3倍洞徑，則功能函數(shù)為：

(17)

首先假設(shè)c和φ之間的相關(guān)結(jié)構(gòu)為高斯相關(guān)結(jié)構(gòu)即Normal Copula，分別建立2～4階隧道開挖塑性區(qū)半徑的隨機響應(yīng)面，然后基于該隨機響應(yīng)面通過蒙特卡洛模擬(N=105)得出塑性區(qū)半徑的分布。圖2表示的是塑性區(qū)半徑的累積概率函數(shù)(圖中：y為允許的相對塑性區(qū)半徑)，為便于驗證本方法，同時將基于Nataf變換并直接采用蒙特卡洛方法(N=105)得出的結(jié)果進行比較。注意建立2～4階隨機響應(yīng)面時采用原解析模型進行隧道塑性區(qū)半徑計算的次數(shù)分別為6，10，15次，其后蒙特卡洛模擬只是簡單的代數(shù)計算，而直接采用蒙特卡洛模擬需要調(diào)用原模型105次。

圖2 隧道開挖塑性區(qū)半徑的累計概率函數(shù)

表4比較了不同方法下基于式17的隧道開挖失效概率。若以Nataf變換加直接采用蒙特卡洛方法得出的結(jié)果作為精確解，則3階隨機響應(yīng)面所對應(yīng)的結(jié)果已具備較高精度。

表4 高斯相關(guān)結(jié)構(gòu)下可靠性分析結(jié)果對比

上述結(jié)果是假設(shè)c和φ之間的相關(guān)結(jié)構(gòu)為高斯相關(guān)結(jié)構(gòu)，然而實際上該假設(shè)并不一定成立，為檢驗非高斯相關(guān)結(jié)構(gòu)是否對隧道開挖可靠性存在較大影響，本文分別研究了當相關(guān)結(jié)構(gòu)用t Copula (λ=2)和Frank Copula表示時的隧道開挖失效概率。

圖3對比了基于3階隨機響應(yīng)面的不同相關(guān)結(jié)構(gòu)下隧道塑性區(qū)半徑分布的累積概率函數(shù)曲線?？梢园l(fā)現(xiàn)，不同相關(guān)結(jié)構(gòu)對可靠性分析結(jié)果有著顯著影響。表5比較了具體的失效概率以及相較于高斯相關(guān)結(jié)構(gòu)時非高斯相關(guān)結(jié)構(gòu)對應(yīng)失效概率的偏差。可以看出當分別假設(shè)相關(guān)結(jié)構(gòu)需要用t Copula (λ=2)和Frank Copula來表示時，基于高斯相關(guān)結(jié)構(gòu)假設(shè)得出的隧道失效概率會明顯高于或低于實際的失效概率，相對偏差分別為-10.61%和6.76%。因此，基于高斯相關(guān)結(jié)構(gòu)假設(shè)有可能導(dǎo)致過于保守或過于冒進的施工及設(shè)計方案。

圖3 隧道開挖塑性區(qū)半徑的累計概率函數(shù)

Copula(相關(guān)結(jié)構(gòu))Normal(高斯)tCopula(λ=2)Frank失效概率p(G<0)0.02810.02540.0300偏差/%—-9.616.76

4 算例2

算例2研究如圖4所示的某1∶1邊坡，邊坡高度10 m。土體抗剪強度參數(shù)即有效粘聚力c和內(nèi)摩擦角φ均為隨機變量，且服從表6給出的分布，相關(guān)系數(shù)ρ=-0.458，土體重度假設(shè)為常數(shù)。

圖4 邊坡穩(wěn)定性示意

參數(shù)分布類型形狀參數(shù)比例參數(shù)c/kPaWeibull1.8658.24φ/(°)Gamma16.981.19

邊坡穩(wěn)定性分析是典型的隱性極限狀態(tài)函數(shù)問題，其穩(wěn)定性通常用安全系數(shù)FS表示，含義為潛在滑動面上的土體抗剪強度與滑動荷載的比值。邊坡穩(wěn)定性分析中極限狀態(tài)函數(shù)為隱性，其原因在于最小安全系數(shù)對應(yīng)的滑動面未知，需要通過一定的搜索算法進行識別，而且土體抗剪強度參數(shù)變化時，該潛在滑動面的位置也可能發(fā)生變化。

這里采用Bishop法計算邊坡安全系數(shù)，滑動面形狀假設(shè)為圓弧形，通過網(wǎng)格搜索確定給定土體抗剪強度下邊坡安全系數(shù)的最小值及其對應(yīng)的滑動面。為驗證所提出方法的有效性，基于Rosenblatt變換通過蒙特卡洛(N=104)產(chǎn)生不同相關(guān)結(jié)構(gòu)下服從已知邊緣分布和相關(guān)系數(shù)的土體抗剪強度參數(shù)樣本，然后直接采用邊坡穩(wěn)定性計算模型算出樣本所對應(yīng)的安全系數(shù)的分布，之后采用隨機響應(yīng)面法計算邊坡的失效概率并與蒙特卡洛方法進行對比。

表7給出了該邊坡可靠性分析的對比結(jié)果，基于隨機響應(yīng)面法得出的失效概率與蒙特卡洛所得出的結(jié)果基本一致。但是隨機響應(yīng)面法對原解析模型調(diào)用的次數(shù)遠低于直接采用蒙特卡洛法調(diào)用原模型的次數(shù)。在本例中，基于5階隨機響應(yīng)面對原模型的調(diào)用次數(shù)僅為21次，之后邊坡安全系數(shù)的計算只需通過Hermite多項式進行簡單的代數(shù)計算即可，總耗時約15 s左右，而采用樣本數(shù)為10000的蒙特卡洛模擬則需要1小時以上。相關(guān)結(jié)構(gòu)對于邊坡失效概率仍有較大影響。當層狀土有效粘聚力相關(guān)結(jié)構(gòu)為Frank時，若仍基于Gauss相關(guān)結(jié)構(gòu)(Nataf變換)進行分析會明顯低估邊坡失效概率，從而造成過于經(jīng)濟的設(shè)計。

表7 邊坡可靠性分析結(jié)果對比

5 討 論

為研究相關(guān)結(jié)構(gòu)如何影響可靠性分析結(jié)果，以案例1為例，圖5a，5b分別顯示了不同相關(guān)結(jié)構(gòu)下3階隨機響應(yīng)面配點在標準隨機空間和原隨機空間的位置(圖中，U1，U2分別表示標準正態(tài)空間中的c和φ)?？梢钥闯鲈跇藴孰S機空間中，候選配點和基于滿秩矩陣實際選取的配點是相同的，不依賴于隨機變量之間的相關(guān)系數(shù)和相關(guān)結(jié)構(gòu)；但是當配點從標準隨機空間向原隨機空間進行映射時，由于c和φ之間存在負相關(guān)關(guān)系，因此在原隨機空間中配點位置相較于標準隨機空間進行了旋轉(zhuǎn)，但不同相關(guān)結(jié)構(gòu)(Copula)對該旋轉(zhuǎn)有著一定影響，除原點外的所有配點在旋轉(zhuǎn)時均有不同程度的偏移，而且配點距離原點越遠，不同相關(guān)結(jié)構(gòu)對于該配點旋轉(zhuǎn)后位置的影響越大。因此，在假設(shè)不同相關(guān)結(jié)構(gòu)時，同樣的配點最后對應(yīng)的系統(tǒng)響應(yīng)值會有所不同，進而導(dǎo)致了隨機響應(yīng)面待定系數(shù)發(fā)生改變，從而得出不同的可靠性分析結(jié)果。

圖5 3階隨機響應(yīng)面配點的位置

由于相關(guān)結(jié)構(gòu)可能對可靠性結(jié)果產(chǎn)生明顯的影響，選擇合適的Copula函數(shù)來表征隨機變量之間的相關(guān)結(jié)構(gòu)顯得尤為重要。然而，在缺少實測數(shù)據(jù)的情況下難以準確確定變量之間的相關(guān)結(jié)構(gòu)。因此，在給定隨機變量的邊緣分布和相關(guān)系數(shù)的條件下，需要補充有關(guān)相關(guān)結(jié)構(gòu)信息。例如，若已知變量之間存在非對稱的下尾相關(guān)性，則可以采用Clayton Copula計算失效概率，而不必受限于高斯相關(guān)結(jié)構(gòu)假設(shè)的限制。

如果缺少相關(guān)結(jié)構(gòu)的信息，可以檢驗不同非高斯相關(guān)結(jié)構(gòu)相比高斯相關(guān)結(jié)構(gòu)對應(yīng)的可靠性結(jié)果偏差。如果偏差較大，再補充數(shù)據(jù)用于準確界定變量之間的相關(guān)結(jié)構(gòu)或者直接選取較為保守的設(shè)計值。

6 結(jié) 論

本文研究了不完備概率信息下可靠性分析的隨機響應(yīng)面方法，基于Copula理論建立了聯(lián)合概率分布模型，通過Rosenblatt變換實現(xiàn)了配點從標準隨機空間到原隨機空間的映射，并以某圓形隧道開挖和某層狀邊坡為例探討了不同相關(guān)結(jié)構(gòu)對隨機響應(yīng)面構(gòu)建及可靠性結(jié)果的影響，得出以下結(jié)論：

(1)當假設(shè)不完備概率信息條件下變量之間的相關(guān)結(jié)構(gòu)為高斯相關(guān)結(jié)構(gòu)時，基于Nataf變換和基于Ronsenblatt變換得出的可靠性結(jié)論是一致的，說明兩種等概率變換形式在假設(shè)Normal Copula時是等效的；

(2)隨機響應(yīng)面可以較為精確地將原有隱性功能函數(shù)顯性化，但不完備概率信息下隨機變量之間的相關(guān)結(jié)構(gòu)未知，而假設(shè)不同相關(guān)結(jié)構(gòu)對于可靠性結(jié)果有著較為顯著的影響，常用的高斯相關(guān)結(jié)構(gòu)會明顯低估或高估系統(tǒng)的失效概率進而導(dǎo)致過于保守或過于經(jīng)濟的設(shè)計和施工方案；

(3)不同相關(guān)結(jié)構(gòu)即Copula選擇對于可靠性分析結(jié)果的影響是通過改變配點在原隨機空間中的位置來實現(xiàn)的，除原點外的所有配點在映射后均有不同程度的偏移，而且配點距離原點越遠，不同相關(guān)結(jié)構(gòu)對于該配點旋轉(zhuǎn)后位置的影響越大，從而導(dǎo)致相同的配點卻對應(yīng)不同的系統(tǒng)響應(yīng)值并影響可靠性分析結(jié)果。

[1] Ditlevsen O, Madsen H O. Structural Reliability Methods[M]. Chichester : John Wiley & Sons Ltd, 2007.

[2] Xu B, Low B K. Probabilistic stability analyses of embankments based on finite element method[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(11): 1444-1454.

[3] Cho S E. Probabilistic stability analyses of slopes using the ANN-based response surface[J]. Computers and Geotechnics, 2009, 36(5): 787-797.

[4] Luo X, Li X, Zhou J, et al. A Kriging-based hybrid optimization algorithm for slope reliability analysis[J]. Structural Safety, 2012, 34(1): 401-406.

[5] Huang S P, Liang B, Phoon K K. Geotechnical probabilistic analysis by collocation-based stochastic response surface method: An Excel add-in implementation[J]. Georisk: Assessment and Management of Risk for Engineered Systems and Geohazards, 2009, 3(2): 75-86.

[6] Tang X S, Li D Q, Rong G, et al. Impact of copula selection on geotechnical reliability under incomplete probability information[J]. Computers and Geotechnics, 2013, 49: 264-278.

[7] Der Kiureghian A, Liu P L. Structural reliability under incomplete probability information[J]. Journal of Engineering Mechanics, 1986, 112(1): 85-104.

[8] Li D, Chen Y, Lu W, et al. Stochastic response surface method for reliability analysis of rock slopes involving correlated non-normal variables[J]. Computers and Geotechnics, 2011, 38(1): 58-68.

[9] Lebrun R, Dutfoy A. An innovating analysis of the Nataf transformation from the copula viewpoint[J]. Probabilistic Engineering Mechanics, 2009, 24(3): 312-320.

[10]Mollon G, Dias D, Soubra A H. Probabilistic analysis of pressurized tunnels against face stability using collocation-based stochastic response surface method[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2011, 137(4): 385-397.