CN115755178B - 一种基于积分地震子波的时间域全波形反演方法 - Google Patents

Abstract

本发明属于海洋地震勘探领域，公开了一种基于积分地震子波的时间域全波形反演方法，包括如下步骤：输入野外观测地震数据和初始速度参数模型；根据野外观测地震数据提取地震子波，对地震子波进行时间积分，构建积分地震子波；以积分地震子波作为震源，数值求解时间域波动方程得到模拟地震数据；将野外观测地震数据与模拟地震数据逐点做差，以残差的时间积分作为震源，再次数值求解时间域波动方程，得到残差波场，同时计算梯度算子；计算速度参数模型更新量迭代反演，终止迭代并输出反演结果。本发明构建积分地震子波生成模拟地震记录，解决了时间域全波形反演对低频数据的压制问题，在梯度中增加了中低波数分量比重，提高了参数反演的精度。

Description

一种基于积分地震子波的时间域全波形反演方法

技术领域

本发明属于海洋地震勘探领域，特别涉及一种基于积分地震子波的时间域全波形反演方法。

背景技术

海洋地震反演通过观测到的地震数据估算相应的地球物理参数，进而反推海底地下的结构形态及物质成分，可以有效地识别地质构造、预测自然灾害和勘探油气藏。地震波的传播速度参数不仅是处理、解释地震数据的主要依据，而且其本身就是反映底下介质构造和岩性的重要数据。因此，求取地震波在地下传播的精确速度显得尤为重要。

全波形反演以波动理论为基础，匹配模拟地震数据与野外观测地震数据的波形信息建立目标函数，利用局部最优化的方法迭代寻找目标函数的全局最小值，可以高分辨率地恢复海底底质的速度参数，是目前海洋地震勘探领域的前沿方向之一。全波形反演可以在时间域实现，也可以在频率域实现，频率域全波形反演对计算机的要求更为苛刻，所以适用性更好的当属时间域全波形反演。但是，时间域全波形反演存在对低频数据的压制问题，不利于中低波数参数分量的恢复，容易使反演陷入局部极小值而导致失败。因此，解决时间域全波形反演这种低频压制问题不仅可以增加反演的容错率，同时可以提高反演的精度。

发明内容

为解决上述技术问题，本发明提供了一种基于积分地震子波的时间域全波形反演方法，以达到提高反演精度的目的。

为达到上述目的，本发明的技术方案如下：

一种基于积分地震子波的时间域全波形反演方法，包括如下步骤：

步骤一，输入数据：

输入野外观测地震数据 d _obs (x, t)和初始速度参数模型m₀；

步骤二，构建积分地震子波：

根据野外观测地震数据 d _obs (x, t)提取地震子波，对地震子波进行时间积分，构建积分地震子波：

                         (1)

式中， t表示时间， f(t)表示根据野外观测地震数据提取的地震子波，表示时间积分后的地震子波，即积分地震子波；

步骤三，正演模拟：

以积分地震子波作为震源，数值求解时间域波动方程得到所有时间点的正演模拟波场 u(x, t)，从正演模拟波场 u(x, t)中得到接收点处的数据，即模拟地震数据 d _cal (x, t)，时间域波动方程表示形式如下：

                               (2)

式中， u(x, t)表示正演模拟波场，x表示空间位置矢量，m表示速度参数模型，初始值为m₀；表示空间导数之和；

步骤四，梯度算子计算：

将野外观测地震数据 d _obs (x, t)与模拟地震数据 d _cal (x, t)逐点做差，则残差表示为 d _obs (x, t)- d _cal (x, t)，以残差的时间积分作为震源替换方程(2)中的积分地震子波，然后再次数值求解时间域波动方程(2)，得到残差波场 q(x, t)，同时根据方程(3)计算梯度算子；

                            (3)

式中， q(x, t)表示残差波场， s和 g分别表示产生野外观测地震数据的震源点和接收点；

步骤五，计算速度参数模型更新量迭代反演：

先给予一个试探步长，再用Armijo条件的一维线搜索方法求取合适的步长作用于梯度算子计算本次迭代的速度参数模型更新量，保证每次迭代的速度参数模型更新量在20m/s-100m/s之间，利用方程(4)计算更新后的速度参数模型m：

                                      (4)

式中，m₀表示初始速度参数模型， λ表示步长，表示本次迭代的速度参数模型更新量；

返回步骤三，将更新后的速度参数模型m代入方程(2)进行正演模拟，重新计算模拟地震数据，以野外观测地震数据与模拟地震数据残差的平方和[ d _obs (x, t)- d _cal (x, t)]²作为判定条件，判断反演是否收敛，进而选择继续迭代更新或是终止迭代并输出反演结果。

上述方案中，步骤四中，当判定条件小于门限值或连续5次增加，则认为反演已收敛。

优选地，所述门限值为0.001。

通过上述技术方案，本发明提供的一种基于积分地震子波的时间域全波形反演方法具有如下有益效果：

常规的时间域全波形反演所用的子波是直接从野外观测地震数据中提取得到的，在计算梯度算子时，地震数据相当于通过了一个高频滤波器，低频信息被压制，进而速度参数的中低波数分量得不到有效的更新。本发明构建积分地震子波生成模拟地震记录，解决了时间域全波形反演对低频数据的压制问题，在梯度中增加了中低波数分量比重，提高了参数反演的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为本发明实施例所公开的一种基于积分地震子波的时间域全波形反演方法流程示意图。

图2为Marmousi2速度参数模型，其中(a)为真实速度参数模型，(b)为初始速度参数模型。

图3 为地震子波与积分地震子波频带对比，其中，黑色实线为期待的频带，左边灰色实线为积分地震子波的频带，右边灰色实线为常规地震子波的频带。

图4 为反演速度参数模型对比，其中(a)为常规地震子波时间域全波形反演结果，(b)为积分地震子波时间域全波形反演结果。

图5 为单道反演速度曲线对比，其中(a)为常规地震子波时间域全波形反演结果，(b)为积分地震子波时间域全波形反演结果，其中黑色实线为真实速度曲线，灰色实线分别为反演速度曲线和初始速度曲线。

图6为梯度对比，其中(a)为常规地震子波时间域全波形反演梯度，(b)为积分地震子波时间域全波形反演梯度。

图7为不同距离处的单道梯度波数谱，其中(a)为距离0.8km处单道梯度波数谱，(b)为距离1.2km处单道梯度波数谱，(c)为距离1.4km处单道梯度波数谱，(d)为距离1.6km处单道梯度波数谱，其中黑色实线为常规子波时间域全波形反演梯度的波数谱，灰色实线为积分地震子波时间域全波形反演梯度的波数谱。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

本发明提供了一种基于积分地震子波的时间域全波形反演方法，如图1所示，包括如下步骤：

步骤一，输入数据：

输入野外观测地震数据 d _obs (x, t)和初始速度参数模型m₀；初始速度参数模型m₀是已知的，其速度值从浅层到深层渐变增加，本实施例中选用Marmousi2速度参数模型，如图2中(b)所示。真实速度参数模型是未知的，如图2中(a)所示，需要通过本发明的方法反演得到。

步骤二，构建积分地震子波：

根据野外观测地震数据 d _obs (x, t)提取地震子波，对地震子波进行时间积分，构建积分地震子波：

                     (1)

式中， t表示时间， f(t)表示根据野外观测地震数据提取的地震子波，表示时间积分后的地震子波，即积分地震子波。

如图3所示，右边灰色实线为根据野外观测地震数据提取地震子波的频带，左边灰色实线为积分地震子波的频带，积分地震子波的频带与期望的数据频带基本一致，而普通地震子波的频带整体偏高，低频部分被压制。

步骤三，正演模拟：

以积分地震子波作为震源，数值求解时间域波动方程得到所有时间点的正演模拟波场 u(x, t)，从正演模拟波场 u(x, t)中得到接收点处的数据，即模拟地震数据 d _cal (x, t)，时间域波动方程表示形式如下：

                               (2)

式中， u(x, t)表示正演模拟波场，x表示空间位置矢量，m表示速度参数模型，初始值为m₀；表示空间导数之和。

步骤四，梯度算子计算：

将野外观测地震数据 d _obs (x, t)与模拟地震数据 d _cal (x, t)逐点做差，则残差表示为 d _obs (x, t)- d _cal (x, t)，以残差的时间积分作为震源替换方程(2)中的积分地震子波，然后再次数值求解时间域波动方程(2)，得到残差波场 q(x, t)，同时根据方程(3)计算梯度算子；

                               (3)

式中， q(x, t)表示残差波场， s和 g分别表示产生野外观测地震数据的震源点和接收点。

步骤五，计算速度参数模型更新量迭代反演：

先给予一个试探步长，再用Armijo条件的一维线搜索方法求取合适的步长作用于梯度算子计算本次迭代的速度参数模型更新量，保证每次迭代的速度参数模型更新量在20m/s-100m/s之间，利用方程(4)计算更新后的速度参数模型m：

                                         (4)

式中，m₀表示初始速度参数模型， λ表示步长，表示本次迭代的速度参数模型更新量。

返回步骤三，将更新后的速度参数模型m代入方程(2)进行正演模拟，重新计算模拟地震数据，以野外观测地震数据与模拟地震数据残差的平方和[ d _obs (x, t)- d _cal (x, t)]²作为判定条件，判断反演是否收敛，当判定条件小于门限值或连续5次增加，则认为反演已收敛，这里门限值选择为0.001；进而选择继续迭代更新或是终止迭代并输出反演结果。

反演结果见图4，其中，图4中(a)为常规子波时间域全波形反演结果，图4中(b)为积分地震子波时间域全波形反演结果，积分地震子波由于解决了低频数据被压制的问题，其反演结果的波数谱更加连续。反演结果的抽道曲线对比见图5，图5中(a)为常规子波时间域全波形反演结果，图5中(b)为积分地震子波时间域全波形反演结果，可以看出当低频数据被压制的问题解决后，反演的结果更加贴近于真实速度模型，波数谱得到连续的恢复，反演精度得到明显改善。

为了从理论上说明本发明的积分地震子波的优势，将方程(3)变换到频率域可得：

                          (5)

式中，U^*(x,s)表示频率域正演模拟波场的共轭， Q(x,r,s)表示频率域残差波场， ω表示圆频率，s和r分别表示产生野外观测地震数据的震源点和接收点的位置矢量。

假设震源子波近似为单位振幅零相位，方程(5)可以写成格林函数的形式：

                        (6)

式中，G₀(x,s)和G₀(x,r)分别为震源点和接收点的格林函数，为野外观测地震数据和模拟地震数据之差。

定义背景速度参数 c ₀，入射单频率平面波传播方向，远场近似的散射单频率平面波方向，则有：

                                          (7)

式中， k ₀= ω/ c ₀。将方程(7)带入方程(6)中得：

          (8)

基于Born近似，可以写为：

                      (9)

式中， δm(x)表示模型参数m的扰动。

将方程(7)带入方程(9)，再将方程(9)带入方程(8)中可得：

        (10)

通过方程(10)可以看出，输入野外观测地震记录计算梯度相当于通过了一个高通滤波器， ω ⁴将压制地震数据的低频部分。

本发明构建积分地震子波方程(1)，变换到频率域：

                                          (11)

采用积分地震子波，此外，残差波场也进行时间积分，梯度算子(3)变换到频率域为：

                           (12)

方程(12)解决了梯度算子对低频地震数据的压制问题。图6和图7对比了不同的地震子波产生的模拟地震数据计算梯度算子波数分量的比重，可以看出若采用野外观测地震数据提取的地震子波计算梯度算子，会存在方程(10)的高通滤波器的问题，梯度算子的低波数分量被压制，见图6中(a)和图7中(a)-(d)中黑色实线；而采用积分地震子波，梯度算子的波数谱很明显的增加了中低波数的比重，见图6中(b)和图7中(a)-(d)中灰色实线。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (3)

1.一种基于积分地震子波的时间域全波形反演方法，其特征在于，包括如下步骤：

步骤一，输入数据：

输入野外观测地震数据d _obs (x,t)和初始速度参数模型m₀；

步骤二，构建积分地震子波：

根据野外观测地震数据d _obs (x,t)提取地震子波，对地震子波进行时间积分，构建积分地震子波：

    (1)

式中，t表示时间，f(t)表示根据野外观测地震数据提取的地震子波， 表示时间积分后的地震子波，即积分地震子波；

步骤三，正演模拟：

以积分地震子波作为震源，数值求解时间域波动方程得到所有时间点的正演模拟波场u(x,t)，从正演模拟波场 u(x,t)中得到接收点处的数据，即模拟地震数据d _cal (x,t)，时间域波动方程表示形式如下：

        (2)

式中， u(x,t)表示正演模拟波场，x表示空间位置矢量，m表示速度参数模型，初始值为m₀；表示空间导数之和；

步骤四，梯度算子计算：

将野外观测地震数据d _obs (x,t)与模拟地震数据d _cal (x,t)逐点做差，则残差表示为d _obs (x,t)-d _cal (x,t)，以残差的时间积分作为震源替换方程(2)中的积分地震子波，然后再次数值求解时间域波动方程(2)，得到残差波场q(x,t)，同时根据方程(3)计算梯度算子；

 (3)

式中，q(x,t)表示残差波场，s和g分别表示产生野外观测地震数据的震源点和接收点；

步骤五，计算速度参数模型更新量迭代反演：

先给予一个试探步长，再用Armijo条件的一维线搜索方法求取合适的步长作用于梯度算子计算本次迭代的速度参数模型更新量，保证每次迭代的速度参数模型更新量在20m/s-100m/s之间，利用方程(4)计算更新后的速度参数模型m：

    (4)

式中，m₀表示初始速度参数模型，λ表示步长，表示本次迭代的速度参数模型更新量；

返回步骤三，将更新后的速度参数模型m代入方程(2)进行正演模拟，重新计算模拟地震数据，以野外观测地震数据与模拟地震数据残差的平方和[d _obs (x,t)-d _cal (x,t)]²作为判定条件，判断反演是否收敛，进而选择继续迭代更新或是终止迭代并输出反演结果。

2.根据权利要求1所述的一种基于积分地震子波的时间域全波形反演方法，其特征在于，步骤四中，当判定条件小于门限值或连续5次增加，则认为反演已收敛。

3.根据权利要求2所述的一种基于积分地震子波的时间域全波形反演方法，其特征在于，所述门限值为0.001。

Images