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

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

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

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

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

1 降雨入滲概化模型

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

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

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

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

圖1 降雨徑流示意圖

t時刻巖體中的含水量可表示為

Vt=Vt-1+ΔV

(1)

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

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

(2)

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

巖體中t時刻整體飽和度定義為

(3)

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

1.2 強度參數(shù)折減

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

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

(4)

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

1.3 強度參數(shù)折減

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

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

(5)

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

1.4 計算流程

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

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

2 遺傳算法

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

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

(6)

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

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

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

(7)

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

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

(8)

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

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

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

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

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

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

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

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

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

4.1 數(shù)值模型

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

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

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

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

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

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

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

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

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

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

圖6 王家臺剖面圖

圖7 塑性剪應(yīng)變云圖

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

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

圖8 現(xiàn)場監(jiān)測與數(shù)值計算位移曲線對比圖

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

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

5 結(jié)論

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

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

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