基于Massflow模擬的察隅縣N047泥石流運(yùn)動特征研究

趙思遠(yuǎn) 吳 琦 柴向陽

（1.華北水利水電大學(xué)地球科學(xué)與工程學(xué)院，河南 鄭州 450003；2.河南省遙感地質(zhì)勘察有限公司，河南 鄭州 450000）

0 引言

通常認(rèn)為降雨是誘發(fā)地質(zhì)災(zāi)害的眾多外部因素之一。雨水源源不斷地匯入溝道，洪流裹挾著溝道內(nèi)、岸坡上的松散固體物質(zhì)沿溝道向下游方向流動，從而誘發(fā)泥石流災(zāi)害。因此本研究將對20 年一遇的暴雨條件下的N047 泥石流災(zāi)害采用Massflow 軟件模擬，以確定其運(yùn)動的基本規(guī)律，并模擬預(yù)測其在50 年、100 年一遇條件下的流量和泥深，為后期防災(zāi)、減災(zāi)提供必要的參考數(shù)據(jù)。

1 泥石流的基本特征

察隅縣竹瓦根鎮(zhèn)知美村N047泥石流地理坐標(biāo)為98°08'12.10″E，28°28'18.13″N，流域面積為1.27 km2，流域內(nèi)最高海拔4 481 m，最低點為3 490 m，相對高差達(dá)到989 m。如圖1 所示，溝槽橫斷面呈“V”型，流域狹長呈帶狀，溝道較筆直，未見支溝，主溝縱坡比489‰。溝域縱向長度約2.1 km，呈對稱分布，為季節(jié)性溝道，溝道內(nèi)可見孤石。溝口處地形寬闊，相對彎曲。溝源周邊群山圍繞，呈“圈椅狀”，植被覆蓋率一般，溝道深切，侵蝕較嚴(yán)重。堆積區(qū)的扇形地完整性為95%，主軸坡降為248‰，扇面發(fā)展呈淤高趨勢。扇長約155 m，扇寬約330 m，擴(kuò)散角為95°，溝口至主河道（日東河）距離約580 m。該泥石流主要威脅對象為知美村居民點，威脅溝口兩側(cè)居民7戶40人，房屋51間，威脅財產(chǎn)約619萬元。

圖1 泥石流全貌

2 泥石流的形成條件

該泥石流具備發(fā)生、發(fā)展所需的必要條件。①溝源處山體陡峭，呈立體的“圈椅狀”結(jié)構(gòu)，即三面環(huán)繞的凹槽地形。一面出口，后緣及兩側(cè)匯水面積較大，均利于匯水。形成區(qū)、流通區(qū)溝道陡峻、狹窄、深切，呈深“V”型，溝床縱坡較陡，水流在此可顯著加速。②坡面出露主要為全新統(tǒng)殘坡積物，較為松散破碎。溝道中游兩側(cè)斜坡侵蝕嚴(yán)重，在水動力作用下易破壞，形成堆積物，進(jìn)一步堆積于坡體下部，成為泥石流的重要物質(zhì)來源。③降雨充沛，且為該泥石流的唯一水源條件。

3 基于mssflow的泥石流運(yùn)動特征

Massflow 軟件是一款由中科院成都山地所歐陽朝軍開發(fā)的具有二階精度與自適應(yīng)求解特點的地表過程動力學(xué)模擬軟件［1］。是通過對三維Naiver-Stokes 方程［2］的深度積分得出質(zhì)量與動量的控制方程，并利用有限差分法（MacCormack-TVD）結(jié)合MPICH 和OpenMP 技術(shù)［3-5］求解方程所開發(fā)的模擬軟件。

3.1 模型建立

首先，采用衛(wèi)星下載知美村N047 泥石流流域的等高線地形數(shù)據(jù)，利用ArcGis軟件的空間分析功能轉(zhuǎn)換成DEM 柵格數(shù)據(jù)。其次，使用ArcGIS 的“Arc Toolbox”工具箱中的“數(shù)據(jù)管理工具”，可以劃分為5 m×5 m 大的計算網(wǎng)格，完成地形數(shù)據(jù)的處理。然后，根據(jù)ArcGIS 軟件的“柵格轉(zhuǎn)ASCII”轉(zhuǎn)換工具，將已劃分網(wǎng)格的DEM 數(shù)據(jù)轉(zhuǎn)化成ASCII格式文檔，可為Massflow 軟件識別。最后，打開Massflow軟件后臺編譯程序Microsoft Visual Studio，在程序代碼中輸入模型的相關(guān)信息，利用程序的“重新生產(chǎn)解決方案”“開始執(zhí)行不調(diào)試”執(zhí)行程序，再根據(jù)Tecplot軟件，輸出計算模型，如圖2所示。

圖2 泥石流三維模型

3.2 參數(shù)選取

3.2.1 重度。知美村N047 泥石流屬溝谷型泥石流，結(jié)合野外調(diào)查成果，確定知美村泥石流重度γ為1.572 t/m3、泥沙修正系數(shù)φ為0.54，以及堵塞系數(shù)Dc為2.2。

3.2.2 流量。泥石流流量的計算設(shè)置在溝口及擬設(shè)工程治理部位的典型斷面上，所有斷面流量的計算皆在假設(shè)未實施工程治理的條件下進(jìn)行。

首先，根據(jù)公式計算暴雨時的最大洪峰流量。該公式適合于全面產(chǎn)流條件下的全面匯流和部分匯流兩種情況，可對橫斷面處的出水點在不同降水頻率下的匯水流量進(jìn)行計算，其計算公式為式（1）。

式中；QB為最大洪峰流量，m3/s；Ψ為洪峰徑流系數(shù)；i為最大平均暴雨強(qiáng)度；F為集水面積，km2。計算結(jié)果見表1。

表1 暴雨洪峰流量計算結(jié)果

其次，采用雨洪法計算各斷面的泥石流峰值流量。其計算公式為式（2）。

式中；Qc為泥石流峰值流量，m3/s；Fc泥沙修正系數(shù)；QP為頻率為P的暴雨洪水流量；Dc為泥石流堵塞系數(shù)。計算結(jié)果見表2。

表2 泥石流峰值流量計算結(jié)果

3.2.3 流量過程線。依據(jù)流量過程線［6］（見圖3）、野外勘察與泥石流多次反演，可知知美村N047泥石流為20年一遇流量所激發(fā)，歷經(jīng)30 min（1 800 s）。

圖3 泥石流過程線

3.2.4 啟動點與摩擦模型。啟動點的選取需具備泥石流成災(zāi)的四個條件，即物源量豐富、降雨集中、溝道地形易于匯水，以及溝道地形能為泥石流啟動提供勢能空間。根據(jù)溝道內(nèi)的監(jiān)測結(jié)合啟動點具備的條件，選取溝道中游處一點，如圖2 泥石流三維模型中所標(biāo)位置，作為模擬的啟動點。Voellmy 模型改良自庫倫模型并且適用于泥漿和泥石流災(zāi)害［7］，具有較好模擬效果，故本次模擬選取Voellmy模型。

3.3 模擬結(jié)果分析

通過對知美村N047 泥石流的反演，得出該泥石流在20 年一遇狀況下60 s、100 s、600 s、900 s、1 500 s、1 800 s的泥深、流速，以及各時刻最高泥深與最大流速分布，如圖4、圖5、圖6所示。

圖4 不同時刻的泥深分布情況

圖5 不同時刻的流速分布情況

各淤積點航拍照片如圖7所示。分析模擬圖可知。

圖7 各淤積點照片

①0～60 s 時，泥石流剛啟動，由于動能較低，流速緩慢，流體由于重力作用，沿著溝道向下運(yùn)動，并不斷侵蝕兩側(cè)溝岸的松散物，此時溝道內(nèi)開始堆積，部分位于溝岸松散物下的溝道堆積較深，出現(xiàn)淤積（淤積點1）；60～100 s 時，由于后續(xù)流體的蠕動，加之溝道狹窄，流體增速明顯，加速通過溝道中游，逐漸開始加速堆積；100～600 s 時，流體持續(xù)加速，后續(xù)流體沖破淤積點1，裹挾溝道內(nèi)的松散物抵達(dá)淤積點2，該點位于溝道寬窄變化的交點且東側(cè)岸坡侵蝕嚴(yán)重，松散物激增，形成淤積，堵塞溝道；600～900 s 時，泥石流持續(xù)加速至溝口轉(zhuǎn)彎處，由于此處相對寬闊加之流體的彎道效應(yīng)，開始出現(xiàn)淤積；900～1 500 s 時，流體流至堆積區(qū)，溝道顯著變寬，溝口處流速達(dá)到最大值，并出現(xiàn)漫流，泥深增大，在轉(zhuǎn)彎處形成淤積（淤積點3），堵塞溝道；1 500～1 800 s 時，流體逐漸減速，向溝口兩側(cè)漫流堆積，泥深不斷增加，直到流體停歇。

②導(dǎo)致泥石流運(yùn)動變化的主要因素為溝道地形的不斷變換，溝道狹窄深切，地形陡峭，流速增大，溝道寬闊平緩，流速減小，泥位增加；溝道內(nèi)出現(xiàn)淤積，是由于流體不斷沖刷侵蝕溝岸，松散物大量涌入溝道造成溝道溝形發(fā)生變化，流體沿著溝道向下流動（即X、Y方向），因溝形空間的變化，流體在Z軸方向發(fā)生了漫流，即泥石流遇到障礙物時沖起爬高，此時動能部分轉(zhuǎn)化為勢能，部分能量有所消耗，導(dǎo)致流速降低，流體出現(xiàn)淤積。

③隨著時間的增加，泥深在溝道內(nèi)堆積深度越來越厚，直到泥石流停止運(yùn)動，泥深不再增加。流速隨著時間的變化，先增大后減小，最后泥石流停止運(yùn)動，流速為0 m/s；流體經(jīng)過了“蠕動-加速-減速-停歇”的過程。泥石流泥深隨著時間的推移，在溝道內(nèi)形成兩處淤積點，泥深最高可達(dá)2.9～4.1 m，在溝道寬緩處和堆積區(qū)漫流堆積，泥深較為平均，深度為0.016～0.54 m，局部最深可達(dá)1.5 m。

3.4 模擬預(yù)測分析

根據(jù)不同降雨頻率，預(yù)測知美村泥石流在50年一遇、100 年一遇的條件下，泥石流的泥深、流速如圖8所示。

圖8 50年、100年一遇泥深與流速分布

通過模擬20 年一遇、50 年一遇、100 年一遇的泥石流流態(tài)，得出1 800 s 時20 年一遇，泥石流最大流速5.134 m/s、最大泥深4.039 m；50 年一遇最大流速5.687 m/s、最大泥深4.950 m；100 年一遇最大流速6.476 m/s、最大泥深5.352 m；由此分析可得。

①在不同頻率（P=5%、P=2%、P=1%）的降雨條件下，相同點泥石流的流速、泥深及流體在堆積區(qū)的范圍均不相同，且這三者的大小均隨流量的增加而增大。

②在不同頻率（P=5%、P=2%、P=1%）的降雨條件下，溝道中的三個淤積點仍有不同程度的淤積，堵塞溝道，說明溝道兩側(cè)岸坡物源豐富。

3.5 模擬結(jié)果驗證

方法一；現(xiàn)場測量泥石流沖出物2.1×103m3，通過ArcGis 對模擬20 年一遇泥石流堆積區(qū)的堆積面積及泥深測量，得出模擬的泥石流沖出方量為1.8×103m3，誤差率-14%（即擬合度為86%）。

方法二；根據(jù)精度計算公式式（3），計算該模型模擬的精度［8］。

式中；Accuracy為精度；S0為實際泥石流堆積范圍和模擬結(jié)果的重疊區(qū)域面積；Sa為實際堆積面積；Sn為模擬堆積面積。

如表3所示，知美村N047泥石流模擬的精度達(dá)84.5%，精度較高。

表3 知美村N047泥石流模擬結(jié)果驗證

通過以上兩種不同方式的驗證，Massflow 模擬的知美村N047 泥石流，精度（擬合度）為84.5%，即模擬結(jié)果較為可靠。

4 結(jié)論

①本研究采用Massflow 數(shù)值仿真模擬軟件對察隅縣知美村N047 泥石流進(jìn)行了不同降雨狀況下的模擬（以20 年一遇狀況為主），得出結(jié)果精度（擬合度）可達(dá)84.5%，即模擬結(jié)果較為可靠。

②該泥石流從開始到結(jié)束流體歷經(jīng)了“蠕動（0～60 s）-加速（60～1 500 s）-減速（1 500～1 800 s）-停歇（＞1 800 s）”的過程。而導(dǎo)致泥石流運(yùn)動變化的原因是溝道地形的變化及堆積物涌入溝道而造成的溝形變化。

③該泥石流在不同頻率的降雨條件下，相同點位流速與泥深并不相同，流體在堆積區(qū)的范圍也不相同，它們均與流速呈正相關(guān)。