CN103105623B - 一种地震勘探中的数据波形处理方法 - Google Patents

Abstract

本发明公开了一种地震勘探中的数据波形处理方法，其具体是基于波动方程的三维自由表面多次波预测方法，应用于地震勘探资料处理中的三维表面多次波预测，该方法利用GPU(图形处理器)加速全三维表面多次波预测算法，即通过GPU和CPU协同并行计算，将计算密集的运算转移到GPU上进行计算，可获得较高的计算效率；可处理来自复杂地下介质的地震资料。该方法考虑了一次反射点与炮点和接收点不在同一条线上的空间效应、优于常规的二维算法，无需对地下介质进行简单近似，因而基于波动方程的全三维的表面多次波预测算法更符合地下介质的真实情况，可使得地震资料多次波的振幅和相位得到较为准确的预测。

Description

一种地震勘探中的数据波形处理方法

技术领域

本发明涉及地球物理勘探领域，具体涉及地球物理勘探数据波形处理方法，特别涉及为基于图形处理器利用反馈迭代法计算三维表面多次波预测方法，该方法应用于地震勘探资料处理中的三维表面多次波预测。

背景技术

丰富的油气资源蕴藏于海域，但在海洋勘探地震资料中常包含有丰富的、强能量的表面多次波。多次波的存在会影响地震资料的偏移成像效果，引起对地震资料解释乃至地下构造的错误认识，并进一步影响探井井位的部署以及勘探的成功率。因此，为了提高海洋地震勘探资料的信噪比，对地震资料中的多次波通常予以压制。基于波动方程的多次波压制方法通常包含两个步骤：即多次波预测和自适应相减。基于波动方程的反馈迭代法是预测表面多次波，尤其是复杂地下介质表面多次波的有力工具；该方法为数据驱动的方法，也更适合于当前复杂地下介质勘探的迫切需要。最初，Anstey and Newman(1966)利用单道的自相关压制表面多次波；Kennett(1979)提出了一维空间表面多次波的正演模拟方法和反演解决方案，但是，此算法需要对数据采集方式和地下介质的性质作许多假设才能应用于实际数据，所以实际应用效果并不好；荷兰Delft大学以Berkhout和Verschuur为首的DELPHI小组在此方面作出了卓越的贡献，先后提出了SRME(surface-related multiple elimination)方法、CFP(common-focus-point)方法、以及逆数据域(inverse data space)方法。Berkhout(1982)成功地提出了数据矩阵和反馈模型的概念，初步奠定了反馈迭代法压制表面多次波的基础，并考虑了数据的采集性质。Verschuur and Berkhout(1992)正式提出了利用级数展开法压制表面多次波。Berkhout(1997)和Verschuur(1997)将级数展开法发展为迭代法压制多次波，并将其成功应用于复杂地下构造的多次波压制，其良好的叠前处理优势，受到了许多地球物理工作者的关注与研究(Dragoset and Jericevic，1998；Wapenaar，1990；Borselen，2005；Pica，2005；牛滨华，2002)，并被广泛应用于工业界。Kelamis and Verschuur(2000)将SRME方法应用于陆上地震数据。Wang(2004，2007)直接通过矩阵求逆方案预测近偏移距缺失地震数据的多次波，巧妙回避常规的近偏移距外推问题，但该方法仍包含预测和相减过程。

上述的二维方法可有效应用于一条测线的表面多次波预测和压制，但由于未考虑地震波在地下介质中传播的空间效应，因此无法有效应用于三维地震勘探资料。

发明内容

针对传统的二维地震资料表面多次波预测方法的不足之处，本发明提供了一种地震勘探中的数据波形处理方法，该方法为基于图形处理器计算三维表面多次波预测方法。

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

一种地震勘探中的数据波形处理方法，其具体为基于图形处理器利用反馈迭代法计算三维表面多次波预测方法，其采用以下步骤：

1)使用海上专用的电缆和检波器，在观测船航行中连续进行地震波的激发和接收，采集含有自由表面多次波的三维地震数据；

2)对采集到的三维地震数据，采用Radon变换的方法实现地震数据规则化，使得规则化后地震数据的炮间距与道间距相等，进而得到规则化后的时空域数据，为反馈迭代法计算准备地震数据；

3)利用傅里叶变换将规则化后的时空域数据变换到频率空间域；

4)循环读取各个频率分量上的数据，在每个频率分量上形成数据矩阵：

假定在x方向放置N_x个检波器，在y方向放置N_y个检波器，则每炮数据有N_xN_y个检波器接收，即每一炮会产生N_xN_y道记录，假定在此矩形区域的N_xN_y个位置都放炮，每个网格点做一次震源位置，每炮记录在N_xN_y个检波器上接收到的信号被记录下来，结果会得到N_xN_yN_xN_y道数据，就是(N_xN_y)²道数据；将所有地震数据道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集信息；在每个震源位置都做上述构建数据矩阵的重复，结果得到大小为N_xN_y的方阵；所形成的数据矩阵中，矩阵的每列为3D共炮点道集，每行为3D共接收点道集；

5)在每个频率分量上对数据矩阵做乘法运算：

在每个频率分量上，应用步骤4)形成的数据矩阵完成矩阵的自乘运算，使用矩阵乘法表述来表示x和y方向二维空间褶积；

6)将预测的多次波数据进行反傅里叶变换到时空域；

上述步骤5)在图形处理器中进行。

其中，反馈迭代法预测多次波的过程是地震数据与地震数据本身在时空域的褶积。

进一步地，针对三维地震数据模型，采用反馈迭代法预测多次波(SRME)，其预测自由表面多次波的方程可表示为：

$M_0^{\left(i\right)}(x_r,y_r,x_s,y_s,f)=\underset{y_k}{\mathrm{\Σ}}\underset{x_k}{\mathrm{\Σ}}{P}_0^{(i-1)}(x_r,y_r,x_k,y_k,f)P(x_k,y_k,x_s,y_s,f)$

式中：是第i次迭代预测的多次波，

是i-1次迭代后估计的有效波，P是含多次波的总波场。

优选地，步骤(2)中采用Radon变换的方法实现地震数据规则化具体步骤为：

对任一CMP道集，Radon变换的最小二乘正变换公式为

(L^HS^TSL+λI)M＝L^HS^TSD

式中S为对角阵，I单位阵，M和D分别是模型空间和数据空间的向量，该式实现了地震数据从数据空间到模型空间的变换；

上式中，左端矩阵(L^HS^TSL+λI)元素可表示为

$a_{\mathrm{ij}}=\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{s}_l^2\exp \left[\left({\mathrm{i\ωx}}_l^2\right)(q_k-q_l)\right]+{\mathrm{\λI}}_{\mathrm{jk}}$

式中n为CMP道集中的地震道数，x_l为各道的偏移距，s_l为各道的权系数，I_jk为单位矩阵的各元素；

与Radon变换的最小二乘正变换公式相应的Radon域采样公式为：

$\mathrm{\Δq}\≤\frac{2\mathrm{\π}}{\mathrm{\ω}(x_{\max }^2-x_{\min }^2+4\mathrm{\Δ}{x}_a^2)}$

$q_{\max }\≤\frac{2\mathrm{\π}}{\mathrm{\ω\Δ}{x}_a^2}$

其中，x_max是最大偏移距，x_min是最小偏移距，是原始地震数据做x平方拉伸后沿x²轴的最大间隔；均匀采样情况， ${\mathrm{\Δx}}_a^2=x_{\max }^2-{(x_{\max }-\mathrm{\Δx})}^2={2x}_{\max }\mathrm{\Δx}-{\left(\mathrm{\Δx}\right)}^2.$

对每个频率成分完成上述步骤后，再按下式进行Radon反变换，即：

D＝LM

可实现地震数据规则化。

另外地，步骤(4)中在每个频率分量上形成数据矩阵具有如图9所示的表达形式，其中N_x和N_y分别是三维数据体的主测线和联络测线数。

假定在此矩形区域的N_xN_y个位置都放炮，每个网格点做一次震源位置，每炮记录在N_xN_y个检波器上接收到的信号被记录下来，结果会得到(N_xN_y)²道数据。所有地震数据道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集信息，在每个震源位置都做这样的重复，结果得到大小为N_xN_y的方阵。用这样的矩阵表示，矩阵的每列为3D共炮点道集，每行为3D共接收点道集。

优选地，步骤(5)中在每个频率分量上对数据矩阵做乘法，利用统一计算设备架构平台CUDA的函数matrix_mul实现。

进一步地，应用CUDA编程语言实现GPU通用计算，其具体的步骤和方法是，在首个低频频率分量上，将数据矩阵由CPU传送到GPU上，在GPU上完成数据矩阵的乘法运算，再将数据矩阵的乘积结果传送回CPU上，完成此频率分量的多次波预测，对其他的频率分量，依次执行此操作，可将所有频率分量上的矩阵乘法执行完毕，从而完成表面多次波预测。

有益效果

本发明所述的全三维地震资料表面多次波预测方法，考虑了地震射线在地下介质中传播的可能射线路径，符合地下介质对地震波反射的真实情况，预测的表面多次波更为真实、可靠。与滤波法比较，基于波动方程的表面多次波预测方法是数据直接生成多次波的方法，利用空间褶积预测多次波，而不经过速度模型等参数模型的估计，可处理复杂地下介质勘探的三维地震资料，由于其适合横向变速，并可作为预测和消除长程层间多次波的基础，得到广泛的研究。为了在更广泛应用时便于理解，SRME就特指利用空间褶积进行多次波预测的方法。而基于信号处理的滤波法仅对简单介质或中等复杂介质的表面多次波处理有效，无法适应“低、深、难、隐”的复杂勘探目标的地震勘探资料。

基于波动方程的表面多次波预测方法，其本质思想来源于数据矩阵和反馈模型，该方法为完全数据驱动的方法，不需要地下介质的先验信息，可应用于复杂地下构造的多次波预测。但是预测过程中，尤其在三维表面多次波预测中，涉及的大量空间褶积运算，使得预测效率较低。图形处理器(GPU)在科学计算方面呈现的高速度和高性能的特点，以及其较高的并行计算效率，为提高SRME方法表面多次波预测效率提供了新的契机，GPU本身具有可编程性，具有强大的浮点运算能力，可同时执行成百上千个进程，其并行度远远超过CPU串行执行程序的能力。因此GPU更加适用于高并行性和高密集度的数值计算，GPU以CUDA(统一计算设备架构)为编程平台，CUDA编程语言是针对通用计算GPU的C语言环境，应用CUDA编程可方便地实现GPU通用计算。利用GPU和CPU协同并行预测表面多次波，可极大地提高计算效率，对基于波动方程的三维表面多次波预测技术的应用有重要意义。

附图说明

附图1为复杂介质中一阶表面多次波的射线路径，其中S为震源位置，R为检波点位置，经过地表反射的反射点、震源以及检波点不在同一条测线上。

附图2为三维地震采集观测系统，若灰色的圈代表检波点位置，那么在每个检波点位置都要放置震源，图中仅用三个震源示意，所有的网格点都是接收点，并且每个网格点做一次震源位置。

附图3为三维水平层状介质表面多次波正演数据的速度模型。

附图4为本发明基于图形处理器计算三维表面多次波预测方法的流程示意图。

附图5为本发明基于图形处理器分块计算数据矩阵乘法。

附图6为理论模拟的某炮数据56条测线的含一阶表面多次波的地震记录。

附图7为利用全三维算法进行多次波预测的结果。

附图8为三维含表面多次波数据的压制效果。

附图9为在每个频率分量上形成数据矩阵的表达形式。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明中，一种地震勘探中的数据波形处理方法具体为基于图形处理器利用反馈迭代法计算三维表面多次波预测方法，其采用以下步骤：

1)使用海上专用的电缆和检波器，在观测船航行中连续进行地震波的激发和接收，采集含有自由表面多次波的三维地震数据；

2)对采集到的三维地震数据，采用Radon变换的方法实现地震数据规则化，使得规则化后地震数据的炮间距与道间距相等，进而得到规则化后的时空域数据，为反馈迭代法计算准备地震数据；

3)利用傅里叶变换将规则化后的时空域数据变换到频率空间域；

4)循环读取各个频率分量上的数据，在每个频率分量上形成数据矩阵：

假定在x方向放置N_x个检波器，在y方向放置N_y个检波器，则每炮数据有N_xN_y个检波器接收，即每一炮会产生N_xN_y道记录，假定在此矩形区域的N_xN_y个位置都放炮，每个网格点做一次震源位置，每炮记录在N_xN_y个检波器上接收到的信号被记录下来，结果会得到N_xN_yN_xN_y道数据，就是(N_xN_y)²道数据；将所有地震数据道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集信息；在每个震源位置都做上述构建数据矩阵的重复，结果得到大小为N_xN_y的方阵；所形成的数据矩阵中，矩阵的每列为3D共炮点道集，每行为3D共接收点道集；

5)在每个频率分量上对数据矩阵做乘法运算：

在每个频率分量上，应用步骤4)形成的数据矩阵完成矩阵的自乘运算，使用矩阵乘法表述来表示x和y方向二维空间褶积；

6)将预测的多次波数据进行反傅里叶变换到时空域；

上述步骤5)在图形处理器中进行。

其中，反馈迭代法预测多次波的过程是地震数据与地震数据本身在时空域的褶积。

进一步地，针对三维地震数据模型，采用反馈迭代法预测多次波(SRME)，其预测自由表面多次波的方程可表示为：

$M_0^{\left(i\right)}(x_r,y_r,x_s,y_s,f)=\underset{y_k}{\mathrm{\Σ}}\underset{x_k}{\mathrm{\Σ}}{P}_0^{(i-1)}(x_r,y_r,x_k,y_k,f)P(x_k,y_k,x_s,y_s,f)$

式中：

是第i次迭代预测的多次波，

是i-1次迭代后估计的有效波，P是含多次波的总波场。

步骤(2)中采用Radon变换的方法实现地震数据规则化具体步骤为：

对任一CMP道集，Radon变换的最小二乘正变换公式为

(L^HS^TSL+λI)M＝L^HS^TSD

式中S为对角阵，I单位阵，M和D分别是模型空间和数据空间的向量，该式实现了地震数据从数据空间到模型空间的变换；

上式中，左端矩阵(L^HS^TSL+λI)元素可表示为

$a_{\mathrm{ij}}=\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{s}_l^2\exp \left[\left({\mathrm{i\ωx}}_l^2\right)(q_k-q_l)\right]+{\mathrm{\λI}}_{\mathrm{jk}}$

式中n为CMP道集中的地震道数，x_l为各道的偏移距，s_l为各道的权系数，I_jk为单位矩阵的各元素；

与Radon变换的最小二乘正变换公式相应的Radon域采样公式为：

$\mathrm{\Δq}\≤\frac{2\mathrm{\π}}{\mathrm{\ω}(x_{\max }^2-x_{\min }^2+4\mathrm{\Δ}{x}_a^2)}$

$q_{\max }\≤\frac{2\mathrm{\π}}{\mathrm{\ω\Δ}{x}_a^2}$

其中，x_max是最大偏移距，x_min是最小偏移距，

是原始地震数据做x平方拉伸后沿x²轴的最大间隔；均匀采样情况， ${\mathrm{\Δx}}_a^2=x_{\max }^2-{(x_{\max }-\mathrm{\Δx})}^2={2x}_{\max }\mathrm{\Δx}-{\left(\mathrm{\Δx}\right)}^2.$

对每个频率成分完成上述步骤后，再按下式进行Radon反变换，即：

D＝LM

可实现地震数据规则化。

另外地，步骤(4)中在每个频率分量上形成数据矩阵具有如图9的表达形式，其中N_x和N_y分别是三维数据体的主测线和联络测线数。

假定在此矩形区域的N_xN_y个位置都放炮，每个网格点做一次震源位置，每炮记录在N_xN_y个检波器上接收到的信号被记录下来，结果会得到(N_xN_y)²道数据。所有地震数据道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集信息，在每个震源位置都做这样的重复，结果得到大小为N_xN_y的方阵。用这样的矩阵表示，矩阵的每列为3D共炮点道集，每行为3D共接收点道集。

步骤(5)中在每个频率分量上对数据矩阵做乘法，利用统一计算设备架构平台CUDA的函数matrix_mul实现。

进一步地，应用CUDA编程语言实现GPU通用计算，其具体的步骤和方法是，在首个低频频率分量上，将数据矩阵由CPU传送到GPU上，在GPU上完成数据矩阵的乘法运算，再将数据矩阵的乘积结果传送回CPU上，完成此频率分量的多次波预测，对其他的频率分量，依次执行此操作，可将所有频率分量上的矩阵乘法执行完毕，从而完成表面多次波预测。

本发明基于以下客观事实：

1)表面多次波通常存在于海洋地震数据中，在海域激发震源，同时进行三维地震数据的采集，采集过程中，在经济限制允许的情况下，尽可能地增加拖缆的条数，以减小后续数据规则化的压力。

2)对实际采集到的三维地震数据，联络测线较稀疏，因此无法直接应用全三维表面多次波预测算法，需要结合实用的地震数据重建技术对原始地震数据进行规则化重建。本发明采用Radon变换的方法实现地震数据重建和规则化处理；

所述的地震数据重建和规则化处理是重建全波场数据，使炮间距与道间距相等；

3)基于波动方程的反馈迭代法预测表面多次波是通过在频域的每个频率分量上进行空间褶积运算完成的，为此，将规则化处理后的时空域数据经Fourier变换到频率空间域；

4)在每个频率分量上形成全三维基于波动方程的SRME方法所要求的数据矩阵，使得矩阵的每行为3D地震数据的共接收点道集，每列为3D地震数据的共炮点道集，此矩阵包含了当前频率分量上的地震数据的单频信息，数据矩阵的正确表示和构建，是正确多次波预测的重要前提和保障；

5)在每个频率分量上对所形成的数据矩阵由CPU传到GPU上，并在GPU上完成乘法运算；

所述的数据矩阵的乘法运算，充分利用了图形处理器(Graphic Processing Unit，GPU)在矩阵乘法方面的运算优势，大幅度地提高算法的计算效率，使得全三维表面多次波预测算法在实际应用成为一种计算速度快、成本低、可应用的技术方案。

所述的数据矩阵乘法运算利用CUDA(统一计算设备架构)为编程平台，其主要思路是利用CUDA语言加速计算三维表面多次波预测算法，突破以往算法在计算量方面的瓶颈。

将预测的多次波数据进行反傅里叶变换到时空域；

6)将加速计算预测的多次波数据，由GPU传到CPU，并在CPU上进行反傅里叶变换计算，得到预测的时空域地震数据；

在上述4)中，在每个频率分量上形成数据矩阵，假定采集是在x-y平面上一矩形区域均匀网格内进行的，在x方向放置N_x个检波器，在y方向放置N_y个检波器，则每一炮会产生N_xN_y道记录，假定在此矩形区域的N_xN_y个位置都放炮，在每个网格点上都做一次炮位置，即有N_xN_y个炮记录，对每炮记录，有N_xN_y道记录能被测量到，结果会得到(N_xN_y)²道数据。同2D方法一样，将数据变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集。

在本发明中，充分利用了GPU近年发展起来的在科学计算方面高速度和高性能的特点，应用CUDA编程语言，实现表面多次波预测的并行计算。GPU的使用可达到优异的加速比并获得较高的并行计算效率。由于GPU具有并行的计算优势，在GPU上实现矩阵相乘以加速运算，提高多次波预测的速度。而多次波压制的其他运算由CPU顺序执行。利用GPU和CPU(CPPC)协同并行计算预测表面多次波。具体来说，协同并行计算就是由CPU负责执行顺序型的代码，而将密集运算的矩阵乘法转移到GPU上进行并行计算。由于GPU的强大运算能力，该方法可在较短的时间内，解决大数据体多次波预测的问题。本发明可节省大量的计算时间。

在5)中，调用统一计算设备架构平台CUDA的函数库matrix_mul计算频率空间域的数据矩阵乘法运算。

本发明的实际应用效率与具体的图形处理器硬件相关，在一定范围内，图形处理器所拥有的流处理器越多、速度越快、速度提升的效果将越好。下面结合附图进一步进行说明：

本发明利用基于图形处理器的三维表面多次波预测的方法对由图3所示速度模型正演的三维水平层状介质表面多次波数据预测多次波，其中图6某点放炮模拟的地震数据；图7全3D算法预测的多次波；图8借助于某多次波相减算法获得的多次波压制结果。

参考附图1，地下介质的真实结构为三维，因此不同地层之间的分界面具有三维结构，这表明在分界面上产生的一次波和多次波具有三维的传播和反射效应。当考虑一阶表面多次波时，自由表面的反射点可能位于自由表面的任何位置，与产生多次波反射面的形状有关。考虑相对较为简单的情况，假定反射面在cross-line方向具有倾角，图1所示为一阶多次波的射线路径图，可观察到表面的反射点并没有位于连接炮点和检波点的直线上，而是沿着倾角方向在横向上有移动。

参考附图2，全三维多次波预测算法描述了地震数据空间2D褶积过程。这个过程也要求在联络测线和主测线方向震源和检波器需要有同样密度的采样，即炮间距等于道间距。图2所示的观测系统中，若灰色的圈代表检波点位置，那么在每个检波点位置都要放置震源，图中仅用三个震源示意。

参考附图3，3D水平层状介质表面多次波正演数据的速度模型，通过这个模型来测试本发明的适用性与高效性。

参考附图4，公开了本发明基于图形处理器的三维表面多次波预测方法的流程图，图4中用设备一词指代图形处理器。包括下述步骤：(1)使用海上专用的电缆和检波器，在观测船航行中连续进行地震波的激发和接收，采集含自由表面多次波的三维地震数据；(2)对采集到的三维地震数据，采用Radon变换的方法实现地震数据规则化，使得规则化后地震数据的炮间距与道间距相等；(3)将时空域数据变换到频率空间域；(4)在每个频率分量上形成数据矩阵；(5)在每个频率分量上对数据矩阵做乘法运算；(6)将预测的多次波数据进行反傅里叶变换到时空域；上述步骤(5)在图形处理器中进行。

参考附图5，为本发明基于图形处理器分块计算矩阵乘法。假定计算维度分别为(WA，HA)和(WB，WA)的矩阵A和矩阵B的乘积C，则每个线程块均负责计算C的一个子方阵C_sub，块内每个线程负责计算C_sub的一个元素。

参考附图6所示，利用本发明对图3所示三维水平层状介质进行含表面多次波波场正演模拟的模型数据，该图为理论模拟的某炮数据56条测线的含一阶表面多次波的地震记录，其中0.55s和0.8s附近为一次波能量，1.0s和1.3s附近为多次波能量。有效的多次波预测方法可以正确地预测1.0s和1.3s附近为多次波。

参考附图7所示，利用本发明对6所示的多次波模型数据进行全三维多次波预测的结果，目前，SRME仍主要局限于2D算法。在地震数据处理中，经常利用2D算法预测表面多次波，然而，在许多情况下，对地下介质和观测系统所做的2D假设在某种程度上是无效的，使得预测的多次波出现了很大的误差，难以在自适应相减中有效地减去。

参考附图8所示，三维含多次波数据的压制效果。为了验证本发明的多次波预测效果，采用伪多道自适应匹配相减的方法，将预测的多次波从含多次波的总波场中减去，可获得理想的多次波压制效果。

显然上述实施例仅为清楚的说明本发明所做的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上，还可以做出其他不同形式的变化或变动，这里无需也无法对所有实施方式予以穷举。由此所引申的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (4)

1.一种地震勘探中的数据波形处理方法，其具体为基于图形处理器利用反馈迭代法计算三维表面多次波预测方法，其采用以下步骤：

1）使用海上专用的电缆和检波器，在观测船航行中连续进行地震波的激发和接收，采集含有自由表面多次波的三维地震数据；

2）对采集到的三维地震数据，采用Radon变换的方法实现地震数据规则化，使得规则化后地震数据的炮间距与道间距相等，进而得到规则化后的时空域数据，为反馈迭代法计算准备地震数据；

3）利用傅里叶变换将规则化后的时空域数据变换到频率空间域；

4）循环读取各个频率分量上的数据，在每个频率分量上形成数据矩阵：

假定在x方向放置N_x个检波器，在y方向放置N_y个检波器，则每炮数据有N_xN_y个检波器接收，即每一炮会产生N_xN_y道记录，假定在此矩形区域的N_xN_y个位置都放炮，每个网格点做一次震源位置，每炮记录在N_xN_y个检波器上接收到的信号被记录下来，结果会得到N_xN_yN_xN_y道数据，就是（N_xN_y）²道数据；将所有地震数据道变换到频域后，可在每个频率分量上构建数据矩阵，每个数据矩阵将包含单频分量上的全3D采集信息；在每个震源位置都做上述构建数据矩阵的重复，结果得到大小为N_xN_y的方阵；所形成的数据矩阵中，矩阵的每列为3D共炮点道集，每行为3D共接收点道集；

5）在每个频率分量上对数据矩阵做乘法运算：

在每个频率分量上，应用步骤4）形成的数据矩阵完成矩阵的自乘运算，使用矩阵乘法表述来表示x和y方向二维空间褶积；

6）将预测的多次波数据进行反傅里叶变换到时空域；

上述步骤5）在图形处理器中进行。

2.根据权利要求1所述的地震勘探中的数据波形处理方法，其中反馈迭代法预测多次波的过程是地震数据与地震数据本身在时空域的褶积。

3.根据权利要求1所述的地震勘探中的数据波形处理方法，其特征在于所述步骤（5）中在每个频率分量上对数据矩阵做乘法，利用统一计算设备架构平台CUDA的函数matrix_mul实现。

4.根据权利要求3所述的地震勘探中的数据波形处理方法，其特征在于，应用统一计算设备架构平台CUDA编程语言实现图形处理器GPU通用计算，其具体的步骤和方法是，在首个低频频率分量上，将数据矩阵由CPU传送到图形处理器GPU上，在图形处理器GPU上完成数据矩阵的乘法运算，再将数据矩阵的乘积结果传送回CPU上，完成此频率分量的多次波预测，对其他的频率分量，依次执行此操作，可将所有频率分量上的矩阵乘法执行完毕，从而完成表面多次波预测。

Images