具体实施方式

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

首先对本说明书实施例中所使用的主频近似方法和流程具体进行解释说明。

时空域中由高斯束表征的正向波场W^(F)(x₀,t；x_s)可写为：

其中，τ为中心射线的走时，x₀＝(x₀,z₀)代表成像点的空间坐标,x_s＝(x_s,0)代表炮点的空间坐标,ω是角频率，ε是初始束参数。P(s)和Q(s)是基于观测数据和地震波在地下介质传播的初始速度v₀，以及动力学射线追踪方程的基本解，n是射线的横向分量。

基于此，可以利用观测数据的主频来简化与旅行时虚部相关的正向波场的计算：

其中，和分别代表正向波场的实部和虚部，≈为本发明中主频近似的体现。

把方程(2)代入到方程(1)，则得到W^(F)(x₀,t；x_s)的近似表达式：

其中，

其中，ω_m是观测数据的主频。

对方程(3)作傅里叶变换得:

然后，我们使用上行射线追踪策略来构建时间反传波场，记录观测数据P_U(x_r,t)从接受点x_r到地下成像点x₀的反向传播可以用Kirchhoff积分实现：

其中，G(x_r,t；x₀,t₀)是格林函数,它是射线偏移方法中震源与接收点之间的重要的波场映射函数。利用一系列高斯束的叠加去近似格林函数：

在高频近似的条件下，格林函数的导数表达式可进一步被简化为：

其中，p_z表示接收点的垂直射线慢度。

把式(7)、(8)代入到式(6)中,我们可以得到反传时间波场的表达式：

我们对式(9)右侧的时间域表示部分通过作傅里叶变换得到频率域部分，W^(R)(x₀,t₀)可以重新表示为：

同样地,利用观测数据的主频来简化与旅行时虚部相关的反向波场的计算：

其中，≈为本发明中主频近似的体现。

把式(11)代入到式(10)中,W^(R)(x₀,t₀)能被进一步简化为：

其中，

对方程(12)作反傅里叶变换，可得：

其中，

g(x_r,t)＝∫ωP_U(x_r,ω)exp(iωt)dω (15)

最后，利用互相关成像条件来计算成像结果：

把式(5)、(14)代入到式(16)得：

对于上行和下行射线的角度间隔，我们借鉴传统频率域高斯束偏移成像方法中心射线的角度间距表达式：

其中，ω_ref是观测数据的参考频率，ω_hig是观测数据的最高频率。

在循环算法中，我们需要计算上行和下行射线的数量：

其中，a_max是最大偏移角，a_min是最小偏移角。通过该方式可以基于观察数据的主频在成像过程中进行更便捷的的迭代求和。

上述部分对于本说明书实施例中所采用的主频近似原理和成像流程进行了具体说明。基于前述内容，本说明书的实施例提供了一种基于主频近似的声波介质时空域高斯束偏移成像方法，如图1所示，该过程具体包括以下步骤：

S101，输入速度参数场v和观测数据P_U(x_r,t)，其中，x_r代表接收点的空间坐标，t是地震波传播的时间；

S103，确定所述观察数据的主频ω_m；

S105，根据所述速度参数场v，通过对主频近似确定正向波场W^(F)(x₀,t；x_s)和时间反传波场W^(R)(x₀,t₀)的虚部，其中，x₀代表成像点的空间坐标,x_s代表炮点的空间坐标，t为时间，t₀为波场传播的初始值。

具体而言，采用如下公式确定正向波场的虚部：

其中，为正向波场的虚部，为空间方位角，Im表示复数的虚部，是正向波场走时的实部，Re表示复数的实部，τ为正向波场走时的初初始值，表示正向波场的振幅，δ表示脉冲函数，s是射线的弧长分量，s₀是射线弧长分量的初始值，ω是角频率，ε是初始束参数，P(s)和Q(s)是动力学射线追踪方程的基本解，v₀是地震波在地下介质传播的初始速度，n是射线的横向分量；

以及，采用如下公式确定时间反传波场的虚部：

其中，是时间反传波场走时的虚部，是时间反传波场走时的实部，是时间反传波场振幅实部的共轭，g(x_r,t)为自定义参数。

S107，根据所述正向波场和时间反传波场的虚部进行互相关生成成像结果。

即采用如下互相关公式生成成像结果：

其中，I(x₀)为在空间点x₀的成像结果，p_z是接收点的垂直射线慢度。

与现有技术相比，本发明既能保持与传统时空域高斯束偏移方法相比拟的成像精度和分辨率，又大大提高了时空域高斯束偏移方法的实现效率。开发基于主频近似的声波介质时空域高斯束偏移成像技术，为复杂构造区域的地震数据处理提供了高精度成像保障，且具有更少的时间成本消耗，提高了后续解释工作的质量与效率。

以下给出本说明是实施例在模型中的实际效果说明。

将本发明所提供的方法首先应用于一个简单的洼陷模型成像，取得了较理想的成像效果。真实速度模型(如图2所示)；建立移动接收的观测系统，输入平滑P波速度场和线性正演模拟得到的观测炮记录；采用互相关成像条件将正传和反传波场互相关，得到传统频率域高斯束偏移成像结果(如图3所示)、传统时空域高斯束偏移成像结果(如图4所示)、本说明书实施例所提供的成像结果(如图5所示)。

在图3中，频率域高斯束偏移方法产生一些偏移假象(红色箭头)。在图4和图5中，我们可以看到，由于时间反传波场采用上行射线追踪策略，所有反射界面都能清晰成像，且精度相当。对比两种时空域高斯束偏移方法的运行时间，本说明书实施例所提供的成像方法可以将洼陷模型(图2)的计算效率提高136.0倍。

将本发明所提供的方法最后应用于国际标准Marmousi模型成像，取得了不错的成像效果。真实速度模型(如图6所示)；建立移动接收的观测系统，输入平滑P波速度场和线性正演模拟得到的观测炮记录；采用互相关成像条件将正传和反传波场互相关，得到传统频率域高斯束偏移成像结果(如图7所示)、传统时空域高斯束偏移成像结果(如图8所示)、本说明书实施例所提供的成像结果(如图9所示)。

可以看出，两种时空域高斯束偏移方法对浅层的成像质量高于频率域高斯束偏移方法(如图7、8、9蓝色椭圆和红色矩形所示部分)。从图8、9中可以看出，两种时空域高斯束偏移方法成像精度大致相同。对比两种时空域高斯束偏移方法的运行时间，本说明书实施例所提供的成像方法可以将Marmousi模型(图7)的效率提高39.9倍。

本说明书实施例所提供的成像方法在保证成像精度和分辨率的同时，大大提高了传统时空域高斯束偏移方法的计算效率，为时空域高斯束偏移成像方法的发展与应用提供了有利的条件，为复杂构造区域的解释工作提供较准确的成像基础，为老油田二次勘探开发提供强有力的技术支持。

对应的，本申请实施例还提供一种计算机设备，所述设备包括包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述程序时实现如前所述的基于主频近似的声波介质时空域高斯束偏移成像所述的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置、设备和介质类实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可，这里就不再一一赘述。

上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤或模块可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。