1 引言

固体潮源于月亮和太阳以及其他天体的引力作用^[1-3].固体潮的影响可以通过安置在地表的高精度仪器观测出来，如用于地震研究的应变仪^[4-5].应变反映了质体内任一点因外力作用引起的形状和尺寸的相对改变，它与地球内部的应力状态密切相关，因此可有效反映地震前后的应力分布和传输.应变固体潮是应变测量时间序列中的最主要的部分，因此从观测资料中有效扣除应变固体潮的影响，并应用于地球内部介质特征和地震的研究具有非常重要的意义.

应变或应力固体潮的理论值通常根据地表的Love数计算获得^[1-3]，而实际上应变仪是安装在地球内部，因此采用地球内部的Love数来计算理论的应变或应力固体潮则更加符合实际.除此之外，地球内部应力固体潮可用于地震的潮汐触发的研究^[6-8]，潮汐应力对地震临界状态的影响有可能是导致地震最终发生的力源.尽管地球内部的应力潮汐已经应用于潮汐触发地震的研究^[6-8]，但具体的关于地球内部应力的计算问题则语焉不详.吴庆鹏曾从点源的角度，根据地球在引潮力作用下的变形讨论了固体潮应力，给出了球谐展开的二阶和三阶的地球内部函数结果；但是并没有给出最终的结果，而是给出了以天体天顶距为参数的公式^[9].因此如需获得地球内部任一点的应力在地球地固系下的应力和应变仍需进行进一步的推导和计算；另外目前的观测精度越来越高，球谐展开的四阶项尽管小，但也是需要考虑的.孙文科也曾对地球内部的重力、倾斜固体潮潮汐因子进行了研究，并且给出了北京地下300 km应力潮汐的计算实例，表明应力潮汐的振幅达到10³ Pa^[10].

为了让读者更加深入详细地了解地球内部应变和应力潮汐的理论和计算问题，以及应力和应变潮汐在地球内部不同区域的变化特征，本文对此进行了详细的介绍与讨论，分析了应变和应力固体潮在地球内部不同深度和不同纬度处的时间变化特征.作为应力固体潮在潮汐触发地震研究中的应用实例，本文也探讨了断层面潮汐应力相位与地震发震时刻的相关性.

2 计算理论 2.1 引潮位及其导数的计算

根据引潮位展开表，引潮位可以表示为^[2]

(1)

其中，G_nm是大地系数，n是球谐展开的阶数，m是展开的次数，其二阶表达式为^[2]

(2)

三阶表达式为

(3)

四阶表达式为

(4)

其中D₂称为Doodson常数，r是地球内部一点到地心的距离，R是地球半径，Φ为纬度，A_ω是频率为ω的潮波的相对幅度，Δ_φω是相应的相位改正数（由于有些正弦波要统一表示为余弦函数的形式），m取0、1、2、3和4对应于长周期、周日、半日、1/3日和1/4日，τ是平地方太阴时，它与世界时t和计算点经度的关系如下：

(5)

（1）中的省略号部分是包含s和h等5个天文参数的线性组合，与τ一起，6个天文参数的系数构成Doodson编码^[1-2].

需要注意的是，传统的Doodson展开方法和现代引潮位展开方法略有不同，如使大地系数保持不变，则现代方法展开的A_ω需要除以一个转换系数（表 1），以便得到正确的结果^[11].

2.2 应力与应变固体潮

通常所使用的Love数是研究地表固体潮现象的重要参数，固体潮产生的地表位移的三个分量可以用Love数h_n和l_n来表示，对于地球内部问题，同样可以采用Love数.我们采用如下的Love数的定义^[12]：

(6)

式中，引潮位的展开系数取为1，而没有选择文献[12]中的Rg₀，这种归一化系数的选取并不会改变最终的结果.那么，地球内部或表面任一点的垂直和水平位移为

(7)

其中，U，V_θ和V_λ分别是垂直位移、余纬向和经向的水平位移，T_n是n阶引潮位，θ是计算点的余纬，g₀是地表重力（加速度）值，N是截断的阶数，目前通常为4阶.由（2）-（4）可知，引潮位中包含（r/R）ⁿ因子，而Love数中含有（R/r）ⁿ因子，因此直接用习惯的y₁和y₃来代替内部Love数则使得这两个因子相互抵消，从而使公式更加简洁，其表达式如下：

(8)

值得注意的是，文献[10]定义的Love数为

(9)

但是在此Love数定义下计算地表位移时，采用了地表的引潮位，所以位移的表达式与（8）式一致.

根据应变与位移的关系，可得应变张量中6个独立分量的表达式为

(10)

再根据应力与应变的关系：

(11)

即可获得地球内部的潮汐应力，其中μ和λ′是r处的拉梅常数， 是应变的体膨胀，即应变张量矩阵的迹，δ_ij是Kronecker符号.

从上述的公式可以看出，计算潮汐应力和应变的关键是求得引潮位及其1阶和2阶偏导数，它们可以根据（2）-（5）式推导获得；另外需要计算y₁和y₃对r的导数，它们是y₁、y₂、y₃和y₄以及地球的拉梅常数的函数，其中y₁、y₂、y₃和y₄分别表示在引潮力作用下地球内部和表面产生的径向位移、径向应力、水平位移和水平应力，它们的数值结果可以通过求解基于某一地球模型的微分方程组的边值问题获得^{[2, 10, 13]}.

2.3 计算步骤

综上所述，可以将计算步骤总结如下：

（i）根据固体潮理论求解微分方程组边值问题获得y₁、y₂、y₃和y₄随r变化的值，并导出y₁和y₃对r的导数值；（ii）根据（1）-（5）式求得引潮位（r=R）对Φ和λ的偏导数，其中对λ的偏导数即为对τ的偏导数；（iii）根据（10）和（11）式分别计算获得固体潮应变和应力.

3 计算结果与讨论 3.1 地球内部Love数

计算中我们采用了球对称、非自转、弹性和各项同性的PREM地球模型^[14]，但对此模型进行了修正，对液态的海洋层和地壳最上层进行了平均^[15]，根据微分方程组的边值问题求解获得了对应于2、3、4阶的y₁、y₂、y₃和y₄随r变化的值，以及y₁和y₃对r的导数值，结果如图 1，Love数可由（6）式获得.

我们仅给出了从地表到800 km深度处的值，不管是对于应力应变固体潮的研究还是与地震相关的研究，这个深度已经足够.y₁~y₄都是r的连续函数，而y₁和y₃对r的导数由于拉梅常数的不连续在间断面处是不连续的.y₂~y₄都随深度的增大而增大；而y₁随着阶数的不同其变化特征也不相同.我们将该结果与文献[10]进行了对比：在地表，本文和文献[10]的h₂值分别为0.5988和0.5810，l₂值分别为0.0846和0.0903，h₂的导数值分别为-0.4229×10^-7和-0.3410×10^-7（转换到国际单位，下同），l₂的导数值分别为-0.8071×10^-7和-0.4910×10^-7.由以上数值可知，本文和文献[10] h₂的结果比较一致，而l₂的差异稍大，但本文l₂结果与文献[13]的约0.0840比较接近.我们知道，改变最表层的地球模型参数，对地表低阶的Love数影响甚微，但是会明显改变h₂和l₂的导数值，因为它们与拉梅常数密切相关.对于h₃和l₃及其导数值以及在地球内部也有类似的结论.

3.2 应力与应变潮汐

图 2给出了经度和纬度都为0的地方应力的6个独立分量对应于不同深度、40h内（中间时刻为2013年1月1日0时）的时间变化，其中引潮位展开表采用Xi的结果^[16].由图可以看出：深度变化对应力潮汐的影响非常显著，随着深度的增加，应力潮汐的振幅也随之增大.rr、rθ和rλ应力分量在浅层的振幅非常小，只有10 Pa的量级，在地球深部其振幅达到10³ Pa；其他分量的振幅在地球浅层达到10² Pa，到地球深部达到10³ Pa；除了rθ应力分量外，其他应力分量的潮汐都表现为明显的半日波特征；三个正应力分量的潮汐时间变化在相位上具有较好的一致性.

图 2同时也给出了相对应的应变潮汐的变化，由于剪切应变与剪切应力只相差一个比例系数，因此它们的变化特征一致，所以图中我们只给出了正应变的结果.从图 2可以看出：正应变的振幅都在10^-8量级，随着深度的增大，应变潮汐的振幅都有减小的趋势，其中rr分量的正应变潮汐振幅对深度的依赖关系与图 1中的dy₁/dr一致，θθ和λλ分量的正应变潮汐振幅对深度的依赖关系也可由图 1中的y₁和y₃获得.

为了更加清晰地了解深度对应力和应变潮汐的影响，我们对深度从10 km到400 km、间隔10 km处的潮汐曲线在与图 2相同时间段内取其最大值（图中实线）和最小值（图中虚线），最大和最小值的差值反映了潮汐的振幅及变化幅度，结果如图 3所示.从图 3中可以看出，应力的6个分量的潮汐变化幅度都随深度增大而增大，其中θθ、λλ和θλ分量由于地球模型参数在内部边界的不连续而产生不连续的变化，其他3个应力分量潮汐的变化幅度与深度近似呈现线性增大的关系；rr分量的正应变随深度减小，在大约250 km深度处，应变增大，然后又随深度减小；θθ分量正应变分量潮汐的变化幅度随深度变化很小，λλ分量正应变分量潮汐的变化幅度随深度增大逐渐减小.

图 4给出了应力潮汐6个独立分量在不同纬度处的时间变化，对于同一纬度，经度的不同只会导致相位的变化，并不会导致振幅的变化，因此我们选择零子午线上纬度从北纬90°到南纬90°，间隔10°，深度为20 km的地点.从图中可以看出：6个应力分量除了在两极外都表现出明显的半日波变化特征，而rθ和rλ应力分量的潮汐在两极表现为周日波特征，θλ分量的潮汐在赤道表现为周日波特征；rθ分量在南北两极的值是有一些差异的，这是因为球谐展开的三阶潮汐项中的周日项是大小相等而符号相反的；rλ分量的潮汐在南北两极是对称的；需要指出的是图中θθ分量在两极处也是半日波特征，只是由于振幅相对其他纬度区域较小，并不是图中显示的长周期特征；其实这些变化特征也不难从大地系数及其导数推断出来.同样地，图 4也给出了3个正应变潮汐随纬度的变化，其变化都比较剧烈，rr分量的正应变潮汐在两极表现为长周期特征，其他的潮汐变化主要表现为半日波特征，正应变潮汐也没有明显的纬度上的对称性.

图 5给出了应力和应变潮汐的变化幅度随纬度的变化关系，图中每间隔1°计算一个值.同样地，实线表示最大值，虚线表示最小值.除了应力的rθ分量表现出关于南北纬度的反对称外，其他的应力和应变潮汐的变化幅度都表现为南北纬度的对称特征，但是从图 4可知它们不是严格对称的，只是在变化趋势上具有明显的这一特性，这一结论也可通过仔细观测图 5中的变化曲线得出.通过最大值和最小值的差值我们可以获得应力和应变潮汐在哪一纬度区变化幅度最大，如rr正应力分量潮汐的最大变化幅度出现在两极附近；对于rr正应变，其值在南北纬20°附近.

3.3 应力潮汐应用

作为一个应用实例，我们研究了潮汐触发地震的问题.选取了我国中部地区（纬度范围25°N-40°N，经度范围100°E-106°E）从1976-2010年发生的共101个地震，地震资料从GCMT获得（http://www.globalcmt.org/[2013-01-11]），该资料提供地震发生的时刻、震源三维坐标及断层方位参数，即走向角、倾角和滑动角.由此，我们计算了震源处的潮汐应力，并根据应力的坐标旋转计算了地震断层面上潮汐的正应力、剪应力、应力张量的第一不变量（张量矩阵的迹）和库伦破坏应力的时间序列.

对于距离海洋较近的地区，海潮负荷甚至比固体潮都要大^[17]，而距离海洋较远的区域，海潮负荷的影响比固体潮小大约一个数量级.为了对海潮负荷如何影响潮汐的振幅与相位有一个具体的认识，以远离海洋的中部地区的汶川地震的震源参数为例，我们首先计算了该地震断层面上的潮汐应力，然后计算离海洋较近的上海台的应力潮汐（采用相同的参数），图 6给出了断层面上的正应力和剪应力潮汐变化曲线，由于上海台（图 6右边一列）的海潮负荷影响较大，因此在表示固体潮和海潮负荷共同影响时，在图的右边采用了不同的虚线标尺.由图可知，对于远离海洋的地区，尽管海潮负荷对潮汐振幅的贡献不大，但是对潮汐的相位还是有较大的影响；而对于离海洋较近的地区，海潮负荷不仅对潮汐振幅而且对潮汐相位都有非常大的影响.因此我们在计算时也加入了海潮负荷的影响.关于如何计算内部的海潮负荷应力和应变潮汐可参考文献[17].

由断层面上的应力潮汐变化的时间序列，根据发震时刻和文献[8]的方法获得发震时刻的潮汐相位和Schuster统计结果的p值，其中p=e^-A/N，A=（Σcosψ_i）² +（Σsinψ_i）²，ψ_i是某一地震发震时刻的潮汐相位.p值反映了拒绝地震随机发生的零假设的显著性水平，其值越小说明地震的发生越不是随机的.在此，我们选取p < 5%，也就是95%的置信水平，来判定地震发生与断层面上应力潮汐的相关性.

图 7给出了地震发生时刻应力潮汐相位的分布特征，区间间隔为30°，统计了在每个应力潮汐相位区间的地震数目的百分比，图中虚线是拟合的余弦函数，由此函数可以看出发震时刻潮汐相位的分布特征，同时代表相关性大小的p值也在图中给出.由图可以看出，正应力潮汐和应力第一不变量潮汐与地震都具有明显的相关性，其p值分别为0.45%和0.55%；而剪应力潮汐、以及库伦破坏应力潮汐都与地震没有显著的相关性，因为其p值分别为58.29%、87.15%.我们一共采用了101个地震，从统计学的角度来说，样本总量可以认为已经足够，因此该统计结果具有统计学意义.

4 结论

本文介绍了计算地球内部应力与应变固体潮的理论公式.基于PREM地球模型，采用高精度的引潮位展开表，根据固体潮理论计算了包含直到球谐展开四阶项的地球内部应力和应变固体潮6个独立分量随时间的潮汐变化，分析了不同深度和纬度处的应力和应变的变化特征，结果可为地球内部的应变和应力测量资料的处理提供固体潮改正值.

作为内部应力固体潮应用于潮汐触发地震研究的实例，本文以发生在我国中部地区的101个地震为例，研究了发震时刻与断层面上的应力潮汐相位之间的相关性.我们发现断层面上的正应力潮汐和应力张量的第一不变量潮汐与地震发震时刻具有明显的相关性，这种相关性可以认为是潮汐触发了地震的证据.该应用实例采用了断层面上的应力潮汐相位作为研究对象，而尽管我们所研究的区域距离海洋较远，但是海潮负荷对潮汐相位仍有较大的影响，因此在计算应力潮汐时必须顾及海潮负荷的贡献，否则将导致统计结果的偏差，从而得出不正确的结论.