2016, Vol. 59
2. 大地测量与地球动力学国家重点实验室, 武汉 430077
2. State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China
多波多分量地震数据采集技术能够获得包含多波信息的矢量场地震数据，不同波型反映了探测目标的不同信息. 多波的传播属性与地下介质物性参数有直接关系，可以利用纵波和横波传播的运动学和动力学属性及其差异，来研究油气藏储层的岩石岩性和物性、岩石裂隙发育状况、流体属性检测等(Li，1997; Li et al.，1999; Qian et al.，2007). 因此，多波数据能够提供仅依靠纵波所不能获得的油藏参数评估(Davis，2001). 然而，当前更大的挑战是多分量数据的准确处理和多波地震信息的有效提取和利用.
常规的处理方法是基于标量波动方程的，主要有两种方法: 一种方法是认为PP波和PS波能量分别集中在垂直分量和水平分量，对单分量数据进行标量处理; 另一种方法是在偏移前进行波型分离，对分离后的数据进行声波偏移(黄中玉等，2007; 许世勇等，1999). 常规处理方法主要利用波的传播和偏振方向等几何属性实现多波型的识别，而没有考虑与弹性参数相关的动力学属性. 但是在多数情况下，多波型之间是相互耦合的，难以彻底分离. 另外，标量波动方程不能准确反映弹性波传播和扩散的动力学性质，不能准确实现弹性波的保幅偏移和反演.
自20世纪80年代，基于弹性波动方程的多分量偏移方法已经出现. 弹性波动方程偏移方法能够准确反映弹性波的波场传播. 该偏移方法主要有两个发展方向：弹性逆时偏移和弹性Kirchhoff偏移.
Sun和McMechan(1986)、Chang和McMechan(1987)利用有限差分实现了弹性逆时偏移，Chang和McMechan(1994)将其推广到3D空间. Mittet(1994)讨论了交错网格有限差分法存在的应力边界条件问题. Zhe和Greenhalgh(1997)在地表及浅层利用弹性波动方程进行波场外推，获得独立的P波和S波位移势能，以位移势能代替位移向下延拓. Sun和McMechan(2001)在浅层利用散度和旋度进行波场分离，之后进行标量偏移，同时考虑了S波极性翻转问题. Sun等(2006)将其推广到三维三分量(3D3C)偏移. Yan和Sava(2008，2009)实现了角度域和各向异性介质的弹性逆时偏移. Ravasi和Curtis(2013)同时利用海底质点速度矢量和声波水压四分量(4C)地震数据实现海底弹性逆时偏移.
Kirchhoff积分方法是最灵活、最高效的偏移方法，并且也是很成熟的保幅偏移方法.Pao和Varatharajulu(1976)提出了非均匀介质中的弹性Kirchhoff-Helmholtz 积分(简称弹性KH积分)，Kuo和Dai(1984)、Dai和Kuo(1986)将其应用于弹 性波场外推. Berkhout和Wapenaar(1989)、 Wapenaar和Haimé(1990)分别推导了声波和弹性波单程波Rayleigh积分，并将Rayleigh积分应用到 波场延拓. Hokstad(2000)基于弹性/黏弹性Kirchhoff 积分和Claerbout提出的沉降观测(survey-sinking)成像条件实现了多分量偏移，并依据单极或双极震源与应力或位移边界条件的等价性处理了边界问题. Schleicher等(2001)讨论了各向异性介质中弹性KH积分的几何射线近似(GRA). Druzhinin(2003)应用稳相性质和高频近似推导了无应力边界条件下的Kirchhoff积分偏移，并利用极化滤波实现波型分离.
当前，基于弹性波动方程的多分量数据处理方法仍处于发展阶段，多分量偏移方法存在多种问题有待讨论. 观测到的多分量地震数据往往是边界点的位移或质点速度矢量，作为弹性波动方程边界条件，这是不完整的(Mittet，1994). 求解弹性波动方程同时需要位移矢量(或质点速度矢量)边界条件和应力边界条件. 在边界条件不明确的情况下，不同波型不能准确分离，偏移成像也会形成明显假象. 为了解决应力边界条件问题，多数弹性逆时偏移方 法(Sun and McMechan，1986; Chang and McMechan，1987; Zhe and Greenhalgh，1997; Sun and McMechan，2001)和部分弹性Kirchhoff积分方法(Kuo and Dai，1984; Sena and Toksoz，1993)将应力边界条件默认假设为自由地表边界条件，即边界上正应力和剪切应力均为零. 为此，Rechtien(1985)曾对Kuo和Dai(1984)所做的设定提出质疑. 这种假设只适用于理想的陆地地表探测，对VSP和OBC地震数据均不符合. Druzhinin(2003)假设P、S波传播相互独立，没有耦合，利用极化滤波实现了弹性Kirchhoff偏移中的波型分离. 现代海底观测系统可以利用4C探测器同时记录海底质点速度矢量和声波水压，Ravasi和Curtis(2013)利用声波水压代替应力边界条件实现了海底弹性逆时偏移，这需要增加数据采集量和计算量.
常规互相关成像条件不适用于矢量波场成像，并且PS转换波存在极性翻转问题. 尤其在3D情况下，S波的椭圆偏振、旋度矢量等属性使矢量波场成像变得更为复杂. P波震源和S波震源成像之间存在相互串扰，当P、S波震源都比较强的时候，相互串扰需要压制. 因此，弹性矢量场成像需要设计适合自身的成像条件.
除此之外，当前偏移方法本身也存在多种问题. Kirchhoff积分偏移存在常规射线理论的多种问题，比如多路径问题、焦散区和阴影区、临界和超临界区等. 逆时偏移计算量和内存需求量过大，并且存在高频频散现象和低频反射干扰问题等.
高斯束偏移具有Kirchhoff偏移的优点，比如计算高效，陡倾界面成像. 高斯束偏移可以自动解决多路径和多走时问题，这是传统Kirchhoff偏移不能解决的问题. Hill(1990，2001)提出了叠后和叠前高斯束偏移. 此后，高斯束偏移被推广到各向异性偏移(Alkhalifah，1995; Zhu et al.，2006)、不同域偏移(Nowack et al.，2003; Gray，2005; 岳玉波，2011; 段鹏飞等，2013)和保幅偏移(Gray and Bleistein，2009). 然而，当前这些偏移方法都是用来处理标量数据，有关多波多分量弹性波场数据的高斯束偏移方法还没有做太多的研究.
本文详细分析了不同边界情况下多波型在弹性各向同性介质中的传播情况(自由空间、海底、自由地表)，确定了不同探测环境的应力边界条件. 通常我们所关注的有效信号是上行波，基于上行波传播假设，建立了平面波场位移与应力关系式，进而确定了模型的应力边界条件(Li and Mao，2015). 在此基础上，准确推导了3D多波多分量高斯束波场延拓和偏移表达式. 多波型波场分离在延拓和偏移过程中自动实现，不需要在偏移前进行单独处理. 为了解决矢量波场成像条件问题，本文采用了散度/旋度互相关成像条件，并且通过约定PS转换波的正向旋转方向解决了3D空间的极性翻转问题. 本文给出了2D和3D模型数据的偏移成像结果，验证了所提供方法在波场延拓、多波型自动分离和偏移成像中的准确性.
2 符号约定由于本文涉及大量公式推导，在此对多种表示形式进行统一约定. 对于矢量和张量的表示，本文引用 Červený(2001)的表示形式. 对于多分量物理量，本文同时使用了分量表示与矢量(张量)表示形式，比如 MIJ为 M 的分量， enν为 e ν 的分量等. 本文中的公式遵守爱因斯坦求和约定，并且该约定只应用于下标符号，上标符号不执行求和约定，比如：
|
为了公式表示和计算执行方便，3D空间统一使用右手正交坐标系，并且约定下标的对应关系为(1，2，3)⇔(z，x，y)，即z轴为第一维，y轴为最后一维，比如：
|
本文使用3D中心射线坐标系(q1，q2，s)(Červený，2001)，对应的坐标基为 e 1，e 2，e 3≡ t，如图 1所示.
|
图 1 3D中心射线坐标系 Fig. 1 3D ray-centered coordinate system |
由于弹性波传播方向和偏振方向在偏移过程中的重要性以及弹性KH积分的复杂性，我们需要对慢度矢量做正向约定. 在没有特殊说明的情况下，KH积分和格林函数中出现的慢度矢量 p 及其分量 均满足图 2所示的正向约定，即均由边界指向介质内侧. 这样约定是为了方便高斯束的应用，该约定与标量波场中常用的方向有所不同，如布莱斯坦等(2004).
|
图 2 正向慢度矢量 Fig. 2 Positive direction of slowness vectors |
符号α、β、υ分别代表P波波速、S波波速和多波型波速的统一表示. 除了在均匀介质中和坐标特殊标示情况下，均表示边界上的物性参数. 符号ν用于表示不同波型，并且作为上标符号显示. 在3D各向同性介质中，S波一般情况下为椭圆偏振，为了表示和计算方便，将其分解为S1和S2两种线性偏振横波. 图 3显示了S波在中心射线垂直截面中的椭圆偏振和S1、S2波的线性偏振. 因此，ν表示P、S1、S2三种波型，并且以 e ν 表示三种波型偏振方向矢量. e ν 均为单位矢量，且与中心射线坐标基有以下对应关系 e S1= e 1，e S2= e 2，e P= t . 因此，文中对矢量 e ν和e n 不做区分.
|
图 3 S波的椭圆偏振 Fig. 3 Elliptic polarization of S waves |
高斯束是波动方程在高频傍轴近似情况下的一种特解. 它假设波场频率较高，并且认为在中心射线附近波场变化缓慢，可用二阶抛物型方程表示(Červený et al.，1982). Červený和PŠenČík(1983)推导了3D弹性非均匀各向同性介质中的高斯束表达式，本文给出了它的归一化多波表示形式：
|
(1) |
其中， u GBv为点 x 的ν 波高斯束位移矢量， r和p ν 分别为高斯束的入射点和入射慢度矢量， q =(q1，q2)T 为垂直射线的平面坐标系. 动力学参数 Q 为2×2矩阵，反映了波场能量的几何扩散情况. a ν为ν 波偏振方向矢量，且有关系式
|
(2) |
其中， MIJ为2×2 矩阵 M 的元素，为波前走时在(q1，q2)坐标系的二阶偏导数. 式(2)包含了与坐标qJ相关的附加项，这是偏振方向矢量的校正项.当 x 离中心射线很近时， qJ 很小，校正项可以忽略不计.
3.2 格林函数张量弹性各向同性介质中的格林函数张量具有高频近似表示形式(Červený，2001)：
|
(3) |
式中 Uiv为ν 波点源波场位移矢量. 格林函数张量 Gin(x，r，ω)表示点 r 沿第n 维方向激起的单位体力震源在点 x 波场位移的第 i 分量，它满足互易关系
|
在均匀介质中 U ν 可以表示为
|
(4) |
它可表示成弹性高斯束的叠加形式(e.g.，Hill，2001)：
|
(5) |
式中， p ν=(pzv，pxv，pyv) 为慢度矢量，其中， px和py 是与波型无关的两个独立变量，为了方便起见，记为 p =(px，py). p ν 的三个分量满足关系式 pzv= (1/υ2－px2－py2)1/2，其中υ 与波型有关. 利用式(3)和(5)可将格林函数张量 Gin(x，r，ω) 表示为多波型高斯束的叠加.
3.3 弹性KH积分弹性KH积分相当于弹性波动方程的积分表示形式(Červený，2001; Sena and Toksoz，1993). 积分方程反向延拓可以表示为
|
(6) |
这里符号*代表复共轭， nj 为边界 Σ 的外法线方向矢量的分量， dSr 为边界点 r 在Σ 的微元面， τij 为应力张量，并且满足广义虎克定律
|
(7) |
|
(8) |
式中， cijkl 为弹性参数张量，在各向同性介质中可以由两个独立的Lamé弹性模量 λ、μ 表示，
|
(9) |
其中 δij 为Kronecker Delta函数. Lamé模量与波速和密度有关系式 λ=ρ(α2－2β2)，μ=ρβ2.
从积分方程(6)可以看出，在物性参数模型确定的情况下，格林函数 Gni和Hnji 均可通过计算求得. 而积分方程波场延拓同时需要应力边界条件 τij(r，ω) 和位移边界条件 ui(r，ω) 作为输入. 当只有 ui(r，ω) 已知时，不能获得波场的唯一解. 只有在应力边界条件同时确定时才能准确求解积分方程. 应力张量 τij(r，ω) 虽然共有9个分量(其中6个独立分量)，但由求和项 τijnj可知，τij 中只有3个独立分量在式(6)中起作用，即边界 Σ 上的正应力和剪切应力(Červený，2001).
为了公式推导方便，将应力位移关系式(7)和(8)代入积分方程(6)，反向延拓积分方程可以写成位移及其偏导数表示形式：
|
(10) |
文中设定 Σ 外法线方向为 n =(－1，0，0)，微元dSr=dxrdyr， 这时式(10)表示上行波的反向延拓，可以展开表示为
|
(11) |
利用射线理论高频近似关系式 
|
(12) |
当同时存在应力和位移边界条件时，波场上下行传播情况可以明确识别. 由于我们只关注上行波的偏移成像，因此在式(12)推导过程中只考虑上行波的传播. 当边界条件不完整，缺少应力边界条件时，会出现上下行传播方向混乱(Mittet，1994). 因此，建立有效的应力边界条件求解方法非常重要.
3.4 弹性高斯束波场延拓首先，通过局部倾斜叠加实现局部平面波分解. Hill(2001)给出了高斯窗口划分归一化表示形式：
|
(13) |
式中， L =(0，Lx，Ly) 为束中心坐标， ωl和wl分别为参考频率和参考有效半宽度，常数 a 为单个束中心网格面积大小.
在局部范围内，假设 λ(r)≈λ(L)，μ(r)≈μ(L)，p ν(r)≈ p ν(L). 点源波场高斯束叠加的相移表示形式为
|
(14) |
将式(13)和(14)代入式(12)，并将多波型分开表示，整理可得
|
(15) |
式中 uGB，nv表示 u GBv的第n 分量.(15)式比声波高斯束方法要复杂得多，被积函数中包含一系列多项式. 为了进一步整理式(15)，对 dSr和dpxdpy 调换积分顺序. 由于在局部范围内，积分多项式中除了位移分量 un 及其空间偏导数之外，其他量均近似为 常数项. 因此，只对多项式中的 un和∂un/∂xlr进行dSr 的积分，即局部倾斜叠加(local slant stack)(Hill，2001)，其中 ∂un/∂xlr 与应力边界条件有关.如果设定函数 f(r，ω) 的局部倾斜叠加为LSS，
|
(16) |
LSS是一种加窗积分变换，它将 (r =(rx，ry)，ω) 域的数据变换到 (p =(px，py)，ω) 域. 这时， un 及其空间偏导数的局部倾斜叠加有关系式
|
(17a) |
|
(17b) |
Dn 为分解后的局部平面波. 位移分量 un 包含了多种波型，因此局部平面波 Dn 也是多种波型的叠加. 由广义虎克定律可知，应力张量可以通过求取位移矢量的空间偏导数获得，并且水平边界面上的正应力和剪切应力均与位移在 z 方向的偏导数有关. 由式(17)可以看出， un在x、y 方向的偏导数很容易与 un 建立关系，在 z 方向的偏导数比较难以确定. 位 移在 z 方向的偏导数同时与波的上下行传播方向、模型边界状况和波型类别有关. 以下将对三种最常见波场传播模型进行具体应力分析和位移空间偏导数求取. 以简单层状模型为例，如图 4所示，当观测边界 Σ(z=0) 以上分别为固体介质、流体介质和真空时，模型空间分别为自由空间、海底、自由地表.
|
图 4 简单层状模型 Fig. 4 Simple layering model |
如图 4所示，当观测边界 Σ 以上为固体介质时，模型称为自由空间模型. 假设该模型边界 Σ 上下波速和密度是连续且局部均匀的，如图 5所示. 实际应用中多种情况符合该假设，比如边界为完全匹配吸收的有限差分波场合成，VSP数据采集以及其他地下信号接收等. 由于边界附近的介质连续且均匀， ν 波入射到边界 Σ 后不会激起反射波、转换波和透射波，因此边界上多波型的传播是相互独立的，没有耦合. 这时，多波型的识别与分离只与波的传播和偏振方向等几何属性有关.
|
图 5 自由空间上行波 Fig. 5 Upgoing waves in free-space model |
下面对多波型局部平面波 Dn(L，p，ω) 进行分析. Dn 是相互独立的 ν 波局部平面波 Dnv 的叠加，即 
|
(18) |
puzv 为上行波慢度 z 分量. 由此可知， puzv 与波型 ν 有关. 由于上行波方向与积分方程(15)慢度矢量方向相反，所以有 puzv=－pzv，其中，pzv 为积分方程(15)慢度z分量.
如果将矢量场 Dnv 写成标量场和偏振方向矢量表示形式
|
(19) |
多波型局部平面波可以表示为
|
(20) |
式中， env为3×3 偏振方向矢量矩阵，是实现自由空间模型多波型分离的关键因子. env 矩阵是通过矢量 p ν 建立的. 对于某一特定的 p =(px，py)，由关系式pzv=(1/υ2－px2－py2)1/2 可得3种波型的慢度矢量 p ν， 由此可确定3条中心射线在边界 Σ 上的坐标基矢量 e 1，e 2，t (注意：3条中心射线中共有9个坐标基矢量，但是其中只有3个矢量与矩阵 env 有关). 根据单位矢量 e n与 e ν 对应关系可知， p P射线的 t 、 p S1射线的 e 1、 p S2射线的 e 2， 分别对应矩阵 env中的 e P、 e S1、 e S2.矩阵 e nv 建立以后，标量场 Dν 的求解可通过矩阵 env的求逆实现. 因此，利用式(20)可以将多波型局部平面波 Dn 分解为相互独立的 Dν.
在上行波假设情况下，分离后的 ν 波局部平面波及其空间偏导数有以下关系式：
|
(21a) |
|
(21b) |
这样，位移在z方向的偏导数也同时获得了. 将式(21)代入式(15)，整理可得
|
(22) |
由此可见，在自由空间模型中，弹性KH积分方程中前后两部分被积项是相等的. 被积函数的多项式可以作进一步整理.
P波：
|
S1波：
|
S2波与S1波类似. 这样，式(22)可以进一步简化为
|
(23) |
式(23)就是所获得的自由空间传播模型的波场延拓积分方程. 该表达式与声波高斯束波场延拓表达式非常相似(Hill，2001)，但是它包含了非常完整的弹性波矢量场信息.
由于本文所做向下反向延拓只对单程的上行波有效，并且实际观测的有效反射地震信号也是上行波，因此本文只对上行波做了分析. 实际观测到的数据中是存在下行波的，但人们往往将其看作噪声，比如多次波等.
3.4.2 海底如图 6所示，当观测边界 Σ 以上为流体介质时，模型称为海底模型. 海底模型的多分量地震数据是三种模型中最难以处理的，它没有固定的应力边界条件，也没有相互独立的多波型波场传播属性. 为了获得准确的边界条件，Ravasi和Curtis(2013)利用了同时记录海底质点速度矢量和声波水压的4C地震数据. 在没有声波水压时，多波型难以准确分离和成像. 假设观测数据只有上行波信号，海底界面上具有精确的应力与位移关系式. 这种假设与我们所关注的地下反射信号是一致的.
|
图 6 海底上行波 Fig. 6 Upgoing waves in ocean-bottom model |
下面我们对应力边界条件进行详细分析. 应力分析只需在2D空间进行，并且很容易推广到3D空间. 对于2D空间的P、SV波场，海底界面满足关系式
|
(24) |
式中上标 w 代表边界 Σ 以上的海水(或流体介质)， σzz和τzx 分别表示正应力和剪切应力， uz 表示位移 z 分量. 如果假设只有上行波入射到边界 Σ， 如图 7所示，在边界上会同时激起反射P波、反射SV波和透射声波. 这里的入射波可以是P波、SV波，也可以是P、SV波的混合.
如果入射波为平面波，则慢度矢量满足Snell定律 pxw=pxP=pxS=pxinc， 其中上标inc代表入射波. 设定边界 Σ(z=0) 上下的平面波位移为
|
(25) |
式中 u w=(uzw，uxw) 只包含上行声波成分， u =(uz，ux) 同时包含上行和下行P、SV波.
流体中位移矢量 u w 具有无旋性质，即
|
(26) |
由(26)式可得
|
(27) |
或者表示为
|
(28) |
由此可知，位移水平分量 uxw 虽然不满足边界连续性，即 ux≠uxw， 但是它可用 uzw 表示.
将正应力 σzzw 表示成位移表示形式，并将式(28)代入，可得
|
(29) |
根据边界条件(24)，可得
|
(30) |
2D空间的SH波具有 uz=0和σzz=0， 因此，SH波也满足式(30). 式(30)也很容易推广到3D空间. 这样，我们获得了海底模型在上行假设情况下的应力和位移关系式. 根据广义虎克定律(7)，整理可得 uz 的空间偏导数表达式
|
(31) |
其中
|
式中， pxw和pyw 为流体上行波慢度水平分量，上标“，”以后的数字表示指数，如 pxw，2=(pxw)2. cw 为流体声波波速. 另外，由于上行波方向与积分方程(15)慢度矢量方向相反，所以有 pxw=－px，pyw=－py，其中px和py 为积分方程(15)慢度分量.
利用式(24)中的零剪切应力边界条件，式(15)可以整理为
|
(32) |
在高频近似情况下，局部平面波Dn的偏导数只需考虑最大奇异项的求导，忽略非奇异项. 因此Dn近似满足偏导数表达式(31). 将式(31)代入式(32)，整理可得
|
(33) |
设定
|
(34) |
其中 Ď ν 为标量局部平面波，是 Dn 的加权求和. Ď ν 添加“ ^ ”符号是为了与式(20)中的 Dν 区分开，两者在推导原理和表示形式上都有很大差别. Wnv为3×3 权重矩阵
|
(35) |
Wnv 是海底模型多波型波场延拓的关键因子，它同时与应力和位移边界条件有关. 这样，式(33)可以简化为
|
(36) |
式(36)即为所获得的海底传播模型的波场延拓积分方程.
因此，在上行波假设情况下，可以获得准确的海底应力与位移关系式，实现准确的海底波场延拓.
3.4.3 自由地表如图 7所示，当观测边界 Σ 以上为真空时，模型称为自由地表模型. 当前，多数多分量偏移方法均假设模型为自由地表模型. 自由地表边界只 存在上行波入射，不存在下行波入射，如图 7所示，并且满足正应力和剪切应力均为零，即 σzz=τzx=τzy=0.
|
图 7 自由地表上行波 Fig. 7 Upgoing waves in free-surface model |
自由地表模型波场延拓与海底模型非常类似，相当于海底模型的一种极限情况. 如果对海底模型参数取 ρw→0，cw→0， 就等价于自由地表模型. 海底模型波场延拓积分方程(36)也适用于自由地表模型，这时的权重矩阵 Wnv 为
|
(37) |
弹性波场为矢量波场，常规的互相关成像条件不适用于矢量波成像，并且PS转换波存在极性翻转问题. 为了解决矢量波成像问题，本文利用散度/旋度互相关成像条件. 该方法借用了弹性逆时偏移常用的Helmholtz分解方法(Yan and Sava，2008)，但目的不是做P、S波分解. 各向同性介质中弹性波位移矢量满足关系式
|
(38) |
其中， Φ和Ψ 分别为标量位场和矢量位场. P、 S 有具体物理意义， P 为体应变， S 为旋转应变. 对于2D空间的SV波， S 为单分量. 由于常规数据采集震源以P波震源为主，因此，本文仅列出了PP和PS偏移成像表达式. 对 P或 S 做互相关，可得成像表达式：
|
(39a) |
|
(39b) |
式中，算子 Δ x 下标表示作用于点 x，Sgn(x) 为一单位矢量函数，并且有表达式
|
(40) |
在弹性各向同性介质中，PS转换波的旋度矢量具有同时垂直于入射波和转换波传播方向的性质. 由此可见，单位矢量 Sgn(x) 代表PS转换波的一种旋转方向，如图 8所示. 本文将其约定为转换波的正向旋转方向. 在3D空间中，矢量波场的旋度仍然是矢量，将其投影到 Sgn(x) 方向可获得标量场，同时解决了3D空间的极性翻转问题. 在2D空间， Sgn(x) 为单分量，起到符号函数的作用. 当慢度矢量 p r(x)和p s(x) 同向时， Sgn(x) 出现0/0无意义现象. 但是此时P波垂直入射地下界面，不会形成转换波，因此可以取其为 0 矢量.
|
图 8 PS转换波的传播方向与旋转方向 Fig. 8 Propagation and rotation directions of PS converted waves |
将互相关成像条件写成频率域表示形式：
|
(41a) |
|
(41b) |
利用点源正向延拓、多分量地震数据反向延拓和矢 量场成像条件，可整理得到偏移成像表达式.
由式(5)可知，P波点源正向延拓可表示为
|
(42) |
式中下标 s 表示震源位置 s .
如果将弹性高斯束写成标量和偏振方向矢量的表示形式：
|
(43) |
则弹性高斯束的散度和旋度有以下关系式：
|
(44a) |
|
(44b) |
|
(44c) |
(44)式将在互相关成像公式推导中应用.
4.1 自由空间将式(23)和(42)代入矢量场互相关成像条件(41)，并调换积分和求导顺序，再将式(44)代入，可得自由空间模型PP反射波和PS转换波成像公式
|
(45) |
|
(46) |
成像公式(46)是PS1、PS2成像的合并表示，即 IPS=IPS1+IPS2， 式中的标量 uGBS 满足关系式 uGBS=uGBS1=uGBS2，D S 为矢量场，
|
(47) |
为了方便公式推导和计算实现，文中将S波的复杂偏振情况分解为S1和S2两种线性偏振. 但是S1与S2是同一种波型，除了偏振方向相互垂直以外，具有相同的波速分布和波场扩散等性质. S1和S2的划分与坐标系的设置有关(Červený，2001)，属于人为行为. 这与SV和SH波的划分并不相同. 在弹性各向同性介质中，PS转换波仅存在P-SV转换波，不存在P-SH转换波，因此在PS成像中真正有效的是P-SV转换波成像. 由于SV转换波中同时包含S1和S2波成分，所以S1和S2的合并成像才能完整地反映P-SV成像效果. 成像公式(45)和(46)与声波高斯束偏移成像公式(Hill，2001)相比可以看出，整理后的弹性波偏移成像并不比声波偏移过于复杂，其整体与声波偏移具有非常相似的表示形式.
4.2 海底与自由空间模型类似，将式(36)和(43)代入矢量场互相关成像条件(41)，并调换积分和求导顺序，再将式(44)代入，可得海底模型PP反射波和PS转换波成像公式：
|
(48) |
|
(49) |
式中， Ď S 为矢量场，
|
(50) |
自由地表模型与海底模型有相同的偏移成像公式，在此不单独列出.
5 数值试验
本文利用2D和3D模型数据的偏移成像结果，来验证所提供方法的准确性. 2D情况相当于3D的一种特殊情况. 使3D空间介质参数沿y轴方向保持不变，并将3D格林函数转换为2D格林函数，即可简化出2D偏移公式. 为简单起见，文中模型仅显示P波波速，并通过取P、S波速比为
简单水平层状介质的P波速度模型如图 4所示. 反射界面在深度z=1000 m，并且界面上下P波波速分别为2000 m·s-1和2400 m·s-1. 炮点和检波点均位于边界面 Σ(z=0)， 并且炮点位于x=2000 m，检波点位于x=0~4000 m，间隔10 m，共401点. 基于此简单模型的参数取值，分别对自由空间、海底和自由地表三种边界模型做了记录合成，即检波点以上分别设为完全匹配吸收、海水和真空三种不同介质.
图 9显示了自由空间和海底模型的记录合成图，自由地表合成记录与海底类似，没有单独显示. 由图可以看出，海底模型记录比自由空间模型记录明显复杂，海底所接收到的上行反射波同时包含PP、PS、SP、SS四种反射波. 在多分量数据的实际应用中，一般认为PP和PS是我们所关注的有效信息，SP和SS为干扰噪声. 并且海底模型记录不具有P直达波z分量极弱的特点，这与海底多波耦合关系有关. 在多分量偏移处理之前，海底直达波和面波利用倾斜滤波器得以压制，而由S波激起的SP和SS波很难与P波激起的PP和PS波彻底分离. 并且在水平层状介质中，SP和PS波同相轴完全重合，更加难以分离.
|
图 9 自由空间和海底模型的记录合成图 (a)和(b)分别为自由空间模型x和z分量；(c)和(d)分别为海底模型x和z分量. Fig. 9 Synthetic records of free-space and ocean-bottom models (a)and(b)are x and z components of free-space model，respectively；(c)and(d)are x and z components of ocean-bottom model，respectively. |
不同模型记录利用不同偏移方法进行了多分量偏移成像. 三种记录利用了相同的偏移模型. 图 10显示了自由空间和海底模型记录分别利用各自的偏 移方法所得到的偏移成像结果. 由图 10a和10b可以看出，自由空间偏移方法能够很准确地将P和S波分离并准确成像，PS成像的极性翻转问题也得到了解决. 图 10d所示的海底模型PS成像说明海底偏移方法也能准确实现海底偏移成像. 图 10c存在SP波对PP成像所形成的明显的串扰. 这是由于本文所提供的多分量高斯束偏移方法，能够准确实现海底边界上的P、S上行波自动分离，但是，其中的P上行波同时包含PP反射波和SP转换波成分. 因此，在PP偏移成像中，P上行波中所包含的PP、SP同时参加互相关成像，由此产生了SP Crosstalk噪声. Zhe和Greenhalgh(1997)研究了PP成像和SS成像，并且讨论了PP、PS、SP、SS波之间的串扰影响.
|
图 10 自由空间和海底模型记录偏移成像结果 (a)和(b)分别为自由空间模型PP和PS成像;(c)和(d)分别为海底模型PP和PS成像. Fig. 10 Migration results of free-space and ocean-bottom records (a)and(b)are PP and PS images of free-space model，respectively;(c)and(d)are PP and PS images of ocean-bottom model，respectively. |
为了比较三种模型偏移成像之间的差别，下面我们对每一种模型记录均做三种偏移方法成像(为了方便，仅对PS转换波成像进行对比). 图 11显示了三种偏移方法对三种模型记录成像效果对比图. 图中共有9张PS波成像图片，以3×3形式排列. 其中，上行、中行和下行分别为自由空间、自由地表和海底模型记录的成像结果，左列、中列和右列分别为自由空间、自由地表和海底偏移方法的成像结果. 因此，只有主对角线上的3个成像准确利用了各自的偏移方法. 由图可以看出，也只有这3个成像得到了非常准确的波型分离和偏移成像，其他成像均存在不同程度的P波成分对PS成像的干扰. 由此可以得到，应力边界条件的分析与应用为多分量偏移成像带来了准确的波型自动分离和偏移成像. 另外，由于自由地表和海底模型记录只做了简单的倾斜滤波，测线两端数据的突然截断导致边界处的滤波处理不准确，使成像结果存在明显划弧现象.
|
图 11 三种模型记录利用三种偏移方法的成像效果对比 Fig. 11 Comparison of imaging results on three model records using three migration methods |
2D弹性SEG/EAGE盐丘模型由原始的声波SEG/EAGE盐丘模型改进和扩展而来(Xie and Wu，2005). 为了使模型适用于弹性波场合成和偏移，对声波模型做了几处修改和扩展. 将声波波速模型作为弹性P波波速模型，并将原始上表面的海水介质波速由海底P波波速延伸替代. 图 12显示了该弹性模型的P波波速模型. 为了计算方便，有限差分记录合成利用了完全匹配吸收方法，相当于自由空间模型. 这种情况下，不存在海底(地表)中的下行S波传播情况，所以在PP成像中不存在SP Crosstalk噪声. 合成的记录共有325炮，炮间距为48.768 m. 单炮记录中利用326个检波点接收信号，检波点间隔为12.192 m，偏移距从-1981.2 m至1981.2 m. 图 13显示了第151炮的单炮记录，该记录炮点位于x=7315.2 m. 由图可以看出，记录中的直达波z分量并不为零，并且随偏移距的增加，直达波z分量逐渐增强. 这是由于在该模型中直达波总是以回折波的形式传播，而不是沿直线方向传播.
|
图 12 2D弹性SEG/EAGE盐丘模型的P波波速 Fig. 12 P wave speed of 2D elastic SEG/EAGE salt model |
|
图 13 2D弹性SEG/EAGE盐丘模型第151炮单炮记录 (a)和(b)分别为记录的x和z分量 Fig. 13 Record of the 151st shot in 2D elastic SEG/EAGE salt model (a)and(b)are x and z components of the record，respectively. |
多分量偏移模型是由弹性P、S波波速模型平滑处理后得到的. 依据记录合成情况，选择自由空间偏移方法成像. 图 14显示了该模型记录的多分量偏移成像结果，由图可见，多分量高斯束偏移方法对PP反射波和PS转换波均能准确成像，对盐丘掩盖区域成像和高度倾斜界面成像都有很好的效果. 比较PP和PS成像效果可以看出，PS成像比PP成像有明显较高的分辨率. 另外，由于S波波速较低，波长较短，PS成像可能会形成一些假频噪声，尤其在浅层成像部分.
|
图 14 2D弹性SEG/EAGE盐丘模型偏移成像结果 (a)和(b)分别为PP和PS成像. Fig. 14 Migration imaging results of 2D elastic SEG/EAGE salt model (a)and(b)are PP and PS images，respectively. |
为了验证该方法在3D空间的适用性，以及3D空间S1和S2合并成像的可行性，本文提供了3D单界面水平层状模型的单炮偏移成像结果. 层状模型x轴和y轴坐标范围均为0~4000 m，界面深度1000 m，界面上下P波波速分别为2000 m·s-1和2200 m·s-1. 炮点坐标为(x，y)=(2000 m，2000 m). 检波点位于x=0~4000 m，y=0~4000 m范围内，间隔dx=dy=10 m，以401×401阵列排列.
图 15为该模型的PP成像效果图，其中图 15c中的L1和L2线即为图 15a和15b的截取位置. 由图 15a可见，该方法准确实现了P、S波型分离. 图 15b中成像位置上下均为测区突然截断导致的划弧现象. 因此，该方法所获得的PP成像是很合理的.
|
图 15 3D简单层状模型PP成像 (a)x=1450 m剖面；(b)y=2000 m剖面；(c)z=1000 m切片. Fig. 15 PP imaging of 3D simple layering model (a)Profile x=1450 m；(b)Profile y=2000 m；(c)Slice z=1000 m. |
图 16显示了3D多分量高斯束偏移PS成像效果图. 该PS成像即为本文所提出的S1和S2合并成像.取值位置与图 15完全相同. 由图 16c可以看出，PS成像的振幅与方位角无关. 图中A所标示的为弱成像区域，位于P波垂直入射点附近，为P波小角度入射区域. 该区域的PS反射系数很小. 因此，该成像效果符合实际情况. 另外，该偏移方法的成像效果与3D中心射线坐标系的设定无关，即修改该坐标系的设定不会影响PS成像. 由此可以证明S1和S2合并成像的可行性.
|
图 16 3D简单层状模型PS成像 (a)x=1450 m剖面;(b)y=2000 m剖面;(c)z=1000 m切片. Fig. 16 PS image of 3D simple layering model (a)Profile x=1450 m;(b)Profile y=2000 m;(c)Slice z=1000 m. |
图中B所指向的多个同相轴是由于测区范围边界的突然截断导致的划弧现象，并非P波串扰噪声. 与2D情况明显不同，3D空间的边界划弧为弧 面，并且沿x轴和y轴方向均会产生划弧面. 而P 波串扰噪声在该PS成像效果中得到了很好的压制. 另外，由图 15和图 16相比可知，PS成像比PP成像具有更大的有效成像范围.
由此可见，本文所提供的3D多分量偏移方法既可实现P、S波的准确自动分离，也可获得PP和PS合理的成像效果.
6 结论我们对自由空间、海底和自由地表三种模型的应力边界条件进行了详细分析. 虽然应力边界条件可能是未知的，但是在上行传播假设情况下，应力与位移有明确的关系式. 在自由空间模型中，多波传播相互独立，没有耦合. 利用波场传播和偏振方向等几何属性可以实现多波型的分离. 在海底模型中，依据应力与位移在海底界面的连续性，可以建立上行波应力位移关系式. 而自由地表模型相当于海底模型的一种特殊情况，海底的各种表达式均适用于自由地表. 由于波场均可分解为平面波的叠加，远源波场也可看作局部平面波，因此我们对独立的(局部)平面波进行应力分析. 在此基础上准确推导了3D多波多分量高斯束波场延拓表达式，同时准确实现了多波型的自动分离.
为了解决常规互相关成像条件不适用于矢量场成像的问题，我们采用了散度/旋度互相关成像条件. 通过约定PS转换波的正向旋转方向，解决了3D空间PS转换波复杂的偏振和旋转情况所带来的成像问题. 由此，获得了准确的适用于3D弹性介质空间中的多波多分量高斯束叠前深度偏移成像公式，其表达形式与声波高斯束偏移公式相似，且包含了非常完整的弹性波信息. 本文利用2D和3D模型数据的成像结果验证了该方法的准确性.
| Alkh alifah T. 1995. Gaussian beam depth migration for anisotropic media. Geophysics , 60(5): 1474–1484. | |
| Berkhout A J, Wapenaar C P A. 1989. One-way versions of the Kirchhoff integral. Geophysics , 54(4): 460–467. | |
| Bleistein N, Cohen J K, Stockwell J W. Multidimensional Seismic Imaging, Migration, and Inversion . (in Chinese) Beijing: Science Press, 2004 . | |
| ČervenýV, PopovM M, PšenčíkL. 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach. Geophys. J. Int. , 70(1): 109–128. | |
| ČervenýV, PšenčíkL. 1983. Gaussian beams and paraxial ray approximation in three-dimensional elastic inhomogeneous media. J. Geophys. , 53: 1–15. | |
| ČervenýV. Seismic Ray Theory. Cambridge University Press, 2001 . | |
| Chang W F, McMechan G A. 1987. Elastic reverse-time migration. Geophysics , 52(10): 1365–1375. | |
| Chang W F, McMechan G A. 1994. 3-D elastic prestack reverse-time depth migration. Geophysics , 59(4): 597–609. | |
| Dai T F, Kuo J T. 1986. Real data results of Kirchhoff elastic wave migration. Geophysics , 51(4): 1006–1011. | |
| Davis T L. 2001. Multicomponent seismology-The next wave. Geophysics , 66(1): 49. | |
| Druzhinin A. 2003. Decoupled elastic prestack depth migration. J. Appl. Geophys. , 54(3-4): 369–389. | |
| Duan P F, Cheng J B, Chen A P, et al. 2013. Local angle-domain Gaussian beam prestack depth migration in a TI medium. Chinese J. Geophys. , 56(12): 4206–4214. doi: 10.6038/cjg20131223. | |
| Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics , 70(4): S71–S77. | |
| Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics , 74(2): S11–S23. | |
| Hill N R. 1990. Gaussian beam migration. Geophysics , 55(11): 1416–1428. | |
| Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics , 66(4): 1240–1250. | |
| Hokstad K. 2000. Multicomponent Kirchhoff migration. Geophysics , 65(3): 861–873. | |
| Huang Z Y, Sun J K, Zhu S J, et al. Multicomponent Seismic Technology . (in Chinese) Beijing: Petroleum Industry Press, 2007 . | |
| Kuo J T, Dai T F. 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver. Geophysics , 49(8): 1223–1238. | |
| Li X L, Mao W J. 2015. Multicomponent Gaussian beam migration in elastic medium.//77th EAGE Conference and Exhibition. EAGE. | |
| Li X Y. 1997. Fractured reservoir delineation using multicomponent seismic data. Geophys. Prospecting , 45(1): 39–64. | |
| Li X Y, Yuan J, Ziokowski A, et al. 1999. Estimating Vp/Vs ratio from converted waves-A 4C case example.//61st EAGE Conference and Exhibition. EAGE. | |
| Mittet R. 1994. Implementation of the Kirchhoff integral for elastic waves in staggered-grid modelling schemes. Geophysics , 59(12): 1894–1901. | |
| Nowack R L, Sen M K, Stoffa P L. 2003. Gaussian beam migration for sparse common-shot and common-receiver data.//73rd Annual International Meeting, SEG Expanded Abstracts. Dallas, Texas, 1114-1117. | |
| Pao Y H, Varatharajulu V. 1976. Huygens' principle, radiation conditions, and integral formulas for the scattering of elastic wave. J. Acoust. Soc. Am. , 59(6): 1361–1371. | |
| Qian Z P, Chapman M, Li X Y, et al. 2007. Use of multicomponent seismic data for oil-water discrimination in fractured reservoirs. The Leading Edge , 26(9): 1176–1184. | |
| Ravasi M, Curtis A. 2013. Elastic imaging with exact wavefield extrapolation for application to ocean-bottom 4C seismic data. Geophysics , 78(6): S265–S284. | |
| Rechtien R D. 1985. On "Kirchhoff elastic wave migration for the case of noncoincident source and receiver," by John R. Kuo and Ting-fan Dai (GEOPHYSICS , 49: 1223–1238. | |
| Schleicher J, Tygel M, Ursin B, et al. 2001. The Kirchhoff-Helmholtz integral for anisotropic elastic media. Wave Motion , 34(4): 353–364. | |
| Sena A G, Toksoz M N. 1993. Kirchhoff migration and velocity analysis for converted and nonconverted waves in anisotropic media. Geophysics , 58(2): 265–276. | |
| Sun R, McMechan G A. 1986. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles. Proceedings of the IEEE , 74(3): 457–465. | |
| Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics , 66(5): 1519–1527. | |
| Sun R, McMechan G A, Lee C S, et al. 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data. Geophysics , 71(5): S199–S207. | |
| Wapenaar C P A, Haimé G C. 1990. Elastic extrapolation of primary seismic P-and S-waves. Geophys. Prospecting , 38(1): 23–60. | |
| Xie X B, Wu R S. 2005. Multicomponent prestack depth migration using the elastic screen method. Geophysics , 70(1): S30–S37. | |
| Xu S Y, Li Y P, Ma Z T. 1999. Separation of P-and S-wave fields via the τ-q transform. China Offshore Oil and Gas (Geology) (in Chinese) , 13(5): 334–337. | |
| Yan J, Sava P. 2008. Isotropic angle domain elastic reverse-time migration. Geophysics , 73(6): S229–S239. | |
| Yan J, Sava P. 2009. Elastic wave-mode separation for VTI media. Geophysics , 74(5): WB19–WB32. | |
| Yue Y B. 2011. Study on Gaussian beam migration methods in complex medium[Ph. D. thesis] (in Chinese). Qingdao:China University of Petroleum (East China). | |
| Zhe J P, Greenhalgh S A. 1997. Prestack multicomponent migration. Geophysics , 62(2): 598–613. | |
| Zhu T F, Gray S, Wang D L. 2006. Prestack Gaussian-beam depth migration in anisotropic media.//76th Annual International Meeting, SEG Expanded Abstracts. New Orleans, Louisiana, 2362-2366. | |
| 布莱斯坦N, 科恩JK, 斯托克韦尔JW. 多维地震成像、偏移和反演中的数学. 北京: 科学出版社, 2004 . | |
| 段鹏飞, 程玖兵, 陈爱萍, 等. 2013. TI介质局部角度域高斯束叠前深度偏移成像. 地球物理学报 , 56(12): 4206–4214. | |
| 黄中玉, 孙建库, 朱仕军, 等. 多分量地震技术. 北京: 石油工业出版社, 2007 . | |
| 许世勇, 李彦鹏, 马在田. 1999. τ-q变换法波场分离. 中国海上油气(地质) , 13(5): 334–337. | |
| 岳玉波. 2011. 复杂介质高斯束偏移成像方法研究[博士论文]. 青岛:中国石油大学(华东). | |