一种三维起伏地表地震资料偏移成像方法、装置及系统

摘要

本申请公开了一种三维起伏地表地震资料偏移成像方法、装置及系统，其中，方法包括：获取三维地形偏导数；获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；根据三维地形偏导数、炮点和接收点分别到成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，根据计算结果、每个炮点和每个接收点分别到地下各成像点的走时、每个炮点和每个接收点分别到地下各成像点的振幅对三维起伏地表地震资料偏移成像。

说明书

技术领域

本申请涉及地震资料处理技术领域，特别涉及一种三维起伏地表地震资料偏移成像方法、装置及系统。

背景技术

地震勘探是一种利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应来推断地下构造形态和岩层性质的方法。地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘察、区域地质研究和地壳研究等方面，地震勘探也得到广泛应用。地震勘探在分层的详细程度和勘查的精度上，都优于其他地球物理勘探方法。

地震勘探通常分为折射波地震勘探和反射波地震勘探，目前在石油和天然气资源勘查领域，主要采用反射波地震勘探。地震勘探过程由地震数据采集、地震数据处理和地震资料解释三个阶段组成。第一阶段是地震数据采集，其主要任务是在油气勘探有利区带，布置二维或三维测线，使用炸药震源或可控震源激发地震波，在测线上等间距布置多个检波器来接收地震波信号，以等时间间隔离散采样地震数据，并以数字形式记录在磁盘上。第二阶段是地震数据处理，其主要任务是利用计算机和软件，加工处理野外采集的地震数据，将地震数据变成能够反映地下构造的地震剖面图和能够反映地下岩性变化的地震波振幅、频率和传播速度等信息图。地震数据处理主要包括噪声压制、反褶积、水平叠加和偏移成像等关键环节。第三个阶段是地震资料解释，其主要任务是分析地震数据处理所提供的资料，确定地下岩层的产状和构造关系，找出有利的含油气目标，提供钻井位置。地震资料解释主要包括构造解释、储层预测和流体识别等。

地震资料偏移成像是地震资料处理的关键环节之一，也是提交处理成果的最后环节。目前工业界常用的深度偏移成像技术是克希霍夫叠前深度偏移成像，其中研究热点之一是真振幅深度偏移成像。1985年Beylkin提出了一种逆散射间断面成像方法，导出了一个关键的变化率函数，它本质上是一种变量之间变换的雅可比，在真振幅深度偏移成像或逆散射成像中起到关键作用，后人称之为Beylkin行列式。Cohen等对Beylkin行列式进行了详细的诠释。Bleistein在研究逆散射成像过程中，将李普曼—许温格散射方程看成是一个傅里叶变换，将波速扰动和散射波场联系起来，再猜想一个傅里叶反变换就可以由散射波场计算波速扰动量。在此过程中，需要求解待定系数，结合δ函数的三维积分公式，用波数矢量代替圆频率和地面炮点接收点坐标变量。在积分变量变换过程中，就需要计算Beylkin行列式。所以说Beylkin行列式是逆散射成像的核心。Hubral和Schleicher等先后分别为叠后地震资料和叠前地震资料介绍了真振幅偏移，认为真振幅偏移的波场应该体现反射系数随入射角变化、消除一次反射波的几何扩散。他们用稳相法对绕射叠加积分进行渐近评估，得到的权函数与Beylkin行列式相当。Beylkin行列式计算一般都是建立在运动学射线追踪和动力学射线追踪基础之上。然而，现有文献介绍的Beylkin行列式计算方法都是针对水平地表的，至今未见到三维起伏地表Beylkin行列式的计算方法，常规的水平地表Beylkin行列式无法满足三维起伏地表地震资料真振幅深度偏移成像或逆散射成像需求。

发明内容

为了满足三维起伏地表地震资料真振幅深度偏移成像或逆散射成像需求，本申请实施例提出一种三维起伏地表地震资料偏移成像方法、装置及系统。

为实现上述目的，本申请实施方式提供一种三维起伏地表地震资料偏移成像方法，包括：

获取三维地形偏导数；

获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；

获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值；

根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

优选地，获取三维地形偏导数的步骤包括：

从地震叠前数据道头字中提取地表高程数据；

对地表高程数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，获取三维地形偏导数的步骤包括：

根据接收点高程数据获取地表高程测量数据；

对地表高程测量数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时通过运动学射线追踪方法获得。

优选地，每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵通过动力学射线追踪方法获得。

优选地，对三维起伏地表Beylkin行列式进行计算的步骤包括：

根据每个炮点和每个接收点分别到地下各成像点的射线方向确定炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系和几何扩散矩阵、炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角对三维起伏地表Beylkin行列式进行计算。

为实现上述目的，本申请实施方式提供一种三维起伏地表地震资料偏移成像装置，包括：

三维地形偏导数单元，用于获取三维地形偏导数；

中心射线坐标系、射线方向和走时单元，用于获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；

振幅和几何扩散矩阵单元，用于获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；

计算单元，用于根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值；

成像单元，用于根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

优选地，所述三维地形偏导数单元包括：

提取模块，用于从地震叠前数据道头字中提取地表高程数据；

第一二阶差分模块，用于对地表高程数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，所述三维地形偏导数单元包括：

地表高程测量数据获取模块，用于根据接收点高程数据获取地表高程测量数据；

第二二阶差分模块，用于对地表高程测量数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，所述中心射线坐标系、射线方向和走时单元通过运动学射线追踪方法获得每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时。

优选地，所述振幅和几何扩散矩阵单元通过动力学射线追踪方法获得每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵。

优选地，所述计算单元包括：

下倾角和方位角获取模块，用于根据每个炮点和每个接收点分别到地下各成像点的射线方向确定炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角；

三维起伏地表Beylkin行列式计算模块，用于根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系和几何扩散矩阵、炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角对三维起伏地表Beylkin行列式进行计算。

为实现上述目的，本申请实施方式提供一种三维起伏地表地震资料偏移成像系统，所述系统包括：存储器和处理器，所述存储器中存储计算机程序，所述计算机程序被所述处理器执行时，实现以下功能：

获取三维地形偏导数；

获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；

获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值；

根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

优选地，获取三维地形偏导数，所述计算机程序被所述处理器执行时，实现以下功能：

从地震叠前数据道头字中提取地表高程数据；

对地表高程数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，获取三维地形偏导数，所述计算机程序被所述处理器执行时，实现以下功能：

根据接收点高程数据获取地表高程测量数据；

对地表高程测量数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

优选地，对三维起伏地表Beylkin行列式进行计算，所述计算机程序被所述处理器执行时，实现以下功能：

根据每个炮点和每个接收点分别到地下各成像点的射线方向确定炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系和几何扩散矩阵、炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角对三维起伏地表Beylkin行列式进行计算。

由上可见，与现有技术相比较，本申请提供的技术方案首先获取三维地形偏导数，再用运动学射线追踪计算每个炮点或每个接收点到地下各成像点的射线方向和走时，并用动力学射线追踪方法计算每个炮点或每个接收点到地下各成像点的振幅和几何扩散矩阵，最后计算三维起伏地表Beylkin行列式，根据计算结果应用于三维起伏地表真振幅深度偏移成像或逆散射成像。与背景技术中提到的适用于水平地表Beylkin行列式计算方法相比，本技术方案考虑到地形的起伏变化，在三维起伏地表Beylkin行列式的计算中引入了地形偏导数，为真振幅深度偏移成像或逆散射成像提供合理的权函数，使之适应于起伏地表区地震资料成像。可以说在几乎没有增加计算成本的基础上，为起伏地表地震资料成像提供了真正的而非近似的Beylkin行列式。

附图说明

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

图1为本申请实施例提出一种三维起伏地表地震资料偏移成像方法流程图；

图2为本申请实施例中常速背景介质中走时图；

图3为本申请实施例中常速背景介质中射线的下倾角图；

图4为本申请实施例中常速背景介质中射线的方位角图；

图5为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₁图；

图6为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₁图；

图7为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₂图；

图8为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₂图；

图9为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₃图；

图10为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₃图；

图11为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₁图；

图12为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₁图；

图13为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₂图；

图14为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₂图；

图15为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₃图；

图16为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₃图；

图17为本申请实施例提出的一种三维起伏地表地震资料偏移成像装置功能框图；

图18为本申请实施例提出的一种三维起伏地表地震资料偏移成像系统示意图。

具体实施方式

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

为了解决背景技术中提到的问题，为保证起伏地表三维地震资料真振幅深度偏移成像或逆散射成像的实施，本申请实施例提出一种三维起伏地表地震资料偏移成像方法，如图1所示。包括：

步骤101)：获取三维地形偏导数。

在本实施例中，从地震叠前数据道头字中提取地表高程数据，计算三维起伏地表地形偏导数：

式中，f_x和f_y分别是x方向和y方向上的地形偏导数，f(x,y)为地表高程数据。由上式可知，可用二阶差分格式计算三维起伏地形偏导数。

另外，还可以考虑到井炮震源有一定激发井深，采用接收点高程数据来形成地表高程测量数据，然后对地表高程测量数据进行二阶差分计算获得三维起伏地形偏导数。

步骤102)：获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时。

在本实施例中，炮点到成像点的中心射线坐标系、接收点到成像点的中心射线坐标系以及每个炮点和每个接收点分别到地下各成像点的射线方向和走时通过运动学射线追踪方法获得。

步骤103)：获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵。

在本实施例中，每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵通过动力学射线追踪方法获得。

步骤104)：根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值。

在本实施例中，三维起伏地表Beylkin行列式如下：

式中，h(x,ξ)表示三维起伏地表Beylkin行列式，x＝(x₁,x₂,x₃)表示地下成像点位置，ξ＝(ξ₁,ξ₂)表示炮点位置或接收点位置。其中，炮点位置可表示为x_s(ξ₁,ξ₂,ξ₃)＝(ξ_1s,ξ_2s,ξ_3s)；接收点位置可表示为x_r(ξ₁,ξ₂,ξ₃)＝(ξ_1r,ξ_2r,ξ_3r)。V(x)为成像点位置上速度。det|·|表示矩阵[·]的行列式。炮点位置x_s到成像点位置x的走时为T(x,x_s)，接收点位置x_r到成像点位置x的走时为T(x,x_r)。a^(s)和b^(s)分别表示炮点到成像点的射线在成像点位置的下倾角和方位角，a^(r)和b^(r)分别表示接收点到成像点的射线在成像点位置的下倾角和方位角。

对矩阵[·]中后两行的元素构成的矩阵进行分解，获得炮点走时场二阶偏导数矩阵及接收点走时场二阶偏导数矩阵。其中，炮点走时场二阶偏导数矩阵表示为：

接收点走时场二阶偏导数表示为：

其中，和分别表示炮点到成像点的射线在炮点位置的下倾角和方位角，和分别表示接收点到成像点的射线在接收点位置的下倾角和方位角。和分别是震源所在位置x方向和y方向的地形偏导数，和分别是接收点所在位置x方向和y方向的地形偏导数。炮点到成像点的中心射线坐标系和以及接收点到成像点的中心射线坐标系和是由运动学射线追踪确定的，几何扩散矩阵Q₂(x,x_s)和Q₂(x,x_r)分别是由炮点和接收点到成像点的动力学射线追踪确定的。

对于多次覆盖的地震资料，地震勘探工区内一个炮点激发的地震波将会由多个接收点来接收。换言之，一个接收点会依次接收来多个炮点的地震波。为了避免重复计算，炮点走时场二阶偏导数矩阵进一步地表示为：

接收点走时场二阶偏导数矩阵进一步地表示为：

因此，获得如下计算公式：

由此可知，通过三维地形偏导数每个炮点和每个接收点分别到地下各成像点的中心射线坐标系炮点到成像点的射线在炮点位置的下倾角和方位角接收点到成像点的射线在接收点位置的下倾角和方位角每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵Q₂(x,x_s)和Q₂(x,x_r)对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值D^(s)以及接收点走时场二阶偏导数矩阵的各元素值D^(r)。

步骤105)：根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

在实施例中，首先输入的背景速度模型为常速模型，速度为1000m/s。背景速度模型大小为50km×20km×6km。x方向起点坐标为0m、离散点个数为10001个、间距为5m。y方向起点坐标为-10000m、离散点个数为5个、间距为5000m。z方向起点坐标为0m、离散点个数为1200个、间距为5m。

为了对比水平地表和起伏地表对Beylkin行列式的影响，设置了三种地表：

(1)水平地表f_x＝0,f_y＝0；

(2)x方向起伏地表f_x＝1,f_y＝0，表示x方向45°斜坡、y方向水平；

(3)y方向起伏地表f_x＝0,f_y＝1，表示x方向水平、y方向45°斜坡。

本实施例中，为反映水平地表和起伏地表的Beylkin行列式的差异，这里仅计算一个点源(可以看成是炮点，也可以看成是接收点)的情况，该点位源于(15000m,0m,1000m)。成像区范围为：x方向从14000m到16000m,间距为50m；y方向从-100m到100m,间距为50m，z方向从0m到2000m,间距为50m。这里显示的剖面为y＝0m的纵剖面。这里显示了输出的走时、成像点上射线的下倾角和方位角以及Beylkin行列式中六个元素(d₂₁、d₂₂、d₂₃、d₃₁、d₃₂、d₃₃)。

如图2所示，为本申请实施例中常速背景介质中走时图。由上文可知，速度为1000m/s，呈等圆形，其中心位于震源位置。无论地表如何起伏，地下走时图保持不变，这意味地形起伏对走时没有影响。如图3所示，为本申请实施例中常速背景介质中射线的下倾角图。由图3可知，下倾角是射线与下垂线的夹角，图中下半部下倾角从0°到90°、上半部下倾角从90°到180°。如图4所示，为本申请实施例中常速背景介质中射线的方位角图。方位角是指射线在水平面上投影与x轴的夹角，图中右半部的方位角为0°、左半部的方位角为180°。图2、图3和图4都是由运动学射线追踪计算的，它们都不受地形起伏影响。

如图5所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₁图；如图6所示，为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₁图；如图7所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₂图；如图8所示，为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₂图；如图9所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₂₃图；如图10所示，为本申请实施例中起伏地表(f_x＝1,f_y＝0)的d₂₃图；如图11所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₁图；如图12所示，为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₁图；如图13所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₂图；如图14所示，为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₂图；如图15所示，为本申请实施例中水平地表(f_x＝0,f_y＝0)的d₃₃图；如图16所示，为本申请实施例中起伏地表(f_x＝0,f_y＝1)的d₃₃图。图5～图16为本申请实施例中常速背景介质中计算的Beylkin行列式中各元素对比图。其中，图5和图6分别为水平地表和起伏地表的d₂₁图，这里水平地形偏导数为f_x＝0,f_y＝0，而起伏地形偏导数取f_x＝1，f_y＝0，意味着起伏地形在x方向地形倾角为45°。对比图5和图6发现，d₂₁值以炮点为对称点，当地形倾斜时，这种对称图形发生了偏转，并且这种偏转角度随着地形倾角的增大而增大。由此看出，x方向地形起伏对Beylkin行列式中元素d₂₁的影响较大。图7和图8分别为水平地表和起伏地表的d₂₂图，这里地形偏导数同前。对比图7和图8发现，x方向地形起伏对d₂₂影响不大。同时也反映了d₂₂独特的花纹分布。图9和图10分别为水平地表和起伏地表的d₂₃图，这里地形偏导数同前。对比图9和图10发现，以炮点为对称点的d₂₃也随地形起伏而偏转，并且数值上随着地形偏导数的增加而明显增大。由此看出，x方向地形起伏对Beylkin行列式中元素d₂₃的影响很明显。为了考虑d₃₁、d₃₂和d₃₃受y方向地形起伏的影响，现在取f_x＝0,f_y＝1，即y方向地形倾角为45°。图11和图12分别为水平地表和起伏地表(f_x＝0,f_y＝1)的d₃₁图，对比发现y方向地形起伏对Beylkin行列式中元素d₃₁的影响也较大，同样导致了d₃₁图形的偏转。图13和图14分别为水平地表和起伏地表(f_x＝0,f_y＝1)的d₃₂图，显示了漂亮的花纹图案，然而两者没有明显差异。图15和图16分别为水平地表和起伏地表(f_x＝0,f_y＝1)的d₃₃图，两者有明显差异，表明y方向地形起伏对Beylkin行列式中元素d₃₃的影响很明显。

如图17所示，为本申请实施例提出的一种三维起伏地表地震资料偏移成像装置功能框图。包括：

三维地形偏导数单元201，用于获取三维地形偏导数；

中心射线坐标系、射线方向和走时单元202，用于获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；

振幅和几何扩散矩阵单元203，用于获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；

计算单元204，用于根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值；

成像单元205，用于根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

在本实施例中，所述三维地形偏导数单元201包括：

提取模块，用于从地震叠前数据道头字中提取地表高程数据；

第一二阶差分模块，用于对地表高程数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

在本实施例中，所述三维地形偏导数单元包括：

地表高程测量数据获取模块，用于根据接收点高程数据获取地表高程测量数据；

第二二阶差分模块，用于对地表高程测量数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

在本实施例中，所述中心射线坐标系、射线方向和走时单元通过运动学射线追踪方法获得每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时。

在本实施例中，所述振幅和几何扩散矩阵单元通过动力学射线追踪方法获得每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵。

在本实施例中，所述计算单元包括：

下倾角和方位角获取模块，用于根据每个炮点和每个接收点分别到地下各成像点的射线方向确定炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角；

三维起伏地表Beylkin行列式计算模块，用于根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系和几何扩散矩阵、炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角对三维起伏地表Beylkin行列式进行计算。

如图18所示，为本申请实施例提出的一种三维起伏地表地震资料偏移成像系统示意图。所述系统包括：存储器和处理器，所述存储器中存储计算机程序，所述计算机程序被所述处理器执行时，实现以下功能：

获取三维地形偏导数；

获取每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、射线方向和走时；

获取每个炮点和每个接收点分别到地下各成像点的振幅和几何扩散矩阵；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系、每个炮点和每个接收点分别到地下各成像点的射线方向、每个炮点和每个接收点分别到地下各成像点的几何扩散矩阵对三维起伏地表Beylkin行列式进行计算，确定炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值；

根据炮点走时场二阶偏导数矩阵的各元素值以及接收点走时场二阶偏导数矩阵的各元素值确定三维起伏地表Beylkin行列式的表达式，根据所述三维起伏地表Beylkin行列式的表达式、每个炮点和每个接收点分别到地下各成像点的走时和振幅对三维起伏地表地震资料偏移成像。

在本实施例中，获取三维地形偏导数，所述计算机程序被所述处理器执行时，实现以下功能：

从地震叠前数据道头字中提取地表高程数据；

对地表高程数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

在本实施例中，获取三维地形偏导数，所述计算机程序被所述处理器执行时，实现以下功能：

根据接收点高程数据获取地表高程测量数据；

对地表高程测量数据采用二阶差分格式计算获得三维起伏地表地形偏导数。

在本实施例中，对三维起伏地表Beylkin行列式进行计算，所述计算机程序被所述处理器执行时，实现以下功能：

根据每个炮点和每个接收点分别到地下各成像点的射线方向确定炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角；

根据三维地形偏导数、每个炮点和每个接收点分别到地下各成像点的中心射线坐标系和几何扩散矩阵、炮点到成像点的射线在成像点位置的下倾角和方位角、接收点到成像点的射线在成像点位置的下倾角和方位角、炮点到成像点的射线在炮点位置的下倾角和方位角、接收点到成像点的射线在接收点位置的下倾角和方位角对三维起伏地表Beylkin行列式进行计算。

在本实施方式中，所述存储器包括但不限于随机存取存储器(Random AccessMemory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard DiskDrive,HDD)或者存储卡(Memory Card)。

在本实施方式中，所述处理器可以按任何适当的方式实现。例如，所述处理器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application SpecificIntegrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式等等。

本说明书实施方式提供的三维起伏地表地震资料偏移成像系统，其存储器和处理器实现的具体功能，可以与本说明书中的前述实施方式相对照解释，并能够达到前述实施方式的技术效果，这里便不再赘述。

本技术方案提出一种新的地球物理油气勘探领域地震数据真振幅深度偏移成像或逆散射成像中Beylkin行列式的计算方案，利用该方案获取的计算结果实现山地、沙漠、黄土塬等三维起伏地表地面地震资料成像处理。

在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device,PLD)(例如现场可编程门阵列(Field Programmable GateArray，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(Altera Hardware DescriptionLanguage)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(RubyHardware Description Language)等，目前最普遍使用的是VHDL(Very-High-SpeedIntegrated Circuit Hardware Description Language)与Verilog2。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。

本领域技术人员也知道，除了以纯计算机可读程序代码方式实现客户端、服务器以外，完全可以通过将方法步骤进行逻辑编程来使得客户端、服务器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种客户端、服务器可以被认为是一种硬件部件，而对其内部包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本申请各个实施方式或者实施方式的某些部分所述的方法。

本说明书中的各个实施方式均采用递进的方式描述，各个实施方式之间相同相似的部分互相参见即可，每个实施方式重点说明的都是与其他实施方式的不同之处。尤其，针对客户端的实施方式来说，均可以参照前述方法的实施方式的介绍对照解释。

本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

虽然通过实施方式描绘了本申请，本领域普通技术人员知道，本申请有许多变形和变化而不脱离本申请的精神，希望所附的权利要求包括这些变形和变化而不脱离本申请的精神。