2017-08-08，四川省阿坝藏族自治州九寨沟发生Ms 7.0级地震，触发了1 883处同震地质灾害^[16-17]，造成九寨沟景区景观损毁严重，景区公路及栈道均不同程度受损。熊猫海及周边区域是受灾最严重的区域之一。因此，本文选取熊猫海周边约52 km²作为研究区域（见图 1）。研究区位于九寨沟县漳扎镇九寨沟国家森林公园内（103.82°E，33.20°N），岩性以生物碎屑灰岩与灰质白云岩为主，具有深切峡谷地貌，因此在震前地质灾害就十分发育^[18-22]。然而，该区域植被覆盖率高达80.39%^[23]，普通光学影像无法辨识植被下的真实地形，造成对该区域震前地质灾害的位置、特征与发育分布规律的认识十分不足，导致无法深入研究九寨沟震区震前与同震及震后地质灾害的关系与时空演化规律。

图  1  研究区地理位置及地震概况

Figure 1.  Geographical Location and Earthquake Survey of the Study Area

本文采用机载LiDAR技术获取九寨沟景区五花海-日则沟区域的点云数据（四川测绘地理信息局提供），生成了高精度DEM（0.5 m），数据获取时间为2018年10月。机载LiDAR技术与传统航测相比，不仅速度与精度更高，而且自动化程度高，受天气影响小。相比传统卫星航测生成的DEM，机载LiDAR技术获取的点云数据处理后得到的DEM精度可达0.1~0.5 m，且去除植被影响后，能反映真实地形，更有利于辨识区域地质构造和地质灾害。研究区震前DEM精度为12.5 m，生成的山体阴影图如图 2所示，震前DEM精度较低，无法用其进行历史地质灾害解译。震后LiDAR DEM精度较高，清晰地反映了研究区的地貌特征。因此，本文主要基于震后LiDAR DEM的三维地形与地貌特征，进行震前滑坡的解译。值得指出的是，九寨沟地震诱发的同震滑坡绝大部分为浅表层滑坡且数量有限（如图 2（a）中的无植被区域），因此不会对震前地形产生明显影响。采用震后高精度LiDAR DEM仍然可以获取震前历史地质灾害，特别是大型地质灾害的地貌特征，帮助识别和解译震前地质灾害。

图  2  影像、12.5 m与0.5 m DEM生成的山体阴影图对比

Figure 2.  Comparison of Optical Image, Maps Generated by 12.5 m and 0.5 m DEM Hillshade

基于DEM的滑坡识别一般是通过生成山体阴影图、坡度图、坡向图以及地形粗糙度图等获取地貌特征，从而进行滑坡解译。使用ArcGIS栅格表面工具生成山体阴影图时，需要输入太阳光入射角与太阳高度角，不同的角度造成山体的阴影在不同方位，这会在一定程度上影响滑坡的识别（阴影区的滑坡常无法识别）。为了更好地识别滑坡，文献[24]提出了RRIM这种基于三维数据的二维可视化方法，其由地形正开度、地形负开度与坡度叠加形成。RRIM方法可补充现有地形可视化手法的不足，并将环境光的影响应用于地形显示，使目视解译能够更清晰地辨认地貌差别，更有利于山区地质灾害的识别与精准解译。

文献[25]提出通过计算一定步长内的最大天顶角与天底角来表现地面的凹凸程度，即地形开度。地形负开度表示曲面的凹度，取值范围一般在0~2之间，值越大说明该处地形越向内凹，在图像上负开度值最大的区域一般为河流、山谷及冲沟等。地形正开度表示曲面的凸度，取值范围一般在0~2之间，值越大说明该处地形越向外凸，在图像上正开度值最大的区域一般为山脊与坡顶。地形开度具有与光源无关的特性，因此可消除其他可视化方法中地形阴影的不利影响，而且计算结果受DEM噪点干扰的影响要比其他参数小。

地形开度以栅格DEM计算，基于一个栅格点与8个方向的目标栅格点以直线相连，以半径L为步长，得到8个方向的地形断面。分别计算8个方向目标栅格点的最大天顶角与天底角。地形正开度为8个方向上最大天顶角的平均值，地形负开度则为8个最大天底角的平均值。地形正开度值越高，说明地形越凸；地形负开度值越高，说明地形越凹。地形开度使地形的凹凸更加明显，突出了微地貌信息^[26]。

根据文献[24]提出的公式计算脊谷指数I：

式中，O_P为地形正开度；O_N为地形负开度。

脊谷指数的取值范围一般为－1~1。当脊谷指数为正时，值越大代表该处地形越凸，在图像上亮度越高，如山脊等区域；当脊谷指数为负时，值越小代表该处地形越凹，在图像上越暗，如道路、湖泊与侵蚀沟等区域。在脊谷指数图上，地质灾害源区与堆积区的明暗界线十分清晰，并且可以通过观察源区的光暗程度，细致地区分出滑坡后缘、崩塌后壁等灾害源区的界线。

将脊谷指数I与地形坡度叠加，且脊谷指数I图层以渐变灰色表示，地形坡度以渐变红色表示。生成RRIM图的流程如图 3所示，先根据DEM计算出地形正开度、地形负开度及坡度，再采用式（1）计算脊谷指数I，将坡度图层叠加在脊谷指数图层之上，调整两个图层的色带以及透明度，最后得到可以清晰辨认地质灾害边界的RRIM图。根据文献[24, 27]的经验，选用红色为底是因为红色对人眼来说较为显眼，能最大限度地展现细微的地貌信息。通过RRIM展现的详细地貌信息，更有利于辨识植被覆盖地区的地质灾害。

图  3  RRIM生成过程图及滑坡解译

Figure 3.  RRIM Process and Landslide Interpretation

为了说明RRIM方法较山体阴影图的优点，本文选取九寨沟漳扎镇某区域进行了对比分析。由图 4可见，由于植被覆盖，通过震前普通光学影像无法辨识地质灾害。通过震后高精度遥感影像可以准确获取同震滑坡范围（植被破坏区域），但依然无法解译震前地质灾害。通过LiDAR DEM可以确定地质灾害的范围，但RRIM方法能够更好地反映地质灾害的地貌特征，实现精细遥感解译。此外还可发现，DEM精度越高，RRIM效果越好。可见，高分辨率LiDAR数据结合RRIM的方法可以很好地帮助精细解译植被茂密地区的地质灾害。

图  4  影像、山体阴影图与RRIM图对比

Figure 4.  Comparison of Image, Hillshade and RRIM

山区地质灾害解译主要是通过遥感影像、山体阴影或其他可视化图像上地质灾害的形态、阴影、纹理以及色调等进行^[28]。历史地质灾害往往被植被覆盖，因此较难通过可见光影像解译。此外，图像的分辨率对于识别山区地质灾害也十分重要。根据山区常见地质灾害——崩塌、滑坡和泥石流的地貌特征，可以进行地质灾害的识别与解译。本文总结了各类地质灾害的地貌解译特征，并对比了RRIM图和山体阴影图，可见RRIM图能够更清晰地展现微地貌特征，有助于地质灾害的精细解译（见图 5）。

图  5  山区地质灾害地貌解译特征示意图

Figure 5.  Characteristics of Geohazards Interpretation in Mountainous Areas

崩塌灾害地貌特征明显，多发育在陡坡或岩体破碎并存在临空面的区域，常成群发育，具有明显的崩塌后壁与崩塌堆积体。崩塌后壁陡峭且可见破碎的岩体结构，色调偏暗，部分崩塌后缘可见参差线状崩塌陡坎。崩塌体通常破碎成碎块堆积于坡脚，形成锥状岩堆（崩塌倒石锥），图像呈粗糙点状，通过高分辨率的DEM数据甚至可辨识堆积碎块石的形态与粒径。在RRIM图中，可明显观察到崩塌后壁粗糙嶙峋的岩体结构普遍呈暗红色，坡体下方崩塌倒石锥图像上呈细点状，颗粒感小且均匀，颜色也较后壁浅，崩塌后壁与崩塌堆积体分界线十分明显。从图 5的崩塌示例图中可见坡面发育多个崩塌灾害。

滑坡灾害通常呈簸箕形或舌形的平面形态，个别滑坡可见滑坡壁、滑坡舌、滑坡台阶、封闭洼地等地貌特征。滑坡按照物质组成可划分为岩质滑坡与土质滑坡^[29]，其中土质滑坡后缘呈圈椅状，具有圆弧形滑面，后缘常见平行于滑坡壁的弧形拉张裂隙，前缘常见由于隆起形成的鼓张裂隙与滑坡体向两侧扩散形成的扇形裂隙；岩质滑坡边界通常由各组结构面控制，呈折线形，滑坡前缘常有局部滑塌，后缘常见由节理发育形成的拉张裂隙^[30]。崩塌与岩质滑坡的区别一般通过后壁的形态区分，滑坡的后壁形态往往呈弧形或折线形，并向后凹陷，而崩塌后壁主要受岩体结构面控制，常表现为陡直的破碎岩壁。

基于DEM生成的山体阴影等可视化图像中，可根据圈椅状后壁、山脊线是否连续、坡体的突然凹陷、中前缘地形隆起等地貌特征进行滑坡辨识。主要识别标志为前缘的较大规模滑塌变形与圈椅状后壁，以及后缘的拉裂缝。在高精度影像上前缘变形呈高亮度显示，后缘拉裂缝表现为线状或条带状阴影或高亮度影像特征，而滑坡后壁色调偏暗。在RRIM图中，滑坡后壁与周围坡体相比，亮度更暗，颜色更深，滑坡堆积体则颜色较浅，亮度变大，图像较滑坡后壁均匀。土质滑坡后壁呈弧形，堆积体可见高亮度条状图像，即堆积体上的地形隆起；岩质滑坡后壁呈折线形，可见控制滑坡的结构面，岩质滑坡堆积体的图像较土质滑坡堆积体粗糙，呈细小点状。由图 5滑坡示例可见滑坡后壁呈明显的圈椅状，滑坡体有明显下错，且下方堆积体呈高亮状态。

泥石流通常发育在地质陡峻、纵坡降大的山区，典型泥石流流域可划分为物源区、流通区与堆积区。泥石流按照流域形态可划分为沟道泥石流与坡面泥石流。其中，沟道泥石流流域呈扇形或狭长条形，沟谷地形，沟长坡缓，规模较大；坡面泥石流流域呈斗状，无明显流通区，沟短坡陡，规模较小。泥石流最典型的识别标志为泥石流的堆积区，常见于沟口，呈圆弧扇形，即泥石流堆积扇，质感细腻且色调较浅，扇体上常见侵蚀沟，表现为线状阴影。通常在沟口发现堆积扇即可说明沟内曾爆发过泥石流，沟内常可发现较多的崩塌、滑坡等灾害形成的堆积体，即泥石流物源区。在RRIM图中，泥石流堆积扇呈高亮度状态，堆积扇颜色一般呈浅红-白色，容易与周围山体区分开。沟口可见明显的弧形扇缘，扇体上可见线形深色侵蚀沟。如图 5中泥石流示例，沟道出口有明显的泥石流堆积扇，并且能清楚识别到扇体上的侵蚀沟。

首先使用SAGA GIS中Topographic Openness工具利用基于LiDAR点云数据生成的研究区DEM分别生成地形正开度与地形负开度栅格文件；再通过ArcGIS栅格表面工具栏中的坡度工具生成地形坡度栅格文件，然后在栅格计算器中根据式（1）计算得到脊谷指数I栅格文件。选择相应色带后，将脊谷指数I与地形坡度进行叠加，并根据研究区地形情况对两个图层的透明度、亮度、对比度进行调整，得到关于研究区的RRIM图（见图 6）。从图 6中可以清晰地看出区域内震前地质灾害的发生范围以及其源区与堆积区的分界线。

图  6  研究区RRIM图

Figure 6.  RRIM Map of the Study Area

采用人工目视解译方法对研究区震前地质灾害进行解译。人工目视解译是利用图像的影像特征与空间特征，与多种非遥感信息资料相结合，根据逻辑推理与相关地质信息分析出地质灾害的范围、位置的常见基础人工解译方法^[31]。解译的同时，根据灾害的特征，参考文献[18, 27, 32]的地质灾害分类方法，对震前地质灾害进行细致分类（见表 1），将地质灾害分为崩塌、岩质滑坡、土质滑坡、沟道泥石流、坡面泥石流5种类型。研究区震前地质灾害总面积约11.8 km²，共311处，解译结果见图 7。数量与规模上岩质崩塌以绝对优势占主导地位，这也符合本文研究区以灰岩为主的岩性情况。滑坡也以岩质滑坡为主，而崩塌泥石流与沟道泥石流数量与面积大致相同。

图  7  震前各类别地质灾害解译

Figure 7.  Interpretation of All Kinds of Pre-earthquake Geohazards

根据文献[18, 19, 22]对九寨沟景区震前灾害的研究，对比解译结果，可以验证解译的准确性。文献[18-19]在研究九寨沟景区地质灾害的研究时，受到数据精度与景区植被丰富、地势险峻的限制，识别出熊猫海及周边范围的地质灾害共计8处，文献[20]则是基于九寨沟景区泥石流治理工程，对九寨沟景区部分泥石流进行了研究，位于本文研究区的仅3条。而本文基于高分辨率去植被DEM与RRIM图，识别出311处震前地质灾害。从图 8中可以看出，基于RRIM图解译的311处灾害包括了文献[18-20]标记的灾害点，且类型相同，并对其范围进行了细化。

图  8  解译所得震前灾害与实地考察对比

Figure 8.  Interpretation of Pre-earthquake Geohazards and Comparison of Field Investigation

九寨沟震区地质灾害主要以岩质崩塌为主，受植被影响较小，因此未考虑植被因素，且震前灾害发育时间跨度较大，结合前人对地质灾害分布规律的研究^{[16-17, 33]}，本文只选取6个地形因子：高程、坡度、坡向、地形粗糙度、曲率以及LF（landform fuzzy）对震前研究区地质灾害的分布特征进行分析，生成的因子栅格图见图 9。其中，LF是根据简单的地貌列表对研究区内的DEM位置进行模糊分类，主要从相对位置和功能上进行分类，本文将地貌单元划分为13种类型，可通过SAGA GIS计算得到。其余因子均通过ArcGIS计算得到。

图  9  研究区因子栅格图

Figure 9.  Factor Grids of the Study Area

通过ArcGIS软件的重分类功能与空间统计分析功能将各因子进行分级，并计算每个因子各分级的面积及分级内灾害的面积，结果见图 10。可以看出，震前灾害集中发育于高程2 600~3 000 m范围内，即坡体中下部，且坡度在40°~45°时灾害集中发育，符合地质灾害发育普遍规律^[33]。研究区震前地震灾害发育的优势坡向在东与南东两个方向，且在地面粗糙度为1~1.4时地质灾害较发育，尤其是粗糙度为1~1.2时，震前灾害发育比较明显。曲率将研究区分为凹（< 0）、凸（> 0）、平（=0）3类区域，且凹值越小，区域凹陷程度越大，凸值越大，区域凸起程度越大。研究发现，曲率为－0.5~0与0.5~1范围内，即在轻微凹陷与一般凸起的区域，震前灾害较发育。同时，从图 10中可以看出，在13个地貌类型中，坡脚洼地和坡肩凸坡震前地质灾害集中发育。

图  10  因子分级统计图

Figure 10.  Factor Grading Statistics

本文利用高精度LiDAR点云数据，采用RRIM数据处理方法，避免了光源、阴影、植被等因素的干扰，提出了山区地质灾害的地貌识别标志与遥感解译方法。本文在九寨沟震区地质灾害集中发育的熊猫海区域进行了震前地质灾害解译，建立了震前地质灾害数据库，以更好地分析该区域地质灾害的时空演化规律，为预测该区域的未来地质灾害提供数据支撑。主要结论如下：

1）高分辨率LiDAR数据结合RRIM的方法在识别植被茂密山区地质灾害方面具有独特的优势，可精准反映地质灾害的地貌特征，有助于地质灾害准确识别与解译。

2）总结了山区常见地质灾害——崩塌、滑坡和泥石流的地貌解译特征：崩塌的主要识别标志为陡峭且可见破碎岩体结构的崩塌后壁与下方的锥形堆积体（崩塌倒石锥），堆积体在图像上表现为粗糙点状或粒状，在高精度图像上甚至可以辨识堆积碎块石的形态与粒径；滑坡的主要识别标志为滑坡前缘的较大规模滑塌变形与圈椅状后壁，在高精度图像上前缘变形呈高亮度显示，凹陷的后缘色调偏暗；泥石流的主要识别标志为泥石流堆积扇，常见于沟口，质感细腻且色调较浅，沟内常可发现较多的崩塌、滑坡等灾害形成的堆积体。

3）在九寨沟熊猫海区域，共解译出震前地质灾害311处，总面积约11.8 km²，并将其分成5类：崩塌、岩质滑坡、土质滑坡、沟道泥石流和坡面泥石流，其中以崩塌为主。震前地质灾害主要分布在高程2 600~3 000 m范围内，坡度40°~45°、坡向东和南东及地面粗糙度1~1.2的区域。

RRIM图凸显了地形变化，不仅对地质灾害的解译具有显著效果，也可应用于侵蚀沟与岩体结构面的识别等。

值得注意的是，RRIM图虽然能展示细微的地貌信息，但其缺陷在于缺少定量化的其他地形信息。为补充其他地形信息，可将RRIM图与其他地形参数图层进行叠加，如与等高线图层进行融合，可补全高程信息。同时，DEM精度越高，RRIM图展现的地貌信息越丰富；当DEM精度较低时，RRIM图无法放大细微地貌特征。而在进行图层的叠加时，脊谷指数图层与坡度图层的透明度、对比度、亮度等参数均可根据不同解译要求进行调整。