CN108181657A - 全波形反演梯度计算中分离偏移和层析成像模式的方法 - Google Patents

Abstract

本发明涉及一种全波形反演梯度计算中分离偏移和层析成像模式的方法，包括以下内容：1)针对目标设计角度滤波器f(α)，初始化角度滤波梯度为0；2)对在某时间片要进行互相关的正传震源波场和反传余量波场分别沿各个空间方向进行傅里叶变换，相应得到波数域波场；3)对两波数域波场分别通过角度滤波器f(α)进行加权计算；4)重复步骤2)和步骤3)通过角度滤波器f(α)完成所有时间片的加权计算，累加得到最后的梯度；5)对步骤4)得到的梯度值沿空间方向进行反傅里叶变换并乘以系数，得到实现偏移成像和层析成像模式分离的梯度值。

Description

全波形反演梯度计算中分离偏移和层析成像模式的方法

技术领域

本发明是关于一种全波形反演梯度计算中分离偏移和层析成像模式的方法，涉及地震数据处理技术领域。

背景技术

全波形反演是利用地震数据来重建地球内部介质性质和结构的有效手段，它是通过优化受波动方程约束的目标函数来实现的。然而目前全波形反演对大、中等尺度的背景速度的重构都需要大偏移距的地震数据，这是因为在大偏移距的条件下折射波很丰富，从而可以在全波形反演中实现层析成像来消除模型中的大、中等尺度的误差。传统的采集观测系统受偏移距的限制，其基本是由反射波占主导的，例如海上拖缆数据。对这种以反射波为主导的地震数据进行全波形反演本质上相当于做非线性最小二乘逆时偏移，因此在全波形反演中偏移成像模式主导了反演，虽然层析成像模式也起了一定的作用，但是比起偏移成像模式在反演中的作用却是很微弱的。为了突出其中某一种模式的作用，非常有必要将全波形反演中的两种模式分离开来。

实现全波形反演中的偏移成像模式和层析成像模式的分离目前有以下几项技术：将模型参数化为大尺度的背景速度和小尺度的波阻抗或者反射系数，从而使背景速度的梯度主要是由层析成像分量组成，而波阻抗或反射系数的梯度主要是由偏移成像分量组成，从而实现全波形反演中这两种模式的分离。但是这种方法不能完全彻底地分离波场中散射波和入射波，因此梯度中的这两种分量也不能完全分离。另一种方法是将波场中的上行波场和下行波场、左行波场和右行波场分开，再分别进行互相关来产生梯度中层析成像分量和偏移成像分量。然而这种方法假设散射波和入射波沿不同的方向传播，比如，入射波向下传播而散射波向上传播。但是只从上下或左右来划分传播方向是不精确的，很多情况下，散射波和入射波会沿着这四个方向中的相同方向传播。因此这种方法在很多情况下也不能正确地分离梯度中的这两种分量。

逆散射成像方法计算出体模量(bulk modulus)的梯度和密度的梯度，然后加权相减得到速度的梯度，加权相加得到波阻抗的梯度。然而，其对速度的梯度沿张开角度进行了(1-cos²θ/2)的加权，对波阻抗梯度沿张开角度进行了cos²θ/2的加权，θ是张开角(其范围为0到360度)。这里，因为梯度中的层析成像分量是在张开角接近180度时形成，而偏移分量在其它张开角形成，所以速度梯度中的层析成像分量得到加强，波阻抗梯度中的偏移分量得到了加强，从而实现这两个分量的分离。但是因为这种余弦平方函数的加权变化缓慢，并包含了所有角度，所以存在两种分量的互相泄露。根据震源波场和余量波场做互相关的入射波和回传波的张开夹角来从常规的梯度中直接提取出这两种分量。当张开夹角接近180度时，由这两个波互相关产生的梯度是由层析成像分量主导的，否则是由偏移分量主导的。因此可以设计一种角度滤波来实现这两种分量的分离，现有技术中存在时间延迟域(Time-lag extendeddomain)的角度滤波、坡印廷矢量(Poynting vector)法以及局部倾角叠加(local slant-stack)。其中，时间延迟域的角度滤波方法中假设速度是常速度，克服这种困难需要通过坐标变换，并且该方法假设介质中的速度是各项同性，因此其计算比较复杂和不精确。坡印廷矢量法中地震波传播的方向可以通过坡印廷矢量来进行计算。因此，可以通过坡印廷矢量来计算入射波和回传波传播的方向，从而来计算这两个波在每个点的夹角，进而实现全波形反演梯度中两种成分分量的分离。但是在一个点上坡印廷矢量是所有波传播方向的一种综合表象，当多个方向的波叠加在一块儿时，坡印廷矢量不能给出每个波的传播方向，因此这种方法对复杂波场的层析成像分量和偏移分量的分离不是有效的。局部倾角叠加通过在空间滑动的窗口内实现波场沿各个传播方向的分离，然后实现入射波和回传波夹角的估计，从而实现角度滤波。但是这种方法假设在有限大小窗口内的波速是常数，并且由于窗口较小，导致分辨率较低，因此这种方法也不是非常精确。

发明内容

针对上述问题，本发明的目的是提供一种全波形反演梯度计算中分离偏移和层析成像模式的方法，能够高效精确区分偏移成像分量和层析成像分量。

为实现上述目的，本发明采取以下技术方案：一种全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于包括以下内容：1)针对目标设计角度滤波器f(α)，初始化角度滤波梯度g(k)为0；2)对在某时间片t_k要进行互相关的正传震源波场和反传余量波场分别沿各个空间方向进行傅里叶变换，相应得到波数域波场；3)对两波数域波场分别通过角度滤波器f(α)进行加权计算；4)重复步骤2)和步骤3)通过角度滤波器f(α)完成所有时间片的加权计算，累加得到最后的梯度；5)对步骤4)得到的梯度值沿空间方向进行反傅里叶变换并乘以系数，得到实现偏移成像和层析成像模式分离的梯度值。

进一步地，所述步骤3)对两波数域波场通过角度滤波器f(α)进行加权计算，具体过程为：

首先，对两波数域的波场P_s(k_s)和P_r(k_r)进行带滤波的乘积运算：

其次，将上式结果累加到g(k)中都得到新的角度滤波梯度：

g(k)＝g(k)+g_k(k)。

为实现上述目的，本发明采取以下技术方案：一种全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于包括以下内容：

1)针对目标设计角度滤波器f(α)，初始化角度滤波梯度g(x)为0；

2)对在时间片t_k要进行互相关的正传震源波场和反传余量波场沿各个空间方向进行傅里叶变换，得到两波数域波场；

3)对两波数域波场通过角度滤波器f(α)进行带滤波的乘积运算，并对运算结果进行沿空间方向的反傅里叶变换；

4)将反傅里叶变换结果进行梯度累加；

5)重复步骤2)到步骤4)，累加得到最后的梯度；

6)对步骤5)得到的最后梯度值乘以系数，得到实现偏移成像和层析成像模式分离的梯度值。

进一步地，步骤3)对两波数域波场通过角度滤波器f(α)进行带滤波的乘积运算：式中，k_si和k_rj分别表示第i个震源波场的波数和第j个余量波场的波数，p_s是正传震源波场，p_r是反传余量波场，i、j是波数的序号。

进一步地，所述步骤4)将反傅里叶变换结果进行梯度累加：

g(x)＝g(x)+g_k(x)。

进一步地，所述步骤1)针对目标设计角度滤波器f(α)具体过程为：

如果要突出偏移分量，则设计为小角度通过的角度滤波器，如果要突出层析成像分量，则设计为大角度通过的角度滤波器，如果要提取角度道集，则设计为很窄的带通的角度滤波器，其中，小角度通过滤波器，其过度带140°到160°，开口角度大于160度被滤掉；大角度通过滤波器，其过度带150°到170°，角度小于150°被滤出；很窄的带通的角度滤波器以60°为中心的，其过度带是50°到70°，其它角度被滤出。

本发明由于采取以上技术方案，其具有以下优点：1、本发明通过傅里叶变换把震源波场和余量波场映射到波数域，然后通过选取不同张开角的角度滤波器来区分偏移成像分量和层析成像分量，它不受模型速度变化、大倾角反射界面、波传播方向以及速度各向异性的影响，可以取得非常好的效果，能够高效精确区分偏移成像分量和层析成像分量。2、本发明的平面波域的角度滤波不仅可以用来提取全波形反演梯度中的层析成像分量和偏移成像分量，还可以用来滤除逆时偏移中的低频噪音，同时可以用于构建共成像点的角度道集。

附图说明

图1是本发明在平面波域产生梯度过程中波数映射的示意图；

图2是本发明全波形反演梯度中偏移成像和层析成像模式分离流程图；

图3是本发明中角度滤波器f(α)示意图，图(a)是大角度通过的角度滤波器，图(b)是小角度通过的角度滤波器，图(c)是很窄的带通的角度滤波器；

图4是根据本发明的一个实施例的两层速度模型，检波器和炮点位于上层速度模型中；

图5是本发明实施例效果示意图，(a)是常规的梯度；(b)是0°to 180°的范围内经过小角度滤波后的梯度；(c)是0°to 180°的范围内大角度滤波后的梯度；(d)是100°to120°范围内经过小角度滤波后的梯度；(e)100°to 120°范围内经过大角度滤波后的梯度；(f)是140°to 170°范围内经过大角度滤波后的梯度，(b)和(d)显示的是梯度中的偏移分量，(c)、(e)和(f)显示的是梯度中的层析成像分量。

具体实施方式

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

波动方程可以表示为

Ap＝s (1)

式中，A是波动方程算子，p是波场，s是震源子波。

全波形反演可以表示为优化(最大或最小化)一个目标函数φ。目标函数φ是一个范数，具有多种形式，例如预测数据和观察数据差的二范数平方。

式中，T表示转置，d₀是在空间和时间离散点观测到的数据，d是预测数据，δd是余量。

目标函数对模型m的导数，即梯度可以表示为：

式中，p_s是正传震源波场，是波动方程算子对模型的导数，p_r＝A^-TD^Tδd是反传余量波场，D是选择矩阵使得d＝Dp。

从公式(3)可以看出，全波形反演的梯度是由正传震源波场与反传余量波场做零延迟互相关并乘上得到的。公式(3)中可以包含算子部分和系数部分。在算子部分作用波场后，目标函数对模型m的梯度可以表示为：

式中，P_s和P_r是经中算子部分作用后的波场，β是中系数。

公式(4)在频率域可以表示为：

在时间域可以表示为：

公式(5)和(6)是常规方法计算全波形反演的梯度，也是地震偏移成像的成像条件。

全波形反演梯度中的层析成像分量其实在地震偏移中被看成低频噪音。因此此处讨论的方法也适合地震偏移成像中低频噪音的去除及共成像点角度道集的生成。因为本发明的全波形反演梯度计算中分离偏移成像和层析成像模式的方法是不依赖时间和频率，所以对于时间域和频域内的角度滤波方法是一样的。

下面通过具体实施例详细推导在时间域通过角度滤波来实现层析成像分量和偏移分量的分离，偏移成像中低频噪音的去除和共成像点角度道集的生成。对频率域的实施只需要将时间域实施中的变量替换为频率。

沿空间方向x(x是空间坐标轴变量，对2维模型其包括(x,y)两个分量，对3维模型其包括(x,y,z)3个分量)的各个分量对波场进行空间域傅里叶变换得到：

式中，*表示两个波场中任意两个元素之间的乘积，k_s、k_r和k分别表示震源波场，余量波场和梯度的波数，i、j是波数的序号，k_si和k_rj分别表示第i个震源波场的波数和第j个余量波场的波数。在傅里叶波数域，每一个元素代表一个有特定波数、振幅和相位的平面波，因此通过空间傅里叶变换，波场被映射到平面波域，因此此处的波数域即一种平面波域。在傅里叶域，梯度是由公式(7)做乘积再求和得到的。例如，把震源波场的一个平面波和反传余量波场的一个平面波做乘积，即

则梯度中的波数满足以下关系式

k＝k_s+k_r (9)

乘积的波数映射如图1所示(k有x和z两个分量),根据公式(8)和(9)可以精确的计算出公式(7)所表达的梯度。

梯度中层析成像分量在当震源波场与余量波场中波的传播方向相反时，即图1所示的张开角α等于180度时形成。考虑到地震数据中频带是有限的，为层析成像部分贡献的波的张开角大约在180度附近。相反，由于偏移部分是由反射波所贡献的，为偏移部分贡献的波的张开角要小得多。因此，可以对公式(7)引入一个角度滤波器f(α)，则带角度滤波的梯度公式可以表示为

式中，角度α可以通过下式求取

因此，通过公式(10)和(11)，全波形反演梯度中偏移成像分量和层析成像分量能被分离，其具体实施将在下面进行详细介绍。在公式(10)中，可以先进行各时间片的计算和累加，然后再进行反傅里叶变换回到空间域，然后再乘以系数β。因此，公式(10)是适合变速度的模型。

如图2所示，基于上述原理，本发明全波形反演梯度计算中分离偏移和层析成像模式的方法，包括以下内容：

1、针对目标设计一个角度滤波器f(α)，初始化角度滤波梯度g(k)为0。

如果要突出偏移分量，则设计为小角度通过的角度滤波器。

如果要突出层析成像分量，则设计为大角度通过的角度滤波器。

如果要提取角度道集，则设计为很窄的带通的角度滤波器。

角度滤波器的设计可以根据实际需要设置，本发明是设计一个随角度变化的线性滤波器，而一对来自震源波场和余量波场的平面波形成的角度可以通过公式(11)计算得到，如图3所示，图(a)是大角度通过滤波器，其过度带150°到170°，角度小于150°被滤出；图(b)是小角度通过滤波器，其过度带140°到160°，开口角度大于160度被滤掉；图(c)是以60°为中心的窄带通过滤波器，其过度带是50°到70°，其它角度被滤出。

2、对第k个时间片t_k，要进行互相关的正传震源波场P_s(x)和反传余量波场P_r(x)沿各个空间方向进行傅里叶变换，得到波数(平面波)域的波场P_s(k_s)和P_r(k_r)。

3、对波数域波场P_s(k)和P_r(k)通过角度滤波器f(α)进行加权计算，具体过程为：

对波数域的波场P_s(k_s)和P_r(k_r)进行带滤波的乘积运算：

将上式结果累加到g(k)中，具体过程：

g(k)＝g(k)+g_k(k) (13)

4、重复步骤2和3通过角度滤波器f(α)完成所有时间片t_k的计算，累加得到最后的梯度g(k)。

5、对步骤4计算得到的梯度g(k)进行沿空间的反傅里叶变换并乘以系数β得到实现梯度中偏移成像和层析成像模式分离的梯度g(x)。

在一个优选的实施例中，上述步骤5中的反傅里叶变换也可以在完成步骤3的对波数域的波场P_s(k_s)和P_r(k_r)进行带滤波的乘积运算后进行，然后将反傅里叶变换结果进行梯度累加：

g(x)＝g(x)+g_k(x)。

对得到的最后梯度值乘以系数β，得到实现偏移成像和层析成像模式分离的梯度值。

下面通过具体实施例对本发明的全波形反演梯度计算中分离偏移和层析成像模式的方法进行详细说明。

实施例：双层模型

本实施例是将本发明的这种偏移和层析成像模式分离应用于一个两层速度模型，如图4所示。震源和检波器位于上层速度模型中，初始速度模型采用一个低于真实速度的两层模型，其中直达波在观测数据和预测数据中均被滤掉。由于初始模型存在反射界面，因此震源波场和余量波场均存在反射波。图5(a)是常规的梯度，其中包含了层析成像分量和偏移分量的全波形反演的梯度。偏移分量呈椭圆形，其上的每一点从震源到检波点具有相同的走时。通过把所有震源和检波点的成分相叠加来构建全波形反演的梯度信息可以矫正界面位置但是无法更新背景速度。其余的梯度部分是层析成分，它提供了波形反演更新背景速度所需要的信息，在地球物理界俗称为“兔子耳朵”。图5(b)展示了过渡带从0°到180°的低角度通过滤波后的偏移成分，过渡带具有余弦平方形状。这个滤波器实质上等价于Laplacian滤波或者逆散射成像条件(Whitmore和Crawley，2012)。图5(c)与图5(b)相对应的高角度通过滤波。这个滤波器与逆散射成像条件对加强层析成分具有相同的作用。可以看到偏移成分被压制了，因此更突显了兔子耳朵部分。然而，由于采用了较宽的过渡带，这种分离是不理想的，可以看到偏移成分中一些被弱化的水平部分，如图5(b)错误地被贡献到层析成分如图5(c)。由于逆散射级数是在空间域实现的，使得对这一问题的改进变得相当困难。但是本发明的方法是在平面波域实施的，通过张开角度来从梯度中筛选地震波就变得非常自然而有效了。图5(d)显示了过渡带从100°到120°的低角度通过滤波后的偏移分量。与图5(b)相比，滤波没有减弱水平成分的振幅。图5(e)显示了过渡带从120°到140°高角度通过滤波后的层析分量，与图5(c)相比，表现出了对两种模式更好的分离结果。图5(f)展示了将过渡带进一步提高到140°到170°的分离结果。

上述各实施例仅用于说明本发明，其中方法的各实施步骤都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (6)

1.一种全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于包括以下内容：

1)针对目标设计角度滤波器f(α)，初始化角度滤波梯度g(k)为0；

2)对在某时间片t_k要进行互相关的正传震源波场和反传余量波场分别沿各个空间方向进行傅里叶变换，相应得到波数域波场；

3)对两波数域波场分别通过角度滤波器f(α)进行加权计算；

4)重复步骤2)和步骤3)通过角度滤波器f(α)完成所有时间片的加权计算，累加得到最后的梯度；

5)对步骤4)得到的梯度值沿空间方向进行反傅里叶变换并乘以系数，得到实现偏移成像和层析成像模式分离的梯度值。

2.如权利要求1所述的全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于，所述步骤3)对两波数域波场通过角度滤波器f(α)进行加权计算，具体过程为：

首先，对两波数域的波场P_s(k_s)和P_r(k_r)进行带滤波的乘积运算：

其次，将上式结果累加到g(k)中都得到新的角度滤波梯度：

g(k)＝g(k)+g_k(k)。

3.一种全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于包括以下内容：

1)针对目标设计角度滤波器f(α)，初始化角度滤波梯度g(x)为0；

2)对在时间片t_k要进行互相关的正传震源波场和反传余量波场沿各个空间方向进行傅里叶变换，得到两波数域波场；

3)对两波数域波场通过角度滤波器f(α)进行带滤波的乘积运算，并对运算结果进行沿空间方向的反傅里叶变换；

4)将反傅里叶变换结果进行梯度累加；

5)重复步骤2)到步骤4)，累加得到最后的梯度；

6)对步骤5)得到的最后梯度值乘以系数，得到实现偏移成像和层析成像模式分离的梯度值。

4.如权利要求3所述的全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于：步骤3)对两波数域波场通过角度滤波器f(α)进行带滤波的乘积运算：式中，k_si和k_rj分别表示第i个震源波场的波数和第j个余量波场的波数，p_s是正传震源波场，p_r是反传余量波场，i、j是波数的序号。

5.如权利要求3所述的全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于：所述步骤4)将反傅里叶变换结果进行梯度累加：

g(x)＝g(x)+g_k(x)。

6.如权利要求1或3所述的全波形反演梯度计算中分离偏移和层析成像模式的方法，其特征在于：所述步骤1)针对目标设计角度滤波器f(α)具体过程为：

如果要突出偏移分量，则设计为小角度通过的角度滤波器，如果要突出层析成像分量，则设计为大角度通过的角度滤波器，如果要提取角度道集，则设计为很窄的带通的角度滤波器，其中，小角度通过滤波器，其过度带140°到160°，开口角度大于160度被滤掉；大角度通过滤波器，其过度带150°到170°，角度小于150°被滤出；很窄的带通的角度滤波器以60°为中心的，其过度带是50°到70°，其它角度被滤出。