断层识别是地震资料解释、储层预测和描述的重要基础。断层识别的精度决定构造成图的精度,因此, 如何有效地从地震数据中进行断层的自动拾取, 一直是地震数据解释研究的热点问题[1-2]。常用的断层识别方法大致分为以下几种: 常规的属性识别法、自动追踪解释法、图像处理识别法和机器学习识别法。属性识别法一般选取对断层属性比较敏感的统计量来识别断层, 如相干体属性、方差、灰度中值和曲率等[3-4]。自动追踪解释法包括蚂蚁追踪法、稀疏脉冲反褶积法以及边缘检测技术类方法[5-6]。图像处理识别法包括基于特征映射法、基于波形特征类方法和图像分割类方法等[7-8]。近年来, 不少学者使用机器学习类的方法(如卷积神经网络、深度学习和聚类网络等), 对储层构造属性进行识别[9-10]。这几类方法有各自的优、缺点和适用条件, 在断层解释工作中, 一般根据实际地震资料的具体情况选取合适的断层识别方法[2,11]。
方向滤波和边界保持滤波是图像增强处理中常用的技术, 其物理过程是使用与图像局部方向一致的滤波器进行平滑处理, 压制噪声并增强图像沿纹线方向的连续性。边界保持的作用是在滤波的同时对图像边缘进行保护, 避免边缘信息因滤波而模糊化。本文提出基于方向滤波和边界保持滤波实现的三维叠后数据增强方法, 通过方向滤波增强地震同相轴的信噪比和连续性, 通过边界保持滤波来保持断层信息, 避免因使用断层两侧的数据点进行平滑导致断层模糊化。在此基础上, 本文对该方法做两点改进: 1) 用多方向断层识别代替常用的单方向断层识别, 可在一定程度上避免将倾斜地层错误地识别为断层, 更好地适用于存在倾斜地层的情况; 2)将二维边界保持滤波方法扩展到三维地震数据, 同时保证沿测线(inline)和垂直于测线(crossline)地震剖面数据的连续性。最后, 将本文方法应用于合成数据和实际工区地震资料, 检验方法的有效性。
1 方法原理
1.1 方向滤波和边界保持滤波
方向滤波和边界保持滤波是广泛应用的图像增强方法。方向滤波是沿着图像方向场的纹线方向进行滤波, 使图像沿纹线方向的数据得到平滑, 垂直于纹线方向的数据得到分离, 达到压制噪声和地层背景干扰的效果[12-13]。方向场代表图像中每个像素点所在的脊线或谷线在该点的切线方向, 是图像处理的重要识别特征之一。将方向场应用于图像信息增强, 可以最大限度地保持图像的分辨率, 并实现压制噪声的效果。方向场提取的基本思想是, 计算灰度图像中每一点在所有方向上的统计量特征,并根据各个方向上统计量的差异确定该点的方向,常用的方法包括基于点方向图、块方向图或连续分布方向图的方向提取法。基于点的方向场提取法计算速度快, 但对噪声比较敏感, 提取的方向场准确性不高。本文使用基于块的方向场提取法, 常用的是梯度法, 即根据数据图像中每一点在各个方向上的梯度值确定该点的方向。基于梯度的方法所需运算量大于基于点的方法, 但梯度的连续性可以保证所提取方向场的连续性。
梯度法实现方向场提取的过程如下: 数据图像I(x,y)的梯度g(x,y)大小表示为
梯度的方向表示为
也就是说, 可以将图像与高斯函数的一阶导数进行卷积, 得到图像的梯度信息。进一步地, 将梯度与其转置相乘, 构造梯度结构张量
其最大特征值对应的特征向量中包含地震图像的方向场信息。得到方向场信息后, 就可以对图像进行方向滤波, 即沿着每个点对应的方向场数据进行平滑处理, 得到压制噪声的效果。
针对断层区域进行方向滤波时, 可能造成断层边缘模糊化。由于断层两侧的地层信息存在不连续性, 如果强行将断层两侧的数据点进行平滑处理,会使得断层两侧的地层变得连续, 断层信息被破坏。本文使用边界保持滤波来避免这一问题。常用的边界保持滤波方法包括 Kuwahara 滤波、扩散滤波和中值滤波等[12,14-15]。本文采用 Kuwahara 滤波,其思路是计算目标点所有邻域内的特征值, 选择特征值中最大的数值作为该点的目标值。Kuwahara滤波可以通过在边界处移动滤波器来避开边界, 对边界起到较好的保护作用。具体应用时, 需要事先识别断层(即边界)位置, 当滤波器包含断层时, 需要移动滤波器来避开断层(图 1)。
图1 边界保持滤波示意图
Fig. 1 Schematic diagram of edge-preserving filtering
本文通过相似度信息 C 来识别断层位置:其中, νi 是将局部图像数据通过奇异值分解得到的特征值, 按降序排列, 因此 ν1 是最大特征值。相似度信息可用于表征图像中是否存在断层, 对于不存在断层的水平连续地层, ν1=1, 其他特征值均为 0,因此 C=1; 对于存在 N 个断层的情况, 将水平地层切割成 N+1 段, 则将奇异值分解后得到的前 N+1 个特征值远大于其他特征值, 此时 C＜1。图 2 为一个简单相似度信息计算示例, 对于水平地层剖面, 将奇异值分解后得到的特征值为[2.7029, 0, 0, …, 0],因此 C=1 (图 2(a)); 对于存在断层的地层剖面, 将奇异值分解后得到的特征值为[2.1250, 1.5249, 0.6815,0, 0, …, 0], 因此 C=0.4906＜1 (图 2(b))。实际应用时, 一般选取目标点邻域的局部图像来计算相似度信息。另外, 考虑到实际地震资料中可能存在倾斜地层和背景噪声干扰, 即使局部图像中不存在断层,C 值也可能略小于 1, 一般选取阈值 λ=0.95, 当 C＜λ时, 则表明目标点邻域的局部图像中存在断层。
图2 相似度信息计算示例
Fig. 2 Examples of continuity calculation
1.2 多方向断层识别和三维断层增强
我们针对边界保持滤波做两点改进——多方向断层识别和三维断层增强, 使得本文方法能够更好地适用于实际叠后地震资料的处理。
在计算相似度信息 C 的时候, 只要地震资料在水平方向不连续, 都可能导致 C＜1, 即被识别为存在断层, 因此倾角较大的倾斜地层可能被错误地识别为断层。为适应倾斜地层的情况, 本文使用多方向断层识别的方法, 与常规的沿水平方向的单方向断层识别方法相比, 多方向断层识别是沿着一定范围倾角的倾斜方向计算 C 值。实现过程如下: 选取一系列倾角 -θ n , -θ n -1 ,...,0,...,θ n -1,θn , 其中最大倾角θn 根据实际地震剖面中观察到的最大地层倾角确定; 针对倾角 θi, 以目标点为中心, 将局部图像进行角度为 θi 的旋转, 对旋转后的图像数据计算 Ci(图 3), 选取 Ci 的最大值作为目标点的 C 值。通过多方向断层识别, 倾角小于 θn 的地层就不会被识别为断层, 增强了本文方法对倾斜地层的适用性。
图3 多方向断层识别示意图
Fig. 3 Schematic diagram of multi-directional fault recognition
进一步地, 本文将二维断层增强方法扩展为三维断层增强方法, 以期更好地适用于叠后三维地震资料的处理。在二维断层增强方法中, 需要对每条inline 的地震剖面进行滤波操作, 相当于沿 inline 方向修改数据点的数值, 但可能造成 crossline 的地震剖面数据不连续, 表现为出现异常噪声或地层不连续, 进而影响断层识别的效果, 因此有必要将二维断层增强方法扩展为三维断层增强方法。严格的三维断层增强方法需要在三维空间进行方向角计算和三维空间的断层识别, 当空间采样不均匀或空间采样间隔较大时, 还面临三维数据空间插值的问题,会产生巨大的运算量。
本文采用简化方式实现伪三维的边界保持滤波, 可以避免严格的三维方法所需的巨大运算量:针对三维数据体的目标点, 分别在该目标点沿 inline 和 crossline 的地震剖面确定滤波器对应的数据点, 将这两组数据点共同用于滤波计算, 得到的数值作为目标点数值, 这样就可以实现简化的三维方向滤波。这种简化的三维方法是直接对二维方法进行扩展, 方向场计算、断层识别和边界保持滤波所用数据点范围的确定均在二维方向进行, 不涉及三维数据运算(仅最后使用多个方向的数据点进行滤波时涉及三维数据运算), 因此运算量比严格的三维方法大大减少。该三维方法既能达到滤波效果,也能保持 inline 和 crossline 的地震剖面包含的断层信息, 避免出现某个方向数据不连续的现象。
1.3 总流程
如图 4 所示, 本文方法分为以下 4 个主要步骤。
图4 三维断层增强方法流程
Fig. 4 Main procedures of three-dimensional fault enhancement technique
1) 利用梯度法计算叠后地震剖面的方向场。这一步涉及将图像数据与高斯函数的一阶导数进行卷积, 经参数测试, 对于采样率为 1 ms, 主频为 30 hz 的地震数据, 标准差大于 4 ms 的高斯函数, 均可保证计算得到的方向信息的准确性, 因此本文所有实例的高斯函数标准差均选取 5 ms。
2) 多方向断层识别, 需要进行参数测试, 确定3 组参数。第一组参数: 计算相似度信息时选取的局部数据时窗大小, 包括空间宽度和时间长度, 时窗大小的选取需要适中。对于空间宽度, 过小则对随机噪声过度敏感, 在断层识别结果中引入噪声;过大则在时窗内引入过多的水平同相轴信息, 导致断层识别模糊。对于时间长度, 过小则局部数据包含的断层信息较少, 降低断层识别的精度, 过大则会引入较大的运算量。第二组参数: 最大倾角数值,一般根据在实际地震剖面中观察到的最大地层倾角来确定。第三组参数: 判断是否存在断层的阈值 λ。
3) 根据前两步得到的方向场信息和断层信息,实现边界保持滤波。滤波前需要通过参数测试确定滤波器的时窗大小, 时窗大小的选取需要适中, 保证既可以达到滤波效果, 又避免数据的过度平滑。
4) 在二维边界保持滤波的基础上, 结合 inline和 crossline 地震剖面数据, 实现三维边界保持滤波。
2 合成数据实例
我们通过合成数据对比, 验证本文方法的有效性。所用数据为半合成数据, 截取实际地震剖面中的连续水平地层数据, 通过对数据图像进行旋转,模拟倾斜地层情况, 通过在选定的空间位置对数据进行时间方向的位移, 模拟断层的分布情况, 最终形成的合成数据接近断层发育区域的实际生产数据。经参数测试, 确定该实例的最优参数如下: 断层识别所用时窗大小为 40 ms×105 m, 地层最大倾向为 40 ms/km, 相似度阈值信息为 0.9, 滤波时窗长度为 5 个数据点。原始地震剖面的断层模糊, 且存在背景噪声干扰(图 5(a)), 经本文方法增强断层后的剖面信噪比显著提高, 且断层边界更清晰(图5(b))。图 5(c)和(d)显示, 断层增强后的相干剖面对断层的刻画更清晰、更准确。由于本文方法使用多方向断层识别, 在断层增强后的剖面中, 原始相干剖面中因存在倾斜地层而产生的虚假断层信息被有效地压制。
图5 断层增强前后的结果对比
Fig. 5 Results before and after fault enhancement
3 实际数据应用
应用本文提出的断层增强方法处理工区 A 的三维叠后地震资料。经参数测试, 确定该实例的最优参数如下: 断层识别所用时窗大小为 30 ms×105 m,地层最大倾向为 25 ms/km, 相似度阈值信息为 0.95,滤波时窗长度为 7 个数据点。图 6 展示 inline 地震剖面的断层增强效果。原始剖面存在明显的随机噪声干扰, 且断层边缘较模糊, 增强断层后, 数据信噪比和地层同相轴的连续性明显提高, 断层边缘更加清晰。图 6(c)和(d)显示, 在应用本文方法进行断层增强后的相干剖面中, 噪声得到显著的压制, 原本模糊不清的断层信息被清晰地刻画出来(虚线椭圆所示区域)。
图6 inline 地震剖面及相干剖面对比
Fig. 6 Seismic profile along inline and comparison of coherence slices
图 7 展示相干切片的对比。原始相干切片受到背景噪声干扰, 相干体属性难以识别(图 7(a)), 应用本文方法进行断层增强后, 相干切片的信噪比明显提高, 相干体属性刻画更加清晰, 原本模糊不清的相干体属性得以清晰地呈现(图 7(c)中虚线椭圆所示区域)。由于工区 A 存在较多的倾斜地层(图 6(a)), 与多方向断层识别结果相比, 单方向断层识别结果中出现较多的噪声干扰, 相干体属性也相对模糊(图 7(b)中虚线椭圆所示区域)。因此, 对于倾斜地层, 与常用的单方向断层识别方法相比, 多方向断层识别方法更适用。
图7 基于单方向识别和多方向识别的方法在相干切片上的断层增强效果对比
Fig. 7 Comparison of the single- and multi- directional methods of fault enhancement on coherent slices
最后, 对比三维增强方法与二维增强方法的处理效果。图 8 展示 crossline 地震剖面的断层增强效果。与原始地震剖面相比, 使用二维方法处理 inline 方向的地震剖面, 信噪比有所提升, 但仍然有较多的背景噪声残留, 主要原因是仅对 inline 的地震剖面进行二维数据增强, 可能引起 crossline 的地震剖面数据不连续。与二维增强结果相比, 在应用本文的三维方法增强断层后的地震剖面中, 背景噪声被显著地压制, 地层同相轴的连续性得到提高。
图8 二维方法和三维方法在crossline 地震剖面上的断层增强效果对比
Fig. 8 Comparison of the two- and three- dimensional methods of fault enhancement on crossline seismic profiles
相干切片的对比结果(图 9)显示, 原始相干切片存在较多的噪声干扰(图 9(a)); 二维增强方法可以在一定程度压制噪声, 对相干体属性的刻画更加清晰, 但图 9(b)中虚线椭圆所示区域的相干体属性仍然模糊; 三维增强技术的处理效果最佳, 相干切片信噪比进一步提高, 相干体属性更加清晰, 图9(c)中虚线椭圆所示区域的两组相干体属性得以清晰地呈现。
图9 二维方法和三维方法在地震剖面和相干切片上的断层增强效果对比
Fig. 9 Comparison of the two- and three- dimensional methods of fault enhancement on seismic profiles and coherent slices
图 9(d)～(f)分别展示图 9(a)～(c)中虚线所示测线对应的地震剖面, 通过地震剖面的对比再次验证本文方法的断层增强效果。从图 9(d)～(f)的地震剖面中可以看到明显的“断层发育”现象, 但以此判断相干切片(图 9(a)～(c))中的相干体属性能够反映断层信息是不够严谨的。相干体属性表征的是地震反射特征突变的位置, 可能由断层、河道边界、火成岩边界以及砂泥分界面所导致。工区 A 位于渤海湾盆地, 该区域断层系统异常发育, 所以相干体属性对应断层信息的概率较大, 但在缺少进一步证据的情况下, 仍无法排除相干体属性对应砂-泥边界或其他情况。判断相干体属性的地质成因并非本文的研究重点, 因此仅做此说明, 不做深入的论证。
4 结论
本文提出一种基于多方向识别的三维断层增强方法。通过方向滤波增强地震同相轴的连续性, 并压制背景噪声; 用多方向断层识别代替常用的单方向断层识别, 提高本文方法对存在倾斜地层情况的适用性; 通过边界保持滤波, 保持地震剖面中的断层信息; 将二维断层增强方法扩展到三维, 实现以较低的运算量达到三维断层增强的效果。合成数据和实际工区数据的应用实例均证明了本文方法的有效性, 能够在有效地压制背景噪声的同时, 得到断层分辨率更好的相干剖面, 有助于进一步的构造解释以及油气输导体系评价。
[1] Gersztenkorn A, Marfurt K J. Eigenstructure-based coherence computations as an aid to 3-D structural and stratigraphic mapping. Geophysics, 1999, 64(5):1468-1479
[2] 李婷婷, 侯思宇, 马世忠, 等. 断层识别方法综述及研究进展. 地球物理学进展, 2018, 33(4): 1507-1514
[3] 郑静静, 印兴耀, 张广智, 等. 基于 Curvelet 变换的多尺度分析技术. 石油地球物理勘探, 2009, 44(5):543-547
[4] 杨国权, 刘延利, 张红文. 曲率属性计算方法研究及效果分析. 地球物理学进展, 2015, 30(5): 308-312
[5] 宋建国, 孙永壮, 任登振. 基于结构导向的梯度属性边缘检测技术. 地球物理学报, 2013, 56(10): 345-355
[6] 刘财, 刘海燕, 彭冲, 等. 基于加权一致性的蚁群算法在断层检测中的应用. 地球物理学报, 2016, 59(10): 3859-3868
[7] 朱剑兵, 赵培坤. 国外地震相划分技术研究新进展.勘探地球物理进展, 2009, 32(3): 167-171
[8] Chopra S, Marfurt K J. Seismic facies classification using some unsupervised machine learning methods //88th Annual Meeting International Conference and Exhibition, SEG Expanded Abstracts. Anaheim, 2018:2056-2060
[9] 付超, 林年添, 张栋, 等. 多波地震深度学习的油气储层分布预测案例. 地球物理学报, 2018, 61(1):293-303
[10] Yuan Zhenyu, Huang Handong, Jiang Yuxin, et al. An enhanced fault-detection method based on adaptive spectral decomposition and super-resolution deep learning. Interpretation, 2019, 7(3): 1-63
[11] 郎晓玲, 彭仕宓, 康洪全. 碳酸盐岩缝洞型储层地球物理响应特征及预测方法研究. 北京大学学报(自然科学版), 2012, 48(5): 775-784
[12] Bakker P. Image structure analysis for seismic interpretation [D]. Delft: Delft University of Technology.2002
[13] 杨培杰, 穆星, 张景涛. 方向性边界保持断层增强技术. 地球物理学报, 2010, 53(12): 2992-2997
[14] 刘财, 李红星, 陶春辉, 等. 模糊嵌套多级中值滤波方法及其在地震数据处理中的应用. 地球物理学报, 2007, 50(5): 1534-1542
[15] Lavialle O, Pop S, Germain C, et al. Seismic fault preserving diffusion. Journal of Applied Geophysics,2006, 61(2): 132-141