CN111308550A - 一种页岩vti储层的各向异性参数多波联合反演方法 - Google Patents

Abstract

本发明公开了一种页岩VTI储层的各向异性参数多波联合反演方法，所述方法包括如下步骤：获取地震子波序列，并建立地震子波矩阵；通过获取的测井数据，建立初始模型；通过初始模型和地震子波矩阵计算PP(反射纵波)和PS(转换横波)合成数据，以及合成记录和实际数据的残差；计算PP和PS数据正演过程对模型参数的正演偏导数；建立反演目标泛函，根据数据残差结果和正演过程偏导数计算目标泛函对模型参数的泛函偏导数；根据泛函偏导数计算初始模型更新梯度，并更新初始模型；通过更新后的初始模型得出的反演结果换算各向异性参数。

Description

一种页岩VTI储层的各向异性参数多波联合反演方法

技术领域

本发明涉及页岩各向异性储层地震勘探技术领域，具体涉及一种页岩VTI储层的各向异性参数多波联合反演方法。

背景技术：

地震反演是利用地下目标地质体的地震波响应信息来获得目标介质构造形态、岩石物理以及含流体属性相关的弹性信息的方法，是地球物理学者挖掘未知储层信息的重要途径。对于以往常规简单构造的油气储层，地震反演一直基于各向同性介质假设，并形成了一套完善的油气储层检测技术体系。然而，随着全球范围内对石油和天然气需求的不断增长，地球物理学家们开始着手探索不同类型的油气藏，油气勘探方向逐步由常规油气储层向一些特殊的非常规油气储层。页岩等非常规储层成为了目前油气勘探的热点和难点。然而对于这类储层的研究继续沿用各向同性假设理论是不合理且行不通的。

各向异性介质中，具有垂向对称轴的横向各向同性介质，即VTI介质是最普遍存在的类型。例如结构各向异性VTI(Transverse isotropy with veritcal axis ofsymmetry)储层，即沉积地层中层状结构所形成的VTI结构；又如岩石各向异性VTI储层，即岩石的矿物颗粒定向排列所形成的VTI结构。页岩属于典型的VTI介质储层。以往用于常规储层的反演方法并不适用于这类VTI介质。因此，开展适用于VTI介质的反演方法研究，对于储层参数的准确估算至关重要。众多学者对VTI介质的正演模拟算法进行了研究，包括反射系数精确公式以及近似表达式。基于VTI介质的反演算法大多选用了近似表达式作为正演算子。但近似式在模拟精度和弱反射等方面限制了相应的反演准确性。因此，基于反射系数精确表达式的反演有助于提高反演精度。但是VTI反演需同时估算五个参数，多参数同时反演需克服反演稳定性差等问题。因此，采用多波联合的反演方法有助于提升反演稳定性，同时提高估算精度。

发明内容

本发明的目的在于提供一种页岩VTI储层的各向异性参数多波联合反演方法，以解决现有技术中导致的各向异性参数反演估算精度差的缺陷。

一种页岩VTI储层的各向异性参数多波联合反演方法，所述方法包括如下步骤：

获取地震子波序列，并建立地震子波矩阵；

通过获取的测井数据，建立初始模型；

通过初始模型和地震子波矩阵计算PP和PS合成数据，以及合成数据残差；

计算PP和PS数据正演过程对模型参数的正演偏导数；

建立反演目标泛函，根据合成数据残差和正演偏导数，计算目标泛函对模型参数的泛函偏导数；

根据泛函偏导数计算初始模型更新梯度，并更新初始模型；

通过更新后的初始模型得出的反演结果换算各向异性参数。

进一步的，获取地震子波序列，并建立地震子波矩阵的方法包括如下步骤：

获取研究区的PP和PS叠前偏移距数据，利用三维速度体，将其转为PP和PS角道集数据；

从PP和PS角道集中提取用于PP和PS数据合成的地震子波序列，建立包含不同角度子波序列的子波矩阵W_PP(θ)和W_PS(θ)。

进一步的，通过获取的测井数据，建立初始模型的方法包括如下步骤：

获取研究区井口处的原始目标参数组合，包括弹性参数：垂向纵波速度V_P0、垂向横波速度V_S0、密度ρ以及各向异性参数ε和δ；

模型参数向量为：

新模型参数组合包括四个刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，模型参数向量为：

新模型参数组合与原模型参数组合之间的具体数学关系如下：

对新模型参数数据进行低通滤波处理，基于层位数据约束，利用插值处理获得新参数模型组合的低频初始模型。

进一步的，通过初始模型和地震子波矩阵计算PP和PS合成数据，以及合成数据残差的方法包括如下步骤：

基于褶积理论，由子波矩阵和反射系数序列褶积获得合成记录：

d_PP＝G_PP(m)＝W_PP·R_PP；

d_PS＝G_PS(m)＝W_PS·R_PS；

其中反射系数矩阵为r＝[R_PP R_PS T_PP T_PS]，可利用VTI等效介质的Graebner反射系数精确方程计算获得：

S·r＝b；

其具体矩阵表达式如下：

其中，上标T和B代表反射界面上下层介质的弹性参数或中间变量；下标P和S参数表示与纵波和横波相关的弹性参数及相关中间变量；此外，p表示纵波水平慢度，s_*为垂直慢度；参数

n、a、b、e均为数学中间变量；以上参数均可由刚度参数c_*和密度ρ计算求得；设定d_PP和d_PS为实际PP和PS角道集数据，则两种反射波的实际与合成数据残差为：

△d_PP＝d_PP-G_PP(m)；

△d_PS＝d_PS-G_PS(m)。

进一步的，计算PP和PS数据正演过程对模型参数的正演偏导数的方法包括如下步骤：

PP和PS反射波正演过程G_PP(m)和G_PS(m)对模型参数m的偏导数如下：

其中，模型参数向量m包括：

m＝[c₁₁ c₁₃ c₃₃ c₅₅ ρ]；

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

其中，

和

为中间变量矩阵对不同模型参数的偏导数：

进一步的，建立反演目标泛函，根据合成数据残差和正演偏导数，计算目标泛函对模型参数的泛函偏导数的方法包括如下步骤：

基于正则化理论，建立目标泛函：

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

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

进一步的，根据泛函偏导数计算初始模型更新梯度，并更新初始模型的方法包括如下步骤：

针对线性反演策略，模型的更新迭代格式如下：

△_km＝m_k+1-m_k＝-β·g_k；

其中，k表示第k次反演迭代，β表示更新步长，g_k为第k次迭代的更新梯度，负梯度方向即为模型的更新方向；选用Guass-Newton优化算法，梯度计算式如下：

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

其中，

进一步的，通过更新后的初始模型得出的反演结果换算各向异性参数的方法包括如下步骤：

反演模型组合为{c₁₁,c₁₃,c₃₃,c₅₅,ρ}，可利用数学关系换算两个Thomsen各向异性参数ε和δ；具体数学表达式如下：

本发明的优点在于：该种具有垂向对称轴的横向各向同性VTI储层的各向异性参数多波联合反演方法，选用了Graebner精确表达式作为反演的正演算子来提高各向异性参数的反演精度，以此来克服以往的近似表达式所带来的存在模拟误差等问题；将PP(纵波-纵波)、PS(纵波-转换横波)单数据反演法拓展为PP和PS多波数据联合反演，选用多波数据可提升参数的敏感性，同时提高反演的稳定性；选用测井模型数据进行方法测试，并对比了PP数据和PP-PS联合数据的反演结果；测试表明联合反演明显提高了刚度参数的反演精度，相应地提高了各向异性参数的估算准确性。

附图说明

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

图2为井口处获得的原始弹性参数组合；

图3为井口处通过数学关系换算得到的新弹性参数组合；

图4为刚性参数C₃₃的PP数据和PP-PS联合数据反演结果；

图5为刚性参数C₅₅的PP数据和PP-PS联合数据反演结果；

图6为刚性参数C₁₁的PP数据和PP-PS联合数据反演结果；

图7为刚性参数C₁₃的PP数据和PP-PS联合数据反演结果；

图8为密度参数ρ的PP数据和PP-PS联合数据反演结果；

图9为PP数据和PP-PS联合数据反演结果换算得到的Thomsen各向异性参数ε；

图10为PP数据和PP-PS联合数据反演结果换算得到的Thomsen各向异性参数δ。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

如图1至图10所示，一种页岩VTI储层的各向异性参数多波联合反演方法，所述方法包括如下步骤：

步骤一：获取地震子波序列，并建立地震子波矩阵：

获取研究区的PP和PS叠前偏移距数据，利用三维速度体，将其转为PP和PS角道集数据；

从PP和PS角道集中提取用于PP和PS数据合成的地震子波序列，建立包含不同角度子波序列的子波矩阵W_PP(θ)和W_PS(θ)；

搜集工区的PP和PS角道集数据d_PP(θ)和d_PS(θ)；从地震角度道集中按照不同角度提取子波序列，建立子波矩阵W_PP(θ)和W_PS(θ)：

该矩阵为一个增广矩阵，PP和PS数据道集均共有m个入射角度。下标(*)_PP和(*)_PS分别表示PP和PS数据对应的子波参数。w(θ_i)为θ_i入射角地震道抽取的地震子波序列建立的矩阵：

其中，PP数据提取的子波序列长度为L1，PS数据子波序列长度为L2。

步骤二：通过获取的测井数据，建立初始模型：

获取研究区井口处的原始目标参数组合，包括弹性参数：垂向纵波速度V_P0、垂向横波速度V_S0、密度ρ以及各向异性参数ε和δ；

模型参数向量为：

新模型参数组合包括四个刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，模型参数向量为：

新模型参数组合与原模型参数组合之间的具体数学关系如下：

对新模型参数数据进行低通滤波处理，基于层位数据约束，利用插值处理获得新参数模型组合的低频初始模型；

搜集井口处的弹性参数纵波速度V_P0、横波速度V_S0、密度ρ，以及弱各向异性参数ε和δ，即参数组合{V_P0,V_S0,ρ,ε,δ}。设定反演参数组合为{c₁₁,c₁₃,c₃₃,c₅₅,ρ}，其中c₁₁、c₁₃、c₃₃和c₅₅均为刚度参数，可以由速度和各向异性参数换算求得，具体表达式如下：

利用数学关系换算得到反演参数组合，进行低通滤波后，插值处理建立低频的初始模型用于反演过程。

步骤三：通过初始模型和地震子波矩阵计算PP和PS合成数据，以及合成数据残差：

基于褶积理论，由子波矩阵和反射系数序列褶积获得合成记录：

其中反射系数矩阵为r＝[R_PP R_PS T_PP T_PS]，可利用VTI等效介质的Graebner反射系数精确方程计算获得：

S·r＝b (8)

其具体矩阵表达式如下：

其中，参数f^T和g^T满足

和

上标T和B代表反射界面上下层介质的弹性参数或中间变量。下标P和S参数表示与纵波和横波相关的弹性参数及相关中间变量。此外，s_x表示水平慢度，s_z为垂直慢度；参数

n、a、b、e均为数学中间变量。以上参数均可由刚度参数c_*和密度ρ计算求得。

具体表达式如下：

其中，

针对横波相关的参数：

其中，

其他中间变量为：

设定d_PP和d_PS为实际PP和PS角道集数据，则两种反射波的实际与合成数据残差为：

步骤四：计算PP和PS数据正演过程对模型参数的正演偏导数；

PP和PS反射波正演过程G_PP(m)和G_PS(m)对模型参数m的偏导数如下：

其中，模型参数向量m包括：

m＝[c₁₁ c₁₃ c₃₃ c₅₅ ρ] (23)

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

其中，

和

为中间变量矩阵对不同模型参数的偏导数：

其中，上标(*)¹和(*)²分别对上层和下层介质参数求导数。矩阵S对上层参数的导数如下：

矩阵S对下层参数的导数如下：

矩阵b对上层参数的导数如下：

矩阵b对下层参数的导数如下：

式中，

和

均是中间参数变量对模型参数的导数。式中，m^*包括有

和ρ^*＝T,B。

步骤五：建立反演目标泛函，根据合成数据残差和正演偏导数，计算目标泛函对模型参数的泛函偏导数：

基于正则化理论，建立目标泛函：

其中，C_m为模型向量的协方差矩阵，u为模型向量的期望。

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

步骤六：根据泛函偏导数计算初始模型更新梯度，并更新初始模型；

针对线性反演策略，模型的更新迭代格式如下：

△_km＝m_k+1-m_k＝-β·g_k (36)

其中，k表示第k次反演迭代，β表示更新步长，g_k为第k次迭代的更新梯度，负梯度方向即为模型的更新方向。选用Guass-Newton优化算法，梯度计算式如下：

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

其中，

步骤七：利用步骤六中更新模型参数计算新的合成记录，并计算记录的误差精度；若误差精度满足要求，则停止迭代，输出模型参数结果；若误差精度不满足要求，则重复步骤三至六，直到满足误差精度或达到最大迭代次数。

步骤八：通过更新后的初始模型得出的反演结果换算各向异性参数：

反演模型组合为{c₁₁,c₁₃,c₃₃,c₅₅,ρ}，可利用数学关系换算两个Thomsen各向异性参数ε和δ。具体数学表达式如下：

在本实施例中，利用反射系数精确式计算正演记录的方法：

基于褶积理论，合成记录可以由子波矩阵和反射系数序列褶积获得：

其中R_PP和R_PS表示为PP和PS反射系数矩阵。式(1)和(2)有如下矩阵形式：

式中，入射角θ_i的合成记录表示为

和

d_PP(θ_i)＝[d_PP(t₁,θ_i)d_PP(t₂,θ_i)…d_PP(t_M,θ_i)]^T (45)

d_PS(θ_i)＝[d_PS(t₁,θ_i)d_PS(t₂,θ_i)…d_PS(t_M,θ_i)]^T (46)

入射角θ_i的反射系数表示为

和

R_PP(θ_i)＝[R_PP(t₁,θ_i)R_PP(t₂,θ_i)…R_PP(t_M,θ_i)]^T (47)

R_PS(θ_i)＝[R_PS(t₁,θ_i)R_PS(t₂,θ_i)…R_PS(t_M,θ_i)]^T (48)

θ_i入射角和t_j时间的反射系数R_PP(t_j,θ_i)和R_PS(t_j,θ_i)可由Graebner反射系数精确方程：

S·r＝b (49)

其中，r＝[R_PP R_PS T_PP T_PS]^T为反射系数和透射系数矩阵。式(10)矩阵表达式如下：

其中，参数f^T和g^T满足

和

上标T和B代表反射界面上下层介质的弹性参数或中间变量。下标P和S参数表示与纵波和横波相关的弹性参数及相关中间变量。此外，s_x表示水平慢度，s_z为垂直慢度；参数

n、a、b、e均为数学中间变量。以上参数均可由刚度参数c_*和密度ρ计算求得。

具体表达式如下：

其中，

针对横波相关的参数：

其中，

其他中间变量为：

搜集到的井口处的弹性参数纵波速度α₀、横波速度β₀、密度ρ和各向异性参数ε和δ。利用数学公式计算刚度参数和密度，下为参数向量具体形式：

m＝[c₁₁,c₁₃,c₃₃,c₅₅,ρ]

[(c₁₁)₁,…,(c₁₁)_n,(c₁₃)₁,…,(c₁₃)_n,(c₃₃)₁,…,(c₃₃)_n,(c₅₅)₁,…,(c₅₅)_n,ρ₁,…,ρ_n]

(60)

利用参数向量可计算反射系数R_PP和R_PS。

在本实施例中，多波联合的各向异性反演方法

基于正则化理论，建立反演的目标泛函：

其中，C_m为模型向量的协方差矩阵，u为模型向量的期望。

针对线性反演策略，模型的更新迭代格式如下：

m_k+1＝m_k+△_km＝-β·g_k (62)

其中，k表示第k次反演迭代，β表示更新步长，g_k为第k次迭代的更新梯度，负梯度方向即为模型的更新方向。选用Guass-Newton优化算法，梯度计算式如下：

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

式中，γ(m_k)为第k次迭代时目标泛函J(m_k)对模型参数的一阶导数，具体表达式如下：

其中，H(m_k)为第k次迭代时目标泛函J(m_k)的伪海森矩阵，其具体表达为:

其中，

和

为PP和PS正演过程的导数。

利用更新模型参数计算新的合成记录，并计算记录的误差精度；若误差精度满足要求，则停止迭代，输出模型参数结果；若误差精度不满足要求，则开始下一迭代，直到满足误差精度或达到最大迭代次数。利用输出的反演结果，换算各向异性参数，具体数学表达式如下：

在本实施例中，正演算子对模型参数的导数的方法：

PP和PS反射波正演过程G_PP(m)和G_PS(m)对模型参数m的偏导数

和

如下：

其中，模型参数向量m包括：

m＝[c₁₁ c₁₃ c₃₃ c₅₅ ρ] (70)

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

其中，

和

为中间变量矩阵对不同模型参数的偏导数：

其中，上标(*)¹和(*)²分别对上层和下层介质参数求导数。矩阵S对上层参数的导数如下：

矩阵S对下层参数的导数如下：

矩阵b对上层参数的导数如下：

矩阵b对下层参数的导数如下：

式中，

和

均是中间参数变量对模型参数的导数。式中，m^*包括有

和ρ^*＝T,B。

图2为井口处搜集的原始弹性参数组合，包括纵波速度、横波速度、密度和两个各向异性参数；

图3为井口处通过数学关系换算得到的新弹性参数组合，包括四个刚度参数和密度参数；

图4为刚性参数C₃₃的PP数据和PP-PS联合数据反演结果；

图5为刚性参数C₅₅的PP数据和PP-PS联合数据反演结果；

图6为刚性参数C₁₁的PP数据和PP-PS联合数据反演结果；

图7为刚性参数C₁₃的PP数据和PP-PS联合数据反演结果；

图8为密度参数ρ的PP数据和PP-PS联合数据反演结果。

图4至8中，黑色虚线为初始模型，黑色实线为真实参数，灰色点线为反演结果。可以看出四个刚度参数采用PP-PS联合数据的反演结果明显优于PP数据反演结果。两种反演方法的密度结果都不理想，但是密度的结果并不影响各向异性参数的估算结果。

图9为PP数据和PP-PS联合数据反演结果换算得到的Thomsen各向异性参数ε；

图10为PP数据和PP-PS联合数据反演结果换算得到的Thomsen各向异性参数δ。黑色实线为真实参数，灰色点线为估算结果。PP数据的估算结果中，参数ε明显优于参数δ，参数δ估算结果并不理想；相比于PP反演结果，PP和PS联合反演明显提高了参数ε和δ的估算精度。测井数据测试结果说明了本发明的多波联合反演方法的有效性。

由技术常识可知，本发明可以通过其它的不脱离其精神实质或必要特征的实施方案来实现。因此，上述公开的实施方案，就各方面而言，都只是举例说明，并不是仅有的。所有在本发明范围内或在等同于本发明的范围内的改变均被本发明包含。

Claims (8)

1.一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于，所述方法包括如下步骤：

获取地震子波序列，并建立地震子波矩阵；

通过获取的测井数据，建立初始模型；

通过初始模型和地震子波矩阵计算PP和PS合成数据，以及合成数据残差；

计算PP和PS数据正演过程对模型参数的正演偏导数；

建立反演目标泛函，根据合成数据残差和正演偏导数，计算目标泛函对模型参数的泛函偏导数；

根据泛函偏导数计算初始模型更新梯度，并更新初始模型；

通过更新后的初始模型得出的反演结果换算各向异性参数。

2.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：获取地震子波序列，并建立地震子波矩阵的方法包括如下步骤：

获取研究区的PP和PS叠前偏移距数据，利用三维速度体，将其转为PP和PS角道集数据；

从PP和PS角道集中提取用于PP和PS数据合成的地震子波序列，建立包含不同角度子波序列的子波矩阵W_PP(θ)和W_PS(θ)。

3.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：通过获取的测井数据，建立初始模型的方法包括如下步骤：

获取研究区井口处的原始目标参数组合，包括弹性参数：垂向纵波速度V_P0、垂向横波速度V_S0、密度ρ以及各向异性参数ε和δ；

模型参数向量为：

新模型参数组合包括四个刚度参数c₁₁、c₁₃、c₃₃、c₅₅和密度ρ，模型参数向量为：

新模型参数组合与原模型参数组合之间的具体数学关系如下：

对新模型参数数据进行低通滤波处理，基于层位数据约束，利用插值处理获得新参数模型组合的低频初始模型。

4.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：通过初始模型和地震子波矩阵计算PP和PS合成数据，以及合成数据残差的方法包括如下步骤：

基于褶积理论，由子波矩阵和反射系数序列褶积获得合成记录：

d_PP＝G_PP(m)＝W_PP·R_PP；

d_PS＝G_PS(m)＝W_PS·R_PS；

其中反射系数矩阵为r＝[R_PP R_PS T_PP T_PS]，可利用VTI等效介质的Graebner反射系数精确方程计算获得：

S·r＝b；

其具体矩阵表达式如下：

其中，上标T和B代表反射界面上下层介质的弹性参数或中间变量；下标P和S参数表示与纵波和横波相关的弹性参数及相关中间变量；此外，p表示纵波水平慢度，s_*为垂直慢度；参数l、n、a、b、e均为数学中间变量；以上参数均可由刚度参数c_*和密度ρ计算求得；设定d_PP和d_PS为实际PP和PS角道集数据，则两种反射波的实际与合成数据残差为：

△d_PP＝d_PP-G_PP(m)；

△d_PS＝d_PS-G_PS(m)。

5.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：计算PP和PS数据正演过程对模型参数的正演偏导数的方法包括如下步骤：

PP和PS反射波正演过程G_PP(m)和G_PS(m)对模型参数m的偏导数如下：

其中，模型参数向量m包括：

m＝[c₁₁ c₁₃ c₃₃ c₅₅ ρ]；

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

其中，

和

为中间变量矩阵对不同模型参数的偏导数：

6.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：建立反演目标泛函，根据合成数据残差和正演偏导数，计算目标泛函对模型参数的泛函偏导数的方法包括如下步骤：

基于正则化理论，建立目标泛函：

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

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

7.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：根据泛函偏导数计算初始模型更新梯度，并更新初始模型的方法包括如下步骤：

针对线性反演策略，模型的更新迭代格式如下：

△_km＝m_k+1-m_k＝-β·g_k；

其中，k表示第k次反演迭代，β表示更新步长，g_k为第k次迭代的更新梯度，负梯度方向即为模型的更新方向；选用Guass-Newton优化算法，梯度计算式如下：

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

其中，

8.根据权利要求1所述的一种页岩VTI储层的各向异性参数多波联合反演方法，其特征在于：通过更新后的初始模型得出的反演结果换算各向异性参数的方法包括如下步骤：

反演模型组合为{c₁₁,c₁₃,c₃₃,c₅₅,ρ}，可利用数学关系换算两个Thomsen各向异性参数ε和δ；具体数学表达式如下：

Images