一种岩层孔隙率分布计算方法
阅读说明：本技术 一种岩层孔隙率分布计算方法 (Rock stratum porosity distribution calculation method ) 是由 王少锋 皮滋滋 尹江江 周子龙 于 2021-08-02 设计创作，主要内容包括：本发明公开了一种岩层孔隙率分布计算方法,包括输出地层运动原始三彩图像；对地层运动原始三彩图像进行预处理,生成对应的三维图像数字矩阵；生成岩层裂隙层间分布图；沿岩层方向遍历岩层裂隙层间分布图,得到岩层裂隙信息提取区域图；并绘制岩层裂隙层孔隙率对应的分布曲线图；对岩层裂隙层间分布图进行全图像范围遍历,生成岩层孔隙率分布云图。本发明提供的这种岩层孔隙率分布计算方法,利用计算机对岩层图像的处理,减少了人工运算出现的遗漏或疏忽,避免了传统孔隙率计算方法中实验造成的有关误差,且在一定程度上加快了运算速度,从数据上直观呈现了岩层孔隙率在一定走向上的分布；同时减少了其他图像干扰,计算更加准确、高效。(The invention discloses a rock stratum porosity distribution calculation method which comprises the steps of outputting original three-color images of stratum movement; preprocessing the original three-color image of the stratum movement to generate a corresponding three-dimensional image digital matrix; generating a stratum fracture interlayer distribution map; traversing the stratum fracture interlayer distribution map along the stratum direction to obtain a stratum fracture information extraction area map; drawing a distribution curve chart corresponding to the porosity of the rock stratum fracture layer; and traversing the stratum fracture interlayer distribution map in a full-image range to generate a stratum porosity distribution cloud map. According to the rock stratum porosity distribution calculation method provided by the invention, the computer is used for processing the rock stratum images, so that omission or negligence caused by manual operation is reduced, related errors caused by experiments in the traditional porosity calculation method are avoided, the operation speed is accelerated to a certain extent, and the distribution of the rock stratum porosity in a certain trend is visually presented from data; meanwhile, the interference of other images is reduced, and the calculation is more accurate and efficient.)
技术领域
本发明属于岩土工程领域，具体涉及一种岩层孔隙率分布计算方法。
背景技术
如今随着基建水平的提高，地质开发工程的不断发展，在岩土工程中面临许多安全问题。例如，矿产资源开采后的踩空区容易因为维护不当而造成塌陷；工程开挖引起透水事故。上述安全问题的原因，一般是由于工程实践中岩层产生了大量裂隙，改变了地层的流体流动性质；裂隙的发育状态一般通过孔隙率进行计算。因此，对岩体工程项目中地层孔隙率分布进行研究及计算具有重要的应用意义。
岩层的孔隙率是地层中各种孔洞和裂隙的总和与岩体总体积的比值，是岩层的重要物理性质，是研究岩石动力学的重要数据资料。国内外对岩层孔隙率分布计算开展了诸多研究。一般认为的岩体孔隙率计算方法为，在忽略岩石内部孔隙率的情况下，对岩体进行冒落实验，冒落的岩石不发生变形，孔隙率可用公式表达：n＝Vf/(Vs+Vf)，n为岩层的孔隙率；Vf为岩体孔隙体积，用冒落前后岩体体积之差表示；Vs为岩体破碎冒落后的岩石总体积。除此之外，在实验室里利用岩石标本或钻井取芯进行测定、采用声波测井和中子密度测井对地层某一位置的孔隙率值进行测定。同时，现有技术中的岩层孔隙分布计算方法如：CN107402176B公开了一种裂缝孔隙度的确定方法和装置，其中，该方法包括：获取目标区域的露头样品、岩心样品；对露头样品进行三维激光扫描，建立目标区域的储层地质模型；通过目标区域的储层地质模型，确定裂缝的延伸长度；根据岩心样品，确定裂缝的开度、裂缝的充填系数、裂缝的贯穿系数；根据裂缝的延伸长度、裂缝的开度、裂缝的充填系数、裂缝的贯穿系数，确定目标区域的裂缝孔隙度。由于该方案利用表征效果较好的露头样品建立储层的地质模型，对裂缝进行定量确定；再结合岩心数据，准确地确定裂缝孔隙度，从而解决了现有方法中存在的确定的裂缝孔隙度精度不高、准确度较差的技术问题，达到了快速、准确地定量确定裂缝孔隙度的技术效果。CN105866002B发明公开了一种精确的含油页岩核磁共振孔隙度测试方法，通过测试不同等待时间、回波间隔饱和煤油页岩核磁共振T2谱分布和孔隙度，并与氦气孔隙度对比标定等待时间和回波间隔，最佳等待时间和回波间隔时核磁共振T2谱分布即可精确表征含油页岩储层物性特征。发明将影响含油页岩核磁共振孔隙度测试精度因素归结于样品饱和流体的类型和测试参数(等待时间和回波间隔)，通过饱和煤油消除页岩样品水化作用的影响；通过分析T2谱分布变化和氦气孔隙度标定最佳等待时间和回波间隔；采用所述方法，通过测试不同等待时间和回波间隔饱和煤油页岩核磁共振T2谱分布和孔隙度，即可获得能够精确表征含油页岩孔隙度的最佳测试参数。
上述现有技术中的裂隙孔隙度计算方法及核磁共振孔隙度测试方法在一定程度上均可以对地层或掩体的裂纹、孔隙度分布进行分析及计算，但技术还存在一些局限性。例如，裂隙孔隙度计算方法通过对目标区域取样进行激光扫描，可确定裂缝的延伸长度、贯穿系数等多项性质参数，但激光扫描的方法缺乏对地层全貌进行精确研究，缺乏综合性的概括，且本身对实验仪器及操作具有较高的要求；核磁共振法通过共振普分布及孔隙度与标定物的对应参数进行比较，得到页岩储层的含油饱和度和孔隙度，但不是通过直接测量手段，存在一定的误差；同时现有的技术方式较为复杂，运算速度低。
发明内容
本发明的目的在于提供一种岩层孔隙率分布计算方法，该方法能够减少计算步骤，并且应用范围广，效率高、可靠性强。
本发明提供的这种岩层孔隙率分布计算方法，包括如下步骤：
S1.对岩层裂隙富集区域进行模拟，输出对应的地层运动原始三彩图像；
S2.对地层运动原始三彩图像进行预处理，生成对应的三维图像数字矩阵；
S3.通过分析三维图像数字矩阵，去除多余岩层间像素信息，得到岩层裂隙层间分布图；
S4.沿岩层方向遍历岩层裂隙层间分布图，得到岩层裂隙代表像素点和岩层裂隙信息提取区域图；并绘制岩层裂隙层孔隙率对应的分布曲线图；
S5.设置全图像范围遍历窗口，对岩层裂隙层间分布图进行全图像范围遍历，全图像范围遍历窗口沿竖直和水平走向，生成岩层孔隙率分布云图。
所述的步骤S1具体包括：利用基于离散单元法，对岩层裂隙富集区域的地层进行模拟，对岩石基质和岩石接触进行参数化处理，生成岩层条件参数，根据岩层条件参数，输出对应的地层运动原始三彩图像image1；原彩图存储了以三种颜色[r,g,b]通道为颜色模式的像素值。
所述的岩层条件参数包括：块体的质量密度、弹性性能、塑性性能用于确定材料的力学特性；刚度、剪切刚度、抗拉强度、摩擦角、内聚力用于确定块体间接触的力学特性；地层块体尺寸。
所述的步骤S2中，预处理包括：将地层运动原始三彩图像image1生成对应的三维图像数字矩阵Array1，数字矩阵Array1采用三维形式反映三种颜色[r,g,b]通道，并分别在其对应通道以二维数组形式[[·]i]j显示像素位置，并通过像素值反映三种原彩色素值；其中，每一个像素位置包含三个颜色通道，红绿蓝三色像素值以[r,g,b]的形式，在像素位置中分布，同时r,g,b∈[0,255]。
所述的步骤S3具体包括如下步骤：
A1.去除步骤S1输出的地层运动原始三彩图像image1中的多余上网格线；
A2.去除步骤S1输出的地层运动原始三彩图像image1中的多余下网格线。所述的步骤A1具体包括：以下部网格线沿上像素行，所有列为基准，设：
w1＝I(i,j,1)/I(i,j,2)
w2＝I(i,j,2)/I(i,j,3)
其中，w1表示像素点红色通道对应的像素值与绿色通道像素值的比值；w2表示像素点绿色通道对应的像素值与蓝色通道像素值的比值；I(i,j,1)表示像素点在红通道的像素值；I(i,j,2)表示像素点在绿通道的像素值；I(i,j,3)表示像素点在蓝通道的像素值；
若满足网格线沿上设定范围，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0。
所述的步骤A2具体包括，以下部网络线沿下像素行，所有列为基准，设：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
满足网格线向下设定范围，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0；I(i,j,1)表示像素点在红通道的像素值；I(i,j,2)表示像素点在绿通道的像素值；I(i,j,3)表示像素点在蓝通道的像素值；sum表示像素点在红、绿、蓝三通道的像素值之和。
所述的步骤S4，岩层裂隙信息提取区域图image3的生成方法包括：遍历岩层裂隙层间分布图image2中的像素，像素二维数组形式为[[·]i]j，对像素点的三颜色通道像素值[r,g,b]进行求和，并设定阈值T，计算：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
其中，I(i,j,1)表示像素点在红通道的像素值；I(i,j,2)表示像素点在绿通道的像素值；I(i,j,3)表示像素点在蓝通道的像素值；sum表示像素点在红、绿、蓝三通道的像素值之和；当sum＞T时，I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝255；当sum≤T时，I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0，经过遍历和计算，像素值为0的像素点所在区域构成的图像为岩层裂隙信息提取区域图image3。
所述的步骤S4，岩层裂隙层孔隙率对应的分布曲线图image4的生成方法包括：对岩层裂隙信息提取区域图image3，设置遍历窗大小为s*k，步长取值为s，s表示s*k窗口的长，k表示s*k窗口的宽，k取定值，代表岩层裂隙信息提取区域图image3图像原宽；在一个s*k遍历窗中，岩层裂隙层孔隙率大小Pr为：
Pr＝n/N
其中，n为s*k遍历窗岩层裂隙像素点数；N为s*k遍历窗总遍历像素点数。
所述的步骤S5包括：窗口扫描岩层裂隙信息提取区域图image3，并设置遍历窗大小为m*h，步长取m；首先沿着水平遍历，再以步长为m，沿着竖直走向对岩层裂隙信息提取区域图image3进行遍历；在一个m*h遍历窗中，岩层裂隙层孔隙率大小Pr'为：
Pr'＝n'/N'
其中，n'为m*h遍历窗岩层裂隙像素点数；N'表示m*h遍历窗总遍历像素点数；将每一个m*h遍历窗中的岩层裂隙率大小Pr'进行计算，对每个窗口的孔隙率大小通过渐进颜色进行可视化，用渐进显示出不同位置遍历窗所代表孔隙率大小，得到在岩层裂隙信息提取区域图image3标准下，岩层孔隙率的分布云图image5。
本发明提供的这种岩层孔隙率分布计算方法，利用计算机对岩层图像的处理，减少了孔隙率人工运算可能出现的遗漏或疏忽，避免了因传统孔隙率计算方法中进行相关实验造成的有关误差，且在一定程度上大大加快了岩层孔隙率分布的运算速度，从数据上直观呈现了岩层孔隙率在一定走向上的分布；同时减少了其他图像干扰，计算更加准确、高效。
附图说明
图1为本发明方法的流程示意图。
图2为本发明实施例的地层运动原始图像示意图。
图3为本发明实施例的岩层裂隙层间分布图示意图。
图4为本发明实施例的岩层裂隙信息提取区域图示意图。
图5为本发明实施例的岩层裂隙分布曲线图。
图6为本发明实施例的岩层孔隙率的分布云图示意图。
具体实施方式
如图1为本发明方法的流程示意图：本发明提供的这种岩层孔隙率分布计算方法，包括如下步骤：
S1.对岩层裂隙富集区域进行模拟，输出对应的地层运动原始三彩图像；
步骤S1具体为，利用基于离散单元法，对岩层裂隙富集区域的地层进行模拟，根据岩层条件参数，输出对应的地层运动原始三彩图像image1；地层运动原始三彩图像image1凸显了岩层裂缝扩张的定性表现，如裂隙发育、分布的位置及裂隙的密度参数等；原彩图存储了以三种颜色[r,g,b]通道为颜色模式的像素值。
其中，在利用基于离散单元法的UDEC数值软件对地层运动进行模拟时，对岩石基质(块)和岩石接触(岩层中的节理和岩层)进行参数化处理，生成岩层条件参数，岩层条件参数包括：块体的质量密度、弹性性能、塑性性能等以确定材料的力学特性；刚度、剪切刚度、抗拉强度、摩擦角、内聚力等以确定块体间接触的力学特性；地层块体尺寸等。
S2.对地层运动原始三彩图像进行预处理；
预处理过程具体包括：将地层运动原始三彩图像image1生成对应的三维图像数字矩阵Array1，数字矩阵Array1采用三维形式反映三种颜色[r,g,b]通道，并分别在其对应通道以二维数组形式[[·]i]j显示像素位置，并通过像素值反映三种原彩色素值。
其中，每一个像素位置包含三个颜色通道，红绿蓝三色像素值以[r,g,b]的形式，在像素位置中分布，同时r,g,b∈[0,255]。
S3.通过分析三维图像数字矩阵Array1，对地层运动原始三彩图像image1进行检测，对岩层裂隙边界位置信息进行提取，去除多余岩层间像素信息，得到岩层裂隙层间分布图image2；
对岩层裂隙边界位置信息进行提取的具体方法包括如下步骤，去除步骤S1输出的地层运动原始三彩图像image1中的多余上下网格线，多余上下网格线表示多余层间信息，使裂隙层间信息更加突出，减少了计算机计算负担，加快运算速度。
A1.以下部网格线沿上像素行，所有列为基准，设：
w1＝I(i,j,1)/I(i,j,2)
w2＝I(i,j,2)/I(i,j,3)
其中，，w1表示像素点I[i,j]红色通道对应的像素值与绿色通道像素值的比值；w2表示像素点I[i,j]绿色通道对应的像素值与蓝色通道像素值的比值。I(i,j,1)表示像素点在红通道的像素值；I(i,j,2)表示像素点在绿通道的像素值；I(i,j,3)表示像素点在蓝通道的像素值。
满足网格线沿上设定范围，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0；在本实施例中网格线沿上设定范围为0.9＜w1,w2＜1.1，即0.9和1.1分别为本实施例中像素通道值对应比值的上下界；
A2.以下部网络线沿下像素行，所有列为基准，设：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
满足网格线向下设定范围，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0；在本实施例中网格线向下设定范围为130＜sum＜170；sum表示像素点在红、绿、蓝三通道的像素值之和。
通过步骤A1和A2实现了层间裂隙特征信息的提取，得到岩层裂隙层间分布图image2。
S4.沿岩层方向遍历岩层裂隙层间分布图image2，设定阈值T，T值根据岩层裂隙层间分布图image2裂隙像素点实际值选取，在本实施例中选取8；对RGB值之和Isum小于设定阈值T的像素点进行计算、识别和提取，得到岩层裂隙代表像素点gray[i,j]，以及岩层裂隙信息提取区域图image3；将岩层裂隙代表像素点gray[i,j]进行求和，得到遍历窗岩层裂隙像素点数n1，并计算岩层裂隙层孔隙率大小；并绘制岩层裂隙层孔隙率对应的分布曲线图image4；
步骤S4中，生成岩层裂隙信息提取区域图image3的具体方法为：遍历岩层裂隙层间分布图image2中的像素，像素二维数组形式为[[·]i]j，对像素点的三颜色通道像素值[r,g,b]进行求和，并设定阈值T，计算：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
当sum＞T时，I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝255；当sum≤T时，I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0，经过遍历和计算，像素值为0的像素点所在区域构成的图像为岩层裂隙信息提取区域图image3。
绘制岩层裂隙层孔隙率对应的分布曲线图image4的具体方法包括：对岩层裂隙信息提取区域图image3，设置遍历窗大小为s*k，步长取值为s，s表示s*k窗口的长，k表示s*k窗口的宽，k取定值，代表岩层裂隙信息提取区域图image3图像原宽；在一个s*k遍历窗中，岩层裂隙层孔隙率大小Pr为：
Pr＝n/N
其中，n为s*k遍历窗岩层裂隙像素点数；N为s*k遍历窗总遍历像素点数。并通过计算机程序进行可视化，由此可得到在水平方向上单位长度的层裂隙分布曲线图image4。
S5.设置全图像遍历窗口，对岩层裂隙层间分布图image2进行全图像范围遍历，全图像范围遍历窗口沿竖直和水平走向，按照设置步长m逐像素移动窗口；对岩层裂隙层间分布图image2完成全图像范围遍历后，可视为岩层孔隙率分布云图image5。
步骤S5具体包括，窗口扫描岩层裂隙信息提取区域图image3，并设置遍历窗大小为m*h；步长取m，在本实施例中m＝1；首先沿着水平遍历，再以步长为m，沿着竖直走向对岩层裂隙信息提取区域图image3进行遍历；在一个m*h遍历窗中，岩层裂隙层孔隙率大小Pr'为：
Pr'＝n'/N'
其中，n'为m*h遍历窗岩层裂隙像素点数；N'表示m*h遍历窗总遍历像素点数；将每一个m*h遍历窗中的岩层裂隙率大小Pr'进行计算，对每个窗口的孔隙率大小通过渐进颜色进行可视化，用渐进显示出不同位置遍历窗所代表孔隙率大小，由此得到在岩层裂隙信息提取区域图image3标准下，岩层孔隙率的分布云图image5。
在本实施例中，以某煤矿岩层作为本发明方法应用实例。山西某煤矿采用地下长壁采煤法对某煤层进行开采，该煤矿集中在最接近地表的4号煤层。4号煤层厚度为6.5m，采煤机截割高度为3.0m，放顶煤厚度为3.5m。长壁开采的宽度为200米。4号煤层及其上覆地层地质条件相对简单，无断层或水文影响，煤层及上覆地层结构基本为水平结构。
方法应用步骤如下：
步骤1、如图2为本发明实施例的地层运动原始图像示意图。利用基于离散单元法的UDEC数值软件对岩层裂隙富集区域的地层运动进行模拟，在UDEC模型中，第2关键层(主关键层)和第1关键层的块体尺寸分别设计为42.6m和23.5m。同时获得了4号煤层以上煤系地层的实测岩性和几何性质，以及地层中完整岩石基质的力学性质。输出对应的地层运动原始图像如图2，图像显示了岩层裂缝发育的定性表现，如裂隙分布的位置和裂隙的密度参数；以及RGB三通道的色素值。
步骤2、对计算机程序输入1步骤中的地层运动原始图像image1，读取并存储原图，生成对应的三维图像数字矩阵Array1，数字矩阵Array1在三维中反映三种颜色通道：r、g、b，并分别在其对应通道以二维数组形式[[·]i]j形式显示像素位置，并通过像素值反映红-绿-蓝三种原彩色素值；
步骤3、如图3为本发明实施例的岩层裂隙层间分布图示意图。通过分析步骤2步骤中地层运动原始图像矩阵Array，对地层运动原始图像及其岩层裂隙边界位置信息进行提取的具体方法为：去除步骤1输出的地层运动原始图像中的多余上下网格线(表示多余层间信息)，目的在于使裂隙层间信息更加突出，且减少计算机计算负担，加快运算速度。以下部网格线沿上像素行，所有列为基准，设：
w1＝I(i,j,1)/I(i,j,2)
w2＝I(i,j,2)/I(i,j,3)
其中，w1表示像素点红色通道对应的像素值与绿色通道像素值的比值；w2表示像素点绿色通道对应的像素值与蓝色通道像素值的比值；I(i,j,1)表示像素点在红通道的像素值；I(i,j,2)表示像素点在绿通道的像素值；I(i,j,3)表示像素点在蓝通道的像素值；
满足0.9＜w1,w2＜1.1，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0；
以下部网络线沿下像素行，所有列为基准，设：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
满足130＜sum＜170，则令I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0。
经过上述算式实现了层间裂隙特征信息的提取，得到岩层裂隙层间分布图。
步骤4、遍历扫描步骤3中的岩层裂隙层间分布图像素，像素二维数组形式为[[·]i]j，并对像素的三颜色通道像素值[r,g,b]进行求和，并设定阈值T，T值根据岩层裂隙层间分布图裂隙像素点实际值选取，本图片取为8：
sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)
其中，当sum＞T时，设置I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝255；当sum≤T时，设置I(i,j,1)＝I(i,j,2)＝I(i,j,3)＝0。如图4为本发明实施例的岩层裂隙信息提取区域图示意图。经过遍历与计算，像素值为0的像素点所在区域构成的图像即为岩层裂隙信息提取区域图。
步骤4中，如图5为本发明实施例的岩层裂隙分布曲线图。绘制岩层裂隙分布曲线的方法为，对岩层裂隙信息提取区域图设置遍历窗大小为2*10，步长取2。在一个遍历窗(2*10)中，岩层裂隙层孔隙率大小Pr为Pr＝n/N；
其中，n为遍历窗(2*10)岩层裂隙像素点数，N为遍历窗(2*10)总遍历像素点数。将数据结果通过计算机程序进行可视化，由此可得到如图5所示的在水平方向上单位长度中岩层裂隙分布曲线图。
步骤5、窗口扫描图4所示的岩层裂隙信息提取区域图，按照步骤4提取岩层裂隙区特征，然后设置遍历窗大小为5*10，步长取1，首先沿着水平，再以步长为1，沿着竖直走向对图4所示的岩层裂隙信息提取区域进行遍历。在一个遍历窗(5*10)中，岩层裂隙层孔隙率大小Pr'为：
Pr'＝n'/N'
其中，n'为遍历窗(5*10)岩层裂隙像素点数，N'为遍历窗(5*10)总遍历像素点数。如图6为本发明实施例的岩层孔隙率的分布云图示意图。将每一个遍历窗中的岩层裂隙层孔隙率大小进行计算，并通过计算机程序，对每个窗口的孔隙率大小通过渐进颜色进行可视化，用渐进显示出不同位置遍历窗所代表孔隙率大小，由此得到在岩层裂隙信息提取区域图标准下，如图6所示的岩层孔隙率的分布云图。
对该煤矿岩层4号煤层运用本方法分析计算岩层孔隙率分布。计算结果如图2～6所示。
从图2～6可知，在该煤矿岩层4号煤层运用一种基于数字图像识别的岩层孔隙率分布计算方法，成功对该煤岩层进行UEDC地层模拟、数字矩阵提取、岩层裂隙边界信息提取及岩层裂隙层间分布图获取、岩层裂隙信息提取区域图获取、孔隙率分布曲线图及云图的绘制。其中，通过两种对孔隙率的计算的可视化结果——孔隙率分布曲线图及云图表明，该方法可以高效快速地解算出岩层孔隙率，且整体运算实施过程符合工作岩体实际情况。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：基于深度学习的色选机分选方法