基于支持向量機和貝葉斯方法的巖體參數反分析

張 雅 賢，侯 中 杰

(1.中央軍委后勤保障部安置住房保障中心，北京 100036； 2.重慶大學 土木工程學院，重慶 400045)

0 引 言

隨著對巖土體認識的不斷深入，人們就巖土體本身的性質具有不確定性這點達成了共識。在巖土體參數不確定性表征方面，國內外學者做了大量研究工作。主要的表征方式為基于統計特征的概率分布形式，包括正態分布[1-2]、對數正態分布[3-4]、Beta分布[5]、極值I型和Weibull分布[6]。早期學者們主要是采用隨機變量模型來對巖土參數的不確定性加以研究。然而，類似的研究在刻畫巖土材料的隨機特性時都是將同層巖土體看作均質材料，即僅在抽樣時對參數取不同數值，而針對同一樣本在不同空間處的參數卻往往取的是同一值，這樣可能導致與實際情況不相符合。近些年來隨著土體空間變異性概念的提出[7]，眾多學者開始采用參數隨機場模型來研究巖體參數并取得眾多研究成果，這意味著對巖體參數的研究實現了從“點特性”到“空間特性”的轉變。Degroot等[8]為了更好地研究土體參數的自相關性，引入了極大似然估計法以及相關函數法。Phoon等[9-10]則采用統計學理論分析了不同實驗參數轉換的相關性、不確定性以及不同巖土材料的變異性等問題。梅松華等[11]結合工程實例，采用不同的方法實現了對土體參數相關距離的求解。隨著隨機場模型的進一步應用，眾多以空間隨機場作為理論基礎的巖土工程可靠性分析方法逐漸被廣泛應用[12-16]。

當前，學術界在巖土工程參數不確定性方面開展的研究非常多。但是截至目前，學術界和工業界仍尚未形成基于整個施工過程的巖土參數不確定性分析方法，而與巖土參數不確定性表征方式相關的研究則較為常見。為此，本文在巖體參數識別的基礎上，融合貝葉斯理論與多輸出支持向量機算法，提出基于位移與松弛深度的巖體參數概率反分析方法，并以白鶴灘水電站左岸壩基邊坡為例，研究了邊坡巖體參數的不確定性特征及其隨施工開挖的動態更新過程。

1 多輸出支持向量機和貝葉斯方法原理

1.1 多輸出支持向量機原理

多輸出最小二乘支持向量機(MSVM)主要源自早期的單輸出支持向量機(SVM)[17-20]。對樣本集為{(x1，y1)，(x2，y2)，…，(xl，yl)}，x∈Rd，y∈Rm，其中，l、d和m分別表示樣本數、輸入變量的維數和輸出變量的維數，多輸出最小二乘支持向量機(MSVM)對應的最小化問題為

(1)

式中：wi為權向量，b為閾值向量，e為誤差，φ1(xi)為核函數，γ為誤差因子。引入Lagrange乘子α，將式(1)轉化為

(2)

以w、b、e和α為自變量對式(2)求偏導，并令偏導等于零，易得：

(3)

i=1，2，3，…，l；j=1，2，3，…，m

從而能形成關于w、b、e和α的線性方程，對該方程求解即可得出多輸出最小二乘支持向量回歸函數表達式：

(4)

式中：K(xi，x)=φ(xi)·φ(x)。

1.2 貝葉斯方法原理

貝葉斯方法(Bayesian method)基于概率理論，在已知事情結果的情況下反推事情發生的原因。假設隨機事件A發生的概率為P(A)，隨機事件B發生的概率為P(B)，那么A在B發生的情況下發生的概率P(A|B)為

(5)

這就是貝葉斯公式的最核心的一般表達式，將上式運用到反分析中時，事件A則是先驗信息，B是觀測結果，運用觀測結果來反推先驗信息就是反分析。式(1)所示的貝葉斯公式是參數反分析的基礎。

(6)

式中：p(b|a)稱作似然函數。此時，根據式(5)所示的貝葉斯理論得到樣本與參數的聯合分布為

d(b，a)=p(b|a)Rp(a)

(7)

在有觀測信息b的情況下，再對參數a進行新的估計。考慮到：

d(b，a)=R(a|b)z(b)

(8)

式中：z(b)是b獨立的邊緣概率密度函數，求解方法為

(9)

由式(8)得到參數a的估計：

(10)

根據統計資料，有先驗分布Rp(a)，觀測資料為b，則代入式(6)就能推斷出參數a新的概率分布。這種新的分布融入了實測資料，將更加準確，稱之為后驗分布。

由以上可知，為推斷出可靠的后驗參數估計，先驗信息尤為關鍵。有限的先驗信息一般通過實地鉆孔測量得到，或者通過大量的文獻調查類比得到。得到先驗信息后，可根據不同的參數選擇不同構建先驗分布的方式，常用高斯曲線來構建概率密度函數，利用均值和變異系數來構造高斯函數作為先驗分布。

2 位移-松動圈與巖體參數間的支持向量機模型

2.1 模型建立

參數的不確定性是直接影響工程設計及安全性評價可靠性的重要因素。采用多輸出支持向量機(MSVM)建立監測位移與巖體參數之間的關系，然后采用貝葉斯方法融合已有的信息對巖土力學參數不確定性分布規律做出合理的估計。MSVM可用于表征監測位移、松動圈與彈性模量等巖土力學參數之間存在的非線性關系，MSVM(X)模型定義如下：

(11)

式中：X=(x1，x2，…，xN)，xi是巖土力學參數(如彈性模量、摩擦角等)組成的向量，N為力學參數的個數。Y=(y1，y2，…，yQ)是監測數據的Q維向量，包括位移和松弛深度。

將貝葉斯理論融合在MSVM模型中，不同測點力學參數與位移松弛深度之間的關系可以通過MSVM模型得出。監測點的預測結果(y1，y2，…，yk)與相應的監測結果(可以用Ymon1，Ymon2，…，Ymonk來表示)相等的可能性是參數θ的條件概率密度函數，即：

(12)

模型參數的先驗分布f(θ)由參數的均值向量μθ=(μθ1，μθ2，…，μθn)和協方差矩陣Cθ給定[21]，其中：

(13)

(14)

式(12) 中的Nk(y)定義如下：

(15)

(16)

后驗分布可以通過優化方法獲取，在求得后驗均值之后，可用下式求得后驗標準差：

(17)

其中，

(18)

(19)

2.2 模型實現基本步驟

B-MSVM采用貝葉斯算法更新巖土力學參數，并且考慮了巖土體內在的不確定性。位移-松弛深度與巖土力學參數的對應關系對更新算法至關重要，本文采用 MSVM來取代傳統的數值分析或模型計算方法，其主要步驟如圖1所示。

圖1 基于支持向量機和貝葉斯方法的反 分析基本步驟Fig.1 Basic steps of back analysis based on support vector machine and Bayesian method

(1) 收集監測面分布、地質條件等實際工程信息，并借助ABAQUS、Phase2等主流有限元軟件建立數值模型。

(2) 根據步驟(1)中所收集到的工程信息，對用于識別的訓練組參數分布類型及其范圍加以確定，并通過抽樣法構建目標樣本。

(3) 借助數值模型計算確定樣本位移。

(4) 構建MSVM樣本。構建巖土力學參數組合，樣本集由采樣點及其數據構成。通過有限元法或其他方法確定每個樣本集的松弛深度、位移并定義MSVM模型。通過對公式(2)求解進而獲得MSVM 模型。

(5) 由公式(12)建立貝葉斯更新模型，并計算識別參數。

(6) 由公式(17)確定參數的后驗標準差，并進一步得到巖土力學參數。

(7) 根據新的開挖步中監測數據，實時更新土體參數和不確定性分布情況。

3 工程實例應用

3.1 工程簡介

白鶴灘水電站左岸壩基邊坡如圖2所示。邊坡主要由砂巖、泥巖、灰巖、凝灰巖組成。裸露的巖石主要有6類：泥巖石灰巖互層、泥巖石灰巖混合巖、層狀凝灰巖、砂巖、砂巖泥巖互層、泥巖。左壩肩屬陡峭的反傾巖質邊坡，開挖高度從635 m降低到215 m，高差達420 m，開挖切角的范圍35°～65°，分布約14個15～20 m高的長椅狀臺階。邊坡上存在50多個大大小小的斷層，較大的斷層對邊坡的穩定性有很大影響。

圖2 白鶴灘水電站左岸壩基邊坡實景圖Fig.2 Realistic view of left bank dam foundation slope of Baihetan Hydropower Station

為記錄開挖期間及開挖完成之后邊坡的變形，共選取了13條監測剖面，位置分布如圖3所示。監測表面位移和深部變形的儀器均在邊坡上沿著監測剖面布置，這些儀器可以監測邊坡的變形情況，現場收集了通過多點伸長計記錄的變形數據。

圖3 監測斷面位置示意Fig.3 Location of monitoring section

柱狀節理玄武巖在邊坡測試區域內集中發育，且柱狀節理被較多的緩傾角結構面和錯動帶切割，除一些受構造影響的區域，柱狀節理玄武巖的柱體軸線傾角基本均為60°～80°。

3.2 反分析過程

選取典型的監測斷面1-1來反演邊坡巖土體的力學參數，沿著1-1地質剖面的巖床包括砂巖、砂巖泥巖互層、泥巖等。此外，考慮了12個重要的斷層建立有限元數值模型。數值計算中破壞準則采用摩爾庫倫準則，斷層采用節理單元模擬。據原位應力測試的結果，采用式(20)、(21)模擬計算原位應力。對于全風化和強風化的巖體：

(20)

式中：σx和σy分別表示與深度方向垂直的x方向和y方向的原位應力，σz表示深度方向上的應力，γ表示巖石密度，μ表示泊松比，h表示巖石上覆層的厚度。

對弱風化和無擾動的巖體：

(21)

式中：Kx和Ky分別表示原位應力在x水平方向和y水平方向的側壓力系數。

通過數值模型、現場試驗或者室內試驗可以計算得出復合巖體的物理力學參數。根據給定的邊坡穩定性要求，反分析計算出符合正態分布的9個巖土力學參數：主要包含7個彈性模量，其中1個是弱風化巖體的彈性模量，另6個是6種巖體的彈性模量；另外2個參數是弱風化巖體和無擾動巖體的原位應力測試的側壓力系數。考慮這9個參數的原因一方面是因為這些參數對邊坡穩定性有較大影響；另一方面是因為如果選擇更多的參數會導致計算過程更為復雜，工程應用較為不便。

利用拉丁超立方體抽樣方法建立了一個總數為50的樣本組。基于前文的數值模型，采用有限元軟件計算每個樣本參數對應的位移值，生成50組含有不同巖土力學參數和位移的訓練樣本。借助MSVM算法來進一步刻畫巖土力學參數與位移之間的聯系。采用同樣方法構建20個MSVM模型測試樣本。通過研究計算位移與MSVM預測位移之間的聯系得出兩者的對比如圖4所示。分析圖4可以得知兩者具有較高一致性，該對比結果也說明在概率反演分析中可用MSVM方法替代數值模型。

圖4 有限元計算位移和MSVM預測位移的對比Fig.4 Comparison of displacement calculated by FEM and predicted by MSVM

3.3 結果比較及驗證

基于本文提出的反分析方法，可以得到14個不同開挖步條件下7種不同巖體的彈性模量和2個原位側壓力系數。為了探討不同開挖條件下B-MSVM方法的適用性，又進行了兩階段分析；第一階段是采用監測點的位移反分析計算了13個開挖步后的巖土力學參數；第二階段是基于第一階段獲取的巖體參數的平均值和標準差，采用第14個開挖步后的位移反分析計算巖土力學參數。

不同開挖步下由計算結果反演的巖土力學參數分布如圖5所示。 隨著開挖的進行，每個開挖步下參數的分布變窄，標準差也隨之降低。這種參數分布的變化表明，所提出的方法可以更加準確地進行反分析，同時所提出的模型可以用于指導工程設計和施工。

通過監測點的附加信息和數據生成的標準差加深了對巖土力學參數不確定性的理解。與傳統的反分析方法相比，概率反分析提供了更多的關于巖土力學參數的信息，并用一種切實可行的方式體現了巖土力學參數的不確定性。這種方法比確定性方法更合理，并能更進一步地匹配地質工程的復雜性和不確定性特征。

采用有限元數值軟件基于已識別的參數計算了監測點的位移，結果如圖6所示。基于參數識別計算的位移和實測位移具有較好的一致性。基于參數識別的第二階段計算位移比第一階段計算位移更接近實測位移，這表明位移的持續更新能提高反分析模型的準確性，概率反分析可以用來進行邊坡的可靠度分析和可靠度設計。

MSVM方法預測的監測點的位移分布如圖7所示。圖中位移分布變窄，巖土體參數分布也變窄，這表明附加的開挖信息降低了預測位移的不確定性。這不僅說明了位移的不確定性，也說明更多的監測信息(先驗信息)可以有效地降低不確定性的影響。

圖5 不同開挖步下不同參數的不確定性對比Fig.5 Comparison of uncertainty of different parameters under different excavation steps

圖6 監測點的計算位移與預測位移對比Fig.6 Comparison of calculated displacement and predicted displacement of monitoring points

圖7 不同高程開挖時B-MSVM法預測的監測點位移 分布對比Fig.7 Comparison of displacement distribution of monitoring points predicted by Bayesian-MSVM method during excavation at different elevations

4 結 論

將多輸出支持向量機方法和貝葉斯理論引入到了邊坡巖體物理力學參數的不確定性反分析過程之中，從而實現參數的不確定性表征，并且實現了不確定性特征隨著施工過程的動態更新。提出基于位移與松弛深度的巖體參數概率反分析方法，將該方法應用于白鶴灘水電站左岸壩基邊坡，主要得出如下結論：

(1) 概率反分析方法提供了更多的關于巖土力學參數的信息，比確定性方法更合理，并能更進一步地匹配地質工程的復雜性和不確定性特征。

(2) 更多的監測信息(先驗信息)可以有效降低不確定性的影響，位移的持續更新能提高概率反分析模型的準確性，相應的反分析的結果可以用來進行邊坡的可靠度分析和可靠度設計。

(3) 本文提出的方法可以實現參數及其不確定性特征隨施工過程的動態更新，可以推廣應用到類似工程中。