2. 中国石油大学 (华东) 地球科学与技术学院, 青岛 266580
2. School of Geosciences, China University of Petroleum, Qingdao 266580, China
随着地震勘探日趋向着深层复杂地区发展,面向复杂地区成像的地震处理成像方法研究越来越受到重视,由于地球广泛存在波动各向异性和黏滞性,地震波的各向异性一般体现在地下传播的地震波会发生横波分裂现象、其传播速度是传播方向的函数以及相互耦合的体波等,黏滞性呈现为地震波振幅随传播距离的增加而呈指数规律衰减.为了获得高精度保幅成像,需要校正各向异性和黏滞性对成像的影响.
各向异性介质的逆时偏移成像主要有两种方式实现:一是多波多分量成像,即求解各向异性介质的Christoffel方程得到P波和S波在频率域的偏振方向,从而求得空间域波场分离算子,进而分离出纯P波进行偏移成像 (Dellinger and Etgen, 1990; Yan and Sava, 2009).但这种方式波场外推过程计算量和内存消耗巨大,而且需要的参数较多.
另一方式是利用各向异性拟声波方程.Alkhalifah (1998)最先提出拟声波方程,他将P-SV波频散关系中对称轴上的SV波速度设为零,然后将它反变换到时间-空间域,得到拟声波方程,但是该方程为四阶微分方程.为了能提高计算效率、减少内存消耗,早期部分学者通过简化频散关系,将四阶微分方程分解成耦合的二阶微分方程组,从而得到计算效率更高的方案 (Zhou et al., 2006; Du et al., 2008).虽然上述声学近似方程能够较为方便的实现波场的数值模拟和逆时偏移,但存在残留横波干扰问题,同时它也破坏了原刚度矩阵的正定性 (Fletcher et al., 2009; Fowler et al., 2010; Zhang et al., 2011).针对各向异性介质逆时偏移不稳定的问题,Fletcher等 (2009)通过引入非零的SV波方式来克服.虽然这种做法可以保持波场正反向传播的稳定,但在偏移过程中易产生一些假象,如非零SV波与P波的串扰假象.Zhang等 (2012)则推导出稳定的TTI声波方程的自共轭微分算子,并通过能量守恒的方式证明了其方程在波场传播过程中是稳定的.Duveneck和Bakker (2011)从弹性波基本理论出发推导出更具有物理意义的各向异性声波方程,并且在TTI介质中考虑了角度 (Thomson参数) 的非均匀性,虽然能在一定程度上克服由角度变化引起的不稳定性,但还是不能完全解决残余横波的问题.程玖兵等 (2012)通过对平面波形式的弹性波方程 (即Christoffel方程) 实施一种代表向波矢量方向投影的相似变换,推导出了一种适应任意各向异性介质、运动学上与原始弹性波方程完全等价,在动力学上突出qP波的新方程,即qP波伪纯模式波动方程.
另一方面,地下介质普遍存在黏滞性,表现为对地震波的吸收衰减作用, 即地震波在介质中传播的同时,其高频成份容易丢失, 而振幅则按指数规律衰减.目前关于黏弹介质的理论模型主要有:Kelvin-Voigt模型 (Madja et al., 1985)、用连续频谱的弛豫时间表示的标准线性固体模型 (Madja et al., 1985)、标准线性固体模型 (Aki and Richards, 1980) 等.早期的标准线性固体模型是从流变学中引进而来 (Emmerich and Korn, 1987).Carcione等 (1988)在此基础上提出了更为完善的基于广义Zener体的线性黏弹理论,并结合流变学理论,提出新的黏弹性介质中的本构方程,即将应力表示为松弛函数与应变的卷积,并将其推广到各向异性黏弹介质中.Robertsson等 (1994)将交错网格有限差分算法应用到基于广义标准线性体模型的一阶速度-应力方程正演模拟中,并得到了不错的效果.
在成像方面,黏介质的偏移成像也随着黏介质理论的完善不断发展.由于它对地震波的衰减与吸收的补偿效应,使其成为实现“保幅偏移”的重要方式.早期业界广泛采用反Q滤波技术进行偏移成像,而这种技术只能进行一维的补偿,并不能准确地恢复波场的本来特征.进而由一维发展到二维,利用伪谱法实现黏声波方程的叠后深度偏移 (Dai and West, 1994).另一种实际生产中常用的方法就是结合射线理论,实现了克希霍夫反Q叠前偏移 (Traynin et al., 2008; Xie et al., 2009).然而上述方法只能对地震波振幅有所补偿,并没有考虑黏滞地层对相位的影响.Fletcher等 (2012)在波场延拓过程中同时对振幅和相位进行校正,实现了各向同性介质黏声波逆时偏移.基于带伪微分算子的常Q频散关系,Zhang等 (2010)导出时间-空间域的各向同性黏声波动方程,该方程能很好的解决相位延迟问题,并将其应用于叠前逆时偏移;Suh等 (2012)将Zhang等 (2010)的理论扩展到各向异性介质中,弥补了黏拟声波的衰减效应.
然而对各向异性黏介质的保幅成像,我们往往只侧重于某一方面的研究,即各向异性或黏滞性,这使得成像效果都不理想.所以要想实现保幅成像必须同时考虑各向异性和黏滞性的影响,而这需要一个能准确描绘地震波在各向异性衰减介质中传播规律的方程.文中通过引入标准线性固体 (SLS) 模型来模拟地震波的吸收衰减规律,并结合各向异性拟声波近似理论,推导出基于SLS模型的各向异性衰减拟声波方程.理论模型的正演模拟和逆时偏移结果表明该方程能够同时兼顾地层的各向异性和黏滞性.
1 基于拟微分算子的各向异性黏声波方程建立由徐等 (2016)推导的各向异性速度-应力拟声波方程组,经简化可以导出耦合二阶的拟声波方程,该方程与Duveneck和Bakker (2011)相一致.公式为
(1) |
其中σH与σV表示垂向与横向波场,vp、vs、ε、δ和γ为Thomson参数,将标准线性固体模型 (SLS) 扩展到各向异性介质中,得到各向异性介质中的黏声波方程为
(2) |
其中:
(3) |
其中H(t) 为Heaviside函数,rH和rV分别为垂向与横向波场的记忆变量,τσ和τε分别代表应力和应变松弛时间,均可由品质因子Q表示,公式为
(4) |
(5) |
其中ω为角频率,且τ=τε/τσ-1,由于实际地层的Q值一般都大于20(Robertsson et al., 1994),所以τ « 1.若将方程2经多位傅里叶变化,可得到频率波数域表达式为
(6) |
其中ωH和ωV为角频率,kx、ky和kz为空间波数.将H(t) 函数分解,可将方程6右边中的第二项分为两项:
(7) |
简化方程7可得到:
(8) |
方程8所示的频散关系与Carcione等 (1988)提出的标准线性固体系列 (SLSs) 模型相一致.对方程8因式分解 (分母化),可得到:
(9) |
有上述方程8和9, 及τ « 1,很容易理解下面两项假设:
(10) |
将方程10带人方程9中,化简可得如下频散关系为
(11) |
根据复数的基本性质,我们可以将方程11转换为
(12) |
因为τ « 1,可将方程11和12化简为
(13) |
(14) |
联立上述两式得到最终的频散关系式为
(15) |
方程15左边第二项与第三项为波场的传播项,分别描述了波沿着水平方向与垂直方向的传播规律.最后一项为衰减项,描述了波传播过程中的衰减规律.可以看出角频率ωH的衰减项除了与水平方向上的波数kx、ky有关,同时也受到垂直方向上波数kz的影响,角频率ωV也是如此,这主要是由于介质的各向异性对波场衰减项的影响.实际处理中,往往人们忽略这种影响, 这样可将方程写成方程16式为
(16) |
将上述方程转换至时间空间域,公式为
(17) |
可以看出方程17两式中的衰减项只针对于各自的波场进行吸收衰减.
2 模型试算由于方程中的拟微分算子是对偏导数的开方,这很难用传统的有限差分法实现数值模拟,故本文采用伪谱法来实现.具体的实现方式为:对时间域导数仍应用有限差分法,而对空间域偏导数使用伪谱法实现.并且本文笔者将伪谱法扩展到能求解空间域分数阶偏导数.设u(x) 为求解的函数,根据傅式微分的性质,可以一阶、二阶的偏导数表示为
(18) |
(19) |
将傅式微分扩展到分数阶,空间域分数阶偏导数则可表示为
(20) |
其中FFT、iFFT分别为傅里叶正、反变换,β为分数.利用上式的方法对方程17式进行求解,其差分格式为
(21) |
24式就构成求解黏声VTI介质伪谱法差分格式,其中am为差分系数,L为差分阶数.
3 二维均匀各向异性介质模型正演模拟为了验证新推导的各向异性衰减拟声波的衰减特性及拟声波的运动学特征,选择简单的二维均匀各向异性介质模型进行正演模拟.模型概述:各向异性参数Vp0=2000 m/s,ε=0.24, δ=0.1;模型大小为600×600,网格大小dx=dz=5 m;震源为主频40 Hz的雷克子波,时间采样间隔为0.5 ms;Q模型:分别为20、50和无穷大 (完全弹性介质, 本文用Q=10000来替代).如图 1为三个不同Q值的波场快照.图 3为二维XOZ平面内绕Z轴旋转45度的TTI介质的正演结果.图 2为模型中坐标为 (1 km,1 km) 点处的波形图和振幅谱.
从图 1与图 3的波场快照可见, 一方面拟声波波前面呈现椭圆状,这符合各向异性介质下准P波的运动学特征,即准P波在各个方向上的相速度是不一致的;另一方面随着Q值的降低,衰减作用越明显,在波场快照中表现为同相轴变模糊,这也符合标准线性固体衰减模型的一般规律.结合图 2中的波形图和振幅谱可知随着Q值的不断减小,吸收衰减越明显,振幅值不断减小,主频不断降低,并向低频靠近,同时频带宽度变窄,高频成分衰减更明显.说明本文的各向异性衰减拟声波方程能较为准确地描述黏声各向异性介质中地震波的运动学特征.
其次,为了验证上文引入的规则化算子在各向异性逆时偏移成像中的作用,笔者应用简单的层状模型,如图 3所示.模型横、纵向均为1.5 km,采样间隔为10 m;其各向异性参数ε和δ分别为0.2、0.1;品质因子为50;上层速度为2 km/s,下层速度为2.5 km/s.为了压制qSV波的干扰,将震源点和接受点位置设置为各向同性介质.
图 4为不同时刻的波场快照,(a) 没有加规则化算子,(b) 添加规则化算子.从图中看出未加规则化算子的波场逆时传播,会产生不稳定现象,而且随着时间的逆推,会越来越不稳定;加入规则化算子后,不稳定现象明显消除.对比图 5各成像结果 (VTI黏声波数据进行VTI声波偏移、VTI黏声波数据进行VTI黏声波偏移以及VTI声波数据进行VTI声波偏移),TI声波逆时偏移没有考虑吸收衰减效应,成像结果保幅性较差;VTI黏声波逆时偏移补偿了反射界面到检波器传播过程中产生的振幅衰减,基本上恢复了由于衰减引起的振幅的降低,验证了通过引入规则化算子,消除逆时传播不稳定,在成像过程中考虑黏滞性,提高了成像的分辨率,振幅的保真度.
为验证文中推导的各向异性衰减标量波方程在逆时偏移中的可行性,文中选VTI_HESS模型进行试算, 如图 6中的 (a)-(d) 分别为P波速度、Q模型、模型以及模型.其中Q值取值为30到500,深层有一个强的衰减体,与其相对应的速度谱上有一个高速体,P波速度取值为1500 m/s到4500 m/s, 和的取值范围分别为0到0.28以及0到0.16,其中和的模型中有个各向异性的突变体 (图中已经标出);采用25 Hz的雷克子波,总共600炮,炮间距为10 m,每炮检波器个数为600,检波器间隔5 m,时间采样长度为4.0 s,时间采样间隔为4 ms,成像条件为震源归一化互相关.
如图 7中的 (a)-(d) 分别是常规声波方程逆时偏移、各向同性黏声波方程逆时偏移、各向异性拟声波方程逆时偏移以及各向异性衰减拟声波方程逆时偏移成像结果.通过对比看出,常规声波方程逆时偏移构造形态成像混乱,构造复杂区域绕射波没有完全收敛;特别是深部的各向异性突变体几乎不能显现;深部能量弱,成像效果差.各向同性黏声波方程逆时偏移只是对构造深部的振幅有所补偿,但没有解决常规声波方程逆时偏移的问题.这主要是因为各向同性介质下的声波方程不能完整的描述各向异性介质下标量波的传播规律,比如P波在各个方向上的传播速度不一致.所以该方法不能校正介质各向异性引起的绕射波不归位、能量不聚焦等现象,尤其是大倾角地区.各向同性拟声波方程逆时偏移对HESS模型断层构造成像清晰,在高陡倾角断层面处也有很好的成像效果,同相轴连续性较好,地下突变体轮廓明显.但由于没有考虑黏滞地层的吸收衰减,所以模型的深层成像效果不是很理想,主要是能量较弱,同相轴不清晰.而各向异性衰减拟声波方程逆时偏移即能很好的描述各向异性地震波传播规律,又考虑了地层的吸收衰减影响.从成像效果 (d) 看出:同相轴清晰、深部能量强,振幅分布更加均衡.
本文将标准线性固体模型引入到各向异性介质中,以各向异性二阶拟声波方程组为出发点,推导出时间空间域的各向异性介质黏声波方程组.通过理论分析和模型试算,得到了如下几点结论与认识:
(1) 方程中用一个伪微分算子来描述地震波的衰减特征,模型的正演模拟表明该算子能很好的描述地震波在黏介质中的衰减规律,能很好的结合介质的各向异性;另一方面,该方程也能准确的描述各向异性介质下的标量波的传播规律.而且该方程中没有记录变量,也就消除了实际地震勘探中没有该记录导致的不稳定问题.
(2) 将本文推导的各向异性介质黏声波方程运用到逆时偏移中.成像结果表明该方程能同时校正介质的各向异性和黏滞性对地震波的影响,成像剖面的分辨率和振幅保真度能得到较大的提高.
(3) 为了解决波场逆传时的高频不稳定现象,文中通过引入规则化算子来构建稳定的逆时传播方程,该算子能很好的压制额外的高频成分,保证波场逆冲中的稳定性.
5.2本文的各向异性黏声波方程是基于声学近似推导而来,虽然能较为准确的表述准P波的运动学信息,但同时存在残留横波,动力学特征不准确等问题.另一方面,本文的黏声波逆时偏移只补偿了地震波在黏介质中的吸收衰减,并没有补偿地震波在透射与反射过程中的能量损失;另外还存在低频噪声、相位难以校正、波场逆时传播高频不稳定、深部幅值弱等问题.故如何实现各向异性黏介质的保幅成像,仍有待于更进一步的研究.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持![] | Aki K, Richards P G. 1980. Quantitative Seismology:Theory and Methods[M]. San Francisco: W.H. Freeman and Company. |
[] | Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media[J]. Geophysics, 63(2): 623–631. DOI:10.1190/1.1444361 |
[] | Carcione J M, Kosloff D, Kosloff R. 1988. Viscoacoustic wave propagation simulation in the earth[J]. Geophysics, 53(6): 769–777. DOI:10.1190/1.1442512 |
[] | Carcione J M. 1990. Wave propagation in anisotropic linear viscoelastic media:Theory and simulated wavefields[J]. Geophysical Journal International, 101(3): 739–750. DOI:10.1111/gji.1990.101.issue-3 |
[] | Cheng J B, Chen M G, Wang T F, et al. 2014. Description of qP-wave propagation in anisotropic media, Part Ⅱ:Separation of pure-mode scalar waves[J]. Chinese Journal of Geophysics (in Chinese), 57(10): 3389–3401. DOI:10.6038/cjg20141025 |
[] | Cheng J B, Kang W. 2012. Propagating pure wave modes in general anisotropic media, Part Ⅰ:P-wave propagator[C].//82nd Ann. Internat Mtg., Soc. Exp. Geophys. Expanded Abstracts, 1-6. |
[] | Cheng J B, Kang W. 2014. Simulating propagation of separated wave modes in general anisotropic media, Part Ⅰ:qP-wave propagators[J]. Geophysics, 79(1): C1–C18. DOI:10.1190/geo2012-0504.1 |
[] | Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part Ⅰ:Pseudo-pure-mode wave equations[J]. Chinese Journal of Geophysics (in Chinese), 56(10): 3474–3486. DOI:10.6038/cjg20131022 |
[] | Chu C L, Macy B K, Anno P D. 2011. An accurate and stable wave equation for pure acoustic TTI modeling[C].//81nd Ann. Internat Mtg., Soc. Exp. Geophys. Expanded Abstracts, 179-184. |
[] | Dai N X, West G F. 1994. Inverse Q migration[C].//64th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 1418-1421. |
[] | Dellinger J, Etgen J. 1990. Wave-field separation in two-dimensional anisotropic media[J]. Geophysics, 55(7): 914–919. DOI:10.1190/1.1442906 |
[] | Du X, Fletcher R P, Fowler P J. 2008. A new pseudo-acoustic wave equation for VTI media[C].//70th EAGE Conference and Exhibition incorporating SPE EUROPEC 2008. Extended Abstracts, H033. |
[] | Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media[J]. Geophysics, 76(2): S65–S75. DOI:10.1190/1.3533964 |
[] | Emmerich H, Korn M. 1987. Incorporation of attenuation into time-domain computations of seismic wave fields[J]. Geophysics, 52(9): 1252–1264. DOI:10.1190/1.1442386 |
[] | Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media[J]. Geophysics, 74(6): WCA179–WCA187. DOI:10.1190/1.3269902 |
[] | Fletcher R P, Nichols D, Cavalca M. 2012. Wavepath-consistent effective Q estimation for Q-compensated reverse-time migration[C].//74th EAGE Conference and Exhibition Incorporating EUROPEC 2012. Extended Abstracts. |
[] | Fowler P J, Du X, Fletcher R P. 2010. Coupled equations for reverse time migration in transversely isotropic media[J]. Geophysics, 75(1): S11–S22. DOI:10.1190/1.3294572 |
[] | Liu F Q, Morton S A, Jiang S S, Net al. 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media[C].//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2844-2848. |
[] | Madja G, Chin R C Y, Followill F E. 1985. A perturbation theory for Love waves in an elastic media[J]. Geophysical Journal International, 80(1): 1–34. DOI:10.1111/gji.1985.80.issue-1 |
[] | Pestana R C, Ursin B, Stoffa P L. 2011. Separate P-and SV-wave equations for VTI media[C].//81st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 163-167. |
[] | Pestana R C, Ursin B, Stoffa P L, et al. 2011. Separate P-and SV-wave equations for VTI media[C].//81st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 163-167. |
[] | Robertsson J, Blanch J O, William W S. 1994. Viscoelastic finite-difference modeling[J]. Geophysics, 59(9): 1444–1456. DOI:10.1190/1.1443701 |
[] | Suh S, Yoon K J, Cai J, et al. 2012. Compensating visco-acoustic effects in anisotropic resverse-time migration[C].//82nd Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 1297-1301. |
[] | Thomsen L. 1986. Weak elastic anisotropy[J]. Geophysics, 51(10): 1954–1966. DOI:10.1190/1.1442051 |
[] | Traynin P J, Liu J, Reilly J M. 2008. Amplitude and bandwidth recovery beneath gas zones using Kirchhoff prestack depth Q-migration[C].//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2412-2416. |
[] | Xie Y, Xin K F, Sun J, et al. 2009. 3D prestack depth migration with compensation for frequency dependent absorption and dispersion[C].//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2919-2923. |
[] | Yan J, Sava P. 2009. Elastic wave-mode separation for VTI media[J]. Geophysics, 74(5): WB19–WB32. DOI:10.1190/1.3184014 |
[] | Zhan G, Pestana R C, Stoffa P L. 2012. Decoupled equations for reverse time migration in tilted transversely isotropic media[J]. Geophysics, 77(2): T37–T45. DOI:10.1190/geo2011-0175.1 |
[] | Zhang Y, Wu G C. 2013. Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration[J]. Chinese Journal of Geophysics (in Chinese), 56(6): 2065–2076. DOI:10.6038/cjg20130627 |
[] | Zhang Y, Zhang H Z, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation[J]. Geophysics, 76(3): WA3–WA11. DOI:10.1190/1.3554411 |
[] | Zhang Y, Zhang P, Zhang H Z. 2010. Compensating for visco-acoustic effects in reverse-time migration[C].//80th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3160-3164. |
[] | Zhou H B, Zhang G Q, Bloor R. 2006. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[C].//76th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 194-198. |
[] | 程玖兵, 陈茂根, 王腾飞, 等. 2014. 各向异性介质qP波传播描述Ⅱ:分离纯模式标量波[J]. 地球物理学报, 57(10): 3389–3401. DOI:10.6038/cjg20141025 |
[] | 程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述Ⅰ:伪纯模式波动方程[J]. 地球物理学报, 56(10): 3474–3486. DOI:10.6038/cjg20131022 |
[] | 张岩, 吴国忱. 2013. TTI介质qP波逆时偏移中伪横波噪声压制方法[J]. 地球物理学报, 56(6): 2065–2076. DOI:10.6038/cjg20130627 |