2. 中国石油大学 (华东) 地球物理系, 青岛 266580
2. Department of Geophysics, China University of Petroleum (East China), Qingdao 266580, China
地震偏移成像技术是地震勘探数据处理的三大基本技术之一,从古典的射线类偏移到现阶段的波动方程偏移,主要围绕着提高成像精度和减少计算成本两个方面发展.虽然基于波动方程数值解的分步傅里叶法、傅里叶有限差分法、广义屏法以及逆时偏移法等具有较高的成像精度,但较低的计算效率和较大的内存开销使它们在现今的计算机水平下还难以满足实际的批量生产需求.Kirchhoff偏移和高斯束偏移等射线类方法以其高效、稳健、灵活、强适应性的独特优势,在实际生产中仍发挥着主要作用 (Gray和方伍宝, 2002;徐升和Lambaré,2006;孙建国,2012;李振春,2014).
Schneider (1978)首先对Kirchhoff偏移理论进行了系统地阐述,之后Carter和Frazer (1984)、Gray (1986)、Gray和May (1994)、Takahashi (1995)、Chang (1998)、Kim和Cook (1999)、Sethian和Popovici (1999)等学者对其在非均匀介质中的应用和走时表的求取、存储以及调用提出了一些建设性的意见,使其在实际生产中得到了迅速发展.由于采用运动学射线追踪计算单值走时表,Kirchhoff偏移存在射线覆盖不到的阴影区、振幅奇异化的焦散区以及不能处理多波至等一系列问题.为克服Kirchhoff方法的缺陷,Hill (1990, 2001)、Nowack等 (2003)和Gray (2005)(Gray and Bleistein, 2009) 等提出了高斯束偏移,该方法不仅保留了射线类偏移的高效、灵活性以及陡倾反射界面准确成像的优势,而且还具有与波动方程类偏移可比拟的成像精度.
复杂构造探区深度域速度型的建立一直勘探地球物理界研究的热点和难点,可行的方法之一是交互式的速度建模,即对不同的速度模型进行快速深度域偏移,使它们在短时间内得到检验和标定,并根据成像效果优化原始的速度模型.为提高偏移速度,国内外许多学者在成像技术上做了大量有益尝试,并取得一定成绩.Wang和Pann (1996)将匹配追踪应用于Kirchhoff偏移中,提出了一种输入型的快速成像方法,使计算速率得到了一定提高.Li等 (1998)进一步优化了匹配追踪算法,提出了双参数控制的快速Kirchhoff深度偏移方法.Hua和McMechan (2001, 2003) 根据射线的出射方向,只在特定的角度进行射线追踪和成像,使传统的Kirchhoff深度偏移速度提高了近十倍.Gao等 (2006, 2007) 采用先拾取同相轴再偏移的策略,提出了一种快速射线束偏移方法,也使偏移速率提高了近一个数量级.Liu等 (2013)根据tau-p域地震数据的稀疏表征方式,提出了特征波分解和成像方法,也使偏移速率得到一定提高.
对上述理论方法进行研究总结后,本文将匹配追踪算法和tau-p域同相轴拾取技术应用于常规高斯束偏移中,提取了一种快速深度域成像方法.该方法首先利用匹配追踪计算出束中心地震道同相轴的位置 (波峰或波谷),然后根据tau-p域地震记录的振幅特征进一步确定出这些同相轴的优势出射方向,并再次利用匹配追踪算法将这些同相轴的优势方向分量提取出来,最后采用高斯束方法将这些提取的同相轴偏移至成像域.相对现有的快速偏移方法,本文方法应用了两次匹配追踪分解,使同相轴的拾取更加准确可靠.此外,该方法在偏移时只对同相轴的优势方向进行波场延拓和成像,减少了成像时的倾角扫描范围,提高了计算效率.同时,这种选择性的成像方式也可以压制相关噪声,消除偏移假象,提高主要构造的成像质量.模型算例和实际资料的处理结果验证了本文方法的有效性和一定适应性.
1 方法原理 1.1 匹配追踪稀疏分解匹配追踪是将输入信号在超完备原子库中进行自适应稀疏分解.超完备原子库可由所选择的母子波经扩张、调制和平移生成,在地震资料处理中,常用的母子波有Morlet子波和Ricker子波等.由于Ricker子波波形简单,延迟时间短,收敛较快,并且与地震信号具有相似的结构特征 (韩海英等,2014),在此采用Ricker子波作为匹配追踪稀疏分解的母子波.设D为超完备子波库,定义D={gi:i=1, 2, 3…},其中gi为由主频fi和中心时间τi确定的Ricker型时频原子 (如图 1所示),具体表达式为
(1) |
设s(t) 为待分解的地震信号,经N次迭代分解后得到:
(2) |
其中gi为第i次迭代得到的最佳匹配子波,ai为对应的振幅,RN[s(t)]为N次分解后的残差,当N足够大时,RN[s(t)]可看成有效信号中的干扰噪声.式 (2) 表明,地震信号可以近似看成由fi、τi和ai双参数控制Ricker子波的线性组合,通过在完备原子库中扫描,能够确定出最佳匹配子波的主频、延迟时间和振幅参量.匹配追踪的实质就是这三个控制参数的匹配寻优过程.
1.2 tau-p域同相轴拾取匹配追踪仅仅将地震数据在时间维上进行了稀疏分解,为了进一步在空间上压缩地震数据,提高偏移速率,可根据tau-p域地震记录的振幅特征对同相轴的优势方向分量的进行动态拾取.
在同相轴拾取时,首先应对束内地震道集 (见图 2a) 进行局部tau-p变换,然后在tau-p域内 (见图 2b) 求得每一个同相轴振幅极值所对应的射线参数P (见图 2b的红色线条),最后根据得到的P值序列再次利用匹配追踪算法将不同时刻的时频原子从tau-p域记录上提取出来.经过在时间维和空间维上的稀疏分解,束内地震道集就由原来成千上万个采样点压缩成了数十个或数百个 (取决于波场复杂性) 时频原子 (见图 2c和2d).
图 3为采用本文方法和传统方法计算同相轴时间延迟的对比结果,可以看出,传统方法仅将地震道按一定的时间间隔划分成不同的时窗,不能保证在每一个时窗内提取的都是同相轴的波峰或波谷 (见图 3b红色椭圆),并且对于简单的地震记录还存在分解的零道 (见图 3b红色箭头)(Hua et al., 2001;Mo et al., 2009),而本文采用了自适应的匹配追算法,可以准确地计算出同相轴的波峰或波谷 (见图 3a红色椭圆),具有较高的提取精度.
对地震道集分别进行tau-p域射线参数的计算和同相轴的拾取,单炮地震记录就可稀疏表征为
(3) |
式中x S、x R和xL分别为炮点、检波点和束中心的坐标,gi(pS, pR, t) 为第i个时频原子,ai为其振幅,pS和pL分别炮点和束中心的射线参数 (水平慢度).
根据Červeny等 (1982, 1984) 的相关文献可知,由高斯束表征的单方向频率域波场为
(4) |
其中i为虚数单位,Δp为射线参数间隔,pz为垂直慢度,AX′和TX′分别为高斯束的复值振幅和复值走时.
根据Gray和Bleistein (2009)提出的真振幅高斯束偏移公式,稀疏分解后地震记录的反向延拓波场和震源激发的正向波场可表示为
(5) |
其中x为成像点坐标,θS和θL分别为炮点和束中心的射线出射角 (见图 4),vS和vL分别为炮点和束中心处的速度,“*”代表复值共轭,gi(pS, pR, w) 为gi(pS, pR, t) 傅里叶变换.
将 (4) 带入 (5) 式,并应用反褶积成像条件:
(6) |
可得到基于匹配追踪分解和tau-p域同相轴拾取的快速高斯束偏移公式为
(7) |
该式表明,本文方法不再采用传统的以采样点为单位的数据映射方式,而是以分解的时频原子为单位进行偏移处理,因此可以进一步提高计算效率.
1.4 计算流程快速高斯束偏移的基本计算流程可分为如下几个步骤:
(1) 地震记录的稀疏表征.首先对输入的叠前共炮点地震记录按经典的5道或7道进行小道集的划分 (每一个小道集可看成一条地震波束),然后对共炮点道集的束内地震道采用tau-p变换和两次匹配追踪算法提取出每一个同相轴的优势方向分量,并存储它们的检波点射线参数,最后在共检波点道集内寻找出对应的炮点-束中心道位置,并根据其tau-p域的振幅特征计算每一个同相轴的炮点射线参数.
(2) 射线追踪.利用运动学射线追踪构建包含空间坐标和走时信息的中心射线场,利用动力学射线追踪沿每条中心射线计算高斯束的动力射线参数,并进一步根据高斯束表达式 (Červeny等,1982) 构建包含复值振幅和复值走时的旁轴射线场.
(3) 正反向波场延拓.对拾取的每一个同相轴,根据其炮点和检波点射线参数从 (2) 获得的旁轴射线场中筛选出合适的角度分量,并利用式 (5) 计算高斯束表征的正反向波场.
(4) 扫描成像.对深度域成像网格点进行扫描,当正向和反向波场均有能量贡献时,采用反褶积成像条件计算此网格点的成像值.
将拾取的所有同相轴均按照 (2)~(4) 步骤进行成像并叠加,就可获得最终的成像结果.
2 模型算例在实现基于匹配追踪和tau-p域同相轴拾取的快速高斯束偏移的基础上,本文首先通过如图 5a所示的洼陷模型对方法的正确性和有效性进行了测试.该模型网格大小1201×301,纵横向间隔分别为10 m和5 m,声波正演记录共201炮,每炮301道,道间距为10 m,采样间隔为1 ms,记录长度2.5 s.
图 5b为使用SU开源软件包中sukdmig2d程序得到的Kirchhoff偏移结果.可以看出,Kirchhoff偏移在对反射界面成像的同时,也引入了大量偏移噪声,主要包括:(1) 由于广角反射与临界角内反射波相位不匹配 (图 6a红色箭头为广角反射波,图 6b蓝色椭圆为广角反射偏移结果) 造成的浅层叠加噪声 (见图 5b蓝色矩形框);(2) 由于界面间多次反射造成的多次波 (见图 5b红色箭头);(3) 由于全方向划弧造成的相关噪声 (见图 5b蓝色椭圆).图 5c和5d为分别采用常规高斯束方法和本文方法得到的偏移结果.可以看出,相对Kirchhoff方法,高斯束偏移能够消除大部分由于全方向划弧造成的相干噪声,但由于理论缺陷仍不能去除浅层叠加噪声 (见图 5c蓝色矩形框) 和多次反射波 (见图 5c红色箭头),而本文方法在对地震数据进行稀疏分解后只在优势方向进行高斯束成像,从根本上避免了三种偏移噪声的出现,大大改善了成像结果的信噪比 (见图 5d).但是,如图 5d中红色椭圆所示,本文方法使信噪比提高的代价是牺牲了一定的成像精度.
为了进一步验证本文方法对复杂地质构造的适应性,在此对Marmous数据集进行测试.该数据集共240炮,每炮96道接收,右侧激发,最小炮检距200 m,最大炮检距2575 m,道间距25 m,采样间隔4 ms,记录长度3 s.图 7a为其速度模型,网格大小为737×750,纵横向间隔分别为4 m和12.5 m,该模型不但具有尖灭、背斜和侵入体等典型的地质构造,还包含复杂的断裂系统,是测试深度域成像方法的经典模型.图 7b-d分别为采用不同偏移方法得到的成像结果,通过仔细观察对比可得到以下认识:(1) Kirchhoff偏移虽然可以对主要构造清晰成像,但是其偏移结果中含有大量相关噪声 (见图 7b红色椭圆),并且由于其不能对多波至成像,导致复杂断裂系统中的陡倾反射界面模糊不清 (见图 7b蓝色箭头);(2) 由于高斯束偏移在具体实现时采用定向射线追踪和逐角度扫描成像技术,因此可以对陡倾断面准确成像 (见图 7c蓝色箭头),并且在一定程度上也减弱了成像结果中的相关噪声 (见图 7c红色椭圆);(3) 基于匹配追踪和tau-p域同相轴拾取技术,本文方法首先对地震数据进行时间和空间上的稀疏分解,然后将主要反射能量用于高斯束成像,在保证对地下主要构造和陡倾断面清晰成像的同时 (见图 7d蓝色箭头),大大改善了成像结果的信噪比 (见图 7d).
图 8为使用不同方法对某探区实际资料进行偏移处理得到的成像结果,可以看出,本文方法不仅可以对地下主要反射界面清晰成像 (见图 8d蓝色箭头),也能够消除了传统成像方法中大量的数值噪声和偏移假象 (见图 8红色椭圆).
表 1为计算时间和压缩前后数据量的对比结果,可以看出,地震数据的压缩比例主要取决于地震波场的复杂性,当波场较为简单时数据压缩比可达到上百倍,比如洼陷模型的压缩比为115.6倍,而当波场较为复杂时,地震数据的压缩比也可控制在10倍左右,比如Marmous数据集的压缩比为13.4倍,实际资料的压缩比为14.1倍.表中不同成像方法的计算时间表明,相对Kirchhoff和常规高斯束偏移,本文方法的计算效率提高了一个数量级以上,对某些数据集可达到数十倍,这主要是因为本文方法采用了匹配追踪、tau-p域拾取和以子波为单位的快速偏移等一系列加速策略.
基于匹配追踪分解和tau-p域拾取算法,本文发展了一种快速叠前深度域高斯束成像方法,并通过对典型模型的测试和对实际资料的处理验证了该方法的有效性和一定适应性.模型测试和实际资料对比分析结果表明:(1) 根据tau-p域地震记录的振幅特征,本文采用两次匹配追踪算法对叠前地震数据进行时间和空间上的稀疏分解,使叠前数据压缩了数十倍,为后续的快速偏移提供了数据基础;(2) 由于本文方法只对同相轴的优势角度分量进行波场的延拓和成像,并采用了以子波为单位的快速偏移策略,可使深度域偏移的计算效率提高一个数量级以上,为今后交互式的偏移速度分析和建模奠定了基础;(3) 由于本文方法只对主要的反射能量进行成像,可以在保证对地下主要构造和陡倾反射面清晰成像的同时,从根本上消除Kirchhoff和传统高斯束偏移中由于广角反射相位不匹配造成的浅层叠加噪声、多次波成像噪声和全方向划弧引起的相干噪声,大大改善了成像结果的信噪比.
3.2虽然本文实现了基于匹配追踪和tau-p域同相轴拾取的快速高斯束偏移,但还需要在如下几个方面做进一步的深入研究:(1) 虽然本文方法可以对地下主要的构造清晰成像,但成像的精度和分辨率稍低,因此需要进一步从稀疏分解和成像条件两个方面探索更佳的偏移策略;(2) 在快速高斯束偏移的基础上,进一步研究相应的偏移速度分析方法和交互的速度建模方法.
致谢 感谢油气重大专项 (2016ZX05002-005-07HZ,2016ZX05014-001-008HZ)、国家重点基础研究发展计划 (973) 子课题 (2014CB239006) 和国家自然基金 (41274124) 联合资助,感谢科罗拉多矿院CWP课题组提供的SU开源软件包.[] | Carter J A, Frazer L N. 1984. Accommodating lateral velocity changes in Kirchhoff migration by means of Fermat's principle[J]. Geophysics, 49(1): 46–53. DOI:10.1190/1.1441560 |
[] | ČervenýV, PopovM M, PšenčíkI. 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach[J]. Geophysical Journal International, 70(1): 109–128. |
[] | Gao F, Zhang P, Wang B, et al. 2007. Interactive seismic imaging by fast beam migration[C].//69th EAGE Conference & Exhibition. |
[] | Gao F C, Zhang P, Wang B, et al. 2006. Fast Beam migration-a step toward interactive imaging[C].//2006 SEG Annual Meeting, Society of Exploration Geophysicists. Expanded Abstracts, 2470-2474. |
[] | Gray S H. 1986. Efficient traveltime calculations for Kirchhoff migration[J]. Geophysics, 51(8): 1685–1688. DOI:10.1190/1.1442217 |
[] | Gray S H. 2005. Gaussian beam migration of common-shot records[J]. Geophysics, 70(4): S71–S77. DOI:10.1190/1.1988186 |
[] | Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration[J]. Geophysics, 74(2): S11–S23. DOI:10.1190/1.3052116 |
[] | Gray S H, FANG Wu-Bao. 2002. Seismic migration problems and solutions[J]. Progress in Exploration Geophysics, 25(2): 44–60. |
[] | Gray S H, May W P. 1994. Kirchhoff migration using eikonal equation traveltimes[J]. Geophysics, 59(5): 810–817. DOI:10.1190/1.1443639 |
[] | HAN Hai-Ying, WANG Zhi-Zhang, WANG Zong-Jun, et al. 2014. Matching pursuit based on Ricker-like seismic wavelet[J]. Petroleum Exploration, 53(1): 17–25. |
[] | Hill N R. 1990. Gaussian beam migration[J]. Geophysics, 55(11): 1416–1428. DOI:10.1190/1.1442788 |
[] | Hill N R. 2001. Prestack Gaussian-beam depth migration[J]. Geophysics, 66(4): 1240–1250. DOI:10.1190/1.1487071 |
[] | Hua B L, McMechan G A. 2001. Parsimonious 2-D poststack Kirchhoff depth migration[J]. Geophysics, 66(5): 1497–1503. DOI:10.1190/1.1487095 |
[] | Hua B L, McMechan G A. 2003. Parsimonious 2D prestack Kirchhoff depth migration[J]. Geophysics, 68(3): 1043–1051. DOI:10.1190/1.1581075 |
[] | Kim S, Cook R. 1999. 3-D traveltime computation using second-order ENO scheme[J]. Geophysics, 64(6): 1867–1876. DOI:10.1190/1.1444693 |
[] | Li X G, Wang B, Pann K, et al. 1998. Fast migration using a matching pursuit algorithm[C].//1998 SEG Annual Meeting, Society of Exploration Geophysicists. Expanded Abstracts, 1732-1735. |
[] | LI Zhen-Chun. 2014. Research status and development trends for seismic migration technology[J]. Oil Geophysical Prospecting, 49(1): 1–21. |
[] | Liu S Y, Feng B, Wang H Z, et al. 2013. Characteristic wave decomposition and imaging[C].//75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013. |
[] | Madariaga R. 1984. Gaussian beam synthetic seismograms in a vertically varying medium[J]. Geophysical Journal International, 79(2): 589–612. DOI:10.1111/j.1365-246X.1984.tb02243.x |
[] | Mo L W, Shurtleff R, Hays D. 2009. 3-D shot-profile Kirchhoff fast beam depth migration[C].//CPS/SEG Beijing 2009 International Geophysical Conference and Exposition. |
[] | Nowack R L, Sen M K, Stoffa P L. 2003. Gaussian beam migration for sparse common-shot and common-receiver data[C].//SEG Int. Exposition and 73rd Annual Meeting. Expanded Abstracts, 1114-1117. |
[] | Schneider W A. 1978. Integral formulation for migration in two and three dimensions[J]. Geophysics, 43(1): 49–76. DOI:10.1190/1.1440828 |
[] | Sethian J A, Popovici A M. 1999. 3-D traveltime computation using the fast marching method[J]. Geophysics, 64(2): 516–523. DOI:10.1190/1.1444558 |
[] | SUN Jian-Guo. 2012. The history, the state of the art and the future trend of the research of Kirchhoff-type migration theory:A comparison with optical diffraction theory, some new results and new understanding, and some problems to be solved[J]. Journal of Jilin University (Earth Science Edition), 42(5): 1521–1552. |
[] | Sun Y H, Qin F H, Checkles S, et al. 2000. 3-D prestack Kirchhoff beam migration for depth imaging[J]. Geophysics, 65(5): 1592–1603. DOI:10.1190/1.1444847 |
[] | Takahashi T. 1995. Prestack migration using arrival angle information[J]. Geophysics, 60(1): 154–163. DOI:10.1190/1.1443742 |
[] | Wang B, Pann K. 1996. Kirchhoff migration of seismic data compressed by matching pursuit decomposition[C].//Expanded Abstracts of The Technical Program, SEG 66th Annual International Meeting. Expanded Abstracts, 1642-1645. |
[] | XU Sheng, Lambaré G. 2006. True amplitude Kirchhoff prestack depth migration in complex media[J]. Chinese Journal of Geophysics, 49(5): 1431–1444. DOI:10.3321/j.issn:0001-5733.2006.05.022 |
[] | GrayS H, 方伍宝. 2002. 地震偏移问题及其解决方案[J]. 勘探地球物理进展, 25(2): 44–60. |
[] | 韩海英, 王志章, 王宗俊, 等. 2014. 基于Ricker类地震子波的匹配追踪[J]. 石油物探, 53(1): 17–25. |
[] | 李振春. 2014. 地震偏移成像技术研究现状与发展趋势[J]. 石油地球物理勘探, 49(1): 1–21. |
[] | 孙建国. 2012. Kirchhoff型偏移理论的研究历史、研究现状与发展趋势展望——与光学绕射理论的类比、若干新结果、新认识以及若干有待于解决的问题[J]. 吉林大学学报 (地球科学版), 42(5): 1521–1552. |
[] | 徐升, LambaréG. 2006. 复杂介质下保真振幅Kirchhoff深度偏移[J]. 地球物理学报, 49(5): 1431–1444. DOI:10.3321/j.issn:0001-5733.2006.05.022 |