一种基于测井电成像图像的砂砾岩岩性自动识别方法

袁 阳， 侯国强， 马涛涛， 周 辉， 董汉强， 潘 涛

(1.中国石油集团测井有限公司吐哈分公司，新疆维吾尔自治区 吐鲁番 838200；2.延长油田油气勘探公司 石油勘探开发部，陕西 延安 716000)

砂砾岩油藏埋藏深，储集性能非均质性强，复杂的岩性导致砂砾岩地层的油气储集空间类型变化多样。复杂的岩石骨架削弱了储集流体的常规测井响应，造成常规测井信息识别、划分砂砾岩岩性、评价储层难度很大。

由于电成像测井资料分辨率高，可得到砂砾岩地层中砾岩颗粒形状、大小、空间位置及砂岩和泥质条带的直观描述，已经应用到砂砾岩储层的岩性识别研究中。但以往的研究主要是岩性的定性识别，表现为基于电成像图像的图像-岩性识别模式，需要解释人员根据模式人工做出岩性的判断、划分，没有实现计算机岩性自动定量识别[1-4]。目前，文献[5]采用数字图像处理技术对砂砾岩图像进行特征提取，再基于人工神经网络自动识别岩性；文献[6]提取了电成像图像中的均值、方差、分形维数、包络和非均质5条特征曲线，分别利用相关性和神经网络两种方法进行岩性自动识别；文献[7]通过图像灰度转化，构建灰度共生矩阵，计算样本的对比度、相关度、熵、均匀度和能量5个特征，利用贝叶斯判别分析法对岩性进行自动判别；文献[8]通过计算图像灰度值的平均值和方差，以及方差曲线的变差方差根，构建岩性的定量判别标准模板，进行岩性定量识别；文献[9]和文献[10]从电成像测井图像中提取定量表征砾石含量的结构指示参数，结合常规测井资料建立交会图，定量自动识别岩性。但这类基于图像模板特征提取的识别方法，未考虑图像的井眼覆盖率，图像中极板间的空白条带未充填，使得样本图像的特征计算结果不准确；并且这类方法依赖建立的图像标准模板，标准模板的有限性导致识别各种多变的未知图像时，出现准确率低、泛化能力差等问题。

本文以图像中砾岩颗粒提取为核心，首先采用MCA(Morphological Component Analysis)算法充填图像中极板间的空白条带，生成全井眼图像,再进行基于边缘流的多尺度砾岩颗粒边缘提取和基于曲线演化的砾岩颗粒区域分割，得到初步分割后的砾岩颗粒区域。但由于电成像测量技术、成像技术的限制和复杂井况，初步分割的颗粒区域与人工视觉判别结果有很大区别，因此在处理井段引入区块地质沉积状况时，特别是在将相应沉积条件下的砾岩颗粒磨圆度这一先验知识作为约束条件时，自顶向下地驱动砂砾岩颗粒区域的合并和分裂。最后，构建砂砾岩粒度谱和岩性精细剖面，得到砂砾岩岩性的精细评价成果。本文方法应用于YJ油田深层砂砾岩难动用油藏的测井评价中，为优选射孔井段和压裂层段提供了准确依据。

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

由于电成像测井在大多数井眼条件下都存在空白条带，达不到全井眼覆盖，为保证砾岩颗粒提取的准确度，采用MCA算法[11]对空白条带进行充填。MCA将图像分解为几何结构un和纹理ut两种形态成分，并假设两个形态成分在某个特定的基库或过完备字典下是稀疏的，且基库或过完备字典之间具有不相干性。MCA模型表示如下[12]：

(1)

图1 MCA算法具体实现步骤

图2 原始成像图及处理结果示例

MCA算法具体实现结果如图2(b)所示，可以看到图2(a)的原始图像经算法充填后，较好地保持了砾岩颗粒边缘的连续性。本文以每次迭代得到的充填图像与原始图像在掩膜模板上的平均像素灰度误差作为迭代结束的条件。表1给出了图2(a)原始图像分别在迭代5，10，20，30，40，50和100次时的平均像素灰度误差值。

由表1看到，算法迭代到20次时，误差急剧下降，随着迭代次数的增加，误差变化不大，说明算法迭代20次以后，充填结果基本相同，实验结果也证明了这一点。因此，通过适当地选择误差阈值(或两次迭代的差异值)，可以自适应地结束迭代，得到可以接受的充填结果，从而避免固定较大迭代次数的方法带来的计算冗余。

表1 不同迭代次数下的平均像素灰度误差值

2 基于边缘流的多尺度砾岩颗粒边缘提取

边缘流算法[13]在初始尺度计算边缘流能量，并将各个像素点的边缘流能量投影到水平、垂直方向，进行累加，得到初始边缘流向量场；随着尺度提升，当初始向量场幅度小于设定幅度阈值时，当前尺度的向量场等效于初始向量场；当初始向量场与当前向量场角度小于设定的角度阈值时，则将当前向量场叠加上初始向量场。在多个尺度重复上述过程，直至到达最大设定尺度，得到最终的多尺度边缘流向量场。最后，根据Possion方程，由最终的多尺度边缘流向量场的散度计算得到图像的边缘流函数，结果如图2(c)所示，边缘流函数通过归一化转换为灰度图像显示，可以看到灰度值较小的砾岩颗粒边缘，但有些边缘并不连续。

3 基于曲线演化的砾岩颗粒区域分割

基于曲线演化的砂砾岩颗粒区域分割算法[14]，是根据上节计算得到的边缘流函数，通过求解水平集方程，计算区域分割单像素曲线，其特点是可以将不闭合的断裂边缘进行封闭。水平集方程如下所示[14]：

(2)

式中，φ为迭代计算曲线演化结果图像；t为时间；S为边缘能量场；▽为取散度；w为平滑权重系数；C为边缘流函数；K为一个表示收缩或扩张的常数。最后对曲线演化结果图像分别以高斯函数水平、垂直方向的一阶导数为卷积核，进行卷积运算，将幅度值大于设定阈值的像素点选作边缘点，结果如图2(d)所示，可以看到蓝色的单像素宽度曲线封闭了各个砾岩颗粒区域。

4 区块地质沉积约束条件下的砂砾岩颗粒区域合并与分裂

从图2(d)中可以看到，各个区域被闭合的蓝色单像素曲线区分开来，但有的砾岩颗粒区域明显与人眼视觉判别结果不同，出现明显的陡峭凹陷。图3为砂砾岩颗粒区域分裂示例。从图3(a)的区域分割结果可以看出，同样各个区域被闭合的蓝色单像素曲线区分开来，但有的砾岩颗粒区域存在明显的粘连情况，两个邻近的砾岩颗粒被划分为一个区域。这些砾岩颗粒的形状特征出现的概率在当前地质沉积条件下很低。

分析当前区块地质沉积状况，由此得到砾岩颗粒磨圆度的统计规律，作为约束条件知识驱动颗粒区域的合并与分裂过程。由颗粒磨圆度判定凹陷的颗粒分割区域是否进行相邻区域的合并，以及明显粘连的颗粒分割区域是否进行区域的分裂，由一个区域分裂为多个区域，并将每个分裂后的区域作为单独的一个砾岩颗粒。具体方法是采用自适应曲率计算技术，定位砾岩颗粒区域中边界上的凹点、凹点对和邻近凸点。由邻近的非砾岩区域的相对位置，合并邻近区域与该砾岩区域。同时，根据定位的凹点对之间的弦长和计算的对应弧长，与弧弦比阈值比较，将该砾岩区域分裂为多个砾岩颗粒，这样就解决了初步分割区域的不完整和粘连情况。

图4为砂砾岩颗粒区域合并示例。图4(a)是图2(d)使用并查集提取的各个分割区域，得到各区域的包含像素点，并以不同颜色表示，根据砾岩灰度阈值，判定是否为砾岩颗粒。提取砾岩颗粒边界像素，并以顺时针排序。自适应计算边界曲率，定位凹点及其邻近凸点，凹点和其邻近凸点形成一黄色矩形区域，如图4(b)所示。如果邻近区域的质心在该矩形区域内，则合并邻近区域；如图4(c)中所示的左下方浅绿色表征的砾岩颗粒区域，经过区域合并处理后，形成了完整的砾岩颗粒区域，符合当前砾岩颗粒磨圆度的设定。图3(b)为图3(a)使用并查集提取的各个分割区域，计算砾岩区域各个凹点对之间的弦长以及对应的弧长，并根据凹点的次数权重，使用权重大于1的凹点对对应弦，切分该砾岩区域。图3(b)中大面积的深蓝色区域表明存在明显的粘连情况；如图3(c)中所示，该区域由红色的凹点对之间相连的弦，分裂为两个不同的区域(分别用深红色与深蓝色表征)，这样，每个区域的形状特征符合当前颗粒磨圆度的知识约束。

图4 砂砾岩颗粒区域合并示例

5 粒度谱构建及精细岩性剖面计算

砂砾岩粒度谱的构建和岩性精细剖面的计算采用滑动窗口方法，如图5所示，统计某深度点对应滑动窗口内，砾岩颗粒区域面积的大小，并由此划分区域面积的BIN值(面积大小的统计分组数目，每个BIN(Binary In Net)代表了各自的面积大小范围)，统计属于每个BIN值范围内的砾岩颗粒的数目，构成粒度谱，作为解释成果直观地展现砂砾岩颗粒的大小分布。如图6所示，精细岩性剖面计算首先根据砾岩颗粒大小阈值，将砾岩颗粒细分为细砾岩、中砾岩和粗砾岩，并根据滑动窗口内各岩性颗粒面积除以滑动窗口面积，得到精细岩性剖面中各细分岩性的百分比，而砂岩岩性百分比由砾岩颗粒占比的残余值计算，并形成曲线，由解释软件进行充填绘制，形成直观的砂砾岩精细岩性剖面。

图5 砂砾岩粒度谱构建原理

图6 砂砾岩精细岩性剖面计算原理

6 应用效果分析

砂砾岩岩性精细评价成果图及与NRA岩性识别方法对比如图7所示。

图7 砂砾岩岩性精细评价成果图及与NRA岩性识别方法对比

图7(b)给出了Y920井3555.517～3557.808 m井段的精细评价综合处理成果图，可以看到构建的粒度谱和精细岩性剖面，与图7(a)所示的NRA(Neutron Resistivity Acoustic)岩性识别方法[15]成果图进行对比可以看到，在解释井段，两者的对应性较好；解释井段上部砾岩总组分偏少，对应NRA方法的砾状砂岩；下部砾岩总组分逐渐增加，对应NRA方法的砾岩岩性。粒度谱直观地展现了相应井段砾岩颗粒的大小分布，岩性剖面相对于传统的砾岩、砾状砂岩和含砾砂岩的岩性划分，更加精细地区分了砾岩颗粒的大小，而不是笼统地根据井段中所有砾岩部分所占面积划分砂砾岩岩性，不考虑砾岩颗粒的独立性和特殊性。

图8为Y936井3676～3677.6 m井段处理成果图，岩性剖面中以细砾岩和中砾岩为主，对比常规解释结论(油层)，表明细砾岩、中砾岩的物性、含油性较好，符合YJ油田砂砾岩体区块的储层“四性”关系，据此优选压裂设计井段，压裂成功实施后，获得初产近20 t的产能，目前日稳产油10 t。

图8 Y936井3676～3677.6 m井段处理成果图

7 结束语

本文的岩性自动识别方法通过图像预处理以及在地质沉积约束条件知识驱动下的精细图像处理技术，有效地提取了图像上的砾岩颗粒，并据此构建粒度谱和精细的岩性剖面，为用户提供了直观、准确的砂砾岩评价成果。

本文方法可在今后的老井复查和新井解释中发挥重要作用，为开拓服务市场和油田的勘探开发提供重要的技术支撑。