一种基于纹理特征提取的地震图像检索方法

技术领域

本发明涉及地震信息处理技术领域，具体涉及一种基于纹理特征提取的地震图像检索方法。

背景技术

图像处理技术和计算机视觉方法在地震解释方面的潜力是无可争议的，并且这方面的研究还在逐渐增加；近年来，图像检索应用越来越普遍，人们可以通过智能手机、可联网数码相机等电子设备，对自己感兴趣的文本或图像信息进行检索；图像检索已在国防建设、医疗、通讯、机器视觉等许多领域得到广泛应用；同样，对地震图像检索的需求也日益增多，所谓地震图像检索，就是从已有地震图像(地震剖面、地震切片等)提取其信息，包括同相轴的形状及其振幅、相位、频率等参数，进而依据这些信息即可在解释专家库的众多经典解释成果中检索到最相近的剖面，从而帮助解释人员得到正确的判断；地震图像检索包括对地震数据的分析，以生成关于地层性质和结构的合理模型和预测；由此可知地震图像检索的重点是地震图像信息的提取；

在数据采集和后期处理之后，地震结构或地震剖面可以通过其纹理特征在视觉上进行区分；基于文本的图像检索依赖于图像元数据的注释质量，其标注往往会受到人为的主观理解片面或检索条件的误差，导致信息不对称；而基于内容的图像检索(Content-based Image Retrieval，CBIR)则使用基于计算机视觉的工具，根据图像的颜色、形状、纹理、边缘等底层视觉特征，建立特征向量，采用高维索引技术，根据距离计算方法计算特征向量之间的距离，进行相似度匹配；目前大量研究人员对基于内容的图像检索展开研究，实现CBIR的关键技术之一是图像特征提取；而在图像的底层特征中，纹理特征由于提取难度较大，一直是研究领域中的热点，纹理特征是地震图像的主要特征之一；

在地震纹理特征分析和提取方面，研究最广泛的是灰度共生矩阵 (GLCM)、局部二值模式(LBP)和Gabor滤波器三种纹理描述子；作为最常用的一种统计分析方法，基于灰度共生矩阵算法除了在提取地震纹理属性的断层、裂缝和储层特征的描述以及研究区域相带的分类中发挥着作用，在图像检索中也有着广泛的研究和应用，其优点是原理简单，易于理解，且能较精确地反映纹理信息，但由于其不具有旋转不变性，而待检索的图像很多是人工拍照上传到数据库，图像在方向上有很大的随意性，因此使用GLCM进行检索会影响检索的准确率；LBP是一种有效的纹理特征提取算法，主要用于解决多分辨率和旋转情况下图像的纹理分类问题，从某种意义上说，该方法将局部纹理信息与全局纹理信息同时融合到一起，但由于LBP对噪声敏感、对光照变化稳定性弱，所以当噪声出现时，会有较多的伪边缘、弱边缘的问题；基于信号处理的Gabor滤波器可以对图像进行多尺度、多方向变换，拥有时域和频域的最佳分辨率，能够将一幅图像从实数域变换到复数域，但是Gabor滤波器的计算性能完全取决于相关参数的选择；因此，对纹理描述符的性能研究已经成为许多领域的研究热点，在地震分析中，纹理描述符被用于计算地震图像的附加属性，从而为解释者提供新的数据可视化，或辅助半自动分类和聚类技术；这时，纹理特征往往被认为是简单的二阶地震属性，并且，在目前公开的文献中，很少有论文纯粹集中在计算机视觉(ComputerVision，CV)处理，也未有相关公开的利用纹理特征提取地震图像的方法来解决上述问题。

发明内容

针对上述存在的问题，本发明旨在提供一种基于纹理特征提取的地震图像检索方法，本方法通过将LBP算子和Gabor过滤器引入地震图像检索领域，对图像纹理特征进行提取，能够有效提取图像中的纹理特征，进而提高地震图像检索过程的准确性，具有检测精确性高和鲁棒性好的特点。

为了实现上述目的，本发明所采用的技术方案如下：

一种基于纹理特征提取的地震图像检索方法，包括步骤

Step1.在图像特征数据库中，对切片图进行分割和建立数据库；

Step2.接收需要查询的地震图像；

Step3.通过提取查询图像和数据库中的纹理特征，根据相似性度量对提取出的特征向量进行比较，计算出视觉相似度；

Step4.在数据集进行样本检索，检索出与查询图像视觉相似性最高的前N个最相似的候选图像块；

Step5.在获得候选图像块之后，采用欧几里得距离进行图像的相似性度量。

优选的，步骤Step1所述的对切片图进行分割和建立数据库的过程包括

Step101.将图像特征数据库中的每张切片图分割为不同大小的图像块；

Step102.根据分割后的图像块像素分类，为每个图像块指定一个唯一标签表明其所属分类，建立数据集。

优选的，在步骤Step2所述的地震图像输入之前，每张灰度图像要经过一定的预处理步骤，将图像强度重新调整为64个灰度级，然后，进行直方图均衡化。

优选的，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用局部二值模式对图像空间域纹理特征进行提取。

优选的，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用基于信号处理的Gabor滤波器纹理提取方法进行纹理特征提取。

优选的，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用基于统计的GLCM纹理结构描述方法进行纹理特征提取。

优选的，步骤Step4所述的样本检索的过程包括

Step401.在检索过程中，输入查询图像块，将每个查询图像块与数据库中所有图像块进行比较，以评估地震图像间的视觉相似度；

Step402.然后，从数据库中返回前N个最相似的候选图像块。

优选的，所述检索方法还包括步骤Step6.计算与查询图像类别相关的检索样本与总检索样本之间的比率，进行检索准确度分析。

优选的，所述检索方法还包括步骤Step7.每分类准确度分析，

在所有查询图像块的检索准确度检索的基础上，进一步分别针对每个类别的图像块进行计算，进行每分类准确度分析，并输出最终检索结果。

本发明的有益效果是：本发明公开了一种基于纹理特征提取的地震图像检索方法，与现有技术相比，本发明的改进之处在于：

本发明提出了一种基于纹理特征提取的地震图像检索方法，本方法通过将LBP算子和Gabor过滤器引入地震图像检索领域，对图像纹理特征进行提取，能够有效提取图像中的纹理特征，进而提高地震图像检索过程的准确性，具有检测精确性高和鲁棒性好的优点；同时为地震图像检索、分割和分类方面的研究提供理论依据。

附图说明

图1为本发明基于纹理特征提取的地震图像检索方法的流程图。

图2为本发明实施例3像素相对分布示意图。

图3为本发明实施例4地震图像3D数据体图。

图4为本发明实施例4Netherlands数据集的示例图。

图5为本发明实施例4提取自Netherlands数据集的横测线的两个不同分类C1类和C2类的图像图。

其中：在图5中，C

具体实施方式

为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。

实施例1：参照附图1-5所示的一种基于纹理特征提取的地震图像检索方法，包括步骤

Step1.图像特征数据库中，对切片图进行分割和建立数据库；

Step101.将图像特征数据库中的每张切片图分割为不同大小的图像块；

图像分割均采用矩阵分割的方式，按照像素尺寸分割成40Pix、 50Pix、60Pix等大小，类似于图5所示；

Step102.根据分割后的图像块像素分类，为每个图像块指定一个唯一标签表明其所属分类，建立数据集；具体来说，若某个图像块至少75％的像素属于分类A，则该图像块将被指定标签A；若某个图像块包含超过25％的噪声，则其将被丢弃，建立数据集；

Step103.使用步骤Step101得到的不同大小的图像块，以估计纹理描述子在处理不同尺寸的图像时的鲁棒性；

Step2.接收需要查询的地震图像；

Step3.通过提取查询图像和数据库中的纹理特征，根据相似性度量对提取出的特征向量进行比较，以计算出视觉相似度；

Step4.在数据集进行样本检索，检索出与查询图像视觉相似性最高的前几幅图像；

Step401.在检索过程中，给定查询图像块，将每个查询图像块与数据库中所有图像块进行比较，以评估地震图像间的视觉相似度；

图像的相似性度量：采用欧几里得距离和曼哈顿距离，在本实施例实现的算法中，欧几里德距离和曼哈顿距离都用于比较、查找、查询图像与数据库中图像之间的相似性；其具体计算过程见公式(1) 欧式距离和曼哈顿距离；需要遍历所有样本；通过实验发现，本实施例采用欧式距离和曼哈顿距离进行相似性度量时，产生的度量结果相似，因此为了简单起见，本实施例实验数据仅展示了欧式距离的度量结果；

Step402.然后，从数据库中返回前N个最相似的候选图像块；

Step5.图像的相似性度量

在经步骤Step402获得候选图像块之后，采用欧几里得距离进行图像的相似性度量，其中：

式中：

d

两数据集(向量)a，b都是n维，其中，a＝[x

Step6.计算与查询图像类别相关的检索样本与总检索样本之间的比率，进行检索准确度分析

由于每个图像块仅与单个类别相关，实际上可确认在检索到的前 N个样本中，通过分析描述子检索相似区域的能力，以及计算特征向量的时间，来分析纹理特征描述方法的性能；由于每个图像块仅与单个类别相关，实际上可确认在检索到的前N个样本中，与查询图像属于相同类别的候选图像块数量；由此，可计算与查询图像类别相关的检索样本与总检索样本之间的比率，作为检索准确度；并记录每次检索任务的执行时间以进行性能分析；

Step7.每分类准确度分析

在所有查询图像块的检索准确度检索的基础上，进一步分别针对每个类别的图像块进行计算，进行每分类准确度分析，并输出最终检索结果；

除了全局实验，针对每个类别的分类图像块还进行了实验，实验对整个数据集进行检索，每种纹理描述符取最佳参数，对不同图块大小、重叠模式和距离度量的值取各自的平均值，并计算每个类别的个体准确性。

优选的，在步骤Step2所述的地震图像输入之前，每张灰度图像要经过一定的预处理步骤(灰度级的设定，将每采样值用6bit量化，即得到64个灰度级)，将图像强度重新调整为64个灰度级，然后，进行直方图均衡化。

优选的，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用局部二值模式(LBP)对图像空间域纹理特征进行提取，其具体过程包括：

LBP用于描述一个像素与其局部圆形邻域之间的强度(灰度)差，表示为(P，R)，其中P定义了均匀分布在半径为R的圆形邻域中的像素数量；为确保对强度变化具有鲁棒性，LBP采用差分的符号(而非差分的实际值)来形成唯一性二进制编码，以描述局部纹理模式；

(1)LBP计算如下：

式中：g

(2)利用均匀描述子(Uniform Local Binary Pattern，LBP

具体为设邻域像素的个数为P，均匀LBP描述子从原始的LBP 的2

(3)进一步考虑纹理模式的旋转效应，将具有相同的循环位移码的二值模式归入一个旋转不变模式，表示为LBP

其中，

当纹理图像旋转时，灰度值g

通过统一化模式，LBP

(4)利用局部半径指数(LRI)的变体LRI-A计算给定方向上的边间距离，表示相邻的平滑区域的宽度，相应的计算规则如下，

对于像素x

式中：Th为阈值，定义一个边；

K为指定纹理元素的半径大小(单位是1个像素间隔)；

通常情况下，纹理元素的LRI探测半径最大为4个像素间隔； Th表示一个边，同时控制着噪声敏感度，而K则决定了计算的复杂度。

实施例2：与上述实施例1不同的是，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用基于信号处理的纹理提取方法进行纹理特征提取，其具体过程包括：

由Gabor函数形成的Gabor滤波器的本质，是一类特殊的加有 Gauss窗的Fourier变换，其中窗函数对于图像则决定了其在空间域中的局限特质，如果移动窗函数窗口中心，就可以获得不同的位置上图像的空间域的信息；除此之外，因为经过Gabor变换后Gauss函数的属性不变，因此不管是在空间域上还是在频率域上Gabor变换都是局部的；如此说来，如果用Gabor变换针对纹理图像来进行纹理分析，就可以既满足空间域上的局部性，又满足频率域上的局部性，用二维 Gabor函数所形成的二维Gabor滤波器来进行图像的纹理特征提取；另外，还有一个重要原因，就是Gabor函数与人类及其他的哺乳动物的视觉表皮皮层细胞的感受视野相似；再加上很好的空间局部特性和方向选择特性，可以捕捉到纹理图像在不同频率、不同方向上的边缘和局部特性；正因为如此，Gabor滤波器的理念是利用不同参数的 Gabor滤波器对图像进行卷积，并利用结果的统计量(通常为均值和方差)作为纹理描述子；

(1)设wavelength为所构造滤波器的大小，θ为滤波器的方向， f为径向频率，将Gabor滤波器分为奇滤波器和偶滤波器，离散Gabor 滤波器(sin和cos)定义如下：

其中，x

x

y

(2)设u，v表示滤波器尺度，W

实施例3：与上述实施例不同的是，步骤Step3所述的提取查询图像和数据库中的纹理特征的过程为利用基于统计的纹理结构描述方法进行纹理特征提取，其具体过程包括：

统计类方法是基于纹理图像局部区域内灰度的统计特性，测量像素的空间分布情况，从纹理图像中提取稳定的统计特征值，来描述图像中不同区域间的相异性，以及同一区域中的相似性，如：均值、中值、方差、能量、熵、惯性矩、平稳性和相关性等；基于灰度共生矩阵的描述子一直被视为纹理分析的有效工具；GLCM矩阵描述了图像中沿某个特定方向的两个相邻像素的灰度之间的共生模式；本质上，其代表了一个二值直方图，近似于相邻灰度值的联合概率分布；

GLCM能够捕捉指定方向上的选定邻域的纹理模式，举例来说，在GLCM中，当灰度共生矩阵的元素主要聚集在对角线附近时，能量值越大即图像整体的纹理较粗，反之能量越小即图像整体的纹理越细；若偏离对角线的元素的数值较大，则表示灰度水平的急剧变化；若偏离对角元素的数值较小，则表示灰度变化很小；给定一个维数M ×N的8位灰度图像I，其中I[m,n]表示位置[m,n]处的灰度值；

GLCM描述了图像中满足一定距离和一定方向的两点灰度出现的概率，每个矩阵由两个参数定义，角度φ和距离d，令C

式中：

其中，

基于GLCM，假定灰度等级的总数量为K，则相应的概率P计算如下：

式中：

P

C

实施例4：与上述实施例不同的是，为验证实施例1-3所述基于纹理特征提取的地震图像检索方法的有效性，本实施例设计以下方案对上述方法进行验证，包括

步骤一：建立实验评价采用的数据集

本方法所处理的地震图像，是在震源产生振动(或脉冲)的过程中获得的，振动(或脉冲)进入地球，穿过具有不同物理性质的地层，并返回地表作为地震数据记录；这些数据经过处理，由研究地下结构和特性的解释人员进行分析；这里使用的地震数据可以被看作是一个由水平的二维地震图像(切片)叠加而成的三维数据体(如图3所示)，这些水平切片可以沿着纵测线和横测线考虑，而垂直轴可以表示深度或时间，这取决于数据是如何处理的；本实施例的实验数据来自 Netherlands F3数据集(图例如图4)的三条纵测线，包含951×462 个像素，以及三条横测线，包含651×462个像素；对于每个实验，选择中央切片作为查询，剩下的两个切片被分配到数据库；专家先前解释了每个切片，因此将每个像素分配到8个类别，图像特征检索数据库是通过将每个地震切片划分成块来建立的；由于单个地震数据量可能达到数十GB，因此必须根据每个应用程序的需求评估在线和离线检索之间的权衡；本实施例中，我们考虑离线的情况，也就是说，对于数据库中的每一个分块，都会计算并存储纹理描述符以进行后验比较；由于数据存储和处理时间可能会决定所采用的策略(在线或离线检索)，因此必须考虑纹理描述符及其参数的选择；本实施例为纹理描述符选用的特征向量如下：

1)LBP：LBP存在多个变体，本实施例考虑了原始算子LBP和 LRI-A局部半径指数模式，两者均为每个图像分块生成256个值；

2)GLCM：特征向量的大小可计算为：特征数×距离×角度；

3)Gabor滤波器：特征向量的大小可计算为：特征数×频率×角度；

步骤二：实验方案确定

采用地震图像检索的目的是确定：

1)不同纹理描述方法对同质区域的影响；

2)不同纹理描述方法参数的选择。

地震图像检索主要包含以下步骤：

1)接收需要查询的地震图像；

2)通过提取查询图像和数据库中的纹理特征，根据相似性度量对提取出的特征向量进行比较，以计算出视觉相似度；

3)检索出与查询图像视觉相似性最高的前几幅图像。

在图像特征数据库中，将每个切片图像分割为图像块，每个图像块均有一个唯一标签表明其所属分类，该分类由地震专家(专业人员) 注释；其后，给定查询图像块B，地震图像检索的目标是返回与B有着相同地震类别的N个图像块；另外，在地震图像输入之前，每张灰度图像(tiff格式)要经过一定的预处理步骤，将图像强度重新调整为64个灰度级(即每个像素的数值范围[0，63])；然后，进行直方图均衡化；

实验环境为Intel(R)Core(TM)i5-5200UCPU@2.2GHz，系统基于Python的框架；具体的操作步骤和参数选择如下：

(1)图像块分割和数据集建立：将每张切片图分割为不同大小的图像块，并根据注释为每个图像块指定一个标签；具体来说，若某个图像块至少75％的像素属于分类A，则该图像块将被指定标签A；若某个图像块包含超过25％的噪声，则其将被丢弃；使用不同大小的图像块，以估计纹理描述子在处理不同尺寸的图像时的鲁棒性；本实施例实验中根据图像分辨率选取图像块尺寸，图5显示了其中两组样例，两个类别C

表1：两种重叠模式下的查询图像和数据集的尺寸和数量

(2)检索样本的数量：将每个查询图像块与数据库中所有图像块进行比较，以评估地震图像间的视觉相似度；然后，从数据库中返回前N个最相似的候选图像块；本实施例实验中测试了N的不同数值，即N＝5、10和15；此外，本实施例发现纹理特征描述子及其参数的行为不受N的数值的影响；因此，为简单起见，给出了N＝ 5时的实验结果；

(3)图像的相似性度量：采用欧几里得距离和曼哈顿距离，在本实施例实现的算法中，欧几里德距离和曼哈顿距离都用于比较、查找、查询图像与数据库中图像之间的相似性；

其中，欧几里得距离，也称为欧氏距离(Euclidean Distance)，在数学上称为L2范数，即向量各个元素差的平方和然后求平方根；其定义为：

式中：

d

曼哈顿距离，在数学上被称为L1范数，它计算从一个数据点到另一个数据点距离的相应样本的差异总和，其定义为：

式中：d

通过实验发现，本实施例采用欧式距离和曼哈顿距离进行相似性度量时，产生的度量结果相似，因此为了简单起见，本实施例实验数据仅展示了欧式距离的度量结果；

(4)检索准确度：通过分析描述子检索相似区域的能力，以及计算特征向量的时间，来分析纹理特征描述方法的性能；由于每个图像块仅与单个类别相关，实际上可确认在检索到的前N个样本中，与查询图像属于相同类别的候选图像块数量；由此，可计算与查询图像类别相关的检索样本与总检索样本之间的比率，作为检索准确度；并记录每次检索任务的执行时间以进行性能分析；

(5)每分类准确度：最后，本实验不但考虑了针对所有查询图像块的检索准确度和性能，而且分别针对每个类别的图像块进行了计算；该分析旨在评估描述子及其相应参数的组合，是否会对具有特定视觉特征的区域产生影响；

实验的另一项主要内容，是各个纹理特征描述方法(算子)的参数设置，我们针对每个描述符测试了不同的参数值对检索质量的影响。评价参数及相应的取值如下：

1)LBP及其改进方法LRI-A

半径：1，2，3，4，5，6；

采样点数量：8个，因为较大数值使得计算耗时较久，内存需求很大；

2)GLCM

特征：对比度、熵、均值、差分熵、差异性、互信息等；

距离：1，2；

角度：[0°，45°,90°,135°]

3)Gabor滤波器

频率：0.25，0.5，0.75和1；

特征：方差；

角度：[0°，45°,90°,135°]；

步骤三：结果与分析

通过大量实验，计算出纹理特征描述子的准确度(相对于前5个最相似的检索样本)；表2给出了Netherlands地震图像数据集纵测线的检索结果，表3给出了Netherlands地震图像数据集横测线的检索结果，两个表格结构相似，对于每个图像块尺寸、图像块重叠模式和描述子，表中每列分别是：最优参数、所得到的最高准确度和相应的检索任务所需时间(单位为秒)；表中展示的数据所用的相似度测量方法是以欧式距离来计算的，实验表明，使用欧几里德距离或曼哈顿距离获得的结果相似；

表2：Netherlands数据集的纵测线检索结果

表3：Netherlands数据集的横测线检索结果

从表2，表3中可以观察到，每个纹理特征描述子的最优参数会受到图像块尺寸大小的影响，但一般而言适用以下规则：1)在 Netherlands数据集上，对于LBP的最优半径值为1(见表2和表3)；具体来说，该数值能够分别得到91％(纵测线)和83％(横测线)的最优准确度；原LBP算子性能大多数优于LRI-A算子的性能，除了一些较大的图块，这可能由于LRI-A算子是在空间域中，使用简单的像素值比较来提取纹理特征所致；GLCM的最优特征为差异性和差分熵，差异性大约是60％的测试纵测线和33％的测试横测线中的最佳参数；而差分熵在大约30％的测试纵测线和33％的测试横测线中性能最优；Gabor滤波器的最优频率在纵测线中均为0.5，在横测线中，最佳频率值有60％选择为1；

除了全局实验，针对每个类别的分类图像块还进行了实验，实验对整个数据集进行检索，每种纹理描述符取最佳参数，对不同图块大小、重叠模式和距离度量的值取各自的平均值，并计算每个类别的个体准确性。表4表示了Netherlands数据集的检索结果，表中的每一行显示的是每个类别的检索平均准确度。

表4：每个类别检索的平均准确度

上述实验表明，在某些情况下，使用重叠模式会稍微提高精度；无论选择怎样的图像块尺寸和重叠模式，在应用最优参数时，LBP、 Gabor和GLCM均能够得到相似或相近的检索精度结果；另一方面， LBP的运行性能优于Gabor和GLCM：LBP比其他两者的速度快得多，平均比Gabor快约2.75倍，比GLCM快约4.34倍；由于LBP 的特征向量非常大(256个位置)，其要求的内存也大得多；在需要对多个地震数据体的特征进行预计算时，该内存需求可能会成为性能瓶颈；尽管如此，虽然在地震分类文献中并未进行探讨，但它仍被证明是一个可靠的描述符；对于Gabor的最优参数，实验结果表明，最优频率值可能取决于图像尺寸；对于LBP的半径参数来说同样如此，图像尺寸越大，最优数值也越高；LRI-A的性能也与图像尺寸有关，这种局部半径指数模式在图像分块较大时其准确度较高。另外，我们注意到，从表4每个类别检索的平均准确度来看，某些类别比其他类别更容易产生良好的检索结果，并且这四种描述符的性能总的来说变化趋势相同。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。