CN107356972B - 一种各向异性介质的成像方法 - Google Patents

Abstract

本发明公开了一种各向异性介质的成像方法，包括如下步骤：S1，获得观测数据，根据观测数据获得速度场和各相异性参数；S2，对观测数据进行偏移得到反射系数；S3，根据速度场和各相异性参数，对反射系数进行反偏移得到模拟数据；S4，计算观测数据和模拟数据的数据残差；S5，根据数据残差获得最终成像结果。通过采用本发明提供的方法，充分考虑了地下介质的各向异性这一实际情况，能够更加精确的模拟地震波在真实地下介质中的传播规律，进一步提高偏移成像结果的分辨率和保真度，使得偏移成像结果中不同深度的同相轴振幅更加均衡，成像结果更为准确。

Description

一种各向异性介质的成像方法

技术领域

本发明涉及一种各向异性介质的成像方法，属于地球物理勘探领域。

背景技术

各向异性介质即物理性质具有方向特性的地球物理介质。通常，地下沉积地层大多成层分布，呈现横向各向同性(VTI，vertical transverse isotropic)介质特征。随着高精度地震勘探技术的发展，能够更加真实且全面的反映地震波在地下传播规律的各向异性介质成像方法越来越引人关注。当前的地下介质反演过程中，通常使用迭代算法，以获得成像结果。

但是，现有技术中的缺陷在于，通常的反演过程中，大多是基地下介质的各向同性假设，忽略介质的各向异性，导致模拟波场和真实波场在传播过程中出现偏差，降低成像时的准确性和精度。

发明内容

本发明的目的在于，提供一种各向异性介质的成像方法，它可以解决当前技术中存在的问题，根据地下介质的各向异性这一实际情况，获得更准确的反演成像结果。

为解决上述技术问题，本发明采用如下的技术方案：一种各向异性介质的成像方法，包括以下步骤：

S1，获得观测数据，根据观测数据获得速度场和各向异性参数；

S2，对观测数据进行偏移得到反射系数；

S3，根据速度场和各向异性参数，对反射系数进行的反偏移得到模拟数据；

S4，计算观测数据和模拟数据的数据残差；

S5，根据数据残差获得最终成像结果。

本发明提供的方法，充分考虑了地下介质的各向异性这一实际情况，能够更加精确的模拟地震波在真实地下介质中的传播规律，基于线性反演思想通过反复迭代进一步提高偏移成像结果的分辨率和保真度，使得偏移成像结果中不同深度的同相轴振幅更加均衡，成像结果更为准确。此外，本发明继承了传统单程波类偏移成像方法高精度、高效率的优势，同时占用较小的存储空间，偏移成像噪音少，对速度误差的敏感性更弱。

附图说明

图1为本发明的一个实施例流程示意图。

图2为本发明本发明的一个实施例中的速度场和各向异性参数图，其中

(a)为速度场，(b)为描述纵波水平速度和垂向速度关系的函数ε的图像，(c)为描述纵波各向异性强度δ；

图3a为对一个模型现有技术的成像结果，图3b为对前述模型本发明的成像结果。

下面结合附图和具体实施方式对本发明作进一步的说明。

具体实施方式

实施例1

一种各向异性介质的成像方法，如图1所示，包括以下步骤：

S1，获得观测数据，根据观测数据获得速度场和各向异性参数。

观测数据d^obs根据各向异性Hess模型获得，模型主要由左侧的高速体岩丘，中间的强各向异性异常体以及右侧的高陡断层组成。模型大小为901-376，纵横向网格间距均为10m。

观测数据共151炮，每炮401道接收，采用中间放炮两边接收的滚动观测系统，炮间隔60m，道间隔10m。记录时间3s，时间采样间隔3ms。根据模型得到的速度和各向异性参数正演使用雷克子波，主频20Hz。

所述的各向异性参数包括描述纵波水平速度和垂向速度关系的函数ε，和，描述纵波各向异性强度的参数δ，如图2所示，

S2，对观测数据进行偏移得到反射系数。

例如可以采用高阶傅里叶有限差分算子(FFD)偏移得到偏移成像结果，这里可以采用的方法很多，属于常规技术，就不再赘述。

S3，根据速度场和各向异性参数，对反射系数进行反偏移得到模拟数据。

具体的说，根据得到的反射系数m_k，根据速度场和各向异性参数进行正演模拟(也叫反偏移)得到模拟数据Lm_k。其中k表示迭代次数，L表示Born线性正演算子。

这里首先需要说明的是，地下波在反射时通常存在多次反射波，一次反射波是指地震波向下传播时遇到地层界面形成的反射波，一次反射波是相对多次波而言的，但是除了一次波以外的其他反射波通常能量较弱被当作噪音进行压制。因此，实际处理中所述的反射波通常指一次反射波。

频率域中的反射波U_s(x_r|x_s，ω)可以表示成以下公式

其中f(ω)是频率域震源子波，G(x'|x_s,ω)是VTI介质中由震源位置x_s传播到散射点x'的Green函数，G(x_r|x',ω)是由散射点x'传播到接收点x_r的Green函数。m(x')是x'点的类反射系数，即最小二乘偏移成像的目标。

需要说明的是，在各向异性介质中m(x')不仅仅与速度有关还受到各向异性参数ε,δ的影响。U₀(x'|x_s,ω)是由震源传播到散射点x'的入射波场，U_s(x_r|x_s，ω)是在检波点x_r接受到的一次反射波。结合公式(1)，具体的反偏移过程如下：

S31，根据速度场和各向异性参数，从地表开始，对震源进行逐层延拓，在地下每一点记录并存储入射波场，直到延拓到最深层。

具体的说，根据速度场和各向异性参数，以及震源(震源可由测数据得到)，从地表开始，利用高阶傅里叶有限差分(即FFD)算子震源进行逐层延拓得到入射波场，在地下每一点记录并存储入射波场，直到延拓到最深层。

在数值模拟的过程中，会对地下介质进行网格剖分，所述的每一点在模拟过程中是指网格的交点，其物理含义为地下介质的每一点。而所述的延拓即指的是模拟地震波在地下的传播情况。所述的高阶FFD波场延拓包括频率波数域的相移，频率空间域时移和频率空间域有限差分补偿，具体方式如下：

I，频率波数域的相移:

U₁(k_x,z_j+1；ω)＝U(k_x,z_j；ω)exp(ik_z0Δz) (2)

其中U(k_x,z_j；ω)是深度z_j处入射的总波场，j表示地层深度序号，k_x表示水平波数，ω表示圆频率，相移算子exp(ik_z0Δz)中k_z0表示背景垂直波数，i表示虚数单位，Δz表示延拓深度间隔。经过相移处理延拓得到的出射波场U₁(k_x,z_j+1；ω)为背景波场。

式(3)中c(z),ε₀(z),δ₀(z)是背景速度和各向异性参数。他们只随深度z发生变化。

II，频率空间域时移:

U₂(x,z_j+1；ω)＝U₁(x,z_j+1；ω)exp(iωΔsΔz), (4)

其中U₁(x,z_j+1；ω)是由U₁(k_x,z_j+1；ω)经过有关x的傅里叶反变换得到的。是介质的慢度扰动。时移校正主要是对小角度入射的波场由于背景速度和真实速度之间的偏差引起的波场传播误差进行补偿。

III，频率空间域有限差分补偿:

上式在计算时按照级联算法分两步进行具体计算

上计算得到的波场U(x,z_j+1；ω)即为z_j+1处最终的出射波场。这里的有限差分补偿主要是对由于各向异性参数偏差引起的高角度传播的波场误差进行校正。其中i是虚数单位，l₁，l₂，l₃是各向异性参数的函数。具体转换公式如下：

ξ＝(1+2ε),η＝2(ε-δ),ξ₀＝(1+2ε₀),η₀＝2(ε₀-δ₀) (7)

换言之，由公式(3)、(4)、(6)联合得到出射波场。需要说明的是，波场从上一层出射，然后入射到下一层，上一层的出射波场就是下一层的入射波场，

S32，根据每一点的入射波场与反射系数，获得二次震源。

即会获得每一点上各自的二次震源，二次震源的作用是形成反射波，因此加载反射波的过程已经包含了二次震源的加载。

S33，在最深层将反射波初始化为0，向地表逐层延拓，在延拓过程中，反射波场叠加来自每一层的反射波，直至地表；

具体的说，在最深层将反射波初始化为0，并利用高阶FFD算子逐层延拓，在延拓过程中，反射波场叠加上来自每一层的反射波，直至地表；这里在采用高阶FFD算子逐层延拓时也包括如步骤S31中所述的频率波数域的相移，频率空间域时移和频率空间域有限差分补偿，这里就不再重复说明。

S34，根据延拓到地表时的反射波场通过傅里叶反变换得到模拟数据。

S4，计算观测数据和模拟数据的数据残差，计算方法为模拟数据Lm_k减去观测数据d^obs即可。

S5，根据数据残差获得最终成像结果，具体的方法如下：

S51，判断数据残差是否小于给定阈值，如果是，则输出反射系数作为成像结果；否者进入步骤S52。

S52，根据数据残差得到梯度，具体的说，包括如下步骤：

根据速度场v和各向异性参数，加载震源和当前数据残差，g_k+1＝L^T(Lm_k-d_obs)，其中L^T表示偏移算子；

利用傅里叶有限差分正向延拓算子逐层延拓入射波场，利用反向延拓算子逐层延拓数据残差；

对延拓到每一层的入射波场取共轭，并与延拓到相应位置的观测记录按照以下公式所示的互相关成像条件提取成像值：

其中*表示取共轭；

最后根据成像值获得梯度值g。

利用傅里叶有限差分进行计算综合了效率和精度两方面的优势，傅里叶有限差分算法比现有技术相比精度和效率更高。

在这个过程中的利用FFD算子逐层延拓也包含如步骤S31和S33中所述的频率波数域的相移，频率空间域时移和频率空间域有限差分补偿，这里就不再重复说明。

S53，根据梯度计算共轭梯度方向的修正参数。

具体的说，根据梯度计算共轭梯度方向的修正参数其中P^-1表示照明预条件算子,T表示转置，共轭梯度是由梯度计算得到的，共轭梯度具有更好的全局收敛性，梯度仅具有较好的局部收敛性。

S54，根据修正参数计算共轭梯度方向和共轭梯度更新步长。

即，计算共轭梯度方向z_k+1＝P^-1g_k+1+βz_k，和共轭梯度更新步长

S55，利用共轭梯度方向和步长更新反射系数，并返回步骤S3。

即，利用共轭梯度方向和步长更新反射系数：m_k+1＝m_k-αz_k+1，并返回步骤S3，重新开始迭代，直至满足S51中的终止条件数据残差小于给定阈值。最终得到的成像结果如图3b所示。

传统的不考虑各向异性的方法所得到的结果图3a所示，对模型中间部分强各向异性区域，偏移结果存在明显的绕射波不收敛现象，强各向异性体的上边界几乎没有对应同相轴。而对比本发明的结果图3b来看，无论是对高陡的岩丘左侧边界，高陡断层，还是岩丘下深部层状构造均达到了比较理想的成像效果。同时子波被压缩，剖面整体的同相轴更加细致清晰，分辨率有明显提高。

通过采用本发明提供的方法，充分考虑了地下介质的各向异性这一实际情况，能够更加精确的模拟地震波在真实地下介质中的传播规律，基于线性反演思想通过反复迭代，本发明进一步提高偏移成像结果的分辨率和保真度，使得偏移成像结果中不同深度的同相轴振幅更加均衡，成像结果更为准确。此外，本发明继承了传统单程波类偏移成像方法的优点，占用较小的存储空间，偏移成像噪音少，对速度误差的敏感性更弱。

当然，上述实施例说明并非是对本发明的限制，本发明不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (3)

1.一种各向异性介质的成像方法，其特征在于，包括如下步骤：

S1，获得观测数据，根据观测数据获得速度场和各向异性参数，其中，所述各向异性参数包括纵波水平速度和垂向速度关系的函数，和，纵波各向异性强度，所述纵波是根据所述观测数据得到的；

S2，对观测数据进行偏移得到反射系数；

S3，根据速度场和各向异性参数，对反射系数进行反偏移得到模拟数据，具体包括：

根据速度场和各向异性参数，从地表开始，对震源进行逐层延拓，在地下每一点记录并存储入射波场，直到延拓到最深层；

根据每一点的入射波场与反射系数，获得二次震源；

在最深层将反射波初始化为0，向地表逐层延拓，在延拓过程中，反射波场叠加来自每一层的反射波，直至地表；

根据延拓到地表时的反射波场得到模拟数据；

其中，在延拓过程中包括频率波数域的相移、频率空间域时移和频率空间域有限差分补偿，具体包括如下方式：

对频率波数域的相移，U₁(k_x,z_j+1；ω)＝U(k_x,z_j；ω)exp(ik_z0Δz)，其中U(k_x,z_j；ω)是深度z_j处入射的总波场，j表示地层深度序号，k_x表示水平波数，ω表示圆频率，相移算子exp(ik_z0Δz)中k_z0表示背景垂直波数，i表示虚数单位，Δz表示延拓深度间隔，其中，c为背景速度场，ε₀为描述纵波水平速度和垂向速度关系的参数，δ₀为描述纵波各向异性强度的参数，c、ε₀和δ₀均只随深度z发生变化；

频率空间域时移，U₂(x,z_j+1；ω)＝U₁(x,z_j+1；ω)exp(iωΔsΔz)，其中U₁(x,z_j+1；ω)是由U₁(k_x,z_j+1；ω)经过有关x的傅里叶反变换得到的波场值，Δs为介质的慢度扰动；

频率空间域有限差分补偿:l₁，l₂，l₃是各向异性参数的函数，U(x,z_j+1；ω)即为z_j+1处最终的出射波场，v为速度场；

S4，计算观测数据和模拟数据的数据残差；

S5，根据数据残差获得最终成像结果。

2.如权利要求1所述的一种各向异性介质的成像方法，其特征在于，所述步骤S5中，根据数据残差获得最终成像结果，具体包括：

S51，判断数据残差是否小于给定阈值，如果是，则输出反射系数作为成像结果；否者进入步骤S52；

S52，根据数据残差得到梯度；

S53，根据梯度计算共轭梯度方向的修正参数；

S54，根据修正参数计算共轭梯度方向和共轭梯度更新步长；

S55，利用共轭梯度方向和步长更新反射系数，并返回步骤S3。

3.如权利要求2所述的一种各向异性介质的成像方法，其特征在于，所述步骤S52中，根据数据残差得到梯度，具体包括：

S61，根据速度场v和各向异性参数，加载震源和当前数据残差；

S62，利用傅里叶有限差分正向延拓算子逐层延拓入射波场，利用反向延拓算子逐层延拓数据残差；

S63，对延拓到每一层的入射波场取共轭，并与延拓到相应位置的观测记录按照如下公式提取成像值；

S64，根据成像值获得梯度值。