CN101545984A - 基于小波变换的地震相干体计算方法 - Google Patents

Abstract

基于小波变换的地震相干体计算方法，步骤如下：1)对地震资料进行预处理，2)对处理过的地震资料进行精细频谱分析，确定地震资料有效频谱范围；3)对地震资料进行层位解释，并进行内插、平滑等处理；4)以地震数据、解释的地震层位为输入，利用高分辨导数小波函数，对输入地震数据的某一个地震道做小波变换，获得小波变换分频结果；5)计算瞬时小波域的瞬时振幅或/和瞬时频率或/和瞬时相位；6)利用基于小波瞬时相位的相干体算法或基于小波变换的相干体算法计算瞬时相干体；7)重构瞬时相干体，获得高频或低频相干数据体；8)对每一个地震道重复上述4－7步骤，得到所有地震道的相干计算结果；9)利用所求取的瞬时相干数据进行制图。

Description

基于小波变换的地震相干体计算方法

技术领域

本发明涉及油气田勘探技术领域，属于地震资料解释范畴，具体地说是一种基于小波变换的地震相干体计算方法。

背景技术

地震相干是指相邻地震道之间地震属性，如：波形、振幅、频率、相位等相似程度的测量。

一般情况下，现在所作的相干都是基于振幅的计算，利用多道相似性将三维振幅数据体经计算转化为相关系数数据体，在显示上强调不相关异常，突出不连续性。它的前提在于：假设地层是连续的，地震波有变化、也是渐变的，因此相邻道、线之间是相似的。当地层连续性遭到破坏而发生变化时，如断层、尖灭、侵入、变形等，导致地震道之间的波形特征发生变化，进而导致局部道与道之间的相关性表现边缘相似性的突变，地层边界、特殊岩性体的不连续性会得到低相关值的轮廓。

通过计算地震相干数据体突出那些不相干的地震数据和提取三维相关属性体，由纵向和横向上局部的波形相似性可以得到三维地震相关性的估计值，从而把三维反射振幅数据体转换成三维相似系数或相关值的数据体。

地震相干的最有效和最重要的应用之一是同相轴的检测，特别是当同相轴的幅值较小并隐藏在噪声中时，除检测外，该测量还对同相轴的强度赋予一个定量的值，若其估值是一个能量归一化的度量值，则很容易转换为信噪比，因此相干可以对数据质量给出评估，它的计算可以在相干较弱或被噪声干扰的情况下提供出数据相似性的定量值，互相关则是它的计算基础。目前主要有如下几种计算方法。

1)C1相干算法

目前的应用软件，如Landmark、GeoFrame等算法大部分都是以经典的归一化互相关为基础的计算(Bahorich和Farmer，1995，1996)，称之为第一代算法C1。它利用主测线和联络测线方向的相关系数合成主联方向的相关系数。C1算法具有计算量小、易于实现的优点，但是，受资料的限制大。

设相邻两地震道为x(n)，y(n)时窗长度为K，那么二者在纵测线l延迟的互相关函数为：

$C_{\mathrm{xy}}\left(l\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}x\left(i\right)\·y(i-l)$

假设p和q分别代表地层的倾角和方位角，那么对三维三道和多道情况分别用如下公式计算其相干值：

三道：

J道(J>3)： $C_1(p,q)={\left|\underset{j=2}{\overset{J}{\mathrm{\Π}}}\frac{C_{1j}}{\sqrt{C_{11}\·C_{\mathrm{jj}}}}\right|}^{1/J-1}$

则C1算法的相干值为： $C_1=\underset{p,q}{\max }{C}_1(p,q)$

2)C2相干算法

该算法和C3算法引入了协方差矩阵，使其可对任意道数进行相似分析，估计其相干性。C2相干算法除了在噪声环境下更稳健地测量相干之外，垂直分析时窗能被限制在只有几个时间样点范围内，能够精确做出薄而小的地层特征图。

假设有J道地震记录在所选的分析时窗内，坐标(x_j，y_j)处的值为u_j。以t＝nΔt为中心的一对视倾角(p，q)的2M+1个采样点，这2M+1个采样点对应着一个J×J的协方差矩阵C：

这里u_jm＝u_j(mΔt-px_i-qy_i)，是地震道沿着t＝mΔt-px_i-qy_i的内插值。那么C2相干算法也可用公式描述如下：

$C_2(p,q)=\frac{{\stackrel{\‾}{a}}^{T}C\stackrel{\‾}{a}}{\mathrm{JTr}\left(C\right)}$

式中向量α＝(1，1，...，1)^T， $\stackrel{\‾}{a}={(\mathrm{1,1},...,1)}^T,\mathrm{Tr}\left(C\right)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{C}_{\mathrm{jj}}$ 表示协方差矩阵的迹。那么，C2相干算法得到的相干值为： $C_2=\underset{p,q}{\max }{C}_2(p,q)$

3)C3相干算法

C3相干算法借助C2相干算法中引入的协方差矩阵来实现，其分辨率较C1、C2算法更高，同时不需要层位的约束。假设λ_j(j＝1，2，...，J)是协方差矩阵C的第j个特征值，其中λ₁是最大值。C3相干算法的实现如下：

$C_3(p,q)=\frac{{\mathrm{\λ}}_1}{\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{\mathrm{\λ}}_j}$

那么C3相干算法得到的相干值为：

$C_3=\underset{p,q}{\max }{C}_3(p,q)$

地震解释过程中最大的问题仍是断层解释，M.Bahorich等于1995年提出地震相干数据的应用方法，从而使得断层的自动解释成为可能。K.J.Marfurt等于1998年提出了基于相似性的相干算法计算地震属性，该方法稳定性强，对断层的划分精度很高。我国一些学者引入该项技术后，在解释断层上取得了经验和成果，并在此基础上提出了改进。

M.Bahorich等提出的方法是在地震振幅数据体上作相干分析，振幅纵横强弱关系，信噪比都会直接影响计算出的相干体的质量，尤其是断层两侧地震信号信噪比低，在计算出的相干体上断点模糊，噪声干扰很大，而且能量也较弱。K.J.Marfurt等提出采用地震信号和地震信号Hilbert变换来计算相干体，增强稳定性，并且取得了明显效果。但是，Hilbert变换的质量在计算相干体时非常重要，该变换对噪声十分敏感，如果原始地震数据中噪声太大，很可能计算出Hilbert变换的误差很大，以致影响K.J.Marfurt等的方法应用效果。

小波变换方法在20世纪80年代兴起后，迅速用于地震资料处理领域。小波变换之所以有优于Fourier变换的特点在于它可以研究信号的局部特征，而Fourier变换是着重研究信号的整体特征。小波函数可以根据信号特征进行构造，高静怀在Morlet小波基础上构造出了适应地震资料处理的改进Morlet小波函数，王西文在改进的Morlet小波基础上构造出了适应高分辨地震资料处理的导数小波函数，这就为小波变换的时间-尺度域分析、分离信号和噪声以及分频处理带来极大的方便。小波变换计算出的Hilbert变换有较好的抗噪性，这为准确地提取瞬时特征参数奠定了基础。

M.Bahorich等和K.J.Marfurt等提出的方法都是在全频带上计算相干体，而未考虑到突出一定频带相干体的问题，尤其是突出用以解释小断层的高频段相干体信息。

发明内容

本发明要解决的技术问题在于提供一种一种基于小波变换的地震相干体计算方法。应用本发明，能提高时频分析精度、提高计算出的相干体质量，突出被忽略的小断层信息，提高地震解释精度。

本发明主要包括如下步骤：

1)对现有技术获得的地震资料进行去噪音、拓宽频谱预处理，以提高地震资料的信噪比，拓宽地震资料的有效频带宽度；

2)对步骤1所处理过的地震资料进行精细频谱分析，确定地震资料有效频谱范围；

3)对地震资料进行层位解释，并对拾取的地震解释层位进行内插、平滑等处理；

4)以地震数据、解释的地震层位为输入，利用高分辨导数小波函数，对输入地震数据的某一个地震道做小波变换，获得小波变换分频结果；其中导数小波函数如式(1)：

$\frac{\∂g}{\∂t}=(\mathrm{im}-c^{2}t)e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}$

式(1)中：m为角频率，c为常数，t为时间；

5)利用式(2)：

$\frac{1}{C_g}\underset{-\∞}{\overset{+\∞}{\∫}}S(b,a)\frac{\mathrm{da}}{a}=S\left(a\right)+\mathrm{iH}\left[S\left(b\right)\right]$

计算瞬时小波域的瞬时振幅或/和瞬时频率或/和瞬时相位；

式(2)中：Cg为常数；S(b，a)为地震道的小波变换；H[S(b)]为S(b)的Hilbert变换；

6)利用基于小波瞬时相位的相干体算法或基于小波变换的相干体算法计算瞬时相干体；其中，基于小波瞬时相位的相干体算法如式(3)：

式(3)中：maxρ_x(b，a，l，x_i，y_i)和maxρ_y(b，a，m，x_i，y_i)分别表示在延迟为l和m时，ρ_x和ρ_y的最大值。a为尺度因子定义在频带[a_i，a_i+1]，i＝1，2...N-1，频率随着序号i的增大频率增高；ρ_x地震道位置(x_i，y_i)和(x_i+1，y_i)，地震道延迟l的互相关系数；ρ_y是在地震道位置在(x_i，y_i)和(x_i，y_i+1)，地震道延迟m的互相关系数；ρ_xy是上述沿纵测线方向即地震道延迟l和沿横测线方向即地震道延迟m的3D相关系数；

基于小波变换的相干体算法如式(4)：

$C(b,a,p,q)$

$=\frac{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\{{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}{J\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\{{\left[S_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[S_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}$

式(4)中：Δb为采样时间间隔；a∈[a_i，a_i+1]；

7)利用式(5)：

${{\mathrm{\ρ}}^{*}}_{\mathrm{xy}}(b,a,l,m)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i{\mathrm{\ρ}}_{\mathrm{xy}}(b,a,l,m)$

$\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i=1$

重构瞬时相干体，获得高频或低频相干数据体；式(5)中：d_i是重构系数，d_i∈[0，1]；

8)对每一个地震道重复上述4—7步骤，得到所有地震道的相干计算结果；

9)利用所求取的瞬时相干数据进行制图，提供给地震资料解释人员，用于断层识别等研究。

本发明采用模拟地震子波的小波函数分频计算瞬时相位，提高了时频分析精度；根据M.Bahorich等方法用分频的瞬时相位计算相干体，通过重构系数，利用基于小波瞬时相位的相干体算法对一定频带相干体放大或减小，来突出特定频段相干体，易于突出被忽略的小断层信息，提高地震解释精度；直接用小波分频变换的实、虚部，根据基于小波变换的相干体算法，克服了直接计算Hilbert变换带来的误差，提高了相干体计算质量。

本发明具体实现原理如下：

1)理论基础

地震信号S(t)的小波变换可以表示为：

$S_R(b,a)=\frac{1}{a}{\∫}_{-\∞}^{+\∞}S\left(t\right)g_R\left(\frac{t-b}{a}\right)\mathrm{dt}$

$S_I(b,a)=\frac{1}{a}{\∫}_{-\∞}^{+\∞}S\left(t\right)g_I\left(\frac{t-b}{a}\right)\mathrm{dt}$

式中：b为时间因子，且t，b∈R；a为尺度因子，且a∈R\{0}，S_R(b，a)为小波变换的实部，S_I(b，a)为小波变换的虚部，

是小波函数。

如果小波函数满足：

g(t)∈L1(R，dt)∩L2(R，dt)

g(ω)∈L1(R\{0}，dω/ω)∩L2{R\{0}，dω/ω}

Cg≠0

则

$\frac{1}{C_g}\underset{-\∞}{\overset{+\∞}{\∫}}S(b,a)\frac{\mathrm{da}}{a}=S\left(a\right)+\mathrm{iH}\left[S\left(b\right)\right]$

式中，C_g为常数；S(b，a)为地震道的小波变换；H[S(b)]为S(b)的Hilbert变换；

可分别对应不同尺度因子a，即相当不同频率下定义瞬时振幅，瞬时相位，瞬时频率。以及不同尺度因子a的阻尼瞬时频率。

小波函数选取如下式所示改进的Morlet小波，该小波函数用常数c控制高斯函数来调制小波函数g(t)的特征，可最大程度地模拟地震子波。

$g\left(t\right)=e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}$

式中，c为常数；m为角频率。

为了更好的模拟高频率地震信号，采用如下的高分辨导数小波函数，主要用于高频段地震信号的分解。

$\frac{\∂g}{\∂t}=(\mathrm{im}-c^{2}t)e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}$

式中：m为角频率，c为常数，t为时间。

2)基于瞬时相位的相干体算法

据此：如果在地震道位置(x_i，y_i)和(x_i+1，y_i)，地震道延迟l的互相关系数ρ_x，在地震道位置在(x_i，y_i)和(x_i，y_i+1)，地震道延迟m的互相关系数ρ_y，那么可将上述沿纵测线方向(地震道延迟l)和沿横测线方向(地震道延迟m)的3D相关系数ρ_xy定义为：

其中：maxρ_x(b，a，l，x_i，y_i)和maxρ_y(b，a，m，x_i，y_i)分别表示在延迟为l和m时，ρ_x和ρ_y的最大值。a为尺度因子定义在频带[a_i，a_i+1]，i＝1，2...N-1，频率随着序号i的增大频率增高。

3)基于小波变换的相干体算法

相似算法是稳定的算法，其先定义一个以分析点为中心的J道椭圆或矩形分析时窗。若取分析点坐标(x，y)为局部中心，则定义相似系数为σ(b，a，p，q)。式中给出的相似估计，但是对一些小的相干地震同相轴是不稳定的。因此，在一个高度为2w(w＝KΔb)ms的垂直时窗上计算一个平均相似系数，该平均相似系数，即相干估计值定义为：

$C(b,a,p,q)$

$=\frac{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\{{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}{J\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\{{\left[S_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[S_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}$

式中，Δb为采样时间间隔；a∈[a_i，a_i+1]。

4)重构相干体

通常依据基于瞬时相位的相干体算法(简称本发明相干算法1)所计算的相干体仅是在一定频带下的相干体，因此可以重构相干体，其定义为：

${{\mathrm{\ρ}}^{*}}_{\mathrm{xy}}(b,a,l,m)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i{\mathrm{\ρ}}_{\mathrm{xy}}(b,a,l,m)$

$\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i=1$

式中，d_i是重构系数，d_i∈[0，1]。

根据基于小波变换的相干体算法(简称本发明相干算法2)及其S_R(b，a)和S_I(b，a)的定义，用计算出的相干体，其重构相干体定义为：

$C^{*}(b,a,p,q)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_{i}C[b,a,p,q]$

$\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i=1$

附图说明

图1(a)是实施例1基于小波瞬时相位的相干体算法计算得到相干切片，

图1(b)是实施例1基于小波变换的相干体算法计算得到相干切片，

图2(a)是实施例2馆陶组底界高频相干，

图2(b)是实施例2馆陶组底界高频相干与断层叠合

图3(a)是实施例2馆陶组底界低频相干，

图3(b)是实施例2馆陶组底界低频相干与断层叠合，

图4(a)是实施例2原馆陶组底界构造图，

图4(b)是实施例2新馆陶组底界构造图，

图5(a)是实施例3周清庄油田的3D地震数据中的一条剖面，

图5(b)是实施例3该地震剖面采用K.J.Marfurt方法计算相干剖面，

图5(c)是实施例3该地震剖面采用K.J.Marfurt方法利用小波计算的相干剖面，

图5(d)是实施例3该地震剖面采用本发明方法计算的高频瞬时相位剖面，

图5(e)是实施例3该地震剖面采用本发明方法计算出低频瞬时相位剖面，

图5(f)是实施例3该地震剖面采用本发明方法重构瞬时相位剖面。

具体实施方式

实施例1

本发明包括如下步骤：

对利用现有技术，通过野外高分辨率地震采集设备所获得的原始地震数据进行去噪音、拓宽频谱预处理，以提高地震资料的信噪比，拓宽地震资料的有效频带宽度；

2、对步骤1所处理过的地震资料进行精细频谱分析，确定地震资料有效频谱范围；

3、对地震资料进行层位解释，并对拾取的地震解释层位进行内插、平滑等处理；

4、以地震数据、解释的地震层位为输入，利用高分辨导数小波函数，对输入地震数据的某一个地震道做小波变换，获得小波变换分频结果；其中导数小波函数如式(1)：

$\frac{\∂g}{\∂t}=(\mathrm{im}-c^{2}t)e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}$

式(1)中：m为角频率，c为常数，t为时间；

由式(2)：

$\frac{1}{C_g}\underset{-\∞}{\overset{+\∞}{\∫}}S(b,a)\frac{\mathrm{da}}{a}=S\left(a\right)+\mathrm{iH}\left[S\left(b\right)\right]$

计算瞬时小波域的瞬时振幅、瞬时频率或瞬时相位；

式(2)中：Cg为常数；S(b，a)地震道的小波变换；H[S(b)]为S(b)的Hilbert变换

利用基于小波瞬时相位的相干体算法或基于小波变换的相干体算法计算瞬时相干体；其中，基于小波瞬时相位的相干体算法如式(3)：

式(3)中：maxρ_x(b，a，l，x_i，y_i)和maxρ_y(b，a，m，x_i，y_i)分别表示在延迟为l和m时，ρ_x和ρ_y的最大值；a为尺度因子定义在频带[a_i，a_i+1]，i＝1，2...N-1，频率随着序号i的增大频率增高；ρx是在地震道位置(x_i，y_i)和(x_i+1，y_i)，地震道延迟l的互相关系数；ρy是在地震道位置(x_i，y_i)和(x_i，y_i+1)，地震道延迟m的互相关系数；ρxy是上述沿纵测线方向和沿横测线方向的3D相关系数；即ρxy是地震道延迟I和地震道延迟m的3D相关系数；

基于小波变换的相干体算法如式(4)：

$C(b,a,p,q)$

$=\frac{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\{{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}{J\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\{{\left[S_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[S_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}$

式(4)中：Δb为采样时间间隔；a∈[a_i，a_i+1]；

利用式(5)：

${{\mathrm{\ρ}}^{*}}_{\mathrm{xy}}(b,a,l,m)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i{\mathrm{\ρ}}_{\mathrm{xy}}(b,a,l,m)$

$\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i=1$

重构造瞬时相干体，获得高频或低频相干数据体；式(5)中：d_i是重构系数，d_i∈[0，1]；

对每一个地震道重复上述4—7步骤，得到所有地震道的相干计算结果；

利用所求取的瞬时相干数据进行制图，提供给地震资料解释人员，用于断层识别等研究。

图1(a)为沿层采用本发明相干算法1即基于小波瞬时相位的相干体算法所获得的高频相干切片；图1(b)为沿层采用本发明相干算法2即基于小波变换的相干体算法所获得的高频相干切片。由图可以得出如下结论，该发明计算相干切片具有识别断层或不均质体分辨率高的特点。

实施例2 图2(a)—图4(b)为本发明在大港油田唐东—马南地区的应用情况。其步骤与实施例1相同。

图2(a)是实施例2馆陶组底界高频相干，该结果可以很好的识别微小断裂；图中非常清楚地显示出主要控制油田构造断裂有两组，一组以F1为代表的北北东向断裂；一组是以F2为代表的西西东向断裂。图中，A处是一组4-5条断层，北北东向小断层交汇在F2断裂上，而且断裂展布清晰，分辩率高；B处是近东西向断裂交汇在F2断裂上，交汇点清晰；C处是北北东向断裂交汇在F2断裂上，断裂点清晰；D处是一组3条雁行式排列北北东向断裂，分辩很清楚；E处是一条北北东向断裂与一条北东向断裂交汇点。

图2(b)实施例2馆陶组底界高频相干与构造图的叠合，叠合结果说明该相干结果与解释结果一致，符合地质认识，而且左侧三维数据解释断裂与右图相干体切片断裂线非常吻合。尤其是A，B，C，D，E处断裂交汇点清晰。

图3(a)实施例2馆陶组底界低频相干，该结果可以很好的识别较大断裂，图中F1和F2断裂非常清楚。但是，A处北北东向这组小断裂较高频相干体分辩不清，B，C处断裂交汇处与高频相干切片相比模糊不清，D处这组雁行排列断裂间距较大，低频相干体仍能分辨清楚，E处断裂交汇点模糊不清。

图3(b)实施例2馆陶组底界低频相干与构造图的叠合，叠合结果说明该相干结果与解释结果一致，大断裂格局非常清楚符合地质认识。

图4(a)和图4(b)分别是原沙三I底构造图和是原沙三IV底构造图。两图相比，断裂特征相同，说明断裂是继承性的，这点与实际不符。图中出现近似直角弧形断裂，在塑性地层中不可能出现这种断裂，说明断裂解释有误。图中很难分析出构造受力方向和主要控制构造断裂体系，说明断裂组合有问题。因此，这种构造图无法为精细储层预测研究提供正确的构造信息；相反，还会引到出错误结论。

实施例3 图5(a)—图5(f)为本发明在大港油田周清庄应用的结果。其步骤与实施例1相同。图5(a)为实施例3周清庄油田的3D地震数据中的一条剖面，图5(a)中所示的S31是一层生物灰岩储层，S31至S34是一套河流相的砂泥岩互层沉积，其中砂岩也是非常好的油气储层。二套油气储层受小断裂控制。该油田已处于开发阶段的中后期，研究受大断裂控制的小断裂对油气储层控制作用十分重要。研究区断裂发育，主要发育—组控制箕状断陷正断层，从地震剖面信息上很难解释出图中所示这组发育的正断层。图中所示断裂是经过精细构造研究后的解释成果。

图5(b)实施例3该地震剖面采用K.J.Marfurt等的方法用Hilbert变换计算相干剖面，由于地震信号信噪比较低，这使Hilbert变换误差增大，以致在计算的相干剖面中断层断点非常模糊，难以解释出图中所示这组发育的正断层。

图5(c)为实施例3该地震剖面采用K.J.Marfurt等的方法利用小波计算的相干剖面，小波变换方法得到Hilbert变换，再计算相干剖面即基于小波变换的相干体算法。很明显，图5(c)中断层断点较明显。但是，该图在用于小断层解释时，分辨率太低。

图5(d)为实施例3该地震剖面采用本发明方法计算的高频瞬时相位剖面，首先利用模拟地震子波的高分辨导数小波函数计算出高频瞬时相位。很明显，图5(d)中所示断层断点清晰，尤其是在CDP1-CDP46；时间2100ms-2500ms之间小断层断点非常清楚。这为精细构造解释和小断层展布的研究提供了基础资料。

图5(e)为实施例3该地震剖面采用本发明方法计算出低频瞬时相位剖面，首先利用本发明方法模拟地震子波的小波函数计算出低频瞬时相位。很明显，图5(e)中所示在CDP1-CDP46；时间2100ms-2900ms之间，控制箕状凹陷大断层断点非常清楚。

图5(f)为实施例3该地震剖面采用本发明方法重构瞬时相位剖面，取高频重构系数d₁＝0.7；低频重构系数d₂＝0.3得到重构相干剖面。图5(f)中所示高频信息增强，一些小断层显示更明显，很容易解释出图5(f)中所示这组发育的正断层。但是，在高频信息放大的同时，相应的高频噪声也变大。

总之，本发明给出基于小波变换的二种地震数据的相干体算法。

基于小波瞬时相位的相干体算法是用模拟地震子波的小波函数或高分辨导数小波函数分频计算瞬时相位，根据分频瞬时相位计算相干体，再进行重构。该算法可以提高断层解释的分辨率，其应用效果是明显的。

基于小波变换的相干体算法是用小波变换得到的实、虚部即Hilbert变换，直接计算相干体。该算法可以克服直接计算Hilbert变换带来的误差，提高计算出的相干体质量。基于小波瞬时相位的相干体算法可以得到高、低频瞬时相位信息，有利于精细构造解释，而且计算出的相干体在小断层分辨程度方面也较基于小波变换的相干体算法效果明显。根据本发明的研究实例，采用分频重构的相干体，易于突出被忽略的小断层信息。实际应用也证明基于小波瞬时相位的相干体算法是十分有效的。

Claims (1)

1、一种基于小波变换的地震相干体计算方法，包括如下步骤：

1)对现有技术获得的地震资料进行去噪音、拓宽频谱预处理，以提高地震资料的信噪比，拓宽地震资料的有效频带宽度；

2)对步骤1所处理过的地震资料进行精细频谱分析，确定地震资料有效频谱范围；

3)对地震资料进行层位解释，并对拾取的地震解释层位进行内插、平滑等处理；

4)以地震数据、解释的地震层位为输入，利用高分辨导数小波函数，对输入地震数据的某一个地震道做小波变换，获得小波变换分频结果；其中导数小波函数如式(1)：

$\frac{\∂g}{\∂t}=(\mathrm{im}-c^{2}t)e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}$

式(1)中：m为角频率，c为常数，t为时间；

5)利用式(2)：

$\frac{1}{C_g}\underset{-\∞}{\overset{+\∞}{\∫}}S(b,a)\frac{\mathrm{da}}{a}=S\left(a\right)+\mathrm{iH}\left[S\left(b\right)\right]$

计算瞬时小波域的瞬时振幅或/和瞬时频率或/和瞬时相位；

式(2)中：Cg为常数；S(b，a)地震到的小波变换；H[S(b)]为S(b)的Hilbert变换；

6)利用基于小波瞬时相位的相干体算法或基于小波变换的相干体算法计算瞬时相干体；

其中，基于小波瞬时相位的相干体算法如式(3)：

式(3)中：max ρ_x(b，a，l，x_i，y_i)和max ρ_y(b，a，m，x_i，y_i)分别表示在延迟为l和m时，ρ_x和ρ_y的最大值；a为尺度因子定义在频带[a_i，a_i+1]，i＝1，2...N-1，频率随着序号i的增大频率增高；ρx是在地震道位置(x_i，y_i)和(x_i+1，y_i)，地震道延迟l的互相关系数；ρy是在地震道位置(x_i，y_i)和(x_i，y_i+1)，地震道延迟m的互相关系数；ρ xy是上述沿纵测线方向和沿横测线方向的3D相关系数；即ρ xy是地震道延迟1和地震道延迟m的3D相关系数；

基于小波变换的相干体算法如式(4)：

$C(b,a,p,q)$

$=\frac{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\{{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}{J\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\{{\left[S_R(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2+{\left[S_I(b+\mathrm{k\Δb}-{\mathrm{px}}_j-{\mathrm{qy}}_j,a,x_j,y_j)\right]}^2\}}$

式(4)中：Δb为采样时间间隔；a∈[a_i，a_i+1]；

7)利用式(5)：

${{\mathrm{\ρ}}^{*}}_{\mathrm{xy}}(b,a,l,m)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i{\mathrm{\ρ}}_{\mathrm{xy}}(b,a,l,m)$

$\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{d}_i=1$

重构瞬时相干体，获得高频或低频相干数据体；式(5)中：d_i是重构系数，d_i∈[0，1]；

8)对每一个地震道重复上述4—7步骤，得到所有地震道的相干计算结果；

9)利用所求取的瞬时相干数据进行制图，提供给地震资料解释人员，用于断层识别等研究。

Images