2018年吉林松原MS5.7地震地磁异常分析

2021-11-10张志宏郭安宁李梦莹黄明威杨牧萍

科学技术与工程 2021年29期

张志宏, 郭安宁, 李梦莹, 黄明威, 杨牧萍

(1.辽宁省地震局，沈阳 110034; 2.中国地震局兰州地震研究所，兰州 730000)

经过几十年探索和实践，中外已发现很多与地震相关的前兆现象[1-3]。地磁学是地震前兆观测的重要手段，其主要研究内容是地球自身及周围空间磁场的起源、变化、分布规律及相关现象。地磁学通过实验、观测和理论研究，认识到了地球磁场的分布及变化规律。地磁背景场分布与变化的时、空都具有延拓性，且可以用物理意义明确的函数来表达。从而可以从空间有限个点的观测值推算其他的点值，用过去的变化预测未来的趋势。近年来，地震与地磁场的关系的研究取得了一些进展。目前中国地磁预测的地震方法有20多个，主要有：加卸载响应比法[4]、地磁逐日比法[5]、日变低点位移法[6]、谐波振幅比法[7]、地磁Z分量日变化空间相关法[8]和地磁每日一值空间相关法[9]。这6种不同的分析方法在预测地震三要素方面都有各自的优缺点：日变低点位移法预测强震发生的时间有较高的报准率，但是预测地震发生的地点区域有其局限性[6]。地磁Z分量日变化空间相关法对数据的质量要求很高；地磁加卸载响应比法和地磁逐日比法在预测地震发生的时间和震级具有独特的优势，震中范围预测存在不足[10-11]；地磁谐波振幅比法预测地震的震中大都在异常台站的周围，但是异常持续时间较长，不利于捕捉震前短期信号[12-13]。因此，对一次地震进行综合预测成为一个新的探索方向。郭安宁等[14]研究表明多方法分析异常可以较好地把握发震短临时间节点。解滔等[15]利用电磁多方法综合提取到了较好的2008年汶川MS8.0地震前电磁异常信号。

基于此，现分析2015—2019年东北地区固定地磁台连续观测波形数据，将地磁加卸载响应比法和地磁谐波振幅比法方法结合各自在预测地震优势，综合研究2018年5月28日松原MS5.7地震前后数据的变化特征，作为地磁多方法预测地震的震例。

1 震中附近地质构造背景

如图1所示为东北地区横跨中朝地台和天山—兴安地槽两大地质构造单元，以赤峰—开原断裂为界，北部的天山—兴安地槽自西向东可分为大兴安岭华力西褶皱带、内蒙华力西褶皱带、松辽拗陷和吉黑华力西褶皱带；南部的中朝地台自西向东可分为内蒙地轴燕山台褶皱带、华北断拗和胶辽台隆[16]。2018年5月28日吉林松原发生MS5.7地震，震源深度为10 km，震中位于松辽构造盆地[17]。

红色实心圆为东北地区公元1100年以来的6级以上地震；蓝色三角为分析地磁异常所使用的台站；黄色五角星为2018年吉林松原MS5.7地震的震中位置图1 东北地区地质构造与地震活动和地磁台站分布情况Fig.1 Geological structure, seismic activity and distribution of geomagnetic stations in Northeast China

2 计算原理与地磁资料处理

2.1 加卸载响应比法计算原理

地震的孕育和发生是一个复杂的非线性过程，因此，20世纪90年代初，尹祥础等[18]提出了定量反映非线性系统不稳定程度的参数——加卸载响应比Y，其可表述为

(1)

式(1)中：“+”和“-”分别为应力加载和卸载，某区域在特定时段内所有正地震(个数为N+)能量E的m次幂的和与负地震(个数为N-)能量的m次幂的和的比值定义为加卸载响应比Y。当孕震介质处于低应力的稳态弹性阶段时，Y接近于1。反之，Y随着失稳的过程逐渐增大。自此，加卸载响应比理论在各种地球物理参量得到了广泛应用[19]。冯志生等[20]研究发现地磁垂直分量Z日变相比其他分量震前异常更为显著。

地磁加卸载响应比计算公式为

(2)

式(2)中：RZ(max)为地磁垂直分量日变化幅度极大值；RZ(min)为出现之后的第一个极小值。

2.2 克里金插值法原理

1951年南非工程师Krige D G首先提出了克里金法(Kriging)，相比三角网线内插法、最小曲率法、反距离加权法和多项式拟合等常规插值法，克里金插值法在保持原场特征的同时更加合理地解释了地球物理异常[21]。

该方法是运用变异函数理论对有限区域内的区域变化的变量无偏最优化估计。其优点在于考虑待插值点与邻近观测数据点的位置的同时，对邻近点之间的空间位置也有所参考，使其误差接近于零[22]。

设定一个研究区D，区域化变量为Z(x,y)[Z(x,y)∈D]在采样i(i=1,2,…,n)点处的属性值为Zi(i=1,2,…,n)，待插值Z0点的插值结果则可以表述为

(3)

式(3)中：λi(i=1,2,…,m)为待确定权重系数；m为已知采样点。

根据Kriging原理，Z(x,y)满足二阶平稳假设，即

Cou{Zi,Zi+h}=E[Zi，Zi+h]-m2=C(h)

(4)

根据无偏性要求可推导得到权重系数方程为

(5)

式(5)中：h为滞后距离；Nh为滞后距离数量；Nk为求取半变异函数值的距离为h的样本对数目。求解权重系数可得λ到待插点估计值Z0。

克里金插值法的半变异函数可表述为

(6)

为了弥补研究区固定地磁台网离散的缺点，对加卸载响应比值进行密集规则网格化插值处理，得到直观的二维等值线空间分布。

2.3 地磁谐波振幅比计算原理

根据电磁感应理论，对于随时间周期变化的不均匀场源，视地球介质为均匀各向同性的平面导体的条件下，地磁谐波振幅比的定义为

(7)

式(7)中：YZHx(NS)和YZHy(EW)分别为南北向和东西向的地磁谐波振幅比；Z(ω)、Hx(ω)、Hy(ω)分别为地磁场的垂直分量、南北向水平分量及东西向水平分量的谱值。

地磁谐波振幅比也可以表达为磁测深视电阻率,即

(8)

式(8)中：ρax和ρay分别为北南向和东西向的深部电阻率；ω为圆频率;μ0为磁导率;λ为常数。

由式(8)可知，地磁谐波振幅比能够反映介质深部电阻率的变化。即地磁谐波振幅比算法具有明确的物理意义。

2.4 地磁数据处理方法

选取震中距700 km内11个地磁固定台的7套FHD磁力仪、2套GM4型磁通门磁力仪、3套M15磁力仪。地磁加卸载响应比异常属于中短期异常，因此，选取了2018年上半年的地磁Z分量日变化分钟值采样数据进行48阶傅氏拟合滤波后计算日变幅度，得到各个台站日变幅度逐日变化序列；其次根据日变幅度逐日变化序列计算各台站的地磁加卸载响应比值，提取所有台站加卸载响应比值，采用Surfer 进行克里金网格化插值处理，得到空间插值等值线数据。地磁谐波振幅比计算时间为2016年6月—2019年6月，谐波振幅比具体计算分为5个步骤：①计算每天地磁三分量采样数据谱幅度；②获得2个方向 5～65 min各周期的地磁谐波振幅比；③按 10 min的频带宽度计算各频带地磁谐波振幅比的频带均值，获得 10、20、30、40、50、60 min的地磁谐波振幅比；④计算 10、20、30、40、50、60 min的逐日变化序列；⑤采用窗长一年的滑动平均法消除各频带地磁谐波比逐日变化序列的年变化。

3 地磁异常分析

3.1 加卸载响应比异常分析

地磁场是地球内部磁性岩石及分布在地球内部和外部电流体系产生的各种磁场成分的叠加。由于磁场起源的不同，各种磁场成分的空间分布和时间变化规律也大不相同，按照场源的划分，此次地震场源为外源磁场。外源场起源于地表以上的空间电流系，这些电流体系随着时间变化较快，变化磁场的畸变具有短时性。因次，地磁加卸载响应比值异常表现为短时且跨越式的特征。

如表1所示计算了2008 年以来东北地区11个固定台站观测连续波形，剔除磁暴引起的日变幅畸变异常数据后得到东北地区加卸载响应比均值P，结果表明东北地区地磁场P值背景场区间是[1.53，1.75]。

表1 东北地区地磁加卸载响应比均值P

图2 东北地区2018年1—6月地磁加卸载响应比时序曲线Fig.2 Time series curve of geomagnetic loading and unloading response ratio in Northeast China from January to June 2018

研究表明，地震对应的地磁加卸载异常表现为震前中短期的高值异常，中国大陆地区异常阈值为3.0对应地震的概率较高[23]。如图2所示为2018 年5月28日MS5.7地震前6个月内的加卸载响应比时序曲线。图3(a)～图3(g)是2018 年5月28日MS5.7地震前6个月内加卸载响应比高于阈值3.0的7次异常的时空分布图。图3(a)中，2018年1月5日距离此次地震震中约137 km的吉林三岗台加卸载响应比超阈值异常，异常值为3.1。经调查异常日为磁静日，三岗台周边环境及观测系统正常。图3(b)中，1月18日内蒙锡林浩特及周边地区及辽宁铁岭地区加卸载响应比超阈值异常，其中辽宁铁岭台为3.1、锡林浩特为3.2。图3(c)中，2月11日德都地区加卸载响应比存在高值区，德都台计算结果为3.1。图3(d)中，2月14日辽宁营口、朝阳地区存在高值区，其中响应比值均为3.0。图3(e)中，2月20日出现了较大面积的高值异常区，东北东部大部分地区处于高值水平，大连、营口及铁岭地区的加卸载响应比值处于3.2～4.0的高值水平。图3(f)中，2月26日黑龙江德都比值为3.0。图3(g)中，4月24日东北西部地区出现了大范围的加卸载超阈值异常，其中距离震源最近的三岗台为3.8、最远的锡林浩特比值达到了4.9的高值，通河的比值为3.3，此次异常是2018年5月28日松原MS5.7地震前最后一次异常。

3.2 地磁谐波振幅比异常分析

选取研究区距离震中最近的三岗和铁岭台的Z、H、D三分量分钟值数据作为计算谐波振幅比。

图3 2018年5月28日松原MS5.7地震前地磁加卸载响应比异常时空演化图Fig.3 Temporal and spatial evolution of the abnormal geomagnetic load-unload response ratio before the Songyuan MS5.7 earthquake on May 28, 2018

如图4所示，铁岭台谐波振幅比分布在(0.150，0.425)区间，YZHx和YZHy两个方向 2016—2017年 6 个周期变化具有明显的同步性。2018年上半年 10、20 、30 min 周期开始呈下降趋势(图中红色箭头)，而 40、50、60 min 正好反向(图中黑色箭头)。2018 年 5 月 28 日松原MS5.7 地震发生在两个方向10、20 、30 min 周期趋势下降阶段，震后两个方向6个周期恢复上升准同步形态。

如图5所示，三岗台计算结果显示谐波振幅比分布在(0.046，0.210)区间内，谐波振幅比YZHx和YZHy两个方向 2016—2017年主要表现为准同步趋势上升特征。YZHx20、50 min 周期2018年初开始表现为趋势下降异常(图5中红色箭头)，YZHy50 min同样表现为趋势下降异常(图中黑色箭头)。2018年5月28日松原MS5.7地震后两个方向50 min 周期异常恢复，YZHx20 min 周期异常恢复后再次出现下降异常。

图4 2016—2019年铁岭台谐波振幅比Fig.4 Tieling station harmonic amplitude ratio from 2016 to 2019

4 讨论

震磁前兆信息的研究主要集中在基本磁场和变化磁场。按照场源的划分，研究区固定地磁台观测数据属于变化磁场。地震的孕育不仅是震源区孕震的过程，同时也和区域应力场及多种环境因素的变化密切相关。变化磁场是联系地震孕育的内、外环境的媒介，其“感应-触发”作用，对地震预报有着重要的意义。地磁谐波振幅比和地磁加卸载响应比可以理解为“感应”和“触发”的关系。

谐波振幅比计算结果表明，三岗台和铁岭台谐波振幅比两个方向的不同周期存在不同步异常，异常特征为下降—转折—恢复，2018年5月28日MS5.7地震发生异常恢复过程阶段，震后基本恢复年变形态。2018年5月28日松原MS5.7 震中位于NE向扶余—肇东断裂带和NW向第二松花江断裂带交汇部位，在研究东北地区 NW 向地质构造时，郯庐断裂带北段和两个分支断裂带是研究者关注的重点。如图6所示，三岗台及铁岭台从广义的地质构造来看都属于郯庐断裂带北段，此断裂带是中国大陆历史强震活跃区和现今地震重点危险区之一，第二松花江断裂 EW 向的东端与依兰—伊通断裂交叉。李君等[24]利用双差定位方法对此次地震进行了重新定位，并使用CAP方法求解了震源机制解，结果表明此次地震的发震断层为松花江第二断裂。松原地区应力场受到大尺度板块构造应力控制，同时受到区域构造活动的影响[25]。台站处于高导异常体边界附近是地磁谐波振幅比异常的机理解释[26]，由此可推断，三岗台和铁岭台可能处于同一高导异常体边界带上。

图5 2016—2019年三岗台谐波振幅比Fig.5 Sangang station harmonic amplitude ratio from 2016 to 2019

f1为松花江第二断裂；f2为扶余—肇东断裂；f3为依兰—伊通断裂；f4为密山—敦化断裂图6 2018年5月28日松原MS5.7地震地磁谐波振幅比异常台站分布图Fig.6 Distribution map of the abnormal stations of the Songyuan MS5.7 earthquake and geomagnetic harmonic amplitude ratio on May 28, 2018

2018年5月28日松原MS5.7地震前地磁加卸载响应比出现了7次超阈值异常。其中，1月5日除三岗台外研究区空间磁场日变处于正常背景。与其类似的还有2月11日和2月26日。1月18日锡林浩特台周围出现约2.5×105km2异常区，同时铁岭台出现了小范围异常。2月14日营口台、朝阳台周边存在小范围异常区。这5次跨越式的短时异常从此次地震的发震区域构造角度可靠性不足。而2月20日和4月26日的异常台站超过了3个，此次地震震中属于两次异常的异常控制区域。因此，2月20日和4月24日可能是“触发”阶段。

综上所述，2018年5月28日松原MS5.7地震区域应力场的“感应”可认为开始于2015年的下半年，“触发”时间节点分别为2018年的2月和4月末。