BSN矿山微震监测技术及其应用

2021-07-21 09:50张达戴锐曾志毅冀虎石雅倩常莹韩鹏

中国地震 2021年2期

张达戴锐曾志毅冀虎石雅倩常莹韩鹏

1)矿冶科技集团有限公司矿山工程研究设计所，北京 102628 2)中国-南非矿产资源开发利用联合研究中心，北京 102628 3)中国-南非矿产资源可持续开发利用“一带一路”联合实验室，北京 102628 4)南方科技大学，深圳市深远海油气勘探技术重点实验室，广东深圳 518055 5)南方科技大学，地球与空间科学系，广东深圳 518055

0 引言

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

随着高精度、高采样率传感器的研发和计算机技术的飞速发展，矿山微震高精度定位、实时监测已成为可能。随着微震定位技术的不断完善，微震监测已能够实现矿山危险源的时空定位及稳定性评估，近年来广泛应用于深井矿山安全监测，已成为当前行业高度认可的地压灾害监测预警手段(Gibowicz et al，1998；田玥等，2001；李庶林等，2005；姜福兴等，2006；唐礼忠等，2006；冯夏庭等，2013；王璐琛等2016；Palgunadi 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)

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

图 2 BSN微震数据采集站

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

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

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

1.2 微震数据处理

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

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

1.2.1 微震信号降噪

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

(1)

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

(2)

(3)

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

本微震监测系统中实现了带通滤波方法和基于同步挤压小波变换(Daubechies et al，1996、2011)的小波系数硬、软阈值滤波技术，滤波效果如图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算法公式为

(4)

CF(j)=Y2(j)

(5)

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

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

AIC(i)=ilg(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)特征作为微震信号初至波拾取的特征函数，利用偏度、峰度极值点位置信息拾取初至波。其计算公式如下

(7)

(8)

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

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

图 7 两种自动拾取方法的初至波拾取结果与人工拾取结果的误差统计分布

1.2.3 微震事件定位

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

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

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

(9)

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

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

(10)

AΔX=r

(11)

其中，系数矩阵A为N×4矩阵，r为N×1矩阵，震源校正量ΔX为4×1矩阵，则

(12)

(13)

(14)

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

1.2.4 地震活动性分析

微震活动性分析是微震监测研究的重要环节，其可以更快速地核定矿区地震危险性，方便矿山管理者选择相应的对策措施，在保证矿山正常生产的情况下，确保采矿工人的安全。地震活动性分析方法包括地震活动的时间和空间分析，较为简单的分析方法是利用震源位置的空间分布了解矿区内岩石破裂的区域，以及用微地震事件位置表示裂缝的空间形态。为了进一步利用微震事件位置信息，许多学者通过聚类方法实现了对微震事件的划分(Wang et al，2019；刘德彪等，2019)，减少了人工主观划分微震事件的不确定性，发现潜在的微震集群，从而有效地分析微震事件分布特征和活动规律。震级-频度关系(lgN=a-bM)是统计地震学研究中最基础的规律之一(Gutenberg et al，1994)，其中a值表征地震活动水平，b值表征地震频次在不同震级上分配的比例。b值是判定地震危险性区域的重要指标，其可以为地震危险性方面的研究提供依据(史海霞等，2018；Xie et al，2019)，目前被广泛应用于地震预测和地震危险性研究中。高b值对应低应力，低b值则对应高应力，因此b值可以间接判断地下应力的状态和变化，进而分析矿区地震活动性(张楚旋等，2016)。传染性余震序列(ETAS)模型(Ogata et al，1988、1989、2001；余娜等，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为垂向。从XZ和YZ面可以看出微震事件主要集中在-5m中段，与该时间段内的矿山开采深度一致。本系统采用盖格定位方法，对定位结果的误差估计由地震波到时的观测值与计算值之间的差，即时间残差来表示；而定位结果在空间坐标XYZ上的误差由时间残差和走时方程系数矩阵的协方差矩阵来估计。微震事件定位误差统计分布结果如图12 所示，可以看到X、Y和Z方向的定位误差均值大致为2.5m，均方根误差在2.5～3.8m之间，总体误差分布较小，表明定位结果较为精确。通过XY面的微震事件空间分布特征可以看到微震活动大致分布在3个区域(图11 中红色虚线椭圆所示)，微震事件较为集中，由此可初步判断这3个区域的地震活动性较高。

图 11 微震事件时空分布

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

图 12 定位结果在X、Y、Z方向上的误差统计分布

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

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

3 结论和讨论

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

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