基于微震监测的页岩气开采区b值特征研究

胡隽 谭玉阳 张海江 曹俊兴 梁春涛 王权锋 徐彬

1)成都理工大学，油气藏地质及开发工程国家重点实验室，成都 610059 2)成都理工大学，数学地质四川省重点实验室，成都 610059 3)中国海洋大学，山东青岛 266100 4)中国科学技术大学，合肥 230026 5)成都理工大学，地球物理学院，成都 610059

0 引言

利用水力压裂技术改造储层来开采页岩气可能会诱发中小级别的地震，其原因是快速扩展的水力裂缝可能延伸到附近处于临界状态的断层并改变滑动面的应力状态，导致断层活化(Ellsworth，2013； Bao et al，2016； Hu et al，2018b； Tan et al，2020)。从全球范围看，虽然只有小比例页岩气开发地发生了疑似被水力压裂所诱发的地震，但是一些诱发地震的频率和震级足以引起政府和民众的重视(Foulger et al，2018)。加拿大阿尔伯塔地区弯曲湖(Crooked lake)附近的诱发地震序列与附近多个水力压裂在时间上高度吻合，震级范围主要在MW1.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年发生的数次ML>4.0和2次ML>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日的震级最大的MS4.8事件，采用Lei等(2017)通过CAP法得到的更可靠的MW4.56矩震级。压裂平台位置数据和生产信息来自相关研究(Lei et al，2017；Hu et al，2018a；Meng et al，2019；Tan et al，2020)。

图 1 研究区位置与地震事件投影(a)绘制了2009—2018年间发生的2级以上地震事件(数据来源于国家地震台网)，颜色对应地震发生时间，红框表示的研究区位于四川盆地南缘，放大后显示在(b)图中；(b)中的地震事件由布设的临时流动台站监测并精定位，圆圈大小对应不同震级

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

lgN(Mc)=a-b×Mc

(1)

式中，Mc表示地震最小完整震级，N(Mc)表示大于等于Mc的地震事件数，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)

(2)

(3)

式中，Mi为第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)。

图 2 TomoDD精定位后的地震事件分布(a)为研究区内部发生的深度10km以内的精定位地震事件平面分布图，不同颜色对应不同深度，AB和CD为两条剖面；(b)、(c)分别为AB剖面和CD剖面在方框内部提取的精定位事件的纵向投影，深度约为-10～0.5km，白色虚线为可能的断层走向；(d)给出了3个方向的精定位相对误差；(e)为精定位采用的地层速度模型(Lei et al，2017)

基于上述综合目录，本文计算了2017年1月至2018年6月间40km×40km范围内的b值。由于估算b值需要采用大于最小完整震级(Mc)的地震事件，因此首先采用最大曲率方法估计出Mc值。图3(a)中的灰色柱状数据指示出近正态分布的震级特征，对应的Mc值仅有0.51，误差为±0.04，其精度比区域台网监测事件所估计的Mc值(大于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年区域台网和国家台网地震目录，基于最小完整震级ML1.0计算出本区域的b值为0.9。Meng等(2019)利用区域临时流动台站在2015年2月至2017年12月期间监测的地震目录，基于最小完整震级ML1.1计算出本区域的b值为0.93。上述b值较本文的计算结果略低，可能的原因有两方面： 一是研究时间段不同，我们研究的是2017年2月至2018年6月的地震目录，晚于前人研究；二是通过近场监测获得了更多ML1.0以下的微震事件，使得估计出的b值更高。一般情况下，大多数水力压裂诱发地震的震级较小(小于ML3.0)，因此更多地震事件集中在小震级区域，拟合b值的直线斜率相应增大，估计出的b值更高。

图 3 根据精定位后的地震事件估计出的b值分布特征(a)为最大似然估计得到的b值、a值和最小完整震级(Mc)，采用bootstrap重采样方法通过1000次检验得到；(b)1000次重采样估计的b值

图 3(a)中给出的拟合直线(红色)与地震数据(蓝色)之间存在较大空隙，置信度不高。这里存在较明显的双线性特征(黑色虚线)，大于和小于ML2.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个地震快速增加的时间段也同步展现出了较大的振荡。

图 4 地震的震级-频度-时间分布及其对应的b值演化特征(a)监测期间的地震事件的震级-时间分布图，不同事件颜色对应不同发震时间，黑色曲线表示累积地震数量；(b)地震频度-时间分布图，红色柱体表示每周的地震数量；(c)监测期间研究区的b值和SCL值演化趋势

前人研究表明，郯庐断裂带南段、龙门山断裂带南段和鲜水河断裂带道孚段等大型断裂构造附近的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软件包(1)https://staff.aist.go.jp/xinglin-lei/，采用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)。

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

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

图 5 3km和10km深度的b值和Mc值的水平向分布特征(a)b值分布特征；(b)Mc值分布特征；(c)不同深度的频度分布

3 结论

本研究基于13个流动台站和国家地震台网固定台站监测到的12000余个地震事件，采用双差成像算法进行精定位，获得10499个重新定位后的地震事件。精定位后的地震剖面能够观察到条带状地震集中区，可能指示出一些产状较为陡峭的隐伏断层被活化。采用最大曲率方法估计出Mc值为 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值减小，Mc值增大，说明更深部的应力累积程度随着结晶基底的出现而增大。

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