CN106154318A - 一种地震频谱重建方法 - Google Patents
一种地震频谱重建方法 Download PDFInfo
- Publication number
- CN106154318A CN106154318A CN201510182064.2A CN201510182064A CN106154318A CN 106154318 A CN106154318 A CN 106154318A CN 201510182064 A CN201510182064 A CN 201510182064A CN 106154318 A CN106154318 A CN 106154318A
- Authority
- CN
- China
- Prior art keywords
- spectrum
- amplitude
- wavelet
- frequency
- formula
- 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.)
- Granted
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种地震频谱重建方法,涉及地球物理勘探领域。本方法包括:(1)获取地震信号x(t),对所述地震信号x(t)进行傅里叶变换,获得该地震信号的振幅谱X(f);(2)选取子波频谱公式,拾取子波频谱的峰值频率m的初始值;(3)计算子波频谱的峰值振幅a;(4)根据最大匹配投影原理对子波频谱的峰值振幅a和子波频谱的峰值频率m进行局部优化得到局部最优的子波频谱的峰值振幅和子波频谱的峰值频率;(5)将子波频谱的峰值振幅a和子波频谱的峰值频率m带入子波频谱公式获得子波的振幅谱;利用局部最优的子波频谱的峰值振幅a和子波频谱的峰值频率m对地震信号的振幅谱X(f)进行重建,获得重建后的地震频谱以及残差。
Description
技术领域
本发明涉及地球物理勘探领域,具体涉及一种地震频谱重建方法。
背景技术
地震频谱成像技术已成为一项储层特色解释技术,已在地震属性分析方面发挥了重要的作用。频谱分解的核心是信号的时频分析方法,时频分析方法一般分为线性时频分析方法、双线性时频分析方法和参数化时频分析方法等,所有的时频分析方法都有其自身固有的优点和缺陷。利用时频分析方法获得的频率道集,在每一个频率上计算的所谓的“单频信息”都是一定频带范围内信息的叠合,而且反应的是所有地震响应在局部频带范围内信息的叠合;同时,在每一个时间深度上,瞬时谱展示的是整个频率范围内所有的频率成分,没有提供给用户感兴趣的成分,不能获得每一个地震响应所对应的频率成分。
分频属性能够更好地刻画储集层,但是解释人员往往趋向于通过带通滤波拾取一定频带的成分。然而地震子波的主瓣的震荡是有限的,转换到频率域后它对应着一个特定频带范围,解释人员很难从频谱中拾取某一个地震响应的频率范围。如果能够获得每一个频率对应的地震响应,那么就能够对感兴趣的特定储层进行更细致的研究。2010年,Tomasso提出了“频谱重构”的概念,通过正演模拟来估计地震信号的频率成分。不同于传统的频谱分解方法,频谱重构是利用模型的频谱来重建信号的频谱,从地震信号中提取出有意义的分量,而不是分解该信号。频谱重构能够只显示某一个特定的频谱成分,这个特定的成分与某一个特定的层相联系,而这个特定的层可能就是用户感兴趣的目的层。然而,Tomasso采取手工拾取频谱成分和振幅,这种做法不准确,而且严重的依赖于个人的经验,而且针对大数据体应用时也受限制。
发明内容
本发明的目的在于解决上述现有技术中存在的难题,提供一种地震频谱重建方法,采用最大匹配投影原理自动对频谱进行重建,重建过程不包含方程求解,稳定性较好,且重建结果不受人为因素影响。
本发明是通过以下技术方案实现的:
一种地震频谱重建方法,包括:
(1)获取地震信号x(t),对所述地震信号x(t)进行傅里叶变换,获得该地震信号的振幅谱X(f);
(2)选取子波频谱公式,拾取子波频谱的峰值频率m的初始值;
(3)计算子波频谱的峰值振幅a;
(4)根据最大匹配投影原理对子波频谱的峰值振幅a和子波频谱的峰值频率m进行局部优化得到局部最优的子波频谱的峰值振幅和子波频谱的峰值频率;
(5)将子波频谱的峰值振幅a和子波频谱的峰值频率m带入子波频谱公式获得子波的振幅谱;利用局部最优的子波频谱的峰值振幅a和子波频谱的峰值频率m对地震信号的振幅谱X(f)进行重建,获得重建后的地震频谱以及残差;
(6)判断所述残差是否满足给定条件,如果是,则转入步骤(7),如果否,则返回步骤(2);
(7)结束。
所述步骤(2)中子波频谱公式如下:
公式(1),
其中a为该子波频谱的峰值振幅、m为该频谱的峰值频率,f为频率。
所述步骤(3)是利用下式计算子波频谱的峰值振幅a:
公式(3)
其中为子波振幅谱集合,<riX,Ri>表示riX和Ri的内积,为第i次重建振幅谱的标准化;如果i=1,ri-1X=r0X即初始振幅谱X;如果i>1,则ri-1X为第i次重建后的残差;argmax表示取最大投影。
所述步骤(4)是这样实现的:
根据公式(3)估算参数ai:给定参数ai的范围,在该范围内取多值作为ai,将各值与峰值频率m的初始值mi一起带入公式(1)获得多个子波的振幅谱,这些子波的振幅谱集合即为;在该子波振幅谱集合中确定满足公式(3)的一个子波振幅谱Ri(mi,f),该振幅谱对应的振幅即为估算的参数ai;其中参数ai的范围为n*Peak~Peak,其中n∈[0.1,0.9],Peak为riX的峰值振幅的值;
在获得初始的估算的参数ai和mi后,在局部范围内搜索参数集ξ,局部范围可表示为[ξ-Δξ,ξ+Δξ],其中Δa为振幅间隔,Δm为峰值频率间隔;在局部范围内,根据公式(3)确定最优的ai和mi,获得最优的参数集ξi,具体实现方法如下:在局部范围内取多个值作为ai和mi,形成多组重建参数ξi={ai,mi},将每一组参数ξi带入公式(1)获得子波的振幅谱,然后寻找一组满足公式(3)参数ξi,即为第i次重建的参数ξi={ai,mi},即得到局部最优的子波频谱的峰值振幅和子波频谱的峰值频率。
所述步骤(5)中的利用局部最优的子波频谱的峰值振幅和子波频谱的峰值频率对地震信号的振幅谱X(f)进行重建,获得残差是这样实现的:
重建后的振幅谱的形式为公式(2)所定义的形式:
公式(2),
其中n为重建次数、ai为与第i次重建的频谱Ri对应的局部最优的子波频谱的峰值振幅、mi为与第i次重建的频谱Ri对应的局部最优的子波频谱的峰值频率,rnX为第n次重建后的残差,且r0X=X,标记i为1到n的正整数, 为第i次重建的子波振幅谱,对于每一次重建,都有
所述步骤(6)中的给定条件为:
所述残差rn+1X的能量应小于等于原始振幅谱能量r0X的5%,r0X=X。
与现有技术相比,本发明提供的地震频谱重建方法采用最大匹配投影原理自动对频谱进行重建,子波频谱的振幅和峰值频率可自适应的调节,获得的重建频谱的和更能接近原始频谱;并且重建过程不包含方程求解,稳定性更好;重建过程不受人为因素影响。
附图说明
图1是根据本发明的一个实施方式的步骤框图。
具体实施方式
下面结合附图对本发明作进一步详细描述。
如图1所示,本发明地震频谱重建方法包括以下步骤:
(1)选取子波频谱形式,该子波频谱形式是由公式(1)定义的Ricker子波振幅谱形式:
公式(1),
其中a为该子波频谱的峰值振幅、m为该频谱的峰值频率,f为频率(公式R(m,f)是计算在频率为f时的振幅值),因此Ricker子波的振幅谱由参数集ξ={a,m}(如果选取不同的子波,将采用不同的参数集,Ricker 子波只是由{a,m}两个参数确定,如果采用其他的子波,将使用该子波所对应的几个参数来形成参数集即可)来描述;选取子波振幅谱形式可以从本领域相关的基础知识中获取,不限于Ricker子波振幅谱,比如常用的Morlet子波振幅谱也可作为子波振幅谱形式,如果采用其他的地震子波,那么子波的振幅谱形式也会变化,不过计算的整个思路是一样的。对于其他子波振幅谱形式就是其对应的公式。
(2)获取地震信号x(t),对所述地震信号x(t)进行傅里叶变换,获得该地震信号的振幅谱X(f);
对地震信号进行傅里叶变换,求取振幅谱都是公知的常用方法,其变换和计算方法在公知教科书中都有详细说明。
(3)对振幅谱X(f)进行多次重建,直到满足重建退出条件(从第一次重建开始,如果第一次重建的结果即满足了退出条件,那么就退出,如果不满足条件,则迭代进行第二次重建,依此类推,直到满足重建条件即停止计算。),重建后的振幅谱的形式为公式(2)所定义的形式:
公式(2),
其中n为重建次数、ai为与第i次重建的频谱Ri对应的振幅、mi为与第i次重建的频谱Ri对应的峰值频率,rnX为第n次重建后的残差,且r0X=X,标记i为1到n的正整数,为第i次重建的子波振幅谱,对于每一次重建,都有
根据最大匹配投影原理计算第i次重建时的参数集ξi={ai,mi};计算第i次重建参数ξi={ai,mi}时包括:拾取振幅谱ri-1X的峰值振幅对应的频率作为初始mi,根据公式(3)估算参数ai:
公式(3)
其中为子波振幅谱集合,<riX,Ri>表示riX和Ri的内积,为第i次重建振幅谱的标准化;如果i=1,ri-1X=r0X即初始振幅谱X;如果i>1,则ri-1X为第i次重建后的残差;argmax表示取最大投影。
具体的,根据公式(3)估算参数ai包括:给定参数ai的范围,在该范围内取多值作为ai,将各值与初始mi一起带入公式(1)获得多个子波的振幅谱,这些小波的振幅谱集合即为;在该子波振幅谱集中确定满足公式(3)的一个子波振幅谱Ri(mi,f),该振幅谱对应的振幅即为估算的参数ai。其中参数ai的范围为n*Peak~Peak,其中n∈[0.1,0.9],Peak为riX的峰值振幅的值。
为了更精确的得到第i次的重建参数ξi={ai,mi},在获得初始的估算的参数ai和mi后,还应在局部范围内搜索参数集ξ,局部范围可表示为[ξ-Δξ,ξ+Δξ],其中Δa为振幅间隔,Δm为峰值频率间隔;在局部范围内,根据公式(3)确定最优的ai和mi,获得最优的参数集ξi。具体实现方法为:在局部范围内取多个值作为ai和mi,形成多组重建参数ξi={ai,mi},将每一组参数ξi带入公式(1)获得子波的振幅谱,然后寻找一组满足公式(3)参数ξi,即为第i次重建的参数ξi={ai,mi}(即最优的参数集ξi,在范围[ξ-Δξ,ξ+Δξ]内生成多个参数ξi={ai,mi},通过公式(3)判断最好的一组作为搜寻结果)。
其所述重建退出条件可以为所述残差rn+1X的能量应小于原始振幅谱能量(即r0X=X)的5%,甚至更小,具体应视实际情况而定。每一次迭代的信号振幅谱的能量即为按频率对各能量(表达式就是公式(1),公式(1)表示的是一个频率f下的能量,信号振幅谱的能量包含了多个频率,比如如果信号采样间隔为2ms,那么最大频率就是250HZ,同时如果频率域采样间隔为1HZ,那么就有251个公式(1),公式(1)中的f就是0-250,将计算得到的251个值相加,再乘以频率域采样间隔,即得到信号振幅谱的能量。)进行积分,即将每个频率下的振幅相加乘以频率域采样间隔即可。
(4)将所述计算出的参数ai和mi(就是步骤(3)最后得到的重建的参数 ξi={ai,mi})带入公式(1)得到第i次重建出的雷克子波振幅谱。
本发明提供的地震频谱重建方法是多次迭代的过程。对于第一次迭代重建,有即此时将X(f)作为初始振幅谱,根据上述方法确定ξ1={a1,m1},由此得到第一次迭代重建的子波振幅谱R1(m1,f)。由于则有在第二次迭代中,将r1X作为初始振幅谱,根据上述方法确定ξ2={a2,m2},由此获得第二次重建的子波振幅谱R2(m2,f)。在第三次迭代中,将r2X作为初始振幅谱。以此类推,在经过N次迭代之后就能得到第N次的残差rnX,当残差满足退出条件后,便不再对残差进行重建,即地震频谱重建结束。重建结束后,残差rnX应被视为噪音数据对其进行处理,如直接滤掉。
总之,本发明提供的地震频谱重建方法采用最大匹配投影原理自动对频谱进行重建,子波频谱的振幅和峰值频率可自适应的调节,获得的重建频谱的和更能接近原始频谱;并且重建过程不包含方程求解,稳定性更好;重建过程不受人为因素影响,可方便的应用于较大的地震数据,且计算效率较高。
上述技术方案只是本发明的一种实施方式,对于本领域内的技术人员而言,在本发明公开了应用方法和原理的基础上,很容易做出各种类型的改进或变形,而不仅限于本发明上述具体实施方式所描述的方法,因此前面描述的方式只是优选的,而并不具有限制性的意义。
Claims (6)
1.一种地震频谱重建方法,其特征在于:所述方法包含以下步骤:
(1)获取地震信号x(t),对所述地震信号x(t)进行傅里叶变换,获得该地震信号的振幅谱X(f);
(2)选取子波频谱公式,拾取子波频谱的峰值频率m的初始值;
(3)计算子波频谱的峰值振幅a;
(4)根据最大匹配投影原理对子波频谱的峰值振幅a和子波频谱的峰值频率m进行局部优化得到局部最优的子波频谱的峰值振幅和子波频谱的峰值频率;
(5)将子波频谱的峰值振幅a和子波频谱的峰值频率m带入子波频谱公式获得子波的振幅谱;利用局部最优的子波频谱的峰值振幅a和子波频谱的峰值频率m对地震信号的振幅谱X(f)进行重建,获得重建后的地震频谱以及残差;
(6)判断所述残差是否满足给定条件,如果是,则转入步骤(7),如果否,则返回步骤(2);
(7)结束。
2.根据权利要求1所述的地震频谱重建方法,其特征在于:所述步骤(2)中子波频谱公式如下:
公式(1),
其中a为该子波频谱的峰值振幅、m为该频谱的峰值频率,f为频率。
3.根据权利要求2所述的地震频谱重建方法,其特征在于:所述步骤(3)是利用下式计算子波频谱的峰值振幅a:
公式(3)
其中为子波振幅谱集合,<riX,Ri>表示riX和Ri的内积,为第i次重建振幅谱的标准化;如果i=1,ri-1X=r0X即初始振幅谱X;如果i>1,则ri-1X为第i次重建后的残差;argmax表示取最大投影。
4.根据权利要求3所述的地震频谱重建方法,其特征在于:所述步骤(4)是这样实现的:
根据公式(3)估算参数ai:给定参数ai的范围,在该范围内取多值作为ai,将各值与峰值频率m的初始值mi一起带入公式(1)获得多个子波的振幅谱,这些子波的振幅谱集合即为在该子波振幅谱集合中确定满足公式(3)的一个子波振幅谱Ri(mi,f),该振幅谱对应的振幅即为估算的参数ai;其中参数ai的范围为n*Peak~Peak,其中n∈[0.1,0.9],Peak为riX的峰值振幅的值;
在获得初始的估算的参数ai和mi后,在局部范围内搜索参数集ξ,局部范围可表示为[ξ-Δξ,ξ+Δξ],其中Δa为振幅间隔,Δm为峰值频率间隔;在局部范围内,根据公式(3)确定最优的ai和mi,获得最优的参数集ξi,具体实现方法如下:在局部范围内取多个值作为ai和mi,形成多组重建参数ξi={ai,mi},将每一组参数ξi带入公式(1)获得子波的振幅谱,然后寻找一组满足公式(3)参数ξi,即为第i次重建的参数ξi={ai,mi},即得到局部最优的子波频谱的峰值振幅和子波频谱的峰值频率。
5.根据权利要求4所述的地震频谱重建方法,其特征在于:所述步骤(5)中的利用局部最优的子波频谱的峰值振幅和子波频谱的峰值频率对地震信号的振幅谱X(f)进行重建,获得残差是这样实现的:
重建后的振幅谱的形式为公式(2)所定义的形式:
公式(2),
其中n为重建次数、ai为与第i次重建的频谱Ri对应的局部最优的子波频谱的峰值振幅、mi为与第i次重建的频谱Ri对应的局部最优的子波频谱的峰值频率,rnX为第n次重建后的残差,r0X=X,标记i为1到n的正整数,为第i次重建的子波振幅谱,对于每一次重建,都有
6.根据权利要求5所述的地震频谱重建方法,其特征在于:所述步骤(6)中的给定条件为:
所述残差rn+1X的能量应小于等于原始振幅谱能量r0X的5%,r0X=X。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510182064.2A CN106154318B (zh) | 2015-04-16 | 一种地震频谱重建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510182064.2A CN106154318B (zh) | 2015-04-16 | 一种地震频谱重建方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106154318A true CN106154318A (zh) | 2016-11-23 |
CN106154318B CN106154318B (zh) | 2018-08-31 |
Family
ID=
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109946740A (zh) * | 2019-03-01 | 2019-06-28 | 成都理工大学 | 一种基于宽平谱地震子波整形的地震分辨率增强技术 |
CN110069836A (zh) * | 2019-04-03 | 2019-07-30 | 河海大学 | 一种高低频段交替与目标谱匹配的改进影响矩阵方法 |
CN112578438A (zh) * | 2019-09-29 | 2021-03-30 | 中国石油化工股份有限公司 | 一种地震子波提取方法及系统 |
CN113703049A (zh) * | 2021-08-27 | 2021-11-26 | 中铁二十局集团安哥拉国际有限责任公司 | 一种基于n次傅里叶谱的子波估计方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008001334A3 (en) * | 2006-06-30 | 2008-02-28 | Spectraseis Ag | Signal integration measure for seismic data |
CN102590858A (zh) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | 基于宽频子波重构的双程波成像方法 |
CN103048690A (zh) * | 2012-12-11 | 2013-04-17 | 成都理工大学 | 基于最优地震子波提取的快速匹配投影分解的地层反射拾取技术 |
CN103913771A (zh) * | 2014-04-01 | 2014-07-09 | 中国石油大学(华东) | 地震数据处理方法、装置和系统 |
CN103995289A (zh) * | 2014-05-19 | 2014-08-20 | 中国石油大学(华东) | 基于时频谱模拟的时变混合相位地震子波提取方法 |
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008001334A3 (en) * | 2006-06-30 | 2008-02-28 | Spectraseis Ag | Signal integration measure for seismic data |
CN102590858A (zh) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | 基于宽频子波重构的双程波成像方法 |
CN103048690A (zh) * | 2012-12-11 | 2013-04-17 | 成都理工大学 | 基于最优地震子波提取的快速匹配投影分解的地层反射拾取技术 |
CN103913771A (zh) * | 2014-04-01 | 2014-07-09 | 中国石油大学(华东) | 地震数据处理方法、装置和系统 |
CN103995289A (zh) * | 2014-05-19 | 2014-08-20 | 中国石油大学(华东) | 基于时频谱模拟的时变混合相位地震子波提取方法 |
Non-Patent Citations (3)
Title |
---|
YIHUA CAI ETL.: "Spectral recomposition using separable nonlinear least squares", 《SEG LAS VEGAS 2012 ANNUAL MEETING》 * |
李娜 等: "频谱重构技术及其在储层预测中的应用", 《中国地球科学联合学术年会 2014》 * |
黄跃 等: "多子波分解与重构中子波的优选", 《石油物探》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109946740A (zh) * | 2019-03-01 | 2019-06-28 | 成都理工大学 | 一种基于宽平谱地震子波整形的地震分辨率增强技术 |
CN109946740B (zh) * | 2019-03-01 | 2020-06-30 | 成都理工大学 | 一种基于宽平谱地震子波整形的地震分辨率增强方法 |
CN110069836A (zh) * | 2019-04-03 | 2019-07-30 | 河海大学 | 一种高低频段交替与目标谱匹配的改进影响矩阵方法 |
CN112578438A (zh) * | 2019-09-29 | 2021-03-30 | 中国石油化工股份有限公司 | 一种地震子波提取方法及系统 |
CN113703049A (zh) * | 2021-08-27 | 2021-11-26 | 中铁二十局集团安哥拉国际有限责任公司 | 一种基于n次傅里叶谱的子波估计方法 |
CN113703049B (zh) * | 2021-08-27 | 2024-04-23 | 中铁二十局集团安哥拉国际有限责任公司 | 一种基于n次傅里叶谱的子波估计方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103995289B (zh) | 基于时频谱模拟的时变混合相位地震子波提取方法 | |
CN107102356B (zh) | 基于ceemd的地震信号高分辨率处理方法 | |
CN101515200B (zh) | 基于瞬态视觉诱发脑电的目标选择方法 | |
CN106405645B (zh) | 一种基于资料品质分析的信噪比可控的地震拓频处理方法 | |
CN106249299A (zh) | 强反射屏蔽下薄层弱反射地震能量恢复方法及装置 | |
CN106291700A (zh) | 基于同步挤压变换的地震加权平均瞬时频率提取方法 | |
CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
CN107272062A (zh) | 一种数据驱动的地下介质q场估计方法 | |
CN106483563B (zh) | 基于互补集合经验模态分解的地震能量补偿方法 | |
CN103728663A (zh) | 一种时频分析方法 | |
CN109946739A (zh) | 一种基于压缩感知理论的地震剖面增强方法 | |
CN106873036A (zh) | 一种基于井震结合的去噪方法 | |
CN107703546A (zh) | 一种基于小波变换的新阈值函数地震资料去噪方法 | |
CN108427146B (zh) | 一种识别沉积地层中的米兰科维奇周期的方法及系统 | |
CN111060961B (zh) | 基于多信息约束反演的品质因子确定方法、装置及系统 | |
Nilsen et al. | Are there multiple scaling regimes in Holocene temperature records? | |
CN110989020B (zh) | 一种音频大地电磁数据噪声干扰的滤波方法及系统 | |
CN104635264B (zh) | 叠前地震数据的处理方法及设备 | |
CN104143031A (zh) | 一种基于小波多尺度分解的植被指数时序数据重构方法 | |
CN104422956B (zh) | 一种基于稀疏脉冲反演的高精度地震谱分解方法 | |
Du et al. | Study on optical fiber gas-holdup meter signal denoising using improved threshold wavelet transform | |
CN107703548B (zh) | 基于沉积物品质因子和回波损失级曲线峰谷的浅地层层界划分方法 | |
CN107728213A (zh) | 一种小波新阈值函数地震资料去噪方法 | |
CN106154318A (zh) | 一种地震频谱重建方法 | |
CN106154318B (zh) | 一种地震频谱重建方法 |
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 |