CN114859410A - 基于差分进化算法的b-cosfire滤波器的地震同相轴提取方法 - Google Patents

Abstract

本发明公开了一种基于差分进化算法的B‑COSFIRE滤波器的地震同相轴提取方法，本发明的方法中，只用两个值来描述每个像素，一个来自条形结构选择滤波器，另一个来自条形结构末端选择滤波器，并通过求和将它们合并。本发明提出的B‑COSFIRE滤波器是通用的，因为它的选择性在实现时不是预定义的，而是由自动配置过程中给定的原型模式确定的。本发明使用的加权几何平均值是一个非线性函数，它只有当所有感兴趣点的滤波器都存在时才会产生响应。这种类型的函数比其他方法对噪声更具鲁棒性。本发明通过差分进化算法可自动选择不同宽度的条形结构，具有高度灵活性，可拾取不同宽度的同相轴。本发明实现了同相轴目标轮廓的显著提取，为后续地震解释、储层预测等做准备。

Description

基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取 方法

技术领域

本发明属于地震勘探技术领域，涉及一种组合移位滤波响应模型(B-COSFIRE)滤波器，并利用差分进化算法确定滤波器子集，以最大限度地提高同相轴提取任务的性能，是一种基于视神经节细胞感受野特性(B-COSFIRE)及差分进化算法的地震数据同相轴检测方法。

技术背景

在地震勘探图像处理过程中，对地层剖面进行构造解释，并且对地震图像的检测点的位置说明是对地震剖面进行构造解释的内容。由于视觉认知过程在地震数据处理中往往遵循地震波反射剖面同相轴拾取的结构影响，因此通过地震剖面的同相轴拾取视觉信息来获取有效油气在地层层位在地震剖面中的位置，然后通过有效的地层层位信息来实现地震数据处理解释，以提高地震勘探的效率。

在地震勘探中采集信号中的信息大部分都包含在同相轴拾取当中，所以我们讨论如何实现视觉特征提取是很有意义的，如:同相轴的自激时间可以在一定程度上反映界面的深度信息，同相轴的形状反映介质中地震波的传播速度信息。由于同相轴的地震子波变化源于地震波传播过程等条件，因此，同相轴的检测和拾取对于地震数据的处理和解释非常重要。通常获取反射地震资料的方式为:地震波由人工地震源(常用爆炸物)激发。在地震波传播期间，当地震波遇到具有不同介质特性的岩石界面时，部分能量将被反射回地面。具备视觉特征的计算模型接收获得多个地震记录(信号)，每道地震信号分别对应一个检波器。地震勘探地震波在地下反射多个反射界面时，每道地震信号就会产生包含多个反射地震子波。此时，连接地震记录上振动剖面相位的极值(通常称为峰或谷)的线称为同相轴。

讨论地震数据中的同相轴信息，当前对地震勘探已开发出各种不同的检测技术，主要包括:

(1)边缘检测法：是将二维地震信号的每个采样点振幅值转换为不同的灰度值，并将地震波道集视为灰度图像，以便在图像处理中使用边缘检测方法来检测图像中的相轴。用这种方法的同相轴拾取信息比较模糊，分辨率受空间和频率干扰大，受噪声因素影响很大，检测结果是灰度突变区域的包络，分辨率很低，往往不能直接用作同相轴检测结果。

(2)神经网络法：是使用已知的同相轴作为标准样本来模拟神经网络，并使用错误返回方法逐渐修改神经元之间的连接权重。神经网络方法的过程需要足够的样本支撑，然而样本的选择更加困难，新样本的采用也会影响神经网络，大量的迭代操作导致更长的处理时间。

(3)互相关法：基本原理是利用地震道之间同相轴的波形相似度特征提取易受噪声影响的同相轴，空间分辨率随着计算迹线的增加而减小，使计算量变大。

(4)链匹配法：在这种方法是每个波形由具有多个特征的峰和谷链接表示，然后将同相轴的拾取转换为链之间的最佳匹配问题，即找到最小的链接路径。这种方法同样受噪声、空间与频率的干扰影响，所以同样难以解决复杂的同相轴处理问题。

(5)Tu等提出一种自动拾取同相轴系统：此方法包括三个步骤：二维匹配滤波器，卡尔曼滤波器和灵活模板。

视觉认知过程在地震数据处理中，根据地震波反射剖面的同相轴拾取对结构解释的影响，通过地震剖面的同相轴拾取来获取有效油气在地层层位中的位置，然后通过有效的地层层位信息来实现地震数据的处理与解释，以提高地震勘探的效率。模拟生物视觉识别机制来实现地震数据目标轮廓提取将会大大提高地震勘探与人工智能的精密链接，这正是当下研究的热点和难点。

发明内容

为了解决现有技术中存在的问题，本发明的目的是提出一种基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，用于地震资料同相轴拾取技术，拾取的同相轴可用于层序地层的研究、储层预测及储藏特征提取。

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

基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，包括以下步骤：

步骤1)：获取地震数据，并对地震数据进行预处理；

步骤2)：用DoG模型的二维表达式模拟感受野模型，即B-COSFIRE滤波器，所述DoG模型能够实现具有兴奋性的中心区域和抑制性环绕的抑制域；

步骤3)：设置两组B-COSFIRE滤波器，将两组B-COSFIRE滤波器的响应组合在一起，用于地震数据的同相轴拾取，所述两组B-COSFIRE滤波器中，一组用于检测同相轴，另一组用于检测同相轴末端；

步骤4)：以预处理后的地震数据作为B-COSFIRE滤波器的输入，对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波器的输出；

步骤5)：针对步骤4)的输出结果，用差分进化算法作为特征选择方法评估对系统性能的影响，确定最大化同相轴描绘性能的最小特征子集，完成地震同相轴提取。

进一步地，所述地震数据为sgy格式。

进一步地，所述预处理具体为：叠前地震数据的包络峰值瞬时频率与倾斜叠加峰值瞬时振幅的小波融合。

进一步地，所述DoG模型的二维表达式如下所示；

其中，具有兴奋性的中心区域的对应值设定成0.5σ，σ为四周抑制域的Gaussian函数的空间分布程度，

为On中心型感受野，即中心为兴奋域，外围为抑制域，

为Off中心型感受野，定义为：

对于给定位置(x,y)和地震数据I的强度分布I(x',y')，DoG的响应c_σ(x,y)计算如下：

其中，*为卷积，|·|⁺为半波整流。

进一步地，所述对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波器的输出，具体为：

用一个长度为ρ_i的向量将每个模糊的DoG响应向支撑区域的中心移动，相关的变化矢量是(Δx_i,Δy_i),Δx_i＝-ρ_icosφ_i，Δy_i＝-ρ_isinφ_i，在φ_i反方向上计算元组(σ_i,ρ_i,φ_i)的模糊和DoG的响应并将其进行移位计算，如下式所示：

式中，σ_i为DoG模型的标准差，(ρ_i,φ_i)为极坐标，G_σ'(x',y')为DoG最大加权响应，-3σ'≤x',y'≤3σ'，其中，σ'＝σ'₀+αρ_i，σ'₀和α的值是常数。

进一步地，所述步骤4)中求取其加权几何平均具体如下所示：

式中，

S为B-COSFIRE滤波器输出响应r_s(x,y)的对应元组，加权几何平均值是一个与型函数，即一个B-COSFIRE滤波器只有当所有传入模糊和移位响应

大于零时才能实现响应。

进一步地，所述差分进化算法以变异操作、交叉操作和选择操作进行迭代搜索。

进一步地，在所述变异操作中，利用局部适应度函数引导各子成分的变异方向。

进一步地，在所述交叉操作中，采用均匀交叉的方式。

与现有技术相比，本发明至少具有以下有益效果：

本发明提出一种视神经节细胞感受野特性(B-COSFIRE)及差分进化算法的地震数据同相轴检测方法，该方法是一种可训练的滤波器方法，采用一种无监督的方式，与基于机器学习技术的监督方法相比，本发明的方法中，我们只用两个值来描述每个像素，一个来自条形结构选择滤波器，另一个来自条形结构末端选择滤波器，并通过求和将它们合并。这提供了一种非常有效的拾取同向轴的方法。而后者随着需要执行更多的操作，时间复杂性随着维度数的增加而增加。本发明提出的B-COSFIRE滤波器是通用的，因为它的选择性在实现时不是预定义的，而是由自动配置过程中给定的原型模式确定的。本发明使用的加权几何平均值是一个非线性函数，它只有当所有感兴趣点的滤波器都存在时才会产生响应。这种类型的函数比其他方法对噪声更具鲁棒性。本发明通过差分进化算法可自动选择不同宽度的条形结构，具有高度灵活性，可拾取不同宽度的同相轴。本发明实现了同相轴目标轮廓的显著提取，为后续地震解释、储层预测等做准备。

进一步的，地震数据先进行预处理后进行B-COSFIRE滤波器处理，所述预处理为叠前地震数据的包络峰值瞬时频率与倾斜叠加峰值瞬时振幅的小波融合，能最大限度地提供最优质量数据。

进一步的，在条形结构末端配置一个新的B-COSFIRE滤波器来实现末端响应。本发明使用位于条形结构末尾的点作为兴趣点。与对称B-COSFIRE滤波器相比，非对称B-COSFIRE滤波器在条形结构末端获得了更强的响应。

进一步的，加权几何平均值是一个与型函数，即一个B-COSFIRE滤波器只有当所有传入模糊和移位响应的值大于零时才能实现响应。模糊响应和移位响应的贡献随着距离B-COSFIRE滤波器支撑中心的增加而减小。

进一步的，B-COSFIRE滤波器的多功能性在于其可训练性，它们可以根据目标为滤波器自动配置最优参数，得到最优拾取结果。

进一步的，本发明提出的B-COSFIRE滤波器是非线性的，因为它通过乘以一组高斯差分(DoG)滤波器的输出来实现方向选择性，这些滤波器的支持以共线的方式对齐。它可以容忍旋转的变化和轻微的变形。此外，COSFIRE是一种可训练的滤波器方法。这意味着滤波器的选择性在实现中没有预定义，而是在自动配置过程中由用户指定的原型模式(例如直容器、分叉或交叉点)确定的。

进一步的，采用差分进化算法来根据B-COSFIRE滤波器的响应构造特征向量，差分进化算法可评估许多特征的组合贡献，可探索更大的解空间。

进一步的，本发明地震数据是sgy格式，其包含大量信息，本发明根据视觉认知方法应用于地震领域的独特性，不是将地震数据转化为灰度图像，而是用数据本身，从数据中可以提取更多属性信息，本发明即用sgy格式的地震数据并从中提取了瞬时频率、瞬时振幅等属性后，进行同相轴拾取的操作，最大限度地提高分辨率、保留信息量。

进一步的，本发明基于叠前及叠后地震数据进行同相轴拾取。叠前地震数据比叠后地震数据含有更丰富的信息量，包括地震采集到的共炮点道集、VSP数据等；并且没有经过NMO拉伸和叠加，振幅和频率信息没有损坏，可以反映一些细微的地层特征，但是其缺点是信噪比低，噪声大。因此本发明针对叠前地震数据信噪比低的特点，提出一种高精度的同相轴拾取方法。当然，本发明的方法同样也适用于叠后地震数据的同相轴拾取。

附图说明

图1为Radon变换示意图。

图2为使用水平合成的条形结构的B-COSFIRE滤波器配置示例，其中(a)为B-COSFIRE滤波器的原理图，(b)为DoG滤波器的响应，(c)为沿着目标点的同心圆的局部最大DoG响应，(d)为最终滤波器的草图。

图3为B-COSFIRE滤波器的示意图。

图4为使用垂直合成的条形结构末端的B-COSFIRE滤波器配置示例，其中(a)为条形结构末端原型，(b)为通过(a)配置的非对称B-COSFIRE滤波器。

图5为实际地震资料，其中(a)为原始叠前CSP资料；(b)为(a)的瞬时振幅剖面。

图6为包络峰值瞬时频率。

图7为包络峰值瞬时频率与倾斜叠加峰值振幅的小波融合剖面。

图8为Radon变换得到的倾角剖面。

图9为B-COSFIRE滤波器响应结果。

图10为B-COSFIRE滤波器拾取同相轴的结果。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

针对地震剖面同相轴轮廓及位置特征的提取，本发明提出一种视神经节细胞感受野特性(B-COSFIRE)及差分进化算法的地震数据同相轴检测方法。该方法用DoG模型的二维表达式来模拟感受野模型，该模型可实现具有兴奋性(即阳性)的中心区域和抑制性环绕(即阴性)的抑制域的功能。对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波响应模型。最后用差分进化法优选最佳特征选择方法，以确定最大化同相轴描绘性能的最小特征子集。本发明提出的模型能够对数据量庞大、噪声干扰繁多的地震数据目标轮廓进行有效的同相轴边缘检测，和其他传统的边缘检测相比，本发明提出的B-COSFIRE滤波器能拾取不同宽度的地震数据的同相轴，更加符合生物视觉机理，该方法拾取的同相轴比现有方法拾取的同相轴更连续和精确，为今后地震数据解释提供新的研究方向。

连接地震记录上振动剖面相位的极值(通常称为峰或谷)的线称为同相轴。在地震勘探中采集信号中的信息大部分都包含在同相轴拾取当中，所以讨论如何实现视觉特征提取是很有意义的，如：同相轴的自激时间可以在一定程度上反映界面的深度信息，同相轴的形状反映介质中地震波的传播速度信息。由于同相轴的地震子波变化源于地震波传播过程等条件，因此，同相轴的检测和拾取对于地震数据的处理和解释非常重要。因此本发明提出一种基于反射地震数据的高精度同相轴拾取方法。

本发明提出一组B-COSFIRE滤波器来解决细长结构的拾取问题，即地震数据的同相轴拾取，该滤波器可选择不同宽度的同相轴。B-COSFIRE滤光片法最初是用于描绘视网膜血管的。这种滤波器也被用于分析计算机断层扫描血管造影(CTA)图像。这表明它们适用于各种细长结构的拾取及分割应用。本发明将两个B-COSFIRE滤波器(一个专用于检测同相轴，另一个专用于检测同相轴末端)的响应组合在一起。同相轴末端滤波器的参数选择应确保两个滤波器的性能最大化。这意味着同相轴末端滤波器的配置依赖于同相轴滤波器。滤波器的参数选择很重要，选择每个滤波器的配置参数是为了在所有条形结构的常用宽度上发挥最佳性能。本发明通过信息论和机器学习来确定B-COSFIRE滤波器的子集，对不同宽度的条形结构进行选择，用差分进化算法作为特征选择方法评估了对系统性能的影响。

首先对地震数据进行预处理，所述预处理为叠前地震数据的包络峰值瞬时频率与倾斜叠加峰值瞬时振幅的小波融合；然后介绍B-COSFIRE滤波器原理和基于差分进化算法的特征选择过程。

下面简单介绍预处理过程。在本发明的实施例中，假设震源子波用如下具有4个参数的常相位子波近似：

式中，

为调制角频率，δ为能量衰减因子，A'和φ分别为振幅和相位，t为时间，

对(1)式两边做Fourier(傅里叶)变换得：

式中，ω为角频率，φ为相位。

根据Barnes的定义，常相位子波的包络峰值瞬时频率(EPIF)为子波振幅谱对Fourier频率加权的平均频率，即

式中，τ为时间，f为频率，S(τ,f)为τ时刻震源子波的傅里叶变换。

将式(2)的模代入式(3)得到零时刻震源子波的EPIF：

倾斜叠加峰值振幅是基于瞬时振幅计算的，瞬时振幅的计算公式如下：

A(t)＝|s(t)+i·H[s(t)]| (5)

其中，

为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为一道地震数据，S(b,a)是地震数据的小波变换，

为小波函数，g(t)为Fourier变换的实部，a为尺度因子，b为平移因子，s(t)关于g(t)的小波变换定义为：

式中，t,b∈R，a>0；g(t)∈L¹(R,dt)∩L²(R,dt)，L¹表示L¹范数，即向量中各个元素绝对值之和，

是g(t)的复共轭。

倾斜叠加即局部Radon变换，本发明采用的是局部线性Radon变换。线性Radon变换的积分路径是线性的，又称为τ-p变换，其频率域的离散形式为：

式中，p为直线的斜率，M(f,p)＝∫m(τ,p)e^-j2πfτdτ，

为瞬时振幅A(t,x_m)的傅里叶变换，x_m为第m个参考道，nx为进行叠加的道数，m(τ,p)是时间域的Radon变换。局部线性Radon变换的计算过程如图1所示：在瞬时振幅剖面上选定参考道，对于参考道上的某个时间截距τ_j，将参考道附近的几道沿n_p个具有不同斜率p_j(j＝1,2,…,n_p)的直线(斜率以Δp为间隔采样)进行叠加，计算该时间截距处瞬时振幅沿不同方向叠加的和，将该叠加值记录在τ-p坐标轴相应的位置(τ_j,pj)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大。将最大叠加值的平均值放置在t-x域(时间-距离域)中的对应位置(τ_j,x_m)上，可以构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面。线性Radon变换的计算方法有很多种，本发明采用高精度频率-空间域矩阵相乘法计算。

至此已求取了包络峰值瞬时频率和倾斜叠加峰值振幅(二者均为剖面形式的二维图像)，接下来将二者进行小波融合。基于小波变换的二维图像融合主要分为三个步骤：

(1)对每张待融合的原图像分别进行小波分解，得到不同方向的分解系数。

(2)对分解得到的不同方向的系数进行融合处理。

(3)对融合后得到的系数进行小波逆变换即可得到最终的结果图像。

接下来求取B-COSFIRE滤波器。人类的视觉神经系统经过进化，可以接收和认知外界场景，具有极高的分辨率。计算机视觉即研究如何仿真实现与此类似的视觉机制。当感受器官接收到光刺激及兴奋时，神经节细胞就会把感觉信息输送给中枢神经，产生响应的区域称为经典感受野。经典感受野是同心圆感受野。在同心圆感受野外面，发现一个对经典感受野具有调制作用的大外周，则经典感受野加上大外周就是非经典感受野。神经生理学将大脑视觉系统分为三级：视网膜部分，外膝体部分，视皮层部分。感受野有视网膜神经节细胞感受野、外侧膝状体神经细胞感受野、视皮层细胞经典感受野等，其中视皮层细胞经典感受野又分为简单细胞感受野、复杂细胞感受野、超复杂细胞经典感受野等。经典感受野的结构是同心圆式的，中心与外周具有相互拮抗作用。这种结构可采用双高斯差(Differenceof Gaussian,DoG)模型来描述。图2(a)为B-COSFIRE滤波器的原理图，配置为一个水平杆，它用作输入中心DoG滤波器在其支撑区域中心的特定位置的响应，图2(b)为DoG滤波器的响应，图2(c)为沿着目标点(由中心的交叉标记标识)的同心圆的局部最大DoG响应，图2(d)为最终滤波器的草图，图中斑点的大小对应于高斯模糊函数的标准偏差。下式给出DoG模型的一维表达式：

其中，A₁,A₂分别为中央兴奋域及其周围抑制域的敏感度，σ₁为中央域Gaussian函数的均方差，σ₂为抑制域Gaussian函数的均方差。

COSFIRE为移位滤波响应模型组合，将滤波响应组合与条形选择相结合的方法叫B-COSFIRE，能实现条状特征的提取。其对方向变化敏感，对转动与微小的形状变化有较强的应激性和适应性。本发明提出的B-COSFIRE滤波器是非线性的，其非线性响应通过乘以一组高斯差分(DoG)滤波器的输出来实现方向选择性，这些滤波器组支持以共线的方式对齐，它可以容忍旋转的变化和轻微的变形。下式给出了DoG模型的二维表达式来模拟感受野模型，该模型可实现具有兴奋性(即阳性)的中心区域和抑制性环绕(即阴性)的抑制域的功能：

其中，中心兴奋域的对应值设定成0.5σ，σ为四周抑制域的Gaussian函数的空间分布程度，即标准差，这种类型的功能是对大脑外侧膝状核(LGN)中的一些细胞的一种公认的计算模型。

为On中心型感受野，即中心为兴奋域，外围为抑制域。

为Off中心型感受野，定义为：

上述DoG模型能够模拟细胞的同心圆拮抗结构。

对于给定位置(x,y)和地震数据I的强度分布I(x',y')，DoG的响应c_σ(x,y)可计算如下：

其中，*为卷积，|·|⁺为半波整流。然后，考虑沿着同心圆在给定的兴趣点周围的DoG响应，并从它们中选择具有局部最大值的DoG。

B-COSFIRE的方向选择与条型组织的方向相同。我们用三个参数来描述第i个点：利用设定好σ的DoG模型，以指定点作为圆心，在N个同心圆范围内分析c_σ(x,y)响应。针对地震数据同相轴双曲线的特点，B-COSFIRE滤波组合的设置如图3。浅灰圆形的圆心为B-COSFIRE滤波组合的中心点，中间竖线与同心圆相交产生的黑点为DoG响应的极值点。该滤波组合参数为：

S＝{(σ_i,ρ_i,φ_i)|i＝1,…,n} (12)

S定义了一个对给定原型具有选择性偏好的B-COSFIRE滤波器。其中，σ_i为DoG滤波器的标准差，(ρ_i,φ_i)为极坐标，表示其位置，n是需要配置的DoG滤波器的数目。下式设定的参数值用来对图2(c)的横线式进行分析，由图2(c)所示的输入模式的自动分析确定，其中ρ和φ由前面得到的点自动分析生成：

首先使用DoG函数卷积输入图像，该DoG函数在集合S的元组中具有指定的标准偏差。接下来进行DoG滤波器模糊移位运算，以便在相关点的优先位置允许一些误差。模糊操作定义为计算一个DoG滤波器的加权阈值响应的最大值。模糊操作计算DoG最大加权响应G_σ'(x',y')的局部邻域中的最大DoG响应以提高其容错性，权值利用DoG滤波响应和高斯函数c_σ'(x',y')的系数相乘取得，其标准差σ'是距离滤波器支撑中心ρ_i的线性函数：σ'＝σ'₀+αρ_i，σ'₀和α的值是常数，我们根据应用程序调整它们，ρ_i表示相关点的优先位置和滤波器支撑中心的距离。

然后，本发明用一个长度为ρ_i的向量将每个模糊的DoG响应向支撑区域的中心移动，这是φ_i的互补角。相关的变化矢量是(Δx_i,Δy_i),Δx_i＝-ρ_icosφ_i，Δy_i＝-ρ_isinφ_i。在φ_i反方向上计算元组(σ_i,ρ_i,φ_i)的模糊和DoG的响应并将其进行移位计算，如下式所示：

式中，-3σ'≤x',y'≤3σ'。

接下来计算B-COSFIRE滤波响应模型。B-COSFIRE滤波器输出响应r_s(x,y)为对应元组S中所有模糊和移位的DoG响应的加权几何平均：

式中，

(ρ_i,φ_i)是第i个DoG响应相对于滤波器支撑中心的极坐标，|·|_τ表示最大响应在分数τ(0≤τ≤1)处的响应的阈值。加权几何平均值是一个与型函数，即一个B-COSFIRE滤波器只有当所有传入模糊和移位响应

大于零时才能实现响应。

上述程序配置了B-COSFIRE滤波器，该滤波器可用于水平定向容器。为了实现多方向选择性，可以通过使用不同方向的原型模式来配置多个B-COSFIRE滤波器。加权几何平均为“与函数”，在全部响应均满足“＞0”时，可以得到响应值，指定点离中心点越远，模糊及移位的作用将越小。由此可见，设定好的条形方向是该滤波器的前提，我们操纵每个元组的参数ψ，并创建一个新的集合R_ψ(S)，这里对集合S进行整合，加入R_ψ(S)实现具备旋转不变形，则有：

它表示一个B-COSFIRE滤波器，其方向偏好为ψ弧度，与原始滤波器S的方向偏好相偏移。可跟踪各方向的条型目标，将所得全部响应求和，计算各像素在(0°，15°，30°，…，165°)方位的响应值，选择最大值作为该点特征。本发明选取一组具有不同方向偏好的B-COSFIRE滤波器，可表示为：

其中，

n_r＝12。

因为地震数据一次反射波及多次反射波不连续、对称特征不明显，同相轴拾取时难以完整检测，而B-COSFIRE为“逻辑与”函数，仅在全部输入均被激活时才会得到响应，原则上其无法在不连续处得到响应。因此，考虑在一次及多次反射波处构建一个新的滤波器，将之前检测连续同相轴及直达波的滤波器称为对称结构的B-COSFIRE，将不连续处及反射波处的滤波器称为非对称结构的B-COSFIRE。该滤波器将在360°的范围内计算响应，也就是说，该滤波器和对称滤波器不一样的是，它仅针对一侧的目标点，实际应用时要把方位数目加倍，得到所有方向的条型结构端点，各个方向均存在是目标点的可能，如图4所示，同心圆的圆心为其响应中心，黑点为响应极值。

地震同相轴的宽度可能从1个像素到若干像素不等，这取决于炮检距，炮检距越大，高频损失越大，低频越多，波形越宽。因此，本发明配置了一组B-COSFIRE滤波器，由21个条形结构{S₁,…,S₂₁}和21个条形结构末端{S₂₂,…,S₄₂}组成，可选择不同厚度的条形结构。大尺度滤波器对宽同相轴具有选择性，并且对背景噪声具有耐受性，但沿窄同相轴的响应小。相反，小尺度滤波器对窄同相轴具有更高的选择性，但对背景噪声的耐受性较差。他们的响应融合在一起，有望在各种尺度上实现更好的同相轴提取性能。我们为每个同相轴位置(x,y)构造一个像素级特征向量v(x,y)，该向量由滤波器组中42个BCOSFIRE滤波器的响应，加上瞬时振幅的强度值g(x,y)组成：

其中，

是BCOSFIRE滤波器S_i的旋转容限响应。

在分类之前，本发明将反双曲正弦变换函数应用于特征向量的每个元素。它减少了数据中的偏斜，定义如下：

对于υ_i和θ>0的较大的值，函数的应用类似于对数变换。令θ→0，f(υ_i,θ)→υ_i。然后计算Z分数，使43个特征值都标准化。本发明分别对每个图像应用Z分数归一化过程，以补偿图像之间的照明变化。

以上设计的滤波器库可能有许多冗余滤波器。本发明研究了各种特征选择方法，以确定最大化同相轴描绘性能的最小特征子集。本发明使用由大小为N×43的矩阵组成的训练数据作为输入，其中N对应于从训练数据中随机选择的像素数，列数对应于滤波器组的大小加上EPIF数据。

反演优化技术已广泛应用于许多领域，在地震勘探中可用于反演地下复杂构造的介质参数。反演方法有很多，如波形反演、层析成像等，波形反演优化技术又有局部优化类反演方法和全局优化类反演方法等。局部优化方法计算速度快，不足是高维参数反演中矩阵求逆困难，容易陷入局部极小值，且反演精度依赖于初值。与局部优化类方法相比，全局优化类方法对初值无依赖性，且不需要计算矩阵的逆，不足是运算量大，计算耗时。近年来由于计算机运算能力显著提高，全局优化方法已得到广泛应用。

全局优化算法有模拟退火法、差分进化算法、遗传算法、协同进化算法，以及改进的协同选择差分进化算法(DE-CCS)、协同变异差分进化算法(DE-CCM)等。进化算法是根据自然界中各物种进化的规律，依据物竞天择、适者生存的原则，利用计算机模拟个体(即物种)达到最优过程的算法，通过群体一代代的迭代，将不利于物种延续的基因不断淘汰，将好的基因传承下去，最终使个体达到最优。本发明使用差分进化算法(DE)，用局部适应度函数引导每一个子成分最优基因的选择，然后进行反演，在大量可能的组合中搜索性能最好的特征子集。

DE方法首先给定一个含M个个体的初始群体X：

X＝[x₁,x₂,…,x_M] (20)

该群体中的第k个个体表示为：

x_k＝[x_k,1,x_k,2,…,x_k,N],k＝1,…,M (21)

式中，M为初始个体数，根据实际情况设置，如M＝400，N为43个待估参数的维度(即EPIF和42个B-COSFIRE滤波器)。适应度函数计算训练数据与选定列的平均精度。用L₁范数定义的第k个个体的局部适应度函数为：

式中，d_o,j(·)为第j列选定的同相轴数据，d_k,j(·)为由第k个个体正演计算的第j列训练数据，NR为观测地震数据的道数，g(l)是高斯窗函数，定义为：

式中，l为变量，σ为标准差，用来控制窗的宽度，高斯窗的中心位于子波包络峰值处。对于叠后数据，时窗的中心应位于子成分的中心，对于叠前数据，随着偏移距的增加，时窗的中心应沿着双曲线滑行。用第k个个体x_k为B-COSFIRE的参数并拾取地震同相轴的位置，每道计算一个局部适应度值，则第k个个体对应的局部适应度值矢量为：

f_k＝[f_k,1,f_k,2,…,f_k,N] (24)

式中，元素f_k,1,f_k,2,…,f_k,N为1～N道的局部适应度值。

因此对于含M个个体的群体X来说，其正演数据对应的局部适应度值为M×N的矩阵：

初始群体X也是M×N的矩阵。

差分进化算法以变异操作、交叉操作和选择操作进行迭代搜索，三个主要操作分别介绍如下。

有效的变异策略可以使个体中的基因(待估参数)快速得到改良，经过多代循环最终得到最优个体。为了寻找基因变异的正确方向，在变异操作中，本发明利用局部适应度函数来引导各子成分的变异方向，从而提出一种新的变异算子。本发明将群体X中第j列的所有变量进行升序排列：

式中，上标s表示排序后的序列。该序列对应的局部适应度值为：

令排序前的序列对应的局部适应度函数为f_j，

由

中的元素对应于f_j中的值所构成，而非将f_j中的元素进行升序排序。从

中拾取该矢量的最小值以及最小值两侧邻域内的NP个值，这(2×NP+1)个值的解析表达式可用具有最小值的凹函数多项式拟合来逼近：

式中，A，B，r₁，r₂为复合方程的4个参数，可通过反演方法求取，x为第j列排序后的变量，y为对应的局部适应度值。该凹函数的极小值点对应的变量可能是一个更好的变异基因。为了确定A，B，r₁，r₂这四个参数的值，定义目标函数为：

式中，

为计算的局部适应度值，

为待拟合的多项式的值。因为有4个待估参数，Hessian矩阵的逆比较容易求取，因此可以用高斯-牛顿法求解该方程组。

然而，全局优化算法需要大量迭代，其计算量本身已比较大，而高斯-牛顿法是一种反演方法，计算速度比较慢。为了提高计算速度，实际中本发明采用如下的直接方法进行变异。令局部适应度函数的最小值对应的变量为x_r,j，在群体的第j列中随机选取除最小值外的两个变量x_r1,j和x_r2,j，则第k个个体的第j个元素的变异策略是：

式中，F是控制参数，其值是0～1之间的随机变量。所有个体进行循环以更新第j列的基因，变异后的基因在最优值点附近扰动，不同个体变异后的值略有差异，这是因为随机选取的变量不同以及控制参数不同，而这也符合自然规律，体现了变异发生的随机性。

将所有个体进行循环后，为了使深层参数反演精确，本发明构造局部循环流程，即进行多次个体循环。该局部迭代循环过程中，我们观察到一种现象，随着迭代次数的增加，局部适应度值(即误差能量)依然很大时不同局部适应度值之间的变化却很小。这是因为随着变异次数的增加，这一列参数逐渐接近最优值，在参数未到达最优值时他们之间的差异已越来越小，变异缺乏突破的动力，导致变异前后的值差异很小，如果不引入新的变异机制，群体中第j列的所有分量将趋于常数。为了改善反演精度，本发明提出如下的加权变异策略。给定误差能量的精度阈值，当误差能量大于阈值而两次个体迭代的误差能量之差小于另一阈值时，用局部适应度值加权的方法对变异扰动量进行增益，公式如下：

式中，λ是控制参数，可根据精度选取，上标M表示变异后的值，权系数p₁为归一化的局部适应度值，计算方法如下：

式中，f_r,j,f_r1,j,f_r2,j是对应于基因x_r,j,x_r1,j,x_r2,j的局部适应度值。该加权变异策略利用了误差能量的信息，使得变异扰动量可根据误差能量自适应地进行调整。

交叉操作：随机地从目标个体

与其相应的变异个体

中抽取基因，组合成一个试验个体

常用的两种交叉方法是均匀交叉和指数交叉，本发明采用均匀交叉的方式，公式如下：

其中CR∈[0,1)是交叉率，n_rand是1到D之间的一随机整数，它保证了试验个体中至少有一个基因是来自于变异个体。

常规选择算子是根据目标个体

和试验个体

的全局适应度值将其中一个个体整体选入下一代。数学表达式为：

在每一次迭代中，我们使用90％的训练集配置一个带有线性核的SVM分类器，并将其应用于剩余的10％。在训练集中全部样本训练一次之后，我们将个体按适应度大小降序排列，只保留居前10％的个体。我们使用这组精英个体，通过对随机选择的精英个体对进行交叉操作，生成其余90％的新个体。新个体中的每一个体都有10％的概率发生变异(即将编码从1变为0或从0变为1)。我们运行这些迭代步骤，直到所有个体达到最优，不再发生变异。最后，我们选择群体中具有最高适应度值的个体的位置所对应的滤波器。

我们使用选定的特征来训练具有线性核的SVM分类器。SVM分类器特别适合于二元分类问题，因为它找到了一个最佳分离超平面，最大化了类之间的边界。

算例分析

本发明将提出的方法用于实际共炮点道集(CSP)，该道集有595道，最小炮检距为90m，检波器之间的距离为10m，1250个采样点，2ms采样。图5(a)是原始叠前地震数据，从图中可以看出，原始数据有明显的同相轴，但是背景噪声同样很大，尤其是近偏移距和浅地表处的同相轴无法分辨，且同相轴的波峰、波谷相互影响，不利于拾取倾角等其他属性信息。图5(b)是瞬时振幅剖面。图6是包络峰值瞬时频率(EPIF)剖面，从图中可以看出，包络峰值瞬时频率减少了波峰、波谷的相互干扰，信息进一步集中，但是深层的噪声依然很大。图7为EPIF与倾斜叠加峰值振幅剖面的小波融合结果，可以看出融合后的图像突显了同相轴，一定程度上弱化了噪声的影响。图8为Radon变换拾取的倾角剖面。图9为B-COSFIRE滤波器的响应。图10为用本发明的方法进行同相轴拾取的结果，可以看出，浅层的同相轴拾取较完整，深层由于噪声较大，信息淹没在噪声中，拾取同相轴较少，但厚地层的同相轴已经完整拾取，说明了本发明的有效性。

本发明用DoG模型的二维表达式来模拟感受野模型，该模型可实现具有兴奋性(即阳性)的中心区域和抑制性环绕(即阴性)的抑制域的功能。对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波响应模型。最后用差分进化法优选最佳特征选择方法，以确定最大化同相轴描绘性能的最小特征子集。本发明在分析生物视觉初级视皮层信息处理机制的基础上，构建了地震资料同相轴拾取方法。本发明提出的模型能够对数据量庞大、噪声干扰繁多的地震数据目标轮廓进行有效的同相轴边缘检测，和其他传统的边缘检测相比，本发明提出的B-COSFIRE滤波器能拾取不同宽度的地震数据的同相轴，更加符合生物视觉机理，该方法拾取的同相轴比现有方法拾取的同相轴更连续和精确，为今后地震数据解释提供新的研究方向。

Claims (9)

1.基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，包括以下步骤：

步骤1)：获取地震数据，并对地震数据进行预处理；

步骤2)：用DoG模型的二维表达式模拟感受野模型，即B-COSFIRE滤波器，所述DoG模型能够实现具有兴奋性的中心区域和抑制性环绕的抑制域；

步骤3)：设置两组B-COSFIRE滤波器，将两组B-COSFIRE滤波器的响应组合在一起，用于地震数据的同相轴拾取，所述两组B-COSFIRE滤波器中，一组用于检测同相轴，另一组用于检测同相轴末端；

步骤4)：以预处理后的地震数据作为B-COSFIRE滤波器的输入，对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波器的输出；

步骤5)：针对步骤4)的输出结果，用差分进化算法作为特征选择方法评估对系统性能的影响，确定最大化同相轴描绘性能的最小特征子集，完成地震同相轴提取。

2.根据权利要求1所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述地震数据为sgy格式。

3.根据权利要求1所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述预处理具体为：叠前地震数据的包络峰值瞬时频率与倾斜叠加峰值瞬时振幅的小波融合。

4.根据权利要求1所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述DoG模型的二维表达式如下所示；

其中，具有兴奋性的中心区域的对应值设定成0.5σ，σ为四周抑制域的Gaussian函数的空间分布程度，

为On中心型感受野，即中心为兴奋域，外围为抑制域，

为Off中心型感受野，定义为：

对于给定位置(x,y)和地震数据I的强度分布I(x',y')，DoG的响应c_σ(x,y)计算如下：

其中，*为卷积，|·|⁺为半波整流。

5.根据权利要求4所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述对DoG响应进行模糊和移位，求取其加权几何平均得到B-COSFIRE滤波器的输出，具体为：

用一个长度为ρ_i的向量将每个模糊的DoG响应向支撑区域的中心移动，相关的变化矢量是(Δx_i,Δy_i),Δx_i＝-ρ_icosφ_i，Δy_i＝-ρ_isinφ_i，在φ_i反方向上计算元组(σ_i,ρ_i,φ_i)的模糊和DoG的响应并将其进行移位计算，如下式所示：

式中，σ_i为DoG模型的标准差，(ρ_i,φ_i)为极坐标，G_σ'(x',y')为DoG最大加权响应，-3σ'≤x',y'≤3σ'，其中，σ'＝σ'₀+αρ_i，σ'₀和α的值是常数。

6.根据权利要求5所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述步骤4)中求取其加权几何平均具体如下所示：

式中，

S为B-COSFIRE滤波器输出响应r_s(x,y)的对应元组，加权几何平均值是一个与型函数，即一个B-COSFIRE滤波器只有当所有传入模糊和移位响应

大于零时才能实现响应。

7.根据权利要求1所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，所述差分进化算法以变异操作、交叉操作和选择操作进行迭代搜索。

8.根据权利要求7所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，在所述变异操作中，利用局部适应度函数引导各子成分的变异方向。

9.根据权利要求7所述的基于差分进化算法的B-COSFIRE滤波器的地震同相轴提取方法，其特征在于，在所述交叉操作中，采用均匀交叉的方式。

Images