2008年汶川地震触发滑坡量大于造山作用生长量

R.N.Parker A.L.Densmore N.J.Rosser M.de Michele Yong Li Runqiu Huang S.Whadcoat D.N.Petley

不言自明，地震通过反复的垂直位移产生地形（Avouac，2008），此外大地震也是滑坡的主要触发因素（Keefer，1994），在构造与地表过程之间驱动山带演化的竞争中起主导作用（Hovius et al，1997；Densmore et al，1998；Whipple，2009；Hovius et al，2000）。 最 近 的 研 究 （Keefer，1994；Malamud et al，2004；Larsen et al，2010；Guzzetti et al，2009）表明，滑坡能够产生持续的高侵蚀率（为1～10mm／a的尺度），这对我们了解山脉地形怎样产生提出了挑战：如果地震触发滑坡产生的可侵蚀沉积体积超过了在造山带同震产生的岩石体积增量，那么假如这种沉积物被其他侵蚀过程清除出造山带，造山带的体积和平均高程就肯定会减少。因此大地震在产生同震岩体上升和促进滑坡侵蚀（Hovius et al，2011）的相对作用对理解地壳平流与剥蚀之间的平衡至关重要。

2008年5月12日发生在中国四川省的汶川MW7.9地震，由于震级大，区域地形陡峭，而且同震滑坡发生范围广泛（Dai et al，2011；Sato and Harp，2009），对研究滑坡与造山带演化之间的关系相当理想。该地震发生在龙门山区域，而该区域之下岩性组合复杂，包括元古代花岗岩、一个古生界被动陆缘序列、较厚的三叠系至始新统（？）前陆盆地相地层和一个微暴露的弱固结新生界沉积物层（Burchfiel et al，1995）。龙门山中的断层形成于晚三叠世（Li et al，2003），并且直到第四纪都作为右旋逆冲斜滑断层仍持续活动（Densmore et al，2007）。在北川和彭灌断层上，该地震产生了大于10m的斜向逆冲走滑地表滑动（Shen et al，2009；Liu－Zeng et al，2009）（图1）。GPS、InSAR数据反演（Shen et al，2009）与实地观察（Liu－Zeng et al，2009）表明，右旋走滑和逆冲倾滑断层位移的尺度和比例沿破裂迹线均显著不同，映秀和北川附近是两个不同的滑动和矩释放集中区（图1）。

图1 汶川地震触发的同震地壳抬升和滑坡。黑色区域代表单体滑坡。粗黑色线条表示地表破裂迹线（Liu－Zeng et al，2009），五角星表示震中。灰色轮廓表示滑坡标绘使用的影像区域。背景是根据de Michele等（2010）修改的基于合成孔径雷达分析得到的同震地壳抬升场地。粗灰色线表示用于结果投影的与破裂平行的剖线

为约束滑坡的侵蚀，使用地震后30天之内的高分辨率卫星影像标绘了龙门山地区面积为13 800km2内的同震和震后短期发生的滑坡图（见方法一节）。我们对这个原始滑坡编目数据进行重新采样，转为滑坡密度Pls：

式中，Als代表一个选定窗口区域At内全部滑坡的面积（Meunier et al，2007）。应用公式（1）得到的Pls值从震中附近的＞60%（At＝1km2）到 低 海 拔 的 四 川 盆 地 的 0%（图1）。Pls值沿断裂也有显著变化，在映秀镇附近沿岷江峡谷为高值（图1），其次是东北部，尤其在主要的横向河谷地区。这部分（而非完全）反应了沿地表破裂走向的差异（Liu－Zeng et al，2009）。Dai等（2011）的研究给出了Pls值在不同岩性中的明显差异，以及Pls值与震源距离的复杂关系。假定滑坡的发生并不完全受同震形变的控制，就可能导致构造地壳抬升与滑坡侵蚀在样式和方量上的不匹配。

要理解汶川地震中构造与滑坡过程之间的平衡，需要一个将个体滑坡面积Ai转换为总体积Vls的标度关系：

图2 滑坡发生和同震位移沿走向的变化。所有数据以1km的宽度为间隔投影到平行于破裂的A—A′线（图1）。（a）每1km宽的条带滑坡总面积。（b）全球基岩滑坡标度关系式（Larsen et al，2010）应用于每1km宽条带内各个滑坡导出的滑坡体积；其他关系式显示出类似的模式。（c）每1km宽条带内的净同震地壳物质改变量（de Michele et al，2010）。（d）同震地壳物质改变量减去滑坡方量得到的净同震物质改变体积。（e）沿同震破裂的卫星影像覆盖的样本区域面积分布。滑坡面积和体积的局部极小值并不对应着小面积的样本区

式中，n表示滑坡的数量，标度参数α和γ是随位置和山坡过程（例如基岩或者浅层滑坡体）而变的常数。我们通过公式（2）使用已出版的标度参数（Larsen et al，2010；Guzzetti et al，2009）以及在该研究区中41个滑坡野外测量得到的参数进行计算。所得结果（表1）现示出惊人的一致，并对滑坡物质的可能方量有一级约束。以γ＝1.332±0.005，应用Larsen等（2010）的适合全部滑坡类型的全球最佳拟合关系式，得到Vls＝5.73＋0.41／－0.38km3。根据Larsen等（2010）的基岩滑坡全球最佳拟合关系式（γ＝1.35±0.01）和根据野外测量导出的关系式（γ＝1.388±0.087）均得到Vls≈9km3的类似值，而根据Guzzetti等（2009）的全球关系式得到Vls＝15.2＋2.0／－1.8km3。表1中的体积预测值是最小值，因为影像覆盖了大部分但并非全部的地表破裂区（见方法一节），但是与13 800km2滑坡标绘区内0.42～1.1m的空间平均剥蚀量一致。将这些估算值换算为滑坡侵蚀率需要知道北川断层上触发滑坡的大地震的复发周期，但这些复发周期由几个宽间距探槽点的有限测年数据（Ran et al，2010；Lin et al，2010）或应变累积推断速率（Shen et al，2009）约束较差。假设合理的地震复发间隔为2000～4000 年 （Shen et al，2009；Ran et al，2010），则得到仅由滑坡引起的长期空间平均侵蚀率为0.1～0.6mm／a，类似于在龙门山东部用宇宙成因核素分析的相似千年时间尺度上震前总侵蚀速率0.2～0.6mm／a（Ouimet et al，2009）。

表1 滑坡标度关系式和体积估计值

这些滑坡方量的估计结果可以与同震地壳抬升造山带的物质结果进行对比。de Michele等（2010）倒转升序和降序模式的合成孔径雷达（SAR）数据（见方法一节），获得了以约350m为间距跨区域三维地表位移矢量数据（图1）。我们根据这些数据（公式3）得到的滑坡分布区域内垂直分量上的净增方量Vt＝2.6±1.2km3。这小于所有滑坡的体积估计结果，而仅与一个标准误差相当（表1），意味着地震导致的龙门山地壳方量增加远小于滑坡释放的潜在方量（图2）。然而，这种直接比较有两个重要的问题。第一，合成孔径雷达数据获得于2006年11月至2008年8月，因此记录着同震与震后的滑坡，以及同震与震后的形变。虽然滑坡仅影响了整个区域（13 800km2）的4%，而其对Vt的影响则截然不同。此外，滑坡造成的地面形变的干扰导致了合成孔径雷达分析的局部不连贯，这些不连贯像素并不用于计算地表位移（de Michele et al，2010）。反演得到的地表位移大小和方向与实地观察密切一致（de Michele et al，2010；Liu－Zeng et al，2009），表明在造山带规模的尺度下，滑坡引起的地表变化对位移估计并没有太强的影响。第二，也是更重要的，估计的滑坡方量并不一定等同于侵蚀量，要转换为造山带尺度的侵蚀率，需要滑坡物质有效地被冲出造山带（Hovius et al，2011）。虽然在地震之前，有一些泥沙沿着龙门山主要河流的峡谷中存储，但是在整体上，该区域内分布的还是裸露的基岩山坡和小于100m厚的沉积体 （Ouimet et al，2009；Kirby et al，2003），表明同震滑坡物质在整个地震复发周期内很有可能被有效地移出山口。但是缺乏震前和震后的输沙量数据使我们无法去量化滑坡物质的移除率（Hovius et al，2011；Dadson et al，2004）。

因此，如果斜坡与河流输沙过程在下一个触发大量滑坡的大地震来临之前能将汶川地震滑坡物质全部移出汶川地震区，那么该地震将可能对造山带造成大的净方量损失。这个不平衡如何影响龙门山地形的生长呢？我们强调，我们的结果是滑坡侵蚀与构造作用之间竞争的瞬时量度，仅间接地适用于定义造山带的长期体积平衡（Densmore et al，1998）。有可能该山脉正如 Godard等（2009）所认为的那样地形正在衰变，侵蚀速率超过了岩石隆起速率，虽然这种模式仍然需要更为仔细的热年代学调查验证。第二种可能是，一些长期地壳隆升通过震间形变（Perfettini et al，2010）或余滑（Freed et al，2006；Hsu et al，2006）累积，尽管后一种机制通常贡献了同震位移的一小部分。换句话说，长期地壳隆升的重要部分可能发生在更频繁的较小或较深的地震，这些地震产生较低的地面加速度峰值（Orphal and Lahoud，1974），并且引发滑坡的数量会低得 多 （Keefer，1994； Malamud et al，2004）。在那种情况下，大型或浅源地震会主要造成由小型或深源地震形成的造山带地形的剥蚀，这样就使斜坡保持在一定的阈值梯度。为了支持这个观点，Ouimet（2010）指出龙门山的短期（1 000年）侵蚀速率为0.2～0.3mm／a，低于百万年时间尺度的侵蚀速率（0.5～0.7mm／a；Godard et al，2009），并认为，一些大地震的侵蚀速率能够赶上较长期的岩石隆起速率。在确定给定地震条件影响下的滑坡精确模式与方量时，气候条件也可能发挥作用；然而，如果我们估计的滑坡侵蚀速率与长期和短期的区域侵蚀速率之间在量级上一致，短期的气候随时间的变化就不太可能对物质平衡发挥重要作用。更可能的是，汶川地震中的地壳抬升与滑坡侵蚀之间的平衡是反常的，不能按多个地震周期推算。看来是，具有较大缩短分量的地震将会导致岩石方量的净增加，而走滑为主的地震事件，因为广泛的滑坡但岩石隆起有限而会引起物质的净减少。汶川地震中的右旋和逆冲滑动在不同断层段是高度分开的（Liu－Zeng et al，2009），并且这些断层段上地壳抬升与走滑的比率在不同地震之间也可能不同（Densmore et al，2007）。因此即使滑坡模式和总体积仍旧相同，相继地震之间这种比率的大差异也会预期在净体积平衡中有大的时间变化。总之，汶川地震涉及的构造与侵蚀体积之间明显而有争议的不匹配，说明还需要进一步了解大地震在背景区域侵蚀速率和造山带长期演化模式中的作用。

方法

滑坡 制图

我们使用EO－1和SPOT 5号卫星的影像开发了进行单体地震滑坡客观制图的半自动探测算法（见附件信息）。用密度阈值和20°的坡度阈值去除峡谷中的假滑坡区，从EO－1号卫星影像中提取滑坡区域。有些工作（Dai et al，2011）表明，坡度小于20°的区域很少发生滑坡。用20°坡度蒙片的非监督分类圈定SPOT 5号卫星影像上的滑坡区域。用一系列面向特征的过滤器去除由公路和田地产生的假滑坡区，并且对该地图进行目视检查与纠正。这得出的滑坡图总面积为13 800km2（图1），覆盖了225km地表破裂（Shen et al，2009；Liu－Zeng et al，2009）中的150km，因此这里计算出的滑坡总面积和体积是下限值。然而，与野外证据（Liu－Zeng et al，2009）、断层模型（Shen et al，2009）、合成孔径雷达分析结果（de Michele et al，2010），以及与通过卫星影像和航空照片手工编制的独立滑坡编目图对比（Dai et al，2011）表明，标绘的滑坡区覆盖了绝大部分同震滑移区，代表了此次地震主要影响区的典型例子。

同震方量估计

通过组合C波段与L波段的空间卫星装载合成孔径雷达波幅数据，de Michele（2010）得出了由汶川地震导致的三分量同震地表位移场。在这里，我们使用向上的或者垂直的分量去计算龙门山体积的净同震变化，忽视掉低起伏的四川盆地的高程改变（图1）。在滑坡分布图覆盖的龙门山地区（图1），我们计算了方量的净变化：

式中A代表格网面积，Ux代表每个格网的垂直位移，n代表所有栅格的数量，得出Vt＝2.6×109m3。位移与地面实测数据之间差异的标准差不适合作为Vt不确定性的数据反映，因为随机（不相关）误差可能会对标绘区域的总方量导致可忽略不计的净贡献。代之以，我们通过评估远离地震断裂的非形变区内Ux值统计变化的大小来估计Vt的不确定性。我们选择了四川盆地36km×36km的区域，其距离地表破裂45km，含有高水平的噪声（平均0m，标准差为1.5m）。在选定的区内，我们提取了30个剖面，每个36km长，并使用最小二乘法通过线性回归拟合每个剖面。因为y轴截距值影响每36km×1个像素面积下的方量估计，我们检查了每个剖面的y轴截距参数并计算了30个y轴截距参数与地面实测数据之间的均方根误差（RMSE）。计算得到的均方根误差为0.10m。当在整个研究区应用该误差值时，则等效于Vt的估计不确定性为1.2×109m3。

附件信息

1 滑坡制图方法

为了对地表破裂区进行快速的地震滑坡制图，使用光谱与地形相结合的模型开发了半自动滑坡探测技术。这些技术的完整描述见Parker（2010）。

1.1 全色影像分类算法（EO－1号卫星影像）

在为了研究而得到的10m分辨率全色EO－1号卫星影像中（附表1），滑坡的区域为亮色区域，暗色的植被已经去掉，岩土都暴露了出来。结果是图像有鲜明的反差，在滑坡影响区像素（数字值，或DN）值高，相邻山坡像素值很低。通过应用合适的像素强度阈值，基于每个影像的最佳可视化分类，滑坡与非滑坡区域提取为二元分类。然而，还有一些数字值高的地区包括城市区、公路和高沉积的河床。为了从分类中去除这些区域，使用90m分辨率的SRTM DEM数据将低于20°坡度的区域删除。这个最优的坡度掩蔽值通过对比从10°到40°（每1°一个间隔）掩蔽值得到的滑坡区和在原始影像中目视可鉴别的滑坡区用试错法估计得到。

1.2 多光谱分类算法（SPOT 5号卫星影像）

在5m和10m分辨率的多光谱SPOT 5号卫星影像（附表1）中，滑坡断崖显示出与周边植被覆盖斜坡有高水平的光谱反差。我们遵循Borghuis等（2007）的方法并应用基于最大似然方法的非监督分类将滑坡与非滑坡区域分开。用非监督分类来对影像进行光谱分类，然后基于每个影像的最佳目视结果选择滑坡类与非滑坡类。再次应用20°的坡度掩蔽（上述提到的最优的）消除类似于滑坡断崖光谱特征的显像区和河流。附图1展示了一个单体滑坡的人工与非监督分类结果的对比。

附表1 滑坡分布制图中用到的卫星影像

1.3 面向对象的滤波

所有覆盖区的分类掩蔽被组合起来并且转化为滑坡多边形。在一些地区，由陡峭斜坡上的公路段、耕地，或者沿由粗糙斜坡掩蔽分辨率产生的粗糙边界产生了额外的假抬升。基于分类特征的尺寸、形状和走向特性，用一系列的面向对象的过滤器消除这些解译错误。在这个过程中，主要考虑的过滤器是考虑滑坡多边形的方向与地形相对一致的方向。这些区域的长轴方向与斜坡的下坡方向一致的被保留下来，而那些长轴方向与斜坡下坡方向的夹角超过40°的则排除掉。相邻像素组合的总面积小于300m2（在10m分辨率的EO－1号卫星和SPOT 5号卫星影像的3个像素，在5m分辨率的SPOT 5号卫星影像的12个像素），起初被认为是滑坡的，在这里也被删除。最终，我们移除了圆比率大于0.7的多边形，因为这些主要是以耕地为主，而那些滑坡长宽比大于7的区域也被移除，因为那些主要是公路和河流。所有的这些过滤器的最优值均基于与在最初影像上目视鉴别滑坡区的试错法的比较予以选择。

最终，全图经过目视检查，并且进行了必要的手动纠正。在这个阶段，震前的陆地卫星7号增强式专题绘图仪＋ （Landsat 7 ETM＋）和陆地卫星5号（Landsat 5）影像用于对比并且删除地震前已经存在的滑坡。使用这个方法的基础是新形成的滑坡在植被覆盖区内产生了与陆地卫星影像上同样可见的清晰滑坡断崖。高程大于3 500m的区域植被覆盖非常少，因此新形成的滑坡无法识别出来。因此，这些区域，还有被云覆盖的区域，从最终的分类覆盖图中去掉了（图1）。

该滑坡分布图与穿越龙门山3个6×6km的确认样本区内人工绘出的滑坡进行了对比。滑坡自动制图的结果相对人工解译的滑坡结果面积低，在这3个区域之内的低估率为6.2%～22.7%。这个低估是由于错误的解译（一些滑坡区未识别出）和错误的删除（一些滑坡区面积识别有误）组合产生的结果。通过自动方法而非人工目视解译方法得到的区域面积为整个样本区的3.2%～5.6%，而被自动制图删除的区域等于人工目视解译区域的7.4%～12.5%。这些错误导致了净解译面自动识别与人工解译的重叠率为58.7%～66.2%，类似于Borghuis等（2007）应用相同的方法得到的结果，比率为53%～66%。作为总的样本面积标定的百分比，这些结果表明总的滑坡密度被低估了1.8%～9.3%。

因为滑坡制图结果依赖于光谱值的设定阈值，无论是在光谱强度（全色影像）还是光谱特征（多光谱影像），图像分辨率对单体滑坡的所得面积是一阶约束的。给定影像的像素是否在滑坡内包括取决于其数字值（DN值）；因此，圈定出亚像素级别的滑坡边界是不可能的。这个限制其他方法同样存在，比如使用航片来进行滑坡制图。

附图1 人工与自动滑坡制图结果的对比。背景影像是SPOT 5号卫星多光谱影像（5 m分辨率）。（a）由人工制图得到的北川县景家山／北川新北中学滑坡的边界。（b）应用坡度掩蔽和面向对象过滤器之前SPOT 5号卫星影像的非监督分类结果。黑色区域为滑坡区。（c）应用20°坡度掩蔽之后非监督分类的结果。黑色区域为滑坡。面向对象的过滤器尚未应用到分类结果

附图2 滑坡的概率密度与滑坡面积的关系。灰色圆圈为本次研究的滑坡，黑色正方形为1994年北岭地震触发的滑坡（Harp and Jibson，1996），白色圆圈为1999年集集地震触发的滑坡（Lin et al，2003；Khazai and Sitar，2003；Chen and Wan，2004）

2 滑坡分布样式和总体统计

总计我们在13 800km2的区域内得出了73 367处滑坡。由我们的制图技术导出的滑坡区域概率密度分布，与其他大地震触发滑坡数据集报道的非常类似 ［如1994年北岭地震（Harp and Jibson，1996）和1999年集集地震（Lin et al，2003；Khazai and Sitar，2003；Chen and Wan，2004）］。附图2说明，我们的数据可用类似的逆伽玛或者双帕雷托分布函数描述（Malamud，2004）。

我们得到的滑坡总数量比Dai等（2011）在41 750km2范围内使用人工目视解译得到的总数56 847处滑坡要多一些。要精确地对比这两个数据集比较困难，但是我们发现我们的总滑坡密度，用（∑Als）／Amap计算（其中Amap表示总研究区面积），为4.1%，而Dai（2011）中得到的为1.94%。这可以解释为我们的研究区更集中在近地表破裂的区域（图1），并不包括远离地表破裂的滑坡发生较少的区域（Malamud et al，2004）。

3 滑坡体积

在滑坡密度高的地区，汇聚到一起的滑坡特征很难或者不可能鉴别和划分出单独的滑坡，无论是使用人工目视解译还是自动或者半自动解译技术均如此。在这些区域，会出现滑坡数量被低估和单体滑坡面积被高估的情况。这些不确定性对滑坡的体积估计具有重大的影响，因为是滑坡面积与体积之间为非线性关系（Larsen et al，2010）。对含有多个滑坡的单一特征估算的体积会大于将这些滑坡单个制图时计算的总体积。因此，计算得到的连接在一起的滑坡的体积会是最大值。这种情况不仅出现在这次调查中，也出现在所有根据滑坡面积标度滑坡体积的研究中（Larsen et al，2010；Guzzetti et al，2009），尤其是在滑坡高密度分布区。如果没有独立的滑坡面积－频率关系数据，就很难评价这种面积高估对我们结果的影响。然而，我们使用Dai等（2011）的近似滑坡累积次数－面积关系去导出滑坡体积的一阶估计，理由是他们的人工制图技术可能比我们的半自动解译方法对面积的过高估计具有（尽管还未知）不同的敏感性。我们导出的总滑坡物质方量与基于我们的数据所得结果（表1）一致，在15%～20%范围内，表明我们证明的滑坡量与构造变动量之间的不匹配并不是由于我们制图方法的结果。此外，由我们的制图导出的滑坡面积概率密度分布（附图2），并未说明我们的数据集较之其他地震事件在大单体滑坡的面积上有大尺度的偏差（从而高估滑坡的方量），这也是如果我们的半自动解译技术对这一问题特别敏感就可以预料到的。