SNREI地球形变理论模拟研究进展

孙和平，周江存，徐建桥

1. 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室，湖北 武汉 430077; 2. 中国科学院大学，北京 100049

我们生活的地球是太阳系中唯一生机盎然的行星，地球内部发生许多动力学过程，如构造活动、地震等，表面受到大气、海洋等的负荷作用，同时也会受到来自太阳系其他天体的引潮力作用，以及太阳的热辐射等。伴随这些过程与力的作用，地球会产生形变。有些形变是能够直观地感受到的，比如地震产生的破坏性错动；有些形变只有通过高精度仪器设备才能探测到，如潮汐形变，尽管它在地面的垂向位移幅度有数十厘米之大[1]，但因为没有相对运动，所以人们站在地表是无法感知到这种形变的。实际上，人们对于地球形变的研究，不仅仅局限于狭义的几何量(或运动学变量)，如地面的位移、倾斜、应变等，而且拓展到了许多物理量，如外部重力场的变化、内部的应力场变化等动力学变量。随着针对不同观测量的仪器观测精度的提高，对于地球形变的研究不仅可以检验形变理论的正确性，还可以进一步完善地球模型与分析和认知不同的动力学过程。实际上，这些不同的动力学过程都可以用等效的力源来表示，如地震的等效体力源、火山的Mogi源、流体运动的膨胀源等。因此，研究地球的形变就与不同力源的作用紧密联系在了一起。

研究地球形变通常基于两种坐标系统：柱坐标系和球坐标系。前者通常适用于研究小范围形变，如中小地震引起的近场形变；后者适用于研究全球范围内的大尺度形变，如大地震引起的中远场形变、海潮负荷引起的全球形变等。关于地球受力形变的研究理论最早可能要追溯到对平面半无限空间负荷问题的系统工作，该理论主要涉及平衡方程与本构关系方程[2]。随着观测精度的提高，对理论的发展也提出了新的要求。如GNSS、InSAR、DAS(分布式声速传感器)可以探测毫米级乃至更微小的形变信息；重力测量的精度也达到了μgal量级，如FG-5重力仪，特别是超导重力仪的探测精度甚至要再高2个数量级。此外，倾斜仪、应变仪等观测精度也越来越高。因此，通常在研究全球性的大尺度形变时，除了需要采用球坐标系外[3-4]，还需要考虑地球的分层结构、自引力及由变形导致的质量调整产生的附加引力位等。因此，在平衡方程中需要考虑重力及引入泊松方程[5-6]。地震反演研究表明，地球更加接近于一个弹性的分层球体[7]。因此，现今地球形变的研究大多以此为基本模型。

为了方便研究，通常假设地球是球对称、不自转、完全弹性和各向同性的，也就是所谓的SNREI(spherically symmetric,non-rotating,elastic and isotropic)地球[8]，如目前常用的PREM地球模型，除了考虑地幔的各向异性外，也给出了一个等效的各向同性模型[7]。此外，采用该类模型研究地球形变常常隐含地包括地球的可压缩性与自引力。采用SNREI模型的一个方便之处在于球型形变与环型形变的解耦，从而简化数值计算。然而对于研究动态的形变问题，特别是地球自由振荡问题，常常需要考虑地球的椭率与自转[8-10]，这两个因素不仅会导致球型形变与环型形变的耦合，还会引起振型的谱峰发生分裂，即所谓的正常谱峰分裂，相关内容可以参考文献[11—12]。

地球表面的动力过程主要体现为全球水循环，它导致水在大气、海洋与地表和地下不停运移与交换。在地球的表面，约70%是海洋，在引潮力作用下，海洋产生潮起潮落现象，即海洋潮汐，简称海潮。而海潮产生的海水质量的时空变化会改变固体地球表面的荷载，由此导致的影响称为海潮负荷效应。此外，大气中的水汽含量、温度、气压等导致其密度产生变化，也就是质量的时空变化，它也会使地球发生形变。同时大气在地表的温度、气压的变化也会使地球发生形变。与海潮负荷类似，地表的河流、地下水等也会产生负荷效应。

不管是地球外部、表面或内部的力源，地球变形过程都可以用基本的物理方程与地球介质的弹性参数来描述，因此，可以建立相应的初边值问题，从而模拟出变形量的大小。地球的形变通常用无量纲的勒夫数来描述[3,13]，它是一个比值，反映了地球在力的作用下的反应强度。勒夫数大，说明单位力作用下的变形量也大。因此，对于地球形变问题的研究，基础的工作就是基于一个给定的地球模型，计算不同的勒夫数与相应的格林函数，并结合实测数据进行模拟计算与分析。鉴于地球固体潮的研究已经非常成熟，相关改正方法和参数可以参考IERS规范。本文着重介绍地表质量和热负荷与地震引起的地球形变方面的理论模拟研究进展情况。

1 地球形变理论初边值问题求解

初边值问题涉及泛定方程与边界条件和初始条件[14]，泛定方程通常是用数学方法对基本的物理定理或定律进行描述所获得的一组微分方程。对于地球而言，主要包括由牛顿第二定律描述的运动方程，对于不考虑惯性力时的静态形变问题也称为平衡方程；联系弹性体(即地球介质)应力和应变的本构关系，即虎克定律；描述扰动位与扰动密度的泊松方程。如果要考虑与温度场有关的热效应，需要考虑热传导方程。对于SNREI地球来说，可以将所有变量展开为球函数级数。在理论研究中，采用点力源是非常方便的，并将这个力源置于球坐标的极轴上，如地表负荷源置于北极点，内部位错源的震中置于北极点。同时，对力源同样做球谐展开，这样就可以得到关于球谐展开系数的微分方程，这些方程最后都可以写成如下统一的形式

(1)

式中，r是球坐标系下的矢径；Y是由不同变量构成的矢量；A是系数矩阵；F是(等效的)单位体积力的球谐展开系数。对于力源在地表的负荷问题，可直接用F建立地表边界条件[6,15-16]；对于力源在内部的地震位错问题，可将F转换为内部震源处的不连续边界条件[5]。因此，最终只需要求解齐次的微分方程组，即式(1)右边不再出现F。在地心处，满足正则条件，即位移场、附加引力位和温度梯度为零。如果地表没有力源，则需要满足自由边界条件，即正应力与剪应力为零。此外，如果方程涉及时间，还需要初始条件。这样通过求解这个初边值问题就可以得到任意位置的Y，从而得到勒夫数

(2)

式中，y1、y3、y5分别表示Y中的垂直位移、水平位移和附加引力位分量，对于有直接引力项的力源δ=1(如地表质量负荷)，其他力源δ=0，n和m分别是球谐展开的阶和次。对于地表负荷问题，m=0；对于位错问题，m=0,1,2。这样，就可得到格林函数

(3)

式中，u、v、w、φ分别为垂直位移、切向位移、横向位移、附加引力位；N为构建边界条件时的规格化系数；Pnm为Legendre函数；ψ是计算点到源点的角距离；φ是方位角；最右边的括号适用于位错问题，kc和ks是某个常数，至于是取余弦还是正弦取决于不同位移场分量与位错模式[17-18]。

对于SNREI模型，传统的研究地球纯弹性形变的方法是采用四阶Runge-Kutta方法[19-20]。需要注意的是，变量之间的相关性会随着球谐阶数的增大而增强，因此，当球谐阶数非常大时，会导致最后确定待定系数的矩阵奇异。文献[20]注意到变量包含因子rn，因此，先从变量中扣除rn再进行计算，再结合积分过程中的归一化方法，可有效改善相关性问题。

对上述的初边值问题，可采用有限元方法的思路直接进行数值求解，或者从最基本的方程出发构造弱解形式。文献[21]在有限元方法的基础上发展了谱有限元方法，在球面上采用球谐函数展开，在矢径方向采用有限元，但该方法没有考虑地球的可压缩性及液态的外核。文献[22—23]将该方法扩展到可压缩地球模型，但同样没有考虑液态的外核。文献[24]借鉴文献[19]的方法，引入势函数并假设液核初始模型满足Adams-Williamson条件，推导出液核中平衡方程及其对应的弱解形式，使谱有限元方法能够应用于更加真实的地球模型。

文献[25]提出了一种近似解析解方法，通过合理假设地球介质分层参数随r的变化关系，得到了常系数矩阵A，从而得到了每一层的解析解，然后通过传播矩阵方法与边界条件得到最终的解。文献[26]曾经也给出一个解析解，但它的假设与真实地球不太符合。不同于文献[26]的假设，文献[25]的假设更加合理，只要层厚足够小，通过这种假设所得模型与最初模型是完全一致的。数值试验表明，实际上层厚可以取得很大，就能得到满意的结果，因此这种方法非常高效，但对球谐展开的高阶项存在数值溢出问题。文献[27]采用传播矩阵中的刚度矩阵方法解决了溢出问题。文献[28]对比了多种高效的传播矩阵方法。关于不同求解方法的讨论可参考文献[29]。

2 研究进展

2.1 地表质量负荷

在年际以下的时间尺度上，地球表面的质量调整主要是全球的水循环，体现在海洋、大气与陆地水的变化中，其中海洋的影响又分为潮汐与非潮汐效应。海洋与陆地水的影响可用式(4)表示(以地表的垂直位移为例)

(4)

式中，a是地球半径；θ和λ分别是计算点的余纬与经度；θ′和λ′分别是源点(即负荷作用点)的余纬与经度；ρ是海水或陆地水密度；H是海潮潮高或海面高变化量或等效水高；u称为垂直位移的格林函数；ψ与源点和计算点的坐标满足式(5)

cosψ=cosθcosθ′+sinθsinθ′cos(λ-λ′)

(5)

除了采用式(4)的数值积分外，还可以采用球谐展开方法。球谐展开方法就是将H展开为球谐级数，u采用式(3)中的形式，根据球谐函数的正交关系，得到式(6)

(6)

式中，ρe是地球平均密度(来自于N)，带上标C和S的Hnm表示H的球谐展开系数，注意这里为了方便计算，采用了大地测量学中使用的完全规格化的球函数(以上横线加以区别)，并且N的值为单位引力位。文献[30]讨论了重力负荷计算中关于直接引力项系数的选取问题。文献[31]提出了一种加速式(6)收敛的方法，并指出由于负荷质量分布不连续引起的吉布斯效应在沿海地区对结果会有很大的影响。文献[32]采用该方法研究了海潮负荷对垂直位移的影响。数值积分方法的计算量大、精度高，而球谐展开方法由于比较简单，计算量小但精度相对低。鉴于此，文献[33]提出了近区采用数值积分、远区采用球谐展开的方法。文献[34]对这两种方法进行了比较，计算了海潮负荷导致的垂直位移，指出在沿海区域由两种方法得到的结果平均差异为3 mm，最大达到10 mm。文献[35]给出了比较完整的计算重力、位移、倾斜和应变的球谐展开方法的计算公式，通过对两种方法计算的结果的比较(表1)，也指出该方法只适合于内陆地区。

表1 不同海潮负荷计算方法在上海站的比较[35]

随着计算机技术的发展，数值计算效率得到了很大的提高，因此应选择用式(4)进行数值计算。目前，有许多开放软件可以计算式(4)，如文献[36]编写的SPOTL包中的nloadf函数，该方法采用了文献[37]提出的积分格林函数方法，有效提高了计算效率，同时考虑了不同海洋区域的海水密度变化影响。文献[38]在GOTIC基础上发布的改进版本GOTIC2；Bos和Scherneck提供的计算海潮负荷的网页版(http:∥holt.oso.chalmers.se/loading/)。文献[39]发布的SGOTL，随着Python的流行，也出现了Python版本的软件[40]。

由式(4)可知，两个关键的量是H和u，前者对应着负荷源模型，如海潮模型，后者对应着地球模型。实际上，H中隐含一个关于海陆边界的信息，由于海潮或海面变化只存在于海洋，而陆地水和大气负荷(考虑反变气压计效应)只存在于陆地。因此，在研究陆地水的负荷效应时，只需要改变海洋和陆地的判断就可以直接使用计算海潮负荷的程序进行计算。STRM地形数据为确定高空间分辨率的海陆边界提供了直接的资料，利用分辨率为90 m的地形数据优化计算点周边的海陆边界。文献[41]计算了沿海台站的海潮负荷效应对重力的影响，对比了由不同海陆边界获得的结果，并指出了高精度海陆边界对于准确计算海潮负荷的重要性(图1)。后来在计算中对中国沿海全域的海陆边界进行了优化，构建了我国沿海和岛屿的海潮负荷模型[42]。

图1 利用不同分辨率的海陆边界计算的不同沿海台站重力海潮负荷[41]

得益于卫星测高技术的发展，流体动力学模拟与验潮站资料的增多，海潮模型的精度得到了很大的提高，特别是近海潮汐模型，已广泛应用于海潮负荷效应的模拟。全球有多家机构从事全球大气、海洋和陆地水模型的构建工作，如FES系列的全球海潮模型[43]，ECCO模型(www.ecco-group.org)、CPC模型(www.cpc.ncep.noaa.gov/products/GODAS)和GLDAS模型(earth.gsfc.nasa.gov)等。通过这些模型，开展了许多负荷效应的模拟应用工作[44-53]。文献[54—55]研究了海潮模型在重力海潮负荷改正中的适定性问题。文献[56]研究了利用不同海潮模型用于重力改正对解算的自由核章动参数的影响。这些研究为海潮模型的选择提供了重要的参考。近年来，InSAR技术广泛应用于地球形变的研究，文献[57]研究了海潮负荷对InSAR视线方向形变的影响。文献[58—59]比较了香港地区部分GPS站得到的海潮负荷和模型的比较，类似的工作也用于新西兰GPS站观测数据的开展。

由不同地球模型计算获得的格林函数是不同的。如果考虑不同区域地壳的参数，如CRUST1.0，则格林函数在近区会有很大的不同[41,60-61]。图2显示了两种地球模型计算的重力负荷格林函数在近场随角距的变化，从中可以看出很明显的差异。

图2 不同地球模型对海潮负荷格林函数计算的影响[41]

在模拟大气对重力观测影响时，不同于研究海洋的影响，由于大气的垂向分布范围较广，不能将其看成单层，而需要考虑大气的垂向分布来计算直接引力。文献[62]采用一个服从指数分布的大气模型计算了直接引力项，并与计算点的局部大气压强进行比较得到了一个大气重力导纳值(直接引力/局部压强)，发现结果比用平面大气模型小7%～8%，说明了考虑地球曲率的必要性。文献[63]指出，大气负荷产生的重力变化达到30 μgal，其中90%由计算点周边50 km的大气引起，同时地表的大气温度发挥很大的作用。文献[64]通过标准大气模型计算了大气重力格林函数，为进行全球大气重力负荷效应的积分计算奠定了基础。文献[16]也计算了大气重力格林函数，但采用不同的方法处理直接引力项，即除了对地表的大气压强产生的正应力负荷作用外，将大气的引力看成是固体地球之外的作用(实际上该力可看成引潮力)，然后根据大气的空间分布进行积分得到直接引力项的格林函数。实际上，由于采用的是基于小形变的线性理论，因此地球对地表质量负荷的响应等价于潮汐响应与地表正压力响应的叠加。因此，文献[16]的结果与以前的结果是一致的。然而实际研究中，由于三维大气观测资料的缺乏，常采用大气重力导纳值进行大气负荷效应的扣除。大气负荷除了引起重力变化外，同样会引起位移、倾斜和应变等变化，而模拟它们的方法与海潮负荷中的方法是一致。

2.2 地表热负荷

由于四季变化与日夜扭转，地表的温度是时刻变化的。热胀冷缩对于地球也不例外，关于温度与应变的关系在材料力学领域研究得比较多，因此在半无限空间体系下比较成熟。虽然人们对地球上典型岩石的热力学参数都有很好的了解[65]，但是整个地球在地表热负荷作用下的形变问题研究并不多。文献[66]研究了分层球体的热弹性变形问题，并用差分方法进行了数值计算，但跟平面半空间模型下的理论一样，没有考虑任何的体力。文献[67]研究了均匀球的热弹性变形问题，同样也没有考虑体力项。并且为了方便，在理论推导过程中作了一定的假设，因此，在球谐阶数超过一定的值之后会有较大误差。采用球谐展开的方法计算了地表温度的周年变化对位移场的影响，发现与采用平面半无限空间模型得到的结果有较大的差异，进一步说明了采用球坐标系的必要性。对于地球来说，自引力是不容忽视的。文献[68]利用全耦合的方程，考虑预应力和自引力，并采用文献[25]的近似解析解方法及DVP传播矩阵方法，求得了任意球谐展开阶数的热负荷勒夫数。并且在不作任何假设的情况下，给出了不考虑自引力的均匀球的严格解析解，其中热负荷勒夫数是以球贝塞尔函数表示的。该理论为研究SNREI地球的地球热源引起的地球形变，特别是重力场变化，奠定了理论基础。

2.3 同震形变

最早关于地球地震位错形变的研究采用的是平面柱坐标，如经典的Okada解析公式[69]。采用柱坐标系的主要原因可能源于材料或工程研究领域的理论传承。文献[70—71]将相关理论拓展到球形地球模型并给出了解析解，但采用的模型没有考虑重力场耦合，即自引力效应。文献[72]采用了1个考虑预应力、自引力、液态外核，径向非均匀等更接近真实地球的模型开展了相关工作。文献[17,73—74]在此基础上作了深入的系统性研究，采用文献[5]的方法模拟了点源位错引起的位移场和重力场等变化。通过将有限断层离散化，将每个子断层认为是点源，得到了有限断层产生的同震形变[75]。相关的理论体系可参考文献[76]。

文献[4]比较了由平面模型和球模型所得结果，说明在近场二者是一致的。但是后续的研究表明，在中远场二者的差异较明显，说明了使用球形地球模型的必要性，这一点也被后来的研究进一步强调[77-78]。

文献[79]提出了计算同震形变的互换定理，即内部的位错源产生的地表形变可以用3种地表力源作用下在位错源所在位置产生的形变的组合来表示，这使得同震形变的计算更加简单，因为这使得位错形变问题只涉及3组在地心正则的解，或者说只需要解齐次的微分方程组而不需要考虑力源项(从力源转化的内部不连续性或特解项)。文献[80]给出了上述3种地表力源引起的地球变形的渐近解。因此利用互换定理，文献[81—83]研究了同震形变的渐近解，该方法也被用来计算断层在地表时的位错勒夫数渐近值[84]。文献[85]也用该方法加速内部同震体应变格林函数计算的收敛，说明渐近值在格林函数计算中的作用。文献[3]给出了均匀球的3组在地心正则的解。因此，可以直接运用互换定理计算均匀球内的位错产生的地表形变，并且采用文献[80]的渐近解可确保格林函数的收敛[86]。但是，对于一个更加真实的分层地球，均匀球的渐近值无法被应用。

文献[18]利用文献[25]的近似解析解方法和DVP传播矩阵法研究了位错勒夫数的计算问题，他们的方法解决了以前的方法在超高阶的计算中的变量相关导致矩阵奇异的问题，因此，可以计算到很高的阶数。由于用不同于文献[74]的位错勒夫数定义，能够通过观察高阶位错勒夫数值和球谐阶数的关系直观地发现二者之间的联系，并确定了拟合公式与位错勒夫数渐近值，从而克服了特殊情形下格林函数不收敛的计算难题[87]。渐近值加速格林函数收敛的作用可从式(7)看出，同样以垂直位移为例

(7)

式中，第1个求和的差值随着阶数n的增大衰减更快，而第2个求和有严格的解析式[87-88]，因此，能够很快得到收敛的格林函数。文献[89]也对地震引起内部变形进行了研究，同样获得了收敛的格林函数。当球谐阶数很高时，重力的影响越来越不明显，因此，通过忽略自引力并运用DVP传播矩阵的特殊性质。文献[90]在理论上获得了位错勒夫数的渐近值，从而彻底解决了位错格林函数的计算问题。

最近十几年，发生了多次震级较大的地震，这些大地震引起了全球尺度的显著变形[91]，同时也对地球自转[92-93]、地球参考框架[94-95]、地球的体积变化[96](图3)，以及地球内部能量的转换产生了较大影响[97-99]。这些研究是对地震形变模拟理论的重要应用。

图3 地震对地球体积大小的累积影响[96]

3 总结与展望

本文回顾了地球弹性变形的理论、地表负荷和地震位错引起形变的理论模拟方面的主要研究进展，地球弹性变形理论发展至今，可以说已经非常完善，相关的研究成果也非常多。然而随着技术的发展，观测精度越来越高，地球变形的理论也需要得到更进一步的发展。

在短时间尺度下，地球表现出显著的弹性特征，因此，地球的弹性形变模拟在目前大地测量和地球物理学的研究中发挥着非常重要的作用。因为地球的形变模拟提供了联系力源与形变的纽带——格林函数，因此一方面可以确定不同力源导致的地球形变，另一方面也可以在反演地球内部结构、物性参数，以及力源参数、甚至动力过程中发挥重要作用。地球是一个复杂的系统，有太多的因素会导致地球的形变。随着高精度观测资料的获取，会发现更多观测与理论不相符合的情形，或用现有理论无法解释的现象。比如越来越多的证据表明，在地热开采区或注水区的地震具有明显的非双力偶特征，并认为这可能与流体有关[100]，然而具体的机制仍不明确，对震后孔隙弹性回弹的认识仍有待加强，因此，需要发展孔隙弹性位错理论；地温观测发现了明显的温度场同震变化[101]和温度的潮汐变化信号[102]，前者与构造活动的关系及后者与引潮力的联系都是值得深入研究的课题，因此，需要发展地球的热弹性变形理论，包括引潮力作用与地震位错作用的相应理论。热弹性理论不仅可用于计算由于温度场变化导致的地球形变，还可以研究变形引起的温度场变化，尽管人们一直认为后者的效应是很小的，通常采用所谓单向耦合的理论，但是地球内部温度的测量目前已经能够达到非常高的精度(0.1 mK)，因此，需要发展双向耦合的热弹性理论及开展相关的应用研究。

在长时间尺度上，地球表现出滞弹性的特征，相关的研究涉及冰后均衡调整与震后黏弹性松弛，根据对应性原理，其在拉普拉斯域的边值问题与弹性边值问题是一致的，因此，相关的弹性理论可以直接应用，只需要最后对结果实施反拉普拉斯变化即可。相关的黏弹性形变的研究进展可参考文献[103]。

我国正大力发展非常规能源的开采(如页岩气和地热等)及水电开发，要评估这类活动对环境的影响就涉及流体与温度场，即孔隙弹性与热弹性的问题，或将二者统一的孔隙热弹性的问题。现在也有一些基于柱坐标系统理论的研究，但是现有理论可能已经不能满足当今观测精度越来越高的需求。因为这些理论没有考虑预应力与地球的自引力，然而这些会影响最终的计算结果[104]。考虑自引力后，由于平面半无限空间模型天然的缺陷，需要作某种假设[105]。因此，作为多学科交叉的地球形变问题，有必要发展球坐标系统下的相应理论对观测结果或现象作客观的描述与解释，同时也有必要开展观测网络的建设获取更加丰富的高精度观测数据，以增强我们对地球动力过程的认知。