李 炫，楊本勇，范建容，張建強(qiáng)，李磊磊

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

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

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

1 研究區(qū)概況

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

2 流域單元?jiǎng)澐?/h2>

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

2.1 閾值確定和流域劃分

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

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

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

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

計(jì)算公式為：

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

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

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

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

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

圖2 s與si差值變化曲線

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

2.2 流域修正

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

圖3 岷江上游地區(qū)流域單元?jiǎng)澐?/p>

3 流域單元泥石流危險(xiǎn)性評(píng)價(jià)

3.1 評(píng)價(jià)因子的選取與分級(jí)

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

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

表1 評(píng)價(jià)因子之間的相關(guān)矩陣

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

表2 評(píng)價(jià)因子等級(jí)劃分

3.2 評(píng)價(jià)因子權(quán)重的確定

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

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

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

（3）計(jì)算主成分的累積方差貢獻(xiàn)率，并確定出累積貢獻(xiàn)率大于85%的主成分?jǐn)?shù)；

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

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

表3 評(píng)價(jià)因子權(quán)重值

3.3 流域單元的危險(xiǎn)度評(píng)價(jià)

4 結(jié)果與驗(yàn)證

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

圖4 岷江上游泥石流危險(xiǎn)分布

5 結(jié) 論

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

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

［1］ 張懷珍，范建容，郭芬芬，等.中國(guó)單溝泥石流危險(xiǎn)度評(píng)價(jià)模型比較研究［J］.水土保持研究，2011，18（4）：20－26.

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

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

［4］ 唐川，朱大奎.基于GIS技術(shù)的泥石流風(fēng)險(xiǎn)評(píng)價(jià)研究［J］.地理科學(xué)，2002，22（3）：300－304.

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

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

［7］ 王曉朋，潘懋，徐岳仁.基于流域單元的泥石流區(qū)域危險(xiǎn)性評(píng)價(jià)［J］.山地學(xué)報(bào)，2006，24（2）：177－180.

［8］ 李麗.舟曲縣泥石流危險(xiǎn)性評(píng)價(jià)［D］.北京：中國(guó)地質(zhì)大學(xué)，2012.

［9］ 王瓊.基于流域尺度的震后汶川縣潛在泥石流危險(xiǎn)性評(píng)價(jià)［D］.成都：成都理工大學(xué)，2012.

［10］ 劉希林，蘇鵬程.四川省泥石流風(fēng)險(xiǎn)評(píng)價(jià)［J］.災(zāi)害學(xué)，2004，19（2）：23－28.

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

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

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

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

［15］ 賴晗，芮小平，梁漢東，等.基于均值變點(diǎn)分析的三峽庫區(qū)河網(wǎng)提取研究［J］.測(cè)繪科學(xué)，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ū)泥石流危險(xiǎn)性評(píng)估：以都江堰—汶川公路為例［J］.中國(guó)水土保持科學(xué)，2012，10（1）：6－11.

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

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