CN106525680B - 岩心孔隙度参数场的获取方法 - Google Patents

Abstract

本发明公开了一种岩心孔隙度参数场的获取方法，属于石油地球物理勘探领域。该方法包括：在k张岩心二值化图像上切割边长为L的方形截面，并将其对应的多个二维矩阵按三维顺序叠加，形成三维矩阵并封装至单元数组的各个元胞数组中，获取每个元胞数组对应的孔隙度，进而获得单元数组的孔隙度参数场。根据单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值，测定岩心样品的平均孔隙度实测值，当平均孔隙度模拟值和平均孔隙度实测值的误差小于5％时，确定L和k的取值为切割参数，并且所获取的单元数组的孔隙度参数场作为岩心样品的岩心孔隙度参数场。该方法能保证获得的岩心孔隙度参数场真实反映岩心样品孔隙度，且减少数据量。

Description

岩心孔隙度参数场的获取方法

技术领域

本发明涉及石油地球物理勘探领域，特别涉及一种岩心孔隙度参数场的获取方法。

背景技术

数字岩心是基于二维扫描电镜图像或三维CT扫描图像，运用计算机图像处理技术，通过一定的算法完成数字岩心重构，其能够最大程度地反映出地层信息，对于岩石物理的研究具有重要的意义。通过数字岩心所建立的孔隙介质模型的孔隙结构最大程度地与实际岩石相同，这不仅可用来精确地预测岩石的物理特性，还可以用来了解不同的物理响应之间的内在关系。所以，提供一种通过构建三维数字岩心来获取岩心孔隙度的方法十分必要。

现有技术(201510526748.x)提供了一种岩心的孔隙度矩阵反演方法，该方法包括以下步骤：利用X射线对岩心样品进行CT扫描，得到岩心CT图像；对所述岩心CT图像依次进行灰度化处理和对比度增强处理，得到岩心灰度图；对所述岩心灰度图进行二值化处理，得到多个二值化矩阵，并将所述多个二值化矩阵按三维顺序叠加成三维矩阵；将所述三维矩阵封装至单元数组的各个元胞数组中，通过对每个元胞数组的孔隙度进行局部反演，随后对所有元胞数组的孔隙度进行反演，实现岩心孔隙度矩阵反演，得到岩心孔隙度参数场。其中，所述的反演即为将图像转化为孔隙度数值矩阵的过程。该方法通过获取能精确表征岩心样品孔隙结构信息的三维矩阵，并将所得到的三维矩阵构建成多个元胞数值矩阵，并对其中的孔隙度进行反演，从而获取岩心孔隙度参数场。

发明人发现现有技术至少存在以下问题：

现有技术提供的方法得到的岩心孔隙度参数场数据量大，在利用其模拟计算岩心样品的平均孔隙度时其计算量大，同时将其作为岩心数值试验研究的基础数据进行相应计算时，其计算量也很大。

发明内容

本发明实施例所要解决的技术问题在于，提供了一种既能保证获得的岩心孔隙度参数场能真实反映岩心样品孔隙度特征，又能减少其平均孔隙度模拟计算量的岩心孔隙度参数场的获取方法。具体技术方案如下：

一种岩心孔隙度参数场的获取方法，包括：利用X射线沿岩心样品的轴向间隔预定距离对所述岩心样品进行CT扫描，得到k张岩心CT图像；对所述k张岩心CT图像依次进行灰度化处理、对比度增强处理、二值化处理，得到k张岩心二值化图像。进一步地，所述方法还包括：

步骤a、在k张所述岩心二值化图像上分别切割边长为L的方形截面，并将k张所述方形截面所对应的多个二维矩阵按三维顺序叠加，形成三维矩阵；

步骤b、将所述三维矩阵封装至单元数组的各个元胞数组中，通过对每个所述元胞数组中的孔隙度参数进行局部反演，获得每个所述元胞数组对应的孔隙度，进而获得所述单元数组的孔隙度参数场；

步骤c、根据所述单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值，测定所述岩心样品的平均孔隙度实测值，当所述平均孔隙度模拟值和所述平均孔隙度实测值之间的误差小于5％时，确定所述L和所述k的取值为切割参数，并且将基于所述切割参数所获取的单元数组的孔隙度参数场作为所述岩心样品的岩心孔隙度参数场。

具体地，所述步骤c还包括：确定所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差是否大于5％；

如果否，则确定所述L和所述k的原值为所述切割参数；

如果是，则增大所述L和所述k的取值，并重复所述步骤a至所述步骤c，直至所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差小于5％，将此时确定的所述L和所述k的新值作为所述切割参数。

具体地，所述方法还包括：待所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差小于5％后，继续增大所述L和所述k的取值，并重复所述步骤a至所述步骤c，检测对应得到的所述平均孔隙度模拟值是否稳定，如果稳定，则确定满足所述稳定条件的所述L和所述k的取值的最小值为所述切割参数。

具体地，将每一张所述岩心二值化图像的边缘位置处作为所述切割边长L的定位起点；

并且，在增大所述切割边长L时，新的切割边长自所述定位起点处开始延伸并覆盖原切割边长。

具体地，通过采用饱和煤油法测定并获取所述岩心样品的平均孔隙度实测值。

具体地，通过采用加权平均法、平均值法或者最大值法对k张所述岩心CT图像进行灰度化处理，使k张所述岩心CT图像中所有像素点的R、G、B分量相同。

具体地，通过采用MATLAB软件，并调用所述MATLAB软件中图像处理工具箱中的imajust()函数，来对经所述灰度化处理后的k张所述岩心CT图像进行对比度增强处理，形成k张岩心灰度图。

具体地，通过调用所述MATLAB软件中的命令函数im2bw()将k张所述岩心灰度图设置成二值图像函数g(x,y)，以实现所述二值化处理，获得k张所述岩心二值化图像以及所述岩心二值化图像对应的多个二维矩阵；

其中，

其中，f(x,y)为岩心灰度图的对应函数，

1代表灰度值为255的岩心孔隙，

0代表灰度值为0的岩心骨架，

t代表预设的灰度阈值。

具体地，所述单元数组中包括900-10000个所述元胞数组，并且每个所述元胞数组中包括10000-90000个像素点。

具体地，所述根据所述单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值包括：

根据每个所述元胞数组对应的孔隙度，通过第一公式计算所述单元数组的平均孔隙度；

根据所述单元数组的平均孔隙度，通过第二公式计算所述三维切割样品的平均孔隙度模拟值；

所述第一公式如下所示：

其中，F为单元数组的平均孔隙度，Y为元胞数组对应的孔隙度；

所述第二公式如下所示：

其中，M为三维切割样品的平均孔隙度模拟值。

本发明实施例提供的技术方案带来的有益效果是：

本发明实施例提供的方法，通过在岩心样品的k张岩心二值化图像上分别截取同等大小，且具有预定切割边长L的方形截面，该k个方形截面按三维顺序叠加后将形成方体形三维切割样品。在此基础上，获取该三维切割样品对应的单元数组的孔隙度参数场，并进一步获取该三维切割样品的平均孔隙度，并将其与实际测量得到的该岩心样品的平均孔隙度实测值进行比较，若两者的误差小于5％时即可确定上述k和L的取值为可行的，此时就可将三维切割样品所对应的孔隙度参数场作为岩心样品的岩心孔隙度参数场，进而作为岩心数值研究的基础数据。可见，本发明实施例提供的方法不仅能较好地保留岩心样品的真实孔隙度特征，并且还能避免对岩心样品的所有岩心二值化图像的三维矩阵进行全部的反演计算，仅仅通过对三维切割样品的三维矩阵进行反演计算即可，有效提高了岩心样品的岩心孔隙度参数场的获取效率。并且，利用本发明实施例提供的方法获取得到的岩心孔隙度参数场可以作为岩心数值试验研究的基础数据使用。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例1提供的岩心样品的一张岩心二值化图像；

图2是本发明实施例1提供的从图1所示的岩心二值化图像中切割得到的方形截面；

图3是本发明实施例1提供的单元数组的局部示意图；

图4是本发明实施例1提供的单元数组的孔隙度参数场的局部示意图；

图5是本发明实施例1提供的L的取值与三维切割样品的平均孔隙度模拟值之间的关系图；

图6是本发明实施例1提供的k的取值与三维切割样品的平均孔隙度模拟值之间的关系图。

具体实施方式

除非另有定义，本发明实施例所用的所有技术术语均具有与本领域技术人员通常理解的相同的含义。为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

本发明实施例提供了一种岩心孔隙度参数场的获取方法，该方法包括：利用X射线沿岩心样品的轴向间隔预定距离对岩心样品进行CT扫描，得到k张岩心CT图像；对k张岩心CT图像依次进行灰度化处理、对比度增强处理、二值化处理，得到k张岩心二值化图像。在此基础上，该方法还包括如下步骤：

步骤101、在k张岩心二值化图像上分别切割边长为L的方形截面，并将k张方形截面所对应的多个二维矩阵按三维顺序叠加，形成三维矩阵。

步骤102、将上述三维矩阵封装至单元数组的各个元胞数组中，通过对每个元胞数组中的孔隙度参数进行局部反演，获得每个元胞数组对应的孔隙度，进而获得单元数组的孔隙度参数场。

步骤103、根据单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值，测定岩心样品的平均孔隙度实测值，当平均孔隙度模拟值和平均孔隙度实测值之间的误差小于5％时，确定L和k的取值为切割参数，并且将基于切割参数所获取的单元数组的孔隙度参数场作为岩心样品的岩心孔隙度参数场。

本发明实施例提供的方法，通过在岩心样品的k张岩心二值化图像上分别截取同等大小，且具有预定切割边长L的方形截面，该k个方形截面按三维顺序叠加后将形成方体形三维切割样品。在此基础上，获取该三维切割样品对应的单元数组的孔隙度参数场，并进一步获取该三维切割样品的平均孔隙度，并将其与实际测量得到的该岩心样品的平均孔隙度实测值进行比较，若两者的误差小于5％时即可确定上述k和L的取值为可行的，此时就可将三维切割样品所对应的孔隙度参数场作为岩心样品的岩心孔隙度参数场，进而作为岩心数值研究的基础数据。可见，本发明实施例提供的方法不仅能较好地保留岩心样品的真实孔隙度特征，并且还能避免对岩心样品的所有岩心二值化图像的三维矩阵进行全部的反演计算，仅仅通过对三维切割样品的三维矩阵进行反演计算即可，有效提高了岩心样品的岩心孔隙度参数场的获取效率。并且，利用本发明实施例提供的方法获取得到的岩心孔隙度参数场可以作为岩心数值试验研究的基础数据使用。

对于岩心样品的平均孔隙度实测值来说，通过采用饱和煤油法测定并获取该岩心样品的平均孔隙度实测值。饱和煤油法测定岩石孔隙度为本领域所常见，的，其操作一般如下所示：首先将岩心样品烘干后称其干重，然后将干燥的岩心样品放入煤油中直至其孔隙内充满煤油(通过将吸有煤油的岩心样品取出，擦干表面的煤油，并称量其湿重，待其湿重稳定后即可说明岩心样品的内部孔隙已被煤油充满饱和)，此时湿重与干重的差值即为孔隙内所充填煤油的质量，并且基于煤油的密度可知，可容易地获得孔隙的体积。

进一步地，步骤103具体还包括：确定平均孔隙度模拟值与平均孔隙度实测值之间的误差是否大于5％；

如果否，则确定L和k的原值为切割参数；

如果是，则增大L和k的取值，并重复步骤101至步骤103，直至平均孔隙度模拟值与平均孔隙度实测值之间的误差小于5％，将此时确定的L和k的新值作为切割参数。

将第一次确定的L和k的取值称为原值，如果利用L和k的原值作为切割参数在岩心样品上切割k张边长为L的方形截面，并且按三维顺序将它们叠加形成三维切割样品。第一种情况：如若基于该三维切割样品所得到的平均孔隙度模拟值与岩心样品的平均孔隙度实测值之间的误差小于5％(当然误差越小越好)时，即可说明该三维切割样品即能够准确表征岩心样品真实的孔隙度特征，表明了基于该三维切割样品所得到的单元数组的孔隙度参数场能够作为该岩心样品的岩心孔隙度参数场来使用，在保证真实度的前提下降低了岩心孔隙度参数场的支持数据量。此时，将确定L和k的原值作为切割参数来对岩心样品进行选择性保留。

第二种情况下，如若基于该三维切割样品所得到的平均孔隙度模拟值与岩心样品的平均孔隙度实测值之间的误差大于或等于5％时，此时可确定该三维切割样品即无法准确地表征岩心样品真实的孔隙度特征，即L和k的原值无法作为切割参数。在这种情况下，需要增大L和k的取值，并重复上述步骤101至步骤103进行验证，直至利用该L和k的某一新值所确定的三维切割样品的平均孔隙度模拟值与岩心样品的平均孔隙度实测值之间的误差小于5％，将此时确定的L和k的新值作为切割参数。可以理解的是，在本发明实施例中，切割参数意味着能使得到的三维切割样品的孔隙度特征真实地反应岩心样品的孔隙度特征。

在增大上述L和k两者的取值时，可以使它们同时增加，也可以先固定其中一个而使另外一个增加，优选采用后者。举例来说，可以先固定住L的取值，随后增加k的取值，待k的取值接近满足切割参数的要求时，再固定住k的取值，随后增加L的取值，直至L也满足切割参数的要求。

进一步地，本发明实施例提供的方法还包括：待平均孔隙度模拟值与平均孔隙度实测值之间的误差小于5％后，继续增大L和k的取值，并重复步骤101至步骤103，检测对应得到的平均孔隙度模拟值是否稳定，如果稳定，则确定满足稳定条件的L和k的取值的最小值为切割参数。

为了保证所得到的三维切割样品的孔隙度参数场能真实反映岩心样品的孔隙度特征，当利用特定的L和k的进行切割后，三维切割样品的平均孔隙度模拟值与平均孔隙度实测值之间的误差小于5％，此时可在此基础上进一步继续增大L和k的取值进行验证，并重复上述步骤101至步骤103，检测对应得到的三维切割样品的平均孔隙度模拟值是否稳定(即与岩心样品的平均孔隙度实测值之间的误差是否稳定在5％以内)，如果稳定，则确定满足稳定条件的L和k的取值中的最小值为切割参数，即还可以理解为将验证之前L和k的取值作为切割参数。

在对切割边长L进行定位时，优选将每一张岩心二值化图像的边缘位置处作为切割边长L的定位起点，这样能够增加L的取值范围，举例来说，当在岩心二值化图像的左上边缘处设定L的定位起点时，在增大L的取值时，其能够沿水平方向向右逐渐延伸以及沿竖直方向向下逐渐延伸，如此获得的方形截面能逐渐靠近并覆盖岩心二值化图像的中部位置，并且还能获得更多个方形截面。进一步地，在增大该切割边长L时，新的切割边长自该定位起点处开始延伸并覆盖原切割边长，通过如此设置，可使后续面积逐渐增大的新方形截面逐渐覆盖前方形截面，在前方形截面的基础上逐渐递进，避免了无规律地设定切割边长，即避免了方形截面在岩心二值化图像上的位置无规律，提高了切割参数的确定效率。可以理解的是，对于k张岩心二值化图像同时增大L的取值，并且保持L的取值相同。

通过本领域常用的X射线计算机层析成像技术即可实现X射线对岩心样品的CT扫描，其可以无损伤地检测待测岩心样品的组成和结构，且能精确地反映出岩心样品的内部结构。在扫描过程中，为了更精确地获知岩心样品的内部结构，一般沿着岩心样品的轴向间隔预定距离对岩心样品进行CT扫描，得到k张岩心CT图像，其中k不仅表示岩心CT图像的张数，张数越多，后续获取得到的单元数组的孔隙度参数场约精确，举例来说，对于一个岩心样品来说，

上述的预定距离可以为5-40个像素点的距离总长(例如5个、10个、15个、20个、25个、30个、35个等)，使岩心样品的轴向长度除以上述预定距离，然后加1即可得到k值。

在对k张岩心CT图像依次进行灰度化处理的过程中，通过采用本领域常见的加权平均法、平均值法或者最大值法对每一张岩心CT图像进行灰度化处理，以使该岩心CT图像中任意像素点的R、G、B分量相同。本领域技术人员可以理解的是，通过加权平均法、平均值法或者最大值法对岩心CT图像进行灰度化处理为本领域常用的现有技术，其可以通过使用商业软件，例如MATLAB软件来实现，本发明实施例在此不对其作具体的说明。

在对灰度化处理后的岩心CT图像进行对比度增强处理的过程中，通过采用MATLAB软件，调用MATLAB软件中图像处理工具箱中的imajust()函数，来实现上述对比度增强处理，得到优化的岩心灰度图。通过对灰度化处理后的图像进行对比度增强处理，使其从比较集中的某个灰度区间变成在全部的灰度范围内均匀分布，得到高对比度的灰度图像，以便于更好的确定灰度阈值t。

在进行二值化处理的过程中，通过调用MATLAB软件中的命令函数im2bw()将k张岩心灰度图设置成二值图像函数g(x,y)，以实现二值化处理，获得k张岩心二值化图像以及该岩心二值化图像对应的多个二维矩阵；

其中，

其中，f(x,y)为岩心灰度图的对应函数，1代表灰度值为255的岩心孔隙，0代表灰度值为0的岩心骨架，t代表预设的灰度阈值，其中，在获取岩心二值化图像的同时，针对每一张岩心二值化图像均能获得其所对应的二维矩阵。其中，灰度阈值t可调用MATLAB工具箱的graythresh函数，采用最大类间方差法得到，基于待处理的灰度图像的实际灰度，利用该graythresh函数可生成一个优选的灰度阈值t，进而获得理想的岩心二值化图像。

对于岩心二值化图像来说，其指的是仅有1值(白色)和0值(黑色)的图像，可参照图1，由于每一张岩心二值化图像均对应有一个二维矩阵，当通过步骤101在k张岩心二值化图像上分别切割边长为L的方形截面后，可以将k张方形截面所对应的多个二维矩阵按三维顺序叠加，形成三维矩阵(k张方形截面按三维顺序叠加将形成方体形状的三维切割样品)。

在得到三维切割样品的三维矩阵后，将上述三维矩阵封装至单元数组的各个元胞数组中，通过对每个元胞数组中的孔隙度参数进行局部反演，获得每个元胞数组对应的孔隙度，进而获得单元数组的孔隙度参数场。可以理解的是，上述的局部反演指的是利用元胞数组中的孔隙度参数0和1进行计算，获得该元胞数组的平均孔隙度的过程。该局部反演可通过Math Works公司提供的MATLAB软件基于如下所示的公式来实现：

其中，Y为元胞数组的平均孔隙度，b为元胞数组对应的二维矩阵中行与列的像素点的数目(像素点即指的是1值和0值)，b²为元胞数组中像素点总数目(即1值和0值的数目总和)，y为元胞数组中1值的总数目。

对于三维切割样品来说，其中所含元胞数组的数目a即为平面精度，可以理解的是，元胞数组的数目越多，其中封装的数据越少，反演的精度越高。为了提高反演精度，本发明实施例中，单元数组中的元胞数组的个数优选为900-10000。且每个元胞数组中包括10000-90000个像素点，例如为10000、20000、30000、40000、50000、60000、70000、80000或90000个像素点。每两个像素点之间的距离为1-4μm，优选为2μm。

进一步地，上述的根据单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值包括：

步骤1031、根据每个元胞数组对应的孔隙度，通过第一公式计算单元数组的平均孔隙度。

其中，第一公式如下所示：

其中，F为单元数组的平均孔隙度，Y为元胞数组对应的孔隙度。

步骤1032、根据单元数组的平均孔隙度，通过第二公式计算三维切割样品的平均孔隙度模拟值。

第二公式如下所示：

其中，M为三维切割样品的平均孔隙度模拟值。

可通过Math Works公司提供的MATLAB软件基于上述第一公式和第二公式来实现上述计算过程。在增加L和k值的前提下，可以通过调用MATLAB软件中的sum和mean函数来进行上述计算。本领域技术人员可以理解的是，在提供了各种基础数据的前提下，利用MATLAB软件进行上述计算为本领域的常规技术手段，本发明实施例在此不作具体说明。

以下将通过具体实施例1进一步地描述本发明。在以下实施例1中采用来自某油田的真实岩心样品依次经灰度化处理、对比度增强处理、二值化处理后，得到k张岩心二值化图像，基于所使用的CT扫描仪，采用MATLAB软件中的imread读取k张岩心二值化图像即可。通过饱和煤油法测得该岩心样品的平均孔隙度实测值为0.19。

实施例1

附图1示出了该岩心样品的其中一张岩心二值化图像，其共计含有567800个像素点，相邻两个像素点之间的距离为2μm。在进行岩心孔隙度参数场的获取时，首先在该岩心二值化图像的左上边缘处设定切割边长L的定位起点(以确保具有最大化的图像切割范围和可用于切割的数据)，以该岩心二值化图像左上角处第一个像素点为原点(以便于对各像素点的查找)，将定位起点设定在(300，300)处，自该定位起点开始水平向右和竖直向下延伸20个像素点，作为L的原值，随后以该切割边长L对该岩心二值化图像进行切割，得到方形截面，其如附图2所示。同样地，对该岩心样品的其他k-1张岩心二值化图像也进行同样的切割，其中相邻两张岩心二值化图像之间的间距取20个像素点的长度和，为40μm。

将k张方形截面所对应的多个二维矩阵按三维顺序叠加，形成三维矩阵。将以上各三维矩阵封装至单元数组的各个元胞数组中(其中，单元数组的局部结构如附图3所示，其包括多个元胞)，通过MATLAB软件对每个元胞数组中的孔隙度参数进行局部反演，获得每个元胞数组对应的孔隙度，进而获得单元数组的孔隙度参数场。其中，单元数组的孔隙度参数场的局部结构如附图4所示。随后，根据该单元数组的孔隙度参数场，利用MATLAB软件计算三维切割样品的平均孔隙度模拟值。其中所用的计算公式如实施方式部分所示，在此不再赘述。

将该三维切割样品的平均孔隙度模拟值和岩心样品的平均孔隙度实测值进行比较，如附图5及附图6所示，当L的取值为17时，并且k的取值为2时，所对应的三维切割样品的平均孔隙度模拟值与岩心样品的平均孔隙度实测值之间存在较大的误差，随后逐渐增大L和k的取值，重新进行上述步骤，直至L的取值等于28时，并且k的取值为20时，所对应的三维切割样品的平均孔隙度模拟值与岩心样品的平均孔隙度实测值之间的误差明显减小至0.3％左右，并且随着L和k的取值的逐渐增大，所对应的三维切割样品的平均孔隙度模拟值也十分稳定，此时，即可确定切割参数如下：L的取值等于28，并且k的取值等于20。利用该切割参数获得的三维切割样品的孔隙度参数场能真实地反应该岩心样品的岩心孔隙度参数场，并且能够作为岩心数值研究的基础数据。

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

Claims (6)

1.一种岩心孔隙度参数场的获取方法，包括：利用X射线沿岩心样品的轴向间隔预定距离对所述岩心样品进行CT扫描，得到k张岩心CT图像；对所述k张岩心CT图像依次进行灰度化处理、对比度增强处理、二值化处理，得到k张岩心二值化图像，其特征在于，所述方法还包括：

步骤a、在k张所述岩心二值化图像上分别切割边长为L的方形截面，并将k张所述方形截面所对应的多个二维矩阵按三维顺序叠加，形成三维矩阵；

步骤b、将所述三维矩阵封装至单元数组的各个元胞数组中，通过对每个所述元胞数组中的孔隙度参数进行局部反演，获得每个所述元胞数组对应的孔隙度，进而获得所述单元数组的孔隙度参数场；

步骤c、根据所述单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值，测定所述岩心样品的平均孔隙度实测值，当所述平均孔隙度模拟值和所述平均孔隙度实测值之间的误差小于5％时，确定所述L和所述k的取值为切割参数，并且将基于所述切割参数所获取的单元数组的孔隙度参数场作为所述岩心样品的岩心孔隙度参数场，

其中，所述根据所述单元数组的孔隙度参数场计算三维切割样品的平均孔隙度模拟值包括：

对每个所述元胞数组中的孔隙度参数进行局部反演，通过第一公式计算所述每个所述元胞数组对应的孔隙度；

根据每个所述元胞数组对应的孔隙度，通过第二公式计算所述单元数组的平均孔隙度；

根据所述单元数组的平均孔隙度，通过第三公式计算所述三维切割样品的平均孔隙度模拟值；

所述第一公式如下所示：

其中，Y为元胞数组的平均孔隙度，b为元胞数组对应的二维矩阵中行与列的像素点的数目，b²为元胞数组中像素点总数目，y为元胞数组中1值的总数目；

所述第二公式如下所示：

其中，F为单元数组的平均孔隙度，Y为元胞数组对应的孔隙度；

所述第三公式如下所示：

其中，M为三维切割样品的平均孔隙度模拟值；

所述步骤c还包括：确定所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差是否大于5％；

如果否，则确定所述L和所述k的原值为所述切割参数；

如果是，则增大所述L和所述k的取值，并重复所述步骤a至所述步骤c，直至所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差小于5％，将此时确定的所述L和所述k的新值作为所述切割参数；

所述方法还包括：待所述平均孔隙度模拟值与所述平均孔隙度实测值之间的误差小于5％后，继续增大所述L和所述k的取值，并重复所述步骤a至所述步骤c，检测对应得到的所述平均孔隙度模拟值是否稳定，如果稳定，则确定满足所述稳定条件的所述L和所述k的取值的最小值为所述切割参数；

将每一张所述岩心二值化图像的边缘位置处作为所述切割边长L的定位起点；

并且，在增大所述切割边长L时，新的切割边长自所述定位起点处开始延伸并覆盖原切割边长。

2.根据权利要求1所述的方法，其特征在于，通过采用饱和煤油法测定并获取所述岩心样品的平均孔隙度实测值。

3.根据权利要求1所述的方法，其特征在于，通过采用加权平均法、平均值法或者最大值法对k张所述岩心CT图像进行灰度化处理，使k张所述岩心CT图像中所有像素点的R、G、B分量相同。

4.根据权利要求3所述的方法，其特征在于，通过采用MATLAB软件，并调用所述MATLAB软件中图像处理工具箱中的imajust()函数，来对经所述灰度化处理后的k张所述岩心CT图像进行对比度增强处理，形成k张岩心灰度图。

5.根据权利要求4所述的方法，其特征在于，通过调用所述MATLAB软件中的命令函数im2bw()将k张所述岩心灰度图设置成二值图像函数g(x,y)，以实现所述二值化处理，获得k张所述岩心二值化图像以及所述岩心二值化图像对应的多个二维矩阵；

其中，

其中，f(x,y)为岩心灰度图的对应函数，

1代表灰度值为255的岩心孔隙，

0代表灰度值为0的岩心骨架，

t代表预设的灰度阈值。

6.根据权利要求1所述的方法，其特征在于，所述单元数组中包括900-10000个所述元胞数组，并且每个所述元胞数组中包括10000-90000个像素点。

Images