CN110542928A - 基于vti各向异性传播矩阵的地震响应模拟方法 - Google Patents

Abstract

本发明的一种基于VTI各向异性传播矩阵的地震响应模拟方法，包括以下步骤:根据测井数据通过岩石物理建模得到VTI介质弹性参数；根据非弹性薄互层反射系数的广义传播矩阵理论，计算相应反射波的频变反射系数R_f；根据频变反射系数R_f求得相应反射波的振幅谱U_f；对U_f做反傅里叶变换得到反射波波形u_t，本发明方法获得的地震合成记录与常规褶积模型获得的地震合成记录相比，提高了对波形的描述精度。

Description

基于VTI各向异性传播矩阵的地震响应模拟方法

技术领域

本发明涉及油气地球物理技术领域，尤其涉及一种基于VTI各向异性传播矩阵的地震响应模拟方法。

背景技术

地震资料研究表明，具有水平对称轴的裂隙诱导的VTI各向异性与具有沉积平面垂直对称轴的沉积或岩性所诱导的各向异性组合广泛存在。另外，储层中两种岩性的薄互层也经常出现。常规地震合成记录是由反射系数与子波褶积得到的，其中反射系数作为可以是由Zoeppritz方程或反射系数近似表达式计算的。常用的反射系数近似表达式局限于描述单一界面的假设，并且基于时间域的褶积计算方法不能够描述由于薄互层的反射地震波呈现复合波即干涉、调谐等现象。Aki和Richards于1980年提出了应用传播矩阵模拟地震波在层状介质中的响应特征，奠定了利用频率域信息描述复杂复合地震波场的基础。随着地震各向异性的测井分析、地震处理与解释的发展，需要能够同时模拟薄互层和各向异性特征的储层地震波的响应特征方法，提高合成地震记录精度。

地震合成记录与实际地震记录进行井震标定的过程中，根据测井数据计算的地震合成记录信息越全面，越接近真实结果，为地震解释及反演提供支撑的井震标定的精度越高。目前，常规的地震合成记录是建立在单一界面的假设理论之下，且很少考虑各向异性。本申请旨在解决或至少一定程度上缓解现有技术中的缺陷。

发明内容

为了解决上述技术问题，本发明提出一种提高地震响应模拟的精度的基于VTI各向异性传播矩阵的地震响应模拟方法。

本发明的一种基于VTI各向异性传播矩阵的地震响应模拟方法，包括以下步骤:

根据测井数据通过岩石物理建模得到VTI介质弹性参数；

根据非弹性薄互层反射系数的广义传播矩阵理论，计算相应反射波的频变反射系数R_f；

根据频变反射系数R_f求得相应反射波的振幅谱U_f；

对U_f做反傅里叶变换得到反射波波形u_t。

进一步地，所述岩石物理建模过程包括：

由HS上下界限平均理论计算泥页岩中非粘土类组分等效颗粒模量，得到由矿物组分计算的各向同性介质弹性模量；

由Backus平均理论计算页岩VTI各向异性固体基质的弹性参数；

应用多尺度裂缝理论考虑水平微裂缝的形状、孔缝系统连通性、流体类型和粘滞性等因素，得到VTI各向异性弹性模量；其中，水平微裂纹为本领域技术人员公知的技术，在此不做赘述。

计算VTI介质横波速度、各向异性参数。

进一步地，所述频变反射系数R_f包括薄互层的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T。

进一步地，

(1)对于P波入射，所述薄互层的反射、透射系数向量r＝[R_PP,R_PS,T_PP,

T_PS]^T由下式计算，

其中，矩阵A₁与A₂分别为与上、下层介质物性参数有关的传播矩阵；

B_α(α＝1,...,N)为具有N层结构的中间薄互层的传播矩阵；i_P为P波入射向量，

其中，

其中，i为虚数单位，ω为入射波频率，为薄互层总厚度，变量β、γ、W和Z的下标P、S分别对应准纵波、准横波，1、2分别对应上、下介质，P1、P2、S1、S2则分别表示上介质的准纵波、下介质的准纵波、上介质的准横波和下介质的准横波，

其中，β、γ、W和Z去掉下标后，简写各自表达式为，

W＝p₅₅(γs_x+βs_z)

Z＝βp₁₃s_x+γp₃₃s_z

其中，p.v.意为取复数的主值，P_ij是VTI介质的刚度，ρ是VTI介质密度。对于γ，符号“+”对应qP波(即准P波)，符号“－”对应qS波(即准S波)，s_x为水平波慢度，

E＝{[(p₃₃-p₅₅)cos²θ-(p₁₁-p₅₅)sin²θ]²+(p₁₃+p₅₅)²sin²2θ}^1/2

s_z为垂直波慢度，

其中，在sz的符号表达式，(+,-)表示向下传播qP波，(+,+)表示向下传播qS波，(-,-)表示向上传播qP波，(-,+)表示向上传播qS波，K1、K2、K3，是为了简化公式所设的中间参数，无实际意义；

Bα＝T(0)T-¹(hα)，

其中，

(2)所述P波入射向量i_P通过下式计算：

i_P＝iω[β_P1,γ_P1,-Z_P1,-W_P1]^T。

进一步地，所述振幅谱U_f通过下式计算得到：

U_f＝W_f×R_f

其中，W_f为频率域的地震子波。

进一步地，所述反射波波形u_t为

其中，f表示角频率，i为虚数单位，t为时间，由于各个频率下的反射、透射系数向量是根据VTI介质岩石物理反演得到的弹性参数计算得到的，因此此时得到的反射波形可以认为是基于VTI各向异性传播矩阵计算得到的。

与现有技术相比，本发明的基于VTI各向异性传播矩阵的地震响应模拟方法利用传播矩阵方法基于层状模型的特征，模拟地震响应时，不仅考虑与频散和衰减有关的地层岩性因此，而且考虑与薄层、薄互层有关的地层结构因素以及介质的各向异性，本发明获得的地震合成记录与常规褶积模型获得的地震合成记录相比，提高了对波形的描述精度，即提高了地震响应模拟的精度，也有助于提高对储集岩分布特征的认识。

上述技术特征可以各种适合的方式组合或由等效的技术特征来替代，只要能达到本发明的目的。

附图说明

在下文中将基于仅为非限定性的实施例并参考附图来对本发明进行更详细的描述。其中：

图1是本发明的基于VTI各向异性传播矩阵的地震响应模拟方法流程图；

图2是本发明实施例某井页岩储层段利用现有方法地震响应模拟结果；

图3是图2所示的某井页岩储层段考虑了VTI各向异性的地震响应模拟结果；

图4是利用本发明计算的PP波叠加数据随VTI介质的厚度变化结果。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

下面以某井页岩储层段基于VTI各向异性传播矩阵的地震响应模拟的具体实例来说明本发明的技术方案。

一种基于VTI各向异性传播矩阵的地震响应模拟方法，该方法包括以下步骤:

步骤1、根据测井数据通过岩石物理建模得到VTI介质弹性参数；首先，根据矿物组分、孔隙度、流体性质、有机质含量、纵波速度、横波速度等测井数据通过岩石物理建模得到VTI介质弹性参数，进而预测横波速度和各向异性参数等作为下一步计算弹性参数的准备条件。

步骤2、根据非弹性薄互层反射系数的广义传播矩阵理论，计算相应反射波的频变反射系数R_f；R_f是指各个频率下的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T。

步骤3、根据频变反射系数R_f求得相应反射波的振幅谱U_f，其中，

U_f＝W_f×R_f。

步骤4，对U_f做反傅里叶变换得到反射波波形u_t，

其中f表示角频率，i为虚数单位，t为时间。由于各个频率下的反射、透射系数向量是根据VTI介质岩石物理反演得到的弹性参数计算得到的，因此此时得到的反射波形可以认为是基于VTI各向异性传播矩阵计算得到的。

进一步地，步骤2中，岩石物理建模过程包括：

首先，由HS(Hashin-Shtrikman)上下界限平均理论计算泥页岩中石英、方解石、白云石以及干酪根等非粘土类组分等效颗粒模量，得到由矿物组分计算的各向同性介质弹性模量。

其次，由Backus平均理论计算页岩VTI各向异性固体基质的弹性参数。

然后，应用Chapman的多尺度裂缝理论考虑水平微裂缝的形状、孔缝系统连通性、流体类型和粘滞性等因素，得到VTI各向异性弹性模量。

最后，计算VTI介质横波速度、各向异性参数。

进一步地，对于P波入射，薄互层的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T由下式计算，

其中，矩阵A₁与A₂分别为与上、下层介质物性参数有关的传播矩阵；B_α(α＝1,...,N)为具有N层结构的中间薄互层的传播矩阵；i_P为P波入射向量，其中，

其中，i为虚数单位，ω为入射波频率，为薄互层总厚度，变量β、γ、W和Z的下标P、S分别对应准纵波、准横波，1、2分别对应上、下介质，P1、P2、S1、S2则分别表示上介质的准纵波、下介质的准纵波、上介质的准横波和下介质的准横波，

其中，β、γ、W和Z去掉下标后，

W＝p₅₅(γs_x+βs_z)

Z＝βp₁₃s_x+γp₃₃s_z

其中，p.v.意为取复数的主值。P_ij是VTI介质的刚度，ρ是VTI介质密度。对于γ，符号“+”对应qP波，符号“－”对应qS波，s_x为水平波慢度，s_z为垂直波慢度，

E＝{[(p₃₃-p₅₅)cos²θ-(p₁₁-p₅₅)sin²θ]²+(p₁₃+p₅₅)²sin²2θ}^1/2

其中，在sz的符号表达式，(+,-)表示向下传播qP波(+,+)表示向下传播qS波，(-,-)表示向上传播qP波，(-,+)表示向上传播qS波，K1、K2、K3，是为了简化公式所设的中间参数，无实际意义；

B_α＝T(0)T^-1(h_α)，

其中，

(2)所述P波入射向量i_P通过下式计算：

i_P＝iω[β_P1,γ_P1,-Z_P1,-W_P1]^T。

进一步地，由非弹性薄互层反射系数的广义传播矩阵理论，可计算各个频率下的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T，即相应反射波的频变反射系数R_f，振幅谱U_f通过频变反射系数与频率域的地震子波W_f相乘计算得到，即：

U_f＝W_f×R_f。

进一步地，通过对U_f做反傅里叶变换则可以获得时间域的反射波波形u_t，获得VTI介质薄互层的高精度地震响应特征。

其模拟结果如图2和图3所示，图2为该井页岩储层段利用现有方法地震响应模拟结果，图3为该井页岩储层段考虑了VTI各向异性的地震响应模拟结果。

如图2所示，图中左边三列分别是纵波速度、横波速度和各向异性参数。右边是各向同性反射率法计算的地震响应特征。如图3所示，左边三列分别是纵波速度、横波速度和各向异性参数。右边是VTI各向异性传播矩阵的地震响应特征。从图3中可以观察到测井数据中地层高、低波阻抗的变化与地震反射波峰、波谷具有良好的对应关系。井震标定结果能够为地震解释与反演提供依据。另外，从图2和图3的对比发现，考虑了VTI各向异性的结果，可以得到更精确的地震AVO响应特征，其主要表现在远偏移距振幅降低和相位的变化。

另外，波的干涉作用、以及页岩地层内部波的各向异性传播使反射波形和振幅都趋于复杂。

本发明的基于VTI各向异性传播矩阵的地震响应模拟方法基于岩石物理建模及反演得到的VTI介质弹性参数，应用各向异性传播矩阵理论计算具有各向异性的油气层地震AVO响应。反射体模型为具有有限厚度、内部呈各向异性的地质体，此时的反射波是如下动力学信息的综合响应，包括界面处波阻抗的差异、地震波的在反射体内的各向异性扩散、以及反射波的调谐与干涉效应。利用传播矩阵方法基于层状模型的特征，模拟地震响应时，不仅考虑与频散和衰减有关的地层岩性因此，而且考虑与薄层、薄互层有关的地层结构因素以及介质的各向异性。因此，利用本发明的方法获得的地震合成记录与常规褶积模型获得的地震合成记录相比，提高了对波形的描述精度，即提高了地震响应模拟的精度。如图4所示，显示了利用本发明计算的PP波叠加数据随VTI介质的厚度从10m增加到70m时的变化情况，从图中可以看出PP波波形随VTI介质厚度的变化而明显变化。

虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的结果下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。

Claims (6)

1.一种基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，该方法包括以下步骤:

根据测井数据通过岩石物理建模得到VTI介质弹性参数；

根据非弹性薄互层反射系数的广义传播矩阵理论，计算相应反射波的频变反射系数R_f；

根据频变反射系数R_f求得相应反射波的振幅谱U_f；

对U_f做反傅里叶变换得到反射波波形u_t。

2.根据权利要求1所述的基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，所述岩石物理建模过程包括：

由HS上下界限平均理论计算泥页岩中非粘土类组分等效颗粒模量，得到由矿物组分计算的各向同性介质弹性模量；

由Backus平均理论计算页岩VTI各向异性固体基质的弹性参数；

应用多尺度裂缝理论考虑水平微裂缝的各因素，得到VTI各向异性弹性模量；

计算VTI介质横波速度、各向异性参数。

3.根据权利要求1所述的基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，所述频变反射系数R_f包括薄互层的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T。

4.根据权利要求3所述的基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，

(1)对于P波入射，所述薄互层的反射、透射系数向量r＝[R_PP,R_PS,T_PP,T_PS]^T由下式计算，

其中，矩阵A₁与A₂分别为与上、下层介质物性参数有关的传播矩阵；B_α(α＝1,...,N)为具有N层结构的中间薄互层的传播矩阵；i_P为P波入射向量，其中，

其中，i为虚数单位，ω为入射波频率，为薄互层总厚度，变量β、γ、W和Z的下标P1、P2、S1、S2分别表示上介质的准纵波、下介质的准纵波、上介质的准横波和下介质的准横波，

其中，β、γ、W和Z去掉下标后，

W＝p₅₅(γs_x+βs_z)

Z＝βp₁₃s_x+γp₃₃s_z

其中，p.v.意为取复数的主值，P_ij是VTI介质的刚度，ρ是VTI介质的密度。对于γ，符号“+”对应qP波，符号“－”对应qS波，s_x为水平波慢度，s_z为垂直波慢度，

E＝{[(p₃₃-p₅₅)cos²θ-(p₁₁-p₅₅)sin²θ]²+(p₁₃+p₅₅)²sin²2θ}^1/2

其中，s_z的符号表达式，+、－表示向下传播qP波，+、+表示向下传播qS波，－、－表示向上传播qP波，－、+表示向上传播qS波，K₁、K₂、K₃，是为了简化公式所设的中间参数，无实际意义；

Bα＝T(0)T^-1(hα)，

其中，

其中，s_p和s_s分别表示纵波和横波的波慢度，

(2)所述P波入射向量i_P通过下式计算：

i_P＝iω[β_P1,γ_P1,-Z_P1,-W_P1]^T。

5.根据权利要求4所述的基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，所述振幅谱U_f通过下式计算得到：

U_f＝W_f×R_f

其中，W_f为频率域的地震子波。

6.根据权利要求5所述的基于VTI各向异性传播矩阵的地震响应模拟方法，其特征在于，所述反射波波形u_t为

其中f表示角频率，i为虚数单位，t为时间。