一种基于P波到时与波形振幅的震源定位方法与流程

本发明涉及一种基于P波到时与波形振幅的震源定位方法，结果明确，适用性强，准确率高。

背景技术：

微震监测作为一种地下工程安全监测的重要技术手段之一，已经广泛被应用于矿山安全、水电工程等领域。而微震震源空间位置的确定以及微震震源定位精度的提高，一直是微震监测技术的重要内容。但是，微震与爆破事件波形较为复杂，S波的首波通常夹杂在P波的尾波中，使得S波到时不仅难以辨识，且易受操作人员专业知识和经验的影响，极大降低了S波拾取准确度及数量，严重制约微震监测实时分析的效果。

震源定位方法到目前为止有十余种之多，主要包括几何方法、物理方法和数学方法等。Geiger定位方法是通过迭代计算使时间残值最小化。在此基础上，Buland对其进行修改并发展出基于经典的最小二乘法定位计算，使其能更好适应地震定位。Tatantola提出了Bayesian定位方法的严格公式和解。Spence针对时差定位技术模型的不确定导致走时异常进行了有意义的研究。Crosson应用联合测定方法(SSH)对局部事件进行了定位研究。以上定位方法，都是以预先测定平均速度或给出平均速度模型为前提，故在现场进行定位时，平均速度测量的准确性直接影响着定位精度。而且，对S波到时的自动准确拾取仍存在难题，需要依靠具有经验的专业工作人员人工拾取，定位准确性与精度存在较大的随机性。

可见现有的微震与爆破事件震源定位方法存在较大的局限，需要研究一种适用性强、准确率高的定位方法。

技术实现要素：

本发明所要解决的技术问题是提供一种基于P波到时与波形振幅的震源定位方法，该定位方法适用性强、准确率高。

发明的技术解决方案如下：

一种基于P波到时与波形振幅的震源定位方法，包括以下步骤：

步骤1：通过微震监测系统获取震源参数；

震源参数包括：各传感器坐标(x_i,y_i,z_i)，其中i为传感器编号，i＝1,2,…,n；震源起振时刻t₀，各传感器的P波到时t_Pi，各传感器的波形振幅趋势线的最大振幅MaxAmp_i及其对应时刻t_Mi；

步骤2：根据步骤1获取的初步求解震源坐标(x,y,z)及P波波速v，得到一次定位的震源坐标(x₀,y₀,z₀)和初步求解的P波波速v⁰；

步骤3：计算各传感器i(i＝1,2,…,n)的波形振幅趋势线的斜率K_i：

步骤4：基于步骤2得到的v⁰和步骤3得到所有传感器的波形振幅趋势线的斜率K_i(i＝1,2,…,n)，通过logistic函数拟合得到传感器的波形振幅趋势线的斜率K与震源的P波波速v^c之间的函数关系式K(v^c)：

其中，r为传感器的波形振幅趋势线的斜率随速度的增长速率(r是logistic函数的一个参数，经数值分析软件拟合数据之后可直接得出)，K_max为所有传感器的波形振幅趋势线的斜率中的最大值，即K_max＝max{K_i(i＝1,2,…,n)}；K₀为所有传感器的波形振幅趋势线的斜率中的最小值，即K_max＝min{K_i(i＝1,2,…,n)}；

步骤5：计算修正后的传感器i对应的P波波速

保持传感器i的波形振幅趋势线的斜率K_i不变，将K_i作为K(v^c)代入公式二，计算相应的v^c值，即为修正后的传感器i对应的P波波速

步骤6：二次定位震源坐标(x,y,z)；

将修正后的传感器i对应的P波波速代入以下目标函数：

其中，min表示求最小值；Δt_ij＝t_Pi-t_Pj；

求解目标函数，其最优解即为二次定位的震源坐标(x,y,z)。

所述步骤2具体为：

建立以下目标函数：

其中，min表示求最小值；Δt_ij＝t_Pi-t_Pj；

由此求得使得全部理论值与实测值Δt_ij的偏离平方和最小时的震源坐标(x,y,z)及P波波速v，即一次定位的震源坐标(x₀,y₀,z₀)和初步求解的P波波速v⁰。

有益效果：

本发明选取震源振动波形图中的最易准确拾取的参数：震源起振时刻t₀，各传感器P波到时t_Pi，最大振幅MaxAmp_i及其对应时刻t_Mi，优化定位所用的P波波速，推导求差式非线性拟合函数用于震源定位，事实证明其结果明确，适用性强，准确率高。

本发明基于P波到时与波形振幅对震源坐标进行确定，减轻了人工拾取S波的数据处理压力，避免了因S波拾取不准确而带来的误差，提高了定位精度与准确性。此方法具有适用性强、准确性高等特点。

附图说明

图1为本发明流程图。

具体实施方式

下面将结合实际算例，对本发明提出的一种基于P波到时与波形振幅的震源定位方法作进一步说明。本发明的思想描述如下：本发明通过微震监测系统获取主要震源参数；利用速度—距离公式和两点间距离公式得到震源与任意传感器i、j之间的等量关系，求得一次震源坐标；再通过优化P波波速提高准确性，推导求差式非线性拟合函数用于震源定位，计算出震源的准确坐标。此方法具有适用性强、准确率高的优点。

步骤1：通过微震监测系统获取震源参数；

震源参数包括：各传感器坐标(x_i,y_i,z_i)，其中i为传感器编号，i＝1,2,…,n；震源起振时刻t₀，各传感器的P波到时t_Pi，各传感器的波形振幅趋势线的最大振幅MaxAmp_i及其对应时刻t_Mi；

步骤2：根据步骤1获取一次定位的震源坐标(x,y,z)及初步求解的P波波速v，分别记为(x₀,y₀,z₀)和v⁰；；

设l_i为传感器i至震源的距离，由两点间距离公式可得：

由速度—距离公式可知，传感器i的P波到时的理论值t_i为：

即：

由式(4)可确定任意2个不同的传感器i和j的P波到时之差的理论值为：

计算传感器i和j的P波到时之差的实测值Δt_ij：

Δt_ij＝t_Pi-t_Pj(5)

全部理论值与实测Δt_ij的偏离平方和可描述全部理论值与实测值的偏离程度，其偏离越小，则相应的一次定位的震源坐标(x₀,y₀,z₀)及初步求解的P波波速v⁰与其真实值的拟合度越好，准确度越高，因此，建立以下目标函数：

其中，min表示求最小值，

由此求得使得全部理论值与实测值Δt_ij的偏离平方和最小时的震源坐标(x,y,z)及P波波速v，即一次定位的震源坐标(x₀,y₀,z₀)和初步求解的P波波速v⁰；；

步骤3：计算传感器i(i＝1,2,…,n)的波形振幅趋势线的斜率K_i：

步骤4：基于步骤2初步求解的P波波速v⁰和步骤3得到所有传感器的波形振幅趋势线的斜率K_i(i＝1,2,…,n)，通过logistic函数拟合得到传感器的波形振幅趋势线的斜率K与震源的P波波速v^c之间的函数关系式K(v^c)：

其中，r为传感器的波形振幅趋势线的斜率随速度的增长速率(r是logistic函数的一个参数，经数值分析软件拟合数据之后可直接得出)，K_max为所有传感器的波形振幅趋势线的斜率中的最大值，即K_max＝max{K_i(i＝1,2,…,n)}；K₀为所有传感器的波形振幅趋势线的斜率中的最小值，即K_max＝min{K_i(i＝1,2,…,n)}；

步骤5：计算修正后的传感器i对应的P波波速

保持传感器i的波形振幅趋势线的斜率K_i不变，将K_i作为K(v^c)代入函数关系式(8)，计算相应的v^c值，即为修正后的传感器i对应的P波波速

步骤6：再次求解震源坐标(x,y,z)；

将修正后的传感器i对应的P波波速作为v，代入目标函数(6)有：

上式必定取得非负值，因此其最小值一定存在。求解(10)，其最优解即为二次定位的震源坐标(x,y,z)。

本发明选取震源振动波形图中的最易准确拾取的参数：震源起振时刻t₀，各传感器P波到时t_Pi，最大振幅MaxAmp_i及其对应时刻t_Mi，优化定位所用的P波波速，推导求差式非线性拟合函数用于震源定位，事实证明其结果明确，适用性强，准确率高。实施实例：

由步骤2可得一次定位震源坐标如下所示：

表1 5组震源的实际坐标和一次定位坐标

由步骤3、步骤4、步骤5可得表2如下所示：

表2 5组震源的K-v函数关系式

由步骤6可求得二次震源坐标如下表所示

表3 5组震源的实际坐标和二次定位坐标

由表3可知，经速度修正后，二次定位坐标精度明显高于一次定位坐标精度。

以上所述仅为本发明的实施例而已，并不用以限制本发明，凡在本发明精神和原则之内，所作任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。