0 引言

地震活动是一种重要的地球动力学现象.地震可以造成巨大的社会和经济损失，同时提供了人们对板块运动、板块内部变形、地球结构和断层活动等现象的理解.地震的发生使地球内部质量、应力和应变重新分布，造成地表位置和重力场的改变.上述现象均可以被现代大地测量技术所检测到，如重力卫星GRACE(Gravity Recovery and Climate Experiment)、全球导航卫星系统GNSS(Global Navigation Satellite System)(Khazaradze and Klotz, 2003; Anzidei et al., 2009; Han et al., 2010; Heki and Matsuo, 2010; Chen, 2019; Gautam et al., 2019).另一方面，根据弹性地球模型位错理论，将震源视为断层两侧的相对运动，同样可以得到同震位移、重力、倾斜和应变等变化.位错理论将地表同震变形和震源机制解之间建立起关系，并用来解释大地测量观测现象，对地震学的发展起到了重要的作用.

同震变形最早由Steketee (1958)提出，他将位错理论应用到固体物理学中解释晶体材料的强度.随后，许多地球物理学家研究了典型位错源的同震变形，包括垂直走滑断层和倾滑断层，同时将同震变形观测资料应用于断层滑动分布的反演过程中来(Chinnery, 1961, 1963; Maruyama, 1964; Press, 1965; 陈运泰等, 1975, 1979; Iwasaki and Sato, 1979; Davis, 1983).Okada(1985, 1992)基于均匀各向同性半无限空间模型，给出了一套完整的剪切位错和引张位错引起的地球表面和内部位移、应变和倾斜的解析表达式，同时考虑了点位错和有限断层位错.Meade (2007)将该解析解推广到有限三角形位错形变中.Okubo(1991, 1992)研究了半无限空间介质内的重力变化问题，同样给出了一套计算同震重力位和重力变化的公式.

事实上，地球并不是均匀的，它的材料组成和物理性质在径向上变化最快.另一方面，由于地球近似球形的结构，在计算远场(距离震中~100 km以外)同震变形时，曲率的影响可以达到25%以上(Sun and Okubo, 2002; Dong et al., 2014).Sun(1992)将位错理论从层状半无限空间地球模型扩展到球形层状地球模型，定义了4种独立点源的地震位错Love数，给出了精确计算SNREI地球模型下Love数和格林函数的数值方法，以及计算任意地震同震变形的有限断层积分算法(Sun and Okubo, 1993, 1998; Sun et al., 2009).考虑到地球内部三维不均匀弹性介质地球模型的影响，付广裕基于空间分辨率为5°×5°(对应于36阶球谐函数)的P波速度模型、S波速度模型和地球介质密度参数，计算了地球横向不均匀构造对2阶重力固体潮因子以及地震变形的影响(Fu, 2007; Fu and Sun, 2007).上述微扰理论只适用于较小区域的研究，对于大范围三维构造模型的地震变形问题不适合.

地震位错Love数是地球介质对震源函数响应的特征数，它的精确计算直接影响地球变形格林函数的精度，以至于影响地震断层反演的可靠性(Sun and Okubo, 1993).值得注意的是，任何计算地震变形的位错理论均基于一维的地球模型，无法反映出局部构造的影响.因而，在震源函数一定的情况下，地球模型将直接影响地震位错Love数、格林函数以及同震变形的计算结果(Okubo et al., 2002).Sun和Okubo(2004a, b)在研究中使用的位错Love数是根据球对称地球模型得到的，如1066A模型(Gilbert et al., 1975)或PREM模型(Dziewonski and Anderson, 1981).然而，这些地球模型是由全球地震波数据得到的，对于局部地区而言，这些模型与真实地球具有一定的差异.这意味着位错Love数和格林函数的精度直接依赖于所采用地球模型的精度.由于大地测量观测数据包含更真实的地球构造信息，可以利用该数据研究地震位错Love数以及格林函数.GPS和DInSAR能够提供精确的地表变形观测数据，尤其是GPS提供的水平方向和垂直方向的位移；GRACE重力卫星可以检测到地震引起的重力场变化，拟合得到同震重力阶跃(Sun and Zhou, 2012)，为地震位错Love数的反演提供了可能.因此，通过GRACE和GNSS技术确定的地震位错Love数以及格林函数更能反映出区域构造特征，在解释大地测量观测数据反演地震断层滑动分布的研究上具有重要的作用.

Sun等(2006b)研究了GRACE的同震变形特征，并提出了利用该数据反演地震位错Love数的基本思想，从理论上论证了该研究的正确性，并以2004年Sumatra地震(M_W9.3)为例作了数值模拟，保证该方法数据计算结果的精度.Yang等(2015)基于2011年Tohoku-Oki地震(M_W9.0)，将上述方法应用到实际观测中，并且考虑了海水质量重新分布的贡献，将GRACE观测到的同震变形与理论同震变形建立观测方程，从而反演得到体现局部构造特征的重力位位错Love数k.结果表明，利用反演得到的位错Love数计算的同震重力变化与GRACE观测值符合得更好.杨君妍等(2017)给出了远场GPS水平位移在位错Love数反演中的变形特征分析，指出远场位移可以用来约束4个独立震源中变化最大的量.

值得注意的是，上述研究仅开展了重力位位错Love数k的反演，事实上，位错Love数还包括垂直位移位错Love数h和水平位移位错Love数l.在实际应用中，我们需要获得一套完整的区域地震位错Love数以及相应的格林函数.本文首先介绍了球形地球模型位错理论的基本原理，包括定义位错Love数和4个独立位错源的格林函数，并通过数值积分给出同震变形的计算公式；其次给出利用大地测量观测数据反演位错Love数以及确定正则化因子的方法；最后以2011年日本东北地震(M_W9.0)为例，反演体现局部构造特征的地震位错Love数和格林函数，以验证该理论方法的可行性.

1 地震位错Love数及格林函数

考虑一个球对称、不旋转、弹性和各向同性的(SNREI)球形地球模型(Dahlen, 1968)，在球心距r_s处定义一个无限小的位错源dS位于坐标系(e₁, e₂, e₃)上，并且滑动矢量为v，法矢量为n，滑动角λ和倾角δ.其中坐标轴e₁和e₂分别表示赤道平面方向φ=0和φ=π/2，e₃沿着北极轴方向.U=(U/2)－(－U/2)定义为两个断层的相对错动.如图 1所示.

根据准静态位错理论，在地球任意观测点(a, θ, φ)的同震位移和重力变化为(Sun and Okubo, 1993; Sun et al., 2006a)：

(1)

(2)

其中，y_{1, m}^{n, ij}, y_{3, m}^{n, ij}, y_{5, m}^{n, ij}和y_{1, m}^{t, n, ij}表示位错理论基本微分方程的解：y_{1, m}^{n, ij}和y_{3, m}^{n, ij}为径向和水平位移的解；y_{5, m}^{n, ij}是引力位的解.a表示地球半径, θ和φ分别表示观测点的余纬和经度.由于i, j=1, 2, 3，它们的组合是9个参量.考虑到震源函数和断层几何模型的对称性，我们将9个震源解减为4个，分别表示为垂直断层水平走滑破裂震源(ij=12)、垂直断层上下倾滑破裂震源(ij=32)、垂直断层水平引张破裂震源(ij=22)和水平断层垂直引张破裂震源(ij=33).通过4个独立震源解的组合，我们可以计算断层在地球任意位置的同震变形.地震位错Love数的定义为如下：

(3)

(4)

(5)

(6)

其中，g₀是万有引力常数.地震位错Love数是无量纲的，表示弹性地球在内部作用力下的变形响应.位错Love数的大小与地震断层参数无关，跟位错理论所使用的地球模型相关.上述变化式可写为：

(7)

(8)

(9)

(10)

式(7)中，Y_n^m(θ, φ)=P_n^m(cosθ)e^imφ，P_n^m(cosθ)是非正规化的连带Legendre函数. 称为“地震位错因子”，用来描述地震的大小.

于是，4个独立解同震重力变化格林函数的表达式：

(11)

(12)

(13)

(14)

垂直位移格林函数的表达式：

(15)

(16)

(17)

(18)

水平位移格林函数的表达式：

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

一旦地震震源参数被确定，我们就可以通过上述4个独立点源的格林函数和地震位错因子的组合计算地球表面和内部的同震变形.另外，由于震源位于北极，m只需取0~2即可.滑动矢量v和法矢量n可由滑动角λ和倾角δ表示.对于一个剪切型震源：

(27)

(28)

根据断层的几何参数，任意剪切型位错源的同震重力位变化、垂直位移变化和水平位移变化分别表示为

(29)

(30)

(31)

(32)

值得注意的是式(29)—(31)中k_n0³³－k_n0²²、h_n0³³－h_n0²²和l_n0³³－l_n0²²，尽管表示两组不同的引张型震源，但在计算中相互耦合在一起，因此作为一个整体来研究.事实上，地震不一定发生在北极轴上，在实际应用中，我们需要将观测点的球面坐标系相应地转换到震源北极坐标系下，通过球面三角关系式进行换算即可.

2 大地测量观测与同震变形的关系及反演方法 2.1 大地测量观测与同震变形的关系

考虑到地震位错Love数是震源深度和球谐方程阶数的函数，我们需要将实际大地测量观测数据转化到谱域中来.GRACE卫星能够检测到地表或地球内部因质量迁移导致的重力变化(National Research Council et al., 1997; Wahr et al., 1998)，用于研究大型俯冲型地震产生的同震和震后变化.Sun和Okubo(2004a, b)认为大于M_W9.0的地震所产生的同震变化可以被GRACE观测到，如2004年苏门答腊地震(M_W9.3)和2011年日本东北地震(M_W9.0).GRACE通过每月提供的Stokes球谐系数计算同震重力场模型.需要注意的是GRACE观测到的大型俯冲型地震都发生于海底俯冲带，地球洋壳的变形会引起海水质量的重新分布，跟固体地球位错理论计算的重力变化相比，该影响不能被忽略(De Linage et al., 2009; Broerse et al., 2011; Sun and Zhou, 2012).这里我们采用一个近似的方法来计算海水改正，选取半无限空间的布格改正层来描述海水重新分布的影响(Heiskanen and Moritz, 1967):

(33)

式中，δg(a, θ, φ)表示海水效应引起的重力变化，ρ_w为海水密度，u_r(a, θ, φ)为球形分层地球模型位错理论的同震垂直位移.O(θ, φ)为海洋函数：

(34)

将其从GRACE重力观测中扣除，这样，我们便得到了GRACE观测值去掉海水效应的同震重力变化在谱域中的计算值.

尽管全球分布了密集的GPS观测网，可以提供丰富的地壳同震位移信息，但大部分地震都发生在洋壳俯冲带地区.海底的GPS数据在地震发生后进行测量，难以分辨同震和震后信息，且包含了震后变形的影响，数据误差较大.但是断层附近的海底GPS观测位移仍然可以成为我们理解地球局部构造信息的约束条件.由于需要计算谱域中的同震垂直位移，局部的GPS观测值无法约束低阶位移信息，我们利用理论垂直位移值替代全球没有GPS观测值的方法来解决此问题.为了减少理论值和观测值之间的误差，这里我们选用球形分层地球模型位错理论作为理论垂直位移的基础.

为了与地震位错Love数建立关系，本研究在全球空间范围Ω内，将GPS垂直位移展开到谱域中，其Stokes系数可写成：

(35)

其中，C_lm^GPS和S_lm^GPS是正规化的系数，P_lm是正规化的勒让德(Legendre)函数.

另一方面，GPS观测值在空间域中各阶的垂直位移变化为：

(36)

其中Ω^r为局部研究区域，(θ′, λ′)∈Ω^r.此时，我们得到了GPS垂直位移变化在谱域中的计算值.大地测量观测技术作为地球内部构造信息的监测手段，有效的处理并得到精确的同震变化信息至关重要.

理论上，在没有任何误差的情况下，地震位错理论同震变化ΔΦ^theory(a, θ, φ)应该与大地测量观测同震变化ΔΦ^observed(a, θ, φ)相等，即：

(37)

实际上，上述观测方程对任意球谐阶数都适用.假设震源类型为剪切型位错源，结合同震变形格林函数与平衡方程，我们得到谱域中每一阶的同震变化观测方程：

(38)

其中lo作为未知的位错Love数(k, h和l)，f₁, f₂和f₃是与地震断层参数有关的量.式(38)中位错Love数跟正演中不一样，正演计算时我们根据地球模型得到地震位错Love数及格林函数，而这里把它们当作未知数.以同震垂直位移为例：

(39)

(40)

(41)

(42)

值得注意的是，由于GRACE数据带有南北条带状噪声，需要高斯滤波来提高信噪比(Wahr et al., 1998).为了使方程平衡，将理论同震重力变化中加入同样滤波半径的高斯平滑.

2.2 反演方法

地球物理学中的反演问题一般可以归结为求解如下方程：

(43)

G为f_i(θ, λ)(i=1, 2, 3)组成的球对称地球模型位错理论同震变化的格林函数，m是未知量地震位错Love数，d为大地测量同震变化，e是模型值与观测值之差.

另一方面，由于观测数据分布的不充足和观测数据误差的影响，方程(43)作为一般线性方程组反演时，难以约束整个位错Love数的分布.为了避免不合理解的出现，需要加入先验模型来约束：

(44)

m₀是先验地球模型(如PREM模型)下得到的位错Love数，ε为先验模型与反演模型之间的差异.结合方程(43)与(44)，反演问题变为求目标函数的最小值，即：

(45)

式中，α为平滑因子或正则化阻尼因子，定义为观测值和先验信息的相对权重(Pollitz et al., 1998).

目前有多种方法可以估计正则化因子α，如L曲线法(Hansen and O′leary, 1993)、广义交叉验证法(Generalized Cross Validation, GCV)(Golub et al., 1979)和赤池贝叶斯信息准则法(Akaike′s Bayesian Information Criterion, ABIC)(Yabuki and Matsu′ura, 1992)等.前两种方法使用较多，L曲线法是寻找正则化解和残差项的平衡点，最终使两部分都比较小，从而确定正则化阻尼因子α.广义交叉验证法与L曲线法相比，容易受到相关性误差的影响.ABIC采用统计学方法，以Bayes算法为基础，寻找失配-粗糙平衡曲线的拐角点，对于简单的问题，ABIC准则与L曲线具有一致性.本研究中选用最常见的L曲线法作为估计α的方法.

通过一维搜索式(45)，得到令该式最小的正则化因子，从而得到最佳位错Love数反演结果：

(46)

为了更有效地描述模型的不确定性，我们通过95%的置信区间来检验反演结果的可靠性(Aster et al., 2018)：

(47)

3 反演实例

结合以上反演理论和方法，以2011年Tohoku-Oki地震(M_W9.0)为例，反演日本岛及其邻近区域的地震位错Love数和格林函数.日本东北大地震发生于太平洋板块和北美板块的俯冲带上，一直俯冲到欧亚板块下，最深可达700 km，属于典型的逆冲型地震.此次地震在日本海沟处引起了明显的同震和震后变形，并被大地测量技术所观测到.本文选用一个典型的地震断层滑动分布模型(Wei et al., 2011)，它是由GSN宽频带地震波形和GPS观测数据联合反演得到，断层分为25×14个子块，尺寸25 km×20 km，走向201°，倾角9°.首先采用了德克萨斯空间研究中心(Center for Space Research at the University of Texas, CSR)每月发布的GRACE L2 RL04数据进行反演，该数据由60阶的Stokes系数组成，时间序列包括了震前96个月和震后36个月的数据.由于地球扁率C₂₀的观测精度较差，我们采用了SLR (Satellite Laser Ranging)观测数据来代替该项(Chen et al., 2007).

为了计算GRACE观测数据的误差分布，将拟合公式中同震阶跃的不确定度(1σ)作为GRACE重力变化的不确定度.图 2所示，误差范围在0.3~0.6 μGal之内，震中附近的不确定度在0.35 μGal左右.该数据结果表明重力卫星GRACE观测得到的同震重力变化可以用于位错Love数的反演.位错理论给出PREM模型下的位错Love数，将其作为先验约束条件，通过0-Tikhonov反演方法得到体现局部构造特征的位错Love数.在选取正则化阻尼因子时，采用L曲线的方法(图 3).

图 4是利用GRACE观测数据反演得到的位错Love数和PREM模型下的位错Love数的对比.每个模型都会得到有限断层层数×阶数n个未知量，这里我们只选取其中一层用来对比分析.从图 4中可以看出，反演值围绕着PREM模型值上下波动，显示研究区域与平均地球模型的区别.图 4显示位错Love数k_n2¹²的低阶反演值大于理论值；当n>20阶时，反演值和理论值在同一水平上；k_n1³²和k_n0³³－k_n0²²的反演结果和理论值都差距较大.结果还表明从浅层到深层两者表现出不同的变化，显示地球结构在深度上存在较大差异.对比3组位错Love数，当j=32，也就是垂直倾滑破裂源时，波动最大.上述反演结果跟我们研究区域的大小和重力卫星GRACE的空间分辨率有关.本研究选取的区域为日本岛及其邻近地区，相对于全球范围来说相对较小，低阶重力场是全球平均信息，不依赖重力卫星观测而改变.从另一个方面说，GRACE卫星的空间分辨率约为300~400 km，并且经过350 km的高斯平滑处理.同样，我们在计算理论同震重力变化时也采用了350 km高斯平滑，所以反演结果不受平滑处理的影响.通过3组反演结果对比分析，对于10~60阶的位错Love数，反演值与PREM模型理论值不同，体现了局部构造信息，达到了本文的研究目的.

图 5给出了同震重力变化的对比结果.图 5d和e分别为GRACE观测值与球形地球模型位错理论同震重力变化和反演结果计算的同震重力变化之间的残差，后者比前者小~80%，说明反演得到的同震重力变化跟实际GRACE卫星观测值符合得更好.本文提出的利用GRACE重力卫星观测数据反演地震位错Love数的方法切实可行，并且得到更精确的格林函数，反映出发震区域局部构造信息，更好地解释大地测量观测现象.

同样，我们选用日本JPL和Caltech的ARIA团队处理的GEONET数据(0.3版本)给出的1232个GPS观测站数据(标准偏差约为1~5 cm)和7个海底GPS观测位移(标准偏差约为0.25 m)，反演2011年日本东北大地震(M_W9.0)垂直位移位错Love数，也得到了类似的反演结果.由于日本地震是一个低角度逆冲型地震，垂直倾滑所占比重较大，所以垂直倾滑破裂源h_n1³²的反演值相对于PREM模型波动最大，h_n2¹²和h_n0³³－h_n0²²变化相对较小.考虑到GPS数据不能约束低阶位移场，当n>10时，反演得到的位错Love数才能体现局部构造模型跟全球平均模型的差别.图 6d显示了理论值与GPS观测值的差距，图 6e是反演值与GPS观测值之间的差别，通过反演得到的位错Love数计算的同震垂直位移比PREM模型下的理论值更加接近于大地测量观测值，并且提高了~60%，可以看出通过本文的研究方法GPS观测值很好的约束了局部地球模型.

4 结论与讨论

现代大地测量技术能够检测到地震引起的地壳变形和地球内部质量迁移，提供精确的同震变化，用于研究地震断层的滑动分布，改进地球物理模型，提高人类对地震孕育发生过程的认知.上述研究是建立在地震位错理论之上，而位错理论中定义的位错Love数的精确度依赖于所采用地球模型的精度，通常情况下地球模型是全球的平均模型(即一维地球模型)，不能体现局部或区域构造特征.利用这样的地球模型计算的理论变形与实际大地测量观测值具有一定的差异.通过实际大地测量观测数据反演位错Love数和4个独立点源的格林函数，其结果包含了地球内部更真实的构造信息.

所以，本文提出了一整套利用GRACE和GPS数据反演局部位错Love数(k, h和l)及格林函数的方法.位错理论只对弹性固体地球有效，而实际大型俯冲型地震发生在海沟俯冲带处，文中使用一个0阶近似布格改正，从GRACE同震重力变化中扣除海水的效应.GRACE数据有南北条带状噪声，在使用高斯平滑时，相应的理论重力变化也要进行平滑处理.局部的GPS值无法约束全球低阶位移场信号，为了得到谱域的同震位移变化，我们利用球形分层地球模型位错理论的同震位移值替代全球没有GPS观测站的数据.以2011年日本东北大地震(M_W9.0)为例，分别利用GRACE重力变化和GPS垂直位移反演局部构造位错Love数和格林函数.结果显示，lo_n1³²相对于PREM模型变化最大，其他两组lo_n2¹²和lo_n0³³－lo_n0²²则变化较小，这与日本地震低角度逆俯冲的滑动方式相吻合.两组实验都证明，通过实际大地测量观测数据反演的位错Love数和格林函数计算的同震变形跟观测值符合得更好，说明反演得到的局部位错Love数和格林函数能体现地球真实的局部构造信息，具备为该地区地震研究提供更加精确的地球模型的潜力.

付广裕(Fu, 2007)研究的三维不均匀地球模型产生的重力变化是基于36阶地球模型计算的，只适用于较小区域三维构造的情况.Aagaard等(2013)、Bagge和Hampel (2017)基于有限元算法建立黏弹性变形的数值计算方法，并且考虑了横向不均匀性.Pan等(2015)通过矢量球谐函数给出了球对称非旋转、弹性和横向各向异性的负荷Love数.基于此，Zhou等(2019)推导出了任意阶数并且高精度的地震位错Love数，从而使高分辨率地球模型广泛应用.以上所有方法都是基于已知的模型以正演的方式建立了位错理论.本研究通过大地测量观测数据反演得到的地震位错Love数和格林函数比正演地球模型更加符合该地区的构造特征，通过此模型反演的断层滑动分布，避免了将地球模型的误差带入到反演结果中.

根据上述方法，我们可以研究地震多发地区的地球构造模型，如日本岛及其附近、川滇地区、智利等区域，给出体现局部三维构造特征的位错Love数和格林函数，建立区域构造模型的数据库.通过反演得到的位错Love数和格林函数调整地球浅层局部构造模型参数，更好地解释大地震再次发生而产生的同震变化.随着GRACE-follow-on的发射和越来越多GPS观测站的建立，空间分辨率和噪声都会得到很大程度的改善，能够确定的重力场和位移场阶数将增加，预期能得到更好的位错Love数反演结果.