一种微震震源自动定位及可靠性综合评价方法与流程

本发明涉及一种微震震源自动定位及定位结果评价方法，尤其是一种适用于矿山、隧道、边坡和油气开发等各种岩土工程的微震震源自动定位及定位可靠性评价的方法，属于微震监测领域。

背景技术：

微震监测技术是利用煤岩材料受载破坏产生的微震波形信号来研究和评价煤岩稳定性的一种地球物理实时监测技术。近年来，微震监测技术在矿山煤岩动力灾害监测预警、油气开发、边坡和堤坝稳定性监测等各类岩土工程中得到了广泛应用。微震震源定位是微震监测技术的重要内容和研究热点，随着科学技术的发展和现场要求的提高，高精度的微震震源自动定位及定位结果可靠性(定位精度和震源求解系统稳定性)评价是微震监测的发展趋势。微震监测台网布设、微震波形识别和到时拾取、波速模型、定位算法是影响微震震源定位的主要因素。其中，P、S波及异常波形识别、高精度到时自动拾取和稳健的定位算法是高精度自动定位关键，同时也需要对自动定位结果可靠性给出正确的评价。

目前微震波形主要通过时频、幅频、持续时间等特征进行分析，对背景噪音的识别和滤除效果较好，然而对于P、S波和异常波形识别难度较大。文献《岩石力学与工程学报》，2014年第8期，李楠，“基于到时差值的微震波异常信号识别方法”，该文献根据到时差值分析对延迟波和外部异常波两种异常波形信号进行识别，没有涉及P、S波和其它异常波形识别。微震波形到时拾取主要采用长短时窗口(STA/LTA)比值法，但当波形信号信噪比较低或初动不明显时到时拾取误差较大，而且该类方法没有利用波形信号的频率特征来拾取到时。在微震定位方面，应用最多的就是基于到时不同原理发展而来的各种震源定位方法，例如Geiger法以及在此基础上改进的各种方法，然而这类算法对初始值依赖较高，当初始值选取不当，定位精度将大大降低，甚至不能收敛。近年来提出的Powell法、粒子群算法等一定程度上降低了对初始值的依赖程度，但是并没有从根本上解决求解系统发散问题。文献《煤炭学报》，2014年第12期，李楠，“基于L1范数统计的单纯形微震震源定位方法”，该文献采用基于L1范数统计准则的单纯性算法进行震源定位，未对到时拾取和异常波形识别进行论述，也没有同时采用P、S波进行定位。另外，目前微震波形到时自动拾取精度和微震震源自动定位精度也不能满足现场需求，很多情况下需要人工进行修正，这不仅大大增加了工作量，费时费力；而且人工定位受个体人主观判断影响较大。目前微震定位结果评价方面主要考虑对定位精度的评价，然而微震震源定位评价不仅应包括定位精度评价，还应包含定位稳定性评价，即定位可靠性综合评价。文献《煤炭学报》，2013年第11期，李楠，“微震震源定位可靠性综合评价模型”，该文献提出了一种能够对定位精度和稳定性进行全面评价的模型，然而该模型没有将涉及微震震源定位方法，而仅仅是对已有的定位结果进行评价。因此，缺少一种同时具有微震震源自动定位和定位可靠性评价功能于一体的综合方法。

因此，发明一种集微震波形到时自动拾取、微震波形自动识别、微震震源自动定位、定位可靠性综合评价于一体的微震震源定位及可靠性综合评价方法是非常有意义的。

技术实现要素：

技术问题：针对上述问题，本发明的目的在于提供一种微震震源自动定位及可靠性综合评价方法，通过建立矿山、隧道、边坡和油气开发等微震监测中微震波形到时自动拾取算法、P、S波及异常波形自动识别方法、稳健的自动定位算法、微震震源定位可靠性综合评价体系，最终形成一种高精度的微震震源自动定位及可靠性综合评价方法，该方法能够大大提高微震震源自动定位精度，并且给出正确的定位可靠性评价结果，从而为微震监测的成功应用奠定基础。

技术方案：一种微震震源自动定位及可靠性综合评价方法，包括基于瞬时频率和可变分辨率的包络函数的微震波形到时自动拾取算法IFEP、基于到时差值和残差分析微震波形自动识别方法APSI、单纯形微震震源自动定位算法Simplex、微震震源定位可靠性综合评价体系SLRE四大部分；

所述IFEP能够自动拾取高精度的微震波形到时数据，为APSI和Simplex奠定基础；

所述APSI在IFEP自动拾取到时的基础上，基于到时差值分析和残差分析，结合微震传感器空间布设形态，对微震波形中的延迟波和外部异常波两类异常波形进行自动识别剔除，同时还能够确定各个传感器拾取微震波的P、S波类型，为Simplex奠定基础；

所述Simplex能够正确采用对应的P、S波到时和波速进行震源自动定位，并且能够任意选择L1范数统计准则和L2范数统计准则计算事件残差；

所述SLRE能够对微震震源定位可靠性(定位精度和稳定性)进行有效评价；

IFEP、APSI、Simplex、SLRE四部分有序结合，互为前提和基础，环环相扣，形成一套完整的微震震源自动定位和可靠性评价方法。

所述基于瞬时频率和可变分辨率的包络函数的微震波形到时自动拾取算法IFEP，具体步骤如下：

(1)将微震波形时序信号进行小波变换，得到波形信号的时频二维联合分布，通过求解波形信号时频二维分布的边缘分布，得到波形信号频率随时间的变化关系；公式如下：

$f\left(t\right)=E(f,t)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{\mathrm{fg}}_{f,t}{d}_f=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}fG(f,t)d_f$

从而计算得到波形信号的任意时刻的瞬时频率；式中：g_f,t代表某一频率成分f在某个时刻t出现的概率；

(2)设定包络分析分辨率d；求取波形信号在分辨率d时的所有极值点，获得该分辨率下的极值点时序数据；采用三次样条插值，获得完整的包络函数序列；改变分辨率d，从而得到不同分辨率的包络函数序列，最终获得不同分辨率信号包络函数序列；

(3)将步骤(1)中得到的微震波形信号的瞬时频率和步骤(2)中得到的可变分辨率的包络函数作为时间序列输入，分别计算包络信号和瞬时频率的长短时窗比值序列R1和R2，令R＝R1×R1，将R值与设定的阈值进行比较，从而自动拾取微震波形到时，将自动拾取的微震波形到时称为观测到时，第i个传感器的观测到时记为t_i。

所述基于到时差值和残差分析微震波形自动识别方法APSI，具体步骤如下：

设任意两传感器T_i和T_j之间的欧式距离为2c_ij，根据权利要求2中所述的IFEP观测到时分别为t_i和t_j，并令t_i＜t_j；现场实测微震监测区域中P波波速为v_P，S波波速为v_S；分别计算传感器之间的观测到时差值OAD_ij和P波到时差值理论极限TLP_ij：

OAD_ij＝t_j-t_i

${\mathrm{TLP}}_{ij}=\frac{2c_{ij}}{v_p}$

设传感器T_i是P波到时，T_j是S波到时，当两传感器和震源位于同一直线上时，它们之间的观测到时差值取得极大值，且该极大值计算公式如下：

${\mathrm{DAD}}_{ij}=\frac{d_i+2c_{ij}}{v_s}-\frac{d_i}{v_p}=d_i\left(\frac{v_p-v_s}{v_p{v}_s}\right)+\frac{2c_{ij}}{v_s}$

式中：DAD_ij记为传感器T_j的延迟波到时差值理论极限，d_i为传感器T_i到测区域边缘的最大距离；

根据观测到时差值OAD_ij、P波到时差值理论极限TLP_ij、延迟波到时差值理论极限DAD_ij，建立到时差值分析表1：

表1到时差值分析表

表中：

微震传感器台站残差采用下式进行表示：

γ_i＝t_i-tt_i

式中：t_i是观测到时，tt_i是计算走时；其中t_i不受计算的影响，因此台站残差γ_i主要受计算走时tt_i的影响；建立残差分析表2：

表2残差分析表

表中：k为残差系数，且k＞0。

根据到时差值分析、残差分析，同时结合微震传感器台网空间布设形态，对微震波形中的延迟波和外部异常波两类异常波形进行识别剔除，同时还能够识别确定各个传感器拾取微震波的P、S波类型，为震源定位奠定基础。

所述单纯形微震震源自动定位算法Simplex，具体步骤如下：

(1)首先将微震发震时间t₀分离出去，采用单纯形微震震源自动定位算法Simplex对震源空间坐标(x₀,y₀,z₀)进行求解，然后再计算发震时间t0；

(2)根据权利要求2所述的IFEP微震波形到时自动拾取结果和权利要求3所述的APSI微震波形自动识别结果，剔除异常波形信号，同时采用正确的、对应的P波和S波到时和波速进行震源定位；

(3)选择采用L1范数统计准则或L2范数统计准则计算事件残差；

基于L1范数统计准则的事件残差计算公式为：

$E_{L1}=\frac{\underset{i=1}{\overset{n}{\Σ}}|t_i-t_{ti}-t_{median}|}{n-m}$

基于L2范数统计准则的事件残差计算公式为：

$E_{L2}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\γ}}_i^2}{n-m}}$

式中：t_median是所有(t_i-tt_i)的中位数；n是有效微震传感器台站的个数；m表示自由度，即未知数的个数。

所述微震震源定位可靠性综合评价体系SLRE，包含事件残差指标R、敏感度指标M和触发序列指标T，并且综合上述三种指标建立了微震震源定位可靠性综合评价指标P，能够对微震震源定位可靠性进行综合评价；具体步骤如下：

(1)利用权利要求4所述的单纯形微震震源自动定位算法Simplex得到微震震源自动定位结果后，SLRE能够自动对定位精度和稳定性进行综合评价；

(2)根据L1范数统计准则计算得到的事件残差大小，确定事件残差指标R的评价标准；微震波波速增加或降低10％后，采用Simplex进行重新定位，根据重定位结果与原来的定位结果之间的误差L的大小，确定敏感度指标M的评价标准；根据观测触发序列和计算触发序列之间的错配程度，同时结合微震有效传感器个数，确定触发序列指标T的评价标准；

(3)引入了微震震源定位综合评价指标P，并有P＝R+M+T，根据P的评分情况，SLRE将震源定位结果划分为5个等级，分别为A、B、C、D、E，实现了微震震源定位可靠性定量化评价。

有益效果：本发明是一种系统的微震震源自动定位以及定位可靠性综合评价方法(AMSLRE)，实现了微震波形到时自动拾取、微震波形自动识别、震源自动定位以及定位可靠性综合评价，其主要优点在于：

(1)IFEP基于瞬时频率和可变分辨率的包络函数进行微震波形到时自动拾取，兼顾了微震波形的幅值和频率信息，对于信噪比较低的微震波形信号，极大的提高了到时自动拾取精度，从而为微震波形自动识别和微震震源自动定位奠定基础。

(2)APSI不需要通过复杂的波形分析对微震波类型进行识别，而是基于到时差值分析、残差分析以及微震传感器台网布设形态对微震波形中的延迟波和外部异常波两种异常波形进行自动识别和剔除，同时还能够识别确定各个传感器拾取微震波的P、S波类型，为采用对应的P、S波到时和波速进行震源定位奠定基础。

(3)单纯形微震震源定位算法(Simplex)不需要求解偏导数和逆矩阵，大大降低了计算量，而且不会出现发散问题，从而使计算过程更加稳定和强健。另外Simplex能够同时采用对应的P、S波到时和波速进行震源定位，而且也能够灵活采用L1范数统计准则或L2范数统计准则，从而大大提高微震震源自动定位精度。

(4)微震震源定位可靠性综合评价体系(SLRE)综合了事件残差指标R、敏感度指标M和触发序列指标T，不仅能够对微震震源定位精度进行评价，还能够对定位求解系统和定位结果稳定性进行有效评价；SLRE给出了综合评价指标P的量化评价标准以及对应的震源定位划分等级，实现了微震震源定位可靠性定量化评价。

(5)AMSLRE是IFEP、APSI、Simplex和SLRE的有序结合，是一套系统完整的微震震源自动定位和可靠性评价方法，IFEP、APSI、Simplex和SLRE环环相扣、互为前提和基础，最终实现高精度微震震源自动定位和定位可靠性综合评价。

附图说明

图1是延迟波到时差值分析示意图。

图2是单纯形在误差空间中的移动实例。

图3是一种微震震源自动定位及可靠性综合评价方法(AMSLRE)主要内容和流程图。

图4是IFEP微震波形到时自动拾取实例。给出了其中两种不同信噪比的现场微震波形采用IFEP自动拾取到时的实例，图4(a)为高信噪比的现场微震波形采用IFEP自动拾取到时的实例图，图4(b)为低信噪比的现场微震波形采用IFEP自动拾取到时的实例图。

具体实施方式

下面结合实施例对本发明作进一步的描述：

本发明的微震震源自动定位及可靠性综合评价方法，具体实施方案如下：

(1)基于瞬时频率和可变分辨率的包络函数的微震波形到时自动拾取算法(IFEP)根据小波分析理论，将微震波形时序信号进行小波变换，得到时频二维联合分布：

g_f,t＝G(f,t) (1)时频图上任一点的值g_f,t代表某一频率成分f在某个时刻t出现的概率。

通过求解波形信号时频二维分布的边缘分布，得到波形信号频率随时间的变化关系。其公式如下：

$f\left(t\right)=E(f,t)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{\mathrm{fg}}_{f,t}{d}_f=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}fG(f,t)d_f---\left(2\right)$

公式(2)即为波形信号的瞬时频率计算公式，根据公式(2)能够准确得到波形信号的任意时刻的瞬时频率，从而为从波形信号的频率特征自动拾取到时奠定基础。

包络分析是工程中常用的一种信号分析方法，形成包络线的基本特征点是波形信号的极值点。可变分辨率的包络函数实现步骤具体如下：

a.求出波形信号的全部极值点；

b.设定包络分析分辨率d；

c.求取波形信号在分辨率d时的所有极值点，获得该分辨率下的极值点时序数据；

d.采用三次样条插值，获得完整的包络函数序列；

e.返回第一步，重复上述步骤，得到不同分辨率的包络函数序列。最终，获得不同分辨率信号包络函数序列。

将微震波形信号的瞬时频率和可变分辨率的包络函数作为时间序列输入，分别计算包络信号和瞬时频率的长短时窗比值序列R1和R2，令R＝R1×R1，将R值与设定的阈值进行比较，从而自动拾取微震波形到时，将自动拾取的微震波形到时称为观测到时，第i个传感器的观测到时记为t_i。从而为微震波形自动识别和微震震源自动定位奠定基础。

(2)基于到时差值和残差分析微震波形自动识别方法(APSI)

设任意两传感器T_i和T_j之间的欧式距离为2c_ij，它们观测到时分别为t_i和t_j，并令t_i＜t_j；现场实测微震监测区域中P波波速为v_P，S波波速为v_S；根据公式(3)和公式(4)分别计算传感器之间的观测到时差值OAD_ij和P波到时差值理论极限TLP_ij：

OAD_ij＝t_j-t_i (3)

${\mathrm{TLP}}_{ij}=\frac{2c_{ij}}{v_p}---\left(4\right)$

如图1所示，设传感器T_i是P波到时，T_j是S波到时，当两传感器和震源位于同一直线上时，它们之间的观测到时差值取得极大值，且该极大值计算公式如下：

${\mathrm{DAD}}_{ij}=\frac{d_i+2c_{ij}}{v_s}-\frac{d_i}{v_p}=d_i\left(\frac{v_p-v_s}{v_p{v}_s}\right)+\frac{2c_{ij}}{v_s}---\left(5\right)$

上式中DAD_ij记为传感器T_j的延迟波到时差值理论极限，d_i为传感器T_i到测区域边缘的最大距离，采用公式(5)计算延迟波到时差值理论极限。

根据观测到时差值OAD_ij、P波到时差值理论极限TLP_ij、延迟波到时差值理论极限DAD_ij，建立到时差值分析表(表1)。

表1到时差值分析表

注：其中

微震传感器台站残差可用下式进行表示：

γ_i＝t_i-tt_i (6)

上式中，t_i是观测到时，tt_i是计算走时。其中t_i不受计算的影响，因此台站残差γ_i主要受计算走时tt_i的影响。根据公式(6)，建立残差分析表(表2)。

表2残差分析表

注：表中k为残差系数，且k＞0，应用中，k根据实际情况选取。

根据到时差值分析、残差分析，同时结合微震传感器台网空间布设形态，对微震波形中的延迟波和外部异常波两类异常波形进行识别剔除，同时还能够识别确定各个传感器拾取微震波的P、S波类型，为震源定位奠定基础。

(3)单纯形微震震源定位算法(Simplex)

分别采用L1范数统计准则和L2范数统计准则计算时间残差，其中基于L1范数统计准则的事件残差计算公式为：

$E_{L1}=\frac{\underset{i=1}{\overset{n}{\Σ}}|t_i-t_{ti}-t_{median}|}{n-m}---\left(7\right)$

式(7)中t_median是所有(t_i-tt_i)的中位数。

基于L2范数统计准则的事件残差计算公式为：

$E_{L2}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\γ}}_i^2}{n-m}}---\left(8\right)$

式(7)和(8)中，n是有效微震传感器台站的个数；m表示自由度，即未知数的个数，在微震震源定位中m＝4。

单纯形微震震源定位算法(Simplex)的实施主要步骤如下：

a.将微震发震时间t₀分离出去，对震源空间坐标(x₀,y₀,z₀)进行求解，然后再计算发震时间t₀；

b.在微震监测区域内任选四点，以这四个点为顶点构造初始单纯形；

c.根据IFEP微震波形到时自动拾取结果和APSI微震波形自动识别结果，剔除异常波形信号，同时选取对应的P、S波到时和波速；

d.采用L1范数统计准则(公式(7))或L2范数统计准则(公式(8))计算单纯形各个顶点的事件残差；

e.分别比较四个点的事件残差，确定最大值和最小值；根据扩展、映射、收缩、压缩四种变换形式在误差空间中寻找新的顶点，组成新的单纯形；

f.重复以上步骤，不断移动单纯形，直到得到最佳震源空间坐标或者满足迭代终止条件时，迭代终止，同时计算微震发震时间，最终得到微震震源定位结果(x₀,y₀,z₀,t₀)。图2表示一个单纯形在误差空间中的迭代过程。

(4)微震震源定位可靠性综合评价体系(SLRE)

事件残差能够在一定程度上反映定位误差的大小，因此提出评价微震震源定位精度的事件残差指标R。采用L1范数统计准则(公式(7))计算事件残差，根据事件残差大小，得到事件残差指标R的评价标准，详见表3。

表3事件残差指标评价标准

注：表中α是事件残差修正系数。

将微震波波速增加或降低10％后，采用Simplex进行重新定位，将重定位结果与原来的定位结果之间的误差记为L，L可采用下式进行计算：

L＝S′(x₀,y0,z₀)|_v′＝v±10_％v-S(x₀,y0,z₀) (9)

上式中，M为敏感度指标，S(x₀,y₀,z₀)为原始定位结果，S′(x₀,y₀,z₀)为微震波速度增加或降低10％后，在采用相同震源定位方法的情况下，得到的重定位结果。

系统越稳定对输入数据的抵抗能力就越强，因此根据公式(9)，提出评价评价微震震源定位稳定性的敏感度指标M。根据重定位误差L的大小，确定敏感度指标M的评价标准，详见表4。

表4敏感度指标评价标准

注：表中β是敏感度修正系数。

根据观测触发序列和计算触发序列之间的错配程度，将触发序列划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ共5个等级，其中等级越高表明观测触发序列和计算触发序列之间的错配程度越小；同时为了考虑微震传感器台网对震源定位影响，结合有效传感器数量，建立能够同时评价微震震源定位精度和定位稳定性的触发序列指标T，具体评价标准如表5所示。

表5触发序列指标评价标准

综上所述，引入微震震源定位可靠性综合评价指标P，并令P＝R+M+T，从而建立微震震源定位可靠性(定位精度和稳定性)综合评价体系(SLRE)，根据SLRE将震源定位结果划分为5个等级，分别为A、B、C、D、E。SLRE的评价标准如表6所示。

表6微震震源定位可靠性综合评价体系

(6)将IFEP、APSI、Simplex、SLRE有序结合，建立微震震源自动定位及可靠性综合评价方法(AMSLRE)，AMSLRE的主要内容和流程见图3。

微震监测系统共有16个传感器，分别记为T1(x₁,y₁,z₁),T2(x₂,y₂,z₂),T3(x₃,y₃,z₃),…,T16(x₁₆,y₁₆,z₁₆)，在微震监测区域中进行已知震源位置的人工爆破实验，爆破事件记为B(x,y,z,t)，其中t为爆破时间；根据现场实测微震波P波波速为4158m/s，S波波速为2483m/s。

2016年5月18日15时37分30秒879毫秒在微震监测区域进行了一次人工爆破试验。精确测量该人工爆破事件的空间坐标为x＝2672.2,y＝4712.8,z＝-132.5，将该人工爆破事件记为B(2672.2,4712.8,-132.5,879)。

微震监测系统接收到了该爆破事件诱发的微震波形，首先采用IFEP自动拾取各个微震传感器接收的微震波形到时，图4给出了其中两种不同信噪比的现场微震波形采用IFEP自动拾取到时的实例，从图中可以看出对于低信噪比的微震波形信号(图4(b))IFEP仍然能够自动拾取高精度的波形到时。采用IFEP自动拾取的各个传感器到时详见表7。

表7 IFEP自动拾取传感器到时

注：“/”表示对应传感器没有接收到微震波形。

根据公式(3)、(4)、(5)，通过计算得到该爆破事件的到时差值分析表，详见表8。

表8爆破事件到时差值分析表

根据到时差值分析表，传感器T16和T1拾取的都是延迟波到时，因此它们不应参与定位计算。

根据公式(6)进行残差分析，本实施例中残差系数k设定为1.2。残差分析结果见表9。

表9爆破事件残差分析表

根据残差分析，将传感器T3拾取到时修正为S波到时。因此，基于到时差值和残差分析微震波形自动识别方法(APSI)得到的爆破事件微震波形识别结果如表10所示。

表10 APSI爆破事件微震波形识别结果

在IFEP自动拾取到时和APSI微震波形识别的基础上，采用单纯形微震震源定位算法(Simplex)对爆破事件B进行震源自动定位，自动定位结果如表11所示。

表11 Simplex震源自动定位结果

件采用微震震源定位可靠性综合评价体系(SLRE)对爆破事件的微震震源自动定位可靠性进行评价，其中事件残差指标R得分为8，敏感度指标M得分为-2，触发误序列指标T得分为4，综合评价指标P得分为：

P＝R+M+T＝8+(-2)+4＝10 (10)

微震震源定位可靠性综合评价结果详见表12所示。

表12 SLRE定位可靠性综合评价结果

从实施例可知，本发明的一种微震震源自动定位及可靠性综合评价方法(AMSLRE)，不仅实现了高精度的微震震源自动定位，而且对自动定位可靠性评价结果符合实际情况。

上述实施例虽然对本发明的技术方案进行了比较详细的说明，但是并不是对本发明的限制，任何对本发明的技术方案进行修改或者同等替换，而不脱离本发明技术方案的实质精神和范围，均应涵盖在本发明的权利要求范围中。