1 引言

在地震数据中，断层主要表现为反射层位的位移或中断、横向上振幅不连续，断层两侧的反射振幅存在明显差异，这些特征在地震剖面上可以被肉眼识别，这也是人工识别断层的基础(严哲，2010)．但是人工识别断层效率低，主观性强，解释能力局限．

所以，学者们针对断裂系统的精确解释，提出并应用了很多断层的自动或半自动检测技术．Bahorich和Farmer(1995)提出了在三维断层解释技术中占据重要地位的基于互相关的第一代相干算法；Marfurt等(1998)提出了基于多道相似性的第二代相干算法；Gersztenkorn和Marfurt(1999)提出了基于特征结构的第三代相干算法．相干体技术在很大程度上提高了断层的解释效率，大部分断层自动检测技术也主要是以这三代相干体技术为核心所进行的改进．然而，传统的基于滑动窗口的相干体属性对处理数据信息中的局部非平稳特征有一定的难度(刘洋等，2014)．为了解决以上问题，学者们陆续开发了一些新的断层自动检测技术，其中Randen等(2001)提出了“人工蚂蚁”断层自动检测技术，并由Pedersen等(2002, 2003)进行了更深层次的分析阐述，该技术利用人工蚂蚁对三维地震数据体生成的突出振幅空间不连续性的断层增强属性体(倾角方位属性、混沌属性、方差属性等)进行处理，从而提取断层信息，该技术在2002年以Ant tracking模块被引入到Petrel软件中，得到广泛商用．Dorn等(2005)提出基于相干体数据的断层自动提取技术，该技术融合了地质先验知识的信号处理技术，后被应用到商业软件Epos中，但该技术只是对相干体断层响应进行了加强，对断层响应的连续性改善作用不大．赵伟(2009)利用蚁群算法对相干切片图像进行聚类处理，实现了“简单”断层线的提取和重构，但结果依赖相干体的断层检测能力，且计算效率低．李香臣(2010)利用蚁群优化聚类算法和基于图搜索式蚁群优化算法对构造特征进行了增强．严哲(2010)提出方向约束的蚁群追踪算法，进一步改进之后，提出了基于蚁群算法的自动断层追踪的方法(Yan et al., 2013)，该方法是将相干体数据作为断层增强属性体，再利用蚁群算法对相干体数据进行处理以达到检测断层的目的，与Petrel软件中Ant tracking模块的断层追踪结果对比，该方法检测到的断层连续性高，但结果好坏也直接依赖于相干体算法对断层的检测能力．

断层增强属性体会沿着断层面产生响应，这是断层自动检测的基础，但也夹杂了噪声、地层残余响应、角度不整合接触等假象．因此，为了提高断层检测方法的信噪比，本文引入了一致性的概念．一致性由Kass和Witkin(1987)提出，是指在对指纹图像边缘检测处理中，对指示局部方向场信息的各向异性强度的一个度量，其计算过程不涉及滑动窗口，计算结果受人为选取参数约束小，且不受类似地震同相轴的指纹脊谷线形状的影响，即若将其引用到地震数据处理，其计算结果将不会受地层倾角的影响．因此，本文将其作为一种新的断层增强属性应用于地震数据断层检测，并提出加权一致性的概念，来提高一致性属性的抗噪性，再结合蚁群算法构成了新的断层检测方法．在理论模型和实际数据试算阶段，对于相同的断层检测过程，本文所提方法在抗噪性、收敛性及压制地层残余响应等方面均明显优于常规的基于方差属性的蚁群算法和基于C3相干属性的蚁群算法．

2 方法原理 2.1 一致性

在指纹图像中，灰度沿指纹脊谷线垂向方向变化最快，沿走向方向变化最慢，这种特征被称为各向异性(Kass and Witkin, 1987)．该特征可用于求取指纹脊谷线的局部方向场信息，Kass和Witkin(1987)将衡量指示局部方向场的各向异性强度的度量称为“一致性”．

地震数据实际是对地下沉积地层的反射波成像，地震记录中的同相轴类似于指纹图像中的脊谷线，因此借鉴一致性算法，对地震数据中指示局部地层走向的各向异性强度作以衡量．一致性的求取过程如下：

在直角坐标系下，二维地震数据体中每个数据点的梯度向量[G_x(x，y)，G_y(x，y)]^T可表示为：

(1)

其中，T表示向量转置，I(x，y)表示点(x，y)的振幅值矩阵．[G_x(x，y)，G_y(x，y) ]^T简记为[G_x，G_y]^T，本文选用Soble算子计算梯度．

在极坐标系下，公式(1)可表示为：

(2)

其中，G_ρ为极径，表示梯度向量的模长，G_φ为极角，表示梯度向量的方位角．则直角坐标系下的梯度向量又可表示为：

(3)

定义数据点的平方梯度向量(Kass and Witkin, 1987) [G_{s, x}，G_{s, y} ]^T表达式为：

(4)

然后，提取一个以点(x，y)为中心的包含n×n个数据点的窗口B，将该数据块的平方梯度向量表示为：

(5)

其中，

类似指纹图像，地震记录中每个数据点的梯度向量总是指向灰度变化最快的方向，并与地震同相轴走向垂直．同一地震同相轴具有两侧边界线，分布在两侧边界线上数据点的梯度向量呈对称分布，即指向相反．平方梯度向量方向角是梯度向量方向角的两倍，这样指向相反的向量变成了指向相同的向量．所以对于同一条同相轴，两侧边界上点的平方梯度向量指向一致．

为了衡量数据块B的各向异性强度，即数据块内各点的平方梯度向量方向的集中程度，Kass和Witkin(1987)提出一致性概念：

(6)

定义Coh_B为数据块B中心点(x，y)的一致性，其中，｜｜表示向量的模．

如果数据块B内各点的平方梯度向量指向都相同，则该局部范围内各向异性指示性很强，与平方梯度向量方向垂直的地震同相轴走向也一致，表明地层连续，此时，所有向量模的和等于所有向量和的模，Coh_B等于1；如果数据块B内各点的平方梯度向量指向在(0-2π)范围内均匀分布，则该局部范围内各向异性特征极不明显，很难获取地震同相轴走向信息，表明地层不连续或存在强随机噪声，此时，所有向量模的和远大于所有向量和的模，Coh_B趋向于0；其他的情况，Coh_B介于0和1之间．本文给出的图件中一致性均用Coh表示．

对图 1a无噪声理论合成地震记录计算一致性，结果如图 1c所示，可以看出，在地层连续处，一致性数值较高，而在断层处及地层不整合接触位置，一致性数值较低，这样就突出了地层的不连续性，该结果表明一致性数据可以作为新的断层增强属性体应用于断层检测.对加入强随机噪声的合成地震记录(图 1b)计算一致性，结果如图 1d所示，可以看出，一致性数据在强噪声干扰下，对地层的横向不连续性识别效果非常差．所以，为了提高抗噪性，本文提出了加权一致性的概念．

2.2 加权一致性

首先将待处理的地震数据分割成N个基础数据块Γ(图 2中对应Ⅰ | Ⅱ | Ⅲ | Ⅳ | Ⅴ | Ⅵ | Ⅶ | Ⅷ | Ⅸ)，每个小数据块由n×n个数据点组成，然后再将相邻的2×2个基础数据块Γ组成新的数据块D，同时将相邻的3×3个基础数据块Γ组成一个新的数据块M，如图 2所示．

接下来计算目标数据块V的一致性．从图 2可以看出，目标数据块V在相邻的四个邻域D中都有被包含，只是在不同的数据块D中所处的位置不同，分别处在四个不同方向，每个基础数据块Γ的平方梯度向量由公式(5)决定，一致性由公式(6)决定．则数据块D_j的平方梯度向量为：

(7)

其中i是数据块D_j中的数据点．因为每个数据块D_j由四个基础数据块Γ组成，所以数据块D_j的平方梯度向量相当于对其包含的四个基础数据块Γ的平方梯度向量求和．例如D₁={Ⅰ，Ⅱ，Ⅳ，Ⅴ}，则有：

(8)

(9)

理想情况下，在数据块M内，如果反映局部地层信息的同相轴连续且相互平行，那么相互重叠的四个数据块D₁，D₂，D₃，D₄就具有非常高的相关性．此时，在这四个数据块内，所有数据点的平方梯度向量指向一致，那么由公式(7)可知，这四个数据块的平方梯度向量指向也一致，衡量数据块各向异性强度的度量值Coh_Dj均相等，并且与目标数据块V的一致性Coh_V相等．在数字图像处理中，利用邻域信息进行冗余计算可压制噪声干扰，并能有效恢复噪声图像(Wang et al., 2005)．类似的，如果目标数据块V内噪声干扰大，那么数据块V内各数据点的平方梯度向量方向就比较散乱，Coh_V值相应会变小，不能指示出真正的地层走向信息．此时，就可利用数据块D₁，D₂，D₃，D₄的平方梯度向量信息来估算目标数据块V的一致性，从而压制噪声干扰，提高信噪比．

将数据块D₁，D₂，D₃，D₄组成新的数据块D^*，那么数据块D^*的平方梯度向量为：

(10)

再结合公式(8)、(9)可知，公式(10)可用基础数据块的平方梯度向量表示：

(11)

其中，M={Ⅰ，Ⅱ，Ⅲ，Ⅳ，Ⅴ，Ⅵ，Ⅶ，Ⅷ，Ⅸ}，r_Γ表示基础数据块Γ参与计算的次数，显而易见，其值等于该数据块Γ被数据块D_j包含的次数．将系数r_Γ进行归一化：

(12)

其中．定义ω_Γ为覆盖率(Wang et al., 2005)，令ω_Γ替换公式(11)中的r_Γ，则有：

(13)

可见，目标数据块V的覆盖率是最高的．覆盖率不同表明每个基础数据块的重要程度不同．将公式(13)定义为数据块M的加权平方梯度向量．如果任意数据块Γ∈ M的平方梯度向量的方向和大小都相同，即

(14)

则有：

(15)

由公式(15)可知，当任意数据块Γ∈M的平方梯度向量的方向和大小都相同时，目标数据块V的平方梯度向量等于数据块M的加权平方梯度向量．因此，目标数据块V的一致性Coh_V就可根据基础数据块Γ的平方梯度向量加权而得：

(16)

当公式(14)成立时，有：

(17)

定义Cohw_V为目标数据块V中心点的加权一致性．

公式(17)表明：如果M内各基础数据块Γ的平方梯度向量的方向和大小都相同，Cohw_V等于1，表明该局部范围内地震同相轴走向一致，说明地层连续．如果M内各基础数据块Γ的平方梯度向量的方向和大小不尽相同，Cohw_V介于0和1之间． Cohw_V接近于0时，表明该局部范围内地震同相轴信息很弱，可能是强噪声的影响，也可能存在不连续地层．本文给出的图件中加权一致性均用Cohw表示．

图 3a和图 3b分别是对图 1a和图 1b计算得到的加权一致性数据．可见，加权一致性数据相较于一致性数据，不仅能在一定程度上压制噪声干扰，还能弱化地层不整合接触和地层残余响应所产生的伪断层信息，使地层的横向不连续性更为突出、准确．

为了验证加权一致性数据的优越性，对实际地震数据(图 4a)计算一致性(图 4b)和加权一致性(图 4c)．由图可见，对于实际地震数据，加权一致性数据相较于一致性数据能更好地压制随机噪声和地层残余响应，突出地层的横向不连续性，提高信噪比．但是作为新的断层增强属性体，无论是一致性数据还是加权一致性数据，都是基于一定窗口范围内数据块的计算，对不连续性存在着水平方向或垂直方向的弥散作用(严哲，2010)，使识别到的地层不连续性表现为一定宽度且不连续的条带，而且加权一致性并不能彻底地消除噪声干扰及地层残余响应．所以，为了进一步提高断层检测的质量，本文将加权一致性属性与蚁群追踪算法相结合，提出了新的断层检测方法——基于加权一致性的蚁群算法．

2.3 基于加权一致性的蚁群算法在断层检测中的应用

基于加权一致性的蚁群追踪算法进行断层检测的过程主要分为以下几步：

(1)  初始化蚂蚁位置

将二维加权一致性数据分为若干个互不相交的m×m的相邻数据块，在数据块中，蚂蚁的初始位置由第i个数据点的被选择概率P_i决定：

(18)

其中，Cohw_i表示数据块内第i个数据点的加权一致性，加权一致性值越低的点被选中的概率越大，可降低蚂蚁在其他区进行无意义追踪的概率．

(2)  断层追踪

蚂蚁所在点(i, j)到其搜索范围R内候选点(r, s)的转移概率为(段海滨，2005)：

(19)

其中，τ_(r，s)^n-1为第n-1次迭代后点(r，s)的信息素浓度，η_(r，s)表示点(r，s)处的启发函数，α、β分别是信息素和启发函数重要程度因子，ω表示选择概率权重(刘闻，2013)．初始时刻，各数据点信息素浓度相同．启发函数η_(r，s)由加权一致性Cohw_(r，s)决定：

(20)

可见，在初始信息素相同的情况下，加权一致性值越低的点被蚂蚁选中的概率越大，可使蚂蚁的搜索路径尽快收敛到断层上．

(3)  判断终止条件

蚂蚁追踪的终止条件是：

(21)

其中，终止条件S(0＜S≤1)表示蚂蚁追踪过程中所被记录的异常步数与总步数的比值(Yan et al., 2013)．S_min、S_max分别是给定的最小和最大终止条件(Yan et al., 2013)．

(4)  更新信息素

若当前迭代中所有蚂蚁的追踪被终止，即按(22)式对数据点(i，j)处的信息素浓度进行实时更新(段海滨，2005)：

(22)

其中，ρ(0≤ρ≤1)表示信息素挥发系数，Δ τ_(i，j)^n-1表示在第n-1次迭代中数据点(i，j)上的信息素增量(段海滨，2005)：

(23)

其中，m_ant是本次迭代中经过数据点(i，j)的蚂蚁数量，Δ τ_(i，j)^n-1(k)表示第k只蚂蚁在第n-1次迭代中在(i，j)处释放的信息素浓度，其取值与蚂蚁追踪的路径长度L(k)成正比(Yan et al., 2013)：

(24)

c为信息素更新系数．

当前迭代结束时，记录所有蚂蚁的当前所在位置，使其作为下一次迭代时的初始位置．

(5)  输出结果

不断重复(2)—(4)步骤，直到所有迭代次数完成．常规蚁群算法是直接输出信息素剖面图，但因蚂蚁追踪的随机性，致使不同参数设置下所输出的信息素浓度范围不同，易导致色标不一致．因此，本文对追踪结果做归一化改进后再输出：

(25)

其中，τ ⁿ表示第n次迭代之后的信息素浓度矩阵，C表示利用蚁群算法进行断层追踪的最终输出结果．

基于加权一致性的蚁群算法进行断层检测的流程如图 5所示．其中k_max表示进行断层追踪的蚂蚁总数，n_max表示最大迭代次数．

3 理论模型分析

为了验证基于加权一致性的蚁群算法进行断层检测的可靠性，选取三维“qdome”模型(Claerbout and Fomel, 2012)中的主测线1.5 km位置的联络测线垂直剖面进行测试(图 1a)，该模型包含顶部的水平地层、中部的高斯型地层、底部的倾斜地层、不同类型地层区域边界的水平不整合面及一个正断层．图 1b是加入正态分布随机噪声以后的合成地震记录．图 1a和图 1b的地震数据的时间深度为500个采样点，采样间隔为0.01 s，测线长度为150道地震道，道间距为0.01 km．

用基于加权一致性的蚁群算法分别对图 1a和图 1b进行断层检测．参数设置为：分割地震数据的小数据块大小为2×2；终止条件S_min=0.2、S_max=0.5；信息素浓度初始值为1，信息素更新系数c=1；信息素和启发函数的重要程度因子分别是α=1、β=3；信息素挥发系数ρ=0.01；迭代次数为5次．检测结果如图 6a和图 6b所示，可见，无论是理论无噪声地震数据(图 1a)，还是加入强随机噪声的地震数据(图 1b)，本文所提方法都能准确、清晰且完整地检测出断层信息，抗噪性强．

选取基于方差属性的蚁群算法(图 6c、图 6d)和基于C3相干属性的蚁群算法(图 6e、图 6f)对图 1a和图 1b进行断层检测．对比可见，对于理论无噪声地震数据(图 1a)，三种方法都能准确、清晰、完整地检测出断层位置，但基于C3相干属性的蚁群算法对大倾角地层很敏感，检测结果伴有明显的地层残余响应．而对于加入强随机噪声的地震数据(图 1b)，后两种方法对噪声压制效果差，严重干扰对断层的识别．图 6断层检测结果表明本文所提方法在理论模型中进行断层检测是可行且有效的，而且本文所提方法的抗噪性最优．

4 实际数据分析

选取墨西哥湾某地区的二维时间偏移剖面对本文所提方法进行测试(Claerbout, 2010)(图 4a)．该地震数据的时间深度范围为0.4~1.8 s，有350个采样点，采样间隔为0.004 s，测线长度为200道地震道．

用基于加权一致性的蚁群算法对图 4a进行断层检测．参数设置为：分割地震数据的小数据块大小为2×2；终止条件S_min=0.2、S_max=0.5；各数据点的信息素浓度初始值为1，信息素更新系数c=0.5；信息素和启发函数的重要程度因子分别是α=1、β=3；信息素挥发系数ρ=0.1；迭代次数为20次．检测结果如图 7a所示，结果能清晰且较完整地指示出断层位置及延伸长度，并较好地压制了噪声和地层残余响应这两类伪断层信息的影响．为了进一步比较，给出由基于方差属性的蚁群算法和基于C3相干属性的蚁群算法对图 4a计算而得到的断层信息(图 7b、图 7c)．因为方差算法和C3相干算法对噪声和地层残余响应的压制效果差，致使蚁群追踪到的断层信息较模糊且不完整，对于地层信息匮乏的断层区域，这两种方法并不能检测出应有的断层信息．图 7的对比结果证明了本文所提方法在实际地震数据中进行断层检测的可行性和优越性．

5 计算效率分析

为了分析本文所提方法进行断层检测的计算效率，本文列出基于加权一致性的蚁群算法、基于方差属性的蚁群算法和基于C3相干属性的蚁群算法对相同地震数据(图 1a，图 1b，图 4a)进行断层检测的运行时间表(表 1)．本文运行程序所用电脑配置为：操作系统：Windows 10专业版64位；处理器： Inter(R) Core(TM) i3-3220 CPU@ 3.30GHz(4CPUs)，~3.3 GHz；内存：4096 MB RAM．三种方法在蚁群追踪环节的参数设置相同，由表 1可见，本文所提方法的计算效率最高．

6 结论

本文提出了新的断层自动检测方法——基于加权一致性的蚁群算法，得出了以下结论：

(1)  一致性是对指纹图像中指示局部方向场信息的各向异性强度的一个度量，将其应用于地震数据处理，可指示地层的不连续性．

(2)  提出加权一致性的概念，与一致性相比，前者具有更高的抗噪性能，更大程度地消除了不整合面和地层残余响应这两类伪断层信息，更能突显地层的横向不连续性．

(3)  对于相同的断层检测过程，对比常规的基于方差属性的蚁群算法及基于C3相干属性的蚁群算法，本文所提方法具有以下优点：抗噪性更高，对伪断层信息压制效果更好，断层边界收敛能力更强；在增加断层连续性的基础之上，对断层的延伸长度展现的更完整；稳定性和计算效率更高．