CN115392090B - 基于三维数字岩心的地震波频散和衰减特征的预测方法 - Google Patents

Abstract

本发明公开基于三维数字岩心的地震波频散和衰减特征的预测方法，涉及油气勘探开发技术领域，包括以下步骤：基于数字岩心技术获得流体分布；获得应力张量与岩石骨架位移场及流体压力的关系；按照研究目标尺寸构建均匀背景模型；使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场。本发明可实现以数字岩心技术表征三维流体分布，绕开传统基于流体分布建模过程，解决了三维模型中大纵横比、小目标体无法建模的困难，精准预测地震波频散和衰减特征。

Description

基于三维数字岩心的地震波频散和衰减特征的预测方法

技术领域

本发明涉及油气勘探开发技术领域，具体为基于三维数字岩心的地震波频散和衰减特征的预测方法。

背景技术

流体分布是描述地下流体特征的重要物性参数，其与地震波频散和衰减的定量关系是地球物理流体反演的重要约束条件。依据流体分布预测纵波模量的基本方法是以Biot框架为理论背景，推导出两种流体分布下的频散和衰减的解析解。如假设流体分布为周期层状和球状两种规则水气分布模式，White及Dutta推导了流体分布差异带来的纵波频散特征。然而，层状和球状水气分布模式与实际的复杂水气分布情况相差很远，更为接近实际情况的方法是假设流体分布具有随机性，使用随机介质模型中统计平滑概念表征复杂的流体分布，同样依据Biot框架推导解析解，实现对速度频散和衰减特征的预测。对于更加复杂的实际流体分布，无法推导解析解，更为广泛的解法是利用有限元方法求解准静态Biot孔弹性方程，实现对孔隙岩石中非均质性(流体分布、裂缝、非均匀孔隙)诱发频散和衰减的预测。该方法广泛的适用性使得输入真正流体分布预测频散和衰减存在可能。随着时间推移，知识积累，逐渐形成了以CT图像所得流体分布作为输入，预测纵波模量频散和衰减的理念。

然而，现有方法及技术在实现上述理念时存在巨大挑战，尤其对于三维岩石样品而言，依据流体分布建模预测其频散和衰减几乎是不可能实现的，原因是：在三维流体分布中，流体分布特征表现出气泡小、分散广及纵横比差异大的特点，从而无法依据CT扫描结果建立三维几何模型并实现网格划分，使得三维地震波频散和衰减的预测存在困难。

发明内容

本发明针对现有技术存在的困境，提供了一种基于三维数字岩心的地震波频散和衰减特征的预测方法。该方法基于数值岩心技术评估三维流体分布，通过间接使用该信息，无需依据流体分布建立几何模型，实现精准预测地震波频散和衰减，解决了常规技术中面对小地质体、大纵横比流体分布特征时无法建立几何模型、划分网格的问题。

为了解决上述技术问题，本发明提供了如下的技术方案：

基于三维数字岩心的地震波频散和衰减特征的预测方法，包括以下步骤：

S1：基于数字岩心技术获得流体分布；

S2：获得应力张量与岩石骨架位移场及流体压力的关系；

S3：按照研究目标尺寸构建均匀背景模型；

S4：使用有限元求解均匀背景模型，获得固体位移场和流体压力场。

作为本发明的进一步技术方案为，所述基于数字岩心技术获得流体分布，具体为：

其中，CT_Water是样品完全饱和水时的灰度值，CT_Air是样品饱和空气时的灰度值，其中，CT_water+air是部分饱和岩石的灰度值，θ是规则化系数，其中，0＜θ≤1。

作为本发明的进一步技术方案为，所述基于数字岩心技术获得流体分布，还包括：

根据流体的体积模量、粘度、密度以及流体饱和度定义关键孔弹性参数，具体为：

K_f＝S_AirK_Air+K_w(1-S_Air)， (2-1)

ρ_f＝ρ_AirS_Air+ρ_w(1-S_Air)， (2-3)

其中，K_Air和K_w分别为空气和水的体积模量，ρ_Air和ρ_w分别为空气和水的密度，η_Air和η_w分别为空气和水的粘度，S_Air为空气饱和度。

作为本发明的进一步技术方案为，所述获得应力张量与岩石骨架位移场及流体压力的关系，具体包括：

多孔材料中固体和流体应满足质量守恒方程(3)和达西定律(4)，其频率域形式为：

其中，

是哈密顿算子，ω是角频率，α是Biot-Willis系数，

ρ_b＝(1-φ)ρ_s+φρ_f， (5)

其中，ρ_f是流体密度，ρ_s是固体骨架密度，φ是孔隙度，

其中，S是孔隙的迂曲因子，η是流体的粘度，κ是渗透率，

是虚数单位，岩石骨架的位移为：

流体压力为：

P_f＝(P_i)， (8)

岩石骨架的应变张量为：

其中i,j＝1,2,3，是欧几里得空间维数；

应力张量σ与岩石骨架位移场及流体压力的关系为：

其中，

K_d是样品的排水体积模量，G_d是样品排水剪切模量，K_f是流体的体积模量，K_g是岩石颗粒的体积模量，δ_ij函数为：

作为本发明的进一步技术方案为，所述按照研究目标尺寸构建均匀背景模型，具体包括：对研究目标尺寸进行常规四面体网格剖分，以及依据步骤S1所获结果寻找均匀背景模型中异常边界，进行网格微调：

其中Mat是与空间变量x，y，z，网格大小变量d有关的位置函数，用于定位异常边界位置；根据Mat值的范围，确定粗网格和细网格的位置FP。

作为本发明的进一步技术方案为，所述按照研究目标尺寸构建均匀背景模型，还包括：定义边界条件：

①在不排水条件下，样品所有边界流量应为0，其频率域形为：

其中n是边界外法向向量，

是虚数单位；

②在笛卡尔坐标系下，根据纵波模量的定义，研究目标区域侧面边界横向应变为0，即：

ε_ij＝0， (15)

其中，i,j＝1,2是欧几里得空间维数，对应x、y方向；

样品下边界为固定边界，即：

u_i＝0， (16)

其中，i＝1,2,3是欧几里得空间维数，对应x、y、z方向；目标区域上边界加载周期震荡应力S0；

σ₃₃＝S0。 (17)

作为本发明的进一步技术方案为，所述使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场，具体为：

纵波模量定义为M为：

其中，σ₃₃为样品垂向应力，ε₃₃为样品垂向应变；则纵波模量衰减为QM：

本发明的有益效果是：

本发明可实现以数字岩心技术表征三维流体分布，并预测地震波频散和衰减特征，其优点是绕开传统方法中基于流体分布的建模过程，解决三维模型中大纵横比、小目标体无法建模的困难，精准预测地震波频散和衰减特征，使用本发明计算时，直接将流体分布信息耦合到物理场中，避开使用复杂几何模型计算时存在的劣质网格，大幅度提高计算速度。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。

图1为本发明提出的基于三维数字岩心的地震波频散和衰减特征的预测方法流程图；

图2为常规计算方法与本发明方法对比图；

图3为本发明实施例提供的汽水分布图；

图4a为本发明实施例提供的流体模型图；

图4b为本发明实施例提供的耦合弹性参数模型图；

图4c为本发明实施例提供的密度模型图；

图4d为本发明实施例提供的流体粘度模型图；

图5为本发明实施例边界加载情况图；

图6为本发明实施例网格划分结果图；

图7a为本发明实施例在频率20000Hz时的流体压力分布图；

图7b为本发明实施例在频率20000Hz时的应变图；

图8a为本发明实施例提供的纵模模量测量结果与本发明预测结果对比图。

图8b为本发明实施例提供的纵波衰减测量结果与本发明预测结果对比图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

参见图1，基于三维数字岩心的地震波频散和衰减特征的预测方法，包括以下步骤：

步骤S1，基于数字岩心技术获得流体分布；

步骤S2，获得应力张量与岩石骨架位移场及流体压力的关系；

步骤S3，按照研究目标尺寸构建均匀背景模型；

步骤S4，使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场。

其中步骤S1中，所述基于数字岩心技术获得流体分布；具体为：

其中，CT_Water是样品完全饱和水时的灰度值，CT_Air是样品饱和空气(干燥样品)的灰度值，其中，CT_water+air是部分饱和岩石的灰度值(即同时包括水和空气)，θ是规则化系数，其中，0＜θ≤1，保证饱和度数值处于0和1之间。

在步骤S1中，基于数字岩心技术获得流体分布：还包括：

根据流体的体积模量、粘度以及密度与流体饱和度定义关键孔弹性参数：流体的体积模量、粘度以及密度与流体饱和度密切相关，包含在方程组(2-1)～(2-4)中:

K_f＝S_AirK_Air+K_w(1-S_Air)， (2-1)

ρ_f＝ρ_AirS_Air+ρ_w(1-S_Air)， (2-3)

其中，K_Air和K_w分别为空气和水的体积模量，ρ_Air和ρ_w分别为空气和水的密度，η_Air和η_w分别为空气和水的粘度，S_Air为空气饱和度。

在步骤S2中，获得应力张量与岩石骨架位移场及流体压力的关系；具体包括：

多孔材料中固体和流体应满足质量守恒方程(3)和达西定律(4)，其频率域形式为：

其中，

是哈密顿算子，ω是角频率，α是Biot-Willis系数，

ρ_b＝(1-φ)ρ_s+φρ_f， (5)

其中，ρ_f是流体密度，ρ_s是固体骨架密度，φ是孔隙度，

其中，S是孔隙的迂曲因子，η是流体的粘度，κ是渗透率，

是虚数单位，岩石骨架的位移张量为：

流体压力为：

P_f＝(P_i)， (8)

岩石骨架的应变张量为：

其中i,j＝1,2,3，是欧几里得空间维数；

应力张量σ与岩石骨架位移场及流体压力的关系为：

其中，

K_d是样品的排水体积模量，G_d是样品排水剪切模量，K_f是流体的体积模量，K_g是岩石颗粒的体积模量，δ_ij函数为：/>

在步骤S3中，所述按照研究目标尺寸构建均匀背景模型，具体包括：对研究目标尺寸进行常规四面体网格剖分，以及依据步骤S1所获结果寻找均匀背景模型中异常边界，进行网格微调：

其中Mat是与空间变量x，y，z，网格大小变量d有关的位置函数。根据Mat值得范围，确定粗网格和细网格的位置FP。

在步骤S3中，所述按照研究目标尺寸构建均匀背景模型，还包括：定义边界条件：

①在不排水条件下，样品所有边界流量应为0，其频率域形为：

其中n是边界外法向向量，

是虚数单位；

②在笛卡尔坐标系下，根据纵波模量的定义，研究目标区域侧面边界横向应变为0，即：

ε_ij＝0， (15)

其中，i,j＝1,2是欧几里得空间维数，对应x、y方向。样品下边界为固定边界，即：

u_i＝0， (16)

其中，i＝1,2,3是欧几里得空间维数，对应x、y、z方向；目标区域上边界加载周期震荡应力S0，

σ₃₃＝S0。 (17)

在步骤S4中，所述使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场，具体为：

纵波模量定义为M为：

其中，σ₃₃为样品垂向(z方向)应力，ε₃₃为样品垂向应变；则纵波模量衰减为Q_M：

本发明所述方法可实现以数字岩心技术表征三维流体分布，并预测地震波频散和衰减特征。其优点是绕开传统方法中基于流体分布的建模过程，从而在保证预测精度的同时，解决了三维模型中大纵横比、小目标体无法建模的困难。

为说明发明效果，以附图2(a)表示三维圆柱形模型为例，该模型内部圆球为空气，外部包裹水。

传统方法是根据附图2(a)建立几何模型，划分网格，计算结果获得的纵向位移场，如图2(c)所示，计算过程总耗时为30分32秒。

使用本发明计算时，不需要建立几何模型，直接将与流体分布信息耦合到孔弹方程中，可以看到，二者计算结果一致，但本发明所使用的计算时间为20分钟10秒，提高五倍之多，且避开了建立几何模型细节的过程。

选取某碳酸盐为实施例，其具体参数如表1所示。利用排水方法饱和该岩石，使其含水饱和度为90％，含气饱和度为10％。为使本发明的目的、技术方案和优点表述更加清楚，依据本发明技术方案，对纵波模量的频散和衰减特征进行预测。将预测结果与测量结果(如图8a和图8b正方形点所示)对比，以展示本发明所得预测结果的可靠性。

表1某碳酸岩、水和空气的特性

具体步骤为：

根据CT扫描图像，利用数字岩心技术，计算汽水分布特征，如图3所示。可以看到，该流体分布中存在很多不规则的气泡，使用常规方法无法建模；

利用上述的流体分布计算各个弹性参数，计算结果如图4所示：其中，图4a为流体模型；图4b为耦合弹性参数；图4c为密度；图4d为流体粘度；

建立均匀背景模型，如图5所示；

划分网格，如图6所示，充满小气泡的位置网格得到加密；

定义边界条件，如图7a和7b所示，圆柱形顶部记载应力S0，底部为固定约束，侧面的偏应变为0；

求解方程，得到流体压力(图7a)，骨架垂向应变(图7b)，据此计算纵向纵波模量和对应的频散衰减关系，如图8a和图8b中的黑色曲线，分别为纵波模量和纵波衰减，可以看到预测结果与实测结果匹配结果良好。证明了本发明技术的有效性。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，包括以下步骤：

S1：基于数字岩心技术获得流体分布；

S2：获得应力张量与岩石骨架位移场及流体压力的关系；

S3：按照研究目标尺寸构建均匀背景模型；

S4：使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场；

所述获得应力张量与岩石骨架位移场及流体压力的关系，具体包括：

多孔材料中固体和流体应满足质量守恒方程(3)和达西定律(4)，其频率域形式为：

其中，

是哈密顿算子，ω是角频率，α是Biot-Willis系数，

ρ_b＝(1-φ)ρ_s+φρ_f， (5)

其中，ρ_f是流体密度，ρ_s是固体骨架密度，φ是孔隙度，

其中，S是孔隙的迂曲因子，η是流体的粘度，κ是渗透率，

是虚数单位，岩石骨架的位移为：

流体压力为：

P_f＝(P_i)； (8)

岩石骨架的应变张量为：

其中i,j＝1,2,3，是欧几里得空间维数，

应力张量σ与岩石骨架位移场及流体压力的关系为：

其中，

K_d是样品的排水体积模量，G_d是样品排水剪切模量，K_f是流体的体积模量，K_g是岩石颗粒的体积模量，δ_ij函数为：

2.根据权利要求1所述的基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，所述基于数字岩心技术获得流体分布，具体为：

/>

其中，CT_Water是样品完全饱和水时的灰度值，CT_Air是样品饱和空气的灰度值，其中，CT_water+air是部分饱和岩石的灰度值，θ是规则化系数，其中，0＜θ≤1。

3.根据权利要求1所述的基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，所述基于数字岩心技术获得流体分布，还包括：

根据流体的体积模量、粘度、密度以及流体饱和度定义关键孔弹性参数，具体为：

K_f＝S_AirK_Air+K_w(1-S_Air)， (2-1)

ρ_f＝ρ_AirS_Air+ρ_w(1-S_Air)， (2-3)

其中，K_Air和K_w分别为空气和水的体积模量，ρ_Air和ρ_w分别为空气和水的密度，η_Air和η_w分别为空气和水的粘度，S_Air为空气饱和度。

4.根据权利要求1所述的基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，所述按照研究目标尺寸构建均匀背景模型，具体包括：对研究目标尺寸进行常规四面体网格剖分，以及依据步骤S1所获结果寻找均匀背景模型中异常边界，进行网格微调：

其中Mat是与空间变量x，y，z，网格大小变量d有关的位置函数，用于定位异常边界位置；根据Mat值的范围，确定粗网格和细网格的位置FP。

5.根据权利要求1所述的基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，所述按照研究目标尺寸构建均匀背景模型，还包括：定义边界条件：

①在不排水条件下，样品所有边界流量应为0，其频率域形为：

其中n是边界外法向向量，

是虚数单位；

②在笛卡尔坐标系下，根据纵波模量的定义，研究目标区域侧面边界横向应变为0，即：

ε_ij＝0， (15)

其中，i,j＝1,2是欧几里得空间维数，对应x、y方向；

样品下边界为固定边界，即：

其中，i＝1,2,3是欧几里得空间维数，对应x、y、z方向；目标区域上边界加载周期震荡应力S0；

σ₃₃＝S0。 (17)

6.根据权利要求1所述的基于三维数字岩心的地震波频散和衰减特征的预测方法，其特征在于，所述使用有限元求解均匀背景模型，获得应变固体位移场和流体压力场，具体为：

纵波模量定义为M为：

其中，σ₃₃为样品垂向应力，ε₃₃为样品垂向应变；则纵波模量衰减为Q_M：

/>

Images