Detection of Rock Mass Defects Using Seismic Tomography

Abstract

Elastic wave tomography is a non-destructive technique used to obtain velocity images of internal structures. Wave velocity tomography has been used to quantitatively describe abnormal structures and identify boundaries in different media and shapes. This research aims to quantitatively describe anomalous structures on a laboratory scale and identify boundaries in different media and shapes. Wave velocity tomography technology was applied to non-empty and empty anomalous areas, and the accuracy of the resulting velocity images was quantified. The reliability of boundary recognition was studied and analyzed. The accuracy of the results was found to relate to the size and shape of the anomalous structures estimated in the prior model, as well as the wave velocity difference between the anomalous and normal regions. In practical engineering, when the presence of abnormal areas cannot be determined, a completely homogeneous prior model is used. The results indicate that, although the inversion velocity may not be accurate, it is sufficient to identify boundaries. The larger the velocity difference, the clearer the boundary to be recognized. Furthermore, to assess the practical efficacy of this physical detection method, the wave velocity tomography technology was used to detect hidden goafs within a tantalum–niobium mine belonging to the Hunan Fuguihengtong Mining Company. Utilizing physical detection, 12 potential goafs and 100 confirmed goafs were identified within the mine, demonstrating the efficacy of the proposed empty anomalous area detection methodology in accurately identifying and delineating the boundaries of hidden goafs.

1. Introduction

Elastic wave tomography is a non-destructive technique used to obtain velocity images of interior structures. This technique is widely used in engineering projects for geotechnical characterization. For example, Leiceaga et al. (2015) [1] used cross-well tomography to verify the integrity of reservoir sites, Olivier et al. (2015) [2] investigated seismic processes using in situ measurements of seismic velocity variations in an underground mine, and James et al. (2017) [3] used distributed acoustic sensing to detect fractures through relative seismic velocity changes.

With the continuous advancement and maturation of technology, the tomography of seismic wave velocity, amplitude, and time characteristics has gradually become an important technique for understanding the internal structure of rock masses or rocks [4], and this technique has also been effectively applied to rocklike materials in the laboratory to study micro-crack development (e.g., [5,6,7,8,9]), identify the location of defects (e.g., [10,11,12,13]), and validate numerical simulations (e.g., [14,15,16]).

Compared to acoustic emission (AE) technology, which is also widely used for the non-destructive testing of internal defects in rock masses, seismic tomography primarily relies on recording the propagation velocity and attenuation characteristics of seismic waves (typically P-waves and S-waves) in geotechnical media. By recording reflected or refracted waves from multiple seismic sources and receivers, it allows for the inversion of subsurface structural features. In contrast, AE is an active monitoring technique, typically used to detect small-scale damage or stress concentration areas within rock masses or structures, focusing on localized dynamic events within materials rather than constructing an image of the entire subsurface structure [17]. Although there are differences in principles and applications, the combination of both methods enables more comprehensive non-destructive testing of internal defects in rock masses in geotechnical engineering. In 1993, Jansen et al. [5] combined tomographic difference imaging and acoustic emission techniques. They proved that the combination of these methods is useful for evaluating induced damage in rock. Stanchits et al. (2011) [8] used this combination to monitor brittle failure initiated by water injection. Velocity imaging and injected volume measurements confirmed the close relationship between the water front and the induced AE cloud. To understand the physical mechanisms governing CO₂–water replacement, Lei and Xue (2009) [10] carried out seismic measurements in porous sandstone under laboratory conditions. They confirmed that, when employed with differential arrival times and relative amplitudes determined from the cross correlation of waveforms, the difference tomography method can efficiently image changes in velocity and attenuation associated with the partial replacement of water. To better understand rock fracture and fault instability driven by high-pressure fluid sources, Lei et al. (2016) [15] conducted rock fracture tests, using granite under a confining pressure with fluid injection in a laboratory setting. Their study, which incorporated laboratory experiments and numerical simulations, indicates that such an approach is helpful for finding a better model in terms of both simulating experimental results and for upscaling purposes.

The studies mentioned above qualitatively show that most rock physical properties can be tomographically imaged, including pore fluid content, pore pressure, rock composition, temperature, elastic moduli, and water saturation, which also demonstrates the broader applications of seismic tomography beyond mining engineering. For example, by analyzing the velocity and attenuation characteristics of rock mass, it is possible to identify groundwater flow in different strata, as well as changes in porosity and permeability within the rock layers, thereby enabling the detection and monitoring of potential water resources (e.g., [18,19]). Additionally, during tunnel excavation and operation, imaging the velocity distribution around the tunnel can reveal potential geological weaknesses and fractures, aiding in the assessment of tunnel construction safety (e.g., [20]).

The purpose of our study is to quantitatively delineate anomalous structures whose physical properties are different from the surrounding rock mass by using seismic tomography on a laboratory scale, and identify the boundaries of anomalous areas, finally achieving non-destructive detection of rock mass defects. The inversion program in our study involves four prerequisites:

• All the receivers are used as active transducers, providing rays with known source positions;
• The receivers can be located in three dimensions;
• The real velocity structure demonstrates strong transverse or orthotropic anisotropy;
• The velocity contrasts are strong and occur over very short distances (millimeters).

2. Methodology

2.1. Theoretical Arrival Time Calculation

The first step was to accurately compute the arrival times in a complex, anisotropic, and strongly heterogeneous medium. The velocity and anisotropic structures were defined along a regular orthonormal grid with relatively coarse spacing, typically larger than the wavelength of the elastic waves used to probe the system. We chose the fast-marching algorithm, which allowed us to compute arrival times at all nodes in a given grid from a source. The implementation used second-order finite differences to solve the anisotropic Eikonal equation. To mitigate inaccuracies regarding the source point, the trial arrival times at all points in a cubic box surrounding the source point were computed using analytical solutions, assuming a constant velocity and an anisotropy parameter equal to the source [21]. The wavefront for an active source (Figure 1, left) and a passive acoustic emission (Figure 1, right) on a 100 mm × 100 mm anisotropic 2D plane were calculated using the above method.

2.2. Ray Tracing

Once the arrival times were computed at all points, it was relatively easy to compute the ray paths a posteriori. In the isotropic case, ray tracing is straightforward because the ray angle at every point is perpendicular to the wavefront. In the anisotropic case, the rays are no longer perpendicular to the wavefront since the phase angle is not the same as that of the group angle. Therefore, at every point, it was necessary to compute the group angle (from the phase angle, which was given by the gradient of the arrival time field) to obtain the ray orientation [21].

2.3. Tomography Procedure

The inverse problem can be expressed as

$$d=g\left(m\right),$$

(1)

where $d$ is a vector of observable quantities, $\mathrm{m}$ is a vector of model parameters, and g is the physical model relating to observable variables (the measured arrival times for all acoustic emissions, as well as active surveys, at a fixed array of receivers).

The model parameters are the velocity and anisotropy parameter at all nodes in a 3D grid and the position of the acoustic emissions. Anisotropy is described using a very simple elliptical vertical transverse isotropic model. The vector of model parameters is formed by the sequence of $V_h$ values of each node; the E values at each node; and quadruplets $x, y, z, t_0$ for each acoustic emission. The prior model parameter vector is denoted mprior. As the wave velocity is Jeffrey’s parameter, the logarithm of $V_h$ is used as the corresponding cartesian parameter [21]. The covariance matrix in the model space, CM, is determined by the variance ${\sigma }_v^2$ and covariances on $\mathrm{l}\mathrm{n}\left(V_h\right)$; the variance ${\sigma }_E^2$ and covariances on $E$; and the variances ${\sigma }_x^2$, ${\sigma }_y^2$, ${\sigma }_z^2$, and ${\sigma }_{t_0}^2$ on source locations. The inverse problem is solved iteratively using the quasi-Newton method, assuming Gaussian errors in all data and model parameters.

3. Synthetic Experiment

Following the workflow illustrated in Figure 2 for identifying rock mass anomalies and their boundaries using seismic tomography at the experimental scale, we carried out wave velocity inversions on a small scale with the purpose of imaging the wave structure and demarcating the boundaries according to the wave velocity difference. Each synthetic test described below was completed in the following six stages:

1.

Prepare the original model, including setting the velocity corresponding to each node in the normal and abnormal zones, arranging the sensor array, and generating random acoustic emission events.

2.

Calculate the arrival times for each AE event source and active ultrasonic source to each sensor using the forward method described in Section 2.1.

3.

Prepare the prior model for the inversion.

4.

Configure the constraints for inversion, including the grid refinement factor; the standard deviation of survey picks and event picks; the maximum number of iterations; residual tolerance to ceasing iteration; the quasi-Newton step size; the standard deviation ${\sigma }_V$ of $\mathrm{l}\mathrm{n}\left(V_h\right)$; the standard deviation ${\sigma }_E$ of $E$; the correlation length $\lambda$; the standard deviations ${\sigma }_x$, ${\sigma }_y$, and ${\sigma }_z$ of the source positions; and the standard deviation ${\sigma }_{t_0}$ of the original time $t_0$ of the sources.

5.

Apply the method described in Section 2.3 for the tomography inversion.

6.

Evaluate the accuracy of the velocity images and identify the boundaries of the abnormal area according to the tomography results. The area defects (e.g., sectional defects in different shapes) and structure defects (e.g., linear defects with different widths) will be tested separately.

Figure 2. The schematic diagram of the workflow for identifying rock mass anomalies and their boundaries using seismic tomography at the experimental scale.

Figure 2. The schematic diagram of the workflow for identifying rock mass anomalies and their boundaries using seismic tomography at the experimental scale.

3.1. Delineation of Non-Empty Anomalous Areas

A total of 12 synthetic tests were carried out for four geometric sizes, multiplying three velocity structures. The geometry size of the original model is provided in Figure 3a, which shows a circular area with a radius of 30 mm in a 100 mm × 100 mm square plane. The abnormal region mentioned in this section refers to the circular region. The area outside the circle is defined as the normal region. The wave velocity of the normal zone is 4.5 km/s. We tested three scenarios with abnormal region wave speeds of 1.5 km/s, 2.5, km/s, and 3.5 km/s. The actual radius of the anomalous zone is considered unknown information in the inversion. When preparing the prior models, we successively set the initial iteration radius to 10 mm, 20 mm, 30 mm, and 40 mm, as shown in Figure 3b. In each test, the wave velocity in the abnormal area in the prior model is the same as in the original model; the only difference is the range of the abnormal area.

Twelve sensors are arranged around the plane at positions of 10 mm, 50 mm, and 90 mm on each side (shown as triangles in Figure 3). In total, in this 100 mm × 100 mm square plane, 100 AE event locations are randomly generated using a uniform distribution (indicated by green crosses in Figure 3), ensuring each position within the given area has an equal probability of event occurrence. After event generation, the distances between adjacent events are calculated to verify whether these distances follow a uniform distribution or meet the expected distribution, thereby confirming the spatial randomness of the events. The plane is divided into 40 × 40 coarse grids, and each coarse cell is divided into a 10 × 10 refined grid. The standard deviation is 0 for both survey picks and event picks. The maximum number of iterations, the residual tolerance to ceasing iteration, and the quasi-Newton step size are set to 10, 1 × 10⁶, and 1.0, respectively. The standard deviations on $\ln \left(V_h\right)$ and $E$ are 0.1. The standard deviations on ${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z,$ and ${\sigma }_{t_0}$ are 0.2. The correlation length $\lambda$ is 5. According to the method described in Section 2 and the parameter configurations introduced above, the wave velocity field in each scenario is inversed.

Figure 3 shows the inversion results for the scenarios investigated using the original model (an abnormal zone radius of 30 mm with wave speeds of 1.5 km/s, 2.5 km/s, and 3.5 km/s) and four prior models (an abnormal zone radius of 10 mm, 20 mm, 30 mm, and 40 mm). The inversion results were examined in two stages: first, the accuracy of the velocity images was quantified; second, the reliability of the boundary recognition was analyze.

To analyze the accuracy of the inversions, we first determined the velocity value of each node in the inversion results and separately calculated the velocity distribution within the actual range of the normal area (inner radius: 30 mm) and the abnormal area (outer radius: 30 mm), as shown in Figure 4a. The inversion accuracy of the normal area is significantly higher than that of the abnormal area. We noted that, except for Tests 4 and 8, the inversion results in the normal area are all distributed around 4.5 km/s. The inversion accuracy of the anomalous zone is affected by the initial iteration size (prior model) and the true normal–abnormal velocity difference. For the abnormal area, only the inversion results of Tests 3, 7, and 11 are close to the actual value, showing that the estimated size in the prior model has a strong influence on the accuracy of the inversion results. In other words, it is necessary to use other methods to estimate the spatial extent of anomalous structures before performing wave velocity inversion. The errors in Tests 9 to 12 are the smallest, followed by Tests 5 to 8, and the errors in Tests 1 to 4 are the largest. The results indicate that the greater the difference in wave velocity between the abnormal and normal areas, the lower the inversion accuracy.

In addition to the statistical comparisons, we compared the calculated P-wave velocity value ${\mathrm{V}}_{\mathrm{t}}$ of each grid node with the actual value $V_o$ of the original model. If the difference between them is less than ±20%, the tomography result is considered valid; see Figure 4b. It can be seen that the inversion of the nodes in the normal area is accurate. However, under the influence of the prior geometry and the actual velocity difference, the error of each node in the abnormal area of Test 1 is more than 20%; that is, the inversion result is unreliable. Tests 2, 4, 5, 6, and 8 have an error of less than 20% of the nodes in the abnormal area and less than 600% of the nodes. Less than 60% of the nodes in the abnormal areas of Test 2, 4, 5, 6, and 8 have an error of less than 20%. In order to ensure the reliability of the inversion results, the velocity difference in the real model cannot be too large.

The main purpose of wave velocity inversion is to identify abnormal structures, especially their boundaries. It can be seen from Figure 3 that the smaller the velocity difference in the real model, the more accurate the velocity in the anomalous area obtained from the inversion, and the greater the velocity difference in the real model, the more obvious the boundary of the abnormal area in the inversion result. Clear boundaries can be seen at a radius of 30 mm in the velocity images for all cases, except for Test 4. However, there are many cases (e.g., Tests 1, 2, 5, and 8) for which there is more than only one single clear boundary in the inversion result. In fact, Figure 3 shows that the largest velocity difference nodes (potential boundaries) in the inversion results do not necessarily appear at the real junction between the normal and abnormal areas, but at the initial boundary established by the prior model. The first two columns and the fourth column in the first row of Figure 3c demonstrate this particularly clearly.

Considering that the abnormal area in our model is a circle, when the source is not inside the abnormal area, the ray traces will pass through the boundary in a tangential manner. When the source is located inside the abnormal area, the ray trace will have a larger corner angle at the boundary nodes. No matter which of the above situations applies, the rays will be gathered in the nodes at the boundary (this is also the principle of visual boundary identification described previously). Therefore, we calculate the number of intersections of the ray paths in the nodes and identify the boundary from their density.

The number of intersections in each node is first counted. The centers of the nodes with the most intersections in each row are selected as the identified boundary points, and these identified points are connected into an enclosed polyline (EP). The coverage inside the EP is identified as the abnormal area, and each node inside is labeled ‘abnormal’. The coverage outside the EP is the identified normal area, and each node is labeled ‘normal.’ As each node in the original model has a real label, the correct rate of the identified boundary can be calculated by comparing the recognized label with the real label:

$$\left\{\begin{array}{c}{cr}_{nor}=N_1/N_2\\ {cr}_{abn}=N_3/N_4,\\ {cr}_{tot}=N_5/N_6\end{array}\right.$$

(2)

where cr_nor, cr_abn, and cr_tot are the correct rates for the normal area, abnormal area, and total area, respectively; N₁ is the number of nodes that were recognized as normal and fall in a truly normal area; N₂ is the actual number of nodes in the original normal area; N₃ is the number of nodes that were recognized as abnormal and fall in the truly abnormal area; N₄ is the actual number of nodes in the original abnormal area; N₅ is the number of nodes whose identified labels are consistent with the original labels; and N₆ is the total number of nodes in the model. The accuracy of the identified boundary for each inversion case is shown in Figure 5. It can be seen that, except for Tests 4, 8, and 9, the accuracy rates of the boundary identification results are all above 90%, regardless of whether they were determined for the normal area, abnormal area, or overall performance. Test 9 has the lowest accuracy of the identified abnormal boundary. By comparing the parameters of Test 9 with the other cases, it can be determined that when the real velocity difference is small and the shape in the prior model is inaccurate, the boundary error obtained via identification is obviously large.

As we do not know the real boundary in practice, the abnormal structure and its border cannot be determined using a single inversion. It is possible to reveal the potential boundaries through multiple inversions, comparing the results each time. Even though we do not know the true shape (a circle with radius of 30 mm) of the abnormal area, we can describe the geometry of this abnormal area by comparing the 12 inversion results.

The above analysis was based on the following assumptions: (1) We know there are abnormal areas. (2) We know the wave velocity of the abnormal area (the wave velocities of the abnormal area in the original and prior models are the same). (3) We do not know the exact range of the abnormal area. Moreover, in practical engineering, we do not know whether or not there are abnormal areas or what the wave velocity of rocks in abnormal areas is, let alone the rough geometry of abnormal areas. Therefore, we included the following tests in our method.

The shape of the original model is the same as in Figure 3a, and the event and sensor locations are consistent with previous tests. The wave velocities in the origin anomalous area are settled as 1.5 km/s, 2.5 km/s, and 3.5 km/s, respectively. However, the prior model is completely homogeneous as we do not know if an abnormal area is present. The results are shown in Figure 6.

The homogeneous prior models used the average velocity of the normal and abnormal areas. From Figure 6, we can infer that although the inversion velocity is not necessarily accurate when using a homogeneous model, it is enough to identify the boundary. A greater velocity difference in the original model results in a clearer boundary identification. For each of the three cases, the accuracy of the intersection density method for identifying boundaries is over 90%.

3.2. Delineation of Empty Anomalous Areas

We tested more extreme cases in which the abnormal area is empty (velocity = 0 km/s). Our approach included tests of an empty circle, square, and triangle. In the original model, the normal area wave velocity is 4.5 km/s, the abnormal area wave speed is 0 km/s, and the prior model is isotropic and homogeneous, with a wave speed of 4.5 km/s.

As shown in Figure 7, the inversion accuracy of the velocity field gradually decreases from the circle to the square and triangle. From the ray distribution and intersection points, the identified boundary of the circular void is the clearest. Only an approximate boundary can be recognized when a rectangular anomalous area forms, for which the identified boundary accuracy rate is only 80%. In the case of the triangular anomalous area, there is a large error in the recognized boundary, and the overall accuracy rate is only 70%. That is, the rounder the abnormal area, the easier it is to identify the boundary, and the more angular the shape, the harder it is to describe clearly.

4. Application

4.1. Study Area

This study was conducted in a tantalum–niobium mine belonging to the Hunan Fuguihengtong Mining Company. The mine is located in Gaolong Town, 50 km northeast of Chaling County, adjacent to the borders of Hunan and Jiangxi. The geographical coordinates are as follows: east longitude, 113°45′29″~113°47′05″; north latitude, 27°02′30″~27°03′28″. Three groups dominate the north of the mining area’s tectonic line (NNE, NNE, and NNE). In this area, folds and faults are developed, the latter of which predominate, and the metallogenic conditions are superior.

According to historical data on the mine and the results of the field investigation, a large number of underground mining tunnels and gob areas were generated in the early underground mining period in what is currently the tantalum and niobium mine (as shown in Figure 8). However, due to the limitations of the previous investigation technology and management level, some hidden gob areas remain undetected and must be identified to ensure the safety of future tantalum and niobium open-pit mining operations. It was necessary to define and identify the underground goafs formed due to historical mining in the mine area. In this study, wave velocity tomography technology was used to identify the boundaries of empty and non-empty anomalous areas to delineate hidden goafs in the mining area.

4.2. Data

In the physical detection of goafs, we employed three-component seismic geophones to capture seismic wave arrival times and waveforms. A total of 27 sensors were deployed in areas with higher microseismic activity and key tunnels within the tantalum–niobium mine, connected via an underground wired network to form a seismic monitoring network. This setup enabled the real-time acquisition and transmission of seismic wave signals, facilitating data analysis and early warning for monitoring personnel. The sensor deployment layout is as illustrated in Figure 9, and their precise deployment coordinates can be found in

Table 1. Coordinates of sensor placement.

Sensor Serial Number	China Geodetic Coordinate System 2000
	X	Y	Z
1	2,993,105.819	38,477,928.200	656.803
2	2,993,170.253	2,993,170.253	648.763
3	2,992,894.087	2,992,894.087	655.071
4	2,992,595.972	2,992,595.972	653.044
5	2,993,121.047	2,993,121.047	716.031
6	2,993,079.673	2,993,079.673	755.657
…	…
23	2,992,987.938	2,992,987.938	894.922
24	2,993,043.968	2,993,043.968	881.048
25	2,993,097.857	2,993,097.857	894.841
26	2,993,363.524	2,993,363.524	881.492
27	2,993,611.819	2,993,611.819	898.803

.

4.3. Results

After obtaining the arrival times of seismic waves from the sensors, noise signals are identified and removed by comparing their waveforms with those of various noise sources. Based on the characteristics of downhole noise sources and their waveforms, noise signals such as drilling impact, blasting vibrations, and electrical interference are accurately recognized and eliminated, thereby improving the quality of the signal set. The target area was divided into a grid, with 5 m × 5 m grids filling six section planes. Based on the parameter settings from Section 3’s synthetic experiment and the specific conditions of the target area, the maximum number of iterations, the residual tolerance to ceasing iteration, and the quasi-Newton step size are set to 20, 1 × 10⁶, and 1.0, respectively. The standard deviations on $\ln \left({\mathrm{V}}_h\right)$ and $\mathrm{E}$ are 0.1. The standard deviations on ${\mathrm{\sigma }}_{\mathrm{x}}$, ${\mathrm{\sigma }}_{\mathrm{y}}$, ${\mathrm{\sigma }}_{\mathrm{z}},$ and ${\mathrm{\sigma }}_{{\mathrm{t}}_0}$ are 0.2. The correlation length $\mathrm{\lambda }$ is 5. Using the goafs’ area and location based on the statistical analysis of historical data as the prior models, the P-wave arrival time data collected from monitoring were used to invert the velocity field for each plane, following the method described in Section 2.

After obtaining the wave velocity imaging results for each section plane, goaf areas are rapidly identified based on the differences in wave velocity between non-empty and empty anomalous areas. Subsequently, smooth curves connect grid units situated at the periphery of these mined-out zones, delineating their boundaries. Figure 10 presents the layout of the results for different sections based on the physical detection results. In these figures, the red area denotes non-extracted rock mass, while the white grid area indicates the goaf regions; the blue line delineates the boundaries of these goaf areas.

Table 2. Goaf area of different stopes derived from the statistical analysis of historical data and physical detection, with comparative analysis of results.

Section	The Serial Numbers of Stopes	The Goaf Area S1 Based on the Statistical Analysis of Historical Data (m2)	The Goaf Area S2 Based on the Results of Physical Detection (m2)	The Rate of Change, R, from S2 to S1 $\mathbf{R}={\displaystyle \frac{{\mathbf{S}}_2-{\mathbf{S}}_1}{{\mathbf{S}}_1}}\times 100\%$
First Section	68160	497.96	520.2328	4.47%
	68132	481.95	270.9463	−43.78%
	68134	507.22	559.0360	10.22%
	68136	553.93	500.1195	−9.71%
	68138	319.20	312.7774	−2.01%
Second Section	68232	821.18	835.1407	1.70%
	68234	633.92	552.5929	−12.83%
	68238W	491.07	549.2915	11.86%
	…	…	…	…
	45261E	365.28	369.5426	1.17%
	63234E	551.01	571.1484	3.65%
Third Section	68334	306.18	312.0997	1.93%
	68336	250.01	252.0154	0.80%
	68336E	462.81	483.1273	4.39%
	…	…	…	…
	45332E	287.97	304.2145	5.64%
	63334	407.58	417.9953	2.56%
Fourth Section	68434	322.47	325.5593	0.96%
	68436	421.5	429.3596	1.86%
	68436N	785.77	882.6283	12.33%
	…	…	…	…
	63442W	719.72	783.9989	8.93%
	63442E	651.69	705.7343	8.29%
Fifth Section	68534	249.84	251.9532	0.85%
	68536	323.14	328.0267	1.51%
	68538W	560.96	581.8787	3.73%
	…	…	…	…
	55540	294.79	301.6728	2.33%
	55540E	501.26	516.3566	3.01%
Sixth Section	68636E	810.43	844.0019	4.14%
	68638	536.79	547.8482	2.06%
	68640	496.57	505.3309	1.76%
	…	…	…	…
	55640	361.91	371.8363	2.74%
	55640E	307.06	314.1288	2.30%

presents the goaf area of different stopes derived from the statistical analysis of historical data and physical detection, along with a comparative analysis of these two results. Owing to spatial constraints, Table 2 presents only a subset of the data. A comprehensive dataset can be found in the Excel file titled ‘The Goaf Area of Various Stopes Derived from Statistical Analysis of Historical Data and Physical Detection.’

• Figure 10a is a planar representation of the physical detection results for the first section of the tantalum–niobium mine. Analysis of Figure 10a reveals five distinct goaf zones within this section. Historical records corroborate this observation, confirming that there are indeed five identified stopes in the first section, specifically numbered 68160, 68132, 68134, 68136, and 68138.
• Figure 10b represents the physical detection results for the second section of the tantalum–niobium mine. Statistical analysis reveals the presence of three mineralized veins in this section: Vein 68, which encompasses stopes 68232, 68234, 68234E, 68238W, 68238E, 68240, 68242W, 68242, 68244, 68246W, 68246, 68261_down, 68263E_down, 68263W_down, 68261^up, 68263E^up, and 68263W^up; Vein 45, which includes stopes such as 45265E, 5263N, 45263, 45263NE, 45261W, and 45261E; and Vein 63, which is represented by stope 63234E. Furthermore, statistical data indicate that there are a total of twenty-four voids distributed throughout this section.
• Figure 10c represents the physical detection results for the third section. The statistical analysis reveals three mineralized veins in this section: Vein 68, encompassing stopes 68334, 68336, 68336E, 68338, 68340, 68342W, 68342E, 68344, and 68346W; Vein 45, which includes stopes such as 45365S^up, 45365N^up, 45365S_down, 45365N_down, 45363S^up, 45363N^up, 45363S_down, 45363N_down, 45361S_down, 45361N_down, 45361S_down, 45361N_down, 45361E^up, 45361E_down, 45360^up, 45360_down, 45332, and 45332E; and, finally, the singularly identified Vein 63, containing stope 63334. Furthermore, statistical data indicate a total of thirty goafs distributed throughout the third section.
• Figure 10d represents the physical detection results for the fourth section. The statistical analysis reveals four distinct mineralized veins within this section: Vein 68, which encompasses stopes 68434, 68436, 68436N, 68438W, 68438, 68440, 68442W, 68442E, 68444, and 68444E; Vein 45, comprising stopes 45465E, 45463, 45461N, 45461S, 45461E, 45460, and 45432; the third vein, identified as the Vein 43, that includes stopes 43438W^up, 43438W_down, 43438E^up, 43438E_down, 43440^up, and 43440_down; and, finally, the designated area known as Vein 63, which consists of stopes 63440, 63442W, and 63442E. Furthermore, statistical data indicate a total of twenty-six goafs distributed throughout the fourth section.
• Figure 10e represents the physical detection results for the fifth section. The statistical analysis reveals three distinct mineralized veins within this section: Vein 68, comprising stopes 68534, 68536, 68538W, 68538, 68540, 68542W, 68542E, 68544, and 6546; Vein 63, including stopes 63538E, 63540, 63542, and 63542W; and Vein 55, encompassing stopes 55540 and 55540E. Furthermore, it is noted that a total of fifteen goafs are distributed throughout the fifth section.
• Figure 10f represents the physical detection results for the sixth section. The statistical analysis reveals the presence of three mineralized veins within this section: Vein 68, encompassing stopes 68636E, 68638, 68640, 68642S, 68642N, 68642E, 68644, 68646, and 68646W; Vein 63, with its associated stope 63642; and Vein 55, which includes stopes 55640 and 55640E. Furthermore, the statistical data indicate a total of twelve goafs distributed throughout the sixth section.

5. Discussion

Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 present the distribution of data, including the goaf area of different stopes derived from the statistical analysis of historical data and physical detection, along with a comparative analysis of these two results. Through a comprehensive analysis of the data and their distribution, the following can be ascertained:

1.

The physical detection result for the goaf area of Stope 68132 in the first section is 43.78% lower than the historical statistical data, while the discrepancy between the physical detection outcomes and historical statistics for goaf areas in the other four stopes is a maximum of 10.22%. Consequently, it can be inferred that Stope 68132 may contain a goaf, whereas it is confirmed that all four other stopes within the first section have established goafs.

2.

The physical detection results for the goaf areas of stopes in the second section—specifically, 68261_down, 68263E_down, 68263W_down, 68261^up, 68263E^up, and 68263W^up, as well as stopes 45265E and 45263N—indicate significant decreases relative to historical statistical data by percentages of 83.69%, 76.13%, 86.05%, 78.29%, 79.85%, 76.00%, and 80.39%, respectively. Conversely, the physical detection results for goaf areas in stopes, such as those designated as 68240 and 68242W, show increases of approximately 67.35% and 58.95%. Additionally, changes in the goaf area results for stopes identified as both 68234E and 68242 were around 30%, while alterations in the remaining twelve stopes did not exceed 12.83%. Consequently, it can be inferred that the eight specified stopes (i.e., 68261_down, 68263E_down, 68263W_down, 68261^up, 68263E^up, 68263W^up, and both 45265E and 45263N) within the second section may contain goafs; conversely, all the other sixteen stopes possess confirmed goafs.

3.

Based on data regarding the underground goaf areas of the third section illustrated in Figure 13, it can be observed that the physical detection result for the goaf area in stope 68342W is 81.95% lower than the historical statistical figure. In contrast, discrepancies between physical detection outcomes and historical statistics for goaf areas in other stopes do not exceed 10%. This suggests that while a goaf may exist within stope 68342W in the third section, other stopes within this section possess such goafs.

4.

The distribution of data from the fourth section indicates that the physical detection results for the goaf area in stope 43438E_down are 84.96% lower than historical statistical averages, while those for stope 68436N exceed historical averages by 12.33%. Furthermore, discrepancies between physical detection results and historical statistics for goaf areas in other stopes do not surpass 10%. Consequently, it can be inferred that a goaf likely exists in stope 43438E_down within the fourth section, and it is confirmed that goafs are present in other stopes of this section.

5.

Figure 15 shows that the discrepancy between the physical detection results and the historical statistical data regarding the goaf area of all 15 stopes in the fifth section is less than 10%. Notably, stope 63542W exhibits the greatest deviation, at 7.11%. This suggests that all 15 stopes within this section contain goafs.

6.

Based on the data distribution illustrated in Figure 16, the physical detection result for the goaf area in stope 68642S in the sixth section is found to be 86.56% lower than the historical statistical result. In contrast, discrepancies between physical detection outcomes and historical statistics for goaf areas in other stopes remain below 10%. Notably, stope 68646 exhibits the most significant change at a rate of 5.80%. Consequently, it can be inferred that a goaf may exist within stope 68642S in the sixth section, and the presence of goafs in the other stopes within this section is confirmed.

In summary, among the sixty-one stopes of Vein 68, nine stopes are suspected to contain goafs, specifically, stopes 68132, 68261_down, 68263E_down, 68263W_down, 68261^up, 68263E^up, 68263W^up, 68342W, and 68642S, and the remaining fifty-two stopes are confirmed to exhibit goafs. Furthermore, within the context of Vein 45, comprising a total of thirty-one stopes, two—stope 45265E and stope 45263N—are also suspected to have goafs; the remaining forty-three stopes are confirmed to exhibit goafs. Additionally, from the analysis of the six stopes associated with the vein under consideration, Vein 43, stope 43438E_down is indicated to potentially contain goafs, while the other five stopes are confirmed to possess them. Lastly, all ten stopes in relation to Vein 63 and four in connection with Vein 55 have been definitively established to contain goafs.

To verify and ensure the accuracy and reliability of the detection results, we calculated the relative error values by comparing the detection results of the goaf from synthetic and field experiments with their respective prior models. We further obtained the mean relative error, standard deviation of relative error, and the percentage of high-accuracy results (relative error < 10%), and plotted the distributions of these metrics across the synthetic experiment and field experiment (as shown in Figure 17).

From Figure 17, it can be observed that, at the experimental scale, both the mean relative error and standard deviation are low. When the wave velocity tomography method is applied at the practical engineering scale, the mean relative error and standard deviation increase slightly but remain within acceptable ranges, with the mean relative error below 10% and the standard deviation below 25%. This indicates that the goaf detection results from the field experiments exhibit low overall error and minimal fluctuation, yielding satisfactory results.

Examining the detection results of individual sections, the accuracy is highest in the third, fourth, and fifth sections, where the mean relative error is generally below 7%, and the percentage of high-accuracy results reaches 91–100%. In contrast, the detection accuracy in the second section is lower, with larger fluctuations in relative error. Based on the factors influencing the accuracy of anomalous boundary identification discussed in Section 3.1 and Section 3.2, combined with the characteristics of the goaf distribution in the second segment, it is inferred that the reduced accuracy may result from the irregular shape, the larger edge angles, and the scattered distribution of the goafs in this section. These factors lower the velocity inversion accuracy, thereby affecting the boundary identification precision of the goaf.

By integrating historical data with results derived from physical detection, a total of 12 potential stopes and 100 confirmed stopes were identified in the tantalum–niobium mine.

6. Conclusions

This study proposed the use of wave velocity tomography to identify abnormal structures and determine the boundaries of abnormal areas. The method was tested on both non-empty and empty areas, and the accuracy of the obtained velocity images was quantified. The method’s boundary recognition reliability was also analyzed. The conclusions are as follows:

• Twelve tests were conducted for four geometric dimensions multiplied by three velocity structures. The inversion of the velocity images showed that the inversion accuracy was significantly higher in normal areas than in abnormal areas. The size of the estimated anomalous structure in the prior model has a significant impact on the accuracy of the results. The larger the wave velocity difference between the anomalous and normal regions, the lower the inversion accuracy of the velocity image.
• By calculating the number of intersections of ray paths in nodes and identifying boundaries based on their density, it was determined that when the actual speed difference is small and the shape in the prior model is inaccurate, the identified boundary error is significantly larger. Among the 12 test results, only three of them had boundary recognition results with accuracy rates below 90%.
• Due to the unknown existence of abnormal areas in practical engineering contexts, we used a completely homogeneous prior model. The results showed that although the inversion velocity may not be accurate, it is sufficient to identify boundaries. The larger the velocity difference in the original model, the clearer the recognized boundary.
• Tests considering circle, square, and triangle areas showed that when the abnormal area is empty, the rounder the abnormal area, the higher the accuracy of the velocity field inversion and the easier it is to identify the boundary; the larger the edge angle of the abnormal area, the lower the inversion accuracy of the velocity field, and the more difficult it is to clearly describe the boundary.
• Physical detection of concealed goafs in the tantalum–niobium mine belonging to the Hunan Fuguihengtong Mining Company was conducted through wave velocity imaging which, when integrated with historical statistical data regarding goaf occurrence in each section of the mine, revealed a total of 112 identified or suspected goafs. Of these, 100 were confirmed as actual goafs and 12 remained suspected goafs. This finding underscores the efficacy of the proposed method for detecting empty anomalous areas, particularly in accurately identifying and delineating the boundaries of hidden goafs.
• By conducting error analysis on the goaf detection results obtained from synthetic and field experiments, and comparing the mean relative error, standard deviation, and the percentage of high-accuracy results across different scenarios, it is observed that the error values of goaf detection remain low and stable at both the experimental scale and the practical engineering scale. This demonstrates that the proposed method for detecting anomalous rock mass regions maintains high accuracy in practical applications. However, the detection accuracy in Section 2 is suboptimal, which is likely due to the irregular shape, large edge angles, and scattered distribution of the goafs in this section.
• Through both simulated experiments and field applications, this study demonstrates that the proposed method of using wave velocity tomography to identify anomalous structures offers advantages such as fast detection, wide coverage, and high accuracy. Compared to conventional methods for detecting rock mass anomalies, this approach has unique benefits: On one hand, as a non-destructive method for rock mass defect detection, it is cost-effective and more efficient than drilling-based techniques for obtaining structural information. On the other hand, unlike methods that rely on electromagnetic or resistivity measurements to infer rock properties and structures, this approach does not require survey line layout, making it more suitable for detecting deep rock mass anomalies and providing detailed geological imaging.