CN109085643A - 早至波的分步联合反演方法 - Google Patents

Abstract

本发明提供一种早至波的分步联合反演方法，包括：输入观测数据和初始速度模型；定义走时层析目标泛函表达式，利用伴随状态法获得目标泛函关于模型的梯度表达式；利用获得的走时层析梯度对模型背景场进行更新；利用获得准确的模型背景场为全波形反演的初始模型；利用得到的速度模型，建立全波形反演的目标泛函，利用伴随状态法获得目标泛函关于模型的梯度表达式；利用获得的全波形反演梯度对模型背景场进行更新；利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果。该早至波的分步联合反演方法获得的速度场无论是在浅层还是深层，都能够为成像提供更为准确的速度模型，验证了联合反演在速度场精度和反演稳定性方面的优势。

Description

早至波的分步联合反演方法

技术领域

本发明涉及地震资料处理技术领域，特别是涉及到一种早至波的分步联合反演方法。

背景技术

初至波层析反演在目前来说是一种近地表成像的标准方法，但对于复杂地质情况比如说隐藏的速度层区不能很好的解决；早至波波形反演是一个稳健的工具对于解决复杂的结构反演。然而因为波形反演的多解性在实际中波形反演的结果并不能与初至时间相吻合，并且波形反演要求要有一个较好的初始模型，因而许多学者提出可初至波层析与波形联合反演算法。

联合反演考虑的是解决波形反演局部极小的问题(Haber and Oldenburg,1997；Juliàet al.,2000；Gallardo and Meju,2003；Colombo and De Stefano,2007)，MAYHALED于2013年提出了使用一种混合的全波形反演，使用波动方程反射波旅行时反演来更新低波数分量，使用FWI来更新高波数分量；Zhang等人于2014年提出了初至波走时与波形联合反演的算法进行近地表结构体的成像；该方法建立了走时层析与全波形反演的基础；Jiao K，Sun D等于2015年提出了一种适应性的FWI算法，使用通过建立走时信息与模型的关系，使用FWI来修正背景速度场避免了陷入局部极值；Sun Mengyao等人提出了一种层析与全波形反演交替算法，该算法减小了计算内存需要，通过对实际资料的测试表明该算法能够明显提高反演结果的效果。但是多解性仍然是旅行时反演和FWI反演中的固有问题，联合反演并不能解决这些问题，约束反演结果。为此我们发明了一种新的早至波的分步联合反演方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种融合走时层析与全波形建模方法的优点，适应陆上早至波信号的特点，获得具有较高精度的全尺度(浅中深)速度模型，破解生产中速度模型建立精度低的难题的早至波的分步联合反演方法。

本发明的目的可通过如下技术措施来实现：早至波的分步联合反演方法，该早至波的分步联合反演方法包括：步骤1：输入观测数据和初始速度模型；步骤2：定义走时层析目标泛函表达式，利用伴随状态法获得目标泛函关于模型的梯度表达式；步骤3：利用获得的走时层析梯度对模型背景场进行更新；步骤4：利用获得准确的模型背景场为全波形反演的初始模型；步骤5：利用步骤4中得到的速度模型，建立全波形反演的目标泛函，利用伴随状态法获得目标泛函关于模型的梯度表达式；步骤6：利用获得的全波形反演梯度对模型背景场进行更新；步骤7：利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果。

本发明的目的还可通过如下技术措施来实现：

在步骤1中，利用快速扫描算法，对简单双层模型和复杂实际模型进行试算，证明快速扫描算法计算精度和效率。

在步骤1中，输入要观测的原始单炮数据，以走时层析反演的到的速度模型作为初始速度模型m₀，对观测数据做切除，获得观测数据的早至波，同时基于初始速度模型通过正演获得模拟的地震数据，同样对模拟地震数据做切除。

在步骤2中，利用快速扫描法计算原始单炮数据的走时；获取目标泛函对观测数据与计算数据的误差进行测量，并且将射线类的解大型的稀疏的走时层析方程组转化为非线性最优化问题，所述基于稀疏台站数据的走时层析的目标泛函为：

其中：E为误差函数；T_obs为观测数据，T为模拟数据残差，Ω指的是整个模型空间，指的是模型的边界也将就是检波点所在的位置，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置；模拟数据根据早至波的梯度相应范围分析选取相应的偏移距和时窗使其与观测数据匹配；利用共轭梯度法更新模型，当误差泛函值小于允许值时，更行停止；使用非线性共轭梯度方法求解目标泛函。

在步骤2中，利用伴随状态法更新梯度，利用抛物线插值法计算步长，目标泛函的梯度表示为：

其中：E为误差函数，Ω指的是整个模型空间，λ为梯度，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置。

在步骤3中，利用步骤2中获得的梯度，对梯度做预条件处理得到d^k，通过求得的d^k计算步长α^k，利用公式m^k＝m^k-1-α^kd^k，更新速度模型，如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

在步骤5中，全波形反演的目标泛函为：

上式中，u(t,X,X_S；v)的意义是模拟得的全波场，d_ObS(t,X_r,X_S)的意义是测量得的炮记录，X_r的意义是与炮点相关的检波点的坐标，X_S的意义是炮点的坐标，v的意义是速度，R的意义是限定检波器位置限定算子，Ru(t,X,X_S；v)＝u(t,X＝X_r,X_S；v)。

在步骤5中，利用伴随状态法得到目标函数对速度模型的梯度，g(v)的表达式如下：

上式中，g(v)是得到的梯度，u是正传波场值，R指的是检波点位置，L^*指的是波场反传算子，式5表达的意思是梯度等于正传波场关于时间的二阶导数乘以检波点位置的残差的反传波场，残差就是Ru-d_ObS。

在步骤5中，为了进一步得到共轭梯度方向，引入共轭梯度方向修正因子β_k，这样得到优化后的方向如下：

β_k是指共轭梯度里面的参数，式6是共轭梯度算法，g_k是梯度，k-1,k分别代表迭代之前与迭代之后；

全波形反演的目标泛函的梯度表达式为：

式中：m为模型参数；v为背景速度；u为正传波场，u^*为残差反传波场。

在步骤6中，获得模型数据的早至波，早至波全波形反演更新梯度对梯度做预条件处理得到d^k，通过以上获得计算步长α^k，利用公式m^k＝m^k-1-α^kd^k，更新速度模型；如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

在步骤6中，利用抛物线插值法计算步长：

α指的是反演用的步长，E₀、E₁、E₂是同一个梯度由不同试验步长α0、α1、α2根据速度模型更新公式(式9)计算出来的目标泛函值；

更新速度模型：

m_k+1＝m_k+α_kd_k(式9)

m_k、m_k+1都表示模型，也就是进行正演使用的速度模型，k、k+1代表迭代之前与迭代之后，d_k是对梯度做预条件处理得到的。

在步骤7中，利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果对真实模型、初始模型、走时层析反演模型以及联合反演模型的不同方位进行了显示对比。

本发明中的早至波的分步联合反演方法，属于地震资料的反演处理，具体是一种可得到高精度的速度场建模结果。首先利用走时层析方法快速地为全波形反演提供相对准确的初始模型，一定程度上避免了周波跳跃问题。其次，全波形反演在走时层析基础上进一步细化模型，揭示复杂和细小构造。最后，将两者加权结合，构建统一的反演框架，形成一套稳定的自动化反演策略。本发明融合两种建模方法的优点，适应陆上早至波信号的特点，获得具有较高精度的全尺度速度模型，破解生产中速度模型建立精度低的难题。本发明的走时层析能够快速地为全波形反演提供相对准确的初始模型，一定程度上避免周波跳跃问题。另一方面，全波形反演能够在走时层析基础上进一步细化模型，揭示复杂和细小构造。最后，将两者加权结合，构建统一的反演框架，能够形成一套稳定的自动化反演策略。利用不同反演结果所得到的速度场进行RTM成像结果进行分析，采用联合反演所获得的速度场无论是在浅层还是深层，都能够为成像提供更为准确的速度模型，验证了联合反演在速度场精度和反演稳定性方面的优势。

附图说明

图1为本发明的走时层析反演方法的一具体实施例的流程图；

图2为本发明的一具体实施例中快速扫描算法计算走时等值线图；

图3为走时层析反演的不同偏移距的目标泛函曲线图；

图4为三点抛物线插值方法求取步长的示意图；

图5为三维Marmousi模型的示意图；

图6为做走时层析反演的初始模型,输入的最初模型的示意图；

图7为走时层析反演得到的模型的示意图；

图8为本发明的全波形反演方法的一具体实施例的流程图；

图9为全波形反演的不同偏移距目标泛函曲线图；

图10为以图7中的走时层析反演模型作为初始模型进行全波形反演得到的联合反演得到的模型的示意图；

图11为图5中真实模型的不同方位显示的示意图；

图12为图6中初始模型的不同方位显示的示意图；

图13为图7中走时层析反演模型的不同方位显示的示意图；

图14为图10中联合反演模型的不同方位显示的示意图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

(1)走时层析反演

图1为本发明的走时层析反演方法的流程图。

首先输入要观测的数据(原始单炮数据)，确定初始速度模型m₀(如图7所示)，利用快速扫描法计算原始数据与模型数据走时(图2)，获取不同偏移距的目标泛函(图3)，然后利用伴随状态法计算梯度，对梯度做预条件处理得到d^k，通过求得的d^k计算步长α^k，利用公式m^k＝m^k-1-α^kd^k，更新速度模型，如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

1)输入观测数据(原始单炮数据)和初始速度模型；需要给定一个初始速度模型如图6所示；利用快速扫描算法，对简单双层模型和复杂实际模型进行了试算，证明了快速扫描算法计算精度和效率，如图2所示。

2)定义走时层析目标泛函表达式，利用伴随状态法获得目标泛函关于模型的梯度表达式。利用快速扫描法计算原始单炮数据的走时；

快速扫描法的主要思想是基于因果关系将走时场传播的方向分成有限个组，对于每一组分别利用高斯—赛德尔迭代方法求解非线性逆风差分格式离散化后的方程组。每一次迭代也称为一次扫描，每次扫描(高斯—赛德尔迭代)按一定的方向求解沿该方向传播的走时场。以二维声波方程为例，扫描方向可以分为四组：右上、左上、左下和右下。初值在这四个方向上依次传播，可以通过四组不同顺序的高斯—赛德尔迭代：

式1中的I、J分别代表模型中扫描的X方向和Z方向。

依据步骤2制定了走时层析的目标泛函表达式，并绘制了不同偏移距的走时层析的目标泛函，如图3所示。获取目标泛函对观测数据与计算数据的误差进行测量，并且将射线类的解大型的稀疏的走时层析方程组转化为非线性最优化问题。所述基于稀疏台站数据的走时层析的目标泛函为：

其中：E为误差函数；T_obs为观测数据，T为模拟数据残差，Ω指的是整个模型空间，指的是模型的边界也将就是检波点所在的位置，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置。在此，模拟数据根据早至波的梯度相应范围分析选取相应的偏移距和时窗使其与观测数据匹配。利用共轭梯度法更新模型，当误差泛函值小于允许值时，更行停止；使用非线性共轭梯度方法求解目标泛函。

3)利用获得的走时层析梯度对模型背景场进行更新。在步骤3中利用伴随状态法更新梯度，利用抛物线插值法计算步长，如图4所示。步骤3中目标泛函的梯度表示为：

其中：E为误差函数，Ω指的是整个模型空间，λ为梯度，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置。

步骤3中得到更新模型作为步骤4中的初始模型。

4)最后利用获得准确的模型背景场为全波形反演的初始模型。进入步骤4更新速度模型，得到最终的走时层析反演模型，如图7所示。

(2)全波形反演

图8为本发明的早至波的全波形反演方法的流程图；

首先输入要观测的数据(原始单炮数据)，以走时层析反演的到的速度模型作为初始速度模型m₀(如图6所示)，对观测数据做切除，获得观测数据的早至波，同时基于初始速度模型通过正演获得模拟的地震数据，同样对模拟地震数做切除。获得模型数据的早至波，早至波全波形反演(EFWI)更新梯度对梯度做预条件处理得到d^k，通过以上获得计算步长α^k，利用公式m^k＝m^k-1-α^kd^k，更新速度模型，如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

1)利用走时层析反演中步骤4得到的速度模型(图7)，建立全波形反演的目标泛函，利用伴随状态法获得目标泛函关于模型的梯度表达式。依据步骤4中得到的反演模型，作为全波形反演的初始模型，建立全波形反演的目标泛函，并绘制出目标泛函曲线(图9)；伴随状态法求解目标泛函，更新梯度并用抛物线插值法计算步长，更新速度模型。全波形反演的目标泛函为：

上式中，u(t,X,X_S；v)的意义是模拟得的全波场，d_ObS(t,X_r,X_S)的意义是测量得的炮记录，X_r的意义是与炮点相关的检波点的坐标，X_S的意义是炮点的坐标，v的意义是速度，R的意义是限定检波器位置限定算子(因为实际不是在整个模型空间观测数据，而是在检波点处)，Ru(t,X,X_S；v)＝u(t,X＝X_r,X_S；v)。

利用伴随状态法得到目标函数对速度模型的梯度，g(v)的表达式如下：

上式中，g(v)是得到的梯度，u是正传波场值，R指的是检波点位置，L^*指的是波场反传算子，式5表达的意思是梯度等于正传波场关于时间的二阶导数乘以检波点位置的残差的反传波场，残差就是Ru-d_ObS。

为了进一步得到共轭梯度方向，引入共轭梯度方向修正因子β_k，这样得到优化后的方向如下：

β_k是指共轭梯度里面的参数，式6是共轭梯度算法，g_k是梯度，k-1,k分别代表迭代之前与迭代之后；

全波形反演的目标泛函的梯度表达式为：

式中：m为模型参数；v为背景速度；u为正传波场，u^*为残差反传波场。

2)利用获得的全波形反演梯度对模型背景场进行更新当误差泛函值小于允许值时，则停止更新模型，输出最终反演模型，如图10所示；

利用抛物线插值法计算步长：

α指的是反演用的步长，E₀、E₁、E₂是同一个梯度由不同试验步长α0、α1、α2根据速度模型更新公式(式9)计算出来的目标泛函值；

更新速度模型：

m_k+1＝m_k+α_kd_k(式9)

m_k、m_k+1都表示模型，也就是进行正演使用的速度模型，k、k+1代表迭代之前与迭代之后，d_k是对梯度做预条件处理得到的。

3)最后利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果对三维Marmousi模型(图5)、初始模型(图6)、走时层析反演模型(图7)以及联合反演模型(图10)的不同方位进行了显示对比，如图11-图14；其中图5中的Marmousi模型是由法国石油研究院所属协会于1988年做出的模型，它是基于非洲西南部安哥拉的地质模型而得出的，是地震物理科研中经常可以用到的标准模型之一，在本实施方式中，利用走时层析与全波形反演分步联合反演的到的模型(图10)与真实的Marmousi模型(图5)作对比。

Claims (12)

1.早至波的分步联合反演方法，其特征在于，该早至波的分步联合反演方法包括：

步骤1：输入观测数据和初始速度模型；

步骤2：定义走时层析目标泛函表达式，利用伴随状态法获得目标泛函关于模型的梯度表达式；

步骤3：利用获得的走时层析梯度对模型背景场进行更新；

步骤4：将获得的模型背景场作为全波形反演的初始模型；

步骤5：利用步骤4中得到的速度模型，建立全波形反演的目标泛函，利用伴随状态法获得目标泛函关于模型的梯度表达式；

步骤6：利用获得的全波形反演梯度对模型背景场进行更新；

步骤7：利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果。

2.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤1中，利用快速扫描算法，对简单双层模型和复杂实际模型进行试算，证明快速扫描算法计算精度和效率。

3.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤1中，输入要观测的原始单炮数据，以走时层析反演的到的速度模型作为初始速度模型m₀，对观测数据做切除，获得观测数据的早至波，同时基于初始速度模型通过正演获得模拟的地震数据，同样对模拟地震数据做切除。

4.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤2中，利用快速扫描法计算原始单炮数据的走时；获取目标泛函对观测数据与计算数据的误差进行测量，并且将射线类的解大型的稀疏的走时层析方程组转化为非线性最优化问题，所述基于稀疏台站数据的走时层析的目标泛函为：

其中：E为误差函数；T_obs为观测数据，T为模拟数据残差，Ω指的是整个模型空间，指的是模型的边界也将就是检波点所在的位置，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置；模拟数据根据早至波的梯度相应范围分析选取相应的偏移距和时窗使其与观测数据匹配；利用共轭梯度法更新模型，当误差泛函值小于允许值时，更行停止；使用非线性共轭梯度方法求解目标泛函。

5.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤2中，利用伴随状态法更新梯度，利用抛物线插值法计算步长，目标泛函的梯度表示为：

其中：E为误差函数，Ω指的是整个模型空间，λ为梯度，c是模型速度，T_obs是观测到的走时，x是向量(x,z)点在模型中的位置。

6.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤3中，利用步骤2中获得的梯度，对梯度做预条件处理得到d^k，通过求得的d^k计算步长α^k，利用公式m^k＝m^{k -1}-α^kd^k，更新速度模型，如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

7.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤5中，全波形反演的目标泛函为：

上式中，u(t,X,X_S；v)的意义是模拟得的全波场，d_ObS(t,X_r,X_S)的意义是测量得的炮记录，X_r的意义是与炮点相关的检波点的坐标，X_S的意义是炮点的坐标，v的意义是速度，R的意义是限定检波器位置限定算子，Ru(t,X,X_S；v)＝u(t,X＝X_r,X_S；v)。

8.根据权利要求7所述的早至波的分步联合反演方法，其特征在于，在步骤5中，利用伴随状态法得到目标函数对速度模型的梯度，g(v)的表达式如下：

上式中，g(v)是得到的梯度，u是正传波场值，R指的是检波点位置，L^*指的是波场反传算子，式5表达的意思是梯度等于正传波场关于时间的二阶导数乘以检波点位置的残差的反传波场，残差就是Ru-d_ObS。

9.根据权利要求8所述的早至波的分步联合反演方法，其特征在于，在步骤5中，为了进一步得到共轭梯度方向，引入共轭梯度方向修正因子β_k，这样得到优化后的方向如下：

β_k是指共轭梯度里面的参数，式6是共轭梯度算法，g_k是梯度，k-1,k分别代表迭代之前与迭代之后；

全波形反演的目标泛函的梯度表达式为：

式中：m为模型参数；v为背景速度；u为正传波场，u^*为残差反传波场。

10.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤6中，获得模型数据的早至波，早至波全波形反演更新梯度对梯度做预条件处理得到d^k，通过以上获得计算步长α^k，利用公式m^k＝m^k-1-α^kd^k，更新速度模型；如果得到的最新速度模型m^k满足：(上一次迭代得到的目标泛函值-当前迭代的目标泛函值)/上一次迭代得到的目标泛函值小于0.000001，则输出m^k最为最终速度模型，否则k＝k+1，循环计算。

11.根据权利要求10所述的早至波的分步联合反演方法，其特征在于，在步骤6中，利用抛物线插值法计算步长：

α指的是反演用的步长，E₀、E₁、E₂是同一个梯度由不同试验步长α0、α1、α2根据速度模型更新公式(式9)计算出来的目标泛函值；

更新速度模型：

m_k+1＝m_k+α_kd_k (式9)

m_k、m_k+1都表示模型，也就是进行正演使用的速度模型，k、k+1代表迭代之前与迭代之后，d_k是对梯度做预条件处理得到的。

12.根据权利要求1所述的早至波的分步联合反演方法，其特征在于，在步骤7中，利用获得准确的模型背景场进行全波形反演，获得高精度的速度场建模结果对真实模型、初始模型、走时层析反演模型以及联合反演模型的不同方位进行了显示对比。