CN112213773B - 一种地震分辨率提高方法及电子设备 - Google Patents

一种地震分辨率提高方法及电子设备 Download PDF

Info

Publication number
CN112213773B
CN112213773B CN201910631191.4A CN201910631191A CN112213773B CN 112213773 B CN112213773 B CN 112213773B CN 201910631191 A CN201910631191 A CN 201910631191A CN 112213773 B CN112213773 B CN 112213773B
Authority
CN
China
Prior art keywords
seismic
spectrum
frequency
seismic data
reflection coefficient
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
CN201910631191.4A
Other languages
English (en)
Other versions
CN112213773A (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.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
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 China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201910631191.4A priority Critical patent/CN112213773B/zh
Publication of CN112213773A publication Critical patent/CN112213773A/zh
Application granted granted Critical
Publication of CN112213773B publication Critical patent/CN112213773B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364Seismic filtering
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms

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,根据获取到的目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率表达式,对目标地震褶积模型的边缘似然函数进行最大化处理,得到第一待反演伪反射系数估计值;
步骤2,将所述第一待反演伪反射系数估计值代入所述目标后验概率表达式进行迭代计算,得到第二待反演伪反射系数估计值和当前迭代变量,并获得当前迭代次数,判断该当前迭代次数是否大于预设迭代阈值,以及判断当前迭代变量的值是否小于预设值;
步骤3,当所述当前迭代次数大于预设迭代阈值且所述当前迭代变量的值小于预设值时,将所述第二待反演伪反射系数估计值作为新的第一待反演伪反射系数估计值,并返回执行步骤2,直至所述当前迭代次数等于预设迭代阈值或所述当前迭代变量的值小于预设值时,将获得的第二待反演伪反射系数估计值作为谱约束伪反射系数估计值;
步骤4,根据所述谱约束伪反射系数估计值对获取到的地震数据进行处理,得到地震数据谱约束伪反射系数频谱,利用测井数据对所述地震数据谱约束伪反射系数频谱进行差异性滤波,以确定地震数据谱约束伪反射系数可信频谱范围;
步骤5,对所述谱约束伪反射系数可信频谱进行拓频处理,以获得经过拓频处理的拓频地震数据,其中,所述拓频地震数据的频带宽度大于拓频处理前的地震数据的频带宽度,从而根据拓频地震数据提高了地震分辨率。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,所述地震数据为叠后地震数据。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,在执行所述步骤1之前,所述方法还包括:
接收用户输入的目标时间段,并根据该目标时间段和地震数据中的子波构建时变地震物理子波库;
根据所述时变地震物理子波库和所述目标时间段,得到地震数据时间域表达式;
对所述地震数据时间域表达式进行傅里叶变换,以得到地震数据频率域表达式;
对所述地震数据频率域表达式在地震数据的指定频段内进行离散化处理,以得到频率域地震褶积模型;
根据所述频率域地震褶积模型和预设的谱约束伪反射系数,得到所述目标地震褶积模型;
当噪音服从零均值高斯分布时,根据贝叶斯原理和所述目标地震褶积模型得到该目标地震褶积模型的似然函数和目标地震褶积模型的边缘似然函数。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,根据该目标时间段和地震数据中的子波构建时变地震物理子波库的步骤包括:
在所述目标时间段内将随时间变化的非稳态地震数据作为稳态地震数据;
对所述稳态地震数据进行开窗处理,并在选定的时窗内采用井震标定的方法对子波进行提取;
对提取后的子波进行汇总,以得到所述时变地震物理子波库。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,所述步骤1包括:
在贝叶斯框架下,根据所述目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率,其中,所述目标后验概率表达式为:
p(m1|d,h,σ2)=p(d|m12)p(m1|h)/p(d|h,σ2)
其中,m1为谱约束伪反射系数序列,d为所述目标地震褶积模型,h为服从零均值高斯分布的谱约束伪反射系数的方差的转置矩阵,h=[h1...hK]T包含了K个独立的参数,hk为第k个谱约束伪反射系数的方差,K为选定时间段的总采样点数,且k为大于0且小于K的自然数,σ2为噪音方差,p(d|m12)为所述地震数据振幅谱的似然函数,p(m1|h)为所述自动相关判别先验信息,p(d|h,σ2)为所述地震数据振幅谱的边缘似然函数;
对所述目标地震褶积模型的边缘似然函数使用最大期望算法,以获得所述转置矩阵的估计值和所述噪音方差的估计值;
根据所述转置矩阵的估计值、噪音方差的估计值和所述目标后验概率表达式,得到所述第一待反演伪反射系数估计值。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,所述地震数据的指定频段为地震数据的主频段,所述地震数据的主频段为所述地震数据振幅谱峰值的0.707倍处的两个频率值之间的频率范围。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,所述当前迭代变量的表达式为:
Figure GDA0003468676060000041
其中,d为所述目标地震褶积模型,G=F,
Figure GDA0003468676060000042
为从时间域的δ函数δ(t-tk)根据傅里叶变换得到的频率域结果,e为自然常数2.71828,i为虚数符号,π为圆周率,t为时间,tk为第k个谱约束伪反射系数的时间位置,k为大于0且小于为选定时窗内的总采样点数的自然数,fm为待处理的第m个频段,其中,f为频率,m为大于0且小于M的自然数,M为当前非零的谱约束伪反射系数的个数,m1为谱约束伪反射系数序列。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,利用测井数据对所述地震数据谱约束伪反射系数频谱进行差异性滤波,以确定地震数据谱约束伪反射系数可信频谱范围的步骤包括:
对测井数据中的纵波速度和密度数据进行方波化处理,以得到测井反射系数,并根据测井数据中的子波和所述测井反射系数,得到测井反射系数频谱;
采用滑动频率域窗口对所述测井反射系数频谱和所述谱约束伪反射系数频谱进行计算,以得到所述测井反射系数频谱与所述谱约束伪反射系数频谱的相关系数;
根据所述相关系数在地震数据中低频和高频的突变点,采用四点梯形滤波算子选取待滤波的频谱范围,对所述待滤波的频谱范围作归一化处理,以得到归一化处理后的待滤波频谱范围,并根据所述相关系数与所述归一化处理后的待滤波频谱范围,得到最终频率域滤波器参数;
根据最终频率域滤波器参数对所述地震数据进行滤波,以得到所述地震数据谱约束伪反射系数可信频谱范围。
在本发明实施例较佳的选择中,在上述地震分辨率提高方法中,对谱约束伪反射系数的可信频谱进行拓频处理的步骤包括:对所述地震数据谱约束伪反射系数可信频谱范围进行反傅里叶变换。
同时,本发明给出了一种电子设备,包括存储器、处理器以及存储在该存储器上并可在其处理器上运行的计算机程序,该计算机程序被处理器执行时,可实现上述方法和实施例中任一项的地震分辨率提高方法。
本发明提供的地震分辨率提高方法及电子设备,通过在频率域中预设谱约束伪反射系数,结合自动相关判别先验信息,获得准确的拓展结果和有效的低频信息,解决了现有技术中对地震数据低频拓展能力低,拓频中无法兼顾信噪比和分辨率的同时提升,拓展地震数据不准确且噪音高的问题。
为使本发明的上述目的、特征和优点能更明显易懂,下文特举较佳实施例,并配合所附附图,作详细说明如下。
附图说明
通过结合附图阅读下文示例性实施例的详细描述可更好地理解本公开的范围。其中所包括的附图是:
图1为本发明实施例提供的地震分辨率提高方法流程图。
图2为本发明实施例提供的不含噪音地震数据拓频实例。
图3为本发明实施例提供的含噪音地震数据拓频实例。
图4为实际资料拓频实例。
图5为实际资料拓频前后频谱对比。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例只是本发明的一部分实施例,而不是全部的实施例。
因此,以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围,而是仅仅表示本发明的选定实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
应注意到:相似的标号和字母在下面的附图中表示类似项,因此,一旦某一项在一个附图中被定义,则在随后的附图中不需要对其进行进一步定义和解释。在本发明的描述中,术语“第一”、“第二”等仅用于区分描述,而不能理解为只是或暗示相对重要性。
请结合参阅图1和图2,本发明实施例提供了一种地震分辨率提高方法,该处理方法包括步骤S110至步骤S170:
步骤S110,接收用户输入的目标时间段,并根据该目标时间段和地震数据中的子波构建时变地震物理子波库。
其中,上述步骤S110可以是,基于高频衰减补偿法构建时变地震物理子波库,基于反Q滤波法构建时变地震物理子波库,基于谱模拟法构建时变地震物理子波库,基于最小相位假设构建时变地震物理子波库,基于分段井震标定构建时变地震物理子波库。
在本实施例中,采用基于分段井震标定构建时变地震物理子波库的方式。
具体的,在本实施例中,上述步骤S110包括步骤S112至S116。
步骤S112:在所述目标时间段内将随时间变化的非稳态地震数据作为稳态地震数据。
步骤S114:对所述稳态地震数据进行开窗处理,并在选定的时窗内采用井震标定的方法对子波进行提取。
步骤S116:对提取后的子波进行汇总,以得到所述时变地震物理子波库。
可以理解,这种分段井震标定估计子波的方法与最小相位假设获得子波的估计相比更加同时符合振幅和相位的特征,且实用性更强,更稳定,获得的子波可以更好地满足该层段的实际地质情况和需要。
同时,步骤S114还包括记录选定时窗内的总采样点数。
在本实施例中,所述地震数据为叠后地震数据。叠前地震数据与叠后地震数据是就地震数据处理中的偏移来说的,在地下介质的产状不是水平时,反射地震的同相轴会发生偏移,不能反映该处地下介质的真实产状,这时要求对地震数据进行偏移归位。偏移有两种方式:叠加之前偏移和叠加之后偏移,前者就叫叠前,后者叫叠后,一般来说,在构造比较简单的地区使用叠后地震数据可以缩短生产运行周期,进而节约开支。
步骤S120,根据所述时变地震物理子波库和预设的谱约束伪反射系数,获得所述目标地震褶积模型。
具体的,在本实施例中,上述步骤S120包括步骤S122至S128。
步骤S122:根据所述时变地震物理子波库和所述目标时间段,得到地震数据时间域表达式。
其中,所述的地震数据时间域表达式为:
Figure GDA0003468676060000071
其中s(t)为地震记录,t为时间,tk为第k个谱约束反射系数的时间位置,k为取值范围从1到K的自然数,K为选定时窗内的总采样点数,w(t)为地震物理子波,rk为第k个反射系数的幅值,δ(t-tk)为δ函数在时间t与tk的差值时间段的表达,tk为第k个反射系数的时间位置,n(t)为存在于地震资料中的噪音信号。
步骤S124:对所述地震数据时间域表达式进行傅里叶变换,以得到地震数据频率域表达式。
其中,所述地震数据频率域表达式为:
Figure GDA0003468676060000072
其中,f是频率,k为取值范围从1到K的自然数,K为选定时间段的总采样点数,rk为第k个反射系数的幅值,i为复数单位,S(f)为地震数据的频谱,W(f)为地震物理子波的频谱,e为自然常数2.71828,tk为第k个谱约束反射系数的时间位置,π为圆周率,N(f)为地震数据中噪音的频谱。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S126:对所述地震数据频率域表达式在地震数据的指定频段内进行离散化处理,以得到频率域地震褶积模型。
其中,选定频率成分fm作为后续处理的频段,对所述地震数据频率域表达式可以被离散化并写成:
Figure GDA0003468676060000073
将离散化的地震数据频率域表达式写成矩阵形式,获得频率域地震褶积模型表达式为:
d=GR+n,
其中,d为地震数据振幅谱,d=S/W,G=F,R为反射系数矩阵,R=[r1...rk]T,rk为第k个反射系数的幅值,n=N/W,S为地震数据频谱矩阵,S(fm)为地震数据在第m个频段的频谱,S=[S(f1)...S(fm)]T,W为地震数据物理子波频谱矩阵,W(fm)为地震物理子波在第m个频段的频谱,W=[W(f1)...W(fm)]T,N为地震数据中噪音频谱矩阵,N(fm)为地震数据中噪音在第m个频段的频谱,N=[N(f1)...N(fm)]T,[·]T代表矩阵的转置,F为δ函数在频率域的矩阵表达式,
Figure GDA0003468676060000081
为从时间域的δ函数δ(t-tk)根据傅里叶变换得到的频率域结果,e为自然常数2.71828,i为虚数符号,π为圆周率,t为时间,tk为第k个谱约束反射系数的时间位置,k为取值范围从1到K的自然数,K为选定时间段的总采样点数,f为频率,m为取值范围从1到选定时间段的总采样点数的自然数,n为噪音。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S128:根据所述频率域地震褶积模型和预设的谱约束伪反射系数,得到所述目标地震褶积模型。
在本实施例中,基于所述目标地震褶积模型的表达式为:
d=Gm1+n,
其中,d为地震数据振幅谱,d=S/W,G=F,n=N/W,S(fm)为地震数据在第m个频段的频谱,S=[S(f1)...S(fm)]T,W(fm)为地震物理子波在第m个频段的频谱,W=[W(f1)...W(fm)]T,N为地震数据中噪音频谱矩阵,N(fm)为地震数据中噪音在第m个频段的频谱,N=[N(f1)...N(fm)]T,[·]T代表矩阵的转置,F为δ函数在频率域的矩阵表达式,
Figure GDA0003468676060000091
为从时间域的δ函数δ(t-tk)通过傅里叶变换到频率域的结果,e为自然常数2.71828,i为虚数符号,π为圆周率,t为时间,tk为第k个谱约束反射系数的时间位置,k为取值范围从1到选定时间段的总采样点数的自然数,f为频率,m1为谱约束伪反射系数序列,n为噪音。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
在本实施例中,所述地震数据的指定频段为地震数据的主频段,所述地震数据的主频段为所述地震数据振幅谱峰值的0.707倍处的两个频率值之间的频率范围。
可以理解,地震数据主频是指地震数据振幅谱的峰值频率,即频谱曲线极大值所对应的频率,而地震数据的主频段为所述地震数据振幅谱峰值的0.707倍处的两个频率值之间的频率范围。在一般的地震信号中,主频段的信噪比较高,低频段和高频段的地震信号的信噪比较低,因此为了避免非主频段信号的低信噪比对后续处理的影响,本实施例选取地震数据的主频段进行利用。
在本实施例中,预设的谱约束伪反射系数近似于真实反射系数,由于在求解过程中会采用稀疏约束和层状介质假设,且真实反射系数中的小幅值或超薄层的反射系数是无法被恢复的,因此可以采用谱约束伪反射系数来代替真实反射系数进行后续地震数据处理。其中,谱约束伪反射系数与真实反射系数具有一样的表达形式。
步骤S130,根据获取到的目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率表达式,对目标地震褶积模型的边缘似然函数进行最大化处理,得到第一待反演伪反射系数估计值。
具体的,在本实施例中,上述步骤S130包括步骤S131至S135。
步骤S131,将自动相关判别先验信息作为谱约束待反演伪反射系数的先验分布。
在本实施例中,所述自动相关判别先验信息表达式,即谱约束待反演伪反射系数的先验分布表达式为:
Figure GDA0003468676060000101
其中,m1为谱约束伪反射系数序列,h为服从零均值高斯分布的谱约束伪反射系数的方差的转置矩阵,h=[h1...hK]T包含了K个独立的参数,这代表每一个谱约束伪反射系数都服从零均值高斯分布,hk为第k个谱约束伪反射系数的方差,其中,K为选定时间段的总采样点数,k为大于0且小于K的自然数,e为自然常数2.71828,i为虚数符号,π为圆周率,rk为第k个反射系数的幅值。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S132,当噪音服从零均值高斯分布时,根据贝叶斯原理和所述目标地震褶积模型得到该目标地震褶积模型的似然函数和目标地震褶积模型的边缘似然函数。
在本实施例中,所述地震数据振幅谱的似然函数表达式为:
Figure GDA0003468676060000102
其中,d为地震数据振幅谱,m1为谱约束伪反射系数序列,σ2为噪音方差,π为圆周率,M为当前非零谱约束伪反射系数的个数,E为噪音方差矩阵,E=σ2I,I为单位矩阵,e为自然常数2.71828,G=F,F为δ函数在频率域的矩阵表达式,[·]T为矩阵的转置。
在本实施例中,所述地震数据振幅谱的边缘似然函数表达式为:
Figure GDA0003468676060000103
其中,d为地震数据振幅谱,h为服从零均值高斯分布的谱约束伪反射系数的方差的转置矩阵,σ2为噪音方差,π为圆周率,M为当前非零谱约束伪反射系数的个数,Q=σ2I+GH-1GH,e为自然常数2.71828,G=F,F为δ函数在频率域的矩阵表达式,H为服从零均值高斯分布的谱约束伪反射系数的方差的对角矩阵,H=diag(h1...hk)包含了k个独立的参数,这代表每一个谱约束伪反射系数都服从零均值高斯分布,hk为第k个谱约束伪反射系数的方差,k为选定时窗内的总采样点数,K为选定时间段的总采样点数,[·]T为矩阵的转置。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S133,在贝叶斯框架下,根据所述目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率,其中,所述目标后验概率表达式为:
Figure GDA0003468676060000111
其中,m1为谱约束伪反射系数序列,d为所述目标地震褶积模型,K为指定的拓频实施的目标时间段内的总采样点数,h为服从零均值高斯分布的谱约束伪反射系数的方差的转置矩阵,h=[h1...hK]T包含了K个独立的参数,这代表每一个谱约束伪反射系数都服从零均值高斯分布,hk为第k个谱约束伪反射系数的方差,其中,K为选定时间段的总采样点数,k为大于0且小于K的自然数,σ2为噪音方差,p(d|m12)为所述地震数据振幅谱的似然函数,p(m1|h)为所述自动相关判别先验信息,p(d|h,σ2)为所述地震数据振幅谱的边缘似然函数,π为圆周率,Σ为协方差矩阵,Σ=(H+σ-2GTG)-1,μ=σ-2ΣGTd,H为服从零均值高斯分布的谱约束伪反射系数的方差的对角矩阵,H=diag(h1…hk)包含了k个独立的参数,这代表每一个谱约束伪反射系数都服从零均值高斯分布,G=F,F为δ函数在频率域的矩阵表达式,e为自然常数2.71828,μ为所述目标后验概率的均值,[·]T为矩阵的转置。
在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S134,对所述目标地震褶积模型的边缘似然函数使用最大期望算法,以获得所述转置矩阵的估计值和所述噪音方差的估计值;
在本实施例中,采用如下公式求解所述转置矩阵的估计值和所述噪音方差的估计值:
Figure GDA0003468676060000121
其中,hk为第k个谱约束伪反射系数的方差,其中,K为选定时间段的总采样点数,k为大于0且小于K的自然数,
Figure GDA0003468676060000122
为第k个目标后验概率的均值的平方,Σ为协方差矩阵,Σk为第k个协方差矩阵,k为大于0且小于K的自然数,σ2为噪音方差,d为所述目标地震褶积模型,G=F,F为δ函数在频率域的矩阵表达式,m1为谱约束伪反射系数序列,hm为第m个谱约束伪反射系数的方差,Σm为第m个协方差矩阵,m为取值范围从1到M的自然数,M为当前非零谱约束伪反射系数的个数。在本实施例中,选定时间段的总采样点数等于选定时窗内的子波总数。
步骤S135,根据所述转置矩阵的估计值、噪音方差的估计值和所述目标后验概率表达式,得到所述第一待反演伪反射系数估计值。
步骤S140,将所述第一待反演伪反射系数估计值代入所述目标后验概率表达式进行迭代计算,得到第二待反演伪反射系数估计值和当前迭代变量,并获得当前迭代次数,判断该当前迭代次数是否大于预设迭代阈值,以及判断当前迭代变量的值是否小于预设值。
当步骤S140执行结果为否时,即当所述当前迭代次数大于预设迭代阈值且所述当前迭代变量的值小于预设值时,将所述第二待反演伪反射系数估计值作为新的第一待反演伪反射系数估计值,并返回执行步骤S130。
在本实施例中,所述预设迭代阈值为可设置的大于0的自然数,所述当前迭代变量的表达式为:
Figure GDA0003468676060000123
其中,d为所述目标地震褶积模型,G=F,F为δ函数在频率域的矩阵表达式,
Figure GDA0003468676060000124
为从时间域的δ函数δ(t-tk)根据傅里叶变换得到的频率域结果,e为自然常数2.71828,i为虚数符号,π为圆周率,t为时间,tk为第k个谱约束伪反射系数的时间位置,k为大于0且小于选定时窗内的总采样点数的自然数,fm为待处理的第m个频段,其中,f为频率,m为取值范围从1到M的自然数,M为当前非零的谱约束伪反射系数的个数,m1为谱约束伪反射系数序列,
Figure GDA0003468676060000131
为二范数的平方。
当步骤S140执行结果为是时,即所述当前迭代次数等于预设迭代阈值或所述当前迭代变量的值小于预设值时,执行步骤S150,步骤S150为将获得的第二待反演伪反射系数估计值作为谱约束伪反射系数估计值。
在本实施例中,所述谱约束伪反射系数估计值用于得到地震数据的低频段、中频段和中高频段的所述谱约束伪反射系数频谱。由于在低频段、中频段和中高频段中,所述谱约束伪反射系数估计值和真实反射系数值具有很高的一致性,因此可以采用所述谱约束伪反射系数估计值推导出谱约束伪反射系数频谱,利用谱约束伪反射系数频谱的低频段、中频段和中高频段代替真实反射系数推导出的反射系数频谱的低频段、中频段和中高频段进行后续滤波处理。
步骤S160,根据所述谱约束伪反射系数估计值对获取到的地震数据进行处理,得到地震数据谱约束伪反射系数频谱,利用测井数据对所述地震数据谱约束伪反射系数频谱进行差异性滤波,以确定地震数据谱约束伪反射系数可信频谱范围。
在本实施例中,步骤S160中的地震数据为时变地震物理子波库中的子波,步骤S160即根据所述时变地震物理子波库中的子波和所述谱约束伪反射系数估计值得到谱约束伪反射系数频谱,步骤S160实现方法为:将所述谱约束伪反射系数估计值和做傅里叶变换后的所述时变地震物理子波库中的子波相乘,即可得到所述谱约束伪反射系数频谱。
在本实施例中,所述测井数据用于校验所述地震数据。
具体的,在本实施例中,上述步骤S160包括步骤S162至S168。
步骤S162,对测井数据中的纵波速度和密度数据进行方波化处理,以得到测井反射系数,并根据测井数据中的子波和所述测井反射系数,得到测井反射系数频谱。
在本实施例中,根据研究区域地质资料和情况,确定参考薄层厚度,选取研究区域内测井数据中纵波速度和密度数据,根据参考薄层厚度进行方波化,利用方波化后的纵波速度和密度数据计算得到测井反射系数。通过对所述测井反射系数和测井数据中的子波做傅里叶变换,即可得到所述测井反射系数频谱。
步骤S164,采用滑动频率域窗口对所述测井反射系数频谱和所述谱约束伪反射系数频谱进行计算,以得到所述测井反射系数频谱与所述谱约束伪反射系数频谱的相关系数。
步骤S166,根据所述相关系数在地震数据中低频和高频的突变点,采用四点梯形滤波算子选取待滤波的频谱范围,对所述待滤波的频谱范围作归一化处理,以得到归一化处理后的待滤波频谱范围,并根据所述相关系数与所述归一化处理后的待滤波频谱范围,得到最终频率域滤波器参数。
其中,本实施例使用四点梯形滤波算子对所述测井反射系数频谱进行筛选,将低频段相关系数的突变点投影在人工选定的所述测井反射系数频谱的低频截至范围内,得到的投影点即为四点梯形滤波算子的起始频率,将高频段相关系数的突变点投影在人工选定的所述测井反射系数频谱的高频截至范围内,得到的投影点即为四点梯形滤波算子的终止频率;将待归一化处理的测井反射系数频谱四点梯形滤波算子进行归一化处理,获得归一化的测井反射系数频谱四点梯形滤波算子。所述归一化的测井反射系数频谱四点梯形滤波算子与所述相关系数的乘积为最终频率域滤波器参数。
步骤S168,根据最终频率域滤波器参数对所述地震数据进行滤波,以得到所述地震数据谱约束伪反射系数可信频谱范围。
可以理解,根据谱约束伪反射系数的估计值和可信频段滤波,避免了由稀疏约束导致的反射系数估计值不准引起的问题。同时,经过差异性滤波的所述谱约束伪反射系数在低频段、中频段和中高频段的信息是准确的,因此可以根据谱约束反推地震数据中主频段外的信息,重点恢复地震数据中的低频段成分,使得地震数据倍频程得到有效提高,同时压制子波旁瓣,达到拓展低频段,补偿高频段的目的。
步骤S170,对所述谱约束伪反射系数可信频谱进行拓频处理,以获得经过拓频处理的拓频地震数据,其中,所述拓频地震数据的频带宽度大于拓频处理前的地震数据的频带宽度,从而根据拓频地震数据提高了地震分辨率。
在本实施例中,对谱约束伪反射系数的可信频谱进行拓频处理的步骤包括:对所述地震数据谱约束伪反射系数可信频谱范围进行反傅里叶变换。
同时,本发明实施例还提供一种电子设备,包括存储器、处理器以及存储在该存储器上并可在其处理器上运行的计算机程序,该计算机程序被处理器执行时,可实现如步骤S110至步骤S170中任一项的地震分辨率提高方法。
可以理解,通过拓频处理的地震数据能够有效消除子波影响,高质量地恢复低频段、中频段和中高频段的信息,而地震数据的频带宽度和质量是决定地震分辨率的重要因素,因此,拓频处理后的地震数据解决了现有技术中对地震数据的低频拓展能力低,拓频中无法兼顾信噪比和分辨率的同时提升,拓展地震数据不准确且噪音高的问题。
结合图2,给出了当地震数据中薄层厚度非常薄且无噪音的情况下,本发明地震分辨率提高方法的拓频能力和可信度。图2中的(a)表示时间域中拓频前后地震数据对比和谱约束反射系数(其中谱约束反射系数为图2中的参考反射系数),其横坐标Time表示时间,时间单位为ms,即毫秒,纵坐标Amplitude表示振幅;图2中的(b)表示频率域中拓频前后地震频谱对比,其横坐标Frequency表示频率,频率单位为Hz,即赫兹,纵坐标Amplitude表示振幅;图2(c)表示谱约束反射系数频谱与真实反射系数频谱对比,其横坐标Frequency表示频率,频率单位为Hz,即赫兹,纵坐标Amplitude表示振幅;图2(d)表示谱约束反射系数与真实反射系数对比,其横坐标Time表示时间,时间单位为ms,即毫秒,纵坐标Amplitude表示振幅。
其中,从图2中的(d)可以看到,谱约束反演伪反射系数并没有与真实反射系数完全一致,由于薄层效应的影响,获得的谱约束伪反射系数出现了偏差,但从图2中的(c)中可以清晰的看到,在地震主频段两边的一定区域内,反演获得的谱约束伪反射系数频谱和真实的反射系数频谱非常一致,因此图2中的(b)中经拓频后的地震数据频谱是可信的,频谱确实得到了有效的拓展,尤其是低频段非常好的得到了恢复。
结合图3,给出了当地震数据中薄层厚度非常薄且含噪音的情况下,地震分辨率提高方法的拓频能力和可信度。图3中的(a)表示含噪音的时间域中拓频前后地震数据对比和谱约束反射系数,其横坐标Time表示时间,时间单位为ms,即毫秒,纵坐标Amplitude表示振幅;图3中的(b)表示含噪音的频率域中拓频前后地震频谱对比,其横坐标Frequency表示频率,频率单位为Hz,即赫兹,纵坐标Amplitude表示振幅;图3中的(c)表示含噪音的谱约束反射系数频谱与真实反射系数频谱对比,其横坐标Frequency表示频率,频率单位为Hz,即赫兹,纵坐标Amplitude表示振幅;图3中的(d)表示含噪音的谱约束反射系数与真实反射系数对比,其横坐标Time表示时间,时间单位为ms,即毫秒,纵坐标Amplitude表示振幅。
从图3可以看出,噪音对谱约束伪反射系数的估计值产生了一定干扰,但反演结果仍然较佳,拓展后的频谱较好地还原了真实反射系数的频谱特征,在拓展频段内的频谱信息是可信性的。因此,处理后的地震数据不但提高了分辨率,还压制了噪音,提高了信噪比。
结合图4和图5,图4和图5为实际地震数据的拓频实例。图4的横坐标Trace表示地震记录,即本文所述地震数据,纵坐标Time表示时间,时间单位为ms,即毫秒;图5横坐标Frequency表示频率,频率单位为Hz,即赫兹,纵坐标Amplitude表示振幅,图5的图例中,Pre-process表示拓频前的地震数据频谱,After-process表示拓频后的地震数据频谱。
从图4中可以清晰地看出,拓频后的地震数据具有更高的分辨率,很多薄层得以显现,从图5中可以看到,地震数据频谱得到了有效的拓展,恢复了有益的信息,且剖面的随机噪音得以衰减,为后续解释和反演提供了坚实的基础。
综上所述,本发明所述地震分辨率提高方法,通过在频率域中预设谱约束伪反射系数,结合自动相关判别先验信息,能够获得准确的拓展结果和有效的低频信息,并结合实验数据进行说明,能够解决现有技术中对地震数据低频拓展能力低,拓频中无法兼顾信噪比和分辨率的同时提升,拓展地震数据不准确且噪音高的问题。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
虽然本发明所公开的实施方式如上,但所述的内容只是为了便于理解本发明而采用的实施方式,并非用以限定本发明。任何本发明所属技术领域内的技术人员,在不脱离本发明所公开的精神和范围的前提下,可以在实施的形式上及细节上作任何的修改与变化,但本发明的保护范围,仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种地震分辨率提高方法,包括:
步骤1,根据获取到的目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率表达式,并对目标地震褶积模型的边缘似然函数进行最大化处理,得到第一待反演伪反射系数估计值;
步骤2,将所述第一待反演伪反射系数估计值代入所述目标后验概率表达式进行迭代计算,得到第二待反演伪反射系数估计值和当前迭代变量,并获得当前迭代次数,判断该当前迭代次数是否大于预设迭代阈值,以及判断当前迭代变量的值是否小于预设值;
步骤3,当所述当前迭代次数大于预设迭代阈值且所述当前迭代变量的值小于预设值时,将所述第二待反演伪反射系数估计值作为新的第一待反演伪反射系数估计值,并返回执行步骤2,直至所述当前迭代次数等于预设迭代阈值或所述当前迭代变量的值小于预设值时,才将获得的第二待反演伪反射系数估计值作为谱约束伪反射系数估计值;
步骤4,根据所述谱约束伪反射系数估计值对获取到的地震数据进行处理,得到地震数据谱约束伪反射系数频谱,利用测井数据对所述地震数据谱约束伪反射系数频谱进行差异性滤波,以确定地震数据谱约束伪反射系数可信频谱范围;
步骤5,对所述地震数据谱约束伪反射系数可信频谱进行拓频处理,以获得拓频地震数据,其中,所述拓频地震数据的频带宽度大于拓频处理前的地震数据的频带宽度,从而提高地震分辨率。
2.根据权利要求1所述的地震分辨率提高方法,其特征在于,所述地震数据为叠后地震数据。
3.根据权利要求1所述的地震分辨率提高方法,其特征在于,在执行所述步骤1之前,所述方法还包括:
接收用户输入的目标时间段,并根据该目标时间段和地震数据中的子波构建时变地震物理子波库;
根据所述时变地震物理子波库和所述目标时间段,得到地震数据时间域表达式;
对所述地震数据时间域表达式进行傅里叶变换,以得到地震数据频率域表达式;
对所述地震数据频率域表达式在地震数据的指定频段内进行离散化处理,以得到频率域地震褶积模型;
根据所述频率域地震褶积模型和预设的谱约束伪反射系数,得到所述目标地震褶积模型;
当噪音服从零均值高斯分布时,根据贝叶斯原理和所述目标地震褶积模型得到该目标地震褶积模型的似然函数和目标地震褶积模型的边缘似然函数。
4.根据权利要求3所述的地震分辨率提高方法,其特征在于,根据该目标时间段和地震数据中的子波构建时变地震物理子波库的步骤包括:
在所述目标时间段内将随时间变化的非稳态地震数据作为稳态地震数据;
对所述稳态地震数据进行开窗处理,并在选定的时窗内采用井震标定的方法对子波进行提取;
对提取后的子波进行汇总,以得到所述时变地震物理子波库。
5.根据权利要求4所述的地震分辨率提高方法,其特征在于,所述步骤1包括:
在贝叶斯框架下,根据所述目标地震褶积模型、目标地震褶积模型的似然函数、目标地震褶积模型的边缘似然函数以及自动相关判别先验信息,得到目标后验概率,其中,所述目标后验概率表达式为:
p(m1|d,h,σ2)=p(d|m12)p(m1|h)/p(d|h,σ2),
其中,m1为谱约束伪反射系数序列,d为所述目标地震褶积模型,h为服从零均值高斯分布的谱约束伪反射系数的方差的转置矩阵,h=[h1 ... hK]T包含了K个独立的参数,hk为第k个谱约束伪反射系数的方差,K为选定时间段的总采样点数,且k为大于0且小于K的自然数,σ2为噪音方差,p(d|m12)为所述地震数据振幅谱的似然函数,p(m1|h)为所述自动相关判别先验信息,p(d|h,σ2)为所述地震数据振幅谱的边缘似然函数;
对所述目标地震褶积模型的边缘似然函数使用最大期望算法,以获得所述转置矩阵的估计值和所述噪音方差的估计值;
根据所述转置矩阵的估计值、噪音方差的估计值和所述目标后验概率表达式,得到所述第一待反演伪反射系数估计值。
6.根据权利要求3所述的地震分辨率提高方法,其特征在于,所述地震数据的指定频段为地震数据的主频段,所述地震数据的主频段为所述地震数据振幅谱峰值的0.707倍处的两个频率值之间的频率范围。
7.根据权利要求1所述的地震分辨率提高方法,其特征在于,所述当前迭代变量的表达式为:
Figure FDA0003468676050000031
其中,d为所述目标地震褶积模型,
Figure FDA0003468676050000032
Figure FDA0003468676050000033
为从时间域的δ函数δ(t-tk)根据傅里叶变换得到的频率域结果,e为自然常数2.71828,i为虚数符号,π为圆周率,t为时间,tk为第k个谱约束伪反射系数的时间位置,k为大于0且小于为选定时窗内的总采样点数的自然数,fm为待处理的第m个频段,其中,f为频率,m为大于0且小于M的自然数,M为当前非零的谱约束伪反射系数的个数,m1为谱约束伪反射系数序列。
8.根据权利要求1所述的地震分辨率提高方法,其特征在于,利用测井数据对所述地震数据谱约束伪反射系数频谱进行差异性滤波,以确定地震数据谱约束伪反射系数可信频谱范围的步骤包括:
对测井数据中的纵波速度和密度数据进行方波化处理,以得到测井反射系数,并根据测井数据中的子波和所述测井反射系数,得到测井反射系数频谱;
采用滑动频率域窗口对所述测井反射系数频谱和所述地震数据谱约束伪反射系数频谱进行计算,以得到所述测井反射系数频谱与所述地震数据谱约束伪反射系数频谱的相关系数;
根据所述相关系数在地震数据中低频和高频的突变点,采用四点梯形滤波算子选取待滤波的频谱范围,对所述待滤波的频谱范围作归一化处理,以得到归一化处理后的待滤波频谱范围,并根据所述相关系数与所述归一化处理后的待滤波频谱范围,得到最终频率域滤波器参数;
根据最终频率域滤波器参数对所述地震数据进行滤波,以得到所述地震数据谱约束伪反射系数可信频谱范围。
9.根据权利要求1所述的地震分辨率提高方法,其特征在于,对所述地震数据谱约束伪反射系数可信频谱进行拓频处理的步骤包括:对所述地震数据谱约束伪反射系数可信频谱范围进行反傅里叶变换。
10.一种电子设备,包括存储器、处理器以及存储在该存储器上并可在其处理器上运行的计算机程序,该计算机程序被处理器执行时,实现如权利要求1至9中任一项所述的地震分辨率提高方法。
CN201910631191.4A 2019-07-12 2019-07-12 一种地震分辨率提高方法及电子设备 Active CN112213773B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910631191.4A CN112213773B (zh) 2019-07-12 2019-07-12 一种地震分辨率提高方法及电子设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910631191.4A CN112213773B (zh) 2019-07-12 2019-07-12 一种地震分辨率提高方法及电子设备

Publications (2)

Publication Number Publication Date
CN112213773A CN112213773A (zh) 2021-01-12
CN112213773B true CN112213773B (zh) 2022-04-22

Family

ID=74048602

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910631191.4A Active CN112213773B (zh) 2019-07-12 2019-07-12 一种地震分辨率提高方法及电子设备

Country Status (1)

Country Link
CN (1) CN112213773B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114002736B (zh) * 2021-09-07 2023-07-07 中国矿业大学 一种基于权重反褶积的地震勘探多频数据融合方法
CN117434604A (zh) * 2023-02-20 2024-01-23 中国石油化工股份有限公司 地震数据处理方法、装置、存储介质及电子设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104280765A (zh) * 2013-07-11 2015-01-14 中国石油化工股份有限公司 基于变子波反射系数反演的地震高分辨处理方法
CN109425889A (zh) * 2017-08-22 2019-03-05 中国石油化工股份有限公司 一种用于刻画古岩溶暗河的方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4747054A (en) * 1984-11-28 1988-05-24 Conoco Inc. Method for non-linear signal matching
US5416750A (en) * 1994-03-25 1995-05-16 Western Atlas International, Inc. Bayesian sequential indicator simulation of lithology from seismic data
US9207356B2 (en) * 2013-07-29 2015-12-08 Chevron U.S.A. Inc. System and method for estimating a reservoir parameter using joint stochastic inversion of multisource geophysical data
CN104181589A (zh) * 2014-08-20 2014-12-03 中国石油集团川庆钻探工程有限公司地球物理勘探公司 一种非线性反褶积方法
CN105467442B (zh) * 2015-12-09 2017-10-03 中国石油大学(北京) 全局优化的时变稀疏反褶积方法及装置
CN105842732B (zh) * 2016-03-16 2018-03-23 中国石油大学(北京) 多道稀疏反射系数的反演方法及系统
CN108663711B (zh) * 2018-04-04 2019-10-01 电子科技大学 一种基于τ分布的贝叶斯地震反演方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104280765A (zh) * 2013-07-11 2015-01-14 中国石油化工股份有限公司 基于变子波反射系数反演的地震高分辨处理方法
CN109425889A (zh) * 2017-08-22 2019-03-05 中国石油化工股份有限公司 一种用于刻画古岩溶暗河的方法

Also Published As

Publication number Publication date
CN112213773A (zh) 2021-01-12

Similar Documents

Publication Publication Date Title
Liu et al. A 1D time-varying median filter for seismic random, spike-like noise elimination
US9268047B2 (en) Geophysical surveying
Gómez et al. A simple method inspired by empirical mode decomposition for denoising seismic data
Wu et al. Noise attenuation for 2-D seismic data by radial-trace time-frequency peak filtering
CN109254324B (zh) 全频保幅地震数据处理方法和装置
CN110187388B (zh) 一种基于变分模态分解的稳定地震品质因子q估计方法
US10895654B2 (en) Method for generating optimized seismic target spectrum
CN112213773B (zh) 一种地震分辨率提高方法及电子设备
CN107179550B (zh) 一种数据驱动的地震信号零相位反褶积方法
Yong-Shou et al. A time-varying wavelet extraction using local similarity
Al-Dossary et al. Lineament-preserving filtering
CN112882099B (zh) 一种地震频带拓宽方法、装置、介质及电子设备
CN112835103A (zh) 自适应去鬼波与宽频准零相位反褶积联合处理方法及系统
CN111399057A (zh) 一种基于非凸稀疏约束的地震资料噪声压制方法
Alsalmi et al. Mask filtering to the Wigner-Ville distribution
CN112817040B (zh) 宽频准零相位反褶积处理方法、装置、电子设备及介质
CN112099083A (zh) 一种基于双谱谱比对数的品质因子估计方法及系统
Sajid et al. Nonstationary differential resolution: An algorithm to improve seismic resolution
CN110568491B (zh) 一种品质因子q的估算方法
CN110749923A (zh) 一种基于范数方程提高分辨率的反褶积方法
CN112764108B (zh) 一种基于改进经验小波变换的新型地震资料噪声压制算法
CN112731526A (zh) 依据地震衰减截距检测油气储层的方法
CN117849875B (zh) 一种地震信号分析方法、系统、装置及存储介质
CN113740909B (zh) 一种基于稀疏s变换和自适应对数谱比法的地震衰减估计方法、系统、设备及存储介质
CN113419276B (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
GR01 Patent grant
GR01 Patent grant