引言

随着勘探技术的发展，多分量地震技术逐渐受到业界的重视(孙歧峰，杜启振，2011)． 为了充分利用多分量资料的纵横波信息，处理流程通常都包括波场分离环节(董敏煜，2002)． 依据不同的处理方法，可将波场分离技术分为两大类(李庆忠，王建花，2000). 第一类是基于标量波场理论的波场分离技术，即对地震记录实现纵波与横波的波场分离． 工业界通常将多分量地震记录的垂直分量视为P波，并利用坐标旋转得到径向(R)和切向(T)分量，然后采用相对成熟的声波处理流程对上述分量分别进行处理． 除此之外，还可利用极化滤波(Etgen，1988； 李志远等，2013)、 f-k滤波(Amundsen et al，1998)、 τ-P变换(许世勇等，1999)、 拉东(R and on)变换(冯晅等，2011)、 平面波分解(Devaney，Oristaglio，1986)等方法实现对地震记录纵横波的波场分离． 第二类是基于矢量波场理论的分离技术，即在波场传播过程中实现对波场快照的波场分离． 在各向同性介质方面，Sun和McMechan(2001)基于亥姆霍兹(Helmholtz)分解理论提出利用散度和旋度算子实现P波与S波的波场分离，但会引入90°相移． 针对这一问题，Sun等(2001)利用希尔伯特(Hilbert)变换实现了相位校正． 相对于各向同性介质，各向异性介质波场分离的难点在于质点振动方向与波场传播方向并不完全平行或正交，因此不能直接利用散度和旋度算子进行波场分离(Yan，Sava， 2009a，b)． Dellinger和Etgen(1990)依据克里斯托弗尔(Christoffel)方程，在波数域中将弹性波场向偏振方向投影，从而得到qP波和qS波． 在此基础上，Yan和Sava(2009a，b)进一步计算得到空间域不同尺寸空间窗的波场分离算子，并利用该算子与波场值卷积得到qP波和qS波． 此外，Zhang和McMechan(2010)还提出了一种在波数域获得每个笛卡尔坐标分量上保幅qP波和qS波的波场分解方法． 与第一类波场分离技术相比，第二类波场分离技术更能够描述波场在地下介质中传播的矢量特性． 因此，针对第二类波场分离技术展开研究具有重要意义．

弹性波场数值模拟是实现第二类波场分离的前提． 有限差分法由于其计算效率高、 占用内存小以及便于实现等优点而成为弹性波场数值模拟的重要工具． 相对于中心规则网格有限差分法，交错网格有限差分法(Madariaga，1976)的精度更高且频散压制效果更好，因此被广泛应用于各种介质的波场数值模拟中. 例如： Virieux(1984，1986)利用交错网格有限差分法模拟了各向同性介质中SH波、 P波及SV波的传播； Carcione和Helle(1999)利用交错网格有限差分法实现了孔隙介质中弹性波场的数值模拟； 裴正林(2004)利用交错网格高阶有限差分法实现了三维各向异性介质弹性波数值模拟． 但交错网格有限差分法的不足之处在于质点振动的速度分量定义在不同的网格点上，利用这种差分格式进行数值模拟得到的波场值不能够直接应用到Yan和Sava(2009a，b)所提出的各向异性介质波场分离方法中． 而旋转交错网格(Saenger et al，2000)将所有的应力分量定义在同一网格点上，将所有的质点振动速度分量定义在另一网格点上，且由旋转交错网格有限差分得到的波场值能够直接应用到各向异性介质波场分离过程中． 同时，采用旋转交错网格有限差分进行数值模拟时，弹性模量无需进行插值处理(严红勇，刘洋，2012)． 综上可见，旋转交错网格有限差分更适于各向异性介质波场数值模拟及其波场分离过程．

本文拟以具有垂直对称轴的横向各向同性(transversely isotropic with a vertical axis，简写为VTI)介质的波场分离为例，从波场分离的实现过程出发，对比分析交错网格和旋转交错网格的参数定义方式对VTI介质波场分离(Yan，Sava，2009a)造成的影响，并通过数值试验对其进行分析验证，以进一步说明旋转交错网格在各向异性波场数值模拟及其波场分离方面的优势.

1.   VTI介质一阶速度-应力方程及两种网格定义方式

二维情况下，VTI介质一阶速度-应力方程(Virieux，1986)为

式中： v_x和v_z分别为x方向和z方向的质点速度分量； τ_xx和τ_zz为正应力，τ_xz为切应力； C₁₁，C₁₃，C₃₃和C₅₅为刚度系数； ρ为密度． 分别采用交错网格有限差分(Madariaga，1976)和旋转交错网格有限差分(Saenger et al，2000)对式(1)进行数值模拟．

图 1给出了交错网格(Virieux，1986)和旋转交错网格(Saenger et al，2000)应力场和速度场的定义方式示意图，其中： 实线表示定义在时间上的主网格，虚线表示在时间t上交错了Δt/2的半网格； O点为参考点，P₁点和P₂点分别代表某一点x方向和z方向上的质点振动速度分量． 比较图 1a与图 1b可知，质点振动的速度分量在交错网格中没有定义在同一个网格节点上，而在旋转交错网格中则定义在同一个网格节点上．

图  1  交错网格(a)和旋转交错网格(b)示意图

Figure  1.  Schematic diagram of staggered-grid (a) and rotated staggered-grid (b)

下载: 全尺寸图片 幻灯片

2.   VTI介质波场分离过程

VTI介质中，qP波与qS波的偏振方向相互垂直，利用偏振矢量将波场投影至偏振方向可得到qP波和qS波(Dellinger，Etgen，1990)． 在波数域中，通过求解克里斯托弗尔方程，可得qP波和qS波的偏振矢量. 二维情况下求解

可得qP波的偏振矢量U_qP为

qSV波的偏振矢量U_qSV为

This page contains the following errors:

Below is a rendering of the page up to the first error.

图  2  空间域波场分离算子. (a) x方向; (b) z方向

Figure  2.  The separation operators of wave-field in space domain (a) x direction; (b) z direction

下载: 全尺寸图片 幻灯片

3.   两种交错网格对波场分离的影响

这里以qP波的波场分离过程为例，详细阐述两种交错网格对VTI介质波场分离所产生的影响．

波数域和空间域中的qP波(Yan，Sava，2009a)分别表示为

This page contains the following errors:

Below is a rendering of the page up to the first error.

波数域和空间域中各速度分量的傅里叶变换对为

如图 1a所示，在交错网格中，O点处的质点振动速度分量可利用P₁和P₂点处的质点振动速度分量表示为

如图 1b所示，在旋转交错网格中，O点处的质点振动速度分量可利用P₁和P₂点处的质点振动速度分量表示为

将式(9)和式(10)带入式(6)，可得

此即为空间域中的qP波.

Yan和Sava(2009a)提出的VTI介质波场分离方法是将V_x(k_x，k_z)和V_z(k_x，k_z)以及v_x(x，z)和v_z(x，z)视为同一个网格节点上的质点振动速度分量进而实现波场分离． 对比式(11)与式(12)可得，由于交错网格质点振动速度分量未定义在同一个网格节点上，导致用于VTI介质波场分离的x方向质点振动速度分量比z方向多出平移算子exp(-i2πk_xΔx/2)，z方向质点振动速度分量比x方向多出平移算子exp(-i2πk_zΔz/2)，从而造成波场分离不彻底. 针对这一问题，Zhang和McMechan(2010)利用两点平滑估算其中间值，从而得到同一个网格节点上的质点振动速度分量，但该方法只是对波场中间值的近似，并不能完全表示两点中间的波场值． Du等(2014)提出在波数域中对波场值进行平移校正，从而实现两点中间值的准确求取； 但鉴于本文采用的VTI介质波场分离方法(Yan，Sava，2009a)的处理流程是在空间域中进行的，这种波数域中的插值方法将会造成额外的计算量，故不适用于本文．

如式(12)所述，由旋转交错网格有限差分得到的用于波场分离的x方向与z方向的质点振动速度分量在网格节点上的对应位置保持一致，不需要进行额外的插值处理，同时可避免因波场插值造成的波场值不准确及效率降低等问题． 另外，在使用旋转交错网格进行波场数值模拟的过程中也不需要对弹性模量进行插值处理.

4.   数值算例

本节将采用均匀VTI介质模型进行弹性波场数值模拟和波场分离，通过模拟测试结果对比分析交错网格和旋转交错网格对VTI介质波场分离(Yan，Sava，2009a)造成的影响. 本文采用的模型网格大小为600×300，水平方向和垂直方向的空间采样间隔均为10 m，时间采样间隔为0.001 s，模型参数如表 1所示． 交错网格有限差分和旋转交错网格有限差分均采用空间8阶、 时间2阶进行弹性波场数值模拟，利用非分裂完全匹配层(non-splitting perfect matched layer，简写为NPML)对边界反射进行吸收衰减． 震源子波选取主频为35 Hz的雷克子波，其解析式为f(t)=[1－2(πf_pt)²]exp[－(πf_pt)²]，式中f_p为子波主频； 震源位置坐标为(2800 m，500 m)．

表  1  均匀VTI介质参数表

Table  1.  The parameters for homogeneous VTI media

δ	ε	vP0/(m·s-1)	vS0/(m·s-1)	ρ/(g·cm -3)
-0.10	0.19	3290	1728	2000

下载: 导出CSV 

| 显示表格

在弹性波场数值模拟过程中，提取第400个时间步长(即第400 ms)的质点振动速度x分量与z分量的波场快照进行波场分离，以阐明这两种交错网格对VTI介质波场分离所造成的影响．

图 3和图 6分别为交错网格和旋转交错网格波场数值模拟得到的波场快照. 对比图 3与图 6可知，采用两种交错网格有限差分进行波场数值模拟均能得到清晰的波场快照． 图 4给出了直接利用交错网格波场数值模拟得到的波场快照进行波场分离得到的qP波和qS波； 图 5给出了利用交错网格波场数值模拟并进行两点平均插值处理后的波场快照进行波场分离得到的qP波和qS波(Zhang，McMechan，2010)． 可以看出，经过两点插值处理后进行波场分离得到的结果要优于未采用波场插值处理直接进行波场分离的结果； 但如图 4和图 5中箭头所指位置所示，qP波中仍然残留有qS波，qS波中同样残留有qP波，且qS波中的残留现象更为严重． 图 7为旋转交错网格波场数值模拟得到的波场快照进行波场分离的结果． 比较图 4、 图 5及图 7中箭头和方框所示可知，使用旋转交错网格数值模拟得到的波场值进行波场分离得到的qP波和qS波较为干净彻底，且qP波的连续性较好．

图  3  交错网格数值模拟得到的波场快照(400 ms). (a) x分量； (b) z分量

Figure  3.  Snapshots by using staggered-grid numerical simulation scheme (400 ms). (a) x-component; (b) z-component

下载: 全尺寸图片 幻灯片

图  4  交错网格数值模拟得到的波场快照直接进行波场分离的结果. (a) qP波; (b) qS波

Figure  4.  Separation results by using the snapshots directly from staggered-grid numerical simulation scheme. (a) qP-wave; (b) qS-wave

下载: 全尺寸图片 幻灯片

图  5  交错网格数值模拟得到波场快照并进行两点平均插值处理后进行 波场分离得到的结果. (a) qP波; (b) qS波

Figure  5.  Separation results by using the snapshots from staggered-grid numerical simulation scheme with two points average interpolation. (a) qP-wave; (b) qS-wave

下载: 全尺寸图片 幻灯片

图  6  旋转交错网格数值模拟得到的波场快照(400 ms). (a) x分量； (b) z分量

Figure  6.  Snapshots by using rotated staggered-grid numerical simulation scheme (400 ms). (a) x-component; (b) z-component

下载: 全尺寸图片 幻灯片

图  7  旋转交错网格数值模拟得到的波场快照直接进行波场分离得到的结果. (a) qP波; (b) qS波

Figure  7.  Separation results by using the snapshots obtained from rotated staggered-grid numerical simulation. (a) qP-wave; (b) qS-wave

下载: 全尺寸图片 幻灯片

为了进一步说明使用旋转交错网格波场数值模拟得到的波场值进行波场分离的优势，本文采用如图 8所示的VTI介质Hess模型进行波场数值模拟和波场分离． 模型网格大小为600×350，水平方向和垂直方向的采样间隔为10 m，时间采样间隔为0.001 s； 采用主频为35 Hz的雷克子波，震源位置坐标为(3200 m，700 m)； 交错网格有限差分和旋转交错网格有限差分均采用空间8阶、 时间2阶进行弹性波场数值模拟； 采用非分裂完全匹配层对边界反射进行吸收衰减，同时选取600 ms的波场快照进行波场分离． 图 9和图 10分别给出了利用交错网格和旋转交错网格进行数值模拟得到的波场快照以及对波场快照进行波场分离的结果，对比其中红色箭头标示处可验证，使用旋转交错网格波场数值模拟得到的波场值进行波场分离得到的qP波与qS波分离效果良好．

图  8  VTI 介质Hess模型. (a) δ; (b) ε; (c) ρ; (d) v_P0; (e) v_S0

Figure  8.  Hess model for VTI media. (a) δ; (b) ε; (c) ρ; (d) v_P0; (e) v_S0

下载: 全尺寸图片 幻灯片

图  9  Hess模型交错网格波场数值模拟结果及波场分离结果 (a) x分量波场快照； (b) z分量波场快照； (c) qP波； (d) qS波

Figure  9.  Snapshots for Hess model by using staggered-grid numerical simulation and the separation results. (a) x-component snapshot; (b) z-component snapshot; (c) qP-wave; (d) qS-wave

下载: 全尺寸图片 幻灯片

图  10  Hess模型旋转交错网格波场数值模拟结果及波场分离结果 (a) x分量波场快照； (b) z分量波场快照； (c) qP波； (d) qS波

Figure  10.  Snapshots for Hess model by using rotated staggered-grid numerical simulation and the separation results. (a) x-component snapshot； (b) z-component snapshot； (c) qP-wave； (d) qS-wave

下载: 全尺寸图片 幻灯片

5.   讨论与结论

本文针对交错网格和旋转交错网格对VTI介质波场分离造成的影响进行了详细分析，并选取简单模型和复杂模型进行了相应的数值实验，得到如下结论：

1)直接采用交错网格有限差分所得波场快照进行波场分离，无法彻底分离qP波与qS波； 尽管通过对模拟结果进行简单的算术平均插值处理后再进行波场分离，能够在一定程度上提高波场分离的效果，但仍不能彻底解决波场无法完全分离的问题．

2)采用旋转交错网格有限差分所得波场快照进行波场分离，能够得到分离较为彻底的qP波和qS波，且无需进行任何插值处理．

因此，当对各向异性介质进行波场数值模拟并进行相应的波场分离时，旋转交错网格有限差分是一种较为理想的选择．

评审专家和编辑部对本文提出了宝贵意见和建议，中国石油大学(华东)EWEP课题组提供了帮助，作者在此表示衷心的感谢.