CN112326524A - 一种基于ct扫描图像的岩石孔渗测量方法 - Google Patents

Abstract

本发明提供一种基于CT扫描图像的岩石孔渗测量方法，包括：对岩样进行CT扫描得到扫描结果并构建常温常压数字岩心模型，并对常温常压数字岩心进行三维重建，得到三维数字岩心模型；采用离散单元对三维数字岩心模型进行填充，得到离散单元位置和体积信息；构建考虑离散单元粘结作用的接触模型，计算高温高压下三维数字岩心模型中离散单元的位移，并进行离散单元位置更新，并得到高温高压三维数字岩心模型；根据常温常压数字岩心模型及高温高压三维数字岩心模型，得到三维孔隙网络模型，计算孔隙结构特征参数以及渗透率。本发明克服了目前CT扫描无法在高温高压条件下开展，高温高压下岩石孔隙结构变化无法分析的困难，简化了计算过程。

Description

一种基于CT扫描图像的岩石孔渗测量方法

技术领域

本发明属于非常规油气开发技术领域，具体涉及一种基于CT扫描图像的高温高压高应力条件下的岩石孔隙结构表征与渗透率测量方法。

背景技术

近年来，我国油气资源消耗量日益增加，石油对外依存度逐年攀升。同时，伴随着中浅层油气田的长期勘探与开采，在盆地中浅层取得大的突破已越来越困难，油气勘探目的层逐渐由中浅层转向深层(3500～4500m)和超深层(大于4500m)。实现深层超深层油气能源的接替对我国的能源安全和可持续发展具有重要意义。我国深层超深层油气藏开采潜力大，但面临“三高”(高温、高压、高应力)难题，必须考虑热流固耦合效应对岩石孔隙结构的影响。认识真实岩心的微观孔隙结构是研究渗流机理和提高开发效果的基础。CT扫描技术通过对实验用岩石进行扫描，获取岩石三维扫描图像，通过图像处理方法构建反映孔隙和固体骨架特征的数字岩心。基于数字岩心提取孔隙网络模型，在孔隙网络模型基础上分析扫描区域孔隙结构性质来对地下岩石孔隙结构特征分析以及渗透率计算。目前该技术受限于CT扫描设备条件，仅能对岩石在常温常压下的孔隙结构以及渗透率进行表征研究，无法直接进行高温高压下的岩心扫描，因此目前的CT扫描技术无法分析高温高压高应力条件下的岩石孔隙结构特征与渗透率。

数值模拟手段是解决高温高压高应力条件下的岩石孔隙结构特征与渗透率变化的有效方法。通过考虑温度-流动-应力影响的热流固耦合数值模拟方法可以充分揭示岩石孔隙结构特征与渗透率变化规律。常见热流固耦合数值模拟方法主要有格子玻尔兹曼方法与离散单元法耦合，计算流体动力学和离散单元法耦合及光滑粒子流体动力学与离散单元法耦合。上述热流固耦合数值模拟方法存在计算量大，耗费时间长的问题，不适用于解决高温高压高应力条件下的岩石孔隙结构特征与渗透率变化模拟计算问题。

发明内容

为解决现有技术存在的问题，本发明提供了一种基于CT扫描图像的高温高压高应力条件下的岩石孔隙结构表征与渗透率测量方法，通过采用离散单元法模拟高温高压高应力条件下岩石形变，采用孔隙网络模型方法快速分析计算岩石形变过程中的孔隙结构与渗透率变化，大大简化了渗透率计算过程的同时保证了计算的准确性。

为实现上述目的，本发明提供了如下方案：本发明提供一种基于CT扫描图像的岩石孔渗测量方法，其特征在于，包括以下步骤：

对岩样进行CT扫描得到扫描结果，根据扫描结果构建常温常压数字岩心模型，并对常温常压数字岩心进行三维重建，得到三维数字岩心模型；

采用离散单元对所述三维数字岩心模型进行填充，得到离散单元位置和体积信息；

构建考虑离散单元粘结作用的接触模型，根据所述接触模型计算高温高压下所述三维数字岩心模型中离散单元的位移，根据所述位移进行所述离散单元位置更新，并根据更新后的位置及所述体积信息得到高温高压三维数字岩心模型；

根据所述常温常压数字岩心模型及所述高温高压三维数字岩心模型，得到三维孔隙网络模型，根据所述三维孔隙网络模型计算孔隙结构特征参数以及渗透率。

优选地，所述常温常压下的数字岩心模型的构建方法为：根据岩样孔径分布数量级选取对应的CT扫描设备分辨率，并对岩样进行CT扫描，得到扫描岩样灰度数据体；计算所述扫描岩样灰度数据体的灰度值分布，根据灰度值分布及实验室测量得到的孔隙度结果选取灰度阈值；根据所述灰度阈值对扫描岩样固体骨架区域及孔隙区域进行阈值分割，得到常温常压下的数字岩心模型。

优选地，所述考虑离散单元粘结作用的接触模型的构建过程为：

根据岩石力学试验测量得到离散单元的物性参数；分别采用三轴压缩力学试验及数值模拟计算测试获得岩石岩样的应力-应变曲线，并根据岩石岩样的应力-应变曲线进行校正单元间粘结作用力的校正，得到离散单元之间的粘结键大小；根据所述离散单元的物性参数及所述离散单元之间的粘结键构建所述考虑离散单元粘结作用的接触模型。

优选地，所述离散单元的位移的计算方法为：根据所述考虑离散单元粘结作用的接触模型计算在高温高压过程中所述离散单元受到外部压力作用时通过挤压碰撞传递的压力，并与重力、所述离散单元之间已有的粘结力形成合力，得到每两个离散单元之间的合外力作用，并根据合外力作用得到所述离散单元的位移。

优选地，所述渗透率的计算过程为：根据所述三维孔隙网络模型分别计算单个圆形孔隙的流体流量、单个正方形孔隙的流体流量及单个三角形孔隙的流体流量并分别得到圆形孔隙、正方形孔隙及三角形孔隙的传导率；根据所述圆形孔隙、正方形孔隙及三角形孔隙的传导率计算所述三维孔隙网络模型上每个孔隙流量和压力分布；根据所述流体流量和所述压力分布计算给定条件下的所述孔隙网络模型的流体流量，并根据所述孔隙网络模型的流体流量计算渗透率。

优选地，所述三维孔隙网络结构模型采用最大球方法构建。

本发明公开了以下技术效果：

(1)本发明通过在常温常压下对岩样进行扫描获得常温常压下数字岩心与离散单元方法结合模拟高温高压下岩石孔隙结构变化，克服了目前CT扫描无法在高温高压条件下开展，高温高压下岩石孔隙结构变化无法分析的困难；

(2)本发明方法在节省了计算时间的同时保证了计算精度，传统计算高压下的数字岩心渗透率基于离散单元法耦合格子玻尔兹曼方法、光滑粒子动力学方法展开，计算时间长，一般需要数周的时间，并且需要巨大的高性能计算机集群作为计算资源支撑；本发明通过引入离散单元法耦合孔隙网络模型方法来准确计算数字岩心渗透率，仅需要几分钟时间进行计算，大大简化了渗透率计算过程的同时保证了计算的准确性；

(3)传统岩心实验方法仅能测量得到不同温度压力下的孔隙空间统计信息(孔隙半径分布)，无法测量得到孔隙形状以及孔隙空间联通性变化。本发明基于不同温度阶段下的数字岩心孔隙空间，通过提取孔隙网络模型能够准确定量计算孔隙空间统计参数(孔隙半径分布、形状因子分布)以及空间结构参数(配位数分布)，克服了目前实验方法的局限性，并且相比实验方法具有时间短、费用低的优点，适合大规模推广应用。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明基于CT扫描图像的岩石孔渗测量方法流程示意图；

图2为本发明实施例高温高压数字岩心的模型构建流程示意图；

图3为本发明实施例离散单元结构示意图；

图4为本发明实施例孔喉单元结构示意图；

图5为本发明是实施例孔隙空间的简化与不同规则几何体横截面对应的形状因子示意图；

图6为本发明实施例CT扫描岩样阈值分割后的数字岩心结构示意图，其中，黑色表示骨架，灰色表示孔隙；其中图6(a)为CT扫描岩样阈值分割后的数字岩心结构示意图；图6(b)为通过离散单元法模拟得到的高温高压中间阶段(上覆岩层压力10MPa)岩石变形后的数字岩心结构示意图；图6(c)为通过离散单元法模拟得到的高温高压最终阶段(上覆岩层压力40MPa)岩石变形后的数字岩心结构示意图；

图7为本发明实施例CT扫描岩样阈值分割后的数字岩心中的孔隙空间结构示意图，其中，黑色表示骨架，灰色表示孔隙；其中：图7(a)为CT扫描岩样阈值分割后的数字岩心中的孔隙空间结构示意图；图7(b)为通过离散单元法模拟得到的高温高压中间阶段(上覆岩层压力10MPa)岩石变形后的数字岩心中的孔隙空间结构示意图；图7(c)为离散单元法模拟得到的高温高压最终阶段(上覆岩层压力40MPa)岩石变形后的数字岩心中的孔隙空间结构示意图；

图8为本发明实施例基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型结构示意图，其中，黑色表示骨架，灰色表示孔隙；其中：图8(a)为基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型结构示意图；图8(b)为基于离散单元法模拟得到的高温高压中间阶段(上覆岩层压力10MPa)岩石变形后的数字岩心提取得到的孔隙网络模型结构示意图；图8(c)为基于离散单元法模拟得到的高温高压最终阶段(上覆岩层压力40MPa)岩石变形后的数字岩心提取得到的孔隙网络模型结构示意图；

图9为本发明实施例基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型孔隙半径分布示意图，其中，黑色表示骨架，灰色表示孔隙；其中：图9(a)为基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型孔隙半径分布示意图；图9(b)为基于离散单元法模拟得到的高温高压中间阶段(上覆岩层压力10MPa)岩石变形后的数字岩心提取得到的孔隙网络模型孔隙半径分布示意图；图9(c)为基于离散单元法模拟得到的高温高压最终阶段(上覆岩层压力40MPa)岩石变形后的数字岩心提取得到的孔隙网络模型孔隙半径分布示意图；

图10为本发明实施例基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型配位数分布示意图，其中，黑色表示骨架，灰色表示孔隙；其中：图10(a)为基于CT扫描岩样阈值分割后的数字岩心提取得到的孔隙网络模型配位数分布示意图；图10(b)为基于离散单元法模拟得到的高温高压中间阶段(上覆岩层压力10MPa)岩石变形后的数字岩心提取得到的孔隙网络模型配位数分布示意图；图10(c)为基于离散单元法模拟得到的高温高压最终阶段(上覆岩层压力40MPa)岩石变形后的数字岩心提取得到的孔隙网络模型配位数分布示意图。

具体实施方式

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

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例1

如图1所示，本发明提供一种基于CT扫描图像的岩石孔隙结构与渗透率测量方法，包括以下步骤：

S1、对岩样进行CT扫描，根据扫描结果，构建常温常压数字岩心模型，并对常温常压数字岩心进行三维重建，得到三维数字岩心模型。

根据岩样孔径分布数量级选取对应的分辨率Ro对岩样进行CT扫描，得到扫描岩样灰度数据体，其中体素尺寸为N_x×N_y×N_z，对应的物理尺寸为N_xRo×N_yRo×N_zRo。

计算扫描岩样灰度数据体的灰度值分布，根据灰度值分布及实验室测量得到的孔隙度结果选取灰度阈值；根据灰度阈值对扫描岩样固体骨架区域、孔隙区域进行阈值分割，得到常温常压数字岩心模型。

本实施例在常温常压下采用分辨率为1μm的CT扫描岩样、阈值分割后的数字岩心，其中体素大小为200×200×200。

通过Avizo软件重建常温常压数字岩心模型，得到三维数字岩心模型。三维数字岩心几何体数据文件包含整个岩心每个像素的三维位置坐标信息。

S2、采用离散单元对三维数字岩心模型进行填充，得到离散单元位置和体积信息。

采用离散单元对三维数字岩心模型进行填充，即在每个像素点生成一个球体颗粒，相互接触的单元体之间设置有碰撞或者粘结作用力，并记录球体颗粒的位置和体积信息。离散单元体的球心为像素中心，离散单元体体积为已知。本实施例构建的离散单元结构如图3所示。

S3、根据岩石物理参数数据计算离散单元填充后的三维数字岩心模型的离散单元之间的粘结键大小和物性参数，根据粘结键大小和物性参数构建考虑粘结作用的离散单元接触模型，根据接触模型计算高温高压下三维数字岩心中离散单元的位移，并进行离散单元位置更新，得到高温高压三维数字岩心模型。

由于步骤S2中得到的离散单元填充后的三维数字岩心模型为一般条件下的通用模型，而对于具体的不同的岩层情况中离散单元的物性参数及单元间实际粘结键大小则需要通过岩石物理参数数据进行校正。本实施例根据岩石力学试验测量得到如杨氏模量，泊松比，导热系数及密度等离散单元的物性参数；采用三轴压缩等力学试验测试获得岩样的应力-应变曲线，采用数值模拟计算上述岩样的应力-应变曲线，多次校正单元间粘结作用力，直到两次应力-应变曲线拟合程度在合理范围，得到离散单元之间的粘结键大小，从而获得考虑离散单元粘结作用的接触模型中的必要参数。

本实施例构建的接触模型的过程为：

计算相邻离散单元之间的法向力F_n。法向力F_n是相邻离散单元之间重叠量δ_n的函数，具体为：

其中：E^*为当量杨氏模型，R^*为当量半径，具体为：

其中：E_i、n_i、R_i和E_j、n_j、R_j分别代表互相接触的单元i和j的杨氏模量、泊松比和球体半径。

计算离散单元的法向阻尼力

具体为：

其中：

是当量质量；

是相对速度的法向分量；S_n为法向刚度；β为常数

e为恢复系数。

计算离散单元的切向力F_t，具体为：

F_t＝-S_tδ_t (5)

其中：δ_t为切向重叠量，S_t为切向刚度，G^*为当量剪切模量。

计算离散单元的切向阻尼力

具体为：

本实施例的离散单元接触模型采用有限大小的粘结键承受切向和法向的位移，直到达到最大的法向和切向剪切应力，即粘结断裂点，具体为：

其中：R_B为粘结半径；S_n、S_t分别是法向和切向刚度；d_t为时间步长；v_n、v_t分别是法向和切向速度；ω_n、ω_t分别是法向切向角速度；dF_n和dF_t为最大法向和切向剪切应力；dM_n和dM_t为法向和切向扭矩；A是有效接触面积；J为转动惯量。

当法向和切向剪切应力超过某个定义的值时，粘结键断裂，具体为：

其中：d_max和t_max分别为法向和切向剪切应力临界值。

根据构建的接触模型计算在高温高压过程中三维数字岩心每一个时刻内的离散单元体积，并记录速度和位置。即当受到外部压力作用时，离散单元之间通过挤压碰撞传递压力，并与重力、离散单元之间已有的粘结力等形成合力，得到每两个离散单元之间的合外力作用，根据牛顿运动力学定律可计算离散单元从t到t+1时间内的加速度、速度及位移。

离散单元i,j的热流量Q_i.j计算：

Q_i,j＝h_cDT_i,j (11)

传热系数h_c为：

离散单元温度随时间更新根据：

其中：r^*为离散单元几何平均半径；E^*是有效杨氏模量；m_i,C_i,T分别是离散单元的质量、比热容和温度；Q_heat代表对流和传导热通量的总和。

根据式(11)-(13)可以计算并更新各个离散单元颗粒的温度，从而得到岩心的温度场分布，这是进行孔隙网络流动模拟时所需要的重要初始边界条件。

在离散单元填充的数字岩心三维空间中，通过编制Python程序遍历每一个三维坐标，根据已知离散单元位置信息、球体中心和半径关系判断每一个三维坐标在球体内还球体外。若三维坐标点距离球心的距离小于离散半径，则代表此坐标点位于离散单元占据的空间内，属于数字岩心的骨架区域，则用骨架符号“1”表示，反之其余空间坐标点用孔隙符号“0”表示。

S5、根据常温常压下的数字岩心及高温高压三维数字岩心，得到不同温度压力阶段的三维孔隙网络模型，计算孔隙结构特征参数以及渗透率，本实施例构建的三维孔隙网络模型如图8所示。

根据常温常压数字岩心模型以及不同温度压力阶段的三维数字岩心模型，采用最大球方法提取得到保留孔隙拓扑结构的孔隙网络模型。每个孔隙和喉道形状根据形状因子G确定，形状因子G通过下式进行计算：

其中：A为孔隙孔喉截面面积；P_d为孔隙孔喉截面形状的周长。

根据孔隙网络模型上每一个孔隙喉道形状因子取值以及不同规则形状孔隙的形状因子取值范围(表1)，对每个不规则孔隙喉道形状进行简化，本实施例的孔喉单元结构如图4所示。

表1

对于截面形状为三角形的孔隙、喉道单元体，其三个内角的确定方法为：

假设三角形的三个内角分别为β₁、β₂和β₃，并有β₁<β₂<β₃，则计算β₁、β₂、β₃的具体步骤为：

根据给定的形状因子数值G，首先计算β₂的取值区间[β_2,min,β_2,max]，具体为：

根据式(15)、(16)的取值区间随机选取β₂的数值。

然后，计算β₁，具体为：

最后，可得到β₃的值：

由于孔隙网络模型中的孔隙、喉道具有和原始数字岩心孔隙空间等效的形状因子，因此孔隙网络模型能够准确反映岩心孔隙结构特征。

根据三维孔隙网络模型每一个孔隙的尺寸、形状、联通性计算整体三维孔隙网络模型的平均孔径、平均配位数、不同形状孔隙所占比例及渗透率。

基于孔隙网络模型的渗透率计算主要过程为：

单个圆形孔隙中流体流量随压差变化关系可表示为：

其中：△P为单个孔隙两端压力降；l为单个孔隙长度；μ为流体粘度，高温高压条件下流体粘度值根据实验室流体粘度测量结果进行计算。

单个正方形孔隙中流体流量q_squ随压差变化关系可表示为；

其中：w为正方形孔隙边长。

单个三角形孔隙中流体流量随压差q_tri变化的变化关系可表示为：

通过传导率的概念来描述流体在单个孔隙中的流动能力，传导率g定义为：

其中：q表示单个孔隙中流体流量；g表示单位压差下流体通过单个孔隙中的流量。

圆形孔隙流体传导率g_cir、正方形孔隙流体传导率g_squ、三角形孔隙流体传导率g_tri可分别表示为：

对于每个孔隙流体流入流出量相等，因此可表示为：

Q_ij＝g_ij(p_i-p_j) (27)

其中：N_i表示孔隙网络模型中所有孔隙个数；p_i和p_j分别表示孔隙i和孔隙j上的压力；Q_ij表示孔隙i与孔隙j之间的流体流量。根据图4，相连孔隙间传导率g_ij可表示为：

其中：L_i，L_j，L_t分别表示孔隙i半径、孔隙j半径、喉道长度；g_i，g_j，g_t分别表示孔隙i传导率、孔隙j传导率、喉道传导率。

给定孔隙压力、地层温度、压力梯度即可对孔隙网络整体求解得到定压差条件下通过孔隙网络模型的流体流量Q_t。

渗透率k基于流体流量Q_t，采用达西公式进行计算，具体为：

其中：N_inlet为入口孔隙个数；q_inlet为每个入口孔隙上的流体流量；L_pnm为孔隙网络模型长度；A_pnm为孔隙网络模型截面面积；△P_pnm为孔隙网络模型上压力降。

本实施例给定地层温度为400K，压力梯度为0.1MPa/m，并计算不同温度压力阶段的渗透率、孔隙半径分布、配位数分布以及孔隙形状因子，计算结果见表2、图9及图10。

表2

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (6)

1.一种基于CT扫描图像的岩石孔渗测量方法，其特征在于，包括以下步骤：

对岩样进行CT扫描得到扫描结果，根据扫描结果构建常温常压数字岩心模型，并对常温常压数字岩心进行三维重建，得到三维数字岩心模型；

采用离散单元对所述三维数字岩心模型进行填充，得到离散单元位置和体积信息；

构建考虑离散单元粘结作用的接触模型，根据所述接触模型计算高温高压下所述三维数字岩心模型中离散单元的位移，根据所述位移进行所述离散单元位置更新，并根据更新后的位置及所述体积信息得到高温高压三维数字岩心模型；

根据所述常温常压数字岩心模型及所述高温高压三维数字岩心模型，得到三维孔隙网络模型，根据所述三维孔隙网络模型计算孔隙结构特征参数以及渗透率。

2.根据权利要求1所述的基于CT扫描图像的岩石孔渗测量方法，其特征在于，所述常温常压下的数字岩心模型的构建方法为：根据岩样孔径分布数量级选取对应的CT扫描设备分辨率，并对岩样进行CT扫描，得到扫描岩样灰度数据体；计算所述扫描岩样灰度数据体的灰度值分布，根据灰度值分布及实验室测量得到的孔隙度结果选取灰度阈值；根据所述灰度阈值对扫描岩样固体骨架区域及孔隙区域进行阈值分割，得到常温常压下的数字岩心模型。

3.根据权利要求1所述的基于CT扫描图像的岩石孔渗测量方法，其特征在于，所述考虑离散单元粘结作用的接触模型的构建过程为：

根据岩石力学试验测量得到离散单元的物性参数；分别采用三轴压缩力学试验及数值模拟计算测试获得岩石岩样的应力-应变曲线，并根据岩石岩样的应力-应变曲线进行校正单元间粘结作用力的校正，得到离散单元之间的粘结键大小；根据所述离散单元的物性参数及所述离散单元之间的粘结键构建所述考虑离散单元粘结作用的接触模型。

4.根据权利要求1所述的基于CT扫描图像的岩石孔渗测量方法，其特征在于，所述离散单元的位移的计算方法为：

根据所述考虑离散单元粘结作用的接触模型计算在高温高压过程中所述离散单元受到外部压力作用时通过挤压碰撞传递的压力，并与重力、所述离散单元之间已有的粘结力形成合力，得到每两个离散单元之间的合外力作用，并根据合外力作用得到所述离散单元的位移。

5.根据权利要求1所述的基于CT扫描图像的岩石孔渗测量方法，其特征在于，所述渗透率的计算过程为：根据所述三维孔隙网络模型分别计算单个圆形孔隙的流体流量、单个正方形孔隙的流体流量及单个三角形孔隙的流体流量并分别得到圆形孔隙、正方形孔隙及三角形孔隙的传导率；根据所述圆形孔隙、正方形孔隙及三角形孔隙的传导率计算所述三维孔隙网络模型上每个孔隙流量和压力分布；根据所述流体流量和所述压力分布计算给定条件下的所述孔隙网络模型的流体流量，并根据所述孔隙网络模型的流体流量计算渗透率。

6.根据权利要求1所述的基于CT扫描图像的岩石孔渗测量方法，其特征在于，所述三维孔隙网络结构模型采用最大球方法构建。

Images