球形地球模型的地震位错理论及其应用＊.docx

球形地球模型的地震位错理论及其应用*孙文科付广裕周新徐长仪唐河董杰周江存杨君妍王武星刘泰1）中国北京100049中国科学院大学地球与行星科学学院2）中国北京100083中国地质大学（北京）地球物理与信息技术学院3）中国北京100085应急管理部国家自然灾害防治研究院4）中国北京100029中国科学院地质与地球物理研究所5）中国北京100036中国测绘科学研究院6）中国武汉430071中国科学院精密测量科学与技术创新研究院7）中国北京100036中国地震局地震预测研究所1球形地球模型地震位错理论的发展地震位错理论是研究震源机制、反演断层、解释观测数据等的必要理论基础，在地震学、地球动力学、大地测量学中占有核心地位,地震位错理论需要基于一定的地球模型来建立，它相当于地震震源与地震变形之间的系统函数，从理论上给出地球在地震力源作用下的变形响应.由于位错理论的计算精度取决于地球模型，地球模型的准确与否将直接影响理论计算精度及其适用范围.所以，基于合理的地球模型来建立位错理论具有重要科学意义.为精确地计算同震和震后变形、合理地解释全球大地测量观测数据、更准确地反演断层滑动分布，本团队在研究中采用了球形地球模型,并据此建立了地震位错理论新体系,新位错理论包含地球曲率、层状结构效应、黏弹性效应、三维非均匀效应等，在物理上更严谨、更合理，同时适用范围更宽泛.而建立球形地球模型的位错理论需要采用与半无限空间位错理论不同的研究思路和数学方法，需要在球坐标系中求解地震变形的偏微分方程组，解决积分变换、级数收敛、积分稳定性等一系列数学问题.根据观测技术以及理论本身的发展，本研究团队先后针对一维分层弹性球形地球模型、三维不均匀弹性球形地球模型、一维分层黏弹性球形地球模型建立了不同的理论分支，形成了较完善的地震位错理论体系.1.1 弹性球形地球模型的地震位错理论球形地球模型的地震位错理论首先采用一维分层球形弹性地球模型，即球对称、不旋转、弹性和横向各向同性的地球模型(SNREI模型，Dahlen,1968).SNREI模型是一大类地球模型，其中1066A模型(Gilbert,Dziewonski,1975)、PREM模型(DZieWonski,Anderson,1981)和IASP91模型(Kennett,Engdahl,1991)等被广泛使用.基于1066A和PREM模型，孙文科等提出了弹性球形地球模型的地震位错理论,采用半解析半数值方法求解地震变形的微分方程组，定义了地震位错勒夫数，克服了勒夫数计算、求和、截断和加速收敛等数值计算困难，给出了四个独立点震源的格林函数(SUn,1992,2022,2022a,b；Sun,Okubo,1993,1998；Sunetal,2022a,2022),并在此基础上进一步给出了计算任意有限断层在任意位置产生变形的数值积分算法(Sun,OkUbO,1998).利用该理论可以解释实际地震产生的位移、应变、重力位、重力等的同震变化，为地表同震变形解释、断层破裂反演等研究奠定了基础，目前也已被广泛使用.更重要的是，这些研究搭建了球形地球模型位错理论的理论框架，为后续发展三维地球变形和黏弹性地球变形理论奠定了基础.团队随后对地表破裂源、地震造成的地球动力学变化等问题进行了研究，扩展了地震位错理论的内涵及其适用范围.例如：孙文科和董杰借助互易定理(Okubo,1993)推导了地表破裂源产生的地表变形的地震位错勒夫数和格林函数，解决了地表震源的奇异解问题(Sun,Dong,2022)；徐长仪和周江存等给出了地震造成的地球体积变化、低阶重力场系数变化、地球自转变化和重力势能变化的计算公式(Zhouetal,2022,2022a,2022；Xu,Sun,2022；Xuetal,2022；Xu,Chao,2022；Xu,2022)；孙文科、徐长仪和周江存等给出了特大地震造成的地心移动解(SUrbDong,2022；Xu,Chao,2022；Zhouetal,2022).近几年周江存等提出了计算位错勒夫数的近似解析解方法，提高了计算效率，同时解决了内部格林函数的收敛问题(ZhOUetal,2022a,b,2022).迄今为止，上述弹性球形地球模型的地震位错理论已经得到了广泛的应用.利用地震变形格林函数，SUn和OkUbo(1998)成功地解释了1964年阿拉斯加MW9.2大地震的同震重力变化，其理论计算结果与观测的重力变化基本吻合，表明新理论可以合理地解释现代大地测量技术观测到的同震变形.此外，为了便于国内科研人员使用该理论，付广裕和孙文科(2022b)撰文详细地介绍了球体位错理论计算程序的总体设计与具体实现过程，并开始在科学社区分享其Fortran程序,得到了学界的认可.更重要的是，该理论因其物理模型的合理性而优于半无限空间位错理论，为地球物理学、大地测量学和地球动力学等应用提供了全新的、可靠的理论保证(孙文科，2022b).该理论已经成功应用于地震断层破裂反演、地震解释以及火山产生的变形，并且可以用来解释同震变化(Okubo,2022).例如,Imanishi等(2022)利用我们团队的球形地球位错理论合理解释了超导重力仪记录到的2022年日本十胜冲(Tokachi-Oki)MW8.3地震产生的亚微伽级重力变化，该研究结果发表在国际著名期刊SCienCe上，属于世界上首次报道亚微伽级同震重力变化.2022年孙文科系统地总结了弹性地震位错理论与部分应用的研究成果，撰写了专著地震位错理论(孙文科，2022a),这是目前国内外关于球形地球模型位错理论的唯一专著.1.2 三维不均匀弹性地球模型的地震位错理论随着现代大地测量观测精度的日益提高，真实三维地球结构对地震变形的影响逐步被观测到(HUetal,2022；Suito,Freymueller,2022；Freedetal,2022；Diaoetal,2022；Tianetal,2022),然而上述基于球对称地球模型的地震位错理论却无法体现局部的三维效应.为了提升地震变形的理论支撑作用、提高对地震变形观测数据的解析能力，在上述一维球形模型的地震位错理论基础上，付广裕等提出了三维不均匀地球模型的地震位错理论(Fu,Sun,2022,2022,2022；Fuetal,2022；付广裕，孙文科，2022a).基于微扰方法的三维地球模型的位错理论，提高了地震变形和重力潮汐的计算精度.具体地说，FU和SUn(2022)推导了密度的横向增量对重力潮汐影响的计算公式，弥补了Molodenskiy(1980)理论的不足，并计算了半日潮重力潮汐因子的全球分布，为三维地球模型的地震变形研究提供了参考.之后，就三维地球模型的地震变形问题，将实际三维不均匀地球分解成球对称地球模型和对应的横向不均匀增量，分别进行计算得到球对称解和三维响应，之后叠加得到最终变形(Fu,Sun,2022,2022,2022；Fuetal,2022).球对称解可以直接利用上述球对称地球模型位错理论(SUn,1992)计算得到.三维响应即地球的横向不均匀结构对同震物理场变化的影响可以分为震源的响应和地球横向不均匀结构的响应，可分别由震源函数的扰动及平衡方程式的变分来求解.经过理论推导，给出了六个特殊点源位错引起的位移、重力变化等计算公式(一个垂直走滑位错，两个相互垂直的倾滑位错，三个引张位错).对这些公式进行线性组合即可计算任意位置、任意类型位错产生的同震变化.依据36阶P波速度模型,FU和SUn(2022,2022)利用岩石试验经验关系式推导出三维S波速度模型和密度模型，并计算出三种典型类型的点源位错产生的同震重力变化，其结果显示三维响应与位错类型、震源深度均相关，其最大响应约占球对称解的0.5%,且S波速度模型的影响最大.数值结果同时表明：三维响应中震源的响应与地球横向不均匀结构的响应处于同一量级，地球横向不均匀结构对同震位移的贡献达到球对称解的1%-2%(Fuetal,2022；付广裕，孙文科，2022a).上述研究成果得到了国内外同行的普遍认可.发表在Science上的Ito和Simons(2022)一文将Fu和Sun(2022)的成果视为潮汐理论研究领域的重要进展加以介绍.国际大地测量联合会2022年度的Bomford奖获得者YoshiyukiTanaka博士在多篇论文中引用上述研究(Tanakaetal,2022,2022,2022),将其视为位错理论研究领域的前沿进展加以介绍,在孙文科教授的专著地震位错理论序言中，许厚泽院士用一句话将上述研究总结如下：”利用变分方法研究了极为困难的三维地球模型的同震变形问题”.这均说明三维不均匀地球模型的地震位错理论虽然难度很大，但是它在实际应用中具有重要意义和价值.1.3 黏弹地球模型的地震位错理论地球介质除具有弹性性质之外，其在较长时间尺度上还具有显著的流变特性(Takeuchi,Hasegawa,1965；Anderson,0,Connell,1967；汪汉胜等，1997；PeItier,2022)；地震变形除了呈现明显的同震变形外，还表现出显著的震后变形(NUr,Mavko,1974；Freed,Burgmann,2022；Wang,2022；Chenetal,2022).大地震的震后变形已经可以被现代大地测量技术清晰地观测到，例如，GNSS.InSAR和GRACE均可检测到震后信号(Panetetal,2022；Ryderetal,2022；Wangetal,2022；刘泰等，2022；Cambiotti,2022)为了解释震后变形观测数据，需要建立基于黏弹性地球模型的地震位错理论.为此，很多学者针对黏弹性地球模型加以研究，提出了相应的计算方法，例如简正模叠加、回路积分和POStTidder等方法(Pollitz,1997；Wang,1999；Tanakaetal,2022,2022；Melinietal,2022).然而，黏弹地球模型的地震位错理论的发展具有一定的挑战性，对于复杂的黏弹性结构地球模型，始终存在数值计算困难.例如：简正模方法无法处理地球的径向连续结构或介质的可压缩性；回路积分方法需要针对特定黏弹性模型，通过数值搜索方法确定积分路径；Post-Widder的实现需要扩展精度函数库.总之，简正模分布复杂、积分核数值震荡、收敛慢、计算量大等问题均未得到完美解决(唐河，孙文科，2022).为此，唐河和孙文科基于上述弹性位错理论的基本思想，发展和完善了基于黏弹地球模型的地震位错理论(Tang,Sun,2022；Tangetal,2022a,b).黏弹性变形理论的核心是在复数域计算复勒夫数，并进行逆拉普拉斯(LaPIaCe)变换以得到时域勒夫数,首先将黏弹性地震变形方程组进行拉普拉斯变换至复数域，转换为拉普拉斯域内的等效弹性问题，从而可以利用上述弹性位错理论的处理方案加以解决，得到复域的地震位错勒夫数和格林函数.该步骤所使用的计算方法已在弹性理论中得到充分验证，可以确保计算的准确性.之后，在将复域解变换回时间域解时，需要处理逆拉普拉斯数值积分，但积分核的震荡特性造成了严重的数值积分困难，极具挑战性.为此，我们对逆拉普拉斯积分核作分式函数近似和级数展开(Valsa,Brancik,1998；Tang,Sun,2022；Tangetal,2022a,b),得到一个交错级数，从而使得时域勒夫数的计算表达为沿虚轴采样的复勒夫数之和.本方法采用的逆拉普拉斯积分方法也已在电气工程等多个领域得到充分验证，具有很高的计算精度.这样，我们构建了一种新的计算黏弹性地震变形的虚轴积分方法，这是黏弹地震变形理论的最新方法.该方法对于任意线性黏弹性地球模型和任意力源类型均适用，且具有一致的计算精度，其计算方法实现简单、收敛快，不需要扩展精度函数库.因此，该方法非常适于黏弹性结构的反演研究.此外，本团队还发展了黏弹性均质球模型的震后变形解析解.例