CN109285222A - 有机页岩高分辨率数字岩心构建与分析方法 - Google Patents
有机页岩高分辨率数字岩心构建与分析方法 Download PDFInfo
- Publication number
- CN109285222A CN109285222A CN201811035256.0A CN201811035256A CN109285222A CN 109285222 A CN109285222 A CN 109285222A CN 201811035256 A CN201811035256 A CN 201811035256A CN 109285222 A CN109285222 A CN 109285222A
- Authority
- CN
- China
- Prior art keywords
- image
- value
- shale
- threshold
- gray
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 239000011435 rock Substances 0.000 title claims abstract description 67
- 238000004458 analytical method Methods 0.000 title claims abstract description 27
- 238000000034 method Methods 0.000 claims abstract description 34
- 230000011218 segmentation Effects 0.000 claims abstract description 26
- 238000002591 computed tomography Methods 0.000 claims abstract description 15
- 238000003384 imaging method Methods 0.000 claims abstract description 14
- 230000005469 synchrotron radiation Effects 0.000 claims abstract description 11
- 229910052500 inorganic mineral Inorganic materials 0.000 claims abstract description 8
- 239000011707 mineral Substances 0.000 claims abstract description 8
- 238000002360 preparation method Methods 0.000 claims abstract description 5
- 238000012545 processing Methods 0.000 claims abstract description 4
- 238000001914 filtration Methods 0.000 claims description 24
- 238000005516 engineering process Methods 0.000 claims description 22
- 239000011159 matrix material Substances 0.000 claims description 9
- 230000000007 visual effect Effects 0.000 claims description 9
- 239000011148 porous material Substances 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 230000001174 ascending effect Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 4
- 235000002918 Fraxinus excelsior Nutrition 0.000 claims 1
- 239000002956 ash Substances 0.000 claims 1
- 210000000988 bone and bone Anatomy 0.000 claims 1
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000002474 experimental method Methods 0.000 description 9
- 238000011160 research Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 229910052683 pyrite Inorganic materials 0.000 description 4
- 239000011028 pyrite Substances 0.000 description 4
- NIFIFKQPDTWWGU-UHFFFAOYSA-N pyrite Chemical compound [Fe+2].[S-][S-] NIFIFKQPDTWWGU-UHFFFAOYSA-N 0.000 description 4
- 238000005481 NMR spectroscopy Methods 0.000 description 3
- 239000002734 clay mineral Substances 0.000 description 2
- 239000000470 constituent Substances 0.000 description 2
- 230000005611 electricity Effects 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 239000003079 shale oil Substances 0.000 description 2
- 238000002922 simulated annealing Methods 0.000 description 2
- 238000009736 wetting Methods 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000010884 ion-beam technique Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003129 oil well Substances 0.000 description 1
- 239000010453 quartz Substances 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N silicon dioxide Inorganic materials O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T19/00—Manipulating 3D models or images for computer graphics
- G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration by the use of local operators
-
- G06T5/70—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10056—Microscopic image
- G06T2207/10061—Microscopic image from scanning electron microscope
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
- G06T2207/20032—Median filtering
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Architecture (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本发明公开了一种有机页岩高分辨率数字岩心构建与分析方法，其特征在于：有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：A.页岩样品制备与同步辐射扫描；B.对页岩岩心灰度图像进行阈值滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于CT扫描重建过程中产生的噪声点；C.对CT图像进行阈值分割，划分岩石的骨架和孔隙。本发明的有益效果是：完整介绍了一种有机页岩高分辨率成像与数字岩心构建的方法与流程。在骨架与孔隙的阈值求取中，采用两组间性质差异与组内差异之比达到最大值确定阈值，准确实现了孔隙组分、矿物、干酪根等组分的分辨。
Description
技术领域
本发明涉及石油测井岩石物理研究领域，具体地说，涉及一种有机页岩高分辨率数字岩心构建与分析方法。
背景技术
页岩油气是一种重要的非常规油气类型，针对有机页岩岩石物理特性研究是储层测井评价与地震预测的基础。岩石物理研究主要包括理论方法研究，岩石物理实验以及近年来发展起来的数字岩石物理分析技术。由于页岩孔隙小，矿物成分复杂，数字岩心分析及其数值模拟技术是一个有效的分析手段。因此，如何结合有机页岩的特点构建有效的接近实际的数字岩心模型则是微观数值模拟的基础与关键。
目前，数字岩心建模方法有两大类：物理实验方法和数值重建方法。物理实验方法借助扫描电镜和CT成像仪等高精度仪器获取岩心二维图像，经三维重建即可得到三维数字岩心(Hazlett et al.,1995)。近年来，更高分辨率的FIB-SEM(Focused Ion BeamCombined with SEM)图像扫描技术，其分辨率高达3nm，能够满足页岩的岩石物理实验。数值重建方法主要是基于岩心切片和重建算法建立数字岩心，如模拟退火法(王晨晨等，2012)、过程法(Jin等，2011)、多点地质统计学以及一些混合方法，如过程法和模拟退火法(刘学峰等，2010)等。对于有机页岩大多采用FIB-SEM方法。
对于有机页岩，目前采用FIB-SEM成像技术比较昂贵，限制了该项技术的使用。而且，该技术不只针对页岩，实验的结果是针对具体的岩心，不能代表形形色色的岩石，不能反映不同润湿性的情况，也不能反映不同饱和度的情形，更不能模拟不同含气饱和度的情形。因此，目前的技术进展限制了后续微观数值模拟中影响因素和微观响应规律的研究。
发明内容
本发明正是为了解决上述技术问题而设计的一种针对页岩油气的数字岩心模型的制作方法，能够用于模拟不同饱和度，不同润湿性的含气或含油页岩的微观数值模拟，能够为后续的电学、声学、核磁共振等物理特性的微观数值模拟提供基础模型。
本发明解决其技术问题所采用的技术方案是：
一种有机页岩高分辨率数字岩心构建与分析方法，有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：
A.页岩样品制备与同步辐射扫描
a1.借助高分辨率电子显微镜选取页岩样品，放入CT仪器载物台进行CT扫描；
a2.扫描图像重构。从扫描图像中剔除背景图像，得到X-Y平面的单帧16位页岩岩心灰度图像；
a3.用图像处理软件把多帧16位图像构建为三维图像，观察岩心样本的整体结构特征和孔隙特征；
B.数字图像滤波。对步骤a2所得页岩岩心灰度图像进行图像滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于CT扫描重建过程中产生的噪声点；
C.CT图像组分分割。对步骤B所得CT图像进行阈值分割，划分岩石的骨架和孔隙：
c1.阈值计算及骨架与孔隙的准确分割
采用阈值分割法，其原理是：岩心CT图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值，通过阈值选取技术选取合理的阈值，将CT图像中的骨架和孔隙准确地分割开来；
c2.骨架矿物的精细分割
对于一般岩石CT图像，只要准确地分割CT图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩CT图像来说，由于其成分复杂，还需要对页岩CT图像的矿物成分进行进一步分割，其主要步骤为：
第一、导入CT图像数据，设置数据格式，将三维CT图像裁剪为需要的尺寸；
第二、依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；
第三、使用灰度阈值分割程序，精细地划分其矿物成分；
第四、导出数据并构建数字岩心，即岩心三维成像图。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤A中同步辐射扫描是指对页岩岩心的高分辨率扫描，其高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm，大视场模式成像视场大小为60μm×60μm。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤B中阈值滤波是指采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理：
二维中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑；图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:
进行滤波处理时，常常选用不同窗口大小[用N表示]，窗口大小的数值为奇数，优先采用窗口大小为3×3的正方形窗口作为中值滤波的模板，其模板定义公式如下：
式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值；具体计算步骤如下:
第一、将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；
第二、读取出二维窗口中各个像素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值，即中值；
第三、用计算得到的中值替代中值滤波窗口中心点的灰度值。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤C中的阈值选取技术包括最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值技术。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤C中的阈值选取技术采用最大类间方差阈值技术，采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过选取不同的阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值；
设Imax，Imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为N1(k)、N2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1 2(k)、σ2 2(k)，全部像素的平均灰度值为μ1，则组内方差σw2(k)：
σw 2(k)＝N1σ1 2+N2σ2 2 [3]
组件方差为：
σB 2(k)＝N1(μ1-μT)2+N2(μ2+μT)2 [4]
则满足的下式的k*即为所求的阈值，
f(k*)＝max{f(k)＝σB 2(k)/σw 2(k)} [5]
式中，
本发明的有益效果是：本发明完整介绍了有机页岩高分辨率成像与数字岩心构建的方法与流程。在骨架与孔隙的阈值求取中，采用两组间性质差异与组内差异之比达到最大值确定阈值，准确实现了孔隙组分、矿物、干酪根等组分的分辨。该方法得到的数字岩心已经在核磁共振、声学孔隙尺度数值模拟与微观特性得到了成功应用。为孔隙结构研究以及后续的微观电学、声学、核磁共振、渗流等岩石物理特性提供了条件。
附图说明
图1为页岩岩心X射线CT扫描单帧灰度图像。
图2为页岩岩心X射线CT扫描3D图像。
图3为页岩岩心X射线CT扫描图未滤波的图像。
图4为页岩岩心X射线CT扫描图滤波后的图像。
图5为页岩CT图像阈值计算界面图。
图6为页岩岩心矿物成分、孔隙流体的分割图像。
图7为重构的页岩三维数字岩心图像。
图8为另一块页岩三维数字岩心图像。
图9为页岩孔隙空间。
图10为页岩干酪根。
图11为页岩岩石骨架。
图12为页岩黄铁矿。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
如图1-12所示，本发明一种有机页岩高分辨率数字岩心构建与分析方法，有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：
A.页岩样品制备与同步辐射扫描
a1.借助高分辨率电子显微镜选取页岩样品，放入CT仪器载物台进行CT扫描；
a2.扫描图像重构。从扫描图像中剔除背景图像，得到X-Y平面的单帧16位页岩岩心灰度图像；
a3.用图像处理软件把多帧16位图像构建为三维图像，观察岩心样本的整体结构特征和孔隙特征；
B.数字图像滤波。对步骤a2所得页岩岩心灰度图像进行图像滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于CT扫描重建过程中产生的噪声点；
C.CT图像组分分割。对步骤B所得CT图像进行阈值分割，划分岩石的骨架和孔隙：
c1.阈值计算及骨架与孔隙的准确分割
采用阈值分割法，其原理是：岩心CT图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值，通过阈值选取技术选取合理的阈值，将CT图像中的骨架和孔隙准确地分割开来；
c2.骨架矿物的精细分割
对于一般岩石CT图像，只要准确地分割CT图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩CT图像来说，由于其成分复杂，还需要对页岩CT图像的矿物成分进行进一步分割，其主要步骤为：
第一、导入CT图像数据，设置数据格式，将三维CT图像裁剪为需要的尺寸；
第二、依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；
第三、使用灰度阈值分割程序，精细地划分其矿物成分；
第四、导出数据并构建数字岩心，即岩心三维成像图。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤A中同步辐射扫描是指对页岩岩心的高分辨率扫描，其高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm，大视场模式成像视场大小为60μm×60μm。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤B中阈值滤波是指采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理：
二维中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑；图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:
进行滤波处理时，常常选用不同窗口大小[用N表示]，窗口大小的数值为奇数，优先采用窗口大小为3×3的正方形窗口作为中值滤波的模板，其模板定义公式如下：
式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值；具体计算步骤如下:
第一、将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；
第二、读取出二维窗口中各个像素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值，即中值；
第三、用计算得到的中值替代中值滤波窗口中心点的灰度值。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤C中的阈值选取技术包括最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值技术。
所述有机页岩高分辨率数字岩心构建与分析方法，其步骤C中的阈值选取技术采用最大类间方差阈值技术，采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过选取不同的阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值；
设Imax，Imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为N1(k)、N2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1 2(k)、σ2 2(k)，全部像素的平均灰度值为μ1，则组内方差σw 2(k)：
σw 2(k)＝N1σ1 2+N2σ2 2 [3]
组件方差为：
σB 2(k)＝N1(μ1-μT)2+N2(μ2+μr)2 [4]
则满足的下式的k*即为所求的阈值，
f(k*)＝max{f(k)＝σB 2(k)/σw 2(k)} [5]
式中，
实施例：
1.页岩样品制备与同步辐射扫描
页岩岩心高分辨率扫描是在中国科学院高能物理研究所北京同步辐射装置试验站下的4W1A-X射线成像实验站进行。实验站所用光源是4W1A光束线是北京同步辐射装置最长的一条光束线，可以提供准直性较好的X射线光束，其中纳米成像设备可获得的空间分辨率最高可达30纳米。X射线纳米分辨三维成像实验可以提供高分辨模式和大视场模式这两种实验模式：高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm(标准样品)，大视场模式成像视场大小为60μm×60μm。因为页岩岩样比较致密，同时页岩对于X射线的吸收较强，采用高分辨模式进行CT图像扫描，主要步骤为：
(1)借助高分辨率电子显微镜选作页岩样品，放入CT仪器载物台进行CT扫描；
(2)扫描图像重构。从扫描图像中剔除背景图像，得到X-Y平面的单帧16位灰度图像，如图1；
(3)用TXM3DViewer软件把多帧16位图像构建为三维图像，以展示其整体特征及形态，如图2。
可以看到，页岩样本的CT扫描图像中小孔隙占整个样本的大多数，多为独立的单独孔隙。
2.阈值滤波
三维页岩岩样切片图像中含有一定的噪声(图3)，需要对图像进行滤波处理。采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理。中值滤波器是一种非线性的滤波技术，它是基于对于信号数值的排序统计理论，对信号进行重新恢复的技术，它可以很好地抑制噪声，同时保留信号中真实的信息和图像的边缘特征。
图像中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑。图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:
进行滤波处理时，常常选用不同窗口大小(用N表示)，窗口大小的数值为奇数，采用窗口大小为3*3的正方形窗口作为中值滤波的模板，其模板定义公式如下：
式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值。具体计算步骤如下:
(1)将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；
(2)读取出二维窗口中各个体素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值(中值)；
(3)用计算得到的中值替代中值滤波窗口中心点的灰度值。
通过编写二维中值滤波程序，对岩心CT图像进行滤波处理，其对比效果如图4所示。可以看出，经过中值滤波处理过后，岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于CT扫描重建过程中产生的噪声点，使得孔隙和骨架的灰度值分布更加明显。
3.CT图像组分分割
岩石CT图像的矿物组分分割是生成三维数字岩心的非常重要的一步，对CT图像进行阈值分割，划分岩石的骨架和孔隙。
1)骨架与孔隙的分割与阈值计算
最常用的方法就是阈值分割法。其原理是，岩心CT图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值。
如何选取合理的阈值就是三维数字岩心能否符合实际岩心的关键。常见的阈值选取技术主要有：最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值。采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过不断选取阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值。
设Imax，Imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为N1(k)、N2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1 2(k)、σ2 2(k)，全部像素的平均灰度值为μ1，则组内方差σw 2(k)：
σw 2(k)＝N1σ1 2+N2σ2 2 (3)
组件方差为：
σB 2(k)＝N1(μ1-μT)2+N2(μ2+μT)2 (4)
则满足的下式的k*即为所求的阈值，
f(k*)＝max{f(k)＝σB 2(k)/σw 2(k)} (5)式中，
通过所求的阈值，就可以将CT图像中的骨架和孔隙准确地分割开来。
根据上述原理，编写了CT图像阈值求取程序，通过该程序可以实现CT图像不同大小的裁剪，对CT图像再用最大类间方差法求取阈值，并将该代码编写为GUI界面，如图5所示。
2)骨架矿物的精细分割
对于一般岩石CT图像，只要准确的分割CT图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩CT图像来说，由于其成分复杂，还需要对页岩CT图像的矿物成分进行进一步分割。进行CT图像处理的主要步骤为：
(1)导入CT图像数据，设置数据格式，将三维CT图像裁剪为需要的尺寸；
(2)依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；
(3)使用灰度阈值分割程序，精细的划分其矿物成分，骨架为红色，孔隙为蓝色，黄色为干酪根，浅绿色为黄铁矿，如图6所示；
(4)导出数据并构建数字岩心,绘制岩心三维成像图，如图7所示。
由X射线CT扫描原理可知，如果岩石组分密度较高，则在灰度图像上显示为高亮度，其灰度值也低，如黄铁矿；而密度较低的组分，在灰度图像上显示为低亮度，其灰度值也高，如孔隙内的流体；密度介于两者之间的组分，则在灰度图像上表现为灰白色，例如脆性矿物和黏土矿物。
图8为另一块有机页岩的数字岩心构建实例。从同步辐射图像扫描图像，到纳米成像、数字岩心构建以及组分分析。图中最终处理结果为：黏土矿物为92.2％，石英为3.1％，两者为骨架如图9所示；孔隙约2.4％(图10)，干酪根为2.0％(图11)，黄铁矿为0.3％(图12)。
本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。
Claims (5)
1.一种有机页岩高分辨率数字岩心构建与分析方法，其特征在于：有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：
A.页岩样品制备与同步辐射扫描
a1.借助高分辨率电子显微镜选取页岩样品，放入CT仪器载物台进行CT扫描；
a2.扫描图像重构；从扫描图像中剔除背景图像，得到X-Y平面的单帧16位页岩岩心灰度图像；
a3.用图像处理软件把多帧16位图像构建为三维图像，观察岩心样本的整体结构特征和孔隙特征；
B.数字图像滤波；对步骤a2所得页岩岩心灰度图像进行图像滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于CT扫描重建过程中产生的噪声点；
C.CT图像组分分割；对步骤B所得CT图像进行阈值分割，划分岩石的骨架和孔隙：
c1.阈值计算及骨架与孔隙的准确分割
采用阈值分割法，其原理是：岩心CT图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值，通过阈值选取技术选取合理的阈值，将CT图像中的骨架和孔隙准确地分割开来；
c2.骨架矿物的精细分割
对于一般岩石CT图像，只要准确地分割CT图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩CT图像来说，由于其成分复杂，还需要对页岩CT图像的矿物成分进行进一步分割，其主要步骤为：
第一、导入CT图像数据，设置数据格式，将三维CT图像裁剪为需要的尺寸；
第二、依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；
第三、使用灰度阈值分割程序，精细地划分其矿物成分；
第四、导出数据并构建数字岩心，即岩心三维成像图。
2.根据权利要求1所述的有机页岩高分辨率数字岩心构建与分析方法，其特征在于：步骤A中同步辐射扫描是指对页岩岩心的高分辨率扫描，其高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm，大视场模式成像视场大小为60μm×60μm。
3.根据权利要求1所述的有机页岩高分辨率数字岩心构建与分析方法，其特征在于：步骤B中阈值滤波是指采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理：
二维中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑；图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:
进行滤波处理时，常常选用不同窗口大小[用N表示]，窗口大小的数值为奇数，优先采用窗口大小为3×3的正方形窗口作为中值滤波的模板，其模板定义公式如下：
式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值；具体计算步骤如下:
第一、将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；
第二、读取出二维窗口中各个像素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值，即中值；
第三、用计算得到的中值替代中值滤波窗口中心点的灰度值。
4.根据权利要求1所述的有机页岩高分辨率数字岩心构建与分析方法，其特征在于：步骤C中的阈值选取技术包括最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值技术。
5.根据权利要求1所述的有机页岩高分辨率数字岩心构建与分析方法，其特征在于：步骤C中的阈值选取技术采用最大类间方差阈值技术，采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过选取不同的阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值；
设Imax，Imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为N1(k)、N2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1 2(k)、σ2 2(k)，全部像素的平均灰度值为μ1，则组内方差σw 2(k):
σw 2(k)＝N1σ1 2+N2σ2 2 [3]
组件方差为：
σB 2(k)＝N1(μ1-μT)2+N2(μ2+μT)2 [4]
则满足的下式的k*即为所求的阈值，
f(k*)＝max{f(k)＝σB 2(k)/σB 2(k)/σw 2(k)} [5]
式中，
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811035256.0A CN109285222A (zh) | 2018-09-06 | 2018-09-06 | 有机页岩高分辨率数字岩心构建与分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811035256.0A CN109285222A (zh) | 2018-09-06 | 2018-09-06 | 有机页岩高分辨率数字岩心构建与分析方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109285222A true CN109285222A (zh) | 2019-01-29 |
Family
ID=65184175
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811035256.0A Pending CN109285222A (zh) | 2018-09-06 | 2018-09-06 | 有机页岩高分辨率数字岩心构建与分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109285222A (zh) |
Cited By (23)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110095388A (zh) * | 2019-04-18 | 2019-08-06 | 中国石油大学(北京) | 碎屑岩颗粒结构的确定方法及装置 |
| CN110441204A (zh) * | 2019-07-09 | 2019-11-12 | 大庆油田有限责任公司 | 一种基于数字岩心模拟的致密储层压裂液伤害数字化评价方法 |
| CN110599473A (zh) * | 2019-09-06 | 2019-12-20 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| CN110873724A (zh) * | 2018-09-03 | 2020-03-10 | 中国石油化工股份有限公司 | 一种页岩岩心图像采集方法 |
| CN110910411A (zh) * | 2019-11-14 | 2020-03-24 | 中国科学院力学研究所 | 一种具有尺寸自适应的页岩裂缝自动提取方法 |
| CN110927194A (zh) * | 2019-12-11 | 2020-03-27 | 中国科学院地质与地球物理研究所 | 确定泥页岩有机孔含量和孔径分布的方法 |
| CN111179296A (zh) * | 2020-01-08 | 2020-05-19 | 青岛科技大学 | 一种基于数字岩心技术研究岩石导热特性的新方法 |
| CN111680440A (zh) * | 2020-06-03 | 2020-09-18 | 中国地质大学（北京） | 不同温度下数字岩心声学特性分析方法 |
| CN112070774A (zh) * | 2020-09-16 | 2020-12-11 | 西南石油大学 | 一种用于页岩数字岩心图像分割的神经网络优化方法 |
| CN112131776A (zh) * | 2020-08-26 | 2020-12-25 | 中国地质大学（北京） | 裂缝性三维数字岩心油水或气水分布模型的构建方法 |
| CN112132965A (zh) * | 2020-09-25 | 2020-12-25 | 中国矿业大学 | 一种岩土体孔裂隙结构多尺度表征方法 |
| CN112801205A (zh) * | 2021-02-20 | 2021-05-14 | 中国科学院地质与地球物理研究所 | 岩石纹层结构的获取方法及装置 |
| CN112816388A (zh) * | 2020-12-31 | 2021-05-18 | 中国石油大学(北京) | 基于ct和数字岩心三维重构的油砂渗流性能测试方法 |
| CN113515847A (zh) * | 2021-05-12 | 2021-10-19 | 中国矿业大学 | 一种基于K-means聚类算法的非均质岩石数字岩芯建模方法 |
| CN113607618A (zh) * | 2021-07-02 | 2021-11-05 | 中铁二十局集团有限公司 | 基于工业ct扫描的岩石损伤分析方法、装置及设备 |
| CN113609696A (zh) * | 2021-08-16 | 2021-11-05 | 中国地质大学（北京） | 基于图像融合的多尺度多组分数字岩心构建方法及系统 |
| CN113916916A (zh) * | 2021-09-29 | 2022-01-11 | 西南石油大学 | 一种用于页岩数字岩心三维渗流-颗粒流动耦合的模拟方法 |
| CN114092489A (zh) * | 2021-11-02 | 2022-02-25 | 清华大学 | 多孔介质渗流通道提取、模型训练方法、装置及设备 |
| CN114518309A (zh) * | 2020-11-20 | 2022-05-20 | 中国石油化工股份有限公司 | 一种页岩油储层中夹层孔隙有效三维表征的方法及装置 |
| CN116402675A (zh) * | 2023-03-23 | 2023-07-07 | 中国地质科学院地质力学研究所 | 一种基于页岩组分标定的图像配准方法 |
| CN116858630A (zh) * | 2023-06-14 | 2023-10-10 | 中国石油大学(华东) | 一种真实页岩干酪根孔隙模型的构建方法及系统 |
| CN117036635A (zh) * | 2023-10-10 | 2023-11-10 | 中国石油大学(华东) | 一种基于图像纹理分类的页岩多尺度数字岩心构建方法 |
| CN113607618B (zh) * | 2021-07-02 | 2024-05-10 | 中铁二十局集团有限公司 | 基于工业ct扫描的岩石损伤分析方法、装置及设备 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105115874A (zh) * | 2015-08-18 | 2015-12-02 | 中国石油天然气股份有限公司 | 多源信息融合的多组分三维数字岩心构建方法 |
| CN105487121A (zh) * | 2015-12-03 | 2016-04-13 | 长江大学 | 基于ct扫描图像与电成像图像融合构建多尺度数字岩心方法 |
| CN105551018A (zh) * | 2015-12-03 | 2016-05-04 | 长江大学 | 面向碳酸盐岩复杂储层数字岩心的各组分提取方法 |
| CN106127816A (zh) * | 2016-03-08 | 2016-11-16 | 中国石油大学(华东) | 一种页岩基质储层孔隙空间表征方法 |
-
2018
- 2018-09-06 CN CN201811035256.0A patent/CN109285222A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105115874A (zh) * | 2015-08-18 | 2015-12-02 | 中国石油天然气股份有限公司 | 多源信息融合的多组分三维数字岩心构建方法 |
| CN105487121A (zh) * | 2015-12-03 | 2016-04-13 | 长江大学 | 基于ct扫描图像与电成像图像融合构建多尺度数字岩心方法 |
| CN105551018A (zh) * | 2015-12-03 | 2016-05-04 | 长江大学 | 面向碳酸盐岩复杂储层数字岩心的各组分提取方法 |
| CN106127816A (zh) * | 2016-03-08 | 2016-11-16 | 中国石油大学(华东) | 一种页岩基质储层孔隙空间表征方法 |
Non-Patent Citations (1)
| Title |
|---|
| NOBUYUKI OSTU: "A Tlreshold Selection Method from Dray-level Histograms", 《IEEE》 * |
Cited By (37)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110873724A (zh) * | 2018-09-03 | 2020-03-10 | 中国石油化工股份有限公司 | 一种页岩岩心图像采集方法 |
| CN110095388A (zh) * | 2019-04-18 | 2019-08-06 | 中国石油大学(北京) | 碎屑岩颗粒结构的确定方法及装置 |
| CN110441204A (zh) * | 2019-07-09 | 2019-11-12 | 大庆油田有限责任公司 | 一种基于数字岩心模拟的致密储层压裂液伤害数字化评价方法 |
| CN110441204B (zh) * | 2019-07-09 | 2022-06-17 | 大庆油田有限责任公司 | 一种基于数字岩心模拟的致密储层压裂液伤害数字化评价方法 |
| CN110599473B (zh) * | 2019-09-06 | 2022-04-15 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| CN110599473A (zh) * | 2019-09-06 | 2019-12-20 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| CN110910411A (zh) * | 2019-11-14 | 2020-03-24 | 中国科学院力学研究所 | 一种具有尺寸自适应的页岩裂缝自动提取方法 |
| CN110910411B (zh) * | 2019-11-14 | 2022-12-09 | 中国科学院力学研究所 | 一种具有尺寸自适应的页岩裂缝自动提取方法 |
| CN110927194A (zh) * | 2019-12-11 | 2020-03-27 | 中国科学院地质与地球物理研究所 | 确定泥页岩有机孔含量和孔径分布的方法 |
| CN111179296A (zh) * | 2020-01-08 | 2020-05-19 | 青岛科技大学 | 一种基于数字岩心技术研究岩石导热特性的新方法 |
| CN111680440A (zh) * | 2020-06-03 | 2020-09-18 | 中国地质大学（北京） | 不同温度下数字岩心声学特性分析方法 |
| CN111680440B (zh) * | 2020-06-03 | 2024-03-19 | 中国地质大学（北京） | 不同温度下数字岩心声学特性分析方法 |
| CN112131776A (zh) * | 2020-08-26 | 2020-12-25 | 中国地质大学（北京） | 裂缝性三维数字岩心油水或气水分布模型的构建方法 |
| CN112131776B (zh) * | 2020-08-26 | 2023-05-05 | 中国地质大学（北京） | 裂缝性三维数字岩心油水或气水分布模型的构建方法 |
| CN112070774A (zh) * | 2020-09-16 | 2020-12-11 | 西南石油大学 | 一种用于页岩数字岩心图像分割的神经网络优化方法 |
| CN112132965B (zh) * | 2020-09-25 | 2024-03-26 | 中国矿业大学 | 一种岩土体孔裂隙结构多尺度表征方法 |
| CN112132965A (zh) * | 2020-09-25 | 2020-12-25 | 中国矿业大学 | 一种岩土体孔裂隙结构多尺度表征方法 |
| CN114518309A (zh) * | 2020-11-20 | 2022-05-20 | 中国石油化工股份有限公司 | 一种页岩油储层中夹层孔隙有效三维表征的方法及装置 |
| CN112816388A (zh) * | 2020-12-31 | 2021-05-18 | 中国石油大学(北京) | 基于ct和数字岩心三维重构的油砂渗流性能测试方法 |
| CN112801205B (zh) * | 2021-02-20 | 2022-04-15 | 中国科学院地质与地球物理研究所 | 岩石纹层结构的获取方法及装置 |
| CN112801205A (zh) * | 2021-02-20 | 2021-05-14 | 中国科学院地质与地球物理研究所 | 岩石纹层结构的获取方法及装置 |
| CN113515847B (zh) * | 2021-05-12 | 2023-09-05 | 中国矿业大学 | 一种基于K-means聚类算法的非均质岩石数字岩芯建模方法 |
| CN113515847A (zh) * | 2021-05-12 | 2021-10-19 | 中国矿业大学 | 一种基于K-means聚类算法的非均质岩石数字岩芯建模方法 |
| CN113607618A (zh) * | 2021-07-02 | 2021-11-05 | 中铁二十局集团有限公司 | 基于工业ct扫描的岩石损伤分析方法、装置及设备 |
| CN113607618B (zh) * | 2021-07-02 | 2024-05-10 | 中铁二十局集团有限公司 | 基于工业ct扫描的岩石损伤分析方法、装置及设备 |
| CN113609696B (zh) * | 2021-08-16 | 2024-03-08 | 中国地质大学（北京） | 基于图像融合的多尺度多组分数字岩心构建方法及系统 |
| CN113609696A (zh) * | 2021-08-16 | 2021-11-05 | 中国地质大学（北京） | 基于图像融合的多尺度多组分数字岩心构建方法及系统 |
| CN113916916B (zh) * | 2021-09-29 | 2024-03-19 | 西南石油大学 | 一种用于页岩数字岩心三维渗流-颗粒流动耦合的模拟方法 |
| CN113916916A (zh) * | 2021-09-29 | 2022-01-11 | 西南石油大学 | 一种用于页岩数字岩心三维渗流-颗粒流动耦合的模拟方法 |
| CN114092489B (zh) * | 2021-11-02 | 2023-08-29 | 清华大学 | 多孔介质渗流通道提取、模型训练方法、装置及设备 |
| CN114092489A (zh) * | 2021-11-02 | 2022-02-25 | 清华大学 | 多孔介质渗流通道提取、模型训练方法、装置及设备 |
| CN116402675B (zh) * | 2023-03-23 | 2023-11-28 | 中国地质科学院地质力学研究所 | 一种基于页岩组分标定的图像配准方法 |
| CN116402675A (zh) * | 2023-03-23 | 2023-07-07 | 中国地质科学院地质力学研究所 | 一种基于页岩组分标定的图像配准方法 |
| CN116858630A (zh) * | 2023-06-14 | 2023-10-10 | 中国石油大学(华东) | 一种真实页岩干酪根孔隙模型的构建方法及系统 |
| CN116858630B (zh) * | 2023-06-14 | 2024-04-09 | 中国石油大学(华东) | 一种真实页岩干酪根孔隙模型的构建方法及系统 |
| CN117036635A (zh) * | 2023-10-10 | 2023-11-10 | 中国石油大学(华东) | 一种基于图像纹理分类的页岩多尺度数字岩心构建方法 |
| CN117036635B (zh) * | 2023-10-10 | 2024-01-05 | 中国石油大学(华东) | 一种基于图像纹理分类的页岩多尺度数字岩心构建方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109285222A (zh) | 有机页岩高分辨率数字岩心构建与分析方法 | |
| CN105487121B (zh) | 基于ct扫描图像与电成像图像融合构建多尺度数字岩心方法 | |
| CN105115874B (zh) | 多源信息融合的多组分三维数字岩心构建方法 | |
| US8699769B2 (en) | Generating artificial hyperspectral images using correlated analysis of co-registered images | |
| US8909508B2 (en) | Petrographic image analysis for determining capillary pressure in porous media | |
| CN107709699B (zh) | 生成多孔岩石样本的三维微型模型 | |
| US20210321963A1 (en) | Systems and methods for enhanced imaging and analysis | |
| CN103698803A (zh) | 一种岩石孔隙结构表征方法及装置 | |
| DE102010039652A1 (de) | Mosaikaufnahmeerzeugung | |
| DE112019005143T5 (de) | System zur co-registrierung medizinischer bilder unter verwendung eines klassifikators | |
| CN106650811A (zh) | 一种基于近邻协同增强的高光谱混合像元分类方法 | |
| CN110264472B (zh) | 一种基于sd-oct的水果近皮下细胞无损成像方法 | |
| CN114565658A (zh) | 基于ct技术的孔隙尺寸计算方法及装置 | |
| CN113989343A (zh) | 基于注意力机制的传感器融合深度重建数据驱动方法 | |
| CN113075731B (zh) | 深层储层连续性井筒数字建模方法及装置 | |
| CN107622476B (zh) | 基于概率生成模型的图像超分辨处理方法 | |
| CN103778615B (zh) | 基于区域相似性的多聚焦图像融合方法 | |
| CN104732190B (zh) | 一种基于正交纹理相关分析的合成孔径声纳目标检测方法 | |
| DE112019004112T5 (de) | System und verfahren zur analyse mikroskopischer bilddaten und zur erzeugung eines annotierten datensatzes zum training der klassifikatoren | |
| CN113610863B (zh) | 多曝光图像融合质量评估方法 | |
| CN103955915A (zh) | 基于稀疏表示和多字典的sar图像分割 | |
| CN114705606A (zh) | 一种基于网络化分析的岩石内部关键渗流节点的封堵方法 | |
| KR102224466B1 (ko) | 기계학습을 이용한 페로브스카이트의 구조 분석 방법 | |
| Eitel et al. | Feature visualization for convolutional neural network models trained on neuroimaging data | |
| CN109325449A (zh) | 基于样本更新的卷积神经网络目标检测框架 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190129 |