不规则地震数据的五维插值处理方法及装置的制造方法
【技术领域】
[0001] 本发明涉及地震勘探技术领域，特别涉及一种不规则地震数据的五维插值处理方 法及装置。
【背景技术】
[0002] 地震勘探主要分为地震资料采集、处理与解释三个阶段。在地震资料采集阶段，无 论是陆地勘探还是海洋勘探，受各种复杂因素的影响，实际采集资料的空间采样只能是近 似规则的，甚至是不规则的。例如，对于陆地勘探，在城区、公路、河网等地表条件复杂的地 区，无法布置规则的观测系统，必须进行不规则采集;对于海洋勘探，受洋流的影响，实际电 缆测线存在羽状漂移的现象，同样难以保证规则的空间采样。在地震资料处理阶段，许多数 据处理算法(如偏移算法和多次波去除算法)都需要、或者至少得益于规则空间采样的地震 数据。此外，在进行连片或者时移地震数据处理时，不同批次采集的地震数据的观测系统参 数一般是不同的，这也涉及到地震数据的空间采样规则化问题。因此，利用特殊的插值算法 解决地震数据空间采样不规则问题是十分必要的。
[0003] 目前，大多数地震数据插值算法还仅局限于三维插值，但野外采集的地震数据本 质上是五维坐标的函数，两维用于确定炮点空间位置，两维用于确定检波点空间位置，还有 一维用于确定采样点时间，发展五维插值算法可以更充分地利用采集的地震数据，进而获 得更好插值结果。与三维插值算法相比，五维插值算法面临的首要难题就是巨大的数据量 和计算量。Trad(2009)将Liu和Sacchi(2004)的三维最小加权范数插值算法推广至五维，算 法实现时，所有的矩阵与向量乘积运算都可以利用快速傅里叶变换(FFT)完成，因而保证了 计算效率，但算法要求地震数据的空间采样是规则的，无法用于不规则采集地震数据，具有 很大的局限性。Jin(2010)基于不规则空间采样假设，提出了基于衰减最小范数傅里叶反演 的五维地震数据插值算法，算法采用不等间隔快速傅里叶变换(NFFT)实现频繁的矩阵与向 量乘积运算，在一定程度上改善了计算效率，但由于高维数据NFFT的计算效率远不及FFT， 并且NFFT本身是一种近似算法，Jin的方法在计算效率方面仍有待改善。
【发明内容】
[0004] 本发明实施例提供了一种不规则地震数据的五维插值处理方法，用以提高地震数 据处理的效率，该方法包括：
[0005] 将不规则地震数据变换到频率域和空间域；
[0006] 对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换 算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；
[0007] 将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间 域。
[0008] 本发明实施例提供了一种不规则地震数据的五维插值处理装置，用以提高地震数 据处理的效率，该装置包括：
[0009] 地震数据变换模块，用于将不规则地震数据变换到频率域和空间域；
[0010] 五维插值处理模块，用于对于每一个变换到频率域和空间域的地震数据的频率切 片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的 频率波数域数据；
[0011]五维插值数据变换模块，用于将五维插值后的频率波数域数据变换到频率域和空 间域，再变换到时间域和空间域。
[0012] 与现有技术相比较，本发明实施例提供的技术方案，通过将不规则地震数据变换 到频率域和空间域;对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速 傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域 数据;将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域， 实现了对不规则地震数据的快速处理，整个计算过程仅涉及少量NFFT运算，绝大部分的矩 阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，提高了地震数据处理的效 率，具有重要的实际应用价值。
【附图说明】
[0013] 此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不 构成对本发明的限定。在附图中：
[0014] 图1是本发明实施例中不规则地震数据的五维插值处理方法的流程示意图；
[0015]图2是本发明实施例中插值前的地震数据观测系统示意图；
[0016] 图3是本发明实施例中插值前的同一个CMP面元的地震数据观测系统示意图；
[0017] 图4是本发明实施例中与图3对应的插值后的的同一个CMP面元的地震数据观测系 统示意图；
[0018] 图5是本发明实施例中采用不同快速算法的运算时间对比图表示意图；
[0019]图6是本发明实施例中与图3对应的插值前的CMP道集地震数据示意图；
[0020] 图7是本发明实施例中与图4对应的插值后的CMP道集地震数据；
[0021] 图8是本发明实施例中不规则地震数据的五维插值处理装置的结构示意图。
【具体实施方式】
[0022] 为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对 本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并 不作为对本发明的限定。
[0023]本发明提出了一种新的迭代方案，整个计算过程仅涉及少量NFFT运算，绝大部分 的矩阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，具有重要的实际应用价 值。下面结合附图1至8对该方案进行详细介绍如下。
[0024] 图1是本发明实施例中不规则地震数据的五维插值处理方法的流程示意图，如图1 所示，该处理方法包括如下步骤：
[0025] 步骤101:将不规则地震数据变换到频率域和空间域；
[0026] 步骤102:对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅 里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数 据；
[0027] 步骤103:将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间 域和空间域。
[0028] 与现有技术相比较，本发明实施例提供的技术方案，通过将不规则地震数据变换 到频率域和空间域;对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速 傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域 数据;将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域， 实现了对不规则地震数据的快速处理，整个计算过程仅涉及少量NFFT运算，绝大部分的矩 阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，提高了地震数据处理的效 率，具有重要的实际应用价值。
[0029] 在一个实施例中，上述步骤101之前可以包括：
[0030] 不规则地震数据进行去噪、静校正和线性动校正处理；
[0031] 从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据；
[0032] 具体实施时，不规则地震数据进行去噪、静校正和线性动校正处理，可以包括:对 采集的地震数据进行前期处理操作，包括去噪、静校正、线性动校正等。
[0033] 具体实施时，从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数 据，可以包括:抽取想要进行五维插值的地震数据為。
[0034]上述tit为时间采样坐标，下标i t = {0，1，2，…，Nt_ 1}，即时间方向有Nt个采样点， 设时间采样间隔为At,则tit = it· At;
[0035]上述^为空间采样坐标，下标4={0，1，2，-_，仏-1}，即空间方向有仏个采样点， I=氐,λ' ,毛,),t = ((i" )/f, (qf, , (? )")是四维空间中的一个点，其中元，.》A 一种一般表示形式，根据实际需要，可以分别代表炮点x坐标、炮点y坐标、检波点x坐标、检 波点y坐标，也可以分别代表中心点X坐标、中心点y坐标，炮检距，方位角。
[0036]上述步骤101可以包括:将抽取的待五维插值处理的不规则地震数据变换到频率 域和空间域。
[0037]具体实施时，将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域 可以包括：
[0038]对^ = ((? )" ,(.? )" ,(.?丄,(毛),、)中的每一维坐标进行归一化变换，变换结果记为 Xix =(( X0 ) ix，（ XI) ix，（ X2 ) ix，（ X3 ) ix )，变换后的地震数据表示为 J(元。
[0039]
其中，η={0,1,2,3}表示不 同的维，对于任意一维坐标(尽?Λ，ix= {〇, 1，2,…，ΝΧ-1}，取其最大值和最小值，分别记为 vmaxn和vminn
，nn是根据实际需 要而设定的插值后的坐标点数，插值后各个空间维的采样间隔为
[0040]利用快速傅里叶变换将地震数据d(Xix，tit)变换到频率域d(Xix，ω ?ω)。
[0041 ] 上述d(xix, ω?ω)与d(xix，tit)满足如下关系式，
,? , ω?ω = ?ω · Δ ω ,
，i ω = {〇，1，2，···，Nt_l}。这里需要 保证地震数据d(Xlx，tlt)满足采样定理，设地震数据d(Xlx，t lt)的有效信号能量刚好处于 0 ~ 的频带内，其中则需要满足条件i comax < fl〇or(Nt/2) 〇
[0042] 生成不等间隔离散逆傅里叶变换矩阵F(Xix，kik)，其中，丨1={0，1，2，-_，队_1}为矩 阵的行索弓丨，ik={0，l，2，···，η〇η?Π 2η3-1}为矩阵的列索引。
[0043]
，其中，kik=(ko(iQ)，ki(ii)，k2(i2)，k3(i3))ik是四维空间 中的一个点，kn(in) = in-Cn，Cn=fl〇〇r((nn-l)/2)，in= {〇, 1，2,…,ηη-1}，n = 〇, 1，2,3,函数 fl〇〇r()表示向下取整，kikXix表示四维向量kik与Xix的内积。并且ik与i〇、ii、i2、i3保持关系 式ik=((i〇 · m+ii) · n2+i2) · M+i3,其中，η〇、ηι、Π 2、η3的定义见步骤(3)。
[0044] 在一个实施例中，上述步骤102可以包括：
[0045] 对每一个频率切片i ω={〇，1，2，···4 comax}，按照如下公式进行处理：
[0046] biM = FHdiM；
[0047] yiGj = CG(bicj，Wicj);
[0048] Uicj = conj(yicj) · *yicj;
[0049] wicj+i = (sqrt(UiGj/sum(UiGj) )+wicj)/2;
[0050] 其中，i comax是根据实际需要选取的频率上限索引，选取原则是使得地震数据d (Xlx，t lt)的有效信号能量刚好处于?的频带内；
[0051 ] wo为一个包含nonimm个元素的列向量，且wo的每一个元素都等 中，wo是Wiu中下角标i ω =〇时的情况；
[0052] F=F(xix，kik)，为一个Νχ行、norumm列的矩阵;其中，F(xix，kik)为不等间隔离散逆 傅里叶变换矩阵，ix={〇，l，2，···，NX-1}为矩阵的行索弓丨，ik={0，l，2，···，n〇mn2n3-l}为矩
其中，kik= (k〇(i〇)，ki(ii)，k2(i2)，k3(i3) )ik是四维空间 中的一个点，kn(in) = in-Cn，Cn=fl〇〇r((nn-l)/2)，in= {〇, 1，2, ··· ,ηη-1}，n = 〇, 1，2,3,函数 fl〇〇r()表示向下取整，kikXix表示四维向量kik与Xix的内积，ik与i〇、ii、i2、i3保持关系式ik = ((?〇 * η?+??) · Π 2+?2) * Π3+?3；
[0053] 4& = [^/(\。,％.,).办％丨,1)乂办:,以,,..>^ 向量；
[0054] 函数conj〇表示对输入向量的每一个元素求复共辄，并输出一个向量;运算符· * 表示两个向量的对应元素相乘，并得到一个向量；函数sqrt〇表示对输入向量的每一个元 素开平方，并输出一个向量;函数sum()表示对输入向量的所有元素求和，并输出一个标量； 运算符H表示求矩阵的共辄转置；
[0055] 函数y = CG(b，w)表示如下预条件共辄梯度算法，①设定常量，包括:迭代终止相对 容差tol = l(T4,最大迭代次数maxit = 30，A = FHF ;②计算迭代初始变