鲁甸地震的滑坡物质运移规律与地形特征

陈晓利 刘春国 传一健 兰 剑 魏延坤

1)中国地震局地质研究所，北京 100029 2)中国地震台网中心，北京 100045 3)北京大学数学科学学院，北京 100871

0 引言

滑坡在构造与侵蚀的相互作用过程中起着重要作用(Montgomery，2001；Larsenetal.，2010)。构造运动使山脉隆升、海拔高程及局部高差增大；侵蚀作用通过使山坡物质从高处向低处运移来限制山脉的高度，并调节或平衡了构造隆升和河流下切作用。在构造活动区，滑坡使山坡变陡的速率可能超过河流下切作用的影响，它引起的物质运移对局部地貌演化有着重要影响(Malamudetal.，2004；Fanetal.，2012)。与长时间缓慢的侵蚀过程不同，地震滑坡是瞬间的突发性事件，一次大地震能产生成百上千的滑坡，所造成的剥蚀、堆积等将使震区地貌发生显著改变。一些规模较大的滑坡会对局部地形产生重要影响，而地形条件是地震滑坡的一个重要控制因素(Hsü，1975；Ashfordetal.，1997；Chigiraetal.，2006；Meunieretal.，2008；Huangetal.，2009)。

地震滑坡涉及从山坡临界状态的形成、地震事件的触发到滑坡体运动、堆积和震后滑坡体演变的过程(Montgomery，2001)。每一次地震滑坡之后，山坡物质的分布状态就会重新进行调整：一方面，临界失稳状态的物质以 “滑坡”的方式运移进入一个相对稳定的状态；另一方面，物质运移的过程形成了一些新的不稳定区域，如滑坡体的后缘将成为潜在滑坡危险区。总体上，地震诱发的滑坡是对地表物质均衡状态的调整。

Keefer(2002)及Malamud等(2004)基于大量实例研究了地震滑坡引起的侵蚀速率。根据震级与诱发滑坡数量的关系，通过样本统计分析获得一定震级下相应的地震诱发滑坡总体积的计算公式，即地震事件造成的物质剥蚀量。这类关于地震滑坡剥蚀量的计算一般将整体的滑坡体积(运移物质)平均到一定范围，获得的是该区域内的平均剥蚀量。例如，汶川地震后，Parker等(2011)依据目前广泛应用的地震滑坡体积计算公式和震后SAR数据反映出的地形变化，计算了同震滑坡产生的物质剥蚀量和构造隆升形成的物质净增量，认为汶川地震在地震复发周期内将导致山体物质的减少。实际上，地震诱发滑坡的空间分布受到震级、断裂、河流、地形等的影响，使滑坡相对集中地分布在一些区域，并造成这些区域的地形变化。例如，汶川地震诱发了100多个规模>50i000m3的大型滑坡，这些滑坡受到发震断裂的控制，主要分布在地表破裂带附近5km范围内(许强等，2010)；由于滑坡规模大，地表破裂带附近的地貌受到较为明显的影响(王运生等，2008；祁生文等，2009；张永双等，2009；吴树仁等，2010)。除此之外，此次地震还形成了828个堰塞湖，这些堰塞湖规模不等，泄洪时间长短不一，无疑会影响河流附近的地貌形态(Fanetal.，2012)。

从汶川地震诱发的滑坡可以看出，发震断裂和河流附近是大型滑坡的聚集区，也是地表物质侵蚀的主要场所。中强地震所诱发滑坡的规模和分布范围相对于高震级地震都会有所降低，但在一定条件下也可能触发较大规模的滑坡，造成显著的地貌改观，如2014年鲁甸MS6.5地震触发的红石岩滑坡形成了高达120m的堰塞坝(Changetal.，2016)。研究表明，该次地震触发滑坡的空间分布没有明显受震中或断裂的控制，一些规模较大的滑坡发生在牛栏江及其支流等具有较陡边坡和较大地形高差的河谷两岸，显示出局部地形对滑坡分布的重要控制作用(Chenetal.，2015)。而鲁甸震区现有的地形条件是该区长期地貌演化的结果。

2014年鲁甸地震的震区位于青藏高原东南部边缘与华南块体之间，地貌上属于云贵高原的一部分。中新世以前，这里形成了较完整统一的准平原。由于经历了中新世末及以后的多次构造运动的强烈抬升，以及伴生的侵蚀、溶蚀等外动力作用，原来的准平原被雕塑、改造成现代的高山峡谷相间的山间凹地景观(Changetal.，2016)。河流下切使河岸变陡、岸坡岩体强度降低，为沿岸滑坡的发生创造了条件。鲁甸地震触发的红石岩大滑坡的数值模拟结果表明，地震动的触发起到了重要的作用(陈晓利等，2015)。然而，地震作用仅是大型滑坡的一种触发条件，并不是必要条件。2017年茂县新磨村滑坡(Suetal.，2017；Huetal.，2018；Wangetal.，2018)、2018年10月发生在金沙江上游地区的滑坡堵江等事件都不是由地震直接触发的，而是其他多种因素长期效应的累积结果，其中，地形可能是较重要的条件(Yangetal.，2019)。

本文以鲁甸地震影响区域为研究区，以网格为基本分析单元对研究区的地貌特点进行定量描述和统计分析，在此基础上开展地貌参数与鲁甸地震滑坡的相关性分析，揭示滑坡分布与地貌特点的耦合关系。与以往从滑坡体本身所处的地质地貌特征分析研究滑坡与影响因素之间的相互关系不同，本文着重从区域的地貌特点出发，通过地貌特征参数自相关性的研究与比较判断地表物质的分布特征，寻找其中的分异性单元(滑坡易发区)，从而认识滑坡发生的地貌特点。之后，进一步基于研究区地貌分布特点对未来研究区内可能发生物质运移的位置进行了预测。

1 数据与方法

1.1 数据

本文研究区是围绕鲁甸地震震中的一个矩形范围，几乎覆盖了此次地震诱发的全部滑坡(图1)，整体面积约为263km2。将研究区按照0.01°(经度)×0.1°(纬度)大小进行网格划分(每个网格的面积为1.104i1km2)，得到238个单元。每个单元具有4个基本属性，即表示滑坡严重程度的滑坡发生率和3个表示地形地貌特征的参数：平均坡度、平均高程和地形高差。各个参数的定义如下：

图1 鲁甸地震诱发的滑坡分布Fig.1 Distribution of landslides triggered by 2014 Ludian MS6.5 earthquake.数字表示单元内的滑坡发生面积比率(LAR)

(1)滑坡面积发生率

计算每个网格内的滑坡面积比(Landslide Area Ratio，LAR)，即滑坡面积占网格单元面积的百分比：LAR=滑坡面积/单元面积(%)，代表滑坡发生率。图1 为2014年鲁甸地震诱发滑坡面积的发生率分布。研究区共有111个网格单元受到影响，这些单元的滑坡发生率变化范围为0.03%～36.1%，对应的滑坡平面面积为300～398i544m2。其中20个单元的LAR>10% ，是受到鲁甸地震诱发滑坡影响较严重的地区，这些区域几乎都沿河谷两岸分布。如SE向距离震中约8km的龙泉河、沙坝河沿岸有2处的滑坡发生率达到36%；在震中西部震中距为4km处，沿牛栏江河谷也有一处的滑坡发生率达到26.4%。

(2)地形参数

基于研究区震前的SRTM 30m数字高程模型(DEM)(Shuttle Radar Topography Mission 30)计算得到每个网格单元的平均坡度(S)、平均高程(H)与地形高差(R)。这些地形参数所表达的含义如图2 所示。

图2 网格单元地形参数释义Fig.2 Topographic parameters of grid cell.S1—S5为通过DEM数据计算出的坡度；L为单元边长，R为局部地形高差，θ为期望坡度

每个网格单元的平均坡度(S)定义为该网格所包含的从DEM数据点中获得坡度的算数平均值，即

(1)

式中，n为网格内的坡度数据点的个数，Si为网格内各数据点的坡度。研究区的平均坡度分布如图3 所示。

图3 鲁甸地震诱发的滑坡分布与平均坡度Fig.3 Landslides and average slope angles in the Ludian seismic area.数字表示单元内的平均坡度

与平均坡度的计算类似，平均高程表示为

(2)

式中，n为网格内的高程数据点的个数，Hi为网格内各个点的高程。

局部地形高差(R)是单元格内最大高程与最小高程的差值，即

R=max(elevation)-min(elevation)

(3)

1.2 方法

一定区域范围内，地形条件对形成滑坡具有重要的作用。在对研究区进行网格划分的基础上，可以通过统计分析获得这些子区域地貌特征的相关性认识。如果某些子区域的地貌参数特点与区域整体地貌特点之间存在较大差异，则认为这些子区域的地貌特征与区域特征不符，是需要进行物质调整的单元，即易于发生滑坡的单元。

如何寻找可以衡量某一类地貌参数变化的物理量则需要根据不同的研究目的设定。杨晓平等(2019)通过分析南迦巴瓦地区的地形坡度和高程变异系数对活动断裂进行了识别。本文提出了一个新的网格单元坡度的度量值：期望坡度(Expected Slope Angle，简称ES)，表示数据的差异程度。对于每个单元，根据其地形高差(R)和单元的边长(L)计算出坡度角θ，

θ=tan-1(R/L)

(4)

该坡度角θ即为期望坡度ES。

从图2 及其计算方式可以看出，期望坡度(θ)和算数平均坡度(S)之间在物理意义和计算方法上均存在差异。期望坡度是以网格单元为度量尺度的描述地形特点的一种表示方法，仅与单元的局部地形高差和单元边长密切相关；而算数平均坡度与单元内坡度计算点本身的坡度数值及计算点的数量相关。可以认为，期望坡度是对一个分析单元实际地形的一种平滑。研究区期望坡度的分布见图4。

图4 鲁甸地震诱发的滑坡分布与期望坡度Fig.4 Landslides and expected slope angles in the Ludian seismic area.数字表示单元内的期望坡度

这里，引入统计学中 “变异系数”的概念来解释期望坡度和平均坡度的关系。变异系数又称 “标准差率”，是衡量资料中各观测值变异程度的另一个统计量。具体到本文所研究的问题，即通过局部地形高差获得的ES就是一个表示差异程度的量，而平均坡度又是一个表示平均数的量。从这种角度上，可以定义一个用期望坡度/平均坡度的指标，该指标能够在去除量纲的情况下分析不同地形的差异。

依据其统计学意义，可以用期望坡度和平均坡度这二者的相互关系来刻画地貌特征。具体而言，依据这二者之间的相对大小，可以将地貌的类型分为3类，如图5 所示(即后文图7 中A、B、C 3块区域的划分依据)。

图5 平均坡度和期望坡度的相互关系代表的不同地貌特征示意图Fig.5 Sketch showing topographic features derived from relationship between average and expected slope angles.θ为期望坡度角，A、B、C曲线为简化的地貌形态，其含义见正文。S1—S5是网格单元内数据点的分区，为便于理解，按照海拔由高到低的顺序排列数据点S1—S5

图6 坡度、高程、地形高差的关系Fig.6 Relationships between slope angles，elevations and topographic reliefs in the study area.

图7 研究区网格单元的平均坡度和期望坡度的关系及滑坡发生率Fig.7 Average slope angles(S)versus expected slope angles(ES)of grid cells and landslide occurrence rates denoted by landslide area ratios(LAR)in the study area.橙色点代表2014年鲁甸地震滑坡发生率LAR>10% 的单元，数字是其LAR值。A、B、C代表的区域地貌特征见正文解释

(1)A型地貌：地形单元的期望坡度>平均坡度，单元的局部地形高差较大，大部分地区的坡度较为缓和；

(2)B型地貌：地形单元的期望坡度<平均坡度，单元内的坡度普遍较陡；

(3)C型地貌：地形单元的期望坡度几乎与平均坡度相等，表明单元内的坡度分布或梯度变化较为均匀。

在以上地貌类型划分依据基础上，可进一步分析地貌类型与同震滑坡的关系。

2 结果

(1)总体地形特征

统计分析表明，研究区全部分析单元的平均坡度、地形高差和高程三者两两之间存在不同的相关性(图6)：1)平均坡度与局部地形高差呈正相关，局部地形高差值越大则平均坡度值就越大；2)平均坡度和高程呈负相关，高程越大，则平均坡度值越小；3)地形高差与高程呈负相关，随着高程值的逐渐增大，局部地形高差则变小。

Liu-Zeng等(2008)在研究青藏高原北部、中部、东南部地貌特征时提出了正地形和负地形的概念：如果高程与地形高差、坡度正相关，则为正地形；如果高程与地形高差、坡度负相关，则为负地形。其研究表明，青藏高原北部和中部以正地形为特征，东南部以负地形为特征，反映了水系侵蚀对于控制高原地貌的重要性。本文的研究区位于青藏高原东南部，高程与坡度和地形高差呈负相关，故属于负地形地貌，这与Liu-Zeng等(2008)的认识是一致的，说明研究区内河流侵蚀是影响地貌特征的重要因素。此次鲁甸地震诱发滑坡主要沿河谷分布，也是河流侵蚀作用的一种间接表现，说明局部地形条件对地震滑坡具有重要控制作用。

(2)滑坡分布特征

研究区各个单元的期望坡度和平均坡度的相关性分析结果表明，二者之间存在很好的正线性相关关系(图7)，其拟合公式为

y=0.940i5x+5.099i9 (R2=0.822i9)

(5)

R2>0.8，表明具有显著的回归效果。

然而，除了大量样本点沿拟合线密集分布外，还有一些样本点呈现出较为明显的离散性，接下来将对这部分样本点进行分析。为了分析这些离散程度较大单元的特点，首先需要确定哪些单元点是离散点。为此，采用

y=0.940i5x+5.099i9±2.5

(6)

作为上限和下限筛选出与回归拟合线相差较大的点(即ES与S的差异较大，位于图7 中的矩形框外)。在这样的筛选准则下，R2=0.822i9这一拟合结果能够保证62%的单元位于拟合线±2.5的范围内(即图7 中的矩形框内)，这部分样本点也就是该线性关系的主体数据。之后，再将剩余的38%的离散点根据其位于矩形框的左上侧(18%)与右下侧(20%)分为2类(图7 中矩形框外的A区和B区)。

在这样的划分原则下，A区共有43个单元，其中有7个滑坡面积 ≥10% 的单元；B区有49个单元，其中也有7个滑坡面积 ≥10% 的单元。因此，在研究区内20个LAR≥10% 的单元中，共有14个单元的期望坡度与平均坡度有较大的差异，位于图7 中矩形框外的A区和B区内，占总数的70%(14/20=70%)。该结果表明，不同的地貌环境与滑坡空间分布之间具有较为密切的关系，滑坡发生密度较大的单元其地貌特点与研究区整体的地貌特点具有较为明显的差异。因此，通过分析研究区整体的地貌特点，可以较好地识别出潜在的滑坡发生单元。我们的分析结果也认为鲁甸震区的地貌主要以C类型为主，即单元内以坡度分布较为均匀的地貌类型占主导。

3 讨论

对于一定范围的区域，高低起伏的地形是该区长期物质运移、分配的结果。不同地形具有不同的重力势能。在地貌演化的各阶段中，无论局部地形高差大或小、地形坡度陡或缓，每个阶段都存在一个相对稳定的状态，它反映了构造隆升与侵蚀之间的平衡(Montgomery，2002)。如果用地貌参数的相互关系来体现这个稳定状态，则应存在某种相对稳定的关系，反映出地貌侵蚀的主要机制。因此，在地貌参数数据的分析上，这种相对稳定的关系可表现为一定范围内各个单元都遵循某种规则，存在一定的自相关性。

本文研究区的高程与坡度和地形高差呈负相关，表明研究区内河流侵蚀是影响地貌特征的重要因素。在这样的地貌侵蚀机制下，尽管本研究所提出的期望坡度与平均坡度是从不同尺度对地形特点进行描述，但二者之间表现出较高的相关性(R2=0.822i9)，表明2种衡量地形变化程度的参数之间存在相对稳定的关系，代表一种占主导的地貌形态所具有的特征。因此，可利用期望坡度作为衡量标准来剔除一些 “特殊”的单元，这些特殊单元与Blöthe等(2015)提出的过量载荷地形具有一定的相似之处。Blöthe等(2015)认为，若一定区域范围内达到稳态的地貌形态，坡体应同时保持一个特征性的临界坡度，超过这个临界角度的地形就称为过量地形。过量地形代表了一种不稳定状态的地形，其在地形特征分析上应该表现出一定的 “特殊性”。本文中与稳定状态不符的单元，则可能是存在过量地形载荷分布的区域。这些单元将最终通过滑坡等物质运移方式进行地形高差或坡度的调整，从而靠近该区的分布规律(如拟合线)以达到与整体的和谐。

这种地形地貌分布上的差异性与滑坡空间分布及滑坡规模等之间的关系值得进行深入探讨。2017年8月8日发生于四川省阿坝州的九寨沟MS7.0地震触发了大量滑坡崩塌。与鲁甸地震诱发的滑坡相比，九寨沟地震所诱发滑坡的规模均较小，以浅表型滑坡为主。这2次地震发生处的地质构造背景和地貌环境也不相同。九寨沟地震发生在青藏高原中部巴颜喀拉次级块体的东缘地区，南北地震带中段的东昆仑断裂带东端和岷江断裂带的交会区域。地貌上，九寨沟地震震区位于青藏高原与四川盆地两大地貌单元过渡的、具有河流深切割特征的高山峡谷地带，地势由西向东快速降低，形成较陡的地形，最高海拔>4i500m，最低海拔约为1km，坡度一般>30°。采用本文提出的期望坡度和平均坡度的分析方法对九寨沟震区的地貌特征进行分析的结果如图8 所示，与鲁甸地震震区的分析结果具有显著不同。九寨沟地震震区的平均高程和平均坡度之间、平均高程与地形高差之间没有明显的相关性，难以用正地形或负地形进行描述；而期望坡度与平均坡度间的关系更为复杂。如果鲁甸地震震区的地貌参数特征可以代表以河流侵蚀为主的地貌类型，且河谷是大型滑坡主要的发育场所，那么九寨沟地震震区的地貌参数分析结果所代表的地貌类型与其众多的浅表性滑坡分布之间还需进行进一步研究。

图8 九寨沟地震震区平均高程与平均坡度和地形高差的关系(a)、平均坡度与期望坡度的关系(b)Fig.8 The relationship between average elevation and average slope and topographic relief(a)and the relationship between average slope and expected slope(b)in the Jiuzhaigou earthquake affected area.

4 结论

地貌演化具有阶段性和动态性的特点，一定区域内，各个地貌单元具有向相对稳定状态演化的趋势。滑坡作为山坡物质运移的方式之一，在地貌演化过程中起到了重要的作用。2014年鲁甸MS6.5地震诱发滑坡是该区在短时间内发生的地貌物质运移的集中调整过程。

这些地震触发的滑坡主要沿河流分布，表明河流侵蚀使河岸地形变陡、强度降低，形成发生物质运移的有利条件，也增加了地震滑坡的易发性。本文基于SRTM 30m数字高程模型(DEM)，利用网格方法对鲁甸地震区的地形特征进行定量分析，结果表明，研究区的高程与坡度、地形高差呈负相关，反映出显著的河流侵蚀效应；地形特征在期望坡度与平均坡度2个不同尺度下表现出很高的一致性，可能代表着地貌在稳定演化阶段所具有的一种特征。而与该特征不符的地貌单元，则是可能发生滑坡进行物质调整的区域。2014年鲁甸地震触发的大部分滑坡发生在期望坡度与平均坡度差异较大的区域，表明了区域地貌特点对滑坡总体上的控制作用。作为对比的九寨沟地震震区的地貌参数分析结果则表现出不同的特点。这种对地貌特点与滑坡分布及规模之间相互关系的认识有助于提高地震诱发滑坡预测的准确性及对潜在较大规模的滑坡场所进行早期识别。

致谢审稿人的意见和建议使本文得以完善，在此表示衷心感谢!