CN104914464B - 基于相空间滤波策略的地震瞬时属性提取方法 - Google Patents

Abstract

本发明公开了一种基于相空间滤波策略的地震瞬时属性提取方法，在小波变换域自适应确定有效信号能量分布空间，在此基础上，利用小波变换与Hilbert变换的关系，提出了含噪信号瞬时属性分析的方法，本发明提出衡量瞬时属性抗噪性能的定量指标，通过计算不同信噪比的合成地震记录的瞬时属性，对比不同方法的抗噪性能，本文方法的计算结果具有良好的抗噪性能和精度，将该方法应用于实际资料，能够更加清晰地刻画储层的特征。

Description

基于相空间滤波策略的地震瞬时属性提取方法

技术领域

本发明属于地球物理勘探中的信号处理领域，涉及地震资料的瞬时属性提取，尤其涉及一种基于相空间滤波策略的地震瞬时属性提取方法。

背景技术

地震属性是指由地震数据经过数学变换导出的有关地震波的几何形态、运动学特征、动力学特征和统计学特征，这些特征能够从视觉上反映储层的形态及其含油气性。20世纪70年代，地震属性分析开始被引入地震解释中。20世纪90年代以来，由于储层描述和三维数据体解释的需要，地震属性分析技术迅速发展，地震属性在地层构造解释、储层岩性和物性特征描述、油气藏预测与动态监视方面得到了广泛应用。目前还没有一个公认的地震属性分类，按照物理意义，地震属性可以分为时间、振幅、频率、相干、衰减等几大类，按照属性拾取方法，每类又可以分为瞬时属性、单道时窗属性、多道时窗属性。其中瞬时属性是在地震波到达的位置上拾取的属性，包括瞬时振幅、瞬时相位、瞬时频率、瞬时带宽等等。近年来，地震属性的数量突增，然而，瞬时属性仍然是地震数据地质解释的支柱，也是目前商业处理软件必备的技术模块。本发明主要研究地震资料瞬时属性的提取。

在信号处理领域，以瞬时振幅、瞬时相位、瞬时频率为代表的信号瞬时属性的概念由来已久，在Gabor和Ville的工作之后，大量学者在相关领域做了研究，Boashash对这些工作做了综述。Taner等人于1979年引入复地震道分析，提出通过复地震道求取瞬时属性的方法，并给出了这些属性的物理意义以及在地震解释中的应用。瞬时振幅与相邻层的岩性变化及油气聚集有关。瞬时相位反映界面的不连续性、断层、不整合面和层序边界等。瞬时频率的变化可以有效刻画地层的厚度和岩性变化，指示油气的分布等。Robertson和Nogami将瞬时属性用于孔隙砂岩薄层厚度估计，Chopra和Marfurt利用瞬时属性进行不连续性、断层和横向不连续性检测。Liu和Marfurt利用瞬时属性检测和刻画曲流河的分布并确定其厚度。Zeng利用瞬时频率异常来指示薄层。Gao等人利用瞬时频率进行地震资料Q值估计。

常用的提取地震瞬时属性的方法可以分为以下三类：

(1)Hilbert变换方法。对地震道做Hilbert变换，转化为复地震道，其中实部为原地震道数据，虚部为其Hilbert变换。得到复地震道之后，就可以在每个采样点计算振幅、相位和频率等属性，即瞬时属性。在Taner的工作之后，Hilbert变换法被广泛用于计算地震瞬时属性，至今大部分商业软件仍采用该方法。很多学者对该方法进行了发展和改进，Barnes于1996年提出二维复地震道分析的概念。Luo等人于2003年提出广义Hilbert变换并给出其在地球物理方面的应用。Barnes于2007年提出加权瞬时频率的概念。Lu和Zhang将Hilbert变换方法进行推广，提出基于自适应滤波器计算瞬时频率的方法。

(2)基于时频分析的方法。时频分析方法利用信号的时频分布来求取瞬时属性。Boashash等人提出基于时频分析的自适应瞬时频率估计方法。Stankovic等人利用自适应窗时频分布计算瞬时频率。高静怀等人提出了在相空间计算地震瞬时属性的方法。Marfurt等人提出基于窄带谱分析的地震资料瞬时属性提取方法。Steeghs和Drijkoningen提出基于二次时频分布的地震层序分析以及属性提取方法。Huang等人提出基于经验模式分解(EMD)计算瞬时频率的方法。Han等人将完全总体经验模式分解(CEEMD)方法用于提取地震资料的瞬时频率。

(3)基于反演的方法。Fomel等人提出利用反演的方法来获得局部属性，相比于瞬时频率，局部频率物理意义更加明确，且应用效果明显。Liu等人提出基于反演的方法，计算地震资料的瞬时频率。

常用的基于Hilbert变换估计瞬时参数的方法对噪声很敏感，且由于滤波器的截断效应，使得计算出的瞬时属性精度低。

发明内容

本发明目的在于克服现有技术的不足，提供了一种具有良好的抗噪性能和精度的基于相空间滤波策略的地震瞬时属性提取方法。

为达到上述目的，本发明采用以下技术方案：

基于相空间滤波策略的地震瞬时属性提取方法，包括以下步骤：

步骤1：由目标区时间和频率范围确定m和n的区间；

小波原子由基本小波进行伸缩和平移得到，即

在进行数值计算时，需要对尺度因子a和平移因子t进行离散化，对尺度因子进行指数化离散记为其中a₀为离散化步长且a₀＞1，平移因子t的离散化步长为t₀，离散化后的小波族表示为

平移因子t所对应的n的范围由待分析数据的时间范围确定；根据广义Morse小波尺度和频率的关系，确定尺度因子a的范围，进而确定m的范围；

步骤2：由公式(14)求变换域系数；

C_m,n＝〈s,ψ_m,n〉, (13)

其中C_m,n为标架系数，s为离散化的信号构成的向量，{Ψ_m,n}为小波紧标架，〈s,ψ_m,n〉表示二者的内积；

步骤3：阈值操作，得到有效信号对应系数；

对系数运用阈值策略，软阈值和硬阈值操作分别定义为

和

其中λ为阈值参数；

步骤3：由公式(2)计算解析信号；

其中，s(t)是信号，h(t)是s(t)的Hilbert变换，c(t)是s(t)对应的解析信号；

步骤4：根据公式(3)-(5)计算瞬时属性；

利用c(t)来计算瞬时振幅、瞬时相位和瞬时频率：

其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。

进一步，步骤1中广义Morse小波尺度a和频率f的关系由下式确定：

进一步，步骤3中阈值参数λ由下式计算得到：

其中p为百分比。

本发明提出了基于相空间滤波策略的地震瞬时属性提取方法，在小波变换域自适应确定有效信号能量分布空间，在此基础上，利用小波变换与Hilbert变换的关系，提出了含噪信号瞬时属性分析的方法。本发明提出衡量瞬时属性抗噪性能的定量指标，通过计算不同信噪比的合成地震记录的瞬时属性，对比不同方法的抗噪性能，本文方法的计算结果具有良好的抗噪性能和精度。将该方法应用于实际资料，能够更加清晰地刻画储层的特征。

此外，在数值计算时，需要对连续小波变换中的平移因子和尺度因子进行离散化以得到小波族，本发明给出一种离散化方式，推导出基于该离散化方式小波标架的上下界估计及标架判定准则，证明相比常用的Morlet小波，离散化的广义Morse小波更易构成紧标架。

附图说明

图1 50Hz的Ricker子波的时间尺度谱图；

图2合成地震记录和采用不同方法计算出来的瞬时频率；

图3不同信噪比的合成地震记录；

图4采用不同方法计算的含噪合成记录的瞬时频率；

图5用不同的小波函数计算的10dB含噪合成记录的瞬时频率；

图6实际资料算例。

具体实施方式

以下通过具体实施例和附图对本发明方案做具体说明：

小波变换法计算瞬时属性

实地震道s(t)的小波变换定义为

当基本小波ψ(t)为解析小波时，高静怀等证明了以下定理：

如果ψ(t)为解析小波，其实部ψ_R(t)为偶函数且其傅里叶变换Ψ_R(ω)满足如下条件：则对任意的信号有下列等式成立

其中h(t)是s(t)的Hilbert变换，c(t)是s(t)对应的解析信号。该定理建立了小波变换和Hilbert变换之间的关系。因此，本发明可以利用c(t)来计算瞬时振幅、瞬时相位和瞬时频率：

其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。

小波函数的选取

广义Morse小波是一类双参数解析小波族，最早由Daubechies等人研究时频局域化算子时导出，是时频联合局域化问题的解。其频域表达式为

其中U(ω)为单位阶跃函数，β和γ为小波的参数，且β＞0，γ＞0。α_β,γ为归一化常量：

α_β,γ≡2(eγ/β)^β/γ. (27)

本发明选用广义Morse小波计算信号的瞬时属性，主要原因有：

(1)广义Morse小波是严格解析的，而常用的Morlet小波只是在调制频率足够大的时候近似解析，在时间局域化要求比较高的时候，Morlet小波会出现负频率泄露。更重要的是，上面定理中要求基本小波为解析小波，小波偏离解析性会造成该定理不再成立，影响所求瞬时属性的准确性。而广义Morse小波严格解析，完全满足定理要求；

(2)广义Morse小波具有更高的自由度，可以通过调整参数，表现出不同的形态和特点，以匹配不同的地震信号；

(3)离散后的广义Morse小波更容易构成紧标架，使得小波变换的运算量大大降低，计算小波变换的系数和从小波系数重构信号变得很方便，同样也有利于后面的滤波策略中对系数的处理。另外，利用冗余的标架对于噪声的压制也非常有利。

小波函数的离散化方法

小波原子由基本小波进行伸缩和平移得到，即

在进行数值计算时，需要对尺度因子a和平移因子t进行离散化，对尺度因子进行指数化离散，记为其中a₀为离散化步长且a₀＞1，平移因子t的离散化步长为t₀。这样，离散化后的小波族可以表示为

在本发明的离散化方法中，t₀和待分析信号的采样间隔相同，因此离散化的小波族满足平移不变性，可以看做由字典平移产生。对于具备平移不变性的字典，有以下定理：

如果存在B≥A＞0，使得对于有

则该平移不变字典构成一个标架，A和B分别为标架的上界和下界。该定理将字典构成标架的条件转化为φ_m(x)的傅里叶变换Φ_m(ω)应满足的条件。

对于小波函数，基于本发明的离散化方法，有则上述定理的条件等价为

本发明采用舒适度δ刻画标架紧的程度：

舒适度越小，说明标架的下界A和上界B越接近，标架越接近于紧标架。当δ＝0时，上下界相等，标架为紧标架。

记

然后计算在不同的尺度离散化步长a₀下，离散化后的Morlet小波和广义Morse小波(GMW)构成的标架的上下界和舒适度，如表1所示。本发明选取广义Morse小波的参数为β＝1,γ＝3，为了确保Morlet小波的近似解析性，参数取σ＝6。从表中可以看出，基于本发明的离散化方式，Morlet小波族和广义Morse小波族均可以构成标架，但是广义Morse小波的舒适度更小，更容易构成紧标架，特别是在采样步长a₀比较大的时候，广义Morse小波族已经可以被看做紧标架，这可以大大降低运算量。

表1 离散化后的小波族构成的标架的上下界和舒适度

对于可以构成紧标架的小波族，和公式(1)相对应的离散化后的小波变换和信号重构可以表示为：

C_m,n＝〈s,ψ_m,n〉, (34)

其中C_m,n为标架系数，s为离散化的信号构成的向量，{Ψ_m,n}为小波紧标架，〈s,ψ_m,n〉表示二者的内积，A是标架的界(上界与下界相等)。离散化的小波族构成紧标架，计算小波变换的系数和从小波系数重构信号变得很方便，同样也保证了后面求解优化问题时算法的收敛性。

有效信号能量分布空间的确定

当选用合适的小波函数对含噪信号进行小波变换，将其投影到时间尺度域的时候，有效信号的能量将会分布在较小的子空间V，噪声的能量会扩散到比较大的子空间V′，甚至整个时间-尺度域。换句话说，有效信号和少数的系数相对应，而噪声几乎分布在全部的系数上。如果确定了有效信号对应的子空间V，得到有效信号对应的系数，将其它系数置零，则噪声就会在变换域被压制，信噪比得到提高。

为了确定有效信号的能量分布空间，本发明提出了以下滤波策略：

(1)通过公式(14)，将含噪数据投影到小波标架上。首先确定公式中m和n的范围。平移因子t所对应的n的范围由待分析数据的时间范围确定。根据广义Morse小波尺度和频率的关系，确定尺度因子a的范围，进而确定m的范围；

(2)对系数运用百分比阈值策略。软阈值和硬阈值操作分别定义为

和

其中阈值参数λ由下式计算：

其中p为百分比。

通过该滤波策略，得到和有效信号对应的系数，这些系数构成噪声被压制以后的信号对应的时间-尺度谱图。图1(a)画出了不含噪的50Hz的Ricker子波的谱图，图1(b)为含噪的50Hz的Ricker子波(信噪比SNR＝5dB)的谱图，可见它已经被噪声污染，图1(c)为运用滤波策略后得到的谱图，可以看出噪声被压制，有效信号被清晰的表现出来。

确定有效信号对应的能量分布空间后，就可以利用公式(2)求取相应的解析信号，然后计算各种瞬时属性。

图1为50Hz的Ricker子波的时间尺度谱图，在含噪信号的小波变换谱图中，有效信号能量分布受到噪声的影响，经过滤波策略以后，噪声得到压制，谱图清晰地反映出有效信号的能量分布。

本发明的物质基础是地震数据体，采用的逐道处理办法。具体步骤为：

步骤1：由目标区时间和频率范围确定m和n的区间；

步骤2：由公式(14)求变换域系数；

步骤3：阈值操作，得到有效信号对应系数；

步骤3：由公式(2)计算解析信号；

步骤4：根据公式(3)-(5)计算瞬时属性。

效果分析

合成信号算例

本发明对50Hz的Ricker子波同反射系数序列卷积生成的合成地震记录进行研究。图2(a)为不含噪的合成记录，图2(b)为加上高斯白噪声的含噪记录(信噪比SNR＝5dB)，图2(c)和(e)为分别采用Hilbert变换和本发明提出的基于滤波策略的方法计算的瞬时频率，其中广义Morse小波的参数选为β＝1,γ＝3。需要指出的是，对于不含噪信号，在滤波策略中没有进行阈值操作。从图2可以看出两种方法都能得到准确的结果。然而，当信号含有噪声时，图2(d)和(f)给出了分别采用两种方法计算的瞬时频率。在采用本发明提出的方法时，m根据合成记录的频率范围1Hz－150Hz确定，n根据时间范围确定。在滤波策略中，选取百分比p＝10％。从图2(d)可以看出，通过传统的Hilbert变换方法得到的瞬时频率结果中，有效信号的瞬时频率因为噪声的存在被完全掩盖。然而，通过本发明方法计算的结果(图2(f))仍然能够清晰地显示有效信号的瞬时频率。

图2合成地震记录和采用不同方法计算出来的瞬时频率，可以看出Hilbert变换法对噪声非常敏感，计算结果已经无法反映有效信号的瞬时频率，而本发明方法具有良好的抗噪性能和精度，在噪声较强时仍能够准确地求得有效信号的瞬时频率

接下来，本发明分别给图2(a)的合成记录加上不同程度的噪声，其信噪比分别为10dB、5dB和1dB，如图3所示。

图3不同信噪比的合成地震记录

使用本发明提出的方法计算不同信噪比合成记录的瞬时频率。在采用滤波策略时，分别使用硬阈值和软阈值。两种阈值方法的计算结果和Hilbert变换方法计算的结果如图4所示。

另外，为了定量衡量不同方法计算出来的瞬时频率的抗噪性能，本发明定义了瞬时频率信噪比：

其中f₀是不含噪信号的瞬时频率，f_n为含噪信号的瞬时频率。

各种方法的瞬时频率信噪比FSNR见表2。和传统的Hilbert变换方法相比，本发明的方法显著提高了瞬时频率的抗噪性能，即使数据被严重的噪声干扰，仍能得到有效信号对应的瞬时频率。另外，在阈值策略中使用软阈值的效果要优于硬阈值。

图4采用不同方法计算的含噪合成记录的瞬时频率。从上到下，信噪比分别为10dB、5dB和1dB，从左到右，蓝色线表示依次采用传统的Hilbert方法，本发明提出的方法(软阈值)和本发明提出的方法(硬阈值)的计算结果，红色的线为不含噪合成记录的瞬时频率。可以看出和传统的Hilbert变换方法相比，本发明的方法具有良好的抗噪性能和精度，同时，在滤波策略中采取软阈值的结果要好于硬阈值

表2 不同方法计算的不同信噪比的合成记录的瞬时频率信噪比FSNR

本发明继续检验选取不同的小波函数给结果带来的影响。在采用本发明提出的方法时，分别将Morlet小波和广义Morse小波作为基本小波，并运用同样的阈值策略，计算信噪比SNR为10dB的合成记录的瞬时频率，结果分别画在图5(a)和(b)，对应的FSNR分别为23.4937和13.8837，可见采用广义Morse小波计算的瞬时频率抗噪性能更强，结果更准确。

图5用不同的小波函数计算的10dB含噪合成记录的瞬时频率。蓝线为采用不同小波函数的计算结果，红线为不含噪合成记录的瞬时频率。可以看出采用广义Morse小波的效果要好于Morlet小波

实际地震资料算例

为了进一步检验本发明方法的有效性，本发明将其用于实际地震资料。图6(a)是三维数据体的一个过井剖面，图6(b)和(c)分别是采用Hilbert变换和本发明方法计算的瞬时频率。

图6实际资料算例，该数据为一个连过井剖面，井A是一个高产井，分别用Hilbert法和本发明方法计算的瞬时频率剖面，箭头指示的位置为含气层造成的频率异常区域

图6(b)的结果完全被噪声影响，而图6(c)的结果有效地反映了地质结构和储层分布。在1.75s、2s和2.25s处，可以看出由于薄层干涉带来的高频特征。在1.90s到1.95s之间存在着一个低频异常，图中用箭头指示。这一类低频阴影是由含气储层中高频能量的迅速衰减造成的，可以用来指示储层的分布。井A在该异常区域所在的层钻得优质的工业气流，和瞬时频率剖面上表现出来的异常吻合。由于本发明方法良好的抗噪性能图4-6(c)的结果更清晰地刻画了这一含气储层的频率异常。

Claims (3)

1.基于相空间滤波策略的地震瞬时属性提取方法，其特征在于包括以下步骤：

步骤1：由目标区时间和频率范围确定m和n的区间；

小波原子由基本小波进行伸缩和平移得到，即

${\mathrm{\ψ}}_{t,a}\left(\mathrm{\τ}\right)=\frac{1}{a}\mathrm{\ψ}\left(\frac{\mathrm{\τ}-t}{a}\right).---\left(1\right)$

在进行数值计算时，需要对尺度因子a和平移因子t进行离散化，对尺度因子进行指数化离散记为其中a₀为离散化步长且a₀＞1，m为整数，即m＝…,-3,-2,-1,0,1,2,3,…,平移因子t的离散化步长为t₀，离散化后的小波族表示为

${\mathrm{\ψ}}_{m,n}\left(x\right)=a_0^{-m}\mathrm{\ψ}\[a_0^{-m}(x-{\mathrm{nt}}_0)\]---\left(2\right)$

平移因子t所对应的n的为整数，即n＝…,-3,-2,-1,0,1,2,3,…，范围由待分析数据的时间范围确定；根据广义Morse小波尺度和频率的关系确定尺度因子a的范围，进而确定m的范围，其中f为频率(单位为Hz)，β和γ为广义Morse小波的参数，均为大于0的实数；

步骤2：由公式(3)求变换域系数；

C_m,n＝<s,ψ_m,n>, (3)

其中C_m,n为标架系数，s为离散化的信号构成的向量，{Ψ_m,n}为步骤1中离散化后的小波族构成的小波紧标架，<s,ψ_m,n>表示二者的内积；

步骤3：阈值操作，得到有效信号对应系数；

对系数运用阈值策略，软阈值和硬阈值操作分别定义为

$S_{\mathrm{\&lambda;}}\left(u\right)=\left\{\begin{array}{cc}u-\mathrm{\&lambda;}\frac{u}{\left|u\right|}& \left|u\right|\&GreaterEqual;\mathrm{\&lambda;},\\ 0& \left|u\right|<\mathrm{\&lambda;},\end{array}\right.---\left(4\right)$

和

其中λ为阈值参数，u为采取阈值操作的变量，S_λ(u)表示通过对变量u采取门限为λ的软阈值操作之后的结果，H_λ(u)表示通过对变量u采取门限为λ的硬阈值操作之后的结果；

步骤4：由公式(6)计算解析信号；

$c\left(t\right)\&equiv;\frac{1}{K_{\mathrm{\&psi;}}}{\&Integral;}_0^{\mathrm{\&infin;}}{W}_{\mathrm{\&psi;}}(t,a)\frac{da}{a}=s\left(t\right)+ih\left(t\right),---\left(6\right)$

其中，s(t)是信号，h(t)是s(t)的Hilbert变换，c(t)是s(t)对应的解析信号，为信号s(t)的小波变换，ψ(t)为基本小波，*代表复共轭，其中Ψ_R(ω)为基本小波ψ(t)实部的Fourier变换；

步骤5：根据公式(7)-(9)计算瞬时属性；

利用c(t)来计算瞬时振幅、瞬时相位和瞬时频率：

$\hat{e}\left(t\right)=\sqrt{{\mathrm{Re}}^2\&lsqb;c\left(t\right)\&rsqb;+{\mathrm{Im}}^2\&lsqb;c\left(t\right)\&rsqb;},---\left(7\right)$

$\widehat{\mathrm{\&theta;}}\left(t\right)=arctan\&lsqb;\frac{\mathrm{Im}\&lsqb;c\left(t\right)\&rsqb;}{\mathrm{Re}\&lsqb;c\left(t\right)\&rsqb;}\&rsqb;,---\left(8\right)$

$\hat{f}\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d}{dt}\&lsqb;arctan\&lsqb;\frac{\mathrm{Im}\&lsqb;c\left(t\right)\&rsqb;}{\mathrm{Re}\&lsqb;c\left(t\right)\&rsqb;}\&rsqb;\&rsqb;,---\left(9\right)$

其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。

2.根据权利要求1所述的基于相空间滤波策略的地震瞬时属性提取方法，其特征在于：步骤1中广义Morse小波尺度a和频率f的关系由下式确定：

$a=\frac{{(\mathrm{\&beta;}/\mathrm{\&gamma;})}^{1/\mathrm{\&gamma;}}}{2\mathrm{\&pi;}f}$

其中，β和γ为广义Morse小波的参数，均为大于0的实数。

3.根据权利要求1所述的基于相空间滤波策略的地震瞬时属性提取方法，其特征在于：步骤3中阈值参数λ由下式计算得到：

$\mathrm{\&lambda;}=p\&CenterDot;\underset{m,n}{max}\left\{\right|C_{m,n}\left|\right\},---\left(10\right)$

其中p为百分比。