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

(中國(guó)計(jì)量大學(xué) 災(zāi)害監(jiān)測(cè)技術(shù)與儀器國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，浙江 杭州 310018)

0 引言

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

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

1 分析方案

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

圖1 分析方案設(shè)計(jì)Fig.1 Design of analysis scheme

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

2.1 ABAQUS強(qiáng)度折減法

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

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

cm=c/Fr

(1)

φm=arctan(tanφ/Fr)

(2)

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

使用強(qiáng)度減法所選取的屈服準(zhǔn)則[8]是Mohr-Coulomb：

(3)

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

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

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

(4)

隱含層到輸出層映射為：

(5)

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

(6)

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

(7)

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

2.3 PSO算法

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

(8)

(9)

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

(10)

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

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

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

2)計(jì)算每個(gè)粒子的適應(yīng)度，如式(11)：

(11)

式中：yi和fi表示實(shí)測(cè)值和預(yù)測(cè)值。

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

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

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

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

2.4 Monte Carlo模擬方法

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

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

(12)

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

(13)

(14)

用β表示可靠度指標(biāo)，則β可表示為：

β=μz/σz

(15)

則失效概率：

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

(16)

式中：Φ(β)為標(biāo)準(zhǔn)正態(tài)分布。

3 實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

圖6 FV1隨U1的變化關(guān)系Fig.6 Variation relationship of FV1 with U1

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

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

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

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

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

3.2 訓(xùn)練結(jié)果分析

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

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

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

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

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

3.3 蒙特卡羅模擬計(jì)算

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

X=μ+σx′

(17)

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

(18)

(19)

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

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

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

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

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

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

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

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

3.4 不同方法對(duì)比分析

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

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

4 結(jié)論

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

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

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