CN105093325B - 一种定量的提频方法 - Google Patents
一种定量的提频方法 Download PDFInfo
- Publication number
- CN105093325B CN105093325B CN201410218603.9A CN201410218603A CN105093325B CN 105093325 B CN105093325 B CN 105093325B CN 201410218603 A CN201410218603 A CN 201410218603A CN 105093325 B CN105093325 B CN 105093325B
- Authority
- CN
- China
- Prior art keywords
- mrow
- frequency
- analytical function
- mover
- function
- 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
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
Abstract
本发明提供了一种定量的提频方法，属于地震资料信号处理领域。本方法包括：(1)利用原始信号道集f(t)与原始信号道集f(t)的希尔伯特变换构建解析函数Z1(t)；(2)在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数为Z2(t)；(3)由解析函数Z1(t)、Z2(t)得到新的解析函数Z(t)；(4)取新的解析函数Z(t)的实部作为提高频率后信号的输出；(5)在给定的频率范围(fmin，fmax)Hz内，Δf从最小频率值开始，以给定的步长循环迭代，并将每个频率的信号输出并累加，直到最大频率值为止即完成定量提频过程。
Description
技术领域
本发明属于地震资料信号处理领域，具体涉及一种定量的提频方法。
背景技术
频谱分析是现代信号处理与分析中的一个重要手段，将信号源发出的信号强度按频率顺序展开，使其成为频率的函数，并考察变化规律，称为频谱分析，它是将时域信号变换至频域，其目的是把复杂的时间历程波形经过傅里叶变换分解为若干单一的谐波分量来研究，从而获得信号的频率结构以及各谐波和相位信息。若信号频率过低时，地震资料的分辨率则较低，可以通过一定的数字信号处理手段对信号进行提频处理，使得信号的主频变高及频带变宽，以此来提高分辨率。
现有提频技术主要有以下几种方法：
(1)谱白化法。谱白化是一种展宽信号振幅谱的方法，来达到补偿频率衰减的目的。谱白化方法主要存在以下问题：分频后各个频带数据随时间增益无法满足时间方向保持谱白化前的相对振幅和波形关系；分频后的各个频带数据增益后再叠加时，无法确定各个频率振幅的比例成分；
(2)反Q滤波。反Q滤波是地震资料处理的常规处理技术，反Q滤波的目的是消除地震波在传播过程中的振幅衰减和相位畸变，使信号高频成分得到恢复，从而提高地震资料的分辨率。但是反Q滤波的准确程度与品质因子Q值求取的正确性密切相关，由于地下构造复杂，影响地震波衰减的因素非常多，Q值往往难以准确求取，从而导致反Q滤波不准确。
发明内容
本发明的目的在于解决上述现有技术中存在的难题，提供一种定量的提频方法，以解析函数与三角函数为基础，实现对信号定量频率的提高，从而提高分辨率。
本发明是通过以下技术方案实现的：
一种定量的提频方法，包括：
(1)利用原始信号道集f(t)与原始信号道集f(t)的希尔伯特变换构建解析函数Z1(t)；
(2)在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数为Z2(t)；
(3)由解析函数Z1(t)、Z2(t)得到新的解析函数Z(t)；
(4)取新的解析函数Z(t)的实部作为提高频率后信号的输出；
(5)在给定的频率范围(fmin，fmax)Hz内，Δf从最小频率值开始，以给定的步长循环迭代，并将每个频率的信号输出并累加，直到最大频率值为止即完成定量提频过程。
所述步骤(1)是这样实现的：
希尔伯特变换表达式为：
利用原始信号道集f(t)与原始信号道集f(t)的希尔伯特变换构建解析函数Z1(t)为：
其中f(t)为解析函数Z1(t)的实部，为解析函数Z1(t)的虚部。
所述步骤(2)是这样实现的：
在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数Z2(t)为：
Z2(t)＝cos(2πΔft)+jsin(2πΔft) (4)
其中cos(2πΔft)为解析函数Z2(t)的实部，sin(2πΔft)为解析函数Z2(t)的虚部。
所述步骤(3)是这样实现的：
由解析函数Z1(t)、Z2(t)得到新的解析函数：
Z(t)＝Re(Z1)Re(Z2)-Im(Z1)Im(Z2)+j{Re(Z1)Im(Z2)+Im(Z1)Re(Z2)}
代入Z1(t)、Z2(t)的实部和虚部得到解析函数Z(t)的表达式为：
本方法的优点是可以进行有方向性(频率减小方向、频率增大方向、频率减小和增大双方向)的频率拓展，至于三个方向中是进行的哪一个方向的处理，取决于给定的频率范围(fmin，fmax)的取值，如(-10，0)是频率减小方向，(0，10)是频率增大方向，(-10，10)是频率减小和增大两个方向同时。因此，一旦给定了频率范围(fmin，fmax)就相当于明确了频率拓展的方向，完成第(4)步后，只是完成了(fmin，fmax)范围内的一个频率Δf，整个方法的完成必须从最小频率值开始，以给定的步长循环迭代，直到最大频率值为止，即若步长为Δk时，则Δf＝fmin+n*Δk，n＝0，1，...。
与现有技术相比，本发明的有益效果是：本发明所描述的可以对信号频谱进行定量的频率拓展，并且可以向频率增大与频率减小中任意方向或同时两个方向的拓展，使得频谱在一定可控的范围内得到拓展来提高频率，从而提高了分辨率。
附图说明
图1是本发明方法的步骤框图。
图2a-1是未做提频处理的雷克子波波形图。
图2a-2是对应图2a-1的频谱。
图2b-1是fmin＝-5Hz，fmax＝0Hz处理后的子波。
图2b-2是对应图2b-1的频谱。
图3a-1是未做提频处理子波。
图3a-2是对应图3a-1的频谱。
图3b-1是fmin＝0Hz，fmax＝18Hz处理后的子波。
图3b-2是对应图3b-1的频谱。
图4a-1是未做提频处理的子波。
图4a-2是对应图4a-1的频谱。
图4b-1是fmin＝-5Hz，fmax＝18Hz处理后的子波。
图4b-2是对应图4b-1的频谱。
图5a是含三层水平同相轴模型及叠加，上到下同相轴主频依次为40Hz，30Hz，20Hz。
图5b是fmin＝-10Hz，fmax＝20Hz提频处理后的模型及叠加。
图6a是图5a中某单道频谱。
图6b是fmin＝-10Hz，fmax＝20Hz提频处理后同一单道频谱。
具体实施方式
下面结合附图对本发明作进一步详细描述：
本发明为提高分辨率提供了一种定量的提频方法。具体实现步骤是：
众所周知，通过傅里叶正反变换可以实现信号时间域和频率域的转变，傅氏正反变换如公式(1)
在频率域，傅里叶变换的典型用途是将信号分解成幅值谱-显示与频率对应的幅值大小，横轴一般为频率，纵轴为振幅。
当原始地震信号低频化时，在频率域表现为子波主频小、频带较窄，因而具有较低的分辨率，将不利于信号的分析和进一步的地震解释工作。
为了能够定量提高地震剖面分辨率，利用原始信号及其希尔伯特变换结果构建解析函数，即首先对信号进行希尔伯特变换，并与原始输入信号构建解析函数，然后以需要提高的频率值的三角函数分别构建实部和虚部，得到另一解析函数，最后利用以上两个解析函数得到频率提高后信号的解析函数，并取其实部作为输出信号。具体见流程图，如图1所示。
设实函数道f(t)的希尔伯特变换为则希尔伯特变换表达式为：
由原始信号构建的解析函数Z1(t)为：
其中f(t)为解析函数的实部即为原始输入信号，为解析函数的虚部即为输入信号f(t)的希尔伯特变换结果。
在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数Z2(t)为：
Z2(t)＝cos(2πΔft)+jsin(2πΔft) (4)
其中cos(2πΔft)为解析函数的实部，sin(2πΔft)为解析函数的虚部。
最后，由解析函数Z1(t)、Z2(t)得到新的解析函数表达式为：
Z(t)＝Re(Z1)Re(Z2)-Im(Z1)Im(Z2)+j{Re(Z1)Im(Z2)+Im(Z1)Re(Z2)}
代入Z1(t)、Z2(t)的实部和虚部得到解析函数Z(t)的表达式为：
取Z(t)的实部作为提高频率后信号的输出。
程序设计中，从最小频率值到最大频率值循环迭代，(Δf为(fmin，fmax)Hz内从小到大的取值，若步长为Δk，则Δf＝fmin+n*Δk，n＝0，1，...)，实现频率增大(起始值fmin＝0，最大值为fmax)或频率减小(起始值为fmin，最大值fmax＝0)或同时两个方向的拓展(起始值为fmin，最大值为fmax)。
为了验证本方法的正确性，做如下实验。图2a-1是雷克子波波形图和其对应的频谱(图2a-2)，图2b-1和图2b-2是提频范围为(-5，0)Hz提频处理后的子波波形图和频谱图对比结果；图3a-1和图3a-2、图3b-1和图3b-2和图4a-1和图4a-2，图4b-1和图4b-2则分别是提频范围为(0，18)Hz、(-5，18)Hz处理后的子波波形图和频谱图对比结果。从以上对比图可以看到，取不同的提频范围即可实现不同方向的定量频率的拓展。
图5a为含有三层水平同相轴模型及叠加，同相轴从上至下主频依次为40Hz，30Hz，20Hz；图5b为其提频范围为(-10，20)Hz处理后的同相轴模型及叠加，图6a、图6b为图5a、图5b中某一道对应的频谱，可知同相轴变“细”的同时，频带变宽，分辨率提高。
上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
Claims (4)
1.一种定量的提频方法，其特征在于：所述方法包括：
(1)利用原始信号道集f(t)与原始信号道集f(t)的希尔伯特变换构建解析函数Z1(t)；
(2)在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数为Z2(t)；
(3)通过以下公式由解析函数Z1(t)、Z2(t)得到新的解析函数Z(t)：
Z(t)＝Re(Z1)Re(Z2)-Im(Z1)Im(Z2)+j{Re(Z1)Im(Z2)+Im(Z1)Re(Z2)}；
(4)取新的解析函数Z(t)的实部作为提高频率后信号的输出；
(5)在给定的频率范围(fmin，fmax)Hz内，Δf从最小频率值开始，以给定的步长循环迭代，并将每个频率的信号输出并累加，直到最大频率值为止即完成定量提频过程。
2.根据权利要求1所述的定量的提频方法，其特征在于：所述步骤(1)是这样实现的：
希尔伯特变换表达式为：
<mrow>
<mover>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>^</mo>
</mover>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>&pi;</mi>
</mfrac>
<msubsup>
<mo>&Integral;</mo>
<mrow>
<mo>-</mo>
<mi>&infin;</mi>
</mrow>
<mrow>
<mo>+</mo>
<mi>&infin;</mi>
</mrow>
</msubsup>
<mfrac>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>t</mi>
<mo>-</mo>
<mi>&tau;</mi>
</mrow>
</mfrac>
<mi>d</mi>
<mi>&tau;</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
利用原始信号道集f(t)与原始信号道集f(t)的希尔伯特变换构建解析函数Z1(t)为：
<mrow>
<msub>
<mi>Z</mi>
<mn>1</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>j</mi>
<mover>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>^</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
其中f(t)为解析函数Z1(t)的实部，为解析函数Z1(t)的虚部。
3.根据权利要求1所述的定量的提频方法，其特征在于：所述步骤(2)是这样实现的：
在提频范围(fmin，fmax)Hz内，假设定量提高的频率值为Δf，利用三角函数构建的另一解析函数Z2(t)为：
Z2(t)＝cos(2πΔft)+jsin(2πΔft) (4)
其中cos(2πΔft)为解析函数Z2(t)的实部，sin(2πΔft)为解析函数Z2(t)的虚部。
4.根据权利要求1所述的定量的提频方法，其特征在于：所述步骤(3)包括：
代入Z1(t)、Z2(t)的实部和虚部得到解析函数Z(t)的表达式为：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>Z</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>{</mo>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
<mi>&Delta;</mi>
<mi>f</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>^</mo>
</mover>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
<mi>&Delta;</mi>
<mi>f</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
<mo>}</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mi>j</mi>
<mrow>
<mo>{</mo>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
<mi>&Delta;</mi>
<mi>f</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mover>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>^</mo>
</mover>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
<mi>&Delta;</mi>
<mi>f</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
<mo>}</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
1
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410218603.9A CN105093325B (zh) | 2014-05-22 | 2014-05-22 | 一种定量的提频方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410218603.9A CN105093325B (zh) | 2014-05-22 | 2014-05-22 | 一种定量的提频方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105093325A CN105093325A (zh) | 2015-11-25 |
| CN105093325B true CN105093325B (zh) | 2017-08-18 |
Family
ID=54574179
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410218603.9A Active CN105093325B (zh) | 2014-05-22 | 2014-05-22 | 一种定量的提频方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105093325B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105700014B (zh) * | 2016-01-26 | 2018-05-15 | 电子科技大学 | 一种基于频域显著性检测的地震属性分析方法 |
| CN112987092A (zh) * | 2021-02-05 | 2021-06-18 | 中国石油天然气股份有限公司 | 一种台盆区针对缝洞储层的地震资料预处理方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5850622A (en) * | 1996-11-08 | 1998-12-15 | Amoco Corporation | Time-frequency processing and analysis of seismic data using very short-time fourier transforms |
| CN102053276A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种地震数字信号的复数道集二维滤波方法 |
| CN102466819A (zh) * | 2010-11-03 | 2012-05-23 | 中国石油天然气集团公司 | 地震信号的频谱分析方法及装置 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6963804B2 (en) * | 2002-11-19 | 2005-11-08 | Saudi Arabian Oil Company | Seismic data processing method to enhance fault and channel display |
-
2014
- 2014-05-22 CN CN201410218603.9A patent/CN105093325B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5850622A (en) * | 1996-11-08 | 1998-12-15 | Amoco Corporation | Time-frequency processing and analysis of seismic data using very short-time fourier transforms |
| CN102053276A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种地震数字信号的复数道集二维滤波方法 |
| CN102466819A (zh) * | 2010-11-03 | 2012-05-23 | 中国石油天然气集团公司 | 地震信号的频谱分析方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105093325A (zh) | 2015-11-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102116868B (zh) | 一种地震波分解方法 | |
| CN105760347A (zh) | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 | |
| CN103412329B (zh) | 一种提高地震数据分辨率的方法 | |
| CN103941091A (zh) | 基于改进emd端点效应的电力系统hht谐波检测方法 | |
| CN103837884B (zh) | 基于时域分析的数字核脉冲信号梯形成形算法 | |
| CN105093325B (zh) | 一种定量的提频方法 | |
| CN107480391A (zh) | 基于数据驱动的近断层非平稳地震动模拟方法 | |
| GB2447235A (en) | Method for estimating absorption parameter Q(T) | |
| CN112084732B (zh) | 一种基于fpga的谐波补偿方法 | |
| CN103927761A (zh) | 一种基于稀疏表示的故障微弱信号特征提取方法 | |
| CN105277973A (zh) | 一种基于匹配追踪的子波分解优化方法 | |
| CN105044769B (zh) | 提高地震信号的分辨率的方法 | |
| CN107219551A (zh) | 拓宽地震数据频带的方法及装置 | |
| Nam et al. | A super-resolution spectrogram using coupled PLCA. | |
| CN103543331B (zh) | 一种计算电信号谐波和间谐波的方法 | |
| CN103604404B (zh) | 基于数值积分的加速度信号测取位移方法 | |
| CN106249282A (zh) | 一种适用于avaf反演的频率域地震道集生成方法 | |
| CN106483555B (zh) | Enpemf数据的格林函数-spwvd时频分析方法 | |
| CN108806721A (zh) | 信号处理器 | |
| CN102543091B (zh) | 一种模拟音效的生成系统及方法 | |
| CN103941280A (zh) | 基于冲激响应不变法的数字核脉冲高斯成形方法 | |
| CN104467836A (zh) | 一种时钟信号发生方法及系统 | |
| CN104156509A (zh) | 一种噪声合成方法 | |
| CN108507743A (zh) | 一种基于lmd的冲击响应谱时域信号合成的方法及装置 | |
| CN103344988A (zh) | 基于k-l分解的可控震源信号相位检测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |