CN109188519B - 一种极坐标下的弹性波纵横波速度反演系统及方法 - Google Patents

Abstract

本发明提出了一种极坐标下的弹性波纵横波速度反演系统及方法，属于石油勘探领域，该方法包括如下步骤：输入初始纵横波速度模型、多分量地震记录；极坐标系下的网格剖分；将纵横波速度模型变换到极坐标系下；求取极坐标系下的纵横波梯度公式；求取加权速度更新步长；更新极坐标系下的纵横波速度，判断是否满足误差条件；将反演的纵横波速度变换到笛卡尔坐标系下；输出反演的纵横波速度。本发明能够很好地解决一些特殊地质环境中，如在井孔环境、隧道空间等的速度反演难题；同时将多分量方法引入到本发明中，能够同时对纵横波速度进行反演，得到更真实、准确的纵横波速度，为井孔环境、隧道空间等特殊地质构造的成像及解释提供基础。

Description

一种极坐标下的弹性波纵横波速度反演系统及方法

技术领域

本发明属于石油勘探领域，具体涉及一种极坐标下的弹性波纵横波速度反演系统及方法。

背景技术

反射地震勘探以提取并应用地震波所携带的信息获取高精度地下介质参数分布为主要目标。随着地震勘探程度的加深，需要获取精度更高的地下介质参数分布，常规的地震数据处理方法由于理论上的限制已不能完全满足地震勘探的需求。波形反演是一种建立在完备理论基础上的反演方法，以正演数据与观测数据的波形匹配为准则通过非线性反演获取地下介质参数分布。由于充分利用了数据中的波形信息(包含旅行时、振幅等)并将完整的波动方程作为正演算子，波形反演在理论上具有更高分辨率。此外，基于弹性介质假设而非声介质假设能够更真实的描述地震波在实际介质中的传播，同时为多参数反演提供了前提。在理论上，弹性介质假设下的波形反演方法能够在一定程度上弥补常规地震数据处理方法的不足，提供精度更高的地下介质参数估计，为油气藏的精细描述提供更多的地球物理证据。

在一些特殊的环境(如井孔环境、隧道空间等)中，传统的速度反演方法无法对以井孔为辐射中心的速度进行准确的反演。

发明内容

基于上述技术问题，本发明提供一种极坐标下的弹性波纵横波速度反演系统及方法。

本发明所采用的技术解决方案是：

一种极坐标下的弹性波纵横波速度反演系统，包括顺序连接的输入模块、网格生成模块、坐标变换模块、梯度求取模块、步长求取模块、速度更新模块、坐标反变换模块和输出模块；

输入模块，被配置为用于输入初始纵横波速度模型、多分量地震记录；

网格生成模块，被配置为用于极坐标系下的网格剖分；

坐标变换模块，被配置为用于将纵横波速度模型变换到极坐标系下；

梯度求取模块，被配置为用于求取极坐标系下的纵横波梯度公式；

步长求取模块，被配置为用于求取加权速度更新步长；

速度更新模块，用于更新极坐标系下的纵横波速度；

坐标反变换模块，用于判断是否满足误差条件，将反演的纵横波速度变换到笛卡尔坐标系下；

输出模块，输出反演的纵横波速度。

一种极坐标下的弹性波纵横波速度反演方法，采用上述的一种极坐标下的弹性波纵横波速度反演系统，包括以下步骤：

步骤1：采用输入模块，输入初始纵横波速度模型、多分量地震记录，地表高程，震源和建波点位置；

步骤2：采用网格生成模块，进行极坐标系下的网格剖分；

步骤3：采用坐标变换模块，将纵横波速度模型变换到极坐标系下；将笛卡尔坐标系(x,z)下的初始纵波速度、横波速度模型变换到极坐标系下(θ,r)；

步骤4：采用梯度求取模块，求取极坐标系下的纵横波梯度公式；极坐标下关于拉梅常数λ和μ的梯度公式由下式求得：

其中，g(λ)为关于λ的梯度公式，g(μ)为关于μ的梯度公式，t为计算时间，T_max为最大计算时间，θ为角度变量，r为半径变量，x_s为震源位置，v_r为r方向的速度分量，v_θ为θ方向的速度分量，τ_θθ和τ_rr分别表示θ方向和r方向的正应力，τ_rθ表示切应力；上标‘*’表示伴随变量；

为r方向速度分量的伴随变量，为θ方向速度分量的伴随变量，

和

分别表示θ方向和r方向正应力的伴随变量，

表示切应力的伴随变量。

再通过如下所示的链式法则求取纵横波速度的梯度：

g(V_p)＝2ρV_pg(λ), (3)

g(V_s)＝-4ρV_sg(λ)+2ρV_sg(μ)， (4)

其中，V_p为纵波速度，V_s为横波速度，ρ为密度；

步骤5：采用步长求取模块，求取加权速度更新步长；通过如下方程计算模型更新量

其中，m表示纵横波速度模型向量，p_i表示第i个参数的更新量，E表示目标函数，n为反演参数的个数，i,j为计数变量，σ_j为不同参数之间的权重，H为Hessian矩阵，‘<>’表示矢量叉乘公式；

步骤6：采用速度更新模块，更新极坐标系下的纵横波速度，判断是否满足误差条件；采用如下所示的公式更新纵横波速度：

其中，k为迭代次数；

采用如下所示的公式判断是否满足误差条件，如果满足了条件，则进入步骤7，如果不满足条件，则进入步骤4：

步骤7：采用坐标反变换模块，将反演的纵横波速度变换到笛卡尔坐标系下；

步骤8：采用输出模块，输出反演的纵横波速度。

优选的，在步骤4中，将数据残差的L2范数作为目标函数：

其中，u(t,x_r,x_s)表示地震波场；R表示将地震波场限制到检波器位置的算子；d_obs(t,x_r,x_s)表示输入的多分量地震数据，x_s和x_r分别为震源位置和检波器位置的坐标；使用变分法则推导梯度公式；目标函数的变分为：

其中，δ表示变化量；模型参数的扰动δm(λ,μ)会产生波场的扰动量δu；极坐标系下二维一阶速度-应力方程为

将m+δm和u+δu带入方程(11)可得：

方程(12)减去方程(11)可得：

根据方程(13)可得δu＝LF'，其中L表示正演模拟算子，F’是新构造的震源项，由下式求得：

因此，

δE(m)＝<R(LF',(Ru-d_obs)>＝<F',L^*R^*(Ru-d_obs)>， (15)

其中，R^*表示将检波器记录的波场扩展到整个模型空间，L^*表示波场的逆时传播；因此，L^*R^*(Ru-d_obs)表示模型空间残差波场反传，这里用U^*代替L^*R^*(Ru-d_obs)，用伴随状态理论确定目标函数的梯度公式；伴随波场U^*可由下式求得

L^*U^*(m)＝d-d_obs, (16)

其中，L^*表示伴随算子；

则

通过对λ和μ求导可得关于λ和μ的梯度公式：

利用链式法则，可以求得最终V_p和V_s的梯度公式：

g(V_p)＝2ρV_pg(λ), (20)

g(V_s)＝-4ρV_sg(λ)+2ρV_sg(μ). (21)。

优选的，在步骤5中，对于多参数反演，δm可以通过线性加权求和得到多参数更新值：

目标函数的二阶展开为：

令可得

定义<Hp_i,p_i>＝h_i,j，

方程(22)可以改写为

hσ＝-b, (25)

其中σ＝(σ₁,σ₂,···)^T，只需要得到h^-1就可以求得不同参数的步长。

本发明的有益技术效果是：

本发明开创性地实现了极坐标系下的弹性波多分量纵横波速度反演方法，该方法能够很好地解决一些特殊地质环境中的速度反演难题，如在井孔环境、隧道空间等。同时将多分量方法引入到本发明中，能够同时对纵横波速度进行反演，得到更真实、准确的纵横波速度，为井孔环境、隧道空间等特殊地质构造的成像及解释提供基础，有助于提高井孔环境、隧道空间等特殊地质构造的成像及反演精度，服务于复杂地质构造的油气勘探。

附图说明

图1为本发明的一种极坐标下的弹性波纵横波速度反演方法的流程图；

图2示出井孔环境下真实的Marmousi模型；其中，(a)，(b)为Vp模型；(c)，(d)为Vs模型；(a)，(c)为笛卡尔坐标系；(b)，(d)为极坐标系；

图3示出初始Marmousi模型；其中，(a)为Vp模型；(b)为Vs模型；

图4示出反演Marmousi模型；其中，(a)为Vp模型；(b)为Vs模型；

图5示出真实模型和反演模拟的差；其中(a)为Vp模型；(b)为Vs模型；

图6主要示出从模型得到的地震记录，其中(a)，(d)为从真实Marmousi模型得到的地震记录；(b)，(e)为从反演Marmousi模型得到的地震记录；(c)，(f)为两者的差；(a)-(c)为θ方向；(d)-(f)为r方向；

图7为半无限空间下真实Marmousi模型；其中，(a)，(b)为Vp模型；(c)，(d)为Vs模型；(a)，(c)为笛卡尔坐标系；(b)，(d)为极坐标系；

图8主要示出半无限空间下的模型；其中，(a)，(b)为半无限空间下初始Marmousi模型；(c)，(d)为半无限空间下反演的Marmousi模型；(a)，(b)为Vp模型；(c)，(d)为Vs模型；

图9为本发明中一种极坐标下的弹性波纵横波速度反演系统的结构示意图。

具体实施方式

在一些特殊的环境(如井孔环境、隧道空间等)中，传统的速度反演方法无法对以井孔为辐射中心的速度进行准确的反演。为此，本发明提出了一种极坐标下的弹性波纵横波速度反演方法，开发极坐标下的弹性波纵横波速度反演，提高井孔环境、隧道空间等特殊地质构造的成像及反演精度。

下面结合附图以及具体实施例对本发明作进一步详细说明：

实施例1

一种极坐标下的弹性波纵横波速度反演系统，其结构如图9所示，包括输入模块、网格生成模块、坐标变换模块、梯度求取模块、步长求取模块、速度更新模块、坐标反变换模块、输出模块；以上各个模块顺序连接；

输入模块，被配置为用于输入初始纵横波速度模型、多分量地震记录；

网格生成模块，被配置为用于极坐标系下的网格剖分；

坐标变换模块，被配置为用于将纵横波速度模型变换到极坐标系下；

梯度求取模块，被配置为用于求取极坐标系下的纵横波梯度公式；

步长求取模块，被配置为用于求取加权速度更新步长；

速度更新模块，用于更新极坐标系下的纵横波速度；

坐标反变换模块，用于判断是否满足误差条件；将反演的纵横波速度变换到笛卡尔坐标系下；

输出模块，用于输出反演的纵横波速度。

实施例2

在上述实施例的基础上，本发明还提出一种基于黏声拟微分方程的黏声起伏地表正演模拟方法，即极坐标下的弹性波纵横波速度反演方法，其流程如图1所示，具体包括以下步骤：

步骤1：输入初始纵横波速度模型、多分量地震记录，地表高程，震源和建波点位置；

步骤2：极坐标系下的网格剖分；

步骤3：将纵横波速度模型变换到极坐标系下；将笛卡尔坐标系(x,z)下的初始纵波速度、横波速度模型变换到极坐标系下(θ,r)；

步骤4：求取极坐标系下的纵横波梯度公式；将数据残差的L2范数作为目标函数：

其中，u(t,x_r,x_s)表示地震波场；R表示将地震波场限制到检波器位置的算子；d_obs(t,x_r,x_s)表示输入的多分量地震数据，x_s和x_r分别为震源位置和检波器位置的坐标，t为计算时间。使用变分法则推导梯度公式。目标函数的变分为：

其中，E表示目标函数，m表示纵横波速度模型向量，δ表示变化量。模型参数的扰动δm(λ,μ)会产生波场的扰动量δu。极坐标系下二维一阶速度-应力方程为

其中，θ为角度变量，r为半径变量，v_θ为θ方向的速度分量，τ_θθ和τ_rr分别表示θ方向和r方向的正应力，τ_rθ表示切应力，x表示坐标位置，λ和μ为拉梅常数，ρ为密度，将m+δm和u+δu带入方程(28)可得：

方程(29)减去方程(28)可得：

根据方程(30)可得δu＝LF'，其中L表示正演模拟算子，F’是新构造的震源项，由下式求得：

因此，

δE(m)＝<R(LF',(Ru-d_obs)>＝<F',L^*R^*(Ru-d_obs)>， (32)

其中，R^*表示将检波器记录的波场扩展到整个模型空间，L^*表示波场的逆时传播。因此，L^*R^*(Ru-d_obs)表示模型空间残差波场反传，用U^*代替L^*R^*(Ru-d_obs)。用伴随状态理论确定目标函数的梯度公式。伴随波场U^*可由下式求得

L^*U^*(m)＝d-d_obs, (33)

其中L^*表示伴随算子。

则

其中，上标‘*’表示变量的伴随。通过对λ和μ求导可得关于λ和μ的梯度公式：

其中，g(λ)为关于λ的梯度公式，g(μ)为关于μ的梯度公式，T_max为最大计算时间。利用链式法则，可以求得最终V_p和V_s的梯度公式：

g(V_p)＝2ρV_pg(λ), (37)

g(V_s)＝-4ρV_sg(λ)+2ρV_sg(μ). (38)

其中，V_p为纵波速度，V_s为横波速度。

步骤5：求取加权速度更新步长；

对于多参数反演，δm可以通过线性加权求和得到多参数更新值：

其中，p_i表示第i个参数的更新量，n为反演参数的个数，i为计数变量，σ为不同参数之间的权重。

目标函数的二阶展开为：

其中，H为Hessian矩阵，‘<>’表示矢量叉乘公式。令

可得

定义<Hp_i,p_i>＝h_i,j，

方程(39)可以改写为

hσ＝-b, (42)

其中σ＝(σ₁,σ₂,···)^T。只需要得到h-¹就可以求得不同参数的步长。

步骤6：更新极坐标系下的纵横波速度，判断是否满足误差条件；采用如下所示的公式更新纵横波速度：

其中，k为迭代次数。

采用如下所示的公式判断是否满足误差条件，如果满足了条件，则进入步骤7，如果不满足条件，则进入步骤4：

步骤7：将反演的纵横波速度变换到笛卡尔坐标系下；

步骤8：输出反演的纵横波速度。

本发明开创性地实现了极坐标系下的弹性波多分量纵横波速度反演方法，该方法能够很好地解决一些特殊地质环境中的速度反演难题，如在井孔环境、隧道空间等。同时，将多分量方法引入到本发明中，能够同时对纵横波速度进行反演，得到更真实、准确的纵横波速度，为井孔环境、隧道空间等特殊地质构造的成像及解释提供基础，有助于提高井孔环境、隧道空间等特殊地质构造的成像及反演精度，服务于复杂地质构造的油气勘探。

应用实验

本发明提出一种基于黏声拟微分方程的黏声起伏地表正演模拟方法，即极坐标下的弹性波纵横波速度反演方法，应用于井孔环境Marmousi模型和半无限空间Marmousi模型数据，取得了理想的计算效果。该方法首先对井孔环境下的Marmousi模型(图2)进行试算。输入初始纵横波速度模型(图3)、多分量地震记录；极坐标系下的网格剖分；将纵横波速度模型变换到极坐标系下；求取极坐标系下的纵横波梯度公式；求取加权速度更新步长；更新极坐标系下的纵横波速度，判断是否满足误差条件；将反演的纵横波速度变换到笛卡尔坐标系下；输出反演的纵横波速度(图4)。为了论证本发明方法的正确性，将反演的结果与真实结果求差得到的结果如图5所示，可以看出，反演的结果非常接近于真实的模型。

接下来对半无限空间的Marmousi模型(图7)进行试算，以论证本发明方法应用到半无限空间环境中也可以取得较好的效果。输入初始速度模型，采用本发明方法得到的反演结果如图8所示，可以看出，本发明的反演方法能够对复杂的地质构造具有较好的反演结果。两个不同空间的模型试算证明了本发明可以在井孔环境及半无限环境中都可以得到很好的反演结果。

上述方式中未述及的部分采取或借鉴已有技术即可实现。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (3)

1.一种极坐标下的弹性波纵横波速度反演方法，包括极坐标下的弹性波纵横波速度反演系统，系统包括顺序连接的输入模块、网格生成模块、坐标变换模块、梯度求取模块、步长求取模块、速度更新模块、坐标反变换模块和输出模块；

输入模块，被配置为用于输入初始纵横波速度模型、多分量地震记录；

网格生成模块，被配置为用于极坐标系下的网格剖分；

坐标变换模块，被配置为用于将纵横波速度模型变换到极坐标系下；

梯度求取模块，被配置为用于求取极坐标系下的纵横波梯度公式；

步长求取模块，被配置为用于求取加权速度更新步长；

速度更新模块，用于更新极坐标系下的纵横波速度；

坐标反变换模块，用于判断是否满足误差条件，将反演的纵横波速度变换到笛卡尔坐标系下；

输出模块，输出反演的纵横波速度；

其特征在于，方法包括以下步骤：

步骤1：采用输入模块，输入初始纵横波速度模型、多分量地震记录，地表高程，震源和建波点位置；

步骤2：采用网格生成模块，进行极坐标系下的网格剖分；

步骤3：采用坐标变换模块，将纵横波速度模型变换到极坐标系下；将笛卡尔坐标系(x,z)下的初始纵波速度、横波速度模型变换到极坐标系下(θ,r)；

步骤4：采用梯度求取模块，求取极坐标系下的纵横波梯度公式；极坐标下关于拉梅常数λ和μ的梯度公式由下式求得：

其中，g(λ)为关于λ的梯度公式，g(μ)为关于μ的梯度公式，t为计算时间，T_max为最大计算时间，θ为角度变量，r为半径变量，x_s为震源位置，v_r为r方向的速度分量，v_θ为θ方向的速度分量，τ_θθ和τ_rr分别表示θ方向和r方向的正应力，τ_rθ表示切应力；

为r方向速度分量的伴随变量，

为θ方向速度分量的伴随变量，

和

分别表示θ方向和r方向正应力的伴随变量，

表示切应力的伴随变量；

再通过如下所示的链式法则求取纵横波速度的梯度：

g(V_p)＝2ρV_pg(λ), (3)

g(V_s)＝-4ρV_sg(λ)+2ρV_sg(μ)， (4)

其中，V_p为纵波速度，V_s为横波速度，ρ为密度；

步骤5：采用步长求取模块，求取加权速度更新步长；通过如下方程计算模型更新量

其中，m表示纵横波速度模型向量，p_i表示第i个参数的更新量，E表示目标函数，n为反演参数的个数，i,j为计数变量，σ_j为不同参数之间的权重，H为Hessian矩阵，‘<>’表示矢量叉乘公式；

步骤6：采用速度更新模块，更新极坐标系下的纵横波速度，判断是否满足误差条件；采用如下所示的公式更新纵横波速度：

其中，k为迭代次数；

采用如下所示的公式判断是否满足误差条件，如果满足了条件，则进入步骤7，如果不满足条件，则进入步骤4：

步骤7：采用坐标反变换模块，将反演的纵横波速度变换到笛卡尔坐标系下；

步骤8：采用输出模块，输出反演的纵横波速度。

2.根据权利要求1所述的一种极坐标下的弹性波纵横波速度反演方法，其特征在于，在步骤4中，将数据残差的L2范数作为目标函数：

其中，u(t,x_r,x_s)表示地震波场；R表示将地震波场限制到检波器位置的算子；d_obs(t,x_r,x_s)表示输入的多分量地震数据，x_s和x_r分别为震源位置和检波器位置的坐标；使用变分法则推导梯度公式；目标函数的变分为：

其中，δ表示变化量；模型参数的扰动δm(λ,μ)会产生波场的扰动量δu；极坐标系下二维一阶速度-应力方程为

将m+δm和u+δu带入方程(11)可得：

方程(12)减去方程(11)可得：

根据方程(13)可得δu＝LF'，其中L表示正演模拟算子，F’是新构造的震源项，由下式求得：

因此，

δE(m)＝<R(LF',(Ru-d_obs)>＝<F',L^*R^*(Ru-d_obs)>， (15)

其中，R^*表示将检波器记录的波场扩展到整个模型空间，L^*表示波场的逆时传播；因此，L^*R^*(Ru-d_obs)表示模型空间残差波场反传，这里用U^*代替L^*R^*(Ru-d_obs)，用伴随状态理论确定目标函数的梯度公式；伴随波场U^*由下式求得

L^*U^*(m)＝d-d_obs, (16)

其中，L^*表示伴随算子；

则

通过对λ和μ求导可得关于λ和μ的梯度公式：

利用链式法则，可以求得最终V_p和V_s的梯度公式：

g(V_p)＝2ρV_pg(λ), (20)

g(V_s)＝-4ρV_sg(λ)+2ρV_sg(μ). (21)。

3.根据权利要求2所述的一种极坐标下的弹性波纵横波速度反演方法，其特征在于，在步骤5中，对于多参数反演，δm通过线性加权求和得到多参数更新值：

目标函数的二阶展开为：

令

可得

定义<Hp_i,p_i>＝h_i,j，

方程(22)可以改写为

hσ＝-b, (25)

其中σ＝(σ₁,σ₂,…)^T，只需要得到h^-1就可以求得不同参数的步长。

Images