CN110780348B - 基于立体成像条件的一次波和多次波联合成像方法及系统 - Google Patents

Abstract

本发明公开了一种基于立体成像条件的一次波和多次波联合成像方法及系统，包括：获取地表观测点的地震数据，利用表层多次波去除方法从地震数据中提取震源子波；根据地震数据和震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，根据地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，利用局部倾斜叠加方法分别对延拓后震源波场和延拓后检波波场进行分解，根据设定的速度场，利用平面波解构滤波器确定分解震源波场和检波波场局部斜率的取值范围；在局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。本发明通过考虑时间和局部斜率对成像互相关的约束，提高成像质量。

Description

基于立体成像条件的一次波和多次波联合成像方法及系统

技术领域

本发明涉及地震成像技术领域，特别是涉及基于立体成像条件的一次波和多次波联合成像方法及系统。

背景技术

长期以来，地震数据中发育的多次波都被视为噪音，在偏移成像前需对其进行压制去除。然而多次波也是来自地下界面的真实反射，具有更广的照明范围，可以提供更丰富的地球内部精细结构信息。传统的多次波成像方法需要首先进行多次波分离，这是个耗时耗力的过程。为了节省计算成本，部分学者提出了一种一次波和多次波联合成像的思路。

目前发展的方法多通过选取脉冲震源与原始记录之和作为震源波场，利用单程波偏移结合二维反褶积成像条件实现一次波和多次波的联合成像。但是该方法中脉冲震源子波的选择存在一定的问题，子波能量过小会影响一次反射波的成像效果，而能量过大则会造成多次波成像信息被一次波成像结果所掩盖，且子波选取不当会带来成像位置不准确等问题；此外，成像过程中不同阶多次波引起的不相干波场互相关会引入大量的串扰假象，严重影响成像的质量，二维反褶积成像条件可以压制一定的串扰噪声，但是对复杂构造的适应性较差，降低了传统方法的应用价值。

发明内容

基于此，本发明提供一种基于立体成像条件的一次波和多次波联合成像方法及系统，从而提高一次波和多次波联合成像的质量。

为实现上述目的，本发明提供了如下方案：

一种基于立体成像条件的一次波和多次波联合成像方法，所述成像方法包括：

获取地表观测点的地震数据，所述地震数据包括一次波和多次波；

利用表层多次波去除方法从所述地震数据中提取震源子波；

根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场；

根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场；

利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场；

根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围；

在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

可选的，所述利用表层多次波去除方法从所述地震数据中提取震源子波，具体包括：

根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子；

根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

可选的，所述根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子，具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

可选的，所述单程波波动方程具体为：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

可选的，所述在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

一种基于立体成像条件的一次波和多次波联合成像系统，所述成像系统包括：

地震数据获取模块，用于获取地表观测点的地震数据，所述地震数据包括一次波和多次波；

震源子波提取模块，用于利用表层多次波去除方法从所述地震数据中提取震源子波；

延拓后震源波场确定模块，用于根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场；

延拓后检波波场确定模块，用于根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场；

波场分解模块，用于利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场；

局部斜率取值范围确定模块，根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围；

成像结果确定模块，用于在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

可选的，所述震源子波提取模块，具体包括：

地表加权算子确定单元，根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子；

震源子波确定单元，用于根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

可选的，所述地表加权算子确定单元，具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

可选的，所述单程波波动方程具体为：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

可选的，所述成像结果确定模块，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供的基于立体成像条件的一次波和多次波联合成像方法及系统，利用表层多次波去除方法从地震数据中提取震源子波，为高质量成像提供了高精度的子波；在估计的地层局部斜率基础上，同时利用延拓波场的传播时间和传播方向来约束成像过程，基于多参数约束的立体成像条件来实现一次波和多次波的联合成像，提高一次波和多次波联合成像的质量。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例一种基于立体成像条件的一次波和多次波联合成像方法的流程图；

图2为本发明实施例简单层状模型速度场的示意图；

图3为本发明实施例提取的震源子波与真实子波对比图；

图4为本发明实施例简单层状模型单炮地震成像结果图；

图5为本发明实施例复杂模型速度场的示意图；

图6为本发明实施例复杂模型多炮地震成像结果图；

图7为本发明实施例一种基于立体成像条件的一次波和多次波联合成像系统的结构示意图。

具体实施方式

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

本发明的目的是提供一种基于立体成像条件的一次波和多次波联合成像方法及系统，从而提高一次波和多次波联合成像的质量。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明实施例一种基于立体成像条件的一次波和多次波联合成像方法的流程图。参见图1，一种基于立体成像条件的一次波和多次波联合成像方法，所述成像方法包括：

S1：获取地表观测点的地震数据，所述地震数据包括一次波和多次波。

S2：利用表层多次波去除方法从所述地震数据中提取震源子波。

S2具体包括：

S201：根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子。

S202：根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

具体的，利用表层多次波去除(SRME，Surface-Related Multiple Elimination)方法预测滤波因子，即地表加权算子，利用傅里叶变换提取震源子波。

假设地震数据已消除震源和检波器效应，且地表反射系数为-1,则基于反馈迭代理论的SRME方法确定地表加权算子的具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

A(w)＝S^-1(w)，S(w)为震源子波的傅里叶变换，假设其仅与时间有关，则A(w)即表示为逆震源子波的傅里叶变换。通过反馈迭代法可以在未知震源子波的情况下通过最小二乘反演理论迭代实现多次波的预测，同时多次迭代以后可以求得逆震源子波A(w)，对其进行求逆并求取傅里叶反变换即可获取震源子波，表示如下：

s(t)＝IFFT(A^-1(w))。

目前，多数子波提取方法都是在某些假设条件下基于信号分析理论实现的，且要求处理的数据不含多次波。本发明基于波动理论利用叠前数据中的多次波实现震源子波的提取更为可靠，且适用于含多次波数据的子波提取。

S3：根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场。

S4：根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场。

所述S3和S4具体为：

将提取的震源子波与采集的地震数据同时作为虚拟震源，将地震数据作为接收记录，利用单程波波动方程理论进行深度波场延拓。

传统一次波和多次波联合成像方法多通过将任选的脉冲震源与原始记录之和作为虚拟震源，基于单程波理论实现一次波和多次波的联合成像，但是该方法中脉冲震源子波的选择存在一定的问题，震源子波能量过小会影响一次反射波的成像效果，而能量过大则会造成多次波成像信息被一次波成像结果所掩盖，且震源子波选取不当会造成错误的成像位置等问题。

因此，本发明提出一种更为合理的一次波和多次波联合成像方法，首先从地震数据中提取真实的震源子波，将提取的子波和地震数据(同时包含一次波和多次波)代替任选的脉冲震源作为虚拟震源即输入震源波场，将记录的地震数据作为接收波场即检波波场，震源波场和检波波场表示如下：

S(x_r,z₀,t)＝s(x_s,z₀,t)+d(x_r,z₀,t)

D(x_r,z₀,t)＝d(x_r,z₀,t)

其中，S(x_r,z₀,t)表示距离地面z₀深度处位于x＝x_r处的震源波场，D(x_r,z₀,t)表示距离地面z₀深度处位于x＝x_r处的检波波场，s(x_s,z₀,t)表示距离地面z₀深度处位于x＝x_s处的震源子波，d(x_r,z₀,t)表示距离地面z₀深度处位于x＝x_r处的地震数据。

本发明以单程波方程作为波场延拓算子，采用有限差分算法来进行深度延拓，本发明基于以二维声波方程为例

将上式变换到频率域，可得单程波传播Helmhotz方程：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

将上式进行根式展开，采用有限差分法对其进行求解即可实现上下行波场的深度延拓。

S5：利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场。

具体的，将延拓后的震源波场和检波波场利用局部倾斜叠加进行分解，将其分解为关于局部斜率p的一系列分量，震源波场和检波波场表示如下：

S(x,z,t)＝∫W_S(x,z,p,t)dp,

D(x,z,t)＝∫W_D(x,z,p,t)dp,

其中，p表示关于空间和时间的局部斜率三维矢量，S(x,z,t)表示延拓深度为z处的震源波场，D(x,z,t)表示延拓深度z处的检波波场，W_S(x,z,p,t)表示分解的震源波场即表示与局部斜率相关的一系列震源波场，W_D(x,z,p,t)表示分解的检波波场即表示与局部斜率相关的一系列检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

传统成像条件仅通过正反向延拓波场的传播时间进行匹配成像，这也是导致成像假象的弊端所在。考虑到不相关波场互相关成像时具有不同的传播方向，因此，有必要通过在应用成像条件时利用多参数(时间和局部斜率)同时约束成像过程来进行成像串扰的压制。

S6：根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围。

具体的，设定的速度场是人为设定的。

为了提高成像精度和计算效率，本发明引入平面波解构滤波器来进行地层倾角的估计，并根据估计的地层倾角选择一定范围内的局部斜率来约束成像条件中不同斜率的叠加。

分解的震源波场和检波波场在成像过程中仅在局部斜率p与成像地层界面斜率一致时才会产生真实的成像值。因此，局部斜率参数的选取对于成像的精度至关重要。通常，为了提高成像精度，参与叠加的局部斜率的范围应足够大，采样间隔应足够小，而与之伴随的则是高计算成本。为了兼顾成像精度和计算成本，本发明引入平面解构滤波器从速度场中估计反射层的斜率，以此作为成像过程中局部斜率p的约束。通过在反射层斜率附近选取一定范围的数值作为局部斜率的取值范围，该取值范围选取时应视地下构造的复杂情况进行选取，构造越复杂选取范围越大，简单近似平层的情况下根据反射层角度在附近(-π/18,π/18)角度范围内选取即可，这样可以避免了大量局部斜率p循环所导致的高计算成本，同时也保证了成像的精度。

从局部平面波物理模型出发，基于下式的局部平面波微分方程利用平面解构滤波器即可进行局部斜率的估计：

其中，v表示以空间位置为函数的速度场，p_z表示速度场的局部斜率。

S7：在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

与传统成像方法中的成像互相关仅由时间约束是不同的，本发明中成像互相关由时间和局部斜率同时约束，对应成像条件称为立体成像条件。在估计的地层局部斜率基础上，同时利用延拓波场的传播时间和传播方向来约束成像过程，基于多参数约束的立体成像条件来实现一次波和多次波的联合成像，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

该立体成像条件对应成像过程具体描述为：首先基于估计的地层局部斜率，确定局部斜率p的取值范围，然后在此范围内对分解的震源波场和检波波场进行互相关并求和实现最终的成像。该过程中不同斜率成分的分解波场互相关增加了传播方向的限制，有效的避免了串扰噪声的产生，同时在估计的较为准确的地层倾角约束下保证了成像的精度，避免了大量p循环造成的高计算成本。

图2为本发明实施例简单层状模型速度场的示意图，在此模型横向中间位置深度5m处设置一震源，震源子波设定为雷克子波，主频为25Hz，采用中间放炮两边接收观测系统，101道接收，道间距为20m。图3为本发明实施例提取的震源子波与真实子波对比图，参见图3可知，采用本发明方法基于波动理论利用多次波提取的子波与实际子波基本完全吻合，这一结果验证了本发明中子波提取方法的可行性。图4为本发明实施例简单层状模型单炮地震成像结果图，参见图4，是采用不同方法对单炮地震记录成像结果的对比，图4中(a)部分表示任选脉冲震源利用传统互相关成像方法得到的一次波和多次波联合成像结果，选取子波与实际子波的相位差异造成了成像界面不准确，子波能量过强造成多次波成像信息能量较弱，成像结果无法体现利用多次波获取更丰富地下构造信息的优势。同时，不相关波场互相关还引入了大量的串扰噪声。图4中(b)部分应用本发明方法提取震源子波并利用传统互相关成像条件得到的一次波和多次波联合成像结果，一次波和多次波联合成像的优势显而易见，多次波在成像过程中提供了更为丰富的地下构造信息，这也是联合成像相比传统一次波成像的优势，但是成像假象严重降低了成像的质量。图4中(c)部分应用本发明方法基于立体成像条件得到的一次波和多次波联合成像结果中不仅一次波和多次波联合成像优势得到了体现，同时通过时间和局部斜率多参数约束，成像噪声得到了有效的压制，有效的验证了本发明的可行性。

图5为本发明实施例复杂模型速度场的示意图，参见图5，该模型包括平层、背斜、断层等复杂构造，在此模型上布置100炮，起始震源位于(1000m，10m)处，炮间距为20m，震源子波设定为雷克子波，主频为25Hz，采用中间放炮两边接收观测系统，101道接收，道间距为20m。图6为本发明实施例复杂模型多炮地震成像结果图，参见图6，是采用不同方法对所有炮地震数据成像结果的对比，图6中(a)部分表示任选脉冲震源子波利用传统互相关成像方法得到的一次波和多次波联合成像结果，图6中(b)部分应用本发明方法提取子波并利用传统互相关成像条件得到的一次波和多次波联合成像结果，图6中(c)部分应用本发明方法基于立体成像条件得到的一次波和多次波联合成像结果。图6(b)的成像结果和图6(a)成像结果相比，采用真实子波能获取相比传统一次波和多次波联合成像更多的地下构造信息(如图白色竖线右侧成像范围所示)，同时传统互相关成像条件仅通过时间约束成像过程，两图中都存在严重的串扰噪声(如图中白色箭头处所示)，这些假象严重降低了成像的可靠性。图6(c)的成像结果不仅一次波和多次波联合成像优势得到了体现，同时通过在成像过程中增加了局部斜率约束，成像噪声得到了有效的压制(白色箭头所示)，本发明的有效性再次得到了验证。

图7为本发明实施例一种基于立体成像条件的一次波和多次波联合成像系统的结构示意图，参见图7，一种基于立体成像条件的一次波和多次波联合成像系统，所述成像系统包括：

地震数据获取模块201，用于获取地表观测点的地震数据，所述地震数据包括一次波和多次波。

震源子波提取模块202，用于利用表层多次波去除方法从所述地震数据中提取震源子波。

延拓后震源波场确定模块203，用于根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场。

延拓后检波波场确定模块204，用于根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场。

波场分解模块205，用于利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场。

局部斜率取值范围确定模块206，根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围。

成像结果确定模块207，用于在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

所述震源子波提取模块202，具体包括：

地表加权算子确定单元，根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子。

震源子波确定单元，用于根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

所述地表加权算子确定单元，具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

所述单程波波动方程具体为：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

所述成像结果确定模块207，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

本发明提出了一种基于立体成像条件的一次波和多次波联合成像方法及系统，首先基于波动理论表层多次波去除的子波提取方法以获取高精度的震源子波，准确的震源子波保证了成像位置的准确性，能更好的平衡一次波和多次波成像的贡献，然后将提取的震源子波与原始地震数据同时作为一次波和多次波联合成像的虚拟震源，最后在成像过程中引入立体成像条件通过时间和局部斜率多参数控制实现成像串扰的压制，提高了一次波和多次波联合成像的质量，达到了同时利用一次波和多次波获取更为可靠、丰富的构造信息的目的。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于立体成像条件的一次波和多次波联合成像方法，其特征在于，所述成像方法包括：

获取地表观测点的地震数据，所述地震数据包括一次波和多次波；

利用表层多次波去除方法从所述地震数据中提取震源子波；

根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场；

根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场；

利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场；

根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围；

在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

2.根据权利要求1所述的基于立体成像条件的一次波和多次波联合成像方法，其特征在于，所述利用表层多次波去除方法从所述地震数据中提取震源子波，具体包括：

根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子；

根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

3.根据权利要求2所述的基于立体成像条件的一次波和多次波联合成像方法，其特征在于，所述根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子，具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

4.根据权利要求1所述的基于立体成像条件的一次波和多次波联合成像方法，其特征在于，所述单程波波动方程具体为：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

5.根据权利要求1所述的基于立体成像条件的一次波和多次波联合成像方法，其特征在于，所述在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

6.一种基于立体成像条件的一次波和多次波联合成像系统，其特征在于，所述成像系统包括：

地震数据获取模块，用于获取地表观测点的地震数据，所述地震数据包括一次波和多次波；

震源子波提取模块，用于利用表层多次波去除方法从所述地震数据中提取震源子波；

延拓后震源波场确定模块，用于根据所述地震数据和所述震源子波，采用单程波波动方程沿深度轴进行正向波场延拓，确定深度为z处的波场，记为延拓后震源波场；

延拓后检波波场确定模块，用于根据所述地震数据，采用单程波波动方程沿深度轴进行反向波场延拓，确定深度为z处的波场，记为延拓后检波波场；

波场分解模块，用于利用局部倾斜叠加方法分别对所述延拓后震源波场和所述延拓后检波波场进行分解，得到不同局部斜率的震源波场和不同局部斜率的检波波场；

局部斜率取值范围确定模块，根据设定的速度场，利用平面波解构滤波器确定震源波场和检波波场局部斜率的取值范围；

成像结果确定模块，用于在所述局部斜率的取值范围内，将分解的震源波场和检波波场进行互相关确定成像结果。

7.根据权利要求6所述的基于立体成像条件的一次波和多次波联合成像系统，其特征在于，所述震源子波提取模块，具体包括：

地表加权算子确定单元，根据所述地震数据，基于反馈迭代的表层多次波去除方法确定地表加权算子；

震源子波确定单元，用于根据所述地表加权算子，采用傅里叶反变换确定所述震源子波。

8.根据权利要求7所述的基于立体成像条件的一次波和多次波联合成像系统，其特征在于，所述地表加权算子确定单元，具体公式为：

ΔP⁽ⁿ⁺¹⁾(z₀)＝P(z₀)-A⁽ⁿ⁺¹⁾(w)ΔP⁽ⁿ⁾(z₀)P(z₀)

其中，n表示迭代次数，n＝0,1,2...，z₀表示地表观测点的深度，ΔP⁽ⁿ⁺¹⁾(z₀)表示第n次迭代后一次波数据，ΔP⁽⁰⁾(z₀)＝P(z₀)，P(z₀)表示地震数据，A⁽ⁿ⁺¹⁾(w)表示第n次迭代后地表加权算子，w表示角频率。

9.根据权利要求6所述的基于立体成像条件的一次波和多次波联合成像系统，其特征在于，所述单程波波动方程具体为：

其中，

u(x,z,t)表示波场，v(x,z,w)表示地震波速度，+符号表示上行波场，－符号表示下行波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间，w表示角频率。

10.根据权利要求6所述的基于立体成像条件的一次波和多次波联合成像系统，其特征在于，所述成像结果确定模块，具体成像公式为：

I(x,z)＝∫∫W_S(x,z,p,t)W_D(x,z,p,t)dpdt

其中，p表示空间和时间的局部斜率三维矢量，W_S(x,z,p,t)表示分解的震源波场，W_D(x,z,p,t)表示分解的检波波场，x表示水平方向坐标，z表示垂直方向坐标，t表示地震波传播时间。

Images