CN111736208B - 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质 - Google Patents

Abstract

本发明公开了一种变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质，该方法首先将微震信号从NEZ坐标系旋转至RTZ坐标系；再分别从Z方向和T方向拾取P波和S波初至时刻；提出变权重联合P波和S波初至数据的微震定位目标函数；进而构建Bayes后验概率密度函数，并采用马尔科夫链蒙特卡洛(MCMC)采样计算模型的后验概率密度函数；判断是否接受新模型，并进入下一次迭代，直至达到设定迭代次数；以MCMC采样后期序列的均值作为微震事件位置。该方法具有震源定位约束性好、易于得到全局最优和自动调整联合权重等特点。

Description

变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质

技术领域

本发明属于微震监测领域，尤其是涉及一种变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质。

背景技术

微震监测技术是一种能够进行实时动态监测的地球物理技术，其广泛用于矿山开采和隧道工程。微震监测系统通过布设宽频带高灵敏度传感器采集微震信号，分析表征这些动力学灾害的发生时间，位置、震级和震源机制等特征参数，以推断岩体的力学状态、进而采取有效的防控措施。其中，微震事件定位能反映动力灾害发生的位置，且其是计算震级、反演震源机制和进行矿山灾害风险性评估的核心基础，由此高精度微震定位至关重要。

目前，微震事件震源定位常用基于走时的射线追踪方法，它由传感器观测的特征震相走时数据和通过射线追踪计算的理论走时的差异建立目标函数，Li et al.(2016)给出了4种常用的定位目标函数。经典的Geiger定位法(Geiger,1912)在地震学领域应用广泛，其将时差方程转换为线性方程组，再通过迭代求解震源位置，一些学者从时差方程目标函数和求解方式等方面对Geiger法进行了改进(Buland,1976；Thurber,1985)。另一种常用的时差定位方法是Waldhauser and Ellsworth(2000)提出的双差定位法，其假定两个相近的地震事件激发波场的传播路径相似，有效的降低了相近地震传播到台站时共同路径上结构异常对走时的影响，Li et al.(2013)将微分方位角信息引入至双差定位法中，进一步提高了定位精度。反演问题目标函数的收敛性常和选择的最优化反演算法息息相关，为此粒子群算法、遗传算法、梯度下降法、网格法和模拟退火算法等方法被用于求解上述定位问题的目标函数。

射线走时定位方法通常使用P波，主要借助传感器P波初至到时来确定震源位置，然而微震信号P波能量通常较弱，容易受背景噪音、邻近事件微震信号错误划分和波场传播的衰减和波前愈合效应等因素的影响而难以准确拾取，导致单一P波初至数据无法保障射线走时方法的定位精度和稳定性。而微震信号S波能量通常较强，初动振幅不易被衰减掉，Li et al.(2016)的结果显示S波微震定位也能取得较好的定位效果，但实际监测中S波初至可能受P波信号影响，拾取准确度较低。由此，如何建立一种考虑P、S波数据质量的联合方法对提高射线走时定位精度具有重要意义。

发明内容

本发明提出了一种变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质，该方法具有震源定位约束性好、易于得到全局最优和自动调整联合权重等特点。

本发明的技术方案如下：

一方面，一种变权重联合P波和S波初至数据的微震事件Bayes定位方法，包括以下步骤：

步骤1：将微震信号从NEZ坐标系旋转至RTZ坐标系；

步骤2：分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

步骤3：建立变权重联合P波和S波初至数据的微震定位目标函数；

步骤4：利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

步骤5：基于MCMC对参数模型进行参数采样，利用采样后的参数计算所述Bayes后验概率密度函数值，判断是否接受基于MCMC采样获得的参数模型，并进入下一次迭代，直至达到设定迭代次数；

步骤6：取MCMC迭代采样后期序列的震源位置参数的均值作为微震事件位置。

进一步地，所述变权重联合P波和S波初至数据的微震定位目标函数为：

其中，

和

分别为Z方向第i个P波和T方向第i个S波初至观测到时，

和

分别为第i个P波和第i个S波理论到时，

M₁和M₂分别为P波和S波初至数据的个数，w表示联合权重，w∈(0,1)。

进一步地，所述基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数如下：

其中，

为关于震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w的参数模型，即

为考虑系数的P波初至和S波初至到时观测数据依次排开所组成的一个列向量，

为考虑系数的P波初至和S波初至到时理论数据依次排开所组成的一个列向量，且

为对角线元素为

大小为(M₁+M₂)×(M₁+M₂)的矩阵，M₁和M₂分别为

包含的

和

的个数，且

为矩阵

的协方差，

为矩阵

的协方差，w∈(0,1)。

进一步地，所述基于MCMC对参数模型进行参数采样，是指利用MCMC对震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w均进行采样更参数新模型，每次随机更新上述参数中的一个参数，记更新后的参数模型为

判断是否接受基于MCMC采样获得的参数模型是指：当

时，接受

当

时，拒绝

其中，u为在[0,1]内均匀分布随机数据，在每次判断时，均重新随机生成具体取值。

在具体的求解过程中，更新参数模型后可由新参数模型的震源位置(x₀′,y₀′,z₀′)计算P波和S波理论传播时间，理论传播时间分别与新模型发震时刻t₀′相加得到当前参数模型下P波和S波理论到时

进一步地，可得到更新后的

由此得到

的表达式，获得更新参数模型后的

并判断是否接受MCMC采样新模型

当

时，接受

当

时，拒绝

u'∈[0,1]为在每次判断过程中重新生成的均匀分布的随机数据；

进一步地，计算

时，采用对

取对数进行计算，当

时，

当

时，

其中，

为参数模型更新后，对应更新

得到的矩阵。

进一步地，将微震信号从NEZ坐标系旋转至RTZ坐标系，是指按照以下公式进行转换：

其中，N、E、Z分别表示微震信号在正北方向，正东方向、垂直向上的振幅值，R为震源位置与台站连线方向的振幅值，T为垂直于R和Z方向所构成平面的振幅值，α为R方向与N方向的夹角。

α可根据专利201510973875.4《一种用于均匀速度场的信号源定位方法》初定位结果得到；

进一步地，从Z方向和T方向微震信号拾取P波和S波初至时刻，是指首先采用专利《一种矿山微震信号P波初至时刻联合拾取方法》分别从RTZ坐标系的Z方向和T方向确定P波和S波初至时刻，对该方法无初至拾取波形借助人工方法进行拾取。

进一步地，取MCMC采样后期序列5000～20000个采样点的均值作为微震反演结果，并得到微震事件位置。

随着MCMC采样迭代次数的增加，各个参数序列趋于收敛；

另一方面，一种变权重联合P波和S波初至数据的微震事件Bayes定位系统，包括：

微震信号坐标系转换单元：将微震信号从NEZ坐标系旋转至RTZ坐标系；

P波和S波初至时刻拾取单元：用于分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

目标函数建立单元：用于建立变权重联合P波和S波初至数据的微震定位目标函数；

Bayes后验概率密度函数及求解单元：用于利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

MCMC求解单元：利用MCMC对参数模型进行采样，求解Bayes后验概率目标函数，并判断是否接受新模型；

微震事件定位单元：以MCMC采样后期序列的均值作为微震事件最终定位结果。

再一方面，一种可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行上述的一种变权重联合P波和S波初至数据的微震事件Bayes定位方法。

有益效果

本发明提供了一种变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质，主要用于解决单一震相初至数据较少使得微震定位约束有限、单一震相初至数据总可能含有大拾取误差导致微震定位稳定性不足的问题。本方法包括如下步骤：首先将微震信号从NEZ坐标系旋转至RTZ坐标系；再分别从Z方向和T方向拾取P波和S波初至时刻；提出变权重联合P波和S波初至数据的微震定位目标函数；进而构建Bayes后验概率密度函数，并采用MCMC采样计算模型的后验概率密度函数；判断是否接受新模型，并进入下一次迭代，直至达到设定迭代次数；以MCMC采样后期序列的均值作为微震事件位置。该方法充分利用了微震信号的P波和S波初至数据，增强了微震事件定位的约束性。Bayes方法以一定概率接受新模型，可跳出局部最优解，大大降低初始迭代值的影响。此外，变权重方式能根据P波和S波数据的质量自动调整P波和S波部分的权重，从而P波和S波数据中大误差拾取对微震事件定位精度的影响。

附图说明

图1是本发明所述方法流程图；

图2是微震信号旋转至RTZ坐标系的波形图；

图3是测试事件及各传感器位置图；

图4是不同数据集和数据质量下Bayes定位MCMC迭代过程图，其中，(a)为对各P和S波初至理论数据加上-1ms～1ms均匀分布随机误差后MCMC迭代过程图；(b)为在(a)之上对2个P波初至数据加上100ms大误差的MCMC迭代过程图；(c)为在(a)之上对2个S波初至数据加上100ms大误差的MCMC迭代过程图；(d)为仅使用P波初至数据，且对各P初至理论数据加上-1ms～1ms均匀分布随机误差后MCMC迭代过程图；(e)为在(d)之上对2个P波初至数据加上100ms大误差的MCMC迭代过程图。

具体实施方式

下面将结合附图1～4，对本发明做进一步地说明。

本发明所述的方法思想如下：针对常用的P波拾取方法定位结果约束不足、可能受大拾取误差影响等问题，提出一种变权重联合P波和S波初至数据的微震事件Bayes定位方法。该方法利用微震信号的P波和S波数据，增强微震事件定位的约束性，借助Bayes方法和MCMC采样求解函数目标值，并以一定概率接受新模型，跳出局部最优解，大大降低初始迭代值的影响，进而利用变权重方法根据P波和S波数据的质量自动调整P波和S波部分的权重，从而P波和S波数据中大误差拾取对微震事件定位精度的影响。

如图1所示，一种变权重联合P波和S波初至数据的微震事件Bayes定位方法，包括以下步骤：

步骤1：将微震信号从NEZ坐标系旋转至RTZ坐标系；

将微震信号从NEZ坐标系旋转至RTZ坐标系，是指按照以下公式进行转换：

其中，N、E、Z分别表示微震信号在正北方向，正东方向、垂直向上的振幅值，R为震源位置与台站连线方向的振幅值，T为垂直于R和Z方向所构成平面的振幅值，α为R方向与N方向的夹角。

α可根据专利201510973875.4《一种用于均匀速度场的信号源定位方法》初定位结果得到；

步骤2：分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

从Z方向和T方向微震信号拾取P波和S波初至时刻，是指首先采用专利《一种矿山微震信号P波初至时刻联合拾取方法》分别从RTZ坐标系的Z方向和T方向确定P波和S波初至时刻，对该方法无初至拾取波形借助人工方法进行拾取。

步骤3：建立变权重联合P波和S波初至数据的微震定位目标函数；

所述变权重联合P波和S波初至数据的微震定位目标函数为：

其中，

和

分别为Z方向第i个P波和T方向第i个S波初至观测到时，

和

分别为第i个P波和第i个S波理论到时，

M₁和M₂分别为P波和S波初至数据的个数，w表示联合权重，w∈(0,1)。

步骤4：利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

所述基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数如下：

其中，

为关于震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w的参数模型，即

为考虑系数的P波初至和S波初至到时观测数据依次排开所组成的一个列向量，

为考虑系数的P波初至和S波初至到时理论数据依次排开所组成的一个列向量，且

为对角线元素为

大小为(M₁+M₂)×(M₁+M₂)的矩阵，M₁和M₂分别为

包含的

和

的个数，且

为矩阵

的协方差，

为矩阵

的协方差，w∈(0,1)。

步骤5：基于MCMC对参数模型进行参数采样，利用采样后的参数计算所述Bayes后验概率密度函数值，判断是否接受基于MCMC采样获得的参数模型，并进入下一次迭代，直至达到设定迭代次数；

所述基于MCMC对参数模型进行参数采样，是指利用MCMC对震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w均进行采样更参数新模型，每次随机更新上述参数中的一个参数，记更新后的参数模型为

判断是否接受基于MCMC采样获得的参数模型是指：当

时，接受

当

时，拒绝

其中，u为在[0,1]内均匀分布随机数据，在每次判断时，均重新随机生成具体取值。

在具体的求解过程中，更新参数模型后可由新参数模型的震源位置(x₀′,y₀′,z₀′)计算P波和S波理论传播时间，理论传播时间分别与新模型发震时刻t₀′相加得到当前参数模型下P波和S波理论到时

进一步地，可得到更新后的

由此得到

的表达式，获得更新参数模型后的

并判断是否接受MCMC采样新模型

当

时，接受

当

时，拒绝

u'∈[0,1]为在每次判断过程中重新生成的均匀分布的随机数据；

计算

时，采用对

取对数进行计算，当

时，

当

时，

其中，

为参数模型更新后，对应更新

得到的矩阵。

步骤6：取MCMC迭代采样后期序列的震源位置参数的均值作为微震事件位置。

随着MCMC采样迭代次数的增加，各个参数序列趋于收敛，取MCMC采样后期序列5000～20000个采样点的均值作为微震反演结果，并得到微震事件位置。

实施例

图2是微震信号从NEZ坐标系旋转至RTZ坐标系的波形图。由图2知，P波初至在Z方向最为清晰，而S波在R和T方向更为清晰，且T方向的SH波起振更为明显。由此，T方向的S波初至和Z方向的P波初至数据被用于Bayes联合定位。

图3是测试事件及传感器位置图。图中三角形代表传感器；五角星为爆破试验测试事件位置，其(X,Y,Z)坐标为(381400,2997000,1000)m。设定P波传播速度为5600m/s，S波的传播速度为3000m/s，由传播距离和速度易得各传感器接收信号的P波和S波传播时间。

图4是不同数据集和数据质量下Bayes定位MCMC迭代过程图。其中，(a)为对各P和S波初至理论数据加上-1ms～1ms均匀分布随机误差后MCMC迭代过程图；(b)为在(a)之上对2个P波初至数据加上100ms大误差的MCMC迭代过程图；(c)为在(a)之上对2个S波初至数据加上100ms大误差的MCMC迭代过程图；(d)为仅使用P波初至数据，且对各P初至理论数据加上-1ms～1ms均匀分布随机误差后MCMC迭代过程图；(e)为在(d)之上对2个P波初至数据加上100ms大误差的MCMC迭代过程图。其中，单一P波初至数据时，其权重设为固定值1。由图可以得出：无大拾取误差时，单一P波初至数据、P波和S波联合数据Bayes反演方法定位稳定性很好，且联合方法优于单一特征，表明联合方法能更好地约束震源位置；对P波初至数据加了大误差后，单一P波初至数据Bayes方法定位不收敛，而P波和S波联合数据Bayes反演结果仍然很好，证明了变权重方式能根据P波和S波初至数据质量来自动调整联合权重，提高定位精度。

一种变权重联合P波和S波初至数据的微震事件Bayes定位系统，包括：

微震信号坐标系转换单元：将微震信号从NEZ坐标系旋转至RTZ坐标系；

P波和S波初至时刻拾取单元：用于分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

目标函数建立单元：用于建立变权重联合P波和S波初至数据的微震定位目标函数；

Bayes后验概率密度函数及求解单元：用于利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

MCMC求解单元：利用MCMC对参数模型进行采样，求解Bayes后验概率目标函数，并判断是否接受新模型；

微震事件定位单元：以MCMC采样后期序列的均值作为微震事件最终定位结果。

应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。

一种可读存储介质，包括计算机程序指令，其特征在于：所述计算机程序指令被处理终端执行时使所述处理终端执行一种变权重联合P波和S波初至数据的微震事件Bayes定位方法。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详尽的说明，所属领域的普通技术人员应当理解，上述实施例仅仅是对本发明的示意性实现方式的解释，实施例中的细节并不构成对本发明范围的限制，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案的等效变换、简单替换等显而易见的改变，均落在本发明保护范围之内。

Claims (9)

1.一种变权重联合P波和S波初至数据的微震事件Bayes定位方法，其特征在于，包括以下步骤：

步骤1：将微震信号从NEZ坐标系旋转至RTZ坐标系；

步骤2：分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

步骤3：建立变权重联合P波和S波初至数据的微震定位目标函数；

步骤4：利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

步骤5：基于MCMC对参数模型进行参数采样，利用采样后的参数计算所述Bayes后验概率密度函数值，判断是否接受基于MCMC采样获得的参数模型，并进入下一次迭代，直至达到设定迭代次数；

步骤6：取MCMC迭代采样后期序列的震源位置参数的均值作为微震事件位置。

2.根据权利要求1所述的方法，其特征在于，所述变权重联合P波和S波初至数据的微震定位目标函数为：

其中，

和

分别为Z方向第i个P波和T方向第i个S波初至观测到时，

和

分别为第i个P波和第i个S波理论到时，

M₁和M₂分别为P波和S波初至数据的个数，w表示联合权重，w∈(0,1)。

3.根据权利要求1所述的方法，其特征在于，所述基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数如下：

其中，

为关于震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w的参数模型，即

为考虑系数的P波初至和S波初至到时观测数据依次排开所组成的一个列向量，

为考虑系数的P波初至和S波初至到时理论数据依次排开所组成的一个列向量，且

为对角线元素为

大小为(M₁+M₂)×(M₁+M₂)的矩阵，M₁和M₂分别为

包含的

和

的个数，且

为矩阵

的协方差，

为矩阵

的协方差，w∈(0,1)。

4.根据权利要求1所述的方法，其特征在于，所述基于MCMC对参数模型进行参数采样，是指利用MCMC对震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w均进行采样更参数新模型，每次随机更新上述参数中的一个参数，记更新后的参数模型为

判断是否接受基于MCMC采样获得的参数模型是指：当

时，接受

当

时，拒绝

其中，u为在[0,1]内均匀分布随机数据，在每次判断时，均重新随机生成具体取值。

5.根据权利要求4所述的方法，其特征在于，计算

时，采用对

取对数进行计算，当

时，

当

时，

6.根据权利要求1所述的方法，其特征在于，将微震信号从NEZ坐标系旋转至RTZ坐标系，是指按照以下公式进行转换：

其中，N、E、Z分别表示微震信号在正北方向，正东方向、垂直向上的振幅值，R为震源位置与台站连线方向的振幅值，T为垂直于R和Z方向所构成平面的振幅值，α为R方向与N方向的夹角。

7.根据权利要求1所述的方法，其特征在于，取MCMC采样后期序列5000～20000个采样点的均值作为微震反演结果，并得到微震事件位置。

8.一种变权重联合P波和S波初至数据的微震事件Bayes定位系统，其特征在于，包括：

微震信号坐标系转换单元：将微震信号从NEZ坐标系旋转至RTZ坐标系；

P波和S波初至时刻拾取单元：用于分别从转换后的RTZ坐标系的Z方向和T方向拾取P波和S波初至时刻；

目标函数建立单元：用于建立变权重联合P波和S波初至数据的微震定位目标函数；

Bayes后验概率密度函数及求解单元：用于利用震源位置(x₀,y₀,z₀)、发震时刻t₀和联合权重w组建参数模型，并基于参数模型构建基于P波初至和S波初至的观测数据和理论数据之差的Bayes后验概率密度函数，求解微震定位目标函数；

MCMC求解单元：利用MCMC对参数模型进行采样，求解Bayes后验概率目标函数，并判断是否接受新模型；

微震事件定位单元：以MCMC采样后期序列的均值作为微震事件最终定位结果。

9.一种可读存储介质，包括计算机程序指令，其特征在于：所述计算机程序指令被处理终端执行时使所述处理终端执行权利要求1-7任一项所述的一种变权重联合P波和S波初至数据的微震事件Bayes定位方法。

Images