0 引言

近几年来，四川盆地西南部的长宁地区地震活动十分频繁。尤其是2018年底以来，该地区地震活动强度明显增强，2018年12月16日和2019年1月3日建武向斜区域先后发生兴文M_S5.7和珙县M_S5.3地震，2次地震相距不到10km。几个月后，2019年6月17日北侧的长宁背斜区域发生了长宁M_S6.0地震。研究发现，研究区内建武向斜地区的地震活动和页岩气水力压裂具有较强的时空相关性(Lei et al，2017；Meng et al，2019)，而长宁背斜地区的地震活动和盐矿注水密切相关(阮祥等，2008；Sun et al，2017)。与兴文M_S5.7和珙县M_S5.3地震具有较少的余震活动不同，长宁M_S6.0地震的余震活动十分强烈，仅5级以上的余震就发生了4次，反映了研究区内复杂的地震活动背景。Zuo等(2020)对长宁地区地壳三维速度结构和地震活动特征的研究结果也显示，长宁背斜地区和建武向斜地区的地下介质结构和地震活动特征均存在明显差异。

地震是地壳岩石在应力作用下突然发生破裂或错动的结果，因此，探测和研究地壳介质和应力状态的时空演化过程是分析地震活动成因和机理的主要途径。一个地区不同阶段发生地震的震源参数在一定程度上反映了该区域构造应力场的变化(周少辉等，2018)，研究应力降的时空分布特征有助于认识地震孕育发生的物理过程。此外，有研究(Abercrombie et al，1993；华卫等，2010、2012；Hua et al，2013a、2013b、2015)认为对于诱发地震，受流体参与的影响，同等震级地震的应力降比构造地震小一个数量级，是识别诱发地震的指标之一。目前对于该认识仍然存在争议，认为其可能是由于不合适的衰减校正和假设引起的(Tomic et al，2009)。Pennington等(2021)通过使用多种方法对俄克拉马州2011年M_W5.7地震序列进行研究发现，对于高质量的数据，虽然不同方法得到的应力降其相对分布具有一致性，但是其绝对值存在系统性的偏差。而且，地质构造、震源模型、地震矩以及台站的选取等多种因素均可能对应力降的计算结果产生影响(Yu et al，2020)。通过使用相同方法对同一时间段、同一地区、不同类型地震活动的震源参数进行研究，可以最大限度地排除方法和数据选取的干扰，分析不同类型地震震源参数的差异。长宁地区近年来台站覆盖条件较好，地震活动频繁，且大量地震呈丛集状分布，有利于寻找波形相似的地震对，因此本文选择使用谱比法来计算长宁地区2018年12月1日至2019年7月21日期间地震的震源参数，据此分析该地区震源参数之间的相互关系，并探讨应力降的时空演化特征。

1 数据

我们收集了2018年12月1日至2019年7月21日长宁地区固定台站和流动台站记录的波形数据，台站仪器型号为BBVS-60、CMG-60T、CMG-3ESP-60、KS-2000M-60等，采样率均为100Hz。在此期间，中国地震台网共记录了研究区内(104.5°E~105.1°E，28.0°N~28.6°N) 13565个地震事件，其中M_L > 1.0地震5975个。本文使用Zuo等(2020)基于波速比一致性约束的双差层析成像方法(Guo et al，2018)联合反演得到的长宁地区地震精定位结果作为震源位置，研究区内台站和地震分布如图 1所示。

2 方法及流程

地震的震源参数包括地震矩、应力降、震源破裂半径等，这些参数描述了震源的静力学特征。在一定的震源模型下，震源参数是由地震观测谱拟合理论震源谱和震源谱参数而计算得到的(赵翠萍等，2011)。震源谱参数包括零频极限Ω₀和拐角频率f_c，其中Ω₀为震源谱的低频渐近线值，也称为零频极限值或震源谱振幅；f_c为震源谱的低频渐近线与高频渐近线交点处的频率。Ω₀主要反映地震的大小，f_c与地震破裂尺度有关。地震波能量在传播过程中不但随着传播距离出现几何扩散，还受到传播路径上介质的吸收和散射等，在到达台站接收仪器前，地表下方浅层的介质也会对地震波产生影响，因此观测波形是震源激发的信息经过上述各种过程的产物。在频率域，台站j记录到的地震i的观测位移谱O_ij(f)可以表示为

$ O_{i j}(f)=S_{i}(f) P_{i j}(f) G_{j}(f) I_{j}(f) S u r_{j} $

(1)

式中，f为频率；S_i(f)为地震i的震源谱；P_ij(f)为地震波从震源i到台站j的传播路径效应项，其描述地震波在传播过程中的衰减，包括地震波的几何扩散和非弹性衰减；G_j(f)为台站j的局部场地效应，其描述台站附近近地表地层介质对地震波的放大作用；I_j(f)为台站j的仪器响应函数；Sur_j为地表自由表面效应，可以根据不同体波的位移反射系数与入射角的关系，消除地表自由表面效应的影响。经过对仪器项I_j(f)和地表自由表面效应Sur_j的处理后，观测位移谱O_ij(f)可表述为

$ O_{i j}(f)=S_{i}(f) P_{i j}(f) G_{j}(f) $

(2)

由式(2)可知，由地震记录获得震源谱S_i(f)，必须消除表达式右边其他各项的影响，而P_ij(f)项与台站-震源的传播路径及频率有关，G_j(f)则与各个台站具体的场地及频率有关。在由观测数据恢复震源谱并开展震源参数的各项研究中，由于路径、场地等影响之间的相互耦合及不确定，许多研究或是将介质品质因子Q值设定为常数，或是不考虑场地的影响。

地震的理论震源谱模型也是震源参数研究中首先要考虑的问题。理论震源模型可以表示为(Boatwright，1978)

$ \varOmega(f)=\varOmega_{0} \frac{\mathrm{e}^{-\left(\frac{{\rm{ \mathsf{ π} }} f}{\varphi}\right)}}{\left(1+\left(f / f_{c}\right)^{\gamma n}\right)^{1 / \gamma}} $

(3)

其中，t为地震波从震源到台站的走时，n为高频衰减系数，γ控制震源谱的拐角形状。实际研究中，根据观测位移谱得到理论震源谱的研究方法主要有2种，第一种为直接对观测谱进行衰减和各台站场地响应校正的方法(刘杰等，2003)，第二种为经验格林函数谱比法。谱比法将目标地震附近震级相差大于1的小地震信号作为经验格林函数(EGF)，由于其到同一台站具有近似相同的传播路径，可以通过目标事件和EGF的频谱比来消除所有的衰减效应(Abercrombie，2015)

$ \varOmega(f)=\frac{\varOmega_{0}^{M}}{\varOmega_{0}^{E}}\left[\frac{1+\left(f / f_{c}^{E}\right)^{\gamma n}}{1+\left(f / f_{c}^{M}\right)^{\gamma n}}\right]^{1 / \gamma} $

(4)

其中，上标M和E分别代表主事件和经验格林函数。可见EGF方法可以准确获取目标地震的f_c。

使用多窗口谱估计方法(Thomson，1982；Prieto et al，2007)和台站3个分量的波形计算地震的震源谱。根据式(3)拟合单个地震事件的频谱得到该地震的零频极限Ω₀，利用式(4)拟合主事件和EGF的频谱比获得主事件的拐角频率f_c。对于4级以下和4级以上的地震分别使用S波震相2s和4s的时间窗，要求每个地震的信噪比不小于2。在计算Ω₀时，对观测波形扣除了仪器响应和地表自由表面效应。根据吴微微等(2016)计算得到的长宁地区台站的场地响应来扣除场地对Ω₀的影响。前人研究结果显示，式(3)中Q值的选取对Ω₀的影响较小(Clerc et al，2016；Holmgren et al，2019)，根据赵翠萍等(2011)和吴微微等(2016)对四川地区介质品质因子的研究结果，将Q值取为500。对于震源谱高频衰减系数n及拐角形状参数γ的选择，赵翠萍等(2011)认为ω^-2模型(Brune，1970)适合中国大陆的中小地震，但也有研究认为Boatwright模型更适合谱比法(Holmgren et al，2019)，且考虑到诱发地震与构造地震的震源谱可能存在差别，本文参考Onwuemeka等(2018)的做法，将n和γ分别确定在2~3和1~2范围内。为主事件选择EGF时，要求满足2个地震之间的距离不超过5km和所使用震相的波形互相关系数大于0.7，且至少有3条满足条件的震相外，还需要震级相差大于1。我们人工检查了所有的频谱拟合结果，并去除其中较差的结果。当一个地震有多个满足条件的EGF时，根据方差大小使用加权平均法(Abercrombie，2014)确定主事件的拐角频率。图 2为使用谱比法得到的一个地震的拐角频率。

得到零频极限和拐角频率之后，根据圆盘破裂模型(Eshelby，1957)并利用式(5)~(7)即可得到地震矩M₀、震源半径r和应力降Δσ等震源参数

$ M_{0}=\frac{4 {\rm{ \mathsf{ π} }} R \rho \beta^{3} \varOmega_{0}}{U_{\theta \varphi}} $

(5)

$ r=\frac{2.34 v_{\mathrm{s}}}{2 {\rm{ \mathsf{ π} }} f_{c}} $

(6)

$ \Delta \sigma=\frac{7 M_{0}}{16 r^{3}} $

(7)

其中，ρ为岩石密度，取2700kg/m³；β为震源深度处S波速度；R为震源距；U_φθ为辐射花样因子，对于S波取$\sqrt {2/5} $将主事件在每个台站的记录按式(3)和式(5)得到的地震矩的中值，作为该地震最终的地震矩。

3 结果与讨论 3.1 震源参数之间的关系

通过计算得到了长宁地区2018年12月1日至2019年7月21日442个M_L1.3~4.7地震的震源参数，其应力降位于0.02~7.26MPa范围内，超过90%的地震应力降小于2MPa。与赵翠萍等(2011)得到的中国大陆典型地区构造地震的应力降进行对比(图 3(a))，可以看出本研究提供了2.5级以下地震的震源参数样本，长宁地区M_L2.5~4.7的地震应力降在中国大陆典型地区构造地震的应力降范围内。图 3(b)~3(d)给出了长宁地区不同震源参数之间的关系，从图 3(b)可以看出，地震矩取对数后与M_L具有较好的线性比例关系，反映了本文计算得到的地震矩可靠性较高。式(8)和式(9)分别为长宁背斜地区和建武向斜地区地震矩和震级关系的最小二乘拟合结果。结果显示长宁背斜地区同等震级地震的地震矩整体略高于建武向斜地区，但差异较小，远小于同一区域内同等震级地震的地震矩变化范围(图 3(b))，该差异可能与两区域的构造环境有关。

$ \begin{aligned} \lg M_{0} =0.94 M_{\mathrm{L}}+10.15 \end{aligned} $

(8)

$ \lg M_{0} =0.91 M_{\mathrm{L}}+10.10 $

(9)

图 3(c)显示研究区内地震拐角频率介于2~21Hz之间，随着震级增大，拐角频率整体呈下降趋势。不同震级地震的应力降虽然变化范围较大，整体仍然呈现随震级增大而升高的特征(图 3(d))。发生在长宁背斜地区与建武向斜地区同等震级地震的应力降和拐角频率未表现出显著的差异性特征。式(10)和式(11)分别为长宁背斜地区和建武向斜地区拐角频率和震级关系的拟合结果，式(12)和式(13)分别为2个地区应力降和震级关系的拟合结果。拟合结果也显示长宁背斜和建武向斜2个区域地震的拐角频率、应力降与震级的关系差异较小。对于3级以下地震，建武向斜地区地震的应力降略低于长宁背斜地区，而对于3级以上地震，建武向斜地区可以计算应力降的地震数量较少，拟合结果并不可靠。且拐角频率与应力降取对数后与震级的线性关系也相对较差，因此长宁背斜与建武向斜地区应力降的差异还需进一步讨论。

$ \lg f_{c}=-0.15 M_{\mathrm{L}}+1.29 $

(10)

$ \lg f_{c}=-0.12 M_{\mathrm{L}}+1.26 $

(11)

$ \lg \Delta \sigma=0.50 M_{\mathrm{L}}-1.57 $

(12)

$ \lg \Delta \sigma=0.56 M_{\mathrm{L}}-1.72 $

(13)

一些研究认为应力降与震源深度存在相关性(Satoh，2006；Huang et al，2017)。图 3(d)展示了长宁地区地震应力降与震源深度(Zuo et al，2020)的关系，结果显示，除个别震级小且深度小于5km的地震具有异常高的应力降外，具有较大应力降的地震主要集中在5km以下深度，同时这些地震也基本具有较高的震级。对于相同震级的地震，应力降大小与震源深度的关系并不明显，与赵翠萍等(2011)的中国大陆主要地震活动区的震源参数研究结果一致。

3.2 长宁地区应力降的时空分布特征

考虑到研究区内长宁背斜地区和建武向斜地区的地震活动特征和发震成因存在明显区别，图 4(a)、4(b)分别给出建武向斜和长宁背斜区的地震震级和应力降的时间序列。由于较难获取页岩气开采时间和压裂过程的有关数据，故难以进行深入细致的分析。由图 4(a)可以看出，研究时段内建武向斜地区地震持续活跃，地震频次显著较高，该地区为长宁页岩气开采区，地震应力降多数为小于1MPa的较低水平，且在2次M_S > 5.0地震前后地震活动水平和应力降均无显著变化。研究期间仅有3次较高应力降(3MPa左右)的地震活动，其中一个高应力降地震位于研究区东南侧的一条已知断层上(图 5)，距离页岩气井相对较远；另外2个地震均位于建武向斜西北侧的一个地震丛集区，2019年4月之后高应力降地震较集中，可能反映了该处具有较高的应力积累。虽然该处也分布着多个压裂井，但更细致的分析需要搜集压裂作业的详细时间和数据。Yu等(2020)通过对加拿大水力压裂区地震的震源参数研究发现，靠近注入井的地震应力降(0.1~1MPa)比远距离地震的应力降(1~10MPa)低，可能是因为井附近较高的孔隙压力和裂缝密度使地壳岩石无法储存和释放更大的应力。Clerc等(2016)对水力压裂诱发地震的研究结果则显示其应力降为0.1~149MPa，处于典型的构造事件范围内，未发现地震应力降与井的距离之间具有相关性。Wang等(2020)的研究结果也显示水力压裂诱发地震序列的应力降处于构造事件通常观测到的范围内。本文的结果显示，长宁页岩气开采区的地震应力降总体较低，但是在开采井附近也有少量高应力降地震活动(图 5)。结合建武向斜地区2018年12月16日和2019年1月3日接连发生2次5级以上地震活动，这些高应力降地震可能反映了该地区具有较高的应力积累。

在长宁M_S6.0地震发生之前，长宁背斜地区可以计算出应力降的地震较少，应力降也不高，这些地震均位于长宁M_S6.0主震东侧。2019年3月在其北东侧出现了一个应力降相对较高的地震活动，可能反映了震源区的异常应力增强，是主震可能的前兆信号。长宁M_S6.0地震发生之后，短期内余震序列的应力降明显高于正常水平且快速衰减，与有关研究认为余震序列的应力降较大相一致。较高应力降的地震空间上集中在余震序列的西北段，主震周围和东侧地震的应力降相对较低。

为了排除震级的影响，图 6仅给出了M_L2.0~2.5地震应力降的空间分布，结果同样显示应力降高值异常主要分布于长宁背斜西北侧的珙县附近，反映了西北段为高应力区，揭示了几次5级左右强余震的高应力背景。Li等(2021)研究表明这种高应力是由于长宁M_S6.0主震引起的静态库伦应力增加所致。此外，地震精定位结果(Zuo et al，2020)显示在该位置存在与整体NW向展布的余震序列呈共轭分布的NE向的地震条带，说明该区域具有复杂的断层结构，这种交汇构造更容易产生应力集中(Andrews，1989；左可桢等，2018；张博等，2019)，导致发生高强度、高频次的余震活动。

图 6显示对于M_L2.0~2.5范围内的地震，在建武向斜地区仅有一个应力降相对较大的地震(接近1.5MPa)，其位于兴文M_S5.7地震的北侧，且距离页岩气开采井较远，其余地震的应力降普遍较低，该地区应力降平均值为0.40MPa；长宁背斜地区相同震级范围的地震应力降的平均值为0.54MPa，稍高于建武向斜地区地震的平均应力降值。长宁背斜地区计算出应力降的地震绝大多数为长宁M_S6.0地震的余震，总体距离双河盐矿较远，反映了构造地震的震源参数特征。建武向斜地区的地震活动则与页岩气水力压裂有关(雷兴林等，2020)。因此，建武向斜地区与长宁背斜地区应力降的差异可能反映了两地区地震活动成因和机理的差异。此外，位于长宁背斜区双河盐矿周围地震的应力降相比于背斜西侧的地震活动也处于相对低的范围，可能与盐矿长期注水导致的高孔隙压力降低了断层强度有关(Sumy et al，2017)。

本文结果显示，建武向斜地区地震的应力降相对长宁背斜地区要低，但是差异并不显著。Huang等(2017)研究认为诱发地震的应力降受构造应力控制。长宁页岩气开采区的诱发地震是由于流体扩散导致已有断层的活化引起的(Tan et al，2020)，因此诱发地震的应力降与发震断层上的应力积累状态有关。研究区内中强地震的连续发生表明该地区已经积累了足够高的应力，这可能是造成长宁背斜地区和建武向斜地区地震活动特征和发震成因存在明显区别，但是应力降差异较小的原因。

4 结论

本文使用谱比法计算得到长宁地区442个M_L1.3~4.7地震的震源参数，分析了不同震源参数之间的相互关系，并讨论了研究区内长宁背斜和建武向斜2个区域地震应力降的时空分布特征。结果显示：

(1) 长宁地区应力降位于0.02~7.26MPa范围内，平均值为0.69MPa。长宁M_S6.0地震后短期内部分余震的应力降显著较高，除此之外大多数地震的应力降均小于2MPa。

(2) 取对数后的地震矩与震级之间有良好的线性关系。随着震级增大，拐角频率整体呈下降趋势，应力降呈现随震级增大而增大的特征。应力降与震源深度之间无明显的相关性。

(3) 研究区内长宁M_S6.0地震发生后，短期内余震的应力降明显升高，随后快速衰减。地震应力降的空间分布特征显示，具有高应力降的地震主要集中在长宁背斜的西侧、M_S6.0地震的NW向余震条带的西北段，揭示出主震破裂造成西北段应力增强是几次M_S＞5.0强余震的主要成因。

(4) 对于M_L2.0~2.5地震，建武向斜页岩气开采区地震的应力降略低于长宁背斜地区，但差异并不显著。

致谢: 四川省地震局提供波形数据，文中图件使用GMT软件绘制，在此一并表示感谢。