自重可压缩层状球体地球模型下地震变形算法的优化研究

刘泰， 唐河， 佘雅文， 付广裕， 孟国杰

1 中国地震局地震预测研究所， 北京 100036 2 中国科学院大学计算地球动力学重点实验室， 北京 100049 3 中国地质大学(北京)地球物理与信息技术学院， 北京 100083

0 引言

随着大地测量技术的发展，GNSS(Global Navigation Satellite System)和GRACE(Gravity Recovery And Climate Experiment)等手段观测到的同震及震后变形数据日益丰富，为研究地震破裂过程和震后变形机制提供了数据支持(刁法启等, 2012; 刘泰等, 2017; 梁明等, 2018; Wang et al., 2022).其中，近场地震变形研究能够加深人们对地震震后余滑时空分布和震源区地幔黏滞性构造横向特征的认识(Diao et al., 2014; Hu et al., 2016)，而引入远场地震变形数据能够约束地震震级和深部地幔黏滞性结构(陈飞等, 2020; Wang et al., 2014).地震位错理论是研究地震变形的重要理论工具，其基于的地球模型直接影响地震变形的计算精度，如地球模型分层结构和不连续性影响近场地震变形计算(Dong et al., 2014, 2021)，地球曲率会影响远场地震变形计算(Dong et al., 2016)，地球自重效应会影响震后变形计算(Pollitz, 1997)，长时间震后变形与地球的可压缩相关(Tanaka et al., 2006).可见，基于接近真实的地球模型，完善地震变形计算方法具有重要的理论意义和应用价值.

数十年来，学者们基于不同几何特征的地球模型发展了多种地震变形计算方法.例如，Okada(1985, 1992)通过总结整理前人的工作，给出了一套完整、简洁、实用的同震形变计算公式，成为均匀半无限空间地震位错理论的经典表达式.Wang等(2003, 2006)发展了半无限空间模型的地震位错理论，考虑了地球的层状结构，提高了近场地震变形的计算精度.Sun(1992)和Pollitz(1996)发展了层状球对称地球模型的地震位错理论，考虑了地球曲率，提高了远场地震变形的计算精度.

相比于外部力源，如负荷、潮汐和剪切，产生的地球变形问题，地震变形问题略微复杂，前者需要求解齐次微分方程组，后者需要求解非齐次微分方程组.地震变形的微分方程组含有力源项(Takeuchi and Saito, 1972)，为非齐次问题.Okubo(1993)提出了地球内部点源位错在地表产生的变形与外部力源在震源处产生的形变存在互换关系(互换定理)，从而简化了地表地震变形的计算过程.当研究地震引起的地球内部变形时，微分方程组的解需满足震源处不连续和地表自由边界条件.为此，Takeuchi和Saito(1972)和Sun(1992)提出了两次数值积分、组合通解和特解的求解方法：从地心积分至地表得到通解，由震源处积分至地表得到特解，再利用地表边界条件确定通解的组合系数，进而得到微分方程组解矢量.该求解方法逐步成为了求解地震变形的常用方法之一.

当观测点靠近震源时，地球内部变形的格林函数需要通过高阶解矢量来得到.但对于高阶微分方程组求解，上述方法在数值积分的过程中会出现数值溢出和精度丢失现象.数值溢出问题是由向上积分过程中积分区间过长引起，出现在负荷变形和地震变形的计算中.最初提出的解决方法是调整积分起始半径，使其随阶次的增加而增大(Sun and Okubo, 1993; Martens et al., 2019).更为有效的方法是从数值放大因子(rn，这里n表示阶数，r表示半径)出发，在积分过程中对该因子进行分离(汪汉胜等, 1996; Liu et al., 2020).精度丢失现象是由于震源上方通解和特解组合引起，在微分方程组的高阶求解中尤为显著.一个直接的解决方法是使用第三方扩展函数库来进行高存储精度的计算，该方法已在黏弹性变形计算中得到应用(Melini et al., 2008).但使用该方法编写的计算程序依赖于第三方数据库，给用户带来了不便.更为有效的方法是调整数值积分策略、改变积分路径，对地震变形微分方程的齐次部分进行两次相向积分，即由地心积分至震源和由地表积分至震源，再利用震源间断条件求解震源上方和下方通解组合系数(Pan, 1997; Wang, 1999; Fukahata and Matsu′ura, 2005; Zhou et al., 2021).

为此，本文使用两次相向积分的方法来求解地震变形的微分方程组，以解决高阶微分方程组求解的数值不稳定问题，从而得到微分方程组解矢量.首先介绍了微分方程组的两种求解方法，接着与均质地球模型解析解进行比较来验证方法正确性，最后研究层状地球模型下震源附近地震变形.

1 地震变形理论

1.1 地震变形的基本方程

假设地球内部r=rs处存在点源位错，如图1a所示，其激发的位移场、应力场和引力位变化同时满足平衡方程、应力-应变关系和泊松方程(Takeuchi and Saito, 1972).利用球谐函数可以将以上物理量展开，整理可得两组微分方程组，分别对应球型变形和环型变形.平衡方程最初的非齐次项是体力，而位错源的等效体力包含了δ函数和δ函数的导数项.由于点位错源也可以用双力偶表达，Takeuchi和Saito(1972)利用震源处位移和应力间断特征推导了震源函数表达式，实际计算中只要把上述体力用双力偶代替即可.

dYs/dr=AsYs+δ(r-rs)ss,

(1)

dYt/dr=AtYt+δ(r-rs)st,

(2)

如图1b所示，断层上下盘相对运动矢量为Uv，坐标轴e1和e2分别表示纬度φ=0°和φ=90°的方向，坐标轴e3对应极轴方向，位错滑动矢量v和断层面法线矢量n可以表示为

(3)

式中v1,v2,v3和n1,n2,n3是v和n在球坐标系下的三个分量，即vinj(i,j=1,2,3)，为了方便后面位错源的表示，用角标ij代替vinj.Sun(1992)定义了四种独立震源，即垂直断层水平走滑震源(ij=12)、垂直断层上下倾滑震源(ij=32)、垂直断层水平引张震源(ij=22)和水平断层垂直引张震源(ij=33)，如图1c所示.通过以上四种独立震源的线性组合可以得到任意类型的位错源引起的变形.

图1 球坐标系下位错模型 (a,b)； (rs,0,0)为震源所在位置，(r,θ,φ)为球内任意一点，断层面的倾角为δ，滑动角为λ，上下盘断层面的相对错动量为U，n和v分别为断层面法线矢量和位错滑动矢量； (c) 表示四种独立震源

1.2 数值积分求解微分方程组

Sun(1992)提出了球体地球模型下点源位错引起的地表形变的计算方法，通过两次积分求解非齐次微分方程组，如图2a 所示，本文将该方法称为两次同向积分法.以球型变形为例，此时该非齐次方程组的解矢量分为两部分：

图2 两种数值积分路径示意图

Y=X·β+Z,

(4)

式中，X为对应齐次方程组的通解矩阵(6×3)，Z为非齐次方程组的特解矩阵(6×1)，β为三组通解的组合系数.其中特解部分由地心附近半径很小的球面(r=r0)积分至地表(r=a)，Love(1911)最早给出了通解部分的积分初值，特解部分由震源函数从震源处积分至地表得到，以满足震源处间断条件.当变形深度在震源上方时，解矢量由通解和特解组成，而震源下方的解矢量完全由通解组成：

(5)

满足地表边界条件，即地球表面应力分量为零和引力位连续，以及引力位导数穿过单层的不连续性(Sun and Sjöberg, 1999; Sun, 1992)，有y2(a)=y4(a)=y6(a)=0，可以得到三组通解的组合系数βi(i=1,2,3).然后带入公式(5)，得到任意深度处的解矢量(yj(r),j=1,2,…,6).

由于通解部分积分区间较长，在求解高阶微分方程组时，为了避免积分过程中出现数值溢出，积分起始半径需随着阶次的增加而增加(Sun, 1992).Liu等(2020)将该地表变形理论推广至地球内部，利用正规化方法避免地球内部变形计算过程中数值溢出的问题.但是当观测点位于震源上方时，解矢量由通解和特解组合得到，阶次足够高时会出现精度丢失现象，本研究通过调整积分路径解决该问题.具体操作是对微分方程的齐次部分进行两次积分，由地心积分至震源和由地表积分至震源(图2b)，本文将该方法称为两次相向积分法.此时解矢量完全由通解组成：

(6)

(7)

(8)

满足震源处间断条件后(Y2(rs)-Y1(rs)=s)，同时得到震源下方和上方通解的组合系数αi和βi(i=1,2,3).

(9)

可以发现本研究两次相向积分法和传统的两次同向积分法的差别在于处理震源间断和地表边界条件的先后顺序，优先处理地表自由边界可以避免微分方程组非齐次部分的求解，微分方程组解矢量完全由通解组成，此时震源上方和震源下方具有不同的组合系数.另外，两者方法在震源下方的积分过程完全一致，孙文科(2012)中详细介绍了内外核边界和核幔边界的处理过程，这里不在赘述.

2 地震变形计算结果

基于上述基本解法，给定不同的源函数(ss,st)求解通解组合系数，可得到对应震源的解矢量(Ys,Yt).本研究将基于Sun(1992)定义的四种独立震源进行求解，地球分层介质参数选自PREM地球模型(Dziewonski and Anderson 1981)，如图3所示，重力值由密度对地球半径积分得到.

图3 地球模型参数随深度的变化

本研究将从两个方面展示微分方程组的求解结果：(i) 忽略自重的均质球体地球模型下，与解析解(Dong et al., 2021)比较，验证本研究求解方法的正确性；(ii) 层状球体地球模型下，与两次同向积分法(Sun, 1992; Liu et al., 2020)比较，展示本研究求解方法的优势.

2.1 均质地球模型下解矢量

对于均质地球模型和层状地球模型的积分过程完全一致，只是每个子层上的介质参数不同，为此，本研究只需验证均质地球模型下求解的正确性，即可保证层状地球模型下的正确性，均质地球模型的介质参数设置为PREM模型的顶层参数值(λ0=34.2 GPa,μ0=26.6 GPa,ρ0=2.6 cm3·g-1).忽略自重的均质地球模型，微分方程组存在解析解(Dong et al., 2021)，为了与该解析解比较，本研究将均质地球模型的密度减少为原来的1%，以达到忽略自重的效果.

图4展示了均质地球模型下四种独立震源的解矢量(ds=20 km;dh=10 km，ds为震源深度，dh为变形观测深度).由于忽略自重的均质地球模型下解析解没有y5和y6，并且环型变形不受自重的影响，这里仅仅展示了变形解y1,y2,y3,y4.可以发现当密度减少为原来的1%时，数值积分解与解析解基本一致，验证了本研究求解方法的正确性.这里与Liu等(2020)的无自重处理方式不同，后者只是简单的将模型重力值减少为1%，没有改变密度值，可见通过减小模型的密度值以达到忽略自重的效果更为合理.另外，水平走滑震源受自重的影响较大，垂直倾滑受自重的影响较小，而自重主要影响垂直引张震源和水平引张震源的解y4，其中自重对垂直引张震源的影响超过100阶.整体而言，自重主要影响微分方程组的低阶求解，随着阶次增加，数值积分解与解析解趋于一致.

图4 均质地球模型下四种独立震源的解矢量(ds=20 km;dh=10 km)

2.2 层状地球模型下解矢量

这里主要比较了1.2部分中两种积分路径下的解矢量.介质参数来源于真实分层地球模型(图3)，由于地核主要影响10阶以内的微分方程组求解(Liu et al., 2020).在保证精度前提下节省计算时间，本研究将1000阶以上的积分起始半径设置为核幔边界，由于在每层积分过程中采用正规化处理，高阶微分方程组的求解不会出现数值溢出现象.

图5展示了层状地球模型下水平走滑震源的解y1(ds=20 km;dh=0 km,5 km,10 km,15 km).图5a1中虚线表示该震源深度ds和变形深度dh下解的截断阶次(Nmax=10a/|ds-dh|)，可以发现截断阶次以内，两种积分路径下的求解结果基本一致，超过截断阶次后，传统的求解方法出现数值不稳定现象.这是由于随着计算阶次的增加，传统求解方法中通解(图c1中实线)和特解(图c1中虚线)的数值量级逐渐增加，但两者符号相反，绝对值逐渐接近.在计算机的正常存储精度下，当通解和特解的有效数字完全一样时，组合得到的解矢量精度丢失，呈现数值震荡现象(图b1中红色圆点)，图(b1)和(c1)的数值计算结果展示在表1中.另外可以发现超过截断阶次后，该现象愈加明显.这种数值不稳定现象是组合通解和特解过程中精度丢失引起的(公式(5)).前面也已经介绍通过第三方扩展函数库可以进行高存储精度的计算，但该方法增加了通解和特解的有效数字，能够一定程度上提高计算阶次，但无法从根本上解决该问题.但本研究通过调整积分路径，震源上方和下方的解矢量均由通解组成(公式(6)和(7))，从根本上解决了该精度丢失的问题，得到的微分方程组的解矢量比较稳定，如图5b1中黑色曲线所示.

图5 层状地球模型下水平走滑震源的解y1(ds=20 km)

表1 层状地球模型下水平走滑震源的解y1(ds=20 km;dh=0 km)

3 震源附近地震变形研究

地震引起的地表变形研究中，渐近解结果显示层状地球模型和均质地球模型下的高阶解矢量基本一致，说明了浅源地震引起的地表变形受层状结构的影响很小.因此，可以利用均质地球模型下地表破裂源引起的近场地表变形对层状地球模型进行替代(Sun and Dong, 2013).本研究将进一步探讨层状结构对地球内部地震变形的影响，计算了均质地球模型和层状地球模型下微分方程组解的比值：

(10)

其中，yinhomo.为层状地球模型下的变形解，yhomo.为均质地球模型下的变形解，均质球参数来源于分层模型顶层(λ0=34.2 GPa,μ0=26.6 GPa,ρ0=2.6 cm3·g-1)，可以发现随着阶次的增加，两者的比值趋于常数，如图6中红线所示.另外，当均质球的介质参数设置为震源附近参数(λ0=45.9 GPa,μ0=44.1 GPa,ρ0=2.9 cm3·g-1)时，两种地球模型下的高阶变形解比值趋于1，如图6中蓝线所示，此时近场地震变形受层状结构的影响较小.这里只展示了走滑震源的计算结果，该特征对其余三个独立震源也适用，为节省篇幅，不逐一展示.

图6 均质地球模型和层状地球模型下走滑震源变形解的比值，纵坐标为比值的绝对值(ds=20 km;dh=19 km)

该部分将主要关注震源附近的地震变形计算，当变形计算深度接近震源深度时，需要对更高阶次的解矢量求和，甚至需要计算无穷阶次的解矢量，数值解无法实现该计算要求.Liu等(2020)计算了20 km震源在20 km处引起的变形，当对105的解矢量进行求和时，0.03°以内的格林函数不收敛.Takagi和Okubo(2017)推导了均质地球模型下地球内部变形的渐近解，解决了均质地球模型下震源附近地震变形不收敛的问题.本研究将进一步解决层状地球模型下震源附近地震变形不收敛的问题.

具体思路是对两种地球模型下解矢量的差异进行求和得到层状结构的影响，再与均质地球模型的渐近解叠加，从而得到层状地球模型下格林函数.例如，计算层状地球模型下垂直位移格林函数时，可以做一个变换，即加上和减去一个均质地球模型下的变形解，而等式关系保持不变：

(12)

图7展示了层状结构对水平走滑震源附近变形的影响，均质球的介质参数设置为震源附近参数.该结果对应公式(11)的第二项，可以发现ur,er θ,er φ分量的层状结构影响趋于稳定，而对于其余6个分量，越靠近震源，受层状结构的影响越小.将其与对应分量的渐近解叠加即可得到层状地球模型下格林函数.

图7 层状结构对水平走滑震源的影响(ds=20 km)，纵轴为均质地球模型和层状地球模型下格林函数的差异，不同线型代表不同的变形深度

4 结论

本文基于自重、可压缩和连续分层球体地球模型，利用龙格库塔数值积分方法求解微分方程组.对传统求解方法的积分路径进行调整，对微分方程的齐次部分进行两次相向积分，由地心积分至震源和由地表积分至震源，利用震源间断条件求解震源上方和震源下方通解的组合系数，避免高阶微分方程组求解中精度丢失的问题.

通过与均质地球模型下微分方程组的解析解比较，验证了求解方法的正确性.并发现严格的无自重处理方式是缩小模型的密度值，而不是简单的缩小重力值.另外，水平走滑震源受自重的影响较大，垂直倾滑受自重的影响较小，而自重主要影响垂直引张震源和水平引张震源的变形解y4，其中自重对垂直引张震源的影响超过100阶.

最后比较了均质地球模型和层状地球模型下高阶解矢量，发现两者的比值趋于稳定，当均质地球模型参数设置为震源所在层参数时，并且变形深度和震源处于模型同一层，两者的比值趋于1.即与层状结构对浅源地震地表变形的影响类似，对于地球内部地震变形，层状结构对震源附近变形的影响较小.基于该发现，结合均质地球模型下渐近解，可解决分层地球模型下震源附近地震变形发散的问题.

致谢感谢中国科学院大学孙文科教授提供的球体位错理论配套程序.感谢两位匿名审稿专家对本文的修改和完善提出了宝贵的意见.