Abstract
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.
1. 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. The environment experienced in mines is complex, where the vibration waves experienced are produced by processes such as rock fracture in deep formations. These usually show interference when monitored, due to the presence of other signals, the so-called signal clutter. This 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. The 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. Thus, the filtering of the actual from the abnormal wave signals can be realized.
The 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 high-precision 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–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. The prediction model created was established with the support vector machine to classify both the valid and invalid P wave results picked up. The 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. The 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.
The velocity model is the key factor that affects the accuracy of the microseismic source location algorithm [5]. The 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. The 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]. The 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. The 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. The application of the isotropic velocity model to the actual anisotropic structure will lead to unrealistic velocity values. Therefore, 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–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. The 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. The 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].
The 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. The P wave, S wave, and abnormal wave are distinguished according to the time when they are picked up. The 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. Thus, the required accurate wave velocity data are difficult to obtain. The 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 different areas of the same mine can be significantly different. Then, 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. Therefore, 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.
2. Discrimination of the Type of First Arrival Wave
There 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. Therefore, 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(x0, y0, z0, t0), where its spatial coordinate is (x0, y0, z0), and the time of occurrence is t0. The two sensors are positioned at S1 and Si, and the times of receiving the vibration signal are t1 and ti, respectively, where t1 < ti.
For a microseismic event, the channel with the largest signal amplitude is chosen in which the time of receiving the vibration signal is t1. 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
The above relationship holds when the wave velocity is given by . The theoretical limit of the time difference between the two sensors (for the P wave) is given by TLP1i [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 TLS1i; thus,where 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 TLR1i; that is,where is the velocity of the R wave, when (ti − t1) > TLR1i, the value of ti 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 ti is associated with the R wave signal and t1 for the P wave signal; then, when the source and two sensors are in the same straight line,
Here, TLPR1i is the theoretical limit in the value of the R wave to the time difference at sensor Si. r1 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 r1’ from sensor S1 to the boundary of the monitoring area is usually used (instead of r1); thus,
Therefore, there must be a different discrimination relationship to the different wave velocity types, as shown in 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. Thus, when (ti − t1) ≤ TLP1i, the signals received by both sensors are P wave types. The P wave velocity is used to calculate the distance difference between the two sensors and the source. When TLP1i < (ti − t1) ≤ TLS1i, it can be determined that the Si sensor receives an S wave, while the S1 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 S1 does not affect the subsequent calculation. Knowing that TLS1i < (ti − t1) ≤ TLR1i, it can be determined that the Si sensor receives the R wave, while the S1 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 S1 does not affect the subsequent calculation. As TLR1i < (ti − t1) ≤ TLPR1i, the S1 sensor receives a P wave or an S wave, while the Si 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 (ti − t1) > TLPR1i, 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.
3. 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. The 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.
The 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/s2, 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.
(a)
(b)
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. This 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.
4. 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 (ti − t1) > TLPR1i, 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.
4.1. 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 t1 and t2, is determined. The correlation coefficient R of x and y is given bywhere 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. The 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
Then, the correlation number function is given bywhere ρ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. Thus, 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. The cross-correlation function of the finite sampling discrete signal is given as follows:
The 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 Fs, the delay or advance time, td, between the two signals in the two channels can be obtained. Thus, td >0 indicates that the x signal appears and the delay time td occurs before the y signal appears, and td <0 indicates that the y signal appears first.
4.2. Simulation Testing
In order to test the recognition ability of signal cross-correlation to the signal delay, the following signal simulation is carried out. The sampling frequency is set to be Fs = 500 Hz; the sampling number is given by N = 1024; and n = (−1) N/2 : 1: N/2−1. The frequency of the signal 1(x1) is f1 = 70 Hz, and the amplitude a1 = 10. The frequency of the signal 2 (x2) f2 = 20 Hz, and the amplitude a2 = 5. The delay number of signal 2 compared with signal 1 is given by d = 20, where the corresponding delay time is 0.04 s. The formula for the generation of signal 1 is then given by
Signal 2 is generated in the same way, as shown in Figure 4.
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 td = 2.088 − N/Fs = 2.088 − 1024/500 = 0.04 s.
Therefore, 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 td = 2.088 − N/Fs = 2.088 − 1024/500 = 0.04 s.
Therefore, 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 cross-correlation method, adding random noise, represented by n1 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. The delay time is then given by td = 2.09 − N/Fs = 2.09 − 1024/500 = 0.044 s.
The added noise results in a signal delay error occurring. Therefore, 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.
5. Velocity Correction of the First Arrival Wave
5.1. Velocity Correction Method
The 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. The corresponding wave velocity, , 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 Si and the S1 values is recalculated by the use of the signal cross-correlation algorithm.
When the signal velocity received by the two sensors is given by , the theoretical limit of the time difference between the two sensors is recorded as TLD1i, and so
When td > TLD1i, it is evident that the sensor, Si, receives abnormal signals, such as delayed waves that arrive at the sensors. Therefore, the wave velocity corresponding to different types of waves discussed in Table 1 can be modified, as is shown in Table 3. The procedure figure of velocity correction of the first arrival wave is shown in Figure 9.
5.2. 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.
The 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. The 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. Therefore, 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 represents the other 4 sensors except 8#; 1 in the first column on the left in Table 4 represents 8#. The correction results for the first break types of events 1 and 2 are shown in Table 4. After field calibration, is set to be 3520.25 m/s.
It can be seen from Table 4 that the time difference between the signals recorded by all the sensors, except sensor 2# and 8#, is less than TLP1i, and therefore, the first arrival wave they receive is the P wave. The 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 t1i is all less than TLD1i, where the wave velocity 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.
6. 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. Thus, it can be judged that the first arrival wave may be P wave, S wave, R wave or abnormal wave. When the 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. The 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. The first arrival wave velocity is simplified to be , , 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. The 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. The field test results obtained and reported here show clearly that the method proposed in this paper is not only feasible but an effective tool.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the Ph.D. Foundation of Natural Science Foundation of Shandong Province (Grant no. ZR2019BEE019). Sun and Grattan were pleased to be supported by the Royal Academy of Engineering.