基于聚類和最大似然法的汶川災區(qū)泥石流滑坡易發(fā)性評價

胡凱衡，崔鵬，韓用順，游勇

(中國科學院山地災害與地表過程重點實驗室;中國科學院水利部成都山地災害與環(huán)境研究所:610041，成都)

2008 年5 月12 日，四川省龍門山地區(qū)發(fā)生了Ms8.0 的強烈地震，震中位于汶川縣映秀鎮(zhèn)，距成都約80 km(圖1)。地震引發(fā)了大量的滑坡、崩塌、滾石、堰塞湖［1-2］，這些同震的地質災害為震后泥石流提供了豐富的松散物源。居民點和基礎設施，如公路、橋梁和水電站都受到了嚴重威脅。震后的地質災害造成了更多的人員傷亡和財產損失，如2008 年9 月24 日，北川縣一場24 h 降雨量為272.4 mm 的暴雨激發(fā)了大面積的泥石流，僅北川老縣城附近就有42 人死亡，幾十間安置房被摧毀。Cui 等［1］認為，震后地質災害的活躍期最長可以持續(xù)20 年左右。為了指導震后災害評估、搬遷安置和災害防治，有必要對震區(qū)開展災害易發(fā)性評價工作。

災害的易發(fā)性可以看成某一地學單元發(fā)生災害的空間可能性，通常被描述為一些影響因子的函數(shù)，如地層、斷裂、坡度、土地利用等［3-6］。目前，已經發(fā)展了許多的易發(fā)性評估方法［3-14］，He 等［14］把這些方法大致劃分為定性/半定量方法和定量方法。定量方法包括統(tǒng)計和機制性方法。后驗統(tǒng)計方法，如判別分析和邏輯回歸方法，可以從災害分布數(shù)據建立易發(fā)性和影響因子之間的統(tǒng)計關系;但如果沒有足夠的災害分布數(shù)據，就無法應用這些方法，而地震發(fā)生后短時間內無法了解災害的分布情況，因此，急需研究在缺少或者只知道部分災害分布的情況下對地質災害的易發(fā)性進行快速評估的方法。

為此，筆者擬采用聚類分析和最大似然法對強震后泥石流滑坡的易發(fā)性進行快速評估。聚類分析法是多元統(tǒng)計分析中的一種常用非監(jiān)督分類法，在數(shù)據挖掘、模式識別、圖像分析和機器學習等研究中得到了廣泛應用。與后驗方法相比，非監(jiān)督分類方法有2 個優(yōu)點:可以不用任何訓練數(shù)據將一個多元數(shù)據集分成幾個子類;可以避免在賦予影響因子的權重或分類閾值時的主觀不確定性。將這一方法應用于汶川地震極重災區(qū)泥石流滑坡的易發(fā)性評價中具有其他方法不可比擬的快速評價的優(yōu)勢。

1 研究區(qū)概況

研究區(qū)位于青藏高原東緣到四川盆地的過渡帶龍門山地區(qū)，是受此次地震影響最大的區(qū)域。龍門山位于四川盆地西北緣，長約500 km，寬約30 ～50 km，呈北東—南西方向延伸，西北與青藏高原相鄰，東南與四川盆地相接。地勢陡峻，一般高度在3 km以上，最高處超過4 km，而其東南的四川盆地平均高度約500 m。龍門山高高地突起在四川盆地西緣，構成了中國西部與東南部的主要地形分界線。這一區(qū)域主要有3 條大的斷層，分別是龍門山后山斷裂、中央斷裂和前山斷裂。后山和中央斷裂之間主要由志留系—泥盆系淺變質巖和前寒武系雜巖構成，中央和前山斷裂之間則主要由上古生界—三疊系沉積巖構成，包括古生代灰?guī)r和含煤巖系［15］。在地貌上，這一區(qū)域也可劃分為3 個單元:中央斷裂以西海拔高于2 000 m 的高山區(qū);中央和前山斷裂之間海拔在1 000 ～2 000 m 之間的丘陵區(qū);山前斷裂帶以東南的平原地貌區(qū)(海拔400 ～700 m)。這一區(qū)域的氣候在夏季主要為太平洋高壓帶來的亞熱帶暖濕季風控制，來自東南的暖濕氣流受龍門山脈的阻擋，使當?shù)爻蔀榇ㄎ鞯貐^(qū)的暴雨中心區(qū)，山脈東側的年平均降水1 100 ～1 500 mm。據l954—1987 年綿竹漢旺氣象站資料，年均降水量1 503 mm，最大年降水量2 082 mm(1961 年)，多年平均最大1 d 降雨量149 mm，并多次出現(xiàn)大暴雨或特大暴雨，而北川縣城觀測到的最大年降水量為2 340 mm。山脈西側主要由暖溫帶半干旱的大陸季風控制，年均降水量在500 mm 左右。這些自然條件使得這一區(qū)域成為中國最容易發(fā)生地質災害的區(qū)域之一。研究區(qū)域為汶川地震極重災區(qū)四川省，總面積2 萬5 800 km2，包括汶川縣、北川縣、綿竹市、什邡市、青川縣、茂縣、安縣、都江堰市、平武縣、彭州市10 個縣市(圖1)，為國家災后恢復重建的重點規(guī)劃區(qū)域。

圖1 研究區(qū)和龍門山斷裂Fig.1 Location of study area and Longmenshan faults

2 研究方法

2.1 聚類分析和最大似然法

災害易發(fā)性的聚類分析方法是基于這樣的假設:易發(fā)性相同的空間單元其影響因子的值也應該比較接近，假設影響災害易發(fā)性的因子有n 個，每個因子在空間上的分布值離散后就是一個柵格，n 個因子的值就組成一個柵格列(圖2)，在同一個空間單元上就有n 個值，構成一個向量(X1，X2，…，Xn)，而一個向量可以看成n 維屬性空間的一個屬性點，這樣，n 個柵格就對應屬性空間的一個點集。按照這個假設，易發(fā)性比較接近的點應該集中在某個范圍內，或者說一個子集里。聚類分析的目的就是將屬性空間的一個點集劃分為不同的子類或子集，并計算每個子類的統(tǒng)計特征值，如均值和方差。在實際操作時，不是所有的空間單元都參與到子類劃分中，只需要選擇研究區(qū)域一個樣板區(qū)來進行聚類分析，獲取不同子類的特征值。進行聚類分析前只需要給出子類的個數(shù)。

圖2 易發(fā)性聚類分析示意圖(假設影響因子有3 個)Fig.2 Schematics of basic idea of susceptibility mapping with cluster analysis(supposed that there are three effective factors X1，X2，X3，n=3)

聚類分析的算法采用的是k-平均聚類法。計算步驟如下:1)假設k 個子類的初始質心(即均值)均勻或者隨機分布在屬性空間(圖3(a));2)計算點集的每一個點與k 個子類質心的距離，將每一個點歸到與其距離最小的那個子類(圖3(b));3)計算第2 步得到的新子類的質心(圖3(c));4)重復第1步直到每個子類的點數(shù)變動少于2%。

圖3 k-平均聚類算法示意圖Fig.3 Schematics of k-mean clustering algorithm

聚類分析完成后，可以采用最大似然法MLC(maximum likelihood classification)判斷研究區(qū)域的任意一個空間單元(網格)屬于哪一個子類。根據每個子類的均值和協(xié)方差矩陣，基于Bayesian 判別準則和多元屬性空間高斯分布的假設，MLC 方法能計算每一個單元屬于不同子類的概率。MLC 方法中假設每個子類的先驗概率都是相同的，當然，也可以賦給每個子類不同的先驗概率(圖4)。對某一空間單元來說，它所屬的子類就是概率最大的那個子類，落入同一個子類的那些空間單元可以認為易發(fā)性相等。

2.2 影響因子

地質災害發(fā)生的3 個充分條件是勢能、水源和物源。根據這3 個條件，震后易發(fā)性評價的影響因子選擇了坡度、相對高差、地表徑流深和地震烈度。坡度反映了松散物質起動的動力條件，坡度越陡，物質就越容易起動，但僅用坡度一個因子不能完全反映勢能的影響，如泥石流溝口的堆積扇區(qū)域比較平坦，坡度比較小，而泥石流堆積扇卻是最容易遭受泥石流危害的地方;因此，相對高差也是比較重要的。將相對高差定義為

圖4 易發(fā)性評價流程圖Fig.4 Flow diagram of susceptibility mapping by the cluster method

式中:R 為相對高差，m;Hmax、Hmin分別為面積等于ΔA 的一個空間鄰域內單元高程的最大值和最小值，m。地表徑流深可以反映降水和土壤類型等，是一個綜合因子。地震影響著震后的地質災害的坡體破壞程度或物質來源。如果純粹考慮因地震而產生的地質災害，可以認為坡體破壞程度和松散物質量與地震能量，即地震烈度大致呈正比關系;所以，用地震烈度來表示物質這一因素。

2.3 數(shù)據和分類

坡度和相對高差柵格數(shù)據來自SRTM(Shuttle Radar Topography Mission)提供的DEM 數(shù)據，空間分辨率為74 m×74 m(圖5(a)和(b))。坡度用ArcGIS 中的表面分析工具計算，相對高差用ArcGIS的鄰域統(tǒng)計函數(shù)計算。鄰域的搜索半徑設為10 個網格單元，即ΔA≈0.74×0.74π≈1.72 km2。地震烈度柵格圖層來自于美國地調局網站上公布的MMI(Modified Mercalli Intensity)矢量數(shù)據(圖5(c))。地表徑流深柵格數(shù)據來自四川省水文局制作的四川省1956—2000 年平均年徑流深等值線圖(圖5(d))。

因為聚類分析中采用的距離是歐幾里得距離，所以，聚類的結果對柵格圖層的取值范圍比較敏感。為了消除這一影響，每個柵格圖層都用它的最大值歸一化，使其值在0 到1.0 之間，如地震烈度的最大值為9.9，烈度柵格層除以9.9 后，取值在0.545 到1.0 之間。歸一化后，采用k-平均聚類方法將柵格集分為極低、低、中等、高、極高5 個子類，并輸出每個子類的統(tǒng)計特征值(表1)。聚類分析的樣本集為10×10 的子塊中取一個單元，迭代到100 次之后5個子類的均值和方差都沒有變化為止。最大似然估計用這些子類的最終統(tǒng)計特征值來判斷每個單元屬于哪一個子類。所有的這些操作，包括數(shù)據轉換、柵格計算和聚類計算都在ArcGIS 中實現(xiàn)。

表1 k-平均聚類法得到的5 個子類的歸一化均值Tab.1 Means of five clusters obtained by k-mean clustering

3 結果與分析

采用聚類和MLC 方法最終得到5 種級別的泥石流滑坡易發(fā)性評價結果，但是，為了知道哪個子類對應哪種易發(fā)性，還需要進一步做出判斷。基于一些常識或專家的經驗，很容易判斷出第1 個子類為極高易發(fā)性區(qū)，第2 個子類為高易發(fā)區(qū)，第5 子類為極低易發(fā)區(qū);但是，不太容易判斷第3 和第4 個子類哪個是中等易發(fā)區(qū)，哪個是低易發(fā)區(qū)。為此，用地質災害的點密度(個/km2)來判別這2 個子類。一般來說，災害的點密度值越大，易發(fā)性也就越高，災害點密度的數(shù)值來自從遙感圖像和應急調查中獲取的有限地質災害分布數(shù)據。這個分布數(shù)據只有3 967 個災害點，是地震發(fā)生后短時間內獲得的。計算結果(表2)表明，第3 個子類的點密度為0.123，而第4 個子類的點密度為0.082。由此可以判斷出第3 個子類為中等易發(fā)區(qū)，第4 個子類為低易發(fā)區(qū)。另外，第1、2 和5 子類的點密度分別為0.391、0.198和0.071，說明對這3 個子類的易發(fā)性分析是正確的。最終的易發(fā)性分區(qū)圖見圖6。

圖5 影響因子的柵格圖Fig.5 Grid maps of four effective factors

表2 各子類的面積和災害的點密度Tab.2 Area and geo-hazard density of five clusters

從表1 可以觀察到，極高子類的地震烈度、地表徑流深度值都是最大，而它的坡度和相對高差并不是最大。與此相反，極低子類具有最小的坡度和最小的相對高差值。從圖6 可以發(fā)現(xiàn)，即使在高易發(fā)區(qū)域(圖6 中的粉紅色區(qū)域)，放大之后可以看到有極低易發(fā)區(qū)域(綠色區(qū)域)存在，如茂縣、青川和平武縣縣城附近，這幾個縣城沒有直接受到震后地質災害的威脅，泥石流滑坡的危害也不明顯;因此，這一分區(qū)圖可以用來尋找危險區(qū)中的安全島(即小塊的低易發(fā)區(qū))。從分區(qū)圖中還可以發(fā)現(xiàn)一點，在極高易發(fā)區(qū)和高易發(fā)區(qū)以及高易發(fā)區(qū)與極低易發(fā)區(qū)分界線附近，地質災害點分布比較密集，如北川老縣城就位于極高易發(fā)區(qū)和高易發(fā)區(qū)交界處，就在地震后不到1 年的時間里，滑坡和泥石流就毀壞了大部分的縣城(圖7)。

圖6 地質災害分布和易發(fā)性評價圖Fig.6 Maps of geo-hazards and susceptibility

圖7 2009 年4 月16 日的北川縣城遺址Fig.7 City of Beichuan on 16 April，2009 inundated by post-quake rainstorm-triggered debris flows

4 結論與討論

在沒有足夠訓練數(shù)據的情況下，采用聚類分析和最大似然方法可以快速對一個大區(qū)域的地質災害易發(fā)性進行評估，且不需要給定各個影響因子的權重。

研究選取了坡度、相對高差、地震烈度和地表徑流深4 個影響因子，雖然沒有考慮更多的因子，如斷裂距離、巖性等，但最終得到的地質災害易發(fā)性評價圖比較合理，與地質災害的點密度值吻合。從易發(fā)性評價圖可以發(fā)現(xiàn)，在大塊的高易發(fā)區(qū)內存在一些小的極低易發(fā)區(qū)，這些極低易發(fā)區(qū)可以作為安全的居民點。另外，在極高易發(fā)區(qū)和高易發(fā)區(qū)以及高易發(fā)區(qū)和極低易發(fā)區(qū)分界線附近，地質災害點分布比較密集。

［1］ Cui P，Chen X Q，Zhu Y Y，et al.The Wenchuan Earthquake(May 12，2008)，Sichuan Province，China，and resulting geohazards［J］.Nat Hazards，2011，56(1):19-36

［2］ 黃潤秋，李為樂.“5·12”汶川大地震觸發(fā)地質災害的發(fā)育分布規(guī)律研究［J］.巖石力學與工程學報，2008，27(12):2585-2592

［3］ Dai Fuchu，Lee C F.Landslide characteristics and slope instability modeling using GIS，Lantau Island，Hong Kong［J］.Geomorphology，2002，42:213-228

［4］ Lin M L，Tung C C.A GIS-based potential analysis of the landslides induced by the Chi-Chi earthquake［J］.Eng Geol，2004，71:63-77

［5］ Kamp U，Growley B J，Khattak G A，et al.GIS-based landslide susceptibility mapping for the 2005 Kashmir earthquake region［J］.Geomorphology，2008，101: 631-642

［6］ Carrara A，Crosta G，F(xiàn)rattini P.Comparing models of debris-flow susceptibility in the alpine environment［J］.Geomorphology，2008，94:353-378

［7］ Carrara A，Cardinali M，Detti R，et al.Gis techniques and statistical-models in evaluating landslide hazard［J］.Earth Surf Proc Land，1991，16:427-445

［8］ Jade S，Sarkar S.Statistical-models for slope instability classification［J］.Eng Geol，1993，36:91-98

［9］ Jiang H，Eastman J R.Application of fuzzy measures in multi-criteria evaluation in GIS［J］.Int J Geogr Inf Sci，2000，14:173-184

［10］Van Westen C J，Rengers N，Soeters R.Use of geomorphological information in indirect landslide susceptibility assessment［J］.Nat Hazards，2003，30:399-419

［11］Lee S，Ryu J H，Min K D，et al.Landslide susceptibility analysis using GIS and artificial neural network［J］.Earth Surf Proc Land，2003，28:1361-1376

［12］Lan H X，Zhou C H，Wang L J，et al.Landslide hazard spatial analysis and prediction using GIS in the Xiaojiang watershed，Yunnan，China［J］.Eng Geol，2004，76:109-128

［13］Ercanoglu M.Landslide susceptibility assessment of SE Bartin(West Black Sea region，Turkey)by artificial neural networks［J］.Nat Hazard Earth Sys，2005，5: 979-992

［14］He Y P，Beighley R E.GIS-based regional landslide susceptibility mapping: a case study in southern California［J］.Earth Surf Proc Land，2008，33:380-393

［15］李勇，黃潤秋，周榮軍，等.龍門山地震帶的地質背景與汶川地震的地表破裂［J］.工程地質學報，2009，17(1):3-18