基于Voellmy模型的哀牢山地區(qū)滑坡泥石流災害運動過程數(shù)值模擬

陳亮+宗士昌+李祥龍+陳紅旗

摘要：我國西南山區(qū)滑坡泥石流災害常發(fā)，滑坡泥石流災害體運移堆積的范圍是開展區(qū)域地質(zhì)災害危險性評價的重要指標之一。以云南省哀牢山地區(qū)水塘鎮(zhèn)芭蕉樹滑坡泥石流災害為研究對象，利用Voellmy流體運動模型進行災害體運移堆積全過程數(shù)值模擬反分析，獲取數(shù)值模擬模型參數(shù)。利用反分析得到的模型參數(shù)建立與芭蕉樹滑坡泥石流所處地質(zhì)環(huán)境條件相似的戛灑鎮(zhèn)南恩小學滑坡數(shù)值模型，進行滑坡泥石流災害運移堆積過程模擬分析，模擬結果與野外地質(zhì)災害危險性評估結論符合程度較高。

關鍵詞：滑坡；泥石流；Voellmy模型；數(shù)值模擬

中圖分類號：P642 文獻標識碼：A 文章編號：1672-1683（2017）03-0113-07

Abstract：Landslide and debris flow are two common geological hazard types in the South-western mountainous area in China.One of the key factors for hazard assessment of landslides and debris flows is their post-initiation mobility.First，the runout process of Bajiaoshu landslide/debris flow was simulated with the Voellmy frictional model，and parameters were obtained.Then the simulation parameters were used to simulate the runout process of Nan′en School landslide for runout prediction and hazard assessment.The results proved to be consistent with field investigation results and earlier hazard assessment conclusions.

Key words：landslide；debris flow；Voellmy model；numerical simulation

滑坡泥石流災害是我國西南山區(qū)常見的地質(zhì)災害類型，其高速遠距離運移堆積能力使其具有較強的破壞力和很大的影響范圍。對其危險性的評價，除了需要計算評估滑坡體穩(wěn)定性和泥石流發(fā)生頻率之外，還需要評估滑坡泥石流發(fā)生后災害體的運動速度和堆積范圍。

21世紀以來，國內(nèi)外對滑坡-碎屑流運動過程的研究主要體現(xiàn)在經(jīng)驗統(tǒng)計[1-4]、力學模型計算、物理模型試驗及數(shù)值模擬[5-7]4個方面。A.Heim[8]首先提出滑坡架空等效摩擦因數(shù)ef 概念；A.E.Scheidegger等[9-10]在此基礎上提出一系列基于體積-等效摩擦因數(shù)的滑距預測公式；王念秦等[11-12]也在我國大型滑坡研究的基礎上提出了滑距預測經(jīng)驗公式。力學計算模型中，滑坡塊體模型應用最為廣泛，包括單一質(zhì)點塊體模型[13-14]、獨立單元模型[15]以及可變塊體運動二維模型[16]等；物理模型試驗主要通過室內(nèi)模型模擬碎屑流物質(zhì)的滑動和擴散過程，分析碎屑流鏟刮滑面、裹挾堆積體等機制，并為滑距計算提供試驗參數(shù)[17-18]。由于滑坡在運動過程中屬于非牛頓流體，在流動時不同高度的流體速度不同，因此層間會有摩擦所產(chǎn)生的剪應力[19-20]。漿體的流變特性是指其剪切變形時的剪應力與剪切率的關系。在描述這種非牛頓流變特性時，通常選取摩擦模型、Voellmy模型和Bingham模型[19-20]。

本文首先以云南省哀牢山地區(qū)水塘鎮(zhèn)芭蕉樹滑坡泥石流災害為研究對象，利用Voellmy模型進行災害體運移堆積全過程數(shù)值模擬反分析，測試該模型模擬泥石流運移堆積過程的適用性，并通過試錯方法獲取數(shù)值模擬模型參數(shù)；利用反分析得到的模型參數(shù)建立與芭蕉樹滑坡泥石流所處地質(zhì)環(huán)境條件相似的戛灑鎮(zhèn)南恩小學滑坡數(shù)值模型，進行滑坡泥石流災害運移堆積過程模擬分析，并與現(xiàn)場調(diào)查結果和危險性評估結論進行對比，進一步確定了Voellmy模型的適用范圍和模擬參數(shù)選取方法。

1 Voellmy流體運動學模型原理

在連續(xù)介質(zhì)模型中，泥石流流體被假設為一種非穩(wěn)定以及非均質(zhì)的流體，可以被2個主要流體參數(shù)來刻畫，即流體高度H（x，y，t）（m）；平均流速U（x，y，t）（m/s）

2 水塘鎮(zhèn)芭蕉樹滑坡泥石流運動過程模擬反分析

2.1 滑坡基本特征

芭蕉樹滑坡泥石流災害2002年發(fā)生于云南省哀牢山地區(qū)水塘鎮(zhèn)鄉(xiāng)邦邁村原芭蕉樹村民小組上方坡體。滑坡位于哀牢山中部，為構造剝蝕高中山地貌?；掳l(fā)生在哀牢山群下亞組（Ptaa）地層之中，滑體巖性為黑云斜長片麻巖，片理層較薄，最薄片理厚度僅為2 cm，片理面可見白色的拉伸線理和眼球狀構造，片理產(chǎn)狀在50°～65°∠36°～45°之間，巖層順傾坡外。片麻巖底部為一層厚度不足5 cm的白云母片巖夾層，抗剪強度和摩擦系數(shù)均較小，該夾層成為順巖體片理面的易滑面，成為坡體失穩(wěn)的優(yōu)勢面。作為紅河深大斷裂的一支，哀牢山斷裂橫穿滑坡體?；略磪^(qū)位于哀牢山斷裂下盤，滑坡堆積區(qū)則覆蓋于哀牢山斷裂上盤[21]。

該滑坡的平面形態(tài)呈“長橢圓形”（見圖1），滑坡最大縱向長度為540 m，滑坡后緣陡坎至前緣剪出口處最大長度為126 m，地形呈圈椅狀，滑體分為兩部分，主堆積區(qū)平均寬度120 m，其前緣左側為平均長度50 m，寬26 m的次級堆積區(qū)，整個滑體呈長橢圓形，中部較寬，兩頭較窄，滑坡面積為5.36萬m2?；麦w后緣高程1 170 m，前緣高程1 002 m，整體落差168 m?；缕骄L度157 m，平均寬度120 m，整體平均厚度6～8 m，體積約13.2萬m3。通過對滑坡結構特征的分析，可將滑坡分為5個區(qū)域（見圖2）：滑源區(qū)、滑體飛行區(qū)、碎屑流堆積區(qū)、彈射體堆積區(qū)以及邊緣氣浪影響帶。

圖1 芭蕉樹滑坡地質(zhì)平面

Fig.1 Geological plan of Bajiaoshu landslide

圖2 芭蕉樹滑坡地質(zhì)剖面

Fig.2 Bajiaoshu Landslide profile

通過對滑坡附近居民的走訪及查問，滑坡所處山體斜坡近年來一直比較穩(wěn)定，并未發(fā)現(xiàn)貫通裂縫及局部剝、墜落等斜坡變形跡象。但2002年8月，滑坡所處區(qū)域發(fā)生了一次罕見的極端降雨過程，從8月13日開始，滑坡所處的山坡及坡腳出現(xiàn)洪水，洪水深度接近成年人膝蓋，山坡表面的殘坡積物等松散堆積層被洪水沖刷帶走形成較小規(guī)模的坡面泥石流。13日上午8時30分，該滑坡開始啟動，伴隨著轟鳴的爆裂聲，滑坡體整體從前緣剪出口沖出，向前下方飛行了50 m左右后墜地解體，滑坡主體分解成為碎屑流繼續(xù)向前下方滑動364 m后停止運動；滑坡左側一小部分滑體墜地后發(fā)生彈射，再次飛行170 m后墜地分解；滑坡主體分解成為碎屑流運動中，其前緣和側緣產(chǎn)生近30秒的強大氣浪，這股氣浪將附近的芭蕉樹小組摧毀。

2.2 運動過程數(shù)值模擬反分析

2.2.1 模型建立與模擬參數(shù)取值

數(shù)字高程模型是滑坡泥石流數(shù)值模擬最基礎、最重要的步驟，這直接決定了模擬結果的精度以及正確性。通過對工作區(qū)地圖進行矢量化后，利用ArcGIS軟件，建立TIN（Triangulated Irregular Network）最優(yōu)化不規(guī)則三角網(wǎng)模型，通過對TIN模型數(shù)字高程點的采集，獲取模型網(wǎng)格點數(shù)字高程，并進行內(nèi)插計算，最終獲得研究區(qū)DEM。通過GIS模擬系統(tǒng)內(nèi)置分析代碼，可以得到流體的流速高度過流量等關鍵參數(shù)，并將計算所得結果以GIS格式數(shù)據(jù)的形式進行三維效果展示。在前述連續(xù)介質(zhì)模型原理基礎上，利用動力數(shù)值模擬計算軟件，對該滑坡泥石流進行實例研究。

模型長度=752 m，寬度=650 m，共計48 995個節(jié)點，48 552個單元根據(jù)現(xiàn)場地質(zhì)災害調(diào)查結果，估算滑坡初始固體物質(zhì)總量為121 469 m3，建立的三維模型見圖3，為滑坡滑動前平面位置，圖中彩色部分為滑坡體，灰色為地表。

滑坡體的勘察作業(yè)中進行了對滑體物質(zhì)的鉆探取樣和室內(nèi)土體剪切試驗，試驗結果見表1，包括土樣的快剪強度和殘余強度。由于此模擬關注于當滑坡啟動之后發(fā)生運移和堆積的過程，在模擬中將滑坡體定義為已經(jīng)產(chǎn)生滑動破壞的黏性巖土體，因此模擬反分析的參數(shù)取值區(qū)間以殘余剪切強度為依據(jù)，而滑體密度采用了土工試驗值，滑體厚度采用工程地質(zhì)剖面圖中滑體的平均厚度10 m。地表參數(shù)表示滑坡泥石流物質(zhì)的運動路徑上的地表的摩阻力，由于地表的復雜和不均勻性，無法直接用試驗獲取該參數(shù)，因此作為變量，利用迭代反分析獲得。滑坡初始速度設置為零。經(jīng)過多次反分析模擬，通過對滑坡泥石流模擬運動過程和堆積形態(tài)與實際滑坡泥石流堆積形態(tài)進行比對，最終芭蕉樹滑坡泥石流模擬參數(shù)取值見表1。

2.2.2 模擬反分析結果

數(shù)值模擬反分析結果見圖4、圖5，其中圖4為芭蕉樹滑坡泥石流運動堆積模擬過程。從圖5中可以看出，模擬的滑坡經(jīng)歷了沿滑動面下滑、沿平緩地表運移和完成堆積的過程，在整個運動過程中滑坡體產(chǎn)生了很大的變形，不再保持滑動前的形態(tài)，巖土體流體化，形成了沿一條較窄的路徑堆積的形態(tài)。

圖5為滑坡沿運移堆積方向的堆積體縱剖面圖，圖中綠色曲線為地形剖面線，紅色曲線為滑坡泥石流堆積體表面剖面線，其縱坐標讀取右側坐標軸高程值；灰色區(qū)域為滑坡堆積體剖面，其縱坐標讀取左側坐標軸堆積厚度值。將圖5結果與圖2實測堆積剖面做對比，可以看出滑坡物質(zhì)在下滑后形成的運移堆積體堆積形態(tài)與圖中實際現(xiàn)場調(diào)查得出的形態(tài)吻合度較高，可以看出數(shù)值模擬得到的堆積剖面基本吻合，說明模擬建模和參數(shù)取值是合理的。

圖4 芭蕉樹滑坡泥石流運動堆積模擬過程

Fig.4 Flow simulation process of Bajiaoshu landslide and debris flow

圖5 芭蕉樹滑坡堆積體縱剖面

Fig.5 Longitudinal profile of Bajiaoshu landslide

3 南恩小學滑坡泥石流運動過程模擬

3.1 滑坡基本特征

該滑坡體位于云南省哀牢山地區(qū)戛灑鎮(zhèn)南恩河東南側的斜坡體之上，為一由片麻巖大塊石與坡積土混合組成的崩塌堆積體滑坡，見圖6?；潞缶壈l(fā)育在老滑坡堆積體中，即老滑坡堆積體前緣過渡部位，呈半圓弧狀，覆蓋層由第四系含碎石粉質(zhì)黏土和層內(nèi)夾的崩塌片麻巖孤石組成，下伏基巖為片麻巖，無基巖出露，地面標高1 906～1 908 m，坡度約40°，地形較陡?；轮骰较驗?0°，后緣高程1 770 m，前緣高程1 604 m，高差166 m，縱向長571 m，橫向?qū)?60～227 m，面積為6.5萬m2，滑體厚度一般18～26 m，平均厚度約22 m，規(guī)模約143.0萬m3，為大型中層牽引式土層滑坡。

滑坡體的發(fā)展分為兩個階段，第一階段是20世紀90年代初S307省道的修建和1999年7月南恩小學的修建，開挖坡腳減小了坡體的阻滑段，改變了坡體內(nèi)部的應力狀態(tài)，隨著坡體內(nèi)的應力調(diào)整，堆積體蠕滑變形逐漸加大；第二階段是暴雨觸發(fā)階段，2002年8月11日-14日的暴雨是這次滑坡形成的最直接的誘因，據(jù)調(diào)查數(shù)據(jù)顯示，2002年8月11日-14日，境內(nèi)發(fā)生強降雨，8月11日至12日，降雨量為4 mm，13日升至27 mm，到14日降雨量猛升至101 mm，在連續(xù)降雨，特別是14日的強降雨情況下，大量的雨水下滲，導致滑坡區(qū)土體飽和，滑體重量加大，基覆界面處土體的抗剪強度急劇降低，滑坡從前緣開始逐漸滑動，南恩小學操場出現(xiàn)張裂縫、教室開裂、擋墻拉裂及S307省道擋墻拉裂，后緣的居民房屋墻體位移明顯。斜坡體前部的強烈變形導致老滑坡堆積體后緣出現(xiàn)下錯，并在隨后的自身調(diào)整中達到了自我平衡，目前處于基本穩(wěn)定狀態(tài)。但考慮到滑坡體在強降雨條件下可能會產(chǎn)生滑動并演化為泥石流，對下游南恩河沿岸多個村寨存在很大威脅，因此需要進行滑坡泥石流啟動后運動堆積距離和破壞范圍進行模擬預測，以進行危險性評價。

3.2 運動過程數(shù)值模擬

利用南恩小學滑坡勘察資料，建立三維地形模型和滑坡體初始位置模型見圖7。由于南恩小學滑坡所處位置地質(zhì)環(huán)境條件與芭蕉樹滑坡泥石流基本一致，因此直接應用芭蕉樹滑坡泥石流數(shù)值模擬反分析得到的模擬參數(shù)（表1），進行南恩小學滑坡運動過程模擬。

假設滑坡體已啟動，初始速度為零。數(shù)值模擬得到的滑坡泥石流運移堆積模擬過程見圖8。從圖8南恩小學滑坡運移堆積模擬過程平面圖中可以看出，滑坡運移堆積產(chǎn)生了分叉，分別沿北側和南側兩條山谷進行由西向東的運動。圖9為勘察結果與模擬結果對比圖，底圖為勘察所評估的南恩小學滑坡泥石流運動路徑，底圖西南角覆蓋層為數(shù)值模擬滑坡泥石流堆積范圍，可以看出模擬結果與評估結果符合程度很高。圖10和圖11分別為北側和南側山谷內(nèi)部滑坡泥石流堆積體剖面圖。對比模擬結果和地質(zhì)災害調(diào)查結果，可以發(fā)現(xiàn)模擬中滑坡泥石流的運移路徑存在分叉，分為南北兩條，北部堆積體又存在二級分叉。這一模擬結果體現(xiàn)出了以GIS為平臺運用Voellmy模型的優(yōu)勢，即以數(shù)字柵格為模擬單元，進行單元內(nèi)質(zhì)量守恒和整體動量守恒的數(shù)值計算，不受傳統(tǒng)有限體積法模擬中滑坡模型的網(wǎng)格剖分限制，可以實現(xiàn)滑坡體的離散性模擬。

4 結論

（1）Voellmy流體運動模型遵守動量守恒定律，可以較好得模擬滑坡泥石流啟動后的運移和堆積過程，從而得到滑坡泥石流災害的破壞范圍和破壞強度。

（2）在相似的地質(zhì)環(huán)境條件下，對同一類型的地質(zhì)災害進行破壞運動過程數(shù)值模擬所需的計算參數(shù)可以取同一組數(shù)值。

（3）Voellmy模型參數(shù)的取值會顯著影響模擬結果，由于模型參數(shù)不易通過實驗獲取，因此合理的方法是選取所需模擬的滑坡泥石流所在區(qū)域內(nèi)已發(fā)生的滑坡泥石流災害進行數(shù)值模擬反分析，通過迭代法得到Voellmy模型參數(shù)取值，作為模擬參數(shù)。

參考文獻（References）：

[1] 〖ZK（#]Prochaska A B，Santi P M，Higgins J D，et al.Debris-flow runout predictions based on the average channel slope （ACS）.Eng Geol[J].Engineering Geology，2008，98（1）：29-40.DOI：10.1016/j.enggeo.2008.01.011

[2] Berti M，Simoni A.DFLOWZ：A free program to evaluate the area potentially inundated by a debris flow[J].Computers & Geosciences，2014，67（2）：14-23.DOI：10.1016/j.cageo.2014.02.002

[3] Gartner J E，Cannon S H，Santi P M，et al.Empirical models to predict the volumes of debris flows generated by recently burned basins in the western U.S.[J].Geomorphology，2008，96（s 3-4）：339-354.DOI：10.1016/j.geomorph.2007.02.033

[4] Dahl M P J，Mortensen L E，Jensen N H，et al.Magnitude-frequency characteristics and preparatory factors for spatial debris-slide distribution in the northern Faroe Islands[J].Geomorphology，2013，188（436）：3-11.

[5] Chen H X，Zhang S，Peng M，et al.A physically-based multi-hazard risk assessment platform for regional rainfall-induced slope failures and debris flows[J].Engineering Geology，2015.DOI：10.1016/j.enggeo.2015.12.009

[6] Calvo B，Savi F.A real-world application of Monte Carlo procedure for debris flow risk assessment[J].Computers & Geosciences，2009，35（5）：967-977.DOI：10.1016/j.cageo.2008.04.002

[7] Aronica G T，Biondi G，Brigandì G，et al.Assessment and mapping of debris-flow risk in a small catchment in eastern Sicily through integrated numerical simulations and GIS[J].Physics & Chemistry of the Earth Parts A/b/c，2012，49（211）：52-63.DOI：10.1016/j.pce.2012.04.002

[8] HEIM A.Landslide and human lives [M].Vancouver：Bitech Press，1989：57-79.

[9] SCHEIDEGGER A E.On the prediction of the reach and velocity of catastrophic landslides[J].Rock Mechanics，1973，5（4）：231-236.DOI：10.1016/0148-9062（74）91709-4

[10] COROMINAS J.The angle of reach as a obility index for small and large landslide[J].Canadian Geotechnical Journal，1996，33（2）：260-271.

[11] 王念秦，張倬元.一種典型黃土滑坡的滑距預測方法[J].西北大學學報，2003，33（1）：111-114.（WANG Nian-qin ，ZHANG Zhuo-yuan.A forecasting method of sliding distance on typical loess landslides[J].Journal of Northw est University，2003，33（1） ：111-114.（in Chinese）） DOI：10.3321/j.issn：1000-274X.2003.01.030

[12] 王思敬，王效寧.大型高速滑坡的能量分析及其災害預測[C].滑坡論文集.成都：四川科技出版社，1989：117-124.（WANG Si-jing，WANG Xiao-ning.Energy analysis and disaster prediction of large - scale high - speed landslide[C].Landslide Proceedings.Chendu：Sichuan Science and Technology Press，1989：117-124.（in Chinese））

[13] FANNIN R J，WISE M P.A method for calculation of debris flow travel distance[C]// Proceedings of the 48th Canadian Geotechnical Conference.[S.l.]：[s.n.]，1996：643-650.

[14] MCLELLAN P J，KAISER P K.Application of two-parameter model to rock avalanches of the Mackenzie mountains[C]// Proceedings of the 4th International Symposium on Landslide.[S.l.]：[s.n.]，1984：135-140.

[15] CUNDALL P A，STRACK O D L.A discrete numerical model for granular assemblies[J].Geotechnique，1979，29（1）：47-65.

[16] 劉忠玉，馬崇武，苗天德，等.高速滑坡遠程預測的塊體運動模型[J].巖石力學與工程學報，2000，19（6）：742-746.（LIU Zhong-yu，MA Chong-wu，MIAO Tian-de，et al.Kinematic block model of long run-out prediction for high-speed landslides[J].Chinese Journal of Rock Mechanics and Engineering，2000，19（6）：742-746.（in Chinese）） DOI：10.3321/j.issn：1000-6915.2000.06.012

[17] EGASHIRA S，HONDAB N，ITOHC T.Experimental study on the entrainment of bed material into debris flow[C]// Physics and Chemistry of the Earth，Part C：Solar，Terrestrial and Planetary Science.[S.l.]：[s.n.]，2001：645-650.DOI：10.1016/S1464-1917（01）00062-9

[18] MANZELLA I，LABIOUSE V.Qualitative analysis of rock avalanches propagation by means of physical modelling of non-constrained gravel flows[J].Rock Mechanics and Rock Engineering，2008，41（1）：133-151.DOI：10.1007/s00603-007-0134-y

[19] 馬宗源，廖紅建，張駿.Bingham型黏性泥石流流體的三維數(shù)值模擬[J].西安交通大學學報，2008，42（9）：1146-1150.（MA Zong-yuan，LIAO Hong-jian，ZHANG Jun.Three dimensional numerical simulation of bingham viscous debris flow fluid[J].Journol of Xi′an jiaotong university，2008，42（9）：1146-1150.（ in Chinese）） DOI：10.3321/j.issn：0253-987X.2008.09.018

[20] 齊超，邢愛國，殷躍平，等.東河口高速遠程滑坡-碎屑流全程動力特性模擬[J].工程地質(zhì)學報，2012，20（3）：334-339.（QI Chao，XING Ai-guo，YIN Yue-ping，et al.Numerical simulation of dynamic behavior of Donghekou rockslide-debris avalanche[J].Journal of Engineering Geology，2012，20（3）：334-339.（ in Chinese）） DOI：10.3969/ j.issn.1004-9665.2012.03.005

[21] 中國地質(zhì)環(huán)境監(jiān)測院.云南新平縣地質(zhì)災害詳細調(diào)查成果報告[R].2007，106-111.（ China Institute of Geo-Environment Monitoring.Report on the findings of geological hazard in Xinping county，Yunnan province[R].2007，106-111.（in Chinese））