CN108535773B - 一种微地震事件检测方法及系统 - Google Patents
一种微地震事件检测方法及系统 Download PDFInfo
- Publication number
- CN108535773B CN108535773B CN201810768507.XA CN201810768507A CN108535773B CN 108535773 B CN108535773 B CN 108535773B CN 201810768507 A CN201810768507 A CN 201810768507A CN 108535773 B CN108535773 B CN 108535773B
- Authority
- CN
- China
- Prior art keywords
- data
- microseism
- time window
- sequence
- value
- 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 74
- 238000000034 method Methods 0.000 claims abstract description 45
- 238000005457 optimization Methods 0.000 claims abstract description 21
- 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 19
- 238000012549 training Methods 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 8
- 230000001502 supplementing effect Effects 0.000 claims description 5
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 235000013399 edible fruits Nutrition 0.000 claims 1
- 238000007689 inspection Methods 0.000 claims 1
- 238000009826 distribution Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000001419 dependent effect Effects 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
- 238000004458 analytical method Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 210000003739 neck Anatomy 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 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/288—Event detection in seismic signals, e.g. microseismics
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Emergency Management (AREA)
- Business, Economics & Management (AREA)
- Acoustics & Sound (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开一种微地震事件检测方法及系统。该方法包括:获取利用一阶线性回归函数和克里金插值方法对已标记微地震波动的微震检测序列进行训练得到的时间窗最优长度L;获取按时间顺序的待检测微地震数据序列;建立第一时间窗和紧随第一时间窗的第二时间窗,两个时间窗的长度均为L;对微地震数据序列的左侧和右侧进行数据填充;从第一时间窗的中心对准微地震数据序列的第一个数据滑动两个时间窗,直至第一时间窗的中心对准微地震数据序列的最后一个数据,每滑动一个数据的步长,对两个时间窗内的数据进行卡方检验确定微地震数据序列中的各个数据处是否发生微地震波动。本发明的方法及系统,能有效避免噪声对检测精度的影响,提高检测准确度。
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:记录所有发生微地震波动的数据在所述微地震数据序列中的序号。
微地震事件检测所面临的最大问题是噪声对检测精度的不良影响。噪声超过一定程度,会造成检测精度极速下降。而本发明则考察微地震数据落入区间的次数,这种统计方法可以在一定程度上避免噪声的影响。
因此,本发明是对相邻的两个数据序列(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 |
---|---|---|---|
CN201810768507.XA CN108535773B (zh) | 2018-07-13 | 2018-07-13 | 一种微地震事件检测方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810768507.XA CN108535773B (zh) | 2018-07-13 | 2018-07-13 | 一种微地震事件检测方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108535773A CN108535773A (zh) | 2018-09-14 |
CN108535773B true CN108535773B (zh) | 2019-07-05 |
Family
ID=63488246
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810768507.XA Expired - Fee Related CN108535773B (zh) | 2018-07-13 | 2018-07-13 | 一种微地震事件检测方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108535773B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110146919B (zh) * | 2019-06-25 | 2020-12-04 | 广东石油化工学院 | 基于正交投影的微震事件检测方法和系统 |
CN110333530B (zh) * | 2019-06-25 | 2020-11-10 | 广东石油化工学院 | 一种微震事件检测方法和系统 |
CN110146921B (zh) * | 2019-06-28 | 2021-06-11 | 广东石油化工学院 | 基于狄拉克分布概率的微震事件检测方法和系统 |
CN110414723B (zh) * | 2019-07-09 | 2022-05-17 | 中国石油大学(北京) | 基于微震事件的裂缝油气藏历史拟合的方法、装置及系统 |
CN110361779B (zh) * | 2019-07-14 | 2021-10-01 | 广东石油化工学院 | 一种基于卡方分布的微震事件检测方法和系统 |
CN112180430B (zh) * | 2020-09-23 | 2021-08-20 | 中国矿业大学 | 一种存在干扰信号下的矿震p波初至识别方法 |
CN112578442B (zh) * | 2020-11-30 | 2021-12-28 | 青岛海洋地质研究所 | 一种去除海洋地震勘探尾流噪音的方法 |
CN112965104B (zh) * | 2021-02-24 | 2023-03-28 | 中海石油(中国)有限公司 | 一种智能油气丛式井网井下微地震监控方法 |
CN115016001B (zh) * | 2022-08-09 | 2023-05-16 | 北京科技大学 | 一种基于维度变换的微震信号到时自动拾取方法 |
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 CN201810768507.XA patent/CN108535773B/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 |
---|---|
CN108535773A (zh) | 2018-09-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108535773B (zh) | 一种微地震事件检测方法及系统 | |
CN104297785B (zh) | 岩相约束储层物性参数反演方法及装置 | |
Strebelle et al. | Modeling of a deepwater turbidite reservoir conditional to seismic data using principal component analysis and multiple-point geostatistics | |
AU2011356685B2 (en) | Methods and systems regarding models of underground formations | |
CN104636980B (zh) | 针对河道砂油藏类型油气汇集条件的地球物理表征方法 | |
CN108897038B (zh) | 一种微地震事件检测方法及系统 | |
CN105425292A (zh) | 一种油气预测方法及装置 | |
WO2009076149A2 (en) | Systems and methods for utilizing cell based flow simulation results to calculate streamline trajectories | |
CN111399048B (zh) | 一种对断溶体计算相关属性及数据加权重构的方法 | |
CN109655903A (zh) | 页岩层横波速度预测方法及系统 | |
CN108535777B (zh) | 一种地震波初至检测方法及系统 | |
CN112147684A (zh) | 表征同沉积断层活动强度的方法和装置 | |
CN113552621B (zh) | 页岩气地应力确定方法和装置 | |
CN104316958A (zh) | 一种识别不同尺度地层断裂的相干处理方法 | |
CN104181598A (zh) | 计算地层不连续性属性值的方法及装置 | |
CN110389382A (zh) | 一种基于卷积神经网络的油气藏储层表征方法 | |
CN105043390A (zh) | 基于泛克里金法的重力场插值方法 | |
RU2289829C1 (ru) | Способ геофизической разведки для выявления нефтегазовых объектов | |
Deutsch et al. | Challenges in reservoir forecasting | |
Grimstad et al. | Identification of unknown permeability trends from history matching of production data | |
Ganguly et al. | Spatiotemporal variations in the connectivity of wetting and nonwetting phases during water alternating gas injection | |
CN108535778B (zh) | 一种地震波初至检测方法及系统 | |
CN105572732A (zh) | 一种属性变化率逐次晋级的裂缝发育带检测方法 | |
CN114153002A (zh) | 储层天然裂缝三维地质建模方法、装置、电子设备及介质 | |
CN114859414B (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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190705 Termination date: 20200713 |