地球物理学报  2013, Vol. 56 Issue (5): 1592-1608   PDF    
基于地震图像校准的共成像点道集增强技术
李钟晓 , 陆文凯 , 王季     
清华大学自动化系 智能技术与系统国家重点实验室 清华信息科学与技术国家实验室(筹), 北京 100084
摘要: 由共成像点道集抽取的共偏移距道集可以当成相同地下成像的多次观测.由于偏移误差的影响, 在不同共偏移距道集上, 同一采样点存在水平和垂直方向上的错位.本文提出一种基于地震图像校准的共成像点道集增强技术, 实现了不同偏移距道集在时间和空间上的逐点匹配对齐.在本文中以2D局部归一化互相关来表征共偏移距道集和叠加道集在相同时空位置的相似度, 假设在不同水平和垂直移动量时的互相关满足连续凸函数, 利用求导方法估计共偏移距道集在该位置处的非整数校正量, 最后采用双线性内插方法估计成像振幅.传统道集拉平技术在垂直方向进行校正量消除, 本文方法能有效估计共偏移距道集中的水平和垂直校正量, 并在亚像素域估计正确的成像振幅.模型数据和实际数据的处理结果表明, 本文方法能有效增强共成像点道集中同相轴的一致性, 提高叠加结果的分辨率.
关键词: 地震图像校准      水平和垂直校正量      求导方法      双线性内插     
Common image gathers enhancement based on seismic image warping
LI Zhong-Xiao, LU Wen-Kai, WANG Ji     
State Key Laboratory of Intelligent Technology and Systems, Tsinghua National Laboratory for Information Science and Technology, Department of Automation, Tsinghua University, Beijing 100084, China
Abstract: The common offset gathers extracted from common image gathers can be seen as multiple observations about the same underground image. Due to migration errors, the real image positions deviate from the true positions with lateral and vertical misalignments. This paper proposes a seismic image warping technique to enhance the common image gathers. This method aligns different common offset gathers point to point in both the lateral and vertical direction. In this paper the 2D local normalized cross-correlation can be used to measure the similarity of the input gather and the model gather at the same temporal and spatial position. By assuming the 2D cross-correlation computed at various lateral and vertical relative lags satisfies a continuous convex function, the derivation method can be used to estimate the non-integer correction values corresponding to the maximum correlation at the given position. After that the accurate image amplitude can be estimated by the bilinear interpolation. Compared with the traditional gather flattening method, the proposed method can estimate the lateral and vertical correction values of the common offset gathers and obtain the accurate imaging amplitude in the sub-pixel domain. The processing results of the 2D synthetic and real dataset demonstrate the proposed method can enhance the event consistency of common image gathers and improve the resolution of the stack..
Key words: Seismic image warping      Lateral and vertical correction values      Derivation method      Bilinear interpolation     
1 引言

在理想情况下经2D叠前偏移生成的共成像点道集(Common Image Gather,CIG道集)对应于相同的地质地层位置.我们期望CIG道集中的各道保持高度的一致性,然而由于偏移速度误差、各向异性等的影响,地质构造并没有正确的成像.实际成像点在水平和垂直方向上会偏离真实成像点,导致CIG道集中的地震同相轴没有得到有效拉平.由共成像点道集抽取的共偏移距道集(Common Offset Gather,COG道集)可当成相同地下成像的多个观测.COG道集上各实际成像点存在的水平和垂直错位,会降低叠加剖面的分辨率,并导致地震属性分析和地震反演精度的降低.

到目前为止,已经提出许多道集拉平方法[1-13]来消除叠前地震道集(如共中心点道集(Common Middle Point,CMP道集)或CIG道集)中各点的垂直校正量(剩余时差),提高叠加剖面的分辨率.基于互相关的方法[1-9]通过计算待校正地震道和模型地震道的互相关来自适应估计和消除垂直方向的校正量;以相位替换为基础的方法[10-13]是根据信号到达时完全包含在信号相位谱的理论,首先计算模型道的相位谱,然后对待校正道的相位谱进行替换来完成垂直校正量消除.这些方法能有效消除垂直方向各个成像点与真实成像点的错位,但忽略了水平错位对叠加质量的影响.另外,动校正后的CMP道集或叠前偏移后的CIG道集中的局部相关值作为叠加道的加权系数,能有效提高地震数据叠加的信噪比和成像精度[14-15].Perez等将图像校准技术应用在3D叠前偏移数据的时间切片上,来消除叠前偏移数据中的水平校正量[16].但经过高精度偏移速度分析后,CIG道集仍然会存在未拉平的同相轴而产生垂直方向的剩余时差[2].如何有效地同时估计并消除COG道集中存在的垂直和水平校正量,提高CIG道集同相轴的一致性,成为值得研究的问题.

本文利用地震图像校准技术对2D叠前偏移后的COG道集各点进行校正,同时消除COG道集的水平和垂直扰动量,增强CIG道集同相轴的一致性和叠加结果的分辨率.图像校准技术[17]是一种用来进行图像变换的数学工具,被广泛应用在医学成像和计算机图形学等领域[18-19],用来拾取和消除不同图像之间的扰动量.在石油物探领域,图像校准技术用来监测时移油藏变化[20-25]和评价AVO分析中的不确定性[26]等.本文方法将叠加剖面作为模板道集,计算各个COG道集和模板道集的局部2D互相关值,同时假设互相关值随2D相对移动量的变化满足连续凸函数,通过计算互相关函数的最大值来估计COG道集中各点的水平和垂直校正量.为消除计算误差的影响,需要对估计的2D校正量进行加权均值滤波处理.与传统的基于互相关的剩余时差消除方法[1-9]或水平校正量消除方法不同[16],本文在水平和垂直方向计算2D互相关最大值对应的水平和垂直校正量,并在亚像素域估计正确的成像点幅值.基于互相关的自动统计剩余静校正方法应用在野外静校正之后,消除由于地形突变、风化层基底和风化层速度变化所造成的走时畸变.所估计校正量与采样时间无关,只与空间位置有关,每道估计一个时延量并进行剩余静校正量消除[27].本文方法应用在2D叠前偏移之后,消除由于偏移速度误差、偏移算法等造成的水平和垂直错位.所估计的校正量与COG道集各点的空间和时间位置有关.

下面首先介绍图像校准技术的基本原理,然后给出叠前道集水平和垂直校正量消除的算法流程图,最后采用所提方法处理模型数据和实际数据,并与传统基于互相关的垂直校正量[2-4]消除的效果进行对比.

2 基本原理

由于偏移误差的影响,COG道集中实际成像点在水平和垂直方向偏离真实位置,如图 1所示.图像校准技术将待校正道集(共偏移距道集)与模型道集进行匹配,来逐点估计水平和垂直校正量.考虑到信噪比要求,模型道集通常选择为CIG道集中近偏移距到中偏移距的部分叠加[4].由于对COG道集各点估计的水平和垂直校正量往往位于非整数位置,因此需要对COG道集进行重采样估计真实的成像点振幅.

图 1 实际成像点和理想成像点的关系 Fig. 1 The relation ship between the real image point and the ideal image point

在本文方法中[dt(rs),dx(rs)]代表所估计的垂直校正量和水平校正量,即图 1中实际成像点偏离真实成像点的垂直和水平位移,rs分别代表实际成像点的时间和空间位置.如果gb(rs)表示待校正道集,校正后道集ga(rs)可以表示为:

(1)

垂直校正量dt和水平校正量dx与输入道集中成像点的位置rs有关.输入道集和模型道集的局部互相关可被用来衡量两个道集在该点的相似性.图 2描述了地震图像校准的基本原理.a图表示待校正道集(上)和模型道集(下),以待校正道集中的某点为中心开窗(图b),计算待校正道集和模型道集在该点不同移动量处的2D互相关值(图c),搜索最大互相关值对应的水平和垂直校正量,并对待校正道集进行时空变换和重采样(图d).

图 2 地震图像校准示意图[1622] (a)待校正道集(上)模型道集(下);(b)计算2D互相关提取的局部窗口数据;(c)计算的2D互相关(横轴:空间相对移动量纵轴:时间相对移动量),互相关最大值对应的移动量为所要估计的水平和垂直校正量;(d)在整个道集上逐点估计2D校正量,并对待校正道集进行时空变换和重采样. Fig. 2 The diagram of seismic image warping[16, 22] (a)The gather to be corrected (top) and the model gather (bottom); (b)The extracted local window data used for 2D cross-correlation computation; (c) The computed 2D cross-correlation (lateralaxis: spatial relative lag, vertical axis: temporal relative lag). The lateral and vertical corrected values correspond to the maximum cross-correlation; (d) The corrected values are computed point by point in the whole gather. Then the gather to be corrected is transformed in the space and time domain and resampled.

图像校准分为校正量估计和重采样两部分.本文方法首先对COG道集逐点计算水平和垂直校正量,并估计成像点位置,然后对COG道集重采样估计成像点振幅.

2.1 水平和垂直校正量计算

在本文方法中,为消除不同时间和空间相对移动量时窗口能量差异对互相关计算的影响,度量两个道集的相似性时采用2D归一化互相关[28]

(2)

其中,

C(rsl1l2)代表待校正道集g和模板道集f之间的2D归一化互相关值,d0(rsl1l2)代表待校正道集g和模板道集f之间的2D互相关值,d1(rsl1l2)代表模板道集f的自相关值,d2(rsl1l2)代表待校正道集g的自相关值,l1l2表示时间和空间相对移动量,[-pp]、[-qq]分别表示时间和空间相对移动量的范围,[-ww]、[-vv]分别表示归一化互相关窗口的时间和空间范围.

为计算2D局部归一化互相关值,在道集f中包含点(rs)开一个窗口f(t1x1),其中t1=(r-w,…,r-1,rr+1,…,r+w),x1=(s-v,…,s-1,ss+1,…,s+v),在道集g中包含点(rs)开一个窗口g(t2x2),其中,t2=(r-w,…,r-1,rr+1,…,r+w),x2=(s-v,…,s-1,ss+1,…,s+v).起初选择t1=t2x1=x2,然后以不同的2D相对移动量l=(l1l2)改变(t2x2),计算在该相对移动量时的局部归一化互相关值C(rsl1l2).整个计算过程如图 3所示.窗口的选择根据地震子波长度和待校正道集上的同相轴倾角决定[17].

图 3 2D局部互相关计算示意图[17] Fig. 3 The diagram of 2D local cross-correlation computation[17]

传统方法在道集的稀疏点上按图 3所示计算局部归一化互相关值,并通过搜索最大互相关来估计校正量,然后进行插值计算其它点的校正量[29].虽然这种方法提高了计算效率,但降低了道集上其它各点校正量估计的精度.为了均衡互相关和校正量计算的精度和效率,首先计算待校正道集g在相对移动量l=(l1l2)时,与模板道集f在每一个点的乘积,然后对乘积结果进行快速均值滤波[30-31]以计算局部互相关值d0(rsl1l2):

(3)

其中,FMFT表示快速均值滤波,时间和空间窗口大小分别为2*w+1,2*v+1,.*表示点乘运算.同理:

(4)

最后采用公式(2)计算每一点的归一化互相关值C(rsl1l2).整个计算过程如图 4所示.这种归一化互相关计算方法简单易行,并能够避免在对相邻位置点计算局部互相关时,产生的重复相乘和相加运算.重复该过程计算在其它相对移动量的局部归一化互相关值.对各个待校正道集g计算互相关时采用相同的模板道集f,因此道集f的自相关值首先计算并进行保存,以进一步降低各道集归一化互相关的总体计算量.传统方法[2-4]在垂直方向计算1D归一化互相关值,因此本文方法在归一化互相关的计算量大约是传统方法的(2q+1)(2v+1)倍.

图 4 2D归一化互相关计算示意图 Fig. 4 Illustration for computing the 2D normalized cross correlation

在计算不同2D移动量对应的局部互相关值后,遍历搜索局部互相关的最大值,对应的校正量位于整数采样点处.然而,即使小于整数采样点的校正量会产生足够大的地震同相轴扰动[22].为避免地震数据的离散采样特性,可以假设在不同整数相对移动量时的互相关值满足一个连续凸函数[32-34].在此基础上求解互相关函数的最大值来估计非整数的校正量.

Newton-Raphson方法[16, 24]已经被用来估计最大互相关值对应的非整数校正量.我们采用求导方法来估计采样点(rs)处最大互相关对应的校正量,式(2)定义的2D归一化互相关函数在最大值附近区域内可近似满足如下的二元抛物面函数[32-34]

(5)

计算得到整数相对移动量(l1l2)下的归一化互相关后,未知参数向量(abcdef)可以利用多项式拟合方法进行估计.令归一化互相关函数C(rsl1l2)的一阶导数为零,可以得到采样点(rs)处最大互相关对应的时间校正量dt和空间校正量dx如下式:

(6)

图 5a给出了校正量估计的示意图,计算各个整数移动量对应的互相关,搜索所计算互相关最大值对应的整数移动量(左图),并取其周围的整数移动量及其互相关值对多项式进行拟合(中图),然后利用求导方法估计非整数校正量(右图),叉号表示所估计的校正量位置.

图 5 成像点位置和振幅计算示意图[16, 24] (a)不同移动量下的归一化互相关值(左图),用于拟合多项式的互相关值(中图),叉号为求导方法估计的校正量位置(右图);(b)校正前成像点位置(左图)和校正后成像点位置(右图). Fig. 5 The diagram of estimating image point position and amplitude (a) The normalized cross-correlation of different relative lags (left). The normalized cross-correlation chosen for polynomial fitting (middle). The cross corresponds to the corrected values position by the derivation method (right); (b) The image point position before correction (left) and after correction (right).

地震数据的局部连续性决定了所估计的水平和垂直校正量也是局部连续的[16].而地震道集中的噪声会影响校正量计算的准确性,校正量估计的野值或奇异值会导致输出的校正道集产生不合理的突变或不连续现象[16].沿主要地质构造特征移动的校正量会计算得到大的互相关值,因此各点的最大互相关值可用来衡量该点处所估计校正量的相对重要性[16].本文对所估计的水平校正量和垂直校正量分别进行加权均值滤波处理,加权系数选择为各点所对应的最大互相关值.而且,最大互相关对应的时间和空间位置不能偏离原始成像点位置太远,因此采用简单的阈值处理将所估计的时间校正量和空间校正量限制在时间相对移动量[-pp]和空间相对移动量[-qq]的范围之内.

上述所计算的校正量表征待校正道集和输出道集的映射关系.基于这种时间-空间映射关系所估计的成像点位置在大部分情况下位于非整数位点.因此,需要在亚像素域对待校正道集进行重采样.

2.2 基于双线性内插的重采样

为同时消除2D叠前偏移COG道集的水平和垂直校正量,应该对抽取的每一个COG道集逐点进行图像校准.利用式(1)表示的输入和输出道集的映射关系,对待校正道集上的每一点估计正确的成像点位置,其中dt、dx为求导方法计算的时间和空间方向的校正量.图 5b为对COG道集上某点进行校正量消除,并估计正确成像点位置的示意图.图 5b右图中白点表示输出道集中的成像点位置,图 5b左图中白点表示该成像点在待校正道集的位置.

图 1图 5b可知,地震数据的离散采样特性决定了所估计的成像点位置位于亚像素域.Wolberg利用多种内插算子对待校准图像进行重采样[18].在确定亚像素位置周围最近的四个离散采样点后,如图 5b左图中的黑点位置,本文采用双线性内插方法对所估计各点的成像振幅进行重建.局部的内插算子能够尽量减小对地震数据中断层等非连续性的模糊效应,并能提高重采样的计算效率[16].

2.3 算法流程图

综上所述,对叠前地震道集(如偏移后的CIG道集或动校正后的CMP道集)进行水平和垂直校正量消除的算法流程图如图 6.

图 6 基于地震图像校准的水平和垂直校正量消除算法流程图 Fig. 6 The diagram of removing lateral and verical corrected values based on seismic image warping

本文方法计算不同时间和空间相对移动量下的2D归一化互相关,并估计互相关最大值对应的时间和空间校正量.传统方法计算时间相对移动量下的1D互相关最大值来估计时间校正量,当将本文方法表征2D局部归一化互相关的空间窗口范围v设为0,空间相对移动量的范围q设为0,即可得到传统方法[2-4]消除垂直校正量的处理结果.

3 应用实例

为验证本文方法的有效性,对存在水平和垂直校正量的人工模拟数据和实际数据进行处理,测试该方法消除水平和垂直校正量的效果,并与传统方法[2-4]进行垂直校正量消除的结果进行对比.在该应用实例中,[-pp]、[-qq]分别表示时间和空间相对移动量的范围,[-ww]、[-vv]分别表示归一化互相关窗口的时间和空间范围,mn分别表示对整个道集估计的水平和垂直校正量进行加权均值滤波的时间和空间大小.

3.1 合成数据

首先人工合成CIG道集进行测试.选取一个叠加剖面,采样间隔0.002s,道间距12.5m,图 7a所示的叠加剖面有1000道,每道有1000个采样点,对黑色框图内各点随机施加水平和垂直扰动量并进行重采样后的结果如图 7b所示,图 8a为各点包含的水平校正量,在-4道到4道范围内随机选取,图 8c为垂直校正量,在-4个采样点到4个采样点范围内随机选取.图 7c图 7d图 7a图 7b上面白色方框的放大结果,图 7e图 7f图 7a图 7b下面白色方框的放大结果,对比图 7c图 7d图 7e7f箭头位置可发现有明显的水平和垂直校正量存在.对同样区域随机施加水平和垂直扰动量模拟生成50个COG道集,抽取第500个CIG道集如图 9b所示,对比图 9a未施加校正量的CIG道集,可知水平和垂直校正量的存在使CIG道集同相轴的一致性变差.

图 7 模拟的COG道集 (a)黑色框图中不包含水平和垂直校正量的COG道集;(b)黑色框图中包含水平和垂直校正量的COG道集;(c)图(a)上面白色方框的放大结果;(d)图(b)上面白色方框的放大结果;(e)图(a)下面白色方框的放大结果;(f)图(b)下面白色方框的放大结果. Fig. 7 The simulated COG gathers (a) The COG gather without lateral and vertical corrected values in the black rectangle; (b) The gather with lateral and vertical corrected values in the black rectangle; (c) The zoomed area of upper white rectangle in Fig.(a); (d) The zoomed area of upper white rectangle in Fig.(b); (e) The zoomed area of lower white rectangle in Fig.(a); (f) The zoomed area of lower white rectangle in Fig.(b).
图 8图 7b估计的校正量和真实校正量 (a)真实的水平校正量;(b)采用本文方法估计的水平校正量;(c)真实的垂直校正量;(d)采用本文方法估计的垂直校正量;(e)采用传统方法估计的垂直校正量 Fig. 8 The estimated corrected values and true corrected values for Fig.7b (a) True lateral corrected values; (b) The estimated lateral corrected values; (c) True vertical corrected values; (d) The estimated vertical corrected values by the method in this paper; (e) The estimated vertical corrected values by the traditional method.
图 9 模拟的CIG道集 (a)不包含水平和垂直校正量的CIG道集;(b)包含水平和垂直校正量的CIG道集;(c)同时消除水平和垂直校正量的CIG道集;(d)只消除垂直校正量的CIG道集. Fig. 9 The simulated CIGs (a) The CIG without lateral and vertical corrected values; (b) The CIG with lateral and vertical corrected values; (c) The CIG after eliminating the lateral and vertical corrected values simultaneously; (d) The CIG after removing the vertical corrected values only.

采用图像校准的方法,p=4采样点,q=4道,w=2采样点,v=2道,m=5采样点,n=5道,同时消除水平和垂直校正量.图 8b图 8d为对图 7b的COG道集所估计的水平和垂直校正量,逐个COG道集处理后的第500个CIG道集如图 9c所示.同时对比传统方法消除垂直校正量的效果,p=4采样点,q=0道,w=2采样点,v=0道,m=5采样点,n=5道,图 8e为对图 7b所估计的垂直校正量,不估计水平校正量.对比图 8c图 8d图 8e框图可知,本文方法在估计水平校正量的同时,所估计的垂直校正量更准确.图 9d为消除垂直校正量后抽取的第500个CIG道集.对比图 9c图 9d箭头位置表明本文算法能有效地同时消除水平和垂直校正量.

3.2 实际数据

将本文方法应用到实际数据处理中.该数据包含叠前时间偏移生成的1120个CIG道集,生成该CIG道集前已经进行了多次波消除,每个道集包含61道,道间距50m,每道2501个采样点,采样间隔2ms.图 10a为前50个偏移距道集叠加生成的模板道集,图 10b为第10个偏移距道集,只显示250到750道,200到3000ms的剖面,图 10c图 10d分别为图 10a图 10b框图的放大显示,对比箭头可知该偏移距道集相对模板道集存在水平和垂直校正量.

图 10 叠加模板道集和第10个COG道集 (a)模板道集;(b)第10个COG道集;(c)图(a)黑色框图的放大显示;(d)图(b)黑色框图的放大显示. Fig. 10 The stack model gather and the tenth COG (a)The model gather; (b) The tenth COG; (c) The zoomed area of the black rectangle in Fig.(a); (d) The zoomed area of the black rectangle in Fig.(b).

选择处理参数,p=3采样点,q=2道,w=30采样点,v=2道,m=5采样点,n=5道,同时消除水平和垂直校正量.图 11(a,b)为对第10个COG道集同时进行水平和垂直校正量消除,所估计的水平校正量和垂直校正量.对比传统方法进行垂直校正量消除的效果,选择处理参数,p=3采样点,q=0道,w=30采样点,v=0道,m=5采样点,n=5道.图 11c为传统方法对第10个COG道集所估计的垂直校正量,不进行水平校正量估计.下面从CIG道集和叠加剖面上,将本文方法同时进行水平和垂直校正量消除的效果和传统方法进行垂直校正量消除的效果进行对比.

图 11 第10个COG道集估计的校正量 (a)本文方法估计的水平校正量;(b)本文方法估计的垂直校正量;(c)传统方法估计的垂直校正量. Fig. 11 The estimated corrected values of the tenth COG (a) The estimated lateral corrected values by the method in this paper; (b) The estimated vertical corrected values by the method in this paper; (c) The estimated vertical corrected values by the traditional method.

图 12a为处理前第401个CIG道集,图 12b为对401个CIG道集同时进行水平和垂直校正量消除的结果,为对比显示处理效果,将图 12(a、b)方框内的区域进行放大,图 13(a,b)分别对应图 12(a,b) 中框图的放大结果.图 14(a,b)分别为第582个CIG道集处理前后的结果,图 15(a,b)分别为图 14(a,b)中框图的放大结果.对比图 13(a,b)图 15(a,b)框图内的同相轴可知,经过处理后的成像点道集中同相轴的一致性增强.

图 12 第401个CIG道集 (a)处理前;(b)本文方法同时进行水平和垂直校正量消除;(c)传统方法进行垂直校正量消除. Fig. 12 The 401 CIG (a) The result before process; (b) The result after removing the lateral and vertical corrected values with this paper′s method; (c) The result after removing the vertical corrected values with the traditional method.
图 13 图 12a的局部放大图 (a)图 12a中框图的放大结果;(b)图 12b中框图的放大结果;(c)图 12c中框图的放大结果. Fig. 13 The zoomed area of Fig.12 (a) The zoomed area of the rectangle in Fig. 12a; (b) The zoomed area of the rectangle in Fig. 12b; (c) The zoomed area of the rectangle in Fig. 12c.
图 14 第582个CIG道集 (a)处理前;(b)本文方法同时进行水平和垂直校正量消除;(c)传统方法进行垂直校正量消除. Fig. 14 The 582 CIG (a)The result before process; (b) The result after removing the lateral and vertical corrected values with this paper′s method; (c) The result after removing the vertical corrected values with the traditional method.
图 15 图 14的局部放大 图(a)图 14a中框图的放大结果;(b)图 14b中框图的放大结果;(c)图 14c中框图的放大结果. Fig. 15 The zoomed area of Fig.14 (a) The zoomed area of the rectangle in Fig. 14a; (b) The zoomed area of the rectangle in Fig. 14b; (c) The zoomed area of the rectangle in Fig. 14c.

图 12c为传统方法对第401个CIG道集进行垂直校正量消除的结果,图 13c图 12c中框图的放大结果;图 14c为传统方法在第582个CIG道集的处理结果,图 15c图 14c中框图的放大结果.对比图 13b图 13c图 15b图 15c框图内的同相轴可知,同时进行水平和垂直校正量消除,能得到同相轴一致性更好的CIG道集.

下面在叠加剖面上对比校正量消除的效果.图 16a为处理前的叠加道集,图 16b为本文方法同时进行水平和垂直校正量消除的叠加剖面,其中只显示出250到750道,200ms到3000ms的剖面,将其中的部分区域进行放大,图 17(a,b)分别给出了图 16(a,b)右边框图的放大结果,图 17(d,e)图 16(a,b)左边框图的放大结果.黑色箭头表明经水平和垂直校正量消除后叠加剖面的同相轴变清晰,连续性变强.图 16左边竖线表示第401个CIG道集的叠加道位置,右边竖线表示第582个CIG道集的叠加道位置,图 12-15已经给出了第401个和582个CIG道集的处理前后结果,经过水平和垂直校正量消除后的CIG道集同相轴的扰动减小,同相轴的一致性增强,这样所得到的叠加结果的分辨率更高,有利于进一步的地震反演和解释.

图 16 叠加剖面 (a)处理前;(b)本文方法同时进行水平和垂直校正量消除;(c)传统方法进行垂直校正量消除. Fig. 16 The stack profile (a) The result before process; (b) The result after removing the lateral and vertical corrected values with this paper′s method; (c) The result after removing the vertical corrected values with the traditional method.

图 16c为传统方法进行垂直校正量消除的叠加结果,图 17(c,f)分别为图 16c右边框图和左边框图的放大结果.将图 17(b,c,e,f)中的框图部分放大,分别显示在图 18(a,b,c,d).图 18(a,b)图 18(c,d)中的黑色箭头表明,同时进行水平和垂直校正量消除后,叠加结果的同相轴更加清晰连续,分辨率更高,进一步验证了本文算法同时消除水平和垂直校正量的有效性.

图 17 图 16的局部放大图 (a)图 16a右边框图放大结果; (b)图 16b右边框图放大结果; (c)图 16c右边框图放大结果; (d)图 16a左边框图放大结果; (e)图 16b左边框图放大结果; (f)图 16c左边框图放大结果. Fig. 17 The zoomed area of the Fig. 16 (a)The zoomed area of the right rectangle in Fig. 16a; (b)The zoomed area of the right rectangle in Fig. 16b; (c)The zoomed area of the right rectangle in Fig. 16c; (d)The zoomed area of the left rectangle in Fig. 16a; (e)The zoomed area of the left rectangle in Fig. 16b; (f)The zoomed area of the left rectangle in Fig. 16c.
图 18 图 17框图的局部放大图 (a)图 17b框图放大结果; (b)图 17c框图放大结果; (c)图 17e框图放大结果; (d)图 17f框图放大结果. Fig. 18 The zoomed area of the Fig. 17 (a) The zoomed area of the rectangle in Fig. 17b; (b)The zoomed area of the rectangle in Fig. 17c; (c)The zoomed area of the rectangle in Fig. 17e; (d)The zoomed area of the rectangle in Fig. 17f.

在本文方法中时间和空间相对移动量范围过大,易导致所估计的成像点位置偏离真实位置,校准后产生串层现象;时间和空间相对移动量范围过小,会导致较大校正量的估计误差.过大的互相关时间和空间窗口导致校正量估计的分辨率降低,窗口过小,互相关计算易受噪声影响.同样对所估计水平和垂直校正量进行加权均值滤波的参数应试验选取,在滤除噪声和野值影响的同时,尽量保持所估计校正量的局部细节.

4 结论

基于地震图像校准的水平和垂直校正量消除方法应用在经过动校正的CMP道集或叠前偏移后的CIG道集,处理之前需要进行多次波消除.本文方法以叠加剖面为模板,逐个COG道集消除水平和垂直校正量.以不同的水平和垂直相对移动量计算2D归一化互相关函数,采用求导方法估计非整数校正量,并在亚像素域估计正确的成像点振幅.模型数据和实际数据的处理结果表明本文方法可以同时消除水平和垂直方向上的校正量,增强CIG道集中同相轴的一致性,提高叠加剖面的分辨率.

本文方法在COG道集上逐点计算2D归一化互相关函数并估计校正量,可以提高整个道集上校正量消除的准确性,计算量主要集中在互相关计算和校正量估计上.为均衡计算效率和计算精度,本文采用先点乘再进行快速均值滤波的方法,来计算待匹配道集和模板道集的自相关和互相关,最后计算2D归一化互相关.2D互相关计算的窗口大小需要通过根据地震子波长度和同相轴的倾角试验获得.对所估计的校正量进行2D加权均值滤波能有效消除计算误差影响,提高成像点位置和振幅估计的准确性.

参考文献
[1] 夏洪瑞, 陈德刚, 周开明. 剩余动校正量的拾取与消除. 石油地球物理勘探 , 1997, 32(6): 872–877. Xia H R, Chen D G, Zhou K M. Pickup and removal of residual normal moveout. Oil Geophysical Prospecting (in Chinese) , 1997, 32(6): 872-877.
[2] Hinkley D, Bear G W, Dawson C. Prestack gather flattening for AVO. 74th Annual International Meeting, SEG, Expanded Abstracts, 2004:271-273. http://www.oalib.com/references/18986256
[3] Duveneck E, Traub B. Automatic moveout correction by local event correlations on coherency-enhanced gathers. 76th Annual International Meeting, SEG, Expanded Abstracts, 2006:3036-3040. http://www.oalib.com/references/18986257
[4] Gulunay N, Gamar F, Hoeber H, et al. Robust residual gather flattening. 77th Annual International Meeting, SEG, Expanded Abstracts, 2007:229-233.
[5] 慎国强, 王玉梅, 孟宪军, 等. 基于时频分析的地震道集校平技术应用. 中国石油大学学报(自然科学版) , 2010, 34(1): 34–36. Shen G Q, Wang Y M, Meng X J, et al. Application of seismic gather flattening technique based on time-frequency analysis. Journal of China University of Petroleum (Edition of Natural Science) (in Chinese) , 2010, 34(1): 34-36.
[6] Lu W K, Liu J. Estimation and removal of the residual normal moveout based on the S-transform. Congress on Image and Signal Processing , 2008: 573-576.
[7] Gulunay N, Magesan M, Roende H. Gather flattening based on event tracking for each time sample. 70th EAGE Conference and Exhibition, 2008. http://www.oalib.com/references/18986279
[8] Gulunay N, Magesan M, Roende H. Gather flattening. The Leading Edge , 2007, 26(12): 1538-1543. DOI:10.1190/1.2821939
[9] Gulunay N. Polarity blind and polarity sensitive gather flattening methods. 74th EAGE Conference and Exhibition, 2012. http://www.oalib.com/references/18986260
[10] 林伯香, 孙建国. 相位替换法剩余时差校正. 石油物探 , 2001, 40(3): 15–22. Lin B X, Sun J G. Residual moveout correction by using phase replacement. Geophysical Prospecting for Petroleum (in Chinese) , 2001, 40(3): 15-22.
[11] 王开燕, 李慧, 官波, 等. 相位剩余时差校正方法在提高分辨率中的应用. 大庆石油学院学报 , 2007, 31(6): 14–17. Wang K Y, Li H, Guan B, et al. Application of residual moveout correction to the improvement of resolution. Journal of Daqing Petroleum Institute (in Chinese) , 2007, 31(6): 14-17.
[12] 王云专, 杨立伟, 李素华. 剩余时差校正及泊松比反演. 地球物理学进展 , 2006, 21(1): 214–218. Wang Y Z, Yang L W, Li S H. Residual moveout correction and Poisson's ratio inversion. Progress in Geophysics (in Chinese) , 2006, 21(1): 214-218.
[13] Lichman E. Automated phase-based moveout correction. 69th Annual International Meeting, SEG, Expanded Abstracts, 1999:1150-1153. http://www.oalib.com/references/18986264
[14] Liu G C, Fomel S, Jin L, et al. Stacking seismic data using local correlation. Geophysics , 2009, 74(3): V43-V48. DOI:10.1190/1.3085643
[15] Liu G C, Fomel S, Chen X H. Stacking angle-domain common-image gathers for normalization of illumination. Geophysical Prospecting , 2011, 59(2): 244-255. DOI:10.1111/gpr.2011.59.issue-2
[16] Perez G, Marfurt K J. Warping prestack imaged data to improve stack quality and resolution. Geophysics , 2008, 73(2): P1-P7. DOI:10.1190/1.2829986
[17] Wolberg G. Digital image warping. IEEE Computer Society Press, 1990. http://www.oalib.com/references/18986268
[18] Lee S Y, Chwa K Y, Hahn J, et al. Image morphing using deformation techniques. Journal of Visualization and Computer Animation , 1996, 7(1): 3-23. DOI:10.1002/(ISSN)1099-1778
[19] 郑亚琴, 田心. 医学图像配准技术研究进展. 国际生物医学工程杂志 , 2006, 29(2): 88–92. Zheng Y Q, Tian X. Development of medical image registration. International Journal of Biomedical Engineering (in Chinese) , 2006, 29(2): 88-92.
[20] Nickel M, Sonneland L. Non-rigid matching of migrated time-lapse seismic. 69th Annual International Meeting, SEG, Expanded Abstracts, 1999:872-875. http://www.oalib.com/references/18986271
[21] Rickett J E, Lumley D E. Cross-equalization data processing for time-lapse seismic reservoir monitoring:A case study from the Gulf of Mexico. Geophysics , 2001, 66(4): 1015-1205. DOI:10.1190/1.1487049
[22] Hall S A, MacBeth C, Barkved O I, et al. Cross-matching with interpreted warping of 3D streamer and 3D ocean-bottom-cable data at Valhall for time-lapse assessment. Geophysical Prospecting , 2005, 53(2): 283-297. DOI:10.1111/gpr.2005.53.issue-2
[23] Hall S A. A methodology for 7D warping and deformation monitoring using time-lapse seismic data. Geophysics , 2006, 71(4): 021-031.
[24] 郭念民, 吴国忱. 非重复采集时移地震正演模拟及可行性分析. 地球物理学进展 , 2012, 27(1): 232–245. Guo N M, Wu G C. Forward-modeling and feasibility study of non-repeating acquired time-lapse seismic exploration. Progress in Geophys. (in Chinese) , 2012, 27(1): 232-245.
[25] 赵伟. 中国海上时移地震技术应用的可行性研究. 勘探地球物理进展 , 2003, 26(1): 30–34. Zhao W. Feasibility study on time-lapse seismic offshore China. Progress in Exploration Geophysics (in Chinese) , 2003, 26(1): 30-34.
[26] Grubb H, Tura A, Hanitzsch C. Estimating and interpreting velocity uncertainty in migrated images and AVO attributes. Geophysics , 2001, 66(4): 1208-1216. DOI:10.1190/1.1487067
[27] 李振春, 张军华. 地震数据处理方法. 东营: 中国石油大学出版社, 2006 . Li Z C, Zhang J H. The Method of Seismic Data Processing (in Chinese). Dongying: China University of Petroleum Press, 2006 .
[28] Barnea D I, Silverman H F. A class of algorithms for fast digital image registration. IEEE Trans. Computers , 1972, 21(2): 179-186.
[29] Perez G, Marfurt K J. Fine-tuning your seismic image:prestack data warping to improve stack quality and resolution. 76th Annual International Meeting, SEG, Expanded Abstracts , 2006: 2475-2479.
[30] Pan J J, Tang Y Y, Pan B C. The algorithm of fast mean filtering. Wavelet Analysis and Pattern Recognition Conference , 2007: 244-248.
[31] Rakshit S, Ghosh A, Shankar B U. Fast mean filtering technique (FMFT). Pattern Recognition , 2007, 40(3): 890-897. DOI:10.1016/j.patcog.2006.02.008
[32] Steger C. An unbiased detector of curvilinear structures. IEEE Transactions on Pattern Analysis and Machine Intelligence , 1998, 20(2): 113-125. DOI:10.1109/34.659930
[33] Steger C. Unbiased extraction of curvilinear structures from 2D and 3D images. München:Technical University of München, 1998. http://www.oalib.com/references/18986287
[34] 雷鸣, 张广军. 基于互相关的图像匹配亚像素定位. 光电工程 , 2008, 35(5): 108–113. Lei M, Zhang G J. Image orientation algorithm with subpixel accuracy based on correlative matching method. Opto-Electronic Engineering (in Chinese) , 2008, 35(5): 108-113.