2. 南方科技大学地球与空间科学系, 深圳 518055;
3. 中海油研究总院有限责任公司, 北京 100028;
4. 西安交通大学电子与信息工程学院波动与信息研究所, 西安 710049;
5. 海洋石油勘探国家工程实验室, 西安 710049
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
离散数字信号处理是现代地震数据处理的基础.而采样定理,又称香农采样定律或奈奎斯特采样定律又是连续信号离散化的基本依据,是实际连续信号和可供计算机处理的离散数字信号之间的桥梁.采样定理指出,对于一维时间信号,如果带宽小于采样频率的一半(即奎斯特频率),则采样信号可通过重建公式完全表示原信号(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) |
该表达式中
算子反假频在Kirchhoff偏移中得到了较为充分的发展.最简单高效的假频压制方法即是限制偏移孔径(王华忠等,2010)或者直接限制偏移算子的大角度分量,该策略对小角度反射层结构成像改善效果明显,但并不能完全消除偏移假频,而且对成像陡倾角结构损伤巨大(Abma et al., 1999).原则上,利用公式(1)反假频需要计算走时对空间的导数
射线束成像在计算效率、陡倾角成像、信噪比、算法稳定性以及灵活性等方面均有一定优势.该方法将邻近道数据同时向下偏移,对宽方位、宽频带和高密度地震数据采集偏移成像优势明显.尤其是共偏移距射线束偏移成像,地震信号线性特征强,利于射线束形成以及输出偏移距道集.波束形成是射线束偏移中关键步骤(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) |
后文为表示简单,在不引起混淆的情况下,将地震数据表示为
(3) |
其中(psx, psy)和(prx, pry)分别表示源和检波器射线参数,τ表示倾斜叠加在时间轴上的截距.射线束方法对每个超道集偏移时,可计算以超道集中心位置
积分类偏移方法相较于波动方程类偏移优势之一即是可以在更为灵活的地震数据集进行偏移,如基于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)处t(Δhx, Δhy)时刻波场平移至
(8a) |
(8b) |
(8c) |
(8d) |
图 3所示为超道集NMO前后波场图,该超道集源-检波器位置如图 2所示,为五维数据体.为了将该五维数据体显示为二维,源-检波器坐标排列顺序由慢至快为ys, xs, yr, xr.图 3a和3b分别为原始和NMO后超道集,二者在大时间处只有微小差异,但是在浅层小时间处(小于1 s)图 3a中略有倾斜的反射信号在图 3b中被校平.
经过部分NMO校正之后的超道集已经没有了局部偏移距的差异(Δhx=Δhy=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时的滤波器频谱图,由图中可见三角滤波器可有效衰减数据中的高频成分,从而达到滤波的目的.但其对有效频带(低于第一个陷波频率部分)高频数据有一定衰减,因此会对没有发生频散数据频率部分能量有所损伤.即使三角滤波器的第一个陷波频率在整个数据有效频带之外,应用该滤波器后仍然会降低数据高频能量,因此三角滤波后的数据成像结果表现为陡倾角结构分辨率有所下降.
三维反假频准则公式(10)是在满覆盖均匀采样下的准则,实际有限孔径以及非规则采样数据不同超道集中的等效采样间隔各有不同,如本次测试所用数据源-检波器crossline和inline方向均为100 m的超道集,其最少包含1道数据,最多为107道数据.即使对满覆盖均匀采样等理想情况下,由于三角滤波反假频操作对高频数据“过份”压制,准则(10)亦不可百分之百执行.鉴于该问题的重要性和复杂性,我们修改公式(10)在其中加入权重系数来对反假频的程度有所控制,式(10)修改如下:
(16) |
其中wx和wy是引入的两个权重,分别控制crossline和inline方向上反假频量.
反假频倾斜叠加算法步骤总结如下:
(1) 对于超道集中给定地震道
(2) 对给定的(pmx, pmy),计算倾斜叠加时移量t_shift=pmxΔmx+pmyΔmy;
(3) 对给定的(pmx, pmy),计算反假频滤波算子时移量
(4) 根据公式(13)累加至倾斜叠加结果
在反假频倾斜叠加计算时,时移量t_shift和k有可能不是采样点的整数倍,若取最近的整数点,会引起较大误差.本文通过内插来减小地震数据的时间采样间隔,以此提高倾斜叠加和反假频精度.在内插数据的基础上执行完反假频倾斜叠加后,再将倾斜叠加后数据减采样至原始间隔.
图 5所示为图 3中超道集数据倾斜叠加结果,其中crossline和inline角度范围均为-65°~65°,角度间隔为5°.图 5a为无反假频倾斜叠加结果,为三维数据体.为了显示成二维图像,设置慢度坐标pmx为快维,pmy为慢维,可以看出明显的随机噪声分布在所有角度分量上.图 5b和图 5c为反假频权重系数全部取1.0和1.5的倾斜叠加结果,应用三角反假频后随机噪声得到明显压制.图 5c相对图 5b在大角度能量上衰减更强,反假频后大角度倾斜叠加数据分辨率明显低于小角度.为了显示细节,图 5d, 5e和5f分别选取图 5a, 5b和5c中的部分慢度倾斜叠加结果放大显示.由于本超道集数据采样稀疏,倾斜叠加不仅有表现为随机噪声的算子假频,还有在大角度出现的数据假频(如图 5d中椭圆标识部分).反假频操作不仅压制了偏移假频,而且还能压制部分数据假频.
本测试所用工区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).
图 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深度剖面上信噪比较高,偏移噪声得到有效压制.
本文提出了数据域共偏移距波束形成中的反假频策略,并将其应用于三维海上实际数据偏移成像,在叠加成像和偏移距道集上均取得了明显效果.该三角反假频滤波算子直接在时域实时完成,在对地震数据进行因果和反因果积分后,可以直接与时域倾斜叠加波束形成结合完成,避免了将数据变换至频域进行滤波操作.然而,在观测系统源和检波器采样非常稀疏,反假频操作对陡倾角结构的影响表现为分辨率降低、能量减弱.因此,本文提出在反假频公式中加入权重系数,以对反假频程度有所控制.最佳反假频权重系数可以通过选择测试线进行多次测试决定.本文反假频方法是在叠前数据域完成,根据地表速度决定的慢度分量在数据域滤波,而在偏移过程中的反假频可结合算子在实际地下传播倾角进行更精确滤波,可以减弱数据域反假频中噪声压制和分辨率降低之间的折衷效应,进一步提升反假频精度.
致谢谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.非常感谢朱兆林在程序编写上的大力帮助,感谢陈文超关于地震信号噪声压制的讨论.感谢中海油研究总院有限责任公司提供三维实际测试数据.
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. 地震资料噪声衰减与几何属性分析方法研究[博士论文]. 西安: 西安交通大学.
|