CN104216017A - 空间相关的非平稳地震信号拓频方法 - Google Patents

Abstract

本发明公开了一种空间相关的非平稳地震信号拓频方法。其包括以下步骤：地震数据输入、空变处理、时变处理、求取频率补偿系数、平滑补偿系数和S变换振幅拓频处理。本发明的有益效果是：本发明的空间相关的非平稳地震信号拓频方法通过将S变换、时变处理、空变处理相结合，求得每一道地震记录的补偿系数，补偿其高频成分，拓宽频带宽度，实现地震记录信号的高保真性和高分辨率。

Description

空间相关的非平稳地震信号拓频方法

技术领域

本发明属于非平稳信号拓频技术领域，尤其涉及一种空间相关的非平稳地震信号拓频方法。

背景技术

自20世纪50年代起，石油、天然气等矿产资源消耗总量逐渐增大，石油供应不仅关系到民计民生，也关系到国家经济和军事保障。但随着地震勘探的不断深入，寻找地下油气藏的难度越来越大。近年来，复杂隐蔽油气藏的勘探成为主要地质勘探目标，这就要求地震资料具有高分辨率。因此，提高地震信号分辨率是保证地震解释准确的关键，在岩性划分和地震层位分析方面上都很有意义。根据信号随时间的变化规律，可将自然界中的信号简单地分为平稳信号与非平稳信号。地震信号作为一种短时、突变的非平稳随机信号，利用傅里叶变换处理的话，由于傅里叶变换自身的局限性，将导致处理结果精度很低。时频分析方法将一维时域平面和频域平面的信号二维联合表示，能够突出非平稳时间信号的局部性能，从而更加准确细致地描述信号特征。目前常用的时频分析方法有：短时傅里叶变换(Short-Time FourierTransform，STFT)、Gabor变换、小波变换(Wavelet Transform，WT)等。时频分辨率是衡量时频分析方法优劣和时频聚集性能的重要指标，这些传统的时频分析方法相比之下，都各有优劣：(1)STFT的优点是算法可直接利用实现容易，物理意义明确，反映整体时频趋势；缺点是其窗函数固定不变，导致时频分辨率也固定不变，低频端和高频端的时频分辨率相同，不能同时考虑地震资料中高频端和低频端的需求；(2)WT继承了STFT的优势，实现了变窗处理的要求，具有良好的时频局域化特性，用来处理地震信号有很好的效果；但WT的母小波有多种且没有选择的准则，其有限长的小波基容易造成能量泄露，并且尺度因子与频率的关系也较模糊，只是一种时间-尺度域变换方法。Stockwell(1996)等人提出的S变换吸收并发展了STFT和WT，它采用宽度随频率变化的高斯窗函数，其时频分辨率也随着频率变化，且基本小波不必满足容许性条件，具有类似多分辨率的特性。它的运算可直接利用FFT实现，大大降低了运算复杂性。S变换(ST)是近几年发展起来的一种时频分析算法，是短时傅里叶变换(STFT)和连续小波(WT)变换思想的结合，可以看作以高斯窗为窗函数，且随频率变化的短时傅里叶变换。其高斯窗的具体表现为在高频处较窄，在低频处较宽。与短时傅里叶变换和连续小波相比，S变换在时间域和频率域都具有更好的分辨率，且与傅里叶变换有直接关系，具有无损可逆性，在多个应用领域中发挥积极作用。地震信号是典型的非平稳信号，其频谱随时间变化而变化。如果对一道地震记录仅采用一个时窗进行频谱分析，求取补偿系数，则结果很可能存在浅层修正过量，或深层补偿不足。浅层修正过量会导致高频噪声放大，降低了信号的保真性；深层补偿不足，地震记录分辨率仍达不到目标要求。目前，地震拓频处理技术在不断发展完善当中，其中常见的方法主要有反Q滤波、时变谱白化、反褶积等。反Q滤波是基于吸收衰减模型设计反滤波因子，以达到补偿大地吸收滤波的目的。品质因子Q的准确程度将直接影响反Q滤波的整体效果，然而，对品质因子Q的准确求取往往很难做到，限制了该方法的广泛应用。时变谱白化实质是通过对分频信号的增益，实现原信号中不同频带数据的均衡，来达到提高数据纵向分辨率的目的。时变谱白化仅对信号的振幅谱进行处理，将振幅谱展平以实现频带的拓宽。但时变谱白化的基本假设前提是反射系数序列的谱是白的，即在全频带内其谱值都为1，但真实的地震记录却不满足该假设。反褶积是通过压缩子波，去除子波影响，以提高地震资料的时间分辨率的过程。但是该方法都必须基于一定的假设条件，如假定地震子波已知，而实际情况并不满足，这样便限制了反褶积方法的应用。且由于地下复杂的地质结构，常常导致反射地震记录的褶积模型不可靠。因此，虽然反褶积对于提高分辨率具有很好的效果，但也有许多局限性和缺陷。

发明内容

为了解决以上问题，本发明提出了一种空间相关的非平稳地震信号拓频方法。

本发明的技术方案是：一种空间相关的非平稳地震信号拓频方法，包括以下步骤：

S1.将地震记录数据进行输入；

S2.对每一道地震记录数据进行空变处理，以当前道为中心，对其两边取相邻N道，则将空间窗窗宽表示为：

spacewin＝2N+1，

将空间窗内的当前道表示为：

wintrace(N+1)，

其中，(N+1)表示窗内道序号；

S3.对每一道地震记录数据进行相同的时变处理，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，具体包括以下步骤：

S31.对每一道地震记录数据加入滑动时窗，时窗窗长为ΔT，求取时窗内地震记录数据的振幅谱；

S32.将步骤S31中得到的时窗内地震记录数据的振幅谱作为相应时窗中心点T_j处的时间点谱smx(f,T_j)，表示为：

smx(f,T_j)＝abs(FFT(trace(T_j:T_j+ΔT,i)))，

其中，abs表示取绝对值，FFT表示快速傅氏变换，trace表示当前道，f表示频率方向，i＝1,2,...inline，inline表示道号；

S33.根据步骤S32中得到的时间点谱，在频率方向f上利用最小二乘法插值，得到与地震记录数据变换振幅谱同维度的矩阵数据；

S34.将步骤S33中得到的与地震记录数据变换振幅谱同维度的矩阵数据进行叠加，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，表示为：

$\mathrm{Sumsmx}=\underset{i=1}{\overset{\mathrm{spacewin}}{\mathrm{\Σ}}}\mathrm{wintrace}\left(i\right).\mathrm{smx},$

其中，i表示道序号；

S4.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，求取频率补偿系数；

S5.对步骤S4中得到的频率补偿系数，进行曲面平滑处理；

S6.根据步骤S5中处理后的频率补偿系数，对地震记录数据进行S变换振幅拓频处理，具体包括以下步骤：

S61.对当前道作S变换，得到变换矩阵S(f,T)，

S62.将步骤S61中得到的S(f,T)与步骤S5中处理后的频率补偿系数相乘，得到拓频后的时频谱S′(f,T)；

S63.对步骤S62中得到的时频谱S′(f,T)进行逆S变换，得到拓频后的地震信号，完成拓频处理。

进一步地，上述步骤S4根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，求取频率补偿系数，具体包括以下步骤：

S41.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，沿时间方向T进行遍历；

S42.对每一道时间点谱自峰值频率起，对其衰减规律依据指数函数exp(-βf)进行指数拟合，得到相应的吸收参数-β；

S43.将β作为频率补偿系数，得到频率补偿系数矩阵gain，表示为：

gain＝exp(βf)。

进一步地，上述步骤S5对步骤S4中得到的频率补偿系数，进行曲面平滑处理，具体包括以下步骤：

S51.以步骤S4中得到的频率补偿系数矩阵内的各点为中心点进行遍历；

S52.取3×3移动窗口进行加权计算；

S53.根据步骤S52中得到的计算结果取权重均值代替原中心点。

本发明的有益效果是：本发明的空间相关的非平稳地震信号拓频方法通过将S变换、时变处理、空变处理相结合，求得每一道地震记录的补偿系数，补偿其高频成分，拓宽频带宽度，实现地震记录信号的高保真性和高分辨率。

附图说明

图1是本发明的空间相关的非平稳地震信号拓频方法的流程示意图。

图2是本发明的空变处理示意图。

图3是本发明的曲面平滑处理示意图。

图4是运用本发明的空间相关的非平稳地震信号拓频方法前后地震记录对比图。

图5是本发明实施例中原始剖面图。

图6是本发明实施例中运用本发明的空间相关的非平稳地震信号拓频方法后的剖面图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，为本发明的空间相关的非平稳地震信号拓频方法的流程示意图。一种空间相关的非平稳地震信号拓频方法，包括以下步骤：

S1.将地震记录数据进行输入。

S2.对每一道地震记录数据进行空变处理，以当前道为中心，对其两边取相邻N道，则将空间窗窗宽表示为：

spacewin＝2N+1，

将空间窗内的当前道表示为：

wintrace(N+1)，

其中，(N+1)表示窗内道序号。

如图2所示，为本发明的空变处理的示意图。我们在时间方向上采用滑动窗口的同时，在空间方向上也进行加窗处理，在横向上多道统计，以保证道间相对稳定性。即对每一道地震记录处理时，以当前道为中心，对其两边取相邻N道。空间窗包含2N+1道地震记录，将空间窗窗宽表示为spacewin＝2N+1，将当前道表示为wintrace(N+1)。

S3.对每一道地震记录数据进行相同的时变处理，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，具体包括以下步骤：

S31.对每一道地震记录数据加入滑动时窗，时窗窗长为ΔT，求取时窗内地震记录数据的振幅谱；

S32.将步骤S31中得到的时窗内地震记录数据的振幅谱作为相应时窗中心点T_j处的时间点谱smx(f,T_j)，表示为：

smx(f,T_j)＝abs(FFT(trace(T_j:T_j+ΔT,i)))，

其中，abs表示取绝对值，FFT表示快速傅氏变换，trace表示当前道，f表示频率方向，i＝1,2,...inline，inline表示道号；

S33.根据步骤S32中得到的时间点谱，在频率方向f上利用最小二乘法插值，得到与地震记录数据变换振幅谱同维度的矩阵数据；

S34.将步骤S33中得到的与地震记录数据变换振幅谱同维度的矩阵数据进行叠加，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，表示为：

$\mathrm{Sumsmx}=\underset{i=1}{\overset{\mathrm{spacewin}}{\mathrm{\Σ}}}\mathrm{wintrace}\left(i\right).\mathrm{smx},$

其中，i表示道序号。

通过对每一道地震记录数据加入滑动时窗，可以解决地震数据浅层修正过量，深层补偿不足的问题，提高信号的保真性和分辨率。

S4.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，求取频率补偿系数，具体包括以下步骤：

S41.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，沿时间方向T进行遍历；

S42.对每一道时间点谱自峰值频率起，对其衰减规律依据指数函数exp(-βf)进行指数拟合，得到相应的吸收参数-β；

S43.将β作为频率补偿系数，得到频率补偿系数矩阵gain，表示为：

gain＝exp(βf)。

S5.对步骤S4中得到的频率补偿系数，进行曲面平滑处理，具体包括以下步骤：

S51.以步骤S4中得到的频率补偿系数矩阵内的各点为中心点进行遍历；

S52.取3×3移动窗口进行加权计算；

S53.根据步骤S52中得到的计算结果取权重均值代替原中心点。

如图3所示，为本发明的曲面平滑处理的示意图。中心点(Point of Interest)的权重高于窗内其他点，若超出矩阵边界，则等价于以0值填充。取中心点权重为2，窗内其余点权重为1，则计算结果表示为：

$\begin{array}{c}\mathrm{gain}(i,j)=(\underset{\mathrm{jj}=j-1}{\overset{j+1}{\mathrm{\Σ}}}\mathrm{gain}(i-1,\mathrm{jj})+\underset{\mathrm{jj}=j-1}{\overset{j+1}{\mathrm{\Σ}}}\mathrm{gain}(i+1,\mathrm{jj})\\ +\mathrm{gain}(i,j-1)+\mathrm{gain}(i,j+1)\\ +2*\mathrm{gain}(i,j))\÷\underset{i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{j-1}{\overset{j+1}{\mathrm{\Σ}}}{\mathrm{weight}}_{\mathrm{ij}}\end{array},$

其中，i,j表示频率补偿系数矩阵gain的行列号。

S6.根据步骤S5中处理后的频率补偿系数，对地震记录数据进行S变换振幅拓频处理，具体包括以下步骤：

S61.对当前道作S变换，得到变换矩阵S(f,T)，

S62.将步骤S61中得到的S(f,T)与步骤S5中处理后的频率补偿系数相乘，得到拓频后的时频谱S′(f,T)；

S63.对步骤S62中得到的时频谱S′(f,T)进行逆S变换，得到拓频后的地震信号，完成拓频处理。

这里将当前道表示为trace(t,i)，t表示时域上的时间方向。当i为定值时，s(t)＝trace(t,i)表示一道地震记录数据。通过S变换，将地震记录数据变换到时频域上，得到变换矩阵S(f,T)。S变换是一种常用的时频分析算法，这里我们不作详细叙述。通过将变换矩阵S(f,T)与频率补偿系数矩阵gain相乘，得到拓频后的时频谱S′(f,T)，这里仅对振幅谱进行修正，保留原有的相位信息。对时频谱S′(f,T)进行逆S变换得到拓频后的地震信号，完成拓频处理。如图4所示，为运用本发明的空间相关的非平稳地震信号拓频方法前后地震记录对比图。如图5所示，为本发明实施例中原始剖面图。如图6所示，为本发明实施例中运用本发明的空间相关的非平稳地震信号拓频方法后的剖面图。由图4可以看出，运用本发明的空间相关的非平稳地震信号拓频方法可以使地震记录信号的高频成分得到有效补偿；由图5和6可以看出，地震记录信号分辨率得到显著提高。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种空间相关的非平稳地震信号拓频方法，其特征在于，包括以下步骤：

S1.将地震记录数据进行输入；

S2.对每一道地震记录数据进行空变处理，以当前道为中心，对其两边取相邻N道，则将空间窗窗宽表示为：

spacewin＝2N+1，

将空间窗内的当前道表示为：

wintrace(N+1)，

其中，(N+1)表示窗内道序号；

S3.对每一道地震记录数据进行相同的时变处理，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，具体包括以下步骤：

S31.对每一道地震记录数据加入滑动时窗，时窗窗长为ΔT，求取时窗内地震记录数据的振幅谱；

S32.将步骤S31中得到的时窗内地震记录数据的振幅谱作为相应时窗中心点T_j处的时间点谱smx(f,T_j)，表示为：

smx(f,T_j)＝abs(FFT(trace(T_j:T_j+ΔT,i)))，

其中，abs表示取绝对值，FFT表示快速傅氏变换，trace表示当前道，f表示频率方向，i＝1,2,...inline，inline表示道号；

S33.根据步骤S32中得到的时间点谱，在频率方向f上利用最小二乘法插值，得到与地震记录数据变换振幅谱同维度的矩阵数据；

S34.将步骤S33中得到的与地震记录数据变换振幅谱同维度的矩阵数据进行叠加，得到与地震记录数据变换矩阵维度相同的叠加振幅谱，表示为：

$\mathrm{Sumsmx}=\underset{i=1}{\overset{\mathrm{spacewin}}{\mathrm{\Σ}}}\mathrm{wintrace}\left(i\right).\mathrm{smx},$

其中，i表示道序号；

S4.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，求取频率补偿系数；

S5.对步骤S4中得到的频率补偿系数，进行曲面平滑处理；

S6.根据步骤S5中处理后的频率补偿系数，对地震记录数据进行S变换振幅拓频处理，具体包括以下步骤：

S61.对当前道作S变换，得到变换矩阵S(f,T)，

S62.将步骤S61中得到的S(f,T)与步骤S5中处理后的频率补偿系数相乘，得到拓频后的时频谱S′(f,T)；

S63.对步骤S62中得到的时频谱S′(f,T)进行逆S变换，得到拓频后的地震信号，完成拓频处理。

2.如权利要求1所述的空间相关的非平稳地震信号拓频方法，其特征在于：所述步骤S4根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，求取频率补偿系数，具体包括以下步骤：

S41.根据步骤S3中得到的与地震记录数据变换矩阵维度相同的叠加振幅谱，沿时间方向T进行遍历；

S42.对每一道时间点谱自峰值频率起，对其衰减规律依据指数函数exp(-βf)进行指数拟合，得到相应的吸收参数-β；

S43.将β作为频率补偿系数，得到频率补偿系数矩阵gain，表示为：

gain＝exp(βf)。

3.如权利要求1所述的空间相关的非平稳地震信号拓频方法，其特征在于：所述步骤S5对步骤S4中得到的频率补偿系数，进行曲面平滑处理，具体包括以下步骤：

S51.以步骤S4中得到的频率补偿系数矩阵内的各点为中心点进行遍历；

S52.取3×3移动窗口进行加权计算；

S53.根据步骤S52中得到的计算结果取权重均值代替原中心点。