CN115308798B - 一种高速低内存消耗的预测储层岩石弹性波速度的方法 - Google Patents

Abstract

本发明公开了一种高速低内存消耗的预测储层岩石弹性波速度的方法，该方法不仅能够精确的预测实验结果，且具有内存消耗小，计算时间短的特点。其成功的解决了无法使用传统孔弹性方程预测实验室三维实验样品弹性波速度的问题，在岩石力学、岩石物理学，油气定量地震学等解释等领域具有重要的研究意义。

Description

一种高速低内存消耗的预测储层岩石弹性波速度的方法

技术领域

本发明涉及油气勘探开发的技术领域，具体涉及一种高速低内存消耗的预测储层岩石弹性波速度的方法。

背景技术

地震波在岩石储层中的传播速度是地震勘探的重要基础。地下岩石由孔隙、岩石骨架以及孔隙中的流体构成，称为孔弹性介质。当地震波通过孔弹性介质时，诱发孔隙中的流体和岩石骨架发生相互作用，岩石表现出不同的孔弹性参数。实验是获得岩石弹性参数最直接的方法，也是定量地震学中重要的约束条件。然而，实验成本昂贵且需要大量的技术成本，作为替代方法，建立数值模型预测饱和流体岩石的弹性波速度不仅是解释实验室测量数据的重要理论依据，更是一种预测孔弹性介质声波速度的重要手段。

对于地下储层，流体和固体的相互作用规律满足动量守恒定律，与此同时，地层内流体流动受到流体压力和岩石应力影响，其流动过程遵守达西定律。上述两组定律结合即为Biot孔弹性方程，其描述了弹性波在地下岩石中传播时的基本规律。然而上诉规律描述的是储层的宏观规律，无法直接用于解释实验室测量数据，为此Rubino在2016年于JGR-Solid earth期刊上发表了题为《Numerical upscaling in 2-D heterogeneousporoelastic rocks:Anisotropic attenuation and dispersion of seismic waves》的文章，通过对二维模型设置纵波边界条件，实现以孔弹性理论解释介观裂缝诱发地震波的频散和衰减的机理。虽然该方法提出了使用基础定律预测介观介质弹性特征的可能，但并不能用于解释实验室测量结果，究其原因：首先，实验室岩石样品为圆柱形，且测量过程不满足所谓的纵波边界条件，因此无法直接使用上述方法预测实验室岩石样品的弹性参数；其次，上述方法在求解过程中使用了孔弹性方程的简化版本，使得其在理论上存在瑕疵，限制了模型的计算精度；此外，上述方法只能应用于二维模型，无法建立三维模型模拟实际样品，其本质原因是：有限元方法求解时通常使用牛顿迭代法和LU矩阵分解法。牛顿迭代法计算需要内存相对较低，计算效率相对高，但需要雅可比矩阵为正定矩阵；LU矩阵分解法则计算精度高，可直接求解，但涉及到矩阵分解，因此需要大量的内存，且计算计算时间长。对于孔弹性方程而言，其刚度矩阵受多种因素的影响，如流体分布，孔隙分布特征等，无法确定雅可比矩阵的正定性，保证牛顿迭代法顺利求解。因此，LU矩阵分解法成为求解孔弹性方程的必要选择。对于三维模型而言，使用LU矩阵分解求解需要大量的计算时间和内存，随着离散化自由度的提升，内存增长速度和计算时间陡然上升。如对于一个直径40，长度80里面的样品，如考虑有效描述样品中流体分布特征或其孔隙造成的非均质性，计算内存甚至需要2-3TB左右，这远超现有计算机的配置极限，因此建立3D模型描述实际样品存在挑战。综上所述，目前缺少一种建立三维模型，并能正确求解、预测实验室弹性波速度的方法。

发明内容

鉴于现有技术的不足，本发明旨在于提供一种高速低内存消耗的预测储层岩石弹性波速度的方法，本发明不仅能够精确的预测实验结果，且具有内存消耗小，计算时间短的特点。其成功的解决了无法使用传统孔弹性方程预测实验室三维实验样品弹性波速度的问题，在岩石力学、岩石物理学，油气定量地震学等解释等领域具有重要的研究意义。

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

一种高速低内存消耗的预测储层岩石弹性波速度的方法，所述方法包括：

S1获取建立模型所需要的参数；其中，需要岩石的物性参数作为输入，可通过实验室测量获取，

S2依据岩石样品尺寸，按照岩石样品R1上端黏贴铝块A1、下端黏贴铝块A2的方法构建模型；

S3设置铝块样品支配方程及边界条件；

S4设置岩石样品的支配方程及边界条件；

S5网格剖分；其中，应变片S1和S2进行细化网格剖分；上部铝块A1和岩石样品R1进行常规四面体网格剖分；根据流体的分布特征及样品的非均质性特征，进行精细化网格的剖分；下部铝块A2按照样品R1下表面网格进行延拓，进行三维空间的剖分；

S6依据网格离散方程组，并利用混合求解方法求解；

S7计算纵横波速度。

需要说明的是，所述步骤S1中，包括：

S1.1参考铝制样，长度为40mm，直径为40mm，密度为：

ρ_al＝2700Kg*m^-3；

S1.2使用游标卡尺测量岩石样品的直径DR，长度LenR，利用天平测量干燥岩石质量MR，计算干燥岩石密度：

S1.3依据SY/T 6385-2016《覆压下岩石孔隙度和渗透率测定方法》标准，测定岩石样品的孔隙度φ和渗透率κ；

S1.4依据分析标准SY/T5163-2018，对岩石进行矿物分析；通过矿物分析和VRH公式，计算岩石颗粒体积模量K_g；

S1.5依据SY/T7410.1-2018分析标准对岩石孔隙结构和流体分布进行扫；描和重建，进而获得岩石的孔隙结构和流体分布特征；

S1.6依据流体分布特征计算流体的体积模量K_f，粘度η；

S1.7依据行波法，测得样品的干燥体积模量K_d和剪切模量μ。

3、根据权利要求1所述的速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S2中，铝块A1的中心和岩石样品R1的中心添加沿着纵向和横向分别添加两个线段(S1和S2)，用于表示测量轴向应变和径向应变的应变片。

需要说明的是，所述步骤S3中，包括：

铝块A1和A2是各向同性均匀材料，其支配方程为应力-应变方程，为提高计算效率，采用频率域形式:

其中ρ_al是铝制材料的密度，u^s＝[U，V，W]是铝制材料的位移场，σ是应力场张量。对于铝块A1而言，其上表面加载周期震荡震源(频率域为常数)，即轴向加载应力σ_zz＝0.1[MPa]，对于铝块A2而言，其下表面为固定位置，其位移张量u^s＝0。

需要说明的是，所述步骤S4中，包括：

岩石样品中流固耦合应该符合孔弹性方程：

其中

是哈密顿算子，ω是角频率；ρ_f是流体密度，ρ_s是固体骨架密度，ρ_b＝ρ_s+φρ_f，

φ是孔隙度，S是孔隙的曲折因子，η是流体的粘度，κ是渗透率，

是虚数单；P_f是流体压力，u^s＝[U，V，W]是固体骨架位移，σ是应力张量，与流体压力P_f和固体骨架位移u^s，干燥岩石体积模量K_d，干燥岩石剪切模量μ有关。其中α是Biot-Willis系数，M为耦合弹性模量：

岩石样品R1上下表面被铝块A1、A2封闭，侧面被胶套封闭，因此岩石各边界处于不排水条件，所有流量应为0，具体为：

其中n是R1表面外法向向量，

是虚数单位。

需要说明的是，所述步骤S5中，包括：

利用网格划分结果，通过离散化转为线性方程组：

其中

采用混合方法，即牛顿迭代法和LU矩阵分解法混合的方式求解线性方程组，具体步骤为：

A设置初始条件X^i＝0＝0；

B计算L₁和雅可比矩阵

则

C基于

计算L₂和雅可比矩阵

则

D基于

计算L₃和雅可比矩阵

则

E基于

利用LU矩阵分解法求解

F如果|Xⁱ⁺¹-Xⁱ|＜∈，则停止计算，否则回到步骤A。∈是相对误差。

需要说明的是，所述步骤S7中，包括：

依据应变计算杨氏模量E和泊松比ν：

其中杨氏模量E_al为铝块杨氏模量72GPa，

为S1所在位置轴向应变的平均值，

为S2所在位置轴向应变的平均值，

为S3位置径向应变的平均值。

纵波速度计算结果为：

横波速度计算结果为：

本发明有益效果在于：

1、在建模理论方面：传统建模方法使用纵波边界条件，与实验室测量样品的实际边界条件不同，无法正确解释实验结果。本发明从基础孔弹性力学理论出发，依据实验室测量样品实际情况设置边界条件，因此可正确解释实验室测量结果，并做出预测。

2、在适用范围方面：本发明适用于实验室测试样品尺度，通常在厘米量级。与传统方法相比，本发明考虑了尺度效应。

3、在计算效率方面：传统牛顿迭代方法无法求解该模型，而LU矩阵分解法求解则需要耗费大量的时间和内存。本发明提出的混合求解方法需要低内存和少量的时间，具有高速性和低内存消耗性特征，具体效果可参照表1。

4、在计算精度方面：传统方法因为尺度和边界条件不同，预测结果与实验结果存在很大差异。本发明基于实际情况，采用完整的孔弹性方程，选用合适的边界条件和求解方法，获得的预测结果于实验结果一致(图4)，证明了该方法的高精度性。

附图说明

图1为本发明的技术方案流程图；

图2为本发明的实施例建模示意图；

图3为本发明的实施例模型网格划分图；

图4为本发明的实施例模型求解过程各物理场收敛速度图；

图5为本发明的实施例预测结果与实际结果对比图。

具体实施方式

下将结合附图对本发明作进一步的描述，需要说明的是，本实施例以本技术方案为前提，给出了详细的实施方式和具体的操作过程，但本发明的保护范围并不限于本实施例。

如图1所示，本发明为一种高速低内存消耗的预测储层岩石弹性波速度的方法，所述方法包括：

S1获取建立模型所需要的参数；其中，需要岩石的物性参数作为输入，可通过实验室测量获取；

S2依据岩石样品尺寸，按照岩石样品R1上端黏贴铝块A1、下端黏贴铝块A2的方法构建模型；

S3设置铝块样品支配方程及边界条件；

S4设置岩石样品的支配方程及边界条件；

S5网格剖分；其中，应变片S1和S2进行细化网格剖分；上部铝块A1和岩石样品R1进行常规四面体网格剖分；根据流体的分布特征及样品的非均质性特征，进行精细化网格的剖分；下部铝块A2按照样品R1下表面网格进行延拓，进行三维空间的剖分；

S6依据网格离散方程组，并利用混合求解方法求解；

S7计算纵横波速度。

需要说明的是，所述步骤S1中，包括：

S1.1参考铝制样，长度为40mm，直径为40mm，密度为：

ρ_al＝2700Kg*m^-3；

S1.2使用游标卡尺测量岩石样品的直径DR，长度LenR，利用天平测量干燥岩石质量MR，计算干燥岩石密度：

S1.3依据SY/T 6385-2016《覆压下岩石孔隙度和渗透率测定方法》标准，测定岩石样品的孔隙度φ和渗透率κ；

S1.4依据分析标准SY/T5163-2018，对岩石进行矿物分析；通过矿物分析和VRH公式，计算岩石颗粒体积模量K_g；

S1.5依据SY/T7410.1-2018分析标准对岩石孔隙结构和流体分布进行扫；描和重建，进而获得岩石的孔隙结构和流体分布特征；

S1.6依据流体分布特征计算流体的体积模量K_f，粘度η；

S1.7依据行波法，测得样品的干燥体积模量K_d和剪切模量μ。

3、根据权利要求1所述的速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S2中，铝块A1的中心和岩石样品R1的中心添加沿着纵向和横向分别添加两个线段(S1和S2)，用于表示测量轴向应变和径向应变的应变片。

需要说明的是，所述步骤S3中，包括：

铝块A1和A2是各向同性均匀材料，其支配方程为应力-应变方程，为提高计算效率，采用频率域形式:

其中ρ_al是铝制材料的密度，u^s＝[U，V，W]是铝制材料的位移场，σ是应力场张量。对于铝块A1而言，其上表面加载周期震荡震源(频率域为常数)，即轴向加载应力σ_zz＝0.1[MPa]，对于铝块A2而言，其下表面为固定位置，其位移张量u^s＝0。

需要说明的是，所述步骤S4中，包括：

岩石样品中流固耦合应该符合孔弹性方程：

其中

是哈密顿算子，ω是角频率；ρ_f是流体密度，ρ_s是固体骨架密度，ρ_b＝ρ_s+φρ_f，

φ是孔隙度，S是孔隙的曲折因子，η是流体的粘度，κ是渗透率，

是虚数单；P_f是流体压力，u^s＝[U，V，W]是固体骨架位移，σ是应力张量，与流体压力P_f和固体骨架位移u^s，干燥岩石体积模量K_d，干燥岩石剪切模量μ有关。其中α是Biot-Willis系数，M为耦合弹性模量：

岩石样品R1上下表面被铝块A1、A2封闭，侧面被胶套封闭，因此岩石各边界处于不排水条件，所有流量应为0，具体为：

其中n是R1表面外法向向量，

是虚数单位。

需要说明的是，所述步骤S5中，包括：

利用网格划分结果，通过离散化转为线性方程组：

其中

采用混合方法，即牛顿迭代法和LU矩阵分解法混合的方式求解线性方程组，具体步骤为：

A设置初始条件X^i＝0＝0；

B计算L₁和雅可比矩阵

则

C基于

计算L₂和雅可比矩阵

则

D基于

计算L₃和雅可比矩阵

则

E基于

利用LU矩阵分解法求解

F如果|Xⁱ⁺¹-Xⁱ|＜∈，则停止计算，否则回到步骤A。∈是相对误差。

需要说明的是，所述步骤S7中，包括：

依据应变计算杨氏模量E和泊松比v：

其中杨氏模量E_al为铝块杨氏模量72GPa，

为S1所在位置轴向应变的平均值，

为S2所在位置轴向应变的平均值，

为S3位置径向应变的平均值。

纵波速度计算结果为：

横波速度计算结果为：

实施例

为使本发明的目的、技术方案和优点表述更加清楚,依据本发明技术方案流程，以实际岩石样品为实施例，预测其在部分饱和情况下纵横波速度。具体步骤为：

①依据技术方案步骤(1)所示，获得建立模型所需参数；经过测量样品尺寸直径为40mm，高为38mm，质量为113.05克(g)，密度为2367.4Kg/m³。岩石孔隙度为10.8％，渗透率为0.02毫达西(mD)，岩石颗粒体积模量K_g＝78GPa。样品饱和采用气、水混合，样品内的气水分布通过CT扫描获得，其中水的体积模量是2.25GPa，气体的体积模量是1Bar，水的密度是1000Kg/m³，空气的密度是140Kg/m³，水的粘度是0.001Pa·s，气体的粘度是1E-5Pa·s。排水体积模量为K_d＝24GPa，排水剪切模量为μ＝15.2GPa，

②依据技术方案步骤(2)所示，建立圆柱形岩石样品模型，如图2所示。样品(R1)上端黏贴铝块A1、下端黏贴铝块A2。岩石样品R1的中心轴向添加长度5mm线段S1，表示测量样品轴向应变的应变片，径向添加长度5mm线段S2，表示测量样品径向应变的应变片。铝块A1偏上部位添加长度为5mm的线段S3，表示测量铝块轴向应变的应变片。

③依据技术方案步骤(3)所示，对铝块A1和A2添加支配物理场和边界条件，如图1所示，对于铝块A1而言，其上表面加载周期震荡轴向应力0.1[MPa]，对于铝块A2而言，其下表面为固定位置，采用固定约束边界条件，即u^s＝[U，V，W]＝0；

③依据技术方案步骤(4)所示，对岩石样品添加支配物理场，如方程(2)和(3)所示，频率按照以10为底的指数变化，初始值为-1，间隔为0.1，最大值为3，因此频率范围为10^-2：0.1：3Hz，角频率为ω＝2πf,样品内部流体分布由CT扫描获得，流体体积模量、粘度和密度与流体分布密切相关。样品处于不排水边界条件，整个样品表面的流量均为0。

④依据技术方案步骤(5)所示，划分网格，铝块A1为均匀四面体网格，应变片S1、S2和S3位置按线性网格加密；岩石样品R1网格按照样品非均质程度自适应加密；样品下部铝块A2为固定样品的作用，不参与纵横波速度的计算，因此为粗网格，最终结果如图3所示；

⑤按照技术方案步骤(6)所示，依据网格离散方程组，并利用混合求解方法求解，各物理场迭代求解次数不超过50次，有效收敛求解时间约7个小时，各物理场收敛曲线如图4所示；

⑥依据技术方案步骤(7)所示，计算纵横波速度，图5为预测结果与实测结果对比图，看到二者匹配良好，证明了本发明的有效性。

对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。

Claims (7)

1.一种高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述方法包括：

S1获取建立模型所需要的参数；其中，需要岩石的物性参数作为输入，可通过实验室测量获取；

S2依据岩石样品尺寸，按照岩石样品R1上端黏贴铝块A1、下端黏贴铝块A2的方法构建模型；

S3设置铝块样品支配方程及边界条件；

S4设置岩石样品的支配方程及边界条件；

S5网格剖分；其中，应变片S1和S2进行细化网格剖分；上部铝块A1和岩石样品R1进行常规四面体网格剖分；根据流体的分布特征及样品的非均质性特征，进行精细化网格的剖分；下部铝块A2按照样品R1下表面网格进行延拓，进行三维空间的剖分；

S6依据网格离散方程组，并利用混合求解方法求解；

S7计算纵横波速度。

2.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S1中，包括：

S1.1参考铝制样，长度为40mm，直径为40mm，密度为：

ρ_al＝2700Kg*m^-3；

S1.2使用游标卡尺测量岩石样品的直径DR，长度LenR，利用天平测量干燥岩石质量MR，计算干燥岩石密度：

S1.3依据SY/T 6385-2016《覆压下岩石孔隙度和渗透率测定方法》标准，测定岩石样品的孔隙度φ和渗透率k；

S1.4依据分析标准SY/T 5163-2018，对岩石进行矿物分析；通过矿物分析和VRH公式，计算岩石颗粒体积模量K_g；

S1.5依据SY/T 7410.1-2018分析标准对岩石孔隙结构和流体分布进行扫描和重建，进而获得岩石的孔隙结构和流体分布特征；

S1.6依据流体分布特征计算流体的体积模量K_f，粘度η；

S1.7依据行波法，测得样品的干燥体积模量K_d和剪切模量μ。

3.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S2中，铝块A1的中心和岩石样品R1的中心添加沿着纵向和横向分别添加两个线段S1和S2，用于表示测量轴向应变和径向应变的应变片。

4.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S3中，包括：

铝块A1和A2是各向同性均匀材料，其支配方程为应力-应变方程，为提高计算效率，采用频率域形式:

其中ρ_al是铝制材料的密度，u^s＝[U，V，W]是铝制材料的位移场，σ是应力场张量；对于铝块A1而言，其上表面加载周期震荡震源，频率域为常数，即轴向加载应力σ_zz＝0.1[MPa]，对于铝块A2而言，其下表面为固定位置，其位移张量u^s＝0。

5.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S4中，包括：

岩石样品中流固耦合应该符合孔弹性方程：

其中

是哈密顿算子，ω是角频率；ρ_f是流体密度，ρ_s是固体骨架密度，ρ_b＝ρ_s+φρ_f，

φ是孔隙度，S是孔隙的曲折因子，η是流体的粘度，κ是渗透率，

是虚数单位；P_f是流体压力，u^s＝[U，V，W]是固体骨架位移，σ是应力张量，与流体压力P_f和固体骨架位移u^s，干燥岩石体积模量K_d，干燥岩石剪切模量μ有关；其中α是Biot-Willis系数，M为耦合弹性模量：

岩石样品R1上下表面被铝块A1、A2封闭，侧面被胶套封闭，因此岩石各边界处于不排水条件，所有流量应为0，具体为：

其中n是R1表面外法向向量，

是虚数单位。

6.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S5中，包括：

利用网格划分结果，通过离散化转为线性方程组：

其中

采用混合方法，即牛顿迭代法和LU矩阵分解法混合的方式求解线性方程组，具体步骤为：

A设置初始条件X^i＝0＝0；

B计算L₁和雅可比矩阵

则

C基于

计算L₂和雅可比矩阵

则

D基于

计算L₃和雅可比矩阵

则

E基于

利用LU矩阵分解法求解

F如果|Xⁱ⁺¹-Xⁱ|＜∈，则停止计算，否则回到步骤A；∈是相对误差。

7.根据权利要求1所述的高速低内存消耗的预测储层岩石弹性波速度的方法，其特征在于，所述步骤S7中，包括：

依据应变计算杨氏模量E和泊松比v：

其中杨氏模量E_al为铝块杨氏模量72GPa，

为S1所在位置轴向应变的平均值，

为S2所在位置轴向应变的平均值，

为S3位置径向应变的平均值；

纵波速度计算结果为：

横波速度计算结果为：

Images