CN109188511A - 一种砂泥岩薄互层介质多波avo联合反演方法 - Google Patents

Abstract

本发明公开了一种砂泥岩薄互层介质多波AVO联合反演方法，针对复杂薄互层介质，利用Ruger提出的VTI介质反射系数方程，首先利用叠前AVA角道集反演出等效VTI介质的5个参数，然后进一步利用VTI介质参数与其内部单层参数和N/G的关系反演出薄互层内部砂岩和泥岩的参数以及互层结构N/G，结果表明，无论是对VTI介质还是薄互层的内部参数，反演均能取得令人满意的效果。

Description

一种砂泥岩薄互层介质多波AVO联合反演方法

技术领域

本发明涉及地震勘探技术领域，具体涉及一种砂泥岩薄互层介质多波AVO联合反演方法。

背景技术

随着油气勘探的不断深入，简单的构造型油气藏已经开发殆尽，未来的勘探重点是隐蔽的岩性油气藏。在我国东部和西部的一些盆地，绝大多数的中、新生代含油气盆地均以薄层砂岩泥岩沉积为主^[1]。薄互层储层内部砂体通常呈条带状，数量多，厚度远远小于地震的勘探分辨率，并且岩性和厚度横向上变化快，针对薄互层储层的勘探存在很大的问题。由于地震频带内难以区分薄互层，地震波的波长远远大于薄互层，因此国际上常把薄互层整体等效为各向异性VTI(Vertical Transversely Isotropic)介质进行研究^[2-4]。薄互层的互层结构可以用N/G(net-to-gross)来表示，其中N表示砂岩的总厚度，G表示薄互层介质的总厚度^[5-6]。Thomsen用P波和SV波的速度以及两个各向异性参数来描述VTI介质的性质^[7]，所以薄互层可以通过P波和SV波的速度、两个各向异性参数和N/G来描述。

AVO(Amplitude versus offset)技术是目前进行油气预测的主要方式，它通过分析振幅随着炮检距的变化能反映储层的岩性和流体信息。AVO分析的理论技术是Zoeppritz方程，它精确的反映了反射波和透射波与入射角的关系。为了研究各向异性介质的AVO，许多学者在不同的假设条件下提出了不同的公式。Banik把Thomsen参数用到横向各项同性介质的PP波的反射系数公式中^[8]。提出在弱各向异性和弱对比度界面下的近似位移和能量场的PP波和PS波反射/透射系数^[9]。Véronique和推导了弱各向异性介质中任意界面对比的平面波位移反射/透射系数的近似公式^[10]。Rüger提出了有关VTI和HTI(垂向各向异性)介质的PP波和PS波反射系数公式^[11]。这些公式有更明确的物理意义，且形式上更加简单，对进行岩性参数反演非常有利。

目前对于AVO参数反演大都是基于各向同性理论或者P-P波单分量的反演。张丰麟等采用广义线性反演算法结合贝叶斯理论，使用精确Zoeppritz方程进行反演，使得反演结果更加准确且稳定^[12]。Lu等提出了基于精确Zoeppritz方程的纵横波联合反演方法，大大提高了反演精度且能适用于较大角度的反演^[13]。徐善辉等将遗传算法引入到TTI介质的各向异性参数反演中，利用反射系数反演出各向异性参数和对称轴倾角、方位角等信息^[14]。这样的反演方法无论在反演原理上还是精度上都无法适用于薄互层介质的反演且无法反演出薄互层的内部参数和互层结构N/G。

参考文献

[1]王赟,杨春,芦俊.薄互层弹性波反演面临的困境[J].地球物理学报,2018,61(03):1118–1135.

[2]Backus G E.Long-wave elastic anisotropy produced by horizontallayering[J].Journal of Geophysical Research,1962,67(11):4427–4440.

[3]Postma G W.Wave propagation in a stratified medium[J].Geophysics,1955,20(4):780–806.

[4]Sayers C M.Long-wave seismic anisotropy of heterogeneousreservoirs[J].Geophys.J.Int.,1998,132(132):667–673.

[5]A.Stovas,M.and P.Avseth,“AVO attribute inversion forfinely layered reservoirs,”Geophysics,vol.71,no.3,pp.C25–C36,2006.

[6]Vernik L,D Fisher,S Bahret.Estimation of net-to-gross from P and Simpedance in deepwater turbidite[J].The Leading Edge,2002,21(4):380–387.

[7]Thomsen L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954–1966.

[8]Banik,N.C.An effective anisotropy parameter in transverselyisotropic media[J],Geophysics,1987，52：1654–1664.

[9]V.Weak-contrast reflection/transmission coefficients inweakly anisotropic elastic media:P-wave incidence[J],Geophys.J.Int.,1999，138：553–562.

[10]Véronique,F.&I.First-order reflection/transmissioncoefficients for unconverted plane P waves in weakly anisotropic media[J],Geophys.J.Int.,2010，183：1443–1454.

[11]Rüger,A.Reflection coefficients and azimuthal AVO analysis inanisotropic media[D].1996,Ph.D.thesis,Colorado School of Mines.

[12]张丰麒，魏福吉，王彦春等.基于精确Zoeppritz方程三变量柯西分布先验约束的广义线性AVO反演[J].地球物理学报，2013，56(06)：2098–2115.

[13]Lu,J.,Z.Yang,Y.Wang,and Y.Shi,2015,Joint PP and PS AVA seismicinversion using exact Zoeppritz equations:Geophysics,80,no.5,R239–R250.

[14]徐善辉,韩立国,郭建.TTI介质各向异性参数多波反演与PS波AVO分析[J].地球物理学报,2012,55(02)：569–576.

[15]Levenberg,K.,1944,A Method for the Solution of Certain Non-LinearProblems in Least Squares:Quarterly of Applied Mathematics,2,164–168.

[16]Marquardt,D.W.,1963,An Algorithm for Least-Squares Estimation ofNonlinear Inequalities:Journal of the Society for Industrial and AppliedMathematics,11,431–441.

发明内容

针对现有技术的不足，本发明旨在提供一种一种砂泥岩薄互层介质多波AVO联合反演方法，无论是对VTI介质还是薄互层的内部参数，反演均能取得令人满意的效果。

为了实现上述目的，本发明采用如下技术方案：

一种砂泥岩薄互层介质多波AVO联合反演方法，包括如下步骤：

S1、抽取PP波AVA道集和PS波AVA道集，并将PS波资料压缩到PP波的旅行时上，使PS波与PP波时间剖面相匹配；

S2、建立VTI介质参数反演的初始模型：

2.1)利用测井资料建立纵横速度和密度的线性初始模型；

2.2)根据岩石物理测井信息或者速度分析，得到初始各向异性参数模型；

S3、建立薄互层内部参数反演的初始模型：

3.1)根据测井资料得到的砂泥岩纵横波速度、密度以及厚度资料；

3.2)利用步骤3.1)得到的资料建立砂泥岩纵横波速度、密度和N/G的线性初始模型；

S4、多波AVA反演

S4.1、VTI介质参数反演：

4.1.1)将步骤S2中建立得到的纵横速度和密度的线性初始模型以及初始各向异性参数模型转化为以Δε,Δδ,为模型参数如下：

M＝[A₁,B₁,C₁,D₁,E₁,…A_m,B_m,C_m,D_m,E_m]^T,

其中，D＝Δε,E＝Δδ,；m表示反演剖面中每一道的采样点数，则模型参数的维数为5*m；记初始模型参数为M₀；

4.1.2)利用式(1)计算模型参数的反射系数序列，并与提取出的子波进行褶积，得到PP波和PS波的合成记录；

其中，

V_P0和V_S0分别表示P波和SV波的垂直速度，ρ是密度，ε和δ是VTI介质的各向异性参数，i和j分别为P波的入射角和透射角；下标1和2分别表示反射界面上下层的弹性参数；

4.1.3)根据式(2)、(3)计算出模型增量ΔM和目标函数Q₁，如果Q₁不满足精度要求，则ΔM叠加到初始参数模型上得到新的模型参数，并返回步骤4.1.2)；直到Q₁满足精度要求，迭代停止；此时得到一个达到目标分辨率的纵波速度、横波速度、密度和各向异性参数的数据体；

ΔM＝(J^TJ)^-1J^T[Y(M,θ,t)-Y(M₀,θ,t)], (2)

Q₁＝||(J^TJ+λI)·ΔM-J^T[Y(M,θ,t)-Y(M_inv,θ,t)]||². (3)

设实测的地震记录Y(M,θ,t)＝[y(M,θ₁,t),y(M,θ₂,t),…,y(M,θ_n,t)]^T，其中n表示CMP的道数，θ表示入射角，t表示时间，地震记录等于反射系数序列和子波的褶积，即：

Y＝ω(t)*R(θ),

式中，ω(t)是地震子波的序列，R(θ)是反射系数序列，具体表示为：

对地震记录Y在M₀处进行Taylor展开并略去高阶项，得到：

ΔY＝JΔM,

其中，ΔM＝M-M₀代表模型参数的增量，ΔY＝Y(M,θ,t)-Y(M₀,θ,t)为实测数据地震记录与初始模型参数地震记录计算得到的地震记录之差；J是Jacobian矩阵，是地震记录Y(M,θ,t)对于模型参数M的一阶导数；J表示为：

S4.2、薄互层内部参数反演：

4.2.1)薄互层内部参数的模型参数定义为：

P＝[V_Pa1,V_Sa1,ρ_a1,V_Pb1,V_Sb1,ρ_b1,r₁,…V_Pax,V_Sax,ρ_ax,V_Pbx,V_Sbx,ρ_bx,r_x]；

r表示N/G；记薄互层内部参数的初始模型参数为P₀；

4.2.2)用初始模型参数为P₀计算得到VTI介质模型参数S：

S＝[V_P1,V_S1,ρ₁,ε₁,δ₁,…V_Px,V_Sx,ρ_x,ε_x,δ_x]^T；

式中x表示VTI介质的个数；

4.2.2)求解下式，得到模型增量ΔP和目标函数Q₂；如果Q₂不满足精度要求，则ΔP叠加到初始模型参数上得到新的模型参数，并用新的模型参数计算VTI介质模型参数S，并再次计算下式，直到Q₂满足精度要求，迭代停止；此时获得一个达到目标分辨率的砂泥岩纵波速度、横波速度和密度以及砂岩与储层厚度比N/G的数据体；

ΔS＝GΔP；

Q₂＝||(G^TG+λI)·ΔP-G^T(S-S_inv)||²；

其中，ΔP＝P-P₀代表模型参数的增量，ΔS＝S-S₀是实际VTI介质参数与通过初始模型参数计算得到的VTI参数之差；G是Jacobian矩阵，是VTI介质模型参数S对于模型参数P的一阶导数,G表示为：

需要说明的是，步骤S1的具体流程如下：

S1.1、抽取PP波AVA道集

1.1.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

1.1.2)确定反演所需要的入射角，抽取其振幅参数得到PP波的AVA道集；

S1.2、抽取PS波AVA道集

1.2.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

1.2.2)确定反演所需要的入射角，抽取其振幅参数得到PS波的AVA道集；

1.2.3)通过纵横波速度比利用下式将PS波资料压缩到PP波的旅行时上，使PS波与PP波时间剖面相匹配，便于多波层位的对比：

本发明的有益效果在于：本发明针对复杂薄互层介质，利用Ruger提出的VTI介质反射系数方程，首先利用叠前AVA角道集反演出等效VTI介质的5个参数，然后进一步利用VTI介质参数与其内部单层参数和N/G的关系反演出薄互层内部砂岩和泥岩的参数以及互层结构N/G，结果表明，无论是对VTI介质还是薄互层的内部参数，反演均能取得令人满意的效果。

具体实施方式

以下将对本发明作进一步的描述，需要说明的是，本实施例以本技术方案的前提，给出了详细的实施方式和具体的操作过程，但本发明的保护范围并不限于本实施例。

以下将对本实施例所提供的一种砂泥岩薄互层介质多波AVO联合反演方法的基本原理作具体描述。

一、VTI介质反射系数方程

Zoeppritz方程是基于各向同性理论，针对VTI介质，本实施例中把Ruger提出的各向异性方程^[11]重新整理，将其转化为对波阻抗差、密度差和各向异性参数差的反射系数方程如公式(1)所示，为后续进行VTI介质的参数反演做基础。

其中，

V_P0和V_S0分别表示P波和SV波的垂直速度，ρ是密度，ε和δ是VTI介质的各向异性参数，i和j分别为P波的入射角和透射角。下标1和2分别表示反射界面上下层的弹性参数。

二、VTI介质参数与N/G的关系

Thomsen提出可以用P波和SV波的垂直速度和各向异性参数来表示VTI介质^[7]，可表示为：

其中V_P0和V_S0分别表示P波和SV波的垂直速度，ε和δ是VTI介质的各向异性参数。c_ij(其中，i与j表示公式(3)中符号c对应的下标)表示VTI介质的弹性常数，它与N/G的关系可以表示为^[5]：

其中，VTI介质的密度可以表示为：

VTI介质中砂岩和泥岩的弹性参数差和密度差可以表示为：

下标a和b分别表示VTI内部单层砂岩和泥岩的参数，V_Pa、V_Sa和ρ_a分别为砂岩的纵波速度、横波速度和密度，V_Pb、V_Sb和ρ_b分别为泥岩的纵波速度、横波速度和密度，r是砂岩和泥岩平均横波速度和平均纵波速度的比值，则可以得到VTI介质参数与VTI介质内部砂泥岩参数和N/G的关系：

根据公式(6)中VTI介质参数与其内部砂泥岩参数和N/G的关系，可以通过VTI介质的参数进行反演，得到薄互层内部砂岩参数、泥岩参数和N/G的值。

三、多波AVA反演方法

3.1、VTI介质参数反演

由公式(1)的反射系数方程可以知道，VTI介质的AVA曲线主要取决于五个参数：

那么利用叠前5道或者更多道的纵波或者转换P-SV波就可以反演出这5个参数。当上层岩性参数已知时，可以进一步的得到V_P0，V_S0，ρ，ε和δ等5个参数。为了表示方便，令：

则反演的模型参数可以表示为：

M＝[A₁,B₁,C₁,D₁,E₁,…A_m,B_m,C_m,D_m,E_m]^T, (9)

其中m表示反演剖面中每一道的采样点数，则模型参量的维数为5*m。设实测的地震记录Y(M,θ,t)＝[y(M,θ₁,t),y(M,θ₂,t),…,y(M,θ_n,t)]^T，其中n表示CMP的道数，θ表示入射角，t表示时间。地震记录等于反射系数序列和子波的褶积，即：

Y＝ω(t)*R(θ), (10)

式中，ω(t)是地震子波的序列，R(θ)是反射系数序列，具体表示为：

假设初始模型参数为M₀，对地震记录Y在M₀处进行Taylor展开并略去高阶项，可以得到：

ΔY＝JΔM, (12)

其中，ΔM＝M-M₀代表模型参数的增量，ΔY＝Y(M,θ,t)-Y(M₀,θ,t)实测数据地震记录与初始模型参数地震记录计算得到的地震记录之差。J是Jacobian矩阵，是地震记录Y(M,θ,t)对于模型参数M的一阶导数。此次进行纵横波联合反演，因此Y(M,θ,t)中包含了P-P和P-SV两种道集信息，所以J可以表示为：

式中反射系数对于模型参数的偏导数可以通过中心差分来近似计算，以P波的反射系数对A的一阶导数为例，中心差分可以表示为:

3.2、薄互层内部参数反演方法

与VTI介质的参数反演类似，对于薄互层内部参数的反演，模型参数定义为：

P＝[V_Pa1,V_Sa1,ρ_a1,V_Pb1,V_Sb1,ρ_b1,r₁,…V_Pax,V_Sax,ρ_ax,V_Pbx,V_Sbx,ρ_bx,r_x] (15)

r表示N/G。VTI介质参数可以表示为：

S＝[V_P1,V_S1,ρ₁,ε₁,δ₁,…V_Px,V_Sx,ρ_x,ε_x,δ_x]^T (16)

式中x表示VTI介质的个数。假设初始模型参数为P₀，对VTI介质参数S在P₀处进行Taylor展开并略去高阶项，可以得到：

ΔS＝GΔP (17)

其中，ΔP＝P-P₀代表模型参数的增量，ΔS＝S-S₀是实际VTI介质参数与通过初始模型计算得到的VTI参数之差。G是Jacobian矩阵，是VTI介质参数S对于模型参数P的一阶导数,G可以表示为：

G中VTI介质参数对于模型参数的偏导可以根据公式(6)直接求出。以第k层VTI介质参数[V_Pk,V_Sk,ρ_k,ε_k,δ_k]对其内部参数[V_Pak,V_Sak,ρ_ak,V_Pbk,V_Sbk,ρ_bk,r_k]的求导为例，给出具体求导后的结果：

3.3、反演方程的求解

两步反演中反演方程的解法类似，以第一步VTI介质参数的反演方程为例，详细给出求解方法。公式(12)可转化为：

J^TΔY＝J^TJΔM, (20)

所以可以得到：

ΔM＝(J^TJ)^-1J^T[Y(M,θ,t)-Y(M₀,θ,t)], (21)

由于J^TJ有时会出现不可逆的情况，此时利用LM算法[14-15]将公式(21)转化为：

ΔM＝(J^TJ+λI)^-1J^T[Y(M,θ,t)-Y(M₀,θ,t)], (22)

其中，I是单位矩阵，λ是一个正实数。此时，可以构建出反演的目标函数：

Q₁＝||(J^TJ+λI)·ΔM-J^T[Y(M,θ,t)-Y(M_inv,θ,t)]||². (23)

在反演迭代中，模型增量ΔM不断的加到M上进行计算，式中Y(M_inv,θ,t)代表模型反演出的地震记录。当Y(M,θ,t)和Y(M_inv,θ,t)非常接近时，迭代停止，得到反演出的模型参数。

四、AVO转化为AVA

动校正之后的道集Y(M,x,t)是炮检距的函数，也就是AVO。在反演中用到的是AVA，因此需要把AVO转化为AVA。通常采用采用直射线法，根据炮检距X和纵横波速度V_P和V_S，可以得到P波的偏移距X_P。对于P-P波，对于P-SV波，然后根据目的层的深度H和炮检距X，根据公式(24)可以计算出P-P波和P-SV波的入射角θ。

θ＝arctan(X_P/H). (24)

五、多波层位匹配

通常P-SV波的速度小，旅行时长，所以P-SV波的剖面在纵向上大于P-P波的剖面。进行多波联合反演的关键是多波的标定和压缩匹配。对于P-SV的压缩主要靠纵横波的联合标定和精细解释来实现，这其中最关键的是确定横波和纵波的速度γ，可以通过标志层来求取主要反射界面的速度比，并进一步平滑和差值得到整个工区的速度比剖面。利用地震反射标志层作为约束，根据公式(25)将P-SV的时间剖面进行压缩，使其与P-P波时间剖面相匹配。

六、提取多波地震子波

通常地震子波的提取方法主要有两大类：一种确定性子波，一种是统计性子波。确定性子波是利用测井资料提取的，这使得它的应用很有局限性，而统计性子波不需要测井资料也可以得到比较准确的子波。因此在进行多波层为匹配之后，用多道统计方法提取纵波子波和转换波子波。由于干扰和反射系数序列不完全非相关性的影响，需要对子波的振幅谱和相位谱做整形处理。

本实施例提供一种砂泥岩薄互层介质多波AVO联合反演方法，具体流程包括如下步骤：

一、抽取AVA道集

1、抽取PP波AVA道集

1.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

1.2)确定反演所需要的入射角，抽取其振幅参数得到PP波的AVA道集；

2、抽取PS波AVA道集

2.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

2.2)确定反演所需要的入射角，抽取其振幅参数得到PS波的AVA道集；

2.3)通过纵横波速度比利用下式将PS波资料压缩到PP波的旅行时上，使PS波与PP波时间剖面相匹配，便于多波层位的对比：

二、建立初始模型

1、VTI介质参数反演

1.1)利用测井资料建立高分辨率的纵横速度和密度的线性初始模型。

1.2)根据岩石物理测井信息或者速度分析，得到初始各向异性参数模型。

2、薄互层内部参数反演

2.1)根据测井资料得到更精细的砂泥岩纵横波速度、密度以及厚度资料。

2.2)利用上述资料建立砂泥岩纵横波速度、密度和N/G的线性模型。

三、多波AVA反演

1、VTI介质参数反演

1.1)将初始模型转化为以Δε,Δδ,为参数的初始模型。

1.2)利用公式(1)计算初始模型的反射系数序列，并与提取出的子波进行褶积，得到PP波和PS波的合成记录。

1.3)根据公式(22)计算出模型增量ΔM和目标函数Q₁。如果Q₁不满足精度要求，则ΔM叠加到初始模型上得到新的模型参数，并返回2)。直到满足条件，迭代停止。

此时可以得到一个达到目标分辨率的纵、横波速度、密度和各向异性参数的数据体。

2、薄互层内部参数反演

2.1)用初始模型计算得到VTI参数序列S。

2.2)类比于公式(20)–(23)求解公式(17)，得到模型增量ΔP和目标函数Q₂。如果Q₂不满足精度要求，则ΔP叠加到初始模型上得到新的模型参数，并返回1)。直到满足条件，迭代停止。

此时获得一个达到目标分辨率的砂泥岩纵、横波速度和密度以及砂岩与储层厚度比N/G的数据体。

以下将通过测试对本实施例方法的性能作进一步的描述。

通过建立薄互层模型来验证本实施例的反演方法，模型参数如表1所示。通过公式(2)–(5)可以计算出该薄互层地层的等效VTI参数如表2所示。

表1砂泥岩薄互层地质模型参数

表2等效VTI介质参数

1、VTI介质参数的反演

反演的第一步是针对VTI介质参数的反演。根据表1的模型数据，采用起始角度为0，终止角度为60度生成反射系数进行反演，可以反演出Δε,Δδ这五个参数。当已知上层VTI介质参数时，可以算出下层VTI介质的参数如表3所示，可以看出反演结果比较精确。

表3 VTI介质的反演结果

二、薄互层内部参数反演

通过第一步的反演得到VTI介质参数后，可以进一步反演VTI介质的内部参数，反演结果如表4所示。可以看出，反演值与真实值比较接近，反演效果较好。

表4薄互层内部参数的反演结果

从上述数据可知，本实施例提供的反演方法，不仅能反演出薄互层地层的等效VTI参数，还能进一步反演出薄互层内部砂岩和泥岩的参数以及薄互层的互层结构N/G。通过模型测试，可以看出本实施例的反演方法具有较好的反演效果，能反演出比较精确的参数值，为薄互层的勘探和储层预测提供更精确的信息。

对于本领域的技术人员来说，可以根据以上的技术方案和构思，作出各种相应的改变和变形，而所有的这些改变和变形都应该包括在本发明权利要求的保护范围之内。

Claims (2)

1.一种砂泥岩薄互层介质多波AVO联合反演方法，其特征在于，包括如下步骤：

S1、抽取PP波AVA道集和PS波AVA道集，并将PS波资料压缩到PP波的旅行时上，使PS波与PP波时间剖面相匹配；

S2、建立VTI介质参数反演的初始模型：

2.1)利用测井资料建立纵横速度和密度的线性初始模型；

2.2)根据岩石物理测井信息或者速度分析，得到初始各向异性参数模型；

S3、建立薄互层内部参数反演的初始模型：

3.1)根据测井资料得到的砂泥岩纵横波速度、密度以及厚度资料；

3.2)利用步骤3.1)得到的资料建立砂泥岩纵横波速度、密度和N/G的线性初始模型；

S4、多波AVA反演

S4.1、VTI介质参数反演：

4.1.1)将步骤S2中建立得到的纵横速度和密度的线性初始模型以及初始各向异性参数模型转化为以Δε,Δδ,为模型参数如下：

M＝[A₁,B₁,C₁,D₁,E₁,…A_m,B_m,C_m,D_m,E_m]^T,

其中，D＝Δε,E＝Δδ,；m表示反演剖面中每一道的采样点数，则模型参数的维数为5*m；记初始模型参数为M₀；

4.1.2)利用式(1)计算模型参数的反射系数序列，并与提取出的子波进行褶积，得到PP波和PS波的合成记录；

其中，

V_P0和V_S0分别表示P波和SV波的垂直速度，ρ是密度，ε和δ是VTI介质的各向异性参数，i和j分别为P波的入射角和透射角；下标1和2分别表示反射界面上下层的弹性参数；

4.1.3)根据式(2)、(3)计算出模型增量ΔM和目标函数Q₁，如果Q₁不满足精度要求，则ΔM叠加到初始参数模型上得到新的模型参数，并返回步骤4.1.2)；直到Q₁满足精度要求，迭代停止；此时得到一个达到目标分辨率的纵波速度、横波速度、密度和各向异性参数的数据体；

ΔM＝(J^TJ)^-1J^T[Y(M,θ,t)-Y(M₀,θ,t)], (2)

Q₁＝||(J^TJ+λI)·ΔM-J^T[Y(M,θ,t)-Y(M_inv,θ,t)]||². (3)

设实测的地震记录Y(M,θ,t)＝[y(M,θ₁,t),y(M,θ₂,t),…,y(M,θ_n,t)]^T，其中n表示CMP的道数，θ表示入射角，t表示时间，地震记录等于反射系数序列和子波的褶积，即：

Y＝ω(t)*R(θ),

式中，ω(t)是地震子波的序列，R(θ)是反射系数序列，具体表示为：

对地震记录Y在M₀处进行Taylor展开并略去高阶项，得到：

ΔY＝JΔM,

其中，ΔM＝M-M₀代表模型参数的增量，ΔY＝Y(M,θ,t)-Y(M₀,θ,t)为实测数据地震记录与初始模型参数地震记录计算得到的地震记录之差；J是Jacobian矩阵，是地震记录Y(M,θ,t)对于模型参数M的一阶导数；J表示为：

S4.2、薄互层内部参数反演：

4.2.1)薄互层内部参数的模型参数定义为：

P＝[V_Pa1,V_Sa1,ρ_a1,V_Pb1,V_Sb1,ρ_b1,r₁,…V_Pax,V_Sax,ρ_ax,V_Pbx,V_Sbx,ρ_bx,r_x]；

r表示N/G；记薄互层内部参数的初始模型参数为P₀；

4.2.2)用初始模型参数为P₀计算得到VTI介质模型参数S：

S＝[V_P1,V_S1,ρ₁,ε₁,δ₁,…V_Px,V_Sx,ρ_x,ε_x,δ_x]^T；

式中x表示VTI介质的个数；

4.2.2)求解下式，得到模型增量ΔP和目标函数Q₂；如果Q₂不满足精度要求，则ΔP叠加到初始模型参数上得到新的模型参数，并用新的模型参数计算VTI介质模型参数S，并再次计算下式，直到Q₂满足精度要求，迭代停止；此时获得一个达到目标分辨率的砂泥岩纵波速度、横波速度和密度以及砂岩与储层厚度比N/G的数据体；

ΔS＝GΔP；

Q₂＝||(G^TG+λI)·ΔP-G^T(S-S_inv)||²；

其中，ΔP＝P-P₀代表模型参数的增量，ΔS＝S-S₀是实际VTI介质参数与通过初始模型参数计算得到的VTI参数之差；G是Jacobian矩阵，是VTI介质模型参数S对于模型参数P的一阶导数,G表示为：

2.根据权利要求1所述的砂泥岩薄互层介质多波AVO联合反演方法，其特征在于，步骤S1的具体流程如下：

S1.1、抽取PP波AVA道集

1.1.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

1.1.2)确定反演所需要的入射角，抽取其振幅参数得到PP波的AVA道集；

S1.2、抽取PS波AVA道集

1.2.1)根据地震道的道头信息，确定目的层的深度H和不同炮检距X，得到对应的入射角θ；

1.2.2)确定反演所需要的入射角，抽取其振幅参数得到PS波的AVA道集；

1.2.3)通过纵横波速度比利用下式将PS波资料压缩到PP波的旅行时上，使PS波与PP波时间剖面相匹配，便于多波层位的对比：