一种提高地震数据分辨率增强有效弱信号能量的方法

摘要

本发明涉及一种地震数据高分辨率有效弱信号能量增强的方法，所述方法需将叠后或偏移后的地震数据逐道进行处理，对单道地震数据执行以下步骤：1)对单道地震数据进行Hilbert变换，获取单道地震数据的包络、时移量及瞬时属性，联合井震标定提取的平均子波，建立Morlet子波库；2)基于第1)步建立的子波库，利用匹配追踪算法分解出组成单道地震数据的各有效子波；3)改进谱白化算法，对分解得到的各有效子波进行自适应拓频和能量补偿处理；4)将第3)步处理后的有效子波叠加，得到提高分辨率增强有效弱信号能量后的单道地震数据。该方法综合测井、地震资料，将谱白化算法改进并与匹配追踪算法结合，达到拓宽资料频带增强有效弱信号能量的效果。

说明书

技术领域

本发明属于油气勘探中地震资料分析与处理领域，涉及一种地震数据处理 方法，具体地说，是关于一种匹配追踪多子波分解与改进谱白化相结合的提高 地震数据分辨率增强有效弱信号能量的方法。

背景技术

高分辨率地震数据是油藏描述的先决条件，精细油气描述、空间非均质性 的判断以及监测油气田开发过程中储层流体变换等都需要高分率地震数据作为 基础。

传统地震勘探中，地震记录可以看成是子波与地层反射系数的褶积，常规 提高分辨率的处理方法有反褶积、谱白化，传统地震模型的建立基于若干基本 假设，其中之一就是假设子波是平稳不变的，在此模型的基础上，Rayleigh (1945)通过对两个相同反射系数和零相位子波的研究，定义地震纵向分辨率 的极限为1/4波长；Ricker(1953)将零相位子波混叠，通过研究其混叠后波形的 二阶倒数来定义纵向分辨率的极限，Widess(1973)在Rayleigh的基础上取极 性相反的反射系数和零相位子波进行研究，利用时间一阶导数将纵向分辨率极 限定义为1/8波长；Schoenberger(1974)在固定振幅谱的基础上通过改变子波 来研究纵向分辨率，上述作者对纵向分辨率的研究主要集中在对子波研究的基 础上，此时提高分辨率的关键在于压缩地震子波。Kallweit和Wood(1982)在 同向系数和反向反射系数的楔状模型的基础上，揭示了调谐振幅和地层厚度之 间的变化规律；Juhlin(1993)和Liu(2003)讨论了1/8波长薄层和AVO的关 系，至此提高分辨率的研究正在从单一的子波压缩，逐渐发展到子波压缩、薄 层振幅联合考虑的情况。

实际上，地震子波在传播过程中，在遇到反射界面发生反射时，会形成振 幅和形状随时间变化的子波，由于不同深度的反射对应的传播路程长短不一， 子波变化的程度也不尽相同，当子波传回地面被检波器接受并记录下来时，形 成的地震道与传统的地震道褶积模型相差较大。由于地震反射的不同特征(频 率、振幅、相位)可以用于表征不同的地下地层岩性、物性、及含流体性，因 此出现了多子波分解及重构的方法。运用多子波分解重构技术能够更准确的了 解地面接收的地震信号，为地震信号处理提供了一个新的思路。

Mallat等(1993)首次提出了匹配追踪算法，其基于Gabor时频原子库， 能够从冗余原子库中选取时频原子将信号进行多子波分解。匹配追踪能够将地 震道分解成主频不同、时移不同、相位不同以及尺度不同的子波集合，从而实 现地震道在时间、频率内的精确分解。然后可按照需求选择不同频率的单一地 震子波或某一频段的子波集，从而进一步重构生成新的数据体。

尽管匹配追踪分解能够得到较好的信号表示，但是该方法计算量过于庞大， 因此，在保证其计算精度的同时，若能进一步降低匹配追踪算法的计算量是目 前面临的急需解决的问题之一。单一应用匹配追踪分解虽然能够将地震资料中 的信号和噪声分离，通过剔除噪声残差来提高地震资料的信噪比，但其并未考 虑地震资料分辨率的提高以及有效弱信号能量增强的问题，传统谱白化虽然是 一种展宽频谱的基本方法，但其并不能针对子波的特点，对子波振幅进行自适 应的能量增强。

发明内容

基于现有技术中存在的上述问题，本发明的目的在于提供一种提高地震数 据分辨率增强有效弱信号能量的方法。

为实现上述发明目的，本发明采用了以下技术方案：一种提高地震数据分 辨率增强有效弱信号能量的方法，所述方法需将叠后或偏移后的地震数据逐道 进行处理，每个单道地震数据的执行，其特征在于，包括以下步骤：

1)在对单道地震数据进行Hilbert变换的基础上，获取单道地震数据的包 络、频率、相位、时移属性，联合井震标定过程中提取的平均子波的宽度，来 共同建立Morlet子波库；

2)在第1)步建立的Morlet子波库的基础上，利用匹配追踪多子波分解 算法分解出M个组成单道地震数据的有效子波；

3)采用改进的谱白化算法，对第2)步分解得到的各有效子波进行自适 应拓频和能量补偿处理；

4)将第3)步处理后的有效子波进行叠加，得到提高分辨率增强有效弱 信号能量后的单道地震数据。

所述步骤1)中，Morlet子波库的建立步骤包括：

①输入N道叠后或偏移后的地震数据，对第j道地震数据道进行Hilbert变 换，构建复地震道p_j(t)

$p_j\left(t\right)=s_j\left(t\right)+i{\tilde{s}}_j\left(t\right)---\left(1\right)$

在等式(1)中，s_j(t)为第j道地震数据(j＝1，2…N)，为s_j(t)的Hilbert 变换结果；

②通过p_j(t)，计算第j道地震信号包络的最大值及包络最大值所对应的时 间

$a_0=\underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(2\right)$

$u_0=\arg \underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(3\right)$

等式(2)中a₀为第j道地震信号包络的最大值，等式(3)中u₀为第j道 地震信号包络最大值所对应的时间；

③通过p_j(t)，获取第j道地震信号的瞬时频率和瞬时相位，同时记录下u₀处的瞬时频率和瞬时相位

${\omega }_0=\frac{d}{\mathrm{dt}}\left(\arctan \frac{{\tilde{s}}_j\left(t\right)}{s_j\left(t\right)}\right),t=u_0---\left(5\right)$

等式(4)中为瞬时相位，等式(5)中ω₀为瞬时频率；

④结合测井资料，进行井震标定，并提取平均子波，确定平均子波的宽度 σ₀，并以σ₀为中心确定σ的搜多范围，其可避免大量公式计算，提高计算效率；

⑤以先验信息点(u₀，φ₀，ω₀，σ₀)为中心，创建子波中心位置u、瞬时频率 ω、瞬时相位φ以及子波宽度σ的参数集R，u的变化步长为时间采样间隔， 瞬时频率ω的变化步长为频率采样间隔，瞬时相位φ的变化步长为5度，子波 宽度σ的变化步长为0.5；

⑥利用第⑤步构建的参数集R，构建Morlet子波库，其中Morlet小波的表 达式为

将参数集R中的不同参数带入(6)式，即可获得不同的Morlet子波。

所述步骤2)中，匹配追踪多子波分解的具体步骤包括：

在步骤1)获得的Morlet子波库的基础上，第j道地震信号s_j通过正交分 解可以表示为：

s_j＝＜s_j，m_r0＞m_r0+R⁽¹⁾s_j     (7)

等式(7)中R⁽¹⁾s_j为信号s_j分解得到m_r0子波后的残差，且R⁽⁰⁾s_j＝s_j，由 于m_r0和R⁽¹⁾s_j正交，则有：

||s_j||²＝|＜s_j，m_r0＞|²+||R⁽¹⁾s_j||²     (8)

假设已迭代m次(m≥0)，其残差为R^(m)s_j，在子波库中按照内积 ＜R^(m)，m_rm＞最大准则来选取最佳小波基m_rm，则有：

$m_{\mathrm{rm}}\left(t\right)=\arg \underset{m_{\mathrm{rm}}\in D}{\max }\frac{<R^{\left(m\right)}{s}_j,m_{\mathrm{rm}}>}{\left|\right|m_{\mathrm{rm}}\left|\right|}---\left(9\right)$

将残差进行分解，可得：

R^(m)s_j＝＜R^(m)s_j，m_rm＞m_rm+R^(m+1)s_j       (10)

第j道地震信号s_j共进行M次分解后，该地震道信号可以表示为：

$s_j={\Sigma }_{m=0}^{M-1}{a}_m{m}_{\mathrm{rm}}+R^{\left(M\right)}{s}_j---\left(11\right)$

等式(11)中a_m为m_rm(t)的振幅，R^(M)s_j是残差，当进行M次迭代之后 R^(M)s_j可视为噪声，振幅a_m为：

$a_m=\frac{|<R^m{s}_j,m_{\mathrm{rm}}>|}{{\left|\right|m_{\mathrm{rm}}\left|\right|}^2}---\left(12\right)$

等式(12)中是小波基的标准化，令g_rm＝a_m× m_rm为每次分解后得到的最佳匹配子波。

所述步骤3)中，对步骤2)中分解得到的各有效子波进行改进谱白化计 算，具体步骤包括：

①利用步骤1)中构建的复地震道p_j(t)来获取s_j(t)地震道的瞬时包络e_j(t)

$e_j\left(t\right)=\underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(13\right)$

②通过s_j(t)的瞬时包络来求取该地震道第m个最佳匹配子波的白化滤波器

$f_m\left(t\right)=\frac{v}{e_j\left(t\right)+{ϵ}_j\left(m\right)v}---\left(14\right)$

等式(14)中v为地震信号s_j(t)的振幅包络e_j(t)的最大值，ε_j(m)为第j道 的第m个子波的白噪因子，ε_j(m)由组成地震信号的各最佳匹配子波的振幅强 弱决定，则有：

ε_j(m)＝ε_base+ε_wave×E_j(m)      (15)

等式(15)中ε_base和ε_wave为常量，E_j(m)为第j道的第m个子波的能量比 重，则E_j(m)为：

$E_j\left(m\right)=\frac{a_m}{{\Sigma }_0^{M-1}{a}_m}---\left(16\right)$

③计算子波在改进谱白化后瞬时振幅

通过步骤3)可以将s_j(t)分解变为M个最佳匹配子波g_rm的叠加和，对于 任一最佳匹配子波g_rm(t)，m＝0，1，2...M，其振幅和瞬时频率分别为a_m(t)和 ω_m(t)，子波白化后的瞬时振幅为：

${\tilde{a}}_m\left(t\right)=a_m\left(t\right)f_m\left(t\right)---\left(17\right)$

所述步骤4)中，结合步骤3)分解得到的最佳小波基m_rm和步骤4)算得 的补偿能量后的瞬时振幅，最终获得第j道地震信号处理后的结果S_j：

$S_j={\Sigma }_{m=0}^{M-1}{\tilde{a}}_m{m}_{\mathrm{rm}}---\left(18\right)$

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于采用匹 配追踪与改进谱白化相结合的方法对地震道的时间域信号s_j(t)进行处理，因 此本发明能够在剔除噪声残差提高地震资料信噪比的同时，自适应的增强有效 弱信号的能量，提高资料的分辨率。2、本发明改变了传统获取子波尺度参数σ 的方法，从井震标定获取的平均子波出发，来确定平均子波宽度σ₀，创建了以 σ₀为中心的σ参数集，因此本发明能够充分结合测井数据，进一步减少匹配追 踪算法中迭代寻找σ₀的计算量，从而提高匹配追踪的计算效率。3、本发明结 合子波能量占地震道总能量的比重，改进了白化滤波器的表达式，因此本发明 能够对匹配追踪分解得到的每个子波进行改进谱白化操作，最终实现能量较弱 的有效子波其能量得到补偿的效果。基于以上优点，本发明可以广泛应用于地 震数据的处理中，处理后的地震资料其分辨率有所提高，且有效弱信号的能量 有所加强。

附图说明

图1是本发明的流程示意图；

图2(a)是本发明所涉及的理论含噪地震道；

图2(b)是本发明对理论含噪模型道进行匹配追踪分解后筛选出来的有效 子波；

图2(c)是本发明将匹配追踪分解后筛选出来的有效子波直接进行重构的 结果；

图2(d)是本发明将匹配追踪分解后筛选出来的有效子波先进行改进谱白 化处理后进行重构的结果；

图3(a)是本发明对图2(a)所示的地震道进行时频分析的结果；

图3(b)是本发明对图2(c)所示的地震道进行时频分析的结果；

图3(c)是本发明对图2(d)所示的地震道进行时频分析的结果；

图4(a)是本发明所要处理的实际含噪地震剖面；

图4(b)是本发明对图4(a)应用地震数据分辨率增强有效弱信号能量 处理后的地震剖面；

图5(a)是本发明对图4(a)所示的实际含噪地震剖面求取的振幅谱示意 图；

图5(b)是本发明对图4(b)中应用地震数据高分辨率有效弱信号能量 增强方法处理后的地震剖面求取的振幅谱示意图；

图6(a)图6(b)图6(c)是本发明对图4(a)中实际含噪地震剖面的 分频结果示意图；

图7(a)图7(b)图7(c)是本发明对图4(b)中应用地震数据高分辨 率有效弱信号能量增强方法处理后的地震剖面分频结果示意图。

具体实施方式

下面结合附图和实例对本发明进行详细的描述。

本发明从地震道的组合形式出发，野外采集得到的地震记录，经过预处理 后变为叠后地震数据或偏移数据，根据本发明方法需逐道对地震数据进行处理。

图1示出本发明提高地震数据分辨率增强有效弱信号能量的方法流程，包 括以下步骤：

1)对单道地震数据进行Hilbert变换，构建复地震道；从复地震道中获取 单道地震数据的包络、频率、相位、时移属性；利用已知的测井曲线进行井震 标定，并提取井震标定过程中平均子波的宽度；结合上述信息来共同建立 Morlet子波库：

①输入N道叠后或偏移后的地震数据，对第j道地震数据道进行Hilbert变 换，构建复地震道p_j(t)

$p_j\left(t\right)=s_j\left(t\right)+i{\tilde{s}}_j\left(t\right)---\left(1\right)$

在等式(1)中，s_j(t)为第j道地震数据(j＝1，2…N)，为s_j(t)的Hilbert 变换结果；

②通过p_j(t)，计算第j道地震信号包络的最大值及包络最大值所对应的时 间

$a_0=\underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(2\right)$

$u_0=\arg \underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(3\right)$

等式(2)中a₀为第j道地震信号包络的最大值，等式(3)中u₀为第j道 地震信号包络最大值所对应的时间；

③通过p_j(t)，获取第j道地震信号的瞬时频率和瞬时相位，同时记录下u₀处的瞬时频率和瞬时相位

${\omega }_0=\frac{d}{\mathrm{dt}}\left(\arctan \frac{{\tilde{s}}_j\left(t\right)}{s_j\left(t\right)}\right),t=u_0---\left(5\right)$

等式(4)中为瞬时相位，等式(5)中ω₀为瞬时频率；

④结合已有的测井资料，进行井震标定，并提取平均子波，从而可以根据 地震数据的采样率来进一步确定平均子波的宽度σ₀，由于地震信号中各子波的 宽度差距不大，若确定平均子波宽度σ₀，并以σ₀为中心确定σ的搜多范围，则 可以避免大量公式计算，最终可以提高计算效率；

⑤以先验信息点(u₀，φ₀，ω₀，σ₀)为中心，创建子波中心位置u、瞬时频率 ω、瞬时相位φ以及子波宽度σ的参数集R，并根据实际特点确定各参数的变 化步长，本步骤取u的变化步长为时间采样间隔，瞬时频率ω的变化步长为频 率采样间隔，瞬时相位φ的变化步长为5度，子波宽度σ的变化步长为0.5；

⑥利用第⑤步构建的参数集R，构建Morlet子波库，其中Morlet小波的表 达式为

参数集将R中的不同参数带入(6)式，即可获得不同的Morlet子波。

2)在步骤1)建立的Morlet子波库的基础上，利用匹配追踪多子波分解 算法分解出M个组成单道地震数据的有效子波

匹配追踪算法是一个从冗余子波库中迭代寻找最佳匹配子波的过程。每次 迭代可以在子波库中获取一个与地震道信号s_j或地震信号与已分解子波的残差 信号的最佳小波基m_rm，然后从信号中剔除最佳匹配子波，以获得此次迭代计 算后的残差信号；在进行M次迭代后，则地震信号可以由分解得到的M个最 佳小波基的线性组合得到，此处通过对资料的前期分析，以分解M次作为分解 结束的条件。

在步骤1)获得的Morlet子波库的基础上，第j道地震信号sj通过正交分 解可以表示为：

s_j＝＜s_j，m_r0＞m_r0+R⁽¹⁾s_j      (7)

等式(7)中R⁽¹⁾s_j为信号s_j分解得到m_r0子波后的残差，且R⁽⁰⁾s_j＝s_j，由 于m_r0和R⁽¹⁾s_j正交，则有：

||s_j||²＝|＜s_j，m_r0＞|²+||R⁽¹⁾s_j||²       (8)

若要使残差R⁽¹⁾s_j达到最小，则从子波库中选取的子波m_r0需要使内积＜ s_j，m_r0＞最大，假设已迭代m次(m≥0)，其残差为R^(m)s_j，在子波库中按照 内积＜R^(m)，m_rm＞最大准则来选取最佳小波基m_rm，则有：

$m_{\mathrm{rm}}\left(t\right)=\arg \underset{m_{\mathrm{rm}}\in D}{\max }\frac{<R^{\left(m\right)}{s}_j,m_{\mathrm{rm}}>}{\left|\right|m_{\mathrm{rm}}\left|\right|}---\left(9\right)$

将残差进行分解，可得：

R^(m)s_j＝＜R^(m)s_j，m_rm＞m_rm+R^(m+1)s_j      (10)

当第j道地震信号s_j进行M次分解后，则该地震道信号可以表示为：

$s_j={\Sigma }_{m=0}^{M-1}{a}_m{m}_{\mathrm{rm}}+R^{\left(M\right)}{s}_j---\left(11\right)$

等式(11)中a_m为m_rm(t)的振幅，R^(M)s_j是残差，当进行M次迭代之后， R^(M)s_j可视为噪声，振幅a_m为：

$a_m=\frac{|<R^m{s}_j,m_{\mathrm{rm}}>|}{{\left|\right|m_{\mathrm{rm}}\left|\right|}^2}---\left(12\right)$

等式(12)中是小波基的标准化，令g_rm＝a_m× m_rm为每次分解得到的最佳匹配子波。

3)采用改进的谱白化算法，对步骤2)分解得到的各有效子波进行自适应 拓频和能量补偿处理

①利用步骤1)中构建的复地震道p_j(t)来获取s_j(t)地震道的瞬时包络e_j(t)

$e_j\left(t\right)=\underset{t}{\max }\left(\sqrt{s_j{\left(t\right)}^2+{\tilde{s}}_j{\left(t\right)}^2}\right)---\left(13\right)$

②通过s_j(t)的瞬时包络来求取该地震道第m个最佳匹配子波的白化滤波器

$f_m\left(t\right)=\frac{v}{e_j\left(t\right)+{ϵ}_j\left(m\right)v}---\left(14\right)$

等式(14)中v为地震信号s_j(t)振幅包络e_j(t)的最大值，ε_j(m)为第j道的 第m个子波的白噪因子，ε_j(m)由组成地震信号的各最佳匹配子波的振幅强弱 决定，则有：

ε_j(m)＝ε_base+ε_wave×E_j(m)      (15)

等式(15)中ε_base和ε_wave为常量，E_j(m)为第j道的第m个子波的能量比 重，则E_j(m)为：

$E_j\left(m\right)=\frac{a_m}{{\Sigma }_0^{M-1}{a}_m}---\left(16\right)$

③计算子波在改进谱白化后瞬时振幅

通过步骤3)可以将s_j(t)分解变为M个最佳匹配子波g_rm的叠加和，对于 任一最佳匹配子波g_rm(t)，m＝0，1，2…M，其振幅和瞬时频率分别为a_m(t)和 ω_m(t)，a_m(t)在子波延续时间点上与白化滤波器相乘，从而得到子波白化后的 瞬时振幅为：

${\tilde{a}}_m\left(t\right)=a_m\left(t\right)f_m\left(t\right)---\left(17\right)$

4)结合步骤2)分解得到的最佳小波基m_rm和步骤3)算得的补偿能量后的瞬 时振幅，最终获得第j道地震信号处理后的结果S_j：

$S_j{=\Sigma }_{m=0}^{M-1}{\tilde{a}}_m{m}_{\mathrm{rm}}---\left(18\right)$

图2给出了单道合成地震记录算例图，其中图2(a)是含噪合成地震记录； (b)是利用匹配追踪多子波分解所获得的8个有效子波；(c)是将图(b)的 8个有效子波直接叠加后的结果；(d)是对(b)中的8个有效子波先进行改进 谱白化计算然后进行叠加的结果。对比图2(a)和(c)可见，匹配追踪多子 波分解可以有效剔除白噪声，将分解得到的有效子波进行叠加其结果较原始数 据信噪比有显著的提高，对比图2(c)和(d)可见，改进谱白化算法可以有 效加强有效弱信号的能量。

图3给出了图2(a)、图2(c)、图2(d)的时频分析结果，其中图3(a) 是图2(a)原始含噪合成地震记录所对应的时频分析结果，由图可见原始地震 道中噪声对时频分析的影响较为严重，背景中含有杂质噪声，部分有效弱信号 时频特征显示不清；(b)是图2(c)所对应的时频分析结果，由图可见匹配追 踪多子波分解对噪声有良好的剔除效果，从时频分析图中可以看出其背景噪声 得到了较好的剔除；(c)是图2(d)所对应的时频分析结果，由图可见改进谱 白化算法可以有效加强有效弱信号的能量。

图4给出了应用地震数据高分辨率有效弱信号能量增强方法处理前后的实 际地震数据对比图，图4(a)是原始实际地震资料；(b)是应用地震数据高分 辨率有效弱信号能量增强方法处理后的地震记录，对比图4(a)和(b)可见 地震数据高分辨率有效弱信号能量增强方法处理后的地震记录，其纵横向分辨 率均有所提高，有效弱信号能量得以加强，其同相轴变得清晰连续。

图5给出了图4(a)、图4(b)的频谱图，图5(a)是原始实际地震资料 的频谱图；(b)是应用地震数据高分辨率有效弱信号能量增强方法处理后的地 震记录的频谱图，对比图5(a)和(b)可见应用地震数据高分辨率有效弱信 号能量增强方法处理后的地震记录其频带较原始实际地震资料有明显的拓宽。

图6给出了图4(a)的分频示意图，图6(a)是原始实际地震资料25Hz 的分频结果；图6(b)是原始实际地震资料35Hz的分频结果；图6(c)是原 始实际地震资料45Hz的分频结果。图7给出了图4(b)的分频示意图，图7 (a)是应用地震数据高分辨率有效弱信号能量增强方法处理后的地震记录的 25Hz的分频结果；图7(b)是应用地震数据高分辨率有效弱信号能量增强方 法处理后的地震记录的35Hz的分频结果；图7(c)是应用地震数据高分辨率 有效弱信号能量增强方法处理后的地震记录的45Hz的分频结果。对比图6(a) 和图7(a)可见，应用地震数据高分辨率有效弱信号能量增强方法处理后的地 震记录在25Hz分频效果图中，其3400ms左右的薄储层分辨率更高，连续性更 好。对比图6(b)和图7(b)可见，应用地震数据高分辨率有效弱信号能量 增强方法处理后的地震记录在35Hz分频效果图中，剖面的信噪比、储层的分 辨率有所提高。对比图6(c)和图7(c)可见，应用地震数据高分辨率有效弱 信号能量增强方法处理后的地震记录在45Hz分频效果图中，剖面的信噪比有 所提高，3400ms左右的储层连续性更好，分辨率更高。