(1)
北京大学学报(自然科学版) 第61卷 第6期 2025年11月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 61, No. 6 (Nov. 2025)
doi: 10.13209/j.0479-8023.2025.085
摘要 针对球状分层地球模型, 通过假设每层中的介质参数具有特定的变化规律并引入展平变换, 将球状模型中变系数的常微分方程组转化为可以使用广义反射–透射系数法(GRTCM)求解的常系数的常微分方程组, 并计算理论地震图。计算结果表明, 对于计算震中距在 100°以内的理论地震图, 广义反射–透射系数法是正确且高效的。随着球状模型中层厚减小, 该方法的误差得到明显改善。
关键词 合成理论地震图; 广义反射–透射系数法(GRTCM); 展平变换
水平层状均匀介质中地震波的传播是地震学中一个经典的课题。在处理近场问题时(D<3°), 由于震中距较小, 地震波频率较高, 可以忽略地球的曲率, 采用水平地球模型来近似真实的地球模型, 从而计算地震波的传播速度。在均匀的水平层状介质系统中, 可以通过传输矩阵法[1]合成完整的理论地震图。基于传输矩阵法计算理论地震图的研究结 果[2–3], Chen[4]引入离散波数积分法[5], 建立一套能系统地模拟层状均匀半空间中完整弹性波场的方法——广义反射–透射系数方法(generalized reflection-transmission coefficient method, GRTCM)。利用广义反射–透射系数, 可以从根本上排除数值计算中不稳定的指数增长项。在高频、厚层情况下, 该算法具有比其他方法更好的精确性和稳定性。经过一系列的完善和优化研究[6–10], 该方法已广泛应用于分层半空间介质中理论地震图的合成和面波频散曲线的计算。
对于远震问题, 通常需要在球状地球模型中使用球坐标来求解波动方程[11]。球坐标系下, 求解地震的波动方程比直角坐标系(或柱坐标系)下困难。在地震学研究中, 应用较广泛的方法是简正振型叠加法[12], MINEOS 程序[13]便是这一方法的代表。MINEOS 程序首先计算球对称地球模型中的本征频率和本征函数, 然后进行简正振型叠加, 得到低频(<0.15Hz)下的理论地震图。但是, 当频率稍高时, 简正振型叠加法的计算就会变得极为耗时。
本研究致力于将广义反射–透射系数方法推广到球状地球模型中, 提供一种效率更高的球状地球模型中理论地震图的计算方法。展平变换[14–17]可将球坐标系下地震波的方程转换到直角坐标系下, 然后用关于水平地球模型的成熟理论方法计算球状地球模型的波动问题。对于 SH 波, Biswas 等[14]给出精确的展平变换公式。对于 P-SV 波, Biswas[16]忽略与地球半径成反比的高阶小项, 提出 P-SV 波展平变换的近似求解公式。Arora 等[18]在球状分层介质模型中, 假设每一层介质的波速和拉梅系数随半径 r 变化, 提出适用于 P-SV 波的精确的展平变换公式。Bhattacharya[19]利用精确的展平变换公式, 修改 Wang[20]程序中的传输矩阵, 得到理论地震图。同时, 精确的展平变换还被 Bhattacharya[21]引入面波速度的计算中。但是, 上述工作基于矩阵方法, 未能解决远场情况下的数值积分问题以及数值搜索中漏根、耗时巨大等一系列问题, 未能广泛地应用到地震学研究中。本研究依照精确的展平变换的基本思路, 将球状地球模型下的弹性动力学方程变换为可以使用广义反射–透射系数方法求解的形式, 然后利用成熟的广义反射–透射系数方法, 为球状地球模型中合成理论地震图提供一条新的路径。
不均匀无重力各向同性介质中的弹性动力学方程[22]为
(1)
式中, l 和m 为拉梅系数; r 为密度; u为质点的位移量。我们在图 1 所示球状分层模型中求解该方程。引入球坐标系(r, q, f), 假设任意一层中的介质参数l, m 和r 都是关于 r 的函数, 在频率域中由式(1)可以推导出
(2)
其中, F 表示频率域中的点源。第j 层中, 利用式(3)所示矢量基函数
和
[11]展开位移场U:
(3)
(4)
其中, 
。假设每层中介质参数满足如下关系:
(5)
其中,l0,m0,r0 和 p 为任意常数, a 为地球半径。将式(4)和式(5)代入式(2), 得到如下常微分方程组:
(6)
通过解耦, 可以将式(6)分为环状振型和球状振型分别来求解。为方便理解, 我们首先对环状振型(即定义位移应力矢量
)进行推导。
式(6)中关于Wlm的方程为特殊结构的变系数线性常微分方程——欧拉方程。引入展平变换:
(7)
将式(7)与变换后的位移应力矢量
代入式(6)中关于 Wlm的方程中, 在任意层 j 中, 介质参数为l j,m j 和r j, 可得一阶线性常微分方程组:
图1 球状分层介质模型
Fig. 1 Spherical layered model
(8)
(9)
(10)
(11)
通过这样的数学处理, 
不再包含自变量r 或 z, 这就将球状地球模型下的弹性动力学方程转化成关于
和
的一阶线性常微分方程组, 该方程组适合用广义反射–透射系数方法进行处理。
我们仅考虑地表到核–幔边界(core mantle boun-dary, CMB)的介质模型(由式(7)可知展平变换在 r= 0 时存在奇异性), 使用广义反射–透射系数方法[2]对式(8)进行处理。仍然以环状振型为例, 式(8)的通解为
(12)
(13)
(14)
(15)
其中, 
。对第 j 层和第j+1 层之间的界面应用连续性条件, 可得
(16)
为方便推导, 我们引入
和
来代表下行波与上行波相关的项。值得一提的是, 经展平变换后z 为负值, 因此式(14)的定义与常规定义不同。
参照 Wu 等[9]的推导, 我们定义下行波相关的广义反射–透射系数:
(17)
定义上行波相关的广义反射–透射系数:
(18)
分别代入式(16)中, 可得递推关系式
(19)
(20)
其中, Fij 是
的元素, Gij 是 G=(Ej+1)–1Ej的元素。
考虑到第 jcmb 层为核–幔边界的固体层, 我们设定的边界条件除自由表面边界条件外, 还有核–幔边界牵引力为零的边界条件。由两个边界条件可以推得
(21)
对于球状振型, 直接给出矩阵 Ej:
(22)
(23)
将
中
替换为
, 即可得到
; 将
中
替换为
, 即可得到
(i=1, 2, 3, 4)。
我们通过特殊的介质结构(式(5))和展平变换(式(7)), 将球状地球模型下的弹性动力学方程变换为可以使用广义反射–透射系数方法求解的形式。因此, 广义反射–透射系数方法中源项
与
的推导[2]在本文方法中同样适用。
根据广义反射–透射系数方法, 我们需要得到ω-k 域中解的显式表达式来合成理论地震图, 因此需要推导展平变换后的震源项。
1.3.1 单力点源
单力点源在球坐标系下形式为
(24)
使用式(3)和(4), 对式(24)进行分解, 得到 fr, fs 和 ft的表达式:
(25)
(26)
根据广义反射–透射系数方法的公式[2], 通过式(25)和(26)求得广义反射–透射系数方法中的源项
和
。
1.3.2 双力偶点源
双力偶点源在球坐标系下的形式为
(27)
使用式(3)和(4), 对式(27)进行分解, 得到 fr, fs和 ft的表达式:
(28)
根据广义反射–透射系数方法的公式[2], 通过式(27)和(28)求得广义反射–透射系数方法中的源项和。
因为广义反射–透射系数法中的数值积分方法较为成熟, 于是我们将式(4)的求和转化为积分来合成理论地震图。
我们通过泊松求和, 将关于角度级数 l 的求和转换为波数 k 的积分。这一变换基于沃森变换(Wat-son transformation), 计算公式如下:
(29)
式(29)是一个恒等式, 对在 k 的实轴附近解析的任何函数 f (k)都有效。积分沿着图 2 所示复平面中的封闭回路 C 进行, 被积函数在 k=1/2, 3/2, 5/2, …处有简单极点。
灰色圆点为 Watson 回路 C=C++C–; k=l+1/2 处简单极点的留数是 ip–1f(l+1/2)
图2 复波数平面的示意图
Fig. 2 Schematic diagram of the complex wavenumber plane
泊松求和公式是沃森公式的变体, 容易从方程(29)得到
(30)
球状地球模型中的格林函数张量如下:
(31)
代入泊松求和公式(式(30)), 则有
(32)
其中, 在点
处
。利用 Legen-dre 函数的性质, 可以得到
(33)
当 l → ∞时, 有
(34)
其中, Δ 表示震中距。因此, Legendre 函数的求和可以转化为贝塞尔函数的积分。从式(33)的结构可以知道, G 是多项积分的求和, 每一项对应着图 3 所示圆弧路径, 我们取长度最短的圆弧路径 1 来合成理论地震图。
图3 波在地球中传播时的路径编号
Fig. 3 Orbit numbering for waves traveling around the Earth
因此, 勒让德函数的求和可以转化为贝塞尔函数的积分。从式(33)的结构可以知道, G 是多项积分的求和, 每一项对应着图 3 所示的圆弧路径, 我们取长度最短的圆弧路径 1 来合成理论地震图。
本文的理论推导中, 每一层的介质变化如式(5)所示。用这种特殊的介质模型来近似地球的介质结构, 需要分较多的层(即 N 取较大值)来最小化层内的不均匀性。这里以无海洋的 PREM 模型为例, 对其地表到核–幔边界之间的介质模型进行拟合。
假设给定模型中第 j 层的 P 波和 S 波速度以及密度的平均值分别为
, 介质参数在第 j 层中的变化分别如式(5)所示, 则波速与密度在 p=2 时可以写成
(35)
选择适当的常数
,
和
, 使得层内𝛼(𝑟), 𝛽(𝑟)和𝜌(𝑟)的平均值分别等于。
用广义反射–透射系数方法合成全球模型下的理论地震图时, 在震中距较大的时候, 离散波数法的数值计算效率很低。这是因为用离散波数法计算波数域积分时, 积分的步长实际上是由时间窗长度和震中距共同决定的距离参数 L 确定的(离散波数的积分步长为 2π/L)。当震中距较大时, 相应的时间窗就会很长, 距离参数 L 会进一步变大, 积分步长变小, 数值计算效率就会变得很低。
针对这一数值计算中的困难, Chen 等[23]提出一种数值积分方法——自适应 Filon 积分方法, 有效地提高了远场情况下理论地震图的计算效率。本文引入这个方法来计算远场理论地震图。
震中距D=60°, 震源深度为 200km, 震源时间函数为用于模拟物理脉冲的平方半正弦波
图4 单力源情况下的合成理论地震图
Fig. 4 Synthetic theoretical seismograms for the case of a single force source
图5 单力点源的合成理论地震图的震相标注(以Z分量为例)
Fig. 5 Phase annotations of synthetic theoretical seismograms for a single force point source (using the Z-component as an example)
基于 PREM 模型, 在单力点源情况下计算震中距D=60°时的理论地震图。为模仿物理脉冲, 我们将震源时间函数设定为平方半正弦波, 即
sin2(pt/ t), 其中 τ 取 50ms。
图 4 为计算得到的理论地震图。将其中 Z 分量单独绘图, 得到图 5。使用 TauP 计算走时, 结果如表1 所示。合成理论地震图与 TauP 算法得到的走时之间的误差小于地震图时间域的采样间隔，因此两者的结果完全吻合。
QSSP 是一款使用简正振型理论计算球状地球模型中完整合成地震图的程序[24]。该程序采用混合算法, 对低频(小谐波次数)进行数值积分, 对高频(大谐波次数)采用解析的传播矩阵算法, 可以计算体波、面波、自由振荡、均匀海洋中的海啸、标准大气层中的次声波和静态变形。我们在震中距为60°时选取爆炸源, 使用 QSSP 合成理论地震图, 用0.005～0.1Hz 的带通滤波器进行滤波, 使用不同的最小层厚(dz)进行 PREM 模型的拟合, 结果如图 6中黑色曲线所示。当dz 为 20km 时, GRTCM 的计算结果的波形相对于 QSSP 表现为提前, 相关系数为0.864。这一问题在最小层厚为 2km 时得到明显的改善(相关系数提升到 0.972)。由此可见, 球状介质中 GRTCM 合成理论地震图的一个主要误差来源是对地球介质模型拟合的误差。
表1 单力点源情况下计算得到的理论地震图走时
Table 1 Travel time table in the case of a single force point source
算法走时/sP震相pP震相PcP震相S震相ScS震相 TauP602.57611.74648.011094.041191.08 地震图602.60611.60648.001094.101191.10
说明: 震中距为 60°, 理论地震图的震源深度为 200km。
仍然考虑 PREM 模型, 在双力偶情况下, 计算理论地震图的震源参数取为 strike = 150°, dip = 40°, rake = −27°, 震源深度为 200km。我们取多个震中距, 绘制理论地震图(图 7), 从其中能够看到面波及其频散现象。由于展平变换数学上的缺陷, 当震中距较大时, 无法观察到内核震相。由于我们使用泊松公式进行近似, 只计算了一条路径, 所以无法观察到从地球另一侧传播过来的地震波。
震源深度为 200km, 介质模型为 PREM, 均使用 0.005～0.1Hz 的带通滤波器进行滤波
图6 不同的分层厚度下使用GRTCM合成的理论地震图
Fig. 6 Theoretical seismograms synthesized using GRTCM with different layer thicknesses
震源参数: strike = 150°, dip = 40°, rake = −27°; 震源深度为200 km
图7 双力偶时多个震中距的理论地震图
Fig. 7 Theoretical seismograms at multiple epicentral distances for a double-couple source
表2 台站位于地表的理论地震图计算时间
Table 2 Calculation time of theoretical seismograms for stations located on the Earth’s surface
程序理论地震图计算时间/s采样频率为0.2 Hz采样频率为0.4 Hz采样频率为0.6 Hz采样频率为0.8 Hz采样频率为1.0 Hz GRTCM (本文)3661563193225803006 QSSP7322086254732063986
说明: 震中距为 60°, 震源深度为 30km。
表2 对比本文 GRTCM 程序与 QSSP 程序的计算效率。可以看出, 当震中距约为 60°时, 本文程序具有更好的计算效率。但是, 在接收点深度和震源深度未发生变化的情况下, QSSP 程序只需要计算一次格林函数, 就能合成多个震中距的理论地震图, 而广义反射–透射系数方法对不同的震中距需要进行不同的数值积分。
本文发展了一套适用于球状地球模型的合成理论地震图的半解析–半数值计算方法。对于球状分层地球模型, 用正交完备的矢量基函数展开波场。假设每层中的介质参数具有特定的变化规律, 从而将球状模型中的弹性动力学方程转化为可以使用广义反射–透射系数法求解的形式。对于广义反射–透射系数方法在球状分层地球模型中的数值问题, 本文也提出一系列解决方案。本文将这套方法的计算结果与 QSSP 等主流合成理论地震图程序的计算结果进行比对, 当 GRTCM 的震中距为 60°时, 在与QSSP 相同精度情况下具有更高的计算效率。本文的研究结果拓宽了广义反射–透射系数方法的适用范围, 不仅可以提高正演计算的效率和精度, 而且可以更深刻地理解地震波的物理性质。
本文计算结果中远震(100°以上)误差的震相缺失主要源于使用泊松求和公式展开后, 质点的位移会变成多个路径求和的形式, 我们使用第一项来近似完整的理论地震图, 导致远震震相的缺失。后续研究中, 我们将考虑增加对其他传播路径的积分来解决这个问题。
参考文献
[1] Thomson W T. Transmission of elastic waves through a stratified solid medium. Journal of Applied Physics, 1950, 21(2): 89–93
[2] Apsel R J, Luco J E. On the Green’s functions for a layered half-space, part II. Bulletin of the Seismolo-gical Society of America, 1983, 73(4): 931–951
[3] Yao Z X, Harkrider D. A generalized reflection-trans-mission coefficient matrix and discrete wavenumber method for synthetic seismograms. Bulletin of the Seis-mological Society of America, 1983, 73: 1685–1699
[4] Chen X F. Seismogram synthesis in multi-layered half-space, part I. theoretical formulations. Earthquake Re-search in China, 1999, 13(2): 149–174
[5] Bouchon M. Discrete wave number representation of elastic wave fields in three‐space dimensions. Journal of Geophysical Research: Solid Earth, 1979, 84(B7): 3609–3614
[6] Zhang H M, Chen X F. Self-adaptive Filon’s integra-tion method and its application to computing synthetic seismograms. Chinese Physics Letters, 2001, 18(3): 313–315
[7] 张海明, 陈晓非, 张似洪. 峰谷平均法及其在计算浅源合成地震图中的应用. 地球物理学报, 2001, 44(6): 805–813
[8] 何耀锋, 陈蔚天, 陈晓非. 利用广义反射–透射系 数方法求解含低速层水平层状介质模型中面波频散曲线问题. 地球物理学报, 2006, 49(4): 1074–1081
[9] Wu B, Chen X F. Stable, accurate and efficient compu-tation of normal modes for horizontal stratified mo-dels. Geophysical Journal International, 2016, 206(2): 1281–1300
[10] Wu B, Chen X F. Accurate computation of leaky modes for anomalous layered models. Annals of geophysics, 2017, 60(6): S0663
[11] Dahlen F, Tromp J. Theoretical global seismology. Prin-ceton: Princeton University Press, 1998
[12] Yang H Y, Zhao L, Hung S H. Synthetic seismograms by normal-mode summation: a new derivation and nu-merical examples. Geophysical Journal International, 2010, 183(3): 1613–1632
[13] Masters G, Woodhouse J H, Freeman G. Mineos v1.0.2, computational infrastructure for geodynamics [EB/OL]. (2011) [2024–02–17]. https://geodynamics.org/resour ces/mineos
[14] Biswas N, Knopoff L. Exact earth-flattening calcula-tion for Love waves. Bulletin of the Seismological So-ciety of America, 1970, 60(4): 1123–1137
[15] Müller G. Approximate treatment of elastic body waves in media with spherical symmetry. Geophysical Jour-nal International, 1971, 23(4): 435–449
[16] Biswas N. Earth-flattening procedure for the propaga-tion of Rayleigh wave. Pure Appl Geophys, 1972, 96: 61–74
[17] Bhattacharya S. Earth-flattening transformation for P-SV waves. Bulletin of the Seismological Society of America, 1996, 86(6): 1979–1982
[18] Arora S, Bhattacharya S, Gogna M. Rayleigh wave dis-persion equation for a layered spherical earth with ex-ponential function solutions in each shell. Pure Appl Geophys, 1996, 147(3): 515–536
[19] Bhattacharya S N. Synthetic seismograms in a sphe-rical earth using exact flattening transformation. Geo-physical Research Letters, 2005, 32(21): doi: 10.102 9/2005GL024152
[20] Wang R. A simple orthonormalization method for sta-ble and efficientcomputation of Green’s functions. Bull Seismol Soc Am, 1999, 89(3): 733–741
[21] Bhattacharya S N. Computation of surface-wave velo-cities in a nongravitating elastic spherical earth using an exact flattening transformation. Bulletin of the Seis-mological Society of America, 2009, 99(6): 3525–3528
[22] Singh S J, Ben‐Menahem A. Decoupling of the vector wave equation of elasticity for radially heterogeneous media. The Journal of the Acoustical Society of Ame-rica, 1969, 46(3B): 655–660
[23] Chen X F, Zhang H M. An efficient method for compu-ting Green’s functions for a layered half-space at large epicentral distances. Bulletin of the Seismological So-ciety of America, 2001, 91(4): 858–869
[24] Wang R, Heimann S, Zhang Y, et al. Complete synthe-tic seismograms based on a spherical self-gravitating earth model with an atmosphere-ocean-mantle-core structure. Geophysical Journal International, 2017, 210 (3): 1739–1764
Generalized Reflection-Transmission Coefficient Method for Synthesizing Theoretical Seismograms in a Spherical Earth Model
Abstract For a spherical layered Earth model, by assuming that the medium parameters in each layer follow specific variation laws and introducing earth flattening transformation, the ordinary differential equations (ODEs) with variable coefficients in the spherical model are transformed into ODEs with constant coefficients that can be solved using the generalized reflection-transmission coefficient method (GRTCM), and the theoretical seismograms are calculated. Numerical results show that the generalized reflection-transmission coefficient method is accurate and efficient for computing theoretical seismograms at epicentral distances within 100 degrees. The error of the method will be significantly reduced as the layer thickness of spherical model is reduced.
Key words Synthetic theoretical seismograms; Generalized reflection-transmission coefficient method (GRTCM); earth flattening transformation
国家自然科学基金(92255307)资助
收稿日期: 2024–09–27;
修回日期: 2024–10–06