土質邊坡位移概率反分析與失穩概率預測

姜 廣 倫，孔 存 芝，南 驍 聰，王 升

(1.山東高速云南發展有限公司，云南 昆明 650206； 2.山東高速四川產業發展有限公司，四川 成都 610041；3.山東高速工程檢測有限公司，山東 濟南 250002； 4.山東省高速公路技術和安全評估重點實驗室，山東 濟南 250002； 5.成都理工大學 地質災害防治與地質環境保護國家重點實驗室，四川 成都 610059)

0 引 言

受地質條件及地形地貌等因素的影響，通過原位試驗獲取巖土體物理力學參數的成本高昂，甚至難以實現。因此，反分析方法已成為獲取巖土體力學參數的重要手段之一。

巖土體自身的非連續性、非均勻性和各向異性使其表現出明顯的不確定性[1]，因此巖土體參數的不確定性是參數反演中十分重要的部分。概率反分析能通過概率分布來科學地量化參數的不確定性。相比確定性反分析[2]，概率反分析需要進行大量的數據統計，實施起來更有難度。貝葉斯概率反分析方法在邊坡風險管理、公路設計、地震模擬等方面得到了不同程度的應用，取得了大量有價值的成果[3-4]。Zhang等[5]基于貝葉斯理論進行了邊坡穩定性的概率反分析研究，提出了更新參數的2種方法：基于回歸模型的優化方法和邊坡本構模型線性簡化模型方法。Zhang等[6]采用馬爾科夫鏈蒙特卡羅(MCMC)方法對多層土體參數進行反分析，有效解決了邊坡多目標參數反演的問題，他還針對MCMC方法效率低下的問題，采用二階響應面法建立了輸入土體參數與輸出穩定性系數的顯式函數關系。左自波等[7]和鄭亞飛等[8]基于監測的時變降雨量數據，通過貝葉斯理論對非飽和土的滲流參數進行了隨機反分析，采用自適應差分演化Metropolis算法在一定程度上解決了通過MCMC得到后驗樣本計算效率低下的問題。然而，目前常用的基于貝葉斯的概率反分析方法一般需要通過MCMC得到數量龐大的隨機樣本。這些隨機樣本都需代入邊坡數值模型計算變形或穩定性系數，這對于實際工程中的復雜邊坡而言，其計算量十分巨大。

同時，在獲得巖土體物理力學參數概率反演結果的基礎上，現有研究大多關注邊坡的持續變形行為，而往往忽略了邊坡的穩定性特征，未能完整構建從邊坡巖土體物理力學參數概率反分析到失穩概率預測的一個完整模型。邊坡失穩概率預測一般可結合極限平衡法[9]或有限元強度折減法[10]。相比極限平衡法，強度折減法可同時考慮受力平衡和變形協調條件，有諸多優勢。Ma等[11]認為強度折減法不僅可以計算失穩概率，還可以結合可靠度方法識別多條代表性滑面，并且計算結果與極限平衡法較為吻合。然而強度折減法的主要缺陷是在結合可靠度方法時計算量過大。響應面法是一種提高計算效率的主流手段，Li等[12]系統地研究了響應面法在邊坡可靠度分析中的應用。此外，近年來通過提高代理模型精度和收斂速度的主動學習代理模型策略在可靠度分析領域得到了快速發展[13-14]。Liu等[15]在邊坡可靠度分析中引入主動學習克里金代理模型，基于圓弧形滑面的極限平衡法提高了可靠度分析的效率。張天龍等[16]采用主動學習徑向基函數代理模型對邊坡進行了系統失穩概率分析。近年來，人工神經網絡[17-18]、支持向量機[19-20]、多變量輸出支持向量機[21-22]等機器學習的方法也在巖土工程領域獲得了廣泛應用，為高效開展邊坡可靠度分析提供了基礎。

本文基于有限差分數值模擬方法，分別構建邊坡位移概率反分析模型和失穩概率預測模型。針對兩種模型計算量大的問題，引入克里金代理模型方法，建立輸出參數(即邊坡穩定性系數和變形等)與輸入參數(如彈性模量、黏聚力和內摩擦角等)間的顯式函數關系式代替邊坡數值分析模型。并以滑坡為例，利用坡表監測位移值作為觀測信息，結合貝葉斯反分析理論，利用克里金代理模型和馬爾科夫鏈蒙特卡羅方法獲得隨機變量的后驗分布信息，并結合多點位移觀測值對反分析結果進行驗證。最后，基于參數反分析結果，結合強度折減法和克里金代理模型，預測邊坡失穩概率，從而構建一套從位移概率反分析到失穩概率預測的邊坡穩定性評價方法，并驗證該方法的適用性與有效性。

1 基于監測位移的參數概率反分析方法

1.1 數值模型與本構模型

本文采用理想彈塑性本構模型，即假定屈服條件與破壞條件相同，選用摩爾-庫倫破壞準則，在彈性階段，彈性模量E、泊松比v滿足：

(1)

(2)

式中：λ，μ為拉梅常數。

在塑性階段，當剪應力與正應力滿足以下條件時發生屈服：

τ=c+σtanφ

(3)

式中：c是黏聚力，φ是內摩擦角。

采用有限差分數值模擬軟件FLAC3D(6.0)[23]進行數值計算，計算收斂條件為不平衡力率(即體系最大不平衡力與典型內力的比率)小于10-5。

1.2 貝葉斯反分析基本理論

通過數值模擬可以動態地展現巖土體的應力應變行為。但由于巖土體自身的非連續性和非均勻性使其表現出高度不確定性，導致巖土體物理力學參數也存在較高的不確定性，因此必須通過概率理論量化參數的不確定性。此外，數值模型也包含不同類型的不確定性，均需進行定量分析。綜上，采用貝葉斯反分析校準土體參數，用先驗分布量化參數的不確定性，同時使用模型誤差因子來表征模型的不確定性[24]，用以下公式來校準土體參數：

y=ε·D(θ)

(4)

式中：y是現場觀測的位移值，ε是模型誤差，θ是隨機變量向量，D(θ)是數值模型的預測值。

研究表明[18，25]，彈性模量E、黏聚力c和內摩擦角φ為影響土體變形最顯著的物理力學參數。因此本文將對這3個參數進行概率反分析。

在給定θ值時，觀測結果為y的可能性被稱為似然函數，即L(θ|y)，它可以表示為θ的條件概率密度函數：

(5)

fε(·)為模型誤差的概率密度函數，在反分析之前把關于θ的信息進行擬合得到概率密度函數，該函數稱為先驗分布，即f(θ)。根據貝葉斯定理，進行反分析更新后得到的θ的后驗概率密度函數為

(6)

式中：k為概率密度函數的歸一化系數。公式(6)是一個貝葉斯問題，可以通過MCMC方法來求解[6]。

1.3 馬爾科夫鏈蒙特卡羅(MCMC)方法

MCMC方法的基本思想是先為目標后驗分布建立一個平穩分布，即馬爾科夫鏈，再通過馬爾科夫鏈產生后驗分布的隨機樣本，從而用蒙特卡羅法求后驗分布對應的期望值、標準差以及極大后驗概率密度值等。MCMC方法允許從任意的分布中抽取樣本，然后校準這些樣本以更好地近似并最終收斂到目標后驗分布。常用的MCMC方法有Gibbs抽樣法、Metropolis抽樣法等[26]。本文采用Metropolis抽樣法進行10萬次MCMC模擬，對邊坡進行位移概率反分析。

1.4 克里金(Kriging)代理模型

概率反分析涉及到數以萬計的重復數值模擬。對于復雜且計算要求高的動態數值模型，例如本次研究中使用的有限差分數值模擬模型，直接進行10萬次模擬需耗時1 042 d，顯然是不可行的。因此，本文利用多響應克里金代理模型近似表達隨機輸入變量與不同位置位移輸出之間的關系。克里金模型是一種基于統計假設的插值方法，已廣泛用于逼近數值計算模型[27]。在克里金模型中，未采樣點θ處的響應被預測為回歸模型和高斯模型的組合：

G(θ)=F(θ)+z(θ)

(7)

式中：F(θ)是通過回歸分析得到的趨勢函數；z(θ)是一個平穩的高斯過程，均值為0，協方差使用高斯相關計算函數[28]。克里金代理模型的詳細構建方法可參考文獻[29]。

本文利用Lophaven等[30]開發的DACE工具包來構建克里金代理模型。代理模型的目標是用最少的訓練樣本來構建最準確的代理，而代理模型的準確性不僅與代理模型的類型有關，還與選擇采樣點的方法有關。因此需要一種有效的樣本設計策略，目前雖然普遍使用隨機抽樣法，但它不能保證目標區域內樣本的均勻分布。為了確保隨機樣本充分覆蓋整個目標區域，通常需要額外的采樣點，這就意味著計算成本的增加。本文采用一種混合自適應采樣策略[31]，可保證目標區域中的采樣點盡可能均勻分布，用于選擇采樣點以構建克里金代理模型。

2 邊坡失穩概率預測方法

2.1 邊坡穩定性功能函數

為了計算邊坡的失穩概率，首先需要構建邊坡穩定性功能函數：

g(θ)=FS(θ)-1

(8)

式中：FS(·)為基于FLAC3D使用強度折減法求解得到的邊坡穩定性系數。當功能函數g(·)大于0，則認為邊坡處于穩定狀態；反之，則不穩定。因此，邊坡的失穩概率可表達成如下的積分形式：

(9)

2.2 強度折減法及主動學習代理模型

強度折減法是一種通過不斷按比例折減土體強度參數來求解邊坡穩定性系數的方法。相比極限平衡法，強度折減法不需要事先假定滑面形狀與位置，同時可以對邊坡失穩過程進行動力分析，計算應力、應變等指標[32]。

由于強度折減法計算比較耗時，采用FLAC3D計算一次簡單邊坡往往都需要幾分鐘，若利用概率反分析得到的10萬個后驗樣本來進行邊坡失穩概率計算，這就意味著需要模擬10萬次，其計算時間是不可接受的。因此，本文通過高效的主動學習克里金代理模型[33]來構建邊坡隨機變量(即黏聚力和內摩擦角)與功能函數的顯式代數表達式。影響代理模型精度和收斂速度的一個關鍵問題就是如何確定訓練樣本集中訓練點的個數和位置。本文通過主動學習函數[33]，選出新的最優樣本點加入到樣本集中，并更新代理模型，如此往復循環迭代。選取最優樣本點的原則為：① 位于極限狀態面附近區域，以提供最有效的信息來提高模型的計算精度；② 遠離訓練樣本集中已有點，避免冗余采樣，以便加快收斂速度。為此可以設置一個最小距離限值d，可以通過公式(10)計算[14]：

(10)

式中：λ為一常數項，依據經驗在0.1～0.5之間取值，本文取0.2[16]。常數項之后是一個距離值，表示當前樣本集中最稀疏的采樣區域處兩個樣本點ui、uj之間的距離。該距離值可以有效地避免局部采樣過密的情況，以便更快地遍布極限狀態面附近區域，使得訓練過程更加高效穩健。

在代理模型達到收斂條件后，便可用代理模型來代替實際功能函數。根據最后建立的克里金代理模型，便可以結合蒙特卡羅模擬快速計算出滑坡失穩概率。

3 方法實施流程

土質邊坡位移概率反分析與失穩概率預測的整體方法流程如圖1所示。方法實現步驟如下：

(1) 構建邊坡有限差分強度折減法數值模型，并確定隨機變量先驗信息和其他基本計算參數。

(2) 利用Kriging算法構建數值模型的代理模型，建立輸入隨機變量與輸出位移的顯示表達式。

(3) 利用實際監測的特征點位移值構建似然函數。

(4) 通過貝葉斯理論更新隨機變量的先驗信息，并結合MCMC方法產生隨機變量的后驗分布。

(5) 在第4步得到的后驗樣本中選取少量樣本點構成初始訓練樣本帶入數值模型，計算獲得穩定性功能函數值。

(6) 通過訓練樣本訓練穩定性功能函數的Kriging代理模型。

(7) 利用第6步中獲得的Kriging代理模型預測后驗樣本中所有點的功能函數值，同時計算當前迭代步驟下的失穩概率。

(8) 判斷Kriging代理模型訓練是否收斂。若收斂，則把最后一次迭代過程中得到的失穩概率作為邊坡的失穩概率預測結果輸出；若不收斂，則利用主動學習函數在后驗樣本中選出最優點添加到訓練樣本中，并返回至第6步重新訓練Kriging代理模型，直至模型訓練收斂。

4 工程算例

4.1 算例概況

為了驗證本文所提方法的有效性，采用吳敏之等[34]介紹的邊坡案例進行驗證，如圖2所示。該邊坡可分為上下2層，上層為碎石土，下層為微風化花崗巖。共有2個位移監測點，經現場監測數據顯示，邊坡發生失穩時A監測點X方向位移大小為23 mm，B監測點X方向位移大小為18 mm[34]。本文將用A監測點的位移監測數據作為觀測信息進行貝葉斯概率反分析，預測B監測點的位移，并進行對比，從而驗證反分析的有效性。

4.2 參數先驗信息及模型誤差

用于FLAC3D計算的數值模型參數如表1所列。根據該邊坡土體的試驗結果[34]，碎石土彈性模量為30 MPa，黏聚力為29.2 kPa，內摩擦角為15.5°，故將試驗結果的值作為待反演參數的先驗均值，同時按照變異系數0.20[35-36]設計標準差，各參數的先驗信息如表2所列。

表1 數值模型計算參數Tab.1 Input parameters of numerical model

表2 待反演參數的先驗信息Tab.2 Priori information of the parameters

根據Kung等[37]的研究表明，模型誤差可以近似認為服從ε～N(1，0.252)的隨機變量。

4.3 數值模型與代理模型

對該邊坡的幾何模型全部采用結構化進行劃分，該模型由33 258個節點數和21 600個單元數構成，如圖3所示。

把彈性模量E、黏聚力c和內摩擦角φ作為輸入，以A監測點的X方向位移為輸出，通過自適應采樣算法克里金模型構建FLAC3D數值模型的代理模型，如圖4所示。采樣次數為500，利用工作站(Intel(R) Core(TM)i9-10900k CPU @ 3.7 GHz 32 GB RAM)進行計算，共計12.5 h。代理模型的確定性系數R2為0.997，足以滿足計算要求。

4.4 概率反分析結果

在利用MCMC求解貝葉斯方程時，馬爾科夫鏈總樣本數量為11萬個，剔除最開始的波動段數目1萬個，用于產生后驗的樣本數目為10萬個。馬爾科夫鏈樣本的接受率隨著縮減因子的增加而降低，當接受率為20%～40%時，馬爾科夫鏈是最有效的[38]。因此本文調整縮減因子為1.2，使得對應的接受率為29%，在有效范圍內。由于使用了代理模型，每一次的計算時間縮短至0.003 s，因此通過MCMC采樣得到的10萬次后驗分布(見圖5～7)，反分析總計算時間約為5 min。

先驗分布與后驗分布的統計特征對比如下：① 彈性模量的均值經反分析后變大，說明先驗分布低估了彈性模量的值。黏聚力與內摩擦角的均值經反分析后變小，說明先驗分布高估了碎石土的抗剪強度。② 彈性模量的標準差由6 MPa降低至4.69 MPa，黏聚力的標準差由5.84 kPa降低至1.22 kPa，內摩擦角的標準差由3.10°降低至1.82°。從后驗分布中可以看出，后驗分布較先驗分布明顯收窄，各參數的標準差都明顯降低，說明貝葉斯反分析大大降低了參數的不確定性。

彈性模量與內摩擦角的后驗分布接近廣義極大值分布，黏聚力的后驗分布接近正態分布。彈性模量后驗分布95%的置信區間為[28.76 MPa，46.78 MPa]；黏聚力后驗分布95%的置信區間為[21.7 kPa，25.32 kPa]；內摩擦角后驗分布95%的置信區間為[9.82°，16.10°]。此結果可為邊坡失穩概率預測及邊坡穩定性定量風險評估提供數據支撐。

4.5 邊坡位移預測及驗證

為了驗證反分析結果的有效性，通常需要將結果進行反饋分析，并與實際觀測值進行比較。本文用3個參數的后驗概率密度極大值(彈性模量34.7 MPa，黏聚力22.97 kPa，內摩擦角10.8°)對該邊坡進行位移預測，X方向的位移云圖如圖8所示。預測出X方向最大位移為23.6 mm，發生在坡腳處。B監測點X方向的位移大小為20.2 mm，與觀測值18.0 mm相差2.2 mm。這種偏差主要是后驗參數的不確定性造成的。土體參數的后驗概率密度極大值代表土體參數最可能的取值，但這并不代表土體參數的真實值。因為從后驗分布中可以看出(見圖5～7)，土體參數的組合是非唯一的，使用參數的后驗概率密度極大值預測的位移跟實際監測位移相比存在較小的偏差是正常的。因此預測結果是合理的，由此證明了此方法的有效性。

5 邊坡穩定性分析及失穩概率預測

利用彈性模量、黏聚力、內摩擦角的后驗概率密度極大值以及其他參數進行FLAC3D強度折減計算。計算獲得邊坡的穩定系數為0.80，說明邊坡處于不穩定狀態，極易發生破壞失穩，這與實際觀察結果一致。最終形成貫通的最大剪應變增量帶，即潛在滑帶，如圖9所示。

通過主動學習克里金代理模型策略建立的邊坡失穩概率分析模型如圖10所示，訓練樣本總數為100個。收斂條件為：最后5次迭代過程所計算的失穩概率指標變異系數不大于0.001[16]。根據計算，得出該邊坡的失穩概率為86.77%，說明發生失穩的可能性極高。

本文失穩概率預測總計算時間為25 h。由于FLAC3D進行強度折減計算一次需要約15 min，若直接對10萬個后驗樣本進行數值計算則需要1 042 d，這顯然是不現實的。由此可以看出利用主動學習克里金代理模型策略構建邊坡失穩概率預測模型是極其高效的。

6 結 論

本文基于邊坡坡表位移監測數據，通過貝葉斯理論對土體物理力學參數進行概率反分析，并將反分析結果進行反饋分析，驗證了該方法的可靠性。同時利用反分析得到的隨機變量后驗分布對滑坡失穩概率做出了預測，主要得出以下結論。

(1) 貝葉斯概率反分析方法可充分利用邊坡坡表變形監測數據作為觀察信息，結合先驗信息推斷土體物理力學參數后驗信息。通過參數更新，可有效校準土體參數的平均值并極大地降低其不確定性。

(2) 利用反分析獲得的參數極大概率密度值對坡表B點位移進行預測，預測結果與實際監測結果吻合較好，驗證了貝葉斯反分析方法的適用性。

(3) 針對邊坡概率反分析計算量大的問題，引入自適應采樣多響應克里金代理模型，將計算時間降低到12.5 h；而針對邊坡失穩概率預測計算量大的問題，引入主動學習克里金代理模型，將計算時間降低到25 h，極大地提高了邊坡概率分析的計算效率。