CN109001800B - 一种基于地震数据的时频分解与气藏检测方法及系统 - Google Patents
一种基于地震数据的时频分解与气藏检测方法及系统 Download PDFInfo
- Publication number
- CN109001800B CN109001800B CN201810801169.5A CN201810801169A CN109001800B CN 109001800 B CN109001800 B CN 109001800B CN 201810801169 A CN201810801169 A CN 201810801169A CN 109001800 B CN109001800 B CN 109001800B
- Authority
- CN
- China
- Prior art keywords
- time
- seismic data
- data
- frequency spectrum
- generating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 50
- 238000001514 detection method Methods 0.000 title claims abstract description 48
- 238000001228 spectrum Methods 0.000 claims abstract description 204
- 238000000034 method Methods 0.000 claims abstract description 25
- 230000009466 transformation Effects 0.000 claims abstract description 21
- 238000010521 absorption reaction Methods 0.000 claims abstract description 8
- 230000006870 function Effects 0.000 claims description 22
- 238000004458 analytical method Methods 0.000 claims description 20
- 238000004590 computer program Methods 0.000 claims description 11
- 238000003860 storage Methods 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 6
- 230000009286 beneficial effect Effects 0.000 abstract description 6
- 238000012545 processing Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 10
- 230000003595 spectral effect Effects 0.000 description 7
- 230000004044 response Effects 0.000 description 6
- 230000002776 aggregation Effects 0.000 description 5
- 238000004220 aggregation Methods 0.000 description 5
- 230000036039 immunity Effects 0.000 description 5
- 238000012937 correction Methods 0.000 description 4
- 238000009826 distribution Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 239000004215 Carbon black (E152) Substances 0.000 description 2
- 229920003776 Reny® Polymers 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 229930195733 hydrocarbon Natural products 0.000 description 2
- 150000002430 hydrocarbons Chemical class 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000007493 shaping process Methods 0.000 description 1
- 239000010902 straw Substances 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种基于地震数据的时频分解与气藏检测方法及系统，包括：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；根据所述最优时频谱，生成衰减梯度数据及分频数据；根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测。本申请具有提高分频数据与能量吸收衰减特征的属性剖面数据的精度的有益效果。
Description
技术领域
本发明涉及石油地球物理勘探技术领域，尤其涉及一种基于地震数据的时频分解与气藏检测方法及系统。
背景技术
利用地震数据进行储层预测的技术领域主要包括属性和反演两大部分，反演主要是叠后波阻抗与反射系数反演、叠前AVO(Amplitude Varies with Offset)属性与弹性阻抗反演等，主要反映地层弹性性质与储层流体特征等；属性则是地震数据一系列数学变换的总和，包括时间、振幅、频率、相位、波形等等，反映储层厚度、几何形态、空间展布和烃类赋存特征等。其中，频率域和时频域属性常常依赖于Fourier变换和诸多时频分析方法的理论支持。
目前地震数据频谱分解常用的方法包括短时傅里叶变换(Short Time FourierTranform，STFT)，小波变换(Wavelet Transform,WT)及时频连续小波变换(Time-Frequency Continous Wavelet Transform，TFCWT))，S变换类(包括：ST(S Transform)及GST(Generalized S Transform))，Wigner-Ville分布(Wigner-Ville Distributioin)类(包括：WVD，PWVD，SPWVD及RSPWVD)，匹配追踪分解(MPD)，经验模态分解系列(包括：EMD，EEMD及CEEMD)，反演谱与基追踪分解(BP)，同步挤压变换(包括：WSST及FSST)等等，尽管新方法层出不穷，目前工业界主流频谱分解技术仍以CWT和GST为主，这是其保守性、滞后性和基于成本考量的结果。若对精度有要求但计算成本可放宽限制，则以MPD居多。一般而言，地震解释中涉及三维数据体的计算，MPD频谱分解的计算效率常常让业界难以接受。事实上，如果CWT的参数选择恰当，仍然可以获得高于预期的高分辨率时频谱，从而为后续的频谱分解与衰减特征分析提供理想的属性剖面。
WT广义上是时频分析的一种方法，也是信号处理中一种常用分析工具，因其具备多尺度分辨特性的优势，常常用于时间序列的局部非平稳特性分析。一般而言，WT根据其定义公式中对时间变量、平移因子和尺度因子三个参数离散方式的不同，可以分为连续小波变换(Continuous Wavelet Transform,CWT)和离散小波变换两大类(Discrete WaveletTransform,DWT)。DWT基于框架理论和紧凑支撑的正交小波基，常用于计算机编码、语音处理与图像重构领域；CWT则常采用非正交连续小波函数作为核函数或者母小波，常用于多学科信号处理中时间序列非平稳与局部频率变化特征分析。在地震勘探领域，CWT得到了极其广泛的应用，包括地震数据的去噪、频散衰减分析与烃类检测等。
CWT在实际应用时，必须要面对的一个重要问题就是母小波的选取，不同的母小波对CWT的应用效果具有重大的影响，母小波选择不当甚至会对分析结果造成难以解释的困境。地震勘探的数据应用中，由于物理意义明确，且表达式简单易懂，Ricker子波与Morlet子波是目前最为常用的两种选择。Ricker子波在信号处理领域，也被称作墨西哥草帽小波，本质上是Gauss小波关于时间的负二阶导数；Morlet子波本质上是加入频率调制项的Gauss小波。Ricker子波主要用于井震标定和波场模拟，Morlet子波则主要应用于地震信号稀疏分解与衰减分析等方面。然而，在实际模拟与分析中发现，将一维时间信号映射到二维时频平面上时，基于Ricker子波和Morlet子波的CWT时频表征并不能获得聚集性最高的时频谱，一方面两类小波并不能拟合任意类型的波形信号，另一方面，对小波参数的控制缺乏针对待分析信号的优化标准。
为减少主观调参的影响，信息熵理论为时频分析的能量聚集性优劣提供了定量判据。
目前工业界主流频谱分解技术仍以连续小波变换(CWT)和广义S变换(GST)为主，这是其保守性、滞后性和基于成本考量的结果。若对精度有要求但计算成本可放宽限制，则以匹配追踪(MP)算法居多。一般而言，地震解释中涉及三维数据体的计算，MP算法频谱分解的计算效率常常让业界难以接受。而CWT和GST方法的参数选择是一个主观性较强的问题，其分辨率和抗噪性也有进一步提升的空间，因此，如何减少主观调参的影响，从而提高地震数据频谱分解与衰减属性分析分辨率及抗噪性，是当前亟待解决的技术问题。
发明内容
为了解决现有技术中的缺陷，本发明提供了一种基于地震数据的时频分解与气藏检测方法及系统，利用雷尼(Rényi)熵判据获取最优化的STFT第一时频谱，利用不同阶数的连续小波变换生成第二时频谱，通过将第一时频谱与第二时频谱进行互相关，获取最优化的Paul-CWT最优时频谱，并根据最优时频谱得到分频数据，同时利用EAA方法获取衰减梯度数据，进行高精度的气藏检测，具有提高分频数据与能量吸收衰减特征的属性剖面数据的精度、客观性和鲁棒性的有益效果。
为了实现上述目的，本发明提供了一种基于地震数据的时频分解与气藏检测方法，该方法包括：
根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；
根据所述最优时频谱，生成衰减梯度数据及分频数据；
根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测。
本发明还提供了一种基于地震数据的时频分解与气藏检测系统，该系统包括：
第一生成单元，用于根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
第二生成单元，用于根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
互相关单元，用于将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；
第三生成单元，用于根据所述最优时频谱，生成衰减梯度数据及分频数据；
气藏检测单元，用于根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测。
本发明还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；
根据所述最优时频谱，生成衰减梯度数据及分频数据；
根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测。
本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；
根据所述最优时频谱，生成衰减梯度数据及分频数据；
根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测。
本发明提供的一种基于地震数据的时频分解与气藏检测方法及系统，包括：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；根据所述最优时频谱，生成衰减梯度数据及分频数据；根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测，本申请具有提高地震数据频谱分解与衰减属性分析分辨率及抗噪性的有益效果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本申请的一种基于地震数据的时频分解与气藏检测方法的流程图；
图2是本申请一实施例中的基于地震数据的时频分解与气藏检测方法的流程图；
图3是本申请一实施例中的步骤S202的流程图；
图4是本申请一实施例中的步骤S203的流程图；
图5a是本申请一实施例中的不同阶数(n＝2～5)Paul小波的时间域响应图；
图5b是本申请一实施例中的不同阶数(n＝2～5)Paul小波的频率域响应图；
图6是本申请一实施例中的步骤S205的流程图；
图7是本申请一实施例中的衰减梯度属性提取示意图；
图8a是本申请一实施例中的尺度0.034下4.04阶时间域Paul小波与Ricker子波关系图；
图8b是本申请一实施例中的尺度0.034下4.04阶频率域Paul小波与Ricker子波关系图；
图9a是本申请一实施例中的尺度0.2下24.6阶时间域Paul小波与Ricker子波关系图；
图9b是本申请一实施例中的尺度0.2下24.6阶频率域Paul小波与Ricker子波关系图；
图10是本申请一实施例中的时频谱分辨率对比图；
图11是本申请的一种基于地震数据的时频分解与气藏检测系统的结构示意图；
图12是本申请一实施例中的第一生成单元的结构示意图；
图13是本申请一实施例中的第二生成单元的结构示意图；
图14是本申请一实施例中的互相关单元的结构示意图；
图15是本申请一实施例中的第三生成单元的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
关于本文中所使用的“第一”、“第二”、……等，并非特别指称次序或顺位的意思，亦非用以限定本发明，其仅为了区别以相同技术用语描述的元件或操作。
关于本文中所使用的“包含”、“包括”、“具有”、“含有”等等，均为开放性的用语，即意指包含但不限于。
关于本文中所使用的“及/或”，包括所述事物的任一或全部组合。
针对现有技术中存在的缺陷，本发明提供了一种基于地震数据的时频分解与气藏检测方法，其流程图如图1所示，该方法包括以下步骤：
S101：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
S102：根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
S103：将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱；
S104：根据最优时频谱，生成衰减梯度数据及分频数据；
S105：根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。
由图1所示的流程可知，本申请利用雷尼(Rényi)熵判据获取最优化的STFT参考时频谱，对实数阶Paul小波，通过与参考谱基础互相关获取最优化的Paul-CWT时频谱，根据最优化的时频谱得到分频数据，并利用EAA方法获取稳健的衰减梯度属性剖面，进行高精度的气藏检测，具有提高地震数据频谱分解与衰减属性分析分辨率及抗噪性的有益效果。
为了使本领域的技术人员更好的了解本申请，下面列举一个更为详细的实施例，如图2所示，本实施例提供了一种基于地震数据的时频分解与气藏检测方法，该方法包括以下步骤：
S201：读取地震数据及层位数据。
具体实施时，激发并记录地震波，按常规地震资料进行高保真处理，形成叠后地震数据剖面或三维数据体S(x,y,t)，这里地震资料高保真处理包括几何扩散补偿、地表一致性振幅恢复、地表一致性子波整形反褶积、叠前保幅去噪、折射波静校正、反射波地表一致性剩余静校正、道集规则化处理、动校正、倾角时差校正、并经过合理的偏移叠加技术实现地震资料的偏移归位获得叠后地震数据S(x,y,t)。其中叠后地震数据S(x,y,t)包括若干道地震数据，各道地震数据s(t)相互独立。
对工区地震资料进行基本的层位解释，追踪并闭合储层顶底层位，获得相应层位数据H(x,y,T)。
其中，x代表纵测线坐标，y代表联络测线坐标，t代表时间深度，T代表目的层时间深度坐标。
以下以一道地震数据s(t)的处理过程为例。
S202：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱。
如图3所示，步骤S202包括以下实现步骤：
S301：根据读入的地震数据，利用短时傅里叶变换生成每道地震数据对应的第三时频谱。
具体实施时，以一道地震数据s(t)的处理过程为例，对地震数据s(t)进行短时傅里叶变换，生成地震数据对应的初始时频谱STFT(t,f)。
根据选取的窗函数及时窗长度，对每道地震数据进行短时傅里叶变换生成各第三时频谱。对叠后地震数据S(x,y,t)进行逐道处理，以其中一道地震数据s(t)为例进行阐述，选定STFT的窗函数h(t)类型(一般选Gauss窗)和时窗长度Lh，例如Lh＝60:2:120，根据时窗长度Lh的各个值对地震数据s(t)进行STFT变换，生成各时窗长度Lh值的地震数据s(t)对应的参考时频谱STFT(t,f)，地震数据s(t)进行STFT生成初始时频谱STFT(t,f)的公式如公式(1)所示：
STFT(t,f)＝∫s(τ)h*(τ-t)e-i2πfτdτ (1)
其中，τ的取值范围是实数域(即负无穷至正无穷之间)，f为频率，t为时间，h*(τ-t)e-i2πfτ为核函数，STFT(t,f)为短时傅里叶变换后的初始时频谱。
根据公式(1)生成地震数据s(t)的不同时窗长度Lh值的初始时频谱STFT(t,f)，由各初始时频谱STFT(t,f)组成地震数据s(t)对应的一系列初始时频谱STFT(t,f,Lh)。其中地震数据s(t)的初始时频谱STFT(t,f)即地震数据s(t)对应的第三时频谱TFR(t,f)。
S302：根据各第三时频谱，利用雷尼熵定义生成各第三时频谱对应的雷尼熵，将最小的雷尼熵对应的第三时频谱作为第一时频谱。
具体实施时，Baraniuk等人提出利用广义的Rényi熵估计信号所含信息与时频谱复杂程度的思路。地震数据s(t)对应的第三时频谱TFR(t,f)的Rényi熵定义如公式(2)所示：
其中，α为Rényi熵的阶数，一般取3；底数b一般取2即可。Rényi熵是Shannon熵的广义形式，Shannon熵是Rényi熵在α→1时的极限。相比Shannon熵，具有更高的普适性。Rényi熵值越小，说明信号在时频域的能量复杂度越低，时频聚集性越高；反之，则复杂度越高，能量分布越均匀，聚集性越差。根据公式(3)将聚集性最高的第三时频谱TFR(t,f)作为第一时频谱TFR0(t,f)即参考时频谱：
S203：根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱。
如图4所示，步骤S203具体包括以下步骤：
S401：对每道地震数据进行不同阶数的连续小波变换，生成每道地震数据对应的不同阶数的小波函数。
具体实施时，选取Paul小波，Paul小波的定义与相应的Fourier变换如公式(4)所示：
其中n∈N为正整数，ψn(t)是Paul小波时间域振幅；
是Paul小波频率域响应。
用Gamma函数Γ(n)代替阶乘可以获得正实数域连续阶数变化的Paul小波如公式(5)所示，从而在波形变化上更具一般性和灵活性：
其中n∈R+的正实数。
公式(4)及(5)中，ω是角频率，n是Paul小波的阶数，j是虚数单位，H(ω)是单位阶跃函数。Paul小波阶数从2变化到20，Paul小波阶数变换的步长根据经验值可选取值为0.5或者0.25等。
当n＝2～5时，不同阶数(n＝2～5)Paul小波的时间域与频率域响应，如图5a及5b所示。图5a为不同阶数(n＝2～5)的Paul小波时间域响应图，图5b为不同阶数(n＝2～5)Paul小波的频率域响应图。
对各道地震数据s(t)根据公式(6)进行连续小波变换CWT：
其中，ψu,s(t)是小波核函数，是由Paul小波ψ(t)经过时间平移u与尺度变换s，生成每个道地震数据s(t)的一族不同阶数n的小波函数Wx(u,s)，Wx(u,s)是CWT的计算结果，二维矩阵形式表达的小波系数。
S402：根据不同阶数的小波函数及预设的采样率与记录长度，生成每道地震数据的各阶数对应的第二时频谱。
具体实施时，构建连续小波变换CWT的频率轴：为保证从小波变换中采用的尺度因子对应的尺度坐标轴转换过来的频率轴可以线性变化，读取叠后地震数据S(x,y,t)，同时获取对应的采样率为dt，采样点数为Nt，读取叠后地震数据的记录长度为Nt·dt，则尺度轴Scale如公式(7)所示：
S204：将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱。
具体实施时，首先将地震数据s(t)的第一时频谱TFR0(t,f)与该道地震数据s(t)的每个阶数对应的第二时频谱TFR(t,f,n)做互相关运算，输出地震数据s(t)的各第二时频谱TFR(t,f,n)对应的相关值。
其次，从该道地震数据s(t)的每个阶数对应的互相关值中选取出最大互相关值，并将该道地震数据s(t)的最大值互相关值对应的第二时频谱TFR(t,f,n)作为该道地震数据s(t)的最优时频谱TFR(t,f)。
具体计算公式如公式(8)所示：
TFR(t,f)＝TFR(t,f,n)|cov<TFR(t,f,n),TFR0(t,f)>→max (8)
S205：根据最优时频谱，生成衰减梯度数据及分频数据。
如图6所示，步骤S205具体包括：
S501：根据每道地震数据对应的最优时频谱，利用吸能分析法生成衰减梯度数据。
具体实施时，对每道地震数据s(t)的最优时频谱TFR(t,f)采用吸能分析(Energyabsorption analysis，EAA)方法，即认为地震数据t0时刻瞬时最优时频谱TFR(t,f)具有指数衰减的形式，如公式(9)所示：
TFR(t0,f)＝A·e-αf (9)
其中，α是衰减梯度数据，A是常数。上式两边取对数，则为公式(10)：
lnTFR(t0,f)＝lnA-αf (10)
如图7所示，根据拟合瞬时振幅谱能量占比60％的频点E60％和90％的频点E90％两个位置的斜率，从而求得衰减梯度数据α，如公式(11)：
S502：根据每道地震数据对应的最优时频谱，生成一分频数据。
根据步骤S201到步骤S203，处理叠后地震数据S(x,y,t)的各道地震数据s(t)，出于效率考虑，由于各道地震数据s(t)处理相互独立，可以采用并行处理，提高处理效率。选定一系列频点f，根据各道地震数据s(t)的最优时频谱TFR(t,f)生成最终的一个分频数据TFR(x,y,t,f)，其中f是选定的频点向量，按这些离散频率点f提取分频数据TFR(x,y,t,f)，在程序中通过频率轴网格插值完成。
S206：根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。
具体实施时，分频数据TFR(x,y,t,f)结合步骤S201读取的层位数据H(x,y,T)，可以提取不同频段分频振幅体和沿层切片，利用沿层切片进行气藏检测，提高检测精度。
本申请提供的一种基于地震数据的时频分解与气藏检测方法，利用雷尼(Rényi)熵判据获取最优化的STFT第一时频谱，利用不同阶数的连续小波变换生成第二时频谱，通过将第一时频谱与第二时频谱进行互相关，获取最优化的Paul-CWT最优时频谱，并根据最优时频谱得到分频数据，同时利用EAA方法获取衰减梯度数据，进行高精度的气藏检测，具有提高分频数据与能量吸收衰减特征的属性剖面数据的精度的有益效果。即实数阶Paul小波的连续小波变换可以任意拟合Ricker子波与Morlet小波，其时频谱在低频端具有更高的频率分辨率，在高频端具有更高的时间分辨率，同时基于Rényi熵判据的阶数优化避免了阶数选择的主观偏差，以此为基础的衰减特征剖面数据则具有更高的抗噪性能。因此，本申请具有提高地震数据频谱分解与衰减属性分析分辨率及抗噪性的有益效果。
为了更加清晰的体现本申请的效果，本申请提供了具体实施中的部分结果：
图8a及8b分别从时间域与频率域展示Paul小波与地震勘探中常用的Ricker子波之间的关系，尺度0.034下4.04阶Paul小波几乎相当于尺度为1峰值频率20Hz的Ricker子波，图8a是尺度0.034下4.04阶时间域Paul小波与Ricker子波关系图，图8b是尺度0.034下4.04阶频率域Paul小波与Ricker子波关系图。
图9a及9b分别从时间域与频率域展示Paul小波与地震勘探中常用的Morlet子波之间的关系，尺度0.2下24.6阶的Paul小波几乎相当于尺度0.04峰值频率20Hz的Morlet子波，图9a是尺度0.2下24.6阶时间域Paul小波与Ricker子波关系图，图9b是尺度0.2下24.6阶频率域Paul小波与Ricker子波关系图。
图10中对不同时频分析方法的时频谱分辨率进行了对比，其中(a)Ricker子波道集；(b)道集叠加合成信号；(c)STFT(100ms时窗)；(d)基于32阶Gauss导数小波的CWT；(e)基于Bump小波(μ＝5,σ＝1.2)的CWT；(f)基于Morlet小波(ω0＝5.8)的CWT；(g)基于6阶Paul小波的CWT。
基于与上述基于地震数据的时频分解与气藏检测方法相同的申请构思，本发明还提供了一种基于地震数据的时频分解与气藏检测系统，如下面实施例所述。由于该系统解决问题的原理与基于地震数据的时频分解与气藏检测方法相似，因此该基于地震数据的时频分解与气藏检测系统的实施可以参见基于地震数据的时频分解与气藏检测方法的实施，重复之处不再赘述。
图11为本申请实施例的一种基于地震数据的时频分解与气藏检测系统的结构示意图，如图11所示，该基于地震数据的时频分解与气藏检测系统包括：第一生成单元101、第二生成单元102、互相关单元103、第三生成单元104及气藏检测单元105。
第一生成单元101，用于根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱。其中由若干道地震数据组成叠后地震数据。
第二生成单元102，用于根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱。
互相关单元103，用于将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱。
第三生成单元104，用于根据最优时频谱，生成衰减梯度数据及分频数据。
气藏检测单元105，用于根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。
在一个实施例中，如图12所示，第一生成单元101包括：傅里叶变换模块201及第一时频谱输出模块202。
傅里叶变换模块201，用于根据读入的地震数据，利用短时傅里叶变换生成每道地震数据对应的第三时频谱。
第一时频谱输出模块202，用于根据各第三时频谱，利用雷尼熵定义生成各第三时频谱对应的雷尼熵，将最小的雷尼熵对应的第三时频谱作为第一时频谱。
在一个实施例中，傅里叶变换模块201，具体用于根据选取的窗函数及时窗长度，对每道地震数据进行短时傅里叶变换生成各第三时频谱。
在一个实施例中，如图13所示，第二生成单元102包括：连续小波变换模块301及第二时频谱生成模块302。
连续小波变换模块301，用于对每道地震数据进行不同阶数的连续小波变换，生成每道地震数据对应的不同阶数的小波函数。
第二时频谱生成模块302，用于根据不同阶数的小波函数及预设的采样率与记录长度，生成每道地震数据的各阶数对应的第二时频谱。
在一个实施例中，如图14所示，互相关单元103包括：相关值输出模块401及最优时频谱输出模块402。
相关值输出模块401，用于将每道地震数据的对应的第一时频谱分别与每道地震数据的对应的各第二时频谱进行互相关运算，输出每道地震数据的各第二时频谱对应的相关值。
最优时频谱输出模块402，用于将每道地震数据的最大的相关值对应的第二时频谱作为最优时频谱，输出每道地震数据对应的最优时频谱。
在一个实施例中，如图15所示，第三生成单元104包括：衰减梯度数据生成模块501及分频数据生成模块502。
衰减梯度数据生成模块501，用于根据每道地震数据对应的最优时频谱，利用吸能分析法生成衰减梯度数据。
分频数据生成模块502，用于根据每道地震数据对应的最优时频谱，生成一分频数据。
基于与上述基于地震数据的时频分解与气藏检测方法相同的申请构思，本申请提供一种计算机设备，如下面实施例所述。由于该计算机设备解决问题的原理与基于地震数据的时频分解与气藏检测方法相似，因此该计算机设备的实施可以参见基于地震数据的时频分解与气藏检测方法的实施，重复之处不再赘述。
在一个实施例中，计算机设备包括：存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，如图1所示，所述处理器执行所述计算机程序时实现以下步骤：
S101：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
S102：根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
S103：将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱；
S104：根据最优时频谱，生成衰减梯度数据及分频数据；
S105：根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。
基于与上述基于地震数据的时频分解与气藏检测方法相同的申请构思，本申请提供一种计算机可读存储介质，如下面实施例所述。由于该计算机可读存储介质解决问题的原理与基于地震数据的时频分解与气藏检测方法相似，因此该计算机可读存储介质的实施可以参见基于地震数据的时频分解与气藏检测方法的实施，重复之处不再赘述。
在一个实施例中，计算机可读存储介质上存储有计算机程序，如图1所示，该计算机程序被处理器执行时实现以下步骤：
S101：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
S102：根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
S103：将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱；
S104：根据最优时频谱，生成衰减梯度数据及分频数据；
S105：根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。
本发明提供的一种基于地震数据的时频分解与气藏检测方法及系统，包括：根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；将第一时频谱分别与各第二时频谱进行互相关运算，生成最优时频谱；根据最优时频谱，生成衰减梯度数据及分频数据；根据读入的层位数据、衰减梯度数据及分频数据，生成沿层切片并进行气藏检测。本申请利用雷尼(Rényi)熵判据获取最优化的STFT第一时频谱，利用不同阶数的连续小波变换生成第二时频谱，通过将第一时频谱与第二时频谱进行互相关，获取最优化的Paul-CWT最优时频谱，并根据最优时频谱得到分频数据，同时利用EAA方法获取衰减梯度数据，进行高精度的气藏检测，具有提高地震数据频谱分解与衰减属性分析分辨率及抗噪性的有益效果。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
Claims (8)
1.一种基于地震数据的时频分解与气藏检测方法，其特征在于，包括：
根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱，具体包括：
对每道地震数据进行不同阶数的连续小波变换，生成每道地震数据对应的不同阶数的小波函数；
根据所述不同阶数的小波函数及预设的采样率与记录长度，生成每道地震数据的各阶数对应的第二时频谱；
将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱，具体包括：
将每道地震数据的对应的所述第一时频谱分别与每道地震数据的对应的各所述第二时频谱进行互相关运算，输出每道地震数据的各所述第二时频谱对应的相关值；
将每道地震数据的最大的相关值对应的第二时频谱作为最优时频谱，输出每道地震数据对应的最优时频谱；
根据所述最优时频谱，生成衰减梯度数据及分频数据，具体包括：
根据每道地震数据对应的最优时频谱，利用吸能分析法生成衰减梯度数据；
根据每道地震数据对应的最优时频谱，生成分频数据；
根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测；
其中，根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱，具体包括：
根据读入的地震数据，利用短时傅里叶变换生成每道地震数据对应的第三时频谱；根据各所述第三时频谱，利用雷尼熵定义生成各所述第三时频谱对应的雷尼熵，将最小的雷尼熵对应的第三时频谱作为所述第一时频谱。
2.根据权利要求1所述的基于地震数据的时频分解与气藏检测方法，其特征在于，由若干道所述地震数据组成叠后地震数据。
3.根据权利要求2所述的基于地震数据的时频分解与气藏检测方法，其特征在于，所述利用短时傅里叶变换生成每道地震数据对应的第三时频谱，包括：根据选取的窗函数及时窗长度，对每道地震数据进行短时傅里叶变换生成各所述第三时频谱。
4.一种基于地震数据的时频分解与气藏检测系统，其特征在于，包括：
第一生成单元，用于根据读入的地震数据，利用短时傅里叶变换及雷尼熵定义生成第一时频谱；
第二生成单元，用于根据读入的地震数据，利用不同阶数的连续小波变换生成各阶数对应的第二时频谱；
所述第二生成单元包括：
连续小波变换模块，用于对每道地震数据进行不同阶数的连续小波变换，生成每道地震数据对应的不同阶数的小波函数；
第二时频谱生成模块，用于根据所述不同阶数的小波函数及预设的采样率与记录长度，生成每道地震数据的各阶数对应的第二时频谱；
互相关单元，用于将所述第一时频谱分别与各所述第二时频谱进行互相关运算，生成最优时频谱；
所述互相关单元包括：
相关值输出模块，用于将每道地震数据的对应的所述第一时频谱分别与每道地震数据的对应的各所述第二时频谱进行互相关运算，输出每道地震数据的各所述第二时频谱对应的相关值；
最优时频谱输出模块，用于将每道地震数据的最大的相关值对应的第二时频谱作为最优时频谱，输出每道地震数据对应的最优时频谱；
第三生成单元，用于根据所述最优时频谱，生成衰减梯度数据及分频数据；
所述第三生成单元包括：
衰减梯度数据生成模块，用于根据每道地震数据对应的最优时频谱，利用吸能分析法生成衰减梯度数据；
分频数据生成模块，用于根据每道地震数据对应的最优时频谱，生成分频数据；
气藏检测单元，用于根据读入的层位数据、所述衰减梯度数据及所述分频数据，生成沿层切片并进行气藏检测；
其中，所述第一生成单元具体包括：
傅里叶变换模块，用于根据读入的地震数据，利用短时傅里叶变换生成每道地震数据对应的第三时频谱；
第一时频谱输出模块，用于根据各所述第三时频谱，利用雷尼熵定义生成各所述第三时频谱对应的雷尼熵，将最小的雷尼熵对应的第三时频谱作为所述第一时频谱。
5.根据权利要求4所述的基于地震数据的时频分解与气藏检测系统，其特征在于，由若干道所述地震数据组成叠后地震数据。
6.根据权利要求5所述的基于地震数据的时频分解与气藏检测系统，其特征在于，所述傅里叶变换模块，具体用于根据选取的窗函数及时窗长度，对每道地震数据进行短时傅里叶变换生成各所述第三时频谱。
7.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1所述的基于地震数据的时频分解与气藏检测方法的步骤。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1所述的基于地震数据的时频分解与气藏检测方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810801169.5A CN109001800B (zh) | 2018-07-20 | 2018-07-20 | 一种基于地震数据的时频分解与气藏检测方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810801169.5A CN109001800B (zh) | 2018-07-20 | 2018-07-20 | 一种基于地震数据的时频分解与气藏检测方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109001800A CN109001800A (zh) | 2018-12-14 |
| CN109001800B true CN109001800B (zh) | 2020-03-10 |
Family
ID=64596504
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810801169.5A Active CN109001800B (zh) | 2018-07-20 | 2018-07-20 | 一种基于地震数据的时频分解与气藏检测方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109001800B (zh) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11269101B2 (en) * | 2019-04-16 | 2022-03-08 | Saudi Arabian Oil Company | Method and system of direct gas reservoir detection using frequency slope |
| CN112305612B (zh) * | 2019-07-23 | 2022-05-20 | 中国海洋石油集团有限公司 | 高分辨率复谱分解时频空间域振幅随偏移距变化校正方法 |
| CN112711070B (zh) * | 2019-10-24 | 2024-02-20 | 中国石油化工股份有限公司 | 一种基于地震信号分解的油气检测方法及装置 |
| CN113589364B (zh) * | 2020-04-30 | 2023-04-28 | 中国石油化工股份有限公司 | 基于佐布里兹方程约束的地震数据规则化处理方法 |
| US11977198B2 (en) | 2020-10-06 | 2024-05-07 | Saudi Arabian Oil Company | Isofrequency volumes ratio workflow to detect gas reservoirs in 3D domain |
| US11333780B2 (en) | 2020-10-09 | 2022-05-17 | Saudi Arabian Oil Company | Method and system for processing a three-dimensional (3D) seismic dataset |
| US11592589B2 (en) | 2021-01-14 | 2023-02-28 | Saudi Arabian Oil Company | Seismic attribute map for gas detection |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2003900418A0 (en) * | 2003-01-30 | 2003-02-13 | Qr Sciences Limited | Improvements in Signal Processing For Detection Of NQR Signals |
| US7751277B2 (en) * | 2008-03-17 | 2010-07-06 | Pgs Geophysical As | Method for interpolating seismic data by anti-alias, anti-leakage Fourier transform |
| US9507042B2 (en) * | 2012-08-31 | 2016-11-29 | Lumina Geophysical LLC | System and method for constrained least-squares spectral processing and analysis of seismic data |
| CN103235339A (zh) * | 2013-04-09 | 2013-08-07 | 中国石油大学(北京) | 一种时频分解地震流体识别方法 |
| CN104458170B (zh) * | 2014-11-07 | 2017-01-11 | 桂林电子科技大学 | 机械装备监测振动信号的时频图处理方法及系统 |
| CN106405639B (zh) * | 2015-07-30 | 2018-11-23 | 中国石油化工股份有限公司 | 一种叠前地震储层岩性参数的反演方法 |
-
2018
- 2018-07-20 CN CN201810801169.5A patent/CN109001800B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109001800A (zh) | 2018-12-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109001800B (zh) | 一种基于地震数据的时频分解与气藏检测方法及系统 | |
| Chen et al. | Empirical low-rank approximation for seismic noise attenuation | |
| Gómez et al. | A simple method inspired by empirical mode decomposition for denoising seismic data | |
| Wang et al. | Seismic sparse-spike deconvolution via Toeplitz-sparse matrix factorization | |
| Margrave et al. | Gabor deconvolution: Estimating reflectivity by nonstationary deconvolution of seismic data | |
| US8280695B2 (en) | Method to adapt a template dataset to a target dataset by using curvelet representations | |
| Yuan et al. | Stable inversion-based multitrace deabsorption method for spatial continuity preservation and weak signal compensation | |
| Van der Baan et al. | Robust wavelet estimation and blind deconvolution of noisy surface seismics | |
| US8103453B2 (en) | Method of seismic data interpolation by projection on convex sets | |
| CN104849756B (zh) | 一种提高地震数据分辨率增强有效弱信号能量的方法 | |
| Wang et al. | Enhancing resolution of nonstationary seismic data by molecular-Gabor transform | |
| US20180306940A1 (en) | Method for Generating Optimized Seismic Target Spectrum | |
| WO2000042448A1 (en) | Method of attenuating noise in three dimensional seismic data using a projection filter | |
| Zhou et al. | Sparse dictionary learning for seismic noise attenuation using a fast orthogonal matching pursuit algorithm | |
| Chen et al. | Five-dimensional seismic data reconstruction using the optimally damped rank-reduction method | |
| Iqbal et al. | Detection and denoising of microseismic events using time–frequency representation and tensor decomposition | |
| WO2017136133A1 (en) | Efficient seismic attribute gather generation with data synthesis and expectation method | |
| Wu et al. | Fast principal component analysis for stacking seismic data | |
| Zhou et al. | A hybrid method for noise suppression using variational mode decomposition and singular spectrum analysis | |
| Bai et al. | Seismic signal enhancement based on the low‐rank methods | |
| Aeron et al. | Broadband dispersion extraction using simultaneous sparse penalization | |
| Yang et al. | A seismic interpolation and denoising method with curvelet transform matching filter | |
| Ridsdill-Smith | The application of the wavelet transform to the processing of aeromagnetic data | |
| Xin et al. | Seismic high-resolution processing method based on spectral simulation and total variation regularization constraints | |
| Cheng et al. | Multi-trace nonstationary sparse inversion with structural constraints |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |