CN116152259B - 一种基于图神经网络的储层渗透率计算方法 - Google Patents
一种基于图神经网络的储层渗透率计算方法 Download PDFInfo
- Publication number
- CN116152259B CN116152259B CN202310437515.7A CN202310437515A CN116152259B CN 116152259 B CN116152259 B CN 116152259B CN 202310437515 A CN202310437515 A CN 202310437515A CN 116152259 B CN116152259 B CN 116152259B
- Authority
- CN
- China
- Prior art keywords
- node
- connected domain
- value
- crack
- pore
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 43
- 230000035699 permeability Effects 0.000 title claims abstract description 37
- 238000013528 artificial neural network Methods 0.000 title claims abstract description 18
- 239000011148 porous material Substances 0.000 claims abstract description 53
- 238000004458 analytical method Methods 0.000 claims abstract description 38
- 239000011435 rock Substances 0.000 claims abstract description 35
- 230000009977 dual effect Effects 0.000 claims abstract description 25
- 238000000034 method Methods 0.000 claims abstract description 17
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 230000011218 segmentation Effects 0.000 claims abstract description 4
- 239000004973 liquid crystal related substance Substances 0.000 claims description 36
- 239000013598 vector Substances 0.000 claims description 27
- 238000011176 pooling Methods 0.000 claims description 23
- 238000004422 calculation algorithm Methods 0.000 claims description 21
- 239000011159 matrix material Substances 0.000 claims description 18
- 230000006870 function Effects 0.000 claims description 13
- 238000012549 training Methods 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 6
- 238000012360 testing method Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 3
- 238000003709 image segmentation Methods 0.000 claims description 3
- 238000012856 packing Methods 0.000 claims description 3
- 230000000149 penetrating effect Effects 0.000 claims description 3
- 238000004804 winding Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 2
- 238000005096 rolling process Methods 0.000 abstract 1
- 238000013135 deep learning Methods 0.000 description 3
- 238000013507 mapping Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000007637 random forest analysis Methods 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 238000012706 support-vector machine Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/40—Image enhancement or restoration by the use of histogram techniques
-
- G06T5/70—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Abstract
本发明公开了一种基于图神经网络的储层渗透率计算方法，包括以下步骤：对岩石样本进行CT扫描，得到岩石序列图像样本；对岩石序列图像样本进行预处理，并进行孔缝自动阈值分割，得到孔缝序列图像切片；将孔缝序列图像切片叠加并进行三维连通域分析，得到连通域分析后的轮廓编号，并按照物理特性区分连通域中的孔隙和裂缝；根据区分后的孔隙和裂缝以及连通域分析后的轮廓编号，构建孔隙‑裂缝双重网络模型；提取多个岩石序列图像样本的孔隙‑裂缝双重网络模型，并各自提取一个图结构，将提取的图结构在图卷积网络GCN中训练回归模型，使用训练好的模型实现渗透率的自动计算，本发明相较现有方法综合了孔隙度分析，扩展性更强。
Description
技术领域
本发明涉及渗透率的计算方法，尤其是一种基于图神经网络的储层渗透率计算方法。
背景技术
近年来，深度学习由于其能够自动提取数据集中丰富的特征关系，在地质和岩石领域也有了广泛的应用。渗透率是反映岩石储流能力的重要数字表征，渗透率高精度的自动计算对于石油，煤层气，页岩气的运移以及勘探开发有重要作用，深度学习计算渗透率的方式是根据储层孔隙度，孔隙结构和孔隙连通性等参数学习一种映射关系，利用学习的映射函数高效计算结果。
目前自动计算渗透率的主要方式有两种，一类是基于规则编程的常规算法，这类算法需要将孔隙和裂缝节点化，构建链接节点的喉道，搭建出能反映孔隙-裂缝介质真实结构的数字岩心模型，之后利用泊肃叶定律，达西定律等规则计算渗透率，这类算法局限性在于定律需要遵循一定假设条件，且单一规则不具备很好的扩展性，只适用于特定样本；另一类算法是基于深度学习的自动预测算法，利用实验获得的岩心分析数据以及测井数据构建多维特征向量，使用支持向量机，随机森林等机器学习算法建立数学模型来预测渗透率，这类算法的缺点在于数据的获取难度较大，算法不能很好的结合数字岩心中孔隙结构信息。
发明内容
针对现有技术中的上述不足，本发明提供的一种基于图神经网络的储层渗透率计算方法解决了现有技术未考虑孔隙度对渗透率影响较大的问题。
为了达到上述发明目的，本发明采用的技术方案为一种基于图神经网络的储层渗透率计算方法，包括以下步骤：
S1、对岩石样本进行CT扫描，得到岩石序列图像样本；
S2、对岩石序列图像样本进行预处理，并进行孔缝自动阈值分割，得到孔缝序列图像切片；
S3、将孔缝序列图像切片叠加并进行三维连通域分析，得到连通域分析后的轮廓编号，并按照物理特性区分连通域中的孔隙和裂缝；
S4、根据区分后的孔隙和裂缝以及连通域分析后的轮廓编号，构建孔隙网络模型和裂缝网络模型，并将孔隙网络模型和裂缝网络模型融合构建孔隙-裂缝双重网络模型；
S5、计算孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度；
S6、提取多个岩石序列图像样本的孔隙-裂缝双重网络模型，并利用孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度训练回归模型，得到训练好的模型，使用训练好的模型实现渗透率的自动计算。
进一步地：所述步骤S2包括以下分步骤：
S21、对岩石序列图像样本进行高斯滤波去噪，得到图像灰度直方图分布Hist；
其中，高斯滤波的卷积核大小为5*5；
S22、对图像灰度直方图分布Hist使用0.25[1，2，1]模板平滑法处理得到连续可导函数
；
S23、对函数
计算一阶导数得到/>
，以/>
中第一个导数为零的点为起点，在查找范围T∈[0，100]内寻找导数值变化小于等于7个阈值的连续点对应的灰度值集合，以灰度值集合最右边的值作为最佳灰度阈值th，对图像分割，将低于最佳灰度阈值th的部分作为孔缝序列图像切片。
进一步地：所述步骤S3包括以下分步骤：
S31、将孔缝序列图像切片叠加得到三维数据体D；
其中，D[x，y，z]=0为背景像素，D[x，y，z]=255为孔隙像素；
S32、设置编号
和并查集字典U，遍历三维数据体D，得到修改后的三维数据体D’；
S33、遍历修改后的三维数据体D’，当前节点为非零像素点时，查询并查集字典U，将根节点的值赋予当前节点，完成遍历则完成三维连通域分析，并得到连通域分析后的轮廓编号；
当当前节点为非零像素点时，以当前节点的像素值为查找键，查询并查集字典U中对应的值；
若查找键和对应的值不相等，则将查找到的对应的值作为查找键继续迭代查询，直至键和值相等，并将相等的值赋值给当前节点，继续遍历下一个节点，完成遍历则完成三维连通域分析，得到连通域分析的结果；
其中，每个轮廓对应一个轮廓编号，轮廓内像素点的像素值为轮廓编号；
S34、根据连通域分析的结果，获取连通域坐标（x min，x max），（y min，y max），（z min，z max），计算连通域外接长方体的最长边与最短边比值ratio，根据ratio判断连通域孔缝类型；
ratio的计算公式为：
其中，若ratio大于等于5，将该连通域孔缝类型标记为裂缝，裂缝区域中的像素点为裂缝像素点；
若当ratio小于5，则将该连通域孔缝类型标记为孔隙，孔隙区域中的像素点为孔隙像素点；
x min代表连通域在x轴最小值，x max代表连通域在x轴最大值，ymin代表连通域在y轴最小值，ymax代表连通域在y轴最大值，zmin代表连通域在z轴最小值，zmax代表连通域在z轴最大值。
进一步地：所述步骤S32中，得到修改后的三维数据体D’的具体步骤为：
当当前节点为孔隙像素点时，查找孔隙像素点空间上的十三邻域孔隙像素集合S；
若S不为空，则将S中最小值min赋值给当前节点，遍历集合S，将每次查找到的像素值作为键，min作为值组合成一个键值对，到并查集字典U中查找是否存在该键值对，若不存在则将该键值对加入并查集字典U ；
若S为空，则表示十三邻域中没有孔隙像素，将当前节点赋值为
的值；之后，/>
的值加1，并在并查集字典U中加入一个键值对；
其中，该键值对的键和值均为当前节点的值，
的初始值为1。
进一步地：所述步骤S4包括以下分步骤：
S41、对连通域分析后的三维数据体D’ ’中所有的孔隙像素点采用扩张算法找到以对应点为中心的内切球，并根据连通域分析后的轮廓编号建立链表；
所述链表的节点表示为（轮廓编号，x坐标，y坐标，z坐标，内切球半径）；
S42、对链表进行简化，简化方式为：若链表中一个节点的内切球是另一个节点内切球的子集，则从链表中删除该节点，得到简化后的链表；
S43、对每个简化后的链表采用成簇算法将节点分为不同的簇，不同簇间的相连节点为喉道节点，内切球为主球节点，保留喉道节点和主球节点，将其余节点从链表删除，根据喉道节点与主球节点的链路构建孔隙网络模型；
S44、使用基于距离变化的细化算法，计算裂缝像素点到岩石基质的最短距离并将其作为像素值，查找裂缝像素点二十六邻域内对立方向像素值比中心点值小的数量对A，若A大于7，则将该裂缝像素点作为中心点，并提取出该裂缝的中心线；
S45、对中心线上的像素点使用最大球填充法找到最大内切球计算半径，并采用链表的方式存储，基于链表中节点间的链路构建裂缝网络模型；
S46、从连通域分析后的三维数据体D’ ’中找到有接触面的孔隙区域和裂缝区域，即找到了孔隙网络模型和裂缝网络模型中对应的链表，找到两个链表对应的最大球节点，并将两个最大球节点链接起来融合成一个新的链表，构建出孔隙-裂缝双重网络模型。
进一步地：所述步骤S5包括以下分步骤：
S51、统计孔缝序列图像切片中的孔缝像素数和背景像素数，并计算孔缝序列图像切片的孔隙度，孔缝序列图像切片的孔隙度
计算公式如下：
其中，P pore（.）为孔缝像素数，P skeleton（.）为像素数，i为切片计数标识，
为孔缝序列图像切片的孔隙度，n为图像切片数量；
S52、对连通域分析后的三维数据体D’ ’中所有孔缝轮廓提取中心线链表，并根据提取的中心线链表计算迂曲度
，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与提取的中心线的孔缝轮廓的轮廓编号相同的链表节点的迂曲度/>
，迂曲度/>
的计算公式如下：
其中，
为迂曲度，m为中心线链表中节点个数，j为节点计数标识，/>
为第j个节点的位置坐标，L m为第m个节点的位置坐标，/>
为第j个节点与第j+1个节点间的欧氏距离，/>
为第1个节点与第m个节点间的欧氏距离；
S53、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，计算各个轮廓编号对应的粗糙度，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度。
进一步地：所述步骤S53包括以下分步骤：
S53-1、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，统计z=h平面上和z=h平面下的像素点个数，以z=h平面上的像素点个数和z=h平面下的像素点个数数量值最接近的平面作为平均面；
其中，h∈[0，z max]为z轴方向的一个平面；
S53-2、遍历孔隙区域和裂缝区域中所有孔缝像素点，统计每个像素点周围二十六个邻域的数据；
当邻域内不全是孔隙像素，则判定该孔隙像素点为轮廓面像素点；
S53-3、计算轮廓面像素点到平均面的距离之和，得到孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度S a，其公式为：
其中，S a为粗糙度，l x为平均面的宽度，l y为平均面的长度，E(x，y)为被测量的轮廓面和建立的平均面之间的z坐标距离的算术平均。
进一步地：所述步骤S6包括以下分步骤：
S61、提取孔隙-裂缝双重网络模型中贯穿整个网络模型的图结构，并用邻接矩阵存储图信息，试验渗透率作为图结构的标签信息；
其中，图结构中主球节点的特征向量表示为(最大球半径，粗糙度S a，迂曲度
)；
S62、将从多个岩石样本提取到的图结构打包成小批量，生成一张含r个不相连小图的大图，其中，大图的邻接矩阵对应为在对角线上把r张小图的邻接矩阵拼接在一起；
S63、使用图神经网络架构选取图卷积网络GCN对图结构进行特征更新，获得每个小图中的各个节点；
S64、对更新后的每个小图中的各个节点的特征向量的集合依次进行平均池化、最大池化和最小池化的处理，并将处理后的节点向量拼接成3*C的特征矩阵；
其中C为图结构更新后的节点维度；
S65、将节点更新后的每个小图使用社区检测的图分解算法，将每个小图分解为四个子图和八个子图，分别对每个子图进行平均池化、最大池化和最小池化的处理得到12*C和24*C的特征矩阵，并将所有特征矩阵拼接为39*C的大特征矩阵；
S66、计算大特征矩阵中39个信息向量的注意力分布，并计算其加权平均后与小图对应的孔隙度特征
拼接得到图水平向量表示G w，其计算公式如下：
其中，G w为第w个小图的向量特征表示，
为第w个小图的孔隙度，G为count个小图的向量特征表示集合，count为岩石样本的总个数，X i为第i个输入向量，，/>
为第i个输入向量X i的权重，concat(，)为将第二个参数拼接到第一个参数尾部的拼接函数，q为查询值，exp(.)为取自然常数e的x次幂运算，dot(.)为点积运算函数；
S67、将G作为全连接层的输入，选择交叉熵损失函数计算模型的损失，并使用向后传播的方式更新参数，将多个岩石序列图像样本中的80%作为训练样本，其余作为测试样本，并将迭代次数设置为50，采用Adam梯度下降学习的速率调整方式，并将学习率设置为0.001，在图卷积网络GCN中训练回归模型，得到训练好的模型，用于储层渗透率的自动计算。
本发明的有益效果为：
1.本发明设计的岩石储层渗透率自动计算流程，相较于常规算法不仅考虑到了岩石孔隙结构信息，还综合了孔隙度对渗透率影响较大的特征进行分析，能处理多种地质样本，扩展性更强；
2.利用空间金字塔池化的思想，用平均池化、最大池化和最小池化的处理方式取代单一的全局平均池化方式，减少信息的损失，加入注意力机制对不同信息向量加权，能够提升模型效果；
3用神经网络学习到了孔缝网络结构和渗透率，孔隙度和渗透率的映射关系，可提高渗透率的计算精度，计算的效率也得到了极大的提升。
附图说明
图1为本发明所述的基于图神经网络的储层渗透率计算方法的流程图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所示，在本发明的一个实施例中，提供了一种基于图神经网络的储层渗透率计算方法，包括以下步骤：
S1、对岩石样本进行CT扫描，得到岩石序列图像样本；
S2、对岩石序列图像样本进行预处理，并进行孔缝自动阈值分割，得到孔缝序列图像切片；
S3、将孔缝序列图像切片叠加并进行三维连通域分析，得到连通域分析后的轮廓编号，并按照物理特性区分连通域中的孔隙和裂缝；
S4、根据区分后的孔隙和裂缝以及连通域分析后的轮廓编号，构建孔隙网络模型和裂缝网络模型，并将孔隙网络模型和裂缝网络模型融合构建孔隙-裂缝双重网络模型；
S5、计算孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度；
S6、提取多个岩石序列图像样本的孔隙-裂缝双重网络模型，并利用孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度训练回归模型，得到训练好的模型，使用训练好的模型实现渗透率的自动计算。
在本实施例中，所述步骤S2包括以下分步骤：
S21、对岩石序列图像样本进行高斯滤波去噪，得到图像灰度直方图分布Hist；
其中，高斯滤波的卷积核大小为5*5；
S22、对图像灰度直方图分布Hist使用0.25[1，2，1]模板平滑法处理得到连续可导函数
；
S23、对函数
计算一阶导数得到/>
，以/>
中第一个导数为零的点为起点，在查找范围T∈[0，100]内寻找导数值变化小于等于7个阈值的连续点对应的灰度值集合，以灰度值集合最右边的值作为最佳灰度阈值th，对图像分割，将低于最佳灰度阈值th的部分作为孔缝序列图像切片。
在本实施例中，所述步骤S3包括以下分步骤：
S31、将孔缝序列图像切片叠加得到三维数据体D；
其中，D[x，y，z]=0为背景像素，D[x，y，z]=255为孔隙像素；
S32、设置编号
和并查集字典U，遍历三维数据体D，得到修改后的三维数据体D’；
S33、遍历修改后的三维数据体D’，当前节点为非零像素点时，查询并查集字典U，将根节点的值赋予当前节点，完成遍历则完成三维连通域分析，并得到连通域分析后的轮廓编号；
当当前节点为非零像素点时，以当前节点的像素值为查找键，查询并查集字典U中对应的值；
若查找键和对应的值不相等，则将查找到的对应的值作为查找键继续迭代查询，直至键和值相等，并将相等的值赋值给当前节点，继续遍历下一个节点，完成遍历则完成三维连通域分析，得到连通域分析的结果；
其中，每个轮廓对应一个轮廓编号，轮廓内像素点的像素值为轮廓编号；
S34、根据连通域分析的结果，获取连通域坐标（x min，x max），（y min，y max），（z min，z max），计算连通域外接长方体的最长边与最短边比值ratio，根据ratio判断连通域孔缝类型；
ratio的计算公式为：
其中，若ratio大于等于5，将该连通域孔缝类型标记为裂缝，裂缝区域中的像素点为裂缝像素点；
若当ratio小于5，则将该连通域孔缝类型标记为孔隙，孔隙区域中的像素点为孔隙像素点；
x min代表连通域在x轴最小值，x max代表连通域在x轴最大值，ymin代表连通域在y轴最小值，ymax代表连通域在y轴最大值，zmin代表连通域在z轴最小值，zmax代表连通域在z轴最大值。
在本实施例中，所述步骤S32中，得到修改后的三维数据体D’的具体步骤为：
当当前节点为孔隙像素点时，查找孔隙像素点空间上的十三邻域孔隙像素集合S；
若S不为空，则将S中最小值min赋值给当前节点，遍历集合S，将每次查找到的像素值作为键，min作为值组合成一个键值对，到并查集字典U中查找是否存在该键值对，若不存在则将该键值对加入并查集字典U ；
若S为空，则表示十三邻域中没有孔隙像素，将当前节点赋值为
的值；之后，/>
的值加1，并在并查集字典U中加入一个键值对；
其中，该键值对的键和值均为当前节点的值，
的初始值为1。
在本实施例中，所述步骤S4包括以下分步骤：
S41、对连通域分析后的三维数据体D’ ’中所有的孔隙像素点采用扩张算法找到以对应点为中心的内切球，并根据连通域分析后的轮廓编号建立链表；
所述链表的节点表示为（轮廓编号，x坐标，y坐标，z坐标，内切球半径）；
S42、对链表进行简化，简化方式为：若链表中一个节点的内切球是另一个节点内切球的子集，则从链表中删除该节点，得到简化后的链表；
S43、对每个简化后的链表采用成簇算法将节点分为不同的簇，不同簇间的相连节点为喉道节点，内切球为主球节点，保留喉道节点和主球节点，将其余节点从链表删除，根据喉道节点与主球节点的链路构建孔隙网络模型；
S44、使用基于距离变化的细化算法，计算裂缝像素点到岩石基质的最短距离并将其作为像素值，查找裂缝像素点二十六邻域内对立方向像素值比中心点值小的数量对A，若A大于7，则将该裂缝像素点作为中心点，并提取出该裂缝的中心线；
S45、对中心线上的像素点使用最大球填充法找到最大内切球计算半径，并采用链表的方式存储，基于链表中节点间的链路构建裂缝网络模型；
S46、从连通域分析后的三维数据体D’ ’中找到有接触面的孔隙区域和裂缝区域，即找到了孔隙网络模型和裂缝网络模型中对应的链表，找到两个链表对应的最大球节点，并将两个最大球节点链接起来融合成一个新的链表，构建出孔隙-裂缝双重网络模型。
在本实施例中，所述步骤S5包括以下分步骤：
S51、统计孔缝序列图像切片中的孔缝像素数和背景像素数，并计算孔缝序列图像切片的孔隙度，孔缝序列图像切片的孔隙度
计算公式如下：
其中，P pore（.）为孔缝像素数，P skeleton（.）为像素数，i为切片计数标识，
为孔缝序列图像切片的孔隙度，n为图像切片数量；
S52、对连通域分析后的三维数据体D’ ’中所有孔缝轮廓提取中心线链表，并根据提取的中心线链表计算迂曲度
，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与提取的中心线的孔缝轮廓的轮廓编号相同的链表节点的迂曲度/>
，迂曲度/>
的计算公式如下：
其中，
为迂曲度，m为中心线链表中节点个数，j为节点计数标识，/>
为第j个节点的位置坐标，L m为第m个节点的位置坐标，/>
为第j个节点与第j+1个节点间的欧氏距离，/>
为第1个节点与第m个节点间的欧氏距离；
S53、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，计算各个轮廓编号对应的粗糙度，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度。
所述步骤S53包括以下分步骤：
S53-1、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，统计z=h平面上和z=h平面下的像素点个数，以z=h平面上的像素点个数和z=h平面下的像素点个数数量值最接近的平面作为平均面；
其中，h∈[0，z max]为z轴方向的一个平面；
S53-2、遍历孔隙区域和裂缝区域中所有孔缝像素点，统计每个像素点周围二十六个邻域的数据；
当邻域内不全是孔隙像素，则判定该孔隙像素点为轮廓面像素点；
S53-3、计算轮廓面像素点到平均面的距离之和，得到孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度S a，其公式为：
其中，S a为粗糙度，l x为平均面的宽度，l y为平均面的长度，E(x，y)为被测量的轮廓面和建立的平均面之间的z坐标距离的算术平均。
在本实施例中，所述步骤S6包括以下分步骤：
S61、提取孔隙-裂缝双重网络模型中贯穿整个网络模型的图结构，并用邻接矩阵存储图信息，试验渗透率作为图结构的标签信息；
其中，图结构中主球节点的特征向量表示为(最大球半径，粗糙度S a，迂曲度
)；
S62、将从多个岩石样本提取到的图结构打包成小批量，生成一张含r个不相连小图的大图，其中，大图的邻接矩阵对应为在对角线上把r张小图的邻接矩阵拼接在一起；
S63、使用图神经网络架构选取图卷积网络GCN对图结构进行特征更新，获得每个小图中的各个节点；
S64、对更新后的每个小图中的各个节点的特征向量的集合依次进行平均池化、最大池化和最小池化的处理，并将处理后的节点向量拼接成3*C的特征矩阵；
其中C为图结构更新后的节点维度；
S65、将节点更新后的每个小图使用社区检测的图分解算法，将每个小图分解为四个子图和八个子图，分别对每个子图进行平均池化、最大池化和最小池化的处理得到12*C和24*C的特征矩阵，并将所有特征矩阵拼接为39*C的大特征矩阵；
S66、计算大特征矩阵中39个信息向量的注意力分布，并计算其加权平均后与小图对应的孔隙度特征
拼接得到图水平向量表示G w，其计算公式如下：
其中，G w为第w个小图的向量特征表示，
为第w个小图的孔隙度，G为count个小图的向量特征表示集合，count为岩石样本的总个数，X i为第i个输入向量，/>
为第i个输入向量X i的权重，concat(，)为将第二个参数拼接到第一个参数尾部的拼接函数，q为查询值，exp(.)为取自然常数e的x次幂运算，dot(.)为点积运算函数；
S67、将G作为全连接层的输入，选择交叉熵损失函数计算模型的损失，并使用向后传播的方式更新参数，将多个岩石序列图像样本中的80%作为训练样本，其余作为测试样本，并将迭代次数设置为50，采用Adam梯度下降学习的速率调整方式，并将学习率设置为0.001，在图卷积网络GCN中训练回归模型，得到训练好的模型，用于储层渗透率的自动计算。
在本发明的描述中，需要理解的是，术语“中心”、“厚度”、“上”、“下”、“水平”、“顶”、“底”、“内”、“外”、“径向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的设备或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性或隐含指明的技术特征的数量。因此，限定由“第一”、“第二”、“第三”的特征可以明示或隐含地包括一个或者更多个该特征。
Claims (6)
1.一种基于图神经网络的储层渗透率计算方法，其特征在于，包括以下步骤：
S1、对岩石样本进行CT扫描，得到岩石序列图像样本；
S2、对岩石序列图像样本进行预处理，并进行孔缝自动阈值分割，得到孔缝序列图像切片；
S3、将孔缝序列图像切片叠加并进行三维连通域分析，得到连通域分析后的轮廓编号，并按照物理特性区分连通域中的孔隙和裂缝；
S4、根据区分后的孔隙和裂缝以及连通域分析后的轮廓编号，构建孔隙网络模型和裂缝网络模型，并将孔隙网络模型和裂缝网络模型融合构建孔隙-裂缝双重网络模型；
S5、计算孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度；
所述步骤S5包括以下分步骤：
S51、统计孔缝序列图像切片中的孔缝像素数和背景像素数，并计算孔缝序列图像切片的孔隙度，孔缝序列图像切片的孔隙度
计算公式如下：
其中，P pore（.）为孔缝像素数，P skeleton（.）为像素数，i为切片计数标识，
为孔缝序列图像切片的孔隙度，n为图像切片数量；
S52、对连通域分析后的三维数据体D’ ’中所有孔缝轮廓提取中心线链表，并根据提取的中心线链表计算迂曲度
，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与提取的中心线的孔缝轮廓的轮廓编号相同的链表节点的迂曲度/>
，迂曲度/>
的计算公式如下：
其中，
为迂曲度，m为中心线链表中节点个数，j为节点计数标识，/>
为第j个节点的位置坐标，L m为第m个节点的位置坐标，/>
为第j个节点与第j+1个节点间的欧氏距离，/>
为第1个节点与第m个节点间的欧氏距离；
S53、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，计算各个轮廓编号对应的粗糙度，即为孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度；
S6、提取多个岩石序列图像样本的孔隙-裂缝双重网络模型，并利用孔缝序列图像切片的孔隙度以及孔隙-裂缝双重网络模型中每个节点的迂曲度和粗糙度训练回归模型，得到训练好的模型，使用训练好的模型实现渗透率的自动计算；
所述步骤S6包括以下分步骤：
S61、提取孔隙-裂缝双重网络模型中贯穿整个网络模型的图结构，并用邻接矩阵存储图信息，试验渗透率作为图结构的标签信息；
其中，图结构中主球节点的特征向量表示为(最大球半径，粗糙度S a，迂曲度
)；
S62、将从多个岩石样本提取到的图结构打包成小批量，生成一张含r个不相连小图的大图，其中，大图的邻接矩阵对应为在对角线上把r张小图的邻接矩阵拼接在一起；
S63、使用图神经网络架构选取图卷积网络GCN对图结构进行特征更新，获得每个小图中的各个节点；
S64、对更新后的每个小图中的各个节点的特征向量的集合依次进行平均池化、最大池化和最小池化的处理，并将处理后的节点向量拼接成3*C的特征矩阵；
其中C为图结构更新后的节点维度；
S65、将节点更新后的每个小图使用社区检测的图分解算法，将每个小图分解为四个子图和八个子图，分别对每个子图进行平均池化、最大池化和最小池化的处理得到12*C和24*C的特征矩阵，并将所有特征矩阵拼接为39*C的大特征矩阵；
S66、计算大特征矩阵中39个信息向量的注意力分布，并计算其加权平均后与小图对应的孔隙度特征
拼接得到图水平向量表示G w，其计算公式如下：
其中，G w为第w个小图的向量特征表示，
为第w个小图的孔隙度，G为count个小图的向量特征表示集合，count为岩石样本的总个数，X i为第i个输入向量，/>
为第i个输入向量X i的权重，concat(，)为将第二个参数拼接到第一个参数尾部的拼接函数，q为查询值，exp(.)为取自然常数e的x次幂运算，dot(.)为点积运算函数；
S67、将G作为全连接层的输入，选择交叉熵损失函数计算模型的损失，并使用向后传播的方式更新参数，将多个岩石序列图像样本中的80%作为训练样本，其余作为测试样本，并将迭代次数设置为50，采用Adam梯度下降学习的速率调整方式，并将学习率设置为0.001，在图卷积网络GCN中训练回归模型，得到训练好的模型，用于储层渗透率的自动计算。
2.根据权利要求1所述的基于图神经网络的储层渗透率计算方法，其特征在于，所述步骤S2包括以下分步骤：
S21、对岩石序列图像样本进行高斯滤波去噪，得到图像灰度直方图分布Hist；
其中，高斯滤波的卷积核大小为5*5；
S22、对图像灰度直方图分布Hist使用0.25[1，2，1]模板平滑法处理得到连续可导函数
；
S23、对函数
计算一阶导数得到/>
，以/>
中第一个导数为零的点为起点，在查找范围T∈[0，100]内寻找导数值变化小于等于7个阈值的连续点对应的灰度值集合，以灰度值集合最右边的值作为最佳灰度阈值th，对图像分割，将低于最佳灰度阈值th的部分作为孔缝序列图像切片。
3.根据权利要求2所述的基于图神经网络的储层渗透率计算方法，其特征在于，所述步骤S3包括以下分步骤：
S31、将孔缝序列图像切片叠加得到三维数据体D；
其中，D[x，y，z]=0为背景像素，D[x，y，z]=255为孔隙像素；
S32、设置编号
和并查集字典U，遍历三维数据体D，得到修改后的三维数据体D’；
S33、遍历修改后的三维数据体D’，当前节点为非零像素点时，查询并查集字典U，将根节点的值赋予当前节点，完成遍历则完成三维连通域分析，并得到连通域分析后的轮廓编号；
当当前节点为非零像素点时，以当前节点的像素值为查找键，查询并查集字典U中对应的值；
若查找键和对应的值不相等，则将查找到的对应的值作为查找键继续迭代查询，直至键和值相等，并将相等的值赋值给当前节点，继续遍历下一个节点，完成遍历则完成三维连通域分析，得到连通域分析的结果；
其中，每个轮廓对应一个轮廓编号，轮廓内像素点的像素值为轮廓编号；
S34、根据连通域分析的结果，获取连通域坐标（x min，x max），（y min，y max），（z min，z max），计算连通域外接长方体的最长边与最短边比值ratio，根据ratio判断连通域孔缝类型；
ratio的计算公式为：
其中，若ratio大于等于5，将该连通域孔缝类型标记为裂缝，裂缝区域中的像素点为裂缝像素点；
若当ratio小于5，则将该连通域孔缝类型标记为孔隙，孔隙区域中的像素点为孔隙像素点；
x min代表连通域在x轴最小值，x max代表连通域在x轴最大值，ymin代表连通域在y轴最小值，ymax代表连通域在y轴最大值，zmin代表连通域在z轴最小值，zmax代表连通域在z轴最大值。
4.根据权利要求3所述的基于图神经网络的储层渗透率计算方法，其特征在于：所述步骤S32中，得到修改后的三维数据体D’的具体步骤为：
当当前节点为孔隙像素点时，查找孔隙像素点空间上的十三邻域孔隙像素集合S；
若S不为空，则将S中最小值min赋值给当前节点，遍历集合S，将每次查找到的像素值作为键，min作为值组合成一个键值对，到并查集字典U中查找是否存在该键值对，若不存在则将该键值对加入并查集字典U ；
若S为空，则表示十三邻域中没有孔隙像素，将当前节点赋值为
的值；之后，/>
的值加1，并在并查集字典U中加入一个键值对；
其中，该键值对的键和值均为当前节点的值，
的初始值为1。
5.根据权利要求4所述的基于图神经网络的储层渗透率计算方法，其特征在于，所述步骤S4包括以下分步骤：
S41、对连通域分析后的三维数据体D’ ’中所有的孔隙像素点采用扩张算法找到以对应点为中心的内切球，并根据连通域分析后的轮廓编号建立链表；
所述链表的节点表示为（轮廓编号，x坐标，y坐标，z坐标，内切球半径）；
S42、对链表进行简化，简化方式为：若链表中一个节点的内切球是另一个节点内切球的子集，则从链表中删除该节点，得到简化后的链表；
S43、对每个简化后的链表采用成簇算法将节点分为不同的簇，不同簇间的相连节点为喉道节点，内切球为主球节点，保留喉道节点和主球节点，将其余节点从链表删除，根据喉道节点与主球节点的链路构建孔隙网络模型；
S44、使用基于距离变化的细化算法，计算裂缝像素点到岩石基质的最短距离并将其作为像素值，查找裂缝像素点二十六邻域内对立方向像素值比中心点值小的数量对A，若A大于7，则将该裂缝像素点作为中心点，并提取出该裂缝的中心线；
S45、对中心线上的像素点使用最大球填充法找到最大内切球计算半径，并采用链表的方式存储，基于链表中节点间的链路构建裂缝网络模型；
S46、从连通域分析后的三维数据体D’ ’中找到有接触面的孔隙区域和裂缝区域，即找到了孔隙网络模型和裂缝网络模型中对应的链表，找到两个链表对应的最大球节点，并将两个最大球节点链接起来融合成一个新的链表，构建出孔隙-裂缝双重网络模型。
6.根据权利要求5所述的基于图神经网络的储层渗透率计算方法，其特征在于，所述步骤S53包括以下分步骤：
S53-1、获取连通域分析后的三维数据体D’ ’中每个轮廓编号对应的轮廓的指定裂缝的外接长方体，统计z=h平面上和z=h平面下的像素点个数，以z=h平面上的像素点个数和z=h平面下的像素点个数数量值最接近的平面作为平均面；
其中，h∈[0，z max]为z轴方向的一个平面；
S53-2、遍历孔隙区域和裂缝区域中所有孔缝像素点，统计每个像素点周围二十六个邻域的数据；
当邻域内不全是孔隙像素，则判定该孔隙像素点为轮廓面像素点；
S53-3、计算轮廓面像素点到平均面的距离之和，得到孔隙-裂缝双重网络模型中链表节点中的轮廓编号与指定点裂缝的外接长方体的轮廓编号相同的链表节点的粗糙度S a，其公式为：
其中，S a为粗糙度，l x为平均面的宽度，l y为平均面的长度，E(x，y)为被测量的轮廓面和建立的平均面之间的z坐标距离的算术平均。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310437515.7A CN116152259B (zh) | 2023-04-23 | 2023-04-23 | 一种基于图神经网络的储层渗透率计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310437515.7A CN116152259B (zh) | 2023-04-23 | 2023-04-23 | 一种基于图神经网络的储层渗透率计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116152259A CN116152259A (zh) | 2023-05-23 |
| CN116152259B true CN116152259B (zh) | 2023-07-04 |
Family
ID=86352873
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310437515.7A Active CN116152259B (zh) | 2023-04-23 | 2023-04-23 | 一种基于图神经网络的储层渗透率计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116152259B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117132870B (zh) * | 2023-10-25 | 2024-01-26 | 西南石油大学 | 一种CenterNet与混合注意力相融合的机翼结冰检测方法 |
| CN117611485B (zh) * | 2024-01-24 | 2024-04-02 | 西南石油大学 | 一种基于时空图神经网络的三维岩心渗透率预测方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104471614A (zh) * | 2012-05-18 | 2015-03-25 | 领英股份有限公司 | 使用数字岩石物理成像由岩石样品来评估岩石属性的方法和系统 |
| CN105487121A (zh) * | 2015-12-03 | 2016-04-13 | 长江大学 | 基于ct扫描图像与电成像图像融合构建多尺度数字岩心方法 |
| CN112419250A (zh) * | 2020-11-13 | 2021-02-26 | 湖北工业大学 | 路面裂缝数字图像提取、裂缝修补及裂缝参数计算方法 |
| CN113963130A (zh) * | 2021-10-25 | 2022-01-21 | 中国石油大学(华东) | 一种针对岩心裂隙的裂隙网络模型的构建方法 |
Family Cites Families (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9135502B2 (en) * | 2009-05-11 | 2015-09-15 | Universitat Zu Lubeck | Method for the real-time-capable, computer-assisted analysis of an image sequence containing a variable pose |
| CN109389128B (zh) * | 2018-08-24 | 2021-08-27 | 中国石油天然气股份有限公司 | 电成像测井图像特征自动提取方法及装置 |
| CN109887083A (zh) * | 2019-01-29 | 2019-06-14 | 中国石油集团测井有限公司西南分公司 | 一种裂缝-孔隙双重介质耦合渗透率模型的建立方法 |
| CN109900617B (zh) * | 2019-03-21 | 2022-06-07 | 西南石油大学 | 一种基于声电成像测井图的裂缝性地层渗透率曲线计算方法 |
| US20210097390A1 (en) * | 2019-10-01 | 2021-04-01 | Chevron U.S.A. Inc. | Artificial learning fracture system and method for predicting permeability of hydrocarbon reservoirs |
| CN110853138B (zh) * | 2019-11-21 | 2023-08-18 | 科吉思石油技术咨询(北京)有限公司 | 双重介质碳酸盐岩孔隙-裂缝双重网络模型构建方法 |
| CN111583148A (zh) * | 2020-05-07 | 2020-08-25 | 苏州闪掣智能科技有限公司 | 基于生成对抗网络的岩心图像重构方法 |
| CN111815583B (zh) * | 2020-06-29 | 2022-08-05 | 苏州润迈德医疗科技有限公司 | 基于ct序列图像获取主动脉中心线的方法和系统 |
| CN112102487B (zh) * | 2020-09-17 | 2022-02-15 | 西南石油大学 | 一种基于多重混合分形的非常规储层三维渗透率确定方法 |
| CN113359212B (zh) * | 2021-06-22 | 2024-03-15 | 中国石油天然气股份有限公司 | 一种基于深度学习的储层特征预测方法及模型 |
| CN113947041A (zh) * | 2021-10-14 | 2022-01-18 | 长江大学 | 一种页岩储层压裂缝网络扩展流动一体化模拟方法 |
-
2023
- 2023-04-23 CN CN202310437515.7A patent/CN116152259B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104471614A (zh) * | 2012-05-18 | 2015-03-25 | 领英股份有限公司 | 使用数字岩石物理成像由岩石样品来评估岩石属性的方法和系统 |
| CN105487121A (zh) * | 2015-12-03 | 2016-04-13 | 长江大学 | 基于ct扫描图像与电成像图像融合构建多尺度数字岩心方法 |
| CN112419250A (zh) * | 2020-11-13 | 2021-02-26 | 湖北工业大学 | 路面裂缝数字图像提取、裂缝修补及裂缝参数计算方法 |
| CN113963130A (zh) * | 2021-10-25 | 2022-01-21 | 中国石油大学(华东) | 一种针对岩心裂隙的裂隙网络模型的构建方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116152259A (zh) | 2023-05-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116152259B (zh) | 一种基于图神经网络的储层渗透率计算方法 | |
| CN110853138B (zh) | 双重介质碳酸盐岩孔隙-裂缝双重网络模型构建方法 | |
| CN106981093B (zh) | 一种分区约束耦合的三维地层并行建模方法 | |
| CN110599506B (zh) | 一种复杂异形曲面机器人三维测量的点云分割方法 | |
| CN105844602A (zh) | 一种基于体元的机载lidar点云三维滤波方法 | |
| CN109145129B (zh) | 基于层次三元组损失函数的深度度量学习方法及其装置 | |
| CN109859114A (zh) | 基于局域平滑性和非局域相似性的三维点云修复方法 | |
| CN107330734A (zh) | 基于Co‑location模式和本体的商业地址选择方法 | |
| CN104038792B (zh) | 用于iptv监管的视频内容分析方法及设备 | |
| CN113255892B (zh) | 一种解耦合的网络结构搜索方法、设备及可读存储介质 | |
| CN114998890B (zh) | 一种基于图神经网络的三维点云目标检测算法 | |
| CN107818338B (zh) | 一种面向地图综合的建筑物群组模式识别的方法及系统 | |
| CN113255677B (zh) | 一种岩体结构面及产状信息快速提取方法、设备及介质 | |
| CN113128518B (zh) | 基于孪生卷积网络和特征混合的sift误匹配检测方法 | |
| CN115309846B (zh) | 一种基于平行系数的道路网结构识别方法 | |
| CN110472543B (zh) | 一种基于局部连接特征匹配的机械图纸对比方法 | |
| Kieler et al. | Matching river datasets of different scales | |
| Zhang | Classification of Urban Land Use Based on Graph Theory and Geographic Information System. | |
| CN116452826A (zh) | 基于机器视觉的遮挡情况下煤矸石轮廓估计方法 | |
| CN115375925A (zh) | 一种基于相位信息和深度学习的水下声呐图像匹配算法 | |
| CN114937215A (zh) | 一种城市功能区的识别方法及装置 | |
| CN112818551A (zh) | 一种有效提高地理探测器模型精度的新算法 | |
| CN112183879A (zh) | 一种城市功能区的分类方法及装置、电子设备和存储介质 | |
| CN113033599B (zh) | 基于空间随机森林算法的露头地质体岩层分层方法 | |
| CN107341151A (zh) | 图像检索数据库生成方法、增强现实的方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |