CN105527650B - Microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick - Google Patents

Abstract

The invention discloses microseismic signals under a kind of engineering yardstick and p ripple first arrival automatic identification algorithms, comprise the following steps, reading sets sample point data in timing window；Calculate characteristic function, the long short time average of sampled point；Carry out microseismic signals and its P ripple first arrivals automatic Picking judges；Parameter value caching is reset；Data in window when reading next section.The present invention is directed to microseismic signals feature under engineering yardstick, improve the recognizer of seismic field, the present invention can be accurately under automatic Picking engineering yardstick the relatively low weak signal of microseismic signals particularly to-noise ratio, the accuracy rate of automatic Picking P ripple first arrivals is improved simultaneously, and then improves the promptness and accuracy of the geo-hazard early-warnings such as rock burst, ore deposit shake, landslide.

Description

An automatic recognition algorithm for microseismic signals and p-wave first arrivals at an engineering scale

technical field

The invention relates to the technical field of microseismic monitoring. It specifically relates to an automatic identification algorithm for microseismic signals and p-wave first arrivals on an engineering scale, and the method can be widely used in mining engineering, water conservancy and hydropower engineering, petroleum engineering, geotechnical engineering, and underground engineering.

Background technique

Microseismic source location is an important part of microseismic monitoring and disaster early warning, and the identification of microseismic signals and P-wave first arrival picking are the keys to microseismic source location. Most of the existing microseismic signal identification methods come from the field of seismicity, and there are many methods, mainly including: STA/LTA algorithm in the time domain to construct characteristic functions based on energy and energy changes in the time domain; Differences, such as the AIC algorithm; in addition, there are neural network methods, high-order statistical methods, and wavelet transform methods based on wavelet theory. STA/LTA and its improved algorithm are widely used in the seismic field because of its simple algorithm, high computational efficiency, and suitable for real-time processing. These methods have good applicability for the processing of seismic signals, but the rock fracture signals generated at the engineering scale are different from natural seismic signals:

1) The frequency of seismic signals is generally below tens of Hz, while the frequency of microseismic signals is generally tens to hundreds of Hz, and some are as high as several thousand Hz; 2) The duration of seismic signals is longer, generally greater than 1 second, and the frequency of microseismic signals lasts The time is short, generally no more than 0.1 second; 3) Compared with natural earthquakes, microseismic events are generally small in magnitude, and smaller microseismic events are not easy to be picked up; 4) In engineering environments, the background is caused by factors such as mechanical construction, electrical interference, and vehicle driving. The noise is complex, making most of the detected microseismic signals have a low signal-to-noise ratio.

Therefore, it is inappropriate to use the signal identification and P-wave arrival time picking methods in the seismic field for engineering-scale microseismic signals, mainly as follows: 1) It is difficult to identify weak signals with low signal-to-noise ratio; The error is large.

To this end, in view of the above shortcomings, a method for automatic identification of microseismic signals and their first p-wave arrivals at an engineering scale is invented, which improves the picking accuracy of microseismic signals and their p-wave first arrivals at an engineering scale, thereby improving rock The timeliness and accuracy of early warning of geological disasters such as explosions, mine earthquakes, and landslides are very necessary.

Contents of the invention

The purpose of the present invention is to overcome the problems existing in the prior art, to provide an automatic recognition algorithm for microseismic signals and p-wave first arrivals at engineering scales, to improve the picking accuracy of microseismic signals and first arrivals of p waves at engineering scales, and to further improve The timeliness and accuracy of early warning of geological disasters such as rock bursts, mine earthquakes, and landslides.

An algorithm for automatic identification of microseismic signals and p-wave first arrivals on an engineering scale, comprising the following steps:

Step 1. Read the sampling point data of the waveform in the set time window monitored by the vibration sensor in real time, and establish a plane Cartesian coordinate system, in which the X-axis corresponds to the sampling point number of the waveform, and the Y-axis corresponds to the amplitude of the waveform. In the plane Cartesian coordinates Symmetrically calibrate the sampling points of the waveform in the time window in the system;

Step 2. Calculate the number N of sampling points, set the upper limit M1 of the number of zero-crossing points, and the lower limit M2 of the number of zero-crossing points; set the time judgment threshold T of the microseismic signal;

Step 3, setting the weighting factor K _(A) and characteristic function CF _(A) of the Ath sampling point;

Set the Ath sampling point short-term average value STA _(A) and long-term average value LTA _(A) ;

Set the dynamic threshold r(A) of the Ath sampling point;

Where A is the serial number of the number of sampling points;

Step 4, taking the sampling point with the smallest number in the plane Cartesian coordinate system as the current sampling point, and the number of the current sampling point is i;

Step 5, setting algorithm parameter variable M, S, L, t, the initial value of M, S, L, t is all set to 0, and wherein M is the number of zero-crossing points; S is a counter; The calculation formula of L is L= 3+M/3; t is the duration of the microseismic signal picked up;

Step 6, compare the size of the short-term average STA _(i) of the current sampling point and the dynamic threshold r (i);

Step 7. If the short-term average value STA _(i) of the current sampling point is less than the dynamic threshold r(i) of the current sampling point, then the current sampling point is not the first arrival point of the P wave. At this time, the next sampling point of the current sampling point point as the current sampling point, return to step (6); if the short-term average STA _(i) of the current sampling point is greater than or equal to the dynamic threshold r(i) of the current point, assign the number i of the current sampling point to k , calculate the amplitude P1 of the current sampling point, assign the amplitude value P1 of the current sampling point to the cache maximum amplitude P0, and then proceed to step (8);

Step 8, calculating the amplitude P2 of the next sampling point k+1 of the sampling point numbered k;

Step 9. If the next sampling point k+1 of the sampling point whose number is k is not a zero-crossing point, assign the larger value of the buffered maximum amplitude P0 and the amplitude P2 of the sampling point k+1 to the cached maximum Amplitude P0, assign k+1 to k, when the value of k is less than the number N of waveform sampling points, return to step 8; when the value of k is greater than or equal to the number N of waveform sampling points, the next One sampling point is used as the current sampling point, return to step 5;

If the number of the sampling point is k and the next sampling point k+1 is a zero-crossing point, the value M of the number of zero-crossing points is increased by 1, and k+1 is assigned to k, and then step 10 is performed;

Step 10. If the value M of the number of zero-crossing points is greater than the upper limit M1 value of the number of zero-crossing points, the current sampling point is not the first arrival point of the P wave, and the next sampling point of the current sampling point is taken as the current sampling point, and return to step 5;

If the zero-cross point number M value is less than or equal to the zero-cross point number upper limit M1 value, then calculate the judgment parameter δ=LTA(i)×M of the current sampling point, and then proceed to step 11;

Step 11, compare the size of the short-term average value STA _(k) of the sampling point k with the numbering and the judgment parameter δ of the current sampling point:

If the short-term average value STA _(k) of the sampling point that is numbered k is less than or equal to the determination parameter δ of the current sampling point, the counter S value is reset to zero, and P2 is assigned a value to P0, and k is added by 1 at this moment, and returns to step 8;

If numbering is the short-term average value STA _(k) of the sampling point of k greater than the determination parameter δ of the current sampling point, the counter S value adds 1, calculates L=3+M/3, enters step 12;

Step 12, compare the counter S value and the size of L:

If the counter S value is less than or equal to L, assign P2 to P0, add 1 to k, and return to step 8;

If the counter S value is greater than L, calculate the duration t of this section of waveform from the current sampling point i to the sampling point numbered k, t=(k-i)/f, where f is the sampling frequency;

Step 13. When t is less than the microseismic signal time judgment threshold T or the number M of zero crossing points is less than the lower limit M2 of the number of zero crossing points, the sampling point between the current sampling point i and the sampling point numbered k is not a microseismic signal , take the next sampling point of the sampling point numbered k as the current sampling point, clear the values of M, P0, P1, P2, S, and L, and proceed to step 5; when t is greater than or equal to the time threshold T and zero crosses When the number of points M is greater than or equal to the lower limit of the number of zero-crossing points M2, the sampling point between the current sampling point i and the sampling point numbered k is a microseismic signal, and the number i of the current sampling point is assigned to T1, and k is assigned to T2, wherein the sampling point numbered T1 is the first arrival point of the P wave of the microseismic signal, and the sampling point numbered T2 is the end point of the signal, and enter step 14;

Step 14. Use the next sampling point of the sampling point numbered k as the current sampling point, clear the values of M, P0, P1, P2, S, and L, and proceed to step 5 until the waveform in the set time window is Sampling point data analysis is complete.

Symmetry calibration as described above involves the following steps:

Select a set number of consecutive sampling points in the time window as calibration samples,

If the ratio of the number of sampling points with positive amplitude to the number of sampling points with negative amplitude in the calibration sample is not within the set ratio range, the sampling point data of the waveform in the set time window needs to be calibrated symmetrically, and the calibration sample is calculated The average value of the amplitude, adding the amplitude of the sampling points of the waveform in the set time window and the average value of the corrected sample amplitude to complete the symmetrical calibration;

If the ratio of the number of sampling points with positive amplitude to the number of sampling points with negative amplitude in the calibration sample is not within the set ratio range, the sampling point data of the waveform within the set time window does not need to be symmetrically calibrated.

The weighting factor K _(A) of the Ath sampling point as described above is based on the following formula:

The feature function CF _(A) of the A sampling point is based on the following formula:

Among them, A is the serial number of the number of sampling points in the waveform after symmetrical calibration, and y _(A) is the amplitude of the Ath sampling point of the waveform after symmetrical calibration, is the first order difference of y _(A) ;

The short-term average STA _(A) of the A sampling point is based on the following formula:

STA _(A) ＝STA _(A-1) +C ₃ ×[CF _(A) -STA _(A-1) ]

The long-term average LTA _(A) of the A sampling point is based on the following formula:

LTA _(A) ＝LTA _(A-1) +C ₄ ×[CF _(A) -LTA _(A-1) ]

The dynamic threshold r(A) of the Ath sampling point is based on the following formula:

r(A)＝C ₅ ×LTA _(A)

Among them, C ₃ is the short-term average coefficient, C ₄ is the long-term average coefficient, and C ₅ is the trigger threshold.

When the P wave arrives, it will cause changes in the waveform characteristics of the microseismic signal such as amplitude or frequency. If the signal-to-noise ratio of the signal is lower, the first arrival of the P wave will be less obvious. Therefore, improving the sensitivity of the picking algorithm to such changes will help to improve the first-arrival picking ability of P waves for microseismic signals on an engineering scale. In the algorithm, the weighting factor K, the characteristic function CF, and the dynamic long-short-time average ratio ε can all reflect the characteristic changes of the waveform and directly affect the picking accuracy. The dynamic long-short-time average ratio ε is the long-term average LTA _(A) Ratio of mean STA _(A) :

Analyze the sensitivity of the weighting factor K, characteristic function CF, and dynamic long-short-time average ratio ε to characteristic changes such as signal amplitude or frequency waveform of the present invention, as shown in Figures 2, 3, and 4, respectively. According to the comparison results, the three variables of the present invention change significantly at the place where the frequency or amplitude changes, and reach the peak quickly. Therefore, in theory, the present invention has an obvious effect on picking up signals with low signal-to-noise ratio and the first arrival of P waves.

Randomly select 555 microseismic signals as samples from the real-time monitoring data of Jinping Water Diversion Tunnel, and define two criteria for judging the picking results: 1) Successful signal picking: the sample signal can be automatically identified without leakage and misjudgment. This will facilitate manual further processing of the automatically picked-up signal parts. The more signals are successfully picked up, the higher the signal pick-up rate. 2) Accurate P-wave first-arrival picking: Comparing the P-wave first-arrival results picked up by the automatic algorithm with the manual picking results, if the picking error is not greater than 3 sampling points, the automatic picking result of the algorithm is considered accurate. The signal pickup rate of the present invention is 94.23%, the P wave pickup accuracy rate is 73.51%, and the pickup rate fully satisfies the current positioning error in the microseismic field.

Analyzed the microseismic monitoring data of Jinping's deep underground laboratory from June 15, 2015 to July 15, 2015, a total of 524 microseismic signals and 109 events were used for positioning by manually picking up P waves and the picking results of the present invention. As shown in Figures 6 and 7, respectively, the 7# laboratory is a high-risk area for rock fracture. It can be seen that the positioning results of the present invention have a relatively high concentration of microseismic events, and the spatial distribution of microseismic events is close to the manual picking results, which is consistent with the actual project. The situation matches. 13 events with very low signal-to-noise ratio, inconspicuous P waves, and obviously abnormal positioning results in the manual picking process, as shown in the triangle in Fig. All can be located successfully, as shown by the triangle in Fig. 7; 2) The positioning effect of the 13 events picked up by the present invention is significantly improved compared with manual picking, and most of them are within the allowable range of positioning errors.

Beneficial effect:

The present invention aims at the characteristics of microseismic signals at the engineering scale, and improves the recognition algorithm in the field of earthquakes. The present invention can accurately and automatically pick up the microseismic signals at the engineering scale, especially weak signals with low sex-to-noise ratio, and at the same time improve the efficiency of automatically picking up the first arrival of P waves. Accuracy, thereby improving the timeliness and accuracy of early warning of geological disasters such as rock bursts, mine earthquakes, and landslides.

Description of drawings

Fig. 1 is a flowchart of the present invention;

Fig. 2 is the sensitivity analysis diagram of weighting factor K of the present invention to amplitude or frequency change; As Fig. 2 (a) when amplitude changes, the weighting factor amplitude of the present invention increases approximately 2 times; As Fig. 2 (b) when frequency When a change occurs, the variation range of the K value of the present invention is more obvious.

Fig. 3 is the sensitivity analysis diagram of characteristic function CF of the present invention to amplitude or frequency variation; The amplitude at the sudden change point is from 1 to 16, and the waveform symmetry axis changes from Y=1 to Y=16 before and after the sudden change; the present invention has obvious amplitude changes and symmetrical axis deviations before and after the sudden change, and the characteristic function curve changes rapidly, as shown in Fig. 3 ( b) When the frequency of the analog signal changes, the amplitude of the characteristic function of the present invention at the mutation point ranges from 1 to 15.6, the amplitude changes obviously, and the amplitude increases rapidly.

Fig. 4 is the epsilon ratio of the present invention sensitivity analysis figure to amplitude or frequency variation; Shown in Fig. 4 (a): the epsilon variation curve of the present invention reaches peak value in 3 samplings after waveform amplitude changes, and curve variation is obvious; 4(b) When the frequency changes, the frequency change position of the present invention can be picked up manually manually, and the curve changes obviously.

Fig. 5 is the variation figure of S-L value of picking up process of example 1 of the present invention;

Figure 6. The left view of the P-wave positioning effect of artificially picked up microseismic signals;

Fig. 7 is the left side view of the P wave positioning effect of the automatic pickup of the microseismic signal by the present invention.

detailed description

In order to make the technical means, creative features, work flow, and use methods of the present invention achieve the purpose and effect easily understood, the present invention will be further described below in conjunction with specific embodiments. The scope of protection of the present invention is not limited by the following examples.

The vertical rock coverage of Jinping Underground Laboratory reaches 2,400 meters. It is currently the laboratory with the deepest rock coverage in the world. The maximum principal stress in the project area can reach 63MPa, which is a typical high-stress area. Among them, 7# and 8# laboratories are under construction. Strong to extremely strong rockbursts have been encountered, causing serious threats and losses to the safety of on-site personnel and equipment. During the on-site microseismic monitoring process, it was found that because the sensor is far away from the 7-8# laboratory construction section with a high risk of rockburst, about 400-500m, the microseismic signal is severely attenuated during the propagation process, and the background noise of the on-site monitoring is complicated. Most of the monitored microseismic signals have low signal-to-noise ratios, small energy amplitudes, and the first arrival of P waves is not obvious, making it difficult to identify microseismic signals and pick up P waves. Using the engineering microseismic sample signal, compare and analyze the results of the present invention’s picking and manual picking, and use the particle swarm algorithm to optimize the best parameter set of the present invention suitable for this project: C ₃ =0.25, C ₄ =0.08, C ₅ =2 , upper limit M1=300, lower limit M2=15, T=0.01. This example takes the analysis of real-time monitoring sampling point data in a certain time window of the project as an example to illustrate.

Example 1:

An algorithm for automatic identification of microseismic signals and p-wave first arrivals at an engineering scale, the process flow is shown in Figure 1, and the steps are as follows:

(1) Read the sampling point data of the waveform in the set time window monitored by the vibration sensor in real time, and establish a plane Cartesian coordinate system X0Y, where the X-axis corresponds to the sampling point number of the waveform, from the origin of the plane Cartesian coordinate system to the X-axis The number of the positive axis starts from 1, and the sampling points are arranged from small to large according to the number. The sampling point number is a positive integer, and the Y axis corresponds to the amplitude of the waveform. The sampling points of the waveform in the time window are symmetrically calibrated in the plane Cartesian coordinate system. Go to step (2).

Waveform symmetry calibration: select a set number (1000 in this example) of continuous sampling points in the time window as calibration samples,

If the ratio of the number of sampling points with positive amplitude to the number of sampling points with negative amplitude in the calibration sample is not within the set ratio range (in this example, it is less than 0.95 or greater than 1.05), then within the selected time window The sampling point data of the waveform needs to be calibrated symmetrically, and the average value of the corrected sample amplitude is calculated, and the sample point amplitude of the waveform in the time window set in step (1) is added to the average value of the corrected sample amplitude to complete the symmetrical calibration. This is conducive to accurately counting the number M of zero-crossing points in the waveform, and reducing the algorithm picking error.

If the ratio of the number of sampling points whose amplitude is positive to the number of sampling points whose amplitude is negative in the calibration sample is within the set ratio range (in this example, it is greater than or equal to 0.95 and less than or equal to 1.05), then step (1) The sampling point data of the waveform in the medium-set time window does not need to be symmetrically calibrated.

(2) Calculate the number N of sampling points of the waveform in the plane Cartesian coordinate system after symmetrical calibration. Set the upper limit of the number of zero crossing points M1 = 300, the lower limit of the number of zero crossing points M2 = 15; set the microseismic signal time judgment threshold T = 0.01, the unit is second. The duration of the microseismic signal is short, and setting the upper limit M1 of the number of zero crossing points and the lower limit M2 of the number of zero crossing points can effectively avoid electrical spike signal interference with a short duration and misjudgment of most non-microseismic signals with a long duration.

(3), utilize the sampling point data of the waveform in the plane Cartesian coordinate system after symmetric calibration to calculate the weighting factor K _(A) of the Ath sampling point of the present invention:

And the characteristic function CF _(A) of the Ath sampling point:

Among them, A is the number of sampling points in the waveform after symmetrical calibration, and y _(A) is the amplitude of the Ath sampling point of the waveform after symmetrical calibration, is the first difference of y _(A) .

Calculate the short-term average STA _(A) and long-term average LTA _(A) :

STA _(A) ＝STA _(A-1) +C ₃ ×[CF _(A) -STA _(A-1) ]

LTA _(A) ＝LTA _(A-1) +C ₄ ×[CF _(A) -LTA _(A-1) ]

Wherein STA _(A) is the short-term average value of the A sampling point; C ₃ is the short-term average coefficient, and C ₃ takes a value of 0.25 in the example 1; LTA _(A) is the long-term average value of the A sampling point; _C4 is the long-term average coefficient, and the value of _C4 in Example 1 is 0.08.

Calculate r(A)=C ₅ ×LTA _(A) , where r(A) is the dynamic threshold of the Ath sampling point, C ₅ is the trigger threshold, and the value of C ₅ in Example 1 is 2.

(4) Take the sampling point with the smallest number in the planar Cartesian coordinate system as the current sampling point, and the number of the current sampling point is i.

(5), setting algorithm parameter variable M, S, L, t, initial value all is set to 0, and wherein M is the number of zero-crossing points; S is counter, plays counting effect; The calculation formula of L is L=3+M /3; t is the duration of the picked-up microseismic signal.

(6) Compare the short-term average STA _(i) value of the current sampling point with the dynamic threshold value r(i).

(7), if the short-term average STA _(i) of the current sampling point is less than the dynamic threshold r(i) of the current sampling point, then the current sampling point is not the first arrival point of the P wave, and the next Take the sampling point as the current sampling point and return to step (6); if the short-term average value STA _(i) of the current sampling point is greater than or equal to the dynamic threshold r(i) of the current point, then the current sampling point may be the first arrival of P wave point, assign the number i of the current sampling point to k, which is k=i, and calculate the magnitude P1 of the current sampling point. In order to find out the first peak value after the current sampling point, the magnitude of the current sampling point P1 Assign a value to the cache maximum amplitude value P0 and temporarily store it as the maximum peak value, and then proceed to step (8).

(8) Calculate the amplitude P2 of the next sampling point k+1 whose number is k in the waveform of the coordinate system after the symmetrical calibration.

(9), if the next sampling point k+1 whose sampling point number is k in the waveform after symmetrical calibration is not a zero-crossing point, compare the maximum amplitude P0 of the cache and the amplitude P2 of the sampling point k+1, which will be compared Assign a large value to the maximum amplitude of the buffer P0, assign k+1 to k, when the value of k is less than the number N of waveform sampling points, return to step (8); when the value of k is greater than or equal to the number N of waveform sampling points , then take the next sampling point of the current sampling point as the current sampling point, that is, i=i+1, and return to step (5).

If the next sampling point k+1 of the sampling point numbered k in the coordinate system waveform after symmetrical calibration is a zero-crossing point, then the number M of the zero-crossing point is increased by 1, that is, M=M+1, and k+1 is assigned a value Give k, then go to step (10).

(10), judge if the zero-cross point number M value is greater than the zero-cross point number upper limit M1 value, then the current sampling point is not the first arrival point of the P wave, then the next sampling point of the current sampling point is used as the current sampling point, That is, i=i+1, return to step (5).

If the value M of the number of zero crossing points is less than or equal to the upper limit M1 value of the number of zero crossing points, then calculate the judgment parameter δ=LTA(i)×M of the current sampling point, and then proceed to step (11).

(11), comparing numbering is the size of the short-term average value STA _(k) of the sampling point k and the determination parameter δ of the current sampling point:

If the short-term average STA _(k) of the sampling point numbered k is less than or equal to the determination parameter δ of the current sampling point, the counter S value is reset to zero, and P2 is assigned to P0, and k is added by 1 at this time, that is, k=k+ 1. Return to step (8);

If the sampling point short-term average value STA _(k) that numbering is k is greater than the judgment parameter δ of current sampling point, counter S value adds 1, promptly S=S+1, calculates L=3+M/3, enters step (12 ).

(12), compare the size of the counter S value and L: if the counter S value is less than or equal to L, then P2 assignment is given to P0, now k is added 1, promptly k=k+1, returns to step (8); if the counter S value is greater than L, calculate the duration t of the waveform from the current sampling point to the sampling point numbered k in the coordinate system waveform after symmetrical calibration, t=(k-i)/f, where f is the sampling frequency.

(13), when t is less than the microseismic signal time judgment threshold T, T is 0.01, the unit is second, or when the number M of zero crossing points is less than the lower limit M2 of the number of zero crossing points, M2 is 15, then judge the current sampling point i The sampling point between the sampling point and the sampling point k is not a microseismic signal, and it is necessary to continue to pick up the signal after the sampling point k at this time. At this time, the next sampling point of the sampling point k in the waveform coordinate system after symmetrical calibration point as the current sampling point, clear the values of M, P0, P1, P2, S, and L, and proceed to step (5); when t is greater than or equal to the microseismic signal time threshold T and the number of zero-crossing points M is greater than or equal to the zero-crossing point When the lower limit of the number is M2, it is judged that the sampling point between the current sampling point i and the sampling point numbered k is a microseismic signal, and i is assigned to T1, and k is assigned to T2, wherein the sampling point numbered T1 is the microseismic signal At the first arrival point of the P wave, the sampling point numbered T2 is the signal end point, and the step (14) is entered.

(14) Continue to automatically analyze the next signal. After symmetrical calibration, the next sampling point of the sampling point numbered k in the waveform coordinate system is used as the current sampling point, and the values of M, P0, P1, P2, S, and L are cleared. value, proceed to step (5) until the analysis of the sampling point data of the waveform in the set time window is completed.

The basic principles and main features of the present invention and the advantages of the present invention have been shown and described above. Those skilled in the industry should understand that the present invention is not limited by the above-mentioned embodiments. What are described in the above-mentioned embodiments and the description only illustrate the principle of the present invention. Without departing from the spirit and scope of the present invention, the present invention will also have Variations and improvements are possible, which fall within the scope of the claimed invention. The protection scope of the present invention is defined by the appended claims and their equivalents.

Claims (3)

1. microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick, it is characterised in that comprise the following steps：

The sample point data for setting waveform in timing window that step 1, reading vibrating sensor are real-time monitored, sets up flat square seat The sampling point number of mark system, wherein X-axis correspondence waveform, the amplitude of Y-axis correspondence waveform, when in plane right-angle coordinate pair in window The sampled point of waveform is symmetrically calibrated；

Step 2, calculating sampled point number N, set zero cross point number upper limit M1, zero cross point number lower limit M2；Microseism is set Signal time decision threshold T；

Step 3, the weighted factor K for setting the A sampled point_(A)With characteristic function CF_(A)；

Set the A sampled point short time average STA_(A)With it is long when average value LTA_(A)；

Set the dynamic threshold r (A) of the A sampled point；

Wherein A is the numbering of sampled point number；

Step 4, the sampled point of minimum is numbered using in plane right-angle coordinate as current sampling point, the numbering of current sampling point is i；

Step 5, setting algorithm parameter variable M, S, L, t, M, S, L, t initial value are disposed as 0, and wherein M is zero cross point number； S is counter；L calculation formula is L=3+M/3；T is the duration of the microseismic signals picked up；

Step 6, the short time average STA for comparing current sampling point_(i)With dynamic threshold r (i) size；

If the short time average STA of step 7, current sampling point_(i)Less than the dynamic thresholding r (i) of current sampling point, then currently adopt Sampling point is not P ripple Onset points, now using next sampled point of current sampling point as current sampling point, returns to step (6)； If the short time average STA of current sampling point_(i)More than or equal to the dynamic thresholding r (i) of current point, then by the volume of current sampling point Number i is assigned to k, calculates the amplitude P1 of current sampling point, the amplitude size P1 of current sampling point is assigned into caching maximum amplitude P0, then carries out step (8)；

Step 8, the amplitude P2 for calculating the next sampled point k+1 for numbering the sampled point for being k；

It is if the next sampled point k+1 for the sampled point that the numbering of step 9, sampled point is k is not zero cross point, caching is maximum Higher value is assigned to caching maximum amplitude P0 in amplitude P0 and sampled point k+1 amplitude P2, k+1 is assigned into k, when k value is small When waveform sampling point number N, return to step 8；When k value is more than or equal to waveform sampling point number N, then by current sampling point Next sampled point be used as current sampling point, return to step 5；

If next sampled point k+1 that the numbering of sampled point is k is zero cross point, zero cross point number M values Jia 1, by k+1 K is assigned to, step 10 is then carried out；

If step 10, zero cross point number M values are more than zero cross point number upper limit M1 values, current sampling point is not P ripple first arrivals Point, then regard next sampled point of current sampling point as current sampling point, return to step 5；

If zero cross point number M values are less than or equal to zero cross point number upper limit M1 values, the critical parameter δ of current sampling point is calculated =LTA (i) × M, then carries out step 11；

Step 11, the short time average STA for comparing the sampled point that numbering is k_(k)With the critical parameter δ of current sampling point size：

If the short time average STA for the sampled point that numbering is k_(k)Less than or equal to the critical parameter δ of current sampling point, counter S values Zero, is assigned to P0 by P2, now adds 1 by k, return to step 8；

If the short time average STA for the sampled point that numbering is k_(k)More than the critical parameter δ of current sampling point, counter S values Jia 1, L=3+M/3 is calculated, into step 12；

Step 12, compare counter S values and L size：

If counter S values are less than or equal to L, P2 is assigned to P0, k plus 1, return to step 8；

If counter S values are more than L, the duration t of this section of waveform from the current sampling point i sampled points for being k to numbering is calculated, T=(k-i)/f, wherein f are 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 When, then current sampling point i and numbering be k sampled point between sampled point be not microseismic signals, will number be k sampled point Next sampled point resets M, P0, P1, P2, S, L value as current sampling point, carries out step 5；When t is more than or equal to the time When threshold values T and zero cross point number M are more than or equal to zero cross point number lower limit M2, then the sampling that current sampling point i is k with numbering Sampled point between point is microseismic signals, and the numbering i of current sampling point is assigned into T1, k is assigned into T2, wherein numbering is T1 Sampled point be microseismic signals P ripple Onset points, the sampled point that numbering is T2 is signal end point, into step 14；

Step 14, the next sampled point for the sampled point for being k using numbering reset M, P0, P1, P2, S, L as current sampling point Value, carries out step 5, until the sample point data for setting waveform in timing window has been analyzed.

2. microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick according to claim 1, its feature exist In described symmetrical calibration comprises the following steps：

The continuous sampled point of number is set during selection in window as calibration samples,

If the ratio that amplitude is positive sampled point number in calibration samples and amplitude is negative sampled point number is not in setting ratio In the range of, then setting the sample point data of waveform in timing window need to symmetrically calibrate, and calculate the average value of calibration samples amplitude, will set When window in the sampled point amplitude of waveform be added with the average value of calibration samples amplitude, complete symmetrical calibration；

If the ratio that amplitude is positive sampled point number in calibration samples and amplitude is negative sampled point number is not in setting ratio In the range of, then the sample point data for setting waveform in timing window is not required to symmetrical calibration.

3. microseismic signals and p ripple first arrival automatic identification algorithms under a kind of engineering yardstick according to claim 1, its feature exist In,

The weighted factor K of the A described sampled point_(A)Based on below equation：

<mrow> <msub> <mi>K</mi> <mrow> <mo>(</mo> <mi>A</mi> <mo>)</mo> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <munderover> <mo>&amp;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>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>A</mi> </munderover> <msup> <msub> <mover> <mi>y</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </msub> <mn>2</mn> </msup> </mrow> </mfrac> </mrow>

The characteristic function CF of the A described sampled point_(A)Based on below equation：

<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>&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>A</mi> <mo>)</mo> </mrow> </msub> <mn>2</mn> </msup> <mo>&amp;times;</mo> <mfrac> <mrow> <munderover> <mo>&amp;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>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>A</mi> </munderover> <msup> <msub> <mover> <mi>y</mi> <mo>&amp;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>

Wherein, A is the numbering of sampled point number in waveform after symmetrical calibration, y_(A)For the A sampled point of waveform after symmetrical calibration Amplitude,For y_(A)First-order difference；

The short time average STA of the A described sampled point_(A)Based on below equation：

STA_(A)=STA_(A-1)+C₃×[CF_(A)-STA_(A-1)]

The A described sampled point it is long when average value LTA_(A)Based on below equation：

LTA_(A)=LTA_(A-1)+C₄×[CF_(A)-LTA_(A-1)]

The dynamic threshold r (A) of the A described sampled point is based on below equation：

R (A)=C₅×LTA_(A)

Wherein, C₃For short-time average coefficient, C₄For it is long when mean coefficient, C₅For triggering threshold values.