2022, Vol. 42
2. 成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都市东三路1号，610059;
3. 成都理工大学数学地质四川省重点实验室，成都市东三路1号，610059;
4. 成都理工大学计算机与网络安全学院(牛津布鲁克斯学院)，成都市东三路1号，610059
2019-06-17四川长宁地区发生MS6.0地震，震中位于长宁县双河镇，震源深度16 km，极震区烈度达Ⅷ度，四川、云南、贵州、重庆等多地震感明显。此次MS6.0地震是四川盆地边缘地区历史最大地震[1]，使得长宁地区地震灾害研究再次回到地震学术界视野[2-3]。
2015年以来，长宁地区中等强度地震频发。部分学者认为这一系列地震与工业作业、流体注入有关[4-9]，也有学者认为与区域复杂的地质构造相关[1, 10]。关于此次长宁MS6.0地震的孕震机理目前仍存在较大争议。地震序列精确定位后的空间展布特征和震源机制反演特征是确定序列发震构造和判定孕震机理的有效方法，而长宁MS6.0地震发生后，国内外相关研究机构给出的地震震源机制解差异巨大。为此，本文采用宽频带波形资料和震相报告数据，使用Geiger重定位方法对长宁MS6.0地震序列空间分布特征展开研究，然后使用ISOLA近震全波形反演方法对序列中MS≥5.0地震震源机制进行反演，讨论本次长宁地震发震成因，以期研究该区域地震触发机制和孕震机理，为后期长宁地区工业注水活动提供参考。
1 区域地质概况长宁地区位于四川盆地南缘，扬子板块西缘，地震活动频繁。长宁MS6.0地震震中位于长宁-双河复式大背斜，背斜轴迹呈NW-SE方向展布，其中西北区域狭窄、东南区域宽广，背斜从西北端向西南端弯曲成近东西方向(图 1)。震中区以西为华蓥山断层带，由若干大小不等的断层组成；震中区以东发育多条平行褶皱构造带，对区域内沉积构造和构造变形作用明显[10]。震中区内褶皱断裂走向较为复杂，主要为NE、NW、EW向，可能与扬子板块周围构造带的相互挤压[11]以及多期构造作用有关[12]。
|
图 1 长宁区域地质构造[10] Fig. 1 Regional geological structure of Changning |
采用Geiger定位法对长宁地震序列(2019-06-17~07-04)进行震源重定位。Geiger定位法是通过地震波的走时和空间关系建立目标函数：
| $ \frac{\partial T_n}{\partial x_0} \Delta x+\frac{\partial T_n}{\partial y_0} \Delta y+\frac{\partial T_n}{\partial z_0} \Delta z+\Delta t=b_n $ | (1) |
式中，n=1，2，…，k，Δx、Δy、Δz分别为震源坐标x、y、z校正量，Δt为发震时间校正量，
选取迭代初值m0=(x0, y0, z0, t0)求导获得校正量，使用奇异值分解(SVD)最小二乘法求解，使时间残差逐渐最小化，从而获得最优定位解[13]，重定位流程见图 2。
|
图 2 地震震源重定位流程 Fig. 2 Flow chart of hypocenter relocation |
截至2019-07-04 11:00，共记录到MS2.0及以上余震563次。其中，5.0~5.9级地震4次，4.0~4.9级地震21次，3.0~3.9级地震83次，2.0~2.9级地震455次。根据Wen等[14]对四川区域地下一维层状速度模型的研究成果(表 1，表中M为莫霍面深度)，对长宁MS6.0地震及余震进行重新定位。
|
表 1 一维速度模型 Tab. 1 1-D velocity model |
为确保得到精确的重定位结果，对余震进行严格筛选，标准为：1)震级≥2.0；2)至少被6个地震台站所接收；3)将研究区划分成0.5 km×0.5 km×0.5 km小格，区块中如果含有地震事件，则只选择最大数量的台站所记录的地震事件，防止地震过度集中和维持相对统一的空间分布；4)每次地震至少包含4次P波或S波初至时间。剔除不符合条件的地震震中数据，最后得到研究区域463次地震的震中位置。
通过Geiger定位法得到长宁地震主震及余震震中分布(图 3)。从图中可以看出，地震序列总体沿NW65°方向展布，长约25 km，宽约10 km，呈条带状分布。地震序列震源深度表现为西深东浅、北深南浅，即沿着NW-SE方向为北西深、南东浅。深度统计结果(图 4)显示，重定位前余震震源深度在0~15 km之间，中小地震发震层位深度在15 km以上，主要分布于中上地壳，绝大多数地震深度分布在3~8 km，重定位后线性集中趋势更加显著，震中多集中分布在2~6 km之间，5次MS≥5.0地震的发震深度均在6~10 km之间。
|
图 3 长宁地震序列震源重定位分布 Fig. 3 Hypocenter relocation distribution of Changning earthquake sequence |
|
图 4 重定位前、后震源深度统计对比 Fig. 4 Comparison of hypocenter depth statistics before and after relocation |
反演数据来源于四川区域地震台网提供的长宁MS6.0地震以及4次MS≥5.0余震震相报告。采用ISOLA全波形反演方法进行震源机制反演[15]。ISOLA方法是将时间域迭代反卷积多点源矩张量反演方法发展到全波形矩张量反演中, 其中任意矩张量产生的地震记录均是基本地震记录的线性组合, 即
| $ \boldsymbol{d}(x, y, t)=\sum\limits_1^6 m_i \boldsymbol{G}_i(x, y, t), i=1, 2, \cdots, 6 $ | (2) |
| $ \mathit{\boldsymbol{d}}{\rm{ = }}\mathit{\boldsymbol{Gm}} $ | (3) |
式中，d为N个波形值的列向量(每个三分量台站包含的所有时间样本)，G为N×6系数矩阵，每列为G1、G2、…、G6(每列都与d以相同方式排序)，m为mi系数的列向量。
由于波形反演得到的超定方程d = Gm无精确解，故采用最小二乘法得到近似解。将方程写为：
| $ \boldsymbol{m}=\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{G}\right)^{-1} \boldsymbol{G}^{\mathrm{T}} \boldsymbol{d} $ | (4) |
经过最小二乘法和网格搜索法等对地震最优矩张量解进行反演。反演得到的矩张量m为对称方阵，可进一步将其分解为ISO(isotropic)、DC(double couple)以及CLVD(compensated linear vector dipole)[16]，即m = mISO+ mDC+ mCLVD。最后利用约化方差(VR)和条件数(CN)来快速评估理论波形和观测波形之间的相关性。
3.2 矩张量反演结果在利用ISOLA方法时，首先基于台站对地震事件包围相对完整、地震波形信噪比高、P波初动明显的原则，对长宁MS6.0地震区域50~200 km之间的台站进行筛选(图 5)。将震中坐标固定，在震源深度方向(Z方向)进行搜索，搜索深度范围为1~10 km，搜索步长设为1 km，每个台站的滤波频带范围设定为0.04~0.09 Hz。采用Liu等[2](图 6)的速度模型和离散波速法[17]计算全场格林函数，按照台站震中距合成理论地震图，在计算格林函数时，矩心时间在发震时刻5 s范围内以0.25 s作为搜索步长，并采用全矩张量求解模式进行计算，最终反演得到长宁MS6.0地震震源机制解(图 7~9)。
|
图 5 主震震中与地震台站分布 Fig. 5 Distribution of epicenter and seismic stations |
|
图 6 长宁地区地壳一维速度模型 Fig. 6 1-D crustal velocity model in Changning area |
|
图 7 矩心深度与互相关系数 Fig. 7 Centroid depth and cross-correlation coefficient |
|
图 8 长宁MS6.0主震矩心深度和波形互相关系数 Fig. 8 Centroid depth and cross-correlation coefficient of waveforms of Changning MS6.0 main shock |
|
图 9 观测波形与合成波形对比 Fig. 9 Comparison between observed and synthetic waveform |
从图 7、8可以看出，当矩心时间和矩心深度发生变化时，最终解的观测波形与合成波形的互相关系数也随之发生变化，当矩心时间为2.25 s、矩心深度为4 km时，互相关系数达到最高值0.916。图 9为最终解的观测波形与合成理论地震图。10个台站中，剔除在反演过程中波形拟合较差的三分量波形数据后，平均约化方差VR=84%。长宁MS6.0地震震源参数显示(表 2)，主震节面Ⅰ走向/倾角/滑动角分别为20°/72°/151°，节面Ⅱ走向/倾角/滑动角分别为120°/62°/21°，P轴方位角、倾角分别为71°、6°，T轴方位角、倾角分别为337°、33°。
|
|
表 2 长宁地震序列震源机制解 Tab. 2 Focal mechanism solutions of Changning earthquake sequence |
使用刀切法[18]评估矩张量反演参数的稳定性，对长宁地震序列震源机制解展开研究。首先加入反演所需的地震台站，根据Liu等[2]的一维速度模型，使用每次增减一个台站或增减一个三分量台站通道的方法，使得每个地震事件的所有断层均向对应的震源球进行反演，图 10为5次地震的所有解节线图。从图中可以看出，5次地震断层面节面线分布均较为集中，说明多次反演结果稳定，震源机制反演参数值可靠。
|
图 10 采用刀切法分析震源参数不确定性 Fig. 10 Uncertainty of source parameters by Jackknifing method |
对长宁地震序列重新定位后，发现地震序列总体沿NW65°方向呈长条状分布，震中多集中在2~6 km深度，长宁主震震中分布在余震区东南端，绝大多数余震位于MS6.0地震北西侧，主震附近及南东侧余震较少，呈现单侧破裂特征。研究地震序列震源深度整体分布趋势发现，地震序列表现为西深东浅、北深南浅特征，即沿着NW-SE方向为北西深、南东浅。根据长宁地震重定位后的震源深度特征，沿长宁地震破裂方向AA′作垂直剖面图[19](图 11)。
|
图 11 P波和S波速度垂直剖面图 Fig. 11 P and S wave velocity profiles |
图 11中色度代表相对速度扰动值，从图中可以看出，长宁MS6.0地震震中处于P波和S波高低速异常体的过渡带，P波和S波在横向上均表现出较强的不均一性，图中长宁MS6.0地震附近P波和S波速度高于左侧，说明研究区地下东南部比西北部岩层结构更加致密，力学性质更强，使得岩层不易发生破裂和错动。主震产生的地震波能量向西北部扩展，导致余震在主震西北方向密集分布。矩张量反演结果显示，节面Ⅰ走向/倾角/滑动角分别为20°/72°/151°；节面Ⅱ走向/倾角/滑动角分别为120°/62°/21°，震源深度为6 km。本研究得到的震源机制解与中国科学院地质与地球物理研究所给出的结果较为一致，结合重新定位后的地震分布特征，认为节面Ⅱ为长宁地震主震断层面，为发育在白象岩-狮子滩背斜的高倾角左旋走滑型地震。主震震源矩心深度为4 km，主压应力轴方位角、倾角分别为71°、6°，主张应力轴方位角、倾角分别为337°、33°，与余震展布特征一致。通过分析每次地震的断层面倾角特征可知，2019-06-17珙县MS5.1地震发震断层为逆倾滑动断层，其余3次余震发震断层均为逆冲型断层。矩张量反演结果如表 3所示，其中约化方差VR、条件数CN用来评判反演结果的正确性。从表中可以看出，5次地震DC含量均小于70%，其中2019-06-17珙县MS5.1地震DC含量更是低至22.3%，地震中存在较大的非双力偶成分，表明此次长宁MS6.0地震并非简单的构造地震。长宁地区注水活动频繁，而此次长宁地震序列震源深度较浅，震中多集中分布在2~6 km深度，5次MS≥5.0地震矩心深度均小于10 km，同时表 3中震源机制数据显示非双力偶含量(ISO+CLVD)占比很高，推测长宁地震的发生与工业注水活动息息相关。该结论与Wang等[8]和Li等[20]的研究结果一致。此外，长宁MS6.0地震序列震源机制解显示，发震断层走向以NW向和近NS向为主，P轴方位角以NE-近EW方向为主(表 2)，长宁地震序列总体表现出NE-SW向挤压，少量NW-SE向张拉特征。该特点反映出本次地震序列震源区域的主压应力方向应该为NE-SW向。但发生在长宁地震南侧的2018-12-16兴文MS5.7和2019-01-03珙县MS5.3地震呈现出NW-SE向挤压和NE-SW向拉张特征[1]，相隔不远却呈现出完全相反的应力场特征，说明四川盆地南部边缘地区构造变形特征具有显著差异，该区应力场转变复杂。
|
|
表 3 ISOLA矩张量反演结果 Tab. 3 The result of ISOLA moment tensor inversion |
1) 由地震序列重定位得到的463条震源数据可知，长宁MS6.0地震余震总体呈条带状分布在白象岩-狮子滩背斜。主震位于地震序列东南端，而绝大多数余震位于主震西北端，呈现单侧破裂特征。分析P波和S波速度剖面可知，长宁地震序列东南端地下岩层结构更加致密，力学性质更强，使得岩层不易发生破裂和错动。余震震源深度在0~15 km之间，绝大多数地震震中分布在2~6 km深度，沿NW-SE方向深度特征表现为北西深、南东浅。
2) 主震矩震级为MW5.78，震源矩心深度为4 km，断层面走向/倾角/滑动角为120°/62°/21°，为高倾角左旋走滑型地震。主压应力轴方位角、倾角分别为71°、6°，主张应力轴方位角、倾角分别为337°、33°，与余震展布特征相符。2019-06-17珙县MS5.1地震发震断层为逆倾滑动断层，2019-06-18长宁MS5.3、2019-06-22珙县MS5.4以及2019-07-04珙县MS5.6地震发震断层均为逆冲型断层。
3) 从长宁MS6.0地震主震和4次MS≥5.0余震震源机制反演结果来看，双力偶(DC)成分占比均不高，非双力偶(ISO+CLVD)成分占比很高，且5次MS≥5.0地震矩心深度均小于10 km，同时层析成像结果表明长宁地震处于P波和S波高低速异常体的过渡带，表明此次地震并非构造地震。此次长宁MS6.0地震序列应力场表现为以NE-SW向挤压为主，而发生在长宁区域南侧的2018-12-16兴文MS5.7地震和2019-01-03珙县MS5.3地震均表现为以NW-SE向挤压为主，两者应力场方向差别明显，表明长宁区域应力场较为复杂。综合分析认为，此次长宁地震为流体注入以及地质构造活动共同引发。
致谢: 感谢四川省地震局提供相关地震数据。
| [1] |
易桂喜, 龙锋, 梁明剑, 等. 2019年6月17日四川长宁MS6.0地震序列震源机制解与发震构造分析[J]. 地球物理学报, 2019, 62(9): 3432-3447 (Yi Guixi, Long Feng, Liang Mingjian, et al. Focal Mechanism Solutions and Seismogenic Structure of the 17 June 2019 MS6.0 Sichuan Changning Earthquake Sequence[J]. Chinese Journal of Geophysics, 2019, 62(9): 3432-3447)
(  0) |
| [2] |
Liu J Q, Zahradník J. The 2019 MW5.7 Changning Earthquake, Sichuan Basin, China: A Shallow Doublet with Different Faulting Styles[J]. Geophysical Research Letters, 2020, 47(4)
( 0) |
| [3] |
Li T, Sun J B, Bao Y X, et al. The 2019 MW5.8 Changning, China Earthquake: A Cascade Rupture of Fold-Accommodation Faults Induced by Fluid Injection[J]. Tectonophysics, 2021, 801
( 0) |
| [4] |
阮祥, 程万正, 张永久, 等. 四川长宁盐矿井注水诱发地震研究[J]. 中国地震, 2008, 24(3): 226-234 (Ruan Xiang, Cheng Wanzheng, Zhang Yongjiu, et al. Research of the Earthquakes Induced by Water Injections in Salt Mines in Changning, Sichuan[J]. Earthquake Research in China, 2008, 24(3): 226-234 DOI:10.3969/j.issn.1001-4683.2008.03.004)
( 0) |
| [5] |
Sun X L, Yang P T, Zhang Z W. A Study of Earthquakes Induced by Water Injection in the Changning Salt Mine Area, SW China[J]. Journal of Asian Earth Sciences, 2017, 136: 102-109 DOI:10.1016/j.jseaes.2017.01.030
( 0) |
| [6] |
Lei X L, Wang Z W, Su J R. The December 2018 ML5.7 and January 2019 ML5.3 Earthquakes in South Sichuan Basin Induced by Shale Gas Hydraulic Fracturing[J]. Seismological Research Letters, 2019, 90(3): 1099-1110 DOI:10.1785/0220190029
( 0) |
| [7] |
Meng L Y, McGarr A, Zhou L Q, et al. An Investigation of Seismicity Induced by Hydraulic Fracturing in the Sichuan Basin of China Based on Data from a Temporary Seismic Network[J]. Bulletin of the Seismological Society of America, 2019, 109(1): 348-357 DOI:10.1785/0120180310
( 0) |
| [8] |
Wang S, Jiang G Y, Weingarten M, et al. InSAR Evidence Indicates a Link between Fluid Injection for Salt Mining and the 2019 Changning(China) Earthquake Sequence[J]. Geophysical Research Letters, 2020, 47(16)
( 0) |
| [9] |
Jia K, Zhou S Y, Zhuang J C, et al. Nonstationary Background Seismicity Rate and Evolution of Stress Changes in the Changning Salt Mining and Shale-Gas Hydraulic Fracturing Region, Sichuan Basin, China[J]. Seismological Research Letters, 2020, 91(4): 2170-2181 DOI:10.1785/0220200092
( 0) |
| [10] |
常祖峰, 张艳凤, 王光明, 等. 2019年四川长宁MS6.0地震的地质构造成因——区域性构造节理贯通、破裂结果[J]. 地球学报, 2020, 41(4): 469-480 (Chang Zufeng, Zhang Yanfeng, Wang Guangming, et al. The Geological Genesis of the 2019 Changning MS6.0 Earthquake in Sichuan: Connecting and Rupturing of Regional Structural Joints[J]. Acta Geoscientica Sinica, 2020, 41(4): 469-480)
( 0) |
| [11] |
张岳桥, 董树文, 李建华, 等. 中生代多向挤压构造作用与四川盆地的形成和改造[J]. 中国地质, 2011, 38(2): 233-250 (Zhang Yueqiao, Dong Shuwen, Li Jianhua, et al. Mesozoic Multi-Directional Compressional Tectonics and Formation-Reformation of Sichuan Basin[J]. Geology in China, 2011, 38(2): 233-250 DOI:10.3969/j.issn.1000-3657.2011.02.001)
( 0) |
| [12] |
覃作鹏, 刘树根, 邓宾, 等. 川东南构造带中新生代多期构造特征及演化[J]. 成都理工大学学报: 自然科学版, 2013, 40(6): 703-711 (Qin Zuopeng, Liu Shugen, Deng Bin, et al. Multiphase Structural Features and Evolution of Southeast Sichuan Tectonic Belt in China[J]. Journal of Chengdu University of Technology: Science and Technology Edition, 2013, 40(6): 703-711)
( 0) |
| [13] |
张凌凡, 陈忠辉, 秦凡, 等. 岩石断裂混合矩张量反演与数值分析[J]. 矿业科学学报, 2019, 4(5): 394-402 (Zhang Lingfan, Chen Zhonghui, Qin Fan, et al. Hybrid Moment Tensor Inversion and Numerical Analysis of Rock Fracture[J]. Journal of Mining Science and Technology, 2019, 4(5): 394-402)
( 0) |
| [14] |
Wen B H, Yang H, Wei Q. Statistical Characteristics and Regression Model of Deep and Superdeep Velocity and Density Structure of Sichuan Basin, SW China[C]//CPS/SEG Beijing 2018 International Geophysical Conference and Exposition, Beijing, 2018
( 0) |
| [15] |
Sokos E N, Zahradník J. ISOLA a Fortran Code and a Matlab GUI to Perform Multiple-Point Source Inversion of Seismic Data[J]. Computers and Geosciences, 2008, 34(8): 967-977 DOI:10.1016/j.cageo.2007.07.005
( 0) |
| [16] |
Vavryčuk V. Moment Tensor Decompositions Revisited[J]. Journal of Seismology, 2015, 19(1): 231-252 DOI:10.1007/s10950-014-9463-y
( 0) |
| [17] |
Bouchon M. A Simple Method to Calculate Green's Functions for Elastic Layered Media[J]. Bulletin of the Seismological Society of America, 1981, 71(4): 959-971 DOI:10.1785/BSSA0710040959
( 0) |
| [18] |
Boyd O S, Dreger D S, Lai V H, et al. A Systematic Analysis of Seismic Moment Tensor at the Geysers Geothermal Field, California[J]. Bulletin of the Seismological Society of America, 2015, 105(6): 2969-2986 DOI:10.1785/0120140285
( 0) |
| [19] |
孙权, 裴顺平, 苏金蓉, 等. 2019年6月17日四川长宁MS6.0地震震源区三维速度结构[J]. 地球物理学报, 2021, 64(1): 36-53 (Sun Quan, Pei Shunping, Su Jinrong, et al. Three-Dimensional Seismic Velocity Structure across the 17 June 2019 Changning MS6.0 Earthquake, Sichuan, China[J]. Chinese Journal of Geophysics, 2021, 64(1): 36-53)
( 0) |
| [20] |
Li W, Ni S D, Zang C, et al. Rupture Directivity of the 2019 MW5.8 Changning, Sichuan, China, Earthquake and Implication for Induced Seismicity[J]. Bulletin of the Seismological Society of America, 2020, 110(5): 2138-2153 DOI:10.1785/0120200013
( 0) |
2. State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, 1 Dongsan Road, Chengdu 610059, China;
3. Geomathematics Key Laboratory of Sichuan Province, Chengdu University of Technology, 1 Dongsan Road, Chengdu 610059, China;
4. College of Computer Science and Cyber Security(Oxford Brookes College), Chengdu University of Technology, 1 Dongsan Road, Chengdu 610059, China