徐奕梓， 樊曉一,2， 張友誼， 田述軍， 溫 翔， 劉浩南， 鄭榆楓， 廖洪陽

(1.西南石油大學(xué)土木工程與測繪學(xué)院，四川 成都 610500； 2.西南石油大學(xué)工程安全評估與防護(hù)研究院，四川 成都 610500； 3.西南科技大學(xué)土木工程與建筑學(xué)院，四川 綿陽 621010)

0 引言

四川省漢源縣地處橫斷山脈北段東緣，為川西高原與四川盆地之間的過渡地帶，區(qū)內(nèi)地形起伏大，地質(zhì)構(gòu)造復(fù)雜且地震活動強(qiáng)烈，是滑坡災(zāi)害的頻發(fā)區(qū)[1-2]。近年來，發(fā)生了多起地質(zhì)災(zāi)害，導(dǎo)致嚴(yán)重的人員傷亡和財產(chǎn)損失[3-5]。

災(zāi)難性滑坡常具有隱蔽性、突發(fā)性特征，現(xiàn)場監(jiān)測和實時獲取災(zāi)害參數(shù)存在較大的難度，模型試驗和數(shù)值模擬成為了理解和分析這類災(zāi)害的有效方法。滑坡運(yùn)動距離和堆積特征是滑坡防災(zāi)減災(zāi)和工程防護(hù)的主要技術(shù)參數(shù)，前人利用室內(nèi)實驗?zāi)Ｐ吞骄苛嘶逻\(yùn)動距離與滑坡體積之間的關(guān)系以及滑坡運(yùn)動距離的影響因素[6-7]，總結(jié)了體積、垂直高度、基底摩擦角、地形條件和不同偏轉(zhuǎn)角度等因素對堆積特性的影響[8-9]。計算機(jī)的發(fā)展為滑坡運(yùn)動過程的模擬帶來更多的可能性，DAN-W模型、離散元模型、DAN3D模型和Massflow模型是4種基于不同理論方法和本構(gòu)模型建立起來的模型，均能對滑坡失穩(wěn)破壞和運(yùn)動堆積過程進(jìn)行演化分析。國內(nèi)學(xué)者分別采用這4種模型對舟曲泄流滑坡、三溪村滑坡、水城滑坡、甘肅黑方臺滑坡的運(yùn)動特征、堆積特征、能量轉(zhuǎn)化特征和沖擊鏟刮效應(yīng)展開數(shù)值模擬，取得了較好的研究成果[10-13]。

降雨、地震、水庫蓄水、人類工程活動、斜坡結(jié)構(gòu)、特殊地層結(jié)構(gòu)、地質(zhì)構(gòu)造等是滑坡發(fā)生的主要因素[3,13-15]，地形條件、巖土體組成則影響了滑坡的運(yùn)動過程和致災(zāi)區(qū)域。2020年8月21日發(fā)生的中海村滑坡在地形條件、巖土體組成上不同于該區(qū)域的其他滑坡，本文通過現(xiàn)場調(diào)查并結(jié)合滑坡區(qū)域地形地貌、地層巖性和氣候水文條件，分析滑坡災(zāi)后特征，在此基礎(chǔ)上建立Massflow數(shù)值模型，反演整個滑坡動態(tài)過程并對滑坡動力學(xué)參數(shù)和運(yùn)動堆積形態(tài)進(jìn)行深入探討，旨在揭示影響滑坡啟動的關(guān)鍵因素和動力學(xué)特征，為該類滑坡運(yùn)動機(jī)理的研究提供參考和借鑒。

1 滑坡區(qū)地質(zhì)環(huán)境條件

1.1 地理位置及地形地貌

滑坡位于距漢源縣城約5 km的富泉鎮(zhèn)中海村，地處大渡河左岸，地貌類型屬于侵蝕構(gòu)造中低山地貌，海拔高程820～1 060 m，斜坡坡向N15°E，滑坡區(qū)兩側(cè)各發(fā)育1條沖溝，沖溝在坡體下部省道S306處匯聚。斜坡平均坡度17°，滑坡區(qū)分別在海拔高程940 m、1 000 m處分布2級緩坡臺地(圖1)。

圖1 研究區(qū)地理位置及地形地貌

1.2 地層巖性

滑坡區(qū)地層巖性單一，主要以第四系昔格達(dá)組半成巖為主(圖2)。半成巖是介于軟巖和松散土體之間的特殊巖體，含泥量較高，具有明顯的水平層理構(gòu)造。根據(jù)現(xiàn)場鉆孔結(jié)果顯示，中海村滑坡山體巖性在垂直方向自上而下依次為強(qiáng)風(fēng)化層—泥巖半成巖—泥巖(圖3)，強(qiáng)風(fēng)化層和泥巖半成巖主要呈黃褐色，泥巖呈灰白色。出露在滑坡區(qū)域周圍的地層由老至新依次有寒武系川盆西南筇竹寺組—婁山關(guān)組、中上二疊統(tǒng)喀大崩組+奔子欄組、上二疊統(tǒng)揚(yáng)子西南緣峨眉山玄武巖組+宣威組/黑泥哨組。

圖2 昔格達(dá)組分布示意圖[16]

圖3 滑坡現(xiàn)場巖石特征

1.3 降雨特征

漢源縣屬于亞熱帶季風(fēng)濕潤氣候，冬暖夏涼，四季分明。2010—2020年期間，研究區(qū)年平均氣溫16 ℃，年降雨量600 mm，其中2020年降雨量為858.9 mm (圖4a)。2020年全年降雨量不均勻分布，7—8月為雨季，降雨量分別為167.5 mm和351.5 mm，約占全年降雨量的60%(圖4b)。8月強(qiáng)降雨主要集中在滑坡發(fā)生的前7天，16日出現(xiàn)最大降雨，日降雨量達(dá)56.5 mm (圖4c)。經(jīng)實地調(diào)查，滑坡災(zāi)后現(xiàn)場表面出現(xiàn)了流水特征，認(rèn)為滑坡發(fā)生前受連續(xù)降雨的影響，雨水沿裂縫結(jié)構(gòu)面下滲使得動靜水壓力增加，飽和巖土體強(qiáng)度降低并引起應(yīng)力的變化，導(dǎo)致天然巖體結(jié)構(gòu)遭到破壞產(chǎn)生變形，而半成巖因成巖度不高且具遇水易軟化膨脹的特性，對于強(qiáng)風(fēng)化層來說滲透性較大，因此，強(qiáng)降雨和半成巖特性是導(dǎo)致滑坡災(zāi)害發(fā)生的關(guān)鍵因素。

(a) 年平均氣溫及降雨量 (b) 月降雨量 (c) 日降雨量

2 滑坡基本特征

2.1 滑坡規(guī)模及運(yùn)動演變特征

滑坡形態(tài)整體呈長條形狀，最大水平運(yùn)動距離約600 m，最大寬度約為260 m，前后高差約180 m(后緣高程1 035 m，前緣高程855 m)，滑坡運(yùn)動等效摩擦系數(shù)為0.3?；聠訒r，約15×104m3的山體在前期連續(xù)強(qiáng)降雨的影響下失穩(wěn)并快速滑出，隨即破碎解體，在上部緩坡平臺的下部，滑坡體向下加速運(yùn)動并鏟刮坡體上表面的松散土體，使體積逐漸增加至20×104m3； 運(yùn)動過程中受地形影響，運(yùn)動方向多次發(fā)生改變?；麦w沖擊掩埋了山坡平臺上和山腳的民房，部分滑坡邊緣區(qū)房屋嚴(yán)重受損，最終導(dǎo)致S360省道道路中斷，造成7人死亡2人失聯(lián)，8棟民房被毀(圖5)。

(a) 滑坡前 (b) 滑坡后

2.2 滑坡分區(qū)及運(yùn)動堆積特征

該滑坡發(fā)生在古滑體的殘坡積層[13]。對災(zāi)后現(xiàn)場的勘察顯示，滑坡區(qū)域呈階梯狀分布，通過當(dāng)?shù)鼐用瘾@悉，1989年滑坡源區(qū)后緣曾發(fā)生滑坡。在滑坡現(xiàn)場，利用無人機(jī)遙感技術(shù)生成滑后DEM數(shù)據(jù)，對滑坡前后高程變化進(jìn)行了差分計算(圖6)。

圖6 滑坡前后的高程變化

結(jié)合實地調(diào)查結(jié)果和高程變化特征，將致災(zāi)范圍沿滑坡運(yùn)動方向分為滑源區(qū)、堆積區(qū)Ⅰ、鏟刮區(qū)和堆積區(qū)Ⅱ共4個部分。

2.2.1 滑源區(qū)及堆積區(qū)Ⅰ

滑源區(qū)位于S360省道北側(cè)，主滑方向為N8°E，據(jù)現(xiàn)場測量，滑源區(qū)寬度約為165 m，滑坡后壁的最大下錯距離為14 m，且壁面完整平滑，形成了典型的“圈椅狀”滑坡壁，具有明顯的剪切摩擦痕跡(圖7)。野外滑坡痕跡及特點顯示，滑源區(qū)土體失穩(wěn)向下滑動，滑源區(qū)后部瞬間缺少支撐，受到巖土體向下運(yùn)動時的張拉和剪切作用，使得滑源區(qū)后緣部出現(xiàn)2條明顯的裂縫。第一條為直線形裂縫，裂縫高程為1 051 m，長度為25 m，裂縫發(fā)展方向為E40°S和W40°N，屬于張拉裂縫。第二條裂縫為弧形裂縫，裂縫高程為1 043 m，長度為11 m，裂縫發(fā)展方向為東南E30°S和W40°N。由于第二條裂縫的位置臨近滑源區(qū)，受到張拉力的同時也受到了滑源區(qū)坡體向下運(yùn)動時的剪切力，使得第二條裂縫出現(xiàn)下錯，最大下錯距離約為1 m，屬于剪切裂縫。

堆積區(qū)Ⅰ位于斜坡3級平臺上，堆積方向保持不變，區(qū)域內(nèi)最大堆積厚度約8 m。在滑坡發(fā)生后平臺上的一棟居民房屋完全被掩埋，導(dǎo)致2位老人失聯(lián)。

圖7 滑坡裂縫

2.2.2 鏟刮區(qū)

鏟刮區(qū)位于2級平臺上。部分滑坡體在堆積區(qū)Ⅰ開始堆積后，另一部分滑坡體沿著下墊面向下運(yùn)動，運(yùn)動方向發(fā)生改變，主滑方向為N24°E。該段滑坡下墊面場地為梯田，表層土體為強(qiáng)風(fēng)化層，由于常年受到人類活動的干預(yù)，土體較為松軟，且殘積坡積物較厚，滑坡體向下運(yùn)動時不斷鏟刮表面的松散土體使滑坡體體積增加，鏟刮表層土體的同時攜卷滑坡區(qū)兩側(cè)土體?；麦w向下鏟刮一段距離后由于地形的阻擋作用運(yùn)動方向再次偏轉(zhuǎn)，發(fā)生第二次堆積。

2.2.3 堆積區(qū)Ⅱ

第二次運(yùn)動方向偏轉(zhuǎn)后，即滑坡體進(jìn)入最終堆積階段(1級平臺)。主堆積方向為N22°W，沿堆積方向的最大水平距離約為120 m，最大堆積厚度為5 m。在堆積區(qū)Ⅱ存在坡度較大的陡坎(圖8)，滑坡體拋射加速后沖毀并掩埋了公路內(nèi)外側(cè)的房屋。根據(jù)滑坡的原始地形和堆積特征，中海村滑坡運(yùn)動地形為階梯型地形。階梯型滑坡相對于凹面型滑坡和坡腳型滑坡而言，最大水平運(yùn)動距離最小，但在坡腳位置具有最大沖擊力[17]。因此，雖然滑坡在堆積區(qū)Ⅱ運(yùn)動距離短，但導(dǎo)致了嚴(yán)重的房屋損壞和掩埋。

圖8 滑坡現(xiàn)場

3 滑坡數(shù)值模擬

3.1 Massflow介紹

動態(tài)反演分析是在獲取滑坡基本參數(shù)和地形數(shù)據(jù)條件下，借助數(shù)值分析工具進(jìn)行的動態(tài)特征分析。Massflow從流體力學(xué)的N-S方程出發(fā)，基于深度積分的連續(xù)介質(zhì)力學(xué)理論，提出具有二階精度的高效能數(shù)值仿真模型，其高效且擴(kuò)展自如的特點能讓軟件使用者根據(jù)自身需求使用Fortran語言對源代碼進(jìn)行添加和修改[18]。該軟件可用來模擬如山體滑坡、碎屑流、泥石流、堰塞湖、潰壩、雪崩等多種山區(qū)地質(zhì)災(zāi)害[19]，不僅可有效反演災(zāi)害動態(tài)過程，模擬結(jié)果還能為工程實踐提供參考。近年來，該方法在地質(zhì)災(zāi)害研究中取得了豐碩的成果，已應(yīng)用在如香港滑坡[20]、北京房山滑坡[21]、國外的諾拉泥石流[20]、國內(nèi)的紅椿溝泥石流[22]以及冷漬溝泥石流[23]的工程治理效果評價中。

3.2 模型參數(shù)選取

在Massflow已構(gòu)建的Coulomb-viscous模型、Voellmy模型、Manning模型這3種物理摩擦模型中，選取適用于滑坡、碎屑流和巖崩等地質(zhì)災(zāi)害的Coulomb-viscous模型，該摩擦模型需考慮的物理力學(xué)參數(shù)有： 密度ρ、黏聚力c、基底摩擦系數(shù)μ、孔隙水壓力系數(shù)ru、基底摩擦角δ和內(nèi)摩擦角φ。以上模擬參數(shù)主要根據(jù)文獻(xiàn)調(diào)研[24-25]和模擬校準(zhǔn)確定： 泥巖半成巖的密度范圍為1.8～2.3 g/cm3，取2.05 g/cm3，黏聚力c為28 kPa，內(nèi)摩擦角φ為25°，基底摩擦系數(shù)采用等效摩擦系數(shù)(f=0.3)，基底摩擦角則利用等效摩擦系數(shù)的反正切函數(shù)計算(δ約為17°)，孔隙水壓力系數(shù)控制在0.6～0.9范圍?？紤]到體積膨脹，將初始體積乘以系數(shù)1.2。

在孔隙水壓力系數(shù)分別為0.6、0.7、0.8和0.9情況下進(jìn)行試錯試驗，通過反向分析得到在孔隙水壓力系數(shù)為0.8時，模擬結(jié)果與實際情況更加吻合。為了量化分析數(shù)值模擬堆積厚度分布，將模擬結(jié)果在滑坡主滑方向、滑源區(qū)、堆積區(qū)Ⅰ和堆積區(qū)Ⅱ進(jìn)行分析(圖6)，使其與航空圖像和野外測量數(shù)據(jù)進(jìn)行對比驗證，對比結(jié)果如圖9所示。由于模型中沒有考慮表面微觀地貌的變化和房屋對滑坡體運(yùn)動的影響，剖面3-3′和4-4′的研究結(jié)果存在一定差異，但出現(xiàn)了較為一致的變化特色，因此，數(shù)值模擬結(jié)果具有可行性。

(a) 剖面1-1′ (b) 剖面2-2′

(c) 剖面3-3′ (d) 剖面4-4′

3.3 滑坡全運(yùn)動過程模擬

3.3.1 滑坡體運(yùn)動形態(tài)變化特征

采用Massflow模型的仿真模擬結(jié)果如圖10所示，以10 s間隔繪制不同時間步長的滑坡運(yùn)動形態(tài)，滑坡運(yùn)動過程共持續(xù)約70 s。為了更直觀地觀察不同時間厚度的變化特征，在主滑方向進(jìn)行剖切，從整個堆積過程可見，滑坡最大運(yùn)動距離約600 m，滑坡體主要堆積在3個平臺上，最大堆積厚度為8 m(圖11)。

圖10 中海村滑坡不同時刻運(yùn)動厚度分布

圖11 中海村滑坡不同時刻厚度分布

3.3.2 運(yùn)動速度特征

滑坡運(yùn)動速度的大小對整個滑坡致災(zāi)強(qiáng)度有著至關(guān)重要的影響。圖12中顯示了不同時間段滑坡體的速度分布特征。10 s時最大運(yùn)動速度位于滑坡體前部，20 s時最大運(yùn)動速度位于滑坡體中部，30 s時最大滑坡運(yùn)動速度位于滑坡體中后部； 在此期間最大運(yùn)動速度為16.2 m/s，均出現(xiàn)在3級平臺上。40～60 s階段，當(dāng)滑坡體運(yùn)動至1級平臺時滑坡體前緣速度開始減小，由于慣性作用，后緣顆粒的能量傳遞到前緣顆粒，使得前緣顆粒繼續(xù)向前運(yùn)動，掩埋了對側(cè)房屋，而此時滑源區(qū)物質(zhì)向3級平臺上繼續(xù)堆積，直到70 s時完全停止運(yùn)動。

圖12 中海村滑坡不同時刻運(yùn)動速度分布

在目前已有的眾多滑坡各區(qū)段運(yùn)動速度估算方法中，最為常見且有效的方法是謝德格爾(A. E. Scheidegger)法[26]，其計算公式如下：

。

(1)

式中：V為運(yùn)動速度，m/s；g為重力加速度,通常取9.8 m/s2；H為研究點距滑坡最高點的垂直距離,m；L為研究點距滑坡后緣的水平距離,m；f為等效摩擦系數(shù)。

經(jīng)計算，整個滑坡運(yùn)動速度曲線如圖13所示。對比2種方法分析的滑坡運(yùn)動速度，在整個運(yùn)動速度的演化過程中，滑源區(qū)以外的滑坡體運(yùn)動速度規(guī)律基本一致且出現(xiàn)3個峰值，峰值速度均出現(xiàn)在3個緩坡平臺的后部，其原因在于坡度由小到大的顯著變化致使運(yùn)動速度增加，從大到小的變化因滑體撞擊地面而速度減小。在滑坡整體的速度演化分布上，由于沒有考慮地形因素對速度的影響，導(dǎo)致計算值大于數(shù)值模擬值。

圖13 運(yùn)動速度對比

4 結(jié)論

(1)野外調(diào)查和地形高程變化特征顯示，中海村滑坡前后緣高差180 m，最大水平運(yùn)動距離600 m，總堆積體積為20×104m3。滑坡區(qū)域總體可分為滑源區(qū)、堆積區(qū)Ⅰ、鏟刮區(qū)和堆積區(qū)Ⅱ這4個區(qū)域。

(2)降雨和半成巖特性是導(dǎo)致滑坡災(zāi)害發(fā)生的關(guān)鍵因素?；滤趨^(qū)域地層巖性為第四系昔格達(dá)組泥巖半成巖，該巖層具有明顯的水平層理構(gòu)造，透水性強(qiáng)，遇水易軟化膨脹； 滑坡發(fā)生前的連續(xù)強(qiáng)降雨增加了地下水的入滲，誘發(fā)了滑坡。

(3)數(shù)值模擬結(jié)果表明，孔隙水壓力系數(shù)為0.8時模擬結(jié)果與實際情況吻合?；逻\(yùn)動持續(xù)時間約為70 s，最大堆積厚度為8 m，最大運(yùn)動速度為16.2 m/s。

(4)中海村滑坡速度在運(yùn)程上的分布演化存在3個峰值，峰值速度和主堆積區(qū)均分布在3個緩坡平臺上。