負(fù)荷勒夫數(shù)漸近值解析表達(dá)式

周江存，潘爾年

1 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院，大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430077 2 College of Engineering, Cleveland State University, Cleveland, OH 44115, USA

0 引言

全球水循環(huán)是地球表層質(zhì)量遷移的重要形式，它引起全球水質(zhì)量在時(shí)間和空間的重新分布，會(huì)進(jìn)一步導(dǎo)致地球的變形.該變形主要有三種形式：(1)大氣負(fù)荷，主要是大氣質(zhì)量的變化對(duì)地球的引力作用與壓強(qiáng)變化對(duì)地表的負(fù)載作用(Merriam，1992；Van Dam et al.，1994；Sun and Luo，1998；Guo et al.，2004)；(2)潮汐和非潮汐海洋負(fù)荷，主要是海水質(zhì)量的時(shí)空變化對(duì)地球的引力作用和對(duì)地表的負(fù)載作用(Longman，1962，1963；Farrell，1972；Goad，1980；許厚澤等，1982)；(3)陸地水負(fù)荷，主要是陸地水的質(zhì)量變化對(duì)地球的引力作用和對(duì)地表的負(fù)載作用(Zhou et al.，2009).這三種效應(yīng)可以引起地表幾厘米(cm)的位移和數(shù)十微伽(10-8m·s-2)的重力變化.在目前高精度的大地測(cè)量學(xué)研究中，這三種效應(yīng)都是必須要加以考慮的.對(duì)這些效應(yīng)的模擬，目前已有比較成熟的理論，通常采用負(fù)荷質(zhì)量與負(fù)荷格林函數(shù)的褶積積分計(jì)算獲得，其中格林函數(shù)是負(fù)荷勒夫數(shù)或其組合與Legendre函數(shù)(或稱為球函數(shù))或其導(dǎo)數(shù)的乘積的無(wú)窮級(jí)數(shù)求和(Longman，1963；Farrell，1972；Guo et al.，2004).鑒于球函數(shù)有解析的遞推公式，因此計(jì)算格林函數(shù)的基礎(chǔ)就是計(jì)算出負(fù)荷勒夫數(shù).

由于格林函數(shù)是一個(gè)無(wú)窮級(jí)數(shù)的求和，對(duì)于一個(gè)無(wú)窮級(jí)數(shù)求和的問(wèn)題來(lái)說(shuō)，高階勒夫數(shù)的高效計(jì)算以及格林函數(shù)的收斂性如何是數(shù)值計(jì)算的關(guān)鍵.為了計(jì)算高階的負(fù)荷勒夫數(shù)，Longman(1963)提出了對(duì)變量進(jìn)行無(wú)量綱化的處理；汪漢勝等(1996)通過(guò)觀察各變量的變化特征提出了r冪因子法；特別地，Pan等(2015)基于合理的假設(shè)提出了計(jì)算負(fù)荷勒夫數(shù)的解析解方法；Chen等(2018)基于解析解方法提出了一種稱為剛度矩陣的傳播矩陣方法來(lái)計(jì)算超高階負(fù)荷勒夫數(shù)，克服傳統(tǒng)傳播矩陣方法的溢出問(wèn)題；Zhou等(2019a)運(yùn)用DVP(dual variable and position)傳播矩陣法獲得了任意球諧階數(shù)的勒夫數(shù)，并提供了相應(yīng)的MATLAB程序計(jì)算不同類型的勒夫數(shù)(包括潮汐、負(fù)荷、剪切和位錯(cuò)).實(shí)際上，這些表征地球在不同力源(如引潮力、地表正應(yīng)力或剪切應(yīng)力、以及內(nèi)部斷層的等效體力等)作用下產(chǎn)生變形的勒夫數(shù)之間并不完全獨(dú)立，它們之間的關(guān)系可參考潘爾年和丁中一(1986)和Okubo(1993).

為了保證格林函數(shù)的收斂，F(xiàn)arrell(1972)給出了負(fù)荷勒夫數(shù)的漸近值，即當(dāng)球諧階數(shù)n趨于無(wú)窮大時(shí)勒夫數(shù)的值，然后通過(guò)Kummer變換可有效降低截?cái)嚯A數(shù)(約10000階)并加速格林函數(shù)的收斂，從而為能夠準(zhǔn)確地模擬負(fù)荷引起的地球形變奠定了基礎(chǔ).可見(jiàn)，負(fù)荷勒夫數(shù)漸近值對(duì)負(fù)荷格林函數(shù)來(lái)說(shuō)是至關(guān)重要的.除了Farrell(1972)之外，Okubo和Endo(1985)、Guo等(2004)也對(duì)如何計(jì)算負(fù)荷勒夫數(shù)漸近值給出了不一樣的理論計(jì)算方法.我們最近提出了計(jì)算位錯(cuò)勒夫數(shù)的漸近值的方法，在忽略自引力的情形下，漸近值可根據(jù)傳播矩陣的特性很容易地求出(Zhou et al.，2021a).位錯(cuò)問(wèn)題和負(fù)荷問(wèn)題的差別在于邊界條件的不同，因此該方法可直接應(yīng)用于求取負(fù)荷勒夫數(shù)的漸近值.

本文首先對(duì)目前已有的上述三種計(jì)算漸近值的方法進(jìn)行簡(jiǎn)單的回顧，然后著重介紹一種新的能夠獲取負(fù)荷勒夫數(shù)漸近值解析表達(dá)式的方法.最后將我們的方法與已有的計(jì)算方法進(jìn)行了比較，并給出了定量的數(shù)值結(jié)果展示新的解析表達(dá)式在表達(dá)負(fù)荷勒夫數(shù)漸近特征方面的優(yōu)勢(shì).

1 已有方法介紹

對(duì)于地球受力變形，我們采用如下的描述地球變形的微分方程組邊值問(wèn)題.其中，微分方程組為

(1)

其中λ，μ是拉梅常數(shù)，ρ是密度，g是重力加速度(習(xí)慣上稱之為重力)，它們都是球坐標(biāo)系的徑向分量r的函數(shù)，G是萬(wàn)有引力常數(shù)，UL、UM、Φ、TL、TM分別是垂向位移、水平位移、附加引力位、垂向正應(yīng)力、垂向剪應(yīng)力的球諧展開(kāi)的第n階的系數(shù)(Pan et al.，2015)，Q滿足

(2)

因此，這6個(gè)變量要分別對(duì)不同的球諧階數(shù)n求解.該方程是研究地球在不同力源作用下發(fā)生靜態(tài)變形的基礎(chǔ)微分方程組，對(duì)于不同的力源需采用相應(yīng)的邊界條件(Love，1911；Farrell，1972；Okubo and Endo，1985；Sun and Okubo，1993).

對(duì)于地表質(zhì)量負(fù)荷問(wèn)題，通常將負(fù)荷看成點(diǎn)源，并且將該點(diǎn)源放置于北極點(diǎn).由于對(duì)稱性，在負(fù)荷作用下的變形與經(jīng)度無(wú)關(guān)，即球諧展開(kāi)的次數(shù)m=0，此時(shí)球函數(shù)退化為L(zhǎng)egendre多項(xiàng)式.

下面介紹邊界條件.為方便起見(jiàn)，我們定義兩個(gè)矢量

(3)

在地心處正則，即

U(r=0)=0.

(4)

在內(nèi)部邊界連續(xù)，以及在地表

(5)

其中，ga是地表重力，a是地球半徑.在推導(dǎo)該式中的邊界值時(shí)，認(rèn)為負(fù)荷點(diǎn)源的質(zhì)量M滿足單位引力位假定(即GM/a=1，可參考Sun and Sj?berg，1999)，也就是說(shuō)負(fù)荷質(zhì)量M=a/G=me/aga，其中me為地球質(zhì)量.也有一些研究采用地球質(zhì)量假定，即M=me(如Guo et al.，2004).不管用哪種定義，在計(jì)算格林函數(shù)時(shí)都要除以負(fù)荷質(zhì)量以得到單位質(zhì)量負(fù)荷的結(jié)果，因此只要推導(dǎo)過(guò)程自洽都不影響最終的結(jié)果.

負(fù)荷勒夫數(shù)定義為

(6)

其中，帶下標(biāo)n的負(fù)荷勒夫數(shù)hn、ln、kn表示對(duì)應(yīng)于球諧展開(kāi)的n階的地表結(jié)果.因?yàn)?1)式是關(guān)于各個(gè)物理量的球諧展開(kāi)系數(shù)的，所以不同的球諧階數(shù)對(duì)應(yīng)不同的勒夫數(shù).相應(yīng)的負(fù)荷格林函數(shù)為

(7)

其中，u，v和φ分別為垂直位移、切向位移和附加引力位，Pn是Legendre多項(xiàng)式，θ是計(jì)算點(diǎn)余緯.可以看到格林函數(shù)是一個(gè)無(wú)窮級(jí)數(shù)求和，需要截?cái)嗵幚?但是對(duì)于負(fù)荷問(wèn)題需要截?cái)嗟胶芨叩碾A數(shù)，這給計(jì)算帶來(lái)很大的時(shí)間成本.為此Farrell(1972)根據(jù)Boussinesq問(wèn)題的解(見(jiàn)下文)給出了負(fù)荷勒夫數(shù)的漸近值，從而可有效降低截?cái)嗟碾A數(shù)，并加速格林函數(shù)計(jì)算的收斂.

應(yīng)用Kummer變換，(7)式改寫(xiě)為

(8)

其中，帶下標(biāo)∞的負(fù)荷勒夫數(shù)h∞、l∞、k∞表示相應(yīng)的漸近值.由于負(fù)荷勒夫數(shù)漸近值具有嚴(yán)格解析表達(dá)式，右端的第二項(xiàng)通常具有用三角函數(shù)表示的嚴(yán)格解析形式(Farrell，1972；Zhou et al.，2019b)，而第一項(xiàng)求和的收斂要比(7)式中的求和快得多，因此截?cái)嚯A數(shù)N可以取得足夠小(對(duì)于位移與附加引力位，N=3000；對(duì)于傾斜或應(yīng)變N=10000足以保證計(jì)算的收斂).可見(jiàn)負(fù)荷勒夫數(shù)漸近值在計(jì)算負(fù)荷格林函數(shù)時(shí)發(fā)揮著至關(guān)重要的作用.

對(duì)于當(dāng)n很大時(shí)的情形，我們可以忽略液核的影響，因?yàn)橐汉说拇嬖谥挥绊戄^低階的負(fù)荷勒夫數(shù)，所以可以采用純固體的地球模型.實(shí)際上，正如下文將要討論的那樣，只有地球最外層的參數(shù)才會(huì)影響負(fù)荷勒夫數(shù)的漸近值，這也是為什么可以直接用均勻球(采用地球最外層結(jié)構(gòu)參數(shù)和地表重力值，如Okubo，1988)或僅僅采用最外層參數(shù)(如郭俊義，2000；Guo et al.，2004)來(lái)解算負(fù)荷勒夫數(shù)的漸近值.

下面分別介紹目前已有的三種求解負(fù)荷勒夫數(shù)漸近值的方法.

1.1 Farrell方法

在球坐標(biāo)系下，當(dāng)角距離(θ)非常小的時(shí)候，問(wèn)題就等效于平面的Boussinesq問(wèn)題，因?yàn)榇藭r(shí)地球的曲率效應(yīng)可忽略(Farrell，1972).對(duì)于Boussinesq問(wèn)題，采用柱坐標(biāo)及柱函數(shù)(Bessel函數(shù))系(用z，r，φ表示，將負(fù)荷點(diǎn)源放置于坐標(biāo)原點(diǎn)，由于對(duì)稱性，變形與φ無(wú)關(guān).注意此時(shí)的r表示距離z軸的距離)，所有變量都利用Hankel變換(潘爾年等，1986)，即

(9)

其中，ξ表示波數(shù).為了簡(jiǎn)便，仍用與球諧展開(kāi)系數(shù)相同的符號(hào)表示變換后的系數(shù).在此情形下，球函數(shù)趨近于柱函數(shù)，即

(10)

并且

(11)

因此，(7)式的求和可以表示成如下的積分：

(12)

根據(jù)(11)式，(12)式變成

(13)

對(duì)于Boussinesq問(wèn)題，一個(gè)點(diǎn)負(fù)荷作用下的解為(在Hankel變換域)

(14)

其中ρa(bǔ)表示地球模型最上層的密度，即地表密度，〈ρ〉表示地球平均密度.將(14)式代入(9)式，并與(13)式進(jìn)行比較，可得

(15)

在此過(guò)程中，要運(yùn)用到(11)式中的n與ξ的近似關(guān)系n=aξ，因此其結(jié)果精度為O(1/n).我們稱(15)式的結(jié)果為一級(jí)漸近值，以區(qū)別于Guo等(2004)給出的更高一級(jí)的漸近值.

1.2 Okubo和Endo方法

Okubo和Endo(1985)、Okubo(1988)推導(dǎo)負(fù)荷勒夫數(shù)漸近值的方法是基于Love(1911)給出的均勻球變形的通解表達(dá)式，用X=[ULUMΦTLTMQ]T表示由6個(gè)變量表示的這個(gè)通解矢量，其中上標(biāo)T表示轉(zhuǎn)置，其解可表示為

(16)

其中，c是由6個(gè)常數(shù)標(biāo)量組成的待定系數(shù)矢量，且

[D(r)]=

(17)

(18)

(19)

式中

+[2f±+n(n+1)]zn(k±r)},

(20)

(21)

(22)

(23)

α和β分別為P波和S波速度，jn(x)是球Bessel函數(shù)，n為球諧階數(shù).利用邊界條件(5)式可確定c，從而求得地表的位移與附加引力位，即

U(a)=D(a)[E(a)]-1T(a).

(24)

由于高階負(fù)荷勒夫數(shù)只與最上層的結(jié)構(gòu)有關(guān)，而如果采用密度為地表密度的均勻球，則該均勻球表面的重力將不等于真實(shí)地球的地表重力.Okubo(1988)指出，可用真實(shí)地球的地表重力值代替均勻球表面的重力，這樣得到的負(fù)荷勒夫數(shù)漸近值就是真實(shí)地球的結(jié)果.因此地球模型的參數(shù)代入(24)式，并將各量展開(kāi)成n的級(jí)數(shù)后，就可得到如(15)式所示的負(fù)荷勒夫數(shù)漸近值.

1.3 Guo方法

郭俊義等(郭俊義，2000；Guo et al.，2004)根據(jù)地表邊界條件考察了每個(gè)變量與球諧階數(shù)n的量級(jí)關(guān)系，得到

(25)

例如，TL隨著n的增大而增大，與n呈線性關(guān)系，而UL不隨n的增大而增大.然后通過(guò)變量代換

(26)

使得所有的新變量zi具有關(guān)于n相同的量級(jí)，并在微分方程中忽略高階小量，從而得到

(27)

該微分方程組存在解析解，可以求出負(fù)荷勒夫數(shù)h和l的漸近值.但是k的漸近值無(wú)法求出，原因是舍去的項(xiàng)太多了.因此他們繼續(xù)保留高一階小量，并求解一個(gè)非齊次的微分方程組，從而求出了k的漸近值.

在此過(guò)程中，Guo等(2004)未介紹如何求得(27)式的解析解，而是直接給出了結(jié)果.因?yàn)樵摲匠探M實(shí)際上是變系數(shù)的，盡管將n/r移到了方程組的左邊.所以其是否存在解析解，并不是那么直觀地能夠看出來(lái).實(shí)際上，因?yàn)?27)式是歐拉型，所以通過(guò)變量代換，

r=eξ,

(28)

可將該方程組化成常系數(shù)的微分方程組，從而有解析解.Guo等(2004)以一級(jí)漸近值為基礎(chǔ)，求解關(guān)于高一級(jí)小量的非齊次微分方程組，給出了更高一級(jí)的漸近值.但是對(duì)于負(fù)荷問(wèn)題，更高一級(jí)的漸近值在計(jì)算負(fù)荷格林函數(shù)時(shí)是沒(méi)有必要的，因?yàn)槔?15)式的漸近值已經(jīng)能夠很高效地獲得收斂的格林函數(shù)(參考下文的圖2).

圖2 負(fù)荷格林函數(shù)隨角距離的變化(規(guī)格化系數(shù)為1012aθ)“asym.with g”表示采用考慮自引力的50000階數(shù)值結(jié)果作為漸近值，“asym.without g”表示采用不考慮自引力的50000階數(shù)值結(jié)果作為漸近值.

2 解析解方法

本文介紹的推導(dǎo)負(fù)荷勒夫數(shù)漸近值的新方法是以Pan等(2015)提出的解析法和DVP傳播矩陣方法為基礎(chǔ).其中解析法的思想是使得地球變形的通解具有解析的形式，而DVP傳播矩陣法是建立高階地表負(fù)荷勒夫數(shù)與地表邊界條件之間的直接聯(lián)系.下面首先簡(jiǎn)單介紹解析法和DVP傳播矩陣法，具體內(nèi)容可參考(Pan et al.，2015)，然后給出推導(dǎo)負(fù)荷勒夫數(shù)漸近值的過(guò)程.

2.1 解析解構(gòu)造

Pan等(2015)提出了求解(1)式的解析法.對(duì)于一個(gè)分層地球模型，(1)式是沒(méi)有解析解的，但是通過(guò)合理的假設(shè)可以構(gòu)造出近似的解析解.首先，在地核中假設(shè)每層密度是常數(shù)，重力是r的線性函數(shù); 在地幔和地殼中，重力是常數(shù)，密度是r的反函數(shù)，如在第j層中，重力和密度滿足

(29)

和

(30)

其中，b和c是常數(shù)，rj-1表示第j層的下邊界，rj表示第j層的上邊界，rc表示核幔邊界的半徑，ρ0j是第j層的平均密度，r0j是使該層質(zhì)量守恒的等效半徑，它滿足

(31)

以上的假設(shè)是否合理，一方面可以考察假設(shè)后的模型與原模型的近似程度，另一方面取決于計(jì)算結(jié)果是否準(zhǔn)確.Pan等(2015)的數(shù)值結(jié)果已經(jīng)證明了結(jié)果的可靠性以及假設(shè)的合理性，實(shí)際上上述假設(shè)正是基于PREM模型的參數(shù)分布特征作出的，假設(shè)后的模型與原模型的符合程度是非常高的，差異僅為0.6%.并且這種方法也成功運(yùn)用到求解位錯(cuò)以及地表熱負(fù)荷引起的變形問(wèn)題中，也被證明是非常高效的(Zhou et al.，2019a，b，2020，2021a,b).

此外，作一個(gè)變量代換

(32)

顯然，對(duì)于j層的下邊界有ξ=0.再以rTL，rTM和rQ為新變量，則(1)式變?yōu)?/p>

(33)

如果再假設(shè)第j層內(nèi)的彈性模量為常數(shù)(目前的研究基本都是這么做的)，那么可以看出上述微分方程組的系數(shù)都是常數(shù)，這樣就存在解析解.為簡(jiǎn)便計(jì)，把(33)式表示成矢量和矩陣的形式，即

(34)

其中，

(35)

為了方便，我們用球諧階數(shù)n對(duì)變量作了規(guī)格化處理使各個(gè)變量量級(jí)相當(dāng).相應(yīng)的矩陣A可從(33)式導(dǎo)出(系數(shù)作相應(yīng)的調(diào)整，乘上與n有關(guān)的系數(shù)).此時(shí)，通解可以表示成(如在地球模型的第j層)

(36)

(36)式右邊第一個(gè)矩陣是系數(shù)矩陣A的特征向量組成的矩陣，下標(biāo)j表示它由第j層的參數(shù)計(jì)算獲得，第二個(gè)矩陣中的s1和s2分別是由具有正實(shí)部(上標(biāo)+)和負(fù)實(shí)部(上標(biāo)-)的三個(gè)特征值組成，即

(37)

從而，由傳播矩陣法建立第j層上下邊界處(ξ=ξj和ξ=0)的兩組解之間的關(guān)系：

(38)

其中上標(biāo)u和l分別表示第j層的上下邊界，上標(biāo)-1表示矩陣求逆.然而，Pan等(2015)的數(shù)值結(jié)果表明，由于數(shù)值溢出，普通的傳播矩陣方法，即(38)式，無(wú)法正確算出高階(如n>6000)結(jié)果.下面要介紹的DVP傳播矩陣法可以很好地解決這個(gè)問(wèn)題.

2.2 DVP傳播矩陣法

該方法的核心思想是避免導(dǎo)致溢出的大數(shù)產(chǎn)生，其技巧僅在于對(duì)普通的傳播矩陣法進(jìn)行小的改動(dòng)，即

1)將(36)式改寫(xiě)為

由于ξj是常數(shù)，所以最終只是改變了待定系數(shù)，從(36)式中的c1變成(39)式的c3，并不會(huì)改變最終的解.

2)在第j層中，由上下邊界的解我們可以得到與(38)式不一樣的關(guān)系：

(40)

其中，

(41)

需要指出的是，這種傳播矩陣方法的特別之處在于變量的位置關(guān)系與傳統(tǒng)的傳播矩陣方法是不一樣的(對(duì)比(38)與(40)兩式).從(41)式可以看出，在傳播矩陣的計(jì)算中不涉及大的正數(shù)的指數(shù)，從而避免大數(shù)的產(chǎn)生，解決了溢出的問(wèn)題.

2.3 漸近值推導(dǎo)

基于上面的方法原理，我們就可以推導(dǎo)出負(fù)荷勒夫數(shù)的漸近值.特別地，當(dāng)球諧階數(shù)n很大時(shí)，由于特征值接近于±n，所以(41)式中的〈-s1ξj〉和〈s2ξj〉都接近于0，當(dāng)n為無(wú)窮大時(shí)，它們就等于0.因此(41)式變成

(42)

由(40)式進(jìn)一步得到

(43)

這說(shuō)明，當(dāng)球諧階數(shù)很大時(shí)，U和T是相關(guān)的，這種相關(guān)性會(huì)導(dǎo)致傳統(tǒng)的龍格-庫(kù)塔數(shù)值積分法無(wú)法準(zhǔn)確計(jì)算出非常高階的負(fù)荷勒夫數(shù); 然而這種相關(guān)性可以使我們能夠直接求得表面的負(fù)荷勒夫數(shù)漸近值.對(duì)于地球的最外層(設(shè)為第p層)，我們可由(43)式直接得出

(44)

其中，方括號(hào)的下標(biāo)p表示方括號(hào)內(nèi)的矩陣由第p層的參數(shù)求得.這說(shuō)明，只要我們知道最外層的系數(shù)矩陣A的特征向量，就可以得到當(dāng)n趨于無(wú)窮大時(shí)的U，從而得到負(fù)荷勒夫數(shù)的漸近值.

通常情況下，矩陣A的特征向量是無(wú)法解析地表達(dá)出來(lái)的.但是我們知道，對(duì)于地表質(zhì)量負(fù)荷問(wèn)題，在計(jì)算高階的負(fù)荷勒夫數(shù)時(shí)，地球的自引力效應(yīng)是可以忽略的，球諧階數(shù)越高，自引力效應(yīng)的影響越不明顯(見(jiàn)下文圖1).實(shí)際上Farrell(1972)從Boussinesq近似得到的漸近值也是忽略了自引力影響的.

圖1 負(fù)荷勒夫數(shù)隨球諧階數(shù)的變化

對(duì)于忽略自引力的情形，在(33)式關(guān)于應(yīng)力的微分方程中就可以舍去與密度或重力有關(guān)的項(xiàng)，從而位移和應(yīng)力是與附加引力位無(wú)關(guān)的，因此可以首先解出(如Wason and Singh，1972；Fang et al.，2014)，然而附加引力位是與位移相關(guān)的，所以繼續(xù)求附加引力位時(shí)需要求解一個(gè)非齊次的微分方程組，這也是1.3節(jié)中Guo等(2004)求負(fù)荷勒夫數(shù)k的漸近值采用的方法，但這反而增加了難度.實(shí)際上，我們可以直接運(yùn)用2.1節(jié)的方法，此時(shí)矩陣A的特征向量可以很容易地求出(如通過(guò)MATLAB的符號(hào)運(yùn)算，具體表達(dá)式可參考Zhou et al.，2021a).

因此根據(jù)矩陣A的特征向量，以及(5)式和(44)式可以很容易求得

(45)

(45)式是高階負(fù)荷勒夫數(shù)的嚴(yán)格解析表達(dá)式，給出了負(fù)荷勒夫數(shù)與球諧階數(shù)n的關(guān)系.特別地，當(dāng)n趨于無(wú)窮大時(shí), 取極限就得到

(46)

由于我們是求負(fù)荷勒夫數(shù)在n趨于無(wú)窮大時(shí)的漸近值，我們可以將第p層設(shè)得足夠薄，因此有

ρ0pr0p=ρa(bǔ)a,

(47)

即地表密度與地球半徑的乘積.這樣就可以得到與(15)式完全一致的結(jié)果.實(shí)際上對(duì)于最外層來(lái)說(shuō)，地球模型通常是均勻的，因此ρ0p就是這個(gè)平均值，而由(31)式計(jì)算得到的r0p實(shí)際已經(jīng)非常接近地球半徑a了(對(duì)于PREM模型來(lái)說(shuō)，誤差小于萬(wàn)分之三)，就算直接使用r0p也可以確保k∞的準(zhǔn)確性.

2.4 數(shù)值結(jié)果

圖1給出了根據(jù)(45)式確定的負(fù)荷勒夫數(shù)隨球諧階數(shù)的變化, 并與Guo等(2004)考慮高一階小量的漸近值進(jìn)行了比較.其中地球模型采用PREM(Dziewonski and Anderson, 1981)，其最外層的海洋層替換為固體層，其參數(shù)設(shè)置如下：λ=3.43×1010Pa，μ=2.66×1010Pa，ρ=2600 kg·m-3.從圖中可以看出，隨著球諧階數(shù)n的增大，負(fù)荷勒夫數(shù)趨近于一個(gè)常數(shù).由于我們推導(dǎo)高階負(fù)荷勒夫數(shù)解析公式時(shí)未考慮自引力的影響，因此所得結(jié)果在高階部分與不考慮自引力的數(shù)值結(jié)果比較接近，在n=5000左右時(shí)二者就已經(jīng)比較一致；但是在低階部分與考慮自引力的數(shù)值結(jié)果仍然存在很大的差異.通常計(jì)算負(fù)荷格林函數(shù)時(shí)只要考慮10000階以下的負(fù)荷勒夫數(shù)，因此在研究負(fù)荷形變時(shí)自引力的影響一般是需要加以考慮的.但是從圖中可以看出，隨著n的增大，自引力的影響越來(lái)越小，因此在研究負(fù)荷勒夫數(shù)漸近特征時(shí)，自引力是可以忽略的，這可從圖1中不同結(jié)果在n很大時(shí)具有的一致性得到驗(yàn)證.

此外，我們也將本文給出的漸近值與Guo等(2004)的結(jié)果進(jìn)行了比較.結(jié)果說(shuō)明，在n<20000左右時(shí)，Guo等(2004)的漸近值與考慮自引力的數(shù)值結(jié)果更加接近，這是因?yàn)樗麄兺茖?dǎo)漸近值時(shí)部分地考慮了自引力的影響；但是對(duì)于n>20000的部分，我們的結(jié)果不管是和考慮自引力還是和不考慮自引力的數(shù)值結(jié)果都更加接近，對(duì)于位移的勒夫數(shù)，這種情況更加明顯.

為了顯示漸近值在計(jì)算格林函數(shù)中的作用，作為例子，我們根據(jù)(8)式計(jì)算了垂直位移、水平位移和附加引力位的格林函數(shù)隨角距離θ的變化關(guān)系，結(jié)果如圖2所示.為了更好的顯示結(jié)果，格林函數(shù)值都乘上了一個(gè)規(guī)格化系數(shù)1012aθ.首先，通過(guò)對(duì)比采用不同漸近值獲得的格林函數(shù)，發(fā)現(xiàn)漸近值之間的微小差異不會(huì)影響格林函數(shù)的計(jì)算.并且我們發(fā)現(xiàn)負(fù)荷勒夫數(shù)漸近值在附加引力位的計(jì)算中并不是必需的，不采用Kummer變換同樣可以獲得收斂的格林函數(shù).圖2也說(shuō)明了如果不采用Kummer變換，即使截?cái)嚯A數(shù)N=50000，位移格林函數(shù)仍然是不收斂的，說(shuō)明漸近值對(duì)于位移格林函數(shù)的計(jì)算來(lái)說(shuō)是必不可少的，這從另一個(gè)側(cè)面證明了漸近值的重要性.

3 討論

從2.3節(jié)的推導(dǎo)我們可以看出，負(fù)荷勒夫數(shù)的漸近值只和地球模型地表的參數(shù)有關(guān)，而不論地球內(nèi)部的結(jié)構(gòu)如何.其實(shí)，這一點(diǎn)早就為我們所熟知，即高階的負(fù)荷勒夫數(shù)與地殼和上地幔介質(zhì)密切相關(guān)，因此采用龍格-庫(kù)塔方法計(jì)算高階負(fù)荷勒夫數(shù)時(shí)可以直接從地幔開(kāi)始積分.然而，通過(guò)我們的推導(dǎo)可以從本質(zhì)上了解其原因：地球每個(gè)分層內(nèi)的U和T之間的相關(guān)性隨著球諧階數(shù)的增大而增加，從而地表的負(fù)荷勒夫數(shù)可由地表的邊界條件直接確定.

相比于Okubo和Endo方法，我們給出的新方法不需要作均勻球的假設(shè)(或基于均勻球模型開(kāi)展相關(guān)推導(dǎo))，不需要用真實(shí)地球的地表重力代替均勻球表面的重力.目前基于均勻球的地球受力變形的數(shù)值模擬必須要替換地表重力才能得到正確的結(jié)果，如Sun(2004)對(duì)位錯(cuò)引起的地表位移場(chǎng)變化的模擬.與Guo等(2004)的方法一樣，我們的方法實(shí)際上也只需要考慮地球模型的最外層即可.

從Guo等(2004)給出的負(fù)荷勒夫數(shù)漸近值的更高階的結(jié)果可以看出，在位移(h和l)結(jié)果中耦合進(jìn)了重力的影響.但是將(44)式展開(kāi)成n的羅朗級(jí)數(shù)后高階解中沒(méi)有重力項(xiàng)，原因在于我們忽略了自引力的影響.對(duì)于地表負(fù)荷問(wèn)題，忽略自引力不會(huì)影響負(fù)荷勒夫數(shù)漸近值的一級(jí)結(jié)果，即O(1)量級(jí)的常數(shù)項(xiàng)，但是會(huì)影響高一級(jí)的結(jié)果，但是這種影響隨著球諧階數(shù)的增大越來(lái)越弱，并且高一級(jí)的結(jié)果在格林函數(shù)的計(jì)算中是不必要的.在2.3節(jié)中我們忽略自引力的目的是為了獲得矩陣A的特征向量的解析形式，從而獲得漸近值的解析形式，這樣就可以運(yùn)用Kummer變換.Kummer變換的關(guān)鍵有兩點(diǎn)：一方面，漸近值必須與實(shí)際的勒夫數(shù)比較接近，從而它們的差是一個(gè)高階小量，并且這個(gè)差異隨著球諧階數(shù)的增大而減小，這樣(8)式右端的第一個(gè)級(jí)數(shù)求和能夠很快地收斂；另一方面，漸近值需要具有嚴(yán)格解析式，從而(8)式右端的第二個(gè)級(jí)數(shù)求和有簡(jiǎn)單的解析式，能夠用簡(jiǎn)單的函數(shù)表示.Chen等(2018)采用的剛度矩陣法或者2.2節(jié)介紹的DVP都能夠計(jì)算任意球諧階數(shù)的負(fù)荷勒夫數(shù)，但是要獲得收斂的負(fù)荷格林函數(shù)，具有解析表達(dá)式的漸近值對(duì)于簡(jiǎn)化負(fù)荷格林函數(shù)的計(jì)算仍是非常必要的，因?yàn)榻財(cái)嗟胶芨叩那蛑C階數(shù)是很費(fèi)時(shí)低效的.采用擬合的方法確定當(dāng)n非常大時(shí)的勒夫數(shù)隨n變化的解析式(即漸近值)也是一個(gè)非常有效的方法(如Pan et al.，2015；Zhou et al.，2019b).

本文的方法除了采用Pan等(2015)的解析解方法獲得常微分方程組的系數(shù)矩陣A的特征向量解析表達(dá)式外，另一個(gè)關(guān)鍵之處是采用了DVP傳播矩陣方法，(42)式顯示了其獨(dú)特的性質(zhì)，階數(shù)n越大，其計(jì)算方面的優(yōu)越性越明顯.DVP方法只是在傳統(tǒng)的傳播矩陣方法基礎(chǔ)上作了小的改動(dòng)，變換了變量的位置，就解決了數(shù)值溢出的問(wèn)題，并且計(jì)算也非常的高效.

最后需要強(qiáng)調(diào)的是，從表面上看，本文的方法涉及構(gòu)造解析解的理論方法和DVP傳播矩陣方法.實(shí)際上推導(dǎo)負(fù)荷勒夫數(shù)的過(guò)程非常簡(jiǎn)單，可以分成以下兩個(gè)步驟：首先計(jì)算常微分方程組的系數(shù)矩陣A的特征向量從而得到(44)式右端方括號(hào)中的矩陣；然后根據(jù)地表邊界條件由(44)式通過(guò)矩陣運(yùn)算直接得到需要的漸近值的解析表達(dá)式.該過(guò)程只涉及矩陣特征向量的計(jì)算及相關(guān)的矩陣運(yùn)算，原理簡(jiǎn)單，也非常容易編程實(shí)現(xiàn)(Pan，2019；Zhou et al.，2019a).

4 結(jié)論

本文基于求解地球變形的解析解方法和DVP傳播矩陣方法，導(dǎo)出了高階負(fù)荷勒夫數(shù)的解析表達(dá)式，進(jìn)而確定了負(fù)荷勒夫數(shù)的漸近值.漸近值的推導(dǎo)只與地球模型的最外層有關(guān)，推導(dǎo)過(guò)程只涉及矩陣特征向量計(jì)算和簡(jiǎn)單的矩陣運(yùn)算.本文所得結(jié)果與已有方法給出的結(jié)果是一致的，并且更加準(zhǔn)確地描述了負(fù)荷勒夫數(shù)的漸近特征.

利用DVP傳播矩陣法對(duì)于n非常大時(shí)所具備的特點(diǎn)，我們不但可以看出變量U和T的關(guān)聯(lián)性，這也是導(dǎo)致傳統(tǒng)的龍格-庫(kù)塔數(shù)值積分法失效的本質(zhì)原因，還可以直接根據(jù)地表邊界條件求出負(fù)荷勒夫數(shù)的漸近值.這說(shuō)明DVP傳播矩陣法在數(shù)值計(jì)算方面有其獨(dú)特的優(yōu)勢(shì)，可以廣泛應(yīng)用于分層地球模型的變形研究中，如負(fù)荷形變、地震形變等靜態(tài)變形問(wèn)題，以及地球簡(jiǎn)正?；蚝铣傻卣饒D等動(dòng)態(tài)變形問(wèn)題中.

致謝本文第二作者于1982年至1984年在北京大學(xué)地質(zhì)系師承王仁教授攻讀碩士學(xué)位，學(xué)位論文題目為《成層地球模型對(duì)體勢(shì)力及表載的響應(yīng)》.該論文利用球諧函數(shù)與傳播矩陣方法研究分層球形地球在日月引潮力和地表負(fù)載作用下的位移場(chǎng)和應(yīng)力場(chǎng)變化，并采用數(shù)值積分方法計(jì)算勒夫數(shù).在王仁教授的教誨下，科研能力得到了很大的提升，為后來(lái)從事力學(xué)和數(shù)學(xué)相關(guān)的研究奠定了堅(jiān)實(shí)的基礎(chǔ)，對(duì)現(xiàn)在的研究工作仍有很大的幫助.