中国地震  2021, Vol. 37 Issue (2): 332-348
BSN矿山微震监测技术及其应用
张达1,2,3, 戴锐1,2,3, 曾志毅4,5, 冀虎1,2,3, 石雅倩1,2,3, 常莹1,2,3, 韩鹏4,5     
1. 矿冶科技集团有限公司矿山工程研究设计所, 北京 102628;
2. 中国-南非矿产资源开发利用联合研究中心, 北京 102628;
3. 中国-南非矿产资源可持续开发利用"一带一路"联合实验室, 北京 102628;
4. 南方科技大学, 深圳市深远海油气勘探技术重点实验室, 广东深圳 518055;
5. 南方科技大学, 地球与空间科学系, 广东深圳 518055
摘要:深井高应力、高岩压诱发的冒顶、片帮、塌方、岩爆等地压灾害严重影响了矿山安全生产。随着高精度、高采样率传感器的研发和计算机技术的飞速发展,微震高精度实时监测技术不断完善,微震监测已成为矿山地压灾害监测预警的重要手段。通过设计、研发完全自主知识产权的BSN(BGRIMM Seismic Network)微震监测系统,实现了微震信号数据采集、多通道时钟同步、噪声压制、到时自动拾取、震源定位等功能,并结合统计地震学的相关方法,初步实现了基于微震活动性分析的岩体危险性评价。本文主要从矿山微震监测数据采集、预处理、定位、活动性分析等模块介绍BSN微震监测系统全流程,并对监测系统在广西某矿山的实际应用案例进行分析。
关键词深井开采    矿山微震监测    微震去噪    震源定位    地震活动性分析    
Technology and Application of BSN Microseismic Monitoring in Mines
Zhang Da1,2,3, Dai Rui1,2,3, Zeng Zhiyi4,5, Ji Hu1,2,3, Shi Yaqian1,2,3, Chang Ying1,2,3, Han Peng4,5     
1. Institute of Mining Engineering, BGRIMM Technology Group Co., Ltd., Beijing 102628, China;
2. China-South Africa Joint Research Center for Exploitation and Utilization of Mineral Resources, Beijing 102628, China;
3. China-South Africa BRI Joint Laboratory for Sustainable Development and Utilization of Mineral Resources, Beijing 102628, China;
4. Shenzhen Key Laboratory of Deep Offshore Oil and Gas Exploration Technology, Southern University of Science and Technology, Shenzhen 518055, Guangdong, China;
5. Department of Earth and Space Sciences, Southern University of Science and Technology, Shenzhen 518055, Guangdong, China
Abstract: The mining safety is seriously effected by the ground disasters such as roof caving, banding, collapse, and rock burst induced by high stress in deep well. With the development of high-precision, high-sampling rate sensors and computer technology, real-time microseismic monitoring with high-precision was achieved and became an important means of ground pressure disaster monitoring and early warning. The BSN(BGRIMM seismic network)microseismic monitoring system with completely independent intellectual property rights was developed by BGRIMM Technology Group. The BSN system has realized microseismic signal acquisition, multi-channel time synchronization, denoising, automatic onset time picking, and source location. By employing the seismic activity analysis methods, the BSN system could provide preliminary results of the risk assessment of the rock mass. In this paper we introduce the BSN microseismic monitoring products and show the details of the whole process including data collection, processing, positioning, and risk analysis. An actual application case of the monitoring system in a mine in Guangxi is further demonstrated.
Key words: Deep mining     Microseismic monitoring     Denoising     Microseismic positioning     Seismic activity analysis    
0 引言

我国矿产资源正逐步进入深井开采阶段,深部高应力、高岩压诱发的冒顶、片帮、塌方、岩爆等地压灾害严重影响了矿山安全生产(Liu et al,2016Salvoni et al,2016)。根据国家矿山安全监察局最新公布的统计数据,2020年全国共发生矿山事故434起、死亡573人。因此,矿山安全监测具有重大社会需求和实用价值。

随着高精度、高采样率传感器的研发和计算机技术的飞速发展,矿山微震高精度定位、实时监测已成为可能。随着微震定位技术的不断完善,微震监测已能够实现矿山危险源的时空定位及稳定性评估,近年来广泛应用于深井矿山安全监测,已成为当前行业高度认可的地压灾害监测预警手段(Gibowicz et al,1998田玥等,2001李庶林等,2005姜福兴等,2006唐礼忠等,2006冯夏庭等,2013王璐琛等2016Palgunadi et al,2020)。在矿山微震监测领域,国外起步较早,南非矿震技术研究院(Institute of Mine Seismology)和加拿大工程地震组织ESG(Engineering Seismology Group)作为该领域研究的先行者,积累了丰富的地压灾害防控经验,其开发的微震监测系统占据了绝大多数矿山市场。国内相关研究起步较晚,成熟的监测系统相对较少。针对这一现状,矿冶科技集团有限公司自2014年开始逐步研发了完全自主知识产权的BSN(BGRIMM Seismic Network)微震监测系统,实现了多通道时钟同步、噪声压制、到时自动拾取、震源定位等功能。并结合统计地震学相关方法,初步实现了基于微震活动性分析的岩体危险性评价。本文主要从矿山微震监测数据采集、预处理、定位、活动性分析等模块介绍BSN微震监测系统全流程,并对监测系统在广西某矿的实际应用案例进行分析。

1 BSN微震监测系统 1.1 微震数据采集

微震传感器作为采集微震波形信号的换能单元,将被测信号由震动信号转换为可被微震监测系统识别的电信号。由于爆破震动信号能量在频带375~500Hz表现较强,而矿山岩体破裂信号(微震)能量主要集中于低频带0~250Hz(朱权洁等,2012),低频成分的监测能够更好地记录到矿山微震信号,了解矿区内岩石破裂的区域,有利于矿山的安全开采。因此,矿冶科技集团有限公司自主研发了具有较好低频响应特性、可实现硬/软岩下微震信号高保真拾取的矿山单轴微震监测检波器(图 1(a)),其幅频响应曲线如图 2(b)所示,可以看到频响曲线在8~1000Hz频带内较为稳定,具有良好的响应特征。仪器在10Hz频点的自噪音小于6×10-7m/s · Hz-1/2。检波器支持安装姿态的一体化自测量、检波器编号自辨识等功能。安装姿态一体化自测量具有传感器标识、指标、属性、类别等关键参数的在线识别功能,同时,为减少人工测量检波器姿态的繁杂工作及测量误差,在传感器内部集成姿态测量单元,可实现倾角等姿态的动态获取和上传。

图 1 BSN微震检波器外观(a)及其幅频响应曲线(归算到输入振幅为1m/s)(b)

图 2 BSN微震数据采集站

BSN微震数据采集站(图 2) 作为硬件平台的核心数据管控装置,主要由数采装置、时间同步装置、不间断供电装置等构成,可以分别实现数据的高精度、高同步、不间断的并行采集、处理、存储及传输等功能,同时也支持有效数据预触发、可靠数据存储、网管型程控管理等功能。图 3为安装在实际矿山中的BSN微震数据采集站,具体参数如表 1所示。

图 3 BSN微震数据采集站现场安装

表 1 检波器与数据采集站主要性能参数

BSN-Trigger微震监测系统数据采集软件主要实现硬件与软件的配置、数据解析、数据存储和显示等功能。该软件由多个管理逻辑通讯协同运行的软件模块组成,整体通讯由连接到数据采集控制器的通讯子系统来实现,具有较高的数据传输效率。综合数据管理工具支持微震事件数据查询,数据可存储于任何一台具有相应数据库配置的本地或网络计算机,以便进行微震数据的本地及远程处理和分析。

1.2 微震数据处理

矿山微震信号的特点是频率范围广、信噪比弱,这使得对信号的处理变得比较困难,采矿现场的背景噪声、工频干扰、人工活动等信号会严重干扰有效的微震事件波形信息,极大地降低微震波形的信噪比,导致P波到时拾取难度增大、可靠性降低。矿山微震监测系统每天监测到海量的微震事件,虽然可以通过数据处理员人工操作直观地拾取P波到时,但人工手动拾取波形到时易受数据处理员的经验和主观性影响,且耗时耗力。因此,矿山微震弱信号降噪与到时自动拾取对微震实时监测技术具有重要意义。

BSN-Processor微震数据处理软件由矿冶科技集团有限公司自主研发,运行在远程网络终端,主要实现现场数据的滤波、震相到时识取、震相分类、波速校准、震源定位、震源参数计算等数据处理功能,并将处理及分析结果通过网络推送至安全监测服务平台,为矿山安全监管员提供可靠的实时矿山安全开采信息,辅助其准确做出矿山安全开采决策。

1.2.1 微震信号降噪

在地震数据处理中,带通滤波是一种常用方法,但当信号频带和噪音频带有重合时去噪效果不佳。此外,带通滤波方法须预先知道有效的滤波频带,易产生波形失真,影响到时拾取的精度。为了提升去噪效果,一些学者开发出了一系列去噪方法,如主成分分析(Principal Component Analysis,PCA)(Hagen,1982Bekara et al,2009朱鹤文等,2020)、奇异值分解(Singular Value Decomposition,SVD)和经验模态分解(Empirical Mode Decomposition,EMD)(Ursin et al,1985刘志鹏等,2012)等滤波方法,这些滤波方法均在时间域(一维)进行。也有学者使用时频变换(Time Frequency Transform,TFT),将一维时间域的地震波信号变换到二维时间-频率域内,在时-频图谱上较为清楚地识别出微震信号位置,进而进行信号处理和去噪。时频域小波系数阈值去噪方法因原理简单(Donoho et al,19941995),去噪效果显著而受到广泛关注(Chang et al,2000Mousavi et al,20162017程浩等,2018)。阈值法利用数据本身来决定哪些系数是重要的,具有较好的自适应性和稳健性。Donoho等(1994)研究表明,对于方差为σ2的高斯白噪声来说,小波系数的通用阈值可表达为

$ \lambda=\sigma \sqrt{2 \lg N} $ (1)

其中,N为信号长度。通用阈值方法假设小波系数大部分是由噪声造成的,而信号主要体现在能量较大的小波系数中。因此,硬、软阈值去噪函数可分别表示为(曹静杰等,2015)

$ W_{H}^{\prime}=\left\{\begin{array}{ll} |W|, & |W| \geqslant \lambda \\ 0, & |W| <\lambda \end{array}\right. $ (2)
$ W_{S}^{\prime}=\left\{\begin{array}{ll} \operatorname{sgn}(W)(|W|-\lambda), & |W| \geqslant \lambda \\ 0, & |W| <\lambda \end{array}\right. $ (3)

式中,W表示小波系数,$ W_{H}^{\prime}$$W_{S}^{\prime} $表示经过硬、软阈值函数处理后的小波系数,sgn(W)表示符号函数。当小波系数小于阈值λ时,则去除该小波系数,否则将其保留。其中硬阈值函数是直接保留大于阈值λ的小波系数,而软阈值是对大于阈值λ的小波系数进行平滑减弱,使滤波后的信号更加平滑。

本微震监测系统中实现了带通滤波方法和基于同步挤压小波变换(Daubechies et al,19962011)的小波系数硬、软阈值滤波技术,滤波效果如图 4所示,由图可见原始微震信号信噪比较弱,初至信息被噪音严重干扰。利用不同滤波方法去除背景噪音后,小波系数软阈值去噪方法的效果最好,信噪比由原始波形的3.4837提升到最高的11.5499(信噪比为以初至时刻为中心向前、向后各延拓一个短时窗内波形的振幅均方根之比),并且能够较好地去除与信号同频带内的背景噪音,有利于初至波识别和精确拾取。

图 4 不同滤波方法的微震信号去噪效果对比 (a)原始微震信号记录及其时频图;(b)带通滤波后微震信号记录及其时频图;(c)硬阈值滤波后微震信号记录及其时频图;(d)软阈值滤波后微震信号记录及其时频图
1.2.2 微震到时自动拾取

微震信号到时拾取的准确性直接影响震源定位的精度(刘翰林等,2017)。到时拾取主要分为以下几类:第一类是长短时窗比法(STA/LTA)(Allen,1978刘晗等,2014),以长、短2个时窗内波形振幅、能量或者波形包络函数等特征函数的比值作为初至波拾取的判断准则;第二类是假设噪声分量和微地震信号分量具有不同的统计性质,利用高阶统计量(峰度与偏度)、熵等特征函数实现微震信号的自动识别和拾取(Saragiotis et al,2002郭铁龙等,2017);第三类运用互相关理论对目标波形与模版波形进行匹配,得到相对初至到时差(Senkaya et al,2014),然后与模板波形的绝对到时相加,得到目标波形的初至到时。上述方法虽然已经有效地实现了到时自动拾取,但拾取精度仍不能满足微地震/声发射震源定位需要。因此,一些学者综合利用微震信号与环境噪音的多种特征差异,提高微震信号识别和拾取效果。例如,宋维琪等(2011)利用小波分解与Akaike信息准则相结合,提高了微震信号的拾取精度;谭玉阳等(2016)综合利用地震信号与噪音在能量、偏振和统计方面存在的差异,有效地提高了初至波拾取方法的抗噪能力和拾取精度。

在本微震监测系统中,已实现了基于STA/LTA-AIC和基于峰度、偏度(K-S)两种自动初至波拾取方法。其中,STA/LTA算法公式为

$ \begin{array}{c} \frac{\mathrm{STA}}{\mathrm{LTA}}(i)=\frac{\frac{1}{W_{\mathrm{S}}} \sum _{j=i-W_{\mathrm{S}}}^{i} \mathrm{CF}(j)}{\frac{1}{W_{\mathrm{L}}} \sum _{j=i-W_{\mathrm{L}}}^{i} \mathrm{CF}(j)} \end{array} $ (4)
$ \mathrm{CF}(j)=Y^{2}(j) $ (5)

式中,i为采样时刻;WSWL分别为短时窗和长时窗的长度;Y为波形振幅值;CF(j)为在j时刻的微震信号的特征函数值,表征微震数据振幅、能量及频率的变化。该算法的主要影响因素包括STA和LTA的时窗长度、触发阈值及特征函数的选取。在本微震监测系统中,选择振幅平方作为对微震事件自动检测和拾取的特征函数。通常,短时窗长度设置为微震信号主周期的2~3倍,长时窗为短时窗长度的5~10倍(刘晗等,2014)。经过反复测试对比分析,本微震系统的短时窗长度为100个采样点,长时窗长度设置为8倍的短时窗长度。

AIC是建立在熵概念的基础上,从信息论和极大似然原理推导出来,用来估计模型的复杂度和衡量此模型拟合优良性的一种标准。由于噪声分量和微地震信号分量的统计特性差别较大,导致在微震信号和噪声交界点处,两种信号的拟合度最差,此交界点的AIC值最小,即为微地震事件的初至时刻。给定一个长度为N的微震数据Y,假设第i点为噪声分量和微地震信号分量的最佳分界处,则信号在第i点被分成两段,对应的AIC值可近似表示为(张唤兰等,2013)

$ \operatorname{AIC}(i)=i \lg ({var}(Y[1, i]))+(N-i-1) \lg ({var}(Y[i+1, N]))+C $ (6)

其中,var表示数据的方差,Y[i+1,N]为信号在窗长为N-(i+1)内的振幅,C为常数。

偏度、峰度是度量微震记录高斯性特征的重要方法之一,基于微震信号为非高斯分布、背景噪音记录为高斯分布的假设(Saragiotis et al,2002),采用信号振幅的偏度(Skewness)与峰度(Kurtosis)特征作为微震信号初至波拾取的特征函数,利用偏度、峰度极值点位置信息拾取初至波。其计算公式如下

$ \begin{array}{c} S(i)=\frac{1}{W} \sum\limits_{j=i}^{i+W}\left[\frac{Y(j)-\bar{Y}}{\sigma}\right]^{3} \end{array} $ (7)
$ K(i)=\frac{1}{W} \sum\limits_{j=i}^{i+W}\left[\frac{Y(j)-\bar{Y}}{\sigma}\right]^{4}-3 $ (8)

其中,W为时窗长度,$ \bar{Y}$σ为时窗内振幅均值和标准差。时窗长度W不宜过大或过小,应尽量刚好包含微震信号。时窗太长导致敏感性降低,反之使得一些高振幅噪音被误识,降低拾取精度。经过反复测试,时窗长度W为400个采样点时,拾取结果较为理想。

利用上述STA/LTA-AIC和K-S两种自动拾取方法来处理单道的实际微震信号,处理结果如图 5图 6所示。长短时窗(STA/LTA)比值法对信号的识别具有较好的效果,但是初至波拾取精度不够高,选择准确的阈值较为困难,而AIC拾取方法具有较高的局部拾取精度,因此选择联合STA/LTA与AIC作为本微震监测系统中的一种拾取方法。另一种基于峰度、偏度(K-S)的自动拾取方法同样具有较高的拾取精度,当微震记录的峰度、偏度同时满足判别阈值时,可实现微震信号的精确拾取。从图 6拾取结果来看,两种不同自动拾取方法的拾取结果均具有更高的精度,满足定位精度要求,而且可以节省人工成本,提高效率。为了更好地展现自动拾取方法的精度,将多个微震信号的自动拾取结果与人工拾取结果进行对比,得到两种自动拾取方法的结果与人工拾取结果的误差统计分布(自动拾取结果减去人工拾取结果),如图 7所示,由图可见两种拾取方法的总体拾取误差较小,主要分布在-20~0个采样点区间,标准差约为10个采样点。考虑到采样率为6000Hz,对应的到时拾取误差为0~3.3ms,标准差约为1.6ms。自动拾取方法的初至波结果相对于人工拾取结果总体上偏小(时间提前),该结果是合理的,因为人工拾取一般是根据振幅极大值点位置,然后在该位置附近拾取波,因此可能会造成初至波拾取结果偏大。人工拾取到的初至波偏差会随着信噪比的减小而增大。

图 5 基于STA/LTA-AIC和K-S特征函数的自动拾取方法

图 6 不同拾取方法的初至波拾取结果

图 7 两种自动拾取方法的初至波拾取结果与人工拾取结果的误差统计分布
1.2.3 微震事件定位

矿山微震具有频带宽、发生频次高等特点。利用微震监测技术可以分析岩体损伤破裂过程中产生的微震信号,并对微震事件实现定位。目前,微震监测系统普遍采用单事件绝对定位方法来高效率地进行定位,其原理和方法可以借鉴传统地震学的研究(田玥等,2002)。双差定位方法可以提高地震事件相对位置的定位精度(Waldhauser et al,2000),前人在此基础上陆续开发出了精度更高的震源定位方法及能够同时反演精细速度结构和地震位置的双差层析成像方法(Zhang et al,2003Guo et al,2017)。

本微震监测系统中的定位模块基于Geiger(1912)的定位方法,该方法已经广泛应用于微震监测中(林峰等,2010毛庆辉等,2015)。假定微震事件发生在空间位置(xyz)和时间t0,地震波经过矿山介质的传播,被空间位置为(xiyizi)的微震监测系统传感器i记录下到达时刻ti,定义走时残差为ζi=ti-(Ti+t0)。其中,Ti为地震波从发震位置到传感器i的旅行时间,是关于发震位置(xyz)的函数。

对于单一速度模型,假定地震波传播速度为v,则

$ T_{i}=\frac{\sqrt{\left(x-x_{i}\right)^{2}+\left(y-y_{i}\right)^{2}+\left(z-z_{i}\right)^{2}}}{v} $ (9)

实际观测区域的介质模型并非单一速度模型,传感器记录到的地震波到达时刻也存在误差,因此存在地震波走时残差,利用观测到时与计算到时的残差(ri=(Ti+t0)-ti)最小的方法来估计发震位置(xyz)和发震时刻t0

对残差ri进行一阶泰勒展开

$ r_{i}=\left(T_{i}+t_{0}\right)-t_{i}=d t+\left(\frac{\partial T_{i}}{\partial x} d x+\frac{\partial T_{i}}{\partial y} d y+\frac{\partial T_{i}}{\partial z} d z\right) $ (10)

其中,i=1,2,…,NN为观测到时数目;dxdydz为震源坐标校正量;dt为发震时刻校正量;$ \frac{\partial T_{i}}{\partial x}、\frac{\partial T_{i}}{\partial y}、\frac{\partial T_{i}}{\partial z}$为走时空间偏导数。将式(10)写为矢量形式

$ \boldsymbol{A} \boldsymbol{\Delta} \boldsymbol{X}=\boldsymbol{r} $ (11)

其中,系数矩阵AN×4矩阵,rN×1矩阵,震源校正量$\boldsymbol{\Delta} \boldsymbol{X} $为4×1矩阵,则

$ \boldsymbol{A}=\left[\begin{array}{cccc} 1 & \frac{\partial T_{1}}{\partial x} & \frac{\partial T_{1}}{\partial y} & \frac{\partial T_{1}}{\partial z} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \frac{\partial T_{N}}{\partial x} & \frac{\partial T_{N}}{\partial y} & \frac{\partial T_{N}}{\partial z} \end{array}\right] $ (12)
$ \boldsymbol{r}=\left[\begin{array}{c} r_{1} \\ \vdots \\ r_{N} \end{array}\right] $ (13)
$ \boldsymbol{\Delta} \boldsymbol{X}=\left[\begin{array}{l} d t \\ d x \\ d y \\ d z \end{array}\right] $ (14)

以最早接收到地震波信号的传感器坐标和观测时刻作为发震参数的起始值,带入式(11)进行迭代计算,直至残差r小于阈值或达到最大迭代次数则停止迭代,得到震源定位结果。

1.2.4 地震活动性分析

微震活动性分析是微震监测研究的重要环节,其可以更快速地核定矿区地震危险性,方便矿山管理者选择相应的对策措施,在保证矿山正常生产的情况下,确保采矿工人的安全。地震活动性分析方法包括地震活动的时间和空间分析,较为简单的分析方法是利用震源位置的空间分布了解矿区内岩石破裂的区域,以及用微地震事件位置表示裂缝的空间形态。为了进一步利用微震事件位置信息,许多学者通过聚类方法实现了对微震事件的划分(Wang et al,2019刘德彪等,2019),减少了人工主观划分微震事件的不确定性,发现潜在的微震集群,从而有效地分析微震事件分布特征和活动规律。震级-频度关系(lgN=a-bM)是统计地震学研究中最基础的规律之一(Gutenberg et al,1994),其中a值表征地震活动水平,b值表征地震频次在不同震级上分配的比例。b值是判定地震危险性区域的重要指标,其可以为地震危险性方面的研究提供依据(史海霞等,2018Xie et al,2019),目前被广泛应用于地震预测和地震危险性研究中。高b值对应低应力,低b值则对应高应力,因此b值可以间接判断地下应力的状态和变化,进而分析矿区地震活动性(张楚旋等,2016)。传染性余震序列(ETAS)模型(Ogata et al,198819892001余娜等,2018)也是分析地震活动性的重要方法之一,该模型通过分析余震的衰减速率和背景地震活动性变化来研究地震风险程度,其中背景地震活动性与地区的应力状态和应力加载速率相关。通过ETAS模型可以客观、定量地分析矿区地震活动性,提升矿山安全监测水平。

微震震级受辐射花样的方位性、传播路径衰减、传感器响应等影响,会存在一定误差。而BSN微震监测系统采用的是单轴传感器,近场的辐射花样方向性响应引起的误差不能校正,所以目前微震监测系统中的震级计算误差较大,而统计地震学中的G-R关系和ETAS模型两种地震活动性分析模型建立在地震目录完备的前提下,对震级准确性要求较高,故G-R关系和ETAS模型还未应用于本微震监测系统中。因此,目前BSN微震监测系统中仅开发了利用震源位置时空分布对地震活动性进行分析的功能模块。

1.3 微震数据分析与发布

BSN-3D Analysis微震数据分析软件运行于远程网络终端,其主要在微震监测数据处理的基础上,通过采用二维图表或三维图形,完成震级/数量-时间统计、时空分布、震级分布、能量-地震矩、累积视体积-能量指数等专业分析,从而结合矿山采矿作业,进一步完成数据深入分析及可视化。

基于矿山安全监测分析云服务平台,建立数据报表、分析报告、预警信息等模板库,实现异地数据同步与前台应用开发,最终达到分析—统计—报表—报告—预警过程的在线发布,使用户能够快速获得有关事件的参数及其位置信息。该服务平台提供便捷的缩放、平移和旋转的功能,实现微震数据三维发布和交互,微震监测安全分析结果的综合展示如图 8所示。通过云服务平台,用户可远程、实时监测矿山的开采情况。

图 8 矿山安全监测分析云服务平台
2 工程应用

目前,矿冶科技集团有限公司自主研发的BSN微震监测系统已成功在广西某矿山实现工程应用,该矿矿体分布广、矿体薄,且属急倾斜矿体,岩石结构性差,矿石开采过后留下了大量的采空区,存在重大安全隐患。虽然前期矿山也采用了钢筋混凝土进行支护,但现场观察到这些支护结构大部分均有变形,最大的裂缝宽度已超过4cm,并有错位情况出现,地压问题突出,如图 9所示。故对该矿山安装了BSN微震监测系统,实时监测矿山微震情况。

图 9 矿山地压安全隐患点

采用微震监测技术实时监测采空区及开采区域岩体破裂发生的微震事件,研究微震事件在时间、空间上的分布特征,分析监测区域地压活动规律,为矿山安全、合理生产提供参考。微震监测系统在75m中段、35m中段与-5m中段各穿脉采空区或采场周边共布置了40个微震监测点,其中75m中段共布置16个监测点,分别布设在0#穿脉、2#穿脉、3#穿脉、4#穿脉、5#穿脉、7#穿脉;35m中段共布置8个监测点,分别布设在0#穿脉、1#穿脉、3#穿脉;-5m中段共布置16个监测点,分别布设在0#穿脉、2#穿脉、3#穿脉、4#穿脉、5#穿脉、7#穿脉,具体位置如图 10所示。

图 10 各中段微震监测点布置图 三角形为检波器;灰色线段为采掘巷道

该矿2020年11月的微震事件定位结果如图 11所示,其中X为EW向,Y为SN向,Z为垂向。从XZYZ面可以看出微震事件主要集中在-5m中段,与该时间段内的矿山开采深度一致。本系统采用盖格定位方法,对定位结果的误差估计由地震波到时的观测值与计算值之间的差,即时间残差来表示;而定位结果在空间坐标XYZ上的误差由时间残差和走时方程系数矩阵的协方差矩阵来估计。微震事件定位误差统计分布结果如图 12所示,可以看到XYZ方向的定位误差均值大致为2.5m,均方根误差在2.5~3.8m之间,总体误差分布较小,表明定位结果较为精确。通过XY面的微震事件空间分布特征可以看到微震活动大致分布在3个区域(图 11中红色虚线椭圆所示),微震事件较为集中,由此可初步判断这3个区域的地震活动性较高。

图 11 微震事件时空分布

图 12 定位结果在XYZ方向上的误差统计分布

图 11中的微震事件位置信息进行进一步的地震活动性分析,将矿山监测区域离散为10m×10m×5m的空间网格(水平向为10m,垂向为5m),求取每个网格内的地震累计数目,结果如图 13所示。从图中可以看出矿区内大致有3处地震活动性较为活跃的区域(图 13中红色虚线椭圆),微震事件较为集中,表明这3个区域的地震危险程度较高。图 14为矿山安全员在11月28日巡查时发现的3处围岩破坏及支护明显变形区域(图 13中五角星所示位置),与定位结果主要集中在-5m中段的情况一致,其中2处地压变化区域与微震事件聚集性分布区域较为接近(图 13中A和B所示的五角星),验证了地震活动性分析结果的正确性。图 15展示了矿区微震事件空间累积数目随时间(周)的变化,颜色越深,代表该区域的微震事件数目越多,呈集中现象。研究表明矿山地震的可能机制可分为两类,一类与开采活动直接相关,由高应力岩石突然破裂失稳造成的较小震级矿震;另一类与开采活动间接相关,由较大尺度采掘空间或整个矿区尺度引起的区域应力重分布所致的较大震级矿震(李铁等,2006王泽伟等,2017)。矿山开采过程中地压是随时间动态变化的,随着矿山开采时间的推移,微震活动性较为活跃的区域逐渐凸显出来,微震事件空间累积数目的时空变化较好地反映了矿山潜在危险区域的时空演化特征。实际应用案例表明,微震事件整体时空特征及演化趋势对现场采矿及安全生产有一定的指导作用,便于一线采矿作业部门及时调整工作计划。基于BSN微震监测结果及危险性分析,已建议矿山技术人员加强对应力集中区和微震事件聚集区域开展现场巡查,严格要求现场工作人员注意施工安全,及时加强支护等保护措施。

图 13 矿区微震事件的空间累积数目 黑色三角形为检波器;红色虚线椭圆为地震活动性强的区域;实心五角星为矿山安全员实地巡查发现地压显现较为明显的区域

图 14 实际矿山围岩变形破坏情况 (a)、(b)、(c)分别对应图 13中实心五角星A、B、C所标注的矿区位置

图 15 矿区微震事件的空间累积数目随时间的变化 (a)2020年11月1日—2020年11月07日;(b)2020年11月1日—2020年11月14日;(c)2020年11月1日—2020年11月21日;(d)2020年11月1日—2020年11月30日;黑色三角形为检波器
3 结论和讨论

BSN微震监测技术及产品已实现微震信号数据采集、多通道时间同步、信号噪声压制、到时自动拾取、震源定位、地震活动性分析等功能,并在实际矿山安全监测中得到了良好的应用。本系统在矿山稳定运行,在复杂工况环境下具备良好的抗干扰能力,且具备可拓展的功能以及可持续的后处理分析服务能力。基于云服务模式,该系统实现了数据的采集、传输、分析、地压管控的一体化,并结合现场-远程的联合分析,强化对矿山安全生产的指导作用,形成高效监测、自动处理、深入分析、智能预警及灾害防控等综合性地压监测解决方案,具备较高的自动化分析处理能力。

基于小波系数阈值的去噪方法虽然具有较好的去噪效果,但是去噪后的波形仍存在干扰初至波拾取的少量背景噪音,因此需要研发去噪效果更好的信号滤波方法。本微震监测系统的定位方法对初至波拾取精度和速度模型准确度敏感性较高,而矿山岩体存在各向异性、传播介质波速不均匀等,会带来基于均匀波速模型的微震定位误差。目前,已有学者发展了震源位置与速度结构联合反演策略,并在矿山监测中取得了良好的应用效果(Zhang et al,2015Qian et al,2018Wang et al,2018)。因此,后续需要对BSN系统升级,将目前较为先进的定位方法进行集成,以提高矿山微地震事件的定位精度和矿山的精细速度结构。此外,目前系统采用的是基于到时拾取的定位方法,今后拟将无需初至波拾取的基于偏移成像的定位方法植入BSN系统,进一步提升定位精度和稳定性(Kao et al,2004Li et al,2020曾志毅等,2020)。同时,还需要在矿山布设更多三轴微震传感器,并对微震震级进行计算和标定,为地震活动性的进一步分析提供依据。

参考文献
曹静杰、王本锋, 2015, 基于一种改进凸集投影方法的地震数据同时插值和去噪, 地球物理学报, 58(8): 2935-2947.
程浩、袁月、王恩德等, 2018, 基于小波变换的自适应阈值微震信号去噪研究, 东北大学学报(自然科学版), 39(9): 1332-1336.
冯夏庭、陈炳瑞、张传庆等, 2013, 岩爆孕育过程的机制、预警与动态调控, 北京: 科学出版社.
Gibowicz S J, Kijko A, 1998. 矿山地震学引论. 修济刚, 译. 北京: 地震出版社.
郭铁龙、张雪梅、邹立晔, 2017, STA/LTA—AIC算法对地震P波震相拾取稳定性影响, 地震地磁观测与研究, 38(3): 13-17. DOI:10.3969/j.issn.1003-3246.2017.03.003
韩骏、姚令侃, 2015, 基于地震活动性参数b值的地应力评估方法研究, 铁道标准设计, 59(7): 36-39, 127.
姜福兴、杨淑华、成云海等, 2006, 煤矿冲击地压的微地震监测研究, 地球物理学报, 49(5): 1511-1516. DOI:10.3321/j.issn:0001-5733.2006.05.032
李纪汉、刘晓红、郝晋升等, 1987, 岩石摩擦滑动的声发射b值, 西北地震学报, 9(4): 36-38.
李庶林、尹贤刚、郑文达等, 2005, 凡口铅锌矿多通道微震监测系统及其应用研究, 岩石力学与工程学报, 24(12): 2048-2053. DOI:10.3321/j.issn:1000-6915.2005.12.008
李铁、蔡美峰、蔡明, 2006, 采矿诱发地震分类的探讨, 岩石力学与工程学报, 25(增刊Ⅱ): 3679-3686.
林峰、李庶林、薛云亮等, 2010, 基于不同初值的微震源定位方法, 岩石力学与工程学报, 29(5): 996-1002.
刘德彪、李夕兵、李响等, 2019, 基于LOF的K-means聚类方法及其在微震监测中的应用, 中国安全生产科学技术, 15(6): 81-87.
刘晗、张建中, 2014, 微震信号自动检测的STA/LTA算法及其改进分析, 地球物理学进展, 29(4): 1708-1714.
刘翰林、吴庆举, 2017, 地震自动识别及震相自动拾取方法研究进展, 地球物理学进展, 32(3): 1000-1007.
刘志鹏、赵伟、陈小宏等, 2012, 局部频率域SVD压制随机噪声方法, 石油地球物理勘探, 47(2): 202-206.
毛庆辉、王彦春、王鹏等, 2015, 改进的微震事件反演重定位方法及其应用, 石油物探, 54(3): 359-366. DOI:10.3969/j.issn.1000-1441.2015.03.016
史海霞、孟令媛、张雪梅等, 2018, 汶川地震前的b值变化, 地球物理学报, 61(5): 1874-1882.
宋维琪、吕世超, 2011, 基于小波分解与Akaike信息准则的微地震初至拾取方法, 石油物探, 50(1): 14-21. DOI:10.3969/j.issn.1000-1441.2011.01.002
谭玉阳、于静、冯刚等, 2016, 微地震事件初至拾取SLPEA算法, 地球物理学报, 59(1): 185-196.
唐礼忠、杨承祥、潘长良, 2006, 大规模深井开采微震监测系统站网布置优化, 岩石力学与工程学报, 25(10): 2036-2042. DOI:10.3321/j.issn:1000-6915.2006.10.014
田玥, 陈晓非, 2001. 关于一种新的地震定位方法的探讨. 见: 2001年中国地球物理学会年刊——中国地球物理学会第十七届年会论文集. 昆明: 云南科技出版社.
田玥、陈晓非, 2002, 地震定位研究综述, 地球物理学进展, 17(1): 147-155. DOI:10.3969/j.issn.1004-2903.2002.01.022
王璐琛、常旭、王一博, 2016, 干涉走时微地震震源定位方法, 地球物理学报, 59(8): 3037-3045.
王泽伟、李夕兵、尚雪义等, 2017, 基于VFOM的矿山微震震源定位及近震震级标定, 岩土工程学报, 39(8): 1408-1415.
余娜、张晓清、袁伏全, 2018, 青海地区地震序列参数的早期特征研究, 中国地震, 34(4): 695-703. DOI:10.3969/j.issn.1001-4683.2018.04.009
曾志毅、张建中, 2020, 利用微地震记录互相关成像的震源定位方法, 石油地球物理勘探, 55(2): 360-372, 388.
张楚旋、李夕兵、董陇军等, 2016, 顶板冒落前后微震活动性参数分析及预警, 岩石力学与工程学报, 35(增刊Ⅰ): 3214-3221.
张唤兰、朱光明、王云宏, 2013, 基于时窗能量比和AIC的两步法微震初至自动拾取, 物探与化探, 37(2): 269-273.
朱鹤文、韩立国、陈瑞鼎, 2020, 主成分分析和字典学习联合地震数据去噪方法, 世界地质, 39(3): 656-663. DOI:10.3969/j.issn.1004-5589.2020.03.015
朱权洁、姜福兴、于正兴等, 2012, 爆破震动与岩石破裂微震信号能量分布特征研究, 岩石力学与工程学报, 31(4): 723-730. DOI:10.3969/j.issn.1000-6915.2012.04.011
Allen R V, 1978, Automatic earthquake recognition and timing from single traces, Bull Seism Soc Am, 68(5): 1521-1532.
Bekara M, van der Baan M, 2009, Random and coherent noise attenuation by empirical mode decomposition, Geophysics, 74(5): V89-V98. DOI:10.1190/1.3157244
Chang S G, Yu B, Vetterli M, 2000, Adaptive wavelet thresholding for image denoising and compression, IEEE Trans Image Process, 9(9): 1532-1546. DOI:10.1109/83.862633
Daubechies I, Lu J F, Wu H T, 2011, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl Comput Harmon Anal, 30(2): 243-261. DOI:10.1016/j.acha.2010.08.002
Daubechies I, Maes S, 1996, A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models. In: Aldroubi A, Unser M. Wavelets in Medicine and Biology. Boca Raton: CRC Press, 527~546.
Donoho D L, Johnstone I M, 1994, Ideal spatial adaptation by wavelet shrinkage, Biometrika, 81(3): 425-455. DOI:10.1093/biomet/81.3.425
Donoho D L, Johnstone I M, 1995, Adapting to unknown smoothness via wavelet shrinkage, J Am Stat Assoc, 90(432): 1200-1224. DOI:10.1080/01621459.1995.10476626
Geiger L, 1912, Probability method for the determination of earthquake epicenters from the arrival time only, Bull St Louis Univ, 8: 60-71.
Guo H, Zhang H J, 2017, Development of double-pair double difference earthquake location algorithm for improving earthquake locations, Geophys J Int, 208(1): 333-348. DOI:10.1093/gji/ggw397
Gutenberg B, Richter C F, 1994, Frequency of earthquakes in California, Bull Seism Soc Am, 34(4): 185-188.
Hagen D C, 1982, The application of principal components analysis to seismic data sets, Geoexploration, 20(1-2): 93-111. DOI:10.1016/0016-7142(82)90009-6
Kao H, Shan S J, 2004, The source-scanning algorithm: mapping the distribution of seismic sources in time and space, Geophys J Int, 157(2): 589-594. DOI:10.1111/j.1365-246X.2004.02276.x
Li L, Tan J Q, Schwarz B, et al, 2020, Recent advances and challenges of waveform-ased seismic location methods at multiple scales, Rev Geophy, 58(1): e2019RG000667.
Liu X Z, Tang C A, Li L C, et al, 2018, Microseismic monitoring and stability analysis of the right bank slope at Dagangshan hydropower station after the initial impoundment, Int J Rock Mech Min Sci, 108: 128-141. DOI:10.1016/j.ijrmms.2018.06.012
Mousavi S M, Langston C A, 2017, Automatic noise-removal/signal-removal based on general cross-validation thresholding in synchrosqueezed domain and its application on earthquake data, Geophysics, 82(4): V211-V227. DOI:10.1190/geo2016-0433.1
Mousavi S M, Langston C A, Horton S P, 2016, Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform, Geophysics, 81(4): V341-V355. DOI:10.1190/geo2015-0598.1
Ogata Y, 1988, Statistical models for earthquake occurrences and residual analysis for point processes, J Am Stat Assoc, 83(401): 9-27. DOI:10.1080/01621459.1988.10478560
Ogata Y, 1989, Statistical model for standard seismicity and detection of anomalies by residual analysis, Tectonophysics, 169(1-3): 159-174. DOI:10.1016/0040-1951(89)90191-1
Ogata Y, 2001, Increased probability of large earthquakes near aftershock regions with relative quiescence, J Geophys Res: Solid Earth, 106(B5): 8729-8744. DOI:10.1029/2000JB900400
Palgunadi K H, Poiata N, Kinscher J, et al, 2020, Methodology for full waveform near real-time automatic detection and localization of microseismic events using high(8 kHz)sampling rate records in mines: application to the garpenberg mine(Sweden), Seism Res Lett, 91(1): 399-414. DOI:10.1785/0220190074
Qian J W, Zhang H J, Westman E, 2018, New time-lapse seismic tomographic scheme based on double-difference tomography and its application in monitoring temporal velocity variations caused by underground coal mining, Geophys J Int, 215(3): 2093-2104. DOI:10.1093/gji/ggy404
Salvoni M, Dight P M, 2016, Rock damage assessment in a large unstable slope from microseismic monitoring-MMG Century mine(Queensland, Australia)case study, Eng Geol, 210: 45-56. DOI:10.1016/j.enggeo.2016.06.002
Saragiotis C D, Hadjileontiadis L J, Panas S M, 2002, PAI-S/K: a robust automatic seismic P phase arrival identification scheme, IEEE Trans Geoscience Remote Sens, 40(6): 1395-1404. DOI:10.1109/TGRS.2002.800438
Senkaya M, Karsli H, 2014, A semi-automatic approach to identify first arrival time: the cross-correlation technique(CCT), Earth Sci Re J, 18(2): 107-113.
Ursin B, Zheng Y, 1985, Identification of seismic reflections using singular value decomposition, Geophys Prospect, 33(6): 773-799. DOI:10.1111/j.1365-2478.1985.tb00778.x
Waldhauser F, Ellsworth W L, 2000, A double-difference earthquake location algorithm: method and application to the northern Hayward Fault, California, Bull Seism Soc Am, 90(6): 1353-1368. DOI:10.1785/0120000006
Wang Z W, Li X B, Shang X Y, 2019, Distribution characteristics of mining-induced seismicity revealed by 3-D ray-tracing relocation and the FCM clustering method, Rock Mech Rock Eng, 52(1): 183-197. DOI:10.1007/s00603-018-1585-z
Wang Z W, Li X B, Zhao D P, et al, 2018, Time-lapse seismic tomography of an underground mining zone, Int J Rock Mech Min Sci, 107: 136-149. DOI:10.1016/j.ijrmms.2018.04.038
Xie W Y, Hattori K, Han P, 2019, Temporal variation and statistical assessment of the b value off the pacific coast of Tokachi, Hokkaido, Japan, Entropy, 21(3): 249. DOI:10.3390/e21030249
Zhang H, Thurber C, 2003, Double-difference tomography: the method and its application to the Hayward fault, California, Bull Seism Soc Am, 93(5): 1875-1889. DOI:10.1785/0120020190
Zhang X, Zhang H J, 2015, Wavelet-based time-dependent travel time tomography method and its application in imaging the Etna volcano in Italy, J Geophys Res: Solid Earth, 120(10): 7068-7084. DOI:10.1002/2015JB012182