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

黏弹地球地震变形理论研究进展和展望

唐河 孙文科

引用本文:
Citation:

黏弹地球地震变形理论研究进展和展望

    通讯作者: 孙文科, sunw@ucas.ac.cn
  • 中图分类号: P315.01, P315.4, P315.2

Progress and prospect of deformation theory in the viscoelastic earth

    Corresponding author: Sun Wenke, sunw@ucas.ac.cn
  • CLC number: P315.01, P315.4, P315.2

  • 摘要:Love(1911)研究了自重球体的弹性变形后,基于不同的黏弹性地球模型,许多科学家都对地震变形问题进行了深入研究,主要发展了基于半无限空间和球形地球模型的黏弹地球地震变形理论. 地震变形问题通常经积分变换、基函数展开等技术处理后,简化为求解满足特定震源和地表边界条件的常微分方程组问题. 针对这一独特的数学物理边值问题,本文以全解析、半解析和数值积分解等求解形式概述了近几十年发展的基于规则几何形态地球模型的黏弹地球地震变形理论,并讨论了各种方法的特点. 此外,针对三维地球模型,本文也简单回顾了目前的研究进展和存在问题. 总之,本文综述了过去半个多世纪以来的黏弹地球地震变形理论的发展历史、研究现状和最新进展,并讨论目前存在的问题和未来的发展趋势.
  • 图 1  分层球形地球模型中黏弹地震变形转化为复域等效弹性变形(一致性原理)示意图. (a)为Burgers体黏弹性变形. (b)为复域内等效弹性变形. Burgers体由Maxwell体和Kelvin体串联得到,弹簧原件表示弹性变形,阻尼器表示黏性变形,$\rho (r)$为密度,$\lambda {\rm{(}}r{\rm{)}}$$\mu (r)$为Lamé参数,$\eta (r)$为黏度,$\lambda {\rm{(}}s{\rm{)}}$$\mu (s)$为复域Lamé参数,它们均是半径r的函数,下标“K”和“M”分别表示Kelvin体和Maxwell体. 地球模型中:内核为弹性固体,外核为液态,地壳和地幔为黏弹性介质

    Figure 1.  The diagram of transforming a viscoelastic seismic deformation problem into an equivalent elastic deformation problem in complex domain (the correspondence principle) in a layered spherical Earth model. (a) the viscoelastic deformation of the Burgers body. (b) the equivalent elastic deformation in complex domain. The Burgers body is obtained by a series connection of a Maxwell body and a Kelvin body. Spring element represents elastic deformation; damper represents viscous deformation; $\rho (r)$ is density, $\lambda {\rm{(}}r{\rm{)}}$ and $\mu (r)$ are Lamé parameters, $\eta (r)$ is viscosity; $\lambda {\rm{(}}s{\rm{)}}$ and $\mu (s)$ is the complex Lamé parameters. Subscripts "K" and "M" represent Kelvin body and Maxwell body, respectively. All parameters are functions of radius r. In the Earth model: the inner core is elastic solid; the outer core is liquid; the crust and mantle are viscoelastic media

    图 2  由复Love数或Green函数,经逆Laplace变换计算时间域内Love数或Green函数的方法示意图. 负半实轴上的灰色点表示奇异点,“…”表示省略的奇异点.(a)中的黑色圆圈表示用Bromwich回路积分计算该奇异点处的留数,向上的有向直线表示逆Laplace积分定义式,对应路径${\rm{Re}} {\rm{(}}s{\rm{) = }}c$;(b)中阴影中的五角星表示Post-Widder方法计算点;(c)中有向矩形回路表示回路积分路径;(d)中半圆形表示积分核近似方法的路径,其内虚线表示纵向积分的实际路径,黑色点表示采样位置,$a$为引入参数,$t$为时间

    Figure 2.  The method of calculating Love numbers or Green function in time domain by inverse Laplace transform from complex Love numbers or Green function. The grey points on the negative semireal axis denote singular points. “…” denotes omitted singular points. The black circle in (a) represents the residue at the singular point calculated by Bromwich contour integral. The upward directed straight line ${\rm{Re}} {\rm{(}}s{\rm{) = }}c$ represents the inverse Laplace transform by the definition. In (b) the pentagram in the shadow represents the calculation point of the Post-Widder method. The directed rectangular in (c) represents the contour integration path. The semicircle in (d) represents the path of the integral kernel approximation method. Here the dotted line inside (d) represents the actual path of the vertical integration, the black dot represents the sampling position, $a$ is the introduced parameter, and $t$ is the time

  • [1]

    Aki K, Richards P G. Quantitative Seismology: Theory and Methods[M]. W. H. Freeman, 1980.
    [2]

    Bekaert D P S, Segall P, Wright T J, et al. A network inversionfilter combining GNSS and InSAR for tectonic slip modeling [J]. Journal of Geophysical Research-Solid Earth, 2016, 121(3): 2069-2086. doi: 10.1002/2015JB012638
    [3]

    Bonafede M, Ferrari C. Analytical models of deformation and residual gravity changes due to a Mogi source in a viscoelastic medium [J]. Tectonophysics, 2009, 471(1-2): 4-13. doi: 10.1016/j.tecto.2008.10.006
    [4]

    Broerse T, Riva R, Simons W, et al. Postseismic GRACE and GPS observations indicate a rheology contrast above and below the Sumatra slab [J]. Journal of Geophysical Research-Solid Earth, 2015, 120(7): 5343-5361. doi: 10.1002/2015JB011951
    [5]

    Cambiotti G. Joint estimate of the coseismic 2011 Tohoku earthquake fault slip and post-seismic viscoelastic relaxation by GRACE data inversion[J]. Geophysical Journal International, 2020, 220(2): 1012-1022.
    [6]

    Cambiotti G, Barletta V R, Bordoni A, et al. A comparative analysis of the solutions for a Maxwell Earth: the role of the advection and buoyancy force[J]. Geophysical Journal International, 2009, 176(3): 995-1006. doi: 10.1111/j.1365-246X.2008.04034.x
    [7]

    Cambiotti G, Klemann V, Sabadini R. Compressible viscoelastodynamics of a spherical body at long timescales and its isostatic equilibrium [J]. Geophysical Journal International, 2013, 193(3): 1071-1082. doi: 10.1093/gji/ggt026
    [8]

    Cambiotti G, Sabadini R. The compressional and compositional stratifications in Maxwell earth models: the gravitational overturning and the long-period tangential flux [J]. Geophysical Journal International, 2010a, 180(2): 475-500. doi: 10.1111/j.1365-246X.2009.04434.x
    [9]

    Cambiotti G, Sabadini R, Bordoni A. The continuous relaxation spectrum of Maxwell Earth models: a new method [J]. Egu General Assembly, 2010b, 12:9288.
    [10]

    Cambiotti G, Sabadini R, Yuen D A. Time-dependent geoid anomalies at subduction zones due to the seismic cycle[J]. Geophysical Journal International, 2018, 212(1): 139-150. doi: 10.1093/gji/ggx421
    [11]

    Cannelli V, Melini D, Piersanti A, et al. Application of the Post-Widder Laplace inversion algorithm to postseismic rebound models [J]. Nuovo Cimento Della Societa Italiana Di Fisica C-Colloquia on Physics, 2009, 32(2): 197-200.
    [12]

    Chao B F, Ding H. Spherical harmonic stacking for the singlets of Earth's normal modes of free oscillation [J]. Geophysical Ressearch Letters, 2014, 41(15): 5428-5435. doi: 10.1002/2014GL060700
    [13]

    Chao B F, Gross R S. Changes in the Earth's rotation and low-degree gravitational field induced by earthquakes [J]. Geophysical Journal International, 1987, 91(3): 569-596. doi: 10.1111/j.1365-246X.1987.tb01659.x
    [14]

    Chen W, Ray J, Li J C, et al. Polar motion excitations for an Earth model with frequency-dependent responses: 1. A refined theory with insight into the Earth's rheology and core-mantle coupling [J]. Journal of Geophysical Research: Solid Earth, 2013, 118(9): 4975-4994. doi: 10.1002/jgrb.50314
    [15]

    Chen W, Shen W B, Han J C, et al. Free wobble of the triaxial earth: Theory and comparisons with international earth rotation service (IERS) data [J]. Surveys in Geophysics, 2009, 30(1): 39-49. doi: 10.1007/s10712-009-9057-3
    [16] 程威,胡小刚. 南极上地幔结构异常对长周期自由振荡的影响 [J]. 地球物理学报, 2018, 61(8): 3211-3218. doi: 10.6038/cjg2018K0057

    Cheng W, Hu X G. The effect of upper mantle structure on long period Earth's free oscillation at the South Polar[J]. Chinese Journal of Geophysics, 61(8): 3211-3218 (in Chinese). doi: 10.6038/cjg2018K0057
    [17]

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

    Davis J L, Elósegui P, Mitrovica J X, et al. Climate‐driven deformation of the solid Earth from GRACE and GPS [J]. Geophysical Ressearch Letters, 2004, 31(24): 24605-24608. doi: 10.1029/2004GL021435
    [19]

    Ding H. Attenuation and excitation of the ∼6 year oscillation in the length-of-day variation [J]. Earth and Planetary Science Letters, 2019, 507(131-139. doi: 10.1016/j.jpgl.2018.12.003
    [20]

    Ding H, Chao B F. Data stacking methods for isolation of normal-mode singlets of Earth's free oscillation: Extensions, comparisons, and applications [J]. Journal of Geophysical Research: Solid Earth, 2015, 120(7): 5034-5050. doi: 10.1002/2015JB012025
    [21]

    Dziewonski A M, Anderson D L. Preliminary reference Earth model [J]. Physics of the Earth and Planetary Interiors, 1981, 25(4): 297-356. doi: 10.1016/0031-9201(81)90046-7
    [22]

    Fang M, Hager B H. A singularity-free approach to post glacial rebound calculations [J]. Geophysical Ressearch Letters, 1994, 21(19): 2131-2134. doi: 10.1029/94GL01886
    [23]

    Fang M, Hager B H. The singularity mystery associated with a radially continuous Maxwell viscoelastic structure [J]. Geophysical Journal International, 1995, 123(3): 849-865. doi: 10.1111/j.1365-246X.1995.tb06894.x
    [24]

    Farrell W E. Deformation of the Earth by surface loads [J]. Reviews of Geophysics & Space Physics, 1972, 10(3): 761-797.
    [25]

    Feng M X, Bie L D, Rietbrock A. Probing the rheology of continental faults: decade of post-seismic InSAR time-series following the 1997 Manyi (Tibet) earthquake [J]. Geophysical Journal International, 2018, 215(1): 600-613. doi: 10.1093/gji/ggy300
    [26]

    Fernández J, Rundle J B. Gravity changes and deformation due to a magmatic intrusion in a two-layered crustal model [J]. Journal of Geophysical Research: Solid Earth, 1994, 99(B2): 2737-2746. doi: 10.1029/93JB02449
    [27]

    Fernández J, Yu T T, Rundle J B. Deformation produced by a rectangular dipping fault in a viscoelastic-gravitational layered earth model. Part I: Thrust fault—FLTGRV and FLTGRH FORTRAN programs [J]. Computers & Geosciences, 1996, 22(7): 735-750.
    [28]

    Fu G Y, Sun W K. Effects of spatial distribution of fault slip on calculating co-seismic displacement: Case studies of the Chi-Chi earthquake (M.*?>=>W7.6) and the Kunlun earthquake (M.*?>=>W7.8) [J]. Geophysical Ressearch Letters, 2004, 31(21):L21601.
    [29]

    Fu G Y, Sun W K. Surface coseismic gravity changes caused by dislocations in a 3-D heterogeneous earth [J]. Geophysical Journal International, 2008, 172(2): 479-503. doi: 10.1111/j.1365-246X.2007.03684.x
    [30]

    Fu G Y, Sun W K, Fukuda Y, et al. Effects of Earth’s curvature and radial heterogeneity in dislocation studies: Case studies of the 2008 Wenchuan earthquake and the 2004 Sumatra earthquake [J]. Earthquake Science, 2010a, 23(4): 301-308. doi: 10.1007/s11589-010-0727-5
    [31]

    Fu G Y, Sun W K, Fukuda Y C, et al. Coseismic displacements caused by point dislocations in a three-dimensional heterogeneous, spherical earth model [J]. Geophysical Journal International, 2010b, 183(2): 706-726. doi: 10.1111/j.1365-246X.2010.04757.x
    [32]

    Gaver D P. Observing stochastic processes and approximate transform inversion [J]. Operations Research, 1966, 14(3): 444-459. doi: 10.1287/opre.14.3.444
    [33]

    Gharti H N, Langer L, Tromp J. Spectral-infinite-element simulations of coseismic and post-earthquake deformation [J]. Geophysical Journal International, 2019a, 216(2): 1364-1393. doi: 10.1093/gji/ggy495
    [34]

    Gharti H N, Langer L, Tromp J. Spectral-infinite-element simulations of earthquake-induced gravity perturbations [J]. Geophysical Journal International, 2019b, 217(1): 451-468. doi: 10.1093/gji/ggz028
    [35]

    Gilbert F. Excitation of the normal modes of the earth by earthquake sources [J]. Geophysical Journal International, 1971, 22(2): 223-226. doi: 10.1111/j.1365-246X.1971.tb03593.x
    [36]

    Gilbert F, Backus G. A computational proplem encountered in a study of the earth's normal modes[C]// Proceedings of the December 9-11, 1968, fall joint computer conference, part II on-AFIPS '68 (Fall, part II). San Francisco, California; ACM. 1968: 1273-1277.
    [37]

    Han S C, Sauber J, Pollitz F. Postseismic gravity change after the 2006-2007 great earthquake doublet and constraints on the asthenosphere structure in the central Kuril Islands [J]. Geophysical Ressearch Letters, 2016, 43(7): 3169-3177. doi: 10.1002/2016GL068167
    [38]

    Hanyk L, Matyska C, Yuen D A. Initial-value approach for viscoelastic responses of the earth’s mantle. In: Wu Patrick (ed) Dynamics of the Ice Age Earth: A Modern Perspective[M]. Trans Tech Publication Ltd, Switzerland, 1998: 135-154.
    [39]

    Hanyk L, Moser J, Yuen D A, et al. Time-domain approach for the transient responses in stratified viscoelastic earth models [J]. Geophysical Ressearch Letters, 1995, 22(10): 1285-1288. doi: 10.1029/95GL01087
    [40]

    Hanyk L. Viscoelastic response of the Earth: Initial-value approach[D]. Charles University, 1999.
    [41]

    Hanyk L, Yuen D A, Matyska C. Initial-value and modal approaches for transient viscoelastic responses with complex viscosity profiles [J]. Geophysical Journal International, 1996, 127(2): 348-362. doi: 10.1111/j.1365-246X.1996.tb04725.x
    [42]

    Hashima A, Takada Y, Fukahata Y, et al. General expressions for internal deformation due to a moment tensor in an elastic/viscoelastic multilayered half-space [J]. Geophysical Journal International, 2008, 175(3): 992-1012. doi: 10.1111/j.1365-246X.2008.03837.x
    [43]

    Haskell N A. The dispersion of surface waves in multilayered media [J]. Bullseismol Socam, 1951, 43:1-2.
    [44]

    Heki K, Miyazaki S, Tsuji H. Silent fault slip following an interplate thrust earthquake at the Japan Trench[J]. Nature, 1997, 386(6625): 595-598. doi: 10.1038/386595a0
    [45]

    Hetland E A, Hager B H. Postseismic and interseismic displacements near a strike‐slip fault: A two‐dimensional theory for general linear viscoelastic rheologies [J]. Journal of Geophysical Research: Solid Earth, 2005, 110:B10401. doi: 10.1029/2005JB003689
    [46]

    Hu Y. Three-dimensional viscoelastic finite element model for postseismic deformation of the great 1960 Chile earthquake [J]. Journal of Geophysical Research, 2004, 109:B12403. doi: 10.1029/2004JB003163
    [47] 胡小刚,薛秀秀, 郝晓光. 利用地球简正模耦合研究上地幔过渡区方位各向异性 [J]. 地球物理学报, 2012, 55(6): 1903-1911. doi: 10.6038/j.issn.0001-5733.2012.06.011

    Hu X G, Xue X X, Hao X G. Study of azimuthal anisotropy in the transition zone of the Earth’s upper mantle with the coupling of normal modes [J]. Chinese Journal of Geophysics. 2012, 55(6): 1903-1911 (in Chinese). doi: 10.6038/j.issn.0001-5733.2012.06.011
    [48]

    Imanishi Y, Sato T, Higashi T, et al. A network of superconducting gravimeters detects submicrogal coseismic gravity changes [J]. Science, 2004, 306(5695): 476-478. doi: 10.1126/science.1101875
    [49] 江颖,徐建桥, 孙和平. 深内部地球结构对内核平动振荡本征周期的影响[J]. 地球物理学报, 2014, 57(4): 1041-1048. doi: 10.6038/cjg20140403

    Jiang Y, Xu J Q, Sun H P. The influence of deep interior structure on the eigenperiod of inner core’s translational oscillations[J]. Chinese Journal of Geophysics. 2014, 57(4): 1041-1048 (in Chinese). doi: 10.6038/cjg20140403
    [50]

    Jonsson S, Segall P, Pedersen R, et al. Post-earthquake ground movements correlated to pore-pressure transients [J]. Nature, 2003, 424(6945): 179-183. doi: 10.1038/nature01776
    [51]

    Kanamori H, Anderson D L. Importance of physical dispersion in surface wave and free oscillation problems: Review [J]. Reviews of Geophysics, 1977, 15(1): 105-112. doi: 10.1029/RG015i001p00105
    [52]

    Langer L, Gharti H N, Tromp J. Impact of topography and three-dimensional heterogeneity on coseismic deformation [J]. Geophysical Journal International, 2019, 217(2): 866-878. doi: 10.1093/gji/ggz060
    [53]

    Liu T, Fu G Y, She Y W, et al. Green’s functions for post-seismic strain changes in a realistic earth model and their application to the Tohoku-Oki M.*?>=>W 9.0 earthquake [J]. Pure Appl Geophys, 2018, 176(9): 3929-3949.
    [54]

    Longman I M. A Green's function for determining the deformation of the Earth under surface mass loads: 1. Theory [J]. Journal of Geophysical Research, 1962, 67(2): 845-850. doi: 10.1029/JZ067i002p00845
    [55]

    Longman I M. A Green's function for determining the deformation of the Earth under surface mass loads: 2. Computations and numerical results [J]. Journal of Geophysical Research, 1963, 68(2): 485-496. doi: 10.1029/JZ068i002p00485
    [56]

    Love A E H. Some Problems of Geodynamics[M]. Cambridge: Cambridge University Press, 1911.
    [57] 鹿璐, 李小凡, 陈楠, 等. 核幔边界地幔侧不同尺度横向非均匀对地球自由振荡的影响 [J]. 地球物理学进展, 2017, 32(04): 1465-1473. doi: 10.6038/pg20170407

    Lu L, Li X F, Chen N, et al. Influence of different scales of lateral inhomogeneous structures on the mantle side of the CMB on the Earth's free oscillation [J]. Progress in Geophysics, 32(4): 1465-1473 (in Chinese). doi: 10.6038/pg20170407
    [58]

    Matsu'ura M, Tanimoto T. Quasi-static deformations due to an inclined, rectangular fault in a viscoelastic half-space [J]. Journal of Physics of the Earth, 1980, 28(1): 103-118. doi: 10.4294/jpe1952.28.103
    [59]

    Mcconnell R K. Isostatic adjustment in a layered earth [J]. Journal of Geophysical Research, 1965, 70(20): 5171-5188. doi: 10.1029/JZ070i020p05171
    [60]

    Melini D, Cannelli V, Piersanti A, et al. Post-seismic rebound of a spherical Earth: new insights from the application of the Post-Widder inversion formula [J]. Geophysical Journal International, 2008, 174(2): 672-695. doi: 10.1111/j.1365-246X.2008.03847.x
    [61]

    Metois M, Socquet A, Vigny C. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone [J]. Journal of Geophysical Research-Solid Earth, 2012, 117(B3): 406-421.
    [62]

    Molodenskiy S M. The influence of horizontal inhomogeneities in the mantle on the amplitude of tidal oscillations [J]. Physics of the Solid Earth, 1977, 13(2): 77-80.
    [63]

    Molodenskiy S M. The effect of lateral heterogeneities upon the tides [J]. BIM Fevrier, 1980, 80:4833-4850.
    [64]

    Moore J D, Yu H, Tang C H, et al. Imaging the distribution of transient viscosity after the 2016 M.*?>=>W7.1 Kumamoto earthquake [J]. Science, 2017, 356(6334): 163-167. doi: 10.1126/science.aal3422
    [65]

    Okada Y. Surface deformation due to shear and tensile faults in a half-space [J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154.
    [66]

    Okada Y. Internal deformation due to shear and tensile faults in a half-space [J]. Bulletin of the Seismological Society of America,1992, 82(2): 1018-1040.
    [67]

    Okubo S. Asymptotic solutions to the static deformation of the Earth - I. Spheroidal mode [J]. Geophysical Journal International, 1988, 92(1): 39–51. doi: 10.1111/j.1365-246X.1988.tb01119.x
    [68]

    Okubo S. Potential and gravity changes raised by point dislocations [J]. Geophysical Journal International, 1991, 105(3): 573-586. doi: 10.1111/j.1365-246X.1991.tb00797.x
    [69]

    Okubo S. Gravity and potential changes due to shear and tensile faults in a half-space [J]. Journal of Geophysical Research: Solid Earth, 1992, 97(B5): 7137-7144. doi: 10.1029/92JB00178
    [70]

    Piersanti A, Spada G, Sabadini R. Global postseismic rebound of a viscoelastic Earth: Theory for finite faults and application to the 1964 Alaska earthquake [J]. Journal of Geophysical Research-Solid Earth, 1997, 102(B1): 477-492. doi: 10.1029/96JB01909
    [71]

    Piersanti A, Spada G, Sabadini R, et al. Global post-seismic deformation [J]. Geophysical Journal International, 1995, 120(3): 544-566. doi: 10.1111/j.1365-246X.1995.tb01838.x
    [72]

    Pollitz F F. Postseismic relaxation theory on the spherical earth [J]. Bulletin of the Seismological Society of America, 1992, 82(1): 422-453.
    [73]

    Pollitz F F. Gravitational viscoelastic postseismic relaxation on a layered spherical Earth [J]. Journal of Geophysical Research-Solid Earth, 1997, 102(B8): 17921-17941. doi: 10.1029/97JB01277
    [74]

    Pollitz F F. Transient rheology of the uppermost mantle beneath the Mojave Desert, California [J]. Earth and Planetary Science Letters, 2003a, 215(1-2): 89-104. doi: 10.1016/S0012-821X(03)00432-1
    [75]

    Pollitz F F. Post-seismic relaxation theory on a laterally heterogeneous viscoelastic model [J]. Geophysical Journal International, 2003b, 155(1): 57-78. doi: 10.1046/j.1365-246X.2003.01980.x
    [76]

    Pollitz F F, Burgmann R, Banerjee P. Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth [J]. Geophysical Journal International, 2006, 167(1): 397-420. doi: 10.1111/j.1365-246X.2006.03018.x
    [77]

    Post E L. Generalized differentiation [J]. Transactions of the American Mathematical Society, 1930, 32(1-4): 723-781.
    [78]

    Qiu Q, Moore J D P, Barbot S, et al. Transient rheology of the Sumatran mantle wedge revealed by a decade of great earthquakes [J]. Nature Communications, 2018, 9(1): 995-1007. doi: 10.1038/s41467-018-03298-6
    [79]

    Rodkin M V, Kaftan V I. Post-seismic relaxation from geodetic and seismic data [J]. Geodesy and Geodynamics, 2017, 8(1): 13-16. doi: 10.1016/j.geog.2017.01.001
    [80]

    Rogister Y, Rochester M G. Normal-mode theory of a rotating earth model using a Lagrangian perturbation of a spherical model of reference [J]. Geophysical Journal International, 2004, 159(3): 874-908. doi: 10.1111/j.1365-246X.2004.02447.x
    [81]

    Rundle J B. Static elastic-gravitational deformation of a layered half-space by point couple sources [J]. Journal of Geophysical Research, 1980, 85(Nb10): 5355-5363. doi: 10.1029/JB085iB10p05355
    [82]

    Rundle J B. Viscoelastic-gravitational deformation by a rectangular thrust-fault in a layered earth [J]. Journal of Geophysical Research, 1982, 87(Nb9): 7787-7796. doi: 10.1029/JB087iB09p07787
    [83]

    Sabadini R, Yuen D A, Boschi E. The effects of post-seismic motions on the moment of inertia of a stratified viscoelastic earth with an asthenosphere[J]. Geophysical Journal of the Royal Astronomical Society, 1984, 79(3): 727-745. doi: 10.1111/j.1365-246X.1984.tb02865.x
    [84] 单新建, 屈春燕, 龚文瑜, 等. 2017年8月8日四川九寨沟7.0级地震InSAR同震形变场及断层滑动分布反演 [J]. 地球物理学报, 2017, 60(12): 4527-4536. doi: 10.6038/cjg20171201

    Shan X J, Qu C Y, Gong W Y, et al.Coseismic deformation field of the Jiuzhaigou M.*?>=>S7.0 earthquake from Sentinel-1A InSAR data and fault slip inversion. Chinese Journal of Geophysics,60(12): 4527-4536 (in Chinese). doi: 10.6038/cjg20171201
    [85]

    Spada G, Boschi L. Using the Post-Widder formula to compute the Earth's viscoelastic Love numbers [J]. Geophysical Journal International, 2006, 166(1): 309-321. doi: 10.1111/j.1365-246X.2006.02995.x
    [86]

    Steketee J A. On Volterra's dislocations in a semi-infinite elastic medium [J]. Canadian Journal of Physics, 1958a, 36(2): 192-205. doi: 10.1139/p58-024
    [87]

    Steketee J A. Some geophysical application of the elasticity theory of dislocations [J]. Canadian Journal of Physics, 1958b, 36(9): 1168-1198. doi: 10.1139/p58-123
    [88]

    Suito H, Freymueller J T. A viscoelastic and afterslip postseismic deformation model for the 1964 Alaska earthquake [J]. Journal of Geophysical Research-Solid Earth, 2009, 114(B11): 404-426.
    [89]

    Sun T H Z, Wang K L. Viscoelastic relaxation following subduction earthquakes and its effects on afterslip determination [J]. Journal of Geophysical Research-Solid Earth, 2015, 120(2): 1329-1344. doi: 10.1002/2014JB011707
    [90]

    Sun W K. Asymptotic theory for calculating deformations caused by dislocations buried in a spherical earth: geoid change [J]. Journal of Geodesy, 2003, 77(7-8): 381-387. doi: 10.1007/s00190-003-0335-4
    [91]

    Sun W K. Asymptotic solution of static displacements caused by dislocations in a spherically symmetric Earth [J]. Journal of Geophysical Research-Solid Earth, 2004a, 109(B5): 402-419.
    [92]

    Sun W K. Short Note: Asymptotic theory for calculating deformations caused by dislocations buried in a spherical earth — gravity change [J]. Journal of Geodesy, 2004b, 78(1): 76-81.
    [93] 孙文科. 地震位错理论[M].北京: 科学出版社, 2012a.

    Sun W K. Theory on Earthquake Dislocation[M]. Beijing: Science Press, 2012a (in Chinese).
    [94] 孙文科. 地震位错理论在地震学研究中的作用与存在的问题[J]. 国际地震动态, 2012b, 000(6): 17.

    Sun W K. The theory of earthquake dislocation in seismic research functions and problems [J]. International Earthquake Dynamics.2012b, 000(6): 17 (in Chinese).
    [95]

    Sun W K, Dong J. Geo-center movement caused by huge earthquakes [J]. Journal of Geodynamics, 2014, 76(1):67-73.
    [96]

    Sun W K, Okubo S. Surface potential and gravity changes due to internal dislocations in a spherical earth-I. Theory for a point dislocation[J]. Geophysical Journal International, 1993, 114(3): 569-592. doi: 10.1111/j.1365-246X.1993.tb06988.x
    [97]

    Sun W K, Okubo S. Surface potential and gravity changes due to internal dislocations in a spherical earth-II. Application to a finite fault [J]. Geophysical Journal International, 2002, 132(1): 79-88. doi: 10.1046/j.1365-246x.1998.00400.x
    [98]

    Sun W K, Zhou X. Coseismic deflection change of the vertical caused by the 2011 Tohoku-Oki earthquake (M.*?>=>W 9.0) [J]. Geophysical Journal International, 2012, 189(2): 937-955. doi: 10.1111/j.1365-246X.2012.05434.x
    [99]

    Takeuchi H, Saito M. Seismic surface waves [J]. Methods in Computational Physics Advances in Research & Applications, 1972, 11(1): 217-295.
    [100]

    Tanaka Y, Hasegawa T, Tsuruoka H, et al. Spectral-finite element approach to post-seismic relaxation in a spherical compressible Earth: application to gravity changes due to the 2004 Sumatra–Andaman earthquake [J]. Geophysical Journal International, 2015, 200(1): 299-321. doi: 10.1093/gji/ggu391
    [101]

    Tanaka Y, Klemann V, Fleming K, et al. Spectral finite element approach to postseismic deformation in a viscoelastic self-gravitating spherical Earth [J]. Geophysical Journal International, 2009, 176(3): 715-739. doi: 10.1111/j.1365-246X.2008.04015.x
    [102]

    Tanaka Y, Okuno J, Okubo S. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (I) - vertical displacement and gravity variation [J]. Geophysical Journal International, 2006, 164(2): 273-289. doi: 10.1111/j.1365-246X.2005.02821.x
    [103]

    Tanaka Y, Okuno J, Okubo S. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (II)-horizontal displacement [J]. Geophysical Journal International, 2007, 170(3): 1031-1052. doi: 10.1111/j.1365-246X.2007.03486.x
    [104]

    Tang H, Sun W K. Asymptotic expressions for changes in the surface co-seismic strain on a homogeneous sphere [J]. Geophysical Journal International, 2017, 209(1): 202-225.
    [105]

    Tang H, Sun W K. Asymptotic co- and post-seismic displacements in a homogeneous Maxwell sphere [J]. Geophysical Journal International, 2018a, 214(1): 731-750. doi: 10.1093/gji/ggy174
    [106]

    Tang H, Sun W K. Closed-form expressions of seismic deformation in a homogeneous Maxwell Earth model [J]. Journal of Geophysical Research-Solid Earth, 2018b, 123(7): 6033-6051. doi: 10.1029/2018JB015594
    [107]

    Tang H, Sun W K. New method for computing postseismic deformations in a realistic gravitational viscoelastic Earth model [J]. Journal of Geophysical Research-Solid Earth, 2019, 124(5): 5060-5080. doi: 10.1029/2019JB017368
    [108]

    Tang H, Dong J, Sun W K. An approximate method to simulate post-seismic deformations in a realistic earth model[M]. International Association of Geodesy Symposia, 2020a.
    [109]

    Tang H, Zhang L, Chang L, et al. Optimized approximate inverse Laplace transform for geo-deformation computation in viscoelastic Earth model [J]. Geophysical Journal International, 2020b, 223(1): 444–453. doi: 10.1093/gji/ggaa322
    [110]

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

    Valsa J, Brancik L. Approximate formulae for numerical inversion of Laplace transforms [J]. International Journal of Numerical Modelling-Electronic Networks Devices and Fields, 1998, 11(3): 153-166. doi: 10.1002/(SICI)1099-1204(199805/06)11:3<153::AID-JNM299>3.0.CO;2-C
    [112]

    Vermeersen L L A, Sabadini R. A new class of stratified viscoelastic models by analytical techniques [J]. Geophysical Journal International, 1997, 129(3): 531-570. doi: 10.1111/j.1365-246X.1997.tb04492.x
    [113]

    Vermeersen L L A, Sabadini R, Spada G. Analytical visco-elastic relaxation models [J]. Geophysical Ressearch Letters, 1996, 23(7): 697-700. doi: 10.1029/96GL00620
    [114]

    Wang H S. Surface vertical displacements, potential perturbations andgravity changes of a viscoelastic earth model induced by internal point dislocations[J]. Geophysical Journal International, 1999a, 137(2): 429-440.
    [115]

    Wang H, Wu P. Effects of lateral variations in lithospheric thickness and mantle viscosity on glacially induced surface motion on a spherical, self-gravitating Maxwell Earth [J]. Earth and Planetary Science Letters, 2006a, 244(3-4): 576-589. doi: 10.1016/j.jpgl.2006.02.026
    [116]

    Wang H S, Wu P. Effects of lateral variations in lithospheric thickness and mantle viscosity on glacially induced relative sea levels and long wavelength gravity field in a spherical, self-gravitating Maxwell Earth [J]. Earth & Planetary Ence Letters, 2006b, 249(3-4): 368-383.
    [117]

    Wang H S, Wu P. Analytical approach for the toroidal relaxation of viscoelastic earth [J]. Geophysical Journal International, 2006c, 167(1): 1-19. doi: 10.1111/j.1365-246X.2006.02980.x
    [118] 汪汉胜, Patrick Wu, Hugo Schotman, 等. 横向非均匀地球负荷问题的CLFE有限元算法的有效性 [J]. 地球物理学报, 2006, 49(6): 1657-1664. doi: 10.3321/j.issn:0001-5733.2006.06.012

    Wang H S, Wu P, Schotman H, et al. Validation of the coupled Laplace-Finite-Element method for the loading problem of a laterally heterogeneous sphencal Earth[J]. Chinese Journal of Geophysics, 2006, 49(6): 1657-1664 (in Chinese). doi: 10.3321/j.issn:0001-5733.2006.06.012
    [119]

    Wang K. The seismogenic zone of subduction thrust faults—17.Elastic and viscoelastic models of crustal deformation in subduction earthquake cycles[M].Dixon T, Moore J C(ed.). New York: Columbia University Press, 2007a: 540-575.
    [120]

    Wang K, Hu Y, He J. Deformation cycles of subduction earthquakes in a viscoelastic Earth [J]. Nature, 2012a, 484(7394): 327-332. doi: 10.1038/nature11032
    [121]

    Wang L, Shum C K, Simons F J, et al. Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry [J]. Geophysical Ressearch Letters, 2012b, 39(7): L07301.
    [122] 王启欣, 江在森, 武艳强, 等. 不同模型下地震位错理论的对比及其应用进展综述 [J]. 地震学报, 2015, (4): 690-704. doi: 10.11939/jass.2015.04.014

    Wang Q X, Jiang Z S, Wu Y Q, et al. A review on comparison and progress in applications of earthquake dislocation theories based on different models[J]. Acta Seismologica Sinica, 2015, 37(4): 690-704 (in Chinese). doi: 10.11939/jass.2015.04.014
    [123]

    Wang R, Wang H S. A fast converging and anti-aliasing algorithm for Green's functions in terms of spherical or cylindrical harmonics [J]. Geophysical Journal International, 2007, 170(1): 239-248. doi: 10.1111/j.1365-246X.2007.03385.x
    [124]

    Wang R J. A simple orthonormalization method for stable and efficient computation of Green's functions [J]. Bulletin of the Seismological Society of America, 1999b, 89(3): 733-741.
    [125]

    Wang R J. The dislocation theory: a consistent way for including the gravity effect in (visco)elastic plane-earth models [J]. Geophysical Journal International, 2005a, 161(1): 191-196. doi: 10.1111/j.1365-246X.2005.02614.x
    [126]

    Wang R J. On the singularity problem of the elastic-gravitational dislocation theory applied to plane-Earth models [J]. Geophysical Ressearch Letters, 2005b, 32(6):L06307
    [127]

    Wang R J, Lorenzo-Martín F, Roth F. PSGRN/PSCMP—a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory [J]. Computers & Geosciences, 2006, 32(4): 527-541.
    [128]

    Wang W Z, Wang H, Liu Y C, et al. A comparative study of the methods for calculation of surface elastic deformation [C]//Proceedings of the Institution of Mechanical Engineers Part J-Journal of Engineering Tribology, 2003, 217(J2): 145-153.
    [129]

    Wen Y M, Li Z H, Xu C J, et al. Postseismic motion after the 2001 M.*?>=>W 7.8 Kokoxili earthquake in Tibet observed by InSAR time series [J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B8): 405-409.
    [130]

    Widder D V. The inversion of the laplace integral and the related moment problem [J]. Transactions of the American Mathematical Society, 1934, 36(1-4): 107-200.
    [131]

    Wong M C, Wu P. Using commercial finite-element packages for the study of Glacial Isostatic Adjustment on a compressible self-gravitating spherical earth – 1: Harmonic loads [J]. Geophysical Journal International, 2019, 217(3): 1798-1820. doi: 10.1093/gji/ggz108
    [132]

    Wu P. Using commercial finite element packages for the study of earth deformations, sea levels and the state of stress [J]. Geophysical Journal International, 2004, 158(2): 401-408. doi: 10.1111/j.1365-246X.2004.02338.x
    [133]

    Xu C J, Gong Z, Niu J M. Recent developments in seismological geodesy [J]. Geodesy and Geodynamics, 2016, 7(3): 157-164. doi: 10.1016/j.geog.2016.04.009
    [134]

    Xu C Y, Sun W K. Earthquake-origin expansion of the Earth inferred from a spherical-Earth elastic dislocation theory [J]. Geophysical Journal International, 2014, 199(3): 1655-1661. doi: 10.1093/gji/ggu364
    [135]

    Yang H Y, Tromp J. Synthetic free-oscillation spectra: an appraisal of various mode-coupling methods [J]. Geophysical Journal International, 2015a, 203(2): 1179-1192. doi: 10.1093/gji/ggv349
    [136]

    Yang J Y, Zhou X, Yi S, et al. Determining dislocation love numbers using GRACE satellite mission gravity data [J]. Geophysical Journal International, 2015b, 203(1): 257-269. doi: 10.1093/gji/ggv265
    [137]

    Yu T T, Rundle J B, Fernández J. Deformation produced by a rectangular dipping fault in a viscoelastic-gravitational layered earth model. Part II: Strike-slip fault—STRGRV and STRGRH FORTRAN programs [J]. Computers & Geosciences, 1996, 22(7): 751-764.
    [138]

    Zhang G Q, Shen W B, Xu C Y, et al. Coseismic gravity and displacement signatures induced by the 2013 Okhotsk M.*?>=>W8.3 earthquake [J]. Sensors, 2016, 16(9):1410. doi: 10.3390/s16091410
    [139] 钟敏, 罗少聪, 申文斌, 等. 地球主惯性矩的精密确定及三轴分层地球自转动力学研究 [J]. 科技资讯, 2016, 14(18): 181-181. doi: 10.3969/j.issn.1672-3791.2016.18.102

    Zhong M, Luo S C, Shen W B, et al. Precise determination of the Earth's principal moment intertia and study of tri-axial layered earth rotation dynamics[J]. Science & Technology Information, 2016, 14 (18): 181-181 (in Chinese). doi: 10.3969/j.issn.1672-3791.2016.18.102
    [140]

    Zhou J, Pan E, Bevis M. A point dislocation in a layered, transversely isotropic and self-gravitating Earth. Part I: analytical dislocation Love numbers [J]. Geophysical Journal International, 2019a, 217(3): 1681-1705. doi: 10.1093/gji/ggz110
    [141]

    Zhou J, Pan E, Bevis M. A point dislocation in a layered, transversely isotropic and self-gravitating Earth — Part II: accurate Green's functions [J]. Geophysical Journal International, 2019b, 219(3): 1717-1728. doi: 10.1093/gji/ggz392
    [142]

    Zhou J, Pan E, Bevis M. A point dislocation in a layered, transversely isotropic and self-gravitating Earth – Part III: Internal deformation [J]. Geophysical Journal International, 2020, (3):1681-1705.
    [143]

    Zhou X, Cambiotti G, Sun W, et al. The coseismic slip distribution of a shallow subduction fault constrained by prior information: the example of 2011 Tohoku (M.*?>=>W 9.0) megathrust earthquake [J]. Geophysical Journal International, 2014, 199(2): 981-995. doi: 10.1093/gji/ggu310
    [144]

    Zhou X, Sun W K, Zhao B, et al. Geodetic observations detecting coseismic displacements and gravity changes caused by the M.*?>=>W= 9.0 Tohoku-Oki earthquake [J]. Journal of Geophysical Research: Solid Earth, 2012, 117: B05408.
  • [1] 吴忠良王龙李丽张晓东邵志刚李营孙珂车时 . 中国地震科学实验场:地震预测与系统设计. 地球与行星物理论评, 2021, 52(): 1-5. doi: 10.16738/j.dqyxx.2021-028
    [2] 孙伟家王一博魏勇赵亮 . 火星地震学与内部结构研究. 地球与行星物理论评, 2021, 52(): 1-13. doi: 10.16738/j.dqyxx.2021-016
  • 加载中
图(2)
计量
  • 文章访问数:  1316
  • HTML全文浏览量:  615
  • PDF下载量:  67
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-07-02
  • 网络出版日期:  2020-08-27
  • 刊出日期:  2021-01-01

黏弹地球地震变形理论研究进展和展望

摘要: Love(1911)研究了自重球体的弹性变形后,基于不同的黏弹性地球模型,许多科学家都对地震变形问题进行了深入研究,主要发展了基于半无限空间和球形地球模型的黏弹地球地震变形理论. 地震变形问题通常经积分变换、基函数展开等技术处理后,简化为求解满足特定震源和地表边界条件的常微分方程组问题. 针对这一独特的数学物理边值问题,本文以全解析、半解析和数值积分解等求解形式概述了近几十年发展的基于规则几何形态地球模型的黏弹地球地震变形理论,并讨论了各种方法的特点. 此外,针对三维地球模型,本文也简单回顾了目前的研究进展和存在问题. 总之,本文综述了过去半个多世纪以来的黏弹地球地震变形理论的发展历史、研究现状和最新进展,并讨论目前存在的问题和未来的发展趋势.

English Abstract

    • 近几十年来,以地震为研究对象、大地测量学为观测手段的新兴交叉学科——地震大地测量学不断发展,有力地促进了人们对地震过程及地震变形的认识,同时也推动了地震变形问题的理论研究. 通过不同时空尺度的大地测量观测数据,人们已经充分了解以断层闭锁、同震破裂、孔隙弹性变形、震后黏性松弛等为具体内容的地震孕育、破裂及震后变形过程,同时也伴随着一系列相关地震位错理论的发展与进步.

      地震周期通常是一个复杂的过程,大体上分为同震、震后和震间变形几个主要子过程,持续时间从几秒到几十年不等. 研究震间—同震—震后的物理过程,将帮助我们估计应变能的积累状况和地壳—地幔的物性参数,理解长期构造过程和地震循环周期之间的关系,提高对大地震特别是俯冲带地震的预警能力(Wang, 2007a; Cambiotti et al., 2018). 此外,理解地震变形可以帮助我们了解地下的应力状态、地壳—地幔的力学特性和地下介质参数,加深对地下结构的认知(Wang, 2007a; Wang et al., 2012a). 总之,现代大地测量技术在长时间尺度的变形观测,可以从同震变形、震后—震间变形角度,加深我们对地震孕育—发生机制的理解以及地壳—地幔物态性质的认知.

      同震形变是指在几分钟内的快速应变能释放,产生显著的各种地球物理场变化,这些变化都可以被现代大地测量观测技术观测和检测到,如GNSS、地表重力仪、合成孔径雷达、空间重力卫星等;这些变化也可以被弹性地震位错理论很好地解释(Sun, 1993).

      震间和震后变形比较复杂,包含了多种不同的物理机制,反映了包括断层闭锁耦合、余震、慢滑动、孔隙弹性响应和黏弹性松弛的复杂过程(Heki et al., 1997; Piersanti et al., 1997; Jonsson et al., 2003; Pollitz et al., 2006; Metois et al., 2012; Wang et al., 2012b; Tanaka et al., 2015). 仅通过观测资料来区分某个地震后的变形物理机制是非常困难的,还需要相应的理论模拟进行解释、补充或约束. 近期的研究发现,如果忽略黏性地幔的蠕动,则所有地表变形都必须由断层的闭锁和滑动来解释,从而导致对闭锁状态的错误描述(Wang et al., 2012a). 在一次大的俯冲地震后,地壳形变持续进行,变形模式演化过程非常复杂(Wang et al., 2012a),因为这种震后变形主要是由地震破裂引起的黏弹性应力松弛和断层不同部位的持续滑动(余滑)或再锁(Suito and Freymueller, 2009)造成的. 人们通常认为余滑主要引起破裂带附近的短期(几年)变形,孔隙弹性变形则集中于断层源附近且持续时间更短,而黏弹性松弛主要产生长期和大范围的变形(Wang et al., 2012a). 然而,目前很难用常规大地测量数据来检验这一假设的有效性.

      总之,震后变形包含丰富的地球物理信息,可以用来确定构造区域的流变结构,估算历史大地震的震源参数,研究应力松弛对其他地震的触发等问题. 震间变形和震后变形在物理上虽然有不同的地球物理机制,但在观测上却很难将其区分,即,如何把这几种不同物理机制的变形有效分离仍然是一个尚未解决的科学问题.

      随着大地测量学观测技术的发展和高质量观测数据的积累,震后变形问题受到了越来越多的关注(Wen et al., 2012; Broerse et al., 2015; Han et al., 2016; Rodkin and Kaftan, 2017). 这些现代大地测量为科学家研究地震变形和分离物理信号提供了一定可能性(Davis et al., 2004; Bekaert et al., 2016; Xu et al., 2016). 然而,某些过程,如余滑和黏弹性松弛产生的变形,其在时间和空间域具有复杂的特征(Sun and Wang, 2015; Moore et al., 2017; Qiu et al., 2018),仅从观测数据中分离它们仍然是相当困难的. 为了解决这个问题,一个较为有效的方法是,首先从理论上计算地震后的黏弹性松弛变形,再从观测值中减去该松弛变形,从而分离出震后余滑变形;另一种方法可以是,利用GNSS等观测确定余滑效应,在观测之中扣除它,从而得到黏弹松弛变形. 无论如何,建立和发展一个完善的黏弹地球模型的地震变形理论(位错理论)是重要的理论基础.

      几十年来,地震位错理论不断发展和完善,用来描述球对称、非自转、黏弹性、各向同性和分层黏弹地球的地震准静态变形过程. 位错理论的核心内容是应用多种数学物理方程的求解技巧,解算平衡方程、应力-应变关系、位场泊松方程等组成的定解边值问题,给出地球内部地震震源所产生的位移场、应力场、重力场等变化. 位错理论经过几十年的发展,所采用的地球模型逐渐由半无限空间发展到逼近三维真实地球,微分方程求解过程逐渐简化,数值奇异、震荡等数学问题逐渐被解决. 这些进展按照数学解算方法划分为多种表现形式,即,全解析、半解析和全数值形式. 本文在介绍黏弹地球地震变形问题的基础理论(第1节)后,主要以这三种表现震后变形形式的顺序介绍黏弹地球变形问题的主要进展(第2节),并简单介绍三维地球模型变形理论(第3节),以及目前存在的问题和未来发展展望.

    • 地震断层破裂对震源附近的介质产生了强烈的应力扰动. 对于较大震级的地震,特别是发生于地壳深部的大地震,强大的应力扰动会导致上地幔介质发生持续的变形响应,表现为流变性介质发生持续流动,内部应力不断发生调整. 由此产生的地球物理变形场是时间的非线性函数,其求解通常比较困难. 为了简化问题,一般假定地幔介质为线性流变介质,例如,Kelvin体、Burgers体、Maxwell体或广义Maxwell体等,利用这些介质模型的线性应力-应变关系来约束地球的黏弹变形. 地震产生的黏弹变形场通常由平衡方程、泊松方程、应力-应变关系来表述,由于应力-应变关系依赖于时间,其变形场的求解为四维问题. 近几十年来,地震变形理论的发展大致可以归结为:(1)介质模型由弹性过渡为黏弹性;(2)几何模型由均匀球或半无限空间过渡为分层球;(3)介质的可压缩性、自重效应等逐渐加以考虑和解决;(4)计算复杂性逐渐被简化.

      黏弹地球模型中流变介质的变形,本质上是通过求解数学物理方程来确定的. 对于黏弹性的四维变形问题,根据对时间的依赖关系不同,可以分为两种主要方法:第一种方法是将应力-应变关系在时间域内进行离散化处理,得到递推形式的方程,从而可以先求解某个时刻点的变形,之后递推得到后续时刻的变形(Hanyk et al., 1995);第二种方法是将应力-应变关系进行Laplace变换,得到同弹性介质形式相同的复域应力-应变关系,再求解Laplace域内的“等效”弹性变形的复域解,最后进行逆Laplace变换得到时间域内的解,这称之为一致性原理(Mcconnell et al., 1965). 图1给出一致性原理的示意图.

      图  1  分层球形地球模型中黏弹地震变形转化为复域等效弹性变形(一致性原理)示意图. (a)为Burgers体黏弹性变形. (b)为复域内等效弹性变形. Burgers体由Maxwell体和Kelvin体串联得到,弹簧原件表示弹性变形,阻尼器表示黏性变形,$\rho (r)$为密度,$\lambda {\rm{(}}r{\rm{)}}$$\mu (r)$为Lamé参数,$\eta (r)$为黏度,$\lambda {\rm{(}}s{\rm{)}}$$\mu (s)$为复域Lamé参数,它们均是半径r的函数,下标“K”和“M”分别表示Kelvin体和Maxwell体. 地球模型中:内核为弹性固体,外核为液态,地壳和地幔为黏弹性介质

      Figure 1.  The diagram of transforming a viscoelastic seismic deformation problem into an equivalent elastic deformation problem in complex domain (the correspondence principle) in a layered spherical Earth model. (a) the viscoelastic deformation of the Burgers body. (b) the equivalent elastic deformation in complex domain. The Burgers body is obtained by a series connection of a Maxwell body and a Kelvin body. Spring element represents elastic deformation; damper represents viscous deformation; $\rho (r)$ is density, $\lambda {\rm{(}}r{\rm{)}}$ and $\mu (r)$ are Lamé parameters, $\eta (r)$ is viscosity; $\lambda {\rm{(}}s{\rm{)}}$ and $\mu (s)$ is the complex Lamé parameters. Subscripts "K" and "M" represent Kelvin body and Maxwell body, respectively. All parameters are functions of radius r. In the Earth model: the inner core is elastic solid; the outer core is liquid; the crust and mantle are viscoelastic media

      第一种方法,主要是对于不同的流变性模型,需要将时间域应力-应变关系离散化,得到特定递推形式的应力-应变关系,之后将其与平衡方程、泊松方程一起做矢量球谐展开,得到一系列时间域内离散化的非齐次形式的常微分方程组,通过求解弹性变形的常微分方程组,并进行递推,可以计算后续时刻的变形(Hanyk et al., 1995; Hanyk et al., 1996; Hanyk et al., 1998; Hanyk, 1999). 因为该方法是时间域内的递推算法,计算比较复杂,相邻时间节点间计算误差的累积效应会导致长时间尺度的变形求解算法不稳定,近些年来并未有很多后续研究.

      目前应用比较广泛的方法是借助一致性原理(Mcconnell et al., 1965),将四维的黏弹变形问题转换为复域内的三维等效弹性变形问题来求解. 具体的求解思路为:(1)特定黏弹性介质的应力-应变关系作Laplace变换,转化为与弹性介质变形相一致的形式,此时黏弹Lamé参数转化为依赖于复频率的函数,黏弹问题转换为复域弹性变形问题;(2)在Laplace域中求解等效的弹性问题获得复域解,再将复域解作逆Laplace变换回到时间域中,得到黏弹变形问题的解.

      黏弹介质的变形问题和弹性介质的变形问题非常相似,二者仅在介质的本构(应力-应变)关系上有区别. 弹性问题的应力仅随应变变化,而黏弹性问题的应力既依赖于应变也依赖于应变率,它们在本质上有显著区别;但对于平衡微分方程和泊松方程,它们的形式是完全一样的.

      下面以常用的Maxwell介质和Burgers介质(如图1a中所示)为例,具体说明如何使用一致性原理. 三维介质中,Maxwell体的本构关系中应力随着应变率变化,具体形式为:

      ${\dot{ \tau }} + \dfrac{{{\mu _{\rm{M}}}}}{{{\eta _{\rm{M}}}}}\left({{{\tau }} - \dfrac{1}{3}{\rm{Tr}}({{\tau }}){{I}}} \right) = \lambda {{I}}\nabla \cdot {\dot{ u}} + {\mu _{\rm{M}}}\left({\nabla {\dot{ u}} + {{(\nabla {\dot{ u}})}^{\rm{T}}}} \right)$

      (1)

      式中,${{{\tau}}}$为应力,${{u}}$为位移,符号上的“˙”表示对时间$t$求导,$\lambda $${\mu _{\rm{M}}}$为Lamé参数,${\eta _{\rm{M}}}$为黏滞参数(黏度),所有地球介质参数均是半径$r$的函数,为了方便简介而省略了(下同),${{I}}$为单位矩阵,Tr为矩阵的迹. 引入新的复变参数$s$,作Laplace变换得到式(1)在复数域中的表达形式为:

      ${{\tau }}{\rm{(}}s{\rm{)}} = \tilde \lambda {\rm{(}}s{\rm{)}}{{I}}\nabla \cdot {{u}}{\rm{(}}s{\rm{)}} + \tilde \mu {\rm{(}}s{\rm{)}}\left[ {\nabla {{u}}{\rm{(}}s{\rm{)}} + {{(\nabla {{u}}{\rm{(}}s{\rm{)}})}^{\rm{T}}}} \right]$

      (2)

      它和弹性介质的本构关系形式完全一致,这为求解黏弹性变形问题带来了极大的便利. 需要强调的是,此时Lamé参数$\tilde \lambda {\rm{(}}s{\rm{)}}$$ \tilde \mu \left(s \right) $是复函数:

      $ \left\{ {\begin{array}{*{20}{l}} {\tilde \mu (s) = \dfrac{{{\mu _{\rm{M}}}s}}{{s + {\mu _{\rm{M}}}/{\eta _{\rm{M}}}}}}\\ {\tilde \lambda (s) = \dfrac{{\lambda s + k{\mu _{\rm{M}}}/{\eta _{\rm{M}}}}}{{s + {\mu _{\rm{M}}}/{\eta _{\rm{M}}}}}} \end{array}} \right.,\;\;k = \lambda + 2{\mu _{\rm{M}}}/3 $

      (3)

      式中,$s$为Laplace变量,量纲为时间的倒数;${\eta _{\rm{M}}}$为Maxwell黏度;$k$为体变模量,它不随时间发生松弛变化.

      Burgers介质体的本构关系更为复杂,但经Laplace变换后,仍然可以写成与弹性介质一样的形式. 和Maxwell介质不同,Burgers介质体对应的复Lamé参数较为复杂,可以表达为:

      $\left\{ {\begin{array}{*{20}{l}} {\tilde \lambda (s) = \lambda + {{2{\mu_{\rm{M}}}} / 3} - \tilde \mu (s)} \\ {\tilde \mu (s) = \dfrac{{s\left({s + \dfrac{{{\mu_{\rm{K}}}}}{{{\eta_{\rm{K}}}}}} \right){\mu_{\rm{M}}}}}{{{s^2} + s\left({\dfrac{{{\mu_{\rm{M}}} + {\mu_{\rm{K}}}}}{{{\eta_{\rm{K}}}}} + \dfrac{{{\mu_{\rm{M}}}}}{{{\eta_{\rm{M}}}}}} \right) + \dfrac{{{\mu_{\rm{M}}}{\mu_{\rm{K}}}}}{{{\eta_{\rm{M}}}{\eta_{\rm{K}}}}}}}} \end{array}} \right.$

      (4)

      式中,带下标“K”和“M”的参数分别对应Kelwin体和Maxwell体,依次表示快速蠕变和稳态蠕变过程涉及的物理量. 实际上,Burgers体是比Maxwell体更一般的流变介质模型,可以蜕化为Maxwell体和Kelvin体:当${\eta_{\rm{K}}} \to \infty$时,Burgers退化为Maxwell体;当${\eta_{\rm{M}}} \to \infty$时,Burgers退化为Kelvin体.

      借助一致性原理,我们可以看到黏弹介质变形的基本方程在Laplace域中与弹性介质中的形式完全一致,唯一的区别是黏弹性介质中的Lamé参数不再是实数,而是一个随流变介质模型不同而异、并且依赖于Laplace变量$s$的复函数. 从而,先前不同学者提出的弹性变形问题的微分方程、边界条件、求解方法均可以在Laplace域中直接使用.

      在Laplace域求解了复域微分方程后,可得到复位错Love数和Green函数,它们可以统一记为$G(s)$. 获得Laplace域内的解$G(s)$仅仅是完成了黏弹性变形计算的一半工作,我们最终目的是获得时间域的解$G(t)$. 为此,我们需要计算如下逆Laplace变换:

      $G{\rm{(}}t{\rm{) = }}\dfrac{1}{{2\text{π} i}}\int_{c - i\infty }^{c + i\infty } {\dfrac{{G{\rm{(}}s{\rm{)}}}}{s}{{\rm{e}}^{st}}{\rm{d}}s} $

      (5)

      式中,$i$为虚数单位,$c$(收敛横坐标)需要大于被积函数的所有奇异点.

      上述公式(5)在黏弹变形问题研究中是非常重要的公式之一,数学上称为黎曼—梅林反演公式. 在地球模型接近于真实地球时,Laplace域内的解$G{\rm{(}}s{\rm{)}}$的性质非常复杂,存在诸多“奇异点”(如图2中各子图坐标轴上的小圆圈),这些奇异点的密集程度和分布模式取决于地球模型的分层结构以及黏弹参数的复杂性;更重要的是积分核中含有指数函数${{\rm{e}}^{st}}$,它随着时间$t$的增加而呈现明显的振荡现象,亦存在数值向上溢出的风险. 所以,实际计算中,上述逆Laplace变换的计算非常困难.

      图  2  由复Love数或Green函数,经逆Laplace变换计算时间域内Love数或Green函数的方法示意图. 负半实轴上的灰色点表示奇异点,“…”表示省略的奇异点.(a)中的黑色圆圈表示用Bromwich回路积分计算该奇异点处的留数,向上的有向直线表示逆Laplace积分定义式,对应路径${\rm{Re}} {\rm{(}}s{\rm{) = }}c$;(b)中阴影中的五角星表示Post-Widder方法计算点;(c)中有向矩形回路表示回路积分路径;(d)中半圆形表示积分核近似方法的路径,其内虚线表示纵向积分的实际路径,黑色点表示采样位置,$a$为引入参数,$t$为时间

      Figure 2.  The method of calculating Love numbers or Green function in time domain by inverse Laplace transform from complex Love numbers or Green function. The grey points on the negative semireal axis denote singular points. “…” denotes omitted singular points. The black circle in (a) represents the residue at the singular point calculated by Bromwich contour integral. The upward directed straight line ${\rm{Re}} {\rm{(}}s{\rm{) = }}c$ represents the inverse Laplace transform by the definition. In (b) the pentagram in the shadow represents the calculation point of the Post-Widder method. The directed rectangular in (c) represents the contour integration path. The semicircle in (d) represents the path of the integral kernel approximation method. Here the dotted line inside (d) represents the actual path of the vertical integration, the black dot represents the sampling position, $a$ is the introduced parameter, and $t$ is the time

      为了完成公式(5)的积分计算,很多科学家做了大量研究工作,其目的就是如何处理这一无穷震荡积分. 图2汇总了目前为止典型的处理积分公式(5)的主要方法,代表了目前为止的研究黏弹变形问题的主要历史和进展,后面第2节还将稍微详细介绍和讨论. 图2a给出最早提出并且最直观的解决方法,即,通过特征方程求根找到Laplace域中$G{\rm{(}}s{\rm{)}}$的所有奇异点${s_i}$,用Bromwich回路积分方法计算奇异点对应的变形,再将每个奇异点对应的变形即某个模态相互累加,得到最终的黏滞性变形(Piersanti et al., 1995; Pollitz, 1997; Vermeersen and Sabadini, 1997). 图2b是一种特殊的回避极点的计算方法Post-Widder方法(Post, 1930; Widder, 1934; Gaver, 1966). 该方法的核心思想是借用离散的Post公式,在五角星位置进行复域采样并加权求和,来计算逆Laplace变换,避免了Bromwich回路积分的计算,从而绕过了留数定理和求根过程. 图2c给出了复平面中作闭合回路积分的方法,它可以一次性考虑所有模态对最终变形的贡献(Tanaka et al., 2006; Tanaka et al., 2007);一旦确定了合适的积分路径,就可以沿矩形路径进行逆Laplace积分,从复Love数中获得真实的时间域变形. 图2d是一种最近提出的比较新颖且有效的近似办法——纵轴积分方法(Tang and Sun, 2019; Tang et al., 2020b). 这个方法本质上是一个逆Laplace变换的近似算法,其核心思想是将积分核近似为一个分式函数,在解析地应用留数定理后,将逆Laplace积分转化为在图2d虚线上采样的一个交错级数求和问题,从而彻底回避了寻找复域Love数或者处理Green函数的奇异点的麻烦,并且因为对积分核的特殊处理,也不存在数值溢出问题.

      总之,在处理逆Laplace积分公式(5)以及求解微分方程过程中科学家们发展了丰富有效的数值技巧、成熟完善的处理策略,最终形成了不同的理论分支,用于计算时间域Green函数. 获得Green函数后,对有限断层做二维数值积分可得到任意断层在任意位置产生的地震变形场. 以下我们从全解析解、半解析解,以及数值解等形式具体介绍各种位错理论(地震变形)的研究进展.

    • 黏弹性地球模型中地震变形的求解,需要设定具体的地球几何模型和物理参数模型,而不同的物理模型会导致变形问题求解难度不同,目前已经发展了全解析的、半解析的以及数值的方法来求解点震源产生的变形(Green函数). 目前仅对分层半无限空间和球形地球模型等几何上规则的模型,发展了一系列较成熟的理论方法. 对于半无限空间和球形地球模型,其数学处理方式比较类似,均使用了基函数展开或者积分变换等常用方法,但为了行文方便,我们下面分别阐述.(1)对于分层半无限空间地球模型,一般将复域三维空间的变形方程做Hankel变换(Aki and Richards, 1980),将极坐标中的微分方程转化为波数域内的常微分方程组,该常微分方程组仅与深度、波数、介质的物性参数相关,且存在解析的通解. 通过Thomson-Haskell传播矩阵的方法(Thomson, 1950; Haskell, 1951),可以得到在多分层模型中,满足特定震源和边界条件的解析解. 之后将波数域内的解作逆Laplace变换和逆Hankel变换得到极坐标下的时域Green函数.(2)对于分层球形地球模型,一般将复域三维空间的变形方程做矢量球谐函数展开,得到一组仅与地球半径、球谐阶次、介质的物性参数相关,以展开系数为未知函数的复域线性微分方程组(Tanaka et al., 2006; Tanaka et al., 2007; Tang and Sun, 2019);通过解析或者数值的方法求解后,得到地表和内部的球谐展开系数,它们称为复域位错Love数;对复域Love数作逆Laplace变换以得到时间域内的位错Love数(Tanaka et al., 2006; Melini et al., 2008; Cambiotti et al., 2013; Tang and Sun, 2019);之后将Love数与球谐函数相乘并通过适当的级数加速方法(Farrell, 1972; Wang and Wang, 2007b)得到变形Green函数. 根据不同的物理模型中位错Love数和Green函数是否可以解析求解,下面我们以全解析方法、半解析方法、数值方法的形式分别介绍黏弹地球变形问题的研究进展.

    • 无论是数学上还是实际应用上,人们都渴望得到全解析解,因为解析解数学上严谨,计算速度快,计算精度高. 此外,通过全解析的Green函数可以直观地理解物理模型和地震变形场之间的函数关系,深化对变形问题的认知. 然而,全解析变形Green函数仅存在于比较简单几何形状的地球模型,如均匀半无限空间和均质球体.

    • 对于半无限空间,Steketee(1958a, 1958b)将位错概念引入到地震学中,并给出了位错所产生位移的一般表达式. 此后,许多学者都讨论了半无限弹性空间中剪切位错造成的弹性变形,特别是Okada(1985, 1992),在系统总结前人成果基础上给出了半无限空间弹性模型中内部任意矩形断层运动产生的地表位移、倾斜和应变等解析表达式. Okubo(1991, 1992)给出了大地水准面和重力变化的闭合形式的Green函数. 借助一致性原理,通过解析逆Laplace变换可以便捷地得到黏弹性介质如Maxwell体和Burgers体中的对应Green函数. 例如,假定半无限空间填充Maxwell介质,Matsu'ura和Tanimoto(1980)给出了由矩形断层产生的地表位移、应变、倾斜等变形的解析Green函数. 此外,Hetland和Hager(2005)考虑了一个弹性层覆盖于黏弹性半无限空间上的平层地球模型,给出了无限走滑断层产生变形的解析Green函数. Bonafede和Ferrari(2009)提出了在具有Maxwell流变介质的半无限空间地球模型中,由Mogi源引起的位移变形和重力变化的解析Green函数.

    • 对于均质球,早在100多年前,Love(1911)就导出了一般弹性自重球体变形的解析解. Okubo(1988)进一步将该解析解渐近展开为球谐阶次的幂级数,从而推导出其渐近解,为此后各种地震变形的渐近解的推导提供了基础. 结合地震变形问题中特定的地表和内部位错源处的边界条件(Takeuchi and Saito, 1972),求解地震变形的非齐次微分方程组,可得到渐近地震位错Love数,将位错Love数与球谐函数相乘并对阶次解析求和,便可以得到全解析的Green函数. Sun(2003, 2004a, 2004b)以及Tang和Sun(2017)具体得到了所有弹性位错Love数的渐近级数形式,以及同震位移、重力、大地水准面、应变等物理量的渐近Green函数,将渐近解方法发展为一套成熟完备的理论. 需要指出的是,幂级数形式的位错Love数仅是准确解的级数近似,在表示高阶变形上有效,阶次越高,近似效果越好(Okubo, 1988).

      此后,Tang和Sun(2018a)将该渐近方法与一致性原理(Mcconnell, 1965)相结合,先求解复域等效弹性变形问题,再通过解析逆Laplace变换,得到了均质Maxwell黏弹自重球体中的震后变形的渐近解. 他们发现微分方程组渐近解在Maxwell均质球中存在复域二阶极点,使得均质Maxwell球体的地震变形出现线性发散项,但该线性项在特征松弛时间范围内相比于复域一阶极点的效应可以忽略,表明这些黏弹解析Green函数在数学上是有效的.

      对于黏弹均质球,若忽略介质的自重效应,则变形微分方程的求解会进一步简化,Pollitz(1992)给出了简化情形下的一般解. Tang和Sun(2018b)进一步求解了该模型下的地震变形问题,给出了解析形式的位错Love数及其渐近展开式,得到了完全解析形式的震后变形Green函数. 此外,均质球或者分层球,其环型变形的微分方程比较简单,可以得到多分层模型中的解析解. Wang和Wu(2006c)得到了2~5层分层球形模型中的环型震后Love数,分析了变形模态和分层结构的关系. 然而,由于球型变形的复杂性,目前还没有解析解.

      应该指出的是,黏弹均质球的全解析Green函数有着计算简便的优势,可以作为理论参考,也可以用于定性分析,但是在解释观测数据上时应该考虑弹性地壳的影响.

    • 对于较为复杂的地球模型,例如,分层半无限空间或者层状球等介质模型,一般很难给出完全解析解. 为了克服这个困难,通常可以采用半解析半数值的处理方法,简称半解析解方法.

    • 针对分层半无限空间模型中的黏弹性地震变形,许多学者提出了不同的求解方法. Rundle(1980, 1982)提出了计算分层黏弹自重介质中点源和有限矩形源产生的同震和震后变形的数值方法. 随后Fernández和Rundle(1994)讨论了同样模型下地下膨胀源的弹性变形问题,Fernández等(1996)以及Yu等(1996)给出了两层介质中计算地表位移变形的Fortran 77程序. 但是Wang(2005)发现Rundle等(Rundle, 1980; Rundle, 1982; Fernández and Rundle, 1994; Fernández et al., 1996; Yu et al., 1996)错误地处理了重力效应,他指出,若考虑半无限空间中的重力效应和可压缩性,当波数小于某个特定值时,半无限空间模型的地震变形问题的解无法确定. 该困难可以通过多种近似处理来回避,如不考虑半空间介质的可压缩性或重力效应. 忽略介质的重力效应,Hashima等(2008)得到了黏弹性多分层半空间中由任意矩张量引起的内部变形场的一般表达式. 近似处理后,黏弹性变形问题是可解的,但随着震后时间增长,剪切模量随着时间松弛变小时,其数值解不稳定并造成结果的不准确. Wang(2005a, 2005b)沿用Longman(1962, 1963)的方法,假定分层介质中的密度梯度满足初始的静压力平衡条件,即Adams-Williamson条件,回避了数值不稳定性问题,并且在该近似下任意波数对应的变形都是可解的. 此外Wang(1999a)Wang等(2003)强调了Thomson-Haskell传播矩阵在直接使用时,通解在向下传递过程中会因为“大数吃小数”而失去线性独立性,使得特定的边界条件无法满足,为此,他们提出了将通解在每一层做简单的正交归一化处理,从而安全地将解向下传递. 之后,Wang等(2006)给出了完整的Fortran计算程序PSGRN/PSCMP,用于模拟分层半空间自重黏弹性模型的地表位移、应变、倾斜、重力和大地水准面等物理量的变化. 至此,平面分层黏弹性介质的地震变形理论趋于完善.

      应当指出,Wang等(2006)发表的分层数量不限的平层地球模型下的黏弹变形计算方法和程序,由于忽略了地球的曲率,在计算中比较简单,可以方便地应用于近场地震变形的计算和解释,但它在计算远场变形时存在忽略地球曲率所产生的误差.

    • Gilbert(1971)指出地震产生的静态位移场可以用自由振荡简正模的叠加来表示(参考图2a),主要需要两个步骤:(1)求解弹性—重力(或者黏弹性—重力)运动方程得到简正模解;(2)根据震源机制叠加简正模得到地震变形Green函数. 步骤二的过程相对简单,关键在于第一步. 目前简正模的求解主要有三个发展阶段. 首先是基于完全弹性均质球形地球模型的半解析半数值积分方法;其次是基于黏弹性(滞弹性)分层球形地球模型的求解;再次是考虑带椭球率和旋转的黏弹性(滞弹性)分层地球模型的求解(Dahlen and Tromp, 1998). 需要明确的是,滞弹性变形和黏弹性变形比较类似. 针对滞弹性问题,其Lamé参数是以赖于频率和品质因子的复变量,其对应的简正模态也是复域内的;而黏弹性问题其本构关系依赖于时间,经Laplace变换后本构关系和弹性情形形式相同,但其Lamé参数变为依赖于黏度的复变量,需要在Laplace域内求解. 有多种观测技术可以观测到地球的自由振荡,其在地球的内部结构研究中有重要应用(Chao and Ding, 2014; Ding and Chao, 2015; Ding, 2019);自由振荡简正模方法在求解滞弹性和黏弹性问题上有很长的发展历史,以下我们仅概述与黏弹性变形紧密相关的研究.

      针对分层球形地球模型,已有很多学者采用简正模叠加方法研究了地震黏弹性变形. 例如,Chao和Gross(1987)最早用简正模叠加的方法计算了弹性自重地球模型中地震产生的地球极移、日长变化和低阶重力场变化. 自Sabadini等(1984)开始,在假定介质不可压缩条件下,Vermeersen等(1996)Vermeersen和Sabadini(1997)通过简正模方法,得到了分层球形模型中的同震变形的半解析Love数. Pollitz(1992)Sabadini等(1984)Piersanti等(1995)使用传播矩阵和简正模叠加的方法,计算了分层、自重、不可压缩的黏弹地球模型中点源产生的地震变形. Wang(1999b)利用简正模方法半解析地计算了分层、自重、可压缩的黏弹地球模型中的地震变形. Cambiotti等(2009)Cambiotti和Sabadini(2010a)以及Cambiotti等(2013)对简正模方法深入研究后发现,如果假定地幔中的密度满足G-B假设(Gilbert and Backus, 1968),即每个均质层中的重力和半径成正比,则自重可压缩介质的黏弹性变形存在近似的解析解. 他们对简正模方法的有效拓展,使得可以考虑自重和可压缩性介质的黏弹性变形. 周江存团队(Zhou et al., 2019a; Zhou et al., 2019b; Zhou et al., 2020)发现,对G-B假设进行进一步扩展,即假定地核任何一层中密度是恒定的、重力随半径线性变化,而地幔任何一层中密度与半径成反比、重力恒定,则可以将分层自重可压缩球体中弹性变形微分方程转换为常系数微分方程组,从而得到解析形式的传播矩阵,用于快速计算位错Love数和Green函数. 这种方法可以和一致性原理结合使用,通过解析的逆Laplace变换得到半解析形式的震后位错Love数和Green函数.

      但以上Cambiotti团队和周江存团队的研究目前尚停留在简正模求解的第二阶段. 其实,考虑椭球率和旋转效应的简正模方法也有了很大的发展. 例如,Chen等(2009)钟敏等(2016)讨论了三轴椭球地球的自由振荡,Chen等(2013)研究了黏弹性和核幔耦合对地球自转轴定向的影响;Rogister和Rochester(2004)使用微扰方法,通过对静压平衡的球形地球模型做两次拉格朗日扰动,得到了旋转椭球地球模型的自有震荡简正模解;不同于微扰方法,Yang和Tromp(2015a)等基于复球谐函数和实球谐函数,提出了一种数学上更严密的方法用于求解非球状自转滞弹性地球模型的自由振荡简正模解. 这些研究考虑了三维几何形态或滞弹性结构对自由振荡的影响(Kanamori and Anderson,1977胡小刚等,2012江颖等,2014Ding and Chao,2015鹿璐等,2017程威和胡小刚,2018),但尚未应用到黏弹性地球的变形研究中.

      需要指出,由简正模方法来表示地震的静态变形比较复杂,因为它需要非常多的特征频率的变形,即简正模态,叠加求和才能表达一个静态位移场. 这类似于信号处理中,用一系列三角函数的组合来表达阶跃函数. 更重要的是,简正模计算在高阶变形和低频变形时会遇到数值困难,其高阶截断误差会造成近场变形计算的收敛问题. 而确定黏弹介质的变形问题变得更加繁杂,其主要的困难在于给定某一特征频率,我们需要通过数值寻根的办法找到所有满足复域特征方程的极点,并计算对应极点的各阶次变形. 根据Vermeersen等(1996)的研究,对于某个阶次的变形,即使在一个简单的分层地球模型中也存在无限多的极点. 在连续分层可压缩的地球模型中,寻找极点会更加困难,计算量也非常大. 因此,用简正模叠加的方法求解黏弹介质中的变形并非是一个好的选择.

    • Spada和Boschi(2006)Melini等(2008)的研究表明,通过应用Post-Widder形式的逆Laplace变换公式(Gaver, 1966),适当修改简正模计算静态变形的公式,可以半解析地解决震后变形问题. 该方法的主要优点是可以绕过数值办法寻找复数域极点的求根过程和逆Laplace回路积分(参考图2b). Post-Widder方法(Post, 1930Widder, 1934Gaver, 1966)的核心思想是借用离散的Post公式:

      $G(t,M) = \dfrac{{\ln 2}}{t}\sum\limits_{n = 1}^{2M} {{\zeta _n}} {G_h}(n{{\ln 2} / t})\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$

      $\begin{split} {\zeta _n} =& {(- 1)^{M + n}}\sum\limits_{j = [(n + 1)/2]}^{\min (n,M)} {\dfrac{{{j^{M + 1}}}}{{M!}}} \left({\begin{array}{*{20}{c}} M \\ j \end{array}} \right)\left({\begin{array}{*{20}{c}} {2j} \\ j \end{array}} \right)\left({\begin{array}{*{20}{c}} j \\ {n - j} \end{array}} \right)\\{G_h}(s) =& \dfrac{{G(s)}}{s} \\[-10pt]\end{split}$

      (6)

      (式中,$2M$为采样点个数,“[n]”为向下取整函数,$\left({{array}{*{20}{c}} M \\ j {array}} \right)$为组合数)来进行逆Laplace变换,避免了Bromwich回路积分的计算,从而绕过了留数定理和求根过程. 通过该方法,可以将时间域中的实Green函数表示为Laplace域中沿实轴采样的复Green函数的加权和(Spada and Boschi, 2006; Melini et al., 2008). 使用该方法,可以方便地使用简正模叠加和传播矩阵的方法,将已有的弹性变形计算程序改写为黏弹性变形计算程序,用于Maxwell和Burgers流变学地球模型中的黏弹变形的计算. 但该方法需要非常高的数值计算精度,支持IEEE扩展精度格式的30位有效数字的商业Fortran编译器也不能保证Post-Widder方法的成功应用(Melini et al., 2008),必须使用第三方扩展函数库来进行高精度的数值计算. Post-Widder方法本质上需要一个高精度的计算机系统,实施起来非常不便.

    • 对于复杂地球模型且考虑各种介质参数的情况,数值积分是一个重要途径. 尽管数值解不如解析解或半解析解简单,但是它可以给出最准确的计算结果,在实际应用时,高精度的数值解往往是最佳的选择. 在研究弹性地球同震变形时,孙文科团队(Sun and Okubo, 1993; Sun and Okubo, 2002; Sun and Zhou, 2012; Sun and Dong, 2014; Xu and Sun, 2014)考虑了球形分层地球模型以及介质的自重、可压缩性、径向分层结构,给出了同震物理场变化,如位移、重力、大地水准面等. 该团队还进一步发展了同震和震后变形的渐近解理论,考虑地球曲率效应并又保证了数学上的解析性,比Okada(1985, 1992)Okubo(1991, 1992)基于半无限空间的表达式物理上更加合理,可以作为分层球地震变形的高阶近似,来代替高阶数值Love,对于快速计算地震变形有实际应用价值.

    • 在考虑固态内核—液态外核—固体地幔等多物态分层结构,顾及球形地球的曲率效应,并严格考虑地球介质的可压缩性、自重效应、流变性等情形下,复域变形微分方程组需要严格的数值积分,如Runge-Kutta方法来由地心向上积分求解,得到复域数值的位错Love数,然后通过数值逆Laplace变换得到时域Love数,再通过对各阶次的Love数求和得到数值Green函数. 而数值逆Laplace变换计算,则是其中最关键的一步.

      根据Fang和Hager(1994)的研究,一般Maxwell黏弹性分层模型中的复域地球变形,其内在的奇异点分布复杂、存在分支点,简正模方法在数学上的严密性有待验证. 此外,Fang和Hager(1995)分析了奇异点的分布范围,指出黏度剖面复杂的地球模型,其复域变形Love数奇异点分布非常复杂. Tanaka等(2006, 2007)证明了一维分层黏弹性模型中,所有奇异点均分布于复域实轴上,克服了简正模方法无法同时考虑可压缩性和地球径向连续流变结构的固有困难,提出在复平面内作回路积分的方法,来同时考虑复杂模型下无数变形模态${s_i}$的贡献,得到了数值的Green函数,计算公式为:

      $G(t) = \sum\limits_i {{\rm{Res}} \left[ {{G_h}(s){{\rm{e}}^{st}},{s_i}} \right]}$

      (7)

      式中,$i$表示第$i$个奇异点${s_i}$,Res表示求留数,${G_h}(s)$${{G(s)} / s}$.

      该方法的主要思想是要确定一个适当的数值积分回路,以便包含复域所有奇异点. 实际上,积分路径是由以下权衡条件而给出的(参考图2c):(1)路径上复域Love数的变化比较平滑;(2)较弱的积分核振荡性. 前者要求路径远离实轴(较大的${s_2}$),后者则要求积分路径更接近原点(较小的$\left| {{s_1}} \right|$${s_2}$${s_3}$). 此外,矩形路径的左右边界,即参数${s_1}$${s_3}$,应通过数值寻根方案确定最大和最小的极点后给定. 如果对PREM模型(Dziewonski and Anderson, 1981)进行轻微修改,消除不稳定的分层结构(Cambiotti and Sabadini, 2010a),即没有大密度层覆盖于低密度层之上,则最大的极点应该为零;那么,只需要寻找最小的极点即可. 然而,Tanaka等(2006, 2007)声称在没有正极点的情况下计算的震后Love数在t=0时的结果与使用纯弹性理论(Sun and Okubo, 1993)计算的不一致. 因此,如果应用回路积分来获得PREM地球模型的同震和震后变形,必须要确定最大和最小的极点.

    • 为了更好的解决复杂地球模型下复域内有无数个奇异点对应多个变形模态的问题,并克服上述矩形积分路径不易确定的问题,Tang和Sun(2019)借用Valsa和Brancik(1998)提出的逆Laplace变换的近似计算方法,提出了从复域Love数计算时域Love数的新方法——纵轴积分方法.

      该方法的主要步骤为:(1)根据Valsa和Brancik(1998)的研究,引入参数$a$,逆Laplace变换中的积分核可以近似地展开为一个收敛的泰勒级数,即:

      ${{\rm{e}}^{st}} \approx \text{π} {{\rm{e}}^a}\sum\limits_{n = 0}^\infty {\dfrac{{{{(- 1)}^n}(n + {1 / 2})}}{{{{(n + {1 / 2})}^2}{{\text{π}} ^2} + {{(a - st)}^2}}}} $

      (8)

      式中,$n$为级数项数,$s$为Laplace变量,$t$为时间. 将该泰勒级数代入黎曼—梅林反演公式(5),可以将逆Laplace变换公式近似为一个交错级数求和的形式,而级数的每一项可以用非振荡的无穷积分来表达:

      $G(t) \!\approx\! \dfrac{{{{\rm{e}}^a}}}{{2i}}\!\sum\limits_{n = 0}^\infty\! {\left[ {{{(- 1)}^n}\left({n + \dfrac{1}{2}} \right)\!\int_{c - i\infty }^{c + i\infty } \!{\dfrac{{{G_h}(s)}}{{{{(n + {1 / 2})}^2}{{\text{π}} ^2} + {{(a - st)}^2}}}{\rm{d}}s} }\! \right]} $

      (9)

      式中,$a$为参数,$i$为虚数单位,$n$为级数项数,$c$为积分收敛横坐标,$s$为Laplace变量,$t$为时间,${G_h}(s)$${{G(s)} / s}$.

      (2)对上述的无穷积分补充一个半圆形的线积分,形成一个新的回路积分(参考图2d);(3)对新的回路积分解析地使用留数定理,将其简化为两个奇异点的留数之和;(4)将逆Laplace变化解析地表达为一个交错级数求和,并使用适当的加速技术(如欧拉变换)来加快级数的收敛,其最终的表达公式为:

      $G(t) \approx \dfrac{{{{\rm{e}}^a}}}{t}\sum\limits_{n = 1}^\infty {{{(- 1)}^n}{\rm{Im}} \left\{ {{G_h}\left[ {\dfrac{a}{t} + i\left({n - \dfrac{1}{2}} \right)\dfrac{{\text{π}} }{t}} \right]} \right\}} $

      (10)

      式中,各参数意义与公式(9)类似,${\rm{Im}} $表示取复数的虚部,${G_h}(s) $${{G(s)} / s} $.

      Tang和Sun(2019)发现时域Love数可以由复域Love数在与虚轴平行的直线上采样,并加权求和给出,同时该方法彻底回避了寻找奇异点和振荡积分的问题,可以应用于复杂地球模型中地震黏弹性变形的计算. 数值计算表明,新方法收敛快、数值稳定,仅需要大约40个复域Love数的采样便可以计算给定0~100年的震后变形. 该方法有如下优点:(1)可以方便地控制计算精度;(2)数值稳定性好、精度高;(3)不依赖于地球的黏弹性参数、适用性强.

    • 在这节里我们简单介绍一下三维地球模型位错理论的进展. 将分层球形地球的变形理论发展到三维黏弹性地球的变形理论是比较困难的,目前正处于一个逐步发展、不断探索的阶段. 目前为止的研究不是很多,但是主要有以下几项研究和进展.

      三维地球的黏弹性变形研究,主要有两种途径:(1)利用三维地球的自由振荡理论,由简正模叠加方法求解地震的静态变形场;(2)直接求解在零频下的三维地球变形方程. 第一种方法,可以充分利用目前三维地球的自由振荡理论. 地震学中已经发展了多种有效方法来计算横向不均匀性的三维地球模型中的地震波场. 通过简正模叠加,理论上可以得到地震的静态变形场,但需要求解非常高频的自由振荡简正模和高球谐阶次的变形,将面临数值震荡、计算精度问题,这可能导致经逆Laplace变换后的时域变形的精度不够. 此外,在近场Green函数的计算中可能需要很高的截断阶次,其逆Laplace的变换将难以保证精度. 这或许正是三维弹性地球自由振荡理论尚未在三维黏弹地球的地震变形中广泛使用的一个原因.

      第二种方法是在零频的条件下直接求解三维地球的静态变形. 需要强调的是,在零频条件下地球的变形方程需要重新调整,求解过程中需要特殊处理以便构建微分方程的独立解. 虽然这种方法需要更多的理论工作,但它比简正模方法更直接有效,其需要处理的数值问题也更少,计算效率高. 例如,付广裕和孙文科等便提出了以微扰理论为基础的三维地球模型地震变形的理论和方法(Fu and Sun, 2004; Fu and Sun, 2008; Fu et al., 2010a; Fu et al., 2010b). 其实,Molodenskiy(1977, 1980)最早利用微扰理论研究了地球固体潮汐变形问题,但是,Molodenskiy(1977, 1980)的微扰动理论需要假定“陆地—海洋”模型,这与真实的地球三维结构有较大偏离,限制了其在实际地震中的应用. Fu和Sun(2008)以及Fu等(2010)基于Molodenskiy(1977, 1980)的微扰动理论思想,探讨了三维非均匀结构地球模型的地震变形问题,其主要思路是将三维地球变形的解表达为球对称的解与横向非均匀性扰动解的叠加. 此时,地表的一般变形需要由6个独立震源产生变形的组合来计算:一个垂直断层走滑破裂、两个正交垂直断层倾滑破裂、一个水平断层垂直引张、两个垂直断层水平引张. 扰动解分为震源介质处的横向不均性和全球横向非均匀性构造产生影响,可以分别独立处理. 前者可以通过改变球对称地球模型中的震源参数来直接求解,而后者需要通过对平衡微分方程的变分法来求解. 他们进一步利用地震层析成像所得到的三维地球构造结果和岩石实验经验公式,实际计算了2011年日本东北大地震(M.*?>=>W9.0)所产生的同震位移和重力变化,探讨了三维构造对地震变形的影响. 该理论目前只适用于计算三维模型中的同震变形,尚未应用到黏弹地球介质的变形研究中. 此外,Pollitz(2003a, 2003b)利用了三维横向非均匀性的简正模计算结果,给出了震后变形的半解析半数值的计算方法,讨论了局部构造的三维效应,然而他的方法不适用于一般的全球三维模型. Wu(2004)就横向非均匀、自重、不可压缩黏弹地球模型,提出了耦合位扰动Laplace方程的三维有限元方法用于求解三维模型中的负荷问题,而最近Wong和Wu(2019)提出了迭代体力法突破了不可压缩性的限制. Wang和Wu(2006a, 2006b)通过实际算例详细讨论了横向非均匀性对冰后回弹和海平面变化的影响,指出冰后回弹的GNSS观测可以体现地幔的横向非均匀性.

      三维地球的变形研究,归根结底是要提出基于三维地球模型中的变形计算理论和方法. 然而,目前尚无一般化的理论能考虑内部分层界面的不规则起伏、局部构造、地表的地形效应等复杂模型中的同震和震后变形. 从一维到三维地球模型,随着地球模型越来越复杂,其计算过程也越发复杂. 目前由地震学层析成像等方法构建的横向不均匀性的三维模型不断更新,它们为三维地震变形理论提供了一个适当的地球模型,但基于这些模型的黏弹性固体潮变形和地震变形等变形理论和方法却有待进一步发展.

    • 一个多世纪以来,基于层状半无限空间和层状球形黏弹地球模型的变形理论不断发展完善,针对不同地球模型的几何形状和物理特性,涌现了诸多理论分支. 在半无限空间和均质球等简单地球模型中存在地震变形的全解析解,它可以作为其他数值计算的验证标准,并可以揭示地球模型的参数和地震变形场的直接关系,对深入认识地震变形机理有重要意义. 对于分层半无限空间和球形地球模型,借助一致性原理、基函数展开、积分变换等有效技巧,地震变形问题最终被转化为在波数域或者球谐域内的常微分方程组,之后通过数值方法求解复域内的边值问题. 对不同的地球物理模型,科学家们提出了不同的逆Laplace变换求解策略,如简正模叠加方法(Gilbert, 1971; Cambiotti et al., 2010b; Cambiotti et al., 2013)、Post-Widder方法(Spada and Boschi, 2006; Melini et al., 2008; Cannelli et al., 2009)、回路积分方法(Tanaka et al., 2006; Tanaka et al., 2007)和纵轴积分方法(Tang and Sun, 2019; Tang et al., 2020a; Tang et al., 2020b)等从复域内恢复时域内的变形. 总之,基于规则几何地球模型的黏弹地球变形理论已经趋于成熟,目前的理论可以同时考虑流变介质的自重、可压缩性、连续性径向结构等物理因素,准确模拟地震产生的位移、重力、应变等物理场,为地震变形的物理解释和地震循环过程的研究提供了理论基础.

      对于不规则三维地球模型,目前除了有限元等纯数值计算方法,尚无一般化的理论能考虑内部分层界面的不规则起伏、局部构造、不规则地形等复杂地球模型中的地震变形. 但已经有较多的研究者(Hu,2004汪汉胜等,2006Tanaka et al.,2009Gharti et al.,2019aGharti et al.,2019bLanger et al.,2019)用有限元等纯数值方法来模拟地震变形或负荷变形,并应用于实际的观测解释中. 有限元等数值方法的确可以模拟三维地球模型中的地震变形,但它一般需要大规模高性能计算机,普通用户很难实际进行. 纯数值的地震变形计算方法,超出了本文的研究范围,在此我们不做更多讨论.

      总之,规则几何形状的地球模型中的黏弹性地球变形理论,已经趋于成熟. 三维横向非均匀地球的准静态变形的理论的发展是非常困难的,但或许可以借鉴非规则分层的半无限空间地球模型中求解理论地震图的思路,如扩展广义反射、透射矩阵方法,之后引入介质的自重效应考虑重力对同震变形的影响,得到顾及地形、非均匀性和重力的同震—震后变形计算方法. 将层状球形地球模型的地震变形理论发展为三维地球的准静态变形理论,将是一个具有挑战的科学课题,也是位错理论发展的方向.

      地震变形理论,特别是黏弹性地震变形理论,在大地测量学和地球物理学中有广泛的实际应用. 例如,利用地震变形理论,可以理论上计算出大地震造成的地壳变形和重力场变化,从而来解释GNSS、InSAR、地表重力仪和空间重力卫星等大地测量数据(Imanishi et al., 2004; Zhou et al., 2012; Zhang et al., 2016; Feng and Bie, 2018; Liu et al., 2018). 此外,以地震位错理论为桥梁,结合这些变形观测数据,可以反演地震断层的破裂过程、滑动模型,推断地下介质的介质参数特别是上地幔的黏度结构(Zhou et al.,2014Yang et al.,2015b单新建等,2017Cambiotti,2020). 鉴于本文主题是介绍黏弹地球地震变形理论的现状和进展,我们对该理论的实际应用部分不做过多讨论,感兴趣的读者可以进一步参考相关论文和专著(孙文科,2012a2012b王启欣等,2015).

      致谢

      非常感谢两位匿名审稿专家的建设性意见以及编辑部对本文的仔细校对和编辑.

参考文献 (144)

目录

    /

    返回文章
    返回