2.1 Forward modeling

In recent years, numerous forward modeling methods based on eikonal solvers have emerged [41], with FSM being one of the most popular. FSM was first proposed by ZHAO [42] based on rectangular grids as a tracing method for solving numerical problems, such as computing the distance function. QIAN et al [43] subsequently proposed a novel ordering strategy to enhance FSM for unstructured grids. The forward modeling based on FSM can compute first arrivals with more efficient and stable performance than other ray tracing methods, making it suitable for non-uniform velocity fields with drastic variations, such as potential hazard region detection. Therefore, we apply the unstructured FSM [43, 44] as the forward modeling method to calculate the first arrivals for irregular opencast mining structures.

The interested region is preprocessed by an open-source software Gmsh, and triangular mesh nodes are obtained, noted as . Next, the source nodes need to be located, and the shooting method obtains traveltimes of surrounding nodes, which remain unchanged in subsequent iterations. The sweeping order of FSM follows a specific direction when solving the eikonal equation and matters in the convergence of FSM. The sweeping direction for the two-dimensional rectangular coordinate system can be divided into four parts: up-right, up-left, down-right, and down-left. However, the order of the unstructured meshes is different [44]. It needs to set noncollinear points in the region as the reference points [45], and other points are sorted in ascending and descending orders based on their distances. Then, traveltimes are iteratively updated based on defined orders, and the minimum value is retained among the original and updated values. The iteration process continues until the convergence requirement is satisfied.

2.2 Backward inversion

As the actual model differs from the prior model, the further optimization of the prior model is needed toward the actual model through traveltime differences at the sensor nodes. The mine environment is often accompanied by blasting, rock drilling, and microseismic signals during mining operations. It provides ample data for identifying unknown underground structures. However, data on the opencast mining area is predominantly composed of microseismic signals, and the relatively small amount of data exacerbates inversion ill-posedness, posing challenges to opencast mining structure detection. To mitigate the ill-posedness resulting from insufficient data and source-sensor pairs, the damping regularization parameter and smoothing regularization parameter are introduced when solving the inversion equation by the least square method, expressed as follows:

(1)

where L represents the sparse raytracing path containing the grid distance between sources and sensors; and are the damping regularization parameter and smoothing regularization parameter, respectively; S represents the smoothing relationship between model disturbance and its surrounding model disturbance; represents the traveltime differences between calculated data and observed data. In practice, it is necessary to continuously attempt to select acceptable regularization parameters and . After obtaining the model disturbance, it needs to update the prior model. However, the updated model may not necessarily fit the observed data well. The updated model must be set as the prior model and iterated repeatedly until it satisfies the convergence criterion.

3.1 Numerical experiments

The numerical experiments can simulate the propagation process of acoustic waves between rock, soil, and potential hazard regions (such as cavities or water). This helps us understand the interaction between acoustic waves and the opencast mining area medium and the influence of various factors on the acoustic imaging results. Additionally, based on the numerical experiment results, we can optimize parameter settings of the forward modeling method and regularization operators to enhance the accuracy and resolution of detecting potential hazard regions.

The numerical experiments are based on two-dimensional unstructured models, which size 50 cm×50 cm and are arranged with different circular empty abnormal regions (R=5 cm, 10 cm, 15 cm, 20 cm). The real models are depicted in Figures 1 and 2, with left and right white nodes representing sources and sensors, respectively. The background velocity is set as the initial model 4000 m/s and the empty abnormal region velocity is set 340 m/s. Empty abnormal regions are located at the center of these models. The tomography inversion process involves 15 events and 225 rays. The amount of data is far less than the number of velocity nodes in the model. The limited number and distribution of active sources in opencast mining area monitoring result in insufficient rays. This exacerbates the ill-posedness of the inversion equation, presenting a challenge for detecting abnormal regions in the opencast mining area. Inversion results are compared to validate the performance of the proposed method for various potential abnormal regions.

Figure 1全屏查看

True model and inversion result of circular empty abnormal regions R=5 cm (a, b), and R=10 cm (c, d)

Figure 2全屏查看

True model and inversion result of circular empty abnormal regions R=15 cm (a, b), and R=20 cm (c, d)

Figure 1(b) shows inversion results with the empty abnormal region R=5 cm. A region can be observed forming and clustering in the center of the background velocity, with a velocity value relatively smaller than other regions. This corresponds to the empty abnormal region marked by a black dotted line. However, the value of the formed low-velocity region does not match the empty abnormal region velocity. The tomography inversion refers to solving the undetermined problem, where differences between observed data and model parameters result in multiple solutions to the inversion problem. Thus, the tomography technique aims to approximate the real model, and the obtained inversion results satisfy the requirement for abnormal region identification.

Furthermore, the slight decrease in velocity in the marked regions may be attributed to small differences in data originating from the abnormal region. The abnormal regions only affect rays passing through them in the prior model, thus inducing receiver data differences. The insufficient differences between the prior and real models lead to a slight degree of model update. However, this demonstrates that the current prior model has satisfied the preliminary identification of the abnormal region, and there is no need for further iterative updates. Thus, the geometric differences of abnormal regions between the prior and real models are crucial factors affecting the identification performance. The inversion result of the empty abnormal region R=10 cm is shown in Figure 1(d). The inversion result reconstructs the distribution of the abnormal region, which is close to the marked region. Compared with the empty abnormal region R=5 cm, there are more obvious diffraction phenomena around the empty abnormal region R=10 cm, and it also represents more data differences. Besides, a dark blue region is in the center of the low-velocity region, representing a sharper velocity decrease compared to the surrounding light blue regions. These blue regions enable the reconstruction of the marked empty abnormal region, and these unstructured meshes accurately identify the circular range. Figure 2 demonstrates the actual model and inversion result of circular empty abnormal regions R=15 cm and R=20 cm. As the size of the empty abnormal region changes, the scale of the low-velocity region in the inversion results becomes more extensive. This maintains consistency with the unstructured marked region. There are also some velocity decreases around sources and sensors. These low-velocity artifacts around sources and sensors can be caused by dense rays, which are sensitive to the model perturbation and tend to accumulate errors during the iteration. In general, their size directly influences the accuracy of identifying abnormal regions. If the anomalous region is too small, there will not be enough rays passing through it to adequately solve the inversion equation for the traveltime perturbation and reconstruct the distribution of the anomalous region. Conversely, when the abnormal region is too large, causing almost all rays to be affected and diffract, it will introduce more perturbation from arrival differences. As the equations are solved under-determined, errors are inevitably introduced, leading to artifacts around the inversion results. Thus, the design of the sensor layout must ensure an appropriate amount of traveltime perturbation. The inversion results have validated the feasibility of the proposed method for identifying unstructured abnormal regions under the theoretical two-dimensional condition.

3.2 Laboratory experiments

The laboratory experiments allow for the direct observation of the interaction between acoustic waves and rock samples in a controlled environment, facilitating the verification of the proposed method’s reliability. In practice, the propagation of acoustic waves is inevitably influenced by the coupling between sources and the medium, between different components of the medium, and between the receivers and the medium. These interactions introduce errors and effects. Therefore, applying the proposed method in laboratory experiments and comparing it with the results of numerical simulations can verify its accuracy. This verification can guide the design and interpretation of subsequent field experiments.

Unlike the theoretical condition in the numerical experiments, error factors in traveltime errors, the prior model, and sensor coupling are unavoidable in laboratory experiments. It is essential to verify whether the proposed method can identify unstructured abnormal regions in laboratory settings. The laboratory experiments comprise two parts: the first part (referred to as A) aims to verify the performance in identifying single empty abnormal regions, and the second part (referred to as B1 and B2) focuses on identifying complex multiple empty abnormal regions.

The initial model and its model mesh of experiment A are shown in Figure 3. The background velocity V=4500 m/s is set as the initial model, and sensors are placed on both sides as sources and receivers. Figure 4 shows the experiment layouts and inverted models of circular empty abnormal regions R=5 cm, and R=10 cm. For the circular empty abnormal region R=5 cm, some low-velocity cells are gathering around the central region and sensors. The velocity of the center region is lower than others, and it corresponds to the actual distribution of the empty abnormal region. The dense rays similar to the situation in the numerical experiments, might cause the reason for the form of low-velocity regions around sensors. The insufficient data differences for the circular empty abnormal region R=5 cm can hardly provide an identification result with higher accuracy.

Figure 3全屏查看

Initial model (a) and model mesh (b) of experiment A. Red nodes and black nodes represent sources and receivers, respectively.

Figure 4全屏查看

The experiment layouts and inverted models of circular empty abnormal regions R=5 cm (a, b), and R=10 cm (c, d)

Besides, the unavoidable error factors in the laboratory situation also increase the difficulty of tomography inversion. The inverted model of the circular empty abnormal region R=10 cm shows a more obvious distribution of abnormal regions. As the size of the actual empty abnormal region becomes more extensive, the low-velocity region from the inverted models in Figure 5 becomes larger accordingly. Though the low-velocity region is elliptical due to the insufficient ray coverage on the top and bottom sides, the lateral distribution of the low-velocity region is precisely located at the marked regions.

Figure 5全屏查看

Experiment layouts and inverted models of circular empty abnormal regions R=15 cm (a, b), and R=20 cm (c, d)

The buried unstructured abnormal regions in opencast mining areas might be discrete rather than clustered, and it is necessary to validate the identification capability of the proposed method for multiple abnormal regions. The experiment layout of B1 is shown in Figure 6(a), and there are two different shapes of cavities representing potential low-velocity regions.

Figure 6全屏查看

Experiment layout (a) and inverted model (b) of experiment B1

The sensors at both sides are excited in turn as the source, and the corresponding sensors on the other side are used as the receiving ends. The amount of travel time data collected is 288. The inverted model in Figure 6(b) shows that the empty regions’ approximate distribution areas are consistent with low-velocity regions. Although there are some low-velocity regions around the actual cavity area due to the uncertainty of the inversion process, the velocity values inside cavities are lower than those around surrounding cavities. Therefore, the distribution of the dark blue regions with lower velocity is obviously closer to the actual cavity area. It shows that in practical applications, a specific speed threshold can be set to better delineate potential abnormal areas. Figures 7(a) and (b) demonstrate the experiment layout and the model mesh of B2. There are five equal-sized circular empty regions in the marble blocks, and sensors are placed on the front and back to form a two-dimensional monitoring range. Sensors on the front utilize the pulse signals as active sources for receivers on the back. Amount 144 rays are participating in the identification of the multiple abnormal regions. The background velocity v=4500 m/s is set as the initial model, as shown in Figure 7(c). Moreover, it can be seen from the inverted model in Figure 7(d) that five empty abnormal regions are reconstructed clearly by the proposed method from the initial model. The distribution of empty abnormal regions accurately matches the marked regions. Though there are some artifacts in the inverted model, the value of empty abnormal regions is lower than those of artifacts.

Figure 7全屏查看

Experiment layout (a), model mesh (b), initial model (c), and inverted model (d) of experiment B2. Red nodes and black nodes represent sources and receivers, respectively.

The laboratory experiments are inevitably influenced by error factors, posing significant challenges to the application of the tomography technique. Low-velocity regions around the source and receiver regions may affect the identification results. This is caused by intense rays traveling from sources to receivers, inevitably bringing a large number of model perturbations. Although these model perturbations are numerically small, they gradually accumulate with iterations, eventually aggregating to form low-speed artifact regions. The results have validated the feasibility of the proposed method for identifying unstructured abnormal regions at the laboratory scale. The proposed method demonstrates excellent performance in reconstructing abnormal regions using sparse rays across a large number of meshes.

3.3 Field experiments

The field experiments can verify the applicability of the results of numerical and laboratory experiments in real scenarios and identify potential anomalous areas. Field experiments were conducted on an unmined rare earth mine slope in Zhongshan-Fuchuan County. According to the circumstances of the rare earth mine slope, 4 monitoring stations and 16 sensors are installed, as shown in Figures 8(a) and (b). Each sensor is equipped with a wave-guide rod for receiving wave propagation better. The rare earth mine slope sizes 45 m×50 m×30 m, and its meshes are generated by Gmsh in Figure 8(c). The hammering points were scattered on the rare earth mine slope, and the waves emitted were used as microseismic sources. According to the preliminary geological data survey, the rare earth mine slope belongs to the fully covered slope. There are three strata from top to bottom: clay bed, completely decomposed granite, and semi-weathered granite, as shown in Figure 8(d). Considering the velocity differences between strata, the wave path tends to pass through with higher velocity under the diffraction effect. Thus, the feasibility of the proposed method can be validated by reconstructing the geological structure of the rare earth mine slope from the initial model. The installed microseismic monitoring system and hammering operations work as receivers and active sources for the potential abnormal region identification and stability analysis.

Figure 8全屏查看

(a) Layout of the microseismic monitoring system; Distribution of sensors and sources (b), unstructured meshes (c) and strata (d) of the rare earth mine slope

The inverted results of the rare earth mine slope are shown in Figure 9. Different from the compact structure of the rock mass, the soft soil occupies the majority of the mine slope. The wave propagation needs to pass across air, oil, gravel, and the interface between them. The diffraction and refraction phenomena greatly increase the complexity of the ray path and the time required, which leads to the calculated average velocity based on received arrivals being much lower than the velocity in the complete medium. Thus, the tomography inversion for the mine slope is based on the relative velocity changes. The identification region mainly focuses on the center of the mine slope, corresponding to the layout of the microseismic monitoring system. The scale of the relative low-velocity region reaches the maximum when Y=31 m. It is inclined along the mine slope, and the distribution extends approximately from the middle crest to the left side. Furthermore, from the velocity slice Y=31 m in Figure 9(c), meshes can be divided into three strata according to the velocity value. The A strata with the lowest relative velocity is mainly distributed along the surface, corresponding to the clay bed of the rare earth mine slope. The B strata has a relatively higher velocity and distribute closely along the clay bed, which can be considered completely decomposed granite. The deep C strata cannot be deduced constrained by the shallow rays, and it might be improved when provided with the crosshole signals. However, the rare earth mine slope is still unexcavated, and the crosshole blasting operation cannot be carried out.

Figure 9全屏查看

Inverted result (a) and rays (b); the slice Y=31 m of the inverted result (c) and rays (d); the slice X=20 m of the inverted result (e) and rays (f)

When the distances between sources and receivers are far enough, the diffraction will happen that rays pass the B strata from the A strata to receivers, as shown in Figure 9(d). Figures 9(e) and (f) demonstrate the velocity slice and rays of X=20 m. With the decrease in height, the relative propagation velocity decreases gradually and then increases. The distribution angle of low-velocity regions is consistent with the inclination angle of the mine slope and is on the left side of the mine slope. These regions are close to the shallow layer, mainly between the clay bed and part of the completely decomposed granite, consistent with the actual situation. Due to the soil of the clay bed and the completely decomposed granite being relatively loose, the velocity between the two layers is relatively small. It also explains that the wave velocity in the mine slope is lower than in the general rock mass.

During the field experiments, the proposed method is implemented using hammering points as active sources on the surface. Results show that the velocity distribution aligns with the strata of the rare earth mine slope. Although field results may not precisely define the actual distribution of each stratum due to the absorption of acoustic waves by the soil medium, this paper supplements the technological capabilities of mine slope monitoring systems with acoustic tomography. It provides the wave velocity field as an identification index for potential hazards in the traditional monitoring of mine slopes. The application of this method is of significant importance for ensuring the safety of mine slopes.