Crustal density variation and its relationship with the seismogenic structure of strong earthquakes in Sichuan-Yunnan region
-
摘要: 为研究川滇地区地壳密度的时空演化特征与强震孕育的关系，本文基于2011—2014年川滇地区的重力复测资料，利用阻尼最小二乘反演算法，获得了川滇地区0—50 km深度范围内分辨率为55 km×55 km×10 km （长×宽×高）的三维动态密度变化模型。以所获取的动态密度变化为依据，分析了川滇地区三维密度变化特征与2013年四川芦山MS7.0、2014年云南鲁甸MS6.5和四川康定MS6.3地震的关系，并由此对强震重点构造部位的深部地壳结构特征、孕震背景及区域动力学过程进行了深入分析。结果显示：川滇地区出现多个与主要活动断裂带展布方向基本一致的密度变化高梯度带，在三次地震的震中区及其附近观测到明显的区域性密度变化异常。15—35 km深度范围内的密度变化水平剖面显示：强震容易发生在上地壳密度变化正、负异常过渡的高梯度带和密度变化四象限分布的中心；中地壳深度密度变化低异常是强震孕育的主要介质条件；下地壳深度密度变化低异常或密度变化高梯度带均有可能是孕育地震的主要介质结构。0—50 km深度的垂直剖面上的密度变化结果显示，地震震中区及附近浅部、深部地壳呈现解耦变化。壳内垂向正负密度变化过渡带可能是强震孕育的又一个主要特征构造。Abstract: In order to study the relationship between the spatio-temporal evolution characteris-tics of crustal density and the seismogenic process of strong earthquakes in Sichuan-Yunnan region, this paper obtained the three-dimensional dynamic density variation model of 0−60 km depth with the resolution of 55 km×55 km×10 km （length, width and height） based on the re-measurement gravity data in the Sichuan-Yunnan region during the period 2011−2014 by the damped least squares inversion algorithm. Then we analyzed the relationship between the three-dimensional dynamic density structure characteristics and the seismogenic characteristics of the LushanMS7.0 in 2013, the Ludian MS6.5 in 2014 and the Kangding MS6.3 in 2014. Furthermore we studied the deep crustal structure, seismogenic background, and regional dynamics of the strong earthquakes. The results show that density variation model exhibited multiple density-change high-gradient zones, which is consistent with the extending direction of the main active fault zones in the Sichuan-Yunnan region. There were obvious regional density variation anomalies in the epicentral area of the three earthquakes mentioned above. The density change horizontal profiles in the depth of 15−35 km showed that strong earthquakes were likely to occur in the center of the high-gradient transition of the positive and negative density variation as well as the center of the four-quadrant distribution of the density variation in the upper-crust, and low density-change anomalies in the mid-crust is required for strong earthquakes. In the lower-crust low anomalies and the high gradient zone of the density change may be the main medium structure for the strong earthquakes. The density variation along the vertical section in the depth of 0−50 km shows that decoupling changes appeared in the shallow crust and deep crust of the epicenters and their adjacent areas. The transitional zone between positive and negative density variation in vertical direction in the crust maybe the main structure for occurrence of strong earthquakes.
-
引言
地震孕育及发生过程中伴随着不同程度的电磁辐射异常（Molchanov et al，1992；Karakelian et al，2002；Matsushima et al，2002；汤吉等，2008，2010），天然地震发生前，也常常观测到电磁辐射异常信号（Han et al，2015；Wang et al，2018；Chang et al，2017；Nagao，2002；Uyeda et al，2002；Hattori，2004）。依据信号波长与传播距离的关系，周期介于几百秒到几百赫兹的信号不满足电磁波要求，地震研究人员将其称为地震电磁扰动（姚休义，冯志生，2018）。
在地震电磁扰动研究中，通过有效的解析方法能够从较强的电磁干扰背景中提取出相对较弱的地震电磁扰动信号。近年来，应用较为广泛的方法有极化法、主成分分析法、分形分析法和梯度定向法等。随着区域台站的增加和仪器采样精度的提高，各种解析方法均积累了很多成功的震例（杜爱民等，2004；Cianchini et al，2012；Hattori et al，2013；潘晖等，2014；李琪等，2015；袁桂平等，2015），其中能反映地震地磁扰动异常的产生、传播和分布特征的方法主要有地磁垂直强度极化法和水平椭圆极化法（Hayakawa et al，2000；Gladychev et al，2001；Molchanov et al，2003；Prattes et al，2008；冯志生等，2010；Li et al，2011）。
Hayakawa等（1996）利用震中距为65 km的日本关岛地磁台的磁通门磁力仪观测资料，分析了1993年8月8日关岛M8.0地震前后地磁场极化值YZH的变化特征，结果显示震前两个月0.01—0.05 Hz （100—20 s）频段的YZH逐渐增大，发震时达到最大，震后逐渐恢复。此后的众多震例研究验证了地震前数天至三个月内地磁场存在高极化值异常现象（Hayakawa et al，2000；Ismaguilov et al，2001；Molchanov et al，2003；Hattori，2004；Prattes et al，2008；李琪等，2015；何畅，冯志生，2017；姚休义，滕云田，2017）。冯志生等（2010）的研究表明，垂直强度极化存在季节变化，提取异常前需予以剔除。上述大部分研究仅涉及震中周围部分台站的观测数据。Hattori （2004）的研究涉及两个台站的数据，最大震中距为118 km；李琪等（2015）利用2009年9月至12月云南西部地磁台阵的连续观测数据，分析了距离2009年11月2日宾川M5.0地震震中143 km范围内的四个地磁台资料，由于台站少，空间特征研究仅限于异常幅值的线性拟合分析。许多研究甚至仅采用单台数据进行异常分析，无法给出异常的空间分布特征。
2015年前后，中国地震局国家地磁台网中心新增了大量磁通门观测仪器，不仅大幅度提升了南北地震带地磁场监测的密度，填补了青藏高原及新疆部分地区地磁监测的空白，而且将采样频率从1 min提高至1 s，使得整个中国大陆西部强震多发地区的地磁场观测水平得到了大幅提升。上述磁通门秒采样观测资料为开展不高于1 Hz的地磁场信号分析工作提供了数据支持。
中国大陆的强震主要分布于西部，对西部地区开展极化异常时空演化特征分析具有十分重要的现实意义。鉴于此，本文拟利用极化法对110°E以西区域的磁通门秒采样观测数据进行处理分析，获取震前地震磁扰动异常信号，并分析其时空演化特征，旨在为该地区后续震情跟踪分析提供参考。
1. 数据及分析方法
1.1 地磁秒采样观测数据
中国地震局国家地磁台网中心总计约150个地磁观测站，观测仪器主要有磁通门磁力仪、质子矢量磁力仪以及Overhauser磁力仪，大部分台站架设了两台以上的观测仪器，以保证数据记录的完整性。其中，磁通门磁力仪产出的秒采样数据可用于地磁超低频（ultra-low frequency，缩写为ULF）异常分析。2011年，国家地磁台网中心仅有磁通门磁力仪112套；2015年前后，随着大规模的台网建设，磁通门磁力仪增至191套。截至2019年，国家地磁台网中心配备磁通门磁力仪189套（图1）。
![]()
图 1 2011年至2019年国家地磁台网中心磁通门磁力仪数量Figure 1. The number of fluxgate magnetometers in National Geomagnetic Networks Center of China from 2011 to 2019目前，中国地震局国家地磁台网中心在110°E以西区域的76个台点共安装120套磁通门磁力仪。为了研究大范围地磁场ULF信号特征，本文采用台站分布相对合理（图2）、观测资料较为丰富的2015—2019年数据进行分析。
![]()
图 2 2019年中国西部磁通门磁力仪（秒采样）分布Figure 2. Distribution of fluxgate magnetometers （second sampling） in western Chinese mainland in 20191.2 地震事件
根据中国地震台网地震目录，2015年至2020年1月，中国大陆西部共发生M6.0以上强震11次。除此之外，本文将中国大陆边界附近发生的4次M6.0以上强震也纳入分析范围，但不包括其强余震（表1）。
表 1 2015年至2020年1月中国大陆西部及周边M≥6.0强震信息Table 1. Information of strong earthquakes （M≥6.0） during the period from 2015 to January of 2020 in western Chinese mainland and its surrounding areas发震时刻
年-月-日 时:分:秒M 北纬/° 东经/° 震源深度/km 参考位置 2020-01-19 21:27:55 6.4 39.83 77.21 16 新疆伽师县 2019-11-21 07:50:45 6.0 19.50 101.35 10 老挝 2019-06-17 22:55:43 6.0 28.34 104.9 16 四川长宁县 2019-04-24 04:15:48 6.3 28.40 94.61 10 西藏墨脱县 2017-11-18 06:34:19 6.9 29.75 95.02 10 西藏米林县 2017-08-09 07:27:52 6.6 44.27 82.89 11 新疆精河县 2017-08-08 21:19:46 7.0 33.2 103.82 20 四川九寨沟县 2016-12-08 13:15:03 6.2 43.83 86.35 6 新疆呼图壁县 2016-11-25 22:24:30 6.7 39.27 74.04 10 新疆阿克陶县 2016-10-17 15:14:49 6.2 32.81 94.93 9 青海杂多县 2016-06-26 19:17:11 6.7 39.43 73.4 10 吉尔吉斯斯坦 2016-01-21 01:13:13 6.4 37.68 101.62 10 青海门源县 2015-12-07 15:50:03 7.4 38.2 72.9 30 塔吉克斯坦 2015-07-03 09:07:46 6.5 37.6 78.2 10 新疆皮山县 2015-04-25 14:11:26 8.1 28.2 84.7 20 尼泊尔 1.3 极化分析
提取地震磁扰动信号的关键是突出源于震源区的异常信号，压制源于外源场的信号。Hayakawa等（1996）由此发展出了地磁垂直分量与水平分量幅值比分析方法：
$$ {Y}_{ZH}{\text{＝}}\left|\frac{Z{\text{（}}\omega {\text{）}}}{H{\text{（}}\omega {\text{）}}}\right|{\text{，}}H{\text{（}}\omega {\text{）}}{\text{＝}}\sqrt{{H}_{X}^{2}{(\omega){\text{＋}}H}_{Y}^{2}{\text{（}}\omega {\text{）}}} {\text{，}}$$ (1) 式中，Z（ω）为地磁垂直分量的傅里叶谱值，H（ω）为地磁水平分量全矢量的傅里叶谱值，HX（ω）为地磁水平分量南北向傅里叶谱值，HY（ω）为地磁水平分量东西向傅里叶谱值，地磁各分量单位为nT，ω为圆频率。由于该方法反映了磁场强度或能量在垂直方向的分配比例，将其简称为地磁垂直强度极化法，将地磁垂直分量与水平分量的幅值比YZH称为地磁垂直强度极化值。
1.4 信号处理
根据趋肤效应，可以推算源于震源处的信号频率范围。趋肤深度定义为
$$ \delta {\text{＝}}\sqrt{\frac{2\rho }{\mu \omega }} {\text{，}}$$ (2) 式中δ为趋肤深度，ρ为电阻率，μ为磁导率，μ＝μ0＝4π×10−7 H/m，μ0为真空中的磁导率，ω为角频率。由于不能准确地测得震源处的电阻率，我们用沉积岩的电阻率ρ＝100 Ω·m代替。根据表2所示的2015年以来中国大陆西部YZH高值异常信息，本研究的地震震源深度处于6—30 km范围内。将上述参数代入式（2），计算得到ω的范围为4.42—0.18 Hz。考虑到资料采样率的局限性以及电阻率的不确定性，并综合前人的研究成果（冯志生等，2010；李琪等，2015），我们将频带范围限定到0.01—0.2 Hz。数据处理的步骤如下：① 将磁通门三分量Z，HX和HY每天的秒资料分为96段，每段15 min，分段计算垂直分量和水平分量的傅里叶谱值；② 计算每天5—100 s （0.01—0.2 Hz）周期内各频点的极化值均值，获得极化值的逐日变化序列；③ 对极化值逐日变化序列进行傅里叶拟合（数据一般不短于半年），获得周期大于半年的傅里叶拟合变化曲线，即为极化值的年变化曲线；④ 将5—100 s周期内各频点极化值减去年变化后，剔除低于二倍均方差的部分，仅保留高值部分；⑤ 对筛选出的极化高值的日均值逐日变化序列再次进行傅里叶拟合，如果仍然存在年变化成分，则再次扣除极化值高值周期大于半年的成分，所得残差为最终处理得到的极化高值序列，即YZH逐日值。
表 2 2015年至2019年中国大陆西部$Y_{ZH2}$ 高值异常信息Table 2. Information of high YZH2 value anomalies in western Chinese mainland from 2015 to 2019序号 异常出现时间
年-月-日异常持续时间/d 异常台站最多日（异常日）
年-月-日异常台站数 总台数 YZH2 1 2015-01-16 6 2015-01-18 28 36 0.73 2 2015-02-09 9 2015-02-13 14 36 0.39 3 2015-07-17 5 2015-07-18 8 36 0.29 4 2015-09-26 5 2015-09-29 13 27 0.47 5 2015-11-18 11 2015-11-23 21 27 0.92 6 2016-01-27 8 2016-01-30 17 41 0.75 7 2016-08-14 9 2016-08-17 15 33 0.92 8 2016-09-11 9 2016-09-15 14 33 0.88 9 2016-12-01 6 2016-12-02 14 37 0.75 10 2017-01-13 5 2017-01-15 8 38 0.41 11 2017-02-13 5 2017-02-14 10 37 0.49 12 2017-03-11 9 2017-03-18 14 34 0.66 13 2017-08-24 6 2017-08-28 22 40 0.65 14 2017-10-29 6 2017-10-30 14 40 0.70 15 2018-04-16 3 2018-04-17 11 43 0.62 16 2018-06-10 6 2018-06-12 17 43 0.44 17 2019-09-22 3 2019-09-24 13 54 1.15 18 2019-11-12 5 2019-11-12 26 56 0.85 以往震例（姚休义，冯志生，2018）表明，震前的超低频异常往往持续数天，且不同台站出现时间可能略有差异。虽然这种差异并不会影响到时间演化过程，但在研究异常空间特征时，可能会因为前后出现时间的差异而导致异常的空间范围与实际不相符，因此对逐日值进行5日滑动平均处理，将高值异常进行时间上的延拓，得到YZH1。
为便于异常空间特征的分析，利用各台YZH1的二倍均方差对其进行归一化，对归一化后的残差曲线减1，得到YZH2，该值被称为归一置零极化值。具体算法为：
$$ {Y}_{ZH2}{\text{＝}} \frac{{Y}_{ZH1}}{2\sqrt{{\displaystyle\sum \limits_{j {\text{＝}} 1}^n \Bigg[ {{{{\text{（}} {{Y_{ZH1}}} {\text{）}}}_j} {\text{－}} {{ {\overline {{{{\text{（}} {{Y_{ZH1}}} {\text{）}}}_1} {{{\text{（}} {{Y_{ZH1}}} {\text{）}}}_2} \cdots {{{\text{（}} {{Y_{ZH1}}} {\text{）}}}_j}\cdots {{{\text{（}} {{Y_{ZH1}}} {\text{）}}}_n}} } }}} \Bigg]^2}}}{\text{－}}1 {\text{.}}$$ (3) YZH2高于0，等价于YZH1高于2倍方差。图3中展示了喀什台2017年YZH1与YZH2的变化曲线，可以看到两者具有相同的变化趋势。 利用现有台站的归一置零极化值，在空间上进行简单的数学插值（克里金法），即可得到大陆西部地区每日的归一置零极化值的空间分布。
![]()
图 3 喀什台2017年YZH1与YZH2的对比Figure 3. Comparison of YZH1 with YZH2 results of Kashi station in 20172. 结果分析
2.1 极化法计算结果
对2015年以来的76个台点120套观测仪器所产出的数据进行上述处理后，每套观测资料计算得到一个YZH1序列。因观测环境变化或仪器故障等原因，每年能够正常产出数据的台站数出现小幅波动。
YZH序列与YZH1序列在中国大陆西部各台具有较好的同步性。以2017年为例（图4），由南向北分布的通海台（THA）、西昌台（XIC）、江油台（JOY）、天水台（TYS），以及由西向东分布的喀什台（KSH）、且末台（QMO）、都兰台（DUL）、英鸽台（YGE），其YZH与YZH1序列均具有相同的变化趋势，但变化幅度未呈现随经度或纬度规律变化的特征。
![]()
图 4 沿经线和纬线分布的部分台站及其YZH和YZH1 2017年结果Figure 4. Location of some stations and their YZH and YZH1 results in 2017基于以上规律，对YZH1序列进一步处理，并再一次筛选参与分析的台站。因原始数据存在一些干扰、突跳甚至错误，YZH1时序曲线出现显著偏离背景的情况，会出现与周边台站出现显著的不一致现象，因此在确认原始数据出错的情况下，删除显著偏离背景时段的YZH1结果。对于有两套以上仪器的台站，仅选择数据连续性更好、与周边台站同步性更高的一套参与后续的空间特征分析。
2.2 极化高值与地磁扰动
为了研究极化值与外源场扰动之间的关系，我们在东京地磁国际数据中心网站http://wdc.kugi.kyoto-u.ac.jp下载了2015年至2019年间的地磁Kp指数和Dst指数，并将两者与YZH （逐日值）进行相关性分析。计算结果显示每天的Kp指数最大值与YZH的相关系数为−0.449，表明Kp指数与YZH存在弱的负相关，即Kp指数高值与YZH低值可能存在弱的相关性，这与冯志生等（2010）的研究结果一致，即“当K指数为高值时，极化结果为低值”。这一结果同时也符合该方法的预期，即来自电离层外空磁场一次源的垂直分量幅值与水平分量幅值的比值较小。此外，Dst指数最大值与YZH的相关系数为0.176，表明两者不相关。
从2017年Kp指数、Dst指数与喀什台YZH1结果（图5）来看，极化高值往往出现在Kp指数小于3的时段和Dst指数处于±20之间的时段，即地磁扰动较弱的时段。
![]()
图 5 喀什台2017年YZH1序列与Kp （a）和Dst （b）指数曲线的对比Figure 5. Comparison of YZH1 series with Kp （a） and Dst （b） curves in Kashgar station in 20172.3 中国大陆西部超低频高值异常
YZH1单台高值现象较多，7个及以下台站出现高值时，高值台站也往往孤立分布，其可靠性有待验证。因此本文仅讨论2015年至2019年中国大陆西部出现的8个以上台点同步YZH1
高值现象。 表2给出了2015年至2019年中国大陆西部共出现的18次8个以上台点同步YZH2高值现象，图6对每月出现高值异常的频次进行了统计，结果表明：异常出现频次最高的月份为1月和9月，均为3次；其次是2月、8月、11月，分别出现2次；除5月外，其余各月均有1次高值异常。总体呈春夏低秋冬高的趋势，但并未呈明显的优势分布。
![]()
图 6 YZH1高值异常在12个月中的频次分布 Figure 6. The frequency distribution of YZH1 high value anomalies in 12 months图7给出了2015年至2019年期间YZH1高值频次与地震频次的对比，可见在YZH1高值频次较低的2018年，后续无M6.0以上强震发生。
![]()
图 7 YZH1高值异常频次与地震频次图Figure 7. The frequency of YZH1 high value anomalies and earthquake frequency为便于后续描述，我们按时间顺序分别将18次高值现象编号为1—18，并取高值台站最多日作为研究日，对上述高值异常的空间分布特征进行研究。
图8给出了2015年至2019年期间部分台站的YZH1时序图，图9和10给出了18次高值异常的空间分布图。从空间等值线图不难看出，每次高值现象出现时，高值区域的位置和面积均不相同：研究区域东部的高值区分布特征较为复杂，总体来说高值区域面积较小；西部地区由于台站稀疏，更容易出现大面积高值区，这种大面积高值可能与插值有关，但也从一定程度上反映了相关区域实际YZH2偏高的情况。从时间变化过程来看，2015年1月、11月和2016年9月高值区面积较大（图9），2018年的两次高值区面积较小（图10），其余高值区的面积相差不明显。
![]()
![]()
图 9 2015年至2016年磁场高值日YZH2的空间分布图Figure 9. Spatial distribution of YZH2 on the magnetic field high value day from 2015 to 2016![]()
图 10 2017年至2019年磁场高值日YZH2的空间分布图 Figure 10. Spatial distribution of YZH2 on the magnetic field high value day from 2017 to 2019为了进一步验证高值现象与强震相关的可信度，本文选择中国东部的静海台、高邮台与西部的喀什台进行对比。以2017年为例（图11），喀什台的YZH1随时间变化的趋势与高邮台大体相同，但不同时段的相对变化幅值差异较大。图中红色虚线标示出YZH1出现高值异常的时间，喀什台YZH1序列高值时间与其基本吻合，而东部的高邮台YZH1序列在上述时间并无明显的相对高值出现。计算上述两台站的YZH1序列相关系数为0.26，表明两台站资料存在弱相关性；而静海台与喀什台的曲线差别较大，YZH1序列相关系数为0.05，表明两者不相关。
![]()
图 11 喀什台、静海台和高邮台的空间分布（a）以及2017年极化YZH1时序曲线对比图（b）图中红色竖线标示中国大陆西部YZH1出现高值异常的时间Figure 11. The locations of the stations Kashi，Jinghai and Gaoyou （a） and YZH1 time series curves in 2017 （b）The red vertical line indicates the time when the YZH1 high value anomaly appeared in western Chinese mainland2.4 强震与YZH高值的时空关系
2.4.1 时空关系
2015—2019年中国大陆西部共出现18次8个以上台站YZH1同步高值现象，共发生15次M6.0以上强震；2015—2017年内共出现14次多台同步高值现象，发生11次M6.0以上强震；2018—2019年共出现4次高值现象，较2015—2017年频次明显降低，强震频次也有所降低，共发生4次强震。YZH高值频次与地震频次的相关系数为0.47，即二者存在弱相关。
对2015年1月至2020年1月中国大陆西部及周边15次强震发生前半年的YZH1高值现象进行回溯性研究，结果显示，除2019年4月24日墨脱M6.3地震和2019年6月17日长宁M6.0地震以外，其余13次强震发生前半年内，YZH1均出现多台同步高值。上述13次强震中，除2016年11月25日新疆阿克陶M6.7地震外，其余12次强震震中及周边均存在高值区。
值得注意的是，在2015年7月3日皮山M6.5地震后15天、2016年1月21日门源M6.4地震后9天、2016年12月8日呼图壁M6.2地震后37天、2017年8月9日精河M6.6地震后19天以及2017年8月8日九寨沟M7.0地震后20天内，出现了YZH1多台同步高值现象且地震震中位于高值区（表3）。
表 3 单次高值异常及短期内多次高值异常叠加后发生强震的高值区面积Table 3. The high value areas of single high-value anomalies and the superposition of several high-value anomalies in a short period which followed with strong earthquakes序号 出现高值日期
年-月-日震后出现高值区的
地震（提前的天数）半年内高值区强震 高值区面积①
/（104 km2）高值区总面积②
/（104 km2）1 2015-01-18 无 2015-07-03皮山M6.5 251.6 313 2 2015-02-13 2015-04-25尼泊尔M8.1 3 2015-07-18 2015-07-03皮山M6.5
（15）2015-12-07塔吉克斯坦M7.4 13 20.6 4 2015-09-29 无 2016-01-21门源M6.4 20 144 5 2015-11-23 2015-12-07塔吉克斯坦M7.4 123.7 6 2016-01-30 2016-01-21门源M6.4 （9） 2016-06-26吉尔吉斯斯坦M6.7 17.5 84.7 7 2016-08-17 无 2016-12-08呼图壁M6.2 24.1 82.2 8 2016-09-15 2016-10-17杂多M6.2 58 9 2016-12-02 无 2016-12-08呼图壁M6.2 18.2 80.8 10 2017-01-13 2016-12-08呼图壁M6.2 （37） 无 47.2 47.2 11 2017-02-14 无 2017-08-08九寨沟M7.0 41.6 84.4 12 2017-03-18 2017-08-09精河M6.6 30.5 13 2017-08-28 2017-08-08九寨沟M7.0 （19）
2017-08-09精河M6.6 （20）无 39.5 14 2017-10-31 无 16.1 15 2018-04-17 无 无 32.4 16 2018-06-12 17 2019-09-24 无 2019-11-21老挝M6.0 19.2 72.3 18 2019-11-12 注：① 指地震震中所在高值区的面积；② 指一次高值现象出现的多个高值区面积之和。 2.4.2 短期内反复出现高值的区域易发强震
上述研究还存在一个现象，YZH2多台高值现象可能在短期内多次出现，且后续强震的震中位置往往为多次出现高值的区域。为了进一步分析这种现象，我们将短期（两个月）内出现且高值区有重叠的多次YZH2高值进行简单叠加，具体做法为：对同一点的N次多台高值日YZH2求平均，经过上述处理，多次均为高值的区域仍为大于0的高值区。但为了排除可能的震后效应产生的影响，强震发生后震中附近区域出现的高值不叠加。
经过筛选，短期内（两个月）出现且高值区有重叠的YZH2多台同步高值现象共有6组，序号分别为①-②，④-⑤，⑦-⑧，⑩-⑪，⑮-⑯，⑰-⑱，叠加结果如图12所示。叠加过程使得多次高值的区域突出，高值区面积缩小。不难看出，后续强震震中大部分位于多次叠加后的高值区域及其边缘，2016年1月25日阿克陶M6.7地震和2019年11月21日伽师M6.0地震除外。其中阿克陶地震前出现了两次高值现象，但在震中附近均未出现高值区，而伽师地震前58天在震中出现了高值区。
![]()
图 12 2015—2019年6组高值异常叠加前、后的YZH2与后续6个月内地震的空间关系Figure 12. Spatial relationship of six groups of high YZH2 anomalies before and after superposition as well as the distribution of earthquakes in the following six months from 2015 to 20192.4.3 一次或一组多台高值现象中多个高值区发生强震
图12显示，一次高值现象中会出现多个YZH2高值区域。多次高值现象累加后，也产生多个高值区域。从上述研究中我们不难发现，一次多台高值现象中的高值区或一组高值现象叠加后出现的多个高值区均有可能发生强震，一个高值区内也可能发生多次强震。这种现象在叠加后的结果中表现得更明显，例如：①与②的叠加结果中，尼泊尔地震和皮山地震位于同一个高值区内；④与⑤的叠加结果中，塔吉克斯坦地震和门源地震的震中分别位于两个高值区内；⑦与⑧的叠加结果中，呼图壁地震和杂多地震的震中分别位于两个高值区内；⑩与⑪的叠加结果中，九寨沟地震和精河地震的震中分别位于两个高值区内。
![]()
图 12 2015—2019年6组高值异常叠加前、后的YZH2与后续6个月内地震的空间关系Figure 12. Spatial relationship of six groups of high YZH2 anomalies before and after superposition as well as the distribution of earthquakes in the following six months from 2015 to 20192.4.4 YZH2高值区面积与后续强震震级的关系
我们基于叠加结果，将18次高值异常重新分组，产生12组高值，每组均有一幅YZH2空间分布图。统计单张YZH2空间分布图中高值区的总面积及后续半年内发生强震的高值区域面积（表3），对上述两个参量与后续地震震级的关系进行了分析，结果如图13所示。可以看出，高值区面积有随震级增大而增大的趋势，相关系数为0.89。
![]()
图 13 叠加后的YZH2高值区面积与后续6个月内地震震级的关系Figure 13. The relationship between YZH2 high-value areas of the superposition of several high-value anomalies and the magnitude of strong earthquakes in six months3. 讨论与结论
本文基于2015至2019年期间中国大陆西部磁通门秒采样观测资料，开展了5—100 s频段地磁垂直强度的极化分析。分析结果表明，中国大陆西部各台站的极化结果具有较好的一致性，沿经向和纬向均无明显的形态及幅值变化，往往出现多台同步高值现象。虽然出现同步高值的台站较多、空间面积较大，但对极化结果与Kp指数及Dst指数的互相关系数分析结果表明，极化高值与外源场扰动的关系不大。对YZH2高值与后续强震的时空关系研究显示：YZH2高值现象出现后的半年内，高值区域及周边容易发生强震；有部分地震后短时间内，震中周边也出现了极化高值；高值现象出现后，多个高值区均有可能发生强震；后续强震的震级与高值区面积呈正相关。值得注意的是，上述结果仅仅是基于5年的数据所得，因统计样本有限，可能不能完全反映客观规律。随着资料的进一步积累，将进一步加深对地磁垂直强度极化高值异常与地震关系的认识。
尽管地磁垂直强度极化法经过了很多震例的验证，但仍有研究人员质疑该方法提取震磁异常信息的可靠性（Thomas et al，2009；Masci，2011），主要存在二点疑问：第一，因磁层电流体系的影响无法排除，认为极化值升高很可能与外源场变化有关；第二，以往研究中震前异常通常仅存在于震中距最小的单一台站，在无法排除数据干扰的情况下，异常的可靠性也常被质疑（姚休义，冯志生，2018）。
关于第一点，本文中进行了简单的讨论，认为极化高值与外源场扰动的关系不大。冯志生等（2010）也曾就此问题展开研究，其结果表明地磁垂直强度极化高值与地磁活动性指数之间存在较好的负相关性。Currie和Waters （2014） 对各频段极化值与Kp指数、AE指数及SYM-H指数的相关性分析结果表明，各频段极化值与上述各种地磁活动性指数之间的相关系数均低于0.3，即两者并不相关。
关于第二点，本文开展了大区域多台站极化高值与地震相关性的研究，结果显示高值往往在较大范围内的多个台站同时出现，与台站附近的干扰无关。李琪等（2015）利用极化法在多个台站均提取到了震前极化异常，且异常特征符合地震电磁扰动信号的衰减特征（Huang，Ikeya，1998）。上述研究均表明极化值升高的结果真实可信。
但与以往研究结果不同的是，本研究显示震前出现高值的空间范围不局限于数十千米，时间也不局限于数日之内，有可能达到半年左右。梅世蓉（1996）在华北地区开展的研究表明，一次M6.0以上地震前常出现持续时间长达数十年的椭圆形活动增强区，其长度至少长达160 km，M8.0以上强震前甚至能够达到1 000 km以上。而活动增强区的微震破裂很可能产生电磁辐射并出现极化高值。
同时出现的多个异常区有可能发生多次强震。根据马宗晋（1980）的研究结果，震前区域应力会在多点集中，跨区域的多次强震之间很可能存在联系而出现强震连发的情况。如果这种联系是存在的，那么在某种条件的激发下，就有可能在多个区域同时出现极化高值。而同时出现的多个极化高值区后续是否发震，取决于每个孕震体的发育情况。震级与高值区面积呈正相关的结果也支持上述推论。
此外，震后一个月左右震中附近有可能出现极化值再次升高的情况。关于这一现象的机理尚需要进一步研究。
在数据下载及处理过程中，上海市地震局朱培育高级工程师提供了软件支持；文中插图主要使用GMT （Wessel，Smith，1995）绘制，作者在此表示衷心的感谢。
-
![]()
图 1 研究区域地质构造及强震分布
F1：龙门山断裂带；F2：鲜水河断裂带；F3：小江断裂带；F4：金沙江断裂带；F5：丽江—小金河断裂带；F6：楚雄—建水断裂带
Figure 1. The regional geological tectonics and earthquake distribution in the studied area
F1：Longmenshan fault zone；F2：Xianshuihe fault zone；F3：Xiaojiang fault zone；F4：Jinshajiang fault zone；F5：Lijiang-Xiaojinhe fault zone；F6：Chuxiong-Jianshui fault zone
![]()
图 2 2011—2014年的动态重力变化分布图
（a） 2011—2012年；（b） 2012—2013年；（c） 2013—2014年
Figure 2. Distribution of dynamic gravity variation from 2011 to 2014
（a） 2011−2012；（b） 2012−2013；（c） 2013−2014
![]()
图 3 15，25和35 km深度下检测板分辨率测试结果，格网剖分单元为55 km×55 km×10 km
（a） 2011—2012年；（b） 2012—2013年；（c） 2013—2014年
Figure 3. The results of the checkboard resolution test with the depths of 15，25，and 35 km， where 55 km×55 km×10 km combination abnormal block was adopted
（a） 2011−2012；（b） 2012−2013；（c） 2013−2014
![]()
图 4 基于2011—2012 （左）、2012—2013 （中）、2013—2014年 （右）重力数据的阻尼最小二乘反演的关键参数及迭代收敛过程
（a） 重力变化残差方差的对数
${{\rm{ lg}}\ \left\|\ { d}-{ {G\rho}}\ \right\|}$ （横轴）和密度变化模型估计误差方差的对数${{\rm{ lg}}\ \left\|\ { \rho}\ \right\|}$ （纵轴）的折中曲线；（b） 重力变化数据拟合残差随迭代次数的变化；（c） 最终密度结果正演所得异常值与实测之差的分布直方图Figure 4. Key parameters and iterative process of the damped least squares inversion based on the gravity variation from 2011 to 2012 （left），from 2012 to 2013 （middle） and，from 2013 to 2014 （right）
（a） The tradeoff curve combined by the logarithmic
${{\rm{ lg}}\left\|{ d}-{ {G\rho}}\right\| }$ of gravity changes residual variances （horizontal axis） and the logarithmic$ {{\rm{ lg}}\left\|{ \rho}\right\|}$ of density changes model variances （vertical axis）；（b） Variation of the gravity changes residual variances with the number of iterations；（c） Histogram of the differences between observed gravity changes and the calculated changes from forward modeling![]()
图 5 15，25和35 km深度下川滇地区地壳2011—2014年三维动态密度变化分布图
（a） 2011—2012年；（b） 2012—2013年；（c） 2013—2014年
Figure 5. Distribution of 3-D dynamic crustal density variation in Sichuan-Yunnan region from 2011 to 2014 at the depths of 15，25，and 35 km
（a） 2011−2012；（b） 2012−2013；（c） 2013−2014
![]()
图 6 2011—2014年川滇地区沿27.25°N （上）和32.25°N （下）垂直剖面的地壳密度变化图
（a） 2011—2012年；（b） 2012—2013年；（c） 2013—2014年
Figure 6. Crustal density change profiles along the latitude 27.25°N （upper panels） and 32.25°N（lower panels） in Sichuan-Yunnan region from 2011 to 2014
（a） 2011−2012；（b） 2012−2013；（c） 2013−2014
-
陈石,王青华,王谦身,王岩,卢红艳,徐伟民,石磊,郭凤义. 2014. 云南鲁甸MS6.5地震震源区和周边三维密度结构及重力场变化[J]. 地球物理学报,57(9):3080–3090. doi: 10.6038/cjg20140933 Chen S,Wang Q H,Wang Q S,Wang Y,Lu H Y,Xu W M,Shi L,Guo F Y. 2014. The 3D density structure and gravity change of Ludian MS6.5 Yunnan epicenter and surrounding region[J]. Chinese Journal of Geophysics,57(9):3080–3090 (in Chinese). doi: 10.6038/cjg20140933（inChinese）
郝明. 2012. 基于精密水准数据的青藏高原东缘现今地壳垂直运动与典型地震同震及震后垂直形变研究[D]. 北京: 中国地震局地质研究所: 56−57. Hao M. 2012. Present Crustal Vertical Movement of Eastern Tibetan Plateau and Coseismic and Postseismic Vertical Deformation of Two Typical Earthquakes[D]. Beijing: Institute of Geology, China Earthquake Administration: 56−57 （in Chinese）.
侯强,邹文远,欧明霖,丁小军. 2018. 青藏高原东南缘壳幔力学耦合及其动力学意义[J]. 大地测量与地球动力学,38(10):991–1000. Hou Q,Zou W Y,Ou M L,Ding X J. 2018. Mechanical coupling in southeastern Qinghai-Tibetan Plateau and its geodynamic implications[J]. Journal of Geodesy and Geodynamics,38(10):991–1000 (in Chinese).
李红蕾,方剑,王新胜,刘杰,崔荣花,陈铭. 2017. 重力及重力梯度联合反演青藏高原及邻区岩石圈三维密度结构[J]. 地球物理学报,60(6):2469–2479. doi: 10.6038/cjg20170634 Li H L,Fang J,Wang X S,Liu J,Cui R H,Chen M. 2017. Lithospheric 3-D density structure beneath the Tibetan Plateau and adjacent areas derived from joint inversion of satellite gravity and gravity-gradient data[J]. Chinese Journal of Geophysics,60(6):2469–2479 (in Chinese). doi: 10.6038/cjg20170634（inChinese）
梁伟锋,刘芳,徐云马,祝意青,郭树松,赵云峰,刘练. 2013. 青藏高原东缘重力观测及对芦山M7.0地震的反映[J]. 地震工程学报,35(2):266–271. doi: 10.3969/j.issn.1000-0844.2013.02.0266 Liang W F,Liu F,Xu Y M,Zhu Y Q,Guo S S,Zhao Y F,Liu L. 2013. Gravity observation in the eastern margin of Qinghai-Tibetan Plateau and reflect to the Lushan M7.0 earthquake[J]. China Earthquake Engineering Journal,35(2):266–271 (in Chinese).
楼海,王椿镛. 2005. 川滇地区重力异常的小波分解与解释[J]. 地震学报,27(5):515–523. doi: 10.3321/j.issn:0253-3782.2005.05.006 Lou H,Wang C Y. 2005. Wavelet analysis and interpretation of gravity data in Sichuan-Yunnan region,China[J]. Acta Seismologica Sinica,27(5):515–523 (in Chinese).
孟小红,石磊,郭良辉,佟拓,张盛. 2012. 青藏高原东北缘重力异常多尺度横向构造分析[J]. 地球物理学报,55(12):3933–3941. doi: 10.6038/j.issn.0001-5733.2012.12.006 Meng X H,Shi L,Guo L H,Tong T,Zhang S. 2012. Multi-scale analyses of transverse structures based on gravity anomalies in the northeastern margin of the Tibetan Plateau[J]. Chinese Journal of Geophysics,55(12):3933–3941 (in Chinese). doi: 10.6038/j.issn.0001-5733.2012.12.006（inChinese）
申重阳. 2005. 地壳形变与密度变化耦合运动探析[J]. 大地测量与地球动力学,25(3):7–12. Shen C Y. 2005. Preliminary analysis of coupling movement between crustal deformation and density change[J]. Journal of Geodesy and Geodynamics,25(3):7–12 (in Chinese).
申重阳,李辉,孙少安,刘少明,玄松柏,谈洪波. 2009. 重力场动态变化与汶川MS8.0地震孕育过程[J]. 地球物理学报,52(10):2547–2557. doi: 10.3969/j.issn.0001-5733.2009.10.013 Shen C Y,Li H,Sun S A,Liu S M,Xuan S B,Tan H B. 2009. Dynamic variations of gravity and the preparation process of the Wenchuan MS8.0 earthquake[J]. Chinese Journal of Geophysics,52(10):2547–2557 (in Chinese).
王椿镛,李永华,楼海. 2016. 与青藏高原东北部地球动力学相关的深部构造问题[J]. 科学通报,61(20):2239–2263. Wang C Y,Li Y H,Lou H. 2016. Issues on crustal and upper-mantle structures associated with geodynamics in the northeastern Tibetan Plateau[J]. Chinese Science Bulletin,61(20):2239–2263 (in Chinese). doi: 10.1360/N972016-00160
王嘉沛,申重阳,玄松柏,杨光亮,郝洪涛. 2018. 2017年四川九寨沟MS7.0地震震前地壳密度变化特征[J]. 大地测量与地球动力学,38(8):771–776. Wang J P,Shen C Y,Xuan S B,Yang G L,Hao H T. 2018. Density variation in the crust before the Jiuzhaigou,Sichuan,MS7.0 earthquake of 2017[J]. Journal of Geodesy and Geodynamics,38(8):771–776 (in Chinese).
王绪本,高原,王志,余嘉顺. 2017. 青藏高原东缘深部地球物理与大陆动力学研究进展[J]. 地球物理学报,60(6):2030–2037. doi: 10.6038/cjg20170601 Wang X B,Gao Y,Wang Z,Yu J S. 2017. Research progress on deep geophysics and continental dynamics in eastern Tibetan Plateau[J]. Chinese Journal of Geophysics,60(6):2030–2037 (in Chinese). doi: 10.6038/cjg20170601（inChinese）
闻学泽,杜方,易桂喜,龙锋,范军,杨攀新,熊仁伟,刘晓霞,刘琦. 2013. 川滇交界东段昭通、莲峰断裂带的地震危险背景[J]. 地球物理学报,56(10):3361–3372. doi: 10.6083/cjg20131012 Wen X Z,Du F,Yi G X,Long F,Fan J,Yang P X,Xiong R W,Liu X X,Liu Q. 2013. Earthquake potential of the Zhaotong and Lianfeng fault zones of the eastern Sichuan-Yunnan border region[J]. Chinese Journal of Geophysics,56(10):3361–3372 (in Chinese). doi: 10.6083/cjg20131012（inChinese）
吴微微,吴朋,魏娅玲,孙玮. 2017. 川滇活动块体中—北部主要活动断裂带现今应力状态的分区特征[J]. 地球物理学报,60(5):1735–1745. doi: 10.6038/cjg20170511 Wu W W,Wu P,Wei Y L,Sun W. 2017. Regional characteristics of stress state of main seismic active faults in mid-northern part of Sichuan-Yunnan block[J]. Chinese Journal of Geophysics,60(5):1735–1745 (in Chinese). doi: 10.6038/cjg20170511（inChinese）
许昭永,文丽敏,石宝文,许峻,樊俊屹,李正媛,苏有锦. 2015. 川滇菱形地块的应变能积累释放周期和强震预测[J]. 地震学报,37(5):774–786. doi: 10.11939/jass.2015.05.006 Xu Z Y,Wen L M,Shi B W,Xu J,Fan J Y,Li Z Y,Su Y J. 2015. The strain energy accumulating-releasing period and strong earthquake prediction in Sichuan-Yunnan rhombic block[J]. Acta Seismologica Sinica,37(5):774–786 (in Chinese). doi: 10.11939/jass.2015.05.006（inChinese）
杨海燕,胡家富,赵宏,苏有锦. 2009. 川西地区壳幔结构与汶川MS8.0级地震的孕震背景[J]. 地球物理学报,52(2):356–364. Yang H Y,Hu J F,Zhao H,Su Y J. 2009. Crust-mantle structure and seismogenic background of Wenchuan MS8.0 earthquake in western Sichuan area[J]. Chinese Journal of Geophysics,52(2):356–364 (in Chinese).
祝意青,梁伟锋,徐云马. 2008. 重力资料对2008年汶川MS 8.0地震的中期预测[J]. 国际地震动态,(7):36–39. doi: 10.3969/j.issn.0253-4975.2008.07.006 Zhu Y Q,Liang W F,Xu Y M. 2008. Medium-term prediction of MS8.0 earthquake in Wenchuan,Sichuan by mobile gravity[J]. Recent Developments in World Seismology,(7):36–39 (in Chinese).
祝意青,徐云马,吕弋培,李铁明. 2009. 龙门山断裂带重力变化与汶川8.0级地震关系研究[J]. 地球物理学报,52(10):2538–2546. doi: 10.3969/j.issn.0001-5733.2009.10.012 Zhu Y Q,Xu Y M,Lü Y P,Li T M. 2009. Relations between gravity variation of Longmenshan fault zone and Wenchuan MS8.0 earthquake[J]. Chinese Journal of Geophysics,52(10):2538–2546 (in Chinese). doi: 10.3969/j.issn.0001-5733.2009.10.012（inChinese）
祝意青,刘芳,李铁明,郑兵,王青华. 2015a. 川滇地区重力场动态变化及其强震危险含义[J]. 地球物理学报,58(11):4187–4196. doi: 10.6038/cjg20151125 Zhu Y Q,Liu F,Li T M,Zheng B,Wang Q H. 2015a. Dynamic variation of the gravity field in the Sichuan-Yunnan region and its implication for seismic risk[J]. Chinese Journal of Geophysics,58(11):4187–4196 (in Chinese). doi: 10.6038/cjg20151125（inChinese）
祝意青,付广裕,梁伟锋,徐云马. 2015b. 鲁甸 MS6.5、芦山MS7.0、汶川MS8.0地震前区域重力场时变[J]. 地震地质,37(1):319–330. doi: 10.3969/j.issn.0253-4967.2015.01.025 Zhu Y Q,Fu G Y,Liang W F,Xu Y M. 2015b. Earthquake predictions: Spatial-temporal gravity changes before the Ludian MS6.5,Lushan MS7.0 and Wenchuan MS8.0 earthquakes[J]. Seismology and Geology,37(1):319–330 (in Chinese).
Backus G E,Gilbert J F. 1967. Numerical applications of a formalism for geophysical inverse problems[J]. Geophys J Int,13(1/2/3):247–276.
Bai D H,Unsworth M J,Meju M A,Ma X B,Teng J W,Kong X R,Sun Y,Sun J,Wang L F,Liang C S,Zhao C P,Xiao P F,Liu M. 2010. Crustal deformation of the eastern Tibetan Plateau revealed by magnetotelluric imaging[J]. Nat Geosci,3(5):358–362. doi: 10.1038/ngeo830
Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve[J]. SIAM Rev,34(4):561–580. doi: 10.1137/1034115
Jiang W L,Zhang J F,Tian T,Wang T X. 2012. Crustal structure of Chuan-Dian region derived from gravity data and its tectonic implications[J]. Phys Earth Planet Inter,212/213:76–87. doi: 10.1016/j.pepi.2012.07.001
Lei J S,Zhao D P. 2016. Teleseismic P-wave tomography and mantle dynamics beneath eastern Tibet[J]. Geochem Geophys Geosyst,17(5):1861–1884. doi: 10.1002/2016GC006262
Li Y G,Oldenburg D W. 1998. 3-D inversion of gravity data[J]. Geophysics,63(1):109–119. doi: 10.1190/1.1444302
Li Y G. 2001. 3-D inversion of gravity gradiometer data[C]//2001 SEG Annual Meeting. San Antonio, Texas: Society of Exploration Geophysicists: 1470−1473.
Merriam D F. 2006. The Conrad discontinuity in the Midcontinent (USA)[J]. Trans Kansas Acad Sci,109(3):125–130.
Zhang G,Wang X B,Fang H,Guo Z M,Zhang Z B,Luo W,Cai X L,Li J,Li Z,Wu X. 2015. Crust and upper mantle electrical resistivity structure in the Panxi region of the eastern Tibetan Plateau and its significance[J]. Acta Geologica Sinica,89(2):531–541. doi: 10.1111/1755-6724.12445
Zhang J,Li W Y,Tang X C,Tian J,Wang Y C,Guo Q,Pang Z H. 2017. Geothermal data analysis at the high-temperature hydrothermal area in western Sichuan[J]. Science China Earth Science,60(8):1507–1521. doi: 10.1007/s11430-016-9053-2
-
期刊类型引用(2)
1. 韩晨，梁锋，李海龙，王涛，张敏，唐灵. 基于背景噪声的自相关与HVSR联合成像在广东惠州城市浅层结构探测中的应用. 地质论评. 2022(04): 1413-1422 . 
百度学术
2. 高峰. 基于数学形态滤波算子的三维时域Green函数求解方法研究. 电子设计工程. 2020(24): 24-28 . 百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 1191
- HTML全文浏览量: 644
- PDF下载量: 126
- 被引次数: 3
