CN101630018B - 一种控制全声波方程反演过程的地震勘探数据处理方法 - Google Patents

Abstract

本发明涉及一种控制全声波方程反演过程的地震勘探数据处理方法，采集理论模型地震数据，获得深度域体积模量和密度初始模型；用伪谱法模拟声学介质中地震波场；计算理论模型与初始模型模拟数据的差；给定δ＝1.0e-6，计算误差能量E，当E≤δ，停止反演，输出反演结果，如果E＞δ，继续下列步骤；计算残差数据；计算体积模量和密度模型的共轭修改量；计算修改步长，获得梯度；对初始模型进行修改；将修改后的模型数据作为新的初始模型，通过对反演过程进行自动控制，使声波方程反演稳定、快速收敛，反演收敛速度比常规方法提高了近5倍。

Description

一种控制全声波方程反演过程的地震勘探数据处理方法

技术领域

本发明涉及一种控制地震数据全声波方程反演过程，提高地震勘探全声波方程反演精度和对地质目标(尤其是岩性地层油气藏)的勘探能力的地震勘探数据处理方法。

背景技术

地震反演是利用野外观测的地震数据，通过数学方法和计算技术求取地层物理参数，如速度、密度和阻抗等。全声波方程反演直接采用了声学介质中地震波的传播方程，与波阻抗反演、弹性阻抗反演和AVO反演相比，它能处理任意复杂的非均匀介质情况，提高地震反演的精度和对地质目标(尤其是岩性地层油气藏)的勘探能力。

全声波方程反演计算费用昂贵，限制了它在实际生产中的应用。解决这一问题的关键：一是提高声波方程计算波场的效率，另一是提高反演的收敛速度。Tarantola等(1988)提出的深度域延拓计算地震波场和Pratt等(1999)提出的通过计算部分频率的波场值实现全波场的模拟都显著提高了求解波动方程的速度。这两种方法减少了单次迭代反演所需的时间，但由于地震反演基本上都是迭代进行的，迭代次数越多，反演的时间越长，费用就越高。

对于波动方程这类巨大计算量的反演，目前采用的最好的算法是最速下降法或共轭梯度法。这两种算法收敛速度比较快，且比较稳定，但仍需要一定的迭代次数。要将波动方程反演用于实际地震数据处理，常规的最速下降法或共轭梯度法所需的迭代次数是不能被允许的，也就是说迭代次数是目前波动方程反演工业化应用的一个瓶颈。

现有最速下降法或共轭梯度法中，人们采用线性搜索法或抛物线拟合法可以在一定程度上降低迭代次数，但在每一迭代中至少增加了两次额外的计算地震波场工作，显著增加了每一步的计算成本。

发明内容：

本发明的目的是利用一种控制方法对野外实测地震数据进行声波方程反演过程的自动控制，降低波动方程反演的计算成本，提高地震反演的精度和对地层油气藏的判断能力。

本发明是一种基于声波理论对野外实测地震数据进行反演，并通过对反演过程的控制，提高反演的速度和精度。具体步骤包括：

(1).野外激发、接收获得地震波数据，用常规处理方法对该数据进行静校正、地表一致性振幅补偿和叠前去除噪音，获得炮集数据；

(2).从步骤(1)获得的炮集数据中抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；

①步骤(2)中速度分析为常规地震速度分析；

②步骤(2)中时间域层速度用速度分析获得的均方根速度和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)$

其中t_i为第i层双程旅行时，v_i为第i层的速度；

③步骤(2)中时间-深度转换是按下式将时间域的层速度转换成深度域层速度，

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

其中t_i为第i层双程旅行时，h_i为第i层的厚度，v_i为第i层的速度；

④步骤(2)中密度(模型参数1)初始模型用Gardner公式计算，

ρ＝0.31v^0.25    (3)

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

⑤步骤(2)中体积模量(模型参数2)初始模型用下列方法计算

K＝ρv²    (4)

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

(3).利用步骤2建立的初始体积模量和密度模型，进行地震波场正演模拟；

①步骤(3)中地震波场正演模拟利用的是声学介质中的波动方程；

②步骤(3)中用常规的伪谱法进行波场模拟，即对空间变量进行付氏变换，在付氏域中计算声波方程中波场对空间变量的偏导数；声波方程中波场对时间变量的偏导数用二阶中心差分法计算；

③步骤(3)中震源函数采用零相位雷克子波；

(4).按下式计算步骤(3)模拟的炮集数据与步骤1获得测量数据的差(残差数据)：

Δd＝d_obs-d_cal    (5)

其中Δd为残差数据，d_obs为测量地震数据，d_cal为计算地震数据。

(5).按下式计算误差能量，

$E=\frac{1}{2}{\mathrm{\Δd}}^t{C}_D^{-1}\mathrm{\Δd}---\left(6\right)$

式中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵；

当E＜δ时，输出反演结果，并停止反演，否则，继续步骤6～步骤9，

直到E＜δ，δ为任意给的一个非常小的数，一般取δ＝1.0e-4～1.0e-6；

步骤(5)中数据协方差矩阵用下列方法计算

$C_D=\frac{1}{X}\·\frac{{\mathrm{\σ}}_d^2}{t^{2p}}---\left(7\right)$

其中X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；

指数p，对于二维问题，取p≥0.5，对于三维情况，取p≥1；

(6).将步骤(4)获得的残差数据作为源数据，用正演模拟类似的方法进行逆时传播，获得残差数据的剩余地震波场；

①步骤(6)中的逆时传播与步骤(3)的波场正演过程相反，即解波动方程，按反时间方向进行波场传播；

②与骤(3)正演模拟所用的震源相对应，步骤(6)中进行波场模拟时，震源用的是步骤(4)计算的残差数据；

(7).将步骤3获得的正演地震波场和步骤(6)获得剩余地震波场分别对时间变量进行求导，进行零延迟互相关，然后对时间变量进行积分，对激发炮数进行累加，获得初始密度参数的共轭修改量；将步骤(3)获得的正演地震波场和步骤(6)获得剩余地震波场分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，获得初始体积参数的共轭修改量，即按下式计算模型参数的共轭修改量；

$\mathrm{\δ}\hat{K}\left(x\right)=\frac{1}{K^2\left(x\right)}\underset{s}{\mathrm{\Σ}}{\∫}_0^T\mathrm{dt}\stackrel{\·}{P}(x,t;x_s)\stackrel{\·}{\mathrm{\ψ}}(x,t;x_s)---\left(8\right)$

$\mathrm{\δ}\widehat{\mathrm{\ρ}}\left(x\right)=\frac{1}{{\mathrm{\ρ}}^2\left(x\right)}\underset{s}{\mathrm{\Σ}}{\∫}_0^T\mathrm{dt}\mathrm{grad}\mathrm{\ψ}(x,t;x_s)\·\mathrm{grad}P(x,t;x_s)$

其中P为正演波场，ψ为回传的剩余波场，或表示波场对时间的一阶偏导数，grad(·)为梯度运算，T为地震记录长度，K和ρ为模型参数，

和

为模型参数共轭量；

(8).利用步骤(7)获得的模型参数的共轭修改量，按下列方法的迭代修改模型参数，

m_n+1＝m_n-α_ng_n    (9)

其中，m＝(K，ρ)为模型，m_n+1和m_n分别为第n次修改后的模型和当前模型；α_n为第n次迭代修改步长；g_n为第n次迭代的梯度；

或①步骤(8)中梯度g_n用下列方法获得，

$g_n=(\mathrm{\δ}{\hat{K}}_n,\mathrm{\δ}{\hat{p}}_n)---\left(10\right)$

其中

和

是第n步迭代的模型参数的共轭修改量，用骤(7)公式计算。

②步骤(8)中修改步长α_n为反演过程的控制，用下列方法计算

其中，

为映射函数，k为反馈增益系数，J为地震数据对模型参数的导数，E为误差能量；

③步骤(8)②中映射函数用下列方法计算

当n＞1时，函数

当n＝1时，

④步骤(8)②中反馈增益系数取满足k＞(10/t_s)的任一数(t_s为系统固有时间，t_s取5～20范围任一值)；

(9).将修改后的模型数据作为新的初始模型，重复步骤(3)～步骤(5)。

发明的效果：

本发明通过对反演过程进行自动控制，使声波方程反演稳定、快速收敛。理论模型测试结果表明控制反演收敛速度比常规方法显著提高了近5倍。实例应用为我国西部某气田实际数据反演，经过两次迭代后，获得的反演结果与实测数据基本一致。

附图说明：

图1是可控声波方程反演流程图；

图2是理论模型合成的地震数据；

图3是反演结果与真实数据图；

图4是反演误差能量随迭代次数变化曲线；

图5是野外测量地震数据图(左为频谱，右为炮集数据)；

图6是经过两次迭代反演的结果与实测数据图；

实施例1：

理论模型为一个三层介质模型，图2中黑色曲线为三层介质模型参数随深度变化曲线。图3为理论模型合成地震数据，对该数据进行反演，具体实现步骤为：

1.采集理论模型地震数据，抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；

2.基于初始模型，用伪谱法模拟声学介质中地震波场；

3.按式(5)计算理论模型与初始模型模拟数据的差；

4.给定δ＝1.0e-6，按式(6)计算误差能量E，当E≤δ，停止反演，输出反演结果，如果E＞δ，继续下列步骤；

5.步骤(3)计算的残差数据作为输入数据，用伪谱法计算声波方程的波场，并按反时间方向传播；

6.按式(8)计算体积模量和密度模型的共轭修改量；

7.按式(11)计算修改步长，并用(10)式获得梯度；

8.用(9)式对初始模型进行修改；

9.将修改后的模型数据作为新的初始模型，重复步骤2～步骤4。

图2灰线显示了5次迭代后的反演结果，与真实模型数据吻合的很好。图4显示了分别用常规共轭梯度法和可控反演方法反演的误差能量随迭代次数变化曲线，控制反演方法收敛速度快，约为常规方法的5倍，且稳定。

实施例2：

实例2是位于我国西部的某气田。该气田为陆相湖泊沉积，水下分流河道发育。目标层埋藏深度大、埋藏时间长、成岩演化程度高、物性差，属于典型的低孔隙度、低渗透率气田。储层纵向上厚度大、平面上大面积分布，但有效砂岩单层厚度薄、横向变化大，非均质性强。沉积微相、物性分析和试气成果证实，有效储层主要发育在三角洲平原和三角洲前缘亚相中的高能分流河道中。

由于有效储层和围岩纵波速度和阻抗差异小，常规地震反演结果达不到勘探的要求。2006年我们用可控声波方程叠前反演方法对该区地震数据进行了反演。图5为野外测量地震数据(左为其频谱，右为炮集数据)，用于地震反演的输入，图6显示了两次迭代反演后钻井处的反演结果(灰线)和实际测井数据(黑线)，反演获得数据与实测数据基本一致。

具体实现步骤为：

1.采集地震数据，用常规方法对地震进行静校正、地表一致性振幅补偿和叠前去噪，形成炮集数据；

2.用常规方法抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；

3.～10.同“实施例1”中步骤2～9。

Claims (1)

1.一种控制全声波方程反演过程的地震勘探数据处理方法，其特征在于：

具体步骤包括：

(1).野外激发、接收获得地震波数据，用常规处理方法对该数据进行静校正、地表一致性振幅补偿和叠前去除噪音，获得炮集数据；

(2).从步骤(1)获得的炮集数据中抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；

①步骤(2)中速度分析为常规地震速度分析；

②步骤(2)中时间域层速度用速度分析获得的均方根速度和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)$

其中t_i为第i层双程旅行时，v_i为第i层的速度，V_r，i为第i层均方根速度，t_i-1为第i-1层双程旅行时，V_r，i-1为第i-1层均方根速度；

③步骤(2)中时间-深度转换是按下式将时间域的层速度转换成深度域层速度：

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

其中t_i为第i层双程旅行时，t_i-1为第i-1层双程旅行时，h_i为第i层的厚度，v_i为第i层的速度；

④步骤(2)中密度初始模型用Gardner公式计算：

ρ＝0.31v^0.25       (3)

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

⑤步骤(2)中体积模量初始模型用下列方法计算：

K＝ρv²            (4)

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

(3).利用步骤(2)建立的初始体积模量和密度模型，进行地震波场正演模拟；

①步骤(3)中地震波场正演模拟利用的是声学介质中的波动方程；

②步骤(3)中用常规的伪谱法进行波场模拟，即对空间变量进行付氏变换，在付氏域中计算声波方程中波场对空间变量的偏导数；声波方程中波场对时间变量的偏导数用二阶中心差分法计算；

③步骤(3)中震源函数采用零相位雷克子波；

(4).按下式计算步骤(3)模拟的炮集数据与步骤(1)获得测量数据的残差数据：

Δd＝d_obs-d_cal        (5)

其中Δd为残差数据，d_obs为测量地震数据，d_cal为计算地震数据；

(5).按下式计算误差能量：

$E=\frac{1}{2}\mathrm{\Δ}{d}^t{C}_D^{-1}\mathrm{\Δd}---\left(6\right)$

式中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵；

当E＜δ时，输出反演结果，并停止反演，否则，继续步骤(6)～步骤(9)，直到E＜δ，δ为任意给的一个非常小的数，一般取δ＝1.0e-4～1.0e-6；

步骤(5)中数据协方差矩阵用下列方法计算：

$C_D=\frac{1}{X}\·\frac{{\mathrm{\σ}}_d^2}{t^{2p}}---\left(7\right)$

其中X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；

指数p，对于二维问题，取p≥0.5，对于三维情况，取p≥1；

(6).将步骤(4)获得的残差数据作为源数据，用正演模拟类似的方法进行逆时传播，获得残差数据的剩余地震波场；

①步骤(6)中的逆时传播与步骤(3)的波场正演过程相反，即解波动方程，按反时间方向进行波场传播；

②与步骤(3)正演模拟所用的震源相对应，步骤(6)中进行波场模拟时，震源用的是步骤(4)计算的残差数据；

(7).将步骤(3)获得的正演地震波场和步骤(6)获得剩余地震波场分别对时间变量进行求导，进行零延迟互相关，然后对时间变量进行积分，对激发炮数进行累加，获得初始密度参数的共轭修改量；将步骤(3)获得的正演地震波场和步骤(6)获得剩余地震波场分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，获得初始体积参数的共轭修改量，即按下式计算模型参数的共轭修改量：

$\mathrm{\δ}\hat{K}\left(x\right)=\frac{1}{K^2\left(x\right)}\underset{s}{\mathrm{\Σ}}{\∫}_0^T\mathrm{dt}\stackrel{\·}{P}(x.t;x_s)\stackrel{\·}{\mathrm{\ψ}}(x,t;x_s)---\left(8\right)$

$\mathrm{\δ}\widehat{\mathrm{\ρ}}\left(x\right)=\frac{1}{{\mathrm{\ρ}}^2\left(x\right)}\underset{s}{\mathrm{\Σ}}{\∫}_0^T\mathrm{dt}\mathrm{grad}\mathrm{\ψ}(x,t;x_s)\·\mathrm{grad}P(x,t;x_s)$

其中P(X，t；Xs)为正演波场，ψ(X，t；Xs)为回传的剩余波场，或

表示波场对时间的一阶偏导数，grad(·)为梯度运算，T为地震记录长度，K(X)和ρ(X)为模型参数，

和为模型参数

共轭量；

(8).利用步骤(7)获得的模型参数的共轭修改量，按下列方法的迭代修改模型参数：

m_n+1＝m_n-α_ng_n    (9)

其中，m＝(K，ρ)为模型，m_n+1和m_n分别为第n次修改后的模型和当前

模型；α_n为第n次迭代修改步长；g_n为第n次迭代的梯度；

①步骤(8)中梯度g_n用下列方法获得：

$g_n=(\mathrm{\δ}{\hat{K}}_n,\mathrm{\δ}{\widehat{\mathrm{\ρ}}}_n)---\left(10\right)$

其中

和是第n步迭代的模型参数的共轭修改量，用步骤(7)公式计算；

②步骤(8)中修改步长α_n为反演过程的控制，用下列方法计算：

其中，

为映射函数，k为反馈增益系数，J为地震数据对模型参数的导数，E为误差能量；

③步骤(8)的②中映射函数用下列方法计算：

当n＞1时，函数

当n＝1时，

④步骤(8)的②中反馈增益系数取满足k＞(10/t_s)的任一数，t_s为系统固有时间，t_s取5～20范围任一值；

(9).将修改后的模型数据作为新的初始模型，重复步骤(3)～步骤(5)。

Images