基于岩心样品二维孔隙图像的渗透率预测方法与流程

本发明属于石油天然气勘探开发技术，特别涉及一种通过岩心样品二维孔隙图像预测其渗透率的方法。

背景技术：

对油田开发人员而言，渗透率是一项需重点关注的储层参数。这是因为渗透率是决定一口井是否进行完井和投产的重要依据；同时，它也是钻井时油层保护、完井时射孔方案选择、最佳排液位置确定和生产速率及三次采油措施制定的重要参考标准。一直以来，渗透率主要通过物理实验测定。并且，通过物理实验测定的渗透率值与实际地层渗透率的吻合程度将直接决定油田开发人员制定的开采方案的准确性。但由于地下储层岩石结构的复杂性，导致其渗透性往往表现出非线性特征，因此常规研究很难获得其渗透率的准确解析解；而进行物理实验需要的经济成本较高，耗费的时间较长，而且对地下储层岩石某些细微孔隙的渗透特征也难以定量表征。鉴于此，探讨预测地下储层岩石渗透率的新技术、新方法就成为目前油田科研人员的当务之急。要求这些新技术、新方法既能够提高地下储层岩石渗透率预测的准确性，又要缩短预测过程所耗费的时间，并减少投资。

技术实现要素：

为了克服上述现有技术的缺陷，本发明的目的在于提供基于岩心样品二维孔隙图像的渗透率预测方法，只要孔隙和吼道与周围的岩石颗粒可以分隔开，那么二维孔隙图像就可以作为最基本的输入来预测相应岩心样品的渗透率，通过统计岩心样品二维孔隙图像中孔隙的个数和每个孔隙的面积来预测其渗透率值，本发明可提高渗透率预测的准确性，降低物理实验预测渗透率所耗费的经济成本，缩短预测时间。

为了达到上述目的，本发明的技术方案为：

基于岩心样品二维孔隙图像的渗透率预测方法，包括以下步骤：

步骤一、岩心样品二维孔隙图像的预处理：利用Photoshop软件对已有岩心样品二维孔隙图像进行裁切，获得孔隙和吼道的有效研究区域，并把非bmp格式图像全部转化为bmp格式的灰度图像，之后利用小波变换技术去除其噪声，增强其对比度；

步骤二、岩心样品二维孔隙图像的二值化：采用迭代阈值法对岩心样品二维孔隙图像进行二值化处理，完成之后整幅岩心样品二维孔隙图像由两种颜色构成：黑色和白色，黑色表示岩心样品二维孔隙图像中的孔隙与吼道，白色表示孔隙与吼道周围的岩石颗粒；

步骤三、岩心样品二维孔隙图像中孔隙个数与每个孔隙面积的统计：首先给出二值化后岩心样品二维孔隙图像中孔隙与吼道的定义，在此基础上利用图像腐蚀操作消除其中的吼道，接着通过轮廓跟踪算法对岩心样品二维孔隙图像中的孔隙轮廓进行跟踪，统计其中孔隙的个数num_pore；最后利用扫描线种子填充算法，计算出每个孔隙的面积；

步骤四：岩心样品二维孔隙图像中全部孔隙的等效半径求解：采用等效圆方法求取岩心样品二维孔隙图像中每一个孔隙的半径，即：把步骤三计算获得的岩心样品二维孔隙图像中每一个孔隙的面积等效成一个圆的面积，然后根据圆的面积计算公式S＝πr²求出每一个孔隙的等效半径r；

步骤五、岩心样品二维孔隙图像孔隙传导系数C的求取：利用公式(1)求取岩心样品二维孔隙图像的孔隙传导系数C，具体公式如下：

上式中，N为岩心样品二维孔隙图像中孔隙的总个数，它的值等于步骤三中的num_pore，i为岩心样品二维孔隙图像中孔隙的编号，c_i为岩心样品二维孔隙图像中第i个孔隙的传导系数，r_i为步骤四中求取的岩心样品二维孔隙图像中第i个孔隙的等效半径；

步骤六、基于岩心样品二维孔隙图像的渗透率预测：基于岩心样品二维孔隙图像的渗透率预测公式如下所示：

上式中，k为与岩心样品二维孔隙图像相对应的岩心样品在该截面处的渗透率值，C为该岩心样品二维孔隙图像的孔隙传导系数，A为该岩心样品二维孔隙图像的面积；

步骤七、岩心样品渗透率的计算：上述步骤只是预测了岩心样品在某一截面处的渗透率，要计算整块岩心样品的渗透率，还需对许多幅岩心样品的二维孔隙图像进行处理，求得许多个渗透率的预测值，然后取它们的最大值、最小值、平均值和标准差等，通过这些参数来表征整块岩心样品的渗透率。

所述的步骤三具体为：

(1)岩心样品二维孔隙图像中孔隙与吼道的定义：把岩心样品二维孔隙图像中由三个或三个以上岩石颗粒所包围的空间称为孔隙，把连接相邻两孔隙的最小通道称为吼道，即两颗岩石颗粒之间的空间；

(2)对岩心样品二维孔隙图像进行腐蚀操作，消除其中的吼道；

(3)基于轮廓跟踪算法岩心样品二维孔隙图像中孔隙个数统计的具体过程为：

①对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，搜索黑色像素点；

②在获取第一个黑色像素点之后，将其作为孔隙轮廓跟踪的起始点，并将该点进入方向的反方向作为起始方向链码，在确定跟踪的起始点和起始方向链码之后，完成对整个孔隙轮廓的跟踪；

③完成第一个孔隙轮廓的跟踪之后，把参数num_pore的初始值设为0，用来表示该岩心样品二维孔隙图像中孔隙的总个数，参数num_pore的值加1，并把该孔隙用颜色RGB((255-1),255,255)填充，用RGB颜色模式表示不同孔隙的颜色，RGB颜色模式可在计算机屏幕上重现16777216种不同的颜色；

④对岩心样品二维孔隙图像继续进行由左到右、由上到下的扫描，获取第二个黑色像素点，将其作为第二个孔隙跟踪的起始点，并将该点进入方向的反方向作为起始方向链码，这样确定了第二个孔隙跟踪的起始点和起始方向链码，在此基础上完成对该孔隙整个轮廓的跟踪，并把参数num_pore的值再加1，把该孔隙用颜色RGB((255-2),255,255)填充；

⑤对岩心样品二维孔隙图像继续进行由左到右、由上到下的扫描，确定第i个孔隙跟踪的起始点和起始方向链码，完成对该孔隙整个轮廓的跟踪，把参数num_pore的值再加1，并把该孔隙用颜色RGB((255-i),255,255)填充；

⑥继续对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，完成岩心样品二维孔隙图像中所有孔隙的轮廓跟踪与颜色填充和个数统计，最终实现不同的孔隙用不同的颜色表示；用num_pore的最终值表示该岩心样品二维孔隙图像中孔隙的总个数；

(4)基于扫描线种子填充算法岩心样品二维孔隙图像中每个孔隙面积计算。

所述的基于扫描线种子填充算法岩心样品二维孔隙图像中每个孔隙面积计算具体过程为：

①对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，找到第一个非白色点(x_first,y_first)，初始化一个空栈并将该点入栈；

②判断该栈是否为空，如果为空则结束该循环，否则取出栈顶元素将其作为当前扫描线的种子点(x,y)，这里x即为x_first，y即为y_first，y为当前的扫描线；从种子点(x,y)出发，沿当前扫描线向左、右两个方向填充至边界，并分别标记填充过区段的左、右边界端点为x_left和x_right，分别检查与当前扫描线相邻的y-1和y+1两条扫描线在区间中的像素，从x_left开始向x_right方向搜索，若存在非边界且尚未填充的像素点，则找出这些相邻像素点中最右边的一个，将其作为种子点压入栈中，重复上述填充过程；

③在填充过程中，用参数number₁表示填充过的所有像素的总数，同时在填充结束之后把填充区域内所有非白色像素点转化为白色像素点；

④继续对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，把扫描过的第二个填充区域内的所有像素总数记为number₂，把扫描过的第N个填充区域内的所有像素总数记为number_N；

⑤最终，由单个像素的面积乘上number₁就可以得到扫描过的第一个孔隙的面积，由单个像素的面积乘上number_N就可以得到扫描过的第N个孔隙的面积，就获得了该岩心样品二维孔隙图像中全部孔隙的面积。

本发明基于计算机图像处理技术，将岩心样品二维孔隙灰度图像的像素处理与岩心样品的渗透率预测相结合，在一定程度上可提高渗透率预测的准确性，并降低通过物理实验预测渗透率所耗费的经济成本，缩短渗透率预测所需要的时间。

附图说明

图1为本发明中基于岩心样品二维孔隙图像预测渗透率的流程图。

图2为本发明中基于轮廓跟踪算法岩心样品二维孔隙图像中统计孔隙个数的流程图。

图3为本发明中基于扫描线种子填充算法岩心样品二维孔隙图像中计算每个孔隙面积的流程图。

具体实施方式

下面结合附图对本发明的技术方案做详细叙述。

参照图1，基于岩心样品二维孔隙图像的渗透率预测方法，包括以下步骤：

步骤一、岩心样品二维孔隙图像的预处理：利用Photoshop软件对已有岩心样品二维孔隙图像进行处理，删除无效区域，获得孔隙和吼道的有效研究区域，并把非bmp格式的图像全部转化为bmp格式。之后把彩色岩心样品二维孔隙图像转化为灰度图像(如已经是灰度图像则无需转化)，最后利用小波变换技术去除其噪声，增强其对比度；

步骤二、岩心样品二维孔隙图像的二值化：利用迭代阈值法对岩心样品二维孔隙图像进行二值化处理，具体步骤为：①求取整幅岩心样品二维孔隙图像灰度的算术平均值Z_average，将其设置为初始阈值，即T⁰＝Z_average；②根据阈值T^k(k＝0,1,2,3,…)将岩心样品二维孔隙图像分割为孔隙吼道和岩石颗粒两部分，以黑色代表孔隙与吼道，白色代表岩石颗粒，再分别求出这两部分的平均灰度值Z_pore和Z_rock；③使得T^k+1＝(Z_pore+Z_rock)/2(k＝0,1,2,3,…)；④如果T^k＝T^k+1(k＝0,1,2,3,…)，则停止迭代，否则使k＝k+1，重复②-④；⑤通过上述迭代计算，求得最终阈值T_final，之后把所有灰度大于T_final的像素点的灰度值重新设置为1，把所有灰度小于T_final的像素点的灰度值重新设置为0，至此完成了岩心样品二维孔隙图像的二值化；

步骤三、岩心样品二维孔隙图像中孔隙个数与每个孔隙面积的统计：首先给出二值化后岩心样品二维孔隙图像中孔隙与吼道的定义，在此基础上利用图像腐蚀操作消除其中的吼道(本发明中所有讨论过程都是在图像的3*3领域内运算执行，后面不再重复说明)，接着利用轮廓跟踪算法对岩心样品二维孔隙图像中的孔隙轮廓进行跟踪，统计图像中孔隙的个数num_pore；最后利用扫描线种子填充算法，计算出每个孔隙的面积。具体步骤为：(1)岩心样品二维孔隙图像中孔隙与吼道的定义：把岩心样品二维孔隙图像中由三个或三个以上岩石颗粒所包围的空间称为孔隙，把连接相邻两孔隙的最小通道称为吼道，即两颗岩石颗粒之间的空间；(2)对岩心样品二维孔隙图像进行腐蚀操作，消除其中的吼道；(3)基于轮廓跟踪算法岩心样品二维孔隙图像中孔隙个数统计的基本流程见图2，具体过程为：①对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，搜索黑色像素点；②在获取第一个黑色像素点之后，将其作为孔隙轮廓跟踪的起始点，并将该点进入方向的反方向作为起始方向链码，在确定跟踪的起始点和起始方向链码之后，完成对整个孔隙轮廓的跟踪；③完成第一个孔隙轮廓的跟踪之后，把参数num_pore(其初始值设为0，用来表示该岩心样品二维孔隙图像中孔隙的总个数)的值加1，并把该孔隙用颜色RGB((255-1),255,255)填充(本专利中用RGB颜色模式表示不同孔隙的颜色，RGB颜色模式可在计算机屏幕上重现16777216种不同的颜色，在这里完全够用)；④对岩心样品二维孔隙图像继续进行由左到右、由上到下的扫描，获取第二个黑色像素点，将其作为第二个孔隙跟踪的起始点，并将该点进入方向的反方向作为起始方向链码，这样确定了第二个孔隙跟踪的起始点和起始方向链码，在此基础上完成对该孔隙整个轮廓的跟踪，并把参数num_pore的值再加1，把该孔隙用颜色RGB((255-2),255,255)填充；⑤对岩心样品二维孔隙图像继续进行由左到右、由上到下的扫描，确定第i个孔隙跟踪的起始点和起始方向链码，完成对该孔隙整个轮廓的跟踪，把参数num_pore的值再加1，并把该孔隙用颜色RGB((255-i),255,255)填充；⑥继续对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，完成岩心样品二维孔隙图像中所有孔隙的轮廓跟踪与颜色填充(最终实现不同的孔隙用不同的颜色表示)和个数统计(用num_pore的最终值表示该岩心样品二维孔隙图像中孔隙的总个数)；(4)基于扫描线种子填充算法岩心样品二维孔隙图像中每个孔隙面积计算的基本流程见图3，具体过程为：①对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，找到第一个非白色点(x_first,y_first)，初始化一个空栈并将该点入栈；②判断该栈是否为空，如果为空则结束该循环，否则取出栈顶元素将其作为当前扫描线的种子点(x,y)(这里x即为x_first，y即为y_first)，y为当前的扫描线。从种子点(x,y)出发，沿当前扫描线向左、右两个方向填充至边界，并分别标记填充过区段的左、右边界端点为x_left和x_right，分别检查与当前扫描线相邻的y-1和y+1两条扫描线在区间中的像素，从x_left开始向x_right方向搜索，若存在非边界且尚未填充的像素点，则找出这些相邻像素点中最右边的一个，将其作为种子点压入栈中，重复上述填充过程；③在填充过程中，用参数number₁表示填充过的所有像素的总数，同时在填充结束之后把填充区域内所有非白色像素点转化为白色像素点；④继续对岩心样品二维孔隙图像进行由左到右、由上到下的扫描，把扫描过的第二个填充区域内的所有像素总数记为number₂，把扫描过的第N个填充区域内的所有像素总数记为number_N；⑤最终，由单个像素的面积乘上number₁就可以得到扫描过的第一个孔隙的面积，由单个像素的面积乘上number_N就可以得到扫描过的第N个孔隙的面积，如此，就获得了该岩心样品二维孔隙图像中全部孔隙的面积；

步骤四：岩心样品二维孔隙图像中全部孔隙的等效半径求解：采用等效圆方法求取岩心样品二维孔隙图像中每一个孔隙的等效半径，即：把前面计算得到的岩心样品二维孔隙图像中每一个孔隙的面积等效成一个圆的面积，根据圆的面积计算公式S＝πr²求出每一个孔隙的等效半径r；

步骤五、岩心样品二维孔隙图像孔隙传导系数C的求取：利用公式(1)来求取岩心样品二维孔隙图像的孔隙传导系数C，具体公式如下：

上式中，N为岩心样品二维孔隙图像中孔隙的总个数，它的值等于步骤三中的num_pore，i为岩心样品二维孔隙图像中孔隙的编号，c_i为岩心样品二维孔隙图像中第i个孔隙的传导系数，r_i为步骤四中求取的岩心样品二维孔隙图像中第i个孔隙的等效半径；

步骤六、基于岩心样品二维孔隙图像的渗透率预测：Mathieu等(2007)提出了基于岩心样品二维孔隙图像的渗透率预测公式，如下所示：

上式中，k为与岩心样品二维孔隙图像相对应的岩心样品在该截面处的渗透率值，N为岩心样品二维孔隙图像中孔隙的总个数，A为该岩心样品二维孔隙图像的面积，c_eff为岩心样品二维孔隙图像中单个孔隙的传导系数。但Mathieu等并没有给出c_eff的详细求解公式，因此本发明中采用步骤五中求取的岩心样品二维孔隙图像的孔隙传导系数C来代替Nc_eff，这样得到新的渗透率预测公式(3)，如下所示：

上式中，k为与岩心样品二维孔隙图像相对应的岩心样品在该截面处的渗透率值，C为该岩心样品二维孔隙图像的孔隙传导系数，A为该岩心样品二维孔隙图像的面积；

步骤七、岩心样品渗透率的计算：上述步骤只是预测了岩心样品在某一截面处的渗透率值，要计算整块岩心样品的渗透率，还需对许多幅岩心样品的二维孔隙图像进行处理，求得许多个渗透率的预测值，然后取它们的最大值、最小值、平均值和标准差等，通过这些参数来表征整块岩心样品的渗透率。