黃土溝壑區(qū)泥石流易發(fā)性分析
——以蘭州老狼溝為例

侯云龍，胡帥軍，彭啟園，許俊凱，吳秀剛

（1.甘肅工程地質(zhì)研究院，蘭州 730030；2.蘭州大學(xué) 土木工程與力學(xué)學(xué)院，蘭州 730000）

隨著西部地區(qū)社會(huì)經(jīng)濟(jì)的快速發(fā)展，大量農(nóng)業(yè)、水利等工程活動(dòng)在黃土高原地區(qū)快速開展。工程活動(dòng)的開展對(duì)地形地貌、地表植被覆蓋情況有著不同程度的破壞，造成水土流失、土地退化，甚者誘發(fā)滑坡、崩塌、泥石流等地質(zhì)災(zāi)害。誘發(fā)的地質(zhì)災(zāi)害對(duì)經(jīng)濟(jì)社會(huì)、人類財(cái)產(chǎn)安全帶來巨大威脅和損失。泥石流易發(fā)性是泥石流災(zāi)害危險(xiǎn)性評(píng)價(jià)工作中的有效指標(biāo)，對(duì)泥石流災(zāi)害預(yù)測(cè)和設(shè)計(jì)治理措施有著重要的指導(dǎo)意義。

目前有關(guān)泥石流易發(fā)性評(píng)價(jià)的研究可以歸納為定性評(píng)價(jià)、半定量評(píng)價(jià)和定量評(píng)價(jià)三類方法。最早采用的是定性評(píng)價(jià)方法，此方法主要用于單條泥石流溝道的分析，是從實(shí)踐中發(fā)展起來的[1-2]。隨著研究的深入化，模糊數(shù)學(xué)法、層次分析法、主成分分析法、二元統(tǒng)計(jì)[3-6]等大量數(shù)學(xué)模型被引入泥石流易發(fā)性評(píng)價(jià)工作中。隨著各種數(shù)學(xué)模型的引入，定量、半定量的泥石流易發(fā)性評(píng)價(jià)方法得到空前發(fā)展。近年來，國(guó)際上各種用于易發(fā)性評(píng)價(jià)的泥石流模擬方法也被廣泛應(yīng)用[7-8]，其中經(jīng)驗(yàn)?zāi)Ｐ途哂袩o需大量數(shù)據(jù)支持、推廣性強(qiáng)的特點(diǎn)，被廣泛應(yīng)用，如Flow-R[9]。經(jīng)驗(yàn)?zāi)Ｐ拖噍^確定性模型對(duì)基礎(chǔ)數(shù)據(jù)的數(shù)量和質(zhì)量要求較低，可以用于缺乏大量復(fù)雜數(shù)據(jù)的研究區(qū)。與統(tǒng)計(jì)模型相比，由于經(jīng)驗(yàn)?zāi)Ｐ椭械臄?shù)據(jù)與物理特性具有一定的相關(guān)性，因此它們具有相當(dāng)好的可復(fù)制性。

首先介紹Flow-R 經(jīng)驗(yàn)?zāi)Ｐ偷脑砑胺椒?，選擇甘肅省蘭州市老狼溝流域作為研究對(duì)象，在研究區(qū)的數(shù)字高程模型（Digital Elevation Model，DEM）數(shù)據(jù)基礎(chǔ)上進(jìn)行泥石流易發(fā)性評(píng)價(jià)，進(jìn)一步分析不同DEM 分辨率、不同臨界坡度閾值、不同運(yùn)動(dòng)擴(kuò)散方式對(duì)泥石流易發(fā)性概率的影響。以期為后續(xù)的危險(xiǎn)性、風(fēng)險(xiǎn)評(píng)價(jià)和治理工作帶來一定的指導(dǎo)。同時(shí)對(duì)其他相似地質(zhì)條件下的泥石流易發(fā)性評(píng)價(jià)提供一定的借鑒意義。

1 研究區(qū)概況

老狼溝位于蘭州市城區(qū)南側(cè)，皋蘭山以東，城關(guān)區(qū)皋蘭山鄉(xiāng)民族村以北。研究區(qū)距離蘭州火車站950 m。老狼溝具有典型的黃土溝壑區(qū)地貌，溝谷為南北走向，南北方向長(zhǎng)約2.3 km，東西寬約0.9 km，面積約0.84 km2。溝的西側(cè)發(fā)育兩條支溝，均呈北東—南西走向，長(zhǎng)度分別為0.5 km 和0.4 km，如圖1所示。研究區(qū)整體呈南高北低的特征，高程范圍1 720～2 090 m，巨大的高差使得研究區(qū)的平均坡度達(dá)到了31.7°，甚至在溝口處達(dá)到了58°，為泥石流物源的擴(kuò)散運(yùn)動(dòng)提供了良好的溝道條件。

現(xiàn)場(chǎng)調(diào)查及相關(guān)資料顯示老狼溝地層較為簡(jiǎn)單。黃土層分布廣而深，構(gòu)成了研究區(qū)地層的主體，下部的第三系泥巖和砂礫巖僅在溝口處有出露，研究區(qū)溝口出露的泥巖產(chǎn)狀為傾向175°，傾角21°[10]。研究區(qū)地處盆地東部南側(cè)的中新生代中部斷陷帶內(nèi)。強(qiáng)烈的剪切擠壓作用使斷陷帶的更新統(tǒng)先褶皺后斷裂，導(dǎo)致該地區(qū)成為東部斷塊并且一直抬升[11]。同時(shí)，黃土地區(qū)本身生態(tài)脆弱，研究區(qū)植被覆蓋較少，水土流失嚴(yán)重。地震、暴雨的發(fā)生經(jīng)常為泥石流的形成帶來大量松散堆積物源[10]。溝口附近有大量農(nóng)家樂、公園等基礎(chǔ)設(shè)施，以及住宅建筑，受到泥石流災(zāi)害的嚴(yán)重威脅。

2 研究方法

Flow-R 是用于泥石流數(shù)值模型的開源軟件。該軟件的特點(diǎn)：對(duì)數(shù)據(jù)集的要求較低，可通過圖形用戶界面定制數(shù)據(jù)輸入、算法和參數(shù)[9-12]。該模型允許在用戶給定條件下自動(dòng)劃分源區(qū)，并基于各種傳播算法和簡(jiǎn)單的摩擦定律對(duì)傳播范圍進(jìn)行評(píng)估。DEM 是所有數(shù)據(jù)中貫穿整個(gè)模擬過程的基礎(chǔ)數(shù)據(jù)資料，其數(shù)據(jù)質(zhì)量對(duì)結(jié)果的準(zhǔn)確性至關(guān)重要。該模型的模擬過程可分為兩步：首先根據(jù)形態(tài)學(xué)和用戶自定義的準(zhǔn)則對(duì)源區(qū)域進(jìn)行識(shí)別，其次根據(jù)摩擦規(guī)律和流動(dòng)方向算法對(duì)這些源區(qū)進(jìn)行泥石流的運(yùn)動(dòng)模擬，如圖2所示[9]。由于大尺度下常受侵蝕和沉積的影響，對(duì)于準(zhǔn)確判定泥石流的質(zhì)量和體積存在較大難度，因而，該模型的計(jì)算過程中沒有考慮泥石流體積和質(zhì)量這一參數(shù)。模型模擬為區(qū)域范圍的易發(fā)性初步評(píng)估提供了一個(gè)實(shí)質(zhì)性的基礎(chǔ)。該模型還用于對(duì)巖崩、雪崩和洪水等其他自然災(zāi)害的易發(fā)性評(píng)估。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

2.1 模型原理

2.1.1 源區(qū)識(shí)別

泥石流發(fā)生的關(guān)鍵因素是松散物源（產(chǎn)砂量）、水源（入水量）和坡度[13]。松散物源和坡度是主要控制因素，降水帶來的水源是觸發(fā)因素[14]。在區(qū)域尺度下，主要選用以下四個(gè)基礎(chǔ)評(píng)價(jià)因子反映上述泥石流發(fā)生的三大關(guān)鍵因素：巖性分布圖、平面曲率、坡度和匯流累積量。因?yàn)閹r性和地貌形態(tài)分別影響碎屑物的形成和沉積，因而可以用巖性和曲率這兩種評(píng)價(jià)因子反映松散物源。匯流累積量可以間接反映降雨和融雪帶來的入滲量。第三個(gè)因素是地形坡度，可以一定程度反映土體和巖屑的抗剪強(qiáng)度。此外，還可以添加其他評(píng)價(jià)因子，比如地質(zhì)圖、土地利用圖來更精確地識(shí)別源區(qū)。

圖2 Flow-R模擬原理示意圖Fig.2 Schematic diagram of Flow-R simulation

不同的評(píng)價(jià)因子可以代表不同類型的空間信息，并使用用戶定義的參數(shù)進(jìn)行處理[12]。因此,基于索引的方法，每個(gè)評(píng)價(jià)因子的柵格單元分為：（1）有利的，有可能啟動(dòng)；（2）排除的，不能啟動(dòng)；（3）忽略的，啟動(dòng)判別時(shí)可以忽略這個(gè)參數(shù)（圖3）。各評(píng)價(jià)因子按照如下規(guī)則來確定源區(qū)：如果某柵格單元至少一次被選為有利的，但從未被認(rèn)定為排除的，則該單元可被視為源區(qū)單元（圖3）[9]。

圖3 源區(qū)識(shí)別方法示意圖Fig.3 Schematic diagram of identification method for source area

2.1.2 擴(kuò)散運(yùn)動(dòng)模擬

擴(kuò)散區(qū)域的評(píng)估一方面基于擴(kuò)散的流動(dòng)方向權(quán)重算法計(jì)算泥石流的運(yùn)動(dòng)路徑，另一方面基于基本的能量平衡計(jì)算泥石流的最大運(yùn)動(dòng)距離。兩者都可以根據(jù)需要選擇組合使用。從確定每個(gè)源區(qū)的潛在泥石流路徑、泥石流方向（以概率為單位）、計(jì)算泥石流的流動(dòng)距離入手，對(duì)泥石流的流動(dòng)和擴(kuò)散進(jìn)行評(píng)價(jià)。

（1）擴(kuò)散評(píng)價(jià)方程

擴(kuò)散方向算法定義了從物源柵格開始，到周圍不同方向柵格的傳播概率[15]。有多種算法可用：

Holmgren 算法[16]在多流向算法中增加了一個(gè)指數(shù)x參數(shù)，控制擴(kuò)散度：

式中：i，j為流動(dòng)方向?yàn)閕方向的易發(fā)性占比，tanβi為中心柵格到i方向柵格的坡度斜率，x為變量指數(shù)。當(dāng)x增大時(shí)，擴(kuò)散度減小。在野外和實(shí)驗(yàn)室測(cè)量的基礎(chǔ)上，Claessens 等[17]提出泥石流的指數(shù)值為4。

持久性函數(shù)可以很好地反映物源運(yùn)動(dòng)過程的慣性行為，以及因運(yùn)動(dòng)方向相對(duì)于前一個(gè)方向的變化對(duì)總體流動(dòng)方向改變的權(quán)重[18]。

式中：是慣性在i方向流動(dòng)的易發(fā)性比例，a(i)為前一方向與從中央單元到下一單元方向中間的夾角。

持久性函數(shù)中的wa(i)有三種取值方法（表1）：第一種為正切，第二種為余弦，第三種基于Gamma’s 模型[18]。在每個(gè)慣性分布中，與流動(dòng)方向相反的慣性權(quán)重值為零（W180=0），以避免最終的向后傳播，從而節(jié)省計(jì)算時(shí)間。

表1 持久性函數(shù)在運(yùn)動(dòng)方向評(píng)估中的權(quán)重Tab.1 Weights of persistence equation in motion direction evaluation

根據(jù)式（3）將擴(kuò)散方向算法得到的值與持久性函數(shù)的權(quán)重結(jié)合起來，即為單元的總體易發(fā)性概率：

式中：i，j為流動(dòng)方向；pi為方向i的易發(fā)性概率；為擴(kuò)散方向算法計(jì)算的流動(dòng)易發(fā)性比例；為持久性計(jì)算的流動(dòng)易發(fā)性比例；p0為中心單元先前確定的易發(fā)性概率值。

（2）運(yùn)動(dòng)距離方程

運(yùn)動(dòng)距離是利用能量平衡、摩擦損失函數(shù)和最大速度閾值來確定的?；诮y(tǒng)一的能量平衡方程和簡(jiǎn)單的摩擦定律[式（4）]控制泥石流的運(yùn)動(dòng)距離，以及向兩側(cè)的橫向擴(kuò)散范圍。

式（4）中的摩擦損失可以用兩種算法進(jìn)行評(píng)估：Perla’s 模型[19]的兩參數(shù)摩擦模型和SFLM 模型（The simplified friction-limited model）。根據(jù)參數(shù)的選擇，這兩種方法都可以產(chǎn)生相似的傳播區(qū)域[20]。

2.2 評(píng)價(jià)因子及控制參數(shù)

在ArcGIS 中數(shù)字化處理老狼溝流域范圍的1∶10 000地形圖，計(jì)算出分辨率為10 m的DEM。研究區(qū)整體呈南高北低的特征，高程范圍1 720～2 090 m，巨大的高差提供了很好的匯水條件。通過空間分析計(jì)算出流域范圍的坡向、坡度等地形信息。對(duì)DEM 進(jìn)行平滑處理，填充匯流槽，得到了累積匯流量和平面曲率。

研究區(qū)的平均坡度達(dá)到了31.7°，甚至在溝口處達(dá)到了58°，坡度是決定泥石流觸發(fā)的主要因素。大多數(shù)泥石流發(fā)生在坡度大于15°的地區(qū)[13，21]。因此，認(rèn)為所有潛在物源的位置在坡度大于15°的區(qū)域。大量調(diào)查發(fā)現(xiàn)，在平面曲率為凹的地方更易發(fā)現(xiàn)泥石流源區(qū)[20]。根據(jù)Bavman 等人[15]在瑞士德沃州，Delmonaco與Fischer等人[16，22]在挪威的經(jīng)驗(yàn)，泥石流源區(qū)的平面曲率閾值被設(shè)定為-2/100 m-1。研究區(qū)內(nèi)主要以第四系松散堆積為主，因此認(rèn)為在巖性方面全區(qū)均具有成為潛在源區(qū)的可能。

3 結(jié)果與討論

利用Flow-R 模型在研究區(qū)范圍內(nèi)模擬出天然狀態(tài)僅受降雨影響的洪水災(zāi)害、受崩塌堆積體的泥石流災(zāi)害、受已有滑坡影響的泥石流災(zāi)害以及受崩塌滑坡災(zāi)害共同影響的泥石流災(zāi)害易發(fā)性評(píng)價(jià)結(jié)果。

圖4（a）為不考慮流域內(nèi)崩塌滑坡提供的松散物源，對(duì)老狼溝流域中的洪水現(xiàn)象進(jìn)行模擬分析。識(shí)別出源區(qū)主要分布在主溝溝頭處各支溝中。流域范圍內(nèi)的降雨經(jīng)過坡面徑流匯聚至各溝道內(nèi)，順著支溝流入主溝，一路向北流至溝口，在溝口平坦地形條件下擴(kuò)散形成扇狀。圖4（b）為受崩塌影響的泥石流易發(fā)性概率圖，崩塌形成的碎屑物質(zhì)對(duì)局部地區(qū)的地形產(chǎn)生一定程度的改變，這些碎屑物質(zhì)可以視為很好的泥石流災(zāi)害物源，受溝底溪流的淘蝕作用，溝頭的崩塌堆積物作為物源發(fā)生泥石流，流域內(nèi)的泥石流易發(fā)區(qū)向溝頭延伸，但仍然集中在溝道內(nèi)。圖4（c）為受滑坡影響的泥石流易發(fā)性概率圖，滑坡的發(fā)生常伴有大量裂縫的形成，同時(shí)對(duì)坡面的局部坡度、匯水條件有一定的改變作用。此外，溝道兩側(cè)的滑坡堆積體可為泥石流提供大量物源，不難發(fā)現(xiàn)溝道兩側(cè)坡面的松散物源發(fā)生泥石流，研究區(qū)內(nèi)的泥石流運(yùn)動(dòng)擴(kuò)散范圍有了明顯的擴(kuò)大?；潞捅浪鸀?zāi)害為泥石流的形成提供了大量物源和動(dòng)能。該工況下的泥石流易發(fā)區(qū)范圍兼顧前兩種工況的特點(diǎn)，相較僅考慮受滑坡影響的泥石流運(yùn)動(dòng)過程，這種極端條件下的泥石流運(yùn)動(dòng)過程在溝口處仍具有較大能量，模擬結(jié)果如圖4（d）所示。

圖4 不同條件下的模擬結(jié)果Fig.4 Simulation results under different conditions

3.1 不同地形分辨率對(duì)泥石流易發(fā)性的影響

選用不同分辨率的DEM數(shù)據(jù)，分析DEM分辨率對(duì)潛在源區(qū)的識(shí)別能力。選用老狼溝流域范圍5、10、20 m 分辨率的DEM 數(shù)據(jù)，采用相同的源區(qū)識(shí)別標(biāo)準(zhǔn)和運(yùn)動(dòng)擴(kuò)散參數(shù)。分析不同地形分辨率下的崩塌滑坡同時(shí)影響下的泥石流災(zāi)害發(fā)生的易發(fā)性評(píng)價(jià)。通過模擬泥石流的運(yùn)動(dòng)過程，發(fā)現(xiàn)DEM 的精度不僅影響源區(qū)的識(shí)別，同時(shí)由于數(shù)據(jù)質(zhì)量的不同，對(duì)擴(kuò)散運(yùn)動(dòng)過程也有一定影響，如圖5。5 m分辨率的DEM最精細(xì)，誤識(shí)別出的坡面微起伏地貌作為潛在源區(qū)發(fā)生類似地表徑流的擴(kuò)散運(yùn)動(dòng)，模擬出的運(yùn)動(dòng)擴(kuò)散范圍較實(shí)際情況有點(diǎn)偏大，造成運(yùn)動(dòng)范圍過于分散，結(jié)果與實(shí)際不符。同時(shí)，5 m分辨率的數(shù)據(jù)計(jì)算耗時(shí)較長(zhǎng)。20 m精度的DEM 數(shù)據(jù)過于粗糙，很多小規(guī)模的運(yùn)動(dòng)過程不能很好地反映出來。10 m 分辨率的DEM 模擬出的結(jié)果最為吻合。

圖5 不同分辨率DEM下同時(shí)受崩塌與滑坡影響的泥石流擴(kuò)散運(yùn)動(dòng)結(jié)果Fig.5 Simulated diffusion movements of debris flow affected by collapse and landslide simultaneously under different DEM resolutions

3.2 不同臨界坡度閾值對(duì)泥石流易發(fā)性的影響

分析15°、25°和35°三種臨界地形坡度閾值，崩塌與滑坡共同影響下的泥石流易發(fā)性結(jié)果。隨著地形坡度閾值的增大, 坡面上識(shí)別處的物源范圍減小，流域內(nèi)的泥石流易發(fā)區(qū)域也隨著坡度閾值的增大而明顯縮?。▓D6）。當(dāng)坡度閾值為35°時(shí)，易發(fā)區(qū)域主要集中在主溝中；當(dāng)坡度閾值為25°時(shí)，易發(fā)區(qū)域主要分布在各溝道內(nèi)和坡面上的大型洼地區(qū)域，當(dāng)坡度閾值為15°時(shí)，流域范圍內(nèi)約2/3區(qū)域?yàn)槟嗍饕装l(fā)區(qū)。因此，需要根據(jù)現(xiàn)場(chǎng)調(diào)查的潛在物源情況，選擇最合適的地形坡度閾值作為確定泥石流易發(fā)性的參數(shù)，識(shí)別出最為接近的潛在物源堆積體。物源的減少造成溝道內(nèi)的泥石流物質(zhì)的運(yùn)動(dòng)動(dòng)能較小，到達(dá)溝口處的能量迅速消失，運(yùn)動(dòng)距離變短。

圖6 不同源區(qū)啟動(dòng)坡度閾值的泥石流運(yùn)動(dòng)結(jié)果Fig.6 Movements of debris flow under different starting slope thresholds

3.3 不同延伸角對(duì)泥石流易發(fā)性的影響

對(duì)10 m 分辨率下識(shí)別出的洪水源區(qū)，通過改變SFLM模型中的延伸角參數(shù)模擬其運(yùn)動(dòng)過程的差異，分析該參數(shù)對(duì)運(yùn)動(dòng)距離的影響。因?yàn)閲?guó)際上大量研究認(rèn)為該參數(shù)的取值多在7°～11°，因此選取7°、9°和11°三個(gè)角度進(jìn)行對(duì)比分析（圖7），模擬出考慮崩塌滑坡影響下的泥石流易發(fā)性評(píng)價(jià)結(jié)果。如圖7 所示，大量的物源提供了充足的能量，即使當(dāng)延伸角設(shè)為11°時(shí)，泥石流依然流出溝口。但對(duì)比發(fā)現(xiàn)，溝口處的能量相較7°和9°的能量有所減少，可推斷出流出溝口后，其運(yùn)動(dòng)距離最短。由于SFLM 模型是建立在簡(jiǎn)化能量平衡的基礎(chǔ)上，該角度越大代表著因摩擦造成的能量損耗也嚴(yán)重。因此，該參數(shù)具有一定的物理基礎(chǔ)，具有較好的實(shí)用性。當(dāng)碎屑物質(zhì)顆粒比較粗糙，顆粒之間的摩擦增強(qiáng)時(shí)，該角度往往較大，運(yùn)動(dòng)距離較短。

圖7 不同延伸角對(duì)泥石流運(yùn)動(dòng)距離的影響Fig.7 Influence of different reach angles on movement distance of debris flow

4 結(jié)論

（1）Flow-R 是基于高程模型數(shù)據(jù)的泥石流經(jīng)驗(yàn)性模型，該模型只需少量的空間數(shù)據(jù)和經(jīng)驗(yàn)性閾值參數(shù)。采用數(shù)據(jù)索引法處理不同的空間數(shù)據(jù)進(jìn)行源區(qū)識(shí)別，基于簡(jiǎn)化的能量守恒定律對(duì)識(shí)別出的源區(qū)進(jìn)行擴(kuò)散運(yùn)動(dòng)模擬。模擬結(jié)果可以較好地為更精確的泥石流研究提供指導(dǎo)。

（2）對(duì)僅受降雨影響的洪水易發(fā)性概率、分別受崩塌和滑坡影響及其共同影響的泥石流易發(fā)性概率進(jìn)行模擬。結(jié)果顯示，作為擴(kuò)散運(yùn)動(dòng)過程的物質(zhì)基礎(chǔ)，源區(qū)識(shí)別是正確進(jìn)行易發(fā)性評(píng)價(jià)的前提和關(guān)鍵。崩塌、滑坡形成的松散堆積體為泥石流運(yùn)動(dòng)提供了良好的物質(zhì)基礎(chǔ)，極大地提高了泥石流的易發(fā)性概率。

（3）選用不同分辨率的DEM 數(shù)據(jù)，分析DEM數(shù)據(jù)分辨率對(duì)泥石流易發(fā)性評(píng)價(jià)的影響，發(fā)現(xiàn)DEM 分辨率過大或過小時(shí)，均會(huì)造成泥石流源區(qū)的誤判，反映出的泥石流擴(kuò)散運(yùn)動(dòng)范圍不是很合理。對(duì)研究區(qū)域，10 m 分辨率的DEM 是最佳的基礎(chǔ)數(shù)據(jù)。

（4）設(shè)定不同的泥石流源區(qū)啟動(dòng)坡度閾值，通過將模擬結(jié)果與現(xiàn)場(chǎng)實(shí)際松散堆積體的對(duì)比，找出合理的泥石流發(fā)生地形條件，發(fā)現(xiàn)泥石流發(fā)生的最低坡度值越小，可作為物源的松散物質(zhì)越多，地形坡度影響沉積物的堆積。

（5）不同的運(yùn)動(dòng)擴(kuò)散方式影響泥石流的擴(kuò)散范圍和運(yùn)動(dòng)距離，基于簡(jiǎn)化能量平衡的SFLM模型可以反映出，碎屑物質(zhì)顆粒比較粗糙，顆粒之間的摩擦增強(qiáng)時(shí)，運(yùn)動(dòng)距離較短。