0 引 言

谱元法(Spectral Element Method，SEM)属于广义有限元法，是建立在波动方程的变分或弱形式理论基础上，求解偏微分方程的一种有效的数值计算方法(严珍珍等，2009).该法最早由Patera在求解Navier-Stokes方程数值解时提出(Patera，1984).它将有限元法和伪谱法相结合，不仅具有有限元法处理复杂结构的几何灵活性，还具有伪谱法高精度和快速收敛的特性.谱元法的基本思想是首先将计算区域划分成有限个单元，然后在每个单元上采用伪谱法，通过在每个单元中配置不均匀的分布节点，把单元的近似解表示成截断的正交多项式展开式，最后采用Galerkin方法求解得到全局的近似解，从而完成对区域内波传播过程的数值模拟(王秀明等，2007；彭海阔，2007；严珍珍等，2009；刘启方等，2013).

采用不同的正交多项式，就可以得到不同形式的谱元法，常见的有Chebyshev谱元法和Legendre谱元法.Priolo和Seriani(1991)首先在声波的传播过程中采用了Chebyshev谱元法，随后该法才被应用于地震波传播问题的数值模拟(Tessmer and Kosloff， 1994；Laurenzano et al., 2008).Maday和Patera(1989)最早将Lagrange插值基函数引入到谱元法中，并与GLL(Gauss-Lobatto-Legendre)积分相结合，将积分节点取为插值点，从而形成了Legendre谱元法.Komatitsch和Vilotte(1998)通过模拟地震波在复杂地质体中的传播过程，证明了Legendre谱元法在地震动模拟中的有效性.此后，科研工作者在此基础上开展了大量的研究工作，使谱元法不仅能用于非均匀、各向异性以及复杂地质构造区域二维或三维地震动场的模拟，还能用于全球地震波的传播模拟，并认为在达到同等精度的前提下，该法能够通过采用较稀疏单元划分，很好的解决有限元方法因单元内高阶插值带来的Runge现象，且易于实现并行计算(Komatitsch and Tromp， 1999; Priolo，2001; Komatitsch and Tromp， 2002; Komatitsch et al., 2004; Lee et al., 2008; Tromp et al., 2008；Stupazzini et al., 2009; 严珍珍等，2009; Smerzini and Villani， 2012).

与传统方法相比，谱元法是一种新的数值计算方法.虽然其在理论上还有待进一步完善，但由于结合了有限元法处理复杂结构的几何灵活性和伪谱法高精度、快速收敛的特性，目前已成为进行大尺度、复杂地质结构模型地震动数值模拟分析的重要工具.然而，根据目前掌握的资料，在中国采用谱元法进行地震动模拟分析的研究并不多，对其理论进行详细分析和总结的更少.与Chebyshev谱元法相比，Legendre谱元法由于将GLL积分节点取为插值点，生成的质量矩阵为对角矩阵，数值计算为简单的显示时间积分过程，具有计算实现过程简单、计算效率高等特点.因此，本文将以Legendre谱元法为例，首先从多个方面对其理论进行阐述，然后对其在国内外的应用现状进行分析和总结，最后采用该法简单分析了埋藏基岩地形对地表地震动的影响作用.本文对推动谱元法在地震动模拟中的应用具有重要的理论价值.

1 谱元法理论 1.1 弹性波动方程及其弱形式

假定有一非均匀各向异性的介质空间，如图 1所示.图中模型的区域为

分别为模型的吸收边界和自由边界，边界的外法线单位矢量为n；x _S为震源位置.若将地震波近似为弹性波，则 Ω 空间内的线弹性波动方程可表示为(安艺敬一和理查兹P G，1986)

该式称为波动方程的强形式，式中ρ为密度，u为位移矢量，σ为应力张量，f 为震源项，Δ 为梯度算子，点号为张量点积的运算符.位移矢量 u的初始条件为：

，其中 u · 为速度矢量.震源项f为一双力偶力源，一般可以表示为(Komatitsch and Tromp， 1999)

式中，M 为地震矩张量，δ(x－x _S)为震源的Dirac Delta分布，S(t)为震源时间函数.谱元法计算中，一般输入(2)式中的点源地震矩张量参数.对于大尺度的、更复杂的震源则可通过多个点源模型的叠加，用矩张量表达式来模拟(Komatitsch and Tromp，1999; 严珍珍等，2009).

在弹性各向异性介质中，应力张量 σ 可按Hooke定律表示如下：

式中，C 为刚度张量，在三维各向异性介质中有21个独立分量，在二维完全各向异性介质中有6个完全独立的分量； ε 为应变张量，ε = 1 2[Δu+(Δu)^T] ；T为矩阵转置符号.应力张量 σ 在自由边界δ Ω 上的外法线分量为 0，即

在谱元法中，为了对波动方程进行数值求解，需采用加权余量原理在方程(1)的两端同时点乘一个任意的测试函数 w，然后在整个计算域 Ω 内进行积分，即：

对(4)式右边第一项在模型域内和边界处进行分部积分，即可得到波动方程的弱形式(Komatitsch and Tromp， 1999; Peter et al.，2011; De Martin，2011)

式中，各参数的含义同上.

理论上，对任何的测试函数 w，波动方程的强形式和弱形式都是一样的.弱形式最明显的特征是它自然满足σ·n= 0 的边界条件.故采用弱形式不仅能很好的模拟波动在线弹性介质中的传播，而且对非弹性介质、表面地形和面波传播的模拟也更为精确(Komatitsch and Viotte， 1998; Komatitsch and Tromp， 1999).但是，弱形式的计算量较大，对计算机的内存空间有一定要求.

1.2 空间离散

为了求解方程(5)的弱形式解，与有限元法类似，首先将整个计算区域 Ω 分解成N_e个不重叠的单元，如图 2所示.若将第 e个单元记为Ω _e，则有：

同时，吸收边界 Γ 也被分解成了N_b个互不重叠的边界单元 Γ _e，即：

谱元法中，常将计算单元划分为结构网格，如二维模拟中的四边形单元和三维模拟中的六面体单元.这样不仅具有通过张量积构造基函数的方便性，还能保证模拟结果具有谱方法的高阶精度性，且易于实现并行运算(Komatitsch and Tromp， 1999; Komatitsch et al.，2001).但是，对于不能全划分为结构网格的复杂模型，也可以划分为非结构网格，如二维模拟中的三角形单元和三维模拟中的四面体单元.不同的是，此时需要经过一些特殊的处理，如通过集合投影方式构建Fekete内插节点进行精确近似时，就需要特殊的投影算子等( Taylor & Wingate，2000; Taylor et al.，2000).为了探讨采用非结构网格划分单元时的模拟精度，Komatitsch等(2001)分别采用结构网格、非结构网格和混合网格对均匀介质中弹性波的传播特性进行了数值模拟研究.研究结果表明，在相同单元尺寸下，非结构网格上的谱元法精度比结构网格上的低.

为了保证叙述的简洁性，文中以下部分的论述以二维谱元法为例，三维谱元法的情况依此类推.

在将二维模型域剖分成N_e个四边形单元后，为了便于计算，常引入一个正方形标准参考单元区域 Λ =[－1，1]×[－1，1]，其上任意一点 ψ = ξ，η 都满足－1≤ξ，η≤1的条件.其中，ψ 中的元素也常用ξ_i(i=1，2)来表示，并计

ξ₁=ξ，ξ₂=η .

现假定对任意单元 Ω _e，在单元局部坐标和标准参考单元间存在一个可逆的映射函数F_e: Λ  Ω _e，通过它可进行实际物理区域与参考单元之间的转换，从而使计算统一在参考单元中进行.因此，对于单元 Ω _e上的任意一点 x = x，y，在参考单元中都有一点 ψ = ξ，η 与其相对应，即有

式中，x _a= x ξ_a，η_a，a=1，…，n_a；n_a为定义四边形单元几何形状的节点总数，n_a=4或9，如图 3所示；ξ，η为标准参考单元中的局部坐标；N_a(ξ，η)为四边形单元中第a个节点的形函数.

在谱元法中，形函数N_a(ξ，η)常由1阶或2阶的Lagrange插值基函数构成.式(9)为n阶Lagrange插值基函数的表达式，即：

从(9)式中可以看出，对于任意一点lⁿ_k x 在任意一点x_l处的值满足关系式为

式中，δ_kl为克朗内克(Kronecker)符号.当k=l时，δ_kl=1；当k≠l时，δ_kl=0.

根据式(9)和(10)，很容易得出1阶Lagrange插值基函数的两个插值点：x₀=－1，x₁=1，其对应的多项式分别为：l¹₀ x = 1－x /2，l¹₁ x = 1+x /2；同理，2阶Lagrange插值基函数的三个插值点：x₀=－1，x₁=0，x₂=1，其对应的多项式分别为：l²₀ x =x x－1 /2，l²₁ x =1－x²，l²₂ x =x x+1 /2.这样一来，将这些不同的插值基函数进行不同的组合就能得到一系列不同形式的形函数N_a(ξ，η).

在二维谱元法中，形函数的个数与插值基函数的阶数是相关的，并满足如下关系式为

式中，n_a为形函数的个数；n为Lagrange插值基函数的阶数.

根据式(11)可知，图 3中4节点四边形单元的形函数可由1阶Lagrange插值基函数的多项式组合而成，即N₁(ξ，η)=l¹₀(ξ)l¹₀(η)，N₂(ξ，η)=l¹₁(ξ)l¹₀(η)，N₃(ξ，η)=l¹₀(ξ)l¹₁(η)，N₄(ξ，η)=l¹₁(ξ)l¹₁(η).同理，9节点四边形单元的形函数则可由2阶Lagrange插值基函数的多项式组合而成.

由于进行了映射变换，任意一个物理单元与标准参考单元间的面积关系可表示为

式中，J_e为单元的雅可比行列式，即

.为了求得J_e的值，需对(8)式中的 x 求偏导数，即：

从上式中可以看出，对雅可比行列式J_e的求值实际上就是对形函数N_a偏导数的求解，也就是对Lagrange多项式的组合求偏导数.因此，为了保证J_e值的存在，在单元划分好后，一定要对质量较差的网格进行平滑处理.只有确保了J_e值的存在，才能保证物理单元与标准参考单元间映射的唯一性和可逆性(Komatitsch and Tromp， 1999; Tromp et al.， 2008).

1.3 单元内插函数及积分

在完成区域剖分后，为了求解方程(5)的解，还需确定计算单元的具体函数表达式.在传统的有限单元法中，表征单元特性的函数表达式一般都选用低阶的Lagrange插值多项式做基函数，而在谱元法中则常选用高阶的Lagrange插值多项式做基函数，这是谱元法与有限单元法的主要区别.一般认为选用4-8阶Lagrange多项式作为基函数能很好地保证计算的精度和稳定性(Seriani and Priolo， 1994; Komatitsch and Vilotte， 1998).此外，在谱元法中，Lagrange多项式的插值节点常取为GLL积分节点.GLL积分节点定义为方程

的n+1个根，记作ξ_i∈ －1，1，i=0，1，…，n.式中，P′ n(ξ)为n阶Legendre多项式的一阶导数.图 4是5个4阶Lagrange插值多项式的GLL积分点分布图，从图中可以看出，这些点并不均匀分布，而是一种特殊的分布，Lagrange插值多项式在GLL积分点处的值为0或1.

GLL积分节点确定后，按照1.2节中的方法构造Lagrange插值多项式，得到每个单元的形函数，然后再按照一定的映射规则就可以得到每个单元 Ω _e上函数的多项式逼近函数.如对二维空间内的任意函数f(u或w 的分量)可近似的表示为

式中，l_α ξ 、l_β η 为Lagrange插值基函数，其性质如前所述；f^αβ为函数f在点 x(ξ_α，η_β)处的函数值，即：

同理，f(x(ξ，η))的梯度可近似地表示为

其中，f(x(ξ，η))为f(x(ξ，η))在参考单元上的偏微分，即：

式中， l′ _α(ξ)、 l′ _β(η) 为Lagrange插值基函数的微分；偏微分记号_i=_xi，i=1，2，且x₁=x，x₂=y； ψ / x 的值由对 x / ψ 求逆得到，x / ψ 的值由式(13)确定.

f(x(ξ，η))在单元 Ω _e上的数值积分，则可通过GLL积分进行计算，即：

式中，ω_α，ω_β，α，β=0，…，n为GLL积分的权系数，是大于零且与每个单元无关的常数(Komatitsch et al.， 2001).

1.4 全局聚合计算

根据以上对离散单元 Ω _e内插函数地分析，易得：

式中，ρ^αβ_e=ρ_e(x(ξ_α，η_β))；w^αβ_e=w_e(x(ξ_α，η_β))；C^αβ_e=C_e(x(ξ_α，η_β))；ε^αβ_e=ε_e(x(ξ_α，η_β))；其它符号含义同上.

对于震源项，根据式(2)并结合Dirac Delta分布的性质，可得到

若震源尺度较大且滑动机制复杂，则将(23)式变换为(Komatitsch and Tromp， 1999)

式中，S_s为有限断层面； m(x _s，t)为地震矩密度张量；其它符号含义同上.

最后，整理以上公式，并按照一定的顺序进行叠加，就可得到单元 Ω 在时间域上的二阶线性常微分方程组，即：

式中，M为全局质量矩阵； U 为全局位移向量；K为全局刚度矩阵；F为震源项.通过以上分析，可以发现与普通有限单元建立过程不同的是，谱单元的建立过程中不仅选择了特殊的插值基函数，还将GLL积分节点作为Lagrange插值多项式的插值节点，从而使式(25)中的全局质量矩阵变成对角矩阵，极大地减少了数值模拟过程中的计算时间，提高了计算效率(De Martin，2011).

1.5 时间离散

为了求解(25)式所示的线性方程组，谱元法中采用的方法有经典Newmark方法、4阶6步龙格库塔法(RK4-6)和4阶4步龙格库塔法(RK4).其中经典Newmark方法为最常采用的方法，该法是1959年由N.M.Newark发展的一类时间步进法，它基于两个公式(Newmark，1959)：

式中，U _n，U_n和U _n分别为与时间t_n对应的位移、速度和加速度矢量；Δt为时间步长，若把整个时间域 0，T 离散成N_T个时间点，则有Δt=T/N_T，t_n=nΔt，0≤n≤N_T；γ和β为Newmark参数，都为实数且有γ∈[0，1]，β∈[0，1/2]，它们对算法的稳定性、精度和耗散特性起控制作用.例如，当γ≠1/2时，Newmark法具有一阶精度；当γ=1/2时，对于任意的β值，都具有二阶精度；当β=0时，该法则是显式的(邢誉峰和郭静，2012；汪文帅等，2012).

在谱元法中，常取γ=，β=0来进行计算，此时式(26)、(27)变为

根据1.4节的内容，由于全局质量矩阵M是对角矩阵，故 U_n+1可由式(25)计算得出，即：

综合(28)、(29)、(30)三式，易知其计算过程是显示的，2阶精度的，且条件稳定.目前，经典Newmark方法已广泛应用在地震波传播的谱元法模拟中，这样不仅能极大节约计算的内存需求和时间，还特别易于实现并行运算(Komatitsch et al., 2000; Chaljub et al.， 2007).

1.6 算法的精度与稳定性

研究结果表明，采用谱元法进行数值计算时，为了使地震波场能得到合理准确的模拟，保证计算结果的精度，要求在每个波长内至少包含5个GLL积分点(Faccioli et al., 1997; Komatitsch and Vilotte， 1998).这实质上就是要求模型的单元尺寸、Lagrange多项式的阶数和传播介质的最短波长之间需满足一定的关系(Cupillard et al., 2012)，即：

式中，d为模型单元尺寸；n为Lagrange多项式的阶数，4≤n≤8；λ_min为传播介质的最短波长.若选Ricker子波为震源时间函数，则有λ_min=V_s/2.5f_p，其中V_s为介质中S波波速，f_p为Ricker子波的峰值频率.

此外，与其它数值计算方法一样，谱元法中也需要保证计算的稳定性.若采用Newmark积分方法，由于该法为条件稳定，为了保证整个计算过程的稳定，则要求时间步长Δt的取值必须满足CFL(Courant-Friedrichs-Lewy)条件(Chaljub et al., 2007; Cupillard et al., 2012)，即：

式中，Δx为两个相邻GLL积分点间的距离；V_p为介质中P波波速；CFCFL条件数，其取值详见表 1，该值具有随着波速比的增大而逐渐减小的趋势( DeBasabe and Sen， 2010).

2 谱元法的应用

Seriani 等(1992)最先将谱元法引入弹性波动方程的数值求解中，此后科研工作者在此基础上开展了大量的研究工作，使其不仅适用于非均匀、各向异性以及复杂地质构造区域的二维或三维地震动场模拟，还能用于全球地震波的传播模拟.例如，Komatitsch和Vilotte(1998)用谱元法对弹性波在复杂地质体中的传播进行了模拟，证明该法是一种进行二维或三维地震动模拟的有效方法，特别是能较准确的模拟瑞利波的传播、地震波的表面衍射以及不规则表面瑞利波与体波的转换等；Komatitsch和Tromp(1999)分别采用谱元法和离散波数法/ 反射率法模拟了三维层状模型的地震动响应，

结果表明无论是采用单点力震源，还是采用地震矩张量震源，谱元法对体波和面波的传播过程都能较好的模拟，且还能考虑粘弹性衰减、地形等因素对地震动的影响；Priolo(1999，2001)从网格划分、震源定义、模拟精度以及计算效率等方面详细讨论了谱元法在工程地震中的应用，并根据对意大利西西里岛东岸卡塔里亚市地震动反应的模拟结果，开展了地震危险性评价工作.Komatitsch 等(2004)用谱元法模拟了洛杉矶盆地在Hollywood(2001，M_W 4.2)地震和Yorba Linda(2002，M_W 4.2)地震作用下的三维地震动反应，结果表明该法能较好地模拟盆地模型中周期大于2 s、区域模型中周期大于6 s的地震动.Laurenzano等(2008)采用二维谱元法研究了Marche地区Treia和Cagli市在同一参考地震作用下的地震动反应，结果表明震源机制、地质构造以及场地条件是影响该区地震动反应大小的主要因素，并认为大面积深厚冲积物沉积层是导致Cagli市出现部分烈度异常现象的主要原因.Smerzini和Villani(2012)对2009年L’Aquila地震(M_W 6.3)的地震动反应进行了三维谱元法模拟，通过对模拟结果与实际震害分布情况的对比分析，认为该法能很好地模拟震中区频率不超过2.5 Hz的地震动，所合成的地震动波形与L’Aquila镇周围地区实际的强震动记录波形吻合良好.

此外，Lee 等(2008，2009)基于谱元法分别研究了台北盆地内和阳明山地区的地震动反应，结果表明盆地深度和浅层剪切波速是控制盆地内地震动放大的主要因素，山顶或山脊处的PGA值相对水平地表的变化一般在50%左右，而与山谷之间的差值最大可达2倍以上.严珍珍等(2009)根据弹性波理论，在考虑地表地形、地球介质衰减及地球椭率等特性的基础上，采用AK135理论地球模型，利用谱元法分别模拟了由人工点源和复合源激发的汶川大地震地震波全球传播过程，结果不仅显示了地震波在地球表面的传播形态，还证明了汶川大地震破裂过程为一个由多阶段破裂组成的复合破裂过程.周红(2010)采用谱元法模拟了点源激发下SH波的地表地震动反应，研究了地形几何形态、地震震源特性对地表地震动的影响，结果表明地震动的放大作用不仅受震源位置的控制，还受局部的地表张角控制，张角越小放大作用越强.胡元鑫等(2011)在同时考虑三维真实地形和介质衰减特性的条件下，结合ASTER DEM模型，采用谱元法对汶川地震的地形效应进行了模拟，结果较好地表现出了地震波波前形态在龙门山陡峻复杂地形影响下发生变形或扭曲，以及水平分量变化幅度大于竖向分量的特征，并且与平坦地形相比，震中附近大部分区域地形的PGA相对放大系数约为10%，部分区域地形对PGA分布及放大效应的影响强于震中距.刘启方等(2013)在考虑地形影响的基础上，采用谱元法模拟了施甸盆地在中等设定地震作用下的地震动反应，结果表明小型盆地内的地震动在这种地震作用下具有很大的空间差异，较深的小型凹陷和覆盖层厚度剧烈变化处是地震动强度增大的多发区，并认为盆地内地震动持续时间较基岩地震动明显增加，也是造成盆地震害严重的重要原因之一.

通过以上对谱元法在国内外应用现状的总结，可以看出谱元法的确是一种进行地震动数值模拟的有效方法.它不仅能模拟出各种地震波的传播过程，还能较好地反映出地形条件、场地土特性、土层结构以及断裂构造等多个因素对地震动的影响.

3 算 例

通过以上对谱元法理论和应用现状的介绍，接下来将采用该法对埋藏基岩凸起地形下地表地震动的变化情况进行简单分析.众所周知，研究埋藏基岩地形对地震动的影响，首先应研究同一震源发出的地震波在基岩面上不同部位的变化趋势，然后再研究地震波经过基岩面到达覆盖层表面后的变化趋势.因此，下面的分析中，将建立两个模型，一个为只有基岩层的基岩模型，另一个为在基岩模型上增加一覆盖层的组合模型.为了较好地体现出基岩地形对地震动的影响作用，所建模型的介质都为均匀且各向同性的.

3.1 模型建立

图 5为建立的基岩凸起地形模型，其中凸起基岩面的形状由标准椭圆方程x²/a²+z²/b²=1决定.若取椭圆地形宽度的一半为w=400 m，最大凸起高度为d=200 m，则基岩凸起地形的陡度为d/w=0.5.规定模型的左下角顶点为坐标原点，模型的横向(或水平向)为X方向，纵向(或垂直向)为Z方向，其中X方向取向右为正，Z方向取向上为正.

在图 5中，基岩模型为局部凸起的不规则模型，其横向尺寸为4000 m，纵向尺寸最小为1000 m，最大为1200 m；组合模型为规则的四边形模型，其横向尺寸为4000 m，纵向尺寸为1250 m.两个模型的表面都设置有49个横向间隔为50 m的计算点，每个计算点的具体位置如图中黑色小三角形所示，其中点17与点33之间的点位于基岩凸起界面或与凸起区相对应的地表面上.模型中各层介质的物理力学性质参数值和尺寸大小如图所示.

采用二维谱元法进行模拟，每个模型中的震源均设定为位于模型(300 m，300 m)处(图中五角星标记处)的双力偶点源.震源子波选用Ricker子波，取其峰值频率为f_p=6 Hz，截止频率为f_c=15 Hz.输入地震波的类型为P-SV波，模拟过程中同时记录X、Z方向上的加速度值.

为了满足计算的精度和稳定性要求，模型中划分单元的尺寸和时间步长都按1.6节中的要求选取，每个方向上划分单元的多少由该方向的最短尺寸确定.对于基岩模型，横向划分单元210个，纵向划分单元40个，总单元8400个；对于组合模型，横向划分单元210个，基岩层纵向划分单元40个，覆盖层纵向划分单元13个，总单元11130个.每个单元均为具有9个节点的四边形单元，并选用4阶Lagrange插值多项式为其函数表达式的基函数.两个模型设定一样的边界条件，即上边界为自由边界，左、右和下边界均为完美匹配层吸收边界，其中完美匹配层的厚度为8个单元层厚.每个模型的计算时间步长均为0.0003 s，计算总时长为15 s.

3.2 时程分析

图 6为基岩凸起地形模型在P-SV波作用下的加速度时程图.为了更好地体现出各计算点波形的变化特征，图中的波动传播时长均取为8 s，且各计算点的时程曲线按震源距的大小从下到上依次排列.

从图中可以看出，虽然每个图的加速度时程在整体上都具有随震源距的增大而逐渐减小的变化趋势，但在它们之间还是表现出了一些明显的差别.例如，基岩模型中，水平分量的加速度值在点17与点25之间明显偏小，而这种现象在垂直分量上却表现的不够明显，垂直分量上的加速度值除在点17附近有一定突变外，其余各点的变化都较为平缓.与基岩模型相比，组合模型中各点的加速度时程则表现出了不一样的变化特点，其内各点的波形不仅持续时间较长，还包含有较为丰富的面波成份.此外，从图中还可以看出，组合模型中水平、垂直两个方向上的加速度时程值，在与基岩凸起地形(点17与点33之间)相对应的区域都明显偏小，且偏小的范围较基岩凸起区大.

通过以上对基岩凸起地形模型中各计算点加速度时程图的定性分析，可以得出：

(1)基岩凸起地形的存在，对基岩面和地表面上的地震动强度都具有一定的影响，但影响程度各不相同.

(2)基岩凸起地形的存在，常会导致基岩面和地表面上部分计算点的地震动强度降低.其中，基岩面上地震动强度降低的点一般位于凸起基岩地形的上坡段，并且水平分量的变化幅度大于垂直分量；地表面上地震动强度降低的点则主要位于与基岩凸起地形相对应的区域，并且降低的范围较凸起区域大.

(3)覆盖层的存在，不仅延长了其上各计算点波形的持续时间，还使得它们的波形变得复杂，包含了较为丰富的面波成份.

3.3 PGA值分析

为了更好的反映出模型中各计算点地震动强度的差异，图 7给出了基岩凸起地形模型中各计算点PGA值的变化关系.可以看出，两个模型的PGA值在整体上都较好地表现出了随震源距增大而逐渐减少的变化趋势.

在基岩模型中，两个方向上的PGA值在基岩面起伏处(点17、点33附近区域)的变化幅度都较大，而在其它地方则变化不大.对于水平分量，PGA值曲线在接近基岩凸起处(点14、15、16)开始发生变化，到点17处出现陡降现象，尔后(点18、19、20)经过小幅振动后又逐渐增加，直到在点33附近区域又出现较大幅度的振动，最后才逐渐趋于平稳.而对于垂直分量，虽然其PGA值曲线在基岩面凸起处也有一定的变化，但变化的幅度不大.总之，从图中可明显看出垂直分量受基岩面起伏的影响比水平分量要小.在组合模型中，由于覆盖层的存在，其上各点PGA值的变化情况与基岩模型相比具有一定的差别.两个方向PGA值曲线的变化趋势基本相同，其中PGA值的突变区都在点17和点38附近，突变区之间则都为变化平缓且分布范围较广的PGA值减弱区.水平方向上PGA值减弱区的范围为点16与点38之间的区域，垂直方向上则为点17与点38之间的区域，两者之间相差不大，并且都大于基岩地形凸起区.

此外，根据组合模型与基岩模型之间各点PGA比值的变化关系也可以反映出基岩地形对地震动的影响.图 8为基岩凸起地形模型中各计算点PGA比值的变化关系曲线.从图中可以看出，PGA比值即放大系数，在水平和垂直两个方向上都具有相同的变化趋势，其中垂直方向的放大系数值整体上偏大.与图 7中出现PGA值减弱区一样，图 8中也出现了放大系数减弱区，并且两者的范围基本相等.这说明，在与基岩凸起地形相对应的区域内，地表面地震动不仅强度较小，其相对基岩面地震动的放大倍数也较小.

同理，通过以上对基岩凸起地形模型中各计算点PGA值变化关系的分析，可以发现同一震源发出的地震波，由于基岩地形的不同，其上接受波的形态和能量大小也不一样.主要表现为：

(1)基岩面上地震动强度的分布受基岩地形和地震波入射方向的控制.一般来说，地震波入射方向越接近基岩面的切线方向，其上接收的能量越小；越接近基岩面的法线方向，其上接收的能量越多.

(2)基岩地形的起伏对基岩面地震动水平分量的影响大于垂直分量.水平分量与垂直分量相比，不仅地震动强度大，在基岩面起伏处的变化幅度也大.

(3)基岩凸起地形的存在，使与凸起地形相对应区域的地表面地震动强度和放大系数都相对较小，形成地震动强度相对减弱区.与基岩凸起区相比，地震动强度相对减弱区的范围较大.

综上所述，可以看出对基岩凸起地形模型的谱元法研究，无论是时程分析还是PGA分析，都较好地体现出了基岩凸起地形对地表地震动的影响作用，即基岩凸起地形的存在，减少了地面地震动的强度，形成了地震动强度减弱区.这与Dezfulian和Seed(Dezfulian and Seed， 1970; 1971)采用有限元法得出的结论是一致的.究其原因，主要是由于凸起地形的存在，改变了地震波的射线路径，使其在基岩与覆盖层的交界面上出现多次折射、反射和透射现象，造成地震波能量在地表的发散，从而导致与凸起地形相对应区域处的地震动强度减小，形成地震动强度减弱区.

4 讨论与结论

通过本文的分析，可以看出谱元法是一种进行地震动数值模拟的有效方法.它不仅能模拟出各种地震波的传播过程，还能较好地反映出地形条件、场地土特性、土层结构以及断裂构造等多个因素对地震动的影响.然而，与国外相比，国内采用谱元法进行地震动反应分析的明显偏少，其发展是明显滞后的，并且对谱元法理论进行详细总结的也很少，这很不利于该法在地震动模拟中的推广应用.因此，本文首先从弹性波动方程及其弱形式、空间离散、单元内插函数及积分、全局聚合计算、时间离散以及算法的精度和稳定性等多个方面，对谱元法的理论进行了详细的分析和总结，然后对其在地震动模拟中的应用现状进行了归纳，最后通过一个简单算例的分析，在证明采用谱元法进行地震动反应分析有效性的同时，也较好地体现出了基岩凸起地形对地表地震动的影响作用.本文不仅对推广谱元法在国内的应用具有重要的理论意义，还对开展全面的地震动模拟工作具有重要的参考价值.

当前，在地震工程领域内，科研工作者已不再满足对一些局部、简单模型的数值模拟分析，而逐渐趋向于一些大尺度、复杂的区域性甚至全球性模型的数值模拟分析.因为，这不仅能加深对研究问题的认识，探究问题的本质规律，还能发现一些新的科学问题，从而得出更加完善、合理、科学的解释.例如，若要对沉积盆地、烈度异常区等一些尺度较大、地质构造和局部场地条件比较复杂地区的地震动反应进行全面的研究，都需要建立大尺度、复杂的数值模型，这不仅能保证模拟结果的相对真实性，还能减少由模型简化、尺寸缩小等带来的变异性和不确定性.然而，建立大尺度、复杂的模型就会出现计算量增大、精度和稳定性不易保证等问题，即使在计算机技术空前发展的今天，这也对所采用的数值计算方法提出了较高的要求.谱元法的提出正好解决了这些问题，它不仅具有高精度和快速收敛的特性，还易于实现并行运算，用其进行数值模拟计算能大大减少内存需求和计算时间.

总的来说，只要合理考虑介质的各向异性、地质构造背景、场地条件、地形地貌以及地震动衰减特征等多种影响因素，建立正确的计算模型，在对盆地效应、地形效应、烈度异常、地震危险性评价以及区域或全球地震波传播进行研究的过程中，采用谱元法进行数值模拟分析是可行的.