CN105527650B - 一种工程尺度下微震信号及p波初至自动识别算法 - Google Patents
一种工程尺度下微震信号及p波初至自动识别算法 Download PDFInfo
- Publication number
- CN105527650B CN105527650B CN201610090278.1A CN201610090278A CN105527650B CN 105527650 B CN105527650 B CN 105527650B CN 201610090278 A CN201610090278 A CN 201610090278A CN 105527650 B CN105527650 B CN 105527650B
- Authority
- CN
- China
- Prior art keywords
- point
- mrow
- sampled point
- current sampling
- amplitude
- 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
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 21
- 238000005070 sampling Methods 0.000 claims description 95
- 238000004364 calculation method Methods 0.000 claims description 3
- 239000011435 rock Substances 0.000 abstract description 9
- 238000000034 method Methods 0.000 description 11
- 230000008859 change Effects 0.000 description 8
- 238000012544 monitoring process Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 5
- 230000035772 mutation Effects 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 238000010276 construction Methods 0.000 description 3
- 238000010206 sensitivity analysis Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000002045 lasting effect Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000006854 communication Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000007774 longterm Effects 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
- 238000005457 optimization Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012421 spiking 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/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
本发明公开了一种工程尺度下微震信号及p波初至自动识别算法，包括以下步骤，读取设定时窗内采样点数据；计算采样点的特征函数、长短时平均值；进行微震信号及其P波初至自动拾取判断；参数值缓存清零；读取下一段时窗内数据。本发明针对工程尺度下微震信号特点，改进地震领域的识别算法，本发明能够准确的自动拾取工程尺度下的微震信号特别是性噪比较低的弱信号，同时提高了自动拾取P波初至的准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性。
Description
技术领域
本发明涉及微震监测技术领域。具体涉及一种工程尺度下微震信号及p波初至自动识别算法，该方法可广泛用于矿业工程、水利水电工程、石油工程、岩土工程以及地下工程。
背景技术
微震源定位是微震监测与灾害预警的重要组成部分，而微震信号的识别及P波初至拾取是微震源定位的关键。现有的微震信号识别方法大部分来源于地震领域，方法很多，主要有：根据在时间域能量和能量变化构建特征函数时间域的STA/LTA算法；根据地震信号到达前后地震波形数据统计性质的差别，如AIC算法；此外还有神经网络法、高阶统计法、以及基于小波理论的小波变换法等。STA/LTA及其改进算法，由于算法简单、计算效率高、适于实时处理，在地震领域被广泛应用。这些方法用于地震信号的处理具有较好的适用性，但工程尺度下产生的岩石破裂信号与天然地震信号不同：
1)地震信号频率一般为几十赫兹以下，而微震信号的频率一般为几十到几百赫兹，有的高达几千赫兹；2)地震信号持续时间较长，一般大于1秒，微震信号持续时间较短，一般不超过0.1秒；3)微震与天然地震相比，一般震级小，较小的微震事件不易被拾取；4)工程环境下由于机械施工、电气干扰、车辆行驶等因素造成背景噪音复杂，使得监测到的大多数微震信号信噪比较低。
因此，将地震领域的信号识别及P波到时拾取方法用于工程尺度微震信号是不合适的，主要表现为：1)难以识别低信噪比的微弱信号；2)P波初至自动拾取误差较大。
为此，针对上述不足，发明一种可以有效识别工程尺度下微震信号及其p波初至的自动识别的方法，提高工程尺度下微震信号及其P波初至的拾取准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性，是非常必要的。
发明内容
本发明的目的在于克服现有技术存在的问题，提供了一种工程尺度下微震信号及p波初至自动识别算法，提高工程尺度下微震信号及其P波初至的拾取准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性。
一种工程尺度下微震信号及p波初至自动识别算法，包括以下步骤：
步骤1、读取振动传感器实时监测到的设定时窗内波形的采样点数据，建立平面直角坐标系，其中X轴对应波形的采样点编号，Y轴对应波形的幅值，在平面直角坐标系中对时窗内波形的采样点进行对称校准；
步骤2、计算采样点个数N，设置零交叉点个数上限M1，零交叉点个数下限M2；设置微震信号时间判断阀值T；
步骤3、设定第A个采样点的加权因子K(A)和特征函数CF(A)；
设定第A个采样点短时平均值STA(A)和长时平均值LTA(A)；
设定第A个采样点的动态阈值r(A)；
其中A为采样点个数的编号；
步骤4、以平面直角坐标系中编号最小的采样点为当前采样点，当前采样点的编号为i；
步骤5、设置算法参数变量M、S、L、t，M、S、L、t的初值均设置为0，其中M为零交叉点个数；S为计数器；L的计算公式为L＝3+M/3；t为拾取出的微震信号的持续时间；
步骤6、比较当前采样点的短时平均值STA(i)和动态阈值r(i)的大小；
步骤7、若当前采样点的短时平均值STA(i)小于当前采样点的动态阀值r(i)，则当前采样点不是P波初至点，此时将当前采样点的下一个采样点作为当前采样点，返回到步骤(6)；若当前采样点的短时平均值STA(i)大于等于当前点的动态阀值r(i)，则将当前采样点的编号i赋值给k，计算当前采样点的幅值P1，将当前采样点的幅值大小P1赋值给缓存最大幅值P0，然后进行步骤(8)；
步骤8、计算编号为k的采样点的下一个采样点k+1的幅值P2；
步骤9、如果采样点的编号为k的采样点的下一个采样点k+1不是零交叉点，将缓存最大幅值P0和采样点k+1的幅值P2中较大值赋值给缓存最大幅值P0，将k+1赋值给k，当k的值小于波形采样点个数N时，返回步骤8；当k的值大于等于波形采样点个数N时，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；
如果采样点的编号为k的下一个采样点k+1是零交叉点，则零交叉点个数M值加1，将k+1赋值给k，然后进行步骤10；
步骤10、若零交叉点个数M值大于零交叉点个数上限M1值，则当前采样点不为P波初至点，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；
若零交叉点个数M值小于等于零交叉点个数上限M1值，则计算当前采样点的判定参数δ＝LTA(i)×M，然后进行步骤11；
步骤11、比较编号为k的采样点的短时平均值STA(k)与当前采样点的判定参数δ的大小：
若编号为k的采样点的短时平均值STA(k)小于等于当前采样点的判定参数δ，计数器S值归零，将P2赋值给P0，此时将k加1，返回步骤8；
若编号为k的采样点的短时平均值STA(k)大于当前采样点的判定参数δ，计数器S值加1，计算L＝3+M/3，进入步骤12；
步骤12、比较计数器S值与L的大小：
若计数器S值小于等于L，则将P2赋值给P0，将k加1，返回步骤8；
若计数器S值大于L，计算从当前采样点i到编号为k的采样点的这段波形的持续时间t，t＝(k-i)/f，其中f为采样频率；
步骤13、当t小于微震信号时间判断阀值T或零交叉点个数M小于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点不是微震信号，将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5；当t大于等于时间阀值T且零交叉点个数M大于等于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点为微震信号，将当前采样点的编号i赋值给T1，将k赋值给T2，其中编号为T1的采样点为微震信号的P波初至点，编号为T2的采样点为信号结束点，进入步骤14；
步骤14、将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5，直至将设定时窗内波形的采样点数据分析完。
如上所述的对称校准包括以下步骤：
选取时窗中设定个数的连续的采样点作为校正样本，
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据需对称校准，计算校正样本幅值的平均值，将设定时窗内波形的采样点幅值与校正样本幅值的平均值相加，完成对称校准；
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据不需对称校准。
如上所述的第A个采样点的加权因子K(A)基于以下公式：
所述的第A个采样点的特征函数CF(A)基于以下公式：
其中，A为对称校准后波形中采样点个数的编号，y(A)为对称校准后波形的第A个采样点的幅值，为y(A)一阶差分；
所述的第A个采样点的短时平均值STA(A)基于以下公式：
STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]
所述的第A个采样点的长时平均值LTA(A)基于以下公式：
LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]
所述的第A个采样点的动态阈值r(A)基于以下公式：
r(A)＝C5×LTA(A)
其中，C3为短时平均系数，C4为长时平均系数，C5为触发阀值。
当P波到来时，将引起微震信号振幅或者频率等波形特征变化，若信号的信噪比越低，则P波初至越不明显。因而提高拾取算法对该类变化的敏感性，将能有利于提高工程尺度下微震信号的P波初至拾取能力。而算法中的加权因子K、特征函数CF、动态长短时平均比值ε均能反应波形的特征变化，直接影响拾取准确率，动态长短时平均比值ε为长时平均值LTA(A)与短时平均值STA(A)的比值：
分析本发明的加权因子K、特征函数CF、动态长短时平均比值ε对信号振幅或者频率波形等特征变化的敏感性，分别如图2、3、4。对比结果分析本发明三种变量在频率或振幅变化处数值改变明显，达到峰值的迅速，因此在理论上本发明对信噪比低的信号及其P波初至的拾取效果明显。
从锦屏引水隧洞实时监测的数据中随机选取555个微震信号为样本，定义两个拾取结果判断准则：1)信号成功拾取：能将样本信号自动识别出来而不产生漏波与误判。这样将有利于人工对自动拾取出来的信号部分进行进一步处理。信号成功拾取个数越多，信号拾取率越高。2)P波初至拾取准确：将自动算法拾取P波初至结果与人工拾取结果相比，若拾取误差不大于3个采样点，则认为该算法自动拾取结果准确。本发明信号拾取率为94.23％，P波拾取准确率为73.51％，拾取率完全满足微震领域目前定位误差。
分析2015年6月15日至2015年7月15日锦屏极深地下实验室微震监测数据，一共524个微震信号，109个事件，将人工拾取P波、本发明拾取结果用于定位，结果分别如图6、7所示，其中7#实验室是岩石破裂高风险区域，可以看出本发明定位结果微震事件的聚集度较高，微震事件空间分布规律与人工拾取结果接近，与工程实际情况相符。将人工拾取过程中信噪比非常低、P波不明显、定位结果明显异常的13个事件，如图6中的三角形所示，用本发明拾取定位分析：1)本发明拾取的P波结果均能成功定位，如图7中的三角形所示；2)本发明该13个事件拾取定位效果相对于人工拾取改进明显，绝大部分在定位误差允许的范围内。
有益效果：
本发明针对工程尺度下微震信号特点，改进地震领域的识别算法，本发明能够准确的自动拾取工程尺度下的微震信号特别是性噪比较低的弱信号，同时提高了自动拾取P波初至的准确率，进而提高岩爆、矿震、塌方等地质灾害预警的及时性与准确性。
附图说明
图1为本发明流程图；
图2为本发明加权因子K对振幅或频率变化敏感性分析图；如图2(a)当振幅发生变化时，本发明的加权因子振幅约增大2倍；如图2(b)当频率发生变化时，本发明K值变化幅度较为明显。
图3为本发明特征函数CF对振幅或频率变化敏感性分析图；图3(a)模拟背景噪音中出现幅值较小的微震信号，即幅值由1突变到2时，本发明特征函数在突变点的幅值由1到16，突变前后波形对称轴由Y＝1变为Y＝16；本发明在突变前后幅值变化和对称轴偏移明显，特征函数曲线突变迅速，图3(b)模拟信号出现频率变化时，本发明特征函数在突变点的幅值由1到15.6，幅值变化明显，幅值大小增长迅速。
图4为本发明ε比值对振幅或频率变化敏感性分析图；图4(a)所示：本发明ε变化曲线在波形振幅变化后3个采样内即达到峰值，曲线变化明显，；如图4(b)当频率变化时，人工能准确拾取出本发明频率改变位置，曲线变化明显。
图5为本发明实例1拾取过程S-L值变化图；
图6人工拾取微震信号P波定位效果左视图；
图7本发明自动拾取微震信号P波定位效果左视图。
具体实施方式
为了使本发明的技术手段、创作特征、工作流程、使用方法达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。本发明的保护范围不受以下实例的限制。
锦屏地下实验室垂直岩石覆盖达2400米，是目前世界岩石覆盖最深的实验室，工程区的最大主应力可达63MPa，属典型的高应力区，其中7#、8#实验室工程施工时已遭遇到强～极强岩爆，给现场人员和设备安全造成了严重威胁和损失。现场微震监测过程中发现，由于传感器距离岩爆风险较高的7-8#实验室施工段距离较远，大概为400～500m，微震信号在传播过程中衰减严重，且现场监测背景噪音复杂，监测到的微震信号大多数信噪比低，能量幅值小，P波初至不明显，致使微震信号不易识别，P波不易拾取。利用该工程微震样本信号，将本发明拾取和人工拾取结果对比分析，利用粒子群算法优化出适于该工程的本发明最佳参数组：C3＝0.25，C4＝0.08，C5＝2，上限M1＝300，下限M2＝15，T＝0.01。本实例以分析该工程某时窗内实时监测采样点数据为例加以说明。
实例1：
一种工程尺度下微震信号及p波初至自动识别算法，流程如图1所示，其步骤如下：
(1)、读取振动传感器实时监测到的设定时窗内波形的采样点数据，建立平面直角坐标系X0Y，其中X轴对应波形的采样点编号，从平面直角坐标系的原点到X轴正轴方向编号从1开始，按照编号从小到大排列采样点，采样点编号为正整数，Y轴对应波形的幅值，在平面直角坐标系中对时窗内波形的采样点进行对称校准，进入步骤(2)。
波形对称校准：选取时窗中设定个数(在本例中为1000个)连续的采样点作为校正样本，
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内(在本例中为小于0.95或大于1.05)，则选取的时窗内波形的采样点数据需对称校准，计算校正样本幅值的平均值，将步骤(1)中设定时窗内波形的采样点幅值与校正样本幅值的平均值相加，完成对称校准，这样有利于准确统计波形中的零交叉点个数M，减小算法拾取误差。
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值在设定比值范围内(在本例中为大于等于0.95小于等于1.05)，则步骤(1)中设定时窗内波形的采样点数据不需对称校准。
(2)计算对称校准后平面直角坐标系中波形的采样点个数N。设置零交叉点个数上限M1＝300，零交叉点个数下限M2＝15；设置微震信号时间判断阀值T＝0.01，单位为秒。微震信号持续时间短，设置零交叉点个数上限M1和零交叉点个数下限M2，能有效避免持续时间较短的电气尖峰信号干扰和对大部分持续时间较长的非微震信号误判。
(3)、利用对称校准后平面直角坐标系中的波形的采样点数据计算本发明的第A个采样点的加权因子K(A)：
和第A个采样点的特征函数CF(A)：
其中A为对称校准后波形中采样点个数的编号，y(A)为对称校准后波形的第A个采样点的幅值，为y(A)一阶差分。
计算短时平均值STA(A)和长时平均值LTA(A)：
STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]
LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]
其中STA(A)为第A个采样点的短时平均值；C3为短时平均系数，实例1中C3取值为0.25；LTA(A)为第A个采样点长时平均值；C4为长时平均系数，实例1中C4取值为0.08。
计算r(A)＝C5×LTA(A)，r(A)为第A个采样点的动态阀值，C5为触发阀值，实例1中C5值为取2。
(4)、以平面直角坐标系中编号最小的采样点为当前采样点，当前采样点的编号为i。
(5)、设置算法参数变量M、S、L、t，初值均设置为0，其中M为零交叉点个数；S为计数器，起计数作用；L的计算公式为L＝3+M/3；t为拾取出的微震信号的持续时间。
(6)、比较当前采样点的短时平均STA(i)值和动态阀值r(i)的大小。
(7)、若当前采样点的短时平均值STA(i)小于当前采样点的动态阀值r(i)，则当前采样点不是P波初至点，此时将当前采样点的下一个采样点作为当前采样点，返回到步骤(6)；若当前采样点的短时平均值STA(i)大于等于当前点的动态阀值r(i)，则当前采样点可能为P波初至点，将当前采样点的编号i赋值给k，即为k＝i，计算当前采样点的幅值大小P1，为找出当前采样点后第一个峰值，将当前采样点的幅值大小P1赋值给缓存最大幅值P0暂存为最大峰值，然后进行步骤(8)。
(8)、计算对称校准后坐标系波形中采样点的编号为k的下一个采样点k+1的幅值P2。
(9)、如果对称校准后波形中采样点的编号为k的下一个采样点k+1不是零交叉点，比较缓存最大幅值P0、采样点k+1的幅值P2的大小，将较大值赋值给缓存最大幅值P0，将k+1赋值给k，当k的值小于波形采样点个数N时，返回步骤(8)；当k的值大于等于波形采样点个数N时，则将当前采样点的下一个采样点作为当前采样点，即i＝i+1，返回步骤(5)。
如果对称校准后坐标系波形中采样点的编号为k的下一个采样点k+1是零交叉点，则零交叉点个数M值加1，即M＝M+1，将k+1赋值给k，然后进行步骤(10)。
(10)、判断若零交叉点个数M值大于零交叉点个数上限M1值，则当前采样点不为P波初至点，则将当前采样点的下一个采样点作为当前采样点，即i＝i+1，返回步骤(5)。
若零交叉点个数M值小于等于零交叉点个数上限M1值，则计算当前采样点的判定参数δ＝LTA(i)×M，然后进行步骤(11)。
(11)、比较编号为k的采样点的短时平均值STA(k)与当前采样点的判定参数δ的大小：
若编号为k的采样点的短时平均值STA(k)小于等于当前采样点的判定参数δ，计数器S值归零，将P2赋值给P0，此时将k加1，即k＝k+1，返回步骤(8)；
若编号为k的采样点短时平均值STA(k)大于当前采样点的判定参数δ，计数器S值加1，即S＝S+1，计算L＝3+M/3，进入步骤(12)。
(12)、比较计数器S值与L的大小：若计数器S值小于等于L，则将P2赋值给P0，此时将k加1，即k＝k+1，返回步骤(8)；若计数器S值大于L，计算对称校准后坐标系波形中从当前采样点到编号为k的采样点的这段波形的持续时间t，t＝(k-i)/f，其中f为采样频率。
(13)、当t小于微震信号时间判断阀值T，T为0.01，单位为秒，或零交叉点个数M小于零交叉点个数下限M2时，M2为15，则判断当前采样点i与编号为k的采样点之间的采样点不是微震信号，需对此时的采样点k后的信号继续拾取，此时将对称校准后波形坐标系中编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤(5)；当t大于等于微震信号时间阀值T且零交叉点个数M大于等于零交叉点个数下限M2时，则判断当前采样点i与编号为k的采样点之间的采样点为微震信号，将i赋值给T1，k赋值给T2，其中编号为T1的采样点为微震信号的P波初至点，编号为T2的采样点为信号结束点，进入步骤(14)。
(14)、继续对下一个信号自动分析，对称校准后波形坐标系中编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤(5)，直至将设定时窗内波形的采样点数据分析完。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (3)
1.一种工程尺度下微震信号及p波初至自动识别算法，其特征在于，包括以下步骤：
步骤1、读取振动传感器实时监测到的设定时窗内波形的采样点数据，建立平面直角坐标系，其中X轴对应波形的采样点编号，Y轴对应波形的幅值，在平面直角坐标系中对时窗内波形的采样点进行对称校准；
步骤2、计算采样点个数N，设置零交叉点个数上限M1，零交叉点个数下限M2；设置微震信号时间判断阀值T；
步骤3、设定第A个采样点的加权因子K(A)和特征函数CF(A)；
设定第A个采样点短时平均值STA(A)和长时平均值LTA(A)；
设定第A个采样点的动态阈值r(A)；
其中A为采样点个数的编号；
步骤4、以平面直角坐标系中编号最小的采样点为当前采样点，当前采样点的编号为i；
步骤5、设置算法参数变量M、S、L、t，M、S、L、t的初值均设置为0，其中M为零交叉点个数；S为计数器；L的计算公式为L＝3+M/3；t为拾取出的微震信号的持续时间；
步骤6、比较当前采样点的短时平均值STA(i)和动态阈值r(i)的大小；
步骤7、若当前采样点的短时平均值STA(i)小于当前采样点的动态阀值r(i)，则当前采样点不是P波初至点，此时将当前采样点的下一个采样点作为当前采样点，返回到步骤(6)；若当前采样点的短时平均值STA(i)大于等于当前点的动态阀值r(i)，则将当前采样点的编号i赋值给k，计算当前采样点的幅值P1，将当前采样点的幅值大小P1赋值给缓存最大幅值P0，然后进行步骤(8)；
步骤8、计算编号为k的采样点的下一个采样点k+1的幅值P2；
步骤9、如果采样点的编号为k的采样点的下一个采样点k+1不是零交叉点，将缓存最大幅值P0和采样点k+1的幅值P2中较大值赋值给缓存最大幅值P0，将k+1赋值给k，当k的值小于波形采样点个数N时，返回步骤8；当k的值大于等于波形采样点个数N时，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；
如果采样点的编号为k的下一个采样点k+1是零交叉点，则零交叉点个数M值加1，将k+1赋值给k，然后进行步骤10；
步骤10、若零交叉点个数M值大于零交叉点个数上限M1值，则当前采样点不为P波初至点，则将当前采样点的下一个采样点作为当前采样点，返回步骤5；
若零交叉点个数M值小于等于零交叉点个数上限M1值，则计算当前采样点的判定参数δ＝LTA(i)×M，然后进行步骤11；
步骤11、比较编号为k的采样点的短时平均值STA(k)与当前采样点的判定参数δ的大小：
若编号为k的采样点的短时平均值STA(k)小于等于当前采样点的判定参数δ，计数器S值归零，将P2赋值给P0，此时将k加1，返回步骤8；
若编号为k的采样点的短时平均值STA(k)大于当前采样点的判定参数δ，计数器S值加1，计算L＝3+M/3，进入步骤12；
步骤12、比较计数器S值与L的大小：
若计数器S值小于等于L，则将P2赋值给P0，将k加1，返回步骤8；
若计数器S值大于L，计算从当前采样点i到编号为k的采样点的这段波形的持续时间t，t＝(k-i)/f，其中f为采样频率；
步骤13、当t小于微震信号时间判断阀值T或零交叉点个数M小于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点不是微震信号，将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5；当t大于等于时间阀值T且零交叉点个数M大于等于零交叉点个数下限M2时，则当前采样点i与编号为k的采样点之间的采样点为微震信号，将当前采样点的编号i赋值给T1，将k赋值给T2，其中编号为T1的采样点为微震信号的P波初至点，编号为T2的采样点为信号结束点，进入步骤14；
步骤14、将编号为k的采样点的下一个采样点作为当前采样点，清零M、P0、P1、P2、S、L的取值，进行步骤5，直至将设定时窗内波形的采样点数据分析完。
2.根据权利要求1所述的一种工程尺度下微震信号及p波初至自动识别算法，其特征在于，所述的对称校准包括以下步骤：
选取时窗中设定个数的连续的采样点作为校正样本，
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据需对称校准，计算校正样本幅值的平均值，将设定时窗内波形的采样点幅值与校正样本幅值的平均值相加，完成对称校准；
若校正样本中幅值为正的采样点个数与幅值为负的采样点个数的比值不在设定比值范围内，则设定时窗内波形的采样点数据不需对称校准。
3.根据权利要求1所述的一种工程尺度下微震信号及p波初至自动识别算法，其特征在于，
所述的第A个采样点的加权因子K(A)基于以下公式：
<mrow>
<msub>
<mi>K</mi>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mo>)</mo>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>A</mi>
</munderover>
<msup>
<msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>A</mi>
</munderover>
<msup>
<msub>
<mover>
<mi>y</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
所述的第A个采样点的特征函数CF(A)基于以下公式：
<mrow>
<msub>
<mi>CF</mi>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mo>)</mo>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<msup>
<msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<msub>
<mover>
<mi>y</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>A</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
<mo>&times;</mo>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>A</mi>
</munderover>
<msup>
<msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>A</mi>
</munderover>
<msup>
<msub>
<mover>
<mi>y</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
其中，A为对称校准后波形中采样点个数的编号，y(A)为对称校准后波形的第A个采样点的幅值，为y(A)一阶差分；
所述的第A个采样点的短时平均值STA(A)基于以下公式：
STA(A)＝STA(A-1)+C3×[CF(A)-STA(A-1)]
所述的第A个采样点的长时平均值LTA(A)基于以下公式：
LTA(A)＝LTA(A-1)+C4×[CF(A)-LTA(A-1)]
所述的第A个采样点的动态阈值r(A)基于以下公式：
r(A)＝C5×LTA(A)
其中，C3为短时平均系数，C4为长时平均系数，C5为触发阀值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610090278.1A CN105527650B (zh) | 2016-02-17 | 2016-02-17 | 一种工程尺度下微震信号及p波初至自动识别算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610090278.1A CN105527650B (zh) | 2016-02-17 | 2016-02-17 | 一种工程尺度下微震信号及p波初至自动识别算法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105527650A CN105527650A (zh) | 2016-04-27 |
| CN105527650B true CN105527650B (zh) | 2017-10-17 |
Family
ID=55769986
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610090278.1A Expired - Fee Related CN105527650B (zh) | 2016-02-17 | 2016-02-17 | 一种工程尺度下微震信号及p波初至自动识别算法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105527650B (zh) |
Families Citing this family (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106199703B (zh) * | 2016-08-26 | 2018-11-16 | 中国矿业大学 | 一种微震震源自动定位及可靠性综合评价方法 |
| CN107479094B (zh) * | 2017-09-18 | 2018-11-30 | 辽宁工程技术大学 | 一种实现地震预警的方法 |
| CN107992802B (zh) * | 2017-11-10 | 2021-12-21 | 桂林电子科技大学 | 一种基于nmf的微震弱信号识别方法 |
| CN108319915B (zh) * | 2018-01-31 | 2020-06-19 | 中国科学院武汉岩土力学研究所 | 一种岩爆信号阈值动态调整的多时窗简化形式识别方法 |
| CN108254781B (zh) * | 2018-01-31 | 2019-07-16 | 中国科学院武汉岩土力学研究所 | 一种岩爆信号阈值动态调整的多时窗完整形式识别方法 |
| CN108376245B (zh) * | 2018-02-02 | 2022-02-11 | 广西师范大学 | 基于ud通道的时空序列图像震源识别方法 |
| CN108426949A (zh) * | 2018-02-14 | 2018-08-21 | 国家***第二海洋研究所 | 一种海底沉积物声学原位数据初至识别拾取方法 |
| CN110146920B (zh) * | 2019-06-26 | 2020-11-10 | 广东石油化工学院 | 基于幅值相对变化的微震事件检测方法和*** |
| CN110297271B (zh) * | 2019-06-26 | 2020-09-11 | 中国矿业大学 | 一种用于矿震报警的单分量探头p波初至到时修正方法 |
| CN112526602B (zh) * | 2020-11-16 | 2023-10-20 | 重庆大学 | 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 |
| CN112364296B (zh) * | 2020-11-17 | 2023-08-04 | 东北大学 | 一种基于深度学习的p波到时自动拾取方法 |
| CN113239620B (zh) * | 2021-05-11 | 2023-03-28 | 中国电建集团华东勘测设计研究院有限公司 | 一种基于gpu加速的用于岩土材料本构模型参数识别的改进粒子群方法 |
| CN114333193B (zh) * | 2022-01-11 | 2024-04-16 | 广西北投交通养护科技集团有限公司 | 道路交通便携式应急声光警报及施工边坡监测装置 |
| CN115146678B (zh) * | 2022-07-04 | 2023-05-09 | 长江水利委员会长江科学院 | 一种***振动信号的p波振相初至识别方法和*** |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103995290A (zh) * | 2014-06-03 | 2014-08-20 | 山东科技大学 | 一种高精度微震p波震相初至自动拾取方法 |
| CN104062677A (zh) * | 2014-07-03 | 2014-09-24 | 中国科学院武汉岩土力学研究所 | 一种多功能综合集成高精度智能微震监测*** |
| CA2515345C (en) * | 2003-02-08 | 2014-12-09 | Robert Hughes Jones | Estimating the time of arrival of a seismic wave |
| CN104777512A (zh) * | 2014-12-31 | 2015-07-15 | 安徽万泰地球物理技术有限公司 | 一种地震波s波自动拾取的算法 |
| CN104914468A (zh) * | 2015-06-09 | 2015-09-16 | 中南大学 | 一种矿山微震信号p波初至时刻联合拾取方法 |
| CN105223614A (zh) * | 2015-09-23 | 2016-01-06 | 中南大学 | 一种基于dwt_sta/lta的含噪信号p波初至峰度拾取方法 |
-
2016
- 2016-02-17 CN CN201610090278.1A patent/CN105527650B/zh not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2515345C (en) * | 2003-02-08 | 2014-12-09 | Robert Hughes Jones | Estimating the time of arrival of a seismic wave |
| CN103995290A (zh) * | 2014-06-03 | 2014-08-20 | 山东科技大学 | 一种高精度微震p波震相初至自动拾取方法 |
| CN104062677A (zh) * | 2014-07-03 | 2014-09-24 | 中国科学院武汉岩土力学研究所 | 一种多功能综合集成高精度智能微震监测*** |
| CN104777512A (zh) * | 2014-12-31 | 2015-07-15 | 安徽万泰地球物理技术有限公司 | 一种地震波s波自动拾取的算法 |
| CN104914468A (zh) * | 2015-06-09 | 2015-09-16 | 中南大学 | 一种矿山微震信号p波初至时刻联合拾取方法 |
| CN105223614A (zh) * | 2015-09-23 | 2016-01-06 | 中南大学 | 一种基于dwt_sta/lta的含噪信号p波初至峰度拾取方法 |
Non-Patent Citations (2)
| Title |
|---|
| STA/LTA算法拾取微地震事件P波到时对比研究;吴治涛 等;《地球物理学进展》;20101031;第25卷(第5期);第1577-1582页 * |
| 微震信号自动检测的STA／LTA算法及其改进分析;刘晗 等;《地球物理学进展》;20140429;第29卷(第4期);第1708-1714页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105527650A (zh) | 2016-04-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105527650B (zh) | 一种工程尺度下微震信号及p波初至自动识别算法 | |
| CN110032975B (zh) | 一种地震震相的拾取方法 | |
| CN108919353B (zh) | 一种微震波形初至到时的自动分级拾取与优选方法 | |
| CN105261136B (zh) | 一种光纤监测报警***中屏蔽天气干扰的方法及装置 | |
| CN106382981B (zh) | 一种单站次声波信号识别提取方法 | |
| CN103776480A (zh) | 基于多次移动平均的微小故障检测方法和装置 | |
| CN103994817A (zh) | 一种基于长距离光纤多发事件的振动振源的识别方法 | |
| CN111126434B (zh) | 基于随机森林的微震初至波到时自动拾取方法及*** | |
| EP2132720A1 (en) | Method and apparatus for monitoring a structure | |
| US20240078413A1 (en) | Massive data-driven method for automatically locating mine microseismic source | |
| CN103994815A (zh) | 一种识别光纤振源的方法 | |
| CN106646610B (zh) | 一种利用偏振约束aic算法自动拾取微震初至的算法 | |
| CN104902509A (zh) | 基于top-k(σ)算法的异常数据检测方法 | |
| CN110398775B (zh) | 隧道突涌水灾害微震事件信号波动初至拾取方法及*** | |
| CN114330120B (zh) | 一种基于深度神经网络预测24小时pm2.5浓度的方法 | |
| CN104977602B (zh) | 一种地震数据采集施工的控制方法及装置 | |
| CN116699724B (zh) | 一种时间域激发极化数据质量评价方法、体系及*** | |
| CN105222885A (zh) | 一种光纤振动检测方法及装置 | |
| CN103994816A (zh) | 一种基于光纤多发事件的识别方法 | |
| CN113065388A (zh) | 一种实时土体类别识别方法、***及一种挖掘机 | |
| CN109212601B (zh) | 一种地震数据异常测点检测方法 | |
| CN104819382B (zh) | 一种用于光纤预警***的自适应恒虚警率振源检测方法 | |
| CN106771928A (zh) | 一种局部放电脉冲初至时刻在线拾取方法 | |
| CN110632191A (zh) | 一种基于决策树算法的变压器色谱峰定性方法和*** | |
| CN111681124B (zh) | 一种深部砂岩型铀矿化信息三维氡异常识别方法及*** |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into 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: 20171017 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |