引言

长宁—兴文地区位于青藏高原东缘比较稳定的四川盆地东南缘，处于华蓥山断裂带东南侧多个褶皱带的过渡区（图1），是历史上的弱震区、少震区（易桂喜等，2019；Lei et al，2019a；解孟雨，孟令媛，2021）。据中国地震台网测定，2019年6月17日22时55分，四川省宜宾市长宁县发生M_S6.0地震，震中位于（28.34°N，104.90°E），震源深度为16 km。随后23时36分，在西北侧距离震中约11 km处发生了珙县M_S5.1地震。6月18日5时3分在M_S6.0地震震中附近发生长宁M_S5.3地震，6月22日22时29分珙县发生M_S5.4地震，7月4日在余震区西北端再次发生珙县M_S5.6地震，这几次地震显示出典型的震群特征。此次强震震群活动打破了长宁—兴文地区M6.0地震的空白，也给当地居民的生命财产安全造成了重大危害，引起了广泛关注。

图  1  双差定位中使用的台站、主要活动断层及研究区地震分布

Figure  1.  Distribution of the stations used in double difference location，active faults and earthquakes in the studied area

下载: 全尺寸图片 幻灯片

地震重定位结果能为地下断层分布和地震活动规律提供依据，是判定发震构造、确定震源特征的重要基础资料。应力降是指地震发生前后位错面的应力变化，一定程度上反映了该区域的应力状态，故应力降的时空演化可为震源区的应力场变化、强余震预测等研究提供参考依据。 长宁—兴文地区近年来发生数次大震，许多研究人员对该地区地震活动性变化的情况和地震危险性进行了深入探讨（Lei et al，2017，2019a；何登发等，2019；Tan et al，2020），但其结论并不一致。为了探究该地区在2019年6月17日M_S6.0地震之后的地震活动性特征及区域应力特征，本文首先利用双差定位方法（Waldhauser，Ellsworth，2000）对2019年至2020年长宁—兴文地区的地震进行重定位，以获得较为精确地震空间展布特征，然后利用研究区2019年6月17日至2020年12月31日发生的208个M_L3.0以上地震的波形资料计算拐角频率、地震矩、应力降等震源参数，并通过对区域内地震的震源参数与震级的定标关系和应力水平评估，探讨2019年6月17日M_S6.0之后震源区应力状态的改变及其与研究区内显著地震的关系，以期为进一步分析该地区的强震危险性提供参考。

1.   余震重定位

利用全国编目网提供的2019年1月1日至2020年12月31日长宁—兴文地区M_L0.0以上地震的震相报告，采用双差定位方法（Waldhauser，Ellsworth，2000）对小震进行重定位。双差定位方法的基本原理是，当事件i和j到达台站k和I的距离远远大于两事件震源之间的距离时，可将事件i和j组成一个地震对，此时可以认为台站k和I记录到该地震对（i和j）的到时差来源于事件i和j震源位置的差异。利用射线理论即可组建该地震对相对于台站k和I的到时差方程如下：

式中，Δt_jk和Δt_ik分别表示地震j和i到达台站k的走时残差，X_i和X_j表示震源位置矢量，S_jk和S_ik表示地震j和i到台站k的地震射线慢度矢量，ε_jk和ε_ik表示到时拾取误差。反演所有台站、所有地震对的到时差方程组即可得到地震丛集的每个事件相对于丛集矩心的定位，这种定位方法能够有效地减小地壳速度结构和速度模型对定位结果的误差影响。

本文选择震中距170 km以内的64个台站（图1）所记录到的Pg波、Sg波震相资料，其中震相到时数据存在明显误差的事件需人工删除。如图2所示，Pg和Sg震相的震源距-到时曲线基本为两条直线（黑色虚线），其中Sg波和Pg波的震源距-到时曲线的拟合斜率分别为0.299和0.176，之后删除偏离拟合直线±2 s的结果。分别赋予P波和S波的先验权重为1.0和0.5，初始速度模型采用Lei 等（2017）的一维速度模型，v_P/v_S设为1.73，选用最小二乘（least squares QR，缩写为LSQR）法进行重定位 ，最终获得1万9 667个地震事件的重定位结果。

图  2  相对定位前震源距-到时曲线及震相数据筛选

Figure  2.  Focal distance versus travel time curves before relative location and seismic phase data screening

下载: 全尺寸图片 幻灯片

长宁—兴文地区的地震活动主要受控于区域内的构造环境，研究区内主要的地质构造为白象岩—狮子滩背斜、长宁背斜和建武向斜（图3a），出露地表的断层较少，且规模、走向、倾角不一（阮祥等，2008；朱利锋等，2016；何登发等，2019；易桂喜等，2019）。从重定位后研究区M_L≥2.0地震相对定位的空间分布（图3a）中可以看出长宁地区的地震活动具有明显的分区特征，大致以28.3°N为界，分成南北两个区域：在北部背斜区，地震序列呈NW−SE向展布，长宁M_S6.0地震位于序列的东南端；南区地震分布相较于北部较为分散，主要集群分布在建武向斜东、西两端。图3b和3c绘制了重定位后沿经度、纬度方向的深度剖面，可以看出研究区的震源深度在18 km以内，优势分布处于4—6 km深度（图3d）。已有研究表明长宁背斜区内2019年长宁M_S6.0地震的震源机制解P轴主要沿ENE−WSW向，整体呈逆冲兼走滑的构造特征（易桂喜等，2019），而长宁南部的建武向斜区内，2018年兴文M_S5.7和2019年珙县M_S5.1两次强震的震源机制解P轴均沿WNW−ESE，呈走滑型的构造特征（Lei et al，2017，2019a）。易桂喜等（2019）计算长宁地区的中强地震应变结果显示，长宁M_S6.0地震序列震源区（北区）的发震断层整体以NS−SW向挤压为主，兼具一定的NW−SE向拉张分量，与Wu等（2015）基于GPS测量所得结果一致，但是南区的2018年12月16日M_S5.7和2019年1月3日M_S5.3地震呈现NW−SE向挤压、NE−SW向拉张，与北区有很大不同。从构造环境上来看，长宁地区位于川东南的盆山结合带上，处于构造变形模式的转换区，多个构造带叠加使得该区的构造环境十分复杂（何登发等，2019；易桂喜等，2019；解孟雨，孟令媛，2021）。胡幸平等（2021）反演了长宁地区地壳应力场的精细结构，其结果显示长宁地区地壳应力场的最大主应力轴基本呈近水平状态，但北部背斜区表现为逆冲型应力结构，南部向斜区表现为走滑型应力结构，且应力结构类型与对应区域内中强地震震源机制解的吻合度较低，在力学上甚至是相抵触的。结合上述研究成果可知，长宁地区南、北区的构造环境和应力场均有较大差异，也因此造成了地震活动分布和震源机制的不同。

图  3  重定位后的震中分布

（a） 重定位震中分布；（b，c） 重定位后沿经度、纬度方向深度分布；（d） 重定位后深度分布统计

Figure  3.  Distribution of relocated earthquake events

（a） Distribution of relocated epicenter；（b，c） Depth distribution along longitude and latitude，respectively；（d） Statistics of depth distribution after relocation

下载: 全尺寸图片 幻灯片

2.   应力降

应力降Δσ表示地震发生时破裂面上剪切应力的变化量，表达式为

式中，Δσ为静态应力降，σ₀和σ₁为地震发生前后断层面上的平均应力。应力降的大小与震源深度、区域应力条件、构造环境等诸多地震活动因素有关（Allmann，Shearer，2009；Uchide et al，2014）。长宁—兴文地区北部背斜区是我国重要的盐矿开采区，有多处注、取水采盐点（阮祥等，2008；张致伟等，2015），南部背斜区是重要的页岩气开采区，现今通过水力压裂技术进行页岩气工业开采（Lei et al，2019a，b；Meng et al，2019）。这些工业活动引发当地地震活动性的改变，属于诱发地震的范畴（阮祥等，2008；张致伟等，2015；Lei et al，2019a，b）。有研究认为，注水诱发地震的应力降均值低于构造地震应力降均值（Hough，2015），还有一些认为诱发地震和构造地震的应力降并无显著差异（Tomic et al，2009；Clerc et al，2016；Huang et al，2017）。由于研究区域、计算方法不同，加上应力降的计算也存在一定误差，对于诱发地震和构造地震应力降的异同点目前尚无定论，但是应力降的时空变化能够给出区域应力相对大小及变化的信息，这对于大震震源区的余震预测具有指示意义。

2.1   资料

本文利用中国地震台网中心提供的地震波形资料，计算了2019年6月17日至2020年12月31日长宁—兴文地区M_L≥3.0地震的地震矩和应力降等震源参数。这些地震中，M_L3.0—3.9地震251次，M_L4.0—4.9地震45次，M_L≥5.0地震7次。在近震源条件下，选取震中距小于300 km的台站记录到的波形数据，筛选出信噪比高的波形计算震源参数。这样，最终得到208个事件的震源参数。

2.2   方法原理

根据Wyss （1970），地震应力降表示为

式中r为震源半径，M₀为地震矩，这两个参数可通过地震波形分析和反演得到（李艳娥等，2012）。

对地震波形进行傅立叶变换得到观测谱，鉴于震中距较小不考虑介质衰减的影响，去除仪器响应、几何扩散校正，即可得到震源谱，表示为：

式中，Ω₀为震源谱零频极限，f_c为拐角频率。震源谱零频极限和拐角频率采用以下方式来确定。从三分向数字地震观测仪器记录中得到三个方向相应的震源谱，即分量震源谱。各分量震源谱矢量合成为位移谱Ω（f ）：

速度谱为

本文使用的数字化波形记录为速度记录，根据上式可以先计算速度谱，再得到位移谱。合成位移谱为

根据式（4）可知，低频段时，${ \varOmega } （ f ） $≈Ω₀，合成位移谱平方的积分为

则

根据实际波形数据得到位移谱之后，选择低频段水平部分（图4中f₁－f₂）进行积分，再根据式（9）得到零频极限值Ω₀。之后截取图（4）中f₁－f₃频段范围，将得到的Ω₀带入式（4），得到该频段范围内的理论值Ω_ti。根据相同频段范围内观测到的频谱数据Ω_i，利用

图  4  拐角频率拟合过程示例

L₁表示将Ω₀和f_c代入式（4）得到的频谱曲线，L₂为实际观测频谱曲线

Figure  4.  The example of corner frequency fitting process

L₁ represents the spectrum curve obtained by substituting Ω₀ and f_c into Equation （4），and L₂ is the observation spectrum curve

下载: 全尺寸图片 幻灯片

可以得到Ω_i和Ω_ti之间距离的平方和$ {R}_{\varOmega } $。结合式（4）和式（10）可以看出，R_Ω值将随频率f的变化而变化，当R_Ω达到最小值时，理论值与实际观测值的偏差达到最小，此时的频率称为拐角频率f_c。

实际计算时，首先挑选台站记录波形，再进一步挑选震源谱形状较好的数据进行反演计算。如图4所示，f₁和f₂可在低频段的平坦部分选取，f₃在衰减段中选取。只要震源谱的形状好，f₁，f₂和f₃的任意选取对计算结果产生的误差均不会对计算结果有本质影响（李艳娥等，2012）。

根据Aki和Richards （1980）可得

式中，M₀为地震矩，ρ为地壳介质密度（可取2.71 g/cm³），υ为波速（v_S取3.5 km/s），d为震源距，R为辐射因子（可用其均方根替代），r为震源破裂尺度。此时已知M₀和r，可根据式（3）计算应力降。

对某次地震的多个台站记录的波形数据，根据上述方法分别求得每个台站的震源参数，再获得多台站的平均值。为了消除个别台站异常值对平均值的影响，在计算多台平均值时采用如下方法（Archuleta et al，1982）：

式中，x_i为各台站的震源参数，N为台站数，Δx为误差因子。

3.   结果

3.1   震源参数的基本特征和关系

震源参数之间并非完全独立，震源参数之间的定标关系与地震的自相似性有关，自相似性也叫地震现象的标度不变性，因此震源参数的研究是数字地震学的重要部分，对于解释地震的物理意义具有重要意义（华卫等，2010，2012；赵翠萍等，2011；周少辉等，2018）。本文根据计算得到的长宁—兴文地区地震震源参数，分别得到研究区地震矩、震源半径、应力降、拐角频率与震级的关系，如图5所示。可见：应力降（对数）和拐角频率（对数）与震级的线性相关性较弱，但二者也呈现较明显的趋势变化，即随着震级的增大，应力降对数随震级增大而整体增大，拐角频率随震级增大而整体减小；地震矩（对数）、震源半径（对数）与震级的相关性较好。采用最小二乘拟合，得到如下关系式：

图  5  研究区震源参数与震级M_L的关系

Δσ为应力降，f_c为拐角频率，M₀为地震矩，r为震源半径，M_L为震级，R为相关系数，下同

Figure  5.  The relationship between source parameters and magnitude M_L in the studied area

Δσ is stress drop，f_c is corner frequency，M₀ is seismic moment，r is source radius，M_L is magnitude，R is correlation coefficient，the same below

下载: 全尺寸图片 幻灯片

从数值统计结果上看，研究区应力降范围在0.061—5.285 MPa，应力降均值为0.369 MPa，表现出变化幅度大、均值低的特征，拐角频率分布于0.36—4.41 Hz之间。

根据地震定位结果，长宁—兴文地区的地震活动具有明显的空间成丛现象，且南北两区强震的震源机制解也存在较大差异（易桂喜等，2019），为了便于进一步讨论震源参数，将研究区以28.3°N为界分为南北两个区域，通过对比两个区域震源参数的异同来讨论该地区的构造特点和发震原因。从数值上看，北区应力降处于0.061—5.285 MPa之间，均值为0.397 MPa；南区应力降处于0.071—1.075 MPa之间，均值为0.307 MPa。北部背斜区应力降均值略大于南部向斜区，与左可桢和赵翠萍（2021）的结果一致，但应力降数值整体上小于左可桢和赵翠萍（2021）的计算结果。这是由于应力降的计算结果取决于拟合震源谱的拐角频率，但其绝对大小有一定的误差，因此借由不同方法计算得到的应力降绝对值会有一定的差异。此外，研究区应力降整体偏低，90%以上的地震应力降小于1 MPa。已有研究也表明该地区多数地震的应力降、视应力低于四川盆地其它区域（程万正等，2006；张致伟等，2015；左可桢，赵翠萍，2021）。对于这一现象的解释需要结合研究区的工业开采背景。作为全国重要的盐矿、页岩气开采区，长宁地区已有多年的注水开采历史（Lei et al，2017，2019a，2020；Tan et al，2020），工业开采中向地下岩层注入流体会导致岩石孔隙压升高，从而降低断层强度，导致应力降、视应力等震源参数低于四川盆地其它地区（Houston，2001；Hough，2015；Huang et al，2017）。

采用最小二乘拟合，分别得到南、北两区震级（图6a，c）、震源半径（图6b，d）与震级的定标关系：

图  6  北区（a，b）、南区（c，d）震源参数地震矩M₀、震源半径r与震级M_L的关系

Figure  6.  The relationship between seismic moment M₀，source radius r and magnitude M_L in the northern region （a，b） and southern region （c，d）

下载: 全尺寸图片 幻灯片

从图6a和6c可看出南北两区的地震矩（对数）与震级之间均呈较强的线性关系，从拟合方程来看，南北两区地震矩与震级的关系差别小、相关系数高，说明计算的地震矩结果可靠。

3.2   应力降时空演变

从应力降与震级的关系（图7b）中可以看出，南区和北区的应力降整体都随震级的增大而增大，因此在讨论应力降演变中需要考虑震级的影响。本文通过对比同一震级档地震的应力降演化来减小震级的影响，以探究区域内应力降演变与地震活动之间的关系。图7a展示了M_L3.0—3.9地震应力降随时间的变化，可见：2019年6月17日主震发生之后，短期内应力降值明显高于平均值，并在一周之内发生九次M_L4.5以上地震，震源区应力得到释放，应力降也随之下降；10月后北区应力降再次升高，且明显高于南区；次年3月后南区应力降也逐渐升高。根据中国地震台网目录，2019年10月至2020年12月期间研究区共发生八次M_L4.3以上余震。

图  7  应力降∆σ随时间（a）和震级M_L（b）的变化

Figure  7.  Variation of the stress drop ∆σ with time （a） and magnitude M_L （b）

下载: 全尺寸图片 幻灯片

图8a标注了长宁M_S6.0主震和2019年10月至2020年12月期间的九次高应力降地震和八次M_L4.3以上地震。这里的高应力降地震是指应力降值大于等于研究区M_L3.0—3.9地震应力降均值2倍（即Δσ≥2$\overline{\mathrm{\Delta }\mathrm{\sigma }}$＝0.532 MPa）的事件。图8b绘制了研究区高应力降地震和显著地震的震中分布。由于重定位后部分事件丢失，编号a，b，d三次地震根据重定位前的目录绘制。在时间序列上，长宁M_S6.0主震后应力降明显高于背景值并快速下降，2019年10月后应力降再次升高，之后发生的较大震级事件与高应力降事件在时间上对应得较好：11月21日发生M_L3.4事件，应力降达到0.618 MPa，北区发生11月27日M_L4.6和12月29日M_L4.7；2020年3月7日发生M_L3.2事件，应力降达到0.696 MPa，北区发生4月23日M_L4.7；2020年9月6日发生M_L3.5事件，应力降达到0.671 MPa，2020年9月6日发生M_L4.3 （图8a）。南区自2019年6月17日长宁主震之后到次年3月，应力降保持在较低水平，2020年3月之后，南区地震活动性增强、应力降升高，之后发生两次M_L4.5以上地震，总体上看南区高值应力降与强余震在时间上未呈明显的对应关系（图8a），但是在空间分布上，从图8b可以看到编号为h和f的两次地震与高应力降地震的震源位置相近，而且研究区内高值应力降地震与区域内显著地震震中的距离在20 km以内。

图  8  （a） 北区和南区M_L3.0—3.9应力降随时间变化及强余震；（b） 2019年10月至2020年12月高应力降地震及其强余震的空间分布

Figure  8.  （a） Stress drop changes of M_L3.0—3.9 earthquakes and strong aftershocks with time in the northern region and the southern region；（b） Spatial distribution of high stress drop earthquakes and strong aftershocks from October 2019 to December 2020

下载: 全尺寸图片 幻灯片

进一步对M_L3.0—3.9地震的应力降值进行空间网格化，设定网格大小为0.05°×0.05°，且确保每个网格内至少有三个地震事件参与计算，求得网格内的应力降均值，如图9左所示。可见北区应力降高值分布于长宁M_S6.0地震序列的两端，南区的应力降高值分布在建武向斜两端。图9右中白色圆圈为研究区2019年6月17日至2020年12月31日内发生的M_L4.0以上地震震中分布。

图  9  M_L3.0—3.9应力降网格均值分布（网格大小0.05°×0.05°）

黑点表示本文用于计算应力降均值的地震，圆圈表示2019年6月17日至2020年12月31日期间长宁地区M_L4.0以上地震

Figure  9.  Mean grid distribution of stress drop of M_L3.0—3.9 （grid size 0.05°×0.05°） Black dots denote the earthquakes used for calculating mean stress drop，white circles denote the earthquakes above M_L4.0 in Changning area from 17 June 2019 to 31 December 2020

下载: 全尺寸图片 幻灯片

应力降能够反映震源区在地震中的应力释放水平，但地震释放的应力远小于背景应力。根据应力降的定义，震后应力降高值表明该区域整体应力水平较高。长宁M_S6.0地震后研究区M_L3.0以上地震应力降的计算结果显示，主震发生后震源区应力水平未恢复至背景值，应力降高值区仍存在较高的地震危险性，且M_L3.0—3.9地震应力降的演变与后续中强地震（图9圆圈所示）在时间、空间上都有较好的对应关系，这一结果表明去除震级影响所得的应力降变化对震源区的强震预测具有一定的指示意义。此外，本文得到应力降的地震事件数目较少，图8也显示部分高应力降地震与后续显著地震的时空关联不统一，目前工作在对后续强震的准确预测方面尚存不足。

4.   讨论与结论

本文研究区位于四川盆地内部，构造活动相对稳定，虽然区内没有大规模的活动构造，但2018年以来的数次强震和较高的地震活动性均表明长宁地区仍存在发生强震的可能，因此需要重新评估区域内的地震危险性（易桂喜等，2019；解孟雨，孟令媛，2021）。本文对长宁地区2019年6月17日至2020年12月31日M_L≥2.0地震进行重定位，并计算了在此期间208个M_L≥3.0地震事件的震源参数，还对研究区南区和北区在不同构造背景下的震源参数差异进行了对比，特别是重点从应力降时空演变的角度分析了研究区的应力特征及其与显著地震之间的关系，得出以下结论：

1） 长宁地区地震活动具有明显分区特征，大致以28.3°N为界分为南、北两个部分。北区的地震活动沿长宁背斜至白象岩—狮子滩背斜沿NW−SE方向呈条带状分布；南区的地震分布较为分散，呈现为分布在建武向斜两端的两个地震集群。

2） 长宁地区应力降处于0.061—5.285 MPa范围内，应力降均值为0.369 MPa，90%以上的地震事件应力降值小于1 MPa，北区的应力降均值略微大于南区。

3） 地震矩对数、震源尺度与震级之间具有良好的线性关系，应力降整体随震级的增大而增大，而拐角频率随震级的增大，整体上有减小的趋势。去除震级影响，M_L3.0—3.9应力降的时空演变显示，在长宁M_S6.0地震发生之后，该震级档地震的应力降高值与后续M_L4.0以上地震在时间序列上有较好的对应，且空间上相距20 km以内。此外，对M_L3.0—3.9地震应力降的空间分布网格化求得的网格均值显示，高应力降区与长宁主震震源区M_L4.0以上地震集聚区在空间上密切相关。结合以上结果，本研究认为M_L3.0—3.9地震应力降高值区仍存在较大的地震危险性，并且同震级档应力降时空演变对区域内中强地震的预测具有一定的指示意义。

感谢中国地震台网中心提供的计算数据，感谢中国地震局地震预测研究所周连庆研究员及匿名审稿人对本文提出的修改意见。