(1)
北京大学学报(自然科学版) 第61卷 第6期 2025年11月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 61, No. 6 (Nov. 2025)
doi: 10.13209/j.0479-8023.2025.084
摘要 针对精确展平变换方法计算地球的振荡模式时, 在搜根阶段容易出现漏根的问题, 在广义反射–透射系数法(GRTCM)框架中引入精确展平变换, 将广义反射–透射系数法推广到环状模式球状无重力弹性地球模型, 然后利用久期函数族计算该模型中的环形振荡模式, 计算结果与 MINEOS 的精确结果非常吻合。
关键词 环形振荡模式; 久期函数族; 精确展平变换; 广义反射–透射系数法(GRTCM)
地球内部的物理特性和结构(包括密度、弹性和不连续性)决定其简正振型的特征频率和模式, 因此地球的自由振荡可以用来反演地球结构的横向变化, 并表征激发它们的源[1–2]。同时, 简正振型的数值计算也是简正振型叠加法合成理论地震图的重要环节[3]。因此, 在理论地震学的研究中, 简正振型的计算一直是非常重要的课题。
在理论地震学研究中, 弹性动力学的偏微分方程在引入正交基函数后呈现为常微分方程组的形式, 特征函数与特征值由这个常微分方程组控制。对这个常微分方程组的求解, 传统的方法主要有直接数值积分或矩阵法两种, 这些半解析方法都在高频下表现出数值不稳定性。Gilbert 等[4]曾提出“次级向量(minor vector)”方法, 通过搜寻解集的子行列式来解决这个问题。Woodhouse[5]详细地总结了计算简正振型的数值方法, 使得简正振型在数值上的有效计算成为可能。基于文献[5]中的算法, Masters等[3]编写了程序 MINEOS, 实现简正振型的数值计算, 克服了简正振型叠加法计算理论地震图的主要困难。但是, 在实际应用中, 该方法仍然可能面临数值精度下降的问题, 特别是当频率超过某一阈值时。MINEOS 在低频时表现良好, 但在周期为数秒时, 不能很好地工作。在实际使用中, MINEOS 的数值精度会在大约 50mHz 之后下降。因此, 提高高频球状振型计算精度的问题仍然有待解决。Zhao 等[6–7]用 Wentzell-Kramers-Brillouin-Jeffreys(WKBJ)近似法, 在具有有限不连续点的地球模型中, 基本上解决了 MINEOS 的上述问题。但是, 因为 WKBJ方法在转折点附近存在固有缺陷, 导致与 MINEOS的结果相比, 出现模式缺失和明显的数值差异。
与球状地球模型不同, 在平层的地球模型中, 面波的上述数值不稳定性问题通过 Chen[8]提出的广义反射–透射系数法(generalized reflection-transmis-sion coefficient method, GRTCM)简洁而高效地解决了。该方法基于 Luco 等[9]和 Apsel 等[10]的技术, 通过引入重整化和虚频率[11], 大幅度地提升数值稳定性。同时, 这套方法非常适合迁移到更复杂的模型中处理地球物理问题, 例如拓展到流体饱和层状孔隙介质中震电波场[12]的计算。此外, 该方法还被不断地优化。Pei 等[13]指出, 广义反射–透射系数可以直接计算, 进一步提升该算法的计算效率。袁腊梅 等[14]在广义反射–透射系数的算法基础上, 对解的形式稍做变形, 使得本征模求解涉及的矩阵无量纲化, 矩阵各元素在量级上不至于差别太大, 提升了计算精度。
在面波频散曲线的计算中, 广义反射–透射系数法发展出类似 Kennett[15]提出的高效、稳定计算面波问题的方法, 并得到广泛的应用。但是, 在模型含低速层的情况下, 该方法需要将搜根步长分割得很细, 从而导致计算速度很慢, 不然就可能漏掉能量主要陷在低速层中的那些振型。为了解决这一问题, 何耀锋[16]引入久期函数族的概念, 发现不同的久期函数都能快速地搜到一部分根, 启用整个久期函数族就能保证不漏根。Chen 等[17]在后续研究中应用久期函数族, 成功地将海洋地球模型产生的Stoneley modes 计算到很高的频率并保证精度。Fan等[18]根据高频情况下频散方程的近似分解, 对频散曲线上的振型进行分类, 并导出给定频率下振型数目的预测公式和实际的搜根策略, 对频散方程的加速求解有重要指导意义。对于球状介质, Wu 等[19]采用 Langer 近似法[20], 为不均匀球形层提供统一的渐近解, 使用广义反射–透射系数法, 有效地计算地球的环状振型。
展平变换是一种将球坐标系下的地震波方程转换到直角坐标系下, 然后基于水平地球模型的成熟理论来计算球状地球模型的波动问题。对于 SH 波, Biswas 等[21]给出精确的展平变换公式。对于 P-SV波, Biswas[22]忽略与地球半径成反比的高阶小项, 提出 P-SV 波展平变换的近似求解公式。Bhattachar-ya[23]在球状分层介质模型中, 假设每一层介质的波速和拉梅系数随球体半径 r 变化, 提出适用于 P-SV 波的精确展平变换公式。Bhattacharya[24]利用精确的展平变换公式, 简单地修改 Wang[25]的程序中的传输矩阵, 得到理论地震图。精确的展平变换还被 Bhattacharya[26]引入面波速度的计算中, 并证明行之有效。但是, 由于没有使用鲁棒性强的数值方法, 只简单地使用 Chen[8]1993 年提出的久期函数进行计算, 搜根效率很低, 并且容易存在漏掉基阶根的缺陷。
受 Wu 等[19]采用 Langer 近似来求解环状振型的思路启发, 我们合理地推测, 采用精确展平变换, 可以将球状地球模型中求解简正振型的问题转化成广义反射–透射系数法体系可以解决的问题, 从而有效地计算地球的简正振型。本研究依照精确的展平变换基本思路, 将球状地球模型下的弹性动力学方程变换为可以使用广义反射–透射系数方法求解的形式, 利用成熟的广义反射–透射系数方法, 为球状地球模型中面波频散的计算提供一条新的路径。
考虑分层的无重力横向各向同性介质中的无源的弹性动力学方程[27]:
(1)
式中, l 和m 为拉梅系数, r 为密度, u 为时间域中质点的位移量。我们在图 1 所示球状分层模型中求解该方程。引入球坐标系(r, q, f), 假设任意一层中的介质参数l,m 和r 都是关于 r 的函数, 在频率域中由式(1)可以推导出
(2)
其中, U 为频率域中指点的位移量。假设第 j 层中介质参数(图 1)满足如下关系:
图1 p=2时球状分层介质模型的示意图
Fig. 1 Schematic diagram of a spherically stratified medium model with p=2
(3)
其中, l0, m0, r0 和 p 为任意常数, a 为地球半径。我们利用矢量基函数
,
和
[28]展开位移场 U:
(4)
将式(4)代入式(2), 通过解耦得到环状振型的控制方程:
(5)
从式(5)可以看出, p 取值为−2 时, 式(4)会得到极大的简化。
本文主要研究环状振型, 故定义位移应力矢量
进行推导。式(5)中, 关于 Wlm 的方程为特殊结构的变系数线性常微分方程, 称为欧拉方程。引入变换
(6)
式(6)在地球物理学中称为精确展平变换[29]。值得注意的是, 由于 r≤a, 所以 z 为负值, 与平层中的坐标轴 z[8]不同, 此处 z 向上表示增加。将展平变换(式(6))和变换后的位移应力矢量
代入关于 Wlm 的方程(式(5))中, 在任意层j中, 其介质参数为lj, mj 和rj, 可得一阶线性常微分方程组
(7)
(8)
(9)
通过上述变换, 
中不再含自变量 r 或 z, 我们达成了前面所述目标: 将球状地球模型下的环状振型的弹性动力学方程转化为关于
的一阶线性常微分方程组, 该方程组适合用广义反射–透射系数方法进行处理。
仅考虑地表到核–幔边界(core mantle boundary, CMB)的介质模型(精确展平变换在 r=0 存在奇异性), 我们使用广义反射–透射系数方法[11]对式(7)进行处理。方程(7)的通解为
(10)
(11)
(12)
其中, 
。根据第 j 层与第 j+1 层之间的界面连续性条件, 可得
(13)
为了推导的方便, 我们引入
和
来代表下行波与上行波相关的项。值得一提的是, 展平变换后
为负值, 因此式(12)的定义与常规定义不同。
参照广义反射–透射系数[30]的推导, 我们定义下行波相关的广义反射–透射系数:
(14)
定义上行波相关的广义反射–透射系数:
(15)
分别代入式(13)中, 可得递推关系式:
(16)
(17)
其中, Fij 是
的元素, Gij 是
的元素。
考虑到 jcmb 层为核–幔边界的固体层, 我们除设定自由表面边界条件外, 还设定核–幔边界牵引力为零的边界条件。由两个边界条件可以推得
(18)
久期函数族的概念最早由何耀锋等[31]于 2006年提出, 是伴随广义反射–透射系数法发展出来的一套用于研究面波频散的理论。在广义反射–透射系数法中, 我们求得广义反射–透射系数后, 就可以定义频散方程来求取给定频率处的根(即相速度)。首先考虑固体层中的环状振型, 在图 1 所示的第 j层中, 有如下关系:
(19)
根据式(18), 可以得到
(20)
根据广义反射–透射系数法求解面波频散曲线的研究结果[8], 式(19)具有非零解就意味着系数行列式为零, 即
(21)
所以, 介质层数与频散方程数相等。另外, 根据自由表面应力为零的边界条件, 可以得到第一层的频散方程:
(22)
在文献[31]中, 式(20)的左边为久期函数, 第 j层介质中的久期函数记为 Sj; 式(21)是与自由表面相关的频散方程, 其左边记为久期函数 S0。于是, 就有 N+1 个久期函数, 统称为久期函数族:
(23)
从理论上讲, 由久期函数族中的久期函数构成的任意一个频散方程都包含给定频率对应的所有本征模的信息。然而, 实际上每一个久期函数搜根都有各自的优势。一般来说, 每一层的久期函数都比较容易搜到大于该层 S 波速度的相速度, 对于相速度小于该层 S 速度的振型, 则较难搜到。因此, 如果利用所有低速层对应的久期函数, 加上直接根据自由表面边界条件导出的久期函数, 就能够高效地搜出任意频率下所有实际存在的本征模。
在我们的理论推导中, 每一层的介质变化如图1 所示, 用这种特殊的介质模型来近似地球的介质变化, 需要分非常多的层(即 N 取较大值)来最小化层内的不均匀性。我们以无海洋的 PREM 模型为例, 对该模型中地表到核–幔边界之间的介质模型进行拟合, 结果如图 2 所示。
为了最小化层内不均匀性的影响, 每层保持较小的厚度。假设给定模型中第 j层的 P 波和 S 波速度以及密度的平均值分别为
, 介质参数在第 j 层中的变化分别如式(3)所示, 由此可得 p=2时的波速与密度:
(24)
其中, 常数a0j, b0j 和r0j 的选择方式是使得层内 𝛼(𝑟), 𝛽(𝑟)和𝜌(𝑟)的平均值分别等于。
图2 无海洋的PREM模型
Fig. 2 PREM model without oceans
图3 GRTM对PREM模型计算的环形振荡的本征频率(红色加号)与MINEOS的精确计算结果(蓝色圆点)比较
Fig. 3 Comparison of the eigenfrequencies of toroidal oscillations calculated for the PREM model (indicated by red plus signs) with the precise calculation results from MINEOS (indicated by blue dots)
图4 广义反射–透射系数法得出的特征频率与MINEOS结果之间的相对误差
Fig. 4 Relative errors of the GRTM eigenfrequencies with respect to the MINEOS results for the PREM model
我们将久期函数组方法应用于图 2 所示的地球模型中, 计算环状振型, 并与 MINEOS 的计算结果进行比对(图 3)。本文方法计算出的环形特征频率与 MINEOS 计算的精确结果极为吻合, 并且没有漏根现象。由于该方法对地球模型的连续性没有要求, 故不论实际模型存在多少个间断面, 都能作为程序的输入。
通过计算广义反射–透射系数法中利用久期函数族得出的特征频率与 MINEOS计算结果之间的相对误差来评估算法的准确度。相对误差定义为
(25)
计算结果如图 4 所示。对于图 4 中环形振荡模式的特征频率, 久期函数族的计算结果具有很高的准确性, 相对误差均低于 0.5%, 其中 69.2%的相对误差低于 0.3%。相对误差未呈现规律性, 根据对算法的分析, 我们认为这是使用幂函数对地球介质模型进行拟合导致的误差。地球的介质参数l, m和r 在 核–幔边界以上随着 r 增大而不断减小。但是, 精确展平变换需要的特殊介质参数(式(23))是幂函数, α(r)和 β(r)在每一层中随着 r 增大而不断增大, 因此在拟合的地球介质中存在非常多的振荡现象。针对这一问题, 可以通过增加分层模型的层数来减小误差。
本文利用精确展平变换的基本思路, 通过设定球状层中的特殊介质参数, 将球状模型中变系数的常微分方程组转化为数学意义上的欧拉方程组, 再通过变换展平变换, 将方程组变成可以使用广义反射–透射系数法求解的常微分方程组。在此基础上, 通过引入久期函数族这一有效的工具, 对球状无重力弹性地球模型中的环形振荡模式进行计算。这种新的计算方法解决了展平变换求取频散信息时容易出现的漏根问题, 提升了搜根效率。并且, 相对误差的计算结果证明, 这种方法具有较高的精确性。基于久期函数族的搜根策略, 也让程序在求解阶段有非常高的计算效率。相对于 MINEOS, 该方法在高频时仍然具有良好的计算效率。
在未来的工作中, 如果在考虑引入重力项的前提下使用广义反射–透射系数方法计算环状振型, 需要重新建立弹性动力学方程。重力加速度满足如下微分关系:
(26)
可以看到, 式(26)仍然适用于欧拉方程理论, 但需要重新定义位移应力矢量来完成考虑重力的球状弹性地球模型中环形振荡模式的计算。
致谢 研究工作中得到南方科技大学李正波博士和鹏城实验室吴勃博士的帮助, 在此表示衷心的感谢。
参考文献
[1] Gilbert F, Dziewonski A M, Bullard E C. An appli-cation of normal mode theory to the retrieval of struc-tural parameters and source mechanisms from seismic spectra. Philosophical Transactions of the Royal Socie-ty of London Series A: Mathematical and Physical Sci-ences, 1975, 278: 187–269
[2] Dziewonski A M. Studies of the seismic source using normal mode theory. Proc Enrico Fermi Sch Phys, 1983, 85: 45–137
[3] Masters G, Woodhouse J H, Freeman G. Mineos v1.0.2, computational infrastructure for geodynamics [EB/ OL]. (2011) [2024–02–17]. https://geodynamics.org/ resources/mineos
[4] Gilbert F, Backus G E. Propagator matrices in elastic wave and vibration problems. Geophysics, 1966, 31 (2): 326–332
[5] Woodhouse J H. The calculation of the eigenfrequ-encies and eigenfunctions of the free oscillations of the earth and sun. Seismological Algorithms: Computatio-nal Methods and Computer Programs, 1988, 230: 321–370
[6] Zhao L, Dahlen F A. A symptotic normal modes of the earth — II. eigenfunctions. Geophysical Journal Inter-national, 1995, 121(2): 585–626
[7] Zhao L, Dahlen F A. A symptotic eigenfrequencies of the Earth’s normal modes. Geophysical Journal Inter-national, 1993, 115(3): 729–758
[8] Chen X. A systematic and efficient method of compu-ting normal modes for multilayered half-space. Geo-physical Journal International, 1993, 115(2): 391–409
[9] Luco J E, Apsel R J. On the Green’s functions for a layered half-space, part I. Bulletin of the Seismological Society of America, 1983, 73(4): 909–929
[10] 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): 930–951
[11] Chen X. Seismogram synthesis in multi-layered half-space, part I. theoretical formulations. Earthquake Re-search in China, 1999, 13(2): 149–174
[12] Ren H, Chen X, Huang Q. Numerical simulation of co-seismic electromagnetic fields associated with seismic waves due to finite faulting in porous media. Geophy-sical Journal International, 2012, 188(3): 925–944
[13] Pei D, Louie J N, Pullammanappallil S K. Improve-ments on computation of phase velocities of Rayleigh waves based on the generalized R/T coefficient me-thod. Bulletin of the Seismological Society of America, 2008, 98(1): 280–287
[14] 袁腊梅, 凡友华, 孙书荣. 广义反射–透射系数算法的无量纲化. 地震学报, 2009, 31(4): 377–384
[15] Kennett B. Seismic wave propagation in stratified me-dia. Cambridge: Cambridge University Press, 2009
[16] 何耀锋. 面波频散反演问题的初步研究[D]. 北京: 北京大学, 2005
[17] Chen W, Chen X. Modal solutions in stratified multi-layered fluid-solid half-space. Science in China Series D: Earth Sciences, 2002, 45(4): 358–365
[18] Fan Y H, Chen X F, Liu X F, et al. Approximate decomposition of the dispersion equation at high fre-quencies and the number of multimodes for Rayleigh waves. Chinese Journal of Geophysics, 2007, 50(1): 222–230
[19] Wu B, Zhao L, Chen X. Uniformly asymptotic eigen-solutions of the Earth’s toroidal modes. Geophysical Journal International, 2022, 228(1): 250–258
[20] Woodhouse J. Asymptotic results for elastodynamic propagator matrices in plane-stratified and spherically-stratified earth models. Geophysical Journal Interna-tional, 1978, 54(2): 263–280
[21] Biswas N, Knopoff L. Exact earth-flattening calcu-lation for Love waves. Bulletin of the Seismological Society of America, 1970, 60(4): 1123–1137
[22] Biswas N. Earth-flattening procedure for the propaga-tion of Rayleigh wave. Pure Appl Geophys, 1972, 96: 61–74
[23] Bhattacharya S N. Earth-flattening transformation for P-SV waves. Bulletin of the Seismological Society of America, 1996, 86(6): 1979–1982
[24] Bhattacharya S N. Synthetic seismograms in a spheri-cal earth using exact flattening transformation. Geo-physical Research Letters, 2005, 32(21): doi: 10.1029/ 2005GL024152
[25] Wang R. A simple orthonormalization method for stable and efficient computation of Green’s functions. Bulletin of the Seismological Society of America, 1999, 89(3): 733–741
[26] Bhattacharya S N. Computation of surface-wave velo-cities in a nongravitating elastic spherical earth using an exact flattening transformation. Bulletin of the Seismological Society of America, 2009, 99(6): 3525–3528
[27] Shapiro N, Campillo M, Margerin L, et al. The energy partitioning and the diffusive character of the seismic coda. Bulletin of the Seismological Society of Ame-rica, 2000, 90(3): 655–665
[28] Dahlen F, Tromp J. Theoretical global seismology. Princeton: Princeton University Press, 1998
[29] Arora S, Bhattacharya S, Gogna M. Rayleigh wave dis-persion equation for a layered spherical earth with exponential function solutions in each shell. Pure Appl Geophys, 1996, 147(3): 515–536
[30] Wu B, Chen X. Stable, accurate and efficient com-putation of normal modes for horizontal stratified mo-dels. Geophysical Journal International, 2016, 206(2): 1281–1300
[31] 何耀锋, 陈蔚天, 陈晓非. 利用广义反射–透射系数方法求解含低速层水平层状介质模型中面波频散曲线问题. 地球物理学报, 2006, 49(4): 1074–1081
Computation of Toroidal Modes in a Spherical Non-Gravity Elastic Earth Model Using a Family of Secular Functions
Abstract When using the exact flattening transformation method to compute the oscillation modes, root-leakage issues can easily arise during the root-search phase. By introducing the exact flattening transformation within the framework of the generalized reflection-transmission coefficient method (GRTCM), this method can be extended to the spherical non-gravity elastic Earth model. Subsequently, using the duration function families, we calculate the toroidal modes in this model. The results of the proposed method align closely with the precise outcomes from MINEOS.
Key words toroidal modes; family of secular functions; exact earth flattening transformation; generalized reflection-transmission coefficient method (GRTCM)
国家自然科学基金(92255307)资助
收稿日期: 2024–09–27;
修回日期: 2024–10–01