地震数据高信噪比高分辨率时频谱的生成方法与流程

本发明涉及地球物理信号分析领域，特别是涉及一种地震数据高信噪比高分辨率时频谱的生成方法。

背景技术：

地震数据的时频谱分析是油藏地球物理研究的一个关键手段。直观地表述，时频谱分析是对任意地震单道按照其时间采样顺序逐点生成相应的频率谱，是将随着时间变化的一维地震道转换成随着时间和频率变化的二维“谱”。这个二维谱称之为时频谱。

维威（wigner-ville）时频分布方法是计算地震道时频谱的有效方法之一。其实现过程可分为三个基本步骤：计算地震道的解析信号，计算解析信号的瞬时自相关函数，计算自相关函数的傅立叶变换。因为维威时频分布的计算方法中采用了解析信号同时还采用了瞬时自相关函数，所以维-威时频分布具有高分辨率的特点。但是，正是因为瞬时自相关函数的计算却又导致了信号中所含各分量之间严重的交叉干扰。

为了抑制上述信号分量之间的交叉干扰，研究人员提出一种叫做平-虚维威分布的计算方法。它包含两层意义，一是“平滑”指的是沿时间方向对干扰的抑制，一是“虚拟”其实指的是沿频率方向对干扰的抑制。可是，无论是时间方向还是频率方向的抑制处理本质上是低通滤波，其效果是损失了高频成分，降低了时频谱的分辨率。

由于上述时频谱计算方法是单道分析方法，所生成的时频谱受地震数据的信噪比影响，不但时频谱的信噪比较低，而且道与道之间的时频谱缺乏空间一致性，对地震时频谱的解释应用造成了困难。

综合上述分析，标准维威时频谱分辨率高但有严重的交叉干扰，而平-虚维威时频谱虽抑制了交叉干扰但分辨率降低，加之这些时频谱的计算都是单道运算，所生成的地震时频谱信噪比差，受地震数据噪声的影响严重，对时频谱的应用和地球物理解释造成了困难。

技术实现要素：

本发明目的在于克服现有方法的上述缺陷，提出一种地震数据高信噪比高分辨率时频谱的生成方法。该地震时频谱的生成方法兼具高信噪比和高分辨率的两高优点。

为了实现上述目的，本发明提出的地震数据高信噪比高分辨率时频谱生成方法通过沿地层倾角的多道加权合成，利用地震反射同相轴的空间连续性，提高地震时频谱的信噪比；通过构建自适应蒙版滤波器的方法，增强地震时频谱的分辨率。该方法兼具高分辨率和高信噪比的两高优点。

作为优化，地震数据高信噪比高分辨率时频谱生成方法具体步骤是：（1）估算地震剖面上各时间-空间位置的倾角；（2）针对每一时间位置构建多道加权的合成道，并计算其解析信号的瞬时自相关函数；（3）计算标准维威（wigner-ville）时频分布及其相应的平-虚维威分布；（4）构建时频谱自适应的蒙版滤波器；（5）实现对标准维威时频分布的蒙版滤波，获得优化的地震道时频谱。

作为优化，步骤（1）中从实际地震剖面上估算各时间-空间位置的倾角是采用离散平面波析构方程的方式构建非线性方程组，通过求解非线性问题的计算方法估算时变和空变的地层倾角。

步骤（1）估算地震剖面上各时间-空间位置的倾角,采用平面波析构方程：

,

式中t是地震道的时间样点，x是地震道的空间位置，u(t,x)是地震波场，σ是线性同相轴的倾角。将上述平面波析构方程离散化并以矢量矩阵形式表示为

f(σ)=e^t[c(σ)○u]e=0,

式中u是数据矩阵，c(σ)是由二维滤波器的系数按照逆序组成的矩阵，[c(σ)○u]是按元素相乘或两个矩阵的哈达玛（hadamard）乘积，e是全1矢量。上式中引入全1矢量e实现对矩阵[c(σ)○u]中所有元素的求和。通过求解非线性方程组f(σ)≈0估算时变和空变的地层倾角。

作为优化，步骤（2）在每一时间位置实现以下三步：依据局部倾角进行多道加权构建合成道x(t)；计算合成道的解析信号z(t)；计算解析信号的瞬时自相关函数：k(t,τ)=z(t+τ/2)z*(t-τ/2),式中τ是瞬时自相关的延滞时间，z*(t)是解析信号z(t)的复数共轭。

作为优化，步骤（3）利用上述步骤（2）中得到的瞬时自相关函数，计算标准维威时频分布如下:

,

并计算相应的平-虚维威分布如下:

式中h(τ)是以时间为变量的对称函数，其作用是沿频率轴的低通滤波，g(t)是另一个以时间为变量的函数，其作用是沿时间轴的低通滤波,且g(0)=h(0)。

作为优化，步骤（4）蒙版滤波器构建式：

,

式中比值|wsp(t,f)|²/|w(t,f)|²是两个功率谱的相似度，ε是控制蒙版滤波器边界位置的调节参数，η则是控制蒙版滤波器边界平滑度的参数，设置η≥2。为了估算的稳定性起见，在（t,f）点处的功率谱估算是以该点为中心的二维小窗口内功率谱之和。

作为优化，构建的蒙版滤波器具有最平通带特征，其调节参数ε有效控制着滤波器的边界位置；设置边界平滑度参数η=4调节参数ε=0.05时，能够从时频谱中滤除任何功率谱比值小于20%的成像谱点。

作为优化，步骤（5）实现过程分为两步：对维威时频分布的振幅谱|w(t,f)|进行滤波，；重构维威时频谱，，式中ɵ(t,f)是维威时频分布的原始相位谱，最终获得高信噪比高分辨率的地震道优化时频谱。

本发明采用上述实施方案之后，具有如下有益效果：（1）本发明的有益效果在于其高分辨率特征，该方法生成的时频谱具有与标准维威分布一样的高分辨率，同时该方法具有与平-虚维威分布一样的清除交叉干扰的能力。（2）本发明的有益效果还在于其高信噪比特征，在运用实际三维地震数据进行储层预测时，该方法能够清楚地分开砂岩储层的空间分布和含煤层的空间分布。

附图说明

图1是本发明地震数据高信噪比和高分辨率时频谱生成方法的基本实现流程图；图2是本发明实施例的实际地震剖面（上图）及其估算的倾角（下图）图。

图3是本发明提出的蒙版滤波器，具有最平通带特征。图中横轴是两个功率谱的比值|wsp(t,f)|²/|w(t,f)|²，计算所用边界平滑度参数η值取为4，三条蒙版滤波器系数曲线对应的调节参数ε值分别为0.1,0.05,0.01。

图4是本发明应用于一个合成道的实施例，合成道含有4个子波，图中对比合成道的标准维威分布，合成道的平-虚维威分布，以及本发明方法生成的时频谱。图4显示本发明方法生成的优化时频谱具有与标准维威分布一样的高分辨率，同时该方法具有与平-虚维威分布一样的清除交叉干扰的能力。

图5是本发明应用于实际三维地震记录进行储层预测的实施例。上图是三维地震数据体的时间切片。中图的频谱切片反映该区煤系地层的空间分布，该区煤系地层的反射地震主频为15赫兹。下图的频谱切片则是该区砂岩储层的空间分布，该区砂岩储层反射地震主频为20赫兹。该图显示本发明方法的高信噪比特征，能够将砂岩储层和煤系地层的空间分布清楚地分开。

具体实施方式

本发明地震数据高信噪比高分辨率时频谱的生成方法是通过沿地层倾角的多道加权合成，利用地震反射同相轴的空间连续性，提高地震时频谱的信噪比；通过构建自适应蒙版滤波器的方法，增强地震时频谱的分辨率。具体步骤是：（1）估算地震剖面上各时间-空间位置的倾角；（2）针对每一时间位置构建多道加权的合成道，并计算其解析信号的瞬时自相关函数；（3）计算标准维威（wigner-ville）时频分布及其相应的平-虚维威分布；（4）构建时频谱自适应的蒙版滤波器；（5）实现对标准维威时频分布的蒙版滤波，获得优化的地震道时频谱。

步骤（1）中从实际地震剖面上估算各时间-空间位置的倾角是采用对平面波析构方程离散化的方式构建非线性方程组，通过求解非线性问题的计算方法估算时变和空变的地层倾角。步骤（1）估算地震剖面上各时间-空间位置的倾角,采用平面波析构方程：

式中t是地震道的时间点，x是地震道的空间位置，u(t,x)是地震波场，σ是线性同相轴的倾角。将上述平面波析构方程离散化并以矢量矩阵形式表示为f(σ)=e^t[c(σ)○u]e=0,式中u是数据矩阵，c(σ)是由二维滤波器的系数按照逆序组成的矩阵，其维数与数据矩阵u相同，[c(σ)○u]是按元素相乘或两个矩阵的哈达玛（hadamard）乘积，e是全“1”矢量。上式中引入全“1”矢量e实现对矩阵[c(σ)○u]中所有元素的求和。通过求解非线性方程组f(σ)≈0估算时变和空变的地层倾角。

步骤（2）在每一时间位置实现以下三步：依据局部倾角进行多道加权构建合成道x(t)；计算合成道的解析信号z(t)；计算解析信号的瞬时自相关函数：k(t,τ)=z(t+τ/2)z*(t-τ/2),式中τ是自相关函数的延滞时间，z*(t)是解析信号z(t)的复数共轭。

步骤（3）利用步骤（2）中自相关函数，计算标准维威时频分布如下:

,

并计算相应的平-虚维威分布如下:

,

式中h(τ)是以时间为变量的对称函数，其作用是沿频率轴的低通滤波，g(t)是另一个以时间为变量的函数，其作用是沿时间轴的低通滤波,且g(0)=h(0)。

步骤（4）蒙版滤波器构建式如下：

,

式中比值|wsp(t,f)|²/|w(t,f)|²是两个功率谱的相似度，ε是控制蒙版滤波器边界位置的调节参数，η则是控制蒙版滤波器边界平滑度的参数，设置η≥2。为了稳定起见，在(t,f)点处的振幅估算是以该点为中心的二维小窗口内振幅之和。构建的蒙版滤波器具有最平通带特征，其调节参数ε有效控制着滤波器的边界位置；设置边界平滑度参数η值为4调节参数ε值取为0.05时，能够从时频谱中滤除任何功率谱比值小于20%的成像谱点。

步骤（5）实现过程分为两步：对维威时频分布的振幅谱|w(t,f)|进行滤波，；重构维威时频谱，

，式中ɵ(t,f)是维威时频分布的原始相位谱，最终获得高信噪比高分辨率的地震道优化时频谱。

为使本发明的技术方案和优点更加清楚明白，下面结合附图对本发明实施作详细具体说明。本发明提供了一种地震数据高信噪比高分辨率时频谱的生成方法。图1是本发明方法的基本实现流程图。

步骤s101，估算地震剖面上各时间-空间位置的倾角。采用平面波析构方程：

式中t是地震道的时间点，x是地震道的空间位置，u(t,x)是地震波场，σ是线性同相轴的倾角。将上述平面波析构方程离散化并以矢量矩阵形式表示为

f(σ)=e^t[c(σ)○u]e=0

式中u是数据矩阵，c(σ)是由二维滤波器的系数按照逆序组成的矩阵，[c(σ)○u]是按元素相乘或两个矩阵的哈达玛（hadamard）乘积，e是全1矢量，上式中引入全1矢量e实现对矩阵[c(σ)○u]中所有元素的求和。通过求解非线性方程组f(σ)≈0估算时变和空变的地层倾角。

图2是本发明实施例的实际地震剖面（上图）及其估算的倾角（下图）。图中箭头是估算的局部倾角方向，与背景的实际地震实际的同相轴吻合。

步骤s102，在每一时间位置实现以下三步：依据局部倾角进行多道加权构建合成道x(t)；计算合成道的解析信号z(t)；计算解析信号的瞬时自相关函数：k(t,τ)=z(t+τ/2)z*(t-τ/2),式中τ是自相关函数的延滞时间，z*(t)是解析信号z(t)的复数共轭。

步骤s103，利用步骤s102中得到的瞬时自相关函数，计算标准维威时频分布如下:

,

并计算相应的平-虚维威分布如下:

式中h(τ)是以时间为变量的对称函数，其作用是沿频率轴的低通滤波，g(t)是另一个以时间为变量的函数，其作用是沿时间轴的低通滤波,且g(0)=h(0)。

步骤s104，蒙版滤波器构建式如下：

,

式中比值|wsp(t,f)|²/|w(t,f)|²是两个功率谱的相似度，ε是控制蒙版滤波器边界位置的调节参数，η则是控制蒙版滤波器边界平滑度的参数，设置η≥2。为了稳定起见，在(t,f)点处的振幅估算是以该点为中心的二维小窗口内振幅之和。如图3所示，构建的蒙版滤波器具有最平通带特征，其调节参数ε有效控制着滤波器的边界；设置边界平滑度参数η值为4调节参数ε值取为0.05时，能够从时频谱中滤除任何功率谱比值小于20%的成像谱点。

步骤s105，实现对标准维威时频分布的蒙版滤波。实现过程分为两步：对维威时频分布的振幅谱|w(t,f)|进行滤波，；重构维威时频谱，，式中ɵ(t,f)是维威时频分布的原始相位谱，最终获得高信噪比高分辨率的地震道优化时频谱。

本发明还提供两个实施例[图4、图5]，说明本发明具有如下有益效果：

（1）如图4所示，一个合成道含有4个子波，而在其标准维威分布上，除了对应4个子波的能量团之外，子波之间的交叉干扰非常明显；在合成道的平-虚维威分布上，虽然这些交叉干扰得以完全清除，但是4个有效能量团却明显扩散，因而时频谱的分辨率变低；最终，本发明方法生成的时频谱具有与标准维威分布一样的高分辨率，同时该方法具有与平-虚维威分布一样的清除交叉干扰的能力。

（2）图5是本发明应用于实际三维地震记录进行储层预测的实施例。上图是三维地震数据体的时间切片。中图的频谱切片所反映的是该区煤系地层的空间分布情况，该区煤系地层的反射地震主频为15赫兹。注意中图的煤系地层空间变化与上图的有些部位有相似之处。下图的频谱切片则是该区砂岩储层的空间分布情况，该区砂岩储层的反射地震主频为20赫兹。如图5所示，本发明的有益效果还在于其高信噪比特征，在运用实际三维地震数据进行储层预测时，该方法能够清楚地分开砂岩储层的空间分布和煤系地层的空间分布。