刘旭跃， 周 巍， 张 兵， 黄 骏

(中国石化石油 物探技术研究院，南京 211103)

一种基于图像学的地震层位自动追踪方法

刘旭跃， 周 巍， 张 兵， 黄 骏

(中国石化石油 物探技术研究院，南京 211103)

构造复杂地区地震数据受信噪比低、噪声干扰等影响，层位追踪结果精度不高，为解决该问题，提出了一种适用于低信噪比，同相轴不稳定的基于图像学的地震层位自动追踪方法。该方法首先采用图像学的全变分方法对地震数据进行保护边缘的降噪处理，然后采用基于相似系数的相干图像处理技术对地震剖面进行计算获得相干剖面，同时提取出构造的倾角，接着对同相轴进行自动追踪，程序中采用链式结构存储数据，便于检通过实际数据测试，该方法能有效去除噪声，增强同相轴连续性，提高自动追踪效果。

全变分; 相干图像处理; 同相轴倾角; 层位自动追踪

0 引言

地震野外获得的原始资料，经过室内处理后，得到可供解释的地震剖面和其他成果图件。对这些图件进行分析研究，可以推断出地震剖面上各反射层所相当的地质层位，从而了解地下地质的情况。其中层位地识别和追踪是地震层位解释的重要环节，它对后续的处理和解释工作有很大影响，是地震层位反射同相轴空间分布地直观体现[1]。

目前常规应用的层位追踪方法主要包含：①人工拾取层位；②计算机自动追踪[3]。人工拾取层位是在二维地震剖面上，利用波形相似性对地层连续反射同相轴进行手动追踪得到层位线[4]。它的工作量大，效率低，准确性差。随着计算机和勘探技术的发展，层位自动追踪技术得到了快速发展，它通过计算机计算，自动拾取出同相轴的位置。随着自动追踪效率的不断提高，逐渐取代了人工拾取层位。国、内外对层位追踪的方法主要分为两类:①基于波形特征的自动追踪；②基于相关的自动追踪[5]。波形特征追踪方法只寻找搜索时窗内特征点的像素波形结构形态，不考虑道与道之间的相关性，准确度不高，(边缘检测法[6])。基于相关的自动追踪是根据在搜索时窗内，经过随机过程的互相关分析，计算相邻地震道的互相关来反映同相轴的不连续性，又称为相干图像处理技术，主要应用于描述地震数据的空间连续性，检测相干性变差位置处的异常地质体(如断层、河道等)。Bahorich, M.等[7]提出了第一代相干性计算方法即互相关法，其缺点是利用数据较少，分辨率不高且对噪音的压制不好；Marfurt等[8]提出了第二代相干性计算方法，相似系数法，利用了多道的信息进行相干性的计算，提高了相干剖面的分辨率，而且对噪音的压制有着明显的增强；Gersztenkorn等[9]提出了第三代相干性计算方法即本征结构法，该方法在稳定性方面以及噪音压制方面均优于前两代相干性算法，但其计算量也远远大于前两代算法。此外，还有其他一些算法，如局部结构熵、高阶统计量等。利用相干体技术检测同相轴、提取同相轴倾角的原理是利用同相轴的横向连续性较好，相干性较强；而地质异常体、噪音的横向连续性较差，相干性较弱。同相轴在相干计算后的相干剖面上更为明显。同时由于在计算相干性的过程中对倾角进行扫描，使用与界面倾角相一致的倾角进行相干性计算时得到的相干性最大，因此可以获得界面的倾角，并且能够用于同相轴追踪等后续工作。

为了提高复杂地区低信噪比资料层位自动追踪的精确度，这里提出了一种基于图像学的地震层位自动追踪方法。图像学去噪的方法很多，常见有小波变换[9]、全变分方法等。我们先应用图像学的全变分去噪技术压制噪声，然后基于相似系数的相干性算法对倾角进行扫描，把最大相干值对应的倾角做为界面倾角，最后进行同相轴追踪。

1 方法原理

我们主要应用了图像学中的全变分(TV)去噪技术和基于相似系数的相干成像技术。

1.1 全变分去噪技术

全变分去噪方法可以在保护边缘信息的同时，降低图像的噪声。考虑一个加噪的模型，令μ为原始的清晰图像，μ0为被噪声污染的图像，即：

μ0(x,y)=μ(x,y)+n(x,y)

(1)

其中：n是具有零均值；方差为σ2的随机噪声；如果Ω表示图像的定义域，像素点(x，y)∈Ω。

设Ω是Rn中的有界开子集，μ为局部可积函数，则其全变分定义为式(2)。

(2)

定义的全变分去噪能量泛函为式(3)。

(3)

通常有噪声的全变分比没有噪声的全变分明显大，最小化全变分可以消除噪声，因此基于全变分的降噪可以归结为如下最小化问题:

(4)

满足约束条件：

(5)

(6)

如果一味地对图像的全变分进行极小化，那么代表细节与纹理的许多自身特征也会被一并抹掉。为此，ROF模型用于图像降噪的出发点是最小化“能量泛函”，即

(7)

它的前一项要求输出μ的全变分尽可能小，称为平滑项；后一项则要求μ与μ0尽可能相近，称为数据保真项，它主要起保留图像特征和降低图像失真度的作用。参数λ用来平衡这两个相互冲突的要求。

由梯度下降法，可以得到TV平滑模型：

(8)

在全变分中引入一个小的正数β，对它进行正则化，即用

(9)

(10)

最终的正则化模型为式(11)。

(11)

求解 公式(7)的偏微分方程，利用边界约束条件推出：

(12)

(13)

=>

λ(u-u0)-Δu=0

(14)

离散迭代格式为：

(15)

正则化参数λ的取值非常关键，算法步骤如下：

1)读入带有噪声的数据u0。

2)初始化参数：n=0，Δt=0.25，u0=u0，divp=0。

3)当n小于最大迭代次数时，循环执行如下操作：n=n+1，根据离散迭代公式计算下一步的un,计算扩散项的值。

4)结束迭代，最后一次为去噪的结果。

1.2 相似系数的相干性计算

Chopra, S.等[10]提出基于相似系数的相干性计算，可以利用与中心道相邻的多道信息进行计算，如图1所示。

图1 相似系数法参与计算的地震道Fig.1 The similarity coefficient method is involved in the calculation of seismic trace

在进行相干性计算时，首先定义一个分析窗口(椭圆窗或者矩形窗)，在分析窗之内以分析点为中心包含J个地震道，如图2所示。方位角Φ和长短轴a、b决定着分析窗的方向和大小，当Φ0、Φa、Φ、a、b确定后分析窗内包含的J的地震道就确定了相干时窗的大小一般根据地震反射波的视周期T而定，通常取T/2～T/3。当计算的相干时窗小于T/2时，不能描述一个完整的波峰或波谷，噪音与层位的影响突出；当计算时窗大于T/2时，包含了多个反射波同相轴，可能会降低分辨率。就参加运算的道数来说，从断层成像清晰度和随机噪声压制程度看，一般参与相干计算的道数越多，平均效应越大，对断层的分辨率越低，这时突出的主要是大断层；相反，相干道数少，平均效应小，就会提高对地层边界、断层、特别是突出了对小断层的分辨率。

图2 相似系数法的计算窗口选取[3]Fig.2 The calculation window selection of similarity coefficient method(a)Φ0 Inine方向的方位角;(b)Φa 分析窗长轴的方位角

(16)

其中：uj为分析窗内第j道；xj、yj分别为该道距离中心道的距离；p、q分别表示x、y方向的倾角。图2是当p、q为0时的情况，当p、q分别为p1、q1时，如图3所示，p1、q1分别表示分析窗沿x、y方向的倾角。

Cs(t,p,q)=

(17)

用不同的倾角进行扫描可以获得使Cs(t,p,q)最大的p和q，即可以获得同相轴的倾角，然后可以进行同相轴追踪等后续工作。

在二维情况下，相干性Cs可以表示为

(18)

基于相似系数的相干性计算利用了多道的信息，考虑同相轴整体的连续性，因此,分辨率优于互相关算法，并且对噪音的压制也较互相关算法好。同时在计算相干性的时候,可以根据文献[11]计算得到同相轴上每一点的倾角，因此我们采用这种方法进行倾角的拾取与同相轴的追踪。

2 自动追踪技术流程及数据结构

自动追踪技术实现流程如图3所示。

在相干计算后，得到的构造倾角可以进行同相轴自动追踪。通过计算机确定拾取层位起始点坐标，自动计算出参考道的时间深度h，检索出该点相干性最大的对应倾角，并记录。然后以下一道为参考道，将时延t赋给h,循环计算，直至计算到相干时窗的最后一道。实现的追踪流程为图4所示。

计算得到的一个同相轴的追踪结果采用链式存储方式保存，其优点是存储量较小，相对顺序存储方式而言，它的插入和删除元素速度较快，方便编辑追踪结果。而且它没有空间限制，存储元素的个数无上限，基本只与内存空间大小有关。给每一个链条分配一个编号用以区别剖面上其他同相轴追踪结果，方便实现层位追踪结果检索的功能。Patel等[12]提出将识别层位转化为网格模型的方法，分离错误层位连接。具体数据结构如图5所示，每条链表存储一条同相轴，包含一个头结点和若干子节点。头结点中存储同相轴的编号、节点个数等，每一个节点存储每一个追踪到的该同相轴上的点坐标和局部倾角。最后所有链表的头结点组成一个头结点数组，便于检索，数据结构如图5所示。

3 应用实例

在基于linux平台的中石化iCluster系统中开发了层位自动追踪功能，进行实际数据测试。选用塔河资料为例，主测线范围在450～1100，CDP范围450～1500，工区内地层较为清晰。

3.1 降噪效果分析

这里采用峰值信噪比(PSNR)[12]和EPI进行评价。PSNR值大时，降噪处理效果就好。EPI是边缘细节保留情况的主要评价标准。

图3 倾角为p1、q1时的计算窗口选取Fig.3 The calculation window selection of the inclination of p1 and q1(a)椭圆形分析窗的情况;(b)矩形分析窗的情况

图4 技术实现流程Fig.4 Technical realization process

图5 数据结构设计Fig.5 Design of data structure

(19)

(20)

图6是原始数据的部分剖面图，图7是采用全变分方法，λ=1.0时经处理得到的部分剖面图。处理后的地震剖面同相轴连续性更好，层间结构也更清晰。

表1 降噪效果Tab.1 Noise reduction effect

图6 原始的剖面图Fig.6 The original profile

图7 全变分处理后的剖面图Fig.7 Profile of total variation after treatmene

3.2 层位追踪结果

都采用基于相似系数的相干计算，分两种情况进行对比：①是不经过全变分降噪处理;②是经过全变分降噪处理，图7是不经过全变分降噪处理的追踪结果，绿色是追踪轨迹。图8是经过全变分降噪处理的追踪结果，绿色是追踪轨迹。比较图8、图9发现，不用全变分处理的追踪结果，跳动大，局部有异常，经过全变分降噪处理后追踪结果有改善，符合较好。

图8 不经过全变分处理的追踪结果Fig.8 Tracking results without total variation

图9 经过全变分处理的追踪结果Fig.9 The tracking results of the total variation

4 结论

全变分方法较好地实现了去噪和保护边缘的统一，保护了同相轴。基于相似系数的相干算法分辨率高，计算后得到的同相轴倾角可用于后续的同相轴自动追踪。该方法在实际资料测试中得到光滑连续的同相轴，与地质层位相符，取得良好的效果。

[1] 谢辉.基于地震层位趋势约束的小层系地质层位随机模拟[J].科技创新与应用,2013.31:83-85.XH.Stochasticsimulationofgeologicalhorizonofsmallstratabasedonthetrendofseismichorizon[J].ScienceandtechnologyinnovationandApplication, 2013,31:83-85.(InChinese)

[2]HOYESJ,CHERETT．Areviewofglobalinterpretationmethodsforautomated3Dhorizonpicking[J].TheLeadingEdge，2011,30(1):38-47.

[3] 彭文,熊晓军,韩小俊.基于高阶统计量的层位自动追踪方法[J].新疆石油地质，2006,6(27):743-745.PENGW,XIONGXJ,HANXJ.Horizonautomatictrackingmethodbasedonhigherorderstatistics[J].XinjiangPetroleumGeology,2006,6(27):743-745.(InChinese)

[4] 郑公营,曾婷婷.地震层位自动追踪技术研究[J].物探化探计算技术,2013,35(6):711-716.ZHENGGY,ZENGTT.Researchonautomatictrackingtechnologyofseismichorizon[J].Geophysicalprospectingandgeochemicalexploration, 2013,35(6):711-716.(InChinese)

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

[6]BAHORICHMS,FARMERSL.3Dseismicdiscontinuityforfaultsandstratigraphicfeatures[J].TheLeadingEdge,1995,14(10):1053-1058.

[7]MARFURTKJ,KIRLINR,FARMERSL,etal.3DSeismicAttributesUsingaSemblance-basedCoherencyAlgorithm[J].Geophysics,1998,63(4):1150-1165.

[8]GERSZTENKORNADAM,MarfurtKJ.Eigenstructure-basedCoherenceComputationsasanAidto3DStructuralandStratigraphicMapping[J].Geophysics,1999,64(5): 1468-1479.

[9]GALIANA-MERINOJJ，ROSA-HERRANZJL，RosaCintasS,eta1.SeismicWaveTool：continuousanddiscretewaveletanalysisandfilteringformultichannelseismicdata[J].ComputerPhysicsCommunications,2013,184(1):162-171.

[10]CHOPRAS,MARFURTKJ.Seismicattributeexpressionofdifferentialcompaction[J].TheLeadingEdge, 2012, 31(31): 1418-1422.

[11]MARFURTKJ．Robustestimatesof3Dreflectordipandazimuth[J].Geophysics, 2006,71(4)：29-40.

[12]PATELD，BRUCKNERS，VIOLAI,eta1.Seismicvolumevisualizationforhorizonextraction[C].PacificVisualizationSymposium(PacificVis),2010IEEE.IEEE,Taipei,2010:73-80.

[13]公成敏. 全变分原理在地震数据去噪中的应用[J].计算机与数字工程,2014,297(42):1271-1274.GONGCM.Thevariationalprincipleanditsapplicationindenoisingofseismicdata[J].Compuuteranddigitalengineering,2014,297(42):1271-1274.(InChinese)

An automatic tracking method for seismic horizons based on image theory

LIU Xuyue, ZHOU Wei, ZHANG Bing, HUANG Jun

(Sinopec Geophysical Research Institute,Nanjing 211103,China)

Seismic data in complex area is affected by low signal-to-noise ratio, noise interference, horizon tracking results are not high, this paper proposes a method to solve the problem. The seismic horizon automatic tracking method based on image theory, the whole variation method to calculate the seismic profile, and then the coherent image processing technology to calculate the seismic profile. And then the phase axis automatic tracking. It is proved that the method can achieve good results using the actual data test.

total variation; coherent image processing; the angle of the same phase; automatic tracing of horizons

2015-12-28 改回日期：2016-03-01

中石化科技部项目(P14152)

刘旭跃(1985-)，男，工程师，现从事物探软件工作， E-mail：Liuxy@sinopec.com。

1001-1749(2017)01-0064-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.01.10