CN111175821B - 一种vti介质的各向异性参数分步反演方法 - Google Patents

Abstract

本发明公开了一种VTI介质的各向异性参数分步反演方法，包括：获取研究工区的叠前地震道集数据；利用获取的道集数据，建立目标参数刚度参数的初始模型；基于初始模型，利用VTI等效介质反射系数精确式模拟合成记录；基于测井初始模型的参数值，计算VTI介质精确方程正演算子的偏导数；建立反演目标泛函，同时计算目标函数对模型参数；基于正演算子偏导数和目标函数一阶导数，计算模型更新梯度，并更新模型；重复进行反演迭代，直到模型误差降到预设范围，迭代停止，输出参数的反演结果；利用反演结果向量，通过数学关系换算各向异性参数。

Description

一种VTI介质的各向异性参数分步反演方法

技术领域

本发明属于非常规储层地震勘探技术领域，属于地震资料多参数反演范畴，具体是一种VTI介质的各向异性参数分步反演方法，用于提高各向异性参数的反演精度。

背景技术

随着对地下介质认知的深入，地质勘探人员开始意识到地下介质并不是各向同性的，而往往是各向异性的。在多种各向异性类型中，具有垂直对称轴的横向各向同性介质(Transverse isotropy with vertical axis of symmetry,VTI)是地下介质中存在最普遍的各向异性。常见的VTI介质包括有：薄的层状结构引起的结构各向异性VTI储层，以及矿物颗粒定向排列引起的岩石各向异性VTI储层。页岩等非常规储层为强VTI(垂直各向同性)各向异性特征。因此，开展适用于VTI介质的反演方法研究，用以挖掘页岩等非常规储层的弹性及各向异性信息来准确表征储层特征非常重要。

虽然，针对VTI介质的正演研究较多(如文献1～4公开的现有技术)，但是针对VTI介质的反演研究起步较晚。目前已存在的VTI介质反演主要是基于Rüger近似表达式开展研究(如文献5、6公开的现有技术)。侯栋甲等(2014)将Rüger近似式作为反演的正演算子，并采用贝叶斯方法建立目标函数开展了五参数的多波联合各向异性同时反演(如文献7公开的现有技术)。Zhang等(2019)基于Rüger近似式提出了两部反演策略(如文献8公开的现有技术)。虽然近似表达式数学形式简单物理意义明确，其在中大角度存在误差，而各向异性参数主要对中大角度数据敏感，因此选用近似式会增加各向异性参数估算的误差。此外，近似式还受到其弱反射假设的限制。VTI等效介质反射系数精确表达式已经被多位学者推导给出(如文献9、10公开的现有技术)，虽然数学表达式复杂，但是基于波动方程的精确解析解，可以精确模拟VTI介质的反射情况。正演方法还是各类反演算法的基础，反演算法中选用的正演算子模拟程度越高，越能获得准确的反演结果。因此，开展基于反射系数精确式的反演方法研究，以提高参数估计的精度，对页岩等非常规储层的参数估算和属性预测至关重要。

上述引用非专利文献如下：

[1]Wright J.The effect of transverse isotropy on reflection amplitudeversus offset[J].Geophysics,1987,52(4):564-567.

[2]Hron F,Daley PF.Reflection and transmission coefficients fortransversely isotropic media[J].Bulletin of the Seismological Society ofAmerica,1977,67(3):661-675.

[3]Daley P F,Hron F.Reflection and transmission coefficients forseismic waves in ellipsoidally anisotropic media[J].Geophysics,1979,44(1):27-38.

[4]Thomsen L A.Weak anisotropic reflections,Offset-DependentReflectivity:Theory and Practice of AVO Analysis[M].Library.seg.org,1993:103–114

[5]Rüger A.P-wave reflection coefficients for transversely isotropicmodels with vertical and horizontal axis of symmetry[J].Geophysics,1997,62(3):713-722.

[6]Rüger A.Variation of P-wave reflectivity with offset and azimuthin anisotropic media[J].Geophysics,1998,63(3):935-947.

[7]侯栋甲,刘洋,任志明,等.基于贝叶斯理论的VTI介质多波叠前联合反演[J].石油物探,2014,53(3):294-303.

[8]Zhang,F.,Zhang,T.,andLi,X.Seismic amplitude inversion for thetransversely isotropic media with vertical axis of symmetry:GeophysicalProspecting,2019,67,2368-2385.

[9]Graebner M.Plane-wave reflection and transmission coefficients fora transversely isotropic solid[J].Geophysics,1992,57(11):1512-1519.

[10]寻浩，董敏煜，牟永光.横向各向同性介质中的AVO[J].石油地球物理勘探，1997，32(1):45-56.

发明内容

本发明的目的是针对现有技术存在的问题，提供一种VTI介质的各向异性参数分步反演方法，可以提高各向异性参数的反演精度。

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

一种VTI介质的各向异性参数分步反演方法，包括以下步骤：

步骤一，获取研究工区的叠前地震道集数据，并从所述叠前地震道集数据中提取不同角度的地震子波序列；

获取工区PP数据的叠前偏移距道集，基于3D速度数据，将叠前偏移距数据转换为叠前角度道集；

从叠前道集中提取地震子波序列，建立不同角度下的子波矩阵W(θ)，其为一个增广矩阵：

其中，w(θ_i)为角度为θ_i时的子波矩阵，由子波序列组成：

步骤二，获取包含各向异性参数的测井数据，将各向异性参数转换为刚度参数，建立刚度参数的初始模型；

获取井口处的弹性参数和各向异性参数曲线，包括垂直方向纵波速度α₀、垂直方向横波速度β₀、密度ρ、Thomsen各向异性参数ε和Thomsen各向异性参数δ；这五参数为反演的目标参数，目标参数向量m_o为：

m_o＝[α¹ … αⁿ β¹ … βⁿ ρ¹ … ρⁿ ε¹ … εⁿ δ¹ … δⁿ] (3)

然而目标五参数敏感性差异过大，因此反演参数转为新参数组合，包括刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ。新参数向量m，具体形式为：

利用数学关系，将上述原始目标五参数组合转为新参数组合，具体表达式如下：

将新目标参数曲线通过Backus平均算法获取低频数据信息；

基于井口低频数据，利用层位解释结果做约束，插值获得2D/3D的低频初始模型。

步骤三，利用VTI等效介质反射系数精确式模拟合成记录，计算合成记录和实际数据的残差；

首先计算等效VTI介质的反射系数序列，精确方程表达式如下：

r＝S^-1·b (7)

其中，S和b矩阵表达式如下：

其中，上标(*)^T和(*)^B分别表示了上层和下层介质的属性参数。下标(*)_P和(*)_S分别表示与纵波和横波相关的属性参数。式中，符号l、n、a、b、d和e均为中间变量参数，且均为刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ的函数；具体表达式如下：

式中，l_P、l_S、n_P和n_S为极化向量的方向余弦，表达式如下：

式(9)至(17)中，q_P和q_S分别代表P波和S波的垂直慢度，为刚度参数和密度的函数。

VTI等效介质精确公式模拟为非线性正演方法，基于褶积理论，利用子波矩阵W和反射系数序列r，可计算合成记录d_syn：

d_syn＝G(m)＝W·r (18)

其中，G表示精确公式的正演过程，模型参数向量m的函数；d_syn为正演合成数据，W为子波矩阵，式(18)具有如下矩阵形式：

式中，d(θ_i)为入射角θ_i对应的合成数据，如下：

d_syn(θ_i)＝[d_syn(t₁,θ_i) d_syn(t₂,θ_i) … d_syn(t_n,θ_i)]^Τ (20)

从入射角θ_i的实际地震道集中抽取子波，构建子波矩阵w(θ_i)，如式(2)所示。r(θ_i)为入射角θ_i时计算得到的反射系数序列，通过式(19)计算得到：

r(θ_i)＝[r(t₁,θ_i) r(t₂,θ_i) … r(t_n,θ_i)]^Τ (21)

计算实际道集数据d_obs和正演模拟记录d_syn的残差Δd：

Δd＝d_obs-d_syn＝d_obs-G(m) (22)

步骤四，基于测井初始模型的参考值，计算VTI介质精确方程正演算子对模型参数的偏导数；

精确式正演算子G(m)对模型参数m的导数为：

反射系数r对模型参数的导数为：

其中，

和

均为中间变量矩阵对模型参数的偏导数；

其中，针对上层模型参数的导数如下：

其中，

针对下层介质模型参数，其偏导数为：

其中，

步骤五，建立反演目标函数，计算目标函数的一阶导数；

建立反演目标函数，具体表达式为：

其中，C为模型参数的协方差矩阵，u为模型参数的期望；

目标函数J(m)对模型参数m的一阶导数具体表达式如下：

步骤六，基于正演算子偏导数和目标函数一阶导数，计算模型更新梯度，并利用迭代获取模型更新结果；

选用高斯牛顿的优化算法，伪海森矩阵的形式为：

针对第k次迭代，模型更新梯度表达式为：

g_k＝H(m_k)^-1·γ(m_k) (38)

设定负梯度为模型的更新方向，模型参数的迭代更新公式为：

m_k+1＝m_k-α·g_k (39)

步骤七，重复步骤三至步骤六，进行下一次反演迭代，直到模型的数据残差降低到预设范围内，迭代停止，输出刚度参数的反演结果；

步骤八，利用刚度参数的反演结果，获取最终各向异性参数的估算结果；

反演目标参数组合包括刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，利用数学关系换算各向异性参数，包括Thomsen各向异性参数ε和Thomsen各向异性参数δ；具体数学表达式如下：

ε＝(c₁₁-c₃₃)/2c₃₃ (40)

与现有技术相比，本发明的有益效果是：本发明通过将VTI介质精确反射系数引入各向异性参数反演中，推导了精确公式对于模型参数的偏导数，选用了线性反演策略以提高反演时效，相较于传统的基于Rüger近似式的各向异性参数反演方法，本发明的基于反射系数精确式的VTI介质反演方法提高了各向异性参数的反演精度。

附图说明

图1为本发明一种VTI介质的各向异性参数分步反演方法的流程示意图；

图2为本发明实施例中测井中获得的弹性参数和各向异性参数曲线图；

图3为本发明实施例中根据测井数据计算得到的新参数曲线图；

图4为本发明实施例中基于反射系数精确式的叠前反演方法获得的新目标参数反演结果示意图；

图5为本发明实施例中通过数学公式换算得到的Thomsen各向异性参数曲线图。

具体实施方式

下面将结合本发明中的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例，都属于本发明保护的范围。

如图1所示，本实施例提供了一种基于等效VTI介质反射系数精确式的各向异性参数反演方法，具体包括以下步骤：

一、基于VTI介质精确反射系数公式的正演：

获取井口处的弹性参数和各向异性参数曲线，包括垂直方向纵波速度α₀、垂直方向横波速度β₀、密度ρ、Thomsen各向异性参数ε和Thomsen各向异性参数δ。这五参数为反演的目标参数，目标参数向量m_o为：

m_o＝[α¹ … αⁿ β¹ … βⁿ ρ¹ … ρⁿ ε¹ … εⁿ δ¹ … δⁿ] (1)

然而目标五参数敏感性差异过大，因此反演参数转为新参数组合，包括刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ。新参数向量m，具体形式为：

利用数学关系，将上述原始目标五参数组合转为新参数组合，具体表达式如下：

将新目标参数曲线通过Backus平均算法获取低频数据信息；利用层位解释结果做约束，插值获得低频初始模型。

从叠前道集中提取地震子波序列，建立不同角度下的子波矩阵W(θ)，其为一个增广矩阵：

其中，w(θ_i)为角度为θ_i时的子波矩阵，由子波序列组成：

计算得到正演合成数据d：

d＝W·r (7)

其有如下矩阵形式：

式中，入射角θ_i所对应的数据为d(θ_i)，具体数学表达式为：

d(θ_i)＝[d(t₁,θ_i) d(t₂,θ_i) … d(t_M,θ_i)]^Τ (9)

对应的反射系数序列为：

r(θ_i)＝[r(t₁,θ_i) r(t₂,θ_i) … r(t_M,θ_i)]^Τ (10)

等效VTI介质的反射系数序列可由如下的精确方程表达式计算得到：

r＝S^-1·b (11)

其中，S和b具体表达式为：

其中，上标(*)^T和(*)^B分别表示了上层和下层介质的属性参数。下标(*)_P和(*)_S分别表示与纵波和横波相关的属性参数。式中，符号l、n、a、b、d和e均为中间变量参数，且均为刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ的函数。具体表达式如下：

式中，l_P、l_S、n_P和n_S为极化向量的方向余弦，表达式如下：

式(14)至(21)中，q_P和q_S分别代表P波和S波的垂直慢度，为刚度参数和密度的函数。

二、基于VTI介质精确反射系数公式的叠前反演方法：

建立反演目标函数，具体表达式为：

其中，C为模型参数的协方差矩阵，u为模型参数的期望。

目标函数J(m)对模型参数m的一阶导数具体表达式如下：

设定负梯度为模型的更新方向，模型参数的迭代更新公式为：

m_k+1＝m_k-α·g_k (24)

其中，针对第k次迭代，模型更新梯度表达式为：

g_k＝H(m_k)^-1·γ(m_k) (25)

本发明选用高斯牛顿的优化算法，其伪海森矩阵的形式为：

其中，

为正演算子对模型的偏导数。

当模型更新使得误差达到要求，则可停止迭代，输出模型结果。利用反演结果向量，通过数学关系换算各向异性参数，包括Thomsen各向异性参数ε和Thomsen各向异性参数δ。具体数学表达式如下：

ε＝(c₁₁-c₃₃)/2c₃₃ (27)

δ＝[(c₁₃+c₅₅)²-(c₃₃-c₅₅)²]/[2c₃₃(c₃₃-c₅₅)] (28)

三、精确反射系数公式对模型的偏导数：

精确式正演算子G(m)对模型参数m的导数

可由反射系数公式对模型导数求取：

反射系数r对模型参数的导数：

其中，

和

均为中间变量矩阵对模型参数的偏导数。其中，针对上层模型参数的导数如下：

其中，

针对下层介质模型参数，其偏导数为：

其中，

图2(a～e)为本实施例测井中获得的弹性参数和各向异性参数数据，包括纵波速度、横波速度、密度、Thomsen各向异性参数ε和各向异性参数δ。图3(a～e)是计算获得的新参数组合，包括有刚度参数C₃₃、C₅₅、C₁₁、C₁₃和密度。图4(a～e)为基于精确反射系数公式的五参数同时反演结果，其中虚线为反演初始模型，实线为真实模型，点线为反演结果；从图4中可以看出，C₁₁、C₃₃、C₅₅反演结果比较好，C₁₃反演效果较差，密度反演结果最差。图5(a、b)为利用图4反演结果换算而得到的Thomsen各向异性参数结果，实线为真实模型，点线为反演结果。虽然图4中密度结果很差，但是并不影响各向异性参数的估算结果，由于C₁₁和C₃₃反演结果较理想，则可获得较好的各向异性参数ε结果。而C₁₃效果不佳则影响了各向异性参数δ。说明该方法可以获得较好的Thomsen各向异性参数ε的估算结果。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (8)

1.一种VTI介质的各向异性参数分步反演方法，其特征在于，包括以下步骤：

步骤一，获取研究工区的叠前地震道集数据，并从所述叠前地震道集数据中提取不同角度的地震子波序列；

步骤二，获取包含各向异性参数的测井数据，将各向异性参数转换为刚度参数，建立刚度参数的初始模型；

步骤三，利用VTI等效介质反射系数精确式模拟合成记录，计算合成记录和实际数据的残差；

步骤四，基于初始模型参数，计算精确方程正演算子对模型参数的偏导数；

步骤五，建立反演目标函数，计算目标函数的一阶导数；

步骤六，计算模型更新梯度，利用迭代获取模型更新结果；

步骤七，重复步骤三至步骤六，直到模型的数据残差降低到预设范围内，迭代停止，输出刚度参数的反演结果；

步骤八，利用刚度参数的反演结果，获取最终各向异性参数的估算结果。

2.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤一具体包括：

获取工区PP数据的叠前偏移距数据，基于速度数据，将所述叠前偏移距数据转换为叠前角道集数据；

从所述叠前角道集数据中提取地震子波序列，建立不同角度下的子波矩阵。

3.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤二具体包括：

获取井口处的弹性参数和各向异性参数，所述各向异性参数包括垂直方向纵波速度α₀、垂直方向横波速度β₀、密度ρ、Thomsen各向异性参数ε和Thomsen各向异性参数δ；

将所述各向异性参数的组合转换为包含刚度参数的新数据组合，所述新数据组合包括刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，具体表达式如下：

通过Backus平均算法获取新目标参数曲线的低频数据信息；

基于井口低频数据，利用层位解释结果做约束，插值获得2D/3D的低频初始模型。

4.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤三具体包括：

利用VTI等效介质反射系数精确方程计算反射系数序列：

r＝S^-1·b

其中，S和b矩阵表达式如下：

其中，上标(*)^T和(*)^B分别表示了上层和下层介质的属性参数，下标(*)_P和(*)_S分别表示与纵波和横波相关的属性参数；式中，符号l、n、a、b、e均为中间变量参数，且均为刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ的函数；

利用子波矩阵W和反射系数序列r，计算合成记录d_syn：

d_syn＝G(m)＝W·r

计算实际道集数据d_obs和正演模拟记录d_syn的残差Δd：

Δd＝d_obs-d_syn＝d_obs-G(m)

其中，G(m)为精确式正演算子，m为目标模型参数向量。

5.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤四中，计算精确方程正演算子对模型参数的偏导数的方法为：

目标模型参数矩阵m具有如下形式：

精确式正演算子G(m)对模型参数m的导数为：

反射系数r对模型参数的导数为：

其中，

和

均为中间变量矩阵对模型参数的偏导数；c₁₁、c₁₃、c₃₃、c₅₅均为刚度参数；ρ为密度参数；W为子波矩阵；b、S为计算反射系数序列r的中间过程矩阵。

6.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤五具体包括：

建立反演目标函数，具体表达式为：

其中，C为模型参数的协方差矩阵，u为模型参数的期望；Δd为实际道集数据d_obs和正演模拟记录d_syn的残差；d_obs为实际地震道集数据；

目标函数J(m)对模型参数m的一阶导数具体表达式如下：

其中，λ为正则化参数。

7.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤六具体包括：

选用高斯牛顿的优化算法，伪海森矩阵的形式为：

针对第k次迭代，模型更新梯度表达式为：

g_k＝H(m_k)^-1·γ(m_k)

设定负梯度为模型的更新方向，模型参数的迭代更新公式为：

m_k+1＝m_k-α·g_k

其中，J(m)为反演目标函数；G(m)为精确式正演算子；λ为正则化参数；C为模型参数的协方差矩阵；γ(m_k)为第k次迭代时，反演目标函数J(m)对目标模型参数m的一阶导数；α·g_k为模型的更新量，其中α为模型更新步长，g_k为模型更新梯度。

8.根据权利要求1所述的一种VTI介质的各向异性参数分步反演方法，其特征在于，步骤八具体包括：

反演目标参数组合包括刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，利用数学关系换算各向异性参数，包括Thomsen各向异性参数ε和Thomsen各向异性参数δ；具体数学表达式如下：

ε＝(c₁₁-c₃₃)/2c₃₃

Images