緱變彩, 王 帆

(1. 武漢科技大學(xué)城市學(xué)院 城建學(xué)部, 湖北 武漢 430083； 2. 華中科技大學(xué) 土木與水利工程學(xué)院， 湖北 武漢 430074)

土體作為典型的天然材料具有非勻質(zhì)性、各向異性，呈現(xiàn)出空間變異特征，而土的工程特性依賴于土體的空間平均特性，因此在進(jìn)行巖土結(jié)構(gòu)變形與穩(wěn)定性分析時(shí)需要考慮土體參數(shù)的空間變異性對(duì)分析結(jié)果的影響[1]。邊坡作為典型的巖土結(jié)構(gòu)在進(jìn)行穩(wěn)定性分析時(shí)也需要對(duì)土體參數(shù)的空間變異性進(jìn)行充分模擬。

當(dāng)前針對(duì)空間變異土體的邊坡穩(wěn)定性問題國(guó)內(nèi)外已有諸多研究成果[2～5]，其中分析的重點(diǎn)和難點(diǎn)主要體現(xiàn)在兩個(gè)方面，即土體參數(shù)空間變異性的模擬和邊坡失穩(wěn)概率的高效計(jì)算。土體參數(shù)的空間變異模擬通?；陔S機(jī)場(chǎng)理論，通過均值函數(shù)、方差函數(shù)等描述拓?fù)淇臻g中任意一點(diǎn)的隨機(jī)特性，通過自相關(guān)函數(shù)描述拓?fù)淇臻g中任意兩點(diǎn)之間的相似性[6]。隨機(jī)場(chǎng)離散是實(shí)現(xiàn)土體參數(shù)空間變異性模擬的關(guān)鍵步驟之一，目前常用的離散方法可分為空間離散法和級(jí)數(shù)展開法[7]。空間離散法包括中點(diǎn)法、形函數(shù)法、局部平均法等，空間離散法的誤差精度與隨機(jī)場(chǎng)網(wǎng)格劃分有關(guān)，對(duì)于自相關(guān)長(zhǎng)度較小的情況需要較為精細(xì)的隨機(jī)場(chǎng)網(wǎng)格，導(dǎo)致計(jì)算效率較低。級(jí)數(shù)展開法是將隨機(jī)場(chǎng)表示為級(jí)數(shù)展開形式，通過截?cái)嗾归_項(xiàng)近似表達(dá)原有隨機(jī)場(chǎng)。Karhunen-Loève(K-L)級(jí)數(shù)展開是其中較常用的一種，相較于其他方法，K-L法從隨機(jī)場(chǎng)離散誤差和離散后輸入變量空間維度的角度講是最優(yōu)的[8,9]，因此在邊坡可靠性分析中廣泛用于土性空間變異性模擬。

邊坡穩(wěn)定性分析方法包括極限平衡法、有限元法和極限分析法等，采用一階/二階可靠度、點(diǎn)估計(jì)、蒙特卡洛等方法估計(jì)邊坡失穩(wěn)概率[10～12]。然而這些方法一般需要進(jìn)行多次穩(wěn)定性分析，尤其是穩(wěn)定性分析模型較復(fù)雜(如有限元)或小概率問題時(shí)計(jì)算效率較低。針對(duì)這一問題，采用代理模型是一種有效的解決辦法，即通過實(shí)驗(yàn)設(shè)計(jì)進(jìn)行有限次數(shù)的分析獲取輸入輸出樣本，再通過訓(xùn)練學(xué)習(xí)得到新模型并用于可靠度分析。常用代理模型包括二次響應(yīng)面、克里金法、人工神經(jīng)網(wǎng)絡(luò)等[13,14]。其中多項(xiàng)式混沌展開(Polynomial Chaos Expansion，PCE)因其嚴(yán)格的數(shù)學(xué)背景，可表達(dá)任意有限方差隨機(jī)響應(yīng)，且對(duì)于平滑的輸入輸出關(guān)系收斂迅速，因此得到廣泛應(yīng)用。該方法的應(yīng)用可分為侵入式和非侵入式兩種，如采用隨機(jī)伽遼金法求解PCE系數(shù)的譜隨機(jī)有限元屬于侵入式方法[15]，而采用隨機(jī)配點(diǎn)法求解PCE系數(shù)則屬于非侵入式方法[16]。通常而言，侵入式方法(如隨機(jī)伽遼金法)相比非侵入式方法(如隨機(jī)配點(diǎn)法)精度更高，但由于需要對(duì)原計(jì)算模型進(jìn)行修改，使用過程較為復(fù)雜，甚至對(duì)于某些復(fù)雜非線性系統(tǒng)其對(duì)應(yīng)的PCE方程難以推導(dǎo)；而非侵入式PCE模型的構(gòu)建與求解與原計(jì)算模型是解耦合的，因此適用范圍更廣?？紤]土性空間變異的邊坡穩(wěn)定性分析屬于典型的非線性問題，且輸入空間維度較高，而PCE展開項(xiàng)和未知系數(shù)的個(gè)數(shù)會(huì)隨PCE階數(shù)和輸入空間維度的增加迅速增長(zhǎng)，導(dǎo)致系數(shù)求解所需要的原有模型計(jì)算成本過高。

為解決這一問題，Blatman[17]提出了一種基于稀疏多項(xiàng)式混沌展開(sparse PCE，sPCE)的代理模型方法。該方法利用稀疏效應(yīng)準(zhǔn)則，采用壓縮感知算法重構(gòu)稀疏的展開式，即非零系數(shù)項(xiàng)，運(yùn)用基于留一誤差的回歸算法求解有限的展開式系數(shù)，實(shí)現(xiàn)了求解精度與計(jì)算效率之間的平衡。本研究采用sPCE方法評(píng)估土性空間變異邊坡的失穩(wěn)概率，分析不同輸入變量的全局敏感度。為實(shí)現(xiàn)這一目標(biāo)，采用PCE逼近理論構(gòu)造邊坡穩(wěn)定性代理模型，利用正交匹配追蹤[18]實(shí)現(xiàn)PCE稀疏化重構(gòu)，結(jié)合拉丁超立方抽樣建立目標(biāo)精度下的PCE系數(shù)自適應(yīng)求解[19]。該方法分別用于分析不同隨機(jī)場(chǎng)模型和邊坡穩(wěn)定性計(jì)算模型的邊坡失穩(wěn)概率，為高效評(píng)估邊坡可靠性提供思路。

1 土性隨機(jī)場(chǎng)離散的K-L法

設(shè)H(z)是一隨機(jī)場(chǎng)，z∈Ω是d維拓?fù)淇臻gΩ?Rd中的一點(diǎn)。另設(shè)隨機(jī)場(chǎng)均值函數(shù)為μ(z)，自協(xié)方差函數(shù)為C(z,z′)，則該隨機(jī)場(chǎng)可展開為如下表達(dá)式：

(1)

式中：ξi為一系列零均值、單位方差且不相關(guān)的隨機(jī)變量；λi和φi(z)分別為自協(xié)方差函數(shù)的特征值和特征函數(shù)，即求解如下廣義特征值問題：

(2)

為便于計(jì)算，需要對(duì)K-L級(jí)數(shù)進(jìn)行截?cái)啵瑑H保留有限項(xiàng)特征值及特征函數(shù)來近似表達(dá)原隨機(jī)場(chǎng)，其對(duì)應(yīng)的均方誤差為：

(3)

通常需要對(duì)特征值從大到小進(jìn)行排序，取前m項(xiàng)從而使得在允許誤差下K-L級(jí)數(shù)展開項(xiàng)數(shù)最少。

對(duì)于高斯隨機(jī)場(chǎng)，ξi為獨(dú)立標(biāo)準(zhǔn)高斯隨機(jī)變量；而對(duì)于非高斯隨機(jī)場(chǎng)，ξi通常難以獲得。因此K-L級(jí)數(shù)展開主要用于表達(dá)高斯隨機(jī)場(chǎng)，而對(duì)于非高斯隨機(jī)場(chǎng)則采用場(chǎng)轉(zhuǎn)換理論將K-L級(jí)數(shù)展開表達(dá)的高斯場(chǎng)轉(zhuǎn)換為對(duì)應(yīng)的非高斯場(chǎng)。由于土性參數(shù)如抗剪強(qiáng)度通常是非負(fù)的，理論上不適合采用高斯場(chǎng)模擬，因此實(shí)際中常假設(shè)為對(duì)數(shù)正態(tài)隨機(jī)場(chǎng)，并采用下式進(jìn)行轉(zhuǎn)換：

HLognormal(z)=exp(HGauss(z))

(4)

2 稀疏多項(xiàng)式混沌展開方法

2.1 多項(xiàng)式混沌展開

PCE逼近理論的基本思想是用一組正交多項(xiàng)式之和近似表達(dá)概率空間中的隱式函數(shù)。設(shè)響應(yīng)量y與n維獨(dú)立隨機(jī)變量x=[x1,x2…,xn]之間存在某函數(shù)關(guān)系y=f(x)，則y可以展開為一組正交基函數(shù)：

(5)

=γαδαβ

(6)

2.2 展開式系數(shù)求解

建立PCE代理模型的關(guān)鍵在于展開式系數(shù)求解。為便于計(jì)算，需要對(duì)展開式進(jìn)行截?cái)啵∽畲箅A數(shù)p以內(nèi)的展開項(xiàng)，即|α|≤p對(duì)應(yīng)的基函數(shù)。對(duì)于非侵入式模型構(gòu)建，常用的系數(shù)求解方法有投影法、配點(diǎn)法、回歸法等。相較而言，在調(diào)用相同原計(jì)算模型次數(shù)的情況下，回歸法的收斂速度更快，因此這里采用回歸法進(jìn)行求解。首先通過實(shí)驗(yàn)設(shè)計(jì)得到輸入空間內(nèi)的一組樣本X=[x(1),x(2),…,x(N)]T，然后調(diào)用原計(jì)算模型得到對(duì)應(yīng)的響應(yīng)量Y=[y(1),y(2),…,y(N)]T，則展開式系數(shù)為：

ηα=(ΨTΨ)-1ΨTY

(7)

若對(duì)于n維輸入變量采用p階PCE擬合其響應(yīng)關(guān)系，則展開項(xiàng)總項(xiàng)數(shù)P為：

(8)

可以看到PCE展開總項(xiàng)數(shù)隨著輸入空間維度和PCE階數(shù)快速增長(zhǎng)，而實(shí)際中為確保系數(shù)求解的精度，通常實(shí)驗(yàn)設(shè)計(jì)采樣次數(shù)需要k≥2P，因此對(duì)于高維問題或高階PCE展開則會(huì)導(dǎo)致原模型計(jì)算次數(shù)顯著增加。在某些情況下，即使采用K-L級(jí)數(shù)展開表達(dá)隨機(jī)場(chǎng)，其輸入變量維度也較高，例如土性空間變異性通常采用在原點(diǎn)處不可微的指數(shù)型自相關(guān)函數(shù)，且在垂直方向上的自相關(guān)長(zhǎng)度較小，從而導(dǎo)致離散誤差收斂速度較慢，需要較多的獨(dú)立標(biāo)準(zhǔn)高斯隨機(jī)變量去近似表達(dá)離散隨機(jī)場(chǎng)模型。因此需要采用更有效的方法估計(jì)PCE展開式系數(shù)。

2.3 正交多項(xiàng)式基函數(shù)的稀疏化選取

根據(jù)稀疏效應(yīng)準(zhǔn)則，實(shí)際問題中占主導(dǎo)的展開式通常是僅包含單獨(dú)輸入變量的多項(xiàng)式基函數(shù)和變量之間的低階交叉項(xiàng)，而其他多項(xiàng)式對(duì)應(yīng)的系數(shù)趨近于零，對(duì)結(jié)果影響較小，可以忽略。因此，可以通過合理選取多項(xiàng)式基函數(shù)提高PCE模型建立的效率。

基函數(shù)選取的關(guān)鍵在于如何確保所構(gòu)建的代理模型其泛化誤差是足夠小的。采用訓(xùn)練誤差判斷模型精度容易造成過擬合現(xiàn)象，因此不適宜作為選擇標(biāo)準(zhǔn)。為解決這一問題，通常采用交叉檢驗(yàn)的方法進(jìn)行誤差判斷。特別地，當(dāng)交叉檢驗(yàn)每次僅選取一個(gè)樣本作為驗(yàn)證集則對(duì)應(yīng)留一法(Leave-one-out，LOO)，相應(yīng)的LOO誤差為：

(9)

(10)

式中：hi是矩陣Ψ(ΨTΨ)-1ΨT的第i個(gè)對(duì)角元素。

上述LOO誤差為基函數(shù)的選取提供了一種有效的評(píng)價(jià)指標(biāo)。常用的基函數(shù)選取方法包括貪心算法和凸松弛算法等，其中貪心算法給出的是L0范數(shù)(非零PCE系數(shù)的個(gè)數(shù))最小化問題的近似解，而凸松弛算法則是將上述L0范數(shù)最小化問題轉(zhuǎn)化為L(zhǎng)1范數(shù)(各PCE系數(shù)絕對(duì)值之和)最小化問題后求解[20]。本文采用貪心算法中的正交匹配追蹤算法進(jìn)行重要基函數(shù)的選取[21]，該算法通過不斷驗(yàn)證殘差和備選基函數(shù)之間的正交性確定最似非零系數(shù)的基函數(shù)，然后再以式(7)求解當(dāng)前選取的基函數(shù)對(duì)應(yīng)的展開項(xiàng)系數(shù)。具體算法如下：

(3)更新已選基函數(shù)集Ak=Ak-1∪jk和候選基函數(shù)集Ck=Ck-1jk。

(4)用式(7)求解當(dāng)前已選基函數(shù)Ψα,k={Ψαj}j∈Ak的系數(shù)ηα,k，并更新殘差rk=Y-Ψα,kηα,k。

(5)用式(9)計(jì)算當(dāng)前PCE代理模型的LOO誤差εLOO,k。

(6)重復(fù)(2)～(5)直至k=min(N,P)。

3 考慮土性空間變異的邊坡可靠度與全局敏感性分析

基于上述K-L級(jí)數(shù)展開和sPCE代理模型可實(shí)現(xiàn)土性空間變異的邊坡可靠度高效評(píng)估，并分析不同變量的全局敏感性。具體流程如圖1所示。

圖1 分析流程

(2)對(duì)土性隨機(jī)場(chǎng)進(jìn)行K-L分解，求得截?cái)嗾`差內(nèi)的K-L表達(dá)。

(3)采用拉丁超立方抽樣獲取實(shí)驗(yàn)設(shè)計(jì)樣本。

(4)建立邊坡穩(wěn)定性計(jì)算模型，計(jì)算樣本對(duì)應(yīng)的邊坡安全系數(shù)。

(5)通過正交匹配追蹤選取多項(xiàng)式基函數(shù)。

(7)依據(jù)構(gòu)建的sPCE代理模型，采用蒙特卡洛法估計(jì)邊坡失穩(wěn)概率；計(jì)算不同輸入變量的Sobol指數(shù)STi：

(11)

4 算例分析

4.1 一維隨機(jī)場(chǎng)

圖2 邊坡示意/m

圖2為某1∶2邊坡，設(shè)土體不排水抗剪強(qiáng)度隨深度方向的空間變異性可用平穩(wěn)對(duì)數(shù)正態(tài)隨機(jī)場(chǎng)Su(z)模擬，均值為40 kPa，標(biāo)準(zhǔn)差為10 kPa，自相關(guān)函數(shù)為ρ(z,z′)=exp(-|z-z′|/l)，l=3 m為自相關(guān)長(zhǎng)度，土體重度為γ=18 kN/m3。邊坡穩(wěn)定性計(jì)算常采用極限平衡法或有限元法?？紤]到邊坡失效概率較小，采用有限元法驗(yàn)證時(shí)計(jì)算成本較高，故此處邊坡穩(wěn)定性計(jì)算采用極限平衡法中的Bishop法，滑移面假設(shè)為圓弧型。由于本文提出的sPCE代理模型構(gòu)建屬于非侵入式方法，因此將極限平衡法替換為有限元法構(gòu)建sPCE代理模型也是完全可行的。

設(shè)K-L截?cái)嘣试S誤差為5%，則需要前21階特征值和特征函數(shù)。允許的LOO誤差設(shè)為1×10-8，初始實(shí)驗(yàn)設(shè)計(jì)樣本大小和每次擴(kuò)充樣本大小均為20。

圖3給出了sPCE代理模型構(gòu)建的LOO誤差演化，經(jīng)過11次迭代降低到允許誤差范圍內(nèi)。圖4給出了該誤差下構(gòu)建的sPCE模型在訓(xùn)練集上的結(jié)果對(duì)比。

圖3 LOO誤差變化

圖4 sPCE模型在訓(xùn)練集上的結(jié)果對(duì)比

為驗(yàn)證模型的有效性，以蒙特卡洛仿真(Monte Carlo Simulation，MCS)50000次的原模型計(jì)算結(jié)果作為參考，失效概率為9×10-4，估計(jì)誤差為±14.9%。圖5給出了基于構(gòu)建的sPCE模型預(yù)測(cè)結(jié)果與實(shí)際計(jì)算值的對(duì)比。圖6給出了兩種方法求得的邊坡安全系數(shù)的概率密度函數(shù)，基于sPCE模型求得的邊坡安全系數(shù)的分布與MCS方法一致。然而MCS方法共調(diào)用原邊坡穩(wěn)定性分析模型5萬次，耗時(shí)約4 h，而采用sPCE方法僅調(diào)用了220次。另外，建立的代理模型僅包含218項(xiàng)多項(xiàng)式基函數(shù)，僅為完整PCE展開項(xiàng)數(shù)的0.07%。

圖5 sPCE模型在測(cè)試集上的結(jié)果對(duì)比

圖6 邊坡安全系數(shù)概率密度函數(shù)對(duì)比(算例1)

基于sPCE模型可以高效地計(jì)算安全系數(shù)的統(tǒng)計(jì)矩和失效概率。表1給出了不同允許誤差下安全系數(shù)的均值與方差，可以看到，隨著允許誤差的降低，需要更高階的PCE展開項(xiàng)來減小擬合誤差，相應(yīng)地模型給出的統(tǒng)計(jì)量和失效概率逐漸逼近MCS方法得到的結(jié)果，模型精度隨LOO允許誤差的降低而上升。

表1 不同允許誤差下的可靠性分析結(jié)果(算例1)

圖7給出了不同輸入變量的全局敏感性，且Sobol指數(shù)的總和略大于1，可見不同特征模式的交互項(xiàng)對(duì)全局敏感性影響較小。影響邊坡安全系數(shù)變化的主要是前4階特征模式，圖8給出了前4階模式的特征值及其特征函數(shù)。可以看到，由于邊坡安全系數(shù)取決于土體的空間平均特性，因此對(duì)于不排水抗剪強(qiáng)度在局部的空間變化(即高階特征模式)并不敏感。

圖7 不排水抗剪強(qiáng)度參數(shù)空間變異特征模式的敏感性

圖8 前4階特征值和特征函數(shù)

4.2 二維隨機(jī)場(chǎng)

設(shè)K-L截?cái)嘣试S誤差為5%，則需要前82階特征值和特征函數(shù)。允許的LOO誤差設(shè)為1×10-8，初始實(shí)驗(yàn)設(shè)計(jì)樣本大小和每次擴(kuò)充樣本大小均為10。該問題屬于典型的高維建模，對(duì)于4階sPCE模型，備選基函數(shù)的數(shù)目高達(dá)212萬，因此從完整的sPCE展開基函數(shù)中進(jìn)行選取效率十分低下。因此基于稀疏效應(yīng)原則，僅將包含單變量的基函數(shù)及包含兩個(gè)不同變量的低階交互項(xiàng)作為備選基函數(shù)，使得相應(yīng)的備選基函數(shù)數(shù)量降低至3650項(xiàng)。

表2給出了不同允許誤差下的分析結(jié)果，在允許誤差為1×10-8時(shí)需要構(gòu)建4階模型來達(dá)到設(shè)定的誤差范圍，失效概率估計(jì)為5.8×10-2。為驗(yàn)證模型的準(zhǔn)確性，將代理模型得到的邊坡安全系數(shù)的分布與1000次蒙特卡洛仿真結(jié)果進(jìn)行比較，從圖10可以看到兩者基本一致，蒙特卡洛方法估計(jì)的失效概率為6.0×10-2，估計(jì)誤差為±12.5%。

圖10 邊坡安全系數(shù)概率密度函數(shù)對(duì)比(算例2)

表2 不同允許誤差下的可靠性分析結(jié)果(算例2)

圖11給出了不同輸入變量的全局敏感性， Sobol指數(shù)的總和約等于1，說明不同特征模式的交互項(xiàng)對(duì)全局敏感性影響可以忽略。影響邊坡安全系數(shù)變化的主要是前2階特征模式，而對(duì)高階特征模式并不敏感，圖12給出了對(duì)應(yīng)的特征值及其特征函數(shù)。

圖11 不排水抗剪強(qiáng)度參數(shù)空間變異特征模式的敏感性

圖12 前2階特征值和特征函數(shù)

5 結(jié) 論

考慮邊坡土性空間變異性時(shí)輸入空間維度通常較高，采用多項(xiàng)式混沌展開逼近理論擬合響應(yīng)量其多項(xiàng)式基函數(shù)會(huì)快速增長(zhǎng)，給代理模型的建立提出了挑戰(zhàn)。采用正交匹配追蹤算法可以自動(dòng)選擇重要的基函數(shù)，實(shí)現(xiàn)模型稀疏化構(gòu)建，減少原計(jì)算模型的調(diào)用次數(shù)。

邊坡可靠性分析算例表明，對(duì)于K-L離散土性隨機(jī)場(chǎng)帶來的高維輸入空間問題，僅需較少次調(diào)用原邊坡穩(wěn)定性分析模型，即可得到較為精確的邊坡安全系數(shù)概率分布及其統(tǒng)計(jì)矩，構(gòu)建的代理模型中非零系數(shù)多項(xiàng)式基函數(shù)遠(yuǎn)小于完整多項(xiàng)式混沌展開的項(xiàng)數(shù)，且主要是僅包含單變量的基函數(shù)及包含兩個(gè)不同變量的低階交互項(xiàng)，大量的高階交互項(xiàng)對(duì)擬合邊坡安全系數(shù)并沒起到作用。

全局敏感性分析表明，由于邊坡安全系數(shù)取決于土體的空間平均特性，因此影響邊坡安全系數(shù)變化的主要是低階特征模式，而對(duì)于高階特征模式并不敏感，因此在構(gòu)建邊坡穩(wěn)定性分析的多項(xiàng)式混沌展開代理模型時(shí)還可以進(jìn)一步簡(jiǎn)化，對(duì)于K-L展開僅截取其低階特征值和特征函數(shù)。