CN114859404A - 超采样地震波形匹配方法及装置 - Google Patents
超采样地震波形匹配方法及装置 Download PDFInfo
- Publication number
- CN114859404A CN114859404A CN202110074803.1A CN202110074803A CN114859404A CN 114859404 A CN114859404 A CN 114859404A CN 202110074803 A CN202110074803 A CN 202110074803A CN 114859404 A CN114859404 A CN 114859404A
- Authority
- CN
- China
- Prior art keywords
- time shift
- waveform
- shift amount
- point
- seismic waveform
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 93
- 238000005070 sampling Methods 0.000 title claims abstract description 69
- 230000009466 transformation Effects 0.000 claims abstract description 77
- 238000010219 correlation analysis Methods 0.000 claims abstract description 42
- 230000006870 function Effects 0.000 claims description 25
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 238000004590 computer program Methods 0.000 claims description 15
- 238000000605 extraction Methods 0.000 claims description 9
- 238000003860 storage Methods 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 5
- 238000010586 diagram Methods 0.000 description 22
- 230000008569 process Effects 0.000 description 11
- 238000009826 distribution Methods 0.000 description 10
- 238000012545 processing Methods 0.000 description 10
- 238000004458 analytical method Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 7
- 238000001228 spectrum Methods 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 5
- 238000006243 chemical reaction Methods 0.000 description 5
- 238000005755 formation reaction Methods 0.000 description 5
- 238000000354 decomposition reaction Methods 0.000 description 4
- 238000012952 Resampling Methods 0.000 description 3
- 238000011065 in-situ storage Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 208000010392 Bone Fractures Diseases 0.000 description 1
- 206010017076 Fracture Diseases 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
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/30—Analysis
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
技术领域
本发明涉及地球物理勘探技术领域，尤其涉及一种超采样地震波形匹配方法及装置。
背景技术
本部分旨在为权利要求书中陈述的本发明实施例提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。
地震勘探是一种利用地下介质弹性和密度差异，通过观测和分析地层对人工激发的地震波的响应来推断地下岩层的性质和形态的地球物理勘探方法。地震勘探是钻探前勘测石油天然气资源和固体矿产资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面也得到广泛应用。
利用地震勘探可以对地下地质构造进行成像。通过地震采集、叠加、偏移等复杂的处理过程，可形成地震数据体。地震数据体通常由一系列规则排列的地震道组成，每一个地震道可视为一个延续数秒的地震波形。地震波形是地下介质的地震响应，随地震波传播时间变化。在沿某一方向所抽取的地震剖面上，由相邻地震波形的同相性所形成的地震同相轴可指示地层界面的位置和形态。
利用地震波形及地震同相轴可对地下地质构造进行分析，获取构造形态信息，如地层的界面空间形态、地层厚度、地层中的断裂发育特征等。在这一过程中，地震同相轴的倾角信息(因地震同相轴倾角主要反映地层的倾斜形态，通常将其称之为地层倾角)扮演着关键的角色。要利用地震波形特征求取地层倾角，对相邻地震波形进行时间域匹配是一个有效途径。
多年来，针对波形类信号的匹配方法得到了深入研究，出现了众多的波形匹配技术，但归纳起来主要分为三大类：特征提取法、预测法和距离法。特征提取法的主要代表是奇异值分解，预测法的主要代表是平面波分解法，距离法的主要代表是欧几里得距离法和相关法。在这些代表性算法中，奇异值分解法相当复杂，且波形匹配的精度依赖于特征提取的类型及数量；平面波分解法运算量很大，欧几里得距离法虽然运算简单但当信噪比较低时匹配精度无法满足要求。与以上其它算法相比，相关法波形匹配技术，具有良好的抗干扰性和稳定性，在波形匹配中得到广泛应用且常被视作是波形匹配精度的衡量标准。
相关法波形匹配的技术核心是相关分析。相关分析(Correlation analysis)是指研究两个随机变量或现象之间是否存在某种依存关系，并对有依存关系的变量或现象探讨它们的相关程度或方向的一种方法或手段。相关分析的结果是相关(Correlation)或相关系数(Correlation Coefficients)，二者基本是等同的，主要用于确定两个随机变量之间线性关系的强度和方向。
利用相关分析进行时移量估计的基本原理是，在可能的时移范围内，将目标波形与某一时移量的参考波形进行互相关而获得相关系数。相关系数大表示波形相似度大(匹配程度高)，反之则相似度小(匹配程度低)。一般地，假设两个地震波形连续且呈线性关系，在最佳匹配点的两侧，相关系数会逐渐减小，呈现高斯分布特征。当求取了所有时移量所对应的相关系数后，从中搜索相关系数最大的匹配点的时移量作为目标波形相对于参考波形的倾角时差。
常用的互相关算法有：标准互相关和广义互相关。无论是标准互相关，还是广义互相关，无论是在时间域还是频率域进行互相关，其最高分辨率均无法超出采样间隔的限制，也就是说，利用互相关所求取的时间延迟精度必然小于采样间隔，这一特点会严重影响地震波形匹配的精度。因此，通常情况下要想实现超采样地震波形匹配，必须对地震波形进行加密采样，这显然会导致波形匹配效率严重降低。
为了提高波形匹配精度，现有技术中，利用多项式拟合对相关系数序列进行加密采样，然后从中搜索相关系数最大值所对应的时移量。这一技术可提高时移量估计的精度，但受地层旋回现象和地震数据信噪比的影响，存在两个显著的缺点：首先，由于地层旋回的存在，窄频带的地震波形存在着“特征循环”现象，即来自不同地层的地震响应会出现相同或相似的特征，如果时移量搜索范围过大，这一特点常常会导致地震波形匹配错误；其次，受噪音影响，相关系数变化剧烈，导致相关系数曲线的高斯分布特征发生畸变，多项式拟合会偏离正态分布，因而会造成时移量估算错误。
此外，由于标准傅里叶快速变换(FFT)的频谱存在栅栏效应，在现有技术中，采用Z变换对频谱进行细化、基于高阶累积量对频谱进行变换等算法试图提高时延估算精度，尽管这些算法在一定程度上提高了时移量估算精度，但并未解决由于窄频带地震波形“特征循环”现象所带来的波形匹配错误这一问题，因而其应用效果不够理想。
针对上述问题，目前尚未提出有效的解决方案。
发明内容
本发明实施例中提供了一种超采样地震波形匹配方法，用以解决现有技术中利用相关法进行地震波形匹配方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题。该方法包括：获取待匹配的目标地震波形和参考地震波形；循环执行如下步骤，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差：根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；根据相关系数序列确定初始时移点和初始时移量；以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；计算拟合时移量与初始时移量的拟合误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量。
本发明实施例中还提供了一种超采样地震波形匹配装置，用以解决现有技术中利用相关法进行地震波形匹配的方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题。该装置包括：波形提取模块，用于获取待匹配的目标地震波形和参考地震波形；波形匹配时差确定模块，用于循环执行波形尺度变换模块、相关系数序列生成模块、初始时移量确定模块、相关系数拐点确定模块、高斯曲线拟合模块和参数调整模块的功能，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差；波形尺度变换模块，用于根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；相关系数序列生成模块，用于遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；初始时移量确定模块，用于根据相关系数序列确定初始时移点和初始时移量；相关系数拐点确定模块，用于以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；高斯曲线拟合模块，用于对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；参数调整模块，用于计算拟合时移量与初始时移量的拟合误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量。
本发明实施例中还提供了一种计算机设备，用以解决现有技术中利用相关法进行地震波形匹配的方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题。该计算设备包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，该处理器执行所述计算机程序时实现上述超采样地震波形匹配方法。
本发明实施例中还提供了一种计算机可读存储介质，用以解决现有技术中利用相关法进行地震波形匹配的方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题。该计算机可读存储介质存储有执行上述超采样地震波形匹配方法的计算机程序。
本发明实施例中，在获取到待匹配的目标地震波形和参考地震波形后，根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；遍历时移量搜索范围内的多个时移量，对参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列，以确定初始时移点和初始时移量；进而以初始时移点为中心，从相关系数序列中识别出两个拐点，并对两个拐点之间的相关系数序列进行高斯曲线拟合，以确定拟合时移量，最后判断拟合时移量与初始时移量的误差绝对值是否小于预设阈值，当拟合时移量与初始时移量的误差绝对值大于或等于预设阈值的情况下，重新配置对参考地震波形和目标地震波形进行尺度变换的波形尺度变换指数，重复执行上述步骤，直到拟合时移量与初始时移量的误差绝对值小于预设阈值的拟合时移量，将此事的拟合时移量作为目标地震波形相对于参考地震波形的波形匹配时差。
本发明实施例中，将地震波形多尺度变换、自适应时移量搜索范围约束和高斯曲线拟合引入地震波形匹配中，突破了利用相关分析所估算的时移量精度无法突破采样间隔的限制，显著提高了地震波形匹配的精度，增强了地震波形匹配的稳定性，能够满足基于地震层位自动追踪、地震地层体追踪、地震属性分析等地震数据处理解释要求。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本发明实施例中提供的一种超采样地震波形匹配方法流程图；
图2为本发明实施例中提供的一种超采样地震波形匹配方法的具体实现流程图；
图3为本发明实施例中提供的目标地震波形与参考地震波形的对比示意图；
图4为现有技术中提供的地震波形匹配原理示意图；
图5为本发明实施例中提供的超采样地震波形匹配原理示意图；
图6为现有技术中提供的地震波形匹配方法的地震层位自动追踪效果示意图；
图7为本发明实施例中提供的超采样地震波形匹配方法的地震层位自动追踪效果示意图；
图8为本发明实施例中提供的一种超采样地震波形匹配装置示意图；
图9为本发明实施例中提供的一种可选的超采样地震波形匹配装置示意图；
图10为本发明实施例中提供的一种计算机设备示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
本发明实施例中提供了一种超采样地震波形匹配方法，图1为本发明实施例中提供的一种超采样地震波形匹配方法流程图，如图1所示，该方法包括如下步骤：
S101，获取待匹配的目标地震波形和参考地震波形。
在具体实施时，可以根据待匹配点位置、匹配半径和时移量搜索范围从参考地震道和目标地震道中分别抽取参考地震波形和目标地震波形。在一个实施例中，上述S101可以通过如下步骤来实现：获取配置的待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2；根据待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2，从参考地震道中提取坐标范围为t0-r+l1～t0+r+l2、波形长度为2r+1+|l1|+|l2|的参考地震波形，并从目标地震道中提取坐标范围为t0-r～t0+r、波形长度为2r+1的目标地震波形。
S102，循环执行S1021～S1026，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差。
需要说明的是，上述S102中的初始时移量是根据目标地震波形与参考地震波形在时移量搜索范围内的相关系数序列确定的一个时移量(通常为相关系数最大值对应的时移量)。
需要注意的是，现有技术中基于相关法的地震波形匹配方法，直接根据初始时移量来确定目标地震模型相对于参考地震模型的波形匹配误差，导致时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误。本发明实施例中提供的超采样地震波形匹配方法，通过执行S1021～S1026，能够将地震波形多尺度变换、自适应时移量搜索范围约束和高斯曲线拟合引入地震波形匹配中，能够显著提高地震波形匹配的精度，增强地震波形匹配的稳定性。
由于拟合时移量与初始时移量的拟合误差绝对值大于或等于预设阈值，表明高斯拟合结果是不合理的，因而，需要对目标地震波形和参考地震模型进行尺度变换，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值。在具体实施时，本发明实施例中的预设阈值可以是1。
S1021，根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；
在一个实施例中，上述S1021可以通过如下步骤来实现：采用正弦插值算法，对参考地震波形和目标地震波形进行尺度变换，使得尺度变换后的参考地震波形和目标地震波形的波形长度为变换前波形长度的2k倍，其中，k表示波形尺度变换指数。
S1022，遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列。
在一个实施例中，上述S1022可以通过如下步骤来实现：遍历时移量搜索范围内的每个时移量，从参考地震波形上提取与目标地震波形的波形长度相等的地震波形作为待匹配波形，对待匹配波形和目标地震波形进行互相关运算，得到一个时移量对应的相关系数；根据时移量搜索范围内所有时移量对应的相关系数，生成一个相关系数序列。
S1023，根据相关系数序列确定初始时移点和初始时移量。
在一个实施例中，上述S1023可以通过如下步骤来实现：将相关系数序列中相关系数最大值对应的时移点确定为初始时移点；将初始时移点对应的时移量确定为初始时移量。
S1024，以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；
具体地，上述S1024可以通过如下步骤来实现：获取相关系数序列对应的曲线段。根据二阶导数的变化特征，从曲线段中识别出第一拐点和第二拐点，其中，第一拐点位于初始时移点的前面，第二拐点位于初始时移点的后面。
S1025，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量。
具体地，上述S1024可以通过如下步骤来实现：采用最小二乘法，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定高斯函数的特征参数，其中，特征参数包括：峰值位置、峰值和半宽度；根据高斯函数的特征参数，将峰值位置与中心位置相减，得到拟合时移点；将拟合时移点对应的时移量确定为拟合时移量。
S1026，计算拟合时移量与初始时移量的拟合误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量。
在具体实施时，本发明实施例中提供的超采样地震波形匹配方法，可以将波形尺度变换指数初始化为0，当拟合时移量与初始时移量的拟合误差绝对值大于或等于预设阈值时，将波形尺度变换指数从0开始依次增加，每次增加1，实现对参考地震波形和目标地震波形的多尺度变换。
由上可知，本发明实施例中提供的超采样地震波形匹配方法，在获取到待匹配的目标地震波形和参考地震波形后，根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；遍历时移量搜索范围内的多个时移量，对参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列，以确定初始时移点和初始时移量；进而以初始时移点为中心，从相关系数序列中识别出两个拐点，并对两个拐点之间的相关系数序列进行高斯曲线拟合，以确定拟合时移量，最后判断拟合时移量与初始时移量的误差绝对值是否小于预设阈值，当拟合时移量与初始时移量的误差绝对值大于或等于预设阈值的情况下，重新配置对参考地震波形和目标地震波形进行尺度变换的波形尺度变换指数，重复执行上述步骤，直到拟合时移量与初始时移量的误差绝对值小于预设阈值的拟合时移量，将此事的拟合时移量作为目标地震波形相对于参考地震波形的波形匹配时差。
通过本发明实施例中提供的超采样地震波形匹配方法，能够将地震波形多尺度变换、自适应时移量搜索范围约束和高斯曲线拟合引入地震波形匹配中，突破了利用相关分析所估算的时移量精度无法突破采样间隔的限制，显著提高了地震波形匹配的精度，增强了地震波形匹配的稳定性，能够满足基于地震层位自动追踪、地震地层体追踪、地震属性分析等地震数据处理解释要求。
图2为本发明实施例中提供的一种超采样地震波形匹配方法的具体实现流程图，如图2所示，具体包括：
S201，从地震道中提取目标地震波形和参考地震波形。其中，参考地震道是指作为匹配基准的地震道；目标地震道是指待匹配的地震道；地震波形匹配，是指计算或求取目标地震道相对于参考地震道在某一时间点处(即待匹配点)的时移量。地震波形匹配的衡量标准是，当将目标地震波形按该时移量反向位移后，目标地震波形与参考地震波形的相似程度最高(在视觉上表现为两个地震波形沿时间轴对齐)；待匹配点是该地震道上的一个采样点；匹配半径是以待匹配点为中心的一个以采样点为单位的数值，用于确定参与互相关分析的地震波形长度(设匹配半径为r，则参与互相关分析的地震波形长度为2r+1)；时移量搜索范围，是指目标地震道相对于参考地震道的最佳时移量的最大分布范围；根据待匹配点位置、匹配半径和时移量搜索范围从参考地震道和目标地震道中分别抽取两个地震波形，是指分别从的两个地震道中提取出两个长度不一样的地震波形。设匹配点为t0，匹配半径为r，时移搜索范围为l1～l2，则从参考地震道中所提取的地震波形坐标范围为t0-r+l1～t0+r+l2，其长度为2r+1+|l1|+|l2|；从目标地震道中所提取的地震波形坐标范围为t0-r～t0+r，其长度为2r+1。
S202，设置初始地震波形变换指数。其中，地震波形尺度变换指数，是指对地震波形进行尺度变换时的一个地震道重采样指标。地震波形尺度变换指数是以2为底的指数，设地震波形尺度变换指数为k，则尺度变换后的地震波形的长度是原地震波形长度的2k倍；将地震波形尺度变换指数初始化为0，是指设置初始的地震波形尺度变换指数为0，即不进行尺度变换，或者说尺度变换后的地震波形长度为原地震波形长度的2°倍(其值为1，即与原地震波形长度一致)。
S203，对地震波形进行尺度变换。其中，尺度变换，是指采用正弦插值算法对地震波形进行插值。所谓正弦插值，是一种经典的用于波形类信号的重采样算法，可表示为：
其中，x(t)为原始信号，y(t)为插值后的信号，M为参与插值的采样点数，T为采样间隔。
需要注意的是，当尺度变换指数不为0时，要同时对目标地震波形和参考地震波形进行尺度变换。
S204，对待匹配的两个地震波形进行互相关分析。其中，互相关分析，既可以是标准互相关，也可以是广义互相关，既可以在时间域进行计算，也可以在频率域计算。
对于连续信号，标准互相关可表示为
其中，x和y为待匹配的地震波形，τ为时间延迟量，N为采样点数。
对于离散信号，标准互相关可表示为：
其中，x和y为待匹配的地震波形，
和
分别为x和y的平均值，N为采样点数。
广义互相关又称归一化互相关，是指在频域中利用窗函数对信号x和y的傅里叶变换进行加权而得到的，即：
其中，Gxx(τ)为信号x的自功率谱，Gyy(τ)为信号y的自功率谱，Gxy(τ)为信号x和y的互功率谱。
其中，根据时移搜索范围内所有的时移量，对两个地震波形进行互相关分析，形成相关系数序列，包含3个步骤：(1)针对时移搜索范围内的某一个时移量，根据该时移量从参考地震波形中抽取一个与目标地震波形长度相等的地震波形，形成待匹配波形；(2)对目标地震波形和待匹配波形进行互相关运算，获得一个相关系数；(3)将所有时移量所对应的相关系数组成相关系数序列。若时移搜索范围为l1～l2，则相关系数序列的长度为|l1|+|l2|+1。
S205，确定初始时移量。即从相关系数序列中搜索相关系数最大值及其对应的时移点，并将该时移点所对应的时移量标记为初始时移量。
受相关分析算法的限制，所获得的时移量必然是信号采样间隔的整数倍，因而其最高分辨率无法超过采样间隔。
在传统的基于相关分析的时移量估算中，到这一步就已经达到目的，但在本专利中，为了提高时移量估算精度，此时所获得的时移量只是一个初始的时移量，后续还要进一步处理。
S206，搜索相关系数曲线的拐点。其中，拐点，又称反曲点，是指曲率发生变化的点(即连续曲线的凹弧与凸弧的分界点)。若该曲线所对应的函数在拐点处有二阶导数，则二阶导数在拐点处异号(由正变负或由负变正)；以初始时移点为中心，从相关系数系列中寻找两个拐点，是指将相关系数序列视作一个曲线段，根据二阶导数的变化特征从该曲线段中识别出两个拐点，其中一个拐点位于初始时移点的前面，另一个拐点位于初始时移点的后面。
S207，对相关系数序列进行高斯曲线拟合。其中，对两个拐点间的相关系数序列进行高斯曲线拟合，是指将两个拐点间的相关系数序列视作一个高斯曲线，利用最小二乘法对其进行高斯曲线拟合，获得该高斯曲线的拟合参数(峰振幅、峰位置和半宽度)。
需要说明的是，本发明实施例中采用高斯曲线拟合算法对相关系数序列进行拟合，是因为相关系数序列分布具有典型的高斯分布特征，但受地震信号信噪比的影响，高斯分布特征常常会受到污染和扭曲变形。采用高斯曲线拟合对相关系数序列进行拟合的目的是将相关系数序列拟合为标准正态分布的高斯曲线。
高斯拟合，是指一种以高斯函数为目标对曲线型数据点集进行函数逼近的数据拟合方法。设有一组数据(xi，yi)，i＝0，1，2...n,将其用高斯函数描述：
其中，xp、yp和S分别为待估算的高斯曲线的峰值位置、峰值和半宽度(尺度)信息。
对上式两边取自然对数，可表示为
令
lnyi＝zi (7)
并考虑所有待拟合的数据点，则公式(5)可表示为：
进一步地，上式可简记为
Z＝XB (12)
根据最小二乘法原理，公式(11)的解为
B＝(XTX)-1XTZ (13)
根据公式(7)～公式(10)，可计算出高斯函数的特征参数xp、yp和S。
S208，确定拟合时移量和拟合时移点。根据公式(4)，xp即为待估算的高斯曲线的峰值位置，将其与高斯函数的中心位置相减，即可获得拟合时移点，并将拟合时移点所对应的时移量标记为拟合时移量。
S209，计算时移量拟合误差。拟合误差为拟合时移量与初始时移量之差的绝对值。
S210，判断拟合误差是否符合要求。若拟合误差小于1.0，表明高斯拟合结果是合理的，执行S211。若拟合误差大于或等于1.0，将地震波形尺度变换指数增加1，根据拐点更新时移搜索范围，重复S203～S209，是指当拟合时移量与初始时移量的误差绝对值大于或等于1.0时，表明高斯拟合结果是不合理的，此时需要对目标波形和参考波形进行尺度变换，然后重新进行分析和计算，直至拟合时移量与初始时移量的误差绝对值符合要求并退出循环。
S211，将拟合时移点记为波形匹配点，根据拟合时移量确定波形匹配时差，将所对应的拟合时移量可作为目标波形相对于参考波形的匹配结果输出。
下面结合图3～图5，对本发明实施例中提供的超采样地震波形匹配方法的具体应用进行详细说明。
图3示出了待匹配的两个地震波形，其中，图3中(a)所示为目标地震波形，图3中(b)所示为参考地震波形。在具体实施时，采用的匹配半径为6、时移搜索半径为6。其中，图3中(a)所示的目标地震波形包含13个采样点(其采样点序号为7-19)，其中，第13号采样点为该参考地震波形的中心点，且被视作是该波形的基准匹配点；图3中(a)所示的参考地震波形包含25个采样点(其采样点序号为1-25)，其中，第13号采样点为该参考地震波形的中心点，且被视作是该波形的基准匹配点。
图4示出了现有技术中提供的基于相关分析算法的地震波形匹配原理图。该地震波形匹配的目标是，以参考地震波形作为模板，从目标地震波形中寻找与之相关性最大且长度相同的波形段(该波形段与参考地震波形具有最大相关系数)，并将其称之为匹配波形。当确定了匹配波形后，该匹配波形的中心点得以确定。将该匹配波形的中心点序号减去目标地震波形的中心点序号，即可获得目标地震波形和参考地震波形之间的时移量。要实现利用模板从地震波形中寻找最佳匹配波形的目的，需要沿时间方向向前和向后滑动参考地震波形，并计算滑动后的目标地震波形与参考地震波形的相关系数。显然，在滑动目标地震波形时，滑动步长不能小于采样间隔，这就是为什么相关分析法时移量估计精度无法超出采样间隔的原因。
在相关分析中，每一次滑动进行一次相关分析并得到一个相关系数，因此当时移搜索半径为6时，所得到的相关系数数量为13。将这些相关系数数量组合为一个相关系数序列，可形成一条如图4中(a)所示的曲线。一般地，该曲线呈正态分布形态，可用高斯函数进行描述。
从相关系数序列中搜索相关系数最大值，可获得最佳匹配点，如图4中(a)所示的t′0。如果两个地震波形完全相同，则时移量必然为0，则其最佳匹配点为时移量为0的点，如图4中(a)所示的t0，如果两个波形存在时差，则时移量不为0。基于这一特点可利用相关分析获得两个地震波形之间的时移量(t′0-t0)。在本实施例中，目标地震波形相对于参考地震波形的时移量为-1.0，也就是说，当将目标地震波形向上移动1个采样点的距离时，目标地震波形与参考地震波形的相关性最大，如图4中(b)所示。
在时移相关分析中，在滑动参考地震波形时，滑动步长不能小于采样间隔，这就是为什么相关分析法时移量估计精度无法超出采样间隔的原因。由于这一原因，相关分析法时移量估计精度无法满足要求，要提高时移量估计精度，只能是对地震信号进行加密采样，但这显然会导致地震波形匹配效率降低。
为了实现超采样地震波形匹配，在本发明实施例中，采用高斯曲线拟合算法对相关系数曲线进行拟合，并根据高斯曲线的性质求取高斯曲线的极值点及其对应的时移量，并将这种方式所获得的时移量称为拟合时移量。图5示出了利用高斯拟合求取拟合时移量的原理，横坐标轴为时移量，纵坐标轴为相关系数，对13个点的原始相关系数曲线进行高斯拟合，可获得一条光滑的高斯曲线(拟合相关系数曲线)。根据公式(7)-公式(10)，可求出该高斯曲线的峰值位置(拟合时移量，图5中的Δt2)。在本实施例中，拟合时移量为-1.392339，而通过相关分析所获得的初始时移量为-1.0(图5中的Δt1)，显然，拟合时移量的精度显著高于初始时移量。
需要说明的是，高斯曲线拟合算法依赖于相关系数曲线的正态分布特征。当时移搜索半径过大、采样精度不足或地震信号信噪比降低到一定程度时，相关系数曲线会出现多个峰值或偏离正态分布形态，此时高斯曲线拟合会出现严重的误差。为避免这一现象的出现，采用拟相关系数拐点搜索模块和合误差计算与分析模块对地震波形匹配过程进行控制。首先，利用以初始时移点为中心，从相关系数序列中寻找2个分别位于时移点两侧的拐点，仅对两个拐点间的相关系数序列进行高斯曲线拟合，这会避免相关系数序列出现多个峰值的现象；其次，计算拟合误差(拟合时移量与初始时移量之差的绝对值)，然后对拟合误差进行分析。如果拟合误差小于1.0，表明相关分析结果是稳定的，此时可转入时移量输出模块并终止地震波形匹配过程；如果拟合误差大于或等于1.0，表明受地震波形采样间隔过大影响而导致相关系数序列出现相关系数失真现象，此时需转入尺度变换指数初始化与调整模块，并开始实施一个新的地震波形匹配过程。在本实施例中，由于时移搜索半径合理、采样精度较高、地震信号信噪比较高，拟合误差为0.392339，满足地震匹配过程终止条件，因而无需修改波形尺度变换指数，可将拟合时移量-1.392339输出并退出地震波形匹配过程。
图6为现有技术中提供的地震波形匹配方法的地震层位自动追踪效果示意图，如图6所示，采用经典的基于相关分析的时移量估计方法，以一个种子点为起点向两侧追踪，所获得的地震层位在某些位置上显著地偏离了地震同相轴所指示的位置，且离种子点越远，这种偏离现象越加严重。
图7为本发明实施例中提供的超采样地震波形匹配方法的地震层位自动追踪效果示意图，如图7所示，采用超采样地震波形匹配方法，以与图6所示的同样的种子点为起点向两侧追踪，所获得的地震层位与地震同相轴所指示的形态一致，层位追踪效果非常理想。
基于同一发明构思，本发明实施例中还提供了一种超采样地震波形匹配装置，如下面的实施例所述。由于该装置解决问题的原理与超采样地震波形匹配方法相似，因此该装置的实施可以参见超采样地震波形匹配方法的实施，重复之处不再赘述。
图8为本发明实施例中提供的一种超采样地震波形匹配装置示意图，如图8所示，该装置包括：波形提取模块81、波形匹配时差确定模块82、波形尺度变换模块821、相关系数序列生成模块822、初始时移量确定模块823、相关系数拐点确定模块824、高斯曲线拟合模块825和参数调整模块的功能826；
其中，波形提取模块81，用于获取待匹配的目标地震波形和参考地震波形；波形匹配时差确定模块82，用于循环执行波形尺度变换模块821、相关系数序列生成模块822、初始时移量确定模块823、相关系数拐点确定模块824、高斯曲线拟合模块825和参数调整模块的功能826，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差；波形尺度变换模块821，用于根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；相关系数序列生成模块822，用于遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；初始时移量确定模块823，用于根据相关系数序列确定初始时移点和初始时移量；相关系数拐点确定模块824，用于以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；高斯曲线拟合模块825，用于对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；参数调整模块826，用于计算拟合时移量与初始时移量的拟合误差绝对值，当拟合误差绝对值大于或等于预设阈值的情况下，重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围。
由上可知，本发明实施例中提供的超采样地震波形匹配装置，通过波形提取模块81获取待匹配的目标地震波形和参考地震波形；通过波形尺度变换模块821根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；通过相关系数序列生成模块822遍历时移量搜索范围内的多个时移量，对参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；通过初始时移量确定模块823根据相关系数序列确定初始时移点和初始时移量；通过相关系数拐点确定模块824以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；通过高斯曲线拟合模块825对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；通过参数调整模块826计算拟合时移量与初始时移量的拟合误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量；通过波形匹配时差确定模块82循环执行波形尺度变换模块821、相关系数序列生成模块822、初始时移量确定模块823、相关系数拐点确定模块824、高斯曲线拟合模块825和参数调整模块的功能826，直到拟合时移量与初始时移量的拟合误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差。
通过本发明实施例中提供的超采样地震波形匹配装置，能够将地震波形多尺度变换、自适应时移量搜索范围约束和高斯曲线拟合引入地震波形匹配中，突破了利用相关分析所估算的时移量精度无法突破采样间隔的限制，显著提高了地震波形匹配的精度，增强了地震波形匹配的稳定性，能够满足基于地震层位自动追踪、地震地层体追踪、地震属性分析等地震数据处理解释要求。
在一个实施例中，本发明实施例中提供的超采样地震波形匹配装置还可以包括：参数配置模块83，用于获取配置的待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2；其中，波形提取模块81还用于根据待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2，从参考地震道中提取坐标范围为t0-r+l1～t0+r+l2、波形长度为2r+1+|l1|+|l2|的参考地震波形，并从目标地震道中提取坐标范围为t0-r～t0+r、波形长度为2r+1的目标地震波形。
在一个实施例中，上述相关系数序列生成模块822还用于遍历时移量搜索范围内的每个时移量，从参考地震波形上提取与目标地震波形的波形长度相等的地震波形作为待匹配波形，对待匹配波形和目标地震波形进行互相关运算，得到一个时移量对应的相关系数；以及根据时移量搜索范围内所有时移量对应的相关系数，生成一个相关系数序列。
在一个实施例中，上述初始时移量确定模块823还用于将相关系数序列中相关系数最大值对应的时移点确定为初始时移点；以及将初始时移点对应的时移量确定为初始时移量。
在一个实施例中，上述波形尺度变换模块821还用于采用正弦插值算法，对参考地震波形和目标地震波形进行尺度变换，使得尺度变换后的参考地震波形和目标地震波形的波形长度为变换前波形长度的2k倍，其中，k表示波形尺度变换指数。
在一个实施例中，上述相关系数拐点确定模块824用于获取相关系数序列对应的曲线段；以及根据二阶导数的变化特征，从曲线段中识别出第一拐点和第二拐点，其中，第一拐点位于初始时移点的前面，第二拐点位于初始时移点的后面。
在一个实施例中，上述高斯曲线拟合模块825还用于采用最小二乘法，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定高斯函数的特征参数；根据高斯函数的特征参数，将峰值位置与中心位置相减，得到拟合时移点；将拟合时移点对应的时移量确定为拟合时移量，其中，特征参数包括：峰值位置、峰值和半宽度。
图9为本发明实施例中提供的一种可选的超采样地震波形匹配装置示意图，该装置与图2所示的方法对应，如图9所示，该装置包括：地震波形输入模块901、尺度变换指数初始化与调整模块902、地震波形尺度变换模块903、地震波形相关分析模块904、初始时移量估算模块905、相关系数拐点搜索模块906、高斯曲线拟合模块907、拟合时移量确定模块908、拟合误差计算与分析模块909和时移量输出模块910。
其中，地震波形输入模块901，用于从计算机文件或数据库中读入参考地震道和目标地震道的地震波形数据，形成两个地震波形序列；地震波形尺度变换指数初始化与调整模块902，用于对地震波形尺度变换中的尺度变换指数进行初始化，并在每一次循环中将初始化指数递增，递增量为1；地震波形尺度变换模块903，用于根据尺度变换指数对地震波形进行重采样，形成新的地震波形序列；地震波形相关分析模块904，用于对参考地震道和目标地震道进行滑动相关分析，形成相关系数序列；初始时移量估算模块905，用于从相关系数序列中搜索相关系数最大处的时移点，并将该时移点命名为初始时移点，将其所对应的时移量命名为初始时移量；相关系数拐点搜索模块906，用于以初始时移点为中心，从相关系数序列中寻找两个分别位于时移点两侧的拐点(即第一拐点和第二拐点)；高斯曲线拟合模块907，用于对两个拐点间的相关系数序列进行高斯曲线拟合，形成高斯拟合系数序列；拟合时移量确定模块908，用于根据高斯拟合系数序列求取相关系数曲线的极大值点及其对应的时移量，并将其命名为拟合时移量；拟合误差计算与分析模块909，用于计算拟合误差，并根据拟合误差的大小确定是否需要提高地震波形尺度变换指数并开始一个新的地震波形匹配过程；时移量输出模块910，用于输出符合精度要求的时移量，并终止地震波形匹配过程。
需要说明的是，本发明实施例中拟合误差，是拟合时移量与初始时移量之差的绝对值。如果拟合误差小于1.0，表明相关分析结果是稳定的，此时可转入时移量输出模块并终止地震波形匹配过程；如果拟合误差大于或等于1.0，表明受地震波形采样间隔过大影响而导致相关系数序列出现相关系数失真现象，此时需转入尺度变换指数初始化与调整模块，并开始实施一个新的地震波形匹配过程。
基于同一发明构思，本发明实施例中还提供了一种计算机设备，用以解决现有技术中利用相关法进行地震波形匹配的方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题，图10为本发明实施例中提供的一种计算机设备示意图，如图10所示，该计算设备100包括存储器1001、处理器1002及存储在存储器1001上并可在处理器1002上运行的计算机程序，所述处理器1002执行所述计算机程序时实现上述超采样地震波形匹配方法。
基于同一发明构思，本发明实施例中还提供了一种计算机可读存储介质，用以解决现有技术中利用相关法进行地震波形匹配的方法所导致的时移量估算精度无法超出采样间隔的限制，容易出现波形匹配错误的技术问题，该计算机可读存储介质存储有执行上述超采样地震波形匹配方法的计算机程序。
综上所述，本发明实施例中提供了一种超采样地震波形匹配方法、装置、计算机设备及计算机可读存储介质，能够将地震波形多尺度变换、自适应时移量搜索范围约束和高斯曲线拟合引入地震波形匹配中，突破了利用相关分析所估算的时移量精度无法突破采样间隔的限制，显著提高了地震波形匹配的精度，增强了地震波形匹配的稳定性，能够满足基于地震层位自动追踪、地震地层体追踪、地震属性分析等地震数据处理解释要求。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (16)
1.一种超采样地震波形匹配方法，其特征在于，包括：
获取待匹配的目标地震波形和参考地震波形；
循环执行如下步骤，直到拟合时移量与初始时移量的误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差：
根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；
遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；
根据相关系数序列确定初始时移点和初始时移量；
以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；
对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；
计算拟合时移量与初始时移量的误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量。
2.如权利要求1所述的方法，其特征在于，获取待匹配的目标地震波形和参考地震波形，包括：
获取配置的待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2；
根据待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2，从参考地震道中提取坐标范围为t0-r+l1～t0+r+l2、波形长度为2r+1+|l1|+|l2|的参考地震波形，并从目标地震道中提取坐标范围为t0-r～t0+r、波形长度为2r+1的目标地震波形。
3.如权利要求2所述的方法，其特征在于，遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列，包括：
遍历时移量搜索范围内的每个时移量，从参考地震波形上提取与目标地震波形的波形长度相等的地震波形作为待匹配波形，对待匹配波形和目标地震波形进行互相关运算，得到一个时移量对应的相关系数；
根据时移量搜索范围内所有时移量对应的相关系数，生成一个相关系数序列。
4.如权利要求1所述的方法，其特征在于，根据相关系数序列确定初始时移点和初始时移量，包括：
将相关系数序列中相关系数最大值对应的时移点确定为初始时移点；
将初始时移点对应的时移量确定为初始时移量。
5.如权利要求1所述的方法，其特征在于，根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换，包括：
采用正弦插值算法，对参考地震波形和目标地震波形进行尺度变换，使得尺度变换后的参考地震波形和目标地震波形的波形长度为变换前波形长度的2k倍，其中，k表示波形尺度变换指数。
6.如权利要求1所述的方法，其特征在于，以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点，包括：
获取相关系数序列对应的曲线段；
根据二阶导数的变化特征，从所述曲线段中识别出第一拐点和第二拐点，其中，所述第一拐点位于初始时移点的前面，所述第二拐点位于初始时移点的后面。
7.如权利要求1所述的方法，其特征在于，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量，包括：
采用最小二乘法，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定高斯函数的特征参数，其中，所述特征参数包括：峰值位置、峰值和半宽度；
根据高斯函数的特征参数，将峰值位置与中心位置相减，得到拟合时移点；
将拟合时移点对应的时移量确定为拟合时移量。
8.一种超采样地震波形匹配装置，其特征在于，包括：
波形提取模块，用于获取待匹配的目标地震波形和参考地震波形；
波形匹配时差确定模块，用于循环执行波形尺度变换模块、相关系数序列生成模块、初始时移量确定模块、相关系数拐点确定模块、高斯曲线拟合模块和参数调整模块的功能，直到拟合时移量与初始时移量的误差绝对值小于预设阈值，根据拟合时移量确定目标地震波形相对于参考地震波形的波形匹配时差；
波形尺度变换模块，用于根据配置的波形尺度变换指数，对参考地震波形和目标地震波形进行尺度变换；
相关系数序列生成模块，用于遍历时移量搜索范围内的每个时移量，对尺度变换后的参考地震波形和目标地震波形进行互相关分析，生成一个相关系数序列；
初始时移量确定模块，用于根据相关系数序列确定初始时移点和初始时移量；
相关系数拐点确定模块，用于以初始时移点为中心，从相关系数序列中识别出第一拐点和第二拐点；
高斯曲线拟合模块，用于对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定拟合时移点和拟合时移量；
参数调整模块，用于计算拟合时移量与初始时移量的误差绝对值，若误差绝对值小于预设阈值则退出循环；若误差绝对值大于或等于预设阈值，则重新配置一个波形尺度变换指数，根据第一拐点和第二拐点更新时移量搜索范围，重新计算拟合时移量。
9.如权利要求8所述的装置，其特征在于，所述装置还包括：参数配置模块，用于获取配置的待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2；
其中，所述波形提取模块还用于根据待匹配点位置t0、匹配半径r和时移量搜索范围l1～l2，从参考地震道中提取坐标范围为t0-r+l1～t0+r+l2、波形长度为2r+1+|l1|+|l2|的参考地震波形，并从目标地震道中提取坐标范围为t0-r～t0+r、波形长度为2r+1的目标地震波形。
10.如权利要求8所述的装置，其特征在于，所述相关系数序列生成模块还用于遍历时移量搜索范围内的每个时移量，从参考地震波形上提取与目标地震波形的波形长度相等的地震波形作为待匹配波形，对待匹配波形和目标地震波形进行互相关运算，得到一个时移量对应的相关系数；以及根据时移量搜索范围内所有时移量对应的相关系数，生成一个相关系数序列。
11.如权利要求8所述的装置，其特征在于，所述初始时移量确定模块还用于将相关系数序列中相关系数最大值对应的时移点确定为初始时移点；以及将初始时移点对应的时移量确定为初始时移量。
12.如权利要求8所述的装置，其特征在于，所述波形尺度变换模块还用于采用正弦插值算法，对参考地震波形和目标地震波形进行尺度变换，使得尺度变换后的参考地震波形和目标地震波形的波形长度为变换前波形长度的2k倍，其中，k表示波形尺度变换指数。
13.如权利要求8所述的装置，其特征在于，所述相关系数拐点确定模块用于获取相关系数序列对应的曲线段；以及根据二阶导数的变化特征，从所述曲线段中识别出第一拐点和第二拐点，其中，所述第一拐点位于初始时移点的前面，所述第二拐点位于初始时移点的后面。
14.如权利要求8所述的装置，其特征在于，所述高斯曲线拟合模块还用于采用最小二乘法，对第一拐点和第二拐点之间的相关系数序列进行高斯曲线拟合，确定高斯函数的特征参数；根据高斯函数的特征参数，将峰值位置与中心位置相减，得到拟合时移点；将拟合时移点对应的时移量确定为拟合时移量，其中，所述特征参数包括：峰值位置、峰值和半宽度。
15.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7任一项所述超采样地震波形匹配方法。
16.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至7任一项所述超采样地震波形匹配方法的计算机程序。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110074803.1A CN114859404A (zh) | 2021-01-20 | 2021-01-20 | 超采样地震波形匹配方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110074803.1A CN114859404A (zh) | 2021-01-20 | 2021-01-20 | 超采样地震波形匹配方法及装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114859404A true CN114859404A (zh) | 2022-08-05 |
Family
ID=82622868
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110074803.1A Pending CN114859404A (zh) | 2021-01-20 | 2021-01-20 | 超采样地震波形匹配方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114859404A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116187064A (zh) * | 2023-02-14 | 2023-05-30 | 中国科学院国家空间科学中心 | 一种连续信号时间序列二阶导数的数值仿真方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011005353A1 (en) * | 2009-07-06 | 2011-01-13 | Exxonmobil Upstream Research Company | Method for seismic interpretation using seismic texture attributes |
| US20120250748A1 (en) * | 2011-04-04 | 2012-10-04 | U.S. Government As Represented By The Secretary Of The Army | Apparatus and method for sampling and reconstruction of wide bandwidth signals below nyquist rate |
| US20160231357A1 (en) * | 2015-01-26 | 2016-08-11 | Guzik Technical Enterprises | Method and apparatus for data acquisition with waveform trigger |
| CN107329170A (zh) * | 2017-08-07 | 2017-11-07 | 中国石油天然气集团公司 | 一种地震波形匹配方法及装置 |
| CN107563297A (zh) * | 2017-08-07 | 2018-01-09 | 中国石油天然气集团公司 | 一种波形匹配方法及装置 |
| CN108957536A (zh) * | 2018-06-14 | 2018-12-07 | 中国石油天然气股份有限公司 | 基于波形匹配反演的地层倾角预测方法及装置 |
-
2021
- 2021-01-20 CN CN202110074803.1A patent/CN114859404A/zh active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011005353A1 (en) * | 2009-07-06 | 2011-01-13 | Exxonmobil Upstream Research Company | Method for seismic interpretation using seismic texture attributes |
| US20120250748A1 (en) * | 2011-04-04 | 2012-10-04 | U.S. Government As Represented By The Secretary Of The Army | Apparatus and method for sampling and reconstruction of wide bandwidth signals below nyquist rate |
| US20160231357A1 (en) * | 2015-01-26 | 2016-08-11 | Guzik Technical Enterprises | Method and apparatus for data acquisition with waveform trigger |
| CN107329170A (zh) * | 2017-08-07 | 2017-11-07 | 中国石油天然气集团公司 | 一种地震波形匹配方法及装置 |
| CN107563297A (zh) * | 2017-08-07 | 2018-01-09 | 中国石油天然气集团公司 | 一种波形匹配方法及装置 |
| CN108957536A (zh) * | 2018-06-14 | 2018-12-07 | 中国石油天然气股份有限公司 | 基于波形匹配反演的地层倾角预测方法及装置 |
Non-Patent Citations (4)
| Title |
|---|
| ANDREAS KJELSRUD EVENSEN 等: "Time-lapse Tomographic Inversion Using a Gaussian Parameterization of the Velocity Changes", GEOPHYSICS, vol. 75, no. 4, 10 September 2010 (2010-09-10), pages 29 * |
| 刘瑞林 等: "动态波形匹配预测火山岩地层的横向变化", 地球物理学进展, vol. 19, no. 04, 30 December 2004 (2004-12-30), pages 946 - 952 * |
| 杨学博 等: "大光斑LiDAR全波形数据小波变换的高斯递进分解", 红外与毫米波学报, vol. 36, no. 06, 31 December 2017 (2017-12-31), pages 749 - 755 * |
| 龚雪萍 等: "改进纵波与转换波时间匹配方法研究", 石油地球物理勘探, vol. 47, no. 05, 31 October 2012 (2012-10-31), pages 698 - 704 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116187064A (zh) * | 2023-02-14 | 2023-05-30 | 中国科学院国家空间科学中心 | 一种连续信号时间序列二阶导数的数值仿真方法 |
| CN116187064B (zh) * | 2023-02-14 | 2024-03-12 | 中国科学院国家空间科学中心 | 一种连续信号时间序列二阶导数的数值仿真方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10180515B2 (en) | Trace downsampling of distributed acoustic sensor data | |
| US11181653B2 (en) | Reservoir characterization utilizing ReSampled seismic data | |
| WO2004061727A1 (en) | Method of conditioning a random field to have directionally varying anisotropic continuity | |
| CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
| CN102707314A (zh) | 一种多路径双谱域混合相位子波反褶积方法 | |
| CN109143331B (zh) | 地震子波提取方法 | |
| CN105093315B (zh) | 一种去除煤层强反射信号的方法 | |
| WO2005017564A1 (en) | Dip value in seismic images | |
| CN104730576A (zh) | 基于Curvelet变换的地震信号去噪方法 | |
| CN109541691B (zh) | 一种地震速度反演方法 | |
| CN114859404A (zh) | 超采样地震波形匹配方法及装置 | |
| CN114966837A (zh) | 基于卷积型w2距离目标函数的全波形反演方法及装置 | |
| CN101825722B (zh) | 一种鲁棒的地震信号瞬时频率的估计方法 | |
| CN108169795A (zh) | 基于随机采样的数据规则化方法 | |
| CN114114421B (zh) | 基于深度学习的导向自学习地震数据去噪方法及装置 | |
| CN115685318A (zh) | 一种基于动态匹配的抗假频地震数据插值方法、电子设备及存储介质 | |
| AU2015309050B2 (en) | Joint estimation of electromagnetic earth responses and ambient noise | |
| CN110865409B (zh) | 一种基于波阻抗低秩正则化的地震波阻抗反演方法 | |
| CN111929726B (zh) | 地震相干数据体处理方法及装置 | |
| CN111352158B (zh) | 地震信号增强方法及装置 | |
| CN111965706B (zh) | 地震反演方法及装置 | |
| CN113126163B (zh) | 五维地震数据噪声衰减方法及装置 | |
| CN112558147B (zh) | 一种井中微地震数据的偏振分析方法及系统 | |
| CN113534250A (zh) | 一种基于快速匹配追踪的多尺度地震反演方法 | |
| Cai et al. | Multiscale dilated denoising convolution with channel attention mechanism for micro-seismic signal denoising |
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 |