基于倾角校正的地震层位追踪算法

刘 鑫， 车翔玖， 林森乔

（吉林大学计算机科学与技术学院，吉林 长春 130012）

基于倾角校正的地震层位追踪算法

刘 鑫， 车翔玖， 林森乔

（吉林大学计算机科学与技术学院，吉林 长春 130012）

由于地下构造复杂且原始采集地震数据的信噪比低、噪声干扰强等原因，常无法正确拾取地震层位，如出现层位跳变等问题。为了解决该问题，本文提出了一种抗干扰的地震层位追踪算法，首先基于Morlet小波进行滤波去噪，并进行相干计算以获得倾角控制体；然后使用倾角信息进行初始追踪，加快增长速度；若初始追踪失败，根据同相轴类型，在搜索时窗内继续追踪候选点，追踪过程利用倾角信息获得正确的邻域进行增长判断；本文还使用了多个容许偏差值进行逐步追踪的同时，通过权值参数控制倾角在断层处停止增长。实验表明本文算法识别的层位面与原始数据的地质层位吻合良好，且断点位置较准确。

层位识别；地震解释；相干体；小波变换

地震解释的主要工作是对地震资料进行处理、分析和研究，进而推断地下构造和岩层性质，寻找可能的含油构造，如断层。地下的岩层在地壳运动中，受力达到一定强度时会发生破裂、错位，从而形成了断层。层位是不同地层之间的分界面，因此断层的识别可通过判断其与层位的交点来进行。最开始的断层解释方法是在层位识别的基础上进行的，通过检测各层位的不连续性，对断层进行解释和组合。随着地震相干体技术的发展，使得基于相干体的断层自动识别成为可能。但层位追踪仍然是地震解释中基础而关键的一步，因为将层位追踪结果与断层识别结果放到三维地震体内结合可以对断层进行验证。

国内外已有相关层位识别技术的文献，主要可分为以下几种：①边缘检测法[1-2]，将地震体数据看作二维灰度图像序列，通过图像处理中常用的边缘检测方式来提取同相轴，常用的有：Canny之类的边缘检测算子、导数、数学形态学等。由于地下构造复杂且信噪比低，边缘检测结果含噪声大，往往很难直接用于地震解释；且该方法属于二维地震解释技术，得到二维识别结果序列后，地质解释人员据此来推测、想象地下岩层空间构造，然而相邻两片之间的信息很难建立联系，因此基于二维切片的层位追踪往往需要反复修改和验证，导致整个地震层位的解释效率很低。②插值法[3]，是手动拾取的简化，该方法不需要逐点拾取整个层位，而是离散地手动拾取部分层位点作为控制点，然后通过插值法对控制点进行插值得到完整层位。该方法结果的准确性受控制点集规模大小的影响较大，控制点集过大使得算法趋近手动拾取过程，人工任务繁重、效率低，且过多的人为干预使得结果随机性增大；控制点集过小，则大部分原始地震信息没得到利用，造成地层结构的三维细节遗失，拾取结果往往不可信。③神经网络法[4-5]，其工作过程包括学习（训练）和使用（预测）两部分。首先选取个别地震道进行层位解释，将其解释结果作为标准样本来训练网络。训练过程耗时较长，但训练过程结束后则可使用该网络进行层位的自动识别。由于神经网络学习过程需要有足够多的分布均匀的学习样本，加大了样本选择的难度，且当使用过程中发现预测结果较偏离真实结果时还要再拾取新样本返回训练过程，反复的迭代运算使得整个层位识别效率降低。④区域增长法，不同于二维层位解释技术，利用区域连通性实现三维地震层位的自动识别。该算法原理简单，无需训练过程，只需从若干种子点出发，根据预先定义的连接标准（一般为地质特征或局部波形的相似程度）自动增长至所有满足连接标准的像素集合，从而完成层位的追踪。这些优点区域增长法是国内外研究学者和商业地震解释软件中常用层位追踪方式。

然而基于种子点的区域增长算法仍有些亟待解决的问题，例如种子点的选取问题、噪声干扰问题、以及连接标准选取的好坏对结果影响很大等。为了解决这些问题，国内外研究者提出了不同的改进算法。针对连接准则选取问题，Revol和Jourlin[6]提出了基于最小偏差的区域增长法(R-RGA)。这种方法用局部数据直方图的紧缩或者扩增使得已分割的区域符合一致性准则Gmax要求。并在此基础上进一步提出了基于评价函数的区域增长算法(RM-RGA)[7]，以解决之前的 Gmax的自动选取问题。由于Gmax是随机选取的且种子点选择是随机的，而种子点选择不当会导致错误分割的问题，针对这一问题马仁安等[8]提出基于局部信息的区域增长算法(LI-RGA)，该算法根据原始选择的种子点的邻域信息来确定真正的种子点和Gmax。针对层位实质上是三维空间中的曲面结构而非三维体结构这一特征，Faraklioti和Petrou[9]基于面检测技术，同时结合区域增长的连通性实现了单体素厚度层位的识别。针对噪声干扰和层位性质相近造成的层位窜位问题，Patel等[10]提出将识别层位转化为三维网格模型的方法，通过分割层位网格模型，地震解释专家借助于相关专业知识能够分离错误层位连接。Hoyes和Cheret[11]对Petrel、Opendtect等商业地震解释软件中的层位解释技术进行了综述，其中Opendtect主要是基于倾角控制的区域增长层位识别方法；Petrel是基于局域波形相似性的区域增长算法，生成层位序列后，再由解释人员组合层位片段。波形的相似性通常由互相关计算进行度量，为了降低噪声对相似值计算的影响，利用噪声的高阶累积量恒为零的特点，近年来出现了以高阶累积量代替互相关计算[12]。

上述算法中无论是基于阈值连接还是波形相似度，若采用“单连接”形式，即在生长过程中从某个已确定为层位点的体素出发，其邻域内的所有点都与该点进行连接判断。这种“单连接”方式由于只有两道数据参与判断，受噪声影响较大；若采用“全局连接”形式，将单像素的值与其邻域内的平均值进行比较，算法的稳定性得以提高。区域生长的原理是从已确定的区域根据连接准则向未确定的区域进行扩张生成目标物体，若采用“全局连接”形式，由于地震层位往往是非规则的，为了解决邻域确认这个问题，本文提出利用相干体算法计算的倾角对邻域进行校正。本文还提出了使用多个容许偏差值进行逐步追踪的同时，通过权值参数控制倾角在断层处停止增长。

传统区域增长方法向邻道扩张时，未考虑层位倾角，当搜索时窗内有相似同相轴时容易发生层位窜层的错误。本文采用国外著名的Opendtect软件提供的荷兰 F3演示数据，该数据最明显的特征是底部有一个大型S形层面。图1为传统区域增长算法对该S层面的追踪结果，图2为放大“大断层”处的某一切片。图2中正确的层位应该为红色箭头所指的层面，但由于搜索时窗内黄色箭头所指的层面与目标层位相似，发生了窜层，便产生了图方框中的伪断层的现象。为了解决这一问题，本文对传统相干体算法进行了改进，通过小波域能量分布自适应确定时窗长度并计算得到的倾角。

图1 传统算法下的层位追踪整体效果(蓝色框内为两个可能的大断层，上方蓝色框处的层位追踪走向错误)

图2 “大断层”处的某一切片(红色箭头所指为正确层位走向，黄色箭头所指为与正确层位相似的层面)

1 本文算法描述

1.1 预处理

区域增长算法对噪声较敏感，所以需要先抑制噪声，提高数据的信噪比再进行层位解释工作。现代商业地震解释软件中最常用的是倾角中值滤波，除此之外邓小英等[13]提出了基于卡尔曼滤波的方法以提高地震资料的信噪比进行层位追踪。但该算法需要对噪声和信号的统计特性做一定的假定，然后通过合适的数学方式提供信噪比。地震波在地层介质传播的过程中，由于复杂的地下结构以及地震波的散射和吸收衰减等因素，造成地震信号通常是以非平稳信号形式呈现的，时频分析恰恰是这类非平稳信号分析的有效手段。常用的时频分析方法有短时傅里叶变换[14]、Gabor变换[15]、Wigner-Ville时频分布函数[16-17]、小波变换[18-20]等方法。小波变换不仅是图像处理领域常用的去噪滤波方式[21-22]，因其多分辨分析的优势近年来成为了新一代相干体技术的代表之一。因此本文选择小波阈值法进行去噪滤波，利用小波变换将地震数据分成几个频带进行滤波处理，有选择地突出特定频带，在去噪的同时保留了有效频段，能够提取隐藏在地震资料中的有用信息，提高地震解释的准确度。

1.2 小波阈值去噪

小波变换具有很强的去数据相关性。在小波域内，信号的能量集中在大的小波系数中，而噪声的能量却分布于整个小波域内，因此，经小波分解后，幅值比较大的小波系数一般以信号为主，而比较小的系数在很大程度上是噪声。因此，采用阈值法可以保留信号系数，减少大部分噪声系数。设地震道信号为f (t)，其连续小波变换W(a, b)定义为：

其中，a为尺度因子，对应频率信息；b为平移因子，对应时空信息；是母小波ψ的复共轭。

傅里叶变换F(ω)为：

结合式(1)和式(2)，则式(1)可重新写成：

由式(3)可知信号f (t)的小波系数求取可通过两步完成：①计算信号f (t)的傅里叶频谱F(ω) 及母小波的复共轭函数(t)的频谱(aω)；②计算F(ω)(aω)的傅里叶逆变换。

将信号在不同尺度上进行小波分解并得到相应的小波系数后，保留所有低频系数，对高频系数｛ωa,b｝采取硬阈值法进行处理，选取阈值其中N是信号的长度，则处理后的高频系数｛a,b｝为：

图3 原始地震信号与小波去噪后效果对比

小波逆变换可写为：

图3(a)为原始地震信号，包含有大量噪声，经小波阈值法去噪后效果如图3(b)所示，地震数据中噪声明显降低，有利于进一步层位追踪。

1.3 自适应时窗相干体

高怀静等[23]针对地震资料处理提出了改进的Morlet小波，能够适应于各种地震子波的分析，且morlet是复小波，能够有效地提取地震信号中的幅值和相位信息。其改进的Morlet小波的定义为：

其中，fc为中心频率，fb为带宽。不同的fc，fb可以得到不同形态的小波函数，fc值随着尺度因子a的增大向低频方向移动。实际应用时需要根据不同的地震子波记录，选择适当的fc，fb，得到最接近地震子波的小波函数，用其对原始地震数据做小波分频处理。小波相干体的具体计算步骤如下：

步骤 1. 原始三维数据体中每道都用式(5)作为母小波做小波变换，并根据小波域计算瞬时特征参数，再按数据体排序，即可得到多个分频数据体V(ai)。

步骤 2. 对每个分频数据体用二代或三代相干体算法计算相干值，得到分频相干体C(ai)。

步骤 3. 根据式(6)重构相干体，通过重构系数di，对一定频带相干体进行放大或减小，突出特定频段的相干体。

在步骤 2中计算相干体时以计算点为中心选取固定长度时窗很难兼顾浅层和深层，即在分析浅层时出现跨多个同相轴，在分析深层时因时窗长度不够，在同相轴零点附近相干值易受随机噪声影响。

传统相干体算法采取固定时窗在全频带上计算相干体，无法兼顾一道信号低频到高频的波动，即分析高频段时出现跨多个地震子波，而在分析低频段时，无法显示一个完整的波峰或波谷。为了解决固定时窗带来的问题，本文提出在利用小波去噪的同时，根据小波域能量集中区域对应的频率计算对应的窗口值，可随波形的变化选择不同大小的时窗。设某一地震道上某一采样点经小波变换至不同频率的小波域中，该样点在各个分频数据体中具有不同的系数，其中最大的系数对应的频率为fmax，则该样点相应的时窗长度L为：

图4中红线代表地震道对应的小波能量图如图4(a)所示，逐点根据小波域能量集中区域对应的频率计算对应的窗口值如图4(b)所示。图中可以看出层位均匀区域能量集中在高频部分，对应小时窗；嘈杂区域能量集中在低频区域，对应大时窗，可降低噪声对相干计算的干扰。利用自适应时窗方法求取相干体后，根据文献[24]，可以进一步计算各点倾角，得到倾角控制体。

图4 自适应时窗确定

1.4 改进后的层位追踪算法

通过鼠标与地震剖面交互设置识别层位的种子点后，在搜索时窗内查找离交互点最近的满足搜索特征的正确种子点，其中满足追踪特征的种子点有波峰、波谷和零交叉点。本文的追踪算法具体实现如下：

步骤 1. 算法使用倾角控制体中的倾角信息进行初始追踪，以加快增长速度，若初始追踪满足追踪特征则作为候选点；

步骤 2. 若初始追踪失败，即没有在邻域内找到满足追踪特征的点，则进一步判定是否处于同一同相轴，在搜索时窗内继续追踪候选点，追踪过程中利用倾角信息获得正确的邻域进行增长判断，即若追踪的点满足同相轴特征并且倾角大小在倾角信息允许范围之内，则将该点加入种子点并标记其追踪过，否则直接将该点标记为追踪过；

步骤 3. 追踪过程中，使用多个表征倾角范围和振幅范围的容许偏差值进行逐步追踪。层位较为均匀时，倾角较小且相邻点振幅相差小，因此最初选择较小的容许偏差值进行追踪，使较为平整均匀的层位部分先被识别出来，然后在后续中逐步增大容许偏差值并以前一轮偏差值所追踪到的结果作为种子点进行追踪，直至如果所有的地震道都被追踪到或超过了最大容许偏差值则追踪结束；

步骤 4. 为了更好地对断层结果进行验证，希望层位的增长能够在断层处停止，断层两侧层位倾角较大且振幅相差很大，因此，在增长时，通过式(7)对增长进行控制，其中，vcur为当前增长容许偏差值，vmax为最大容许偏差值，NumOfCandidate为候选点个数，NumOfNeiborhood为邻域内未被标记追踪点的个数。

2 实验结果分析和比较

本文实验环境如下：CPU：Inter(双核，主频2.0 GHz), 内存：1 G, GPU：NVIDIA GeForce 8400 M GS(2个SM，每个含8个SP，显存128 MB)。编程环境：Visual Studio2005+OpenGL+CUDA。本文算法的实验数据为 Netherlands Offshore F3 Block，该数据总共有631×951道，每道有462个采样点。

图5是传统算法下层位追踪的结果，因层位倾角过大，领域信息不准确，追踪过程中，在搜索时窗内遇到与正确层位相似的同相轴，产生窜层的现象，导致层位追踪不完整；图6是本文算法下层位追踪的结果，进行追踪时，首先利用倾角计算每一点的领域，在遇到相似同相轴时，没有直接停止追踪，而是根据倾角信息，判断相似同相轴处的点是否满足倾角范围，增强领域确定的准确性，最终层位追踪在断点处停止，结果显示，本文算法下的层位追踪结果避免了受搜索时窗内错误层位的相似同相轴的误导，比传统算法下的层位追踪结果更完整、准确。图7是本文算法下的层位追踪整体效果图，与图1相比，错误的层位得到纠正，层位追踪的准确度极大提高。图6与图7共同验证了本文算法比传统算法更适合层位追踪。

图5 传统算法下层位追踪效果(红圈处为实际倾角较大的层位，未使用倾角校正，层位追踪在搜索时窗内倾角较小的相似同相轴处停止，造成错误追踪结果)

图6 本文算法下层位追踪效果(红圈处为实际倾角较大层位，因使用倾角校正，层位追踪在正确的同相轴处停止，追踪结果正确)

图7 本文算法下的层位追踪整体效果

3 总结与展望

为了解决地震数据体信噪比低造成的层位错误追踪问题，本文提出了一种抗干扰的地震层位追踪算法，首先利用小波变换滤波去噪，同时利用倾角信息计算初始追踪点，若初始追踪满足追踪特征则作为候选点，否则在搜索时窗内进一步使用倾角信息追踪候选点，最终实现层位的准确追踪。但是，本文算法在追踪速度上仍需要继续改进，下一步工作重点考虑使用并行化等方法提高断层追踪速度。

[1] 李红星, 刘 财, 陶春辉. 图像边缘检测方法在地震剖面同相轴自动检测中的应用研究[J]. 地球物理学进展, 2007, 22(5): 1607-1610.

[2] Li Lili, Ma Guoqing, Du Xiaojuan. New method of horizon recognition in seismic data [J]. Geoscience and Remote Sensing Letters, IEEE, 2012, 9(6): 1066-1068.

[3] 李学森. 地震反射层位构造信息三维可视化显示技术研究[J]. 地球物理学进展, 2005, 20(3): 735-740.

[4] 姚 姚. 用人工神经网络实现同相轴自动拾取[J]. 石油地球物理勘探, 1994, 29(1): 111-116.

[5] Leggett M, Sandham W A, Durrani T S. Geophysical applications of artificial neural networks and fuzzy logic [M]. Springer: Netherlands, 2003: 31-44.

[6] Revol C, Jourlin M. A new minimum variance region growing algorithm for image segmentation [J]. Pattern Recognition Letters, 1997, 18(3): 249-258.

[7] Revol-Muller C, Peyrin F, Carrillon Y, et al. Automated 3D region growing algorithm based on an assessment function [J]. Pattern Recognition Letters, 2002, 23(1): 137-150.

[8] 马仁安, 张二华, 杨静宇, 等. 不规则地质体的分割与体绘制方法研究[J]. 计算机研究与发展, 2005, 42(5): 853-887.

[9] Faraklioti M, Petrou M. Horizon picking in 3D seismic data volumes [J]. Machine Vision and Applications, 2004, 15(4): 216-219.

[10] Patel D, Bruckner S, Viola I, et al. Seismic volume visualization for horizon extraction [C]//Pacific Visualization Symposium (PacificVis), 2010 IEEE. IEEE, Taipei, 2010: 73-80.

[11] Hoyes J, Cheret T. A review of global interpretation methods for automated 3D horizon picking [J]. The Leading Edge, 2011, 30(1): 38-47.

[12] 冯智慧, 刘 财, 冯 晅, 等. 基于互四阶累积量一维切片的地震层位自动拾取方法[J]. 石油地球物理勘探, 2011, 46(1): 58-63.

[13] 邓小英, 胡 健, 李 月, 等. 一种新的基于卡尔曼滤波的地震记录同相轴跟踪方法及性能分析[J]. 地球物理学报, 2014, 57(1): 270-279.

[14] Lu Wenkai, Li Fangyu. Seismic spectral decomposition using deconvolutive short-time Fourier transform spectrogram [J]. Geophysics, 2013, 78(2): 43-51.

[15] 陈 红, 彭真明, 王 峻, 等. 地震信号分数阶Gabor变换谱分解方法及应用[J]. 地球物理学报, 2011, 54(3): 867-873.

[16] Wu Xiaoyang, Liu Tianyou. Spectral decomposition of seismic data with reassigned smoothed pseudo Wigner-Ville distribution [J]. Journal of Applied Geophysics, 2009, 68(3): 386-393.

[17] 赵迎月, 顾汉明, 李宗杰, 等. Wigner-Ville高阶时频谱及其在塔中奥陶系缝洞型储层预测中的应用[J]. 石油地球物理勘探, 2010, 45(5): 688-694.

[18] 岳文正, 陶 果. 小波变换在识别储层流体性质中的应用[J]. 地球物理学报, 2003, 46(6): 863-869.

[19] Satish S, Partha S R, Phil D A, et al. Spectral decomposition of seismic data with continuous-wavelet transform [J]. Geophysics, 2005, 70(6): 19-25.

[20] Galiana-Merino J J, Rosa-Herranz J L, Rosa-Cintas S, et al. SeismicWaveTool: continuous and discrete wavelet analysis and filtering for multichannel seismic data [J]. Computer Physics Communications, 2013, 184(1): 162-171.

[21] 朱晓临, 李雪艳, 邢 燕, 等. 基于小波和奇异值分解的图像边缘检测[J]. 图学学报, 2014, 35(4): 563-570.

[22] 王 凤, 万智萍. 基于阈值与人眼特性的小波图像压缩算法[J]. 图学学报, 2013, 34(6): 80-86.

[23] 高静怀, 汪文秉, 朱光明, 等. 地震资料中小波函数的选取研究[J]. 地球物理学报, 1996, 39(3): 392-400.

[24] Marfurt K J. Robust estimates of 3D reflector dip and azimuth [J]. Geophysics, 2006, 71(4): 29-40.

Seismic Horizon Extraction Based on Dip Correction

Liu Xin, Che Xiangjiu, Lin Senqiao

(College of Computer Science and Technology, Jilin University, Changchun Jilin 130012, China)

To solve the problem of incorrect horizons′ extraction caused by the low signal-to-noise ratio of original seismic data sets, we present an anti-noise seismic horizon tracking algorithm, which includes following steps. Firstly, use wavelet transform for de-noising and the calculation of dip angle information. Then, estimate the initial tracking position based on the dip information. Experiments show that horizon tracking results using our method match well with the geological horizons.

horizon extraction; seismic interpretation; coherence technique; wavelet transform

P 631

A

2095-302X(2015)03-0418-07

2014-10-08；定稿日期：2014-10-24

国家自然科学基金资助项目(61170005，60905022)

刘 鑫(1989-)，女，安徽淮北人，硕士研究生。主要研究方向为地震数据处理。E-mail：15104459360@163.com

车翔玖(1969-)，男，吉林吉林人，教授，博士，博士生导师。主要研究方向为医学图像分割、图像传输与信息隐藏、大数据三维可视化及其在地学与医学中的应用。E-mail：chexj@jlu.edu.cn