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

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

(中國科學(xué)院山地災(zāi)害與地表過程重點(diǎn)實(shí)驗(yàn)室;中國科學(xué)院水利部成都山地災(zāi)害與環(huán)境研究所:610041，成都)

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

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

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

1 研究區(qū)概況

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

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

2 研究方法

2.1 聚類分析和最大似然法

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

圖2 易發(fā)性聚類分析示意圖(假設(shè)影響因子有3 個(gè))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-平均聚類法。計(jì)算步驟如下:1)假設(shè)k 個(gè)子類的初始質(zhì)心(即均值)均勻或者隨機(jī)分布在屬性空間(圖3(a));2)計(jì)算點(diǎn)集的每一個(gè)點(diǎn)與k 個(gè)子類質(zhì)心的距離，將每一個(gè)點(diǎn)歸到與其距離最小的那個(gè)子類(圖3(b));3)計(jì)算第2 步得到的新子類的質(zhì)心(圖3(c));4)重復(fù)第1步直到每個(gè)子類的點(diǎn)數(shù)變動(dòng)少于2%。

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

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

2.2 影響因子

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

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

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

2.3 數(shù)據(jù)和分類

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

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

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

3 結(jié)果與分析

采用聚類和MLC 方法最終得到5 種級(jí)別的泥石流滑坡易發(fā)性評(píng)價(jià)結(jié)果，但是，為了知道哪個(gè)子類對(duì)應(yīng)哪種易發(fā)性，還需要進(jìn)一步做出判斷?；谝恍┏ＷR(shí)或?qū)＜业慕?jīng)驗(yàn)，很容易判斷出第1 個(gè)子類為極高易發(fā)性區(qū)，第2 個(gè)子類為高易發(fā)區(qū)，第5 子類為極低易發(fā)區(qū);但是，不太容易判斷第3 和第4 個(gè)子類哪個(gè)是中等易發(fā)區(qū)，哪個(gè)是低易發(fā)區(qū)。為此，用地質(zhì)災(zāi)害的點(diǎn)密度(個(gè)/km2)來判別這2 個(gè)子類。一般來說，災(zāi)害的點(diǎn)密度值越大，易發(fā)性也就越高，災(zāi)害點(diǎn)密度的數(shù)值來自從遙感圖像和應(yīng)急調(diào)查中獲取的有限地質(zhì)災(zāi)害分布數(shù)據(jù)。這個(gè)分布數(shù)據(jù)只有3 967 個(gè)災(zāi)害點(diǎn)，是地震發(fā)生后短時(shí)間內(nèi)獲得的。計(jì)算結(jié)果(表2)表明，第3 個(gè)子類的點(diǎn)密度為0.123，而第4 個(gè)子類的點(diǎn)密度為0.082。由此可以判斷出第3 個(gè)子類為中等易發(fā)區(qū)，第4 個(gè)子類為低易發(fā)區(qū)。另外，第1、2 和5 子類的點(diǎn)密度分別為0.391、0.198和0.071，說明對(duì)這3 個(gè)子類的易發(fā)性分析是正確的。最終的易發(fā)性分區(qū)圖見圖6。

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

表2 各子類的面積和災(zāi)害的點(diǎn)密度Tab.2 Area and geo-hazard density of five clusters

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

圖6 地質(zhì)災(zāi)害分布和易發(fā)性評(píng)價(jià)圖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 結(jié)論與討論

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

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

［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ā)地質(zhì)災(zāi)害的發(fā)育分布規(guī)律研究［J］.巖石力學(xué)與工程學(xué)報(bào)，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］李勇，黃潤秋，周榮軍，等.龍門山地震帶的地質(zhì)背景與汶川地震的地表破裂［J］.工程地質(zhì)學(xué)報(bào)，2009，17(1):3-18