CN111722280A - 一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质 - Google Patents
一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质 Download PDFInfo
- Publication number
- CN111722280A CN111722280A CN202010601258.2A CN202010601258A CN111722280A CN 111722280 A CN111722280 A CN 111722280A CN 202010601258 A CN202010601258 A CN 202010601258A CN 111722280 A CN111722280 A CN 111722280A
- Authority
- CN
- China
- Prior art keywords
- wave
- arrival
- acoustic emission
- time
- bayes
- 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 80
- 238000005070 sampling Methods 0.000 claims abstract description 21
- 238000012614 Monte-Carlo sampling Methods 0.000 claims abstract 2
- 230000004807 localization Effects 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 11
- 238000001615 p wave Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 238000003860 storage Methods 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000006641 stabilisation Effects 0.000 claims description 3
- 238000011105 stabilization Methods 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 230000006870 function Effects 0.000 description 35
- 238000010586 diagram Methods 0.000 description 23
- 230000000694 effects Effects 0.000 description 14
- 238000004364 calculation method Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 238000011156 evaluation Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000000875 corresponding effect Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 240000008574 Capsicum frutescens Species 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000005669 field effect Effects 0.000 description 1
- 238000011478 gradient descent method Methods 0.000 description 1
- 230000035876 healing Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000006467 substitution reaction Methods 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/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- 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)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Data Mining & Analysis (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Computational Biology (AREA)
- Algebra (AREA)
- Evolutionary Biology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Operations Research (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明公开了一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质，该方法：首先采用低通滤波器对声发射信号高频噪音降噪；拾取确定声发射信号P波初至；建立基于P波初至理论数据与观测数据差的Bayes定位后验概率密度函数；进而采用马尔科夫链蒙特卡罗采样方法迭代求解P波初至系统观测误差；利用P波初至拾取到时tij加上P波初至系统观测误差Tj，得到校正后的P波初至数据；进而利用Bayes定位后验概率密度函数、MCMC采样方法对校正后的P波初至数据开展定位。该方法拾取P波初至数据速度较快、精度高，可去除P波初至系统观测误差，提高了P波初至数据的准确性，Bayes定位方法易于得到全局最优解。
Description
技术领域
本发明属于声发射领域，尤其是涉及一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质。
背景技术
P波初至到时拾取质量直接关系到声发射事件定位精度。目前常用绝对到时法确定P波初至，例如长短时窗均值比法、自回归池赤准则法、峰度法、偏度法、分形维法、卓越周期法和人工神经网络法等。上述P波初至绝对到时拾取方法利用单个波形振幅的起跳点，然而受背景噪声和波场传播的波前愈合效应等影响，P波初至绝对振幅可能很不清晰或淹没于背景噪声之中，由此造成P波初至绝对到时拾取困难。
为解决上述问题，一些研究者开始采用波形互相关法测量相对到时来替代传统的P波初至绝对到时拾取，其对加窗波形段进行互相关操作，取互相关最大值作为对比信号的到时差。因此，波形互相关法主要取决于加窗波形段优势频段主峰的位置，受背景噪声的影响较P波初至绝对到时拾取要小很多，且研究表明波形互相关法对于地震波相对到时拾取非常有效。
上述P波到时拾取方法均建立在波形不存在P波初至系统观测误差的假定上，然而传感器场地效应和中转设备处理响应时间等可能导致稳定的P波初至系统观测误差；那么观测波形的P波到时应该等于震源实际发生的时刻、P波由震源到传感器的传播时间和P波初至系统观测误差之和，由此声发射事件定位前需校正P波初至系统观测误差。
震源定位目标函数的收敛性常和选择的最优化反演算法息息相关。Geiger定位法采用了一阶梯度下降法，计算快捷、但易受初始值选择的影响，而Thurber和Li et al分别采用包含二阶偏导数的牛顿法和高斯牛顿法对目标函数进行求解，提高了算法的稳定性，但二阶hessian矩阵计算量巨大；Prugger and Gendzwill和Li et al将单纯形法引入到震源定位中，具有较好的计算速度和定位精度。上述梯度基算法对于初始点的依赖性很强，容易陷入局部最优解，定位稳定性较差。
发明内容
本发明提出了一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质，该方法针对现有P波初至拾取方法不稳定，提出用《一种矿山微震信号P波初至时刻联合拾取方法》自动拾取P波初至，并采用人工方法校正拾取不好和无P波初至拾取的信号，降低了错误拾取率，增强了自动拾取P波初至的稳定性。利用Bayes建立基于P波初至理论与观测数据差的P波初至系统观测误差反演模型，有助于高维数据反演。最后，利用Bayes、MCMC方法对校正后的P波初至数据展开定位，极大的提高了定位精度。
本发明提供了以下技术方案：
一方面，一种去除P波初至系统观测误差的声发射事件Bayes定位方法，其特征在于，包括以下步骤：
步骤1：采用低通滤波器对声发射信号高频噪音降噪；
步骤2：拾取声发射信号P波初至观测到时；
步骤3：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
步骤4：构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
所述参数模型包括的参数为声发射事件位置、发震时刻以及P波初至系统观测误差；
所述利用MCMC对参数模型中的参数采样，是指改变参数模型中的任意一个参数，获得新的参数模型；
步骤5：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
步骤6：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
采用Bayes定位方法易于得到全局最优解，能够给出声发射事件定位结果的可靠性评价，迭代后期各参数趋于稳定，则结果可靠；若迭代后期各参数震荡范围较大，则结果可靠性有限。
进一步地，使用《一种矿山微震信号P波初至时刻联合拾取方法》自动拾取P波初至，并采用人工方法校正拾取不好和无P波初至拾取的信号；
其中，《一种矿山微震信号P波初至时刻联合拾取方法》具体过程如下：提取微震信号波形数据；最大峰度值到时点(K1)拾取；最小自相关值到时点(K2)拾取；判定∣K1-K2∣是否大于阈值L，若是则舍弃该段微震信号，反之则输出P波到时K＝K2；L取值范围为3-5ms。
进一步地，所述基于声发射事件位置和P波初至系统观测误差的理论到时与观测到时差的目标函数如下：
其中，i为声发射事件个数，j为传感器编号，
为第i个声发射事件的发震时刻，
为第i个声发射事件震源到j号传感器的传播时间，Tj为j号传感器P波初至系统观测误差，tij为观测数据中的第i个声发射事件j号传感器的P波初至波形拾取时刻，参与P波初至系统观测误差反演的传感器数为n1个、声发射事件数为n2个。
假定IMS监测系统一直处于稳定工作状态，即声发射信号P波初至系统观测误差在研究时段保持不变。
设定参与P波初至系统观测误差反演的传感器数为n1个、声发射事件数为n2个、P波初至数据为n3个，每个事件的三维坐标
和发震时刻
未知，则总未知数为4n2+n1个。又因求解方程时总未知数的个数需小于或等于P波初至数据个数，即4n2+n1≤n3。
上述目标函数求解的方式多种多样，然而该问题未知数维数达到了4n2+n1个，为此本申请将建立基于公式(1)的高维全局最优解求解的Bayes定位方法。具体建模过程如下：记观测数据中的P波初至到时dij为第i个声发射事件j号传感器的P波初至波形拾取时刻tij(dij＝tij)，理论模型数据中的P波初至理论到时
为第i个声发射事件j号传感器的P波初至理论到时为
进一步地，所述基于理论数据与观测数据之差的Bayes定位后验概率密度函数：
其中，
为P波初至到时观测数据dij依次排开所组成的一个列向量，
为P波初至理论到时
依次排开所组成的一个列向量，
为表征
拟合程度的协方差矩阵；
为关于声发射事件位置、发震时刻以及P波初至系统观测误差的参数模型。
对于某一具体参数模型
时，可由该模型的声发射事件震源位置
计算第i个声发射事件震源到j号传感器的传播时间
再结合该模型的
和Tj参数，可计算得到
的值，所有声发射事件和传感器的理论到时构成理论模型数据
又
为已知的观测数据，由此可计算得到公式(2)中的
进一步地，由射线理论计算所述Bayes定位后验概率密度函数中的第i个声发射事件震源到j号传感器的传播时间
再用马尔科夫链蒙特卡罗(Markov Chain MonteCarlo,MCMC)采样并结合
来计算声发射事件震源位置和P波初至系统观测误差；
当
时，接受迭代计算中更新参数后的新模型
当
时，拒绝该新模型，并保持为模型
其中，u为0～1的均匀分布随机数据，每次迭代时均需重新生成一个u值。
在每次迭代时，若更新参数后得到的新模型不被接受时，利用前一个被接受的参数模型作为当前迭代获得的参数模型；
进一步地，MCMC采样时每次随机更新声发射事件震源位置、发震时刻和P波初至系统观测误差中的一个参数：
当选中的参数为声发射事件震源位置中的一个参数时，更新速度为M*g；当选中的参数为发震时刻或P波初至系统观测误差时，更新速度为N*g；
其中，g指均值为0、方差为1的标准正态分布，M取值范围为1～3mm，N为1us。
MCMC采样更新某个参数后，其余参数保持当前模型参数，计算这些参数对应的
进而得到
和
由公式(3)判断是否接受新模型。
所述更新速度是指选中的参数更新后的取值，在每次迭代时，g都会生成新的取值；
进一步地，取各参数最后5000～10000次迭代的均值作为各参数反演结果。
MCMC迭代的后期各参数趋于稳定；
进一步地，使用Matlab导入声发射信号，利用经验模态分解将声发射信号自适应分解为多个模态分量，去除前两个模态分量，实现声发射信号高频噪音降噪。
再一方面，一种去除P波初至系统观测误差的声发射事件Bayes定位系统，包括：
声发射信号去噪单元：采用低通滤波器对声发射信号进行去噪；
P波初至时刻拾取单元：用于拾取声发射信号P波初至观测到时；
目标函数构建单元：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
P波初至系统观测误差获取单元：通过构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
单个声发射事件校正单元：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
单个声发射事件定位单元：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
再一方面，一种可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行一种去除P波初至系统观测误差的声发射事件Bayes定位的方法。
有益效果
本发明提供了一种去除P波初至系统观测误差的声发射事件Bayes定位方法、系统及介质，主要用于去除声发射事件定位中P波初至系统观测误差，解决传统定位中受P波初至系统观测误差的影响，Bayes反演有利于实现高维数据反演和全局最优解，并给出定位结果可靠性评价。本方法包括如下步骤：首先采用低通滤波器对声发射信号高频噪音降噪；再用自动法+人工拾取确定声发射信号P波初至；建立基于P波初至理论与观测数据差的Bayes定位后验概率密度函数；进而采用马尔科夫链蒙特卡罗(MCMC)采样方法来迭代求解P波初至系统观测误差；利用P波初至拾取到时tij加上P波初至系统观测误差Tj，得到校正后的P波初至数据；进而利用Bayes定位后验概率密度函数、MCMC方法对校正后的P波初至数据开展定位；最后输出定位结果，并评价定位结果的可靠性。该方法拾取P波初至数据较快、精度高，可去除P波初至系统观测误差，提高了P波初至数据的准确性，Bayes定位方法易于得到全局最优解，并且能够给出声发射事件定位结果的可靠性评价。
附图说明
图1是本发明实例所述方法流程图；
图2是传感器和断铅测试事件位置及坐标图；
图3是1-8号传感器接收断铅实验得到的声发射信号波形图；
图4是1-8号传感器声发射信号波形局部放大后的P波初至拾取图；
图5是校正系统观测误差前后Bayes方法定位效果图，其中，(a)为包含8号传感器且未校正P波初至系统观测误差的Bayes定位效果图，(b)为包含8号传感器且校正P波初至系统观测误差的Bayes定位效果图，(c)为不包含8号传感器且未校正P波初至系统观测误差的Bayes定位效果图，(d)为不包含8号传感器且校正P波初至系统观测误差的Bayes定位效果图。
具体实施方式
本发明实施例所述的方法思想如下：声发射定位时需要考虑P波初至系统观测误差，倘若将其直接作为未知数对单个事件开展震源定位，那么再加上未知数(x0,y0,z0,t0)，总未知数的个数将多余该事件P波初至数据个数，为此需要预先去除P波初至系统观测误差再开展定位。首先随机选取多个事件反演得到这些事件的震源位置、发生时刻和传感器P波初至系统观测误差。进而利用Bayes方法对去除P波初至系统观测误差的单个待定位新声发射事件进行定位。
如图1所示，一种去除P波初至系统观测误差的声发射事件Bayes定位方法，包括以下步骤：
步骤1：采用低通滤波器对声发射信号高频噪音降噪；
使用Matlab导入声发射信号，利用经验模态分解将声发射信号自适应分解为多个模态分量，去除前两个模态分量，实现声发射信号高频噪音降噪。
步骤2：拾取声发射信号P波初至观测到时；
步骤3：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
步骤4：构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
所述参数模型包括的参数为声发射事件位置、发震时刻以及P波初至系统观测误差；
所述利用MCMC对参数模型中的参数采样，是指改变参数模型中的任意参数，获得新的参数模型；
步骤5：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
步骤6：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
采用Bayes定位方法易于得到全局最优解，能够给出声发射事件定位结果的可靠性评价，迭代后期各参数趋于稳定，则结果可靠；若迭代后期各参数震荡范围较大，则结果可靠性有限。
使用《一种矿山微震信号P波初至时刻联合拾取方法》自动拾取P波初至，并采用人工方法校正拾取不好和无P波初至拾取的信号；
《一种矿山微震信号P波初至时刻联合拾取方法》具体过程如下：提取微震信号波形数据；最大峰度值到时点(K1)拾取；最小自相关值到时点(K2)拾取；判定∣K1-K2∣是否大于阈值L，若是则舍弃该段微震信号，反之则输出P波到时K＝K2；L取值范围为3-5ms。
所述基于声发射事件位置和P波初至系统观测误差的理论到时与观测到时差的目标函数如下：
其中，i为声发射事件个数，j为传感器编号，
为第i个声发射事件的发震时刻，
为第i个声发射事件震源到j号传感器的传播时间，Tj为j号传感器P波初至系统观测误差，tij为观测数据中的第i个声发射事件j号传感器的P波初至波形拾取时刻，参与P波初至系统观测误差反演的传感器数为n1个、声发射事件数为n2个。
假定IMS监测系统一直处于稳定工作状态，即声发射信号P波初至系统观测误差在研究时段保持不变。
设定参与P波初至系统观测误差反演的传感器数为n1个、声发射事件数为n2个、P波初至数据为n3个，每个事件的三维坐标
和发震时刻
未知，则总未知数为4n2+n1个，又因求解方程时总未知数的个数需小于或等于P波初至数据个数，即4n2+n1≤n3。
上述目标函数求解的方式多种多样，然而该问题未知数维数达到了4n2+n1个，为此本申请将建立基于公式(1)的高维全局最优解求解的Bayes定位方法，具体建模过程如下：记观测数据中的P波初至到时dij为第i个声发射事件j号传感器的P波初至波形拾取时刻tij(dij＝tij)，理论模型数据中的P波初至理论到时
为第i个声发射事件j号传感器的P波初至理论到时为
所述基于理论数据与观测数据之差的Bayes定位后验概率密度函数：
其中，
为P波初至到时观测数据dij依次排开所组成的一个列向量，
为P波初至理论到时
依次排开所组成的一个列向量，
为表征
拟合程度的协方差矩阵；
为关于声发射事件位置、发震时刻以及P波初至系统观测误差的参数模型。
对于某一具体参数模型
时，可由该模型的声发射事件震源位置
计算第i个声发射事件震源到j号传感器的传播时间
再结合该模型的
和Tj参数，可计算得到
的值，所有声发射事件和传感器的理论到时构成理论模型数据
又
为已知的观测数据，由此可计算得到公式(2)中的
由射线理论计算所述Bayes定位后验概率密度函数中的第i个声发射事件震源到j号传感器的传播时间
再用马尔科夫链蒙特卡罗(Markov Chain Monte Carlo,MCMC)采样并结合
和来计算声发射事件震源位置和P波初至系统观测误差；
当
时，接受迭代计算中更新参数后的新模型
当
时，拒绝该新模型，并保持为模型
其中，u为0～1的均匀分布随机数据，每次迭代时均需重新生成一个u值。
MCMC采样时每次随机更新声发射事件震源位置、发震时刻和P波初至系统观测误差中的一个参数：
当选中的参数为声发射事件震源位置中的一个参数时，更新速度为M*g；当选中的参数为发震时刻或P波初至系统观测误差时，更新速度为N*g；
其中，g指均值为0、方差为1的标准正态分布，M取值范围为1～3mm，N为1us。
所述更新速度是指选中的参数更新后的取值，在每次迭代时，g都会生成新的取值；
MCMC采样更新某个参数后，其余参数保持当前模型参数，计算这些参数对应的
进而得到
和
由公式(3)判断是否接受新模型。
各参数最后5000～10000次迭代的均值作为各参数反演结果，MCMC迭代的后期各参数趋于稳定。
实施例
图2是声发射实验各传感器和断铅测试事件的位置图，图中三角形代表传感器，五角星为断铅实验位置。各传感器的具体位置坐标如下表1所示，断铅实验(X,Y,Z)的位置坐标为(50,10,40)mm。
表1 声发射实验传感器布设位置坐标
图3是各传感器接收断铅实验得到的声发射信号波形图。由图可知：在观测时间1024us内，各传感器接收的声发射信号波形相似，另外实现发现8号传感器总是早约260us于其它传感器接收到波形，基于上述分析8号声发射信号与信号为同一声发射事件产生，而50mm×50mm×100mm的砂岩事件声发射信号P波传播时间不会超过50us，有理由相信声发射实验观测到了明显的P波初至系统观测误差，而其他传感器是否存在P波初至系统观测误差需要通过计算验证。
图4是图3中各传感器声发射信号P波初至段波形放大后的P波初至拾取结果图，进而利用公式(1)和(2)的Bayes后验概率目标函数和MCMC采样计算得到各传感器P波初至系统观测误差如表2所示，由表知各传感器确实存在P波初至系统观测误差，进一步地，可通过断铅实验确定去除P波初至系统观测误差后的定位效果。
表2 Bayes反演得到的各传感器P波初至系统观测误差
| 传感器编号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| P波初至系统观测误差(us) | -6.5 | -0.2 | 1.1 | 0.2 | -3.3 | 0.7 | -0.8 | 254.0 |
图5是校正P波初至系统观测误差前后Bayes方法定位效果图，其中，(a)为包含8号传感器且未校正P波初至系统观测误差的Bayes定位效果图，(b)为包含8号传感器且校正P波初至系统观测误差的Bayes定位效果图，(c)为不包含8号传感器且未校正P波初至系统观测误差的Bayes定位效果图，(d)为不包含8号传感器且校正P波初至系统观测误差的Bayes定位效果图。由图5(a)～(d)最后10000次MCMC迭代的平均值得到Bayes反演的定位结果如表3所示。由图表可知，包含8号传感器且未校正P波初至系统观测误差时，Bayes定位结果不收敛，而不包含8号传感器且未校正P波初至系统观测误差时，Bayes定位结果收敛较差，表明8号传感器的P波初至系统观测误差对定位结果影响很大和P波初至系统观测误差对定位结果具有较大影响，而校正P波初至系统观测误差后的Bayes定位效果很好，再次证明P波初至系统观测误差校正对声发射事件定位具有重要意义。
表3 校正P波初至系统观测误差前后Bayes方法定位效果图
基于上述方法，本发明实例还提供了一种去除P波初至系统观测误差的声发射事件Bayes定位系统，包括：
声发射信号去噪单元：采用低通滤波器对声发射信号进行去噪；
P波初至时刻拾取单元：用于拾取声发射信号P波初至观测到时；
目标函数构建单元：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
P波初至系统观测误差获取单元：通过构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
单个声发射事件校正单元：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
单个声发射事件定位单元：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。
一种可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行所述的一种去除P波初至系统观测误差的声发射事件Bayes定位方法。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详尽的说明，所属领域的普通技术人员应当理解，上述实施例仅仅是对本发明的示意性实现方式的解释，实施例中的细节并不构成对本发明范围的限制，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案的等效变换、简单替换等显而易见的改变，均落在本发明保护范围之内。
Claims (10)
1.一种去除P波初至系统观测误差的声发射事件Bayes定位方法，其特征在于，包括以下步骤：
步骤1：采用低通滤波器对声发射信号高频噪音降噪；
步骤2：拾取声发射信号P波初至观测到时；
步骤3：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
步骤4：构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
所述参数模型包括的参数为声发射事件位置、发震时刻以及P波初至系统观测误差；
步骤5：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
步骤6：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
2.根据权利要求1所述的方法，其特征在于，使用《一种矿山微震信号P波初至时刻联合拾取方法》自动拾取P波初至，并采用人工方法校正拾取不好和无P波初至拾取的信号；
其中，《一种矿山微震信号P波初至时刻联合拾取方法》具体过程如下：提取微震信号波形数据；最大峰度值到时点K1拾取；最小自相关值到时点K2拾取；判定∣K1-K2∣是否大于阈值L，若是则舍弃该段微震信号，反之则输出P波到时K＝K2；L取值范围为3-5ms。
3.根据权利要求1所述的方法，其特征在于，所述基于声发射事件位置和P波初至系统观测误差的理论到时与观测到时差的目标函数如下：
其中，i为声发射事件个数，j为传感器编号，
为第i个声发射事件的发震时刻，
为第i个声发射事件震源到j号传感器的传播时间，Tj为j号传感器P波初至系统观测误差，tij为观测数据中的第i个声发射事件j号传感器的P波初至波形拾取时刻，参与P波初至系统观测误差反演的传感器数为n1个、声发射事件数为n2个。
4.根据权利要求1所述的方法，其特征在于，所述基于理论数据与观测数据之差的Bayes定位后验概率密度函数：
其中，
为P波初至到时观测数据dij依次排开所组成的一个列向量，
为P波初至理论到时
依次排开所组成的一个列向量，
为表征
拟合程度的协方差矩阵；
为关于声发射事件位置、发震时刻以及P波初至系统观测误差的参数模型。
5.根据权利要求4所述的方法，其特征在于，由射线理论计算所述Bayes定位后验概率密度函数中的第i个声发射事件震源到j号传感器的传播时间
再利用马尔科夫链蒙特卡罗采样并结合
来计算声发射事件震源位置和P波初至系统观测误差；
当
时，接受迭代计算中更新参数后的新模型
当
时，拒绝该新模型，并保持为模型
其中，u为0～1的均匀分布随机数据，每次迭代时均需重新生成一个u值。
6.根据权利要求1所述的方法，其特征在于，MCMC采样时每次随机更新声发射事件震源位置、发震时刻和P波初至系统观测误差中的一个参数：
当选中的参数为声发射事件震源位置中的一个参数时，更新速度为M*g；当选中的参数为发震时刻或P波初至系统观测误差时，更新速度为N*g；
其中，g指均值为0、方差为1的标准正态分布，M取值范围为1～3mm，N为1us。
7.根据权利要求1所述的方法，其特征在于，取各参数最后5000～10000次迭代的均值作为各参数反演结果。
8.根据权利要求1-7任一项所述的方法，其特征在于，使用Matlab导入声发射信号，利用经验模态分解将声发射信号自适应分解为多个模态分量，去除前两个模态分量，实现声发射信号高频噪音降噪。
9.一种去除P波初至系统观测误差的声发射事件Bayes定位系统，其特征在于，包括：
声发射信号去噪单元：采用低通滤波器对声发射信号进行去噪；
P波初至时刻拾取单元：用于拾取声发射信号P波初至观测到时；
目标函数构建单元：建立基于声发射事件位置和P波初至系统观测误差的理论到时与P波初至观测到时之差的目标函数；
P波初至系统观测误差获取单元：通过构建参数模型与P波初至系统观测误差的Bayes定位后验概率密度函数，并利用MCMC对参数模型中的任意参数采样，生成新的参数模型后，计算随机获取的声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，获得P波初至系统观测误差；
单个声发射事件校正单元：对待求的单个声发射事件的P波初至拾取到时加上P波初至系统观测误差，得到校正后的P波初至到时；
单个声发射事件定位单元：将校正后的P波初至到时代入参数模型中，生成新的参数模型后，利用MCMC对参数模型中的声发射事件位置、发震时刻采样，计算待求解的单个声发射事件的所述Bayes定位后验概率密度函数，判断采样更新后的参数模型是否可接受，通过迭代，利用参数模型在迭代稳定期的均值，定位待求解的单个声发射事件位置。
10.一种可读存储介质，包括计算机程序指令，其特征在于：所述计算机程序指令被处理终端执行时使所述处理终端执行一种去除P波初至系统观测误差的声发射事件Bayes定位的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010601258.2A CN111722280B (zh) | 2020-06-29 | 2020-06-29 | 一种去除p波初至系统观测误差的声发射事件定位方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010601258.2A CN111722280B (zh) | 2020-06-29 | 2020-06-29 | 一种去除p波初至系统观测误差的声发射事件定位方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111722280A true CN111722280A (zh) | 2020-09-29 |
| CN111722280B CN111722280B (zh) | 2021-09-07 |
Family
ID=72569538
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010601258.2A Active CN111722280B (zh) | 2020-06-29 | 2020-06-29 | 一种去除p波初至系统观测误差的声发射事件定位方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111722280B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112731523A (zh) * | 2020-12-22 | 2021-04-30 | 北京环境特性研究所 | 一种震源定位方法、计算机设备及计算机可读存储介质 |
| CN116466390A (zh) * | 2023-02-17 | 2023-07-21 | 南方海洋科学与工程广东省实验室(广州) | 一种大型水库诱发地震实时监测定位方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2001096904A2 (en) * | 2000-06-12 | 2001-12-20 | Schlumberger Technology Corporation | Method for downhole waveform tracking and sonic labeling |
| CN102879801A (zh) * | 2012-08-30 | 2013-01-16 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 一种基于射孔约束的EnKF微地震事件位置反演方法 |
| CN104914468A (zh) * | 2015-06-09 | 2015-09-16 | 中南大学 | 一种矿山微震信号p波初至时刻联合拾取方法 |
| CN107861149A (zh) * | 2017-11-03 | 2018-03-30 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于波形驱动下的叠前纵横波速度比模拟方法 |
| CN108663711A (zh) * | 2018-04-04 | 2018-10-16 | 电子科技大学 | 一种基于τ分布的贝叶斯地震反演方法 |
| WO2020081156A1 (en) * | 2018-10-16 | 2020-04-23 | Halliburton Energy Services, Inc. | Method to improve das channel location accuracy using global inversion |
| CN111175815A (zh) * | 2020-01-06 | 2020-05-19 | 中国石油化工股份有限公司 | 油藏改造微地震监测裂缝震源机制求解方法及系统 |
-
2020
- 2020-06-29 CN CN202010601258.2A patent/CN111722280B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2001096904A2 (en) * | 2000-06-12 | 2001-12-20 | Schlumberger Technology Corporation | Method for downhole waveform tracking and sonic labeling |
| CN102879801A (zh) * | 2012-08-30 | 2013-01-16 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 一种基于射孔约束的EnKF微地震事件位置反演方法 |
| CN104914468A (zh) * | 2015-06-09 | 2015-09-16 | 中南大学 | 一种矿山微震信号p波初至时刻联合拾取方法 |
| CN107861149A (zh) * | 2017-11-03 | 2018-03-30 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于波形驱动下的叠前纵横波速度比模拟方法 |
| CN108663711A (zh) * | 2018-04-04 | 2018-10-16 | 电子科技大学 | 一种基于τ分布的贝叶斯地震反演方法 |
| WO2020081156A1 (en) * | 2018-10-16 | 2020-04-23 | Halliburton Energy Services, Inc. | Method to improve das channel location accuracy using global inversion |
| CN111175815A (zh) * | 2020-01-06 | 2020-05-19 | 中国石油化工股份有限公司 | 油藏改造微地震监测裂缝震源机制求解方法及系统 |
Non-Patent Citations (2)
| Title |
|---|
| JIULONG CHENG 等: "Research Developments and Prospects on Microseismic Source Location in Mines", 《ENGINEERING》 * |
| 谭玉阳 等: "一种改进的基于网格搜索的微地震震源定位方法", 《地球物理学报》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112731523A (zh) * | 2020-12-22 | 2021-04-30 | 北京环境特性研究所 | 一种震源定位方法、计算机设备及计算机可读存储介质 |
| CN112731523B (zh) * | 2020-12-22 | 2022-09-16 | 北京环境特性研究所 | 一种震源定位方法、计算机设备及计算机可读存储介质 |
| CN116466390A (zh) * | 2023-02-17 | 2023-07-21 | 南方海洋科学与工程广东省实验室(广州) | 一种大型水库诱发地震实时监测定位方法 |
| CN116466390B (zh) * | 2023-02-17 | 2023-11-03 | 南方海洋科学与工程广东省实验室(广州) | 一种大型水库诱发地震实时监测定位方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111722280B (zh) | 2021-09-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111722280B (zh) | 一种去除p波初至系统观测误差的声发射事件定位方法 | |
| CN110907991B (zh) | 基于数据场势值的震源定位方法、系统及可读存储介质 | |
| CN107272066B (zh) | 一种含噪地震信号初至走时拾取方法及装置 | |
| CN109375252B (zh) | 考虑不同发震构造最大可信地震的地震动参数评价方法 | |
| CN112014883B (zh) | 一种基于Log-Cosh函数的微震震源定位方法、系统、装置及可读存储介质 | |
| CN114994754B (zh) | 基于直达波和深度震相初动极性的震源机制联合反演方法 | |
| CN110703313B (zh) | 考虑传感器感度的声发射事件震级获取方法、系统及可读存储介质 | |
| CN112731534B (zh) | 一种考虑p波初至系统误差的双声发射事件联合定位方法、系统、电子设备及可读介质 | |
| CN113239022A (zh) | 医疗诊断缺失数据补全方法及补全装置、电子设备、介质 | |
| CN112526602A (zh) | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 | |
| CN108469609B (zh) | 一种用于雷达目标跟踪的检测信息滤波方法 | |
| JP2007081455A (ja) | 音源位置・受音位置推定方法、その装置、そのプログラム、およびその記録媒体 | |
| CN112183407B (zh) | 一种基于时频域谱减法的隧道地震波数据去噪方法及系统 | |
| CN111736208B (zh) | 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质 | |
| CN113866827B (zh) | 一种解释性速度建模地震成像方法、系统、介质和设备 | |
| CN111863016B (zh) | 一种天文时序信号的噪声估计方法 | |
| CN111624655B (zh) | 初至波剩余静校正量的确定方法及装置 | |
| CN112068089B (zh) | 一种基于粒子滤波的序列检索方法 | |
| Li et al. | Three dimensional acoustic source localization and tracking using statistically weighted hybrid particle filtering algorithm | |
| CN111308423A (zh) | 鲁棒声源定位系统及其方法 | |
| CN111596357B (zh) | 一种分析海底采集节点工作状态的方法及装置 | |
| CN112731523B (zh) | 一种震源定位方法、计算机设备及计算机可读存储介质 | |
| CN113219534B (zh) | 一种叠前深度偏移速度质控方法、装置、介质及电子设备 | |
| CN112083480B (zh) | 同步激发采集地震数据的实时监控方法及装置 | |
| CN115184858A (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 |