随着国内经济的快速发展，浅部资源的开发趋于饱满，使得挖掘深度越来越大，而在深部开采过程中遇到的巷道变形剧烈、瓦斯突出、岩爆与冲击地压频发等灾害问题与地应力增大息息相关^[1]。地应力是由自身重力和构造应力共同作用形成，随时间、空间变化的非稳定场^[2]，它会随着采掘工作面形成的超前应力与构造应力的叠加造成局部应力集中从而发生实时变化，同时复杂的地应力场会对工程建设造成严重的影响^[3-5]，尤其是煤层顶板高地应力区域。因此，准确、快速预测出地应力场的分布对开采资源和灾害预测等问题都显得尤为重要。

地应力测量的方法有很多，大体可以分为绝对地应力测量和相对地应力测量。绝对地应力测量主要就是测量地应力的大小和方向，例如水压致裂法、声发射法、钻孔崩落法、套芯应力解除法、应变恢复法等^[6-7]；相对地应力测量则是观测应力随时间变化的动态规律，例如刚性空心包体应力计、钻孔分量应力计、基于微震监测、基于光纤光栅等方法^[8-9]。

其中，微震监测技术具有实时、动态、连续监测煤矿井下动态活动的优点，数据记录中包含有地应力场变化过程的丰富信息^[10-11]。目前，通过层析成像技术探测地应力的方法，可分为主动层析成像和被动层析成像。采用主动源的主动层析成像技术已经得到了较多的应用^[12-13]，2012年王文书等^[14]基于地震CT技术初步建立了以地震波速度异常系数和波速梯度系数为主要因素的冲击地压危险性评价模型。同样在2012年，CAO Anye等^[15]通过澳大利亚南方煤矿LW704煤岩体现场爆破试验和微震监测，研究了地震波在煤岩体中传播过程中的微震效应，但基于微震的被动层析成像技术还没有得到足够的重视和发展。近年来，随着微震监测技术的快速发展以及实时监测的需求，被动层析成像技术逐渐得到关注^[11]。2016年蔡武等^[16]将微震实时监测与层析成像相结合来反演工作面的纵波速度，并与传统监测手段相对比，发现强矿震多发生于高波速区域。郭来功等^[17]在2018年采用微震走时成像技术，对采煤工作面的应力异常进行了监测。相对于主动层析成像，基于微震的被动层析成像监测范围广、安全性高、成本低、数据量大，可以满足长期连续监测的需要，是未来国内外监测冲击危险等地质灾害的一个发展趋势^[11]。

但是由于煤矿微震事件空间分布、观测系统布置、煤层位置等因素的不稳定性，被动层析成像技术的反演结果会受到严重影响。本次研究以三维地质模型为基础，采用快速三维层析成像算法对煤层顶板高地应力区域进行探测，以评估微震事件空间分布的不同以及高地应力区域与煤层顶板的间距对反演结果的影响，以期探测煤层顶板高地应力区域提供理论依据。

1   地震波速度与应力的关系

理论研究和大量的试验已证明，地震波的传播速度与岩石弹性常数、岩性、密度、孔隙率、应力等都有关系。同时，受地形、地质构造等影响，应力的离散度随深度的增加逐渐变大^[18-20]。近些年国内多位学者对此问题进行了实验，2012年巩思园等^[21]在单轴条件下对煤岩样的应力与纵波速度耦合关系进行了实验研究，结果表明应力与纵波速度之间存在幂函数关系，且岩样纵波速度对应力的敏感度要大于煤样纵波速度对应力的敏感度。薛雅荣等^[22]在2019年对突出型煤试样进行单轴加载破坏实验，结果表明纵波速度与突出煤样应力之间存在正相关关系。邓志刚等^[23]在2020年选取不同冲击倾向性煤岩测试其单轴加载和循环加载条件下各个应力阶段纵波速度的变化，发现无论是单轴加载还是循环加载，波速和应力同增同减，符合幂函数关系(图1)。因此，可以根据地震波速度大小，间接判断煤岩层所受应力状况，分析试验区域的应力异常。

图  1  煤岩应力−波速曲线(据文献[23]，修改)

Figure  1.  Stress vs. wave velocity curves of coals (Modified according to reference No.23)

下载: 全尺寸图片 幻灯片

2   快速三维层析反演原理

利用微震监测数据进行层析反演过程中，首先需要对微震事件进行定位，大致分为3个步骤：

(1) 微震事件的绝对位置定位。

(2) 通过双差算法对微震事件进行相对重定位。

(3) 在层析反演过程中同时校正微震事件位置参数和速度参数。

绝对定位微震事件的误差主要原因之一来源于初始速度模型不精确，如果真实速度模型和初始速度模型相差很大，那么微震事件位置相对于真实位置的偏移也可能很大。因此，为了减少初始速度模型不精确对定位造成的影响，算法将绝对定位中搜索的目标函数值(G)改成了一种由残差、距离和震相共同决定的特殊形式：

其中，

式中：M为所记录的地震事件总数；A为走时残差；${{\tau }}_{1}$和${{\tau }}_{2}$为预先定义的残差极限，如果残差大于${{\tau }}_{2}$，则认为这个值偏大或偏小，不进行计算，如果残差在${{\tau }}_{1}$和${{\tau }}_{2}$之间，则认为这个值可能处于速度异常区，进行计算；B为距离，路程较长的射线通常会有更大的残差，因此，其权重比路程较短的射线小；C_PS为震相，纵波速度扰动的影响小于横波，因此，纵波的残差具有更大的权重；d_i为射线i的路径长度。

在进行层析反演之前，首先要建立三维速度网格模型，通过设定平滑界面分割地层单元，利用三阶多项式进行层间网格节点的线性插值，并将速度场进行平滑处理。三维空间任意两点间的射线追踪快速修正方法是在起点设定一系列初始射线，并利用基于射线方程的单点射线的求解数值算法^[24]：

式中：P为三维空间X、Y、Z三个方向的慢度；v(X,Y,Z)为速度分布；dS为积分步长。穿过界面后的射线路径是根据Snell定律得到射线i的折射路径计算为：

反射路径计算为：

式中：${{P}}^{0}$为初始慢度，P为折射或反射后的慢度；N为曲面与射线交点处的法向量；v₁和v₂为地层界面两边的速度。

为满足实时监测的需求，通过以下3种方法来优化算法以提高层析反演的速度：(1) 在计算复杂三维模型时，根据当前点附近的速度变化特征来确定积分步长dS。如果射线通过速度变化很大的界面，就会自动增加dS的步长，以便下一个点尽可能靠近界面。(2) 在计算指定位置的剖面和平面时，节点根据每条垂直线上的射线密度来布设，如果没有射线，则不布设节点。(3) 为了避免节点分布对反演结果的影响，使用多个方向的网格(例如0、22°、45°、66°)进行反演，最后将其合并到一个模型中，以减少由网格方向所导致的节点分布不均匀而造成的误差。

3   模型验证及影响因素

建立X轴为3000 m、Y轴为3000 m、Z轴为1000 m的三维地质模型，在深度600~610 m之间插入厚度为10 m煤层，在深度400~600 m之间插入形状为圆柱体的高速异常体，代表高地应力区域(图2)，并设定模型中地层的物性参数(表1)。

图  2  三维地质模型

Figure  2.  3D geological model

下载: 全尺寸图片 幻灯片

表  1  三维地质模型物性参数

Table  1.  Physical parameters of the 3D geological model

岩性	纵波速度 vP/(m·s−1)	横波速度 vS/(m·s−1)
上方围岩	3000	1700
煤层	1800	1100
下方围岩	3000	1700
高地应力区	4200	2400

下载: 导出CSV 

| 显示表格

3.1   观测系统的建立

目前，煤矿微震数据采集方式可分为井中观测和地面观测^[25]。井中观测系统将震动传感器安置在巷道煤层或顶底板，虽然噪声等干扰小，采集数据分辨率高，但受煤层与围岩高波阻抗的影响，当地震波进入到煤层中后，大部分能量会被禁锢，使得围岩中信号能量变弱，难以识别^[26]。相比较而言，地面观测系统采集数据分辨率虽然较低，但覆盖范围广，成本低且安全性高，采集的数据也容易实时无线传输，便于全方位地实时监测地应力场变化。本次研究将使用地面观测系统对试验区进行监测。

地震层析成像反演的结果除了受初始模型的影响外，很大程度上还受数据完备程度的影响，如微震观测系统、射线覆盖密度以及模型网格尺寸。

将检波器以500 m的间距在3 km×3 km的平面内均匀布置，微震事件随机分布在所建立的三维地质模型中。通过前述射线追踪算法进行正演，得到水平面与垂直面的射线路径(图3)。从水平方向来看，更多的检波器和微震事件可以使射线路径有效的覆盖整个区域，但是从垂直面来看，边缘且较深的区域射线覆盖密度较为稀薄。因此，在观测系统的设计上，为能使结果更加的精确，应该使检波器边界的位置稍大于试验区的位置。

图  3  在500 m间距检波器下的射线覆盖密度

Figure  3.  Density of ray coverage under a geophone spacing of 500 m

下载: 全尺寸图片 幻灯片

在实际生产中，由于微震事件的位置有一定的聚集性，这可能会导致射线覆盖密度不均匀而影响反演结果。因此，通过设置一个简单的模型来对比微震事件聚集分布与均匀分布的反演结果，从图4中可以看出，均匀分布的反演效果较好，而聚集分布的反演结果中波速偏大且会产生低速异常区的假象。因此，在实际应用中，需要对微震事件进行筛选使其均匀分布。

图  4  微震事件聚集与均匀分布的反演结果对比

Figure  4.  Comparison of inversion results under uniform and concentrated distributions of microseismic events

下载: 全尺寸图片 幻灯片

3.2   层析算法稳定性测试

检测板测试，又称棋盘分辨率试验法，是一种常用的层析算法检测方法。它通过在初始速度模型上添加大小相同、正负相间的扰动值，形成棋盘格分布，然后进行层析正反演，最后对比分析反演结果和棋盘格模型，这种方法可以检验反演系统的稳定性。采用尺寸500 m×500 m的正方形网格，速度以3 km/s为基础，以20%上下浮动。由图5可知，反演出的结果几乎与检测板模型一致，但是四周由于射线覆盖密度较低而导致反演效果不够好。

图  5  检测板模型及反演结果

Figure  5.  Detection board model and inversion results

下载: 全尺寸图片 幻灯片

3.3   层析反演结果影响因素分析

由于煤层对地震波的屏蔽效应，使得地震波在地下的传播变得复杂，因此，不同位置的微震事件所产生的地震波在煤岩层中的传播规律会有一定的差异。图6中路线①和路线③都是没有经过煤层直接被检波器所接收的地震波，因此，信噪比会相对较高，而路线②和路线④在经过煤层时，由于地震波会在煤层中多次传播而导致地震波的能量衰减，从而被检波器接收后信号质量会有所降低。因此，微震事件空间分布的不同可能会对反演结果造成影响。

图  6  地震波在煤岩层中的传播路径

Figure  6.  Propagation patterns of seismic waves in coal seams and rocks

下载: 全尺寸图片 幻灯片

为准确探测煤层顶板高地应力区域，基于三维地质模型，将微震事件分别均匀分布在煤岩层中(图7a)、集中分布在煤层上方(图7b)以及集中分布在煤层下方(图7c)三种情况来评估对层析反演结果的影响。

图  7  微震事件空间位置及反演结果

Figure  7.  Spatial locations of microseismic events and inversion results

下载: 全尺寸图片 幻灯片

层析反演结果分成三部分来表现，水平截面、X-Z垂直截面和Y-Z垂直截面。为了突显速度异常区域以及增加算法的稳定性，在算法中加入了一些限制，它可以压制或者忽略非异常区域的波速值，而主要计算异常区域，因此，在反演结果中，非异常区域的速度值会接近于0。由图7a中可知，微震事件均匀分布时，煤层顶板高地应力区域较符合三维地质模型，但是波速偏低。图7b和图7c则是将微震事件空间位置分别集中分布在煤层上方与煤层下方，当微震事件集中分布在煤层上方时，几乎无法在反演结果中看到煤层的位置，而煤层顶板高地应力区域反演效果较好；当微震事件集中分布在煤层下方时，反演结果中高地应力区域位置和波速都被削弱的十分明显，而煤层区域的效果较好。微震事件分布的不同对反演结果的影响较大，一方面可能是因为所处位置的射线路径不能覆盖所有区域；另一方面可能是因为煤层对地震波能量的衰减。

冲击地压的能量大多来源于煤层顶板高地应力区域。通过分析图7a—图7c发现，无论微震事件位置如何分布，煤层顶板高地应力区域的位置和能量都会受到削弱，因此，考虑到高地应力区域与煤层的相对位置可能是导致反演结果精度不高的原因。在微震事件均匀分布的基础上，如图7d所示，将高地应力区域与煤层间距增加了200 m。由图7d可知，高地应力区反演结果确实有了较好的效果，位置与波速都较为精确。煤层对高地应力区反演结果的影响一方面可能是高波速区和低波速区附近都会产生波速梯度，互相产生了影响；另一方面可能是因为间距过近导致煤层对高地应力区域地震波能量的削弱。

4   试验数据分析与应用

4.1   微震观测系统及数据采集

试验区内总体构造形态走向为北北东—南北向西缓倾的单斜，在此基础上发育方向比较单一的宽缓褶曲，褶曲为北北东—南北向。主采二叠系下统3号煤，平均埋藏深度550~600 m，最大埋深深度大于750 m，煤层厚度5.00~7.25 m，平均6.00 m，煤层稳定。

本次研究微震观测系统采用地面观测方法，通过对试验区地质概况的了解，观测系统按照以下2个原则布置：(1) 要覆盖整个试验区区域，且观测系统的边界应大于试验区域。(2) 为了使试验区的微震事件记录误差相对较小，台站应均匀分布在试验区。依据上述基本原则和试验区地质状况，观测系统布置如图8所示。

图  8  试验区台站布置及微震事件分布

Figure  8.  Map showing the location of the test site and the distributions of microseismic stations and microseismic events

下载: 全尺寸图片 幻灯片

在试验区所布置的微震传感器的监测参数为：采样频率500 Hz，不断采集缓存；频率范围为0.1~150 Hz。采集的微震数据通过4G网络实时传回处理中心，图9展示了试验区某一时间段所接收到的20 s微震原始信号，可以看出起振点较为明显，其中N01，N10等代表检波器序号。在对原始信号进行频谱分析后，可以从图10中看出试验区微震信号能量主要分布在5~20 Hz之间，最后采用STA/LTA法与AIC(Akaike Information Criterion)法联合对信号初至自动、准确拾取。

图  9  试验区微震监测的原始信号(部分台站)

Figure  9.  Raw microseismic monitoring signals at the test site (partial stations)

下载: 全尺寸图片 幻灯片

图  10  微震监测原始信号频谱分析

Figure  10.  Frequency spectrum analysis of raw microseismic monitoring signals

下载: 全尺寸图片 幻灯片

微震事件定位分为3个步骤进行，首先通过一维模型实现微震事件的绝对定位，之后使用双差法对微震事件进行相对重定位，双差法是将一定空间范围内两个相邻的地震事件组对，根据组对事件在同一个检波器记录地震波到时差值来获得地震事件的相对位置^[27]。最后在层析反演过程中对微震事件参数和速度参数进行同时校正。为了获得较精确的反演结果，对所有的微震事件进行了筛选，手动选取了少量聚集区域的微震事件，大量选取了其余区域的微震事件，尽量使射线密度较均匀的覆盖整个试验区。从微震事件深度分布图11中可以看出微震事件大部分在深度300~600 m之间，处于煤层上方，依据在三维地质模型中的分析，反演结果会受到一定程度的影响。

图  11  微震事件深度分布

Figure  11.  Depth distributions of microseismic events

下载: 全尺寸图片 幻灯片

4.2   层析反演技术流程

通过对微震监测数据的分析处理为后续反演工作提供了较为准确可靠的前期资料，继而开展快速三维层析成像反演，其主要流程包括：(1) 首先为了保证三维模型的精确度，要对初始一维模型进行数据优化选择，尽可能地选择均匀分布的事件；(2) 定义好初始模型后，根据地震波到达的时间，采用复杂三维介质的射线追踪法重新定位微震事件，之后再使用双差分算法来同时校正微震事件参数和速度参数；(3) 层析成像的速度场参数化会将网格中的节点连接起来形成一个内部速度线性近似的四面体网格，有效的防止速度异常的混淆；(4) 利用三维模型中微震事件定位所得到的射线路径来计算一阶导数矩阵，之后使用LSQR方法对得到的矩阵进行反演。矩阵反演获得的波速异常和初始速度模型叠加，并用作下一次迭代的基本速度模型，当下一次的速度模型与上一次速度模型的差值小于初始设定的阈值时，迭代停止。

4.3   层析成像结果

在通过上述微震数据预处理方法获得微震事件位置和旅行时间后，使用快速三维层析成像算法对试验区进行波速反演，程序运行时间仅需要10 s。通过对试验区地质概况的了解，煤层深度在550~600 m范围，因此，着重对试验区550 m左右深度的平面进行层析成像反演，分别对400、500、550、600 m的4个水平截面进行反演，结果如图12所示。在400 m深度的反演结果中，高地应力区域较小，随着深度的增加，高地应力区域逐渐增大，直到深度550 m时高地应力区域最大，之后逐渐缩小，可以判断高地应力区域大致是中间宽上下窄的形状，并且在550 m深度的反演结果中有较大范围的低速异常区，可能是煤层导致。

图  12  试验区不同深度反演结果

Figure  12.  Inversion results under different depths at the test site

下载: 全尺寸图片 幻灯片

为了对高地应力区域有更深入的了解，使用三维快速层析成像算法对试验区X-Z垂直截面进行波速反演，结果如图13a所示，可以看出高地应力区域的位置基本可以与水平面的反演结果相对应，但是垂直截面反演结果中的波速有一定的衰减，且煤层的位置也不明显，可能是因为微震事件大多数处于煤层的上方，导致反演结果受到影响。试验区地质概况如图13b所示，可以看出高地应力区域处于试验区断层覆盖区域，初步判断可能是断层末端扰动所致，总体来说，采用快速三维层析成像算法反演试验区波速可以对试验区高地应力区域进行有效地判断。

图  13  试验区垂直截面反演结果和地质概况

Figure  13.  Inversion results of a vertical section and the geological conditions of the test site

下载: 全尺寸图片 幻灯片

5   结 论

a. 当微震事件均匀分布在煤岩层空间时，煤层顶板高地应力区域位置较准确，但波速偏小，会影响对地应力大小的判断；当微震事件集中分布在煤层上方或者下方，煤层顶板高地应力区域位置和波速都会受到一定程度的影响；当高地应力区域与煤层顶板的间距过小时，反演结果中位置和波速都会受到影响。

b. 根据地质条件合理地布置微震地面观测系统以及筛选微震事件使其均匀分布可以有效的提高射线覆盖密度，保证微震数据的完整性，为快速三维层析反演提供数据基础。

c. 试验区实际应用表明，快速三维层析成像算法和微震实时监测技术相结合可以有效地对煤层顶板高地应力区域进行探测，为监测冲击地压等动力灾害提供理论基础。