基于谱元法的地震动放大效应研究——以青川县斜坡为例

万子轩

(成都理工大学地球物理学院，四川 成都 610059)

地形的放大效应通常和场地地形及近地表结构息息相关。许多研究(Davisetal，1973；Wongetal，1975；Celebi, 1985)表明，山顶、山脊和河谷处的地震动要明显大于平坦地形。而近地表结构例如山坡上的松软的土层和强风化层能够吸收大量地震波能量，从而导致地震动增强。在汶川大地震中，龙门山断裂所延伸的地方均发生了不同程度的山体滑坡现象。本文所研究的青川县区域，虽然在地震中未发生大型滑坡，但主城区乔庄镇附近的东山、狮子粱和桅杆梁等斜坡均发生了地表裂缝和岩体破裂，对主城区居民的生命财产安全造成了巨大威胁。斜坡的地震动力响应表现为山体震裂、崩塌、滑坡等地质灾害，地质灾害的发育程度与微地貌具有很好的对应关系。研究斜坡地形效应最直接有效的方式就是进行实地监测。然而因为地理条件和成本的约束，往往很难对一个区域内的各处斜坡进行全面的监测。在过去的十几年间，数值模拟方法被广泛地运用于地震动的模拟(Ma et al，2007；Komatitschetal, 2008; Maufroyetal, 2017)，但鲜有针对于小地形体的地震动研究。本文首次利用谱元法(SEM)模拟斜坡地震动及其放大效应，综合分析了该方法在斜坡灾害防治中的应用前景。

1 区域地质条件

青川县位于四川盆地北部边缘，川、甘、陕交界处，县城乔庄镇位于龙门山造山带与秦岭造山带的交汇部位，展布于乔庄河深切沟谷形成的近南北向狭长形二级阶地之上。县城北部属摩天岭构造带高中山区，南部属龙门山，为中高山深切峡谷地貌区。龙门山断裂带中的后山断裂平武—青川断裂恰好穿过这个区域(见图1)，断裂带总体上沿N70°E走向穿过县城主城区。汶川大地震中，映秀—北川断裂破裂带向东北延伸直插平武—青川断裂，一定程度上激活了平武—青川断裂。余震重定位结果(易桂喜等，2012；陈九辉等，2009)表明，破裂带北端是余震的多发地，且余震主要分布在青川区域西西南—北北东方向，震源类型从早期的走滑模式演变到后期的走滑兼逆冲模式。在汶川大地震中，该地虽未发生大型滑坡，但主城区乔庄镇附近的东山、狮子粱和桅杆梁均出现了山体裂缝和岩体破裂，对主城区居民的生命财产安全造成了巨大威胁。因此该地也成为防震减灾的重点区域之一。

东山斜坡位于县城城区东部，乔庄河左岸，坡脚高程约788 m，附近山峰高程大于 1 078 m，地形高差约为290 m，坡度一般为35°～60°，局部为陡崖。斜坡南侧为断层崖，断面产状出露N72°E/NW∠79°，其上盘出露碧口群(AnDbi)厚层状白云质灰岩(N70°E/NW∠45°) 逆冲于下盘茂县群(Sh)千枚状片岩之上(N75°E/SE∠21°)；狮子梁山位于县城城区西部，乔庄河右岸，整个山体东西长约1450 m，南北宽约400 m，最低高程约805 m，最高高程约962 m，高差150余米，斜坡总体坡度为30°～50°，狮子梁斜坡前缘坡角60°～50°，出露地层主要有古生界志留系茂县群黄坪组下段(Sh1)：粉砂质千枚岩夹变质砂岩板岩，元古界震旦系水晶组(Zs)白云质灰岩等。

图1 (左图)小方框为青川县所在地，线段代表断层分布(Xu et al，2009)(右图)为青川县高程图，黑色线段a-b和c-d分别代表东山—狮子粱剖面侧线和桅杆梁剖面测线。

2 基于谱元法的地震动模拟

2.1 谱元法简介

谱元法是一种融合了有限元和伪谱法的广义高阶有限元方法(Komatitsch et al，1998)。该方法能够很好解决有限差分法在高阶差分的情况下无法适应复杂地形和有限元法对于计算机性能限制的问题，非常适用于复杂地形下的地震动模拟。采用了Gauss-Lobatto-Legendre(GLL)积分以得到对角质量矩阵，大大提升了计算效率。其有效性得到充分验证。为保证模拟数值的稳定，该方法要求时间步长满足：

(1)

其中Δt为时间步长，c为库朗常数，取值在0.3～0.4，h为谱单元长度，v为介质的波度。同时，为了达到特定的模拟精度。所用网格的单元长度还需满足：

(2)

其中d为谱单元长度，λ为地震波波长，N为插值阶数，通常取5。因此上式意味着单元长度要小于一个地震波长，才能达到模拟相应频率地震波的要求。本文所有模拟结果均通过基于谱元法的三维数值模拟程序SPECFEM3D计算得出。

2.2 计算模型

利用无人机航空测距所得的高程数据构建高精度地形表面。高程数据水平分辨率约为6 m，能够很好的刻画出青川区域的地表特征(见图2)。在Cubit/Trelis软件中构建出包含地形表面的六面体网格模型(见图3)，采用了三倍加密层使得单元格长度在地表附近足够小，以满足较高频率的地震波模拟需要。整个模型体包含42 480个谱单元，位于模型表面的点之间的平均距离为3 m。Leeetal(2008)对台北阳明山地形效应的对比研究显示，层状速度结果和均匀速度结构模拟出的结果非常类似，且近震模拟中衰减非常微弱，可忽略不计。因此在缺少该区域三维速度信息的情况下，本模型采用均匀介质，根据PREM标准模型，P波和S波的波速分别取5.3 km/s和3.1 km/s。该模型对于11 Hz以下的地震波有着较好的模拟精度。因研究区域内无信息完备的地震可用，本文采用虚拟地震，震源位置位于研究区域的西侧(见图1)，走向为45°，倾角为90°，滑动角为0°，符合该区域断层的高倾角走滑特征和汶川地震余震的分布方位。地震震源深度设为6 km，震源时间函数选择为高斯函数。为了模拟较高频率的地震波，高斯函数的带宽设为0.1 s。根据实际的检测点位分布，在研究区域设置6个测点用来记录合成的加速度波形，分别为位于东山的Q1、Q2、Q3测点，位于狮子粱的Q4、Q5测点和位于桅杆粱的Q6测点。

图2 三维网格模型

图3 高精度数字高程构建的地形和东山、狮子粱斜坡真实地貌

2.3 地表峰值加速度放大效应

Specfem3D能够计算地震波传播至地表时三分量的加速度值，并自动计算地表最大加速度，取其作为该点的地表峰值加速度(PGA)。对于近震来说，地表面的峰值加速度主要受震源的辐射模式控制，从而不能很好地突出该地的地形效应。为此另外构建了不包含真实地形的模型体II。利用下式计算出地表峰值加速度：

(3)

PGAⅠ代表起伏地形计算出的地表峰值加速度，PGAⅡ代表平坦地形计算出的地表峰值加速度。

图4 青川县区域地表放大系数分布

区域的三维放大效应分布如图4所示，山体本身固有的地形效应在模拟中有很好的体现，即：斜坡的凸出部分和转折部分、山体的背处和山顶尖锥处存在地震动放大效应，而凹型的斜坡、斜坡的马鞍处存在地震动减小的情况。可以看出，东山斜坡、狮子山斜坡、桅杆梁斜坡对地震波均存在不同程度的放大效应，但其放大效应有所不同。沿着测线的PGA放大系数分布(见图5)也显示，狮子粱斜坡和东山斜坡的放大效应并不随高程线性变化。狮子粱山体狭长，放大效应主要沿着山体走向分布，放大系数随着其中狮子头处(Q4)尤为明显。东山山体较为浑厚，明显的放大效应出现在约960m高程的坡折处(Q1)，还可以看出，山脚位置(Q3)的放大系数微弱，可以作为基准台站的布设位置。而主城区的平坦地段反而产生了一定的放大效应，幅值约为0.2。这可能是由于在两侧山体散射的地震波在主城区重新汇聚，从而一定程度地加大了该处的地震动。桅杆粱的突出放大效应主要出现在山顶附近，局部峰值位于880 m高程处。可见，对于圆锥形的地形体，其放大效应整体上是随着高程增加而增大的，但突出放大效应不一定出现在山体最顶端部位。此外，桅杆粱放大效应整体沿着西西南—北北东走向分布，即地震波传播的背坡面一侧放大效应较迎坡面一侧更为明显。根据振动波在半无限空间中的传播理论，当振动波垂直入射到背波坡面时，将产生全反射，入射波和反射波相叠加，使坡面处的振动加强，从而导致背坡面灾害发育较应迎坡面更明显(许强等，2010)。这种现象同样出现在了A、B处山体。考虑到背坡面效应较为复杂，不能单由波传播的方向性解释，还受斜坡的陡峭程度、临空关系影响。因此该现象是否与背坡面直接相关还需进一步讨论。

图5 沿测线的PGA放大系数分布

2.4 地震动参数计算

地震动参数与地形效应有着密切的联系，是评价斜坡稳定性的重要手段。Anderson(2004)提出将常用的地震动参数如峰值加速度、峰值速度、加速度反应谱、阿里亚斯烈度、地震动持时等用于评价模拟结果在工程角度上具有一定意义。这里选取峰值加速度、阿里亚斯烈度和重要地震动持时这三个参数用于合成加速度记录。阿里亚斯烈度Ia计算公式如下：

(4)

其中，g表示重力加速度，a(t)为加速度记录，t0为整个加速度记录的持续时间。通常只计算水平合成向的Ia。阿里亚斯烈度定义为单位质量所吸收的地震波能量的总和，能够体现出地震动的强度及其持续时间。重要持时的定义基于阿里亚斯烈度，通常定义为阿里亚斯烈度达在5%Ia和95%Ia之间的时间间隔。

测点在三维地形模型I中的数值和该点在平坦表面模型II中的数值的比值即为该点的放大效应值。图6为各测点三分量的峰值加速度放大系数。可见，东山、狮子粱以水平南北向放大为主，东西向放大效应相对微弱。垂向的放大效应普遍小于水平向，且有随高程增加而增大的趋势。Q3因处于山脚位置，几乎没有放大效应。桅杆粱Q6点以垂向和南北向放大为主。图7给出了各测点阿里亚斯烈度和重要地震动持时的放大系数。作为表示地震动强度的参数，阿里亚斯烈度的放大系数与峰值加速度放大系数的变化趋势相近，但放大系数明显要更大一些。重要地震动持时并未显示出随高程变化的趋势，放大系数更为稳定，主要分布在1.2～1.5。模拟结果能够明显体现出地形使得地震波发生散射，导致地震动持续的时间显著变长，从而积累了更大的地震波能量。这也解释了峰值加速放大不明显的Q3点却在阿里亚斯烈度和地震动持时上有明显的放大效应。本文单独研究了表面基岩地形对强震参数的影响，并未考虑到场地条件。

图6 三分量峰值加速度放大系数

图7 阿里亚斯烈度、地震动持时放大系数

3 对比讨论

Luo (2014)在2009～2010年间对青川县乔庄镇各处斜坡进行了多组平硐检测，分别为位于东山斜坡的Q1、Q2、Q3点、狮子粱斜坡的Q4、Q5点和桅杆粱的Q6点，获得了宝贵的斜坡地震动数据。其对于十多个有感地震的检测分析结果与本文模拟有很大程度上的一致，但也有区别之处，主要在于：(1)根据现场调查，东山871 m突出山嘴，四面临空，动力效应十分强烈。但在模拟结果中，突出放大效应出现在了940 m高程的坡折处(Q1)。事实上，DEM模型中东山870 m坡体平滑，未见突出，可能是该点局布地形尺度过小，目前所用数据无法分辨。(2)桅杆粱的放大效应与监测结果一致， Q6点放大效应明显，山体的地形效应呈现出随高程增加而变大的线性趋势。(3)实际监测中，于狮子粱890 m高程处检测点显示该点位放大效应微弱，而模拟结果中该高程处放大效应强烈。这可能是由于受场地条件约束，该点并未布置在狮子粱主山头上。模拟结果也显示出，狮子粱主山头地形效应强烈，而侧翼的浑厚坡体地形效应微弱。(4)模拟计算得出的峰值加速度和持时的放大效应一般不超过2倍，阿里亚斯烈度不超过2.5倍。而在实际监测中，斜坡的放大系数往往会超过1.0，甚至可达4～5。这也表明，斜坡的放大效应不仅仅与其地形相关，还和场地条件(如松软土层、弱风化层、破碎岩层等)及近地表速度结构息息相关。前者是地震波数值模拟办法无法体现的，而后者则是今后的研究重点。值得说明的是，余震的方位角和震中距皆对放大效应产生一定影响，本文所对比的监测结果是该地最具有代表性的放大效应分布，因此模拟结果无法对应所有余震的分析结果，更不能和汶川大地震所产生的放大效应作直接对比。今后的研究中会进一步研究震源方位和震源机制对地形放大效应的影响。

4 结语

本文利用三维谱元法地震波正演程序SPECFEM3D模拟了虚拟震源作用下的强地面震动，旨在获得一个较为全面的青川区域地形放大效应分布。模拟结果显示，青川县乔庄处附近的东山、狮子粱、桅杆粱等斜坡均存在不同程度的放大效应，其中桅杆粱山顶附近放大效应明显，且主要沿山体走向分布，东北面斜坡放大效应微弱。东山和狮子粱斜坡中高程突出地形处放大效应显著，是地质灾害防治的重点。B、C处地形效应微弱，在汶川地震后的实地调查中也并未见到明显的地质灾害发生，是否存在可能的灾害风险还有待进一步模拟探究。基于强震参数的放大效应显示，东山、狮子粱、桅杆粱的峰值加速度以水平南北向为主，东西向和垂向微弱。阿里亚斯烈度在东山高高程、狮子粱中高程以上和桅杆粱顶部都显示出明显的放大，放大规律与峰值加速度类似。重要持时的放大效应普遍存在于各点，但没有明显的变化规律。综合前文的讨论，基于谱元法的数值模拟能够较好体现出地形体对地震波的放大作用。作为研究斜坡地形效应的一种新的技术手段，可以同场地检测、结构动力分析相结合，以期更加全面地评估斜坡的稳定性。