基于流域單元的岷江上游泥石流危險性評價

李 炫，楊本勇，范建容，張建強，李磊磊

（1.中國科學院 水利部 成都山地災(zāi)害與環(huán)境研究所，成都610041；2.中國科學院大學，北京100049；3.國家測繪地理信息局 第三地理信息制圖院，成都610100）

泥石流是在暴雨、融雪等外界條件的激發(fā)下產(chǎn)生的，來勢兇猛并且過程短暫，具有運動速度快、暴發(fā)突然、能量巨大等特征，其攜帶的大石塊具有沖擊力強和破壞性大的特點［1］，泥石流通常給山區(qū)人民帶來較大的經(jīng)濟財產(chǎn)損失和人員傷亡，如2010年8月7日甘肅省舟曲縣城后山三眼溝及羅家溝突然暴發(fā)大規(guī)模泥石流，造成了1 744人死亡和失蹤［2］。因此，做好泥石流的評估、預(yù)警工作，對防范泥石流災(zāi)害的發(fā)生，減少災(zāi)害帶來的損失，以及區(qū)域的規(guī)劃建設(shè)都具有重要的意義。

泥石流危險性評價是災(zāi)情評估、預(yù)測、防災(zāi)救災(zāi)的基礎(chǔ)［3］。在國內(nèi)好多區(qū)域泥石流危險性評價工作中大都采用柵格單元［4－6］，其結(jié)構(gòu)簡單計算高效，被國內(nèi)學者廣泛的應(yīng)用于泥石流危險性評價中。但這種評價單元與地質(zhì)地貌以及泥石流發(fā)生的其他環(huán)境條件缺乏聯(lián)系，無法體現(xiàn)泥石流發(fā)生的實際情況。而泥石流溝谷的地形地貌特征、豐富的物質(zhì)條件以及降雨是泥石流發(fā)生的必要條件［7］，基于流域單元的泥石流評價能充分考慮到與泥石流形成相關(guān)的內(nèi)在因子的影響。目前流域劃分方法大都是基于DEM數(shù)據(jù)［8－9］，提取流域時閾值的確定大都通過人工嘗試，沒有一個較為精確合理的方法，自動化程度不高，會產(chǎn)生的大量非閉合子流域，劃分出的流域與實際的泥石流流域有所差別。本文選取岷江上游區(qū)域作為研究區(qū)，采用均值變點法和河網(wǎng)密度法探究研究區(qū)流域劃分的合理閾值，并以劃分的流域單元為基礎(chǔ)，從泥石流的形成條件出發(fā)，選取地層巖性、斷裂密度、平均坡度、形狀系數(shù)、溝道比降、24h最大降雨量等六個指標，對岷江上游地區(qū)泥石流的危險性進行綜合的評估。

1 研究區(qū)概況

岷江上游地區(qū)轄四川省阿壩藏族羌族自治州的汶川、茂縣、理縣、黑水和松潘五縣的全部和都江堰小部分區(qū)域，流域面積為23 037km2，流域內(nèi)地勢起伏較大，最高海拔6 246m，最低海拔455m，該區(qū)地貌類型以山原和高山峽谷為主，山地災(zāi)害嚴重，泥石流災(zāi)害處于高風險區(qū)，是四川省乃至是全國泥石流的重災(zāi)區(qū)［10］，地質(zhì)結(jié)構(gòu)復雜，斷裂發(fā)育，龍門山構(gòu)造體系、岷江縱向構(gòu)造體系共同作用，茂汶斷裂、北川—映秀斷裂、岷江斷裂等深大斷裂也分布于此，加上松潘地震帶和龍門山地震帶等活動地震帶的影響，歷史上發(fā)生過許多強大的地震如1976年的松潘地震和2008年的汶川地震等。受地震等災(zāi)害的影響，研究區(qū)內(nèi)地質(zhì)、地形條件都發(fā)生了不同程度的改變，地表植被損毀，松散固體物質(zhì)增多，使得區(qū)域內(nèi)地表抗蝕能力減弱，易損性變大，這些因素都為泥石流的形成創(chuàng)造了良好的條件。同時，該區(qū)在南北走向的特殊地貌和西南季風的共同作用下，焚風效應(yīng)顯著［11］，降雨季節(jié)性明顯，5—8月集中了全年80%～90%的雨量，汛期時常會有暴雨出現(xiàn)，山洪和泥石流災(zāi)害頻發(fā)。對岷江上游地區(qū)進行泥石流危險性評價，能為區(qū)域內(nèi)泥石流災(zāi)害的預(yù)防提供參考，對減少人民經(jīng)濟財產(chǎn)損失和區(qū)域的長遠發(fā)展都具有非常重要的意義。

2 流域單元劃分

數(shù)字高程模型（Digital Elevation Model，DEM）是反映流域地形特征的數(shù)字模型，基于DEM的流域劃分，大大節(jié)省了的人力、物力，從提取的效率和精度來看都是切實可行的［12］。基于DEM的流域劃分需要對數(shù)據(jù)進行填洼處理，用網(wǎng)格坡度確定每個單元格的水流方向，得到對應(yīng)的流向矩陣，計算匯流累積量，再給定一個匯水閾值，凡是匯水大于該閾值的網(wǎng)格，均為河網(wǎng)內(nèi)的網(wǎng)格，將這些網(wǎng)格連接后，得出整個流域的排水網(wǎng)絡(luò)。當流域的排水網(wǎng)絡(luò)生成以后，就可以確定整個流域的界限。因此，閾值的確定是河網(wǎng)提取和流域劃分的關(guān)鍵。

2.1 閾值確定和流域劃分

閾值的確定方法有多種，主要有水系分形分維法［13］、河網(wǎng)密度法［14］、均值變點法［15］等等。本文以空間分辨率為25m×25m的DEM數(shù)據(jù)作為數(shù)據(jù)源，綜合運用均值變點法和河網(wǎng)密度法推求閾值，避免閾值選取的隨意性。河網(wǎng)密度法認為閾值與河網(wǎng)密度（或河源密度）關(guān)系曲線趨于平緩時所對應(yīng)的點為合理的閾值。如圖1所示，為河網(wǎng)密度與閾值之間的變化關(guān)系。由圖1可以看出：隨著閾值的不斷增大，總河網(wǎng)長度會逐漸減少，河網(wǎng)密度也逐漸減少；起初河網(wǎng)密度大小減少比較劇烈，隨著閾值的慢慢增大，河網(wǎng)密度的減少變得越來越緩，當閾值為4 000的時候，二階導數(shù)趨近于0，河網(wǎng)密度變化趨于穩(wěn)定。

圖1 河網(wǎng)密度隨閾值變化趨勢

均值變點法的離散數(shù)據(jù)模型為［15］：

a1＝…am1－1＝b1，am1＝…am2－1＝b2，amq＝…＝aN＝bq＋1，此處1＜m1＜m2＜…＜mq≤N，如果bj＋1≠bj，則mj就是一個變點。

計算公式為：

（1）令i＝2，…，N對每個i將樣本分為2段：x1，x2，…，xi－1和xi，xi＋1，…，xN。

式中：Si——各段樣本的方差和；Xi1，Xi2——每段樣本的算術(shù)平均值。

（2）計算統(tǒng)計量。

式中：xm——樣本均值；s——總體的方差；xt——單個樣本值。

選取平均河網(wǎng)密度作為樣本序列X，按照公式（3）計算樣本總體平均值和方差由公式（2）計算出分段樣本的算術(shù)平均值與方差，從而得到S與Si差值變化曲線如圖2所示。變點的存在會使S和Si的差值增大，這個方法對恰有一個變點的檢驗最為有效［11］。由圖2可以看出第5個點對應(yīng)的差值最大，說明第5個點就是我們要找的變點，第5個點對應(yīng)的閾值為4 000，這與之前用河網(wǎng)密度法所求的閾值相同。

圖2 s與si差值變化曲線

根據(jù)計算的閾值提取出河網(wǎng)，利用GIS工具Hydrology Tools劃分出子流域單元。將劃分出的流域單元與實際的流域單元進行比較發(fā)現(xiàn)，大部分提取的流域單元與實際的流域單元都是比較接近的，但提取出的流域中包含大量的非閉合子流域，這些非閉合的子流域會對泥石流的危險性評價造成了干擾。同時，流域單元劃分時只基于地形數(shù)據(jù)，提取出的流域單元在某些地勢平坦的區(qū)域，流域邊界并不準確，會出現(xiàn)整個流域被刨分成多個流域，某些細小的流域被劃分成一個流域等情況，所以，應(yīng)根據(jù)實際地形對流域單元進行修正。

2.2 流域修正

針對上述提取的流域邊界與實際流域單元不吻合的情況，本文運用DEM生成的山體陰影圖像結(jié)合高分辨率遙感影像，對部分集水單元進行手動修改，去除某些非閉合子流域，按照流域的實際形態(tài)合并或修正子流域邊界，最終得到岷江上游地區(qū)的小流域分布（圖3）。

圖3 岷江上游地區(qū)流域單元劃分

3 流域單元泥石流危險性評價

3.1 評價因子的選取與分級

影響泥石流形成的因素較多，根據(jù)形成機制來看，影響因素主要可以分為：地形地貌條件、物源條件、降雨等誘發(fā)條件。地形地貌條件制約著泥石流的形成、運動和規(guī)模等特征，地形因素包括流域面積、縱比降、坡度、形狀因子、發(fā)育程度等［16］。韓用順等人就選取了巖性、平均坡度、流域面積、相對高差、主溝長度等因子實現(xiàn)了汶川震后災(zāi)區(qū)泥石流危險度評價［17］。譚萬沛等在研究區(qū)域暴雨泥石流預(yù)測時表明形狀系數(shù)和溝道比降對泥石流的發(fā)生具有很大的制約作用，是影響泥石流流速和輸送能力大小的重要條件，是評價泥石流流域危險度的重要因子［18］。孟國才等人在研究岷江上游泥石流特征時指出岷江上游地區(qū)地質(zhì)構(gòu)造復雜，斷裂發(fā)育，新構(gòu)造運動隆升強烈，茂汶斷裂、北川—映秀斷裂、岷江斷裂等深大斷裂分布于此，地層巖性和地質(zhì)構(gòu)造是岷江上游地區(qū)泥石流危險性評價的重要因素［11］。降雨是誘發(fā)泥石流災(zāi)害最直接的因素，對泥石流的形成起關(guān)鍵作用的是1h雨強和10min雨強，但是小時雨強和10min雨強不易觀測和記錄，而日降雨量的相對容易些，且三者之間存在正相關(guān)關(guān)系，一定程度上可以用日最大降水量來同時表示雨量和雨強。

本文在上述分析和前人研究的基礎(chǔ)上，從影響泥石流的三大因素出發(fā)，篩選出對岷江上游地區(qū)泥石流發(fā)生起一定主導作用的地層巖性、斷裂密度、平均坡度、形狀系數(shù)、溝道比降、24h最大降雨量等六個指標建立泥石流危險性的評價模型。通過對研究區(qū)978個流域的6個指標進行量化，并進行相關(guān)分析，以避免影響因子的重復引入，相關(guān)矩陣如表1所示。從相關(guān)矩陣可以看出各個因子之間沒有明顯的相關(guān)關(guān)系，均從不同角度反映了流域的特征。

表1 評價因子之間的相關(guān)矩陣

由于各個評價因子的計量單位不同，取值范圍變幅很大，會影響評價的結(jié)果，所以需要對上述指標進行歸一化處理。本文對各因子采用5級量化指標，通過不同級別反映各因子對泥石流發(fā)育影響程度的差異，等級越高，說明其對泥石流的影響程度越大，泥石流危險性越高。對于地層巖性因子，由獲得的研究區(qū)地質(zhì)巖性圖，綜合考慮本區(qū)巖土體類型、物理力學性質(zhì)，并結(jié)合巖土體結(jié)構(gòu)特征，參考工程巖體分級標準，將研究區(qū)內(nèi)巖組分為堅硬巖組（花崗巖、正長巖、石英等）、較堅硬巖組（白云巖、大理巖等）、軟硬相間巖組（千枚巖、砂質(zhì)泥巖等）、軟弱巖組（含煤層、半固結(jié)巖）、極軟巖組（泥巖和煤系地層）五類，危險度分別賦予1，2，3，4，5。對于斷裂密度、平均坡度、形狀系數(shù)、溝道比降、24h最大降雨量等評價因子，采用自然間斷點分級法，劃分為五個等級，并按照各指標值的大小分別賦值為5，4，3，2，1。具體歸一化方案見表2。

表2 評價因子等級劃分

3.2 評價因子權(quán)重的確定

在評價指標體系中，不同影響因子對泥石流災(zāi)害的貢獻程度是不同的［6］，這種貢獻程度的不同通過各評價因子的權(quán)重系數(shù)加以體現(xiàn)。目前泥石流評價中權(quán)重確定的方法一般有兩種，主觀賦權(quán)法，和客觀賦權(quán)法。主觀賦權(quán)法如層次分析法、專家打分法等，權(quán)重的結(jié)果受到人為因素影響較大。而客觀賦權(quán)法如主成分分析法和因子分析法等是根據(jù)各指標間的相關(guān)關(guān)系或各項指標值的變異程度來確定權(quán)重系數(shù)，可以避免人為因素的影響。本文采用主成分分析法來確定個因子的權(quán)重，計算方法如下［19］：

（1）對各因子進行相關(guān)分析，確定相關(guān)系數(shù)矩陣R（又稱協(xié)方差矩陣）；

（2）計算矩陣R的特征值和特征向量，即通過正交變換將矩陣化為對角矩陣，對矩陣中的對角元素即為所求特征值，按其大小排列，分別稱為第一、二特征值（λ1＞λ2＞…λn）；

（3）計算主成分的累積方差貢獻率，并確定出累積貢獻率大于85%的主成分數(shù)；

（4）計算主成分中載荷系數(shù)，即特征值的方根與相應(yīng)的特征向量的乘積，并將其歸一化，得到各參評因子的權(quán)重。

將量化的978個流域數(shù)據(jù)作主成分分析，以累積方差大于93%確定主成分數(shù)（4個），確定出各個評價因子的權(quán)重值，見表3。

表3 評價因子權(quán)重值

3.3 流域單元的危險度評價

4 結(jié)果與驗證

通過加入研究區(qū)內(nèi)已查明的205條泥石流溝，與評價結(jié)果相疊加，結(jié)果如圖4所示。研究區(qū)內(nèi)有9條泥石流溝位于極低危險區(qū)內(nèi)，占總數(shù)的4.39%，21條泥石流溝位于低度危險區(qū)域內(nèi)，占總數(shù)的10.24%，52條泥石流溝位于中度危險區(qū)域內(nèi)，占總數(shù)的25.36%，84條泥石流溝位于高度危險區(qū)內(nèi)，占總數(shù)的40.98%，39條溝位于極度危險區(qū)域內(nèi)，占總數(shù)的19.03%。已查明的泥石流溝一共有175條位于中度危險、高度危險、極度危險區(qū)域內(nèi)，占總數(shù)的85.37%。結(jié)果分析表明，評價的結(jié)果能較好的反映出岷江上游的泥石流分布情況。

圖4 岷江上游泥石流危險分布

5 結(jié) 論

本文以DEM數(shù)據(jù)為基礎(chǔ)進行流域劃分，綜合運用均值變點法和河網(wǎng)密度法確定出岷江上游地區(qū)流域劃分的合理閾值為4 000，將自動劃分的流域單元，參考實際的遙感影像，對子流域單元進行合并或修正，最終得到岷江上游地區(qū)大小流域共978個。通過遙感影像的修正，劃分的流域與實際的流域比較吻合。

從影響泥石流形成的三大條件出發(fā)，選取地層巖性、斷裂密度、平均坡度、形狀系數(shù)、溝道比降、24h最大降雨量等六個評價指標，運用主成分分析法確定各影響因子的權(quán)重值，對流域的量化數(shù)據(jù)進行歸一化、分級處理，在ARCGIS支持下完成了岷江上游地區(qū)泥石流危險性評價，得到危險性評價圖，其中極度危險度流域97個，占總流域的9.92%，高度危險流域205個，占總流域的20.96%，中度危險流域286個，占總流域的29.24%，低度危險流域240個，占總流域的24.54%，極低危險流域150個，占總流域的15.34%。用已有災(zāi)害數(shù)據(jù)驗證表明，評價結(jié)果能較好的反映出岷江上游地區(qū)泥石流的分布規(guī)律，為該區(qū)泥石流災(zāi)害的防治提供了參考和依據(jù)。

［1］ 張懷珍，范建容，郭芬芬，等.中國單溝泥石流危險度評價模型比較研究［J］.水土保持研究，2011，18（4）：20－26.

［2］ 余斌，楊永紅，蘇永超，等.甘肅省舟曲8.7特大泥石流調(diào)查研究［J］.工程地質(zhì)學報，2010，18（4）：437－444.

［3］ 孫紹騁.災(zāi)害評估研究內(nèi)容與方法探討［J］.地理科學進展，2001，20（2）：122－130.

［4］ 唐川，朱大奎.基于GIS技術(shù)的泥石流風險評價研究［J］.地理科學，2002，22（3）：300－304.

［5］ 寧娜，馬金珠，張鵬，等.基于GIS和信息量法的甘肅南部白龍江流域泥石流災(zāi)害危險性評價［J］.資源科學，2013，35（4）：892－899.

［6］ 王歡，丁明濤，陳廷方.基于GIS的三江并流區(qū)泥石流危險性評價［J］.水土保持通報，2011，31（5）：167－170.

［7］ 王曉朋，潘懋，徐岳仁.基于流域單元的泥石流區(qū)域危險性評價［J］.山地學報，2006，24（2）：177－180.

［8］ 李麗.舟曲縣泥石流危險性評價［D］.北京：中國地質(zhì)大學，2012.

［9］ 王瓊.基于流域尺度的震后汶川縣潛在泥石流危險性評價［D］.成都：成都理工大學，2012.

［10］ 劉希林，蘇鵬程.四川省泥石流風險評價［J］.災(zāi)害學，2004，19（2）：23－28.

［11］ 孟國才，王士革，謝洪，等.岷江上游泥石流災(zāi)害特征分析［J］.災(zāi)害學，2005，20（3）：94－98.

［12］ 陳加兵，勵惠國，鄭達賢，等.基于DEM的福建省小流域劃分研究［J］.地球信息科學，2007，9（2）：74－77.

［13］ 楊邦，任立良.集水面積閾值確定方法的比較研究［J］.水電能源科學，2009，27（5）：11－14.

［14］ 孔凡哲，李莉莉.利用DEM提取河網(wǎng)時集水面積閾值的確定［J］.水電能源科學，2006，23（4）：65－67.

［15］ 賴晗，芮小平，梁漢東，等.基于均值變點分析的三峽庫區(qū)河網(wǎng)提取研究［J］.測繪科學，2012，5（37）：056.

［16］ Liu C N，Dong J J，Peng Y F，et al.Effects of strong ground motion on the susceptibility of gully type debris flows［J］.Engineering Geology，2009，104（3）：241－253.

［17］ 韓用順，李龍偉，朱穎彥，等.汶川地震災(zāi)區(qū)泥石流危險性評估：以都江堰—汶川公路為例［J］.中國水土保持科學，2012，10（1）：6－11.

［18］ 譚萬沛，王成華，姚令侃，等.暴雨泥石流滑波的區(qū)域預(yù)測與預(yù)報：以攀西地區(qū)為例［M］.成都：四川科學技術(shù)出版社，1994.

［19］ 韓小孩，張耀輝，孫福軍，等.基于主成分分析的指標權(quán)重確定方法［J］.四川兵工學報，2012，33（10）：124－126.