CN108072894A - 一种地震信号处理方法 - Google Patents
一种地震信号处理方法 Download PDFInfo
- Publication number
- CN108072894A CN108072894A CN201611071443.5A CN201611071443A CN108072894A CN 108072894 A CN108072894 A CN 108072894A CN 201611071443 A CN201611071443 A CN 201611071443A CN 108072894 A CN108072894 A CN 108072894A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- signal
- seismic
- amplitude
- 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.)
- Pending
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 17
- 238000000034 method Methods 0.000 claims abstract description 18
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 8
- 230000002708 enhancing effect Effects 0.000 claims abstract description 5
- 238000013461 design Methods 0.000 claims abstract description 4
- 230000002087 whitening effect Effects 0.000 claims description 23
- 230000003595 spectral effect Effects 0.000 claims description 14
- 238000001228 spectrum Methods 0.000 claims description 8
- 238000011160 research Methods 0.000 description 7
- 238000012545 processing Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 2
- 241000287196 Asthenes Species 0.000 description 1
- 238000012935 Averaging Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000005303 weighing Methods 0.000 description 1
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
本发明涉及地震解释性处理领域，针对地震资料分辨率低、同相轴连续性差、能量不均衡的问题，本发明提出了一种地震信号处理方法，该方法根据原始的单道地震信号设计出白化滤波器，以希尔伯特‑黄变换为基础，将原始的地震信号通过EMD得到IMF分量，使用白化滤波器分别对每个IMF分量进行振幅增强处理，获得时间域IMF分量信号，对时间域IMF分量信号进行叠加，然后进行多道均衡处理获得最终的地震信号。本发明提出的方法，能够保持EMD算法的分解效率及IMF的精度，能够在合理提高分辨率的同时，保持较好的同相轴连续性和信噪比。
Description
技术领域
本发明涉及地震解释性处理领域，具体涉及一种地震数据处理方法，该方法以希尔伯特-黄变换的算法为基础，得到最终的地震信号。
背景技术
针对地震资料分辨率低的问题，很多学者进行了深入的研究和讨论。常规的方法对于平稳信号有较好的效果，但是对于非平稳信号(如地震信号)而言，常规方法不能满足实际需求。希尔伯特-黄变换(HHT)是一种适用于非平稳、非线性数据的分析研究方法，已在诸多领域如信号处理、地震工程学、生物医学等方面取得了广泛应用。本发明以HHT为技术手段，通过经验模态分解(EMD)得到多个本征模式函数(IMF)，然后分别对单个IMF分量进行提频处理，最终对提频的IMF分量进行叠加，得到提频后的地震信号，达到提高地震资料分辨率的目的。目前关于这方面研究的文献极少，实际应用过程中存在同相轴连续性差、能量不均衡等问题，制约了该方法的有效应用，所以需要进行进一步探索研究。
发明内容
针对现有技术中存在的地震资料分辨率低、同相轴连续性差、能量不均衡问题，本发明提出了一种地震信号处理方法。
本发明提出了一种地震信号处理方法，该方法以希尔伯特-黄变换为基础，将原始的地震信号通过经验模态分解得到本征模式函数分量，分别对每个本征模式函数分量进行振幅增强处理，获得时域信号，对时域信号进行叠加，然后进行多道均衡处理获得最终的地震信号。
通过实际资料的应用，本发明的提出的地震信号处理方法能较好地提高信号的分辨率，同时也具有较好的同相轴连续性和信噪比。
作为对本发明的进一步改进，本发明提出的方法包括如下步骤：
步骤一：获得白化滤波器表达式，如下
其中，为外包络振幅值，v是外包络的最大值即ε为白噪因子，ω表示频率。
作为对步骤一的进一步改进，对于同一道地震数据来说，利用傅里叶变换获得原始地震信号的振幅谱X(ω)，通过求取振幅谱X(ω)的绝对值|X(ω)|，获得所述外包络振幅函数从而获得所述白化滤波器表达式，同一道地震数据，采用统一的白化滤波器f(ω)。
进一步，通过试验表明，当白噪因子ε取值为0.1、0.5、0.8中的一个，尤其当白噪因子ε取值为0.5时，既能有效地提高地震信号的分辨率，也不至于由于分辨率过高造成地震信号的失真。
步骤二：对第j个本征模式函数分量的瞬时振幅进行振幅增强，获得白化后的瞬时振幅，所述白化后的瞬时振幅表达式如下
其中，aj(t)为第j个本征模式函数分量的瞬时振幅，ωj(t)为第j个本征模式函数分量的瞬时频率，j＝1,2,3…n。
作为对步骤二的进一步改进，将原始地震信号x(t)通过经验模态分解获得多个本征模式函数分量，对于其中第j个本征模式函数分量cj(t)来说，其瞬时振幅和瞬时频率分别为aj(t)和ωj(t)，为了进行振幅增强，对于每个时间点的aj(t)，使其与白化滤波器相乘进行振幅增强。
步骤三：获得重构后时间域本征模式函数分量，表达式如下
其中，e为自然常数2.71828，i为虚数符号，为第j个本征模式函数分量谱白化后的瞬时振幅，θj(t)是对应的瞬时相位。
作为对步骤三的进一步改进，在获得所述时域信号表达式过程中，地震信号的相位信息保持不变，利用希尔伯特变换获得信号的瞬时相位，利用谱白化后的瞬时振幅和瞬时相位相结合的方式获得所述时域信号表达式。
步骤四：获得谱白化后的地震信号，表达式为
如果获得的地震信号y(t)存在较为明显的道间能量不一致的问题，可通过多道均衡处理表达式对地震信号y(t)进行处理，可以提高同相轴横向连续性和能量一致性。
步骤五：对谱白化后的单道地震信号进行多道均衡处理，获得最终的地震信号，多道均衡处理表达式为
总之，本发明提出的地震信号处理方法，能够保持经验模态分解算法的分解效率及本征模式函数的精度，而且在合理提高地震信号分辨率的同时，使地震信号保持了较好的同相轴连续性和信噪比。
附图说明
在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：
图1为本发明提出的地震信号处理方法流程图；
图2为单道地震信号；
图3为图2中单道地震信号振幅谱绝对值及外包络；
图4为白化滤波器曲线；
图5为原始地震信号与经过本发明方法处理后地震信号剖面对比图；
在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例。
具体实施方式
在研究区开展了基于HHT的高分辨率地震数据处理方法研究和应用。数据处理是单道进行的，单道地震信号如图2所示。首先从图2中获得地震信号的振幅谱X(ω)，获得振幅谱X(ω)的绝对值|X(ω)|，得到振幅谱外包络曲线如图3所示，设外包络曲线函数为然后通过得到的振幅谱包络曲线，设计白化滤波器，得到白化滤波器表达式如下：
其中，为外包络振幅值，v是外包络的最大值即ε为白噪因子，ω表示频率。
白化滤波器的设计根据白化滤波器表达式计算得到，其中变量参数为白噪因子ε。白噪因子ε的取值直接影响了白化滤波器的设计。为了研究分析不同白噪因子取值情况下的白化滤波器形态，图4显示了ε取值分别为0.1、0.5和0.8时对应的白化滤波器曲线。前人研究认为，参数ε一般选用0.5，这既能有效地提高分辨率，也不至于分辨率过高造成地震信号的失真，因此选用ε＝0.5的白化滤波器，得到ε＝0.5的白化滤波器曲线，如图4。
在确定了白化滤波器之后，将单道地震信号x(t)通过EMD分解，得到n个IMF分量。然后利用ε＝0.5的白化滤波器对每个IMF分量的瞬时振幅进行能量增强，即利用谱白化后的瞬时振幅表达式对瞬时振幅能量进行增强，谱白化后的瞬时振幅表达式如下：
其中，aj(t)为第j个本征模式函数分量的瞬时振幅，ωj(t)为第j个本征模式函数分量的瞬时频率，j＝1,2,3…n。
由于瞬时振幅增强前后数据的瞬时相位保持不变，利用希尔伯特变换求得信号的瞬时相位，通过时间域本征模式函数分量表达式将谱白化后的瞬时振幅转化为时间域IMF信号，重构后时间域本征模式函数分量表达式如下：
其中，e为自然常数2.71828，i为虚数符号，为第j个本征模式函数分量谱白化后的瞬时振幅，θj(t)是对应的瞬时相位。
最后通过谱白化后的地震信号表达式将时间域IMF信号叠加，得到处理后的单道信号y(t)，谱白化后的地震信号，表达式为
由于得到的剖面存在较强的能量不一致性，所以利用多道均衡处理表达式进行道间均衡处理，多道均衡处理表达式为
进行多道均衡处理后，使得数据在具有较高分辨率和信噪比的同时，也具有很好的同相轴连续性和能量一致性。
处理前的地震信号和使用本发明提出的方法处理后剖面对比如图5所示，见图中椭圆框内对比，可以看出使用本发明提出的方法处理后的地震信号分辨率明显提高，即通过本发明提出的方法对地震信号进行处理后，剖面分辨率明显提高，原先叠合的同相轴清晰地分成了两套同相轴，这将大大有利于薄层砂体的识别与分析。
最后说明的是，以上实施例仅用于说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。
Claims (8)
1.一种地震信号处理方法，所述方法以希尔伯特-黄变换为基础，将原始的地震信号通过经验模态分解得到本征模式函数分量，分别对每个本征模式函数分量进行振幅增强处理，获得时间域本征模式函数分量信号，对时间域本征模式函数分量信号进行叠加，然后进行多道均衡处理获得最终的地震信号。
2.根据权利要求1所述的地震信号处理方法，其特征在于，所述方法包括如下步骤：
2.1设计白化滤波器，白化滤波器表达式如下
其中，为外包络振幅值，v是外包络的最大值即ε为白噪因子，ω表示频率；
2.2对第j个本征模式函数分量的瞬时振幅进行振幅增强，获得白化后的瞬时振幅，所述白化后的瞬时振幅表达式如下
<mrow>
<msub>
<mover>
<mi>a</mi>
<mo>~</mo>
</mover>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>a</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&omega;</mi>
<mi>j</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
其中，aj(t)为第j个本征模式函数分量的瞬时振幅，ωj(t)为第j个本征模式函数分量的瞬时频率，j＝1,2,3…n；
2.3获得重构后时间域本征模式函数分量，表达式如下
<mrow>
<msub>
<mover>
<mi>s</mi>
<mo>~</mo>
</mover>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>Re</mi>
<mo>&lsqb;</mo>
<msub>
<mover>
<mi>a</mi>
<mo>~</mo>
</mover>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<msub>
<mi>i&theta;</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>&rsqb;</mo>
</mrow>
其中，e为自然常数2.71828，i为虚数符号，为第j个本征模式函数分量谱白化后的瞬时振幅，θj(t)是对应的瞬时相位；
2.4将时间域本征模式函数分量叠加，获得谱白化后的地震信号，表达式为
<mrow>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<msub>
<mover>
<mi>s</mi>
<mo>~</mo>
</mover>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
2.5对谱白化后的单道地震信号进行多道均衡处理，获得最终的地震信号，多道均衡处理表达式为
<mrow>
<msub>
<mi>z</mi>
<mi>m</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mi>m</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>m</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</munderover>
<msub>
<mi>y</mi>
<mi>j</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mn>3</mn>
<mo>,</mo>
<mi>m</mi>
<mo>=</mo>
<mn>2</mn>
<mo>,</mo>
<mn>3</mn>
<mo>,</mo>
<mn>4</mn>
<mo>...</mo>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
3.根据权利要求2所述的地震信号处理方法，其特征在于，在步骤2.1中，利用傅里叶变换获得原始地震信号的振幅谱X(ω)，通过求取振幅谱X(ω)的绝对值|X(ω)|，获得所述外包络振幅函数从而获得所述白化滤波器表达式。
4.根据权利要求3所述的地震信号处理方法，其特征在于，在步骤2.1中，所述白噪因子ε取值为0.1、0.5、0.8中的一个。
5.根据权利要求4所述的地震信号处理方法，其特征在于，所述白噪因子ε取值为0.5。
6.根据权利要求2所述的地震信号处理方法，其特征在于，在步骤2.2中，将原始地震信号x(t)通过经验模态分解获得所述本征模式函数分量。
7.根据权利要求2所述的地震信号处理方法，其特征在于，在步骤2.3中，利用谱白化后的瞬时振幅和瞬时相位相结合的方式获得所述时域信号表达式。
8.根据权利要求7所述的地震信号处理方法，其特征在于，在获得所述时域信号表达式过程中，地震信号的相位信息保持不变，利用希尔伯特变换获得信号的瞬时相位。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611015681 | 2016-11-18 | ||
| CN2016110156814 | 2016-11-18 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108072894A true CN108072894A (zh) | 2018-05-25 |
Family
ID=62161672
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201611071443.5A Pending CN108072894A (zh) | 2016-11-18 | 2016-11-29 | 一种地震信号处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108072894A (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104849756A (zh) * | 2015-03-31 | 2015-08-19 | 中国地质大学（北京） | 一种提高地震数据分辨率增强有效弱信号能量的方法 |
| CN105044769A (zh) * | 2015-06-10 | 2015-11-11 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 提高地震信号的分辨率的方法 |
| CN105116442A (zh) * | 2015-07-24 | 2015-12-02 | 长江大学 | 岩性油气藏弱反射地震信号的重构方法 |
-
2016
- 2016-11-29 CN CN201611071443.5A patent/CN108072894A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104849756A (zh) * | 2015-03-31 | 2015-08-19 | 中国地质大学（北京） | 一种提高地震数据分辨率增强有效弱信号能量的方法 |
| CN105044769A (zh) * | 2015-06-10 | 2015-11-11 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 提高地震信号的分辨率的方法 |
| CN105116442A (zh) * | 2015-07-24 | 2015-12-02 | 长江大学 | 岩性油气藏弱反射地震信号的重构方法 |
Non-Patent Citations (1)
| Title |
|---|
| 王季: "基于Hilbert谱白化的高分辨率地震资料处理", 《煤炭学报》 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yinfeng et al. | Analysis of earthquake ground motions using an improved Hilbert–Huang transform | |
| CN105760347A (zh) | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 | |
| Du et al. | Source separation of diesel engine vibration based on the empirical mode decomposition and independent component analysis | |
| Zheng et al. | Extreme-point weighted mode decomposition | |
| CN103699513A (zh) | 一种基于多尺度噪声调节的随机共振方法 | |
| CN106886044A (zh) | 一种基于剪切波与Akaike信息准则的微地震初至拾取方法 | |
| Rui et al. | Online wavelet denoising via a moving window | |
| CN102988041B (zh) | 心磁信号噪声抑制中的信号选择性平均方法 | |
| CN104408288A (zh) | 基于小波和参数补偿的多稳态随机共振微弱信号检测方法 | |
| CN105044769B (zh) | 提高地震信号的分辨率的方法 | |
| Zhou et al. | Geomagnetic sensor noise reduction for improving calibration compensation accuracy based on improved HHT algorithm | |
| Zhidong et al. | A new method for processing end effect in empirical mode decomposition | |
| Li et al. | Magnetotelluric signal-noise separation method based on SVM–CEEMDWT | |
| CN108072894A (zh) | 一种地震信号处理方法 | |
| Wang et al. | Desert seismic noise suppression based on multimodal residual convolutional neural network | |
| CN107831535A (zh) | 时变相位分解与重构方法 | |
| Liu et al. | Speech enhancement based on Hilbert-Huang transform | |
| Wang et al. | Multicomponent seismic noise attenuation with multivariate order statistic filters | |
| CN106125148B (zh) | 一种针对有源周期电磁信号的降噪方法及装置 | |
| Sun et al. | Application of adaptive iterative low-rank algorithm based on transform domain in desert seismic signal analysis | |
| CN103577877A (zh) | 一种基于时频分析和bp神经网络的船舶运动预报方法 | |
| Xu et al. | Analysis of the time-frequency characteristics of internal combustion engine vibration signal based on Hilbert-Huang transform | |
| Yang et al. | Robustness analysis of synchrosqueezed transforms | |
| CN107590788A (zh) | 一种遥感图像处理方法 | |
| CN106443771A (zh) | 提高转换波地震数据分辨率的方法及其速度反演方法 |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20180525 |