負荷勒夫數漸近值解析表達式

周江存，潘爾年

1 中國科學院精密測量科學與技術創新研究院，大地測量與地球動力學國家重點實驗室，武漢 430077 2 College of Engineering, Cleveland State University, Cleveland, OH 44115, USA

0 引言

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

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

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

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

1 已有方法介紹

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

(1)

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

(2)

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

對于地表質量負荷問題，通常將負荷看成點源，并且將該點源放置于北極點.由于對稱性，在負荷作用下的變形與經度無關，即球諧展開的次數m=0，此時球函數退化為Legendre多項式.

下面介紹邊界條件.為方便起見，我們定義兩個矢量

(3)

在地心處正則，即

U(r=0)=0.

(4)

在內部邊界連續，以及在地表

(5)

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

負荷勒夫數定義為

(6)

其中，帶下標n的負荷勒夫數hn、ln、kn表示對應于球諧展開的n階的地表結果.因為(1)式是關于各個物理量的球諧展開系數的，所以不同的球諧階數對應不同的勒夫數.相應的負荷格林函數為

(7)

其中，u，v和φ分別為垂直位移、切向位移和附加引力位，Pn是Legendre多項式，θ是計算點余緯.可以看到格林函數是一個無窮級數求和，需要截斷處理.但是對于負荷問題需要截斷到很高的階數，這給計算帶來很大的時間成本.為此Farrell(1972)根據Boussinesq問題的解(見下文)給出了負荷勒夫數的漸近值，從而可有效降低截斷的階數，并加速格林函數計算的收斂.

應用Kummer變換，(7)式改寫為

(8)

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

對于當n很大時的情形，我們可以忽略液核的影響，因為液核的存在只影響較低階的負荷勒夫數，所以可以采用純固體的地球模型.實際上，正如下文將要討論的那樣，只有地球最外層的參數才會影響負荷勒夫數的漸近值，這也是為什么可以直接用均勻球(采用地球最外層結構參數和地表重力值，如Okubo，1988)或僅僅采用最外層參數(如郭俊義，2000；Guo et al.，2004)來解算負荷勒夫數的漸近值.

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

1.1 Farrell方法

在球坐標系下，當角距離(θ)非常小的時候，問題就等效于平面的Boussinesq問題，因為此時地球的曲率效應可忽略(Farrell，1972).對于Boussinesq問題，采用柱坐標及柱函數(Bessel函數)系(用z，r，φ表示，將負荷點源放置于坐標原點，由于對稱性，變形與φ無關.注意此時的r表示距離z軸的距離)，所有變量都利用Hankel變換(潘爾年等，1986)，即

(9)

其中，ξ表示波數.為了簡便，仍用與球諧展開系數相同的符號表示變換后的系數.在此情形下，球函數趨近于柱函數，即

(10)

并且

(11)

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

(12)

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

(13)

對于Boussinesq問題，一個點負荷作用下的解為(在Hankel變換域)

(14)

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

(15)

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

1.2 Okubo和Endo方法

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

(16)

其中，c是由6個常數標量組成的待定系數矢量，且

[D(r)]=

(17)

(18)

(19)

式中

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

(20)

(21)

(22)

(23)

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

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

(24)

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

1.3 Guo方法

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

(25)

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

(26)

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

(27)

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

在此過程中，Guo等(2004)未介紹如何求得(27)式的解析解，而是直接給出了結果.因為該方程組實際上是變系數的，盡管將n/r移到了方程組的左邊.所以其是否存在解析解，并不是那么直觀地能夠看出來.實際上，因為(27)式是歐拉型，所以通過變量代換，

r=eξ,

(28)

可將該方程組化成常系數的微分方程組，從而有解析解.Guo等(2004)以一級漸近值為基礎，求解關于高一級小量的非齊次微分方程組，給出了更高一級的漸近值.但是對于負荷問題，更高一級的漸近值在計算負荷格林函數時是沒有必要的，因為利用(15)式的漸近值已經能夠很高效地獲得收斂的格林函數(參考下文的圖2).

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

2 解析解方法

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

2.1 解析解構造

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

(29)

和

(30)

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

(31)

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

此外，作一個變量代換

(32)

顯然，對于j層的下邊界有ξ=0.再以rTL，rTM和rQ為新變量，則(1)式變為

(33)

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

(34)

其中，

(35)

為了方便，我們用球諧階數n對變量作了規格化處理使各個變量量級相當.相應的矩陣A可從(33)式導出(系數作相應的調整，乘上與n有關的系數).此時，通解可以表示成(如在地球模型的第j層)

(36)

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

(37)

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

(38)

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

2.2 DVP傳播矩陣法

該方法的核心思想是避免導致溢出的大數產生，其技巧僅在于對普通的傳播矩陣法進行小的改動，即

1)將(36)式改寫為

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

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

(40)

其中，

(41)

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

2.3 漸近值推導

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

(42)

由(40)式進一步得到

(43)

這說明，當球諧階數很大時，U和T是相關的，這種相關性會導致傳統的龍格-庫塔數值積分法無法準確計算出非常高階的負荷勒夫數; 然而這種相關性可以使我們能夠直接求得表面的負荷勒夫數漸近值.對于地球的最外層(設為第p層)，我們可由(43)式直接得出

(44)

其中，方括號的下標p表示方括號內的矩陣由第p層的參數求得.這說明，只要我們知道最外層的系數矩陣A的特征向量，就可以得到當n趨于無窮大時的U，從而得到負荷勒夫數的漸近值.

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

圖1 負荷勒夫數隨球諧階數的變化

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

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

(45)

(45)式是高階負荷勒夫數的嚴格解析表達式，給出了負荷勒夫數與球諧階數n的關系.特別地，當n趨于無窮大時, 取極限就得到

(46)

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

ρ0pr0p=ρaa,

(47)

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

2.4 數值結果

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

此外，我們也將本文給出的漸近值與Guo等(2004)的結果進行了比較.結果說明，在n<20000左右時，Guo等(2004)的漸近值與考慮自引力的數值結果更加接近，這是因為他們推導漸近值時部分地考慮了自引力的影響；但是對于n>20000的部分，我們的結果不管是和考慮自引力還是和不考慮自引力的數值結果都更加接近，對于位移的勒夫數，這種情況更加明顯.

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

3 討論

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

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

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

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

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

4 結論

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

利用DVP傳播矩陣法對于n非常大時所具備的特點，我們不但可以看出變量U和T的關聯性，這也是導致傳統的龍格-庫塔數值積分法失效的本質原因，還可以直接根據地表邊界條件求出負荷勒夫數的漸近值.這說明DVP傳播矩陣法在數值計算方面有其獨特的優勢，可以廣泛應用于分層地球模型的變形研究中，如負荷形變、地震形變等靜態變形問題，以及地球簡正模或合成地震圖等動態變形問題中.

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