唐岳灝 姜清輝

摘要：

為評估2018年金沙江白格村兩次特大型滑坡堰塞堵江事件后所形成潛在不穩(wěn)定塊體（殘留體）的穩(wěn)定性及堵江風險，在應用嚴格三維極限平衡分析白格滑坡殘留體穩(wěn)定性的基礎上，結(jié)合應急處置方案，基于SPH-DEM流固耦合模型，對殘留體滑坡范圍和堆積體厚度進行了預測，定量評估崩滑堵江風險。結(jié)果表明：① K1-Ⅰ和K2-Ⅰ滑塊處于臨界狀態(tài)，失穩(wěn)風險較大，存在再次發(fā)生滑坡堵江的可能；② 河道清淤方法可降低堰塞體整體高度，增大滑坡物質(zhì)堆積空間，有效降低再次滑坡堵江的風險。研究成果可為高位滑坡堵江風險預測提供一定參考。

關鍵詞：

白格殘留體； 滑坡堵江； SPH-DEM； 數(shù)值模擬； 風險評估

中圖法分類號：P642.22

文獻標志碼：A

DOI：10.15974/j.cnki.slsdkb.2023.05.006

文章編號：1006-0081（2023）05-0038-07

0 引 言

2018年10月10日和11月3日，在西藏自治區(qū)江達縣波羅鄉(xiāng)白格村與四川省白玉縣絨蓋鄉(xiāng)則巴村交界處的金沙江西藏岸，先后兩次發(fā)生大規(guī)模高位滑坡［1-4］，堵塞金沙江，形成堰塞湖?；卵呷蜎]了上游的村莊和各種生產(chǎn)生活設施，同時潰壩洪水嚴重沖毀滑坡下游的村莊、農(nóng)田和公路、橋梁等基礎設施［5-6］。白格滑坡經(jīng)兩次滑動后，松散物質(zhì)已大為減少，但滑坡體對山體擾動作用巨大（主要表現(xiàn)為拉、拽、刮、鏟、刷），滑坡三面由于臨空卸荷作用不斷加強，變形跡象十分明顯，后緣及兩側(cè)發(fā)育多條深大裂縫，存在大量滑坡殘留體。無人機航空攝影測量及工程地質(zhì)測繪、鉆探、物探等詳細勘察工作的成果表明：前兩次白格滑坡周界后緣有3個主要殘留變形體，其滑坡周界明顯，已經(jīng)發(fā)生了明顯的下滑、錯動，仍具備進一步下滑的趨勢。同時，地表宏觀變形跡象和監(jiān)測成果表明殘留體一直存在前緣溜滑和蠕滑變形，其穩(wěn)定安全問題不容樂觀。因此，開展白格滑坡殘留體的防災減災研究，分析殘留體的穩(wěn)定性和變形趨勢，可為滑坡應急處置提供支撐，并為白格滑坡第三次堵江的可能性、危害性預測提供基本數(shù)據(jù)。

目前，已有一些學者利用數(shù)值方法對白格殘留體堵江風險進行風險評價。曹水合等［7］利用數(shù)值模擬軟件Massflow對白格殘留體堵江范圍和高度進行了風險預測研究，趙程等［8］利用MassMov2D碎屑流模擬軟件對白格滑坡滑源區(qū)的3處潛在不穩(wěn)定巖體進行了預測，周禮等［9］利用PFC3D軟件對白格滑坡滑源區(qū)殘留不穩(wěn)定部分可能失穩(wěn)的運動路徑和堆積范圍進行預測。蔡耀軍等［10］利用PFC3D軟件對殘留體不同失穩(wěn)規(guī)模進行了堰塞體堆積形態(tài)預測。盡管這些數(shù)值仿真取得了一定進展，但也存在一些不足：這些軟件大多使用固體力學的方法模擬滑坡動力學過程，無法處理滑坡堵江產(chǎn)生的流固耦合問題；有的將滑坡過程簡化為流體流動過程，導致模擬結(jié)果與實際情況不完全對應。此外，上述大部分軟件都是基于國外商用程序，后處理的二次開發(fā)功能比較有限，在精細化展示堆積體的三維幾何形態(tài)方面缺乏可擴展性。

本文首先采用嚴格三維極限平衡分析方法分析計算滑體安全系數(shù)，同時確定滑體主滑方向，建立三維滑坡模型。然后使用基于拉格朗日力學體系的離散元（DEM）方法和光滑粒子流體動力學（SPH）方法完全自主開發(fā)的SPH-DEM耦合程序——利用DEM模擬大變形的滑坡動力過程、利用SPH模擬碎屑體入江產(chǎn)生的水動力學問題，以更好地反映滑坡（固體）與水體（流體）間的相互作用。通過SPH-DEM耦合程序可以實現(xiàn)真實場景下的滑坡堵江動力學過程模擬，對強變形區(qū)里不穩(wěn)定體崩落后的發(fā)展過程和堆積形態(tài)進行預測，定量評估崩落后再次堵江的風險，為潰壩及洪水災害分析做參考。

1 滑坡體特征

白格滑坡地質(zhì)剖面如圖1所示。根據(jù)勘察，按殘留體變形特征與空間位置關系，將殘留體劃分為3個區(qū)塊（圖2）：滑坡后緣為K1，滑坡左側(cè)（金沙江下游側(cè)）為K2，滑坡右側(cè)（金沙江上游）為K3。K1，K2和K3都是由滑坡牽引卸荷形成。各殘留體面積較大，根據(jù)不同變形邊界特征又可分為多個子塊，其中，K1殘留體細分為K1-Ⅰ，K1-Ⅱ和K1-Ⅲ，K2殘留體細分為K2-Ⅰ，K2-Ⅱ。

根據(jù)《西藏自治區(qū)昌都市江達縣波羅鄉(xiāng)金沙江白格滑坡應急勘查報告》：K1-Ⅰ和K1-Ⅲ殘留體潛在滑帶為元古界雄松群千枚巖組，巖性以絹云母石英千枚巖為主，為深灰色、灰黑色，變余泥質(zhì)結(jié)構(gòu)，千枚狀構(gòu)造，礦物成分以絹云母為主，含少量石墨、石英；K1-Ⅱ潛在滑帶為片麻巖與侵入蛇紋巖接觸面，巖性以全/強風化蛇綠巖為主，深灰綠色，風化呈淺綠-淺白色，碎塊石土狀；K2-Ⅰ殘留體潛在滑帶為元古界雄松群千枚巖組，巖性以絹云母石英千枚巖為主，深灰色、灰黑色，變余泥質(zhì)結(jié)構(gòu)，千枚狀構(gòu)造，礦物成分以絹云母為主，含少量石墨、石英；K2-Ⅱ殘留體潛在滑面不明顯；K3潛在滑帶為元古界雄松群片麻巖強風化層。典型巖石樣品如圖3所示。

2 滑坡后緣殘留體穩(wěn)定性評價

2.1 計算模型和計算方法

根據(jù)白格滑坡殘留體的地形等高線、滑坡邊界以及地質(zhì)鉆孔、深部測斜儀監(jiān)測確定的滑面深度，建立白格滑坡殘留體（K1-Ⅰ，K1-Ⅱ，K1，K2-Ⅰ，K2-Ⅱ，K3）的三維穩(wěn)定分析計算模型。計算方法采用滑源區(qū)嚴格三維極限平衡分析法，在對現(xiàn)有的邊坡穩(wěn)定性定量評價方法進行總結(jié)分析的基礎上，克服了現(xiàn)有滑坡穩(wěn)定性三維分析需要預先給定主滑方向以及滑體只滿足部分平衡條件的局限性，適用于多結(jié)構(gòu)面組合形成、具有復雜空間形態(tài)滑裂面的滑坡嚴格三維穩(wěn)定分析理論，可以同時計算得到滑體的安全系數(shù)和主滑方向，且嚴格滿足空間力系的6個平衡方程，能夠?qū)崿F(xiàn)對邊坡/滑坡在自然、降雨、地震不同工況條件下的穩(wěn)定狀態(tài)或加固后斜坡穩(wěn)定性的準確評估。圖4為計算所得的K1-Ⅰ、K1-Ⅱ和K2-Ⅰ滑面。

2.2 計算工況和參數(shù)取值

根據(jù)滑坡殘留體的實際情況，可能作用在殘留體上的荷載有：坡體自重、水壓力和地震力。① 坡體自重。在天然狀態(tài)下，坡體自重按天然重度計算。據(jù)當?shù)厮臍庀筚Y料及斜坡巖性組成的實際情況，在連續(xù)降雨或暴雨條件下，降水入滲深度按全入滲計，入滲范圍內(nèi)的土體重度按飽和重度計。② 水壓力。地下水產(chǎn)生的荷載主要是降雨誘發(fā)坡面暫態(tài)飽和區(qū)產(chǎn)生的地下水壓力等。③ 地震力。工區(qū)地震設防烈度為7度，設計地震分組為第二組，設計基本地震加速度值為0.20 g 。計算工況考慮3種：正常工況（自重荷載）、降雨工況（自重荷載+降雨水壓力）、地震工況（自重荷載+地震荷載）。殘留體滑帶力學計算參數(shù)主要依據(jù)巖土力學試驗和工程經(jīng)驗選取，如表1所示。

2.3 穩(wěn)定性評估與穩(wěn)定狀態(tài)分級

采用三維嚴格極限平衡法計算得到的白格滑坡后緣殘留體K1-Ⅰ，K1-Ⅱ，K1，K2-Ⅰ，K2-Ⅱ和K3在不同工況條件下穩(wěn)定性安全系數(shù)以及各個滑塊的體積、主滑方向和穩(wěn)定狀態(tài)，如表2所示，其中主滑方向角為與正北方向夾角，順時針為正。穩(wěn)定狀態(tài)根據(jù)《青藏高原重大滑坡風險防控指南》中滑坡穩(wěn)定狀態(tài)分為不穩(wěn)定 （K<1.0）；臨界穩(wěn)定（1.0≤K≤1.05）；弱穩(wěn)定（1.05

根據(jù)殘留體在正常、降雨、地震工況下的穩(wěn)定性分區(qū)、穩(wěn)定狀態(tài)分級和相應的預警等級，得到殘留體的風險等級（圖2），為下一步進行滑坡殘留體應急處置和堵江風險評估提供了直觀的參考。圖5為計算得到的各滑塊主滑方向與地表外觀測點位移矢量的對比，從中可以看出，主滑方向與外觀測點變形方向吻合，進一步論證了所使用的三維嚴格分析方法的有效性和可靠性。

基于對滑源區(qū)殘留體的穩(wěn)定性評估，可以看出K1-Ⅰ和K2-Ⅰ滑塊處于臨界狀態(tài)，失穩(wěn)風險較大，因此，首先考慮的失穩(wěn)模式為K1-Ⅰ+K2-Ⅰ發(fā)生崩滑（模式1）；其次考慮處于臨界狀態(tài)和弱穩(wěn)定狀態(tài)的滑塊發(fā)生失穩(wěn)，即K1-Ⅰ+K2-Ⅰ+K1-Ⅱ+K3發(fā)生崩滑（模式2）。

3 滑坡體應急處置

由于白格滑坡后緣殘留體K1-Ⅰ和K2-Ⅰ處于臨界穩(wěn)定狀態(tài)，存在再次下滑和堵江風險，因此，需要采取應急處置措施。對白格滑坡后緣殘留變形體的處置原則采用滑源區(qū)—運移區(qū)—堆積區(qū)全過程處理中的“顧兩頭”，分別為滑源區(qū)的處置和金沙江河道堰塞體的處置?；磪^(qū)應急處置采用的方案是削方減載，主要對K1-Ⅰ殘留體進行削坡，其目的是提高殘留體的穩(wěn)定性，降低滑坡風險，如圖6所示。削坡坡比根據(jù)GB 50330-2013《建筑邊坡工程技術(shù)規(guī)范》坡率法，按碎石土邊坡坡率允許值，結(jié)合深大裂縫分布位置，確定為1∶1.7～1∶1.3，分1～8級，坡高8 m，馬道寬5～6 m；總開挖體積約48萬m3。堰塞堆積區(qū)采用的方案是河道清淤，主要針對四川省一側(cè)開展，主要目的是降低堰塞體整體高度，提高水流下泄能力，增大滑坡物質(zhì)堆積空間，降低再次滑坡堵江的風險。具體方案為將堰塞體殘體中靠近河道側(cè)不小于40.0 m范圍內(nèi)、2 920 m高程以上的殘體全部清除。開挖設計方案共有4級馬道，馬道寬度0.2 m，放坡坡比1∶2.0，每一級臺階高度20 m，總開挖體積約245萬m3，如圖7所示。

4 滑坡殘留體堵江評估

4.1 光滑粒子流體動力學（SPH）

4.4 模型計算參數(shù)選擇

從DEM建模角度，需要將滑坡碎屑體概化為等直徑的小球。由于潛在滑坡體的體積超過500萬m3（模式2），如果采用小顆粒構(gòu)建滑坡體會導致顆粒數(shù)量巨大，嚴重影響DEM和SPH模型的計算效率，但過大的顆粒會影響模型計算精度。經(jīng)反復驗算，最終確定采用直徑2.4 m的球形顆粒。模型計算參數(shù)見表3。

4.5 堵江風險評估

采用滑坡啟動—運移—堵江災害鏈全過程模擬的SPH-DEM耦合方法，對白格滑坡后緣殘留體不同塊體組合形成的失穩(wěn)模式在失穩(wěn)觸發(fā)后的動力學過程及其入江后與水流相互作用的過程進行了模擬。根據(jù)滑坡最后形成的堰塞堆積體形態(tài)對滑坡堵江的危害性進行了評估。計算工況：① K1-Ⅰ+K2-Ⅰ（模式1）；② K1-Ⅰ+K2-Ⅰ+K1-Ⅱ+K3（模式2）。圖8～9為兩種失穩(wěn)模式最后形成的堰塞體堆積體厚度和橫剖面、縱剖面。為了對比分析，進行源區(qū)和河道工程應急處置后，相應失穩(wěn)模式造成的堰塞體堆積體厚度和橫剖面如圖10～11所示。表4為應急處置前后滑坡堰塞壩幾何特征統(tǒng)計。

根據(jù)圖8～11和表4分析所得結(jié)果如下。

（1） 應急處置前，兩種失穩(wěn)模式堆積體方量分別為440萬m3和527萬m3，最后堆積形成的堰塞壩高分別為60.2 m和63.0 m。相對模式1，模式2的滑塊方量增加了87萬m3，但最后形成堰塞壩的高度僅增加了約3 m。主要原因是模式2中的K1-Ⅱ滑塊和K3滑塊主要位于上游側(cè)坡，根據(jù)其失穩(wěn)觸發(fā)后運移軌跡，大部分滑坡物質(zhì)堆積到靠金沙江上游側(cè)河道，因此，模式2形成堰塞壩的寬度（順河向）與模式1相比增加了446 m，但壩高增加較少。

（2） 采用應急處置措施后，金沙江河道清淤245萬m3，為后緣殘留體失穩(wěn)后形成的滑坡壩提供了堆積空間，因此，模式1和模式2形成的堰塞壩高度分別為47.0 m和50.5 m，相對應急處置前壩高減少了約13 m，減少了堰塞湖庫容，有效降低了滑坡堰塞湖的危害性。

（3） 根據(jù)計算結(jié)果，按照最危險的情況（模式2），白格后緣殘留體處于臨界穩(wěn)定和弱穩(wěn)定的子滑塊觸發(fā)失穩(wěn)后，其堆積形成的最大壩高為50.5 m，與白格“10·10”滑坡堰塞壩高程大致相同?？紤]滑坡區(qū)交通條件大為改善、人工干預簡單易行等因素，裂縫區(qū)滑坡堵江及其災害鏈風險處于可控范圍。

5 結(jié) 論

（1） 本文采用嚴格三維穩(wěn)定性分析方法對白格后緣殘留體6個滑塊（K1-Ⅰ、K1-Ⅱ、K1、K2-Ⅰ、K2-Ⅱ、K3）在不同工況條件下的穩(wěn)定性進行評估，得到滑塊的安全系數(shù)和主滑方向，計算結(jié)果與現(xiàn)場宏觀定性分析和外觀監(jiān)測得到的各子區(qū)域變形方向吻合。根據(jù)滑塊穩(wěn)定性分析結(jié)果，對白格后緣殘留體進行了穩(wěn)定性分區(qū)、分級，并給出了分級監(jiān)測預警建議。

（2） 針對堆積區(qū)的河道清淤應急處置，采用滑坡災害鏈全過程模擬SPH-DEM方法，對后緣殘留體不同失穩(wěn)模式可能導致的滑坡堵江效果進行了分析。數(shù)值模擬結(jié)果表明：金沙江河道清淤為后緣殘留體失穩(wěn)后形成的滑坡壩提供了堆積空間，減少了堰塞湖庫容，有效降低了滑坡堰塞湖的危害性。

（3） 通過定量分析白格殘留體滑坡穩(wěn)定性與堵江風險，驗證了本文提出的方法可以用于碎屑流堵江的動力學過程及風險評價研究。

參考文獻：

［1］ 王立朝，溫銘生，馮振，等.中國西藏金沙江白格滑坡災害研究［J］.中國地質(zhì)災害與防治學報，2019（1）：1-9.

［2］ 馮文凱，張國強，白慧林，等.金沙江“10·11”白格特大型滑坡形成機制及發(fā)展趨勢初步分析［J］.工程地質(zhì)學報，2019（2）：415-425

［3］ 鄧建輝，高云建，余志球，等.堰塞金沙江上游的白格滑坡形成機制與過程分析［J］.工程科學與技術(shù)，2019，51（1）：9-16.

［4］ 王塞，鄧建輝，陳菲，等.“10·10”西藏白格滑坡運動特征反演分析［J］.哈爾濱工業(yè)大學學報，2022，54（2）：162-170.

［5］ 羅再均，次仁卓瑪，德吉白姆.2018年西藏堰塞湖搶險救災實踐及思考［J］.水利水電快報，2019，40（3）：18-21.

［6］ 謝忱，陳骎，侯奇東，等.金沙江白格堰塞壩自然泄流沖刷潰決過程數(shù)值模擬［J］.人民長江，2021，52（8）：22-29.

［7］ 曹水合，吳新明，鐘東.基于Massflow的金沙江白格滑坡數(shù)值模擬及堵江風險預測［J］.地質(zhì)災害與環(huán)境保護，2021（4）：3-7.

［8］ 趙程，范宣梅，楊帆，等.金沙江白格滑坡運動過程分析及潛在不穩(wěn)定巖體預測［J］.科學技術(shù)與工程，2020，20（10）：3860-3867.

［9］ 周禮，范宣梅，許強，等.金沙江白格滑坡運動過程特征數(shù)值模擬與危險性預測研究［J］.工程地質(zhì)學報，2019，27（6）：1395-1404.

［10］ 蔡耀軍，欒約生，楊啟貴，等.金沙江白格堰塞體結(jié)構(gòu)形態(tài)與潰決特征研究［J］.人民長江，2019（3）：15-22.

［11］ TANG Y，JIANG Q，ZHOU C.A Lagrangian-based SPH-DEM model for fluid-solid interaction with free surface flow in two dimensions［J］.Applied Mathematical Modelling，2018，62（10）：436-460.

［12］ 徐文杰.滑坡涌浪流–固耦合分析方法與應用［J］.巖石力學與工程學報，2020，39（7）：1420-1433.

（編輯：江 燾，高小雲(yún)）

Abstract：

In order to assess the stability and risk assessmentof river blockage for potentially unstable blocks （residual bodies） formed after two super-large landslides and river blockage occurred in Baige village，Jinsha River in 2018，the strict three-dimensional limit equilibrium method was applied to analyze the stability status of the residual bodies.Combing with the emergency response plans，the SPH-DEM fluid-solid coupling model was deployed to predict the deposit range and thickness，and to assess the risk of river blockage quantitatively.The results showed that：① The K1-I and K2-I blocks were in a critical state，with relatively large instability，leading to a high possibility of blocking Jinsha River again.② The clearance of river could decrease the overall height of the landslide damand increase the accumulation space for landslide debris consequently，which could effectively reduce the risk of river blockage.The research result could provide references for predicting the risk of landslide-caused river blockage.

Key words：

Baige residual body； river blockage； SPH-DEM； numerical simulation； risk assessment