基于色差识别的岩石微观结构射线成像探测法

1.本发明属于岩石微观结构图像处理技术，具体涉及一种基于色差识别的岩石微观结构射线成像探测法。


背景技术：

2.岩石微观结构对各种岩石工程灾害预测和防控有重要影响，而微观破岩 (岩石在不同加载条件下的岩石损伤破裂)机理和流体在岩石内部微观结构的流动规律对深部能源开采和地下空间开发具有十分重要的用途。在实际工程应用中(如页岩气、石油开发项目)，现有的岩石微观结构探测方法主要有流体示踪探测法、激光和红外线探测法和射线成像探测法等几类。
3.属于射线成像探测法的x-射线层析成像方法，它的微观结构探测过程包括以下步骤：步骤1、选取材料制作岩石标准试件，并运用x-射线成像设备获取岩石多微观结构特征信息。步骤2、选取适当分割算法，提取并分割岩石内部各个微观结构相；如文献“a threshold section method from gray-levelhistogram”，otsu n，1979.iee trans smc.9(1),62-66.(“一种基于灰度直方图的阈值分割方法”，大津，1979.ieee系统、人与控制处理,9(1),62-66) 公开的两相分割算法，它通过类间方差计算得到一个分割阈值来分割岩石。步骤3、确定各个微观结构相赋值区间，运用计算机进行数据分析和结果可视化。
4.现有的岩石射线成像探测法对不同岩石微观结构都有较好的识别分割效果，但存在以下的问题：
5.1、岩石微观结构分割缺乏统一标准，不同分割算法针对同一岩石的微观结构分割结果可能相差甚远。
6.2、岩石微观结构分割结果易受到成像因素(如分辨率等)和分割算法的影响，容易造成微观结构欠分割或者过度分割，微观结构分割精度受到一定限制。
7.3、针对带裂隙的岩石，无法准确分类识别裂纹和孔隙，易将裂纹和孔隙均当作孔隙，导致岩石微观结构分类识别精度低，例如两相分割法用于复杂岩石或包含其他流体相的岩石，其分割精度较低。
8.本发明针对现有的岩石微观结构射线成像探测方法的分割算法进行改进，以解决上述的问题。


技术实现要素：

9.针对现有的岩石微观结构探测方法的岩石多微观结构分类识别不清，无法准确分类识别裂纹和孔隙，微观结构过度分割或欠分割等问题，本发明要解决的问题是提供一种基于色差识别的岩石微观结构射线成像探测法，它能够精确识别岩石裂纹和孔隙及其他固体基质结构相，提高岩石微观结构识别分割精度，避免了由岩石微观结构相过度分割或欠分割导致预测的岩石物理力学特性误差较高的问题。
10.本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：
11.步骤1、选择指定的砂岩块体，制作具有一定尺寸的长方体砂岩标准试件；在不同加载状态下，采用x-射线成像设备获取试砂岩试件破坏过程的分层图像；将砂岩分层图像离散并转化生成二维(2d)标量矩阵；
12.步骤2、导入转换的砂岩2d标量矩阵，采用三通道演示转换函数将矩阵各个单元添加颜色标签，生成色差矩阵单元；运用色差单元高斯分布初步估算矩阵单元识别阈值；基于各个色差分块矩阵边界单元色差，获得划分砂岩内部各个结构相的色差识别阈值；
13.步骤3、导入步骤2建立的色差矩阵，提取剔除颜色标签之后的矩阵单元，建立矩阵单元的灰阶矩阵，计算灰阶矩阵分块区间的识别阈值，
14.步骤4、将步骤2和步骤3得到的各个结构相的色差识别阈值和灰阶识别阈值取平均值，获得最优的砂岩各个微观结构相的最优化分割阈值；
15.步骤5、运用获得的多个最优化分割阈值，精确快速分割砂岩内部各个微观结构相；通过将各个微观结构相的矩阵单元赋予0-1之间不同的值，进行数字化表征并实现结果可视化。
16.本发明根据色差识别不同岩石成分，划分不同单元区间，从而统一分割标准如岩石中的石英、空隙等均有固定的阈值区间。根据添加的单元色差，能够准确区分并且提取不同岩石成分，从而避免欠分割或过度分割。
17.本发明的技术效果是：
18.本发明在步骤2增设了矩阵单元的色差识别分割算法，该算法基于像素单元，引入颜色标签组成色差单元矩阵，通过色差分块矩阵边界单元色差，计算岩石内部各微观结构相色差识别阈值，并结合步骤3计算的灰阶识别阈值，获得岩石微观结构相最优分割阈值，达到精确分割岩石内部多类微观结构相的目的。本发明实现了精确识别区分岩石裂纹和孔隙及其他固体基质结构相，提高了岩石微观结构识别分割精度，提高了岩石物理力学参数的预测精度。
附图说明
19.本发明的附图说明如下：
20.图1为本发明的最优阈值计算图；
21.图2为一个实施例分割和提取单独相的成像图。
具体实施方式
22.下面结合附图和实施例对本发明作进一步说明：
23.本实施例选单轴压缩破坏时的重庆砂岩为对象，其岩石内部各微观结构相的分割按以下步骤：
24.步骤1、选择指定的砂岩块体，制作具有一定尺寸的长方体砂岩标准试件，当不同加载状态下导致砂岩试件破坏时，采用x-射线成像设备获取试砂岩试件破坏过程的分层图像；将砂岩分层图像离散并转化生成二维(2d)标量矩阵，它的表达式为：
[0025][0026]
式(1)中：i
dm
为包含砂岩几何拓扑结构特征的2d离散数字化矩阵；f为数字化矩阵体素单元，所述体素单元等同于像素单元，均代表成像坐标系中的最小单元；gl为矩阵单元强度值，它由像素灰度值gi决定(0-255)；x、y为成像坐标系中在x轴和y轴方向的坐标，n为矩阵维度。
[0027]
x-射线分层图形是通过x-射线接收平面形成的三维图像矩阵，这个三维图像矩阵通常由一层一层的二维分层图像组成。因此，x-射限分层图形是一个三维立体图，任意选取一层二维分层图像，就可以用公式(1)表示包含三维特征信息的二维分层图像。式(1)仅为组成三维图像矩阵中的一个二维矩阵的数学表达式。
[0028]
步骤2、矩阵单元色差识别
[0029]
1)、导入转换的砂岩2d标量矩阵，采用三通道演示转换函数将矩阵各个单元添加颜色标签，生成色差矩阵单元。
[0030]
文献“digitalmeasurementof2dand3dcracksinsandstonesthroughimprovedpseudocolorimageenhancementand3dreconstructionmethod.”zhaozhi,zhouxiaoping，internationaljournalfornumericalandanalyticalmethodsingeomechanics.2019,43,(16):2565-2584.(“基于非彩色图像增强和三维重构算法的砂岩二维和三维裂纹识别”，赵智，周小平.地质力学数值与分析方法国际期刊.2019,43,(16):2565-2584.)记载了三通道演示转换函数的表达式为：
[0031][0032]
式(2)中：为包含砂岩各微观结构相的2d色差矩阵；f
cd
为带颜色标签的矩阵色差单元；γ
gi
→
rgb
为三色彩通道转换函数，γ
gi
→r为红色通道转换函数；γ
gi
→g绿红色通道转换函数；γ
gi
→b为蓝色通道转换函数。转换函数将岩gi转换成不同通道颜色值。
[0033]
通过式(2)的非线性计算得到的红、绿、蓝三通道颜色值，再赋予矩阵各个单元的色差值，具体过程为：将红、绿和蓝三通道颜色值分别按照一定比例(这个比例可根据编程语言中自带的调色模板设定)分配，如0.5,0.2,0.3进行线性叠加，然后向上取整数获得一
个标签值，这个标签值与单元强度组组合形成色差值，例如153|223,122,87，其中，153为单元强度值，233,122,87分别为红、绿和蓝色通道值，该单元的标签值即为233
×
0.5 122
×
0.2 87
×
0.3＝167，色差值则为153|167；色差值即是矩阵相应单元的强度值和标签值的组合，用于在图像可视化中突出显示各个微观结构岩色差异。
[0034]
2)、计算色差矩阵单元强度高斯分布，运用色差单元高斯分布初步估算矩阵单元识别阈值，它的表达式为：
[0035][0036]
式(3)中：f为数字化矩阵体素单元，为第i个子矩阵分块与第i 1个子矩阵分块之间的识别阈值，用于划分两个子矩阵。d是微分，i为子矩阵分块的序号，i

,i-分别代表微分右极限和左极限；pn为矩阵单元识别阈值个数。
[0037]
3)、基于各个色差分块矩阵边界单元色差，获得划分砂岩内部各个结构相的色差识别阈值，它的表达式为：
[0038][0039]
式(4)中：fk、f
k 1
分别为色差识别过程中搜寻到的第k个和k 1个单元；分别为带颜色标签的第i个和第i 1个矩阵色差单元；为第i个子矩阵相与第i 1个子矩阵相之间的色差识别阈值；p
nn
为色差识别阈值数目，代表可识别p
nn
1个子矩阵相。
[0040]
说明：1、矩阵分块是指纯数学矩阵，矩阵相是指表征不同结构相的矩阵单元形成的矩阵，矩阵分块可包含矩阵相，且各矩阵相之间是由矩阵边界单元区分的。
[0041]
2、为第i个子矩阵分块与第i 1个子矩阵分块之间的识别阈值，用于划分两个子矩阵；为第i个子矩阵相与第i 1个子矩阵相之间的色差识别阈值。简单地说，它们的本质不同，大小略有差异。当矩阵分块中只有一个矩阵相时，与的值可能接近，但前者代表是一个粗糙的识别阈值(即识别精度较低)，后者代表一个相对中等的识别阈值(即识别精度相对有所提高)。
[0042]
步骤3、导入步骤2建立的色差矩阵提取剔除颜色标签之后的矩阵单元，建立灰阶矩阵，计算灰阶矩阵分块区间的识别阈值。
[0043]
剔除颜色标签的举例：色差矩阵中的某一单元，它可以表达为153|223,122,87，153表示矩阵单元强度值，223,122,87表示组成颜色标签的红、绿、蓝三色通道值。剔除颜色标签即取这4个值的平均，并向上取整数，代替矩阵单元的强度值。
[0044]
文献“fracturecharacterizationandpermeabilitypredictionbyporescalevariablesextractedfromx-rayctimagesofporousgeomaterials”，zhao
zhi, zhou xiao-ping,qian qi-hu，science china technological science. 2020,63(5):755-767(“基于x-ray ct图像提取孔隙尺度变量的多孔岩土材料裂纹特征描述和渗透性预测”，赵智，周小平，钱七虎，中国科学-技术科学， 2020，63(5)：755-767)记载：灰阶矩阵分块区间的识别阈值的表达式为：
[0045][0046]
式(5)中：
[0047][0048][0049]
式(5)-(7)中，为第i个子灰阶矩阵与第i 1个子灰阶矩阵之间的识别阈值；p
ns
为灰阶矩阵识别阈值个数；fi为第i次计算灰阶阈值过程中将灰阶矩阵划分为两个子灰阶矩阵的灰阶；为搜寻算符；λ
max
为最大值算法；λ
abs
为绝对值算符；χ为方差；
[0050]ccount
为技术算符，表示统计某一项的个数；下标i→1和i
→
2是指当第i次计算灰阶阈值过程中将灰阶矩阵划分为两个子灰阶矩阵之后分别对应于第1个和第2个矩阵；为移除算符，表示去除某个子矩阵中的其他矩阵中的单元。
[0051]
说明：f为数字化矩阵体素单元，它表示任意一个初始图像矩阵单元； f
cd
为将f经过公式(2)色彩转换后图像矩阵中的任意一个色彩单元，两者取值不同。本步骤的fi为第i次计算灰阶阈值过程中将灰阶矩阵划分为两个子灰阶矩阵的灰阶，它只是一个单元值。
[0052]
步骤4、将步骤2和步骤3得到的各个结构相的色差识别阈值和灰阶识别阈值取平均值，就可以获得最优的砂岩各个微观结构相的最优化分割阈值，它的表达式为：
[0053][0054]
式(8)中，为第i个微观结构相与第i 1个微观结构相之间的最优化分割阈值；λ
average
为平均值算符。
[0055]
图1为最优阈值计算示意图，包括计算得到的最优阈值和单元区间划分示意图。纵轴表为矩阵单元强度相对频率(相对频率指矩阵中具有相同值的各类单元数目与矩阵单元总数的比例)与矩阵单元强度的分布，曲线经过一般的频率直方图统计得到。图2中，由公式(8)计算得到的分割阈值，用于划分裂纹相、空隙相和固体基质相。本实例
中，只有裂纹相、空隙相和固体基质相存在，因此只有和对于其它实例，在存在更多微观结构相时，本发明方法仍可计算出其他相分割阈值，则可能会有和甚至更多分割阈值ts。
[0056]
步骤5、将岩石微观结构相矩阵各个单元赋予新值0～1，进行数字化表征，实现结果可视化。
[0057]
运用式(8)获得的岩石微观结构相(指组成岩石微观几何拓扑结构的不同结构相，一般包括孔隙、裂纹和固体基质相等)各最优分割阈值，将矩阵单元(该矩阵是与岩石图像矩阵具有相同维度的空矩阵)依次赋值，通过依次对不同微观结构相赋值而生成岩石微观结构相矩阵。岩石内部各微观结构相矩阵的表达式为：
[0058][0059]
式(9)中：为分割之后的砂岩微观结构相矩阵，f
mul-phase
为各微观结构单元属性值；i为识别的最优阈值序号，它取值1～p
no
；p
no
为识别的最优阈值个数。
[0060]
本发明通过上述步骤，将i
dm
数字化矩阵、色差矩阵和步骤3的灰阶矩阵，优化生成岩石微观结构相矩阵进行可视化。
[0061]
应用式(9)例如：若i＝1，识别阈值个数为p
no
＝1，可划分两个子分块矩阵，则两个矩阵单元分别赋值0和1/p
no
＝1，i＝2,3
……
时，依次类推。采用式(9)赋值的效果是：1、矩阵单元强度值归一化，将不同的微观结构相归一化后赋予不同值，同一微观结构相的单元取值与其他微观结构相的差异增大，因此可突出显示，能够突出显示所需分析对象可视化结果；2、归一化处理后，同一微观结构相单元取值相同，可单独提取某一微观结构相，按不同微观相提取对应的特征值，便于快速进行特征分析。相对应地，未经归一化处理的各类微观结构相单元取值相近，且同一微观结构相单元取值又不相同，区分难度大。
[0062]
运用步骤4获得的多个最优化分割阈值，就可以精确快速分割砂岩内部各个微观结构相。本发明方法通过归一化，加强微观结构相的特征，可突出识别分割所需研究对象，在本实例中，能精确分割提取裂纹相和孔隙相。如需要研究其他微观结构相时，本发明还可自动计算出其他相分割阈值，提高分割精度。而现有的两相分割算法则无法单独识别并分别裂纹和孔隙相，且对多相分割识别精度低。
[0063]
在图2中的三图为本发明单独提取某一相的特征图。从图2可以看出：单独提取并分割孔隙相和裂纹相，相对于现有技术，本发明方法能够单独提取某一相特征进行分析，且分割精度更高。