0 引言

利用水力压裂技术改造储层来开采页岩气可能会诱发中小级别的地震，其原因是快速扩展的水力裂缝可能延伸到附近处于临界状态的断层并改变滑动面的应力状态，导致断层活化(Ellsworth，2013；Bao et al，2016；Hu et al，2018b；Tan et al，2020)。从全球范围看，虽然只有小比例页岩气开发地发生了疑似被水力压裂所诱发的地震，但是一些诱发地震的频率和震级足以引起政府和民众的重视(Foulger et al，2018)。加拿大阿尔伯塔地区弯曲湖(Crooked lake)附近的诱发地震序列与附近多个水力压裂在时间上高度吻合，震级范围主要在M_W1.7~3.9之间(Schultz et al，2015；Atkinson et al，2016)。同一地区靠近Fox Greek地区的Duvern地层中的水力压裂作业与附近发生的异常地震丛集也存在强时空关联，地震诱发强度与注水量之间存在吻合较好的线性关系(Schultz et al，2017、2018)。美国俄克拉荷马州和俄亥俄州在进行页岩水力压裂后发生了一系列2级以上的地震事件，诱发关系显著(Holland，2013；Darold et al，2014；Skoumal et al，2015)。中国四川省宜宾地区自2014年进行压裂开采以来，珙县、筠连县、长宁县附近的地震事件数量明显增加，诱发地震风险逐年上升，2017年至2019年发生的数次M_L＞4.0和2次M_L＞5.0地震，均可能与页岩气开发区的水力压裂活动有关(Lei et al，2017、2019；Meng et al，2019)。

上述现象表明，很可能因为水力压裂岩层产生的巨大弹塑性应力和高速扩展的流体压力使附近的既有断层被再次活化，显著增加区域地震风险。在此背景下，科学评估和识别地震风险更高的区域，把握风险演化的动态过程，成为一个极其重要的研究方向。一些学者分别采用单元细胞分割估算方法(Ghofrani et al，2016)、改进的贝叶斯网络法(Hincks et al，2018)、蒙特卡罗方法(Shabarchin et al，2017)、数据统计法(Schultz et al，2018)、模糊层次评价法(Hu et al，2018a)、机器学习算法(Pawley et al，2018)对加拿大西部地区、美国得克萨斯州、俄克拉荷马州、堪萨斯州和中国四川盆地南部的诱发地震风险进行评价，得到区域地震风险概率估计分布图。上述方法对诱发地震风险的估计具有科学价值，但主要采用数据统计方法，并未直接针对指针构造风险的地震学参数进行讨论。另一些相关研究计算了较大时空范围内的b值(Bao et al，2016；Schultz et al，2018；Lei et al，2017；Atkinson et al，2016)，发现诱发地震的b值存在双线性特征(Igonin et al，2018)，但仍然缺乏基于页岩气开发区微震监测数据的b值演化分析研究。

本研究针对中国四川省宜宾市长宁页岩气开发区，联合使用专门布设的13个流动台站监测获得的微震数据和国家地震台网的较大震级数据，对2017年2月至2018年7月的区域b值及其时空演化过程进行详细分析，讨论了b值振荡和演化的规律和可能成因，为页岩气开发区的地震风险评估提供了新思路。

1 资料与方法

位于中国四川省宜宾市珙县上罗镇的长宁页岩气开发区块，大地构造上位于青藏高原东缘和四川盆地南缘的盆山结合带，虽远离龙门山断裂带、华蓥山断裂带等大型构造单元，但是局部上位于长宁大背斜以南，隐伏断裂发育，构造环境仍较为复杂(Lei et al，2017)。

长宁区块设计了密集的水力压裂平台，自2014年起开展了总体上由南至北的规模性压裂开采，页岩气产量和地震数量一度同时攀高。2017年2月在长宁页岩气开发区内布设了13台高灵敏度检波器(图 1(b)中的白色倒三角)，获取了一年半时间的近场微地震监测数据。由于检波器限幅原因，临时台阵对于较大震级的地震波形记录并不完整。基于数据的完整性考虑，另收集了中国地震数据共享中心提供的国家地震台网固定台站记录到的远场地震目录。分别对近场和远场两类地震数据进行双差精定位，删去重复数据(基于精度考虑，小震级范围删去远场数据，保留近场数据)，得到本研究使用的综合地震目录。本研究的震级均采用里氏震级，只有发生在2017年5月4日的震级最大的M_S4.8事件，采用Lei等(2017)通过CAP法得到的更可靠的M_W4.56矩震级。压裂平台位置数据和生产信息来自相关研究(Lei et al，2017；Hu et al，2018a；Meng et al，2019；Tan et al，2020)。

b值是地震学中的重要参数，是描述地震频度和震级关系的系数，常用于表示区域内不同震级地震的相对分布特征，其概念由Gutenberg等(1955)提出

$ \lg N\left(M_{\mathrm{c}}\right)=a-b \times M_{\mathrm{c}} $

(1)

式中，M_c表示地震最小完整震级，N(M_c)表示大于等于M_c的地震事件数，a为历史地震活动水平的常数。

岩石实验结果和一些较大震级实例表明b值大小与环境应力累积程度成反比，低b值往往对应高地震风险(Scholz，1968；Fiedler，1974；Sahu et al，1994；Guha，1979；Molchan et al，1990、1999；Imoto，1991；Amitrano，2003；Schorlemmer et al，2005)，其是一个能够衡量区域应力累积状态的重要指标(Chan et al，2001)。本文采用更适用于分析小震活动性的最大似然法对b值进行估计(王辉等，2012)

$ b=\frac{\lg (\mathrm{e})}{\bar{M}-M_{\mathrm{c}}} $

(2)

式中，M_c为最小完整震级，也是本研究采用的起算震级，M表示平均震级，lg(e)值取0.4343。b值的标准偏差δb为(Shi et al，1982)

$ \delta b=2.3 b^{2} \sqrt{\frac{\sum\left(M_{i}-\bar{M}\right)^{2}}{n(n-1)}} $

(3)

式中，M_i为第i个地震事件的震级，n为地震事件总数。

本文就长宁页岩气开发区块在近场监测期间的地震b值进行时间和空间维度的综合分析，结合空间相关距离(SCL)、地震频度和丛集性特征探讨b值的演化规律和可能成因。

2 结果与讨论 2.1 精定位结果和b值

采用双差成像算法TomoDD(Zhang et al，2003、2006)和四川盆地南部的速度模型(Lei et al，2017)对13个流动台站监测到的10499个地震事件进行了精定位，获得10287个重新定位后的地震事件(图 2(a))。计算获得了较高的定位精度，X、Y、Z方向的相对误差均仅有0~20m(图 2(d))，其中3个方向相对定位误差的均值分别为7.10m、7.04m和8.69m，标准差分别为4.29m、4.12m和5.85m，可见Z方向具有相对更大的误差均值和方差。为提高流动台站监测3级以下的微震事件的效率，对监测设备进行了限幅，因此可能存在较大震级地震事件的遗漏。为了保证目录的完整性，将国家地震台网固定台站测定的3级以上地震事件纳入地震目录，形成用于本研究的综合目录。

图 2(a)中的精定位地震事件呈现出一定的区域集中性，可能与附近的页岩气开发平台的压裂开采作业有关。AB和CD两条剖面从EW向和SN向给出了垂向地震分布(图 2(b)、2(c))，剖面上能够观察到条带状地震集中区，预示着该区域可能存在一些产状较为陡峭的隐伏断层(图中的白色虚线)。这些疑似断层的产状基本均大于60°，与相关研究中的震源机制解吻合(Lei et al，2017；Tan et al，2020)。

基于上述综合目录，本文计算了2017年1月至2018年6月间40km×40km范围内的b值。由于估算b值需要采用大于最小完整震级(M_c)的地震事件，因此首先采用最大曲率方法估计出M_c值。图 3(a)中的灰色柱状数据指示出近正态分布的震级特征，对应的M_c值仅有0.51，误差为±0.04，其精度比区域台网监测事件所估计的M_c值(大于1)有较大提升，这一结果显然得益于13个流动台站的近场微震监测。然后，将研究区分割为0.05°×0.05°的等间距网格，设定每个网格的最少地震数目为50条，采用最大似然法(Aki，1965)计算每个网格的b值。由于网格尺度较小，绝大部分网格估计出的b值误差能够控制在0.01~0.05之间。基于上述方法，得出2017年2月至2018年7月间长宁页岩气开发区的b值为0.98，误差为±0.02，a值为4.31，误差为±0.01(图 3(a))。这一b值表示本文讨论的全时空平均水平。Lei等(2017)利用2014年至2016年区域台网和国家台网地震目录，基于最小完整震级M_L1.0计算出本区域的b值为0.9。Meng等(2019)利用区域临时流动台站在2015年2月至2017年12月期间监测的地震目录，基于最小完整震级M_L1.1计算出本区域的b值为0.93。上述b值较本文的计算结果略低，可能的原因有两方面：一是研究时间段不同，我们研究的是2017年2月至2018年6月的地震目录，晚于前人研究；二是通过近场监测获得了更多M_L1.0以下的微震事件，使得估计出的b值更高。一般情况下，大多数水力压裂诱发地震的震级较小(小于M_L3.0)，因此更多地震事件集中在小震级区域，拟合b值的直线斜率相应增大，估计出的b值更高。

图 3(a)中给出的拟合直线(红色)与地震数据(蓝色)之间存在较大空隙，置信度不高。这里存在较明显的双线性特征(黑色虚线)，大于和小于M_L2.5的b值拟合斜率明显不同，与加拿大西部盆地相关研究中类似(Igonin et al，2018)。这一现象的根本原因在于部分震级较大的诱发地震的余震震级较低，导致2.5~4.0级的中等震级地震事件数非常少。5个大于3.0级的地震(其中1个大于4.0级)的余震震级普遍较低，此特征不完全满足描述天然地震余震序列特征的大森定律，而在一定程度上更符合诱发地震的余震特征(Lei et al，2019)。这表明本区域地震的成因有很大一部分与外界应力的扰动有关，最大的怀疑对象即是附近的页岩气压裂开采。另外，采用的1000次Bootstrap重采样结果显示出b值的错位分布(图 3(b))，大部分b值为灰色部分，位于[0.95，1.00]区间，但是同样存在1.05和1.12左右的另2个b值丛集(图 3(b)中蓝色和红色)，这一现象表明b值存在时空差异性。

2.2 b值时空演化及分布特征

为了更全面和深入地探讨b值特征，首先从时间和空间的演化角度进行进一步分析。地震震级、累积数量、地震频率、b值和SCL值随时间演化的计算结果，如图 4所示。计算中采用固定样本数目50个事件进行扫描，即相邻统计点之间均发生了50次地震。图 4(a)中地震事件的颜色对应于不同的地震发生时间，地震事件更为集中的4个时间段用红色圆圈标识，其位置正好对应了累计地震数量曲线的快速增长部分，即曲线斜率突然上升的位置。图 4(b)中用红色柱体表示的地震频度也能提取出相同的地震丛集特征。图 4(c)中的b值在上述5个地震快速增加的时间段也同步展现出了较大的振荡。

前人研究表明，郯庐断裂带南段、龙门山断裂带南段和鲜水河断裂带道孚段等大型断裂构造附近的b值均体现出一定程度的低值特征(王辉等，2012；易桂喜等，2013；吴萍萍等，2015)。因此，b值的分布依赖于地壳介质的不均匀性和大型断裂带。b值的定义(公式(1))表明其值直接反映地震频度和震级之间的相关性，可以认为在远离大型构造活动带的、地震频度少而稳定的地区，b值一般不会出现较大波动。长宁页岩气开发区块距离最近的大型断裂华蓥山断裂有40余千米，在2014年11月进行大规模压裂开采页岩气之前，地震频度很低，属于稳定地区。然而，自该时间点以后区域地震数量和震级均呈现出截然不同的陡增态势(Lei et al，2017)。

图 4(a)、4(c)中的5个红圈标注出了b值振荡的主要部分，其刚好对应地震频率显著增高的5个时区(图 4(a)中的红圈和图 4(b)中的红色柱体)，并且能够观测到一些竖直条带状的地震丛集。上述现象说明某些时间点地震事件数的激增会使得b值出现明显波动。另外，借助Geotaos Map软件包^①，采用Single-link Cluster方法(Bruce et al，1989)对丛集地震的空间相关距离SCL值(Lei et al，2007)进行同步分析，其结果用蓝色线条描绘在图 4(c)中。总体上讲，所有40km×40km范围内的地震事件的SCL值均小于10km，直接阐明地震事件随时间发生存在时空丛集性的显著特点。同时，不难观察到每当密集地震丛集出现(图 4(a)、4(c)中红圈所示)，即b值开始振荡的时域，SCL值均呈现出较为稳定的低值特征(图 4(c)中绿色方框)。这说明5个地震丛集对应的空间相关距离仅有3~5km，基本符合水力压裂诱发地震的5km控制范围(Schultz et al，2017)。

① https://staff.aist.go.jp/xinglin-lei/

其次，进一步探讨地震事件激增的原因。结合相关研究数据(Lei et al，2019；Tan et al，2020)和遥感卫星图片，能够大致了解在b值振荡的对应时间段均存在正在进行的水力压裂作业活动，出现的密集地震丛集可能正是由这类人类活动所诱发。这一结论具有较好的灾害预警意义，因为b值的波动对应该地区的应力状态的变化，如果实时监控到大量的诱发微震事件伴随着b值的急剧变化，可能识别出继续压裂活动的地震高风险和低风险特征。值得注意的是，这并不意味着只要进行水力压裂就一定会触发大量微地震，因为很多压裂时间段内并未出现显著的地震丛集，例如2017年11月至2018年1月H7平台在压裂期间区域地震活动性较低。

另外，b值的演化过程体现出另一个规律，即4个较大震级的中型地震发生之前b值均存在下降的趋势。图 4(a)中用十字符号标注了监测时间段内区域内发生的4个M_L＞3.0的有感地震，虚线顶端标注了对应的具体震级。20世纪60年代的岩石物理实验(Scholz，1968；Mogi，1962)认为低b值对应高应力状态，是岩石破裂的前兆，也是强震孕育的重要条件，这一观点在诸如汶川8.0级地震、川滇地区若干6.0级以上地震、苏门答腊地震等事件中均得到了较好的验证(史海霞等，2018；王辉等，2012；Nanjo et al，2012)。本文讨论的大多数地震事件并不属于典型的天然地震事件，而更可能是与人类活动相关的诱发地震事件。震前低b值现象说明这些诱发事件仍然具有与天然地震类似的应力累积-释放特征，然而应力累积过程一定程度上被人为缩短了。图 4中表征4个3级以上地震的虚线左侧的灰色条带直接对应b值下降区间，其中震级最大的M_W4.56事件在发生后b值发生了显著的反弹回升，与大震级天然地震发生前后的b值特征变化模式最为吻合。

最后，继续分析b值的空间分布和不同深度的演化规律。图 5 (c)显示了地震事件数量随深度的分布，显然在2~4km深处达到最大值，更深和更浅的两侧逐渐下降，呈现出类正态分布特征。因此，选择了3km和10km为目标深度，并选取上、下2km范围内的地震事件扫描b值的空间分布。能够估计出b值在3km深度的投影范围大于10km，这与地震事件密度有关。中部区域的b值较高，外围则较低。并且随着深度增大，b值呈现出减小的趋势。这说明更深部的应力累积程度随着结晶基底的出现而增大。同时，扫描了M_c值的空间分布，总体趋势为中部有效震级较小，外围较大，深度越深M_c值越大。这一特征仍然与地震事件的空间分布密度相关，3km深处地震密集的中间区域能够获得更小的最小完整震级。

3 结论

本研究基于13个流动台站和国家地震台网固定台站监测到的12000余个地震事件，采用双差成像算法进行精定位，获得10499个重新定位后的地震事件。精定位后的地震剖面能够观察到条带状地震集中区，可能指示出一些产状较为陡峭的隐伏断层被活化。采用最大曲率方法估计出M_c值为0.51±0.04，b值为0.98±0.02，a值为4.31±0.01，b值的估计结果略高于已有研究。b值拟合直线存在双线性特征，重采样进一步发现其分布的差异性，指示出b值演化存在明显的时空差异特征。通过分析b值的时空演化过程，发现5个b值振荡的主要部分刚好对应地震频率显著增高的5个时区，且能观测到一些竖直条带状的地震丛集，说明地震事件数的激增会使得b值出现明显波动。同时，b值振荡位置的SCL值均呈现出较为稳定的低值特征，说明5个地震丛集仅约束于3~5km的空间相关距离，符合水力压裂诱发地震的距离控制范围。另外，发现4个较大震级的中型地震发生之前b值均呈现出下降的趋势，说明这些诱发事件仍然具有与天然地震类似的应力累积-释放特征，但其应力累积过程被人为缩短了。最后分析发现随着深度增大，b值减小，M_c值增大，说明更深部的应力累积程度随着结晶基底的出现而增大。

本文选取描述地震风险的重要指标b值，对压裂过程中的长宁页岩气开发区的风险演化特征进行了深入研究，得到了一些初步的结果。这些结论可能具有较好的灾害预警意义，由于b值的波动对应该地区的应力状态的变化，如果实时监控到大量的诱发微震事件伴随着b值的急剧变化，可能识别出继续压裂活动的地震高风险和低风险特征。当然，由于监测时间、研究手段和生产数据的限制，本文的部分观点仅为科学推测，更加准确的结论仍有待进一步讨论。