基于邏輯回歸模型的泥石流易發(fā)性評價與檢驗：以金沙江上游奔子欄—昌波河段為例

吳賽兒，陳 劍，ZHOU Wendy，高玉欣，徐能雄

(1.中國地質(zhì)大學(xué)(北京) 工程技術(shù)學(xué)院，北京 100083; 2.Colorado School of Mines, Department of Geology and Geological Engineering, Colorado Denver 80401)

0 引 言

泥石流作為我國最頻發(fā)的地質(zhì)災(zāi)害之一，多發(fā)生于半干旱山區(qū)以及高原冰川區(qū)，其分布遍及我國31個省、市和自治區(qū)，全國有1 583個縣(市)長期受到泥石流災(zāi)害的困擾，每年要暴發(fā)數(shù)千起泥石流災(zāi)害，其活動區(qū)域面積達480萬km2，共計8萬多處，其中活動強烈的泥石流有8 500多處，面積達130萬km2[1-2]，多爆發(fā)于西南、西北的廣大山區(qū)，在空間和時間上有分布廣泛但相對集中的特征。近年來，頻發(fā)的泥石流災(zāi)害直接造成了平均每年近千的人員傷亡和>10億元的經(jīng)濟損失[3]。因此，以探究區(qū)域泥石流特點、分析在特定環(huán)境下的致災(zāi)因子、進行泥石流預(yù)測與評價等為主的研究工作備受關(guān)注。作為風(fēng)險評價的基礎(chǔ)，通過泥石流的易發(fā)性評價來直觀反映泥石流災(zāi)害的空間分布，不僅可以完善當(dāng)前災(zāi)害程度的評估，同時也為預(yù)防潛在災(zāi)害的戰(zhàn)略準(zhǔn)備提供參考，對于預(yù)測預(yù)警、防災(zāi)減災(zāi)工作意義重大。目前，國內(nèi)外不少學(xué)者開展了廣泛的研究，在泥石流的形成條件、堆積特征、影響因素、評價方法方面取得了重要成果[4-20]。

本文應(yīng)用多元統(tǒng)計分析方法中的邏輯回歸(LR)模型，以GIS為操作平臺，基于小流域為評價單元開展干熱河谷泥石流易發(fā)性評價與檢驗。LR模型的自變量既可以是連續(xù)變量，也可以是二分類或多分類變量，易于用來預(yù)測具有二項特點的概率問題，并且與區(qū)域遙感影像資料的結(jié)合有很好的適應(yīng)性和應(yīng)用價值。隨著LR模型和GIS技術(shù)應(yīng)用于地質(zhì)災(zāi)害空間預(yù)測的理念被國外學(xué)者Gorsevski提出后[21]，LR模型的可操作性在地質(zhì)災(zāi)害空間預(yù)測方面已得到很好的證實，尤其在滑坡災(zāi)害敏感性評價中得到較廣泛的應(yīng)用[21-27]。近年來該方法在泥石流評價方面的應(yīng)用逐漸引起重視[28-30]。

本文研究區(qū)為金沙江上游奔子欄—昌波河段一帶，覆蓋面積3 105 km2，地形上屬于高山峽谷地區(qū)，切割強烈[31]。作為我國“西部大開發(fā)”戰(zhàn)略的重要地區(qū)之一，該區(qū)域內(nèi)泥石流等地質(zhì)災(zāi)害的頻發(fā)對周邊的水利水電工程和公路工程建設(shè)構(gòu)成嚴(yán)重威脅。目前，在研究區(qū)泥石流堆積扇成因機制和泥石流易發(fā)性評價方面已經(jīng)取得一些顯著成果[32-36]，但就泥石流預(yù)測評價而言，在評價指標(biāo)的選取、標(biāo)準(zhǔn)化處理、評價結(jié)果精度及檢驗方法等方面尚有待改進。為彌補前期評價工作的不足，探討如何有效地對評價結(jié)果進一步驗證，本文采用邏輯回歸模型對金沙江上游奔子欄—昌波河段的泥石流易發(fā)性進行了評價和檢驗。

1 區(qū)域地質(zhì)概況及泥石流分布

1.1 研究區(qū)地質(zhì)概況

研究區(qū)分布于四川、云南、西藏交接處的金沙江上游徐龍—奔子欄一帶(圖1)，地處我國西南橫斷山區(qū)；在地貌單元上，位于以“世界屋脊”著稱的青藏高原南東緣，橫斷山脈中段；地理坐標(biāo)為98°55′—99°25′E，28°00′—29°20′N。

圖1 研究區(qū)位置與地貌圖Fig.1 Location and geomorphic map of the study area

地理位置對年降水量分布的影響主要體現(xiàn)在谷地對外來水汽的遮蔽。研究區(qū)峽谷深切、高山對峙，河谷相對高差達1～1.5 km[31]，受西南季風(fēng)和東南暖濕氣流影響水汽難以進入，研究區(qū)處于“雨影區(qū)”，表現(xiàn)為典型干熱河谷型氣候。在金沙江河谷的奔子欄到四川得榮一帶，年降水量只有300 mm左右，巴塘至得榮間金沙江河谷年降水量少于400 mm[36]。降水季節(jié)分布極不均勻，6月至9月雨量非常集中，一般占80%以上，雨季旱季分化明顯。同時，研究區(qū)氣溫在時間和空間上，呈現(xiàn)極大的垂直分布溫差變化。每年6月份河谷最高氣溫可達30 ℃以上，每年1月至2月高山區(qū)最低氣溫則降至-30 ℃。該區(qū)特殊的自然條件造成植被不發(fā)育、地表風(fēng)化剝蝕強烈[31]。

研究區(qū)位于瀕臨歐亞板塊與印度板塊相互碰撞帶的東緣，沉積構(gòu)造環(huán)境復(fù)雜多樣，火山活動極其強烈、頻繁，分布多條構(gòu)造巖漿巖帶、變質(zhì)帶、緊密褶皺沖斷層帶，但地震活動性相對較弱，震級均小于6級。區(qū)內(nèi)巖性復(fù)雜?？v貫?zāi)媳钡慕鹕辰瓟嗔褜⒃搮^(qū)分為東、西兩個地層區(qū)，東區(qū)出露地層為石炭系、上二疊統(tǒng)和下二疊統(tǒng)、中—下三疊統(tǒng)，巖性主要為灰?guī)r、板巖、砂巖、片巖。二疊系地層在西區(qū)分布較廣，中—下三疊統(tǒng)、二疊系下統(tǒng)下部地層主要分布在西區(qū)東北部以及金沙江沿岸，主要出露巖性為板巖、砂巖、火山巖、灰?guī)r[23,25]。研究區(qū)內(nèi)滑坡、邊坡失穩(wěn)、崩塌等不良地質(zhì)現(xiàn)象發(fā)育，地質(zhì)構(gòu)造帶破碎，巖石風(fēng)化剝蝕強烈，全區(qū)物理風(fēng)化現(xiàn)象較為強烈(圖2)，為泥石流提供了充足的松散破碎固體物源。

1.2 泥石流分布

通過遙感解譯和現(xiàn)場調(diào)查資料分析，在研究區(qū)范圍內(nèi)，直接進入金沙江干流河谷和主要支流的泥石流共有91條(圖2)，共涉及流域面積約1 274 km2，占研究區(qū)總面積的44.5%。其中溝谷型泥石流66條，廣泛發(fā)育于全區(qū)，涉及流域面積約1 242 km2，平均流域面積18.8 km2；坡面泥石流25條，分布于金沙江沿岸的斜坡地帶，涉及流域面積32 km2，平均流域面積為1.28 km2[35]。

圖2 研究區(qū)泥石流流域分布圖Fig.2 Distribution map of debris flow catchment in the study area

2 數(shù)據(jù)準(zhǔn)備與研究方法

2.1 基礎(chǔ)數(shù)據(jù)與評價單元

本文使用的數(shù)據(jù)主要包括30 m×30 m分辨率的DEM、Landsat 8遙感影像、1∶5萬區(qū)域地形圖和地質(zhì)圖、降水資料等。研究區(qū)被30 m×30 m的方格分為3 472 074個柵格單元，其中泥石流柵格單元1 545 073個。由于泥石流發(fā)生區(qū)域的平面形狀多呈狹長展布，以往普遍以柵格或規(guī)則網(wǎng)格為單元的泥石流評價方法即便在GIS軟件中易于劃分與計算，但其存在割裂流域整體性的弊端，因此，本文以不規(guī)則小流域作為單元開展評價，以提高易發(fā)性評價區(qū)劃和檢驗的精度。

2.2 方法與步驟

研究工作分為5步進行：

(1)通過GIS平臺，基于研究區(qū)30 m×30 m分辨率的數(shù)字高程模型(DEM)，綜合考慮野外調(diào)查、航空影像、交通圖中河流分布現(xiàn)狀，采用GIS水文分析模塊將研究區(qū)劃分為220個小流域，根據(jù)現(xiàn)場踏勘和遙感影像解譯的泥石流分布點，進而分為91個泥石流樣本和129個非泥石流樣本；

(2)選擇可能促使歷史泥石流發(fā)生的評價指標(biāo)并借助軟件(ArcGIS10.0和ENVI 5.0)，通過DEM數(shù)據(jù)、數(shù)字地質(zhì)圖、降水資料、Landsat 8遙感影像得到指標(biāo)圖層并柵格化；

(3)對各項指標(biāo)標(biāo)準(zhǔn)化處理后的流域樣本通過SPSS軟件進行75%的隨機抽樣，創(chuàng)建訓(xùn)練組和檢驗組數(shù)據(jù)，前者用來構(gòu)建LR模型，后者用來檢驗?zāi)Ｐ托Ч?/p>

(4)嘗試指標(biāo)全部組合并利用ROC曲線測試模型效果，擬合最佳的LR模型被采用；

(5)利用模型得到易發(fā)性評價分區(qū)圖，最后通過檢驗測試評價結(jié)果的精度與可靠性。

2.3 評價指標(biāo)選取

為保證選取指標(biāo)的有效性，基于LR模型的泥石流易發(fā)性評價在兩個假設(shè)前提下進行：(1)預(yù)測的潛在泥石流將發(fā)生在與歷史泥石流相同或近似環(huán)境中；(2)在分析中使用的指標(biāo)在未來相當(dāng)長的時間內(nèi)不會發(fā)生明顯改變。

根據(jù)專家打分法和前人泥石流評價資料，初步選取了泥石流危險性評價中常用的評價指標(biāo)，并進行了相關(guān)性分析以確認(rèn)被選指標(biāo)對泥石流發(fā)生起到了控制作用。初步選擇坡度(F1)、坡向(F2)、流域相對高差(F3)、溝壑密度(F4)、巖性(F5)、植被歸一化指數(shù)(F6)、距斷裂帶距離(F7)、距主干道路距離(F8)、地表粗糙度(F9)、雨季月平均降雨量(F10)，共10項作為評價指標(biāo)。

2.4 指標(biāo)的分類與標(biāo)準(zhǔn)化處理

將選取的10項評價指標(biāo)作為一級指標(biāo)然后再進行分類，得到的二級指標(biāo)使用了不同的量綱，有必要進行標(biāo)準(zhǔn)化處理后得到統(tǒng)一尺度的標(biāo)準(zhǔn)指標(biāo)值，再作為回歸變量導(dǎo)入模型以發(fā)揮作用。指標(biāo)標(biāo)準(zhǔn)化消除了連續(xù)變量(坡度、流域相對高差等)和離散變量(巖性)的數(shù)據(jù)類型差異，簡化了在建立模型時對于離散變量進行啞變量設(shè)置這一步驟，同時，也為后續(xù)指標(biāo)的評價分析提供直觀的數(shù)據(jù)。

LR模型指標(biāo)量化方法多歸于兩種，第一種按照公式(1)計算標(biāo)準(zhǔn)化指標(biāo)值Xi,j：

(1)

其中：i=1，2，3，…，n，為一級指標(biāo)編號；j=1，2，3，…，m，為二級指標(biāo)序號；si,j是各二級指標(biāo)下已發(fā)泥石流面積，km2。在某些研究中，si,j被替換為ni,j，表示各二級指標(biāo)下已發(fā)泥石流溝谷數(shù)。

第二種指標(biāo)量化方法[30]中，引入各二級指標(biāo)分類后面積Si,j以及初始指標(biāo)值Mi,j，計算方法如下：

(2)

(3)

第一種方法僅僅反映了各二級分類指標(biāo)對災(zāi)害的影響程度；第二種方法不僅考慮了各二級指標(biāo)下災(zāi)害的分布情況，也體現(xiàn)了二級指標(biāo)本身的分布情況，本文采用第二種方法進行指標(biāo)標(biāo)準(zhǔn)化。

2.5 邏輯回歸模型

二元邏輯回歸模型是一種反映二值分類響應(yīng)變量(因變量)對一個或多個回歸變量 (自變量) 進行分析的對數(shù)模型，其結(jié)果反映某種事件的發(fā)生概率。在泥石流易發(fā)性評價中，二值響應(yīng)變量代表泥石流事件發(fā)生與否，即1代表泥石流發(fā)生，0代表不發(fā)生[30]；回歸變量為各評價指標(biāo)數(shù)據(jù)，可以是連續(xù)變量，抑或是分類變量，以此構(gòu)建評價模型。

假設(shè)泥石流災(zāi)害發(fā)生概率記作P，值域分布[0，1]，取P/(1-P)的自然對數(shù)，即對P作Logit轉(zhuǎn)換。記評價指標(biāo)集[x1，x2，…，xn]為自變量，建立LR方程：

lgP=Z=C+B1x1+B2x2+…+Bnxn

(4)

(5)

其中：C為常數(shù)，表示無災(zāi)害隱患條件下響應(yīng)變量事件發(fā)生概率與不發(fā)生概率之比的對數(shù)值；B1，B2，…，Bn為邏輯回歸系數(shù)，表示某個指標(biāo)改變一個單位量時，響應(yīng)變量事件發(fā)生概率與不發(fā)生概率之比的對數(shù)值，通過對訓(xùn)練樣本的求解得到其值。

已知邏輯回歸系數(shù)，利用指標(biāo)值x1，x2，…，xn，便可計算各流域單元災(zāi)害發(fā)生概率P，進而根據(jù)P值劃分災(zāi)害易發(fā)性等級。

3 奔子欄—昌波河段泥石流易發(fā)性評價分區(qū)與檢驗

3.1 評價指標(biāo)分析

將初步選取的10項一級評價指標(biāo)通過等間隔劃分，得到各一級指標(biāo)下的二級指標(biāo)分區(qū)。隨后依據(jù)公式(2)和公式(3)進行標(biāo)準(zhǔn)化處理，并對計算獲得的各標(biāo)準(zhǔn)化指標(biāo)值重賦值。賦值越大表示二級指標(biāo)對泥石流災(zāi)害事件發(fā)生的貢獻越大，賦值范圍[1，j]。指標(biāo)標(biāo)準(zhǔn)化及重賦值結(jié)果如表1所示。由表1可見各二級指標(biāo)對研究區(qū)地形地貌的控制程度(Mi,j)、各二級指標(biāo)對研究區(qū)泥石流發(fā)育的影響程度(Xi,j)，以及標(biāo)準(zhǔn)化指標(biāo)與泥石流發(fā)育的相關(guān)性等級(重賦值)。由此對單一指標(biāo)進行致災(zāi)影響性分析。

(1)坡度：地面坡度在地表徑流、地下水補給與排泄、松散物源堆積等方面起到重要的控制作用，從而間接影響泥石流的發(fā)育。研究區(qū)坡度分布區(qū)間為[0°，67°]，泥石流發(fā)育的優(yōu)勢坡度為20°～30°。當(dāng)坡度<20°不利于泥石流的啟動；坡度>40°時，不利于碎屑物源的積累。

(2)坡向：坡向控制著山坡小氣候和水熱比的分布規(guī)律。陽坡和陰坡相比，陽坡受到日照強度和輻射較多，冰雪消融快，巖石風(fēng)化速度快、程度高，巖土結(jié)構(gòu)較為松散，因此泥石流物源較為豐富。造成陽坡泥石流一般比陰坡發(fā)育的現(xiàn)象。由表1可見，研究區(qū)內(nèi)各坡向均有發(fā)育，但陽坡泥石流發(fā)育略具優(yōu)勢。

(3)流域相對高差：高差變化控制著地表植被發(fā)育和基巖裸露情況，研究區(qū)位于金沙江上游，山勢陡峻，河流下切強烈，高山峽谷地形十分發(fā)育，為泥石流啟動提供動力條件。其中，相對高差在600～900 m范圍內(nèi)泥石流發(fā)育比較集中。

(4)溝壑密度：該指標(biāo)表征研究區(qū)的水系發(fā)育與分布。溝壑在地表水的下蝕、搬運和沉積的交替作用下形成，同時也受地質(zhì)構(gòu)造、地層巖性、氣候等因素的影響，受到內(nèi)外地質(zhì)營力的共同作用。研究區(qū)內(nèi)，隨溝壑密度增大，泥石流發(fā)育的優(yōu)勢成“U”形分布，即溝壑密度較大或較小時均有利于泥石流發(fā)生。

(5)巖性：根據(jù)風(fēng)化嚴(yán)重程度和巖性破碎程度將研究區(qū)地層巖性分類。區(qū)內(nèi)出露的板巖、片巖、砂巖、灰?guī)r風(fēng)化嚴(yán)重，形成的坡積物為泥石流提供主要的物源。由表1可見，處于強烈風(fēng)化等級地區(qū)的泥石流最為發(fā)育。本區(qū)第四系松散堆積物零星分布，所占比例較少，不足以成為泥石流物質(zhì)的主要來源。

(6)植被歸一化指數(shù)(NDVI)：植被增強了巖土體抗剪強度和土壤固結(jié)程度，防止水土流失，是影響并判斷泥石流分布的重要因素之一。全區(qū)NDVI值域為[0，0.7]，在該值處于0.1～0.3之間時，地表多為基巖出露與裸土地帶，有利于泥石流發(fā)育。

(7)距斷裂帶距離：研究區(qū)構(gòu)造形跡以南北向為主，其次為北西向與北北西向構(gòu)造，褶皺、斷裂構(gòu)造發(fā)育。斷裂構(gòu)造的交叉復(fù)合部位巖土體尤其破碎，產(chǎn)生的松散碎屑物質(zhì)為泥石流發(fā)育提供了物源儲備。距離斷裂構(gòu)造500～1 000 m處為研究區(qū)泥石流發(fā)育的優(yōu)勢區(qū)間。

(8)距主要道路距離：這一項指標(biāo)在一定程度上反映了人類活動作用。研究區(qū)主要道路靠近金沙江干流，海拔較低，屬于干熱河谷氣候。主要道路地帶植被不發(fā)育，巖土物理風(fēng)化嚴(yán)重，加上人類活動如修路切坡、伐木毀林等現(xiàn)象在這一帶集中，導(dǎo)致近公路段成為泥石流敏感區(qū)。

表1 指標(biāo)分級與重賦值

(續(xù))表1 指標(biāo)分級與重賦值

(9)地表粗糙度：反映地表微地形，即地表起伏變化與侵蝕程度。地表粗糙度對近地面風(fēng)速等風(fēng)沙活動造成影響，同時也反映了地表的抗風(fēng)蝕能力。該指標(biāo)被定義為地表單元曲面面積與投影面積之比。從表1分析結(jié)果來看，該指標(biāo)與泥石流發(fā)育的相關(guān)性不顯著，需要在后續(xù)建立模型過程中進一步觀察，再決定該指標(biāo)的保留與否。

(10)雨季月平均降雨量：降雨是泥石流主要觸發(fā)條件之一。研究區(qū)為干熱河谷，年均降雨量<400 mm，時空分布不均，雨季旱季分化明顯。每年6月至9月為雨季，降雨強度達到峰值，雨量集中。表1的分析結(jié)果顯示，隨降雨量增加，泥石流發(fā)育程度越高。

3.2 評價模型建立與ROC檢驗

對研究區(qū)的220個小流域樣本進行75%隨機抽樣，將得到的流域樣本數(shù)據(jù)組成訓(xùn)練組來建立LR模型；余下25%的流域樣本組成檢驗組，用來檢驗評價效果的可靠性。

訓(xùn)練組樣本數(shù)據(jù)如果不能覆蓋各二級分類指標(biāo)區(qū)間，意味著抽樣不全面，建模效果將會受到影響，因此有必要對抽樣建組結(jié)果進行復(fù)核。如圖3所示，經(jīng)隨機抽樣得到的全體樣本數(shù)據(jù)、訓(xùn)練組樣本數(shù)據(jù)、檢驗組樣本數(shù)據(jù)在各個二級指標(biāo)中的分布趨勢基本吻合，且各區(qū)間不存在抽樣空白，證明抽樣結(jié)果均衡有效。

將訓(xùn)練組樣本數(shù)據(jù)的各評價指標(biāo)的標(biāo)準(zhǔn)化指標(biāo)值Xi,j作為模型的回歸變量，泥石流事件Yi作為響應(yīng)變量，基于SPSS10.0軟件進行回歸分析。

在LR模型建立過程中，當(dāng)選取的評價指標(biāo)全部參與建模時可能導(dǎo)致評價結(jié)果不顯著。為了明確對研究區(qū)泥石流事件的發(fā)生貢獻突出的指標(biāo)，找到預(yù)測模型的最優(yōu)指標(biāo)組合，在選擇參與擬合的指標(biāo)時，本文采用了評價指標(biāo)體系的全部組合方式來進行建模，再利用受試者工作特征曲線(ROC曲線)一一檢測模型擬合效果。

如表2所示，加*的AUC值所對應(yīng)的指標(biāo)在變量數(shù)相同的模型組中，表現(xiàn)了最優(yōu)的建模擬合效果，因此，該指標(biāo)在模型組中被保留下來進而引入下一個變量。以此類推，當(dāng)模型中引入變量9時，AUC已達到最高值0.827，說明繼續(xù)加入第10個變量對模型的擬合效果并沒有良性影響。因此，選擇AUC為0.827時的模型作為最終評價模型，保留參與該建模過程中的9項指標(biāo)，并剔除坡向指標(biāo)F9。最優(yōu)模型的ROC檢驗結(jié)果如圖4所示。研究區(qū)泥石流易發(fā)性評價模型表達式如下：

Z=-10.396+4.957F1-2.278F2+4.860F3+
4.619F4+8.274F5+2.374F6+
7.687F7+16.349F8+6.975F10

(6)

表3和表4分別顯示了LR模型擬合統(tǒng)計量值和邏輯回歸系數(shù)及相關(guān)統(tǒng)計參數(shù)，在邏輯回歸系數(shù)B的分析結(jié)果中，除坡向F2、溝壑密度F4兩項指標(biāo)Sig值分別為0.065、0.055，略大于0.05，其他7項指標(biāo)Sig值均小于0.05，本模型基本通過了5%顯著性水平檢驗，可用于進一步評價。

表2 55組LR模型AUC值測試結(jié)果

注：在建模過程中引入相同數(shù)量的變量時，測試效果最佳的模型AUC值以“*”的形式來表示。

表3 LR模型擬合統(tǒng)計量

注：df為自由度；Sig為顯著性；-2LL為似然函數(shù)的自然對數(shù)的-2倍；Cox & SnellR2和NagelkerkeR2為偽決定系數(shù)。

表4 邏輯回歸系數(shù)及相關(guān)統(tǒng)計參數(shù)

注：B為模型中各個參數(shù)的回歸系數(shù)；S.E.(Std.Error)為標(biāo)準(zhǔn)誤差；Wals為Wald檢驗得到的統(tǒng)計量；df為自由度；Sig為顯著性。

3.3 評價結(jié)果分析與總體檢驗

通過ArcGIS10.0將建立的LR模型運用到指標(biāo)分級重賦值處理后的圖層組進行合成，得到研究區(qū)泥石流易發(fā)性概率圖。認(rèn)為易發(fā)性概率P大于0.5時，泥石流易發(fā)性等級為高或極高；概率小于0.1時，認(rèn)為易發(fā)性等級為極低?；诖?，將研究區(qū)劃分為5個等級區(qū)，即極低易發(fā)區(qū)(P=[0，0.1])，低易發(fā)區(qū)(P=[0.1，0.3])，中等易發(fā)區(qū)(P=[0.3，0.5])，高易發(fā)區(qū)(P=[0.5，0.7])，極高易發(fā)區(qū)(P=[0.7，1])。研究區(qū)泥石流易發(fā)性等級區(qū)劃圖、評價結(jié)果分區(qū)統(tǒng)計如圖5和表5所示。

評價結(jié)果顯示，研究區(qū)泥石流集中發(fā)育在昌波—貢波、徐龍—古學(xué)的金沙江干流沿岸一帶，局部高發(fā)區(qū)分布在瓦卡流域處。中等易發(fā)區(qū)與極低易發(fā)區(qū)的分布范圍分別占全區(qū)的29.78%、12.79%，同比其他3區(qū)分別處于較大和較小兩個極端，但差異并不懸殊，總體上看來，5個等級分區(qū)面積比較均衡。在極高易發(fā)區(qū)，實際發(fā)育泥石流面積465.46 km2，共計32個泥石流流域，泥石流災(zāi)害占該分區(qū)面積的75.95%，占研究區(qū)災(zāi)害總面積的37.15%；高易發(fā)區(qū)的泥石流面積349.46 km2，共發(fā)育泥石流流域26處，占分區(qū)面積的69.24%，占全區(qū)災(zāi)害總面積的27.89%；中等易發(fā)區(qū)發(fā)育泥石流面積334.07 km2，包含24個泥石流流域，占分區(qū)面積的36.12%，占全區(qū)災(zāi)害總面積的26.66%；低易發(fā)區(qū)和極低易發(fā)區(qū)共計存在泥石流流域9處，泥石流流域面積合占全區(qū)災(zāi)害總面積的8.30%。在極高易發(fā)區(qū)和高易發(fā)區(qū)中，已發(fā)泥石流流域共計58處(全區(qū)泥石流流域共91處)，占全區(qū)泥石流流域數(shù)量的63.7%，占全區(qū)災(zāi)害面積的65.0%；在極高易發(fā)區(qū)、高易發(fā)區(qū)以及中等易發(fā)區(qū)中，已發(fā)泥石流流域共計82處，占全區(qū)泥石流流域數(shù)量的90.1%，占全區(qū)災(zāi)害面積的91.7%。由此可見，高易發(fā)區(qū)與低易發(fā)區(qū)的分區(qū)界限清晰且災(zāi)害易發(fā)性區(qū)別顯著，泥石流主要發(fā)育在中等以及中等以上易發(fā)性分區(qū)內(nèi)。

圖3 全體樣本、訓(xùn)練組樣本、檢驗組樣本在各指標(biāo)中的分布頻率統(tǒng)計圖Fig.3 Histograms for the training and testing data subsets

從極高易發(fā)區(qū)至極低易發(fā)區(qū)，分區(qū)泥石流面積與分區(qū)總面積比值、分區(qū)災(zāi)害面積與全區(qū)災(zāi)害總面積比值的分布呈現(xiàn)不同幅度的下降趨勢，說明各分區(qū)之間的災(zāi)害易發(fā)性具有明顯的等級變化特征。從總體檢驗來看，泥石流易發(fā)性評價和分區(qū)結(jié)果與實際情況吻合較好。

3.4 樣本個案檢驗

在SPSS 25%隨機抽樣操作下的檢驗組樣本數(shù)據(jù)共計55個，包含24個泥石流流域和31個非泥石流流域。利用檢驗組樣本對本次易發(fā)性評價結(jié)果進行驗證與檢測，結(jié)果如表6和表7所示。表6顯示了訓(xùn)練組樣本總體預(yù)測精度檢驗結(jié)果，P為預(yù)測發(fā)生災(zāi)害概率，當(dāng)P>0.5即判斷為泥石流流域，得到的預(yù)測命中率為70.8%；判斷非泥石流流域時，在P值分別<0.3和<0.5的水平下，非泥石流流域判斷命中率分別為64.5%、81.6%，兩種情況下，樣本預(yù)測整體命中率為67.7%和76.2%。

圖4 泥石流易發(fā)性評價模型ROC曲線檢驗Fig.4 ROC curve of the LR-based debris-flow susceptibility assessment model

圖5 研究區(qū)泥石流易發(fā)性等級分區(qū)圖Fig.5 Debris-flow susceptibility map of the study area

表7顯示了檢驗組55個樣本流域在各區(qū)的分布，對比該樣本中實際泥石流流域分布情況可見，本次LR模型在各分區(qū)的抽樣調(diào)查中也有較好的表現(xiàn)。隨著易發(fā)性等級降低，檢驗組樣本中實發(fā)泥石流流域數(shù)量顯著減低。個體預(yù)測準(zhǔn)確性方面，極高易發(fā)區(qū)個體樣本預(yù)測一致率為91.7%；其次為高易發(fā)區(qū)，達到75.0%。二者較高的預(yù)測正確率也說明，泥石流集中發(fā)生于極高易發(fā)區(qū)與高易發(fā)區(qū)。評價結(jié)果經(jīng)個案檢驗證明預(yù)測效果良好。

表5泥石流易發(fā)性評價結(jié)果統(tǒng)計

Table5Statisticsofdebris-flowsusceptibilityassessmentresults

易發(fā)性分區(qū)分區(qū)面積/km2(分區(qū)/全區(qū))/%泥石流面積/km2(泥石流面積/分區(qū)面積)/%(分區(qū)災(zāi)害面積/全區(qū)災(zāi)害面積)/%極高612.8919.73465.4675.9537.15高504.6716.25349.4669.2427.89中等924.9929.78334.0736.1226.66低665.6921.44101.2115.208.08極低397.3312.792.760.690.22

表6檢驗組樣本總體預(yù)測精度檢驗結(jié)果

Table6Resultsoftheoverallforecastingaccuracyforthetestingdatasets

流域預(yù)測情況/個實際情況/個命中率/%整體命中率/%泥石流流域P>0.5(極高、高易發(fā)區(qū))172038—非泥石流流域P<0.3(極低、低易發(fā)區(qū))24313167.7非泥石流流域P<0.570.864.581.676.2

表7 檢驗組樣本分區(qū)分布情況與精度

4 討論與結(jié)論

本文對金沙江上游干熱河谷區(qū)的奔子欄—昌波河段進行了泥石流的易發(fā)性研究，并對LR評價模型的建模過程進行優(yōu)化，強化了評價指標(biāo)的篩選過程，在評價結(jié)果的檢驗方法方面做出改進，提高了研究區(qū)泥石流易發(fā)性區(qū)劃的準(zhǔn)確性與評價結(jié)果的可靠度，取得主要認(rèn)識如下：

(1)在指標(biāo)分析篩選和LR建模方面，初步選取10項，在指標(biāo)敏感性分析中標(biāo)定了敏感度較低的指標(biāo)F9。在建立LR模型時關(guān)注了指標(biāo)篩選與組合，使建立的模型最大化地適用于本區(qū)干熱河谷，測得模型預(yù)測成功率為0.827(AUC)；同時，在建模過程中再次驗證了引入指標(biāo)F9對模型意義不大，因此將指標(biāo)“地表粗糙度”一項剔除。通過55組模型的預(yù)測效果對比發(fā)現(xiàn)，同一指標(biāo)在獨立參與評價、與其他指標(biāo)組合參與評價的兩種情況下，發(fā)揮作用的權(quán)重不盡相同，從而體現(xiàn)了指標(biāo)評價體系中各致災(zāi)因素在泥石流發(fā)育過程中的協(xié)同與牽制作用關(guān)系。

(2)研究區(qū)泥石流易發(fā)性評價結(jié)果表明了本區(qū)致災(zāi)主導(dǎo)因素的指標(biāo)排序為：距主干公路距離>巖性>距斷裂帶距離>雨季月平均降雨量，坡度、坡向、流域相對高差、溝壑密度、植被歸一化指數(shù)5項指標(biāo)發(fā)揮輔助作用。研究區(qū)的主干道路比鄰金沙江干流，切坡、修路、砍樹伐林等人類活動比較強烈，導(dǎo)致干流兩岸的崩塌滑坡以及傾倒變形體較發(fā)育，形成了泥石流重要物質(zhì)來源。由此可判斷，在季節(jié)性暴雨提供動力條件的影響下，強烈的人類活動因素促進了本區(qū)干熱河谷泥石流的發(fā)育。

(3)對評價結(jié)果的可靠性先后進行了易發(fā)性區(qū)劃的總體檢驗與流域樣本的隨機個案檢驗?？傮w檢驗表明，以流域數(shù)量和流域面積為單位統(tǒng)計，極高易發(fā)區(qū)和高易發(fā)區(qū)兩處的實發(fā)泥石流與全區(qū)泥石流比例之和，分別為63.7%和65.0%；在極高易發(fā)區(qū)、高易發(fā)區(qū)以及中等易發(fā)區(qū)三區(qū)中，比例分別為90.1%和91.7%，證明了各分區(qū)之間的泥石流易發(fā)性具有明確的等級變化特點，且泥石流易發(fā)區(qū)與低發(fā)區(qū)區(qū)別顯著，易發(fā)性分區(qū)結(jié)果具有實際意義。隨機個案檢驗結(jié)果表明，檢驗組小流域樣本中，被判為泥石流易發(fā)流域的樣本與實發(fā)泥石流流域分布一致的比率隨分區(qū)易發(fā)性等級降低依次為91.7%、75.0%、36.4%、16.7%、0，從流域個案角度證明了各分區(qū)泥石流易發(fā)性呈現(xiàn)逐級減弱的特征?？傮w檢驗與個案檢驗均證明了泥石流易發(fā)性評價結(jié)果可靠性良好。

參考文獻：

[1] 唐川, 朱靜. 云南滑坡、泥石流研究[M]. 北京:商務(wù)印書館, 2003: 76-115.

[2] 謝洪, 鐘敦倫, 韋方強, 等. 我國山區(qū)城鎮(zhèn)泥石流災(zāi)害及其成因[J].山地學(xué)報, 2006,24(1):79-87.

[3] 羅元華. 論泥石流災(zāi)害風(fēng)險評估方法[J]. 中國礦業(yè), 2000, 9 (6): 70-72.

[4] 李闊, 唐川. 泥石流危險性評價研究進展[J]. 災(zāi)害學(xué)報, 2007, 22(1): 106-110.

[5] 叢威青, 潘懋, 李鐵峰, 等. 基于GIS的滑坡、泥石流災(zāi)害危險性區(qū)劃關(guān)鍵問題研究[J]. 地學(xué)前緣, 2006,13(2):185-190.

[6] 羅元華, 陳崇希. 泥石流堆積過程數(shù)值模擬及防災(zāi)效益評估方法[J]. 現(xiàn)代地質(zhì), 2000,14(4): 484-488.

[7] 鐘敦倫, 韋方強, 謝洪. 長江上游泥石流危險度區(qū)劃的原則與指標(biāo)[J]. 山地研究, 1994, 12(2):78-83.

[8] 劉希林, 陳宜娟. 泥石流風(fēng)險區(qū)劃方法及其應(yīng)用:以四川西部地區(qū)為例[J]. 地理科學(xué), 2010,30(4): 558-565.

[9] 崔鵬, 鄒強. 山洪泥石流風(fēng)險評估與風(fēng)險管理理論與方法[J]. 地理科學(xué)進展, 2016, 35(2):137-147.

[10] 王維早, 雷霆, 許強. 基于層次分析法的河北省太行山區(qū)泥石流災(zāi)害危險性評估研究[J]. 地球與環(huán)境, 2010, 38(3): 357-362.

[11] 寧娜, 馬金珠, 張鵬, 等. 基于GIS和信息量法的甘肅南部白龍江流域泥石流災(zāi)害危險性評價[J]. 資源科學(xué), 2013,35(4):892-899.

[12] 陳曉清, 謝洪. 基于GIS的泥石流危險度區(qū)劃研究—以擬建向家壩、溪洛鍍水電工程庫區(qū)為例[J]．土壤侵蝕與水土保持學(xué)報, 1999,5(6):46-50.

[13] 柳金峰, 黃江成, 歐國強, 等.甘肅隴南武都區(qū)泥石流易發(fā)性評價[J]. 中國地質(zhì)災(zāi)害與防治學(xué)報, 2010,21(4):8-12.

[14] 李為樂, 唐川, 楊武年, 等. RS和GIS技術(shù)在區(qū)域泥石流危險性區(qū)劃中的應(yīng)用研究[J]. 災(zāi)害學(xué), 2008,23(2):71-75.

[15] 唐堯, 楊武年. 基于GIS的震后汶川潛在泥石流危險性評價[J]. 上海國土資源, 2012,33(3):57-60.

[16] 吳雨夫,余斌, 亓星, 等. 岷江上游大白楊溝泥石流特征及成因[J]. 現(xiàn)代地質(zhì), 2012,26(1):107-113.

[17] 李雅暉. 基于流域尺度的泥石流危險性評價:以岷江流域汶川縣為例[D]. 成都: 成都理工大學(xué),2011:21-45.

[18] 戴福初, 李焯芬, 黃志全, 等. 火山巖坡殘積土地區(qū)暴雨滑坡泥石流的形成機理[J]. 工程地質(zhì)學(xué)報, 1999,7(2):147-154.

[19] SONG S J, ZHANG B L, FENG W L, et al. Using fuzzy relations and GIS method to evaluate debris flow hazard[J]. Wuhan University Journal of Natural Sciences, 2006, 11(4): 875-881.

[20] YUAN L F, ZHANG Y S. Debris flow hazard assessment based on support vector machine[J]. Wuhan University Journal of Na-tural Sciences, 2006, 11(4): 897-900.

[21] GORSEVSKI P V,GESSLER P E,FOLTZ R B. Spatial prediction of landslide hazard using logistic regression and GIS[J].Transactions in GIS, 2000,10 (3):395-415.

[22] 邢秋菊, 趙純勇, 高克昌. 基于GIS的滑坡危險性邏輯回歸評價研究[J]. 地理與地理信息科學(xué), 2004,20(3):49-51.

[23] 劉藝梁, 殷坤龍, 劉斌. 邏輯回歸和人工神經(jīng)網(wǎng)絡(luò)模型在滑坡災(zāi)害空間預(yù)測中的應(yīng)用[J]. 水文地質(zhì)工程地質(zhì)，2010, 37(5):92-96.

[24] 許沖, 徐錫偉. 邏輯回歸模型在玉樹地震滑坡危險性評價中的應(yīng)用與檢驗[J]. 工程地質(zhì)學(xué)報, 2012, 20(3): 326-333.

[25] 李雪平, 唐輝明. 基于GIS的分組數(shù)據(jù)Logistic模型在斜坡穩(wěn)定性評價中的應(yīng)用[J]. 吉林大學(xué)學(xué)報, 2005,35(3):361-365.

[26] OHLMACHER G C, DAVIS J C. Using multiple logistic regression and GIS technology to predict landslide hazard in northeast Kansas, USA[J]. Engineering Geology,2003,69:331-343.

[27] 王志旺, 廖勇龍, 李端有. 基于邏輯回歸法的滑坡危險度區(qū)劃研究[J]. 地下空間與工程學(xué)報,2006,2(2): 1451-1454.

[28] 潘赟, 叢威青, 潘懋. 基于GIS的遼寧省岫巖縣泥石流災(zāi)害危險性區(qū)劃研究[J]. 北京大學(xué)學(xué)報,2010,46(4): 601-606.

[29] NETRA R R, JOHN R G, ERIC V M, et al. A comparison of logistic regression-based models of susceptibility to landslides in western Colorado, USA[J]. Landslides,2014, 11:247-262.

[30] 鄒強,王青,劉延國. 基于GIS與Logistic模型的公路泥石流易發(fā)性分析[J]. 水土保持通報, 2014,34(3):185-189.

[31] 戴福初. 金沙江-奔子欄水電站庫區(qū)地質(zhì)構(gòu)造、環(huán)境地質(zhì)遙感研究研究報告[R].北京：中國科學(xué)院地質(zhì)與地球物理研究所，2005.

[32] 陳劍, 崔之久, 戴福初, 等. 金沙江奔子欄-達日河段大型泥石流堆積扇的成因機制[J]. 山地學(xué)報, 2011, 29(3) : 312-319.

[33] 姚鑫, 戴福初, 陳劍. 金沙江干熱河谷區(qū)地質(zhì)災(zāi)害遙感研究[J]. 長江流域資源與環(huán)境,2007,16(5)：655-660.

[34] CHEN J, DAI F C, YAO X. Holocene debris-flow deposits and their implications on climate in the upper Jinsha River valley,China[J]. Geomorphology, 2008, 93: 493-500.

[35] 黎艷,陳劍,許沖. 基于AHP的半干旱區(qū)泥石流易發(fā)性評價：以金沙江上游奔子欄—昌波河段為例[J]. 現(xiàn)代地質(zhì), 2015, 29(4):975-982.

[36] CHEN Jian, LI Yan, ZHOU Wendy, et al. Debris-flow susceptibility assessment model and its application in semiarid mountainous areas of the southeastern semiarid mountainous areas of the southeastern Tibetan Plateau[J]. Natural Hazards Review, 2016.doi:10.1061/ (ASCE) NH. 1527-6996.0000229.