CN108363739B - 一种基于稀疏采集的地震资料高低频拓展方法 - Google Patents

一种基于稀疏采集的地震资料高低频拓展方法 Download PDF

Info

Publication number
CN108363739B
CN108363739B CN201810057585.9A CN201810057585A CN108363739B CN 108363739 B CN108363739 B CN 108363739B CN 201810057585 A CN201810057585 A CN 201810057585A CN 108363739 B CN108363739 B CN 108363739B
Authority
CN
China
Prior art keywords
frequency
low
data
detector
spectral density
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
CN201810057585.9A
Other languages
English (en)
Other versions
CN108363739A (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.)
Ocean University of China
Original Assignee
Ocean University of China
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 Ocean University of China filed Critical Ocean University of China
Priority to CN201810057585.9A priority Critical patent/CN108363739B/zh
Publication of CN108363739A publication Critical patent/CN108363739A/zh
Application granted granted Critical
Publication of CN108363739B publication Critical patent/CN108363739B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases

Landscapes

  • Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Data Mining & Analysis (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于稀疏采集的地震资料高低频拓展方法，具体为：稀疏布设低频检波器和高频检波器，获得常规、低频以及高频地震数据；进行匹配，获得匹配后的低频以及高频地震数据；计算同一低频与常规数据功率谱密度比以及高频与常规数据功率谱密度比；分析功率谱密度比曲线的稳定范围，确定可拓展的最低以及最高有效频率；设计频率域拓频算子；求取时域拓频算子；计算所有稀疏点位的时域拓频算子；进行空间插值获得工区所有点位时域拓频算子；将常规地震数据与插值后的时域拓频算子进行褶积，获得拓频后的地震数据，实现常规地震数据高低频拓展。本发明利用稀疏布设获得频带丰富的原始地震资料，通过后续拓频处理实现高分辨地震勘探。

Description

一种基于稀疏采集的地震资料高低频拓展方法
技术领域
本发明属于地震资料采集处理技术领域，具体地说，涉及一种基于稀疏采集的地震资料高低频拓展方法。
背景技术
我国开发地震勘探技术总体上处于发展阶段，国外开发地震勘探技术起步早、探索较多，逐渐成为领先技术高地，但对外技术封锁；国内也开展了大量的探索研究和区块实践，仍然是在传统三维地震基础上缩小面元、增加覆盖次数、增大接收道数的高精度技术系列。
提高高精度勘探纵向分辨率的核心是获得频带宽、主频的的信号，既要高频成分丰富，也要低频信息，那种缺失低频成分的高分辨率是假高分辨率。国外研究成功了低频震源和低频检波器，与高频检波器结合，实现真正意义含有低频信息的宽高频的高分辨率勘探，但是，一是这些设备昂贵，二是国外技术封锁。我国在这方面做了一些探索，尚没有形成完整产品及技术。
现有技术面临的难题：第一、常规检波器低频信息缺失，富含低频信息的宽频信号获取困难；第二、能接收低频信息的低频传感器价格昂贵；第三、单一接收方式对高频的保护不够。
发明内容
有鉴于此，本发明提供了一种基于稀疏采集的地震资料高低频拓展方法，本发明利用稀疏布设低频检波器与高频检波器获得频带丰富的原始地震资料，补充常规地震资料的不足，通过后续拓频处理实现高分辨地震勘探。
为了解决上述技术问题，本发明公开了一种基于稀疏采集的地震资料高低频拓展方法，包括以下步骤：
步骤1、在野外常规三维地震采集的基础上，稀疏布设低频检波器以及高频检波器，同时采集地震信号，获得常规地震数据、低频地震数据以及高频地震数据；
步骤2、利用步骤1获得的低频地震数据以及高频地震数据向常规地震数据进行匹配，包括互相关走时匹配以及振幅匹配，从而获得匹配之后的低频地震数据以及高频地震数据；
步骤3、利用步骤2得到的匹配后数据分别计算同一点位不同炮记录的低频数据与常规数据功率谱密度比PSDRdf以及高频数据与常规数据功率谱密度比PSDRgf；分析同一点位不同炮两种功率谱密度比曲线的稳定范围作为常规地震资料频带可拓展范围，确定可拓展的最低有效频率flow以及最高有效频率fhigh；
步骤4、利用步骤3计算获得的低频数据与常规数据功率谱密度比 PSDRdf、高频数据与常规数据功率谱密度比PSDRgf以及常规地震资料频带最低有效频率flow以及最高有效频率fhigh，设计频率域拓频算子；
步骤5、利用步骤4获得的频率域拓频算子H(k)求取时域拓频算子；
步骤6、重复步骤2-步骤5的方法计算所有稀疏点位的时域拓频算子；
步骤7、利用步骤6获得的时域拓频算子进行空间插值获得工区所有点位时域拓频算子；
步骤8、将步骤1常规地震数据f(t)与步骤7中插值后对应点位的时域拓频算子h(t)进行褶积，获得拓频后的地震数据fdg(t)，从而实现常规地震数据高低频拓展。
进一步地，所述稀疏布设低频检波器以及高频检波器具体为：每 0.5-1km2布设一组低频检波器及高频检波器，与常规检波器埋置于同一点位，间距0.5m以内，保证三种检波器接收端信号相同。
进一步地，所述步骤2中的互相关走时匹配方法具体为：
设同一炮同一道的低频检波器数据为d(t)，常规检波器数据为f(t)，计算两者的互相关函数Rdf(τ)：
Figure BDA0001554282530000031
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围 [0，M-1]；
求取Rdf(τ)最大值对应的时刻τdf，进而对低频检波器数据d(t)进行匹配，互相关走时匹配后的低频数据ds(t)：
ds(t)＝d(t+τdf) (2)
其中，τdf是自相关函数Rdf(τ)最大值对应的时刻，t为采样时间，范围 [0，M-1]；
同理，设高频检波器数据为g(t)，则高频检波器与常规检波器的互相关函数Rgf(τ)为：
Figure BDA0001554282530000032
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围 [0，M-1]；
那么互相关走时匹配后的高频数据为：
dg(t)＝g(t+τgf)
其中，τgf是自相关函数Rgf(τ)最大值对应的时刻，t为采样时间，范围 [0，M-1]。
进一步地，所述步骤2中的振幅匹配方法具体为：
设互相关走时匹配之后的低频检波器数据为ds(t)，常规检波器数据为 f(t)，其中t为采样时间，范围[0，M-1]；求取ds(t)振幅最大值max(ds(t))，以及f(t)振幅最大值max(f(t))，那么低频检波器振幅匹配结果dsf(t)有：
Figure BDA0001554282530000033
同理，设互相关走时匹配之后的高频检波器数据为dg(t)，高频检波器振幅匹配结果dgf(t)为：。
Figure BDA0001554282530000041
其中，max()表示求取函数最大值。
进一步地，所述步骤3中功率谱密度(Power spectral density，PSD)计算方法采用周期图法，设x(n)是一个长度为M的有限长序列，则定义x(n)的 N点离散傅里叶变换为
Figure BDA0001554282530000042
式中，k＝0，1，…，N-1，
Figure BDA0001554282530000043
j2＝-1，N是变换区间长度，N≥M；
其功率谱密度为：
Figure BDA0001554282530000044
那么，设步骤2获得匹配后的低频地震数据为dsf(t)，其功率谱密度为：
Figure BDA0001554282530000045
设步骤2获得匹配后的低频地震数据为gsf(t)，其功率谱密度为：
Figure BDA0001554282530000046
设常规地震数据为f(t)，其功率谱密度为：
Figure BDA0001554282530000047
其中，j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M；
那么，低频数据与常规数据功率谱密度比为：
Figure BDA0001554282530000051
高频数据与常规数据功率谱密度比为：
Figure BDA0001554282530000052
进一步地，所述分析同一点位不同炮两种功率谱密度比曲线的稳定范围作为常规地震资料频带可拓展范围，确定可拓展的最低有效频率flow以及最高有效频率fhigh具体为：功率谱密度比反应了两种检波器的频率响应特征，通过分析功率谱密度比曲线获得常规检波器数据的频率拓展范围，功率谱密度比曲线包括低频检波器的功率谱密度比曲线和高频检波器的功率谱密度曲线，其中，低频检波器的功率谱密度比曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿降频方向呈曲线提升，表示在此频段范围低频检波器频率响应好于常规检波器，低频检波器记录到常规检波器中被压制的低频成分，此区域作为常规检波器可拓展的低频范围，区域2的低频端边界称为最低有效频率 flow；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展；高频检波器的功率谱密度曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿升频方向呈曲线提升，表示在此频段范围高频检波器频率响应好于常规检波器，高频检波器记录到常规检波器中被压制的高频成分，此区域作为常规检波器可拓展的高频范围，区域2的高频端边界称为最高有效频率fhigh；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展。
进一步地，步骤4中的频率域拓频算子计算公式为：
Figure BDA0001554282530000061
其中，N是同一点位选取不同炮记录的数量，N≥10；
Figure BDA0001554282530000064
是第N 道低频数据与常规数据功率谱密度比；PSDRgfn(k)是第n道高频数据与常规数据功率谱密度比，n范围为[1，N]；flow是最低有效频率，即低频检波器的功率谱密度比曲线区域2的左边界；f1是低频检波器的功率谱密度比曲线区域2 的右边界；fhigh是最低有效频率，即高频检波器的功率谱密度曲线区域2的右边界；f2是高频检波器的功率谱密度曲线区域2的左边界。
进一步地，步骤5中的时域拓频算子h(t)表示为：
Figure BDA0001554282530000062
其中j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M。
进一步地，步骤7中的利用步骤6获得的时域拓频算子进行空间插值获得工区所有点位时域拓频算子具体为：
设步骤6获得的稀疏点位时域拓频算子为hi(t)，i是稀疏点位个数，范围[1，N]；xi、yi分别是第i个点的横纵坐标；h0(t)是所求点位的时域拓频算子， x0、y0分别是所求点位横纵坐标；
那么已知拓频算子点位与未知点位的距离：
Figure BDA0001554282530000063
每个已知点的权重系数：
Figure BDA0001554282530000071
未知点位时域拓频算子空间差值结果：
Figure BDA0001554282530000072
根据求取的所有未知点位的时域拓频算子，即获得区所有点位时域拓频算子。
进一步地，步骤8中的实现常规地震数据高低频拓展的过程如下：
Figure BDA0001554282530000073
其中t为样点数，范围[0，M-1]。
与现有技术相比，本发明可以获得包括以下技术效果：
1)本发明借助稀疏的灵敏度高的低频检波器获取低频信号，和常规检波器记录进行调谐，开发提取低频补偿因子技术和实现补偿的技术，以确保常规检波器记录恢复低频信息；研究获取足够高的高频信息，与常规检波器记录进行调谐，获得高频补偿因子，进而实现拓展高频，本发明使用的高性能检波器(包括reftek125a低频检波器、陆用压电高频检波器)，能将常规检波器资料的频宽由8Hz-100Hz提升到4.5Hz-130Hz，实现常规地震资料的高低频拓展。
2)本发明使用的高性能检波器能够获得频带丰富的地震信号，但是其价格昂贵不适用于大规模生产，本发明在常规检波器勘探的基础之上，每 0.5km2-1km2选择一个接受点位，同时布设高性能检波器，实现稀疏采集，从而减少的高性能检波器的使用数量，节约勘探成本。
3)由于稀疏布设高性能检波器，只有较少的点位能够获得高低频补偿因子，因此通过空间插值方法，将若干点位的补偿因子推广到工区全部点位，从而实现整个工区常规检波器资料的高低频拓展。
当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是本发明地震仪造价与频谱响应响特征关系图；
图2是本发明步骤1工区高性能检波器稀疏布设点位图；
图3是本发明步骤1常规检波器、reftek125a低频检波器以及陆用压电检波器的频谱分析图；
图4是本发明步骤2低频检波器互相关走时匹配与振幅匹配前后效果图；其中，a代表常规检波器，b代表reftek125a低频检波器匹配前的记录，c代表低频检波器走时匹配后的记录；
图5是本发明步骤3同一点位10炮记录的低频检波器数据与常规检波器数据功率谱密度比曲线图；
图6是本发明步骤3同一点位10炮记录的高频检波器数据与常规检波器数据功率谱密度比曲线图；
图7是本发明步骤4频率域拓频算子图；
图8是本发明步骤5时域拓频算子图；
图9是本发明步骤7原始拓频算子与差值后拓频算子对比图；其中，左图第8、9道是通过步骤2-5计算的时域拓频算子，右图第8、9道是通过步骤7计算的时域拓频算子；
图10是本发明步骤8常规检波器资料高低频拓展前叠加剖面图；
图11是本发明步骤8常规检波器资料高低频拓展后叠加剖面图。
具体实施方式
以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。
本发明公开了一种基于稀疏采集的地震资料高低频拓展方法，包括以下步骤：
步骤1、在野外常规三维地震采集的基础上，稀疏布设低频检波器以及高频检波器，同时采集地震信号，获得常规地震数据、低频地震数据以及高频地震数据；稀疏布设方案为每0.5-1km2布设一组低频检波器及高频检波器，与常规检波器埋置于同一点位，间距0.5m以内，保证三种检波器接收端信号相同；
如图1所示，本发明每0.5km2-1km2选择一个接受点位，同时布设高性能检波器，实现稀疏采集，从而减少的高性能检波器的使用数量，节约勘探成本；图2是工区高性能检波器稀疏布设点位图，图钉标志布设点位；按照步骤1布设常规检波器以及高性能检波器，采集的地震资料频谱如图3所示， reftek125a低频检波器的低频成分(1Hz-10Hz)比常规检波器丰富，陆用压电的高频成分(30Hz-100Hz)比常规检波器丰富，高性能检波器能够采集常规检波器缺失的频率成分。
步骤2、利用步骤1获得的低频地震数据以及高频地震数据向常规地震数据进行匹配，包括互相关走时匹配以及振幅匹配，从而获得匹配之后的低频地震数据以及高频地震数据；
其中，由于检波器记录方式不同，同一点位不同检波器采集数据在振幅和走时上有差异，如图4低频检波器走时匹配前后效果图，图4a是常规检波器数据，图4b是相同炮reftek125a低频检波器数据，两者振幅走时差异较大，因此需要将低频检波器数据以及高频检波器数据向常规检波器数据匹配，包括互相关走时匹配和振幅匹配先进行走时匹配再进行振幅匹配。
互相关走时匹配方法：以低频检波器数据匹配为例，高频检波器数据匹配与之相同。设同一炮同一道的低频检波器数据为d(t)，常规检波器数据为 f(t)，计算两者的互相关函数Rdf(τ)：
Figure BDA0001554282530000091
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围 [0，M-1]。
求取Rdf(τ)最大值对应的时刻τdf，进而对低频检波器数据d(t)进行匹配，互相关走时匹配后的低频数据ds(t)：
ds(t)＝d(t+τdf) (2)
其中，τdf是自相关函数Rdf(τ)最大值对应的时刻，t为采样时间，范围 [0，M-1]。
同理，设高频检波器数据为g(t)，则高频检波器与常规检波器的互相关函数Rgf(τ)为：
Figure BDA0001554282530000101
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围 [0，M-1]。
那么互相关走时匹配后的高频数据为：
dg(t)＝g(t+τgf)
其中，τgf是自相关函数Rgf(τ)最大值对应的时刻，t为采样时间，范围 [0，M-1]。
振幅匹配方法：以低频检波器数据匹配为例，高频检波器数据匹配与之相同。
设互相关走时匹配之后的低频检波器数据为ds(t)，常规检波器数据为 f(t)，其中t为采样时间，范围[0，M-1]。求取ds(t)振幅最大值max(ds(t))，以及f(t)振幅最大值max(f(t))，那么低频检波器数据振幅匹配结果dsf(t)有：
Figure BDA0001554282530000102
同理，设互相关走时匹配之后的高频检波器数据为dg(t)，高频检波器振幅匹配结果dgf(t)为：。
Figure BDA0001554282530000111
其中，max()表示求取函数最大值。
图4c是reftek125a低频检波器数据进行互相关走时匹配和振幅匹配后的结果，与图4a常规检波器数据对比，低频检波器数据振幅与走时与常规检波器一致。同理，若设高频检波器数据为g(t)，匹配后的数据为gsf(t)。
步骤3、利用步骤2得到的匹配后数据分别计算同一点位不同炮记录的低频数据与常规数据功率谱密度比PSDRdf以及高频数据与常规数据功率谱密度比PSDRgf；分析同一点位不同炮两种功率谱密度比曲线的稳定范围作为常规地震资料频带可拓展范围，确定可拓展的最低有效频率flow以及最高有效频率fhigh。
其中，本发明通过功率谱密度比(Power spectral density ratio，PSDR) 分析不同检波器之间的频谱响应；功率谱密度(Power spectral density，PSD) 计算方法采用周期图法，设x(n)是一个长度为M的有限长序列，则定义x(n) 的N点离散傅里叶变换为
Figure BDA0001554282530000112
式中，k＝0，1，…，N-1，
Figure BDA0001554282530000113
j2＝-1，N是变换区间长度，N≥M。
其功率谱密度为：
Figure BDA0001554282530000114
那么，设步骤2获得匹配后的低频地震数据为dsf(t)，其功率谱密度为：
Figure BDA0001554282530000121
设步骤2获得匹配后的低频地震数据为gsf(t)，其功率谱密度为：
Figure BDA0001554282530000122
设常规地震数据为f(t)，其功率谱密度为：
Figure BDA0001554282530000123
其中，j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M。
那么，低频数据与常规数据功率谱密度比为：
Figure BDA0001554282530000124
高频数据与常规数据功率谱密度比为：
Figure BDA0001554282530000125
功率谱密度比反应了两种检波器的频率响应特征，通过分析功率谱密度比曲线获得常规检波器数据的频率拓展范围。图5是同一点位10炮记录的低频检波器数据与常规检波器数据功率谱密度比曲线。可将曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿降频方向呈曲线提升，表示在此频段范围低频检波器频率响应好于常规检波器，低频检波器记录到常规检波器中被压制的低频成分，此区域作为常规检波器可拓展的低频范围，区域2的低频端边界称为最低有效频率flow；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展。图6是同一点位10炮记录的高频检波器数据与常规检波器数据功率谱密度比曲线。可将曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿升频方向呈曲线提升，表示在此频段范围高频检波器频率响应好于常规检波器，高频检波器记录到常规检波器中被压制的高频成分，此区域作为常规检波器可拓展的高频范围，区域2的高频端边界称为最高有效频率fhigh；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展。
步骤4、利用步骤3计算获得的低频数据与常规数据功率谱密度比
Figure BDA0001554282530000132
高频数据与常规数据功率谱密度比PSDRgf以及常规地震资料频带最低有效频率flow以及最高有效频率fhigh，设计频率域拓频算子：
频率域拓频算子计算公式为：
Figure BDA0001554282530000131
其中，N是同一点位选取不同炮记录的数量，N≥10；
Figure BDA0001554282530000133
是第N 道低频数据与常规数据功率谱密度比；PSDRgfn(k)是第n道高频数据与常规数据功率谱密度比，n范围为[1，N]；flow是最低有效频率，即图5区域2的左边界；f1是图5区域2的右边界；fhigh是最低有效频率，即图6区域2的右边界；f2是图6区域2的左边界；图7是应用以上方法设计的频率域拓频算子。
步骤5、利用步骤4获得的频率域拓频算子H(k)求取时域拓频算子，本方法通过傅里叶逆变换实现；时域拓频算子h(t)表示为：
Figure BDA0001554282530000141
其中j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M。图8是由图7频率域拓频算子H(k)求取的时域拓频算子h(t)。
步骤6、重复步骤2-步骤5的方法计算所有稀疏点位的时域拓频算子。
步骤7、利用步骤6获得的时域拓频算子进行空间插值获得工区所有点位时域拓频算子。由于稀疏布设低频以及高频检波器，只有较少的点位能够获得对应的匹配因子，对于没有布设低、高频检波器的点位就需要通过空间插值技术实现算子求取。
本方法采用反距离权重差值：
设步骤6获得的稀疏点位时域拓频算子为hi(t)，i是稀疏点位个数，范围[1，N]；xi、yi分别是第i个点的横纵坐标；h0(t)是所求点位的时域拓频算子， x0、y0分别是所求点位横纵坐标。
那么已知拓频算子点位与未知点位的距离：
Figure BDA0001554282530000142
每个已知点的权重系数：
Figure BDA0001554282530000143
未知点位时域拓频算子空间差值结果：
Figure BDA0001554282530000151
按步骤7求取所有未知点位的时域拓频算子，即获得区所有点位时域拓频算子。图9是原始拓频算子与差值后拓频算子对比，左图第8、9道是通过步骤2-5计算的时域拓频算子，右图第8、9道是通过步骤7计算的时域拓频算子，两者计算结果几乎相同，证明该空间差值方法的有效性。
步骤8、将步骤1常规地震数据f(t)与步骤7中插值后对应点位的时域拓频算子h(t)进行褶积，获得拓频后的地震数据fdg(t)，从而实现常规地震数据高低频拓展。实现过程如下：
Figure BDA0001554282530000152
其中t为样点数，范围[0，M-1]。
图10是常规检波器资料高低频拓展前叠加剖面，图11是常规检波器资料高低频拓展后叠加剖面，从图中可以看出高低频拓展后，剖面浅层分辨率提高，深层同相轴更加连续，证明本发明的有效性。
上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。

Claims (10)

1.一种基于稀疏采集的地震资料高低频拓展方法，其特征在于，包括以下步骤：
步骤1、在野外常规三维地震采集的基础上，稀疏布设低频检波器以及高频检波器，同时采集地震信号，获得常规地震数据、低频地震数据以及高频地震数据；
步骤2、利用步骤1获得的低频地震数据以及高频地震数据向常规地震数据进行匹配，包括互相关走时匹配以及振幅匹配，从而获得匹配之后的低频地震数据以及高频地震数据；
步骤3、利用步骤2得到的匹配后数据分别计算同一点位不同炮记录的低频数据与常规数据功率谱密度比PSDRdf以及高频数据与常规数据功率谱密度比PSDRgf；分析同一点位不同炮两种功率谱密度比曲线的稳定范围作为常规地震资料频带可拓展范围，确定可拓展的最低有效频率flow以及最高有效频率fhigh；
步骤4、利用步骤3计算获得的低频数据与常规数据功率谱密度比PSDRdf、高频数据与常规数据功率谱密度比PSDRgf以及常规地震资料频带最低有效频率flow以及最高有效频率fhigh，设计频率域拓频算子；
步骤5、利用步骤4获得的频率域拓频算子H(k)求取时域拓频算子；
步骤6、重复步骤2-步骤5的方法计算所有稀疏点位的时域拓频算子；
步骤7、利用步骤6获得的时域拓频算子进行空间插值获得工区所有点位时域拓频算子；
步骤8、将步骤1常规地震数据f(t)与步骤7中插值后对应点位的时域拓频算子h(t)进行褶积，获得拓频后的地震数据fdg(t)，从而实现常规地震数据高低频拓展。
2.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，所述稀疏布设低频检波器以及高频检波器具体为：每0.5-1km2布设一组低频检波器及高频检波器，与常规检波器埋置于同一点位，间距0.5m以内，保证三种检波器接收端信号相同。
3.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，所述步骤2中的互相关走时匹配方法具体为：
设同一炮同一道的低频检波器数据为d(t)，常规检波器数据为f(t)，计算两者的互相关函数Rdf(τ)：
Figure FDA0001554282520000021
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围[0，M-1]；
求取Rdf(τ)最大值对应的时刻τdf，进而对低频检波器数据d(t)进行匹配，互相关走时匹配后的低频数据ds(t)：
ds(t)＝d(t+τdf) (2)其中，τdf是自相关函数Rdf(τ)最大值对应的时刻，t为采样时间，范围[0，M-1]；
同理，设高频检波器数据为g(t)，则高频检波器与常规检波器的互相关函数Rgf(τ)为：
Figure FDA0001554282520000022
其中t为采样时间，范围[0，M-1]，τ是互相关函数自变量，范围[0，M-1]；
那么互相关走时匹配后的高频数据为：
dg(t)＝g(t+τgf)
其中，τgf是自相关函数Rgf(τ)最大值对应的时刻，t为采样时间，范围[0，M-1]。
4.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，所述步骤2中的振幅匹配方法具体为：
设互相关走时匹配之后的低频检波器数据为ds(t)，常规检波器数据为f(t)，其中t为采样时间，范围[0，M-1]；求取ds(t)振幅最大值max(ds(t))，以及f(t)振幅最大值max(f(t))，那么低频检波器振幅匹配结果dsf(t)有：
Figure FDA0001554282520000031
同理，设互相关走时匹配之后的高频检波器数据为dg(t)，高频检波器振幅匹配结果dgf(t)为：
Figure FDA0001554282520000032
其中，max()表示求取函数最大值。
5.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，所述步骤3中功率谱密度(Power spectral density，PSD)计算方法采用周期图法，设x(n)是一个长度为M的有限长序列，则定义x(n)的N点离散傅里叶变换为
Figure FDA0001554282520000033
式中，k＝0，1，…，N-1，
Figure FDA0001554282520000034
j2＝-1，N是变换区间长度，N≥M；
其功率谱密度为：
Figure FDA0001554282520000035
那么，设步骤2获得匹配后的低频地震数据为dsf(t)，其功率谱密度为：
Figure FDA0001554282520000036
设步骤2获得匹配后的低频地震数据为gsf(t)，其功率谱密度为：
Figure FDA0001554282520000037
设常规地震数据为f(t)，其功率谱密度为：
Figure FDA0001554282520000041
其中，j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M；
那么，低频数据与常规数据功率谱密度比为：
Figure FDA0001554282520000042
高频数据与常规数据功率谱密度比为：
Figure FDA0001554282520000043
6.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，所述分析同一点位不同炮两种功率谱密度比曲线的稳定范围作为常规地震资料频带可拓展范围，确定可拓展的最低有效频率flow以及最高有效频率fhigh具体为：功率谱密度比反应了两种检波器的频率响应特征，通过分析功率谱密度比曲线获得常规检波器数据的频率拓展范围，功率谱密度比曲线包括低频检波器的功率谱密度比曲线和高频检波器的功率谱密度曲线，其中，低频检波器的功率谱密度比曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿降频方向呈曲线提升，表示在此频段范围低频检波器频率响应好于常规检波器，低频检波器记录到常规检波器中被压制的低频成分，此区域作为常规检波器可拓展的低频范围，区域2的低频端边界称为最低有效频率flow；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展；高频检波器的功率谱密度曲线分为三个区域，区域1功率谱密度比曲线平稳，表示在此频段两种检波器频率响应相同；区域2功率谱密度比曲线沿升频方向呈曲线提升，表示在此频段范围高频检波器频率响应好于常规检波器，高频检波器记录到常规检波器中被压制的高频成分，此区域作为常规检波器可拓展的高频范围，区域2的高频端边界称为最高有效频率fhigh；区域3功率谱密度比曲线呈无规则变化，表示该区域两种检波器自噪声水平超过有效信号，此频段常规检波器数据的有效信号不能拓展。
7.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，步骤4中的频率域拓频算子计算公式为：
Figure FDA0001554282520000051
其中，N是同一点位选取不同炮记录的数量，N≥10；PSDRafN(k)是第N道低频数据与常规数据功率谱密度比；PSDRgfn(k)是第n道高频数据与常规数据功率谱密度比，n范围为[1，N]；flow是最低有效频率，即低频检波器的功率谱密度比曲线区域2的左边界；f1是低频检波器的功率谱密度比曲线区域2的右边界；fhigh是最低有效频率，即高频检波器的功率谱密度曲线区域2的右边界；f2是高频检波器的功率谱密度曲线区域2的左边界。
8.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，步骤5中的时域拓频算子h(t)表示为：
Figure FDA0001554282520000052
其中j2＝-1，t为样点数，范围[0，M-1]，N是变换区间长度，N≥M。
9.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，步骤7中的利用步骤6获得的时域拓频算子进行空间插值获得工区所有点位时域拓频算子具体为：
设步骤6获得的稀疏点位时域拓频算子为hi(t)，i是稀疏点位个数，范围[1，N]；xi、yi分别是第i个点的横纵坐标；h0(t)是所求点位的时域拓频算子，x0、y0分别是所求点位横纵坐标；
那么已知拓频算子点位与未知点位的距离：
Figure FDA0001554282520000061
每个已知点的权重系数：
Figure FDA0001554282520000062
未知点位时域拓频算子空间差值结果：
Figure FDA0001554282520000063
根据求取的所有未知点位的时域拓频算子，即获得区所有点位时域拓频算子。
10.根据权利要求1所述的基于稀疏采集的地震资料高低频拓展方法，其特征在于，步骤8中的实现常规地震数据高低频拓展的过程如下：
Figure FDA0001554282520000064
其中t为样点数，范围[0，M-1]。
CN201810057585.9A 2018-01-22 2018-01-22 一种基于稀疏采集的地震资料高低频拓展方法 Active CN108363739B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810057585.9A CN108363739B (zh) 2018-01-22 2018-01-22 一种基于稀疏采集的地震资料高低频拓展方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810057585.9A CN108363739B (zh) 2018-01-22 2018-01-22 一种基于稀疏采集的地震资料高低频拓展方法

Publications (2)

Publication Number Publication Date
CN108363739A CN108363739A (zh) 2018-08-03
CN108363739B true CN108363739B (zh) 2021-07-30

Family

ID=63006549

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810057585.9A Active CN108363739B (zh) 2018-01-22 2018-01-22 一种基于稀疏采集的地震资料高低频拓展方法

Country Status (1)

Country Link
CN (1) CN108363739B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111665542B (zh) * 2019-03-06 2023-04-25 中国石油天然气集团有限公司 地震数据拓频方法及系统
CN114609664B (zh) * 2020-12-08 2024-05-28 中国石油天然气集团有限公司 地震数据低频信号补偿方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105259580A (zh) * 2015-10-30 2016-01-20 中国石油大学(北京) 一种可控震源信号低频拓展方法
CN106405645A (zh) * 2016-08-30 2017-02-15 英得赛斯科技（北京）有限公司 一种基于资料品质分析的信噪比可控的地震拓频处理方法
CN106597532A (zh) * 2016-11-14 2017-04-26 中国石油化工股份有限公司 一种结合井资料与层位资料的叠前地震数据频带拓展方法
CN107219551A (zh) * 2016-03-22 2017-09-29 中国石油化工股份有限公司 拓宽地震数据频带的方法及装置
CN107436450A (zh) * 2017-07-26 2017-12-05 西安交通大学 一种基于连续小波变换的地震信号带宽拓展方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8737165B2 (en) * 2010-10-01 2014-05-27 Westerngeco L.L.C. Interferometric seismic data processing for a towed marine survey

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105259580A (zh) * 2015-10-30 2016-01-20 中国石油大学(北京) 一种可控震源信号低频拓展方法
CN107219551A (zh) * 2016-03-22 2017-09-29 中国石油化工股份有限公司 拓宽地震数据频带的方法及装置
CN106405645A (zh) * 2016-08-30 2017-02-15 英得赛斯科技（北京）有限公司 一种基于资料品质分析的信噪比可控的地震拓频处理方法
CN106597532A (zh) * 2016-11-14 2017-04-26 中国石油化工股份有限公司 一种结合井资料与层位资料的叠前地震数据频带拓展方法
CN107436450A (zh) * 2017-07-26 2017-12-05 西安交通大学 一种基于连续小波变换的地震信号带宽拓展方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于高性能稀疏采集的低频恢复技术在胜利东风港探区的应用;于富文等;《海洋地质前沿》;20170331;第33卷(第3期);全文 *

Also Published As

Publication number Publication date
CN108363739A (zh) 2018-08-03

Similar Documents

Publication Publication Date Title
CN111158049B (zh) 一种基于散射积分法的地震逆时偏移成像方法
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
CN108845317B (zh) 一种基于分层介质格林函数的频域逆时偏移方法
CN108845357B (zh) 一种基于同步挤压小波变换估计地层等效品质因子的方法
CN108363739B (zh) 一种基于稀疏采集的地震资料高低频拓展方法
CN105572723B (zh) 基于自相关子波的可控震源扫描信号的设计方法
CN109884709B (zh) 一种基于面波旅行时层析的转换波静校正方法
CN105242309A (zh) 一种地震数据中规则干扰的压制方法及装置
CN109633749B (zh) 基于散射积分法的非线性菲涅尔体地震走时层析成像方法
Zhang et al. Receiver function imaging with reconstructed wavefields from sparsely scattered stations
RU2412454C2 (ru) Способ обработки сейсмических данных с использованием дискретного вейвлет-преобразования
CN105510975B (zh) 提高地震数据信噪比的方法及装置
CN107605470A (zh) 一种纵横波径向速度变化成像方法
CN104793237A (zh) 一种获得宽频可控震源扫描信号的方法和装置
CN105676292A (zh) 一种基于二维曲波变换的三维地震数据去噪方法
CN109425892B (zh) 地震子波的估计方法及系统
Shen et al. Partial frequency band match filtering based on high-sensitivity data: method and applications
CN110895349B (zh) 基于中值滤波处理需求的独立扫描信号的设计方法
CN111899749A (zh) 一种监听混凝土泵送管道作业声音的降噪方法
CN113721297B (zh) 利用速度特征压制沙丘鸣震的方法
Parolai Determination of dispersive phase velocities by complex seismic trace analysis of surface waves (CASW)
CN106368689B (zh) 一种快速生成测井成像图的方法
CN112526611A (zh) 表层地震波品质因子的提取方法及装置
CN107402406B (zh) 一种压制地震数据中大钻噪声的方法
CN112444881A (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