CN104965222B - 三维纵波阻抗全波形反演方法及装置 - Google Patents

Abstract

本发明公开一种三维纵波阻抗全波形反演方法及装置，该方法包括：确定时间域三维一阶速度‑应力变密度声波方程相应的伴随方程、纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；根据时间域三维一阶速度‑应力变密度声波方程及相应的伴随方程，采用高阶交错网格有限差分法，确定正向传播波场和残差逆时反传波场；根据正向传播波场、残差逆时反传波场和纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波速度梯度、密度梯度和纵波阻抗梯度；根据纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。本发明可提高反演过程的可行性，并达到良好的数值频散抑制能力。

Description

三维纵波阻抗全波形反演方法及装置

技术领域

本发明涉及石油勘探地震弹性参数反演技术领域，尤其涉及三维纵波阻抗全波形反演方法及装置。

背景技术

传统的弹性参数反演方法，如基于Zeoppritz方程的弹性阻抗反演方法和AVO(Amplitude Versus Offset，振幅随偏移距的变化)类反演方法，主要存在以下问题：一是受到平面波入射假设的限制，基于褶积模型的叠前地震数据反射振幅的假设不适应大角度入射的情况；二是井稀疏的情况下，很难得到可靠的低频模型。这些问题使得弹性阻抗反演方法和AVO类反演方法不适于对复杂岩性和油气藏的精细描述。上世纪80年代，Tarantola和Mora等提出了基于波动方程的纵横波速度和密度等的反演方法，即全波形反演方法，该方法自动考虑了地震波的反射、透射及波形转换等问题，充分利用地震波的运动学和动力学特征，能够反演出比较可靠的速度以及反映储层物性的弹性参数和密度等信息。

目前，尽管大部分全波形反演实例都是为深度域成像提供高精度的速度模型，也有一些用于储层定量表征的例子。密度是岩石物性参数的重要组成部分，如果能通过反演同时得到比较可靠的纵横波速度、密度资料，这对储层评价、岩性解释和油藏描述等都具有重要作用。弹性波多参数全波形反演研究表明，各参数之间的相互耦合，会增加多参数全波形反演的非线性程度。为降低反演问题的复杂性以及解的稳定性和可靠性，可以减少反演参数，利用声波方程进行双参数全波形反演，得到高精度纵波速度和密度，接着就可以得到纵波阻抗。

Jeong等(2012)提出了一种频率域密度反演策略，首先将密度固定为任意常数，通过单参数反演得到弹性模量，进而将换算得到的纵横波速度作为初始模型，再同时反演速度和密度；Prieux等(2013)先利用大偏移距数据单参数反演速度参数，再利用全偏移距数据同时反演速度、密度参数，有效提高了密度反演的稳定性，但密度反演结果仍然不够理想。他们同时指出，利用速度-密度参数化方式反演得到的速度、密度和阻抗结果要好于利用速度-阻抗参数化方式得到的反演结果。杨积忠等(2014)从变密度声波方程出发，首先研究了密度在速度反演中的重要作用，然后分析了速度对密度反演的影响程度，进而提出了一种有利于速度、密度分步联合反演的策略：第一步，利用给定的初始模型对速度、密度进行同时反演，得到比较可靠的速度反演结果；第二步，利用第一步反演得到的速度和给定的初始密度作为初始模型，继续进行双参数同时反演，这样可以同时得到比较可靠的速度、密度反演结果。

可见由于纵波速度和密度是相互耦合的，同时进行纵波速度和密度反演是比较可行的反演策略，但是，Jeong等、Prieux等和杨积忠等采用的是二阶标量频率域变密度声波方程，而频率域正演内存需求巨大，特别是三维情况下LU矩阵(线性代数中可以分解为一个下三角矩阵和一个上三角矩阵乘积的矩阵，其中L和U分别是下三角和上三角矩阵)的建立对于内存需求极高，甚至难以实现；同时，二阶标量声波方程通常利用高阶中心网格有限差分法求解，要降低数值频散的影响就需要增加有限差分阶数，从而会增加计算量。

发明内容

本发明实施例提供一种三维纵波阻抗全波形反演方法，用以降低对内存的需求，提高反演过程的可行性，并达到良好的数值频散抑制能力，该方法包括：

确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

确定所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

根据所述时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

根据所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

根据所述正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；

根据所述纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

根据所述纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。

本发明实施例还提供一种三维纵波阻抗全波形反演装置，用以降低对内存的需求，提高反演过程的可行性，并达到良好的数值频散抑制能力，该装置包括：

伴随方程确定模块，用于确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

计算公式确定模块，用于确定所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

正向波场确定模块，用于根据所述时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

残差波场确定模块，用于根据所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

梯度计算模块，用于根据所述正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；根据所述纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

参数模型更新模块，用于根据所述纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。

本发明实施例中，基于时间域三维一阶速度-应力变密度声波方程，在时间域进行全波形反演，对于内存需求较小，可行性高；并且，采用时间域三维一阶速度-应力变密度声波方程，可以利用高阶交错网格有限差分法求解，在阶数相同时，能够比中心网格有限差分法有更好的数值频散抑制能力。

附图说明

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

图1为本发明实施例中三维纵波阻抗全波形反演方法的示意图；

图2为本发明实施例中三维纵波阻抗全波形反演装置的示意图；

图3为本发明实施例中真实纵波阻抗切片的示意图；

图4为本发明实施例中初始纵波阻抗切片的示意图；

图5为本发明实施例中反演得到的纵波阻抗切片的示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

为了降低对内存的需求，提高反演过程的可行性，并达到良好的数值频散抑制能力，本发明实施例提供一种三维纵波阻抗全波形反演方法。图1为本发明实施例中三维纵波阻抗全波形反演方法的示意图，如图1所示，该方法可以包括：

步骤101、确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

步骤102、确定时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

步骤103、根据时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

步骤104、根据时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

步骤105、根据正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；根据纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

步骤106、根据纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。

由图1所示流程可以得知，本发明实施例与现有技术采用二阶标量频率域变密度声波方程不同的是，基于时间域三维一阶速度-应力变密度声波方程，在时间域进行全波形反演，对于内存需求较小，可行性高；并且，采用时间域三维一阶速度-应力变密度声波方程，可以利用高阶交错网格有限差分法求解，在阶数相同时，能够比中心网格有限差分法有更好的数值频散抑制能力。

具体实施时，本发明实施例充分利用地震波的运动学和动力学特征，克服传统的基于Zeoppritz方程的弹性阻抗反演方法和AVO类反演方法的缺点，得到比较可靠的反映储层物性的纵波阻抗信息，基于时间域三维一阶速度-应力变密度声波方程进行纵波阻抗的全波形反演，通过采用局部最优化算法，在反演过程中不断更新纵波阻抗模型，进而得到较高精度的纵波阻抗反演结果，为岩性解释和储层含油气性检测提供数据支持。

实施例中，基于时间域三维一阶速度-应力变密度声波方程，采用基于摄动理论的伴随状态法，可以将全波形反演问题转化为包含伴随波场变量的最优化问题。具体的，可以根据包含伴随波场变量的最优化问题，确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程。例如，可以通过对观测波场变量求偏导数，并令这些偏导数为零，结合边界条件及终止条件，得到时间域三维一阶速度-应力变密度声波方程相应的伴随方程。还可以根据包含伴随波场变量的最优化问题，确定时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式。例如，通过对模型参数求偏导数，结合边界条件及终止条件，得到纵波速度梯度、密度梯度及纵波阻抗梯度计算公式。

具体实施时，在确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程、纵波速度梯度、密度梯度及纵波阻抗梯度计算公式之后，根据时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场。实施例中，可以利用真实纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解时间域三维一阶速度-应力变密度声波方程，得到观测炮记录；利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解时间域三维一阶速度-应力变密度声波方程，得到正向传播波场各时刻的快照及合成炮记录。

确定正向传播波场后，再确定残差逆时反传波场。实施例中，根据时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场，可以包括：利用观测炮记录和合成炮记录，得到残差记录；利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法逆时外推时间域三维一阶速度-应力变密度声波方程相应的伴随方程，得到残差逆时反传波场各时刻的快照。

具体实施时，根据正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度，可以包括：利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度。例如，先利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，得到纵波速度梯度；再利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，得到密度梯度；后续可以利用纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，得到纵波阻抗梯度。

实施例中，根据纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。其中，在每次迭代中，可以采用共轭梯度法更新纵波速度和密度模型，共轭方向采用弗莱彻和李维斯提出的FR(Fletcher-Reeves)方法来计算；在每次迭代中，可以采用共轭梯度法更新纵波阻抗模型，共轭方向采用FR方法来计算。这样，基于时间域三维一阶速度-应力变密度声波方程进行纵波阻抗的全波形反演，通过采用局部最优化算法，在反演过程中不断更新纵波阻抗模型，进而得到较高精度的纵波阻抗反演结果，为岩性解释和储层含油气性检测提供高精度纵波阻抗数据体。

实施例中，上述时间域三维一阶速度-应力变密度声波方程可以如下：

其中，p为流体平均压力，v_x为x方向波场速度分量，v_y为y方向波场速度分量，v_z为z方向波场速度分量，t为时间，x、y、z为波场方向，ρ为介质密度，f_p为震源项，K为弹性模量v_p为纵波速度；

上述时间域三维一阶速度-应力变密度声波方程相应的伴随方程可以如下：

其中，λ_x、λ_y、λ_z分别为x、y、z方向上的伴随速度分量，λ_p为伴随压力波场，p^obs为观测压力波场，x为任意点坐标，x_r为接收点坐标；

上述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度计算公式可以如下：

其中，G_ρ为纵波速度梯度，T为波场传播时间；

上述时间域三维一阶速度-应力变密度声波方程相应的密度梯度计算公式可以如下：

其中，为密度梯度；

上述时间域三维一阶速度-应力变密度声波方程相应的纵波阻抗梯度计算公式可以如下：

其中，为纵波阻抗梯度。

基于同一发明构思，本发明实施例中还提供了一种三维纵波阻抗全波形反演装置，如下面的实施例所述。由于该装置解决问题的原理与三维纵波阻抗全波形反演方法相似，因此该装置的实施可以参见三维纵波阻抗全波形反演方法的实施，重复之处不再赘述。

图2为本发明实施例中三维纵波阻抗全波形反演装置的示意图。如图2所示，本发明实施例中三维纵波阻抗全波形反演装置可以包括：

伴随方程确定模块201，用于确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

计算公式确定模块202，用于确定时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

正向波场确定模块203，用于根据时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

残差波场确定模块204，用于根据时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

梯度计算模块205，用于根据正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；根据纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

参数模型更新模块206，用于根据纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果。

具体实施时，伴随方程确定模块201具体可以用于：

通过对观测波场变量求偏导数，并令这些偏导数为零，结合边界条件及终止条件，得到时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

计算公式确定模块202具体可以用于：

通过对模型参数求偏导数，结合边界条件及终止条件，得到纵波速度梯度、密度梯度及纵波阻抗梯度计算公式。

具体实施时，正向波场确定模块203具体可以用于：

利用真实纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解时间域三维一阶速度-应力变密度声波方程，得到观测炮记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解时间域三维一阶速度-应力变密度声波方程，得到正向传播波场各时刻的快照及合成炮记录；

残差波场确定模块204具体可以用于：

利用观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法逆时外推时间域三维一阶速度-应力变密度声波方程相应的伴随方程，得到残差逆时反传波场各时刻的快照。

具体实施时，梯度计算模块205具体可以用于：

利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度。

具体实施时，参数模型更新模块206具体可以用于：

在每次迭代中，采用共轭梯度法更新纵波速度和密度模型，共轭方向采用FR方法来计算；

在每次迭代中，采用共轭梯度法更新纵波阻抗模型，共轭方向采用FR方法来计算。

具体实施时，上述时间域三维一阶速度-应力变密度声波方程可以如下：

其中，p为流体平均压力，v_x为x方向波场速度分量，v_y为y方向波场速度分量，v_z为z方向波场速度分量，t为时间，x、y、z为波场方向，ρ为介质密度，f_p为震源项，K为弹性模量v_p为纵波速度；

上述时间域三维一阶速度-应力变密度声波方程相应的伴随方程可以如下：

其中，λ_x、λ_y、λ_z分别为x、y、z方向上的伴随速度分量，λ_p为伴随压力波场，p^obs为观测压力波场，x为任意点坐标，x_r为接收点坐标；

上述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度计算公式可以如下：

其中，G_ρ为纵波速度梯度，T为波场传播时间；

上述时间域三维一阶速度-应力变密度声波方程相应的密度梯度计算公式可以如下：

其中，为密度梯度；

上述时间域三维一阶速度-应力变密度声波方程相应的纵波阻抗梯度计算公式可以如下：

其中，为纵波阻抗梯度。

由上述实施例可知，本发明实施例中基于时间域一阶速度-应力变密度声波方程的三维纵波阻抗全波形反演方法及装置，由于充分利用地震波的运动学和动力学特征，能够反演出比较可靠的反映储层物性的纵波阻抗信息；由于采用高阶交错网格有限差分法，比基于二阶方程的声波反演方法在采用相同差分阶数时具有更好地数值频散抑制能力；由于在时间域实现，因此波场正演和残差逆时反传比较直接快速。

图3为本发明实施例中真实纵波阻抗切片的示意图；图4为本发明实施例中初始纵波阻抗切片的示意图；图5为本发明实施例中反演得到的纵波阻抗切片的示意图。如图3、图4、图5所示，与真实纵波阻抗切片和初始纵波阻抗切片比较可见，本发明实施例的反演方法具有较强的细节刻画能力，能够得到较为准确的纵波阻抗反演结果。

综上所述，本发明实施例中基于时间域一阶速度-应力变密度声波方程的三维纵波阻抗全波形反演方法及装置，由于采用高阶交错网格有限差分法，比基于二阶方程的声波反演方法在采用相同差分阶数时具有更好地数值频散抑制能力；由于在时间域实现，因此波场正演和残差逆时反传比较直接快速；由于充分利用地震波的运动学和动力学特征，能够反演出比较可靠的反映储层物性的纵波阻抗信息。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种三维纵波阻抗全波形反演方法，其特征在于，包括：

确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

确定所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

根据所述时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

根据所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

根据所述正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；

根据所述纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

根据所述纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果；

根据所述时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场，包括：

利用真实纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解所述时间域三维一阶速度-应力变密度声波方程，得到观测炮记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解所述时间域三维一阶速度-应力变密度声波方程，得到正向传播波场各时刻的快照及合成炮记录；

根据所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场，包括：

利用所述观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法逆时外推所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，得到残差逆时反传波场各时刻的快照。

2.如权利要求1所述的方法，其特征在于，确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程，包括：

通过对观测波场变量求偏导数，并令这些偏导数为零，结合边界条件及终止条件，得到时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

确定所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，包括：

通过对模型参数求偏导数，结合边界条件及终止条件，得到纵波速度梯度、密度梯度及纵波阻抗梯度计算公式。

3.如权利要求1所述的方法，其特征在于，根据所述正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度，包括：

利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度。

4.如权利要求1所述的方法，其特征在于，根据所述纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果，包括：

在每次迭代中，采用共轭梯度法更新纵波速度和密度模型，共轭方向采用弗莱彻-李维斯FR方法来计算；

在每次迭代中，采用共轭梯度法更新纵波阻抗模型，共轭方向采用FR方法来计算。

5.如权利要求1至4任一项所述的方法，其特征在于，所述时间域三维一阶速度-应力变密度声波方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}+\frac{\∂p}{\∂x}=0\\ \mathrm{\ρ}\frac{\∂v_y}{\∂t}+\frac{\∂p}{\∂y}=0\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}+\frac{\∂p}{\∂z}=0\\ \frac{\∂p}{\∂t}+K(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})-f_p=0\end{array}\right.;$

其中，p为流体平均压力，v_x为x方向波场速度分量，v_y为y方向波场速度分量，v_z为z方向波场速度分量，t为时间，x、y、z为波场方向，ρ为介质密度，f_p为震源项，K为弹性模量v_p为纵波速度；

所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_x}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂x}=0\\ \mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_y}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂y}=0\\ \mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_z}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂z}=0\\ \frac{\∂{\mathrm{\λ}}_p}{\∂t}+(\frac{\∂{\mathrm{\λ}}_x}{\∂x}+\frac{\∂{\mathrm{\λ}}_y}{\∂y}+\frac{\∂{\mathrm{\λ}}_z}{\∂z})=(p-p^{obs})\mathrm{\δ}(x-x_r)\end{array}\right.;$

其中，λ_x、λ_y、λ_z分别为x、y、z方向上的伴随速度分量，λ_p为伴随压力波场，p^obs为观测压力波场，x为任意点坐标，x_r为接收点坐标；

所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度计算公式如下：

$G_{\mathrm{\ρ}}={\∫}_0^T({\mathrm{\λ}}_x\frac{\∂v_x}{\∂t}+{\mathrm{\λ}}_y\frac{\∂v_y}{\∂t}+{\mathrm{\λ}}_z\frac{\∂v_z}{\∂t})dt+v_p^2{\∫}_0^T{\mathrm{\λ}}_p(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})dt;$

其中，G_ρ为纵波速度梯度，T为波场传播时间；

所述时间域三维一阶速度-应力变密度声波方程相应的密度梯度计算公式如下：

$G_{v_p}=2{\mathrm{\ρv}}_p{\∫}_0^T{\mathrm{\λ}}_p(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})dt;$

其中，为密度梯度；

所述时间域三维一阶速度-应力变密度声波方程相应的纵波阻抗梯度计算公式如下：

$G_{I_p}=G_{v_p}\·G_{\mathrm{\ρ}};$

其中，为纵波阻抗梯度。

6.一种三维纵波阻抗全波形反演装置，其特征在于，包括：

伴随方程确定模块，用于确定时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

计算公式确定模块，用于确定所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度、密度梯度及纵波阻抗梯度计算公式；

正向波场确定模块，用于根据所述时间域三维一阶速度-应力变密度声波方程，采用高阶交错网格有限差分法，确定正向传播波场；

残差波场确定模块，用于根据所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，采用高阶交错网格有限差分法，确定残差逆时反传波场；

梯度计算模块，用于根据所述正向传播波场、残差逆时反传波场及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度；根据所述纵波速度梯度、密度梯度及纵波阻抗梯度计算公式，确定纵波阻抗梯度；

参数模型更新模块，用于根据所述纵波速度梯度、密度梯度和纵波阻抗梯度，采用局部最优化算法对参数模型进行更新，输出满足收敛条件的三维纵波阻抗全波形反演结果；

所述正向波场确定模块具体用于：

利用真实纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解所述时间域三维一阶速度-应力变密度声波方程，得到观测炮记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法数值求解所述时间域三维一阶速度-应力变密度声波方程，得到正向传播波场各时刻的快照及合成炮记录；

所述残差波场确定模块具体用于：

利用所述观测炮记录和合成炮记录，得到残差记录；

利用初始纵波速度模型和密度模型，采用高阶交错网格有限差分法逆时外推所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程，得到残差逆时反传波场各时刻的快照。

7.如权利要求6所述的装置，其特征在于，所述伴随方程确定模块具体用于：

通过对观测波场变量求偏导数，并令这些偏导数为零，结合边界条件及终止条件，得到时间域三维一阶速度-应力变密度声波方程相应的伴随方程；

所述计算公式确定模块具体用于：

通过对模型参数求偏导数，结合边界条件及终止条件，得到纵波速度梯度、密度梯度及纵波阻抗梯度计算公式。

8.如权利要求6所述的装置，其特征在于，所述梯度计算模块具体用于：

利用正向传播波场各时刻的快照及残差逆时反传波场各时刻的快照，及纵波速度梯度和密度梯度计算公式，确定纵波速度梯度和密度梯度。

9.如权利要求6所述的装置，其特征在于，所述参数模型更新模块具体用于：

在每次迭代中，采用共轭梯度法更新纵波速度和密度模型，共轭方向采用FR方法来计算；

在每次迭代中，采用共轭梯度法更新纵波阻抗模型，共轭方向采用FR方法来计算。

10.如权利要求6至9任一项所述的装置，其特征在于，所述时间域三维一阶速度-应力变密度声波方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂v_x}{\∂t}+\frac{\∂p}{\∂x}=0\\ \mathrm{\ρ}\frac{\∂v_y}{\∂t}+\frac{\∂p}{\∂y}=0\\ \mathrm{\ρ}\frac{\∂v_z}{\∂t}+\frac{\∂p}{\∂z}=0\\ \frac{\∂p}{\∂t}+K(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})-f_p=0\end{array}\right.;$

其中，p为流体平均压力，v_x为x方向波场速度分量，v_y为y方向波场速度分量，v_z为z方向波场速度分量，t为时间，x、y、z为波场方向，ρ为介质密度，f_p为震源项，K为弹性模量v_p为纵波速度；

所述时间域三维一阶速度-应力变密度声波方程相应的伴随方程如下：

$\left\{\begin{array}{c}\mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_x}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂x}=0\\ \mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_y}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂y}=0\\ \mathrm{\ρ}\frac{\∂{\mathrm{\λ}}_z}{\∂t}+K\frac{\∂{\mathrm{\λ}}_p}{\∂z}=0\\ \frac{\∂{\mathrm{\λ}}_p}{\∂t}+(\frac{\∂{\mathrm{\λ}}_x}{\∂x}+\frac{\∂{\mathrm{\λ}}_y}{\∂y}+\frac{\∂{\mathrm{\λ}}_z}{\∂z})=(p-p^{obs})\mathrm{\δ}(x-x_r)\end{array}\right.;$

其中，λ_x、λ_y、λ_z分别为x、y、z方向上的伴随速度分量，λ_p为伴随压力波场，p^obs为观测压力波场，x为任意点坐标，x_r为接收点坐标；

所述时间域三维一阶速度-应力变密度声波方程相应的纵波速度梯度计算公式如下：

$G_{\mathrm{\ρ}}={\∫}_0^T({\mathrm{\λ}}_x\frac{\∂v_x}{\∂t}+{\mathrm{\λ}}_y\frac{\∂v_y}{\∂t}+{\mathrm{\λ}}_z\frac{\∂v_z}{\∂t})dt+v_p^2{\∫}_0^T{\mathrm{\λ}}_p(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})dt;$

其中，G_ρ为纵波速度梯度，T为波场传播时间；

所述时间域三维一阶速度-应力变密度声波方程相应的密度梯度计算公式如下：

$G_{v_p}=2{\mathrm{\ρv}}_p{\∫}_0^T{\mathrm{\λ}}_p(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y}+\frac{\∂v_z}{\∂z})dt;$

其中，为密度梯度；

所述时间域三维一阶速度-应力变密度声波方程相应的纵波阻抗梯度计算公式如下：

$G_{I_p}=G_{v_p}\·G_{\mathrm{\ρ}};$

其中，为纵波阻抗梯度。