地球變形近似解析解的構建及應用

周江存, 潘爾年, 孫和平,3, 徐建橋, 陳曉東, 崔小明

1 中國科學院精密測量科學與技術創新研究院， 大地測量與地球動力學國家重點實驗室， 武漢 430077 2 臺灣陽明交通大學， 災害預防與水環境研究中心， 臺灣新竹 300 3 中國科學院大學， 北京 100049

0 引言

地球在內力和外力作用下都會產生變形，盡管人們無法感知有些力源產生的變形，如日、月等天體的引潮力導致的潮汐變形(Melchior，1983)，但是安置在地表或衛星中的高精度大地測量觀測儀器能夠清楚地識別出這些信號(Sun et al.，2019).因此基于高精度觀測的一些科學研究，如地球參考框架的建立與維持，需要扣除由于地球變形導致的干擾信號，如前述的固體潮信號，以及地表質量負荷、地震等(Altamimi et al.，2016).

描述地球變形的原理是基于無限小形變的假設，僅保留微小形變的一次項，而舍去高次項.通過數學手段建立的物理模型就包括三個基本方程：即(1)描述物質受力平衡的方程，對于靜態的問題可以不考慮慣性項——位移對時間的二次導數；(2)描述彈性或非彈性介質應力應變狀態的本構方程；(3)描述擾動位與擾動密度關系的泊松方程.將上述方程中的變量進行球諧展開后得到關于展開系數的一階常微分方程組，其包含八個變量的八個常微分方程(可參考Alterman et al.，1959；Takeuchi and Saito，1972；郭俊義，2001；徐建橋和孫和平，2003；孫文科，2012).其中，水平位移和剪切應力又有球型和環型之分，因此這八個變量中分別有兩個表示水平位移和剪切應力.

通常用勒夫數表示地球對力源的變形響應特征，它們分別對應著上述方程組中的位移和附加引力位的解.勒夫數是Love(1911)在研究地球的潮汐形變時定義的用來描述地球受引潮力作用而變形的兩個無量綱的量，用h和k表示，分別描述對應于某球諧展開階數的垂直位移與平衡潮潮高(引潮力位/地表重力)的比值和附加引力位與天體對地球引潮位的比值.后由Shida(志田)定義l描述水平位移與平衡潮水平分量的比值(Shida, 1912)，這三個數一般統稱為勒夫數(其中l也稱志田數).后來在研究地球受表面負載作用的變形也采用類似的三個無量綱的數，稱為“負荷勒夫數”(Farrell，1972)；Saito(1978)定義了描述地球受地表剪切力而變形的“剪切勒夫數”.特別地，潘爾年和丁中一(1986)揭示了上述三組勒夫數之間的關系.此外，研究地球受點源位錯(地震)作用而變形時也有類似的三個“位錯勒夫數”(Sun and Okubo，1993).Okubo(1993)也揭示了地表位錯勒夫數與震源處上述三組勒夫數之間的關系.Sun和Dong(2013)利用這種關系計算了地表位錯勒夫數及淺源地震的位錯勒夫數漸近值.

如果我們考慮一個SNREI(Spherically symmetric, Non-Rotating, Elastic and Isotropic，球對稱、非自轉、彈性各向同性)且包含自重的分層地球模型，那么地球變形的球型場和環型場是解耦的，所以由六個常微分方程共同描述了球型場的變形，由兩個常微分方程共同描述了環型場的變形.再基于地球所受力源的性質所確定的邊界條件，求解這個邊值問題，地球的受力變形就可以在理論上模擬出來.由于環型場的問題只包含兩個微分方程，因此求解相對比較簡單，主要的難點體現在包含六個微分方程的球型場的解算中.

由于常微分方程的系數與地球的參數，如密度、重力、拉梅常數，以及球坐標系的矢徑r有關，因此在計算上存在一些困難.當然，如果是一個簡單的地球模型，如均勻球，求解就非常簡單，這得益于Love(1911)從復雜的公式中發現了簡潔的嚴格解析解這一偉大的工作，他給出的結果到現在仍得到廣泛的應用(如Sun，2003，2004a，b；Takagi and Okubo，2017；Tang and Sun，2017).然而，對于一個分層的地球，沒有嚴格意義上的解析解.Gilbert和Backus (1968)通過假設每層中重力是r的線性函數，其他參數是常數，得到了一組用球貝塞爾函數表示的解析解.然而，對于一個接近真實地球的地球模型來說，在地核中重力近似是r的線性函數，因此關于重力的假設在地核中是合理的，但是在地幔中重力不是r的線性函數，該假設是不合理的，其結果必然帶來很大的誤差.因此，一直以來數值計算方法仍是求解地球變形時常用的方法：如龍格-庫塔數值積分法(如汪漢勝等，1996)、差分法(Poulsen，2009)、譜有限元方法(如Martinec，2000；廖彬彬等，2019).其中，四階龍格-庫塔方法是模擬SNREI地球形變最為廣泛使用的方法.

在用龍格-庫塔方法模擬計算中，第一個出現的困難是計算溢出的問題，這受限于計算機的數據表示及存儲能力.這個問題從模擬海潮負荷形變時得到關注，因為涉及的常微分方程組是與球諧階數相關的.不同于地球固體潮的研究只需要球諧展開的2階和3階的結果(Melchior，1983)，至多到4階(Xi，1989)，海潮負荷需要至少數千階的結果并結合漸近值以便采用Kummer變換(原理見下文)獲得收斂的格林函數(Farrell，1972).因此為解決這個問題提出了許多方法，如Longman的無量綱化方法(Longman，1963).該方法使得各個變量的數值相差不是太大(如幾個量級的差異)，計算更加穩定一些，但是這種方法只是讓能夠計算的最大球諧階數提高一些而已.汪漢勝等(1996)注意到每個變量都包含一個rn因子(其中，n是球諧階數)，提出了r冪因子法，即從每個變量中扣除rn因子.單純采用r冪因子法，也不能徹底解決溢出的問題，需要在積分過程中的每一步進行歸一化(如使所有變量的最大值為1).這實際上是每一步都從所有的變量中提取一個共同的常數，從而避免了溢出.鑒于所處理的微分方程組是線性的，所以可以通過地表的邊界條件加以約束.實際上所有積分過程中提取的常數的乘積都被最后利用邊界條件確定的未知常數吸收.因此，如果只計算地表的值，可以無需考慮這些提取的常數大小，然而如果要計算地球內部的形變，需要將這些提取的常數存儲并代回.

解決了溢出問題，對于大部分問題來說，數值模擬都能夠很好地開展.但是在模擬由地震引起的變形問題中，如果震源與場點(計算變形的點)距離比較近，那么所需要的球諧階數n就非常高.Sun和Okubo(1993)給出了一個經驗公式：所需計算的最大階數為10a/d，其中a是地球半徑，d是震源與場點的距離(r方向的投影).我們對此截斷公式也進行了深入的研究，考慮到每個變量的不同漸近性質給出了新的截斷公式(Zhou et al.，2019b).一般情況下，如果d=1 km，需要n=6萬階，如果要模擬地下爆炸對地表的變形影響，需要數十萬階.龍格-庫塔方法對于高階甚至超高階的計算是無能為力的.當n很大時，從地心(實際計算中從中心均勻小球表面)開始，不管初始值如何設定，積分幾步后幾組獨立解都變得線性相關，甚至一樣.這會導致最后在地表確定未知系數時的矩陣產生奇異.我們發現這個情況的出現是因為六個變量不再獨立(詳見下文)，這是由地球變形的物理問題所決定的，屬于內在的屬性.只要采用六個常微分方程，這種情況就不可避免.這種情況在有限元方法中就不會出現，因為有限元方法僅僅以位移和附加位為變量，從而不存在變量相關的問題(廖彬彬等，2019).

計算中最重要的一點是效率問題，在龍格-庫塔方法中，需要積分步長足夠小以便每一層的各個物理量可以合理地假設為常數.我們對位錯勒夫數的計算表明，積分步長不能超過100 m，否則球諧展開的高階結果就會出錯(Zhou et al.，2019a).這樣就需要將地球分成非常多的層，導致計算速度非常慢.在有限元中，使用插值函數同樣也需要節點之間的距離合適，對于球諧階數較高的情形，最終的剛度矩陣也非常大，也會影響計算效率.實際上，有限元方法的優勢主要體現在對三維地球模型的處理.

在地球變形的理論模擬中，如何使計算速度很快，又能夠避免出現計算的溢出問題，并且當球諧階數很高時也能輕松計算，這正是本文將要介紹的解析解方法結合DVP(Dual Variable and Position)傳播矩陣法所具備的特征.

1 地球變形的解析解

由于環型變形求解非常簡單，因此僅討論球型變形.求解球型變形的常微分方程組為(Takeuchi and Saito,1972；Pan et al.，2015)

(1)

其中，UL、UM、φ、TL、TM、Q分別是垂向位移、水平位移、附加引力位、垂向正應力、垂向剪應力和引力位梯度(參考第三個等式)，r為球坐標系的徑向分量，λ，μ是拉梅常數，ρ是密度，g是重力，n是球諧階數.值得注意的是式(1)為齊次微分方程組，對于一個有內部力源的地球變形問題，微分方程組本來應該是非齊次的，但是如果我們假設力源是點源，那么可以將非齊次項轉變成力源處的邊界條件(參考Takeuchi and Saito，1972).因此只需要求解如式(1)的齊次方程組.

歸根結底，求解微分方程組(1)的困難在于系數矩陣是變量，從而解不能簡潔地表示出來.如果基于一些合理的假設，能夠使得系數矩陣變成常數矩陣，那問題就非常簡單了.解析解方法正是基于這樣的思路，首先將地球分成若干層，然后對每一層作一些近似(Pan et al.，2015)：在液態地核的分層中假設密度為常數，重力是r的線性函數；在地幔和地殼的分層中，重力是常數而密度是r的反函數.此外，假設拉梅常數在每一層中為定值，并且定義一個新變量

ξ=ln(r/rj-1),rj-1≤r≤rj,

(2)

式中，rj和rj-1分別是第j層的上下邊界.將新變量和上述的假設代入式(1)，我們就可以得到如下用矢量形式表示的微分方程.

dY/dξ=AY,

(3)

其中Y=(UL,UM,φ,rTL,rTM,rQ)T，最重要的是此時A為6×6的常數矩陣，因此方程(3)存在解析解，該解由系數矩陣A的特征值和特征向量表征，即

(4)

其中，所有的子矩陣都是3×3的，c1和c2是待定系數矢量，各包含3個未知數，

(5)

該方法的關鍵是對密度和重力在地球內部的分布作了假設，從而系數矩陣變成了常數矩陣.那么這些假設是否合理，是否會對結果產生很大的影響是我們首先需要回答的問題.通過與PREM(Preliminary Reference Earth Model，初步參考地球模型)(Dziewonski and Anderson, 1980)比較發現，由上述假設生成的模型與PREM模型的差異非常小，因此這種假設是合理的(Pan et al.，2015).實際上，上述的假設也是在充分考慮PREM模型的參數分布特性的基礎上做出的.最重要的一點是：我們可以在顧及計算效率的情況下，減低層厚使得二者之間的差異更小.實際上只要層厚足夠小，可以用任意的函數擬合模型參數且能夠較好地控制誤差.通過計算發現，模型的厚度可以設成數十公里(Pan et al.，2015；Chen et al.，2018；Zhou et al.，2019a).由此可見，這里的假設要比Gilbert和Backus(1968)合理得多.需要特別指出的是，如果層厚足夠小，那么得到的結果就是真實的解，而非近似解.

盡管得到了理論上的解析解，但是對于PREM這樣的地球模型，特征值和特征向量是不太可能有簡單的解析形式或由普通函數直接表示出來.因此我們需要采用數值計算的方法計算特征值和特征向量，可以理解為是一種關于拉梅常數和球諧階數的特殊函數.如同均勻球的解用球貝賽爾函數表示，雖然球貝賽爾函數是解析的，但仍然需要采用遞推公式及適當的技巧以便能夠準確地將它數值地計算出來.

然而，直接利用式(4)并采用普通的傳播矩陣方法，在球諧階數很高的時候，也會導致數值計算的溢出(Pan et al.，2015).因此，我們介紹一種新的傳播矩陣方法，稱為DVP方法(Pan，2019；Zhou et al.，2019a)，使得傳播穩定而高效.

2 DVP傳播矩陣法

普通傳播矩陣法溢出的原因是：當球諧階數很高時(如n=10000)，矩陣的6個特征值接近于n+1、n、n-1、-n、-(n+1)、-(n+2)，其中，如式(5)所示，對于正的特征值，指數eτ ξ就會是一個更大的數，超出了計算機的存儲能力，從而導致向上溢出.雖然可以把層厚取得足夠小(ξ很小)，從而某一層內不會溢出，但是傳播數層之后，溢出還是會出現.DVP傳播矩陣法就是避免在傳播過程中直接使用含有大的正實部的特征值參與計算，原理如下：

首先將式(4)改成

(6)

其中ξj即ln(rj/rj-1)為j層的上邊界的值.實際上這種改變是從c1中提取了某一個常數，我們用c3表示新的待定常數.顯然有如下的關系：

〈-s1ξj〉c3=c1.

(7)

在某種程度上，這種改變類似于龍格-庫塔方法每一步積分過程中實施的歸一化.為了方便，定義

(8)

式中，我們用n對變量進行了規格化使得每個變量具有相同的量級，這只是使矩陣A稍作改變而已，從而特征向量也作相應調整.但這并不影響問題的本質(特征值不變)以及我們對方法的描述.這樣一來，在第j層的上邊界(ξ=ξj)和下邊界(ξ=0)，分別有

(9a)

和

(9b)

其中E為3×3的單位陣，上標u和l表示上邊界和下邊界.通過將上面的(9a)和(9b)兩式寫成4個由子矩陣表示的等式，調整變量順序，并消去系數后再表示成矩陣形式，有

(10)

其中

(11)

需要注意式(10)中上下邊界處變量的位置，普通的傳播矩陣法的形式為(Pan et al.，2015)

(12)

式中，傳播矩陣P′的具體形式可以很容易從式(4)或(9)導出，剛度矩陣法的形式為(Chen et al.，2018)

(13)

式中，傳播矩陣P″也可以很容易從式(9)導出.可以看出式(10)，(12)和(13)三種方法中上下邊界處變量的位置是明顯不同的，而正是這種簡單的位置調換產生了不一樣的計算效果.

從式(11)可以看出，雖然s1包含具有正實部的特征值，但是前面有一個負號，因此避免了大數的產生，從而防止了溢出.特別地，當n非常大時，可以得到

(14)

因此

(15)

這說明，當球諧階數很大時，U和T是相關的，這也是為什么龍格-庫塔方法無法計算高階解的原因.然而這種相關性可以使我們能夠直接求得位錯勒夫數在震源處的值，從而求出當n很大時候的位錯勒夫數漸近值，并可實施Kummer變換獲得收斂的位錯格林函數(Zhou et al.，2021).此外，我們還發現，

(16)

其中a是地球半徑.也就是說，對于地表負荷這類在地表有非零值邊界條件的問題來說，地表的位移和附加引力位可直接由最上層的特征向量求出，而不管地球內部的結構如何，這也解釋了為什么高階負荷勒夫數只與地殼最外層結構有關.這實際上也是求當n很大時負荷勒夫數漸近值的新方法(見周江存和潘爾年，2021).

有了相鄰兩層的傳播矩陣，就可以將解進行傳播，j層和j+1層之間的傳播矩陣為

(17)

其中

(18)

重復迭代公式(18)，就可以得到初始層和目標層之間的傳播矩陣，從而得到兩個界面上的解之間的關系，由邊界條件即可求得最終的解(Zhou et al.，2019a).

3 爆炸源格林函數

不同的勒夫數是構建不同格林函數以及計算相應力源引起的地面和內部變形導致的位移場、重力場和應變應力場變化的基礎(Pan，2019).需要指出的是，我們定義的位錯勒夫數與以往如Sun和Okubo(1993)略有不同，采用了不同的規格化(詳情可參考Zhou et al.，2019a)，新定義有助于我們直觀地發現它們的漸近性質，即隨著n的增加，它們具有怎樣的變化特征，從而為能夠擬合出漸近值給出直觀的感知，為能夠實施Kummer變換從而加速格林函數的收斂創造了條件(Zhou et al.，2019b).此外，我們最近也提出了另一種方法獲得位錯勒夫數漸近值，即根據傳播矩陣的特點直接求解位錯勒夫數的漸近值(Zhou et al.，2021).漸近值的成功計算為我們能夠模擬淺源爆炸引起的地表形變奠定了基礎.對于像淺斷層或淺源爆炸的情形，我們再無需疊加像上文提及的數十萬那么高階的位錯勒夫數，只需實施Kummer變換疊加少量低階的位錯勒夫數即可.

格林函數涉及球函數的求和，其系數就是勒夫數或其組合.對于球對稱的地球模型，任意的斷層都可以表示成四種獨立的斷層模式的組合(Sun and Okubo，1993；Zhou et al.，2019a)，它們是垂直走滑、垂直傾滑、水平拉張和垂直拉張，分別用兩個數字組合表示，即分別為12，32，22，和33.爆炸源的矩張量可以表示為(Kumagai et al.，2014)

M0=3KsΔVc=(3λs+2μs)ΔVc,

(19)

其中，K表示體膨脹系數，下標s表示震源處的值，ΔVc是由任一個正交的拉張斷裂(11，22或33，由于對稱性，11可通過22的旋轉獲得)產生的無應力體積變化(Aki and Richards，2002).因此爆炸源的格林函數就可用位錯格林函數表示為(也可參考孫文科，2012)

(20)

其中，(r,θ,φ)構成球坐標系，

(21)

(22)

其中，a1、b1、c1和d1均為常數，可由高階位錯勒夫數擬合而得(Zhou et al.，2019b)，或根據分層地球模型的參數直接計算出來(Zhou et al.，2021).(22)式中的兩個等式右端的第一個求和部分比(21)式要收斂得快得多，因此可截斷到nMax，而第二個求和具有嚴格的解析表達式(參考Zhou et al.，2020；Tang and Sun，2021).

4 數值結果與討論

我們以地下100 m、200 m和500 m深處的爆炸為例，模擬爆炸產生的永久變形.首先計算對應不同震源深度的地表位錯勒夫數.采用解析解的方法避免了傳統的龍格-庫塔方法中的變量相關問題，可以很容易地計算超高階的位錯勒夫數.圖1和圖2分別顯示了震源深度為100 m時的球諧階數n從0到107階的位錯勒夫數h和nl的計算結果.由于我們考慮的問題只涉及兩個獨立的拉張型位錯模式，所以圖中只給出這兩個位錯模式的結果.

圖1 地表位錯勒夫數hn(震源深度100 m. 左側兩圖是0～100階的結果，中間兩圖是100～104階的結果，右側兩圖是104～107階的結果)Fig.1 Dislocation Love number hn on the Earth′s surface for source depth at 100 m (the left two plates are for degrees 0～100, the middle two plates are for degrees 100～104 and the right two plates are for degrees 104～107)

圖2 地表位錯勒夫數ln (震源深度100 m. 左側兩圖是0～100階的結果，中間兩圖是100～104階的結果，右側兩圖是104～107階的結果)Fig.2 Dislocation Love number ln on the Earth′s surface for source depth at 100 m (the left two plates are for degrees 0～100, the middle two plates are for degrees 100～104 and the right two plates are for degrees 104～107)

從圖中也可以很容易地看出位錯勒夫數的漸近特征，兩種位錯勒夫數在高階表現出與n呈線性的關系，這在垂直拉張的結果中更加明顯，從低階的結果中就開始顯現，其變化看上去就是一條直線.通過這種漸近特征我們可以擬合出位錯勒夫數漸近值的線性表達式，其表達式的系數列于表1中，這些系數對我們能夠使用Kummer變換加速位錯格林函數的收斂至關重要.從表中可看出，這些系數的線性項具有明顯的特點，h和nl的線性系數大小相等，符號相反；垂直拉張位錯勒夫數的線性系數正好是水平拉張結果的-2倍.

表1 位錯勒夫數漸近值Table 1 Asymptotic dislocation Love numbers

圖3和圖4分別給出了不同深度處的無應力單位體積(即假定M0=3K)爆炸源垂直位移和水平位移的格林函數.不同于地震斷層需要分成若干的子斷層，每個子斷層認為是點源，再疊加所有子斷層對位移的貢獻，爆炸源本身就是一個點源事件，所以從其格林函數中就可看出其對位移場的影響大小.并且位移場是以爆炸源與地心連線為軸旋轉對稱的，即位移場只與球面角距離有關，而與方位無關.

圖3 地下不同深度處無應力單位體積爆炸源引起的地表垂直位移(單位：m)Fig.3 Vertical displacements (unit: m) on the Earth′s surface due to underground explosions of unit stress-free volume at different depths

圖4 地下不同深度處無應力單位體積爆炸源引起的地表水平位移(單位：m)Fig.4 Horizontal displacements (unit: m) on the Earth′s surface due to underground explosions of unit stress-free volume at different depths

為了顯示爆炸源引起的全球變化特征，我們將角距離分成不同的區域以顯示位移場變化的詳細特征.為了顯示方便，在每個區域都從格林函數中提取出一個因子，標注在圖中的右上角.盡管對于爆炸源來說，遠場的位移是非常小的，這里將遠場的結果也顯示出來是為了(1)保持格林函數的完整性；(2)顯示格林函數的收斂性以說明位錯勒夫數漸近值在格林函數計算中的重要性.從圖中可以看出，隨著深度的增加，地表的變形會逐漸減小，顯著的位移場只發生在離爆心大約100～200 m的范圍之內，垂直位移的幅度在6 km之外降低5個數量級，水平位移的幅度在6 km之外降低4個數量級，在50 km之外降低5個數量級.因此，如果爆炸當量比較小，那么顯著的永久形變只會發生在離爆炸中心不遠的地方.當然，這個距離會隨著爆炸當量的增大而增大，同時也隨著爆炸源深度的減小而增大.

5 結論

我們介紹了一種計算地球受力變形勒夫數的解析解方法.該方法通過合理假設地球內部結構參數，使得變系數的常微分方程組變成常系數的常微分方程組，從而其解可以用系數矩陣的特征值和特征向量解析地表達出來.

通過引入一種特殊的DVP傳播矩陣方法，解決了計算的溢出和球諧展開高階項的變量相關問題，從而可以穩定和高效地計算出超高階的位錯勒夫數，并且利用該方法可直接求得震源處的位錯勒夫數漸近值，從而確定任意位置的位錯勒夫數漸近值.特別地，從計算得到的勒夫數隨階數的變化中能夠很直觀地發現它的漸近特征，從而擬合出漸近值.采用Kummer變化加速了位錯格林函數的收斂，由此計算了淺源爆炸的格林函數.

我們模擬分析了淺源爆炸產生的地球表面的永久形變特征.結果表明，采用解析解方法結合DVP傳播矩陣方法可以解決位錯勒夫數和格林函數的計算問題.最后需要說明的是，本文的數值結果是基于靜態的變形理論，未考慮爆炸產生的波及其引起的破壞性變形.

致謝特別感謝我們前期系列工作的合作者美國俄亥俄州立大學Michael Bevis教授.