自重可壓縮層狀球體地球模型下地震變形算法的優(yōu)化研究

劉泰， 唐河， 佘雅文， 付廣裕， 孟國杰

1 中國地震局地震預(yù)測研究所， 北京 100036 2 中國科學(xué)院大學(xué)計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049 3 中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京 100083

0 引言

隨著大地測量技術(shù)的發(fā)展，GNSS(Global Navigation Satellite System)和GRACE(Gravity Recovery And Climate Experiment)等手段觀測到的同震及震后變形數(shù)據(jù)日益豐富，為研究地震破裂過程和震后變形機(jī)制提供了數(shù)據(jù)支持(刁法啟等, 2012; 劉泰等, 2017; 梁明等, 2018; Wang et al., 2022).其中，近場地震變形研究能夠加深人們對地震震后余滑時(shí)空分布和震源區(qū)地幔黏滯性構(gòu)造橫向特征的認(rèn)識(Diao et al., 2014; Hu et al., 2016)，而引入遠(yuǎn)場地震變形數(shù)據(jù)能夠約束地震震級和深部地幔黏滯性結(jié)構(gòu)(陳飛等, 2020; Wang et al., 2014).地震位錯(cuò)理論是研究地震變形的重要理論工具，其基于的地球模型直接影響地震變形的計(jì)算精度，如地球模型分層結(jié)構(gòu)和不連續(xù)性影響近場地震變形計(jì)算(Dong et al., 2014, 2021)，地球曲率會(huì)影響遠(yuǎn)場地震變形計(jì)算(Dong et al., 2016)，地球自重效應(yīng)會(huì)影響震后變形計(jì)算(Pollitz, 1997)，長時(shí)間震后變形與地球的可壓縮相關(guān)(Tanaka et al., 2006).可見，基于接近真實(shí)的地球模型，完善地震變形計(jì)算方法具有重要的理論意義和應(yīng)用價(jià)值.

數(shù)十年來，學(xué)者們基于不同幾何特征的地球模型發(fā)展了多種地震變形計(jì)算方法.例如，Okada(1985, 1992)通過總結(jié)整理前人的工作，給出了一套完整、簡潔、實(shí)用的同震形變計(jì)算公式，成為均勻半無限空間地震位錯(cuò)理論的經(jīng)典表達(dá)式.Wang等(2003, 2006)發(fā)展了半無限空間模型的地震位錯(cuò)理論，考慮了地球的層狀結(jié)構(gòu)，提高了近場地震變形的計(jì)算精度.Sun(1992)和Pollitz(1996)發(fā)展了層狀球?qū)ΨQ地球模型的地震位錯(cuò)理論，考慮了地球曲率，提高了遠(yuǎn)場地震變形的計(jì)算精度.

相比于外部力源，如負(fù)荷、潮汐和剪切，產(chǎn)生的地球變形問題，地震變形問題略微復(fù)雜，前者需要求解齊次微分方程組，后者需要求解非齊次微分方程組.地震變形的微分方程組含有力源項(xiàng)(Takeuchi and Saito, 1972)，為非齊次問題.Okubo(1993)提出了地球內(nèi)部點(diǎn)源位錯(cuò)在地表產(chǎn)生的變形與外部力源在震源處產(chǎn)生的形變存在互換關(guān)系(互換定理)，從而簡化了地表地震變形的計(jì)算過程.當(dāng)研究地震引起的地球內(nèi)部變形時(shí)，微分方程組的解需滿足震源處不連續(xù)和地表自由邊界條件.為此，Takeuchi和Saito(1972)和Sun(1992)提出了兩次數(shù)值積分、組合通解和特解的求解方法：從地心積分至地表得到通解，由震源處積分至地表得到特解，再利用地表邊界條件確定通解的組合系數(shù)，進(jìn)而得到微分方程組解矢量.該求解方法逐步成為了求解地震變形的常用方法之一.

當(dāng)觀測點(diǎn)靠近震源時(shí)，地球內(nèi)部變形的格林函數(shù)需要通過高階解矢量來得到.但對于高階微分方程組求解，上述方法在數(shù)值積分的過程中會(huì)出現(xiàn)數(shù)值溢出和精度丟失現(xiàn)象.數(shù)值溢出問題是由向上積分過程中積分區(qū)間過長引起，出現(xiàn)在負(fù)荷變形和地震變形的計(jì)算中.最初提出的解決方法是調(diào)整積分起始半徑，使其隨階次的增加而增大(Sun and Okubo, 1993; Martens et al., 2019).更為有效的方法是從數(shù)值放大因子(rn，這里n表示階數(shù)，r表示半徑)出發(fā)，在積分過程中對該因子進(jìn)行分離(汪漢勝等, 1996; Liu et al., 2020).精度丟失現(xiàn)象是由于震源上方通解和特解組合引起，在微分方程組的高階求解中尤為顯著.一個(gè)直接的解決方法是使用第三方擴(kuò)展函數(shù)庫來進(jìn)行高存儲精度的計(jì)算，該方法已在黏彈性變形計(jì)算中得到應(yīng)用(Melini et al., 2008).但使用該方法編寫的計(jì)算程序依賴于第三方數(shù)據(jù)庫，給用戶帶來了不便.更為有效的方法是調(diào)整數(shù)值積分策略、改變積分路徑，對地震變形微分方程的齊次部分進(jìn)行兩次相向積分，即由地心積分至震源和由地表積分至震源，再利用震源間斷條件求解震源上方和下方通解組合系數(shù)(Pan, 1997; Wang, 1999; Fukahata and Matsu′ura, 2005; Zhou et al., 2021).

為此，本文使用兩次相向積分的方法來求解地震變形的微分方程組，以解決高階微分方程組求解的數(shù)值不穩(wěn)定問題，從而得到微分方程組解矢量.首先介紹了微分方程組的兩種求解方法，接著與均質(zhì)地球模型解析解進(jìn)行比較來驗(yàn)證方法正確性，最后研究層狀地球模型下震源附近地震變形.

1 地震變形理論

1.1 地震變形的基本方程

假設(shè)地球內(nèi)部r=rs處存在點(diǎn)源位錯(cuò)，如圖1a所示，其激發(fā)的位移場、應(yīng)力場和引力位變化同時(shí)滿足平衡方程、應(yīng)力-應(yīng)變關(guān)系和泊松方程(Takeuchi and Saito, 1972).利用球諧函數(shù)可以將以上物理量展開，整理可得兩組微分方程組，分別對應(yīng)球型變形和環(huán)型變形.平衡方程最初的非齊次項(xiàng)是體力，而位錯(cuò)源的等效體力包含了δ函數(shù)和δ函數(shù)的導(dǎo)數(shù)項(xiàng).由于點(diǎn)位錯(cuò)源也可以用雙力偶表達(dá)，Takeuchi和Saito(1972)利用震源處位移和應(yīng)力間斷特征推導(dǎo)了震源函數(shù)表達(dá)式，實(shí)際計(jì)算中只要把上述體力用雙力偶代替即可.

dYs/dr=AsYs+δ(r-rs)ss,

(1)

dYt/dr=AtYt+δ(r-rs)st,

(2)

如圖1b所示，斷層上下盤相對運(yùn)動(dòng)矢量為Uv，坐標(biāo)軸e1和e2分別表示緯度φ=0°和φ=90°的方向，坐標(biāo)軸e3對應(yīng)極軸方向，位錯(cuò)滑動(dòng)矢量v和斷層面法線矢量n可以表示為

(3)

式中v1,v2,v3和n1,n2,n3是v和n在球坐標(biāo)系下的三個(gè)分量，即vinj(i,j=1,2,3)，為了方便后面位錯(cuò)源的表示，用角標(biāo)ij代替vinj.Sun(1992)定義了四種獨(dú)立震源，即垂直斷層水平走滑震源(ij=12)、垂直斷層上下傾滑震源(ij=32)、垂直斷層水平引張震源(ij=22)和水平斷層垂直引張震源(ij=33)，如圖1c所示.通過以上四種獨(dú)立震源的線性組合可以得到任意類型的位錯(cuò)源引起的變形.

圖1 球坐標(biāo)系下位錯(cuò)模型 (a,b)； (rs,0,0)為震源所在位置，(r,θ,φ)為球內(nèi)任意一點(diǎn)，斷層面的傾角為δ，滑動(dòng)角為λ，上下盤斷層面的相對錯(cuò)動(dòng)量為U，n和v分別為斷層面法線矢量和位錯(cuò)滑動(dòng)矢量； (c) 表示四種獨(dú)立震源

1.2 數(shù)值積分求解微分方程組

Sun(1992)提出了球體地球模型下點(diǎn)源位錯(cuò)引起的地表形變的計(jì)算方法，通過兩次積分求解非齊次微分方程組，如圖2a 所示，本文將該方法稱為兩次同向積分法.以球型變形為例，此時(shí)該非齊次方程組的解矢量分為兩部分：

圖2 兩種數(shù)值積分路徑示意圖

Y=X·β+Z,

(4)

式中，X為對應(yīng)齊次方程組的通解矩陣(6×3)，Z為非齊次方程組的特解矩陣(6×1)，β為三組通解的組合系數(shù).其中特解部分由地心附近半徑很小的球面(r=r0)積分至地表(r=a)，Love(1911)最早給出了通解部分的積分初值，特解部分由震源函數(shù)從震源處積分至地表得到，以滿足震源處間斷條件.當(dāng)變形深度在震源上方時(shí)，解矢量由通解和特解組成，而震源下方的解矢量完全由通解組成：

(5)

滿足地表邊界條件，即地球表面應(yīng)力分量為零和引力位連續(xù)，以及引力位導(dǎo)數(shù)穿過單層的不連續(xù)性(Sun and Sj?berg, 1999; Sun, 1992)，有y2(a)=y4(a)=y6(a)=0，可以得到三組通解的組合系數(shù)βi(i=1,2,3).然后帶入公式(5)，得到任意深度處的解矢量(yj(r),j=1,2,…,6).

由于通解部分積分區(qū)間較長，在求解高階微分方程組時(shí)，為了避免積分過程中出現(xiàn)數(shù)值溢出，積分起始半徑需隨著階次的增加而增加(Sun, 1992).Liu等(2020)將該地表變形理論推廣至地球內(nèi)部，利用正規(guī)化方法避免地球內(nèi)部變形計(jì)算過程中數(shù)值溢出的問題.但是當(dāng)觀測點(diǎn)位于震源上方時(shí)，解矢量由通解和特解組合得到，階次足夠高時(shí)會(huì)出現(xiàn)精度丟失現(xiàn)象，本研究通過調(diào)整積分路徑解決該問題.具體操作是對微分方程的齊次部分進(jìn)行兩次積分，由地心積分至震源和由地表積分至震源(圖2b)，本文將該方法稱為兩次相向積分法.此時(shí)解矢量完全由通解組成：

(6)

(7)

(8)

滿足震源處間斷條件后(Y2(rs)-Y1(rs)=s)，同時(shí)得到震源下方和上方通解的組合系數(shù)αi和βi(i=1,2,3).

(9)

可以發(fā)現(xiàn)本研究兩次相向積分法和傳統(tǒng)的兩次同向積分法的差別在于處理震源間斷和地表邊界條件的先后順序，優(yōu)先處理地表自由邊界可以避免微分方程組非齊次部分的求解，微分方程組解矢量完全由通解組成，此時(shí)震源上方和震源下方具有不同的組合系數(shù).另外，兩者方法在震源下方的積分過程完全一致，孫文科(2012)中詳細(xì)介紹了內(nèi)外核邊界和核幔邊界的處理過程，這里不在贅述.

2 地震變形計(jì)算結(jié)果

基于上述基本解法，給定不同的源函數(shù)(ss,st)求解通解組合系數(shù)，可得到對應(yīng)震源的解矢量(Ys,Yt).本研究將基于Sun(1992)定義的四種獨(dú)立震源進(jìn)行求解，地球分層介質(zhì)參數(shù)選自PREM地球模型(Dziewonski and Anderson 1981)，如圖3所示，重力值由密度對地球半徑積分得到.

圖3 地球模型參數(shù)隨深度的變化

本研究將從兩個(gè)方面展示微分方程組的求解結(jié)果：(i) 忽略自重的均質(zhì)球體地球模型下，與解析解(Dong et al., 2021)比較，驗(yàn)證本研究求解方法的正確性；(ii) 層狀球體地球模型下，與兩次同向積分法(Sun, 1992; Liu et al., 2020)比較，展示本研究求解方法的優(yōu)勢.

2.1 均質(zhì)地球模型下解矢量

對于均質(zhì)地球模型和層狀地球模型的積分過程完全一致，只是每個(gè)子層上的介質(zhì)參數(shù)不同，為此，本研究只需驗(yàn)證均質(zhì)地球模型下求解的正確性，即可保證層狀地球模型下的正確性，均質(zhì)地球模型的介質(zhì)參數(shù)設(shè)置為PREM模型的頂層參數(shù)值(λ0=34.2 GPa,μ0=26.6 GPa,ρ0=2.6 cm3·g-1).忽略自重的均質(zhì)地球模型，微分方程組存在解析解(Dong et al., 2021)，為了與該解析解比較，本研究將均質(zhì)地球模型的密度減少為原來的1%，以達(dá)到忽略自重的效果.

圖4展示了均質(zhì)地球模型下四種獨(dú)立震源的解矢量(ds=20 km;dh=10 km，ds為震源深度，dh為變形觀測深度).由于忽略自重的均質(zhì)地球模型下解析解沒有y5和y6，并且環(huán)型變形不受自重的影響，這里僅僅展示了變形解y1,y2,y3,y4.可以發(fā)現(xiàn)當(dāng)密度減少為原來的1%時(shí)，數(shù)值積分解與解析解基本一致，驗(yàn)證了本研究求解方法的正確性.這里與Liu等(2020)的無自重處理方式不同，后者只是簡單的將模型重力值減少為1%，沒有改變密度值，可見通過減小模型的密度值以達(dá)到忽略自重的效果更為合理.另外，水平走滑震源受自重的影響較大，垂直傾滑受自重的影響較小，而自重主要影響垂直引張震源和水平引張震源的解y4，其中自重對垂直引張震源的影響超過100階.整體而言，自重主要影響微分方程組的低階求解，隨著階次增加，數(shù)值積分解與解析解趨于一致.

圖4 均質(zhì)地球模型下四種獨(dú)立震源的解矢量(ds=20 km;dh=10 km)

2.2 層狀地球模型下解矢量

這里主要比較了1.2部分中兩種積分路徑下的解矢量.介質(zhì)參數(shù)來源于真實(shí)分層地球模型(圖3)，由于地核主要影響10階以內(nèi)的微分方程組求解(Liu et al., 2020).在保證精度前提下節(jié)省計(jì)算時(shí)間，本研究將1000階以上的積分起始半徑設(shè)置為核幔邊界，由于在每層積分過程中采用正規(guī)化處理，高階微分方程組的求解不會(huì)出現(xiàn)數(shù)值溢出現(xiàn)象.

圖5展示了層狀地球模型下水平走滑震源的解y1(ds=20 km;dh=0 km,5 km,10 km,15 km).圖5a1中虛線表示該震源深度ds和變形深度dh下解的截?cái)嚯A次(Nmax=10a/|ds-dh|)，可以發(fā)現(xiàn)截?cái)嚯A次以內(nèi)，兩種積分路徑下的求解結(jié)果基本一致，超過截?cái)嚯A次后，傳統(tǒng)的求解方法出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象.這是由于隨著計(jì)算階次的增加，傳統(tǒng)求解方法中通解(圖c1中實(shí)線)和特解(圖c1中虛線)的數(shù)值量級逐漸增加，但兩者符號相反，絕對值逐漸接近.在計(jì)算機(jī)的正常存儲精度下，當(dāng)通解和特解的有效數(shù)字完全一樣時(shí)，組合得到的解矢量精度丟失，呈現(xiàn)數(shù)值震蕩現(xiàn)象(圖b1中紅色圓點(diǎn))，圖(b1)和(c1)的數(shù)值計(jì)算結(jié)果展示在表1中.另外可以發(fā)現(xiàn)超過截?cái)嚯A次后，該現(xiàn)象愈加明顯.這種數(shù)值不穩(wěn)定現(xiàn)象是組合通解和特解過程中精度丟失引起的(公式(5)).前面也已經(jīng)介紹通過第三方擴(kuò)展函數(shù)庫可以進(jìn)行高存儲精度的計(jì)算，但該方法增加了通解和特解的有效數(shù)字，能夠一定程度上提高計(jì)算階次，但無法從根本上解決該問題.但本研究通過調(diào)整積分路徑，震源上方和下方的解矢量均由通解組成(公式(6)和(7))，從根本上解決了該精度丟失的問題，得到的微分方程組的解矢量比較穩(wěn)定，如圖5b1中黑色曲線所示.

圖5 層狀地球模型下水平走滑震源的解y1(ds=20 km)

表1 層狀地球模型下水平走滑震源的解y1(ds=20 km;dh=0 km)

3 震源附近地震變形研究

地震引起的地表變形研究中，漸近解結(jié)果顯示層狀地球模型和均質(zhì)地球模型下的高階解矢量基本一致，說明了淺源地震引起的地表變形受層狀結(jié)構(gòu)的影響很小.因此，可以利用均質(zhì)地球模型下地表破裂源引起的近場地表變形對層狀地球模型進(jìn)行替代(Sun and Dong, 2013).本研究將進(jìn)一步探討層狀結(jié)構(gòu)對地球內(nèi)部地震變形的影響，計(jì)算了均質(zhì)地球模型和層狀地球模型下微分方程組解的比值：

(10)

其中，yinhomo.為層狀地球模型下的變形解，yhomo.為均質(zhì)地球模型下的變形解，均質(zhì)球參數(shù)來源于分層模型頂層(λ0=34.2 GPa,μ0=26.6 GPa,ρ0=2.6 cm3·g-1)，可以發(fā)現(xiàn)隨著階次的增加，兩者的比值趨于常數(shù)，如圖6中紅線所示.另外，當(dāng)均質(zhì)球的介質(zhì)參數(shù)設(shè)置為震源附近參數(shù)(λ0=45.9 GPa,μ0=44.1 GPa,ρ0=2.9 cm3·g-1)時(shí)，兩種地球模型下的高階變形解比值趨于1，如圖6中藍(lán)線所示，此時(shí)近場地震變形受層狀結(jié)構(gòu)的影響較小.這里只展示了走滑震源的計(jì)算結(jié)果，該特征對其余三個(gè)獨(dú)立震源也適用，為節(jié)省篇幅，不逐一展示.

圖6 均質(zhì)地球模型和層狀地球模型下走滑震源變形解的比值，縱坐標(biāo)為比值的絕對值(ds=20 km;dh=19 km)

該部分將主要關(guān)注震源附近的地震變形計(jì)算，當(dāng)變形計(jì)算深度接近震源深度時(shí)，需要對更高階次的解矢量求和，甚至需要計(jì)算無窮階次的解矢量，數(shù)值解無法實(shí)現(xiàn)該計(jì)算要求.Liu等(2020)計(jì)算了20 km震源在20 km處引起的變形，當(dāng)對105的解矢量進(jìn)行求和時(shí)，0.03°以內(nèi)的格林函數(shù)不收斂.Takagi和Okubo(2017)推導(dǎo)了均質(zhì)地球模型下地球內(nèi)部變形的漸近解，解決了均質(zhì)地球模型下震源附近地震變形不收斂的問題.本研究將進(jìn)一步解決層狀地球模型下震源附近地震變形不收斂的問題.

具體思路是對兩種地球模型下解矢量的差異進(jìn)行求和得到層狀結(jié)構(gòu)的影響，再與均質(zhì)地球模型的漸近解疊加，從而得到層狀地球模型下格林函數(shù).例如，計(jì)算層狀地球模型下垂直位移格林函數(shù)時(shí)，可以做一個(gè)變換，即加上和減去一個(gè)均質(zhì)地球模型下的變形解，而等式關(guān)系保持不變：

(12)

圖7展示了層狀結(jié)構(gòu)對水平走滑震源附近變形的影響，均質(zhì)球的介質(zhì)參數(shù)設(shè)置為震源附近參數(shù).該結(jié)果對應(yīng)公式(11)的第二項(xiàng)，可以發(fā)現(xiàn)ur,er θ,er φ分量的層狀結(jié)構(gòu)影響趨于穩(wěn)定，而對于其余6個(gè)分量，越靠近震源，受層狀結(jié)構(gòu)的影響越小.將其與對應(yīng)分量的漸近解疊加即可得到層狀地球模型下格林函數(shù).

圖7 層狀結(jié)構(gòu)對水平走滑震源的影響(ds=20 km)，縱軸為均質(zhì)地球模型和層狀地球模型下格林函數(shù)的差異，不同線型代表不同的變形深度

4 結(jié)論

本文基于自重、可壓縮和連續(xù)分層球體地球模型，利用龍格庫塔數(shù)值積分方法求解微分方程組.對傳統(tǒng)求解方法的積分路徑進(jìn)行調(diào)整，對微分方程的齊次部分進(jìn)行兩次相向積分，由地心積分至震源和由地表積分至震源，利用震源間斷條件求解震源上方和震源下方通解的組合系數(shù)，避免高階微分方程組求解中精度丟失的問題.

通過與均質(zhì)地球模型下微分方程組的解析解比較，驗(yàn)證了求解方法的正確性.并發(fā)現(xiàn)嚴(yán)格的無自重處理方式是縮小模型的密度值，而不是簡單的縮小重力值.另外，水平走滑震源受自重的影響較大，垂直傾滑受自重的影響較小，而自重主要影響垂直引張震源和水平引張震源的變形解y4，其中自重對垂直引張震源的影響超過100階.

最后比較了均質(zhì)地球模型和層狀地球模型下高階解矢量，發(fā)現(xiàn)兩者的比值趨于穩(wěn)定，當(dāng)均質(zhì)地球模型參數(shù)設(shè)置為震源所在層參數(shù)時(shí)，并且變形深度和震源處于模型同一層，兩者的比值趨于1.即與層狀結(jié)構(gòu)對淺源地震地表變形的影響類似，對于地球內(nèi)部地震變形，層狀結(jié)構(gòu)對震源附近變形的影響較小.基于該發(fā)現(xiàn)，結(jié)合均質(zhì)地球模型下漸近解，可解決分層地球模型下震源附近地震變形發(fā)散的問題.

致謝感謝中國科學(xué)院大學(xué)孫文科教授提供的球體位錯(cuò)理論配套程序.感謝兩位匿名審稿專家對本文的修改和完善提出了寶貴的意見.