CN105005076A - 基于最小二乘梯度更新速度模型的地震波全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，包括以下步骤:1)获得时间域内地震波场的声波方程；2)构建初始速度模型，设置速度模型更新迭代次数N与允许最小误差值ε；3)构造观测波场数据和计算波场数据的波场误差矢量；4)构造目标函数；5)对目标函数计算获得地震波全波形反演的标准方程；6)引入速度模型更新梯度方向g_k和更新量步长α；7)利用最小二乘方法求解速度模型更新梯度方向g_k；8)对更新量步长α进行插值计算；9)对速度模型进行更新得：m_k＝m_k-1+αg_k；当|αg_k|＜ε时或速度模型更新次数达到速度模型更新迭代次数N时，速度模型更新结束；否则，进入步骤3)。本发明可快速完成速度模型更新，广泛应用于地震波全波形反演方法中。

Description

基于最小二乘梯度更新速度模型的地震波全波形反演方法

技术领域

本发明涉及一种地震波全波形反演方法，特别是关于一种基于最小二乘梯度更新速度模型的地震波全波形反演方法。

背景技术

地震波场反演方法很多，如相位反演成像、振幅反演成像以及全波形反演成像。相位反演和振幅反演利用的是波场运动学信息，反演出的速度分布是一种平滑估计值，它是真实速度模型的低频分量，反映的是大尺度的速度结构。基于波动方程理论的全波形反演综合利用了地震记录中振幅、走时和相位等完整的波场信息，通过拟合，全波形反演通常利用给定时窗内的理论波场记录与实测波场记录差值的二次泛函以及额外的地质约束条件建立优化目标函数，利用迭代线性化反演方法或完全非线性反演法求解符合条件的模型，定量提取地下介质的弹性参数，进而为深部大尺度构造演化分析，为勘探地震成像及速度建模等方面提供可靠依据。

目前，全波形反演的方式大致可以分为时间域全波形反演与频率域全波形反演。频率域全波形反演理论与时间域全波形反演理论有着良好的对应关系。两种反演方法的观测数据表现形式不同造成正演波场计算方式的差别，而反演思想完全一致，与时间域相比，频率域波场是相互解耦的。不同频率数据对异常体反映能力不一，根据需要既可以利用部分频段数据反演，也可以使用全频段数据同时进行反演。

全波形反演可利用逆时偏移产生的成像结果近似来更新速度模型，基于逆时偏移的全波形反演获得低频成像不再需要炮集的低频信息且具有很高的保真度。但是，反演使用大量低频的信息来更新背景速度，而偏移使用高频信息来勾画地下构造的比较精确的边界，由于反演和偏移的目标在不同的频带，因此，在实际处理中，由于噪声和模型假象，高保真度很难达到，且常规全波形反演流程中，每次的全波形反演的非线性迭代过程中，只利用一次的逆时偏移的成像结果作为速度模型更新梯度的方向，这样容易导致更新梯度计算的更新量不足，且收敛速度较慢。

发明内容

针对上述问题，本发明的目的是提供一种可有效提高速度模型更新梯度的准确性且可快速完成速度模型更新的基于最小二乘梯度更新速度模型的地震波全波形反演方法。

为实现上述目的，本发明采取以下技术方案：一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，包括以下步骤：

1)在时间域内，地震波场的声波方程为：

$M\left(x\right)\frac{d^{2}u(x,t)}{d{t}^2}=A\left(x\right)u(x,t)+s(x,t);$

其中，u是波场向量，s是震源向量，x和t分别是空间和时间坐标参数，M为质量矩阵，A为刚度矩阵；

2)构建初始速度模型，设置速度模型更新迭代次数N与速度模型允许的最小误差值ε；

3)构造观测波场数据和计算波场数据的波场误差矢量Δd；

4)利用最小二乘流程和波场误差矢量Δd，构造目标函数C(m_k)，目标函数C(m_k)的表达式为：

C(m_k)＝1/2Δd^TΔd；

式中，T表示转置；

m_k为第k次更新的速度模型，更新公式为：

m_k＝m_k-1+Δm；

式中，m_k-1为第k-1次更新的速度模型；Δm为扰动模型；

5)对目标函数进行计算获得地震波全波形反演的标准方程为：

$\mathrm{\Δm}=-{\left[\frac{{\∂}^{2}C\left(m_{k-1}\right)}{\∂{m_n}^2}\right]}^{-1}\frac{\∂C\left(m_{k-1}\right)}{\∂m_n};$

式中，m_n为模型参数；

6)引入速度模型更新梯度方向和更新量步长，扰动模型Δm表示为：

Δm＝αg_k；

式中，g_k为速度模型更新梯度方向，α为更新量步长；

7)利用最小二乘的方法求解地震波全波形反演的速度模型更新梯度方向g_k，其公式为：

$\left[\frac{{\∂}^{2}C\left(m_{k-1}\right)}{\∂{m_n}^2}\right]g_k=-\frac{\∂C\left(m_{k-1}\right)}{\∂m_n};$

8)在速度模型更新梯度方向g_k上对更新量步长α进行插值计算求取更新量步长α；

9)对速度模型进行更新得：

m_k＝m_k-1+αg_k；

当|αg_k|＜ε时或速度模型更新次数达到速度模型更新迭代次数N时，速度模型更新结束；否则，进入步骤3)。

所述步骤3)中，波场误差矢量Δd的计算公式如下：

Δd＝d_obs-d_cal；

式中，d_obs为地震观测波场数据，d_cal为在迭代过程中的计算波场数据；

计算波场数据d_cal的计算公式如下：

式中，为计算波场数据d_cal和波场向量u之间的关系函数。

所述步骤5)中，地震波全波形反演的标准方程具体计算过程包括以下步骤：

(1)对目标函数C(m_k)进行二阶泰勒-拉格朗日展开，得到以下形式：

$C(m_{k-1}+\mathrm{\Δm})=C\left(m_{k-1}\right)+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\frac{\∂C\left(m_{k-1}\right)}{\∂m_j}\mathrm{\Δ}{m}_j+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{l=1}{\overset{M}{\mathrm{\Σ}}}\frac{{\∂}^{2}C\left(m_{l-1}\right)}{\∂m_j\∂m_l}\mathrm{\Δ}{m}_j\mathrm{\Δ}{m}_l+o^2\left(m_{k-1}\right);$

式中，j、l分别为二维数据的位置参数；M为正整数；m_j、m_l分别为更新模型变量m_k-1在二维参数模型中某一个方向的分量；

(2)对步骤(1)中目标函数C(m_k-1+Δm)的展开式以模型参数m_n为变量进行求导：

$\frac{\∂C\left(m_k\right)}{\∂m_n}=\frac{\∂C\left(m_{k-1}\right)}{\∂m_n}+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\frac{{\∂}^{2}C\left(m_{k-1}\right)}{\∂m_j\∂m_n}\mathrm{\Δ}{m}_j;$

当目标函数C(m_k)的导数为零时，目标函数C(m_k)得到极值，此时，扰动模型Δm的表达式为：

$\mathrm{\Δm}=-{\left[\frac{{\∂}^{2}C\left(m_{k-1}\right)}{\∂{m_n}^2}\right]}^{-1}\frac{\∂C\left(m_{k-1}\right)}{\∂m_n};$

扰动模型Δm的表达式即为地震波全波形反演的标准方程。

所述步骤8)中，更新量步长α满足伍尔夫条件：

$C(m_k+\mathrm{\&alpha;}{g}_k)<C\left(m_{k-1}\right)+\mathrm{\&alpha;}{c}_1{g}_k^T\&dtri;C\left(m_k\right);$

$|d^T\&dtri;C(m_k+\mathrm{\&alpha;}{g}_k)/g_k^T\&dtri;C\left(m_{k-1}\right)|<c_2;$

式中，c₁和c₂均为固定参数，其中c₁∈(0，1)，c₂∈(c₁，1)，▽C为目标函数C的梯度。

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于采用迭代非线性反演方法求解符合条件的模型，为深部大尺度构造演化分析、勘探地震成像及速度建模等方面提供可靠依据。2、本发明由于采用在传统全波形反演的计算流程基础上，在每次的全波形反演的非线性迭代过程中，引入最小二乘逆时偏移的流程，通过最小二乘的约束和多次迭代修正，有效地提高了速度模型更新梯度准确性，可快速完成速度模型的更新。综上所述，本发明可广泛应用于地震波全波形反演方法中。

附图说明

图1是本发明的流程示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提供一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，包括以下步骤：

1)在时间域内，地震波场的声波方程为：

$M\left(x\right)\frac{d^{2}u(x,t)}{d{t}^2}=A\left(x\right)u(x,t)+s(x,t)---\left(1\right)$

其中，u是波场向量，s是震源向量，x和t分别是空间和时间坐标参数，M为质量矩阵，A为刚度矩阵。

2)构建初始速度模型，设置速度模型更新迭代次数N与速度模型允许的最小误差值ε。

3)构造观测波场数据和计算波场数据的波场误差矢量Δd，波场误差矢量Δd的计算公式如下：

Δd＝d_obs-d_cal      (2)

式中，d_obs为地震观测波场数据，d_cal为在迭代过程中的计算波场数据。

其中，根据波场向量u计算波场数据d_cal：

式中，为计算波场数据d_cal和波场向量u之间的关系函数。

4)利用最小二乘流程和波场误差矢量Δd构造目标函数C(m_k)如下：

C(m_k)＝1/2Δd^TΔd      (4)

式中，m_k为第k次更新的速度模型，T表示转置。

速度模型的求取过程是对步骤2)中构建的初始速度模型m₀的不断更新过程，第k次更新的速度模型m_k表示为第k-1次更新的速度模型m_k-1与扰动模型Δm的和：

m_k＝m_k-1+Δm      (5)

5)对目标函数C(m_k)进行计算获得地震波全波形反演的标准方程：

$\mathrm{\&Delta;m}=-{\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]}^{-1}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n}---\left(6\right)$

式中，m_n为模型参数。

具体计算过程包括以下步骤：

(1)对目标函数C(m_k)进行二阶Taylo-Lagrange(泰勒-拉格朗日)展开，得到以下形式：

$C(m_{k-1}+\mathrm{\&Delta;m})=C\left(m_{k-1}\right)+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_j}\mathrm{\&Delta;}{m}_j+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;}^{2}C\left(m_{l-1}\right)}{\&PartialD;m_j\&PartialD;m_l}\mathrm{\&Delta;}{m}_j\mathrm{\&Delta;}{m}_l+o^2\left(m_{k-1}\right)---\left(7\right)$

式中，j、l分别为二维数据的位置参数；M为正整数；m_j、m_l分别为更新模型变量m_k-1在二维参数模型中某一个方向的分量。

(2)对步骤(1)中目标函数C(m_k-1+Δm)的展开式以模型参数m_n为变量进行求导：

$\frac{\&PartialD;C\left(m_k\right)}{\&PartialD;m_n}=\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n}+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;m_j\&PartialD;m_n}\mathrm{\&Delta;}{m}_j---\left(8\right)$

当目标函数C(m_k)的导数为零时，目标函数C(m_k)得到极值，此时，扰动模型Δm的表达式为：

$\mathrm{\&Delta;m}=-{\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]}^{-1}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n}---\left(6\right)$

公式(6)为地震波全波形反演的标准方程。

6)引入速度模型更新梯度方向和更新量步长，扰动模型Δm表示为：

Δm＝αg_k        (9)

式中，g_k为速度模型更新梯度方向，α为更新量步长。

7)利用最小二乘的方法求解地震波全波形反演的速度模型更新梯度方向g_k，其公式为：

$\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]g_k=-\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n}---\left(10\right)$

8)在速度模型更新梯度方向g_k上对更新量步长α进行插值计算求取更新量步长α。其中，更新量步长α满足伍尔夫(Wolfe)条件：

$C(m_k+\mathrm{\&alpha;}{g}_k)<C\left(m_{k-1}\right)+\mathrm{\&alpha;}{c}_1{g}_k^T\&dtri;C\left(m_k\right)---\left(11\right)$

$|d^T\&dtri;C(m_k+\mathrm{\&alpha;}{g}_k)/g_k^T\&dtri;C\left(m_{k-1}\right)|<c_2---\left(12\right)$

式中，c₁和c₂均为固定参数，其中c₁∈(0，1)，c₂∈(c₁，1)，▽C为目标函数C的梯度；

9)对速度模型进行更新得：

m_k＝m_k-1+αg_k       (13)

当|αg_k|＜ε时或速度模型更新次数达到速度模型更新迭代次数N时，速度模型更新结束；否则，进入步骤3)。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，流程中涉及各步骤的数据表达方式、数据的连接方式和程序实现方式等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (5)

1.一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，包括以下步骤：

1)在时间域内，地震波场的声波方程为：

$M\left(x\right)\frac{d^{2}u(x,t)}{d{t}^2}=A\left(x\right)u(x,t)+s(x,t);$

其中，u是波场向量，s是震源向量，x和t分别是空间和时间坐标参数，M为质量矩阵，A为刚度矩阵；

2)构建初始速度模型，设置速度模型更新迭代次数N与速度模型允许的最小误差值ε；

3)构造观测波场数据和计算波场数据的波场误差矢量Δd；

4)利用最小二乘流程和波场误差矢量Δd，构造目标函数C(m_k)，目标函数C(m_k)的表达式为：

C(m_k)＝1/2Δd^TΔd；

式中，T表示转置；

m_k为第k次更新的速度模型，更新公式为：

m_k＝m_k-1+Δm；

式中，m_k-1为第k-1次更新的速度模型；Δm为扰动模型；

5)对目标函数进行计算获得地震波全波形反演的标准方程为：

$\mathrm{\&Delta;m}=-{\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]}^{-1}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n};$

式中，m_n为模型参数；

6)引入速度模型更新梯度方向和更新量步长，扰动模型Δm表示为：

Δm＝αg_k；

式中，g_k为速度模型更新梯度方向，α为更新量步长；

7)利用最小二乘的方法求解地震波全波形反演的速度模型更新梯度方向g_k，其公式为：

$\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]g_k=-\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n};$

8)在速度模型更新梯度方向g_k上对更新量步长α进行插值计算求取更新量步长α；

9)对速度模型进行更新得：

m_k＝m_k-1+αg_k；

当|αg_k|＜ε时或速度模型更新次数达到速度模型更新迭代次数N时，速度模型更新结束；否则，进入步骤3)。

2.如权利要求1所述的一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，其特征在于：所述步骤3)中，波场误差矢量Δd的计算公式如下：

Δd＝d_obs-d_cal；

式中，d_obs为地震观测波场数据，d_cal为在迭代过程中的计算波场数据；

计算波场数据d_cal的计算公式如下：

式中，为计算波场数据d_cal和波场向量u之间的关系函数。

3.如权利要求1或2所述的一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，其特征在于：所述步骤5)中，地震波全波形反演的标准方程具体计算过程包括以下步骤：

(1)对目标函数C(m_k)进行二阶泰勒-拉格朗日展开，得到以下形式：

$C(m_{k-1}+\mathrm{\&Delta;m})=C\left(m_{k-1}\right)+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_j}\mathrm{\&Delta;}{m}_j+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;}^{2}C\left(m_{l-1}\right)}{\&PartialD;m_j\&PartialD;m_l}\mathrm{\&Delta;}{m}_j\mathrm{\&Delta;}{m}_l+o^2\left(m_{k-1}\right);$

式中，j、l分别为二维数据的位置参数；M为正整数；m_j、m_l分别为更新模型变量m_k-1在二维参数模型中某一个方向的分量；

(2)对步骤(1)中目标函数C(m_k-1+Δm)的展开式以模型参数m_n为变量进行求导：

$\frac{\&PartialD;C\left(m_k\right)}{\&PartialD;m_n}=\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n}+\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;m_j\&PartialD;m_n}\mathrm{\&Delta;}{m}_j;$

当目标函数C(m_k)的导数为零时，目标函数C(m_k)得到极值，此时，扰动模型Δm的表达式为：

$\mathrm{\&Delta;m}=-{\left[\frac{{\&PartialD;}^{2}C\left(m_{k-1}\right)}{\&PartialD;{m_n}^2}\right]}^{-1}\frac{\&PartialD;C\left(m_{k-1}\right)}{\&PartialD;m_n};$

扰动模型Δm的表达式即为地震波全波形反演的标准方程。

4.如权利要求1或2所述的一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，其特征在于：所述步骤8)中，更新量步长α满足伍尔夫条件：

$C(m_k+\mathrm{\&alpha;}{g}_k)<C\left(m_{k-1}\right)+\mathrm{\&alpha;}{c}_1{g}_k^T\&dtri;C\left(m_k\right);$

$|d^T\&dtri;C(m_k+\mathrm{\&alpha;}{g}_k)/g_k^T\&dtri;C\left(m_{k-1}\right)|<c_2;$

式中，c₁和c₂均为固定参数，其中c₁∈(0，1)，c₂∈(c₁，1)，▽C为目标函数C的梯度。

5.如权利要求3所述的一种基于最小二乘梯度更新速度模型的地震波全波形反演方法，其特征在于：所述步骤8)中，更新量步长α满足伍尔夫条件：

$C(m_k+\mathrm{\&alpha;}{g}_k)<C\left(m_{k-1}\right)+\mathrm{\&alpha;}{c}_1{g}_k^T\&dtri;C\left(m_k\right);$

$|d^T\&dtri;C(m_k+\mathrm{\&alpha;}{g}_k)/g_k^T\&dtri;C\left(m_{k-1}\right)|<c_2;$

式中，c₁和c₂均为固定参数，其中c₁∈(0，1)，c₂∈(c₁，1)，▽C为目标函数C的梯度。