基于降雨入滲概化模型的王家臺(tái)滑坡巖土參數(shù)反演分析

張新偉，朱心廣，馮春*，王心泉，劉新明

(1.北京市地質(zhì)工程勘察院，北京 100048；2.中國(guó)科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100190；3.中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院，北京 100049)

中國(guó)是滑坡災(zāi)害的多發(fā)國(guó)家，降雨降水、地震活動(dòng)、水庫(kù)蓄水等因素是誘發(fā)滑坡災(zāi)害的主要原因。其中，以降雨為誘因誘發(fā)的滑坡災(zāi)害占中國(guó)滑坡災(zāi)害總量的70%～90%[1-4]。降雨型滑坡主要有3個(gè)階段：第一階段是在降雨降水過(guò)程中雨水持續(xù)入滲邊坡體，邊坡體的自重不斷增加且孔隙水壓力持續(xù)升高，直至破壞邊坡體的極限平衡狀態(tài)并產(chǎn)生滑動(dòng)；第二階段是由于降雨降水產(chǎn)生的地下水位升高軟化了邊坡體的潛在滑動(dòng)面，并降低了邊坡體的抗剪強(qiáng)度；第三階段是在入滲和蒸發(fā)的反復(fù)作用下使邊坡體產(chǎn)生了節(jié)理裂隙，從而導(dǎo)致了邊坡體的失穩(wěn)滑動(dòng)[5-6]。

在降雨型滑坡中，國(guó)內(nèi)外學(xué)者的主要關(guān)注點(diǎn)是降雨入滲分析力學(xué)模型和坡體穩(wěn)定性評(píng)價(jià)方法。Hurley等[7]通過(guò)算例描述了地下水位在邊坡中的時(shí)空演化規(guī)律；李寧等[8]將非飽和土視覺(jué)幾何(visual geometry，VG)模型與修正后的Green-Ampt入滲模型相結(jié)合，提出了一個(gè)降雨誘發(fā)淺層滑坡的簡(jiǎn)化計(jì)算模型。陳曉東[9]將剛體極限平衡法與Mein-Larson入滲模型相結(jié)合，建立了一種降雨型邊坡穩(wěn)定性分析的通用分析模式；李龍起等[10]利用有PFC2D研究了竹林溝滑坡整體的失穩(wěn)模式，得出在持續(xù)降雨作用下邊坡中后緣表現(xiàn)為蠕滑-拉裂模式，應(yīng)力和位移整體上從前緣到后緣呈遞減趨勢(shì)；雷遠(yuǎn)見(jiàn)等[11]通過(guò)離散元方法和強(qiáng)度折減法的結(jié)合，討論了含多結(jié)構(gòu)面情況下巖質(zhì)邊坡的穩(wěn)定性。李世海等[12]采用連續(xù)-非連續(xù)單元法(continuum-discontinum element method，CDEM)，將滑坡體的地表裂縫作為分析的參數(shù)基本量，討論了滑坡體的運(yùn)動(dòng)性破壞。韓廷文等[13]在考慮地下水和地震力作用下，得到降雨條件下邊坡穩(wěn)定性及表面水平位移變化規(guī)律。王如賓[14]等對(duì)不同工況下降雨入滲引起的滑坡堆積體的滲流穩(wěn)定性進(jìn)行了數(shù)值分析。田東方等[15]以 Richards 方程和有限元法為基礎(chǔ)，通過(guò)數(shù)值算例討論了徑流補(bǔ)給對(duì)降雨型滑坡穩(wěn)定性的影響。童欣等[16]采用有限元方法，分析了降雨入滲對(duì)含軟弱夾層滑坡的穩(wěn)定性的影響。葉唐進(jìn)等[17]進(jìn)行了降雨條件下滾石滑坡失穩(wěn)機(jī)理的試驗(yàn)研究，推導(dǎo)出了適用于滾石斜坡降雨入滲條件的簡(jiǎn)步法穩(wěn)定性計(jì)算公式；賀可強(qiáng)等[18]以彈塑性理論和損傷力學(xué)為基礎(chǔ)，提出了降雨動(dòng)力增載位移響應(yīng)比的評(píng)價(jià)參數(shù)、預(yù)測(cè)模型和穩(wěn)定性判據(jù)。關(guān)曉迪等[19]通過(guò)數(shù)值計(jì)算和實(shí)驗(yàn)方法研究了不同坡比條件下黃土邊坡降雨入滲規(guī)律和機(jī)理。高謙等[20]利用GPS監(jiān)測(cè)和數(shù)值模擬對(duì)干坡和降雨工況下的邊坡穩(wěn)定性進(jìn)行了測(cè)試。王成華等[21]分析了雨強(qiáng)、坡角和孔隙比對(duì)非飽和砂土的非正交入滲規(guī)律。曾鈴等[22]利用有限元方法對(duì)降雨型邊坡的一維滲流及二維滲流過(guò)程進(jìn)行了計(jì)算。熊勇林等[23]基于有限元-有限差分理論，建立了水-土-氣三相滲流-變形耦合有限元數(shù)值分析。

由于地質(zhì)體具有非連續(xù)、非均勻、流固耦合、未知初始狀態(tài)等特點(diǎn)，僅依靠有限的鉆孔數(shù)據(jù)無(wú)法準(zhǔn)確給出滑坡體的滲流參數(shù)及巖土力學(xué)參數(shù)。為了準(zhǔn)確評(píng)價(jià)降雨條件下王家臺(tái)滑坡的穩(wěn)定性及成災(zāi)風(fēng)險(xiǎn)，現(xiàn)提出一種降雨入滲坡體的概化模型，通過(guò)將該模型內(nèi)嵌于連續(xù)-非連續(xù)單元法(CDEM)中，對(duì)降雨誘發(fā)的滑坡災(zāi)害進(jìn)行快速模擬；借助遺傳算法，通過(guò)監(jiān)測(cè)數(shù)據(jù)與數(shù)值模擬耦合的優(yōu)化分析，即可優(yōu)選出與現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)相匹配的王家臺(tái)滑坡當(dāng)前巖土參數(shù)。

1 降雨入滲概化模型

通常情況下，坡體的降雨計(jì)算涉及固體應(yīng)力場(chǎng)及孔隙滲流場(chǎng)的耦合計(jì)算，計(jì)算過(guò)程中需要采用真實(shí)時(shí)間模擬降雨過(guò)程，數(shù)值模型的計(jì)算時(shí)步受制于模型的滲流參數(shù)、網(wǎng)格尺寸和總計(jì)算時(shí)長(zhǎng)。為了精確描述坡體的流固耦合過(guò)程，重點(diǎn)區(qū)域的網(wǎng)格尺寸大都要求在米量級(jí)以下，且坡體中往往存在高滲透地層。因此，采用顯式耦合算法的求解步長(zhǎng)為秒量級(jí)甚至毫秒量級(jí)，而降雨分析的總時(shí)長(zhǎng)往往在年量級(jí)，從而導(dǎo)致計(jì)算耗時(shí)巨大。

基于上述問(wèn)題，提出了降雨入滲概化模型。該模型重點(diǎn)考慮降雨對(duì)地質(zhì)體強(qiáng)度的影響，忽略降雨引起的滲流壓力及降雨誘發(fā)入滲誘發(fā)地質(zhì)體破裂對(duì)坡體穩(wěn)定性的影響。該模型中需事先假定坡體的滲透層，而后將降雨量轉(zhuǎn)化為可滲透地層的飽和度，當(dāng)飽和度達(dá)到1后，開(kāi)始對(duì)可滲透地層的強(qiáng)度參數(shù)(黏聚力、抗拉強(qiáng)度、內(nèi)摩擦角)進(jìn)行折減，以達(dá)到降雨弱化巖體的目的。該簡(jiǎn)化模型中僅包含固體應(yīng)力場(chǎng)的計(jì)算，該部分計(jì)算可采用準(zhǔn)靜態(tài)算法，因此計(jì)算效率可以大幅度提高。

1.1 地質(zhì)體中儲(chǔ)水量的計(jì)算

假定地質(zhì)體中的雨水流入量為外部降雨，流出量分為地表徑流和地下徑流(圖1)。

圖1 降雨徑流示意圖

t時(shí)刻巖體中的含水量可表示為

Vt=Vt-1+ΔV

(1)

式(1)中：ΔV為t時(shí)刻巖體中雨水的增量，可表示為

ΔV=Rt-αVt-1-βVt-1

(2)

式(2)中：Rt為t時(shí)刻巖降雨量，m3/s；α為地表徑流系數(shù)；β為地下徑流系數(shù)。

巖體中t時(shí)刻整體飽和度定義為

(3)

式(3)中：n為巖體孔隙率，%；V為可透水巖體體積。

1.2 強(qiáng)度參數(shù)折減

當(dāng)St=1時(shí)，開(kāi)始強(qiáng)度參數(shù)折減，折減公式為

T=(T0-1-γDTmax)e-t1/ΔT+γDTmax

(4)

式(4)中：γD為恢復(fù)強(qiáng)度系數(shù)；Tmax為巖體原始強(qiáng)度值，Pa；T0-1為St=1時(shí)的初始時(shí)刻強(qiáng)度值，Pa；t1為巖體內(nèi)部飽和度為1時(shí)的持續(xù)時(shí)間，s；ΔT為強(qiáng)度改變所對(duì)應(yīng)的特征時(shí)間，s。

1.3 強(qiáng)度參數(shù)折減

當(dāng)St=0時(shí)，開(kāi)始強(qiáng)度參數(shù)增加，公式為

T=(T1-0-γUTmax)e-t0/ΔT+γUTmax

(5)

式(5)中：γU為恢復(fù)強(qiáng)度系數(shù)；T1-0為St=0時(shí)的初始時(shí)刻強(qiáng)度值，Pa；t0為巖體內(nèi)部飽和度為0時(shí)的持續(xù)時(shí)間，s。

1.4 計(jì)算流程

降雨快速計(jì)算模型的執(zhí)行流程如圖2所示。

圖2 降雨入滲概化模型流程圖

2 遺傳算法

采用遺傳算法進(jìn)行降雨入滲概化模型關(guān)鍵參數(shù)的反演優(yōu)化。遺傳算法利用簡(jiǎn)單的編碼技術(shù)和進(jìn)化繁殖機(jī)制解決復(fù)雜的多極值優(yōu)化問(wèn)題。

進(jìn)行降雨參數(shù)識(shí)別的第一步就是定義一個(gè)目標(biāo)函數(shù)，然后采用某些優(yōu)化方法求解。采用的目標(biāo)函數(shù)為

(6)

式(6)中：um為地表觀測(cè)位移；uc為數(shù)值計(jì)算位移；α、β、n、γD、γU、ΔT為降雨參數(shù)。

采用改進(jìn)型的十進(jìn)制遺傳算法進(jìn)行動(dòng)態(tài)尋優(yōu)，該算法包括重組、變異、選擇和精英搜索4個(gè)步驟。

步驟1重組。重組是最主要的遺傳運(yùn)算，它同時(shí)對(duì)兩個(gè)個(gè)體操作，組合二者的特性產(chǎn)生新的后代。采用如下所示的凸交叉，在種群中隨機(jī)選擇兩個(gè)個(gè)體xα和xβ，經(jīng)過(guò)交叉后產(chǎn)生兩個(gè)新個(gè)體，表達(dá)式為

(7)

式(7)中：μ為服從[0，1]的均勻分布。

步驟2變異。變異是另一種最基本的運(yùn)算，它在種群中個(gè)體上自發(fā)地產(chǎn)生隨機(jī)的變化。一種簡(jiǎn)單的變異方法是隨機(jī)替換一個(gè)或多個(gè)個(gè)體。在遺傳算法中，變異可以提供初始值其中不含有的個(gè)體，或找回在重組和選擇中丟失的個(gè)體，為種群提供多樣性。

(8)

步驟3選擇。選擇的作用是優(yōu)勝劣汰、適者生存。將多個(gè)個(gè)體值代入式(6)，將值最大的個(gè)體淘汰。

步驟4精英搜索。將較差的個(gè)體淘汰后剩下的即為精英個(gè)體，可進(jìn)行下一輪的遺傳算法尋優(yōu)。

3 王家臺(tái)滑坡的地質(zhì)特征

王家臺(tái)滑坡位于北京市房山區(qū)霞云嶺鄉(xiāng)王家臺(tái)村西南方向130 m左右處(圖3)。所屬地質(zhì)災(zāi)害類(lèi)型為不穩(wěn)定斜坡，長(zhǎng)約143 m，寬為30～50 m，面積約6 000 m2，坡高約50 m，坡度約40°，分布高程516～531 m。地勢(shì)為西部高，東部低。

圖3 王家臺(tái)滑坡航拍圖(1∶1 300)

不穩(wěn)定斜坡西部植被覆蓋一般，主要以灌木為主，植被覆蓋率約60%。斜坡東部呈階梯狀，平臺(tái)與平臺(tái)之間為近直立的陡坎，平臺(tái)處種有玉米、南瓜和大樹(shù)。每個(gè)陡坎處均砌有高1.0～2.5 m的干砌石擋土墻，大約有20級(jí)，結(jié)構(gòu)松散、無(wú)定向性。斜坡底部為108國(guó)道，一旦發(fā)生滑塌，將會(huì)危及108國(guó)道。

該滑坡周邊出露地層為上元古界薊縣系、更新統(tǒng)坡洪積物(Qp3dpl)和全新統(tǒng)沖洪堆積物(Qhapl)，巖層產(chǎn)狀變化不大。大部分山體坡面較為穩(wěn)定。沉積物為礫石、砂、砂質(zhì)黏土、亞砂土，以前兩者為主。區(qū)域大地構(gòu)造位置處于華北板塊中部，燕山臺(tái)褶帶中的次一級(jí)構(gòu)造單元西山迭拗褶帶西部。區(qū)域地質(zhì)構(gòu)造復(fù)雜，處于后呂梁-印支期地臺(tái)蓋層褶皺中的第二亞構(gòu)造層，其西北部有一條張剪性斷裂，東南部位有一條推覆構(gòu)造。

據(jù)斜坡的現(xiàn)狀空間幾何特征、物質(zhì)組成及所處的地質(zhì)環(huán)境分析，不穩(wěn)定斜坡坡體破壞形式主要為滑塌，斜坡主要成分是第四系人工梯田和坡積碎石覆蓋于青白口系下馬嶺組(QNx)泥質(zhì)粉砂巖之上，斜坡因開(kāi)挖形成陡坎，并且斜坡整體植被保持水土能力較差，如遇暴雨天氣，斜坡土體受雨水侵蝕，易發(fā)生滑塌。斜坡主要影響因素包括內(nèi)部條件和外部條件兩類(lèi)，內(nèi)部條件包括地形地貌、坡體結(jié)構(gòu)；外部條件包括降雨、風(fēng)化、地震、人類(lèi)工程活動(dòng)等。

4 數(shù)值模型與計(jì)算工況

4.1 數(shù)值模型

建立如圖4所示的王家臺(tái)數(shù)值計(jì)算模型。模型后緣高度67 m，前緣高度30 m，西向長(zhǎng)度76 m，北向長(zhǎng)度115 m。模型共包含兩類(lèi)地層，其中基巖為泥質(zhì)粉砂巖，覆巖為破碎積石。前緣坡腳處為108國(guó)道。1-1′、2-2′、3-3′為截取的3個(gè)橫截面。

4.2 計(jì)算工況及參數(shù)

反演目標(biāo)是王家臺(tái)滑坡的巖土參數(shù)，具體包括地表徑流系數(shù)、地下徑流系數(shù)、巖體孔隙率、殘余強(qiáng)度系數(shù)、恢復(fù)強(qiáng)度系數(shù)、強(qiáng)度變化特征時(shí)間；收斂條件是0092及0093兩個(gè)地表位移監(jiān)測(cè)點(diǎn)(圖4)的實(shí)測(cè)曲線與數(shù)值模擬曲線的誤差小于20%。

圖4 王家臺(tái)滑坡地質(zhì)模型

在2019年3月—2020年7月王家臺(tái)滑坡降雨觀測(cè)站記錄的460 d降雨數(shù)據(jù)如圖5所示。

圖5 王家臺(tái)滑坡降雨數(shù)據(jù)

地層基本參數(shù)如表1所示。

表1 王家臺(tái)滑坡地層基本參數(shù)

4.3 正常工況(非降雨)下安全系數(shù)計(jì)算

如圖6所示，截取圖4中3個(gè)橫截面，利用GDEM軟件進(jìn)行非降雨工況下王家臺(tái)滑坡災(zāi)害點(diǎn)的安全系數(shù)分析。圖7所示為塑性剪應(yīng)變?cè)茍D，針對(duì)截面1-1′，塑性區(qū)域位于地層分界處，處于邊坡中部坡度較陡處，邊坡安全系數(shù)為1.46；截面2-2′中存在邊坡上部與邊坡中部?jī)蓚€(gè)塑性區(qū)，邊坡安全系數(shù)為1.682；截面3-3′的塑性區(qū)位置與截面2-2′相同，安全系數(shù)為1.45。塑性區(qū)所在位置處為潛在滑動(dòng)危險(xiǎn)區(qū)，結(jié)合3個(gè)截面計(jì)算情況，王家臺(tái)滑坡災(zāi)害點(diǎn)綜合安全系數(shù)為1.53，現(xiàn)場(chǎng)監(jiān)測(cè)安全系數(shù)為1.50，計(jì)算誤差2%。

圖6 王家臺(tái)剖面圖

圖7 塑性剪應(yīng)變?cè)茍D

4.4 降雨工況下參數(shù)反演計(jì)算結(jié)果

基于提出的降雨入滲概化模型，以現(xiàn)場(chǎng)監(jiān)測(cè)得出的0092和0093兩個(gè)監(jiān)測(cè)點(diǎn)的地表位移曲線為優(yōu)化目標(biāo)，結(jié)合遺傳算法進(jìn)行動(dòng)態(tài)調(diào)參，最終數(shù)值計(jì)算得出的位移曲線和現(xiàn)場(chǎng)數(shù)據(jù)對(duì)比曲線如圖8所示。

圖8 現(xiàn)場(chǎng)監(jiān)測(cè)與數(shù)值計(jì)算位移曲線對(duì)比圖

分析曲線可知，0092測(cè)點(diǎn)上升期(前150 d)數(shù)值計(jì)算結(jié)果相比監(jiān)測(cè)結(jié)果偏低，后期與監(jiān)測(cè)結(jié)果基本一致。0093測(cè)點(diǎn)上升期(前150 d)數(shù)值計(jì)算結(jié)果與監(jiān)測(cè)結(jié)果吻合較好，后期相比監(jiān)測(cè)結(jié)果略偏大?？傮w而言，數(shù)值計(jì)算結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果總體趨勢(shì)一致，吻合較好，這也說(shuō)明了提出的降雨入滲概化模型的合理性和有效性。

在給定的初始巖土參數(shù)，結(jié)合降雨入滲概化模型和遺傳算法進(jìn)行王家臺(tái)滑坡當(dāng)前巖土參數(shù)的尋優(yōu)計(jì)算；最終給出與降雨入滲概化模型相匹配的巖土參數(shù)：地表徑流系數(shù)α=0.1，地下徑流系數(shù)β=0.05，巖體孔隙率n=0.01，殘余強(qiáng)度系數(shù)γD=0.2，恢復(fù)強(qiáng)度系數(shù)γU=0.8，強(qiáng)度改變系數(shù)所對(duì)應(yīng)的時(shí)間ΔT=5×106s。

5 結(jié)論

(1)提出了可快速分析降雨入滲誘發(fā)滑坡災(zāi)害的力學(xué)模型。該模型根據(jù)降雨量、地表徑流量及地下徑流量綜合確定滑體中的飽和程度，進(jìn)而根據(jù)飽和度來(lái)折減或增加滑體的強(qiáng)度。由于不需要進(jìn)行滲流應(yīng)力耦合的計(jì)算，因此可大幅提升計(jì)算效率。

(2)運(yùn)用GDEM軟件，截取王家臺(tái)滑坡災(zāi)害點(diǎn)特定剖面進(jìn)行正常工況下安全系數(shù)分析，得出綜合安全系數(shù)為1.53，相比于現(xiàn)場(chǎng)監(jiān)測(cè)誤差為2%。

(3)基于提出的降雨入滲概化模型，結(jié)合遺傳算法，實(shí)現(xiàn)了王家臺(tái)滑坡巖土參數(shù)的反演分析。經(jīng)過(guò)監(jiān)測(cè)數(shù)據(jù)與數(shù)值模擬的優(yōu)化對(duì)比，給出了適用于降雨入滲概化模型的王家臺(tái)滑坡當(dāng)前巖土參數(shù)為：地表徑流系數(shù)α=0.1，地下徑流系數(shù)β=0.05，巖體孔隙率n=0.01，殘余強(qiáng)度系數(shù)γD=0.2，恢復(fù)強(qiáng)度系數(shù)γU=0.8，強(qiáng)度改變系數(shù)所對(duì)應(yīng)的時(shí)間ΔT=5×106s?；谏鲜鰠?shù)，可進(jìn)一步開(kāi)展王家臺(tái)滑坡在不同降雨條件下的成災(zāi)風(fēng)險(xiǎn)分析。