一种基于级联式挤压时频变换的地震资料分析方法及系统

1.本发明属于勘探地球物理技术领域，具体涉及一种基于级联式挤压时频变换的地震资料分析方法及系统。


背景技术：

2.地震资料分析不仅可用于油气勘探领域中的储层预测，而且在描述薄互层的空间连续性等方面也具有广泛应用。时频分析作为一种强有力的工具在地震资料分析中尤为重要，其在实际地震资料分析中已取得广泛应用。未来时频分析技术应用于地震资料分析方面的研究仍具有很大的潜力，而如何获得高精度的地震记录时频谱，是后期进行地震资料分析与处理的首要前提。
3.目前已有的分析地震资料的时频分析方法主要包含：
4.现有技术1：通过平滑伪wvd分布和cohen类方法来分析地震资料。该技术直接对地震资料计算相应的时频谱，从而根据时频谱特征来对相应的地震资料属性进行分析。该技术虽然计算简单，但获得的时频谱在根本上很难消除交叉项干扰，导致时频分辨率降低。
5.现有技术2：以挤压时频分析为代表的地震资料分析技术。该类技术主要对原始的时频变换结果在频率方向上进行重排，能获得高分辨率的时频谱以分析地震资料。该方法对于频率相距较近、变化较快的多频率分量记录获得的时频分辨率有限，因此影响了后续的地震资料分析结果。


技术实现要素：

6.本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于级联式挤压时频变换的地震资料分析方法及系统，用于解决传统方法存在的交叉项干扰及时间分辨率和频率分辨率不可兼得的技术问题。
7.本发明采用以下技术方案：
8.一种基于级联式挤压时频变换的地震资料分析方法，包括以下步骤：
9.s1、根据单道地震记录的多个窗函数下短时傅里叶变换renyi熵最小选择最优大时窗；
10.s2、根据步骤s1选取的最优大时窗，进行最优大时窗下的二阶挤压时频变换；
11.s3、根据感兴趣目标区域确定时频范围，对步骤s2得到的最优大时窗下的二阶挤压时频变换进行重构，得到重构信号；
12.s4、对步骤s3得到的重构信号进行一系列较小时窗下的二阶挤压时频变换，并选取重构误差最小对应的窗函数为最优小时窗；
13.s5、根据步骤s4选取的最优小时窗进行二阶挤压时频变换，对地震数据进行时频分析。
14.具体的，步骤s1中，最优大时窗的宽度σ
lopt
：
15.16.其中，v(x,g
σ
；t,η)为单道地震数据对应的短时傅里叶变换，g
σ
为窗函数，τ为时间变量，η为频率变量。
17.进一步的，单道地震数据对应的短时傅里叶变换v(x,g
σ
；t,η)为：
[0018][0019]
其中，τ为临时变量，x(τ)为单道地震信号，i为虚数单位。
[0020]
具体的，步骤s2中，二阶高阶挤压时频变换为：
[0021][0022]
其中，为t和η位置处的二阶瞬时频率，δ为脉冲函数，ω为频率变量，为最优大时窗。
[0023]
进一步的，t和η位置处的二阶瞬时频率为：
[0024][0025]
其中，表示取复数的实部，表示窗函数为表示窗函数为表示对信号x(t)采用窗函数做短时傅里叶变换的结果，为t和η位置处的一阶瞬时频率，q(t,η)为频率调制因子。
[0026]
具体的，步骤s3中，重构信号xr(t)为：
[0027][0028]
其中，η0为重构频率的下限，η1为频率的上限，为最优大时窗函数在0时刻的取值，为二阶高阶挤压时频变换。
[0029]
具体的，步骤s4中，重构误差error(xr,σ)最小时对应的最优小时窗宽度σ
sopt
为：
[0030][0031]
进一步的，重构误差error(xr,σ)为：
[0032][0033]
其中，δω为一个极小的频率量，xr(t)为重构信号，g
σ
(0)为窗函数在零时刻的取值，为重构信号小时窗下的二阶挤压时频变换，
[0034]
具体的，步骤s5中，针对步骤s3中的重构信号xr(t)，采用步骤s4中所获得的最优小时窗运用步骤s2得到重构信号最优小时窗下的二阶挤压时频变换
从时频域精细分析地震资料。
[0035]
第二方面，本发明实施例提供了一种基于级联式挤压时频变换的地震资料分析系统，包括：
[0036]
选择模块，根据单道地震记录的多个窗函数下短时傅里叶变换renyi熵最小选择最优大时窗；
[0037]
变换模块，根据选择模块选取的最优大时窗，进行最优大时窗下的二阶挤压时频变换；
[0038]
重构模块，根据感兴趣目标区域确定时频范围，对变换模块得到的最优大时窗下的二阶挤压时频变换进行重构，得到重构信号；
[0039]
优化模块，对重构模块得到的重构信号进行一系列较小时窗下的二阶挤压时频变换，并选取重构误差最小对应的窗函数为最优小时窗；
[0040]
分析模块，根据优化模块选取的最优小时窗进行二阶挤压时频变换，对地震数据进行时频分析。
[0041]
与现有技术相比，本发明至少具有以下有益效果：
[0042]
一种基于级联式挤压时频变换的地震资料分析方法，首先对单道地震记录进行较大时窗下的二阶挤压时频变换，然后对感兴趣的目标区域进行重构，接着对上述重构结果进行较小时窗下的二阶挤压时频变换，从而可依据较小时窗下的二阶挤压分析结果来对地震资料进行分析。相比较于已有地震资料的时频分析方法，本发明可提供一种具有极高时频分辨率的时频谱，进而可用于目标层附近的薄层分析等地震处理解释环节，实现地震资料的高精度时频分析。
[0043]
进一步的，根据单道地震记录在不同大时窗下的短时傅里叶变换时频的renyi熵最小原则选取最优大时窗，有利于在众多不同宽度的窗函数中优选出最佳的大时窗函数。
[0044]
进一步的，对单道地震数据计算短时傅里叶变换，有利于获得地震信号时频谱轮廓，并为后续挤压时频变换提供基础的时频变换结果。
[0045]
进一步的，对单道地震记录采用最优大时窗下进行二阶挤压变换，有利于获得高时频分辨率时频谱及不同频段的时频特征。
[0046]
进一步的，计算t和η位置处的二阶瞬时频率，有利于后续挤压操作能够将时频变换系数挤压到正确的位置以获得高精度的时频变换结果。
[0047]
进一步的，对大时窗下的二阶挤压变换目标区域进行重构，有利于提取出特定的时频特征和减小特定频段的重构误差。
[0048]
进一步的，对重构结果在不同时窗下二阶挤压短时傅里叶变换重构误差最小原则选取最优小时窗，有利于针对重构结果在众多不同宽度的窗函数中优选出最佳的窗函数。
[0049]
进一步的，重构误差error(xr,σ)为每个小时窗下重构误差，为后续后续优选最优小时窗提供了一种代价函数。
[0050]
进一步的，对重构结果进行较小时窗下的二阶挤压时频变换，有利于获得时频分辨率更高的时频谱，根据较小时窗下的二阶挤压时频变换结果来获得地震属性信息，有利于对地震资料进行更为精细的时频分析。
[0051]
可以理解的是，上述第二方面的有益效果可以参见上述第一方面中的相关描述，在此不再赘述。
[0052]
综上所述，本发明可以有效、准确地实现提取地震资料的时频谱，采用级联式挤压时频变换，具有高分辨率和快速精确等特点，可用于地震资料的高精度时频分析。
[0053]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0054]
图1为本发明流程图；
[0055]
图2为实际二维地震资料记录及两口井的位置图示；
[0056]
图3为本方法得到的井1周边单道数据目标层附近的时频谱；
[0057]
图4为本方法得到的井2周边单道数据目标层附近的时频谱。
具体实施方式
[0058]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0059]
在本发明的描述中，需要理解的是，术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
[0060]
还应当理解，在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
[0061]
还应当进一步理解，在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合，例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
[0062]
应当理解，尽管在本发明实施例中可能采用术语第一、第二、第三等来描述预设范围等，但这些预设范围不应限于这些术语。这些术语仅用来将预设范围彼此区分开。例如，在不脱离本发明实施例范围的情况下，第一预设范围也可以被称为第二预设范围，类似地，第二预设范围也可以被称为第一预设范围。
[0063]
取决于语境，如在此所使用的词语“如果”可以被解释成为“在
……
时”或“当
……
时”或“响应于确定”或“响应于检测”。类似地，取决于语境，短语“如果确定”或“如果检测(陈述的条件或事件)”可以被解释成为“当确定时”或“响应于确定”或“当检测(陈述的条件或事件)时”或“响应于检测(陈述的条件或事件)”。
[0064]
在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的，其中为了清楚表达的目的，放大了某些细节，并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的，实际中可能由于制造公差或技术限制而有所偏差，并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。
[0065]
本发明提供了一种基于级联式挤压时频变换的地震资料分析方法，根据单道地震
记录的时频谱特征选取最优大时窗，对单道地震记录进行最优大时窗下的二阶挤压时频分析，并将目标区域分析结果进行重构，最终对上述重构结果选择最优小时窗下并进行最优小时窗下的二阶挤压时频分析，从而依据最优小时窗下的二阶挤压分析结果对地震资料进行分析，实现地震资料的高精度时频分析。
[0066]
请参阅图1，本发明一种基于级联式挤压时频变换的地震资料分析方法，以基于级联式挤压时频变换为基础，首先对地震资料在较大时窗下进行二阶挤压时频变换并针对感兴趣区域进行重构分量，其次对上述重构分量进行较小时窗下的二阶挤压时频变换来精细分析地震资料，用于目标层附近的薄层分析等地震处理解释环节，具体步骤如下：
[0067]
s1、计算地震道不同时窗下的短时傅里叶变换谱，选取renyi熵最小时对应的时窗作为最优大时窗；
[0068]
单道地震数据为x(t)，窗函数为σ为窗函数的窗宽。
[0069]
则x(t)对应的短时傅里叶变换为v(x,g
σ
；t,η)：
[0070][0071]
其中，τ为临时积分变量，η为频率变量。
[0072]
然后给定大时窗的变化范围为σ＝[σ
lmin
ꢀ…ꢀ
σ
lmax
]，σ
lmin
为大时窗变化范围的下限，σ
lmax
为大时窗变化范围的上限；计算不同时窗参数下的短时傅里叶变换结果v(x,g
σ
；t,η)，并计算每一个窗函数σ下的v(x,g
σ
；t,η)的renyi熵：
[0073]
re nyi[v(x,g
σ
；t,η)]＝∫∫|v(x,g
σ
；t,η)|2dtdη
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0074]
然后选择renyi熵最小时对应的时窗宽度为最优大时窗的宽度σ
lopt
：
[0075][0076]
s2、对选取的单道地震记录进行最优大时窗下的二阶挤压时频变换；
[0077]
信号为x(t)，选取最优窗函数为x(t)的短时傅里叶变换为则其对应的二阶高阶挤压时频变换
[0078][0079]
其中，为t和η位置处的二阶瞬时频率，通过如下方式获得：
[0080][0081]
其中，表示取复数的实部，表示窗函数为表示窗函数为表示对信号x(t)采用窗函数做短时傅里叶变换的结果；
[0082]
t和η位置处的一阶瞬时频率通过如下方式获得：
[0083]
[0084]
其中，表示在t方向的偏导数；
[0085]
频率调制因子q(t,η)通过如下方式得到：
[0086][0087]
s3、对二阶挤压时频分析获得的结果进行特定频段重构；
[0088]
信号x(t)对应的二阶挤压时频变换为则可针对特定频段重构信号xr(t)：
[0089][0090]
其中，η0为重构频率的下限，η1为频率的上限，为最优大时窗函数在0时刻的取值。
[0091]
s4、计算重构信号不同小时窗下的二阶挤压傅里叶变换，选取二阶挤压傅里叶变换重构误差最小时对应的时窗作为最优小时窗；
[0092]
给定小时窗的变化范围为σ＝[σ
smin
ꢀ…ꢀ
σ
smax
]，其中σ
smin
为时窗下限，σ
smax
为时窗上限；利用步骤s2中相关流程计算不同小时窗参数下的二阶挤压傅里叶变换然后寻找时频谱中的脊线位置ridge(xr,g
σ
；t)：
[0093][0094]
然后计算每个小时窗下二阶挤压傅里叶变换沿脊线ridge(xr,g
σ
；t)的重构误差error(xr,σ)：
[0095][0096]
其中，δω为一个极小的频率量(典型取值小于0.1hz)。
[0097]
最终选择重构误差error(xr,σ)最小时对应的时窗宽度为最优小时窗的宽度σ
sopt
：
[0098][0099]
s5、对重构结果进行最优小时窗下的二阶挤压时频变换，并进行地震资料的时频分析；
[0100]
针对步骤s3中的重构信号xr(t)，采用步骤s4中所获得的最优小时窗运用步骤s2得到重构信号最优小时窗下的二阶挤压时频变换依据此时频变换中的可判定相应的信号特征，从而可从时频域精细分析地震资料。
[0101]
本发明再一个实施例中，提供一种基于级联式挤压时频变换的地震资料分析系统，该系统能够用于实现上述基于级联式挤压时频变换的地震资料分析方法，具体的，该基于级联式挤压时频变换的地震资料分析系统包括选择模块、变换模块、重构模块、优化模块
以及分析模块。
[0102]
其中，选择模块，根据单道地震记录的多个窗函数下短时傅里叶变换renyi熵最小选择最优大时窗；
[0103]
变换模块，根据选择模块选取的最优大时窗，进行最优大时窗下的二阶挤压时频变换；
[0104]
重构模块，根据感兴趣目标区域确定时频范围，对变换模块得到的最优大时窗下的二阶挤压时频变换进行重构，得到重构信号；
[0105]
优化模块，对重构模块得到的重构信号进行一系列较小时窗下的二阶挤压时频变换，并选取重构误差最小对应的窗函数为最优小时窗；
[0106]
分析模块，根据优化模块选取的最优小时窗进行二阶挤压时频变换，对地震数据进行时频分析。
[0107]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0108]
以某油田采样间隔为1ms的包含1200道记录的地震资料为例。
[0109]
请参阅图2，图2为实际二维地震资料记录及测井1(图中用well1标识出位置)和测井2((图中用well2标识出位置))图示，采样间隔为1ms，共有1200道地震记录。
[0110]
请参阅图3，为测井1附近一道地震数据(第75道)的基于级联式挤压时频变换的时频分析结果。
[0111]
从图3可以看到多道记录的时频脊线在1300ms处具有倾斜向下的变化趋势，表明地层厚度随深度的增加而增大。同时，测井1在此位置的测井曲线也证实地层厚度沿时间方向逐渐增大。
[0112]
请参阅图4，为测井2附近一道地震数据(第760道)的基于级联式挤压时频变换的时频分析结果。
[0113]
从图4看到多道记录的时频脊线在1325ms处具有倾斜向上的变化趋势，表明地层厚度随深度的增加而减小。
[0114]
相应地，测井2在此位置的测井曲线也表明层厚逐渐减小。因此上述两个目标区域的时频分析验证了本方法在地震资料分析中的有效性。
[0115]
综上所述，本发明一种基于级联式挤压时频变换的地震资料分析方法及系统，采用大时窗挤压时频变换以获得高频率分辨率时频变换结果，并以此重构感兴趣时频区域的信号，采用小时窗挤压时频变换以获得高时间分辨率时频变换结果，最终相比于传统的地震资料时频分析方法，可以获得同时兼具高频率分辨率和高时间分辨率的时频谱，有利于对地震资料进行高精度时频分析。
[0116]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实
施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0117]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0118]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0119]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0120]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。