干涉测量在地震勘探上的应用始于Claerbout (1968),他首次给出了水平层状介质模型地震透射响应和反射响应之间的关系,即地震透射记录进行自相关等价于其自激自收的模拟记录的推论.之后将此推论应用到三维介质,并命名为“日光成像”(Daylight Imaging).Zhang等 (1989)运用分解平面波的方法证明了“日光成像”在三维均匀介质情况下的可行性.在稳相分析的基础上,Schuster等 (2000)对互相关地震记录进行偏移处理,完成了地下介质的反射波成像, 并在之后发表的文章中首次提出地震干涉测量这个名词.Wapenaar等 (2002)用格林定理证明了Claerbout的推论,为地震干涉测量提供了坚实的数学基础.Calvert等 (2004)提出了一种名为“虚震源”法,指出该方法可以有效避免近地表复杂速度建模过程,并强调“虚震源”法还可以应用于时移地震,监测油藏动态变化.针对相关型地震干涉测量存在的虚假同向轴、震源分布不均匀等问题,Wapenaar和van der Neut (2010)提出了多维反褶积地震干涉测量,该方法在弥补低频和高频能量的损失、改善恢复波场的稳定性和分辨率方面有良好效果.Bellezza和Poletto (2014)等改进了多维反褶积地震干涉测量,通过互褶积将多维反褶积地震干涉测量运用到虚反射方法上,应用Arctic实际资料验证了方法的可行性.综上,近些年国外学者对勘探地球物理地震干涉测量的研究主要集中在波场重建、静校正、裂缝预测和深层勘探等方面 (Duguid et al., 2011; Liu et al., 2013; Liu et al., 2014).
地震干涉测量在国内的研究起步较晚.陶毅等 (2010)对地震干涉测量的提出、发展和完善及其应用进行了总结归纳.罗桂纯等 (2011)利用地震干涉测量对地铁震动产生的环境噪声对钢筋混凝土建筑结构进行了重构分析.黄伟传等 (2012)结合胜利油田Fan158井随钻地震数据,开展了地震干涉测量偏移成像的研究,取得了较好的效果.王晨龙等 (2013)提出了地面与井中观测条件下的微地震干涉逆时定位算法.针对西部山地地震资料低信噪比的问题,徐基祥等 (2014)提出了地震干涉测量法近地表散射波分离技术.乔宝平等 (2014)提出了基于逆虚折射干涉法,有效提取近地表弱信号,成功恢复了折射波信息.叶月明等 (2015)提出了基于地震干涉测量的海底相关层间多次波衰减方法和实现流程,并以模型试算和实际资料验证了该方法的可行性.
稳相分析能够形象解释地震干涉测量的物理含义,因此被不少勘探地球物理专家使用.Snieder (2004)(Snieder et al., 2006) 利用稳相分析探讨了地面地震剖面 (SSP,Surface Seismic Profile) 观测系统地震干涉测量的物理含义,同时分析了地震干涉测量中虚假多次波问题.首先,本文结合Snieder的思路,首次探讨VSP观测系统的稳相分析,解释其地震干涉测量的物理含义,分析不同交叉项的含义,探讨虚假同向轴产生的原因;其次,利用地震干涉测量实现VSP波场向SWP波场转化,探测诸如断层这样的高陡倾角地质构造,偏移成像结果证明该方法的可行性;最后,将该成像结果与传统地面地震勘探在高陡地层成像质量方面进行对比,得出该方法可以有效提高高陡地层成像精度的结论.
1 地震干涉测量原理如图 1所示,在三维声波介质中,由格林定理、Rayleigh定理得到频率域内格林积分表达式为 (Wapenaar et al., 2005):
(1) |
式 (1) 中ℜ表示实部,G(XB, XA, ω) 表示A点激发、B点接收的波场,G(XA, X, ω) 表示X点激发、A点接收的波场,G(XB, X, ω) 表示X点激发、B点接收的波场,ω为圆频率;等式右边第一个积分符号是面积分,积分区域为震源所包围的区域S,n表示外法线方向;j表示单位虚部,ρ(X) 声波介质密度,∂i表示偏导数,i取1、2、3,*表示复共轭.式 (1) 右边积分式子表示频率域内一格林函数与另一格林函数的共轭相乘,经傅里叶反变换可以得到时间域两检波点波场互相关的形式,因此上式也称为互相关型地震干涉测量.这里的波场指全波场,即包含直达波、反射波以及散射波等.式 (1) 表明,对A、B两接收点波场进行相关处理,可以得到以其中一接收点作为新震源另外一点作为接收点的新的地震波场.
式 (1) 包含对空间方向导数的求导,实际资料很难获取.在远场近似情况下,上式可以写为 (Schuster,2009):
(2) |
近似式 (2) 可能会因为省略了部分项引起稳相点的丢失从而导致虚假同向轴产生,但当震源和检波点之间距离足够大时,上式可认为是准确的.由于式 (2) 不存在空间方向导数的求取,远场表达式为地震干涉测量在实际资料中的应用提供了便利.
2 稳相分析稳相分析是一种近似的方法,解决下面的积分函数 (Bleistein,1984),公式为
(3) |
式 (3) 积分为实数域积分,g(x) 为一缓慢变化的函数,ϕ(x) 是实函数,称为相函数,λ表示渐近系数.若λϕ(x) 关于x变化剧烈,那么eiλϕ(x)就是一个快速震荡的函数,式 (3) 在实数域内的积分等于0.然而,当存在x=x*时相位函数满足ϕ(x)′x=x*=0,上式的积分主要受 (x*-ε, x*+ε) 邻域控制,此时x*被称为稳相点,相应的邻域被称为稳相区.积分式 (3) 主要贡献来自稳相点,所以f(λ) 关于x*一阶稳相泰勒近似可以表达为
(4) |
式 (4) 表明,积分表达式的主要贡献来自稳相点.
3 VSP模型的稳相分析本文建立如图 2所示的模型,来探讨VSP波场向SWP波场转化的物理含义.震源位于地表,在均匀介质中,有一断层反射面 (或者其他高陡波阻抗界面),断面倾角为θ0,反射系数为r,井中放置A、B两个检波器,道间距z0,检波器A到断面的距离为d,边界为吸收边界.
频率域内地震波场的传播可以表示为
(5) |
式中,ρ表示介质密度,W(ω) 表示震源子波,ω表示圆频率,k表示波数,R表示地震波传播的距离.定义LAS、LBS为震源S到A、B两个检波器的距离,RA、RB表示两反射点,那么检波器A、B的全波场可以表示为
(6) |
(7) |
那么频率域互相关检波器A、B波场,得:
(8) |
互相关结果包含四项,即T1、T2、T3、T4四项.鉴于稳相分析的直观性,下面利用该方法对对每一项进行解释.
如图 3所示,第一项T1表示A、B两直达波互相关, 公式为
(9) |
由稳相分析可知,积分表达式的主要贡献取决于稳相点,当x满足:
(10) |
即存在θA=θB,当检波器A接收到的直达波场入射角与检波器B接收到直达波场入射角相同时,存在稳相点.此时A、B两波场直达波有一部分重合,互相关后相互重合的路径抵消,得到的结果正好是以检波器B作为新震源、A作为新的接收点接收到的直达波波场.如图 3b所示,稳相点位于井口附近,对于零偏移距VSP资料,A、B直达波波场互相关后存在稳相点,能够准确重构新的直达波场;对于非零偏移距VSP资料,A、B直达波波场互相关后不存在稳相点.因此,T1项存在稳相点,在零偏移距情况下能够精确重构直达波场,理论上该项不会产生虚假同向轴.然而,重构的新的直达波场在资料处理过程中用处不大.
如图 4所示,第二项T2表示A点的反射波与B点的直达波互相关,公式为
(11) |
积分表达式的主要贡献取决于稳相点,当x满足:
(12) |
即存在θA=θB,当检波器A接收到的反射波场入射角与检波器B接收到直达波场入射角相同时,存在稳相点.此时A、B两点波场路径一部分重合,互相关后相互重合的路径抵消,得到的结果正好是以检波器B作为新震源、A接收到来自断面的反射波波场.重构的新的反射波场是以B点为激发点,A点作为检波点,相当于把观测系统重构到井中,实现了震检基准面的转化,可以有效避免因地表起伏引起的震源静校正问题.一方面,重构的新震源不是物理存在的,不会对井产生任何破坏,因此有的专家学者称之为“虚震源”.另一方面,虚震源相比于地面震源而言离井附近的目标体更近,可以有效探测地质体信息,提高成像精度,尤其对于高陡倾角的地质体.
如图 5所示,第三项T3表示A点的直达波与B点的反射波互相关,公式为
(13) |
由稳相分析可知,积分表达式的主要贡献取决于稳相点,当x满足:
(14) |
即不存在θA=θB,在地面找不到这样的震源满足当检波器A接收到的直达波场入射角与检波器B接收到反射波场入射角相同的情况.因检波器A位于检波器B的下方,A接收到地面震源的直达波场不可能和B接收到来自断面反射波场有相同的路径,因此不存在稳相点.数值模拟表明,尽管T3不存在稳相点,但在波场重构中不会产生虚假同向轴.
如图 6所示,第四项T4表示A点的反射波与B点的反射波互相关,公式为
(15) |
由稳相分析可知,积分表达式的主要贡献取决于稳相点,当x满足:
(16) |
即不存在θA=θB,在地面找不到这样的震源,满足当检波器A接收到的反射波场入射角与检波器B接收到反射波场入射角相同的情况.因检波器A位于检波器B的下方,A接收到来自断面反射波场不可能和B接收到来自断面反射波场有相同的路径,因此不存在稳相点.由于检波器A的反射波时间延迟往往大于检波器B反射波的时间延迟,互相关后积分T4不为0,会对整个积分结果产生影响.因此,T4不存在稳相点,但在波场重构中会产生虚假同向轴.在应用地震干涉测量处理时,这项应尽量压制,避免虚假同向轴带来的干扰.
通过以上分析可知,VSP波场向SWP波场转化的准确性主要取决于稳相点.积分项T1、T2存在稳相点,是地震干涉测量的主要贡献项;在稳相点附近通过T1可以准确重构直达波场,T2可以准确重构反射波场.T3、T4不存在稳相点,是地震干涉测量的干扰项.实际应用中,全波场地震干涉测量非常便捷,但全波场地震干涉测量处理后会产生干扰项,形成虚假同向轴,影响波场转化的准确性,对于复杂的地质构造,全波场地震干涉测量后干扰项会增多.目前,为压制VSP波场向SWP波场转化过程中产生的虚假同向轴,在互相关处理之前可进行波场分离,用中值滤波、F-K滤波和τ-p变换等方法实现VSP波场的上下波场分离,针对不同的地质体利用不同的波场进行地震干涉测量处理,可以大大提高重构的SWP波场的质量,有效压制虚假同向轴.
4 模型测试波场重构是地震干涉测量的一个重要应用,为验证地震干涉测量能够有效实现VSP波场和SWP波场的转化,建立如图 7a所示的二维声波模型.模型包括两水平地层,速度分别为1500 m/s、2000 m/s,高陡隆起速度为3500 m/s.地表设置101个震源,炮间距14 m,井中101个检波器,道间距20 m,正演模拟结果如图 9所示.
地震干涉测量实现VSP波场向SWP波场转化的过程分为以下四步:(1) 波场分离,利用中值滤波或者τ-p变换提取VSP的上行直达波;(2) 互相关处理,选取一炮,将分离出第一道的直达波场与第二道的反射波场做互相关;(3) 叠加处理,将每炮互相关处理的结果进行叠加,得到以第一个检波器为震源、第二个检波器为接收点的地震记录;(4) 重复1~3过程,实现不同道之间的处理,得到重构的SWP剖面.图 8a表示通过地震干涉处理重构的SWP剖面图,震源以第40个检波器为“新”的震源,即虚震源,其他检波器为“新”的接收点所接收到的地震波场,波场剖面包括直达波和来自断面的反射波.图 8b表示直接将震源放置在第40个检波器位置所接收到的理论地震波场.对比可知,无论是直达波场还是反射波场,重构的SWP波场与理论波场基本一致,能够有效实现VSP波场向SWP波场的转化,验证了该方法的可行性.
重构的SWP波场有效恢复了来自高陡隆起界面的反射波,对其建立局部速度模型,进行偏移成像,结果如图 9b.成像结果真实反应了高陡隆起界面,高陡界面的形态和位置都能清晰地展现出来,说明该方法通过波场转化,将位于地表的震源转换到井中,能够有效探测高陡波阻抗界面.然而偏移结果不能正确成像水平层,说明对水平层状介质成像效果不好,但这不是该方法关注的重点,因为传统的地面地震勘探或垂直地震勘探已经能够很好的探测小角度地层.因此直接利用重构的SWP波场,可以准确成像高陡反射界面,提高了横向勘探的成像精度.
建立常规地面地震勘探模型,如图 10a所示,用传统地面地震勘探成像结果与该方法进行对比,验证该方法可以有效提高高陡构造的成像精度.模型中,地表设置101个震源和101个检波器.从图 10b的偏移成像剖面可知,水平反射界面成像质量好,轮廓清晰;然而右侧的高陡反射界面成像却非常差,出现虚假同向轴,界面能量弱,不能准确的反映高陡界面的基本形态.对比图 9b与10b的偏移结果,我们可以得出这样的结论:常规地面地震勘探对水平层的分辨率高,对高陡地层的分辨率低;通过地震干涉测量实现VSP波场向SWP波场转化,直接利用SWP波场成像,可以有效提高高陡构造的分辨能力,提高成像精度.
本文利用稳相分析解释了地震干涉测量实现VSP波场向SWP波场转化的物理含义,并通过二维声波模型实现了VSP波场向SWP波场的转化,实现了高陡构造的成像,并将成像结果与常规地面地震勘探成像结果进行对比.本文研究得出如下结论:
(1) VSP波场向SWP波场转化的准确性主要取决于稳相点.积分项T1、T2存在稳相点,是地震干涉测量的主要贡献项;在稳相点附近通过T1可以准确重构直达波场,T2可以准确重构反射波场.T3、T4不存在稳相点,是地震干涉测量的干扰项;T3不存在稳相点,理论上不会产生虚假同向轴,T4不存在稳相点,会虚假同向轴,处理过程中应尽量压制.
(2) VSP波场向SWP波场的转化,实现了震检基准面变换.传统VSP波场是在地表激发井中接收,重构的新的SWP波场是在井中激发井中接收,这样转化将震源重构到井中,可以有效避免第一个检波器以上复杂盖层速度模型的影响,自动消除因地表起伏引起的震源静校正问题.
(3) 重构的SWP波场离高陡目标体更近,建立局部速度模型,直接利用重构的SWP波场对深部高陡目标成像,可以有效提高高陡构造横向勘探精度.与常规地面地震勘探相比,该方法具有提高横向探测精度的优势,可以成像岩丘侧翼、高陡隆起、大角度断层等高陡构造,为长期难以解决的复杂高陡构造成像问题提供了一种新的思路.
致谢 感谢国家油气重大专项 (2016ZX05003) 对本文研究的资助,感谢审稿专家提出的宝贵意见和编辑部的大力支持![] | Bellezza C, Poletto F. 2014. Multidimensional deconvolution and processing of seismic-interferometry Arctic data[J]. Geophysics, 79(3): WA25–WA38. DOI:10.1190/geo2013-0297.1 |
[] | Bleistein N. 1984. Mathematical methods for wave phenomena[M]. New York: Academic Press: 77-82. |
[] | Calvert R W, Bakulin A, Jones T C. 2004. Virtual sources, a new way to remove overburden problems[C].//Proceedings of the 66th Mtg., European Association of Geoscientists and Engineers.Extended Abstracts, 234. |
[] | Claerbout J F. 1968. Synthesis of a layered medium from its acoustic transmission response[J]. Geophysics, 33(2): 264–269. DOI:10.1190/1.1439927 |
[] | Duguid C, Halliday D, Curtis A. 2011. Source-receiver interferometry for seismic wavefield construction and ground-roll removal[J]. The Leading Edge, 30(8): 838–843. DOI:10.1190/1.3626489 |
[] | Huang W C, Ge H K, Wu H Z, et al. 2012. Seismic interferometriy in seismic whiling drilling data processing[J]. Oil Geophysical Prospecting (in Chinese), 47(1): 32–36. |
[] | Liu L B, Zhao Z, Zhao J G. 2013. Characterization of fractures using interferometry techniques with multi-hole, multi-polarization borehole radar[C].//SEG Annual Meeting, Society of Exploration Geophysicists. Expanded Abstracts, 4508-4515. |
[] | Liu Y, Arntsen B, Wapenaar K, et al. 2014. Combining inter-source seismic interferometry and source-receiver interferometry for deep local imaging[C].//SEG Annual Meeting, Society of Exploration Geophysicists. Expanded Abstracts, 5107-5112. |
[] | Luo G C, Liu L B, Qi C, et al. 2011. Structural response analysis of a reinforced concrete building based on excitation of microtremors and passing subway trains[J]. Chinese J. Geophys.(in Chinese), 54(10): 2708–2715. DOI:10.3969/j.issn.0001-5733.2011.10.028 |
[] | Qiao B P, Guo P, Wang P, et al. 2014. Effectively picking weak seismic signal near the surface based on reverse virtual refraction interferometry[J]. Chinese J. Geophys.(in Chinese), 57(6): 1900–1909. DOI:10.6038/cjg20140621 |
[] | Schuster G T, Rickett J. 2000. Daylight imaging in V (x, y, z) media[R]. Stanford Exploration Project, Report105. 55-56. |
[] | Schuster G T. 2009. Seismic Interferometry[M]. New York: Cambridge University Press: 38-40. |
[] | Snieder R. 2004. Extracting the Green's function from the correlation of coda waves:A derivation based on stationary phase[J]. Phys. Rev. E, 69(4): 046610. DOI:10.1103/PhysRevE.69.046610 |
[] | Snieder R, Wapenaar K, Larner K. 2006. Spurious multiples in seismic interferometry of primaries[J]. Geophysics, 71(4): SI111–SI124. DOI:10.1190/1.2211507 |
[] | Tao Y, Fu L Y, Sun W J, et al. 2010. A review of seismic interferometry[J]. Progress in Geophysics (in Chinese), 25(5): 1775–1784. DOI:10.3969/j.issn.1004-2903.2010.05.035 |
[] | Wang C L, Cheng J B, Yin C, et al. 2013. Microseismic events location of surface and borehole observation with reverse-time focusing using interferometry technique[J]. Chinese J. Geophys.(in Chinese), 56(9): 3184–3196. DOI:10.6038/cjg20130931 |
[] | Wapenaar K, Draganov D, Thorbecke J, et al. 2002. Theory of acoustic daylight imaging revisited[C].//Proceedings of the 72nd Annul International Meeting of the Society of Exploration Geophysicists.Expanded Abstracts, 69-72. |
[] | Wapenaar K, Fokkema J, Snieder R. 2005. Retrieving the Green's function in an open system by crosscorrelation:A comparison of approaches[J]. Acoustical Society of America, 118(5): 2783–2786. DOI:10.1121/1.2046847 |
[] | Wapenaar K, van der Neut J. 2010. A representation for Green's function retrieval by multidimensional deconvolution[J]. The Journal of the Acoustical Society of America, 128(6): EL366–EL371. DOI:10.1121/1.3509797 |
[] | Xu J X. 2014. Separating the near-surface seismic scattered wave using seismic interferometry method[J]. Chinese J. Geophys.(in Chinese), 57(6): 1910–1923. DOI:10.6038/cjg20140622 |
[] | Xu J X, McLean B F, Song X J. 2014. A near-surface seismic scattered wave separation method[J]. Petroleum Exploration and Development (in Chinese), 41(6): 705–711. |
[] | Ye Y M, Yao G S, Zhao C L, et al. 2015. Sea bottom peg-leg multiple suppression based on seismic interferometry[J]. Oil Geophysical Prospecting (in Chinese), 50(2): 225–231. |
[] | Zhang L. 1989. Reflectivity estimation from passive seismic data[J]. Stanford Exploration Project-Annual Report, 60: 85–96. |
[] | 黄伟传, 葛洪魁, 吴何珍, 等. 2012. 地震干涉处理方法在随钻地震资料处理中的应用[J]. 石油地球物理勘探, 47(1): 32–36. |
[] | 罗桂纯, 刘澜波, 齐诚, 等. 2011. 基于地脉动和地铁振动的钢筋混凝土建筑结构响应分析[J]. 地球物理学报, 54(10): 2708–2715. DOI:10.3969/j.issn.0001-5733.2011.10.028 |
[] | 乔宝平, 郭平, 王璞, 等. 2014. 基于逆虚折射干涉法有效提取近地表弱地震信号[J]. 地球物理学报, 57(6): 1900–1909. DOI:10.6038/cjg20140621 |
[] | 陶毅, 符力耘, 孙伟家, 等. 2010. 地震波干涉法研究进展综述[J]. 地球物理学进展, 25(5): 1775–1784. DOI:10.3969/j.issn.1004-2903.2010.05.035 |
[] | 王晨龙, 程玖兵, 尹陈, 等. 2013. 地面与井中观测条件下的微地震干涉逆时定位算法[J]. 地球物理学报, 56(9): 3184–3196. DOI:10.6038/cjg20130931 |
[] | 徐基祥. 2014. 地震干涉测量法近地表散射波分离技术[J]. 地球物理学报, 57(6): 1910–1923. DOI:10.6038/cjg20140622 |
[] | 徐基祥, McLeanB F, 宋雪娟. 2014. 近地表地震散射波分离方法[J]. 石油勘探与开发, 41(6): 705–711. |
[] | 叶月明, 姚根顺, 赵昌垒, 等. 2015. 利用地震干涉法衰减海底相关层间多次波[J]. 石油地球物理勘探, 50(2): 225–231. |