CN104570082B - 一种基于格林函数表征的全波形反演梯度算子的提取方法 - Google Patents

Abstract

本发明公开了一种基于格林函数表征的全波形反演梯度算子的提取方法，通过建立全波形反演参数模型；采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；将脉冲点源置于地表原接收点处，首先将脉冲点源置于地表原接收点处，计算该点源产生的格林函数波场；利用优化有限差分正演模拟计算地下各点的虚震源波场；利用格林函数和虚震源波场褶积，计算波场对模型参数的导数，得到波场逆时传播算子；将波场逆时传播算子作用于波场误差，得到全波形反演的梯度算子。本发明的有益效果是，对格林函数波场、虚震源波场和波场误差的模拟精度高，数值频散小，可以有效处理边界反射，效率较高。

Description

一种基于格林函数表征的全波形反演梯度算子的提取方法

技术领域

本发明涉及石油地震勘探速度建模技术领域，尤其涉及一种基于格林函数表征的全波形反演梯度算子的提取方法

背景技术

随着油气勘探的不断深入，人们越来越希望由地震资料获得更多的岩石物性信息，而地震波形反演的目的正是获取一个预测地震记录与实测地震记录拟合最佳的地质模型。与其他反演方法相比，波形反演利用的是整个或部分地震记录波形，此波形不仅包含地震波记录震相的运动学特征(如旅行时、波速等)，还包含地震波动力学特征(如振动方向、振幅、周期、相位等)，可以有效的排除偶然因素的影响，从而提高了计算的可靠性和稳定性，有效地提高反演的精度和分辨率。地震波形反演利用地震波场反演介质的弹性性质，即在利用波动方程描述地震波传播过程的基础上，从地震波形包含的信息中获取介质的性质。全波形反演方法能利用叠前地震波场的运动学和动力学信息，重建地下速度结构，具有揭示复杂地质背景下构造与岩性细节信息的潜力。波形反演参数模型更新量的梯度可以用梯度算子来表示，它是波形反演的关键技术之一，利用梯度算子和海瑟矩阵得到参数模型修正量，通过迭代来修正模型，从而建立满足观测数据和模型数据的最佳匹配的地下地质模型。

上世纪80年代Tarantola等人提出了基于广义最小二乘反演理论的时间域全波形反演方法，对近20多年多维地震反演理论的发展产生了深远的影响。该方法使用理论波场与实际波场误差的L2泛函作为反演的目标函数，对模型进行网格离散化并利用双程波动方程数值模拟方法（如有限差分法、伪谱法等）模拟波场传播过程，通过梯度寻优实现模型参数迭代更新。该方法可适用于井间地震、VSP、广角地震、反射地震等多种观测系统，反演时能够充分利用多种波形信息，而不局限于只利用反射波信息。为提高计算效率，80年代末90年代初Pratt等人将全波形反演理论推广到频率域，形成了频率域全波形反演方法，也称波形层析成像方法。为了避免反演陷入局部极小，全波形反演方法对初始模型精度要求严格。针对地震数据频带宽度有限、反演初始速度模型获取困难等问题，Shin提出利用阻尼波场零频分量反演低频模型作为频率域波形反演的初始模型，即Laplace域全波形反演方法，为全波形反演理论与应用研究注入了新的活力。近几年来国内外全波形反演方法应用研究发展迅速，例如Bunks et al.(1995)研究了时间域内的多尺度波形反演。Pratt et al.(1990)发展了全波形反演成像的并进行油气开发CO2流体注入前面破裂的监测应用，Smithyman et al.(2009)使用波形层析成像进行了近地表探测，Operto et al.(2006)进行全波形反演处理多次覆盖的海洋数据得到地壳尺度的结构。这些研究为区域深部构造及演化分析、浅表层环境调查、宏观速度场建模与成像、岩性参数反演提供了新的有力手段。

发明内容

为了提高全波形反演梯度算子的计算精度，本发明提供了基于格林函数表征的全波形反演梯度算子的提取方法，并利用优化有限差分数值计算虚震源波场和波场误差，该方法可以提高全波形反演的计算精度。

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

基于格林函数表征的全波形反演梯度算子的提取方法，其特征在于包括如下步骤：

1）建立全波形反演参数模型；

2）采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；

3）将脉冲点源置于地表原接收点处，利用下式表示优化有限差分算法计算该点源产生的非均匀介质格林函数波场：

其中G为非均匀介质格林函数波场，c_m为优化有限差分系数，V为介质速度，△x和△z为空间采样间隔，△t为时间采样间隔，i,k为x和z空间坐标，n为时间坐标，δ⁰为脉冲点源，is,ks为点源空间坐标；

4）利用下式表示优化有限差分正演模拟算法计算地下各点的虚震源波场：

其中c_m为优化有限差分系数，U为根据当前模型正演模拟的波场，利用该式就可以用有限差分方法计算虚震源；

5）将步骤3）得到的格林函数波场和步骤4）得到的虚震源波场按照下式进行褶积，计算波场逆时传播算子：

其中J为波场逆时传播算子，G为非均匀介质格林函数波场，f′为虚震源波场，U为根据当前模型正演模拟的波场，V为介质速度；

6）利用步骤4)中所述优化有限差分正演模拟算法公式得到叠前单炮合成数据，并其与观测地震数据求差得到波场残差；

7）根据下式提取格林函数表征的用于全波形反演的梯度算子：

其中：RU-d_obs是波场残差。

本发明的有益效果是，对格林函数波场、虚震源波场和波场误差的模拟精度高，可以有效处理边界反射，效率较高。

附图说明

图1是全波形反演梯度算子的格林函数表征的工作流程图。

图2是均匀介质格林函数的振幅图。其中波速v=2000m/s，网格大小为101×101，脉冲点源位于模型中心处，振幅强度为1。图中表明均匀介质格林函数的振幅等值线为同心圆，并且随着与震源距离的增大，其振幅是逐渐减小的。

图3是真实参数模型图。

图4是初始参数模型图。

图5是观测记录。

图6是合成记录和观测记录求差得到波场误差。

图7a时刻一波场误差逆时传播波场。

图7b时刻二波场误差逆时传播波场。

图7c时刻三波场误差逆时传播波场。

图7d时刻四波场误差逆时传播波场。

图8复杂模型全波形反演的格林函数表征的梯度算子。

图9最终反演结果。

具体实施方式

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

实施例1。一种基于格林函数表征的全波形反演梯度算子的提取方法，其特征在于包括如下步骤：

建立全波形反演参数模型；

采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；

首先将脉冲点源置于地表原接收点处，计算该点源产生的格林函数波场；

利用优化有限差分正演模拟计算地下各点的虚震源波场；

利用格林函数和虚震源波场褶积，计算波场对模型参数的导数，得到波场逆时传播算子；

将波场逆时传播算子作用于波场误差，得到全波形反演的梯度算子。

其具体步骤为：

根据定量地震学，对于声波方程相应格林函数满足

其中G=G(x，z，t，x′，z′，t′)为格林函数，δ＝δ(x′，z′，t′)为时空δ函数，为空间Laplace算子，V为空间速度分布。格林函数可理解为位于(x′，z′)的t′时刻的集中体力δ(x′z′，t′)在速度介质产生的波场函数。在均匀无限各向同性介质中，其解为

而从弹性波方程出发，均匀、各向同性、无限介质弹性动力学格林函数为

利用格林函数在时间和空间上的积分可以求得地震位移。

波形反演的参数模型更新量表达式为

Δm=-H·g (4)

其中:Δm是参数模型的更新量；g＝J·(Ru-d_obs)是参数模型更新量的梯度算子；H=J·J^T是Hessian矩阵；Rn-d_obs是波场残差；是地震波场对模型参数的导数，利用格林函数进行表征。

二维声波波动方程可以表示为：

其中U为波场函数，f为震源函数。根据格林公式，可将声波方程（5）的解表示为格林函数积分形式

U(x，z，t)=∫dt∫∫dx′dz′G(x，z，t，x′，z′，t′)f(x_s，z_s，t′) (6)

假设速度扰动δV引起的散射场为δU则(V+δV)和(U+δU)同样满足波动方程（5）：

利用近似可得

同理，根据格林公式并只考虑单点的速度扰动，方程（7）的格林函数表征为

两边同除以δV则可以得出位移场对模型参数导数的表达式：

即波场对介质速度的偏导数可表示为格林函数场与的褶积，其中被称为虚震源。

因此，波形反演梯度算子的格林函数表征为：

考虑互换原理，将脉冲点源置于地表激发，将由此产生的格林函数波场与地下介质各点的虚震源波场褶积，可以减少正演计算次数，提高计算效率。

利用优化有限差分算法，计算非均匀介质的格林函数和虚震源波场。

1)格林函数计算

二维格林函数满足如下的方程：

其中G(x，z，t)为点源的场，即格林函数；V(x，z)为介质速度。

均匀介质的格林函数可采用解析解进行计算，而非均匀介质的格林函数可采用数值方法进行计算。本发明采用优化有限差分算法计算非均匀介质的格林函数。

利用Taylor展开进行时间和空间差分，可得时间离散形式：

再根据空间差分的形式

可以得出非均匀介质格林函数数值计算公式

其中G_m为优化有限差分系数。利用该式就可以用优化有限差分方法计算非均匀介质格林函数波场。

2）虚震源计算

虚震源可表示为其中是波场对时间的二阶偏导数，因为当前时刻波场u已知，则根据波动方程可以在空间域进行计算。

当前波场满足波动方程为：

其中U(x，z，t)为波场；V(x，z)为介质速度。因此波场对时间的二阶偏导数可表示为：

根据空间差分的形式

可以得出波场对时间的二阶偏导数的计算公式

因此，虚震源计算公式为：

其中c_m为优化有限差分系数。利用该式就可以用有限差分方法计算虚震源。

3）格林函数表征的梯度算子

波场对介质速度的偏导数的格林函数表征为：

因此，利用计算的格林函数场和虚震源波场进行褶积，则可得到波场对介质速度的偏导数，即波场逆时传播算子，再将该算子作用于波场误差上实现波场误差的逆时传播，最终得到全波形反演的格林函数表征的梯度算子：

实验例1。应用本发明全波形反演的格林函数表征的梯度算子，应用于胜利油田某块区的地震数据，取得了理想的成像效果。图3为真实参数模型；图4为初始参数模型；一方面，在当前参数模型基础上利用优化有限差分算法计算合成地震记录，再与观测记录图5求差得到波场误差图6；另一方面利用优化有限差分算法计算格林函数波场，再与计算的虚震源波场褶积，得到波场逆时传播算子；将波场逆时传播算子作用于波场误差，实现波场误差的逆时传播(图7a到图7d)，得到复杂模型全波形反演的格林函数表征的梯度算子（图8），通过多次迭代，获得最终反演结果（图9）。可以看出模拟记录和逆时传播波场数值频散小，模拟精度高，人为边界反射压制较好，计算的梯度算子精度高，反演结果准确。

Claims (1)

1.基于格林函数表征的全波形反演梯度算子的提取方法，其特征在于包括如下步骤：

1)建立全波形反演参数模型；

2)采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；

3)将脉冲点源置于地表原接收点处，利用下式表示优化有限差分算法计算该点源产生的非均匀介质格林函数波场：

$\begin{array}{c}{G}_{i,k}^{n+1}=2G_{i,k}^n-G_{i,k}^{n-1}+\frac{1}{2}{\left(\frac{V\mathrm{\Δ}t}{\mathrm{\Δ}x}\right)}^2\[c_0{G}_{i,k}+\underset{m=1}{\overset{M}{\Σ}}{c}_m(G_{i+m,k}+G_{i-m,k})\]\\ +\frac{1}{2}{\left(\frac{V\mathrm{\Δ}t}{\mathrm{\Δ}z}\right)}^2\[c_0{G}_{i,k}+\underset{m=1}{\overset{M}{\Σ}}{c}_m(G_{i,k+m}+G_{i,k-m})\]+{\mathrm{\δ}}_{is,ks}^0\end{array}$

其中G为非均匀介质格林函数波场，c_m为优化有限差分系数，V为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标，δ⁰为脉冲点源，is,ks为点源空间坐标；

4)利用下式表示优化有限差分正演模拟算法计算地下各点的虚震源波场：

$\begin{array}{c}{f}^{\′}=\frac{2}{V^3}\frac{{\∂}^{2}U}{\∂t^2}{|}_{i,k}^n=\frac{1}{{\mathrm{V\Δx}}^2}\[c_0{U}_{i,k}+\underset{m=1}{\overset{M}{\Σ}}{c}_m(U_{i+m,k}+U_{i-m,k})\]\\ +\frac{1}{{\mathrm{V\Δz}}^2}\[c_0{U}_{i,k}+\underset{m=1}{\overset{N}{\Σ}}{c}_m(U_{i,k+m}+U_{i,k-m})\]\end{array}$

其中c_m为优化有限差分系数，U为根据当前模型正演模拟的波场，利用该式就可以用有限差分方法计算虚震源；

5)将步骤3)得到的格林函数波场和步骤4)得到的虚震源波场按照下式进行褶积，计算波场逆时传播算子：

$J=\frac{\∂U}{\∂V}=G*f^{\′}=G*\frac{2}{V^3}\frac{{\∂}^{2}U}{\∂t^2}$

其中J为波场逆时传播算子，G为非均匀介质格林函数波场，f′为虚震源波场，U为根据当前模型正演模拟的波场，V为介质速度；

6)利用步骤4)中所述优化有限差分正演模拟算法公式得到叠前单炮合成数据，并其与观测地震数据求差得到波场残差；

7)根据下式提取格林函数表征的用于全波形反演的梯度算子：

$g=(G*\frac{2}{V^3}\frac{{\∂}^{2}U}{\∂t^2})\·(RU-d_{obs})$

其中：RU-d_obs是波场残差。