趙 程，范宣梅，楊 帆，周 禮，郭 晨

(成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點實驗室，成都 610059)

圖1 白格滑坡區(qū)域地質(zhì)構(gòu)造圖

白格滑坡發(fā)生于西藏江達(dá)縣波羅鄉(xiāng)白格村與四川省白玉縣絨蓋鄉(xiāng)則巴村交界處，金沙江右岸，分別在2018年10月11日和11月3日發(fā)生兩次滑坡事件，堵塞金沙江，形成堰塞湖，在人工干預(yù)下于11月13日泄流完成，險情解除[1-2]。白格滑坡并非個例。近年來，中國西南地區(qū)地質(zhì)災(zāi)害頻發(fā)，其中具有破壞性大、高速和攜帶物質(zhì)多特點的滑坡碎屑流發(fā)生多次[3-4]。2009年6月5日，重慶市武隆縣鐵礦鄉(xiāng)雞尾山山體發(fā)生大規(guī)模碎屑流，造成10人死亡，64人失蹤，8人受傷[5]。2017年6月24日，四川省茂縣疊溪鎮(zhèn)新磨村后山山體發(fā)生碎屑流，掩埋了坡腳整個新磨村，致使10人死亡，73人失蹤[6-7]。郭晨等[8]通過無人機(jī)航測以及現(xiàn)場和室內(nèi)試驗，研究蔣劉4#滑坡呈近流體狀高速遠(yuǎn)程運動堆積特征。曾琇舒等[9]在野外調(diào)查和無人機(jī)航攝數(shù)據(jù)基礎(chǔ)上分析雷波縣白沙村滑坡-碎屑流的運動過程。目前常用于滑坡碎屑流的數(shù)值模擬軟件有DAN-W、DAN3D、Massflow、MassMov2D等。張睿曉等[10]以滑坡碎屑流為研究對象，通過模型試驗和數(shù)值模擬結(jié)果分析攔擋結(jié)構(gòu)效果。齊超等[11]通過滑坡動力分析軟件DAN-W建立了Voellmy模型、摩擦模型和F-V等3種模型對東河口滑坡-碎屑流進(jìn)行了模擬，得到了較好的模擬結(jié)果。高楊等[12]采用DAN3D數(shù)值方法對深圳人工堆填體滑坡運動過程進(jìn)行了模擬研究，堆積范圍、堆積厚度和運動速度與滑坡實際值基本吻合。夏式偉等[13]通過DAN3D軟件，運用Friction-Voellmy模型模擬了湯家溝滑坡-碎屑流，分析了湯家溝滑坡碎屑流運動過程。Salvatici等[14]通過DAN3D軟件，采用Voellmy模型得到碎屑流的流速和深度和實際情況有較好擬合度。Ouyang等[15]通過Massflow軟件對金沙江白格滑坡進(jìn)行了數(shù)值模擬，其運動過程和堆積形態(tài)與實際情況基本吻合。MassMov2D的模型是基于流體動力學(xué)方程的深度平均形式的二維有限差分解，流體被視為單相材料，其運動受流變學(xué)控制(即通過應(yīng)力與應(yīng)變之間的函數(shù)關(guān)系)，可以通過數(shù)字高程模型(DEM)建立復(fù)雜的地形條件用以實際事件的模擬[16]。通過數(shù)字代碼MassMov2D的開發(fā)和實現(xiàn)，可以在PCRaster GIS環(huán)境下模擬復(fù)雜地形條件下的泥石流事件[17-18]。Scaringi通過MassMov2D模型，對新磨村滑坡-碎屑流進(jìn)行了反演分析，重建了滑坡的運動特征和堆積形態(tài)，對災(zāi)害風(fēng)險評估有參考價值[19]。

通過使用MassMov2D模型對金沙江白格兩次滑坡事件進(jìn)行分析，對滑坡運動過程進(jìn)行了研究，基于對已發(fā)生的2次白格滑坡的模擬，采用參數(shù)反演，并對未來最不利情況，即滑坡物源區(qū)三處潛在不穩(wěn)定巖體同時失穩(wěn)進(jìn)行了模擬預(yù)測，獲得了運動過程與堆積范圍參數(shù)，結(jié)果表明三處潛在不穩(wěn)定巖體同時失穩(wěn)，會再次堵塞金沙江，形成堰塞湖，以期為滑坡的風(fēng)險評價與防災(zāi)減災(zāi)措施制定提供理論指導(dǎo)。

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

白格滑坡發(fā)生于金沙江右岸的山脊處，高程3 720～2 880 m，相對高差達(dá)840 m，河谷呈“V”形發(fā)育，該區(qū)域地貌類型為構(gòu)造侵蝕高山強(qiáng)烈寒凍風(fēng)化地貌，氣候類型為大陸性季風(fēng)高原型氣候，年平均氣溫7.7 ℃，平均降水量600 mm。

白格滑坡所處區(qū)域為金沙江構(gòu)造帶(圖1)，整個構(gòu)造帶呈NNW向展布。區(qū)域內(nèi)主要分布有雪青-龍崗斷裂(F1)、竹英-貢達(dá)斷裂(F2)、則巴-協(xié)塘斷裂(F3)、波羅-木協(xié)斷裂(F4)、貢達(dá)-迪中斷裂(F5)等斷裂和沙東-巴巴背形(M1)，其走向均為NW。其中波羅-木協(xié)斷裂(F4)的總體走向為330°左右，斷裂面向西傾斜，傾角50°～70°，從滑坡物源區(qū)頂部經(jīng)過，對滑坡的影響較大[1]。

2 白格滑坡基本特征

2.1 滑坡基本特征

白格滑坡兩次滑坡事件發(fā)生后，地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點實驗室在第一時間用無人機(jī)航拍獲得2018年10月12日和11月5日兩期影像數(shù)據(jù)，以及由四川省繪地理信息局提供的滑前1∶10 000數(shù)字高程模型(DEM)數(shù)據(jù)。通過以上數(shù)據(jù)處理后獲得了滑坡區(qū)數(shù)字地表模型(DSM)、數(shù)字正射影像圖(DOM)。通過ARCGIS從影像分析結(jié)合現(xiàn)場實地調(diào)查，將滑坡區(qū)域總體上分為三個區(qū)域：滑源區(qū)(Ⅰ)、流通區(qū)(Ⅱ)及堆積區(qū)(Ⅲ)；滑坡區(qū)域后緣及兩側(cè)巖土體在滑坡發(fā)生的影響下，具有明顯裂縫的三處潛在不穩(wěn)定巖體K1、K2、K3(圖2)[1]。

圖2 白格滑坡分區(qū)圖

2.2 滑坡運動特征

金沙江白格于2018年10月11日發(fā)生第一次滑坡事件，通過滑坡前后DEM差值計算其總方量約為1 960×104m3，第一次滑坡發(fā)生時，在滑源區(qū)約1 960×104m3的巖土體失穩(wěn)并向下運動，其運動過程中對坡面的鏟刮作用攜帶坡面巖土體下滑，形成碎屑流，到達(dá)金沙江后激起涌浪沖刷對岸坡面，在涌浪作用下對岸坡體部分物質(zhì)被沖到堆積區(qū)，造成金沙江堵塞形成堰塞湖，經(jīng)計算被鏟刮巖土體體積約240×104m3，第一次滑坡總方量為2 200×104m3，堆積體的最大厚度約為81.4 m[圖3(a)]。堰塞壩于2018年10月12日自然溢流，堰塞湖水位逐漸降低，13日險情解除。

第二次滑坡事件在2018年11月3日發(fā)生，在第一次滑坡的影響下，滑坡后緣K1區(qū)失穩(wěn)形成第二次滑坡，方量約為368×104m3的巖土體沿第一次滑坡形成的溝槽向下滑動，其運動過程中對坡面的鏟刮作用攜帶坡面巖土體形成碎屑流，再次堵塞金沙江，經(jīng)計算被鏟刮巖土體體積約312×104m3，包括滑源區(qū)巖土體、被鏟刮攜帶的巖土體以及第一次滑坡殘留在坡體上的松散堆積體，第二次滑坡的總方量為850×104m3，堆積體的最大厚度約為78.2 m[圖3(b)]。堰塞壩于11月9日在人工干預(yù)下泄洪。

圖3 滑坡堆積厚度圖

2.3 滑坡潛在不穩(wěn)定巖體特征

在第一次滑坡事件作用下，滑坡后緣及兩側(cè)三處潛在不穩(wěn)定巖體[圖2(a)K1、K2、K3區(qū)域][1]出現(xiàn)臨空面為第二次滑坡事件提供了有利條件，第二次滑坡事件對其擾動后，根據(jù)現(xiàn)場監(jiān)測結(jié)果顯示三處潛在不穩(wěn)定巖體仍然在持續(xù)變形[圖2(b)]，存在再次大規(guī)模下滑堵江的可能，因此后續(xù)通過數(shù)值模擬，對三處潛在不穩(wěn)定巖體失穩(wěn)做了預(yù)測研究。

3 白格滑坡數(shù)值模擬分析

3.1 MassMov2D模型描述

該模型的建立是基于Savage-Hutter理論，該理論假定流體為一種具有流變特性的單相均質(zhì)材料。使用流動力學(xué)方程的深度積分的近似，將流動建模為2-D連續(xù)介質(zhì)，這是泥石流數(shù)值模擬的常用方法[20-21]。與使用地形相關(guān)的局部參考系的深度積分模型不同，MassMov2D的控制方程在笛卡爾坐標(biāo)x、y的二維歐幾里得空間中被參考[22]。質(zhì)量和動量平衡方程為

(1)

(2)

(3)

(4)

(5)

根據(jù)朗肯理論中主動和被動兩個極值之間的范圍可知ka≤1≤kp，其取值由碎屑流混合物的內(nèi)摩擦角δ決定。

模型使用Voellmy模型，其庫侖摩擦阻力增加了湍流分量，該湍流分量與湍流因子(ξ)和流動質(zhì)量的厚度成反比，并且與流動速度的平方成正比[23]。模型通過以下關(guān)系描述：

(6)

式(6)中：τ為抗剪切應(yīng)力；g為重力加速度；v為滑體的運動速度；ρ為滑體的密度；f為摩擦系數(shù)；ξ為湍流系數(shù)。

3.2 滑坡模型參數(shù)取值

模擬采用Voellmy模型，當(dāng)湍流系數(shù)ξ的值很高時，湍流分量相對非常小，流變模型就縮減為庫倫模型[24]。模擬采用的參數(shù)在表1中給出，模擬使用5 m精度的數(shù)字高程模型，模型運行時間步長約束為Δt≤1 s。

3.3 模擬結(jié)果分析

3.3.1 10月11日第一次滑坡事件模擬

由表1的參數(shù)對金沙江第一次滑坡事件進(jìn)行模擬，分別選取10、20、30、60、80、100 s的運動過程的速度圖(圖4)、滑坡堆積厚度(圖5)及滑坡堆積厚度曲線(圖6)用來解釋其運動過程。

圖4 第一次滑坡事件模擬運動過程速度

表1 模型參數(shù)取值表

圖5 第一次滑坡事件堆積厚度

圖6 第一次滑坡事件堆積厚度曲線圖

第一次滑坡滑源區(qū)以及運動過程中鏟刮作用攜帶的巖土體總量約2 200×104m3，設(shè)置好滑源區(qū)碎屑流體進(jìn)行模擬。從圖5中可以看出，滑坡形成的碎屑流體很快的向下運動，于30 s時運動速度達(dá)到最高值為69 m/s，30～60 s滑坡速度逐漸減小，100 s時第一次滑坡事件結(jié)束。

第一次滑坡事件堆積體形態(tài)、堆積厚度集中區(qū)域均與實際情況較為相似(圖6)，第一次滑坡堆積情況在1-1′剖面中680 m處為最大堆積厚度達(dá)79.8 m，模擬堆積情況在1-1′剖面中560 m處為最大堆積厚度達(dá)86 m，堆積曲線整體上趨勢較為一致，模擬結(jié)果結(jié)合實際情況是在滑坡發(fā)生過程中有部分巖土體在江水沖刷帶到下游，模擬過程沒有考慮這一因素，故堆積厚度大于實際情況。

3.3.2 11月3日第二次滑坡事件模擬

由表1 的參數(shù)對金沙江第二次滑坡事件進(jìn)行模擬，分別選取10、20、30、60、80、100 s的運動過程的速度圖(圖7)、滑坡堆積厚度(圖8)及滑坡堆積厚度曲線(圖9)用來解釋其運動過程。

圖7 第二次滑坡事件模擬運動過程速度圖

圖8 第二次滑坡事件堆積厚度

圖9 第二次滑坡事件堆積厚度曲線圖

第二次滑坡是由于K1區(qū)巖土體失穩(wěn)，滑源區(qū)巖土體加上滑動過程中鏟刮巖土體總方量為850×104m3進(jìn)行模擬。從圖7中可以看出，滑坡形成的碎屑流體沿第一次滑動形成的溝槽向下運動，速度慢于第一次滑坡事件，于30 s時運動速度達(dá)到最高值為58 m/s，30～60 s滑坡速度逐漸減小，100 s時第二次滑坡事件基本結(jié)束。

第二次滑坡事件堆積體形態(tài)、堆積厚度集中區(qū)域均與實際情況較為相似(圖5)，第二次滑坡堆積情況在2-2′剖面中525 m處為最大堆積厚度達(dá)75.7 m，模擬堆積情況在2-2′剖面中530 m處為最大堆積厚度達(dá)76 m，堆積曲線整體上趨勢較為一致，堆積厚度最大區(qū)域匹配較好。

4 潛在不穩(wěn)定巖體預(yù)測

前文中提到兩次滑坡事件后，根據(jù)現(xiàn)場監(jiān)測結(jié)果顯示滑坡后緣及兩側(cè)有三處潛在不穩(wěn)定巖體[圖2(b)]，存在再次發(fā)生滑坡的可能[1]，因此對這三處潛在不穩(wěn)定巖體進(jìn)行模擬，使用第二次滑坡模擬反演得到的參數(shù)，預(yù)測其發(fā)生滑坡的運動速度、路徑以及堆積范圍和堆積厚度。

按照上述方法建立模型進(jìn)行模擬，考慮三處潛在不穩(wěn)定巖體同時滑動，整個運動和堆積過程持續(xù)約120 s。在30 s時滑坡速度達(dá)到最大值57 m/s，60～80 s滑坡速度逐漸降低，120 s時滑坡運動堆積過程基本結(jié)束(圖10)。通過120 s時滑坡堆積厚度最大為64 m。金沙江白格兩次滑坡事件最大堆積厚度分別為81.4 m和78.2 m，兩次滑坡事件均造成堵江。通過數(shù)值模擬得到的堆積厚度和形態(tài)來看，如果三處潛在不穩(wěn)定巖體同時下滑會再次形成堰塞體堵塞金沙江，此次模擬結(jié)果可為白格滑坡的防災(zāi)提供參考依據(jù)。

圖10 潛在不穩(wěn)定巖體模擬運動過程速度圖

圖11 潛在不穩(wěn)定巖體預(yù)測滑坡堆積厚度

5 結(jié) 論

2018年10月11日和11月3日在金沙江發(fā)生兩次滑坡事件，通過在PCRaster環(huán)境下使用MassMov2D模型對兩次滑坡進(jìn)行了模擬，得到以下結(jié)論。

(1)通過MassMov2D模型模擬金沙江白格兩次滑坡運動過程和堆積特征，結(jié)果顯示，第一次滑坡事件實際堆積最大厚度為81.4 m，模擬堆積最大厚度86 m；第二次滑坡事件實際堆積最大厚度78.2 m，模擬堆積最大厚度76 m，結(jié)果同真實情況較為接近。

(2)通過兩次滑坡反演得到的參數(shù)，用于對三處潛在不穩(wěn)定巖體的失穩(wěn)預(yù)測，預(yù)測結(jié)果顯示，三處潛在不穩(wěn)定巖體同時失穩(wěn)會形成滑坡，再次堵塞金沙江，堆積厚度達(dá)到64 m，因此應(yīng)在滑坡區(qū)布設(shè)監(jiān)測設(shè)備，實時監(jiān)測三處潛在不穩(wěn)定巖體變形情況，同時采取相應(yīng)的防治措施治理。

(3)通過對兩次滑坡事件數(shù)值模擬分析和潛在不穩(wěn)定巖體預(yù)測，說明MassMov2D可以較好地用于碎屑流滑坡的研究和分析，研究結(jié)果可為研究區(qū)災(zāi)害防治提供一定的參考依據(jù)。