CN114814944B - 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 - Google Patents

Abstract

本发明公开了属于地震波信息处理技术领域，具体公开了一种基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法，包括：步骤一、根据地震波正演模拟获取总矢量地震波场U；步骤二、求解经振幅相位校正以后的Helmholtz势纵波P^cor和横波S^cor；以及求解矢量纵波U^P和矢量横波U^S；还包括步骤三、使用标量成像条件和基于振幅和符号的矢量成像条件来分别生成PP反射率图像和PS反射率图像。本发明解决了现有技术中的逆时偏移成像叠加剖面无法反映真实的地下反射率信息的技术问题。

Description

基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法

技术领域

本发明涉及地震波信息处理技术领域，具体涉及一种基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法。

背景技术

当今油气勘探正在由构造油气藏向岩性和隐蔽油气藏、由常规油气开发向页岩气和致密油等非常规油气开发转变，在此趋势下，多波地震勘探的优势意义重大。多波地震勘探技术能够为复杂储层油气藏勘探提供有利的解决方案，其优势在于有效的转换波成像。相比于常规纵波勘探可以获得更多的地下地质信息。多波勘探可以有效减小地球物理反演的多解性，不仅可以改善构造成像，在岩性描述、流体预测、裂缝检测、各向异性检测等油气藏勘探开发关键问题上，也能够提供更佳的解决方案。其中，弹性逆时偏移(elasticreverse-time migration，ERTM)是专门为多分量地震数据而发展的偏移成像技术。弹性逆时偏移直接以对多分量数据进行偏移处理，使用弹性波动方程进行波场延拓，能够保持波场的弹性特征和矢量特征，更为真实地模拟地震波场在地下介质中的传播过程；弹性逆时偏移可同时得到PP、PS、SP和SS等多种成像结果，为地震解释提供更可靠的依据。而充分体现这些特点和优势的关键在于提供有效的转换波成像。

波场分离是进行弹性逆时偏移的前提条件，无论哪种纵横波分离方法，最终弹性波场会分离为标量波场或者矢量波场。针对标量波场，ERTM往往使用标量成像条件。只要标量P波含有正确的相位振幅信息，那么PP标量成像条件就可以近似真反射率信息。但是 PS标量成像条件存在两个问题，一个是3D情况下经Helmholtz分解得到的S波是矢量，往往需要额外的转换计算才能获得一个标量成像结果；另一个是S波在法向入射两侧存在极性相反问题，致使在多炮叠加时转换波成像剖面中的反射同相轴因极性不同相互抵消而受到破坏，最终叠加剖面无法反映转换波的反射率信息，因此需要引入额外的极性校正。极性校正是个复杂费时的过程，往往也需要引入诸如Poynting矢量、界面信息等其他参数，如果这些参数求取不准确，也会影响最终的叠加剖面。

针对矢量波场，ERTM往往采用点积矢量成像条件。点积矢量成像条件直接利用震源和检波点矢量波场点乘相加，计算简便但是成像结果是反射率与入射和反射角的余弦或正弦的综合效应，无法反映真反射率结果。后来基于振幅符号的成像条件将成像条件分成两部分，振幅部分由入射波和反射波振幅的乘积来确定，符号则通过入射波和反射波的矢量内积的符号来确定。这种矢量成像条件对于PP成像也存在一定问题，在入射和反射夹角90°以外的区域仍旧存在极性反转问题。

发明内容

本发明目的在于提供一种基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法，以克服现有技术中的逆时偏移成像叠加剖面无法反映真实的地下反射率信息的技术问题。

本发明所提供的基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法，包括：

步骤一、根据地震波正演模拟获取总矢量地震波场U；

步骤二、根据以下关系求解经振幅相位校正以后的Helmholtz势纵波P^cor和横波S^cor：

其中，α和β分别为纵波和横波速度；

根据以下关系求解矢量纵波U^P和矢量横波U^S：

步骤三、使用如下的标量成像条件和基于振幅和符号的矢量成像条件来分别生成PP 反射率图像和PS反射率图像：

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

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

进一步的，利用有限差分算法计算离散化的震源侧及检波点侧的总矢量地震波场

u、w分别为x和z方向的质点位移，并计算震源侧离散化的Helmholtz势纵波

以及检波点侧的离散化的Helmholtz势纵波

和横波

如下：

其中，i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数，Δt为所述有限差分算法在时间上的采样间隔，Δx为所述有限差分算法在空间上的采样间隔。

进一步的，通过以下的时间迭代计算获得离散化的震源侧矢量纵波

进一步的，通过以下的时间迭代计算获得离散化的检波点侧矢量横波

进一步的，离散化的矢量地震波场

依照如下的时间迭代计算得到:

本发明的原理和有益效果在于，构建了的基于散度和旋度的弹性波场分离技术，既可以得到标量波场，也可以得到矢量波场，且分离的波场与原始波场保持一致的振幅相位信息，基于此波场分离方法的弹性逆时偏移技术可适用于二维和三维情况，利用标量和矢量成像条件分别产生PP和PS成像剖面，由于校正后的Helmholtz势P波保持了振幅和极性信息，在利用标量成像条件生成PP成像时，在入射和反射夹角90°以外的区域也不会不存极性反转问题，另外矢量P和S波在利用基于振幅符号的矢量成像条件生成PS成像时不会受入射反射角影响，因此本发明在多炮叠加以后不会对叠加剖面造成破坏，也无需要引入额外的极性校正，叠加剖面能够真实的反映地下反射率信息。

附图说明

图1为本发明实施例中的基于散度和旋度的弹性波场分离和逆时偏移成像方法的流程图。

图2为本发明实施例中的第一个模拟实验所采用的模型的示意图。

图3为本发明实施例中的第一个模拟实验中0.7s时刻的波场快照示意图，其中，a)和b)子图分别是原始总矢量波场的水平方向位移u和垂直方向位移w，c)和d)子图分别是振幅和相位校正后的Helmholtz势纵波波场P^cor和横波波场S^cor，e)～h)子图分别是矢量纵波和横波的纵向和横向分量。

图4为本发明实施例中的第一个模拟实验的逆时偏移成像结果示意图，其中a)和b)子图是采用标量成像条件，c)和d)子图是采用点积矢量成像条件，e)和f)子图是采用基于振幅符号的矢量成像条件，a)、c)和e)子图是PP单炮成像结果，b)、d)和f)子图是PS单炮成像结果。

图5为本发明实施例中的第二个模拟实验的模型参数以及相应的初始平滑模型的参数示意图，其中，子图a)和b)分别是Marmousi2模型的P波速度和密度，子图c)和d)分别为初始平滑模型的P波速度和密度。

图6为本发明实施例中的第二个模拟实验中Marmousi2模型的PP和PS逆时偏移叠加剖面图。

具体实施方式

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

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

本实施例中的基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法包括:

根据地震记录获取总矢量地震波场U；

于是，本实施例中首先基于二维弹性波动方程在时间空间域进行有限差分正演模拟，计算出矢量地震波场U＝(u,w)。

以二阶位移方程为例，有：

式中，u、w分别为x和z方向的质点位移，α、β分别为地下介质P波和S波速度。

于是，根据有限差分数值计算原理进行离散化，得到u、w分别的时间递推方程：

式中，Δx和Δt分别是有限差分算法在空间和时间上的采样间隔，i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数。

接下来，本实施例分别对震源侧波场和检波点侧波场进行逆时延拓，并进行纵横波场分离。

在纵横波场的分离中，首先需要根据以下关系求解经振幅相位校正以后的Helmholtz 势纵波P^cor和横波S^cor：

其中，α和β分别为纵波和横波速度；

而后根据以下关系求解矢量纵波U^P和矢量横波U^S：

本实施例中，震源侧和检波点侧波场的逆时延拓需要同时进行，其中，震源侧波场逆时延拓时需要将上一步保存的边界波场值作为边界条件，最后时刻的波场值作为初始条件，进行逆时外推；检波点侧波场逆时延拓时需要将去除直达波以后的地震记录作为边界条件，进行逆时外推。震源侧和检波点侧波场均利用前述的u、w分别的时间迭代公式进行逆时外推，然后，利用前述关系在逆时外推过程中均需要进行纵横波场分离。

于是，综合上述计算关系，用于计算出震源侧波场标量纵波

和矢量纵波

的递推方程分别为：

用于计算出检波点侧波场标量纵波

和矢量横波

且矢量横波的计算需要先求取中间量Helmholtz势横波

因此检波点侧波场的递推方程为：

第三步，执行逆时偏移成像条件得到偏移剖面，本发明采用以下的标量和矢量成像条件分别生成PP和PS反射率图：

其中，符号函数的定义为：

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

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

s1.设定有限差分算法相关参数，包括空间采样间隔Δx、时间采样间隔Δt、水平和垂直网格点数nx和nz，模拟时间长度，读取速度模型包括模型大小、纵波速度α和横波速度β和密度ρ，并将模型参数网格离散化，设定炮点个数和炮点位置，设定检波点个数和位置，设定输出参数等。

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

s3.加载震源，根据震源类型在震源位置(nsx,nxz)处加载震源项f_s，本实施例中，设定震源为爆炸源，有：

u_nsx,nxz＝u_nsx,nxz+f_s

w_nsx,nxz＝w_nsx,nxz+f_s。

s4.波场迭代：根据波场有限差分时间迭代公式，利用前一时刻波场值

当前时刻波场值

纵波速度α、横波速度β、密度ρ以及有限差分空间和时间采样间隔Δx、Δt，计算得到下一时刻的波场值

同时为防止模型四周出现强反射，需要加载PML吸收边界。

s5.存储边界波场值。

s6.判断是否达到时间循环总次数nt，未达到则返回s3，否则时间循环终止。

s7.存储最后两个时刻的波场值。

s8.读取已经存储的震源最后两个时刻的波场值作为震源反传的初始条件。

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

s10.震源侧波场逆时延拓，读取已经存储的震源侧波场边界波场值，并利用后时刻震源侧波场值

当前时刻震源波场值

纵波速度α、横波速度β、密度ρ以及有限差分空间和时间采样间隔Δx、Δt，计算得到前时刻的震源侧波场值

s10.震源波场分离。需要计算出震源波场的标量纵波

和矢量纵波

s11.检波点波场逆时延拓。读取检波点记录的波场值，并利用后时刻检波点波场值、当前时刻检波点波场值以及模型和有限差分相关参数，计算得到前时刻的检波点波场值。

s12.检波点波场分离。需要计算出标量纵波

和矢量横波

s13.执行震源归一化的互相关成像条件。利用

和

直接进行互相关计算，

和

进行振幅互相关和符号计算，并且将每一时刻得到的互相关成像图进行叠加。

s14.判断是否达到有限差分时间循环总次数nt，未达到则返回s11，否则时间循环终止。

s15.时间循环结束得到此炮点上的PP和PS成像剖面。

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

s17.判断是否达到总炮点个数nshot，未达到则返回s2，否则炮点循环终止。

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

s19.输出最终的PP和PS成像剖面。

本实施例中的第一个模拟实验采用以下相关参数设定：

1.模型参数：如图2所示，模型大小为5600×1600m²，采用了双层模型，其中第一层P波和S波速度分别是α＝2000m/s、β＝1154m/s，第二层P/S波速分别为α＝2500m/s、β＝1443m/s。

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

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

图3为本次模拟实验中0.7s时刻的波场快照，展示了纵横波分离的结果，其中，a)和b)子图分别是原始总矢量波场的水平方向位移u和垂直方向位移w，c)和d)子图分别是振幅和相位校正后的Helmholtz势纵波波场P^cor和横波波场S^cor，e)～h)子图分别是矢量纵波和横波的纵向和横向分量。

本次模拟实验的逆时偏移成像结果如图4所示，其中a)和b)子图是采用标量成像条件， c)和d)子图是采用点积矢量成像条件，e)和f)子图是采用基于振幅符号的矢量成像条件， a)、c)和e)子图是PP单炮成像结果，b)、d)和f)子图是PS单炮成像结果。从c)和e)子图中可看出，对于PP成像，两种矢量成像结果对于大角度反射波存在极性反转问题，只有标量成像结果(子图a))能够正确指示真实反射率信息。而对于PS成像，标量成像结果(子图b))由于在垂直入射两侧也存在极性反转问题，只有基于振幅符号的矢量成像结果(子图f))才能够指示真实反射率信息。

本实施例中的第二个模拟实验采用以下相关参数设定：

1.模型参数：采用Marmousi2模型，模型参数以及相应的初始平滑模型的参数如图5 所示，其中，Marmousi2模型的P波速度和密度如子图a)和b)所示，初始平滑模型的P波速度和密度如子图c)和d)所述；S波速度设定为

模型大小为4000×2240m²。

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

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

本次模拟实验的逆时偏移成像结果如图6所示，其中a)和b)子图是采用标量成像条件，c)和d)子图是采用点积矢量成像条件，e)和f)子图是采用基于振幅符号的矢量成像条件， a)、c)和e)子图是PP单炮成像结果，b)、d)和f)子图是PS单炮成像结果。与图4中的真实模型进行对比，从图6中可看出PP标量成像结果(子图a)和PS基于振幅符号的矢量成像结果(子图f)基本能够正确恢复层位，尤其针对模型细节的成像更清晰，因此可用来正确指示地下真实反射率信息。

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

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

Claims (6)

1.基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法，包括：

步骤一、根据地震波正演模拟获取总矢量地震波场U；

步骤二、根据以下关系求解经振幅相位校正以后的Helmholtz势纵波P^cor和横波S^cor：

其中，α和β分别为纵波和横波速度；

根据以下关系求解矢量纵波U^P和矢量横波U^S：

其特征在于，

还包括步骤三、使用如下的标量成像条件和基于振幅和符号的矢量成像条件来分别生成PP反射率图像和PS反射率图像：

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

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

3.根据权利要求1所述的方法，其特征在于，利用有限差分算法计算离散化的震源侧及检波点侧的总矢量地震波场

u、w分别为x和z方向的质点位移，并计算震源侧、检波点侧的离散化的Helmholtz势纵波

以及检波点侧的离散化的Helmholtz势横波

如下：

其中，i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数，Δt为所述有限差分算法在时间上的采样间隔，Δx为所述有限差分算法在空间上的采样间隔。

4.根据权利要求3所述的方法，其特征在于，通过以下的时间迭代计算获得离散化的震源侧矢量纵波

5.根据权利要求3所述的方法，其特征在于，通过以下的时间迭代计算获得离散化的检波点侧矢量横波

6.根据权利要求3-5中任一所述的方法，其特征在于，离散化的矢量地震波场

依照如下的时间迭代计算得到:

Images