CN104570076A - 一种基于二分法的地震波初至自动拾取方法 - Google Patents

Abstract

本发明提供了一种基于二分法的地震波初至自动拾取方法，属于石油地球物理勘探领域。本方法包括：获取参考初至：分别计算地震波中直达波、反射波和折射波在水平层状介质中的时距曲线，根据这三个时距曲线得到三个理论初至，然后取这三个理论初至中最小的一个作为参考初至；利用所述参考初至，采用二分法获取初至点；确定波峰，波谷和起跳点：根据拾取到的初至点确定波峰，波谷和起跳点。

Description

一种基于二分法的地震波初至自动拾取方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种基于二分法的地震波初至自动拾取方法。

背景技术

地震波初至拾取是地震资料处理中的一个基础又重要的环节，然而初至拾取工作单一且耗费时间很长，极大地影响震资料处理的进度和效果。近些年来，地球物理研究人员针对此初至拾取已经提出上百种自动拾取算法，归纳为能量比算法、分形维算法、相关法、神经网络法、边缘检测法等几大类。这些算法各有特点，但都具有“速率快的精度较低，精度高的速率较低”的不足之处。目前很多资料处理软件中都具有初至自动拾取这一模块，为了提高拾取准确性，都需要设置一些参数，例如Omega中关于pickworks就设置了很多参数项，初至状态，最大振幅限制，速度等；在Opencps中，参数设置的更多，速度，阈值，选择时窗，基础时窗等。

现阶段，随着勘探区域的复杂以及勘探技术的提高，野外采集的数据量激增，现有的初至自动拾取已经跟不上实际生产需求的需要，需要进一步的提高初至自动拾取的效率。相对于其他类的算法，能量比算法能够清楚的描述初至的主要特点--能量瞬时增大。实现上相对简单，借鉴其他初至自动拾取软件的特点，通过设计参数可以显著地提高拾取的准确性。最主要的特点是速度较快，这适用于数据激增的地震数据。

初至广泛的应用于折射波静校正，VSP解释，地震层析成像中，准确的拾取地震波的初至可以准确的求取静校正量，计算出可靠的平均速度和层速度，关系着下覆地层或结构几何参数和速度参数的正确反演。因此，初至自动拾取算法的好坏直接关系到后续处理的效率和结果。

同时由于数据量的激增，初至拾取任务劳动强度非常大，这对处理人员是非常不利的，而好的自动拾取算法可以帮助处理人员解决这一问题，大大减少劳动消耗，就有更多的时间和精力用于其他处理工作。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于二分法的地震波初至自动拾取方法，准确的拾取初至，减轻处理人员的劳动强度，给予其他处理过程更多的时间和精力，最后得到比较完美的处理结果。

本发明是通过以下技术方案实现的：

一种基于二分法的地震波初至自动拾取方法，包括：

获取参考初至：分别计算地震波中直达波、反射波和折射波在水平层状介质中的时距曲线，根据这三个时距曲线得到三个理论初至，然后取这三个理论初至中最小的一个作为参考初至；

利用所述参考初至，采用二分法获取初至点；

确定波峰，波谷和起跳点：根据拾取到的初至点确定波峰，波谷和起跳点。

所述获取参考初至具体如下：

第i道的参考初至为：

t_i=min{t_i1，t_i2，t_i3}，

其中直达波的时距曲线为：

反射波的时距曲线为：

折射波的时距曲线为： $t_{i3}=\frac{\sqrt{{(Y_s+Y_i)}^2+K_i^2}}{V_1}+\frac{|X_i-K_i|}{V_2},$

其中V₁为估计的地表速度，V₂表示估计的地表的下一层的速度，Y_s表示炮点高程，Y_i表示第i道检波点高程，X_i表示第i道的偏移距。

所述利用所述参考初至，采用二分法获取初至点包括：

(1)通过二分法得到不同的拾取范围：

对于第i道地震数据，第一次拾取范围的起点为起始位置到参考初至的中点，即其中t_i表示第i道的参考初至，Δt为采样间隔，[]为取整运算符；

第一次拾取范围的终点为参考初至到终止位置的中点，即 其中N为单道地震数据的采样点数；

第二次拾取范围的起点为第一次拾取范围的起点到参考初至的中点，即 $N_{\mathrm{ib}2}=[\frac{1}{2}+{\left(\frac{1}{2}\right)}^2)\frac{t_i}{\mathrm{\Δt}}];$

第二次拾取范围的终点为参考初至到第一次拾取范围的终点的中点，即 $N_{\mathrm{ie}2}=[\frac{t_i}{\mathrm{\Δt}}+{\left(\frac{1}{2}\right)}^2*\frac{\mathrm{\Δt}*N-t_i}{\mathrm{\Δt}}];$

依此类推，第j次拾取范围的起点为：第j次拾取范围的终点为 $N_{\mathrm{ie}2}=[\frac{t_i}{\mathrm{\Δt}}+{\left(\frac{1}{2}\right)}^j*\frac{\mathrm{\Δt}*N-t_i}{\mathrm{\Δt}}].$

(2)在各个拾取范围内分别获取初至点：

对于第i道，参考初至时刻为t_i，在第一次迭代过程中，前一个拾取范围为[N_ib1,N_ie1]，后一个拾取范围为[N_ib2,N_ie2]，若在这两个拾取范围内寻找到的初至点是同一个，则说明在第一迭代过程中找到了该道的初至，如果不相同，则进入第二次迭代过程；

在第二迭代过程中，前一个拾取范围为[N_ib2,N_ie2]，后一个拾取范围为[N_ib3,N_ie3]，若在这两个范围内寻找到的初至点是同一个，则说明在第二迭代过程中找到了该道的初至，如果不相同，则进入第三次迭代过程；

依此类推，直到前后两次拾取范围内的初至点是同一个为止，即找到该道的初至点；

(3)对于每一炮集，如果在第j个迭代过程中具有稳健性初至的道集占整炮道集的比例大于80％，则停止迭代，对于具有稳健性初至的道集，选择其各自的第j个迭代过程中拾取的初至点为其各自的初至点，对于不具有稳健性初至的道集，选择其各自的前一个范围内拾取的初至点为其各自的初至点。

所述具有稳健性初至是指在前后两次拾取范围内拾取的初至点是同一个。

在每个拾取范围内获取初至点是通过在该拾取范围内采用滑动时窗能量比算法得到初至的。

所述根据拾取到的初至点确定波峰，波谷和起跳点是采用可变时移量相位推算方法实现的。

与现有技术相比，本发明的有益效果是：为了考察和检验本发明初至自动拾取算法的有效性和精度，在相同共炮点地震数据条件下，采用当前主流商用软件(OMEGA)的初至自动拾取功能进行测试盒对比分析。图7是商用软件(OMEGA)初至自动拾取的成果展示，图8是本发明初至自动拾取的，效果展示，可以看出，本发明初至自动拾取的效果优于当前同类的主流商用软件。

附图说明

图1为直达波理论初至示意图；

图2为反射波理论初至示意图；

图3为折射波理论初至示意图；

图4为参考初至示意图，其中地表速度估计为1900m/s，地表的下一层的速度估计为2800m/s；

图5为高信噪比资料采用本算法自动拾取效果示意图；

图6为一般信噪比资料采用本算法自动拾取效果示意图；

图7为低信噪比资料采用Omega自动拾取效果示意图；

图8为低信噪比资料采用本算法自动拾取效果示意图。

图9为本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

在进行初至自动拾取时主要分为以下三步，首先确定初至大概位置，然后基于二分法得到初至，最后利用相位时移技术精确地得到波峰、波谷和起跳点。这样可以快速准确的进行初至拾取。

如图9所示，本方法包括：

(1)粗略确定初至的位置

通过估计表层速度计算参考初至，得到初至所在的大概位置。

首先给出参考初至的概念，对于水平层状介质，通过估计地下速度，计算地震波中直达波、反射波、折射波的时距曲线，得到理论上检波点首次接收到的波的时刻定义为参考初至，具体如下：

分别计算这三种波在水平层状介质中的时距曲线，其中表层的深度用炮点高程与检波点高程中数值大的一个近似替代，然后取这三个时间中的最小的一个作为所需要的参考初至。也就是说，分别根据直达波的时距曲线、反射波的时距曲线、折射波的时距曲线得到反射波、折射波的理论初至，再取这三个理论初至中的最小的一个作为所需要的参考初至。那么第i道的参考初至为：t_i=min{t_i1，t_i2，t_i3}，其中直达波：反射波：折射波： $t_{i3}=\frac{\sqrt{{(Y_s+Y_i)}^2+K_i^2}}{V_1}+\frac{|X_i-K_i|}{V_2},$ 其中 $K_i=a_0*Y_s/(1+\frac{Y_i}{Y_s}),a_0=\frac{V_1}{\sqrt{V_2^2-V_1^2}},$ V₁为估计的地表速度，V₂表示估计的地表的下一层的速度，Y_s表示炮点高程，Y_i表示第i道检波点高程，X_i表示第i道的偏移距。

理论初至表示理论上检波点首次接收到的波的时刻，这里的波可以表示任何一种波，因此下面的图1图2图3图4中可以用到理论初至这个概念。

其次，地震初至波中，大部分都是直达波、反射波、折射波，利用上述三种波的时距曲线既可以涵盖大多数的初至波，又可以节约计算时间，提高效率。最后，考虑到用户的操作性，可以通过简单地操作获得速度，而不需要外部输入，具体如下：

在炮集数据剖面上，在某一道(除炮点位置的任意一道均可)初至的位置点左击一下，获得此处的初至时间，通过直达波的时距曲线可以计算出表层速度，对表层速度增加一个数值得到表层的下一层速度，本算法中数值取为900。

在介绍图像之前，首先给出理论初至的概念：通过估计地下速度，根据波的时距曲线得到理论上检波点首次接收到的波的信号，称为理论初至。图1表示根据直达波得到的理论初至；图2表示根据反射波得到的理论初至；图3表示根据折射波得到的理论初至；图4表示参考初至，即取这三种波得到的理论初至的最小值。可以看出选择直达波、反射波、折射波作为初至波，得到的参考初至与初至最为接近。

优化地，当用户通过单击初至位置得到参考初至后，用户根据在界面上显示的参考初至是否合理(即初至是否是最贴近参考初至的)来选择是否重新调整参考初至。

(2)通过二分法获取初至

能量比算法的一个不足之处在于当拾取范围不同时得到的结果有很大的区别，即得到的初至没有稳健性。鉴于此，本发明提出，为了提高能量比算法的稳健性，通过二分范围得到不同的拾取范围，分别得到初至点，若前后两次范围内初至一样，则说明在这一迭代范围内此道的初至(每一道只有一个初至)已经找到，否则，继续缩小范围进行查找，直到满足条件为止。此段说明了对于炮集中任意的地震道i，如何确定迭代过程。

首先，给出二分范围的具体操作过程：

对于第i道地震数据，取第一次拾取范围的起点为其中t_i表示第i道的参考初至，Δt为采样间隔，[]为取整运算符，即为起始位置到参考初至的中点；第一次拾取范围的终点其中N为单道地震数据的采样点数，即为参考初至到终止位置的一半作为第一次范围的终点；取第二次范围的起点为即为第一次范围的起点到参考初至的中点；取第二次范围的终点为即为参考初至到第一次范围的终点的中点；依次类推，可以得到第j次范围的起点为 $N_{\mathrm{ibj}}=[\frac{1}{2}+{\left(\frac{1}{2}\right)}^2+...+$ ${\left(\frac{1}{2}\right)}^j\left)\frac{t_i}{\mathrm{\Δt}}\right],$ 终点为 $N_{\mathrm{ie}2}=[\frac{t_i}{\mathrm{\Δt}}+{\left(\frac{1}{2}\right)}^j*\frac{\mathrm{\Δt}*N-t_i}{\mathrm{\Δt}}].$

其次，对于每一炮集，如果在第j个迭代过程中具有稳健性初至的道集占整炮道集的比例大于80％，则选择第j个迭代过程中拾取的初至点，此处不具有稳健性初至的道集选择的是前一个范围内拾取的初至点；另一个准则是迭代过程不能超过13个，这是根据2¹³=8129>8001，单道中采样点个数一般不会超过8001个，在实际资料中，一般选择迭代过程为2个，太多会影响拾取速度。

具体来说，在不同的范围内，采用滑动时窗能量比算法得到初至，并比较在前后两次范围中拾取的具有稳健性初至的道集(即在前后两次范围内拾取的初至点为同一个初至点，)占整个炮集的比例，如果大于80％，则不再继续寻找下一范围(寻找初至时是每道都要找到初至，寻找稳健性初至时是以一炮为单位，查看稳健性初至所在的道集所占的比例)，否则，继续寻找下一范围，并在前后两个范围拾取初至并比较是否满足给出的比例条件。最后，为了避免寻找范围次数过多造成算法的速率降低，根据比较次数(用户通过交互操作给出，一般情况下选择2次)(保证初至点位于拾取范围内)作为结束条件。此段说明了在整个炮集中如何根据迭代过程中的稳健性初至点确定使用第几个迭代过程。

滑动时窗能量比算法是能量比算法的一种，能量比算法以地震记录初至前后时窗内的能量差异特征为基础，参考文献为：《高精度初至自动拾取综合方法研究》(《物探与化探》，34(5)，595-599，2010)，简述如下：

首先，按步长将地震道分为S份；

其次，以公式(1)计算前后两个时窗内的能量比值：

$R=\frac{{\left[{\mathrm{\Σ}}_{t=N0}^{N2}{x}^2\left(t\right)\right]}^{1/2}+\mathrm{\αA}}{{\left[{\mathrm{\Σ}}_{t=N1}^{N0}{x}^2\left(t\right)\right]}^{1/2}+\mathrm{\αA}}---\left(1\right)$

其中，R为能量比，A为一个地震道的相对能量，x(t)为地震记录振幅值，N1为时窗起点，N0为时窗中点，N2为时窗终点，N为所用地震道的样点数，α为稳定系数(这个值是用户通过交互界面从外界输入的，给出的默认值为0.2)；

最后，找到能量比值最大的时窗，在此时窗内寻找具有最大振幅的采样点，即为初至。

通过不断的变换范围，可以更加准确的描述初至特点：第一次能量瞬时增大。原因如下：一方面拾取范围一直包含初至点，另一方面拾取范围的起点和终点距离初至点越来越接近，这样就可以排除其他位置的能量瞬时增大，达到准确描述初至的目的。

(3)确定波峰，波谷，起跳点

由于处理人员的不同需求，有必要确定初至的状态，本发明采用的是可变时移量相位推算技术，根据拾取到的初至确定波峰，波谷，起跳点，这里借用的文献《初至智能拾取技术》(《石油物探》，46(5)，521-530，2007)所给出的方法，简述如下：

对于第i道地震数据，首先，根据上一步拾取的初至，找到这个初至的上下共M个采样点，其中N表示单道地震数据的采样点数；

其次，根据公式(2)求这些采样点的希尔伯特变换：

$h\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}---\left(2\right)$

其中，x(t)为t时刻地震道的振幅值，h(t)为x(t)对应的希尔伯特变换；

然后寻找满足公式(3)的点，即为波峰N_peak：

$\cos \mathrm{\θ}\left(t\right)=\frac{x\left(t\right)}{a\left(t\right)}=1---\left(3\right)$

其中 $a\left(t\right)=\sqrt{x^2\left(t\right)+h^2\left(t\right)},$

a(t)为瞬时包络，θ(t)为瞬时相位，cosθ(t)为相位函数；

最后依据公式(4)和公式(5)分别求取波谷N_trough和起跳点N_offset：

N_trough＝N_peak-τ₁         (4)

N_offset=N_peak-τ₂        (5)

其中 ${\mathrm{\τ}}_1=\frac{2}{m}\left({\mathrm{\Σ}}_{j=i-\frac{m}{2}}^{i+\frac{m}{2}}\mathrm{\Δ}{\mathrm{\τ}}_j\right),{\mathrm{\τ}}_2=\frac{3}{m}\left({\mathrm{\Σ}}_{j=i-\frac{m}{2}}^{i+\frac{m}{2}}\mathrm{\Δ}{\mathrm{\τ}}_j\right),$ m为当前道i的前后相邻的m道，Δτ_j为第j道初至波波峰与正零点之差，满足的点为正零点。

本发明是一种提高初至自动拾取精度和速度的方法，下面以几个不同工区的资料来说明本发明的效果。首先给出信噪比高的资料，每炮含有280道，每道4001个采样点，采样间隔为2ms，从图5中可以看出，当选择Peak时，利用本发明可以完全正确的拾取到初至波峰位置。其次给出具有一般信噪比的资料，每炮含有260道，每道1001个采样点，采样间隔为2ms，图6中结果显示在空道处拾取的不准确，其他位置都可以准确的拾取到波峰。最后给出信噪比较低的资料，每炮含有480道，每道8001个采样点，采样间隔为0.5ms，图7为Omega自动拾取的效果，从图中可以看出，自331道开始至480道，拾取的初至不准确，图8为本发明的自动拾取效果，可以看出本发明在远偏移距处可以拾取到相对正确的结果。

本发明操作过程分为以下几步：

第一步：通过交互界面对话框，填写所需要的两个速度(地表速度和地表的下一层的速度)或者在剖面上初至位置处单击(单击初至位置后，程序中首先获取了初至的时刻，然后根据时距曲线得到地表速度，并根据程序中地表的下一层速度比上一层速度大900的数值得到地表的下一层速度，然后再根据直达波、反射波和折射波的时距曲线计算理论初至，得到参考初至；

第二步：依据处理需要，选择波峰、波谷、起跳点中的一项；

第三步：填写时窗长度、步长、比例系数以及最大迭代次数(其中时窗长度、步长、比例系数是根据滑动时窗能量比算法，需从外界交互得到，最大迭代次数与二分法有关，相当于对炮集地震数据求取了最大迭代次数次通过二分法得到的范围，这个最大迭代次数的作用是防止程序进入死循环)，得到初至(通过点击界面上的Apply键即可在界面上显示拾取到的初至)；

第四步：调整第三步中的参数(如果用户在界面上看到拾取的初至结果不理想，可以通过改变上述几个参数的输入值来修改拾取的结果)，得到理想的结果(对于一般的资料，自动拾取结果达到90％以上即可，但最终仍由用户决定)。

本发明提出了一种地震道初至自动拾取算法，通过引入参考初至，得到初至所在的范围，然后利用二分法准确的获得初至，并利用相位时移技术进行微调，得到波峰，波谷，起跳点。最后对不同信噪比的地震资料进行拾取，结果表明该算法具有快速，准确的特点。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (6)

1.一种基于二分法的地震波初至自动拾取方法，其特征在于：所述方法包括：

获取参考初至：分别计算地震波中直达波、反射波和折射波在水平层状介质中的时距曲线，根据这三个时距曲线得到三个理论初至，然后取这三个理论初至中最小的一个作为参考初至；

利用所述参考初至，采用二分法获取初至点；

确定波峰，波谷和起跳点：根据拾取到的初至点确定波峰，波谷和起跳点。

2.根据权利要求1所述的基于二分法的地震波初至自动拾取方法，其特征在于：所述获取参考初至具体如下：

第i道的参考初至为：

t_i=min{t_i1，t_i2，t_i3}，

其中直达波的时距曲线为：

反射波的时距曲线为：

折射波的时距曲线为： $t_{i3}=\frac{\sqrt{{(Y_s+Y_i)}^2+K_i^2}}{V_1}+\frac{|X_i-K_i|}{V_2},$

其中V₁为估计的地表速度，V₂表示估计的地表的下一层的速度，Y_s表示炮点高程，Y_i表示第i道检波点高程，X_i表示第i道的偏移距。

3.根据权利要求2所述的基于二分法的地震波初至自动拾取方法，其特征在于：所述利用所述参考初至，采用二分法获取初至点包括：

(1)通过二分法得到不同的拾取范围：

对于第i道地震数据，第一次拾取范围的起点为起始位置到参考初至的中点，即其中t_i表示第i道的参考初至，Δt为采样间隔，[]为取整运算符；

第一次拾取范围的终点为参考初至到终止位置的中点，即 其中N为单道地震数据的采样点数；

第二次拾取范围的起点为第一次拾取范围的起点到参考初至的中点，即 $N_{\mathrm{ib}2}=[\frac{1}{2}+{\left(\frac{1}{2}\right)}^2)\frac{t_i}{\mathrm{\Δt}}];$

第二次拾取范围的终点为参考初至到第一次拾取范围的终点的中点，即 $N_{\mathrm{ie}2}=[\frac{t_i}{\mathrm{\Δt}}+{\left(\frac{1}{2}\right)}^2*\frac{\mathrm{\Δt}*N-t_i}{\mathrm{\Δt}}];$

依此类推，第j次拾取范围的起点为：第j次拾取范围的终点为 $N_{\mathrm{ie}2}=[\frac{t_i}{\mathrm{\Δt}}+{\left(\frac{1}{2}\right)}^j*\frac{\mathrm{\Δt}*N-t_i}{\mathrm{\Δt}}].$

(2)在各个拾取范围内分别获取初至点：

对于第i道，参考初至时刻为t_i，在第一次迭代过程中，前一个拾取范围为[N_ib1,N_ie1]，后一个拾取范围为[N_ib2,N_ie2]，若在这两个拾取范围内寻找到的初至点是同一个，则说明在第一迭代过程中找到了该道的初至，如果不相同，则进入第二次迭代过程；

在第二迭代过程中，前一个拾取范围为[N_ib2,N_ie2]，后一个拾取范围为[N_ib3,N_ie3]，若在这两个范围内寻找到的初至点是同一个，则说明在第二迭代过程中找到了该道的初至，如果不相同，则进入第三次迭代过程；

依此类推，直到前后两次拾取范围内的初至点是同一个为止，即找到该道的初至点；

(3)对于每一炮集，如果在第j个迭代过程中具有稳健性初至的道集占整炮道集的比例大于80％，则停止迭代，对于具有稳健性初至的道集，选择其各自的第j个迭代过程中拾取的初至点为其各自的初至点，对于不具有稳健性初至的道集，选择其各自的前一个范围内拾取的初至点为其各自的初至点。

4.根据权利要求3所述的基于二分法的地震波初至自动拾取方法，其特征在于：所述具有稳健性初至是指在前后两次拾取范围内拾取的初至点是同一个。

5.根据权利要求4所述的基于二分法的地震波初至自动拾取方法，其特征在于：在每个拾取范围内获取初至点是通过在该拾取范围内采用滑动时窗能量比算法得到初至的。

6.根据权利要求1所述的基于二分法的地震波初至自动拾取方法，其特征在于：所述根据拾取到的初至点确定波峰，波谷和起跳点是采用可变时移量相位推算方法实现的。