山脊线与坡度和峰值速度放大系数的相关性研究

蒋涵，周红，高孟潭

中国地震局地球物理研究所，北京 100081

1 引言

人们很早就认识到地表的起伏能够对地震动产生很大的影响.通常山顶或者山脊能造成地震动的振幅增加，而山谷则会造成地震动的振幅减小.地震灾害调查也发现在山顶上的建筑往往会遭受到比周围位于平地或山谷中的建筑更严重的损害.不仅如此，山体的顶部在地震发生后也更容易发生崩塌和滑坡（黄润秋和李为乐，2008；王运生等，2009）.因此起伏的地形对地震动的影响的研究，在地震学和工程学上都具有重要意义.

通过震害调查和数值模拟的研究，学者发现地形效应存在着一定规律.比如当入射的波长和起伏的几何特征接近的时候，地形效应会非常明显.而在入射波的波长较大的情况下，地形效应则可以被忽略（Ashford et al.，1997；Geli et al.，1988）.研究也发现，许多因素如入射波场，地下介质，以及局部地形几何形状等都会影响地形效应.Bouchon（1973）分别研究了P，SH和SV波垂直入射到二维峡谷形态的地形效应.他指出在山谷两侧的顶部，入射的SV波所产生的放大倍数最高.Bard（1982）设定了一个正弦形态的山脊，并用SH作为入射波.他发现向前传播的散射波振幅会随入射角度的增加而增大.Geli等（1988）总结了前人对二维的三角与弧形地形的研究成果，他从中得到了随山体高度与基底宽度的比值增加，对应的山顶和基底的地震动幅值之比也会增加的结论，作者还就介质起伏对地形效应的影响进行了探讨.周红等（2010）针对地形拐角角度变化对放大的影响展开研究，结果发现当台阶地形的张角小于180°时，地震动放大会随着角度的减小而增加.

以上的研究成果采用的是简单的二维模型，其计算结果与三维模型的结果相比存在着不小的差异.在三维地形效应中，地表的放大倍数更高，放大模式也更复杂.Liao等（1981）采用了三维截圆柱的模型，模拟得到的放大倍数是相应二维模型的两倍.Spudich等（1996）通过分析Northridge earthquake的强地震动记录，观测到垂直于山脊轴线方向的强烈的震动，其强度大大超过了二维模型的计算值，由此作者推测这是由三维的地形所导致的.该推测被随后三维模型的数值模拟（Bouchon et al.，1996）所证实.

上述几个三维数值模拟所研究的是孤立的山体，真实的地形则往往是多个山脊相连，因此局部的地形效应除了受山脊本身的影响以外，还会受到周围山脊的影响，比如不同走向的山脊的交汇处会有所谓的角落效应（Wang et al.，2006）.对这种真实条件下的三维地形已开展不少数值研究，比如San Andreas Fault地区（Ma et al.，2007），汶川地区（Zhang et al.，2008），台北地区（Lee et al.，2009a，2009b）的研究.学者通过将含真实起伏的地形和水平地形的模拟结果进行比较，从而得到该地区的地震动放大系数的分布特征.然而和简单地形模型相比，真实地貌特征和放大系数的关系还缺乏深入讨论.这当中有个重要原因就是三维地形特征的复杂性.

本文围绕着真实地貌特征与放大系数的特征展开研究工作.为了更好地分析研究地形特征，我们引入地理中的地形分析法，将山脊线和坡度作为地形特征因子来描述山体的特征.这样就将对三维地形特征和放大系数的相关性的研究，转换为对山脊线与坡度和放大系数的相关性的研究.不仅让分析的难度得以降低，也使得分析能够定量化.基于这个思路，本文选取了地形起伏剧烈的芦山地区作为研究对象，用谱元法模拟了由爆炸源产生的地形放大系数.通过分析地形因子和地震动放大系数分布特征，得到了它们之间存在的相关性.

2 芦山地区的三维地形效应模拟

2.1 数值方法

数值模拟是研究地形效应的有力工具，本文中所采用的数值方法为谱元法（SEM）.它最早应用于流体力学的计算中（Patera，1984）.后来由Komatitsch和 Vilotte（1998），Komatitsch和 Tromp（1999）引入到地震波场模拟计算中.谱元法结合了有限元的灵活划分网格特点和谱方法高精度的优势，因而在计算地震学中被广泛使用.

对于在地下介质中传播的地震波，其位移场满足的力学方程如下：

其中ρ是密度，u是位移向量，σ是应力张量，f则表示外力.

谱元法所采用的是（1）式的积分形式.我们将（1）乘以一个权函数w，并在整个计算区域进行分步积分，从而得到如下的积分方程：其中Ω是物理计算区域，Γ是区域的边界，t为边界上的力.

谱元法把计算区域离散为一个个单元.每个单元 由 Gauss－Lobatto－Ledendre（GLL）点 构 成 的Lagrange插值多项式表示.有限元中通常使用低阶次的单元，而在谱元法中则使用高阶次的单元.本文所使用的是6阶次的单元.将这些单元组装起来，并用U代表全局的位移，那么方程（2）可以重新写作

其中M代表质量矩阵，K代表刚度矩阵，F代表外力.

谱元法的积分点和单元节点重合，并且在数值积分的时候利用了GLL点的正交性，从而使得质量矩阵M成为对角的矩阵.M矩阵的对角化是谱元法的一个优势，它能够避免复杂的矩阵运算，提高了计算效率.下一节我们详细介绍基于谱元法的芦山地区的三维计算模型.

2.2 芦山地区的计算模型

我们选择的研究区域位于四川省境内的芦山地区（图1a）.它处在四川盆地和青藏高原东边缘之间.该地区附近地震活动频繁.从图1可以看出，芦山地区的地形复杂且山体起伏剧烈.起伏的地形会对传至地表的地震波产生很大的影响（周红和陈晓非，2006；Zhang et al.，2012）.为了研究该地区由地形造成的地震动的改变，我们考虑了两个计算模型.一个模型加入了芦山的起伏地形，地表的数据来源为gtopo30.另一个是对照模型，其地表为海拔统一为1km的水平地表.除去地表形态的区别，两个模型在其他参数上保持一致.计算区域为长方体，南北向为120km，东西向为110km，垂向上为50km.谱元法网格的单元长度为2.0km，最小的网格间隔为0.17km.

为了减少介质和震源复杂效应对地震动的影响，从而突出由地形起伏造成的地震动放大，我们选择相对简单的速度模型和震源模型.速度模型（见图1b）综合了前人正演（曾祥方等，2013）和反演（张冬丽等，2013）的数据，从地表往下P波速度依次为4.0km·s－1，6.1km·s－1和6.3km·s－1.速度的分界面为水平，目的是避免由介质界面起伏引起的地震动放大对结果的干扰；这样还能让速度界面形态在起伏地表模型和水平地表模型中保持一致.我们在地表下设置了一个爆炸点源作为震源.爆炸源的辐射花样不随着角度发生变化.这样的简化也是为了突出地形效应.震源的位置和2013年4月20日的芦山地震的位置相同.对震源的矩张量采用归一化的处理.震源时间函数则选择了Ricker子波来表示.

地形效应会随地震波的频谱变化而变化（周红和陈晓非，2006），本文选择了三种主频分别为0.8Hz，1.2Hz，1.5Hz的 Ricker子波来研究不同频谱的地形效应.主频1.5Hz的Ricker子波的有效频段最宽.其最大频率为4.0Hz，小于网格截止频率.模拟的整个地震波的传播模拟时长为50s.

3 两个模型的地表峰值速度

通过并行计算，我们模拟了三种主频下的地震波在芦山地区的传播过程，并得到了起伏地形和平层地形模型的地表速度时程.进一步，我们计算了各自的地表峰值速度（PGV）.图2为主频1.5Hz的两个模型在地表PGV的分布特征：在左图的平层模型中，其PGV的分布相对集中，呈圆形展布，越大的PGV值则越靠近震源；而在右图中，PGV的分布随着山体的起伏而呈现出不规则的复杂形状，在山脊上比周围地区的PGV幅值更大.为了更好地量化地形对于PGV的影响，我们将地表每个点的起伏地形的PGV和水平地形的PGV之比定义为PGV放大系数（PGV amplification factor）.图3为1.5Hz主频Ricker子波的PGV放大系数分布.从图中我们发现，在山脊或者山顶附近，通常对应了PGV放大系数大于1.这意味着山脊地形的加入，让地表的PGV幅值增加.而在山谷附近，PGV放大系数小于1.这表明山谷地形使得PGV幅值减小.PGV系数大于1的部分意味着地震动的破坏性增强了，是我们本文的研究重点.

然而PGV放大系数的分布特征是非常复杂的.这是由复杂的地形起伏所导致的.这种情况下，山体本身形状并不规则，并且许多山体并非孤立存在，它们组成复杂的山脊，从而导致了放大系数分布变得复杂.那么如何在复杂的地形上讨论地形特征与放大系数分布的关系成为了一个难题.

4 局部地形和PGV放大系数分布的关系

图1 （a）芦山地区的地形图，蓝色的矩形框为计算区域的地表投影，红色的五角星表示震中位置.（b）计算区域的三维图.南北长为120km，东西宽为110km，深度为50km.速度介质模型分为三层.分界面保持水平Fig.1 （a）The topography of the Lushan region.The computation area of the surface is indicated by the blue rectangle.The red star indicates the epicenter.（b）The 3Dcomputation area.The area is 110km in W－E direction，120km in N－S direction and 50km in depth.The wave－speed model has three horizontal layers

图2 主频1.5Hz的地表峰值速度PGV的分布情况：（a）平层模型；（b）带有地形的模型.PGV值由速度矢量的模长计算得到.由于震源矩张量经过了归一化，因此PGV幅值是无量纲的Fig.2 The PGV distribution of 1.5Hz dominant frequency of models （a）with flat surface and （b）with topographic surface.The PGV is obtained based on the length of velocity vector.Because the seismic moment tensor is normalized，the PGV value is dimensionless

为了更好地分析三维地形特征，我们采取了量化地形特征的方式，用一系列地形特征因子来描述地形特征（汤国安等，2003）.这些因子构成了地形变化起伏的骨架（李志林和朱庆，2001；Hengl and Reuter，2009）.只要使用相应的算法，就可以将这些地形因子从数字高程模型（DEM）中提取出来.从而将对地形和放大系数关系讨论化解为对这些的特征因子和放大系数的关系的讨论.不仅降低了讨论的难度，也能得到量化的分析结果.

4.1 山脊线和PGV放大系数的关系

山脊线组成了地形起伏变化的分界，它可以用来表示山脊的位置.而要从地形数据中获得山脊线，一种常用的算法为水系分析法（O′Callaghan and Mark，1984；Martz and Jong，1988）.它利用了山脊的分水性、山谷的汇水性的特点（汤国安等，2003）.该算法通过模拟水在重力作用下沿着地表坡面流动的过程，得到每个地表点的汇流累积值（flow accumulate）.接着对汇流累积值设定一个阈值，把超过阈值的点认作在山谷线上的点（Martz and Jong，1988）.同样地，我们只需先把地形的海拔值取负，从该“负地形”中便可获得山脊线.我们将芦山的地形数据输入到山脊线计算程序CATCH中，设定相应阈值，便得到了该地区的山脊线.

为了考察PGV放大系数在山脊线上的分布，我们从芦山地区中选出两个区域A和B（见图3中的黑框），将山脊线和三个主频的PGV放大系数的分布情况分别绘制在两个区域中（图4）.从图中我们观测发现当山脊顶部比较宽阔时，对应的山脊线呈现较大的网状；当山脊顶部比较狭窄时，山脊线则相应地变窄甚至成为线状.这表明利用算法得到的山脊线能够帮助我们识别出山脊的位置和形态.

图3 主频为1.5Hz的PGV放大系数分布图.暖色调表示PGV放大系数超过1，冷色调表示PGV放大系数小于1.图中的黑框表示被选择的两个研究区域位置.Fig.3 The PGV amplification factor of 1.5Hz peak frequency in the Lushan region.The warm color denotes that the factor is greater than 1.The cold color denotes that the factor is lower than 1.The black windows in the small picture on left side denote positions of the two selected areas.

从图4中可以看出，三种主频的PGV放大系数的幅值在山脊线周围的分布都是不均匀的.放大系数较大的部分（对应红色区域），均位于以下几种山脊线形态周围：第一种为山脊线出现分支时，通常是由较粗的山脊线分化为几条较细的山脊线，对应的山脊表现为由较粗的山脊分支为几个较窄细的山脊；第二种是山脊线发生汇聚的情况，这种情况下山脊线起始较细，延伸为较粗的、帽形的网状，对应的则多是孤立，或多面临空的山体.第三种情况中的山脊线形态较第一、二种情况复杂，它的宽度发生了变化，伴随有分支，或在延伸方向上发生了较大曲折.

上述分析说明PGV放大系数分布特征很复杂.放大系数较大的区域，往往是处于山脊线在延伸形态发生较大曲折或者宽度发生较大变化的地方.这些山脊形态也和地震引发次生灾害的分布有关系.黄润秋和李为乐（2008）调查了汶川地震灾害情况，文中指出，“单薄的山脊以及孤立或者多面临空的山体，对地震波最为敏感，具有显著的放大效应，这些部位的崩坍滑坡最为发育”.这里提到的“单薄的山脊”和我们第一种或第三种情况中的山脊线形态比较接近；而另一个“孤立或者多面临空的山体”则类似我们的第二种情况.

4.2 坡度和PGV放大系数的关系

本小节我们将研究另一个地形因子坡度和放大系数的相关性.坡度反映的是山体高程的最大变化率（李志林和朱庆，2001；Hengl and Reuter，2009），可以间接地表征地形的起伏形态和结构.地表某点的坡度定义如下：其中fx为x方向上的高程变化率.fy是y方向上的高程变化率.

图4 三种主频的PGV放大系数在区域A和B中山脊线上的分布情况.位于山脊线上的点用灰色的圆点表示，相邻的点用直线相连Fig.4 The ridge lines and PGV amplification distribution of three different peak frequency in the two selected areas A and B.Points on the ridge lines are connected by lines

在数字高程模型中，地表是由一个个离散的网格点所表示.坡度的计算是由某点和它周围的8个点所组成的窗口中进行的.提取坡度的算法一般是数值微分法或者局部曲面拟合法（刘敏等，2007）.本文中使用的是三阶差分算法（Zhou and Liu，2004），该算法也被ArcGIS采用.芦山地区的坡度计算结果见图5b.而图5a为1.5Hz的PGV放大系数图.通过观测，我们发现两者的分布有正相关性：在四个被选区域中，黑色区域的坡度值较大，对应的PGV放大系数较大；灰色区域内的坡度值较低，对应的放大系数值也较低.

进一步，我们用统计方法研究三个主频的PGV放大系数在坡度上的分布特征，我们使用了滑动平均坡度（moving average slope）.我们通过对滑动窗口内全部点的原始坡度求平均，把均值赋给中心点，从而计算出平均坡度.滑动窗口的长度为3km.为了统计出滑动平均坡度和PGV放大系数超过1部分的相关性，我们先从芦山地区的PGV放大系数中选出放大系数大于1的点；再将平均坡度划分为以0.05rad为间隔的区间，并求得各区间内的数据点的放大系数的均值.统计所用的数据分布和结果见图6.图6中的灰色圆点所构成的曲线定量地展示了滑动平均坡度和放大系数均值的关系.考虑到数据离散程度大，为了展示它的分布，我们用方框表示由滑动平均坡度和PGV放大系数组成的区间.方框内的颜色则表示了属于该方框内的数据点的个数.当颜色越接近红色，对应的点的数目越多；越接近蓝色表示数目越少.

图5 （a）1.5Hz的PGV放大系数平面图；（b）芦山地区的坡度分布图，坡度使用了弧度制rad表示Fig.5 （a）The distribution of PGV amplification factor generated by 1.5Hz dominant frequency；（b）The distribution of slope.The unit of the slope is radian

图6 滑动平均坡度和不同主频的PGV放大倍数的关系曲线.背景为统计所用的数据分布情况.对每个方框内的点的数目取了以10为底的对数，颜色表示不同的数目Fig.6 The correlation between the moving average slope and the PGV amplification factor of different dominant frequencies.The background denotes the distribution of the data for calculation.Numbers in small box is in base－10 logarithm.Different colors indicate different numbers

从三个主频的放大系数均值曲线中，我们发现它们都和平均坡度呈现出正相关性.其中1.5Hz的主频的曲线表现出的正相关性最强，随着平均坡度的增加，它的均值PGV放大系数从1.03单调递增到了1.38.而0.8Hz主频的曲线正相关性最弱.以上分析说明了不同的主频所对应的正相关性程度不同.

5 结论

本文用谱元法模拟了主频为0.8Hz，1.2Hz和1.5Hz的爆炸源在芦山地区所产生地震动.从起伏地表模型和水平地表的两个模型的峰值速度PGV中，我们得到了三个主频的峰值速度放大系数.芦山地区的地形复杂使得对它的几何特征与放大系数关系的分析变得困难.我们使用了量化地形特征的方式作为解决这个问题的办法.本文选择了山脊线和坡度这两个地形因子来描述芦山地形，用算法将它们从地形中提取出来.将对三维地形和放大系数关系的研究，转换到对山脊线和坡度与放大系数的关系研究上来.结果表明两个参数和放大系数都有很强的相关性.

我们分析了PGV放大系数在山脊线的分布，发现其幅值存在分布不均的特点.放大系数较大的区域通常为如下三种：（1）山脊线发生了分支，由较粗的山脊线分化为几条较细的山脊线.（2）山脊线发生了汇聚.（3）山脊线形态比（1）和（2）情况中的更复杂，它的宽度发生变化，并且伴随有分支，或在延伸方向上发生曲折.我们又分析了放大系数和坡度的关系，观测发现两者存在正相关.为了量化这种关系，我们计算了坡度的滑动平均，再将滑动平均坡度划分为等距间隔，统计了各自区间上的PGV放大系数均值，从而得到了滑动平均坡度和PGV放大系数均值的相关性曲线.结果表明三个主频的PGV放大系数和坡度呈正相关关系，不同的主频所对应的正相关性程度不同，其中1.5Hz的主频对应的正相关系最强.

以上的研究表明在大规模复杂的地形中，地形特征对于峰值速度放大系数的分布影响很大.考虑真实地形起伏对地震动的影响对地震危险性评估和防灾减灾非常重要.要深入地探究地形效应的背后机制，还需要我们使用更好的模型，模拟更多的区域，并且还应对其他的地形参数与放大系数的关系进行研究.

致谢感谢地球物理研究所二室的陈鲲研究员，九室的研究生贾晓辉对作者的研究给予的有益指导.感谢两位匿名审稿专家为完善本文所提出的宝贵意见.

Ashford S A，Sitar N，Lysmer J，et al.1997.Topographic effects on the seismic response of steep slopes.Bull.Seismol.Soc.Amer.，87（3）：701－709.

Bard P Y.1982.Diffracted waves and displacement field over twodimensional elevated topographies.Geophys.J.Int.，71（3）：731－760.

Bouchon M.1973.Effect of topography on surface motion.Bull.Seismol.Soc.Amer.，63（2）：615－632.

Bouchon M，Schultz C A，Toksöz M N.1996.Effect of threedimensional topography on seismic motion.J.Geophys.Res.：SolidEarth，101（B3）：5835－5846.

Geli L，Bard P Y，Jullien B.1988.The effect of topography on earthquake ground motion：a review and new results.Bull.Seismol.Soc.Amer.，78（1）：42－63.

Hengl T，Reuter H I.2009.Geomorphometry：concepts，software，applications（Vol.33）.Amsterdam：Elsevier.

Huang R Q，Li W L.2008.Research on development and distribution rules of geohazards induced by Wenchuan Earthquake on 12th May，2008.ChineseJournalofRockMechanicsandEngineering（in Chinese），27（12）：2585－2592.

Komatitsch D，Vilotte J P.1998.The spectral element method：an efficient tool to simulate the seismic response of 2Dand 3D geological structures.Bull.Seismol.Soc.Amer.，88（2）：368－392.

Komatitsch D，Tromp J.1999.Introduction to the spectral element method for three－dimensional seismic wave propagation.Geophys.J.Int.，139（3）：806－822.

Lee S J，Komatitsch D，Huang B S，et al.2009a.Effects of topography on seismic－wave propagation：an example from Northern Taiwan.Bull.Seismol.Soc.Amer.，99（1）：314－325.

Lee S J，Chan Y C，Komatitsch D，et al.2009b.Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral－element method and LiDAR DTM.Bull.Seismol.Soc.Amer.，99（2A）：681－693.

Li Z L，Zhu Q.2001.Digital Elevation Model（in Chinese）.Wuhang：Wuhang University Press.

Liao Z P，Yang B P，Yuan Y F.1981.Effect of three－dimensional topography on earthquake ground motion.Earthq.Eng.Eng.Vib.，1（1）：56－77.

Liu M，Tang G A，Wang C，et al.2007.Analysis of the slope uncertainty derived from DEMs.Geo－InformationScience（in Chinese），9（2）：65－69.

Ma S，Archuleta R J，Page M T.2007.Effects of large－scale surface topography on ground motions，as demonstrated by a study of the San Gabriel Mountains，Los Angeles，California.Bull.Seismol.Soc.Amer.，97（6）：2066－2079.

Martz L W，Jong E D.1988.CATCH：a FORTRAN program for measuring catchment area from digital elevation models.Computers＆Geosciences，14（5）：627－640.

O′Callaghan J F，Mark D M.1984.The extraction of drainage networks from digital elevation data.ComputerVision，Graphics，andImageProcessing，28（3）：323－344.

Patera A T.1984.A spectral element method for fluid dynamics：laminar flow in a channel expansion.JournalofComputational Physics，54（3）：468－488.

Spudich P，Hellweg M，Lee W H K.1996.Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge，California，earthquake：implications for mainshock motions.Bull.Seismol.Soc.Amer.，86（1B）：S193－S208.

Tang G A，Yang W Y，Yang X，et al.2003.Some key points in Terrain Variables Deriving from DEMs.ScienceofSurveying andMapping（in Chinese），28（1）：28－32.

Wang G Q，Tang G Q，Jackson C R，et al.2006.Strong ground motions observed at the UPSAR during the 2003M6.5San Simeon and 2004M6.0Parkfield，California，Earthquakes.Bull.Seismol.Soc.Amer.，96（4B）：S159－S182.

Wang Y S，Xu H B，Luo Y H，et al.2009.Study of formation conditions and toss motion program of high landslides induced by earthquake.ChineseJournalofRockMechanicsand Engineering（in Chinese），28（11）：2360－2368.

Zeng X F，Luo Y，Han L B，et al.2013.The LushanMs7.0 earthquake on 20April 2013：A high－angle thrust event.ChineseJ.Geophys.（in Chinese），56（4）：1418－1424，doi：10.6038／cjg20130437.

Zhang D L，Huang P，Zhang X B，et al.2013.Strong ground motion distribution and simulation based on finite fault model of Lushan 7.0earthquake on April 20，2013.Seismologyand Geology（in Chinese），35（2）：423－435.

Zhang W，Shen Y，Chen X F.2008.Numerical simulation of strong ground motion for theMs8.0Wenchuan earthquake of 12May 2008.ScienceinChinaSeriesD：EarthSciences，51（12）：1673－1682.

Zhang W，Zhang Z G，Chen X F.2012.Three－dimensional elastic wave numerical modelling in the presence of surface topography by a collocated－grid finite－difference method on curvilinear grids.Geophys.J.Int.，190（1）：358－378.

Zhou H，Chen X F.2006.A study on frequency responses of topography with different scales due to incident SH wave.ChineseJ.Geophys.（in Chinese），49（1）：205－211.

Zhou H，Gao M T，Yu Y X.2010.A study of topographical effect on SH waves.ProgressinGeophys.（in Chinese），25（3）：775－782.

Zhou Q M，Liu X J.2004.Analysis of errors of derived slope and aspect related to DEM data properties.Computers＆Geosciences，30（4）：369－378.

附中文参考文献

黄润秋，李为乐.2008.“5.12”汶川大地震触发地质灾害的发育分布规律研究.岩石力学与工程学报，27（12）：2585－2592.

李志林，朱庆.2001.数字高程模型.武汉：武汉大学出版社.

刘敏，汤国安，王春等.2007.DEM提取坡度信息的不确定性分析.地球信息科学，9（2）：65－69.

汤国安，杨玮莹，杨昕等.2003.对DEM地形定量因子挖掘中若干问题的探讨.测绘科学，28（1）：28－32.

王运生，徐鸿彪，罗永红等.2009.地震高位滑坡形成条件及抛射运动程式研究.岩石力学与工程学报，28（11）：2360－2368.

曾祥方，罗艳，韩立波等.2013.2013年4月20日四川芦山MS7.0地震：一个高角度逆冲地震.地球物理学报，56（4）：1418－1424，doi：10.6038／cjg20130437.

张冬丽，黄蓓，张献兵等.2013.基于有限断层模型的芦山“4·20”7.0级强烈地震强地面运动特征.地震地质，35（2）：423－435.

周红，陈晓非.2006.不同尺度地形的SH波频率域响应特征研究.地球物理学报，49（1）：205－211.

周红，高孟潭，俞言祥.2010.SH波地形效应特征的研究.地球物理学进展，25（3）：775－782.