CN110907991A - 基于数据场势值的震源定位方法、系统及可读存储介质 - Google Patents
基于数据场势值的震源定位方法、系统及可读存储介质 Download PDFInfo
- Publication number
- CN110907991A CN110907991A CN201911262949.8A CN201911262949A CN110907991A CN 110907991 A CN110907991 A CN 110907991A CN 201911262949 A CN201911262949 A CN 201911262949A CN 110907991 A CN110907991 A CN 110907991A
- Authority
- CN
- China
- Prior art keywords
- microseismic
- wave
- signal
- arrival
- data set
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 107
- 238000003860 storage Methods 0.000 title claims abstract description 12
- 238000005070 sampling Methods 0.000 claims abstract description 22
- 238000004590 computer program Methods 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 11
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 230000006870 function Effects 0.000 description 18
- 238000010586 diagram Methods 0.000 description 15
- 230000004807 localization Effects 0.000 description 7
- 230000000694 effects Effects 0.000 description 6
- 238000001615 p wave Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000010356 wave oscillation Effects 0.000 description 2
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000005290 field theory Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012502 risk assessment Methods 0.000 description 1
- 230000001960 triggered 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
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- 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/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
-
- 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/65—Source localisation, e.g. faults, hypocenters or reservoirs
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Environmental & Geological Engineering (AREA)
- Acoustics & Sound (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Business, Economics & Management (AREA)
- Emergency Management (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于数据场势值的震源定位方法、系统及可读存储介质,该方法首先进行微震信号去噪;再用自动和人工相结合拾取P波初至;采用Bootstrap方法采样,得到多个P波初至子数据集;进而采用双差法和数据场势值获取微震事件初定位结果;利用初定位结果计算震源到传感器的理论传播时间,去除P波初至大误差拾取信号;以去除大误差拾取的P波初至作为新的数据,开展Bootstrap采样和数据场势值定位得到最终的震源位置。该方法可去除微震信号高频和工频噪音,提高了信噪比,使得自动方法拾取的P波初至更准确;基于数据场势值得到定位密度较大的点,能够降低初始迭代中心的影响,从而使微震定位结果更为稳定。
Description
技术领域
本发明属于微震监测领域,尤其是涉及一种基于数据场势值的震源定位方法、系统及可读存储介质。
背景技术
微震监测技术作为一种在岩土工程,尤其是采矿工程和隧道工程中广泛应用的地球物理技术,其通过布设传感器接收断层滑移产生的微震信号,进而分析微震发生位置、震级和能量变化等特征,从而推断岩体的力学状态,并采取有效的防控措施。其中,微震定位可以确定岩体的破裂位置,是后续震级计算、震源机制反演和风险分析的基础,因此高精度震源定位至关重要。
目前微震信号降噪主要基于低通滤波、小波阀值去噪和经验模态分解(EMD)去噪。低通滤波只能滤除特定频率的噪音;小波分析如何选取小波基和分解尺度是应用中的难题,且其将信号按尺度固定频段分解,每一频段只与信号的采样频率有关,而与信号本身信息无关,很容易受到信号中相邻谐波成分的交叠影响。EMD存在频率混淆、端点效应明显等问题,降低了信号的保真度。
微震定位主要借助传感器间P波理论到时差与监测到时差或P波理论传播时间与传感器监测传播时间来确定震源位置,其都是基于P波初至时间或P波相对到时差。国内外提出了长短时窗平均值比法(STA/LTA法)、池赤准则法(AIC法)和高阶统计量法(PAI-S/K法)、分形维数法和神经网络法等用于P波初至拾取,上述方法均以P波起振点作为P波初至时刻,然而其可能受传播而衰减掉,即拾取到的P波初至时刻并非真实的P波初至时刻。有鉴于此,一些研究采用互相关法拾取P波相对到时,而非起振时刻,其在天然地震领域取得了较好了效果。然而,矿山微震信号传播距离较近,尾波非常复杂,难以使用互相关法拾取P波相对到时差。此外,由于噪音、相邻事件信号错误划分等原因,P波初至本身可能存在大误差。
微震定位大多采用到时差的L1和L2范数来定义目标函数,但传统的L1和L2范数定位易受大拾取误差影响,导致定位结果不太稳定。为解决上述问题,许多学者提出了很多新的微震定位方法。董陇军等(2019)就深部复杂采场中波速传播、信号采集等问题做了研究,提出了一种能够提高微震事件定位精度的方法,该方法是基于理论分析和迭代求解的协同定位方法。李夕兵等(2016)提出了一种虚拟场优化方法,该方法利用三维空间中站台确定的双曲面公共交点来定位,其采用了指数衰减函数来降低大拾取误差引起的震源定位不稳定的问题。
发明内容
本发明提出了一种基于数据场势值的震源定位方法、系统及可读存储介质,该方法针对现有P波初至拾取方法不稳定,提出自动与人工相结合拾取P波初至时刻,可有效提高拾取精度且保证较快的数据处理速度。从P波到时数据中有放回抽样展开定位,降低大误差拾取的影响。基于数据场势值得到定位密度较大的点,够降低初始迭代中心的影响,从而使微震定位结果更为稳定。此外,初定位后可去除大误差拾取信号,进而开展高精度定位。
一种基于数据场势值的微震信号大误差拾取去除及震源定位方法,包括以下步骤:
步骤1:对微震信号的高频噪音和工频噪音进行去噪;
步骤2:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
步骤3:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
步骤4:双差法和数据场势值实现微震事件初定位;
采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
步骤5:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
步骤6:以更新后的微震数据集,重复步骤3-4实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
进一步地,对微震信号的高频噪音和工频噪音进行去噪的过程如下:
首先,采用集合经验模态分解将微震信号x(m)自适应分解为多个本征模态分量IMFi和1个剩余分量rk,i=1,2,…,k,k为经验模态分量的个数,对前两个本征模态分量IMF采用软阀值去除高频噪音,将去噪后IMF1和IMF2与其他IMF分量IMF3,IMF4,…,IMFk和剩余分量rk相加,得到去除高频噪音的微震信号;
其次,采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪,得到去除高频和工频后的微震信号;
x(m)为微震信号,m=1,2,…,M,M为微震信号的总采样点数。
进一步地,对高信噪比的去噪后的微震信号,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1,再采用池赤准则法AIC在[T1-10ms,T1+10ms]范围内确定P波初至的准确位置ti;
其中,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1时,使用P波初至后、前0.1s内微震信号振幅均值比,且触发阀值取3。
进一步地,采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集的具体过程如下:
使用Bootstrap采样法对微震数据集S1,S2,…,Sn有放回采样n次,得到一个新的数据集D1,D2,…,Dn;从新的数据集D1,D2,…,Dn中选出所有不相同的数据集,记为{U1,U2,…,Un1},n1≤n;重复上述步骤N次,得到N个数据集{U1,U2,…,Un1},6≤n1,即N个P波初至子数据集,N的取值范围为1000~3000。
进一步地,每个P波初至子数据集对应的定位位置采用以下公式进行计算:
其中,x0,y0,z0,t0分别对应微震事件的三维位置和发生时刻,li为各台站到震源之间的距离,ti为P波到达各台站的时刻,i=1,2,3…,n1,vp为P波传播速度。
上述公式采用Nelder-Mead法快速求解;
每个P波初至子数据集在对应定位位置的数据场势值按照以下公式计算:
其中,(xi,yi,zi)和(xj,yj,zj)分别为第i和j个P波初至子数据集得到的定位位置,σ为距离影响因子,10≤σ≤200。
进一步地,根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集的具体过程如下:
计算传感器位置(xi,yi,zi)与微震事件初定位事件位置(x′,y′,z′)之间P波理论传播时间Ti=li/vp,再计算Ti与基于观测到时的传播时间(ti-t0)的差值,当差值小于等于T时,0.01s≤T≤0.03s,认为该P波初至拾取是可靠的,进入下一步;当差值大于T时,认为该拾取具有较大的误差,删除该P波初至拾取对应的数据集Si。
进一步地,对前两个本征模态分量IMF采用软阀值去除高频噪音时,所采用的软阀值的计算方法包括无偏风险估计阀值、固定阀值、启发式阀值和极大极小阀值中的一种;
采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪时,正弦函数拟合时首先截取一段P波初至前的微震信号(200us~500us),从而得到截取微震信号的局部最大值和最小值点,再根据两个相邻的局部最大值点和局部最小值点确定正弦函数的周期,并依据首个最大值点确定相位,进而将正弦函数延拓至长度等于微震信号长度M的时间序列。
进一步地,微震信号的信噪比计算方式为SNR(dB)=20lg(Asignal/Anoise),Asignal和Anoise分别为信号和噪音的平均振幅;
高信噪比取SNR≥6,低信噪比取SNR<6。
一种基于数据场势值的微震信号大误差拾取去除及震源定位系统,包括:
微震信号去噪单元:对微震信号的高频噪音和工频噪音进行去噪;
P波初至时刻拾取单元:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
P波初至子数据集提取单元:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
微震事件初定位单元:采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
微震数据集更新单元:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
微震事件最终定位单元:以更新后的微震数据集,重复利用P波初至子数据集提取单元、微震事件初定位单元以及微震数据集更新单元,实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
一种可读存储介质,包括计算机程序指令,其特征在于:所述计算机程序指令被处理终端执行时使所述处理终端执行一种基于数据场势值的微震信号大误差拾取去除及震源定位方法。
有益效果
本发明提供一种基于数据场势值的震源定位方法、系统及可读存储介质,主要用于减小微震事件定位中P波大拾取误差的影响,解决传统定位方法受初始迭代中心的影响较大以及微震事件定位不稳定等问题。本方法包括如下步骤:首先采用集合经验模态分解(EEMD)软阈值与正弦函数相结合去除微震信号高频和工频噪音;再采用信噪比、长短时窗均值比(STA/LTA)法、池赤准则(AIC)法实现高信噪比微震信号P波初至的自动拾取,并结合人工拾取低信噪比信号;采用Bootstrap方法对P波初至数据集进行有放回采样,得到多个P波初至子数据集;进而采用双差法对每一个P波初至子数据集展开震源定位,并取数据场势值最大的30~60个点位置均值作为微震事件初定位结果;利用初定位结果计算震源到传感器的理论传播时间,再与监测到的传播时间比较,去除P波初至大误差拾取信号;以去除大误差拾取的P波初至作为新的数据,开展Bootstrap采样和数据场势值定位得到最终的震源位置。该方法可去除微震信号高频和工频噪音,提高信噪比,使得自动方法拾取的P波初至更多、更准确,降低了P波初至的整体拾取时间;从P波到时数据中有放回抽样展开定位,可降低大误差拾取的影响;基于数据场势值得到定位密度较大的点,能够降低初始迭代中心的影响,从而使微震定位结果更为稳定;此外,初定位后可去除大误差拾取信号,进而开展高精度定位。
附图说明
图1是本发明所述方法流程图;
图2是微震信号去噪及自动拾取效果图,其中,(a)为原始微震信号,(b)为去噪微震信号及其AIC序列拾取;
图3是生成理论数据的传感器与待测试事件位置图;
图4是数据场势值初定位和重定位案例图,其中,(a)为数据场初定位案例,(b)为去噪后重定位案例;
图5是三种不同方法的定位误差箱型图,其中,(a)为事件1双差法定位误差箱型图,(b)为事件2双差法定位误差箱型图,(c)为事件1数据场定位误差箱型图,(d)为事件2数据场定位误差箱型图,(e)为事件1去噪后重定位误差箱型图,(f)为事件2去噪后重定位误差箱型图。
具体实施方式
下面将结合附图1~5,对本发明做进一步地说明。
本发明所述的方法思想如下:针对常用的L1和L2范数方法定位结果不稳定、可能受大拾取误差影响等问题,提出一种基于数据场势值的微震信号大误差拾取去除及震源定位方法。该方法借助Bootstrap采样方法和数据场理论,得到定位密度较大的点,降低初始迭代中心的影响,得到较好的初定位结果。基于初定位利用理论传播时间与监测传播时间差去除大误差拾取信号,进而重复上述步骤实现高精度震源定位。
如图1所示,一种基于数据场的微震大误差拾取去除及震源定位方法,包括以下步骤:
步骤1:对微震信号的高频噪音和工频噪音进行去噪;
对微震信号的高频噪音和工频噪音进行去噪的过程如下:
首先,采用集合经验模态分解将微震信号x(m)自适应分解为多个本征模态分量IMFi和1个剩余分量rk,i=1,2,…,k,k为经验模态分量的个数,对前两个本征模态分量IMF采用软阀值去除高频噪音,将去噪后IMF1和IMF2与其他IMF分量IMF3,IMF4,…,IMFk和剩余分量rk相加,得到去除高频噪音的微震信号;
其次,采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪,得到去除高频和工频后的微震信号;
x(m)为微震信号,m=1,2,…,M,M为微震信号的总采样点数。
对前两个本征模态分量IMF采用软阀值去除高频噪音时,所采用的软阀值的计算方法包括无偏风险估计阀值、固定阀值、启发式阀值和极大极小阀值中的一种;
采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪时,正弦函数拟合时首先截取一段P波初至前的微震信号(200us~500us),从而得到截取微震信号的局部最大值和最小值点,再根据两个相邻的局部最大值点和局部最小值点确定正弦函数的周期,并依据首个最大值点确定相位,进而将正弦函数延拓至长度等于微震信号长度M的时间序列。
步骤2:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
对高信噪比的去噪后的微震信号,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1,再采用池赤准则法AIC在[T1-10ms,T1+10ms]范围内确定P波初至的准确位置ti;
其中,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1时,使用P波初至后、前0.1s内微震信号振幅均值比,且触发阀值取3。
微震信号的信噪比计算方式为SNR(dB)=20lg(Asignal/Anoise),Asignal和Anoise分别为信号和噪音的平均振幅;
高信噪比取SNR≥6,低信噪比取SNR<6。
步骤3:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集的具体过程如下:
使用Bootstrap采样法对微震数据集S1,S2,…,Sn有放回采样n次,得到一个新的数据集D1,D2,…,Dn;从新的数据集D1,D2,…,Dn中选出所有不相同的数据集,记为{U1,U2,…,Un1},n1≤n;重复上述步骤N次,得到N个数据集{U1,U2,…,Un1},6≤n1,即N个P波初至子数据集,N的取值范围为1000~3000。
步骤4:双差法和数据场势值实现微震事件初定位;
采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
每个P波初至子数据集对应的定位位置采用以下公式进行计算:
其中,x0,y0,z0,t0分别对应微震事件的三维位置和发生时刻,li为各台站到震源之间的距离,ti为P波到达各台站的时刻,i=1,2,3…,n1,vp为P波传播速度。
上述公式采用Nelder-Mead法快速求解;
每个P波初至子数据集在对应定位位置的数据场势值按照以下公式计算:
其中,(xi,yi,zi)和(xj,yj,zj)分别为第i和j个P波初至子数据集得到的定位位置,σ为距离影响因子,10≤σ≤200。
步骤5:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集的具体过程如下:
计算传感器位置(xi,yi,zi)与微震事件初定位事件位置(x′,y′,z′)之间P波理论传播时间Ti=li/vp,再计算Ti与基于观测到时的传播时间(ti-t0)的差值,当差值小于等于T时,0.01s≤T≤0.03s,认为该P波初至拾取是可靠的,进入下一步;当差值大于T时,认为该拾取具有较大的误差,删除该P波初至拾取对应的数据集Si。
步骤6:以更新后的微震数据集,重复步骤3-4实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
实施例
图2是微震信号去噪及自动拾取效果图。由图2(a)知,原始微震信号存在少量的高频信号和较大的工频信号,使得P波拾取较为困难。采用EEMD软阀值+正弦函数去噪的微震信号信噪比明显提升,P波初至更易于拾取(图2b)。图2b中垂直虚线为STA/LTA法拾取结果,垂直点线为AIC法局部最小值可取的范围。与人工拾取相比,AIC法能取得较好的拾取效果。
图3是生成理论数据的传感器与待测试事件位置图。微震传感器分布较为宽广,其中930m和1080m高度附近分别布置了12个传感器,1120m高度附近布置了4个传感器,传感器的采样频率均为6000Hz。为验证方法的有效性,选取两个典型事件测试定位方法:待测试事件1位于传感器阵列内部,待测试事件2位于传感器阵列外部。考虑到现实微震中,传感器通常不能全部触发,我们从这28个传感器中选取15个展开研究。由P波速度、震源位置和传感器位置生成理论传播时间,再对P波速度加上5%的高斯噪音,对每个传感器接收的理论传播时间加上2ms的高斯噪音。此外,我们还对上述噪音基础上,其中一个传感器加上大误差拾取展开了定位研究。
图4是数据场势值初定位和重定位案例图。选取事件2作为测试,并在高斯噪音的基础上,对其中一个传感器P波初至加上大误差噪音,传统定位2000次的结果如图4(a)中灰色小圆点,灰色越深,数据场势值越大,显然传统定位方法结果很不稳定。基于数据场势值的震源初定位可得到较稳定的定位结果,然而其定位误差仍较大。基于震源初定位去除大误差拾取信号的震源重定位如图4(b)所示,可知去噪后重定位更靠近事件真实位置。
图5是三种不同方法的定位误差箱型图。其中,Ex、Ey、Ez和ET分别为X、Y、Z方向绝对误差和总误差。由图可知,事件位于传感器阵列内部的定位误差整体上小于事件位于传感器阵列外部的定位误差,基于数据场势值的初定位结果明显优于传统的定位结果。当P波初至存在大误差拾取时,去噪后重定位结果优于初定位结果,表明大误差去噪的必要性和该专利定位的有效性。
一种基于数据场势值的微震信号大误差拾取去除及震源定位系统,包括:
微震信号去噪单元:对微震信号的高频噪音和工频噪音进行去噪;
P波初至时刻拾取单元:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
P波初至子数据集提取单元:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
微震事件初定位单元:采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
微震数据集更新单元:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
微震事件最终定位单元:以更新后的微震数据集,重复利用P波初至子数据集提取单元、微震事件初定位单元以及微震数据集更新单元,实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
应当理解,本发明各个实施例中的功能单元模块可以集中在一个处理单元中,也可以是各个单元模块单独物理存在,也可以是两个或两个以上的单元模块集成在一个单元模块中,可以采用硬件或软件的形式来实现。
一种可读存储介质,包括计算机程序指令,其特征在于:所述计算机程序指令被处理终端执行时使所述处理终端执行一种基于数据场势值的微震信号大误差拾取去除及震源定位方法。
本领域内的技术人员应明白,本申请的实施例可提供为方法、系统、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是:以上实施例仅用以说明本发明的技术方案而非对其限制,尽管参照上述实施例对本发明进行了详细的说明,所属领域的普通技术人员应当理解:依然可以对本发明的具体实施方式进行修改或者等同替换,而未脱离本发明精神和范围的任何修改或者等同替换,其均应涵盖在本发明的权利要求保护范围之内。
Claims (10)
1.一种基于数据场势值的微震信号大误差拾取去除及震源定位方法,其特征在于,包括以下步骤:
步骤1:对微震信号的高频噪音和工频噪音进行去噪;
步骤2:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
步骤3:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
步骤4:双差法和数据场势值实现微震事件初定位;
采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
步骤5:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
步骤6:以更新后的微震数据集,重复步骤3-4实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
2.根据权利要求1所述的方法,其特征在于,对微震信号的高频噪音和工频噪音进行去噪的过程如下:
首先,采用集合经验模态分解将微震信号x(m)自适应分解为多个本征模态分量IMFi和1个剩余分量rk,i=1,2,…,k,k为经验模态分量的个数,对前两个本征模态分量IMF采用软阀值去除高频噪音,将去噪后IMF1和IMF2与其他IMF分量IMF3,IMF4,…,IMFk和剩余分量rk相加,得到去除高频噪音的微震信号;
其次,采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪,得到去除高频和工频后的微震信号;
x(m)为微震信号,m=1,2,…,M,M为微震信号的总采样点数。
3.根据权利要求1所述的方法,其特征在于,对高信噪比的去噪后的微震信号,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1,再采用池赤准则法AIC在[T1-10ms,T1+10ms]范围内确定P波初至的准确位置ti;
其中,采用长短时窗均值比法拾取高信噪比信号的P波大致初至T1时,使用P波初至后、前0.1s内微震信号振幅均值比,且触发阀值取3。
4.根据权利要求1所述的方法,其特征在于,采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集的具体过程如下:
使用Bootstrap采样法对微震数据集S1,S2,…,Sn有放回采样n次,得到一个新的数据集D1,D2,…,Dn;从新的数据集D1,D2,…,Dn中选出所有不相同的数据集,记为{U1,U2,…,Un1},n1≤n;重复上述步骤N次,得到N个数据集{U1,U2,…,Un1},6≤n1,即N个P波初至子数据集,N的取值范围为1000~3000。
6.根据权利要求1所述的方法,其特征在于,根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集的具体过程如下:
计算传感器位置(xi,yi,zi)与微震事件初定位事件位置(x′,y′,z′)之间P波理论传播时间Ti=li/vp,再计算Ti与基于观测到时的传播时间(ti-t0)的差值,当差值小于等于T时,0.01s≤T≤0.03s,认为该P波初至拾取是可靠的,进入下一步;当差值大于T时,认为该拾取具有较大的误差,删除该P波初至拾取对应的数据集Si。
7.根据权利要求2所述的方法,其特征在于,对前两个本征模态分量IMF采用软阀值去除高频噪音时,所采用的软阀值的计算方法包括无偏风险估计阀值、固定阀值、启发式阀值和极大极小阀值中的一种;
采用正弦函数拟合的方式对去除高频噪音的微震信号的工频噪音降噪时,正弦函数拟合时首先截取一段P波初至前的微震信号(200us~500us),从而得到截取微震信号的局部最大值和最小值点,再根据两个相邻的局部最大值点和局部最小值点确定正弦函数的周期,并依据首个最大值点确定相位,进而将正弦函数延拓至长度等于微震信号长度M的时间序列。
8.根据权利要求1所述的方法,其特征在于,微震信号的信噪比计算方式为SNR(dB)=20lg(Asignal/Anoise),Asignal和Anoise分别为信号和噪音的平均振幅;
高信噪比取SNR≥6,低信噪比取SNR<6。
9.一种基于数据场势值的微震信号大误差拾取去除及震源定位系统,其特征在于,包括:
微震信号去噪单元:对微震信号的高频噪音和工频噪音进行去噪;
P波初至时刻拾取单元:对高信噪比和低信噪比的去噪后的微震信号,分别采用自动和人工拾取P波初至时刻;
P波初至子数据集提取单元:采用Bootstrap法对微震数据集S1,S2,…,Sn进行采样,得到多个P波初至子数据集;
Si=(xi,yi,zi,ti),i=1,2,…,n,n为有拾取微震信号的数目,(xi,yi,zi)为传感器的三维坐标位置,ti为步骤2拾取的P波初至时刻;
微震事件初定位单元:采用双差法对每一个P波初至子数据集展开震源定位,得到每个P波初至子数据集对应的定位位置;
计算每个P波初至子数据集在对应定位位置的数据场势值,取数据场势值较大的前K个数据场势值对应的定位位置均值作为微震事件初定位位置,K的个数取值为30~60个;
微震数据集更新单元:根据P波理论传播时间与基于观测到时的传播时间之间的差值,去除P波初至大误差拾取对应的微震数据集,得到更新后的微震数据集;
微震事件最终定位单元:以更新后的微震数据集,重复利用P波初至子数据集提取单元、微震事件初定位单元以及微震数据集更新单元,实现微震事件重定位,得到重定位的事件位置(x,y,z),并作为最终定位结果。
10.一种可读存储介质,包括计算机程序指令,其特征在于:所述计算机程序指令被处理终端执行时使所述处理终端执行权利要求1-8任一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911262949.8A CN110907991B (zh) | 2019-12-11 | 2019-12-11 | 基于数据场势值的震源定位方法、系统及可读存储介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911262949.8A CN110907991B (zh) | 2019-12-11 | 2019-12-11 | 基于数据场势值的震源定位方法、系统及可读存储介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110907991A true CN110907991A (zh) | 2020-03-24 |
CN110907991B CN110907991B (zh) | 2021-03-16 |
Family
ID=69824247
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911262949.8A Active CN110907991B (zh) | 2019-12-11 | 2019-12-11 | 基于数据场势值的震源定位方法、系统及可读存储介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110907991B (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112083480A (zh) * | 2020-08-24 | 2020-12-15 | 中国石油天然气集团有限公司 | 同步激发采集地震数据的实时监控方法及装置 |
CN112232321A (zh) * | 2020-12-14 | 2021-01-15 | 西南交通大学 | 一种振动数据干扰降噪方法、装置、设备及可读存储介质 |
CN112230270A (zh) * | 2020-12-14 | 2021-01-15 | 西南交通大学 | 一种地震预警方法、装置、设备及可读存储介质 |
CN112526602A (zh) * | 2020-11-16 | 2021-03-19 | 重庆大学 | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 |
CN112748465A (zh) * | 2020-12-30 | 2021-05-04 | 中国矿业大学(北京) | 基于岩石特征的震源机制反演方法及装置 |
CN113219068A (zh) * | 2021-05-14 | 2021-08-06 | 重庆大学 | 一种基于分组传感器解析解的圆筒声发射定位方法、系统、终端及可读存储介质 |
Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6965849B1 (en) * | 2000-02-10 | 2005-11-15 | Schlumberger Technology Corporation | Method of designing geophysical surveys |
US20140104980A1 (en) * | 2012-10-16 | 2014-04-17 | Colorado School Of Mines | Virtual electrode current injection using seismic focusing and seismoelectric conversion |
CN104133246A (zh) * | 2014-07-31 | 2014-11-05 | 中国石油天然气集团公司 | 一种微地震事件扫描定位方法及装置 |
CN105487114A (zh) * | 2015-12-08 | 2016-04-13 | 中南大学 | 一种微震信号p波初至点综合拾取方法 |
CN106094021A (zh) * | 2016-06-01 | 2016-11-09 | 北京科技大学 | 一种基于到时差数据库的微地震震源快速定位方法 |
CN106990435A (zh) * | 2017-06-07 | 2017-07-28 | 中煤科工集团西安研究院有限公司 | 一种减弱依赖初至拾取精度的微震定位方法及装置 |
CN107045141A (zh) * | 2017-02-24 | 2017-08-15 | 北京科技大学 | 基于反时到时差数据库的微地震/地震震源快速定位方法 |
CN107450098A (zh) * | 2017-08-03 | 2017-12-08 | 中煤科工集团西安研究院有限公司 | 一种煤层底板隐伏突水陷落柱动态定位方法 |
US9880305B2 (en) * | 2013-10-08 | 2018-01-30 | Altan Turgut | Method of passive acoustic depth determination in shallow water |
CN109655918A (zh) * | 2017-10-11 | 2019-04-19 | 中国石油化工股份有限公司 | 地面浅井微地震监测观测台站位置确定方法及系统 |
-
2019
- 2019-12-11 CN CN201911262949.8A patent/CN110907991B/zh active Active
Patent Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6965849B1 (en) * | 2000-02-10 | 2005-11-15 | Schlumberger Technology Corporation | Method of designing geophysical surveys |
US20140104980A1 (en) * | 2012-10-16 | 2014-04-17 | Colorado School Of Mines | Virtual electrode current injection using seismic focusing and seismoelectric conversion |
US9880305B2 (en) * | 2013-10-08 | 2018-01-30 | Altan Turgut | Method of passive acoustic depth determination in shallow water |
CN104133246A (zh) * | 2014-07-31 | 2014-11-05 | 中国石油天然气集团公司 | 一种微地震事件扫描定位方法及装置 |
CN105487114A (zh) * | 2015-12-08 | 2016-04-13 | 中南大学 | 一种微震信号p波初至点综合拾取方法 |
CN106094021A (zh) * | 2016-06-01 | 2016-11-09 | 北京科技大学 | 一种基于到时差数据库的微地震震源快速定位方法 |
CN107045141A (zh) * | 2017-02-24 | 2017-08-15 | 北京科技大学 | 基于反时到时差数据库的微地震/地震震源快速定位方法 |
CN106990435A (zh) * | 2017-06-07 | 2017-07-28 | 中煤科工集团西安研究院有限公司 | 一种减弱依赖初至拾取精度的微震定位方法及装置 |
CN107450098A (zh) * | 2017-08-03 | 2017-12-08 | 中煤科工集团西安研究院有限公司 | 一种煤层底板隐伏突水陷落柱动态定位方法 |
CN109655918A (zh) * | 2017-10-11 | 2019-04-19 | 中国石油化工股份有限公司 | 地面浅井微地震监测观测台站位置确定方法及系统 |
Non-Patent Citations (1)
Title |
---|
史鹏程 等: ""利用到时差分布图方法进行近场地表震源定位"", 《地球物理学报》 * |
Cited By (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112083480A (zh) * | 2020-08-24 | 2020-12-15 | 中国石油天然气集团有限公司 | 同步激发采集地震数据的实时监控方法及装置 |
CN112526602A (zh) * | 2020-11-16 | 2021-03-19 | 重庆大学 | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 |
CN112526602B (zh) * | 2020-11-16 | 2023-10-20 | 重庆大学 | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 |
CN112232321A (zh) * | 2020-12-14 | 2021-01-15 | 西南交通大学 | 一种振动数据干扰降噪方法、装置、设备及可读存储介质 |
CN112230270A (zh) * | 2020-12-14 | 2021-01-15 | 西南交通大学 | 一种地震预警方法、装置、设备及可读存储介质 |
CN112230270B (zh) * | 2020-12-14 | 2021-03-09 | 西南交通大学 | 一种地震预警方法、装置、设备及可读存储介质 |
CN112232321B (zh) * | 2020-12-14 | 2021-03-19 | 西南交通大学 | 一种振动数据干扰降噪方法、装置、设备及可读存储介质 |
CN112748465A (zh) * | 2020-12-30 | 2021-05-04 | 中国矿业大学(北京) | 基于岩石特征的震源机制反演方法及装置 |
CN112748465B (zh) * | 2020-12-30 | 2021-12-10 | 中国矿业大学(北京) | 基于岩石特征的震源机制反演方法及装置 |
CN113219068A (zh) * | 2021-05-14 | 2021-08-06 | 重庆大学 | 一种基于分组传感器解析解的圆筒声发射定位方法、系统、终端及可读存储介质 |
CN113219068B (zh) * | 2021-05-14 | 2022-03-22 | 重庆大学 | 一种基于分组传感器解析解的圆筒声发射定位方法、系统、终端及可读存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN110907991B (zh) | 2021-03-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110907991B (zh) | 基于数据场势值的震源定位方法、系统及可读存储介质 | |
CN106353792B (zh) | 一种适用于水力压裂微震震源定位的方法 | |
CN103837891B (zh) | 微地震初至的高精度拾取方法 | |
CN113189644B (zh) | 一种微震震源定位方法及系统 | |
Dong et al. | Some developments and new insights for microseismic/acoustic emission source localization | |
CN107479094B (zh) | 一种实现地震预警的方法 | |
CN112014883B (zh) | 一种基于Log-Cosh函数的微震震源定位方法、系统、装置及可读存储介质 | |
CN113238280B (zh) | 一种基于格林函数的地震监测方法 | |
CN106646598A (zh) | 一种fast‑aic法微地震信号拾取方法 | |
CN109375253A (zh) | 基于全部发震构造最大可信地震的地震动参数评价方法 | |
CN103913768A (zh) | 基于地震波资料对地表中浅层进行建模的方法及装置 | |
CN109375252A (zh) | 考虑不同发震构造最大可信地震的地震动参数评价方法 | |
CN112183407B (zh) | 一种基于时频域谱减法的隧道地震波数据去噪方法及系统 | |
CN114839673A (zh) | 多震源高效采集波场分离方法、分离系统及计算机设备 | |
CN108802814B (zh) | 一种隧道围岩微震波速的获取方法 | |
CN112180429B (zh) | 利用隧道爆破震动反演的不良地质构造探测系统及方法 | |
CN110703313B (zh) | 考虑传感器感度的声发射事件震级获取方法、系统及可读存储介质 | |
CN106646610B (zh) | 一种利用偏振约束aic算法自动拾取微震初至的算法 | |
CN112526602A (zh) | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 | |
CN111736208B (zh) | 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质 | |
CN109884691B (zh) | 用于随采地震信号的强单频和随机噪声压制方法及系统 | |
CN112731534A (zh) | 一种考虑p波初至系统误差的双声发射事件联合定位方法、系统、电子设备及可读介质 | |
CN111208568B (zh) | 一种时间域多尺度全波形反演方法及系统 | |
Sun et al. | The adaptive particle swarm optimization technique for solving microseismic source location parameters | |
CN113219526A (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 |