• ISSN 2096-8957
  • CN 10-1702/P

远震波场正演模拟方法及应用

桑莹泉 刘有山 徐涛 白志明 解桐桐

引用本文:
Citation:

远震波场正演模拟方法及应用

    作者简介: 桑莹泉(1998-),男,硕士生,主要从事地震学研究. E-mail:sangyq@mail.iggcas.ac.cn.
  • 中图分类号: P315

Forward modeling method and application of teleseismic wavefield

  • CLC number: P315

  • 摘要: 地震波场数值模拟是壳幔结构成像和深部探测的重要基础. 经典的远震波场模拟主要基于一维地球模型,包括解析法、半解析法和数值法等. 这些算法能够高效地计算理论地震图,但由于将地球假设为一维层状介质,难以考虑介质的横向非均匀性对地震波场的影响. 近年来随着计算机性能的不断提高,三维波场数值模拟方法得到快速发展,并被广泛用于局部/区域地震波场模拟及壳幔结构成像. 然而由于计算成本过高,实现全球尺度模型的高频地震波场数值模拟存在较大挑战,因而基于远震波场的混合数值模拟方法逐渐得到关注和应用. 混合法将计算区域分解为全球区域和局部区域,在全球区域中采用一维地球模型近似,应用快速算法计算全球高频理论地震图;在局部目标区域内采用三维数值方法(谱元法、有限差分法等)和注入技术,模拟地震波在非均匀介质中的传播,从而在波场计算的效率和精度之间达到平衡. 随着密集台阵观测的普及,对地下结构成像的分辨率提出了更高的要求,准确高效的地震波场混合模拟方法在高分辨率地震成像领域将发挥日益重要的作用. 本文系统地总结了远震波场数值模拟的一维模拟方法,并在此基础上重点介绍混合模拟方法的原理及应用.
  • 图 1  p平面上的Cagniard路径

    Figure 1.  Cagniard-de Hoop contour

    图 2  DWN方法的物理解释. 将单个源替换为以相等间隔L水平分布的无限多个源阵列. 对于与特定激励频率相对应的给定辐射波长λ,弹性能量仅在离散方向θ上辐射(修改自Bouchon, 2003).

    Figure 2.  Physical interpretation of the DWN method. The single source is replaced by an infinite array of sources distributed horizontally at equal interval L. For a given radiation wavelength k corresponding to a specific frequency of excitation, the elastic energy is radiated in discrete directions $ \theta $ only(modified from Bouchon, 2003

    图 3  层析成像获得的从南美到南非的二维速度剖面(修改自Ritsema et al., 1999; Ni et al., 2003

    Figure 3.  2D velocity section from South America to South Africa obtained from tomography (modified from Ritsema et al., 1999; Ni et al., 2003)

    图 4  WKM法(左)和伪谱法(右)生成的SV波合成地震图对比(修改自Ni et al., 2003

    Figure 4.  Comparison of SV synthetics generated by the WKM method against the pseudo-spectral method (on the right)(modified from Ni et al., 2003)

    图 5  底部半空间上的N层组成的分层半空间

    Figure 5.  A layered half-space consists of N layers over a half-space at the bottom.

    图 6  震源分解图样(从左到右为单极源、偶极源、四极源)(修改自Nissen-Meyer et al., 2014

    Figure 6.  Radiation patterns for monopole, dipole, and quadrupole from left to right (modified from Nissen-Meyer et al., 2014)

    图 7  GRT-FD混合方法原理示意图(修改自Zhao et al., 2008

    Figure 7.  Schematic diagram of the principle of hybrid method (modified from Zhao et al., 2008)

    图 8  用于全波形反演的P(上)和S(下)波速度模型(修改自Monteiller et al., 2015

    Figure 8.  Synthetic model of P (top) and S (bottom) velocities for full waveform inversion (modified from Monteiller et al., 2015)

    图 9  由垂直分量直达P波的(a)全波形反演和(b)伴随层析成像得到的模型,在伴随层析成像情况下L-BFG迭代15次,在全波形反演情况下迭代5次(修改自Monteiller et al., 2015

    Figure 9.  Models obtained by (a) full waveform inversion and (b) adjoint tomography from vertical-component direct P waves, after 15 iterations of the L-BFGS in the case of adjoint tomography and 5 iterations in the case of full waveform inversion (modified from Monteiller et al., 2015)

    图 10  (a)在西藏中部部署的Hi-CLIMB台站(三角形)的地理分布. 红色三角形表示用于观测和合成RF剖面的台站. $ \mathrm{B}{\mathrm{B}}^{'} $Tseng等(2009)估算地壳厚度的剖面. (b)Hung等(2011)绘制的沿剖面$ \mathrm{B}{\mathrm{B}}^{'} $$ {V}_{\mathrm{P}} $和(c)$ {V}_{\mathrm{S}} $扰动. (d)Tseng等(2009)根据$ {T}_{\mathrm{SsPmp}}-{T}_{\mathrm{Ss}} $估计的地壳厚度(蓝点). (e)2005年发生在北纬5.32°、东经123.34°、深度522 km的地震的SV波地震反射剖面(垂直分量速度). 各种散射/转换的相位,用紫色箭头表示. (f)基于局部地壳模型的三维SEM-FK混合方法计算相应的SV波合成地震剖面,该模型综合了CRUST1.0、估计的莫霍剖面、三维$ {V}_{\mathrm{P}} $$ {V}_{\mathrm{S}} $变化. 图4e图4f中的所有地震记录都在Ss震相上对齐(修改自Tong et al., 2014b

    Figure 10.  (a) Geographic distributions of Hi-CLIMB stations (triangles) previously deployed in central Tibet. Red triangles indicate stations used to generate the observed and synthetic RF profiles. $ \mathrm{B}{\mathrm{B}}^{'} $ is the profile along which the crust thickness was estimated by Tseng et al. (2009). (b)VP and (c) VS perturbations along Profile $ \mathrm{B}{\mathrm{B}}^{'} $ mapped by Hung et al. (2011) superimposed onto the estimated Moho. (d) Estimated crust thickness (blue dots) based on $ {T}_{\mathrm{SsPmp}}-{T}_{\mathrm{Ss}} $ by Tseng et al. (2009). (e) Observed SV wave seismic reflection profile (vertical-component velocity) for a 2005 earthquake occurred at $ 5.3{2}^{\circ } $ N, $ 123.3{4}^{\circ } $ E, and a depth of 522 km. Various scattered/converted phases, including SsPmp, are indicated by purple arrows. (f) Corresponding synthetic SV wave seismic profile computed by the 3-D SEM-FK hybrid method based on a local crustal model that incorporates CRUST1.0, the estimated Moho profile, and the 3-D VP and VS variations. All seismograms in Figures 4e and 4f are aligned on the Ss phase (modified from Tong et al., 2014b)

  • [1]

    Aki K, Larner K L. 1970. Surface motion of a layered medium having an irregular interface due to incident plane SH waves[J]. Journal of Geophysical Research, 75(5): 933-954. doi: 10.1029/JB075i005p00933
    [2]

    Al-Attar D, Woodhouse J. 2008. Calculation of seismic displacement fields in self-gravitating earth models—applications of minors vectors and symplectic structure[J]. Geophysical Journal International, 175: 1176-1208. doi: 10.1111/j.1365-246X.2008.03961.x
    [3]

    Alterman Z, Karal Jr F. 1968. Propagation of elastic waves in layered media by finite difference methods[J]. Bulletin of the Seismological Society of America, 58(1): 367-398.
    [4]

    Baker B, Roecker S. 2014. A full waveform tomography algorithm for teleseismic body and surface waves in 2.5 dimensions[J]. Geophysical Journal International, 198(3): 1775-1794. doi: 10.1093/gji/ggu236
    [5]

    Beller S, Monteiller V, Combe L, et al. 2018a. On the sensitivity of teleseismic full-waveform inversion to earth parametrization, initial model and acquisition design[J]. Geophysical Journal International, 212(2): 1344-1368. doi: 10.1093/gji/ggx480
    [6]

    Beller S, Monteiller V, Operto S, et al. 2018b. Lithospheric architecture of the South-Western Alps revealed by multiparameter teleseismic full-waveform inversion[J]. Geophysical Journal International, 212(2): 1369-1388. doi: 10.1093/gji/ggx216
    [7]

    Bielak J, Christiano P. 1984. On the effective seismic input for non-linear soil-structure interaction systems[J]. Earthquake Engineering & Structural Dynamics, 12(1): 107-119.
    [8]

    Bielak J, MacCamy R C, McGhee D S, et al. 1991. Unified symmetric BEM-FEM for site effects on ground motion—SH waves[J]. Journal of Engineering Mechanics, 117(10): 2265-2285. doi: 10.1061/(ASCE)0733-9399(1991)117:10(2265)
    [9]

    Bielak J, Loukakis K, Hisada Y, et al. 2003. Domain reduction method for three-dimensional earthquake modeling in localized regions, Part I: Theory[J]. Bulletin of the Seismological Society of America, 93(2): 817-824. doi: 10.1785/0120010251
    [10]

    Bouchon M, Aki K. 1977. Discrete wave-number representation of seismic-source wave fields[J]. Bulletin of the Seismological Society of America, 67(2): 259-277.
    [11]

    Bouchon M. 1979. Discrete Wave number representation of elastic wave fields in 3-space dimensions[J]. Journal of Geophysical Research, 84(Nb7): 3609-3614. doi: 10.1029/JB084iB07p03609
    [12]

    Bouchon M. 1981. A simple method to calculate Green's functions for elastic layered media[J]. Bulletin of the Seismological Society of America, 71(4): 959-971.
    [13]

    Bouchon M. 2003. A review of the discrete wavenumber method[J]. Pure and Applied Geophysics, 160(3-4): 445-465.
    [14]

    Capdeville Y, Chaljub E, Vilotte J P, et al. 2003a. Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models[J]. Geophysical Journal International, 152(1): 34-67. doi: 10.1046/j.1365-246X.2003.01808.x
    [15]

    Capdeville Y, To A, Romanowicz B. 2003b. Coupling spectral elements and modes in a spherical Earth: an extension to the 'sandwich' case[J]. Geophysical Journal International, 154(1): 44-57. doi: 10.1046/j.1365-246X.2003.01959.x
    [16]

    Carcione J M, Kosloff D, Kosloff R. 1988. Wave-Propagation Simulation in a Linear Viscoacoustic Medium[J]. Geophysical Journal, 93(2): 393-407. doi: 10.1111/j.1365-246X.1988.tb02010.x
    [17]

    Chapman C H. 1976. Exact and approximate generalized ray theory in vertically inhomogeneous-media[J]. Geophysical Journal of the Royal Astronomical Society, 46(2): 201-233. doi: 10.1111/j.1365-246X.1976.tb04154.x
    [18]

    Chapman C H. 1978. A new method for computing synthetic seismograms[J]. Geophysical Journal International, 54(3): 481-518. doi: 10.1111/j.1365-246X.1978.tb05491.x
    [19]

    Chapman C H, Drummond R. 1982. Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory[J]. Bulletin of the Seismological Society of America, 72(6B): S277-S317.
    [20]

    Chapman C H, Orcutt J A. 1985. The computation of body wave synthetic Seismograms in laterally homogeneous media[J]. Reviews of Geophysics, 23(2): 105-163. doi: 10.1029/RG023i002p00105
    [21]

    Chen L, Wen L X, Zheng T Y. 2005. A wave equation migration method for receiver function imaging: 2. Application to the Japan subduction zone[J]. Journal of Geophysical Research: Solid Earth, 110(B11): B11309.
    [22]

    Chen X F. 1990. Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method. I. Theory of two-dimensional SH case[J]. Bulletin of the Seismological Society of America, 80(6A): 1696-1724.
    [23]

    Chen X F, Zhang H M. 2001. An efficient method for computing Green's functions for a layered half-space at large epicentral distances[J]. Bulletin of the Seismological Society of America, 91(4): 858-869.
    [24]

    Chevrot S, Zhao L. 2007. Multiscale finite-frequency Rayleigh wave tomography of the Kaapvaal craton[J]. Geophysical Journal International, 169(1): 201-215. doi: 10.1111/j.1365-246X.2006.03289.x
    [25]

    Clouzet P, Masson Y, Romanowicz B. 2018. Box tomography: First application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[J]. Geophysical Journal International, 213(3): 1849-1875. doi: 10.1093/gji/ggy078
    [26]

    Cummins P R, Geller R J, Takeuchi N. 1994. Dsm complete synthetic seismograms-P-Sv, spherically symmetrical, case[J]. Geophysical Research Letters, 21(15): 1663-1666. doi: 10.1029/94GL01281
    [27]

    Dahlen F, Tromp J. 1998. Theoretical Global Seismology[M]. New Jersey: Princeton University Press.
    [28]

    Dey-Sarkar S, Chapman C. 1978. A simple method for computation of body-wave seismograms[J]. Bulletin of the Seismological Society of America, 68(6): 1577-1593.
    [29]

    Dmitry B, Singh S C, Nobuaki F. 2015. An efficient method of 3-D elastic full waveform inversion using a finite-difference injection method for time-lapse imaging[J]. Geophysical Journal International, 202(3), 1908-1922. doi: 10.1093/gji/ggv268
    [30]

    Fäh D, Suhadolc P, Panza G F. 1993. Variability of seismic ground motion in complex media: the case of a sedimentary basin in the Friuli (Italy) area[J]. Journal of Applied Geophysics, 30(1-2): 131-148. doi: 10.1016/0926-9851(93)90022-Q
    [31]

    Fäh D, Suhadolc P, Mueller S, et al. 1994. A hybrid method for the estimation of ground motion in sedimentary basins: quantitative modeling for Mexico City[J]. Bulletin of the Seismological Society of America, 84(2): 383-399.
    [32]

    Florsch N, Fäh D, Suhadolc P, et al. 1991. Complete synthetic seismograms for high-frequency multimode SH-waves[J]. Pure and Applied Geophysics, 136(4): 529-560. doi: 10.1007/BF00878586
    [33]

    Fu L-Y, Bouchon M. 2004. Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media-I. Theory of two-dimensional SH case[J]. Geophysical Journal International, 157(2): 481-498. doi: 10.1111/j.1365-246X.2004.02135.x
    [34]

    Geller R J, Ohminato T. 1994. Computation of synthetic seismograms and their partial derivatives for heterogeneous media with arbitrary natural boundary-conditions using the direct solution method[J]. Geophysical Journal International, 116(2): 421-446. doi: 10.1111/j.1365-246X.1994.tb01807.x
    [35]

    Geller R J, Takeuchi N. 1995. A new method for computing highly accurate Dsm synthetic seismograms[J]. Geophysical Journal International, 123(2): 449-470. doi: 10.1111/j.1365-246X.1995.tb06865.x
    [36]

    Gilbert F. 1971. Excitation of the normal modes of the earth by earthquake sources[J]. Geophysical Journal of the Royal Astronomical Society, 22: 223-226. doi: 10.1111/j.1365-246X.1971.tb03593.x
    [37]

    Grand S P. 1994. Mantle shear structure beneath the Americas and surrounding oceans[J]. Journal of Geophysical Research Solid Earth, 99: 11591-11622. doi: 10.1029/94JB00042
    [38]

    Graves R W, Helmberger D V. 1988. Upper mantle cross section from Tonga to Newfoundland[J]. Journal of Geophysical Research Atmospheres, 93(B5): 4701-4711. doi: 10.1029/JB093iB05p04701
    [39]

    Gualtieri L, Stutzmann E, Capdeville Y, et al. 2013. Modelling secondary microseismic noise by normal mode summation[J]. Geophysical Journal International, 193(3): 1732-1745. doi: 10.1093/gji/ggt090
    [40]

    Haskell N A. 1953. The dispersion of surface waves on multilayered media[J]. Bulletin of the Seismological Society of America, 43(1): 17-34.
    [41]

    Haskell N A. 1964. Total energy and energy spectral density of elastic wave radiation from propagating faults[J]. Bulletin of the Seismological Society of America, 54(6A): 1811-1841.
    [42]

    Helmberger D V. 1968. The crust-mantle transition in Bering Sea[J]. Bulletin of the Seismological Society of America, 58(1): 179-214.
    [43]

    Helmberger D V. 1974. Generalized ray theory for shear dislocations[J]. Bulletin of the Seismological Society of America, 64(1): 45-64.
    [44]

    Helmberger D V, Harkrider D G. 1978. Modeling earthquakes with generalized ray theory[J]. Air Force Office of Scientific Research, Bolling AFB, DC, 21:174-222.
    [45]

    Helmberger D V. 1983. Theory and application of synthetic seismograms[J]. Earthquakes: Observation, Theory and Interpretation, 85: 174-222.
    [46]

    Hung S H, Chen W P, Chiao L Y. 2011. A data-adaptive, multiscale approach of finite-frequency, traveltime tomography with special reference to P and S wave data from central Tibet[J]. Journal of Geophysical Research: Solid Earth, 116(B6): B06307.
    [47]

    Ivo O, Johana B, Donat F, et al. 2002. 3D hybrid Ray-FD and DWN-FD seismic modeling for simple models containing complex local structures[J]. Studia Geophysica et Geodaetica, 46(4): 711-730. doi: 10.1023/A:1021181422709
    [48]

    Kawai K, Takeuchi N, Geller R J. 2006. Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media[J]. Geophysical Journal International, 164(2): 411-424. doi: 10.1111/j.1365-246X.2005.02829.x
    [49]

    Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation[J]. Geophysical Journal International, 139(3): 806-822. doi: 10.1046/j.1365-246x.1999.00967.x
    [50]

    Komatitsch D, Liu Q Y, Tromp J, et al. 2004. Simulation of ground motion in the Los Angeles Basin base upon the Spectral-Element Method[J]. Bulletin of Seismological Society of America, 94(1): 187-206 doi: 10.1785/0120030077
    [51]

    Kummer B, Behle A, Dorau F. 1987. Hybrid modeling of elastic-wave propagation in two-dimensional laterally inhomogeneous media[J]. Geophysics, 52(6), 765-771. doi: 10.1190/1.1442343
    [52]

    Li D Z, Helmberger D, Clayton R W, et al. 2014. Global synthetic seismograms using a 2-D finite-difference method[J]. Geophysical Journal International, 197(2): 1166-1183.
    [53] 李旭, 陈运泰. 1996. 合成地震图的广义反射透射系数矩阵方法[J]. 地震地磁观测与研究, 17(3): 1-20.

    Li X, Chen Y T. 1996. The generalized reflection-transmission coefficient matrix method for synthetic seismograms[J]. Seismological and Geomagnetic Observation and Research, 17(3): 1-20(in Chinese).
    [54] 李幼铭, 孙永智, 姚振兴. 1980. 关于“反射法”计算综合地震图的实施方案[J]. 地球物理学报, 23(4): 389-395. doi: 10.3321/j.issn:0001-5733.1980.04.005

    Li Y M, Sun Y Z, Yao Z X. 1980. A practical scheme on the computation of synthetic seismograms using “reflectivity method” [J]. Chinese Journal of Geophysics, 23(4): 389-395(in Chinese). doi: 10.3321/j.issn:0001-5733.1980.04.005
    [55]

    Lin C X, Monteiller V, Wang K, et al. 2019. High-frequency seismic wave modelling of the deep Earth based on hybrid methods and spectral-element simulations: a conceptual study[J]. Geophysical Journal International, 219(3): 1948-1969. doi: 10.1093/gji/ggz413
    [56]

    Liu Q, Gu Y J. 2012. Seismic imaging: From classical to adjoint tomography[J]. Tectonophysics, 566-567: 31-66. doi: 10.1016/j.tecto.2012.07.006
    [57]

    Liu S L, Yang D H, Dong X P, et al. 2017a. Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling[J]. Solid Earth, 8(5): 969-986. doi: 10.5194/se-8-969-2017
    [58]

    Liu S L, Yang D H, Ma J. 2017b. A modified symplectic PRK scheme for seismic wave modeling[J]. Computers & Geosciences, 99: 28-36.
    [59]

    Liu T, Zhang H. 2017. Synthetic seismograms for finite sources in spherically symmetric Earth using normal-mode summation[J]. Earthquake Science, 30(3): 125-133. doi: 10.1007/s11589-017-0188-1
    [60] 刘有山, 滕吉文, 刘少林, 徐涛. 2013. 稀疏存储的显式有限元三角网格地震波数值模拟及其PML吸收边界条件. 地球物理学报, 56(9): 3085-3099. doi: 10.6038/cjg20130921

    Liu Y S, Teng J W, Liu S L, Xu T. 2013. Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition[J]. Chinese Journal of Geophysics, 56(9): 3085-3099(in Chinese). doi: 10.6038/cjg20130921
    [61] 刘有山, 滕吉文, 徐涛, 等. 2014. 三角网格谱元法地震波场数值模拟[J]. 地球物理学进展, 29(4): 1715-1726. doi: 10.6038/pg20140430

    Liu Y S, Teng J W, Xu T, et al. 2014. Numerical modeling of seismic wavefield with the SEM based on Triangles[J]. Progress in Geophysics, 29(4): 1715-1726(in Chinese). doi: 10.6038/pg20140430
    [62]

    Liu Y S, Teng J W, Xu T, Badal J. 2017c. Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling[J]. Journal of Computational Physics, 336: 458-480. doi: 10.1016/j.jcp.2017.01.069
    [63]

    Lysmer J, Drake L A. 1972. A Finite Element Method for Seismology[M]//Bolt B A. Methods in Computational Physics: Advances in Research and Applications, 11: 181-216.
    [64] 栾威, 申文斌, 丁浩. 2021. 地球自由振荡弹性简正模研究进展与展望[J]. 地球与行星物理论评, 52(3): 308-325.

    Luan W, Shen W B, Ding H. 2021. Progress and prospect of studies on elastic normal modes of Earth's free oscillation[J]. Reviews of Geophysics and Planetary Physics, 52(3): 308-325(in Chinese).
    [65]

    Masson Y, Cupillard P, Capdeville Y, et al. 2014. On the numerical implementation of time-reversal mirrors for tomographic imaging[J]. Geophysical Journal International, 196(3): 1580-1599.
    [66]

    Masson Y, Romanowicz B. 2017a. Box tomography: localized imaging of remote targets buried in an unknown medium, a step forward for understanding key structures in the deep Earth[J]. Geophysical Journal International, 211(1): 141-163. doi: 10.1093/gji/ggx141
    [67]

    Masson Y, Romanowicz B. 2017b. Fast computation of synthetic seismograms within a medium containing remote localized perturbations: a numerical solution to the scattering problem[J]. Geophysical Journal International, 208(2): 674-692. doi: 10.1093/gji/ggw412
    [68]

    Mita A, Luco J E. 1987. Dynamic-response of embedded foundations-a hybrid approach[J]. Computer Methods in Applied Mechanics and Engineering, 63(3): 233-259. doi: 10.1016/0045-7825(87)90071-5
    [69]

    Moczo P, Bystrický E, Kristek J, et al. 1997. Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[J]. Bulletin of the Seismological Society of America, 87(5): 1305-1323.
    [70]

    Monteiller V, Chevrot S, Komatitsch D, et al. 2013. A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[J]. Geophysical Journal International, 192(1): 230-247. doi: 10.1093/gji/ggs006
    [71]

    Monteiller V, Chevrot S, Komatitsch D, et al. 2015. Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM–DSM hybrid method[J]. Geophysical Journal International, 202(2): 811-827. doi: 10.1093/gji/ggv189
    [72]

    Monteiller V, Beller S, Plazolles B, et al. 2020. On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[J]. Geophysical Journal International, 224(3): 2060-2076.
    [73]

    Ni S D, Ding X M, Helmberger D V. 2000. Constructing synthetics from deep earth tomographic models[J]. Geophysical Journal International, 140(1): 71-82. doi: 10.1046/j.1365-246x.2000.00982.x
    [74]

    Ni S D, Cormier V F, Helmberger D V. 2003. A comparison of synthetic seismograms for 2D structures: Semianalytical versus numerical[J]. Bulletin of the Seismological Society of America, 93(6): 2752-2757. doi: 10.1785/0120030011
    [75]

    Nissen-Meyer T, Dahlen F A, Fournier A. 2007a. Spherical-earth frechet sensitivity kernels[J]. Geophysical Journal International, 168(3): 1051-1066. doi: 10.1111/j.1365-246X.2006.03123.x
    [76]

    Nissen-Meyer T, Fournier A, Dahlen F A. 2007b. A two-dimensional spectral-element method for computing spherical-earth seismograms - I. Moment-tensor source[J]. Geophysical Journal International, 168(3): 1067-1092. doi: 10.1111/j.1365-246X.2006.03121.x
    [77]

    Nissen-Meyer T, Fournier A, Dahlen F A. 2008. A 2-D spectral-element method for computing spherical-earth seismograms-II. Waves in solid-fluid media[J]. Geophysical Journal International, 174(3): 873-888. doi: 10.1111/j.1365-246X.2008.03813.x
    [78]

    Nissen-Meyer T, van Driel M, Stahler S C, et al. 2014. AxiSEM: broadband 3-D seismic wavefields in axisymmetric media[J]. Solid Earth, 5(1): 425-445. doi: 10.5194/se-5-425-2014
    [79]

    Nolet G. 2008, A Breviary of Seismic Tomography[M]. Cambridge: Cambridge University Press.
    [80]

    Opršal I, Matyska C, Irikura K. 2010. The source-box wave propagation hybrid methods: general formulation and implementation[J]. Geophysical Journal of the Royal Astronomical Society, 176(2): 555-564.
    [81]

    Paulssen H. 1988. Evidence For Small Scale Structure of The Upper Mantle[M]. Instituut voor Aardwetenschappen der Rijksunversiteit te Utrecht: Geologica Ultraiectina.
    [82]

    Rayleigh Baron J W S. 1896. The Theory of Sound (Vol. 2) [M]. Macmillan.
    [83]

    Rayleigh L. 1907. On the dynamical theory of gratings[J]. Proceedings of the Royal Society of London Series a-Containing Papers of a Mathematical and Physical Character, 79(532): 399-416.
    [84]

    Ritsema J, Van H H J, Woodhouse J. 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland[J]. Science 286: 1925–1928. doi: 10.1126/science.286.5446.1925
    [85]

    Ritsema J, Deuss A, Van H H J, et al. 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements[J]. Geophysical Journal International, 184(3): 1223-1236. doi: 10.1111/j.1365-246X.2010.04884.x
    [86]

    Robertsson J O A, Chapman C H. 2000. An efficient method for calculating finite-difference seismograms after model alterations[J]. Geophysics, 65(3): 907-918. doi: 10.1190/1.1444787
    [87]

    Roecker S, Baker B, McLaughlin J. 2010. A finite-difference algorithm for full waveform teleseismic tomography[J]. Geophysical Journal International, 181(2): 1017-1040.
    [88]

    Schmidt H, Tango G. 1986. Efficient global matrix approach to the computation of synthetic seismograms[J]. Geophysical Journal of the Royal Astronomical Society, 84(2): 331-359. doi: 10.1111/j.1365-246X.1986.tb04359.x
    [89]

    Shtivelman V. 1984. A hybrid method for wave field computation[J]. Geophysical prospecting, 32(2): 236-257. doi: 10.1111/j.1365-2478.1984.tb00730.x
    [90]

    Shtivelman V. 1985. Two-dimensional acoustic modeling by a hybrid method[J]. Geophysics, 50(8): 1273-1284. doi: 10.1190/1.1441998
    [91]

    Singh S J, Ben-Menahem A. 1969a. Eigenvibrations of the earth excited by finite dislocations-II Spheroidal oscillations[J]. Geophysical Journal of the Royal Astronomical Society, 17(3): 333-350. doi: 10.1111/j.1365-246X.1969.tb00243.x
    [92]

    Singh S J, Ben-Menahem A. 1969b. Eigenvibrations of the earth excited by finite dislocations—I Toroidal oscillations[J]. Geophysical Journal International,17(2): 151-177. doi: 10.1111/j.1365-246X.1969.tb02317.x
    [93]

    Stead R J, Helmberger D V. 1988. Numerical-analytical Interfacing in Two Dimensions with Applications to Modeling NTS Seismograms[M]//Aki K, Wu R-S. Scattering and Attenuations of Seismic Waves, Part I. Basel: Birkhäuser Basel, 157-193.
    [94]

    Taflove A, Hagness S. 2005. Computational electrodynamics: The finite-difference time-domain method[J]. Artech house, 5(15): 629–670.
    [95]

    Takeuchi H, Saito M. 1972. Seismic Surface Waves[M]//Bolt B A. Methods in Computational Physics: Advances in Research and Applications. Elsevier, 11: 217-295.
    [96]

    Takeuchi N, Geller R J, Cummins P R. 1996. Highly accurate P-SV complete synthetic seismograms using modified DSM operators[J]. Geophysical Research Letters, 23(10): 1175-1178. doi: 10.1029/96GL00973
    [97]

    Takeuchi N, Geller R. 2003. Accurate numerical methods for solving the elastic equation of motion for arbitrary source locations[J]. Geophysical Journal International, 154: 852-866. doi: 10.1046/j.1365-246X.2003.02009.x
    [98]

    Tanimoto T. 1984. A simple derivation of the formula to calculate synthetic long-period seismograms in a heterogeneous earth by normal mode summation[J]. Geophysical Journal of the Royal Astronomical Society, 77(1): 275-278. doi: 10.1111/j.1365-246X.1984.tb01934.x
    [99]

    Thomson W T. 1950. Transmission of elastic waves through a stratified solid medium[J]. Journal of Applied Physics, 21(2): 89-93. doi: 10.1063/1.1699629
    [100]

    Tong P, Chen C-W, Komatitsch D, et al. 2014a. High-resolution seismic array imaging based on an SEM-FK hybrid method[J]. Geophysical Journal International, 197(1): 369-395. doi: 10.1093/gji/ggt508
    [101]

    Tong P, Komatitsch D, Tseng T-L, et al. 2014b. A 3-D spectral-element and frequency-wave number hybrid method for high-resolution seismic array imaging[J]. Geophysical Research Letters, 41(20): 7025-7034. doi: 10.1002/2014GL061644
    [102]

    Tseng T L, Chen W P, Nowack R L. 2009. Northward thinning of Tibetan crust revealed by virtual seismic profiles[J]. Geophysical Research Letters, 36: L24304. doi: 10.1029/2009GL040457
    [103]

    Turcotte D L, Schubert G. 2002. Geodynamics[M]. Cambridge: Cambridge University Press.
    [104]

    Van Manen D J, Robertsson J O A, Curtis A. 2007. Exact wave field simulation for finite-volume scattering problems[J]. The Journal of the Acoustical Society of America, 122(4): EL115-EL121. doi: 10.1121/1.2771371
    [105] 王椿镛. 1982. 层状不均匀介质中合成地震图的反射率法[J]. 地球物理学报, 25(5): 424-433. doi: 10.3321/j.issn:0001-5733.1982.05.005

    Wang C Y. 1982. The reflectivity method of synthetic seismogram calculation for layered inhomogeneous media[J]. Chinese Journal of Geophysics, 25(5): 424-433(in Chinese). doi: 10.3321/j.issn:0001-5733.1982.05.005
    [106]

    Wang Y, Chevrot S, Monteiller V, et al. 2016. The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P waves[J]. Geology, 44(6): 475-478. doi: 10.1130/G37812.1
    [107] 温联星, 姚振兴. 1994. 用广义射线和有限差分计算近场理论地震图的混合方法[J]. 地球物理学报, 37(2): 211-219. doi: 10.3321/j.issn:0001-5733.1994.02.009

    Wen L X, Yao Z X. 1994. Theory and application of hybrid method seismograms[J]. Chinese Journal of Geophysics, 37(2): 211-219(in Chinese). doi: 10.3321/j.issn:0001-5733.1994.02.009
    [108]

    Wen L X, Helmberger D V. 1998a. A two-dimensional P-SV hybrid method and its application to modeling localized structures near the core-mantle boundary[J]. Journal of Geophysical Research: Solid Earth, 103(B8): 17901-17918. doi: 10.1029/98JB01276
    [109]

    Wen L X, Helmberger D V. 1998b. Ultra-low velocity zones near the core-mantle boundary from broadband PKP precursors[J]. Science, 279(5357): 1701-1703. doi: 10.1126/science.279.5357.1701
    [110]

    Wen L X. 2002. An SH hybrid method and shear velocity structures in the lowermost mantle beneath the central Pacific and South Atlantic Oceans[J]. Journal of Geophysical Research-Solid Earth, 107(B3): 2055. doi: 10.1029/2001JB000499
    [111]

    Wen L X, Niu F L. 2002. Seismic velocity and attenuation structures in the top of the Earth's inner core[J]. Journal of Geophysical Research, 107(B11): ESE 2-1-ESE 2-13. doi: 10.1029/2001JB000170
    [112]

    Wiggins R A. 1976. A fast, new computational algorithm for free oscillations and surface waves[J]. Geophysical Journal International, 47(1): 135-150. doi: 10.1111/j.1365-246X.1976.tb01266.x
    [113]

    Woodhouse J H. 1988. The calculation of the eigenfrequencies and eigenfunctions of the free oscillations of the Earth and the Sun[J]. Seismological Algorithms, 230: 321-370.
    [114]

    Wu S L, Nozu A, Nagasaka Y. 2020. Accuracy of near-fault fling-step displacements estimated using the discrete wavenumber method[J]. Bulletin of the Seismological Society of America, doi: 10.1785/0120190257.
    [115]

    Wu W, Ni S, Zhan Z, et al. 2018. An SEM-DSM three-dimensional hybrid method for modelling teleseismic waves with complicated source-side structures[J]. Geophysical Journal International, 215(1): 133-154. doi: 10.1093/gji/ggy273
    [116] 谢小碧, 郑天愉, 姚振兴. 1992. 理论地震图计算方法[J].地球物理学报, 35(6): 790-801. doi: 10.3321/j.issn:0001-5733.1992.06.014

    Xie X B, Zheng T Y, Yao Z X. 1992. Methods of synthetic seismograms——a review[J]. Chinese Journal of Geophysics, 35(6): 790-801(in Chinese). doi: 10.3321/j.issn:0001-5733.1992.06.014
    [117]

    Yang H-Y, Zhao L, Hung S-H. 2010. Synthetic seismograms by normal-mode summation: a new derivation and numerical examples[J]. Geophysical Journal International, 183(3): 1613-1632. doi: 10.1111/j.1365-246X.2010.04820.x
    [118] 姚振兴. 1979. 层状介质、非轴对称震源情况下的反射法[J]. 地球物理学报, 22(2): 181-194. doi: 10.3321/j.issn:0001-5733.1979.02.006

    Yao Z X. 1979. Generalized reflection coefficients for a layered medium and asymmetrical source[J]. Chinese Journal of Geophysics, 22(2): 181-194(in Chinese). doi: 10.3321/j.issn:0001-5733.1979.02.006
    [119] 姚振兴. 1980. 剪切位错内源的反射法[J]. 地球物理学报, 23(3): 306-311, 343-344. doi: 10.3321/j.issn:0001-5733.1980.03.008

    Yao Z X. 1980. The reflectivity method for a buried shear dislocation source[J]. Chinese Journal of Geophysics, 23(3): 306-311, 343-344(in Chinese). doi: 10.3321/j.issn:0001-5733.1980.03.008
    [120] 姚振兴, 郑天愉. 1994. 地震波传播理论和应用研究[J]. 地球物理学报, 37(A01): 160-171.

    Yao Z X, Zheng T Y. 1994. Theory of seismic wave propagation and its application[J]. Chinese Journal of Geophysics, 37(A01): 160-171(in Chinese).
    [121]

    Yoshimura C, Bielak J, Hisada Y, et al. 2003. Domain reduction method for three-dimensional earthquake modeling in localized regions, part II: Verification and applications[J]. Bulletin of the Seismological Society of America, 93(2): 825-841. doi: 10.1785/0120010252
    [122]

    Zahradnik J, Moczo P. 1996. Hybrid seismic modeling based on discrete-wave number and finite-difference methods[J]. Pure and Applied Geophysics, 148(1-2): 21-38. doi: 10.1007/BF00882053
    [123]

    Zhang H M, Chen X F, Chang S. 2003. An efficient numerical method for computing synthetic seismograms for a layered half-space with sources and receivers at close or same depths[J]. Pure Applied Geophysics, 160(3): 467-486.
    [124] 张明辉, 刘有山, 侯爵, 等. 2019. 近地表地震层析成像方法综述[J]. 地球物理学进展, 34(1): 48-63. doi: 10.6038/pg2019CC0534

    Zhang M H, Liu Y S, Hou J, et al. 2019. Review of seismic tomography methods in near surface structures reconstruction[J]. Progress in Geophysics, 34(1): 48-63(in Chinese). doi: 10.6038/pg2019CC0534
    [125]

    Zhao L, Wen L X, Chen L, et al. 2008. A two-dimensional hybrid method for modeling seismic wave propagation in anisotropic media[J]. Journal of Geophysical Research, 113(B12): B12307. doi: 10.1029/2008JB005733
    [126]

    Zhou H, Chen X F. 2006. A new approach to simulate scattering of SH waves by an irregular topography[J]. Geophysical Journal International, 164(2): 449-459. doi: 10.1111/j.1365-246X.2005.02670.x
    [127]

    Zhou H, Chen X F. 2008. The localized boundary integral equation discrete wavenumber method for simulating P-SV wave scattering by an irregular topography[J]. Bulletin of the Seismological Society of America, 98(1): 265-279. doi: 10.1785/0120060249
    [128] 朱良保. 1989. 计算横向非均匀介质体波理论地震图的Maslov渐近法[J]. 地球物理学报, 32(2): 163-172. doi: 10.3321/j.issn:0001-5733.1989.02.005

    Zhu L B. 1989. Maslov's asymptotic method of calculating the theoretical seismogram of body waves in the lateral inhomogeneous medium[J]. Chinese Journal of Geophysics, 32(02): 163-172(in Chinese). doi: 10.3321/j.issn:0001-5733.1989.02.005
    [129]

    Zhu L P, Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophysical Journal International, 148(3): 619-627. doi: 10.1046/j.1365-246X.2002.01610.x
    [130] 朱仁益. 1990. 用 WKBJ 法计算衰减介质中的理论地震图[J]. 地球物理学报, (增刊): 212-224.

    Zhu R Y. 1990. WKBJ method of theoretical seismogram calculation for attenuation media[J]. Chinese Journal of Geophysics, (Sup.): 212-224 (in Chinese).
  • [1] 吴忠良王龙李丽张晓东邵志刚李营孙珂车时 . 中国地震科学实验场:地震预测与系统设计. 地球与行星物理论评, doi: 10.16738/j.dqyxx.2021-028
    [2] 孙伟家王一博魏勇赵亮 . 火星地震学与内部结构研究. 地球与行星物理论评, doi: 10.16738/j.dqyxx.2021-016
  • 加载中
图(10)
计量
  • 文章访问数:  96
  • HTML全文浏览量:  64
  • PDF下载量:  14
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-03-09
  • 网络出版日期:  2021-04-19

远震波场正演模拟方法及应用

    作者简介: 桑莹泉(1998-),男,硕士生,主要从事地震学研究. E-mail:sangyq@mail.iggcas.ac.cn

摘要: 地震波场数值模拟是壳幔结构成像和深部探测的重要基础. 经典的远震波场模拟主要基于一维地球模型,包括解析法、半解析法和数值法等. 这些算法能够高效地计算理论地震图,但由于将地球假设为一维层状介质,难以考虑介质的横向非均匀性对地震波场的影响. 近年来随着计算机性能的不断提高,三维波场数值模拟方法得到快速发展,并被广泛用于局部/区域地震波场模拟及壳幔结构成像. 然而由于计算成本过高,实现全球尺度模型的高频地震波场数值模拟存在较大挑战,因而基于远震波场的混合数值模拟方法逐渐得到关注和应用. 混合法将计算区域分解为全球区域和局部区域,在全球区域中采用一维地球模型近似,应用快速算法计算全球高频理论地震图;在局部目标区域内采用三维数值方法(谱元法、有限差分法等)和注入技术,模拟地震波在非均匀介质中的传播,从而在波场计算的效率和精度之间达到平衡. 随着密集台阵观测的普及,对地下结构成像的分辨率提出了更高的要求,准确高效的地震波场混合模拟方法在高分辨率地震成像领域将发挥日益重要的作用. 本文系统地总结了远震波场数值模拟的一维模拟方法,并在此基础上重点介绍混合模拟方法的原理及应用.

English Abstract

    • 解析法是基于高频近似合成地震记录非常有效的手段,主要包括广义射线法、WKBJ和简正模法等. 解析法计算效率较高,但对于有限频率,解析法不够精确,无法提供相比于波动方程的精确或完备解. 解析法的精度足以解决许多三维地震波传播的实际问题,而且这些问题很难用其他方法高效解决.

      (1)广义射线法

      广义射线理论(GRT)是基于Laplace变换的Cagniard-de Hoop方法,由于远震波可以通过初动近似来计算,因此分层介质中每一层中的射线初动(包括反射波和首波)在时间域都可以得到解析解(Helmberger, 1974; Helmberger and Harkrider, 1978). 广义射线法适用于计算水平层状或近似水平层状介质模型情况下远场、中低频等某些震相,对于层数较多或近场高频情况,由于转换波和多次反射、折射波等,计算量将大大增加(李旭和陈运泰,1996).

      z轴朝下的柱坐标系$ \left(r,\theta,z\right) $中,点源在拉普拉斯域的波动方程可以写为:

      $ \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \phi }{\partial r}\right)+\frac{1}{{r}^{2}}\frac{{\partial }^{2}\phi }{\partial {\theta }^{2}}+\frac{{\partial }^{2}\phi }{\partial {z}^{2}}-\frac{{s}^{2}}{{c}^{2}}\phi =0 $

      (1)

      利用分离变量法,可以得到点源情况下全空间下的通解:

      $ \phi (r,\theta,z,s)=\sum\nolimits_{n=0,\pm 1,\cdots }{\mathrm{e}}^{-in\theta }{\int }_{0}^{\infty }{A}_{n}\left(k\right){J}_{n}\left(kr\right){\mathrm{e}}^{-\nu \left|z\right|}{\rm{d}}k $

      (2)

      式中:

      $ \nu =\sqrt{{k}^{2}+\frac{{s}^{2}}{{c}^{2}}} $

      (3)

      式中,$ {A}_{n} $由初始条件决定.

      对于液体全空间中的爆炸源,为得到压力场$ \phi \left(t\right) $的时间域解,需要两次积分,一个关于k(不含$ -\dfrac{{M}_{0}}{4{\text{π}}\rho {c}^{2}} $):

      $ \phi \left(s\right)=\frac{1}{s}{\int }_{0}^{\infty }\frac{k}{\nu }{J}_{0}\left(kr\right){\mathrm{e}}^{-\nu \left|z\right|}{\rm{d}}k $

      (4)

      另一个是对$ \phi \left(s\right) $进行拉普拉斯逆变换.

      Cagniard-de Hoop方法可以同时进行双重积分. 通过变量替换$ k=-isp $,利用:

      $ {J}_{0}(-ispr)=\frac{i}{\text{π}}\left({K}_{0}\left(spr\right)-{K}_{0}(-spr)\right) $

      (5)

      式(4)变为:

      $\begin{array}{l} \phi \left(s\right)=-\dfrac{i}{\text{π}}{\displaystyle\int }_{-i\infty }^{i\infty }\dfrac{p}{\eta }{K}_{0}\left(spr\right){\mathrm{e}}^{-s\eta \left|z\right|}{\rm{d}}p \\ \quad\quad =\dfrac{2}{\text{π}}\mathfrak{I}{\displaystyle\int }_{0}^{i\infty }\dfrac{p}{\eta }{K}_{0}\left(spr\right){\mathrm{e}}^{-s\eta \left|z\right|}{\rm{d}}p \end{array}$

      (6)

      式中$ \Im$表示取虚部,且:

      $ \eta =\sqrt{\frac{1}{{c}^{2}}-{p}^{2}} $

      (7)

      x较大时,式(6)变为(为了简便只取第一项):

      $ \phi \left(s\right)=\sqrt{\frac{2}{{ \text{π} }rs}}\mathfrak{I}{\int }_{0}^{i\infty }\frac{\sqrt{p}}{\eta }{\mathrm{e}}^{-s(pr+\eta |z\left|\right)}{\rm{d}}p $

      (8)

      为了写成拉普拉斯变换的形式,令:

      $ t=pr+\eta \left|z\right| $

      (9)

      式中,t为正实数. 因而问题变为,在式(9)中$ p=p\left(t\right) $应有图1所示的积分路径C,此积分路径称为Cagniard路径(图1).

      图  1  复p平面上的Cagniard路径

      Figure 1.  Cagniard-de Hoop contour

      将正实数t作为一个参数,它可从零变至无穷大,那么由式(9),可解得:

      $ p\left(t\right)=\left\{\begin{array}{c}\dfrac{rt}{{R}^{2}}-\dfrac{\sqrt{{t}_{0}^{2}-{t}^{2}}}{{R}^{2}}\left|z\right|,t<{t}_{0},\\ \dfrac{rt}{{R}^{2}}+i\dfrac{\sqrt{{t}^{2}-{t}_{0}^{2}}}{{R}^{2}}\left|z\right|,t>{t}_{0},\end{array}\right. $

      (10)

      式中,$ R=\sqrt{{r}^{2}+{z}^{2}} $$ {t}_{0}=\dfrac{R}{c} $. 由于除了沿实数p$ p>\dfrac{1}{c} $,积分函数是解析的,因此积分路径$ [0,i\infty ] $可以变成沿着Cagniard路径:

      $ J\left(t\right)={\mathcal{L}}^{-1}\left\{\mathfrak{I}{\int }_{0}^{i\infty }\frac{\sqrt{p}}{\eta }{\mathrm{e}}^{-s(pr+\eta |z\left|\right)}{\rm{d}}p\right\}=\mathfrak{I}\left\{\frac{\sqrt{p}}{\eta }\frac{\mathrm{d}p}{\text{d}t}\right\} $

      (11)

      因此压力场的解为:

      $ \phi \left(t\right)=\frac{1}{\text{π}}\sqrt{\frac{2}{r}}\frac{1}{\sqrt{t}}*J\left(t\right) $

      (12)

      由于$ \dfrac{\mathrm{d}p}{\mathrm{d}t} $$ t={t}_{0} $时对$ J\left(t\right) $的贡献最大,也就是初至. 利用其在$ {t}_{0} $处的近似,可以得到初动近似的解:

      $ \phi \left(t\right)=\frac{1}{{\text{π}}R}\frac{1}{\sqrt{t}}*\frac{H\left(t-{t}_{0}\right)}{\sqrt{t-{t}_{0}}}=\frac{H\left(t-{t}_{0}\right)}{R} $

      (13)

      这与精确解相同.

      (2)WKBJ法

      WKBJ方法起源于WKB理论,该理论为线性微分方程提供了一个近似解. WKB是由三位著名学者Wentzel、Kramer和Brillouin的姓名首字母命名. 在地震学中,由于Jeffreys(还有Rayleigh)对它的早期发展做出了贡献,通常被称为WKBJ方法. 最先由Chapman(1976)以及Wiggins(1976)用不同的方法得到了WKBJ地震记录. 后来,Chapman(1978)推广了这一结果. Dey-Sarka和Chapman(1978)给出了该方法的一些数值例子. Chapman和Orutt(1985)讨论了该方法在深部地球模型中的应用,并与FK方法的一维合成结果进行了比较. Chapman和Drummond(1982)将WKBJ方法扩展到二维,称为Maslov理论. Graves和Helmberger(1988)应用这种方法成功模拟了多个S波震相(S、SS等). 与广义射线理论不同,WKBJ法先进行频率反变换,然后在实路径上进行波数积分. WKBJ方法的优点在于只要采用已有的常微分方程近似解法(WKBJ法),使中间结果直观且容易解释. 然而,WKBJ地震图并不能较好地处理介质速度梯度问题,因为它无法模拟来自梯度带的低频反射,这是WKBJ近似所固有的不足. 因为在解表达式中其频率和空间是解耦的,解在深度z上的相对振幅仅取决于该深度的模型参数(Paulssen, 1988).

      从频率域二维标量波动方程开始简单介绍其原理(Chapman, 1978):

      $ {\nabla }^{2}\phi +{\omega }^{2}{v}^{-2}\left(z\right)\phi =0 $

      (14)

      为得到方程(14)的近似解,要求ω足够大. 试验解的形式如下:

      $ \phi (\omega,{{x}})={\phi }_{o}\left(\omega \right)\sum\nolimits_{n=0}^{\infty }\frac{{A}^{\left(n\right)}\left({{x}}\right)}{(-i\omega {)}^{n}}{\mathrm{e}}^{i\omega T\left({{x}}\right)} $

      (15)

      将展开式代入方程(14),并令$ {\omega }^{2} $项的系数之和为零,得到程函方程:

      $ (\nabla T{)}^{2}-{v}^{-2}(z)=0 $

      (16)

      ω任意次幂项之和为零,得到振幅系数$ {A}^{\left(n\right)}\left(x\right) $.

      在高频情况下,只有WKBJ渐近展开的前导项比较重要,从而得到WKBJ的近似解:

      $ \phi (\omega,{{x}})={\phi }_{o}\left(\omega \right){A}^{\left(0\right)}\left({{x}}\right){\mathrm{e}}^{i\omega T\left({{x}}\right)} $

      (17)

      这就是我们熟知的几何射线理论的解.

      对于两点边值问题,通常存在几种$ T\left(x\right) $$ {A}^{n}\left(x\right) $的解. 因此,在零阶WKBJ近似中可以写成:

      $ \phi (\omega,{{x}})={\phi }_{o}\left(\omega \right)\sum\nolimits_{j=1}^{\mathrm{rays}}{A}_{j}^{\left(0\right)}\left({{x}}\right){\mathrm{e}}^{i\omega {T}_{f}\left(x\right)} $

      (18)

      式中,“rays”代表解的数量.

      需要注意的是,只有在速度分布连续的情况下,零阶近似才有效. 或者,更一般地说,为了确定渐近展开中的第n项,模型参数的第n阶导数必须是连续的. 因此,如果遇到速度不连续,WKBJ近似将会失效. 这种情况下,需要利用Snell定律和适当的反射和透射系数得到界面处的解.

      上述几何射线理论在WKBJ地震记录中的推广涉及空间坐标到慢度域的额外傅里叶变换. WKBJ地震记录的关键是逆变换的顺序:在进行空间坐标变换之前,首先将频率慢度域中的解变换到时域. 对慢度的逆变换为(Dey-Sarkar and Chapman, 1978):

      $\begin{array}{l} \phi \left(t,{{x}}\right)\\ =\dfrac{-1}{{\text{π}}\sqrt{2}}{\phi }_{o}\left(t\right)*{\partial }_{t}{\rm{Im}}\left[\varLambda \left(t\right)*{\displaystyle\int }_{-\infty }^{\infty }{A}^{\left(0\right)}(p,z)\delta (t-\theta (p,{{x}}\left)\right)\mathrm{d}p\right] \end{array}$

      (19)

      式中:

      $ \varLambda \left(t\right)=\lambda \left(t\right)+i \bar{\lambda }\left(t\right) $

      (20)

      其中,$ \lambda \left(t\right)=H\left(t\right){t}^{-1/2} $$ \bar{\lambda }\left(t\right)=H(-t)(-t{)}^{-1/2} $.

      $ \theta (p,{{x}})=\tau (p,z)+px $

      (21)

      对于$ t=\theta (p,{{x}}) $,得到:

      $ \phi (t,{{x}})=\frac{-1}{{\text{π}}{\sqrt{2}}}{\phi }_{o}\left(t\right)*{\partial }_{t}{\rm{Im}}\left[\varLambda \left(t\right)*\sum\nolimits_{t=\theta }\frac{{A}^{\left(0\right)}(p,z)}{\left|{\partial }_{p}\theta \right|}\right] $

      (22)

      对于实$ t=\theta $,逆变换是精确的. 然而,对于某些慢度p$ \tau (p,z) $是复数,会产生指数衰减相位项. 在WKBJ地震图方法中,反变换只在实p等值线上进行.

      (3)简正模方法(normal modes)

      简正模方法是远震波场数值模拟中应用最广泛的方法之一. 1960年智利大地震后,地震学家对简正模进行了大量的研究,并逐渐发展出简正模态叠加(normal mode summation)(Florsch et al., 1991; Yang et al., 2010; Gualtieri et al., 2013; Liu and Zhang, 2017). Gilbert(1971)首先提出了用模态叠加的形式表示位移,并给出了相关的振型特征函数的显式表示. Singh和Ben-Menahem(1969a, 1969b)计算了球对称介质中点源激发的各个模式振型. Takeuchi和Saito(1972)推导了简正模的方程. Tanimoto(1984)得到了用模态叠加计算长周期合成地震图的公式. Woodhouse(1988)提出了一种计算径向特征函数的数值方法,这使得使用模态叠加计算合成地震记录成为可能. Dahlen和Tromp(1998)整合了前人的工作,构建了一个全球地震学理论框架,包括简正模理论、简正模态叠加方法等相关方法.

      所谓简正模态叠加法,就是将地球自由振荡的各个模式按一定系数叠加在一起,构成位移场(栾威等,2021). 地球存在两种自由振荡模式:球形振荡和环形振荡,分别对应于P-SV波系和SH波系. 简正模提供了一组完备的、线性无关的波动方程通解,因此任何特定的解,即由给定源激发的位移场,都可以通过简正波特征函数的线性叠加得到. 在球极坐标系中,震源和台站坐标分别由$ ({r}_{S},{\theta }_{S},{\phi }_{S}) $$ ({r}_{R},{\theta }_{R},{\phi }_{R}) $给出,其中$\mathbf {r}_\mathrm{R} $通常是地球的半径. 对于$ \mathbf {r}_\mathbf{S} $处的脉冲点源的地震,它可由力矩张量M表示,因此接收器$ \mathbf {r}_\mathbf{R} $处每个振型的位移场是力矩张量与震源处模态应变$ {\varepsilon }_{j} $的张量积和台站处模态位移$\mathbf {s}_{j} $的复共轭所给出的模态激励的乘积:

      $ \begin{array}{l} \mathbf{u}\left(\mathbf{r}_\mathrm{R},t;\mathbf{r}_\mathrm{S}\right)\\ =\mathrm{Re}\left[\displaystyle\sum\nolimits_{j=k,l,m,n}{v}_{j}^{-2}\left[\mathbf{M}:{\varepsilon} _{j}^{*}\left(\mathbf{r}_\mathrm{S}\right)\right]\mathbf {s}_{j}\left(\mathbf{r}_\mathrm{R}\right)\left[1-\mathrm{exp}(i{v}_{j}t)\right]\right] \end{array} $

      (23)

      式中,模态应变为:

      $ {\varepsilon }_{j}\left(\bf{r}\right)=\frac{1}{2}\left[\nabla \mathbf{s}_{j}+(\nabla \mathbf{s}_{j}{)}^{\mathrm{T}}\right] $

      (24)

      式中,上标T表示矩阵转置.

      对于具有滞弹性衰减的真实地球模型,简正模的本征频率$ {v}_{j} $是复数,实部和虚部$ {\omega }_{j} $$ {\sigma }_{j} $分别表征模式的时间振荡模式和衰减. 因此:

      $ {v}_{j}={\omega }_{j}+i{\sigma }_{j},{\sigma }_{j}=\frac{1}{2}{\omega }_{j}{Q}_{j}^{-1} $

      (25)

      式中,$ {Q}_{j} $是相应的模态品质因数.

      根据Rayleigh原理,弹性/滞弹性性质的微小扰动对简正模特征函数的影响相比其对特征频率的影响是二阶的,因此可以忽略. 所以若考虑滞弹性,由于忽略了特征函数的影响,在引入与频率相关的体积模量和剪切模量时,会导致所谓的物理频散(Dahlen and Tromp, 1998). 此外,假设地球略有滞弹性,也就是说,$ {\sigma }_{j}<<{\omega }_{j} $,方程(23)变为:

      $ \begin{array}{l} \mathbf{u}(\mathbf{r}_\mathrm{R},t;\mathbf{r}_\mathrm{S})\\ =\displaystyle\sum\nolimits_{j=k,l,m,n}\left[\mathbf{M}:{\varepsilon }_{j}^{*}\left(\mathbf{r}_\mathrm{S}\right)\right]\mathbf{s}_{j}\left(\mathbf{r}_\mathrm{R}\right)\frac{1-\mathrm{exp}(-{\sigma }_{j}t)\mathrm{cos}{\omega }_{j}t}{{\omega }_{j}^{2}} \end{array} $

      (26)

      这是用简正模态叠加表示位移的常见形式.

      简振模方法仍然在全球和区域尺度的面波层析成像中广泛使用(Chevrot and Zhao, 2007; Ritsema et al., 2011),主要优点是以很小的计算成本提供精确的合成地震图,并能包括滞弹性、地球自转和重力的影响;主要缺点是在8 s以下的周期内计算球形振荡很困难(Al-Attar and Woodhouse, 2008),并且要求简振模求和的模数随着频率的增加而急剧增加. 因此,简正模态叠加方法非常适合于计算长周期的地震记录,不适合模拟短周期远震体波.

    • 对于横向均匀的一维地球模型,可以通过半解析法将波动方程分解为水平和垂直分量,水平分量是偏微分方程,其解析解为平面分层模型的矢面调和或球对称模型的矢面球谐. 垂直分量为常微分方程组(ODE)的数值解,通过求解该常微分方程,从而达到精度和效率的折衷(Zhang et al., 2003; Yang et al., 2010).

      目前已发展出各种方法来求解常微分方程组,例如,离散波数法(DWN)、WKM法、频率波数域法(FK)等. 其中WKM法是WKBJ方法的改良版. 这些方法给出了球对称地球模型波动方程的精确解,即忽略了地球结构的横向变化,却充分而精确地考虑了地球球状的曲率效应和有限频率地震波的散射、波前愈合及射线路径在低速区发散、在高速区聚集而引起射线分布非均匀等复杂的波现象问题(Dahlen and Tromp, 1998; 张明辉等,2019).

      (1)离散波数法(DWN)

      离散波数法由Bouchon和Aki(1977)提出,由于其能够精确地求解出完备的格林函数,在许多弹性动力学计算问题中得到了广泛的应用(Bouchon and Aki, 1977; Bouchon, 1979, 1981, 2003; Chen, 1990; Chen and Zhang, 2001; Fu and Bouchon, 2004; Zhou and Chen, 2006, 2008; Wu et al., 2020). 该方法的基本思想是将震源表示为以离散角度传播的均匀和非均匀平面波的叠加(Bouchon and Aki, 1977). 离散波数法的原理可以追溯到Rayleigh,他证明了由正弦波纹界面反射的波只能以离散的角度传播,被称为频谱阶数(Rayleigh and Baron, 1896; Rayleigh, 1907). 水平波数谱中离散阶的存在直接导致了反射界面的周期性. Aki和Larner(1970)利用复频率扩展了Rayleigh方法来研究周期性不规则表面附近的平面波散射. 同样,离散波数方法引入震源的空间周期性来离散辐射波场,并依靠复频域中的傅里叶变换来计算格林函数.

      离散波数法的基本原理如下:无限均匀介质中线源的稳态辐射可以表示为柱面波,或者等价地表示为均匀和非均匀平面波的连续叠加. 因此,用xz表示垂直于线源的平面中的水平轴和垂直轴,任何可观测数据(如位移或应力)都可以写成以下形式(Bouchon, 2003):

      $ F(x,z;\omega)={\mathrm{e}}^{i\omega t}{\int }_{-\infty }^{\infty }f(k,z){\mathrm{e}}^{-ikx}{\rm{d}}k $

      (27)

      当介质有限或垂直非均匀时,积分核具有极点和奇点,水平波数上的积分在数学上和数值上都变得复杂起来. 绕过这些困难的一种简单方法是将其解由公式(27)表示的单源问题替换为沿x轴周期性分布震源的多源问题(图2). 然后,公式(27)被替换为:

      图  2  DWN方法的物理解释. 将单个源替换为以相等间隔L水平分布的无限多个源阵列. 对于与特定激励频率相对应的给定辐射波长λ,弹性能量仅在离散方向θ上辐射(修改自Bouchon, 2003).

      Figure 2.  Physical interpretation of the DWN method. The single source is replaced by an infinite array of sources distributed horizontally at equal interval L. For a given radiation wavelength k corresponding to a specific frequency of excitation, the elastic energy is radiated in discrete directions $ \theta $ only(modified from Bouchon, 2003

      $ G(x,z;\omega)={\int }_{-\infty }^{\infty }f(k,z){\mathrm{e}}^{-ikx}\sum\nolimits_{m=-\infty }^{\infty }{\mathrm{e}}^{ikmL}{\rm{d}}k $

      (28)

      式中,L是周期性震源间隔,方程(28)简化为:

      $ G(x,z;\omega)=\frac{2{\text{π} }}{L}\sum\nolimits_{n=-\infty }^{\infty }f({k}_{n},z){\mathrm{e}}^{-{i}{k}_{n}x} $

      (29)

      式中,$ {k}_{n}=\dfrac{2{\text{π}}}{L}n $,如果级数收敛,可用有限求和方程来近似:

      $ G(x,z;\omega)=\frac{2{\text{π} }}{L}\sum\nolimits_{n=-N}^{N}f({k}_{n},z){\mathrm{e}}^{-i{k}_{n}x} $

      (30)

      通过上面的假设和变换,计算量大幅减小. 已经将单源问题变成了含有周期性震源分布的多源问题(图2). DWN方法需要计算方程(30),也就是$ G(x,z;\omega) $.

      第二步,要从频域中求解的多源问题中提取单源解. 可以通过在复频域中进行傅里叶变换来实现:

      $ g(x,z;t)={\int }_{-\infty +i{\omega }_{I}}^{\infty +i{\omega }_{I}}G(x,z;\omega){\mathrm{e}}^{i\omega t}{\rm{d}}\omega $

      (31)

      式中,$ {\omega }_{I} $表示频率的常量虚部,并且:

      $ {\mathrm{e}}^{{\omega }_{I}T}<<1 $

      (32)

      方程(32)确保了先前无限时间响应解在时间窗口T上的衰减. 因此,假设我们选择了足够大的L,使得最接近的震源到达接收点(x, z)在感兴趣时间窗口T中不受干扰,则可以从多源频率域解计算$ G(x,z;\omega) $得到时域单源解$ f(x,z;t) $

      $ f(x,z;t)={\mathrm{e}}^{-{\omega }_{I}t}{\int }_{-\infty }^{\infty }G(x,z;\omega){\mathrm{e}}^{i{\omega }_{R}t}{\rm{d}}{\omega }_{R} $

      (33)

      可以通过FFT计算积分.

      (2)WKM法

      为了得到更高精度的波形来反映地幔深部结构的非均匀性,Ni等(2000)发展了一种新方法,该方法能直接从块状的层析模型中计算二维合成地震图. 该方法是WKBJ方法的改良版,因此叫做WKM方法. 它跟WKBJ近似不同,WKBJ近似是利用越过或未到达台站的射线,而WKM近似则是和GRT方法类似,利用到达台站的射线. WKM利用一维层状参考模型的射线路径被用来定位每个射线段,通过不断调整$ {p}_{i}\left({t}_{i}\right) $$ {p}_{i} $为射线参数,$ {t}_{i} $为走时)以满足Snell定律及其数值导数$ (\partial p/\partial t) $. Ni等(2000)基于Grand的层析成像模型(Grand et al., 1994)生成了中美洲下方高速$ {D}^{{'}{'}} $区域的WKM合成地震图. Ni等(2003)基于与南非热柱相关的S速度异常的二维模型(图3),分别用WKM与伪谱法合成地震图,并进行比较(图4). 结果显示,两种方法合成的SV波场和SH波场基本吻合(Ni et al., 2003).

      图  3  层析成像获得的从南美到南非的二维速度剖面(修改自Ritsema et al., 1999; Ni et al., 2003

      Figure 3.  2D velocity section from South America to South Africa obtained from tomography (modified from Ritsema et al., 1999; Ni et al., 2003)

      图  4  WKM法(左)和伪谱法(右)生成的SV波合成地震图对比(修改自Ni et al., 2003

      Figure 4.  Comparison of SV synthetics generated by the WKM method against the pseudo-spectral method (on the right)(modified from Ni et al., 2003)

      WKM近似对于合成二维理论地震图非常有效. 但是WKM不能用于地球内部结构有特别尖锐变化的地方,这是因为尖锐介质处速度差异都很大,这样产生的衍射作用很强,而WKM近似不能有效考虑这种衍射作用(Ni et al., 2000).

      (3)频率波数域法(FK)

      频率波数域法(FK)是计算远震波场最重要的方法之一(Thomson, 1950; Haskell, 1953, 1964; Takeuchi and Saito, 1972; Schmidt and Tango, 1986; Zhu and Rivera, 2002),该方法在频率和波数域采用传播矩阵计算地震的位移场分布. 由于该方法理论清晰,计算准确,因此在计算分层模型的理论合成地震图时,应用非常广泛(Tong et al., 2014a).

      已知在各向同性介质中,频率域的弹性波方程可以写成:

      $ -\rho {\omega }^{2}{{u}}=\nabla \cdot \left[\lambda (\nabla \cdot {{u}})+\mu \left(\nabla {{u}}+\nabla {{{u}}}^{\mathrm{T}}\right)\right] $

      (34)

      式中,u代表位移矢量,ρλμ表示介质中密度和拉梅参数的空间分布. 在水平层状介质中,位移场u仅依赖于水平和垂直坐标,即u = u(x, z, ω). 对u做傅里叶变换,将空间域x变换为水平波数域k,将时间域t变换为频率域ω,可得:

      $ {{u}}=\left[-i{y}_{1}(z,k;\omega),{y}_{2}(z,k;\omega),{y}_{3}(z,k;\omega)\right] $

      (35)

      于是方程(34)可以简化成两组一阶常微分方程(Takeuchi and Saito, 1972; Zhu and Rivera, 2002; Liu et al., 2017a). 通过计算特征值和特征向量,得到两个常微分方程组的通解:

      $ \begin{array}{l} \left[\begin{array}{c}{y}_{1}\\ {y}_{3}\\ {y}_{4}\\ {y}_{6}\end{array}\right]=\left[\begin{array}{cccc}-i\dfrac{{v}_{{S}}}{k}& i\dfrac{{v}_{S}}{k}& 1& 1\\ 1& 1& -i\dfrac{{v}_{{P}}}{k}& i\dfrac{{v}_{{P}}}{k}\\ 2k\mu {\gamma }_{1}& 2k\mu {\gamma }_{1}& -2i\mu {v}_{{P}}& 2i\mu {v}_{{P}}\\ -2i\mu {v}_{{S}}& 2i\mu {v}_{{S}}& 2k\mu {\gamma }_{1}& 2k\mu {\gamma }_{1}\end{array}\right]\cdot \\ \left[\begin{array}{cccc}{\mathrm{e}}^{-i{v}_{{S}}z}& 0& 0& 0\\ 0& {\mathrm{e}}^{i{v}_{S}z}& 0& 0\\ 0& 0& {\mathrm{e}}^{-i{v}_{{P}}z}& 0\\ 0& 0& 0& {\mathrm{e}}^{i{v}_{P}z}\end{array}\right]\left[\begin{array}{c}{C}_{1}\\ {C}_{3}\\ {C}_{4}\\ {C}_{6}\end{array}\right] \end{array} $

      (36)

      $ \left[\begin{array}{c}{y}_{2}\\ {y}_{5}\end{array}\right]=\left[\begin{array}{cc}1& -1\\ -i\mu {v}_{{S}}& -i\mu {v}_{{S}}\end{array}\right]\left[\begin{array}{cc}{\mathrm{e}}^{-i{v}_{{S}}z}& 0\\ 0& {\mathrm{e}}^{i{v}_{{S}}z}\end{array}\right]\left[\begin{array}{c}{C}_{2}\\ {C}_{5}\end{array}\right] $

      (37)

      式中,$ {C}_{1}{\text{、}}{C}_{3}{\text{、}}{C}_{4} $$ {C}_{6} $分别对应上行SV波、下行SV波、上行P波和下行P波的振幅;而$ {C}_{2}{\text{、}}{C}_{5} $对应的是上行SH波、下行SH波的振幅. $ p=\dfrac{k}{\omega } $为射线参数,αβ分别表示该层介质中P波和S波的波速,则:

      $ {v}_{{P}}=\omega \sqrt{\frac{1}{{\alpha }^{2}}-{p}^{2}},{v}_{{S}}=\omega \sqrt{\frac{1}{{\beta }^{2}}-{p}^{2}},{\gamma }_{1}=1-\frac{1}{2{p}^{2}{\beta }^{2}} $

      (38)

      $ {v}_{{P}} $$ {v}_{{S}} $分别为P波和S波的垂直波数,$ {\gamma }_{1} $为辅助变量.

      当平面波入射N层层状半空间介质时(图5),对于第m层,上层$ z={z}_{m}{\text{与下层}}z={z}_{m-1} $的变量y可通过传播矩阵$ \mathbf{P}_{m} $联系起来,即:

      图  5  底部半空间上的N层组成的分层半空间

      Figure 5.  A layered half-space consists of N layers over a half-space at the bottom.

      $ {\mathbf{y}}_{m}={\mathbf{P}}_{m}{\mathbf{y}}_{m-1} $

      (39)

      对于P波平面波入射情况,我们将底层设为$ z=0 $,并假设从底部第1层到顶部第n层的传播矩阵分别为$ \mathbf{P}_{1},\mathbf{P}_{2},...,\mathbf{P}_{n} $. 由方程(36),可得到在$ z=0 $处的波场:

      $ \begin{array}{l} {\left[\begin{array}{c}{y}_{1}\\ {y}_{3}\\ {y}_{4}\\ {y}_{6}\end{array}\right]}_{z=0}={\left[\begin{array}{cccc}-i\dfrac{{v}_{{S}}}{k}& i\dfrac{{v}_{{S}}}{k}& 1& 1\\ 1& 1& -i\dfrac{{v}_{{P}}}{k}& i\dfrac{{v}_{{P}}}{k}\\ 2k\mu {\gamma }_{1}& 2k\mu {\gamma }_{1}& -2i\mu {v}_{{P}}& 2i\mu {v}_{{P}}\\ -2i\mu {v}_{{S}}& 2i\mu {v}_{{S}}& 2k\mu {\gamma }_{1}& 2k\mu {\gamma }_{1}\end{array}\right]}_{0} \cdot \\ \left[\begin{array}{c}{C}_{1}\\ {C}_{3}\\ {C}_{4}\\ {C}_{6}\end{array}\right]={\mathbf{E}}_{0}(\mu,\xi ;\omega,p)\left[\begin{array}{c}{C}_{1}\\ {C}_{3}\\ {C}_{4}\\ {C}_{6}\end{array}\right] \end{array} $

      (40)

      地表为:

      $ \left[\begin{array}{c}{y}_{1}\\ {y}_{3}\\ {y}_{4}\\ {y}_{6}\end{array}\right]={\mathbf{P}}_{n}\cdots {\mathbf{P}}_{1}{\mathbf{E}}_{0}\left[\begin{array}{c}{C}_{1}\\ {C}_{3}\\ {C}_{4}\\ {C}_{6}\end{array}\right]={\mathbf{N}}\left[\begin{array}{c}{C}_{1}\\ {C}_{3}\\ {C}_{4}\\ {C}_{6}\end{array}\right] $

      (41)

      式中:

      $ {{N}}={\left[{N}_{ij}\right]}_{4\times 4}={\mathbf{P}}_{n}\cdots {\mathbf{P}}_{1}{\mathbf{E}}_{0} $

      (42)

      上式中$ {\mathbf{E}}_{0} $是由半空间($ z=0 $)模型参数计算得到的矩阵. 而$ {\mathbf{P}}_{i} $中的下角标$ i $表示在第$ i $层介质中的传输矩阵.

    • (1)直接法(DSM)

      直接法是一种Galerkin方法,它在频率域求解波动方程的弱形式(Cummins et al., 1994; Geller and Ohminato, 1994; Geller and Takeuchi, 1995; Takeuchi et al., 1996). 之所以称为直接解法,是因为该解是通过直接求解线性方程组得到的,而非首先计算自由振荡模式,然后通过模态叠加计算合成地震图(Geller and Ohminato, 1994).

      具有自由表面边界条件介质的DSM运动方程为(Takeuchi and Saito, 1972; Kawai et al., 2006):

      $ \left({\omega }^{2}{{T}}-{{H}}+{{R}}\right){{c}}=-{{g}} $

      (43)

      式中,$ \omega $是频率,T是质量(动能)矩阵,H是刚度(势能)矩阵,R是对应于自由边界条件的矩阵算子. c是测试函数的展开系数向量,g是力矢量. 矩阵和向量元素的显式形式为:

      $ \!\!\!\!\!\! {T}_{mn}={\int }_{V}{\left({\phi }_{i}^{\left(m\right)}\right)}^{*}\rho {\phi }_{i}^{\left(n\right)}{\rm{d}}V $

      (44)

      $ \;\; {H}_{mn}={\int }_{V}{\left({\phi }_{i,j}^{\left(m\right)}\right)}^{*}{C}_{ijkl}{\phi }_{k,l}^{\left(n\right)}{\rm{d}}V $

      (45)

      $ {R}_{mn}={\int }_{S}{\left({\phi }_{i}^{\left(m\right)}\right)}^{*}{S}_{ij}{\phi }_{j}^{\left(n\right)}{\rm{d}}S $

      (46)

      $\!\!\!\!\!\!\!\!\!\!\!\!\! {g}_{m}={\int }_{V}{\left({\phi }_{i}^{\left(m\right)}\right)}^{*}{f}_{i}{\rm{d}}V $

      (47)

      位移矢量由N个测试向量测试函数$ {\phi }_{i}\left(n\right) $的线性组合给出:

      $ {u}_{i}\left({{r}}\right)=\sum\nolimits_{n=1}^{N}{c}_{n}{\phi }_{i}\left(n\right)\left({{r}}\right) $

      (48)

      用线性插值函数$ {X}_{k}\left(r\right) $来描述位移的垂直变化:

      $ {X}_{k}\left(r\right)=\left\{\begin{array}{c}(r-{r}_{k-1})/({r}_{k}-{r}_{k-1}), {r}_{k-1}<r\le {r}_{k}\\ ({r}_{k+1}-r)/({r}_{k+1}-{r}_{k}), {r}_{k}\le r<{r}_{k+1}\\ \quad\quad \text{0} \quad\quad\quad\qquad {\text{其他}}\end{array}\right. $

      (49)

      用复矢量球谐$ {{S}}_{lm}^{1}{\text{、}}{{S}}_{lm}^{2}{\text{、}}{{{T}}}_{lm} $来描述位移的横向变化:

      $ \begin{array}{c} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! {{S}}_{lm}^{1}(\theta,\phi)=\left({Y}_{lm}(\theta,\phi),\mathrm{0,0}\right)\\ {{S}}_{lm}^{2}(\theta,\phi)=\left(0,\dfrac{1}{\mathcal{L}}\dfrac{\partial {Y}_{lm}(\theta,\phi)}{\partial \theta },\dfrac{1}{\mathcal{L}\sin \theta }\dfrac{\partial {Y}_{lm}(\theta,\phi)}{\partial \phi }\right)\\ \;\;\, {{{T}}}_{lm}(\theta,\phi)=\left(0,\dfrac{1}{\mathcal{L}\sin \theta }\dfrac{\partial {Y}_{lm}(\theta,\phi)}{\partial \phi },-\dfrac{1}{\mathcal{L}}\dfrac{\partial {Y}_{lm}(\theta,\phi)}{\partial \theta }\right)\end{array} $

      (50)

      式中,$ {Y}_{lm}(\theta,\phi) $是完全归一化的球谐函数.

      在矢量球谐基函数中,位移矢量由下式给出:

      $ {{u}}(r,\theta,\phi)=\sum\nolimits_{lm}{U}_{lm}\left(r\right){{S}}_{lm}^{1}+{V}_{lm}\left(r\right){{S}}_{lm}^{2}+{W}_{lm}\left(r\right){{{T}}}_{lm} $

      (51)

      式中,径向函数与径向测试函数的关系:

      $ {U}_{lm}\left(r\right)=\sum\nolimits_{k}{c}_{lmk}^{1}{X}_{k}\left(r\right) $

      (52)

      $ {V}_{lm}\left(r\right)=\sum\nolimits_{k}{c}_{lmk}^{2}{X}_{k}\left(r\right) $

      (53)

      $ {W}_{lm}\left(r\right)=\sum\nolimits_{k}{c}_{lmk}^{3}{X}_{k}\left(r\right) $

      (54)

      一旦得到位移,就可以计算位移相对于r$ \theta $$ \phi $的偏导数. 测试函数$ {{S}}_{lm}^{1}{\text{、}}{{S}}_{lm}^{2}{\text{、}}{{{T}}}_{lm} $关于$ \phi $$ \theta $的导数可以得到解析表达式:

      $ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \frac{\partial {Y}_{lm}}{\partial \phi }=im{Y}_{lm} $

      (55)

      $ \frac{\mathrm{d}{P}_{lm}}{\mathrm{d}\theta }=m\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }{P}_{lm}-{P}_{lm-1} $

      (56)

      在DSM中,位移在频域通过一系列基函数展开,这些基函数由垂直方向的低阶多项式和角方向的球谐函数的乘积构成. 然后用伽辽金法求解这些基函数的系数. 由于在垂直方向上进行离散化处理,点源通常不会位于垂直节点之间. 在这些情况下,Takeuchi和Geller(2003)提出的震源表示方法可用于合成计算,其精度与震源位于节点处的解相同. 通过调整垂直网格间距、最大角度阶数和截止深度,即使频率在高达2 Hz的情况下,DSM可以非常有效和准确地计算理论地震图(Kawai et al., 2006). 因此,该方法特别适用于短周期远震波场数值模拟.

      (2)AxiSEM

      轴对称谱元法(Axisymmetric Spectral Element Method, AxiSEM)最初由Nissen-Meyer于2007年提出(Nissen-Meyer et al., 2007a, 2007b),是基于球对称的一维地球模型(例如PREM、IASP91)解决三维全球声波和弹性波传播的并行谱元方法,该方法的创新之处在于考虑了球面地球模型的轴对称性,使地震速度和地球密度随地球极轴的旋转而保持不变(Monteiller et al., 2020). 通过考虑对称性,首先对二维大尺度背景下的半圆区域进行常规谱元法数值模拟,然后经过方位角旋转实现降维处理,获取第三个维度任意方位角台站处的波场信息. 通过将震源矩张量$ {{M}} $的6个分量(其响应从震源出发到接收器,为正向波场)以及地震台站的三分量点源$ \widehat{p} $(其响应为从接收器到震源,为反向波场)的全波场响应按照单极源、偶极源、四极源模式分解为独立的二维问题(图6),从而实现三维问题向二维问题的降维. 作为全球地球模型地震图计算方法的一种,轴对称谱元法能够准确地计算介质的质点速度以及应变响应,兼顾了四边形网格谱元法的计算精度,能够实现全球范围内的高频地震波数值模拟.

      图  6  震源分解图样(从左到右为单极源、偶极源、四极源)(修改自Nissen-Meyer et al., 2014

      Figure 6.  Radiation patterns for monopole, dipole, and quadrupole from left to right (modified from Nissen-Meyer et al., 2014)

      为了描述矩张量源的辐射方向图,需要进行多个二维模拟通过线性组合以重建完整的三维波场. 根据矩张量震源产生的地震波场方位依赖性,可以将其表示为具有单极(m=0)、偶极(m=1)和四极(m=2)源项的有限个多极展开. 然后,与给定极阶数m相关的波场可以用柱坐标表示为:

      $ {{{u}}}_{m}\left({{x}}\right)=\left(\begin{array}{c}{u}_{s}\left(\widetilde {{x}}\right)\mathrm{cos}m\varPhi \\ {u}_{\varPhi }\left(\widetilde {{x}}\right)\mathrm{sin}m\varPhi \\ {u}_{z}\left(\widetilde {{x}}\right)\mathrm{cos}m\varPhi \end{array}\right) $

      (57)

      式中,$ {\mathbf{x}}=(s,\varPhi,z) $分别是沿径向轴、方位向轴和对称轴的坐标,$ \widetilde {\mathbf{x}}=(s,z) $与前者相同,用于二维AxiSEM计算区域.

      地震矩张量的全部响应最终表示为对应于四个基本震源(两个单极、一个偶极和一个四极)的格林函数的线性组合:

      $ {\mathbf{u}}({\mathbf{x}},t)=\sum\nolimits_{k=1}^{4}{\widetilde {M}}_{k}\left(\varPhi \right){\widetilde {\mathbf{G}}}_{k}\left(\widetilde {\mathbf{x}},t;{z}_{s}\right)S\left(t\right) $

      (58)

      式中,S(t)为震源时间函数,$ {\widetilde {M}}_{k} $是基本源项的辐射方向图,$ {\widetilde {\mathbf{G}}}_{k} $是与基本源k对应的二维格林函数(位于深度为$ {z}_{s} $的对称轴上). 这四个辐射方向图本身是矩张量解的分量的线性组合:

      $ \begin{array}{c}{\widetilde {M}}_{1}\left(\varPhi \right)={M}_{rr}\\ {\widetilde {M}}_{2}\left(\varPhi \right)=\left({M}_{\theta \theta }+{M}_{\varPhi \varPhi }\right)/2\\ {\widetilde {M}}_{3}\left(\varPhi \right)={M}_{r\theta }\cos \varPhi +{M}_{r\varPhi }\sin\varPhi \\ {\widetilde {M}}_{4}\left(\varPhi \right)=\left({M}_{\theta \theta }-{M}_{\varPhi \varPhi }\right)\mathrm{}\cos 2\varPhi +{M}_{\varPhi \theta }\sin 2{\varPhi }\end{array} $

      (59)

      式中,$ {\widetilde {M}}_{1} $$ {\widetilde {M}}_{2} $是两个单极源,$ {\widetilde {M}}_{3} $是偶极源,$ {\widetilde {M}}_{4} $是四极源.

      一旦计算得到4个基本波场,并且用方程(59)重建地震波场,并被投影回球面坐标系. 最后,将震源绕对称轴旋转,可以得到在地理参考系中任意位置的波场.

    • 尽管基于一维地球模型的快速算法具有计算频率高、运行速度快的优点,但难以处理具有横向非均匀的三维区域模型. 而具有高精度的三维数值方法由于受到硬件计算能力和存储能力的限制,目前仍然无法高效地模拟高频全球尺度的地震波场(例如全球尺度周期为1~2 s的P波,周期3~6 s的S波)(Liu and Gu 2012; Monteiller et al., 2013). 这也限制了利用高频地震资料进行区域乃至全球尺度的高分辨率成像研究. 为了降低远震波正演模拟的计算成本,很多学者已经开发出注入方法来计算区域三维模型中的理论地震图(Bielak and Christiano, 1984; Roecker et al., 2010; Monteiller et al., 2013; Tong et al., 2014a, 2014b; Masson and Romanowicz, 2017a, 2017b; Beller et al., 2018a; Lin et al., 2019). 在勘探地震学(Robertsson and Chapman, 2000; Ivo et al., 2002; Van Manen et al., 2007)和电磁波传播(Taflove and Hagness, 2005)等领域也有类似的技术. 该方法的关键思想是,目标区域通常远离震源区域(也就是远震震源的情况下),因此,波场传播问题通常简化考虑为两部分:全球尺度地震波传播和区域尺度地震波传播问题. 对于全球范围内的地震波场计算,可以采用前面介绍的基于一维地球模型的快速算法. 较为耗时的三维非均匀区域的地震波场计算被限制在较小的区域范围内,并且可以使用合理的计算资源来有效地计算短周期的三维波场. 一般采用三维数值方法来计算目标区域,如有限差分法(FD)(Alterman and Karal, 1968)和谱元法(SEM)(Komatitsch and Tromp, 1999; Liu et al., 2017b)等,能够处理影响地震波传播的多种复杂因素,例如自由表面和内部不连续性的界面起伏,各向同性和各向异性的非均匀介质以及衰减问题等.

    • 早在20世纪80年代,学者们就提出了利用有限差分法与其他各种方法结合的地震波场数值模拟方法. Shtivelman(1984)利用有限差分技术实现声波在结构非均匀部分的传播. 用改进的射线法计算均匀介质中的背景波场(Shtivelman, 1984, 1985),以及将有限差分和边界积分法相结合(Kummer et al., 1987; Stead and Helmberger, 1988). 温联星和姚振兴(1994)提出广义射线和有限差分理论相结合的混合方法,很好地应用了有限差分方法和广义射线法的优点,既能计算非均匀区域的地震波传播,也提高了计算效率. 但基于有限差分法的混合法也存在缺点,如由于离散方法中固有的数值频散经过一段传播距离之后会使波形发生畸变,且难以考虑地球曲率问题,因此它不适合计算大震中距的地震波场(温联星和姚振兴, 1994). Wen和Helmberger(1998a)以及Wen(2002)提出了有效的P-SV和SH混合方法,利用该方法解决地幔边界附近的局部非均匀各向同性结构的波场传播问题. 该混合方法已被用于研究核幔边界附近的各向异性结构(Wen and Helmberger, 1998a, 1998b; Wen and Niu, 2002; Wen, 2002)以及岩石圈和上地幔结构(Chen et al., 2005). Zhao等(2008)提出了用GRT-FD混合法计算了二维各向异性非均匀介质的理论波形.

      混合方法的概念如图7a所示(Zhao et al., 2008). 计算包括广义射线理论(GRT)、GRT-FD界面耦合和FD计算三部分. 假设方形中的局部区域是横向变化的各向异性介质,在局部区域采用FD算法. GRT(Helmberger, 1983)用于计算从震源到FD区域底部的地震波的传播. 将GRT的解与(图7b)中阴影区域的FD计算相结合. 在FD计算的顶部网格点上得到地球表面的合成地震图,在这些计算中采用了地球展平近似(Zhao et al., 2008).

      图  7  GRT-FD混合方法原理示意图(修改自Zhao et al., 2008

      Figure 7.  Schematic diagram of the principle of hybrid method (modified from Zhao et al., 2008)

    • (1)DSM-SEM

      Capdeville等(2003a)提出耦合法,通过区域分解,沿耦合界面匹配位移和牵引力,将谱元法与简正叠加方法连接起来,并将其扩展到“夹层”域(sandwich domain)的情况(Capdeville et al., 2003a, 2003b). 但是,该方法要求目标区域必须位于接收器所在的地球表面的正下方,因此只能用于接收阵列下方的精细结构正反演(如接收函数成像等).

      Monteiller等(2013)结合直接解法(DSM)与谱元法(SEM)模拟远震波形数据,在全球一维球状分层地球模型下采用DSM方法进行全球高频理论地震图计算,存储目标区域边界处的波场;利用注入方法,在预设好的局部目标区域利用三维谱元法求解该区域的非均匀介质波场传播;在区域分界面上采用吸收边界条件耦合两种正演方法,形成一个高效、高精度三维混合全波场计算方法. 该混合方法将计算量较大的三维波场数值模拟限制在较小的局部目标区域内,极大降低了三维非均匀介质中地震波场正演计算的消耗,尤其适用于区域阵列成像应用,使得利用高频远震波形进行全波形反演成为可能. Monteiller等(2015)随后用该混合方法实现了三维全波形反演,并展示了全波形反演较好的分辨率,即使仅在四个远震事件情况下,可以成像具有尖锐莫霍面的简单地壳模型(图89).

      图  8  用于全波形反演的P(上)和S(下)波速度模型(修改自Monteiller et al., 2015

      Figure 8.  Synthetic model of P (top) and S (bottom) velocities for full waveform inversion (modified from Monteiller et al., 2015)

      图  9  由垂直分量直达P波的(a)全波形反演和(b)伴随层析成像得到的模型,在伴随层析成像情况下L-BFG迭代15次,在全波形反演情况下迭代5次(修改自Monteiller et al., 2015

      Figure 9.  Models obtained by (a) full waveform inversion and (b) adjoint tomography from vertical-component direct P waves, after 15 iterations of the L-BFGS in the case of adjoint tomography and 5 iterations in the case of full waveform inversion (modified from Monteiller et al., 2015)

      Masson等(2014)综述了基于注入技术混合法的理论,并将其命名为“箱式”技术,随后Masson和Romanowicz(2017a, 2017b)将混合法应用到波形反演中并将其称为箱式层析成像(box tomography),实现了深部地球成像问题. Wu等(2018)将谱元法模拟限制在一个较小的震源侧区域,有效地利用DSM-SEM混合方法计算具有3-D震源侧结构的远震波. Wang等(2016)利用29个台站记录到的5个远震波形记录,基于混合法开展了远震P波波形反演研究,得到比利牛斯山西部下方的P波和S波速度模型,结果显示波形反演成像得到的速度结构与接收函数偏移图像之间具有很好的一致性.

      这些工作对远震体波全波形反演(FWI)的发展具有重要意义,但是DSM涉及存储和处理大量的球谐系数,这可能会造成严重的I/O瓶颈问题. 对于浅源来说这个问题尤为严重,因为浅源地震的DSM模拟计算需要非常大的角阶数(angular orders)(Monteiller et al., 2020).

      (2)SEM-FK

      Tong等(2014a)提出了三维SEM-FK混合方法,将远震入射波看成平面波,利用计算消耗更小的FK算子作为谱元法的边界输入,并通过理论合成数据验证了方法的可靠性. 通过在区域场外使用FK算子计算一维介质的平面波响应,在区域场内使用谱元法模拟远震体波在三维复杂介质中的传播,在保持高计算效率的同时还能精确地捕捉到线性台站下方非均匀结构体对远震体波的影响,同时混合方法能产生和莫霍面相关的各种震相(及尾波),对速度间断面也有很好地约束(图10). 然而,FK方法忽略了波前和地球的曲率,并且无绝对到时等,在接收函数正演模拟方面有一定的应用前景,但是在远震波形层析成像研究中存在一定的困难和不足.

      图  10  (a)在西藏中部部署的Hi-CLIMB台站(三角形)的地理分布. 红色三角形表示用于观测和合成RF剖面的台站. $ \mathrm{B}{\mathrm{B}}^{'} $Tseng等(2009)估算地壳厚度的剖面. (b)Hung等(2011)绘制的沿剖面$ \mathrm{B}{\mathrm{B}}^{'} $$ {V}_{\mathrm{P}} $和(c)$ {V}_{\mathrm{S}} $扰动. (d)Tseng等(2009)根据$ {T}_{\mathrm{SsPmp}}-{T}_{\mathrm{Ss}} $估计的地壳厚度(蓝点). (e)2005年发生在北纬5.32°、东经123.34°、深度522 km的地震的SV波地震反射剖面(垂直分量速度). 各种散射/转换的相位,用紫色箭头表示. (f)基于局部地壳模型的三维SEM-FK混合方法计算相应的SV波合成地震剖面,该模型综合了CRUST1.0、估计的莫霍剖面、三维$ {V}_{\mathrm{P}} $$ {V}_{\mathrm{S}} $变化. 图4e图4f中的所有地震记录都在Ss震相上对齐(修改自Tong et al., 2014b

      Figure 10.  (a) Geographic distributions of Hi-CLIMB stations (triangles) previously deployed in central Tibet. Red triangles indicate stations used to generate the observed and synthetic RF profiles. $ \mathrm{B}{\mathrm{B}}^{'} $ is the profile along which the crust thickness was estimated by Tseng et al. (2009). (b)VP and (c) VS perturbations along Profile $ \mathrm{B}{\mathrm{B}}^{'} $ mapped by Hung et al. (2011) superimposed onto the estimated Moho. (d) Estimated crust thickness (blue dots) based on $ {T}_{\mathrm{SsPmp}}-{T}_{\mathrm{Ss}} $ by Tseng et al. (2009). (e) Observed SV wave seismic reflection profile (vertical-component velocity) for a 2005 earthquake occurred at $ 5.3{2}^{\circ } $ N, $ 123.3{4}^{\circ } $ E, and a depth of 522 km. Various scattered/converted phases, including SsPmp, are indicated by purple arrows. (f) Corresponding synthetic SV wave seismic profile computed by the 3-D SEM-FK hybrid method based on a local crustal model that incorporates CRUST1.0, the estimated Moho profile, and the 3-D VP and VS variations. All seismograms in Figures 4e and 4f are aligned on the Ss phase (modified from Tong et al., 2014b)

      (3)AxiSEM-SEM

      除了DSM-SEM和FK-SEM,AxiSEM-SEM同样被证明是一种有效的模拟远震波场的混合数值模拟方法(Beller et al., 2018a, 2018b; Lin et al., 2019; Monteiller et al., 2020). 该方法利用AxiSEM计算地震波在一维地球模型中传播的入射波场(Nissen-Meyer et al., 2014),用三维谱元法模拟地震波在岩石圈目标区域的传播. Beller等(2018b)使用该方法研究了西南阿尔卑斯山脉的岩石圈结构,通过对CIF ALPS实验记录的9个远震事件进行3-D各向同性全波形反演(FWI),给出了阿尔卑斯山脉岩石圈的三维密度、P波和S波速度模型.

      由于AxiSEM通过模拟四个二维全矩张量源的波场合成任意台站处的位移波场(Nissen-Meyer et al., 2007a, 2007b, 2008),相比于DSM的三维计算,大大提高了计算效率. 相比于FK方法,由于无平面波假设,AxiSEM模拟出更加接近真实的地震波场. AxiSEM方法的另一个优点是采用了与岩石圈目标区域中相同GLL节点的离散网格,避免了任何边界插值引入的数值噪声(Monteiller et al., 2020).

    • 本文系统总结了远震波场模拟的常用方法. 传统的远震波场模拟方法主要基于一维地球模型,根据其特点可以分为解析法、半解析法和数值法. 解析法主要采用积分变换方法,在高频条件下得到波传播问题的解析解. 半解析法主要借助于数值方法求解波传播问题的积分解. 数值法则主要基于一维球对称模型,对波动方程数值求解. 这些方法具有各自的特点,运算效率较高,但无法模拟三维非均匀介质中地震波的传播.

      随着软、硬件的快速发展,以及越来越多的高质量宽频带和短周期数据,人们对地下结构成像的精度要求也越来越高,仅利用传统方法已无法满足这一要求. 因此,研究人员提出了将全球一维地球模型的快速算法和局部高精度、高效率三维数值方法相结合的混合方法,实现了远震波场数值模拟的高效计算. 根据两种算法的特点可以有多种结合,得到不同的混合方法. 其中,GRT-FD混合法由于采用有限差分法来计算局部区域的三维非均匀地震波场,因而难以适应剧烈的地表起伏并处理自由边界条件,并且难以考虑地球的曲率变化因而难以适应区域尺度的地震波场高精度计算;FK-SEM混合法由于采用FK方法计算全球一维地球模型中的地震波场,无法考虑地球曲率且无绝对到时;DSM-SEM和AxiSEM-SEM混合方法由于采用SEM实现局域区域的三维非均匀地震波场数值模拟,通过设计自适应网格能够有效处理地表起伏并适应地球的曲率变化等,两者作为更为精确的远震波场混合模拟方法,有望在3D区域尺度波形反演、偏移成像中得到更广泛的应用.

      致谢

      感谢中国科学技术大学孙道远教授、美国密歇根州立大学(Michigan State University)陈敏(Chen Min)博士有益的建议. 感谢两位匿名审稿人的宝贵建议,对稿件提升很有帮助.

参考文献 (130)

目录

    /

    返回文章
    返回