《1. Introduction》

1. Introduction

With the gradual depletion of land resources, the focus of geotechnical engineering has gradually moved toward underground construction [1,2]. However, buried abnormal regions in the rock-mass bring significant challenges to underground geotechnical engineering. The safety of underground geotechnical engineering is threatened by various potential underground disasters [3,4] such as rockburst [5,6], roof collapse [7], and water inrush [8]. Such potential hazard areas often differ from the rockmass in stress, causing great damage and casualties during underground excavation and exploration [9,10]. Since the stress field in geotechnical engineering is in constant redistribution [11–14], the spatiotemporal evolution pattern remains to be revealed. It is challenging to use current stress-field identification technology to detect potential hazard areas in complex rock-mass structures.

When the rock-mass is subjected to external stresses as well as high temperature and failures with extensive pre-failure damage, a significant velocity reduction exists [15,16]. Velocity field imaging technology can be utilized to detect the internal damage evolution of the rock-mass during the stable and unstable stages of crack propagation [17]. This technology can be further applied to obtain early warning of potential hazard sources in mining engineering, which is urgently needed in mining practices, in order to provide sufficient time for evacuation from potentially threatened areas [18] and guidance for underground autonomous driverless vehicles [19]. Velocity field imaging technology plays a significant role in the safety and efficiency of geotechnical engineering activities.

Identifying abnormal regions is one of the most classical and fundamental problems in underground geotechnical engineering [20,21]. In the 1970s, Aki et al. [22,23] pioneered the introduction of medical computed tomography into geophysical science. Tomography technology has been developed to explore the stratigraphic structure and geological disasters. The two-point raytracing method [24–26] used to be widely applied for natural seismic tomography. However, due to the low computational efficiency, tendency to track blind areas, and ease of falling into a local optimum, the method was phased out with the development of ray-tracing technology. Nakanishi and Yamaguchi [27] proposed the shortest-path method (SPM), inspired by the path problem of graph theory. In SPM, the region is discretized into units; then, travel time on the unit boundary is calculated, where the node with the minimum travel time is connected as the ray path. Nasr et al. [28] proposed a dynamic SPM (DSPM) and further improved the accuracy and efficiency by adding and removing nodes in the SPM topology. Dong et al. [29,30] developed an improved A* shortest-path algorithm combined with the match method to identify the empty region in a two-dimensional (2D) structure. Apart from the ray-tracing theory, the eikonal equation is an essential principle of tomography. The fast marching method (FMM) combined with narrowband technology describes the propagation mode of the wavefront, effectively determines the update sequence of grid points [31–33], and calculates the travel time of the region. Dong et al. [34] adopted FMM combining active and passive sources in synthesis experiments to investigate the influencing factors of acoustic tomography for complex structures. Jiang et al. [35] proposed an improved FMM method and realized more accurate location of rock fractures, thereby facilitating the detection of damaged regions in the rock-mass. Brantut [36] proposed an active–passive joint tomography inversion algorithm using FMM and applied it to sandstone compaction. However, the proposed active–passive joint tomography inversion algorithm failed to detect empty regions due to the excessive velocity gradient between the rock-mass and empty regions. The fast sweeping method (FSM) obtains a fast numerical solution to the eikonal equation through Gauss–Seidel iterative alternating scanning [37–39]. This method has efficient forward modeling and holds great significance for the engineering application of tomography. As a practical technique, tomography can quantitatively reveal the geological structure during large-scale seismic and oil–gas storage exploration [40–43]. The fracturing of intact rock accompanies microseismic activity as a result of the sudden release of accumulated strain energy around underground openings, and the redistributed stress and formed microfractures can be detected using microseismic signals [44–46]. Tomography is essential for detecting buried abnormal regions in complex mining environments.

This paper proposes the use of traveltime tomography combined with the damped least-squares QR-factorization (DLSQR) and Gaussian filtering (DLSQR–GF) to identify abnormal regions in complex rock-mass structures. The proposed method not only overcomes the limitation of velocity difference between rock and empty regions but also mitigates the impact from isolated velocity mutation in the iteration. SPM, DSPM, and FSM are used as the forward modeling in the traveltime tomography to obtain the arrivals in computational regions. Numerical and laboratory experiments have been conducted to evaluate the identification accuracy and computational efficiency of the proposed method. The proposed method has also been applied in the field experiments to verify its feasibility.

《2. Methods》

2. Methods

《2.1. Forward modeling》

2.1. Forward modeling

The forward modeling of tomography requires precise and fast calculation of arrivals. It is an essential requirement for the early identification of abnormal regions. Typically used methods include ray-tracing technology based on ray theory and wavefront-tracing technology using the eikonal equation. For efficiency, we adopted SPM, DSPM, and FSM as the forward modeling methods in the paper. The basic principles of these methods are briefly introduced below.

2.1.1. Shortest-path method

SPM is derived from the shortest-path problem in graph theory. The region is discretized into units, and nodes are set on the boundary of each unit. The adjacent nodes are calculated with the arrival from the source node and then set as the secondary sources for other nodes. Based on the minimum calculated arrival and the corresponding secondary source node, the connection between nodes is linked and set as the ray-tracing path.

2.1.2. Dynamic shortest-path method

To eliminate errors from the uniform slowness in units and the plane wave assumption, Nasr et al. [28] proposed DSPM to define the slowness in nodes. DSPM first calculates arrivals for all nodes by means of improved SPM and dynamically changes the slowness resolution by adding and removing secondary nodes in the iteration topology. The steepest traveltime gradient is used for each ray pair in the ray path-tracing process. The slowness values at the intersection areas of the ray path and cells are then used to recompute arrivals at each receiver. DSPM improves the computation accuracy without increasing mesh resolution, thereby reducing memory and computation cost; it can obtain better forward modeling of arrivals under an equal SPM computation cost.

2.1.3. Fast sweeping method

The FSM calculates the arrivals of region nodes through Gauss– Seidel iterative alternating scanning. First, the arrival at the source node is initialized to zero, and others are set at an extreme value. Then, the Gauss–Seidel iterative algorithm is used to calculate the solution of the eikonal equation, and the minimum value among the original value and updated value is retained. The scanning process continues until it satisfies the convergence condition. The FSM has a considerably higher computing efficiency than the FMM and far more than the ray-tracing approach. The forward calculation is more efficient with FSM. In tomography, the forward calculation is always the most time-consuming part; thus, FSM is of great significance to the application of tomography.

《2.2. Inversion algorithm》

2.2. Inversion algorithm

The inversion problem of tomography can be solved as the optimization problem composed by observation arrivals and calculated arrivals from the forward modeling, which is expressed in Eqs. (1)–(3): 

where dobs represents the observed arrivals, dcal is the calculated arrivals, mi is the slowness field at i iteration, and dm is the slowness field difference between mi and mi+1. The ray-tracing path L is a nonlinear sparse operator containing the grid distance between the sources and the sensors; it is calculated from the forward modeling calculation, and then combined with the slowness field m to obtain the calculated arrivals dcal in the inversed process. It establishes the relationship between the observed arrivals d and the slowness field m in the inversion region. Since the misfit between the observed and calculated arrivals is mainly caused by the ray paths in different slowness fields, the prior slowness field can be optimized toward the actual slowness field by minimizing the misfit function.

Eq. (3) is an underdetermined and contradictory equation, which can be solved by iterative inversion methods such as the back projection method, the Newton method, and DLSQR. The latter has good convergence and reliability in solving ill-posed linear equations with a low computation cost and can thus be used to avoid the calculation of gradient matrix; as a result, it is only necessary to calculate the non-zero elements of the large sparse matrix. DLSQR simplifies the calculation with the sparsity of the matrix and demands less storage space and computation cost. It is suitable for solving a nonlinear sparse operator L composed of ray paths. DLSQR can be divided into several steps: transforming L into a tridiagonal matrix by a Lanczos vector, constructing the least-squares equation, generating the lower triangular matrix, and solving it. However, DLSQR cannot avoid isolated velocity mutation, which may accumulate errors in the inversion process. We employ a Gaussian filter to mitigate salt-and-pepper noise from the isolated velocity mutation. 

The traveltime tomography inversion is composed of the following steps:

(1) Establishing the initial model and dividing mesh grids. The size and number of grids affect the resolution of the results; increased computation costs and a greater number of signals are required with finer grids.

(2) Collecting signal data, including the arrivals and coordinates of source-receiver sensors and lead breaking points. The collected signal date are performed with the noise reduction algorithm to ensure the accuracy of arrivals.

(3) Executing tomography combing forward modeling and DLSQR–GF. The calculated arrivals from the forward modeling part are optimized with the observed arrivals by DLSQR–GF to obtain the velocity field during the iteration.

(4) Stopping the iteration when the convergence criterion is achieved. The extended information criterion (EIC) can determine the optimum number of iterations using bootstrap statistics [47], avoiding over-iteration caused by using the residual error as the convergence criterion. The EIC and the maximum iteration can be applied as the convergence criterion.

(5) Identifying abnormal regions from the traveltime tomography results. The abnormal regions are clustered by the obviously different velocity cells compared with the background velocity. The ability of identifying abnormal regions can be evaluated by the statistical method.

The traveltime tomography procedure is shown in Fig. 1.

《Fig. 1》

Fig. 1. Flowchart of the tomography inversion algorithm. SPM, DSPM, and FSM are adopted as forward modeling methods, and DLSQR–GF is proposed as the inversion algorithm.

《3. Experiments》

3. Experiments

In order to quantitatively evaluate the performance of traveltime tomography in the early identification of abnormal regions in complex rock-mass structures, numerical, laboratory, and field experiments were conducted. The numerical experiment was based on a 50 m × 40 m × 20 m excavation model containing two elliptic low-velocity abnormal regions. Corners were respectively arranged with receivers with equal intervals. The laboratory experiment consists of 2D and three-dimensional (3D) experiments. The region of interest is a 500 mm × 200 mm × 160 mm cuboid granite sample. There are five boreholes with a diameter (D) of 50 mm at the top, representing low-velocity abnormal regions. The parameters of the numerical and laboratory experiments are shown in Table 1. The positions of the sensors and abnormal regions are shown in Figs. 2(a) and (b). The configuration of sensors in the laboratory experiments was arranged to increase the ray coverage in the rock-mass as much as possible. Disturbance of the slowness field origin from the ray-tracing path and sufficient ray coverage can improve the identification accuracy of abnormal regions in the rock-mass. Figs. 2(c) and (d) depicts the grid meshes in the laboratory experiments. Pulse signals and lead break experiments were conducted in the laboratory experiments as active sources. The basic frequency of the pulse signals is around 160 kHz. By establishing an analytical equation with coordinates, the occurrence time of the lead breaking points was determined using the maximum likelihood estimation method. AMSY-6 multi-channel equipment and VS45-H sensors were used in the laboratory experiment to sample signals with a high standard. The frequency range of the sensors is 20–450 kHz, and the frequency response reaches a peak at 280 kHz.

《Table 1 》

Table 1 Parameters of the numerical and laboratory experiments.

h: height.

《Fig. 2》

Fig. 2. Configuration of sensors and grid meshes in laboratory experiments. (a) 2D sensor network; (b) 3D sensor network; (c) 2D grid meshes; (d) 3D grid meshes.

For early identification, it is necessary to rapidly distinguish abnormal regions from rock-mass. We use the identify rate (IR) and cost rate (CR) to quantitatively evaluate the identification efficiency of traveltime tomography. Effective identification can be said to be achieved when the velocity of abnormal regions is significantly lower than the background velocity. We set the IR as the proportion of the velocity in abnormal regions that is lower than the lower quartile of the background velocity. The lower quartile value is calculated from the mesh grids in non-cavity regions and is then compared with the mesh grid in cavity regions. The IR value can be obtained from the ratio of the number of mesh grids with higher velocity than the lower quartile velocity value in non-cavity regions to all mesh grids in cavity regions. A higher IR value represents the higher identification degree of abnormal regions.

《3.1. Numerical experiments》

3.1. Numerical experiments

To evaluate the practical usefulness of the proposed method under theoretical conditions, we constructed numerical experiments, as shown in Fig. 3. The numerical model is designed to simulate actual abnormal region detection during excavation. The target region has the dimensions 50 m × 40 m × 20 m and is surrounded by other rock-masses. Two different abnormal regions are buried in the target rock, which may affect the excavation process. To detect the abnormal region distribution, there are boreholes set on corners of the target rock-mass, and sensors are arranged at 5 m intervals along the boreholes. The sensor arrangement constructs a monitoring network for the target rock-mass. The parameterization of the target rock-mass is shown in Fig. 3(b). The background velocity is set as Vout = –4500 m·s–1, and the abnormal low-velocity regions are respectively set as 1500 and 340 m·s–1. The velocity of the excavation face region is set as 340 m·s–1. Imaging processes and results using the SPM, DPSM, and FSM of numerical experiments are shown in Fig. 4. The initial model is composed of the target regions with the background velocity and the excavation face region, as shown in Fig. 4(a).

《Fig. 3》

Fig. 3. (a) Layout of abnormal region detection during excavation and (b) its velocity parameterization model. (c–e) Velocity parameterization models of the imaging results from (c) SPM, (d) DSPM, and (e) FSM.

《Fig. 4》

Fig. 4. (a) Initial model and (b) real model: (c–h) Imaging results: (c, e, g) imaging results of (c) SPM, (e) DSPM, and (g) FSM in the 3D numerical experiments; (d, f, h) imaging results of (d) SPM, (f) DSPM, and (h) FSM with 5% arrival errors.

The tomography results of SPM, DSPM, and FSM are shown in Figs. 3(c)–(e). It can be seen from the horizontal velocity profile results that the low-velocity cells are located around the actual abnormal regions. The lower velocity abnormal regions are near the boundary of the target rock-mass, and the other abnormal regions are located near the tunneling face. The inversed model is closed to the designed numerical model. However, the tomography inversion must solve an undetermined, multi-solution problem, whose degree of underdeterminedness is directly related to the signals and grids of the slowness field. Thus, the inversed model may differ from the real model in terms of the shape and values of abnormal regions. The velocity outside of the abnormal region also shows various degrees of change in the tomography results. The cells around low-velocity regions are assigned higher velocity values than the background velocity in order to minimize the arrival differences from the different shape and velocity values of abnormal regions.

Figs. 4(c), (e), and (g) show the perspective views of the initial model and the tomography results. For a better demonstration of abnormal region distribution, the inversed models are sliced along the height. The tomography process begins from the initial model, which only contains the excavation face. The tomography results are obtained by combining forward modeling with DLSQR–GF. The tomography results clearly show the distribution of the excavation face and abnormal regions; thus, the proposed methods realize the detection of abnormal regions. The IR values of SPM, DSPM, and FSM are 93.96%, 94.02%, and 94.16%, respectively. Although there is little difference between the IR values in the numerical experiments using these three methods, FSM exhibits a relatively better performance. With the computation cost of SPM as the evaluation criteria, the CR values of DSPM and FSM are 28.92% and 18.13%, respectively. The computation cost of FSM is significantly lower than those of the others, indicating higher computational efficiency. To verify the robustness of the proposed method against noise, the average 5% arrival errors are added to the theoretical travel times. The tomography results are shown in Figs. 4(d), (f), and (h). The results show that, although the tomography accuracy is affected to some degree in this case, the proposed method can still realize abnormal region detection.

In the 3D tomography numerical experiments, DSPM and FSM demonstrated excellent identification ability and computation efficiency in abnormal region detection. Although SPM satisfied the requirement of identifying abnormal regions, its computation cost is significantly higher than those of other methods.

《3.2. Laboratory experiments》

3.2. Laboratory experiments

In contrast to the ideal conditions in the numerical experiment, error factors are unavoidable in the tomography laboratory experiments, including arrival errors, uniform background velocity, and sensor–media coupling. The accuracy of the tomography is disturbed by these factors. Therefore, it is necessary to verify the application of traveltime tomography in laboratory experiments. We carried out a quantitative evaluation of the identification ability of traveltime tomography for abnormal cavity regions in a rock sample.

3.2.1. 2D tomography laboratory experiments

The configuration of the sensors is shown in Fig. 2(a). Sensors on both sides are used as sources and receivers, respectively. The background velocity Vout = 4500 m·s–1 is set as the initial model, and the tomography results of SPM, DSPM, and FSM are shown in Figs. 5(a)–(d). SPM shows undesirable imaging results in the laboratory experiments. This is due to the more complex cavity structure in the laboratory experiments in comparison with the numerical experiments, which increases the difficulty of the tomography; moreover, the result is affected by the uniform velocity distribution and outliers in the imaging process. Artifacts appear around the sources, receivers, and boundary regions in the SPM tomography result. Thus, it is difficult to effectively identify the distribution of low-velocity abnormal regions in SPM. The results from DSPM and FSM are less impacted by error factors in the laboratory experiments. There are fewer artifacts around the sources and receivers and more obvious diffraction phenomena than in SPM. An accuracy comparison of the initial model and the methods for identifying abnormal cavity regions is provided in Figs. 5(e)–(h). The imaging result from SPM is undesirable, and the surrounding artifacts significantly affect the identification. The identification abilities of DSPM and FSM are relatively similar. FSM has a better distribution in the middle low-velocity abnormal region than DSPM. The identification of the low-velocity region from DSPM and FSM matches the actual distribution of the abnormal region. The IR values of SPM, DSPM, and FSM are 65.83%, 74.59%, and 72.86%, respectively. The IR value of SPM is affected by artifacts outside of the actual cavity regions in the tomography results. There is little difference between the IR values of DSPM and FSM.

《Fig. 5》

Fig. 5. (a–d) Imaging and ray-tracing results of (a) the initial model, (b) SPM, (c) DSPM, and (d) FSM. (e–h) Comparison of the abnormal region distribution between (e) the initial model, (f) SPM, (g) DSPM, and (h) FSM.

The results show a low-velocity band between two left abnormal cavity regions and extended low-velocity zones around other abnormal cavity regions. This may be caused by improper operation in the drilling process. Although the internal damage cannot be seen from the visible break on the surface of the rock sample, it affects the velocity structure, as shown in tomography results. The low-velocity artifacts near the boundaries may be due to insufficient ray coverage and the stability of the forward modeling method. The inversion of the wave velocity field is an underdetermination problem, which may induce artifacts when there is insufficient information or errors in the inversion process.

The tomography laboratory experiments are affected by error factors, including arrival errors, uniform background velocity, and sensor–media coupling. These error factors increased the application difficulty of the early identification of abnormal regions. SPM is severely affected, and the inversed result is undesirable compared with those of the other methods. The identification abilities of DSPM and FSM are similar, with DSPM exhibiting a better performance in terms of computation efficiency.

《Fig. 6》

Fig. 6. Imaging results in 3D tomography laboratory experiments. (a) Initial model, (b) SPM, (c) DSPM, and (d) FSM: (i) are horizontal slices and (ii) are vertical slices.

3.2.2. 3D tomography laboratory experiments

The configuration of the sensors is shown in Fig. 2(b). Sensors are used as the receivers, while lead breaking points are used as the sources. Horizontal and vertical slices of the SPM, DSPM, and FSM tomography results are shown in Fig. 6. The low-velocity regions differ significantly from the background velocity and correspond to the location of cavities. However, 3D tomography experiments have a higher demand for signal quantity and are more vulnerable to error factors than 2D tomography experiments. This aggravates the difficulty of detecting 3D abnormal regions and induces artifacts in regions that the ray path cannot cover. A large low-velocity area appears on the top horizontal slice of the velocity field, and the actual distribution of cavities cannot be identified accurately in SPM. The tomography results of DSPM and FSM have few artifacts, while only the shallow horizontal slices at the top can be updated.

The increasing number of 3D grids demands more effective signals than 2D grids. However, the rays from broken leads cannot fully cover the cavities, so it is difficult to identify cavities effectively under a limited signal. Moreover, it is because the ray path changes in the 3D structure are more complex, increasing the difficulty of the tomography process. Due to insufficient coverage of abnormal regions, the low-velocity band fuzzes the area between each abnormal region. In addition, the concentrated low-velocity regions between holes in the middle and on the left sides may be caused by internal damage to the rock caused by an improper drilling process. The IRs of these methods are all about 0.32, subject to limited signal quantity.

Compared with 2D tomography laboratory experiments, 3D tomography experiments can provide a more comprehensive perspective for analyzing the velocity field. However, due to the more complex ray-tracing paths, the requirements for signals for 3D tomography are more stringent.

《3.3. Field experiments》

3.3. Field experiments

To verify the feasibility of the proposed method in engineering, the proposed method was applied in the Shaanxi Zhen’ao mine. As shown in Fig. 7, the installed microseismic monitoring system is distributed from level 795 to level 860, and is composed of 26 single-component sensors. The signals date used for tomography are multi-source acoustic events from April 1 to 4 in 2022. The target regions can be divided into two parts: The first part includes level 795-1 and level 795-2-810 (level A), and the second includes level 795-4-842 and level 860-1 (level B). The size of the established model is approximately 300 m × 520 m × 200 m, according to the blast location results and the sensor distribution. Based on a comparison of the tomography performance from the numerical and laboratory experiments, DSPM and FSM are applied in the field experiments. The double-difference method is used in the inversion part in order to eliminate the error impact of the occurrence.

《Fig. 7》

Fig. 7. The underground roadway distribution and the microseismic monitoring system in the Shaanxi Zhen’ao mine.

The initial models of levels A and B shown in Figs. 8(a) and (b) are both set as 4500 m·s–1 , corresponding to the velocity of the rock-mass. The imaging results from DSPM and FSM are demonstrated in Figs. 8(c)–(f), where the low-velocity regions correspond to the underground roadway distribution. The results demonstrate that the proposed method can differentiate the underground roadway from the rock-mass, which is equivalent to abnormal region detection during excavation. The imaging results have highvelocity abnormal regions, which are mainly distributed near the sensors and blasting points. On the one hand, these high-velocity regions might be caused by the large velocity disturbance raised by massive rays. The velocity near sensors and blasting points tends to change significantly. On the other hand, the highvelocity regions in the excavation face might be due to deformation and velocity variation caused by stress change. The imaging results from DSPM and FSM show little difference from the results of the field experiments.

Although field experiments are affected by the complex rockmass structure and error factors, the proposed methods still realize the reconstruction of the underground roadway from within the rock-mass. If more diverse data and a more realistic initial model are provided, the inversed results can be further improved.

《Fig. 8》

Fig. 8. Horizontal slices of the initial model at (a) level A and (b) level B; horizontal slices of imaging results from DSPM at (c) level A and (d) level B; horizontal slices of imaging results from FSM at (e) level A and (f) level B.

《4. Conclusion》

4. Conclusion

This paper proposes a traveltime tomography method to identify buried abnormal regions in complex rock-mass structures. The proposed method obtains arrivals in computational regions by means of forward modeling, including SPM, DSPM, and FSM. Next, residuals between observation and calculation arrivals are iteratively optimized with DLSQR–GF, which overcomes the limitation of velocity difference in empty-region detection and mitigates the impact from isolated velocity mutation during the optimization. Various numerical and laboratory experiments have been carried out to quantitatively evaluate the identification accuracy and computational efficiency of forward modeling methods with DLSQR-GF. The reliability and effectiveness of the proposed method with DSPM and FSM are demonstrated by IR and CR results on the detected area with complex abnormal regions. Field experiments were conducted in the Shaanxi Zhen’ao mine to verify the feasibility of the proposed method. The results show that this method can realize the reconstruction of an underground roadway and high-stress regions in the rock-mass. It not only utilizes characteristics of classified events from active blasting, drilling, and other mining activities, but also combines wide distribution of passive events induced by geo-stress changes. The proposed velocity structure imaging method using multi-source acoustic events ensures the huge and diverse data requirements for the highprecision mine structure identification. Therefore, the proposed method meets the requirements of accuracy and efficiency in the early identification of abnormal regions.

The distribution of sensors and sources affectes the accuracy of abnormal region detection in the experiments. The use of passive sources would be more reliable and effective for 3D abnormal region identification. In addition, the identification results can be further improved with a more reasonable sensor layout and initial model. The proposed method provides a theoretical foundation and technical support for application in detecting material heterogeneity and potential abnormal areas in underground geotechnical engineering, such as tunneling, mining, and so on.

《Acknowledgments》

Acknowledgments

We are grateful for the financial support from National Key Research and Development Program of China (2021YFC2900500) and Funds for International Cooperation and Exchange of the National Natural Science Foundation of China (52161135301).

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Longjun Dong, Zhongwei Pei, Xin Xie, Yihan Zhang, and Xianhang Yan declare that they have no conflict of interest or financial conflicts to disclose.