CN114428282B - 一种基于去尺度s变换的地震信号时频变换方法 - Google Patents
一种基于去尺度s变换的地震信号时频变换方法 Download PDFInfo
- Publication number
- CN114428282B CN114428282B CN202210093730.5A CN202210093730A CN114428282B CN 114428282 B CN114428282 B CN 114428282B CN 202210093730 A CN202210093730 A CN 202210093730A CN 114428282 B CN114428282 B CN 114428282B
- Authority
- CN
- China
- Prior art keywords
- time
- frequency
- seismic
- downscaling
- transform
- 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
- 238000000034 method Methods 0.000 title claims abstract description 29
- 238000006243 chemical reaction Methods 0.000 title claims abstract description 28
- 230000009466 transformation Effects 0.000 claims abstract description 22
- 230000001131 transforming effect Effects 0.000 claims abstract description 8
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 238000012937 correction Methods 0.000 claims description 4
- 238000009792 diffusion process Methods 0.000 claims description 2
- 230000003068 static effect Effects 0.000 claims description 2
- 230000009897 systematic effect Effects 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 abstract description 25
- 238000004458 analytical method Methods 0.000 abstract description 17
- 230000000694 effects Effects 0.000 abstract description 9
- 238000011426 transformation method Methods 0.000 abstract description 4
- 238000000354 decomposition reaction Methods 0.000 description 6
- 230000003595 spectral effect Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- YBJHBAHKTGYVGT-ZKWXMUAHSA-N (+)-Biotin Chemical compound N1C(=O)N[C@@H]2[C@H](CCCCC(=O)O)SC[C@@H]21 YBJHBAHKTGYVGT-ZKWXMUAHSA-N 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000001125 extrusion Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- FEPMHVLSLDOMQC-UHFFFAOYSA-N virginiamycin-S1 Natural products CC1OC(=O)C(C=2C=CC=CC=2)NC(=O)C2CC(=O)CCN2C(=O)C(CC=2C=CC=CC=2)N(C)C(=O)C2CCCN2C(=O)C(CC)NC(=O)C1NC(=O)C1=NC=CC=C1O FEPMHVLSLDOMQC-UHFFFAOYSA-N 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
本发明提供一种基于去尺度S变换的地震信号时频变换方法，包括：获取地震信号数据，并进行预处理，获取地震时间序列数据；根据去尺度S变换方法对地震时间序列数据进行变换，获取地震信号变换后的时频表示结果；根据时频表示结果计算获取地震属性参数。该方法针对S变换存在的问题，通过去除尺度因子进行修正，消除了S变换对宽频带地震信号时频分析的误导作用，同时保留S变换的多分辨率特征和与傅里叶频谱的直接关系。
Description
技术领域
本发明涉及地震勘探技术领域，具体涉及一种基于去尺度S变换的地震信号时频变换方法。
背景技术
地震信号本质上是非平稳信号，其频率成分含量随时间记录的变化而变化，蕴藏着丰富的地质信息。谱分解(也称为时频分析)是处理非平稳信号的强有力工具。它可以将一维地震道转换到二维的时间—频率域，以描述地震信号频率分量，揭示地震频率的时变特性与地震资料的特征。时频分析能够突出宽频带地震信号中隐藏的与频率相关的地球物理响应，而这些响应通常与频率异常现象直接相关，因此在实际地震数据处理和解释中得到了广泛的应用。时频分析方法计算出来的属性可用于描述地质体结构，如薄互层、河道、断层和构造异常等，其在地震资料去噪、直接油气指示，地震衰减测量，薄层反射率反演，地层非均质性测定，油藏相照明，孔隙和渗透率分布映射等方面也有着广泛的应用。
由于谱分解本质上是一个非唯一的过程，因此在实际地震勘探中使用了多种时频分析的方法：包括短时傅里叶变换，连续小波变换，S变换，匹配追踪分解，二次型时频分布方法魏格纳—威尔分布，经验模态分解，同步挤压变换，和正则化谱反演方法等。这些不同的谱分解方法基于不同的数学原理，产生不同的时频表示，适用于特定的应用。
其中，S变换结合了连续小波变换(CWT)和短时傅里叶变换(STFT)的特点，以基本小波为固定窗函数，并引入信号频率这一参数来控制窗口大小，即实现了分辨率自适应调节，又增强了高频弱振幅信号的能量，时频分析效果优于两者，且反变换无损可逆，具有良好的计算速度，并且易于实现。S 变换已经广泛应用于不同的实际应用，如时频滤波、噪声衰减、薄层估计和衰减介质地震反褶积等。
然而，Li等人(2016)提出，用信号的傅里叶变换频率响应来解释S 变换的振幅谱可能并不容易，并且提出通过将频率f替换成线性函数或其他函数来提高分辨率。以此为突破口，通过研究发现，S变换通过乘以一个尺度因子，即频率的绝对值，将与频率相关的高斯窗函数归一化处理。这种归一化处理使S变换保持了时间上的能量守恒，却引起了频谱的改变，使得振幅峰值向更高的频率移动。在给定的时间内，S变换的振幅谱与输入信号波形导数的傅里叶变换振幅谱具有相似的特征。对于带宽窄的信号，这种归一化对时频分析的峰值频率影响不大。然而，对于宽频带地震信号，如Ricker 子波，其S变换振幅谱峰值频率明显高于其傅里叶变换峰值频率，因而具有误导性。经进一步研究发现，各种波形的峰值频率数值比较结果支持了S变换峰值频率等于或大于傅里叶变换峰值频率的一般规律。实际地震资料的时间—频率域分析结果表明，在解释S变换频谱分解时应考虑这种效应。
为此，本发明提出了一种新的基于去尺度S变换的地震信号时频变换方法。
发明内容
为解决上述问题，本发明的目的在于提供一种基于去尺度S变换的地震信号时频变换方法，该方法针对S变换存在的问题，通过去除尺度因子进行修正，消除了S变换对宽频带地震信号时频分析的误导作用，同时保留S变换的多分辨率特征和与傅里叶频谱的直接关系。
为实现上述目的，本发明提供了如下的技术方案。
一种基于去尺度S变换的地震信号时频变换方法，包括以下步骤：
获取地震时间序列数据；
根据去尺度S变换对地震时间序列数据进行变换，获取地震信号变换后的时频表示结果；所述去尺度S变换的表达式：
其中，h(t)为时间序列数据；为高斯窗函数；e-i2πft为核函数；τ为高斯窗函数的中间时刻；f为频率。
优选地，还包括：对地震信号数据进行预处理；所述地震信号数据预处理包括以下步骤：对述地震信号数据进行数据解编、格式转换、编辑、几何扩散校正、建立野外观测系统和静校正。
优选地，还包括：根据时频表示结果计算获取地震属性参数。
优选地，所述地震属性参数包括平均频率g(t)：
式中p(t,f)对应时频分布。
优选地，所述地震属性参数还包括瞬时带宽B(t)：
优选地，所述地震属性参数还包括扭度St：
优选地，所述地震属性参数还包括峰度kt：
本发明有益效果：
本发明提出一种基于去尺度S变换的地震信号时频变换方法，该方法针对S变换存在的问题，通过去除尺度因子进行修正，消除了S变换对宽频带地震信号时频分析的误导作用，同时保留S变换的多分辨率特征和与傅里叶频谱的直接关系。相比之下，去尺度S变换分辨率有所提高，能更准确的刻画信号的低频振幅分量，获得更高分辨率的时频谱。这种对S变换方法的有效改进能够适应地震勘探数据处理和定量解释的需要。
附图说明
图1是本发明实施例的方法流程图；
图2a是本发明实施例1的主频为80HZ Ricker子波信号图；
图2b是本发明实施例1的Ricker子波S变换的时频分布图；
图2c是本发明实施例1的Ricker子波去尺度S变换的时频分布图；
图2d是本发明实施例1的t＝0.5s时刻，Ricker子波S变换振幅谱(ST)、去尺度S变换振幅谱(UST)、傅里叶变换振幅谱(FT)和Ricker 子波导数的傅氏变换振幅谱(FTD)；
图3a是本发明实施例1的主频为80HZMorlet小波信号图；
图3b是本发明实施例1的Morlet小波S变换的时频分布图；
图3c是本发明实施例1的Morlet小波去尺度S变换的时频分布图；
图3d是本发明实施例1的t＝0.5s时刻，Morlet小波S变换振幅谱 (ST)、去尺度S变换振幅谱(UST)、傅里叶变换振幅谱(FT)和Morlet 小波导数的傅氏变换振幅谱(FTD)；
图4是本发明实施例2的墨西哥湾某工区2D地震剖面数据；
图5a是本发明实施例2的墨西哥湾近海某工区三维勘探中储层上方记录的地震道；
图5b是本发明实施例2的采用长度200ms，半极值处宽度为20(1/e采样间隔)的高斯窗函数的短时傅里叶变换(STFT)的时频表示图；
图5c是本发明实施例2的S变换的时频表示图；
图5d是本发明实施例2的去尺度S变换的时频表示图；
图5e是本发明实施例2的当地震振幅大于0.02的最大值时，S变换 (ST)、去尺度S变换(UST)和短时傅里叶变换(STFT)的峰值频率对比图；
图5f是本发明实施例2的t＝0.8s时，对应的储层时间位置的S变换 (ST)、去尺度S变换(UST)和短时傅里叶变换(STFT)的振幅谱的对比图；
图6a是本发明实施例2通过S变换计算的瞬时频率属性剖面图；
图6b是本发明实施例2通过去尺度S变换计算的瞬时频率属性剖面图；
图7a是本发明实施例2通过S变换计算的带宽属性剖面图；
图7b是本发明实施例2通过去尺度S变换计算的带宽属性剖面图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例1
一种基于去尺度S变换的地震信号时频变换方法，方法流程图如图1所示，包括以下步骤：
获取地震信号数据，并进行预处理，获取地震时间序列数据；
根据去尺度S变换方法对地震时间序列数据进行变换，获取地震信号变换后的时频表示结果；所述去尺度S变换的表达式：
其中，h(t)为时间序列数据；为高斯窗函数；e-i2πft为核函数；
根据时频表示结果计算获取时频地震属性参数。
本实施例计算了有意义的小波如雷克子波和Morlet小波的S变换和去尺度S变换结果，比较了不同子波中心时间的时频谱与傅里叶频谱之间的差异；
从图2、图3结果可以看出，S变换中的尺度因子将频谱向更高的频率倾斜，从而对信号的时频分析结果产生误导。对于带宽窄的信号，如正弦信号或Morlet小波，S变换的影响很小或可忽略。然而，对于宽频带信号，与地震勘探应用相关的，如Ricker子波，其影响是相当大的。S变换尺度因子是一种低切滤波器，它缩小了信号带宽，并且使峰值频率往高处移动，会被误解为提高了频率分辨率，并将会干扰时频地震属性的计算，如峰值频率、带宽、扭度和峰度。
实施例2
采用长度200ms，半宽20(1/e采样间隔)的高斯窗函数的短时傅里叶变换 (STFT)，S变换(ST)，去尺度因子S变换(UST)，对来自墨西哥海湾的地震数据进行时频分析，用于地震数据中出现的“亮点”和低频阴影的识别方面，如图4所示，并进行定性分辨率比较。
从地震数据时频分析(图5)的结果可以看出：S变换时频谱在低频时能量偏低，峰值频率往高频移动。在t＝0.8s,低阻抗层的情况下(出现强波谷和强波峰)，比较S变换(ST)，去尺度S变换(UST)和短时傅里叶变换 (STFT)方法的振幅谱，可以看到去尺度S变换(UST)和短时傅里叶变换 (STFT)时频谱具有相似的带宽和峰值频率，而S变换(ST)峰值频率偏高。通过去除尺度因子，发现在解析波形和真实地震数据上，去尺度因子的S变换将移向真实的频率并具有类似的带宽。
在对实际地震数据进行时频分析的基础上，结合峰值频率、带宽等时频地震属性，测试去尺度S变换方法的有效性，及其相对于S变换的分辨率提高和改进程度。
对墨西哥湾工区的2D地震资料(图4)，计算频率属性剖面(图6)和带宽属性剖面(图7)，用以比较S变换(ST)和去尺度S变换的频率属性。从两种时频地震属性剖面图可以看出，S变换(ST)和去尺度S变换 (UST)的频率属性剖面有很大的不同。去尺度S变换(UST)峰值频率的频率始终较低，而频谱带宽始终较大。对于同一反射层位(白色虚线部分)，两种时频分析方法的频率属性显示不同的横向变化，这会导致对岩性、流体和地层参数的不同解释。如果作为解释属性，不需要正确地确定峰值频率，但是如果峰值频率被用于定量地震解释，如储层厚度的预测或指导地震反演，则必须考虑S变换对宽频带地震信号的误导作用。
以上仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于去尺度S变换的地震信号时频变换方法，其特征在于，包括以下步骤：
获取地震时间序列数据；
根据去尺度S变换对地震时间序列数据进行变换，获取地震信号变换后的时频表示结果；所述去尺度S变换的表达式：
其中，h(t)为时间序列数据；
为高斯窗函数；e-i2πft为核函数；τ为高斯窗函数的中间时刻；f为频率。
2.根据权利要求1所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，还包括：对地震信号数据进行预处理；所述地震信号数据预处理包括以下步骤：对述地震信号数据进行数据解编、格式转换、编辑、几何扩散校正、建立野外观测系统和静校正。
3.根据权利要求1所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，还包括：根据时频表示结果计算获取地震属性参数。
4.根据权利要求3所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，所述地震属性参数包括平均频率g(t)：
式中p(t,f)对应时频分布。
5.根据权利要求4所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，所述地震属性参数还包括瞬时带宽B(t)：
6.根据权利要求5所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，所述地震属性参数还包括扭度St：
7.根据权利要求5所述的基于去尺度S变换的地震信号时频变换方法，其特征在于，所述地震属性参数还包括峰度kt：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210093730.5A CN114428282B (zh) | 2022-01-26 | 2022-01-26 | 一种基于去尺度s变换的地震信号时频变换方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210093730.5A CN114428282B (zh) | 2022-01-26 | 2022-01-26 | 一种基于去尺度s变换的地震信号时频变换方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114428282A CN114428282A (zh) | 2022-05-03 |
| CN114428282B true CN114428282B (zh) | 2023-04-25 |
Family
ID=81313239
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210093730.5A Active CN114428282B (zh) | 2022-01-26 | 2022-01-26 | 一种基于去尺度s变换的地震信号时频变换方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114428282B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101939763A (zh) * | 2007-12-10 | 2011-01-05 | 爱克发医疗保健公司 | 生成多尺度对比度增强的图像的方法 |
| WO2013107076A1 (zh) * | 2012-01-19 | 2013-07-25 | 东南大学 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
| CN107272063A (zh) * | 2017-07-05 | 2017-10-20 | 西安交通大学 | 基于高分辨率时频分析和一致性度量的非均质性刻画方法 |
| CN107436450A (zh) * | 2017-07-26 | 2017-12-05 | 西安交通大学 | 一种基于连续小波变换的地震信号带宽拓展方法 |
| CN109239780A (zh) * | 2017-07-10 | 2019-01-18 | 中国石油化工股份有限公司 | 基于同步挤压小波变换去除面波的方法 |
| CN109343020A (zh) * | 2018-11-16 | 2019-02-15 | 南京理工大学 | 一种基于改进窗函数的s变换时频分析方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR2933499B1 (fr) * | 2008-07-03 | 2010-08-20 | Inst Francais Du Petrole | Methode d'inversion conjointe de donnees sismiques representees sur des echelles de temps differentes |
| 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 |
-
2022
- 2022-01-26 CN CN202210093730.5A patent/CN114428282B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101939763A (zh) * | 2007-12-10 | 2011-01-05 | 爱克发医疗保健公司 | 生成多尺度对比度增强的图像的方法 |
| WO2013107076A1 (zh) * | 2012-01-19 | 2013-07-25 | 东南大学 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
| CN107272063A (zh) * | 2017-07-05 | 2017-10-20 | 西安交通大学 | 基于高分辨率时频分析和一致性度量的非均质性刻画方法 |
| CN109239780A (zh) * | 2017-07-10 | 2019-01-18 | 中国石油化工股份有限公司 | 基于同步挤压小波变换去除面波的方法 |
| CN107436450A (zh) * | 2017-07-26 | 2017-12-05 | 西安交通大学 | 一种基于连续小波变换的地震信号带宽拓展方法 |
| CN109343020A (zh) * | 2018-11-16 | 2019-02-15 | 南京理工大学 | 一种基于改进窗函数的s变换时频分析方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114428282A (zh) | 2022-05-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111208561B (zh) | 基于时变子波与曲波变换约束的地震声波阻抗反演方法 | |
| Chen | Fast dictionary learning for noise attenuation of multidimensional seismic data | |
| US8538702B2 (en) | Geologic features from curvelet based seismic attributes | |
| Gómez et al. | A simple method inspired by empirical mode decomposition for denoising seismic data | |
| Liu et al. | A 1D time-varying median filter for seismic random, spike-like noise elimination | |
| Li et al. | Investigation of generalized S-transform analysis windows for time-frequency analysis of seismic reflection data | |
| CN109001800B (zh) | 一种基于地震数据的时频分解与气藏检测方法及系统 | |
| CN108020863A (zh) | 一种基于地震奇偶函数的碳酸盐岩薄储层孔隙度预测方法 | |
| Cai et al. | Seismic data denoising based on mixed time-frequency methods | |
| CN113777650B (zh) | 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 | |
| Zoukaneri et al. | A combined Wigner-Ville and maximum entropy method for high-resolution time-frequency analysis of seismic data | |
| Zhou et al. | Sparse dictionary learning for seismic noise attenuation using a fast orthogonal matching pursuit algorithm | |
| Zhong‐Lai et al. | A SECOND‐ORDER SYNCHROSQUEEZING S‐TRANSFORM AND ITS APPLICATION IN SEISMIC SPECTRAL DECOMPOSITION | |
| Aeron et al. | Broadband dispersion extraction using simultaneous sparse penalization | |
| CN113281809B (zh) | 一种地震信号的谱分析方法 | |
| CN111505709B (zh) | 一种基于稀疏谱分解的衰减定性分析的方法 | |
| Khakhki et al. | High-resolution time-frequency hilbert transform using sparsity-aware weighting function | |
| CN114428282B (zh) | 一种基于去尺度s变换的地震信号时频变换方法 | |
| Kazemnia Kakhki et al. | Seismic attributes via robust and high-resolution seismic complex trace analysis | |
| CN111366978A (zh) | 一种基于多次挤压小波变换的地震时频分析方法及系统 | |
| CN107678065B (zh) | 提高地震分辨率的保构造井控空间反褶积方法和装置 | |
| CN110542921B (zh) | 一种基于主分量分析的振幅横向变化刻画方法 | |
| Wang et al. | Structure-oriented edge-preserving smoothing based on accurate estimation of orientation and edges | |
| CN114152981B (zh) | 一种势能波函数域地震资料品质因子估计方法 | |
| Bui et al. | Time-frequency sparse Gabor transform for detecting microseismic events |
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 |