工业开采诱发地震的应力降研究现状和主要科学议题

2022-02-03姜丛尹欣欣蒋长胜翟鸿宇王骞张延保来贵娟尹凤玲

中国地震 2022年4期

姜丛尹欣欣蒋长胜翟鸿宇王骞张延保来贵娟尹凤玲

1)中国地震局地球物理研究所，北京 100081 2)甘肃省地震局，兰州 730000

0 引言

在地壳浅部采用水力压裂等方式进行的页岩气和干热岩等非常规能源开采、废水回注以及抽取方式进行的矿山开采，均可造成地壳应力状态的扰动并产生诱发地震，给这些工业活动的顺利实施以及地表社区与重大基础设施的地震安全带来较大挑战。应力降Δσ是描述地震震源特征的重要参量，被定义为地震破裂前后断层上的平均应力差(Hanks，1979)。由于应力降可被用来考察高频地震动特征(Boore，1983)、区分所在地区的构造背景(Allmann et al，2009；Boyd et al，2017)、根据其自相似性给出关于地震危险性的认识(Abercrombie，1995；Oth et al，2010)，因此其在天然地震的地球动力学研究和减轻地震灾害风险研究中广受青睐。在工业开采诱发地震活动中，应力降还可用来了解水力压裂等工业活动与诱发地震震源性质之间因果关系(Yu et al，2020)，以及作为重要参量来揭示断层活化现象(Chen et al，2020)，这使其成为目前诱发地震研究中的重点内容之一。尤其是开采活动中布设大量的密集地震观测网络，为更深层次利用应力降提供了新机遇(Shearer et al，2019)。

(1)

Δσ=cM0/S/L

(2)

从而实现了直接由地震波测量应力降的过程。由上述定义可见：①受到震源模型假设的影响，Δσ的具体表达形式依赖于震源模型，例如半径为R的圆形断层上的应力降为Δσ=7M0/16R3，长为L宽为w的矩形断层上的应力降为Δσ=2M0/πw2L；②由于断层面积S≈L2，应力降Δσ与1/L3相关，其测定结果的不确定性很大(Allmann et al，2007)；③由于应力降的概念基于静态位错理论，仅考虑了断层的相对移动，没有考虑断层的破裂，因此通常所说的静态应力降要比实际的应力降小。

鉴于上述关于应力降Δσ的定义和物理含义，以及其在新型工业活动和诱发地震灾害风险管控中的重要价值，有必要系统全面地总结目前关于诱发地震应力降研究的新进展、技术发展趋势和已获得的科学认识。而在人为工业活动中的这些进展和认识，对天然地震的研究也同样有重要启发意义。本文尝试对近20年的相关研究进行总结，以期为新型工业活动伴随的诱发地震的减灾研究以及天然地震的研究提供参考。

1 工业开采诱发地震应力降的计算

包括应力降在内的震源参数的测定，需要分离或同时测定几何扩散、介质衰减、场地响应等波形信号的耦合信息。新型工业活动诱发地震多为震级较小的微震事件，波形信号的低信噪比给测定工作带来较大难度。目前已发展的测定方法有如下三类：

(1)谱拟合方法是计算诱发地震的应力降的最为直接的方法。一般首先利用多锥形窗方法(Percival et al，1993)分别对三分量的地动速度谱进行平滑处理，并利用Boatwright点源模型拟合地动速度谱(Boatwright，1980)，即

(3)

(2)单事件的经验格林函数方法(EGF)是诱发地震应力降计算的重要方法，其优点是可一定程度提高拐角频率fc和标量地震矩M0的测量精度(Kwiatek et al，2015)。EGF方法的原理是，对拟求解震源参数的地震，找到位置接近的小震级事件(震级差ΔM>1.0)作为经验格林函数，利用谱阈的反卷积方法，消除待求解事件与EGF事件相同的场地响应、传播路径衰减项，从而获得震源谱。其中拟求解事件i与EGF事件j的地动速度谱的比值为

(4)

对式(4)一般可采用最小二乘法拟合理论谱与实际的观测比值谱的L2范数最小值来获得参数结果(Holmgren et al，2019；Wang et al，2020)。在实际应用中，能否找到严格吻合EGF定义的小地震是制约经验格林函数方法的重要因素(Huang et al，2019)，此外对EGF方法计算结果分辨率的限度也有较多讨论(Abercrombie，2015；Shearer et al，2019)。

μ 为重力参数，ω =[00 0]探测器在轨道坐标系下相对于惯性坐标系的角速度，θ0为探测器的真近点角。设rc=[x y z]T探测器在轨道坐标系与小行星的相对位置，a=[axayaz]T是探测器各轴的控制加速度，是d=[dxdydz]T外部环境的扰动。

(3)利用谱叠加和广义反演方法(Castro et al，1990；Shearer et al，2006；Oth et al，2011)获得更具普遍意义的经验格林函数(EGF)是计算应力降的重要方法。例如，Warren等(2002)以及Shearer等(2006)发展了一种P波谱叠加方法，可以将震源谱与几何扩散、传播路径的衰减、场地响应分离出来。该方法通过将同一台站大量的地震记录进行叠加，迭代后消除几何扩散、路径衰减等差异项，再利用震级分组的震源项与理论震源谱的拟合获得修正后的经验格林函数，从而更为完整地去除几何扩散、路径衰减和场地响应，获得每个地震的震源谱，并利用Brune圆盘震源模型计算应力降等震源参数。Trugman等(2017)在该方法的基础上，发展了非参数重采样方法来估计震源参数的不确定性。为提高谱叠加和广义反演方法中的震源谱参数反演的可靠性，Picozzi等(2017)还发展了一种基于数据驱动策略的无参数计算方法，其核心是设计了震源谱lg10(S(f))对各个震源参数的无量纲的偏导数作为频域相依的敏感度函数，作为反演的权重函数。

2 工业开采诱发地震应力降的平均值

目前在应力降的数值上已取得一定共识。工业开采诱发地震的应力降通常较震源深度更深的天然地震的应力降值偏小(Fehler et al，1991；Abercrombie et al，1993；Hua et al，2013；Hough，2014)。已有的对全球天然地震的研究表明，应力降至少在三个数量级上变化(例如，Shearer et al，2006；Allmann et al，2007、2009)，数值通常在1～100MPa之间(Abercrombie，1995)，且符合正态分布(Allmann et al，2009)。而对工业开采诱发地震，应力降数值则通常在0.01～10MPa之间。本文罗列了一些工业开采区诱发地震应力降的测定结果，以便进行比较研究。

在北美有大量的诱发地震应力降测定案例。例如，Wu等(2019)对2015年美国俄克拉荷马州加斯里(Guthrie)MW4.0诱发地震应力降的测定结果为3.4MPa，这一数值接近Wu等(2018)测定的俄克拉荷马州其他诱发地震应力降的中值(约2MPa)。Kwiatek等(2015)采用Madariaga震源模型观测到美国加利福尼亚州盖尔瑟斯地热田Prati-9井附近的一组诱发地震，震级分布在MW2.0～3.8之间，这些地震的平均应力降约为7MPa，但明显低于Viegas等(2011)采用经验格林函数方法对同一水力压裂注入井测定的应力降中值(28MPa)。Jeong等(2020)对美国得克萨斯州沃思堡盆地的90次诱发地震(ML2.0～3.9)的应力降测定结果显示，平均应力降约为4.46MPa。在加拿大不列颠哥伦比亚省东北部的蒙特尼区块，通过P波估算的应力降中值为0.08MPa，对数标准差为0.8(以lg为单位)，同时利用S波估算的应力降中值为0.03MPa，对数标准差为0.7(Yu et al，2020)。在加拿大西部沉积盆地(WCSB)，Holmgren等(2019)利用S波测定了87次震级为M2.3～4.2的诱发地震的应力降，发现中小震级事件的平均应力降为 7.5±0.5MPa，但测定结果在0.2～98MPa之间变化，且有一个极端事件的应力降测定值为370MPa。

在欧洲和亚洲，部分页岩气和干热岩开采诱发地震的应力降得到了测定和研究。其中，Kettlety等(2020)测定了英国兰开夏郡普雷斯顿新路PNR-1z井2018年发生的震级为ML1.5的诱发地震，平均应力降值为1.0MPa。对于瑞士巴塞尔地热田附近2006年发生的MW3.2±0.1诱发地震，以及圣加仑地区2013年发生的MW3.4±0.1诱发地震事件，Edwards等(2015)测得的平均应力降分别为3.5MPa和3.0MPa。在芬兰赫尔辛基，Kwiatek等(2019)对OTN-3井附近震级为MW2.0的诱发地震测得的平均应力降为1.6MPa。在亚洲的韩国浦项地区，Song等(2019)对PH02井附近的诱发地震计算获得的平均应力降约为2MPa，与Chai等(2020)得到的1.92MPa结果一致。

表1给出了这些诱发地震应力降测定结果。由上述案例可见，对于多数震级在1.5～4.0的诱发地震，平均应力降分布在0.03～7MPa之间，且多数主要分布在1.0～7MPa之间。目前认为诱发地震应力降的数值相比天然地震更小的现象，主要是由诱发地震发生在地壳更浅部造成的。

表1 全球部分地区工业开采诱发地震的平均应力降测定结果

3 影响应力降测定数值的因素

诱发地震的发生机制与天然地震明显不同，至少包括注入流体影响岩石孔隙压力、注入的较冷流体与高温岩石相互作用引起热弹性应变和裂缝滑动(Rawal et al，2014；尹欣欣等，2021)等较为独特的机制，因此影响应力降数值的因素可能也有所不同。下面分别从主观因素造成测定结果的不可靠性以及客观因素造成测定结果的复杂性做出简要总结。

震源模型和常数参数的设置，是应力降计算的重要假设，可直接影响应力降的计算结果。例如，Yu等(2020)使用圆盘震源模型(Brune，1970)和取值S波的k=0.37，相比于使用Madariaga模型(Madariaga，1976)以及取值k=0.21，其计算得到的应力降较后者小约5.5MPa。Yamada等(2007)发现使用圆盘震源模型(Brune，1970)要比使用Madariaga模型(Madariaga，1976)获得的应力降小，并分析认为使用Madariaga模型更为合理。破裂速度的取值也是应力降计算中的重要影响因素，较慢的破裂速度将导致震源半径的估计值较小、应力降的估计值较大。而破裂速度的取值通常带有明显的主观性，例如在Madariaga模型中常取剪切波速 3.9km/s的90%，这在采用中值计算平均应力降的研究中影响更为明显(Allmann et al，2009)。地震波传播介质的品质因子Q值的计算和取值，也会显著影响应力降测定的可靠性(Ko et al，2012)。Q值对岩石的力学性质非常敏感，其空间变化可以作为水力压裂影响局部应力状态的指标(Yu et al，2020)，但当对衰减校正不充分时常常可得到较低的应力降测定结果(Tomic et al，2009)。由于上述震源模型和参数设置的人为主观因素影响，Holmgren等(2019)、Yu等(2020)和Wang等(2020)建议，应力降的绝对数值只有在使用相同的震源模型和参数设置的条件下才可以进行比较。

一些客观因素也显著影响了应力降的分布差异，包括构造环境的差异、存在的水量、破裂速度的变化、断层面正应力的变化等(Tomic et al，2009)。这些因素导致应力降数值时空分布出现差异，另一方面也使得根据应力降数值反推地下施工区或储层应力状态变得更为复杂。这些影响因素至少包括：①水热型地热田和增强型地热田(干热岩)的高温特征，对地壳的摩擦强度等力学行为有重要作用，例如Oth(2013)对日本地区的应力降的研究表明，低应力降的数值在空间上与高热流区有关；②诱发地震在深度分布上的差别以及流体介入地震破裂过程的程度也能明显影响应力降的数值，这是由于震源深度较浅和流体注入时，可能会降低有效应力，从而导致剪切应力更容易克服莫尔-库仑破坏准则进而产生诱发地震(Goertz-Allmann et al，2013)；③诱发地震的应力降在空间上具有方向性差异，且对中等大小的地震(Lui et al，2019)和微震(Folesky et al，2016)均可观测到。而当地震记录的台站覆盖较为稀疏时，应力降测量方式不当则可能带来结果的较大不确定性，例如仅从断层的前向或后向测量，可能会将拐角频率fc的测定结果带来2倍误差，并导致应力降结果达到8倍偏差(Holmgren et al，2019)。

4 与应力降有关的重要科学议题

4.1 诱发地震应力降的空间相关性

探讨诱发地震应力降的空间相关性问题具有重要的科学价值。这是由于断层带和裂隙上可能存在应力的不均匀性(Chen et al，2020；姜丛等，2022)，应力降的空间变化可以揭示断层强度的非均质性，以及反映断层带的水文地质性质，这有助于理解诱发地震的发生机制。

应力降是否具有“近低、远高”分布现象近年来在诱发地震研究中成为热点问题。Yu等(2020)在对加拿大不列颠哥伦比亚省东北部的蒙特尼区块水力压裂井周边3km范围内的484次地震(M-1.0～3.0)的研究中发现，近井端的应力降(0.1～1MPa)要低于远井端的应力降(1～10MPa)，且约小一个数量级。类似的现象也在加拿大阿尔伯塔省克鲁克湖(Crooked Lake)附近的注水诱发地震(Clerc et al，2016)，以及瑞士巴塞尔和萨尔瓦多柏林地区的增强型地热系统(EGS)开采区诱发地震研究(Goertz-Allmann et al，2011；Kwiatek et al，2014)中被证实。上述例子中应力降数值与井的距离相关的现象，在统计学意义上也得到了证明(Yu et al，2020)。Yu等(2020)认为，应力降的这种空间分布特点可能是由于井筒附近水力压裂受损岩石的破裂速度较慢，裂隙密度高或孔隙压力较高，阻止了井筒附近地壳岩石存储和释放较大应力，而远端地震事件的应力扰动受水力压裂的影响则相对有限或几乎可以忽略。

应力降是否具有空间分布的恒定性也是诱发地震研究的关注点之一。在对加拿大不列颠哥伦比亚省东北部的蒙特尼区块水力压裂诱发地震的研究中，Yu等(2020)还发现无论是发生在距离井筒近处还是远处的地震，应力降的数值随着时间变化而大致不变。Kwiatek等(2015)和Picozzi等(2017)分别采用不同的应力降计算方法，在对美国Geysers地热田Prati-9井诱发地震的分析中，也表示并未观察到静态应力降随时间的明显变化。这种空间分布的恒定性，可能表明具有相似流变特性的岩石中的动态破裂通常具有相对恒定的应力降(Yu et al，2020)。还有一些研究发现，地震序列的应力降具有随时间变化的特点，主要表现在主震后应力降随时间增加，在余震序列的后期观察到的应力降总体较高(Sumy et al，2017；Chen et al，2020；Goertz-Allmann et al，2011)，Goertz-Allmann等(2011)认为这是孔隙压力过高导致的。但这些序列的应力降时序变化是否在同一位置、余震序列活动的短期调整变化是否能影响空间分布的恒定性，还缺乏更深入的研究。

4.2 诱发地震应力降与震源机制类型的相依性

即使在天然构造地震研究中，应力降与震源机制类型的依赖性也是重要的研究课题。根据安德森断层理论，逆断层作用的剪应力最高，正断层作用的剪应力最低，走滑型断层的剪应力介于前两者之间，尽管应力降和绝对应力之间的关系尚不清楚，但传统观点认为应力降和绝对应力在与震源机制类型的关系上相同(McGarr，1984；McGarr et al，2002)。但也正因为应力降和绝对应力之间的关系不明确，一些新的观测现象在逐渐揭示其复杂性。例如Allmann等(2009)在对地震应力降的全球变化的研究中观察到了应力降中值对震源机制类型有依赖性，其中走滑事件的应力降最高，中值约为10MPa，而正断层和逆断层事件的应力降较低，中值为2～3MPa，并未显示逆冲型地震具有更高的应力降。

对于诱发地震，Kwiatek等(2015)对美国盖尔瑟斯地热田Prati-9井水力压裂诱发地震的研究表明，尽管正断层型地震的应力降略低于走滑型地震的应力降，但由于走滑型地震在储层深部分布的数量更多，因此这种应力降与震源机制类型的相关性难于验证。研究这一问题的难点在于影响应力降数值大小的因素，例如在一些应力降的计算中假设断层破裂速度是恒定的，而Allmann等(2009)采用可变的断层破裂速度来计算应力降时，可以解释走滑型地震具有较高的应力降可能是由更快的断层破裂速度造成的。

4.3 诱发地震应力降的非自相似性

应力降是否具有自相似性在物理上具有特殊意义，其实质是大地震和小地震是否具有相似的动态破裂过程。在实际应用中，同一地区小地震的应力降空间分布是否与中等以上地震存在相关性、能否用平均的应力降代表这一地区的总体应力降水平，直接关系到是否可用数量更多的小震级的诱发地震来预测未来破坏性诱发地震的应力降，并减少概率性地震危险性评估和地震预警中的地面运动不确定性(Hardebeck，2020)。

对应力降的自相似性/非自相似性研究，最重要的研究方式是考察应力降是否与震级或标量地震矩存在相关性。Holmgren等(2019)在对加拿大西部沉积盆地87次M2.3～4.2级诱发地震应力降的测定中，发现应力降随震级增大呈上升趋势。Novakovic等(2018)和Chen等(2020)在对美国俄克拉荷马州以及Picozzi等(2017)对加利福尼亚州盖尔瑟斯地热田进行诱发地震研究中，也均发现较大震级的诱发地震具有相对更高的应力降。Chai等(2020)对韩国浦项诱发地震序列的应力降研究中也获得了同样结论。其他案例还包括Mayeda等(2005)和Venkataraman等(2006)等的研究。概括起来，目前支持工业开采诱发地震应力降具有自相似性的案例总体上较少，仅有Yamada等(2007)对南非金矿开采诱发的20次微震的静态应力降的研究等少数案例，此外，尽管Picozzi等(2017)在美国加州盖尔瑟斯地热田诱发地震的研究中观察到了大地震群体具有自相似趋势，但对单个地震序列的研究往往发现强烈的非自相似行为。

4.4 诱发地震应力降的深度依赖性

在天然地震的应力降研究中，许多案例显示出应力降与震源深度存在依赖性(Hardebeck et al，2009；Allmann et al，2007；Baltay et al，2017)，并认为这可能与断层的摩擦强度沿深度的变化有关。工业开采诱发地震多发生在地下3～10km深度的地壳浅部，表1给出的大量案例也表明诱发地震的应力降总体相比天然地震偏小，因此目前普遍持有的观点是造成诱发地震应力降比天然地震偏小的原因也在于震源深度的影响(Satoh，2006；Hough，2014；Yenier et al，2015；Zhang et al，2016；Atkinson et al，2017)。

但对无论是诱发地震还是天然地震是否存在应力降大小的深度依赖性目前还有较大争议，例如Wang等(2020)认为在计算应力降时假设剪切波速恒定可能导致人为的深度依赖性，特别是发生在地壳浅部5km以上的地震，Allmann等(2007)也持有类似观点。Holmgren等(2019)进一步解释认为，由于应力降与速度的立方呈正相关，也认为由于剪切波速随深度增加，应力降的估计对假定的震源深度很敏感。在具体的应力降测定研究中，Allmann等(2009)对全球范围的天然地震研究表明，尽管存在较大的分散性，但平均应力降估计值随深度的变化很小。进一步从不同震源机制类型的角度分析，Allmann等(2009)还发现正断层事件和走滑事件的平均应力降均没有深度依赖性，而逆断层作用下平均应力降随深度的增加与Bilek等(1998、1999)在俯冲带地震中发现的结果一致，可以用深度相关的刚度变化来解释。左可桢等(2021)使用谱比法分析四川长宁地区微震震源参数之间的相互关系以及应力降的时空分布特征，也得到应力降总体与震源深度关系并不明显的结果。

尽管正反两种观点的争议较大，表1给出的大量案例至少给出了诱发地震的应力降数值分布范围更窄、总体比天然地震偏小的初步认识。由于应力降测定结果的主观可靠性、客观影响因素较为复杂，这也使得对这一问题的争论和研究在未来仍具有开放性。

5 结论和讨论

鉴于应力降在研究工业开采诱发地震的震源特征、高频地震动特征、区分构造背景和地震危险性分析等诸多领域的重要价值，本文系统地总结了近20年来相关方向的研究进展。分别从计算所用理论方法、获得的平均应力降数值、影响应力降数值的主客观因素、围绕应力降的重要科学议题4个方面做了总结。获得的主要认识如下：

(1)在工业开采诱发地震的理论计算方法上，针对诱发地震微震较多、信噪比较低、实时快速计算需求明显等特点，目前主要发展了三类方法，包括直接采用谱拟合方法求取震源参数、利用小震级事件作为经验格林函数(EGF)提高拐角频率fc和标量地震矩M0测量精度以及利用谱叠加和广义反演方法对多台站多事件同时计算震源参数。不同方法尚存在明显的不足之处，这也限制了对应力降数值本身和分布规律探索的过高期待。

(2)对于工业开采诱发地震的应力降平均值，本文收集了相关案例，结果显示多数震级在1.5～4.0的诱发地震平均应力降分布在0.03～7MPa之间，且多数主要分布在1.0～7MPa之间，应力降数值相比天然地震分布的范围更窄、总体偏小。对这一现象的主要解释集中在诱发地震往往发生在地壳更浅部，以及流体更多地参与诱发地震发生过程等因素上。

(3)影响工业开采诱发地震测定数值的因素，在测定方法的可靠性上，主要受到采用相对简单的圆盘模型等震源模型假定，以及对波速、断层破裂速度、地壳介质品质因子Q值的具体取值差异等因素影响；此外，构造环境的差异、存在的水量、破裂速度的变化、断层面正应力的变化等客观因素也影响了应力降的数值本身。

(4)围绕工业开采诱发地震应力降的科学问题，总结了应力降是否具有“近低、远高”分布特征以及是否存在空间分布恒定性等空间分布依赖性问题，以及诱发地震的应力降是否与震源机制类型存在相依性、是否存在相似性或非自相似性、是否存在深度依赖性等4个问题。其中对应力降存在距离依赖性、存在非自相似性有一定程度的共识，但鉴于应力降测定可靠性和多种客观因素的影响，这些科学问题未来均有较大的开放性、需要持续探索。

工业开采区通常布设更为密集、覆盖程度更好的现代化地震台网，这为诱发地震获得更为可靠的包括应力降在内的震源参数提供了便利。但通过本文的总结可见，受观测技术和测定理论方法的局限，诱发地震的应力降测定的可靠性问题并未获得实质性突破。对诱发地震的应力降很难精确测量似乎成了为数不多的共识性问题(Shearer et al，2019)，Stein等(2009)甚至认为当前的静态应力降研究仅反映了观测谱或震源谱的一般特征、无法直接反映震源的物理特性，甚至会误导对真实的断层活动的认识。由于应力降是由拐角频率通过假设的震源模型和参数间接计算得到的，对此Holmgren等(2019)和Shearer等(2019)建议应该转而关注拐角频率fc，这是相比应力降更基本、有用且更接近实际数据的震源参数。