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

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院， 大地測量與地球動力學(xué)國家重點實驗室， 武漢 430077 2 臺灣陽明交通大學(xué)， 災(zāi)害預(yù)防與水環(huán)境研究中心， 臺灣新竹 300 3 中國科學(xué)院大學(xué)， 北京 100049

0 引言

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

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

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

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

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

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

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

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

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

1 地球變形的解析解

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

(1)

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

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

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

(2)

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

dY/dξ=AY,

(3)

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

(4)

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

(5)

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

盡管得到了理論上的解析解，但是對于PREM這樣的地球模型，特征值和特征向量是不太可能有簡單的解析形式或由普通函數(shù)直接表示出來.因此我們需要采用數(shù)值計算的方法計算特征值和特征向量，可以理解為是一種關(guān)于拉梅常數(shù)和球諧階數(shù)的特殊函數(shù).如同均勻球的解用球貝賽爾函數(shù)表示，雖然球貝賽爾函數(shù)是解析的，但仍然需要采用遞推公式及適當(dāng)?shù)募记梢员隳軌驕?zhǔn)確地將它數(shù)值地計算出來.

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

2 DVP傳播矩陣法

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

首先將式(4)改成

(6)

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

〈-s1ξj〉c3=c1.

(7)

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

(8)

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

(9a)

和

(9b)

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

(10)

其中

(11)

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

(12)

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

(13)

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

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

(14)

因此

(15)

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

(16)

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

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

(17)

其中

(18)

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

3 爆炸源格林函數(shù)

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

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

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

(19)

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

(20)

其中，(r,θ,φ)構(gòu)成球坐標(biāo)系，

(21)

(22)

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

4 數(shù)值結(jié)果與討論

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

圖1 地表位錯勒夫數(shù)hn(震源深度100 m. 左側(cè)兩圖是0～100階的結(jié)果，中間兩圖是100～104階的結(jié)果，右側(cè)兩圖是104～107階的結(jié)果)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 地表位錯勒夫數(shù)ln (震源深度100 m. 左側(cè)兩圖是0～100階的結(jié)果，中間兩圖是100～104階的結(jié)果，右側(cè)兩圖是104～107階的結(jié)果)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)

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

表1 位錯勒夫數(shù)漸近值Table 1 Asymptotic dislocation Love numbers

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

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

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

5 結(jié)論

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

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

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

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