数字岩⼼重建与孔隙参数分析_吴运强
第13卷第5期太赫兹科学与电⼦信息学报Vo1．13，No．5 2015年10⽉ Journal of Terahertz Science and Electronic Information Technology Oct．，2015
⽂章编号：2095-4980(2015)05-0788-07
数字岩⼼重建与孔隙参数分析
吴运强1，赵增义1，王⼦强1，左琛2，滕奇志2
(1.中⽯油新疆油⽥分公司实验检测研究院，新疆维吾尔⾃治区克拉玛依 834000；2.四川⼤学电⼦信息学院，四川成都610065)
摘要：针对岩⼼物性实验周期长且难以复现等问题，利⽤数字岩⼼图像分析技术进⾏岩⼼物性研究。以岩⼼铸体薄⽚图像为训练图像，根据模拟退⽕算法进⾏三维重建，获得岩⼼三维重
建图像。通过计算机断层扫描(CT)技术，获得同组岩⼼序列图像。在2种⽅式获得的图像中利⽤最
⼤球模型提取岩⼼孔隙结构，计算孔隙与喉道尺⼨及分布，将结果与压汞法测定的岩⼼⽑管分布
进⾏⽐较，分析了当前数字岩⼼技术与物性实验之间的区别与原因。结果表明，铸体薄⽚图像和
CT图像可以有效表征岩⼼⼤尺⼨孔隙特征，但物性实验⽅法可以测定微⼩孔隙尺⼨，由于成像设
备所限，数字岩⼼不能刻画该类孔隙特征。
关键词：三维重建；模拟退⽕算法；最⼤球模型；计算机断层扫描；岩⼼物性实验；岩⼼孔隙结构
中图分类号：TN911.73 ⽂献标识码：A doi：10.11805/TKYDA201505.0788
Reconstruction of digital core and analysis of pore parameters
WU Yunqiang1，ZHAO Zengyi1，WANG Ziqiang1，ZUO Chen2，TENG Qizhi2
(1.Institute of Xinjiang Oilfield Detection，PetroChina，Kelamayi Xinjiang 834000，China；
2.College of Electronics and Information Engineering，Sichuan University，Chengdu Sichuan 610065，China) Abstract：Digital core image analysis technique is utilized in the core physical study in order to solve the problems such as time-consuming and difficult to reproduce of the core physical properties
experiments. Taking the two-dimensional(2D) slices of the core as the training images, a three-
dimensional(3D) image is reconstructed by the simulated annealing algorithm. A 2D image quence of the
有理想same ries is acquired by Computed Tomography(CT) technique. The pore structures are extracted from
the two ts of images. Pore-throat sizes and their distribution are calculated by maximal balls model. By
comparing the results with that of the mercury intrusion method experiments, distinctions and reasons
between the current digital core analysis technique and traditional core physical properties experiments
are analyzed.Results show that the large pores could be depicted in 2D slices and CT images. However,
the core physical properties experiments are able to measure the size of micro pores, while the digital core reconstruction method could not due to the limitations of imaging equipments.
Key words：three-dimensional reconstruction；simulated annealing algorithm；maximal balls model；
Computed Tomography；core physical properties experiments；core pore structure
岩⼼物理特性实验，如压汞法测定⽑管压⼒实验、稳态法测定油⽔相对渗透率实验等，都是研究储集层物理特性，特别是渗流特性的有效⼿段，但物理实验周期长、过程繁琐，同时也存在着实验结果难以复现等问题。近年来，数字岩⼼成为了热门研究课题。以⼏张岩⼼⼆维图像为训练图像，通过数学建模⽅法构建三维图像，分析储集层岩⼼物理特性。数字岩⼼分析⽅法运算周期短，可以重复多次对同⼀块岩⼼进⾏实验，且实验结果可以进⾏动态显⽰。
企业金融模拟退⽕重建[1-3]是构建数字岩⼼的经典算法，本⽂将岩⼼铸体薄⽚图像作为训练图像，利⽤模拟退⽕算法重建岩⼼三维结构。同时利⽤计算机断层扫描技术(CT)，扫描同组岩⼼样本，获得序列图像。利⽤最⼤球模型[4]提取岩⼼⼆维图像和三维图像的拓扑结构，在链路的基础上划分孔隙与喉道，计算孔隙结构参数[5]。将得到的孔收稿⽇期：2014-10-21；修回⽇期：2014-12-09
基⾦项⽬：国家⾃然科学基⾦资助项⽬(No.61372174；60972130)
第5期吴运强等：数字岩⼼重建与孔隙参数分析 789
隙尺⼨分布状态与压汞法测定的⽑管束分布状态进⾏⽐较，分析了数字岩⼼分析⽅法与岩⼼物性实验的区别。
1 实验⽅法
1.1 压汞法测定⽑管压⼒曲线
储层岩⽯的⽑管压⼒和湿相饱和度的关系称为⽑管压⼒曲线，随着相关研究⼯作的开展，有关储层的⼏乎全部参数，如束缚⽔饱和度、残余油饱和度、孔隙度、孔喉⼤⼩分布等，在某种程度上都可以利⽤⽑管压⼒资料来确定[6]。⽬前主要有3种⽅法测定岩⼼⽑管压⼒曲线：半渗透隔板法、压汞法、离⼼法。压汞法测量速度快，对样品的形状、⼤⼩并⽆特殊要求，因⽽在⽯油勘探领域获得⼴泛应⽤。
汞对于绝⼤部分岩⼼都是⾮润湿相，如果对汞施加的压⼒⼤于或等于孔隙喉道的⽑管压⼒，汞就能够克服⽑管阻⼒进⼊孔隙。根据进汞的孔隙体积分数和对应压⼒，就能得到⽑管压⼒与岩⼼含汞饱和度的关系，⽣成压汞
法⽑管压⼒曲线[7]。假设孔隙空间是由不同半径的⽑管束组成，则⽑管压⼒与⽑管半径之间的关系如式(1)所⽰：
c c
2cos p r σθ
=
(1) 式中：c p 为⽑管压⼒；σ为表⾯张⼒；θ为润湿接触⾓；c r 为⽑管半径。
利⽤式(1)，可以获得岩⼼⽑管束分布。实验室条件下，σ=0.48 N/m ，θ=140°，则：
c c
0.735
p r =
(2) 1.2 ⼆维图像处理
利⽤岩⼼⼆维图像可以有效刻画岩⼼孔隙结构，为了⽅便对⽐研究，本⽂选择铸体薄⽚图像和CT 序列图像作为研究对象。
1.2.1铸体薄⽚图像
将玻璃单晶通过加压注⼊到岩⽯孔隙空间中，待玻璃单晶固化后磨制成为岩⽯薄⽚，由于孔隙空间中充满了有⾊胶体，在显微镜下可以⽅便地观察孔隙空间形态。图1(a)为某岩⼼样本的铸体薄⽚图像，暗灰⾊部分为玻璃单晶，表征孔隙相，其余部分表征岩⽯相。根据颜⾊或灰度
值对图像进⾏⼆值化处理，就可以得到图1(b)，图中⽤⽩⾊表⽰孔隙相，⿊⾊表⽰岩⽯相。 1.2.2CT 序列图像
CT 技术利⽤X 射线围绕岩⼼进⾏断⾯成像，可以在⽆损条件下直观地展⽰岩⼼内部的结构[8]。在图像中，由于孔隙相与岩⽯相对X 射线的穿透性不同，通过像素灰度值区分孔隙相与岩⽯相，获得岩⼼孔隙空间结构。
从CT 序列图像中抽出其中⼀张，如图2(a)所⽰，由于部分容积效应，圆柱形岩⼼的表⾯部分成像并不清晰，选择⽩⾊⽅框内的部分为处理区域，依据灰度值进⾏⼆值化处理，得到图2(b)。
1.3 岩⼼三维重建
在得到了岩⼼⼆维孔隙结构分布后，利⽤模拟退⽕算法重建岩⼼三维结构。算法从训练图像中提取统计特征函数，如孔隙度、两点相关函数、线性路径函数[1]等作为约束，通过迭代逼近最优结果。模拟退⽕算法本质上是⼀种优化算法，可以整合多张图像、多种统计函数，增强算法的灵活性与通⽤性，能够保证在训练图像并不理想的情况下依然得到合格的结果。模拟退⽕算法主要包括以下3个环节[3]：
(a)
before procesd (b) after procesd
Fig.1 Two-dimensional slice images of core
图1 岩⼼⼆维铸体薄⽚
Fig.3 Reconstruction results of simulated
annealing algorithm
图3 模拟退⽕算法重建结果
(a) before procesd (b) after procesd
Fig.2 CT images of core
图2 岩⼼CT 图像
790
太赫兹科学与电⼦信息学报第
13卷
1) 建⽴数学模型，确定解空间、⽬标函数和初始解。解空间就是重建结果可能的解的集合，⽬标函数就是统计特征函数，初始解是算法开始迭代的初始状态。
2) 确定状态产⽣函数，关键是确定候选解的产⽣⽅式和概率分布。在三维重建中根据概率，通过交换2个点的位置产⽣候选解，候选解产⽣的概率分布可以是均匀分布、正态分布、泊松分布等。
3) 确定初始温度和温度更新函数。
很棒
“温度”变量⽤于控制迭代的速度，⼀般来说，初始温度越⾼，获得最优解的概率越⼤，但所需时间也越长。因此，初始温度的确定应平衡考虑优化效果和优化效率的关系。温度更新函数给出了温度的下降⽅式，下降速度必须与状态产⽣函数相匹配。初始温度和温度更新函数构成冷却进度表，控制算法进程，使模拟退⽕算法逼近渐进收敛性，在有限的迭代次数内返回⼀个近似的最优解。
利⽤图1(b)作为训练图像，通过模拟退⽕算法构建的岩⼼三维结构如图3所⽰。其中，孔隙相使⽤灰⾊显⽰，岩⽯相则不显⽰。利⽤两点相关函数曲线和线性路径函数曲线衡量建模质量，图4列出了重建结果在X ,Y ,Z ⽅向
以及训练图像在X ,Y ⽅向上的两点相关函数曲线和线性路径函数曲线。
重建结果和训练图像的统计特性曲线基本相符，从统计学⾓度验证了重建的有效性。
1.4 最⼤球模型
获得岩⼼⼆维与三维图像后，为了计算岩⼼孔隙参数，需要提取孔隙空间拓扑结构并构建数学模型。最⼤球模型是提取孔隙拓扑结构的有效⽅法，以孔隙中的某点为中⼼⽣成⼀个球，不断增加球的半径，直到球⾯接触到岩⽯点，就得到了该点的最⼤内切球。遍历所有孔隙点，获得每个孔隙点的最⼤内切球，随后去除掉冗余球，孔隙空间中就充满了彼此重叠的球，组合这些球就构成了链路，这些互联互通的链路就反映了孔隙空间的拓扑结构。在链路的基础上对孔隙与喉道进⾏划分，以局部最⼤球为孔隙的起点，局部最⼩球为喉道的起点，以两者之间的某⼀球作为分隔球，将链路划分为孔隙与喉道。孔喉划分结束后，通过计算链路长度和内切球半径，就可以得到孔隙⼤⼩分布、平均孔喉⽐、平均孔隙半径等参数。该部分具体算法可见参考⽂献[4,9]。
2 实验数据
利⽤本⽂所述的⽅法，对多组岩⼼样本进⾏实验，均获得了类似的实验结论，由于⽂章篇幅所限，在这⾥列出2组实验结果。
2.1 实验1
2.1.1图像分析
图1(a)是⼀张铸体薄⽚图像，图像⼤⼩为1 290×972像素，点长度1.09 µm/pixel 。对该张图像进⾏⼆值化处理，可得到孔隙空间的分布状态如图1(b)所⽰。以图1(b)为训练图像，利⽤模拟退⽕算法进⾏重建，可得到三维图像，即图3，该图像⼤⼩为128×128×128像素，点长度为8.28 µm/pixel 。利⽤最⼤球模型分别对这2张图像进⾏孔喉尺⼨计算，得到结果如表1所⽰。
对岩⼼进⾏CT 扫描，可获得序列图像，如图2所⽰。图像⾼980个像素，宽1 013个像素，点长度3.41 µm/pixel 。完成整套CT 扫描，共获得990张图像。采⽤章节1.2.2中所述⽅法进⾏批处理，并去除序列图像⾸尾效果较差的部分，得到800张⼆维⼆值图像，组成了⼀个长宽⾼均为800个像素、共512×106个体素点的三维⽴⽅体。⾸先从800张图像中沿竖直⽅向每隔10张抽取1张图像，共80张⼆维图像，在⼆维平⾯内进⾏孔喉划分与参数计
Fig.4 Statistical function curves of reconstruction and training image
图4 重建结果与训练图像统计函数曲线
(a) autocorrelation function (b) linear path function
0 10 20 30 40 50 60
r /pixel
1.0
0.80.60.40.20-0.2
SA reconstruction result X Y Z training image X Y
l i n e a r p a t h f u n c t i o n
第5期吴运强等：数字岩⼼重建与孔隙参数分析 791
算，并对结果求平均，获得⼆维分析结果。随后，直接将三维图像输⼊到最⼤球模型中进⾏分析，获得三维分析结果。两项结果均列在表1当中。
岩⼼孔隙⼤⼩的分布状态对于油层储集与渗流特性有重要影响，实验通过绘制直⽅图的⽅式，可以直观地呈现孔隙⼤⼩分布状态，由表1中4组分析所得的结果如图5所⽰。
2.1.2压汞法测量⽑管压⼒曲线
女子防身术踢裆
对同组岩⼼，利⽤压汞法测定⽑管压⼒曲线，换算得到岩⼼⽑管半径⼤⼩及分布状态，实验结果如图6，左侧图像为⽑管压⼒曲线图，其中位于上⽅的曲线为进汞曲线，位于下⽅的曲线为退汞曲线，右侧直⽅图为⽑管分布直⽅图。实验数据见表2。
表1 实验1计算机图像分析结果
Table1 Results of computer image analysis in experiment 1
porosity/% average pore radius/µm average throat diameter/µm pore throat ratio
slice image 13.57 12.95 12.34 3.23 reconstruction result 13.57 15.94 16.06 3.39 2D CT image 11.90 27.59 28.05 3.37 3D CT image 11.42 20.91 20.03
4.31
表2 实验1压汞法测量结果
Table2 Results of mercury porosimetry in experiment 1
商品谷物农业
parameters value parameters value
porosity 20.2% average capillary
radius
7.16 µm
displacement
pressure 0.02 MPa max capillary radius 31.02 µm
median saturation pressure 0.43 MPa median saturation
radius
1.70 µm
Fig.6 Curve of mercury porosimetry in experiment 1
图6 实验1压汞法测量结果
100 80 60 40 20 0 mercury saturation/% 40 20 0 mercury saturation increment/%
1001010.10.01
0.001
Fig.5 Histogram of the pore radius in experiment 1
图5 实验1孔隙半径分布图
(c) two-dimensional CT quence images (d) three-dimensional CT quence images
2-1010-1717-2525-3232-4040-4747-5555-6363-7070-7878-8585-9393-100100-108108-115115-123123-130130-138138-145145-1535040
302010
羞涩反义词pore radius/µm
p o r e s n u m b e r
(a) two-dimensional cast thin ction images (b) reconstruction results of simulated annealing
0.700.600.500.400.300.200.100
牛腩肉怎么做好吃0-10 10-20 20-30 30-40 40-50 ＞50
pore radius/µm
pore radius/µm
70605040
302010
大德如阳
792
太赫兹科学与电⼦信息学报第
13卷
2.2 实验2
选择另⼀组岩⼼进⾏实验，采⽤与实验1相同的⽅法进⾏分析，将数字岩⼼分析结果与压汞实验结果分别呈现在表3、图7、表4、图8当中。
