地震灾害直接与地震引起的地面运动有关． 及时掌握震中区地面运动的时空分布，无疑对了解灾情，进而采取救援措施至关重要． 理想的作法是利用高密度观测站与高速畅通的通讯，实时监测和传递地震现场的地面运动图像． 但实际上却存在难以克服的困难，高密度的台站分布几乎不可能实现，通讯高速畅通也不可能保证，尤其在地震发生的时候． 所以，权宜之计是利用遥远而稀疏的台站观测资料获得地震震源的时空过程，然后利用震源的时空过程再现地震现场的地面运动．

经过近30年的努力，利用远场观测资料获取地震震源的时空破裂过程已获得了长足的进展，震后5小时左右获取震源的破裂图像已成为可能(张勇，2008; 张勇等， 2008，2010． 震源的时空破裂图像在地震应急响应中发挥的重要作用也日益得到体现． 例如，2008年四川汶川M_S8.0地震(张勇等，2008)，2010年青海玉树M_S7.1地震(张勇等，2010)等． 然而，地震现场的地面运动图像，包括地面位移图像、 地面速度图像或地面加速度图像，对于应急决策更直观． 因此，利用地震模型快速计算近源的地面运动已成为现今社会发展的必然要求，也符合地震学学科发展的方向．

计算理论地震图的经典方法主要有基于射线理论的广义射线方法(Helmberger，1968； Gilbert，Helmberger，1972； Chapman，1974)和基于波动理论的反射率方法(Kennett，Kerry，1979; Bouchon，1981; Kennett， 1980，1983； Chen，1999)． 广义射线方法需要指定感兴趣的射线路径，而且随着射线路径的复杂化，计算成本也逐渐增加，因此多用于远场较低频的地震波计算． 由于需要指定具体的射线路径，所以，通常的计算结果并不是介质的完全响应． 比较而言，反射率方法却不需要指定具体的传播路径，计算结果自动包含介质的完全响应，因此更适合计算近断层的地面运动．

本文从分层均匀介质中地震波传播的基本理论出发(Kennett，1983)，参考已有的利用广义反射透射系数矩阵方法计算分层均匀介质中点源引起的地面运动的方法和程序(李旭，1993； 李旭，陈运泰，1996； 刘超，2011)，针对任意几何特征的有限震源，考虑震源机制在震源过程中的时空变化，顺应并行计算的发展方向，对计算流程进行优化，重新设计编写计算程序GRTMatSyn，并利用数值试验对该程序的可靠性和计算效率给出评价．

关于分层均匀介质中点源引起的地震波传播理论已经相当成熟(Kennett，Kerry，1979; Kennett， 1980，1983; Yao，Harkrider，1983； Chen，1999). 下面只简要给出与计算地面运动直接有关的核心公式．

如图 1所示，在分层均匀介质中，震源辐射的地震波总是可以分成两部分： 上行波(自震源向上传播的波)和下行波(自震源向下传播的波)． 为了描述问题的方便，在震源深度位置人为设置一个界面(这个界面其实不存在)，使其平行于其它介质界面，并以经过震源垂直分层界面向下的方向为柱坐标的中心轴建立柱坐标系．

图 1 分层均匀介质及坐标构架(Kennett，1983) z0表示自由表面，zL表示介质底界面，zS表示震源处的虚界面． U和D分别表示上行波场和下行波场 Fig. 1 Stratified media and coordinate frame(Kennett，1983) z0 denotes free surface，zL，bottom interface， and zS，imaginary interface of source. U and D represent up-going and down-going wave-fields，respectively

在柱坐标系中，分层均匀介质中点源引起的地面运动的面谐矢量位移解w₀可以表示为(Kennett，1983； 李旭，1993； 李旭，陈运泰，1996)

图 2 广义反射透射系数矩阵(R0SD，T 0SU，RSLD和R0SU)几何解释 z0表示自由表面，zL表示介质底界面，zS表示震源界面 Fig. 2 Geometrical illustration of the generalized reflection and transmission matrices R0SD，T 0SU，RSLD and R0SU z0 denotes free surface，zL，bottom interface， and zS，imaginary interface of source

由式(1)可以看出，面谐矢量位移解w₀依赖于自由表面的接收系数矩阵、 反射系数矩阵、 分层介质中上行波和下行波的广义反透射系数，以及上行波与下行波地震波矢量间断． 如果通过这些量得到面谐矢量位移解w₀，那么很容易通过Fourier-Hankel反变换得到时空域点源引起的地面运动解． 可见，计算地面运动的关键在于计算面谐矢量位移解w₀，而计算面谐矢量位移解w₀的关键则在于计算广义反透射系数矩阵和表征震源作用的地震波矢量间断． 因此，提高计算有限震源引起的地面运动效率的途径之一是优化广义反透射系数矩阵和表示震源作用的地震波矢量间断的计算流程．

如图 3所示，用界面A和界面B表示相邻的两个界面，用界面B和C表示任意多层介质之间的界面，用 R D(z^-_A，z⁺_A)与 T D(z^-_A，z⁺_A)分别表示地震波矢量通过界面A的下行波的反射系数和透射系数矩阵，用 R U(z^-_A，z⁺_A)与 T U(z^-_A，z⁺_A)分别表示地震波矢量通过界面A的上行波的反射系数和透射系数矩阵. 具体某界面的反透射系数的计算比较简单(详见附录B)． 相应地，界面B与C之间的广义反射透射系数矩阵用 RBC_D，T BC_D，RBC_U和T BC_U来表示，那么z^-_A与z^-_C之间的广义反射透射系数矩阵 RAC_D，T AC_D，RAC_U和T AC_U由下式确定.

图 3 地震波矢量在分层介质中传播的几何解释 以水平实线表示文中考虑的3个分界面zA，zB和zC； 虚线表示任意界面； 上标“-”与“+”分别表示分界面的上下边界； 带箭头的实线表示地震波矢量 Fig. 3 Geometrical illustration of the wave vectors in stratified media Solid horizontal lines represent the three interfaces zA，zB， and zC， involved in the text. Discontinued lines denote other interfaces， and superscript “-” and “+” mean upper and lower boundaries of interface. Solid line with arrow represents wave vector

在式(1)中，震源作用以地震波矢量间断∑_U(z_S)和∑_D(z_S)来描述(Kennett，1983)，下列各式把描述一般点源的地震矩张量和地震波矢量间断联系起来．

对于P-SV波

对于SH波

1)当m=0时

2)当m=±1时

3)当m=±2时

有限震源引起的地面运动，总是可以通过多个点源引起的地面运动的适当叠加得到(姚振兴，郑天愉，1985)． 在计算有限震源的地面运动时，我们必须考虑有限震源的一般性． 因为，实际地震的震源几何形状可能是多样的，地震过程中震源机制可能随时间和空间发生变化． 考虑到几何形状的一般性，我们并不把震源局限在一个平面内； 考虑到震源机制的变化，我们使用地震矩张量表示震源．

由式(2)和式(3)可以看出，介质内部任意层之间的广义反透射系数依赖于其介质结构． 如果它们的几何分层和各层的物理属性不发生变化，那么对应的广义反透射系数矩阵不会发生变化． 对于有限震源情形，通常只有震源层界面发生变化，上行波和下行波的广义反透射系数矩阵的变化仅发生在震源层与其邻层之间． 所以，对于多个震源界面的情况，其它层之间的广义反透射系数矩阵无需重复计算．

由于各反透射系数矩阵的计算都类同，所以，我们借助于图 4和图 5仅阐述 R^SL_D的计算过程． 为了阐述问题的方便，如图 4所示，我们假设仅有2个不同深度的震源，震源界面分别为z_S1和z_S2，两者之间可以存在多个介质层．

图 4 不同深度的两个震源与介质底界面之间下行波的 广义反射路径及其系数矩阵 RS1LD和RS2LD的几何解释 zS1与zS2分别表示不同深度的两个震源界面，RS1LD(绿线)与 RS2LD(红线)分别表示 两个震源界面与介质底界面之间下行波的广义反射系数矩阵． zL表示介质底界面， zS1与zL之间具有M层介质，zS2与zL之间具有N层介质，且M>N Fig. 4 Geometrical illustration of generalized reflection traces and their matrices RS1LD and RS 2 L D of down-going wave between the bottom interface and two source-interfaces zS1 and zS2 denote interfaces of the two depth-different sources， and RS1LD and RS2LD denote generalized reflection matrices of the down-going waves between bottom interface and two source-interfaces，respectively. zL denotes bottom interface. There are M layers between zS1 and zL， and N layers between zS2 and zL， and M>N

图 5 两个震源情况下广义反射系数矩阵 RS1DL和S2DL的计算流程示意图 (a)表示已有程序的计算过程；(b)表示本文程序的计算过程． 深度为zS1的震源(水平绿线)与底界面zL之间存在M层介质，绿色箭头表示 RS1LD的计算过程； 深度为zS2的震源(水平红线)与底界面zL之间存在N层介质，红色箭头表示 RS2LD的计算过程 Fig. 5 Computing flow of the generalized reflection matrices RS1LD and RS2LD (a)Computing flow of the referenced program;(b)Computing flow of GRTMatSyn. Horizontal green line is the interface of zS1 with M interfaces above the bottom interface， and green line with arrow shows the computing flow of RS1LD; Horizontal red line is the interface of zS2 with N interfaces above the bottom interface， and red line with arrow shows the computing flow of RS2LD

根据上面提到的广义反透射系数矩阵的计算方法，R^SL_D依赖于此震源分界面与介质底界面之间的介质结构． 具体计算需要自介质底界面z=z_L开始，逐层计算，直至z=z_S的震源界面(Kennett，1983； 李旭，1993； 李旭，陈运泰，1996)． 如图 5所示，假设每个界面的计算耗时为1个时间单位，那么，单独计算 R^S₁L_D需要M个时间单位(绿线)，R^S₂L_D需要N个时间单位(红线)，分别计算 R^S₁L_D与 R^S₂L_D则需要(M+N)个时间单位． 但是，分析两者的迭代过程可以发现，计算 R^S₁L_D或者 R^S₂L_D的差别仅发生在(N+1)层与M层之间，前面的过程完全相同． 因此，在同一循环过程中，可以同时计算 R^S₁L_D与 R^S₂L_D，这样可以节省N个时间单位． 不难想象，随着不同深度的震源数目的增加，节省的时间会更多．

式(4)和式(5)建立了地震波矢量间断与描述一般点源的地震矩张量之间的联系． 我们知道，任意矩张量解或者任何一种震源机制解都可以表示为6个基本矩张量的加权和，即

因此，我们可以预先计算6个基本矩张量对应的地震波矢量间断(∑_U(z_S)和∑_D(z_S))，任意机制的震源引起的地震波矢量间断都可以由此基本的地震波矢量间断合成． 很容易想到，如果有限震源的震源机制少于6种，采取这样的技术路线是不划算的； 如果多于6种，则是节省时间的． 然而实际情况是，一次破坏性地震的有限震源往往由远大于6的点源构成． 所以，在一般情况下，采用上面的技术路线是经济实用的．

从上面讨论可以看出，6个基本矩张量对应6组地震波矢量间断(∑_Uj(z_S)和∑_Dj(z_S)，j=1，6)， 每个震源深度对应一套广义反射透射系数矩阵 R^0S_D，T^0S_U， R^SL_D和R^0S_U． 将其代入式(1)，即可得到6种基本震源引起的面谐矢量位移解w_0j(j=1，6)，则任意矩张量引起的w₀可以表示为

有限的震源可以看作多个点源的组合，每个点源有它自己的深度和震源机制． 通过上述技术途径可以很快得到每个点源的面谐矢量位移解，然后将这些面谐矢量位移解经Fourier-Hankel反变换得到时间-空间域的位移，即地震位移，最后将所有点源在观测点的位移叠加即可得到有限震源引起的地面运动． 概而言之，计算有限震源的地面运动的计算流程如下：

1)计算层状介质中各界面的反射透射系数矩阵、 自由表面的反射系数矩阵和接收系数矩阵.

2)计算6个基本矩张量引起的地震波矢量间断.

3)计算不同震源深度点源对应的广义反射透射系数矩阵.

4)将地震波矢量间断、 广义反射透射系数矩阵、 自由表面的反射系数矩阵和接收系数矩阵组合，得到不同震源深度点源的6个基本矩张量引起的面谐矢量位移解.

5)根据给定的有限震源的点源的深度和震源机制，调用4)的计算结果，生成此点源引起的面谐矢量位移解，并对其进行Fourier-Hankel反变换，得到时-空域的介质响应或格林函数.

6)将5)生成的介质响应与给定的震源时间函数褶积，生成考虑时间过程的点源引起的地面运动.

7)叠加所有点源在给定观测点引起的地面运动，生成有限震源引起的地面运动．

为进一步提高程序的计算效率，我们对计算过程中的两个环节采用了并行计算设计． 第一个环节是广义反透射系数矩阵 R^0S_D，T^0S_U，R^SL_D和R^0S_U的计算． 该环节包括自由表面与震源之间的广义反透射系数矩阵 R^0S_D，T^0S_U和R^0S_U的计算及震源与介质底界面之间广义反射系数矩阵 R^SL_D 的计算两部分. 二者并不发生联系，因此，采用了并行计算． 第二个环节是有限震源的地面运动计算． 在这个环节上有两种可行的并行方案: 第一种是把有限的震源分成多组而观测点为一组； 第二种是把观测点分为多组而震源为一组． 考虑到在通常情况下观测点数目远多于震源的数目，我们采用了后者，即并行计算同一有限震源在不同观测点的地面运动．

为了测试程序的可靠性并检验其计算效率，我们选择相同的介质模型和震源模型，分别利用GRTMatSyn和被广泛使用的程序CompSyn(Spudich，Xu，2003)计算给定点的地面运动，对二者的计算结果和计算时间进行对比分析． CompSyn采用Olson等(1984)的离散波数有限元方法计算格林函数，并采用Spudich和Archuleta(1987)提出的数字技术实现断层面上的位移积分． CompSyn的特点是计算的格林函数反映了介质的完全响应(包括P波、 S波、 面波以及近场项)，并适用于有限的平面断层． 然而，由于采用了有限元算法，所以计算量随着观测点离震源的距离增加而增加； 由于时代的局限性，并没有考虑并行计算； 另外，由于局限于平面断层，所以组成有限断层的点源的走向和倾角不能变化． 我们 之所以使用CompSyn与GRTMatSyn进行比较，是因为二者计算的格林函数都包含了介质的完全响应，且都可以计算有限平面断层引起的地面运动．

由于CompSyn采用震源坐标，所以我们不防定义Z轴垂直向下为正，N轴向北为正，E轴向东为正，坐标系遵循右手法则． 测试断层如图 6所示，走向为0°，倾角为90°，滑动角为0°； 断层埋深5 km； 断层走向方向长6 km，倾向方向宽3 km； 滑动量均为0.75 m，震源时间函数为持续时间1 s的箱型函数． 介质模型选用两层半介质模型，具体参数如表 1所示． 另外，采用的频率范围为0—0.5 Hz，慢度区间0.001—1 s/km，采样间隔为0.25 s，波形时间长度为128 s． 我们计算了如图 6和表 2所示的6个观测点的地面运动．

表 1(Table 1)

表 1 介质模型参数 Table 1 Parameters of the stratified media 注：vP表示纵波波速，vS表示横波波速，ρ表示介质层密度，h表示介质层厚度，Qα表示纵波品质因子，Qβ 表示横波品质因子．

表 1 介质模型参数 Table 1 Parameters of the stratified media

图 6 有限断层模型与观测台站分布 Fig. 6 The model of finite source and observation stations

表 2(Table 2)

表 2 6个观测台站的N-E坐标 Table 2 Coordinates of observation stations

表 2 6个观测台站的N-E坐标 Table 2 Coordinates of observation stations

为了同时测试程序的可靠性和计算效率，如图 7所示，我们设置了4种情形，分别利用1个点源、 3个点源、 6个点源和9个点源来表征有限断层模型． 为了定量比较两组结果之间的一致性或者差异，计算和评价了对应波形的相关系数和最大振幅差． 最大振幅差由计算得到，这里 P₀表示 CompSyn计算结果的最大振幅，P表示GRTMatSyn计算结果的最大振幅． 另外，为了节省篇幅，在此我们只展示南北分向的波形比较． 更多的测试参见邸海滨(2011)文章．

图 7 利用不同数量点源表征的有限震源模型 (a)1个点源；(b)3个点源；(c)6个点源；(d)9个点源 Fig. 7 Four cases of finite source approximated with point-sources (a)1 point-source;(b)3 point-sources;(c)6 point-sources;(d)9 point-sources

图 8展示了4种情况下GRTMatSyn与CompSyn计算结果的比较． 在每个子图中，粗线表示CompSyn计算的地面运动，细线表示GRTMatSyn计算的地面运动． 子图左边的数字从上而下依次为： CompSyn计算结果的最大振幅、 两个结果之间的相关系数和GRTMatSyn计算结果的最大振幅； 右边的字符自上而下依次为观测台站编号与南北向分量； 最下方给出波形对应的时间刻度． 从图 8可以看出，仅用一个点源代替有限断层时，距震源 较近的观测点的波形的一致性较差，相关系数较小，最大振幅差较大. 例如，最近的观测点(No.1)的相关系数仅为0.72，最大振幅差达40%． 随着震中距的增大，二者之间的相关系数增大，最大振幅差减小. 例如，最远的观测点(No.6)的相关系数达到0.92，振幅差减小至24.6%． 不过，随着点源数目的增加(从1到9)，两组结果之间的差异减小，即使在最近的观测点，两个波形之间的相关系数也可以达到0.96，最大振幅差不超过10%． 我们把点源数目增加到更多(例如15)，但相关系数和最大振幅差也基本不变． 我们认为，两组结果之间的这种差异很可能与采用的数值算法有关. 比如，在计算贝塞尔函数时，GRTMatSyn采用精确计算的方法，而CompSyn采用查表并插值的方法．

图 8 4种情况下CompSyn(粗线)与本文程序(细线)计算南北向波形的比较 (a)采用1个点源表示有限震源模型；(b)采用3个点源表示有限震源模型； (c)采用6个点源表示有限震源模型；(d)采用9个点源表示有限震源模型 Fig. 8 Comparison of the waveforms computed with CompSyn(bold line) and GRTMatSyn(light line)in four cases (a)With 1 point-source;(b)With 3 point-sources; (c)With 6 point-sources;(d)With 9 point-sources

为了解程序的计算效率，我们统计了相同情况下，GRTMatSyn与CompSyn的计算时间． 我们使用的计算机主要参数为四核CPU、 4G内存、 500G硬盘、 512M独立显卡． 表 3仅给出了3个点源情况下两个程序的计算耗时． 不难看出，计算同一观测点的地面运动，GRTMatSyn所用时间较少，而CompSyn所用的时间较多； 计算不同的观测点的地面运动，GRTMatSyn所用时间不随震中距的变化而变化，而CompSyn所用的时间随着震中距的增大而增大．

表 3(Table 3)

表 3 本文程序(GRTMatSyn)与CompSyn计算耗时的比较 Table 3 Comparison of the computing time between our program and CompSyn 观测台站 No.1 No.2 No.3 No.4 No.5 No.6 本文程序 12.6 s 12.6 s 12.6 s 12.6 s 12.6 s 12.6s CompSyn 90 s 93 s 101 s 109 s 116 s 131 s

表 3 本文程序(GRTMatSyn)与CompSyn计算耗时的比较 Table 3 Comparison of the computing time between our program and CompSyn

在程序的适用性方面，GRTMatSyn的设计考虑了任意几何特征的有限震源情况以及有限震源的震源机制随时间和空间发生变化的情况． 考虑到大地震断层可能不是一个平面，所以，GRTMatSyn的设计没有局限于平面断层； 考虑到震源机制在震源过程中可能会发生变化，所以，GRTMatSyn从底层设计时就考虑了任意震源机制的需求． 即便是一组和多组爆炸源或其它类型的震源引起的地面运动，GRTMatSyn也能胜任．

在程序的计算效率方面，GRTMatSyn的设计考虑了计算流程的优化和并行计算． 计算流程的优化主要出现在通过一次循环计算所有不同深度的震源对应的广义反透射系数矩阵的环节上，避开早期的针对不同震源深度的多次循环． 并行计算主要出现在两个环节： 第一个环节是广义反透射系数矩阵的计算． 广义反透射系数矩阵的计算分上行波和下行波两种，二者并不发生联系，因此，设计了并行计算； 第二个环节是有限震源的地面运动计算． 在这个环节上，把观测点分为多组进行并行计算．

由于对适用性和计算效率的考虑，GRTMatSyn不但可以计算任意有限震源引起的地面运动，同样也能够计算点源引起的地面运动． GRTMatSyn的计算结果与CompSyn的计算结果的对比分析表明，GRTMatSyn的计算结果是可靠的，计算效率也得到了提高．

除计算点源和有限震源引起的地面运动外，GRTMatSyn在基础研究方面也将发挥积极的作用． 例如，介质结构的波形反演问题． 我们知道，介质结构的波形反演问题是一个非线性反演问题，因此，反演效率直接取决于波形的计算效率． 如果波形计算的效率高，介质结构的波形反演就可行且效率高； 反之，即不可行或效率不高． 又如，区域地震或地方地震的震源机制反演问题，不同于远场地震震源机制的波形反演. 区域地震或地方地震的震源机制反演，由于介质响应的复杂性和震源深度对波形的影响，也通常采用非线性反演，需要反复计算不同震源深度、 基本震源机制的震源产生的波形． 因此，反演效率也直接取决于波形的计算效率．

当然，该程序也有其局限性，仅适用于计算分层均匀介质中震源引起的地面运动． 在计算机技术高速发展的今天，三维介质中地震波的计算已成为可能，或将很快应用于实际地震的应急响应．

自由表面的反射系数矩阵和接收系数矩阵表达式(Kennett，1983)为

对于P-SV波

对于SH波，反射系数矩阵=1，接收系数矩阵=2．

已知地震波在传播过程中，遇到地层分界面会产生反射与透射作用，反射与透射系数矩阵被用于描述这一作用． 假设在深度z处存在一地层分界面，上下层面分别以z^-和z⁺表示． 根据广义反射透射系数矩阵方法的基本理论(Kennett，1983)，以 R D(z^-，z⁺)和T D(z^-，z⁺)表示下行波的反射系数和透射系数矩阵，R U(z^-，z⁺)和T U(z^-，z⁺)表示上行波的反射系数和透射系数矩阵，则

已知地层分界面处的传递矩阵 Q (z^-，z⁺)具有如下表达式(Kennett，1983):

对于P-SV波，

对于SH波，
式中，i 表示虚数单位，α与β分别表示P波与S波的波速，p表示慢度，ρ表示介质密度，q_α=(α^-2-p²)^1/2，q_β=(β^-2-p²)^1/2，且Im ωq_α≥0，Im ωq_β≥0，ε_α=(2ρq_α)^-1/2，ε_β=(2ρq_β)^-1/2.

将公式(B-3)代入(B-1)和(B-2)，即可完成地层分界面处反射透射系数矩阵的计算．