2013年芦山 MS7.0地震同震地磁变化分析
宋成科1), 张海洋2)
1)中国地震局第一监测中心, 天津 300180
2)河北省地震局保定中心台, 保定 071000

〔作者简介〕 宋成科, 男, 1989年生, 2013年于中国地震局地壳应力研究所获固体地球物理专业硕士学位, 工程师, 现主要研究方向为地磁观测及其在地震预报中的应用, 电话: 18502288550, E-mail: songchk@126.com

摘要

同震地壳应力变化能够引起地磁场的变化。 文中利用成都地磁台连续观测资料分析了芦山 MS7.0地震前后的地磁场变化特征, 并根据芦山 MS7.0地震的同震破裂模型和岩石磁化强度-应力的线性模型, 计算了芦山 MS7.0地震同震引起的稳定地磁场变化。 分均值分析结果显示成都地磁台观测的同震稳定地磁变化非常微弱, 以武汉台和兰州台为参考, 同震稳定地磁变化分别为+0.09nT和+0.04nT; 秒观测值分析结果显示同震稳定地磁场变化为+0.02nT。 压磁效应计算结果显示, 断裂带下盘磁场垂直分量以向下变化为主, 断裂带上盘磁场垂直分量以向上变化为主, 水平分量的变化则比较散乱。 当应力响应系数为1×10-3MPa-1、 岩石磁化强度为1A/m时, 计算得到震中稳定地磁场总强度变化达5nT, 而距离震中约15km处为1nT, 在距离震中99km的成都地磁台处为+0.007nT。 以上结果为同震地磁变化研究提供了新的数据。

关键词: 芦山 MS7.0地震; 地磁观测; 压磁效应; 同震磁场变化
中图分类号:P315.72+1 文献标志码:A 文章编号:0253-4967(2020)06-1301-15
ANALYSIS OF COSEISMIC VARIATIONS IN MAGNETIC FIELD OF THE LUSHAN MS7.0 EARTHQUAKE IN 2013
SONG Cheng-ke1), ZHANG Hai-yang2)
1)The First Monitoring and Application Center, CEA, Tianjin 300180, China
2)Baoding Central Seismic Station of Hebei Earthquake Agency, Baoding 071000, China
Abstract

Stress perturbation related to earthquakes could cause the anomalous changes of magnetic field. Numerous observation data from high-precision magnetometers before and after earthquakes have revealed obvious coseismic changes in the geomagnetic field during earthquakes. In order to assess the ability of observing geomagnetic variations related to earthquakes with present geomagnetic stations, the coseismic geomagnetic variations corresponding to the 2013 MS7.0 Lushan earthquake are presented. Precise measurements of geomagnetic fields have been obtained at a continuous geomagnetic station, Chengdu(CDP), which is approximately 99km from the epicenter of the MS7.0 earthquake. The 12th Generation International Geomagnetic Reference Field is used to reduce the secular variation of geomagnetic field. The geomagnetic station, Wuhan(WHN)and Lanzhou(LZH), which are approximately 1 100km and 650km from the epicenter of MS7.0 earthquake respectively, are regarded as reference stations to reduce the diurnal variation of geomagnetic field. A corrected method based on observed data from CDP, WHN and LZH is also used to correct the effect of short-period geomagnetic variations of external origin. Minute average values of CDP from April 16, 2013 to April 23, 2013 is used to analyze the coseismic geomagnetic changes because the Σ Kp indices in this period are less than 16 which indicates that there are no strong magnetic disturbances, solar flares, or solar storms. The result reveals that coseismic geomagnetic variation at CDP is little, which is about +0.09nT and +0.04nT, with WHN and LZH as reference stations. Meanwhile, the coseismic geomagnetic variation is +0.02nT using the geomagnetic values sampled per second observed in 20s before and after earthquake from CDP.
The observed coseismic geomagnetic variations are generally interpreted in terms of the piezomagnetic effect in rocks, which are most probably generated by sudden stress changes resulting from earthquake rupturing. According to slip model of 2013 MS7.0 Lushan earthquake which is constrained by GPS data, the coseismic stress changes are calculated. Then, piezomagnetic magnetic field in the epicenter and surrounding area of the 2013 MS7.0 Lushan earthquake is calculated based on linear model of magnetization and stress. In the model, we assume that the initial magnetization is consistent with the present magnetic field, which has an inclination and declination of 47.2° and -1.7°(westward). The in-situ Curie depth is chosen as 30km because the magnetized rocks beneath 30km contribute negligibly to the ground geomagnetic field. The calculated piezomagnetic effect shows that the vertical component is downward and horizontal component is mainly southward resulting from stress release in footwall of faults. The vertical component is upward and the horizontal component is mainly northward in hanging wall of faults. When the stress sensitivity is 1×10-3MPa-1 and magnetization is 1A/m,the calculated piezomagnetic field is about 5nT, 1nT and +0.007nT in the epicenter, 15km away from epicenter and at CDP, respectively. According to coseismic geomagnetic variations acquired from observed data and calculated data, it is obvious that stations (sites) located close enough to the epicenter could record coseismic signals above several nano Tesla. This implies that it would be difficult to observe geomagnetic signals associated with moderate earthquakes with present geomagnetic stations in this area and we need more improved observations. The piezomagnetic model tended to underestimate the observed coseismic geomagnetic variations with stress sensitivity of 1×10-3MPa-1 and magnetization of 1A/m. In this study,a uniform initial magnetization of 1A/m is expected to be reasonable according to aeromagnetic data and ground-based geomagnetic data, which indicates a relatively small magnetic anomaly. In this case, the stress sensitivity of 2×10-3~3×10-3MPa-1 is suitable for explaining the observed coseismic geomagnetic variations. This result provides new data for coseismic geomagnetic variations. We can evaluate the initial magnetization, the stress sensitivity and slip model of earthquake if more observed data in near-field are available.

Keyword: Lushan earthquake; geomagnetic field observation; piezomagnetic effect; coseismic geomagnetic variations
0 引言

地磁场的变化与构造运动有关(Hao et al., 1982; Yamazaki, 2013), 随着构造应力的变化, 地壳介质的磁性也会产生变化, 从而引起地磁场的变化。 强震的发生会引起地壳变形和应力的显著变化, 进而产生地磁场变化(Yamazaki, 2013), 该现象可由压磁效应(即应力扰动引起的磁场变化)来解释(Bhattacharyya, 1964; Stacy, 1964)。

Sasai(1980, 1991, 1994)根据大量岩石压磁实验结果提出了岩石磁化强度与应力的线性关系, 并结合位错理论发展了压磁效应定量计算模型。 结合压磁效应计算结果与地面地磁观测结果开展分析, 可解释地震地磁异常及构造运动引起的地磁场变化(Johnston et al., 1994; Nishida et al., 2004)。 然而并非在所有中强地震中均能观测到震磁异常, 能否观测到震磁异常与震级及震中距有关。 Hattori等(2004)对震例进行了分析, 认为可观测到震磁异常的震中距与震级的关系可大致表示为0.025R< M-4.5(R为震中距, M为震级), 即M6.0和M7.0地震能够观测到异常的地磁台站的最大震中距分别为60km和100km, 对于一些震源深度较浅的地震, 地磁异常的分布范围可能更大。 在美国圣安德烈斯断层系和日本等地已观测到多个同震地磁异常, 如: Johnston等(1987)发现距离1986年North Palm Springs ML5.9地震震中3km的观测仪器记录到的同震地磁变化量为1.2nT; Mueller等(1990)发现距离1989年Loma Prieta ML7.1地震30km和50km的观测仪器记录到的同震地磁变化量分别为1.1nT和0.6nT; Johnston等(1994)发现距离1992年Landers MW7.3地震震中17km和24km处的观测仪器记录到的同震地磁变化量分别为1.2nT和0.7nT; Sasai等(1997)分析了日本伊豆半岛震群的多次地震同震地磁变化, 其中在1978年近伊豆大岛M7.0地震中所观测的同震变化量可达4.3nT; Utada等(2011)发现在距离2011年Tohoku M9.0地震震中200~500km范围内的台站均可识别出同震地磁变化, 最大变化量可达0.8nT, 最小变化量为0.01nT。 这些地面观测的同震地磁异常均能够用压磁效应解释, 这为衡量区域地壳介质磁性、 评估区域监测能力和认识同震磁场变化分析等提供了重要帮助(Johnston et al., 2006; Okubo et al., 2011; Utada et al., 2011; Yamazaki, 2013)。

青藏高原东缘是中国大陆构造运动较活跃的地区, 地壳物质的E向逃逸受到坚硬的华南地块阻挡, 在青藏高原东缘积累了较高的地壳应力, 强震危险性较高(黄禄渊等, 2017; 庞亚瑾等, 2019)。 近年来该区域先后发生了2008年5月12日汶川MS8.0地震(以下简称汶川地震)、 2013年4月20日芦山MS7.0地震(以下简称芦山地震)和2017年8月8日九寨沟MS7.0地震(以下简称九寨沟地震)。 目前该区域的地磁观测资料较多(主要包括地磁台站观测和流动地磁观测), 距离汶川地震和芦山地震震中100km范围内的成都地磁台先后观测到与汶川地震、 芦山地震和九寨沟地震相关的地磁异常(张凌等, 2016; 何畅等, 2017; 廖晓峰等, 2018; 解滔等, 2018)。 本文以芦山地震为例, 利用连续观测资料分析地磁台站观测的芦山地震前后的磁场变化, 并结合芦山地震的破裂模型计算了同震压磁效应, 最后将计算结果与观测结果进行对比, 分析该地震的震磁异常及地壳介质磁性参数。

图 1 区域构造背景和地磁观测点分布Fig. 1 Tectonic background and distribution of geomagnetic sites.

1 磁场观测数据分析

本文选取距离芦山地震震中99km, 观测质量较高的成都地磁基准台(以下简称CDP)的地磁场观测值分析台站观测到的同震磁场变化特征, 地震前后观测仪器的工作状态稳定, 未出现缺数现象。 CDP台的观测环境良好, 观测室内磁场梯度< 1nT/m, 观测室内日温差< 0.3℃, 年温差< 10.0℃, 相对湿度 < 85% , 均满足地磁观测规范要求。 每周进行2次基线值观测。 连续观测数据采样率为1s, 对秒观测值进行高斯滤波可获得分均值(廖绍欢等, 2016)。 选择兰州地磁基准台(以下简称LZH)和武汉地磁基准台(以下简称WHN)的同步观测资料作为对比数据。 连续观测数据来自国家地磁台网中心。

地面磁场测量值主要包含稳定变化的主磁场、 地壳磁场和瞬变的日变场(包括平静变化和扰动变化), 在分析地震地磁异常时, 需消除主磁场和日变场等有规律的变化成分的影响。 这里选择第12代国际地磁参考场(IGRF-12)作为主磁场模型(Thebault et al., 2015a, b), 用原始观测数据减去IGRF-12主磁场以消除地磁场稳定变化。 日变场与地磁扰动的强度有关, 这里使用ΣKp描述地磁场扰动, ΣKp为每天的8个Kp指数(3h磁情指数)之和, 数据来自世界地磁数据中心(① http://wdc.kugi.kyoto-u.ac.jp。)。

地磁场分均值的连续变化如图2a所示。 3个台站的地磁场变化特征基本一致, 最明显的特征为地磁场的规则日变化及扰动变化, ΣKp越大的时段地磁场扰动越明显。 通常认为午夜时段外源及人为磁场干扰较少, 从午夜均值(00—03时)的变化曲线(图2b)可以看出, 3个台站的磁场变化趋势较为一致。 受强烈外源扰动的影响, 部分时段的磁场变化幅度略有差别, 如0430日和0501日(图2b中的红色箭头所示), CDP、 WHN和LZH的变化量分别为28.9nT、 28.2nT和20.7nT。 地震前后, CDP、 WHN和LZH所有观测日午夜均值的平均值分别为-5.0nT、 -4.9nT、 -3.7nT和-2.8nT、 -2.7nT和-0.6nT。

图 2 芦山地震前后的地磁场变化
a 3个台站的地磁场总强度观测分均值曲线; b 3个台站的地磁场总强度观测午夜均值(00—03时)曲线; c ΣKp, 蓝色箭头为0424日的地磁场扰动, 红色箭头为0501日的地磁场扰动
Fig. 2 Variations of magnetic field before and after Lushan earthquake.

为识别微弱的震磁异常, 需最大限度消除外源场的影响, 常用的技术手段为计算2个台站同步观测的差值(Johnston et al., 2006; Takla et al., 2011; Utada et al., 2011)。 该方法的基本假设是目标台站和参考台站的地磁场变化一致, 且参考台站的地磁场不受地震影响, 因此这里选用距离芦山地震震中分别为650km和1 100km的LZH和WHN作为参考台, 分别计算CDP和LZH、 WHN午夜均值的差值。 由于3个台站距离较远, 外源场变化并不完全相同(Johnston et al., 2006), 这里也采用适用于距离较远台站的地磁数据校正方法对分均值数据进行校正。 该方法使用参考台站的观测值重构目标台站的地磁场(即重构值), 认为目标台站的实际观测值与重构值的残差即是构造运动产生的磁场变化(Utada et al., 2011)。 对于有限长度的离散时间序列, 任意台站的观测值可用参考台站的水平分量和垂直分量的卷积表示:

Fk(i)=j=-LMpkjH0i-j+qk(j)Z0(i-j)(1)

式中, Fk为目标台站总强度重构值, H0为参考台水平分量观测值, Z0为参考台垂直分量观测值, pkqk为校正系数, 可通过最小二乘法求解, LM表示重构一个目标台站数据所使用的参考台站的数据量, 以不少于1个观测日的数据量为宜。 假设数据序列包含n个数据, 则i为该数据序列的第i个位置, i的取值范围为[M+1, n-L], j为第i个位置前或后的第j个位置, j的取值范围为[-L, M]。

图3为直接差值方法和和校正后的地震前后的CDP稳定地磁场变化。 从图3a和3b可看出, 经校正后地磁场变化更加平稳。 以LZH台作为参考, CDP稳定地磁场变化的标准差分别为2.8nT和1.0nT。 直接差值结果表明, ΣKp指数较大的时段地磁场的变化明显较大; 而校正后的计算结果显示较大的地磁场扰动得到抑制。 标准差越小, 则说明该时段内外源场扰动产生的影响越小, 因此认为通过校正能够有效降低外源场的影响。 下文的分析将仅使用校正后的结果。

图 3 芦山地震前后CDP与WHN、 LZH总强度午夜均值(00—03时)的差值
a 20130101—20130731时段以LZH作为参考台站的差值和校正结果; b 20130101—20130731时段以WHN作为参考台站的差值和校正结果
Fig. 3 Difference of the average value of total intensity during the period from 00:00 to 03:00 hours between CDP and WHN before and after Lushan earthquake.

由于0424日的地磁场扰动强烈(如图2b中蓝色箭头所示), 因此仅使用0420—0423时段(ΣKp< 16)午夜值的均值表示震后的地磁场水平。 与之对应, 用0416—0419时段(ΣKp< 16)午夜值的均值表示震前的地磁场水平, 二者的差值为同震稳定磁场变化, 结果见图 4 和表1。 以WHN作为参考台, CDP的同震稳定地磁变化为+0.09nT; 以LZH作为参考台, CDP的同震稳定地磁变化为+0.04nT。 由于使用分均值获得的CDP同震稳定地磁场变化较小, 为了确定该变化是否准确, 下面将对地震前后短时间内的秒观测值进行分析。

图 4 20130416—20130423时段午夜均值的差值
a 以LZH为参考台; b 以WHN为参考台。粗虚线表示地震前后午夜均值差值的平均值, 箭头指示芦山地震
Fig. 4 Difference of midnight-mean value from April 16, 2013 to April 23, 2013 with reference LZH(a) and WHN(b).

表1 芦山地震前后的地磁场总强度变化 Table1 Variations of total intensity of magnetic field before and after Lushan earthquake

根据Okubo等(2011)的研究, 地震地磁变化可先于地震波被观测到, 随着地震破裂的发展和地震波的传播, 地磁场出现动态变化, 地震波到达台站后地磁场趋于稳定, 因此需使用地震波到达后的稳定观测数据进行同震稳定磁场变化分析。 图 5 为CDP在芦山地震前、 后共60s内观测的地磁场变化。 地震发生后16~20s, 由地震波到达所致, 地磁观测数据出现剧烈变化。 地震后20s地磁场变化趋于稳定, 故用地震20s后的观测数据代表震后稳定地磁场变化。 若假设极短时间内地磁场恒定或呈线性稳定变化, 对比地震前、 后短时间内的地磁场即可获得同震稳定地磁变化, 如图 5 所示。 地震前后, CDP处的地磁场变化非常微弱(< 0.1nT)。 用地震波到达之后10s和20s的观测值表示震后地磁场, 并用地震前10s和20s的观测值表示震前地磁场, 获得的同震稳定地磁变化分别为+0.02nT和+0.02nT。

图 5 芦山地震发生前、 后60s内CDP磁场总强度的观测变化曲线
蓝色虚线为用10s观测数据表示的地震前、 后的磁场
Fig. 5 Variations of the total intensity in sixty seconds before and after Lushan earthquake.

根据前文对分均值的分析, 以WHN作为参考台, CDP的同震稳定地磁变化为+0.09nT; 以LZH作为参考台, CDP的同震稳定地磁变化为+0.04nT。 秒观测值分析结果显示, CDP的同震稳定地磁变化为+0.02nT。 2种分析结果在符号和量级上均一致。 综合上述使用不同数据得到的同震稳定磁场的变化结果可知, CDP观测到芦山地震同震产生的稳定地磁变化非常微弱(< 0.1nT)。

2 压磁效应的计算方法及结果
2.1 同震压磁效应的计算方法

地壳应力变化会改变岩石的磁化率和磁化强度(Nagata, 1969; 郝锦绮等, 1989), 进而引起地磁场的变化。 结合压磁理论及芦山地震的同震破裂模型, 可计算地震引起的应力变化及地磁场变化。

本文采用PSGRN/PSCMP开源程序(Wang et al., 2006)计算芦山MS7.0地震的同震应力, 其中地震破裂滑动分布选择Jiang等(2014)利用GPS观测资料反演获得的同震断层破裂滑动分布模型。 将龙门山断裂带及周边地区复杂的岩石圈结构简化为均匀分层结构, 岩石圈各分层的力学参数参考Crust1.0。 之后, 进一步计算芦山地震的同震破裂引起的区域岩石圈磁化强度的变化及磁场变化。

根据岩石磁学实验结果, 应力与磁化强度关系可近似表示为(Sasai, 1980, 1991; Hao et al., 1982)

ΔJ=32βS·J(2)

其中, Δ J为磁化强度的变化, J为岩石初始磁化强度, β 为应力响应系数, δ kl为克罗内克函数(k=l时, δ =1; kl时, δ =0), kl分别为坐标轴, S为偏应力张量。

Skl=σkl-13δkl(σ11+σ22+σ33)(3)

根据胡克定律, 应力与位移的关系可表示为

σkl=λδkl·u+μukxl+ulxk(4)

其中, σ kl为应力, λ μ为拉梅常数, u为位移。

·u=ux+uy+uz=13λ+2μ(σ11+σ22+σ33)(5)

则偏应力张量与位移关系可表示为

Skl=λδkl·u+μukxl+ulxk-13δkl(3λ+2μ)·u=μukxl+ulxk-23δklμ·u(6)

将式(6)代入式(2), 即可通过位移求得磁化强度变化(Hao et al., 1982; Okubo et al., 2004):

ΔJkl=βμ32ukxl+ulxk-δkl·uJk(7)

磁场强度H和标量磁位W的关系为

H=-W(8)

考虑在直角坐标系中的一个微元dv(或dxdydz)(图 6), 其磁化强度变化为Δ J, 则其在任意一点P(x0, y0, z0)处的标量磁位dW可用式(9)表示(Sasai et al., 1978):

d(W)= ΔJ·ρρ3dv=|Δ Jlx0-x+my0-y+nz0-zρ3dxdydz (9)

其中, (l, m, n)为Jx, y, z轴的方向余弦, ρ 为微元到P点的矢量, |ρ |为矢量ρ 的模, |Δ J|为矢量Δ J的模, |ρ |= r0-r= x0-x2+y0-y2+z0-z2

图 6 直角坐标系中的微元Fig. 6 Element body in rectangular coordinate system.

整个磁化体在P点的标量磁位可表示为(Sasai et al., 1978; Nishida et al., 2004; Okubo et al., 2004)

W=|ΔJ|·lx0-x+my0-y+nz0-zρ3dxdydz(10)

对式(10)右边求导数即可获得空间任意点的变化磁场N向分量Δ X、 E向分量Δ Y和垂向分量Δ Z:

ΔX=0H-+-+|ΔJ|(lU1+mU4+nU5)dxdydzΔY=0H-+-+|ΔJ|(lU4+mU2+nU6)dxdydzΔZ=0H-+-+|ΔJ|(lU5+mU6+nU3)dxdydz(11)

其中, (l, m, n)为Jx, y, z轴的方向余弦, H为计算深度。

U1=2x0-x2-y0-y2-z0-z2ρ5U2=2y0-y2-x0-x2-z0-z2ρ5U3=2z0-z2-x0-x2-y0-y2ρ5U4=3(x0-x)(y0-y)ρ5U5=3(x0-x)(z0-z)ρ5U6=3(y0-y)(z0-z)ρ5(12)

基于以上方法可计算出同震稳定磁场变化。 参考实验室结果和理论计算的应力响应系数的量级为10-3MPa-1(Nagata, 1969; Revol et al., 1977; Martin, 1980), 本文沿用诸多震例计算使用的应力响应系数数值1×10-3MPa-1。地表岩石的磁化强度为0.1~2A/m, 随着深度增加磁化强度也相应增加(Johnston et al., 2006), 较多震例研究选用的磁化强度为1~2A/m, 本文计算采用的磁化强度为1A/m。根据熊盛青等(2016)给出的大地热流分析结果, 该区域居里等温面深度约为30km, 且30km以下岩石的磁性变化并不会对地表产生明显的压磁效应(Carmichael, 1977), 因此模型的计算深度选取为30km。 磁偏角和磁倾角根据该区域多个地磁测点的绝对观测结果确定。 计算采用的模型介质磁性参数见表2

表2 压磁模型介质参数 Table2 Parameters in piezomagnetic model
2.2 芦山地震的同震磁场变化结果

根据表2给出的参数建立压磁模型(正方向分别为N向、 E向和垂直向下), 得到的芦山地震压磁效应计算结果如图 7 所示。 计算时采用的区域岩石磁化强度的方向与现今的地磁场方向一致, 即垂直分量向下, 水平分量以N向分量为主, W向分量较小。 从计算结果可以看出, 断裂带NW侧(即断层上盘一侧)地磁场的垂直分量变化以向上为主, 断裂带SE侧(即断层下盘一侧)产生的地磁场垂直分量变化以向下为主(图7b); 震中区同震磁场的水平矢量变化比较散乱, 既有较大的S向变化, 也有较大的N向变化(图7c)。 随着震中距的增加, 同震磁场的变化逐渐变弱。

图 7 芦山地震的同震压磁效应计算结果
a 芦山地震的同震滑动在水平面投影的分布图(Jiang et al., 2014); b 地表磁场垂直分量的同震变化; c 地表磁场水平分量的同震变化; d 地表磁场总强度的同震变化。 箭头代表磁场水平分量的变化方向, 矩形框表示断层滑动分布的范围, LC表示流动地磁观测点
Fig. 7 Calculated piezomagnetic field of Lushan earthquake.

从计算的地磁场总强度变化结果看(图7d), 由于震中的滑动量最大, 震中附近地磁场总强度变化量最大, 高达5nT, 距离震中15km处可达约1nT。 然而, 计算获得的CDP台站处的总强度变化较微弱, 仅为+0.007nT; 在距离震中26km的流动地磁测点处(LC)的总强度变化为-0.25nT。 显然, 明显的压磁效应仅局限在震中周边区域, 这主要是由于芦山地震震级和破裂尺度较小所致。

3 讨论

在计算地震压磁效应时, 必须选取合适的岩石磁性参数。 许多震例研究选取的应力响应系数(1×10-3MPa-1)和磁化强度(1~2A/m)能够解释同震地磁变化。中国地震局流动地磁监测技术团队基于2008—2009年和2013—2015年中国近千个流动地磁的观测结果和几十个地磁台站的观测数据, 在对其进行日变通化改正、长期变化改正并消除主磁场后, 获得了2010.0年代(2010年1月1日0时0分0秒)中国大陆岩石圈磁异常和2015.0年代(2015年1月1日0时0分0秒)中国大陆岩石圈磁异常。从2010.0年代和2015.0年代的岩石圈磁异常可以看出, 研究区的磁异常< 40nT, 比其东部区域的岩石圈磁异常小很多(图8), 说明该区域的岩石磁化强度并不高。实验室相关研究结果表明, 岩石应力响应系数随岩石孔隙率、 颗粒大小等的增加而增大(Kean et al., 1976; Hamano, 1983), 岩石应力响应系数有可能> 1×10-3MPa-1

图 8 区域岩石圈磁场异常
a 2010.0年代区域岩石圈磁异常; b 2015.0年代区域岩石圈磁异常。 黑色五角星为芦山地震震中; 灰色矩形框表示断层滑动分布范围
Fig. 8 Lithosphere geomagnetic anomaly in 2010.0 epoch(a)and 2015.0 epoch(b).

地磁观测结果和压磁效应计算结果的对比分析可为岩石磁性参数的选取提供约束和参考。 对于地磁场的N向分量Δ X, 将式(2)代入式(11), 可得:

ΔX=0H-+-+32βJ|S·e|lU1+mU4+nU5dxdydz(13)

其中, ‖ 表示向量的模, eJ的方向向量。 岩石磁性参数β J在研究区是均匀的, 即β J为常数, 故由式(13)可知Δ Xβ J成正比, 则Δ Y和Δ Z均与β J成正比。 磁场总强度FXYZ 3个分量的关系可用三角函数表示, 因此Fβ J成正比, 比值与偏应力S有关。 由于通过破裂模型计算的同震应力可认为是实际应力变化, 则实际观测的磁场变化Fobserved与实际的(β J)actual之比等于压磁效应计算的磁场变化Fcalculted与计算时假设的(β J)assumed之比, 即(Currenti et al., 2007; Yamazaki, 2013):

(βJ)actual=FobservedFcalculted×(βJ)assumed(14)

观测结果与模型计算结果列于表3。 压磁效应计算所得CDP处的磁场变化约为+0.007nT, 观测所得磁场变化为+0.02~+0.09nT, 则实际的β J可能是压磁效应计算时输入的参数的3倍甚至更大。 在流动地磁观测点处的压磁效应计算磁场变化为-0.25nT, 而流动观测磁场变化为-0.60nT, 则实际的β J可能是计算模型输入的参数的2.4倍。 综合以上2个结果, 为了解释CDP和流动地磁测点观测的同震变化, 实际的β J需要比本文模型使用的数值更大, 即若磁化强度选用1A/m时, 区域应力响应系数为2×10-3~3×10-3MPa-1比较合理, 对应于该岩石磁性参数, 震中周围区域压磁变化将为本文图7 中所示计算结果的2~3倍, 即震中附近的变化值最大将超过10nT。

表3 模型计算和观测的同震地磁场总强度变化 Table3 Model calculated and observed coseismic variations in the geomagnetic total intensity

本文计算压磁效应时使用了均匀的磁化强度和应力响应系数, 这种简化必然导致计算结果与实际情况存在偏差。 非均匀的磁化强度能够增强地磁变化(Oshiman, 1990), 计算时使用尽可能详细的区域磁化强度结果可使计算结果与实际变化结果更为接近, 但目前很难准确获得精细的磁化结构。

4 结论

本文以芦山地震为研究对象, 利用青藏高原东缘的地磁场连续观测资料分析了地震前后的磁场变化特征, 结合芦山地震的破裂模型定量计算了区域同震压磁效应, 并对同震稳定地磁场变化与计算结果进行了对比分析, 主要结论如下:

(1)地磁台站分均值观测资料分析结果显示成都地磁台观测到的芦山地震同震稳定磁场变化非常微弱, 以武汉台和兰州台作为参考, 同震稳定地磁变化分别为+0.09nT和+0.04nT; 使用秒观测值分析结果显示成都地磁台观测到的芦山地震同震稳定磁场变化仅为+0.02nT。

(2)同震压磁计算结果分析显示芦山地震在震中区域产生的压磁效应高达5nT, 距离震中15km处压磁效应也达1nT, 而成都地磁台处的压磁效应仅为+0.007nT, 说明芦山地震的同震压磁变化主要集中于震中附近, 距离震中较远的区域地磁变化非常微弱。

(3)计算得到的压磁场变化量与实际观测的磁场变化量不一致, 这主要与模型计算时选取的区域磁性参数有关, 当区域应力响应系数取为2×10-3~3×10-3MPa-1 时, 压磁效应能够解释目前观测到的芦山地震同震地磁异常。 此外, 在台站观测基础上, 结合压磁效应的定量计算可以帮助研究者由点到面地认识区域同震压磁异常。

致谢 本研究使用了国家地磁台网中心的地磁观测数据和中国地震局流动地磁监测技术团队的数据结果; 审稿专家对文章提出了宝贵的建议。 在此一并表示感谢!

参考文献
[1] 陈斌, 袁洁浩, 王璨, . 2017. 流动地磁监测数据处理流程[J]. 地震研究, 40(3): 335339.
CHEN Bin, YUAN Jie-hao, WANG Can, et al. 2017. Data processing flowchart of Chinese mobile geomagnet monitoring array[J]. Journal of Seismological Research, 40(3): 335339(in Chinese). [本文引用:1]
[2] 郝锦绮, 黄平章, 张天中, . 1989. 岩石剩余磁化强度的应力效应[J]. 地震学报, 11(4): 381391.
HAO Jin-qi, HUANG Ping-zhang, ZHANG Tian-zhong, et al. 1989. The stress effect on remanent magnetization of rocks[J]. Acta Seismologica Sinica, 11(4): 381391(in Chinese). [本文引用:1]
[3] 何畅, 廖晓峰, 祁玉萍, . 2017. 2017年8月8日九寨沟 MS7. 0地震前成都台地磁谐波振幅比异常分析[J]. 中国地震, 33(4): 575581.
HE Chang, LIAO Xiao-feng, QI Yu-ping, et al. 2017. Study on anomalous variation of the geomagnetic harmonic wave amplitude ratios of the Chengdu seismic station before the MS7. 0 Jiuzhaigou earthquake on August 8, 2017[J]. Earthquake Research in China, 33(4): 575581(in Chinese). [本文引用:1]
[4] 黄禄渊, 张贝, 瞿武林, . 2017. 2010智利Maule特大地震的同震效应[J]. 地球物理学报, 60(3): 972984.
HUANG Lu-yuan, ZHANG Bei, QU Wu-lin, et al. 2017. The co-seismic effects of 2010 Maule earthquake[J]. Chinese Journal of Geophysics, 60(3): 972984(in Chinese). [本文引用:1]
[5] 廖绍欢, 李雪浩, 胡俊明, . 2016. 成都地震基准台FHDZ-M15仪与GM4仪观测资料对比分析[J]. 四川地震, (4): 1520.
LIAO Shao-huan, LI Xue-hao, HU Jun-ming, et al. 2016. Comparative analysis of observed data recorded by both geomagnetic instrument FHDZ-M15 and GM4 at Chengdu seismic stand ard station[J]. Earthquake Research in Sichuan, (4): 1520(in Chinese). [本文引用:1]
[6] 廖晓峰, 何康, 张明东, . 2018. 基于平滑伪魏格纳-维勒时频分析的地磁数据研究[J]. 地震, 38(1): 107116.
LIAO Xiao-feng, HE Kang, ZHANG Ming-dong, et al. 2018. Time-frequency analysis of geomagnetic data based on smooth pseudo-Wigner-Ville distribution[J]. Earthquake, 38(1): 107116(in Chinese). [本文引用:1]
[7] 庞亚瑾, 程惠红, 张怀, . 2019. 青藏高原东缘岩石圈结构对现今地表垂向运动影响的数值分析[J]. 地球物理学报, 62(4): 12561267.
PANG Ya-jin, CHENG Hui-hong, ZHANG Huai, et al. 2019. Numerical analysis of the influence of lithospheric structure on surface vertical movements in eastern Tibet[J]. Chinese Journal of Geophysics, 62(4): 12561267(in Chinese). [本文引用:1]
[8] 解滔, 刘杰, 卢军, . 2018. 2008年汶川 MS8. 0地震前定点观测电磁异常回溯性分析[J]. 地球物理学报, 61(5): 19221937.
XIE Tao, LIU Jie, LU Jun, et al. 2018. Retrospective analysis on electromagnetic anomalies observed by ground fixed station before the 2008 Wenchuan MS8. 0 earthquake[J]. Chinese Journal of Geophysics, 61(5): 19221937(in Chinese). [本文引用:1]
[9] 熊盛青, 杨海, 丁燕云, . 2016. 中国陆域居里等温面深度特征[J]. 地球物理学报, 59(10): 36043617.
XIONG Sheng-qing, YANG Hai, DING Yan-yun, et al. 2016. Characteristics of Chinese continent Curie point isotherm[J]. Chinese Journal of Geophysics, 59(10): 36043617(in Chinese). [本文引用:1]
[10] 张凌, 杜爱民, 郎雪, . 2016. 汶川地震前的地磁场日变化特征分析[J]. 地球物理学报, 59(3): 952958.
ZHANG Ling, DU Ai-min, LANG Xue, et al. 2016. Characteristics of geomagnetic regular diurnal variation before Wenchuan earthquake[J]. Chinese Journal of Geophysics, 59(3): 952958(in Chinese). [本文引用:1]
[11] Bhattacharyya B K. 1964. Magnetic anomalies due to prism-shaped bodies with arbitrary polarization[J]. Geophysics, 29(4): 517531. [本文引用:1]
[12] Carmichael R S. 1977. Depth calculation of piezomagnetic effect for earthquake prediction[J]. Earth and Planetary Science Letters, 36(2): 309316. [本文引用:1]
[13] Currenti G, Del Negro C, Johnston M J S, et al. 2007. Close temporal correspondence between geomagnetic anomalies and earthquakes during the 2002—2003 eruption of Etna volcano[J]. Journal of Geophysical Research: Solid Earth, 112(B9): B09103. [本文引用:1]
[14] Hamano Y. 1983. Experiments on the stress sensitivity of natural remanent magnetization[J]. Journal of Geomagnetism and Geoelectricity, 35(5): 155172. [本文引用:1]
[15] Hao J Q, Hastie L M, Stacey F D. 1982. Theory of the seismomagnetic effect: A reassessment[J]. Physics of the Earth and Planetary Interiors, 28(2): 129140. [本文引用:3]
[16] Hattori K, Takahashi I, Yoshino C, et al. 2004. ULF geomagnetic field measurements in Japan and some recent results associated with Iwateken Nairiku Hokubu earthquake in 1998[J]. Physics and Chemistry of the Earth, 29(4-9): 481494. [本文引用:1]
[17] Jiang Z, Wang M, Wang Y, et al. 2014. GPS constrained coseismic source and slip distribution of the 2013 MW6. 6 Lushan, China, earthquake and its tectonic implications[J]. Geophysical Research Letters, 41(2): 407413. [本文引用:1]
[18] Johnston M J S, Mueller R J. 1987. Seismomagnetic observation during the 8 July 1986 magnitude 5. 9 North Palm Springs earthquake[J]. Science, 237(4819): 12011203. [本文引用:1]
[19] Johnston M J S, Mueller R J, Sasai Y. 1994. Magnetic field observations in the near-field the 28 June 1992 MW7. 3 Land ers, California, earthquake[J]. Bulletin of the Seismological Society of America, 84(3): 792798. [本文引用:2]
[20] Johnston M J S, Sasai Y, Egbert G D, et al. 2006. Seismomagnetic effects from the long-awaited 28 September 2004 M6. 0 Parkfield earthquake[J]. Bulletin of the Seismological Society of America, 96(4B): S206S220. [本文引用:4]
[21] Kean W F, Day R, Fuller M, et al. 1976. The effect of uniaxial compression on the initial susceptibility of rocks as a function of grain size and composition of their constituent titanomagnetite[J]. Journal of Geophysical Research, 81(5): 861872. [本文引用:1]
[22] Martin R J. 1980. Is piezomagnetism influenced by microcracks during cyclic loading?[J]. Journal of Geomagnetism and Geoelectricity, 32(12): 741755. [本文引用:1]
[23] Mueller R J, Johnston M J S. 1990. Seismomagnetic effect generated by the October 18, 1989, ML7. 1 Loma Prieta, California, earthquake[J]. Geophysical Research Letters, 17(8): 12311234. [本文引用:1]
[24] Nagata T. 1969. Basic magnetic properties of rocks under the effects of mechanical stresses[J]. Tectonophysics, 9(2-3): 167195. [本文引用:2]
[25] Nishida Y, Sugisaki Y, Takahashi K, et al. 2004. Tectonomagnetic study in the eastern part of Hokkaido, NE Japan: Discrepancy between observed and calculated results[J]. Earth Planets and Space, 56(11): 10491058. [本文引用:2]
[26] Okubo A, Oshiman N. 2004. Piezomagnetic field associated with a numerical solution of the Mogi model in a non-uniform elastic medium[J]. Geophysical Journal International, 159(2): 509520. [本文引用:2]
[27] Okubo K, Takeuchi N, Utsugi M, et al. 2011. Direct magnetic signals from earthquake rupturing: Iwate-Miyagi earthquake of M7. 2, Japan[J]. Earth and Planetary Science Letters, 305(1-2): 6572. [本文引用:2]
[28] Oshiman N. 1990. Enhancement of tectonomagnetic change due to non-uniform magnetization in the Earth's crust: Two-dimensional case studies[J]. Journal of Geomagnetism and Geoelectricity, 42(5): 607619. [本文引用:1]
[29] Revol J, Day R, Fuller M D. 1977. Magnetic behavior of magnetite and rocks stressed to failure: Relation to earthquake prediction[J]. Earth and Planetary Science Letters, 37(2): 296306. [本文引用:1]
[30] Sasai Y. 1980. Application of the elasticity theory of dislocations to tectonomagnetic modelling[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 55(2): 387447. [本文引用:2]
[31] Sasai Y. 1991. Tectonomagnetic modeling on the bases of the linear piezomagnetic effect[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 66(4): 585722. [本文引用:2]
[32] Sasai Y. 1994. Piezomagnetic fields produced by dislocation sources[J]. Surveys in Geophysics, 15(4): 363382. [本文引用:1]
[33] Sasai Y, Ishikawa Y. 1978. Changes in the geomagnetic total force intensity associated with the anomalous crustal activity in the eastern part of the Izu Peninsula(2): The Izu-Oshima-Kinkai earthquake of 1978[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 53(2): 893923(in Japanese). [本文引用:2]
[34] Sasai Y, Ishikawa Y. 1997. Seismomagnetic models for earthquakes in the eastern part of Izu Peninsula, central Japan[J]. Annals of Geophysics, 40(2): 463478. [本文引用:1]
[35] Stacy F D. 1964. The seismomagnetic effect[J]. Pure and Applied Geophysics, 58(1): 522. [本文引用:1]
[36] Takla E M, Yomoto K, Liu J Y, et al. 2011. Anomalous geomagnetic variations possibly linked with the Taiwan earthquake( MW=6. 4)on 19 December 2009[J]. International Journal of Geophysics, 211: 110. [本文引用:1]
[37] Thebault E, Finlay C C, Alken P, et al. 2015a. Evaluation of cand idate geomagnetic field models for IGRF -12[J]. Earth, Planets and Space, 67(1): 787804. [本文引用:1]
[38] Thebault E, Finlay C C, Beggan C D, et al. 2015b. International geomagnetic reference field: The 12th generation[J]. Earth, Planets and Space, 67(1): 7998. [本文引用:1]
[39] Utada H, Shimizu H, Ogawa T, et al. 2011. Geomagnetic field changes in response to the 2011 off the Pacific coast of Tohoku earthquake and tsunami[J]. Earth and Planetary Science Letters, 311(1-2): 1127. [本文引用:4]
[40] Wang R, Lorenzo-Martín F, Roth F. 2006. PSGRN/PSCMP: A new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory[J]. Computers and Geosciences, 32(4): 527541. [本文引用:1]
[41] Yamazaki K. 2013. Improved models of the piezomagnetic field for the 2011 MW9. 0 Tohoku-oki earthquake[J]. Earth and Planetary Science Letters, 363: 915. [本文引用:4]