2018, Vol. 61
2. 中国科学院地球科学研究院, 北京 100029;
3. 中国科学院大学, 北京 100049;
4. 中国科学院地质与地球物理研究所 地球与行星物理重点实验室, 北京 100029
2. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China;
4. Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
弹性波逆时偏移采用弹性波方程进行波场模拟计算，在波场的传播过程中，介质不再是声波介质，纵波入射不仅会产生透射和反射纵波，同时也会产生透射和反射横波，纵波和横波的耦合，使整个波场十分复杂(Chang and McMechan, 1986; 尧德中和周熙襄，1993).如果不将正向延拓的正演波场和逆时延拓波场中的纵横波进行分离，那么得到的偏移剖面会因纵横波相互干扰，导致偏移成像失败(Yan and Sava, 2008).目前，大家对弹性波逆时偏移中纵横波分离的必要性达成了统一认识，但是大多数学者所采用的方法仍然为传统Helmholtz分解方法(Dellinger and Etgen, 1990; Sun et al, 2004；Sun and McMechan, 2001b).Helmholtz分解假设波场由一个无旋场的梯度和一个无散场的旋度组成，这个无旋场就是纵波的势函数，无散场就是横波的势函数.因此可将波场做旋度运算得到横波波场，将波场做散度运算得到纵波波场.由于波场经过旋度和散度运算，故分离波场的物理意义和原波场不再一致，若原波场为位移，那么经过Helmholtz分离后的波场为与质点速度相关的一个量(Li et al., 2016).在求解旋度和散度运算时做了空间一阶偏导数，分离波场的振幅和相位会遭到改变.为了恢复相位和振幅，Sun等(2001a)通过Hilbert变换矫正相位的改变；Sun等(2011)用纵横波速度恢复P波和S波的相对振幅比，通过地表P波速度和频率恢复P和S波的真正振幅值.这样虽然能恢复振幅与相位，但需要存储整个时间上的波场.对逆时偏移来说，存储太大而不易实现.
马德堂和朱光明(2003)给出了纵横波解耦的二阶位移方程，采用伪谱法求解，得到了分离的纵横波位移矢量场，该方法得到的纵横波波场物理意义与总波场一致，且相位和振幅不会受到改变.李振春等(2007)在此基础上，采用P波为无旋场、S波为无散场的思想，给出了纵横波解耦的一阶速度-应力弹性波方程，将应力场和速度场都分成纵波和横波两部分.但是分离的横波应力场不仅包含横波，还有较强的纵波能量(Du et al., 2017).Zhang等(2007)同样推导了纵横波解耦的一阶速度-应力方程，该方程将速度项分为纵波速度和横波速度，二者之和为总速度，通过总速度场求取总应力场.Xiao和Leaney(2010)在一阶速度-应力弹性波交错网格方程中，引入P波应力场，通过该应力场求取P波速度矢量场，总速度矢量场减去P波速度矢量场就得到了S波的速度矢量场，并将该方法用于VSP数据的干涉法逆时偏移.Wang等(2015a)证明了Zhang等(2007)、Xiao和Leaney(2010)的公式是等价的.Yu等(2016)给出了声波-弹性波耦合的方程，从方程中能得到标量压力场.Du等(2017)也对纵横波解耦方法进行了讨论，给出了只包含横波的拟S波应力场.Gu等(2015)和Zhou等(2016)将Zhang等(2007)的公式拓展到了3D介质.本文采用纵横波解耦的一阶速度应力方程，得到了纵横波分离的波场，并与常规Helmholtz分离结果进行对比，分析了纵横波解耦方法在相位、振幅及物理意义上的优势，认为该方法适用于弹性波逆时偏移，能得到相位和振幅不受改变的偏移结果.
目前很多弹性波逆时偏移仍然在2D介质中进行.常规2D勘探假设地质构造在与测线垂直的横向上没有有变化，但是很多地质体在此方向上存在倾角，从而产生侧面反射波，二维偏移结果偏离测线方向，这种偏差给识别勘探目标带来困难.同时，对于小构造产生的绕射波，二维勘探往往不能使其归位.三维地震勘探是在空间上进行的，偏移时按照真倾角方向偏移，可以使回转波、侧面波和绕射波回归到地下真实位置，从而能识别小构造，提高对圈闭探测的成功率(French, 1974).因此我们有必要进行3D弹性波逆时偏移，得到地下准确的构造信息.纵横波解耦的波场分离方法能够保持波场分离结果的相位和振幅不受改变，因此我们采用该方法实现弹性波逆时偏移过程中的波场分离.
逆时偏移采用双程波方程，根据正演波场和逆时传播波场的互相关，就能得到地下高精度的成像结果，针对逆时偏移中的一些问题，如计算量、偏移带来的假象等均得到了很好的解决(杨仁虎等，2010；Yang et al., 2014; 王一博等，2016)，因此逆时偏移已经得到了广泛应用(Chang and McMechan, 1986; 常旭等，2006；Shabelansky et al., 2015).在弹性波逆时偏移中，需要将正演波场和逆时传播波场中的纵横波分离，然后选择正向延拓的P波和逆时延拓的P波进行成像就能得到无串扰的PP偏移结果；用正向延拓的P波和逆时延拓的S波进行成像，就能得到无串扰的PS成像结果.但即使将正演波场和逆时延拓波场中的纵横波分离，弹性波逆时偏移还需要考虑波场为矢量场的问题，因此弹性波逆时偏移的成像条件也相对声波方程更为复杂(Du et al., 2012; Wang and McMechan, 2015b).常规Helmholtz分解得到的P波为标量，2D介质中的S波为平行于y轴的矢量，可视为标量，因此不用考虑矢量波场的成像问题.但是S波在物理上为矢量，存在传播方向，经过Helmholtz变换后，S波波场在过炮点的界面法向两侧极性相反，从而导致PS成像结果在界面法向两侧存在相位反转，在叠加过程中会导致同相轴因极性不同而相互抵消，不能反映纵波转换横波的反射率信息(Duan and Sava, 2015；王维红等, 2017).在3D情况下，Helmholtz分解的P波虽为标量场，但是S波却为有三个分量的矢量场，因此对于转换波的偏移也要考虑矢量波场的成像条件问题(Du et al., 2014).
利用纵横波解耦分离的纵横波波场能很好的保持相位和振幅不受改变，但分离的纵横波都为矢量场，这给PP和PS的成像带来了新的挑战.Gu等(2015)和Zhou等(2016)直接用各个分量分别做互相关，这样每个分量都会得到一个成像结果，不便于解释.Wang和McMechan(2015b；2016)利用纵横波解耦的矢量场进行弹性波逆时偏移，在成像时，采用激发振幅成像条件(张智等，2013)，并将成像结果分为符号项和能量项，能量项根据波场振幅求取；符号项通过波场在界面法向和倾斜方向的投影来确定符号的正负，其中反射界面法向和倾斜方向通过Poynting矢量求取.Shabelansky等(2015)采用纵横波解耦的二阶方程，将检波器记录的波场逆时反传，在逆时反传过程中利用逆时反传P波矢量和S波矢量进行矢量点积，得到界面的成像结果，实现了如微震等数据在未知源的情况进行成像和速度分析.Xiao和Leaney(2010)利用矢量波场的点积实现VSP数据成像.Du等(2017)采用正演P波速度矢量场分别与逆时反传P波和S波速度矢量场的矢量点积互相关，得到了PP和PS的标量成像结果，同时利用辅助P波标量应力场的互相关得到了另一个标量PP成像结果.矢量点积成像方法实现简单，计算量小，但是因为纵横波矢量场包含波传播信息，因此矢量点积成像结果实为反射率和传播角度的综合结果，要想得到真实的反射信息，还需要根据角度做进一步矫正.
本文通过Poynting矢量来对矢量纵横波场进行标量化，对标量化后的结果采用互相关成像条件进行成像.Poynting矢量是速度和应力的张量内积，表征能流密度信息，可以用来代表波场传播方向.但是若总速度场和应力场直接进行张量内积，那么求得的Poynting矢量为纵横波综合结果，不能代表纵横波各自的传播方向，因此需要分别求取纵横波的Poynting矢量(Wang and McMechan, 2015b).对于P波Poynting矢量，P波的速度和P波应力场都能很好的求取，因此P波Poynting矢量也很方便获得.而以往给出的S波的应力场中含有纵波成分(李振春等, 2007；Du et al., 2017)，因此求取的S波Poynting矢量不准确.本文采用Du等(2017)的思想，对于S波应力场，采用纯净的S波速度场进行求取，这样得到只包含S波的拟应力场，采用该应力场求取S波Poynting矢量将更准确.
得到纵横波的Poynting矢量后，利用Poynting矢量对纵横波进行标量化.P波的传播方向和质点振动方向平行，因此通过归一化的P波Poynting矢量能很方便的实现P波场的标量化，将标量化的正演和逆时传播P波场采用常规互相关成像条件，就能得到PP成像结果；而S波的质点振动方向和波的传播方向垂直，不能很好的进行标量化，本文采用正演P波场和逆时反传S波波场的几何关系，得到逆时反传S波波场的参考振动方向，根据该方向对S波场进行标量化，通过标量化后的逆时反传S波场和正演标量P波场的互相关得到PS成像结果.
本文采用纵横波解耦方程得到了3D介质中纵横波分离的波场，并与常规Helmholtz分解结果进行了对比，发现纵横波解耦分离方法具有保相位和振幅的优势.将纵横波解耦的波场分离方法运用于逆时偏移中正演和逆时反传波场的纵横波分离，对于分离的矢量场，本文根据Poynting矢量求取了质点参考振动方向，从而进行标量化，得到了标量PP和PS成像结果.通过部分SEAM模型的偏移结果可以看出该成像条件是有效的.
1 纵横波解耦的3D矢量波场分离弹性波波场分离方法中，传统的Helmholtz分离方法基于纵波为无旋场，横波为无散场的思路，通过对波场求取散度和旋度算子得到P和S波.这种方法虽能很好的将纵横波进行分离，但是从频率-波数域的公式(1)(2)可以看出，该方法乘虚数单位i时会将相位改变90°，乘波数矢量k时会改变振幅.同时可以看出，P波为两矢量的点积，因此为标量；而S波在2D情况下(x-z平面)的叉乘为平行于y轴的矢量，可视为标量，在3D情况下为一个矢量场.若原波场为位移场，那么经过Helmholtz分解后的波场变为与速度相关的波场，物理意义和矢量特性均与原波场不一致.
|
(1) |
|
(2) |
其中，
为了得到不改变相位和振幅，且保持与原波场一致物理意义的纵横波分离方法，本文采用纵横波解耦波场分离公式，并应用于3D介质，通过模型计算观察其分离效果.
1.1 基于速度-应力方程的纵横波解耦波场分离公式弹性波方程(Aki and Richards, 2002)：
|
(3) |
其中，
|
(4) |
|
(5) |
ρ为介质密度，ui为质点位移场，σij为质点的应力，εkl为质点的应变，Cijkl为与介质有关的刚度矩阵，该矩阵有81个参数，在极端各向异性中有21个独立参数，在均匀各向同性介质中，刚度矩阵可由2个独立变量，Lame系数λ和μ表示：
|
(6) |
式中，δij为狄利克雷函数，θ为体应变.将公式(6)代入方程(3)得到：
|
(7) |
将位移u和应变ε的方程组合，得到公式(7)，其中的第一部分为纵波的传播项，该项包含的是纵波的速度，第二部分为横波传播项，包含横波速度.
在2D情况下i, j=1, 2，有如下的应力-应变关系：
|
(8) |
将2D应力-应变关系公式(8)，代入公式(7)，将方程拆为两项并将位移场转化为速度场，就得到纵横波解耦的一阶速度-应力方程：
P波：
|
(9) |
S波：
|
(10) |
总速度场为P和S波场的和：
|
(11) |
应力场与常规交错网格一致：
|
(12) |
对于3D情况，有公式(13)的应变-应力关系：
|
(13) |
将公式(7)拆为两部分，并将3D情况下的公式(13)应力-应变关系代入公式(7)得到：
P波：
|
(14) |
S波：
|
(15) |
总速度场和应力场均与2D一致：
|
(16) |
|
(17) |
在交错网格差分中，P波和S波速度场的和为总的速度场，因此在原有差分格式基础之上(Graves, 1996; Virieux, 1986)，将vPx, vSx放在vx位置，将vPy, vSy放在vy位置，将vPz, vSz放在vz位置，应力项和常规的交错网格方程一致.从公式推导过程中可以看出，该纵横波分离公式是将速度-应力方程中的质点速度项根据纵横波速度而解耦为两项，两项的和为总波场，这种方程得到的纵横波波场和原波场物理意义一致，不会改变相位和振幅信息，而且该公式中的P波项和S波项均保持了与输入波场一致的矢量特征.
1.2 3D纵横波解耦波场分离数值计算及对比分析为了验证3D波场分离的效果，采用两层模型为例，对比观察总波场、Helmholtz分离结果和本文方法分离效果.如图 1所示，给出了纵波模型，该模型大小为nx=400网格点、ny=400网格点、nz=200网格点，空间采样间隔为dx=dy=dz=10 m，上层纵波速度为vP=3000 m·s-1，横波速度为vS=2000 m·s-1，下层纵波速度为VP=4000 m·s-1，横波速度为VS=3000 m·s-1，密度为常数2.56 km·m-3.采用时间二阶，空间八阶的一阶速度-应力交错网格进行正演模拟，在地表x=2000 m, y=2000 m处加载爆炸震源.为了方便对比，在空间上，我们给出了波场快照，通过波场快照能看出本文分离效果的优势；在时间上，我们提取了地表某一位置处的正演波形记录，通过波形能更好的比较分离效果.
|
图 1 3D两层介质的P波速度模型 Fig. 1 P velocity model of 3D 2-layer media |
(1) 波场快照对比
从震源激发出P波波场，P波在界面处发生反射、透射和转换，因此在界面处存在反射P、反射S、透射P、透射S这四种波场.在正演过程中不仅能得到总波场，还能得到P波和S波分离的矢量场，对总波场求取散度和旋度运算，就能得到Helmholtz分离的P和S波.这些波场在x, y, z三个分量上都有值，由于介质在横向上没有变化，因此x方向的波场和y方向的波场振幅和相位相似，因此在这里只用x分量和z分量的波场加以对比说明.图 2中给出了在x=1500 m, y=1500 m, z=1000 m处的切片.由于采用了爆炸震源作为子波激发，因此直达波中只有P波成分，在模型中S波速度小于P波速度，所以在图 2(a，b)中，离界面最近的反射和透射波场为S波，其他波场为P波，因此从波场快照可以看出Helmholtz分离方法和本文的方法都很好的将P波和S波分离.从这几张图的对比中可以较明显的发现，Helmholtz分离的S波z分量为0，图 2i.这是因为模型反射界面水平，S波经过旋度算子后的波场矢量在xy平面上，z方向波场为0.从图 2c中可以看出，Helmholtz分离的P波为一个标量场，其相位和振幅与原始波场中x分量和z分量中的P波分量均不一致.图 2f中的Helmholtz分离的S波x分量也和原始波场x分量中的S波不一致.但从图 2(d，e，g，h)中可以看出，用本文波场分离方法的分离结果在相位和振幅上与原始波场一致，分离的P波也为矢量场，分离的结果没有改变原始波场的物理意义.
|
图 2 3D两层介质波场快照图 (a)总波场的x分量；(b)总波场的z分量；(c) Helmholtz分解得到的P波波场；(d)本方法分离得到的P波x分量；(e)本方法分离得到的P波z分量；(f) Helmholtz分解得到的S波x分量；(g)本方法分离得到的S波x分量；(h)本方法分离得到的S波z分量；(i) Helmholtz分解得到的S波z分量. Fig. 2 Particle-velocity snapshot of 3D 2-layer model (a) Snapshot of total wavefield x-component; (b) Snapshot of total wavefield z-component; (c) Decomposed P-wave by Helmholtz decomposition; (d) P-wave x-component by our method; (e) P-wave z-component by our method; (f) Decomposed S-wave x component by Helmholtz decomposition; (g) S-wave x-component by our method; (h) S-wave z-component by our method; (i) Decomposed S-wave z component by Helmholtz decomposition. |
(2) 波形对比
波场快照在空间上展示了分离效果，可以看出本文的分离方法在振幅和相位方面的优势.为了更好的观察方法的分离效果，提取了模型表面x=1500 m, y=1500 m处的地震正演波形记录，如图 3所示，在时间上观察分离效果.图 3从上到下分别展示了x, y, z三个分量的总波场的地震记录、纵横波解耦方法分离的P及S记录和Helmholtz分离的P及S地震记录.黄线为总波场地震记录，绿线和红线分别为本文方法分离得到的P波和S波记录，蓝线和紫色线分别为Helmholtz分离得到的S波和P波.Helmholtz分离的P波为标量，为了方便对比，在x, y, z三个分量上都画出了该波形.从x, y, z三个分量可以看出，由于介质在横向均匀，提取的地震波形在x=y处，因此x和y分量的地震记录完全一致.从三个分量可以看出，Helmholtz分离的P波和S波的振幅明显和原始地震记录不一致，相位也有错动.但是本文的分离方法不仅能很好的将纵横波分离，而且分离的波形能与原始波场中的对应成分完全重合，分离的P波和S波也都为矢量场.
|
图 3 两层3D介质的在x=1500 m, y=1500 m处的地震记录 Fig. 3 Seismic records at x=1500 m, y=1500 m of 3D 2 layers model |
从上述数值计算可以看出，纵横波解耦的分离方法能使分离的纵横波保持原有的矢量特征，分离波场的物理意义与原波场一致，分离结果能保证相位和振幅不受改变.
2 基于纵横波解耦的3D弹性波逆时偏移基于纵横波解耦的矢量波场分离方法能很好的实现正演波场和逆时反传波场中的纵横波分离，将正演P波和逆时反传P波进行互相关成像就能得到PP的成像结果，将正演P波和逆时反传S波进行成像就能得到PS成像结果.但是该方法分离得到的纵横波波场均为矢量场，矢量场的成像给弹性波逆时偏移带来了新的挑战.Gu等(2015)、Zhou等(2016)直接利用各个对应分量进行互相关成像，这样对于每种成像均有三个成像结果(PxPx、PyPy、PzPz和PxSx、PySy、PzSz).该成像条件得到的成像结果太多不便于解释，在倾斜介质中，也不能保证成像结果极性一致，会存在相位反转现象，影响最终的叠加效果.为了得到一个标量的成像结果，矢量点积成像条件时将正演P波矢量场和逆时反传P波矢量场进行点积成像得到PP成像结果，正演P波矢量场和逆时反传S波矢量场进行点积得到PS成像结果(Xiao and Leaney 2010; Shabelansky et al., 2015; Du et al., 2017).
|
(18) |
|
(19) |
其中，IPP和IPS分别为PP成像结果和PS成像结果，VP, forward和VP, inverse分别为正演和逆时反传P波矢量场.这样的成像结果能很好的解决矢量场的标量成像问题，方便解释.但是该方法中正演波场和逆时反传波矢量场均带有波的传播方向信息，因此其成像结果是反射率和入射角与反射角的综合效应.
|
(20) |
|
(21) |
其中，I′PP和I′PS代表真实的反射率，α和β分别为P波入射角和S波反射角.因此Du等(2017)提出在分析角道集的时候需要根据角度对成像结果进行校正.但是在普通叠加成像，不分析角道集的情况下，这种叠加结果反映的不是真实的反射率信息，不能保持波的动力学特征.而且对于PP波的成像结果中，入射角为45°时，成像结果为0，入射角在45°附近时成像能量也非常小，因此对这附近的成像结果的真实振幅信息也难以恢复.
为了更好的解决矢量波场分离的成像条件问题，本文引入了Poynting矢量.利用Poynting矢量代表波场的传播方向，根据波场传播的几何特征对正演波场和逆时反传波场进行标量化，对标量化之后的波场再进行点积成像，就能得到标量的成像结果.
2.1 Poynting矢量总波场的Poynting矢量为速度和应力的张量内积，本文在波场正演和逆时反传中采用一阶速度-应力交错网格方程，因此Poynting矢量方便求取：
|
(22) |
ni为总波场的Poynting矢量.但总波场的Poynting矢量中纵横波耦合在一起，导致不能代表真正的纵横波传播方向(Wang和McMechan, 2015b)，需要将纵横波的Poynting矢量分离.从纵横波解耦的一阶速度-应力方程能得到纵横波分离的速度场，因此若能得到纵横波分离的应力场，那么就能得到各自的Poynting矢量.李振春等(2007)给出了如下的纵横波应力公式：
P波应力：
|
(23) |
S波应力：
|
(24) |
利用上述应力场，能得到如下纵横波Poynting矢量公式(Wang and McMechan, 2015b)：
P波的Poynting矢量：
|
(25) |
S波的Poynting矢量：
|
(26) |
nP, i和nS, i分别代表P波和S波Poynting矢量.我们仍然采用图 1所示的速度模型进行正演，给出了图 4所示，600 ms处的Poynting矢量和应力场的波场快照.如图 4a所示的P波Poynting矢量的z分量，可以看出该矢量仅包含P波成分，求得的矢量能很好的代表P波传播方向.而图 4(b, c)所示的S波Poynting矢量x和z分量可以看出，求得的S波Poynting矢量不准确，尤其是在界面附近不能很好的反应传播方向.分析原因为S波应力场包含较强能量的P波成分，如图 4(b，c)的S波应力场σS, xz和σS, zz所示，该应力场不仅有S波成分，还有P波成分，因此导致S波的Poynting矢量求取不准确.Du等(2017)在方程(24)的基础上，给出了拟S波应力场.该方法将方程(24)中的速度项，由原来的总速度场替换为S波速度场，得到拟S波应力场：
|
图 4 3D两层介质Poynting矢量和应力场快照图 (a) P波Poynting矢量z分量；(b) S波Poynting矢量x分量；(c) S波Poynting矢量z分量；(d) S波σS, xz应力场；(e) S波σS, zz应力场；(f)拟S波σqS, xz应力场；(g)拟S波σqS, zz应力场；(h)用拟S波应力场得到的Poynting矢量x分量；(i)用拟S波应力场得到的Poynting矢量z分量. Fig. 4 Poynting vector and stress wavefield snapshot of 3D 2-layer model (a) P-wave Poynting vector of x-component; (b) S-wave Poynting vector of x-component; (c) S-wave Poynting vector of z-component; (d) S-wave σS, xz stress; (e) S-wave σS, zz stress; (f) S-wave σqS, xz quasi-stress; (g) S-wave σqS, zz quasi-stress; (h) S-wave Poynting vector of x-component calculated by quasi-stress; (i) S-wave Poynting vector of z-component calculated by quasi-stress. |
|
(27) |
根据此公式得到的拟S波应力场如图 4(f, g)可以看出，σqS, xz和σqS, zz这两个拟S波应力场都只包含S波成分，σqS, xz在界面附近有一点假象，这可以通过平滑速度模型减弱，也可以通过中值滤波消除(Du et al., 2017).但在本文中，拟S波应力场并没有直接使用，只用于求解S波Poynting矢量，而求解的S波Poynting矢量还需要做平滑，因此可以不进行处理.用纯净的拟S波应力场和S波速度场可以得到准确的S波Poynting矢量：
|
(28) |
nqS, i为利用拟S波应力场求得的Poynting矢量.如图 4(h, i)所示为用拟S波应力场求得的Poynting矢量x分量和z分量，用拟S波应力场求得的Poynting矢量能更准确的描述S波的传播方向.
2.2 基于Poynting矢量的矢量波场标量成像条件 2.2.1 PP成像条件对于P波而言，质点的振动方向和波的传播方向平行，可以用Poynting矢量来标量化P波场，但在这之前需要对正演和逆时反传的P波Poynting矢量进行归一化：
|
(29) |
|
(30) |
|
(31) |
该公式把正演波场和逆时反传波场投影到各自的传播方向上，这样就能得到标量的正演波场和逆时反传波场.然后对标量波场采用类似声波的震源归一化成像条件进行成像.
2.2.2 PS成像条件P波场可直接利用其Poynting矢量和质点振动矢量的点乘来进行标量化，但是反射S波场的质点振动方向和S波的传播方向垂直，因此不能直接利用传播方向来标量化S波场.为了得到与S波质点振动方向平行的参考方向，我们引入了入射P波的传播方向，如图 5所示，通过正演P波Poynting矢量叉乘逆时反传S波Poynting矢量，可以得到与入射波和反射波构成的界面垂直的矢量场，对波前面上每个质点都做该运算，就得到如图 5中圆圈所示的逆时针旋转波场.该场只与传播方向有关，不受质点振动方向影响，将该场再与逆时反传S波进行叉乘就能得到与质点振动方向平行的参考方向.从图 5可以看出，左右两边的参考方向相对，而这两边的S波质点振动方向也相对，所以用该场进行标量化，得到的标量场不会存在极性反转.另外，在倾斜地层下，传播方向会随着地层倾角而发生改变，得到的参考方向也会随着界面发生变化，故在倾斜介质中仍能有效的作为S波参考振动方向.
|
图 5 波场传播示意图 Fig. 5 Schematic diagram of wave propagation |
|
(32) |
并将参考方向归一化：
|
(33) |
nqS, inverse为拟S波应力场计算得到的逆时反传S波Poynting矢量，
|
(34) |
P波Poynting矢量和S波参考方向均与对应的质点振动方向平行，在标量化之前，对这些矢量均做了归一化处理，因此标量化的过程不会改变波场的振幅和相位，保持了波的动力学特征.而且这样的成像结果不会存在极性反转，不会引入传播角度的信息，更有利于叠加成像.
2.3 3D弹性波逆时偏移数值计算 2.3.1 层状模型为了验证对纵横波解耦的矢量场采用Poynting矢量进行标量成像的正确性，我们采用图 6所示的层状模型验证偏移结果.该模型大小为nx=600网格点，ny=600网格点，nz=200网格点，空间采样间隔为dx=dy=dz=10 m.从上到下，三层模型的纵横波速度和密度分别为(2600 m·s-1, 1200 m·s-1, 2.3 kg·m-3), (2700 m·s-1, 1300 m·s-1, 2.4 kg·m-3), (2800 m·s-1, 1400 m·s-1, 2.5 kg·m-3).震源位置位于模型表面(3000 m, 3000 m)处，我们给出了单炮偏移结果，如图 7所示，给出了在x=3000 m, y=3000 m, z=1330 m处的切片成像结果，其中z方向的切片位置刚好为最下面的界面位置.如前面的分析所知，矢量点积成像条件会引入传播角度的信息.对PP成像结果而言，在入射角为45°时，成像结果非常小.如图 7a的红色标识区所示.而基于Poynting矢量进行标量化的成像条件不会引入角度信息.如图 7b所示，红色区域的同相轴明显，没有受到角度影响，能很好的反应该位置的反射信息.同样，对于PS成像而言，矢量点积成像条件也会引入角度的影响，在小入射角和较大入射角时，成像结果偏小.如图 7c的y方向切片红色区域所示，该区域为小入射角，该图的成像结果明显小于图 7d的相同区域.同时，从图 7c和图 7d的z方向切片也可以看出，在小角度区域，矢量点积成像结果明细偏小，不能反应正常的反射信息.而基于Poynting矢量的成像结果没有受到传播角度影响，保持了波的动力学特征，能很好的对界面进行成像.
|
图 6 三层介质的P波速度模型 Fig. 6 P velocity of 3 layers model |
|
图 7 水平层状介质偏移结果 (a)矢量点积的PP偏移结果；(b)基于Poynting矢量的PP标量成像结果；(c)矢量点积的PS偏移结果；(d)基于Poynting矢量的PS标量成像结果. Fig. 7 Migration results of layer model (a) PP image result by vector dot product imaging condition; (b) Scalar PP image based on Poynting vector; (c) PS image result by vector dot product imaging condition; (d) Scalar PS image based on Poynting vector. |
本文还选取了3D SEAM模型的一部分来进行3D弹性波逆时偏移(Fehler and Larner, 2008).截取模型大小为nx=600网格点、ny=600网格点、nz=200网格点，空间采样间隔为dx=dy=dz=10 m，采用时间8阶，空间二阶的速度-应力差分格式进行波场正向延拓和逆时延拓.如图 8所示，我们给出了在x=3 km, y=3 km, z=1 km处的模型切片，图 8a, 8b, 8c分别为真实的P波速度、S波速度和密度模型，图 8d, 8e, 8f为平滑后用于偏移的P波速度、S波速度和密度模型.时间采样间隔为0.9 ms，主频为25 Hz.我们用图 8左侧的数据进行正演，得到模型表面三分量地震记录.然后采用平滑后的模型进行偏移.在x和y方向各有13炮均匀分布于模型表面，共为169炮，x和y方向的炮间隔均为50个网格点.每炮检波点均匀覆盖在每个网格上，单炮x和y方向观测系统为300网格点.如图 9所示，利用交错网格将震源正向延拓、检波点数据逆时延拓，利用公式(31)得到PP成像结果，如图 9a所示；利用公式(34)得到PS成像结果，如图 9b所示.为了更好的观察成像效果，将图 9a所示的红色区域对应的PP、PS成像结果放大显示(图 10a和图 10b).从图 9和局部放大图 10可以看出，该方法能很好的对模型结构进行偏移成像.
|
图 8 3D部分SEAM真实和平滑后的模型图 (a) P波真实速度模型; (b) S波真实速度模型; (c)密度真实模型; (d)平滑后的P波速度模型; (e)平滑后的S波速度模型; (f)平滑后的密度模型. Fig. 8 Partial of 3D SEAM model (a) P velocity model; (b) S velocity model; (c) Density model; (d) Smoothed P velocity model; (e) Smoothed S velocity model; (f) Smoothed density model. |
|
图 9 3D SEAM部分模型的偏移结果 (a)基于Poynting矢量的PP成像结果；(b)基于Poynting矢量的PS成像结果. Fig. 9 Migration results of 3D SEAM partial model (a) Scalar PP image based on Poynting vector; (b) Scalar PS image based on Poynting vector. |
|
图 10 SEAM部分模型偏移局部放大图 (a) PP成像结果局部方法图；(b) PS成像结果局部放大图. Fig. 10 Migration results of 3D SEAM partial model (a) Partial magnification of PP image result; (b) Partial magnification of PS image result. |
本文采用纵横波解耦的弹性波方程，利用速度-应力交错网格差分格式，得到了波场传播中纵横波分离的波场.与Helmholtz纵横波分离方法对比，发现纵横波解耦的波场分离方法在分离结果的振幅和相位上与原始波场能保持一致，分离波场的物理意义也与输入波场相同.将纵横波解耦的分离方法用于弹性波逆时偏移中的正演波场和逆时反传波场中的纵横波分离，得到分离的矢量纵横波场.针对矢量波场的成像问题，引入了Poynting矢量，通过该矢量来对纵横波波场进行标量化.但是S波的Poynting矢量的求取时，会有P波能量干扰，针对这一问题，采用了拟S波应力场，该应力场只包含S波而没有P波能量，利用该S波应力场和S波速度场就能得到准确的S波Poynting矢量，该矢量不仅可以用于此处转换波成像条件的求取，还可以用于准确的求取转换波角道集.在成像时，PP波的正演和逆时传播波场均用归一化后的P波Poynting矢量直接标量化，PS波中的逆时传播波场，利用正演P波Poynting矢量和逆时传播S波Poynting矢量来求取参考振动方向，通过该方向来标量化逆时传播的S波波场.得到正演P波的标量场和逆时传播P和S波的标量场后，采用类似声波方程的成像条件，就能得到PP和PS的成像结果，这样的成像结果不会存在相位反转，不会引入传播角度的信息，不需要根据角度对成像结果进行修正，保持了波的动力学特征.
致谢谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.
Aki K, Richards P G. 2002. Quantitative seismology. 2nd ed. Sausalito, CA:University Science Books.
|
|
Chang W F, McMechan G A.
1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition. Geophysics, 51(1): 67-84.
DOI:10.1190/1.1442041 |
|
Chang X, Liu Y K. Gui Z X.
2006. Zero-offset reverse time migration for prediction ahead of tunnel face. Chinese Journal of Geophysics (in Chinese), 49(5): 1482-1488.
|
|
Dellinger J, Etgen J.
1990. Wave-field separation in two-dimensional anisotropic media. Geophysics, 55(7): 914-919.
DOI:10.1190/1.1442906 |
|
Du Q Z, Zhu Y T, Ba J.
2012. Polarity reversal correction for elastic reverse time migration. Geophysics, 77(2): S31-S41.
DOI:10.1190/GEO2011-0348.1 |
|
Du Q Z, Gong X F, Zhang M Q, et al.
2014. 3D PS-wave imaging with elastic reverse-time migration. Geophysics, 79(5): S173-184.
DOI:10.1190/GEO2013-0253.1 |
|
Du Q Z, Guo C F, Zhao Q, et al.
2017. Vector-based elastic reverse time migration based on scalar imaging condition. Geophysics, 82(2): S111-S127.
DOI:10.1190/geo2016-0146.1 |
|
Duan Y T, Sava P.
2015. Scalar imaging condition for elastic reverse time migration. Geophysics, 80(4): S127-S136.
DOI:10.1190/GEO2014-0453.1 |
|
Fehler M, Larner K.
2008. SEG advanced modeling (SEAM):Phase Ⅰ first year update. The Leading Edge, 27(8): 1006-1007.
DOI:10.1190/1.2967551 |
|
French W S.
1974. Two-dimensional and three-dimensional migration of model-experiment reflection profiles. Geophysics, 39(3): 265-277.
DOI:10.1190/1.1440426 |
|
Graves R W.
1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite difference. Bulletin of the Seismological Society of America, 86(4): 1091-1106.
|
|
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. Journal of Applied Geophysics, 112: 62-78.
DOI:10.1016/j.jappgeo.2014.11.008 |
|
Li Z C, Zhang H, Liu Q M, Han W G.
2007. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm. Oil Geophysical Prospecting (in Chinese), 42(5): 510-515.
|
|
Li Z Y, Ma X N, Fu C, et al.
2016. Wavefield separation and polarity reversal correction in elastic reverse time migration. Journal of Applied Geophysics, 127: 56-67.
DOI:10.1016/j.jappgeo.2016.02.012 |
|
Ma D T, Zhu G M.
2003. Numerical modeling of P-wave and S-wave separation in elastic wavefield. Oil Geophysical Prospecting (in Chinese), 38(5): 482-486.
|
|
Shabelansky A H, Malcolm A E, Fehler M C, et al.
2015. Source-independent full wavefield converted-phase elastic migration velocity analysis. Geophysical Journal International, 200(2): 952-966.
DOI:10.1093/gji/ggu450 |
|
Sun R, Chow J, Chen K J.
2001a. Phase correction in separating P-and S-waves in elastic data. Geophysics, 66(5): 1515-1518.
DOI:10.1190/1.1487097 |
|
Sun R, McMechan G A.
2001b. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics, 66(5): 1519-1527.
DOI:10.1190/1.1487098 |
|
Sun R, McMechan G A, Hsiao H H, et al.
2004. Separating P-and S-waves in prestack 3D elastic seismograms using divergence and curl. Geophysics, 69(1): 286-297.
DOI:10.1190/1.1649396 |
|
Sun R, McMechan G A, Chuang H H.
2011. Amplitude balancing in separating P and S-waves in 2D and 3D elastic seismic data. Geophysics, 76(3): S103-S113.
DOI:10.1190/1.3555529 |
|
Virieux J.
1986. P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method. Geophysics, 51(4): 889-901.
DOI:10.1190/1.1442147 |
|
Wang Y B, Zheng Y K, Xue Q F, et al.
2016. Reverse time migration with Hilbert transform based full wavefield decomposition. Chinese Journal of Geophysics (in Chinese), 59(11): 4200-4211.
DOI:10.6038/cjg20161122 |
|
Wang W H, Zhang W, Shi Y, et al.
2017. Elastic reverse time migration based on wavefield separation. Chinese Journal of Geophysics (in Chinese), 60(7): 2813-2824.
DOI:10.6038/cjg20170726 |
|
Wang W L, McMechan G A, Zhang Q S.
2015a. Comparison of two algorithms for isotropic elastic P and S vector decomposition. Geophysics, 80(4): T147-T160.
DOI:10.1190/geo2014-0563.1 |
|
Wang W L, McMechan G A.
2015b. Vector-based elastic reverse time migration. Geophysics, 80(6): S245-S258.
DOI:10.1190/GEO2014-0620.1 |
|
Wang W L, McMechan G A. 2016. Vector-based image condition for 3D elastic reverse time migrations.//86th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, doi:10.1190/segam2016-13512297.1.
|
|
Xiao X, Leaney W S.
2010. Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution:Salt-flank imaging with transmitted P-to-S waves. Geophysics, 75(2): S35-S49.
DOI:10.1190/1.3309460 |
|
Yan J, Sava P.
2008. Isotropic angle-domain elastic reverse-time migration. Geophysics, 73(6): S229-S239.
DOI:10.1190/1.2981241 |
|
Yang P L, Gao J H, Wang B L.
2014. RTM using effective boundary saving:A staggered grid GPU implementation. Computers & Geosciences, 68: 64-72.
|
|
Yang R H, Chang X, Liu Y K.
2012. The influence factors analyses of imaging precision in pre-stack reverse time migration. Chinese Journal of Geophysics (in Chinese), 53(8): 1902-1913.
DOI:10.3969/j.issn.0001-5733.2010.08.016 |
|
Yao D Z, Zhou X R.
1993. Elastic wave reverse time migration of longitudinal and transverse wave separate imagery. Acta Geophysica Sinica (in Chinese), 36(5): 665-673.
|
|
Yu P F, Geng J H, Li X B, et al.
2016. Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration. Geophysics, 81(5): S333-S345.
DOI:10.1190/geo2015-0535.1 |
|
Zhang J L, Tian Z P, Wang C X. 2007. P-and S-wave separated elastic wave equation numerical modeling using 2D staggered-grid.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2104-2109.
|
|
Zhang Z, Liu Y S, Xu T, et al.
2013. A stable excitation amplitude imaging condition for reverse time migration in elastic wave equation. Chinese Journal of Geophysics (in Chinese), 56(10): 3523-3533.
DOI:10.6038/cjg20131027 |
|
Zhou X, Chang X, Wang Y, et al. 2017. Elastic up/down going wavefields separation and imaging condition analysis.//79th Annual International Conference and Exhibition. EAGE, 12-15, doi:10.3997/2214-4609.201701038.
|
|
Zhou X Y, Chang X, Wang Y B, et al. 2016. 3D P-and S-wave separation and elastic reverse time migration.//86th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2304-2308, doi:10.1190/segam2016-13854721.1.
|
|
常旭, 刘伊克, 桂志先.
2006. 反射地震零偏移距逆时偏移方法用于隧道超前预报. 地球物理学报, 49(5): 1482–1488.
|
|
李振春, 张华, 刘庆敏, 等.
2007. 弹性波交错网格高阶有限差分法波场分离数值模拟. 石油地球物理勘探, 42(5): 510–515.
|
|
马德堂, 朱光明.
2003. 弹性波波场P波和S波分解的数值模拟. 石油地球物理勘探, 38(5): 482–486.
|
|
王维红, 张伟, 石颖, 等.
2017. 基于波场分离的弹性波逆时偏移. 地球物理学报, 60(7): 2813–2824.
DOI:10.6038/cjg20170726 |
|
王一博, 郑忆康, 薛清峰, 等.
2016. 基于Hilbert变换的全波场分离逆时偏移成像. 地球物理学报, 59(11): 4200–4211.
DOI:10.6038/cjg20161122 |
|
杨仁虎, 常旭, 刘伊克.
2012. 叠前逆时偏移影响因素分析. 地球物理学报, 53(8): 1902–1913.
DOI:10.3969/j.issn.0001-5733.2010.08.016 |
|
尧德中, 周熙襄.
1993. 纵横波独立成像的弹性波逆时偏移方法. 地球物理学报, 36(5): 665–673.
|
|
张智, 刘有山, 徐涛, 等.
2013. 弹性波逆时偏移中的稳定激发振幅成像条件. 地球物理学报, 56(10): 3523–3533.
DOI:10.6038/cjg20131027 |
|