CN102385066B - 一种叠前地震定量成像方法 - Google Patents

Abstract

本发明公开了一种叠前地震定量成像方法，该方法包括：采集待成像区域的地震数据，根据地震数据获得测量的多炮炮集数据，并建立包含密度和体积模量参数的当前模型；测量待成像区域内一井孔处的密度和纵波纵波速度，根据井孔处的密度和纵波纵波速度获得井孔处的体积模量，建立先验模型；对当前模型进行地震波场正演模拟生成正演地震波场数据，并进行反付氏变换生成模拟炮集数据；根据模拟炮集数据和测量的炮集数据，获得残差数据；分别计算测量的炮集数据和模型数据的协方差矩阵，计算目标函数值，判断目标函数值小于预设值时，输出成像结果，判断目标函数值大于或等于预设值时，进行迭代修改，直至目标函数值小于所述预设值时，输出成像结果。

Description

一种叠前地震定量成像方法

技术领域

本发明关于地球物理地震成像，具体的讲是一种叠前地震定量成像方法。

背景技术

弄清地下地质体的结构是地质、地球物理工作者的首要任务，地震成像是实现这一目标的主要手段。

地震成像方法有两类：射线偏移和波动方程偏移。射线偏移出现在上世纪70年代初，主要用于相对简单的地质体构造形态成像。波动方程偏移最初由J.F.Claerbout提出和实现，随后发展了一系列波动方程偏移方法，如有限差分偏移、频率-波数域偏移、克希霍夫积分偏移等。克希霍夫积分偏移是目前普遍使用的方法。波动方程偏移可用于构造非常复杂的地质体界面成像，同时能通过地震波振幅信息显示界面两侧地质体差异相对大小。

但是，现有的地震成像方法只限于提供地质体构造形态，即界面形态成像，不能提供地质体内部的物理参数数据，即地质体的定量成像，不能直接提供更多反映地质体内部特征的数据，如地质体的岩性、孔隙、含油气性等。

发明内容

本发明实施例提供了一种叠前地震定量成像方法，提供一种高效、高分辨率的叠前地震地质体构造形态及其内部密度和体积模量参数的成像方法，为油气藏勘探提供定量的地球物理数据。

本发明的目的之一，提供一种叠前地震定量成像方法，该方法具体步骤包括：（1）采集待成像区域的地震数据，根据地震数据获得测量的多炮炮集数据；（2）根据测量的炮集数据建立包含密度和体积模量参数的当前模型；（3）测量待成像区域内一井孔处的密度和纵波速度，根据井孔处的密度和纵波速度获得井孔处的体积模量，根据井孔处的密度和体积模量建立先验模型；（4）对当前模型进行地震波场正演模拟生成正演地震波场数据，对正演地震波场数据进行反付氏变换生成模拟炮集数据；（5）根据模拟炮集数据和测量的炮集数据，获得残差数据；（6）计算测量的炮集数据的协方差矩阵和模型数据的协方差矩阵，根据残差数据、当前模型、先验模型、多炮炮集数据协方差矩阵、模型协方差矩阵计算目标函数值；（7）比较目标函数值与预设值的大小：判断目标函数值大于或等于所述预设值时：根据获得的残差数据，获得频率域剩余地震波场数据；根据正演地震波场数据和剩余地震波场数据获得密度和体积模量的共轭修改量；根据密度和体积模量的共轭修改量对当前模型进行迭代修改，然后进行步骤（4）～步骤（7）；判断目标函数值小于所述预设值时，根据当前模型的密度和体积模量输出定量成像结果。

本发明利用频率域单程声波方程延拓和全频带炮集数据，提供地下地质体定量成像，成像精度高，计算效率高。为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1本发明叠前地震数据成像流程图；

图2为当前模型的生成步骤；

图3为待成像区域先验模型的生成步骤；

图4采集数据处理后用于成像的炮集数据图；

图5炮集数据频谱图；

图6初始密度模型图；

图7初始体积模量模型图；

图8地层密度定量成像剖面图；

图9地层体积模量定量成像剖面图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供了一种叠前地震定量成像方法，该方法是一种基于频率域单程声波方程延拓的地震定量成像方法，从采集的地震炮集数据中获取地下构造形态和地质体密度及体积模量参数。

如图1所示，为本发明叠前地震定量成像方法的流程图

步骤S101，采集待成像区域的地震数据，根据所述的地震数据获得测量的多炮炮集数据；

步骤S102，根据所述的测量的炮集数据建立包含密度和体积模量参数的当前模型；

步骤S103，测量待成像区域内一井孔处的密度和纵波速度，根据所述的井孔处的密度和纵波速度获得所述井孔处的体积模量，根据所述井孔处的密度和体积模量建立先验模型；

步骤S104，对所述的当前模型进行地震波场正演模拟生成正演地震波场数据，对所述的正演地震波场数据进行反付氏变换生成模拟炮集数据；

步骤S105，根据所述的模拟炮集数据和测量的炮集数据，获得残差数据；

步骤S106，计算测量的炮集数据的协方差矩阵和模型数据的协方差矩阵，根据所述的残差数据、当前模型、先验模型、多炮炮集数据协方差矩阵、模型协方差矩阵计算目标函数值；

步骤S107，判断目标函数值是否小于预设值：

判断是，执行步骤S111，根据当前模型的密度和体积模量输出定量成像结果；

判断否，进行下述步骤：

步骤S108，根据残差数据获得频率域剩余地震波场数据；

步骤S109，根据正演地震波场数据和剩余地震波场数据获得的密度参数和体积模量参数的共轭修改量；

步骤S110，根据密度参数和体积模量参数的共轭修改量对当前模型进行迭代修改，然后重新执行步骤S104，步骤S105，步骤S106和步骤S107，直至目标函数值小于预设值时，输出成像结果。

为更好理解本发明，下面结合上述步骤对本发明实施例做进一步详细描述。本发明是一种基于频率域单程声波方程延拓的叠前地震定量成像方法，从采集的地震炮集数据中获取地下构造形态和地质体密度及体积模量参数。

本发明公开的叠前地震定量成像方法在本实施例中的具体步骤如下：

1.采集地震数据，对地震数据进行数据处理和频谱分析，获得炮集数据和频谱数据。在该步骤中的数据处理包括解编、坏炮、坏道剔除、静校正、地表一致性振幅补偿和叠前去除噪音。频谱分析采用常规地震数据频谱分析方法。

2.如图2所示为当前模型的生成步骤，从获得的炮集数据中抽取共中心点道集数据，进行纵波速度分析，层纵波速度计算和时间-深度转换，获得包含深度域体积模量和密度参数的当前模型。其中，该步骤中纵波速度分析为常规地震纵波速度分析。

上述的层纵波速度计算采用纵波速度分析获得的均方根纵波速度和如公式（1）所示Dix公式进行计算：

$v_i^2=\frac{v_{r,i}^2{t}_i-v_{r,i-1}^2{t}_{i-1}}{t_i-t_{i-1}}---\left(1\right)$

公式（1）中的t_i为第i层双程旅行时间，v_i为第i层的纵波速度，其中，v_r,i为i层的均方根纵波速度，Δt_i为地震波通过第i层的时间。

上述的时间-深度转换按下式将时间域的层纵波速度转换成深度域层纵波速度，h_i为第i层的厚度。

$h_i=\frac{1}{2}(t_i-t_{i-1})v_i---\left(2\right)$

根据上述的速度v_i计算当前模型的密度参数，用如公式（3）所示的Gardner公式计算：

ρ＝0.31v^0.25（3）

其中，ρ为密度，v为纵波速度。

上述的当前模型的体积模量用下列方法计算

K＝ρv²（4）

其中，K为体积模量，ρ为密度，v为纵波速度。

3.如图3所示为先验模型的生成步骤，分别利用声波测井和密度测井获得井孔处的纵波速度和密度，根据井孔处的纵波速度和密度获得井孔处的体积模量，用层位数据对获得井孔处的密度和体积模量进行控制，对井孔处的密度和体积模量进行内插和外推，获得密度和体积模量先验模型。层位数据为常规地震解释获得，井孔处体积模量用上述的公式（4）计算，内插和外推为常规线性方法。

4.利用步骤2建立的当前模型，进行地震波场正演模拟，获得正演地震波场数据和正演地面接收的模拟炮集数据。其中，地震波场正演模拟利用的是声学介质中的单程波动方程。

用常规的频率域单程声波方程延拓进行波场模拟，即对时间变量进行付氏变换，在频率域中计算声波方程中波场对时间变量的偏导数，沿深度方向延拓，计算正演地震波场数据。

步骤4中震源函数采用零相位雷克子波，震源主频为步骤1中采用常规地震数据频谱分析获得的主频数据。

5.对步骤4中获得模拟炮集数据进行反付氏变换，并按下式（5）计算与步骤1获得与测量数据的差，即残差数据：

Δd＝d_obs-d_cal（5）

其中，Δd为残差数据，d_obs为测量获得的炮集数据，d_cal为模拟获得的炮集数据。

6.计算目标函数，公式如下所示：

$S=\frac{1}{2}[\mathrm{\Δ}{d}^t{C}_D^{-1}\mathrm{\Δd}+{(m-m_{\mathrm{pior}})}^t{C}_M^{-1}(m-m_{\mathrm{pior}})]---\left(6\right)$

（6）式中S为目标函数；Δd^t为Δd的转置；C_D为测量炮集数据协方差矩阵；m=(ρ，K)为当前模型，第一次迭代修改时为步骤2获得的初始模型；m_pior=(ρ_pior，K_poir)为步骤3获得的先验模型，ρ_pior为先验密度，K_poir为先验体积模量；C_M为当前模型协方差矩阵。

当S<ε时，根据当前模型的密度和体积模量输出定量成像结果，并停止迭代修改；否则，继续步骤7～步骤10，直到S<ε。预设值ε为任意给的一个非常小的数，一般取ε=1.0e^-4～1.0e^-6。

在公式（6）中的测量的炮集数据协方差矩阵C_D用下列方法计算

$C_D=\frac{1}{X_D}\&CenterDot;\frac{{\mathrm{\&sigma;}}_d^2}{t^{2p}}---\left(7\right)$

其中，X_D为地震波从炮点经过地下界面反射到达测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，对于二维问题，取p≥0.5，对于三维情况，取p≥1。

（2）步骤6中当前模型协方差矩阵用下列方法计算

$C_M=\frac{1}{X_M}\&CenterDot;{\mathrm{\&sigma;}}_M^2---\left(8\right)$

其中，X_M为模型上的任一点到模型中点的距离，σ_M为当前模型方差。

7.将步骤5获得的残差数据进行付氏变换，并作为源数据，用正演模拟类似的方法，在频率域进行回传延拓，获得频率域残差数据的剩余地震波场。

在步骤7中的回传延拓与步骤4的波场正演过程相反，即在频率域解单程声波方程，按反时间方向进行波场延拓。

与步骤4正演模拟所用的震源相对应，在步骤7中进行波场模拟时，震源用的是步骤5计算获得的残差数据。

8.在频率域分别计算步骤4获得的正演地震波场的全频带数据对时间变量的一阶导数和步骤7获得剩余地震波场全频带数据对时间变量的一阶导数，然后将这两个一阶导数相乘，对激发炮数进行累加，进行反付氏变换，获得体积模量的共轭修改量，通过下式获得体积模量的共轭修改量：

$\mathrm{\&delta;}\hat{K}\left(x\right)=\frac{1}{K^2\left(x\right)}\underset{s}{\mathrm{\&Sigma;}}{\&Integral;}_0^T\mathrm{dt}\stackrel{\&CenterDot;}{P}(x,t;x_s)\stackrel{\&CenterDot;}{\mathrm{\&psi;}}(x,t;x_s)$

其中P为正演地震波场，ψ为剩余地震波场，T为记录地震时间长度，K为体积模量，为体积模量的共轭修改量；

在频率域计算步骤4获得的正演地震波场的全频带数据对空间变量的梯度和步骤7获得剩余地震波场的全频带数据对空间变量的梯度，将获得的两个梯度相乘，对激发炮数进行累加，反付氏变换，按下式对激发炮数进行累加，获得密度的共轭修改量：

$\mathrm{\&delta;}\widehat{\mathrm{\&rho;}}\left(x\right)=\frac{1}{{\mathrm{\&rho;}}^2\left(x\right)}\underset{s}{\mathrm{\&Sigma;}}{\&Integral;}_0^T\mathrm{dt}\mathrm{grad}\mathrm{\&psi;}(x,t;x_s)\&CenterDot;\mathrm{grad}P(x,t;x_s)$

其中P为正演地震波场，ψ为剩余地震波场，T为地震记录长度，grad(·)为梯度运算，ρ为密度，为密度共轭修改量。

9.利用步骤8获得的体积模量和密度的共轭修改量，对当前模型的密度和体积模量进行迭代修改，获得体积模量和密度参数，迭代修改采用常规共轭梯度法。

10.将修改后的模型数据作为新的当前模型，重复步骤4～步骤6。

本发明提出了利用全频带地震炮集数据和频率域单程声波方程延拓反演获取地下地质体定量成像，具有成像精度高，计算效率高的优点。

位于我国西部的某气田。该气田为陆相湖泊沉积，水下分流河道发育。目标层埋藏深度大、埋藏时间长、成岩演化程度高、物性差，属于典型的低孔隙度、低渗透率气田。储层纵向上厚度大、平面上大面积分布，但含气有效砂岩单层厚度薄、横向变化大，非均质性强，与围岩差异小，反射地震特征模糊，常规方法难辨识。

采用本发明提出的方法对研究区的2D地震数据进行处理，获得地层的密度和体积模量定量成像。图4为采集的地震数据经过解编、静校正、地表一致性振幅补偿和部分去噪等处理后的炮集数据，用于定量成像的输入数据。图5为炮集数据频谱，其主频用于正演模拟炮集数据震源的主频；图6和图7分别是初始密度和体积模量模型，用于正演模拟和迭代修改的初始模型；图8是经过两次迭代修改后所得的地层密度成像剖面；图9是经过两次迭代修改后所得的体积模量成像剖面。密度和体积模量分布也成层状，大套地层得到很好成像，大套地层内部弱差异地质体（即有效储层）成像清晰，连续性比较差，与该区的沉积微相和物性分析一致。该测线上有两口井，均为产气井，利用成像数据预测的气藏分布与钻井结果基本吻合。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (6)

1.一种叠前地震定量成像方法，其特征在于，所述的方法包括：

(1)采集待成像区域的地震数据，根据所述的地震数据获得测量的多炮炮集数据；

(2)根据所述的测量的炮集数据建立包含密度和体积模量参数的当前模型；

(3)采用声波测井和密度测井获得井孔处的纵波速度和密度，采用常规地震解释获得的层位数据对所述的井孔处的密度和体积模量进行控制，并对所述的井孔处的密度和体积模量进行内插和外推，获得待成像区域包含密度和体积模量参数的先验模型；

(4)对所述的当前模型采用频率域单程声波方程延拓进行地震波场正演模拟生成正演地震波场数据，对所述的正演地震波场数据进行反付氏变换生成模拟炮集数据；

(5)根据所述的模拟炮集数据和测量的炮集数据，获得残差数据；

(6)计算测量的炮集数据的协方差矩阵和当前模型协方差矩阵，根据所述的残差数据、当前模型、先验模型、测量的炮集数据协方差矩阵、当前模型协方差矩阵计算目标函数值，其中，采用下式计算目标函数：

$S=\frac{1}{2}\&lsqb;{\mathrm{\&Delta;d}}^t{C}_D^{-1}\mathrm{\&Delta;}d+{(m-m_{pior})}^t{C}_M^{-1}(m-m_{pior})\&rsqb;$

其中，S为目标函数；Δd为残差数据；Δd^t为Δd的转置，m＝(ρ，K)为当前模型，ρ为所述的当前模型的密度，K为所述的当前模型的体积模量；m_pior＝(ρ_pior，K_poir)为所述的先验模型，ρ_pior为所述的先验模型的密度，K_poir为所述的先验模型的体积模量；

当前模型协方差矩阵由下式获得：

$C_M=\frac{1}{X_M}\&CenterDot;{\mathrm{\&sigma;}}_M^2$

其中X_M为模型上的任一点到模型中点的距离，σ_M为当前模型方差；

(7)比较所述的目标函数值与预设值的大小：

判断所述的目标函数值大于或等于所述预设值时：根据获得的残差数据，获得频率域剩余地震波场数据；根据所述的正演地震波场数据和剩余地震波场数据获得密度和体积模量的共轭修改量；根据所述的密度和体积模量的共轭修改量对所述的当前模型进行迭代修改，然后进行步骤(4)～步骤(7)；判断所述的目标函数值小于所述预设值时，根据当前模型的密度和体积模量输出定量成像结果；其中，

根据获得的残差数据，获得频率域残差数据的剩余地震波场数据包括：

将获得的残差数据进行付氏变换，将付氏变换后的残差数据作为源数据，采用正演模拟的方法，在频率域进行回传延拓，获得所述的频率域残差数据的剩余地震波场数据。

2.如权利要求1所述的方法，其特征在于，所述的步骤(2)包括：

从获得的测量炮集数据中抽取共中心点道集数据；

对所述的共中心点道集数据进行速度分析、层速度计算和时间-深度转换，获得包含深度域体积模量和密度的当前模型。

3.如权利要求1所述的方法，其特征在于，所述的步骤(5)根据所述的模拟炮集数据和测量的炮集数据获得残差数据包括：采用下式计算残差数据：

Δd＝d_obs-d_cal

其中，Δd为残差数据，d_obs为测量的炮集数据，d_cal为模拟炮集数据。

4.如权利要求1所述的方法，其特征在于，所述的测量的炮集数据协方差矩阵由下式获得：

$C_D=\frac{1}{X_D}\&CenterDot;\frac{{\mathrm{\&sigma;}}_d^2}{t^{2p}}$

其中，X_D为地震波从炮点经过地下界面反射到达测量点的距离；t为地震波从炮点经过地下界面反射到达测量点的时间；σ_d为测量的炮集数据方差；指数p，二维情况下，取p≥0.5，对于三维情况，取p≥1。

5.如权利要求1所述的方法，其特征在于，所述的步骤(7)中的根据所述的正演地震波场数据和剩余地震波场数据获得密度参数和体积模量参数的共轭修改量包括：

在频率域分别计算所述的正演地震波场数据和剩余地震波场数据的全频带数据对时间变量的一阶导数，将所述的两个一阶导数相乘，进行反付氏变换，按下式对激发炮数进行累加，获得体积模量的共轭修改量：

$\mathrm{\&delta;}\hat{K}\left(x\right)=\frac{1}{K^2\left(x\right)}\underset{s}{\&Sigma;}{\&Integral;}_0^T\mathrm{dt}\stackrel{\&CenterDot;}{P}(x,t;x_s)\stackrel{\&CenterDot;}{\mathrm{\&psi;}}(x,t;x_s)$

其中P为正演地震波场，ψ为剩余地震波场，T为地震记录时间长度，K为体积模量，为体积模量的共轭修改量；

在频率域分别计算所述的正演地震波场和剩余地震波场的全频带数据对空间变量的梯度，将所述的两个梯度相乘，进行反付氏变换，按下式对激发炮数进行累加，获得密度的共轭修改量：

$\mathrm{\&delta;}\widehat{\mathrm{\&rho;}}\left(x\right)=\frac{1}{{\mathrm{\&rho;}}^2\left(x\right)}\underset{s}{\&Sigma;}{\&Integral;}_0^{T}dtgrad\mathrm{\&psi;}(x,t;x_s)\&CenterDot;gradP(x,t;x_s)$

其中P为正演地震波场数据，ψ为剩余地震波场数据，T为地震记录时间长度，ρ为密度，为密度共轭修改量。

6.如权利要求1所述的方法，其特征在于，所述的步骤(7)中的根据所述的密度参数和体积模量参数的共轭修改量对所述的当前模型进行迭代修改包括：

采用常规共轭梯度法对所述的当前模型进行迭代修改。