CN111257928B

CN111257928B - 一种震源的定位方法及定位系统

Info

Publication number: CN111257928B
Application number: CN202010088520.8A
Authority: CN
Inventors: 缪思钰; 张海江; 巫芙蓉
Original assignee: University of Science and Technology of China USTC
Current assignee: University of Science and Technology of China USTC
Priority date: 2020-02-12
Filing date: 2020-02-12
Publication date: 2021-08-13
Anticipated expiration: 2040-02-12
Also published as: CN111257928A

Abstract

本申请公开了一种震源的定位方法及定位系统，其中，所述震源的定位方法在震源定位过程中通过台站检测地震波的波形序列，无需人工手动拾取地震波走时，降低了震源的定位方法的人工成本，从而降低了震源的定位方法的整体成本；此外，所述震源的定位方法采用了全新的分组成像条件，该分组成像条件将传统成像条件中的相加变为分组相乘，并加入距离、质量因子控制提高每组成像质量，使得所述震源的定位方法有效地在空间上和时间上对震源进行成像，提高了震源的定位方法的分辨率，有助于计算储层体积，研究裂缝分布和评价压裂效果。

Description

一种震源的定位方法及定位系统

技术领域

本申请涉及震源探测技术领域，更具体地说，涉及一种震源的定位方法及定位系统。

背景技术

微地震(Mine Tremor or Micro Seismic)是一种小型的地震，多由于地下资源的开采过程诱发，多可定义为在开采坑道附近的岩体内因应力场变化导致岩石破坏而引起的那些地震事件。

目前，微地震检测广泛应用于非常规油气田开采和地热资源开发过程。在这些过程中，微地震事件的定位可用于刻画断裂分布，估计储层容量，指导下级压裂以及识别活动断层。因此，微地震事件的检测和定位是微地震监测中十分重要的环节。

传统的地震定位方法无法实时定位，并且无法适应数据量越来越大的趋势，另外由于微地震事件的信噪比低，传统的地震定位方法手动拾取地震波中的P波(PrimaryWave，纵波或称胀缩波)和S波(Transverse Waves，横波或称剪切波)到时的方法定位准确度低。因此基于偏移的定位方法应运而生。这类方法通常需要首先设置多个台站拾取地震波，然后根据不同台站拾取的地震波走时差校正波形，最后根据不同成像条件对震源进行定位，但现有技术中基于偏移的定位方法仍然存在分辨率较低的问题。

发明内容

为解决上述技术问题，本申请提供了一种震源的定位方法及定位系统，以实现在无需手动拾取地震波走时的前提下，提高震源的定位方法的分辨率的目的。

为实现上述技术目的，本申请实施例提供了如下技术方案：

一种震源的定位方法，包括：

将多个台站进行分组，获得多组检测组合；

获取每个所述台站检测的地震波的波形序列；

计算每个所述波形序列的特征函数，所述特征函数包含所述波形序列的振幅信息和相位信息；

将每组所述检测组合中的台站检测的波形序列的特征函数分别代入分组成像条件中，计算获得每组所述检测组合的叠加值；

根据多组所述检测组合的叠加值，确定震源位置；

所述分组成像条件为：

其中，I(τ,s)表示所述检测组合的叠加值，m表示多组所述检测组合的个数，l_s表示地震波中S波的叠加窗口长度，l_p表示地震波中P波的叠加窗口长度，E^S ₍₎表示所述特征函数中S波的波包，E^p ₍₎表示所述特征函数中P波的波包，t_s表示地震波中S波从震源到台站的走时，t_p表示地震波中P波从震源到台站的走时，W表示所述波形序列的质量因子，D表示台站的距离因子，α表示台站和该台站所在检测组合的几何中心的距离倒数的指数。

可选的，所述将多个台站进行分组包括：

按照相邻所述台站之间相距的距离，将多个台站进行分组。

可选的，所述计算每个所述波形序列的特征函数包括：

根据第一预设公式，计算每个所述波形序列的特征函数；

所述第一预设公式为：

其中，C(i)表示第i个波形序列的特征函数，x(i)表示第i个波形序列，x(i-1)表示第i-1个波形序列。

可选的，所述根据多组所述检测组合的叠加值，确定震源位置包括：

比较多组所述检测组合的叠加值，将多组所述叠加值中的最大值对应的检测组合的所在区域确定为震源位置。

一种震源的定位系统，包括：

台站分组模块，用于将多个台站进行分组，获得多组检测组合；

波形获取模块，用于获取每个所述台站检测的地震波的波形序列；

特征计算模块，用于计算每个所述波形序列的特征函数，所述特征函数包含所述波形序列的振幅信息和相位信息；

叠加值计算模块，用于将每组所述检测组合中的台站检测的波形序列的特征函数分别代入分组成像条件中，计算获得每组所述检测组合的叠加值；

震源定位模块，用于根据多组所述检测组合的叠加值，确定震源位置；

所述分组成像条件为：

可选的，所述台站分组模块具体用于，按照相邻所述台站之间相距的距离，将多个台站进行分组。

可选的，所述特征计算模块具体用于，根据第一预设公式，计算每个所述波形序列的特征函数；

所述第一预设公式为：

可选的，所述震源定位模块具体用于，比较多组所述检测组合的叠加值，将多组所述叠加值中的最大值对应的检测组合的所在区域确定为震源位置。

从上述技术方案可以看出，本申请实施例提供了一种震源的定位方法及定位系统，其中，所述震源的定位方法在震源定位过程中通过台站检测地震波的波形序列，无需人工手动拾取地震波走时，降低了震源的定位方法的人工成本，从而降低了震源的定位方法的整体成本；此外，所述震源的定位方法采用了全新的分组成像条件，该分组成像条件将传统成像条件中的相加变为分组相乘，并加入距离、质量因子控制提高每组成像质量，使得所述震源的定位方法有效地在空间上和时间上对震源进行成像，提高了震源的定位方法的分辨率，有助于计算储层体积，研究裂缝分布和评价压裂效果。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1为本申请的一个实施例提供的一种震源的定位方法的流程示意图；

图2为本申请的另一个实施例提供的一种震源的定位方法的流程示意图；

图3为本申请的又一个实施例提供的一种震源的定位方法的流程示意图；

图4为合成测试观测系统的示意图；

图5为合成测试中的P波、S波速度模型示意图；

图6(a)-图6(c)为合成波形分别加上(a)10％噪声，(b)50％噪声，(c)200％噪声后的示意图；

图7(a)-图7(c)为不同信噪比(SNR)下基于四种成像条件的震源成像结果示意图：(a)SNR＝10，(b)SNR＝2，(c)SNR＝0.5；

图8(a)-图8(d)为不同信噪比(SNR)下基于四种成像条件的震源发震时刻反演结果示意图：(a)基于线性叠加成像条件，(b)基于波包叠加成像条件，(c)基于STA/LTA叠加成像条件，(d)基于本申请实施例提供的分组成像条件；

图9为空间域分辨率测试定位结果示意图；

图10为时间域分辨率测试发震时刻反演结果示意图；

图11观测系统(三角)和检测出事件的位置分布(圆点)示意图；

图12为一个典型事件的垂直分量波形，采样率为0.5ms；

图13为四个成像条件采取的原始波形及其不同特征函数示意图；

图14为基于四种成像条件的定位结果示意图；

图15(a)-图15(d)为选取的30个事件基于走时反演的定位结果(和基于全波形偏移定位结果的对比示意图。

具体实施方式

正如背景技术中所述，传统的震源定位方法需要大量人力时间手动拾取地震波中P波到时和S波到时，从而无法实时定位，并且无法适应数据量越来越大的趋势。另一方面，微地震事件信噪比低，很难拾取准确到时，进而影响定位准确性。针对这种情况，基于偏移的定位方法应运而生。这类方法则不需要人工事先拾取到时信息，而是以不同台站走时差校正波形后根据不同成像条件对震源进行定位，根据成像条件算出的最大值点即为震源位置。与传统走时定位方法相比，偏移地震定位方法在节约了大量时间和精力的同时，由于用了除到时之外的波形信息在一定程度上提高了定位精度。

偏移定位方法是通过“叠加”来提高信噪比，从而可以检测出一些微弱地震，并且得到了广泛的应用。这类方法最具有代表性的就是震源扫描算法(SSA)(Kao and Shan,2004)，SSA是同时搜索发震事件和震源位置，叠加波形绝对振幅的成像条件对事件震源在时间空间上成像。Duncan(2005)用叠加原始波形的成像条件，首次将偏移定位方法应用于油气田的微震事件中。为了去除震源极性对波形叠加的相消作用，Kao和Shan(2007)采用了叠加波包的成像条件。同样的，叠加波形STA/LTA的成像条件(Grigoli et al.,2013)联合震源扫描算法(jSSA)(Liang et al.,2016)都可以消除震源极性影响。

但上述基于偏移的定位方法仍然存在分辨率较低的问题。

有鉴于此，本申请实施例提供了一种震源的定位方法，包括：

将多个台站进行分组，获得多组检测组合；

获取每个所述台站检测的地震波的波形序列；

根据多组所述检测组合的叠加值，确定震源位置；

所述分组成像条件为：

其中，I(τ,s)表示所述检测组合的叠加值，m表示多组所述检测组合的个数，l_s表示地震波中S波的叠加窗口长度，l_p表示地震波中P波的叠加窗口长度，E^S ₍₎表示所述特征函数中S波的波包，E^p ₍₎表示所述特征函数中P波的波包，t_s表示地震波中S波从震源到台站的走时，t_p表示地震波中P波从震源到台站的走时，W表示所述波形序列的质量因子，D表示台站的距离因子。

所述震源的定位方法在震源定位过程中通过台站检测地震波的波形序列，无需人工手动拾取地震波走时，降低了震源的定位方法的人工成本，从而降低了震源的定位方法的整体成本；此外，所述震源的定位方法采用了全新的分组成像条件，该分组成像条件将传统成像条件中的相加变为分组相乘，并加入距离、质量因子控制提高每组成像质量，使得所述震源的定位方法有效地在空间上和时间上对震源进行成像，提高了震源的定位方法的分辨率，有助于计算储层体积，研究裂缝分布和评价压裂效果。

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

本申请实施例提供了一种震源的定位方法，如图1所示，包括：

S101：将多个台站进行分组，获得多组检测组合；

可选的，参考图2，所述将多个台站进行分组包括：

S1011：按照相邻所述台站之间相距的距离，将多个台站进行分组。

即可选的，当相邻所述台站之间相距的距离小于预设距离值时，将这些台站分为一组，作为一组检测组合。

另外，在本申请的其他实施例中，还可以根据多个台站的设置顺序或序号等参数进行分组，但通常情况下，采取步骤S1011的方式进行分组有利于提高最终确定的震源位置的精确度。

S102：获取每个所述台站检测的地震波的波形序列；

S103：计算每个所述波形序列的特征函数，所述特征函数包含所述波形序列的振幅信息和相位信息；

可选的，参考图3，所述计算每个所述波形序列的特征函数包括：

S1031：根据第一预设公式，计算每个所述波形序列的特征函数；

所述第一预设公式为：

由第一预设公式可知，所述特征函数不仅包含所述波形序列的振幅信息和相位信息，而且所述特征函数都是正值，有效地避免了震源辐射的不良影响。

S104：将每组所述检测组合中的台站检测的波形序列的特征函数分别代入分组成像条件中，计算获得每组所述检测组合的叠加值；

S105：根据多组所述检测组合的叠加值，确定震源位置；

S1051：比较多组所述检测组合的叠加值，将多组所述叠加值中的最大值对应的检测组合的所在区域确定为震源位置。

所述分组成像条件为：

需要说明的是，地震波中S波从震源到台站的走时具体指S波从震源传播到台站所经历的时间，相应的地震波中P波从震源到台站的走时具体指P波从震源传播到台站所经历的时间。

所述波形序列的质量因子W与所述台站获取的波形序列的信噪比正相关，如果所述波形序列的信噪比(质量)较高，则所述波形序列的质量因子W接近于1，如果所述波形序列的信噪比(质量)较低，则所述波形序列的质量因子接近于0，即0<W<1。

所述台站的距离因子由台站和该台站所在检测组合的几何中心的距离决定，α是表示台站和该台站所在检测组合的几何中心的距离倒数的指数，用来控制距离因子这个变量对最终结果的影响大小，α越大则距离因子对最终结果的影响越大，α越小则距离因子对最终结果的影响越小。

下面将本申请实施例提供的震源的定位方法与基于现有技术中的成像条件的定位方法的对比测试进行说明。

一、首先进行模型测试：

1、模型测试：

我们采用相同的合成数据对本申请实施例提供的震源的定位方法与基于现有技术中的成像条件的定位方法分别进行抗噪性和分辨率测试。采用的台站分布如图4所示，图4中的坐标系为以水平面为XY平面建立的XYZ坐标系，X轴、Y轴和Z轴的单位均为千米(km)，Receiver arrays表示台站，Monitoring Wells表示观测井，YA、YB、YC、YD和YE表示所述观测井的编号，由5口观测井构成，每口观测井中有8个台站，但由于仪器问题，每口观测井中大概只有3-8个台站有效。因此，在合成测试中和实际数据一样只用了27个有效台站。从观测井数据中得到的一维速度模型如图5，图5中Vp表示P波速度，Vs表示S波速度。采用离散波数法(Bouchon,1981)生成合成波形数据。

2、抗噪性测试

真实震源位置设在[x＝0,Y＝0,Z＝1.5]km处，在合成波形上加入不同强度的噪声信号来测试抗噪性。与处理实际数据一样，我们采用10-35Hz相同的滤波频段。图6(a)-图6(c)分布是加上10％、50％、200％噪声的波形，图6(a)-图6(c)中横坐标表示时间(time)t，单位秒(s)，纵坐标表示不同的台站(Receiver)接收到的信号的波形，YA1、YA4、YB1、YB6、YC2、YC5、YC8、YD3和YE3分别为不同台站的编号，图中曲线分别为加上噪声后的波形和经过10-35Hz滤波后的波形。S波后面紧跟的一个强振幅信号是上边界反射波，实际波形则中不会出现，反射波正好可以用来检测当有一个相似的强振幅在真实相位附近该成像条件是否可以准确成像。用以下公式计算加入噪声后的波形信噪比SNR:

其中，U_signal是信号能量，U_noise是噪声能量大小。因此，加入10％、50％、200％噪声后波形的信噪比分别为10、2、0.5，值得注意的是，我们多观测井实际数据的信噪比大约在10左右。

对于本申请实施例提供的震源的定位方法，我们将同一口井中的台站分为一组，27个台站被分为5组。P波、S波的叠加窗长分别是30、40个采样点。由于每个检波器数据质量都相似，质量因子W为1，α为0.2。真实震源位置在(0,0,1.5)，发震时刻为1s。一共有四个搜索变量：震源位置(x,y,z)，发震时刻τ，搜索范围：x为-1.1～1km,y为-2.5～2.5km,z为0～2.5km，网格间隔为0.05km，发震时刻0.5～1.5s间隔为0.1s。图7(a)-图7(c)对比了本申请实施例提供的震源的定位方法(图7中以new method表示)和3个传统方法(线性叠加(图7中以linear表示)、波包叠加(图7中以envelope表示)、STA/LTA(图7中以STA/LTA表示)叠加)震源成像结果，能量值最大的点就是搜索得到的震源位置。随着噪声水平的增加，本申请实施例提供的震源的定位方法始终可以高分辨率的正确定位，而其他传统方法分辨率降低并且定位出现误差(参考表1)。

表1在不同信噪比下，四种成像条件定位结果对比

表1中，Source location表示震源位置，Linear stack表示采用线性叠加的成像条件的传统震源定位方法、Envelope stack表示采用波包叠加的成像条件的传统震源定位方法、STA/LTA stack表示采用长短时窗叠加的成像条件的传统震源定位方法、New method表示本申请实施例提供的震源的定位方法。

图8(a)-图8(d)显示了发震时刻反演结果：图8(a)为基于线性叠加成像条件的发震时刻反演结果，图8(b)为基于波包叠加成像条件的发震时刻反演结果，图8(c)为基于STA/LTA叠加成像条件的发震时刻反演结果，图8(d)为基于本申请实施例提供的分组成像条件的发震时刻反演结果，图8(a)-图8(d)中横坐标为时间Time，单位为s，纵坐标为值栈(Stackvalue)，SNR表示信噪比，从图8(a)-图8(d)中可以得出只有本申请实施例提供的震源的定位方法可以识别出准确的发震时刻。对比传统的成像条件，在时间和空间域，基于分组成像条件的结果抗噪性最好。

3、分辨率测试：

分别从空间分辨率和时间分辨率两方面对本申请实施例提供的震源的定位方法和基于传统成像条件的震源定位方法对比。空间分辨率测试中，设置两个发震时刻相同的相邻震源点(0,0,1.25)和(0,0,1.5)，定位结果见图9，图9中横坐标为Z轴位置，单位km，纵坐标为值栈(Stackvalue)。同样的，图9中本申请实施例提供的震源的定位方法以Newmethod表示，3个传统方法(线性叠加(图9中以linear stack表示)、波包叠加(图9中以Envelope stack表示)、STA/LTA(图9中以STA/LTA stack表示)叠加)。

从图9中可以看出，只有本申请实施例提供的震源的定位方法能够检测出1.25km和1.5km处的事件，传统成像条件不能准确定位。

图10为时间分辨率测试结果，真实发震时刻在1s和1.2s。波包和STA/LTA叠加只能检测出一个事件，线性叠加和基于分组成像条件可以准确检测出两个事件，但是本申请实施例提供的震源的定位方法更加明显。

从合成测试结果可以看出，基于分组的成像条件的高分辨率全波形地震偏移定位算法可以在信噪比较低的情况下稳定正确的定位。

二、多井实际数据地震定位

我们将本申请实施例提供的震源的定位方法和传统成像条件都应用于阿曼一油田的井中检测数据进行对比。5口监测井设置在地震活动性较强的区域，每口井有3-8个三分量检波器，大约在750-1250米深度处(图4)。在记录到的波形数据中，前人检测到了超过15800个地震事件，并对大约5400个事件进行了定位。图11显示了事件的位置(圆点)和台站的分布(三角)，图11中横坐标为东部(East)的位置，单位为km，纵坐标为北部(North)位置，单位为km，Group A和Group B分别表示两组检测组合，图11中实线是两条主要的平行断层，大部分地震事件沿着这两条断层分布。我们从两条断层上各选取了15个事件进行定位，分为A组合B组，总的来说A组数据质量优于B组。用本申请实施例提供的震源的定位方法定位，27个台站分为5组，叠加窗口为200个采样点，由于每组的台站距离较近，距离控制因子α为0.2。

我们先选取了一个事件做不确定性测试，它的垂直分量波形如图12，图12中横坐标为时间t(time)，单位为s，纵坐标为观测井中的台站Receiver。除了观测井YE#，大多数台站数据质量较好，因此，井YE#的质量因子W为0.5，其余均为1。如图13所示，图13中从上到下分别是基于线性叠加成像条件的原始波形，基于波包叠加成像条件的波形包络，基于STA/LTA叠加成像条件的STA/LTA特征函数，基于本申请实施例提供的分组成像条件的特征函数包络，3个传统成像条件(在图13中分别以vertical component、envelope、STA/LTA)和本申请实施例提供的震源的定位方法(图13中以characteristic envelope表示)这四种方法计算出的波形特征对到时灵敏度不同。原始波形和波包相比于其他两种方法到时不是很明显。STA/LTA和特征函数波包都有明显的到时信号，但是STA/LTA有其他信号产生的假象。震源成像结果见图14，相比于其他传统成像条件，新成像条件能够最大程度消除假象，不确定性最低。

我们将本申请实施例提供的震源的定位方法定出的30个事件结果和Zhang(2009)用基于走时的双差定位方法定出的结果进行对比。对比结果见图15(a)-图15(d)，STA/LTA和基于分组成像条件都与走时反演有很好的一致性，但是线性和波包叠加成像条件的结果有一定偏差。值得注意的是，我们是将基于人工拾取到时的走时成像方法结果最为参考，它本身并不是真实事件位置可能存在一定误差。总体来说，我们提出的本申请实施例提供的震源的定位方法的定位结果与它一致说明我们的方法是可接受的，有一定实际意义。

下面对本申请实施例提供的震源的定位系统进行描述，下文描述的震源的定位系统可与上文描述的震源的定位方法相互对应参照。

相应的，本申请实施例提供了一种震源的定位系统，包括：

所述分组成像条件为：

所述第一预设公式为：

综上所述，本申请实施例提供了一种震源的定位方法及定位系统，其中，所述震源的定位方法在震源定位过程中通过台站检测地震波的波形序列，无需人工手动拾取地震波走时，降低了震源的定位方法的人工成本，从而降低了震源的定位方法的整体成本；此外，所述震源的定位方法采用了全新的分组成像条件，该分组成像条件将传统成像条件中的相加变为分组相乘，并加入距离、质量因子控制提高每组成像质量，使得所述震源的定位方法有效地在空间上和时间上对震源进行成像，提高了震源的定位方法的分辨率，有助于计算储层体积，研究裂缝分布和评价压裂效果。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。