芦山地震震后次生滑坡灾害风险评价研究1

刘 睿 苏 伟 张晓东 田丽燕 张群燕

（中国农业大学信息与电气工程学院，北京 100083）

引言

北京时间2013年4月20日8时02分四川省雅安市芦山县（北纬30.3°，东经103.0°）发生了7.0级地震，震源深度13km。据中国地震局网站统计资料，震中芦山县龙门乡99%以上房屋垮塌，截至2013年4月24日上午10时，共发生余震4045次，其中3级以上余震103次，最大余震5.7级。此次地震受灾人口152万，受灾面积12500km2。四川省芦山“4·20”7.0级强烈地震触发了多处山体崩塌、滑坡，造成惨重的人员伤亡和巨大的财产损失，地震灾区山地灾害活动异常活跃，发生次生灾害的危险性非常大。地震次生灾害一般是指强烈地震后，以震动的破坏后果为导因而引起的一系列其它类型的灾害（谭德宝等，2008），包括崩塌、滑坡、塌陷、地裂缝和砂土液化等，是地震灾害的重要组成部分（杜军等，2009）。中国是一个地震多发国家，每年由于地震灾害而造成的损失巨大，其中地震滑坡是一种常见的地震次生灾害形式，它会对人类的生命财产和社会经济活动造成巨大损失，因此人们非常重视对地震次生灾害的研究。目前，对于地震次生滑坡灾害的研究，许多学者采用光学遥感数据、激光雷达数据以及SAR数据等多源遥感数据，针对不同地区，开展了包括地震次生灾害发生范围、强度等影响因子的分析，并提出了地震次生灾害敏感性分析的有效方法（唐川等，2002；李忠生，2003；李为乐等，2008；赵祥等，2009；李秀珍等，2009；祁生文等，2009；张建强等，2009；Daisaku等，2009）。

对于地质灾害的危险性分析，目前广泛应用的模型有：信息量模型、多元统计分析模型、模糊综合评判模型、人工神经网络遗传算法模型等，实现方法是将地理信息系统（Geographic Information System，GIS）技术与相关模型结合起来进行区域地质灾害危险性分析（夏毓超，2010）。对于地质灾害的危险性区划评价的研究（谢全敏等，2007；王萌等，2009；Mandy等，2001），由于各种地质因素本身的不确定性，以及地质因素之间相互作用的复杂性，往往需要以大量地质环境资料为基础，充分利用GIS强大的空间化功能，选取合适的评价预测指标，并运用数学分析模型划分地质灾害危险性等级，从而为地质灾害的预警决策和管理防治提供依据。本文在基于芦山地震受灾区航空影像解译滑坡灾害隐患点的基础上，综合运用GIS技术和层次分析法预测地震引起的次生滑坡灾害的危险性，探讨芦山地震的次生滑坡灾害风险性评价方法。

1 研究区及数据源

1.1 研究区概况（图1）

芦山县位于四川盆地西缘，雅安地区东北部，青衣江上游。北与汶川县连界，东北与崇州市、大邑县、邛崃市毗邻。地跨东经 102°52'至 103°11'，北纬 30°01'至 30°49'，县境南北长86.6km（飞仙关至断头岩）；东西宽：北部为24.4km（芦、崇、大三县（市）交界点至二十四凼），中部为19.42km（芦、大、邛三县（市）交界点至大瓮顶），南部为17.2km（芦、邛、雅三县（市）交界点至六台山）；幅员面积1364.42km2；县城距雅安31km，距成都 156km。芦山县位于龙门山前缘构造带南段，龙门山断裂带位于青藏高原东缘。龙门山断裂带上地震多发，2008年的汶川特大地震就在距离芦山县北部大川镇仅69km处发生。2013年04月20日08时02分在四川省雅安市芦山县（北纬30.3°，东经103.0°）发生7.0级地震，震源深度13km。

图1 研究区（芦山县龙门乡、宝盛乡、太平镇部分地区）（R（Band 3）G（Band 2）B（Band 1））Fig. 1 Location of the study area （R（Band 3）G（Band 2）B（Band 1））

1.2 数据源及数据预处理

本研究使用的数据源有无人机航拍影像、空间分辨率为30m的芦山县DEM、降雨量数据、芦山县行政区划图等。航空影像拍摄于2013年4月20日下午17时，范围包括芦山县核心灾区龙门乡、宝盛乡、太平镇部分地区，以具有地理参考的DEM为参考，对航空影像进行几何精校正。降雨量数据从天气历史气象网站上下载获得，下载地址为：

http://www.wunderground.com/history/station/56287/2013/4/23/DailyHistory.html?req_city=

NA&req_state=NA&req_statename=NA

下载后的气象数据为点位数据，利用Kriging插值方法进行空间插值处理扩展到面上。

2 滑坡遥感解译

滑坡是最为常见的一种地质灾害，由于主要受重力影响，地貌特征明显，可以通过遥感影像解译。滑坡的遥感解译主要是依据形态、色调、纹理特征来判定。在形态上，发生滑坡的区域在遥感影像上主要呈簸箕形、马蹄形、弧形或不规则形状；震后，在地震力的作用下，滑坡以长条形或不规则形状为主。在色调上，震后新形成的滑坡大多以松散堆积物质组成，表面波谱反射能力较强，在影像上以浅色调为主，与周围稳定地形有明显的区别，加之坡体破碎，纹理粗糙。根据这些解译特征，参考地震前视宝卫星影像对滑坡灾害点进行了详细的解译，研究区内共解译滑坡97处，图2分别是解译得到的庐山地震次生滑坡灾害面状解译结果和点状解译结果图。通过遥感解译，芦山地震次生滑坡灾害主要发生在震中的北部和西北部山区等坡位较高的地方，这与汶川地震次生滑坡灾害发生的地质条件相似，在今后开展的地震次生滑坡灾害实地调查中将会得到验证。

图2 芦山地震次生滑坡灾害目视解译结果Fig. 2 Results of the earthquake landslides disaster by visual interpretation

3 评价因子指标选取与敏感性分析

3.1 评价因子指标选取

影响地震次生灾害发生的因素（即灾害致灾因子）众多（夏毓超，2010），其中包括地形地貌、地层岩性、地质构造等影响地质灾害形成的基本地质环境因素，还包括降雨、地震、河流冲刷作用和人类工程活动等的影响，以及诱发地质灾害演化和发生的外在因素。综合国内外研究者所采用的地震次生灾害评价因子，本文确定的地震次生滑坡灾害风险的基本因素评价因子包括：地形地貌因子（坡度、坡向）、环境因子（危险植被指数）和气象因子（震后累积降雨量）。

3.1.1 地形地貌因子

地面坡度对滑坡发育有着重要的影响，总的来说，坡度越大越容易引发滑坡、崩塌等地质灾害。坡向主要反映光照情况，是山坡平面法线在水平面上投影的方向，不同的光照条件可能导致不同的山坡风化程度。一般而言，阳坡植被茂密、冰雪易于融化，岩石的风化速度也快于阴坡，且降水量较多，阳坡上的地质灾害发生的可能性大于阴坡。另外，气流在迎风坡上升的过程中降温降水，而背风坡在下降过程中增温不易于水汽凝结，也导致迎风坡地质灾害发生的可能性大于背风坡。本文采用空间分辨率为30m的DEM数字高程模型数据，并用ARCGIS软件对DEM进行空间运算，计算出研究区每个像素的坡度、坡向，如图3所示。

图3 芦山地震次生滑坡灾害地形地貌因子分析结果Fig. 3 Results of the earthquake landslides disaster by using slope factor

3.1.2 气象因子

降雨量是地震次生滑坡灾害的气象激发因子，当足够的雨水渗入岩土体空隙中，岩土体自重及孔隙水压力均增加，随着雨水的持续渗入，土体达到饱和状态，土石体原本稳定的结构遭受破坏，阻抗力变小，使得松散固体物质变得越来越不稳定，就容易发生滑坡等灾害。在特定的地区，孕育滑坡的客观环境条件一定的情况下，在一定时间内降雨量超过某一临界阈值，就会发生滑坡。所以本文选取2013年4月20日至4月26日的雅安及周边9个气象站（九龙、汉源、峨眉山、康定、乾宁、成都、都江堰、绵阳、宜宾）的气象数据，根据各站点的降雨量，利用克里金（Kringing）插值方法，得到芦山县太平镇、宝盛乡和龙门乡2013年4月20日至4月26日累积降雨量数据。

3.1.3 环境因子

当一个位置发生滑坡现象之后，植被覆盖度降低，植被指数降低，根据滑坡前后植被指数的变化可以分析判断一个位置是否发生滑坡现象。因此，本文利用芦山县宝盛乡、龙门乡、太平镇共3个乡（镇）地震前后的航空遥感影像，通过计算震前、震后的植被指数，进而构建危险植被指数，作为地震次生滑坡灾害的环境因子，结果如图4所示。

收集到的地震前后的遥感影像有震前视宝卫星影像和震后航空影像，但由于航空影像只有红光、绿光、蓝光三个波段，所以我们构建了可见光抗大气指数（Visible Atmospherically Resistant Index，VARI）来计算植被指数。可见光抗大气指数不仅可以提高对土壤背景的鉴别能力，消除地形和群落结构的阴影影响，还极大地削弱了大气对于植被覆盖度检测的干扰，扩展了其灵敏度。可见光抗大气指数的计算公式为：

图4 芦山地震次生滑坡灾害环境因子分析结果Fig. 4 Results of the earthquake landslides disaster by using environmental factor

式中，G为绿波段；R为红波段；B为蓝波段。

根据地震前后植被指数变化计算危险植被指数（Vegetation Damage Index，VDI）的计算公式为：

式中，VARIpre指震前植被指数；VARIpost指震后植被指数。

VDI越高，说明此区域越危险，越容易发生滑坡现象。

3.2 敏感性分析

敏感性分析采用统计学方法，并充分利用 GIS的分析功能（高克昌等，2006），对解译出来的滑坡灾害点与坡度、坡向、震后累积降雨量和危险植被指数4个指标因子的属性进行定量计算。具体方法为：①数据预处理。利用克里金插值法对2013年4月20日至4月26日的累积降水量进行空间插值；②针对每个因子图层，进行属性划分，分级的方法包括分位数分级法、标准差分级法、自定义分级法等；③敏感性分析，利用GIS的空间分析功能统计出各等级滑坡隐患点的个数，具体计算公式为：

式中，Fn为某指标因子的属性n在每100km2内隐患点密度，即滑坡隐患点出现对属性n的敏感性；Nn为指标因子中属性n所包含的隐患点个数；Sn为属性n的面积。

在对各个指标进行敏感性计算之后，需要对各指标进行归一化处理。根据得到的敏感性的值将各指标因子划分为四级，比如，将敏感性最高的1至2个属性类赋值为4，敏感性较高的1至2个属性类赋值为3，敏感性较低的1至2个属性类赋值为2，敏感性最低的属性类赋值为1。对每个指标而言，指标中值越高的属性类，对地震次生滑坡灾害的影响程度越大，它所在区域发生地震次生滑坡灾害的危险性也就越高。

3.3 敏感性分析结果

坡度以1×1网格为分析单元，将坡度划分为<89.9899m、89.9899m—89.9950m、89.9950m—89.9980m、89.9980m—89.9996m共4个等级，敏感性从1到4逐渐增加，分析结果如表1所示；坡向属性以正北方向顺时针旋转，对8个方向进行敏感性分析，分析结果如表2所示，结果表明，北、东北和东方向地震次生滑坡敏感性最高；其次为南、西南和西方向；东南和西北方向敏感性较低。

表1 坡度敏感性分析结果Table 1 Sensitivity analysis of slope factor

表2 坡向敏感性分析结果Table 2 Sensitivity analysis of slope orientation factor

研究区VDI值分为九级，经分析后得到敏感性从4到1依次递减，如表3所示。危险植被指数越高，发生滑坡灾害的可能性越大，统计得到危险植被指数在-0.19到1时，地震次生滑坡敏感性较高，危险植被指数小于-0.41时，地震次生滑坡敏感性较低。

表3 危险植被指数敏感性分析结果Table 3 Sensitivity analysis of the VDI factor

续表

研究区内2013年4月20日至4月26日的累积降雨量也分为九级，其敏感性分析结果如表4所示。

表4 年均降雨量敏感性分析结果Table 4 Sensitivity analysis of rainfall factor

4 地震次生滑坡灾害危险性评价

AHP（Analytic Hierarchy Process）层次分析方法，是20世纪70年代由美国运筹学家Saaty提出的，主要是对一些较为复杂、较为模糊的问题作出决策的多准则决策方法，可用于存在不确定情况及多种评价标准的决策问题，它基于对问题的全面考虑，将定性与定量分析相结合、将决策者的经验予以量化，是一种层次权重决策分析方法（樊晓－，2004；Liang，2009）。本文采用AHP方法对诸多影响因素进行综合分析，根据对影响次生滑坡灾害评价因子的分析，对各因素的重要性进行两两比较，采用 1—9标度法使各因子相对重要性定量化，得出权重值。

4.1 构造判断矩阵

根据层次分析法的基本原理，建立判断矩阵，并将各个指标之间的相对重要程度表示出来，判断矩阵标度（重要性指标）及其含义如表5所示。在研究中主要通过查阅文献和向专家咨询的方式，进行指标因子之间的相互重要性打分，并根据层次分析法原理构建判断矩阵，如表6所示。

表5 判断矩阵标度（重要性指标）及其含义Table 5 The judgment matrix and its concept

表6 指标判断矩阵和权重Table 6 The Factor judgment matrix and weights

计算特征向量并归一化处理，得到权重数W及最大特征根为：

依据以上计算得到的特征向量，对判断矩阵进行一致性检验，公式如下：

式中，CI为判断矩阵一致性指标；RI由大量的实验给出，为判断矩阵的平均随机一致性指标。

根据层次分析法的规定：若 CR<0.1，则认为判断矩阵具有令人满意的一致性，否则就需要调整判断矩阵，直到满意为止。由上述计算结果可知，本次构造的判断矩阵一致性较好。

4.2 地震次生滑坡灾害危险性评价

利用ARCGIS的空间分析功能，采用栅格加权叠加，对4个影响次生滑坡灾害发生的因子进行统计分析，得出各栅格单位的地震次生滑坡灾害危险度值，并将所有栅格单元分为高危险区、中危险区和低危险区三类，如图3所示。

根据芦山地震次生滑坡灾害危险性评价图，高危险区主要分布在坡度较高的山区，该区域断裂构造发达，岩石破碎，且近年来水土破坏严重，植被覆盖低，危险植被指数高，震后累计降雨量也相对较高，各项评价指标危险性都高，因此综合的发生地震次生滑坡灾害危险性也很高；危险程度较低的地区区域地质构造简单，断裂带密度低，且危险植被指数低。

图5 芦山地震次生滑坡灾害风险评价结果Fig. 5 The estimation result of the earthquake secondary landslides disaster

5 结论和分析

（1）本研究在对芦山地震次生滑坡灾害危险性评价的因子选取中引入了危险植被指数、震后累积降雨量、坡度、坡向4个因子，针对不同因子在敏感性计算结果的取值范围上的差异，采取了敏感性归一化的方法进行量化，得到各指标的敏感区间。选取危险植被指数、震后累积降雨量2个因子作为评价因子是本研究的创新之一。由于震后数据是航空影像，只有红、绿、蓝三个波段，所以本研究基于可见光抗大气指数VARI植被指数计算危险植被指数，对于研究区的植被覆盖情况有效地进行了表达。

（2）应用AHP层次分析法对地震滑坡危险等级进行划分是地震滑坡问题研究中的尝试，层次分析法从理论基础角度来讲，是一种较好的确定影响滑坡指标因子权重的方法，并且在一定程度上克服了模糊综合判断多因子和多指标问题，减小了主观因素的影响，也能解决概念外延不明确的问题。

（3）充分利用ARCGIS的空间分析功能，采用栅格加权叠加，通过统计分析，对于芦山地震次生滑坡灾害进行了风险性评价，得到了高危险区、中危险区、低危险区三级危险区，这为灾后重建方案的制定提供了理论基础。

杜军，杨青华，2009. 基于GID和AHP耦合的汶川震后次生地质灾害风险评估. 中国水土保持，（11）：14—16.

樊晓一，乔建平，陈永波等，2004. 层次分析法在典型滑坡危险度评价中的应用. 自然灾害学报，13（1）：72—76.

高克昌，崔鹏，赵纯勇等，2006. 基于地理信息系统和信息量模型的滑坡危险性评价——以重庆万州为例. 岩石力学与工程学报，25（5）：991—993.

李为乐，唐川，杨武年等，2008. RS和GIS技术在县级区域泥石流危险区划中的应用研究——以四川省泸定县为例. 灾害学，23（2）：71—76.

李秀珍，孔纪名，邓红艳等，2009. “5·12”汶川地震滑坡特征及失稳破坏模式分析. 四川大学学报（工程科学版），41（3）：72—77.

李忠生，2003. 国内外地震滑坡灾害研究综述. 灾害学，18（4）：64—70.

祈生文，许强，刘春玲等，2009. 汶川地震极其重灾区地质背景及次生斜坡灾害空间发育. 工程地质学报，17（1）：39—49.

唐川，朱大奎，2002. 基于GIS 技术的泥石流风险评价研究. 地理科学，22（3）：300—304.

谭德宝，申邵洪，陈蓓青，2008. 空间信息技术在地震次生灾害防治中的应用. 人民长江，39（22）：99—102.

王萌，乔建平，2009. 基于GIS的沐川县滑坡危险度区划. 自然灾害学报，18（1）：174—179.

夏毓超，2010. 地震次生滑坡危险性分析. 成都：西南交通大学.

谢全敏，李道明，翟鹏程，2007. 滑坡次生灾害损失评价方法研究. 岩土力学，28（5）：961—970.

赵祥，李长春，苏娜等，2009. 滑坡泥石流的多源遥感提取方法. 自然灾害学报，18（6）：29—32.

张建强，范建容，严冬等，2009. 地震诱发崩塌滑坡敏感性评价——以北川县为例. 四川大学学报（工程科学版），141（32）：140—145.

张友谊，胡卸文，朱海勇等，2007. 滑坡与降雨关系研究展望. 自然灾害学报，16（1）：104—108.

Daisaku Kawabata, Joel Bandibas, 2009. Landslide susceptibility mapping using geological data, a DEM from ASTER image and an artificial Neural Network (ANN). Geomorphology, 113: 97—109.

Liang S.Y., Wang Y., Zhou Z.Q., Yang Y.T., 2009. Geological hazard assessment of small towns based on EAHP.See: Internatinal Conference on Management and Service Science, MASS 2009. IEEE C.S., N.J., United States,ISBN 978-1-4244-4639-1.

Mandy Lineback Gritzner, Andrew Marcus W., Richard Aspinall et al., 2001. Assessing landslide potential using GIS, soil wetness modeling and topographic attributes, Payette River, Idaho. Geomorphology, 37 (1-2): 149－165.