2008年5月12日14时28分四川西部龙门山构造带发生M_W7.9强烈地震，严重震灾波及成都、绵阳、德阳等广大地区，造成了以汶川、映秀为中心及其周边地域建筑物的严重破坏和人员的重大伤亡.地震破坏猛烈且集中，伴随有强烈的突发性，因此地震灾害严重威胁着我国经济发展和人民生命财产安全(滕吉文等，2008).考虑上述情况，利用地震波强地面运动数值模拟来分析地震危险性，并及时了解震情与震中烈度最大区域，可以减少地震造成的人员和经济损失，因此地震波场数值模拟的重要性日益凸显.基于波动方程理论的模拟方法有很多, 如有限差分法、伪谱法、有限元法、谱元法、边界元法等(朱多林和白超英，2011).一般地，影响地震强地面运动数值模拟结果的因素较为复杂，如横向地壳介质分层结构、复杂地形、地震震源输入模型等.多年来关于强地面运动数值模拟的研究包括了对地形、物性参数和震源模型等很多方面的讨论(王真理和李幼铭，1999；张怀等，2009；Douilly et al., 2015)，如，刘礼农等(2004)用单程波方法对非均匀介质与复杂地形模型中的弹性波传播进行计算；张怀等(2009)利用有限元数值模拟方法来模拟受地震激发的福州盆地地震波强地面运动，并考虑沉积层在盆地内部的放大效应，分析盆地结构对强地面运动的影响.Komatitsch和Tromp(2010)结合洛杉矶盆地地形及三维地震波速结构，利用谱元法模拟盆地的地面运动，分析了盆地结构对强地面运动特征的影响；孟令媛和史保平(2011)采用改进的复合震源强地面运动预测模型，根据动态复合震源模型和近断层区域强震模拟的计算方法建立断层破裂运动学模型，结果有利于实现对大震强地面运动的预测.赵宏阳和陈晓非(2017)根据辽宁省海城市在1975年2月4日发生的M_S7.3级地震，结合海城地区的地下速度结构与地形起伏建立震源模型，利用有限差分方法计算了模拟地震的波场传播过程，得到地震烈度分布图.

大量的震害调查、强震观测以及理论研究表明，盆-山结构造成的地形起伏对强地面运动模拟具有较大影响.而龙门山地区的复杂地形对于讨论地震动的影响以及建筑物抗震设防具有重要的理论和工程应用意义(王海云和谢礼立，2010；张建毅等，2012).我们知道，汶川地震位于青藏高原的东缘与四川盆地的西缘，影响区域包含了多类典型的构造单元，如龙门山、四川盆地等；并且，这片地区海拔变化极其剧烈，从四川盆地到龙门山地区约10 km的距离内，海拔从约500 m陡增至约5000 m(邓起东等，1994；许志琴等，2007)，剧烈起伏地形在地震传播过程中产生的影响将不可忽略.已有的汶川地震强地面运动模拟研究，如Bjerrum等(2010)以反演获取的汶川地震震源滑动模型作为输入源，模拟再现了汶川地震的强地面运动；尽管文中建立三种震源破裂对比模型，试图寻求最佳震源参数，却并未考虑地形在数值模拟中的影响.另外，过去多年汶川地震的震源滑动模型有许多更新研究成果，如Zhang等(2012)利用强震动波形数据反演表明，汶川地震过程中存在障碍体延迟破裂现象；如何在强地面运动的模拟中引入障碍体破裂、以及障碍体破裂的波形效应等都值得深入探讨.

图 1

Fig. 1

图 1 龙门山地区汶川地震构造图 红线表示龙门山地区三条主断裂，黑色五角星表示汶川地震震源位置，震源球显示汶川地震震源机制，黑色三角形为波形比对的台站位置. Fig. 1 Seismotectonic map of the Longmen Shan region Red lines are major active faults. Black five-pointed star represents the epicenter of the Wenchuan earthquake. Beach ball shows the mechanism of the Wenchuan earthquake. Stations marked with black triangle are used to illustrate the differences of waveform between the models.

本文基于反演获得的汶川地震震源破裂过程，构造震源模型，重点考虑龙门山地区特殊复杂的地质构造和汶川地震震源破裂过程的障碍体破裂特征；最后给出基于动力学模拟的强地面震动模拟结果及分析，讨论大震作用下复杂地形、运动学反演输入模型对模拟地震动的影响.

本文利用中国科学院力学研究所非连续介质力学与工程灾害联合实验室提供的GDEM软件进行问题解算(Wang et al., 2013；王杰等, 2013, 2015).依据微分方程，即弹性介质的波动方程(张怀等，2009)

(1)

求解复杂介质的波动传播问题，强地面运动特征受震源破裂过程、震源与场点之间复杂的地质条件和地形变化等因素的影响；基本物理方程包括：①弹性动力学方程；②断层滑动模式；③断层的几何形态.其中地震断层在破裂传播过程中，经过复杂介质与地形起伏的波动传播问题，可用弹性动力学方程描述(陈晓非，2005).弹性动力学方程的解包括体力、边界值和应力值的贡献.可以证明，边界值和应力值与体力是等效的.所以，不失一般性，可以将弹性动力学方程的解表示为

(2)

式中u_i(r, t)是t时刻r处的位移在i方向的分量；G_ij是格林函数，它表示t'时刻用于r'处j方向的单位脉冲集中在t时刻r处产生的位移在i方向的分量；f_j(r', t')是物理上真实的等效体力的总和，是t'时刻r'处j方向的力；V表示等效体力的空间.弹性动力学方程能帮助正确地"预测"地面震动情况，从而给出地表受大震影响程度的分布情况(朱耿尚，2014).

"设定地震"模拟结果的质量依赖于断层几何、震源破裂模型及介质参数的选取，通过选定的运动学反演输入模型确定，并直接影响最终模拟结果的正确性，参数选取如表 1所示.在此基础上进行的强地面震动模拟计算，不仅可以评估某一设定的未来地震的地震危险性；也可以"预测"已经发生的地震，如2008年汶川地震，所造成的地面强震动(陈晓非，2005).

表 1 Table 1

表 1 强地面运动模拟参数与物性参数设置 Table 1 Earthquake source parameters used in ground-motion simulations and material parameters 震源 103.4°E/31.06°N/15.6 km 断层长度 213 km 滑移速度 0.5 m·s-1 平均上升时间 1.2 s 拟合频率 0.8 Hz 杨氏模量 7×1010Pa 泊松比 0.28 密度 2.65 g·cm-3

表 1 强地面运动模拟参数与物性参数设置 Table 1 Earthquake source parameters used in ground-motion simulations and material parameters

利用数值计算模拟地震波及强地面运动等问题，因为涉及复杂介质的波动传播问题，需要对有限区域模型进行无反射边界条件的处理(张怀等，2009).本文在模型边界面施加无反射边界条件，主要是为了在动态计算中消除人工边界引起的虚假反射现象，设置黏性边界条件.这种边界条件适用于应力波入射角度较小的情况，特别适用于模型内部爆炸产生的应力波的吸收.

为更好理解地形因素对于强地面数值模拟结果的影响，本文建立具有地形影响因素的模型Ⅰ与无地形影响因素的模型Ⅱ，计算模型尺寸则一致，均为337 km长，290 km宽，深度则为80 km.为获得更精细的强地面运动变化分布，在断层带处对网格进行适当加密.在强地面运动模拟参数与物性参数设置等条件都相同的情况下(表 1)，模拟两类模型同震期破裂过程，得到强地面运动波形模拟结果、地表峰值速度分布图谱.

模型依据汶川地震运动学反演结果(Zhang et al., 2012)中子断层的划分、位错的分布、子断层的震源机制(走向/倾角/滑动角)设置，结合网格化模型，将整个震源模型分割成102个子断层，模型大小以及破裂面的几何构型如图 2所示，并依据反演结果中的震源位置在模型中设定初始破裂点.由于地震破裂在整个断层面上不是同时发生的，而是与空间和时间上的传播过程有关，所以本文从反演结果中反推出断层面上各子断层每一个质点滑动历史的时间导数，即粒子速度.粒子的速度为箱型函数，如图 3所示，x轴表示时间，y轴表示质点滑动速度，t₀表示质点滑动开始时间，t₁表示结束时间.利用基本的箱型函数拟合各子断层的位错分布、破裂量级以及破裂时间过程，并最终获得了汶川地震震源动力学破裂过程(图 4).

图 2

Fig. 2

图 2 汶川地震三维地形模型 (a)三维整体模型；(b)断层示意图；(c)三维地形模型顶视图. Fig. 2 Three-dimensional model of the Wenchuan earthquake (a) The overall three-dimensional model; (b) The diagram of faults; (c) The top view of model.

图 3

Fig. 3

图 3 粒子速度-箱型函数示意图 Fig. 3 Particle velocity-diagram of box function

图 4

Fig. 4

图 4 依据(Zhang et al., 2012)运动学反演模型模拟的汶川地震东西分量动态破裂过程 Fig. 4 Snapshot of east-west dynamic rupture propagation of the 2008 Wenchuan earthquake

可以清晰看到，汶川地震的震源动力学破裂过程遵循脉冲破裂模式；破裂传播主要是沿映秀—北川主断层往北东方向扩展；破裂前缘扩展遇到障碍时，不仅会越过障碍体继续传播，同时还会在回头方向发生微破裂.可以看到，汶川地震的主要能量释放来自于一个形状不规则的障碍体，刚开始的时候障碍体阻碍破裂前缘的扩展，但是随着应力扰动超过这个障碍体的破裂极限时，障碍体迅速发生破裂.利用反演的滑动模型作为输入模型，模拟重现了汶川地震的震源动力学破裂过程，与运动学反演结果存在异曲同工之处(Zhang et al., 2012).图 7—8展示基于汶川地震震源动力学破裂模型获得的强地面运动模拟结果.

图 7

Fig. 7

图 7 两种模型地表位移场东西向传播过程 (a)模型Ⅰ的结果; (b)模型Ⅱ的结果. Fig. 7 Snapshots of east-west ground displacement propagation for models (a) is the result of model Ⅰ, and (b) is the result of model Ⅱ.

图 8

Fig. 8

图 8 两种模型地表速度场东西向传播过程 (a)模型Ⅰ的结果; (b)模型Ⅱ的结果. Fig. 8 Snapshot of east-west ground velocity propagation for model Ⅰ and Model Ⅱ (a) is the result of model Ⅰ, and (b) is the result of model Ⅱ.

图 5、6均为基于汶川地震动力学破裂过程所模拟的强震动加速度波形数据结果.计算输出模拟结果时，参考了实际强震动台站位置，为将来与实际波形的比对提供依据(表 2).尽管影响强震动波形振幅的因素复杂，由于两类模型对比因素只有地形起伏，从中可以发现地形因素对强震动波形具有不可忽略的影响.通过定性比较两个模型的模拟结果，可以发现：

图 5

Fig. 5

图 5 龙门山地区强震动台站的加速度模拟波形比对结果 (a)东西向加速度比对结果; (b)南北向加速度比对结果; (c)上下加速度比对结果.红线表示模型Ⅰ(有地形影响)模拟的强地面震动数据；蓝线表示模型Ⅱ(无地形影响)模拟的强地面震动数据. Fig. 5 The comparison on simulated acceleration of strong motion stations in Longmen Shan (a) is the comparison of E-W component, (b) is the comparison of N-S component, and (c) is the comparison of U-D component. The red line represents the strong ground motion simulated by model Ⅰ. The blue line represents the strong ground motion simulated by model Ⅱ.

图 6

Fig. 6

图 6 四川盆地强震动台站的加速度模拟波形比对结果 (a)东西向加速度比对结果; (b)南北向加速度比对结果; (c)上下加速度比对结果.红线表示模型Ⅰ(有地形影响)模拟的强地面震动数据；蓝线表示模型Ⅱ(无地形影响)模拟的强地面震动数据. Fig. 6 The comparison on simulated acceleration of strong motion stations in Sichuan Basin (a) is the comparisonof E-W component, (b) is the comparison of N-S component, and (c) is the comparison of U-D component. The red line represents the strong ground motion simulated by model Ⅰ. The blue line represents the strong ground motion simulated by model Ⅱ.

表 2 Table 2

表 2 强震动台站位置 Table 2 Sites selected for low-frequency ground motion 经度(°) 纬度(°) 站名 断层距(km) 104.113 30.4483 CDZ 90.6 103.5442 30.5167 DXY 47.3 105.8325 32.0738 GYS 75.9 103.4107 30.154 PJD 66.2 103.3361 31.4767 LXM 47.9 103.4478 31.4642 LXT 38.8 103.6791 31.9469 MXD 57.7 104.5208 32.5158 PWM 45.4

表 2 强震动台站位置 Table 2 Sites selected for low-frequency ground motion

① 一般地，地形在地震动传播过程中具有显著放大作用，尤其在高程变化大的龙门山地区(图 5)；而在地形起伏变化不大的区域，如四川盆地内，则具有明显的抑制作用.

② 强震动台站的断层距对地形效应具有放大或抑制作用.距离断层破裂带越近，地形效应越明显；反之，距离越远，则越不明显.如强震动台DXY，由于断层距较小，使得地形的波形效应被显著放大，尽管其位于四川盆地内；而MXD则由于断层距较大，地形效应被抑制，尽管其位于龙门山山区地形起伏较大的区域.

图 7给出了两类模型地表位移场东西向运动传播过程.图 7左图 0~21 s时间内破裂传播中有一段空区，对应于不规则障碍体位置.21~29 s时间内，随着破裂前缘向北东方向扩展，由于应力扰动超过障碍体的破裂极限时，破裂回溯传播使障碍体迅速发生破裂，此时此地产生最大的地表位移，处于断层北-中段.这也是汶川地震震害的主要分布区域，相对于初始破裂点区域，波场更复杂且幅度更大.同时，从29 s起，破裂沿北东方向进一步传播.由于地形高差起伏的影响(模型Ⅰ)，龙门山地区的地表位移更多是集中分布于断层带附近；当把地形高差影响去除时(模型Ⅱ)，则地表位移分布和扩展范围更大，直至深入到四川盆地内部.

图 8给出两类模型地表速度场东西向的传播过程.16~32.5 s时间内，速度分布特征进一步显示了破裂回溯传播使障碍体快速破裂，并在此处显现出明显的速度高值分布区域.且在49 s之后，破裂传播至远离震源区域，此时速度变化出现高值异常，导致沿同震断层扩展方向的远场速度水平分量衰减趋势较慢，反映汶川地震破裂累积时长较长.从模型对比可以看出模型Ⅰ虽作用范围小但速度场数值大，是地形约束造成能量更集中的缘故.

图 9给出了两个模型地表峰值速度(Peak Ground Velocity，PGV)，图 9a为模型Ⅰ的结果，图 9b为模型Ⅱ的结果，取相同数值范围的颜色栏便于比较.综合来看，映秀—北川断层沿北东走向方向首先出现峰值速度高值区；另一个高值区出现在跨走向拐折区；而震中区域并非地表峰值速度高值区.对比真实地区示意图，也就表明映秀—北川断层的峰值速度特征在汶川、北川附近形成明显的高值分布区域；另外映秀—北川断层后半段相较平均值又出现一个相对高值区，大致对应青川附近，这些结果同野外地质调查大致相符.

图 9

Fig. 9

图 9 两种模型地表峰值速度分布图谱比对 (a)模型Ⅰ的结果; (b)模型Ⅱ的结果. Fig. 9 Simulated PGV distributions in the study area based on the three finite fault-slip models (a) is the result of model Ⅰ, and (b) is the result of model Ⅱ.

将两种模型结果对比，模型Ⅰ中映秀—北川断层的地表峰值速度分布高值区相较模型Ⅱ范围更广，数值更大.另外，模型Ⅰ中灌县—江油断层的高值分布区域处于前半段，而模型Ⅱ中该断层的高值分布区域却处于断层后半段，且不论分布广度还是数值都大于模型Ⅱ中映秀—北川断层上的高值分布区域，猜测可能与龙门山区域特殊地形息息相关.综上所述，模型Ⅰ的最大烈度区域位于汶川与北川附近；而模型Ⅱ的最大烈度区域位于灌县—江油断层的后半段，处于清平、安县附近.

本文基于反演获得的运动学模型，构造震源模型，选用震源附近区域复杂地形数据，通过建立具有地形影响因素的模型Ⅰ与无地形影响因素的模型Ⅱ，模拟两类模型同震破裂过程，计算强地面震动波形模拟结果、地表峰值速度分布图谱，来阐明大震作用下复杂地形因素对于强地面数值模拟结果的影响.且基于峰值速度分布特征，圈定峰值图谱及预测未来大震强地面运动特征分布，进一步推测未来严重震害区.

由于实际的复杂地形、地震波场传播的复杂性、以及一些特定的场地条件，地震的波动传播问题还需要更多的研究.本文通过固定其他变量，只考虑地形因素，定性比较有无地形两种模型，结论如下：

① 一般地，地形在地震动传播过程中具有显著放大作用；而在地形起伏变化不大的区域，则具有明显的抑制作用.

② 强震动台站的断层距对地形效应具有放大或抑制作用.距离断层破裂带越近，地形效应越明显；反之，距离越远，则越不明显.

③ 因为地形高差与障碍体的影响，地震造成的峰值可能出现在震中区域之外.考虑地形影响模型的地表峰值速度(PGV)区域位于汶川与北川附近；而未考虑地形影响模型的PGV区域位于灌县—江油断层的后半段，处于清平、安县附近.

在汶川地震波动传播问题中，近断层效应、方向性效应和地形效应对波场能量传播起到重要的控制作用(陈晓非，2005)，因此进一步的数值模拟工作需要更加准确的震源模型与更加详细的速度介质模型，并综合更多的参数来更好地模拟其强地面运动.

感谢由中国科学院力学研究所非连续介质力学与工程灾害联合实验室提供的GDEM软件.