高位滑坡動力學(xué)特征分析

徐楊青，喬龍騰，楊龍偉，江強強，王孝臣

(1.煤炭科學(xué)研究總院，北京 100013；2.中煤科工集團武漢設(shè)計研究院有限公司，湖北 武漢 430064)

近年來，受到全球極端氣候變化和人類工程活動等影響，高位滑坡地質(zhì)災(zāi)害呈現(xiàn)頻發(fā)狀態(tài)，給人民生命財產(chǎn)造成巨大威脅[1]，如2017年6月24日發(fā)生的四川茂縣新磨村滑坡造成83人死亡[2]，2019年7月23日發(fā)生的貴州省水城縣雞場鎮(zhèn)坪地村滑坡造成43人死亡[3]，2021年2月發(fā)生的印度北阿坎德邦查莫利高位冰川崩裂災(zāi)害造成超過200多人失蹤[4]。由于高位滑坡地質(zhì)具有急劇啟動，破壞性大和運動范圍廣等特點，其防災(zāi)減災(zāi)難度較大，一直是工程地質(zhì)災(zāi)害領(lǐng)域研究的熱點和痛點[5]。

對于高位滑坡災(zāi)害定義目前尚未統(tǒng)一，大量研究成果表明，其主要是指滑坡剪出口與滑坡前緣的垂直高差大于50m，且滑源區(qū)多分布在高山峽谷等地區(qū)，植被茂密，這導(dǎo)致其滑坡演化過程十分隱蔽[6]。目前圍繞高位滑坡災(zāi)害的研究主要聚焦于成災(zāi)機制、滑坡動力學(xué)和遠程運動堆積等方面，其中高位滑坡災(zāi)害動力學(xué)是分析高位滑坡動力災(zāi)害效應(yīng)的重中之重，這是因為滑坡動力學(xué)計算分析得到的滑體速度、堆積體厚度和最遠運動距離等都可為滑坡防災(zāi)減災(zāi)等提供重要參考[7，8]。廣大專家學(xué)者圍繞滑坡動力災(zāi)害學(xué)開展了卓有成效的工作，殷躍平等基于滑塊-彈簧模型的基礎(chǔ)上提出了高位遠程滑坡動力侵蝕犁切模型，并計算了滑體運動過程中受到的犁切阻力[9]。胡偉等通過現(xiàn)場調(diào)查和室內(nèi)試驗發(fā)現(xiàn)雞尾山高位滑坡的底部摩擦溫度高達790℃，提出了高位遠程滑坡滑帶高溫?zé)岱纸鈾C理[10]。高楊等基于Hertz接觸模型和等效流體理論計算分析了動力沖擊鏟刮對高位滑坡體積效應(yīng)等影響，并提出了四種動力沖擊鏟刮模式[11]。嚴炎和崔一飛等通過對高位滑坡災(zāi)害震動信號分析，反演了高位滑坡運動過程[12]。羅剛和程謙恭等總結(jié)了高位高能巖崩災(zāi)害的動力沖擊力及工程防護措施等[13]。邢愛國基于離散元模型方法計算了貴州納雍高位崩塌運動過程及其震動信號特征[14]。Zhang Ming探討了重慶雞尾山滑坡的巖體強度降低可能與礦物質(zhì)分解有關(guān)，進而誘發(fā)了滑坡急劇啟動[15]。Watkins通過遙感技術(shù)分析了火星上的遠程滑坡運動特征，發(fā)現(xiàn)水化硅酸鹽潤滑基底滑動帶這一機制起到重要作用[16]。上述研究為高位滑坡動力學(xué)災(zāi)害研究提供了重要指導(dǎo)方向，也為開展類似研究提供有力支撐。

基于此，本文探討了高位滑坡災(zāi)害運動學(xué)特征，并以菲律賓 Guinsaugon 滑坡為例，通過基底流變模型動力學(xué)方法計算分析滑坡運動速度、堆積體厚度和最遠運動距離等，并與雪橇模型計算結(jié)果進行對比，最后總結(jié)該滑坡的動力學(xué)過程，為類似的高位滑坡地質(zhì)災(zāi)害研究提供借鑒。

1 菲律賓 Guinsaugon 滑坡概況

2006年2月17日，菲律賓萊特島Saint Bemard鎮(zhèn)Guinsaugon村發(fā)生高位滑坡地質(zhì)災(zāi)害，運動過程中掩埋一所小學(xué)，此次災(zāi)害共造成1800人死亡，震驚世界[17-19]。事后調(diào)查發(fā)現(xiàn)，該研究區(qū)地處Philippine 斷層，其活動速率達到2.5cm/a，且該地區(qū)地震活躍，這些都加劇了巖體破碎程度。該地區(qū)降雨豐沛，月均降雨量在600～800mm，滑坡發(fā)生前該地區(qū)出現(xiàn)連續(xù)16d降雨，累計達到700mm。由此可見，地質(zhì)營力活躍、巖體結(jié)構(gòu)破碎和降雨是該滑坡發(fā)生的主要誘發(fā)因素?；磪^(qū)體積約為4×106m3，滑坡后緣高程為800m，剪出口高程約為260m，前緣高程為30m，最遠距離為4100m，呈現(xiàn)出現(xiàn)典型的高位遠程滑坡災(zāi)害特征，詳情如圖1所示[17-19]。根據(jù)Guinsaugon滑坡的運動堆積和空間分布特征，將其劃分為滑源區(qū)、碎屑流區(qū)和堆積區(qū)。

滑源區(qū)：滑源區(qū)海拔高程范圍為380～800m。母巖巖性主要為角礫巖、碎屑巖和沉積巖等為主。巖體發(fā)育多組反傾節(jié)理，這些與巖層面一起將巖體切割成不對稱“楔形體”，高度約為420m。在滑坡發(fā)生前，巖體后緣裂縫已貫通至基巖處，為滑坡發(fā)生提供了優(yōu)勢入滲通道。

碎屑流區(qū)：碎屑流區(qū)位于海拔50～380m，高差約為330m，斜長約為2027m。根據(jù)現(xiàn)場目擊者描述，滑體撞擊坡面并散落堆積，順著原沖溝進行運動，呈現(xiàn)管道型碎屑流。

堆積區(qū)：堆積區(qū)坡度較緩，其水平長度約為1000m，其前緣已接近河流，平面形狀呈現(xiàn)“扇形”。由于其運動路徑上有稻田覆蓋，降低了滑體運動過程中受到的摩擦阻力，加速了滑體的遠程運動。

圖1 Guinsaugon滑坡

2 高位滑坡運動學(xué)特征分析

高位滑坡地質(zhì)災(zāi)害由于其滑源區(qū)多位于脊狀山頂，在降雨和地震等誘發(fā)因素的作用下，巖體結(jié)構(gòu)進一步破碎，加速滑坡失穩(wěn)。高位滑坡急劇啟動后，由于滑坡前后緣間的垂直落差大，經(jīng)過重力勢能轉(zhuǎn)化成動能，導(dǎo)致滑體的運動速度加快，平均速度大于20m/s，也使得滑體運動的最遠距離一般超過1km，形成典型的高位遠程滑坡災(zāi)害，極易造成群死群傷事件。

在運動過程中，根據(jù)滑坡運動過程的空間演化特征及速度變化趨勢，將滑坡運動路徑主要分為滑源區(qū)，碎屑流區(qū)和堆積區(qū)?；磪^(qū)由于植被覆蓋茂密，其巖體裂化演化過程十分隱蔽，這使得滑坡發(fā)生具有突發(fā)性。此外，滑體在剪切啟動過程中釋放的能量也加劇了滑坡運動動能，使得滑體在啟動后，短時間內(nèi)速度增長至2～3m/s。碎屑流區(qū)內(nèi)地表多覆蓋松散堆積物，這使得滑體在運動至該區(qū)域時會鏟刮并裹挾地表松散堆積物，實現(xiàn)了滑坡“體積效應(yīng)”。由于碎屑流區(qū)內(nèi)坡度逐漸由陡變緩，且碎屑流區(qū)內(nèi)地表微地貌起伏變化，這導(dǎo)致滑體在該區(qū)域內(nèi)的速度波動起伏，平均速度在20m/s以上。在堆積區(qū)內(nèi)，由于滑體在運動過程中受到河流、對岸山體等阻礙，其速度逐漸變緩直至停止。從碎屑流區(qū)到堆積區(qū)的過程中，由于地形逐漸開闊，滑體從碎屑流運動狀態(tài)漸變至散落堆積，由于滑體之間的相互擠壓碰撞，使得滑體速度呈現(xiàn)差異化，進而在散落堆積過程中形成定向排列、橫向脊和縱向脊等堆積地貌特征。

綜上所述，高位滑坡孕災(zāi)過程隱蔽，急劇發(fā)生、動力沖擊和遠程成災(zāi)等特點。加強高位滑坡動力學(xué)過程分析，對于高位滑坡風(fēng)險防控顯得尤為重要。

3 高位滑坡的動力學(xué)模型

3.1 模型原理

高位滑坡在快速運動過程中，與地表產(chǎn)生巨大碰撞并粉碎成塊狀，進行形成滑坡碎屑流，鏟刮地表松散堆積物，使得滑體體積增大。本研究主要利用等效流體計算方法，主要是基于拉格朗日解析解來計算組成滑體的一系列塊體在運動過程中受到的基底剪切阻力(圖2和圖3)[20]。在刻畫高位滑坡動力學(xué)過程中，F(xiàn)riction模型(F)和Voellmy模型(V)被廣泛運用，其中Friction模型可以較好地闡釋塊體受到的基底阻力，Voellmy模型適用于濕型滑坡泥石流的流變特性。本文的研究重點是高位滑坡碎屑流的動力學(xué)特征，且受到連續(xù)降雨影響，Guinsaugon滑坡呈現(xiàn)出一定的滑坡泥石流的狀態(tài)，故采用Friction模型和Voellmy模型來進行分析(表1)，其主要是在曲線坐標系中滑體劃分為多個塊體，對各個塊體建立物理平衡方程進行求解，其中各區(qū)域塊體受到的摩擦阻力根據(jù)其Friction模型和Voellmy模型進行選取，進而求解滑體的運動距離和運動速度，具體公式見式(1)—(5)[21]，上述公式已被很多專家學(xué)者編寫成相應(yīng)的應(yīng)用程序，為高位遠程滑坡運動計算提供了支撐。

F=ΥHiBidSsinα+P-T

(1)

Friction模型：T=σ(1-γu)tanφ

(4)

式中，F(xiàn)是高位滑坡下滑力；Υ是滑體重度；H是滑塊的高度；B是滑塊的寬度；P是相鄰塊體間的切向壓力；T是滑體運動過程中受到的阻力；v是滑體速度；M是滑塊的動量通量；S是曲線位移；i是滑塊的編號；j是滑塊間的邊界次序；σ是滑體運動過程中塊體受到的豎向總應(yīng)力；γu是孔隙水壓力系數(shù)；φ是滑塊的動摩擦角。ξ是湍流項系數(shù)。

圖2 菲律賓Guinsaugon滑坡數(shù)值模型

圖3 塊體之間的受力分析

表1 運動路徑模型選擇

3.2 模型建立及參數(shù)選取

根據(jù)高位滑坡動力學(xué)特征，建立高位滑坡動力學(xué)基底流變計算模型(圖2和圖3)，結(jié)合滑體的運動軌跡，將滑體的運動路徑分為滑源區(qū)、碎屑流區(qū)和堆積區(qū)等。由于很難獲取滑體的巖土體強度參數(shù)，根據(jù)已有的大量巖質(zhì)滑坡計算案例和試錯法可知[22-24]，對于高位滑坡遠程動力學(xué)基底流變剪切計算，多采用Friction-Voellmy-Friction計算模型組合。由式(4)可知，F(xiàn)rictional模型主要受到動摩擦角φ和孔隙水壓力γu的影響，根據(jù)已有大量的巖質(zhì)滑坡案例的數(shù)值試驗，確定了模擬參數(shù)見表2。

表2 Guinsaugon滑坡數(shù)值模擬建議參數(shù)

4 計算結(jié)果分析

4.1 滑坡前后緣運動特征

計算結(jié)果顯示(圖4和圖5)，Guinsaugon滑坡運動總時長約為90s，滑體在水平距離X=1896m處時的運動速度為最大，約為52.95m/s，這是由于該處地處坡腳，滑體經(jīng)過重力勢能轉(zhuǎn)化成動能，到達該處時速度最大。

圖4 滑坡前后緣隨時間變化

圖5 滑體前后緣隨滑程變化

由圖4和圖5可見，滑體前緣和后緣的速度顯示滑體運動主要經(jīng)歷三個階段，以滑坡前緣為例：①滑坡急劇加速階段，滑體啟動后，在水平X=960～1000m范圍內(nèi)，滑體速度急劇增加，從0增長到33m/s，這主要是由于滑體急劇啟動，能量處于急劇釋放的過程；②滑體穩(wěn)定加速階段，在X=1000～3000m范圍內(nèi)，滑體的速度保持穩(wěn)定，基本維持在45～53m/s，這主要是由于受地形限制，滑體呈現(xiàn)管道狀碎屑流狀態(tài)，一直保持高速運動狀態(tài)；③滑體減速階段，在X=3000～4100m范圍內(nèi)，滑體速度逐漸降低，直至停止，這是由于滑體運動至平原地區(qū)，坡度急劇變緩，滑體動能逐漸消耗。

4.2 堆積體形態(tài)分析

滑體的堆積體平均厚度約為20～40m，如圖6所示。其中滑源區(qū)只有少量殘余碎屑流堆積，其厚度不足20m；碎屑流區(qū)域由于坡度相對較緩，運動路徑上停積了一定厚度的滑體，其厚度約為30m；在堆積區(qū)內(nèi)，堆積體的形態(tài)呈現(xiàn)“動車”頭型，平均厚度約為15m。此外，滑坡前緣運動至4100m處停止，與實際情況相近。

圖6 堆積體形態(tài)

4.3 典型特征點分析

根據(jù)對于滑坡災(zāi)害分析和已有的現(xiàn)場調(diào)查資料分析，本文選取了水平距離X=1100m為典型特征點，因為該點地處滑坡剪出口附近，能顯示滑坡急劇啟動后速度變化過程和滑源區(qū)堆積體的殘余變化情況。由圖7和圖8可知，該點的速度最大為45.5m/s，最終的堆積體厚度約為27m。

圖7 典型點速度隨時間變化

圖8 典型點堆積體厚度隨時間變化

4.4 模型對比和成災(zāi)模式

4.4.1 模型對比

高位滑坡動力學(xué)方法研究的主要手段包含數(shù)值模擬法和統(tǒng)計學(xué)方法。Scheidegger統(tǒng)計了大量的高位滑坡案例，提出了“等效視摩擦角”概念(圖9)[25]，即滑坡前緣最遠點與滑坡后緣點之間的連線的斜率，其本質(zhì)上就是滑體勢能轉(zhuǎn)化成動能(該方法也稱為“雪橇模型”)，見式(6)。

式中，ν0為滑體的速度；g為重力加速度；H為滑體運動點與滑坡后緣之間的垂直高差；f為等效視摩擦系數(shù)；L為滑坡計算點與滑坡后緣點之間的水平距離。

圖9 等效視摩擦系數(shù)概念示意

通過圖9和式(6)，可以計算Guinsaugon滑坡的滑程、垂直高差與等效視摩擦系數(shù)之間的關(guān)系。由此，計算Guinsaugon滑坡在滑坡剪出口(水平距離約為1000m)處的運動速度：

根據(jù)目擊人描述，Guinsaugon滑坡的平均速度約為27～38m/s[26]，其最大速度預(yù)估為40.5～57m/s。目擊證人的預(yù)估速度、雪橇模型和基底流變模型三者進行比較，雪橇模型計算結(jié)果也顯示滑坡運動過程中速度急劇增加，波動性高速運動和減速三個階段。此外，雪橇模型計算結(jié)果偏大，這是由于雪橇模型在計算過程中，忽略了滑體運動過程中受到的基底剪切阻力，且對滑動路徑微地貌進行了簡化。傳統(tǒng)的雪橇模型適用于應(yīng)急搶險過程中的應(yīng)急估算，而基于基底剪切阻力的計算模型可以精細化刻畫滑坡動力學(xué)運動過程，與實際速度相近。

表3 Guinsaugon滑坡的雪橇模型計算

4.4.2 成災(zāi)模式探討

Guinsaugon滑坡的母巖以角礫巖、碎屑巖等為主，巖體相對破碎，滑源區(qū)坡度較陡，約為48°，滑體呈現(xiàn)楔形體，其高度約為450m，并發(fā)育有多組節(jié)理，使得巖體呈現(xiàn)碎裂結(jié)構(gòu)。在滑坡發(fā)生前16d連續(xù)降雨的影響下，巖土體強度進一步降低，使得楔形體失穩(wěn)滑動，撞擊坡面并呈現(xiàn)碎塊狀，以碎屑流的形式繼續(xù)向坡角運動。隨著坡度變緩，滑體向四周散落，最終直至停止?；诖?，建議將“楔形體”塊體作為野外高位滑坡識別的重要標志之一，這是因為楔形體失穩(wěn)后，滑體撞擊地面呈現(xiàn)出巨大的動力沖擊力，并擾動(動力侵蝕)滑動路徑上松散的堆積物，進而加劇滑體體積效應(yīng)。這種楔形體巖質(zhì)滑坡多發(fā)育在我國西南山區(qū)線型工程兩側(cè)，嚴重威脅交通安全運行。對于楔形體巖質(zhì)滑坡的風(fēng)險防控，多采用清方+支護措施，并加強監(jiān)測。

與Guinsaugon滑坡類似，我國西南山區(qū)的高位遠程滑坡碎屑流災(zāi)害頻發(fā)，諸如2009年發(fā)生的重慶雞尾山滑坡、2017年發(fā)生的貴州納雍滑坡、2022年9月5日四川瀘定發(fā)生6.8級地震誘發(fā)的大量高位遠程滑坡碎屑流災(zāi)害。這種災(zāi)害大多數(shù)在降雨、地震或者人類工程活動誘發(fā)因素的影響，滑源區(qū)巖體強度降低，進而滑坡急劇啟動、動力撞擊坡面后破碎解體，滑坡-碎屑流在向前運動過程中動力侵蝕地表松散堆積體，最后遇到溝谷和河流等阻擋后形成滑坡壩，壩體在堰塞湖的作用下潰決形成洪水災(zāi)害[27]。對于這種高位遠程滑坡鏈式災(zāi)害應(yīng)加強地質(zhì)災(zāi)害風(fēng)險調(diào)查評價，進而減少災(zāi)害損失[28]。

本次高位滑坡災(zāi)害動力學(xué)分析主要以計算滑體運動速度和堆積體厚度等為主，建模過程主要采用條塊法，以等效流體的方式來進行計算，缺乏在建模過程中考慮巖體節(jié)理、且?guī)r土體參數(shù)主要以試錯法來取得，這些與真實情況有一定差距，這也是今后研究的重要方向。

5 結(jié) 論

1)Guinsaugon滑坡屬于典型的高位遠程滑坡災(zāi)害，將其運動路徑劃分為滑源區(qū)、碎屑流區(qū)和堆積區(qū)。滑源區(qū)巖體形狀為楔形體，碎屑流區(qū)呈現(xiàn)管道狀和堆積區(qū)呈扇形。

2)Guinsaugon滑坡的主要誘發(fā)因素為地質(zhì)構(gòu)造運動活躍、巖體結(jié)構(gòu)破碎和降雨，其運動階段主要劃分為急劇啟動-撞擊碎屑化-遠程堆積的模式，其速度運動階段主要為急劇加速-波動性維持-減緩?fù)７e的過程。

3)模型計算結(jié)果顯示，滑坡運動時長約為90s，在水平距離X=1896m處時的運動速度為最大，約為52.95m/s，堆積體的平均厚度約為15～20m，最遠運動距離約為4100m，與實際情況相近。與傳統(tǒng)的雪橇模型計算結(jié)果相比，基于基底剪切阻力計算模型可以更好的精細化刻畫滑體動力學(xué)過程。