典型地震断层错动导致的重力变化及2022门源6.9级地震重力变化分析

钟慧婷, 靳志同,2, 万永革,2

(1. 防灾科技学院, 河北 三河 065201;2. 河北省地震动力学重点实验室, 河北 三河 065201)

0 引言

确定断层的破裂方式对研究区的整个地震孕育过程研究及其震后地震趋势的发展均起到非常重要的作用,而重力场是反映地球介质密度变化和各种环境下动力学特征的最基本、最直接的物理量,重力-地震的联合反演具有显著的互补效果和更完备额的地质学解释[1],所以重力资料在反演震源参数和震源破裂中起着不可或缺的作用。如贾可可等[2]利用同震重力变化反演出日本“3·11”地震的震源破裂参数。郑增记等[3]基于GRACE记录的重力资料反演得到2012年苏门答腊地震断层参数。杨君妍等[4]通过GRACE重力卫星获取的数据反演出了地震位错Love数,并且得到了更准确的格林函数,反演出孕震区域局部构造信息。

除此之外,在地震预报领域,重力资料也发挥了一定的作用。祝意青等指出强震易发生在与构造活动有关联的重力变化正、负异常区过渡的高梯度带上,重力变化等值线的拐弯部位或四象限中心附近[5]。2001年昆仑山口西8.1级地震发生在重力变化高梯度带上,重力负值变化低于-40 μgal[6];2021年云南漾濞6.4级地震震中处于重力变化的高梯度带区域,判断漾濞地区地下活动的能量已经完全释放完,震后短期内无更大地震发生的可能[7]。这些研究表明重力资料可以为地球动力学过程、地震危险性判断提供重要信息。由于重力监测资料中包含了地震断层错动产生的重力变化值,而地震产生的重力变化会对揭示地震孕育等长期动力学行为造成干扰,因此将地震导致的重力变化扣除更有利于揭示地震孕育等地球动力学过程。

2022年1月8日,青海门源发生了6.9级地震,震中位于101.28°E、37.82°N,震源深度为10 km,造成了房屋建筑受损、兰新高铁受损停运等问题。此次地震位于青藏高原东北部边缘的祁连山地震带,由一系列的深大活动断裂组成。作为地震活动主体区域的祁连山次级地块,四周被深大走滑活动断裂包围,形成一个相对独立的较为活动的次级地块[8]。此次6.9级地震的发震构造既不是冷龙岭断裂,也不是托莱山断裂,而是发育在它们之间阶区中的道沟断裂[9-10],主要是由浅层的左旋走滑断层造成的。

(红色沙滩球表示本次地震主震的震源机制解,蓝色三角号表示研究区域的前兆台,红色线代表研究区域内不同的断层,左下角子图表示研究地区所在的位置)图1 研究区地形图Fig.1 Topographic map in the study region

近100年来,该地区已发生高于6.0级的地震20余次,其中发生的距离此次地震最近的一次大地震是2016年门源县6.4级地震。地震断层错动会导致地面重力变化,而中国地震前兆台网的重力仪可以观测到这种变化。这种变化信息对确定断层的错动方式、大小等都起到非常重要的作用。地震前后重力变化的原因一般有两方面:一是震前地震孕育过程中的质量迁移和地表形变引起的重力变化,二是同震位错引起周围物理密度的变化;前一种原因往往是重力变化的主要因素,后一种原因导致的重力变化一般在距离震中几十千米以外处只有约10 μgal量级且随着距离的增加迅速衰减[11]。

2022年青海门源发生6.9级地震,该地震发生在周围地震前兆台布设比较密集的地区,为了甄别重力前兆数据中的该地震错动引起的重力变化,从而为提取有意义的重力地震前兆提供基础数据,并了解该地震对周围物质运移和密度变化的影响,本文采用Okubo地震位错模型,分析了各种类型的地震错动引起的地表重力变化的特征,并计算了2022年青海门源县6.9级地震在周围地区引起的重力变化。

1 Okubo理论及计算公式

国内外许多学者对位错和断层引起的重力变化进行过深入研究。Stekettee[12]最早将位错理论引入地震学,Maruyama[13]给出了垂直和水平伸展断层引起的地表位移场的完整解析解,陈运泰等[14-15]讨论了结合半空间无限空间位错理论和地表形变资料进行震源反演的一般方法,Okada[16]总结了前人的工作,给出了一套完整简洁实用的同震变形计算公式。Okubo[11,17-18]研究了半无限空间介质内剪切和张裂断层的重力变化问题,导出了点源和有限断层的同震重力位和重力变化的解析表达式,Sun[19]解决了球形地球分层模型中地震位错的同震重力变化计算问题。李振洪等利用InSAR技术获取2022年青海门源地震的同震形变场,可以发现本次地震形变仅限局部地区,同时利用了弹性半空间的位错模型确定了本次地震事件的断层几何参数[20],鉴于2022年青海门源地震产生的重力变化仅限于局部地区,本文介绍Okubo的半空间线弹性介质中地震位错产生的重力变化计算方法。

Okubo理论主要计算了有限矩形断层的重力变化,根据最后的推导结果,用双竖约定符号标示为简洁模式[21]:

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

(1)

重力变化可以表示为:

Δg(x1,x2)={ρG[U1Sg(ξ,η)+U2Dg(ξ,η)+
U3·Tg(ξ,η)]+ΔρGU6g(ξ,η)}‖-
βΔh(x1,x2)

(2)

式中:G为万有引力常数;ρ为介质密度;Δρ=ρ′-ρ,ρ′为引张破裂后的填充物质的密度。S、D、T、C则表示走滑、倾滑、引张和引张破裂填充物的贡献,(Sg,Dg,Tg,Cg)是对(S,D,T,C)进行微分得到的,β≈0.309 mgal/m为自由空气重力梯度,Δh为高程变化,它们的具体表达式为:

(3)

(4)

(5)

Cg(ξ,η)=2I2cosδ-sinδ·lg(R+ξ)

(6)

q=x2sinδ-(d-x3)cosδ

(7)

(8)

同样,高程变化可以表示为:

(9)

具体表达式为:

(10)

(11)

(12)

其中:

(13)

(14)

(15)

I5(ξ,η)=2(1-2v)I1secδ

(16)

当cosδ=0时:

(17)

(18)

式中:δ为倾角;v为泊松系数;d为有限矩形断层的深度。

2 典型错动断层在其周围产生的重力变化规律

为了解典型断层引起的重力变化规律,本小节计算了单个左旋走滑断层、正断层、逆断层的重力变化情况,并且进行了分析。走向上单个断层长度为10 km,倾角方向上的断层宽度为10 km,在滑动方向上的位错为5 m,取介质密度为2 670 kg/m3,计算结果分析如下:

2.1 走滑断层重力变化特征

图2为单个走滑断层引起的重力变化分布图,走向为正东,倾角为90°,滑动角为0°,其中AA′为计算图3的剖面位置。图3为破裂中心位于不同深度的地表AA′剖面的重力变化曲线。

图2 走滑断层重力变化(AA′为图3的剖面位置)Fig.2 Gravity changes induced by strike-slip fault (AA′is the position of profile in Fig.3)

图3 同一剖面破裂中心位于不同深度的走滑断层产生的重力变化(黑色水平线代表零重力值)Fig.3 Gravity changes along profile induced by strike-slip fault with fracture center at different depths (The black horizontal line represents the zero gravity value)

由图2可以看出,垂直走滑断层产生的重力变化总体呈现出四象限对称分布;在近处,以断层迹线逆时针算起,一三象限重力变化为正值,二四象限为负值;在远处,一三象限重力变化为负值,二四象限为正值。图3显示,同一剖面重力变化曲线关于断层所在位置对称,重力变化大小随着远离断层而逐渐缩小;对于相同破裂尺度和滑动量的垂直走滑断层,破裂中心深度越深,震中附近的重力变化幅度越小,而自断层位置向远处的重力变化幅度衰减越慢。反之,震源破裂中心越浅,断层产生的重力变化局限在断层附近,峰值较大,但很快衰减。也就是说,相同破裂尺度和破裂量的浅表震源引起的重力变化值局限在震中的局部地区,而相同破裂尺度和破裂量的深部震源引起的重力变化会产生更大范围的重力变化,但重力变化值较小。

2.2 倾滑断层重力变化特征

图4与图7为单个正断层(滑动角为-90°)和逆断层(滑动角为90°)引起的重力变化分布图,走向为正东,倾角为45°。以正断层为例,图5为相同破裂中心位于不同深度的AA′剖面上的重力变化曲线,图6为破裂中心位于6 km处,不同倾角的断层破裂在AA′剖面产生的重力变化。逆断层同理,图像为7～9。

图4 正断层重力变化(AA′为图5的剖面位置)Fig.4 Gravity variation induced by normal fault (AA′ is the position of profile in Fig.5)

图5 同一剖面破裂中心位于不同深度的正断层产生的重力变化Fig.5 Gravity changes along profile induced by normal fault with fracture center at different depths

图6 同一剖面破裂中心位于6 km处不同倾角正断层重力变化Fig.6 Gravity changes along profile induced by normal fault with fracture center at a depth of 6 km and different dip angles

图7 逆断层重力变化(AA′为计算8的剖面位置)Fig.7 Gravity variation induced by reverse fault (AA′ is the position of profile in Fig.8)

由于倾滑断层在震源近处引起的重力变化值过大,为了更好地描述远处的重力变化特征,本文将重力变化大于5 μgal的值不显示,如图4和7。可以发现,在远处,正断层在震源南北方向产生的重力变化为正值,在东西方向产生的重力变化为负值,逆断层反之。图5和图8显示,上盘重力变化的数值总体上比下盘重力变化的数值小;对于相同破裂尺度和滑动量的倾滑断层,破裂中心深度越深,远处的重力变化最大值越小。

图8 同一剖面破裂中心位于不同深度的逆断层的重力变化Fig.8 Gravity changes along profile induced by reverse fault with fracture center at different depths

从图6和图9可看出,在倾角小于45°时,断层上盘产生的重力变化峰值较大,且随着断层的倾角的增大,在断层上盘重力变化波峰越大;在倾角大于45°时,断层下盘产生的重力变化峰值较大,且随着倾角的增大,断层下盘重力变化的峰值变小。当正断层倾角较小时,上盘重力变化主要以负值为主,下盘重力变化主要以正值为主;当倾角较大时,上盘重力变化主要为正值,下盘为负值;逆断层反之。

图9 同一剖面同一深度不同倾角逆断层重力变化Fig.9 Gravity variation along profile induced by reverse fault with fracture center at a depth of 6 km and different dip angles

3 2022年门源地震产生的重力场

本节采取Okubo理论,震中周围介质密度选取为2 670 kg/m3,地震破裂模型采用李振洪等[20]根据InSAR资料联合了升、降轨地表形变场反演的地震破裂模型;反演模型包含走向104°和109°两个断层段,走向104°的断层段长度为10 km,宽度为16 km,倾角为80°,滑动角为0°,最大滑动量为2.5 m;走向为109°的断层段长度为20 km,宽度为16 km,倾角为80°,滑动角为5°,最大滑动量为3 m[20]。将发震断层离散成子断层反演得到:该断层主要受左旋滑动控制,属于高倾角走滑型地震;子断层最大滑动量为3.5 m,出现在地下约4 km处,滑动主要集中在地下2～7 km的区域[20]。本文利用多个子断层叠加的方式,计算出重力变化情况,绘制出图10。

(黑点表示城市,蓝色三角号表示未能监测到重力变化的台站,红色三角号表示能够检测到重力变化的台站,黑色花瓣轮廓线是重力变化为1μgal的分界线,沙滩球是本次地震的震源机制解)图10 2022年1月8日门源地震后重力变化Fig.10 Gravity changes generated by Menyuan earthquake on January 8, 2022

图10显示,本次地震重力增加的区域范围要比重力减小的区域范围小。在数值上,重力变化负值区最大为-8.82 μgal,正值区最大变化为28.28 μgal。震中东北和西南方向的重力减小,震中西北和东南方向的重力增加;即在距离震中远处断层错动前方重力变化为正,后方重力变化为负,与前文叙述单个走滑断层引起的重力变化特征基本一致。

重力变化一般在距离震中几十千米以外处只有约10 μgal的量级,并随着距离的增加迅速减小[11]。本次门源地震断层错动导致的重力变化从震源处向四周衰减得很快,在数百公里之外几乎不发生重力变化。当前CG-5高精度重力仪在重力观测上分辨率可以达到1 μgal[22],图10中,黑色花瓣轮廓线是重力变化为1 μgal的分界线,花瓣内是能够采用重力仪监测到地震错动产生重力变化的地区,震中东北和西南方向、西北和东南方向的花瓣轮廓最远大约在100 km处,超出该花瓣范围的区域将难以监测到重力变化。

图10中有63个用三角号表示的前兆台,红色三角(门源台)标志表示能够监测到此次门源地震重力变化的台站,蓝色为不能够监测到重力变化的台站。需要指出的是,在已知台站中,有一个S9(门源台)前兆台站能够监测到本次门源地震的重力变化,其余台站处重力变化低于1 μgal,超过了重力仪监测的范围,所以未能监测到重力变化。即便在本次门源地震中只有一个前兆台能监测到重力变化,在监测和分析该地区长期重力变化从而揭示构造活动时,也应该扣除该地震造成的重力变化。

重力场的变化,能较好地反映地壳厚度的差异、地壳密度的变化和深部物质迁移等构造活动信息,重力场随时间变化与地震的形成和发展有着内在联系[23],所以重力资料在反演震源参数中也起着重要的作用。当然,单靠重力资料也是远远不够的,需要联合其他多种资料共同约束,从而得到更加准确的震源参数。

4 讨论与结论

为分析地震错动在地表导致的重力变化,本研究根据Okubo位错理论公式,计算分析了典型断层(走滑、倾滑)在地表的重力变化特征:(1)垂直走滑断层产生的重力变化总体呈现出四象限对称分布;对于同一剖面,重力变化大小随着远离断层而逐渐缩小;对于相同破裂尺度和滑动量的垂直走滑断层,深度越深,震中附近的重力变化幅度越小,而自断层位置向远处的重力变化幅度的衰减越慢。反之,震源越浅,断层产生的重力变化局限。(2)倾滑断层上盘重力变化的数值总体上比下盘重力变化的数值小;对于相同破裂尺度和滑动量的倾滑断层,破裂中心深度越深,远处的重力变化最大值越小。同时本文对2022年1月8日门源6.9级地震的破裂分布进行重力变化的计算和分析,得到:本次地震重力增加的区域范围要比重力减小的区域范围小;在数值上,重力变化负值区最大为-8.82 μgal,正值区最大变化为28.28 μgal,震中东北和西南方向重力降低,而震中西北和东南方向的重力增加。位于震中附近的门源台的重力增加3.40 μgal,能够被重力仪器监测到。

本文计算采用是李振洪等[20]的破裂模型,随着资料的公开和共享,更加精确的震源模型会改变本文的计算结果,但本文基于InSAR的数据反演模型抓住了主要的特征,即使后续有更为精确的破裂模型,对我们的计算结果不会较大影响。

余震也会对本文的计算结果也会有一定的影响。万永革等在计算兰德斯地震余震产生的位移场时发现,地震断层面及其附近余震产生了厘米量级的位移变化[24],对计算地表重力变化的数值影响较小;并且目前没有更小地震的震源机制和破裂模型数据,所以本文中没有计算更小地震对地表重力变化的影响。

除了上文所说的震源破裂模型、余震会对计算结果有影响,本文采取震源周围介质为均匀弹性半空间,密度统一为2 670 kg/m3,这也会对计算结果造成一定的影响。在计算门源地震后的重力变化时,同样还采用了密度为2 630 kg/m3、2 650 kg/m3、2 690 kg/m3进行计算,但是最后的计算结果影响都比较小,重力变化特征也基本一致。在不同密度的计算结果中,也是只有一个S9(门源台)前兆台站能够监测到本次门源地震的重力变化,其余台站处重力变化比较小。同样震后余滑或介质黏弹性松弛等产生的影响也较小,所以可以忽略不计。

虽然上述种种问题还有待于今后进一步解决,但本文在一级近似情形下给出了不同类型断层错动在地表产生的重力变化模式,并根据2022年门源地震的破裂模型计算得到了该地震产生的重力变化,这对采用重力资料反演该地震的破裂模型和研究长期重力变化从而揭示其他动力学过程是有益的。