CN117233838B - 一种二维vti介质中的弹性准纵横波场分离和逆时偏移成像方法 - Google Patents

Abstract

本发明涉及地震波资料处理技术领域，具体涉及一种二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法。本发明所提供的二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法，包括：步获取二维VTI介质中的速度波场以及正应力波场；利用VTI介质中的空间域分离算子求得经振幅相位校正以后的Helmholtz势qP标量波场P^cor以及qP和qSV矢量波场；利用前述步骤获得的震源侧和检波点侧分别的qP标量波场以及qP和qSV矢量波场，并生成PP反射性图像和PS反射性图像。该方法实现了在时间空间域进行各向异性波场分离，降低了成像计算量且有利于模型的区域分解和并行加速计算。

Description

一种二维VTI介质中的弹性准纵横波场分离和逆时偏移成像 方法

技术领域

本发明涉及地震波资料处理技术领域，具体涉及一种二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法。

背景技术

深度域地震资料成像是地震资料处理的关键步骤，弹性波逆时偏移技术是当前业界重点发展的多分量地震资料深度域成像技术之一。其中VTI(vertical transverseisotropy)介质模型是最具代表性的地层各向异性模型。

早期的弹性波逆时偏移成像结果是纵波和横波耦合在一起的，其成像结果物理含义不明确，并且不同波场之间存在串扰噪音，影响成像质量。基于纵横波场分离的弹性逆时偏移技术能够合理地利用多波地震数据，得到具有明确物理意义的成像结果，并且能减弱串扰噪音，提高成像质量。同样地，准纵横波场分离是进行各向异性弹性逆时偏移的前提条件。目前波场分离方法均基于准纵、横波(qP、qSV波)的偏振方向，但是由于qP波的偏振方向与传播方向不平行，qSV波的偏振方向与传播方向也不垂直，他们的偏振方向与传播方向的夹角随空间变化，因此，各向同性介质中通过对弹性波场求取散度和旋度的方式无法在各向异性介质中适用。

各向异性介质波场分离算子可由qP或者qSV的偏振向量构建，而偏振向量可通过求解Christoffel方程得到。由于VTI介质的波场分离算子是传播方向的函数，是一个空间波数域混合算子。最直接的方法在波数域中进行波场分离，理论上只能适用于均匀介质，较难扩展到非均匀介质。另外一种分离方法是采用低秩近似方法，但是低秩近似方法需要密集的FFT操作，计算量大，且不利于模型的区域分解和并行加速计算。

另一发面，目前针对VTI介质的弹性逆时偏移成像条件，使用最多的是矢量点积成像条件。矢量点积成像条件直接利用震源和检波点矢量波场点乘相加，计算简便但是成像结果是反射率与入射和反射角的余弦或正弦的综合效应，PP成像会在大偏移距的地方产生极性反转问题，多炮叠加以后不能正确成像，PS成像会在近垂直入射的地方产生较弱的成像剖面，无论是PP还是PS成像均无法反映真反射率结果。

发明内容

本发明目的在于提供一种二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法，该方法能够在时间空间域进行各向异性波场分离。

本发明所提供的二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法，包括：

步骤一，获取二维VTI介质中的速度波场V＝(V_x,V_z)以及正应力波场σ_xx、σ_zz；

步骤二、根据以下关系求得经振幅相位校正以后的Helmholtz势qP标量波场P^cor：

根据以下关系求得qP和qSV矢量波场V_P和V_S：

式中，为VTI介质中的空间域分离算子，表示如下：

和/>是空间域的标量算子，其表达式为：

其中，n＝(n_x,n_z)＝(s_x/|s|,s_z/|s|)为传播方向单位矢量，为坡印廷向量；c_ij为VTI介质的弹性矩阵中的系数，其中角标i、j分别对应系数在弹性矩阵中的行和列；

步骤三、利用前述步骤获得的震源侧和检波点侧分别的qP标量波场以及qP和qSV矢量波场，并生成PP反射性图像和PS反射性图像。

进一步的，步骤二中，以如下简化的方式计算P^cor：

进一步的，步骤三中，使用如下的点积成像条件来生成PP反射性图像和PS反射性图像：

其中，I^PP和I^PS是PP和PS成像剖面，角标src和rev分别代表震源侧和检波点侧，||·||运算是求矢量的模，符号函数sgn^PS(x,t)通过求解震源侧和检波点侧矢量波场的内积来确定如下：

进一步的，还包括，对于每一炮地震记录均采用上述三步进行单炮偏移成像，多炮叠加以后生成最终的成像剖面。

进一步的，步骤1中，在震源处加载震源子波S(t)和检波点记录R(t)，基于VTI介质中的一阶速度-应力波动方程在时间空间域进行采用交错网格的有限差分正演迭代和逆时延拓,计算出震源侧和检波点侧的速度波场和应力波场。

进一步的，步骤1中，还对原始的震源子波S(t)和检波点记录R(t)分别进行-1/ω²滤波，得到滤波后的震源子波S′(t)和检波点记录R′(t)如下：

式中和/>为正反傅里叶变换；

以及，加载滤波后的震源子波S′(t)和检波点记录R′(t)，利用交错网格有限差分正演和逆时延拓，得到震源波场V_s′_r′_c、检波点波场V_r′_e′_c以及相应的应力波场；

在波场逆时延拓的同时利用时间迭代计算进行qP和qSV矢量波场分离。

进一步的，在有限差分正演迭代时，加载PML吸收边界。

进一步的，在有限差分正演迭代时，保存边界处M层的速度和应力波场值，以及最后时刻的速度和应力波场值；

震源和检波点波场的逆时延拓同时进行，其中，震源波场逆时延拓时将保存的边界波场值作为边界条件，最后时刻的波场值作为初始条件。

进一步的，检波点波场逆时延拓时，将去除直达波以后的地震记录R(t)作为边界条件。

进一步的，VTI介质中的一阶速度-应力波动方程表示如下：

式中，V_x、V_z分别为x和z方向的质点速度，σ_xx、σ_zz和σ_xz分别为应力张量中xx方向，zz方向和xz方向的元素，ρ为VTI介质的密度。

本发明的原理和有益效果在于，提出了能够在时间空间域进行各向异性波场分离的方法，该方法无需密集的FFT操作，计算成本低且方便进行GPU加速。

本发明的一些实施例中，采用标量和矢量结合的成像条件，其中PP采用标量成像，而PS采用基于振幅和符号的矢量成像，两者均能得到正确的成像剖面。

附图说明

图1为本发明实施例中的三层VTI模型及0.44s时刻的波场快照图；

图2为本发明实施例中的三层VTI介质中单炮PP和PS逆时偏移成像结果图以及与传统的矢量点积成像条件对比图；

图3为本发明实施例中的复杂VTI模型(Marmousi2模型)的Thomson参数和密度分布图；

图4为本发明实施例中的复杂VTI模型的初始平滑模型及0.88s时刻qP和qSV分离的波场快照图；

图5为本发明实施例中的复杂VTI模型的PP和PS逆时偏移成像叠加剖面图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例中，利用基于波场模型和VTI介质模型的模拟实验示例性的说明本发明的技术方案及技术效果。

本实施例中的逆时偏移成像方法包括:

获取二维VTI介质中的原始混合速度波场V＝(V_x,V_z)以及xx和zz方向的正应力波场σ_xx、σ_zz；

然后，通过应力波场求qP标量波场P^cor用于PP成像，

并进行qP和qSV波的矢量波场分离，求得qP和qSV波的矢量波场V_P和V_S，用于PS成像；

为达到上述目的，核心是建立能够在空间域快速计算的分离算子，往往由qP波或qSV波的偏振矢量产生。通过求解二维VTI介质的Christoffel方程的特征向量可得到qP波和qSV波的偏振矢量在波数域的表达形式分别为：

其中，c_ij是VTI介质的弹性矩阵系数，k_x和k_z是波数矢量的x和z分量。由于/>本发明选择两者之一的qP波的偏振矢量建立分离算子/>但是分离算子无法直接转到空间域进行波场分离，往往只能在波数域进行计算，使其无法有效扩展到非均匀介质。现有技术中，采用低秩近似方法进行波场分离，但是低秩近似的方法要求每个时间迭代步内进行频繁地FFT操作，计算量大且不易进行区域分解和并行加速计算，由此限制了其实际应用。为此本发明借鉴标量算子的思想，令

其中S_k为引入的波数域标量算子：

通过波数矢量与传播方向单位矢量/>的关系/>波数域标量算子S_k可转化为空间域标量算子：

式中，

空间域标量算子S_n的计算则转化为传播方向单位矢量的计算。坡印廷矢量是求解局部传播方向最快速和低成本的一种方式，理论上需要精确地qP波传播方向，但是实际上无法在波分离之前得到分离波的传播方向，因此本例中首先得到一个粗略的分离qP波场，进行传播方向估计，以便用于后续的波场分离。而这个粗略的分离qP波场即是通过各向同性分离算子(散度)得到的分离波场，等价于应力波场(σ_xx+σ_zz)，而后对这个粗略的分离波场求取坡印廷矢量即可，同时由于我们求取的是归一化单位向量，所以可以忽略坡印廷矢量中的时间导数项，即可得到传播方向单位矢量的计算公式：

其中，可以看出标量算子S_n完全可以在空间域完成计算。

那么，将VTI介质中的波数域分离算子转换到空间域，得到VTI介质的空间域分离算子为：

可以看出，当在各向同性介质中时，即令ε＝δ＝0，此时S_n＝S_k＝1，分离算子式(6)则转化为各向同性介质中的波场分离算子：

借鉴各向同性介质中的分离思想，可得到VTI介质中的波场分离理论。假设原始混合速度波场为V＝(V_x,V_z)，那么Helmholtz势qP和qSV波场分别记为P和S：

经振幅和相位校正以后得Helmholtz势qP和qSV波场记为P^cor和S^cor：

其中P^cor即为标量qP波场，可直接用于后期的PP逆时成像。矢量qP和qSV波场分别为V_P和V_S：

式中和/>也是空间域的标量算子，其表达式为：

通过数值实验发现，由于标量算子S_n、是通过计算粗略的分离qP波场的传播方向估计的，其计算具有不精确性，会导致在迭代过程中有较大的累积误差，因此上式是不适合进行时间迭代计算的。

考虑到PP成像只需要标量qP波场P^cor，PS成像需要矢量qP和qSV波场V_P和V_S。针对P^cor的计算可以走一条捷径，直接通过应力求取P^cor，因为他们之间存在以下关系：

针对矢量qP和qSV波的求取，式(10)中对分离的矢量波场V_P和V_S的二阶求导可以转化为对混合波场V的二次积分，因此本实施例中，先对震源子波或检波点记录做-1/ω²滤波，然后进行正演迭代和逆时延拓得到的混合波场为V″，最后利用下式：

可以得到分离的矢量qP波场V^P和矢量qSV波场V^S，它们与原始混合波场保持一致的振幅相位信息。

在此基础上，利用震源和检波点各自的qP标量波场P^cor，以及qP和qSV矢量波场V_P和V_S，可进行后续的弹性逆时偏移成像，后续使用角标src和rec区分震源和检波点波场。

在各向异性弹性逆时偏移过程中，往往使用矢量点积成像条件。但是这种成像结果是反射率与入射和反射角的余弦或正弦的综合效应，不能反映真实的反射率结果，尤其是单炮PP成像会在大偏移距的地方存在极性反转问题，多炮叠加以后不能正确成像。针对以上情况，本发明使用标量和矢量结合的成像条件，其中PP采用标量成像，而PS采用基于振幅和符号的矢量成像条件，即：

式中，I^PP和I^PS是PP和PS成像剖面，角标src和rec分别代表震源和检波点波场，||·||运算是求矢量的模，符号函数通过求解震源和检波点矢量波场的内积来确定：

于是，更为具体的，本实施例中首先在震源处加载震源子波S(t)，基于VTI介质中的一阶速度-应力波动方程在时间空间域进行采用交错网格的有限差分正演和逆时延拓；计算出速度波场V＝(V_x,V_z)和应力波场(σ_xx+σ_zz)。VTI介质中的一阶速度-应力波动方程表示如下：

式中，V_x、V_z分别为x和z方向的质点速度，σ_xx、σ_zz和σ_xz分别为应力张量中xx方向，zz方向和xz方向的元素，ρ为VTI介质的密度。

根据一阶速度-应力波动方程的空间2M阶和时间2阶交错网格有限差分迭代公式，速度和应力分量不仅需要交错的分布在空间网格点上，而且也是交错分布在时间节点上。因此，采用交错网格有限差分方法求解VTI介质中的一阶速度-应力波动方程式(15)，令Δt为有限差分算法在时间上的采样间隔，首先利用t时刻应力波场更新t+1/2Δt时刻速度波场，然后利用t+1/2Δt时刻混合速度波场更新t+Δt应力波场。同时为防止模型四周出现强边界反射，需要加载PML吸收边界。

同理，在加载检波点记录R(t)的基础上，可获得正传的检波点速度波场和应力波场，

存储边界处M层的速度和应力波场值、存储最后时刻的速度和应力波场值。

震源和检波点波场的逆时延拓需要同时进行，其中，震源波场逆时延拓时将保存的边界波场值作为边界条件，最后时刻的波场值作为初始条件，进行逆时延拓。检波点波场逆时延拓时将去除直达波以后的地震记录R(t)作为边界条件，进行逆时延拓。震源和检波点波场具有相同的逆时延拓波动方程式(15)，采用交错网格有限差分方法求解VTI介质中的一阶速度-应力波动方程式(15)，利用t时刻应力波场更新t-1/2Δt时刻速度波场，然后利用t-1/2Δt时刻速度波场更新t-Δt应力波场。最后，分别得到震源和检波点的应力波场。

利用反传的应力波场，值依照公式(12)计算出震源波场和检波点波场的标量qP波场和/>用于PP成像。

另一方面，对原始震源子波S(t)和检波点记录R(t)分别进行-1/ω²滤波，得到滤波后的震源子波S′(t)和检波点记录R′(t)：

式中和/>为正反傅里叶变换。

加载滤波后的震源子波S′(t)和检波点记录R′(t)，并重复上述的交错网格有限差分正演模拟和逆时延拓，得到震源波场V″_src和检波点波场V″_rec；由于进行了滤波处理克服了qP和qSV波不适合进行时间迭代计算的问题，于是，在波场逆时延拓的同时可同时利用时间迭代计算依照公式(11)和(13)进行qP和qSV矢量波场分离。

本实施例中，还对于每一炮地震记录均采用上述步骤进行单炮偏移成像，多炮叠加以后生成最终的成像剖面。

本实施例中的模拟实验的具体实施流程如下：

一、PP成像实施流程

1.读取VTI介质模型参数，包括模型大小、各向异性弹性参数(v_P0、v_S0、ε、δ)以及密度ρ。设定有限差分算法相关参数，包括空间采样间隔Δx、时间采样间隔Δt、水平和垂直网格点数nx和nz、模拟时间长度，并将模型参数网格离散化。设定炮点个数和炮点位置，设定检波点个数和位置，设定输出参数等。

2.进入炮循环，若小于或等于总炮点个数nshot，进行下一步。

3.震源波场正向延拓。读取震源参数包括震源类型、震源位置(nsx,nxz)、震源子波S(t)等参数，

3.1进入时间循环，若小于或等于时间循环总次数nt，进行下一步。

3.2加载震源。根据震源类型在震源位置(nsx,nxz)处加载t时刻震源幅值f_s＝S(t)，若为爆炸源，则在正应力分量上加载σ_xx|_nsx,nxz＝σ_xx|_nsx,nxz+f_s，σ_zz|_nsx,nxz＝σ_zz|_nsx,nxz+f_s。

3.3波场迭代。根据一阶速度-应力波动方程的空间2M阶和时间2阶交错网格有限差分迭代公式，速度和应力分量不仅需要交错的分布在空间网格点上，而且也是交错分布在时间节点上。因此首先利用t时刻应力波场更新t+1/2Δt时刻速度波场，然后利用t+1/2Δt时刻速度波场更新t+Δt应力波场。同时为防止模型四周出现强边界反射，需要加载PML吸收边界。

3.4存储边界处M层的速度和应力波场值

3.5判断是否达到时间循环总次数nt，未达到则返回3.1步，否则时间循环终止。

3.6存储最后时刻的速度和应力波场值

4.震源和检波点波场的逆时延拓

4.1读取已经存储的震源正传波场最后时刻的速度和应力波场值作为震源反传的初始条件。

4.2进入逆时延拓时间循环，若小于或等于时间循环总次数nt，进行下一步。

4.3震源波场逆时延拓。读取已经存储的t时刻震源正传波场边界处M层的速度和应力波场值，利用t时刻应力波场更新t-1/2Δt时刻速度波场，然后利用t-1/2Δt时刻速度波场更新t-Δt应力波场。

4.4震源波场分离。利用震源反传的应力波场值计算出震源波场的标量qP波场

4.5检波点波场逆时延拓。读取t时刻检波点记录的速度波场值，利用t时刻应力波场更新t-1/2Δt时刻速度波场，然后利用t-1/2Δt时刻速度波场更新t-Δt应力波场。

4.6检波点波场分离。利用检波点反传的应力波场值计算出检波点波场的标量qP波场

4.7执行震源归一化的标量互相关成像条件。利用标量波场和/>直接进行互相关计算，并且将每一时刻得到的互相关成像结果进行叠加。

4.8判断是否达到有限差分时间循环总次数nt，未达到则返回第4.2步，否则时间循环终止。

4.9时间循环结束得到此炮点上的PP成像剖面。

5.将每一炮得到的PP成像剖面进行叠加。

6.判断是否达到总炮点个数nshot，未达到则返回第2步，否则炮点循环终止。

7.得到最终的PP叠加剖面，进行Laplace滤波，消除低波数噪音。

8.输出最终的PP成像剖面。

二、PS成像实施流程：

1.读取VTI介质模型参数，包括模型大小、各向异性弹性参数(v_P0、v_S0、ε、δ)以及密度ρ。设定有限差分算法相关参数，包括空间采样间隔Δx、时间采样间隔Δt、水平和垂直网格点数nx和nz、模拟时间长度，并将模型参数网格离散化。设定炮点个数和炮点位置，设定检波点个数和位置，设定输出参数等。

2.读取震源子波S(t)和检波点记录R(x,t)，对两者做-1/ω²滤波得到滤波后的震源子波S′(t)和检波点记录R′(x,t)。

3.进入炮循环，若小于或等于总炮点个数nshot，进行下一步。

4.震源波场正向延拓。

4.1进入时间循环，若小于或等于时间循环总次数nt，进行下一步。

4.2加载震源。根据震源类型在震源位置(nsx,nxz)处加载t时刻震源幅值f_s＝S′(t)，若为爆炸源，则在正应力分量上加载σ_xx|_nsx,nxz＝σ_xx|_nsx,nxz+f_s，σ_zz|_nsx,nxz＝σ_zz|_nsx,nxz+f_s。

4.3波场迭代。根据一阶速度-应力波动方程的空间2M阶和时间2阶交错网格有限差分迭代公式，速度和应力分量不仅需要交错的分布在空间网格点上，而且也是交错分布在时间节点上。因此首先利用t时刻应力波场更新t+1/2Δt时刻速度波场，然后利用t+1/2Δt时刻速度波场更新t+Δt应力波场。同时为防止模型四周出现强边界反射，需要加载PML吸收边界。

4.4存储边界处M层的速度和应力波场值

4.5判断是否达到时间循环总次数nt，未达到则返回4.1步，否则时间循环终止。

4.6存储最后时刻的速度和应力波场值

5.震源和检波点波场的逆时延拓

5.1读取已经存储的震源正传波场最后时刻的速度和应力波场值作为震源反传的初始条件。

5.2进入逆时延拓时间循环，若小于或等于时间循环总次数nt，进行下一步。

5.3震源波场逆时延拓。读取已经存储的t时刻震源正传波场边界处M层的速度和应力波场值，利用t时刻应力波场更新t-1/2Δt时刻速度波场，然后利用t-1/2Δt时刻速度波场更新t-Δt应力波场。

5.4震源波场分离。利用震源反传的速度波场值计算出震源波场的矢量qP波场和qSV波场/>

5.5检波点波场逆时延拓。读取t时刻检波点记录的速度波场值，利用t时刻应力波场更新t-1/2Δt时刻速度波场，然后利用t-1/2Δt时刻速度波场更新t-Δt应力波场。

5.6检波点波场分离。利用检波点反传的应力波场值计算出检波点波场的矢量qP波场和qSV波场/>

5.7执行震源归一化的基于振幅和符号的互相关成像条件。利用震源矢量qP和qSV波场(和/>)与检波点矢量波场矢量qP和qSV波场(/>和/>)进行基于振幅和符号的互相关计算，并且将每一时刻得到的互相关成像结果进行叠加。

5.8判断是否达到有限差分时间循环总次数nt，未达到则返回第5.2步，否则时间循环终止。

5.9时间循环结束得到此炮点上的PS成像剖面。

6.将每一炮得到的PS成像剖面进行叠加。

7.判断是否达到总炮点个数nshot，未达到则返回第3步，否则炮点循环终止。

8.输出最终的PS成像叠加剖面。

本例中，数值模拟实验分为两个，一个为简单层状模型，一个为复杂的Marmousi2模型。

第一个模拟实验采用以下相关参数：

1.模型参数：模型大小为3200×1200m²，三层模型(图1a)，其中第一层Thomson各向异性参数和密度分别为V_P0＝2600m/s,V_S0＝1400m/s,ε＝0,δ＝0,ρ＝1.8kg/m³，第二层参数为V_P0＝2700m/s,V_S0＝1500m/s,ε＝0.15,δ＝-0.1,ρ＝2.0kg/m³，第三层参数为V_P0＝3000m/s,V_S0＝2000m/s,ε＝0.3,δ＝0.25,ρ＝2.2kg/m³。

2.震源参数：震源类型为爆炸源，震源时间函数为25Hz的Ricker子波，震源位于模型中心(1600m,40m)处。

3.有限差分数值算法参数：空间采样间隔Δx＝4m，时间采样间隔Δt＝0.4ms，网格点数801×301，总模拟时长为1.6s，时间迭代次数为nt＝4000。边界条件采用PML吸收边界，厚度为50层。

模型如图1a所示，纵横波分离结果如图1b～1h所示，图中纵轴为深度(Depth)为距离(Distance)，下同；逆时偏移成像结果如图2所示。其中，波场分离是0.44s时刻的波场快照，图1b和1c分别是原始总矢量波场水平(V_x)和垂直方向位移(V_z)，图1d是经振幅和相位校正后的Helmholtz势qP波场(P^cor)，图1e～1h分别是矢量qP波和qSV波的x和z方向分量。图2是单炮PP和PS逆时偏移成像图。

从图1的qP和qS波场分离快照可以看出，本发明的波场分离方法可以较好地将qP和qSV波分离开来，波场分离快照较干净，精度较高。从图2a和b中的逆时偏移成像图可以看出，与图2c和d中传统的矢量点积成像条件(PP/Psbydot-productimaging condition)相比，图2a中的标量PP成像(PPbyscalarimagingcondition)不存在大偏移距极性反转的问题(黑色箭头所指示的地方)，图2b中的基于振幅和符号的PS成像(PSbymagnitude-andsign-basedimagingcondition)在近垂直入射附近(黑色箭头所指示的地方)成像幅值较大，说明无论是PP还是PS成均能反映地下VTI介质的真实反射率信息。

第二个模拟实验采用以下相关参数：

1.模型参数：为各向异性Marmousi2模型，Thomson各向异性参数和密度分布如图3所示，qSV波速度设定为模型大小为4000×2240m²。

2.炮点和检波点参数：震源类型为爆炸源，震源时间函数为25Hz的Ricker子波，炮点深度位于40m，水平位置从40m到3964m总共激发99炮，跑间距为40m。

3.有限差分数值算法参数：空间采样间隔Δx＝4m，时间采样间隔Δt＝0.4ms，网格点数1001×561，总模拟时长为3.2s，时间迭代次数为nt＝8000。边界条件采用PML吸收边界，厚度为100层。

图3是模型参数，图4是初始平滑模型及波场分离快照，图5是各向异性Marmousi2模型的PP和PS逆时偏移成像叠加剖面。

从图4的各向异性Marmousi2模型qP和qS波场分离快照可以看出，本发明的波场分离方法在复杂VTI介质中也可以较好地将qP和qSV波分离开来，波场分离快照较干净，精度较高。图5的复杂VTI介质弹性逆时偏移成像剖面与真实VTI介质模型(图3)进行对比可以看出，无论是PP还是PS成均能对地下VTI介质的真实反射率进行成像。

在模拟计算中，各采样点的初始位移可设定为0，震源位置(nsx,nxz)则加载震源时间函数f_s，于是，矢量地震波场中各采样点在各采样时刻的波场值皆可计算得到。而在实际工作中，震源位置和震源时间函数可以通过对地震的观测或通过观测数据推算得到，进而模拟出总地震矢量波场，但不排除直接测量或推算总地震矢量波场的情况，无论何种方式获得的总地震矢量波场(包括其离散表达)，皆不影响本发明中的分解方法的实施，依旧属于本发明的保护范围。

本发明中未涉及部分均与现有技术相同或可采用现有技术加以实现。尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (10)

1.一种二维VTI介质中的弹性准纵横波场分离和逆时偏移成像方法，其特征在于，包括：

步骤一，获取二维VTI介质中的速度波场V＝(V_x,V_z)以及正应力波场σ_xx、σ_zz；

步骤二、根据以下关系求得经振幅相位校正以后的Helmholtz势qP标量波场P^cor：

根据以下关系求得qP和qSV矢量波场V_P和V_S：

式中，为VTI介质中的空间域分离算子，表示如下：

和/>是空间域的标量算子，其表达式为：

其中，n＝(n_x,n_z)＝(s_x/|s|,s_z/|s|)为传播方向单位矢量，为坡应廷向量；c_ij为VTI介质的弹性矩阵中的系数，其中角标i、j分别对应系数在弹性矩阵中的行和列；ρ为VTI介质的密度；

t为采样时间；

步骤三、利用前述步骤获得的震源侧和检波点侧分别的qP标量波场以及qP和qSV矢量波场，并生成PP反射性图像和PS反射性图像。

2.根据权利要求1所述的方法，其特征在于，步骤二中，以如下简化的方式计算P^cor：

3.根据权利要求1所述的方法，其特征在于，步骤三中，使用如下的点积成像条件来生成PP反射性图像和PS反射性图像：

其中，I^PP和I^PS是PP和PS成像剖面，角标src和rev分别代表震源侧和检波点侧，||·||运算是求矢量的模，符号函数sgn^PS(x，t)通过求解震源侧和检波点侧矢量波场的内积来确定如下：

4.根据权利要求1所述的方法，其特征在于，还包括，对于每一炮地震记录均采用步骤一至三进行单炮偏移成像，多炮叠加以后生成最终的成像剖面。

5.根据权利要求1所述的方法，其特征在于，步骤1中，在震源处加载震源子波S(t)和检波点记录R(t)，基于VTI介质中的一阶速度-应力波动方程在时间空间域进行采用交错网格的有限差分正演迭代和逆时延拓,计算出震源侧和检波点侧的速度波场和应力波场。

6.根据权利要求1所述的方法，其特征在于，步骤5中，还对原始的震源子波S(t)和检波点记录R(t)分别进行-1/ω²滤波，得到滤波后的震源子波S′(t)和检波点记录R′(t)如下：

式中F和F^-1为正反傅里叶变换；

以及，加载滤波后的震源子波S′(t)和检波点记录R′(t)，利用交错网格有限差分正演和逆时延拓，得到震源波场V″_src、检波点波场V″_rec以及相应的应力波场；

在波场逆时延拓的同时利用时间迭代计算进行qP和qSV矢量波场分离。

7.根据权利要求5或6所述的方法，其特征在于，在有限差分正演迭代时，加载PML吸收边界。

8.根据权利要求7所述的方法，其特征在于，在有限差分正演迭代时，保存边界处M层的速度和应力波场值，以及最后时刻的速度和应力波场值；

震源和检波点波场的逆时延拓同时进行，其中，震源波场逆时延拓时将保存的边界波场值作为边界条件，最后时刻的波场值作为初始条件。

9.根据权利要求8所述的方法，其特征在于，检波点波场逆时延拓时，将去除直达波以后的检波点记录R(t)作为边界条件。

10.根据权利要求5所述的方法，其特征在于，VTI介质中的一阶速度-应力波动方程表示如下：

式中，V_x、V_z分别为x和z方向的质点速度，σ_xx、σ_zz和σ_xz分别为应力张量中xx方向，zz方向和xz方向的元素。