一种基于三维成像技术的岩土体渗透率计算方法与流程

本发明涉及渗透率技术领域，尤其涉及一种基于三维成像技术的岩土体渗透率计算方法。

背景技术：

目前，对于工程地质和岩土工程，常规岩土体材料的渗透率测试主要以宏观试验为主，通过对取到的岩芯进行气体稳态法或瞬态法试验，使用达西定律计算得到渗透率大小。这种宏观试验的好处是能够模拟各种应力应变情况，较为真实的模拟各种地质条件下的变化。但缺点也很明显：1、实验需要一整套设备，造价昂贵，操作复杂；2、实验周期长，有的试样渗透率低，实验周期长达几个月；3、实验过程中，由于需要施加一定的围压，存在对试样造成损伤和破坏的风险，对于一些不易取得，或者说在取芯难度和成本较大的情况下，岩芯不易进行宏观实验。即便取得了岩心，宏观实验也会造成一定程度的损伤，无法再次进行后续的分析测试，实验之后岩心便失去了最初的价值。

岩土体材料的渗透率与其内部微观结构息息相关，而岩心内部含有大量跨尺度分布、形态复杂的孔/裂隙，从毫米级的裂缝到纳米级连氮分子都无法通过的微孔隙，使得岩心内部微观结构异常复杂、无序。近年来，随着微细观成像技术的发展，研究人员可以对试样内部的孔隙结构进行直接的观测。通过对微观图像直接的分析，可以将微观岩土结构特征与流体渗流特性之间建立联系，从而来解释宏观渗流实验的一些现象，更加深刻地理解岩土体中流体的迁移规律以及其细观机理。

目前，较为常用的岩土体微观三维成像技术有：计算机断层成像(computedtomography，简称ct)、聚焦离子束(focusedionbeam，简称fib)等。这些成像技术的原理略有不同，但都是通过扫描收集样品的表面测绘学形貌和组成信息的各种信息，创建显示微结构的表面形貌图像。通过对岩土体试样进行一层一层的扫描，创建每一层的岩土体内部微结构形貌图像，使用技术手段将所有断层图像进行处理，三维重构出岩土体的内部孔隙结构立体形貌。

技术实现要素：

发明目的：针对现有岩土体材料渗透率测试手段单一，测试设备昂贵，测试周期长以及无法弄清楚岩土体材料的微观渗流机理等问题，本发明提出一种基于三维成像技术的岩土体渗透率计算方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：

一种基于三维成像技术的岩土体渗透率计算方法，所述方法包括如下步骤：

s1：读取n张岩土体材料的断层微观图像，提取所述n张断层微观图像的孔隙结构，获取所述n张断层微观图像的二值图；

s2：将相邻两个连续的所述断层微观图像的二值图进行合成重构，获取n-1张合成二值图；

s3：根据hagen-poiseuille定律和所述合成二值图，获取每一张合成二值图的渗透率；

s4：根据所述合成二值图的渗透率，确定所述岩土体渗透率，具体为：

其中：n为相邻断层图像的总数目，ki为合成二值图的渗透率。

进一步地讲，所述步骤s1获取n张断层微观图像的二值图，具体如下：

s1.1：读取n张岩土体材料的断层微观图像，获取断层微观图像中每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni；

s1.2：根据所述断层微观图像中每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni，获取岩土体材料灰度图像的灰度直方图曲线上的各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数；

s1.3：根据所述各个p(i)点，拟合获取岩土体材料灰度图像的灰度直方图曲线，提取所述n张断层微观图像的孔隙结构，获取所述n张断层微观图像的二值图。

进一步地讲，所述步骤s2获取n-1张合成二值图，具体如下：

s2.1：读取所述不同断层微观图像的二值图，获取不同断层微观图像的二值图中所有像素点的数值大小；

s2.2：根据所述不同断层微观图像的二值图中所有像素点的数值大小，确定不同断层微观图像的二值图对应的矩阵大小；

s2.3：将所述不同断层微观图像的二值图对应的矩阵大小进行相乘，确定合成矩阵的大小；

s2.4：根据所述合成矩阵的大小，获取n-1张所述合成二值图。

进一步地讲，所述不同断层微观图像的二值图中所有像素点的数值大小具体为：所述不同断层微观图像的二值图中白色像素点的数值大小为0、黑色像素点的数值大小为1。

进一步地讲，所述步骤s3获取每一张合成二值图的渗透率，具体如下：

s3.1：根据所述合成二值图，获取合成二值图中连通孔隙管道的数目，以及所述连通孔隙管道对应的圆管半径大小；

s3.2：根据所述hagen-poiseuille定律、合成二值图中连通孔隙管道的数目和连通孔隙管道对应的圆管半径大小，获取每一张合成二值图的渗透率，具体为：

其中：a为流体的截面面积，ri为孔隙半径，ai为孔隙的面积。

进一步地讲，所述步骤s3.1根据合成二值图，获取合成二值图中连通孔隙管道的数目，以及连通孔隙管道对应的圆管半径大小，具体如下：

s3.1.1：对所述合成二值图中的每一孔隙对象进行计算，获取每一孔隙对象的等距图，根据所述每一孔隙对象的等距图，确定等距图中欧几里德距离最大的像素点；

s3.1.2：以所述欧几里德距离最大的像素点为圆心，以所述像素点到孔隙和基质边界的距离为半径画圆，逐渐缩小欧几里德距离直至为0，并重复步骤s3.1.1至s3.1.2，获取所述孔隙面积随孔径的变化关系，确定所述合成二值图中连通孔隙管道的数目，以及连通孔隙管道对应的圆管半径大小。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

(1)本发明通过对岩土体材料的数字化重构及分析，可以减少由于分析所需的时间过长，而导致方案的延期，也可以避免由于根据经验数据或者类比数据决策，而增加方案的不确定性，进而更加快速的提供渗透率结果；

(2)本发明可以揭示岩土体材料的内部物理现象，从微观角度出发，解释宏观现象，可以从根本上解释并解决岩土体材料的渗流问题；

(3)本发明可以对数字岩心反复利用，测试各种数据，提高了数字岩心的利用率，减少了资源的损失，降低了试验成本。

附图说明

图1是本发明岩土体渗透率计算方法的整体流程示意图；

图2是本发明使用fib成像技术的岩土体微观图像；

图3是本发明第n张fib图像及对应二值化图像；

图4是本发明第n+1张fib图像及对应二值化图像；

图5是本发明二值化图像像素信息提取示意图；

图6是本发明合成二值图；

图7是本发明合成二值图中连通孔隙管道示意图；

图8是本发明相邻图像之间连通孔隙管道的孔径分布及大小；

图9是本发明整个岩土体试样通过断层扫描图像三维重构后示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。

实施例1

参考图1，本实施例提供了一种基于三维成像技术的岩土体渗透率计算方法，该岩土体渗透率计算方法可以应用于三维成像技术得到的任何图像。但是在岩土体渗透率的计算过程中，计算的是与岩土体材料不发生反应的流体渗透率，其中包括气体、水等。更进一步地讲，岩土体材料吸水后的力学性质及孔隙结构会发生改变的岩土体渗透率不适合用本发明的岩土体渗透率计算方法进行求取。在本实施例中，以岩土体材料的气体渗透率的计算方法为例进行说明，其中岩土体渗透率计算方法具体包括如下步骤：

步骤s1：通过matlab代码读取n张岩土体材料的断层微观图像，并提取n张断层微观图像的孔隙结构，获取n张断层微观图像的二值图，具体如下：

参考图2，使用三维成像技术获取岩土体试样大量连续的断层微观图像，在本实施例中，岩土体试样大量连续的断层微观图像为n张，其中每一张图像，均可以直观的看到包含孔隙的微观结构。

步骤s1.1：通过matalb代码读取待测岩土体材料的n张断层微观图像，获取该断层微观图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。

步骤s1.2：根据断层微观图像中每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni，通过如下公式获取待测岩土体材料灰度图像的灰度直方图曲线上的各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数。

步骤s1.3：根据各个p(i)点，拟合获取待测岩土体材料的灰度图像的灰度直方图曲线，并根据该曲线中的峰值数目确定待测岩土体材料的结构，进而提取出断层微观图像的孔隙结构，获取断层微观图像的二值图，该断层微观图像的二值图为n张。具体地讲，在断层微观图像的二值图中，黑色像素表示的是岩土体孔隙结构，白色像素表示的是岩土体表面土骨架。

步骤s2：参考图3和图4，将不同断层微观图像的二值图进行合成重构，获取n-1张合成二值图，其中需要进行合成重构的不同断层微观图像的二值图为相邻的两个连续的断层微观图像的二值图，这是由于两张连续的断层微观图像的二值图之间的距离较小，符合hagen-poiseuille定律的直管流体计算条件。如果是间隔的两张断层微观图像的二值图，其之间的距离将会变大，从而孔隙之间的连通管道将不能近似为直管，从而不符合hagen-poiseuille定律的直管流体计算条件。其具体过程如下：

步骤s2.1：通过matlab代码读取不同断层微观图像的二值图，获取不同断层微观图像的二值图中所有像素点的数值大小，即读取两张连续的断层微观图像的二值图，获取两张连续的断层微观图像的二值图中所有像素点的数值大小。

其中matlab代码读取断层微观图像的二值图，获取的断层微观图像的二值图中所有像素点的数值大小分为两种，分别为：白色像素点的数值大小为0，黑色像素点的数值大小为1。

步骤2.2：参考图5，根据不同断层微观图像的二值图中所有像素点的数值大小，确定不同断层微观图像的二值图对应的矩阵大小。通过获取得到的断层微观图像的二值图中所有像素点的数值大小，可以得到两张连续的断层微观图像的二值图中所有像素点的数值信息，即两个只包含0和1的矩阵，在本实施例中，第一个断层微观图像的二值图中所有像素点的数值信息对应的矩阵为a，连续的第二个断层微观图像的二值图中所有像素点的数值信息对应的矩阵为b。

步骤s2.3：将不同断层微观图像的二值图对应的矩阵大小进行相乘，确定合成矩阵的大小，即将矩阵a和矩阵b进行相乘，获取得到合成矩阵c，具体为：

c＝a·b

其中：c为合成矩阵，a为第一个断层微观图像的二值图中所有像素点的数值信息对应的矩阵，b为第二个断层微观图像的二值图中所有像素点的数值信息对应的矩阵。

其中矩阵a和矩阵b相乘的过程中，如果矩阵a和矩阵b中相同位置处的像素若均为1，则表示该处为相互连通的像素，此时积为1，即在矩阵c中相同位置处的像素也为1。如果矩阵a和矩阵b中相同位置处的像素若含有0，则表示该处不为相互连通的像素，即表示岩土体基质，此时积为0，即在矩阵c中相同位置处的像素也为0。

步骤s2.4：参考图6和图7，根据合成矩阵c的大小，获取n-1张合成二值图。即通过matlab代码将合成矩阵c逆向输出，获取合成二值图，在合成二值图中，黑色像素表示连通孔隙，白色代表不连通区域。

步骤s3：根据hagen-poiseuille定律和合成二值图，获取每一张合成二值图的渗透率，具体过程如下：

步骤s3.1：参考图8，根据合成二值图，获取合成二值图中连通孔隙管道的数目，以及连通孔隙管道对应的圆管半径大小，具体为：

步骤s3.1.1：对图像中的每一孔隙对象进行计算，得到其等距图，即计算每一孔隙内部的每一像素点到孔隙和基质边界的欧几里德距离，对相同的距离采用相同的颜色来表示，最后得到整体图像的等距图。也就是说，对每一张合成二值图中的每一孔隙对象进行计算，即合成二值图中的黑色像素处进行计算，得到其等距图。同时在等距图中找到欧几里德距离最大的像素点。

步骤s3.1.2：以在等距图中找到的欧几里德距离最大的像素点为圆心，并以像素点到孔隙和基质边界的距离为半径画圆，同时该圆所包围的面积即为不小于该孔径的孔隙面积，最后逐渐缩小欧几里德距离直至为0，并重复步骤s3.1.1至s3.1.2，得到累积孔隙面积随孔径的变化关系，从而确定图像中的孔隙管道数量及对应半径。

步骤s3.2：根据hagen-poiseuille定律、合成二值图中连通孔隙管道的数目和连通孔隙管道对应的圆管半径大小，获取合成二值图的渗透率。

在本实施例中，具体地讲，hagen-poiseuille定律是描述流体流经细管所产生的压力损失、压力损失和体积流率、动黏度和管长的乘积成正比，和管径的四次方成反比例。同时此定律适用于不可压缩、不具有加速度、层流稳定且长于管径的牛顿流体。其表达式为：

其中：qi为每个孔隙的流量，r为细管半径，δp为压力损失，l为细管长度，μ为动粘度。

同时流体通过试样微观结构的总流量和通过每个孔隙流量的总和是相等的，也就是说：

q＝∑qi

其中：q为试样微观结构的总流量，qi为每个孔隙的流量。

又因为达西定律：

其中：q为试样微观结构的总流量，k为渗透率，a为流体的截面面积，δp为压力损失，l为细管长度，μ为动粘度。

从而：

其中：k为渗透率，a为流体的截面面积，δp为压力损失，μ为动粘度，l为细管长度，q为试样微观结构的总流量，qi为每个孔隙的流量，ri为孔隙半径。

即合成二值图的渗透率为：

其中：a为流体的截面面积，ri为孔隙半径，ai为孔隙的面积。

步骤s4：参考图9，根据合成二值图的渗透率ki，确定岩土体渗透率，具体为：

其中：n为相邻断层图像的总数目，ki为合成二值图的渗透率。

以上结合附图对本发明的实施方式做出了详细说明，但本发明不局限于所描述的实施方式。对于本领域技术人员而言，在不脱离本发明的原理和精神的情况下，对这些实施方式进行各种变化、修改、替换和变形仍落入本发明的保护范围内。