大地测量学和地球动力学的研究均涉及到地球在外力(如潮汐)或内力(如地震)作用下的形变问题。为了描述地球表面和内部位移场、重力场对这些力的响应，一般采用h、l、k这样一组无量纲的参数表示，统称为勒夫(Love)数。Love数最早由Love^[1]在研究潮汐问题时提出，采用h和k表示与垂直位移和附加引力位相关的无量纲参数；之后Shida(志田)^[2]提出与水平位移相关的参数l，因此l也称为志田数。其中，h表示在引潮力作用下，地球上某点元的径向位移与假定地球表面覆盖一层海水产生的平衡潮相应点元径向位移之比；l表示地球上某点元的切向位移与假定平衡潮相应点元切向位移之比；k表示由于潮汐形变导致地球物质重新分布所引起的附加引力位与引潮位之比。为了研究地表的负荷问题，学者们^[3-5]又提出负荷Love数的概念，其定义与潮汐Love数类似。潮汐和负荷问题由于满足相同的微分方程，只是地表的边界条件不同，因此可以归为同一类边值问题。由于地表有3个边界条件，因此只要获得3组独立的边界条件，那么任意的边界条件都可以用这3组独立边界条件表示出来，从而对于同类问题中任意边界条件的解都可以由这3组独立边界条件确定的解组合得到。基于此，Saito^[6]提出剪切应力Love数用于研究在地球表面作用一个剪切应力时的变形问题。

对于Love数的计算，大多是由运动方程、介质的本构关系以及泊松方程推导获得线性常微分方程组，结合不同的地表和地球内部的边界条件求解出位移场和外部引力场变化，然后根据它们与Love数的关系求得不同力源激发的地球变形问题的Love数^[3-7]。不同于用数值积分方法计算Love数，本文基于简正模理论推导计算潮汐Love数、负荷Love数、正应力Love数和剪应力Love数的公式。简正模是与地球介质相关的本征运动模式，完全决定于地球介质的物理参数，而地球受激发后的运动都可以用简正模表示出来。求解简正模的微分方程与求解Love数的方程是一样的，唯一的不同是前者是自由边界，而后者的边界是受力的，或者说受激发的，这也是可以用简正模来求Love数的理论基础。

为了研究地震激发的地球自由振荡问题，Gilbert^[8]应用简正模理论推导了大地震引起的位移场变化公式；Chao等^[9]在此基础上，应用简正模理论推导并计算了大地震引起的引力场变化的公式，为根据简正模理论研究Love数提供了参考。本文首先简单介绍利用简正模计算位移场和重力场的相关理论，然后推导由不同力源激发的地球位移场和外部重力场扰动的公式，进一步推导得到静态形变下不同Love数在地表的最终计算公式，利用PREM地球模型计算相应的结果，并与常规方法所得结果进行比较。

1 简正模理论

简正模是多自由度线性力学系统中的集体运动模式，一个多自由度系统可以有多个振动模式，每种模式都有一个固定的特征频率和特征函数，这些振动模式互相独立，可以线性叠加。相应地，地球也存在简正模，当系统受到外力作用时，会发生变形，根据运动方程，可以得到位移随时间变化的函数，随着时间的推移，由于衰减，动态位移场将转变成静态位移场^[8]：

$ {\boldsymbol{u}_\alpha } = \mathop \sum \limits_{n = 1}^\infty {s_{\alpha, n}}\left( {\omega _n^{ - 2}\sum\limits_\beta ^N {s_{\beta, n}^*} \cdot {\boldsymbol{f}_\beta }} \right) $

(1)

式中，s_{α, n}表示第α个质点的n个简正模，其中n可以由一组节点数、阶数和次数表示，N表示质点数，s_{β, n}^*表示第β个质点的第n个简正模，上标*表示共轭，ω_n表示简正模的本征频率，f_β表示系统中第β个质点受到的力。

求解受力激发的静态位移场的核心问题是式(1)括号内的力与简正模的乘积，即简正模组合系数的计算。由于地球是连续体，则式(1)中括号部分的求和将转变成体积分^[8]：

$ \begin{array}{l} \boldsymbol{u}\left( \boldsymbol{r} \right) = \mathop \sum \limits_{n = 1}^\infty {s_n}\left( \boldsymbol{r} \right)\omega _n^{ - 2}\mathop \int{s_n^*\left( \boldsymbol{r} \right) \cdot f\left( \boldsymbol{r} \right){\rm{d}}\mathit{\Omega }}\\ \mathop \sum \limits_n^\infty {{\rm{C}}_f}{s_n}\left( \boldsymbol{r} \right) \end{array} $

(2)

式中，u(r)表示静态位移场，s_n^*(r)表示受力点r处的第n个简正模的共轭，Ω表示地球的体积。

外力引起的地球变形通常伴随着体积的变化和内部物质的重新分布，进而外部的引力场会产生扰动。对于地球外部任一点P(r₀, θ₀, φ₀)，引力位的扰动Δφ(r₀)可以表示为^[9]：

$ \left\{ \begin{array}{l} \Delta \varphi \left( {{r_0}} \right) = G\mathop \sum \limits_{l = 0}^\infty \mathop \sum \limits_{m = - l}^l \frac{{4{\rm{ \mathit{ π} }}}}{{2l + 1}}\frac{1}{{r_0^{l + 1}}}{\mathit{\Gamma _{lm}}}Y_{lm}^*\left( {{\theta _0}, {\phi _0}} \right)\\ {\mathit{\Gamma} _{lm}} = \mathop \int{\rho \left( \boldsymbol{r} \right){r^{l - 1}}\boldsymbol{u}\left( \boldsymbol{r} \right)} \cdot \left[{l{Y_{lm}}\left( {\theta, \phi } \right)\boldsymbol{\hat r} + {\nabla _1}{Y_{lm}}\left( {\theta, \phi } \right)} \right]{\rm{d}}\mathit{\Omega } \end{array} \right. $

(3)

式中，G表示万有引力常数，Y_lm表示归一化的广义球谐函数，l、m分别表示球谐函数的阶数和次数，ρ(r)表示地球密度，由式(1)确定，${\nabla _1} = \boldsymbol{\hat \theta} \partial \theta + \boldsymbol{\hat \phi} {\rm{csc}}\theta \partial \phi $，θ和ϕ分别表示地球内部一点的余纬和经度，$\boldsymbol{\hat r}、\boldsymbol{\hat \theta} 、\boldsymbol{\hat \phi} $表示球坐标系中的3个单位基矢量。这里运用变形前密度和位置分量来计算积分，即拉格朗日法，引力位的扰动Δφ可以通过求解Γ_lm得到。因此，只要根据式(2)求得简正模的系数，代入式(1)和式(3)就可以获得地球内部的位移场和地球外部的引力场变化，进而求得不同的Love数。

2 不同类型Love数的计算

在外力作用下，地表上任一点产生的径向位移u_r、切向位移u_θ及附加引力位Δφ都可以通过Love数和引潮位(以潮汐问题为例，其他问题类似)表示：

$ \begin{array}{l} {\boldsymbol{u}_r} = \mathop \sum \limits_{l = 2}^\infty {h_l}\frac{{{W_l}\left( r \right)}}{g}, {\boldsymbol{u}_\theta } = \mathop \sum \limits_{l = 2}^\infty {l_l}\frac{{\partial {W_l}\left( r \right)}}{{g\partial \theta }}\\ {\boldsymbol{u}_\phi } = \mathop \sum \limits_{l = 2}^\infty {l_l}\frac{{\partial {W_l}\left( r \right)}}{{g{\rm{sin \mathit{ θ} }}\partial \phi }}, \Delta \varphi = \mathop \sum \limits_{l = 2}^\infty {k_l}{W_l}\left( r \right) \end{array} $

(4)

式中，g表示地表重力加速度，h_l、l_l和k_l表示Love数, W_l(r)表示引潮位。

2.1 潮汐Love数

潮汐Love数用来表征地球在日月引潮力的作用下所产生的形变特性。对于球对称的地球，Love数的计算与坐标系的选择没有关系。为了简化运算，在球坐标系下，可以将引潮天体置于极轴上(北极上空)，此时引潮力将与经度无关，因此变形也将与经度无关。对于潮汐引起的地球形变，激发的力源是引潮力，引潮力可以通过引潮位的梯度表示：

$ \begin{array}{l} \boldsymbol{f}\left( \boldsymbol{r} \right) = \rho \left( \boldsymbol{r} \right)\nabla {W_l}\left( \boldsymbol{r} \right) = \\ \rho \left( \boldsymbol{r} \right)\nabla \left[{\frac{{GM}}{R}{{\left( {\frac{\boldsymbol{r}}{R}} \right)}^l}\sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} {Y_{l0}}\left( {{\rm{cos}}\theta } \right)} \right] \end{array} $

(5)

为了方便，用0次规格化的广义球谐函数代替习惯采用的勒让德多项式，把简正模和引潮力的表达式代入式(2)，展开可以得到：

$ \begin{array}{l} C_f^{ln} = \frac{{GM}}{{{R^{l + 1}}}}\sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} l\omega _n^{ - 2}\mathop \smallint \limits_0^a \rho \left( r \right){r^{l + 1}} \cdot \\ \left[{{U_{ln}}\left( r \right) + \left( {l + 1} \right){V_{ln}}\left( r \right)} \right]{\rm{d}}r = \frac{{GM}}{{{R^{l + 1}}}}\sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} {\mathit{\Gamma} _{ln}} \end{array} $

(6)

式中，U_ln(r)、V_ln(r)分别表示第l阶、第n个简正模的垂向位移和水平位移(下标m为0，因此略去)，a表示地球半径。结合式(1)、式(4)求得位移和引力位的变化，化简得到潮汐Love数的表达式：

$\left\{ \begin{array}{l} h_l^T = \mathop \sum \limits_{n = 1}^\infty {U_{ln}}\left( a \right)\frac{g}{{{a^l}}}{\mathit{\Gamma} _{ln}}\\ l_l^T = \mathop \sum \limits_{n = 1}^\infty {V_{ln}}\left( a \right)\frac{g}{{{a^l}}}{\mathit{\Gamma} _{ln}}\\ k_l^T = \mathop \sum \limits_{n = 1}^\infty \frac{{4{\rm{ \mathit{ π} }}G}}{{\left( {2l + 1} \right){a^{2l + 1}}}}\mathit{\Gamma} _{ln}^2 \end{array} \right.$

(7)

2.2 正应力Love数和负荷Love数

地球的变形不仅受到引潮力的影响，还受到地球表面物质迁移而产生的负荷变化的影响，如海潮负荷。研究该问题时，通常引入负荷Love数。地面上的质量一方面对地球施加一个引力，其作用与引潮力相同，因此这部分的响应可用潮汐Love数表示；另一方面该质量对地球也有一个垂直方向的压力作用。

对于压力引起的形变问题，根据平衡方程，等效的单位体力可以表示为应力张量的散度，即f(r) =-▽·T，通过分部积分和高斯定理，式(2)简化为^[8]：

$ {C_f} = - \mathop \int {\boldsymbol{n}} \cdot \left( {T \cdot s_{\rm{n}}^*} \right){\rm{d}}{A_0} + \mathop \int {T}:\varepsilon _n^{\rm{*}}{\rm{d}}\mathit{\Omega } $

(8)

式中，n表示由球心向外的方向向量，A₀表示球面，ε_n^*表示应变张量的共轭。等式左边是作用力所作的功，由于仅考虑静态形变，不考虑温度等其他物理量的变化，根据能量守恒，作用力所作的功将全部转变为弹性应变能，也就是说等式左边在数值上等于应变能。注意到等式右边的第2项是应变能的2倍，所以右边第1项面积分(不带负号)在数值上也等于应变能。因此所求的系数即为等式右边的面积分。

若负荷质量作用于地球的北极点，对于负荷压力引起的变形，激发的力源是正应力，则观测点的正应力Love数也可以用一组无量纲的数h_l^P、l_l^P、k_l^P表示。单位点质量引起的正应力的大小等于垂直方向上压强的大小，方向相反(规定张应力为正)^[6]：

${T_{rr}}\left( a \right) = - \mathop \sum \limits_{l = 0}^\infty \sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} \frac{g}{{{a^2}}}{Y_{l0}}$

(9)

代入式(8)展开可以得到系数的表达式：

$ C_f^{ln} = - \sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} {U_{ln}}\left( a \right)g\omega _n^{ - 2} $

(10)

进一步计算可以得到正应力Love数:

$\left\{ \begin{array}{l} h_l^P = - \mathop \sum \limits_{n = 1}^\infty U_{ln}^2\left( a \right)\frac{{\left( {2l + 1} \right){g^2}a\omega _n^{ - 2}}}{{4{\rm{ \mathit{ π} }}G}}\\ l_l^P = - \mathop \sum \limits_{n = 1}^\infty {U_{ln}}\left( a \right){V_{ln}}\left( a \right)\frac{{\left( {2l + 1} \right){g^2}a\omega _n^{ - 2}}}{{4{\rm{ \mathit{ π} }}G}}\\ k_l^P = - \mathop \sum \limits_{n = 1}^\infty {U_{ln}}\left( a \right)\frac{{g\omega _n^{ - 2}}}{{4{\rm{ \mathit{ π} }}{a^l}}}{\mathit{\Gamma} _{ln}} \end{array} \right.$

(11)

负荷引起地球表面的形变和引力位的变化是表面负荷质量的压力和引潮力共同作用的结果，负荷Love数可以通过二者线性相加表示，即

$\left\{ \begin{array}{l} h_l^L = h_l^T + h_l^P\\ l_l^L = l_l^T + l_l^P\\ k_l^L = k_l^T + k_l^P \end{array} \right.$

(12)

2.3 剪应力Love数

地球受到剪切应力的作用会发生剪切形变，剪切形变也可以通过一组无量纲的参数来描述，即剪应力Love数，此时式(8)中的应力即为剪切应力^[6]：

${T_{r\theta }}\left( a \right) = \sqrt {\frac{{\left( {2l + 1} \right)}}{{4\pi }}} \frac{g}{{{a^2}l\left( {l + 1} \right)}}\frac{{{\rm{d}}{Y_{10}}}}{{{\rm{d}}\theta }}$

(13)

代入式(8)展开可以得到系数的表达式：

$C_f^{ln} = \sqrt {\frac{{\left( {2l + 1} \right)}}{{4{\rm{ \mathit{ π} }}}}} g{V_{ln}}\left( a \right)\omega _n^{ - 2}$

(14)

同样，可以推导出地表剪切力引起的位移和引力位变化公式，求得剪应力Love数:

$ \left\{ \begin{array}{l} h_l^S = \mathop \sum \limits_{n = 1}^\infty {U_{ln}}\left( a \right){V_n}\left( a \right)\frac{{\left( {2l + 1} \right){g^2}a{\omega ^{-2}}}}{{4{\rm{ \mathit{ π} }}G}}\\ l_l^S = \mathop \sum \limits_{n = 1}^\infty V_{ln}^2\left( a \right)\frac{{\left( {2l + 1} \right){g^2}a{\omega ^{-2}}}}{{4{\rm{ \mathit{ π} }}G}}\\ k_l^S = \mathop \sum \limits_{n = 1}^\infty {V_{ln}}\left( a \right)\frac{{g{\omega ^{-2}}}}{{4\pi {a^l}}}{\mathit{\Gamma} _{ln}} \end{array} \right. $

(15)

式中，h_l^S、l_l^S、k_l^S表示地表的剪应力Love数。对于剪切力问题，由于没有质量，因此也没有相应的平衡潮，所以Love数的定义只是采用类似的公式，即假设有单位质量存在。

3 计算结果

为了计算不同的Love数，采用PREM地球模型^[11]，但本文对此模型进行少许修正，即对模型的6 368 ~6 371 km的岩石圈和海洋层进行平均以消除液态的海洋层^[12]。计算选取的简正模周期都是大于5 s的2~60阶的球型简正模，简正模通过Mineos软件计算得到，该程序考虑的对象是球对称、无旋转的地球，适用于长周期低频信号的模拟^[13]。由于不可能对所有的无穷多个简正模求和，因此要进行截断，所以存在求和收敛性问题。

在计算中发现，当球谐阶数为60时，大部分Love数只需要叠加频率由低到高的前20个简正模就已经具有很好的收敛性。但是正应力Love数h_l^P需要叠加500个简正模才有很好的收敛性，剪应力Love数l_l^S需要叠加100个简正模，并且随着阶数的增加，h_l^P和l_l^S需要更多简正模的叠加才会收敛。为了保证简正模结果的精度和可靠性，计算的各阶简正模节点数范围为0~500，即每阶501个简正模叠加。

表 1中给出了Farrell^[5]、Okubo等^[14]与本文数值积分(NI)和简正模(NM)方法计算得到的2阶潮汐、负荷和剪应力Love数的结果，数值积分方法和简正模方法所得结果采用相同的地球模型，该模型与Farrell和Okubo等采用的模型有少许的差别。从表 1中可以看到，简正模方法计算得到的潮汐Love数与数值积分计算结果基本一致。Farrell和Okubo等应用的分别是Gutenberg地球模型与1066A地球模型，这2个模型虽然是球对称、非旋转的，但具体的模型参数与PREM地球模型还是有差别的。然而，尽管采用的模型不同，由不同方法获得的Love数数值结果仍然具有很好的一致性。这是因为模型之间的少许差异不会影响低阶的Love数，但对高阶的Love数影响较大。此外，使用PREM地球模型通过数值积分和简正模2种不同的方法计算的Love数数值上也是比较一致的。正应力Love数虽未给出具体数值，但是负荷Love数是潮汐Love数和正应力Love数线性求和的结果，在负荷Love数和潮汐Love数具有一致性的前提下，正应力Love数计算结果也一定是一致的。

图 1给出了用相同地球模型，采用简正模方法和数值积分方法获得的Love数结果的比较。因为对于潮汐问题，通常只关心较低阶数(如 < 5)的结果，所以未给出潮汐Love数结果，只给出了2~60阶的正应力Love数、负荷Love数和剪应力Love数的结果。

从图 1可以看出，数值积分方法和简正模方法计算的Love数曲线基本重合，这也验证了简正模方法计算Love数的有效性。实际上，这2种方法都是对相同的微分方程进行了数值积分，只不过对地表边界条件的处理方式不一样。前者是设定频率为0，直接求满足边界条件的解；后者是根据不同的本征频率求解本征函数，然后根据边界条件确定组合系数。因此从本质上来说，2种方法是一样的。但通过2种方法的比较，揭示了它们之间的内在联系，这有助于更好地了解地球简正模与Love数之间的关系。

4 结语

基于地球简正模理论，研究了不同力源引起的地球位移场和外部引力场的变化；推导了利用地球简正模计算地球的潮汐Love数、负荷Love数、正应力Love数和剪应力Love数的实用公式；采用PREM地球模型计算了相应的Love数，并与其他学者的结果进行了比较。结果表明，各种方法所得结果之间具有很好的一致性，说明了本文公式的正确性和有效性。

本研究使用的是球对称、非自转、弹性、各向同性的地球模型，没有考虑自转、椭率等因素的影响，因为这些因素对地表Love数计算的影响量级不足1%^[15]。虽然采用的是简单的地球模型，但有助于了解简正模与Love数之间的内在联系。

此外，本文所求的Love数均为静态响应，即激发力源的频率设为0，但对于动态响应问题仍有参考意义，因为对于动态问题只需加入共振项即可。

致谢 感谢Masters提供共享软件Mineos用于地球简正模的理论计算。