CN114494637A - 一种基于结构体矩阵的砂岩三维真实模型重构方法 - Google Patents

Abstract

本发明公开了一种基于结构体矩阵的砂岩三维真实模型重构方法，它包括运用X‑射线成像设备获取砂岩标准试件几何拓扑结构图像数据；将图像数据转换成离散三维数字化矩阵；运用Chan‑Vese模型对砂岩数字矩阵进行两相分割；建立结构体矩阵模式库：基于砂岩数据的结构体矩阵确立对应的中心体素单元模式和坐标信息；通过结构体矩阵反演得到几何拓扑结构数据和线性插值，得到砂岩两相微观结构的点坐标和方向矩阵，从而实现砂岩三维模型重构过程和结果可视化。本发明的优点是：运用结构体矩阵完成几何拓扑结构数据的存储、搜寻和恢复，提高了计算效率，降低了数据存储内存消耗，减少了三维重构模型与真实试件之间的误差。

Description

一种基于结构体矩阵的砂岩三维真实模型重构方法

技术领域

本发明属于岩石微观结构图像处理技术领域，具体涉及一种基于数字化结构体矩阵建立砂岩三维真实模型的数字化重构方法。

背景技术

砂岩是常见的建筑材料、地下水资源、矿产资源和石油等能源的储层岩石。这些深部资源的开发与储层岩石的微观几何结构和破坏特征紧密相关，准确建立岩石三维微观结构可以有效还原岩石微观几何结构演化特征，揭示岩石变形破坏和内部流体流动规律，能为深部岩石破裂致灾机理、灾害防控技术及深部资源的开发技术提供一定的理论支持。

现有砂岩三维模型重构主要分为面重构(如移动立方体技术)、体重构技术(如射线追踪技术)和基于树结构模型的重构算法。树结构重构算法主要是通过树形模式结构(一层一层依次二值化向下发展的结构)建立相应数据，然后重构三维模型。这些算法中，移动立方体算法最为简洁，它是一种经典的面重构技术，它的基本思想是在给定的三维数据的每一个体元(一般由八个体素组成)中寻找顶点数值与给定数值的关系，并根据这种分布关系寻找等值面，最后用三角形将等值点连接组成等值面。移动立方体技术的三维模型重构过程包括以下步骤：步骤1、获取研究对象三维数据并读取到内存中；步骤2、提取计算单元8个顶点组成立方体，将立方体顶点值与给定等值面值比较，建立立方体索引；步骤3、在查找表中搜寻索引点，获取三角面片顶点准确位置和法线方向；步骤4、输出三角形面片顶点和法线数据，进行可视化显示。虽然这些技术可以快速建立研究对象的三维模型，但是它们在实际工程中建立岩石三维模型时都具有以下的问题：

1、由于岩石内部空间特征复杂，且岩石内部各类微观结构相之间存在相互影响，部分曲面和空洞三维特征信息识别精度低、容易丢失部分三维模型空间特征信息。

2、三维模型重构过程存储数据多，存储方式未优化，且现有技术数据存储结构冗杂，导致数据搜寻传递过程需要耗费大量时间，因此三维模型重构过程计算效率低、数据存储内存消耗大；

3、由于不同重构算法对岩石内部微观结构识别和分割、空间特征信提取精度较低，导致部分微小孔隙和曲面特征信息等丢失，所以三维重构模型精度与实际研究对象模型存在一定程度差异、特别是微观结构特征。

结构体矩阵是在Python编程语言中定义的结构体矩阵(或结构体数组)，用于存储各种分类特征数据等信息，例如A为一个结构体矩阵，它的形式为 A＝{Ax；Ay；Az；……}，其中每个分类(如Ax)可以是由数字、字符串等标准的向量、数组或者矩阵。结构体矩阵的优点：1、有顺序的数据存储，剔除了现有技术中数据存在中多余空间消耗；2、数据传递、搜寻和恢复过程可精准定位提取特征信息，依次地快速提取有效特征信息(如孔隙特征)。

发明内容

针对现有技术存在的问题，本发明所要解决的技术问题就是提供一种基于结构体矩阵的砂岩三维真实模型重构方法，它能够提高三维模型重构过程的计算效率、降低模型数据存储内存消耗，保留更多的砂岩空间微观结构特性信息，提高模型精度，从而降低重构模型与真实岩石差异。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：

步骤1、选定需建模型的砂岩块体，制作具有一定尺寸的圆柱体或长方体标准试件，运用X-射线成像设备获取砂岩标准试件几何拓扑结构图像数据；

步骤2、将获取的砂岩几何拓扑结构图像数据转换成离散三维(3D)数字化矩阵：

步骤3、运用Chan-Vese模型将步骤2获得的砂岩数字矩阵进行精准的两相分割，计算获得的砂岩两相数字矩阵；所述两相是指岩石内部孔隙相和固体基质相；

步骤4、建立结构体矩阵模式库：采用三维(3D)扫描模板对砂岩两相数字化矩阵进行顺时针数据扫描，提取砂岩微观结构的二进制模式、中心体素单元和像素坐标，并将提取的数据由二进制转化为十进制存储在结构体矩阵中；

所述3D扫描模板相当于一个小型矩阵，如3×3×3形式的模板为 [1,1,1；1,0,1；1,1,1]；所述二进制模式为用模板扫描之后得到的由二进制(0,1) 组成的模式单元，所述中心体素单元为扫描模板中心U对应的单元。

步骤5、计算中心体素单元的分布概率，利用概率密度分布函数确立未知中心体素单元相状态，补充建立的结构体矩阵模式库，确立对应的中心体素单元模式和坐标信息；

步骤6、运用结构体矩阵反演得到的几何拓扑结构数据进行线性插值，获得表征砂岩两相微观结构的点坐标和方向矩阵；

将得到的点坐标和方向矩阵进行可视化显示，获得砂岩两相微观结构的三维数字化重构模型。

本发明的技术效果是：

本方法发明采用结构体矩阵存储、搜寻和恢复岩石微观几何拓扑结构数据；一方面能直接提取岩石微观几何特征信息(包括坐标、几何尺寸等)，采用根据扫描模板中心点的单元相来记录并存储二进制模式库的方式降低数据存储内存消耗，采用十进制代替二级制方式完成存储、搜寻和数据恢复过程，提高了计算效率；另一方面，结构体矩阵剔除了图像分类聚束过程导致的微小特征信息的数据处理误差，能够精确提取微小结构特征信息，因此重构的三维模型误差小。

附图说明

本发明的附图说明如下：

图1为训练数据和模板扫描示意图；

图2为结构体矩阵模式库构建示意图；

图3为数据结构体矩阵示意图；

图4为初始二值训练图像、基于树结构模型重构的二维图像和基于矩阵结构体重构的二维图及三维砂岩模型。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：

本实施例选取重庆砂岩为模型重构对象，其三维模型重构过程包括以下步骤：

步骤1、选择指定的砂岩块体，制作一定尺寸的圆柱体标准试件，运用X-射线成像设备扫描砂岩标准试件，得到包含微观几何拓扑结构的砂岩分层图像数据。

X-射线分层图形是通过X-射线接收平面形成的三维图像矩阵，这个三维图像矩阵通常由一层一层的二维分层图像组成。因此，X-射限分层图形是一个三维立体图。

步骤2、将砂岩微观几何拓扑结构的分层图像数据离散成三维离散标量场数据，构建砂岩三维数字化矩阵:

式(1)中：I_DM为包含砂岩几何拓扑结构特征的3D离散数字化矩阵； f为数字化矩阵体素单元；所述体素单元等同于像素单元，均代表成像坐标系中的最小单元。gl为矩阵单元强度值，它由像素灰度值GI决定(0-255)； x，y，z为成像坐标系中在x轴、y轴和z轴方向的坐标；N为矩阵维度。

步骤3、根据文献“Active contours without edges.”，Chan T F&Vese L A(2001)， IEEE Transactions on Image Processing,10(2),pp266-277(“无边缘活动轮廓”，Chan T F&Vese L A(2001)，IEEE Transactions on Image Processing，第10卷，第2期，第 266-277页)记载的Chan-Vese模型精准识别砂岩孔隙相和固体基质相，识别过程的数学表达式为：

为了简化计算，这里使用水平集函数Φ替换轮廓线函数：Ω_outside代表轮廓线以内的区域。

式(2)、(3)中，E(Ω,f₁,f₂)是能量函数，

是梯度函数算符；I为单位矩阵， L为识别对象的轮廓线长度；H为Heaviside函数；δ为Dirac函数；f₁和f₂分别为轮廓线内单元函数和外单元相函数(内单元相函数和外单元相函数指的是由轮廓线划分的内外区域中的单元，若内单元相为孔隙，则外单元为固体基质，反之亦然)，它们的表达式为：

将式(4)、(5)联立公式(3)进行积分求解运算，可以获得砂岩的两相微观结构数字矩阵，它的表达式为：

式(6)中，

为砂岩两相的3D数字化矩阵，f_b(即f_b(x,y,z)的简写) 为数字化体素单元；f_b＝0时，体素单元代表孔隙或/和空隙，f_b＝1时，体素单元代表固体基质。

步骤4、建立结构体矩阵模式库：

1)、采用3D扫描模板的投影模板对砂岩两相3D数字化矩阵

的投影图进行顺时针数据扫描，依次完成三维空间在各个投影平面(包括XOY、 XOZ和YOZ平面)的数据扫描。

图1演示了模板在XOY平面的扫描过程，图1(a)表示3D训练图像示意图，它由公式(6)计算获得，然后投影在XOY平面形成图1(b)。图 1(d)为建立的三维(3D)扫描模板，尺寸3×3×3(模板尺寸可以任意，它可以根据使用者自行设定，只要不超过训练图像尺寸)，然后投影到XOY平面形成图1(c)。图1(c)作用在图1(b)过程就是结构体矩阵模式库单元采集过程：

扫描模板(如3×3×3模板，在XOY平面模板尺寸为3×3)顺时针扫描砂岩3D数字化矩阵；提取微观结构二进制模式、中心体素单元和像素坐标。中心体素单元是指扫描模板中心点U对应的

体素(值为0或1)，二进制模式是指扫描模板根据模板中心U周围对应点按顺时针方向1-2-3-4(上 -右-下-左)依次提取的体素单元(如1-0-1-1)，即它的二进制模式为1011，像素坐标即为在X轴、Y轴和Z轴对应的值。

由设定的尺寸模板(如3×3×3)依次扫描步骤3得到的三维数字化矩阵，它的表达式为：

式中：f(U)为扫描模板读取矩阵的中心体素单元(即模板中心U对应的单元)的相状态(0或者1)，PL(U)为采集的模式库单元，所述模式库单元是指根据模板识别的二进制模式，包含二进制模式单元和中心体素单元相状态。二进制模式单元f(U+h_α)如图3中列1的1011,1001等；h_α为扫描方向，其中α＝α_x,α_y,α_z分别表示X轴方向、Y轴方向和Z轴方向。

图3中列1的1011,1001是根据公式(7)中的f(U+h_α)获得的二进制模式单元，这些数据在图3中位置“列1”之下列出。f(U+h_α)获得二进制模式的过程结合图1进行说明：当α＝α_z(即

)时，表示扫描模板图1(d)投影到XOY平面形成图1(c)的扫描模板，运用图1(c) 模板扫描图1(a)投影在XOY平面形成的图1(b)，就可以得到相应的f(U) 和另一个f(U+h_α)(此时为

)。如当第一次扫描时，图1(c)中中心点 U对应图1(b)的单元为0(图1(b)中，黑色为0，白色为1。)，即f(U)＝0，而模板扫描位置1,2,3,4依次读取的数据则为1,0,1,1，这些数据构成了二进制模式单元 1011，即

因此，式(7)形成的其中一个模式库单元为0|1011，其中0为中心体素单元相状态，它存储在图3的列2中，分隔符“|”之后的四位数“1011”为二进制模式单元，它储存在图3的列1中。这就是根据扫描模板中心点的单元相来记录并存储二进制模式库。

2)、沿着XOY、XOZ和YOZ平面，分别提取砂岩微观结构的模式库单元、中心体素单元和像素坐标，计算中心体素单元相状态的概率分布。

提取模式库单元的过程如图2所示：

图2为结构体矩阵模式库构建示意图，图中所用模版和训练图像均为 5×5尺寸(模板尺寸可以任意，它可以根据使用者自行设定，只要不超过训练图像尺寸)，将图2(a)的模板作用在图2(b)的训练图像上，按照图2 (a)中顺序顺时针依次读取训练图像图2(b)的单元，则可获得不同的子特征模式库，如图2(c)所示。

计算中心体素单元相状态的概率分布的表达式为：

式(8)中：CDF为条件概率密度函数；prob为条件概率密度算符；PL_j(U) 为不同扫描次数下的总模式库个数；

为第j个模式库单元的中心体素单元相状态(k＝0或1)的次数，Num为扫描次数，j为第j次扫描；k为模式库单元的中心体素单元相状态，取值0或1。

概率计算是为了根据相应的二进制模式库单元恢复未知中心体素单元相，它在步骤5中的子步骤2)和3)中有具体说明。

3)、将提取的数据由二进制转化为十进制存储在结构体矩阵中，如图3 所示。图3为根图2建立的结构体矩阵示意图，用于存储、搜寻和恢复岩石微观特征信息。图3中左边三列为经过图2过程得到的二进制子模式库，列 1为二进制模式，这些二进制模式是根据模板扫描原矩阵的不同位置，即不同的中心点位置得到的，在图3中“列3-坐标信息”一栏标明。列2为中心体素单元相(0或1)的次数，列3为坐标ID；图3中右边三列为经过图3 中坐标三列转换后的十进制子模式库，其中列4为对应于列1转换后的十进制数，列5和列6分别为中心体素单元相(0或1)的次数和坐标，由于次数和坐标信息简单，所以均保持不变。

结构体矩阵的表达式为：

A_struct＝{XOY＝(A_Sub-pattern),XOZ,YOZ} (9)

式(9)中，A_struct为结构体矩阵；XOY＝(A_Sub-pattern)为在XOY平面获取的子结构体矩阵模式库A_sub-pattern，其表达式为：

式(10)中：S为一个完整模式库(PL(U)仅为一个完整模式库中的一个模式库单元)；M为模式库个数，分别有S₁、S₂…S_M。B-pattern、Counters、Coordination ID、D-pattern分别为二进制模式单元、中心体素单元相(为0或1) 的计数次数、坐标信息和十进制模式。B-pattern、Counters、Coordination ID对应于图3中的列1、列2和列3；D-pattern,Counters,Coordination ID对应于图3 中的列4、列5和列6；X0Z、YOZ模式库与XOY结构类型相同，但是取值不同。

步骤5、利用概率密度分布函数确立未知中心体素单元的相状态，补充建立的结构体矩阵模式库，确立对应的中心体素单元模式和坐标信息。

1)、假设未知中心体素单元属于的模式库单元中第1和第3项均为1，根据式(9)和(10)可知，它的模式库单元表达式为：

PL(U)＝ζ|(1,ζ,1,ζ) (11)

式(11)中：ζ为未知中心体素单元的取值。

说明：式(11)与式(7)物理意义相同，式(7)表示扫描三维数字化矩阵提取不同的模式库单元，而式(11)表示一个包含中心体素单元相状态和二进制模式单元的具体模式库单元。

2)、利用式(11)构造的含未知体素单元的模式库单元，搜寻式(9)和(10) 建立的结构体矩阵模式库(图3)，匹配的模式库中单元相状态(0或1)的次数分别为：

式(12)中，读取图3中列1中对应于模式库单元(1，ζ，1，ζ)的二进制模式有1011，1010和1111，对应列2中相0的计数分别为1，2，1；相1的计数为1。

3)、根据式(12)，计算未知中心体素单元为0的概率为0.8，未知中心体素单元为1的概率为0.2。选取概率大的单元相状态(0)替代未知中心体素单元，其结果为：未知中心体素单元相状态为0。然后提出相应模式库单元的坐标信息，则可以确定未知中心体素单元取值和坐标。

重复步骤5，则可快速恢复并重构砂岩数字化矩阵。

步骤5所完成的功能是弥补原矩阵I_DM的不足，这种不足指的是在实际工程中，只获得部分区域的三维数据，对其他区域的三维数据缺失的情况下，可利用步骤5 补足这些缺失的数据，完成三维模型的重构，并最终用于安全性和稳定性分析。

重构第数字化矩阵与原矩阵I_DM中的单元取值不同，且各微小为微观结构相之间存在着一定差异，参见图4中部的黑色小方框标记处。

步骤6、利用结构体矩阵A_struct反演得到的几何拓扑结构数据进行线性插值，获得表征砂岩两相微观结构的点坐标和方向矩阵，

反演过程如下：基于建立的结构体矩阵模式库，通过计算相应单元的概率(参见步骤5中的步骤2)和步骤3)的具体说明)，重复此过程，直到与原矩阵维度相等的新矩阵生成，即完成了图像数据重建过程。从结构体中恢复相应单元的二进制模式和中心体素单元取值(即在步骤5中不断赋值，直到生成的与原矩阵维度相通的新矩阵中的各单元相状态(为0或者1)完全确定，最终完成图像重建过程，反演的结果就是获得新矩阵中各未知单元的相状态。)

线性插值的表达式依次为：

式中：P为计算体素单元的坐标；R为计算体素单元的方向矢量；Γ为等值面值(等值面值取相0和相1的平均值，为0.5)；F₀为体元中被赋值为状态0的顶点单元；F₁为体元中被赋值为状态1的顶点单元，所述体元是指由8个体素单元组成的立方体结构单元；P₀,R₀分别为几何拓扑结构数据中相为0时的任意一点坐标和方向矢量，P₁,R₁分别为几何拓扑结构数据中相为1 时的任意一点坐标和方向矢量。

根据式(13)和式(14)，计算得到表针砂岩几何拓扑结构的点坐标和方向矩阵，对其进行可视化显示，获得砂岩两相微观结构的三维数字化重构模型，如图4所示。图4(a)为原始图像获得的二值化图像(0,1图像，其中0 为黑色，1为白色)，图4(b)基于树结构的生成的二值图像，图4(c)和图4(d) 分别为本发明的二维和三维图像。

从图4(a)，(b)，(c)中部的小黑色方框可以看出：基于本发明方法中的结构体矩阵获得的结果比基于树结构获得结果更加接近原图像，说明重构的图像精度更高，能恢复更多的局部微观结构。此外，根据本方法发明，还可以恢复局部区域的数据，弥补缺失的三维数据，建立精度更好的三维数据模型，如图4(d)。

为进一步验证本方法发明的模型精度、计算效率和数据存储内存消耗，通过上述6个步骤计算了砂岩重构模型(尺寸：300²×100，500²×100，1000²×100)与真实砂岩(尺寸：1000²×100)的孔隙率、孔隙相尺寸、固体基质相尺寸和它们的相对误差见表1：

表1

由表1可以看出：本方法发明重构的模型精度与真实砂岩误差非常小。

为了将现有技术的树结构模型与本方法发明在运算上进行对比，选用了不同模板尺寸(3×3×3，5×5×5，7×7×7，9×9×9，11×11×11，13×13×13，15×15×15，17×17×17， 19×19×19)重构砂岩重构模型(尺寸：300²×100，500²×100，1000²×100)，在不同扫描模板尺寸下的CPU的计算时间和内存消耗见表2：

表2

从表2看出：本发明的CPU计算时间和数据存储消耗大大降低。因此能够在实际的工程中使用本方法发明重构深部砂岩三维模型，研究深部砂岩微观破裂-致灾机理和流体流动规律，促进深部能源开采和地下空间开发技术的发展，极大限度地降低深部岩石工程灾害造成的生命财产损失。

Claims (5)

1.一种基于结构体矩阵的砂岩三维真实模型重构方法，包括以下步骤：

步骤1、选定需建模型的砂岩块体，制作具有一定尺寸的圆柱体或长方体标准试件，运用X-射线成像设备获取砂岩标准试件几何拓扑结构图像数据；

步骤2、将获取的砂岩几何拓扑结构图像数据转换成离散3D数字化矩阵：

其特征是：

步骤3、运用Chan-Vese模型将步骤2获得的砂岩数字矩阵进行精准的两相分割，计算获得的砂岩两相数字矩阵；

步骤4、建立结构体矩阵模式库：采用3D扫描模板对砂岩两相数字化矩阵进行顺时针数据扫描，提取砂岩微观结构的二进制模式、中心体素单元和像素坐标，并将提取的数据由二进制转化为十进制存储在结构体矩阵中；

步骤5、计算中心体素单元的分布概率，利用概率密度分布函数确立未知中心体素单元相状态，补充建立的结构体矩阵模式库，确立对应的中心体素单元模式和坐标信息；

步骤6、运用结构体矩阵反演得到的几何拓扑结构数据进行线性插值，获得表征砂岩两相微观结构的点坐标和方向矩阵；

将得到的点坐标和方向矩阵进行可视化显示，获得砂岩两相微观结构的三维数字化重构模型。

2.根据权利要求1所述的砂岩三维真实模型重构方法，其特征是，在步骤2中，离散的3D数字化矩阵为：

式中：I_DM为包含砂岩几何拓扑结构特征的3D离散数字化矩阵；f为数字化矩阵体素单元；gl为矩阵单元强度值，它由像素灰度值GI决定，x，y，z为成像坐标系中在x轴、y轴和z轴方向的坐标；N为矩阵维度。

3.根据权利要求2所述的砂岩三维真实模型重构方法，其特征是，在步骤3中，所述砂岩两相数字矩阵为：

式中：

为砂岩两相的3D数字化矩阵，f_b为数字化体素单元：f_b＝0时，体素单元代表孔隙或空隙，f_b＝1时，体素单元代表固体基质。

4.根据权利要求3所述的砂岩三维真实模型重构方法，其特征是，在步骤4中，所述结构体矩阵为：

A_struct＝{XOY＝(A_Sub-pattern),XOZ,YOZ}

式中：A_struct为结构体矩阵；XOY＝(A_Sub-pattern)为在XOY平面获取的子结构体矩阵模式库A_sub-pattern，其表达式为：

A_Sub-pattern＝[S₁＝(B-pattern,Counters,Coordination ID,D-pattern,Counters,Coordination ID)；S₂；…S_M]

式中，S为一个完整模式库；M为模式库个数，分别有S₁、S₂…S_M；B-pattern、Counters、Coordination ID、D-pattern分别为二进制模式单元、中心体素单元相的计数次数、坐标信息和十进制模式；X0Z、YOZ模式库与XOY结构类型相同。

5.根据权利要求4所述的砂岩三维真实模型重构方法，其特征是，在步骤6中，所述线性插值表达式为：

式中：P为计算体素单元的坐标；R为计算体素单元的方向矢量；Γ为等值面值；F₀为体元中被赋值为状态0的顶点单元；F₁为体元中被赋值为状态1的顶点单元；P₀,R₀分别为几何拓扑结构数据中相为0时的任意一点坐标和方向矢量，P₁,R₁分别为几何拓扑结构数据中相为1时的任意一点坐标和方向矢量。

Images