CN105549078B - 不规则地震数据的五维插值处理方法及装置 - Google Patents

Abstract

本发明公开了一种不规则地震数据的五维插值处理方法及装置，该方法包括：将不规则地震数据变换到频率域和空间域；对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域。上述技术方案实现了对不规则地震数据的快速处理，整个处理过程绝大部分的矩阵与向量乘积运算都利用快速傅里叶变换FFT算法实现，提高了地震数据处理的效率。

Description

不规则地震数据的五维插值处理方法及装置

技术领域

本发明涉及地震勘探技术领域，特别涉及一种不规则地震数据的五维插值处理方法及装置。

背景技术

地震勘探主要分为地震资料采集、处理与解释三个阶段。在地震资料采集阶段，无论是陆地勘探还是海洋勘探，受各种复杂因素的影响，实际采集资料的空间采样只能是近似规则的，甚至是不规则的。例如，对于陆地勘探，在城区、公路、河网等地表条件复杂的地区，无法布置规则的观测***，必须进行不规则采集；对于海洋勘探，受洋流的影响，实际电缆测线存在羽状漂移的现象，同样难以保证规则的空间采样。在地震资料处理阶段，许多数据处理算法(如偏移算法和多次波去除算法)都需要、或者至少得益于规则空间采样的地震数据。此外，在进行连片或者时移地震数据处理时，不同批次采集的地震数据的观测***参数一般是不同的，这也涉及到地震数据的空间采样规则化问题。因此，利用特殊的插值算法解决地震数据空间采样不规则问题是十分必要的。

目前，大多数地震数据插值算法还仅局限于三维插值，但野外采集的地震数据本质上是五维坐标的函数，两维用于确定炮点空间位置，两维用于确定检波点空间位置，还有一维用于确定采样点时间，发展五维插值算法可以更充分地利用采集的地震数据，进而获得更好插值结果。与三维插值算法相比，五维插值算法面临的首要难题就是巨大的数据量和计算量。Trad(2009)将Liu和Sacchi(2004)的三维最小加权范数插值算法推广至五维，算法实现时，所有的矩阵与向量乘积运算都可以利用快速傅里叶变换(FFT)完成，因而保证了计算效率，但算法要求地震数据的空间采样是规则的，无法用于不规则采集地震数据，具有很大的局限性。Jin(2010)基于不规则空间采样假设，提出了基于衰减最小范数傅里叶反演的五维地震数据插值算法，算法采用不等间隔快速傅里叶变换(NFFT)实现频繁的矩阵与向量乘积运算，在一定程度上改善了计算效率，但由于高维数据NFFT的计算效率远不及FFT，并且NFFT本身是一种近似算法，Jin的方法在计算效率方面仍有待改善。

发明内容

本发明实施例提供了一种不规则地震数据的五维插值处理方法，用以提高地震数据处理的效率，该方法包括：

将不规则地震数据变换到频率域和空间域；

对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域。

本发明实施例提供了一种不规则地震数据的五维插值处理装置，用以提高地震数据处理的效率，该装置包括：

地震数据变换模块，用于将不规则地震数据变换到频率域和空间域；

五维插值处理模块，用于对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

五维插值数据变换模块，用于将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域。

与现有技术相比较，本发明实施例提供的技术方案，通过将不规则地震数据变换到频率域和空间域；对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域，实现了对不规则地震数据的快速处理，整个计算过程仅涉及少量NFFT运算，绝大部分的矩阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，提高了地震数据处理的效率，具有重要的实际应用价值。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1是本发明实施例中不规则地震数据的五维插值处理方法的流程示意图；

图2是本发明实施例中插值前的地震数据观测***示意图；

图3是本发明实施例中插值前的同一个CMP面元的地震数据观测***示意图；

图4是本发明实施例中与图3对应的插值后的的同一个CMP面元的地震数据观测***示意图；

图5是本发明实施例中采用不同快速算法的运算时间对比图表示意图；

图6是本发明实施例中与图3对应的插值前的CMP道集地震数据示意图；

图7是本发明实施例中与图4对应的插值后的CMP道集地震数据；

图8是本发明实施例中不规则地震数据的五维插值处理装置的结构示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

本发明提出了一种新的迭代方案，整个计算过程仅涉及少量NFFT运算，绝大部分的矩阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，具有重要的实际应用价值。下面结合附图1至8对该方案进行详细介绍如下。

图1是本发明实施例中不规则地震数据的五维插值处理方法的流程示意图，如图1所示，该处理方法包括如下步骤：

步骤101：将不规则地震数据变换到频率域和空间域；

步骤102：对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

步骤103：将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域。

与现有技术相比较，本发明实施例提供的技术方案，通过将不规则地震数据变换到频率域和空间域；对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域，实现了对不规则地震数据的快速处理，整个计算过程仅涉及少量NFFT运算，绝大部分的矩阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，提高了地震数据处理的效率，具有重要的实际应用价值。

在一个实施例中，上述步骤101之前可以包括：

不规则地震数据进行去噪、静校正和线性动校正处理；

从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据；

具体实施时，不规则地震数据进行去噪、静校正和线性动校正处理，可以包括：对采集的地震数据进行前期处理操作，包括去噪、静校正、线性动校正等。

具体实施时，从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据，可以包括：抽取想要进行五维插值的地震数据

上述t_it为时间采样坐标，下标it的取值范围是0,1,2,…,N_t-1，即时间方向有N_t个采样点，设时间采样间隔为Δt,则t_it＝it·Δt；

上述为空间采样坐标，下标ix的取值范围是0,1,2,…,N_x-1，即空间方向有N_x个采样点，是四维空间中的一个点，其中是一种一般表示形式，根据实际需要，可以分别代表炮点x坐标、炮点y坐标、检波点x坐标、检波点y坐标，也可以分别代表中心点x坐标、中心点y坐标，炮检距，方位角。

上述步骤101可以包括：将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域。

具体实施时，将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域可以包括：

对中的每一维坐标进行归一化变换，变换结果记为x_ix＝((x₀)_ix,(x₁)_ix,(x₂)_ix,(x₃)_ix)，变换后的地震数据表示为

上述归一化变换方法为其中，η＝{0,1,2,3}表示不同的维，对于任意一维坐标取其最大值和最小值，分别记为vmax_η和vmin_η，具体地和n_η是根据实际需要而设定的插值后的坐标点数，插值后各个空间维的采样间隔为

利用快速傅里叶变换将地震数据d(x_ix,t_it)变换到频率域d(x_ix,ω_iω)。

上述d(x_ix,ω_iω)与d(x_ix,t_it)满足如下关系式，其中，t_it＝it·Δt，ω_iω＝iω·Δω，iω＝{0,1,2,…,N_t-1}。这里需要保证地震数据d(x_ix,t_it)满足采样定理，设地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内，其中则需要满足条件iω_max≤floor(N_t/2)。

生成不等间隔离散逆傅里叶变换矩阵F(x_ix,k_ik)，其中，ix＝{0,1,2,…,N_x-1}为矩阵的行索引，ik＝{0,1,2,…,n₀n₁n₂n₃-1}为矩阵的列索引。

上述其中，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，η＝0,1,2,3，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积。并且ik与 保持关系式其中，n₀、n₁、n₂、n₃的定义见步骤(3)。

在一个实施例中，上述步骤102可以包括：

对每一个频率切片iω＝{0,1,2,…,iω_max}，按照如下公式进行处理：

b_iω＝F^Hd_iω；

y_iω＝CG(b_iω,w_iω)；

u_iω＝conj(y_iω)·*y_iω；

w_iω+1＝(sqrt(u_iω/sum(u_iω))+w_iω)/2；

其中，iω_max是根据实际需要选取的频率上限索引，选取原则是使得地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内；

w₀为一个包含n₀n₁n₂n₃个元素的列向量，且w₀的每一个元素都等于其中，w₀是w_iω中下角标iω＝0时的情况；

F＝F(x_ix,k_ik)，为一个N_x行、n₀n₁n₂n₃列的矩阵；其中，F(x_ix,k_ik)为不等间隔离散逆傅里叶变换矩阵，ix＝{0,1,2,…,N_x-1}为矩阵的行索引，ik＝{0,1,2,…,n₀n₁n₂n₃-1}为矩阵的列索引，其中，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积，ik与保持关系式

为一个包含N_x个元素的列向量；

函数conj()表示对输入向量的每一个元素求复共轭，并输出一个向量；运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和，并输出一个标量；运算符^H表示求矩阵的共轭转置；

函数y＝CG(b,w)表示如下预条件共轭梯度算法，①设定常量，包括：迭代终止相对容差tol＝10^-4，最大迭代次数maxit＝30，A＝F^HF；②计算迭代初始变量，包括：r₀＝b，ρ_-1＝1；③计算迭代终止绝对容差对于i_cg＝{0,1,2,…,maxit-1}，执行如下处理操作：

如果或者i_cg＝maxit，跳出循环，返回否则继续循环:

其中，为标量，表示标量分别与向量的每一个元素相乘，并返回一个向量，表示标量分别与向量的每一个元素相乘，并返回一个向量；

在一个实施例中，所述函数y＝CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q＝Ap，其中A＝F^HF具有多级分块Toeplitz矩阵结构，采用如下快速算法计算q＝Ap：

①设N₀＝n₁n₂n₃、N₁＝n₂n₃、N₂＝n₃、N₃＝1，m_η＝2n_η-1,η＝0,1,2,3，按如下方式生成长度为m_η的向量和其中η＝0,1,2,3：

②设向量是一个四维数组的向量化表示，i_a为向量a的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，此处表示i_a是的函数，以下均采用了类似表达，按如下方式对a的元素赋值，初始化变量

对于执行如下循环(1)至(12)

(1)

(2)

(3)对于执行如下循环(4)至(12)

(4)

(5)

(6)对于执行如下循环(7)至(12)

(7)

(8)

(9)对于执行如下循环(10)至(12)

(10)

(11)

(12)

其中，表示向量的第个元素，η＝0,1,2,3，表示向量a的第个元素，表示矩阵A的第行、列的元素；

③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量

其中，向量是一个四维数组的向量化表示，为向量的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，函数FFT₄()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组；

④设向量是一个四维数组的向量化表示，i_p为向量p的索引，为向量p对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

设向量是一个四维数组的向量化表示，i_p′为向量p′的索引，为向量p′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；按照如下方式对向量p′的各个元素赋值：

⑤计算其中，函数FFT₄()和IFFT₄()分别表示对四维数组做四维快速傅里叶变换和逆变换，并返回一个四维数组；运算符·*表示两个四维数组的对应元素相乘，并得到一个四维数组；

向量是一个四维数组的向量化表示，i_q′为向量q′的索引，为向量q′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；

⑥设向量是一个四维数组的向量化表示，i_q为向量q的索引，为向量q对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，按照如下方式对向量q的各个元素赋值,其中0≤i_η≤n_η-1，η＝0,1,2,3。

具体实施时，整个方法流程中，以上第①～③步只需要计算一次，运算量基本可以忽略，因此，如第⑤步所述，矩阵与向量乘积运算q＝Ap的计算量主要包括一次四维快速傅里叶变换和一次四维快速傅里叶逆变换。如图5所示，本发明人实施例提供的技术方案提高了不规则地震数据处理的效率。

在一个实施例中，还包括：对y_iω执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组y的元素位置进行调整，并返回一个四维数组其中：

是一个四维数组的向量化表示，i_y为向量y的索引，为向量y对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

其中c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，函数mod(v₁,v₂)表示v₁对v₂求0～v₂-1之间的余数，函数floor()表示向下取整。

在一个实施例中，还包括：对执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组做四维快速傅里叶逆变换，并返回一个四维数组s，其中，

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，i_s为向量s的索引，为向量s对应的四维数组索引，i_s的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

在一个实施例中，上述步骤103可以包括：

利用s_iω生成插值后的频率域五维数据其中，为插值后的空间采样坐标，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，iω的取值上限由iω_max变为N_t-1，以满足后续进行快速傅里叶逆变换的需求，

其中，i_η的取值范围是0,1,2,…,n_η-1，vmin_η为插值后各个空间维坐标的最小值，为插值后各个空间维的采样间隔，η＝0,1,2,3，ix与i₀、i₁、i₂、i₃保持关系式ix＝((i₀·n₁+i₁)·n₂+i₂)·n₃+i₃；生成方式如下，用向量代表的频率坐标索引为iω的一部分数据，

则

其中，函数conj()表示对输入向量的每一个元素求复共轭，z表示由n₀n₁n₂n₃个0元素构成的向量；

利用快速傅里叶逆变换将变换的时间域，得到插值后的五维数据

其中，与满足如下关系式，其中，t_it＝it·Δt，it的取值范围是0,1,2,…,N_t-1，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，N_t为采样点个数；Δt为时间采样间隔。

基于同一发明构思，本发明实施例中还提供了一种不规则地震数据的五维插值处理装置，如下面的实施例所述。由于不规则地震数据的五维插值处理装置问题的原理与不规则地震数据的五维插值处理方法相似，因此不规则地震数据的五维插值处理装置的实施可以参见不规则地震数据的五维插值处理方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图8是本发明实施例中不规则地震数据的五维插值处理装置的结构示意图，如图8所示，该装置包括：

地震数据变换模块10，用于将不规则地震数据变换到频率域和空间域；

五维插值处理模块20，用于对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

五维插值数据变换模块30，用于将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域。

在一个实例中，还包括：

地震数据预处理模块，用于对不规则地震数据进行去噪、静校正和线性动校正处理；

地震数据抽取模块，用于从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据；

所述地震数据变换模块具体用于：将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域。

在一个实例中，五维插值处理模块具体用于：

对每一个频率切片iω＝{0,1,2,…,iω_max}，按照如下公式进行处理：

b_iω＝F^Hd_iω；

y_iω＝CG(b_iω,w_iω)；

u_iω＝conj(y_iω)·*y_iω；

w_iω+1＝(sqrt(u_iω/sum(u_iω))+w_iω)/2；

其中，iω_max是根据实际需要选取的频率上限索引，选取原则是使得地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内，其中

w₀为一个包含n₀n₁n₂n₃个元素的列向量，且w₀的每一个元素都等于其中，w₀是w_iω中下角标iω＝0时的情况；

F＝F(x_ix,k_ik)，为一个N_x行、n₀n₁n₂n₃列的矩阵；其中，F(x_ix,k_ik)为不等间隔离散逆傅里叶变换矩阵，ix为矩阵的行索引,ix的取值范围是0,1,2,…,N_x-1，ik为矩阵的列索引，ik的取值范围是0,1,2,…,n₀n₁n₂n₃-1，其中，为虚数单位，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η的取值范围是0,1,2,3，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积，ik与 保持关系式

为一个包含N_x个元素的列向量；

函数conj()表示对输入向量的每一个元素求复共轭，并输出一个向量；运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和，并输出一个标量；运算符^H表示求矩阵的共轭转置；

函数y＝CG(b,w)表示如下预条件共轭梯度算法，①设定常量，包括：迭代终止相对容差tol＝10^-4，最大迭代次数maxit＝30，A＝F^HF；②计算迭代初始变量，包括：r₀＝b，ρ_-1＝1；③计算迭代终止绝对容差对于i_cg＝{0,1,2,…,maxit-1}，执行如下处理操作：

如果或者i_cg＝maxit，跳出循环，返回否则继续循环:

其中，为标量，表示标量分别与向量的每一个元素相乘，并返回一个向量，表示标量分别与向量的每一个元素相乘，并返回一个向量；

在一个实例中，所述函数y＝CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q＝Ap，其中A＝F^HF具有多级分块Toeplitz矩阵结构，采用如下快速算法计算q＝Ap：

①设N₀＝n₁n₂n₃、N₁＝n₂n₃、N₂＝n₃、N₃＝1，m_η＝2n_η-1,η＝0,1,2,3，按如下方式生成长度为m_η的向量和其中η＝0,1,2,3：

②设向量是一个四维数组的向量化表示，i_a为向量a的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，此处表示i_a是的函数，以下均采用了类似表达，按如下方式对a的元素赋值，初始化变量

对于执行如下循环(1)至(12)

(1)

(2)

(3)对于执行如下循环(4)至(12)

(4)

(5)

(6)对于执行如下循环(7)至(12)

(7)

(8)

(9)对于执行如下循环(10)至(12)

(10)

(11)

(12)

其中，表示向量的第个元素，η＝0,1,2,3，表示向量a的第个元素，表示矩阵A的第行、列的元素；

③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量

其中，是一个四维数组的向量化表示，为向量的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，函数FFT₄()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组；

④设向量是一个四维数组的向量化表示，i_p为向量p的索引，为向量p对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

设向量是一个四维数组的向量化表示，i_p′为向量p′的索引，为向量p′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；按照如下方式对向量p′的各个元素赋值：

⑤计算其中，函数FFT₄()和IFFT₄()分别表示对四维数组做四维快速傅里叶变换和逆变换，并返回一个四维数组；运算符·*表示两个四维数组的对应元素相乘，并得到一个四维数组；

其中，向量是一个四维数组的向量化表示，i_q′为向量q′的索引，为向量q′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；

⑥设向量是一个四维数组的向量化表示，i_q为向量q的索引，为向量q对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，按照如下方式对向量q的各个元素赋值,其中0≤i_η≤n_η-1，η＝0,1,2,3。

在一个实例中，还包括：对y_iω执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组y的元素位置进行调整，并返回一个四维数组其中：

是一个四维数组的向量化表示，i_y为向量y的索引，为向量y对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

其中c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，函数mod(v₁,v₂)表示v₁对v₂求0～v₂-1之间的余数，函数floor()表示向下取整。

在一个实例中，还包括：对执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组做四维快速傅里叶逆变换，并返回一个四维数组s，其中，

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，i_s为向量s的索引，为向量s对应的四维数组索引，i_s的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

在一个实例中，所述五维插值数据变换模块具体用于：

利用s_iω生成插值后的频率域五维数据其中，为插值后的空间采样坐标，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，iω的取值上限由iω_max变为N_t-1，以满足后续进行快速傅里叶逆变换的需求，

其中，i_η的取值范围是0,1,2,…,n_η-1，vmin_η为插值后各个空间维坐标的最小值，为插值后各个空间维的采样间隔，η＝0,1,2,3，ix与i₀、i₁、i₂、i₃保持关系式ix＝((i₀·n₁+i₁)·n₂+i₂)·n₃+i₃；生成方式如下，用向量代表的频率坐标索引为iω的一部分数据，

则

其中，函数conj()表示对输入向量的每一个元素求复共轭，z表示由n₀n₁n₂n₃个0元素构成的向量；

利用快速傅里叶逆变换将变换的时间域，得到插值后的五维数据

其中，与满足如下关系式，其中，t_it＝it·Δt，it的取值范围是0,1,2,…,N_t-1，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，N_t为采样点个数；Δt为时间采样间隔。

下面再以实例来进行说明，以便于理解如何实施本发明。

本发明实施例提供不规则地震数据的五维插值处理方法主要包括如下步骤：

(1)、对采集的地震数据进行前期处理操作，包括去噪、静校正、线性动校正等；

(2)、抽取想要进行五维插值的地震数据图2是本发明实施例提供的插值前的地震数据观测***图。

步骤(2)所述的t_it为时间采样坐标，下标it＝{0,1,2,…,N_t-1}，即时间方向有N_t个采样点，设时间采样间隔为Δt,则t_it＝it·Δt。本例中，N_t＝1200，Δt＝0.002秒；

步骤(2)所述的为空间采样坐标，下标ix＝{0,1,2,…,N_x-1}，即空间方向有N_x个采样点，是四维空间中的一个点，其中是一种一般表示形式，根据本例需要，分别代表中心点x坐标、中心点y坐标，炮检距，方位角。本例中，N_x＝7151。

(3)、对中的每一维坐标进行归一化变换，变换结果记为x_ix＝((x₀)_ix,(x₁)_ix,(x₂)_ix,(x₃)_ix)，变换后的地震数据表示为

步骤(3)所述的归一化变换方法为其中，η＝{0,1,2,3}表示不同的维，对于任意一维坐标取其最大值和最小值，分别记为vmax_η和vmin_η，具体地和n_η是根据实际需要而设定的插值后的坐标点数，插值后各个空间维的采样间隔为本例中，vmin₀＝253，vmax₀＝263，vmin₁＝660，vmax₁＝670，vmin₂＝-160，vmax₂＝160，vmin₃＝0，vmax₃＝170，n₀＝11，n₁＝11，n₂＝65，n₃＝18，

(4)、利用快速傅里叶变换将地震数据d(x_ix,t_it)变换到频率域d(x_ix,ω_iω)。

步骤(4)所述的d(x_ix,ω_iω)与d(x_ix,t_it)满足如下关系式，其中，t_it＝it·Δt，ω_iω＝iω·Δω，iω＝{0,1,2,…,N_t-1}。这里需要保证地震数据d(x_ix,t_it)满足采样定理，设地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内，其中则需要满足条件iω_max≤floor(N_t/2)。本例中，Δω≈2.618，iω_max＝175。

(5)、生成不等间隔离散逆傅里叶变换矩阵F(x_ix,k_ik)，其中，ix＝{0,1,2,…,N_x-1}为矩阵的行索引，ik＝{0,1,2,…,n₀n₁n₂n₃-1}为矩阵的列索引。本例中，n₀n₁n₂n₃＝141570。

步骤(5)所述其中，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，η＝0,1,2,3，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积。并且ik与保持关系式其中，n₀、n₁、n₂、n₃的定义见步骤(3)。本例中，c₀＝5，c₁＝5，c₂＝32，c₃＝8。

(6)、对每一个频率切片iω＝{0,1,2,…,iω_max}，执行如下循环，并存储y_iω。

步骤(6)所述的iω_max是根据实际需要选取的频率上限索引，选取原则是使得地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内，其中

步骤(6)所述的w₀为一个包含n₀n₁n₂n₃个元素的列向量，且w₀的每一个元素都等于其中，w₀是w_iω中下角标iω＝0时的情况。

步骤(6)所述的F＝F(x_ix,k_ik)，为一个N_x行、n₀n₁n₂n₃列的矩阵。

步骤(6)所述的为一个包含N_x个元素的列向量。

步骤(6)所述的函数conj()表示对输入向量的每一个元素求复共轭，并输出一个向量；运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和，并输出一个标量；运算符^H表示求矩阵的共轭转置，下同。

步骤(6)所述的函数y＝CG(b,w)表示如下预条件共轭梯度算法，①设定常量。包括迭代终止相对容差tol＝10^-4，最大迭代次数maxit＝30，A＝F^HF；②计算迭代初始变量。包括r₀＝b，ρ_-1＝1；③计算迭代终止绝对容差对于i_cg＝{0,1,2,…,maxit}，执行如下循环

如果或者i_cg＝maxit，跳出循环，返回否则继续循环:

其中，为标量，表示标量分别与向量的每一个元素相乘，并返回一个向量，表示标量分别与向量的每一个元素相乘，并返回一个向量；

步骤(6)所述的函数y＝CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q＝Ap，其中A＝F^HF具有多级分块Toeplitz矩阵结构，可以采用如下快速算法计算q＝Ap。

①设N₀＝n₁n₂n₃、N₁＝n₂n₃、N₂＝n₃、N₃＝1，m_η＝2n_η-1,η＝0,1,2,3，按如下方式生成长度为m_η的向量和其中η＝0,1,2,3：

②设向量是一个四维数组的向量化表示，i_a为向量a的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，此处表示i_a是的函数，以下均采用了类似表达。按如下方式对a的元素赋值。初始化变量

对于执行如下循环(1)至(12)

(1)

(2)

(3)对于执行如下循环(4)至(12)

(4)

(5)

(6)对于执行如下循环(7)至(12)

(7)

(8)

(9)对于执行如下循环(10)至(12)

(10)

(11)

(12)

其中，表示向量的第个元素，η＝0,1,2,3，表示向量a的第个元素，表示矩阵A的第行、列的元素。

③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量其中，向量其中，是一个四维数组的向量化表示，为向量的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3。函数FFT₄()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组。

④设向量是一个四维数组的向量化表示，i_p为向量p的索引，为向量p对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

设向量是一个四维数组的向量化表示，i_p′为向量p′的索引，为向量p′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3。按照如下方式对向量p′的各个元素赋值

⑤计算其中，函数FFT₄()和IFFT₄()分别表示对四维数组做四维快速傅里叶变换和逆变换，并返回一个四维数组；运算符·*表示两个四维数组的对应元素相乘，并得到一个四维数组；其中，向量是一个四维数组的向量化表示，i_q′为向量q′的索引，为向量q′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3。

⑥设向量是一个四维数组的向量化表示，i_q为向量q的索引，为向量q对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。按照如下方式对向量q的各个元素赋值,其中0≤i_η≤n_η-1，η＝0,1,2,3。

⑦整个方法流程中，以上第①～③步只需要计算一次，运算量基本可以忽略，因此，如第⑤步所述，矩阵与向量乘积运算q＝Ap的计算量主要包括一次四维快速傅里叶变换和一次四维快速傅里叶逆变换。

(7)、对y_iω执行函数并存储其中iω的取值范围是0,1,2,…,iω_max。

步骤(7)所述的函数表示对四维数组y的元素位置进行调整，并返回一个四维数组其中，

是一个四维数组的向量化表示，i_y为向量y的索引，为向量y对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

具体调整方法如下，其中c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3,，函数mod(v₁,v₂)表示v₁对v₂求0～v₂-1之间的余数，函数floor()表示向下取整。

(8)、对执行函数并存储s_iω，其中iω的取值范围是0,1,2,…,iω_max。

步骤(8)所述的函数表示对四维数组做四维快速傅里叶逆变换，并返回一个四维数组s。其中，

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，i_s为向量s的索引，为向量s对应的四维数组索引，i_s的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

(9)、利用s_iω生成插值后的频率域五维数据其中，为插值后的空间采样坐标，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，iω的取值上限由iω_max变为N_t-1，以满足后续进行快速傅里叶逆变换的需求。

步骤(9)所述的

其中，i_η的取值范围是0,1,2,…,n_η-1，vmin_η为插值后各个空间维坐标的最小值，定义见步骤(3)。并且ix与i₀、i₁、i₂、i₃保持关系式ix＝((i₀·n₁+i₁)·n₂+i₂)·n₃+i₃。

步骤(9)所述的生成方式如下，为了方便表述，用向量代表的一部分数据，具体地，则

其中，函数conj()表示对输入向量的每一个元素求复共轭。z表示由n₀n₁n₂n₃个0元素构成的向量。

(10)、利用快速傅里叶逆变换将变换的时间域，得到插值后的五维数据

步骤(10)所述的与满足如下关系式，其中，t_it＝it·Δt，it的取值范围是0,1,2,…,N_t-1，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，N_t为采样点个数；Δt为时间采样间隔。

图2是本发明实施例中插值前的地震数据观测***示意图；图3是本发明实施例中插值前的同一个CMP面元的地震数据观测***示意图；图4是本发明实施例中与图3对应的插值后的的同一个CMP面元的地震数据观测***示意图；图5是本发明实施例中采用不同快速算法的运算时间对比图表示意图；图6是本发明实施例中与图3对应的插值前的CMP道集地震数据示意图；图7是本发明实施例中与图4对应的插值后的CMP道集地震数据。

与现有技术相比较，通过图2至图7所示和上述实施例的介绍可知，本发明实施例提供的不规则采集地震数据的五维插值快速算法，核心是提出了一种新的迭代方案,大大提高了计算效率。整个计算过程仅涉及少量NFFT运算，绝大部分的矩阵与向量乘积运算都可以利用FFT实现，大大提高了计算效率，具有重要的实际应用价值。

本领域内的技术人员应明白，本发明的实施例可提供为方法、***、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种不规则地震数据的五维插值处理方法，其特征在于，包括：

将不规则地震数据变换到频率域和空间域；

对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对变换到频率域和空间域的地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域；

对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对变换到频率域和空间域的地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据，包括：

对每一个频率切片iω＝{0,1,2,…,iω_max}，按照如下公式进行处理：

b_iω＝F^Hd_iω；

y_iω＝CG(b_iω,w_iω)；

u_iω＝conj(y_iω)·*y_iω；

w_iω+1＝(sqrt(u_iω/sum(u_iω))+w_iω)/2；

其中，iω_max是根据实际需要选取的频率上限索引，选取原则是使得地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内；

w₀为一个包含n₀n₁n₂n₃个元素的列向量，且w₀的每一个元素都等于其中，w₀是w_iω中下角标iω＝0时的情况；

F＝F(x_ix,k_ik)，为一个N_x行、n₀n₁n₂n₃列的矩阵；其中，F(x_ix,k_ik)为不等间隔离散逆傅里叶变换矩阵，ix为矩阵的行索引,ix的取值范围是0,1,2,…,N_x-1，ik为矩阵的列索引，ik的取值范围是0,1,2,…,n₀n₁n₂n₃-1，其中，为虚数单位，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η的取值范围是0,1,2,3，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积；

为一个包含N_x个元素的列向量；

函数conj()表示对输入向量的每一个元素求复共轭，并输出一个向量；运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和，并输出一个标量；运算符H表示求矩阵的共轭转置；

函数y＝CG(b,w)表示如下预条件共轭梯度算法，①设定常量，包括：迭代终止相对容差tol＝10^-4，最大迭代次数maxit＝30，A＝F^HF；②计算迭代初始变量，包括：r₀＝b，ρ_-1＝1；③计算迭代终止绝对容差对于i_cg＝{0,1,2,…,maxit-1}，执行如下处理操作：

如果或者i_cg＝maxit，跳出循环，返回否则继续循环:

其中，为标量，表示标量分别与向量的每一个元素相乘，并返回一个向量，表示标量分别与向量的每一个元素相乘，并返回一个向量；

x_ix为变换后的地震数据的空间采样坐标，t_it为地震数据的时间采样坐标；

d(x_ix,ω_iω)为d(x_ix,t_it)变换到频率域后的频率为ω_iω的地震数据,ω_iω的下标iω为频率坐标的离散索引，d_iω为频率切片向量；b_iω为频率切片向量d_iω的不等间隔傅里叶变换；y_iω为预条件共轭梯度算法的输出；u_iω为y_iω的二范数距离的平方；w_iω为预条件共轭梯度算法的输入权值；

所述函数y＝CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q＝Ap，其中A＝F^HF具有多级分块Toeplitz矩阵结构，采用如下快速算法计算q＝Ap：

①设N₀＝n₁n₂n₃、N₁＝n₂n₃、N₂＝n₃、N₃＝1，m_η＝2n_η-1,η＝0,1,2,3，按如下方式生成长度为m_η的向量和其中η＝0,1,2,3：

②设向量是一个四维数组的向量化表示，i_a为向量a的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，此处表示i_a是的函数，以下均采用了类似表达，按如下方式对a的元素赋值，初始化变量

对于执行如下循环(1)至(12)

(3)对于执行如下循环(4)至(12)

(6)对于执行如下循环(7)至(12)

(9)对于执行如下循环(10)至(12)

其中，表示向量的第个元素，η＝0,1,2,3，表示向量a的第个元素，表示矩阵A的第行、列的元素；

③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量

其中，向量是一个四维数组的向量化表示，为向量的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，函数FFT₄()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组；

④设向量是一个四维数组的向量化表示，i_p为向量p的索引，为向量p对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

设向量是一个四维数组的向量化表示，i_p′为向量p′的索引，为向量p′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；按照如下方式对向量p′的各个元素赋值：

⑤计算其中，函数FFT₄()和IFFT₄()分别表示对四维数组做四维快速傅里叶变换和逆变换，并返回一个四维数组；运算符·*表示两个四维数组的对应元素相乘，并得到一个四维数组；

向量是一个四维数组的向量化表示，i_q′为向量q′的索引，为向量q′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；

⑥设向量是一个四维数组的向量化表示，i_q为向量q的索引，为向量q对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，按照如下方式对向量q的各个元素赋值,其中i_η的取值范围是0,1,2,…,n_η-1。

2.如权利要求1所述的不规则地震数据的五维插值处理方法，其特征在于，所述将不规则地震数据变换到频率域和空间域之前，包括：

对不规则地震数据进行去噪、静校正和线性动校正处理；

从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据；

所述将不规则地震数据变换到频率域和空间域，包括：将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域。

3.如权利要求1所述的不规则地震数据的五维插值处理方法，其特征在于，还包括：对y_iω执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组y的元素位置进行调整，并返回一个四维数组其中：

是一个四维数组的向量化表示，i_y为向量y的索引，为向量y对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

其中c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，函数mod(v₁,v₂)表示v₁对v₂求0～v₂-1之间的余数，函数floor()表示向下取整。

4.如权利要求3所述的不规则地震数据的五维插值处理方法，其特征在于，还包括：对执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组做四维快速傅里叶逆变换，并返回一个四维数组s，其中，

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，i_s为向量s的索引，为向量s对应的四维数组索引，i_s的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

5.如权利要求4所述的不规则地震数据的五维插值处理方法，其特征在于，将插值后的频率波数域数据变换到频率和空间域，再变换到时间和空间域，完成对不规则地震数据的五维插值处理，包括：

利用s_iω生成插值后的频率域五维数据其中，为插值后的空间采样坐标，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，iω的取值上限由iω_max变为N_t-1，以满足后续进行快速傅里叶逆变换的需求；

其中，i_η的取值范围是0,1,2,…,n_η-1，vmin_η为插值后各个空间维坐标的最小值，为插值后各个空间维的采样间隔，η＝0,1,2,3；生成方式如下，用向量代表的频率坐标索引为iω的一部分数据，则

其中，函数conj()表示对输入向量的每一个元素求复共轭，z表示由n₀n₁n₂n₃个0元素构成的向量；

利用快速傅里叶逆变换将变换的时间域，得到插值后的五维数据

其中，与满足如下关系式，其中，t_it＝it·Δt，it的取值范围是0,1,2,…,N_t-1，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，N_t为采样点个数；Δt为时间采样间隔。

6.一种不规则地震数据的五维插值处理装置，其特征在于，包括：

地震数据变换模块，用于将不规则地震数据变换到频率域和空间域；

五维插值处理模块，用于对于每一个变换到频率域和空间域的地震数据的频率切片，按照快速傅里叶变换算法，对变换到频率域和空间域的地震数据的矩阵和向量进行乘积运算，求取五维插值后的频率波数域数据；

五维插值数据变换模块，用于将五维插值后的频率波数域数据变换到频率域和空间域，再变换到时间域和空间域；

所述五维插值处理模块具体用于：

对每一个频率切片iω＝{0,1,2,…,iω_max}，按照如下公式进行处理：

b_iω＝F^Hd_iω；

y_iω＝CG(b_iω,w_iω)；

u_iω＝conj(y_iω)·*y_iω；

w_iω+1＝(sqrt(u_iω/sum(u_iω))+w_iω)/2；

其中，iω_max是根据实际需要选取的频率上限索引，选取原则是使得地震数据d(x_ix,t_it)的有效信号能量刚好处于的频带内，其中

w₀为一个包含n₀n₁n₂n₃个元素的列向量，且w₀的每一个元素都等于其中，w₀是w_iω中下角标iω＝0时的情况；

F＝F(x_ix,k_ik)，为一个N_x行、n₀n₁n₂n₃列的矩阵；其中，F(x_ix,k_ik)为不等间隔离散逆傅里叶变换矩阵，ix为矩阵的行索引,ix的取值范围是0,1,2,…,N_x-1，ik为矩阵的列索引，ik的取值范围是0,1,2,…,n₀n₁n₂n₃-1，其中，为虚数单位，是四维空间中的一个点，c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η的取值范围是0,1,2,3，函数floor()表示向下取整，k_ikx_ix表示四维向量k_ik与x_ix的内积；

为一个包含N_x个元素的列向量；

函数conj()表示对输入向量的每一个元素求复共轭，并输出一个向量；运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和，并输出一个标量；运算符H表示求矩阵的共轭转置；

函数y＝CG(b,w)表示如下预条件共轭梯度算法，①设定常量，包括：迭代终止相对容差tol＝10^-4，最大迭代次数maxit＝30，A＝F^HF；②计算迭代初始变量，包括：r₀＝b，ρ_-1＝1；③计算迭代终止绝对容差对于i_cg＝{0,1,2,…,maxit-1}，执行如下处理操作：

如果或者i_cg＝maxit，跳出循环，返回否则继续循环:

其中，为标量，表示标量分别与向量的每一个元素相乘，并返回一个向量，表示标量分别与向量的每一个元素相乘，并返回一个向量；

x_ix为变换后的地震数据的空间采样坐标，t_it为地震数据的时间采样坐标；

d(x_ix,ω_iω)为d(x_ix,t_it)变换到频率域后的频率为ω_iω的地震数据,ω_iω的下标iω为频率坐标的离散索引，d_iω为频率切片向量；b_iω为频率切片向量d_iω的不等间隔傅里叶变换；y_iω为预条件共轭梯度算法的输出；u_iω为y_iω的二范数距离的平方；w_iω为预条件共轭梯度算法的输入权值；

所述函数y＝CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q＝Ap，其中A＝F^HF具有多级分块Toeplitz矩阵结构，采用如下快速算法计算q＝Ap：

①设N₀＝n₁n₂n₃、N₁＝n₂n₃、N₂＝n₃、N₃＝1，m_η＝2n_η-1,η＝0,1,2,3，按如下方式生成长度为m_η的向量和其中η＝0,1,2,3：

②设向量是一个四维数组的向量化表示，i_a为向量a的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，此处表示i_a是的函数，以下均采用了类似表达，按如下方式对a的元素赋值，初始化变量

对于执行如下循环(1)至(12)

(3)对于执行如下循环(4)至(12)

(6)对于执行如下循环(7)至(12)

(9)对于执行如下循环(10)至(12)

其中，表示向量的第个元素，η＝0,1,2,3，表示向量a的第个元素，表示矩阵A的第行、列的元素；

③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量

其中，是一个四维数组的向量化表示，为向量的索引，为向量a对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3，函数FFT₄()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组；

④设向量是一个四维数组的向量化表示，i_p为向量p的索引，为向量p对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

设向量是一个四维数组的向量化表示，i_p′为向量p′的索引，为向量p′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；按照如下方式对向量p′的各个元素赋值：

⑤计算其中，函数FFT₄()和IFFT₄()分别表示对四维数组做四维快速傅里叶变换和逆变换，并返回一个四维数组；运算符·*表示两个四维数组的对应元素相乘，并得到一个四维数组；

向量是一个四维数组的向量化表示，i_q′为向量q′的索引，为向量q′对应的四维数组索引，的取值范围是0,1,2,…,m_η-1，η＝0,1,2,3；

⑥设向量是一个四维数组的向量化表示，i_q为向量q的索引，为向量q对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，按照如下方式对向量q的各个元素赋值,其中i_η的取值范围是0,1,2,…,n_η-1。

7.如权利要求6所述的不规则地震数据的五维插值处理装置，其特征在于，还包括：

地震数据预处理模块，用于对不规则地震数据进行去噪、静校正和线性动校正处理；

地震数据抽取模块，用于从处理后的不规则地震数据中抽取待五维插值处理的不规则地震数据；

所述地震数据变换模块具体用于：将抽取的待五维插值处理的不规则地震数据变换到频率域和空间域。

8.如权利要求6所述的不规则地震数据的五维插值处理装置，其特征在于，还包括：对y_iω执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组y的元素位置进行调整，并返回一个四维数组其中：

是一个四维数组的向量化表示，i_y为向量y的索引，为向量y对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

其中c_η＝floor((n_η-1)/2)，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3，函数mod(v₁,v₂)表示v₁对v₂求0～v₂-1之间的余数，函数floor()表示向下取整。

9.如权利要求8所述的不规则地震数据的五维插值处理装置，其特征在于，还包括：

对执行函数其中，iω的取值范围是0,1,2,…,iω_max；

函数表示对四维数组做四维快速傅里叶逆变换，并返回一个四维数组s，其中，

是一个四维数组的向量化表示，为向量的索引，为向量对应的四维数组索引，的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3；

是一个四维数组的向量化表示，i_s为向量s的索引，为向量s对应的四维数组索引，i_s的取值范围是0,1,2,…,n_η-1，η＝0,1,2,3。

10.如权利要求9所述的不规则地震数据的五维插值处理装置，其特征在于，所述五维插值数据变换模块具体用于：

利用s_iω生成插值后的频率域五维数据其中，为插值后的空间采样坐标，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，iω的取值上限由iω_max变为N_t-1，以满足后续进行快速傅里叶逆变换的需求；

其中，i_η的取值范围是0,1,2,…,n_η-1，vmin_η为插值后各个空间维坐标的最小值，为插值后各个空间维的采样间隔，η＝0,1,2,3；生成方式如下，用向量代表的频率坐标索引为iω的一部分数据，

则

其中，函数conj()表示对输入向量的每一个元素求复共轭，z表示由n₀n₁n₂n₃个0元素构成的向量；

利用快速傅里叶逆变换将变换的时间域，得到插值后的五维数据

其中，与满足如下关系式，其中，t_it＝it·Δt，it的取值范围是0,1,2,…,N_t-1，的取值范围是0,1,2,…,n₀n₁n₂n₃-1，N_t为采样点个数；Δt为时间采样间隔。