CN104459770B - 一种高维地震数据规则化方法 - Google Patents

Abstract

本发明提供了一种高维地震数据规则化方法，属于油气勘探开发中地震数据预处理领域。本方法包括：(1)输入原始三维炮道集；(2)利用高维傅里叶变换将时间空间域的所述原始三维炮道集变换到频率空间域中，然后进行分选得到五维炮道集；(3)对一个频率片进行最优化求解，得到无假频频谱；(4)判断是否达到迭代门槛值，如果是，则转入步骤(5)，如果否，则返回步骤(3)；(5)利用高维傅里叶反变换将所述无假频频谱变换到时间空间域，即得到规则炮道集数据；(6)输出所述规则炮道集数据。

Description

一种高维地震数据规则化方法

技术领域

本发明属于油气勘探开发中地震数据预处理领域，具体涉及一种高维地震数据规则化方法。

背景技术

在地震数据采集中，空间采样的不规则以及空间采样间隔过大(采集成本限制)导致地震数据的空间假频经常发生。由于偏移算子不能消除由空间采样引起的假频，以致这些假频能量以假象的形式出现在剖面中，这些假象会误导地震剖面的解释。因此为了实现真振幅成像，必须在偏移前通过地震数据插值等方法除去这些假频能量。在过去的十多年里，地球物理学家发展了很多解决地震数据非规则性和假频问题的方法，大体上可分为四类。描述如下：

第一类，基于褶积理论。这种方法的主要思想是，用地震数据的低频分量构造一个滤波器，并利用这个滤波器实现对地震数据的插值。属于这种类型的方法包括FX域预测滤波插值(Spitz，1991)，TX域预测误差滤波插值(Claerbout，1992)，FK域预测滤波插值(Gulunay，2003)。它们通常要求数据域的同相轴是线性的，只有满足线性假设条件，所构造的滤波器才有意义。这些方法的另一个要求是地震数据的空间采样方式必须是规则的，但是在数据采集过程中这个要求是很难实现的。

第二类，基于偏移算子和反偏移算子。这种方法先应用偏移算子将地震数据投影到成像域中，通过在成像域里完成地震数据的插值工作，然后将成像域的数据反偏移回数据域中。由Stolt(2002)和Trad(2003)介绍的方法都属于这种类型。这种方法的偏移算子种类繁多，包括正常时差校正(NMO)，倾角时差校正(DMO)，零偏移距偏移(MZO)，方位角时差校正(AMO)，叠前时间偏移(PSTM)，叠前深度偏移(PSDM)等。这类方法虽然可以很好的保持地震数据中同相轴的相位信息，但是振幅信息往往是失真的。

第三类，基于FK频谱估计。这种方法可以另外分为两小类：一种基于反演理论，另一种基于相关理论。基于反演理论的主流方法可参考A.J.W.Duijndam(1999)和Sacchi(2004，2007)的文章。这类方法的主要思想是，在反演理论的框架下，应用合适的正演算子和约束条件，估计地震数据的FK频谱。这类反演问题的模型空间和数据空间通常是FX域和FK域的地震数据，而正算子是傅立叶变换。这些方法的约束条件差异很大，不同的约束条件会产生不同的结果，因此约束条件是该方法的关键因素。总所周知，一个反问题经常会对应着多个解，不同的限制条件意味着在不同的解空间寻求结果，因此对于不同的约束条件而言解是完全不同的。另一类基于相关理论的方法包括抗泄露傅立叶变换(ALFT)(Sheng Xu，2005)，凸集投影(POCS)(Ray Abma，2006)，匹配追踪(MP)(Al i Ozbek，)以及其他一些与基于地震数据相关的方法。在相关过程中，这些方法将相关能量最大的基作为地震数据中能量最大的分量，然后选择这个最大能量作为FK谱估计的一个分量。实际上如果将相关理论看作解反演问题的一种方法，那么每次选择的基就是梯度，而相关能量就是这个反演问题的步长。如果数据中的假频不严重，那么这类方法能得到很好的结果。但是如果假频严重，最大相关能量的基可能就不再是地震数据有效信号的分量而是假频，则此方法得到的插值结果就很差。

最后一种，基于矩阵分解。如果将地震数据看作一个高维矩阵，那么在奇异值分解(SVD)中，数据插值可以通过插值奇异向量得到。由Sacchi(2011)提出的方法就属于这一类。如果数据域中的同相轴是线性的，那这种方法将得到更好的结果，因为如果同相轴是线性的，那么大奇异值的所对应的奇异向量的数量将大大减少。但是其缺陷也是显而易见的，因为地震数据的随机性很强，因此在奇异值分解时就会产生奇异值从而导致矩阵误解。

上述四类方法中，最后两类是当今地震数据插值中最为流行的，因为它们都具备高维地震数据插值的能力，并且插值效果也比较理想。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种高维地震数据规则化方法，解决地震数据采集中，空间采样的不规则以及空间采样间隔过大(采集成本限制)导致地震数据的空间假频的问题，这些假频能量以假象的形式出现在剖面中，这些假象会误导地震剖面的解释。本方法是在反演过程中以高维傅里叶变换为正演算子，数据线性拉东变换谱为约束的高维数据规则化方法。

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

一种高维地震数据规则化方法，所述方法包括：

(1)输入原始三维炮道集；

(2)利用高维傅里叶变换将时间空间域的所述原始三维炮道集变换到频率空间域中，然后进行分选得到五维炮道集；

(3)对一个频率片(该频率片是由时间域通过傅里叶变换到频率域，有多少个时间采样点就有多少个频率片。)进行最优化求解，得到无假频频谱；

(4)判断是否达到迭代门槛值，如果是，则转入步骤(5)，如果否，则返回步骤(3)；

(5)利用高维傅里叶反变换将所述无假频频谱变换到时间空间域，即得到规则炮道集数据；

(6)输出所述规则炮道集数据。

所述步骤(3)是采用共轭梯度算法求解下面的公式实现的：

这个公式就是每个频率片要求解的，即公式是频率域的基础上求解。

其中，是无假频的规则数据频谱，即无假频频谱，m_prior是先验模型，C_M为数据的Radon谱，C_D为一般约束条件(最小二乘反演的约束条件，一般给0.01)，d_obs为观测数据(观测数据就是野外得到的炮集数据经过傅里叶变到频率域的数据。)，A表示采样算子矩阵(采样算子就是根据野外观测系统定义的一个矩阵，采集到的地方为1，没有采集到的地方为0.)，A^H为A的共轭转置矩阵，m_prior＝0。

所述数据的Radon谱是这样得到的：

利用步骤(1)中的炮道集数据中同相轴的局部线性，进行线性拉东变换，得到拉东谱，具体实现公式如下：

R(f，p)=∫d(f，x)e^-2πfpxdx

其中，f，p，x分别表示频率、射线参数和空间采样位置，d(f，x)是原始数据的频率空间表示，R(f，p)为所述数据的Radon谱。

所述步骤(4)中的迭代门槛值设置成10^-6(即10的负六次方)。

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

1)本发明方法使用了高精度，高效率的Intel编译库中高维傅里叶正反变换，其优势在于进行同等效果的傅里叶时的速度更加的快捷稳健，提高计算效率；

2)本发明方法利用原始数据同相轴线性拉东谱作为反演约束条件，其优势在于利用同相轴的拉东谱约束算法更加稳健，反演的结果更加的接近真解。

附图说明

图1本发明实施例中的规则采样的五个同相轴理论合成数据；

图2本发明实施例中图1的傅里叶变换规则数据频谱；

图3本发明实施例中图1随机抽取三分之一道后非规则采样数据；

图4本发明实施例中图3的傅里叶变换频谱；

图5本发明实施例中的通过本方法得到的频谱反变换得到的规则数据；

图6本发明实施例中的通过本方法得到的规则数据频谱；

图7本发明实施例中的合成的规则炮道集；

图8本发明实施例中的图7随机抽取三分之一道得到的非规则炮道集；

图9本发明实施例中的通过本方法得到的规则炮道集；

图10本发明实施例中的规则后的炮道集与原始规则炮道集的残差；

图11本发明实施例中随机抽取三维盐丘模型三分之一炮道集以及道集内随机抽取三分之一道得到的非规则炮道集；

图12本发明实施例中通过本方法得到的规则炮道集；

图13本发明实施例中图11局部放大非规则炮道集；

图14本发明实施例中图12局部放大规则炮道集；

图15本发明实施例中某野外实际非规则炮道集；

图16本发明实施例中通过本方法差值充零得到的炮道集；

图17本发明实施例中通过本方法得到的规则炮道集；

图18本发明方法的步骤框图。

具体实施方式

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

本发明通过对具有一定信噪比(即数据中没有强能量的野值和线性面波能量)的炮道集数据进行高维傅里叶变换，然后利用炮道集数据同相轴的局部线性拉东变换频谱作为约束，通过反演将FK谱(频率波数谱)的假频分量减弱得到规则数据的频谱，最后再利用高维傅里叶反变换得到规则数据。

如图18所示，所述方法具体包括以下步骤：

(1)准备具有高分辨率、高信噪比的原始地震炮集数据

(2)利用Intel编译库中高维傅里叶变换将时间空间域的炮集数据变换到频率空间域中，这个过程中将地震道头信息的炮点坐标记录下来形成道头信息索引；野外采集的地震数据是时间和空间的3维数据，而在数据规则化中要把3维数据转换到5维，即频率、x、y、hx和hy，hx和hy分别是x和y方向的偏移距，这样的数据称为高维数据。5维数据是通过高维地震数据插值产生的：输入的是原始3D炮道集，进入程序后进行傅里叶变换到频率空间域，之后再进行分选到5D进行插值，即x、y，hx，hy，t，这样得到的效果更加的接近实际采集的情况；

(3)利用(1)中原始炮道集数据中同相轴的局部线性，进行线性拉东变换，得到拉东谱，具体实现公式如下所示：

R(f，p)=∫d(f，x)e^-2πfpxdx

其中，f，p，x分别表示频率、射线参数和空间采样位置，d(f，x)是原始数据的频率空间表示，R(f，p)为数据的Radon谱；

(4)反演过程中将拉东谱作为约束条件，对每个频率片进行反演得到规则数据的无假频频谱，具体如下：

其中是这个问题的最佳模型即无假频的规则数据频谱，m_prior是先验模型，C_M为数据的Radon谱，C_D为一般约束条件，d_obs为观测数据，A表示采样算子矩阵，A^H为A的共轭转置矩阵。一般而言m_prior=0，解最优化问题用共轭梯度算法求解。

(5)利用Intel编译库中高维傅里叶反变换将(4)中反演得到的无假频频谱变换到时间空间域得到规则采样的炮道集数据，即规则炮道集数据。

图1到图6是2维理论数据验证方法理论；图7到图10是3维理论数据进行算法测试；图11到图17是利用3维盐丘模型进行算法测试的结果。

本发明属于第三种类型中基于反演理论的方法。该方法的正过程即规则且无假频地震数据的傅立叶谱m(f，k)是本发明想要的模型。存在空间假频的且非规则的野外数据d(f，x)是已知的。则正演过程可以表示为d(f，x)=A(x，k)m(f，k)向量d(f，x)，m(f，k)及矩阵A(x，k)的元素可以相应地作如下表达(下面3个公式表示每个矩阵的具体元素的含义，即原始观测数据，反演模型空间矩阵元素和采样矩阵算子的元素)：

d_i(f，x)=d(f，x_i)

m_i(f，k)=m(f，k_i)

a_mn=exp(i2πx_m·k_n)

其中d(f，x_i)是野外数据，x_i是其空间坐标。m(f，k_i)是要估计模型的分量，其坐标是k_i。x_m·k_n是x_m和k_n的内积。f是地震数据的频率。这个正过程与Dui jndam等人(1999)给出的定义相同，反过程利用了地震数据的拉东谱能量分布构造一个椭圆范数做为约束条件，所构造的椭圆范数能提高FK谱的相对能量分布。如果同相轴是线性的，那么拉东变换可以在FK域中有效地聚集这些同相轴的能量，相应地，FK谱中的假频分量将相对地被减弱。因此，在这样的约束条件下插值过程会变得更加稳定，反演问题的解也会变得更加稀疏且更合逻辑。

如果数据采样均匀，每个采样点在数据空间的元素中的地位是等价的，并且是相互独立的。这样数据的协方差矩阵主对角线上的值相等。其余非对角元素为0，即是一个单位阵I。这样起能量分布就是一个球面，实际数据采样是不等间隔的，这样数据的协方差矩阵主对角线值不再相等，也无法用单位阵表示。协方差矩阵不仅要描述各个采样点处能量的强弱关系，还要描述不同采样点之间的相互关系，协方差矩阵的值主要由观测系统的形式决定。首先，协方差矩阵必须能反映各个采样点处如下的能量关系：在采样比较密集的方向上，对应着比较强的能量，因此在协方差矩阵对角线上与之相应的元素的值应该比较大；在采样比较稀疏的方向上，对应着比较弱的能量，因此在协方差矩阵对角线上与之相应的元素的值应该比较小。协方差矩阵主对角线上的元素是采样点的自相关，表示各个采样点处能量的大小；而非主对角线上的元素是不同采样点之间互相关，表示不同采样点之间相关性的大小。在正常情况下，均假设不同点处采集到的地震数据是不相关的，即各个采样点之间是互相独立的。因此，协方差矩阵非主对角线上的元素的值均默认为0。如果每个检波器所接收到的地震数据是相互独立的，那么协方差矩阵非对角线上的元素的值均为0，协方差矩阵变成了一个对角矩阵。即数据空间的单位圆对应的不再是一个球面，而是一个椭球面。

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

Claims (3)

1.一种高维地震数据规则化方法，其特征在于：所述方法包括：

(1)输入原始三维炮道集；

(2)利用高维傅里叶变换将时间空间域的所述原始三维炮道集变换到频率空间域中，然后进行分选得到五维炮道集；

(3)对一个频率片进行最优化求解，得到无假频频谱；

(4)判断是否达到迭代门槛值，如果是，则转入步骤(5)，如果否，则返回步骤(3)；

(5)利用高维傅里叶反变换将所述无假频频谱变换到时间空间域，即得到规则炮道集数据；

(6)输出所述规则炮道集数据,

所述步骤(3)是采用共轭梯度算法求解下面的公式实现的：

$\begin{array}{c}\tilde{m}=m_{prior}+{(A^H{C}_D^{-1}A+C_M^{-1})}^{-1}{A}^H{C}_D^{-1}(d_{obs}-{\mathrm{Am}}_{prior})\\ =m_{prior}+C_M{A}^H{({\mathrm{AC}}_M{A}^H+C_D)}^{-1}(d_{obs}-{\mathrm{Am}}_{prior})\end{array}$

其中，是无假频的规则数据频谱，即无假频频谱，m_prior是先验模型，C_M为数据的Radon谱，C_D为一般约束条件，d_obs为观测数据，A表示采样算子矩阵，A^H为A的共轭转置矩阵，m_prior＝0。

2.根据权利要求1所述的高维地震数据规则化方法，其特征在于：所述数据的Radon谱是这样得到的：

利用步骤(1)中的炮道集数据中同相轴的局部线性，进行线性拉东变换，得到拉东谱，具体实现公式如下：

R(f,p)＝∫d(f,x)e^-2πfpxdx

其中，f,p,x分别表示频率、射线参数和空间采样位置，d(f,x)是原始数据的频率空间表示，R(f,p)为所述数据的Radon谱。

3.根据权利要求1所述的高维地震数据规则化方法，其特征在于：所述步骤(4)中的迭代门槛值设置成10^-6。