CN111665556B - 地层声波传播速度模型构建方法 - Google Patents

Abstract

本发明涉及一种地层声波传播速度模型构建方法，包括：获取探测区域内震源震动的多个单炮数据，单炮数据包括单炮单道数据和单炮多道数据；根据地质任务需求，利用动态时窗对多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取初始速度模型，并基于初始速度模型对震源波形进行正演，得到正演模拟数据；根据观测数据和正演模拟数据计算波场残差；根据波场残差构造误差泛函；利用波场残差反传播到初始速度模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据根据伴随状态法计算误差泛函的梯度；利用误差泛函的梯度更新初始速度模型，得到精确速度模型。

Description

地层声波传播速度模型构建方法

技术领域

本发明涉及地球物理勘探技术领域，尤其涉及一种地层声波传播速度模型构建方法。

背景技术

在勘探地球物理技术领域,浅层建模技术一直是速度建模的难点。目前较为常用的是基于射线类的初至走时层析反演的方法。这类方法计算简单，不要求准确的背景速度场，且应用复杂地表情况，因此一直以来是解决近地表及浅层速度建模的常用方法。但这种方法是基于高频假设的最短路径原理，在低速体发育的情况下，该技术存在建模的‘盲区’。此外，射线类方法对射线密度敏感，当速度变化剧烈，即使是高速异常体，由于全发射现象的发生，将严重影响射线的密度，继而降低反演的精度。更何况，在实际资料应用中，初至的拾取工作量巨大，人工拾取存在误差，且当地表复杂的情况下初至很难拾取准确。

另外一种解决这类问题的新方法就是全波形反演的方法，这种方法基于波动方程，可以较为真实的模拟波的传播波场，因此，不受射线密度影响。但这一新方法的实现还有很多的限制，如要求大偏移距的观测系统、低频信息缺失等。全波形反演的实用化还有很长的路程，尤其陆上资料的实用化。

因此,目前缺少一种针对速度变化剧烈或低速体发育区域构建声波传播速度模型的解决方案。

发明内容

本发明的目的是针对现有技术中存在的缺陷，提供一种地层声波传播速度模型构建方法。

为实现上述目的，第一方面,本发明提供了一种地层声波传播速度模型构建方法，包括：

获取探测区域内震源震动的多个单炮数据，所述单炮数据包括单炮单道数据和单炮多道数据；

利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；

根据所述观测数据获取地层信息和均方根速度信息；

根据所述地层信息和均方根速度信息计算所述地层的各反射层的层速度，得到所述初始速度模型；

基于所述初始速度模型对震源波形进行正演，得到正演模拟数据；

根据所述观测数据和所述正演模拟数据计算波场残差，并根据所述波场残差构造误差泛函；

利用所述波场残差反传播到所述初始速度模型空间，得到残差反传播数据；

利用所述正演模拟数据和所述残差反传播数据根据伴随状态法计算所述误差泛函的梯度；

利用所述误差泛函的梯度更新所述初始速度模型，得到精确速度模型。

进一步的，所述利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据之后，所述方法还包括：

利用小波变换对所述观测数据进行滤波处理，得到处理后的观测数据；

对所述处理后的观测数据进行多次波压制处理。

进一步的，所述基于所述初始速度模型对震源波形进行正演，得到正演模拟数据具体包括:

利用交错网格有限差分法对所述震源波形的波动方程进行时间域离散,得到离散后的波动方程；

根据所述离散后的波动方程和所述初始速度模型确定交错网格各个时刻空间分布的波场值。

进一步的，所述根据所述观测数据和所述正演模拟数据计算波场残差具体包括：

将所述观测数据与所述正演模拟数据做差，得到所述波场残差。

进一步的，所述根据所述波场残差构造误差泛函具体包括：

根据公式

计算误差泛函，其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，d_obs为观测数据，b(m)-d_obs为波场残差，C_D为数据协方差矩阵，C_M为模型的协方差矩阵，m为速度初始模型的模型参数，m_prior为先验信息模型，λ为先验信息比重参数。

进一步的，所述利用所述波场残差反传播到所述初始速度模型空间，得到残差反传播数据具体包括：

将反传算子作用于波场残差,得到所述初始速度模型空间的残差反传播数据。

进一步的，所述利用所述正演模拟数据和所述残差反传播数据根据伴随状态法计算所述误差泛函的梯度具体包括：

根据公式

计算所述误差泛函的梯度；

其中，

k＝ρV_p ²，P_f为正演模拟数据，P_b残差反传播数据，ω为频率，V_P为速度，k、ρ为初始模型参数，E为误差泛函。

进一步的，所述方法还包括：

根据最快下降法、共轭梯度法、拟牛顿法LBFGS确定最优迭代步长。

进一步的，所述利用所述误差泛函的梯度更新所述初始速度模型，得到精确速度模型具体包括：

根据误差泛函的梯度和最优迭代步长对初始速度模型进行更新；

确定迭代终止条件；

当所述初始速度模型满足所述迭代终止条件时，得到精确速度模型。

本发明提供的地层声波传播速度模型构建方法，获取探测区域内震源波形的单炮数据，单炮数据包括单炮单道数据和单炮多道数据；利用时窗对单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取初始速度模型，并基于初始速度模型对震源波形进行正演，得到正演模拟数据；根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；利用波场残差反传播到初始速度模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据计算误差泛函的梯度；利用误差泛函的梯度更新初始速度模型，得到精确速度模型。本发明提供的方法，根据地质任务需求利用初至以后一段时间内波场的运动学和动力学特征完成波形反演，实现浅层的高精度速度建模。

附图说明

图1为本发明实施例提供的地层声波传播速度模型构建方法流程图；

图2为本发明实施例提供的高斯函数图；

图3为本发明实施例提供的二维声波各物理量和介质参数的网格配置示意图；

图4为本发明实施例提供的三维声波各物理量和介质参数的网格配置示意图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述，显然，所描述的实施例仅仅是本发明一部份实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本发明针对浅层建模的难点，尤其在速度变化剧烈，低速绕动体发育的情况下，避开基于射线理论的弱点，综合利用初至以后一段时间内的波场信息实现浅层的高精度建模。尤其解决山前带、黄土塬、沙漠地带等低速体发育的区域，为中深层建模及偏移成像提供可靠的支撑。

图1为本发明实施例一提供的地层声波传播速度模型构建方法流程图。如图1所示，具体包括如下步骤：

步骤101，获取探测区域内震源震动的多个单炮数据；

将地震勘探仪器布设于探测区域内，每放一炮采集一炮的单炮数据，单炮数据包括单炮单道数据和单炮多道数据。根据具体情况和施工需要，发射多炮以采集多组单炮数据。其中，地震勘探仪器具体为408ULS电缆仪器等。

其中，单炮数据为浅剖数据，利用高斯函数将采集到的单炮多道数据处理成单炮单道数据，取高斯函数的一半覆盖整条测线，然后将对应的离散的点，作为权重系数，将权重系数做归一化，使其和为1。整理成一个超道集，如公式(1)所示，从而将单炮多道数据处理成单炮单道数据。

其中，μ为炮点，i为检波点，σ为高斯窗口因子。

选取高斯函数图(如图2所示)中的x轴正半轴部分，横坐标对应检波点，每个点的系数为高斯系数M，如公式(2)所示。

本发明技术方案中的探测区域具体为山前带、黄土塬、沙漠地带等低速体发育的区域。

步骤102，利用时窗对单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；

具体的，在声波传播过程中加以时窗控制，得到近偏移距在近地表以及中浅层传播的波形信息。利用固定时长或非固定时长的移动时窗截取单炮数据中的直达波、浅层反射波、浅层折射波，得到观测数据。

所截取的波场并不特指某一类型的波，它包含了很多波的信息，如直达波、首播、透射波、回折波、折射波等。没有经过反射界面改造，其波形信息中含有该区域内速度异常体的丰富信息，因此本发明利用动态时窗所截取的波场对低速体发育的区域进行高精度速度建模。

在利用时窗对单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据之后，还包括对观测数据进行预处理的步骤，具体包括：利用小波变换对观测数据进行滤波处理，得到处理后的观测数据；对处理后的观测数据进行多次波压制处理。

具体的，利用小波变换对观测数据进行分频去噪，小波变化可无限细分，彼此正交，对含有相干干扰的观测数据进行分频，去噪处理可仅限于很窄的频带中进行，使去噪后对有效波的损失最大限度的减少，不存在傅里叶变换的频泄现象。

采用共中心点叠加法对处理后的观测数据进行多次波压制处理。共中心点叠加法是依据动校正后一次波和多次波之间剩余时差的差异，将不同接收点收到的来自地下同一反射点的不同激发点的信号，经动校正后叠加起来，进行多次波压制。用一次波的速度作动校正，一次波同相轴被校平而多次波仍有剩余时差，通过叠加使一次波得到增强而多次波得到削弱。

另外，还可以采用倾角滤波、速度滤波、扇形滤波等二维滤波法进行多次波压制，滤除多次波保留一次波。

步骤103，获取初始速度模型，并基于初始速度模型对震源波形进行正演，得到正演模拟数据；

具体的,根据观测数据获取地层信息和均方根速度信息；根据地层信息和均方根速度信息计算地层的各反射层的层速度，得到初始速度模型。

对地下进行网格化剖分，模型的大小是m*n，m表示水平方向的网格点数，n表示纵向的网格点数，网格间距为h，初始速度模型的尺寸为水平m*h米，纵向是n*h米。初始速度模型的矩阵形式如公式(3)所示:

声波的波动方程如公式(4)所示:

其中,p为压力场,v_x和v_z分别为横向和纵向速度场；k＝ρv²；t_s为声波传播时间。

在地震勘探领域，声波方程交错网格的思想是把不同的波场分量和地下介质参数配置在不同的网格节点上，时间步进也是采用时间交错步进的。本发明采用的二维声波各物理量和介质参数的网格配置如图3所示，三维声波各物理量和介质参数的网格配置如图4所示。

在声波方程离散前，首先需要采用泰勒展开法推导声波场的空间导数在规则网格和交错网格的高精度近似。

设u(x)有2N+1阶导数，则u(x)在x0+△x,与x0-△x处的2N+1阶泰勒展开式为：

(5)、(6)两式做差得到公式(7)：

同理有：

则任意2N阶精度中心有限差分系数的计算公式为：

令

则有：

化简得：

同理可推导交错网格任意2N阶精度有限差分差分格式和差分系数计算式。

根据泰勒展开式有：

化简得：

其中，a_n为差分系数

因此采用交错网格有限差分，三维一阶速度-应力声波方程可以离散为：

其中，Δx,Δy,Δz,Δt分别为空间和时间采样间隔，p为每个时刻的应力波场，v表示每个时刻的位移波场，x,y,z分别表示不同的方向，ρ表示密度，f表示震源函数。

当给定初值和边值条件时，用上述差分格式可以递推得到各个时刻波场的空间分布。

在有限差分波动方程正演模拟过程中，为了避免数值噪声以及不稳定，对于一个给定的频带宽度，有限差分的网格大小和时间步长需要分别满足频散关系以及稳定性条件。本发明技术方案采用有限差分频散关系需要满足每个最小波长至少需要5个网格点，即为了避免网格色散空间采样间隔需要满足公式(29):

其中Δx为空间网格大小，λ_min是最短波长，

为最小纵波速度，f_max是最大频率。

空间采样网格大小确定之后，还需要选择合适的时间采样大小来满足有限差分数值稳定性条件：

其中Δt为时间采样间隔，

为最大纵波速度。

步骤104，根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；

具体的，将观测数据与正演模拟数据做差，得到波场残差。观测数据与正演模拟数据均为二维的数组，将二维数组内对应的点上的数据相减，得到波场残差。

根据公式(31)构造误差泛函：

其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，d_obs为观测数据，b(m)-d_obs为波场残差，C_D为数据协方差矩阵，C_M为模型的协方差矩阵，m为速度初始模型的模型参数，m_prior为先验信息模型，λ为先验信息比重参数，用来调节模型项和先验信息项的比重。

步骤105，利用波场残差反传播到初始速度模型空间，得到残差反传播数据；

将反传算子作用于波场残差,得到初始速度模型空间的残差反传播数据。

步骤104得到检波点所在位置处的波场残差，将这些点处的波场残差作为震源加载到时间域正演模拟过程中，进行时间反方向传播,得到残差反传播数据。

步骤106，利用正演模拟数据和残差反传播数据计算误差泛函的梯度；

梯度计算是参数反演的关键部分,代表着模型的更新方向，梯度导引类反演方法通过目标泛函对模型参数的导数,来寻找迭代更新方向,以实现对模型的更新。本发明基于伴随状态法进行梯度求取，将正传波场和反传波场的数据残差作为新的震源进行正演,以求取目标函数对模型的梯度。

根据公式(32)计算每个单炮数据的误差泛函的梯度：

其中，m为速度初始模型的模型参数，u为正演模拟播场，B为正演算子，Δd为波场残差，B^-1t为残差反传播数据；将所有单炮数据的梯度进行叠加，得到速度模型空间的全局梯度。

具体的，基于波动方程的特征波全波形反演采用频率域声波波动方程的格林函数进行研究，频率域声波方程的格林函数如公式(33)所示：

选取速度作为模型参数化方式，则常规二范数目标泛函关于速度的梯度为：

其中

k＝ρV_p ²  (37)

其中，P_f为正演模拟数据，P_b残差反传播数据，ω为频率，V_P为速度，k、ρ为初始模型参数，E为误差泛函。

步骤107，利用误差泛函的梯度更新初始速度模型，得到精确速度模型。

具体的，令梯度函数

得到扰动模型，最终的精确速度模型为初始速度模型和扰动模型之和。

根据最快下降法、共轭梯度法、拟牛顿法LBFGS确定最优迭代步长。根据误差泛函的梯度和最优迭代步长对初始速度模型进行更新；确定迭代终止条件；当初始速度模型满足迭代终止条件时，得到精确速度模型。

具体的，利用公式(39)进行模型更新：

m_update＝m_before-α_kd_k  (39)

其中，α_k为第k步的最优迭代步长，d_k为第k步模型的梯度。

迭代终止条件如公式(40)所示：

公式(40)规定了迭代收敛规则，第n次更新的量小于n-1次迭代的一定比例，例如，ε＝0.001，则迭代收敛，反演终止；否则，将更新结果作为输入，进行下一次迭代。

本发明提供的地层声波传播速度模型构建方法，获取探测区域内震源波形的单炮数据，单炮数据包括单炮单道数据和单炮多道数据；利用时窗对单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取初始速度模型，并基于初始速度模型对震源波形进行正演，得到正演模拟数据；根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；利用波场残差反传播到初始速度模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据计算误差泛函的梯度；利用误差泛函的梯度更新初始速度模型，得到精确速度模型。本发明提供的方法，根据地质任务需求利用初至以后一段时间内波场的运动学和动力学特征完成波形反演，实现浅层的高精度速度建模。

专业人员应该还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种地层声波传播速度模型构建方法，其特征在于，所述方法包括：

获取探测区域内震源震动的多个单炮数据，所述单炮数据包括单炮单道数据和单炮多道数据；

利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；

根据所述观测数据获取地层信息和均方根速度信息；

根据所述地层信息和均方根速度信息计算所述地层的各反射层的层速度，得到初始速度模型；

基于所述初始速度模型对震源波形进行正演，得到正演模拟数据；

根据所述观测数据和所述正演模拟数据计算波场残差，并根据所述波场残差构造误差泛函；

利用所述波场残差反传播到所述初始速度模型空间，得到残差反传播数据；

利用所述正演模拟数据和所述残差反传播数据根据伴随状态法计算所述误差泛函的梯度；

利用所述误差泛函的梯度更新所述初始速度模型，得到精确速度模型；

其中，所述基于所述初始速度模型对震源波形进行正演，得到正演模拟数据具体包括:

利用交错网格有限差分法对所述震源波形的波动方程进行时间域离散,得到离散后的波动方程；

根据所述离散后的波动方程和所述初始速度模型确定交错网格各个时刻空间分布的波场值；

所述离散后的波动方程为三维一阶速度-应力声波方程，具体为：

其中，Δx,Δy,Δz,Δt分别为空间和时间采样间隔，p为每个时刻的应力波场，v表示每个时刻的位移波场，x,y,z分别表示不同的方向，ρ表示密度，f表示震源函数，a为差分系数，κ＝ρv²。

2.根据权利要求1所述的方法，其特征在于，所述利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据之后，所述方法还包括：

利用小波变换对所述观测数据进行滤波处理，得到处理后的观测数据；

对所述处理后的观测数据进行多次波压制处理。

3.根据权利要求1所述的方法，其特征在于，所述根据所述观测数据和所述正演模拟数据计算波场残差具体包括：

将所述观测数据与所述正演模拟数据做差，得到所述波场残差。

4.根据权利要求1所述的方法，其特征在于，所述根据所述波场残差构造误差泛函具体包括：

根据公式计算误差泛函，其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，d_obs为观测数据，b(m)-d_obs为波场残差，C_D为数据协方差矩阵，C_M为模型的协方差矩阵，m为速度初始模型的模型参数，m_prior为先验信息模型，λ为先验信息比重参数。

5.根据权利要求1所述的方法，其特征在于，所述利用所述波场残差反传播到所述初始速度模型空间，得到残差反传播数据具体包括：

将反传算子作用于波场残差,得到所述初始速度模型空间的残差反传播数据。

6.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据根据伴随状态法计算所述误差泛函的梯度具体包括：

根据公式计算所述误差泛函的梯度；

其中，k＝ρV_p ²，P_f为正演模拟数据，P_b残差反传播数据，ω为频率，V_P为速度，k、ρ为初始模型参数，E为误差泛函。

7.根据权利要求1所述的方法，其特征在于，所述方法还包括：

根据最快下降法、共轭梯度法、拟牛顿法LBFGS确定最优迭代步长。

8.根据权利要求7所述的方法，其特征在于，所述利用所述误差泛函的梯度更新所述初始速度模型，得到精确速度模型具体包括：

根据误差泛函的梯度和最优迭代步长对初始速度模型进行更新；

确定迭代终止条件；

当所述初始速度模型满足所述迭代终止条件时，得到精确速度模型。

Images