Anomalous Areas Detection in Rocks Using Time-Difference Adjoint Tomography

: Detecting anomalous areas (such as caves, faults, and weathered layers) in rocks is essential for the safety of facilities and personnel in subsurface engineering. Seismic tomography has been proved to be an effective exploration technology in engineering geophysics. However, the complexity, anisotropy, and uncertainty in rock environments pose challenges to the resolution and robustness of tomography methods. Traditional tomography methods have difﬁculty balancing reliability and efﬁciency. Therefore, we developed a time-difference adjoint tomography method combining the arrival-time difference and the adjoint state method. The effectiveness was veriﬁed by numerical experiments and a laboratory-scale acoustic experiment. The effectiveness of the proposed method was demonstrated by the experimental results. The adjoint scheme avoids additional ray tracing and improves the efﬁciency of the inversion, which allows the use of ﬁner forward grids in practice. By considering the differential arrivals of receiver pairs, the proposed method is robust in the face of systematic errors and relatively stable against large random noises. Moreover, the velocity contrast obtained by the proposed method is sharper than for ﬁrst-arrival tomography in the areas where the rays are not dense, resulting in a clearer indication of the anomalous areas in the tomographic image.


Introduction
Subsurface space development is accelerating with the progress of engineering technology and the expansion of resource demand [1,2]. However, anomalous geological structures (such as caves, ground fissures, and faults) may affect the stability of the rock engineering and induce various engineering geological disasters [3][4][5][6][7]. For example, a cavity area hidden in the formation may contain toxic gases and stagnant water [8]. During engineering exploration and surrounding rock excavation, potential hazard sources are likely to cause great damage. Therefore, accurately detecting anomalous areas in rocks is of great significance to the safety of construction facilities and personnel in subsurface engineering [9][10][11].
Geophysical methods play an increasingly important role in geological surveys since they can provide more detailed and valuable structural information [12]. Seismic tomography is a typical representative geophysical method and, to date, it has been successfully applied at various engineering scales (e.g., [13][14][15][16][17][18][19]). In engineering geophysics, artificial seismic waves are generally actively excited at the surface by using blasting, hammering and vibration. Through transmission, reflection, or refraction tomography, the geological structure of the survey area, including the anomalies, can be easily obtained [12]. However,

Methodology
We briefly introduce the workflow of TDAT and clarify some basic tomography principles and methods. For simplicity, all descriptions are based on two-dimensional situations. There is no difficulty in extending this workflow to three-dimensional situations. The schematic workflow is shown in Figure 1.

Model Discretization
Under the condition of high-frequency approximation, the travel time of an elastic wave in an isotropic medium approximately satisfies the eikonal equation with the boundary condition t(x) = 0, x = x s . In the continuous case, solving the eikonal equation is very complex. Therefore, the discretization of the inversion area is necessary. For simplicity, we discretize the model into square meshes. The grid length depends on the size of the inversion area, resolution requirement, and observations. Generally, there are no specific criteria regarding the grid length. The recommended grid length is between 1% and 2% of the maximum scale of the computational area. Since the finite difference method is adopted, the wave field and structural characteristics are described by assigning attribute values to the model. In the process of tomography, the model attributes are always defined on the grid points rather than inside the grids. Therefore, forward and inverse processes in tomography are carried out for grid points.

Model Discretization
Under the condition of high-frequency approximation, the travel time of an elastic wave in an isotropic medium approximately satisfies the eikonal equation with the boundary condition ( ) = 0, = .
In the continuous case, solving the eikonal equation is very complex. Therefore, the discretization of the inversion area is necessary. For simplicity, we discretize the model into square meshes. The grid length depends on the size of the inversion area, resolution requirement, and observations. Generally, there are no specific criteria regarding the grid length. The recommended grid length is between 1% and 2% of the maximum scale of the computational area. Since the finite difference method is adopted, the wave field and structural characteristics are described by assigning attribute values to the model. In the process of tomography, the model attributes are always defined on the grid points rather than inside the grids. Therefore, forward and inverse processes in tomography are carried out for grid points.

Forward Calculation
Before introducing the specific inversion process, we must introduce a method to solve Equation (1). Since we chose the adjoint state method for inversion, it is not neces-

Forward Calculation
Before introducing the specific inversion process, we must introduce a method to solve Equation (1). Since we chose the adjoint state method for inversion, it is not necessary to incur extra cost to trace the ray paths. For this reason, the finite difference method is an attractive option. In this paper, the fast sweeping method (FSM) [36] is used to calculate the travel-time field. The partial differential operators are replaced by the Godunov upwind difference operators [37] in FSM to establish a discrete eikonal equation: where t x min = min t i−1,j , t i+1,j and t y min = min t i,j−1 , t i,j+1 . The unique solution of the equation is determined by repeated Gauss-Seidel iterations. FSM consists of three steps: 1 initializing the model with t s = 0, and assigning large positive values at all other grid points; 2 sweeping the domain via Gauss-Seidel iterations and selecting the smaller of the new solution and the original solution; 3 repeating step 2 until ||t k+1 − t k ||≤ ε . In the two-dimensional case, FSM sweeps the whole domain repeatedly with four alternative orderings, and in the three-dimensional case, the necessary orderings are increased to eight (Figure 2a,b). In practice, we adopted an improved FSM scheme with mixed grids to reduce the spread of error from the source grid point. Specifically, the grids near the source were refined to 1/5 of the other grids, and the refined grids were solved by linear interpolation (Figure 2c). where min = min( −1, , +1, ) and min = min( , −1 , , +1 ). The unique solution of the equation is determined by repeated Gauss-Seidel iterations. FSM consists of three steps: ① initializing the model with = 0, and assigning large positive values at all other grid points; ② sweeping the domain via Gauss-Seidel iterations and selecting the smaller of the new solution and the original solution; ③ repeating step ② until ‖ +1 − ‖ ≤ .
In the two-dimensional case, FSM sweeps the whole domain repeatedly with four alternative orderings, and in the three-dimensional case, the necessary orderings are increased to eight (Figure 2a,b). In practice, we adopted an improved FSM scheme with mixed grids to reduce the spread of error from the source grid point. Specifically, the grids near the source were refined to 1/5 of the other grids, and the refined grids were solved by linear interpolation (Figure 2c).

Misfit Function
The difference between the actual model and the computational model is characterized by a misfit function. The most classical misfit function for FAAT is defined as

Misfit Function
The difference between the actual model and the computational model is characterized by a misfit function. The most classical misfit function for FAAT is defined as where M and R represent the number of sources and receivers, respectively, and t(c, r) and t * (r) represent the values of the calculated first arrivals and observed first arrivals, respectively. TDAT adopts a misfit function in a completely different form from Equation (3), as follows: where d ij = t c, r j − t(c, r i ) − t * r j − t * (r i ) . For FAAT, the seismic origin time and the trigger time at the receiver must be known exactly, while for TDAT, the absolute arrival time or the relative arrival time is available.

Model Update
To minimize J(c), we iteratively invert c by the gradient method. The iteration of the inversion starts from a given initial velocity model c 0 , and the velocity model is updated by continuously applying perturbations until the convergence condition is satisfied. The iterative process is formulated as where,c k represents the exploration direction of the model in a single iteration, and α k is the iteration step size determined by inaccurate linear search based on the Wolfe-Powell conditions. Iterative Equation (5) is realized by the quasi-Newton method, which is used to estimate the feasible perturbation directionc. In the Newton method, the first derivative and the second derivative are considered to calculate the descent direction and to realize the second-order convergence. In the quasi-Newton method, the positive definite matrix is obtained by iterative calculation to approximately replace the Hessian matrix, which not only has superlinear convergence speed but also saves required memory effectively. The estimation ofc synthetically depends on the gradient of the current model, and the gradient and descent direction of the historical iterative model.

Adjoint Variable and Gradient
The nonlinear relationship between the misfit function value J and the velocity model c is implicit, and the conventional method for the gradient requires the computation of Fréchet derivatives. We chose the adjoint state method to calculate the gradient ∇J, which can effectively avoid the computation of Fréchet derivatives [38]. The use of the adjoint scheme can avoid the problem of multi-paths in extreme cases [39] and effectively reduce the inversion cost under dense grid conditions.
The computation of the adjoint field in FAAT has been deduced (see [35,40,41]). In this paper, we derive the adjoint field calculation formulas for TDAT, and the detailed derivation is given in Appendix A. Within the subsurface, the adjoint state variable λ is the solution of ∇·λ∇t = 0, On the surface, we choose the adjoint state variable λ, which is the solution of where n represents the unit vector normal to the surface. Compared with FAAT, solving for the adjoint variable of TDAT only differs in the boundary initialization, and therefore it does not increase any computational complexity. To solve Equations (6) and (7), we adopted a finite difference method based on Gauss-Seidel iterations [35]. Then, λ is normalized by c 3 to obtain the gradient of the misfit function for a single source, and the final gradient used to update the velocity model is the summation of the contributions of all sources [41].
The role of the adjoint system for gradient calculation can be illustrated through a numerical example. Figure 3a shows an initial velocity model with a constant vertical gradient, where the true velocity model contains a positive velocity anomaly (Figure 3b). The synthetic velocity model is obtained by superposition of the initial velocity model and the velocity anomaly. The source is located at the surface, i.e., (0, 0) km. We obtain the gradient for source-receiver pairs by solving the adjoint state variable, as shown in Figure 3c-f. Here, the results clearly show the function of the adjoint system: the adjoint state variable is initialized according to the residuals at the surface, and the initial values then propagate from the receiver to the source in the reverse direction of the ray. In contrast to FAAT, the adjoint source of TDAT contains relative information about the receiver pairs, so the gradient in a local area may be inaccurate. However, as the number of sources and receivers increases, this error can be eliminated by superposition. In the area not covered by the rays, the gradient is 0, which means that the model does not update at these locations. Increasing the number of receivers can improve the ray coverage, and accordingly the gradient calculated by the adjoint state method also covers a larger area.

=1
The role of the adjoint system for gradient calculation can be illustrated through a numerical example. Figure 3a shows an initial velocity model with a constant vertical gradient, where the true velocity model contains a positive velocity anomaly (Figure 3b). The synthetic velocity model is obtained by superposition of the initial velocity model and the velocity anomaly. The source is located at the surface, i.e., (0, 0) km. We obtain the gradient for source-receiver pairs by solving the adjoint state variable, as shown in Figure 3c-f. Here, the results clearly show the function of the adjoint system: the adjoint state variable is initialized according to the residuals at the surface, and the initial values then propagate from the receiver to the source in the reverse direction of the ray. In contrast to FAAT, the adjoint source of TDAT contains relative information about the receiver pairs, so the gradient in a local area may be inaccurate. However, as the number of sources and receivers increases, this error can be eliminated by superposition. In the area not covered by the rays, the gradient is 0, which means that the model does not update at these locations. Increasing the number of receivers can improve the ray coverage, and accordingly the gradient calculated by the adjoint state method also covers a larger area.

Regularization
It can be seen from Figure 3 that the gradient clearly has a sharp spike near the source location. The point-source travel-time field at the source location is non-differentiable due to the source singularity [42], which further leads to a significant increase in the gradient near the source. In addition, the rays highly overlap at the grid points near the receiver locations. The anomalies cause significant interference in the iterations, resulting in the local velocity values being too high or too low. Therefore, proper regularization is necessary for the stability of the inversion.
We used regularization twice in each iteration to help the algorithm find a more appropriate update direction and to obtain a more robust result. The first regularization was used to weaken the anomalies near the sources and receivers. We convoluted the gradient with a larger Gaussian kernel to reduce the weight of the gradient and thus to directly constrain the update near the surface, as shown in Figure 4. The first regularization hardly changed the gradient values in the middle of the computational area. Accordingly, it greatly reduced the gradient near the sources and receivers. This regularization is effective in first-arrival tomography but is only applicable when the sources and receivers are located at the boundary of the computational area. The second regularization is used to smooth the velocity model. Since the velocity model is discrete, proper smoothing is essential to reduce the discontinuity. In this regularization, a small low-pass filter is applied to smooth the velocity model. The filtering range is recommended to cover the whole computational area. In addition, it is necessary to set the boundary of the velocity, including the upper limit and lower limit. Since we used a low-pass filter in each iteration, we only used the cubic interpolation once to optimize the meshes before the final output of the tomography result.
location. The point-source travel-time field at the source location is non-differentiable due to the source singularity [42], which further leads to a significant increase in the gradient near the source. In addition, the rays highly overlap at the grid points near the receiver locations. The anomalies cause significant interference in the iterations, resulting in the local velocity values being too high or too low. Therefore, proper regularization is necessary for the stability of the inversion.
We used regularization twice in each iteration to help the algorithm find a more appropriate update direction and to obtain a more robust result. The first regularization was used to weaken the anomalies near the sources and receivers. We convoluted the gradient with a larger Gaussian kernel to reduce the weight of the gradient and thus to directly constrain the update near the surface, as shown in Figure 4. The first regularization hardly changed the gradient values in the middle of the computational area. Accordingly, it greatly reduced the gradient near the sources and receivers. This regularization is effective in first-arrival tomography but is only applicable when the sources and receivers are located at the boundary of the computational area. The second regularization is used to smooth the velocity model. Since the velocity model is discrete, proper smoothing is essential to reduce the discontinuity. In this regularization, a small low-pass filter is applied to smooth the velocity model. The filtering range is recommended to cover the whole computational area. In addition, it is necessary to set the boundary of the velocity, including the upper limit and lower limit. Since we used a low-pass filter in each iteration, we only used the cubic interpolation once to optimize the meshes before the final output of the tomography result.

Results and Discussion
To evaluate the performance of TDAT for detection of anomalous areas in rocks, numerical experiments and a laboratory-scale acoustic monitoring experiment were conducted.

Numerical Experiments
In the numerical experiments, the velocity anomaly and the corresponding scale that the proposed method can identify in an ideal environment were evaluated. In a constant velocity model with a size of 100 m × 100 m, the background velocity is 5000 m/s. As shown in Figure 5, the first numerical experiment considered the influence of anomalous areas of different sizes. The white dotted line represents the circular

Results and Discussion
To evaluate the performance of TDAT for detection of anomalous areas in rocks, numerical experiments and a laboratory-scale acoustic monitoring experiment were conducted.

Numerical Experiments
In the numerical experiments, the velocity anomaly and the corresponding scale that the proposed method can identify in an ideal environment were evaluated. In a constant velocity model with a size of 100 m × 100 m, the background velocity is 5000 m/s. As shown in Figure 5, the first numerical experiment considered the influence of anomalous areas of different sizes. The white dotted line represents the circular anomalous area in the model, with a radius of 2 m, 4 m, 8 m, and 16 m, respectively. The velocity in the anomalous areas was 6000 m/s. It can be seen from the figure that TDAT accurately captured the high-velocity anomaly in the center of the computational area. In the ideal environment of dense acquisition, TDAT can identify an anomaly with a radius of 2% of the computational area. With the increase in the size of the anomalous area, TDAT gradually reconstructs the wave velocity value.
The second numerical experiment considered the velocity difference between the background and the anomalous area. The radius of the abnormal area was 10 m, and the corresponding wave velocity difference was between −2000 m/s and +2000 m/s. As shown in Figures 6 and 7, TDAT can identify a velocity anomaly exceeding 250 m/s in an ideal environment. Further, it can be observed that TDAT is more sensitive to high-velocity anomalies. This is related to the propagation characteristics of seismic waves. In the inversion process, the rays tend to be close to the high-velocity area. Therefore, the update in the low-velocity area will decrease.
anomalous area in the model, with a radius of 2 m, 4 m, 8 m, and 16 m, respectively. The velocity in the anomalous areas was 6000 m/s. It can be seen from the figure that TDAT accurately captured the high-velocity anomaly in the center of the computational area. In the ideal environment of dense acquisition, TDAT can identify an anomaly with a radius of 2% of the computational area. With the increase in the size of the anomalous area, TDAT gradually reconstructs the wave velocity value. The second numerical experiment considered the velocity difference between the background and the anomalous area. The radius of the abnormal area was 10 m, and the corresponding wave velocity difference was between −2000 m/s and +2000 m/s. As shown in Figures 6 and 7, TDAT can identify a velocity anomaly exceeding 250 m/s in an ideal environment. Further, it can be observed that TDAT is more sensitive to high-velocity anomalies. This is related to the propagation characteristics of seismic waves. In the inversion process, the rays tend to be close to the high-velocity area. Therefore, the update in the low-velocity area will decrease.
The numerical experiments show that TDAT is very effective in an ideal environment. The adjoint system correctly captures the information from the receivers to the sources without ray tracing. The practical application of TDAT may be affected by many factors such as the rock environment and the monitoring environment.     The numerical experiments show that TDAT is very effective in an ideal environment. The adjoint system correctly captures the information from the receivers to the sources without ray tracing. The practical application of TDAT may be affected by many factors such as the rock environment and the monitoring environment.

Experimental Overview
A laboratory-scale homogeneous granite model with five cylindrical holes was used in the experiment (Figure 8a,b). The size of the three-dimensional model was 50 cm × 20 cm × 16 cm. The five anomalous holes were symmetrically located in the model, and their diameter was about 5 cm. Three of the five holes did not completely penetrate the model from top to bottom, while the other two penetrated completely. We filled the three non-penetrating cylindrical holes (A, B, and C) with water, as shown in Figure 8c. From one-dimensional measurement, the propagation velocity of an ultrasonic wave in the granite medium was about 4500 ± 500 m/s, that in the water medium was about 1500 m/s, and that in air was about 340 m/s. Twenty-four acoustic sensors were used in the experiment to monitor the acoustic signals. The sensors have the function of transmitting high-frequency and low-frequency pulse signals, and therefore they act as both sources and receivers. The sensors were symmetrically arranged on the front and rear sides of the model (Figure 8c). They were located in the center of the grid points, and the distance between two adjacent sensors was 4 cm.
that in air was about 340 m/s. Twenty-four acoustic sensors were used in the experiment to monitor the acoustic signals. The sensors have the function of transmitting high-frequency and low-frequency pulse signals, and therefore they act as both sources and receivers. The sensors were symmetrically arranged on the front and rear sides of the model (Figure 8c). They were located in the center of the grid points, and the distance between two adjacent sensors was 4 cm. Half of the 24 sensors were used as ultrasonic sources and the other half as receivers. Therefore, the 24 sensors formed 144 source-receiver pairs. After 180-degree rotation, the number of source-receiver pairs increased to 288. Due to the inhomogeneity and anisotropy of the model and the influence of the experimental environment, the data for the latter 144 were not completely consistent with the data for the former 144. Therefore, we considered the contribution of all 288 source-receiver pairs during inversion. Figure 8d shows the 144 rays in the experiment. It can be seen that the rays cover the middle, the upper side, and the lower side of the inversion area well, while only a few rays pass through the left and right sides (the two sides without sensors).
After the acoustic signal acquisition system was turned on, the 24 sensors transmitted 10 MHz pulses in turn. For each sensor, the pulse was transmitted four times in a row. The acquisition system recorded all acoustic signals received by the 24 sensors during the experiment. We turned on the amplitude limit of the collected signal to filter most of the reflected waves and background noise. To ensure the reliability of the data, three repeated Half of the 24 sensors were used as ultrasonic sources and the other half as receivers. Therefore, the 24 sensors formed 144 source-receiver pairs. After 180-degree rotation, the number of source-receiver pairs increased to 288. Due to the inhomogeneity and anisotropy of the model and the influence of the experimental environment, the data for the latter 144 were not completely consistent with the data for the former 144. Therefore, we considered the contribution of all 288 source-receiver pairs during inversion. Figure 8d shows the 144 rays in the experiment. It can be seen that the rays cover the middle, the upper side, and the lower side of the inversion area well, while only a few rays pass through the left and right sides (the two sides without sensors).
After the acoustic signal acquisition system was turned on, the 24 sensors transmitted 10 MHz pulses in turn. For each sensor, the pulse was transmitted four times in a row. The acquisition system recorded all acoustic signals received by the 24 sensors during the experiment. We turned on the amplitude limit of the collected signal to filter most of the reflected waves and background noise. To ensure the reliability of the data, three repeated experiments were carried out. A total of 7084 events were recorded. These events still contained some noise and reflected signals with high amplitude. After screening, only 6912 first-arrival events were retained as the inversion basis.

Resolution
Checkerboard resolution tests were conducted to determine the resolution of TDAT and to help to select the best grid length matching the observations. The 44 cm × 20 cm area that the first arrivals could cover was selected as the inversion area. The initial model adopted the constant velocity model, and the background velocity was set to 4500 m/s. The checkerboard velocity model c(x) is constructed by applying a velocity perturbation [39] to the initial velocity model c 0 (x): where γ is the maximum relative velocity perturbation and l x , l y represent the lengths of the anomalous areas in the horizontal and vertical axis directions, respectively.
A total of three resolution tests were carried out, considering grid lengths of 1 cm, 0.5 cm, and 0.25 cm (Table 1). Thirty-two square regions with velocity perturbations were evenly distributed in the rectangular inversion area according to Equation (9), as shown in Figure 9. In the tests, we let γ = 8% and l x = l y = 5 cm. The initial environment was completely consistent with the experimental setup, and the parameters required for inversion were exactly the same in each test.   The imaging results after 30 iterations are shown in Figure 10. It can be seen that TDAT and FAAT actually have a similar tomographic resolution. The results also clearly show that the resolution is related to the grid length. The denser grids provide more detail on a small scale, while the coarse grids show more uniform changes in the inversion. In test 1, the adjoint tomography almost successfully inverted the size and position of most of the perturbation regions, while in test 2 and test 3, neither TDAT nor FAAT met the resolution requirements.
As in any other tomography method, the observed data form one of the key factors affecting imaging resolution. The best tomography result requires the number of grids to match the amount of data. Therefore, we adopted an inversion grid length of 1 cm for imaging. In addition, we used a grid length of 0.2 cm in the forward modeling, to obtain a more accurate travel time. The transformation between sparse grids and dense grids was realized by interpolation.

Imaging
TDAT and FAAT were used to reconstruct the velocity structure ( Figure 11). The visualization images show that the tomography results of both TDAT and FAAT accurately reflect the velocity structure characteristics of the hole-containing granite medium model. In the velocity maps, the positions of the five low-velocity areas are basically consistent with the reality. It is very clear that the TDAT image shows a sharper velocity contrast, especially near the sources and receivers. This characteristic of TDAT is in accordance with classical double-difference ray tomography [27]. For FAAT, the background velocity appeared to rise more strongly and deviated from the normal value, and intense diffusion occurred in the low-velocity areas. Perhaps, therefore, more powerful smoothing was required.
From the convergence curves (Figure 12), the value of the misfit function in the two cases decreased by 96.14% and 97.67%, respectively. The quasi-Newton method showed excellent convergence performance. After about 30 iterations, the tomographic result was basically stable, and the subsequent iterations had only a tiny impact. The operation speed of TDAT was 27.7% slower, and the additional computational cost arises from the more complex misfit function.
From the convergence curves (Figure 12), the value of the misfit function in the two cases decreased by 96.14% and 97.67%, respectively. The quasi-Newton method showed excellent convergence performance. After about 30 iterations, the tomographic result was basically stable, and the subsequent iterations had only a tiny impact. The operation speed of TDAT was 27.7% slower, and the additional computational cost arises from the more complex misfit function. Overall, the proposed TDAT method was very effective in the experiment, accurately identifying the five anomalous areas in the granite medium. More importantly, the use of differential arrivals of receivers provided a clearer velocity contrast, which can help to locate the anomalies in rocks better.

Robustness to Observations
In engineering geophysics, the seismic tomography method must be robust enough to deal with anisotropy in rocks, background noises, measurement deviations, source location errors, and system uncertainty. Here, robustness tests were conducted to investigate the influence of errors on the results for TDAT and FAAT.
The physical velocity model we used cannot be accurately measured in practice, therefore in all the robustness tests, the raw observations were assumed to be accurate and error-free. A constant offset was added to the trigger time of sources, to simulate the system errors, and Gaussian random noises were added to the first arrivals of receivers, to simulate the influence of uncertainty and anisotropy. The settings of the errors for the eight groups of robustness tests are listed in Table 2. Based on the synthetic observations, the velocity structure of the granite medium model was reconstructed by TDAT and FAAT (Figures 13 and 14). In the FAAT inversions, a stronger regularization was used to Overall, the proposed TDAT method was very effective in the experiment, accurately identifying the five anomalous areas in the granite medium. More importantly, the use of differential arrivals of receivers provided a clearer velocity contrast, which can help to locate the anomalies in rocks better.

Robustness to Observations
In engineering geophysics, the seismic tomography method must be robust enough to deal with anisotropy in rocks, background noises, measurement deviations, source location errors, and system uncertainty. Here, robustness tests were conducted to investigate the influence of errors on the results for TDAT and FAAT. The physical velocity model we used cannot be accurately measured in practice, therefore in all the robustness tests, the raw observations were assumed to be accurate and error-free. A constant offset was added to the trigger time of sources, to simulate the system errors, and Gaussian random noises were added to the first arrivals of receivers, to simulate the influence of uncertainty and anisotropy. The settings of the errors for the eight groups of robustness tests are listed in Table 2. Based on the synthetic observations, the velocity structure of the granite medium model was reconstructed by TDAT and FAAT (Figures 13 and 14). In the FAAT inversions, a stronger regularization was used to maintain the background velocity within a reasonable range. erroneously identified an extra anomalous low-velocity area (Figure 14a). The results of test 1 suggest that TDAT is far more robust than FAAT in the face of systematic errors. In test 2, a relatively high level of random noise was added, and the imaging result for TDAT showed a significant perturbation. However, the identified anomalous areas were still very clear. The result for FAAT also showed a significant change, and the velocity contrast was relatively less obvious. In test 3 to test 8, two errors were added to the observations, and the constant system errors were higher than the random noises. For TDAT, the effect of the errors mainly originated from the random noises. The irregular  contrast between the background velocity and the velocity of the holes was still very sharp, even in areas where the rays were not dense. For FAAT, the influence of the error spread to the entire inversion area. On the one hand, the velocity contrast was not obvious, and on the other hand, the irregular variation was more significant. In particular, large increases or decreases in background velocity significantly reduced FAAT's robustness and ability to identify anomalous areas. To further compare the performance of the two methods, we considered the parameter IR (identification ratio), defined as For test 1, since only relative information for the arrival times of receivers was considered, the constant offset applied to the source trigger time did not affect the tomography result for TDAT (Figure 13a). However, faced with the interference of systematic errors, FAAT could not correctly identify the anomalous areas in the structure and it erroneously identified an extra anomalous low-velocity area (Figure 14a). The results of test 1 suggest that TDAT is far more robust than FAAT in the face of systematic errors.
In test 2, a relatively high level of random noise was added, and the imaging result for TDAT showed a significant perturbation. However, the identified anomalous areas were still very clear. The result for FAAT also showed a significant change, and the velocity contrast was relatively less obvious. In test 3 to test 8, two errors were added to the observations, and the constant system errors were higher than the random noises. For TDAT, the effect of the errors mainly originated from the random noises. The irregular variation of the images became gradually obvious as the standard deviation of the Gaussian noise increased. Nevertheless, the identification of the five anomalous areas was still quite robust in the imaging results for TDAT. In particular, as the errors increased, the contrast between the background velocity and the velocity of the holes was still very sharp, even in areas where the rays were not dense. For FAAT, the influence of the error spread to the entire inversion area. On the one hand, the velocity contrast was not obvious, and on the other hand, the irregular variation was more significant. In particular, large increases or decreases in background velocity significantly reduced FAAT's robustness and ability to identify anomalous areas.
To further compare the performance of the two methods, we considered the parameter IR (identification ratio), defined as where S overlap represents the overlap of identified anomalous areas and real anomalous areas, and S anomaly represents the total area of a real anomaly. A grid point is considered to be anomalous when the difference between the grid-point velocity and the background velocity exceeds 20%. In other words, in this experiment, the grid point is included in the anomalous area only when the velocity at the grid point is lower than 3600 m/s. Figure 15 shows the values of IR for TDAT and FAAT in eight groups of tests. In all groups of robustness tests, the IR value for TDAT was larger than that for FAAT.  One particular point to note is that TDAT and FAAT should have different ass ment criteria, due to the difference in the misfit functions. In this case, the same arr error had different influence for TDAT and FAAT. Due to the relatively short dista between receivers (the distance between adjacent receivers was 4 cm), the relative arri One particular point to note is that TDAT and FAAT should have different assessment criteria, due to the difference in the misfit functions. In this case, the same arrival error had different influence for TDAT and FAAT. Due to the relatively short distance between receivers (the distance between adjacent receivers was 4 cm), the relative arrivals difference for some receiver pairs was very small (even less than 0.5 µs). Under such conditions, even a small noise can cause a large relative error. At this point, in practice, appropriately increasing the distance between the receivers can effectively reduce the impact of random errors on TDAT.
The comprehensive results reveal that the arrival-time difference can effectively cope with the system errors. Moreover, TDAT can maintain stability in the face of random noises, although the relative errors added are quite large. Nevertheless, TDAT can still obtain a sharper velocity contrast in the face of errors. This feature is essential for detecting anomalous areas in rocks, because usually only the locations and the scales of anomalous areas are important.

Conclusions
A time-difference adjoint tomography method based on the finite difference method and a quasi-Newton method was developed to detect anomalous areas in rocks. The adjoint scheme avoids additional ray tracing, which allows a denser forward mesh to be used in practice. Numerical experiments and a laboratory-scale acoustic monitoring experiment were conducted to verify the effectiveness. The resolution, ability to detect anomalous areas, and robustness of TDAT were comprehensively evaluated and compared with FAAT.
The workflow of TDAT was proven to be very effective, according to the results. Despite the influence of anisotropy and the uncertainty that existed objectively, the five anomalous areas in the rock structure were still accurately identified. When using active seismic sources, TDAT and FAAT had similar resolutions. However, the velocity images of TDAT had a sharper contrast, which is of great significance for the identification of anomalous areas in rocks. The use of arrival-time difference can help to accurately confirm the locations and sizes of anomalous areas.
Since only the relative information between the receivers is considered, TDAT is more robust in the face of systematic errors and remains stable in the face of large random noises. Even if the noises are mixed and increased, the velocity contrast in the result for TDAT is still clear. TDAT is far more effective than FFAT, especially in areas where rays are not dense.
Furthermore, this scheme can be extended to passive seismic sources (such as microseismic and acoustic emission sources) since the workflow does not involve any temporal information about the seismic source.
Author Contributions: Conceptualization: F.W. and L.D.; methodology: X.X. and Z.P.; investigation: X.X. and Z.P.; software: F.W. and X.X.; writing-original draft: X.X.; writing-review and editing: F.W. and L.D. All authors have read and agreed to the published version of the manuscript.

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

Appendix A. Adjoint Scheme for Time-Difference Tomography
We derive the adjoint state scheme of the common-source time-difference misfit function in the discrete case. Without loss of generality, we first consider the case of only one source. Then, the common-source misfit function can be defined as where d ij = t c, r j − t(c, r i ) − t * r j − t * (r i ) . According to the Lagrangian formulation, the extended misfit function becomes L(c, t, λ) = 1 2 During the minimization process, t and λ are independent of c, and the gradient of the misfit function ∂J ∂c is obtained by solving Here, we focus on Equation (A5), which is the key to solving the adjoint system. Notice that ∂d ij ∂t(x) =    0, x / ∈ r i , r j , −1, x = r i , 1, x = r j .
(A6) Therefore, we can choose λ satisfying λ(x)∇t (x) ∂∇t(x) ∂t = 0, x ∈ Ω/r, with the boundary condition Approximately, Equation (A8) can be written as where n is the unit vector normal to the surface. We can solve Equations (A7) and (A9) according to a finite difference method based on Gaussian-Seidel sweeping [35], and the only thing to note is the change of initialization conditions. Finally, the gradient of the misfit function ∂J ∂c is obtained by directly solving Equation (A3). For the case of multiple sources, the contribution of all sources should be considered, so Equation (A9) then becomes