CN105549077A - 基于多级多尺度网格相似性系数计算的微震震源定位方法 - Google Patents

Abstract

本发明公开一种微震震源定位算法，属于地球物理反演算法领域，包括如下步骤：对监测区域进行均匀网格划分；分别计算每个网格分割点对应的地震振幅相似性系数；选取相似性系数最大值处对应的网格位置作为震源定位结果；在定位结果对应的网格处再次进行均匀网格划分、相似性系数计算，选取相似性系数最大值处对应的网格位置更新震源定位结果，直至定位精度满足要求。本发明相对于传统的相似性系数定位方法而言，通过对监测区域进行多级多尺度的网格划分，避免监测区域整体密集网格划分造成的计算时间过长的缺点，也避免监测区域整体稀疏网格划分导致的定位精度过低的缺点。

Description

基于多级多尺度网格相似性系数计算的微震震源定位方法

技术领域

本发明是一种基于多级多尺度网格划分的利用相似性系数的微地震定位反演方法，涉及地球物理探测技术领域，特别涉及微地震监测定位反演方法。

背景技术

地下的微小震动作为地下介质大型动力灾害的前兆，对于微震震源的实时定位，可以有效的监测地下介质的破裂情况以及对于大型动力灾害的预防提供信息。在微地震监测工作中，如何准确的对微震震源进行定位是监测工作成功的关键。微震监测因其可以实时监控介质破裂程度的特点，已大量运用在工程领域中，在煤矿系统利用微震监测判断煤层顶底板破裂情况，在工程物探中利用微震来判断桩基介质的内部情况。

目前成熟的微震震源定位算法大多是基于盖格理论的方法，通过检波器空间位置、地下介质体速度和地震波到达检波器的时间的关系建立方程组，基于最优化理论求得定位结果。基于盖格理论的定位过程一种是利用非线性反演算法中目标函数的几何性质寻求迭代方向寻求最优值。但是这种手段往往使最终结果寻求到局部的最优值导致定位结果不精确。另一种则是利用在搜索空间的随机搜寻来寻找全局的最优值来减低反演算法陷入局部最优的危险。

基于盖格类的定位方法在地震数据信噪比高的情况下可以进行准确的定位，但是微地震监测中得到的数据信噪比较低，利用盖格方法对数据进行定位处理往往得不到理想的定位结果。另外，面对大量地震数据，盖格类方法只是利用了地震波的旅行时间，丢弃了地震数据中的波形数据。

基于相似性系数计算的定位方法是一种同时利用地震波旅行时和地震波形数据计算相似性系数进行定位的方法，在低信噪比或速度模型较复杂时该方法的定位精度较之盖格类方法有所提高，在天然地震定位中得到了很好的应用。微震震源定位借鉴天然地震定位的方法，对基于相似性系数计算的定位方法开展了研究并取得了一定的应用成果。然而微震震源定位在指导资源开采时，对方法的实时性和定位精度都要求较高。为了同时达到较好的定位精度和实时性，本发明提出一种基于多级多尺度网格划分的相似性系数计算方法。首先通过对定位空间进行稀疏且均匀的网格划分得到对应的网格点，计算网格点中相似性系数最高的网格点；然后对相似性系数计算结果最大值的网格点附近的空间进行密集的网格划分，再一次对网格点进行相似性系数的计算，直至定位精度达到要求，最终得到的相似性系数值最高的空间点就可认定是最终的微震震源位置。

本算法的目的在于多级多尺度的网格划分同时提高了定位方法的计算效率和计算精度，利用相似性系数定位的方法通过利用地震波形数据进行精确的微震震源定位。此外，本方法不用拾取地震波从震源到检波器的旅行时间，避免了盖格法在地震信号信噪比不高的情况下，拾取地震波旅行时不准确造成定位不理想。

发明内容

本发明是为了解决传统基于相似性系数计算的震源定位算法计算效率和计算精度无法同时提高，应用在微震震源定位时精度和实时性无法同时满足，导致应用范围受限的问题，通过多级多尺度网格划分，较之传统均匀网格划分时相似性系数计算的定位方法，同时提高了定位方法的计算效率和计算精度，满足微震监测实时性的同时，实现震源的精确定位。

基于多级网格划分和相似性系数的微震震源定位方法的基本实现步骤如下：

步骤一：将微震监测区域均匀离散化为i行×j列的网格模型，从原点处逐行进行编号，将对应的网格分割点坐标位置存储在i行×j列的二维数组L₀中，同时建立i行×j列的二维数组L₁，数组L₁中元素的值为数组L₀中对应网格分割点的平均传播速度，L₁的编号顺序同L₀；

步骤二：根据数组L₀和L₁中的值，计算i行×j列的网格模型中每一个网格位置到每一个检波器的地震波旅行时，并存储在i行×j列×m道的三维数组L₃中，其中m道表示检波器道数；

步骤三：选取长度为N的时窗，根据数组L₃中对应的地震波旅行时对时窗进行移动，选取视窗内的地震振幅参数；

步骤四：根据时窗内选取的振幅信息，将相同行数和列数处对应的所有道的振幅进行相加，计算相似性系数，并将结果存在存储在i行×j列的二维数组L₂中的对应位置；

步骤五：选取步骤四二维数组L₂中相似性系数最大的点，以这个点对应的L₀空间坐标位置作为暂时的定位结果；

步骤六：对步骤五中得到的相似性系数最大的网格进行密集地划分，将其均匀离散化为i行×j列的网格模型，将每个网格分割点对应的空间坐标和空间平均传播速度分别存入步骤一中提到的二维数组L₀和L₁中，重复步骤二至对步骤五，直至定位精度满足要求；

步骤七：将步骤六中计算得到的相似性系数最大的网格点，对应到L₀数组对应的真实物理空间位置处，并将该位置是最终的微震震源定位结果。

进一步的，在所述步骤一中，根据监测区域的大小将其等效成长方形剖面，在对监测区域的大小和网格的大小进行设定后，以坐标轴原点处的网格为起始，编号为1，以测线方向为行，垂直测线方向为列，逐行依次编号，从而将剖面划分为i行和j列，并将对应行数和列数处网格分割点对应的空间坐标存储在数组L₀中，网格内平均速度的值存入二维数组L₁中；若监测空间是三维空间，则可将该区域等效成长方体，二维数组L₀和L₁需要扩充为三维数组。

进一步的，所述步骤二中，每一个网格划分点到每一个检波器的地震波旅行时的具体表示为：

$t_{ij}^{\left(m\right)}=\frac{1}{v_{ij}}\×\sqrt{{(x_i-x_m)}^2+{(y_j-y_m)}^2}$

其中表示三维数组L₁中第i行和第j列所对应的网格点所对应到第m道检波器的地震波旅行时；x_i,y_j表示L₀中第i行和第j列对应空间位置的坐标；x_m，y_m表示第m道检波器的坐标。v_ij为二维数组L₁对应的第i行和第j列网格点的空间物理位置的平均传播速度。

进一步的，所述步骤三中，选取长度为N的时窗，根据数组L₃中对应的地震波旅行时对时窗进行移动，选取时窗内的地震振幅参数，其中时窗的长度N至少要大于权利要求1步骤二中计算得到的数组L₁中的最大值。

进一步的，所述步骤四中，据时窗内选取的振幅信息，将相同行数和列数处对应的所有道的振幅进行相加，计算相似性系数，并将结果存在存储在i行×j列的二维数组L₂中的对应位置，其具体形式为：

$S_{ij}=\frac{{\left(\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{X}_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}{M\×\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{\left(X_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}$

其中S_ij是表示二维数组L₂中，第i行和第j列对应的网格点的相似性系数，N表示时窗的长度，M表示检波器的个数，表示对第m个地震道的数据，将时窗进行长度为的移动后，时窗内第n个点的地震数据幅值。

进一步的，所述步骤五中，选取步骤四二维数组L₂中相似性系数最大的点，以这个点对应的L₀空间坐标位置作为暂时的定位结果；若监测区域出现多个微震震源，则需要建立一个阈值，超过阈值的相似性系数则被认定为暂时的定位结果，对应的空间网格处都需要再次加密细分。

进一步的，所述步骤六中，对步骤五中得到的相似性系数最大的网格进行密集地划分，将每个网格分割点对应的空间坐标和空间平均传播速度分别存入步骤一中提到的二维数组L₀和L₁中，重复步骤二至对步骤五，直至定位精度满足要求。

进一步的，所述步骤七中，其特征在于将步骤六中计算得到的相似性系数最大的网格点，对应到数组L₀中对应的真实物理空间位置处，并将该位置作为最终的微震震源定位结果。

本发明的技术方案，首先通过对监测区域进行稀疏且均匀的网格划分，计算网格点对应的相似性系数，得到一个距离真实位置比较接近的空间位置，再以此位置为中心，以一定的范围进行密集的网格划分，再次通过计算相似性系数，选取相似性系数值最大的点作为最终的定位结果。这样，能够避免监测区域整体的密集网格划分造成的计算时间过长的缺点，也能够避免由于整体的稀疏网格划分导致的定位精度过低的缺点，最终实现高效、高精度的微震震源定位。

附图说明

图1是基于多级多尺度网格相似性系数计算的微震震源定位方法的流程框图。

具体实施方式

基于多级网格划分和相似性系数的微震震源定位方法的基本实现步骤如下：

步骤一：将微震监测区域均匀离散化为i行×j列的网格模型，将对应的网格分割点的坐标位置存储在i行×j列的二维数组L₀中，从原点处逐行进行编号，建立i行×j列的二维数组L₁，数组L₁中元素的值为对应网格空间物理位置的平均传播速度，编号顺序同L₀。

步骤二：根据数组L₁中的值，计算i行×j列的网格模型中每一个网格位置到每一个检波器的地震波旅行时，并存储在i行×j列×m道的三维数组L₃中，其中m道表示检波器道数。具体表示为：

$t_{ij}^{\left(m\right)}=\frac{1}{v_{ij}}\×\sqrt{{(x_i-x_m)}^2+{(y_j-y_m)}^2}$

其中表示三维数组L₁中，第i行和第j列所对应的网格点所对应到第m道检波器的地震波旅行时。x_i,y_j表示L₁中第i行和第j列的坐标，x_m，y_m表示第m道检波器的坐标。v_ij为二维数组L₁对应的第i行和第j列网格点的空间物理位置的平均传播速度。

步骤三：选取长度为N的时窗，根据数组L₃中对应的地震波旅行时对时窗进行移动，选取视窗内的地震振幅参数。

步骤四：根据时窗内选取的振幅信息，将相同行数和列数处对应的所有道的振幅进行相加，计算相似性系数，并将结果存在存储在i行×j列的二维数组L₂中的对应位置。其具体形式为：

$S_{ij}=\frac{{\left(\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{X}_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}{M\×\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{\left(X_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}$

其中S_ij是表示二维数组L₂中，第i行和第j列对应的网格点的相似性系数，N表示时窗的长度，M表示检波器的个数，表示对第m个地震道的数据，将时窗进行长度为的移动后，时窗内第n个点的地震数据幅值。

步骤五：选取步骤四二维数组L₂中相似性系数最大的点，以这个点对应的L₀空间坐标位置作为暂时的定位结果。

步骤六：对步骤五中得到的相似性系数最大的网格进行密集地划分，将其均匀离散化为i行×j列的网格模型，将每个网格分割点对应的空间坐标和空间平均传播速度分别存入步骤一中提到的二维数组L₀和L₁中，重复步骤二至对步骤五，直至定位精度满足要求。

步骤七：将步骤六中计算得到的相似性系数最大的网格点，对应到L₀数组对应的真实物理空间位置处，并将该位置是最终的微震震源定位结果。

下面对基于多级网格划分和相似性系数的微震震源定位方法进行验证

一、正演模型

正演模型是一个长5000米，宽500米的单层均匀介质，速度为500m/s，检波器布设在地表，分布在0米至5000米范围内，检波器间的距离为100米。设计的震源位置的坐标为(2400,1060)处。

二、反演结果

由正演模型参数可得到51道检波器的数据。对于计算相似性系数的时窗，选择长度为30，即时窗内选择30个的地震幅值数据。

第一个将监测区域进行稀疏且均匀的网格划分，具体划分网格空间位置如表1。

表1监测区域网格划分点坐标

划分网格之后，对每一个网格上的点进行对应的相似性系数的计算。得到的结果如表2。

表2网格划分点对应的相似性系数

从表2的结果，可以看出分割点14所对应的相似性系数是最大的，所以认定真实的震源位置在分割点14对应的位置附近，分割点14的空间坐标由表1可知为(2500,1000)。

以分割点14的坐标为中心，进行密集的网格划分，在分割点14横坐标200米的范围，纵坐标100米的范围进行间距为10米的网格划分。这样在一个小范围内得到了21×41个网格划分点。划分点的序号规律为从左向右，从上到下依次增加。再次进行每一个网格划分点的相似性计算。部分结果如表3。

表3部分网格划分点对应的相似性系数

由整体的网格分割点对应的相似性系数，第427个分割点对应的相似性系数是最大的为0.6725。则认为此点所对应的空间位置为最终的定位结果。第427个分割点对应的空间位置为(2400,1600)。此点正是正演时模拟的震源位置。所以验证基于多级网格划分和相似性系数的微震震源定位方法是一种有效的微震震源定位方法。

Claims (8)

1.一种基于多级多尺度网格相似性系数计算的微震震源定位方法，其特征在于利用多级多尺度网格划分的方法计算网格对应分割点的相似性系数进行精确定位，其特征包括：

步骤一：将微震监测区域均匀离散化为i行×j列的网格模型，从原点处逐行进行编号，将对应的网格分割点坐标位置存储在i行×j列的二维数组L₀中，同时建立i行×j列的二维数组L₁，数组L₁中元素的值为数组L₀中对应网格分割点的平均传播速度，L₁的编号顺序同L₀；

步骤二：根据数组L₀和L₁中的值，计算i行×j列的网格模型中每一个网格位置到每一个检波器的地震波旅行时，并存储在i行×j列×m道的三维数组L₃中，其中m道表示检波器道数；

步骤三：选取长度为N的时窗，根据数组L₃中对应的地震波旅行时对时窗进行移动，选取视窗内的地震振幅参数；

步骤四：根据时窗内选取的振幅信息，将相同行数和列数处对应的所有道的振幅进行相加，计算相似性系数，并将结果存在存储在i行×j列的二维数组L₂中的对应位置；

步骤五：选取步骤四二维数组L₂中相似性系数最大的点，以这个点对应的L₀空间坐标位置作为暂时的定位结果；

步骤六：对步骤五中得到的相似性系数最大的网格进行密集地划分，将其均匀离散化为i行×j列的网格模型，将每个网格分割点对应的空间坐标和空间平均传播速度分别存入步骤一中提到的二维数组L₀和L₁中，重复步骤二至对步骤五，直至定位精度满足要求；

步骤七：将步骤六中计算得到的相似性系数最大的网格点，对应到L₀数组对应的真实物理空间位置处，并将该位置是最终的微震震源定位结果。

2.如权利要求1所述的方法，其特征在于，所述步骤一中，其特征在于，在对监测区域的大小、网格的大小进行设定后，得到对应的划分结果，并将网格分割点对应的空间坐标存储在数组L₀中、网格对应地下介质平均速度存储在数组L₁中。

3.如权利要求1所述的方法，其特征在于，所述步骤二中，其特征在于，根据存储网格分割点对应的地下介质平均速度的数组L₁，计算每一个网格划分点到每一个检波器的地震波旅行时，具体表示为：

$t_{ij}^{\left(m\right)}=\frac{1}{v_{ij}}\×\sqrt{{(x_i-x_m)}^2+{(y_j-y_m)}^2}$

其中表示数组L₃中第i行和第j列网格分割点所对应的网格点所对应到第m道检波器的地震波旅行时；x_i,y_j表示L₀中第i行和第j列对应空间位置的坐标；x_m，y_m表示第m道检波器的坐标；v_ij为二维数组L₁对应的第i行和第j列网格点的空间物理位置的平均传播速度。

4.如权利要求1所述的方法，其特征在于，所述步骤三中，其特征在于，选取长度为N的时窗，根据数组L₃中对应的地震波旅行时对时窗进行移动，选取视窗内的地震振幅参数，时窗的长度N至少要大于权利要求1步骤二中计算得到的数组L₁中的最大值。

5.如权利要求1所述的方法，其特征在于，所述步骤四中，其特征在于，根据时窗内选取的振幅信息，将相同行数和列数处对应的所有道的振幅进行相加，计算相似性系数，并将结果存在存储在i行×j列的二维数组L₂中的对应位置，其具体形式为：

$S_{ij}=\frac{{\left(\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{X}_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}{M\×\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{M}{\Σ}}{\left(X_m\right(t_n-t_{ij}^{\left(m\right)}\left)\right)}^2}$

其中S_ij是表示二维数组L₂中，第i行和第j列对应的网格点的相似性系数，N表示时窗的长度，M表示检波器的个数，表示对第m个地震道的数据，将时窗进行长度为的移动后，时窗内第n个点的地震数据幅值。

6.如权利要求1所述的方法，其特征在于，所述步骤五中，选取步骤四二维数组L₂中相似性系数最大的点，以这个点对应的L₀空间坐标位置作为暂时的定位结果；若监测区域出现多个微震震源，则需要建立一个阈值，超过阈值的相似性系数则被认定为暂时的定位结果，对应的空间网格处都需要再次加密细分。

7.如权利要求1所述的方法，其特征在于，所述步骤六中，其特征在于对步骤五中得到的相似性系数最大的网格进行密集地划分，将每个网格分割点对应的空间坐标和空间平均传播速度分别存入步骤一中提到的二维数组L₀和L₁中，重复步骤二至对步骤五，直至定位精度满足要求。

8.如权利要求1所述的方法，其特征在于，所述步骤七中，其特征在于将步骤六中计算得到的相似性系数最大的网格点，对应到数组L₀中对应的真实物理空间位置处，并将该位置作为最终的微震震源定位结果。