基于高分辨率时频分析和一致性度量的非均质性刻画方法与流程

本发明属于地球物理勘探中的信号处理领域，涉及一种基于高分辨率时频分析和一致性度量的非均质性刻画方法。

背景技术：

时频分析可以为地震资料处理和解释提供有用的信息，准确的时频表示对于刻画细小的地层结构，探测和油气储层相关的异常区域至关重要，尤其是在复杂隐蔽油气藏勘探中，对时频分析的准确性提出了更高的要求。传统的时频分析工具，如加窗傅里叶变换(stft)、小波变换(cwt)、s变换(st)等，由于受到窗函数或者小波函数的影响，在非平稳信号的时频表示方面具有一定的局限性，而地震信号恰恰就是典型的非平稳信号。

地震数据的相干体分析技术是地震资料解释的强有力工具，它能清晰地刻画地下断层及岩性变化(即非均质性)等，可提高地震解释的效率。该技术在石油勘探与开发中已得到广泛应用。相干体技术方法主要是利用地震波形的相似性度量。前人已提出了多种相干体技术分析方法。如：第一代相干体(c1)、第二代相干体(c2)、第三代相干体算法(c3)等，这些算法各有其优缺点，但从研究波形的相似性方面看，利用了地震道之间的线性相关性。目前的技术方法针对的对象都是非常规隐蔽性的油气藏，因此对于新的相干体算法，需要更精细刻画地震波形的变化，在有些情况下，常用方法的精度还不能满足要求。事实上，若把地震信号看作随机变量的观察值，两个随机变量之间的关系有多种，线性相关只是其中一种。

技术实现要素：

本发明的目的在于克服上述现有技术的缺点，提供一种基于高分辨率时频分析和一致性度量的非均质性刻画方法，通过计算地震数据的高分辨率的时频分析(tpw同步挤压变换)后，提取单频分量，进行一致性度量的相干体分析，从而进行地震数据的分均质性刻画。将该方法应用于实际地震资料，能够更加清晰地刻画地震数据的非均质性。

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

基于高分辨率时频分析和一致性度量的非均质性刻画方法，包括以下步骤：

步骤1：根据公式(1)计算三参数小波变换系数；

其中，s为原始输入信号，a为尺度因子，b为时移因子，s(t)为实地震道，为ψ(t)的复共轭，为伸缩和平移之后的子波，π为圆周率，s(ω)为原始输入信号的傅里叶变换，ω为角频率，为带有尺度因子的子波的傅里叶变换，e^iωb为平移因子的傅里叶域，i为虚数因子，是基本小波ψ(t)的fourier变换；

步骤2：根据下式计算重排准则；

步骤3：由下式计算重排后的时间-频率域系数；

步骤4：确定有效信号分布空间，提取单频地震数据体；

步骤5：根据下式计算单频数据体的一致性度量相干体，进行地震数据的非均质性刻画；

一致性度量值：

其中，sign表示符号函数，k表示inline方向道号，l表示croosline道号，n表示采样点数，m表示采样点数，n表示总采样点数；cm为分析窗内目标单频地震道分析点的一致性度量，j为一致性度量矩阵的维数，trace(·)表示矩阵的迹，对单频数据体中所有的分析点上述2式计算其一致性度量相干结果，即得到单频数据体的相干估计结果。

本发明进一步的改进在于：

步骤1中，计算三参数小波变换系数的具体方法如下：

式(1)为实地震道s(t)的小波变换定义，根据式(1)可知基本小波平方可积，无直流分量，满足以下容许性条件：

针对薄互层地震信号含有快速变化的振幅和频率分量的特点，选取三参数小波作为小波变换的母小波，三参数小波的定义为：

其中，τ为能量衰减因子，β为能量延迟时间，σ为分析小波调制频率，σ,τ,β∈r且σ,τ＞0；用向量λ＝(σ,τ,β)记参数集合，则ψ(t；σ,τ,β)记为ψ(t；λ)，其他量类似；以向量的形式，上式简写为：

其中，

对(4)式作fourier变换，得到其频域形式为：

步骤2中，计算重排准则的具体方法如下：

假设一个单频余弦信号，即s(t)＝acos(ω0t)，s(t)的fourier变换为根据(1)式，信号s(t)关于三参数小波ψ(t；λ)的小波变换结果为

由于三参数小波没有负频率分量，即ω＜0，且只考虑正尺度，则(7)式简化为

从(8)式看到，假如三参数小波ψ(t；λ)的峰值频率为ξm，则小波变换的结果将在尺度a＝ξm/ω0处取到最大值，并以这个能量最大的尺度为中心形成一个尺度带，造成能量的扩散，为了得到更集中的时频分布，计算重排准则，即瞬时频率：

其中，ws(a,b)≠0；当ws(a,b)＝0时，信号在此处无能量，(9)式将无穷大，所以计算时要求当噪声小于-10db时，选择阈值为10^-8；

如果噪声大于-10db，选取自适应阈值作为阈值median是matlab中的库函数，代表取中值，0.6745是针对高斯噪声的正则化因子；假设有n个时刻，所以对n个时刻求得的ση取均值，能得到最优的阈值；

当参考信号是单频余弦信号时，其小波变换会形成一个尺度带，但是由(9)计算其瞬时频率为：

ωs(a,b)＝ω0(10)。

步骤3中，时间-频率域系数的具体计算方法如下：

根据(9)式计算得到瞬时频率后，将小波变换系数累加到这个频率成分上，即(a,b)→(ωs(a,b),b)，进行能量重排；

对于给定的信号s(t)∈l²(r)的三参数小波变换为ws(a,b)，有如下表达式：

记而是s(t)的解析信号，从而(11)式写为：

因为sa(b)是s(b)的解析信号，故有

将ws(a,b)a^-1都分配给瞬时频率成分上，则就得到挤压后的时频分布，由此得出时间-尺度域到时间-频率域的映射如下：

通过(14)式，在某个固定的时刻b，计算其瞬时频率ωs(a,b)，将所有瞬时频率都为某一频率ω的小波系数通过(14)式累加到一起，这样就完成了能量的重分配，得到了挤压后的时频分布，(14)式的等价形式为：

步骤4中，确定有效信号能量分布空间的具体方法如下：

对重排后的时间频率域系数运用百分比阈值策略进行噪声压制，即软阈值操作，定义为：

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

δ＝p·max{|ts|}(17)

其中p为百分比，ts为本发明变换后的时间频率域系数；

通过滤波策略，得到有效信号对应的时频域系数，这些时频域系数构成噪声被压制以后的信号对应的时间-频率域谱图，根据上述时间-频率域谱图提取单频时频分量，即单频地震数据体，进行地震数据的分均质性刻画。

步骤5的具体方法如下：

基于步骤(4)中提取的单频地震数据体，选取在空间上的inline方向corssline方向分析窗分别记为wⁱ和w^c，在时间方向上分析窗记为w^t，得到一个子数据体；令ni×nc＝j，设wⁱ和w^c中分别有ni和nc道数据，w^t中包含nt个采样点；

取ni＝3，nc＝3，nt＝n；根据单频地震数据的位置标号，按顺序依次将地震到排列为数据矩阵d，数据矩阵d包含n行与j列；数据矩阵d的第i行、第j行元素记为dij，dij是第j个地震道的第i个采样点，表达式如(18)所示

数据矩阵d中的每一列是为分析窗内的一个单频地震道数据，将k列和l列作为两个随机变量，通过式(19)计算得到一致性度量值：

式中，sign表示符号函数；

根据式(19)，对于所有单频数据体中的地震道两两组合，计算其一致性度量的值，也就是说对所有的{(k,l)|k,l＝1,2,…,j}求一致性度量得到一致性度量矩阵c，c为一个正定矩阵，表达式如下：

采用renyi散度度量准则进行计算，其定义为：

式中，j为一致性度量矩阵的维数，trace(·)表示矩阵的迹，cm为分析窗内目标单频地震道分析点的一致性度量；对单频数据体中所有的分析点根据式(19)和(21)计算其一致性度量相干结果，即得到单频数据体的相干估计结果。

与现有技术相比，本发明具有以下有益效果：

本发明提出结合一种高分辨率的时频分析方法(基于三参数小波的同步挤压变换)和一致性度量相干算法的地震数据非均质性刻画方法。时频分析可以为地震资料处理和解释提供有用的信息，在复杂隐蔽油气藏勘探中，精确的时频表示对于指示细小的地层结构和检测含油气储层至关重要。相干体技术是三维地震资料解释中的重要技术，目前常用的相干体技术主要是用线性相关系数作为度量准则，而实际上有些地震道之间的相干性还包含有非线性的相干性。本发明结合一种高分辨率的时频分析技术(tpw同步挤压变换)和基于一致性度量的相干体技术，进行地震数据非均质性刻画。

附图说明

图1为40hz的ricker子波的时频域谱图；其中，(a)无噪的ricker子波，(b)、(c)无噪ricker子波的小波变换和本发明提出变换的时频图；(d)加入信噪比为5db的高斯白噪声后的ricker子波，(e)、(f)含噪ricker子波的小波变换和本发明提出变换的时频图。在无噪和含噪情况下，本发明提出的变换都能得到更稀疏、更稳定的时频结果，谱图能够清晰地反映出有效信号的能量分布；

图2为一致性度量相干体计算中多个地震道的选取示意图；

图3为实际地震数据相干体计算实例；其中，(a)原始实际地震数据剖面，(b)c3相干体结果，(c)一致性度量相干体结果；

图4为原始3d地震数据的沿层切片；该3d数据体crossline道数为1001，inline道数为701；

图5为40hz单频分量数据体的沿层切片；对原始3d地震数据体进行不同的时频变换后，提取40hz的单频数据体的沿层切片，结果分别是基于(a)morlet小波变换，(b)morlet同步挤压变换，(c)tpw小波变换，(d)tpw的同步挤压变换(本发明)。

图6为3d地震数据体的非均质性刻画结果；其中，(a)morlet小波变换和c3相干体，(b)morlet小波变换和一致性度量相干体，(c)morlet同步挤压变换和c3相干体，(d)morlet同步挤压变换和一致性度量相干体，(e)tpw小波变换和c3相干体，(f)tpw小波变换和一致性度量相干体，(g)tpw同步挤压变换和c3相干体，(h)tpw同步挤压变换和一致性度量相干体(本发明)。

具体实施方式

下面结合附图对本发明做进一步详细描述：

参见图1-6，本发明基于高分辨率时频分析和一致性度量的非均质性刻画方法，包括以下步骤：

(1)小波变换定义

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

其中，为ψ(t)的复共轭，a为尺度因子，b为时移因子，是基本小波ψ(t)的fourier变换，是s(t)的fourier变换。基本小波平方可积，无直流分量，满足以下容许性条件：

针对薄互层地震信号含有快速变化的振幅和频率分量的特点，以高静怀提出的最佳匹配地震子波(bmsw)的小波为基础，借鉴harrop等人的工作对传统的morlet小波加以改进，高静怀等人提出了一种小波函数——三参数小波(tpw)。由于具有三个参数，通过调节，它可以灵活的适应我们的需要，将它应用于层序检测和高精度地震资料分辨中都取得了很明显的效果。本发明选取三参数小波作为小波变换的母小波。三参数小波的定义为：

其中，τ为能量衰减因子，β为能量延迟时间，σ为分析小波调制频率，(σ,τ,β∈r且σ,τ＞0)。为了书写方便，用向量λ＝(σ,τ,β)记参数集合，则ψ(t；σ,τ,β)可记为ψ(t；λ)，其他量类似。以向量的形式，上式可简写为：

其中，

对(4)式作fourier变换，得到其频域形式为

(2)计算重排准则(瞬时频率)

假设一个单频余弦信号，即s(t)＝acos(ω0t)，s(t)的fourier变换为根据(1)式，信号s(t)关于三参数小波ψ(t；λ)的小波变换结果为

因为我们选择的三参数小波几乎没有负频率分量，即ω＜0，且只考虑正尺度，则(7)式可以简化为

从(8)式可以看到，假如三参数小波ψ(t；λ)的峰值频率为ξm，则小波变换的结果将在尺度a＝ξm/ω0处取到最大值，并以这个能量最大的尺度为中心形成一个尺度带，造成能量的扩散，为了得到更集中的时频分布，计算重排准则——瞬时频率：

其中，ws(a,b)≠0。当ws(a,b)＝0时，信号在此处无能量，(9)式将无穷大，计算瞬时频率ωs(a,b)没有意义，所以计算时要求当噪声较小时，可以选择阈值为10^-8；如果噪声较大，需要抑制噪声，我们选取自适应阈值作为阈值其中ση是和噪声水平相关的，它是由小波变换系数的前nv个尺度来估计的，median是matlab中的库函数，代表取中值，0.6745是针对高斯噪声的正则化因子。假设有n个时刻，所以可以对n个时刻求得的ση取均值，能得到最优的阈值。

当参考信号是单频余弦信号时，其小波变换会形成一个尺度带，但是可以由(9)计算其瞬时频率为

ωs(a,b)＝ω0(10)

从(10)式可以看到，(9)式定义的瞬时频率就是单频余弦信号的频率，即对于一个余弦信号，它的三参数小波变换得到的时间-尺度域结果会在某个能量最大的尺度附近形成一个尺度带，但是这些尺度对应小波系数通过(10)式计算出来的瞬时频率都为余弦信号的频率ω0，因此我们可以想象将这些尺度的能量都集中到ω0上，从而得到更集中的时频表示。

(3)能量重排——时间尺度域能量映射到时间频率域

根据(9)式计算得到瞬时频率后，我们知道能量应该往这个频率上集中，也就是将小波变换系数累加到这个频率成分上，即(a,b)→(ωs(a,b),b)，进行能量重排。

对于给定的信号s(t)∈l²(r)的三参数小波变换为ws(a,b)，有如下表达式：

记而是s(t)的解析信号，从而(11)式可以写为：

因为sa(b)是s(b)的解析信号，故有

从(12)中可以看到，通过小波变换系数再乘以因子a^-1得到的结果和原信号的解析信号只差一个常数因子，这很容易得到其反变换，如(13)式所示。如果我们将ws(a,b)a^-1都分配给我们上面提到相应的瞬时频率成分上，则就可以得到挤压后的时频分布，由此我们得出时间-尺度域到时间-频率域的映射如下：

通过(14)式，在某个固定的时刻b，计算其瞬时频率ωs(a,b)，将所有瞬时频率都为某一频率ω的小波系数通过(14)累加到一起，这样就完成了能量的重分配，得到了挤压后的时频分布，(14)的等价形式为：

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

当选用合适的小波函数(本发明选取三参数小波)对含噪信号进行小波变换，将其投影到时间尺度域的时候，有效信号的能量将会分布在较小的子空间v，噪声的能量会扩散到比较大的子空间v′，甚至整个时间-尺度域。换句话说，有效信号和少数的系数相对应，而噪声几乎分布在全部的系数上。如果我们确定了有效信号对应的子空间v，得到有效信号对应的系数，将其它系数置零，则噪声就会在变换域被压制，信噪比得到提高。本发明通过同步挤压变换之后，时间尺度域系数经过重排，变换到时间频率域，得到更为稀疏的时频结果，有效信号能量分布空间变小，更有利于噪声压制。

为了更有效地确定有效信号的能量分布空间，我们提出了对重排后的时间频率域系数运用百分比阈值策略进行噪声压制，即软阈值操作，定义为

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

δ＝p·max{|ts|}(17)

其中p为百分比，ts为本发明变换后的时间频率域系数。

通过该滤波策略，得到有效信号对应的时频域系数，这些系数构成噪声被压制以后的信号对应的时间-频率域谱图。图1(a)画出了不含噪的40hz的ricker子波的谱图，图1(b)和(c)为小波变换和本发明变换后的时频谱图，可见本发明变换能够得到更加集中、更加稀疏的时频表示。图1(d)为含噪的40hz的ricker子波(信噪比snr＝5db)的谱图，可见它已经被噪声污染，图1(e)和(f)为运用滤波策略后得到的小波变换和本发明变换的时频谱图，可以看出由于本发明变换更为稀疏，噪声得到明显压制，有效信号被清晰的表现出来。可以证明本发明能够得到分辨率更好的时频分析结果。基于该时频分析方法，可以提取单频时频分量，进行地震数据的分均质性刻画。

(5)基于一致性度量的相干体技术

利用步骤(4)中提取的高分辨率的单频(如40hz)时频分量地震数据体，结合基于一致性度量的相干体技术之后，可以更加精确地刻画地震数据的非均质性。

基于步骤(4)中提取的单频地震数据体，选取在空间上的inline方向corssline方向分析窗分别记为wⁱ和w^c，在时间方向上分析窗记为w^t，可以得到一个子数据体。令ni×nc＝j，设wⁱ和w^c中分别有ni和nc道数据，w^t中包含nt个采样点。

在图2中，选择窗(如矩形窗a所示)，取ni＝3，nc＝3，nt＝n。根据图2中单频地震数据的位置标号，按顺序依次将地震到排列为数据矩阵d(即n行与j列)，矩阵的第i行，第j行元素(记为dij)是第j个地震道的第i个采样点，表达式如(18)所示

当目标道移动到记号为5'或者5”的地震道时，包含目标道的分析窗如图2中虚线矩形框所示，会向inline方向上或者crossline方向上移动(如矩形窗b和窗c)。矩阵d中的每一列是为分析窗内的一个单频地震道数据，将k列和l列作为两个随机变量，通过式(19)可计算得到一致性度量值

式中，sign表示符号函数。

根据式(19)，对于所有单频数据体中的地震道两两组合，计算其一致性度量的值，也就是说对所有的{(k,l)|k,l＝1,2,…,j}求一致性度量得到一致性度量矩阵c，c为一个正定矩阵，表达式如下

对于该正定矩阵c，如果直接计算特征值运算量大。故本发明采用renyi散度度量准则进行计算，其定义为

式中，j为一致性度量矩阵的维数，trace(·)表示矩阵的迹，cm为分析窗内目标单频地震道分析点的一致性度量。对单频数据体中所有的分析点根据式(19)和(21)计算其一致性度量相干结果，即可得到单频数据体的相干估计结果。

实施例：

为了验证本发明方法的有效性，我们将其应用于实际地震资料的非均质性刻画。

图3是利用不同相干体算法进行断层刻画。从实际的资料剖面椭圆和矩形框中所标记的位置处有小的断层。通过相干算法计算得到的相干剖面，如图(b)和(c)，对应位置处，一致性度量相干体技术对于断层的刻画更清晰，与c3算法比较，对相干信息具有一定的增强作用，刻画地质特征信息更明显。

图4是三维地震数据体的一个沿层切片，可见该资料随机噪声十分严重，严重影响地震数据非均质性刻画。我们将不同时频分析方法应用于该地震数据，并提取40hz的单频数据体，进而提取沿层切片，对比不同时频变换方法对该地震数据的刻画能力，结果如图5所示。由图5可见，基于三参数小波(tpw)的同步挤压变换时频分辨率更高，可以更好的刻画地质结构(如红色箭头所示的河道，绿色箭头所示的断层)。

最后，我们对不同的时频分析方法得到的单频数据体(图5结果)，利用不同的相干体算法计算相干体，从而尝试进行地震数据的非均质性刻画，即将图5中单频数据体作为原始数据引入不同的相干体算法。结果如图6所示，图6中结果分别为(a)morlet小波变换和c3相干体，(b)morlet小波变换和一致性度量相干体，(c)morlet同步挤压变换和c3相干体，(d)morlet同步挤压变换和一致性度量相干体，(e)tpw小波变换和c3相干体，(f)tpw小波变换和一致性度量相干体，(g)tpw同步挤压变换和c3相干体，(h)tpw同步挤压变换和一致性度量相干体(本发明)。对比不同的时频分析方法和不同的相干体技术结合后的刻画结果可见，本发明(tpw同步挤压变换和一致性度量相干体方法)，对地震数据的非均质性(河道和断层等地质信息)刻画更为精确和明显。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。