本发明涉及石油地震勘探数据处理和解释领域，是通过一种串联型地震信号优化时频变换方法来解决其他常规时频变换方法存在的分辨率和能量聚集性较低的问题，改善时频分辨率和时频能量在时频平面中的聚集性，有效提高从时频域地震信号中提取地层结构的能力，可作为地震信号处理中可靠的时频变换技术。
背景技术：
时频变换方法是直接应用地震资料进行烃类监测的重要手段之一，也是地震信号处理的重要步骤。时频变换方法旨在通过构造一种时间和频率的密度函数，将一个一维的时间信号以二维的函数形式表示出来，以揭示地震信号中所包含的频率分量，使我们不但能够同时了解到地震信号在时间-频率域的信息，而且可以精准地刻画地震信号随时间变化的特征。因此，时频变换方法对提高地震资料处理和解释的精度以及可靠性有着非常重要的意义。
现有技术中地震信号时频变换方法很多，例如：1)分数阶变换方法，是在短时傅里叶变换等基础上提出的一种对非平稳信号分析的方法，将原来的变换方法从时间-频率域推广到时间-分数阶域，其具有短时傅里叶变换的优点，增强了处理非平稳信号的灵活性。2)同步挤压变换，是在其它变换(如连续小波变换)基础上根据信号各组成分量的时频特征、压缩时频曲线，达到提高时频分辨率的方法。它属于一种时频重排算法，与原始的连续小波变换不同的是，它在提高时频分辨率的同时，也支持信号的重构，主要用于检测地震信号中的低频微弱信号。具体实现的方法有同步挤压小波变换。3)自适应稀疏时频分析方法，它首先通过交替分裂Bregman迭代算法自适应计算给定地震信号的最优高斯窗，使其在时频平面的能量紧凑性达到最大值；其次，利用生成的最优窗在凸约束下进行逆分解，使其混合范数的转置系数最小化，通过以上两个步骤达到提高时频分辨率的目的。
随着地震勘探的不断发展和完善，人们对时频分析的精度要求也大大提高。常规的时频分析方法往往只利用了通频带以内的地震信息，对截频和限频以外的信息很少利用，尤其是地震信号的振幅信息变化较大时，对于那些高、低频弱振幅信息的识别和利用更少。然而，由于传播路径的复杂性及不同时刻不同岩性的地层对地震波吸收、衰减程度的差异使得地震信号呈现明显的时变、非平稳等特征的这类地震信号，需要在时间-频率域同时精确的揭示信号的时频分布特征，这对地震信号的时频分析方法的精度和可靠性提出了更高的要求。
技术实现要素：
本发明的目的是提供一种通过地震信号的自适应固有模态函数分解，并用优化时频变换求解每个固有模态函数的时频谱，通过叠加时频谱来提高时频分辨率的串联型地震信号优化时频变换方法。这种新方法不仅可以提高地震信号的时频分辨率，而且充分利用了原始地震道中的有效信息，同时具有很高的计算效率。
本发明所用的时频分析方法，摒弃了短时傅里叶变换、小波变换等时频变换方法的窗函数在整个时间信号上滑动，汲取了上述时频变换方法以及经验模态分解的优点并进行了改进和发展，可以精确地刻画地震信号中隐藏的高、低频弱振幅信息，其主要优点在于：构建了由信号分量主频自动控制最优参数的频域窗函数，频率覆盖范围广且时频分辨率高，同时具有很高的计算效率，适用于实际地震数据的处理。
本发明充分利用了地震勘探信号中的有效信息，可将其中的有效信息精确地展示在时频平面中，细致刻画了地震信号中不同频率信号分量随时间变化的特征，为后期地震资料的解释提供了有力的保障。
一种串联型地震信号优化时频变换方法，其具有以下优越性：
(1)本发明由于对输入的原始地震道进行完备集经验模态分解(CEEMD)，使输入的原始地震道中各信号分量的频率特征逐渐趋于窄带化，利用寻优算法快速找到对应的最优频域窗函数，以便充分利用原始地震道中的高、低频弱振幅信息。
(2)本发明使用的最优窗函数，消除了其他时频变换的窗函数参数存在的人为主观影响因素，并利用构建的能量最优准则和优化时频分布计算公式，使时频平面中的能量聚集性达到最优，使地震信号的时频分辨率得到了明显提高。
(3)常规的时频变换方法是窗函数沿着整个时间信号滑动，而本发明使最优窗函数与输入信号在频率域快速匹配，以便于最优窗可以在很窄的有限带宽内快速滑动，因此大大减少了时频变换所需的时间，使计算效率得到数倍的提升。
本发明的具体实现原理如下：
首先输入三维地震勘探数据，从其中抽出某一道地震记录进行完备集经验模态分解。假设输入的原始地震道为x(t)，根据经验模态分解原理找到其对应的第一个分量信号h1(t)，并根据经验模态函数(IMF)的判定条件检验h1(t)是否为IMF。
重复上述操作N次，得到N个IMF。对N个IMF求均值，作为完备集经验模态分解的第一个IMF。重复上述操作依次得到第i个IMF，即：
式中，ωk是均值为零且方差为1的高斯白噪声，ε是固定系数，i表示进行完备集经验模态分解的次数。
以每个IMF分量为新的输入信号，此处，将其表示为imti(t)，则输入信号imti(t)的主频fd如下计算：
式中，F[·]表示快速傅里叶变换，max[·]表示求取最大值的运算。
然后，依据主频fd来确定频域窗函数G(f,fd)的位置，因此新的输入信号imti(t)的主频与频域窗函数的主频同步移动，此时的频域窗函数G(f,fd)如下定义：
其中，η和γ为控制G(f,fd)的参数。
利用主频为fd的频域窗函数G(f,fd)和新输入信号imti(t)构建能量最优关系准则的目标函数，得
式中，IMFi(f)表示新输入信号imti(t)的傅里叶变换,即：
IMFi(f)＝F[imti(t)]
利用上述能量最优关系准则目标函数，求取最优窗参数γopt如下：
然后，构建最优频域窗函数如下：
利用上式的最优频域窗函数，建立时频分辨率明显提高的优化时频分布计算公式，则输入信号imti(t)的优化时频分布TFRi(t,f)的计算公式如下：
式中，b＝2fd，a＝fd/5，b-a表示新输入信号imti(t)的带宽。由于该公式在很窄的有限带宽内计算，因此大大减少了时频变换的计算时间，使计算效率得到数倍的提升。
最后，可计算原始地震道x(t)的优化时频分布TFR(t,f)，其表达式如下：
附图说明
图1是合成的Cross-chirp信号的时频谱对比图。其中，横坐标为时间，单位为秒，图1(a)纵坐标为振幅，图1(b)-1(d)纵坐标为频率，单位为赫兹。图1(b)是采用本发明的技术得到的时频谱，图1(c)是采用常规S变换得到的时频谱，图1(d)是采用小波变换得到的时频谱。
图2是某地区采集的一道地震数据的时频谱对比图。其中，横坐标为时间，单位为秒，图2(a)纵坐标为振幅，图2(b)-2(d)纵坐标为频率，单位为赫兹。图2(b)是采用本发明的技术得到的时频谱，图2(c)是采用常规S变换得到的时频谱，图2(d)是采用小波变换得到的时频谱。
图3是某地区采集的二维地震数据的瞬时谱剖面对比图，横坐标为道号，纵坐标为时间，单位为秒。其中，图3(b)是采用本发明的技术得到的25Hz的瞬时谱剖面图，图3(c)是采用常规S变换得到的25Hz瞬时谱剖面图。
具体实施方式
为了实现上述目标，本发明采取以下技术方案：一种串联型地震信号优化时频变换方法，其包括以下步骤：1)输入原始地震道x(t)；2)利用完备集经验模态分解将原始地震道分解成多个不同频段的固有模态函数(IMF)窄带信号分量，并将每个固有模态函数窄带信号分量作为新的输入信号imti(t)；3)求取每个新的输入信号imti(t)的主频fd，并依据能量最优准则目标函数求取每个频域窗函数的位置及其最优参数γopt；4)利用优化时频分布计算公式依次计算每个新输入信号imti(t)的时频分布谱TFRi(t,f)；5)将4)得到的每个时频分布谱进行叠加，得到原始地震道x(t)的时频分布谱，从而得到了时频分辨率明显提高的地震道时频分布谱TFR(t,f)。
本发明的实施实例说明：
图1是合成的Cross-chirp信号的时频谱对比图。在图1(b)中可以明显看出，信号的频率分辨率在低频端和高频端都很高，能量聚集性均很好，表明时频分辨率高。图1(c)是常规S变换的时频谱，虽然较好地反映了信号的时频分布特征，但存在一些缺陷：在低频端，信号的能量聚集性和频率分辨率较好；在高频端，能量聚集性较差，且随着频率的增高，频率分辨率和收敛性变得更差。图1(d)是小波变换的时频谱，但是存在明显的缺陷：整个时频平面的能量聚集性比图1(a)差，且低频端的时频谱出现了畸变。
图2是某地区采集的一道地震数据的时频谱对比图。在图2(c)-(d)的时频谱中，70Hz以上的高频信息消失，另外，图2(d)的其时间-频率关系不能很好地对应，增加了后期解释的不确定性，但在图2(b)中，可以精确描述70-150Hz隐藏的高频信息，较好地反映了原始地震道中的高频信息。因此，图2(b)弥补了图2(c)-(d)的缺点，既挖掘了其他时频变换方法不能刻画的隐藏的高频信息，又显示了更高的时频分辨率。
图3是某地区采集的二维地震数据的瞬时谱剖面对比图。图3(b)是采用本发明得到的25Hz瞬时谱剖面，可以清晰地刻画目的地层的能量分布情况，整个瞬时谱剖面的分辨率较强，在0.60s、0.70s及0.95s处能量连续性较好且比较集中，能更好地对应能量变化与地层的关系。在图3(c)的常规S变换的25Hz瞬时谱剖面中，其能量连续性较差，只能看到0.60s及0.95s处的能量较连续，整个瞬时谱剖面的分辨率和能量聚集性较差，不能清晰地刻画25Hz的信号能量随地层的分布规律。通过比较可以看出：采用本发明得到的瞬时谱剖面，其能量的横向连续性、聚焦性以及分辨率等都得到了较大的提升。
上述各实施例仅用于说明本发明，其中方法的各实施步骤等都是可变化的，凡在本发明技术方法的基础上进行的类似变换和改进，均不排除在本发明的保护范围之外。