運(yùn)用Voellmy模型方法模擬泥石流運(yùn)動(dòng)堆積動(dòng)力過(guò)程

陳 偉

(廣東省環(huán)境地質(zhì)勘查院，廣東 廣州 510000)

0 引言

泥石流災(zāi)害屬于山區(qū)的一種突發(fā)性災(zāi)害，具有流速快、攜帶固體物質(zhì)能量強(qiáng)和破壞力大等特點(diǎn)，是我國(guó)主要地質(zhì)災(zāi)害類(lèi)型之一(謝洪等，2006)。開(kāi)展泥石流危險(xiǎn)范圍的預(yù)測(cè)工作具有重要的現(xiàn)實(shí)意義。

利用新一代的地理信息系統(tǒng)(GIS)，結(jié)合適用于泥石流兩相流體分析的理論模型，綜合地質(zhì)、數(shù)值計(jì)算以及計(jì)算機(jī)科學(xué)等不同學(xué)科，對(duì)泥石流運(yùn)動(dòng)以及堆積過(guò)程進(jìn)行模擬，進(jìn)而為后期泥石流危險(xiǎn)性評(píng)價(jià)區(qū)劃及防治提供關(guān)鍵參數(shù)，是當(dāng)前泥石流理論研究以及工程防治設(shè)計(jì)的研究熱點(diǎn)和難點(diǎn)。從目前國(guó)內(nèi)外泥石流運(yùn)動(dòng)堆積過(guò)程研究現(xiàn)狀出發(fā)，探討Voellmy連續(xù)介質(zhì)模型原理以及GIS在泥石流研究中的應(yīng)用，并以典型泥石流溝為例，通過(guò)模擬泥石流運(yùn)動(dòng)、堆積的動(dòng)力過(guò)程，獲取了泥石流體的最大速度、壓力、過(guò)流量、堆積高度以及堆積面積等關(guān)鍵參數(shù)。

根據(jù)對(duì)國(guó)內(nèi)外文獻(xiàn)的分析，目前對(duì)泥石流堆積扇運(yùn)動(dòng)及堆積預(yù)測(cè)通常采用3種方法：經(jīng)驗(yàn)統(tǒng)計(jì)方法、物理模擬方法、動(dòng)力模型分析。

日本是較早也是較多研究泥石流運(yùn)動(dòng)及堆積預(yù)測(cè)的國(guó)家之一。池谷浩等1979年根據(jù)流域面積推算泥石流沖出量，根據(jù)沖出量推算泥石流堆積長(zhǎng)度和堆積寬度，率先從統(tǒng)計(jì)學(xué)角度探討了這一問(wèn)題(池谷浩，1996);石川芳治(1999)研究了用模型實(shí)驗(yàn)?zāi)M堆積區(qū)具有溝槽和隆起等微地貌時(shí)，泥石流堆積范圍應(yīng)如何修正以及泥石流體中細(xì)粒物質(zhì)對(duì)泥石流堆積擴(kuò)散的影響等。

歐美國(guó)家也很重視泥石流危險(xiǎn)范圍的研究。奧地利很早就進(jìn)行了泥石流危險(xiǎn)范圍的預(yù)測(cè)工作，并引用交通信號(hào)中的紅、黃、藍(lán)3色的特定含義，將泥石流危險(xiǎn)區(qū)分為3類(lèi)(Blaikie et al，1994);Brien等開(kāi)發(fā)了水力模型，模擬泥石流和流域內(nèi)的堆積區(qū);Ghilardi等(2001)開(kāi)發(fā)了1個(gè)數(shù)學(xué)模型，考慮了侵蝕和堆積過(guò)程，以及不同類(lèi)型的堆積物;Pasuto等(2004)采用多學(xué)科方法，從歷史、地形地貌、地質(zhì)結(jié)構(gòu)、氣象、土壤和森林管理等角度，對(duì)泥石流進(jìn)行了評(píng)價(jià)。

唐川等(1992)應(yīng)用泥石流堆積模型實(shí)驗(yàn)，將實(shí)驗(yàn)結(jié)果用多元回歸曲線(xiàn)擬合，得到泥石流堆積形態(tài)函數(shù)與泥石流參數(shù)關(guān)系表達(dá)式;吳積善等(1993)根據(jù)泥石流典型地區(qū)(主要是東川小江流域)的資料，建立了泥石流最大堆積面積、最大堆積長(zhǎng)度、最大堆積寬度和泥石流堆積幅角與流域面積、松散固體物質(zhì)儲(chǔ)量、主溝長(zhǎng)度和流域最大相對(duì)高差的相關(guān)式。

利用動(dòng)力模型進(jìn)行泥石流運(yùn)動(dòng)及堆積預(yù)測(cè)是目前研究的熱點(diǎn)，該模型通常采用數(shù)值方法，利用能量與運(yùn)動(dòng)轉(zhuǎn)換法則對(duì)泥石流運(yùn)移及堆積過(guò)程進(jìn)行模擬。動(dòng)力模型包括離散元模型、質(zhì)量集中模型以及連續(xù)單元模型。其中，連續(xù)介質(zhì)模型是建立在流體力學(xué)之上，利用物質(zhì)－運(yùn)動(dòng)－能量消散的轉(zhuǎn)換方程來(lái)描述泥石流的動(dòng)力過(guò)程。因此，連續(xù)介質(zhì)模型可以更為合理地刻畫(huà)出泥石流固體物質(zhì)的運(yùn)動(dòng)特性，從而被廣泛使用(Hungr et al，1984)。

1 連續(xù)介質(zhì)模型在泥石流中的應(yīng)用

在連續(xù)介質(zhì)模型中，泥石流流體被假設(shè)為一種非穩(wěn)定以及非均質(zhì)的流體，可以被2個(gè)主要流體參數(shù)來(lái)刻畫(huà)(Hungr et al，2009)，即流體高度H(x，y，t)，m;平均流速U(x，y，t)，m/s：

式(1)中，Ux為X方向的速度;Uy為Y方向的速度;T為平均速度的轉(zhuǎn)置矩陣。速度大小可以被定義為：

式(2)中，‖U‖表示對(duì)速度U取絕對(duì)平均值，這就確保U在向量空間中為嚴(yán)格的正速度，流體速度的方向可以用單位向量(nU)來(lái)定義：

Voellmy流變模型用下列質(zhì)量平衡公式來(lái)表達(dá)：

式(4)中，H為流體高度，m;Q x，y，()t，kg/(m2·s)為質(zhì)量源，當(dāng)Q=0時(shí)，表示沒(méi)有物質(zhì)沉積。在X和Y方向上，流體平均深度平衡方程可以分別表示為：

以及：

式(5)、(6)中，Cx以及Cy為剖面系數(shù)，gz為垂直方向的重力加速度。在Voellmy模型中，垂直方向的接觸關(guān)系可以定義為非均質(zhì)的摩爾－庫(kù)倫關(guān)系。ka/p為土壓力系數(shù)，可以用下列關(guān)系式來(lái)表示：

式(7)中，φ為泥石流流體的內(nèi)摩擦角。

通過(guò)以上各式，可以得出Voellmy流變公式：

2 實(shí)例研究

2.1 泥石流溝基本特征

選取地震災(zāi)區(qū)汶川縣某典型泥石流溝作為實(shí)例研究，該溝位于岷江右岸，流域面積5.8 km2，主溝長(zhǎng)約5.2 km，以“V”型谷為主，流域海拔范圍1 320～3 340 m之間，高差 1 020 m，主溝縱比降為200‰。該溝流域內(nèi)出露基巖主要以花崗巖、板巖及千枚巖為主。溝道兩岸多基巖出露，坡度較陡，平均坡度約40°～60°。該溝主要威脅到下游的阿壩師專(zhuān)、村莊、公路及農(nóng)田等承災(zāi)體。通過(guò)對(duì)該溝流域的調(diào)查特征分析，可將該溝可分為清水、物源、流通及堆積4個(gè)區(qū)，各分區(qū)特征見(jiàn)表1。

表1 泥石流流域分區(qū)特征表Table 1 Zonation features of the debris flow basins

通過(guò)實(shí)地調(diào)查，該溝于1958年及1979年8月13日分別暴發(fā)中—大規(guī)模泥石流。其中1979年發(fā)生的泥石流持續(xù)時(shí)間約lh，速度較快，“龍頭”高度為2.0～3.1 m，泥漿飛濺，泥位高度約3.1 m(康景文等，2011)。在“5·12”汶川特大地震的強(qiáng)大作用力之下，溝內(nèi)及兩岸松散固體物質(zhì)激增，在強(qiáng)降雨或綿雨條件下，很容易暴發(fā)嚴(yán)重的泥石流災(zāi)害。

2.2 泥石流特征值計(jì)算

根據(jù)相關(guān)計(jì)算公式，分別計(jì)算出該溝全流域的暴雨歷時(shí)、設(shè)計(jì)暴雨量、徑流深度、設(shè)計(jì)洪水總量、概化矩形歷時(shí)等參數(shù)，該溝的洪水總量、泥石流峰值流量、1次泥石流總量、1次泥石流沖出的固體物質(zhì)總量等計(jì)算結(jié)果(表2、表3、表4、表5)。

表2 暴雨洪水流量計(jì)算結(jié)果Table 2 Calculation of rainstorm and flood flow

表3 泥石流峰值流量計(jì)算結(jié)果Table 3 Calculation of peak flow of debris flow

表4 1次泥石流總量計(jì)算結(jié)果Table 4 Calculation of the total amount of one-time debris flow

表5 1次泥石流固體物質(zhì)總量計(jì)算表Table 5 Calculation of the total solid material of one-time debris flow

根據(jù)計(jì)算結(jié)果，該溝20年、50年及100年一遇的泥石流總量約4 230～9 052 m3。

2.3 泥石流運(yùn)動(dòng)及堆積數(shù)值模型建立

數(shù)字地面模型是泥石流數(shù)值模擬最基礎(chǔ)，也是最重要的步驟，這直接決定了模擬結(jié)果的精度以及正確性。通過(guò)對(duì)工作區(qū)紙質(zhì)地圖進(jìn)行矢量化后，利用ArcGIS軟件，建立TIN(Triangulated Irregular Network)最優(yōu)化不規(guī)則三角網(wǎng)模型，通過(guò)對(duì)TIN模型數(shù)字高程點(diǎn)的采集，獲取模型網(wǎng)格點(diǎn)數(shù)字高程，并進(jìn)行內(nèi)插計(jì)算，最終獲得研究區(qū)DEM。通過(guò)GIS模擬系統(tǒng)內(nèi)置分析代碼，獲取流體的流速、高度、過(guò)流量等關(guān)鍵參數(shù)，并將計(jì)算所得的結(jié)果以GIS格式數(shù)據(jù)的形式以3D效果進(jìn)行展示。

在上述連續(xù)介質(zhì)模型原理的基礎(chǔ)上，采用動(dòng)力數(shù)值模擬計(jì)算軟件，對(duì)該泥石流溝進(jìn)行實(shí)例研究。模型長(zhǎng)度=7 300 m，寬度=5 040 m，共計(jì)50 034個(gè)節(jié)點(diǎn)，49 580個(gè)單元。對(duì)該溝百年一遇的泥石流影響范圍進(jìn)行預(yù)測(cè)，1次泥石流固體物質(zhì)總量為9 052 m3，建立的三維模型如圖1所示。

圖1 泥石流數(shù)值模擬三維模型Fig.1 3D numerical simulation model of debris flow

3 結(jié)果與討論

采用常規(guī)的泥石流流量過(guò)程概化線(xiàn)和洪峰流量進(jìn)行模擬，針對(duì)泥石流的堆積范圍進(jìn)行模擬，并為后期的危險(xiǎn)性評(píng)價(jià)及區(qū)劃提供技術(shù)支撐。通過(guò)模擬計(jì)算，泥石流最大速度、最大壓力、過(guò)流量以及流體堆積高度、堆積特征如圖2、圖3、圖4、圖5、表6所示。

圖2 泥石流最大速度特征圖Fig.2 Characteristic diagram of maximum velocity of debris flow

圖3 泥石流最大壓力特征圖Fig.3 Characteristic diagram of maximum pressure of debris flow

圖4 泥石流最大過(guò)流量特征圖Fig.4 Characteristic diagram of maximum momentum of debris flow

圖5 泥石流最大流體堆積高度特征圖Fig.5 Characteristic diagram of maximum accumulation height of debris flow

表6 泥石流運(yùn)動(dòng)及堆積計(jì)算表Table 6 Calculation of the movement and accumulation of debris flow

從以上分析結(jié)果圖可以看出，該溝泥石流暴發(fā)后，在泥石流的流通區(qū)和物源區(qū)，泥石流流體的最大流速可以達(dá)到6.96 m/s，快速的泥石流流體可以產(chǎn)生很強(qiáng)的側(cè)蝕和底蝕作用，加大了泥石流流體的固體物質(zhì)攜帶能力，并大幅度地增加了固體物質(zhì)運(yùn)移距離，表現(xiàn)出強(qiáng)烈的侵蝕作用。在溝谷的中下段，特別是出山口至主河之間的溝段，泥石流流速已經(jīng)有大幅度的降低，以沉積作用為主。根據(jù)數(shù)值模擬計(jì)算結(jié)果，在1次泥石流固體物質(zhì)總量為9 052 m3的情況下，泥石流流體攜帶固體物質(zhì)沖出溝口后，堆積扇面積可高達(dá)0.147 km2，將對(duì)溝口大量居民、民房、公路等造成嚴(yán)重危害。

另外，根據(jù)模擬，泥石流的堆積扇扇頂剛好沖入岷江，對(duì)岷江主河道造成擠迫，雖未完全阻斷，但會(huì)造成河流的回水，對(duì)位于溝口所處的岷江上游區(qū)段造成潛在威脅。

4 結(jié)論

(1)泥石流危險(xiǎn)范圍預(yù)測(cè)是泥石流危險(xiǎn)性評(píng)價(jià)及區(qū)劃的一項(xiàng)重要研究?jī)?nèi)容，特別是針對(duì)地震災(zāi)區(qū)泥石流災(zāi)害防治有著重要的意義。

(2)在預(yù)測(cè)泥石流危險(xiǎn)范圍時(shí)，動(dòng)力模型分析方法有著經(jīng)驗(yàn)統(tǒng)計(jì)方法以及物理模擬等方法不具備的優(yōu)勢(shì)，可以更為合理地刻畫(huà)出泥石流固體物質(zhì)的運(yùn)動(dòng)特性，從而被廣泛使用。

(3)選取地震災(zāi)區(qū)汶川縣某典型泥石流溝作為實(shí)例研究，該溝流域面積5.8 km2，主溝長(zhǎng)約5.2 km，該溝20年、50年及100年一遇的泥石流總量約為4 230～9 052 m3。

(4)采用動(dòng)力數(shù)值模擬計(jì)算軟件，建立泥石流溝的數(shù)值模型，對(duì)該溝百年一遇的泥石流影響范圍進(jìn)行預(yù)測(cè)，計(jì)算結(jié)果顯示該溝在1次泥石流固體物質(zhì)總量為9 052 m3的情況下，最大速度為6.96 m/s，最大壓力99.96 kPa，最大過(guò)流量為5.55 m2/s，最大流體堆積高度為3.1 m，堆積面積為0.147 km2。不僅對(duì)下游居民造成嚴(yán)重威脅，還有可能阻塞岷江主河道。

(5)根據(jù)野外調(diào)查情況，并結(jié)合數(shù)值模擬研究，證明基于GIS技術(shù)平臺(tái)以及Voellmy連續(xù)介質(zhì)模型可以很好地模擬泥石流運(yùn)動(dòng)、堆積的動(dòng)力過(guò)程，在泥石流危險(xiǎn)范圍預(yù)測(cè)工作中具有較高的準(zhǔn)確性。

池谷浩.1996.云仙水無(wú)川泥砂流失量的推算方法［J］.水土保持科技情報(bào)，(4)：30－34.

康景文，黃建波，周勇，等.2011.汶川縣威州鎮(zhèn)羊嶺溝泥石流特征及成因分析［J］.施工技術(shù)，40(增刊)：180－184.

唐川，劉希林.1992.泥石流動(dòng)力堆積模擬和危險(xiǎn)范圍預(yù)測(cè)模型［M］.北京：科學(xué)出版社.

謝洪，鐘敦倫，韋方強(qiáng)，等.2006.我國(guó)山區(qū)城鎮(zhèn)泥石流災(zāi)害及其成因［J］.山地學(xué)報(bào)，24(1)：79－87.

吳積善，田連權(quán)，康志誠(chéng)，等.1993.泥石流及其綜合治理［M］.北京：科學(xué)出版社.

石川芳治.1999.地震にょる土石流の発生に係わる地形，地質(zhì)條件［J］.砂防學(xué)會(huì)誌，51(5)：35－42.

BLAIKIE P，CANNON T.1994.At Risk：Natural Hazards，Peoples Vulnerability and Disasters［M］.London，UK：Psychology Press.

GHILARDI P，NATALE L，SAVI F.2001.Modelling debris－flow propagation and deposition［J］.Physics and Chemistry of the Earth：Part C，26(9)：651 －656.

HUNGR O，MORGAN G C，KELLERHALS R.1984.Quantitative analysis of debris torrent hazards for design of remedial measures［J］.Canadian Geotechnical Journal，21(4)：663－677.

HUNGR O，MCDOUGALL S.2009.Two numerical models for landslide dynamic analysis［J］.Computers ＆ Geosciences，35(5)：978－992.

JULIEN P Y，O'BRIEN J S.1997.Slected notes on debris flow dynamics［J］.Lecture Notes in Earth Sciences，64：144 －162.

PASUTO A，SOLDATI M.2004.An integrated approach for hazard assessment and mitigation of debris lows in the Italian Dolomites［J］.Geomorphology，61(1/2)：59 － 70.