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

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

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

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

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

1 研究區概況

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

2 流域單元劃分

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

2.1 閾值確定和流域劃分

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

圖1 河網密度隨閾值變化趨勢

均值變點法的離散數據模型為［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——每段樣本的算術平均值。

（2）計算統計量。

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

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

圖2 s與si差值變化曲線

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

2.2 流域修正

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

圖3 岷江上游地區流域單元劃分

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

3.1 評價因子的選取與分級

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

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

表1 評價因子之間的相關矩陣

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

表2 評價因子等級劃分

3.2 評價因子權重的確定

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

（1）對各因子進行相關分析，確定相關系數矩陣R（又稱協方差矩陣）；

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

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

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

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

表3 評價因子權重值

3.3 流域單元的危險度評價

4 結果與驗證

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

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

5 結 論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

［15］ 賴晗，芮小平，梁漢東，等.基于均值變點分析的三峽庫區河網提取研究［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］ 韓用順，李龍偉，朱穎彥，等.汶川地震災區泥石流危險性評估：以都江堰—汶川公路為例［J］.中國水土保持科學，2012，10（1）：6－11.

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

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