2. 黑龙江省普通高校科技创新团队"断层变形、封闭性及流体运移", 大庆 163318
2. Science and Technology Innovation Team on Fault Deformation, Sealing and Fluid Migration, Daqing 163318, China
弹性波逆时偏移是一种以矢量波理论为基础的深度域偏移方法,相比于声波逆时偏移考虑了地震波在地下传播过程中的耦合特征,波场特征更加全面,算法与实际情况更为接近.多分量地震数据的弹性波逆时偏移能够形成多波模式的成像剖面,提供了额外的地下构造信息,从而减少纵波勘探的多解性.假设忽略地层衰减,由于转换波的波长比纵波更短,转换波的成像比反射纵波成像分辨率更高(Yan and Sava, 2008;Zhang and McMechan, 2011;Wang and MeMchan, 2015).因此对多分量地震资料的弹性波逆时偏移成像方法展开研究,更具实际意义.
逆时偏移的思想最早由Whitemore在1983年的SEG年会上提出,此后,一些学者(Baysal et al., 1983;Loewenthal and Mufti, 1983;MeMechan, 1983)先后将它运用到纵波资料的叠后偏移,随着计算机硬件技术的发展,逆时偏移技术已经从叠后发展到叠前,从二维到三维,从声波到弹性波,从各向同性到各向异性(刘国峰等,2009;刘红伟等,2010;王保利等,2013).为了充分利用多波多分量弹性资料的矢量特征,减少纵波勘探的多解性,国内外基于双程弹性波方程的叠前逆时偏移已经进行了一些研究.Chang和McMechan(1987, 1994)基于激发时间成像条件,先后实现了2D和3D弹性波逆时偏移.杜启振和秦童(2009)基于激发时间成像条件实现了横向各向同性介质中的多波多分量叠前逆时偏移.张智等(2013)采用激发振幅成像条件实现了多波多分量资料的弹性波逆时偏移.但是,他们应用矢量波场的分量得到的成像结果,存在未分离纵波和横波场而引起的串扰(crosstalk),物理意义不明确.Sun和McMechan(2001)提出了标量的弹性波逆时偏移方法,该方法先在地面将多分量地震记录分离,然后利用纵波和横波的地震记录分别进行标量波动方程偏移,它解决了应用矢量分量成像产生的串扰问题,但不能准确的处理纵波和横波的能量转换关系,放弃了弹性波场的矢量特征(Sun and McMechan, 2001;Sun et al., 2011;Du et al., 2012).
目前,常用的一种方法是在成像前应用旋度和散度算子分离延拓的矢量波场,然后应用互相关成像条件形成成像结果(Yan and Sava, 2008;Du et al., 2012;Duan and Sava, 2015;Li et al., 2016).但是,这类方法的转换波成像结果会产生极性反转的问题,影响叠加剖面.Du等(2012)利用震源和检波点的Poynting矢量计算符号因子,进一步利用符号因子校正了PS和SP成像剖面极性.Duan和Sava(2015)利用纵横波传播方向、界面方位和横波的偏振之间的关系,校正了PS和SP成像极性反转问题.Choi等(2016)利用震源波场中S波的符号校正了PS和SP成像的极性,并压制了成像噪声.Li等(2016)借助弹性波解耦方程的极化矢量计算符号因子,进一步利用符号因子校正PS和SP的极性反转.另一种基于角度域校正方法更为精确(Rosales and Rickett, 2001;Rosales et al., 2008),但此类方法的计算量较大.由于应用旋度、散度分离波场会产生相位和振幅变化,使波场失去了物理意义,若该方法要实现保幅的偏移还需要额外的校正(Sun and McMechan, 2001;Sun et al., 2011;Li et al., 2016).
最近,矢量分离纵横波场弹性波方程的逆时偏移被提出(Gu et al., 2015;Wang and MeMchan, 2015;Gong et al., 2016;Du et al., 2017),该类方法优于传统散度和旋度算子分离矢量波场的RTM方法,因为矢量分离保持了波场的振幅和相位特征.然而,分离得到的纵波和横波场都是矢量分量形式,若对纵波和横波矢量分量之间直接应用互相关成像(Gu et al., 2015),成像结果与传统意义上的反射系数不同,物理意义不明确给解释带来困难.为了解决该问题,将矢量波场标量化,通过分配合理的极性,能够获得PP和PS反射率成像,然而该方法需要传播角度已知,当涉及多路径时,这类方法变得不合理(Wang and MeMchan, 2015;Gong et al., 2016).
基于上述认识,本文利用矢量分离纵横波场的一阶速度-应力方程组构建矢量的纵波和横波波场,避免了分离波场的振幅和相位变化(马德堂和朱光明,2003).首先通过高阶交错网格对速度-应力方程组离散得到构建弹性波震源和检波点波场的差分算子,然而,利用均匀模型测试验证矢量分离的有效性.对于矢量分离的纵波和横波波场,本文提出应用震源归一化的内积成像条件,避免了纵波和横波矢量分量之间互相关成像结果的物理意义不清晰,而且转换波成像结果自动避免了传统方法的极性反转.
1 方法原理 1.1 弹性波矢量分离波场的速度-应力方程基于矢量分离纵横波的弹性波逆时偏移包括三个部分:1) 构建纵波和横波的震源矢量波场;2) 将多分量地震记录作为速度-应力方程的边界条件,构建纵波和横波的检波点矢量波场;3) 应用内积成像条件.在各向同性介质,构建纵波和横波的矢量波场的速度-应力方程如式(1) 所示,公式为
(1) |
式中,vx和vz分别为质点的水平分量和垂直分量的振动速度;vxp和vzp分别为质点的纵波水平分量和垂直分量的振动速度;vxs和vzs分别为质点的横波水平分量和垂直分量的振动速度;σxx和σzz分别为水平和垂直方向的正应力;τxz为切应力;ρ为介质的密度;λ和μ为介质的拉梅常数;t为时间变量.
1.2 高阶交错网格的有限差分格式根据高阶交错网格有限差分的定义,可以得到速度-应力方程(1) 的离散形式为
(2) |
其中,Dx和Dz为空间高阶差分算子,公式为
(3) |
(4) |
(5) |
式中,Δx和Δz分别表示沿和方向的空间网格大小;N等于空间差分阶数的一半;CnN为差分系数.本文采用时间2阶、空间12阶的差分格式,人工边界条件采用完美匹配层(PML)压制来自边界的反射.
1.3 成像条件与声波RTM不同,在弹性介质中,矢量分离得到的纵波和横波波场均为矢量分量形式,若直接对纵波和横波的矢量分量之间应用互相关成像条件,成像结果的物理意义与传统意义上的反射系数不同.因此,本文提出对分离的矢量纵波和横波波场应用震源归一化的内积成像条件.公式为
(6) |
式中,·代表震源和检波点纵横矢量波场之间的内积;S和R分别代表震源和检波点波场;(x, z)代表空间位置;nt采样时间总数目;下标i, j均代表分离后纵、横波.应用该成像条件,能够得到不同纯波模式之间的成像剖面,如:Ipp成像代表震源波场的矢量P波波场和检波点矢量P波波场的内积并归一化的结果.震源归一化成像条件它有效地消除了震源强度对成像结果的影响,成像振幅更接近于真实的反射系数.另外,还可以有效地补偿由于地震波在地层中传播时产生的能量损失,有利于深部地层成像.
1.4 震源波场逆时重建文中的弹性波逆时偏移应用内积成像条件时,与声波逆时偏移使用互相关成像条件相似,它需要同一时刻的震源波场和检波点矢量波场值,但是二者的传播方向相反,因此最直接的方法就是存储任意一方所有时刻的成像时所需要的波场值,当另一方传播时,提取存储的波场应用成像条件,但是,这会产生巨大存储花费和I/O负担.为了缓解该问题,我们利用额外的一次正演计算代替存储某一方向所有波场的方法(Feng and Wang, 2012;王保利等, 2012),也就是利用保存在边界内若干层的波场逆时重建震源波场的作法应用到本文的方法中.
2 波场分离测试为了测试文中纵横波场矢量分离的有效性,设计了一个均匀各向同性弹性介质模型,模型的大小为400×400网格,两个方向的网格尺寸均为10 m,纵波速度为4000 m/s,横波速度为2312 m/s,密度为1000 kg/m3,四周的PML边界厚度均为50个网格,时间采样间隔为1 ms,在模型中央位置放置混合震源,采用了主频为25 Hz的雷克子波作为激发震源.
图 1中显示了t=0.50 s时的波场快照及波场分离结果.图 1a、e分别为水平和垂直速度波场快照;图 1b、c、e和g分别是根据公式1矢量分离的纵波和横波波场快照;图 1d和h分别是对速度矢量波场应用散度和旋度算子分离的纵波和横波分量波场快照.从图 1a和e可以看出,纵波和横波耦合在速度分量之上,根据Helmholtz定理,应用旋度和散度算子能够很好的分离矢量的波场,如图 1d和1h所示,但是,这样分离的波场使输入波场的振幅和相位均发生了改变.利用文中所述了矢量分离纵横波场方程,纵波和横波场能够很好的分离开来,并且波场的矢量特征得到保存,振幅和相位均未发生改变,如图 1b、c、e和g所示.
对图 1的波场快照抽取单道,抽取的位置如图 1a所示虚线,各波场快照的单道信息如图 2所示.对比传统方法和矢量分离的方法,可以看出矢量分离的纵横波场保存了矢量波场的振幅和相位信息,而应用散度、旋度算子分离波场的传统方法改变了波场输入的矢量信息,使波场失去了物理意义.
为了测试文中所述的基于矢量分离纵横波的弹性波逆时偏移方法的有效性.我们测试了层状介质和岩丘模型.首先利用层状介质说明应用内积成像条件形成的成像结果中,转换波成像不会形成传统方法的极性反转,因此能够有效地避免极性反转对叠加剖面的影响.然后,测试复杂的岩丘模型,进一步验证文中所述方法的有效性.
3.1 层状模型为了清楚地说明内积成像结果不会形成极性反转现象,我们利用传统的旋度、散度分离波场的互相关成像结果进行对比.如图 3为层状介质模型参数,爆炸震源采用主频25 Hz的Ricker子波,共41炮,均匀分布在地面,横向和纵向网格间距均为10 m,采样时间为1 ms,400个检波点均匀分布在地面.图 4a和4b分别为炮点在2 km位置处的单炮成像结果,图 4c和图 4d为41炮叠加的成像结果.图 4a显示了传统方法单炮的PS成像结果,我们能够看见在地层纵波垂直入射点(图 4a中箭头位置)的左右两侧同相轴极性反转,因此,若不校正极性而让多炮成像结果叠加,结果会影响叠加剖面的同相轴,如图 4c所示,为传统方法叠加的PS成像结果,如箭头所示位置同相轴已经发生了扭曲.图 4b为本文内积成像的单炮PS成结果,在垂直入射位置的同相轴并没有出现极性反转,多炮叠加后的结果如图 4d所示,同相轴连续无错段,这说明文中所述的内积成像条件能有效地避免了转换波成像的极性反转.
由于在地层水平时,利用矢量分离得到的震源波场的P波矢量和检波点波场的S波矢量极性分布是相同的,因此应用内积成像条件形成的PS成像结果不会发生极性反转.在2D各向同性介质中,考虑纵波入射时,传统方法应用散度算子分离矢量波场得到的P波分量为标量,旋度算子分离矢量波场得到S波分量的极性在P波垂直入射两侧发生极性反转,因此传统方法直接对分离的P和S波应用互相关成像,PS成像结果会发生极性反转.类似地,当横波为入射波时,旋度算子分离的SP波发生极性反转,传统方法的SP成像结果也会产生极性反转.而文中所述的基于矢量分离纵波和横波波场的弹性波逆时偏移,形成的成像结果有效地避免了极性的反转.
4.3 岩丘模型进一步验证本文方法对复杂构造成像的有效性,我们测试了复杂的岩丘模型.岩丘纵波速度模型如图 5,横波速度
图 6a、b、d和f分别为应用内积成像条件的PP、SS、PS和SP成像结果.图 6c和e分别为传统方法的PS和SP成像结果.观察可知,内积成像条件的PP和PS剖面能够非常清晰地对岩丘和深层的断层成像,同相轴连续无错段.传统成像方法的结果,受极性反转的影响,如图 6e和6c箭头头所示位置,在岩丘底部同相轴错段,能量变弱,叠加剖面受到巨大的损害.应用本文的内积成像条件,PS和SP成像结果同相轴连续,无错段.由于在数值模拟时,采用爆炸源,转换波能量由纵波转换而来,能量相对较弱,所以内积成像条件的SP和SS成像结果不像PP和PS成像那么清晰.
本文提出了一种矢量分离纵横波场的弹性波逆时偏移方法,形成了PP、PS、SP和SS四种模式的高精度纯波成像剖面.相比于应用旋度、散度算子分离矢量波场的传统方法,文中提出的矢量分离方法较好地保持了原始波场的振幅和相位特征,转换波PS和SP成像无极性反转.文中提出的震源归一化内积成像条件形所成的纯波成像剖面能够准确地对复杂地下构造准确成像.本文的方法可进一步推广至三维弹性波逆时偏移高精度成像和多参数的全波形反演.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持![] | Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration[J]. Geophysics, 48(11): 1514–1524. DOI:10.1190/1.1441434 |
[] | Chang W F, McMechan G A. 1987. Elastic reverse-time migration[J]. Geophysics, 52(10): 1365–1375. DOI:10.1190/1.1442249 |
[] | Chang W F, McMechan G A. 1994. 3-D elastic prestack, reverse-time depth migration[J]. Geophysics, 59(4): 597–609. DOI:10.1190/1.1443620 |
[] | Chattopadhyay S, McMechan G A. 2008. Imaging conditions for prestack reverse-time migration[J]. Geophysics, 73(3): S81–S89. DOI:10.1190/1.2903822 |
[] | Choi H, Seol S J, Byun J. 2016. Converted-wave guided imaging condition for elastic reverse time migration with wavefield separation[J]. Exploration Geophysics. DOI:10.1071/EG16003 |
[] | Du Q Z, Guo C F, Zhao Q, et al. 2017. Vector-based elastic reverse time migration based on scalar imaging condition[J]. Geophysics, 82(2): S111–S127. DOI:10.1190/geo2016-0146.1 |
[] | Du Q Z, Qin T. 2009. Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium[J]. Chinese J. Geophys. (in Chinese), 52(3): 801–807. |
[] | Du Q Z, Zhu Y T, Ba J. 2012. Polarity reversal correction for elastic reverse time migration[J]. Geophysics, 77(2): S31–S41. DOI:10.1190/geo2011-0348.1 |
[] | Duan Y T, Sava P. 2015. Scalar imaging condition for elastic reverse time migration[J]. Geophysics, 80(4): S127–S136. DOI:10.1190/geo2014-0453.1 |
[] | Feng B, Wang H Z. 2012. Reverse time migration with source wavefield reconstruction strategy[J]. J. Geophys. Eng., 9(1): 69–74. DOI:10.1088/1742-2132/9/1/008 |
[] | Gong T, Nguyen B D, McMechan G A. 2016. Polarized wavefield magnitudes with optical flow for elastic angle-domain common-image gathers[J]. Geophysics, 81(4): S239–S251. DOI:10.1190/geo2015-0518.1 |
[] | Gu B L, Li Z Y, Ma X N, et al. 2015. Multi-component elastic reverse time migration based on the P-and S-wave separated velocity-stress equations[J]. Journal of Applied Geophysics, 112: 62–78. DOI:10.1016/j.jappgeo.2014.11.008 |
[] | Li B, Liu H W, Liu G F, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J]. Chinese J. Geophys. (in Chinese), 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
[] | Li Z Y, Ma X N, Fu C, et al. 2016. Wavefield separation and polarity reversal correction in elastic reverse time migration[J]. Journal of Applied Geophysics, 127: 56–67. DOI:10.1016/j.jappgeo.2016.02.012 |
[] | Liu G F, Liu H, Li B, et al. 2009. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation[J]. Chinese J. Geophys. (in Chinese), 52(12): 3101–3108. DOI:10.3969/j.issn.0001-5733.2009.12.019 |
[] | Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J]. Chinese J. Geophys. (in Chinese), 53(7): 1725–1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
[] | Loewenthal D, Mufti I R. 1983. Reversed time migration in spatial frequency domain[J]. Geophysics, 48(5): 627–635. DOI:10.1190/1.1441493 |
[] | Ma D T, Zhu G M. 2003. Numerical modeling of P-wave and S-wave separation in elastic wavefield[J]. Oil Geophysical Prospecting (in Chinese), 38(5): 482–486. |
[] | Rosales D, Rickett J. 2001. PS-wave polarity reversal in angle domain common-image gathers[C].//71st Annual International Meeting, SEG. Expanded Abstracts, 1843-1846. |
[] | Rosales D A, Fomel S, Biondi B, et al. 2008. Wave-equation angle-domain common-image gathers for converted waves[J]. Geophysics, 73(1): S17–S26. DOI:10.1190/1.2821193 |
[] | Sun R, McMechan G A. 1986. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles[J]. Proc. IEEE, 74(3): 457–465. DOI:10.1109/PROC.1986.13486 |
[] | Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data[J]. Geophysics, 66(5): 1519–1527. DOI:10.1190/1.1487098 |
[] | Sun R, Chow J, Chen K J. 2001. Phase correction in separating P-and S-waves in elastic data[J]. Geophysics, 66(5): 1515–1518. DOI:10.1190/1.1487097 |
[] | Sun R, McMechan G A, Chuang H H. 2011. Amplitude balancing in separating P-and S-waves in 2D and 3D elastic seismic data[J]. Geophysics, 76(3): S103–S113. DOI:10.1190/1.3555529 |
[] | Virieux J. 1984. SH-wave propagation in heterogeneous media:Velocity-stress finite-difference method[J]. Geophysics, 49(11): 1933–1942. DOI:10.1190/1.1441605 |
[] | Virieux J. 1986. P-SV wave propagation in heterogeneous media:Velocity-stress finite-difference method[J]. Geophysics, 51(4): 889–901. DOI:10.1190/1.1442147 |
[] | Wang B L, Gao J H, Chen W C, et al. 2012. Efficient boundary storage strategies for seismic reverse time migration[J]. Chinese J. Geophys. (in Chinese), 55(7): 2412–2421. DOI:10.6038/j.issn.0001-5733.2012.07.025 |
[] | Wang B L, Gao J H, Chen W C, et al. 2013. Extracting efficiently angle gathers using Poynting vector during reverse time migration[J]. Chinese J. Geophys. (in Chinese), 56(1): 262–268. DOI:10.6038/cjg20130127 |
[] | Whitmore N D. 1983. Iterative depth migration by backward time propagation[C].//53rd Annual International Meeting, SEG. Expanded Abstracts, 382-385. |
[] | Yan J, Sava P. 2008. Isotropic angle-domain elastic reverse-time migration[J]. Geophysics, 73(6): S229–S239. DOI:10.1190/1.2981241 |
[] | Zhang Q S, McMechan G A. 2011. Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse time migration[J]. Geophysics, 76(5): WB135–WB149. DOI:10.1190/geo2010-0314.1 |
[] | Zhang Z, Liu Y S, Xu T, et al. 2013. A stable excitation amplitude imaging condition for reverse time migration in elastic wave equation[J]. Chinese J. Geophys. (in Chinese), 56(10): 3523–3533. DOI:10.6038/cjg20131027 |
[] | 杜启振, 秦童. 2009. 横向各向同性介质弹性波多分量叠前逆时偏移[J]. 地球物理学报, 52(3): 801–807. |
[] | 李博, 刘红伟, 刘国峰, 等. 2010. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
[] | 刘国峰, 刘洪, 李博, 等. 2009. 山地地震资料叠前时间偏移方法及其GPU实现[J]. 地球物理学报, 52(12): 3101–3108. DOI:10.3969/j.issn.0001-5733.2009.12.019 |
[] | 刘红伟, 李博, 刘洪, 等. 2010. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 53(7): 1725–1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
[] | 马德堂, 朱光明. 2003. 弹性波波场P波和S波分解的数值模拟[J]. 石油地球物理勘探, 38(5): 482–486. |
[] | 王保利, 高静怀, 陈文超, 等. 2012. 地震叠前逆时偏移的有效边界存储策略[J]. 地球物理学报, 55(7): 2412–2421. DOI:10.6038/j.issn.0001-5733.2012.07.025 |
[] | 王保利, 高静怀, 陈文超, 等. 2013. 逆时偏移中用Poynting矢量高效地提取角道集[J]. 地球物理学报, 56(1): 262–268. DOI:10.6038/cjg20130127 |
[] | 张智, 刘有山, 徐涛, 等. 2013. 弹性波逆时偏移中的稳定激发振幅成像条件[J]. 地球物理学报, 56(10): 3523–3533. DOI:10.6038/cjg20131027 |