同震重力變化對震源滑動參數的敏感性研究：以2004年蘇門答臘地震為例

穆璇， 白永良*， 單新建， 張國宏， 王振杰

1 中國石油大學(華東)海洋與空間信息學院， 青島 266580 2 中國地震局地質研究所， 北京 100029

0 引言

地震是地殼運動引起的地球表層快速震動，是釋放地球內部能量的過程.快速、準確獲取震源參數，對于揭示發震機制具有重要作用(Mignan et al.，2015).獲取震源滑動參數的主流方法有地震波波形法(Shao et al.，2011；陳俊磊等，2020；何驍慧等，2020；趙翠萍等，2011)和形變觀測反演法(侯麗燕等，2020；邵志剛等，2011，2015；申文豪等，2019；張國宏等，2010).對于海底地震而言，因缺少有效的地震臺站、形變觀測等數據，導致上述方法無法反演海底震源參數；受時間尺度限制，針對地震波動時間很長的大地震或特大地震，這兩類方法也無法準確推測其震源參數(龔正，2016).

由于地震孕育和發生過程中，常伴有地表重力場的局部變化(Han et al.，2010；Heki and Matsuo，2010)，該變化的主分量與地震斷層運動引起的地表形變有關(付廣裕等，2015；賈宇鵬等，2015；談洪波等，2009，2011；燕乃玲等，2003；張國慶等，2015).因此，重力反演是獲取海底大地震初始震源參數的一個有效手段(Zhou et al.，2018；付廣裕等，2018；龔正，2016；劉寧等，2009；鄭增記等，2019).由于地球物理場固有的等效性，反演結果具有非唯一性.如能明確同震重力變化對震源不同滑動參數的敏感性特征，就可在一定程度上對反演結果加以限定，有助于提高震源參數的反演精度.然而，這一特征尚不清晰.

用于分析因變量對多個自變量敏感性的方法包括非參數統計法、曲線斜率法、敏感性系數法、敏感度函數法、正交試驗法、均勻設計法、灰色關聯度分析法以及人工神經網絡方法等(冉濤等，2018).其中，敏感性系數法和敏感度函數法適用于參數較少的情況.以此為基礎，章光和朱維申(1993)及王輝和陳衛忠(2012)針對不同參數具有不同物理量綱時無法直接比較敏感性的問題，提出了無量綱化的敏感性分析方法.在分析同震重力變化對震源滑動參數的敏感性時，各滑動參數也具有不同量綱.因此，本文將引入該方法，評估同震重力變化對震源不同滑動參數的敏感性.

本文將以2004年蘇門答臘地震初始震源參數為基準參數集，基于Okubo位錯理論正演計算地表同震重力變化，再引入無量綱化的敏感性分析方法，分析和比較不同滑動參數影響同震重力變化的敏感性，最終為基于同震重力變化反演初始震源參數提供指導.

1 同震重力變化對震源滑動參數敏感性分析方法

通過正演方法可以計算出震源不同滑動參數下的地表同震重力變化，進而可以分析同震重力變化對各滑動參數的敏感性.

1.1 震源不同滑動參數下的地表同震重力變化模擬

利用Okubo位錯理論可以根據震源滑動模型快速計算出地表同震重力變化.在彈性半無限空間直角坐標系o-x1x2x3中，設斷層寬度為W，斷層長度為L，斷層走向為φ，斷層傾角為δ，斷層質心深度為ξ3，斷層上盤相對于下盤的走滑、傾滑和張裂位錯分量分別為U1、U2、U3.定義滑動角θ為斷層上盤移動方向與斷層走向之間的夾角，滑移量D為滑動角方向的位錯分量.地表任意點(x1,x2,0)處的重力變化Δg可以表示為(Okubo，1992)

Δg(x1,x2)={ρG[U1Sg(ξ,η)+U2Dg(ξ,η)

+U3Tg(ξ,η)]+ΔρGU3Cg(ξ,η)}‖

-βΔh(x1,x2),

(1)

式中ρ、Δρ分別為介質密度、張裂紋內密度與介質密度之差；G為萬有引力常數；β=0.309×10-5m·s-1，為自由空氣重力梯度；ξ、η是為方便對斷層面積分而定義的代換變量；Sg(ξ,η)、Dg(ξ,η)、Tg(ξ,η)、Cg(ξ,η)是分別對走滑、傾滑、張裂和張裂填充物貢獻進行微分得到的系數，表達式由Okubo(1992)給出；Δh表示地表高程變化，對地表任意點(x1,x2,0)的數值可由下式計算:

+U3Th(ξ,η)]‖,

(2)

式中，“‖”為Chinnery記號，表示如下置換關系:

f(ξ,η)‖=f(x1,p)-f(x1,p-W)-f(x1-L,p)

+f(x1-L,p-W),

(3)

其他未提及參數的具體含義和計算過程，請參見Okubo(1992)的定義.

1.2 同震重力變化對震源滑動參數的敏感性分析

當利用αk表示第k個震源滑動參數、dg表示同震重力場振幅時，根據位錯理論式(1)，可以將dg和αk之間的關系抽象為

dg=f(α1,…,αk-1,αk,αk+1,…,αn).

(4)

(5)

根據式(5)繪出特性曲線dg-αk，其斜率可以表征同震重力變化對單一震源滑動參數αk擾動的敏感性.然而，這種由斜率直接表達的敏感性帶有量綱.例如，對斷層深度敏感性的量綱為μGal·km-1，對斷層傾角和滑動角敏感性的量綱為μGal/(°).因此，具有不同量綱各滑動參數的敏感性之間無法直接比較.

k=1,2,3,…,n

(6)

(7)

圖1 蘇門答臘地區構造綱要與地震分布(MW≥7.5)(數據來源：https:∥www.globalcmt.org/CMTsearch. html/[2021-05-14])Fig.1 Plate tectonics and earthquake (MW≥7.5) distribution in the Sumatra region (Data sources：https:∥www.globalcmt. org/CMTsearch.html/[2021-05-14])

2 用于同震重力變化正演的2004年蘇門答臘震源參數

蘇門答臘地區位于印度—澳大利亞板塊和歐亞板塊的俯沖型邊界上(鄧園浩等，2017)(圖1).圖1中震源球展現了全球矩心矩張量解中心(Global CMT)給出的1980年至今，該區域發生的MW≥7.5歷史地震.2004年12月26日，印尼蘇門答臘島以北海域發生了MW9.0級淺源逆沖型大地震，根據Global CMT發布的地震機制顯示，其震源中心位于3.09°N，94.26°E，是繼1960年智利地震后規模第二大的地震，也是GRACE重力衛星第一次記錄到重力變化的地震(劉泰等，2019).

影響地震引發地表同震重力變化的主要震源滑動參數有斷層位置、斷層走向、斷層傾角、滑移量及滑動角(龔正，2016；熊維等，2015).由于斷層平面位置和走向不影響地表同震重力變化幅度(龔正，2016)，所以本文不分析同震重力變化對這兩個參數的敏感性.最終，針對斷層質心深度(簡稱斷層深度)、斷層傾角、滑移量和滑動角這四個相互獨立的震源滑動參數進行地表同震重力變化正演和敏感性分析.

本文設計的方法適用于海底大地震發生后短時間內，基于快速測定的初始震源參數，研究某一震例同震重力變化對震源不同滑動參數的無量綱化敏感性排序，為正式反演該地震震源參數時的參數設置提供依據.因此，考慮應用環境及Okubo位錯理論的模型設置，本文根據2005年東京大學地震研究所(Yamanaka，2005)運用遠震體波波形法給出的2004年蘇門答臘地震震源參數(表1)，設置了單一有限矩形斷層模型.該模型與Okubo位錯理論相匹配，主要滑動參數于震后短時間內獲得，數據完整、模型簡單，適用于快速計算.將這組震源參數設置為基準參數集，統一設定四個參數圍繞基準值上下各波動40%，以此獲取用于敏感性分析的各參數變化范圍(表1).

表1 用于敏感性分析的2004年蘇門答臘地震震源參數基準值(Yamanaka，2005)與變化范圍Table 1 Reference value and variation range of source parameters (Yamanaka，2005) of the 2004 Sumatra earthquake for sensitivity analysis

圖2 斷層深度變化引起地表同震重力變化的正演結果Fig.2 Forward modeling results of surface coseismic gravity changes caused by depth changes

圖3 斷層傾角變化引起地表同震重力變化的正演結果Fig.3 Forward modeling results of surface coseismic gravity changes caused by dip changes

圖4 滑移量變化引起地表同震重力變化的正演結果Fig.4 Forward modeling results of surface coseismic gravity changes caused by slip changes

圖5 滑動角變化引起地表同震重力變化的正演結果Fig.5 Forward modeling results of surface coseismic gravity changes caused by rake changes

圖6 不同斷層深度引起的地表同震重力變化Fig.6 Surface coseismic gravity changes caused by different depths

圖7 不同斷層傾角引起的地表同震重力變化Fig.7 Surface coseismic gravity changes caused by different dips

圖8 不同滑移量引起的地表同震重力變化Fig.8 Surface coseismic gravity changes caused by different slips

圖9 不同滑動角引起的地表同震重力變化Fig.9 Surface coseismic gravity changes caused by different rakes

圖10 同震重力場振幅受斷層深度、斷層傾角、滑移量及滑動角影響的趨勢Fig.10 The variation of amplitude of coseismic gravity field affected by depth, dip, slip and rake

3 正演計算的正確性檢驗

根據2004年蘇門答臘震例斷層面長寬和斷層傾角的取值，以斷層質心的平面位置為中心，建立寬度為600 km、長度為1300 km的地面矩形模擬區.該模擬區域能夠囊括地震引發地表重力場有效變化的空間位置.對待研究的四個獨立震源滑動參數分別進行斷層深度加深、滑移量加大、斷層傾角和滑動角采用常用的特征角度值組合，以檢驗正演計算的正確性.

3.1 斷層深度變化引起地表同震重力變化正演計算的檢驗

為了使斷層深度的檢驗范圍不局限于淺源地震(<70 km)，設置斷層深度分別為10 km、50 km和90 km，其他參數取蘇門答臘地震震源的實際參數固定不變.如圖2所示，淺源地震引起的同震重力變化幅度大但影響范圍小，中源地震引起的同震重力變化幅度小而影響范圍大，符合深度變化導致地表重力變化的一般規律.

3.2 斷層傾角變化引起地表同震重力變化正演計算的檢驗

為了直觀分析斷層傾角變化引起的同震重力變化，同時控制斷層不突破地表，選擇滑動角為0°的走滑型斷層進行檢驗，且固定斷層深度為25 km.設置斷層傾角分別為特征角度值0°、45°和90°，其他參數取蘇門答臘地震震源的實際參數固定不變.如圖3所示，走滑型斷層引發的同震重力變化對稱分布于斷層兩側.隨著斷層傾角不斷增大，引起的同震重力變化范圍擴大.當斷層傾角達到90°時，同震重力變化的幾何分布符合一般規律，即呈內外兩個四象限分布的特征.

3.3 滑移量變化引起地表同震重力變化正演計算的檢驗

固定斷層深度為25 km，設置滑移量分別為20、40、60 m，其他參數取蘇門答臘地震震源的實際參數固定不變.如圖4所示，隨著滑移量呈倍數增加，引起的同震重力變化幅度穩步增大，且保持增幅一致，表明兩者間有明顯的線性相關性，正演結果符合位錯理論式(1)中滑移量與地表重力變化的線性關系.

3.4 滑動角變化引起地表同震重力變化正演計算的檢驗

固定斷層深度為25 km，設置滑動角分別為特征角度值0°、45°和90°，其他參數取蘇門答臘地震震源的實際參數固定不變.如圖5所示，當滑動角為0°時，正、負同震重力變化以垂直走向方向為軸，對稱分布于斷層寬邊兩側.隨著滑動角增大至45°，正、負同震重力變化逐漸向對面插入，幾何分布上出現了明顯的逆時針旋轉.直至滑動角增加到90°時，正、負同震重力變化以平行走向方向為軸，對稱分布于斷層長邊兩側.隨著滑動角從0°旋轉至90°，同震重力變化的幾何位置逆時針旋轉了90°，符合滑動角變化影響同震重力變化幾何分布的理論趨勢.

4 同震重力變化對2004年蘇門答臘震源滑動參數的敏感性

與正演計算的正確性檢驗相同，為囊括地震引發地表重力場有效變化的空間位置，以斷層質心的平面位置為中心，建立寬度為600 km、長度為1300 km的地面矩形模擬區.在敏感性分析的參數取值范圍內(表1)，依次計算地表區域同震重力變化隨斷層深度、斷層傾角、滑移量和滑動角這四個獨立參數變化的情況(圖6a—9a).當單個獨立參數變化時，其他參數取蘇門答臘地震震源的實際參數固定不變.為直觀分析不同滑動參數變化對同震重力變化的影響，在模擬區內劃定一條垂直于、斜交于、近似平行于斷裂走向的AA′、BB′和CC′ 剖面，得到沿三個剖面的同震重力變化隨各參數變化的情況(圖6b—9b).CC′ 剖面沒有劃定為絕對平行于斷裂走向，是為了同時獲取正、負地表同震重力變化在斷層走向首尾處的有效信息.

由圖6—9可知，當僅有斷層深度從6 km增加至14 km時，地表重力場變化范圍擴大(圖6a)，AA′與BB′剖面上同震重力變化振幅減小，CC′剖面上同震重力正變化數值加大、負變化數值減小、振幅變化不明顯(圖6b).當僅有斷層傾角從4.8°增加至11.2°時，地表重力場變化范圍無明顯改變(圖7a)，三條剖面上同震重力正變化數值減小、負變化數值加大、振幅變化不明顯(圖7b).當僅有滑移量從12 m增加至28 m時，地表重力場變化范圍顯著擴大(圖8a)，三條剖面上同震重力變化振幅顯著增加(圖8b).當僅有滑動角增加時，在相對90°幾乎對稱分布的67.2°和112°處，地表重力場變化范圍基本一致(圖9a)，AA′與BB′剖面上同震重力變化分布幾乎重合(圖9b)，但由圖9a中斷層走向的首尾處和圖9b 中CC′剖面上同震重力變化分布可以發現，地表重力場正負變化在空間上有逆時針旋轉趨勢，這是其他三個參數不具備的特性.當滑動角繼續增加至156.8°時，這一逆時針旋轉趨勢依然存在，且整個模擬區域重力場正負最值的空間位置連線已由垂直斷裂走向逐步旋轉至平行斷裂走向(圖9a).AA′與BB′剖面上的同震重力變化振幅明顯減小，CC′剖面上的同震重力變化振幅明顯增大(圖9b).

綜上，斷層深度、斷層傾角、滑移量和滑動角在用于敏感性分析的范圍內變化時，會導致幾條剖面上同震重力變化發生改變，且各自引起的重力變化幅度有所差異.特別地，當滑動角變化時(圖9)，同震重力場幾何分布發生旋轉，沿不同剖面具有不同的數值分布特征.因此，有限個剖面上的同震重力變化不足以說明整個重力場的變化情況.為準確分析2004年蘇門答臘地震同震重力變化對震源滑動參數的敏感性，本文選定該震例在整個模擬區域內的同震重力場振幅dg作為代表區域同震重力變化的系統特性，采用數值法擬合特性曲線(式(5))，大致了解同震重力變化受震源滑動參數擾動的影響，結果如圖10所示.

由系統特性曲線擬合結果(圖10)可知，當斷層深度從6 km增加至14 km的過程中，地震產生的區域同震重力場振幅減小了230.82 μGal，且減小的趨勢逐漸加快(圖10a).當斷層傾角從4.8°增加至11.2°的過程中，地震產生的區域同震重力場振幅增加了56.23 μGal，且增加的趨勢逐漸加快(圖10b).當滑移量從12 m增加至28 m的過程中，地震產生的區域同震重力場振幅增加了1799.29 μGal，且增加的趨勢呈線性(圖10c).當滑動角從67.2°至156.8°的過程中，地震產生的區域同震重力場振幅呈波動形式升降，且于90°處取得極大值、135°處取得極小值，相差了409.59 μGal，越靠近90°和135°處，其變化趨勢越平緩(圖10d).

利用式(7)計算目標滑動參數取值為2004年蘇門答臘地震實際震源參數時的敏感度因子(表2)，得到同震重力變化對該震例各參數的量化敏感性.結果表明，對于2004年蘇門答臘地震，同震重力變化對四個獨立參數的敏感性由高到低的順序為滑移量、滑動角、斷層深度和斷層傾角.其中，斷層傾角的敏感度因子比其他參數的敏感度因子小了一個量級，在四個參數中表現為最不敏感；而滑移量和滑動角表現出的敏感性明顯高于另外兩個參數的敏感性.由于本文方法沒有針對地震類型進行特定約束，是一個泛泛的方法，因此可在各種類型的地震中應用.但著重強調的是，上述結果是基于2004年蘇門答臘地震給出的，對于其他震例，該結論不可直接引用，需根據本文方法重新計算、排序.

表2 2004年蘇門答臘地震各參數的敏感度因子Table 2 Sensitivity factors of the 2004 Sumatra earthquake parameters

5 對震源參數反演的啟示

考慮到在淺源、低傾角、逆沖型的2004年蘇門答臘震例上，同震重力變化對滑移量和滑動角的敏感性比對斷層深度和斷層傾角的敏感性要高得多，建議在基于同震重力變化反演震源參數的過程中，加大對滑移量和滑動角的重視程度.例如，依據先驗信息，精細刻畫高敏感性參數的模擬區間.因對斷層深度和斷層傾角的敏感性較低，僅當滑移量和滑動角反演精度較高時，才能夠保證有效反演出斷層深度和斷層傾角.注意，這一結論僅代表同震重力變化對2004年蘇門答臘震例條件的敏感性，對于其他震例對應的基準參數集，結論會有所不同，但本文方法是適用的.

隨著衛星觀測技術的逐步發展，同震重力變化數據的時間和空間分辨率均會得到不斷提高.當震中位置距離臺站較遠時，可以使用地震臺站的地震波數據反演出初步的震源參數，并以此作為基準參數集.在基準值附近，利用本文方法分析同震重力變化對不同滑動參數的敏感性.依據敏感性差異，對各參數設定不同的反演權重，有助于得到更加準確的震源參數.

6 總結

本文設計了一套定量化表征同震重力變化對震源滑動參數敏感性的方法體系：基于位錯理論和初始震源參數，通過同震重力變化正演得到震源不同滑動參數下的區域地表同震重力場振幅，再利用無量綱化敏感性分析方法，定量評估同震重力變化對震源不同滑動參數的敏感性.上述方法在2004年蘇門答臘地震中的應用結果表明，在該震例上，同震重力變化對滑移量、滑動角、斷層深度和斷層傾角這四個獨立震源滑動參數的敏感性依次降低.其中，滑移量和滑動角為影響該地震同震重力變化的高敏感性參數，在利用同震重力變化反演獲取震源滑動參數時，需著重考慮這兩個參數.另外，此方法同樣適用于其他震例，但在不同震例中，同震重力變化對各滑動參數的敏感性將有所差別.

致謝感謝審稿專家提出的寶貴意見和建議.