单重煤岩孔隙结构模型的构建方法与流程

本发明涉及一种模型构建方法，尤其涉及一种单重煤岩孔隙结构模型的构建方法。

背景技术：

煤岩既是煤层气的源岩，又是煤层气的储集空间，作为储集空间，与常规岩石相比，最重要的区别是煤岩是一种双重孔隙介质，由基质孔隙和切割基质的裂隙组成，二者对煤岩的渗流和储集能力起着不同的作用。

物理、物理化学以及分子物理学的研究证实，介质的孔隙结构决定了孔隙的宏观、微观结构参数。目前煤岩研究的焦点为其微观孔隙结构特征。煤层内流体的储集、渗流能力受其微观孔隙结构影响。流体、特别是多相流体在煤岩层中的渗流特性，煤层微观孔隙结构的影响较温度、压力以及煤岩的润湿性、吸附性等的影响更重要。另外，煤岩的孔隙结构直接影响到煤岩对气体吸附及解吸特性。因此正确认识煤岩的微观孔隙结构对煤层气的勘探开发至关重要。

作者张井利用煤层微观孔隙结构研究成果解释了孔隙结构对煤层气瓦斯突出的影响，作者张玉贵在此基础上进行了危险性预测。作者郁伯铭通过研究煤岩孔隙结构参数与物性关系，为物性参数评价、分类以及预测提供经验公式。总之，煤岩的宏观性质与它的微观孔隙结构有着本质和必然的联系。

我国关于煤岩物性研究工作已开展多年，但对于煤岩微观孔隙结构的研究相对较少，有关煤岩孔隙结构的特征以观察描述法和物理测试法为主，这些方法用于对煤岩的孔隙结构进行分类或者描述其形状特征，如张慧观察大量煤岩扫描电镜图像，对煤岩中的孔隙进行了分类；郝琦利用电子扫描技术煤岩的显微孔隙按成因进行了分类；李明、李前贵、周龙刚、陈贞龙分别使用压汞法，联合电镜扫描的压汞法、扫描电镜法和液氮测法对煤岩进行了测试，对煤岩中大小孔隙的相对比例和孔隙的连通性进行了描述；傅海统计了不同观测尺度下煤岩的宏观和微观裂隙分布情况，基于分形几何理论提出了裂隙长度越短，裂隙分形维数越大的观点。但上述方法，都没有基于煤岩微观孔隙结构定量图像的分析技术，而且也没有进行孔隙结构3d建模及孔隙度计算，另外，我们已知的建模，一般是针对一个煤岩样品建一次模，建模方法复杂，成本高昂，且不具有通用性。

技术实现要素：

本发明的目的就在于提供一种解决上述问题，能对孔隙这种单重结构、孔隙加裂隙的双重结构进行模拟、的煤岩孔隙结构模型的构建方法。

为了实现上述目的，本发明采用的技术方案是这样的：一种单重煤岩孔隙结构模型的构建方法，包括以下步骤：

(1)获取煤岩样品的孔隙结构参数；

选取一大小为(a，b，c)的煤岩样品，将其分为数张薄片，对每张薄片，计算其面孔隙度、孔隙半径均值、孔隙半径标准差；

将所有薄片的面孔隙度求均值，得到煤岩样品的孔隙度p；

将所有薄片的孔隙半径均值再求均值，得到煤岩样品的孔隙半径均值b；

将所有薄片的孔隙半径标准差求均值，得到煤岩样品的孔隙半径标准差c；

(2)构建单重孔隙结构模型；

(21)预设两组圆球个数n1和n2，用蒙特卡罗方法，随机创建2个数据集x1、x2；

其中，x1包含n1个圆球，x2包含n2个圆球，n1初值0，n2＞＞n1，xi表示数据集中第i个圆球半径，x1～n(b，c²)、x2～n(b，c²)，n为正态分布；

(22)分别预设一与x1、x2对应的初始堆积空间，大小为(a，b，∑r1)、(a，b，∑r2)其中，∑r1为x1中n1个圆球半径之和，∑r2为x2中n2个圆球半径之和；

(23)用球体随机堆积算法将n1个圆球、n2个圆球放入各自对应的初始堆积空间中，采用压缩算法将初始堆积空间尺寸压缩至(a，b，c)，得到两个分别与x1、x2对应的压缩空间；

(24)分别计算两个压缩空间的的孔隙度，x1对应的压缩空间孔隙度为p1，x2对应的压缩空间孔隙度为p2；

(25)采用下式进行判断，得出实际煤岩样品中孔隙数量n；

若p1＜p，p2＞p，则调整n1＝(n1+n2)/2，n2不变，重复步骤(22)-(24)；

若p1＞p，p2＞p，则重设n1，n2不变，重复步骤(22)-(24)；

若p1＜p，p2＜p，则n1不变，重设n2，重复步骤(22)-(24)；

若p1＞p，p2＜p，则n1不变，调整n2＝(n1+n2)/2，重复步骤(22)-(24)；

若|p1-p|<0.01p，则n＝n1，将该结构作为单重孔隙结构模型，并记录此时圆球的位置信息和圆球半径；

若|p2-p|<0.01p，则n＝n2，将该结构作为单重孔隙结构模型，并记录此时圆球的位置信息和圆球半径。

作为优选：步骤(1)中，薄片的面孔隙度计算方法为：将煤岩样品的每张薄片分别放置在显微镜下，通过相机拍摄显微镜下的薄片，得到每张薄片对应的图像，用图像二值化方法处理图像，得到灰度图像，再统计黑白两种像素的数量，白色像素点的个数占所有像素点的比例为该薄片的面孔隙度。

作为优选：步骤(1)中，薄片的孔隙半径的均值计算方法为：

(11)利用边缘检测算法得到薄片中每个孔隙的边界，统计出该薄片上孔隙数量和每个孔隙面积，所述孔隙面积为每个孔隙中白色像素点的个数；

(12)假设每个孔隙为圆形，根据孔隙面积求出孔隙半径，再用该薄片上所有孔隙半径之和/孔隙数量，得出孔隙半径的均值。

(13)根据孔隙数量和孔隙半径求得每张薄片的孔隙半径的均值。

本发明中：蒙特卡罗方法是一种采用数学手段产生伪随机数的方法，能依据煤岩样品的孔隙半径均值和标准差随机生成多个不同半径的球状孔隙；本发明需要利用蒙特卡罗方法随机生成两组与圆球对应的数据集；。

球体随机堆积算法是在球体个数、球体半径和堆积空间一定的情况下，将球体随机放入堆积空间中，并保证球体与球体无重叠，本发明利用球体随机堆积方法确定初始堆积空间。

压缩算法是通过随机往某一方向移动球体使得整个堆积空间压缩至一定程度。

二分法能调整圆球或椭球的数量，使模型的孔隙度逐步逼近煤岩样品的孔隙度。

本发明主要分为三部分：

第一部分，通过煤岩样品来获取孔隙结构参数；

第二部分，利用圆球来模拟样品中的圆形的孔隙，根据第一步的参数来建立单重孔隙模型，这一步的模型仅包含圆球形的孔隙，且这种模型可以适用于仅孔隙的研究实验。

第三部分，利用椭球来模拟样品中的圆形的裂隙，改变第二步中部分圆球的纵横比让其变成椭球，建立双重孔隙模型，这一步的模型包含圆球形的孔隙和椭球形裂隙，适用于孔隙、裂隙的研究实验。

在获取煤岩样品的孔隙结构参数，我们采用了以下思路：

(1)化整为零和整合的方法，将整体样品分成多张薄片分别分析，得到每张薄片的参数，再整合在一起得到整体的参数，这种方法获得的参数，准确性高、并且具有一定的代表性，区别于前人的研究以对煤岩整体进行压汞法或者液氮法所测得参数只考虑了相互连通的孔隙相比，这种方法考虑到了独立的孔隙。

(2)采用了煤岩微观孔隙结构定量图像分析技术，例如，通过相机拍摄显微镜下的薄、二值化方法处理图像、统计黑白两种像素的数量等方法，得到薄片的面孔隙度，再通过边缘检测算法得到薄片中每个孔隙的边界等，来计算孔隙的相关参数。

在构建单重孔隙结构模型时，我们采用了以下思路：

(1)通过圆球来模拟模型中的孔隙。

(2)建模时将多种方法如蒙特卡罗方法，球体随机堆积方法、二分法、压缩算法进行了有效的结合，得到两组不同的孔隙度p1、p2。

(3)通过二分法逐步逼近煤岩样品的孔隙度，对两组圆球的数量进行调整，直至满足预设条件，此时能得到确定的圆球数量，而与该圆球数量对应的模型，就是单重孔隙结构模型。

与现有技术相比，本发明的优点在于：

(1)利用本发明模拟的模型，能对孔隙这种单重结构进行模拟，具有建模方式简单、成本低廉、且具有通用性的优点。

(2)蒙特卡洛方法可以用来模拟现实中的随机现象，因此本发明用来模拟非等径球体的半径分布情况；球体随机堆积算法具有简单易行的优点，但是其不能保证所有球体均能放置在一定的堆积空间中，与压缩算法相结合可以解决该缺点；在实际中可以运用二分法思想解决一定约束条件下的求解问题。

(3)通过统计分析可以计算得到煤岩孔隙半径的均值和标准差，以球体表示微孔隙，据此建立煤岩的单重孔隙结构模型，与实际情况吻合较好。

附图说明

图1为本发明的流程图；

图2a为实施例2中采用相机对显微镜下的的第一张薄片进行拍摄的图像；

图2b为实施例2中采用相机对显微镜下的的第二张薄片进行拍摄的图像；

图2c为实施例2中采用相机对显微镜下的的第三张薄片进行拍摄的图像；

图2d为实施例2中采用相机对显微镜下的的第四张薄片进行拍摄的图像；

图2e为实施例2中采用相机对显微镜下的的第五张薄片进行拍摄的图像；

图2f为实施例2中采用相机对显微镜下的的第六张薄片进行拍摄的图像；

图3a为图2a经二值化处理了后及边缘检测算法处理后得到的图像；

图3b为图2b经二值化处理了后及边缘检测算法处理后得到的图像；

图3c为图2c经二值化处理了后及边缘检测算法处理后得到的图像；

图3d为图2d经二值化处理了后及边缘检测算法处理后得到的图像；

图3e为图2e经二值化处理了后及边缘检测算法处理后得到的图像；

图3f为图2f经二值化处理了后及边缘检测算法处理后得到的图像；

图4为实施例2中得到的单重孔隙模型的小球半径分布直方图；

图5为实施例2中得到的单重孔隙模型的小球半径正态分布图。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例1：参见图1，一种单重煤岩孔隙结构模型的构建方法，包括以下步骤：

(1)获取煤岩样品的孔隙结构参数；

选取一大小为(a，b，c)的煤岩样品，将其分为数张薄片，对每张薄片，计算其面孔隙度、孔隙半径均值、孔隙半径标准差；

将所有薄片的面孔隙度求均值，得到煤岩样品的孔隙度p；

将所有薄片的孔隙半径均值再求均值，得到煤岩样品的孔隙半径均值b；

将所有薄片的孔隙半径标准差求均值，得到煤岩样品的孔隙半径标准差c；

薄片的面孔隙度计算方法为：将煤岩样品的每张薄片分别放置在显微镜下，通过相机拍摄显微镜下的薄片，得到每张薄片对应的图像，用图像二值化方法处理图像，得到灰度图像，再统计黑白两种像素的数量，白色像素点的个数占所有像素点的比例为该薄片的面孔隙度。

薄片的孔隙半径的均值计算方法为：

(11)利用边缘检测算法得到薄片中每个孔隙的边界，统计出该薄片上孔隙数量和每个孔隙面积，所述孔隙面积为每个孔隙中白色像素点的个数；

(12)假设每个孔隙为圆形，根据孔隙面积求出孔隙半径，再用该薄片上所有孔隙半径之和/孔隙数量，得出孔隙半径的均值。

(13)根据孔隙数量和孔隙半径求得每张薄片的孔隙半径的均值。

(2)构建单重孔隙结构模型；

(21)预设两组圆球个数n1和n2，用蒙特卡罗方法，随机创建2个数据集x1、x2；

其中，x1包含n1个圆球，x2包含n2个圆球，n1初值0，n2＞＞n1，xi表示数据集中第i个圆球半径，x1～n(b，c²)、x2～n(b，c²)，n为正态分布；

(22)分别预设一与x1、x2对应的初始堆积空间，大小为(a，b，∑r1)、(a，b，∑r2)其中，∑r1为x1中n1个圆球半径之和，∑r2为x2中n2个圆球半径之和；

(23)用球体随机堆积算法将n1个圆球、n2个圆球放入各自对应的初始堆积空间中，采用压缩算法将初始堆积空间尺寸压缩至(a，b，c)，得到两个分别与x1、x2对应的压缩空间；

(24)分别计算两个压缩空间的的孔隙度，x1对应的压缩空间孔隙度为p1，x2对应的压缩空间孔隙度为p2；

(25)采用下式进行判断，得出实际煤岩样品中孔隙数量n；

若p1＜p，p2＞p，则调整n1＝(n1+n2)/2，n2不变，重复步骤(22)-(24)；

若p1＞p，p2＞p，则重设n1，n2不变，重复步骤(22)-(24)；

若p1＜p，p2＜p，则n1不变，重设n2，重复步骤(22)-(24)；

若p1＞p，p2＜p，则n1不变，调整n2＝(n1+n2)/2，重复步骤(22)-(24)；

若|p1-p|<0.01p，则n＝n1，将该结构作为单重孔隙结构模型，并记录此时圆球的位置信息和圆球半径；

若|p2-p|<0.01p，则n＝n2，将该结构作为单重孔隙结构模型，并记录此时圆球的位置信息和圆球半径。

本发明在第一步获取样品参数时，将整体样品分成多张薄片分别分析，得到每张薄片的三个参数，再将所有薄片的三个参数分类，分别求均值，作为煤岩样品的孔隙结构参数。

其次，利用圆球来代替模型中的孔隙，得到单重孔隙模型。

构建模型时，我们采用了蒙特卡罗方法，球体随机堆积方法、二分法、压缩算法，并将这些方法和算法相互结合。

其中，蒙特卡罗方法是为了：在给定了球体半径的均值和标准差后，能得到一定数量的随机数作为球体半径。这样，我们在已知煤岩样品中b、c的值，则以x1～n(b，c²)、x2～n(b，c²)为依据，随机生成两组数据集x1、x2，我们假设x1中包含n1个圆球，那么x1内的数据为n1个圆球的半径，假设x2中包含n2个圆球，那么x2内的数据为n2个圆球的半径；且为了更好的结合二分法，我们设置n1初值0，n2＞＞n1。

球体随机堆积方法，目的是为了建立一个初始堆积空间，为了保证所有球体均能放入堆积空间内，我们将与x1对应的初始堆积空间的高度设置为∑r1，与x2对应的初始堆积空间的高度设置为∑r2，而长和宽，分别同用煤岩样品尺寸(a，b，c)中的a和b，这样，得到两个初始堆积空间，大小分别为(a，b，∑r1)、(a，b，∑r2)。

压缩算法是为了将初始堆积空间尺寸压缩至(a，b，c)，具体方法是让初始堆积空间中的各个小球沿各自随机设定的方向，当移动3个b值的距离或者与其他小球相碰撞时，停止移动；随着小球的不断移动，初始堆积空间逐渐压缩至(a,b,c)，这里的b值，为步骤(1)中得到的煤岩样品的孔隙半径均值b。

二分法在本发明中的作用是：在模型大小和孔隙度一定的情况下，确定模型中孔隙的数量，本发明方案中的步骤(25)、步骤(34)则是二分法的具体体现，在步骤(25)，通过比较p1、p2与p的关系，调整n1或n2的数量，再次计算p1、p2与p比较，再次调整n1或n2的数量，最终逐步与真实值逼近，得到满足二分法条件的n值，此时的结构，则是单重孔隙结构模型。双重孔隙结构模型的方法也同理。

在这里需要注意的地方有：

构建单重孔隙模型时，以圆球代替孔隙。我们先构建两个数据集x1、x2,由于数据集内的元素表示圆球的半径，所以实际上，两个数据集x1、x2各对应一组已知半径的圆球，也就是说，x1中n1个圆球半径为已知，x2中n2个圆球半径也为已知。那么预设初始空间时，(a，b，∑r1)、(a，b，∑r2)尺寸也为已知。

在进行压缩前，通过模拟软件将n1、n2个圆球分别放入大小为(a，b，∑r1)、(a，b，∑r2)的空间中压缩，在压缩到尺寸为(a，b，c)时，得到两个孔隙度p1、p2，比较p1、p2与p的关系，调整n1或n2的数量，接下来就是迭代的过程，调整n1或n2的数量，重新预设初始堆积空间、压缩、求孔隙度p1、p2、比较p1、p2与p的关系，调整n1或n2的数量；直到n1或n2的数量满足预设需求，此时得到最终圆球的数量，也得到了单重孔隙模型。

实施例2：参见图2到图5，一种单重煤岩孔隙结构模型的构建方法，包括以下步骤：

(1)获取煤岩样品的孔隙结构参数；

我们选取一块山西和顺地区太原组煤矿的煤岩样品，分为6张薄片，图片分别编号a，b，c，d，e，f。对每张薄片，计算其面孔隙度、孔隙半径均值、孔隙半径标准差；

具体方法为：

将上述6张薄片分别放置在显微镜下，放大40倍，使用nikon-90相机对显微镜下的薄片进行拍摄，拍摄到的6张薄片的图像如图2所示：

在图2(a)中，煤岩中细小微孔发育，较致密，呈较密集状分布，也见有部分溶蚀孔隙；图2(b)和图2(d)中，煤岩呈多孔状，细小微孔发育，较致密，呈分散状，孔隙呈密集细孔分布于局部层中，也有呈星点状分布于煤岩中；图2(c)中，煤岩发育密集微孔，局部见有少量超大溶蚀孔，呈斑块状，并发育有一条破裂缝；图2(d)和图2(f)中，煤岩中微孔较为发育，且致密，并且发育有4-5条裂隙。通过具体地观察煤岩的镜下图像，可以发现，在煤岩中绝大部分空隙空间为微孔隙，并且其形状近似于圆形，如果以球体代替微孔隙，那么就可以得到煤岩的单重孔隙结构模型。考虑到煤岩中同样地存在着极少数近似于椭球体的裂缝，如果在煤岩的单重孔隙结构模型中加入长轴与短轴差异很大的椭球体来表示裂缝，就能得到了煤岩的双重孔隙结构模型。

我们对上述六张图片进行具体分析：

例如：用图像二值化方法处理图像，得到灰度图像，再统计黑白两种像素的数量，白色像素点的个数占所有像素点的比例为该薄片的面孔隙度，我们得到六张薄片的面孔隙度。

在利用边缘检测算法得到薄片中每个孔隙的边界，得到边界图像，然后根据实施例1中的步骤(11)-(13)，得到每张薄片的孔隙半径均值；

对于每张薄片，利用该薄片上各个孔隙的半径，孔隙半径均值以及孔隙的数量，通过简单的数学标准差计算公式，求得得到每张薄片的孔隙半径标准差；

用图像二值化方法处理图像得到灰度图像、利用边缘检测算法得到薄片中每个孔隙的边界得到的边界图像，参见图3a-图3f，为了方便查看，我们对图片做了对比，图3a-图3f中，左边为二值化处理的灰度图像，右边为边缘检测算法得到的边界图像。

本实施例中，利用matlab图像分析工具将其转化为对应的灰度图像，在灰度图中，白色的像素点为孔隙部分，而黑色的像素点为骨架部分，因此通过统计白色像素点的个数计算了其所占所有像素点的比例得到薄片的面孔隙率，最后，使用matlab中的边缘检测函数得到了所有孔隙的边界，以此统计了每个薄片中孔隙的平均半径及标准差，处理结果如下表1。且由于是图像处理，得到的长度单位均为像素。

表1薄片处理结果表

从表中可以发现，煤岩的孔隙度绝大多数小于10％，而煤岩中孔隙的半径均值相差不大，为1.75像素左右。

根据表1，我们可以计算出煤岩样品孔隙度p＝9.4900％、孔隙半径均值b＝1.7371、孔隙半径标准差c＝0.6394。

(2)构建单重孔隙结构模型；

为了便于理解，我们假设煤岩样品的大小为为50×60×70，单位为像素。然后根据实施例1中的步骤(21)-(25)来构建模型。首先设置不同的圆球数量进行堆积，n1初值0，n2＞＞n1。

每次堆积时，小球半径满足正态分布，其标准差和均值为上述步骤(1)得到的煤岩样品的孔隙半径均值b、孔隙半径标准差c。分别为1.7371，0.6394。至于用于判断循环的孔隙度p也是分析煤岩薄片获得9.4900％，在计算完成后，堆积空间大小为50×60×69.58像素，得到的球体个数为781个，其半径分布情况如图4、图5所示。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。