Microseismic Source Location Based on Improved Artificial Bee Colony Algorithm: Performance Analysis and Case Study

Highly accurate microseismic (MS) localization is the basis for rock damage assessment and disaster warning. The engineering background noise mixed in the MS signal seriously affects the subsequent analysis of the MS signal. A noise reduction method of singular spectral analysis-complementary ensemble empirical mode decomposition-wavelet threshold (CEEMD-WT-SSA) is proposed. The complementary ensemble empirical mode decomposition (CEEMD), complementary ensemble empirical mode decomposition-wavelet threshold (CEEMD-WT) and the proposed method are respectively used for denoise the noisy Ricker wavelet. The signal-to-noise ratio of the denoised signal by the proposed method is 56.77% and 37.88% higher than that of CEEMD and CEEMD-WT methods. Moreover, an adaptive artificial bee colony (ABC) algorithm is applied for MS source location. The time to quantile difference is introduced as the objective function. The blast positioning test results prove that the proposed method improves the positioning accuracy by 44.12% and 47.64% particle swarm optimization (PSO) algorithm and simulated annealing particle swarm optimization (SA-PSO) algorithm, respectively. The MS positions of underground caverns reveal that the calculated cluster of MS events using the adaptive ABC algorithm are more concentrated at the structural plane and appearance deformation failure location and in good agreement with field survey and routine monitoring data.


Introduction
The underground space is in a great demand with the rapid development of economy in China. Myriads of deep underground projects are being built or to be built. Affected by geological structures and high ground stresses [1], deformation, and failure of surrounding rock (slab, collapse, rock burst) [2][3][4] induced by excavation unloading of deep underground caverns are frequently encountered, which seriously threatens people's lives and the safety of production activities.
As an effective monitoring method, the microseismic (MS) monitoring technique could collect elastic waves released by fracturing rocks using embedding sensors in engineering rock mass, and invert MS source information to obtain rock fracture characteristics and delineate surrounding rock damage areas [5]. In addition, early warning for upcoming geological disasters such as rock bursts can also be carried out. At present, MS monitoring, a three-dimensional realtime monitoring technology, has been widely applied to deep underground engineering such as mines, tunnels, and underground powerhouse of hydropower station [6][7][8].
Denoising of MS signal is the basis of MS source location. At present, the mainstream denoising method in China is the improvement of EMD denoising method and wavelet denoising method. Li et al. [9] applied EMD method to rock blasting and impact breaking signals and achieved good results. Han et al. [10] proposed an adaptive threshold EEMD method and applied it to de-noising hydraulic fracturing MS signals. Tang et al. [11] proposed an improved MS denoising method based on different wavelet coefficients of useful signal and noise components, it can retain the MS information well and has a good denoising effect. These methods are of great significance to improve the quality and analysis of MS data.
As a core content of the MS monitoring technology, a high precision of MS location is of great significance for stability analysis [12,13]. Since Geiger [14] firstly proposed the absolute location method, a lot of research on the factors affecting the location accuracy of has been conducted. In recent years, many scholars tried to explore the influence of MS location velocity model on MS location. Lurka et al. [15] used inversion principle and evolutionary algorithm to build a three-dimensional velocity model and applied it to deep mines in the west of the United States. Compared with the fixed velocity model, the three-dimensional velocity model built has improved the positioning accuracy. Based on the assumption of horizontal layered medium condition, Cong et al. [16] proposed a new velocity inversion method, and use DIRECT fast search algorithm to solve the objective function, the results show that this method has good adaptability to different monitoring systems. Wamriew et al. [17] proposed a deep learning method for real-time reconstruction of velocity models from acquired MS signals. Compared with traditional velocity models, this method has higher positioning accuracy. Scholars have also carried out a lot of research on P-wave arrival data. Peng et al. [18] used Log-Cosh function to eliminate the arrival time of remote sensors to improve the location accuracy of mine MS. Rui et al. [19] have proposed a clustering detection method to eliminate the AE positioning method of abnormal arrival time, and have verified that the positioning results of abnormal arrival time elimination have good robustness through simulation experiments. In addition, some scholars have also studied the influence of MS network layout on positioning accuracy. Gong et al. [20] have established a solution model suitable for large-scale network combined planning for MS monitoring system, aiming at the optimization problem of large-scale network, effectively solving the layout and optimization problem of mine seismic network. Li et al. [21] proposed sensor network optimization principles to improve positioning accuracy by studying sensor spatial layout. WardNeale et al. [22] proposed an improved multi-array source inversion positioning method by optimizing station layout. However, due to the complexity of underground engineering geological conditions and construction environment, the precise location of MS source (or rock burst) still faces many challenges. In addition to sensor array layout, velocity model and arrival pick accuracy, sensor position measurement error [23] and optimization algorithm [24] are also factors that affect the location accuracy of seismic source. Among them, sensor array layout, velocity model and time pick accuracy are always the focus of research.
In this study, an improved noise reduction algorithm is introduced to denoise the MS signals. After the MS signals are denoised, the signals are picked up accurately on time by the time-length ratio method -Akaike information criterion (STA/LTA-AIC) method, to achieve the precise pick up of MS signals under the influence of noise. Finally, the adaptive ABC algorithm is used to locate the processed MS signals, the location results can be verified according to the field routine monitoring data and damage results.

Basic Theory
At present, Fourier transform [25], empirical mode decomposition (EMD) and wavelet transform are the main denoising methods commonly used in the field of MS signal denoising. Due to the nonlinear and non-stationary characteristics of MS signals, the fast Fourier transform is not ideal for denoising such signals. The wavelet transform has good local property and high accuracy in the timefrequency domain, but it is more dependent on personal subjective will to select the basis function. Once the selection of the basis function and the noise reduction threshold is different, there will be a big difference [26]. The EMD has advantages in processing non-stationary signals, but it is prone to mode aliasing [27]. In recent years, the improved ensemble empirical mode decomposition (EEMD) based on EMD has made up for the mode aliasing problem existing in EMD to a certain extent. However, intrinsic mode function (IMF) after lumped average cannot accurately maintain the definition of IMF component defined by EMD decomposition, and the signal white noise after reconstruction cannot be completely offset. Therefore, the complementary ensemble empirical mode decomposition (CEEMD) is proposed to solve this problem [28].
Based on the above analysis, the CEEMD method is adopted to decompose noisy MS signals and obtain a series of IMF components from high frequency to low frequency. In addition, the sample entropy is introduced to determine the boundary point between high frequency component and low frequency component by calculating the sample entropy of each IMF components. The wavelet threshold denoising method is used to denoise the high frequency IMF components, and then the high frequency IMF components after denoising and the low frequency components are reconstructed to obtain the denoised MS signal, and then the reconstructed signal is filtered twice by combining with singular spectrum analysis (SSA).

Denoising Principle based on CEEMD-WT-SSA
Sample entropy is an index represents the complexity of time series and its size represents the noise content in the signal. The specific algorithm of sample entropy can be referred to [29], which will not be detailed in this section.
The wavelet threshold denoising method mainly includes three basic steps: (1) Transforming the original signal by wavelet transform. The high frequency wavelet coefficients and low frequency wavelet coefficients are obtained by reasonable wavelet basis function and decomposition layer.
(2) Processing threshold for decomposed wavelet coefficients. By setting thresholds at different scales, the high-frequency wavelet coefficients are treated as thresholds to retain the low-frequency components.
(3) Reconstructing the denoised signal. The processed high frequency and low frequency wavelet coefficients are inverted to obtain the denoised signal.
Wavelet threshold function is generally divided into hard threshold and soft threshold. The hard threshold function is discontinuous, and it is easy to have spikes after denoising. The signal reconstructed by soft threshold function has a good smooth line, but the error is large. The selection of threshold is a key step in wavelet denoising. The commonly used threshold estimation methods include fixed threshold estimation, unbiased likelihood estimation, heuristic estimation and extreme threshold estimation. In this study, hard threshold is selected for denoising, specifically as follows: Where ω is the wavelet coefficient, λ is the threshold value, ζ is the noise standard deviation ， and N is the length of signals.
The CEEMD is based on the EMD decomposition, which eliminates modal aliasing by adding paired, opposite numbers of white noise. The specific steps are as follows: (1) Add I paired white noise (2) The EMD decomposition is used to decompose each group of noise-added signals, and finally the im IMF  and im IMF  components can be decomposed.
(3) The final decomposition result is obtained by averaging the I pieces components. For the decomposed IMF component, calculate the sample entropy value of each component, and finally use the average sample entropy value as the dividing point between the high-frequency component and the low-frequency component to perform wavelet threshold denoising on the high-frequency component. Through simulation experiments, it is concluded that using db4 wavelet basis function for wavelet threshold denoising is better. Reconstruct the highfrequency component after denoising and the low-frequency component without denoising to obtain the denoised signal X(t).

Singular spectrum analysis
The HanKel trajectory matrix of the denoised time series is calculated [30]. Firstly, it is necessary to select an appropriate window length L to arrange the original signals to obtain the trajectory matrix. Generally, it is taken L Singular value decomposition is performed on the trajectory matrix, where the trajectory matrix can be rewritten as Eq. (7).
Where d rand( X )  is the intrinsic component or mode number with nonzero eigenvalues, i λ is the eigenvalue, and each eigenvalue represents the contribution rate of the eigenvector after singular value decomposition to the original signal energy.
Eq. (9) is called the contribution rate of with a length of d N , the time series G is the target signal extracted from the original signal X. The specific reconstruction method is as follows:  , represents the elements in the m-th row and n-th column of the matrix * I X , and each element in G can be calculated using the following equation.
The basic flow chart of the proposed noise removal method is shown in Fig. 1.

Simulation analysis
To verify the effectiveness of the proposed denoising method, Ricker wave of simulated earthquake are introduced for simulation experiment, and the original expression of Ricker wave is as follows: represents amplitude and m f represents center frequency. In this study, m f is taken as 30 Hz and the sampling frequency is 1000 Hz. Fig. 2a and b show the waveform and spectrum diagrams of Ricker wavelet, respectively. Add 22dB Gaussian white noise to the Ricker wavelet to obtain the waveform after noise addition as shown in Fig. 2(c) and (d). In order to illustrate the advantages of the noise reduction method, CEEMD and CEEMD-WT are used to decompose Ricker wavelet with white noise. The results are shown in Fig.  3(a)-(d) respectively. CEEMD removes most of the high frequency noise, but some noise remains.
However, the denoising method based on CEEMD-WT retains more useful information of the original signal and further reduces the noise interference. As can be seen from Fig. 3(e) and (f), the proposed CEEMD-WT-SSA denoising method is more effective, and the waveform after denoising is smoother. To illustrate the effect of the noise reduction method in this study, on the one hand, waveform diagram and spectrum diagram are combined to make a qualitative analysis of the noise reduction effect. On the other hand, the SNR Eq. (13) is introduced to quantitatively describe the noise reduction effect.
Where x( t ) is the signal not polluted by noise; () xt is the signal containing noise; N is the signal sampling point. In the simulation experiment, when the signal is known to be clean without noise pollution, the larger the SNR is, the less noise the signal contains.
Eq. (13) is used to calculate the SNR before and after noise reduction of the above three methods, which was recorded in Table 1. Through the analysis of the calculated results, can be seen that among the three methods, the proposed denoisinh method has a higher SNR after noise reduction than the other two methods, while the CEEMD-WT method has a SNR after noise reduction than the single CEEMD. The calculated results are in agreement with the quantitative analysis of spectral graph, which indicates that the proposed method is more effective. According to Eq. (14), the denoising accuracy of CEEMD-WT-SSA is 56.77% and 37.88% higher than that of CCEMD and CEEMD-WT methods, respectively.
Where Lr stands for improving accuracy; 0 SNR represents the SNR value after CEEMD-WT-SSA denoising, and 1 SNR represents the SNR value after other two methods.

P-wave arrival time picking
In the preprocessing of MS data, accurate arrival picking of waveform is a key step, and the accuracy of picking up is directly related to the positioning accuracy of MS events and subsequent results analysis. Although the arrival time of manually picking up waves is relatively stable, it is easy to introduce human error and has low picking efficiency. Therefore, many scholars at home and abroad have proposed many methods for automatically picking up the arrival time of microseisms [31,32]. proposed STA/LTA to identify signals by using the variation of amplitude energy ratio in the time-domain signal time-short window. This method is simple to calculate and is one of the most commonly used methods at present, but it has a high misjudgment rate. Akaike information criterion (AIC) algorithm [33,34] is one of the most widely used methods, featuring high computational accuracy and good noise robustness. However, the selection of time window is artificially subjective, which may cause large errors. In this study, STA/LTA-AIC algorithm is used to extract the MS signal after noise reduction pre-processing. The idea of this algorithm is as follows: When the MS signal arrives, there will be a sudden change in the time-short mean ratio STA/LTA. When the ratio is greater than a certain threshold δ , the sudden change point is usually regarded as the arrival point of the MS signal [35]. The specific calculation process of STA/LTA is shown in Eqs.
Where n refers to the sampling point, sw is the shorttime window length, lw is the long-time window length, and CF( i ) is the characteristic function value of MS signal at the moment. The characteristic function selected is shown in Eq. (18) below.
AIC was originally used to measure the goodness of statistical model fitting. In the field of MS phase picking, the local minimum value of the model is calculated to pick up the seismic phase. For MS signals x( t ) with n sampling points, the AIC function is calculated as follows:

Construction of MS location objective function
At present, the MS location of underground engineering mainly adopts the more easily identified P-wave seismic phase for positioning, which has less error and is easy to pick up when compared with other seismic phases. Due to the complexity of rock mass media, it is assumed that the rock mass in the MS monitoring area is isotropic homogeneous, that is, the P-wave velocity in the monitoring area remains unchanged. When the time of P wave detected by station i of a MS source O( x ,y ,z ) 0 0 0 is t , then the propagation time of P wave to the i-th station can be expressed as: Where v represents the p-wave velocity, The fixed wave velocity has a large error for complex rock mass media, so the P wave velocity is taken as an unknown quantity v 0 and the difference method is adopted. Then the calculate time c i t is as follows: The arrival time difference between two different sensors i and j is The unknown source seismic time t 0 is eliminated through the difference, and the smaller the time difference is, the less error about MS location, as shown in Eq. (25) below.
The difference method only deals with the difference of time between all stations and a random station. The introduction of the first and third quartiles [36] greatly enhanced the robustness and global of the objective Eq. (25). Last to final MS source objective Eq. (26). The derivation process of Eq. (26) will not be discussed in this paper. Please refer to literature [36] for details.

Improved artificial bee colony algorithm
ABC algorithm [37], proposed by Karaboga in 2005, is a swarm intelligent optimization algorithm inspired by the honey collecting mechanism of bees. The ABC algorithm is used to search for the minimum point of the objective function of MS positioning in Eq. (26). The corresponding position of the minimum point and the focal position of the MS are shown in the following steps. The basic flow chart sees Fig. 5. Step 1: Generate 2N locations (representing potential source locations) in the MS space according to Eq. (27), and calculate the fitness value of the objective function according to Eq. (28): Step 2: Employed bees to search for new nectar sources near nectar sources according to Eq. (29).
Where ij x is the j-th dimension position of nectar source i; kj x is the k-th dimension position of nectar source j, which is not equal to i. ij  is a random number between [-1,1].
Step 3: Compare the nectar sources before and after the search. If the nectar source after the search is better than that before the search, replace the nectar source before the search. The follower bees selected honey sources according to the release information of the employed bees at the honey source in the roulette mode, and updated the selected honey source information according to Eq. (29).
Step 4: Set the parameter limit . If some nectar source remains unchanged after the limit cycles, give up the nectar source. Then the corresponding employed bees are transformed into scout bees, and new honey sources are randomly generated according to Eq. (27).
In the traditional ABC algorithm, the high randomness of roulette selection strategy is easy to lead the algorithm to fall into the local optimal solution. The honey source search equation has a strong random performance. Although it can enhance the search ability, the convergence rate of the algorithm is too slow. To address these problems, the exponential ranking method is introduced to select the nectar source, and proposed a honey source search equation with fast convergence strategy. These two strategies are described in detail below.
The exponential ranking method is a selection strategy based on user-defined probabilities. Firstly, honey sources are ranked according to their fitness. The honey source with the worst fitness is assigned a number 1, and the best honey source is assigned a number N. For the honey source with the number i, the probability of being selected as a follower bee is i p , which is calculated by Eq. (30). In the formula, c is a fixed base, and 0<c<1 needs to be satisfied.
Based on the literature [38], Eq. (29) was improved through many tests, and an adaptive fast convergence Eq. (31) was obtained.
is the global optimal solution under the current loop, t is the current number of cycles and maxcylce is the maximum number of cycles.

Algorithm verification analysis
To verify the positioning accuracy of the MS location method based on the improved artificial bee colony (IABC) algorithm, a three-dimensional model with a size of 200m×200m×200m and a wave velocity of 3800m/s of rock mass is established, as shown in Fig. 6. Eight monitoring stations and five random MS sources were set up in the simulation experiment. The coordinates of the station see Table 2, and the coordinates of the source are shown in Table 3. The time of earthquake onset is 10 seconds. It is assumed that the range of each search dimension of the source is 0~200m. For different optimization algorithms, the positioning accuracy error simulation experiment was carried out, and the results were compared in Fig. 8. The relative error in the figure is the difference between the calculated distance from the source position to the sensor and the actual distance from the source to the sensor. There are three hyperparameters in the IABC algorithm: the number of honey sources N, the total number of iterations Maxcycle, and the upper limit of mining limit . In this study, we set the number of honey sources N=200, the total number of iterations Maxcycle=150, and the upper limit of mining limit =20.  As can be seen from the fitness value changing with the number of iterations in Fig. 7, under the same number of iterations, the convergence speed of PSO algorithm is the slowest, that of SA-PSO algorithm is the second, and that of IABC algorithm is the fastest. Among the three optimization algorithms, the IABC algorithm has the lowest fitness value, indicating that it can find the global optimal value quickly. Fig.  8 clearly shows the positioning errors of the three optimization algorithms for the five focal points in the simulation experiment. It can be seen from the figure that the three optimization algorithms can meet the needs of MS positioning accuracy, but the positioning accuracy of IABC algorithm is higher. For the five simulated focal points, the positioning accuracy errors of IABC algorithm can be kept within 6m, showing good robustness. The positioning error accuracy of SA-PSO algorithm and PSO algorithm is kept within 8m to 10m, which has a certain gap compared with the IABC algorithm.

Project description
The improved optimization algorithm in Section 4 is applied to the location of blasting events and MS events during the excavation of underground powerhouse in a hydropower station in Western China. The underground powerhouse is composed of main powerhouse size of 183.5m×25.8m×65.25m (length × width × height), transformer chamber of 117.2m×16.5m×30.8m (length × width × height), tailrace surge chamber, etc. The lithology of the underground powerhouse section is revealed to be medium-thick sandstone, extremely thin slab and phyllite in the T3Z2 (3) ~ T3Z2 (5) formation, with the maximum vertical burial depth of 250m. The structural plane of the underground workshop area is shown in Fig. 9. The background noise caused by the construction disturbance on multiple working surfaces of underground powerhouse leads to low signal pickup accuracy. The robustness and global performance of the traditional MS location optimization algorithm are poor, and it has fallen into the local optimal solution in the global search, and the convergence rate is slow. Based on previous studies, this paper discusses the location of MS source in the surrounding rock of deep underground cavern by ABC algorithm. A series of results obtained can provide a new idea for MS location of underground powerhouse.

Construction of the MS monitoring system
In order to monitor MS activities in rock mass during excavation disturbance, a high-precision MS monitoring system produced by Engineering Seismology Group (ESG) of Canada was installed in an underground chamber. The system is mainly composed of uniaxial acceleration sensors, Paladin digital signal acquisition instrument, Hyperion digital signal processing system. In this MS monitoring system, a total of 6 acceleration sensors are arranged, and the monitoring mainly covers the side wall between the transformer chamber and the tailrace surge chamber. The sensor layout in this area (see Fig. 9), and Table 4 illustrates the sensor coordinates. Based on the acoustic test of surrounding rock and the adjustment of blasting test parameters, the P-wave velocity is finally determined to be 3831m/s. Six sensors are used to carry out real-time monitoring of the target monitoring area. When the rock breaks and releases an elastic wave, the elastic wave signal will be received by the sensors and converted into an electrical signal. The electrical signal is transmitted via cable to the Paladin acquisition substation and converted into a digital signal (A/D) and transmitted via fiber optic cable to the host computer equipped with Hyperion signal processing.

Location test of blasting events
In order to further verify the effectiveness and accuracy of the proposed MS location algorithm in practical engineering, blasting events during excavation are used to verify the positioning accuracy. The coordinates of the two blasting events were respectively (507299.93, 3495825.32, 2167.6) and (507305.73, 3495801.33, 2161.5). Fig. 10 shows the waveform of blasting event collected and the result of blasting signal picked up by STA/LTA-AIC method in Section 3. IABC, SA-PSO and PSO were used to locate the blasting event, and each method was calculated 25 times and averaged. The location results of three different optimization algorithms for blasting events are presented in Table 4. It can be seen from    In this subsection, we use the proposed signal noise reduction method to denoise the noisy signals (see Fig. 11), and then use the STA/LTA-AIC algorithm to pick up more accurate P-wave arrival times, and finally use the adaptive ABC algorithm to localize the MS signals, and the localization results are shown in Fig. 12(c) and (d). From the comparison of the figures, it can be indicated that the MS events localized by the adaptive ABC algorithm are concentrated in the range of 0+063-0+028m in the middle partition wall between the transformer chamber and the tailrace surge room, and the localization results are more concentrated and there are no abnormal localization results. Fig. 13 shows the multi-point displacement meter monitoring results at the MS event gathering area, from which it can be verified that the deformation appears to be in a fluctuating state to different degrees during June 15-July 15, 2022 by the excavation and blasting. From the construction situation, the tailrace surge chamber is mainly for the second layer of support construction, and the transformer chamber is being excavated for the third layer of blasting. The field investigation found that there was a crack in the initial support concrete and rock fissures revealed by excavation near 0+43m on the downstream side of the transformer chamber. According to the existing geological exploration results, the Fault 20 runs through the transformer chamber and the tailrace surge chamber middle partition wall area, combined with the construction conditions, routine monitoring, site investigation results, the blasting excavation construction, the rock stress in the middle partition wall area again induced self-adjustment, the internal rupture development state intensified, the rupture damage deepened. Fault 20 at the primary rupture surface disturbed by blasting excavation in the extension process the secondary rupture surfaces of different sizes are generated in the extension process after the blasting excavation, and the deterioration of the rock body is deepened, and MS events are gathered near the Fault 20.

Conclusions
Accurate delineation of surrounding rock damage area is of great significance to reduce economic losses and casualties. In this study, an improved CEEMD-WT-SSA method based on sample entropy is proposed to denoise the MS signals, then the STA/LTA-AIC method is used to accurately extract the denoised signals. Finally, the IABC algorithm is introduced to locate the events, and the results are verified in an engineering case. The conclusions are as follows: (1) A high-precision MS monitoring system is successfully constructed in a hydropower station underground cavern group. A CEEMD noise reduction method is proposed by combining wavelet threshold and singular spectrum analysis. Through qualitative and quantitative comparisons with CEEMD and CEEMD-WT, it is demonstrated that the CEEMD-WT-SSA method can remove the background noise more effectively and retain the transient non-stationary characteristics of the MS signal to the maximum extent.
(2) In this study, an objective function based on the difference of P-wave to time-division of stations is introduced and combined with the IABC algorithm for more accurate localization of MS events. In the simulation experiments, the MS localization algorithm proposed in this study has faster convergence speed compared with PSO algorithm and PSO-SA algorithm, and the localization error can be controlled within 6 m. It can be verified by blast data that for blast events M1 and M2, the localization accuracy of the localization algorithm proposed in this study is improved by 44.13% and 53.3% respectively compared with PSO algorithm; compared with the localization accuracy of the PSO-SA algorithm is improved by 47.64% and 42.9%, respectively.
(3) The conventional monitoring data shows that the rock mass near 0+40m downstream of the transformer chamber is significantly affected by excavation unloading. The cracking of the initial sprayed concrete on site and the exposure of excavation further confirm that the MS monitoring system constructed in this paper can effectively delineate the damage area of the surrounding rock. The MS events located by the MS monitoring system and the MS events located by the IABC algorithm have essentially the same orientation on elevation as the gently dipping Fault 20. It is analyzed that the cluster in this area is the result of the joint action of excavation blasting disturbance and Fault 20 weak structural plane.