CN108802813A - 一种多分量地震资料偏移成像方法及系统 - Google Patents

Abstract

本发明公开一种多分量地震资料偏移成像方法及系统。该偏移成像方法包括：获取观测多分量地震记录；获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；获得每一炮对应的多分量震源波场、多分量检波点波场；对多分量震源波场和多分量检波点波场进行纵横波场分离；利用梯度计算公式获得梯度剖面；构建每一炮对应的下降方向剖面；获得每一炮对应的多分量反偏移模拟波场；获得多分量地震记录增量；根据多分量地震记录增量和下降方向剖面，确定最优化步长；根据最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。采用本发明的方法或系统，可获得能够直接反映地下介质的PP、PS、SP和SS反射系数信息的偏移剖面。

Description

一种多分量地震资料偏移成像方法及系统

技术领域

本发明涉及勘探地震领域，特别是涉及一种多分量地震资料偏移成像方法及系统。

背景技术

多分量地震偏移成像技术是地震勘探在油气田开发应用中的关键技术之一，高精度的多分量地震偏移成像方法能够获得高质量的地下介质构造形态，为后续的多分量地震勘探解释提供准确的基础资料。随着多分量地震偏移技术水平的不断提高以及高性能计算技术的迅速发展，多分量地震偏移成像方法逐渐由射线类叠前时间偏移转向波动方程类叠前深度偏移。为了获得具有明确物理意义的偏移剖面，充分发挥横波的优势，业内普遍采用基于纵横波分离的成像条件进行成像。在空间每个成像点上，该成像条件能够获得四个物理意义明确的偏移剖面，分别是PP波、PS波、SP波和SS波剖面，该偏移剖面不受纵横波串扰假象影响，而且能够直接表征地下介质的PP、PS、SP和SS四个反射系数信息，因此能直接用于后续的地下介质属性解释。

实际中，由于诸多问题的存在，比如实际地震数据并不完整且存在噪音、采集孔径有限、介质速度变化剧烈、偏移算子近似误差大等，这些因素导致多分量资料叠前深度偏移剖面存在采集脚印、分辨率低、假象严重、振幅失衡等问题，使偏移剖面难以用于实际生产，增加后续数据处理和解释的难度。多分量地震资料最小二乘偏移通过引入反演理论，能够在一定程度上解决常规偏移剖面分辨率低、假象严重、振幅失衡等问题。但是与常规叠前深度偏移理论不同，最小二乘偏移是基于介质扰动理论进行的公式推导，其偏移剖面与地下介质参数的相对变化有关，该偏移剖面与常规多分量地震资料叠前深度偏移剖面的物理意义完全不同，并不能直接表征地下介质的PP、PS、SP和SS四个反射系数信息，因此该偏移剖面难以直接用于后续的地震资料处理及解释，影响储层预测、流体识别等环节的精度。因此对于实际资料，现有的多分量地震资料最小二乘偏移方法并不能获得具有明确物理意义的高质量偏移剖面。为此，必须建立一套新的、可获得高质量的(高分辨率、高信噪比、振幅均衡)、能够直接表征地下介质的PP、PS、SP和SS四个反射系数信息的多分量地震资料叠前深度偏移方法。

发明内容

本发明的目的是提供一种多分量地震资料偏移成像方法及系统，能够直接表征地下介质的PP、PS、SP和SS四个反射系数信息，以提高偏移剖面的精度和分辨率。

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

一种多分量地震资料偏移成像方法，所述偏移成像方法包括：

获取待进行偏移成像的观测多分量地震记录；

获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场；

获取每一炮对应的多分量检波点波场；

根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场；

根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面；所述梯度剖面为每一炮的单炮梯度剖面的集合；

根据最优化反演算法，构建每一炮对应的下降方向剖面；所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面；

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场；

根据所述多分量反偏移模拟波场，获得多分量地震记录增量；

根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长；

根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。

可选的，所述根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场，具体包括：

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质弹性波方程，对该炮点波场的正向延拓，获得该炮的每一个时刻的多分量震源波场，从而获得每一炮对应的多分量震源波场。

可选的，所述获取每一炮对应的多分量检波点波场，具体包括：

对于每一炮，获取前一次迭代获得的模拟多分量地震记录其中i表示当前的迭代次数；

根据所述观测多分量地震记录和所述模拟多分量地震记录，利用公式确定多分量地震记录残差；其中d^obs为所述观测多分量地震记录，和分别表示观测多分量地震记录d^obs在笛卡尔坐标系x，y和z三个方向的分量；T表示转置；为前一次迭代获得的模拟多分量地震记录；δd_i为多分量地震记录残差，δd_i＝(δd_x,δd_y,δd_z)^T；

根据各向同性介质弹性波动方程，对所述多分量地震记录残差δd_i进行波场逆时延拓，获得该炮的每一个时刻的多分量检波点波场。

可选的，所述根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场，具体包括：

对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；

根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

根据所述下降方向剖面的PP波和SP波分量剖面、所述纵波背景波场、所述纵波偏移速度及所述偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

根据所述下降方向剖面的PS波和SS波分量剖面、所述横波背景波场、所述横波偏移速度及所述偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

在相同的时刻，利用纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场。

可选的，所述根据所述多分量反偏移模拟波场，获得多分量地震记录增量，具体包括：

对所述多分量反偏移模拟波场进行记录采样，获得每一炮对应的多分量地震记录增量。

可选的，根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长，具体包括：

利用公式确定当前迭代的最优化步长；其中，L表示反偏移算子，r_i ^s为第s炮的下降方向剖面，r_i ^s＝(r_i ^PP,s,r_i ^PS,s,r_i ^SP,s,r_i ^SS,s)；r_i ^PP,s，r_i ^PS,s，r_i ^SP,s和r_i ^SS,s分别表示第s炮单炮下降方向剖面对应的PP波，PS波，SP波和SS波分量剖面；为第s炮的多分量地震记录增量，Δd_x，Δd_y和Δd_z分别表示第s炮的多分量地震记录增量在笛卡尔坐标系x，y和z三个方向的分量，α_i为当前迭代的最优化步长。

可选的，所述根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长，之后还包括：

利用公式确定当前第i次迭代后的模拟地震记录；其中，为当前第i次迭代后的模拟地震记录，为前一次迭代后的模拟地震记录，α_i为当前迭代的最优化步长，Δd_i为当前迭代的多分量地震记录增量。

可选的，所述根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面，之后还包括：

根据公式计算当前第i次迭代后的目标函数值；其中，δd_i为当前第i次迭代获得的多分量地震记录残差，Obj_i为当前第i次迭代后的目标函数值；

判断所述目标函数值是否满足收敛标准；

当所述目标函数值满足收敛标准时，将所述偏移剖面确定为最终的偏移剖面；

当所述目标函数不满足收敛标准时，返回“获取每一炮对应的多分量检波点波场”，进入下一次迭代。

一种多分量地震资料偏移成像系统，所述偏移成像系统包括：

观测多分量地震记录获取模块，用于获取待进行偏移成像的观测多分量地震记录；

初始参数获取模块，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；

多分量震源波场获取模块，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场；

多分量检波点波场获取模块，用于获取每一炮对应的多分量检波点波场；

纵横波场分离模块，用于根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场；

梯度剖面获取模块，用于根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面；所述梯度剖面为每一炮的单炮梯度剖面的集合；

下降方向剖面构建模块，用于根据最优化反演算法，构建每一炮对应的下降方向剖面；所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面；

多分量反偏移模拟波场获取模块，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场；

多分量地震记录增量获取模块，用于根据所述多分量反偏移模拟波场，获得多分量地震记录增量；

最优化步长确定模块，用于根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长；

偏移剖面确定模块，用于根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。

可选的，所述多分量反偏移模拟波场获取模块，具体包括：

震源子波设置单元，用于对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；

纵波背景波场获取单元，用于根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

横波背景波场获取单元，用于根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

纵波反偏移算子方程构建单元，用于根据所述下降方向剖面的PP波和SP波分量剖面、所述纵波背景波场、所述纵波偏移速度及所述偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

横波反偏移算子方程构建单元，用于根据所述下降方向剖面的PS波和SS波分量剖面、所述横波背景波场、所述横波偏移速度及所述偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

纵波反偏移模拟波场和横波反偏移模拟波场获取单元，用于基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

纵横波场合成单元，用于在相同的时刻，利用纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场。

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

1)本发明方法为一种多分量地震资料反演偏移方法，与常规叠前深度偏移相比，本发明可以获得高精度、高分辨率、高信噪比、振幅保真的偏移剖面；

2)本发明以地下介质的PP、PS、SP和SS四个反射系数信息为反演目标，通过精确的多分量波场分解及合成算子的应用，可以获得物理意义明确的四个偏移剖面，该偏移剖面直接反映了地下介质的四个反射系数信息，可以直接用于后续的资料解释，极大提高了储层预测及流体识别的精度；

3)本发明可以广泛用于油气勘探领域中，特别是对于复杂构造深部的成像效果更加明显。

附图说明

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

图1为本发明多分量地震资料偏移成像方法的流程示意图；

图2为本发明多分量地震资料偏移成像系统的结构示意图；

图3为本发明实施例1利用有限差分方法模拟的X分量；

图4为本发明实施例1利用有限差分方法模拟的Z分量；

图5为本发明实施例1利用本发明方法分离的P分量；

图6为本发明实施例1中利用本发明方法分离的S分量；

图7为本发明实施例1中利用本发明方法合成的X分量；

图8为本发明实施例1中利用本发明方法合成的Z分量；

图9为本发明实施例2中纵波速度模型；

图10为本发明实施例2中横波速度模型；

图11为本发明实施例2中利用逆时偏移方法所得的PP偏移剖面；

图12为本发明实施例2中利用逆时偏移方法所得的PS偏移剖面；

图13为本发明实施例2中利用本发明方法所得的PP偏移剖面；

图14为本发明实施例2中利用本发明方法所得的PS偏移剖面；

图15为本发明实施例3中纵波速度模型；

图16为本发明实施例3中横波速度模型；

图17为本发明实施例3中利用逆时偏移方法所得的PP偏移剖面；

图18为本发明实施例3中利用逆时偏移方法所得的PS偏移剖面；

图19为本发明实施例3中利用本发明方法所得的PP偏移剖面；

图20为本发明实施例3中利用本发明方法所得的PS偏移剖面。

具体实施方式

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

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

图1为本发明多分量地震资料偏移成像方法的流程示意图。如图1所示，所述偏移成像方法包括以下步骤：

步骤100：获取待进行偏移成像的观测多分量地震记录。根据实际需求，从观测地震记录库中筛选出待进行偏移成像的观测多分量地震记录。待进行偏移成像的观测多分量地震记录用表示，和分别表示观测多分量地震记录d^obs在笛卡尔坐标系x，y和z三个方向的分量，T表示转置。

步骤200：获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数。

步骤300：根据观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场。

具体的，针对每一炮，基于获取的观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波，利用上述纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质弹性波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的多分量震源波场，从而获得每一炮对应的多分量震源波场。

其中所述的各向同性介质弹性波方程为：

所述方程(1)中，V＝(V_x,V_y,V_z)^T表示质点速度场矢量；σ＝(σ_xx,σ_yy,σ_zz,σ_xy,σ_yz,σ_zx)^T表示应力场矢量；t表示时间；x，y和z分别表示笛卡尔坐标的x，y和z三个方向；λ和μ表示介质拉梅常数，ρ表示介质密度；介质纵波偏移速度v_p和横波偏移速度速度v_s可由下式子计算获得：

步骤400：获取每一炮对应的多分量检波点波场。具体过程如下：

设当前迭代次数为i，对于每一炮，获取前一次迭代获得的模拟多分量地震记录其中i表示当前的迭代次数；

根据所述观测多分量地震记录和所述模拟多分量地震记录，利用公式确定多分量地震记录残差；其中d^obs为所述观测多分量地震记录，和分别表示观测多分量地震记录d^obs在笛卡尔坐标系x，y和z三个方向的分量；T表示转置；为前一次迭代获得的模拟多分量地震记录；和分别表示模拟多分量地震记录d^cal在笛卡尔坐标系x，y和z三个方向的分量；δd_i为多分量地震记录残差，δd_i＝(δd_x,δd_y,δd_z)^T；δd_x，δd_y和δd_z分别表示多分量地震记录残差δd_i在笛卡尔坐标系x，y和z三个方向的分量。

根据各向同性介质弹性波动方程，对该炮的多分量地震记录残差δd_i进行波场逆时延拓，获得该炮的每一个时刻的多分量检波点波场。

步骤500：根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场。在相同的时刻，应用精确的多分量波场纵横波分离方程进行纵横波场分离，因此，可以获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场四种波场。

步骤600：根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面。所述梯度剖面为每一炮的单炮梯度剖面的集合。

对每一炮的所述纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场应用梯度计算公式，进而获得每一炮对应的单炮梯度剖面所有的单炮剖面的集合构成了本次迭代的梯度剖面，也即第i次迭代的梯度剖面和分别表示第s炮单炮梯度剖面对应的PP波，PS波，SP波和SS波分量剖面；Ns表示偏移总炮数。具体过程如下：

精确的多分量波场纵横波分离方程为：

所述方程(3)中，P和S分别表示分离的纵波场和横波场；a＝(a_x,a_y,a_z)^T表示横波方向矢量，其计算式为：

所述方程(4)中，E_S和E_R分别表示震源波场的坡印廷矢量和检波点波场的坡印廷矢量，所述的震源波场的坡印廷矢量和检波点波场的坡印廷矢量均由下式计算：

基于方程(5)，利用震源波场的质点速度场矢量和应力场矢量可求出震源波场的波印廷矢量E_S，利用检波点波场的质点速度场矢量和应力场矢量可求出检波点波场的波印廷矢量E_R；

所述应用精确的多分量波场纵横波分离方程对所述多分量震源波场和所述的多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，其特征为：在相同时刻，分别对震源波场和检波点波场应用方程(3)，即可获得该时刻的纵波震源波场P_S、横波震源波场S_S、纵波检波点波场P_R和横波检波点波场S_R；

对于第s炮，应用梯度计算公式获得该炮的单炮梯度剖面，具体为：

步骤700：根据最优化反演算法，构建每一炮对应的下降方向剖面。所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面。构建的每一炮对应的下降方向剖面为r_i ^s＝(r_i ^PP,s,r_i ^PS,s,r_i ^SP,s,r_i ^SS,s)；所有的单炮下降方向剖面的集合构成了第i次迭代的下降方向剖面r_i＝(r_i ¹,r_i ²,...,r_i ^Ns)^T；其中，r_i ^PP,s，r_i ^PS,s，r_i ^SP,s和r_i ^SS,s分别表示第s炮单炮下降方向剖面对应的PP波，PS波，SP波和SS波分量剖面。

步骤800：根据观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场。

对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

同理，根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

根据第i次迭代获得的下降方向剖面r_i的PP波和SP波分量剖面、该炮的纵波背景波场、纵波偏移速度及偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

根据第i次迭代获得的下降方向剖面r_i的PS波和SS波分量剖面、该炮的横波背景波场、横波偏移速度及偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

在相同的时刻，利用精确的纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场δV＝(δV_x,δV_y,δV_z)^T。

各向同性介质声波方程为：

所述方程(7)中，U表示波场；

建立纵波反偏移算子方程，具体为：

所述方程(8)中，U_P表示纵波背景波场，δU_P表示纵波反偏移波场，δV_P＝(δV_Px,δV_Py,δV_Pz)^T表示纵波反偏移质点速度场矢量；所述纵波虚拟震源，具体为：

同理，建立横波反偏移算子方程，具体为：

所述方程(10)中，U_S表示横波背景波场，δU_S表示横波反偏移波场，δV_S＝(δV_Sx,δV_Sy,δV_Sz)^T表示横波反偏移质点速度场矢量；所述横波虚拟震源，具体为：

所述的在相同的时刻，利用精确的纵横波场合成方程对所得的纵波和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场δV＝(δV_x,δV_y,δV_z)^T过程中，所述的精确的纵横波场合成方程具体为：

所述方程(12)中，和为δU_S在笛卡尔坐标的x，y和z三个方向的分量，其计算方程具体为：

所述方程(13)中，cosθ_x，cosθ_y和cosθ_z为横波反偏移波场δP_S在笛卡尔坐标的x，y和z三个方向的方向余弦，其计算方程具体为：

所述方程(14)中，和为所述纵波背景波场的坡印廷矢量E_B在笛卡尔坐标的x，y和z三个方向的分量；和为所述横波反偏移波场的坡印廷矢量E_D在笛卡尔坐标的x，y和z三个方向的分量；其中纵波背景波场的坡印廷矢量E_B和横波反偏移波场的坡印廷矢量E_D均由下式计算：

步骤900：根据多分量反偏移模拟波场，获得多分量地震记录增量。对该炮的多分量反偏移模拟波场进行记录采样，获得第s炮的多分量地震记录增量进而获得每一炮对应的多分量地震记录增量，所有炮的多分量地震记录增量构成了第i次迭代的多分量地震记录增量其中，Δd_x，Δd_y和Δd_z分别表示第s炮的多分量地震记录增量在笛卡尔坐标系x，y和z三个方向的分量。

步骤1000：根据多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长。具体的，利用公式确定当前迭代的最优化步长；其中，L表示反偏移算子，r_i ^s为第s炮的下降方向剖面，r_i ^s＝(r_i ^PP,s,r_i ^PS,s,r_i ^SP,s,r_i ^SS,s)；r_i ^PP,s，r_i ^PS,s，r_i ^SP,s和r_i ^SS,s分别表示第s炮单炮下降方向剖面对应的PP波，PS波，SP波和SS波分量剖面；为第s炮的多分量地震记录增量，Δd_x，Δd_y和Δd_z分别表示第s炮的多分量地震记录增量在笛卡尔坐标系x，y和z三个方向的分量，α_i为当前迭代的最优化步长。

确定最优化步长之后，利用公式确定当前第i次迭代后的模拟地震记录其中，为当前第i次迭代后的模拟地震记录，为前一次迭代后的模拟地震记录，α_i为当前迭代的最优化步长，Δd_i为当前迭代的多分量地震记录增量。

步骤1100：根据最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。利用最优化步长α_i及下降方向剖面r_i，更新第i次迭代的偏移剖面m_i＝m_i-1+α_ir_i；其中，偏移剖面m_i包含和四个分量剖面。

第i次迭代完后，计算第i次迭代的目标函数值Obj_i，判断当前迭代是否满足收敛标准，如果满足则输出最新的偏移剖面为最终的偏移剖面m_i；否则重复步骤400-步骤1100，直至获得最终偏移剖面。

所述计算第i次迭代的目标泛函值Obj_i具体为：δd_i为当前第i次迭代获得的多分量地震记录残差，Obj_i为当前第i次迭代后的目标函数值。

收敛标准具体为：其中，Theshold表示迭代停止的阈值标准，通常选取1.0e-6。

本发明包括：建立精确的多分量波场分解与合成算子；确定多分量波场反偏移算子；确定多分量波场偏移算子；建立梯度计算式；建立目标函数和反演流程；确定最优化步长计算方法。本发明将反演思想与多分量资料叠前深度偏移成像方法相结合，通过建立精确的多分量波场分解与合成算子，可获得高精度、高分辨率、高信噪比、振幅保真的多分量地震资料叠前深度偏移剖面，而且偏移剖面能够直接反映地下介质的PP、PS、SP和SS反射系数信息。该方法可应用到各种复杂地质条件的地震资料偏移中，成像精度和分辨率高，剖面物理意义明确，可直接用于后期地质解译。

本发明还提供一种多分量地震资料偏移成像系统，图2为本发明多分量地震资料偏移成像系统的结构示意图。如图2所示，所述偏移成像系统包括：

观测多分量地震记录获取模块201，用于获取待进行偏移成像的观测多分量地震记录；

初始参数获取模块202，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；

多分量震源波场获取模块203，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场；

多分量检波点波场获取模块204，用于获取每一炮对应的多分量检波点波场；

纵横波场分离模块205，用于根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场；

梯度剖面获取模块206，用于根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面；所述梯度剖面为每一炮的单炮梯度剖面的集合；

下降方向剖面构建模块207，用于根据最优化反演算法，构建每一炮对应的下降方向剖面；所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面；

多分量反偏移模拟波场获取模块208，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场；

多分量地震记录增量获取模块209，用于根据所述多分量反偏移模拟波场，获得多分量地震记录增量；

最优化步长确定模块2010，用于根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长；

偏移剖面确定模块2011，用于根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。

所述多分量反偏移模拟波场获取模块208，具体包括：

震源子波设置单元，用于对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；

纵波背景波场获取单元，用于根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

横波背景波场获取单元，用于根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

纵波反偏移算子方程构建单元，用于根据所述下降方向剖面的PP波和SP波分量剖面、所述纵波背景波场、所述纵波偏移速度及所述偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

横波反偏移算子方程构建单元，用于根据所述下降方向剖面的PS波和SS波分量剖面、所述横波背景波场、所述横波偏移速度及所述偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

纵波反偏移模拟波场和横波反偏移模拟波场获取单元，用于基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

纵横波场合成单元，用于在相同的时刻，利用纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场。

为进一步说明本发明的可行性和有效性，下面举三个实施例：

实施例1：

图3-图8是均匀介质模型350ms时刻地震波场快照。其中，图3为本发明实施例1利用有限差分方法模拟的X分量；图4为本发明实施例1利用有限差分方法模拟的Z分量；图5为本发明实施例1利用本发明方法分离的P分量；图6为本发明实施例1中利用本发明方法分离的S分量；图7为本发明实施例1中利用本发明方法合成的X分量；图8为本发明实施例1中利用本发明方法合成的Z分量。将图3和4与图5和6相比可以看出，纵横波获得了有效的分离，证明了本发明提出的多分量波场分离方法是有效的。将图3和4与图7和8相比可以看出，利用纵横波场有效地合成了多分量波场，证明了本发明提出的多分量波场合成方法是有效的。

实施例2：

图9-图10是二维洼陷偏移速度模型图，图9为本发明实施例2中纵波速度模型；图10为本发明实施例2中横波速度模型。在此模型上设置30个爆炸震源，震源子波设定为雷克子波，主频为20赫兹，起始震源点的位于(150m，100m)处，炮间隔为100m。采用中间放炮两边接收观测系统，单边最大偏移距2200m，最小偏移距为150m，道间距为10m。

图11-图14是图9-10所示二维洼陷模型的多炮叠加偏移剖面。其中，图11为本发明实施例2中利用逆时偏移方法所得的PP偏移剖面；图12为本发明实施例2中利用逆时偏移方法所得的PS偏移剖面；图13为本发明实施例2中利用本发明方法所得的PP偏移剖面；图14为本发明实施例2中利用本发明方法所得的PS偏移剖面。从图11和12中可以看出，剖面存在比较明显的噪声，剖面的分辨率较低，振幅不均衡。从图13和图14中可以看出，本发明方法所得的偏移剖面具有很高的精度、分辨率和信噪比，而且振幅是均衡性很好，也证明了本发明方法的可行性及有效性。

实施例3：

图15-图16是Marmousi-2偏移速度模型，图15为本发明实施例3中纵波速度模型；图16为本发明实施例3中横波速度模型。该模型是验证各种偏移方法成像效果的国际标准模型之一。在此模型上设置87个爆炸震源，震源子波设定为雷克子波，主频为10赫兹，起始震源点的位于(300m，300m)处，炮间隔为300m。采用中间放炮两边接收观测系统，单边最大偏移距5200m，道间距为20m。

图17-图20是图15-16所示Marmousi-2模型的多炮叠加偏移剖面。其中，图17为本发明实施例3中利用逆时偏移方法所得的PP偏移剖面；图18为本发明实施例3中利用逆时偏移方法所得的PS偏移剖面；图19为本发明实施例3中利用本发明方法所得的PP偏移剖面；图20为本发明实施例3中利用本发明方法所得的PS偏移剖面。由图17-图20可以看出，传统方法的成像剖面振幅严重失衡，深部剖面相对振幅并不保真。而本发明方法所得成像剖面，效果更好，分辨率和精度更高，信噪比也更好，振幅均衡性更佳。

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

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

Claims (10)

1.一种多分量地震资料偏移成像方法，其特征在于，所述偏移成像方法包括：

获取待进行偏移成像的观测多分量地震记录；

获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场；

获取每一炮对应的多分量检波点波场；

根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场；

根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面；所述梯度剖面为每一炮的单炮梯度剖面的集合；

根据最优化反演算法，构建每一炮对应的下降方向剖面；所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面；

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场；

根据所述多分量反偏移模拟波场，获得多分量地震记录增量；

根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长；

根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。

2.根据权利要求1所述的偏移成像方法，其特征在于，所述根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场，具体包括：

根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质弹性波方程，对该炮点波场的正向延拓，获得该炮的每一个时刻的多分量震源波场，从而获得每一炮对应的多分量震源波场。

3.根据权利要求1所述的偏移成像方法，其特征在于，所述获取每一炮对应的多分量检波点波场，具体包括：

对于每一炮，获取前一次迭代获得的模拟多分量地震记录其中i表示当前的迭代次数；

根据所述观测多分量地震记录和所述模拟多分量地震记录，利用公式确定多分量地震记录残差；其中d^obs为所述观测多分量地震记录， 和分别表示观测多分量地震记录d^obs在笛卡尔坐标系x，y和z三个方向的分量；T表示转置；为前一次迭代获得的模拟多分量地震记录；δd_i为多分量地震记录残差，δd_i＝(δd_x,δd_y,δd_z)^T；

根据各向同性介质弹性波动方程，对所述多分量地震记录残差δd_i进行波场逆时延拓，获得该炮的每一个时刻的多分量检波点波场。

4.根据权利要求1所述的偏移成像方法，其特征在于，所述根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场，具体包括：

对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；

根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

根据所述下降方向剖面的PP波和SP波分量剖面、所述纵波背景波场、所述纵波偏移速度及所述偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

根据所述下降方向剖面的PS波和SS波分量剖面、所述横波背景波场、所述横波偏移速度及所述偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

在相同的时刻，利用纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场。

5.根据权利要求1所述的偏移成像方法，其特征在于，所述根据所述多分量反偏移模拟波场，获得多分量地震记录增量，具体包括：

对所述多分量反偏移模拟波场进行记录采样，获得每一炮对应的多分量地震记录增量。

6.根据权利要求1所述的偏移成像方法，其特征在于，根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长，具体包括：

利用公式确定当前迭代的最优化步长；其中，L表示反偏移算子，r_i ^s为第s炮的下降方向剖面，r_i ^s＝(r_i ^PP,s,r_i ^PS,s,r_i ^SP,s,r_i ^SS,s)；r_i ^PP,s，r_i ^PS,s，r_i ^SP,s和r_i ^SS,s分别表示第s炮单炮下降方向剖面对应的PP波，PS波，SP波和SS波分量剖面；为第s炮的多分量地震记录增量，Δd_x，Δd_y和Δd_z分别表示第s炮的多分量地震记录增量在笛卡尔坐标系x，y和z三个方向的分量，α_i为当前迭代的最优化步长。

7.根据权利要求1所述的偏移成像方法，其特征在于，所述根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长，之后还包括：

利用公式确定当前第i次迭代后的模拟地震记录；其中，为当前第i次迭代后的模拟地震记录，为前一次迭代后的模拟地震记录，α_i为当前迭代的最优化步长，Δd_i为当前迭代的多分量地震记录增量。

8.根据权利要求3所述的偏移成像方法，其特征在于，所述根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面，之后还包括：

根据公式计算当前第i次迭代后的目标函数值；其中，δd_i为当前第i次迭代获得的多分量地震记录残差，Obj_i为当前第i次迭代后的目标函数值；

判断所述目标函数值是否满足收敛标准；

当所述目标函数值满足收敛标准时，将所述偏移剖面确定为最终的偏移剖面；

当所述目标函数不满足收敛标准时，返回“获取每一炮对应的多分量检波点波场”，进入下一次迭代。

9.一种多分量地震资料偏移成像系统，其特征在于，所述偏移成像系统包括：

观测多分量地震记录获取模块，用于获取待进行偏移成像的观测多分量地震记录；

初始参数获取模块，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数；

多分量震源波场获取模块，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，获得每一炮对应的多分量震源波场；

多分量检波点波场获取模块，用于获取每一炮对应的多分量检波点波场；

纵横波场分离模块，用于根据多分量波场纵横波分离方程对每一炮的多分量震源波场和多分量检波点波场进行纵横波场分离，获得相同时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场；

梯度剖面获取模块，用于根据每一炮同一时刻的纵波震源波场、横波震源波场、纵波检波点波场和横波检波点波场，利用梯度计算公式获得梯度剖面；所述梯度剖面为每一炮的单炮梯度剖面的集合；

下降方向剖面构建模块，用于根据最优化反演算法，构建每一炮对应的下降方向剖面；所述下降方向剖面包括单炮下降方向剖面对应的PP波、PS波、SP波和SS波分量剖面；

多分量反偏移模拟波场获取模块，用于根据所述观测系统参数、纵波偏移速度、横波偏移速度、偏移密度模型及偏移参数，利用多分量波场合成算子，获得每一炮对应的多分量反偏移模拟波场；

多分量地震记录增量获取模块，用于根据所述多分量反偏移模拟波场，获得多分量地震记录增量；

最优化步长确定模块，用于根据所述多分量地震记录增量和每一炮对应的下降方向剖面，确定最优化步长；

偏移剖面确定模块，用于根据所述最优化步长和每一炮对应的下降方向剖面，确定偏移剖面。

10.根据权利要求9所述的偏移成像系统，其特征在于，所述多分量反偏移模拟波场获取模块，具体包括：

震源子波设置单元，用于对于每一炮，根据所述观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波；

纵波背景波场获取单元，用于根据所述纵波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的纵波背景波场；

横波背景波场获取单元，用于根据所述横波偏移速度、偏移密度模型及偏移参数，基于交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得该炮的每一个时刻的横波背景波场；

纵波反偏移算子方程构建单元，用于根据所述下降方向剖面的PP波和SP波分量剖面、所述纵波背景波场、所述纵波偏移速度及所述偏移密度模型，构建该炮的纵波虚拟震源，建立纵波反偏移算子方程；

横波反偏移算子方程构建单元，用于根据所述下降方向剖面的PS波和SS波分量剖面、所述横波背景波场、所述横波偏移速度及所述偏移密度模型，构建该炮的横波虚拟震源，建立横波反偏移算子方程；

纵波反偏移模拟波场和横波反偏移模拟波场获取单元，用于基于所述的纵波反偏移算子方程和横波反偏移算子方程，利用所述纵波偏移速度、所述横波偏移速度及所述偏移密度模型对该炮点进行波场正向延拓，获得该炮的每一个时刻的纵波反偏移模拟波场和横波反偏移模拟波场；

纵横波场合成单元，用于在相同的时刻，利用纵横波场合成方程对纵波反偏移模拟波场和横波反偏移模拟波场进行纵横波场合成，获得该炮的每一个时刻的多分量反偏移模拟波场。