CN102928870B - 基于正则化的非线性地震叠前弹性参数反演方法 - Google Patents

Abstract

本发明提供一种基于正则化的非线性地震叠前弹性参数反演方法，包括以下步骤：根据已知的测井资料与射线弹性阻抗得出多角度的射线弹性阻抗井曲线REI_i；引入已知的层位、断层资料以及上述多角度的射线弹性阻抗井曲线REI_i，以建立整个工区不同角度的射线弹性阻抗模型；使用叠前偏移距道集，计算出上述不同入射角的部分叠加数据体；通过稀疏脉冲反演算法进行反演，得出上述入射角对应的射线弹性阻抗数据体；将上述入射角对应的射线弹性阻抗数据体，导入非线性Landweber-Fridman迭代算法，以得出整个工区的纵、横波速度和岩石密度，进而计算出其它岩性参数。本发明利用多个角度的REI_i数据，并通过非线性Landweber-Fridman迭代法进行反演计算，以获得准确的岩性参数。

Description

基于正则化的非线性地震叠前弹性参数反演方法

技术领域

本发明涉及石油天然气地震勘探领域，尤其是一种基于正则化的非线性地震叠前弹性参数反演方法。

背景技术

传统的叠后波阻抗反演可获得波阻抗剖，运用于储层预测中效果较好。但近年来，随着勘探精度和难度的不断提高，使得叠后波阻抗反演解决地质问题的局限性越来越明显。基于此，叠前地震反演技术受到越来越多的重视，成为地球物理勘探领域一项正在兴起的新技术。

叠前地震反演技术充分利用叠前信息，可以得到除波阻抗之外的很多其他弹性参数信息，大大丰富了储层预测的手段，增强了对复杂储层的描述和流体检测的能力，因此从地震资料中提取介质的弹性参数，并将这些参数与岩性和流体成分联系起来，在油藏描述和油气监测中发挥着重要作用。

现有技术1：

现有的叠前地震反演技术以Connolly提出的弹性阻抗反演(EI)为主。根据Connolly的定义，弹性阻抗是纵横波速度，密度以及入射角的函数。其具体步骤如下：

1、从叠前道集中提取特定入角资料；

2、利用纵、横波速度和密度测井资料计算对应入角的弹性阻抗；

3、以角道集资料代替叠后反演中零炮检距资料，以弹性阻抗井曲线代替传统的波阻抗井曲线，利用测井约束反演软件实现弹性阻抗反演

Connolly弹性阻抗公式：

Connolly根据Richards近似公式和Shuey线性近似公式给出了弹性阻抗(EI)的计算公式：

$EI\left(\mathrm{\θ}\right)={\mathrm{\α}}^{1+{\tan }^2\mathrm{\θ}}{\mathrm{\β}}^{-8{\mathrm{Ksin}}^2\mathrm{\θ}}{\mathrm{\ρ}}^{1-4{\mathrm{Ksin}}^2\mathrm{\θ}}.---\left(1\right)$

该公式的不足在于，随着入射角的变化EI的量纲是变化的，其数值范围也随入射角变化，使得El和AI不容易对比分析。为克服这一弱点，通常选一组比例因子ρ₀对v_p,v_s,ρ进行归一化，这样便得到归一化后的弹性阻抗公式：

$EI\left(\mathrm{\θ}\right)={(v_p/v_{p_0})}^{1+{\tan }^2\mathrm{\θ}}{(v_s/v_{s_0})}^{-8{\mathrm{Ksin}}^2\mathrm{\θ}}{(\mathrm{\ρ}/{\mathrm{\ρ}}_0)}^{1-4{\mathrm{Ksin}}^2\mathrm{\θ}}.---\left(2\right)$

现有技术2：

在弹性阻抗EI的基础上，将射线弹性阻抗与P波反射波之间的关系定义为：

$R_{{\mathrm{pp}}_i}=\frac{{\mathrm{EI}}_{i+1}({\mathrm{\α}}_{i+1},{\mathrm{\β}}_{i+1},{\mathrm{\ρ}}_{i+1},\mathrm{\θ})-{\mathrm{EI}}_i({\mathrm{\α}}_i,{\mathrm{\β}}_i,{\mathrm{\ρ}}_i,\mathrm{\θ})}{{\mathrm{EI}}_{i+1}({\mathrm{\α}}_{i+1},{\mathrm{\β}}_{i+1},{\mathrm{\ρ}}_{i+1},\mathrm{\θ})+{\mathrm{EI}}_i({\mathrm{\α}}_i,{\mathrm{\β}}_i,{\mathrm{\ρ}}_i,\mathrm{\θ})}---\left(3\right)$

在上下岩层弹性参数差异不大的前提假设下，上式可以写成差分的形式，即：

$R_{{\mathrm{pp}}_i}\left(p\right)\≈\frac{1}{2}\mathrm{\Δ}\ln EI(\mathrm{\α},\mathrm{\β},\mathrm{\ρ},\mathrm{\θ})---\left(4\right)$

沿着射线路径R积分，就获得了射线弹性阻抗REI

$REI(\mathrm{\α},\mathrm{\β},\mathrm{\ρ})=\frac{\mathrm{\ρ}\mathrm{\α}}{\cos \mathrm{\θ}}{(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})}^4.---\left(5\right)$

现有技术3：

在射线弹性阻抗REI的基础上，对公式取对数

$ln\left(REI\right(\mathrm{\θ})cos\mathrm{\θ})=ln\mathrm{\ρ}\mathrm{\α}+4ln(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})---\left(6\right)$

如果θ和REI(θ)看作已知，令x＝ρα，则问题转化为求方程

f(x,y)＝ln(REI(θ)cosθ)-lnx-4ln(1-y²sin²θ)＝0(7)

按θ不同写成方程组

$f(x,y)=\left[\begin{array}{c}\ln \left(REI\right({\mathrm{\θ}}_1){\mathrm{cos\θ}}_1)-\ln x-4\ln (1-y^2{\sin }^2{\mathrm{\θ}}_1)\\ \ln \left(REI\right({\mathrm{\θ}}_2){\mathrm{cos\θ}}_2)-\ln x-4\ln (1-y^2{\sin }^2{\mathrm{\θ}}_2)\end{array}\right]---\left(8\right)$

方程组中的第二个ln项都做泰勒展开

$ln(1-y^2{\sin }^2\mathrm{\θ})=-{\sin }^2{\mathrm{\θy}}^2-\frac{1}{2}{(-{\sin }^2{\mathrm{\θy}}^2)}^2+\frac{1}{3}{(-{\sin }^2{\mathrm{\θy}}^2)}^3+\mathrm{...}---\left(9\right)$

0＜sinθ＜1,0＜y＜1(10)

忽略很小的高次项得到 $ln(1-y^2{\sin }^2\mathrm{\θ})=-{\sin }^2{\mathrm{\θy}}^2+\frac{1}{2}{\sin }^4{\mathrm{\θy}}^4---\left(11\right)$

重新得到方程组为

$f(x,y)=\left[\begin{array}{c}\ln \left(REI\right({\mathrm{\θ}}_1){\mathrm{cos\θ}}_1)-\ln x-4{\sin }^2{\mathrm{\θ}}_1{y}^2-2{\sin }^2{\mathrm{\θ}}_1{y}^4\\ \ln \left(REI\right({\mathrm{\θ}}_2){\mathrm{cos\θ}}_2)-\ln x-4{\sin }^2{\mathrm{\θ}}_2{y}^2-2{\sin }^4{\mathrm{\θ}}_2{y}^4\end{array}\right]---\left(12\right)$

f的偏导数矩阵可以写作

$f^{\′}(x,y)=\left[\begin{array}{cc}\frac{\∂f_1}{\∂x}& \frac{\∂f_1}{\∂y}\\ \frac{\∂f_2}{\∂x}& \frac{\∂f_2}{\∂y}\end{array}\right]=\left[\begin{array}{cc}-\frac{1}{x}& 4{\mathrm{ysin}}^2{\mathrm{\θ}}_1-2y^3{\sin }^4{\mathrm{\θ}}_1\\ -\frac{1}{x}& 4{\mathrm{ysin}}^2{\mathrm{\θ}}_2-2y^3{\sin }^4{\mathrm{\θ}}_2\end{array}\right]---\left(13\right)$

用牛顿迭代公式x^(k+1)＝x^(k)-[f'^(k)]^-1f(x^(k))(14)

写成矩阵形式为

$\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\end{array}\right]={\left[\begin{array}{cc}-\frac{1}{x}& 4{\mathrm{ysin}}^2{\mathrm{\θ}}_1-2y^3{\sin }^4{\mathrm{\θ}}_1\\ -\frac{1}{x}& 4{\mathrm{ysin}}^2{\mathrm{\θ}}_2-2y^3{\sin }^4{\mathrm{\θ}}_2\end{array}\right]}^{-1}\left[\begin{array}{c}\ln x-4\ln (1-y^2{\sin }^2{\mathrm{\θ}}_1)-\ln \left(REI\right({\mathrm{\θ}}_1){\mathrm{cos\θ}}_1)\\ \ln x-4\ln (1-y^2{\sin }^2{\mathrm{\θ}}_2)-\ln \left(REI\right({\mathrm{\θ}}_2){\mathrm{cos\θ}}_2)\end{array}\right]---\left(15\right)$

给出初值x⁽⁰⁾，y⁽⁰⁾当||Δx,Δy||²＜δ时停止迭代。

该算法的问题在于，在做泰勒展开后，只有在θ接近零的情况下可以将高次项忽略。而实际情况是，入射角θ可以取到较大的值，使得高次项变得不可忽略。因此该算法的线性化并不完善，使用它对REI数据进行反演，难以获得准确的岩性参数。

上述的三种算法的缺点如下：

1、常规叠后反演不能得到可靠的波阻抗和其他岩性信息；

2、弹性阻抗EI以及归一化得EI都依赖于S波与P波速度比从浅到深保持某一常数不变，这和实际经常不相吻合；

3、弹性阻抗EI以及归一化的EI一般需要三个偏移距数据体才能得到纵横波速度及密度，抗噪能力比较差。

发明内容

针对上述技术的不足之处，本发明提供一种利用多个角度的REI_i数据，并通过基于正则化思想的非线性Landweber-Fridman迭代法进行反演计算，以获得准确的岩性参数的基于正则化的非线性地震叠前弹性参数反演方法。

为实现上述目的，本发明提供一种基于正则化的非线性地震叠前弹性参数反演方法，包括以下步骤：

根据已知的测井资料与射线弹性阻抗REI利用下式得出多角度的射线弹性阻抗井曲线REI_i，其中，所述测井资料中包括纵波速度、横波速度、密度与入射角，

${\mathrm{REI}}_i=REI(\mathrm{\α},\mathrm{\β},\mathrm{\ρ},{\mathrm{\θ}}_i)=\frac{\mathrm{\ρ}\mathrm{\α}}{{\mathrm{cos\θ}}_i}{(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_i)}^4,---\left(16\right)$

在上式中，α为纵波速度，β为横波速度，ρ为岩石密度，θ_i为入射角。

引入已知的层位、断层资料以及上述多角度的射线弹性阻抗井曲线REI_i，以建立整个工区不同角度的射线弹性阻抗模型；

使用叠前偏移距道集，计算出上述不同入射角的部分叠加数据体；

将多角度的部分叠加数据体、连同多角度的射线弹性阻抗模型，通过稀疏脉冲反演算法进行反演，得出上述入射角对应的射线弹性阻抗数据体；

将上述入射角对应的射线弹性阻抗数据体，导入非线性Landweber-Fridman迭代算法，以得出整个工区的纵、横波速度和岩石密度；

由于公式(16)是一个非线性问题，而且针对其进行的线性化改造是困难和不成功的，只有采用非线性反演方法才能较好地解决。

令Z_i＝REI_i，由于REI随θ值的变化而变化，做偏差平方和函数：

$J(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})=\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}{\[Z_i(\mathrm{\ρ},\mathrm{\α},\mathrm{\β};{\mathrm{\θ}}_i)-Z_i^{*}\]}^2---\left(17\right)$

令

$f_i(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})=Z_i(\mathrm{\ρ},\mathrm{\α},\mathrm{\β};{\mathrm{\θ}}_i)-Z_i^{*},(i=1,2,\mathrm{...},N)---\left(18\right)$

则有

$J(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})=\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}{f}_i^2(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})---\left(19\right)$

问题成为求向量(ρ^*,α^*,β^*)，使得

$J({\mathrm{\ρ}}^{*},{\mathrm{\α}}^{*},{\mathrm{\β}}^{*})=\underset{(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})\∈R^3}{min}J(\mathrm{\ρ},\mathrm{\α},\mathrm{\β})---\left(20\right)$

引入Landweber-Fridman迭代法：

设

F(x):X→Y(21)

为非线性映射，X,Y均为Hilbert空间。这时J[x]的梯度为

grad(J[x])＝F′(x)^T(F(x)-y)(22)

其中F′(x)表示泛函F(x)在x点的Frechet导数。沿负梯度方向进行步长为ω(ω||F′(x)||²≤1)的线搜索，则得到标准的非线性Landweber-Fridman迭代公式

$x_{k+1}^{\mathrm{\δ}}=x_k^{\mathrm{\δ}}+{\mathrm{\ωF}}^{\′}{\left(x_k^{\mathrm{\δ}}\right)}^T(y-F(x_k^{\mathrm{\δ}}))---\left(23\right)$

当右端项y不能够精确获得时，如仅仅知道y^δ，同样可以利用带扰动的非线性Landweber-Fridman迭代公式获得准确的x＝{α,β,ρ}

$x_{k+1}^{\mathrm{\δ}}=x_k^{\mathrm{\δ}}+{\mathrm{\ωF}}^{\′}{\left(x_k^{\mathrm{\δ}}\right)}^T(y^{\mathrm{\δ}}-F(x_k^{\mathrm{\δ}})),---\left(24\right)$

根据计算出的纵波速度α、横波速度β和岩石密度ρ，利用下列各式得出整个工区的岩性物理参数，

$\mathrm{\σ}=\frac{0.5*{(\mathrm{\α}/\mathrm{\β})}^2-1}{{(\mathrm{\α}/\mathrm{\β})}^2-1},---\left(25\right)$

$\mathrm{\λ}\mathrm{\ρ}=I_p^2-2*I_s^2,---\left(26\right)$

λ/μ＝(α/β)²-2，(27)

在上式中，I_p为纵波阻抗ρα，I_s为横波阻抗ρβ，σ为泊松比，λ为拉梅常数，μ为岩石的剪切模量。

与现有技术相比，本发明具有以下优点：

1、本发明利用多个角度的REI_i数据，通过非线性Landweber-Fridman迭代法进行反演计算，获得准确的纵横波速比、纵横波阻抗比，以及泊松比等岩性参数。并且由于更大程度地使用了叠前数据中蕴含的信息，整个计算过程相对传统的算法更加稳定可靠；

2、使用多个入射角的弹性阻抗数据，充分利用了叠前数据的信息，且算法对初值不敏感，不必事先构建初始模型；

3、能准确求出纵横波速度比、纵横波阻抗和泊松比等岩性参数，计算过程稳定可靠；

4、在实际资料应用中有较好的抗噪能力和较高的计算效率。

附图说明

图1为本发明反演方法的流程图；

图2为测井流体参数和0°～40°的射线弹性阻抗正演数据示意图；

图3为将各角度的REI_i数据通过非线性Landweber-Fridman迭代法进行反演计算后得到的纵波阻抗I_p的反演结果示意图；

图4为将各角度的REI_i数据通过非线性Landweber-Fridman迭代法进行反演计算后得到的横纵波速比k的反演结果示意图；

图5为通过本发明的反演方法以识别储层的含油气性的示意图；

图6为图5的反演结果剖面图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。

如图1所示，本发明提供一种基于正则化的非线性地震叠前弹性参数反演方法，包括以下步骤：

步骤S1，确定多角度的射线弹性阻抗井曲线REI_i。

具体的，根据已知的测井资料与射线弹性阻抗REI利用下式得出多角度的射线弹性阻抗井曲线REI_i，其中，所述测井资料中包括纵波速度、横波速度、密度与入射角，

${\mathrm{REI}}_i=REI(\mathrm{\α},\mathrm{\β},\mathrm{\ρ},{\mathrm{\θ}}_i)=\frac{\mathrm{\ρ}\mathrm{\α}}{{\mathrm{cos\θ}}_i}{(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_i)}^4,---\left(16\right)$

在上式中，α为纵波速度，β为横波速度，ρ为岩石密度，θ_i为入射角。

步骤S2，建立整个工区不同角度的射线弹性阻抗模型。

具体的，引入已知的层位、断层资料以及上述多角度的射线弹性阻抗井曲线REI_i，以建立整个工区不同角度的射线弹性阻抗模型；

步骤S3，使用叠前角度道集，计算出上述多角度的部分叠加数据体。

步骤S4，将多角度的部分叠加数据体、连同多角度的射线弹性阻抗模型，通过稀疏脉冲反演算法进行反演，得出上述入射角对应的射线弹性阻抗数据体；

具体的，将多角度的部分叠加数据体、连同多角度的射线弹性阻抗模型，通过稀疏脉冲反演算法进行反演，得出多角度的射线弹性阻抗数据体。

步骤S5，将上述入射角对应的射线弹性阻抗数据体，导入非线性Landweber-Fridman迭代算法，以得出整个工区的纵、横波速度和岩石密度。

步骤S6，根据纵波速度α、横波速度β和岩石密度ρ，利用(25)～(27)各式得出整个工区的各类岩性物理参数数据体。

图2为测井流体参数(左图)和0°到40°的射线弹性阻抗(REI)正演数据(右图)。

图3与图4分别是将各角度的射线弹性阻抗(REI)数据同时导入L-F算法进行计算，以得到的纵波阻抗和波速比反演结果。由计算结果可知，L-F算法具有很高的稳定性和准确性，且对初值不敏感，是一种可行的方法。

如图5与图6为本发明反演方法的应用效果，经理论和实践均证实，储层中含有流体会引起纵波速度降低、横波速度不变，使得Vp/Vs值降低，因此可以通过Vp/Vs来识别储层的含油气性。

图5为Vp/Vs反演结果的剖面图，其中，低值区域(红色区域)与测井识别的含气储层一致性较高。

图6为Vp/Vs反演结果切片示意图，其中，低值区域(红色)与测井识别的含气储层一致性较高，且相对关系体现更为明显。

本发明充分地利用叠前数据中蕴含的信息，并使用高精度的非线性反演方法，获得更加准确的纵横波速比、拉梅系数、剪切模量和泊松比等多种参数，提高对岩性和流体识别的置信度。同时改善现有EI和REI反演技术的不稳定性、低抗噪性等问题。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于正则化的非线性地震叠前弹性参数反演方法，包括以下步骤：

根据已知的测井资料与射线弹性阻抗REI利用下式得出多角度的射线弹性阻抗井曲线REI_i，其中，所述测井资料中包括纵波速度、横波速度、密度与入射角，

${\mathrm{REI}}_i=\mathrm{REI}(\mathrm{\α},\mathrm{\β},\mathrm{\ρ},{\mathrm{\θ}}_i)=\frac{\mathrm{\ρ\α}}{\cos {\mathrm{\θ}}_i}{(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_i)}^4,---\left(16\right)$

在上式中，α为纵波速度，β为横波速度，ρ为岩石密度，θ_i为入射角；

引入已知的层位、断层资料以及上述多角度的射线弹性阻抗井曲线REI_i，以建立整个工区不同角度的射线弹性阻抗模型；

使用叠前偏移距道集，计算出上述不同入射角的部分叠加数据体；

将多角度的部分叠加数据体、连同多角度的射线弹性阻抗模型，通过稀疏脉冲反演算法进行反演，得出上述入射角对应的射线弹性阻抗数据体；

将上述入射角对应的射线弹性阻抗数据体，导入非线性Landweber-Fridman迭代算法，以得出整个工区的纵、横波速度和岩石密度；

根据纵波速度α、横波速度β和岩石密度ρ，利用下式得出整个工区的岩性物理参数，

$\mathrm{\σ}=\frac{0.5*{(\mathrm{\α}/\mathrm{\β})}^2-1}{{(\mathrm{\α}/\mathrm{\β})}^2-1},---\left(25\right)$

$\mathrm{\λ\ρ}=I_p^2-2*I_s^2,---\left(26\right)$

λ/μ＝(α/β)²-2，(27)

在上式中，I_p为纵波阻抗ρα，I_s为横波阻抗ρβ，σ为泊松比，λ为拉梅常数，μ为岩石的剪切模量。

2.根据权利要求1所述的基于正则化的非线性地震叠前弹性参数反演方法，其特征在于，所述射线弹性阻抗REI通过下式获得：

$\mathrm{REI}(\mathrm{\α},\mathrm{\β},\mathrm{\ρ},\mathrm{\θ})=\frac{\mathrm{\ρ\α}}{\cos \mathrm{\θ}}{(1-\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})}^4,---\left(5\right)$

在上式中，α为纵波速度，β为横波速度，ρ为岩石密度，θ为入射角。

3.根据权利要求1所述的基于正则化的非线性地震叠前弹性参数反演方法，其特征在于，将多角度的射线弹性阻抗数据体导入非线性Landweber-Fridman迭代算法进行反演计算，以得到下式标准的非线性Landweber-Fridman迭代公式，

$x_{k+1}^{\mathrm{\δ}}=x_k^{\mathrm{\δ}}+\mathrm{\ω}{F}^{\′}{\left(x_k^{\mathrm{\δ}}\right)}^T(y-F\left(x_k^{\mathrm{\δ}}\right)).---\left(23\right)$

4.根据权利要求3所述的基于正则化的非线性地震叠前弹性参数反演方法，其特征在于，当右端项y不能够精确获得，而只是仅仅知道y^δ时，同样可以利用下式带扰动的非线性Landweber-Fridman迭代公式获得准确结果，

$x_{k+1}^{\mathrm{\δ}}=x_k^{\mathrm{\δ}}+\mathrm{\ω}{F}^{\′}{\left(x_k^{\mathrm{\δ}}\right)}^T(y^{\mathrm{\δ}}-F\left(x_k^{\mathrm{\δ}}\right)).---\left(24\right)$