地球物理学报  2018, Vol. 61 Issue (8): 3297-3309   PDF    
三维共偏移距叠前地震数据反假频射线束形成偏移成像
吴帮玉1, 杨辉2, 孙文博3, 姜秀娣3, 高静怀4,5     
1. 西安交通大学数学与统计学院, 西安 710049;
2. 南方科技大学地球与空间科学系, 深圳 518055;
3. 中海油研究总院有限责任公司, 北京 100028;
4. 西安交通大学电子与信息工程学院波动与信息研究所, 西安 710049;
5. 海洋石油勘探国家工程实验室, 西安 710049
摘要:对稀疏/非规则采样或者低信噪比数据,射线束提取困难并伴随有假频产生,对叠加剖面和道集造成严重干扰.为了提升射线束偏移在稀疏和低信噪比地震数据采集中的成像效果,本文提出基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法.射线束偏移首先将地震数据划分为超道集,经过部分NMO后转化为以射线束中心定义的共偏移距数据,倾斜叠加和反假频操作均在局部共中心点坐标上实现.时间域倾斜叠加是对地震数据的时移累加操作,三角低通滤波同样可以在时间域完成,在对地震数据进行因果和反因果积分后,亦为地震数据的时移累加.因此,三角低通滤波与倾斜叠加可在时间域结合同时完成,避免了频域滤波的正反傅里叶变换.本文在反假频公式中加入权重系数,用以对反假频的程度进行控制,达到分辨率和噪声压制的最佳折衷.以某海上三维实际数据为例,文中展示了反假频射线束形成对偏移叠加剖面和共成像点偏移距道集中的噪声进行了有效压制.
关键词: 反假频      三角滤波      局部倾斜叠加      波束形成      三维叠前深度偏移     
Common offset 3D prestack seismic data anti-aliasing beamforming, migration and imaging
WU BangYu1, YANG Hui2, SUN WenBo3, JIANG XiuDi3, GAO JingHuai4,5     
1. School of Mathematics and Statistics, Xi'an Jiaotong University, Xi'an 710049, China;
2. Department of Earth and Space Sciences, Southern University of Science and Technology, Shenzhen 518055, China;
3. CNOOC Research Institute Co., Ltd. Beijing 100028, China;
4. Institute of Wave and Information, School of Electronic and Information Engineering, Xi'an Jiaotong University, Xi'an 710049, China;
5. National Engineering Laboratory for Offshore Oil Exploration, Xi'an 710049, China
Abstract: We propose triangle filter anti-aliasing beam forming to attenuate the migration operator aliasing for the three dimensional beam migration on the real sparse acquisition marine seismic data. Slant stack is an effective way for beam forming, which picks up the local linear seismic events in the seismic profile. The ray parameters in the slant stack is an approximate of the migration operator dip and anti-aliasing is a low-pass frequency filtering operation according to the sampling theory. The frequency band of the filter is determined by the ray parameters of the migration operator dip. Time domain slant stack is the summation of the time shifted seismic traces. The triangle filtering can also be conducted in the time domain. It is also summation of the time shifted data after casual and anti-causal integral. Therefore, triangle filtering and slant stack can be incorporated in the time domain which avoids the transform of data from time to frequency domain. We propose to add weighting factors in the anti-aliasing formula to control the trade off between noise suppression and the resolution decreasing of the steep structures in the image. The three dimension stack image and the common image point offset gathers from the real seismic data show that the proposed method suppresses the operator aliasing effectively.
Key words: Anti-aliasing    Triangle filtering    Local slant stack    Beam forming    3D prestack depth migration    
0 引言

离散数字信号处理是现代地震数据处理的基础.而采样定理,又称香农采样定律或奈奎斯特采样定律又是连续信号离散化的基本依据,是实际连续信号和可供计算机处理的离散数字信号之间的桥梁.采样定理指出,对于一维时间信号,如果带宽小于采样频率的一半(即奎斯特频率),则采样信号可通过重建公式完全表示原信号(Robinson and Treitel, 2008).当采样不满足采样定理时,高于或者处于奈奎斯特频率的部分会导致混叠现象,这部分混入有效频带内的频率成分即称为假频.一般可采取提高采样频率或者引入低通滤波器两种措施避免混叠的发生.对于各维相互独立的多维信号,采样定理以及假频可以参照一维信号在各自维度上单独处理.地震信号处理所面对的是离散的多维数据,尤其是三维情况下,假频是必须面对的问题之一,某些情况甚至直接决定能否继续进行后续的解释或者属性分析等流程.

三维叠前数据本质上为五维信号,即记录时间(t),二维炮点位置坐标(xs, ys)和二维检波点坐标(xr, yr).同样,三维偏移成像结果也可为五维数据,如空间位置坐标(深度(z)和二维水平坐标(x, y))以及地下层位的反射角(θ)和方位角(ϕ).因此,三维叠前偏移成像是多维数据到多维数据的映射,而且每个维度上均为离散量,会涉及多个采样间隔.首先,是观测系统间隔,这其中又涉及到炮点网格和检波器网格,分别在二维空间上采样;其次是地震成像网格,即输出地震图像三维空间位置采样以及反射角和方位角采样间隔.理论上,每个维度上的离散化均需分别满足采样定理.然而,由于受采集成本或者施工限制影响以及综合考虑数据存储量、计算量等因素,无论是数据采集还是后续数据处理环节,均在尽量大的网格上进行操作,不能够总是满足采样定理的要求,即会引入假频.一般称此类由采样引起的假频为数据假频.内插是压制数据假频的有效方法(刘喜武等,2004刘财等,2013),可以解决观测系统不规则及由偏移引起的部分假频问题.但插值重建涉及较大计算以及存储量(刘礼农等,2006),数据量的增大直接降低后续数据处理尤其是偏移成像的计算效率.

数据假频是引起偏移成像中假频噪声的一个原因,另一原因是由叠前偏移操作引起的,称为算子假频(Zhou et al., 2013; Nimsaila et al., 2015).地面接收地震数据是地下介质反射或者散射波场的采样,各种积分类偏移即可视为在叠前五维地震数据中有选择性的“采样累加”,如Kirchhoff偏移(Schneider, 1978)和各类射线束偏移(Gray et al., 2009).图 1所示为二维叠前积分偏移示意图,红色曲线示意偏移算子的积分曲线,地震道等间隔采样Δx.左侧标注的ΔTn(n=1, 2, 3, 4, 5)为偏移算子在各个相邻道之间等效时间采样间隔,可以看出该间隔量不仅与Δx有关,还与偏移算子的斜率相关.若使偏移算子累加量无假频产生,则在偏移曲线上不同道处数据累加所能满足采样定理的频率需满足如下关系(Lumley et al., 1994; Abma et al., 1999):

(1)

图 1 积分法偏移算子假频产生示意图 Fig. 1 Operator aliasing generation for the integral migration method

该表达式中为偏移算子斜率,ΔT为偏移算子在任意相邻道之间的相对采样间隔.由公式(1)可知,无论地震道采样间隔Δx多小,即使时间为连续记录,假频在偏移算子斜率增大到一定程度后仍会出现.在成像结果中,由偏移假频引起的噪声一般表现为随机背景噪声,但有时也会在深部表现为类似地层结构的相干信号(Abma et al., 1999; Claerbout, 2010).公式(1)所示为二维观测系统满覆盖且采样均匀的通用反假频准则,对于三维情况不同研究人员提出了不同的应用公式(Lumley et al., 1994; Claerbout, 1997; Abma et al., 1999Zhang et al., 1999; Biondi, 2001; Gray, 2013).

算子反假频在Kirchhoff偏移中得到了较为充分的发展.最简单高效的假频压制方法即是限制偏移孔径(王华忠等,2010)或者直接限制偏移算子的大角度分量,该策略对小角度反射层结构成像改善效果明显,但并不能完全消除偏移假频,而且对成像陡倾角结构损伤巨大(Abma et al., 1999).原则上,利用公式(1)反假频需要计算走时对空间的导数,该操作在偏移的最内层循环,会严重降低偏移的计算效率(Zhang et al., 2001).为了计算简单以及提高偏移效率,可以用常速(地表速度)来代替.虽然这种近似操作对高倾角结构成像的分辨率和振幅保真度有所下降(Zhou et al., 2013; Gray, 2013),但由于其实用高效,仍然在Kirchhoff偏移反假频中得到广泛应用.Gray(1992)提出与Kirchhoff偏移算子倾角相关的低通滤波反假频算法,该方法利用不同频带的低通滤波器对地震道滤波,产生多个不同带宽的地震数据,偏移时根据算子倾角选择不同数据作为输入.为了实现输入数据在偏移剖面上的平滑过渡,Gray的方法对每道数据需要产生一定数量的截止频率不同的版本,增加了数据预处理计算量和存储量.Claerbout(1992)Lumley等(1994)利用N点三角滤波器,提出了一种实时计算的偏移算子反假频Kirchhoff偏移,以增加一定计算量的代价避免了多频带数据的产生和存储.王华忠等(2010)给出了三维Kirchhoff共偏移距叠前时间偏移中的反假频准则.针对上述反假频中近似带来的误差,Pica(1996)Biondi(2001)以及Zhou等(2013)提出了数据倾角自适应反假频方法,对目标倾角结构成像分辨率和振幅保真均有提升,但牺牲了非目标倾角结构成像质量.Nimsaila等(2015)提出在频率波数域(f-k)或者频率慢度域(f-p)进行炮域数据反假频滤波,用来压制高频逆时偏移成像中的假频噪声.

射线束成像在计算效率、陡倾角成像、信噪比、算法稳定性以及灵活性等方面均有一定优势.该方法将邻近道数据同时向下偏移,对宽方位、宽频带和高密度地震数据采集偏移成像优势明显.尤其是共偏移距射线束偏移成像,地震信号线性特征强,利于射线束形成以及输出偏移距道集.波束形成是射线束偏移中关键步骤(Wu et al., 2015),利用各种局域化分解技术将地震数据分解为不同方向的局部平面波(波束).提取出来的局部平面波是叠前地震数据的一类特征波场(王华忠等, 2015, 2017),其结果直接影响偏移效率和质量.Zhang等(2001)Gray(2013)给出了共炮点波束形成的反假频一般准则.本文旨在针对共偏移距射线束偏移方法,提出三角滤波反假频波束形成算法,并将其应用于Kirchhoff射线束偏移(Sun et al., 2000).后文首先给出射线束偏移对数据的划分,即形成超道集数据,并将地震数据位置表示为超道集中心坐标以及局部坐标,然后将局部源-检波器坐标转换为局部偏移距中心点坐标;其次,通过部分NMO(Normal Moveout),将地震数据转化为以中心点定义的共偏移距数据;然后,给出局部中心点坐标局部倾斜叠加反假频公式.最后,为了适应观测系统有限孔径或采集引起的非规则采样,本文提出在反假频公式中加入权重以对反假频程度进行控制,并以三维海上地震数据为例,展示其应用效果.

1 超道集:束中心坐标+局部相对坐标

射线束偏移将一定数量的相邻道数据同时偏移,这些同时向下传播的地震道组成的集合即可视为超道集(supergather).将多大范围内的相邻数据作为一个超道集来处理,直接影响射线束偏移的效率和质量,是射线束偏移中最为关键的参数之一(Wu et al., 2015).倾斜叠加可以用来提取超道集地震剖面中的局部线性信号,所提取的线性信号即为局部平面波,该平面波的方向为波束地表初始出射方向.若三维叠前地震反射数据表示为u(xs, ys, xr, yr, t),则形成超道集后,地震道位置可写为相对于超道集中心位置坐标的局部坐标

(2)

为超道集中心位置坐标,其中Δxs= 是每一道数据相对于超道集中心位置的坐标,即为局部相对坐标.图 2所示为本文测试所用海上地震数据某一超道集源位置和检波器位置分布示意图.该超道集中共有367道数据,源和检波器位置分别如图 2a2b所示.该数据观测系统crossline和inline采样间隔分别为25 m和12.5 m,由图可以看出源点和检波点位置均未重整化至观测系统采样点上.波束形成可将此类非规则采样数据转化为以超道集中心位置和射线参数规则采样数据,并不影响射线束偏移的应用.由图 2可以看出源的inline和crossline以及检波器的inline采样都较检波器crossline采样稀疏,是产生假频的主要来源.为了避免超道集的划分在偏移过程中引入“脚印”,以及进一步提高倾斜叠加的精度和信噪比,一般划分时超道集之间会有重叠部分,也即有些地震道会属于不同的多个超道集.

图 2 海上地震数据某超道集中源、检波点位置示意图,该超道集中共有367条数据 (a)源点位置;(b)检波点位置. Fig. 2 Demonstration of source and receiver locations from one supergather with 367 traces of the marine seismic data (a) Source locations; (b) Receiver locations.

后文为表示简单,在不引起混淆的情况下,将地震数据表示为,则对其中源和检波点均做倾斜叠加的时间域公式为

(3)

其中(psx, psy)和(prx, pry)分别表示源和检波器射线参数,τ表示倾斜叠加在时间轴上的截距.射线束方法对每个超道集偏移时,可计算以超道集中心位置并以出射角发出射线的走时,地下每个位置的振幅可由抵达该位置处的超道集中心射线参数在数据中内插得到.公式(3)是对超道集的完全倾斜叠加分解,在数据采集较为稀疏时,倾斜叠加后数据量会以数量级的速度急剧增加.如本次测试所用数据,若将超道集中源和检波器的crossline和inline长度均设置为100 m,则超道集中平均道个数为21.2个,仅占观测系统网格数的2.07%.若源和检波器crossline和inline射线角度间隔设置为10°,倾斜叠加角度范围为-50°到+50°,则每个射线参数变量包含11个角度取值,倾斜叠加后每个超道集变为114条数据,增加了690.6倍.为了提高多次重复偏移效率(如速度更新),倾斜叠加的数据可以存储起来以供重复使用,双边倾斜叠加对数据的增加限制了其在大规模数据上的应用.Wu等(2014)利用最小二乘反演方法进行高分辨率波束形成以减少有效射线角度个数,虽然倾斜叠加后数据量大量减小,但对深层弱反射结构成像效果不佳.波束形成可以如公式(3)对源和检波器操作(双边波束形成),也可以对共检波点波场或者共源点波场进行操作(单边波束形成),如各种共炮域或者共检波器波束偏移方法.单边波束形成后的数据通过一系列压缩或者稀疏表示(如匹配追踪),基本能够保持与原始数据量相当或者略有减小.选择共点源或者是共检波器源进行射线束偏移,亦是抑制假频噪声的一种有效策略(Casasanta and Gray, 2015).

积分类偏移方法相较于波动方程类偏移优势之一即是可以在更为灵活的地震数据集进行偏移,如基于Kirchhoff法共偏移距数据的偏移能够方便地输出共偏移距成像点道集(common image offset gather)而无需额外计算量.形成超道集后的地震数据由超道集中心坐标和局部相对坐标表示,可以通过如下公式将局部源点-检波器坐标转换为局部中心点局部偏移距坐标:

(4a)

(4b)

(4c)

(4d)

其中Δhx和Δhy表示局部偏移距,Δmx和Δmy表示局部中心点坐标.经过坐标转换后的超道集 可以表示为 ,则倾斜叠加公式为

(5)

其中的(phx, phy)和(pmx, pmy)分别为偏移距和中心点射线参数,该组参数与源-检波器坐标系下射线参数关系如下:

(6a)

(6b)

(6c)

(6d)

式(5)也是双边波束形成公式,其偏移过程与双边源-检波器波束偏移类似,也存在稀疏采样下射线参数过多而引入数据量增大问题.在偏移距-中心点坐标系下,单边波束形成可针对偏移距坐标或者中心点坐标叠加,即共中心点波场或者共偏移距波场波束形成.共偏移距数据类似一个三维叠后剖面,其中反射信息连续性和相关性好,更易于倾斜叠加提取其中线性信号,因此本文在共偏移距域实现波束形成和偏移.在共偏移距波束形成之前,需将超道集中的局部偏移距(Δhx, Δhy)消除,即将所有道偏移距校正至以超道集中心定义的偏移距上,需要对超道集进行部分NMO.

2 部分正常时差校正

正常时差(NMO)校正是对地震道数据在时间轴上的时移和拉伸,使所有道记录均等效为零偏移距接收数据.对超道集数据进行部分NMO校正公式为(Wu et al., 2014)

(7)

该公式表示需将超道集地震道(Δhx, Δhy)处thx, Δhy)时刻波场平移至,使该地震道数据等效为由定义的偏移距接收数据. 表示中心点为处的RMS速度,t0表示零偏移距时间.其中坐标 是由超道集中心的源-检波器坐标转换为偏移距-中心点坐标,定义如下:

(8a)

(8b)

(8c)

(8d)

图 3所示为超道集NMO前后波场图,该超道集源-检波器位置如图 2所示,为五维数据体.为了将该五维数据体显示为二维,源-检波器坐标排列顺序由慢至快为ys, xs, yr, xr.图 3a3b分别为原始和NMO后超道集,二者在大时间处只有微小差异,但是在浅层小时间处(小于1 s)图 3a中略有倾斜的反射信号在图 3b中被校平.

图 3 超道集数据及部分NMO (a)原始数据;(b)部分NMO后数据.

经过部分NMO校正之后的超道集已经没有了局部偏移距的差异(Δhxhy=0).因此,超道集数据可进一步简化为,对局部中心点坐标倾斜叠加公式为

(9)

公式(9)即为本文采用的无反假频操作的共偏移距倾斜叠加公式.

3 三角低通滤波算子

对三维情况,本文采取如下偏移算子积分反假频准则:

(10)

其中dx和dy分别为观测系统crossline和inline间隔(12.5 m和25 m).若不产生偏移算子假频的最大频率已知,Lumley等(1994)提出可应用三角滤波器在时域对地震数据进行实时滤波,该滤波算子Z域表达式为

(11)

其中α=(k+1)2为归一化因子,k是滤波器算子长度的一半,由公式(10)中的最大频率决定,k= ,其中

(12)

表达式g(Z)可以等效分成三个部分,

(13)

其中每个分式可视为独立部分,三部分可以对地震道按顺序依次独立进行,而1/(1-Z-1)和1/(1-Z)表示对地震道的因果和反因果积分且与滤波器长度k无关,因此可以作为预处理应用于地震数据.

三角滤波器振幅谱为(Lumley et al., 1994)

(14)

由(14)式可知,三角滤波器的陷波频率在如下位置:

(15)

Lumley等(1994)建议将fmax设置为滤波器的第一个陷波频率f1.图 4所示为k取3, 4, 5时的滤波器频谱图,由图中可见三角滤波器可有效衰减数据中的高频成分,从而达到滤波的目的.但其对有效频带(低于第一个陷波频率部分)高频数据有一定衰减,因此会对没有发生频散数据频率部分能量有所损伤.即使三角滤波器的第一个陷波频率在整个数据有效频带之外,应用该滤波器后仍然会降低数据高频能量,因此三角滤波后的数据成像结果表现为陡倾角结构分辨率有所下降.

图 4 三角滤波器频谱图,其中k分别取3, 4, 5 Fig. 4 Spectrum of the triangle filter with k for 3, 4, 5

三维反假频准则公式(10)是在满覆盖均匀采样下的准则,实际有限孔径以及非规则采样数据不同超道集中的等效采样间隔各有不同,如本次测试所用数据源-检波器crossline和inline方向均为100 m的超道集,其最少包含1道数据,最多为107道数据.即使对满覆盖均匀采样等理想情况下,由于三角滤波反假频操作对高频数据“过份”压制,准则(10)亦不可百分之百执行.鉴于该问题的重要性和复杂性,我们修改公式(10)在其中加入权重系数来对反假频的程度有所控制,式(10)修改如下:

(16)

其中wxwy是引入的两个权重,分别控制crossline和inline方向上反假频量.

反假频倾斜叠加算法步骤总结如下:

(1) 对于超道集中给定地震道进行因果积分和反因果积分,得到

(2) 对给定的(pmx, pmy),计算倾斜叠加时移量t_shift=pmxΔmx+pmyΔmy

(3) 对给定的(pmx, pmy),计算反假频滤波算子时移量,以及归一化因子α=(k+1)2,其中最大频率

(4) 根据公式(13)累加至倾斜叠加结果

在反假频倾斜叠加计算时,时移量t_shift和k有可能不是采样点的整数倍,若取最近的整数点,会引起较大误差.本文通过内插来减小地震数据的时间采样间隔,以此提高倾斜叠加和反假频精度.在内插数据的基础上执行完反假频倾斜叠加后,再将倾斜叠加后数据减采样至原始间隔.

图 5所示为图 3中超道集数据倾斜叠加结果,其中crossline和inline角度范围均为-65°~65°,角度间隔为5°.图 5a为无反假频倾斜叠加结果,为三维数据体.为了显示成二维图像,设置慢度坐标pmx为快维,pmy为慢维,可以看出明显的随机噪声分布在所有角度分量上.图 5b图 5c为反假频权重系数全部取1.0和1.5的倾斜叠加结果,应用三角反假频后随机噪声得到明显压制.图 5c相对图 5b在大角度能量上衰减更强,反假频后大角度倾斜叠加数据分辨率明显低于小角度.为了显示细节,图 5d, 5e5f分别选取图 5a, 5b5c中的部分慢度倾斜叠加结果放大显示.由于本超道集数据采样稀疏,倾斜叠加不仅有表现为随机噪声的算子假频,还有在大角度出现的数据假频(如图 5d中椭圆标识部分).反假频操作不仅压制了偏移假频,而且还能压制部分数据假频.

图 5 三维倾斜叠加结果,其中crossline和inline角度范围均为-65°~65°,角度间隔为5° (a)无反假频;(b)反假频权重系数wx=1.0,wy=1.0;(c)反假频权重系数wx=1.5,wy=1.5; (d)选取(a)图中部分慢度系数放大显示;(e)选取(b)图中部分慢度系数放大显示;(f)选取(c)图中部分慢度系数放大显示. Fig. 5 Three dimension slant stack results, the crossline and inline angles are from -65° to 65° with 5° interval (a) Without anti-aliasing; (b) With anti-aliasing weighting factor wx=1.0, wy=1.0; (c) With anti-aliasing weighting factor wx=1.5, wy=1.5; (d) Partial slowness index result demonstration from (a); (e) Partial slowness index demonstration from (b); (f) Partial slowness index demonstration from (c).
4 三维海上实际数据偏移效果试验

本测试所用工区inline和crossline方向长度分别为20.75 km和13.75 m,数据时间采样间隔4 ms,记录长度6 s,最大偏移距6 km.某一超道集源-检波器位置分布如图 2所示,该超道集inline方向偏移距为600 m,crossline偏移距0 m.超道集大小的选择一般根据模型复杂度、射线束参数设置、目标层深度以及成像精度要求和计算时间综合考虑(Sun et al., 2000).选择超道集大小的一个直接有效方法即是应用不同大小超道集划分对同一块测试数据进行多次偏移.本文选择一条inline线作为测试线,所用数据孔径为crossline方向全数据,inline范围为该测试线两侧各4 km内数据,所选择孔径范围内数据总大小为352 G.综合考虑计算能力和成像质量,本次偏移超道集大小crossline和inline均为100 m,角度范围均为-40°到40°,采样间隔10°,倾斜叠加后数据经压缩存储后为432 G,是原始数据的1.22倍.

图 6展示的是对该数据一条inline线共偏移距射线束偏移成像结果,偏移输出网格深度采样间隔为2 m,总深度6 km,crossline采样间隔12.5 m, 输出范围为13.75 km.图 6a偏移基于无反假频波束形成,图 6b图 6c基于不同权重的反假频波束形成,分别为wx=wy=1.0和wx=wy=1.5.无假频成像结果中整个剖面上随机噪声严重,尤其是浅层薄层结构被掩盖难于辨识.反假频后成像结果(图 6b, 6c)随机噪声被明显压制,浅层水平结构辨识度增高.图 7为选择该inline线中部分区域并用波形图放大显示其结构细节.在crossline 1.8 km、深度1.5 km左右的断层(图 7a, 7c, 7e)在反假频结果中可被明显辨识(图 7b).对比图 6b图 6c图 6c中浅部噪声被进一步压制,浅部薄层结构分辨率进一步增强,但一些陡倾角结构明显被削弱(图 7b, 7d, 7f).

图 6 某海上实际数据共偏移距射线束三维成像一条inline线结果 (a)无反假频波束形成;(b)反假频波束形成(反假频系数wx=1.0,wy=1.0);(c)反假频波束形成(反假频系数wx=1.5, wy=1.5). Fig. 6 One inline section from the three dimension common offset beam imaging of the marine field data (a) Beam forming without anti-aliasing; (b) Anti-aliasing beam forming (weighting factor wx=1.0, wy=1.0); (c) Anti-aliasing beam forming (weighting factor wx=1.5, wy=1.5).
图 7 海上地震数据三维成像结果(图 6)部分区域波形图放大显示 (a, b)无反假频波束形成;(c, d)反假频波束形成(反假频系数wx=1.0,wy=1.0);(e, f)反假频波束形成(反假频系数wx=1.5, wy=1.5). Fig. 7 Wiggle plot of partial section from Fig. 6 (a, b) Beam forming without anti-aliasing; (c, d) Anti-aliasing beam forming (weighting factor wx=1.0, wy=1.0); (e, f) Anti-aliasing beam forming (weighting factor wx=1.5, wy=1.5).

图 8所示为对应图 6成像结果的共成像点偏移距道集,每个成像点输出20个偏移距,最小偏移距200 m,偏移距间隔200 m.无假频道集中(图 8a)随机噪声严重且分布于所有偏移距上,尤其是浅部和深部的弱反射信号难于辨识,使AVO分析或基于偏移距道集的速度分析难于进行.经反假频后道集反射信号连续性和相干性明显增强(图 8b, 8c), 反假频后小偏移距分辨率变化较小,大偏移距能量和分辨率均有所降低.

综合考虑偏移效果和效率,我们将反假频系数为(wx=1.0,wy=1.0)应用于体偏移,该套参数压制了部分偏移噪声,而且对陡倾角成像分辨率和能量保持较好,剩余噪声可通过叠后剖面去噪方法进一步去除(王伟,2012).图 9所示为该海洋数据体偏移结果,crossline方向采样如上文所述,体偏选择了以图 6中所示inline为中心的101条线输出,在测试线两侧各50条,总长度2.5 km.所选深度分别为500 m (图 9a), 3 km (图 9b)和5 km (图 9c),可以看到浅层仍保留了一定偏移噪声,在3 km和5 km深度剖面上信噪比较高,偏移噪声得到有效压制.

图 8 共成像点偏移距道集 (a)无反假频波束形成;(b)反假频波束形成(反假频系数wx=1.0,wy=1.0);(c)反假频波束形成(反假频系数wx=1.5, wy=1.5). Fig. 8 Common imaging point offset gather (a) Beam forming without anti-aliasing; (b) Anti-aliasing beam forming (weighting factor wx=1.0, wy=1.0); (c) Anti-aliasing beam forming (weighting factor wx=1.5, wy=1.5).
图 9 三维偏移深度切片(反假频系数wx=1.0,wy=1.0) (a)深度500 m;(b)深度3 km;(c)深度5 km. Fig. 9 Depth slice of the three dimension migration (anti-aliasing weighting factor wx=1.0, wy=1.0) (a) Depth 500 m; (b) Depth 3 km; (c) Depth 5 km.
5 结论与讨论

本文提出了数据域共偏移距波束形成中的反假频策略,并将其应用于三维海上实际数据偏移成像,在叠加成像和偏移距道集上均取得了明显效果.该三角反假频滤波算子直接在时域实时完成,在对地震数据进行因果和反因果积分后,可以直接与时域倾斜叠加波束形成结合完成,避免了将数据变换至频域进行滤波操作.然而,在观测系统源和检波器采样非常稀疏,反假频操作对陡倾角结构的影响表现为分辨率降低、能量减弱.因此,本文提出在反假频公式中加入权重系数,以对反假频程度有所控制.最佳反假频权重系数可以通过选择测试线进行多次测试决定.本文反假频方法是在叠前数据域完成,根据地表速度决定的慢度分量在数据域滤波,而在偏移过程中的反假频可结合算子在实际地下传播倾角进行更精确滤波,可以减弱数据域反假频中噪声压制和分辨率降低之间的折衷效应,进一步提升反假频精度.

致谢

谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.非常感谢朱兆林在程序编写上的大力帮助,感谢陈文超关于地震信号噪声压制的讨论.感谢中海油研究总院有限责任公司提供三维实际测试数据.

References
Abma R, Sun J, Bernitsas N. 1999. Antialiasing methods in Kirchhoff migration. Geophysics, 64(6): 1783-1792. DOI:10.1190/1.1444684
Biondi B. 2001. Kirchhoff imaging beyond aliasing. Geophysics, 66(2): 654-666. DOI:10.1190/1.1444956
Casasanta L, Gray S H. 2015. Converted-wave beam migration with sparse sources or receivers. Geophysical Prospecting, 63(3): 534-551. DOI:10.1111/1365-2478.12226
Claerbout J F. 1992. Earth Soundings Analysis-Processing Versus Inversion. Boston: Blackwell Scientific Publications, Inc.
Claerbout J F. 1997. Anti aliasing. Stanford Exploration Project, Report 73. 395-415.
Claerbout J F. 2010. Basic Earth Imaging. http://sepwww.stanford.edu/sep/prof/.
Gray S H. 1992. Frequency-selective design of the Kirchhoff migration operator. Geophysical Prospecting, 40(5): 356-571.
Gray S H, Xie Y, Notfors C, et al. 2009. Taking apart beam migration. The Leading Edge, 28(9): 1098-1108. DOI:10.1190/1.3236380
Gray S H. 2013. Spatial sampling, migration aliasing, and migrated amplitudes. Geophysics, 78(3): S157-S164. DOI:10.1190/geo2012-0451.1
Liu C, Li P, Liu Y, et al. 2013. Iterative data interpolation beyond aliasing using seislet transform. Chinese Journal of Geophysics, 56(5): 1619-1627. DOI:10.6038/cjg20130519
Liu L N, Chen S M, Zhang J F. 2006. Steep dipping structure depth imaging using wave equation based migration schemes on 3-D sparse dataset. Chinese Journal of Geophysics, 49(5): 1452-1459.
Liu X W, Liu H, Liu B. 2004. A study on algorithm for reconstruction of de-alias uneven seismic data. Chinese Journal of Geophysics, 47(2): 299-305.
Lumley D E, Claerbout J E, Bevc D. 1994. Anti-aliased Kirchhoff migration. //64th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1282-1285.
Nimsaila K, Shao G F, Huang R X, et al. 2015. Premigration data anti-aliasing for reverse time migration. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 4641-4645.
Pica A. 1996. Model-based anti-aliased Kirchhoff 3D PSDM. //58th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
Robinson E A, Treitel S. 2008. Digital Imaging and Deconvolution: The ABCs of Seismic Exploration and Processing. Tulsa, Oklahoma, U. S. A: Society of Exploration Geophysicists.
Schneider W A. 1978. Integral formulation for migration in two and three dimensions. Geophysics, 43(1): 49-76. DOI:10.1190/1.1440828
Sun Y H, Qin F H, Checkles S, et al. 2000. 3-D prestack Kirchhoff beam migration for depth imaging. Geophysics, 65(5): 1592-1603. DOI:10.1190/1.1444847
Wang H Z, Cai J X, Kong X N, et al. 2010. An implementation of Kirchhoff integral prestack migration for large-scale data. Chinese Journal of Geophysics, 53(7): 1699-1709. DOI:10.3969/j.issn.0001-5733.2010.07.021
Wang H Z, Feng B, Liu S Y, et al. 2015. Characteristic wavefield decomposition, imaging and inversion with prestack seismic data. Chinese Journal of Geophysics, 58(6): 2024-2034. DOI:10.6038/cjg20150617
Wang H Z, Feng B, Wang X W, et al. 2017. The theoretical framework of characteristic wave inversion imaging. Geophysical Prospecting for Petroleum, 56(1): 38-49.
Wang W. 2012. Research on seismic data noise attenuation and geometric attributes extraction methods[Ph. D. thesis] (in Chinese). Xi'an: Xi'an Jiaotong University.
Wu B Y, Zhu Z L, Yang H, et al. 2014. High resolution beam forming for 3D common offset Kirchhoff beam migration. //84th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3837-3841.
Wu B Y, Yang H, Zhu Z L, et al. 2015. Local slant stack Kirchhoff common offset beam migration: Application to the SEAM model. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 4283-4287.
Yilmaz Ö. 2001. Seismic data analysis: processing, inversion, and interpretation of seismic data. SEG, Volume (Ⅰ).
Zhang Y, Gray S, Sun J, et al. 2001. Theory of migration anti-aliasing. //71th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 997-1000.
Zhou Y, Etgen J, Whitcombe D, et al. 2013. Dip-adaptive operator anti-aliasing for Kirchhoff migration. //83th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3629-3696.
刘财, 李鹏, 刘洋, 等. 2013. 基于seislet变换的反假频迭代数据插值方法. 地球物理学报, 65(5): 1619-1627. DOI:10.6038/cjg20130519
刘礼农, 陈树民, 张剑锋. 2006. 稀疏采样下陡角度构造的波动方程深度偏移成像. 地球物理学报, 2006: 1452-1459. DOI:10.3321/j.issn:0001-5733.2006.05.024
刘喜武, 刘洪, 刘彬. 2004. 反假频非均匀地震数据重建方法研究. 地球物理学报, 47(2): 299-305.
王华忠, 蔡杰雄, 孔祥宁, 等. 2010. 适于大规模数据的三维Kirchhoff积分法体偏移实现方案. 地球物理学报, 53(7): 1699-1709. DOI:10.3969/j.issn.0001-5733.2010.07.021
王华忠, 冯波, 刘少勇, 等. 2015. 叠前地震数据特征波场分解、偏移成像与层析反演. 地球物理学报, 58(6): 2024-2034. DOI:10.6038/cjg20150617
王华忠, 冯波, 王雄文, 等. 2017. 特征波反演成像理论框架. 石油物探, 56(1): 38-49.
王伟. 2012. 地震资料噪声衰减与几何属性分析方法研究[博士论文]. 西安: 西安交通大学.