CN104570113A - 一种自适应去除地震强反射的方法 - Google Patents

Abstract

本发明提供了一种自适应去除地震强反射的方法，属于油气田勘探开发中地震资料预处理领域。本方法包括：(1)输入三维地震数据体，利用希尔伯特变换计算整个三维地震数据体中每道地震信号的瞬时相位余弦cosθ(t)，得到瞬时相位余弦数据体；(2)针对瞬时相位余弦数据体，采用互相关法自动拾取指定点位置所在的层位，记录每道强反射层位的反射时间t_i；所述的指定点为用户预先设定的强反射在某一道上对应的时间；(3)对输入的一道地震数据，根据该道自动拾取的层位时间，根据最大匹配投影原理匹配出最优小波；然后用该道的地震数据减去匹配的最优小波，得到该道的去除指定强反射后的地震数据。

Description

一种自适应去除地震强反射的方法

技术领域

本发明属于油气田勘探开发中地震资料预处理领域，具体涉及一种自适应去除地震强反射的方法。

背景技术

在地震勘探中，地震资料的主频通常都不高、分辨率较低，而较多的储层都为薄储层，同时受目的层段无效强反射(如煤层、海底强反射等)的影响，目的层的有效反射被淹没在强反射之中，很难从常规的地震资料中识别出目标砂体。

通常层位解释都是针对原始地震剖面的，而原始地震剖面中振幅能量的强弱会影响追踪层位的走向。解释过程中或者参考瞬时相位剖面，而地震剖面的波峰在瞬时相位剖面中对应在0值附近，波谷在瞬时相位剖面中对应-π和π的突变位置，与原始地震数据中同相轴对应关系不直观。而瞬时相位余弦独立于地震振幅强度，它的主波峰、主波谷与地震信号的波峰、波谷对应较好，能够很好的反映地层连续性，可以有效分辨地震波的几何形态，能够用来识别层序的边界以及地震相的外部几何形态。因此，瞬时相位余弦是一个识别地震相几何形态和内部结构较好的几何属性参数。

匹配追踪算法是一种最典型的参数化时频分析方法，它是将信号表示为一系列与信号局部结构特征最佳匹配的时频原子的线性组合，其算法核心为最大匹配投影原理。同时Morlet小波能很好的近似实际地震子波，可以通过对该小波进行伸缩、频移和调制来匹配实际地震信号，得到较好的分解效果，能够更合理的表征信号的局部结构特征。

目前还没有专门去除地震强反射的方法。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种自适应去除地震强反射的方法，采用一种自适应方法去除指定位置的干扰同相轴，突出目标区域的有效反射能量，为精细储层预测服务。本发明充分考虑到层位自动拾取过程中反射能量对拾取结果的影响，采用瞬时相位余弦剖面代替常规的地震反射剖面，能够更准确的拾取指定的反射时间。考虑到Morlet小波能很好的近似实际地震子波，在沿层拾取最优小波过程中采用Morlet小波作为时频原子，保证拾取的地震反射能更合理的表征信号的局部结构特征。

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

一种自适应去除地震强反射的方法，包括：

(1)输入三维地震数据体，利用希尔伯特变换计算整个三维地震数据体中每道地震信号的瞬时相位余弦cosθ(t)，得到瞬时相位余弦数据体；

(2)针对瞬时相位余弦数据体，采用互相关法自动拾取指定点位置所在的层位，记录每道强反射层位的反射时间t_i；所述的指定点为用户预先设定的强反射在某一道上对应的时间；

(3)对输入的一道地震数据，根据该道自动拾取的层位时间，根据最大匹配投影原理匹配出最优小波；然后用该道的地震数据减去匹配的最优小波，得到该道的去除指定强反射后的地震数据；

(4)对整个三维地震数据中的每一道重复步骤(3)，得到每一道的去除指定强反射后的地震数据，所有道去除指定强反射后的地震数据构成去除指定强反射后的地震数据体。

所述步骤(1)具体包含以下步骤：

(11)对每一道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)；

(12)计算地震信号的相位余弦cosθ(t)，如式(1)所示：

$\cos \mathrm{\θ}\left(t\right)=\frac{x_r\left(t\right)}{\sqrt{x_r^2\left(t\right)+x_i^2\left(t\right)}}---\left(1\right)$

所有道地震信号的相位余弦就构成所述瞬时相位余弦数据体。

所述步骤(2)中采用互相关法自动拾取指定点位置所在的层位是这样实现的：确定互相关时窗长度和相对指定位置上下漂移距离，然后使用这两个参数自动拾取得到指定点位置所在的层位；

所述互相关时窗长度为至少包含一个完整的波形。

所述步骤(3)中的根据最大匹配投影原理匹配出最优小波具体包含以下步骤：

(31)对一道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)；计算复信号的瞬时频率和瞬时相位；其中u_n为对应道的时间t_i，f_n为对应时间的瞬时频率，φ_n为对应时间的瞬时相位；

(32)在一组固定的u_n、f_n和φ_n情况下，通过计算最优化公式 $\arg \underset{g_{r_n}\&Element;D}{\max }\frac{|<<R^{\left(n\right)}s,g_{r_n}>>|}{\left|\right|g_{r_n}\left|\right|}$ 得到尺度σ_n；

其中D＝{g_r(t)}_r∈Γ为时频原子字典，是函数R⁽ⁿ⁾s和的内积， $\left|\right|g_{r_n}\left|\right|=\sqrt{<<g_{r_n},g_{r_n}>>};$

(33)利用所述最优化公式在局域内寻找r_n＝{σ_n，f_n，φ_n}的最优值，即公式达到最大时，对应的一组参数便是最优的r_n；搜索范围为[r_n-Δr，r_n+Δr]，其中Δr＝(Δσ，Δf，Δφ)，Δσ为尺度偏移量、Δf为频率偏移量、Δφ为相位偏移量；

(34)通过下式预测到最优的时频原子的振幅为a_n，即 $a_n=\frac{|<<R^{\left(n\right)}s,g_{r_n}>>|}{<<g_{r_n},g_{r_n}>>};$

(35)基于步骤(31)至(34)得到的决定Morlet小波形态的五个参数，根据Morlet小波的表达式(2)得到对应道的最优小波：

$g_{r_n}\left(t\right)=a_n\exp [-4\ln 2\frac{f_n^2{(t-u_n)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left(i(2\mathrm{\&pi;}{f}_n(t-u_n)+{\mathrm{\&phi;}}_n)\right)---\left(2\right).$

与现有技术相比，本发明的有益效果是：

1)本发明方法提供了一种自适应去除地震强反射的方法，根据指定强反射的位置，自适应的去除指定的强反射。

2)本发明方法中拾取层位采用瞬时相位余弦剖面代替常规的地震剖面，避免了地震反射能量对层位追踪的影响，能够更准确的拾取指定强反射的时间。

3)本发明方法采用Morlet小波作为时频原子，在局域内对时频原子进行伸缩、频移和相移得到一个原子库，根据最大匹配投影原理寻找强反射处的最优小波，能够更合理的表征信号的局部结构特征。

附图说明

图1本发明自适应去除地震强反射的步骤框图。

图2本发明实施例中的一道地震信号及其对应的瞬时相位余弦函数曲线。

图3本发明实施例中的一条原始地震剖面，图中1500ms左右显示的黑线为自动拾取的强反射层位；

图4本发明实施例中的一条原始地震剖面的瞬时相位余弦剖面，图中1500ms左右显示的一条线为自动拾取的指定位置的强反射同相轴；

图5本发明实施例中的沿强反射同相轴拾取的最优子波，图中1500ms左右显示的黑线为自动拾取的强反射层位；

图6本发明实施例中的去除强反射后的地震剖面，图中1500ms左右显示的黑线为自动拾取的强反射层位；

具体实施方式

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

一种自适应去除地震强反射的方法，所述方法首先计算整个三维地震数据体瞬时相位余弦；针对瞬时相位余弦数据体，采用互相关法自动拾取指定强反射的反射时间；沿着拾取的反射时间，根据最大匹配投影原理对每道寻找最优小波，每一道的地震数据减去对应的最优小波，便得到去除指定强反射的地震数据体。

如图1所示，本发明自适应去除地震强反射的方法包括以下步骤：

(1)输入三维地震数据体，利用希尔伯特变换计算整个三维地震数据体中每道地震信号的瞬时相位余弦cosθ(t)，得到瞬时相位余弦数据体；

(2)针对瞬时相位余弦数据体，采用互相关法自动拾取指定点位置所在的层位(只需要用户设定一个时间点，本发明就能根据这个时间点自动的拾取出相应的层位出来)，记录每道强反射的反射时间t_i(本发明是用于去除强反射的，所以指定点应该是强反射在某一道上对应的时间，有了这个时间点之后，根据本发明自适应的特点，就能搜索出整个强反射层位的t_i，在指定的道上，指定时间就是t_i，这里相当于是给了个初始值，这个时间是在进行匹配追踪时的u_n，在(31)中有说明)；所述的指定点位置为用户预先设定的想去除的强反射所在的时间；

(3)对输入的一道地震数据，根据该道自动拾取的层位时间，根据最大匹配投影原理匹配出最优小波；然后用该道的地震数据减去匹配的最优小波，得到该道的去除指定强反射后的地震数据(两个信号长度一致，一一对应的减，就是两个等长的信号相减，或者说是两个等长的向量相减。)；

(4)对整个三维地震数据中的每一道重复步骤(3)，便得到去除指定强反射后的地震数据体。

所述步骤(1)具体包含以下步骤：

(11)对每一道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)；

希尔伯特变换是公知的常规算法，本领域的技术人员都可以准确无误的从x_r(t)得到唯一的x_i(t)。

(12)计算地震信号的相位余弦cosθ(t)，如式(1)所示：

$\cos \mathrm{\&theta;}\left(t\right)=\frac{x_r\left(t\right)}{\sqrt{x_r^2\left(t\right)+x_i^2\left(t\right)}}---\left(1\right)$

所有道地震信号的相位余弦就构成所述瞬时相位余弦数据体。

所述步骤(2)中采用互相关法自动拾取指定点位置所在的层位，涉及到几个重要的参数：互相关时窗长度、相对指定位置上下漂移距离(漂移距离就是计算互相关时的上下移动的时间大小，这个需要根据实际的情况设置)；互相关时窗长度应根据需要去除的波形的长度确定，波形长度与地震资料主频或分辨率有关，主频越高时窗越小，主频越低时窗越大；但是互相关时窗长度至少包含一个完整的波形(通常就是比一个完整的波形稍微长一点儿，一般不超过两个波形，太长的话就没有意义了)，即一峰二谷或一谷两峰，视去除同相轴为波峰、波谷而定；

所述步骤(2)中采用互相关法自动拾取指定点位置所在的层位，具体实施如下：如指定的强反射的道所在的线号为m(线范围为1-M，1＜m＜M)，道号为n(道范围为1-N，1＜n＜N)，指定强反射的时间为t₀.首先在该条线m上利用互相关向起始道方向和终止道方向分别采用逐道递推的方式自动追踪层位。如互相关时窗长度为wa，向上漂移距离为a，向下漂移距离为b；计算n-1道的层位时间时，在第n道上取以第n道层位时间为t₀中心，长度为wa的数据作为基准数据；在第n-1道上以t₀为中心取长度为wa的数据与基准数据计算互相关，得到一个互相关系数；以长度wa向上逐点滑动，直到漂移到a长度，向下逐点滑动，直到漂移到b长度，分别计算与基准数据的互相关值，便能得到a+b+1个互相关值，互相关值最大时对应的时窗的中心点便是第n-1道强反射层位的时间。按此方式向两端递推，便能得到指定线上的强反射的时间。然后以这条线上的强反射的层位时间为基准，在各自的道序号上分别向起始线和终止线方向自动搜索，便能得到整个数据体强反射对应的时间，层位与时间是一一对应的，一个层位是由很多个点组成的，每个点都由线号、道号和一个唯一的时间所确定。

所述步骤(3)中的最优小波为Morlet小波，小波的形态由振幅a_n、主频f_n、尺度σ_n、时间u_n和相位φ_n决定，Morlet小波的表达式如式(2)所示：

$g_{r_n}\left(t\right)=a_n\exp [-4\ln 2\frac{f_n^2{(t-u_n)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left(i(2\mathrm{\&pi;}{f}_n(t-u_n)+{\mathrm{\&phi;}}_n)\right)---\left(2\right).$

所述步骤(3)中的根据最大匹配投影原理对每道寻找最优小波具体包含以下步骤：

(31)对每道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)(这个在步骤(11)中是用于计算相位余弦数据体，此处再做一遍，是为了计算瞬时频率和瞬时相位)；计算复信号的瞬时频率和瞬时相位；其中u_n为对应道的时间t_i，f_n为对应时间的瞬时频率(就是上面计算得到的复信号的瞬时频率，就是每一道在时刻u_n，也就是对应的t_i上的瞬时频率和瞬时相位)，φ_n为对应时间的瞬时相位；

对数字信号进行希尔伯特变换、计算瞬时频率和瞬时相位函数，都是公知的常规方法，在三瞬分析中经常使用，其变换和计算方法在公知教科书中都有详细说明。

(32)在一组固定的u_n、f_n和φ_n情况下，通过计算最优化公式得到尺度σ_n；其中D＝{g_r(t)}_r∈Γ为时频原子字典，是函数R⁽ⁿ⁾s和的内积，

(33)利用所述最优化公式在局域内寻找r_n＝{σ_n，f_n，φ_n}的最优值，即公式达到最大时，对应的一组参数便是最优的r_n；搜索范围为[r_n-Δr，r_n+Δr]，其中Δr＝(Δσ，Δf，Δφ)，Δσ为尺度偏移量、Δf为频率偏移量、Δφ为相位偏移量；

(34)在上述两步的基础上，便可以预测到最优的时频原子的振幅为a_n，即 $a_n=\frac{|<<R^{\left(n\right)}s,g_{r_n}>>|}{<<g_{r_n},g_{r_n}>>};$

(35)基于上述步骤便匹配出能够决定Morlet小波形态的五个参数，根据Morlet小波的表达式(2)便可得到对应道的最优小波。

下面通过几个实施例来说明本发明方法的效果：

图2为任意一道地震信号和对应的瞬时相位余弦函数曲线。由于瞬时相位余弦独立于地震振幅，能够很好的反映地层的连续性，从而可以有效的分辨地震波的几何形态，能够用于识别层序的边界或地震相的外部几何形态和内部结构。由图中可看出瞬时相位余弦的主波峰、主波谷与地震信号的波峰、波谷对应较好，如图中短划线箭头所示；同时由于瞬时相位余弦独立于地震振幅，信号本身的能量差异对瞬时相位余弦函数值几乎没有影响，这就避免了地震信号能量对层位追踪的影响；原始地震信号中能量较小的同相轴也能在瞬时相位余弦函数中得到体现，如图中右边点划线箭头所示。

图3为任意的一条线地震剖面。根据工区实际情况，约1.5s处的强反射为目的层上覆煤层和下伏煤层造成的，目的层的有效反射被淹没在强反射之中；由于目的层较薄，地震资料分辨率较低，很难在常规地震资料中识别目的层段的砂体。所以，利用本发明去除由于煤层造成的强反射，突出目的层的有效反射。同时在此说明指定强反射位置为剖面中的强反射对应的时间。

图4为地震剖面对应的瞬时相位余弦剖面，从图中可以看出瞬时相位余弦剖面能够更清楚的显示层位走向，而且避免了地震反射能量对层位追踪的影响。尤其是像图3这类具有强反射的剖面，约1.46s以上的信号相对非常的微弱，很难判断同相轴的走向，而在瞬时相位余弦剖面上能够很清楚的看出1.46s以上地层的反射特征。那么，从瞬时相位余弦剖面中能够更合理的拾取出指定位置的层位走向。

图5为采用本发明去除的指定强反射的反射波，反射波沿自动拾取的地震层位；由于原始地震剖面中该同相轴各道之间能量差异较小，所以匹配最优小波的差异也比较小，对比图3中对应的强反射可知，拾取子波的能量与原始地震剖面相对应，拾取的最优小波是非常合理的。

图6为采用本发明去除煤层强反射后的效果图，从图中可看出，煤层强反射得到压制，目的层段的砂岩反射得到突出。采用去除强反射后的地震剖面进行储层预测能够更容易的圈定有利油气区域，能够够更精确的划分沉积相和沉积微相。

总之，本发明方法中拾取层位采用瞬时相位余弦剖面代替常规的地震剖面，避免了地震反射能量对层位追踪的影响，能够更准确的拾取指定强反射对应的层位。同时在局域内通过对Morlet小波进行伸缩、频移和相移得到一个原子库，根据最大匹配投影原理寻找强反射处的最优小波，能够更合理的拾取指定位置的地震反射和更合理的表征信号的局部结构特征。

本发明提供了一种自适应去除地震强反射的方法，属于油气田勘探开发中地震资料预处理领域。本发明首先计算整个三维地震数据体瞬时相位余弦；针对瞬时相位余弦数据体，采用互相关法自动拾取指定强反射的反射时间；沿着拾取的反射时间，根据最大匹配投影原理对每道寻找最优小波，每一道的地震数据减去对应的最优小波，便得到去除指定强反射的地震数据体。本发明方法用于去除掩盖目的层段储层信息的强反射，突出储层的反射信号，为储层预测提供更可靠的地震数据。

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

Claims (4)

1.一种自适应去除地震强反射的方法，其特征在于：所述方法包括：

(1)输入三维地震数据体，利用希尔伯特变换计算整个三维地震数据体中每道地震信号的瞬时相位余弦cosθ(t)，得到瞬时相位余弦数据体；

(2)针对瞬时相位余弦数据体，采用互相关法自动拾取指定点位置所在的层位，记录每道强反射层位的反射时间t_i；所述的指定点为用户预先设定的强反射在某一道上对应的时间；

(3)对输入的一道地震数据，根据该道自动拾取的层位时间，根据最大匹配投影原理匹配出最优小波；然后用该道的地震数据减去匹配的最优小波，得到该道的去除指定强反射后的地震数据；

(4)对整个三维地震数据中的每一道重复步骤(3)，得到每一道的去除指定强反射后的地震数据，所有道去除指定强反射后的地震数据构成去除指定强反射后的地震数据体。

2.根据权利要求1所述的自适应去除地震强反射的方法，其特征在于：所述步骤(1)具体包含以下步骤：

(11)对每一道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)；

(12)计算地震信号的相位余弦cosθ(t)，如式(1)所示：

$\cos \mathrm{\&theta;}\left(t\right)=\frac{x_r\left(t\right)}{\sqrt{x_r^2\left(t\right)+x_i^2\left(t\right)}}---\left(1\right)$

所有道地震信号的相位余弦就构成所述瞬时相位余弦数据体。

3.根据权利要求2所述的自适应去除地震强反射的方法，其特征在于：所述步骤(2)中采用互相关法自动拾取指定点位置所在的层位是这样实现的：确定互相关时窗长度和相对指定位置上下漂移距离，然后使用这两个参数自动拾取得到指定点位置所在的层位；

所述互相关时窗长度为至少包含一个完整的波形。

4.根据权利要求3所述的自适应去除地震强反射的方法，其特征在于：所述步骤(3)中的根据最大匹配投影原理匹配出最优小波具体包含以下步骤：

(31)对一道地震信号x_r(t)进行希尔伯特变换，得到实信号x_r(t)对应的虚信号x_i(t)；计算复信号的瞬时频率和瞬时相位；其中u_n为对应道的时间t_i，f_n为对应时间的瞬时频率，φ_n为对应时间的瞬时相位；

(32)在一组固定的u_n、f_n和φ_n情况下，通过计算最优化公式 $\arg \underset{g_{r_n}\&Element;D}{\max }\frac{|<<R^{\left(n\right)}s,g_{r_n}>>|}{\left|\right|g_{r_n}\left|\right|}$ 得到尺度σ_n；

其中D＝{g_r(t)}_r∈Γ为时频原子字典，是函数R⁽ⁿ⁾s和的内积， $\left|\right|g_{r_n}\left|\right|=\sqrt{<<g_{r_n},g_{r_n}>>};$

(33)利用所述最优化公式在局域内寻找r_n＝{σ_n，f_n，φ_n}的最优值，即公式达到最大时，对应的一组参数便是最优的r_n；搜索范围为[r_n-Δr，r_n+Δr]，其中Δr＝(Δσ，Δf，Δφ)，Δσ为尺度偏移量、Δf为频率偏移量、Δφ为相位偏移量；

(34)通过下式预测最优的时频原子的振幅为a_n，即

$a_n=\frac{|<<R^{\left(n\right)}s,g_{r_n}>>|}{<<g_{r_n},g_{r_n}>>};$

(35)基于步骤(31)至(34)得到决定Morlet小波形态的五个参数，根据Morlet小波的表达式(2)得到对应道的最优小波：

$g_{r_n}\left(t\right)=a_n\exp [-4\ln 2\frac{f_n^2{(t-u_n)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left(i(2\mathrm{\&pi;}{f}_n(t-u_n)+{\mathrm{\&phi;}}_n)\right)---\left(2\right).$