CN113552632A - 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统 - Google Patents

基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统 Download PDF

Info

Publication number
CN113552632A
CN113552632A CN202010325675.9A CN202010325675A CN113552632A CN 113552632 A CN113552632 A CN 113552632A CN 202010325675 A CN202010325675 A CN 202010325675A CN 113552632 A CN113552632 A CN 113552632A
Authority
CN
China
Prior art keywords
wavelet
seismic
cross
correlation
wave field
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
Application number
CN202010325675.9A
Other languages
English (en)
Other versions
CN113552632B (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 University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
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 University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN202010325675.9A priority Critical patent/CN113552632B/zh
Publication of CN113552632A publication Critical patent/CN113552632A/zh
Application granted granted Critical
Publication of CN113552632B publication Critical patent/CN113552632B/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
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/148—Wavelet transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0248—Filters characterised by a particular frequency response or filtering method
    • H03H17/0255—Filters based on statistics
    • H03H17/0257—KALMAN filters

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明属于地震勘探技术领域，涉及一种基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统，包括以下步骤：S1将各检测点获得的地震信号进行小波转换，并获取小波系数；S2将其中两个检测点获得的地震信号的小波系数进行互相关，获得小波互相关函数；S3根据小波互相关函数，获得小波互相关的图像，并在图像中获取一地震道的目标波场；S4通过卡尔曼滤波，根据地震道的目标波场获取地震道集的目标波场；S5将地震道集的目标波场转换为新的小波互相关函数，并将新的小波系数函数转换至时间偏移距域，获得地震相关信号。其应用了多维度的卡尔曼滤波，能有效压制互相关数据中的噪声，并提高小波域互相关的质量。

Description

基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统
技术领域
本发明是关于一种基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统，属于地震勘探技术领域。
背景技术
地震互相关信号在地球物理学中起着至关重要的作用。例如，在利用旅行时求取速度结构时，可以从观测记录和正演模拟记录的互相关中提取相移信息，以约束地下速度模型。当地震互相关信号波形简单且相似时，可以通过互相关峰值相的滞后时间，来估算两条地震道之间的时差。地震互相关信号另一个重要的应用是地震干涉法：同一震源激发，将两个检波器记录的地震信号互相关并叠加，可以计算出一个检波器激发、另外一个检波器接收的格林函数。该方法既可以应用于勘探中的虚拟源成像，也适用于大地球物理中的背景噪声成像。背景噪声互相关重构的格林函数包括台站之间传播的面波等的多种震相，可用于估算地下不同尺度的介质速度。但是，互相关方法通常假定输入的数据是固定的，然而由于多种波场的存在，实际地震数据并非如此。而且在地震干涉法中，大多数采集的数据无法满足震源封闭的假设条件，因此得到的互相关格林函数往往低信噪比较低。所以，需要对干涉数据进行处理恢复其覆盖范围，使互相关处理有效。
发明内容
针对上述现有技术的不足，本发明的目的是提供了一种基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统，其应用了多维度的卡尔曼滤波，能有效压制互相关数据中的噪声，并提高小波域互相关的质量。
为实现上述目的，本发明提供了一种基于小波域卡尔曼滤波的地震互相关信号拾取方法，包括以下步骤：S1将各检测点获得的地震信号进行小波转换，并获取小波系数；S2将任意两个检测点获得的地震信号的小波系数进行互相关，获得小波互相关函数；S3根据小波互相关函数，获得小波互相关的图像，并在图像中获取一地震道的目标波场；S4通过卡尔曼滤波，根据地震道的目标波场获取地震道集的目标波场；S5将地震道集的目标波场转换为新的小波互相关函数，并将新的小波系数函数转换至时间偏移距域，获得地震相关信号。
进一步，S1中小波系数的公式为：
Figure BDA0002463126710000011
其中，
Figure BDA0002463126710000021
是小波系数，X(t)是t时刻获得的地震信号，
Figure BDA0002463126710000022
是Ricker母小波，符号f和λ分别是小波频率和平移量，T为测试时长。
进一步，S2中的小波域互相关函数WM(f,τ)为：
Figure BDA0002463126710000023
其中，
Figure BDA0002463126710000024
和
Figure BDA0002463126710000025
是分别在检波器A和检波器B处获取的噪声波场的小波系数函数，τ为检波器A和检波器B之间的时间延迟，f和λ分别是小波频率和平移量。
进一步，步骤S3中，在图像中获取一地震道的目标波场的方法为：对小波互相关的图像窗口化，根据每个窗口附近小波系数的局部平均幅度来选择阈值，并保持小波系数的幅度高于阈值的小波系数，将小波系数的幅度低于阈值的小波系数的值设置为零。
进一步，地震道的目标波场对应的互相关函数为：
WM0(f,τ)＝W(f,τ)*f,τWM(f,τ)
Figure BDA0002463126710000026
其中，*为褶积算子，W(f,τ)为窗口长度调整函数，WM(f,τ)为原小波互相关函数，WM0(f,τ)为小波时频域中窗口化之后输出的小波互相关函数，WMi(f,τ)为目标波场对应的小波互相关函数，η为阈值，f和λ分别是小波频率和平移量。
进一步，通过卡尔曼滤波获取地震道集的目标波场的方法为：基于已经获得的地震道的目标波场，通过预测函数预测下一条地震道的目标波场位置和大小，得到下一条地震道的预测目标波场；通过步骤S1-S3，获取下一条地震道的实际目标波场，计算预测目标波场和实际目标波场的协方差矩阵；计算卡尔曼增益矩阵，并根据卡尔曼增益矩阵优化预测函数。
进一步，初始的预测函数的表达式为：
Figure BDA0002463126710000027
初始的协方差矩阵为：
Figure BDA0002463126710000028
其中，k是偏移序列索引，Yk是偏移距k的预测向量，Pk是与Yk有关的预测误差协方差，N和Q分别是预测噪声和预测协方差，变量Ф是模型矩阵，ФT是模型矩阵的转置矩阵。
进一步，卡尔曼增益矩阵的表达式为：
Figure BDA0002463126710000029
其中，G是卡尔曼增益，R是测量协方差，H是模型矩阵，是常量，HT是模型矩阵的转置矩阵。
进一步，更新后的预测函数的表达式为：
Figure BDA0002463126710000031
更新后的协方差矩阵：
Figure BDA0002463126710000032
其中，Zk+1是偏移距k+1的测量向量，I是单位矩阵。
本发明还公开了一种基于小波域卡尔曼滤波的地震互相关信号拾取系统，包括：小波系数获取模块，用于将各检测点获得的地震信号进行小波转换，并获取小波系数；小波互相关函数获取模块，用于将其中两个检测点获得的地震信号的小波系数进行互相关，获得小波互相关函数；目标波场获取模块，用于根据小波互相关函数，获得小波互相关的图像，并在图像中获取一地震道的目标波场；卡尔曼滤波模块，用于通过卡尔曼滤波，根据地震道的目标波场获取地震道集的目标波场；输出模块，用于将地震道集的目标波场转换为新的小波互相关函数，并将小波系数函数转换至时间偏移距域，获得地震相关信号。
本发明由于采取以上技术方案，其具有以下优点：
1、本方法集成了小波变换、互相关、图像目标检测和卡尔曼滤波目标追踪的方法，对目标波场进行有效的测量。
2、本发明能够有效地衰减不需要的波场，并产生高信噪比的互相关函数。
3、本发明提出了一种自动化的地震相关信号处理方法，无需大量人为干预，保证了处理结果的一致性。
附图说明
图1是本发明一实施例中基于小波域卡尔曼滤波的地震互相关信号拾取方法的流程图；
图2是本发明一实施例中基于小波域卡尔曼滤波的地震互相关信号转换结果过程图，图2(a)为经过小波转换的地震信号对应的小波系数图像；图2(b)为局部自适应阈值矩阵图；图2(c)是根据小波系数的幅值将信号分为若干区域的示意图；图2(d)是地震道的目标波场的示意图；图2(e)是经过卡尔曼滤波后的小波系数函数图像；图2(f)是时域中的原始地震信号和经过卡尔曼滤波后的地震信号的对比图，其中浅色的波形为原始地震信号，深色的波形是经过卡尔曼滤波后的地震信号；
图3是2004BP模型的地震信号深度-偏移距图，图3(a)是地震信号的真实速度模型，图3(b)是经过初始平滑获得的地震信号P波速度模型；
图4是2004BP模型的地震信号时间-偏移距图，图4(a)是地震信号的真实速度的地震信号，图4(b)是经过初始平滑获得的地震信号P波速度的地震信号；
图5是2004BP模型的近偏移距地震道的示意图，图5(a)是观测地震信号和模拟地震信号之间的互相关道集的示意图；图5(b)是经过小波转换的地震信号的小波系数；图5(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵，白色虚线框为卡尔曼滤波器的相应预测值；图5(d)是经过卡尔曼滤波的的小波系数函数图像；
图6是2004BP模型的中等偏移距地震道的示意图，图6(a)是观测地震信号和模拟地震信号之间的互相关道集的示意图；图6(b)是经过小波转换的地震信号的小波系数；图6(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵，白色虚线框为卡尔曼滤波器的相应预测值；图6(d)是经过卡尔曼滤波的的小波系数函数图像；
图7是2004BP模型的大等偏移距地震道的示意图，图7(a)是观测地震信号和模拟地震信号之间的互相关道集的示意图；图7(b)是经过小波转换的地震信号的小波系数；图7(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵，白色虚线框为卡尔曼滤波器的相应预测值；图7(d)是经过卡尔曼滤波的的小波系数函数图像；
图8(a)是2004BP模型目标波场的中心位置在测量值和预测值之间的残差在频率方向的分量，图8(b)是2004BP模型目标波场的中心位置在测量值和预测值之间的残差在时间方向的分量。
图9是2004BP模型中的互相关道集信号，其中，图9(a)是常规互相关的互相关道集信号，图9(b)是本发明中方法获得的互相关信号，图中黑线标记为最大峰值给出的时间延迟；
图10(a)是井下的井轨迹(图中线性部分)和检波器(图中的三角形部分)的位置示意图，图10(b)是随机选择由10个检波器串记录的1.2秒长的微地震背景噪声道集。
图11是实施例三中的零偏移距地震道的示意图，图11(a)是根据井下阵列记录的背景噪声重构的互相关道集的示意图；图11(b)是经过小波转换的地震信号的小波系数；图11(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图11(d)是经过卡尔曼滤波的的小波系数函数图像；
图12是实施例三中的大偏移距地震道的示意图，图12(a)是根据井下阵列记录的背景噪声重构的互相关道集的示意图；图12(b)是经过小波转换的地震信号的小波系数；图12(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图12(d)是经过卡尔曼滤波的的小波系数函数图像；
图13是实施例三中的另一条大偏移距地震道的示意图，图13(a)是根据井下阵列记录的背景噪声重构的互相关道集的示意图；图13(b)是经过小波转换的地震信号的小波系数；图13(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图13(d)是经过卡尔曼滤波的的小波系数函数图像；
图14(a)是实施例三中目标波场的中心位置在测量值和预测值之间的残差在频率方向的分量，图14(b)是实施例三中目标波场的中心位置在测量值和预测值之间的残差在时间方向的分量。
图15是实施例三中的互相关道集信号，其中，图15(a)是常规互相关的互相关道集信号，图15(b)是本发明中方法获得的互相关信号。
图16是分布在加勒比海周围的99个宽带地震台站(倒三角)的位置以及通往G.FDF台的射线路径示意图。
图17是实施例四中的互相关道集信号，其中，图17(a)是常规互相关的互相关道集信号，图17(b)是本发明中方法获得的互相关信号。
具体实施方式
为了使本领域技术人员更好的理解本发明的技术方向，通过具体实施例对本发明进行详细的描绘。然而应当理解，具体实施方式的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。在本发明的描述中，需要理解的是，所用到的术语仅仅是用于描述的目的，而不能理解为指示或暗示相对重要性。
互相关是用来度量两个地震道相似性的算法，它是一个地震相对于另一个地震相的时间延迟的函数。在数学上，两条地震地震道之间的互相关函数可以写成：
Figure BDA0002463126710000051
其中XA和XB分别代表观测数据和计算数据。寻找一个τ(小于记录时间T的绝对值)对计算的地震道进行平移，使得互相关函数M(τ)最大化。Δτ＝0表明已知且正确的速度模型，该模型会生成与观测数据同时到达的波。
在地震干涉测量法中，XA和XB分别表示在两个接收器的空间坐标处记录到的被动地震数据集。将A视为虚拟源时在检波器B处得到M(τ)，具有相关时滞长度τ的干涉数据。
在实际应用中，XA和通XB常包含多波模式，会导致数据的不稳定：
Figure BDA0002463126710000052
Figure BDA0002463126710000053
其中下标tar和n分别是在检波器A和检波器B处记录的目标波场和噪声波场。只有
Figure BDA0002463126710000054
与
Figure BDA0002463126710000055
之间的互相关会形成正确的运动事件，而其他项通常会产生伪事件。通常，目标波场是在旅行时层析成像中的传输波，即地震干涉法中在井眼周围传播的P波，或背景噪声干涉法中的瑞利波。
实施例一
本实施例公开了一种基于小波域卡尔曼滤波的地震互相关信号拾取方法，如图1、图2所示，包括以下步骤：
S1将各检测点获得的地震信号进行小波转换，并获取小波系数。
如图2(a)所示，经过小波转换的地震信号对应的小波系数的公式为：
Figure BDA0002463126710000061
其中，
Figure BDA0002463126710000062
是小波系数，X(t)是t时刻获得的地震信号，
Figure BDA0002463126710000063
是Ricker母小波，符号f和λ分别是小波频率和平移量，T为测试时长。该小波系数的函数在小波域中生成的图像为图2(a)。
S2将任意两个检测点获得的地震信号的小波系数进行互相关，获得小波互相关函数，设此处两个检测点为检测点A和检测点B。
小波域互相关函数WM(f,τ)为：
Figure BDA0002463126710000064
其中，
Figure BDA0002463126710000065
和
Figure BDA0002463126710000066
是分别在检波器A和检波器B处获取的噪声波场的小波系数函数，τ为检波器A和检波器B之间的时间延迟，f和λ分别是小波频率和平移量。
S3根据小波互相关函数，获得小波互相关的图像，并在图像中获取一地震道的目标波场。
在图像中获取一地震道的目标波场的方法为局部自适应阈值算法。局部自适应阈值算法的具体过程为：对小波互相关的图像窗口化，根据每个窗口附近小波系数的局部平均幅度来选择阈值，所有的阈值组成的阈值矩阵如图2(b)所示，根据小波系数的幅值将小波系数分为若干区域，并采用不同的颜色深浅将不同区域进行标注，如图2(c)所示，保持小波系数的幅度高于阈值的小波系数，将小波系数的幅度低于阈值的小波系数的值设置为零，最终选出地震道的目标波场如图2(d)所示。
上述局部自适应阈值算法的过程可以采用如下公式表示，即将地震道的目标波场对应的互相关函数表示为：
WM0(f,τ)＝W(f,τ)*f,τWM(f,τ)
Figure BDA0002463126710000067
其中，*为褶积算子，W(f,τ)为窗口长度调整函数，WM(f,τ)为原小波互相关函数，WM0(f,τ)为小波时频域中窗口化之后输出的小波互相关函数，WMi(f,τ)为目标波场对应的小波互相关函数，η为阈值，f和λ分别是小波频率和平移量。
S4通过卡尔曼滤波，根据地震道的目标波场获取地震道集的目标波场。
如图2(e)、(f)所示，通过卡尔曼滤波获取地震道集的目标波场的方法为：基于已经获得的地震道的目标波场，通过预测函数预测下一条地震道的目标波场位置和大小，得到下一条地震道的预测目标波场；通过步骤S1-S3，获取下一条地震道的实际目标波场，计算预测目标波场和实际目标波场的协方差矩阵；计算卡尔曼增益矩阵，并根据卡尔曼增益矩阵优化预测函数。
初始的预测函数的表达式为：
Figure BDA0002463126710000071
初始的协方差矩阵为：
Figure BDA0002463126710000072
其中，k是偏移序列索引，Yk是偏移距k的预测向量，Pk是与Yk有关的预测误差协方差，N和Q分别是预测噪声和预测协方差，变量Ф是模型矩阵，假定为常数，ФT是模型矩阵的转置矩阵。
卡尔曼增益矩阵的表达式为：
Figure BDA0002463126710000073
其中，G是卡尔曼增益，R是测量协方差，H是模型矩阵，是常量，HT是模型矩阵的转置矩阵。
更新后的预测函数的表达式为：
Figure BDA0002463126710000074
更新后的协方差矩阵：
Figure BDA0002463126710000075
其中，Zk+1是偏移距k+1的测量向量，I是单位矩阵。
S5将地震道集的目标波场转换为新的小波互相关函数，并将新的小波系数函数转换至时间偏移距域，获得地震相关信号。
卡尔曼滤波是目前应用最为广泛的滤波方法之一，在通信、导航、图像视频分析等多领域都得到了较好的应用。它具有便于编程实现、实时动态处理数据、计算量小的优点。卡尔曼滤波是一种线性系统状态算法，通过系统输入输出观测数据，对系统状态进行最优估计，在测量方差已知的情况下能够从包含噪声的数据中，估计动态系统的状态。卡尔曼滤波还可以作为强大的空间算法，通过距离测量来进行预测，对一系列距离测量进行处理，来递归更新系统的预测值。每次空间测量后，滤波器都会在测量阶段产生新的偏移距预测。由于卡尔曼滤波使用的测量值中包含随机噪声和其他不确定因素，因此该算法能够得到未知变量的估计值，这比仅基于单一值所得到的预测结果更准确。本方法集成了小波变换、互相关、图像目标检测和卡尔曼滤波目标追踪的方法，对目标波场进行有效的测量，能够有效地衰减不需要的波场，并产生高信噪比的互相关函数，无需大量人为干预，保证了处理结果的一致性。
实施例二
为了进一步说明实施例一中的技术方案，本实施例以勘探地震学的2004BP模型。如图3所示，图3是2004BP模型的地震信号深度-偏移距图，图3(a)是地震信号的真实速度模型，图3(b)是经过初始平滑获得的地震信号P波速度模型。图3(b)是图3(a)中模型进过初始平滑后得到的模型，将图3(b)中模型作为初始地震信号，带入实施例一中技术方案。图4是2004BP模型的地震信号时间-偏移距图，图4(a)是地震信号的真实速度的地震信号，图4(b)是经过初始平滑获得的地震信号P波速度的地震信号。对模型中的数据集进行顶切处理，该模型对应的数据集包含从50km到87km的地震道，将经过地震道估计的零偏移数据用作卡尔曼滤波器追踪系统的初始值。
图5、图6和图7中分别是2004BP模型的近偏移距地震道、中等偏移距地震道和大偏移距地震道的示意图，其中，图(a)是观测地震信号和模拟地震信号之间的互相关道集的示意图；图(b)是经过小波转换的地震信号的小波系数；图(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图(d)是经过卡尔曼滤波的的小波系数函数图像。
通过观测道集(即图3(a)和图4(a)中的地震信号)和模拟道集(即图3(b)和图4(b)中的地震信号)之间的常规互相关来创建互相关道集。选定的近偏移地震道用垂直的虚线在图5(a)中标记，同理，选定的中等偏移地震道用垂直的虚线在图6(a)中标记，选定的大偏移地震道用垂直的虚线在图7(a)中标记。从图5(c)、图6(c)和图7(c)中可以看出通过实施例一中方法计算得到的目标波场与卡尔曼滤波预测的目标波场边界吻合，验证了实施例一中方案的有效性，而且不论是近偏距、中等偏距还是高偏距采用实施例一中方法计算的目标波场均匀测量的范围吻合。图8(a)是目标波场的中心位置在测量值和预测值之间的残差在频率方向的分量，图8(b)是目标波场的中心位置在测量值和预测值之间的残差在时间方向的分量。通过在图8中绘制追踪误差进一步检验测量值和预测值的相关性。图8(a)显示了y轴方向上追踪位置的预测值和检测值并不完全匹配，图8(b)显示了关于时间方向上追踪位置的预测值和检测值也不匹配。观察结果支持了追踪误差的随机性和独立性这一观点，其中总误差始终低于2％，说明测量值和预测值之间具有较好的相关性。图9是2004BP模型中的互相关道集信号，其中，图9(a)是常规互相关的互相关道集信号，图9(b)是本发明中方法获得的互相关信号，图中黑线标记为最大峰值给出的时间延迟。从图9(a)中可以看出至少存在三个较大的时间延迟，而图9(b)中基本不存在明显的实现延迟，这表明通过实施例一中的技术方案获得的互相关信号时间延迟小，相关性高。这些观察结果证实了实施例一中的技术方案具有较好的稳健性。此外，通过实施例一中的技术方案获得的互相关信号具有更好的信噪比，可以较好的应用于其他研究，例如全波形反演。
实施例三
本实施例以四川省的垂直井中的水力压裂监测系统为例对实施例一中的技术方案进行说明。上述水力压裂监测系统的具体结构如图10所示，图10(a)是井下的井轨迹(图中线性部分)和检波器(图中的三角形部分)的位置示意图，其中检波器覆盖的侧面倾角小于2°，图10(b)是随机选择由10个检波器串记录的1.2秒长的微地震背景噪声道集。检波器沿井轨迹分布，微地震背景噪声包含从未知震源到达检波器的波。
图11是实施例三中的零偏移距地震道的示意图，图11(a)是根据井下阵列记录的背景噪声重构的互相关道集的示意图；图11(b)是经过小波转换的地震信号的小波系数；图11(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图11(d)是经过卡尔曼滤波的的小波系数函数图像。图11(a)中信号包括10天的地震信号原始叠加记录，由于震源分布不均匀且在时间上不具有周期性，信噪比较低。将最深的检波器视为源，对其他检波器进行互相关来获得上行P波。产生P波的噪声源位于检波器上方，可能与套管上的活动有关。
图11、图12和图13中分别是本实施中的零偏移距地震道、大偏移距地震道和另一条大偏移距地震道的示意图，其中，图(a)是根据井下阵列记录的背景噪声重构的互相关道集的示意图；图(b)是经过小波转换的地震信号的小波系数；图(c)是经过局部自适应阈值算法得出的两个波场的标记矩阵；图(d)是经过卡尔曼滤波的的小波系数函数图像。将地震道估计后的自相关(零偏移)检测用作卡尔曼滤波器追踪系统的初始值。自相关地震道用垂直的虚线进行标记。从图11(c)、图11(c)和图11(c)中可以看出通过实施例一中方法计算得到的目标波场与卡尔曼滤波预测的目标波场边界吻合，验证了实施例一中方案的有效性。图14(a)是实施例三中目标波场的中心位置在测量值和预测值之间的残差在频率方向的分量，图14(b)是实施例三中目标波场的中心位置在测量值和预测值之间的残差在时间方向的分量。通过在图14中绘制追踪误差进一步检验量值和预测值的相关性。图14(a)显示了y轴方向上追踪位置的预测值和检测值并不完全匹配，图14(b)显示了关于时间方向上追踪位置的预测值和检测值也不匹配。观察结果支持了追踪误差的随机性和独立性这一观点，其中总误差始终低于1％，说明测量值和预测值之间具有较好的相关性。图15是实施例三中的互相关道集信号，其中，图15(a)是常规互相关的互相关道集信号，图15(b)是本发明中方法获得的互相关信号。信号沿图10(a)中底部到顶部的方向传播，每个地震道的最大值由黑线连接，斜率不一致的。预测值在结果峰值上显示出良好的信号连续性，且降低了其他波场的污染。
实施例四
本实施例使用了加勒比海周围部署的99个宽带站垂直部分的数据，其分别方式如图16所示。通过以下方式处理获得的地震连续信号：去仪器响应、去均值、去线性趋势和波形尖灭，将记录降采至1Hz，0.02-0.2Hz带通滤波，按天分割成(84600s)片段。接着应用时域归一化来消除地震信号的干扰，并通过谱白化来增强背景噪声。最后计算出每天的互相关，并根据相位加权叠加。由于信号来源不均匀、叠加时常不足，得到的互相关结果仍会出现低信噪比问题，如图17a所示。而使用本发明的方法得到的目标波场具有高度相关性，如图17b所示。相比之下，图17(b)中黑线斜率明显大于图17(a)，这说明实施例一中方法提高了道集的信噪比。其中图17(b)中黑线斜率表示平均视速度为2.73km/s。
以上三个实施例的结果表明，实施例一中的技术方案可以有效地分离目标波场，并且无需任何其他信息即可产生高信噪比相关道集。在此框架内，通过减少其他波场对目标信号的影响，为勘探和全球地震学的各种应用提供了有价值的互相关输入，证明了本发明的技术价值。
实施例五
基于相同的发明构思，本实施例还公开了一种基于小波域卡尔曼滤波的地震互相关信号拾取系统，包括：
小波系数获取模块，用于将各检测点获得的地震信号进行小波转换，并获取小波系数；
小波互相关函数获取模块，用于将任意两个检测点获得的地震信号的小波系数进行互相关，获得小波互相关函数；
目标波场获取模块，用于根据小波互相关函数，获得小波互相关的图像，并在图像中获取一地震道的目标波场；
卡尔曼滤波模块，用于通过卡尔曼滤波，根据地震道的目标波场获取地震道集的目标波场；
输出模块，用于将地震道集的目标波场转换为新的小波互相关函数，并将小波系数函数转换至时间偏移距域，获得地震相关信号。
上述内容仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围为准。

Claims (10)

1.一种基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，包括以下步骤：
S1将各检测点获得的地震信号进行小波转换，并获取小波系数；
S2将任意两个所述检测点获得的地震信号的所述小波系数进行互相关，获得小波互相关函数；
S3根据所述小波互相关函数，获得小波互相关的图像，并在所述图像中获取一地震道的目标波场；
S4通过卡尔曼滤波，根据所述地震道的目标波场获取地震道集的目标波场；
S5将所述地震道集的目标波场转换为新的小波互相关函数，并将所述新的小波系数函数转换至时间偏移距域，获得地震相关信号。
2.如权利要求1所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，所述S1中小波系数的公式为：
Figure FDA0002463126700000011
其中，
Figure FDA0002463126700000012
是小波系数，X(t)是t时刻获得的地震信号，
Figure FDA0002463126700000013
是Ricker母小波，符号f和λ分别是小波频率和平移量，T为测试时长。
3.如权利要求1所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，所述S2中的小波域互相关函数WM(f,τ)为：
Figure FDA0002463126700000014
其中，
Figure FDA0002463126700000015
和
Figure FDA0002463126700000016
是分别在检波器A和检波器B处获取的噪声波场的小波系数函数，τ为检波器A和检波器B之间的时间延迟，f和λ分别是小波频率和平移量。
4.如权利要求3所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，所述步骤S3中，在所述图像中获取一地震道的目标波场的方法为：对所述小波互相关的图像窗口化，根据每个窗口小波系数的局部平均幅度来选择阈值，并保持小波系数的幅度高于所述阈值的小波系数，将小波系数的幅度低于所述阈值的小波系数的值设置为零。
5.如权利要求4所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，所述地震道的目标波场对应的互相关函数WMT(f,τ)为：
WM0(f,τ)＝W(f,τ)*f,τWM(f,τ)
Figure FDA0002463126700000021
其中，*为褶积算子，W(f,τ)为窗口长度调整函数，WM(f,τ)为原小波互相关函数，WM0(f,τ)为小波时频域中窗口化之后输出的小波互相关函数，WMi(f,τ)为目标波场对应的小波互相关函数，η为阈值，f和λ分别是小波频率和平移量。
6.如权利要求1-5任一项所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，通过卡尔曼滤波获取地震道集的目标波场的方法为：
基于已经获得的地震道的目标波场，通过预测函数预测下一条地震道的目标波场位置和大小，得到下一条地震道的预测目标波场；
通过所述步骤S1-S3，获取下一条地震道的实际目标波场，计算预测目标波场和实际目标波场的协方差矩阵；
计算卡尔曼增益矩阵，并根据所述卡尔曼增益矩阵优化所述预测函数。
7.如权利要求6所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，初始的所述预测函数的表达式为：
Figure FDA0002463126700000022
初始的协方差矩阵为：
Figure FDA0002463126700000023
其中，k是偏移序列索引，Yk是偏移距k的预测向量，Pk是与Yk有关的预测误差协方差，N和Q分别是预测噪声和预测协方差，变量Ф是模型矩阵，ФT是模型矩阵的转置矩阵。
8.如权利要求7所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，所述卡尔曼增益矩阵的表达式为：
Figure FDA0002463126700000024
其中，G是卡尔曼增益，R是测量协方差，H是模型矩阵，HT是模型矩阵的转置矩阵。
9.如权利要求7所述的基于小波域卡尔曼滤波的地震互相关信号拾取方法，其特征在于，更新后的所述预测函数的表达式为：
Figure FDA0002463126700000025
更新后的协方差矩阵：
Figure FDA0002463126700000026
其中，Zk+1是偏移距k+1的测量向量，I是单位矩阵。
10.一种基于小波域卡尔曼滤波的地震互相关信号拾取系统，其特征在于，包括：
小波系数获取模块，用于将各检测点获得的地震信号进行小波转换，并获取小波系数；
小波互相关函数获取模块，用于将任意两个所述检测点获得的地震信号的所述小波系数进行互相关，获得小波互相关函数；
目标波场获取模块，用于根据所述小波互相关函数，获得小波互相关的图像，并在所述图像中获取一地震道的目标波场；
卡尔曼滤波模块，用于通过卡尔曼滤波，根据所述地震道的目标波场获取地震道集的目标波场；
输出模块，用于将所述地震道集的目标波场转换为新的小波互相关函数，并将所述小波系数函数转换至时间偏移距域，获得地震相关信号。
CN202010325675.9A 2020-04-23 2020-04-23 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统 Active CN113552632B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010325675.9A CN113552632B (zh) 2020-04-23 2020-04-23 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010325675.9A CN113552632B (zh) 2020-04-23 2020-04-23 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统

Publications (2)

Publication Number Publication Date
CN113552632A true CN113552632A (zh) 2021-10-26
CN113552632B CN113552632B (zh) 2022-11-01

Family

ID=78129325

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010325675.9A Active CN113552632B (zh) 2020-04-23 2020-04-23 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统

Country Status (1)

Country Link
CN (1) CN113552632B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116881385A (zh) * 2023-09-08 2023-10-13 中国铁塔股份有限公司 轨迹平滑方法、装置、电子设备及可读存储介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103873743A (zh) * 2014-03-24 2014-06-18 中国人民解放军国防科学技术大学 一种基于结构张量和卡尔曼滤波的视频去噪方法
CN105093275A (zh) * 2014-05-09 2015-11-25 中国石油化工股份有限公司 基于速度模型的有效信号提取iss波场分离方法
CN105652325A (zh) * 2016-01-05 2016-06-08 吉林大学 基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法
US9703002B1 (en) * 2003-10-04 2017-07-11 SeeScan, Inc. Utility locator systems and methods
CN109782347A (zh) * 2019-01-18 2019-05-21 南京邮电大学 一种基于小波分析的地震走时反演成像方法
US20190158777A1 (en) * 2015-12-30 2019-05-23 Steve Mann Recompressive sensing, resparsified sampling, and lightspacetimelapse: means, apparatus, and methods for spatiotemporal and spatiotonal timelapse and infinitely long media or multimedia recordings in finite memory
EP3542189A1 (en) * 2016-11-17 2019-09-25 Saudi Arabian Oil Company Use of wavelet cross-correlation for virtual source denoising

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9703002B1 (en) * 2003-10-04 2017-07-11 SeeScan, Inc. Utility locator systems and methods
CN103873743A (zh) * 2014-03-24 2014-06-18 中国人民解放军国防科学技术大学 一种基于结构张量和卡尔曼滤波的视频去噪方法
CN105093275A (zh) * 2014-05-09 2015-11-25 中国石油化工股份有限公司 基于速度模型的有效信号提取iss波场分离方法
US20190158777A1 (en) * 2015-12-30 2019-05-23 Steve Mann Recompressive sensing, resparsified sampling, and lightspacetimelapse: means, apparatus, and methods for spatiotemporal and spatiotonal timelapse and infinitely long media or multimedia recordings in finite memory
CN105652325A (zh) * 2016-01-05 2016-06-08 吉林大学 基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法
EP3542189A1 (en) * 2016-11-17 2019-09-25 Saudi Arabian Oil Company Use of wavelet cross-correlation for virtual source denoising
CN109782347A (zh) * 2019-01-18 2019-05-21 南京邮电大学 一种基于小波分析的地震走时反演成像方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
SZU-WEI CHOU ET AL.: "Wavelet-Based Method for Time-Domain Noise Analysis and Reduction in a Frequency-Scan Ion Trap Mass Spectrometer", 《AMERICAN SOCIETY FOR MASS SPECTROMETRY, 2012》 *
YANG ZHAO ET AL.: "WAVELET-CROSSCORRELATION-BASED INTERFEROMETRIC REDATUMING IN 4D SEISMIC", 《GEOPHYSICS》 *
王文等: "多尺度卡尔曼滤波在图像去噪中的应用研究", 《弹箭与制导学报》 *
谭平等: "基于无抽取Haar算法的实时卡尔曼滤波方法", 《中南大学学报（自然科学版）》 *
赵娟等: "基于小波变换的分形随机信号的卡尔曼滤波", 《电子学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116881385A (zh) * 2023-09-08 2023-10-13 中国铁塔股份有限公司 轨迹平滑方法、装置、电子设备及可读存储介质
CN116881385B (zh) * 2023-09-08 2023-12-01 中国铁塔股份有限公司 轨迹平滑方法、装置、电子设备及可读存储介质

Also Published As

Publication number Publication date
CN113552632B (zh) 2022-11-01

Similar Documents

Publication Publication Date Title
RU2518577C2 (ru) Непрерывный адаптивный анализ поверхностных волн в случае трехмерных сейсмических данных
EA012636B1 (ru) Способ предсказания кратных волн, связанных с поверхностью, на основе данных буксируемой морской сейсмической косы с двумя типами датчиков
GB2444953A (en) Regularisation of irregularly sampled seismic data using sinc functions
EA014282B1 (ru) Геофизический способ разведки
CN111025386B (zh) 一种无分离假象的纵横波分离方法
WO2008053135A2 (en) Method of determining properties of the earth by evaluating the wavefield of turning waves
CN111257939B (zh) 一种时移地震虚拟震源双向波场重构方法和系统
EP2073041A1 (en) Method to estimate a seismic ray parameter for a seismogram
US10379244B2 (en) Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume over a time period
CN107884829A (zh) 一种联合压制浅海obc地震资料多次波的方法
RU2412454C2 (ru) Способ обработки сейсмических данных с использованием дискретного вейвлет-преобразования
Zhao et al. Signal detection and enhancement for seismic crosscorrelation using the wavelet-domain Kalman filter
CN104570116A (zh) 基于地质标志层的时差分析校正方法
CA3224766A1 (en) P/s wave measurement and compensation
Dai et al. Effects due to aliasing on surface-wave extraction and suppression in frequency-velocity domain
CN114415234B (zh) 基于主动源面波频散和h/v确定浅地表横波速度的方法
AU2014101578A4 (en) System and method for seismic adaptive optics
CN113552632B (zh) 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统
CN110261899B (zh) 地震数据z字形干扰波去除方法
CN111257938A (zh) 基于小波互相关时移地震虚拟震源波场重构方法和系统
CN107632321A (zh) 一种偏移成像方法
US11313987B2 (en) Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume over a time period using time-lapse seismic
CN115113277A (zh) 时间域利用海洋地震直达波求远场震源子波的方法
GB2610201A (en) Method of simulating seismic data
Ehsaninezhad et al. Urban subsurface exploration improved by denoising of virtual shot gathers from distributed acoustic sensing ambient noise

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
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Zhao Yang

Inventor after: Zuo Jiahui

Inventor after: Zhou Yangtianli

Inventor after: Chen Haichao

Inventor before: Zhou Yangtianli

Inventor before: Zhao Yang

Inventor before: Zuo Jiahui

Inventor before: Chen Haichao

GR01 Patent grant
GR01 Patent grant