基于降雨入滲概化模型的王家臺滑坡巖土參數反演分析

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

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

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

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

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

1 降雨入滲概化模型

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

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

1.1 地質體中儲水量的計算

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

圖1 降雨徑流示意圖

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

Vt=Vt-1+ΔV

(1)

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

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

(2)

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

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

(3)

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

1.2 強度參數折減

當St=1時，開始強度參數折減，折減公式為

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

(4)

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

1.3 強度參數折減

當St=0時，開始強度參數增加，公式為

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

(5)

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

1.4 計算流程

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

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

2 遺傳算法

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

進行降雨參數識別的第一步就是定義一個目標函數，然后采用某些優化方法求解。采用的目標函數為

(6)

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

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

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

(7)

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

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

(8)

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

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

3 王家臺滑坡的地質特征

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

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

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

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

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

4 數值模型與計算工況

4.1 數值模型

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

4.2 計算工況及參數

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

圖4 王家臺滑坡地質模型

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

圖5 王家臺滑坡降雨數據

地層基本參數如表1所示。

表1 王家臺滑坡地層基本參數

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

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

圖6 王家臺剖面圖

圖7 塑性剪應變云圖

4.4 降雨工況下參數反演計算結果

基于提出的降雨入滲概化模型，以現場監測得出的0092和0093兩個監測點的地表位移曲線為優化目標，結合遺傳算法進行動態調參，最終數值計算得出的位移曲線和現場數據對比曲線如圖8所示。

圖8 現場監測與數值計算位移曲線對比圖

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

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

5 結論

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

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

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