The underground explosion point high precision measurement method based on multidimensional information fusion of vibration sensors

Abstract

Purpose

The purpose of this paper is to solve the problem that the location of the initiation point cannot be measured accurately in the shallow underground space, this paper proposes a method, which is based on fusion of multidimensional vibration sensor information, to locate single shallow underground sources.

Design/methodology/approach

First, in this paper, using the characteristics of low multipath interference and good P-wave polarization in the near field, the adaptive covariance matrix algorithm is used to extract the polarization angle information of the P-wave and the short term averaging/long term averaging algorithm is used to extract the first break travel time information. Second, a hybrid positioning model based on travel time and polarization angle is constructed. Third, the positioning model is taken as the particle update fitness function of quantum-behaved particle swarm optimization and calculation is performed in the hybrid positioning model. Finally, the experiment verification is carried out in the field.

Findings

The experimental results show that, with root mean square error, spherical error probable and fitness value as evaluation indicators, the positioning performance of this method is better than that without speed prediction. And the positioning accuracy of this method has been improved by nearly 30%, giving all of the three tests a positioning error within 0.5 m and a fitness less than 1.

Originality/value

This method provides a new idea for high-precision positioning of shallow underground single source. It has a certain engineering application value in the fields of directional demolition of engineering blasting, water inrush and burst mud prediction, fuze position measurement, underground initiation point positioning of ammunition, mine blasting monitoring and so on.

Citation

Li, J., Yan, X., Zhao, F. and Zhao, X. (2022), "The underground explosion point high precision measurement method based on multidimensional information fusion of vibration sensors", Sensor Review, Vol. 42 No. 3, pp. 281-293. https://doi.org/10.1108/SR-01-2020-0020

Introduction

Distributed-source positioning in shallow underground spaces is a popular subject in the field of underground space positioning (Li et al., 2019, 2015). With this method, many vibration sensors are laid in the near field of the underground source. By extracting the characteristic parameters of seismic phase obtained by each sensor, the source positioning is determined, as shown in Figure 1 (Kwon et al., 2005; Vogl et al., 2009). This is a close-range source measurement method, focusing on local near-field source positioning applications, such as geological monitoring, engineering blasting and antitheft surveillance of cultural relics (Zhang and Zhai, 2019).

With the development of underground location technology, a special type of location event, single-source location (Qin et al., 2019), has emerged. Examples include underground chamber blasting monitoring and ammunition fuze detonation location. For single-source location problems, however, no mature positioning algorithm exists. To address this problem, some researchers have turned to Geiger’s travel time method, which is used in natural seismic positioning (Geiger, 1912). Drawing on Geiger theory, Waldhauser and Ellsworth proposed a double-difference travel time location method to solve the problem of single-source location under the condition of uniform geology and known speed. However, in the case of complex geological conditions such as layered horizons or depressions, the average speed must be measured and the resulting location error becomes larger (Waldhauser and Ellsworth, 2000). Asgari et al. improved the travel time location method to propose a method using joint inversion of source location and velocity, which completes single-source location by updating speed and source position parameters (2015). This method accounts for source without requiring known speed. However, its location accuracy is affected by the initial modeling accuracy of the velocity field and the number of source excitations. Share et al. improved joint inversion of source location and velocity method to propose a single-source location method that does not rely on predicted speed. This method sets the velocity and source coordinates as unknown parameters to be found. A larger number of sensors are used to increase the redundancy of the travel time information, thereby reducing the sparsity of the location equations (Share et al., 2019). Although this method is capable of single-source location, it involves a greater number of deployed sensors, thus increasing the experiment cost and the difficulty of test deployment.

To address the above problems, this study determines how to improve the location accuracy of a single source using a small number of sensors. This paper makes full use of the 3D vector information obtained by 3D vibration sensors. As a shock wave propagation direction has good polarization, the P-wave polarization angle information is included advantageously into the positioning equations (Wang et al., 2010; Khajouei and Goudarzi, 2018). A multiparameter source location model containing polarization angle information and travel time information is established in this way. The quantum-behaved particle swarm optimization (QPSO) algorithm, a swarm intelligence optimization algorithm, is used to solve the source positioning model over a large scale and in a fast scanning mode and ultimately locate the source. This approach can be used for single-source location in applications with a limited number of sensors. The proposed method offers real-time monitoring of underground chamber blasting and rapid determination of the detonation position of ammunition fuzes.

Principle of source positioning based on travel time-polarization angle information

The polarization of a wave is the representation of the spatiotemporal characteristics of the wave field. The P-wave is linearly polarized in the near field of the explosion. The movement direction of the P-wave is consistent with the propagation direction of the wavefront (Li et al., 2017), as shown in Figure 2 (Si et al., 2018).

In a near field, the P-wave has good polarization and strong polarizability (Lu et al., 2010), so the polarization angle is introduced into the travel-time positioning model to build a hybrid positioning model based on travel time-polarization angle, as shown in Figure 3. First, the adaptive covariance matrix (ACM) algorithm is used to extract the P-wave polarization angle of each node and the STA/LTA algorithm is used to extract the arrival time of the first break. Second, with the above angle and time information, the travel time-polarization angle positioning equation is developed. Finally, the equation is evaluated for the source coordinates.

Extraction of positioning parameters

Arrival time of first break

The first break is a special type of wave. As a vibration wave transmitted to the receiving point directly from the source, it is characterized by an early take-off time and strong energy, thus providing an important basis for determining the arrival of the vibration wave (Tariq et al., 2019; Kumar et al., 2018). Based on the short term averaging/long term averaging (STA/LTA) method, this paper introduces higher-order statistics and a general S transform to construct an adaptive time window, removing the need to artificially determine the time window length. This method improves the pick-up accuracy of the first break while guaranteeing the time variance of the instantaneous energy factor.

The STA/LTA method is one of the most popular first break arrival time pick-up algorithms in engineering applications; its principle is shown in Figure 4 (Sabbione and Velis, 2010).

Let the long/short window recognition factor p(i) for the ith time point be:

(1) p(i)=∑j=1Nx(i+j)/N∑j=1Mx(i−j)/M

where x(j) (j = 1,2,…,N) is the vibration data and M and N are the number of samples in the long- and short-time windows, respectively.

To obtain the first break arrival time more accurately, higher-order statistical function is used for an enhanced mutation point and the S transform is performed to adaptively adjust the time window length to the minimum period of the signal in the time window, thereby amplifying the time variance of the first break detection factor. The improved first break detection factor R(i)′ consists of three parts:

(2) R(i)′=P(i)′.Q(i)′.K(i)

where P(i)′ is the instantaneous energy recognition factor under unequal window lengths, Q(i)′ is the instantaneous energy recognition factor under equal window length and K(i) is the kurtosis, a fourth-order statistical function, in the time window at the corresponding time. These terms are written as:

(3) P(i)′ =∑j=−N/2N/2X2(i+j)/N∑j=1MX2(i−j)/M

(4) Q(i)′=∑j=−N/2N/2X2(i+j)/N∑j=−N/2N/2X2(i−j)/N

(5) K(i)=E[(X(i+j)−E(X(i+j)))4]{E[X(i+j)−E(X(i+j))2]4}1/2   (j=−N/2,⋯0,⋯N/2) 

where the long window length M = 5N and the short window length N = k/Ω_i(t), k is the weight and Ω_i(t) is the instantaneous dominant frequency of the ith moment after the generalized S transform of the three-axis vibration signal, which refers to the frequency point with the highest energy among the wide spectrum corresponding to each moment after the generalized S transform. In the near field of the blasting point, when the first break arrives, the instantaneous dominant frequencies of the three axes are consistent (Mousa et al., 2011):

(6) Ωi(t)=Ωix(t)=Ωiy(t)=Ωiz(t)

where X(i) is the resultant of the three energy components:

(7) X(i)=x(i)×x(i)+y(i)×y(i)+z(i)×z(i)

To accent the waveform mutation of the first break, a recognition factor K(i) is constructed using the kurtosis.

P-wave polarization angle information

Here, the ACM algorithm is used to extract the angle information of the P-wave. This algorithm uses the multicomponent signal covariance matrix to determine the main movement direction of the P-wave. As the explosive vibration signal is a nonstationary time-varying signal, the length of the time window is adjusted in real time to suit the instantaneous frequency value of the vibration signal. This method does not require manual selection of the time window length. The specific extraction process using the ACM algorithm is shown in Figure 5.

The Hilbert transform spectrum uih(t) is used to find the analytical signal of the three-component signal u_i(t), (i = x,y,z) and the approximate expression of the multicomponent signal u_i(t) is used to construct the covariance matrix (Diallo et al., 2006):

(8) M¯(t)=|I¯xx(t)I¯xy(t)I¯xz(t)I¯xy(t)I¯yy(t)I¯yz(t)I¯xz(t)I¯yz(t)I¯zz(t)|

T_km(t) is the adaptive window length of M¯(t) at time t and is defined as follows:

(9) Tkm(t)=6πNΩx(t)+Ωy(t)+Ωz(t)=2πNΩavkm(t)

Let the eigenvalues of the covariance matrix M¯(t) be λ_i(i = 1,2,3), with λ₁ ≥ λ₂ ≥ λ₃. v_i is the corresponding characteristic vector of λ_i and the maximum characteristic value λ₁ corresponds to the principal characteristic vector v₁ (v_1x,v_1y,v_1z). Hence, the model for extracting the P-wave angle (Tian and Xu, 2015) is as follows:

(10) θ0(t)=arctan⁡[v1y(t)v1x(t)]

(11) δ0(t)=arctan⁡[v1z(t)v1x(t)2+v1y(t)2]

where θ₀(t) stands for the instantaneous azimuth and δ₀(t) is the instantaneous inclination. The polarization angle of the separated P-wave is extracted by the ACM algorithm, then the extraction is evaluated by the polarization angle.

Polarizability, T, is a metric to evaluate the extraction of the polarization angle of the P-wave. Generally, the polarization is in the range [0, 1] and the closer T is to 1, the better the polarization characteristic of the P-wave and hence the more accurate the extraction. When T is 0, the particle appears spatially as a sphere, that is, there is no polarization. The formula is as follows:

(12) T=(1−ρ2)2+(1−ρ12)2+(ρ2+ρ12)22(1+ρ2+ρ12)2

where ρ is the principal ellipticity polarization parameter and ρ₁ is the minor ellipticity polarization parameter.

Construction of positioning model based on travel time-polarization angle

The travel time-polarization angle positioning model is shown in Figure 6. The time difference between the P-wave reaching two sensor nodes (such as sensors i and j) describes a hyperboloid. The hyperboloid foci are on the two sensor nodes in question and the source is situated on this hyperboloid. Because the travel time error is large, when multiple hyperboloids meet, there must be a false appearance of location (Gambi et al., 2016; Bishop et al., 2008; Wang and Ho, 2017). Based on the principle of multibeam cross positioning, the P-wave polarization angle is introduced into the positioning model to correct the positioning error when the hyperboloids meet.

According to the travel-time positioning principle (Salari et al., 2018) and the geometric relationship between the sensors (Bishop et al., 2009), the following system of equations is established:

(13) {tan⁡γi=z−ziritan⁡βi=x−xiy−yir02=(x−x0)2+(y−y0)2+(z−z0)2ri2=(x−xi)2+(y−yi)2+(z−zi)2roi=r0−ri=vt0i

where the source coordinates are (x,y,z), the sensor coordinates are (x_i,y_i,z_i) (i = 1,2,3,n), r_i (i = 1,2,3,…,n) is the distance from the source to a sensor, (x₀,y₀,z₀) are the reference sensor coordinates, the P-wave speed is v and the difference in times between the reference sensor and each sensor is t_0i. The pitch angle and the azimuth between sensor i and the source are γ_i and β_i, respectively, with i being the serial number of the sensor node.

Solving QPSO-based source positioning model

Principle of QPSO algorithm

QPSO, an improvement on the PSO algorithm, is a random optimization technique based on swarm intelligence (Parvin and Vasanthanayaki, 2019). The most prominent features of the QPSO algorithm are its fast convergence and high global optimization capability. It regards the PSO system as a quantum system and builds an attractive potential field with local attractor p_i,j as the center, so that particles in the bound state of the potential field are able to search, with a probability of 1, for any position within the feasible region. The state model of the QPSO algorithm is shown in Figure 7 (Sun et al., 2012).

In the QPSO algorithm, the speed and position of a particle are all dependent on one parameter. To ensure the convergence of the algorithm, each particle must converge on its own point P, with P = (P_1,P₂..P_d) and P_d being the value of the particle in the dth dimension (Wu et al., 2018):

(14) Pd=(ϕ1∗Pid+ϕ2∗Pgd)(ϕ1+ϕ2)

where φ₁ and φ₂ are random functions between 0 and 1. A median optimal position mbest, defined as the average of the global extreme values of all particles, is introduced to evaluate the variable L in the next iteration of the particle. The formula is as follows:

(15) mbest=1M∑i=1MPi=(1M∑i=1MPi1,…1M∑i=1MPid)

where M is the number of particle swarms and P_i is the global extremum of particle i. Therefore, parameter L is determined by the following equation:

(16) L(t+1)=2*β*|mbest−x(t)|

Hence, the evolution process of the particles can be obtained as follows:

(17) xi,j(t+1)=pi,j±β*|mbest−x(t)|*ln⁡(1/u),u=rand(0,1)

where β is the coefficient creativity. By adjusting this value, it is possible to regulate the convergence speed of the algorithm. Generally, the algorithm can give better results when β decreases linearly from 1.0 to 0.5 (Gan et al., 2018).

Establishing fitness function of the particle swarm using travel time-polarization angle information

The fitness function is used to select and update particles. Here, the source positioning model based on travel time-polarization angle is transformed into a fitness function for the particle swarm, as shown in equation (18). The fitness function is constructed using the polarization angle and travel time information of the four sensors giving the highest polarization, as shown in equation (18):

(18) f(x,y,z,v)=∑i=14((ri−r0−vti0)2+(xi−x+(y−yi)tan⁡βi)2+(z−zi−ritan⁡γi)2)

where (x, y, z) is the target position and v is the propagation speed of the P-wave; the expressions of r_i and r₀ are as shown in equations (19) and (20), respectively:

(19) ri=(x−xi)2+(y−yi)2+(ritan⁡γi)2

(20) r0=(x−x0)2+(y−y0)2+(ritan⁡γi+zi)−z0)2

where λ_i is the pitch of sensor i and the source, β_i is the azimuth of sensor i and the source, (x_i, y_i, z_i) are the position coordinates of the sensor node i = 1,2,3,4, (x₀, y₀, z₀) are the position coordinates of the reference sensor node and t_i0 is the measured time difference between the ith base station and the master station.

Specific implementation process

The QPSO algorithm flow is shown in Figure 8:

Step 1: Set the search range for the source and the velocity, the population size to 40, the number of spatial dimensions to 4 and the number of iterations to 1000 and randomly generate the initial particle swarm;

Step 2: From the objective function (18), find the position with the smallest fitness value in the source group and take it as the optimal source Gbest of the current generation of the source group;

Step 3: From Equation (15), find the average optimal position mbest of 40 sources;

Step 4: Update the location of each source using Equation (17);

Step 5: When the number of iterations is reached, Gbest is output and taken as the optimal source location.

Experimental validation

To verify the advantages of the proposed method, a shallow underground explosion positioning experiment was performed at the test and measuring academy. The proposed method was then compared with source location methods with no predicted speed. Root mean square error (RMSE), spherical error probable (SEP) and fitness were taken as the evaluation indexes of the positioning accuracy. The two methods were evaluated with three location tests. Three rounds of TNT explosives were buried in a shallow space of 20 m × 20 m × 5 m and 10 omnidirectional vibration sensors, made by North University of China, were used to form a vibration sensor array to obtain the 3D vibration signal of the underground explosions. The indices of the sensors are shown in Table 1.

The sensors were laid out as described in Tables 2‐3. The layout, including burst point positions, is shown in Figure 9. Three rounds of TNT explosives were buried in different areas to simulate independent underground initiation point positioning tests. Burst point 1 was in the third quadrant, burst point 2 was in the first quadrant and burst point 3 was on the origin of the coordinate.

A general horizontal datum was set up on the ground according to preset coordinates. The sensor nodes were sequentially installed according to the sensor layout diagram. The sampling rate of the underground burst point positioning system was set to 20 kHz and its sampling time was set to 10 s. The site was allowed to stand for one day and three blasting experiments were carried out successively. The extraction of the positioning parameters is described below, using as an example the vibration signals obtained after the first explosion by sensor nodes 5 and 9.

Figure 10 shows that multiwaveform aliasing exists in the 3D vibration signals acquired by nodes 5 and 9. For the two nodes, recognition factor graphs of the first break arrival time were extracted using the STA/LTA algorithm, as shown in Figure 11.

In extracting the near-field first break arrival time, as shown in Figure 11, there are multiple peak points, so the issue of potential misjudgment arises. We took the time corresponding to the peak as the first break arrival time to determine the first arrival times of the node group, as shown in Table 4.

To simulate three independent positioning tests of the underground burst point, the time interval was set to 30 min. Each test started from time 0. The external trigger mode was used to control the data acquisition system and the trigger signal line was connected to the explosive. When the explosive exploded, the trigger signal was generated, enabling the data acquisition system to start collecting signals. Therefore, the first break arrival time in Table 4 was the relative time that the shock wave propagated from the burst point to the sensor.

After the travel times were obtained, the polarizability of the two nodes was evaluated. The first break arrival times of the node group were all around 3.04 s, so the analysis in the instantaneous polarizability diagram focuses on the polarizability transformation after that moment, as shown in Figure 12.

As shown in Figure 12, in the effective P-wave region, the polarizability of nodes 5 and 9 both exceed 0.97, suggesting good polarization and high extraction accuracy of the polarization angle. The curves of the instantaneous polarization angles of the P-wave beam between the coordinate axes, as obtained with the ACM algorithm in the Cartesian coordinates, are shown in Figure 13.

As shown by Figure 13, the triaxial polarization curve is stable within the effective direct P-wave time period. The abrupt changes of the angle at the aliasing moment of the P-wave final phase and the S-wave first phase are obvious. The polarization angles for all nodes, found with the above method, are shown in Table 5.

As shown in Table 5, the polarizability varies greatly among the sensor node group. The figures in blue indicate sensors with poor polarization, which suggests poor coupling between the sensor and the soil and the resulting amplitude characteristics are therefore unable to characterize the movement direction of the sensor. The figures in red represent groups of sensors with the best polarizability.

The travel time information corresponding to ten sensors was used to construct a travel-time positioning equation with no predicted speed, as shown in equation (21) (Share et al., 2019):

(21) (x−x1)2+(y−y1)2+(z−z1)2−(x−x0)2+(y−y0)2+(z−z0)2=vt10;(x−x2)2+(y−y2)2+(z−z2)2−(x−x0)2+(y−y0)2+(z−z0)2=vt20;(x−x3)2+(y−y3)2+(z−z3)2−(x−x0)2+(y−y0)2+(z−z0)2=vt30;⋯⋯⋯⋯⋯(x−xi)2+(y−yi)2+(z−zi)2−(x−x0)2+(y−y0)2+(z−z0)2=vti0;

where (x, y, z) are the coordinates of the underground source, v is the propagation speed of the direct P-wave throughout the underground local space composed of the five sensor nodes and t_i0 is the time difference between the P-wave reaching the ith sensor node and the reference node, with i = 9.

The two positioning equations were evaluated using QPSO. The specific process was as follows: First, set the search range of the source to x (−50 m, 50 m), y (−50 m, 50 m), z (−50 m, 0 m) and v (100 m/s, 1,000 m/s). The population size was set to 80, the number of spatial dimensions to 4 and the number of iterations to 1,000, to randomly generate the initial particle swarm. Second, with the positioning equation converted into the objective function, the particle with the minimum fitness value in the source swarm was found and taken as the optimal particle Gbest of the current generation of particle swarm. Finally, Gbest was continuously updated to produce the source coordinates corresponding, respectively, to the proposed polarization angle-travel time location method and the positioning method without predicted speed.

The real-time positioning accuracy of the two methods over 1,000 iterations was measured in fitness and SEP. Fitness is a key indicator for evaluating the quality of particles. The smaller the value, the better the particles and the higher the positioning accuracy. SEP is a key indicator to measure positioning stability; it is expressed as follows:

(22) SEP=0.75σx2+σy2+σz2

where σ_x, σ_y and σ_z are the positioning standard deviations along the three axes. The smaller their values, the more stable the positioning model. The search curves corresponding to the two indicators during 1,000 iterations are shown in Figure 14.

Figure 14 shows that the fitness value of the proposed method is always lower than that of the source location method with no predicted speed. With the proposed method, after the first 100 iterations, fitness begins to stabilize below 1.0; the source location method with no predicted speed develops prematurely and its fitness eventually stabilizes at 40. With the SEP as the evaluation index, the results of the proposed method are closer to the true source coordinates and the radius of the SEP is significantly smaller than that in the method with no predicted velocity. The positioning accuracy evaluation indices – the final positioning result, RMSE, maximum SEP and minimum fitness – are given in Table 6.

Table 6 shows that the proposed travel time-polarization angle source positioning method improved the positioning accuracy by almost 30%, giving in all three tests a positioning accuracy within 0.5 m, an SEP within 0.25 m and a fitness less than 1. These indicators are superior to those of the source location method without speed prediction.

In Table 6, the relative error along the x-axis is small, whereas the relative error along the y-axis is large. According to the geometric dilution of precision correlation principle (Zhang and Lu, 2020), this is related to the geometric layout of the sensor array. Combined with the data shown in Figure 9(a), in the x-axis, the sensor array can obtain complete positioning information on both sides of the burst point. However, in the y-axis, the sensor array can only obtain positioning information on one side of the burst point, resulting in a lack of positioning information in the y-axis. Therefore, there exists positioning deviation in the y-axis. To further improve positioning accuracy, future research will determine how to optimize the geometric layout of the sensors.

Conclusions

This paper proposes a method to locate single shallow underground sources based on fusion of multidimensional vibration sensor information. Considering the fact that in the near field of the explosion a shock wave propagation direction has good polarization, the P-wave polarization angle information is included advantageously into the positioning equations. A source location model with a greater degree of accuracy is constructed using both polarization angle information and travel time information. The source positioning model is solved using the QPSO algorithm. The test results suggest that, compared with the source positioning method without predicted speed, the proposed method improved location accuracy by 30%, with all positioning errors less than 50 cm. The proposed method can be used to monitor underground chamber blasting and rapidly determine the detonation position of ammunition fuzes.

This method still has certain shortcomings. As explosion energy increases, the signal obtained by the sensor is saturated, which will cause a significant decrease in the extraction accuracy of the P-wave polarization angle. Further, Table 6 shows that the geometric layout of the sensor array also affects the positioning accuracy of the burst point. Therefore, assuming that the explosive equivalent is known, future research will focus on the optimization of the sensor array layout to eliminate the positioning blind spot and positioning ambiguity.

Figures

Figure 1

Schematic diagram of shallow ground distributed source localization

Figure 2

Polarization characteristics of a shock wave

Figure 3

Overall flowchart for source positioning

Figure 4

Schematic of the sliding window used in the energy factor method

Figure 5

Flowchart of the extraction of the direct P-wave angle information

Figure 6

3D positioning diagram

Figure 7

State model of the QPSO algorithm

Figure 8

QPSO algorithm execution flowchart

Figure 9

Field layout

Figure 10

Time domain waveforms of two nodes

Figure 11

First break recognition factor graph

Figure 12

Instantaneous polarizability diagrams

Figure 13

Triaxial polarization angles of node 6

Figure 14

Evaluation curve of positioning accuracy

Table 1

Indices of the omnidirectional vibration sensors

Indices	Details
Sensor bandwidth (−3 dB)	0–1,800 Hz
Sensor range (g)	±6 g
Sensor accuracy (% FS)	±1
Inter-axis crosstalk (% FS)	1
Physical indices	Volume: ϕ3 cm; density: 1.40 g/cm3; power supply: 5 V ± 5%; total power consumption: <10 mW; output impedance: <100 ohm

Table 2

Coordinates of the sensors

Serial number of sensors	Actual coordinates (m)
	X	Y	Z
0	−4.974	4.931	−1.121
1	−5.183	2.896	−1.687
2	−6.108	3.461	−1.304
3	−6.946	3.964	−1.292
4	−6.767	1.764	−1.096
5	6.048	0.000	−1.700
6	7.055	0.018	−1.403
7	7.753	2.151	−1.110
8	6.931	4.072	−1.090
9	7.751	4.548	−0.802

Table 3

Positions of burst points

Serial no. of	Explosive	Actual coordinates (m)
bust points	quantity (kg)	X	Y	Z
1	4	−8.636	−5.100	−15.000
2	4	−0.567	7.995	−20.000
3	4	0.000	0.000	−10.000

Table 4

First break arrival times of the node group (in s)

Serial no.	First explosion	Second explosion	Third explosion
0	3.03995	2.29690	3.67065
1	3.03885	2.29495	3.66780
2	3.03860	2.29620	3.66950
3	3.03890	2.29700	3.67015
4	3.03995	2.30020	3.67185
5	3.04185	2.29810	3.67315
6	3.04385	2.29895	3.67460
7	3.04530	2.29935	3.67610
8	3.04500	2.29815	3.67080
9	3.04720	2.30095	3.67320

Table 5

P-wave angle information of the node group (in radian)

		First explosion		Second explosion		Third explosion
No.	Polarizability	Pitch angle	Azimuth	Pitch angle	Azimuth	Pitch angle	Azimuth
0	0.7920	0.0848	0.3539	0.0903	1.1107	0.0734	2.3744
1	0.1572	0.0774	0.4620	0.0362	0.7730	0.0530	2.0774
2	0.9539	0.0622	0.3280	0.0753	0.9680	0.1019	2.0857
3	0.9123	0.0580	1.2120	0.0300	2.4860	−0.0388	1.6028
4	0.2553	0.1225	0.2693	0.0990	0.8115	0.1316	1.8607
5	0.9856	0.0767	0.1782	0.0941	1.0723	0.0537	2.0829
6	0.9618	0.0562	1.2341	0.0540	2.4153	0.0985	1.5533
7	0.9913	0.0642	1.1457	0.0923	2.2115	0.1242	1.3014
8	0.9903	0.0619	1.0204	0.1151	2.0834	0.1558	1.0639
9	0.9800	0.0746	1.0062	0.1294	1.9141	0.1736	1.0304

Table 6

Results of underground burst point positioning (m)

Sequence of blasting	Algorithm	X (m)	Y (m)	Z (m)	RMSEa (m)	SEP (m)	Fitness value
1	Actual source coordinates	−8.6360	−5.1000	−15.0000	–	–	–
	Travel time-polarization angle positioning method	−8.9912	−5.2094	−15.1851	0.4152	0.2481	0.4820
	Positioning method without speed prediction	−9.2273	−5.1298	−15.2170	0.6306	0.5942	32.2007
2	Actual source coordinates	−0.5670	7.9950	−20.0000	–	–	–
	Travel time-polarization angle positioning method	−0.7701	8.1813	−20.3354	0.4341	0.2241	0.2863
	Positioning method without speed prediction	−0.8305	7.7905	−19.4042	0.6828	0.5946	34.4301
3	Actual source coordinates	0.0000	0.0000	−10.0000	–	–	–
	Travel time-polarization angle positioning method	−0.0526	−0.1158	−10.3986	0.4184	0.2461	0.0091
	Positioning method without speed prediction	−0.1168	−0.0037	−9.3631	0.6475	0.5506	15.0606

Note:

^aRMSE: root-mean-square error