长宁地区复杂地震活动的力学成因分析

胡幸平， 崔效锋*， 张广伟， 王甘娇， Arno Zang，史丙新， 姜大伟

1 应急管理部国家自然灾害防治研究院， 北京 100085 2 江西省地震局， 南昌 330039 3 GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany 4 四川省地震局， 成都 610041

0 引言

四川长宁地区位于四川盆地南缘，扬子板块西缘，区域构造包括四川盆地与蜀南低陡褶皱带(王适择等, 2013)，是我国主要的井矿盐和页岩气产区之一，存在多处采盐、采气及废水回注井(阮祥等, 2008; 王玉满等, 2016).长宁地区地震频发且近年来增幅显著，据中国地震台网中心的数据，该地区在2018年12月以前并没有5级以上地震，而从2018年12月16日兴文MS5.7地震发生以来，陆续发生了多次5级以上地震，包括2019年1月3日珙县MS5.1地震，以及2019年6月17日长宁MS6.0地震及其后20天内的4次5级以上中强震(图1).由于毗邻盐气开采区，且注水诱发地震又是国内外长期以来的研究热点(Healy et al., 1968; 刁守中等, 1987; Segal, 1989; 张致伟等, 2012; Ellsworth, 2013; Bao and Eaton, 2016; Zang et al., 2019)，因此，长宁地区的地震活动及其成因机理引起了科学界乃至整个社会的广泛关注(阮祥等, 2008; 朱航和何畅, 2014; Sun et al., 2017; Lei et al., 2017, 2019; Meng et al., 2019；何登发等, 2019; 易桂喜等, 2019; 刘敬光等, 2019; 胡晓辉等, 2020; 常祖峰等, 2020).

针对长宁地区地震活动的已有研究指出，2019年6月17日长宁MS6.0地震及其余震序列发生于长宁背斜这一当地主要地质构造的轴部，其震源机制解的P轴主要沿北东东-南西西向，整体呈逆冲兼走滑的构造特征(易桂喜等, 2019)；而2018年12月16日兴文MS5.7地震和2019年1月3日珙县MS5.1地震发生于长宁背斜南边的向斜区，其震源机制解的P轴沿北西西-南东东向，且呈现为走滑型的构造特征(Lei et al., 2017, 2019).由此可见，这两组位于长宁背斜不同构造部位的地震活动，虽然其震源位置相距仅十几公里，但中强震震源机制呈现出明显不同的构造活动特征，表明长宁地区的地震活动具有很强的复杂性(图1).

震源理论认为，地震是在区域应力场作用下岩体沿某一构造面发生破裂的结果(Scholz, 2002)，而震源机制则反映了这一构造活动的具体形态.震源机制解会受到区域应力场的影响，但二者在本质上是不同的，其形态还在很大程度上受到局部发震构造面的影响.那么，长宁地区地震活动的复杂特征是在怎样的力学成因下形成的呢？仅仅是由于发震断层构造不同所造成的，还是孕育地震的应力场本身也存在局部变化？这一问题是研究长宁地区地震成因机理的基础，需要基于区域应力场的可靠结果开展力学分析才能加以解答；而这一应力场结果又需要满足以下两方面的条件.首先，区域应力场反演结果需要有足够的空间精细度，能够对长宁背斜轴部以及南部向斜区这两个相距仅十几公里的不同构造部位之间可能存在的应力场局部变化进行甄别.再者，为了避免应力场反演结果与中强震震源机制解之间的力学对比分析出现自我比较的逻辑问题(例如用中强震震源机制解反演得到的应力场再反过来跟这些中强震震源机制解做力学对比分析)，用于反演区域应力场的数据与这些中强震震源机制解应当在数据方面是相对独立的.目前我们尚未发现满足上述条件的研究成果见刊，而推进这一研究认识正是本文工作的宗旨所在.

图1 长宁地区地质构造及中强震震源机制解 黑色圆点表示2014年5月1日至2019年7月10日0.5级以上地震，红色五角星表示其中5级以上地震，蓝色沙滩球表示本文采用gCAP方法反演得到的震源机制解.图中地震定位结果来自国家地震数据中心的正式观测报告；地层数据引用自 中华人民共和国地质图1∶2500000(中国地质调查局, 2004).Fig.1 Geological map and focal mechanism solutions of moderate to strong earthquakes in Changning area Black dots represent the M≥0.5 earthquakes from May 1, 2014 to July 10, 2019, in which red stars mark those M>5.0 earthquakes. Blue beach balls represent the focal mechanism solutions from gCAP inversion in this paper. All the seismic locations in the figure are obtained from the official observation report of the National Seismic Data Center. Stratigraphic data is referred from the geological map of the People′s Republic of China (1∶2500000) (China Geological Survey, 2004).

1 研究方法

按照上述分析，本文选定区域应力场的反演方法需要考虑空间精细度以及与中强震震源机制相对独立这两方面需求.

利用大量小震P波初动极性求取综合震源机制解是研究区域应力场的重要途径(Aki, 1966;许忠淮等, 1987; 万永革等, 2011a).该方法可以利用大量无法计算出单个震源机制解的小震P波初动极性数据推断区域应力场的主要特征，因此，其所用数据与中强震震源机制解在数据方面是相对独立的.同时，小震P波初动极性在数据量和空间分布上优势显著，且小震震源尺度小，对局部应力场的变化敏感，使得该方法具有研究应力场精细结构的潜力.鉴于这些原因，本文选用小震综合震源机制解方法作为开展区域应力场反演的基础.

然而，如果仅单纯地通过小震综合震源机制解方法进行区域应力场反演，也难以获取能够满足本文区域应力场精细程度的可靠结果，主要体现在以下两个方面.首先，传统定位方法给出的小震定位误差往往有几公里(张国民等, 2002)，这一误差直接影响着能否正确合理地区分出不同构造部位的小震数据来研究区域应力场的精细变化.再者，地震定位和P波速度结构的误差还会造成P波初动数据的射线方向误差，这一射线误差对于近震数据尤为显著(Hardebeck and Shearer, 2002)；而小震综合震源机制解方法所用的数据由于震级小，往往正是以近震数据为主(许忠淮等, 1987; 万永革等, 2011a).虽然有学者认为，在数据十分充足的情况下，数据之间的误差可能会相互抑制，使得应力场反演结果相对稳定(盛书中等, 2015)；但对于精细尺度的应力场反演，需要将反演区域限制在比较小的范围内，此时P波初动极性数据的充足性难以保证，而数据误差的影响就难以忽视(胡幸平, 2018).鉴于上述原因，在利用小震综合震源机制解方法开展区域应力场的精细反演之前，有必要对长宁地区的小震定位以及波速结构进行修正，而双差层析成像方法(Zhang and Thurber, 2003)能够很好地契合这一研究需求.该方法是将双差定位方法(Waldhauser and Ellsworth, 2000)和传统走时层析成像方法相结合，利用邻近地震之间相对走时差和地震到台站之间的绝对走时进行震源位置和波速结构的联合反演，能够得到更为准确的地震定位和波速结构(Zhang and Thurber, 2006; 王长在等, 2018).因此，本文将引入双差层析成像的联合反演，来改进利用小震综合震源机制解方法开展的区域应力场精细研究.

2 双差层析成像联合反演

本文从国家地震数据中心搜集了2014年5月1日至2019年7月10日发生在27.84°N—28.84°N、114.40°E—115.40°E范围内0.5级以上地震事件的正式观测报告，从中提取P波和S波走时，以及P波初动极性.为了保证每个地震事件定位结果的可靠性，我们只选取300 km震中距以内有9条以上走时数据的地震，最终筛选出158个台站记录的19886个地震事件的P波走时数据139109条以及S波走时数据137322条(图2a,b).我们参考Lei等(2017)在长宁地区的研究中所给出的模型设定初始速度模型(图2c)，基于该初始模型利用TauP软件(Crotwell et al., 1999)计算出的理论时距曲线与筛选出的走时数据具有良好的一致性(图2b)，表明本文的初始模型以及筛选出的走时数据均比较合理.利用这些绝对走时数据，并设定邻近地震事件对的最大间隔为15 km，我们又得到了P波和S波相对走时差数据分别为808954条和806272条.

基于上述绝对走时以及相对走时差数据，本文利用Zhang和Thurber(2003)开发的tomoDD程序对长宁地区的地壳波速结构和震源位置进行了联合反演.反演采用带阻尼的LQSR算法，阻尼因子和光滑权重因子通过绘制数据方差和模型方差的均衡曲线来确定其最优值(Aster et al., 2012; 王长在等, 2018)，分别为500和40.此外，通过棋盘格检测板测试(Humphreys and Clayton, 1988)，我们发现长宁背斜区波速结构的横向反演分辨率能够达到0.1°×0.1°，因此将这一核心区域的反演分辨率设定为0.1°×0.1°(图3).

通过反演，本文得到了15762个地震的重定位结果(图4)，其横向及深度的平均相对误差分别为422 m和534 m.相比于初始定位结果(图1)，重定位结果在平面上更为集中，长宁MS6.0地震及其5级以上中强余震更清晰地呈现出沿西北方向线性分布的特征.在震源深度上，重定位后的地震主要集中在2～6 km，与Meng等(2019)对2015年2月至2017年12月发生在28°N—28.25°N、104.6°E—105°E范围内的地震所做的定位结果在总体上是一致的.如果进一步区分震级，本文结果还显示出M<4的小震震源深度仍集中在2～6 km，但M≥4的中强震则更多地分布在6～9 km.对照何登发等(2019)基于反射地震剖面给出的地层结构，可以看出，前者主要分布在基底之上的震旦系及以上地层，而后者主要分布在震旦系以下的基底中.

图2 (a)地震射线分布图; (b)震相走时曲线; (c)初始速度模型Fig.2 (a) Distribution map of seismic rays; (b) Travel-time curve; (c) Initial velocity model

图3 长宁背斜区棋盘格检测板测试结果(预设的速度扰动幅度为±5%)Fig.3 Results of checkerboard test with ±5% velocity perturbation around Changning anticline

与此同时，反演也给出了长宁地区的地壳波速结构.基于检测板测试的恢复分辨率大于0.6的反演结果，本文通过插值得到了长宁背斜区0.01°×0.01°的地壳波速结构(图5).由这一结果我们注意到，长宁背斜区在7 km以上主要表现为明显的S波高速异常.对此，本文整理前人对长宁地区地层构造所做的研究发现，长宁背斜区的地层厚度可达9 km(郭正吾等, 1996; 何登发等, 2019)；在漫长的地质时期中，由于受到来自西南方向的挤压应力，长宁背斜区的地层发生了褶皱、抬升，后又经历地表的风化剥蚀，形成了长宁背斜现今的构造面貌(何登发等, 2019).据此本文推测，正是由于地层的前期抬升和后期剥蚀，造成了长宁背斜区7 km以上的地层与周边区域更深部的地层相对应；由于S波速度一般会随着深度增加，那么这种横向上的地层差异就导致了在同一深度上长宁背斜区的S波速度相比于周边区域呈现为相对的高速异常.虽说由于盐气开采长宁地区存在大量的注水现象(Sun et al., 2017; Lei et al., 2019)，可能会造成S波低速异常(Roland et al., 2012; Guo et al., 2018)，但从本文结果来看，这些注水作用对于当地S波速度结构的影响并不显著，并没有颠覆不同地层S波速度的高低差异.

图4 长宁地区地震重定位平面(a)及深度分布(b)Fig.4 Horizontal (a) and vertical (b) distribution of the relocated earthquakes in Changning area

图5 不同深度反演结果与初始模型的相对差异 黑色圆点表示临近深度内的地震重定位结果，红色五角星表示其中5级以上地震.Fig.5 Relative differences between inversion result and initial model at different depths Black dots indicate the relocated epicenters close to each depth, in which red stars indicate M>5.0 earthquakes.

此外，我们也尝试了其他的初始速度模型(赵珠等, 1997)进行反演，得到了非常相似的地震定位和波速结构的反演结果，表明本文的双差层析成像反演结果对于不同的初始速度模型来说是比较稳定的.

3 区域应力场反演

在上述反演结果的基础上，本文仅选取长宁MS6.0地震发生前的M<4.0小震初动极性数据进行应力场反演，这一方面是为了保证应力场反演所用的数据与中强地震震源机制解的独立性，另一方面则是为了排除长宁MS6.0地震可能造成的局部应力调整.依据双差层析成像得到的P波速度结构及地震重定位结果，本文采用伪弯曲法(Um and Thurber, 1987)计算了这些极性数据的射线方向，获得了6936条修正后的P波初动极性数据.

本文首先采用网格扫描的方法对区域应力场的横向不均匀性进行了扫描.我们将长宁背斜区划分为0.1°×0.1°的网格；对于每一网格点，出于对应力场反演的区域平滑度和局部保真度的双重考虑，选取横向距离(DL)小于20 km的极性数据，并对其设定三次衰减函数的距离权重因子(WD)：

WD=(1-(DL/DL0)3)3,

(1)

式中DL0为横向距离的阈值，等于20 km.此外，考虑到清晰和不清晰的初动极性通常具有不同的拾取误差(Hardebeck and Shearer, 2002)，我们设定其清晰度权重因子分别为1和0.5.本文采用格点尝试方法(许忠淮等, 1983; 俞春泉等, 2009)来求取小震综合震源机制解，搜索步长为2°×2°×2°，选取加权矛盾比最小的200个解作为可选解.通过对这些可选解的三个主应力轴构成的张量进行张量平均，然后求解平均张量的主向量，本文最终确定出网格点上每个主应力轴的平均方向；此外，通过计算网格点上每个主应力轴与可选解的对应主应力轴之间的空间夹角的均方根，本文得出网格点上每个主应力轴的方向离散度.为了保证结果可靠性，本文只对极性数据超过500条的应力反演结果进行展示和分析(图6).

从图6可以看出，除了少数边缘网格点之外，大部分结果的主应力轴的方向离散度都在15°以内，表明这些反演结果是相当稳定的.区域应力场的网格扫描结果显示出，长宁地区地壳应力场的最大主应力轴(σ1)在整个研究区内基本都呈现为低倾角的近水平状态，其方位大体沿着近东西向，但在长宁背斜轴部往北以及南部向斜区这两个区域之间存在不太明显的分叉现象——前者偏向北东东-南西西向，而后者偏向北西西-南东东向.相比之下，这两个局部区域应力场在构造应力类型上的差异更为显著，其中长宁背斜轴部主要呈现为中等主应力轴(σ2)低倾角、最小主应力轴(σ3)高倾角的逆冲型应力结构，而南边的向斜区则主要呈现为σ2高倾角、σ3低倾角的走滑型应力结构.

上述局部应力差异大致以28.3°N的纬度线为界，为了更清晰地对比这两个区域的应力场差异，本文将修正后的P波初动极性数据分成28.3°N以南和以北两组(简称为南组和北组)，对于每个数据组内不考虑距离权重，同样采用上述格点尝试法和参数设定求解小震综合震源机制解(图7).分组结果显示出，南北两组的最大主应力轴(σ1)的平均倾角都较低(南组10.1°，北组12.9°)，接近水平状态，其平均方位虽然从北向南显示有17°左右的顺时针偏转(北组79.4°，南组276.6°)，但从二者可能解的σ1方位分布来看，两组σ1方位差别并不显著，可以说都大致沿着近东西向.与σ1轴相比，σ2和σ3轴呈现出的分组差异明显超出了其可选解的分布范围，其中北组表现为σ2低倾角、σ3高倾角的逆冲型应力结构，而南组则表现为σ2高倾角、σ3低倾角的走滑型应力结构，而且这种构造应力类型的差异是相当显著的.

本文反演得到的南部向斜区内的最大主应力方位与Lei等(2017, 2019)利用该地区震源机制解的反演结果以及王玉满等(2016)整理长宁气田区地应力测量资料得到的最大水平主应力方位比较一致，而长宁背斜轴部的最大主应力方位则与Zheng等(2018)和石玉涛等(2013)给出的该地区北东-南西向快波方向比较接近.

图6 长宁地区网格点上主应力方向及其离散度 黑线代表主应力方向，其长度与主应力轴的倾角成反比，色阶表示每个主应力轴的方向离散度. 黑色圆点表示提供P波初动极性数据的M<4.0小震重定位结果.Fig.6 Orientations of principal stress axes and their direction dispersions determined by grid-scanning in Changning area Black bars indicate the orientations of principal stress axes, and their lengths are inversely proportional with the plunges of the principal stress axes. The direction dispersion of each principal stress axis is color coded. Black dots indicate the relocated M<4.0 earthquakes providing initial motion data of P wave.

图7 长宁地区应力场分组反演结果 左图中圆点表示提供P波初动极性数据的M<4.0小震重定位结果；右图中-、o、+分别表示可选解的σ1、σ2、σ3轴，平均解的 σ1、σ2、σ3轴的方位角和倾角标注在对应应力轴符号之下.北组和南组的地震和应力场反演结果由红色和绿色加以区分.Fig.7 Stress inversion results determined by discriminating seismicity groups in Changning area Dots in the left map indicate the relocated M<4.0 earthquakes providing initial motion data of P wave. -, o, + in the right polar graph indicate optional σ1,σ2, and σ3 axes respectively; the average σ1,σ2, and σ3 axes are marked namely by the symbols of σ1,σ2, and σ3, and their azimuth and plunge are written below. Red/green color discriminate the relocated earthquakes and stress inversion results of north/south seismicity group.

4 区域应力场及中强震震源机制解的力学拟合分析

为了研究区域应力场与中强地震的力学一致性，同时也为了进一步验证局部应力差异，本文将上述两个局部区域的应力结构与中强震震源机制解进行了力学对比分析.对此，本文采用的是计算归一化剪滑分量(Angelier, 2002)的方法.归一化剪滑分量(ω)同时考虑两个参数，一是相对剪应力大小，即区域应力场(应力张量)在震源机制解节面上的计算剪应力除以最大剪应力(胡幸平, 2018; 万永革, 2020)，二是剪滑角，即实际震源机制解节面上的滑动方向和计算得到的剪应力方向之间的角度(Gephart and Forsyth, 1984; 许忠淮, 1985; Michael, 1987).归一化剪滑分量(ω)的变化范围为-1～1，值越小，表明很难用区域应力场解释实际震源机制解所指示的岩体滑动，若为负值，则表明区域应力场与实际震源机制解在力学上是相抵触的.本文之所以使用归一化剪滑分量来衡量区域应力场与震源机制解的一致性，而不是单纯使用剪滑角，是因为后者仅考虑了如果岩石在区域应力场下发生滑动，那么这种滑动与实际震源机制解之间的差异，而前者还同时从剪应力量值的角度考虑了岩石在区域应力场下发生滑动的能力或者说概率.如果区域应力场作用在震源机制解节面上的相对剪应力很小，说明岩石滑动很难被驱动，那么即便是剪滑角很小，也难以认为区域应力场与震源机制解在力学上是一致的.另外，归一化剪滑分量的计算对于震源机制解两个节面中的任意一个，其结果都是相等的，这也避免了从震源机制解两个节面中判断发震断层面的难题.

本文用于对比分析的中强震震源机制解是我们采用国际上常用的gCAP方法(Zhu and Ben-Zion, 2013)反演得到的，所用的波形数据来自中国地震局地球物理研究所国家测震台网数据备份中心.在震源机制反演中，Pnl和S波截取波形窗长分别为35 s和70 s，相应滤波范围为0.02～0.2 Hz和0.02～0.1 Hz；走向、倾角和滑动角的搜索步长为10°，深度为1 km.格林函数采用频率-波数法(FK)来计算(Zhu and Rivera, 2002)，其中采样间隔设为0.1 s，采样点为1024个.震源机制解的反演结果见图1及表1.

表1 长宁地区中强震震源机制解反演结果Table 1 Focal mechanism solutions of moderate to strong earthquakes in Changning area

计算归一化剪滑分量需要用到相对应力大小R；按照Gephart和Forsyth(1984)以及万永革等(2011b)的定义，其表达式为

R=(σ2-σ3)/(σ1-σ3),

(2)

R的变化范围为0～1，但利用小震综合解的反演并不能确定出其具体参数(许忠淮等, 1983; 万永革等, 2011b).对此，本文首先将其设定为0.5这一中间值来进行计算分析(图8).结果显示出，对于对应区域内的应力场和中强震震源机制解(例如同样是长宁背斜轴部的应力场反演结果和中强震震源机制解)，计算得到的归一化剪滑分量都比较大(13个震源机制解中有10个的ω值大于0.6，ω最小值为0.29)，这就表明二者之间在力学上具有良好的一致性，同时这也意味着驱动小震和中强震的应力场是几乎同一的；与之相比，对于非对应区域的应力场和中强震震源机制解(例如长宁背斜轴部的应力场反演结果和南部向斜区内的中强震震源机制解)，计算得到的归一化剪滑分量明显偏小(13个震源机制解中有12个的ω值小于0.3，甚至有6个ω值为负值)，这就表明二者之间在力学上的吻合度较低，甚至是相抵触的.这一强烈的对比结果表明，长宁背斜区的地壳应力场确实存在比较明显的局部差异，而这种应力场的局部改变是造成背斜轴部和南部向斜区内的中强震活动机制存在差异的必要力学基础.

图8 长宁地区中强震震源机制解与 区域应力场的力学一致性Fig.8 Results of mechanical fitting between the stress fields and focal mechanism solutions of moderate to strong earthquakes in Changning area

从图8中也能发现第11号震源机制解在非对应计算中反而会得到更大的归一化剪滑分量.本文认为此个别力学拟合结果的反置可能由于震源机制解反演误差或者其反映的是应力场调整等原因导致的，而这种极个别现象并不能颠覆应力场与中强震震源机制的对比分析认识.

对于长宁地区的相对应力大小R，Lei等(2019)以及刘敬光等(2019)通过震源机制解反演给出了相对较低的取值：前者结果为0.15；后者在R=(σ2-σ1)/(σ3-σ1)的定义下的结果为0.7，如果按照本文R的定义，其结果为0.3.因此，本文也将R设为0.15、0.3等低值，以及试算了0.8这一高值，此外还采用其他学者的震源机制结果(易桂喜等, 2019)进行了多组计算分析.这些设定不同相对应力大小以及采用不同来源的震源机制数据的计算结果都表现出统一的对比趋势，即利用小震综合震源机制解反演得到的应力场和中强震震源机制解在对应区域内一致性较好，但在非对应区域之间一致性明显降低.因此，可以认为本文的对比分析是相当稳定的，而且这一力学分析能够说明长宁背斜轴部和南部向斜区之间确实存在局部应力场的改变.

5 讨论

5.1 长宁地区应力场局部变化的成因

通过长宁背斜区应力场的精细反演及其与中强震震源机制解的力学对比分析，本文发现，长宁地区地壳应力场存在显著的局部变化，长宁背斜轴部表现为最大主应力轴北东东-南西西向的逆冲型应力结构，而南部向斜区则表现为最大主应力轴北西西-南东东向的走滑型应力结构.那么在相距仅十几公里的两个构造部位，应力场为何会表现出这样的差别？

细致观察上述局部应力结构(图7)可以看出，虽然二者在最大主应力方位上也有一定的差异，但更显著的差异在于前者的中间主应力轴近水平、最小主应力轴近直立，而后者与之相反.这就表明相比于南边的向斜区，长宁背斜轴部在横向上受到的挤压作用更强.本文认为这种横向挤压作用的差异可能是由于岩体重力在横向上产生的围压不同而带来的.根据何登发等(2019)对长宁背斜区地层系统的研究，长宁背斜轴部10 km以上的地壳主要由寒武系、震旦系以及更老的结晶基底构成，岩体中包含相当比例的震旦系白云岩和基底变质杂岩；而南部向斜区10 km以上的地壳岩层中有3 km以上厚度的奥陶系至侏罗系沉积岩，岩层组成中页岩所占比例高于前者，而也正是如此，南部向斜区成为长宁气田的集中开采区(王玉满等, 2016).在岩石的力学性质上，一般页岩的泊松比明显低于白云岩和变质杂岩，因此有理由认为，长宁背斜轴部岩体的重力由于泊松效应会在横向上产生更大的围压，而这一横向围压叠加在周边构造活动产生的横向推挤作用后，使得水平方向上的最大和最小主应力(SH、Sh)增加从而成为最大及中等主应力(σ1、σ2)，而垂向应力(SV，可能主要由上覆岩体重力产生)成为最小主应力(σ3)；相比之下，在南部向斜区，岩体重力由于泊松效应产生的横向围压较前者明显更小，垂向应力(SV)仍为三个主应力中的中等主应力(σ2).于是这种由于岩石泊松比差异带来的横向围压差异就导致了长宁背斜轴部和南部向斜区这两个局部区域之间存在逆冲型(SH>Sh>SV)和走滑型(SH>SV>Sh)的应力类型转变.

为了检验上述成因的可能性，本文依据双差层析成像反演的波速结构计算波速比(VP/VS)，然后依据波速比与泊松比(μ)的关系，即

μ=(0.5×(VP/VS)2-1)/((VP/VS)2-1),

(3)

计算了长宁背斜区的地壳泊松比分布(图9).

从图9可以看出，长宁背斜区的岩石泊松比具有明显的横向差异；这种差异在7～10 km的深度上表现得尤为显著，长宁背斜轴部的泊松比明显高于南部向斜区，差值甚至可以达到0.15左右，而这种岩石力学性质的强烈差异很可能就是由地层的横向差异所造成的.另一方面，王玉满等(2016)通过整理长宁气田的地应力数据指出，气田区2300～3200 m深度上的横向差应力值(SH-Sh)为21.4～22.3 MPa，垂向应力(SV)为56～66 MPa.依据上述这些数据，本文对长宁背斜区的地应力量值以及横向围压差异的影响进行估算.

吴珍汉和白加启(1997)根据多种应力资料的统计分析指出，上地壳的最大差应力在25 MPa左右；由此看来，长宁气田区深度3 km左右的横向差应力值(22 MPa)已经处于超高压状态.因此，有理由认为该地区上地壳内的横向差应力不会再随深度明显增加，于是我们将深度7 km左右的横向差应力按25 MPa估计.选择长宁气田所在的南部向斜区的应力结构作为区域应力场的参考基准，根据反演结果，该区应力结构呈现为走滑型，那么这种横向差应力(SH-Sh)与(σ1-σ3)接近.垂向应力(SV)在深度3 km左右约为60 MPa，与上覆岩体的重力相当，据此，我们将深度7 km左右的垂向应力(SV)估算为140 MPa.背斜轴部和南部向斜区的泊松比差值按0.1估算，那么前者由于泊松效应在横向上产生的围压比后者高出14 MPa左右.如果相对应力大小R按0.5估算，则σ2-σ3=(σ1-σ3)×R≈12.5 MPa，可以看出这一由于泊松比不同导致的～14 MPa的横向围压差别具备了造成σ2和σ3在Sh和SV之间转换的量值需求；而如果相对应力大小按照0.15(Lei et al., 2019)或者0.3(刘敬光等, 2019)来估算，则σ2-σ3等于3.75 MPa或者7.5 MPa，那么横向围压差别更是明显超出该值，就更有可能导致σ2和σ3在Sh和SV之间的转换.

图9 长宁背斜区波速比及泊松比分布 黑色圆点表示临近深度内的地震重定位结果，红色五角星表示其中5级以上地震.Fig.9 Poisson′s ratio and velocity ratio at different depths around Changning anticline Black dots indicate the relocated epicenters close to each depth, in which red stars indicate M>5.0 earthquakes.

通过上述应力量值的估算，可以看出，背斜轴部和向斜的岩体在泊松比这一力学性质上的差异，确实能够造成这两个局部应力场具有不同的应力类型.同时，这种岩石泊松比差异造成的围压差在长宁背斜轴部到向斜之间的横向作用，很难说能够与区域应力场在水平面上的任一主应力(SH或Sh)方向一致，于是也会导致局部最大主应力方向发生一定程度的改变(Zoback，1992).前人的应力场数值模拟研究也证明了岩石介质差异具备造成应力类型和应力方向发生改变的这种可能性(刘平江等, 2007).

另外，结合图6和图9还可以发现，南部向斜区内部也依然存在更小尺度的局部应力结构和岩石泊松比的对应差异.例如，兴文MS5.7地震和珙县MS5.1地震附近呈现为走滑型应力结构和低泊松比，而其西南边(104.7°E、28.1°N附近)呈现为逆冲型应力结构和高泊松比.由此可见，在长宁背斜轴部与南部向斜区之间、以及南部向斜区内的更小尺度上，局部应力结构和岩石泊松比之间都表现出很强的对应关系，这就进一步支持岩石泊松比的横向差异很可能就是局部应力结构差异的主要成因.

除了岩石泊松比的横向差异外，其他更为复杂的因素也有可能在一定程度上引起或者说加剧了长宁背斜区的应力场局部变化.例如，长宁背斜区由于盐气开采不断发生的高压流体注入.鉴于背斜轴部采盐和南部向斜区采气的注水工艺可能不同、流体渗透可能受裂隙发育方向的影响、以及注水过程的持续动态性，该地区的流体注入过程造成的地应力变化很可能不是静水压力那样的各向同性，就有可能导致这两个区域的应力方向和结构形态产生差异.再如，现今构造作用以及岩体中残留的古应力场释放在长宁背斜这一古褶皱构造的不同部位是否也存在局部变化(常祖峰等, 2020)，从而导致了现今应力结构的不同，而这种因素的影响也不能完全排除.

上述这些都是造成区域应力场局部改变的可能因素；同时，另一方面，无论是用小震极性还是震源机制解反演得到的都是地震释放的应力，其主应力轴会受局部断层构造影响，可能与实际应力场的主应力轴有所差别(Tian et al., 2013).虽说通过大量不同走向和倾角的断层上的地震数据能够有效地推断出区域应力场的方向特征(许忠淮等, 1983)，且本文所用的小震数据丰富，可以认为上述条件大致成立，但仍难以保证本文反演结果是对实际应力场的“完美”恢复.这种地震释放应力与实际应力场的可能差异也是目前所有基于地震学的应力场反演方法无法完全消除的问题，而这也有可能是造成应力场反演结果出现南北差异的部分原因.

通过上述分析可以得出，岩石泊松比的横向差异具备造成长宁背斜区地壳应力场的局部变化的可能性，因此，这种局部应力变化是可以理解的.但与此同时，这种局部应力变化也有可能是多方面的复杂因素共同作用的产物，甚至还可能包含了少量由于地震学反演造成的“假象”，因此，需要我们结合多种数据和手段、开展更为全面和深入的研究工作才能彻底解答.

5.2 长宁地震发震构造

基于上述双差层析成像联合反演结果、震源机制解、以及区域应力场结果，本文也对长宁地震的发震构造进行了探讨.

从重定位和震源机制解的结果可以发现，长宁地震序列的中强震沿着北西方向线性分布，其震源机制解也几乎都具有北西向的节面，因此，本文选取了2条穿过长宁MS6.0地震的剖面来讨论这一地震的发震断层构造.其中一条剖面沿着长宁地震序列线性分布的北西方向(图4中的AA′，走向NE130°)，另一条与之正交(图4中的BB′，走向NE40°)；开展剖面分析的底图数据包括2个，一是地震相对密度(图10)，另一则是波速结构(图11).

本文地震相对密度的计算方法为：将整个研究区划分为0.01°(经度)×0.01°(纬度)×1 km(深度)的三维网格，对于每个网格点，统计其周边2 km范围内不同震级区间的地震数量，然后通过除以所有网格点周围对应震级区间地震数量的最大值进行归一化，得出各个网格点上的不同震级区间的地震相对密度. 从图10的AA′剖面可以看出，震源深度整体上向北西方向逐渐加深，且不同震级地震的总体震源深度存在分离现象，中强震主要分布在小震底部.对照何登发等(2019)给出的地层结果，可以看出4级以上的中强震主要分布在基底中.由于这些中强震的破裂长度可达数公里(Wells and Coppersmith，1994)，且它们的发震节面和空间位置具有一定的连贯性，因此，本文认为在长宁背斜轴部可能存在北西走向的基底断层，控制着4级以上中强震的活动.这种基底断层可能是一条较大规模(长度超过15 km)的断层，又或是多条中等规模(长度数公里)的断层组合而成；但无论哪种方式，其总体上都基本贯通(或者说组合贯通)了长宁背斜的轴部.依据AA′剖面上的地震相对密度和中强震震源机制解(图10)，并结合该剖面的S波速度结构(图11)，本文勾画了这种基底断层沿着北西向AA′剖面的构造特征，认为其主要分布在6～9 km深度内，且从长宁MS6.0地震处往西北方向略微加深.相比之下，小震几乎都分布在中强震震源以上的地壳浅部，其优势分布深度与中强震存在一定程度的分离，且小震破裂尺度小，因此本文认为这些浅部小震可能并不对应着规模较大的基底断层的活动，而是上覆地层中的裂隙或者微小断层活动而产生的.

图10 垂直剖面上的地震相对密度及震源机制投影 AA′和BB′剖面位置见图4；红色虚线表示推断的断层构造.Fig.10 Relative density of seismicity and projection of focal mechanism solutions on vertical profiles The locations of profiles AA′ and BB′ are sketched in Fig.4. Dashed red lines indicate the inferred seismogenic fault.

图11 垂直剖面上的波速结构及震源机制投影 AA′和BB′剖面位置见图4；黑点表示剖面两侧2 km以内的地震，红点表示4级以上地震；红色虚线表示推断的断层构造.Fig.11 Seismic velocity and projection of focal mechanism solutions on vertical profiles The locations of profiles AA′ and BB′ are sketched in Fig.4. Black dots indicate the relocated earthquakes within 2 km on both sides of each vertical profile, while larger red dots indicates those M≥4.0 events. Dashed red lines indicate the inferred seismogenic fault.

另外，BB′剖面的波速结构似乎呈现出向西南缓倾的低波速异常(图11)，有可能指示着断层附近的岩体破碎.因此，本文猜想长宁背斜轴部的基底断层可能是倾向西南，倾角较缓，在40°左右，而这种倾角的推测也与本文以及易桂喜等(2019)的震源机制解节面的统计特征基本相符.

基于上述断层构造推测(走向130°，倾角40°)，利用长宁背斜轴部的应力场反演结果，分别将相对应力大小R取为0.15、0.3、0.5以及0.8，可以计算出该断层构造上的理论滑动角和相对剪应力大小(胡幸平, 2018; 万永革, 2020).这些计算结果(表2)表明，对于上述R的任一不同取值，计算出的断层面上的相对剪应力大小都超过0.9，表明本文推测的断层构造在区域应力场下具有很强的错动能力.另外，随着R取值的增加，计算出的断层错动的逆冲分量逐渐增加；对于R≤0.5的三个取值，其断层错动性质都是左旋走滑与逆冲比例相当，与实际长宁地震序列中强震震源机制解的总体特征相当吻合；而即使是R取为0.8，计算出的断层错动性质与实际情况相比也没有本质差别.因此，本文认为长宁地震的发生正是由于背斜轴部存在的基底断层在局部区域应力场作用下发生错动的结果.

表2 采用不同相对应力大小计算出的推断断层构造面上的 理论滑动角和相对剪应力大小Table 2 Theoretical rakes and relative shear stresses calculated by using different relative stress ratios on the inferred fault structure

6 结论

本文通过双差层析成像联合反演获取了长宁地区的小震定位和波速结构，在此基础上利用小震综合震源机制解方法反演了长宁地区地壳应力场的精细结构.结果显示出，长宁地区地壳应力场的最大主应力轴在整个研究区内基本都呈现为近水平状态，其方位虽然由北向南发生了一定角度的顺时针旋转，但也基本保持在近东西向；构造应力类型在长宁背斜轴部和南部向斜区之间存在显著差异，前者表现为逆冲型应力结构，而后者表现为走滑型应力结构.

通过上述应力场结果与gCAP反演得到的中强震震源机制解做力学对比分析，本文发现长宁背斜轴部与南部向斜区的局部应力结构与其各自对应区域内的中强震震源机制解的吻合度较高，但与非对应区域内的中强震震源机制解的吻合度较低，甚至在力学上是相抵触的.据此本文认为，区域应力场的局部改变是长宁地区中强震震源机制在十几公里范围内就发生显著变化的必要力学基础.而通过岩石力学估算分析，本文认为长宁地区的局部应力变化很可能主要是由岩石泊松比的横向差异造成的，但也有可能是多方面因素的共同产物.

最后，依据小震精定位和波速结构的剖面投影以及基于区域应力场的估算分析，本文推断了长宁地震的发震构造，认为长宁背斜轴部可能存在基底断层，控制着包括长宁MS6.0地震在内的中强地震的活动；而长宁地区地壳浅部的小震活动与中强震在深度上存在比较明显的分离，可能是由于上覆地层中的裂隙或者微小断层活动而产生的，而与轴部基底的较大规模的断层的关系较弱.

本文对长宁地区的地壳应力场及地震活动的研究分析侧重在空间差异上，而在不同的时间段，例如，在盐气开采的不同时期、以及长宁MS6.0地震等重要地震事件前后，是否也存在变化，由于资料所限并没有展开研究.然而，这种时间上的变化对于该地区地震成因机理研究和活动趋势预判都是相当重要的，值得继续深入研究.

致谢本文研究中所用的地震观测报告来自国家地震数据中心，波形数据来自中国地震局地球物理研究所国家测震台网数据备份中心，绘图通过GMT软件包(Wessel and Luis, 2017)完成；应急管理部国家自然灾害防治研究院的徐锡伟研究员和雷建设研究员对本文工作提供指导意见；匿名专家对本文提出修改意见.在此感谢上述机构和个人对本文工作的支持.