基于稀疏混沌多項式的邊坡可靠度分析方法

緱變彩, 王 帆

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

土體作為典型的天然材料具有非勻質性、各向異性，呈現出空間變異特征，而土的工程特性依賴于土體的空間平均特性，因此在進行巖土結構變形與穩定性分析時需要考慮土體參數的空間變異性對分析結果的影響[1]。邊坡作為典型的巖土結構在進行穩定性分析時也需要對土體參數的空間變異性進行充分模擬。

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

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

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

1 土性隨機場離散的K-L法

設H(z)是一隨機場，z∈Ω是d維拓撲空間Ω?Rd中的一點。另設隨機場均值函數為μ(z)，自協方差函數為C(z,z′)，則該隨機場可展開為如下表達式：

(1)

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

(2)

為便于計算，需要對K-L級數進行截斷，僅保留有限項特征值及特征函數來近似表達原隨機場，其對應的均方誤差為：

(3)

通常需要對特征值從大到小進行排序，取前m項從而使得在允許誤差下K-L級數展開項數最少。

對于高斯隨機場，ξi為獨立標準高斯隨機變量；而對于非高斯隨機場，ξi通常難以獲得。因此K-L級數展開主要用于表達高斯隨機場，而對于非高斯隨機場則采用場轉換理論將K-L級數展開表達的高斯場轉換為對應的非高斯場。由于土性參數如抗剪強度通常是非負的，理論上不適合采用高斯場模擬，因此實際中常假設為對數正態隨機場，并采用下式進行轉換：

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

(4)

2 稀疏多項式混沌展開方法

2.1 多項式混沌展開

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

(5)

=γαδαβ

(6)

2.2 展開式系數求解

建立PCE代理模型的關鍵在于展開式系數求解。為便于計算，需要對展開式進行截斷，取最大階數p以內的展開項，即|α|≤p對應的基函數。對于非侵入式模型構建，常用的系數求解方法有投影法、配點法、回歸法等。相較而言，在調用相同原計算模型次數的情況下，回歸法的收斂速度更快，因此這里采用回歸法進行求解。首先通過實驗設計得到輸入空間內的一組樣本X=[x(1),x(2),…,x(N)]T，然后調用原計算模型得到對應的響應量Y=[y(1),y(2),…,y(N)]T，則展開式系數為：

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

(7)

若對于n維輸入變量采用p階PCE擬合其響應關系，則展開項總項數P為：

(8)

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

2.3 正交多項式基函數的稀疏化選取

根據稀疏效應準則，實際問題中占主導的展開式通常是僅包含單獨輸入變量的多項式基函數和變量之間的低階交叉項，而其他多項式對應的系數趨近于零，對結果影響較小，可以忽略。因此，可以通過合理選取多項式基函數提高PCE模型建立的效率。

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

(9)

(10)

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

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

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

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

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

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

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

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

圖1 分析流程

(2)對土性隨機場進行K-L分解，求得截斷誤差內的K-L表達。

(3)采用拉丁超立方抽樣獲取實驗設計樣本。

(4)建立邊坡穩定性計算模型，計算樣本對應的邊坡安全系數。

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

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

(11)

4 算例分析

4.1 一維隨機場

圖2 邊坡示意/m

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

設K-L截斷允許誤差為5%，則需要前21階特征值和特征函數。允許的LOO誤差設為1×10-8，初始實驗設計樣本大小和每次擴充樣本大小均為20。

圖3給出了sPCE代理模型構建的LOO誤差演化，經過11次迭代降低到允許誤差范圍內。圖4給出了該誤差下構建的sPCE模型在訓練集上的結果對比。

圖3 LOO誤差變化

圖4 sPCE模型在訓練集上的結果對比

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

圖5 sPCE模型在測試集上的結果對比

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

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

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

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

圖7 不排水抗剪強度參數空間變異特征模式的敏感性

圖8 前4階特征值和特征函數

4.2 二維隨機場

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

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

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

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

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

圖11 不排水抗剪強度參數空間變異特征模式的敏感性

圖12 前2階特征值和特征函數

5 結 論

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

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

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