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

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

1 中國地震局地震預測研究所， 北京 100036 2 中國科學院大學計算地球動力學重點實驗室， 北京 100049 3 中國地質大學(北京)地球物理與信息技術學院， 北京 100083

0 引言

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

數十年來，學者們基于不同幾何特征的地球模型發展了多種地震變形計算方法.例如，Okada(1985, 1992)通過總結整理前人的工作，給出了一套完整、簡潔、實用的同震形變計算公式，成為均勻半無限空間地震位錯理論的經典表達式.Wang等(2003, 2006)發展了半無限空間模型的地震位錯理論，考慮了地球的層狀結構，提高了近場地震變形的計算精度.Sun(1992)和Pollitz(1996)發展了層狀球對稱地球模型的地震位錯理論，考慮了地球曲率，提高了遠場地震變形的計算精度.

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

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

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

1 地震變形理論

1.1 地震變形的基本方程

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

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

(1)

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

(2)

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

(3)

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

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

1.2 數值積分求解微分方程組

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

圖2 兩種數值積分路徑示意圖

Y=X·β+Z,

(4)

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

(5)

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

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

(6)

(7)

(8)

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

(9)

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

2 地震變形計算結果

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

圖3 地球模型參數隨深度的變化

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

2.1 均質地球模型下解矢量

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

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

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

2.2 層狀地球模型下解矢量

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

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

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

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

3 震源附近地震變形研究

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

(10)

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

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

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

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

(12)

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

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

4 結論

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

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

最后比較了均質地球模型和層狀地球模型下高階解矢量，發現兩者的比值趨于穩定，當均質地球模型參數設置為震源所在層參數時，并且變形深度和震源處于模型同一層，兩者的比值趨于1.即與層狀結構對淺源地震地表變形的影響類似，對于地球內部地震變形，層狀結構對震源附近變形的影響較小.基于該發現，結合均質地球模型下漸近解，可解決分層地球模型下震源附近地震變形發散的問題.

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