一种微震震源定位方法及系统

本发明涉及微震震源定位方法及系统领域，特别是涉及一种微震震源定位方法及系统。

背景技术：

微震是指外界条件作用下，固体物质内部局部应力集中区，在由初始高能状态向低能状态过渡时，应变能以弹性波的方法快速释放，频率在10-300hz的声发射现象。

微震监测技术源于对地壳自发地震的研究，后应用于矿山地震监测领域，以人工地震和天然地震发生的强度较弱的应力波为监测对象，突破了传统监测设备监测范围在时间与空间上的限制，实现了实时全范围全天候的立体监测。国际上南非、加拿大等国的金属矿山在地压监测应用中取得良好效果。

1、矿山开采过程中，采空区分布未知与矿权越界采取等问题普遍存在，需要对矿山开采过程进行监测。需要借助现代化地震监测技术及地震信号数据处理手段，对矿山井下活动进行准确追踪测定，矿山开采过程中采用爆破方式进行，而爆破会在地壳产生一定的震动。利用微震追踪定位，得到地下开采面的空间位置形态，是一种动态、实时、不间断的监测技术，可用于监测和追踪各种采矿行为。

2、在采矿的过程中，随着采空区的增加和施工影响区域的增多，岩层地质结构会相应的改变，某些采空区域的岩层就可能会用应力失衡的条件下就会发生微小破裂，产生声发射现象，即以声波的形式向周围释放能量，频率在300hz以下的称为微震。这些微小破裂的增多，最终导致大的矿震，危害人们的生命安全。如果能够对地壳微震情况进行监测，及时检测到这些信息，通过计算分析监测的微震信息，确定微震发生的大小及位置，以及微震事件累积数量等，即可以及早发现，提前判断矿山岩体稳定性并提出合理预测，实现高精度、远距离、动态、实时监测，从而起至灾害报警伤脑筋，提前预防灾害发生，降低对人们的危害，遏制重大事故。

无论上述哪种应用，首先要采集到微震事件，通过建立软件，基于各种算法实现所记录信号的分析识别，利用地震定位方法，实现震源定位和追踪。在此过程中，定位算法是定位系统的核心

在微震源定位算法中，全局优化算法对于传统算法有着一定的优势，其对复杂地质的适应性更好，能够建立其更为真实的速度模型，但运算量过于巨大。因此在提高优化运算速度的前提下，如果快速且准确地实现定位是定位算法的最重要目标。由于单一算法的局域性，单一算法往往实用性不高。

技术实现要素：

本发明的目的是提供一种微震震源定位方法及系统，采用以网格搜索法为辅，粒子群算法为主体的算法进行微震震源定位，能够提高定位精度和速度。

为实现上述目的，本发明提供了如下方案：

一种微震震源定位方法，包括：

采用传感器采集微震震源数据；

根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻；

根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置；

根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置。

进一步地，所述根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻，具体包括：

根据所述震源数据，构建特征函数；

根据所述特征函数，采用能量比法构建特征函数比值；

根据所述特征函数比值，基于中值滤波法确定动态阈值；

根据所述动态阈值以及所述特征函数比值确定微震波形的初至时刻。

进一步地，所述特征函数的表达式如下：

fc(i)＝a(i)²-a(i-1)a(i+1)

其中，fc(i)表示特征函数，c表示传感器，a(i)表示微震波形的离散的振幅信号，i＝1,2,...,n，n为记录的离散点数。

进一步地，所述根据所述动态阈值以及所述特征函数比值确定微震波形的初至时刻，具体包括：

当所述动态阈值大于所述特征函数比值时，确定微震波形未至；

当所述动态阈值小于所述特征函数比值时，确定微震波形初至。

进一步地，所述根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置，具体包括：

将震源所在区域作网格划分；

基于所述传感器的坐标位置、所述初至时刻和每个网格点坐标，计算震源的发震时间；

将所述初至时刻与所述发震时间进行对比，基于残差特性分析，取残差最小值对应的网格点空间作为初始位置。

进一步地，所述根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置，具体包括：

基于所述得到的初始位置作为边界值；

根据所述传感器的坐标位置、所述初至时刻以及所述初始位置，确定粒子群内的粒子数量，粒子的位置信息及粒子的速度范围；

根据设定的波速，以传感器的初至时刻为时间原点，代入粒子群算法进行计算，得到每个粒子对应的残差值；

取残差值最小的粒子的坐标和时间为震源的准确位置和发震时间

本发明还提供了一种微震震源定位系统，包括：

数据采集模块，用于采用传感器采集微震震源数据；

初至时刻确定模块，用于根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻；

初始位置确定模块，用于根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置；

震源位置确定模块，用于根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置。

进一步地，所述初始位置确定模块，具体包括：

特征函数构建单元，用于根据所述震源数据，构建特征函数；

比值构建单元，用于根据所述特征函数，采用能量比法构建特征函数比值；

动态阈值构建模块，用于根据所述特征函数比值，基于中值滤波法确定动态阈值；

初至时刻确定单元，用于根据所述动态阈值以及所述特征函数比值确定微震波形的初至时刻。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种微震震源定位方法及系统。方法包括：采用传感器采集微震震源数据；根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻；根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置；根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置。本发明首先利用网格搜索法缩小了粒子群随机搜索的范围，避免了粒子限入局部最小值，引起大的定位误差。因此，粒子群算法只需要在小的搜索范围完成任务，意味着只需要较小的粒子数量，即可达到足够的定位精度和可靠性。从而实现了在较小的运算量的前提下，提高了定位精度和可靠性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例微震震源定位方法的流程图；

图2为传感器信号的识别结果。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种微震震源定位方法及系统，采用以网格搜索法为辅，粒子群算法为主体的算法进行微震震源定位，能够提高定位精度和速度。

在微震震源定位算法中，全局优化算法对于传统算法有着一定的优势，其对复杂地质的适应性更好，能够建立其更为真实的速度模型，但运算量过于巨大。因此在提高优化运算速度的前提下，如果快速且准确地实现定位是定位算法的最重要目标。由于单一算法的局域性，单一算法往往实用性不高。为了实现高精度、高运算速度的定位算法，提出先计算各点的前后窗口能量比的主算法，再结合改进后的中值滤波法对能量比处理，筛选出符合阈值判断的点，实现初至识别。然后根据传感器坐标和初至数据，采用网格搜索法及粒子群法相结合的综合算法，以网格搜索法为辅，粒子群算法为主体的算法，实现对震源的高精度和高速度定位。

网格搜索法的重点在于其对速度模型的建立，例如在均匀介质中速度为常量，震源至监测点的位置为一条直线，而在实际中由于介质的不同，地震波在传播过程中发生折射，因此需利用snell定律对传播路径进行修正。

单一采用网格搜索法时，为提高反演定位精度，划分网格时所用尺度要小，可靠性才高，但运算量十分巨大，每个网格点都需对其路径进行修正计算。特别是对多参数优化问题，计算量会随着网格尺度的减小而呈几何级数增加。

但网格搜索法也是全局优化算法中最基本、最简单、最直观且最有效的方法。虽然对于单次微震定位来说数据量过大，但如果与其它算法结合，将网格搜索法作为基础，实现大尺度的定位范围的确定，可大大缩短定位时间。

此时，本发明引入粒子群优化算法，利用粒子群算法实现在上述定位范围内的高精度、高可靠性的震源定位。粒子群优化算法是通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法。

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。将粒子群算法引入微震定位领域，则转变为在三维空间中寻找走时残差最小的粒子群问题。其位置为一个三维坐标，寻找最小走时残差。

粒子群算法同样属于全局优化算法，当进行运算时，当粒子数量不足够时，很容易使定位陷入局部最优解，导致定结果发散，使定位可靠性降低。此时就需要足够数量的粒子来避免粒子陷入局部最优导致结果发散的情况。而大量的粒子则会大量增加运算量。

本发明的方法首先利用网格搜索法缩小了粒子群随机搜索的范围，避免了粒子限入局部最小值，引起大的定位误差。这样的话，粒子群算法只需要在小的搜索范围完成任务，意味着本发明只需要较小的粒子数量，即可达到足够的定位精度和可靠性。也就是实现了在较小的运算量的前提下，提高了定位精度和可靠性。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1所示，一种微震震源定位方法包括以下步骤：

步骤101：采用传感器采集微震震源数据。

步骤102：根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻。具体包括：根据所述震源数据，构建特征函数；根据所述特征函数，采用能量比法构建特征函数比值；根据所述特征函数比值，基于中值滤波法确定动态阈值；根据所述动态阈值以及所述特征函数比值确定微震波形的初至时刻。

通过传感器收集到的数据，本发明可以从中提取到地震波的初至时刻，具体思路是：利用传感器收集到的源数据构建特征函数，特征函数可用如下表达式确定：fc(i)＝a(i)²-a(i-1)a(i+1)，式中a(i),i＝1,2,...,n，a(i)是微震波形的离散的振幅信号，c代表某一传感器，n为记录的离散点数。该特征函数能同时增强振幅和频率的变化。由于存在噪声，噪声具有一定范围内的随机性，所以本发明将特征函数构建为源数据相邻两项的平方差，对于有效信号，特征函数的数值会随着振幅的增大而增大。

由于传感器存在噪声，所以不能直接通过对数值的测量来判断地震波的初至，但是由于地震波携带着能量，在地震波到来后传感器的波形能量发生明显改变。本发明的方法是采用多算法结合的方式对微震波形的初至时刻进行提取。主算法首先是基于前后窗能量比法，构建特征函数的比值，其表达式为

其中：xi代表特征函数比值，fc(l)为先前所构建的特征函数，c表示某一传感器。i代表离散数据的序数，nlta和nsta分别代表长时窗(lta)和短时窗(sta)所选取的数据点窗口。

而对于噪声，特征函数的数值在一定小范围内进行浮动。构建好特征函数后，再利用这个特征函数构建判断地震波初至的阈值，其目的在于对噪声部分通过阈值来将其覆盖，有效信号部分获得输出，从而达到提取初至的目的。此时用到的算法即为中值滤波法，基于特征函数比值具有的周期特征，构建如下动态阈值表达式：

其中med代表微震数据特征函数比值的所有数据点的中值，std{x1,x2,x3,…xn}表示所有数据点的标准差。stdmax表示最大值，n为构建的窗口长度。利用此方法可以增加噪声状态的阈值，减小有效信号的阈值，从而可以提高提取的初至时刻的精确度。

对于动态阈值，它有以下的特点：在地震波初至未至时，阈值永远大于特征函数比值，在地震波初至后，特征函数比值会迅速增大，此时的动态阈值会小于特征函数比值，刚则大于给定动态阈值的时刻即为微震波的初至时刻，通过步骤102得到了微震波形的初至时刻。

本步骤可以概括为：利用能量比法来构建一个特征函数比值，以便反映出局部数据相比于整体数据的特殊性，用来与阈值比较判断初至时刻。根据特征函数比值，采用改进中值滤波法构建动态阈值，筛选出符合阈值判断的点，实现微震波形初至提取。该初至信息将保存在数据库中，保存提供给下一步使用。

步骤103：根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置。

在步骤102中得到了微震的初至时刻信息，该信息与微震波形数据传感器的坐标位置信息用来实现微震源的准确初值点(即初始位置)的确定。

从步骤102得到微震的初至时刻信息，该信息被保存在后台数据库中，在步骤103中，将步骤102保存的初至信息导入到步骤103中，步骤103还利用了事先输入的传感器坐标位置信息。利用网格搜索算法进行微震源的准确初值点的确定。设在一个均匀单一介质模型中进行定位，设定相应范围，如取介质模型尺寸为30*30*40km，设定微震波的传播波速，如设定波速为5700m/s。将震源所在区域作网格体划分，网格以某一固定距离为间隔，确定每一维度上划分的网格数，在三维空间建立一个一定大小的三维矩阵，每个网格以网格点表示，建立以网格点坐标表示的网格空间；根据介质特性设定微震波的传播速度；基于所述传感器的坐标位置、所述初到时刻和每个网格点坐标，基于走时方程计算震源的发震时间；将所述初至时刻与发震时间进行对比，基于残差特性分析，取残差最小值对应的网格点空间作为初始位置。

在微震事件未发生时，还可利用网格搜索法根据提前建立的速度模型，对各网格点发生微震事件时的初至时刻进行计算，并将各网格点模拟的数据载入数据库，在事件发生时，将与事件信息最匹配的网格点作为初始位置，载入粒子群定位算法即可得到准确位置信息。虽然网格搜索法运算量大，但其大部分运算均可在事件发生前完成，因此仍是一种快速的定位方法，在多事件定位中有着不错的优势。

对于在网格足够细密，速度模型足够准确时，仅使用网格搜索法即可达到准确的定位。而在无法对地质情况有准确的了解时，则需采用较大的网格来抵消各小网格之间的误差。本步骤中采用的为基于走时的网格搜索法，对于不同的粒子群算法均可有较好的兼容。

步骤104：根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置。

将步骤103中得到的初始位置导入到步骤104中，作为本步骤计算的边界值。步骤104将该信息与微震波形数据传感器的坐标信息和步骤102得到的初至时刻实现微震源的准确定位。步骤104中主要采用粒子群算法实现上述过程。粒子群算法依据这些参量，确定需要使用的群内粒子数量，而粒子应在步骤103给定的一个接近真实位置的初始位置内，由于采用网格搜索算法进行了初始位置确定，粒子群距离真实震源位置较近，一方面粒子群内可以设定较少的粒子数量，有效减小计算量；另一方面可以对粒子速度进行适当修正，以保证震源的定位精确度。根据设定的波速，取某一传感器的初至时刻为时间原点，代入粒子群算法的走时公式，并计算残差值，基于粒子群算法规则设定的自适应值，通常设定为时间残差的平方和，将具有不同时间和位置属性的粒子不断的更新迭代，在迭代过程中不断地接近震源的准确位置。在满足残差最小值条件后，迭代自动退出，此时大部分粒子集中于震源真实位置附近，可以取残差值最小的粒子的坐标和时间为震源的准确位置和发震时间，或取其中心位置也可。基于先前的网格搜索算法优化了初始位置，只需要少量的迭代次数和少量的粒子群数量即可以确定震源的准确位置和发震时间。

此时粒子群算法的计算结果，可以有效地减少由于偶然因素所导致的误差，可避免在不断的迭代过程中会有一部分粒子陷入局部最小值，影响最终结果。

网格搜索+粒子群算法通过网格搜索确定了一个较为准确的初始位置值，使粒子群在开始运算时就集中分布在空间距离真实震源较近的某一片区域，经过粒子群算法得到的就是准确度和精确度极高的震源位置。在运算结束后，粒子群集中在真实值附近，从宏观上来看像是粒子群的中心位置移动到了真实值位置。

举例说明：

对某一单通道传感器收集到的一组原始地震波波形信号(图2(c)所示)进行处理：

对图2(c)中地震波信号波形根据定义，计算得到特征函数，如图2(b)所示。利用能量比法计算得到特征函数比值，再采用中值滤波方法得到动态阈值，特征函数比值和动态阈值同步显示在图2(a)中。在图2(a)中，幅值较高的曲线代表中值滤波法得到的各点动态阈值大小，幅值较低的曲线是利用能量比法得到的各处特征函数比值，通过两者比较即可判断微震到达该传感器的初至时刻。

从图2中可以直观地看出，本发明提供的方法处理的阈值一直远远大于能量比，将通过此算法识别出的地震波初至时刻与地震波振幅随时间变化以及地震波波形的变化比较发现，该算法给出的第一个初至时刻点前几乎没有地震波能量，初至点后地震波能量依次增大。可以看出，本发明采用的方法对微震初至识别准确，排除了噪声的干扰。识别速度快。算法依赖能量比法得到特征函数比值和中值滤波法得到的动态阈值作比较，因此在使用前需要选择合适的窗口。该算法可适用于利用某一单通道传感器来检测定位地震波震源的系统之中。

本发明还提供了一种微震震源定位系统，包括：

数据采集模块，用于采用传感器采集微震震源数据；

初至时刻确定模块，用于根据所述震源数据，通过能量比法和中值滤波法确定微震波形的初至时刻；

初始位置确定模块，用于根据所述初至时刻以及所述传感器的坐标位置，采用网格搜索法确定微震震源的初始位置；

震源位置确定模块，用于根据所述初始位置、所述传感器的坐标位置以及所述初至时刻，采用粒子群算法确定微震震源的准确位置。

其中，所述初始位置确定模块，具体包括：

特征函数构建单元，用于根据所述震源数据，构建特征函数；

比值构建单元，用于根据所述特征函数，采用能量比法构建特征函数比值；

动态阈值构建模块，用于根据所述特征函数比值，基于中值滤波法确定动态阈值；

初至时刻确定单元，用于根据所述动态阈值以及所述特征函数比值确定微震波形的初至时刻。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。