叠前高角度快速傅里叶变换地震成像方法、系统、设备

1.本发明属于地球物理勘探、地震成像领域，具体涉及一种叠前高角度快速傅里叶变换地震成像方法、系统、设备。


背景技术：

2.随着勘探的不断深入，极复杂断块油藏、走滑断裂油藏、裂缝性潜山等油气藏成为油气勘探亟需突破的新领域，与此相对应的高陡断块构造、高陡走滑断裂和裂缝性潜山构造，产生地震反射具有高角度、弱振幅等特征。传统基于born近似的一阶逼近分裂步傅里叶变换（ssf）高效地震成像技术仅适用于缓倾角（《15度）地层成像，已无法满足目前勘探地质目标的需求。新世纪以来发展的基于双程波的最小二乘逆时偏移等高精度成像技术，虽然可有效改善高陡构造的成像品质，但存在计算效率极低、高度依赖速度模型精度等缺点而无法实现工业化应用。因此，迫切需要一种能够对高角度等复杂构造进行快速高精度成像的技术。


技术实现要素：

3.为了解决现有技术中的上述问题，即为了解决现有的地震成像方法在对高陡构造（即高角度地层）进行成像时，存在的成像质量、效率低以及过于依赖速度模型精度问题，本发明第一方面，提出了一种叠前高角度快速傅里叶变换地震成像方法，该方法包括以下步骤：s100，获取在地震勘探过程中采集的地震数据；所述地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距；s200，在所述地震数据中，按照炮点和检波点的坐标，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，并将所述时空域共偏移距道集中的时空域共偏移距地震数据利用傅里叶变换由时空域转换至频率波数域，得到频率波数域共偏移距地震数据；s300，基于所述频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度，并根据所述波传播角度，将成像区域分为小角度区域和大角度区域；s400，通过预设的改进的随机游走算法，分别求解小角度区域与大角度区域对应的最优化的常系数；s500，将求解的常系数代入频率波数域波场延拓方程，对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果；s600，将不同深度的不同频率的成像结果进行积分，并利用傅里叶逆变换，将积分后的成像结果转换至空间域，得到共偏移距成像剖面；s700，将所有共偏移距成像剖面进行叠加处理，获得最终的成像结果。
4.在一些优选的实施方式中，根据所述波传播角度，将成像区域分为小角度区域和大角度区域，其方法为：根据地震波波场对、方向的偏导数，获得地震波传播的角度，若
所述地震波传播的角度小于设定角度阈值，则将地震成像区域划分为小角度区域，否则划分为大角度区域。
5.在一些优选的实施方式中，所述预设的改进的随机游走算法其求取有理式系数方法为：其中，表示目标函数，表示地震波传播的角度，表示不同介质的折射率，、、表示有理常系数，表示横向波数，，表示归一化的横向波数，表示背景总波数，表示设定的最大成像角度。
6.在一些优选的实施方式中，所述频率波数域波场延拓方程，其构建方法为：基于真实地下介质的非均匀介质性，构建有理式逼近的波数频散关系；将所述有理式逼近的波数频散关系，代入频率波数域的傅里叶波场延拓基本方程，得到频率波数域波场延拓方程；所述频率波数域波场延拓方程为：其中，表示深度为的波场信息，表示纵向波数，表示深度为的波场信息，，表示中间参数变量，表示虚数，表示深度方向空间采样间隔。
7.在一些优选的实施方式中，所述有理式逼近的波束频散关系为：其中，表示归一化的纵向波数。
8.本发明的第二方面，提出了一种叠前高角度快速傅里叶变换地震成像系统，该系统包括：地震数据获取模块、转换模块、角度区域划分模块、常系数求解模块、波场延拓模块、逆变换模块、剖面叠加模块；所述地震数据获取模块，配置为获取在地震勘探过程中采集的地震数据；所述地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距；所述转换模块，配置为在所述地震数据中，按照炮点和检波点的坐标，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，并将所述时空域共偏移距道集中的时
空域共偏移距地震数据利用傅里叶变换由时空域转换至频率波数域，得到频率波数域共偏移距地震数据；所述角度区域划分模块，配置为基于所述频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度，并根据所述波传播角度，将成像区域分为小角度区域和大角度区域；所述常系数求解模块，配置为通过预设的改进的随机游走算法，分别求解小角度区域与大角度区域对应的最优化的常系数；所述波场延拓模块，配置为将求解的常系数代入频率波数域波场延拓方程，对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果；所述逆变换模块，配置为将不同深度的不同频率的成像结果进行积分，并利用傅里叶逆变换，将积分后的成像结果转换至空间域，得到共偏移距成像剖面；所述剖面叠加模块，配置为将所有共偏移距成像剖面进行叠加处理，获得最终的成像结果。
9.本发明的第三方面，提出了一种设备，至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的叠前高角度快速傅里叶变换地震成像方法。
10.本发明的第四方面，提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的叠前高角度快速傅里叶变换地震成像方法。
11.本发明的有益效果：本发明实现了对高陡构造的高精度地震成像，提升了成像效率。
12.1）本发明通过获得高角度有理式傅里叶变换成像算子（即高角度偏移算子），将其引入到频率波数域波场延拓方程中，对频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果，实现了对高陡构造的高精度成像，其成像角度可达60
°
，远高于传统ssf地震成像方法；2）本发明可对共偏移距道集进行地震成像，计算效率高，可进行工业化实际应用，也可以广泛应用于勘查业相关的地球物理地震成像领域。
附图说明
13.通过阅读参照以下附图所做的对非限制性实施例所做的详细描述，本技术的其他特征、目的和优点将会变得更明显。
14.图1是本发明一种实施例的叠前高角度快速傅里叶变换地震成像方法的流程示意图；图2是本发明一种实施例的叠前高角度快速傅里叶变换地震成像系统的框架示意图；图3是本发明一种实施例的叠前高角度快速傅里叶变换地震成像方法的简略示意图；图4(a)是本发明一种实施例的折射率为0.3的波场快照精度的示意图；
图4(b)是本发明一种实施例的折射率为0.6的波场快照精度的示意图；图4(c)是本发明一种实施例的折射率为0.9的波场快照精度的示意图；图5(a)是本发明一种实施例的seg/eaeg岩丘速度模型的示意图；图5(b)是本发明一种实施例的构造的小角度区域和大角度区域的示意图；图6(a)是传统kirchhoff技术的成像结果示意图；图6(b)是传统ssf技术的成像结果示意图；图6(c)是本发明rfa技术的成像结果示意图；图7(a)是传统kirchhoff成像剖面的示意图；图7(b)是本发明方法成像剖面的示意图；图8是本发明一种实施例的适于用来实现本技术实施例的电子设备的计算机系统的结构示意图。
具体实施方式
15.为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
16.下面结合附图和实施例对本技术作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
17.需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的特征可以相互组合。
18.本发明的一种叠前高角度快速傅里叶变换地震成像方法，如图1所示，该方法包括以下步骤：s100，获取在地震勘探过程中采集的地震数据；所述地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距；s200，在所述地震数据中，按照炮点和检波点的坐标，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，并将所述时空域共偏移距道集中的时空域共偏移距地震数据利用傅里叶变换由时空域转换至频率波数域，得到频率波数域共偏移距地震数据；s300，基于所述频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度，并根据所述波传播角度，将成像区域分为小角度区域和大角度区域；s400，通过预设的改进的随机游走算法，分别求解小角度区域与大角度区域对应的最优化的常系数；s500，将求解的常系数代入频率波数域波场延拓方程，对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果；s600，将不同深度的不同频率的成像结果进行积分，并利用傅里叶逆变换，将积分后的成像结果转换至空间域，得到共偏移距成像剖面；s700，将所有共偏移距成像剖面进行叠加处理，获得最终的成像结果。
19.为了更清晰地对本发明叠前高角度快速傅里叶变换地震成像方法进行说明，下面结合附图3对本发明方法一种实施例中各步骤进行展开详述。
20.本发明提出一种新的基于有理式逼近的双程波快速傅里叶变换地震成像技术，在频率波数域实施，区别于传统基于born近似的单程波快速傅里叶变换地震成像技术，实现了对双程波的高角度近似，进而获得了高角度有理式傅里叶变换成像算子，将其引入到波场延拓公式中，实现了对高陡构造的高精度成像。具体过程如下：s100，获取在地震勘探过程中采集的地震数据；所述地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距；在本实施例中，先获取在地震勘探采集的地震数据，地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距。其中，本发明中设计的速度模型为seg/eaeg岩丘速度模型，如图5(a)所示，最大构造角度约为70
°
，速度模型大小优选设置为15000 m
ꢀ×ꢀ
3700 m，主频为10hz。
21.s200，在所述地震数据中，按照炮点和检波点的坐标，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，并将所述时空域共偏移距道集中的时空域共偏移距地震数据利用傅里叶变换由时空域转换至频率波数域，得到频率波数域共偏移距地震数据；在本实施例中，从检波器阵列采集的地震数据中提取共炮点道集，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，通过傅里叶变换（即图3中的fft）将所述时空域共偏移距道集中的时空域共偏移距地震数据由时空域转换至频率波数域，得到频率波数域共偏移距地震数据，表示截止频率。
22.s300，基于所述频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度，并根据所述波传播角度，将成像区域分为小角度区域和大角度区域；在本实施例中，先根据频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度（即地震波传播的角度），然后根据波传播角度，将成像区域划分为小角度区域和大角度区域，具体为：根据地震波波场对、方向的偏导数，可以获得地震波传播的角度，若所述地震波传播的角度小于设定角度阈值，则将地震成像区域划分为小角度区域，否则划分为大角度区域。
23.图5(b)为构造的小角度区域和大角度区域的示意图，其中，图5(b)中黑色为构造的小角度区域，白色为构造的大角度区域。
24.s400，通过预设的改进的随机游走算法，分别求解小角度区域与大角度区域对应的最优化的常系数；s500，将求解的常系数代入频率波数域波场延拓方程，对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果；在频率波数域，傅里叶波场延拓基本方程为： （1）
其中，，分别表示横向波数、纵向波数，其满足，表示背景总波数，表示深度为的波场信息，表示深度z的傅里叶波场，表示虚数，表示深度方向空间采样间隔。因此，通过公式（1）可以获得不同深度的波场信息。
25.由于真实地下介质为非均匀介质，其波数频散关系可表示为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)其中，为不同介质的折射率，，，其中，表示归一化的横向波数，表示归一化的纵向波数。由上式可以看出，但时，波数频散关系将产生不稳定。其中，不同折射率波场快照精度的对比图，其中，图4(a)、图4(b)、图4(c)分别为折射率为0.3、0.6、0.9的波场快照精度的示意图。为解决该问题，本发明提出以下基于有理式逼近（rational function approximation，简称rfa）的波数频散关系：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（3）将公式（3）带入公式（1）中，可得频率波数域波场延拓的基本方程：（4）其中，、、表示有理常常系数，，表示深度为的波场信息。本发明定义最小二乘函数为：
ꢀꢀ
（5）其中，表示目标函数，为本发明近似频散曲线与公式（2）频散曲线的二范数。通过求取目标函数的最小值可获得最优有理常系数、、，表示地震波传播的角度，表示设定的最大成像角度，由确定。本发明采用改进的随机游走算法进行全局寻优。
26.在本实施例中，通过预设的改进的随机游走算法（即改进的随机游走算法其求取有理式系数方法），分别求解小角度区域与大角度区域对应的最优化的常系数，将求解的常系数代入频率波数域波场延拓方程，获得高角度偏移算子（即公式（4）中大括号的内容），对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像
结果。对所述seg/eaeg 岩丘模型小角度计算区域优化系数分别为、、；大角度区域对应的优化系数分别为、、。
27.s600，将不同深度的不同频率的成像结果进行积分，并利用傅里叶逆变换，将积分后的成像结果转换至空间域，得到共偏移距成像剖面；在本实施例中，将不同频率的成像结果进行积分，并将积分后的成像结果进行傅里叶逆变换（即图3中的ifft），得到共偏移距成像剖面。
28.s700，将所有共偏移距成像剖面进行叠加处理，获得最终的成像结果。
29.在本实施例中，如图6(c)所示，将所有共偏移距成像剖面进行叠加处理，获得本发明的最终成像结果。
30.另外，对岩丘模型开展数值模拟实验，验证了本发明能够在不增加计算量前提下实现对高陡构造的高品质成像，首先以图5(a)所示的seg/eaeg岩丘速度模型为例，图6(a)为传统kirchhoff技术的成像结果，图6(b)为传统ssf技术的成像结果，图6(c)表示本发明rfa技术的成像结果。从图中可以看出，本发明方法能够对高角度实现高分辨率成像（图6(c)箭头所示），传统单程波偏移成像技术无法实现对高陡构造的准确成像。
31.为了进一步说明本专利提出技术的效果，开展了实际资料测试。通过本发明方法成像剖面的示意图（如图7(b)）与传统kirchhoff成像剖面的示意图（如图7(a)）对比可以看出，本发明技术能够对高陡构造实现高精度成像。
32.本发明第二实施例的一种叠前高角度快速傅里叶变换地震成像系统，如图2所示，具体包括以下模块：地震数据获取模块100、转换模块200、角度区域划分模块300、常系数求解模块400、波场延拓模块500、逆变换模块600、剖面叠加模块700；所述地震数据获取模块100，配置为获取在地震勘探过程中采集的地震数据；所述地震数据包括震源类型、震源主频、炮点的位置、最大构造角度、设定的采样时间、给定的速度模型的大小及网格点的间距；所述转换模块200，配置为在所述地震数据中，按照炮点和检波点的坐标，将共炮点道集转换为共偏移距道集，作为时空域共偏移距道集，并将所述时空域共偏移距道集中的时空域共偏移距地震数据利用傅里叶变换由时空域转换至频率波数域，得到频率波数域共偏移距地震数据；所述角度区域划分模块300，配置为基于所述频率波数域共偏移距地震数据，计算不同频率下速度模型所有网格点的波传播角度，并根据所述波传播角度，将成像区域分为小角度区域和大角度区域；所述常系数求解模块400，配置为通过预设的改进的随机游走算法，分别求解小角度区域与大角度区域对应的最优化的常系数；所述波场延拓模块500，配置为将求解的常系数代入频率波数域波场延拓方程，对所述频率波数域共偏移距道集进行分频率逐层波场延拓，获得不同深度的不同频率的成像结果；
所述逆变换模块600，配置为将不同深度的不同频率的成像结果进行积分，并利用傅里叶逆变换，将积分后的成像结果转换至空间域，得到共偏移距成像剖面；所述剖面叠加模块700，配置为将所有共偏移距成像剖面进行叠加处理，获得最终的成像结果。
33.所述技术领域的技术人员可以清楚的了解到，为描述的方便和简洁，上述描述的系统具体的工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
34.需要说明的是，上述实施例提供的叠前高角度快速傅里叶变换地震成像系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。
35.本发明第三实施例的一种设备，至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的叠前高角度快速傅里叶变换地震成像方法。
36.本发明第四实施例的一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的叠前高角度快速傅里叶变换地震成像方法。
37.所述技术领域的技术人员可以清楚的了解到，未描述的方便和简洁，上述描述的设备、计算机可读存储介质的具体工作过程及有关说明，可以参考前述方法实例中的对应过程，在此不再赘述。
38.下面参考图8，其示出了适于用来实现本技术系统、方法、设备实施例的服务器的计算机系统的结构示意图。图8示出的服务器仅仅是一个示例，不应对本技术实施例的功能和使用范围带来任何限制。
39.如图8所示，计算机系统包括中央处理单元（cpu，central processing unit）801，其可以根据存储在只读存储器（rom，read only memory）802中的程序或者从存储部分808加载到随机访问存储器（ram，random access memory）803中的程序而执行各种适当的动作和处理。在ram803中，还存储有系统操作所需的各种程序和数据。cpu801、rom 802以及ram803通过总线804彼此相连。输入/输出（i/o，input/output）接口805也连接至总线804。
40.以下部件连接至i/o接口805：包括键盘、鼠标等的输入部分806；包括诸如阴极射线管、液晶显示器等以及扬声器等的输出部分807；包括硬盘等的存储部分808；以及包括诸如局域网卡、调制解调器等的网络接口卡的通讯部分809。通讯部分809经由诸如因特网的网络执行通讯处理。驱动器810也根据需要连接至i/o接口805。可拆卸介质811，诸如磁盘、光盘、磁光盘、半导体存储器等等，根据需要安装在驱动器810上，以便于从其上读出的计算机程序根据需要被安装入存储部分808。
41.特别地，根据本公开的实施例，上文参考流程图描述的过程可以被实现为计算机软件程序。例如，本公开的实施例包括一种计算机程序产品，其包括承载在计算机可读介质上的计算机程序，该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中，该计算机程序可以通过通讯部分809从网络上被下载和安装，和/或从可拆卸介质
811被安装。在该计算机程序被cpu801执行时，执行本技术的方法中限定的上述功能。需要说明的是，本技术上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是但不限于：电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、ram、rom、可擦式可编程只读存储器（eprom或闪存）、光纤、便携式紧凑磁盘只读存储器（cd-rom）、光存储器件、磁存储器件、或者上述的任意合适的组合。在本技术中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本技术中，计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：无线、电线、光缆等等，或者上述的任意合适的组合。
42.可以以一种或多种程序设计语言或其组合来编写用于执行本技术的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言，如java、smalltalk、c ，还包括常规的过程式程序设计语言，如c语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络，包括局域网或广域网连接到用户计算机，或者可以连接到外部计算机（例如利用因特网服务提供商来通过因特网连接）。
43.附图中的流程图和框图，图示了按照本技术各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
44.术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
45.术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
46.至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些
更改或替换之后的技术方案都将落入本发明的保护范围之内。