一种油藏地质剖面图矢量化方法

段友祥，李宁宁，孙歧峰，李 钰

(中国石油大学(华东) 计算机与通信工程学院，山东 青岛 266580)

0 引 言

地质剖面图是用规定的符号、花纹和颜色,按一定的比例，沿一定的方向，表示一定的距离内地下一定深度内地质现象的图件,一般要表示出地层分层、断层、岩性等地层结构和岩体属性特征。它形象直观地表达了地层的结构构造和地层的沉积规律，是地层在垂向上最直观、有效的表达方式。油藏地质剖面图能为系统分析区域或局部的油藏地质条件、正确指导油藏资源开发利用和优化管理提供依据。

现在很多老油田已经进入了后期开发，小油藏、隐蔽油藏等复杂油藏成为主要的开发对象，运用的技术手段也越来越多样化，要求的基础地质数据资料也越来越多。但是，由于受早期条件的限制，这些老油田的很多地质资料不全或不符合现在开发技术手段的要求。例如，作为重要地质成果图件的地质剖面图，过去大多采用手工绘制(很容易转化为位图图像)或软件绘制位图图像(又称为栅格图像)。它们虽然能够逼真展现地层剖面，但是难以运用现代计算机技术进行处理，图像的数据量大，不能进行保真缩放，难以满足现代的钻井导向、油藏精确描述识别等先进技术的要求。因此，对过去陈旧的地质资料进行深度加工和处理是老油田有效开发的重要基础，矢量化就是其中重要的内容之一。矢量化是将油藏地质剖面图等栅格图像转换为矢量图像的过程，其实质是将图像数据转换为矢量图形数据，从而更方便地对图形进行处理，并精确地表达空间对象的位置、长度等信息[1]。

1 油藏地质剖面图矢量化方法研究

随着计算机技术的发展以及图像处理技术应用的深入，图像矢量化对象由早期的二值位图(即黑白图片)变为现在的彩色图像，矢量化算法也在不断发展。图像矢量化的流程一般包括图像分割、图像轮廓提取、曲线特征点提取、曲线拟合等主要环节。为了更好地实现矢量化，针对各个环节存在的突出问题，人们研究和提出了各种不同的解决思路，从而产生了很多不同的矢量化算法，主要包括基于轮廓的矢量化算法、基于边缘的矢量化算法、基于细化的算法、基于梯度网格的矢量化算法等。

钱晓峰等[2]提出了一种轮廓点列快速矢量化算法,极大地减少了轮廓需存储的像素点个数,节省了内存空间,并为进一步处理如形状匹配、编码等提供了基础。黄雪优[3]采用Canny算法和基于阈值的算法提取边缘,对边缘进行编辑获取闭合的轮廓，根据轮廓点矢量变化率判断轮廓的特征点并完成轮廓分段,对分段轮廓进行拟合初步获得图像矢量化边缘。最后采用矢量图手动编辑的方式对获得的矢量图边缘做出修正,减少特征点个数,提高矢量图的精度。石文莹[4]提出了采用基于细化和自适应网格相结合的矢量化方法，解决了传统基于细化方法容易在交叉处产生畸变的情况。杨云等[5]提出了一种线目标半自动矢量化的方法。该方法通过在用户选择的线目标上加一个滑动窗,采用颜色空间转换、k-均值聚类和区域生长相结合的方法对窗口内的当前线目标进行自适应分割;再通过不断地沿前进方向移动窗口和对窗口内的线目标进行分割、细化及序贯跟踪,来完成线目标的矢量化。王红岩等[1]针对煤矿地质剖面图提出了采用不同的算法删除各种岩性符号，只保留各个地层的边界线，最后自动矢量化的算法。这些算法在实际问题解决中都得到了应用，一些矢量化软件采用了其中的部分算法，但是对于一些特殊问题，这些算法还存在不足。

文中在深入研究已有矢量化算法的基础上，针对油藏地质剖面图的特点，提出了一种基于颜色聚类结合改进的Canny边缘检测的油藏地质剖面图矢量化方法。

1.1 油藏地质剖面图矢量化流程

根据矢量化的一般流程，确定地质剖面图矢量化的流程，如图1所示。

图1 地质剖面图矢量化流程

1.2 颜色聚类分割

油藏地质剖面图的地层层位之间一般采用不同的图例填充，其中图例由颜色和图案组成。对图像首先进行颜色空间转换，然后可以采用颜色聚类消除层位之间的图例。

1.2.1 RGB转Lab

常见的颜色模式有RGB模式、HIS模式、HSV模式和Lab模式等。由于采用RGB模式时，各种颜色在数值上的差距不能很好地反应其色彩差异，于是先将图像由RGB模式转换到Lab模式。Lab模式是国际照明委员会于1976年制定的一种色彩模式，是一种均匀的颜色空间，也是经常用来描述人眼可见的所有颜色的最完备的色彩模型。Lab颜色空间由L、a和b三个坐标轴组成。L表示亮度，值域是从0到100逐渐由黑变白；a和b的取值范围都是-128～+127，a值从小到大的变化是指从绿色变为红色，b值从小到大的变化是从蓝色变为黄色。一般先将RGB转化到XYZ空间，然后再由XYZ空间转化到Lab色彩空间，见式1～3。

(1)

(2)

(3)

设在Lab色彩空间模式下的两种颜色分别为C1和C2，则C1与C2的颜色色差定义为在L、a、b三个通道的欧氏距离，即：

ED(C1,C2)=

(4)

1.2.2 确定初始聚类中心

假设油藏地质剖面图的图例库为Library(其中Library是由图例填充时的颜色值L、a、b组成，L、a、b分别为Lab颜色模型的三个分量)，根据图例库Library分析油藏地质剖面图地层轮廓的填充颜色，从而确定初始聚类中心的颜色值。对于给定的油藏地质剖面图Image(RGB转Lab后的图像)中的每一个像素点p,计算其与每种图例颜色Library(C)(C为图例库中的一种颜色)的欧氏距离，每个像素点取对应最小欧式距离的点的颜色值建立颜色量化图像Image2,即：Image2=Library(C)。其中任意C1≠C，满足：ED(Library(C),Image(p))

1.2.3 颜色聚类

已知k个初始颜色聚类中心,循环执行以下两步：(1)对每个像素点p，计算其与每个初始聚类中心的欧氏距离，并将该像素划分到距离最近的聚类中心上；得到该像素的类别Classp∈{1,2,…,k},即将像素点p划分到第Classp类。(2)重新计算聚类中心。新的聚类中心为图例库Library中与类别平均色最接近的颜色。对于每个类别,先计算该类所包含的所有像素的平均色，然后计算该平均色与Library中的每一种颜色的距离，将距离最小值设为新的聚类中心；直到聚类中心未发生变化或达到最大迭代次数。

1.2.4 平 滑

颜色聚类后图像中还存在一些孤立的颜色块。采用图的最小割算法进行平滑处理，通过最小化能量函数[6]实现，重新分配各像素的类别。最终得到颜色聚类分割后的油藏地质剖面图。

1.3 改进的自适应Canny边缘检测算法

边缘检测是图像处理和计算机视觉中的重要组成部分。边缘检测的目的是识别包括图像特征信息和形状质量的大部分的亮度变化。众所周知的经典边缘检测算子包括Log、Robert、Prewitt、Laplace、Sobel和Canny[7]等。不同的算子有不同的优势。由于Canny算子具有更好的信噪比和检测精度，因此，Canny算子成为评估其他方法的基准，其基本步骤如下：

(1)利用一维高斯函数，对图像进行低通平滑滤波；

(2)用一阶偏导有限差分计算平滑后图像中各点的梯度值和梯度方向；

(3)对梯度幅值进行非极大值抑制；

(4)用双阈值算法检测和连接边缘。

但Canny算子也存在如下缺点：

(1)传统的Canny边缘检测利用高斯函数进行平滑滤波时会过度平滑图像，使图像变得模糊，可能扭曲和丢失弱边缘，特别是边缘的连接和拐角处；同时，高斯函数的方差需要人为设定；

(2)在计算像素点的梯度幅值时，对噪声较敏感，容易检测出假边缘或丢失一些真实边缘的细节[8]；

(3)在高低阈值的选取时，由于高低阈值不是由图像边缘的特征信息决定的，而是需要人为事先设定好的，因此不具有自适应性。

针对以上缺点，Cheng Y[9]用非线性扩散滤波器代替高斯函数平滑图像，求梯度幅值时考虑了像素对角线方向，最后用Otsu算法自动获取阈值,在边缘检测中取得了更好的精度。Zhao M等[10]针对传统Canny方法使用高斯函数去噪时，平滑度取决于高斯核的大小且只对于高斯噪声效果好的缺点，引进了DCT压缩域的概念，提出通过DCT变换将2D离散像素转换为连续DCT域，得到DCT系数和估计噪声，对DCT系数进行修正，然后通过IDCT获得平滑图像的方法，保留了更多有效的边缘信息。Ma X等[11]为解决传统Canny边缘检测算子过度平滑图像和适应性弱的问题，改进了参数Sigma和高低阈值的获取方法。Biswas R等[12]针对Canny算法对于模糊图像高低阈值难于确定的缺点，提出了基于type-2模糊集概念的算法来自动确定阈值。尚长春等[13]采用一种新的四阶偏微分方程的降噪算法对图像去噪,进一步提高降噪效果；采用自适应阈值的方法对图像边缘进行检测,实现了双阈值的自适应提取；基于模糊判决的理论,提出了一种有效的边缘连接方法。温阳东等[14]提出了一种基于改进非极大值抑制(NMS)过程和双阈值求取方法的Canny边缘检测算子，能够获得更好的图像边缘。

在研究以上改进算法的基础上，结合油藏地质剖面图矢量化的特点，即对地层水平层位轮廓矢量化，提出了对Canny算法改进的思路，步骤如下：

1.利用一维高斯函数，对图像进行低通平滑滤波(采用Sigma为1；Gaussian mask size 为5)。

2.用Sobel算子对图像进行卷积，求像素点的梯度，采用水平、45°和135°方向上的一阶梯度模板，如图2所示。

图2 一阶梯度模板

三个方向上的一阶梯度分量分别为Px(i,j)、P45°(i,j)、P135°(i,j)，由图2中对应的一阶梯度模板与图像进行卷积得到(其中i,j分别表示图像像素点的x坐标和y坐标)。每个像素点的梯度幅值和梯度方向计算分别为：

(5)

(6)

但是求完梯度后依然会存在一些噪声，采用如下方法消除部分噪声。

(2)Matrix为梯度图像中每一点的梯度平均值累加和矩阵，其中每一点存放了该点所经过的窗口的梯度平均值的累加和，Matrix1为梯度平均值的累加次数，计算Matrix中每一点的平均值为Matrix2。

(3)判断如果M中某一点的梯度值小于Matrix2中对应点梯度值的b倍(b为倍数因子)，则设该点的梯度值为0，否则该点的梯度值不变。(采用邻域窗口20行20列，step为5，倍数因子为2)。

3.对梯度幅值进行非极大值抑制。

4.用双阈值算法检测和连接边缘，其中的高低阈值选取方法采用改进的Otsu自动阈值选取方法[15]，连接边缘后，删除边缘长度小于等于2的边缘。

1.4 边缘轮廓矢量化

经过边缘轮廓提取后的油藏地质剖面图只包含水平的地层层位轮廓，需要将其矢量化，以便于后期的编辑和应用。采用Two-pass算法先对轮廓信息进行连通区域标记，然后对标记的每个区域采用Douglas-Peucker矢量数据压缩算法提取特征点(文中采用的DP阈值为1)，最后采用三次贝塞尔曲线对提取的特征点进行拟合，从而完成油藏地质剖面图的矢量化。

2 实验结果及分析

2.1 实验1

为了验证提出算法的可行性，实验在MATLAB中以Lena图像为例，并分别添加椒盐噪声(噪声密度为0.05)、高斯噪声(均值为0、方差为0．01)以及泊松噪声(噪声密度为0.05)来进行验证，并与传统经典算法Log、Canny边缘检测算法的检测效果作比较。从定量的角度分析，分别计算各边缘检测算子所对应的峰值信噪比(peak signal-to-noise ratio，PSNR)，如表1和表2所示。

表1 不同噪声环境下各算法对应的PSNR dB

表2 不同噪声密度下各算法对应的PSNR dB

上述实验结果表明，针对不同类型的噪声以及不同程度的高斯噪声，文中算法具有较好的抗噪性和准确性。

2.2 实验2

利用提出的油藏地质剖面图矢量化流程和改进算法，在MATLAB中用实际地质剖面图进行了矢量化实验，如图3和图4所示。

图3 实验结果1

实验过程和结果如下：

(1)采用传统的Canny边缘检测算法，对地质剖面图进行检测，结果见图3(b)和图3(e)。可以看出，图3(b)噪声很多，效果很差。采用改进的Canny边缘检测算法，对地质剖面图进行检测，结果见图3(c)和图3(f)。可以看出噪声少，较好地保留了水平的层位轮廓，效果好。

(2)对原始图先进行颜色聚类，结果见图4(a)和图4(e)；然后再采用改进的Canny边缘检测算法进行边缘检测，得到轮廓检测结果，见图4(b)和图4(f)；最后再进行矢量化，结果见图4(c)和图4(g)。

(3)对图4(b)和图4(f)采用美国斯坦福大学的Vector magic软件进行矢量化，结果见图4(d)和图4(h)。可以看出，文中矢量化结果与Vector magic基本相同，但是文中矢量化方法速度更快，得到的曲线特征点更少，更适合矢量化结果的存储和编辑。与文献[1]中针对每一种图例进行删除的算法相比，文中算法实用性更强。

图4 实验结果2

3 结束语

油藏地质剖面图中的岩性复杂多样，针对其特点提出了油藏地质剖面图的矢量化流程，增加了颜色聚类分割，改进了Canny边缘检测算法。实验结果表明，该方法很好地解决了油藏地质剖面的矢量化问题，矢量化速度快、质量高，为在随钻地质导向综合决策中应用油藏地质剖面奠定了基础。下一步包含断层的油藏地质剖面的矢量化将是研究的重点。