Microseismic Event Location by Considering the Influence of the Empty Area in an Excavated Tunnel

The velocity model is a key factor that affects the accuracy of microseismic event location around tunnels. In this paper, we consider the effect of the empty area on the microseismic event location and present a 3D heterogeneous velocity model for excavated tunnels. The grid-based heterogeneous velocity model can describe a 3D arbitrarily complex velocity model, where the microseismic monitoring areas are divided into many blocks. The residual between the theoretical arrival time calculated by the fast marching method (FMM) and the observed arrival time is used to identify the block with the smallest residual. Particle swarm optimization (PSO) is used to improve the location accuracy in this block. Synthetic tests show that the accuracy of the microseismic event location based on the heterogeneous velocity model was higher than that based on the single velocity model, independent of whether an arrival time error was considered. We used the heterogeneous velocity model to locate 7 blasting events and 44 microseismic events with a good waveform quality in the Qinling No. 4 tunnel of the Yinhanjiwei project from 6 June 2017 to 13 June 2017 and compared the location results of the heterogeneous-velocity model with those of the single-velocity model. The results of this case study show that the events located by the heterogeneous velocity model were concentrated around the working face, which matched the actual conditions of the project, while the events located by the single-velocity model were scattered and far from the working face.


Introduction
During the excavation of a deeply buried tunnel, rockbursts are likely to occur as the working face advances [1]. Rockbursts can have severe consequences, such as casualties, engineering delays, and economic losses [2,3]. To avoid the impact of rockbursts, a technique to provide early warnings is sought. It is determined that many microseismic events occur before rockbursts. Therefore, microseismic monitoring technology is used in many projects to provide rockburst warnings, allowing certain measures to be taken in advance of an event. At present, this technology has achieved some success in avoiding danger from underground excavation [4][5][6][7][8][9].
The accuracy of the microseismic event location (MEL) is the key to microseismic monitoring [10][11][12] because it affects the later analysis and thus the early warnings. The traditional method of MEL is based on a single-velocity model (SVM) [13,14]. This method assumes that the wave velocity in the rock mass is isotropic. However, considerable research and experimental data have shown that the location error of the SVM is very large [3,15]. For example, Jiang et al. found that the P-wave arrival of a receiver far from the source was earlier than that of a receiver near the source in microseismic monitoring of a coal mine [16]. This phenomenon indicates that the SVM is not realistic and the MEL needs to consider the path of wave propagation. Therefore, to improve the accuracy of MEL, many researchers have conducted in-depth research on velocity models [17,18]. Crosson was the first to propose the use of the least squares method to simultaneously invert the velocity and source parameters [19]. Since then, the joint inversion of the velocity and source parameter has become a hot research topic and has been widely used in engineering applications [20][21][22][23]. Konca used joint inversion to invert measured data from the Nias-Simeulue earthquake and obtained good results, emphasizing the importance of multiple data sets in seismic rupture imaging [24]. However, the joint inversion requires many measurements for the solution to be accurate and stable. Therefore, in the absence of a large amount of measurement data, the solution will be unstable due to several unknowns in the equation of the joint inversion, such as anisotropic velocity parameters, source coordinates, and origin time, which makes obtaining an accurate location (AL) difficult. Feng et al. proposed a sectional velocity model to improve the accuracy of MEL and ensure the stability of location results in an excavated tunnel [3]. In the velocity model, the sensors on the same tunnel section are treated as a single group, and the sensors on other tunnel sections are treated as other groups. It is assumed that the velocity of the wave from the microseismic source to the same set of sensors is the same, while the velocity of the wave from the microseismic source to the different sets of sensors varies. This model is particularly useful when the orientations of the strata and the tunnel are nearly perpendicular. However, in some cases, the velocity of the same group of sensors is inconsistent, such as when the orientation of the strata is at an angle to the orientation of the tunnel or the distribution of the strata is more complicated. To generate velocity models for different geological conditions, Ma et al. proposed four different equivalent velocity models [25]. However, these equivalent models cannot produce arbitrarily complex velocity models.
Therefore, establishing a complex velocity model that is consistent with the actual engineering scenario is an important factor for improving the accuracy of MEL in a tunnel. To solve this problem, Peng proposed a mesh-based velocity model [26], which can accurately generate an arbitrarily complex 3D velocity model. In this paper, based on this model, a 3D heterogeneous velocity model (HVM) is proposed for MEL around tunnels. The research methods and innovations of this paper are as follows: First, by using a grid-based modeling method that can express an arbitrary velocity model, a velocity model in the microseismic monitoring area around the excavated tunnel was established. Then, considering the influence of the empty area after excavation on the microseismic location, we proposed the use of the fast marching method (FMM) to accurately calculate the travel time of the microseismic wave from the source to the sensor. Finally, based on the precise travel time calculated using the FMM, we used a new two-step location strategy. The first step was to use block localization (BL) to determine the approximate range of the event. The second step was to use particle swarm optimization (PSO) for AL of the event within the approximate range determined in the first step. The performance of the method was analyzed with synthetic tests and the method was successfully applied to MEL in the Qinling tunnel in China.

Methodologies
In this section, we mainly describe the research methods of this article in the following three parts. First, we detail the method of establishing an arbitrary velocity model. Second, we introduce the BL step, which determines the approximate range of the source. Third, we describe the AL step, which solves the optimal location of the source using PSO.

Establishment of the Velocity Model around Tunnels
With the development of computer science, many 3D modeling software packages have been developed in different fields. These software packages use grid modeling to build complex 3D geological bodies with tetrahedrons or hexahedrons. For microseismic localization, Peng proposed a grid-based velocity model [26]. Below, we briefly introduce the construction of a grid-based velocity model for a tunnel, consisting of the following four steps.
Step 1: According to the layout of the monitoring network, the monitoring area is determined. During tunnel excavation, stress concentration occurs due to the excavation of the tunnel, which disrupts the original stress pattern. When the stress is greater than the maximum stress that the rock mass can withstand, the rock mass is destroyed, and the area near the working surface is the most vulnerable. Therefore, the area near the working face is an important monitoring area in tunnel monitoring.
Step 2: According to the geological conditions in the monitoring area, a geometric model is established. In the absence of a large amount of geological data, it is difficult to obtain a complex and accurate velocity model that is consistent with the actual project [27]. Therefore, we only consider the influence of the empty area in the tunnel after excavation on the velocity model. However, the proposed method can also be applied with a complex velocity model. The wave velocity varies among different media. In a tunnel, the empty area after excavation and the unexcavated rock mass are two distinct media. The propagation velocity of P-waves in a rock mass is generally 5500-7000 m/s. However, the propagation velocity of P-waves in air is approximately 340 m/s. In this paper, these two media are considered, and an HVM for the tunnel is established to improve the MEL accuracy.
Step 3: Monitoring area meshing is then performed. The size of the block is a key parameter that affects the accuracy of the velocity model and the travel time calculated using the FMM [28], thus affecting the accuracy of the MEL. The smaller the size of the block, the higher the location accuracy of the block. However, the computational cost increases sharply as the block size gets smaller. Therefore, the determination of the block size should consider the balance between computational accuracy and computational cost.
Step 4: A geometric model is used as a constraint to assign a velocity value to the block. In this study, the tunnel monitoring area included two different media: the empty area after the excavation and the unexcavated rock mass. A given block belonged to the empty area after excavation or to the unexcavated rock mass, and the corresponding velocity was determined after deciding which medium the block belonged to. The implementation of this method is as follows: a ray is shot from the centroid of a block, and the intersections with all of the polyhedrons are counted. If the number is odd, the point is inside the polyhedron; if it is even, the point is outside. Notably, this method is only applicable to convex polyhedra. Since the tunnel model in this study was a convex polyhedron, this method was applicable. After assigning a velocity to all the blocks, an HVM was established.
Through the above method, the schematic diagram of the establishment of a 2D arbitrary velocity model is shown in Figure 1. Figure 1a shows the microseismic monitoring area, corresponding to step 1 above. Figure 1b shows the monitoring area after meshing, corresponding to step 3 above. Figure 1c shows the geometry of a geological model with three media in the monitoring area, corresponding to step 2 above. Figure 1d shows the figure after assigning a velocity to each block using the geometry of the geological model as a constraint, corresponding to step 4 above.

Block Localization
BL is used to find the block closest to the location of the source (the target block). Microseismic localization based on arrival time theory is the most widely used method [29][30][31]. The basic idea of this method is to calculate the residual between the observed arrival time and the theoretical arrival time and identify the minimum value of the residual in the space as the optimal solution. Based on dividing the monitoring area into the uniform blocks mentioned above, the center of mass of each block is adopted as the BL parameter. Each block is a potential source block. We use (x k , y k , z k ) to represent the coordinates of the center of mass of each block, where k = 1, 2, 3, . . . M and M is the total number of blocks. (x i , y i , z i ) represents the coordinates of the i−th receiver, where i = 1, 2, 3, . . . N and N is the total number of receivers. If the k − th block is the target block, the observation arrival time of the i−th receiver is: where t obs i is the observation arrival time of the i−th receiver, t 0 is the origin time of the source, T obs ik is the travel time from the center of the k−th block to the i−th receiver, and ξ i represents the error of observation time for the i−th receiver. Similarly, the observation arrival time of the j−th receiver is expressed as: By subtracting Equation (2) from Equation (1), the difference ∆t obs ij between the observed arrival times of the i−th and j−th receivers is obtained. ξ i and ξ j are much smaller than T obs ik and T obs jk , respectively; therefore, we assume that ξ i minus ξ j is also equal to 0: The corresponding theoretical arrival time difference ∆T rt ijk is: T rt ik is the theoretical travel time from the k−th block to the i−th receiver, that is, the travel time from point (x k , y k , z k ) to point (x i , y i , z i ). Since the velocity model is heterogeneous, we obtain the theoretical arrival time via the FMM [28].
Therefore, the objective function is expressed as follows: where m is the norm (m ≥ 1). One difficulty with the L2 method for acoustic emission and microseismic source location is that the input errors often do not follow a normal distribution, as is assumed by the method [32]. Therefore, we used m = 1 in this study. According to Equation (5), the objective function value of each block is obtained, and the block with the minimum value is selected as the target block B,

Accurate Location
AL utilizes BL to find the optimal source location in the target block. The location accuracy of the block varies with the size of the block. The smaller the size of the block, the higher the location accuracy of the block. However, the computational cost increases sharply as the block size gets smaller. Therefore, in the first step of the location, we adopted an appropriate mesh size for the BL. Below, we describe the method for further location in the target block after the BL. We assume that P is any point in the target block, with coordinates of x p , y p , z p , and that the velocity in the block is constant, represented by v b . We calculate the theoretical travel time from point P to the i − th receiver using Equation (6) [26]: where ψ is the spatial domain of the targeted block and (p ix , p iy , p iz ) is the gradient vector, defined as follows: This is the final step of the source localization. Then, to solve the origin time t 0 , we use the following objective function G: We use PSO (see Appendix A for the pseudocode) to solve Equation (9) and obtain accurate source coordinates and the origin time.

Results and Discussions
In this section, we first verify the performance of proposed method through synthetic tests and then apply our proposed method to real data.

Synthetic Tests
In this section, we first build a simple tunnel model and then test the accuracy and efficiency of the HVM-based approach with the synthetic microseismic data. Finally, the noise immunity of the HVM-based approach is analyzed by adding the arrival time error, and the SVM-based approach is compared to the HVM-based approach.

Establishment of the Tunnel Engineering Model
As shown in Figure 2, the microseismic monitoring area range was a cube, along x coordinates from 0 m to 200 m, y coordinates from −30 m to 30 m, and z coordinates from −30 m to 30 m. A rectangular section of the tunnel, perpendicular to the x direction, was 5 m × 5 m. A total of six receivers were arranged on both sides of the tunnel, which are represented by the green triangles and denoted as R1, R2,..., R6. The three designed seismic sources, represented using red circles, are denoted S1, S2, and S3. The specific parameters of the receivers and seismic sources are shown in Table 1. The wave velocity in the monitoring area was 5000 m/s. However, since the medium in the tunnel after excavation was air, and the energy of the waves in air is greatly attenuated, the propagation velocity in the air was much lower than the propagation velocity in the rock. Therefore, the velocity of the empty area after the tunnel excavation was set to 340 m/s.
The propagation of waves follow the principle of minimum travel time. As the wave travels slowly in the tunnel after excavation, it bypasses the empty area and travels in the rock. In this example, the rock is homogeneous, and the wave velocity is the same throughout the rock; therefore, the minimum path from the source to the receiver is the shortest path. According to the knowledge of spatial analytic geometry, we can obtain the minimum travel path, namely, the theoretical path, which is represented by a solid black line in Figure 3. The specific calculation method of the theoretical path is as follows. Figure 4a is the tunnel model. The faces ADHE and BFGC are the two sides of the tunnel. R is the receiver that is very close to the face ADHE, and S is the source near the face BFGC. R and S are the vertical projections of R and S on the faces ADHE and BFGC, respectively. We first find the shortest distance from R to S . We expand the tunnel model vertically to obtain Figure 4b. According to the principle that the line segment between two points is the shortest, we obtain the shortest path from S to R , which intersects the line segments AD and BC at points M and N. We then utilize plane DNMR to obtain the shortest path from S to R (green line in Figure 4a). According to the theoretical path, the theoretical arrival time of each source to each receiver can be calculated, as shown in Table 2.  Based on the abovementioned tunnel engineering model, we used an HVM to locate the three sources. First, we meshed the monitoring area with a grid size of 0.5 m × 0.5 m × 0.5 m, thus obtaining 400 × 120 × 120 unit blocks. The tunnel geometry model was used as a constraint to assign velocity values to the unit block. The velocity of the unit block inside the tunnel was 340 m/s, and that of the unit block outside the tunnel was 5000 m/s. BL was carried out using the theoretical arrival time in Table 2 as input parameters. Each block was a potential source, and the residual arrival time of each block was calculated according to Equation (5). The theoretical travel time for the waves between the sources and receivers was calculated using the FMM, as shown in Figure 4. The red lines in the figure represent the calculated path from the sources to the receivers. We used the grid search method to assign the block with the minimum value of Equation (5) as the target block, and the center coordinates of the target block represent the BL results, as shown in Table 3.
AL was then performed in the target block according to Equations (6)- (9). The minimum value of Equation (9) was found using PSO, and the AL results are shown in Table 3. Figure 5 shows that the AL results of the three sources were very close to the theoretical position. The PSO iteration parameters were as follows: The maximum number of iterations was 2000, and the number of seeds was 80. The acceleration parameters of the algorithm were 2 and 2, which affected the local and global optimal values, respectively. The weighted values for the initial and convergence moments were 0.9 and 0.4, respectively. The threshold of the termination algorithm was 1 × 10 −25 , and when the minimum value of the target function was less than this value, the algorithm stopped. The change in the value of the objective function in the PSO iteration is shown in Figure 6, and the three sources converged after 25 iterations.   A description of the computational efficiency is as follows. First, we used the FMM to calculate the travel time of the waves from each receiver to all grids and saved these travel times in a database. As long as the receiver position and velocity model were not changed, each subsequent microseismic event was located via the travel time database, such that the FMM solution needed to be calculated only once. In this paper, the FMM code was based on C++ programming, and the other code was based on MATLAB programming. In this case, it took 51 s to construct the travel time database using the FMM. The BL and AL computation times for the three events are shown in Table 3. After obtaining the travel time database, locating an event took approximately 17 s. All the above programs ran on a 3.6 GHz Intel Core i9-9900k CPU.
As shown in Figure 7, we compared the BL results of the three sources with the AL results. The teal columns in Figure 7 represent the errors of the BL results, and the pink columns represent the errors of the AL results. The location results include the X, Y, and Z coordinate errors, spatial error, origin time error, and minimum value of the target function. The AL errors of S1 and S2 were larger than the BL errors in the Y direction, and the AL errors of S3 were larger than the BL errors in the Z direction. We calculated the spatial location error using Equation (10): where ∆SR represents the spatial location error. ∆XR, ∆YR, and ∆ZR represent the location errors of the X, Y, and Z coordinates, respectively. X c , Y c , and Z c represent the location results of the X, Y, and Z coordinates, respectively. X t , Z t , and Z t represent the theoretical values of the X, Y, and Z coordinates, respectively.
However, the AL results of the three sources were clearly smaller than the BL results in terms of the spatial error, which indicates that the AL results were closer to the theoretical position. In addition, in terms of the error of the origin time and the minimum value of the objective function, the AL results were also better than the BL results.

Comparison and Analysis
Due to the influence of noise in the tunnel, there was a certain error in the arrival picking in actual projects. To further test the practicality of the HVM-based method, we added a certain noise in the theoretical arrival time (see Table 2). We designed three contrasting experiments. First, the theoretical arrival time was used as the location input parameter, and an SVM was adopted to locate the three sources, which was expressed using the SVM. Second, the noisy arrival time was used as the location input parameter, and an HVM was adopted to locate the three sources, which was represented using the HVM(N), where N is stand for noisy arrival time. Third, the noisy arrival time is used as the location input parameter, and an SVM was adopted to locate the three sources, which is represented using the SVM(N). The location results of the three experiments are shown in Table 4, which clearly shows that the origin time obtained by these three experiments was very close to the theoretical origin time. However, the location results of the SVM and SVM(N) were significantly different from those of the HVM(N) in the Y direction. We used these three experiments to compare the errors of the location results of the HVM and SVM in detail, as shown in Figure 8. The origin time error of these four experiments was within 0.6 ms, and the accuracy was very high, which was not used as a criterion. In the case where no noise was added, the location error of the SVM of the three sources in the X direction was smaller than that of the HVM, but the location error of the HVM in the Y direction was smaller than that in the SVM. The location errors of the SVM in the Z direction of S1 and S3 were smaller than those of the HVM, but the opposite was true for S2. Considering the location errors in the X, Y, and Z directions, it was impossible to determine which method had a higher location accuracy. The spatial location error was a comprehensive error in the integrated X, Y, and Z directions. Therefore, the spatial location error calculated using Equation (10) was used as the criterion. Figure 8 clearly shows that the spatial location error of the HVM was smaller than the spatial location error of the SVM. For the average spatial location errors of the three sources, HVM (2.06 m) < SVM (5.51 m). In addition, from the minimum value of the objective function, the value of the HVM-based MEL was smaller. In the case with the added noise, the spatial location error and the residual of the HVM-based method were smaller than those of the SVM-based method. Regarding the average spatial location error of the three sources, HVM (4.95 m) < SVM (6.23 m). In summary, it can be concluded that the location accuracy of the HVM was higher than that of the SVM. Below, we use the proposed method to analyze real data.

Application to Real Data
The Qinling No. 4 tunnel of the Yinhanjiwei project is located south of the Qin Mountains in southern Shaanxi Province, China. The length of this shaft is 5820 m, with a section 6.5 m high and 6.7 m wide. The maximum slope is 11.96%, and the maximum depth is 1600 m. Drilling and blasting are used in this project, resulting in frequent rockbursts that may both damage infrastructure and injure people. To provide safety guidance, a microseismic system was used for 24 h of continuous monitoring. Four accelerometers with a sensitivity of 10 V/g were embedded along two sides of the tunnel. The spatial arrangement of the receiver is shown in Figure 9a. The sampling frequency was set to 10 kHz.
The working face of the actual project is shown in Figure 9d. Due to the poor lighting inside the tunnel, to distinguish the microseismic hole, the sensor installation position was marked with red paint for visibility, as shown in Figure 9c. Some ejected fragments were found at the top of the working face on 15 June 2017, as shown in Figure 9b. According to engineering experience, these fragments formed due to the destruction of the roof rock mass caused by the excavation of the tunnel. Below, we verify this through microseismic monitoring.
During the tunneling process, the microseismic monitoring system monitored a large number of events. The microseismic event and the blasting event could be clearly distinguished by the waveform, as shown in Figure 10. We selected 51 events with better waveforms during the period from 6 June 2017 to 13 June 2017, including 7 blast events and 44 microseismic events. We verified the location effect in two ways: (1) due to the known blasting position of the working face, we verified the location effect based on 7 blasting events; and (2) we verified the location effect based on the spatial relationship between the location results of 44 microseismic events and the working surface.  We used the HVM to locate these 51 events. The monitoring area covered the x coordinates from 3,727,271 m to 3,727,516 m, y coordinates from 502,564 m to 502,761 m, and z coordinates from 558 m to 598 m. The block size was 0.5 m × 0.5 m × 0.5 m, so there were a total of 490 × 394 × 80 unit blocks. There are many ways to calculate the propagation velocity of waves in rock masses. For example, Wang et al. optimized the seismic wave velocity in the deep mining area of a coal mine by using a combination method, residual error optimization method, location error optimization method, location residual optimization method, and combined inversion method [33]. The acquisition of the velocity model was not the focus of this paper. In this application, the wave propagated at a velocity of approximately 6000 m/s in the rock mass. Therefore, we set the velocity of the unexcavated rock mass to 6000 m/s and the velocity of the empty area after excavation to 340 m/s, as shown in Figure 11. The coordinates of the receivers and the arrival time of the seven blasting events are shown in Table 5. The arrival time of the 44 microseismic events is shown in Appendix B.  The location results are shown in Figure 12. In this figure, red indicates the location result of the HVM, and blue indicates the location result of the SVM; five-pointed stars indicate the blast events, and circles indicate the microseismic events. The results of the seven blasting events are shown in Table 6. The results of the 44 microseismic events are shown in Appendix C.
First, we analyzed the results of the HVM. From the location results of the blasting events, it can be clearly seen that the seven blasting events occurred in the vicinity of the working face, which was consistent with the actual engineering excavation. The location results of the microseismic events showed that 44 microseismic events were concentrated around the blasting events. The blasting of the tunnel face caused damage to the surrounding rock mass. The microseismic signal from the rock mass damage could be received by the sensor; therefore, in theory, most of the events occurred near the working face. The location results of the 44 microseismic events were consistent with this theoretical derivation.
Then, we analyzed the location results of the SVM. As shown in Figure 12, the location results of the blasting events suggested that the SVM results were very scattered. Four of the blasting events were located behind the propulsion direction of the working face, and two blasting events were located approximately 85 m in front of the propulsion direction of the working face. From the location results of the microseismic events, the spatial distribution of the 44 microseismic events was consistent with that of the seven blasting events. The spatial distribution was very scattered, irregular, and not concentrated near the working surface.
In summary, we can conclude that the HVM had a high location accuracy and good effect, with a location accuracy that was much higher than that of the SVM.  What caused the location results of the SVM-based method to be so poor? We compared these results with the HVM-based MEL and found that the velocity model error decreased the location accuracy. After excavation, the tunnel was filled with air, and the propagation of the microseismic signal was greatly affected. The wave did not pass directly through the empty zone to the receiver, but rather bypassed the empty zone and reached the receiver by travelling through the unexcavated rock mass. Therefore, in the MEL in the tunnel, we should have fully considered the impact of the empty zone and used the HVM for location.
The HVM-based method proposed in this paper had a high precision and was suitable for MEL during tunnel excavation. The location accuracy of the HVM varied with the size of the mesh. The larger the mesh size was, the lower the accuracy. Conversely, the smaller the mesh size was, the higher the accuracy. However, the smaller the grid was, the lower the computational efficiency. Therefore, in practical engineering applications, we should consider both accuracy and efficiency when determining the appropriate grid size.

Conclusions
We adopted a grid-based modeling method that can express arbitrary velocity models and establish a 3D HVM in the microseismic monitoring area of an excavated tunnel. For the velocity model the tunnel, we fully considered the influence of the empty area on the location result and set the velocity in the tunnel after excavation to 340 m/s. The approximate range of the source was determined using BL, and then the exact position of the source was determined using the AL.
The synthetic tests showed that the location accuracy of the proposed HVM-based method was higher than that of the traditional SVM-based method and that the proposed method had certain anti-interference characteristics. The average spatial location error of the HVM was less than that of the SVM: HVM (2.06 m) < SVM (5.51 m) without adding noise. With added noise, the location accuracy of the HVM-based method was also higher than that of the SVM-based method: HVM (4.95 m) < SVM (6.23 m). Finally, the HVM-based method was applied to the monitoring of the Qinling No. 4 tunnel of the Yinhanjiwei project. The results showed that the event locations of the HVM were concentrated near the working surface, which was in line with observations made during engineering practice. However, the event distribution of the SVM-based method was very scattered and irregular. Therefore, the empty area created by the excavation had a great influence on the microseismic event location around the tunnel. The HVM improved the location accuracy of microseisms around the tunnel and has practical research significance.
Author Contributions: Y.J. and P.P. participated in the data analysis, helped to design the study, drafted the manuscript, and carried out the statistical analyses; Z.H. collected the field data; L.W. conceived the study, lead the design of the study, coordinated the study, and helped draft the manuscript. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A The PSO Pseudocode
Procedure PSO For each particle i Initialize velocity Vi and position Xi for particle i Evaluate particle i and set pBesti = Xi End for gBest = min pBesti While not stop For i = 1 to N Update the velocity and position of particle i Evaluate particle i If f it(Xi) < f it(pBesti) pBesti = Xi If f it(pBesti) < f it(gbest) gBest = pBesti End for End while Print End procedure Appendix B