CN108535777B - 一种地震波初至检测方法及系统 - Google Patents
一种地震波初至检测方法及系统 Download PDFInfo
- Publication number
- CN108535777B CN108535777B CN201810768659.XA CN201810768659A CN108535777B CN 108535777 B CN108535777 B CN 108535777B CN 201810768659 A CN201810768659 A CN 201810768659A CN 108535777 B CN108535777 B CN 108535777B
- Authority
- CN
- China
- Prior art keywords
- time window
- data
- value
- arrival
- earthquake record
- 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.)
- Expired - Fee Related
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 51
- 238000000034 method Methods 0.000 claims abstract description 48
- 238000005457 optimization Methods 0.000 claims abstract description 22
- 238000012417 linear regression Methods 0.000 claims abstract description 18
- 238000000546 chi-square test Methods 0.000 claims abstract description 12
- 239000013589 supplement Substances 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000001502 supplementing effect Effects 0.000 claims description 5
- 235000013399 edible fruits Nutrition 0.000 claims description 3
- 238000009826 distribution Methods 0.000 description 5
- 239000011159 matrix material Substances 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000001427 coherent effect Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
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/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
- G01V1/305—Travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/41—Arrival times, e.g. of P or S wave or first break
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
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
技术领域
本发明涉及地震波初至检测技术领域,特别是涉及一种地震波初至检测方法及系统。
背景技术
地震波初至在矿物资源开发、地震定位以及地震层析成像中都起着十分关键的作用。在石油勘探技术中,地震波信号可以用来了解地质构造,例如石油勘探就需要利用人工地震信号,其中准确检测到地震波到达检波点的初至时刻尤为重要。
从根本上说,地震信号的初至波既包含地震信号又包含噪声信号,是两者的叠加,初至拾取的精确性是地震资料处理中基础而又重要的环节。传统的初至拾取主要依靠手工,但人工拾取耗时且精确度低。为了提高地震波初至的提取效率与提取速度,现在已经出现了自动拾取的方法,例如能量比法、最大振幅法、分形维法等。
初至自动拾取方法中最常用的是基于能量比的方法。根据能量比达到最大值时的时间点就可以估算出地震波初至到来的时刻。但值得注意的是,当信噪比相对较低时这种基于能量比的方法就十分不理想。与能量比法原理相类似,振幅特性也常被用来检测初至,利用这种方法,根据地震波初至前后的振幅比的变化来拾取初至。一般来说,当信噪比足够高时这些方法是有效的。但我国石油勘探环境并不都是理想的,地表复杂的地区噪声干扰非常严重,这导致了现有方法对低信噪比情况下的初至拾取并不理想,精度并不高。
发明内容
本发明的目的是提供一种地震波初至检测方法及系统,有效避免噪声对检测精度的影响,提高检测准确度。
为实现上述目的,本发明提供了如下方案:
一种地震波初至检测方法,包括:
获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
获取检波器采集的待检测的地震记录序列;
建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗;
对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据;
从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达;
记录所有有地震波动到达的数据在所述地震记录序列中的序号,并计算每个有地震波动到达的数据对应的初至。
可选的,所述对所述地震记录序列的左侧和右侧进行数据填充,具体包括:
将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;
将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
可选的,所述每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达,具体包括:
获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;
求取所述第一时间窗内所有数据的最大值和最小值,得到取值区间;
将所述取值区间划分为K段,得到K段小区间;
计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k);其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;k为小区间序号,k=1,2,…,K;
计算数量差值比率:
根据参数K和置信度α,通过查询卡方表获取卡方值χ2(α,K-1);
判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;
若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,则确定在所述地震记录序列的第n个数据处没有地震波动到达。
可选的,利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练的过程为:
获取已确定初至的地震序列的实际初至值;
以子序列个数为横坐标,时间窗长度为纵坐标建立二维坐标系,根据时间窗长度的取值范围和子序列个数的取值范围确定所述二维坐标系内的取值区域;
在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,计算每个插值点所对应的初至预测值;
计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;
将所述取值区域进行网格划分,获取各个网格点作为被插值点;
利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值;
求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
本发明还公开一种地震波初至检测系统,包括:
时间窗最优长度获取模块,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
序列获取模块,用于获取检波器采集的地震记录序列;
时间窗建立模块,用于建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗;
数据填充模块,用于对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据;
地震波初至检测模块,用于从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达;
记录与初至计算模块,用于记录所有有地震波动到达的数据在所述地震记录序列中的序号,并计算每个有地震波动到达的数据对应的初至。
可选的,所述数据填充模块,具体包括:
序列左填充单元,用于将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;
序列右填充单元,用于将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
可选的,所述地震波初至检测模块,具体包括:
子序列最优个数获取单元,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;
区间计算单元,用于求取所述第一时间窗内所有数据的最大值和最小值,得到取值区间;
区间划分单元,用于将所述取值区间划分为K段,得到K段小区间;
区间数据量计算单元,用于计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k);其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;k为小区间序号,k=1,2,…,K;
比率计算单元,用于计算数量差值比率:
卡方值查询单元,用于根据参数K和置信度α,通过查询卡方表获取卡方值χ2(α,K-1);
判断单元,用于判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;
地震初至确定单元,用于若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,则确定在所述地震记录序列的第n个数据处没有地震波动到达。
可选的,该地震波初至检测系统还包括训练模块,所述训练模块用于利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练;所述训练模块包括:
已知序列获取单元,用于获取已确定初至的地震序列的实际初至值;
坐标系建立单元,用于以子序列个数为横坐标,时间窗长度为纵坐标建立二维坐标系,根据时间窗长度的取值范围和子序列个数的取值范围确定所述二维坐标系内的取值区域;
插值点初至预测计算单元,用于在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,计算每个插值点所对应的初至预测值;
预测误差计算单元,用于计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;
网格划分单元,用于将所述取值区域进行网格划分,获取各个网格点作为被插值点;
线性拟合单元,用于利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值;
最优插值点计算单元,用于求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
根据本发明提供的具体实施例,本发明公开了以下技术效果:本发明的地震波初至检测方法及系统,利用卡方检验的方法识别地震波初至,关注地震记录数据落入某个区间的数量,能够避免依赖于地震数据本身的大小,从而有效避免噪声的影响,提高检测准确度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明地震波初至检测方法实施例中训练过程的流程图;
图2为本发明地震波初至检测方法实施例中在坐标系内选取的插值点分布图;
图3为本发明地震波初至检测方法实施例中检测过程的流程图;
图4为本发明地震波初至检测系统实施例的系统结构图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种地震波初至检测方法及系统,在低信噪比的情况下有效避免噪声对检测精度的影响,提高检测准确度。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明将时间窗口内采样值落入某一区间的次数作为指标,统计次数差值之和,并将此差值之和与卡方值相比较,将比较结果作为地震初至是否到达的标志。
该地震波初至检测方法,包括:训练过程和检测过程。
图1为本发明地震波初至检测方法实施例中训练过程的流程图。
参见图1,训练过程如下:
步骤101:获取已确定初至的地震序列的实际初至值。已确定初至的地震序列表示为:NT为已确定初至的地震序列的长度。
步骤102:以子序列个数K′为横坐标,时间窗长度L'为纵坐标建立二维坐标系,根据时间窗长度L'的取值范围和子序列个数K′的取值范围确定所述二维坐标系内的取值区域。
子序列个数K′的取值范围为[1,51],时间窗长度L'的取值范围为[1,201]。且K′和L'均为整数,且L'为奇数。建立坐标系后,子序列个数K′和时间窗长度L'组成参数矢量
步骤103:在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,调用地震波初至检测算法(该地震波初至检测方法与检测过程的检测步骤类似),计算每个插值点所对应的初至预测值;
训练算法的实质是基于Kriging模型的数据插值算法。因此,首选需要明确已知的插值点及其所对应的函数值(本发明以每个插值点所对应的预测误差作为函数值)。以九个插值点为例,图2为本发明地震波初至检测方法实施例中在坐标系内选取的插值点分布图。参见图2,图2中的圆形坐标点即插值点,9个已知的插值点为:
其中为9个插值点的横坐标,即9个插值点对应的时间窗长度;为9个插值点的纵坐标,即9个插值点对应子序列个数。
步骤104:计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;
9个插值点的预测误差集合为:其中分别为9个插值点的预测误差。
步骤105:将所述取值区域进行网格划分,获取各个网格点作为被插值点。
参见图2,图2中方形坐标点为被插值点,图2中仅示意出四个被插值点。本发明的一个实施方式中,将取值区域的横纵坐标范围均划分成10个网格点。划分的点数越多,则最终的结果精确度越高。
步骤106:利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值。
求取拟合值的公式为:
其中为拟合值,C为拟合矩阵。拟合矩阵的计算公式为:
C=R-1[r-Fλ]
其中λ为中间参数,是与拟合变量空间位置有关的一个参数,λ=[FTR-1F]-1[FTR-1r-1],
其中F为回归矩阵,R为相关矩阵,r为另一个中间参数。具体为:
回归矩阵表示线性回归的系数。
相关矩阵R=[Rij]9×9,其中为相关系数,表示Kriging模型中相干函数的形式;表示参数空间中两个参数的举例。Rij和dij均为中间参数。i=1,2,…,9;j=1,2,…,9。
中间参数r=[ri],其中ri表示已知的9个插值点与被插值点(L',K')之间的相干函数,di表示已知的9个插值点与被插值点(L',K')之间的空间距离,i=1,2,…,9。
步骤107:求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
图3为本发明地震波初至检测方法实施例中检测过程的流程图。
参见图3,检测过程如下:
步骤301:获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
步骤302:获取检波器采集的待检测的地震记录序列。
待检测的地震记录序列为p0,p1,…,pN-1,其中N为待检测的地震记录序列的长度,即待检测的地震记录序列所包含的地震数据的数量。
步骤303:建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗。
步骤304:对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据。
其中,对所述地震记录序列的左侧和右侧进行数据填充,具体包括:将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
步骤305:从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达。
所述第一时间窗内的数据序列为其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;pn为位于所述第一时间窗中心的数据;P表示第一时间窗的数据序列,p表示地震记录序列内的各个数据;由于第二时间窗紧随第一时间窗,则所述第二时间窗内的数据序列为
其中,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达,具体包括:
获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;
求取所述第一时间窗Pn内所有数据的最大值和最小值,得到取值区间;其中最小值记为p(n) min=min[Pn],最大值记为p(n) max=max[Pn],则取值区间为[p(n) min,p(n) max]。
将所述取值区间划分为K段,得到K段小区间。其中K的取值一般为10;各个小区间的序号为k,划分之后的每个小区间可以表示为:
k=1,2,…,K
n=1,2,…,N
计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k)。
计算数量差值比率:
根据参数K和置信度α(α一般取为95%),通过查询卡方表获取卡方值χ2(α,K-1);
判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;
若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,即γ(n)≥χ2(α,K-1),则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,即γ(n)<χ2(α,K-1),则确定在所述地震记录序列的第n个数据处没有地震波动到达。
步骤306:记录所有有地震波动到达的数据在所述地震记录序列中的序号nFB,并计算每个有地震波动到达的数据对应的初至。各个有地震波动到达的数据所对应的地震波初至为t=(nFB-1)*T。其中T为检波器采样时间间隔,单位为秒。
地震波初至检测所面临的最大问题是噪声对检测精度的不良影响。噪声超过一定程度,会造成检测精度极速下降。而本发明则考察地震记录数据落入区间的次数,这种统计方法可以在一定程度上避免噪声的影响。
因此,本发明是对相邻的两个数据序列(P和Q)的分布进行统计刻画,如果两者差值之和小于卡方值,则认为这两个数据序列(P和Q)具有相同的分布,没有地震波动到达;如果大于卡方值,则认为这两个数据序列(P和Q)的分布不同,意味着有地震波动到达。利用分布是否相同,可以在一定程度上避免噪声的影响,提高检测精度,并且计算方法十分简单,原理非常明确。
本发明在训练过程中采用Kriging插值模型进行训练。该Kriging插值模型包括两个部分,第一部分利用一阶线性回归函数,可以较好地降低脉冲噪声(尤其是强脉冲噪声)对参数的影响,使得方法具有更好的鲁棒性。第二部分采用指数函数,指数函数强调数据之间的关系,且具有各向同性,与本发明所处理的功率数据正好契合。
综上,本发明的方法具有较高的鲁棒性、实用性和自适应性,可以有效降低噪声对初至拾取的影响。
图4为本发明地震波初至检测系统实施例的系统结构图。
参见图4,该地震波初至检测系统,包括:
时间窗最优长度获取模块401,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
序列获取模块402,用于获取检波器采集的待检测的地震记录序列。
时间窗建立模块403,用于建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗。
数据填充模块404,用于对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据。
所述数据填充模块404,具体包括:序列左填充单元和序列右填充单元。其中序列左填充单元,用于将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;序列右填充单元,用于将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
地震波初至检测模块405,用于从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达。
所述地震波初至检测模块405,具体包括:
子序列最优个数获取单元,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;区间计算单元,用于求取所述第一时间窗内所有数据的最大值和最小值,得到取值区间;区间划分单元,用于将所述取值区间划分为K段,得到K段小区间;区间数据量计算单元,用于计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k);其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;k为小区间序号,k=1,2,…,K;比率计算单元,用于计算数量差值比率:卡方值查询单元,用于根据参数K和置信度α,通过查询卡方表获取卡方值χ2(α,K-1);判断单元,用于判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;地震初至确定单元,用于若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,则确定在所述地震记录序列的第n个数据处没有地震波动到达。
记录与初至计算模块406,用于记录所有有地震波动到达的数据在所述地震记录序列中的序号,并计算每个有地震波动到达的数据对应的初至。
训练模块407,用于利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练;所述训练模块407包括:
已知序列获取单元,用于获取已确定初至的地震序列的实际初至值;坐标系建立单元,用于以子序列个数为横坐标,时间窗长度为纵坐标建立二维坐标系,根据时间窗长度的取值范围和子序列个数的取值范围确定所述二维坐标系内的取值区域;插值点初至预测计算单元,用于在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,计算每个插值点所对应的初至预测值;预测误差计算单元,用于计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;网格划分单元,用于将所述取值区域进行网格划分,获取各个网格点作为被插值点;线性拟合单元,用于利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值;最优插值点计算单元,用于求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
根据本发明提供的具体实施例,本发明公开了以下技术效果:本发明的地震波初至检测方法及系统,利用卡方检验的方法识别地震波初至,关注地震记录数据落入某个区间的数量,能够避免依赖于地震数据本身的大小,从而有效避免噪声的影响,提高提高检测准确度。
对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。
Claims (6)
1.一种地震波初至检测方法,其特征在于,包括:
获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
获取检波器采集的待检测的地震记录序列;
建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗;
对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据;
从所述第一时间窗的中心对准填充前的所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准填充前的所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达;
记录所有有地震波动到达的数据在所述地震记录序列中的序号,并计算每个有地震波动到达的数据对应的初至;
所述对所述地震记录序列的左侧和右侧进行数据填充,具体包括:
将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;
将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
2.根据权利要求1所述的一种地震波初至检测方法,其特征在于,所述每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达,具体包括:
获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;
求取所述第一时间窗内所有数据的最大值和最小值,得到取值区间;
将所述取值区间划分为K段,得到K段小区间;
计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k);其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;k为小区间序号,k=1,2,…,K;
计算数量差值比率:
根据参数K和置信度α,通过查询卡方表获取卡方值χ2(α,K-1);
判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;
若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,则确定在所述地震记录序列的第n个数据处没有地震波动到达。
3.根据权利要求2所述的一种地震波初至检测方法,其特征在于,利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练的过程为:
获取已确定初至的地震序列的实际初至值;
以子序列个数为横坐标,时间窗长度为纵坐标建立二维坐标系,根据时间窗长度的取值范围和子序列个数的取值范围确定所述二维坐标系内的取值区域;
在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,计算每个插值点所对应的初至预测值;
计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;
将所述取值区域进行网格划分,获取各个网格点作为被插值点;
利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值;
求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
4.一种地震波初至检测系统,其特征在于,包括:
时间窗最优长度获取模块,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的时间窗最优长度L;其中L为奇数;
序列获取模块,用于获取检波器采集的待检测的地震记录序列;
时间窗建立模块,用于建立第一时间窗和第二时间窗,所述第一时间窗和所述第二时间窗的长度均为时间窗最优长度L;所述第二时间窗紧随所述第一时间窗;
数据填充模块,用于对所述地震记录序列的左侧和右侧进行数据填充,使所述第一时间窗和所述第二时间窗从所述第一时间窗的中心对准所述地震记录序列的第一个数据滑动至所述第一时间窗的中心对准所述地震记录序列的最后一个数据时,所述第一时间窗和所述第二时间窗内均填充满L个数据;
地震波初至检测模块,用于从所述第一时间窗的中心对准填充前的所述地震记录序列的第一个数据滑动所述第一时间窗和所述第二时间窗,直至所述第一时间窗的中心对准填充前的所述地震记录序列的最后一个数据,每滑动一个数据的步长,对所述第一时间窗内的数据与所述第二时间窗内的数据进行卡方检验,从而确定所述地震记录序列中的各个数据处是否有地震波动到达;
记录与初至计算模块,用于记录所有有地震波动到达的数据在所述地震记录序列中的序号,并计算每个有地震波动到达的数据对应的初至;
所述数据填充模块,具体包括:
序列左填充单元,用于将所述地震记录序列的左侧补充至少个数据;左侧补充的至少个数据均与所述地震记录序列最左侧的数据相同;
序列右填充单元,用于将所述地震记录序列的右侧补充至少个数据;右侧补充的至少个数据均与所述地震记录序列最右侧的数据相同。
5.根据权利要求4所述的一种地震波初至检测系统,其特征在于,所述地震波初至检测模块,具体包括:
子序列最优个数获取单元,用于获取利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练得到的子序列最优个数K;
区间计算单元,用于求取所述第一时间窗内所有数据的最大值和最小值,得到取值区间;
区间划分单元,用于将所述取值区间划分为K段,得到K段小区间;
区间数据量计算单元,用于计算所述第一时间窗内的数据介于每个所述小区间内的数量F(n) P(k)和所述第二时间窗内的数据介于每个所述小区间内的数量F(n) Q(k);其中n表示位于所述第一时间窗中心的数据在所述地震记录序列中的序号;k为小区间序号,k=1,2,…,K;
比率计算单元,用于计算数量差值比率:
卡方值查询单元,用于根据参数K和置信度α,通过查询卡方表获取卡方值χ2(α,K-1);
判断单元,用于判断所述数量差值比率是否大于或等于所述卡方值,得到第一判断结果;
地震初至确定单元,用于若所述第一判断结果表示所述数量差值比率大于或等于所述卡方值,则确定在所述地震记录序列的第n个数据处有地震波动到达;若所述第一判断结果表示所述数量差值比率小于所述卡方值,则确定在所述地震记录序列的第n个数据处没有地震波动到达。
6.根据权利要求5所述的一种地震波初至检测系统,其特征在于,还包括训练模块,所述训练模块用于利用一阶线性回归函数和克里金插值方法对已确定初至的地震序列进行训练;所述训练模块包括:
已知序列获取单元,用于获取已确定初至的地震序列的实际初至值;
坐标系建立单元,用于以子序列个数为横坐标,时间窗长度为纵坐标建立二维坐标系,根据时间窗长度的取值范围和子序列个数的取值范围确定所述二维坐标系内的取值区域;
插值点初至预测计算单元,用于在所述取值区域内随机选取预设数量的点作为插值点,以每个插值点对应的纵坐标作为时间窗建立时第一时间窗和第二时间窗的长度,以每个插值点对应的横坐标作为取值区间的划分个数,计算每个插值点所对应的初至预测值;
预测误差计算单元,用于计算每个插值点所对应的初至预测值相对所述实际初至值的初至拾取误差,得到每个插值点的预测误差;
网格划分单元,用于将所述取值区域进行网格划分,获取各个网格点作为被插值点;
线性拟合单元,用于利用每个插值点的预测误差对每个被插值点进行线性拟合,得到拟合值;
最优插值点计算单元,用于求取所有拟合值中的最小值对应的插值点,得到最优插值点;所述最优插值点对应的横坐标即为子序列最优个数,所述最优插值点对应的纵坐标即为时间窗最优长度。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810768659.XA CN108535777B (zh) | 2018-07-13 | 2018-07-13 | 一种地震波初至检测方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810768659.XA CN108535777B (zh) | 2018-07-13 | 2018-07-13 | 一种地震波初至检测方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108535777A CN108535777A (zh) | 2018-09-14 |
CN108535777B true CN108535777B (zh) | 2019-07-05 |
Family
ID=63488247
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810768659.XA Expired - Fee Related CN108535777B (zh) | 2018-07-13 | 2018-07-13 | 一种地震波初至检测方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108535777B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110954940B (zh) * | 2018-09-26 | 2021-08-24 | 中国石油化工股份有限公司 | 一种基于地表一致性模型估计的初至质量控制方法 |
CN110361779B (zh) * | 2019-07-14 | 2021-10-01 | 广东石油化工学院 | 一种基于卡方分布的微震事件检测方法和系统 |
CN111077569B (zh) * | 2019-12-23 | 2022-05-06 | 中国石油天然气股份有限公司 | 全波形反演中分时窗提取数据的方法及装置 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104076392A (zh) * | 2014-05-28 | 2014-10-01 | 中国矿业大学(北京) | 基于网格搜索和牛顿迭代的微震震源定位联合反演方法 |
CN107526713A (zh) * | 2017-07-04 | 2017-12-29 | 北京航天易联科技发展有限公司 | 一种被动式太赫兹人体安检仪成像时间的确定方法和装置 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6686738B2 (en) * | 2001-04-17 | 2004-02-03 | Baker Hughes Incorporated | Method for determining decay characteristics of multi-component downhole decay data |
-
2018
- 2018-07-13 CN CN201810768659.XA patent/CN108535777B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104076392A (zh) * | 2014-05-28 | 2014-10-01 | 中国矿业大学(北京) | 基于网格搜索和牛顿迭代的微震震源定位联合反演方法 |
CN107526713A (zh) * | 2017-07-04 | 2017-12-29 | 北京航天易联科技发展有限公司 | 一种被动式太赫兹人体安检仪成像时间的确定方法和装置 |
Non-Patent Citations (4)
Title |
---|
A method for microseismic event detection and P-phase picking;G-Akis Tselentis,et al.;《SEG San Antonio 2011 Annual Meeting》;20111231;第1638-1642页 |
A multi-window algorithm for real-time automatic detection and picking of P-phases of microseismic events;Zuolin Chen and Robert R. Stewart;《CREWES Research Report》;20061231;第1-9页 |
Strategy for automated analysis of passive microseismic data based on S-transform, Otsu’s thresholding, and higher order statistics;G-Akis Tselentis,et al.;《GEOPHYSICS》;20121231;第43-54页 |
基于STA/LTA 方法的微地震事件自动识别技术;段建华,等;《煤田地质与勘探 》;20150228;第76-81页 |
Also Published As
Publication number | Publication date |
---|---|
CN108535777A (zh) | 2018-09-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108535773B (zh) | 一种微地震事件检测方法及系统 | |
CN108535777B (zh) | 一种地震波初至检测方法及系统 | |
EP2243043B1 (en) | Subsurface prediction method and system | |
Rentsch et al. | Fast location of seismicity: A migration-type approach with application to hydraulic-fracturing data | |
CN105425292B (zh) | 一种油气预测方法及装置 | |
CN105386756B (zh) | 一种应用应变量计算脆性地层孔隙度的方法 | |
RU2008147704A (ru) | Анализ повторных съемок по данным электромагнитной разведки | |
CN110389382B (zh) | 一种基于卷积神经网络的油气藏储层表征方法 | |
CN112147684B (zh) | 表征同沉积断层活动强度的方法和装置 | |
CN104597494A (zh) | 地震地层体分析方法及装置 | |
CN112180433B (zh) | 地震初至波拾取方法及装置 | |
CN105093318B (zh) | 一种自适应波动方程波场延拓静校正方法 | |
CN106353807B (zh) | 裂缝识别方法和装置 | |
AU2015200555B2 (en) | Correction of sea surface state | |
LIU et al. | Estimation of Near‐Surface Velocity for 3‐D Complicated Topography and Seismic Tomographic Static Corrections | |
GB2525072A (en) | Correction of sea surface state | |
CN108535778B (zh) | 一种地震波初至检测方法及系统 | |
CN112114358B (zh) | 一种基于三维地震资料表征的地下火山通道识别方法 | |
CN102778691B (zh) | 一种计算检波器组内静校正时差的方法 | |
EP3371627B1 (en) | Three-dimensional, stratigraphically-consistent seismic attributes | |
CN105891887B (zh) | 基于叠加数据的速度纵横向高密度分析方法 | |
Finney | The reality of GSSPs | |
CN112241020B (zh) | 稀疏地震数据采集中确定欠采样率的方法和装置 | |
CN110646846B (zh) | Vti介质各向异性参数确定方法、装置和设备 | |
CN104570100B (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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190705 Termination date: 20200713 |
|
CF01 | Termination of patent right due to non-payment of annual fee |