RBF神經(jīng)網(wǎng)絡可靠度分析方法在邊坡穩(wěn)定性研究中的應用*

何永波，李 青，張 寧，李闖將

(中國計量大學 災害監(jiān)測技術(shù)與儀器國家地方聯(lián)合工程實驗室，浙江 杭州 310018)

0 引言

邊坡穩(wěn)定性分析是對滑坡災害預警和治理的關鍵方法之一，目前邊坡穩(wěn)定性的分析方法主要分為2大類：確定性分析法和不確定分析法。確定性分析法目前應用較廣泛的是以極限平衡理論為基礎的安全系數(shù)法，但由于實際邊坡存在大量的不確定性，包括計算參數(shù)的變異性和隨機性[1]，該方法很難反映出邊坡的實際工作狀態(tài)；不確定分析中，最具代表性的是基于概率分析的可靠度分析方法，該理論把各不確定因素看作隨機變量來分析邊坡破壞的可能性，根據(jù)已知的隨機變量統(tǒng)計參數(shù)和概率分布模型以及給定的功能函數(shù)估計邊坡在規(guī)定時間和條件下完成預定功能的概率[2],將不能完成功能的概率稱之為失效概率[3]。

可靠度理論是建立在土體具有的抗力大于荷載效應的概率基礎上進行設計和校核的，因此分析結(jié)果更符合客觀實際。邊坡可靠度的計算常用方法有一次二階矩法(中心點法、設計驗算點法)、蒙特卡羅模擬法、響應面法等。一次二階矩法計算的精度較低，蒙特卡羅法簡單且計算精度高，但是計算量偏大[4]。近年來，隨著非線性理論的發(fā)展，基于替代模型的可靠度分析開始應用起來，利用替代模型來近似構(gòu)造函數(shù)，將復雜的邊坡穩(wěn)定隱式功能函數(shù)顯示化[5]，降低工作量，提升工作效率。常見的替代模型有利用響應面對邊坡可靠度研究，神經(jīng)網(wǎng)絡、支持向量機[6]等智能算法，給邊坡的穩(wěn)定性分析開辟了新的途徑。針對之前各研究方法所遇到的問題，本文提出一種基于粒子群算法(PSO)優(yōu)化徑向基函數(shù)(RBF)神經(jīng)網(wǎng)絡的邊坡可靠度分析方法。

1 分析方案

本文主要根據(jù)地質(zhì)勘探和室內(nèi)外的土工試驗得到邊坡巖土體物理參數(shù)(黏聚力c值和內(nèi)摩擦角φ值)；利用強度折減法計算c,φ對應的安全系數(shù)；將各組數(shù)據(jù)作為模型的訓練樣本帶入RBF網(wǎng)絡模型進行訓練，利用模型強大的數(shù)據(jù)擬合能力，映射出安全系數(shù)和c，φ之間的關系，并使用PSO算法進一步優(yōu)化，構(gòu)建響應面功能函數(shù)；結(jié)合蒙特卡羅法(Monte Carlo)產(chǎn)生大量的隨機數(shù)樣本模擬求解邊坡失穩(wěn)概率，分析方案見圖1。

圖1 分析方案設計Fig.1 Design of analysis scheme

2 基于RBF網(wǎng)絡的可靠度分析

2.1 ABAQUS強度折減法

強度折減法最早由Zienkiewicz等提出，后被許多學者廣泛采用并提出了抗剪強度折減系數(shù)(SSRF)的概念[7]，定義為：假設外載荷保持不變，邊坡內(nèi)土體所能提供的最大抗剪強度與外載荷在邊坡內(nèi)所產(chǎn)生的實際剪應力之比。當邊坡內(nèi)所有土體抗剪強度的發(fā)揮程度相同時，這種抗剪強度系數(shù)相當于穩(wěn)定安全系數(shù)FS。

折減后的抗剪強度參數(shù)可分別表達為:

cm=c/Fr

(1)

φm=arctan(tanφ/Fr)

(2)

式中：c和φ是土體所能提供的抗剪強度參數(shù)；cm和φm是維持平衡所需的抗剪強度參數(shù)；Fr是強度折減系數(shù)。

使用強度減法所選取的屈服準則[8]是Mohr-Coulomb：

(3)

式中：c和φ是土體黏聚力和內(nèi)摩擦角；I1為應力張量的第一不變量；J2是應力偏張量的第二不變量；θ是應力洛德角。因ABAQUS具有求解非線性，處理非均質(zhì)問題以及模擬各種復雜的材料的優(yōu)點，故本文選用此方法進行模擬計算。

2.2 RBF神經(jīng)網(wǎng)絡

RBF神經(jīng)網(wǎng)絡是于1989年由Mooden及Darken提出，是由輸入層、隱含層和輸出層構(gòu)成的三層前饋式神經(jīng)網(wǎng)絡(見圖2)，在參數(shù)選擇適當?shù)那疤嵯拢淠軌蛞越o定的精度逼近任意的連續(xù)函數(shù)[9]。圖中R為網(wǎng)絡輸入數(shù)；P為輸入向量；S為神經(jīng)元數(shù)；a和IW為該層輸出和權(quán)值矩陣；radbas為徑向基函數(shù)；b為隱含層閾值。隱含層常采用K-means算法進行訓練，輸出層采用遞歸最小二乘法進行訓練。它用徑向基函數(shù)作為激勵函數(shù)，由于最常用的基函數(shù)為高斯函數(shù)，定義第i個隱含單元的激活函數(shù)為:

(4)

隱含層到輸出層映射為：

(5)

其中，可通過最小二乘法獲得權(quán)重系數(shù)ω，向量可表示為ω=[ω1,ω2,…,ωH]，一般通過最小化誤差指標函數(shù)訓練網(wǎng)絡,選取均方誤差函數(shù)為:

(6)

式中：N代表訓練樣本點的個數(shù)；ei代表第i個節(jié)點的輸出誤差,可以表示為:

(7)

圖2 RBF徑向基函數(shù)神經(jīng)網(wǎng)絡結(jié)構(gòu)Fig.2 Structure of RBF neural network

2.3 PSO算法

PSO算法[10]是基于群體行為的搜索算法，群體中每個粒子在迭代的過程中不斷改變本身的速度矢量v和位移矢量x，尋找到全局最佳的位置，粒子的迭代過程滿足：

(8)

(9)

ωt=ω2-t(ω2-ω1)/T

(10)

式中：ω1和ω2分別表示初始和最終迭代權(quán)重；t為當前迭代次數(shù)；T為最大迭代次數(shù)。

PSO優(yōu)化RBF模型[11]步驟主要包括：

1)模型初始化，種群的規(guī)模、迭代的次數(shù)和權(quán)重；

2)計算每個粒子的適應度，如式(11)：

(11)

式中：yi和fi表示實測值和預測值。

3)以種群中適應度最小的粒子作為gbest

初始值，將粒子當前位置作為最優(yōu)pbest，找出具有最優(yōu)適應度值的粒子位置作為pbest。

4)比較個體和種群的最優(yōu)解適應度，更小的作為gbest。

5)更新粒子的速度和位置以及迭代的權(quán)重，直到迭代次數(shù)滿足結(jié)束條件，將gbest對應的粒子作為RBF的參數(shù)。

2.4 Monte Carlo模擬方法

Monte Carlo模擬方法是進行可靠度計算的重要手段，又稱為隨機模擬方法或統(tǒng)計試驗方法，由于其限制較小，思路簡單，得到了較為廣泛的應用[12]。可以假設結(jié)構(gòu)的功能函數(shù)已知以及基本隨機變量的概率分布，當選取樣本數(shù)量足夠大時，事件實際發(fā)生的概率可以通過頻率近似得到。

按照滑坡的巖土體性質(zhì)、變形機制及其受力狀態(tài)，首先確定狀態(tài)變量x1,x2,…,xn的參數(shù)統(tǒng)計值及概率分布，參數(shù)分別代表重度，黏聚力，內(nèi)摩擦角等，根據(jù)RBF擬合的響應面函數(shù)確定邊坡的結(jié)構(gòu)功能函數(shù)為：F=g(x1,x2,…xn)，極限狀態(tài)方程可表示為Z=F-1,如此重復N次，得到N個相對獨立的樣本值Z1,Z2,…,Zn，若定Z<0為滑坡失效事件，則在N次抽樣中出現(xiàn)M次，則由伯努利大數(shù)定理可知，失效概率為：

(12)

式(12)為蒙特卡羅法計算出的失效概率，對于得到的N組Z,其均值和標準差分別為：

(13)

(14)

用β表示可靠度指標，則β可表示為：

β=μz/σz

(15)

則失效概率：

Pf=Φ(-β)=1-Φ(β)

(16)

式中：Φ(β)為標準正態(tài)分布。

3 實例分析

半嶺滑坡位于浙江省麗水市遂昌縣黃沙腰鎮(zhèn)村下半嶺自然村，距離遂昌縣城約85 km，半嶺滑坡所在地區(qū)雨季降雨強度大，山區(qū)水流流速大，極易引起滑坡。滑坡平面上總體呈“舌狀”，主軸傾向215°，滑坡堆積體平面投影長540 m，前緣寬約220 m，后緣寬約190 m,高程460～670 m，根據(jù)鉆孔揭露厚度為0～25.9 m，滑坡區(qū)域的現(xiàn)場勘察如圖3所示。

圖3 半嶺滑坡工程現(xiàn)場Fig.3 Site map of Banling landslide project

根據(jù)現(xiàn)場調(diào)查發(fā)現(xiàn)位于滑坡前緣2處小規(guī)模滑坡現(xiàn)象，2018年6月多雨天氣，鉆孔ZK02附近再次發(fā)現(xiàn)長裂縫，該滑坡地下水豐富，軟弱土層浸泡軟化，工程力學性質(zhì)降低，導致堆積體發(fā)生失衡，沿軟弱夾層發(fā)生蠕動，工程地質(zhì)平面圖如圖4，工程地質(zhì)剖面圖如圖5。

根據(jù)半嶺A區(qū)滑坡目前處于臨界狀態(tài)這一現(xiàn)狀，結(jié)合土工實驗以及現(xiàn)場勘察得出最終的計算參數(shù)，取多組黏聚力c值和內(nèi)摩擦角φ值，并計算其均值、標準差、變異系數(shù)，見表1。

表1 滑帶土物理參數(shù)Table 1 Physical parameters of soil in slip zone

3.1 求解邊坡安全系數(shù)

參考Dawson等分析的一個均質(zhì)土坡[7]建立和實際邊坡近似的算例。將飽和狀態(tài)下參數(shù)的均值作為輸入：黏聚力c=13.38 kPa，內(nèi)摩擦角φ=12.26°，重度γ=20 kN/m3，楊氏彈性模量E=20 MPa，變異系數(shù)δ=0.315，泊松比ν=0.3。

首先，在ABAQUS軟件中建立部件，繪制邊坡幾何輪廓，設置材料屬性以及截面特性；其次，裝配部件以及定義分析步之后需要定義載荷以及邊界條件、劃分網(wǎng)格；最后，修改模型的輸入文件，控制場變量變化，進行結(jié)果分析。

圖4 半嶺工程地質(zhì)平面Fig.4 Engineering geological plane of Banling

圖5 半嶺滑坡地質(zhì)剖面Fig.5 Geological profile of Banling landside

當c和φ值取飽和狀態(tài)均值時，選擇場變量FV1(強度折減系數(shù))和x方向的位移U1作為輸出變量，由圖6可見，以位移拐點作為邊坡穩(wěn)定的評價標準，安全系數(shù)為0.81；若以數(shù)值計算不收斂作為評價標準，對應的FV1為0.83。2個數(shù)值與極限平衡法算出的FS=0.82以及傳遞系數(shù)法計算的FS=0.81相比都比較接近，如圖7所示，說明本方法可行且計算效率高。

圖6 FV1隨U1的變化關系Fig.6 Variation relationship of FV1 with U1

圖7 傳遞系數(shù)法計算安全系數(shù)結(jié)果Fig.7 Calculation results of safety coefficient by transfer coefficient method

經(jīng)過多次計算并結(jié)合《滑坡防治工程勘察規(guī)范》(DZ/T 0218-2016)表2可知，天然狀態(tài)下FS為 1.14 左右，基本穩(wěn)定，與現(xiàn)場勘察時地表裂縫現(xiàn)象是一致的；實驗過程中首先從邊坡的前緣發(fā)生塑性應變，在飽和狀態(tài)下FS為0.82左右，處于不穩(wěn)定狀態(tài)，與勘察期間連續(xù)降雨后坡面出現(xiàn)裂縫是一致的，說明滑坡體內(nèi)地下水豐富且降雨對地下水具有一定的影響，穩(wěn)定系數(shù)下降，若持續(xù)強降雨將可能出現(xiàn)滑坡。

根據(jù)麗水市勘察測繪院監(jiān)測結(jié)果顯示，2017年1月—2017年4月，最大位移量為60 mm,最大變化速率值為10～15 mm/d,邊坡達到較危險值，且在連續(xù)強降雨期間邊坡前緣發(fā)生小范圍的崩塌。該滑坡的主要成因機制是粉質(zhì)黏土層長期受地下水浸泡，軟化，在滑坡堆積體的自重影響下沿該層發(fā)生蠕動。

表2 滑坡穩(wěn)定狀態(tài)劃分等級Table 2 Grade division of landslide stability state

本次實驗選取2個邊坡基本物理參數(shù)c值φ值，將土工實驗得到的天然、飽和數(shù)據(jù)總共48組參數(shù)分別通過ABAQUS軟件逐對進行求解安全系數(shù)，將數(shù)據(jù)按比例7∶3劃分后作為網(wǎng)絡模型訓練和檢驗的樣本，如表3～4。

3.2 訓練結(jié)果分析

經(jīng)過多次的計算，取粒子群的規(guī)模為20，粒子維數(shù)12，最大迭代次數(shù)250次，初始權(quán)重為0.9，c1=c2=1.49，當網(wǎng)絡隱含節(jié)點數(shù)為10左右時，網(wǎng)絡的預測性能最好，決定系數(shù)均在99%以上，仿真的誤差最大不超過0.001；而直接使用RBF網(wǎng)絡預測的決定系數(shù)在90%左右，精度明顯低于PSO-RBF預測的結(jié)果,結(jié)果如圖8所示。由此可見PSO優(yōu)化后的RBF神經(jīng)網(wǎng)絡適合于本文的研究，RBF擬合出的響應面函數(shù)S=f(c,φ)如圖9所示。

表3 PSO-RBF網(wǎng)絡訓練樣本Table 3 Training samples of PSO-RBF network

表4 PSO-RBF網(wǎng)絡檢驗樣本Table 4 Testing samples of PSO-RBF network

圖8 PSO-RBF預測結(jié)果Fig.8 PSO-RBF prediction results

圖9 RBF擬合出功能函數(shù)曲面Fig.9 Performance function surface fitted by RBF

3.3 蒙特卡羅模擬計算

隨機變量通常服從的正分布是非標準正態(tài)分布N(0,σ2)，可采用標準正態(tài)分布N(0,1)的隨機變量x′經(jīng)過線性轉(zhuǎn)換[13]得到：

X=μ+σx′

(17)

用二元函數(shù)變換可得到：

(18)

(19)

式中：X1和X2是標準的正態(tài)分布隨機變量，u1和u2為[0,1]區(qū)間均勻隨機數(shù)。

隨機變量X服從對數(shù)正態(tài)分布[14]可通過公式Y(jié)=lnX轉(zhuǎn)換，則Y服從正態(tài)分布。

根據(jù)文獻[15]可知物理參數(shù)的概率分布特性和變異系數(shù)相關，得到天然和飽和的黏聚力均服從對數(shù)正態(tài)分布；同理，2種工況的內(nèi)摩擦角均服從正態(tài)分布，通過對參數(shù)值的處理后，利用matlab自帶函數(shù)隨機獲取15 000組數(shù)據(jù)所得c和φ值概率密度分布如圖10所示。

圖10 c和φ值概率密度分布Fig.10 Distribution of probability densities of and values

將生成的15 000組數(shù)據(jù)帶入訓練好的模型中，依次可求得相應的安全系數(shù),則根據(jù)公式(12)～(16)得到失穩(wěn)概率和可靠度指標,見表5。

表5 蒙特卡羅法計算結(jié)果Table 5 Calculation results of Monte Carlo

模擬過程中嘗試了多次模擬計算，從小到大依次增加模擬次數(shù)，當模擬次數(shù)2 000以內(nèi)時，失穩(wěn)概率值并不穩(wěn)定，波動比較大，但當模擬的次數(shù)超過6 000次后失穩(wěn)概率值逐漸趨于平穩(wěn)，最終模擬15 000次達到平穩(wěn)狀態(tài)，模擬次數(shù)統(tǒng)計見圖11。

圖11 失穩(wěn)概率隨模擬次數(shù)的變化Fig.11 Variation of instability probability with number of simulation

3.4 不同方法對比分析

通過計算可以得到邊坡天然和飽和2種工況的失穩(wěn)概率和可靠度指標，再利用文獻[5，16]中提到的基于BP網(wǎng)絡的可靠度分析法和直接Monte Carlo法得出2種工況邊坡的失穩(wěn)概率和可靠度，計算結(jié)果見表6。以直接Monte Carlo法計算結(jié)果作為對比對象，基于BP網(wǎng)絡可靠度分析法的誤差為5.77%和9.52%；基于RBF網(wǎng)絡可靠度分析法的誤差為2.33%和7.36%；基于PSO-RBF網(wǎng)絡可靠度分析方法的誤差為0.79%和2.59%。綜上所述，本文的方法計算結(jié)果誤差較小，證明了該方法的可行性。

表6 計算結(jié)果對比Table 6 Comparison of calculation results

4 結(jié)論

1)本文從理論分析結(jié)合具體的實驗，驗證了基于ABAQUS和RBF神經(jīng)網(wǎng)絡的邊坡可靠度分析方法的可行性。

2)以具體的邊坡作為實例，利用土工實驗以及野外勘探所得參數(shù)和ABAQUS強度折減法計算結(jié)果構(gòu)造樣本數(shù)據(jù)，利用PSO算法對RBF神經(jīng)網(wǎng)絡模型進行優(yōu)化后擬合出功能函數(shù)；再通過Monte Carlo模擬法得到邊坡的失穩(wěn)概率和可靠指標。

3)對于c值和φ值的分布形式做出了合適的判斷，訓練的樣本數(shù)據(jù)包括了天然和飽和2種狀態(tài)，使得模型的準確度達到了99%以上，該方法和其他相關方法計算的結(jié)果相近，但利用了RBF神經(jīng)網(wǎng)絡訓練簡單，學習收斂速度快，能夠逼近任意非線性函數(shù)的優(yōu)點，構(gòu)造邊坡極限狀態(tài)響應面，擬合精度好，模型簡單直觀，計算效率更高；通過PSO優(yōu)化網(wǎng)絡使得計算精度相對于單一的RBF網(wǎng)絡法有進一步的提高，引入了可靠度分析方法，避免了傳統(tǒng)方法的一些缺點，對日后的滑坡風險評價起到技術(shù)支撐作用，在實際工程中具有一定實用價值。