CN111766624B - 一种地震数据拓频处理方法、装置、储存介质及电子设备 - Google Patents

一种地震数据拓频处理方法、装置、储存介质及电子设备 Download PDF

Info

Publication number
CN111766624B
CN111766624B CN202010575806.9A CN202010575806A CN111766624B CN 111766624 B CN111766624 B CN 111766624B CN 202010575806 A CN202010575806 A CN 202010575806A CN 111766624 B CN111766624 B CN 111766624B
Authority
CN
China
Prior art keywords
wavelet
frequency
data
reflection coefficient
seismic data
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
Application number
CN202010575806.9A
Other languages
English (en)
Other versions
CN111766624A (zh
Inventor
霍守东
舒梦珵
梁瑶
舒国旭
石太昆
徐海
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Garboro Information Technology Co ltd
China University of Petroleum Beijing
Original Assignee
Nanjing Garboro Information Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Nanjing Garboro Information Technology Co ltd filed Critical Nanjing Garboro Information Technology Co ltd
Priority to CN202010575806.9A priority Critical patent/CN111766624B/zh
Publication of CN111766624A publication Critical patent/CN111766624A/zh
Application granted granted Critical
Publication of CN111766624B publication Critical patent/CN111766624B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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
    • G01V1/282—Application of seismic models, synthetic seismograms
    • 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
    • G01V1/30—Analysis
    • 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
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity

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

一种地震数据拓频处理方法、装置、储存介质及电子设备
技术领域
本发明涉及地址勘测技术领域，尤其涉及一种地震数据拓频处理方法、装置、储存介质及电子设备。
背景技术
随着勘探程度的日益加深，对于薄层、精细构造等识别能力的要求日益提高。地震数据是有效描述地下空间介质的手段。但是在地震数据采集过程中，由于受到激发条件和地层吸收的影响，得到的地震主频往往较低；而且在处理过程中，基于褶积模型的反射波处理方法受到其分辨率的影响，分辨能力仅能够达到1/4波长。
随着高分辨处理技术的发展，目前已经发展了测井约束反演、拓频处理、反褶积、小波分析等等方法，还有时频域联合反演反射系数方法、匹配追踪滤波方法、基于盲源反褶积等等技术。
但是这些方法由于未充分利用反射系数的稀疏性和全频段性，对于复杂地区的拓频能力十分有限，且在大多地区难以做到保幅保真处理。提高地震数据分辨率的同时，得到保幅保真真高分辨率的地震数据是亟待解决的问题。
发明内容
(一)发明目的
本发明的目的是提供一种地震数据拓频处理方法、装置、储存介质及电子设备以解决现有技术分辨率低且不保真的问题。
(二)技术方案
为解决上述问题，本发明的第一方面提供了一种地震数据拓频处理方法，包括：
将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
将所述地震数据进行拓展频率处理，同时得到高频子波数据；
将所述高频子波数据与所述反射系数模型褶积得到地震数据拓频结果。
进一步地，所述将地震数据经过薄层匹配追踪处理，得到初始反射系数模型具体为：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)w(ta)-b(t)w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure BDA0002550939140000021
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
进一步地，所述将所述初始反射系数模型进行基于L0范数的稀疏反演处理，得到精准反射系数模型具体为：
目标方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；
通过求解L0范数得到精确反射系数模型R(t)。
进一步地，所述将所述地震数据进行拓展频率处理，同时得到高频子波数据具体包括：
将所述地震数据进行整形正则化处理，得到初低频波数据；
将所述初始子波数据经过连续小波变换和反变换处理，得到高频子波数据。
进一步地，所述将所述地震数据进行整形正则化处理，得到初始子波数据具体为：
利用求解目标函数：
obj＝argmin||d(t)-R(t)w0(t)||2
待求参数初始子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)w0(t)||2s.t.w0(t)＝w测井(t)。
进一步地，所述将所述初始子波数据经过连续小波变换处理，得到第一子波数据具体为：
Figure BDA0002550939140000031
式中，ψ(w)为第一子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
进一步地，所述采用小波母函数在所述小波域对第一子波数据进行子波替换，将替换后的子波进行逆变换处理具体为：
在小波域，利用不同频率的子波W(w1)对第一子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(w1)代表不同尺度c的子波ω(t)在小波域的表达形式，mi为各个子波的系数；
求解得到m向量(即m0,m1,m2......)；
增加m中高频尺度成分的比重，形成新的系数向量m’(即m'0,m'1,m'2......)，最终得到小波域拟合的高频子波ψHF(w)＝m'0W(w0)+m'1W(w1)+m'2W(w2)......；通过逆变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
根据本发明的另一方面，提供一种地震数据拓频处理装置，包括：
薄层匹配追踪模块，用于将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
稀疏反演模块，用于将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
拓展频率模块，用于将所述地震数据进行拓展频率处理，同时得到高频子波数据；
褶积模块，用于将所述高频子波数据与所述反射系数模型褶积得到地震数据拓频结果。
进一步地，所述薄层匹配追踪模块具体执行下述逻辑运算：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)w(ta)-b(t)w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure BDA0002550939140000051
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
进一步地，所述稀疏反演模块具体执行下述逻辑运算：
目标方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；
通过求解L0范数得到精确反射系数模型R(t)。
进一步地，所述拓展频率模块具体包括：
整形正则化单元，用于将所述地震数据进行整形正则化处理，得到低频子波数据；
连续小波变换单元，用于将所述初始子波数据经过连续小波变换处理，得到高频子波数据。
进一步地，所述整形正则化单元具体执行下述逻辑运算：
利用求解目标函数：
obj＝argmin||d(t)-R(t)w0(t)||2
待求参数初始子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)w0(t)||2s.t.w0(t)＝w测井(t)。
进一步地，所述连续小波变换单元具体执行下述逻辑运算：
Figure BDA0002550939140000061
式中，ψ(w)为第一子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
进一步地，所述逆变换模块具体执行下述逻辑运算：
在小波域，利用不同频率的子波W(w1)对第一子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(w1)代表不同尺度c的子波ω(t)在小波域的表达形式，mi为各个子波的系数；
求解得到m向量(即m0,m1,m2......)；
增加m中高频尺度成分的比重，形成新的系数向量m’(即m'0,m'1,m'2......)，最终得到小波域拟合的高频子波ψHF(w)＝m'0W(w0)+m'1W(w1)+m'2W(w2)......；；通过逆变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
根据本发明的又一方面，提供一种计算机存储介质，所述存储介质上存储有计算机程序，所述程序被处理器执行时实现上述技术方案中任意一项所述方法的步骤。
根据本发明的又一方面，提供一种电子设备，包括存储器、显示器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述程序时实现上述技术方案中任意一项所述方法的步骤。
(三)有益效果
本发明的上述技术方案具有如下有益的技术效果：
本发明方法及装置算法稳定性高，得到的反射系数更加稀疏可靠；新剖面频谱拓展而不失自然，在处理过程中不损失数据振幅特征。
附图说明
图1是根据本发明第一实施方式的地震数据拓频处理方法流程图；
图2是根据本发明一具体实施方式的初始的反射系数模型对地震剖面的处理前后对比图；
图3是根据本发明一具体实施方式的基于压缩感知的叠后稀疏反演对地震数据的处理图；
图4是根据本发明一具体实施方式的初始的反射系数模型和精准反射系数模型对地震数据处理前后对比图；
图5是根据本发明一具体实施方式的Morlet小波时间域(左)和频率域(右)示意图；
图6是根据本发明一具体实施方式的Morlet小波替换示意图(上：低频小波；下：高频小波)；
图7是根据本发明一具体实施方式的提频效果及反射系数对比图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
此外，下面所描述的本发明不同实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互结合。
在本发明实施例的第一方面，提供了一种地震数据拓频处理方法，包括：
S1：将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
S2：将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
S3：将所述地震数据进行拓展频率处理，同时得到高频子波数据；
S4：将所述高频子波数据与所述反射系数模型褶积得到地震数据拓频结果。
上述方法稳定性高，得到的反射系数更加稀疏可靠；新剖面频谱拓展不失自然，在处理过程中不损失数据振幅特征。
可选的，所述将地震数据经过薄层匹配追踪处理，得到初始反射系数模型具体为：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)w(ta)-b(t)w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure BDA0002550939140000091
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
可选的，所述将所述初始反射系数模型进行基于L0范数的稀疏反演处理，得到精准反射系数模型具体为：
目标方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；
通过求解L0范数得到精确反射系数模型R(t)。
可选的，所述将所述地震数据进行拓展频率处理，同时得到高频子波数据具体包括：
将所述地震数据进行整形正则化处理，得到低频子波数据；
将所述初始子波数据经过连续小波变换和反变换处理，得到高频子波数据。
可选的，所述将所述地震数据进行整形正则化处理，得到初始子波数据具体为：
利用求解目标函数：
obj＝argmin||d(t)-R(t)w0(t)||2
待求参数初始子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)w0(t)||2s.t.w0(t)＝w测井(t)。
可选的，所述将所述初始子波数据经过连续小波变换处理，得到第一子波数据具体为：
Figure BDA0002550939140000101
式中，ψ(w)为第一子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
可选的，所述采用小波母函数在所述小波域对第一子波数据进行子波替换，将替换后的子波进行逆变换处理具体为：
在小波域，利用不同频率的子波W(w1)对第一子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(w1)代表不同尺度c的子波ω(t)在小波域的表达形式，mi为各个子波的系数；
求解得到m向量(即m0,m1,m2......)；
增加m中高频尺度成分的比重，形成新的系数向量m’(即m'0,m'1,m'2......)，最终得到小波域拟合的高频子波ψHF(w)＝m'0W(w0)+m'1W(w1)+m'2W(w2)......；通过逆变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
根据本发明的另一方面，提供一种地震数据拓频处理装置，包括：
薄层匹配追踪模块，用于将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
稀疏反演模块，用于将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
拓展频率模块，用于将所述地震数据进行拓展频率处理，同时得到高频子波数据；
褶积模块，用于将所述高频子波数据与所述反射系数模型褶积得到地震数据拓频结果。
可选的，所述薄层匹配追踪模块具体执行下述逻辑运算：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)w(ta)-b(t)w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure BDA0002550939140000111
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
可选的，所述稀疏反演模块具体执行下述逻辑运算：
目标方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；
通过求解L0范数得到精确反射系数模型R(t)。
可选的，所述拓展频率模块具体包括：
整形正则化单元，用于将所述地震数据进行整形正则化处理，得到低频子波数据；
连续小波变换单元，用于将所述初始子波数据经过连续小波变换处理，得到高频子波数据。
可选的，所述整形正则化单元具体执行下述逻辑运算：
利用求解目标函数：
obj＝argmin||d(t)-R(t)w0(t)||2
待求参数初始子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)w0(t)||2s.t.w0(t)＝w测井(t)。
可选的，所述连续小波变换单元具体执行下述逻辑运算：
Figure BDA0002550939140000121
式中，ψ(w)为第一子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
可选的，所述逆变换模块具体执行下述逻辑运算：
在小波域，利用不同频率的子波W(w1)对第一子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(w1)代表不同尺度c的子波ω(t)在小波域的表达形式，mi为各个子波的系数；
求解得到m向量(即m0,m1,m2......)；
增加m中高频尺度成分的比重，形成新的系数向量m’(即m'0,m'1,m'2......)，最终得到小波域拟合的高频子波ψHF(w)＝m'0W(w0)+m'1W(w1)+m'2W(w2)......；通过逆变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
在本发明一具体实施例中，提供基于压缩感知的提高分辨率方法，包括的核心技术流程包括：
1.薄层匹配追踪
(1)薄层匹配追踪–双层假设
假设两个单反射层位于时间ta/tb，反射系数为a/b
那么，目标方程为
obj＝∑(d(t)-aw(ta)-bw(tb))2
求极值
Figure BDA0002550939140000131
线性方程的解
Figure BDA0002550939140000141
(2)薄层匹配追踪实例
通过薄层匹配追踪，可以得到较为初始的反射系数模型，如图2所示，左侧为输入地震剖面，右侧为输出反射系数。
2.基于压缩感知的叠后稀疏反演–L0Norm
(1)原理
L0Norm就是计算非零值的个数，目标方程为
obj＝∑(d-w*r)2+λC(r≠0)
L0范数的求解是个非凸函数，只有局部最优解(NPhard)，一般有两种办法，1是L1近似，2是OMP匹配追踪；两者各有优缺点，不同的假设。我们以薄层匹配追踪的结果作为输入，直接求解L0范数问题，得到了很好的结果。
(2)L0Norm示例如图3所示，从左向右依次为输入地震道、L0Norm、L1Norm、L2Norm。通过L0范数反演示例可以看出，L0范数结果更加稀疏，信噪比更高，能够得到频带更宽的地震反射系数模型，为进一步提频处理提供较好的结果。
(3)叠后稀疏反演
将匹配追踪结果作为初始输入，通过L0范数叠后稀疏反演得到更加准确的反射系数模型，如图4所示，从左向右依次为输入地震剖面、匹配追踪结果、稀疏反演结果。
3子波提取-ShapingRegularization
对于线性方程y＝Ax，假设x光滑，我们可以得到
x＝[S(ATA-λ2I)+λ2I]-1SATy
如果遇到了除零问题，我们可以这样解决，由
Figure BDA0002550939140000151
可知
Figure BDA0002550939140000152
而对于简单的地震褶积模型
d(t)＝r(t)*w(t)
在频率域是
D(f)＝R(f)W(f)
子波频谱
Figure BDA0002550939140000153
在子波提取过程中，如果我们假设子波的频谱是光滑的，我们可以采用上述的shapingregularization方法，由
Figure BDA0002550939140000154
可知
Figure BDA0002550939140000155
与此同时，如果有vsp(垂直地震剖面(verticalseismicprofile))地震数据和测井数据，可以将它们得到的地震子波信息，用于约束子波提取。
4.连续小波变换–Morlet小波
(1)Morlet小波
高斯包络下的单频率复正弦函数
Figure BDA0002550939140000161
如图5所示，左侧为Morlet小波时间域右侧为频率域。
(2)Morlet子波替换
我们利用Morlet小波特性，使得宽频子波与原子波的峰值，相位完全吻合，自然的拓展低、高频能量，使得新剖面频谱拓展而不失自然。通过给定的参数(主频，频宽)来替代小波母函数，如图6所示，上方为低频小波，下方为高频小波。
(3)连续小波变换–逆变换
连续小波变换逆变换容许条件(Smith，2008)
Figure BDA0002550939140000162
逆变换公式
Figure BDA0002550939140000163
通过连续小波逆变换，得到最终宽频处理结果。如图7所示，从左向右依次为原始数据；初始反射系数；提频结果和最终反射系数。
上述实施例将薄层匹配追踪算法得到的反射系数作为初始反射系数，用于基于压缩感知的叠后稀疏反演，保证了算法稳定性，得到的反射系数更加稀疏可靠。子波提取和替换过程：首先采用了ShapingRegularization技术进行子波提取，然后利用利用Morlet小波特性，使得宽频子波与原子波的峰值，相位完全吻合，自然的拓展低、高频能量，使得新剖面频谱拓展而不失自然。通过给定的参数(主频，频宽)来替代小波母函数。保证了在处理过程中不损失数据振幅特征。
在本发明实施例的又一方面，提供一种计算机存储介质，所述存储介质上存储有计算机程序，所述程序被处理器执行时实现上述实施例中任意一项所述方法的步骤。
在本发明实施例的又一方面，提供一种电子设备，包括存储器、显示器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述程序时实现上述实施例中任意一项所述方法的步骤。
本发明旨在保护一种地震数据拓频处理方法，包括：将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；将所述地震数据进行拓展频率处理，同时得到高频子波数据；将所述高频子波数据与所述反射系数模型褶积得到地震数据拓频结果。该方法稳定性高，得到的反射系数更加稀疏可靠；新剖面频谱拓展不失自然，在处理过程中不损失数据振幅特征。
应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (14)

1.一种地震数据拓频处理方法，其特征在于，包括：
将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
将所述地震数据进行拓展频率处理，同时得到高频子波数据；
将所述高频子波数据与所述精准反射系数模型褶积得到地震数据拓频结果；
所述将所述地震数据进行拓展频率处理，同时得到高频子波数据具体包括：
将所述地震数据进行整形正则化处理，得到低频子波数据；
将所述低频子波数据经过连续小波变换和反变换处理，得到高频子波数据。
2.根据权利要求1所述的方法，其特征在于，所述将地震数据经过薄层匹配追踪处理，得到初始反射系数模型具体为：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)*w(ta)-b(t)*w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure FDA0003128718110000021
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
3.根据权利要求2所述的方法，其特征在于，所述将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；具体为：将所述初始反射系数模型进行基于L0范数的稀疏反演处理，得到精准反射系数模型；其稀疏反演方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)*w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；w(t)为待求参数低频子波数据；
通过求解L0范数得到精确反射系数模型R(t)。
4.根据权利要求1所述的方法，其特征在于，所述将所述地震数据进行整形正则化处理，得到低频子波数据具体为：
利用求解目标函数：
obj＝argmin||d(t)-R(t)*w0(t)||2
待求参数低频子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)*w0(t)||2 s.t.w0(t)＝w测井(t)。
5.根据权利要求1所述的方法，其特征在于，所述将所述低频子波数据经过连续小波变换和反变换处理，得到高频子波数据具体为：
Figure FDA0003128718110000031
式中，ψ0(w)为高频子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
6.根据权利要求5所述的方法，其特征在于，所述低频子波数据经过连续小波变换和反变换处理，得到高频子波数据，所述反变换处理具体为：
在小波域，利用不同频率的子波W(wi)对高频子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(wi)代表不同尺度c的子波ω(t)在小波域的表达形式，即W(w0)，W(w1)，W(w2)......，mi为各个子波的系数；求解得到m向量，即m0,m1,m2.....；
增加m中高频尺度成分的比重，形成新的系数向量m’，即m′0,m′1,m′2....，最终得到小波域拟合的高频子波：
ψHF(w)＝m'0W(w0)+m′1W(w1)+m'2W(w2)......；
通过反变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
7.一种地震数据拓频处理装置，其特征在于，包括：
薄层匹配追踪模块，用于将地震数据经过薄层匹配追踪处理，得到初始反射系数模型；
稀疏反演模块，用于将所述初始反射系数模型进行稀疏反演处理，得到精准反射系数模型；
拓展频率模块，用于将所述地震数据进行拓展频率处理，同时得到高频子波数据；
褶积模块，用于将所述高频子波数据与所述精准反射系数模型褶积得到地震数据拓频结果；
所述拓展频率模块具体包括：
整形正则化单元，用于将所述地震数据进行整形正则化处理，得到低频子波数据；
连续小波变换单元和反变换模块，用于将所述低频子波数据经过连续小波变换和反变换处理，得到高频子波数据。
8.根据权利要求7所述的装置，其特征在于，所述薄层匹配追踪模块具体执行下述逻辑运算：
假设两个单反射层位于时间ta和tb，反射系数为a和b，其中a,b均为时间t的时间序列，w(ta)和w(tb)分别代表两个基函数，该基函数为小波基或科沃莱特基函数；
从目标方程求解a和b，式中d(t)为地震数据：
obj＝argmin||d(t)-a(t)w(ta)-b(t)w(tb)||2；
argmin||·||2代表求解L2范数；
能够得到解
Figure FDA0003128718110000051
得到初始反射系数模型r(t)＝a(t)+b(t)，其中，t为时间。
9.根据权利要求7所述的装置，其特征在于，所述稀疏反演模块具体执行下述逻辑运算：
目标方程为：
obj＝argmin||R(t)||0 s.t.d(t)＝R(t)*w(t)；
d(t)为输入，求解L0范数，得到精准反射系数模型R(t)；
argmin||·||0代表求解L0范数，s.t.代表约束条件，求解该方程需要迭代算法，第一次迭代时将所述初始反射系数模型r(t)＝a(t)+b(t)作为初始反射系数：R(t)＝r(t)；w(t)为待求参数低频子波数据；
通过求解L0范数得到精确反射系数模型R(t)。
10.根据权利要求9所述的装置，其特征在于，所述整形正则化单元具体执行下述逻辑运算：
利用求解目标函数：
obj＝argmin||d(t)-R(t)*w0(t)||2
待求参数低频子波数据为w0(t)，输入为d(t)和R(t)，*代表褶积运算；
采用整形正则化方法求解，
与此同时，将垂直地震剖面的地震数据或测井数据的地震子波信息w测井(t)约束子波提取
obj＝argmin||d(t)-R(t)*w0(t)||2 s.t.w0(t)＝w测井(t)。
11.根据权利要求9所述的装置，其特征在于，所述连续小波变换单元具体执行下述逻辑运算：
Figure FDA0003128718110000061
式中，ψ0(w)为高频子波数据，c为小波变换的尺度，τ为平移量，t为时间；ω为小波基函数。
12.根据权利要求7所述的装置，其特征在于，所述反变换模块具体执行下述逻辑运算：
在小波域，利用不同频率的子波W(w1)对高频子波数据ψ0(w)进行拟合
argmin||ψ0(w)-m0W(w0)-m1W(w1)-m2W(w2)......||2；
式中W(wi)代表不同尺度c的子波ω(t)在小波域的表达形式，即W(w0)，W(w1)，W(w2)......，mi为各个子波的系数；
求解得到m向量，即m0,m1,m2......；
增加m中高频尺度成分的比重，形成新的系数向量m’，即m'0,m′1,m'2....，最终得到小波域拟合的高频子波：
ψHF(w)＝m'0W(w0)+m′1W(w1)+m'2W(w2)......；
通过反变换，得到拓频后的高频第二子波wHF(t)；最终得到拓频后地震数据：
D(t)＝wHF(t)*R(t)。
13.一种计算机存储介质，其特征在于，所述存储介质上存储有计算机程序，所述程序被处理器执行时实现权利要求1-6中任意一项所述方法的步骤。
14.一种电子设备，其特征在于，包括存储器、显示器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述程序时实现权利要求1-6中任意一项所述方法的步骤。
CN202010575806.9A 2020-06-22 2020-06-22 一种地震数据拓频处理方法、装置、储存介质及电子设备 Active CN111766624B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010575806.9A CN111766624B (zh) 2020-06-22 2020-06-22 一种地震数据拓频处理方法、装置、储存介质及电子设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010575806.9A CN111766624B (zh) 2020-06-22 2020-06-22 一种地震数据拓频处理方法、装置、储存介质及电子设备

Publications (2)

Publication Number Publication Date
CN111766624A CN111766624A (zh) 2020-10-13
CN111766624B true CN111766624B (zh) 2021-10-08

Family

ID=72721596

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010575806.9A Active CN111766624B (zh) 2020-06-22 2020-06-22 一种地震数据拓频处理方法、装置、储存介质及电子设备

Country Status (1)

Country Link
CN (1) CN111766624B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117434604A (zh) * 2023-02-20 2024-01-23 中国石油化工股份有限公司 地震数据处理方法、装置、存储介质及电子设备

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6789018B1 (en) * 2003-08-29 2004-09-07 Nonlinear Seismic Imaging, Inc. Mapping reservoir rocks using frequency spectral broadening and the presence of the slow-wave
CN103336303A (zh) * 2013-06-06 2013-10-02 浙江大学 一种利用声波测井进行地震拓频方法
CN104280765A (zh) * 2013-07-11 2015-01-14 中国石油化工股份有限公司 基于变子波反射系数反演的地震高分辨处理方法
CN110174702A (zh) * 2018-09-30 2019-08-27 中海油海南能源有限公司 一种海上地震数据低频弱信号恢复的方法和系统
CN110568489A (zh) * 2019-07-23 2019-12-13 中国石油化工股份有限公司 一种块状介质的宽频反演方法
CN111025395A (zh) * 2020-01-06 2020-04-17 中国石油化工股份有限公司 基于压缩感知的高斯频率域提高薄互层分辨率的方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6789018B1 (en) * 2003-08-29 2004-09-07 Nonlinear Seismic Imaging, Inc. Mapping reservoir rocks using frequency spectral broadening and the presence of the slow-wave
CN103336303A (zh) * 2013-06-06 2013-10-02 浙江大学 一种利用声波测井进行地震拓频方法
CN104280765A (zh) * 2013-07-11 2015-01-14 中国石油化工股份有限公司 基于变子波反射系数反演的地震高分辨处理方法
CN110174702A (zh) * 2018-09-30 2019-08-27 中海油海南能源有限公司 一种海上地震数据低频弱信号恢复的方法和系统
CN110568489A (zh) * 2019-07-23 2019-12-13 中国石油化工股份有限公司 一种块状介质的宽频反演方法
CN111025395A (zh) * 2020-01-06 2020-04-17 中国石油化工股份有限公司 基于压缩感知的高斯频率域提高薄互层分辨率的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于压缩感知和稀疏反演的地震数据低频补偿;韩立国等;《吉林大学学报(地球科学版)》;20121226;第259-264页 *

Also Published As

Publication number Publication date
CN111766624A (zh) 2020-10-13

Similar Documents

Publication Publication Date Title
Gholami Sparse time–frequency decomposition and some applications
Nazari Siahsar et al. Data-driven multitask sparse dictionary learning for noise attenuation of 3D seismic data
Chen et al. Double-sparsity dictionary for seismic noise attenuation
Aghamiry et al. Compound regularization of full-waveform inversion for imaging piecewise media
Beckouche et al. Simultaneous dictionary learning and denoising for seismic data
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
Lu Adaptive multiple subtraction using independent component analysis
US10393899B2 (en) Automatic tracking of faults by slope decomposition
Dong et al. Signal-to-noise ratio enhancement for 3C downhole microseismic data based on the 3D shearlet transform and improved back-propagation neural networks
Wang et al. High-resolution wave-equation amplitude-variation-with-ray-parameter (AVP) imaging with sparseness constraints
Zhou et al. Robust noise attenuation based on nuclear norm minimization and a trace prediction strategy
CN112882099B (zh) 一种地震频带拓宽方法、装置、介质及电子设备
CN106680874A (zh) 基于波形形态特征稀疏化建模的谐波噪声压制方法
CN105093315B (zh) 一种去除煤层强反射信号的方法
CN110646841B (zh) 时变稀疏反褶积方法及系统
CN111766624B (zh) 一种地震数据拓频处理方法、装置、储存介质及电子设备
Liu et al. Eliminating harmonic noise in vibroseis data through sparsity-promoted waveform modeling
Chen et al. A misfit function based on entropy regularized optimal transport for full-waveform inversion
CN113341463B (zh) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN113109873B (zh) 一种基于秩残差约束的沙漠地区地震信号噪声抑制方法
Wang et al. A Joint Framework for Seismic Signal Denoising Using Total Generalized Variation and Shearlet Transform
Sun et al. Application of adaptive iterative low-rank algorithm based on transform domain in desert seismic signal analysis
CN112363217A (zh) 一种地震数据随机噪声压制方法及系统
Nose-Filho et al. Algorithms for sparse multichannel blind deconvolution

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20231017

Address after: 102249 Beijing city Changping District Road No. 18

Patentee after: China University of Petroleum (Beijing)

Patentee after: Nanjing Garboro Information Technology Co.,Ltd.

Address before: No.1 Dongji Avenue, Jiangning District, Nanjing City, Jiangsu Province (Jiangning Development Zone)

Patentee before: Nanjing Garboro Information Technology Co.,Ltd.