基于时序遥感数据的九寨沟核心景区崩塌发育分布特征

梁靖， 王彦东*， 裴向军

(1.地质灾害防治与地质环境保护国家重点实验室， 成都 610059； 2.中铁二院工程集团有限责任公司， 成都 610031)

2017年8月8日21时19分，世界自然遗产九寨沟风景名胜区(33.20°N，103.82°E)爆发7.0级地震(简称“8·8地震”)。“8·8”地震触发了数以千计的高位或极高位崩塌灾害，是震区发育多、分布广的灾害类型，并造成九核心景区交通旅游设施及景观严重破坏。同时，震后三年的持续关注表明现阶段崩塌活动性仍较高，如景区芦苇海危岩崩塌在汛期降雨时危岩带变形与堆积体失稳频发，此类现象严重影响了景区开放运营、村寨与游客安全以及区域生态环境安全。因而，开展地震前后崩塌灾害的时序解译与分析，对景区地质灾害风险管控与防灾减灾具有重要现实意义。

一直以来，地震地质灾害调查以野外调查为主。随着遥感与地理信息技术的快速发展，基于遥感影像解译的地质灾害调查与分析成为一种更为高效、快捷的手段，并具有成本低、速度快、大范围覆盖等优势[1-3]。得益于此，基于遥感影像的解译识别技术在九寨沟震后抢险救灾与景区灾后恢复重建工作中发挥了重要作用。目前，九寨沟震后灾害遥感解译主要是基于光学遥感检测地震前后地表的各项特征信息变化，再利用目视解译或计算机自动解译等方式实现灾害识别[4]。例如，王绚等[5]便基于九寨沟景区雨季前后3期高分辨率遥感影像解译震后地质灾害，研究了震后地质灾害活动性与地质灾害易发性。董秀军等[6]结合高分辨率光学影像数据和机载激光雷达开展了景区地质灾害解译，并利用光学影像和数字高程模型综合建立了灾害三维解译标志，为景区恢复重建提供了技术支撑。同时，还有学者基于高精度遥感影像进行地震滑坡的分布及范围等信息提取，建立自动提取方法[7-8]。除此之外，国外学者还进行了更早的探索。Mahdavifar等[9]在Iran地震的现场调查的基础上，利用遥感技术分析了地震崩塌分布与形成特征。Dunning等[10]结合遥感影像与实地调查开展了2005年Kashmir地震地质灾害的调查解译工作，分析了地震触发崩塌滑坡的发育形成和失稳特征。

因此，现以“8·8”九寨沟地震前后的崩塌灾害为研究对象，结合多年实地调查、时序遥感影像解译与核密度分析等手段，通过地震前后共计5期高精度遥感影像建立崩塌灾害数据库，进一步分析崩塌空间聚集趋势及发育分布特征，以期为研究区震后重建及生态修复等治理工程提基础资料。

1 研究区概况

九寨沟位于青藏高原与四川盆地大地貌单元过渡的深切高山峡谷地带，地势上南高北低，海拔高度普遍高于2 000 m，平均地形坡约50°，出露地层为寒武系-三叠系地层灰岩、白云质灰岩等。同时，研究区属于高原寒温带-亚寒带季风气候带，垂直分带特征明显，可见年均降雨量为704.3 mm，其中4—10月降雨量占全年的89.61%。此外，由于研究区地处青藏高原东北缘、龙门山断裂西北部的构造活跃区，活跃的构造活动使得地震频发。在20世纪以来共爆发过8次Ms≥6.0的历史地震(图1)，长期强烈的构造活动导致其地质环境脆弱，多期次历史地震也致使山体具有不同程度“内伤”。

“8·8”地震震中位于核心景区以西5 km的比芒村，震源深度达20 km，最大烈度达Ⅸ度，且景区均处于Ⅶ到Ⅸ度的影响范围内。结合现场及同震解译分析表明，九寨沟地震具有如下特点：①地震震级大、震源深、影响范围广；②地震触发的次生地质灾害严重，高位崩塌普遍发育(图2)；③地震对自然景观和生态环境影响严重。

PGA为峰值地面加速度图1 研究区典型历史地震数据[3]Fig.1 Typical historical seismic data in the study area[3]

图2 九寨沟地震区崩塌发育情况Fig.2 Development of collapse in Jiuzhaigou earthquake area

2 数据与方法

收集了九寨沟地震前后共计5期高精度遥感影像和无人机航拍影像(表1、图3)，基于此进行地质灾害解译。为了提高解译的准确性，首先对遥感影像进行几何校正、色彩调整、精确配准、坐标系转换等系列处理。依据《滑坡崩塌泥石流灾害调查规范》及《环境地质遥感监测技术要求》等技术规范，采取人机交互方法详细解译研究区内的地质灾害，结合实地核查及相关报道成果[6]来确定崩塌灾害的解译标志，并利用ARCGIS软件获得区内崩塌灾害的解译成果。

同时，基于L函数与R语言，结合核密度分析来获得来研究崩塌的空间聚集特征。L函数是一种对点数据集进行不同距离的聚类分析方法，可将坡度、岩性等因素引入对空间分布分析中[11]。其是在Ripley’s K函数的改进基础上提出，通过计算随机分布点在半径距离为r下额外随机点的期望数目，从而判断点数据在空间上的分布聚集性[12]。需要

图3 研究区地震前后高清遥感影像Fig.3 High-definition remote sensing images before and after the earthquake

表1 遥感数据来源

指出的是，该函数可以克服自然非恒定强度和事件聚集之间的误差[13]。L函数具体定义为

(1)

式(1)中：K(r)是Ripley’sK函数；r是通过L函数计算得到的理论半径。

再借助R语言平台的Spatstat软件编程计算，可得到L函数绘制的曲线，其中正值代表事件聚集，负值代表事件离散。最后基于ARCGIS中的核密度函数工具绘制不同地质灾害核密度图[14]，该方法广泛运用于滑坡分布研究[15]，从而评价滑坡的空间分布密度。

3 地震前后崩塌编录及空间聚集性

3.1 地震前后崩塌编录

借助ARCGIS软件，对震前(2017-08-05)、同震(2017-09-06)及震后(2018-08-22、2019-08-16、2020-08-28)的遥感影像进行详细解译，每期解译范围达655 km2，得到崩塌的时序解译结果(图4、表2)，可为崩塌灾害发育分布研究提供基础数据。

图4 地震前后崩塌解译结果Fig.4 Interpretation results of collapse before and after the earthquake

对于震前影像(2017-08-05)：研究区在地震前的地质灾害发育程度较弱，主要是由降雨、人类活动和河流下切所引发。并解译出震前崩塌灾害135处，面积约1.72 km2，面积占比0.26%。对于同震影像(2017-09-06)：研究区受强烈的地震力而导致地质灾害激增，总体呈现五花海-箭竹海与丹祖沟强发育、树正沟中等发育、则查洼沟发育程度弱的趋势。共计解译同震崩塌灾害1 219处，面积约6.91 km2，面积占比1.05%。对于震后影像(2018-08-22、2019-08-16、2020-08-28)：可以发现，震后三年内崩塌灾害数量及面积仍呈现缓慢增加趋势，崩塌从2018年的1 430处增加到2020年的1 528处，灾害面积也由7.85 km2扩大至8.43 km2。

综合时序遥感影像的地质灾害编录结果，不难看出研究区易发灾害类型主要为崩塌，并且经过“8·8”地震影响，震后三年内其数量、面积等持续增加，仍具有较高活动性。

3.2 崩塌空间聚集特征

基于R语言得到了地震前后崩塌L函数的计算结果(图5)。由图可知崩塌分布在0～5.6 km范围值大于0，表明其在此范围内具有聚集性。而在分布频率方面，崩塌计算曲线上在1.8 km处表现出最大。需要指出的是，本文中核密度分析对应的搜索半径为L函数分析中获得的最大聚类距离值，因此以1.8 km作为搜索半径，同时对崩塌灾害点赋予面积权重，并得到面积赋权核密度图(图6)。

表2 崩塌解译结果统计表

图5 崩塌灾害的Ripley’s L函数分布Fig.5 Ripley’s L function distribution of collapse

图6 地震前后崩塌面积赋权核密度图Fig.6 The weighted core density map of the collapsed area before and after the earthquake

由图6可知，地震前后崩塌的空间聚集趋势显著不同，且面积赋权核密度的分布位置不尽相同。地震前，崩塌的面积赋权核密度主要集中树正沟及长海流域；而震后则是五花海-箭竹海区域、丹祖沟中部密度等级最高，且沟域内广泛发育崩塌灾害。单独分析震后三年的崩塌聚集趋势，可见最大聚集区域区别并不显著，但聚集范围存在局部扩大，这说明震后的崩塌聚集趋势总体不变，空间上聚集在五花海-箭竹海区域，且每年均会发生不同程度活动。

4 地震前后崩塌发育分布特征分析

基于前述对崩塌空间聚集趋势的分析，此处选取高程、坡度、坡向及岩性四类影响因素，进一步研究时序影像下的崩塌发育分布特征。

图7 不同高程下的崩塌分布Fig.7 Collapse distribution at different elevations

4.1 不同高程下的崩塌发育分布

研究区崩塌高程主要分布于2 200～4 200 m，为便于分析，将高程范围按间隔200 m划分为11个等级。图7(a)显示，地震前大多数崩塌集中在高程3 000～3 600 m内；“8·8”地震后不久，崩塌主要集中在高程2 800～3 400 m，但在0～2 200 m与4 000 m以上并未分布。而在震后三年内，崩塌主要集中2 800～3 800 m内，崩塌面积也逐年缓慢增加，可见震后崩塌仍持续活跃。由图7(b)可知，同震与震后的变化趋势大致相同，80%以上的崩塌集中在2 600～3 800 m内。此外，震后最高面积比在海拔3 200 m左右，震前则在3 600 m左右。由此可见震后崩塌在不同高程等级上仍然处于较高活动性，并无明显衰减趋势。

图8 不同坡度下的崩塌分布Fig.8 Collapse distribution under different slopes

4.2 不同坡度下的崩塌发育分布

为分析崩塌面积在不同坡度的分布，将研究区地形坡度从5°～75°分为15个等级，间隔为5°，绘制的崩塌面积分布见图8(a)。可以发现，地震前后崩塌的地形坡度分布趋势几乎相似，主要集中在30°～55°内；而在震后2018—2020年期间，崩塌面积在所有的坡度等级中都出现了小幅增加，但在大于65°和小于15°地形坡度中分布依旧较小。此外，统计不同坡度等级的崩塌面积比例，如图8(b)所示，崩塌的震后发展趋势大致相似，且85%的崩塌始终处于30°～55°的坡度等级。进一步分析发现，崩塌在震前集中在35°～45°，而震后主要集中在40°～50°。可以认为，崩塌面积及其占比在地形坡度上的分布趋势基本一致，并且大部分崩塌主要集中在陡坡段。

4.3 不同坡向下的崩塌发育分布

将研究区发育的崩塌灾害与地形坡向进行统计，结果如图9所示。可以发现，大多数崩塌均发育在北东-南方向，其中东方向尤为显著；其次为南-北西方向。与高程与坡度相比，地震前后崩塌的发育分布趋势仍大致相同，仅发育面积有所差别。具体而言，震前面积约0.73 km2的崩塌发育于东向，但面积占比却达到43%。震后崩塌面积逐年递增，至2020-08-28已发育2.55 km2，而面积占比相较于震前减少至30%。崩塌在坡向上的发育分布趋势表明，地震触发与震前重力触发的崩塌均发育在北东-南的优势方位。此外，震后崩塌面积及其占比在所有方位均有所增加，表明地震释放的强烈能量还严重影响着非优势方位，导致崩塌大范围发育分布。

图9 不同坡向下的崩塌分布Fig.9 Distribution of collapses on different slopes

4.4 不同岩性下的崩塌发育分布

崩塌编录结果与研究区地层岩性的统计分析表明(图10)，地震前后崩塌的发育及分布趋势基本一致。崩塌面积在地震前后均为石炭系(C)与石炭二叠系(Cp)发育最多[图10(a)]，震前分别为0.7、0.8 km2，在震后的2020年达到2.3、4.9 km2。同时，石炭系(C)与石炭二叠系(Cp)崩塌面积占比也较高[图10(b)]，但震后Cp发育的崩塌占比相较于震前降低了约13%，而石炭系则增加至58%。除此之外，第四系(Q)、三叠系(T)、二叠系(P)、泥盆系下吾拉组(Dx)及泥盆系益硅沟组(Dcy)发育崩塌面积小，占比也低。可以认为，地震前后的崩塌在石炭及二叠系石灰岩地层均表现出高倾向性，显示了厚层坚硬碳酸盐岩地层对崩塌灾害发育分布具有重要控制作用。

图10 不同岩性下的崩塌分布Fig.10 Collapse distribution under different lithologies

5 结论

利用多期次高精度光学遥感影像数据，采用实地调查、遥感解译与核密度分析等方法，分析了崩塌空间聚集趋势及发育分布特征，为九寨沟震区的地震地质灾害时空演化规律与长期效应研究提供了数据支撑。主要结论如下：

(1)对震前(2017-08-05)、同震(2017-09-06)及震后(2018-08-22、2019-08-16、2020-08-28)遥感影像进行时序解译发现，震前崩塌发育极少，同震期间受强烈地震力而导致崩塌激增，共解译1 219处，面积约6.91 km2；震后三年内崩塌数量面积一直缓慢增加，崩塌从2018年的1 430处增加到2020年的1 528处，面积也由7.85 km2扩大至8.43 km2。

(2)崩塌在0～5.6 km范围内聚集性显著，但地震前后聚集地域有所不同。其中，震前集中在树正沟及长海区域，震后则集中于五花海-箭竹海与丹祖沟中部，其面积赋权核密度等级最高。总体而言，2018—2020三年内的总体聚集趋势变化不大，但每年仍有一定程度改变。

(3)地震前后崩塌的发育分布规律分析表明，不同高程、坡度、坡向、岩性及降雨强度的控制性规律为：崩塌显著分布于高程2 800～3 400 m与30°～55°的陡坡段，地震触发与震前重力诱发崩塌多分布在北东-南的优势坡向，岩性上集中分布于石炭系与石炭二叠系的厚层坚硬碳酸盐岩地层。