CN103984018A - 多波联合振幅随入射角变化反演方法 - Google Patents

Abstract

本发明适用于地震波技术领域，提供了一种多波联合振幅随入射角变化反演方法，适于地震波的多波叠前联合反演。所述方法包括：接收PP波与PS波的转换波道集，并提取PP波及压缩后的PS波的不同入射角的子波；根据不同入射角的子波，建立第一与第二矩阵；根据薄互层初始模型，建立第一与第二雅各比矩阵及第一与第二反射系数矩阵；根据所述PP波与PS波的转换波道集，建立第三与第四矩阵，并计算第一与第二目标函数；根据第一、第二、第三与第四矩阵、第一与第二雅各比矩阵、第一与第二反射系数矩阵，得到修正值，并修正薄互层初始模型，且调整所述目标函数。通过应用该方法，解决了现有技术存在的多波叠前反演方法会存在较大的误差的问题。

Description

多波联合振幅随入射角变化反演方法

技术领域

本发明涉及地震波技术领域，尤其涉及一种多波联合振幅随入射角变化反演方法。

背景技术

目前从许多的文献来看，多波振幅随入射角变化(Amplitude VariousAngle,AVA)叠前联合反演方法一般都是针对油气田的单界面反演，采用了阿基-理查德(Aki-Richards)近似公式以及与之相关的改造公式来模拟单界面的反射系数。这类反演方法的特点是假设界面上下的介质的波阻抗是弱反差的。

但对于薄互层来说，薄互层的波阻抗一般与围岩相差太大，并且薄互层是薄层，对于地震波频带来说，薄互层顶、底板的反射是无法区分的，这也使得上述多波叠前AVA联合反演方法用于薄互层会存在较大的误差。

发明内容

本发明的主要目的在于提供一种多波联合振幅随入射角变化反演方法，以解决现有技术存在的多波叠前AVA联合反演方法用于薄互层会存在较大的误差的问题。

为解决上述问题，本发明实施例提供一种目标函数的计算方法，适于地震波的多波叠前联合反演。所述方法包括：1)接收PP波的转换波道集与PS波的转换波道集，并提取所述PP波的不同入射角的子波以及压缩后的所述PS波的不同入射角的子波，其中压缩后的所述PS波为压缩所述PS波至PP波的预设时；2)根据所述不同入射角的子波，建立第一矩阵与第二矩阵；根据一薄互层初始模型，建立第一雅各比矩阵与第二雅各比矩阵及对应所述薄互层初始模型的第一反射系数矩阵与第二反射系数矩阵；3)根据所述PP波的转换波道集与PS波的转换波道集，建立第三矩阵与第四矩阵，并计算第一目标函数与第二目标函数；4)根据所述第一矩阵、所述第二矩阵、所述第三矩阵、所述第四矩阵、所述第一雅各比矩阵、所述第二雅各比矩阵、所述第一反射系数矩阵与所述第二反射系数矩阵，得到修正值，并修正所述薄互层初始模型，且调整所述第一目标函数与所述第二目标函数。

其中，所述第一矩阵、所述第二矩阵、所述第三矩阵、所述第四矩阵、所述第一雅各比矩阵、所述第二雅各比矩阵、所述第一反射系数矩阵、所述第二反射系数矩阵与所述修正值组合成如下公式： $\begin{array}{c}[{(W_{\mathrm{pp}}*G_{\mathrm{pp}})}^T(W_{\mathrm{pp}}*G_{\mathrm{pp}})+{(W_{\mathrm{ps}}*G_{\mathrm{ps}})}^T(W_{\mathrm{ps}}*G_{\mathrm{ps}})]\mathrm{\ΔV}\\ ={(W_{\mathrm{pp}}*G_{\mathrm{pp}})}^T(D_{\mathrm{pp}}-W_{\mathrm{pp}}*R_{\mathrm{pp}0})+{(W_{\mathrm{ps}}*G_{\mathrm{ps}})}^T(D_{\mathrm{ps}}-W_{\mathrm{ps}}*R_{\mathrm{ps}0})\end{array}$ 其中，W_pp为所述第一矩阵，W_ps为所述第二矩阵，G_pp为所述第三矩阵，G_ps为所述第四矩阵，R_pp0为所述第一雅各比矩阵，R_ps0为所述第二雅各比矩阵，(D_pp-W_pp*R_pp0)为所述第一反射系数矩阵，R_ps0为所述第二反射系数矩阵，△V为所述修正值。

其中，所述公式更加入阻尼矩阵；其中所述阻尼矩阵为λI，其中I为单位向量矩阵，λ取极小值。

其中，所述第一目标函数与第二目标函数分别为如下公式：

Q(V)＝‖S_pp-D_pp‖²+‖S_ps-D_ps‖²，

$T=\begin{array}{c}{\mathrm{\λ}}_1\·\left(\right|S_{\mathrm{pp}}-D_{\mathrm{pp}}|+|S_{\mathrm{ps}}-D_{\mathrm{ps}}\left|\right)\\ +{\mathrm{\λ}}_2\·\left(\right|v_{\mathrm{pi}}-v_{\mathrm{pw}}|+|v_{\mathrm{si}}-v_{\mathrm{sw}}|+{|\mathrm{\ρ}}_i-{\mathrm{\ρ}}_w|)\end{array},$

其中，D_pp、D_ps为实际地震记录，S_pp＝W_pp*R_pp和S_ps＝W_ps*R_ps为模型响应，i为反演得到的参数，w为测井模型外推得到的参数，λ₁、λ₂分别对应反演中地震道与测井模型约束作用的权重。

其中，所述调整所述第一目标函数与所述第二目标函数包括迭代所述步骤1)、2)、3)、4)至所述第一目标函数与所述第二目标函数的值最小为止。

其中，所述薄互层初始模型、所述PP波的转换波道集与PS波的转换波道集为标准地震数据记录格式。

其中，所述不同入射角的子波为文本格式。

根据本发明的技术方案，采用多波地震数据以及测井共同约束反演的结果，基于广义反演理论同时反演薄互层的纵波速度、横波速度以及密度。其中，PP波与PS波的反射系数采用策普里兹方程组精确求解。由于反射系数的求解过程无任何近似，所以本发明的目标函数的计算方法即使在无井约束的情况下也能对薄互层的储层参数精确反演。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是根据本发明实施例的模拟得到的PP波的AVA道集的示意图；

图2是根据本发明实施例的模拟得到的PS波的AVA道集的示意图；

图3是根据本发明实施例的多波联合AVA反演方法的流程图。

具体实施方式

本发明的主要思想在于，采用多波地震数据以及测井共同约束反演的结果，基于广义反演理论同时反演薄互层的纵波速度、横波速度以及密度。其中，PP波与PS波的反射系数采用策普里兹(Zoeppritz)方程组精确求解。由于反射系数的求解过程无任何近似，所以本发明的目标函数的计算方法即使在无井约束的情况下也能对薄互层的储层参数精确反演。

为使本发明的目的、技术方案和优点更加清楚，以下结合附图及具体实施例，对本发明作进一步地详细说明。

首先，在各向同性介质中，当一个平面纵波倾斜入射到两种介质分界面时，会产生四种波，即反射纵波、反射横波、透射纵波和透射横波，且满足斯涅尔定律(Snell's Law)，如公式(1.1)所示：

$\frac{V_{p1}}{\sin \mathrm{\α}}=\frac{V_{s1}}{\sin \mathrm{\β}}=\frac{V_{p2}}{\sin {\mathrm{\α}}^{\′}}=\frac{V_{s2}}{\sin {\mathrm{\β}}^{\′}}.---\left(1.1\right)$

根据边界条件解波动方程并引入纵波反射系数R_PP、转换横波反射系数R_PS、纵波透射系数T_PP和转换横波透射T_PS，可以得出四个波的位移振幅应当满足的策普里兹(Zoeppritz)方程组，如公式(1.2)所示：

$\left\{\begin{array}{c}\sin \mathrm{\α}{R}_{\mathrm{PP}}+\cos \mathrm{\β}{R}_{\mathrm{PS}}-\sin {\mathrm{\α}}^{\′}{T}_{\mathrm{PP}}+\cos {\mathrm{\β}}^{\′}{T}_{\mathrm{PS}}=-\sin \mathrm{\α},\\ \cos \mathrm{\α}{R}_{\mathrm{PP}}+\sin \mathrm{\β}{R}_{\mathrm{PS}}-\cos {\mathrm{\α}}^{\′}{T}_{\mathrm{PP}}+\sin {\mathrm{\β}}^{\′}{T}_{\mathrm{PS}}=\cos \mathrm{\α},\\ \cos 2\mathrm{\β}{R}_{\mathrm{PP}}-\frac{V_{S1}}{V_{P1}}\sin 2\mathrm{\β}{R}_{\mathrm{PS}}-\frac{{\mathrm{\ρ}}_2{V}_{P2}}{{\mathrm{\ρ}}_1{V}_{S2}}\cos 2{\mathrm{\β}}^{\′}{T}_{\mathrm{PP}}\frac{{\mathrm{\ρ}}_2{V}_{P2}}{{\mathrm{\ρ}}_1{V}_{S2}}\sin 2-{\mathrm{\β}}^{\′}{T}_{\mathrm{PS}}=-\cos 2\mathrm{\β}\\ \frac{V_{S2}^2}{V_{P1}}\sin 2\mathrm{\α}{R}_{\mathrm{PP}}+V_{S1}\cos 2\mathrm{\β}{R}_{\mathrm{PS}}+\frac{{\mathrm{\ρ}}_2{V}_{S2}^2}{{\mathrm{\ρ}}_1{V}_{P2}}\sin 2\mathrm{\α}{T}_{\mathrm{PP}}-\frac{{\mathrm{\ρ}}_2{V}_{S2}}{{\mathrm{\ρ}}_1}\cos {\mathrm{\β}}^{\′}2T_{\mathrm{PS}}=\frac{V_{S1}^2}{V_{P1}}\sin 2\mathrm{\α}\end{array}\right.---\left(1.2\right)$

在公式(1.2)中，V_P1和V_S1分别为上层介质纵、横波速度；V_P2和V_S2为下层介质的纵、横波速度；ρ₁和ρ₂为上、下层介质的密度；α和α′为纵波入射角和透射角；β和β′为横波反射角和透射角。对上述公式(1.2)求解则可得到P-P波、P-SV波的反射系数R_PP、R_PS。

根据公式(1.2)可以分析各向同性介质条件下薄互层的多波AVA地震响应特征。如表1所示为薄互层模型参数表，从浅至深的地层分别为：第四系地层、砂岩(Sand)、薄互层(Coal)、砂岩、砂岩，其中薄互层厚度(Thickness)为3米。

层	Vp(m/s)	Vs(m/s)	ρ(g/cc)	层(m)
					第四系地层	2800	933	2.28	100
砂岩	3200	1600	2.35	100
					薄互层	2000	1333	2.0	3
砂岩	3300	1650	2.35	80
					砂岩	3400	1700	2.52	70

表1薄互层模型参数表

图1是根据本发明实施例的模拟得到的PP波的AVA道集的示意图。图2是根据本发明实施例的模拟得到的PS波的AVA道集的示意图。如图1及2所示，其中PS波道集已经被压缩至PP波的T0时间，子波为30Hz的雷克子波。

从图1及2可以看出：由于薄互层为薄层，薄互层的PP波与PS波地震反射无法区分顶底板的位置，顶、底板反射形成了复合波形态；PP波与PS波的AVA特征有较大区别，对于第四系地层与砂岩界面来说，PP波随着入射角增大，反射能量减弱，但PS波反之；但对于薄互层反射来说，PP波与PS波的反射能量均随着入射角增大而增大；PP波与PS波的极性大多数相反。

另外，反演采用双目标函数分别为公式(1.3)及(1.4)，如下所示：

Q(V)＝‖S_pp-D_pp‖²+‖S_ps-D_ps‖²→minimum，                  (1.3)

$T=\begin{array}{c}{\mathrm{\λ}}_1\·\left(\right|S_{\mathrm{pp}}-D_{\mathrm{pp}}|+|S_{\mathrm{ps}}-D_{\mathrm{ps}}\left|\right)\\ +{\mathrm{\λ}}_2\·\left(\right|v_{\mathrm{pi}}-v_{\mathrm{pw}}|+|v_{\mathrm{si}}-v_{\mathrm{sw}}|+{|\mathrm{\ρ}}_i-{\mathrm{\ρ}}_w|)\end{array}\→\min \mathrm{imum},---\left(1.4\right)$

在公式(2.3)中，D_pp、D_ps为实际地震记录，S_pp＝W_pp*R_pp和S_ps＝W_ps*R_ps为模型响应；在公式(2.4)中，下标“i”表示反演得到的参数，下标“w”表示测井模型外推得到的参数；λ₁、λ₂分别对应反演中地震道与测井模型约束作用的权重，λ₁+λ₂＝1，当各取50％时，测井数据与地震数据在反演过程中的约束权重各占一半；当λ₁＝1时表明地震数据完全约束反演结果，测井数据不起约束作用；当λ₂＝1时，表明测井数据完全约束反演结果，地震数据不起约束作用，相当于将测井模型进行外推。若某一数据的可靠性较低，可以相应降低其约束的权重。

在公式(1.3)中，对纵波和转换波进行泰勒展开成公式(1.5)，如下所示：

R(V₀+△V₀)＝R(V₀)+G(V₀)△V₀，                   (1.5)

其中，V₀为纵波速度、横波速度和密度的初始猜测值，可表示为向量，如公式(1.6)所示：

(v_p v_s ρ)，                                 (1.6)

△V₀为修正量，如公式(1.7)所示：

(△v_p △v_s △ρ)，                          (1.7)

G(V₀)称为f(V₀)在V₀处的梯度，即雅各比矩阵。

反射系数线性近似公式代入反演目标(公式(1.3))，可得公式(1.8)，如下所示：

Q(V)＝‖W_pp*(R_pp0+G_pp△V)-D_pp‖²+‖W_ps*(R_ps0+G_ps△V)-D_ps‖²＝[W_pp*(R_pp0+G_pp△V)-D_pp]^T[W_pp*(R_pp0+G_pp△V)-D_pp]，    (1.8)+[W_ps*(R_ps0+G_ps△V)-D_ps]^T[W_ps*(R_ps0+G_ps△V)-D_ps]

令可得公式(1.9)，如下所示：

$\begin{array}{c}\∂Q\left(V\right)/\∂\mathrm{\ΔV}={(W_{\mathrm{pp}}*R_{\mathrm{pp}0}-D_{\mathrm{pp}})}^T(W_{\mathrm{pp}}*G_{\mathrm{pp}})+{(W_{\mathrm{pp}}*G_{\mathrm{pp}})}^T(W_{\mathrm{pp}}*R_{\mathrm{pp}0}-D_{\mathrm{pp}})\\ +2{(W_{\mathrm{pp}}*G_{\mathrm{pp}})}^T(W_{\mathrm{pp}}*G_{\mathrm{pp}})\mathrm{\ΔV}+{(W_{\mathrm{ps}}*R_{\mathrm{ps}0}-D_{\mathrm{ps}})}^T(W_{\mathrm{ps}}*G_{\mathrm{ps}})\\ +{(W_{\mathrm{ps}}*G_{\mathrm{ps}})}^T(W_{\mathrm{ps}}*R_{\mathrm{ps}0}-D_{\mathrm{ps}})+2{(W_{\mathrm{ps}}*G_{\mathrm{ps}})}^T(W_{\mathrm{ps}}*G_{\mathrm{ps}})\mathrm{\ΔV}=0\end{array},---\left(1.9\right)$

所以根据公式(1.9)，可以得出公式(1.10)，如下所示：

[(W_pp*G_pp)^T(W_pp*G_pp)+(W_ps*G_ps)^T(W_ps*G_ps)]△V，   (1.10)＝(W_pp*G_pp)^T(D_pp-W_pp*R_pp0)+(W_ps*G_ps)^T(D_ps-W_ps*R_ps0)

其中：

以上，大略说明了本发明之实施例所需要运用到的相关公式，以下将提供对应的实施例来进行说明。根据本发明的实施例，提供了一种目标函数的计算方法，其适于地震波的多波叠前联合反演。

图3是根据本发明实施例的多波联合AVA反演方法的流程图，如图3所示，该方法包括：

步骤S302，接收PP波的转换波道集与PS波的转换波道集，并提取所述PP波的不同入射角的子波以及压缩后的所述PS波的不同入射角的子波，其中压缩后的所述PS波为压缩所述PS波至PP波的预设时。其中，所述PP波的预设时例如为PP波的T0时。

步骤S304，根据所述不同入射角的子波，建立第一矩阵与第二矩阵。进一步来说，例如建立公式(1.10)中的W_pp与W_ps，其中W_pp为所述第一矩阵，W_ps为所述第二矩阵。并且，所述不同入射角的子波例如为文本格式(TXT)。

步骤S306，根据薄互层初始模型，建立第一雅各比矩阵与第二雅各比矩阵及对应所述薄互层初始模型的第一反射系数矩阵与第二反射系数矩阵。进一步来说，例如建立公式(1.10)中的G_pp、G_ps、R_pp0与R_ps0，其中G_pp为所述第一雅各比矩阵，G_ps为所述第二雅各比矩阵，R_pp0为所述第一反射系数矩阵，R_ps0为所述第二反射系数矩阵。并且，所述薄互层初始模型例如为标准地震数据记录格式(SEGY)。

步骤S308，根据所述PP波的转换波道集与所述PS波的转换波道集，建立第三矩阵与第四矩阵，并计算第一目标函数与第二目标函数。进一步来说，例如建立公式(1.10)中的(D_pp-W_pp*R_pp0)与(D_ps-W_ps*R_ps0)，其中(D_pp-W_pp*R_pp0)为所述第三矩阵，(D_ps-W_ps*R_ps0)为所述第四矩阵，而第一目标函数例如为公式(1.3)，第二目标函数例如为公式(1.4)。并且，所述PP波的转换波道集与所述PS波的转换波道集例如也为标准地震数据记录格式(SEGY)。

步骤S310，根据所述第一矩阵、所述第二矩阵、所述第三矩阵、所述第四矩阵、所述第一雅各比矩阵、所述第二雅各比矩阵、所述第一反射系数矩阵与所述第二反射系数矩阵，得到修正值，并修正所述薄互层初始模型，且调整所述第一目标函数与所述第二目标函数。其中，例如根据公式(1.10)求解，可得到△V，其中△V例如为修正值。

另外，[(W_pp*G_pp)^T(W_pp*G_pp)+(W_ps*G_ps)^T(W_ps*G_ps)]在解矩阵方程需要加入阻尼矩阵，以克服奇异，其中阻尼矩阵为λI，I为单位向量矩阵，λ取极小值，例如0.000001等。由于奇异出现的情况不可预料，有时λ取值很难决断。一般来说，若λ取值误差较大，一般不会影响反演结果，但是个别情况会导致整个反演结果发生整体漂移。这个漂移量可以用井数据来校正。

综上所述，根据本发明的技术方案，采用多波地震数据以及测井共同约束反演的结果，基于广义反演理论同时反演薄互层的纵波速度、横波速度以及密度。其中，PP波与PS波的反射系数采用策普里兹(Zoeppritz)方程组精确求解。由于反射系数的求解过程无任何近似，所以本发明的目标函数的计算方法即使在无井约束的情况下也能对薄互层的储层参数精确反演。

以上所述仅为本发明的实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (7)

1.一种多波联合振幅随入射角变化反演方法，适于地震波的多波叠前联合反演，其特征在于，所述目标函数的计算方法包括：

1)接收反射纵波PP波的转换波道集与反射转换横波PS波的转换波道集，并提取所述PP波的不同入射角的子波以及压缩后的所述PS波的不同入射角的子波，其中压缩后的所述PS波为压缩所述PS波至PP波的预设时；

2)根据所述不同入射角的子波，建立第一矩阵与第二矩阵；

根据薄互层初始模型，建立第一雅各比矩阵与第二雅各比矩阵及对应所述薄互层初始模型的第一反射系数矩阵与第二反射系数矩阵；

3)根据所述PP波的转换波道集与所述PS波的转换波道集，建立第三矩阵与第四矩阵，并计算第一目标函数与第二目标函数；

4)根据所述第一矩阵、所述第二矩阵、所述第三矩阵、所述第四矩阵、所述第一雅各比矩阵、所述第二雅各比矩阵、所述第一反射系数矩阵与所述第二反射系数矩阵，得到修正值，并修正所述薄互层初始模型，且调整所述第一目标函数与所述第二目标函数。

2.根据权利要求1所述的目标函数的计算方法，其特征在于，所述第一矩阵、所述第二矩阵、所述第三矩阵、所述第四矩阵、所述第一雅各比矩阵、所述第二雅各比矩阵、所述第一反射系数矩阵、所述第二反射系数矩阵与所述修正值组合成如下公式：

[(W_pp*G_pp)^T(W_pp*G_pp)+(W_ps*G_ps)^T(W_ps*G_ps)]△V，＝(W_pp*G_pp)^T(D_pp-W_pp*R_pp0)+(W_ps*G_ps)^T(D_ps-W_ps*R_ps0)

其中，W_pp为所述第一矩阵，W_ps为所述第二矩阵，G_pp为所述第一雅各比矩阵，G_ps为所述第二雅各比矩阵，R_pp0为所述第一反射系数矩阵，R_ps0为所述第二反射系数矩阵，(D_pp-W_pp*R_pp0)为所述第三矩阵，(D_ps-W_ps*R_ps0)为所述第四矩阵，△V为所述修正值。

3.根据权利要求2所述的目标函数的计算方法，其特征在于，所述公式更加入阻尼矩阵；

其中所述阻尼矩阵为λI，其中I为单位向量矩阵，λ取极小值。

4.根据权利要求2所述的目标函数的计算方法，其特征在于，所述第一目标函数与第二目标函数分别为如下公式：

Q(V)＝‖S_pp-D_pp‖²+‖S_ps-D_ps‖²，

$T=\begin{array}{c}{\mathrm{\λ}}_1\·\left(\right|S_{\mathrm{pp}}-D_{\mathrm{pp}}|+|S_{\mathrm{ps}}-D_{\mathrm{ps}}\left|\right)\\ +{\mathrm{\λ}}_2\·\left(\right|v_{\mathrm{pi}}-v_{\mathrm{pw}}|+|v_{\mathrm{si}}-v_{\mathrm{sw}}|+{|\mathrm{\ρ}}_i-{\mathrm{\ρ}}_w|)\end{array},$

其中，D_pp、D_ps为实际地震记录，S_pp＝W_pp*R_pp和S_ps＝W_ps*R_ps为模型响应，i为反演得到的参数，w为测井模型外推得到的参数，λ₁、λ₂分别对应反演中地震道与测井模型约束作用的权重。

5.根据权利要求1所述的目标函数的计算方法，其特征在于，所述调整所述第一目标函数与所述第二目标函数包括迭代所述步骤1)、2)、3)、4)至所述第一目标函数与所述第二目标函数的值最小为止。

6.根据权利要求1所述的目标函数的计算方法，其特征在于，所述薄互层初始模型、所述PP波的转换波道集与PS波的转换波道集为标准地震数据记录格式。

7.根据权利要求1所述的目标函数的计算方法，其特征在于，所述不同入射角的子波为文本格式。