基于证据权方法的玉树地震滑坡危险性评价

2013-12-14徐锡伟于贵华

地震地质 2013年1期

许冲徐锡伟于贵华

(中国地震局地质研究所，活动构造与火山重点实验室，北京 100029)

0 引言

2010年4月14日，青海省玉树藏族自治州玉树县发生了MS7.1地震。震中位置为33.224°N，96.666°E，震源深度为 17km，震级 MW6.9(U.S.Geological Survey，2010)。据青海省抗震救灾指挥部通报，截至25日17时，青海玉树7.1级地震已造成2 220人死亡，70人失踪，12 135人受伤，其中重伤1 434人。地震触发滑坡对灾区也产生了一定的危害①(殷跃平等，2010a，b;Zhang et al.，2010a)。作者前期主要基于野外考察与遥感影像解译方法，对玉树地震滑坡开展了空间分布调查，并制作了玉树地震滑坡空间分布图。本文拟以玉树地震滑坡空间分布数据为基础，基于GIS平台，采用加法与减法2类证据权方法，考虑影响玉树地震滑坡发生的12个因素，对玉树地震灾区在类似地震条件下可能遭受的滑坡危险性进行评价，为灾区滑坡防灾减灾提供参考。

地震滑坡危险性评价研究是地震事件发生后，关于地震滑坡的一个重要研究方向(许冲等，2009，2010a，b，2011b;许冲等，2012a，b;Xu et al.，2012a，b，c，d)，它需要翔实可靠的滑坡空间分布数据与多种滑坡控制因素专题数据。作者的前期工作(许冲等，2011a，2012;Xu et al.，2012e)为开展本文的地震滑坡危险性评价研究奠定了良好的基础。证据权方法是一种较常见的滑坡危险性评价方法，该方法在国内外得到了较多的应用(van Westen et al.，2003;Lee et al.，2004;Thiery et al.，2004，2007;Zahiri et al.，2006;Neuhauser et al.，2007;Mathew et al.，2007;王志旺等，2007;Dahal et al.，2008a，b;Sharma et al.，2008;Song et al.，2008;van den Eeckhaut et al.，2009;许冲等，2011b;Xu et al.，2012d)。所以，本文基于GIS技术与2类证据权方法，开展玉树地震滑坡危险性评价研究工作。

1 玉树地震及诱发滑坡

青海玉树MS7.1地震发生在印度板块向欧亚大陆挤入形成喜马拉雅碰撞带以北的青藏高原。由于印度板块向青藏高原的推挤作用而导致青藏高原缩短，同时内部的块体沿一些重要的块体边界断裂带发生侧向滑移，造成青藏高原主体向东移动，并在青藏高原内部和块体边界形成不同规模的走滑断裂系与挤出块体。在这种构造背景下，巴颜喀拉(松潘－甘孜)地块向SE方向运动，块体南部的金沙江－鲜水河断裂带全新世以来平均水平滑动速率达5～7mm/a;北部的喀喇昆仑－嘉黎断裂带平均水平滑动速率达10～11mm/a。这2条近EW向的断裂带依次向东推移，使其东部的川、滇地区由于受扬子地块的阻挡而产生复杂的断裂形变与块体转动。玉树地震就是在这种地质背景下发生的，地震发生在鲜水河断裂带西段的甘孜－玉树断裂带上，是甘孜－玉树断裂带突发错动的结果。甘孜－玉树断裂带规模巨大，西起青海治多县那王草曲塘，经当江、玉树、邓柯、玉隆，至四川甘孜县城南，全长约500km。地震产生了走向约300°，65km长的地表破裂带，左旋走滑性质，地表破裂带由一系列挤压鼓包与张裂缝相间排列或雁列式裂缝组成，实测最大水平位错约1.8m(闻学泽等，2003;Xu et al.，2010;Pan et al.，2010;Wu et al.，2010;Bartholomew et al.，2010;Zhang et al.，2010b;顾瑾萍，2010;陈立春等，2010;中国地震局地质研究所，2010;孙鑫喆等，2010;张永双等，2010;张建毅等，2010;张军龙等，2010;石峰等，2010;马寅生等，2010;殷跃平等，2010a;刘云华等，2010;熊探宇等，2010;沈强等，2010;钱晓东等，2010;任俊杰等，2010;李金臣等，2010;Tobita et al.，2011;Lin et al.，2011)。

作者以灾后航空相片与高分辨率遥感影像目视解译为主，辅以野外调查的方法，建立了玉树地震滑坡数据库，库中包含了2 036条用面要素表示的玉树地震滑坡记录(许冲等，2011a，2012;Xu et al.，2012e)，这些滑坡主要分布在一个沿地表断裂带的面积约1 455.3km2的矩形区域内，本文就以这个矩形区域作为研究区。这些滑坡总面积为1.194km2，单体滑坡平均面积为586.4m2。研究区内滑坡面积百分比为0.082%;滑坡点密度为1.4个/km2。图1为玉树地震滑坡与同震地表破裂空间分布图。

2 滑坡影响因子专题图层的制作

结合玉树地震区的地震、地形、地质与其他资料，综合考虑对玉树地震滑坡的多种影响因素，进行滑坡影响因子的选择。综合考虑地震滑坡的外动力因子与内动力因子综合作用，选择地震滑坡的影响因子，地震滑坡的外动力因子为地震因子，选择同震地表破裂与PGA。内动力因子综合作用包括地形因子、地质因子与其他因子，地形影响因子选择高程、坡度、坡向、曲率、坡位、水系;地质影响因子选择地层岩性、与断裂的距离;此外，还考虑了公路与植被对地震滑坡的影响作用，公路影响采用与公路的距离来衡量，植被的影响选择NDVI来衡量。这样，选择了12个对地震滑坡有影响的因子，进行地震滑坡危险性评价。

图1 玉树地震滑坡与同震地表破裂空间分布图Fig.1 Spatial distribution of the 2010 Yushu earthquake triggered landslides and co-seismic surface ruptures.

2.1 地震因子

地震是地震滑坡的触发因素，本文选择地表破裂与PGA作为滑坡发生的地震影响因子。同震地表破裂对地震滑坡控制作用强烈，地震滑坡多沿同震地表破裂密集分布(许冲，2012;Xu et al.，2012)，因此，选择地表破裂作为地震滑坡评价的一个因子，以与地表破裂的距离作为同震地表破裂对滑坡的影响衡量指标。按500m的分级标准，将与地表破裂距离分为21类，分类标准依次为:1)0～500m;2)500～1 000m;3)1 000～1 500m;4)1 500～2 000m;5)2 000～2 500m;6)2 500～3 000m;7)3 000～3 500m;8)3 500～4 000m;9)4 000～4 500m;10)4 500～5 000m;11)5 000～5 500m;12)5 500～6 000m;13)6 000～6 500m;14)6 500～7 000m;15)7 000～7 500m;16)7 500～8 000m;17)8 000～8 500m;18)8 500～9 000m;19)9 000～9 500m;20)9 500～10 000m;21)＞10 000m。地震动峰值加速度(PGA)反映了地震中地面震动的强度，PGA越强，地震滑坡越容易发生，因此，选择PGA作为玉树地震滑坡的另一个外动力影响因子。PGA为来自于U.S.Geological Survey(2010)公布的PGA等值线图，研究区内PGA值的范围是0.12～0.40g。设定PGA分级标准为1)＞0.38g;2)0.34～0.38g;3)0.30～0.34g;4)0.26～0.30g;5)0.22～0.26g;6)0.18～0.22g;7)0.14～0.18g;8)＜0.14g。在ArcGIS软件中，将地表破裂缓冲区文件与PGA文件分别栅格化，空间分辨率均为5m×5m，这样制作出地表破裂距离专题图与PGA专题图。

2.2 地形因子

地形影响因子包括高程、坡度、坡向、曲率、坡位、与水系距离这几个因子。地震滑坡往往在一定的高程上较容易发生(许冲等，2009，2010c)，因此选择高程作为地震滑坡的地形影响因子之一;坡度是影响地震滑坡发生的重要因素，坡度大的区域是地震滑坡的易发区域;坡向、曲率、坡位、与水系距离这些地形因子均对地震滑坡具有一定的影响作用(许冲等，2009)。根据1/5万地形图上的等高线、高程点、水系等信息，制作数字高程模型(DEM)，其空间分辨率为5m×5m。一般情况下，考虑到大部分等高线的水平间距在25m左右，因此1/5万地形图对应着25m×25m分辨率的DEM，以保证不会丢失掉太多的地形信息，然而对于一些地形陡峭的地方，如等高线的水平间距＜25m的区域，会损失掉一些地形信息。而且，玉树地震滑坡的规模一般较小，25m×25m分辨率的DEM难以进行较客观的统计。因此，我们将1/5万比例尺的地形图插值为5m×5m分辨率的DEM，使用ArcGIS中的“topo to raster”功能。这种增加DEM分辨率的方法对于相邻等高线水平间距＞25m的地方，不会改变任何地形信息，但是对于相邻等高线水平间距＜25m的陡峭的地方会增加实际的地形信息。因此，制作了5m×5m分辨率的DEM来进行后续的地震滑坡危险性评价。研究区内高程范围是3 589.7～5 181.4m。高程分类依据200m一个间隔进行分类，分类标准为:1)＜3 800m;2)3800～4 000m;3)4000～4 200m;4)4200～4 400m;5)4400～4 600m;6)4600～4 800m;7)4800～5 000m;8)＞5 000m。

坡度、坡向、曲率、坡位这4个因子图层均是根据DEM在ArcGIS软件中导出的，它们的分辨率与DEM一样，都是5m×5m。坡度基于5°间隔进行分类，分类标准依次为:1)＜5°;2)5°～10°;3)10°～15°;4)15°～20°;5)20°～25°;6)25°～30°;7)30°～35°;8)35°～40°;9)＞40°。坡向分为9类，分类标准如下:1)Flat;2)N;3)NE;4)E;5)SE;6)S;7)SW;8)W;9)NW。将曲率分为12类，分别标准为:1)＜－1;2)－1～－0.1;3)－0.1～－0.05;4)－0.05～－0.01;5)－0.01～－0.005;6)－0.005～－0;7)0～0.005;8)0.005～0.01;9)0.01～0.05;10)0.05～0.1;11)0.1～1;12)＞1。坡位依据的分类标准为山脊、上坡、中坡、平坡、下坡、谷底共6类(Weiss，2006)。水系提取自地形图，将水系按照100m缓冲区进行分类，分类级别为:1)0～100m;2)100～200m;3)200～300m;4)300～400m;5)400～500m;6)500～600m;7)600～700m;8)700～800m;9)800～900m;10)900～1 000m;11)＞1 000m。在ArcGIS软件中将水系缓冲区图层栅格化为5m×5m分辨率的距离水系分级图，这样得到了水系影响因子专题图层。

2.3 地质因子

岩性与断裂为地震滑坡的2个重要的地质影响因子，岩性是地震滑坡发生的物质基础，断裂也与地震滑坡关系密切，断裂穿越区的岩性一般较破碎，进而影响着地震滑坡的发生，因此选择岩性与断裂作为地震滑坡评价的2个地质因子。本文使用的岩性与断裂因子均来自于1/20万地质图。图2为研究区地层岩性图，研究区内包括11类地层，依次为:1)Qh;2)Qhal-pl;3)N;4)T3bt;5)T3kn3;6)T3kn2;7)T3kn1;8)T2jl2;9)T2jl1;10)C—P;11)花岗岩。将岩性区文件栅格化为5m×5m的岩性专题栅格图。将与断裂距离按照如下分级标准进行分类:1)0～50m;2)50～100m;3)100～150m;4)150～200m;5)200～250m;6)250～300m;7)300～350m;8)350～400m;9)400～450m;10)450～500m;11)＞500m。将断裂缓冲区图层栅格化为5m×5m的图层，这样得到了断裂影响因子专题图层。

2.4 其他因子

图2 研究区地层岩性图Fig.2 Geological map of the study area.

除去地震、地质、地形等影响因子外，公路与归一化植被指数(NDVI)也是影响地震滑坡发生的2个因素。因为在公路建设过程中，人工路基的开挖会导致公路上坡的临空面较陡，进而容易在地震条件下发生滑坡，NDVI是归一化植被指数，反映植被覆盖情况，而植被覆盖对滑坡发生有一定的影响。将NDVI分为14类:1)＜0;2)0～0.05;3)0.05～0.1;4)0.1～0.15;5)0.15～0.2;6)0.2～0.25;7)0.25～0.3;8)0.3～0.35;9)0.35～0.4;10)0.4～0.45;11)0.45～0.5;12)0.5～0.55;13)0.55～0.6;14)＞0.6。公路图层来自1/5万地形图与基于SPOT 5数据的解译，将公路线状图按照100m缓冲区进行分类，分类级别为:1)0～100m;2)100～200m;3)200～300m;4)300～400m;5)400～500m;6)500～600m;7)600～700m;8)700～800m;9)80～900m;10)900～1 000m;11)＞1 000m。将公路缓冲区图层栅格化为5m×5m的图层，得到了公路影响因子专题图层。这样，制作了这12个影响因子的栅格格式的专题图，它们的分辨率均是5m×5m。

3 证据权方法及证据权值的计算

证据权方法是一种较成熟的地震滑坡评价方法，该方法以贝叶斯概率统计为基础。该方法的数学表达式为

根据玉树地震滑坡数据，分别计算每个影响因子的滑坡危险性权重。根据式(1)与式(2)，可得到

式中:Npix1表示该影响因子级别内发生滑坡的栅格数;Npix2表示该影响因子级别外发生滑坡的栅格数;Npix3表示该影响因子级别内未发生滑坡的栅格数;Npix4表示该影响因子级别外未发生滑坡的栅格数。

如公式(3)所示，在求取最终权重的步骤中，当前的研究有2种情况。一是将正权重与负权重相加，得到最终权重为

二是将两者相减，得到的最终权重为

本文分别采用这2种证据权方法进行滑坡评价，并对结果进行定量对比，研究哪种最终权重计算方法更合适。图3为12个因子的所有分级的权重计算结果图。可以看出在所有的影响因子分级内，两者表现出了几乎相同的趋势，只存在一些微小的差别。总权重(Wadd与Wsubtract)表现了影响因子对滑坡的重要性，如果总权重为正，表明该级别利于滑坡的发生，如果其为负，表明其不利于滑坡的发生。例如，坡度＞30°的总权重值为正，表明该分级利于滑坡的发生，而0°～30°的坡度范围不利于滑坡的发生。一些因子级别对滑坡的影响作用小，其取值在0附近。在证据权方法中，假设影响因子间是互相独立，不受彼此干扰的。证据权方法与其他二元统计方法(如信息量法、确定性系数法)相比，其优点是不但考虑了由滑坡发生的充分率得到的正权重，而且还考虑了由滑坡发生的必要率得到的负权重，最终权重的计算综合考虑了二者。根据图3的证据权值可以看出因子内部哪些级别的滑坡危险性高，距离地表破裂、水系、公路越近，滑坡危险性越高;坡度越大，滑坡危险性越高;4 000m以下是地震滑坡发生的危险高程;凸凹坡比光滑的斜坡发生滑坡的危险性更高;植被覆盖越低，地震滑坡的危险度越高。

4 2类滑坡索引图的制作与结果评价

得到最终权重值之后，给每个影响因子的分级赋予相对应的权重值，然后通过对各个影响因子层的叠加操作，来获得研究区的地震滑坡危险性索引(Landslide Hazard Index，LHI)图。叠加公式如式(8):

图3 权重计算结果图Fig.3 Weight of evidence values of the twelve factors.

这样，基于2种不同最终权重计算方法(Wadd与Wsubtract)，得到2个地震滑坡危险性评价索引图。运用Wadd方法得到的滑坡危险指数(LSI)范围为－15.805～7.614。运用Wsubtract方法得到的滑坡危险指数(LSI)范围为－17.454～9.312。

运用滑坡面积百分比累积－危险性面积百分比累积曲线下面积方法(Conoscenti et al.，2008;Kamp et al.，2008;Zezere et al.，2008;许冲等，2009，2010b，2011b)对2种结果的正确率进行定量检验。将危险性评价结果按照1%的面积间隔从高到低分为100等分，分别求取这100个级别内的滑坡发生的百分比，构建危险性结果百分比累加与滑坡百分比累积曲线，即结果评价曲线(图4)。基于2种证据权方法得到的结果评价曲线几乎重合。Wadd方法对应的曲线下面积(Area Under Curve，AUC)为80.32%，Wsubtract方法对应的AUC值为80.19%。两者相差很小，Wadd方法的正确率略大于Wsubtract方法，两者得到的结果均用已有滑坡数据得到了很好的验证，表明2种方法得到的结果都是理想的。

分别取“－3，－1，1，3”为阈值，按照证据权值将评价结果分为5类，分别代表滑坡危险性极低、低、中等、高、极高共5个等级，最终评价结果见图5。2种证据权方法得到的结果表现出类似的特点，地震滑坡高危险区分布在地表破裂的东端，这里为断裂拐弯终止的位置，地震在此处释放了大量的能量，造成了大量的坡体裂缝(许冲等，2012)，应加强这些地震滑坡极高危险区内的震后雨季滑坡、泥石流的排查、预防工作。