CN105527650A - Automatic identification algorithm for microseismic signal and p wave first arrival at engineering scale - Google Patents
Automatic identification algorithm for microseismic signal and p wave first arrival at engineering scale Download PDFInfo
- Publication number
- CN105527650A CN105527650A CN201610090278.1A CN201610090278A CN105527650A CN 105527650 A CN105527650 A CN 105527650A CN 201610090278 A CN201610090278 A CN 201610090278A CN 105527650 A CN105527650 A CN 105527650A
- Authority
- CN
- China
- Prior art keywords
- point
- sampled point
- current sampling
- value
- 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.)
- Granted
Links
- 238000001615 p wave Methods 0.000 title abstract 4
- 238000005070 sampling Methods 0.000 claims abstract description 103
- 230000006870 function Effects 0.000 claims abstract description 12
- 238000007405 data analysis Methods 0.000 claims description 3
- 230000008859 change Effects 0.000 description 14
- 238000000034 method Methods 0.000 description 14
- 238000012544 monitoring process Methods 0.000 description 8
- 239000011435 rock Substances 0.000 description 8
- 230000008569 process Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000010276 construction Methods 0.000 description 3
- 238000010206 sensitivity analysis Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 206010034719 Personality change Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000006854 communication Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 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
The invention discloses an automatic identification algorithm for microseismic signals and p wave first arrival at an engineering scale. The automatic identification algorithm comprises the steps of: reading data of sampling points in a set time window; calculating feature functions as well as long and short time average values of the sampling points; carrying out automatic pickup judgment on the microseismic signals and the p wave first arrival; clearing buffer memory of parameter values to zero; and reading data in a time window at a next period. According to the automatic identification algorithm, the identification algorithm in the seismic field is improved aiming at microseismic signal characteristics at the engineering scale, the automatic identification algorithm can automatically pick up microseismic signals especially weak signals with low signal-to-noise ratio at the engineering scale accurately, increases the accuracy of automatically picking up p wave first arrival, and further enhances early warning timeliness and accuracy of geological disasters such as rockburst, mine earthquake and collapse.
Description
Technical field
The present invention relates to On Microseismic Monitoring Technique field.Be specifically related to microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick, the method can be widely used in Mineral Engineering, Hydraulic and Hydro-Power Engineering, petroleum engineering, Geotechnical Engineering and underground works.
Background technology
Microquake sources location is the important component part of micro seismic monitoring and disaster alarm, and the identification of microseismic signals and P ripple first break pickup are the keys of microquake sources location.Existing microseismic signals recognition methods major part derives from seismic field, and method is a lot, mainly contains: according to the STA/LTA algorithm at time domain energy and energy variation construction feature versus time territory; Seismic waveform data statistically qualitative difference before and after arriving according to seismic signal, as AIC algorithm; In addition neural network, high order statistics and the Wavelet Transform etc. based on wavelet theory is also had.STA/LTA and innovatory algorithm thereof, because algorithm is simple, counting yield is high, be suitable for real-time process, extensively used in seismic field.The process that these methods are used for seismic signal has good applicability, but the rock burst signal produced under engineering yardstick is different from earthquake signal:
1) seismic signal frequency is generally less than tens hertz, and the frequency of microseismic signals is generally tens to hundreds of hertz, have up to a few KHz; 2) the seismic signal duration is longer, is generally greater than 1 second, and the microseismic signals duration is shorter, is generally no more than 0.1 second; 3) microseism is compared with earthquake, and general earthquake magnitude is little, and less microseismic event is not easily picked; 4) under engineering-environment because the factors such as mechanical execution, electrical Interference, vehicle traveling cause background noise complicated, make most of microseismic signals signal to noise ratio (S/N ratio)s of monitoring lower.
Therefore, by the Signal analysis of seismic field and P ripple, then pick-up method is used for engineering yardstick microseismic signals is inappropriate, and main manifestations is: 1) be difficult to the feeble signal identifying low signal-to-noise ratio; 2) P ripple first arrival automatic Picking error is larger.
For this reason, for above-mentioned deficiency, invent a kind of can the automatic knowledge method for distinguishing of effectively microseismic signals and the first arrival of p ripple thereof under recognitive engineering yardstick, the pickup accuracy rate of microseismic signals and the first arrival of P ripple thereof under raising engineering yardstick, and then improve promptness and the accuracy of the geo-hazard early-warnings such as rock burst, ore deposit shake, landslide, be very important.
Summary of the invention
The object of the invention is to overcome prior art Problems existing, provide microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick, the pickup accuracy rate of microseismic signals and the first arrival of P ripple thereof under raising engineering yardstick, and then improve promptness and the accuracy of the geo-hazard early-warnings such as rock burst, ore deposit shake, landslide.
Microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick, comprise the following steps:
The sampling number certificate of establishing waveform in timing window that step 1, reading vibration transducer real-time monitor, set up plane right-angle coordinate, the wherein sampled point numbering of the corresponding waveform of X-axis, the amplitude of the corresponding waveform of Y-axis, time in plane right-angle coordinate pair, in window, the sampled point of waveform carries out symmetry and calibrates;
Step 2, calculating sampling point number N, arranges zero cross point number upper limit M1, zero cross point number lower limit M2; Microseismic signals time decision threshold T is set;
Step 3, set the weighting factor K of A sampled point
(A)with fundamental function CF
(A);
Set A sampled point short time average STA
(A)with mean value LTA time long
(A);
Set the dynamic threshold r (A) of A sampled point;
Wherein A is the numbering of sampled point number;
Step 4, to number minimum sampled point in plane right-angle coordinate for current sampling point, current sampling point be numbered i;
Step 5, arrange algorithm parameter variable M, S, L, t, the initial value of M, S, L, t is all set to 0, and wherein M is zero cross point number; S is counter; The computing formula of L is L=3+M/3; T is the duration of the microseismic signals of picking up out;
Step 6, compare the short time average STA of current sampling point
(i)with the size of dynamic threshold r (i);
If the short time average STA of step 7 current sampling point
(i)be less than dynamic thresholding r (i) of current sampling point, then current sampling point is not P ripple Onset point, now using the next sampled point of current sampling point as current sampling point, turn back to step (6); If the short time average STA of current sampling point
(i)be more than or equal to dynamic thresholding r (i) of current point, then by the numbering i assignment of current sampling point to k, calculate the amplitude P1 of current sampling point, by the amplitude size P1 assignment of current sampling point to buffer memory maximum amplitude P0, then carry out step (8);
Step 8, calculating are numbered the amplitude P2 of the next sampled point k+1 of the sampled point of k;
If the next sampled point k+1 being numbered the sampled point of k of step 9 sampled point is not zero cross point, by higher value assignment in the amplitude P2 of buffer memory maximum amplitude P0 and sampled point k+1 to buffer memory maximum amplitude P0, by k+1 assignment to k, when the value of k is less than waveform sampling point number N, return step 8; When the value of k is more than or equal to waveform sampling point number N, then using the next sampled point of current sampling point as current sampling point, return step 5;
If the next sampled point k+1 being numbered k of sampled point is zero cross point, then zero cross point number M value adds 1, by k+1 assignment to k, then carry out step 10;
If step 10 zero cross point number M value is greater than zero cross point number upper limit M1 value, then current sampling point is not P ripple Onset point, then using the next sampled point of current sampling point as current sampling point, return step 5;
If zero cross point number M value is less than or equal to zero cross point number upper limit M1 value, then calculates critical parameter δ=LTA (the i) × M of current sampling point, then carry out step 11;
Step 11, compare the short time average STA of the sampled point being numbered k
(k)with the size of the critical parameter δ of current sampling point:
If be numbered the short time average STA of the sampled point of k
(k)be less than or equal to the critical parameter δ of current sampling point, counter S value makes zero, and by P2 assignment to P0, now k is added 1, returns step 8;
If be numbered the short time average STA of the sampled point of k
(k)be greater than the critical parameter δ of current sampling point, counter S value adds 1, calculates L=3+M/3, enters step 12;
Step 12, compare the size of counter S value and L:
If counter S value is less than or equal to L, then by P2 assignment to P0, k is added 1, returns step 8;
If counter S value is greater than L, calculate from current sampling point i to the duration t of this section of waveform of sampled point being numbered k, t=(k-i)/f, wherein f is sample frequency;
Step 13, when t be less than microseismic signals time decision threshold T or zero cross point number M be less than zero cross point number lower limit M2 time, sampled point then between current sampling point i and the sampled point being numbered k is not microseismic signals, to the next sampled point of the sampled point of k be numbered as current sampling point, reset the value of M, P0, P1, P2, S, L, carry out step 5; When t is more than or equal to time threshold values T and zero cross point number M is more than or equal to zero cross point number lower limit M2, sampled point then between current sampling point i and the sampled point being numbered k is microseismic signals, by the numbering i assignment of current sampling point to T1, by k assignment to T2, the sampled point being wherein numbered T1 is the P ripple Onset point of microseismic signals, the sampled point being numbered T2 is signal ended point, enters step 14;
Step 14, will the next sampled point of the sampled point of k be numbered as current sampling point, reset the value of M, P0, P1, P2, S, L, carry out step 5, until will the sampled point data analysis of waveform in timing window be established complete.
Symmetry calibration as above comprises the following steps:
The continuous print sampled point of number is set in window as calibration samples when choosing,
If the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number is not setting in ratio range in calibration samples, then establish the sampling number of waveform in timing window according to needing symmetrical calibration, calculate the mean value of calibration samples amplitude, to the sampled point amplitude of waveform in timing window be established to be added with the mean value of calibration samples amplitude, complete symmetrical calibration;
If the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number is not setting in ratio range in calibration samples, then establish the sampling number of waveform in timing window according to not needing symmetrical calibration.
The weighting factor K of A sampled point as above
(A)based on following formula:
The fundamental function CF of A described sampled point
(A)based on following formula:
Wherein, A is the numbering of sampled point number in waveform after symmetrical calibration, y
(A)for the amplitude of A sampled point of waveform after symmetry calibration,
for y
(A)first order difference;
The short time average STA of A described sampled point
(A)based on following formula:
STA
(A)＝STA
(A-1)+C
3×[CF
(A)-STA
(A-1)]
Mean value LTA during A described sampled point long
(A)based on following formula:
LTA
(A)＝LTA
(A-1)+C
4×[CF
(A)-LTA
(A-1)]
The dynamic threshold r (A) of A described sampled point is based on following formula:
r(A)＝C
5×LTA
(A)
Wherein, C
3for short-time average coefficient, C
4for mean coefficient time long, C
5for trigger threshold.
When P ripple arrives, will cause the waveform character changes such as microseismic signals amplitude or frequency, if the signal to noise ratio (S/N ratio) of signal is lower, then the first arrival of P ripple is more not obvious.Thus improve picking algorithm susceptibility that such change, the P ripple first break pickup ability of microseismic signals under raising engineering yardstick can be conducive to.And weighting factor K in algorithm, fundamental function CF, ε ratio all can the changing features of response waveform, directly accuracy rate is picked up in impact:
Analyze weighting factor K of the present invention, fundamental function CF, dynamically length time mean ratio ε to the susceptibility of the changing features such as signal amplitude or frequency waveform, respectively as Fig. 2,3,4.Comparing result is analyzed the present invention's three kinds of variablees and is changed obviously at frequency or amplitude variations place numerical value, reaches the rapid of peak value, therefore in theory the present invention to the pickup successful of the low signal of signal to noise ratio (S/N ratio) and the first arrival of P ripple thereof.
From the data of silk screen diversion tunnel Real-Time Monitoring, random selecting 555 microseismic signals are sample, define two pickup result judgment criterion: 1) signal successfully picks up: sample signal can be automatically identified and not produce leaky wave and erroneous judgement.Manually automatic Picking signal section is out further processed being conducive to like this.It is more that signal successfully picks up number, and picking up signal rate is higher.2) P ripple first break pickup is accurate: automatic algorithms is picked up P ripple first arrival result and manually pick up compared with result, if picking errors is not more than 3 sampled points, then thinks that this algorithm automatic Picking result is accurate.Picking up signal rate of the present invention is 94.23%, P ripple pickup accuracy rate is 73.51%, and pickup rate meets the current positioning error in microseism field completely.
Analyze silk screen on July 15 in 15 days to 2015 June profound subterranean laboratory micro seismic monitoring data in 2015,524 microseismic signals altogether, 109 events, by manually picking up P ripple, the present invention picks up result and is used for location, respectively as shown in Figure 6,7, wherein 7# laboratory is rock burst high risk zone to result, can find out that the concentration class of positioning result microseismic event of the present invention is higher, microseismic event space distribution rule is close with artificial result of picking up, and conforms to engineering practice.By very low for signal to noise ratio (S/N ratio) in artificial pick process, that P ripple is not obvious, positioning result is obviously abnormal 13 events, as shown in the triangle in Fig. 6, positioning analysis is picked up: 1) the P ripple result of the present invention's pickup all can successfully be located, as shown in the triangle in Fig. 7 with the present invention; 2) these 13 event pickup locating effects of the present invention improve obviously relative to artificial pickup, and the overwhelming majority is in the scope that positioning error allows.
Beneficial effect:
The present invention is directed to microseismic signals feature under engineering yardstick, improve the recognizer of seismic field, the present invention can the microseismic signals weak signal that particularly to-noise ratio is lower accurately under automatic Picking engineering yardstick, improve the accuracy rate of automatic Picking P ripple first arrival simultaneously, and then improve promptness and the accuracy of the geo-hazard early-warnings such as rock burst, ore deposit shake, landslide.
Accompanying drawing explanation
Fig. 1 is process flow diagram of the present invention;
Fig. 2 is that weighting factor K of the present invention is to amplitude or frequency change sensitivity analysis figure; If Fig. 2 (a) is when amplitude changes, weighting factor amplitude of the present invention about increases 2 times; If Fig. 2 (b) is when frequency changes, K value amplitude of variation of the present invention is comparatively obvious.
Fig. 3 is that fundamental function CF of the present invention is to amplitude or frequency change sensitivity analysis figure; There is the microseismic signals that amplitude is less in Fig. 3 (a) simulation background noise, namely amplitude by 1 sudden change to 2 time, fundamental function of the present invention is in the amplitude of catastrophe point by 1 to 16, and before and after sudden change, waveform symmetry axle becomes Y=16 from Y=1; The present invention amplitude change and axis of symmetry before and after sudden change offset obviously, the sudden change of fundamental function curve rapidly, during the change of Fig. 3 (b) simulating signal frequency of occurrences, fundamental function of the present invention in the amplitude of catastrophe point by 1 to 15.6, amplitude change is obvious, amplitude size rapid development.
Fig. 4 is that ε ratio of the present invention is to amplitude or frequency change sensitivity analysis figure; Shown in Fig. 4 (a): namely ε change curve of the present invention reaches peak value in rear 3 samplings of amplitude of wave form change, and curvilinear motion is obvious; If Fig. 4 (b) is when frequency change, manually accurately can pick up out frequency shift position of the present invention, curvilinear motion is obvious.
Fig. 5 is example 1 pick process S-L value variation diagram of the present invention;
Fig. 6 manually picks up microseismic signals P ripple locating effect left view;
Fig. 7 automatic Picking microseismic signals P of the present invention ripple locating effect left view.
Embodiment
Reaching object to make technological means of the present invention, creation characteristic, workflow, using method and effect is easy to understand, below in conjunction with specific embodiment, setting forth the present invention further.Protection scope of the present invention is not by the restriction of following instance.
The vertical rock of silk screen subterranean laboratory covers and reaches 2400 meters, that current world rock covers the darkest laboratory, the major principal stress of Project Areas can reach 63MPa, belong to the typical region of high stress, wherein suffer from strong ~ strong rock burst on active during the engineering construction of 7#, 8# laboratory, cause serious threat and loss to field staff and device security.Find in on-the-spot micro seismic monitoring process, the 7-8# laboratory construction section higher due to sensor distance rockburst risk is distant, be probably 400 ~ 500m, microseismic signals decays seriously in communication process, and field monitoring background noise is complicated, and the most of signal to noise ratio (S/N ratio) of the microseismic signals monitored is low, energy magnitude is little, the first arrival of P ripple is not obvious, and cause microseismic signals not easy to identify, P ripple not easily picks up.Utilize this project microseism sample signal, by the present invention's pickup and the analysis of artificial pickup Comparative result, utilize particle cluster algorithm optimization to be suitable for the optimal parameter group of the present invention of this project: C
3=0.25, C
4=0.08, C
5=2, upper limit M1=300, lower limit M2=15, T=0.01.When this example is to analyze this project, in window, Real-Time Monitoring sampled point data instance is illustrated.
Example 1:
Microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick, as shown in Figure 1, its step is as follows for flow process:
(1) the sampling number certificate of establishing waveform in timing window that vibration transducer real-time monitors, is read, set up plane right-angle coordinate X0Y, the wherein sampled point numbering of the corresponding waveform of X-axis, from the initial point of plane right-angle coordinate to the positive direction of principal axis numbering of X-axis from 1, arrange sampled point from small to large according to numbering, sampled point is numbered positive integer, the amplitude of the corresponding waveform of Y-axis, time in plane right-angle coordinate pair, in window, the sampled point of waveform carries out symmetry calibration, enters step (2).
Waveform symmetry is calibrated: set number (being 1000 in this example) continuous print sampled point when choosing in window as calibration samples,
If in calibration samples, the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number is not setting in ratio range (in this example for be less than 0.95 or be greater than 1.05), then choose time window in the sampling number of waveform according to needing symmetrical calibration, calculate the mean value of calibration samples amplitude, be added establishing the sampled point amplitude of waveform in timing window in step (1) with the mean value of calibration samples amplitude, complete symmetrical calibration, be conducive to the zero cross point number M in accurate statistics waveform like this, reduce algorithm picking errors.
If the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number (is less than or equal to 1.05 for being more than or equal to 0.95) in this example setting in ratio range in calibration samples, then establish the sampling number of waveform in timing window according to not needing symmetrical calibration in step (1).
(2) the sampled point number N of waveform in symmetrical calibration back plane rectangular coordinate system is calculated.Zero cross point number upper limit M1=300 is set, zero cross point number lower limit M2=15; Arrange microseismic signals time decision threshold T=0.01, unit is second.The microseismic signals duration is short, arranges zero cross point number upper limit M1 and zero cross point number lower limit M2, and the electric spiking effectively avoiding the duration shorter disturbs and judges by accident the long-term non-microseismic signals of major part.
(3), utilize the sampling number of the waveform in symmetrical calibration back plane rectangular coordinate system according to the weighting factor K calculating A sampled point of the present invention
(A):
With the fundamental function CF of A sampled point
(A):
Wherein A is the numbering of sampled point number in waveform after symmetrical calibration, y
(A)for the amplitude of A sampled point of waveform after symmetry calibration,
for y
(A)first order difference.
Calculate short time average STA
(A)with mean value LTA time long
(A):
STA
(A)＝STA
(A-1)+C
3×[CF
(A)-STA
(A-1)]
LTA
(A)＝LTA
(A-1)+C
4×[CF
(A)-LTA
(A-1)]
Wherein STA
(A)it is the short time average of A sampled point; C
3for short-time average coefficient, C in example 1
3value is 0.25; LTA
(A)be A sampled point long time mean value; C
4for mean coefficient time long, C in example 1
4value is 0.08.
Calculate r (A)=C
5× LTA
(A), r (A) is the dynamic thresholding of A sampled point, C
5for trigger threshold, C in example 1
5value is for getting 2.
(4), to number minimum sampled point in plane right-angle coordinate for current sampling point, current sampling point be numbered i.
(5), arrange algorithm parameter variable M, S, L, t, initial value is all set to 0, and wherein M is zero cross point number; S is counter, plays count action; The computing formula of L is L=3+M/3; T is the duration of the microseismic signals of picking up out.
(6) the short-time average STA of current sampling point, is compared
(i)the size of value and dynamic thresholding r (i).
(7) if the short time average STA of current sampling point
(i)be less than dynamic thresholding r (i) of current sampling point, then current sampling point is not P ripple Onset point, now using the next sampled point of current sampling point as current sampling point, turn back to step (6); If the short time average STA of current sampling point
(i)be more than or equal to dynamic thresholding r (i) of current point, then current sampling point may be P ripple Onset point, by the numbering i assignment of current sampling point to k, be k=i, calculate the amplitude size P1 of current sampling point, for finding out first peak value after current sampling point, the amplitude size P1 assignment of current sampling point being kept in as peak-peak to buffer memory maximum amplitude P0, then carries out step (8).
(8) the amplitude P2 being numbered the next sampled point k+1 of k of sampled point in the rear coordinate system waveform of symmetrical calibration, is calculated.
(9) if the next sampled point k+1 being numbered k of sampled point is not zero cross point in the rear waveform of symmetrical calibration, the relatively size of the amplitude P2 of buffer memory maximum amplitude P0, sampled point k+1, by higher value assignment to buffer memory maximum amplitude P0, by k+1 assignment to k, when the value of k is less than waveform sampling point number N, return step (8); When the value of k is more than or equal to waveform sampling point number N, then using the next sampled point of current sampling point as current sampling point, i.e. i=i+1, returns step 0.
If in coordinate system waveform, the next sampled point k+1 being numbered k of sampled point is zero cross point after symmetrical calibration, then zero cross point number M value adds 1, i.e. M=M+1, by k+1 assignment to k, then carries out step (10).
(10) if, judge that zero cross point number M value is greater than zero cross point number upper limit M1 value, then current sampling point is not P ripple Onset point, then using the next sampled point of current sampling point as current sampling point, i.e. i=i+1, returns step (5).
If zero cross point number M value is less than or equal to zero cross point number upper limit M1 value, then calculate critical parameter δ=LTA (the i) × M of current sampling point, then carry out step (11).
(11) the short time average STA of the sampled point being numbered k, is compared
(k)with the size of the critical parameter δ of current sampling point:
If be numbered the short time average STA of the sampled point of k
(k)be less than or equal to the critical parameter δ of current sampling point, counter S value makes zero, and by P2 assignment to P0, now k is added 1, i.e. k=k+1, returns step (8);
If be numbered the sampled point short time average STA of k
(k)be greater than the critical parameter δ of current sampling point, counter S value adds 1, i.e. S=S+1, calculates L=3+M/3, enters step (12).
(12), compare the size of counter S value and L: if counter S value is less than or equal to L, then by P2 assignment to P0, now k is added 1, i.e. k=k+1, returns step (8); If counter S value is greater than L, after calculating symmetrical calibration in coordinate system waveform from current sampling point to the duration t of this section of waveform of sampled point being numbered k, t=(k-i)/f, wherein f is sample frequency.
(13), when t is less than microseismic signals time decision threshold T, T is 0.01, unit is second, or zero cross point number M is when being less than zero cross point number lower limit M2, M2 is 15, then judge that the sampled point between current sampling point i and the sampled point being numbered k is not microseismic signals, need to the signal subsequent pick-up after sampled point k now, now be numbered the next sampled point of the sampled point of k as current sampling point using after symmetry calibration in waveform coordinate system, reset the value of M, P0, P1, P2, S, L, carry out step (5); When t is more than or equal to microseismic signals time threshold values T and zero cross point number M is more than or equal to zero cross point number lower limit M2, then judge that the sampled point between current sampling point i and the sampled point being numbered k is microseismic signals, by i assignment to T1, k assignment is to T2, the sampled point being wherein numbered T1 is the P ripple Onset point of microseismic signals, the sampled point being numbered T2 is signal ended point, enters step (14).
(14), continue next signal automatic analysis, the next sampled point of the sampled point of k is numbered as current sampling point in waveform coordinate system after symmetrical calibration, reset the value of M, P0, P1, P2, S, L, carry out step (5), until will the sampled point data analysis of waveform in timing window be established complete.
More than show and describe ultimate principle of the present invention and principal character and advantage of the present invention.The technician of the industry should understand; the present invention is not restricted to the described embodiments; what describe in above-described embodiment and instructions just illustrates principle of the present invention; without departing from the spirit and scope of the present invention; the present invention also has various changes and modifications, and these changes and improvements all fall in the claimed scope of the invention.Application claims protection domain is defined by appending claims and equivalent thereof.
Claims (3)
1. microseismic signals and a p ripple first arrival automatic identification algorithm under engineering yardstick, is characterized in that, comprise the following steps:
The sampling number certificate of establishing waveform in timing window that step 1, reading vibration transducer real-time monitor, set up plane right-angle coordinate, the wherein sampled point numbering of the corresponding waveform of X-axis, the amplitude of the corresponding waveform of Y-axis, time in plane right-angle coordinate pair, in window, the sampled point of waveform carries out symmetry and calibrates;
Step 2, calculating sampling point number N, arranges zero cross point number upper limit M1, zero cross point number lower limit M2; Microseismic signals time decision threshold T is set;
Step 3, set the weighting factor K of A sampled point
(A)with fundamental function CF
(A);
Set A sampled point short time average STA
(A)with mean value LTA time long
(A);
Set the dynamic threshold r (A) of A sampled point;
Wherein A is the numbering of sampled point number;
Step 4, to number minimum sampled point in plane right-angle coordinate for current sampling point, current sampling point be numbered i;
Step 5, arrange algorithm parameter variable M, S, L, t, the initial value of M, S, L, t is all set to 0, and wherein M is zero cross point number; S is counter; The computing formula of L is L=3+M/3; T is the duration of the microseismic signals of picking up out;
Step 6, compare the short time average STA of current sampling point
(i)with the size of dynamic threshold r (i);
If the short time average STA of step 7 current sampling point
(i)be less than dynamic thresholding r (i) of current sampling point, then current sampling point is not P ripple Onset point, now using the next sampled point of current sampling point as current sampling point, turn back to step (6); If the short time average STA of current sampling point
(i)be more than or equal to dynamic thresholding r (i) of current point, then by the numbering i assignment of current sampling point to k, calculate the amplitude P1 of current sampling point, by the amplitude size P1 assignment of current sampling point to buffer memory maximum amplitude P0, then carry out step (8);
Step 8, calculating are numbered the amplitude P2 of the next sampled point k+1 of the sampled point of k;
If the next sampled point k+1 being numbered the sampled point of k of step 9 sampled point is not zero cross point, by higher value assignment in the amplitude P2 of buffer memory maximum amplitude P0 and sampled point k+1 to buffer memory maximum amplitude P0, by k+1 assignment to k, when the value of k is less than waveform sampling point number N, return step 8; When the value of k is more than or equal to waveform sampling point number N, then using the next sampled point of current sampling point as current sampling point, return step 5;
If the next sampled point k+1 being numbered k of sampled point is zero cross point, then zero cross point number M value adds 1, by k+1 assignment to k, then carry out step 10;
If step 10 zero cross point number M value is greater than zero cross point number upper limit M1 value, then current sampling point is not P ripple Onset point, then using the next sampled point of current sampling point as current sampling point, return step 5;
If zero cross point number M value is less than or equal to zero cross point number upper limit M1 value, then calculates critical parameter δ=LTA (the i) × M of current sampling point, then carry out step 11;
Step 11, compare the short time average STA of the sampled point being numbered k
(k)with the size of the critical parameter δ of current sampling point:
If be numbered the short time average STA of the sampled point of k
(k)be less than or equal to the critical parameter δ of current sampling point, counter S value makes zero, and by P2 assignment to P0, now k is added 1, returns step 8;
If be numbered the short time average STA of the sampled point of k
(k)be greater than the critical parameter δ of current sampling point, counter S value adds 1, calculates L=3+M/3, enters step 12;
Step 12, compare the size of counter S value and L:
If counter S value is less than or equal to L, then by P2 assignment to P0, k is added 1, returns step 8;
If counter S value is greater than L, calculate from current sampling point i to the duration t of this section of waveform of sampled point being numbered k, t=(k-i)/f, wherein f is sample frequency;
Step 13, when t be less than microseismic signals time decision threshold T or zero cross point number M be less than zero cross point number lower limit M2 time, sampled point then between current sampling point i and the sampled point being numbered k is not microseismic signals, to the next sampled point of the sampled point of k be numbered as current sampling point, reset the value of M, P0, P1, P2, S, L, carry out step 5; When t is more than or equal to time threshold values T and zero cross point number M is more than or equal to zero cross point number lower limit M2, sampled point then between current sampling point i and the sampled point being numbered k is microseismic signals, by the numbering i assignment of current sampling point to T1, by k assignment to T2, the sampled point being wherein numbered T1 is the P ripple Onset point of microseismic signals, the sampled point being numbered T2 is signal ended point, enters step 14;
Step 14, will the next sampled point of the sampled point of k be numbered as current sampling point, reset the value of M, P0, P1, P2, S, L, carry out step 5, until will the sampled point data analysis of waveform in timing window be established complete.
2. microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick according to claim 1, is characterized in that, described symmetry calibration comprises the following steps:
The continuous print sampled point of number is set in window as calibration samples when choosing,
If the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number is not setting in ratio range in calibration samples, then establish the sampling number of waveform in timing window according to needing symmetrical calibration, calculate the mean value of calibration samples amplitude, to the sampled point amplitude of waveform in timing window be established to be added with the mean value of calibration samples amplitude, complete symmetrical calibration;
If the ratio of amplitude to be positive sampled point number and amplitude be negative sampled point number is not setting in ratio range in calibration samples, then establish the sampling number of waveform in timing window according to not needing symmetrical calibration.
3. microseismic signals and p ripple first arrival automatic identification algorithm under a kind of engineering yardstick according to claim 1, is characterized in that,
The weighting factor K of A described sampled point
(A)based on following formula:
The fundamental function CF of A described sampled point
(A)based on following formula:
Wherein, A is the numbering of sampled point number in waveform after symmetrical calibration, y
(A)for the amplitude of A sampled point of waveform after symmetry calibration,
for y
(A)first order difference;
The short time average STA of A described sampled point
(A)based on following formula:
STA
(A)＝STA
(A-1)+C
3×[CF
(A)-STA
(A-1)]
Mean value LTA during A described sampled point long
(A)based on following formula:
LTA
(A)＝LTA
(A-1)+C
4×[CF
(A)-LTA
(A-1)]
The dynamic threshold r (A) of A described sampled point is based on following formula:
r(A)＝C
5×LTA
(A)
Wherein, C
3for short-time average coefficient, C
4for mean coefficient time long, C
5for trigger threshold.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610090278.1A CN105527650B (en) | 2016-02-17 | 2016-02-17 | Microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610090278.1A CN105527650B (en) | 2016-02-17 | 2016-02-17 | Microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105527650A true CN105527650A (en) | 2016-04-27 |
| CN105527650B CN105527650B (en) | 2017-10-17 |
Family
ID=55769986
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610090278.1A Expired - Fee Related CN105527650B (en) | 2016-02-17 | 2016-02-17 | Microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105527650B (en) |
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106199703A (en) * | 2016-08-26 | 2016-12-07 | 中国矿业大学 | A kind of microseism focus is automatically positioned and Reliability Synthesis evaluation methodology |
| CN107479094A (en) * | 2017-09-18 | 2017-12-15 | 辽宁工程技术大学 | A kind of method for realizing earthquake pre-warning |
| CN107992802A (en) * | 2017-11-10 | 2018-05-04 | 桂林电子科技大学 | A kind of microseism weak signal recognition methods based on NMF |
| CN108254781A (en) * | 2018-01-31 | 2018-07-06 | 中国科学院武汉岩土力学研究所 | A kind of dynamic adjustment of rock burst signal threshold value it is more when window complete form recognizer |
| CN108319915A (en) * | 2018-01-31 | 2018-07-24 | 中国科学院武汉岩土力学研究所 | A kind of dynamic adjustment of rock burst signal threshold value it is more when window reduced form recognizer |
| CN108376245A (en) * | 2018-02-02 | 2018-08-07 | 广西师范大学 | Time-space serial image focus recognition methods based on the channels UD |
| CN108426949A (en) * | 2018-02-14 | 2018-08-21 | 国家海洋局第二海洋研究所 | A kind of seabed sediment acoustics original position data first arrival identification pick-up method |
| CN110146920A (en) * | 2019-06-26 | 2019-08-20 | 广东石油化工学院 | Microseismic event detection method and system based on the opposite variation of amplitude |
| CN110297271A (en) * | 2019-06-26 | 2019-10-01 | 中国矿业大学 | A kind of simple component probe P wave first arrival-time modification method for mine shake alarm |
| CN112364296A (en) * | 2020-11-17 | 2021-02-12 | 东北大学 | P wave arrival time automatic picking method based on deep learning |
| CN112526602A (en) * | 2020-11-16 | 2021-03-19 | 重庆大学 | P wave arrival time picking method based on long and short time windows and AR model variance surge effect |
| CN113239620A (en) * | 2021-05-11 | 2021-08-10 | 中国电建集团华东勘测设计研究院有限公司 | Improved particle swarm method for parameter identification of geotechnical material constitutive model based on GPU acceleration |
| CN114333193A (en) * | 2022-01-11 | 2022-04-12 | 广西北投交通养护科技集团有限公司 | Portable emergent reputation alarm of road traffic and construction side slope monitoring devices |
| CN115146678A (en) * | 2022-07-04 | 2022-10-04 | 长江水利委员会长江科学院 | P-wave vibration phase first-arrival identification method and system for blasting vibration signal |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103995290A (en) * | 2014-06-03 | 2014-08-20 | 山东科技大学 | High-precision automatic microseism P wave seismic phase first arrival pickup method |
| CN104062677A (en) * | 2014-07-03 | 2014-09-24 | 中国科学院武汉岩土力学研究所 | Multifunctional comprehensive integrated high-precision intelligent micro-seismic monitoring system |
| CA2515345C (en) * | 2003-02-08 | 2014-12-09 | Robert Hughes Jones | Estimating the time of arrival of a seismic wave |
| CN104777512A (en) * | 2014-12-31 | 2015-07-15 | 安徽万泰地球物理技术有限公司 | Seismic S-wave automatic pickup algorithm |
| CN104914468A (en) * | 2015-06-09 | 2015-09-16 | 中南大学 | Mine micro-quake signal P wave first arrival moment joint pickup method |
| CN105223614A (en) * | 2015-09-23 | 2016-01-06 | 中南大学 | A kind of signals and associated noises P ripple first arrival kurtosis pick-up method based on DWT_STA/LTA |
-
2016
- 2016-02-17 CN CN201610090278.1A patent/CN105527650B/en 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 (en) * | 2014-06-03 | 2014-08-20 | 山东科技大学 | High-precision automatic microseism P wave seismic phase first arrival pickup method |
| CN104062677A (en) * | 2014-07-03 | 2014-09-24 | 中国科学院武汉岩土力学研究所 | Multifunctional comprehensive integrated high-precision intelligent micro-seismic monitoring system |
| CN104777512A (en) * | 2014-12-31 | 2015-07-15 | 安徽万泰地球物理技术有限公司 | Seismic S-wave automatic pickup algorithm |
| CN104914468A (en) * | 2015-06-09 | 2015-09-16 | 中南大学 | Mine micro-quake signal P wave first arrival moment joint pickup method |
| CN105223614A (en) * | 2015-09-23 | 2016-01-06 | 中南大学 | A kind of signals and associated noises P ripple first arrival kurtosis pick-up method based on DWT_STA/LTA |
Non-Patent Citations (2)
| Title |
|---|
| 刘晗 等: "微震信号自动检测的STA／LTA算法及其改进分析", 《地球物理学进展》 * |
| 吴治涛 等: "STA/LTA算法拾取微地震事件P波到时对比研究", 《地球物理学进展》 * |
Cited By (22)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106199703B (en) * | 2016-08-26 | 2018-11-16 | 中国矿业大学 | A kind of automatic positioning of microseism focus and Reliability Synthesis evaluation method |
| CN106199703A (en) * | 2016-08-26 | 2016-12-07 | 中国矿业大学 | A kind of microseism focus is automatically positioned and Reliability Synthesis evaluation methodology |
| CN107479094A (en) * | 2017-09-18 | 2017-12-15 | 辽宁工程技术大学 | A kind of method for realizing earthquake pre-warning |
| CN107992802A (en) * | 2017-11-10 | 2018-05-04 | 桂林电子科技大学 | A kind of microseism weak signal recognition methods based on NMF |
| CN107992802B (en) * | 2017-11-10 | 2021-12-21 | 桂林电子科技大学 | NMF-based microseism weak signal identification method |
| CN108319915B (en) * | 2018-01-31 | 2020-06-19 | 中国科学院武汉岩土力学研究所 | Multi-time-window simplified form identification method for dynamically adjusting rock burst signal threshold |
| CN108319915A (en) * | 2018-01-31 | 2018-07-24 | 中国科学院武汉岩土力学研究所 | A kind of dynamic adjustment of rock burst signal threshold value it is more when window reduced form recognizer |
| CN108254781A (en) * | 2018-01-31 | 2018-07-06 | 中国科学院武汉岩土力学研究所 | A kind of dynamic adjustment of rock burst signal threshold value it is more when window complete form recognizer |
| CN108376245A (en) * | 2018-02-02 | 2018-08-07 | 广西师范大学 | Time-space serial image focus recognition methods based on the channels UD |
| CN108426949A (en) * | 2018-02-14 | 2018-08-21 | 国家海洋局第二海洋研究所 | A kind of seabed sediment acoustics original position data first arrival identification pick-up method |
| CN110146920A (en) * | 2019-06-26 | 2019-08-20 | 广东石油化工学院 | Microseismic event detection method and system based on the opposite variation of amplitude |
| CN110297271A (en) * | 2019-06-26 | 2019-10-01 | 中国矿业大学 | A kind of simple component probe P wave first arrival-time modification method for mine shake alarm |
| CN110297271B (en) * | 2019-06-26 | 2020-09-11 | 中国矿业大学 | Single-component probe P wave first arrival time correction method for mine earthquake alarm |
| CN112526602A (en) * | 2020-11-16 | 2021-03-19 | 重庆大学 | P wave arrival time picking method based on long and short time windows and AR model variance surge effect |
| CN112526602B (en) * | 2020-11-16 | 2023-10-20 | 重庆大学 | P-wave arrival time pickup method based on long and short time windows and AR model variance surge effect |
| CN112364296B (en) * | 2020-11-17 | 2023-08-04 | 东北大学 | P-wave arrival time automatic pickup method based on deep learning |
| CN112364296A (en) * | 2020-11-17 | 2021-02-12 | 东北大学 | P wave arrival time automatic picking method based on deep learning |
| CN113239620A (en) * | 2021-05-11 | 2021-08-10 | 中国电建集团华东勘测设计研究院有限公司 | Improved particle swarm method for parameter identification of geotechnical material constitutive model based on GPU acceleration |
| CN114333193A (en) * | 2022-01-11 | 2022-04-12 | 广西北投交通养护科技集团有限公司 | Portable emergent reputation alarm of road traffic and construction side slope monitoring devices |
| CN114333193B (en) * | 2022-01-11 | 2024-04-16 | 广西北投交通养护科技集团有限公司 | Portable emergency acousto-optic alarm for road traffic and construction slope monitoring device |
| CN115146678B (en) * | 2022-07-04 | 2023-05-09 | 长江水利委员会长江科学院 | P-wave vibration phase first arrival identification method and system for blasting vibration signals |
| CN115146678A (en) * | 2022-07-04 | 2022-10-04 | 长江水利委员会长江科学院 | P-wave vibration phase first-arrival identification method and system for blasting vibration signal |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105527650B (en) | 2017-10-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105527650A (en) | Automatic identification algorithm for microseismic signal and p wave first arrival at engineering scale | |
| Ge | Efficient mine microseismic monitoring | |
| CN110109895A (en) | Fender graded unified prediction and application suitable for TBM driving tunnel | |
| CN110823356B (en) | Distributed optical fiber intrusion detection method based on Mel frequency spectrum | |
| US20100073163A1 (en) | Method and apparatus for monitoring a structure | |
| CN105261136B (en) | The method and device of weather interference is shielded in a kind of fiber-optic monitoring warning system | |
| CN106437843B (en) | coal mine bottom plate water guide channel identification method based on microseismic monitoring | |
| CN103135131A (en) | Device for interpreting fractured reservoir prediction | |
| Akhouayri et al. | Automatic detection and picking of P-wave arrival in locally stationary noise using cross-correlation | |
| CN109120336B (en) | False alarm prevention and false alarm prevention method based on phase sensitive optical time domain reflection sensor | |
| CN104459763A (en) | Method and system for detecting position of underground cavity through compactly supported wavelet | |
| CN106382981A (en) | Single station infrasonic wave signal recognition and extraction method | |
| CN106646610B (en) | A kind of algorithm using polarization constraints AIC algorithm automatic Picking microseism first arrivals | |
| CN104536046B (en) | Epicenter excitation signal conformance evaluation methodology based on earthquake record | |
| CN111520193B (en) | Non-contact tunnel engineering construction rock burst real-time forecasting method | |
| CN104977602A (en) | Control method and apparatus for earthquake data acquisition construction | |
| GB2450163A (en) | Detecting the location of seismic events without picking events in received seismic wave data | |
| CN204989215U (en) | Wheel detector signal processing system | |
| CN113050168B (en) | Crack effectiveness evaluation method based on array acoustic logging and acoustic remote detection logging data | |
| JIA et al. | Joint arrival-time picking method of microseismic P-wave and S-wave based on time-frequency analysis | |
| CN114415231A (en) | Microseismic positioning method based on EDT surface probability distribution function of station | |
| CN111123356A (en) | Abnormal track intelligent identification method based on first arrival information | |
| CN104406681B (en) | A kind of method of testing determining microseism velocity of wave in real time | |
| Nabatov | Information entropy as an identifier in rock mass structure determination using low-frequency radars | |
| CN102830426A (en) | Method and device for monitoring tunnel geology |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20171017 |