一种地震反演方法和装置与流程

本发明涉及地质勘探技术领域，特别涉及一种地震反演方法和装置。

背景技术：

地震反演可以指的是利用地表观测到的地震资料，以已知地质规律和钻井、测井资料为约束，推测出地层岩性构造的过程。从地震反演所用的地震资料来分，地震反演可分为：叠前反演和叠后反演；从反演所利用地震的信息来分：地震反演可分为：地震波旅行时反演和地震波振幅反演；从反演的地质结果来分，地震反演可分为：构造反演、波阻抗反演、储层参数反演、地质统计反演等。

目前，一般可以采用以下方法来实现地震储层中具有垂向对称轴的横向各向同性(verticaltransverseisotropy，简称为vti)介质参数的反演：

1)基于各向同性假设的反演方法：

然而，该方法不能反演得到具有各向异性特征的页岩储层的vti参数，且由于在vti各向异性条件下，各项同性反射系数计算公式计算得到的反射系数存在较大的误差，使得基于各向同性假设的反演方法反演得到的常规弹性参数精度不高。

2)全波形反演方法：

虽然该方法可以利用全波场信息来预测储层的vti参数。然而，该方法计算量巨大，在反演尺度和计算效率上不能满足实际油藏储层的要求。

3)基于vti介质的振幅随角度变化的(amplitudeversusangle，简称为ava)反演方法：

该方法和基于各向同性假设的反演方法以及全波形反演方法相比，地震数据解释的可信度得到了提高，主要可以通过结合vti介质的精确反射系数公式以及叠前道集上振幅随偏移距变化等信息，来确定地震储层中的vti参数(弹性参数和各向异性参数)。然而，进行反演的过程中，由于精确反射系数近似公式非常复杂，很难直接用于反演中。同时，由于vti介质的精确反射系数公式计算精度较低，在强各向异性、强阻抗差以及大入射角等情况下，反演得到的vti参数存在较大的误差，不能满足实际油藏储层的精细表征要求。

技术实现要素：

本发明提供了一种地震反演方法和装置，以达到提高vti参数的反演精度，指导油田高效开发的目的。

本发明实施例提供了一种地震反演方法，可以包括：基于待测储层中纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面纵波的垂直慢度关系式；根据待测储层中横波速度在垂直方向的取值、所述纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面横波的垂直慢度关系式；根据所述待测储层平面纵波的垂直慢度关系式、所述待测储层平面横波的垂直慢度关系式，构建新的具有垂向对称轴的横向各向同性介质的反射系数计算方式；利用所述反射系数计算方式，对所述待测储层进行振幅随角度变化的反演。

在一个实施例中，确定所述待测储层平面纵波的垂直慢度，可以包括：计算所述待测储层中所述纵波速度在垂直方向的取值，并计算所述取值的倒数的平方，作为第一储层参数；计算所述水平慢度的平方，并将所述水平慢度的平方与所述纵波速度在近垂直入射时随相位角的变化情况相乘，作为第二储层参数；计算所述纵波速度在横向和垂向的差异值、所述纵波速度在近垂直入射时随相位角的变化情况之差，作为速度差异值，计算所述纵波速度在垂直方向的取值的平方，计算所述水平慢度的四次方，计算所述速度差异值、所述取值的平方、所述水平慢度的四次方之间的差值，并将所得到的差值作为第三储层参数；计算所述第三储层参数的方根，并将所述方根作为所述待测储层平面纵波的垂直慢度。

在一个实施例中，按照以下公式确定所述待测储层平面纵波的垂直慢度：

上式中，qα表示所述待测储层平面纵波的垂直慢度，vp0表示所述纵波速度在垂直方向的取值，δ表示所述纵波速度在近垂直入射时随相位角的变化情况，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值。

在一个实施例中，确定所述待测储层平面横波的垂直慢度，可以包括：计算所述待测储层中横波速度在垂直方向的取值的倒数，并计算所述倒数的平方，作为第四储层参数；计算所述纵波速度在垂直方向的取值的平方与所述水平慢度的四次方之间的乘积，得到第五储层参数；计算所述纵波速度在垂直方向的取值与所述横波速度在垂直方向的取值的比值的平方，并计算所述比值的平方与所述水平慢度的平方的乘积作为水平值，将所述比值的平方与所述水平值相乘，得到第六储层参数；计算所述第五储层参数和所述第六储层参数的差值，作为第七储层参数；计算所述速度差异值和所述第七储层参数的乘积，并将所得到的乘积结果和所述第四储层参数相加，将相加后的结果作为第八储层参数；计算所述第八储层参数、所述水平慢度的平方之间的差值，将所得到的差值作为所述待测储层平面横波的垂直慢度。

在一个实施例中，按照以下公式确定所述待测储层平面横波的垂直慢度：

上式中，qβ表示所述待测储层平面横波的垂直慢度，vs0表示所述横波速度在垂直方向的取值，vp0表示所述待测储层中纵波速度在垂直方向的取值，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值，δ表示所述纵波波速度在近垂直入射时随相位角的变化情况。

本发明实施例还提供了一种地震反演装置，可以包括：纵波确定模块，用于基于待测储层中纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面纵波的垂直慢度关系式；横波确定模块，用于根据待测储层中横波速度在垂直方向的取值、所述纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面横波的垂直慢度关系式；反射系数计算模块，用于根据所述待测储层平面纵波的垂直慢度关系式、所述待测储层平面中横波的垂直慢度关系式，构建新的具有垂向对称轴的横向各向同性介质的反射系数计算方式；反演模块，用于利用所述反射系数计算方式，对所述待测储层进行振幅随角度变化的反演。

在一个实施例中，所述纵波确定模块可以包括：第一储层参数计算单元，用于计算所述待测储层中所述纵波速度在垂直方向的取值，并计算所述取值的倒数的平方，作为第一储层参数；第二储层参数计算单元，用于计算所述水平慢度的平方，并将所述水平慢度的平方与所述纵波速度在近垂直入射时随相位角的变化情况相乘，作为第二储层参数；第三储层参数计算单元，用于计算所述纵波速度在横向和垂向的差异值、所述纵波速度在近垂直入射时随相位角的变化情况之差，作为速度差异值，计算所述纵波速度在垂直方向的取值的平方，计算所述水平慢度的四次方，计算所述速度差异值、所述取值的平方、所述水平慢度的四次方之间的差值，并将所得到的差值作为第三储层参数；纵波的垂直慢度计算单元，用于计算所述第三储层参数的方根，并将所述方根作为所述待测储层平面纵波的垂直慢度。

在一个实施例中，所述纵波确定模块具体可以用于按照以下公式确定所述待测储层平面纵波的垂直慢度：

上式中，qα表示所述待测储层平面纵波的垂直慢度，vp0表示所述纵波速度在垂直方向的取值，δ表示所述纵波速度在近垂直入射时随相位角的变化情况，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值。

在一个实施例中，所述横波确定模块可以包括：第四储层参数计算单元，用于计算所述待测储层中横波速度在垂直方向的取值的倒数，并计算所述倒数的平方，作为第四储层参数；第五储层参数计算单元，用于计算所述纵波速度在垂直方向的取值的平方与所述水平慢度的四次方之间的乘积，得到第五储层参数；第六储层参数计算单元，用于计算所述纵波速度在垂直方向的取值与所述横波速度在垂直方向的取值的比值的平方，并计算所述比值的平方与所述水平慢度的平方的乘积作为水平值，将所述比值的平方与所述水平值相乘，得到第六储层参数；第七储层参数计算单元，用于计算所述第五储层参数和所述第六储层参数的差值，作为第七储层参数；第八储层参数计算单元，用于计算所述速度差异值和所述第七储层参数的乘积，并将所得到的乘积结果和所述第四储层参数相加，将相加后的结果作为第八储层参数；横波的垂直慢度计算单元，用于计算所述第八储层参数、所述水平慢度的平方之间的差值，将所得到的差值作为所述待测储层平面横波的垂直慢度。

在一个实施例中，所述横波确定模块具体可以用于按照以下公式确定所述待测储层平面横波的垂直慢度：

上式中，qβ表示所述待测储层平面横波的垂直慢度，vs0表示所述横波速度在垂直方向的取值，vp0表示所述待测储层中纵波速度在垂直方向的取值，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值，δ表示所述纵波波速度在近垂直入射时随相位角的变化情况。

在本发明实施例中，提供了一种新的所述待测储层平面纵波的垂直慢度公式、所述待测储层平面横波的垂直慢度公式，并利用上述两个公式所构建的vti介质的反射系数计算方式对所述待测储层进行ava反演的方法。上述两个公式不仅具有足够高的计算精度，还具有相对简单的形式和相对较弱的非线性，利用本申请所提出的所述待测储层平面纵波的垂直慢度公式、所述待测储层平面横波的垂直慢度公式对所述vti介质的反射系数计算方式进行简化后，得到的简化后的vti介质的反射系数计算方式进行ava反演时，很好的克服了现有的反演技术中由于假设条件太多而导致的反演精度较低等问题。进一步的，简化后的vti介质的反射系数计算方式在强各向异性、强阻抗差以及大入射角等情况下均具有较高的计算精度。利用本申请所提出的反演结果进行油田勘探时，勘探精度得到了有效的提高。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，

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

图1是本申请提供的一种地震反演方法流程图；

图2是本申请提供的一种地震反演方法反演得到的vti参数，其中，图a表示纵波速度vp，图b表示横波速度vs，图c表示密度，图d表示横向和垂向p波速度的差异结果ε，图e表示p波速度在近垂直入射时随相位角的变化情况δ(e)；

图3是本申请提供的一种地震反演装置的结构框图。

具体实施方式

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

考虑到现有技术中基于vti介质进行ava反演时，由于现有的精确的vti介质反射系数公式非常复杂所导致的油藏勘探精度较低的问题，发明人提出了一种新的vti介质反射系数计算方式的简化方法，即，根据所述待测储层平面p波的垂直慢度公式、所述待测储层平面sv波的垂直慢度公式对vti介质反射系数公式进行简化，再利用简化后的vti介质反射系数计算方式对待测储层进行ava反演。基于此，提出了一种地震反演方法，如图1所示，可以包括以下步骤：

s101：基于待测储层中纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面纵波的垂直慢度关系式。

具体的，所述待测储层可以是页岩储层，也可以是砂岩储层、碳酸盐岩储层等其他各种类型的待测储层，本申请对此不作限定。

其中，在各向异性介质中，不同方向的纵波传播速度是不一样的，所述纵波速度可以是和入射角度有关的。不同入射角度情况下的纵波速度等于垂直入射也就是入射角为0度时的纵波速度乘以一个和角度有关的表达式。

根据波的质点振动方向与传播方向的不同，可以分为横波(s波)以及纵波(p波)。质点振动方向与传播方向相同的波称为p波；质点振动方向与传播方向垂直的波称为s波。质点振动发生在与波的传播面相垂直的面内的波称为横波(sv波)，质点振动发生在与波的传播面相平行的面内的波称为sh波。进一步地，pp波可以是所述p波入射之后的反射波，p波速度可以是纵波速度。

在本申请的一个实施例中，确定所述待测储层平面纵波的垂直慢度，可以包括：计算所述待测储层中所述纵波速度在垂直方向的取值，并计算所述取值的倒数的平方，作为第一储层参数；计算所述水平慢度的平方，并将所述水平慢度的平方与所述纵波速度在近垂直入射时随相位角的变化情况相乘，作为第二储层参数；计算所述纵波速度在横向和垂向的差异值、所述纵波速度在近垂直入射时随相位角的变化情况之差，作为速度差异值，计算所述纵波速度在垂直方向的取值的平方，计算所述水平慢度的四次方，计算所述速度差异值、所述取值的平方、所述水平慢度的四次方之间的差值，并将所得到的差值作为第三储层参数；计算所述第三储层参数的方根，并将所述方根作为所述待测储层平面纵波的垂直慢度。

具体的，可以按照以下公式确定所述待测储层平面纵波的垂直慢度：

上式中，qα表示所述待测储层平面纵波的垂直慢度，vp0表示所述纵波速度在垂直方向的取值，δ表示所述纵波速度在近垂直入射时随相位角的变化情况，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值。

s102：根据待测储层中横波速度在垂直方向的取值、所述纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面横波的垂直慢度关系式。

在本申请的一个实施例中，确定所述待测储层平面横波的垂直慢度，可以包括：计算所述待测储层中横波速度在垂直方向的取值的倒数，并计算所述倒数的平方，作为第四储层参数；计算所述纵波速度在垂直方向的取值的平方与所述水平慢度的四次方之间的乘积，得到第五储层参数；计算所述纵波速度在垂直方向的取值与所述横波速度在垂直方向的取值的比值的平方，并计算所述比值的平方与所述水平慢度的平方的乘积作为水平值，将所述比值的平方与所述水平值相乘，得到第六储层参数；计算所述第五储层参数和所述第六储层参数的差值，作为第七储层参数；计算所述速度差异值和所述第七储层参数的乘积，并将所得到的乘积结果和所述第四储层参数相加，将相加后的结果作为第八储层参数；计算所述第八储层参数、所述水平慢度的平方之间的差值，将所得到的差值作为所述待测储层平面横波的垂直慢度。

具体的，按照以下公式确定所述待测储层平面横波的垂直慢度：

上式中，qβ表示所述待测储层平面横波的垂直慢度，vs0表示所述横波速度在垂直方向的取值，vp0表示所述待测储层中纵波速度在垂直方向的取值，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值，δ表示所述纵波波速度在近垂直入射时随相位角的变化情况。

s103：根据所述待测储层平面纵波的垂直慢度关系式、所述待测储层平面中横波的垂直慢度关系式，构建新的具有垂向对称轴的横向各向同性介质的反射系数计算方式。

s104：利用所述介质反射系数计算方式，对所述待测储层进行振幅随角度变化的反演。

在本申请中，可以基于所述待测储层中的具有垂向对称轴的横向各向同性介质的参数(即vti参数)具有各向异性的特征，采用如下方式对所述待测储层进行ava反演：

在本申请中，所述vti参数可以包括：三个弹性参数和两个各向异性参数。其中，弹性参数可以是：待测储层在垂直方向的纵波速度vp、待测储层在垂直方向的横波速度vs和密度ρ；各向异性参数可以是：横向和垂向p波速度的差异结果ε和p波速度在近垂直入射时随相位角的变化情况δ。当然，也可以包括地层厚度等其他参数，本申请对此不作限定。

s4-1：基于所述待测储层的地震数据提取角度依赖的子波，通过测井数据正演模拟和井旁道地震数据确定振幅缩放因子，并应用于所述子波；确定角度域的pp波与ps波道集。

基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。计算振幅缩放因子，并应用于上述子波，达到模拟记录与实际记录的振幅匹配。其中，当地震数据信噪比较高时，为道集中的每一个地震道使用统一的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比较低时，可分近、中、远偏移距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的影响。

由于单独利用地震数据中的pp波进行反演时，存在较大不确定性，而ps波数据包含丰富的密度和各向异性参数信息，因此，在本申请中通过利用pp波和ps波一起进行ava反演，从而提高反演的稳定性和精度。在本申请中，可以通过结合测井数据以及精确的vti介质反射系数方程正演模拟确定角度域的pp波与ps波道集。

s4-2：利用地震构造解释资料和测井数据，建立初始弹性参数模型。

可以采用三维空间插值法建立弹性参数的模型。首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，从而计算得到地下每个点的弹性参数值，并建立初始弹性参数模型。

s4-3：从vti介质的精确反射系数公式出发，基于弱各向异性假设，略去高阶无穷小量，推导得到新的所述vti介质的反射系数计算方式。基于所述初始弹性参数模型和所述vti介质的反射系数计算方式正演模拟角度域的多波地震叠前道集，并由正演模拟记录和实际记录直接计算pp波与ps波反演残差。

由于本申请所采用的ava反演算法需要求取包括5个vti参数在内的十个未知弹性参数和各向异性参数的一阶偏导数，由于精确的vti反射系数求解表达式非常复杂，因此这个一阶偏导数求解起来非常复杂。通过观察发现，所述待测储层平面p波的垂直慢度公式、所述待测储层平面sv波的垂直慢度公式，这两个慢度复杂的表达形式直接导致了完整表达式的复杂结构，因此，本申请提出了在计算精度损失较小的前提下，对这两个慢度的表达式进行简化，使最终一阶偏导数的求解变得相对简单和稳定。

在精确的vti反射系数方程中，不同波模式的垂直慢度可表示为：

上式中，p＝sin(θ)/vp(θ)；

上式中，qα表示所述待测储层平面p波的垂直慢度，qβ表示所述待测储层平面sv波的垂直慢度公式，cij表示刚度，ρ表示密度，p表示水平慢度，vp(θ)表示p波相速度，θ表示相位角。

thomsen引入了一组适用于vti介质的新的各向异性参数集，即，可以利用两个垂向的速度和三个无量纲的各向异性参数对五个独立的刚度进行替换：

上式中，vp0表示所述待测储层中纵波速度在垂直方向的取值，vs0表示所述待测储层中横波速度在垂直方向的取值，γ表示横向和垂向sh波速度的差异结果，ε表示p波速度在横向和垂向的差异值，δ描述的是p波速度在近垂直入射时随相位角的变化情况。由于在vti介质中sh波从p波和sv波中完全解耦，所以在本申请中只考虑p波和sv波的反射和透射系数。

基于弱各向异性假设，vti介质中的p波相速度vp(θ)和刚度c13可以利用各向异性参数ε和δ线性表示为：

vp(θ)≈vp0(1+δsin²θcos²θ+εsin⁴θ)(3)

c13≈c33(1+δ)-2c55(4)

根据等式(2)，我们可以得到形式相对简单的水平慢度：

将等式(2)和(3)代入等式上述k1、k2、k3的表达式中，并利用thomsen各向异性参数和垂向的速度对独立的刚度进行替换，同时略去高阶项，可得：

将k2和k3代入k1中，可以得到：

略去高阶项无穷小量，可以得到：

令，

则有

将等式(8)和(9)代入等式(1)中，可以得到：

这样，可以得到了形式简单的所述待测储层平面p波的垂直慢度公式、所述待测储层平面sv波的垂直慢度公式。

从现有的vti介质的精确反射系数公式中，选取包含本申请所提出的平面p波的垂直慢度、平面sv波的垂直慢度的所述vti介质的精确反射系数公式，并将形式简单的所述待测储层平面p波的垂直慢度公式、所述待测储层平面sv波的垂直慢度公式，即公式(10)代入选取出的所述vti介质的精确反射系数公式中，从而可以构建出新的vti介质的反射系数计算方式。

下式为所述vti介质的精确反射系数公式，将简化过的水平慢度和垂直慢度代入到下式中，就可以实现对整个方程进行近似简化。

上式中，

上式中，rpp表示pp波反射系数，rps表示ps波反射系数，tpp表示pp波透射系数，tps表示ps波透射系数，(1)表示上层介质，(2)表示下层介质。

再基于所述初始弹性参数模型和所述vti介质的反射系数计算方式，结合上述应用了所述振幅缩放因子的子波，正演模拟角度域的多波地震叠前道集。根据上述正演模拟记录和实际记录直接计算得到pp波与ps波反演残差。

之后，可以利用简化后的所述vti介质的反射系数计算方式，结合所示pp波与ps波反演残差进行ava反演，从而确定vti参数，即：弹性参数和各向异性参数。

s4-4：基于所述待测储层的测井数据确定vti参数的先验分布函数。

基于所述测井数据得到所述待测储层某些位置处的vti参数及其均值，包括：弹性参数和各向异性参数。假设vti参数先验模型服从高斯分布的同时，还包含服从微分拉普拉斯分布的垂向块约束项，通过测井数据分析获得需要的模型参数，并计算确定五个vti参数的方差和协方差，构建协方差矩阵，从而可以形成符合该工区的模型参数先验分布函数r(m)。可以表示为：

上式中，cm为包含弹性参数和各向异性参数相关性的块对角矩阵，μ为vti参数的均值向量(不同的vti参数需要分别求取)，d为一阶微分算子，kl为尺度参数，l＝1,2,3,4,5，m是一个由五个参数向量排成一列组成的待估参数向量。进一步地，上述这些参数对于不同的vti参数可以选取不同的值。

s4-5：利用广义线性反演思想将目标函数改写成关于弹性参数扰动量的函数，然后对扰动量求导并令导数等于零得到扰动量的迭代求解公式，利用迭代重加权最小二乘算法对扰动量进行求解，并利用求解得到的扰动量对模型参数进行更新，重复以上工作步骤直至反演残差达到要求或达到最大迭代次数，输出最终得到的vti参数。

根据bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数。

设vti参数为m^t＝(m1,m2,…,mn)^t，观测地震数据为d^t＝(d1,d2,…,dn)^t，由贝叶斯理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的问题可以归结为求解一个后验概率函数：

上式中，p(d)＝∫p(d|m)p(m)dm为归一化因子，可以看作是常数，p(d|m)为似然函数，p(m)为先验概率分布。

根据上述后验概率确定反演的目标函数j1(m)，具体的，可以表示如下：

上式中，dpp表示实际的pp波观测地震数据，dps表示实际的ps波观测地震数据，gpp(m)表示pp波正演地震记录，gps(m)表示ps波正演地震记录，表示pp波噪声的协方差逆矩阵，表示ps波噪声的协方差逆矩阵，r(m)表示先验约束项。

由于该模型参数不好直接求解，因此可以借助于泰勒展开对上述概率函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量δm的迭代求解公式。

δm^k＝(h^k)^-1γ^k

上式中，

上式中，和表示正演方程关于界面上下层十个参数的一阶偏导数gpp(m^k)和gps(m^k)分别表示第k次迭代之后的正演地震记录。

将所述vti介质反射系数计算方式代入所述迭代求解公式，其中，所述vti介质反射系数计算方程关于vti参数的一阶偏导数可以通过解析或数值方法求得。根据推导，可以得到模型参数的更新迭代公式：

m^k+1＝m^k+η^kδm^k,k＝0,1,2,...

其中，η^k是第k次迭代的步长因子，m表示由所有vti参数组成的参数向量。在本申请中所述步长因子η^k＝1，k＝4。

根据模型参数扰动量的迭代求解公式和由s4-3所确定的pp波与ps波反演残差，采用迭代重加权最小二乘算法计算模型参数的扰动量，重复迭代上述步骤，将本次求解得到的扰动量加在所述初始弹性参数模型上，就可以得到新的vti参数，将这个新的vti参数看作新的初始模型，根据模型参数扰动量的求解表达式又可以求解得到新的扰动量，就这样重复更新和迭代求解，直至最终达到迭代终止条件：由s4-3所确定的pp波与ps波反演残差，而最终得到的vti参数就是最终的反演结果。

如图2所示为本申请提供的一种地震反演方法反演得到的纵波速度vp(a)、横波速度vs(b)、密度ρ(c)、横向和垂向p波速度的差异结果ε(d)和p波速度在近垂直入射时随相位角的变化情况δ(e)，图中纵轴表示时间，单位为：秒，横轴自左至右，从图(a)至图(e)分别表示纵波速度(单位：km/s)、横波速度(单位：km/s)、密度(单位：g/cm3)、ε和δ。基于新近似公式的页岩储层vti参数多波ava反演方法能够较高精度的预测到vti介质的弹性参数和各向异性参数信息，由于引入了包含密度信息的先验分布，并采用ps横波进行联合反演，其对密度模型也预测准确。

基于同一发明构思，本发明实施例中还提供了一种地震反演装置，如下面的实施例所述。由于地震反演装置解决问题的原理与地震反演方法相似，因此地震反演装置的实施可以参见地震反演方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。图3是本发明实施例的地震反演装置的一种结构框图，如图3所示，可以包括：p波确定模块301、sv波确定模块302、vti系数确定模块303、ava反演模块304，下面对该结构进行说明。

纵波确定模块301，可以用于基于待测储层中纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面纵波的垂直慢度关系式；

横波确定模块302，可以用于根据待测储层中横波速度在垂直方向的取值、所述纵波速度在垂直方向的取值、所述纵波速度在近垂直入射时随相位角的变化情况、水平慢度、所述纵波速度在横向和垂向的差异值，确定所述待测储层平面横波的垂直慢度关系式；

反射系数计算模块303，可以用于根据所述待测储层平面纵波的垂直慢度公式、所述待测储层平面中横波的垂直慢度公式，构建新的具有垂向对称轴的横向各向同性介质的反射系数计算方式；

反演模块304，可以用于利用所述反射系数计算方式，对所述待测储层进行振幅随角度变化的反演。

在一个实施例中，所述纵波确定模块可以包括：第一储层参数计算单元，可以用于计算所述待测储层中所述纵波速度在垂直方向的取值，并计算所述取值的倒数的平方，作为第一储层参数；第二储层参数计算单元，可以用于计算所述水平慢度的平方，并将所述水平慢度的平方与所述纵波速度在近垂直入射时随相位角的变化情况相乘，作为第二储层参数；第三储层参数计算单元，可以用于计算所述纵波速度在横向和垂向的差异值、所述纵波速度在近垂直入射时随相位角的变化情况之差，作为速度差异值，计算所述纵波速度在垂直方向的取值的平方，计算所述水平慢度的四次方，计算所述速度差异值、所述取值的平方、所述水平慢度的四次方之间的差值，并将所得到的差值作为第三储层参数；纵波的垂直慢度计算单元，可以用于计算所述第三储层参数的方根，并将所述方根作为所述待测储层平面纵波的垂直慢度。

在一个实施例中，所述纵波确定模块具体可以用于按照以下公式确定所述待测储层平面纵波的垂直慢度：

上式中，qα表示所述待测储层平面纵波的垂直慢度，vp0表示所述纵波速度在垂直方向的取值，δ表示所述纵波速度在近垂直入射时随相位角的变化情况，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值。

在一个实施例中，所述横波确定模块可以包括：第四储层参数计算单元，可以用于计算所述待测储层中横波速度在垂直方向的取值的倒数，并计算所述倒数的平方，作为第四储层参数；第五储层参数计算单元，可以用于计算所述纵波速度在垂直方向的取值的平方与所述水平慢度的四次方之间的乘积，得到第五储层参数；第六储层参数计算单元，可以用于计算所述纵波速度在垂直方向的取值与所述横波速度在垂直方向的取值的比值的平方，并计算所述比值的平方与所述水平慢度的平方的乘积作为水平值，将所述比值的平方与所述水平值相乘，得到第六储层参数；第七储层参数计算单元，可以用于计算所述第五储层参数和所述第六储层参数的差值，作为第七储层参数；第八储层参数计算单元，可以用于计算所述速度差异值和所述第七储层参数的乘积，并将所得到的乘积结果和所述第四储层参数相加，将相加后的结果作为第八储层参数；横波的垂直慢度计算单元，可以用于计算所述第八储层参数、所述水平慢度的平方之间的差值，将所得到的差值作为所述待测储层平面横波的垂直慢度。

在一个实施例中，所述横波确定模块具体可以用于按照以下公式确定所述待测储层平面横波的垂直慢度：

上式中，qβ表示所述待测储层平面横波的垂直慢度，vs0表示所述横波速度在垂直方向的取值，vp0表示所述待测储层中纵波速度在垂直方向的取值，p表示所述水平慢度，ε表示所述纵波速度在横向和垂向的差异值，δ表示所述纵波波速度在近垂直入射时随相位角的变化情况。

尽管本申请内容中提到平面p波的垂直慢度确定方式、平面sv波的垂直慢度确定方式以及ava反演方式等描述，但是，本申请并不局限于必须是本申请实施例所描述的情况。某些行业标准或者使用自定义方式或实施例描述的实施基础上略加修改后的实施方案也可以实现上述实施例相同、等同或相近、或变形后可预料的实施效果。应用这些修改或变形后的平面p波的垂直慢度确定方式、平面sv波的垂直慢度确定方式以及ava反演方式等获取的实施例，仍然可以属于本申请的可选实施方案范围之内。

虽然本申请提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或终端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。

上述实施例阐明的单元、装置或模块等，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申请时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块或子单元的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

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

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

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

本说明书中的各个实施例采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。本申请可用于众多通用或专用的计算机系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的电子设备、网络pc、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。

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