地震勘探方法主要通过地表获取的地震波反射数据来反映地下介质的结构和物性,因此获取叠前地震数据的数量和质量成为能否获取高质量地下信息的关键环节.随着勘探程度的逐渐深入,三维宽方位角地震勘探方法的普遍应用,获取完整地震数据体的任务变得更加困难.由于实际环境和经济因素的限制,观测系统很难记录完整 (尤其近炮检距位置) 的地震波场信息.而近炮检距数据的完整性通常是基于波动理论处理方法应用的前提条件,例如自由表面多次波的预测 (SRME)(van Dedem and Verschuur, 2005),波动方程偏移.近年来,国内外学者就缺失地震数据的重建问题做出了大量的研究并取得了一定的进展.F-X域地震道插值方法 (Spitz,1991),基于Fourier变换 (高建军等,2009;Xu et al., 2005) 以及Radon变换 (张红梅和刘洪,2008) 的数据插值方法具有很好的计算特性,但随着人们对地震波场处理精度的提高,这些方法无法很好地解决复杂地震道插值的问题.分形插值数据重建方法 (李信富和李小凡,2008) 能够很好地恢复地震数据的局部信息,但是有时存在计算效率的问题.基于非一致离散傅里叶变换 (刘喜武等,2004;孟小红等,2008) 以及Radon变换 (Kabir and Verschuur, 1995;王维红等,2007) 的数据插值方法具有很好的计算特性,但仍然不能很好地解决复杂地震波场的数据缺失问题.基于干涉测量的地震数据插值方法 (Wang and Schuster, 2007; Wang et al., 2010) 和构建准一次波的插值方法 (Shan and Guitton, 2004;Curry and Shan, 2010) 以及直接利用含有多次波的地震数据构建拟地震数据的插值方法 (郭书娟等,2013),主要是利用地震数据的动力学可预测特性,将多次波中蕴含的有效信息提取出来,但是对非平稳数据的局部变化处理能力往往较弱.
地震数据的非平稳性已经引起了诸多学者的关注.预测误差滤波器 (PEF) 在数学上是地震有效信息谱的逆,可以有效地表征地震波能量,标准的预测误差滤波器是以地震数据平稳为前提 (Claerbout,1992),利用重叠时窗的方法 (Claerbout,2010) 进行数据分析,地震数据满足局部平稳假设,但是固定时窗的大小,往往限制非平稳信息的表征.非平稳预测误差滤波 (Crawley et al., 1999;Crawley,2000;Guitton et al., 2001) 已经得到了开发和应用,但最终的效果取决于大量人为参数的设定.Fomel (2002)提出利用非平稳平面波分解滤波器的方法对缺失数据进行重建,Liu和Fomel (2011)提出自适应预测误差滤波器 (APEF),这两种方法可以有效地解决带有空间假频的数据重建问题,但仍然存在不准确动力学波场重建的问题,尤其是近炮检距数据插值精度不高.
本文结合地震数据的动力学可预测特性,应用自适应预测误差滤波器从预测虚拟一次波中估算能量表征系数,利用最小二乘反演方法实现地震近炮检距缺失数据的重建.该方法可以有效地恢复缺失的近炮检距地震波场,保证恢复的数据具有准确的动力学特征.
2 理论基础 2.1 虚拟一次波构建在反射地震数据中,多次波通常被视为不需要的噪声,随着地震数据处理技术的完善,研究人员发现多次波也是一种可利用的波场,由于其多次反射的特征,能够反映一次波无法到达的地下信息,逐渐已经成为一次波勘探方法的有效补充 (Curry and Shan, 2010).利用多次波和有效波进行互相关运算,可以构建虚拟一次波.图 1给出了构建虚拟一次波的示意图.
在二维观测系统中 (图 1),横坐标表示检波点,纵坐标表示震源点,其中“
(1) |
其中,p(r1, r2, ω) 表示构建的虚拟一次波频率域数据,ω代表频率域,d(s, r2, ω) 表示s为震源,在r2接收的频率域多次波数据,
在实际地震勘探中,尽管s与r1,s与r2由于观测系统的限制不能过小,但r1与r2之间的距离可以足够小,使得虚拟一次波具有完整的近炮检距记录信息,这部分数据可以成为一次波的有效补充.
构建的虚拟一次波与理想记录在动力学的特征上近似,但并不能完全与原始数据匹配.虚拟一次波和原始数据的振幅会在整体上出现不同,由于互相关的计算,使得子波出现了平方的形式,改变了地震子波的形态.为了解决以上问题,本文利用基于整形正则化的非平稳自回归方法,设计自适应预测误差滤波器,进而表征非平稳虚拟一次波波场的谱特征,由于虚拟一次波与真实一次波的谱具有较好的相似性,因此,虚拟一次波的自适应预测误差滤波器系数可以准确地恢复缺失的地震数据信息,同时可以避免出现额外的噪声.
2.2 自适应预测误差滤波器 (APEF) 估计地震数据的重建可以表述为求解反问题的过程.利用自适应预测误差滤波 (APEF) 解决插值的方法可分为两个步骤,首先,将虚拟一次波作为初始数据,使得虚拟一次波与未知滤波系数的滤波器之间滤波结果最小,估计出滤波系数;然后利用原始数据中的已知数据作为限定条件,使得估计出的滤波器与未知数据 (原始数据中的缺失部分) 的滤波结果最小,完成对缺失数据的重建.
对于非平稳的地震数据,可以假定数据的能量谱是局部平稳的,使得自适应预测误差滤波器在局部具有尺度不变性.估计的非平稳自回归系数Bn可以体现出时间和空间 (t, x) 方向上的变化.例如:带有7个滤波系数的APEF表达形式为:
(2) |
其中,“·”代表零值.为了估计滤波系数Bn(t, x),可以通过带有正则化条件的最小二乘反演方法求取:
(3) |
其中,Sn(t, x)=S(t-i, x-j),表示非平稳的虚拟一次波数据S(t, x) 在时间和空间方向上的移动拷贝,i表示时间变化,j表示空间变化.下角标n表示时间和空间的变化,N表示i和j的总数.D表示正则化算子,ε表示标量正则化参数.所估计出的系数Bn同时反映了时间和空间上的变化.本文应用了Fomel (2009)提出的非平稳自回归,其中的整形正则化 (Fomel,2007) 算子具有较Tikhonov (1963)正则化条件更好的数学特征 (附录A).
构建的虚拟一次波相对于缺失的数据有着近似的动力学能量特征,将虚拟一次波作为输入数据,估算滤波系数,由于自适应预测误差滤波器的非平稳特性,可以很好地表征缺失一次波的非平稳能量谱.
2.3 近炮检距缺失数据重建
利用之前估计的滤波系数
(4) |
(5) |
其中,
本文利用共轭梯度法 (Hestenes and Stiefel, 1952) 计算公式 (3) 和公式 (4) 的最小二乘问题.公式 (5) 作为初始模型 (限定条件),在共轭梯度框架下,利用已知的中、远炮检距地震数据对插值结果进行约束,完成对缺失的近炮检距数据进行重建.
3 理论模型分析为了测试本文方法的有效性,这里采用Sigsbee2B模型进行测试,图 2为Sigsbee2B的速度模型,合成数据共计496炮,每炮348个接收点,时间采样间隔为8 ms,时间长度为8 s,模型中包含一个起伏海底界面,用于产生多次波.图 3为零炮检距数据 (选取空间方向8~24 km,时间方向0~8 s的记录).图 4a为在15 km位置的单炮记录,图 4b为近炮检距缺失15道的单炮记录.从图 4中可以看到,该理论数据表现出复杂的波场特征,存在具有非平稳振幅的一次波、多次波以及散射波的多种波场,常规的插值方法很难同时兼顾波场的动力学特征和非平稳数据特征.
通过公式 (1),可以利用缺失数据的多次波和一次波互相关构建虚拟一次波.为了对比,图 5a为完整数据采样下构建虚拟一次波的零炮检距数据,图 5b为近炮检距缺失条件下构建的虚拟一次波的零炮检距记录,可以看出构建的虚拟一次波和原始数据 (图 3) 近似,近炮检距的缺失对虚拟一次波的构建并未产生较大的影响.仔细对比可以发现,构建的虚拟一次波与原始数据仍然存在着差异,比如虚拟一次波中产生了随机的噪声,而且振幅也不同于原始记录,并且在19 km附近出现了不准确的倾角信息.图 6a为虚拟一次波的单炮记录,可以看出近炮检距同相轴倾角和连续性与原始数据相近,但在远炮检距的预测结果并不理想,可是远炮检距数据并不需要重建,因此该部分数据并不对近炮检距数据重建产生影响.如果将虚拟一次波直接替代缺失数据 (图 6b),可以看出不仅引入了噪声,而且同相轴也不具有连贯性.
为了提高插值结果,利用本文的插值方法,将构建的虚拟一次波作为初始数据,首先利用带有整形正则化的最小二乘反演方法 (公式 (3)) 估算自适应预测误差滤波器的系数,采用17(时间)×7(空间) 为滤波器尺寸参数,20(时间)×5(空间) 的平滑半径参数.然后利用最小二乘反演方法 (公式 (4)) 对缺失的数据进行插值处理.图 7a为本文方法重建之后单炮记录的结果,相比较于常规的预测误差滤波插值方法 (图 7b),同相轴的连续性得到了更好的恢复,尤其在5~8 s的数据更清晰地被体现出来.对重建结果进行验证,图 8a为单道波形对比图,可以看出插值之后的数据在振幅和相位上更接近原始数据.插值误差 (图 8b) 主要集中在散射波的位置.
接下来,测试三维自适应预测误差滤波器的应用,并结合最小二乘反演方法对地震数据进行重建.增加了震源点方向的空间变量,为地震数据的重建提供更多的有效信息.在三维非平稳滤波器情况下,采用7(时间)×7(空间)×7(空间) 为滤波器尺寸参数,20(时间)×5(空间)×10(空间) 的平滑半径参数,利用最小二乘反演方法 (公式 (4)) 对缺失的数据进行插值处理.图 9a为利用三维自适应预测误差滤波器数据重建的结果,共炮检距剖面中同相轴的连续性较好,主要的反射波信息被恢复,在单炮记录中也没有出现多余的同相轴.相比较于二维处理结果 (图 9b),由于噪声引起的抖动明显减少,空间方向的地震信息更平稳,5.5 s附近的同相轴信息在重建结果中得到了更清晰的体现 (多次波信息).
接下来,测试该方法对噪声的适应能力.在理论模型中加入随机噪声,信噪比约为3.7 dB (图 10a),在对近炮检距缺失15道的情况下 (图 10b),利用本文方法对缺失数据进行重建,选取相同的滤波器尺寸参数和平滑半径参数,从插值结果 (图 10c) 上看,与无噪声情况 (图 7a) 相比,在随机噪声没有覆盖有效信息的情况下,重建结果中2.5 s附近同相轴的信息可以得到恢复,并且具有较好的连续性,但5 s以下的同相轴信息在重建结果中受到噪声的影响.理论模型表明该方法对随机噪声具有一定的适用性.
利用墨西哥湾某地区的实际数据来进一步测试本文的方法.图 11为炮集数据体,图中正上部图为3.8 s位置的时间切片,左下图为位置在0.05 km的共炮检距剖面,右下图为位置在5 km的单炮记录.炮间距为0.05 km,最小炮检距0.075 km,最大炮检距3 km,共计250炮.从单炮数据中可以看出,在3.5 s左右存在较明显的表层多次波,但也存在着由于盐丘体形成的复杂散射波.
图 12为利用公式 (2) 计算出的虚拟一次波数据,其中的多次波通过Radon变换 (Foster and Mosher, 1992) 获取.相比于原始记录 (图 11),近炮检距信息得到了恢复,但振幅和相位特征不准确,互相关引起的随机噪声也比较明显,并且由于实际数据中包含复杂的一次波、多次波以及散射波信息,部分海底反射的倾角信息很难被呈现出来.
利用本文提出的插值方法,利用带有整形正则化的最小二乘反演方法估算三维自适应预测误差滤波器的系数,采用15(时间)×5(空间)×10(空间) 为滤波器尺寸参数,25(时间)×10(空间)×10(空间) 的平滑半径参数处理实际地震数据.然后利用最小二乘反演方法 (公式 (4)) 对缺失的数据进行插值处理,如图 13所示,随机噪声得到了明显的压制,浅层海底反射的同相轴得到了较好的插值.尤其是在不同时间深度上具有非平稳振幅关系的波场可以得到合理的恢复.
本文提出利用基于非平稳自回归的自适应预测误差滤波方法对近炮检距缺失的数据进行重建.该方法利用多次波中蕴含的动力学信息预测缺失的地震数据,自适应预测误差滤波器结合非平稳自回归可以在局部对振幅和相位进行校正,三维预测误差滤波器可以有效降低预测过程中引入的噪声,但是不可避免地增加了内存空间和计算时间.通过对模型数据和实际地震数据进行处理,验证了该方法可以合理地对近炮检距缺失数据进行重建.
附录A
假设D为线性算子,最小二乘估计变为线性反演:
(A1) |
其中:
(A2) |
(A3) |
矩阵A的元素为:
(A4) |
整形规则化 (Fomel,2007) 采用一个整形算子 (平滑算子)G代替D,并且相较于Tikhonov正则化 (Tikhonov, 1963) 有着更好的计算性能.利用整形正则化的反演表达式为:
(A5) |
其中
(A6) |
矩阵
(A7) |
其中λ为缩放系数.整形规则化方法的优势在于具有易于控制的参数λ,并且利用一个参数G代替ε和D.G为高斯平滑半径,人为设定参数 (Fomel,2007),而且定义系数λ为Sn(t, x) 的平均值.
致谢感谢两位匿名审稿专家提出的宝贵意见.
Claerbout J F. Earth Soundings Analysis: Processing Versus Inversion. London: Blackwell Scientific Publications, 1992. | |
Claerbout J F. 2010. Image estimation by example: Geophysical soundings image construction-Multidimensional autoregression. Stanford Exploration Project. | |
Crawley S, Clapp R, Claerbout J. 1999. Interpolation with smoothly nonstationary prediction-error filters.//SEG Technical Program Expanded Abstracts. SEG, 1154-1157. | |
Crawley S. 2000. Seismic trace interpolation with nonstationary prediction-error filters [Ph. D. thesis]. Stanford, California: Stanford University. | |
Curry W, Shan G J. 2010. Interpolation of near offsets using multiples and prediction-error filters. Geophysics, 75(6): WB153-WB164. DOI:10.1190/1.3506557 | |
Fomel S. 2002. Applications of plane-wave destruction filters. Geophysics, 67(6): 1946-1960. DOI:10.1190/1.1527095 | |
Fomel S. 2007. Shaping regularization in geophysical-estimation problems. Geophysics, 72(1): R29-R36. | |
Fomel S. 2009. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics, 74(1): V25-V33. DOI:10.1190/1.3043447 | |
Foster D J, Mosher C C. 1992. Suppression of multiple reflections using the Radon transform. Geophysics, 57(3): 386-395. DOI:10.1190/1.1443253 | |
Gao J J, Chen X H, Li J Y, et al. 2009. Study on reconstruction of seismic data based on nonuniform Fourier transform. Progress in Geophys., 24(5): 1741-1747. DOI:10.3969/j.issn.1004-2903.2009.05.026 | |
Guo S J, Fang W B, Xu Z T. 2013. Interpolation of near offset seismic trace using surface-related multiples. Geophysical Prospecting for Petroleum, 52(6): 636-641. | |
Guitton A, Brown M, Rickett J, et al. 2001. Multiple attenuation using a t-x pattern-based subtraction method.//SEG Technical Program Expanded Abstracts. SEG, 1305-1308. | |
Hestenes M R, Stiefel E. 1952. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6): 409-436. DOI:10.6028/jres.049.044 | |
Kabir M M N, Verschuur D J. 1995. Restoration of missing offsets by parabolic Radon transform. Geophysical Prospecting, 43(3): 347-368. DOI:10.1111/gpr.1995.43.issue-3 | |
Li X F, Li X F. 2008. Seismic data reconstruction with fractal interpolation. Chinese J. Geophys., 51(4): 1196-1201. | |
Liu X W, Liu H, Li Y M. 2004. High resolution radon transform and its application in seismic signal processing. Progress in Geophysics, 19(1): 8-15. | |
Liu Y, Fomel S. 2011. Seismic data interpolation beyond aliasing using regularized nonstationary autoregression. Geophysics, 76(5): V69-V77. DOI:10.1190/geo2010-0231.1 | |
Meng X H, Guo L H, Zhang Z F, et al. 2008. Reconstruction of seismic data with least squares inversion based on nonuniform fast Fourier transform. Chinese J. Geophys., 51(1): 235-241. | |
Shan G J, Guitton A. 2004. Migration of surface-related multiples: Tests on the Sigsbee2B datasets.//SEG Technical Program Expanded Abstracts. SEG, 1285-1288. | |
Spitz S. 1991. Seismic trace interpolation in the F-X domain. Geophysics, 56(6): 785-794. DOI:10.1190/1.1443096 | |
Tikhonov A N. 1963. Solution of incorrectly formulated problems and the regularization method. Soviet Math. Dokl., 4: 1035-1038. | |
van Dedem E J, Verschuur D J. 2005. 3D surface-related multiple prediction: A sparse inversion approach. Geophysics, 70(3): V31-V43. DOI:10.1190/1.1925752 | |
Wang W H, Pei J Y, Zhang J F. 2007. Prestack seismic data reconstruction using weighted parabolic Radon transform. Chinese J. Geophys., 50(3): 851-859. | |
Wang Y B, Schuster G T. 2007. Interferometric interpolation of missing seismic data.//SEG Technical Program Expanded Abstracts. SEG, 2688-2692. | |
Wang Y B, Dong S Q, Luo Y. 2010. Model-based interferometric interpolation method. Geophysics, 75(6): WB211-WB217. DOI:10.1190/1.3505816 | |
Xu S, Zhang Y, Pham D, et al. 2005. Antileakage Fourier transform for seismic data regularization. Geophysics, 70(4): V87-V95. DOI:10.1190/1.1993713 | |
Zhang H M, Liu H. 2006. Interpolation of poststack seismic traces based on sparse discrete τ-p transform. Oil Geophysical Prospecting, 41(3): 281-285. | |
高建军, 陈小宏, 李景叶, 等. 2009. 基于非均匀Fourier变换的地震数据重建方法研究. 地球物理学进展, 24(5): 1741–1747. DOI:10.3969/j.issn.1004-2903.2009.05.026 | |
郭书娟, 方伍宝, 徐兆涛. 2013. 基于表层多次波的地震数据插值方法研究. 石油物探, 52(6): 636–641. | |
李信富, 李小凡. 2008. 分形插值地震数据重建方法研究. 地球物理学报, 51(4): 1196–1201. | |
刘喜武, 刘洪, 李幼铭. 2004. 高分辨率radon变换方法及其在地震信号处理中的应用. 地球物理学进展, 19(1): 8–15. | |
孟小红, 郭良辉, 张致付, 等. 2008. 基于非均匀快速傅里叶变换的最小二乘反演地震数据重建. 地球物理学报, 51(1): 235–241. | |
王维红, 裴江云, 张剑锋. 2007. 加权抛物Radon变换叠前地震数据重建. 地球物理学报, 50(3): 851–859. | |
张红梅, 刘洪. 2006. 基于稀疏离散τ-p变换的叠后地震道内插. 石油地球物理勘探, 41(3): 281–285. | |