基于知识的碳酸盐岩电成像测井裂缝自动识别方法

2019-05-28宋锦祥

测井技术 2019年6期

宋锦祥

(中国石油集团测井有限公司大庆分公司,黑龙江大庆163000)

0 引言

通过微电阻率扫描成像测井资料拾取和识别裂缝是碳酸盐岩储层评价中重要的工作之一。目前生产实际中,裂缝拾取和识别主要靠解释人员凭借经验手工完成,存在耗时、效率低下、工作量大和对经验要求较高等问题。现有的裂缝自动拾取和识别方法[1-3]主要是在图像预处理的基础上,进行Hough变换提取裂缝,针对图像环境较为简单,且不能识别裂缝类别。

本文借鉴和使用了人工智能技术,对解释专家裂缝手工拾取和识别经验以及裂缝区块总体知识,进行定量化表征,在信息融合的基础上,形成知识规则库,送入知识推理机,同时电成像图像经过图像预处理后得到的二值化图像以及进一步由图像细化和Hough变换拾取的正弦曲线组,作为推理数据,也送入推理机;由逆向推理机制,得到裂缝拾取和识别结果。

图1 裂缝自动拾取与识别方法处理流程图

1 改进的图像预处理

传统裂缝自动拾取和识别方法图像预处理步骤较为简单,主要有全局图像二值化和图像细化2个步骤。但全局二值化对于低对比度图像和局部明暗变化明显的图像,处理效果并不理想,其核心问题是全局阈值难以选择。因此本文在电成像空白条带充填的基础上,通过直方图匹配方法,提升图像的对比度;采用基于积分图像的自适应二值化步骤,计算自适应的二值化阈值;考虑照明的空间变化,增加局部明暗变化时的二值化处理的鲁棒性。对于二值化图像,本文预处理步骤采用基于曲率的区域分割和基于局部对比度的白色背景区域转换方法,进一步提升图像处理效果。图1为本文提出的裂缝自动拾取和识别方法的处理流程图。

1.1 电成像图像空白条带充填

由于电成像测井在大多数井眼条件下都存在空白条带,达不到全井眼覆盖,为保证裂缝拾取和识别的准确度,需要对电成像图像的空白条带进行充填。采用了基于形态分量分析算法(MCA)的空白条带充填方法[4];MCA方法将图像的结构分量和纹理分量进行迭代分解与重构,具体算法参见文献[4]。

1.2 基于直方图匹配的图像增强

对比度低的图像,其直方图覆盖较宽的灰度级范围并且分布较为均匀。直方图较均匀分布在0～255灰度范围内,此时做图像二值化,灰度阈值选择困难,不得不进行手工试错。而高对比度的图像,其直方图在低灰度范围和高灰度范围呈明显的双峰分布。直方图波峰之间波谷处的灰度可作为理想的图像二值化灰度阈值。因此,为增加图像对比度,需要把灰度级的分布拉开,

使灰度级层次分明。本文采

用了直方图匹配方法,将原图像的灰度直方图规定化到预先设定的双峰高斯函数,进行图像增强。

1.3 基于积分图像的自适应二值化

常规图像二值化采用固定的全局阈值,对于图像局部的明暗变化处理效果不好。采用基于积分图像的自适应阈值方法[5],对图像进行二值化处理。

对于一幅灰度的图像,积分图像中的任意一点(x,y)的值是指从图像的左上角到这个点的所构成的矩形区域内所有的点的灰度值之和,计算公式为

(1)

式中,IntegralImg(x,y)为像素点(x,y)的积分图像值;Img(x,n)为点(x,n)处的灰度值。

对像素点(x,y)二值化时,比较该像素点的灰度值与大小为s×s的邻域块内的积分图像平均值RegionValue

RegionValue=IntegralImg(x2,y2)-

IntegralImg(x2,y1-1)-IntegralImg(x1-

1,y2)+IntegralImg(x1-1,y1-1)

(2)

式中,x1,y1和x2,y2分别为大小为s×s的邻域块左上角和右上角坐标。

RegionValue作为自适应二值化的局部阈值,与像素点(x,y)的灰度值比较,如果(x,y)的灰度值小于平均值RegionValue,则将像素点(x,y)设置为黑色,否则被设置为白色。

1.4 基于曲率的区域分割

自适应二值化后,在局部会出现距离非常近的黑色目标区域断裂,为提升二值化图像处理效果,使其更符合人眼的视觉图像分割结果,采用基于曲率的白色背景区域分割方法,连接断裂的黑色目标区域。具体算法:①使用并查集算法,提取各个白色背景区域,得到各区域的包含像素点及边界像素点;②计算边界像素点的曲率,提取曲率为负值的局部极大值点作为凹点;③根据凹点及曲率计算时用到的2个边界点,计算法向向量方向;④在一定图像距离内,搜索与上述凹点法向向量方向相对的另一凹点,如果存在,则将2个凹点标记为分割凹点对,据此,将凹点对用黑色前景像素连接。

在上述预处理步骤后,还采用了基于局部对比度的白色背景区域转换,将区域面积小,且与包围其区域的局部黑色前景区域对比度相差较小的白色背景区域转换为黑色前景像素;在此基础上,进行图像的细化。在经过图像预处理后,本文采用改进的Hough变换[1]进行裂缝自动拾取。

2 基于知识的裂缝识别方法

2.1 裂缝区块的总体统计知识及定量表征

研究区块属于TH油田奥陶系油藏,成像测量井段主要为石炭系巴楚组和奥陶系一间房组地层。储层总体特征是自然伽马相对低值,电阻率相对低值,部分储层段双侧向呈正幅度差,中子、声波数值变化较小。溶蚀发育,或溶蚀及裂缝均较发育,多数高导缝为溶蚀孔洞沟通形成,部分连通性较差,静态图像暗黄色为主,储层为溶蚀孔隙型储层及裂缝孔隙型为主。诱导缝不发育,井眼崩落不明显。主要岩性为泥晶灰岩、泥质灰岩和灰质泥岩;其中泥晶灰岩段电成像静态图像黄色或暗黄色,高导缝较发育、溶蚀较发育,层理发育较少。泥质灰岩和灰质泥岩段电成像静态图像暗色,层理发育。

因此,本文研究的重点集中在识别、区分高导缝和层理。从0°～90°,以2°为间隔,统计区块高导缝倾角分布频率,归一化后形成高导缝倾角的先验概率统计模型,作为定量化表征的区块总体统计知识。

2.2 解释专家裂缝人工识别知识及定量化表征

2.2.1碳酸盐岩地层电成像裂缝人工拾取知识

在生产实际中,碳酸盐岩电成像解释主要拾取高阻缝、高导缝、诱导缝和层理这4种类型。

(1)正弦线显示为高阻亮黄色、黄色形态,可以直接拾取为高阻缝;形状不连续的也可以拾取为高阻缝。

(2)呈雁列状排列的,或者是180°对称排列的可以拾取为诱导缝;诱导缝中充满低阻钻井液,在电成像图中显示为黑色。

(3)正弦线连续、完整且成组出现,相互平行或接近平行,拾取为层理;一般为低电阻率泥岩夹层,表征为黑色。

(4)正弦线随时中断,由溶蚀现象沟通的,且组内正弦线既不平行,也不规则的,拾取为高导缝,是否为泥质充填或开启缝先不予考虑;高导缝由于泥质充填或裂缝开启,电阻率低,表征为黑色。

(5)低角度正弦线,有溶蚀沟通,且与层理方向不一致的,可以拾取为低角度高导缝,同样为黑色颜色特征。

2.2.2知识的产生式规则化表述及定量表征

裂缝人工拾取知识可以用IF-THEN产生式规则表述,如IF Condition1∧Condition2 THEN Action1;其中Condition1:正弦线连续、完整;Condition2:正弦线成组出现,相互平行或接近平行;Action1:拾取层理。但在实际推理过程中,需要将条件定量表征,以实现计算机自动计算。

总结3个主要条件:正弦线连续性;正弦线与其连通的黑色前景区域上下边界的相似性(判断是否为层理边界或溶蚀孔洞);空间相似度(正弦线组内倾角、倾向相似度)。

正弦线连续性用正弦线在二值图像中黑色前景区域内的点数与正弦线总点数之比计算,其数学表征为

(3)

式中,NumP(SinCurveinBlack)为正弦线在二值图像中黑色前景区域内的像素点数;NumP(SinCurve)为正弦线的总像素点数。

上下边界相似性的数学表征为

BoundSimilarity=

upperBlackBound,scale)+

lowerBlackBound,scale)

(4)

式中,SinCurveinBlack为正弦线在黑色前景区域内的曲线段;upperBlackBound为SinCurveinBlack曲线段所占黑色前景区域的垂直方向上对应的上部边界曲线;lowerBlackBound,scale为SinCurveinBlack曲线段所占黑色前景区域的垂直方向上对应的下部边界曲线;Similarity函数为曲线在每个尺度scale上的相似性度量;Normalize为归一化函数。本文采用多尺度曲线相似度度量方法,在正弦线贯穿的黑色前景区域内,在多个尺度上,通过比较正弦线与上下边界的导数大小,计算曲线相似度。

空间相似度度量的数学表征为

SpatialSimilarity=

(5)

式中,SinCurvei为要度量的第i条正弦线;DASimSinCurvej为不为i的,且与SinCurvej具有相似倾角、倾向的第j条正弦线;n为拾取的正弦线总数。本文通过Distance函数,计算SinCurvej与DASimSinCurvej基线之间的距离,并应用距离越近,空间相似度越大的概念,将距离的倒数之和作为空间相似度度量。

图2 知识规则库框架

2.3 信息融合

裂缝区块的先验概率统计模型与上述3个主要条件都可以作为单独的分类器,对拾取的裂缝进行识别和分类。采用多特征信息融合技术[6],在度量层将4个分类器的不同置信度量综合成1个带有更多信息量的置信度量,来做出更加精确的决策。表示为

(6)

式中,g为融合输出;wk为赋予第k个分类器的权值;fk分别为区块先验概率、正弦线连续性度量值、上下边界相似性度量值和空间相似性度量值。

由于区块的裂缝先验概率统计模型表征了区块的地质构造特征,其高导缝分布符合一定的地质规律,故信息融合中,区块先验概率的融合权重应占主要部分。而由解释专家经验知识可知,正弦线连续性度量、上下边界相似性度量和空间相似性度量具有同等的重要性,故其融合权重应为同等比重。因此,在上述原则下,设定区块先验概率权重w1=0.4,其他3个分类器权重均为0.2;融合权重设定值的准确性在第3节的实验中得到验证。

2.4 规则库构建及逆向推理

产生式规则描述逻辑推理过程是很有效的,但对于大量的实际应用,特别是在测井解释过程中,逻辑推理和数值计算交叉混杂,纯产生式规则不能很好地描述某些领域的专家知识。需要对纯产生式规则作必要的改进,采用关联函数F把规则的前提和结论等2个部分有机地关联起来。例如,产生式规则a,b,c,…⟹r可扩充为AND(a,b,c,…)⟹r;每条扩充的产生式规则中都含有节点(变量)、关联函数和变元,故将这种扩充的产生式规则简称为NFA语言[7-8],即F(a,b,c)⟹r又可以写成r=F(a,b,c)。

关联函数F可以是算术运算函数,或者是逻辑推理函数,也可以是经验公式,甚至是1个内部进行复杂运算的子程序。把表达某个领域知识的全部扩充的产生式规则连接起来,便构成完整的推理网络,形成所要建立的知识规则库(见图2)。

图2中HRFrac,InduceFrac,HCFrac和Layer分别为高阻缝、诱导缝、高导缝和层理的布尔值表示,w1～w4为多特征信息融合计算使用的权重,作为解释参数;DIP为计算的裂缝倾角,BWImg为图像预处理得到的二值图像,FracGroup为图像上拾取的正弦曲线组,这些是在推理之前计算得到的,作为交换信息,放入公共数据黑板[9],供推理过程使用。CalcCon,CalcBsim,CalcSSim分别为正弦线连续性,上下边界相似性和空间相似性计算子例程。PriorProbModel表示在区块先验概率模型中寻找DIP值对应概率的函数。

采用目标制导下的逆向推理机制,其实质是以深度优先为基调的搜索过程。以当前规则结点变量的计算为目标,先假定结论,再依次计算产生式规则的条件是否满足,再修改和确认结论。针对每条正弦线,推理都始于顶结点Frac,依次计算高阻缝、诱导缝、高导缝和层理的布尔值,如果高阻缝和诱导缝不满足条件,则计算是否为高导缝。由推理网络,高导缝的计算向下展开。推理机制模拟人工识别过程,对拾取裂缝进行识别和分类。

3 裂缝拾取和识别结果对比分析

3.1 裂缝区块的先验概率模型

统计了区块内7口井电成像测井段高导缝手工拾取结果,由倾角频数的对应关系,归一化后得到倾角的概率分布图,作为先验概率模型。从模型可知,区块内高导缝的倾角主要集中在70°～80°之间。

3.2 图像预处理结果对比

在测试井中选取5 529.6～5 532.2 m井段的动态电成像图像进行处理[见图3(a)];图3(b)显示空白条带充填后图像;图3(d)为图像增强处理结果;经过本文的图像预处理方法,得到二值化图像,见图3(e);图3(f)为对灰度图像取全局阈值(阈值为108),得到的二值化图像。可以看到,本文方法处理后的二值化图像更加符合人眼视觉的图像分割结果。

3.3 裂缝的自动拾取和识别结果对比

对二值化图像[见图3(e)]细化后,采用文献[1]所述的改进Hough变换方法,首先由细化图像中所有水平距离为T/2的白色像素点配成点对,以图像纵坐标建立一个一维参数空间,并将点对的中点纵坐标在该参数空间投票,检测参数空间的峰值位置,确定正弦线的基线。接着建立正弦线幅度和初始相位的二维Hough变换参数空间,将每个白色像素点的坐标(x,y)映射到该二维空间,同样检测峰值位置,就可得到拾取正弦线的基线、幅度和初始相位。拾取的正弦曲线见图4(b),红色、蓝色和绿色分别表示高角度、中角度和低角度正弦线;规则库推理过程中,多特征信息融合计算结果见表1。

图3 图像预处理效果对比图

由表1的融合计算结果,及规则库推理,可得到正弦线的识别结果[见图4(c)],红色为高导缝,绿色为层理。其中,采用图像后处理方法,对层理的上下边界进行正弦曲线拟合,基线距离大的上下边界拟合为新的层理,同时删除原层理。裂缝自动拾取与识别结果,与图3(d)所示的手工拾取结果比较,基本相符。但手工拾取结果中在5 531 m处未出现层理,对比常规测井数据(见图5),可以看到在5 531 m处,自然伽马相对高值,电阻率值明显下降,结合电成像数据,说明此处应为泥质条带,识别出层理是正确的。

图4 裂缝自动拾取与识别结果对比

表1 多特征信息融合计算结果

正弦线①正弦线②正弦线③正弦线④正弦线⑤正弦线⑥正弦线⑦正弦线⑧倾角/(°)85.95567.96064.84513.34821.5774.52213.34817.556先验概率0.20790.39080.35920.0021750.0173800.0021750.05653连续性0.7180.7720.8861.00.6981.00.8710.985边界相似度0.0690.3870.250.8540.4290.8780.60.9空间相似度0000.2330.30.2880.3820.354融合结果0.5200.5240.5160.1830.3220.1670.2300.175