CN104965224A - 用平均入射角道集进行pp波与ps波联合avo反演方法 - Google Patents

Abstract

本发明公开了一种用平均入射角道集进行PP波与PS波联合AVO反演方法，其中该方法包括：计算PS波压缩至PP波时间域的压缩系数；将所述压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集；将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集；将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波的地震子波；将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。通过本发明，以解决现有技术存在的反演方法出现较大的误差的问题。

Description

用平均入射角道集进行PP波与PS波联合AVO反演方法

技术领域

本发明涉及地震勘探的技术领域，尤其涉及一种用平均入射角道集进行PP波与PS波联合AVO反演方法。

背景技术

随着我国油气田勘探开发的不断深入，构造油气藏越来越少，致密油、页岩气、煤层气等非常规油气资源已经成为勘探的重要方向。在这过程中，地震勘探的技术手段也在不断的发展，以应对日益复杂的油气勘探问题。在诸多地震勘探新技术中，多分量地震技术在解决复杂储层的勘探问题上表现出明显的优势。大量的理论研究表明：利用多分量地震得到的PP波与PS波进行联合AVO反演可以提供可靠的纵、横波速度以及密度信息，并派生出更多的反映地层岩性与含流体性的属性参数，这在很大程度上降低了地震解释的多解性，为复杂油气藏的预测提供了新的思路。

多波AVO反演的基础是通过Zoeppritz方程求解地震反射系数，但该方程的形式比较复杂，没有给出反射系数与地层弹性参数之间的“显式”关系。所以国内外研究AVO反演的一个重要方面是对Zoeppritz方程进行近似，以获得反射系数与地层弹性参数之间的线性关系。这其中，最常用的近似公式是Aki-Richards近似公式，其是在地层界面两侧的弹性参数为弱反差的假设基础上建立的，可用来从反射系数直接反推单界面两侧纵、横波速度以及密度变化量。该理论在纵波上首先得到突破，并在实际勘探中得到充分应用。

虽然，PP波与PS波联合AVO反演在理论上具有明显的先进性，但是在实际应用时见效并不明显。原因主要存在于以下两个方面：(1)首先是PS波入射角道集的信噪比较差，降低了联合反演的可靠性。AVO反演需要获得高质量叠前道集作为保障，但现在处理水平尚难以保证PS波的AVA道集具有较宽的入射角范围，也不能保证每个入射角的覆盖次数均匀、信噪 比较高。(2)另一个原因是反演中用到的反射系数近似公式都是建立在弱反差的假设前提下，无法适用于强反差的地层界面。很多储集层，如含煤地层、非固结砂岩地层、火成岩与页岩地层等，其弹性参数与围岩的差异较大，导致基于Aki-Richards反射系数近似公式的AVO反演方法出现较大的误差。

发明内容

本发明的主要目的在于提供一种用平均入射角道集进行PP波与PS波联合AVO反演方法，以解决现有技术存在的反演方法出现较大的误差的问题。

为解决上述问题，本发明实施例提供一种用平均入射角道集进行PP波与PS波联合AVO反演方法，包括：计算PS波压缩至PP波时间域的压缩系数；将所述压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集；将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集；将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波的地震子波；将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。

其中，所述计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集的步骤之后还包括：基于所述合成地震记录与实际数据的平均入射角道集计算目标函数；判断所述目标函数是否小于预设值，并产生判断结果；当所述判断结果为所述目标函数小于预设值时，输出模型修改量向量；当所述判断结果为所述目标函数未小于预设值时产生初始模型修正量，以修正所述初始模型，并回到计算PS波压缩至PP波时间域的压缩系数的步骤。

其中，所述目标函数满足如下公式：Q(V)＝||W_PP*R_PP-D_PP||²+||W_PS*R_PS-D_PS||²,其中， 分别为实际PP波与PS波在平均入射角分别为θ₁和θ₂的地震记录向量；W_PP、W_PS分别为PP波与PS波的地震子波向量； 分别为PP波与PS波在平均入射角分别为θ₁和θ₂的反射系数。

其中，所述所述模型修改量向量满足如下公式：  $\begin{array}{c}\mathrm{\ΔV}={[{(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}})]}^{-1}\\ \·[{(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},$ 其中，G_PP、G_PS分别为PP波与PS波的雅各比矩阵；R_PP0、R_PS0分别表示R_PP、R_PS在初始模型处进行泰勒展开的零阶项。

其中，所述PP波与PS波的反射系数满足如下公式：

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

其中，R_PP为纵波反射系数；R_PS为转换横波反射系数；T_PP为纵波透射系数；T_PS为转换横波透射；p为射线参数；V_P1和V_S1分别为上层介质纵、横波速度；V_P2和V_S2分别为下层介质的纵、横波速度；ρ₁和ρ₂分别为上、下层介质的密度；α₁和α₂为纵波入射角或反射角和透射角；β₁和β₂为横波反射角和透射角。

根据本发明的技术方案，通过使用的输入数据是平均入射角道集，且利用一个统一的目标函数，将PP波与PS波的正演道集与实际数据平均入射角道集之间的残差结合在一起，以计算精确的PP波与PS波反射系数，进行纵波速度、横波速度以及密度的同时最小二乘反演，使得本发明所获得的参数较为精确。

附图说明

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

图1a与图1b分别为通过Zoeppritz方程计算的PP波反射系数及PS波 反射系数的AVO曲线；

图2是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演装置的方块图；

图3是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演装置的另一方块图；

图4是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演方法的流程图；

图5是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演方法的另一流程图；

图6a、图6c分别是PP波与PS波的AVA道集的波形图；

图6b、图6d分别由5度与20度入射角组成的PP波与PS波道集的波形图；

图7a、图7b、图7c分别为PP波与PS波的5度与20度角道集用本发明的反演装置及方法反演得到的结果的示意图；

图8a、图8b、图8c分别是基于近似反射系数的反演结果；

图9a、图9c分别是加入20％的随机噪音后的PP波与PS波道集的波形图；

图9b、图9d分别由5度与20度入射角组成且加入20％的随机噪音后的的PP波与PS波道集的波形图；

图10a、图10b、图10c分别为在加入20％的随机噪音下PP波与PS波的5度与20度角道集用本发明的反演装置及方法反演得到的结果的示意图。

具体实施方式

本发明的主要思想在于，基于所使用的输入数据是平均入射角道集，且利用一个统一的目标函数，将PP波与PS波的正演道集与实际数据平均入射角道集之间的残差结合在一起，以计算精确的PP波与PS波反射系数，进行纵波速度、横波速度以及密度的同时最小二乘反演，且进一步在反演过程中，通过多次迭代修改初始模型，使得目标函数达到最小，进而使得本发 明所获得的参数较为精确。

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

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

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

其中，p为射线参数；V_P1和V_S1分别为上层介质纵、横波速度；V_P2和V_S2分别为下层介质的纵、横波速度；ρ₁和ρ₂分别为上、下层介质的密度；α₁和α₂为纵波入射角(反射角)和透射角；β₁和β₂为横波反射角和透射角。

令纵波反射系数为R_PP、转换横波反射系数为R_PS、纵波透射系数为T_PP和转换横波透射为T_PS，可以得出四个波的位移振幅应当满足的策普里兹(Zoeppritz)方程组，如公式(1.2)所示：

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

基于公式(1.2)，给出界面两侧介质的弹性参数，可以通过线性方程组的数值解法计算PP波与PS波精确的反射系数。

为了验证通过PP波与PS波AVO曲线上的一个小入射角和一个大入射角的反射系数组合在一起是否可以反推地层的岩性与含油气性，建立了三种等效地层模型，模型参数如表1所示，通过公式(1.2)计算精确的AVO曲线显示在图1a与图1b中。取5度小入射角与20度大入射角，在图1a与图1b中分别拾取了PP波与PS波反射系数曲线上的对应点，如图1a与图1b中虚线所示，并将“两点”进行了直线连接，可以看出：单独用PP波的“两点”反射系数连线(图1a中的虚线)无法区分地层的组合模式，尤其是页岩/水饱和砂岩界面与页岩/气饱和砂岩界面的“两点”反射系数连线非常接近，无 法区分；只有灰岩/气饱和砂岩界面的PP波“两点”反射系数连线与其它两条连线的差异较大。所以单独利用PP波两个角度的平均入射角道集无法保证AVO反演的可靠性。

将图1b与图1a进行对比，可以看出：单独PS波的“两点”反射系数连线(图1b中的虚线)难以区分页岩/气饱和砂岩界面与灰岩/气饱和砂岩界面，但页岩/水饱和界面的PS波“两点”反射系数连线表现出与其它两条连线差异较大。所以单独利用PS波两个角度的平均入射角道集也无法反演地层含油气性的真实情况。但是，如果将PP波与PS波“两点”反射系数连线组合在一起，相当于形成四个角度的入射角道集，则可以把三种界面类型完全区分开。该结论也为PP波与PS波的平均入射角道集的联合反演提供了理论依据。

平均入射角道集相当于把PP波与PS波的AVA道集的角度各自限制成两个角度。其中，较小的平均入射角地震道是由AVA道集中的小入射角范围的地震数据通过加权叠加来获得；同样，较大的平均入射角地震道是通过相对较大的入射角范围的地震数据进行加权叠加来获得。平均入射角的多波局部叠加数据可以用来反演弹性参数。为了使得平均入射角地震道的振幅能够逼近该角度的真实地震振幅，叠前的数据处理要尽量采用保幅算法。平均入射角的选择要根据多分量地震数据采集的实际情况来选择，但为了保持PP波与PS波反射系数的变化趋势，选择的大、小平均入射角要尽量差异大些，并且在实际数据处理中具有可操作性。

在纵波速度与密度之间关系满足Gardner公式以及界面两侧地层的纵、横波速度变化率接近的情况下，PP波单独反演的最佳平均入射角为0度；PS波应选择中远偏移距入射角，考虑到地层的吸收衰减与动校正拉伸率的因素，30度平均入射角比较适合PS波的单独反演。但是在PP波与PS波的联合反演中，PP波叠加的优势振幅集中在中小偏移距，而PS波叠加的优势振幅集中在中大偏移距；所以0度与30度的平均入射角地震道无法同时适用于PP波与PS波的反演。由于PS波在0度入射角时的反射系数为0，联合反演应该选择非零度的小入射角，因此入射角较佳为为5度。对于大入射角，折衷考虑到PP波在大入射角反射系数有可能较小，因此联合反演的大 入射角较佳选择为20度。

表1三种含流体地层模型的弹性参数

另外，为了使得本发明的反演方法能够适用于更为广泛的地层条件，将基于Zoeppritz方程计算精确的反射系数来进行PP波与PS波的联合反演，以提高反演的精度。通过最小二乘法建立的联合反演目标函数Q，如公式(1.3)所示：

Q(V)＝||W_PP*R_PP-D_PP||²+||W_PS*R_PS-D_PS||²    (1.3) 

其中，分别为实际PP波与PS波在平均入射角分别为θ₁和θ₂的地震记录向量；W_PP、W_PS分别为PP波与PS波的地震子波向量，分别为PP波与PS波在平均入射角分别为θ₁和θ₂的反射系数。V＝(V_P,V_S,ρ)为界面两侧的地层弹性参数组成的向量，并且有V_P＝(V_P1,V_P2)、V_S＝(V_S1,V_S2)、ρ＝(ρ₁,ρ₂)。

对PP波和PS波的反射系数向量在初始模型附近进行泰勒展开，统一的矩阵表述形式如公式(1.4)所示：

$R(V_0+\mathrm{\ΔV})=R\left(V_0\right)+\frac{\∂R}{\∂V}\·\mathrm{\ΔV}---\left(1.4\right)$

其中，V₀为纵波速度、横波速度和密度的初始猜测值组成的向量，ΔV为初始模型修正量。令雅各比矩阵将公式(1.4)代入公式(1.3)，可 得公式(1.5)，如下所示：

Q(V)＝||W_PP*(R_PP0+G_PPΔV)-D_PP||²+||W_PS*(R_PS0+G_PSΔV)-D_PS||²    (1.5) 

将公式(1.5)写成矩阵的展开形式，如公式(1.6)所示：

Q(V)＝(W_PP*R_PP0-D_PP)^T(W_PP*R_PP0-D_PP)+(W_PP*R_PP0-D_PP)^T(W_PP*G_PPΔV)

+(W_PP*G_PPΔV)^T(W_PP*R_PP0-D_PP)+(W_PP*G_PPΔV)^T(W_PP*G_PPΔV)

                                                        (1.6) 

+(W_PS*R_PS0-D_PS)^T(W_PS*R_PS0-D_PS)+(W_PS*R_PS0-D_PS)^T(W_PS*G_PSΔV)

+(W_PS*G_PSΔV)^T(W_PS*R_PS0-D_PS)+(W_PS*G_PSΔV)^T(W_PS*G_PSΔV)

为了使得目标函数Q(V)取值最小，要求存在所以公式(1.6)可转换成公式(1.7)，如下所示： 

$\begin{array}{c}\frac{\∂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}}\mathrm{\ΔV})}^T(W_{\mathrm{PP}}*G_{\mathrm{PP}})\\ +[{(W_{\mathrm{PS}}*R_{\mathrm{PS}0}-D_{\mathrm{PS}})}^T(W_{\mathrm{PS}}*G_{\mathrm{PS}})+{(W_{\mathrm{PS}}*G_{\mathrm{PS}}\mathrm{\ΔV})}^T(W_{\mathrm{PS}}*G_{\mathrm{PS}})]=0\end{array}---\left(1.7\right)$

进一步可求解ΔV，如公式(1.8)所示：

$\begin{array}{c}\mathrm{\ΔV}={[{(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}})]}^{-1}\\ \·[{(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}---\left(1.8\right)$

其中，G_PP、G_PS分别为PP波与PS波的雅各比矩阵；R_PP0、R_PS0分别表示R_PP、R_PS在初始模型处进行泰勒展开的零阶项。将上式中的ΔV与初始的弹性参数向量V₀进行向量求和，可得到更新后的初始模型，并进行上述过程的迭代；当目标函数Q(V)取值最小时，反演过程结束。更新后的弹性参数向量V₀即为模型的反演结果。

在实际的联合反演中，输入的是由一系列采样点组成的时间序列，公式(1.8)中的向量需要进行扩展，以适用于时间序列模型的反演。扩展后的雅各比矩阵G_PP与G_PS具有统一形式，如公式(1.9)所示：

$G=\left[\begin{array}{cccc}\frac{{\∂R}_1}{{\∂V}_1}& \frac{{\∂R}_1}{{\∂V}_2}& L& \frac{{\∂R}_1}{{\∂V}_n}\\ \frac{{\∂R}_2}{{\∂V}_1}& \frac{{\∂R}_2}{{\∂V}_2}& L& \frac{{\∂R}_2}{{\∂V}_n}\\ M& M& O& M\\ \frac{{\∂R}_n}{{\∂V}_1}& \frac{{\∂R}_n}{{\∂V}_2}& L& \frac{{\∂R}_n}{{\∂V}_n}\end{array}\right]---\left(1.9\right)$

公式(1.9)中的下标1，2…n表示不同的时间采样点，R_i表示第i个采样点位置的PP波或者PS波反射系数，表示第i个采样点的PP波或者PS波的反射系数随第j个采样点的弹性参数变化量，并且有公式(1.10)，如下所示：

$\frac{{\∂R}_i}{{\∂V}_j}=(\frac{{\∂R}_i}{{\∂V}_{\mathrm{Pj}}},\frac{{\∂R}_i}{{\∂V}_{\mathrm{Sj}}},\frac{{\∂R}_i}{{\∂\mathrm{\ρ}}_j})---\left(1.10\right)$

在上式中，反射系数的求解必须由公式(1.2)进行线性方程组的数值求解；若用Aki-Richards反射系数近似公式建立公式(1.9)的雅各比矩阵，则本发明的反演方法退化为基于Aki-Richards反射系数近似公式的反演方法。

为了在一次迭代中同时修改时间采样点的弹性参数，PS波的平均入射角道集必须和PP波在同一时间尺度上，所以必须把PS波道集压缩至PP波的时间上。该过程可以通过PP波与PS波的叠加剖面进行层位的对比解释，以获得不同层位的压缩系数，从而实现PS波的压缩。若同一地质界面的PP波与PS波反射波同相轴面貌差异太大，则需要工区内有足够多的测井合成记录来辅助完成层位的匹配。另外，地震子波需要在PS波平均入射角道集压缩至PP波时间后再提取，以保证地震子波时间尺度的一致性。由于PP与PS波的子波类型存在差异，所以PP波与PS波必须分别提取子波用于反演。

以上，大略说明了本发明之实施例所需要运用到的相关公式，以下将提供对应的实施例来进行说明。根据本发明的实施例，提供了一种用平均入射角道集进行PP波与PS波联合AVO反演装置。

图2是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演装置的方块图。所述用平均入射角道集进行PP波与PS波联合AVO反演装置200包括第一计算模块210、压缩模块220、选取模块230、叠加及提取模块240、第二计算模块250。

所述第一计算模块210计算PS波压缩至PP波时间域的压缩系数。

所述压缩模块220连接所述第一计算模块210，以将压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集。

所述选取模块230连接所述压缩膜块220，以将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集。

所述叠加及提取模块240连接所述选取模块130，以将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波的地震子波。

所述第二计算模250块连接所述叠加及提取模块240，以将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。

图3是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演装置的另一方块图。

所述用平均入射角道集进行PP波与PS波联合AVO反演装置300包括第一计算模块210、压缩模块220、选取模块230、叠加及提取模块240、第二计算模块250、第三计算模块310、判断模块320与控制模块330。其中，第一计算模块210、压缩模块220、选取模块230、叠加及提取模块240、第二计算模块250的连接关系及操作，可参考图2的实施例，故在此不再赘述。

所述第三计算模块310连接所述所述第二计算模块250，以基于所述合成地震记录与实际数据的平均入射角道集计算目标函数。

所述判断模块320判断所述目标函数是否小于预设值，以产生判断结果。

所述控制模块330连接所述判断模块320、所述第一计算模块210、压缩模块220、选取模块230、叠加及提取模块240、第二计算模块250，以当所述判断结果为所述目标函数小于预设值时输出模型修改量向量。当所述判断结果为所述目标函数未小于预设值时，所述控制模块330产生初始模型修正量，以修正所述初始模型，并控制所述第一计算模块210、压缩模块220、选取模块230、叠加及提取模块240、第二计算模块250重新启动，通过多次迭代修改初始模型，使得目标函数达到最小为止。

另外，根据本发明的实施例，提供了一种用平均入射角道集进行PP波与PS波联合AVO反演方法。

图4是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演方法的流程图。

步骤S402，计算PS波压缩至PP波时间域的压缩系数。

步骤S404，将所述压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集，使得PS波的平均入射角道集和PP波在同一时间尺度上。

步骤S406，将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集。其中，第一入射角范围与第二入射角范围不同，亦即第一入射角范围小于第二入射角范围，或第一入射角范围大于第二入射角范围

步骤S408，将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波的地震子波。

步骤S410，将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。其中，PP波与PS波的反射系数例如根据公式(1.2)进行求解。

图5是根据本发明实施例的用平均入射角道集进行PP波与PS波联合AVO反演方法的另一流程图。

步骤S402，计算PS波压缩至PP波时间域的压缩系数。

步骤S404，将所述压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集，使得PS波的平均入射角道集和PP波在同一时间尺度上。

步骤S406，将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集。其中，第一入射角范围与第二入射角范围不同，亦即第一入射角范围小于第二入射角范围，或第一入射角范围大于第二入射角范围。

步骤S408，将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波 的地震子波。

步骤S410，将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。其中，PP波与PS波的反射系数例如根据公式(1.2)进行求解。

步骤S502，基于所述合成地震记录与实际数据的平均入射角道集计算目标函数。其中，所述目标函数可利用公式(1.3)进行求解。

步骤S504，判断所述目标函数是否小于预设值，并产生判断结果。

步骤S506，当所述判断结果为所述目标函数小于预设值时，输出模型修改量向量。

步骤S508，当所述判断结果为所述目标函数未小于预设值时产生初始模型修正量，以修正所述初始模型，并回到步骤S410，重复上述过程，在目标函数达到最小时停止。其中，所述初始模型修正量为公式(1.8)，可通过公式(1.3)～(1.7)进行求解。

为了验证本发明的反演装置及方法，提供了一种含有薄层的地层模型，其中模型参数如表2所示，第3层的厚度较薄，且与围岩的弹性参数差异较大，其顶、底界面均为强反差界面；其它的地层界面均为弱反差界面。合成记录是用射线追踪基于精确的反射系数模拟得到。PP波与PS波的地震子波分别设置为40Hz以及30Hz的雷克子波，采样率为1ms；PS波的反射时间已经压缩至PP波时间上。

图6a、6c分别是PP波与PS波的AVA道集的波形图，图6b、6d分别由5度与20度入射角组成的PP波与PS波道集的波形图。图7a、图7b、图7c分别为PP波与PS波的5度与20度角道集用本发明的反演装置及方法反演得到的结果的示意图。其中，标号701为初始模型曲线、标号702为真实模型曲线、标号703为平均入射角道集联合反演曲线。

图7a、图7b、图7c可以看出：在初始模型没有任何界面变化的情况下，基于精确反射系数的反演装置及方法能够获得与真实模型非常接近的反演结果；尤其是第3层的薄层，反演曲线在界面位置的突变与真实模型匹配较好。总体上看，与真实模型相比，本发明的反演装置及方法得到模型曲线在地层弹性参数突变的位置略有一些平滑，但是误差不大。

图8a、图8b、图8c是基于近似反射系数的反演结果。其中，标号801为初始模型曲线，标号802为真实模型曲线，标号802为平均入射角道集联 合反演曲线。图8a、图8b、图8c与图7a、图7b、图7c对比可见：基于近似反射系数的反演没能很好的反映地层弹性参数的真实变化规律；即使在弱反差的界面上，其结果误差也较大，尤其是第2层的纵波速度反演结果明显偏小、密度明显偏大。上述数值分析的结果表明：对于一些弹性参数反差较大、厚度较薄的地层，必须采用精确的反射系数进行AVO反演；即使对于弱反差的地层，基于Zoeppritz方程精确反射系数的AVO反演的精度也较高。

在地震记录上加入20％的随机噪音后，测试噪音对本发明的反演装置及方法的影响。如图9b、9d所示，5度与20度平均入射角地震道在加噪的AVA道集(图9a、9c)的基础上用入射角为0-10度的地震道以及入射角为15-25度的地震道分别加权叠加得到。用本发明的反演装置及方法分别对AVA道集以及平均入射角道集进行基于Zoeppritz方程精确反射系数的反演，如图10a、图10b、图10c所示，平均入射角道集反演结1001与AVA道集反演结果1003都很靠近真实模型1002；但是，AVA道集反演结果1003的噪音振幅强于平均入射角道集反演结1001。所以，AVA道集的直接反演相比，用平均入射角道集进行反演具有更强的抗噪性。

综上所述，根据本发明的技术方案，采用的输入数据是平均入射角道集，且利用一个统一的目标函数，将PP波与PS波的正演道集与实际数据平均入射角道集之间的残差结合在一起，以计算精确的PP波与PS波反射系数，进行纵波速度、横波速度以及密度的同时最小二乘反演，且进一步在反演过程中，通过多次迭代修改初始模型，使得目标函数达到最小，进而使得本发明所获得的参数较为精确。

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

Claims (5)

1.一种用平均入射角道集进行PP波与PS波联合AVO反演方法，其特征在于，包括：

计算PS波压缩至PP波时间域的压缩系数；

将所述压缩系数应用至PS波AVA道集上得到PP波时间域的PS波AVA道集；

将PP波与压缩后的PS波AVA道集分选成第一入射角范围和第二入射角范围的AVA道集；

将PP波与PS波的第一、第二入射角范围的AVA道集各自加权叠加，形成平均入射角道集，并在同一时窗内分别提取PP波与PS波的地震子波；

将测井得到的P波速度、S波速度与密度外推成初始模型，以计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集。

2.根据权利要求1所述的用平均入射角道集进行PP波与PS波联合AVO反演方法，其特征在于，所述计算PP波与PS波的反射系数，并与所述地震子波合成平均入射角道集的步骤之后还包括：

基于所述合成地震记录与实际数据的平均入射角道集计算目标函数；

判断所述目标函数是否小于预设值，并产生判断结果；

当所述判断结果为所述目标函数小于预设值时，输出模型修改量向量；

当所述判断结果为所述目标函数未小于预设值时产生初始模型修正量，以修正所述初始模型，并回到计算PS波压缩至PP波时间域的压缩系数的步骤。

3.根据权利要求2所述的用平均入射角道集进行PP波与PS波联合AVO反演方法，其特征在于，所述目标函数满足如下公式：

Q(V)＝||W_PP*R_PP-D_PP||²+||W_PS*R_PS-D_PS||²，

其中， $D_{\mathrm{PP}}=(D_{\mathrm{PP}}^{\left({\mathrm{\θ}}_1\right)},D_{\mathrm{PP}}^{\left({\mathrm{\θ}}_2\right)}),D_{\mathrm{PS}}=(D_{\mathrm{PS}}^{\left({\mathrm{\θ}}_1\right)},D_{\mathrm{PS}}^{\left({\mathrm{\θ}}_2\right)})$ 分别为实际PP波与PS波在平均入射角分别为θ₁和θ₂的地震记录向量；W_PP、W_PS分别为PP波与PS波的地震子波向量； $R_{\mathrm{PP}}=(R_{\mathrm{PP}}^{\left({\mathrm{\θ}}_1\right)},R_{\mathrm{PP}}^{\left({\mathrm{\θ}}_2\right)}),R_{\mathrm{PS}}=(R_{\mathrm{PS}}^{\left({\mathrm{\θ}}_1\right)},R_{\mathrm{PS}}^{\left({\mathrm{\θ}}_2\right)})$ 分别为PP波与PS波在平均入射角分别为θ₁和θ₂的反射系数。

4.根据权利要求3所述的用平均入射角道集进行PP波与PS波联合AVO反演方法，其特征在于，所述模型修改量向量满足如下公式：

ΔV＝[(W_PP*G_PP)^T(W_PP*G_PP)+(W_PS*G_PS)^T(W_PS*G_PS)]^-1

·[(W_PP*G_PP)^T(D_PP-W_PP*R_PP0)+(W_PS*G_PS)^T(D_PS-W_PS*R_PS0)]，

其中，G_PP、G_PS分别为PP波与PS波的雅各比矩阵；R_PP0、R_PS0分别表示R_PP、R_PS在初始模型处进行泰勒展开的零阶项。

5.根据权利要求1所述的用平均入射角道集进行PP波与PS波联合AVO反演方法，其特征在于，所述PP波与PS波的反射系数满足如下公式：

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

其中，R_PP为纵波反射系数；R_PS为转换横波反射系数；T_PP为纵波透射系数；T_PS为转换横波透射；p为射线参数；V_P1和V_S1分别为上层介质纵、横波速度；V_P2和V_S2分别为下层介质的纵、横波速度；ρ₁和ρ₂分别为上、下层介质的密度；α₁和α₂为纵波入射角或反射角和透射角；β₁和β₂为横波反射角和透射角。