利用不同分辨率CT图像构建多尺度三维数字岩心的方法与流程

本发明属于石油测井解释领域，涉及一种利用不同分辨率CT图像构建多尺度三维数字岩心的方法。

背景技术：

测井目前面临着“三低两复杂”和“非常规”储层等地质问题的挑战，进一步加强对储层岩石物理性质的研究是提高测井解释符合率、准确评价油气产能的关键。储层岩石物理性质的深入研究可以基于以下两种手段：一是基于三维数字岩心的数值模拟方法，以数字岩心技术为基础，通过数值模拟获取储层不同尺度上的岩石物理特征参数，从而建立准确、有效的测井解释模型；二是物理模拟方法，研制新型的井下岩石物理测量装置，开展地层原始条件下的岩石物理定点、定向、定域测量，准确获取地层原始参数。由于井下岩石物理测量还处于探索阶段，因此目前最常用的方法还是通过数值模拟手段对储层岩石物理性质进行深入研究。X射线CT是一种无损三维成像方法，通过扫描实验获取岩心的三维灰度图像，灰度值反映了岩石不同组分对X射线吸收系数的差异，可区分岩石骨架和孔隙空间。再采用图像分割技术，将灰度图像转换为二值化图像，建立起三维数字岩心。对于碳酸盐储层非均质性强，孔隙类型多元化(孔、裂缝、洞等)，孔隙大小变化可达好几个数量级，不同尺度上的孔隙对储层的岩石物理特性均有影响，需开发多尺度融合算法解决CT扫描岩心尺寸与分辨率矛盾问题，拓展数字岩心技术在该领域中的应用。

技术实现要素：

本发明的目的在于提供一种利用不同分辨率CT图像构建多尺度三维数字岩心的方法，该方法能够解决CT扫描岩心尺寸与分辨率矛盾问题，拓展数字岩心技术在非均质储层评价中的应用。

为达到上述目的，本发明采用以下技术方案予以实现：

利用不同分辨率CT图像构建多尺度三维数字岩心的方法，包括以下步骤：

1)利用CT扫描对岩心进行不同分辨率扫描，获取高低两种分辨率的三维数字岩心；

2)基于三维数字岩心利用最大球算法提取岩心的三维孔隙空间网络，分别得到粗尺度和细尺度的孔隙网络模型；

3)将利用最大球算法提取的粗尺度孔隙网络模型直接移植过来作为新模型中的粗尺度模型；结合细尺度网络的几何拓扑结构信息随机重建等效细尺度孔隙网络；

4)通过孔隙网络跨尺度耦合实现粗尺度孔隙网络和等效细尺度孔隙网络的叠加。

本发明进一步的改进在于：

所述步骤3)的具体方法为：

3-1)首先利用低分辨率图像灰度与孔隙度概率的分布曲线确定微孔在粗尺度网络中的填充区域；

通过下式(1)求出每个微孔灰度值为孔隙的概率P_φ：

其中，T₁为灰度值阈值下限，T₂为灰度值阈值上限，n为微孔灰度值；

假设多尺度孔隙网络的孔隙空间为Ω，粗尺度孔隙网络孔隙空间为Ω₁，骨架的分布空间为Ω₃，细尺度网络的填充空间域Ω₂为：

Ω₂＝Ω-Ω₁-Ω₃ (2)

3-2)利用尺度变化因子和连通函数确定在粗尺度网络中需要生成的微孔隙数目和需要充填的喉道数目；根据配位数分布，按孔隙半径从大到小依次添加孔隙之间的连通喉道，并保证添加喉道的几何特征同原大孔隙网络一致，进而构建出同原始孔隙网络几何拓扑等效的随机孔隙网络模型。

所述步骤3-2)中，微孔隙数目N'_p由式(3)求得：

N'_p＝ξ³N_p (3)

其中，N_p为原始细尺度网络中的孔隙数目，ξ为尺度变化因子。

计算原细尺度孔隙网络中孔喉参数的相关性，两个随机变量S和T间的相关性用相关系数ρ_(S,T)进行度量，对于S,T的序列{(s₁,t₁),...,(s_n,t_n)}，ρ定义为：

其中，s、t分别表示两个随机变量，n表示序列里面变量个数；相关系数绝对值大表示相关性强；通过计算相关系数，确定各参数之间的从属关系。

所述步骤3-2)中，根据原细尺度网络的连通函数计算需要充填的喉道数目：

描述网络拓扑特征的连通性函数可通过简化的比欧拉示性数来描述孔隙的全局连通性由连通性函数来表示：

式中，N_p(r)为原始孔隙网络孔隙集合中半径大于r的孔隙的数目，N_t(r)为原始孔隙网络喉道集合中半径大于r的喉道数目，V为原始孔隙网络的体积；根据欧拉数定义得到需要填充的喉道数目N_t'(r_k)表示为：

N_t'(r_k)＝N'_p(r_k)-χ_v(r_k)V' (6)

式中，N'_p(r_k)表示从新建孔隙集合p'中孔隙半径大于r_k的孔隙数目，V'为新建网络的体积。

利用微孔隙比例因子提高建模速度；

引入微孔隙比例因子的概念，微孔隙比例因子是指在构建多尺度孔隙网络模型过程中实际填充的微孔隙数目N”_pm占通过尺度变化因子计算所得微孔隙数目N'_pm的比例，定义如下：

α＝N”_pm/N'_pm×100％ (7)

结合微孔隙比例因子的定义，根据原始细尺度孔隙网络模型中孔隙的几何参数变量之间的概率分布函数和变量之间的相关性信息，能够在细尺度孔隙网络模型空间域中生成同原始细尺度孔隙网络几何拓扑结构信息等效的微孔隙，并根据统计分析随机产生各微孔隙单元半径、体积的形状因子，使各变量的概率分布与原始细尺度孔隙网络模型一样。

所述步骤4)的具体方法为：

4-1)确定跨尺度连接的最大孔隙间距：

对于一个喉道连接的两个孔隙，半径为分别为r₁,r₂，定义孔隙间距d_i：

d_i(N₁,N₂)＝d(N₁,N₂)-r₁-r₂ (8)

其中，N₁、N₂分别为任意两个孔隙的编号，d(N₁,N₂)为两个孔隙中心点的距离；遍历原始细尺度孔隙网络的孔隙间距，将max(d_i(N₁,N₂),i＝1,...,n)作为跨尺度连接的孔隙间距的最大值；

4-2)确定跨尺度配位数：

跨尺度配位数定义为一个粗尺度孔隙所能连接的细尺度孔隙的最大数目；对于原始细尺度孔隙网络，其一系列的孔隙{N₁,N₂,…,N_n}和喉道{B₁,B₂,…,B_m}，在区域的中心插入s个同心圆球，同时保证最大的球仍在孔隙网络内；计算数目N_Bj,j＝1,…,s，即穿透相应圆球的喉道数目；也即，如果d(p,N)为同心圆球p和孔隙N间的距离，N_k1和N_k2为喉道B_k连接的两孔隙，并且d(p,N_k1)＜d(p,N_k2)，则：N_Bj＝{B_k|d(p,N_k1)＜r_j＜d(p,N_k2),k＝1,2…,m}；定义跨尺度系数γ_CSC

其中，a_j为第j个球的表面积，s为同心圆的数量，N_Bj为穿透相应圆球的喉道数目；跨尺度系数表示穿透单位孔隙面积的喉道数；

为了确定交叉尺度配位数，考虑粗尺度孔隙网络中的一个孔隙，它的一部分表面积已被粗尺度的喉道占据；对于每一个粗尺度的孔隙N_c，对应半径为r_c，粗尺度配位数z_c，其交叉尺度配位数的计算公式为：

其中，A_k为半径为r_k的粗尺度喉道在球面上所占的面积：

与现有技术相比，本发明具有以下有益效果：

本发明解决了非均质性储层CT扫描岩心尺寸与分辨率矛盾问题。对于国外算法，微孔隙在岩心整个空间中均匀分布(包含骨架)，这与微孔隙的真实分布状态相违背，针对这个问题，发明了微孔非均匀分布多尺度三维数字岩心构建方法。利用本发明构建的数字岩心更符合实际岩心的特征，实现了不同分辨率CT图像的融合构建多尺度三维数字岩心。

【附图说明】

图1是利用CT扫描构建的数字岩心(左低分辨率扫描，右高分辨率扫描)；

图2是基于数字岩心提取的粗尺度孔隙网络模型；

图3是基于数字岩心提取的细尺度孔隙网络模型；

图4是微孔非均匀分布的等效细尺度孔隙网络模型；

图5(a)是低分辨率扫描图像二维切片；

图5(b)是图像灰度-孔隙度概率对应关系；

图6(a)是没有连接之前大小孔隙相互独立，其中N1、N2、N3代表大尺寸孔隙；

图6(b)是大孔隙搜索某一半径内的所有小孔隙(阴影部分)；

图6(c)是大孔隙和搜索半径中的小孔隙用喉道连接起来；

图7是微孔非均匀分布的多尺度孔隙网络模型。

【具体实施方式】

下面结合附图对本发明做进一步详细描述：

本发明利用不同分辨率CT图像构建多尺度三维数字岩心的方法，包括以下步骤：

1)利用CT扫描对岩心进行不同分辨率扫描，获取高低两种分辨率的三维数字岩心。

2)基于三维数字岩心利用最大球算法提取岩心的三维孔隙空间网络，分别得到粗尺度和细尺度的孔隙网络模型。

3)将利用最大球算法提取的粗尺度孔隙网络模型直接移植过来作为新模型中的粗尺度模型；结合细尺度网络的几何拓扑结构信息随机重建等效细尺度孔隙网络，具体步骤如下：

a.确定细尺度孔隙网络填充的空间域。细尺度孔隙网络在填充过程中不能与粗尺度孔隙网络产生交叠，因此细尺度孔隙网络填充的空间域为粗尺度孔隙网络中大孔隙以外的空间。考虑岩石孔隙的真实分布状态与规律，认为低分辨率未识别的微孔隙只存在于尺寸小于粗尺度网络分辨率的微孔区域内。粗尺度孔隙网络由低分辨率扫描得到的数字岩心提取得到，由于受到扫描分辨率的限制，尺寸小于分辨率的孔隙无法精确识别，在成像图上表现为灰度值大于可识别孔隙的灰度值但是低于纯岩石骨架灰度值，通过确定成像灰度与孔隙度的关系可以看出，存在阈值下限T₁和阈值上限T₂使得微孔隙的分布范围满足灰度值T属于区间[T₁,T₂]。同时，通过确定阈值上下限，可以给出微孔隙分布区间每个位置处的孔隙度值：

假设多尺度孔隙网络的孔隙空间为Ω，粗尺度孔隙网络孔隙空间为Ω₁，细尺度网络的填充空间域为：

Ω₂＝Ω-Ω₁-Ω₃

其中Ω₃为骨架的分布空间。

b.等效细尺度网络中需要生成的微孔隙数目由下式求得，

N'_p＝ξ³N_p

其中N_p为原始细尺度网络中的孔隙数目，ξ为尺度变化因子。

c.计算原细尺度孔隙网络中孔喉参数的相关性，两个随机变量S和T间的相关性用相关系数ρ(S,T)是两变量线性关系的度量。对于非线性关系，可以转换为线性关系，对于S,T的序列{(s₁,t₁),...,(s_n,t_n)}，ρ定义为：

相关系数绝对值大表示相关性强。通过计算相关系数，可以确定各参数之间的从属关系。

d.根据原细尺度网络的连通函数计算需要充填的喉道数目。

描述网络拓扑特征的连通性函数可通过简化的比欧拉示性数来描述孔隙的全局连通性由连通性函数来表示：

式中，N_p(r)为原始孔隙网络孔隙集合中半径大于r的孔隙的数目，N_t(r)为原始孔隙网络喉道集合中半径大于r的喉道数目，V为原始孔隙网络的体积。根据欧拉数定义可得需要填充的喉道数目N_t'(r_k)可表示为

N_t'(r_k)＝N'_p(r_k)-χ_v(r_k)V'

式中，N'_p(r_k)表示从新建孔隙集合p'中孔隙半径大于r_k的孔隙数目，V'为新建网络的体积。根据配位数分布，按孔隙半径从大到小依次添加孔隙之间的连通喉道，并保证添加喉道的几何特征同原大孔隙网络一致，进而构建出同原始孔隙网络几何拓扑等效的随机孔隙网络模型。

e.利用微孔隙比例因子提高建模速度。

引入微孔隙比例因子的概念，微孔隙比例因子是指在构建多尺度孔隙网络模型过程中实际填充的微孔隙数目N”_pm占通过尺度变化因子计算所得微孔隙数目N'_pm的比例，定义如下：

α＝N”_pm/N'_pm×100％

结合微孔隙比例因子的定义，根据原始细尺度孔隙网络模型中孔隙的几何参数变量之间的概率分布函数和变量之间的相关性信息，可以在细尺度孔隙网络模型空间域中生成同原始细尺度孔隙网络几何拓扑结构信息等效的微孔隙，并根据统计分析随机产生各微孔隙单元的半径)体积)形状因子等，使各变量的概率分布与原始细尺度孔隙网络模型一样。

4)通过孔隙网络跨尺度耦合实现粗尺度孔隙网络和等效细尺度孔隙网络的叠加。

粗尺度孔隙网络和新生成的细尺度孔隙网络之间是相互不连通的，因此需要添加大孔隙和微孔隙之间的连通喉道来进行两个孔隙网络的整合。具体计算思路如下：在细尺度孔隙网络中心，放置粗尺度孔隙网络中半径最小的孔隙体，计算此孔隙体表面与周围微孔隙喉道的配位数，不断增加孔隙体的半径进而可得到一系列孔隙体表面和周围微孔隙喉道的配位数，进而得到大孔隙和微孔隙之间的配位数分布图，根据大孔隙和微孔隙之间的配位数分布图来添加大孔隙与其相邻微孔隙之间的连接喉道。

a.确定跨尺度连接的最大孔隙间距

对于一个喉道连接的两个孔隙，半径为分别为r₁,r₂，定义孔隙间距：

d_i(N₁,N₂)＝d(N₁,N₂)-r₁-r₂

其中d(N₁,N₂)为两个孔隙中心点的距离。遍历原始细尺度孔隙网络的孔隙间距，将max(d_i(N₁,N₂),i＝1,...,n)作为跨尺度连接的孔隙间距的最大值。

b.确定跨尺度配位数

跨尺度配位数定义为一个粗尺度孔隙所能连接的细尺度孔隙的最大数目。对于原始细尺度孔隙网络，其一系列的孔隙{N₁,N₂,…,N_n}和喉道{B₁,B₂,…,B_m}，在区域的中心插入s个同心圆球，同时保证最大的球仍在孔隙网络内。计算数目N_Bj,j＝1,…,s,即穿透相应圆球的喉道数目。也即，如果d(p,N)为同心圆球p和孔隙N间的距离，N_k1和N_k2为喉道B_k连接的两孔隙，并且d(p,N_k1)＜d(p,N_k2)，则：N_Bj＝{B_k|d(p,N_k1)＜r_j＜d(p,N_k2),k＝1,2…,m}。定义跨尺度系数γ_CSC

其中a_j是第j个球的表面积。跨尺度系数表示穿透单位孔隙面积的喉道数。

为了确定交叉尺度配位数，考虑粗尺度孔隙网络中的一个孔隙，它的一部分表面积已被粗尺度的喉道占据。对于每一个粗尺度的孔隙N_c，对应半径为r_c，粗尺度配位数z_c，其交叉尺度配位数的计算公式为：

其中，A_k为半径为r_k的粗尺度喉道在球面上所占的面积，

实施例：

1.选取四川油田某井岩心为研究对象，首先将待扫描的岩样进行洗油、洗盐处理，目的是去除结晶在岩心喉道里面的盐，然后进行常规物性测量，得到岩心孔隙度和渗透率等参数。

2.岩样的X射线CT扫描

a、仪器预热、载物台位置调节以及X射线参数的设置

打开仪器，对仪器做抽真空预热处理，标准是真空度<10nba,初始化100％完成。v|tome|xs 180的载物台可沿上、下、前、后、左和右六个方向调节。根据X射线CT扫描原理，在扫描过程中载物台要旋转180°或360°。因此，在扫描样品前需通过调整载物台的位置，使样品在载物台旋转过程中不得超出探测器的检测范围，否则无法成像。设置X射线管的电压、电流等参数，依据是保证X射线能穿透样品。

b、选择高低两种分辨率对岩心进行扫描

扫描分辨率的选择是利用CT扫描重构岩心三维图像关键的一步，直接关系在三维图像中能否准确直观反映岩心的孔隙结构。本发明中分别通过比对三维CT图像孔隙度、渗透率与实验孔隙度、渗透率比对，进行扫描分辨率的优选。由于碳酸盐岩的非均质性比较强，首先对柱塞岩样(直径为2.5cm)进行扫描，分辨率为14.4μm,然后对未能识别的孔隙部分钻取小样(直径为5mm)进行高分辨率扫描，分辨率为3.4μm。构建结果如图1所示，其中红色标记部分就是低分辨率未能完全识别的微孔区域。

3.基于三维数字岩心利用最大球算法提取岩心的三维孔隙空间网络，分别得到粗尺度和细尺度的孔隙网络模型，如图2、3和4所示。

4.将利用最大球算法提取的粗尺度孔隙网络模型直接移植过来作为新模型中的粗尺度模型；结合细尺度网络的几何拓扑结构信息随机重建等效细尺度孔隙网络，具体步骤如下：

a.确定细尺度孔隙网络填充的空间域。

在灰度图上能被完全识别孔隙的概率设定为1，对应最大灰度值为T₁，骨架为孔隙的概率设定为0，对应最小灰度值为T₂，微孔隙的分布范围满足灰度值T属于区间[T₁,T₂]，如图5，通过下式可以求出每个微孔灰度值为孔隙的概率。

假设多尺度孔隙网络的孔隙空间为Ω，粗尺度孔隙网络孔隙空间为Ω₁，细尺度网络的填充空间域为：

Ω₂＝Ω-Ω₁-Ω₃

其中Ω₃为骨架的分布空间。

b.等效细尺度网络中需要生成的微孔隙数目由下式求得，

N'_p＝ξ³N_p

其中N_p为原始细尺度网络中的孔隙数目，ξ为尺度变化因子，可以通过原细尺度网络的尺寸和等效细尺度的尺寸相比得到。

c.计算原细尺度孔隙网络中孔喉参数的相关性，两个随机变量S和T间的相关性用相关系数ρ(S,T)是两变量线性关系的度量。对于非线性关系，可以转换为线性关系，对于S,T的序列{(s₁,t₁),...,(s_n,t_n)}，ρ定义为：

相关系数绝对值大表示相关性强。通过计算相关系数，可以确定各参数之间的从属关系。

d.根据原细尺度网络的连通函数计算需要充填的喉道数目。

5.通过孔隙网络跨尺度耦合实现粗尺度孔隙网络和等效细尺度孔隙网络的叠加。在细尺度孔隙网络中心，放置粗尺度孔隙网络中半径最小的孔隙体，计算此孔隙体表面与周围微孔隙喉道的配位数，不断增加孔隙体的半径进而可得到一系列孔隙体表面和周围微孔隙喉道的配位数，进而得到大孔隙和微孔隙之间的配位数分布图，根据大孔隙和微孔隙之间的配位数分布图来添加大孔隙与其相邻微孔隙之间的连接喉道。在确定了最大孔隙间距和跨尺度配位数后，对于给定粗尺度孔隙N_c，标注其邻近细尺度孔隙N_f，计算孔隙间距d_i(N_f,N_c)，统计小于最大孔隙间距d_i,max的细尺度孔隙形成集合{N_f}。计算跨尺度配位数z_cs(N_c)，实际的跨尺度配位数为z_cs(N_c)＝min{z_cs(N_c),N(N_f)}。从细尺度孔隙集合{N_f}内依次选择，生成跨尺度的喉道，根据之前的统计信息产生半径，体积，长度，形状因子等参数。通过上述步骤，可将与原始粗尺度孔隙网络模型和细尺度孔隙网络模型几何拓扑结构信息等效的粗尺度孔隙网络和细尺度孔隙网络进行整合，如图6，构建出同时包含两种尺度孔隙网络模型结构信息的多尺度孔隙网络模型，如图7。

6.基于构建的多尺度网络计算了岩心孔隙度、渗透率、胶结指数以及饱和度指数等参数，并与实验结果进行了比较，如表1所示；结果表明模拟结果和实验结果基本吻合，相对误差在10％左右。

表1计算结果比较

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。