地震发生在断层上，断层的破裂会辐射地震波，并引起地壳变形和地球物质迁移。震源力可以用位错/位移间断来描述。位错理论是现代地震学的重要理论支柱，是模拟同震变形、震后变形以及地震波形的重要物理工具。通过位错理论计算格林函数，并利用大地测量或地震学观测数据反演震源参数和滑移是研究地震破裂、构造地质等方面的重要手段^[1-2]。随着重力卫星GRACE能够观测到全球特大地震，位错理论模型在解释卫星观测信号以及利用同震重力变化确定震源参数的研究中也起着重要作用^[3-7]。

自Steketee^[8]将位错理论引入地震学后，大量学者研究了均匀半无限空间地球模型下的同震变形问题。Okada^[9]总结整理了前人的工作，并给出了点源和矩形断层的解析公式，可以用于计算任意剪切或引张断层在均匀半空间下引起的同震地表位移、应变和倾斜，解析公式简洁实用、计算效率高，广泛应用于地震断层参数和滑移分布的反演问题。陈运泰等^[10-11]讨论了结合半无限空间位错理论和地表形变资料进行震源反演的一般方法，并利用实际观测数据反演了1966年邢台地震和1976年唐山地震的震源过程。Rundle^[12]利用贝塞尔函数得到了无自重的层状半无限空间的同震变形。Okubo^[13]推导并给出了弹性均匀半无限空间的同震重力变化，与Okada相似，该理论模型也是简洁的解析公式。基于Takeuchi等^[14]给出的地震球面波理论，Sun等^[15-20]利用龙格-库塔数值方法解决了在一个球对称的分层地球模型下，由内部位错源引起的地表位移、应变、重力和垂线偏差的变化。研究表明，当震源距超过10°，地球曲率的影响可达10%；震源距大于0.5°时，需要考虑地球径向的不均匀性^[18]。Okubo^[21]提出了互易定理，将位错解与潮汐、剪切力、负荷变形的解联系在一起，即一个内部位错源引起的地表变形可以用在该深度的潮汐解、剪切力解和负荷解线性表示。Pollitz^[22]给出了无自重层状球形地球的同震位移。基于Rundle和Okada的理论模型，Wang等^[23]发展了分层半无限空间的位错理论，并采用传递矩阵的方法得到地表的同震变形。郝金来等^[24]使用了梯形积分和Filon积分相结合的方法得到了分层介质地球模型的同震位移、应力和应变，加快了计算速度。Fu等^[25]在Sun等^[15]的球体位错理论基础上，利用微扰原理研究了球形地球横向不均匀性对同震变形的影响。

将弹性位错理论进行扩展可以得到震后黏弹性变形。对黏弹性本构关系经过Laplace变换后得到与弹性问题相同的形式，按照弹性位错方程的解法得到Laplace域的变形，再对其进行Laplace反变换得到震后黏弹性变形。Rundle^[26]考虑了一个更接近真实地球的自重、分层半无限空间模型，给出了矩形推覆断层引起的黏弹性变形。Pollitz^[27]解决了在一个黏弹性、无自重的球形地球模型下，由位错源引起的区域同震位移场和应变场。Piersanti等^[28]、Sabadini等^[29]给出了自重、不可压缩、层状球形地球下，位错源引起的地表震后黏弹性变形。另外，多位学者给出了计算有限分层数的震后黏弹变形方法^[30-31]。Tanaka等^[32]利用复平面路径积分的方法解决了有限层数问题。Wang等^[33]利用虚轴积分给出了自重、可压缩、分层半无限空间的黏弹性变形，并给出了一套实用的计算程序。

径向一维分层地球模型下的弹性、黏弹性位错理论都较成熟，但以上研究工作都基于前人给出的经典理论方程的解，并未重新给出基本解。为了方便相关领域研究人员理解和使用球体位错理论模型，本文给出不同物理前提条件下的均匀球变形方程的基本解。

1 物理方程

对于无自重的、非自转的地球模型，平衡方程可以写为：

(1)

式中，ρ为介质密度，σ为应力张量，F为外力。如果只考虑弹性地球介质，根据胡克定律，应力应变关系为：

(2)

式中，u为位移矢量，λ和μ为拉梅常数，后者也是剪切模量或刚度，I为单位矩阵，。假设介质为均匀的、各向同性的，则有：

(3)

将式(3)代入式(1)中，并运用：

(4)

可得：

(5)

考虑地球的自重，假设无变形的地球是球对称的，并处于静水平衡状态，初始压力可以写为：

(6)

令r₀和r表示变形体的初始位置和变形后的位置，则与位移有如下关系：

(7)

且初始平衡应力为：

(8)

其散度为：

(9)

另一方面，变形后的密度为：

(10)

故密度增量为：

(11)

变形体的引力位为初始引力位和扰动位之和：

(12)

初始引力位的梯度为正常重力，即。变形单元因重力作用产生的体力增量为：

(13)

将式(9)、(13)代入式(1)中，并考虑，可得自重的平衡方程：

(14)

式中，σ^δ为应力变化，F_D为位错力源。考虑到并不引起混淆，此后用σ替代σ^δ。将式(9)、(13)代入式(5)中，得到平衡方程的另一种形式：

(15)

另外，在球体内部的扰动位满足泊松方程:

(16)

2 球形地球下的常微分方程组

本文考虑的地球模型是均匀的、非自转的和球对称的模型，故需要考虑§1得到的平衡方程(15)和泊松方程(16)在球坐标(r, θ, λ)下的解。首先，将位移、应力和引力位表达为球谐函数形式。对于矢量场，可以分解为球型和环型的多极表达，位移和法向应力矢量可以写为：

(17)

式中，

(18)

式中，是表面梯度算子，Y_lm(θ, λ)=P_lm(cosθ)e^imλ是面谐函数，P_lm(cosθ)为勒让德函数。引力位可以按照标量场展开：

(19)

为了得到方程的解，除了上述未知变量(y₁(r), y₂(r)~y₅(r))之外，还需引入另一个新变量^[14]：

(20)

引力位在地球内部满足泊松方程(16)，在地表满足拉普拉斯方程。地表处的边界条件较为复杂，而y₆(r)在地表为0，该变量的引入将对后面利用边界条件得到方程的通解系数提供方便。同引力位的展开，将体积变化作球谐展开得：

(21)

且χ_l和y₁、y₃有如下关系：

(22)

平衡方程(14)可以分为两个分别关于e_r和e_θ的方程，略去繁琐的数学过程，直接给出这两个方程的化简结果：

(23)

(24)

同理，弹性本构关系方程(2)在球坐标下表示为：

(25)

式中，δ_ij为狄拉克符号，应变分量的球坐标表达为：

(26)

同样，略去数学过程，得到以下两个微分方程：

(27)

(28)

由泊松方程可得一个二阶微分方程：

(29)

联合式(20)~(24)以及式(27)~(29)，并考虑可压缩的(χ_l≠0)一般形式，可得关于y₁(r)~y₆(r)的一阶齐次微分方程组：

(30)

式中，记。式(30)是6个球型的一阶微分方程。相对而言，可以较容易地得到2个关于环型解和的微分方程组：

(31)

需要指出的是，式(30)和(31)给出的都是齐次方程组，不包括位错力项。对于一个非齐次方程组，其解的构成是齐次方程组的通解与非齐次项的特解。非齐次项也可以通过边界条件解决，即位错变形的解可以通过方程组(30)和(31)的通解在地表和位错间断面满足边界条件得到，故这里主要考察齐次方程组的解。

3 均质球下的基本解

本节将在不同的物理简化假设下给出变形方程在均质球的基本解，首先讨论不同地球模型的球型解，最后给出环型解。

3.1 无自重的地球模型

首先，假设地球模型无自重，平衡方程中的G、ψ和g₀均消失，同时泊松方程消失。式(30)中关于y₅和y₆项消失后，6个球型一阶微分方程退化为4个:

(32)

通过数学运算，容易得到上述方程组的解：

(33)

式中，f₁=(l-2)λ+(l-4)μ，f₂=(l+3)λ+(l+5)μ。不考虑自重情况下的解为r^l和r^-l的函数，式(33)中矩阵的逆可以解析得到。通过传递矩阵方法将球心处的解和地表解分别传递到位错间断面的下、上表面，满足位错简单条件和地表处自由表面边界条件得到通解的系数。鉴于解的便利性，该模型在计算同震位移时应用广泛。

3.2 不可压缩的地球模型

在短时间尺度上，地球的体积因外力作用发生变化；在长时间尺度上，变形地球可以近似认为体积不变，即不可压缩介质。不可压缩的地球模型是模拟冰川均衡调整、长时间尺度的震后黏弹性变形等地球物理问题的经典模型。当体积不变时，χ_l=0，密度为常数∂_rρ₀=0，故由式(11)可知密度增量为0，即Δρ=0。此时，平衡方程(14)和泊松方程(16)解耦，且泊松方程退化为拉普拉斯方程。仅考虑齐次方程形式，有：

(34)

(35)

对于不可压缩介质，令

(36)

故在不可压缩条件下，本构关系式(3)为：

(37)

将式(37)代入式(34)中得到：

(38)

式中,

(39)

(40)

将Γ和H在球坐标系展开：

(41)

(42)

将式(41)、(42)作用在式(38)得：

(43)

一方面，对式(38)取散度，并使用=0有：

(44)

式中，，式(44)为拉普拉斯方程，解为：

(45)

同理，方程(35)的解为：

(46)

顾及χ_l=0，由式(22)可得：

(47)

另一方面，结合式(40)可得：

(48)

联合式(43)、(45)和(48)可得y₁的解：

(49)

代入式(47)可得y₃的解：

(50)

将y₁和y₃代入式(28)可得到y₄：

(51)

由应力-应变关系得：

(52)

将式(45)和(49)代入可得y₂的解：

(53)

最后得到y₆：

(54)

至此，得到了自重、不可压缩、球对称的均质地球模型的解，式中(C₁, C₂, C₃, )为待定系数，可以通过地表自由边界条件和位错间断条件求得。

3.3 自重、可压缩的地球模型

当自重效应和压缩性无法忽略时，平衡方程(15)和泊松方程(16)耦合，且可以化为6个一阶微分方程组(30)，其解较无自重模型和不可压缩模型复杂。同上，在均质球下，密度为常量，即∂_rρ₀=0。对平衡方程(15)的齐次方程两边取旋度，有：

(55)

式中，H定义见式(40)。另一方面，对方程(15)两边取散度，并整理后有：

(56)

将球谐展开式代入式(55)、(56)，得:

(57)

(58)

式中，，分别是P波和S波速度；。对式(58)两边作用算子，有：

(59)

式中，。令

(60)

则式(59)可以写为：

(61)

该方程可以化为2个Helmholtz方程，得到4个独立的解：

(62)

式中，j_l(kr)、n_l(kr)(k=k_c±)分别是第一类和第二类球贝塞尔函数，j_l(kr)在r=0处正则，在r=∞处奇异；n_l(kr)与之相反。不失一般性，下面将j_l(kr)和n_l(kr)统一用f_l(kr)表示，求出自重、可压缩的均质球变形方程的基本解。设：

(63)

将其和式(62)一并代入到方程(57)中，可以求得：

(64)

另一方面，顾及χ_l和H_l的定义式(22)以及

(65)

可得：

(66)

(67)

进一步，利用关系rf′_l(kr)=lf_l(kr)-krf_l+1(kr)，式(66)、(67)写为：

(68)

(69)

假设ry_r=C₁f_l(kr)+C₂rf_l+1(kr)，并利用球贝塞尔函数的递推关系：

(70)

可得到y₁的解为：

(71)

式中，

(72)

相似地，可得y₃的解为:

(73)

再将y₁和y₃代入式(27)、(28)中，得到y₂和y₄：

(74)

(75)

式(74)中，β=λ+2μ，与前文定义一致。

泊松方程可以写为：

(76)

可以求得y₅的解：

(77)

最后，由式(20)得到y₆的解：

(78)

由于k=k_c±以及f_l(kr)=(j_l(kr), n_l(kr))，上述过程可以得到4组关于y₁、y₂、y₃、y₄、y₅、y₆的解。另外，可以容易得到二组解：

(79)

3.4 环型解

相对于复杂的球型方程，环型方程(31)较为简单。该方程没有重力项和体积变化项，故在任何物理假设条件下，其解均为：

(80)

4 结语

本文回顾了位错引起的地球变形和引力位变化的方程，并给出了非自转、球对称、各向同性的均质球下的基本解。

1) 对于无自重的地球模型，不考虑引力位的改变，6个一阶微分方程组退化为4个，其球型解为关于r^l和r^-l-1的函数。

2) 对于不可压缩的地球模型，体积不变，平衡方程和泊松方程解耦，其球型解同样为r^l和r^-l-1的函数。

3) 对于顾及自重和可压缩性的模型，其球型解除了2组分别关于r^l和r^-l-1的函数外，还有4组关于特殊函数j_l(kr)和n_l(kr)的解。

4) 环型解的形式为关于r^l和r^-l-1的函数。

在这些解中，r^l和j_l(kr)在球心处正则，无穷远处奇异；r^-l-1和n_l(kr)在球心处奇异，无穷远处正则。

以上3种一维地球模型各有用处，例如，鉴于无自重均质球模型的简易性，可以用于与半无限空间模型比较研究地球的曲率效应^[34]；不可压缩模型可以用于近似模拟长时间尺度的震后变形^[28]；自重、可压缩的球体模型可以用于模拟完整的同震和震后变形、重力场变化^{[19, 32]}。

在实际计算地球因内部位错引起的变形时，需要考虑分层效应。此时，微分方程组(30)和(32)中与地球材料有关的参数不再是常数，而是随着深度改变。无论如何，在某一分层上，这些参数仍然可以视为常数，满足微分方程组，故解的形式为本文给出的基本解。对于无自重和不可压缩的地球模型而言，基本解为简单的解析解，组成的矩阵的逆同样可以解析得到，故可以通过矩阵传递的方法，满足地表自由边界条件和位错间断条件，求出待定系数。对于自重可压缩的地球模型而言，基本解含有球贝塞尔函数，难以解析得到逆矩阵。可以通过龙格-库塔方法将微分方程组(30)的通解从球心处积分到地表，再将位错的特解从位错间断面积分到地表，共同满足地表自由边界条件得到待定系数。

本文着重解决微分方程组(30)的通解，位错特解可以用位移间断表达，其球谱域的展开表达式可以参考文献[14]和[35]。