Determination of First Arrival Wave Type of Microseismic Signals and Approach to Wave Velocity Correction

Given the complex environment experienced in working mines, the vibration waves produced by processes such as rock fracture in deep formations usually show interference effects when monitored due to other signals, the so-called “clutter” in the signal, which are interfered with the clutter. At the same time, owing to the influence of system noise, the first arrival time and the arrival time difference values of the signals obtained cannot easily be determined accurately. The propagation model for the microseismic signals experienced and the discrimination method used to determine the first arrival wave type can be established using knowledge of the spatial geometry between the sensors used and the seismic source. Thus, the filtering of the actual from the abnormal wave signals is possible. Using the theory of signal cross-correlation in this work, a correction method for the arrival velocity of the first microseismic signal has been proposed and evaluated. By calculating the cross-correlation coefficient of the same source vibration signal and finding the position that corresponds to the maximum value of the cross-correlation coefficient, the arrival time difference between the signals seen in the two channels is obtained. Thus, the key conclusions can be drawn from the experiments carried out: when the signal-to-noise ratio of the original signal is low, the time difference can still be determined with high accuracy. Further, a wave velocity correction criterion has also been proposed, where the velocity correction of the S wave or the R wave can be realized by combining the spatial coordinate information on the blasting point and an algorithm representing the signal cross-correlation to arrival time difference is used.


Introduction
Enhancing the technology of microseismic monitoring in mines is very important, to allow for better forecasting of potential problems and thus creating an early warning of a possible rock burst disaster, all with a view to improving mine safety.To create an effective early warning of rock burst events, it is important to understand better both the microseismic and the blasting signals obtained from working mines.e environment experienced in mines is complex, where the vibration waves experienced are produced by processes such as rock fracture in deep formations.ese usually show interference when monitored, due to the presence of other signals, the so-called signal clutter. is noise problem also means that the first arrival time and the arrival time difference values of the signals obtained cannot be easily or accurately picked up.e propagation model of the microseismic signals experienced and the discrimination method to determine the first arrival wave can be established according to the spatial geometry between the sensors used and the seismic source.us, the filtering of the actual from the abnormal wave signals can be realized.
e linear location algorithm created, based on the signal arrival time, is usually sensitive to the arrival time difference of different signals, which then requires the use of a highprecision process to extract uniquely the actual first arrival time signal.When the background noise in the channel is large, thus, the first arrival time information is submerged in the noise, and the first arrival information may be filtered out when a filtering method is used.When only a small amount of microseismic data is present, the first arrival time of the signal can readily be obtained manually.However, at present, with the wider use of microseismic monitoring systems, the effects of the major events generated by multiple systems make it virtually impossible to use the manual approach effectively.Based on the Allen automatic pick-up algorithm (using the ratio of the short-time window average value to the long-time window average value (STA/LTA)) or the improved algorithm [1][2][3][4], when the signal-to-noise ratio of the monitored signal is low, the first arrival point of the primary (P) wave may be obscured by the presence of the noise signal, resulting in a large first arrival signal pick-up error, or even picking-up the wrong signal [5].
Zhu et al. [5] have, in their work, determined the P wave phase of all microseismic signals monitored by using the modified Akaike information criterion and thus calculated the corresponding waveform parameters.
e prediction model created was established with the support vector machine to classify both the valid and invalid P wave results picked up.
e practical application of this method has shown that the invalid P wave signals which have been picked up have been both correctly and effectively identified.Gong et al. [6] have proposed the use of the shearlet-AIC method which is based on the shear wave transformation of the Akaike information criterion.Compared with other methods mentioned in their paper, this method can accurately pick up the initial arrival time of the P wave, even when the signal-to-noise ratio is as low as −8 dB.Guo et al. [7] have introduced a fast scanning method into the field of microseismic monitoring and thus established a three-dimensional fast scanning algorithm in a Cartesian coordinate system and then used the fast scanning algorithm to calculate the first arrival travel time of a signal from a point source in the single velocity model and horizontally layered model, respectively.e method of source location based on the time difference database established with the FSM algorithm has been proposed; all this improves the location accuracy and shortens the time of location.
e velocity model is the key factor that affects the accuracy of the microseismic source location algorithm [5].e first arrival wave types from mine microseismic signals include the P wave, secondary (S) wave, or Rayleigh (R) wave.Generally, it is assumed that the whole rock mass has the same elastic modulus characteristic: therefore, a single velocity model is adopted.In many cases, the first arrival wave type is simplified to the P wave for processing [8,9], and the P wave velocity is used as the basis for the location of the source.e uncertainties in the actual rock mass velocity model and the ignorance of the anisotropy of the medium will lead to systematic-and thus unacceptable-source location errors.Guo et al. have introduced the multitemplate fast marching method (MSFM) to calculate the first arrival time for a complex rock mass containing cavities and different velocity zones in their work [10].e results demonstrate the superior value of the MSFM algorithm as a travel time forward method used in actual applications.
A 5% error in the velocity model will, however, give rise to a large location error in use.e single velocity model is problematic for the treatment of complex rock masses seen frequently in mine engineering, and therefore, it is important to create a complex rock mass velocity model that can be applied practically.e application of the isotropic velocity model to the actual anisotropic structure will lead to unrealistic velocity values.erefore, it is vitally important to consider the anisotropy of the rock mass, measure the anisotropic parameters of the velocity, and then use the anisotropic velocity model to improve the location precision of the actual microseismic source [11][12][13][14].Belayouni et al. have simulated the propagation of a direct wave and a refracted wave by designing a scientific layered wave velocity model [15] and using a ray-tracing algorithm to deal with microseismic events.e establishment of an elastic wave velocity model is a technical problem that must be overcome to realize the location of microseismic source, with high precision.Cong et al. [16] have studied the wave travel time difference between the measurement point and the reference point, using the concept of the isochromatic time surface and based on the assumption of both horizontal layered media conditions and the earliest observation point of the first arrival time in the middle of the observation network.Taking the minimum difference (double time difference) between the measured time difference and the calculated time difference of the first arrival as the constraint condition, the objective function can be constructed to solve the velocity model.Using the microseismic data obtained from known sources, the layered velocity model in horizontally layered media can be obtained by use of this method.Ouyang et al. have used the self-excited microseismic monitoring technology, with the self-excited source used to transmit the vibration signal, to retrieve and monitor the regional wave velocity field.Using the wave velocity after inversion to locate the source, the location error can be reduced to be less than 10 m. e disadvantage of the use of this method is that the self-excited source can only transmit signals over a specific frequency band and thus there is no judgment of the type of first arrival wave [17].
e above methods do not distinguish the different types of the first arrival waves and treat these first arrivals as P waves.Ge and Kaiser [18] have put forward a dynamic velocity model based on events observed.e P wave, S wave, and abnormal wave are distinguished according to the time when they are picked up.
e criterion used is the sequence and time difference of the signals received in each channel.However, when the type of the first arrival wave is determined to be an S wave or an R wave, because the S wave and the R wave of blasting signal in the mine are not obvious, they are easy to be missed due to interference by noise and the P wave follow-up tail wave and the first break time are more difficult to assess.us, the required accurate wave velocity data are difficult to obtain.e S wave or the R wave velocity can be calculated approximately by using the P wave velocity and Poisson's ratio, the modulus of elasticity and lame's coefficient, and so on.Nevertheless, the geological lithology of each mine is different, and thus, the elastic constant values from the different types of mines or indeed 2 Shock and Vibration different areas of the same mine can be significantly different.en, fixed values are used to estimate the S wave and the R wave velocities and the calculated values of S wave and R wave velocities will show huge errors.
In the situation where the P wave velocity is still used to locate the source, a huge location error will occur if there are errors in picking out the different waves and thus the type of wave velocities corresponding to the first arrival of the signal not being a P wave.erefore, it is critically important to study and thus assess the type of the wave velocity at the first arrival point of the signal for the case where the initial judgment is that the signal is not a P wave, taking the targeted calculation method for the arrival time difference to obtain the arrival time difference between the two sensors used and then to establish the wave velocity correction criteria for the S wave or the R wave.Compared with the case where there is no proper determination of the first arrival wave type and velocity correction applied, the algorithm proposed in this paper can be used greatly to improve the accuracy of the calculation of the signal delay time and thus finally to achieve the appropriate velocity correction of the arrival of the first wave.

Discrimination of the Type of First Arrival Wave
ere are different strata involved in the process of mine microseismic signal transmission, and as joints and fractures are developed in the strata, the P wave velocity will be different even in different directions in the same strata.erefore, using a single P wave or S wave velocity to locate the source will give rise to huge errors, as the velocity from the source to each sensor used will then be different.
As shown in Figure 1, it is well known that any three points in the space will lie on the same plane.In this illustration, the source is defined as S(x 0 , y 0 , z 0 , t 0 ), where its spatial coordinate is (x 0 , y 0 , z 0 ), and the time of occurrence is t 0 .e two sensors are positioned at S 1 and S i , and the times of receiving the vibration signal are t 1 and t i , respectively, where t 1 < t i .
For a microseismic event, the channel with the largest signal amplitude is chosen in which the time of receiving the vibration signal is t 1 .If the vibration wave field is an isotropic medium and given that the sum of two sides of the triangle formed (as can be seen from Figure 1) is larger than the third side, then When the source and two sensors lie on the same straight line, equation (1) becomes ( e above relationship holds when the wave velocity is given by v p .e theoretical limit of the time difference between the two sensors (for the P wave) is given by TLP 1i [19]; that is, When both sensors receive S wave signals, the theoretical limit of the time difference between the two sensors is then given by TLS 1i ; thus, where v s is the S wave velocity.When both sensors receive R wave signals, the theoretical limit of time difference between the two sensors is given by TLR 1i ; that is, where v r is the velocity of the R wave, when (t i − t 1 ) > TLR 1i , the value of t i may be affected by abnormal signals such as the R wave signal, a wave with a lower velocity than the R wave, and the delayed wave.First, it is assumed that t i is associated with the R wave signal and t 1 for the P wave signal; then, when the source and two sensors are in the same straight line, Here, TLPR 1i is the theoretical limit in the value of the R wave to the time difference at sensor S i .r 1 is the distance between the sensor and the source corresponding to the minimum value.Although it is unknown, it can be replaced by estimating the maximum value.As shown in Figure 1, the maximum distance r 1 ' from sensor S 1 to the boundary of the monitoring area is usually used (instead of r 1 ); thus, erefore, there must be a different discrimination relationship to the different wave velocity types, as shown in Figure 1: Position relationship between the seismic source and the sensors.
Table 1.It can be seen from Table 1 that, by monitoring the relationship between the time difference and the above parameters, the signal wave type received by the sensor can be determined.us, when (t i − t 1 ) ≤ TLP 1i , the signals received by both sensors are P wave types.e P wave velocity is used to calculate the distance difference between the two sensors and the source.When TLP 1i < (t i − t 1 ) ≤ TLS 1i , it can be determined that the S i sensor receives an S wave, while the S 1 sensor receives a P wave or an S wave.When calculating the distance difference between the two sensors and the source, by using the time difference between them, the S wave velocity is used.Consequently, the P wave or S wave received by S 1 does not affect the subsequent calculation.Knowing that TLS 1i < (t i − t 1 ) ≤ TLR 1i , it can be determined that the S i sensor receives the R wave, while the S 1 sensor may receive the P wave, S wave, or R wave.When calculating the distance of the two sensors from the source by using the time difference between them, the R wave velocity is used; consequently, the P wave, S wave, or R wave received by S 1 does not affect the subsequent calculation.As TLR 1i < (t i − t 1 ) ≤ TLPR 1i , the S 1 sensor receives a P wave or an S wave, while the S i sensor receives an R wave.When the distance difference between the two sensors and the source is calculated by using knowledge of the time difference, the R wave velocity is used.When (t i − t 1 ) > TLPR 1i , it is considered to be the interference of an abnormal wave.Based on the above analysis and according to the arrival time and sensor coordinate information, the wave type information corresponding to the monitoring signal used can be identified, and thus the wave velocity model used in the subsequent source positioning can be established.

Traditional Method of Determining the First Arrival Time of Microseismic Signal
Using the STA/LTA method, the monitoring system grabs the relevant events from the continuous sample points and selects the most effective channels, so as to locate the source according to the effective channel data.e traditional method, which uses the amplitude of the monitoring signal to judge the first arrival of the signal, is affected by any noise that may be present in the system.e signal waveforms obtained from two microseismic events received by a sensor are shown in Figure 2. If the threshold value is set to 0.003 m/s 2 , the signal starting point in Figures 2(a) and 2(b) is the "take-off point 2," while the signal starting point in Figure 2(b) with the same threshold setting is basically in the signal peak area.It is obvious that the error in the signal take-off point is too large, when judged according to the threshold value method.
When using the sliding energy ratio method to determine the signal starting point, the signal time corresponding to the maximum value and the minimum value of the short-and long-time window ratio is taken as the signal starting time and the signal ending time, respectively.According to the STA/LTA method, the take-off points of both events are at the "take-off point 1" in their respective figures, which greatly reduces the error compared with the "take-off point 2." However, the STA/LTA method cannot ensure that the take-off point and signal endpoint can be determined accurately, as is required every time.As shown in Figure 3, for the signal take-off point and end time point determined by the use of the sliding energy ratio method, the red vertical line represents the take-off point, and the blue vertical line represents the signal end line.It is clear from the figure that the picking of the signal end time points of channel 2 and channel 5 is incorrect.
Traditionally, most of the microseismic monitoring systems use a manual method of picking to determine the take-off point of the signal.is is feasible only when the amount of data (that then can be processed manually) is small.Now, with the wider use of microseismic monitoring systems, if the massive data sets generated by multiple systems were to be manually processed, the workload would be greatly increased.However, the error manifested in Figure 3 may appear when STA/LTA is used for automatic picking, which then requires a secondary manual correction to be applied in a large number of cases.
When the signal's first arrival cannot be determined accurately, due to the influence of system noise, the conventional automatic time arrival method or manual method, based on long-and short-time windows, will result in different levels of picking errors in this parameter.

Cross-Correlation Algorithm of Arrival Time Difference of Microseismic Signal
Because the S wave and the R wave of microseismic signals from the mine scale are not so obvious, it is more difficult to determine their first arrival times.In Table 1, when (t i − t 1 ) > TLPR 1i , the signal is preliminarily determined as an abnormal wave; however, if that is determined to be an effective vibration wave after the manual inspection is carried out and the signal cannot be picked up accurately due to the interference seen from system noise, a better, more effective arrival time difference calculation method is needed for verification of the signal.

Signal Cross-Correlation Algorithm.
In the signal analysis process carried out, the degree of correlation between the two event sequences, that is, the degree of correlation between the values of random signals x(t) and y(t) at any two different times, given by t 1 and t 2 , is determined.e correlation coefficient R of x and y is given by where Cov(x, y) is the covariance of x and y, D(x) is the variance of x, and D(y) is the variance of y. e formula for calculating the correlation value of the finite discrete signals is given by If τ is set as the phase difference of the two signals, then 4 Shock and Vibration en, the correlation number function is given by where ρ xy (τ) is a function that changes in value between 0 and 1.In practical applications, x(t) and y(t) need to be discretized to facilitate the processing that needs to be carried out.
us, x(t) � x(n) and y(t) � y(n), are set where n � 1, 2,. ..,N, and N is the total number of samples along the time axis, where R is the intermittent time shift value, r � 1 − N, 2 − N,. ..,0,. ..,N − 2, N − 1. e cross-correlation function of the finite sampling discrete signal is given as follows: e events detected by the microseismic monitoring system used are the responses of each channel to the same source; consequently, they are all related.By calculating the correlation coefficient of the signals between the two channels, determining the value of Rmax corresponding to the maximum value of the cross-correlation coefficient, and dividing by the sampling frequency F s , the delay or advance time, t d , between the two signals in the two channels can be obtained.us, t d >0 indicates that the x signal appears and the delay time t d occurs before the y signal appears, and t d <0 indicates that the y signal appears first.

Simulation Testing.
In order to test the recognition ability of signal cross-correlation to the signal delay, the following signal simulation is carried out.e sampling frequency is set to be F s � 500 Hz; the sampling number is given by N � 1024; and n � (−1) N/2 : 1: N/2−1.e frequency of the signal 1(x 1 ) is f 1 � 70 Hz, and the amplitude a 1 � 10. e frequency of the signal 2 (x 2 ) f 2 � 20 Hz, and the amplitude a 2 � 5. e delay number of signal 2 compared with signal 1 is given by d � 20, where the corresponding delay time is 0.04 s. e formula for the generation of signal 1 is then given by Signal 2 is generated in the same way, as shown in Figure 4.  6

Shock and Vibration
Using the cross-correlation function approach for the two signals, the signal cross-correlation coefficient can be obtained as shown in Figure 5.
As can be seen from Figure 5, the maximum value of the cross-correlation coefficient is 170.1, and the corresponding time is 2.088 s; thus, accordingly, the delay time is erefore, an accurate determination of the delay time between the two signals can be obtained by the use of a signal cross-correlation algorithm.
As can be seen from Figure 5, the maximum value of the cross-correlation coefficient is 170.1, and the corresponding time is 2.088 s; thus, accordingly, the delay time is t d � 2.088 − N/F s � 2.088 − 1024/500 � 0.04 s.
erefore, an accurate determination of the delay time between the two signals can be obtained by the use of a signal cross-correlation algorithm.
In order to verify the adaptability of the signal crosscorrelation method, adding random noise, represented by n 1 with a mean value of 0, and a variance σ^2 � 1, and standard deviation σ � 1, which conforms to a normal distribution, the signal after adding noise is then shown in Figure 6.
After the random noise with standard deviation σ � 1 is added, it can be seen from the signal time-domain diagram that the signal-to-noise ratio is exceedingly low.Because it is random noise, each calculation then arrives at a different result.Figure 7 shows the cross-correlation function diagram, and thus the maximum value of the cross-correlation coefficient is seen to be 228.8, with the corresponding time being 2.09 s. e delay time is then given by t d � 2.09 − N/ F s � 2.09 − 1024/500 � 0.044 s.
e added noise results in a signal delay error occurring.erefore, noise, with a mean value of 0 and σ of 0.1, 0.2, 0.5, and 1 are added, respectively, as discussed below.Each group is measured 10 times, and the signal delay time obtained is illustrated in Table 2.When σ is 0.2, the effect of the added noise signal is demonstrated in Figure 8.
It can be seen that, with the increase in the standard deviation of the noise signal, the amplitude of the noise signal increases.When the standard deviation of the noise signal is given by 0.5, the maximum error of the delay time of signal 2, relative to signal 1, is only 10% of the true value, and when the standard deviation is given by a value of 1, the error increases, but there is also a 70% probability that the delay error of signal is less than 20%.
When the original vibration signal is affected by noise, the signal cross-correlation method can be used to greatly improve the calculation accuracy of the signal delay time, compared with the manual picking method and the STA/ LTA method.

Velocity Correction Method.
e part of the monitoring signal with large energy has an enormous influence on the correlation number; therefore, the signal cross-correlation time difference is seen to be close to the time difference, when both signals are either S waves or R waves.It can be considered that the cross-correlation of the two signals is actually a comparison of the same type wave.e corresponding wave velocity, v d , of the time difference can be calibrated by the blasting process.In Table 1, when the arrival time difference corresponding wave type is not a P wave, the arrival time difference between the S i and the S 1 values is recalculated by the use of the signal cross-correlation algorithm.
When the signal velocity received by the two sensors is given by v d , the theoretical limit of the time difference between the two sensors is recorded as TLD 1i , and so

Shock and Vibration
When t d > TLD 1i , it is evident that the sensor, S i , receives abnormal signals, such as delayed waves that arrive at the sensors.
erefore, the wave velocity corresponding to different types of waves discussed in Table 1 can be modified, as is shown in Table 3.
e procedure figure of velocity correction of the first arrival wave is shown in Figure 9.

Field Testing in a Deep
Mining Environment.In order to verify the approach discussed for the determination of the first arrival wave type and the wave velocity correction, proposed in this paper, a field test has been carried out in a deep mining environment.In this case, eight microseismic sensors were arranged around the test working coalface, designated by numbers of 1#, 2#,...,8#.In Figure 10, the green dots indicate the sensors installed in the roof, and the blue dots indicate the sensors installed in the floor of the mine.
e determination of the first arrival wave type and the wave velocity correction analysis were carried out for the two microseismic events received.Five channels receive the signal from event 1; however, the first arrival of the wave recorded by sensor 2# is not obvious.e signal from four channels was received for event 2, and the wave recorded by the sensor 2# is noisy, so the first arrival time of the P wave is not easily picked out.erefore, a signal cross-correlation algorithm should be used to calculate the time difference between the two events.Sensor 8# receives the signal first; therefore, i in the first column on the left in Table 4     4 that the time difference between the signals recorded by all the sensors, except sensor 2# and 8#, is less than TLP 1i , and therefore, the first arrival wave they receive is the P wave.e waveforms recorded by sensor 2# of the above two events are based on the arrival time difference from the signal cross-correlation algorithm, and the value of t 1i is all less than TLD 1i , where the wave velocity v d is adopted accordingly.When the first arrival    point is not obvious or the signal-to-noise ratio is low, the correction method for the first arrival wave velocity, based on the signal cross-correlation, allows the determination of the first arrival wave type and the wave velocity correction.

Conclusions
A propagation model for the microseismic signal received has been established according to the arrival time information and the spatial coordinate information for the sensors used, and the judgment criterion for the first arrival wave type has been established.us, it can be judged that the first arrival wave may be P wave, S wave, R wave or abnormal wave.When signal-to-noise ratio is low and the type of first arrival wave is not a P type wave, the arrival time difference between the signals received at the channels of the microseismic system could be calculated by use of the signal cross-correlation algorithm.e simulation results illustrate that the use of the signal cross-correlation algorithm can greatly improve the accuracy of the calculation of the signal delay time.Taking into account the cross-correlation algorithm for the signal arrival time difference, a correction table for the first arrival wave velocity has been established.e first arrival wave velocity is simplified to be v p , v d , and the abnormal wave velocity.When the first arrival time of the signal cannot be determined accurately, an accurate calculation of the time difference to the arrival can be achieved.e error triggered by using a single average wave velocity is thus avoided, which makes it feasible to accurately locate the source with high precision in the rock structure.
e field test results obtained and reported here show clearly that the method proposed in this paper is not only feasible but an effective tool.

Figure 2 : 2 )Figure 3 :
Figure 2: First arrival point based on the use of different algorithms.

Figure 7 :
Figure 7: Illustration of the cross-correlation coefficient of signal with noise.

Figure 10 :
Figure 10: Planar graph showing the arrangement of the monitoring points during the field tests.

Figure 9 :
Figure 9: e procedure of velocity correction of the first arrival wave.

Table 1 :
Determination of wave velocity type.Corresponding wave type of r 1 Corresponding wave type of r i Corresponding wave type of arrival time difference ti − t 1 > TLP 1i P P P TLP 1i < t i − t 1 ≤ TLS 1i P or S S S TLS 1i < t i − t 1 ≤ TLR 1i P, S, or R R R TLR 1i < t i − t 1 ≤ TLPR 1i

Table 2 :
Delay time at different noise levels.
represents the other 4 sensors except 8#; 1 in the first column on the left in Table4represents 8#. e correction results for the first break types of events 1 and 2 are shown in Table4.After field calibration, v p is set to be 3520.25 m/s.It can be seen from Table

Table 3 :
Correction table relating to the first arrival wave velocity.Corresponding wave of r 1Corresponding wave of r i Corresponding wave of arrival time differencet i − t 1 ≤ TLP 1i v p v p v p TLP 1i < t i − t 1 � t d ≤ TLD 1i v d v d v d t i − t 1 � t d > TLD 1iAbnormal wave Abnormal wave Abnormal waveObtaining v p and v d by field test 1i < t d ≤ TLD 1i

Table 4 :
Correction of velocity type of first arrival wave.