CN105022091A - 一种无预测速的远场震源快速定位方法 - Google Patents

Abstract

本发明的目的是提供一种无预测速的远场震源快速定位方法，本发明在经典TDOA定位方法的基础上，提出一种新的数学模型，在波速未知的情况下，使用降维的方式，用双曲线的渐近线构建寻优目标函数，得到最优解可以快速稳定的收敛，从而快速得到震源的坐标位置。

Description

一种无预测速的远场震源快速定位方法

技术领域

本发明涉及工业监测定位领域，尤其是一种适用于均匀介质环境下远场微震震源快速定位方法，适用于煤矿、石油井下等以震动波、声波、电磁波作传输媒介的工业监测定位系统。

背景技术

微震监测技术是煤矿灾害事故预测的重要手段。研究发现震动活动的频次、能量分布与井下冲击地压、煤与瓦斯突出等矿井动力灾害关系密切。煤矿微震监测是将微震监测技术应用在煤矿井下，记录震动活动频率、能量、波形等信息，利用震波记录初至时间信息定位煤矿井下微震事件的位置、反演煤岩体的纵波速度、分析应力的分布特征等，为预测冲击危险区域提供支持。

微地震监测的基本做法是通过在井中或地面布置检波器排列接收生产活动所产生或诱导的微小地震事件,并通过对这些事件的反演求取微地震震源位置等参数,然后,应用这些参数对生产活动进行监控或指导。长期以来，对微震震源定位方法、提高对微震源定位的准确性和精度的研究，一直是微震监测技术研究的重要内容。震源定位方法很多，主要包括几何方法、物理方法与数学方法等。目前国内外使用的震动源定位方法，大多是以预先测定平均速度或给出平均速度模型为前提，故在现场进行定位时，平均速度测量的准确性直接影响着定位精度。

目前最为常用无需预先测速的震动源定位的以二维函数做寻优目标的到达时间差定位方法，但是这个方法的求解速度和效率都不够，研究发现震源位于远场时，二维目标函数解值在收敛域附近会发生持续震荡，严重影响了求解的稳定性和速度。

因此，有必要研究出一种无预测速的远场震源快速定位方法，从而解决现有技术当中的上述问题。

发明内容

本发明的目的是提供一种无预测速的远场震源快速定位方法，从而能够快速确定震源的位置。

本发明在经典TDOA定位方法的基础上，提出一种新的数学模型，在波速未知的情况下，使用降维的方式，用双曲线的渐近线构建寻优目标函数，得到最优解可以快速稳定的收敛，从而快速得到震源的坐标位置。

本发明请求保护一种无预测速的远场震源快速定位方法，该方法包括如下步骤：

S101，获取N个检波器的位置坐标(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、...(x_n,y_n)，该N个检波器(N≥4)成等距线型排列在待检测定位的环境中，如图2所示；

S102,在待测区域发生微震活动时，各检波器将接收到震动波信号，连接各检波器的采集处理模块将震动信号集中处理，得到各相邻检波器震波的到达时间差，分别记为t₁₂,t₂₃…,t_(N-1)N；

S103，将每3个检波器定义为一组检波器，如(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为第一组检波器，对应到时差为t₁₂、t₂₃，得待测位置横坐标与波速一组关系式x＝f(v)；

$x=\frac{\frac{x_2+x_3}{2}\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2}-\frac{x_1+x_2}{2}\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2-\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}}---\left(1\right)$

认为(x₄,y₄)、(x₅,y₅)、(x₆,y₆)为第二组检波器，对应到时差为t₄₅、t₅₆，根据(1)式得另一组横坐标与波速关系式x′＝f(v)；

$X^{,}=\frac{\frac{x_5+x_6}{2}\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\frac{x_4+x_5}{2}\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}$

对每两组检波器求分别求|x-x′|的值；

使用遗传算法寻优策略，设定目标函数为所有|x-x′|的和，即为

$\begin{array}{c}\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}{\left|x-x^{\′}\right|}_i=\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}\frac{\frac{x_{6i-4}+x_{6i-3}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\frac{x_{6i-4}+x_{6i-5}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}\\ -\frac{\frac{x_{6i}+x_{6i-1}}{2}\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\frac{x_{6i-1}+x_{6i-2}}{2}\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}\end{array}---\left(2\right)$

找到使目标函数为0的v即可认为求解过程完成。

S104，将求得的v和各已知参数代入式(1)和渐近线公式(3)即可求得目标坐标。

$y=\frac{y_1+y_2}{2}\±\frac{\left(x-\frac{x_1+x_2}{2}\right)\sqrt{{(x_2-x_1)}^2-t_{12}^2{v}^2}}{t_{12}v}---\left(3\right)$

附图说明

图1为本发明方法的流程示意图；

图2为实施例中检波器位置的示意图。

具体实施方式

下面通过具体实施例对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。

以下通过图1对本发明的方法作进一步的介绍。本发明提供一种无预测速的远场震源快速定位方法，该方法包括如下步骤：

S101，获取N个检波器的位置坐标(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、...(x_n,y_n)，该N个检波器(N≥4)成等距线型排列在待检测定位的环境中，如图2所示；

S102,在待测区域发生微震活动时，各检波器将接收到震动波信号，连接各检波器的采集处理模块将震动信号集中处理，得到各相邻检波器震波的到达时间差，分别记为t₁₂,t₂₃…,t_(N-1)N；

S103，将每3个检波器定义为一组检波器，如(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为第一组检波器，对应到时差为t₁₂、t₂₃，得待测位置横坐标与波速一组关系式x＝f(v)；

$x=\frac{\frac{x_2+x_3}{2}\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2}-\frac{x_1+x_2}{2}\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2-\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}}---\left(1\right)$

认为(x₄,y₄)、(x₅,y₅)、(x₆,y₆)为第二组检波器，对应到时差为t₄₅、t₅₆，根据(1)式得另一组横坐标与波速关系式x′＝f(v)；

$X^{,}=\frac{\frac{x_5+x_6}{2}\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\frac{x_4+x_5}{2}\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}$

对每两组检波器求分别求|x-x′|的值；

使用遗传算法寻优策略，设定目标函数为所有|x-x′|的和，即为

$\begin{array}{c}\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}{\left|x-x^{\′}\right|}_i=\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}\frac{\frac{x_{6i-4}+x_{6i-3}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\frac{x_{6i-4}+x_{6i-5}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}\\ -\frac{\frac{x_{6i}+x_{6i-1}}{2}\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\frac{x_{6i-1}+x_{6i-2}}{2}\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}\end{array}---\left(2\right)$

找到使目标函数为0的波速v即可认为求解过程完成。

S104，纵坐标y的计算公式为 $y=\frac{y_1+y_2}{2}\±\frac{\left(x-\frac{x_1+x_2}{2}\right)\sqrt{{(x_2-x_1)}^2-t_{12}^2{v}^2}}{t_{12}v}---\left(3\right)$

将求得的v和各已知参数代入公式(1)和公式(3)即可求得目标坐标(x，y)。

实施例一：

假设检波器坐标为(970,0)、(980,0)、(990,0)、(1000,0)，震源坐标为(600,1000)，波速为3000m/s，通过数值计算可得，各检波器之间的到时差t₂₁＝1.170400ms、t₃₂＝1.197628ms、t₄₃＝1.224583ms。以本例对二维定位问题进行详细说明，实际定位时，已知量为四个检波器的坐标及传感器触发记录的时刻，微震源的位置未知，在这里之所以给出，目的是通过本发明提出的方法进行验证。具体实施步骤如下：

1、列出四个检波器的坐标值：(970,0)、(980,0)、(990,0)、(1000,0)；

2、列出经采集处理后得到的各检波器间的到到时差t₂₁＝1.170400ms、t₃₂＝1.197628ms、t₄₃＝1.224583ms；

3、以标准GA作为寻优方法，以式(2)|x-x′|作为目标函数进行计算。算法运行于MATLAB环境，使用英国谢菲尔德大学开发的GATBX工具箱编程实现，取v∈(1,8165)；

通过分析计算的结果，可以发现本发明的方法收敛速度较快且稳定，一般在迭代20次左右即收敛，并且计算最优的v非常稳定，其值恒定为3002.4，已经较接近真实解。

4、将得到稳定的v值代入式(1)、(3)，即可得到待测震源位置坐标，计算得到坐标为(600.1，998.8)，近似于真实值。

应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (1)

1.一种无预测速的远场震源快速定位方法，其特征在于，该方法包括如下步骤：

S101，获取N个检波器的位置坐标(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、...(x_n,y_n)，该N个检波器(N≥4)成等距线型排列在待检测定位的环境中；

S102,在待测区域发生微震活动时，各检波器将接收到震动波信号，连接各检波器的采集处理模块将震动信号集中处理，得到各相邻检波器震波的到达时间差，分别记为t₁₂,t₂₃…,t_(N-1)N。

S103，将每3个检波器定义为一组检波器，如(x₁,y₁)、(x₂,y₂)、(x₃,y₃)为第一组检波器，对应到时差为t₁₂、t₂₃，得待测位置横坐标与波速一组关系式x＝f(v)；

$x=\frac{\frac{x_2+x_3}{2}\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2}-\frac{x_1+x_2}{2}\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_3-x_2}{t_{23}}\right)}^2-v^2-\sqrt{{\left(\frac{x_2-x_1}{t_{12}}\right)}^2-v^2}}}---\left(1\right)$

定义(x₄,y₄)、(x₅,y₅)、(x₆,y₆)为第二组检波器，对应到时差为t₄₅、t₅₆，根据公式(1)得另一组横坐标与波速关系式x′＝f(v)；

$x^{,}=\frac{\frac{x_5+x_6}{2}\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\frac{x_4+x_5}{2}\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_6-x_5}{t_{56}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_5-x_4}{t_{45}}\right)}^2-v^2}}$

对每两组检波器求分别求|x-x′|的值；

使用遗传算法寻优策略，设定目标函数为所有|x-x′|的和，即为

$\begin{array}{c}\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}{\left|x-x^{\′}\right|}_i=\underset{i=1}{\overset{i=\frac{N}{6}}{\Σ}}\frac{\frac{x_{6i-4}+x_{6i-3}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\frac{x_{6i-4}+x_{6i-5}}{2}\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i-4}-x_{6i-3}}{t_{6i-4,6i-3}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-4}-x_{6i-5}}{t_{6i-4,6i-5}}\right)}^2-v^2}}\\ -\frac{\frac{x_{6i}+x_{6i-1}}{2}\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\frac{x_{6i-1}+x_{6i-2}}{2}\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}{\sqrt{{\left(\frac{x_{6i}-x_{6i-1}}{t_{6i,6i-1}}\right)}^2-v^2}-\sqrt{{\left(\frac{x_{6i-1}-x_{6i-2}}{t_{6i,6i-2}}\right)}^2-v^2}}\end{array}---\left(2\right)$

找到使目标函数为0的波速v即可认为求解过程完成。

S104，纵坐标y的计算公式为 $y=\frac{y_1+y_2}{2}\±\frac{\left(x-\frac{x_1+x_2}{2}\right)\sqrt{{(x_2-x_1)}^2-t_{12}^2{v}^2}}{t_{12}v}---\left(3\right)$

将求得的v和各已知参数代入公式(1)和公式(3)即可求得震源的坐标(x，y)。