同震重力变化对震源滑动参数的敏感性研究：以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年苏门答腊地震中的应用结果表明，在该震例上，同震重力变化对滑移量、滑动角、断层深度和断层倾角这四个独立震源滑动参数的敏感性依次降低.其中，滑移量和滑动角为影响该地震同震重力变化的高敏感性参数，在利用同震重力变化反演获取震源滑动参数时，需着重考虑这两个参数.另外，此方法同样适用于其他震例，但在不同震例中，同震重力变化对各滑动参数的敏感性将有所差别.

致谢感谢审稿专家提出的宝贵意见和建议.