数字岩心多分辨率存储及连通区域分析

彭成

(中国石油化工股份有限公司石油勘探开发研究院， 北京 100083)

随着计算机技术的发展，运用图像技术将电镜或电子计算机断层扫描(computed tomography，CT)扫描后的二维和三维岩心数据进行处理，进而形成数字岩心[1]。数字岩心可用多种计算机算法进行建模、孔隙度分析、孔喉结构分析等，是目前较为常用的岩心分析手段。王云龙等[2]基于综合微计算机断层扫描技术、数字岩心分析技术，实现低渗储层数字岩心渗吸过程的模拟。宋文辉等[3]对数字岩心进行纳米孔隙流动数学模型建立，研究了多尺度孔隙结构特征。

数字岩心的建模方法早期是通过二维图片，通过一些典型算法建模出三维图像，例如蒙特卡洛、模拟退火、过程法等。罗家荣[4]采用过程法模拟数字岩心，通过模拟岩石沉积过程建立数字岩心。

随着计算机处理能力和扫描能力的进步，三维数字岩心图像现在可以直接通过扫描获取。在此基础上，对数字岩心进行相应算法分析，构建孔隙模型、孔喉结构等。王志兵等[5]采用Avizo软件中先进的数字图像处理技术提取孔隙网络模型。柳青兵等[6]采用传统拼接方法，通过相同孔径条件下流过孔隙结构的流体体积变化相等的拼接原理确定不同岩性样品的拼接点位置，进而确定完整孔隙大小。

对于数字岩心存储，目前多采用单一文件存储，数据量大，查看局部数据时速度较慢，且没有区分粗细粒度，应该进行分割存储并存储不同粗细级别的数据来方便查看[7-8]。陈国军等[9]基于分层四叉树对二维数字岩心进行了多分辨率数字岩心体素模型建立，不过没有实现三维数字岩心相关的处理。

对于孔隙的连通性，各个连通区域之间的区分程度不够高，应该将不同连通区域进行分色来帮助显示[10-11]。对于每个连通区域的面积或体积，范围，中心位置以及其他一些信息，应该建立连通区域和其属性的对应关系。

针对上述问题，现通过四叉树和八叉树分块方法对数字岩心进行分布式多级分辨率存取，根据浏览范围选择粒度层级；通过扫描线、三角剖分、四面体剖分和地图染色方法识别和划分连通区域；基于移动立方体算法求取岩心孔隙边界面。通过上述方法，以期为数字岩心的高效存取、连通区域划分和性质分析、孔喉结构的建模和形态分析提供支持。

1 数字岩心多分辨率分块存储

为提升数字岩心的存取效率，采用四叉树和八叉树分别对二维和三维数字岩心进行分块切割存储，使得在获取局部范围内的数字岩心数据时只需读取区域对应的分块即可。将分块数据逐级向上合并形成多级分辨率，并在显示较大范围数字岩心时，根据屏幕像素选择适合的分辨率级别读取对应的分块。

1.1 分块参数配置

数字岩心分为二维三维两种类型，格式主要包括RAW和BMP等。对于RAW格式，除岩心数据本身外不带有其他信息，需要给出岩心的X、Y、Z(长、宽、高)像素数量。对于BMP格式，文件中带有X、Y(长、宽)两个方向的像素值，可以直接获取，Z方向像素数量需要另外给出。接下来指定分块大小，二维岩心设置分块的X、Y像素数量，三维岩心设置X、Y、Z像素数量，后续会将数字岩心按照设定的分块大小切割成块。

分别选取二维和三维RAW格式文件，其中二维RAW文件的X、Y像素值分别为110 342和104 900，文件大小为10 GB，四叉树子块设为长宽1 024，即子块大小为1 MB；三维RAW文件的X、Y、Z像素值分别为1 200、1 200、1 600，文件大小为2 GB，八叉树子块设为长宽高128，即子块大小为2 MB。子块的大小由用户指定，为提升子块分割效率，依照四叉树和八叉树的性质，子块的长宽高大小选取的依据为：尽量保证分割后各个方向上切割出的数量相近，并接近且小于等于一个2的次幂的个数，这样可以保证各个方向上切分后的层级相同，且空块较少。

1.2 分块切割

基于四叉树和八叉树对数字岩心进行分块，首先计算树的层级数，对于二维数字岩心数据，四叉树的层数需同时满足关系式如下。

2(n-1)Xblock≥Xfile

(1)

2(n-1)Yblock≥Yfile

(2)

对于三维数字岩心数据，八叉树的层数还需满足关系式为

2(n-1)Zblock≥Zfile

(3)

式中：n为层数；Xblock、Yblock、Zblock分别为子块X、Y、Z像素值；Xfile、Yfile、Zfile分别为文件的X、Y、Z像素值。

满足这些不等式的最小层数作为树的层级数，依此计算得到本文选取的二维数字岩心数据四叉树的层数为8，对应子块的总个数为16 384个，三维数字岩心数据八叉树的层数为5，按照八叉树规则得到对应子块的总个数为4 096个。基于四叉树和八叉树编码规则对子块赋予编号[12-13]，如图1所示，为八叉树空间结构划分及线性莫顿编码示意图，按照图中左侧所示对三维数字岩心进行切割，每个子块的编码如图1中右侧所示。二维数字岩心按照类似方式基于四叉树进行切割和编码。

图1 八叉树空间结构划分及线性莫顿编码Fig.1 Octree spatial structure partition and linear Morton coding

接下来按子块编号的顺序依次读取对应位置的数字岩心数据，对于靠近边界区域的子块，其长宽高可能达不到设定的子块大小，截取到边界即可，对于边界外的子块，生成空子块。

数字岩心数据量较大，各分块相对独立且处理方式一致，适用于用多线程的方式进行读取，按照四叉树和八叉树的性质，为二维数字岩心设置4的幂级数量的线程，为三维数字岩心设置8的幂级数量的线程，在例子中，设定的二维和三维线程数分别为16和8。为每个线程分配数量相同的子块进行处理。

1.3 分块合并

在取子块数据的过程中，同时进行子块合并工作。对于二维数字岩心数据，每获取完4个子块进行合并，如图2(a)所示，即在X、Y方向间隔获取像素进行粗化，红色格子为所选像素，4个子块间隔抽取合并出右边粗化子块，即在长宽方向上各进行50%的压缩，合并完的子块即上一粒度的粗子块，实际数据效果如图2(b)所示。合并后将这4个块所占内存释放。当上一级子块数量也达到4个时，将它们进行合并存储。依此直到最顶层。

图2 二维数字岩心四叉树分块合并示意图Fig.2 Two dimensional digital core quadtree block merging

对于三维数字岩心数据的合并，与二维类似，每获取8个子块进行合并成上一级粗粒度子块，每级粒度的子块满8个时都进行合并，直到所有子块合并完毕。按照边读取边合并的方式可以大幅度减少内存，正常方式读取所有子块需要的内存为数字岩心文件大小，而边读取边合并的方式占用的内存最多为

M=3Mblocklevel

(4)

式(4)中：level为树层级数；Mblock为单个子块占用内存大小。

对于多线程的模式，各个线程单独完成上面的操作，每个线程做完所有的合并操作后，再将各个线程的块按照同样方法合并起来，直到最顶层。例子中二维线程数为16，即四叉树的第二层，再将它们进行两级合并即可。三维线程数为8，为八叉树的第1层，再进行一级合并即可，最后将这些线程子块合并到最上层即树的第0层。

1.4 读取粒度与区域内分块选择

在数字岩心文件导入并切割存储完成后，可以进行查看，如图3(a)所示，为本文例子中的二维数字岩心图像，图3(b)为三维数字岩心的3个方向切片切出来的立体图像，可以指定任意方向的切片来查看三维数字岩心。在显示方面，由于生成了不同粗细粒度的图像，通过不同的四叉树和八叉树级别表示，越上层的图像越粗。在选择以何种粗细粒度显示时，选择的粒度如果过细，则受屏幕分辨率限制不能完全显示所有像素，如果选择的粒度过粗，则展现的像素量太少。选择粒度的准则是在屏幕分辨率的限制下，尽可能选择更细的粒度。

图3 二维三维数字岩心显示效果Fig.3 2D and 3D digital core display effect

对于数字岩心图像粗细粒度的选择，即四叉树或八叉树级别的选择，二维三维数字岩心所取级别需分别满足以下关系式。

Mapfile/(Pixblock4lev-1)

(5)

Mapfile/(Pixblock8lev-1)

(6)

式中：Mapfile为文件地理坐标范围，式(5)中为面积范围，式(6)中为体积范围，在导入时设定，默认取值为其像素坐标范围；Pixblock在式(5)中为子块像素面积，在式(6)中为体积；lev为所取级别；Mapscr和Pixscr分别为屏幕地理坐标范围和像素范围。式(5)和式(6)分别对应二维和三维数字岩心。

对于局部区域岩心数据查看，根据屏幕所显示的地理坐标范围，以及数字岩心的地理坐标范围的位置区域，得到对应要查看的岩心区域，读取区域内子块数据，子块层级选择按前面所述。

1.5 存取效率及粒度选择实验对比

实验运行的环境为Windows 10，64位操作系统，4核8 GB内存，Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz 2.30 GHz。对1.1节中所述的两个岩心数据进行分布式存储。首先搭建分布式服务器，如图4所示，建立GFS Master服务器负责子块的切分和获取；建立GFS Chunk存储服务器来存放子块数据，最后共建立了3个存储服务器100个虚拟存储节点来实现分布式存储。

图4 分布式服务器架构Fig.4 Distributed server architecture

按照1.1节的切块配置进行切块，分布式存储到各个节点中，各数据的切块结果如表1所示。

表1 实验数据切块结果Table 1 Experimental data segmentation results

与本地文件整体存储相比，分布式分块存储一方面减少本地存储空间开销，另一方面在访问局部区域数据时只需要读取对应子块，减少文件读写开销。表2对比了整体与分块存储在访问一个子块区域对应的数据时，文件读写的跨度。

由表2可以看出，对于整体存储模式数据连续存储，在访问子块区域数据时，为读取子块区域对应的数据，那么会不断跨整个宽度和深度数据去访问下一个位置的数据，造成了文件访问跨度的增加。而分块存储只需要访问对应子块文件，不需要额外的跨度。

天津市高等学校师资培训中心（以下简称天津市中心）成立于1990年，受天津市教委和天津师范大学双重领导。中心下设行政办公室、岗前培训办公室、信息技术部、教师资格认定办公室，当前业务以岗前培训、教师资格认定、高校教师网络培训为主。

表2 访问区域数据时文件读写跨度对比Table 2 Comparison of file read and write span when accessing area data

在访问屏幕范围内数据时，会选择合适的粗细粒度子块数据进行读取，表3显示了与选择原始数据读取相比，选择适合粒度层级进行读取在数据量上的对比，选择的屏幕分辨率为1 024×1 024。

表3 选择合适读取粒度与原始数据对比Table 3 Select the appropriate reading granularity and compare it with the original data

可以看出，多级分辨率的岩性数据，根据浏览范围选择合适的粒度层级展现，依照1.4节中的方法即充分利用屏幕分辨率，又整体上减少了数据读取量，屏幕显示范围越大，减少的数据量越多。

2 数字岩心连通区域识别

数字岩心的孔隙分析、孔喉结构模型计算、孔喉连通区域计算等是数字岩心分析的重要方面。孔隙的识别通过设定的像素阈值进行；对孔喉结构连通性的计算，首先基于扫描线方法，先计算连通线，再将相邻连通线合并成连通区域，最后将各分块连通区域合并成整体连通区域。

2.1 连通线扫描

首先进行扫描参数的设定，包括阈值、扫描线最小长度、连通区域面积或体积的下限等。依次读取最细层级的每个子块进行扫描，在扫描过程中，二维按照X、Y的顺序，三维按照X、Y、Z的顺序依次读取子块的每个点的像素值，对于满足阈值的像素点，相邻的连接起来形成线段，即扫描线，如果扫描线长度小于最小长度则剔除掉。

在例子中，二维数字岩心阈值设定为55，扫描线最小长度为3，连通区域面积下限为22，三维数字岩心阈值设定为55，扫描线最小长度为3，连通区域体积下限为50。即扫描二维数字岩心过程中，像素值小于55的点会记录下来，连续3个或以上这种点可以组成一条扫描线，若干个相邻的扫描线组成的连通区域总面积还要大于等于22；扫描三维数字岩心过程中，像素值小于55的点会记录下来，连续3个或以上这种点可以组成一条扫描线，若干个相邻的扫描线组成的连通区域总体积还要大于等于50。

2.2 分块连通区域生成

生成完子块所有的扫描线后，将扫描线相邻合并形成连通区域，从而得到这个子块的若干个连通区域。对于二维数字岩心，如果扫描线的Y方向距离小于等于1且X方向不相离，则相邻；对于三维数字岩心，如果扫描线的Y、Z方向距离均小于等于1且X方向不相离，则相邻。依照此判断规则将所有扫描线组成一个个连通区域，生成一个连通区域的方法为：通过堆栈的方式，每次在未分组的扫描线中选取一条，判断是否与当前堆栈顶部的扫描线相邻，如果相邻则加入连通区域和堆栈。所有扫描线选取完毕后将堆栈弹出一个对象添加到连通区域。接下来进行第二次循环，操作与第一次相同，以此类推，直到堆栈为空，便生成了一组扫描线表示的连通区域。如图5所示，具有相邻关系的一组连通线组成一个连通区域，共形成了3个连通区域。

图5 扫描线合并连通区域示意Fig.5 Schematic diagram of scan line merging and connecting area

2.3 合并相邻分块连通区域

将相邻分块的连通区域进行合并，对于二维数字岩心，相邻分块为上下左右4个，对于三维数字岩心，相邻分块为上下左右前后6个。对于连通区域存在相邻的两个分块，将它们的连通区域扫描线进行合并得到新的连通区域。

连通区域是否相邻的判断方法为：先判断两个连通区域的范围是否有交集，如果没有交集，则不相邻，如果有交集，则进一步判断：将它们所有的扫描线二维按照Y，三维先按照Z再按照Y进行排序，然后扫描线两两比较是否相邻，如果Y或者Z距离大于1的直接跳过。一种快速的方法是，每次两边只取二维Y方向相邻三行，三维Z方向相邻三行的数据进行比较。对于三维，进一步对这些数据每次取Y方向相邻三行进行比较。如果有相邻的扫描线，则表示两个连通区域相邻。所有分块的相邻连通区域完成合并之后，最后将合并后的各个分块的连通区域信息作为整体的连通区域信息。如图6所示，相邻两个分块的两个连通区域在边界处相邻，子块合并后这两个连通区域合并成一个连通区域。

图6 相邻分块连通区域合并示意Fig.6 Merging of adjacent block connected areas

2.4 合并整体连通区域

3 连通区域分色及边界面计算

为方便查看各个孔喉连通区域形态，需要对其进行染色，基于三角剖分和四面体剖分对二维和三维连通区域分色，并结合地理信息技术将各连通区域作为图元存储形成图层，图元属性包括颜色、面积或体积等，方便查询和管理。基于移动立方体(marching cube，MC)算法[14]对三维数字岩心模型孔隙求取边界面，以便在三维空间中展示孔喉结构。

3.1 三角和四面体剖分连通区域及分色

对各个连通区域的中点进行三角/四面体剖分连接起来。二维采用三角剖分算法，三维采用四面体剖分算法[15]。如图7所示为例子中二维数字岩心连通区域进行三角剖分后的结果。

图7 二维数字岩心三角剖分结果Fig.7 Two dimensional digital core triangulation results

利用三角/四面体剖分的结果，可以得到各个连通区域中心点之间是否有边相连，以及每个点所连的边的数量，如果两个点有边相连，则表示这两个连通区域距离较近，需要染不同颜色。具体每个连通区域染哪种色的计算方法如图8所示。

图8 连通区域染色流程Fig.8 Dyeing process of connected areas

按照地图染色算法，统计各个点连接的边数量，从大到小排序，从边数量最多的点开始染色，每染色一个点，将其四周有边相连的点的边计数减去1，然后将剩余点按照所连接的边的数量从大到小排序，继续从最多边数量点染色，每次染色要与周围有边相连的点的颜色不同。然后依此类推循环下去，直到所有点都被染色。

接下来对于二维数字岩心，将染色结果用地理信息图层表示。每个连通区域作为图层中一个图元对象，并赋予所染颜色。图元对应属性存放这个连通区域的中点坐标、面积、范围、扫描线集合等。在例子中，二维数字岩心生成的地理信息图层如图9所示。

图9 二维数字岩心连通区域识别Fig.9 Recognition of connected area of two-dimensional digital core

对于三维数字岩心，构建三维模型并将各个连通区域作为三维图元对象，存放的属性信息除扫描线集合、中点、范围外，还存储了连通区域的像素体积和地理上的体积。在例子中，三维数字岩心的连通区域染色效果如图10所示，展示了3个方向上的切片。

图10 三维数字岩心连通区域识别Fig.10 Recognition of connected area of 3D digital core

3.2 连通区域边界面计算

三维数字岩心的连通区域表示孔洞，将孔洞剔除掉，剩余的部分代表喉道，构建喉道三维模型的流程如图11所示。

图11 构建喉道三维模型流程Fig.11 Process of constructing laryngeal three-dimensional model

对数字岩心的三维空间范围中每个像素点，如果其属于某个连通区域内，则赋值为0，如果不属于连通区域内，则赋值为1，然后用MC算法，对每个立方体进行处理，得到各个立方体的01分界面，然后合并成整体的01分界面即孔洞分界面[11]。具体的方法为，按Z方向每次取相邻的两层切片，切片中夹了若干个立方体，数量为

n=(Xpix-1)(Ypix-1)

(7)

式(7)中：Xpix为长像素值；Ypix为宽像素值。

对于立方体的8个顶点，01分界面共有256种情况，可以分为15类，如图12所示，列出了不同情况的等值面构建方法。将每个立方体形成的等值面连接起来，即是孔喉的分界面。连接时，与连通区域分组类似，将具有相连关系的分界面划分为一组，属于同一个连通孔喉，互不相连的分界面属于不同孔喉。

图12 移动立方体等势面构建情况组合Fig.12 Combination of construction of equipotential surface of moving cube

例子中将三维数字岩心的等值面生成后的孔喉结构效果如图13所示，可以用两种显示方式，图13(a)～图13(d)是显示等值面，图13(e)～图13(h)是只显示等值面的边即镂空模式，从左到右分别是从小到大设定不同阈值得到的结果，阈值分别为50、70、90、110。阈值设定较小，得到的连通区域小，剔除后剩余的部分较多；阈值设定大，得到的连通区域大，剔除后剩余的部分较少。

图13 三维数字岩心连通区域边界面Fig.13 Edge interface of 3D digital core connected area

4 结论

(1)分别以四叉树和八叉树模式分块存储二维和三维数字岩心数据，生成多级分辨率的岩性数据，根据浏览范围选择合适的粒度层级展现，整体上减少了数据读取量。

(2)通过扫描线合并构建连通区域，基于移动立方体方法获取孔隙分界面，得到了孔喉结构的形态及连通性质。通过三角剖分和四面体剖分对连通区域分色，方便区分显示各连通区域。