Multidirectional Anisotropic Total Variation and Its Application in the Tomography of the Surrounding Rock of Coal Mining Faces

,e computed tomography (CT) reconstruction algorithm is one of the crucial components of the CT system. To date, total variation (TV) has been widely used in CT reconstruction algorithms. Although TV utilizes the a priori information of the longitudinal and lateral gradient sparsity of an image, it introduces some staircase artifacts. To overcome the current limitations of TV and improve imaging quality, we propose a multidirectional anisotropic total variation (MATV) that uses multidirectional gradient information.,e surrounding rock of coal mining faces uses principles of tomography similar to those of medical X-rays. ,e velocity distribution for the surrounding rock can be obtained by the first-arrival traveltime tomography of the transmitted waves in the coal mining face. Combined with the geological data, we can interpret the geological hazards in the coal mining face. To perform traveltime tomography, we first established the objective function of the first-arrival traveltime tomography of the transmitted waves based on theMATV regularization and then used the split Bregmanmethod to solve the objective function.,e simulated data and real data show that the MATV regularization method proposed in this paper can better maintain the boundaries of geological anomalies and reduce the artifacts compared with the isotropic total variation regularization method and the anisotropic total variation regularization method. Furthermore, this approach describes the distribution of geological anomalies more accurately and effectively and improves imaging accuracy.


Introduction
X-ray computed tomography (CT) has been widely used in medical diagnosis [1][2][3][4] since it was invented in 1973. X-ray CT scans automatically explore the internal structure of the human body. One of the key contents of X-ray CT research is the CT reconstruction algorithm [1], which can be divided into an analytical algorithm and an iterative algorithm. e analysis algorithm requires high data completeness. However, in most cases, the data are incomplete, and the image reconstructed by the analysis algorithm will produce more artifacts [3]. Direct iterative algorithms are currently widely used in CT, such as the simultaneous algebraic reconstruction technique (SART) [5], algebraic reconstruction technique (ART) [6], and adaptive algebraic reconstruction technique (AART) [7]. ere is also a type of iterative algorithm that uses the image gradient sparse a priori information, which is known as total variation (TV) regularization. Sidky and Pan [2] introduced TV into medical X-ray CT for the first time. According to the existing research [3,4], the accuracy of this kind of algorithm is higher than the accuracy of a direct iterative method.
TV was proposed by Rudin et al. [8] while studying image denoising. Subsequently, TV has been widely used in image processing [9][10][11], image reconstruction [12,13], and geophysical inversion [14][15][16][17]. However, the results of Chan et al. [18] show that TV produces staircase artifacts. Researchers around the globe have made numerous attempts to suppress staircase artifacts and improve imaging accuracy. ese improvements can be divided into the following categories: (1) generalize the difference of operation in TV to a higher-order difference [18][19][20] or fractional-order difference [21] such as high-order TV [18][19][20][21], total generalized variation [22], and fractional-order TV [23,24]; (2) extend the L1 norm in TV to the Lp quasinorm (0<p < 1) that can obtain sparser solution than the L1 norm, such as TpV [11,17,25,27]; (3) generalize the gradient in TV to the overlapping combination gradient information, such as overlapping group sparsity total variation [28][29][30][31][32]; (4) use Tikhonov regularization and TV regularization to solve the objective function in an alternative way, such as modified TV [33][34][35]; and (5) generalize the two-directional gradient to multidirectional gradient in TV, such as multidirectional TV [3,4,36]. TV can be divided into isotropic total variation (ITV) and anisotropic total variation (ATV) [10,37]. Wang et al. [38] showed that ATV has higher imaging accuracy than ITV. Geological anomalies such as faults, collapse columns, and goafs in coal mining faces can affect production efficiency, which is potentially challenging for safe production [39][40][41]. With the rapid development of comprehensive mechanized mining technology, it is crucial to identify geological anomalies in coal mining faces. Because of the low precision of 3D seismic exploration on the ground, it cannot meet the requirements for safe and efficient coal mine production. Seismic transmission exploration in the coal mining face has the potential to find the geological anomalies in the coal mining face more precisely. e refracted waves that propagate along the surrounding rock will be generated when the seismic waves are excited in the coal mining face. e velocity of the refracted waves that reach the receivers first is the fastest. erefore, the first-arrival traveltime of the transmitted waves in the coal mining face is the propagation time of the refracted waves in the surrounding rock. When the velocity of the surrounding rock changes, it will affect the first-arrival traveltime. rough traveltime tomography, the distribution of the surrounding rock velocity can be obtained. Combined with the geological data from the coal mining face for comprehensive interpretation, we can obtain the geological anomaly distribution map of the coal mining face. To take advantage of the multidirectional gradient information and improve imaging accuracy, we propose a multidirectional anisotropic total variation (MATV) and apply it to the tomography of the surrounding rock of a coal mining face. e results of the two examples show the effectiveness of the proposed method.

Multidirectional Anisotropic Total Variation
As shown in Figure 1(a), the traditional ATV only uses the gradient information of a pixel and its two adjacent pixels, and its expression can be written as follows [10,37]: where P is the pixel matrix, P m,j is the element of the matrix P, and ∇ x P and ∇ y P represent the first-order difference of the pixels in the x direction and y direction, respectively. e matrix P is transformed into a column vector in calculation.
Assuming that the size of the matrix P is r × s, the size of the column vector is (r × s) × 1. ∇ x and ∇ y are different matrices.
, and . Previous studies [3,4] reported that the ITV uses the multidirectional gradient information of pixels and showed that multidirectional ITV could effectively improve imaging accuracy. Inspired by their research, we propose a novel ATV that uses the 8-direction gradient information of a pixel, that is, MATV. For a two-dimensional image, as shown in Figure 1(b), the two-dimensional MATV uses the gradient information of a pixel and its 8 adjacent pixels, and its specific expression is defined as where ∇ 1 P, ∇ 2 P, ∇ 3 P, ∇ 4 P, ∇ 5 P, ∇ 6 P, ∇ 7 P, and ∇ 8 P represent the first-order difference of the pixel in 8 adjacent directions.

Transmitted Waves of a Coal Mining Face.
A typical three-layer symmetric model of rock-coal-rock is shown in Figure 2(a). e velocity of the coal seam is generally lower than the velocity of the surrounding rock. e observation system of seismic transmission exploration in a coal mining face is shown in Figure 2(b). Shot points are placed in one roadway of the coal mining face, and geophones are placed in the other roadway to receive the transmitted waves that pass through the coal mining face. rough transmitted wave analysis, we can obtain the geological anomalies of the coal mining face.
Exciting seismic waves in the coal seam will produce refracted waves that propagate along the surrounding rock, direct waves, and channel waves that propagate in the coal seam. As shown in Figure 2(c), when the seismic waves' incident angle is equal to the critical angle θ, refracted waves that propagate along the surrounding rock will be generated, and its propagation velocity is similar to the velocity of the surrounding rock. erefore, we can obtain the velocity of the surrounding rock through the traveltime tomography of refracted waves. Generally, the P-wave velocity of refracted waves is larger and reaches the geophone earlier. erefore, we can pick up the first-arrival traveltime of the transmitted waves for tomography to acquire the P-wave velocity distribution map of the surrounding rock. e refraction angle θ is 30-45 degrees, and the coal thickness d is typically 2-6 m. It can be estimated that r in Figure 2 en, we can correct the first-arrival traveltime of the transmitted waves.

Establishment of the Objective Function.
In the two-dimensional plane of x and y, if the path of the transmitted wave is L and the velocity of the medium is V(x, y), then the traveltime T of the transmitted waves can be expressed as When the medium velocity changes, the first-arrival traveltime is a nonlinear function of the medium velocity, which can be expressed as follows: where S is the slowness of the model and its value is the reciprocal of the velocity. By using the Taylor series, the first-order derivative of L(S) is written as where S k is the slowness of the model at the k-th iteration.
where A k is the length matrix of the ray on the model grid at the k-th iteration. e backward process of linear traveltime interpolation (LTI) [42] is used to solve A k in this paper. L(S k ) is the traveltime at the k-th iteration. In this paper, the fast sweeping method (FSM) [43] is used to calculate its value. In general, the geological anomalies in the surrounding rock of the coal mining face are relatively isolated, and its velocity distribution is local to the background. We assume that the velocity of the surrounding rock of the coal mining face has sparsity in the gradient domain. In this context, we 4 Shock and Vibration can use the a priori information that velocity has sparsity in the gradient domain to establish an objective function. In the two-dimensional plane, we can establish the constrained optimization problem of the first-arrival traveltime tomography as follows: where d is the first-arrival traveltime, L(S) is the theoretical traveltime calculated by the model's slowness S, λ is the regularization factor, ‖ ‖ 2 represents the errors between d and L(S), and arg min S indicates the value of S when the objective function is minimized.

Solution Method.
e split Bregman algorithm [44] is used to solve equation (8). First, auxiliary variables are introduced, namely,

Shock and Vibration
According to Goldstein and Osher [44], the regularization factor α is set to α � 2λ.
By introducing the intermediate variables of q 1 , q 2 , q 3 , q 4 , q 5 , q 6 , q 7 , and q 8 , equation (9) is transformed into arg min arg min We set the first derivative of the first subproblem (11) to zero [45] and then obtain where A � Equation (13) can be efficiently solved by using the conjugate gradient method. e second subproblem can be solved with a shrinkage operator [46]: where shrink is defined as shrink(θ, α) � sign(θ) max (|θ| − α, 0). e pseudocode of using the split Bregman algorithm to solve equation (8) is shown in Algorithm 1.

Simulated Data Example
When there are fractures in the rock, its velocity will generally decrease. Geological anomalies such as faults, collapse columns, and goafs will lead to the development of rock fractures to varying degrees, which will reduce their velocity. In this section, we design the model with two local lowvelocity blocks on a high-velocity background to test the method proposed in this paper and we compare the results of the MATV regularization with the results of the ITV regularization and ATV regularization. Figure 3(a), the size of the model is 100 m × 100 m, and it contains two local lowvelocity blocks (geological anomalies). e velocities of the background and the two low-velocity blocks are 3500 m/s and 2500 m/s, respectively. As shown in Figure 3(b), by using a 2 m × 2 m grid to divide the model, a total of 2500 grid units can be obtained. In Figure 3(b), receivers are set on the left and top edges of the model and shots are set on the right and bottom edges of the model to form an observation system of seismic transmission exploration. In this observation system, the shot spacing is 2 m, and there are 100 shots in total, the receiver spacing is 2 m, and 100 receivers per shot.

Curved Ray-Tracing. As shown in
In this paper, the FSM is used for traveltime calculation and the backward process of LTI is used for ray-tracing. To verify the accuracy of the traveltime calculation and raytracing, we design the model with the size of 500 m × 500 m and a velocity of 2000 m/s. Figure 4(a) is a comparison chart of the contours of the traveltime calculated by the FSM and theoretically. It can be seen that the two are coincident, which indicates that the error is small. Figure 4(b) shows the error obtained by subtracting the traveltime calculated by the FSM from the theoretical traveltime. It can be seen from Figure 4(b) that the error is on the order of 10 − 5 s, which also shows that the traveltime calculated by the FSM has high accuracy. Figure 5 is a comparison diagram of the ray paths traced by the backward process of LTI and theoretically. Except near the source, the other places are basically coincident, which shows the accuracy of the ray paths.  e traveltime contours curve at the low-velocity blocks, and the traveltime increases due to the existence of the low-velocity blocks. Figure 6(b) shows the traveltime of receivers R1-R51 in Figure 6(a) and more intuitively reflects the changes in the traveltime with each receiver. Figure 6(c) is the ray path diagram of the 25th shot. e ray paths bypass the lowvelocity blocks to find the minimum traveltime path and present a curved shape.

Comparison of Several Tomography Methods.
e firstarrival traveltime data are calculated by the FSM with the model in Figure 3(a). e initial value is the result of the back projection technique (BPT) (Figure 7(a)). en, we compare the results of the ITV regularization (Figure 7(b)), the ATV regularization (Figure 7(c)), and the proposed MATV regularization (Figure 7(d)).
In the TV-based method, it is important to select the appropriate regularization factors. e regularization factor λ in this paper plays a role in balancing the data fitting term and the regularization term. On the one hand, λ that is too small gives a solution that mainly fits the data and has no better representation of the regularization term. On the other hand, a large value of λ results in the poor fitting of data and generates an incorrect solution. In this paper, we choose the regularization factor from several experiments. We notice that the regularization factor λ with 0.1 can obtain a better inversion result and can ensure that the data fitting error is small.
According to the traveltime tomography results of several methods in Figure 7(b)-7(d), all methods can predict the geological anomalies to a certain extent. e results of the BPT (Figure 7(a)) are the worst, which shows the background velocity and can hardly identify the position of low-velocity blocks. e results of the MATV regularization are in good agreement with the low-velocity blocks, followed by the ATV regularization and ITV regularization. In addition, there are obvious "ripple" artifacts in the reconstructed images of the ITV regularization and ATV regularization, which are effectively reduced by the MATV regularization. us, compared with the ITV regularization and ATV regularization, the proposed MATV regularization method has high imaging accuracy and can improve image quality.   Shock and Vibration To intuitively compare the several traveltime tomography results in Figure 7, we extract the results at x � 28 m and x � 68 m for analysis, as shown in Figure 8. e results of the ITV regularization and the ATV regularization deviate from the model in the background and low-velocity regions, while the results of the MATV regularization are relatively close to the model. e tomography errors are obtained by subtracting the model shown in Figure 3(a) from the traveltime tomography results of each method in Figure 7(b)-7(d). Figure 9 shows the histogram curves of the errors for the ITV regularization, ATV regularization, and MATV regularization in which the horizontal axis represents the velocity and the vertical axis represents the number of points. As shown in Figure 9, when the velocity is near 0, the value is significantly higher for the MATV regularization method than for the ITV and the ATV regularization methods. at is to say, compared with the other two methods, the tomography results of the MATV regularization method have more points equal to the model. is indicates that the errors of the MATV regularization method are smaller than the errors of the other two methods at most points.

Quantitative Evaluation.
To evaluate the traveltime tomography results quantitatively, the universal quality index (UQI) [47] and root mean square error (RMSE) are introduced. e RMSE is a common evaluation index of imaging results. When the RMSE is smaller, the imaging quality is better; otherwise, the imaging quality is poor. e RMSE is defined as where V is the model velocity, X is the velocity of the traveltime tomography, and r×s is the model size. e value of UQI is in the range of [− 1, 1], and when the value is closer to 1, the imaging quality will be better. e UQI is defined as where , V is the model velocity, V is the mean velocity of the model, X is the velocity of the tomography results, and r×s is the model size.
e RMSE and UQI of the traveltime tomography results with the ITV regularization, ATV regularization, and MATV regularization are shown in Table 1. Compared with the other methods, the RMSE of the MATV regularization method is smaller, which indicates a better imaging quality than the other methods. e UQI of the MATV regularization method is greater than the UQI of the other methods, which also shows that its imaging quality is better than the other methods.
Two aspects are used to evaluate a tomography method. On the one hand, the tomography results should be basically consistent with the true model. On the other hand, the data fitting errors should be as small as possible. erefore, in addition to the tomography results, the errors between the synthetic data and the observed data in the data domain should also be considered, and their value should be as small as possible.
To illustrate the performance of the MATV regularization method in the data domain, we use the results of the MATV, as shown in Figure 7(d), to perform forward modelling and obtain MATV synthetic data. Next, we compare the MATV synthetic data with the first-arrival traveltime data calculated by the FSM with the model in Figure 3(a). Figure 10(a) shows a comparison between the first-arrival traveltime data (observations) and the MATV synthetic data (synthetic), which shows a small difference. To quantify the difference, the residuals are obtained by subtracting the synthetic data from the observations. e histogram of the residuals (Figure 10(b)) shows that the residuals are distributed in the range of [− 6.0×10 − 4 s, +6.0×10 − 4 s], which suggests that the residuals of the MATV synthetic data and observations are very small.
Compared with the ITV regularization and ATV regularization methods, the advantage of the MATV regularization method proposed in this paper is that it not only better maintains the boundary of the model but also effectively reduces the artifacts in the inversion due to the full use of multidirectional gradient information.

Real Data Example
We apply the proposed MATV regularization method to the first-arrival traveltime tomography of the transmitted waves in a coal mining face. Figure 11 shows the observation system of seismic transmission exploration in the coal mining face with a width of 143 m. Shots are placed in one roadway, and receivers are placed in another roadway to receive the transmitted waves. e shots and receivers are located in the middle of the coal seam. e main purpose of this study is to conduct tomography of the surrounding rock of a coal mining face by picking up the first-arrival of the transmitted waves. In this observational system, the shot spacing is 10 m, and there are 29 shots in total; the receiver spacing is 10 m, and there are 24 receivers per shot. A total of 667 first-arrival data points are selected because the 16th receiver is not working. e first-arrival traveltime of the 12th shot is shown in Figure 12. Figure 13(a)-13(c) shows the results of the first-arrival traveltime tomography with the ITV regularization, ATV regularization, and MATV regularization. e cool color (blue) is the distribution area of the geological anomalies that represent the low-velocity value, while the hot color (red) represents the high-velocity value in Figure 13. F1, F2, F3, F4, and F5 are the exposed faults in the surrounding rock of the coal mining face, and YC is the predicted geological anomaly area. e results of the several methods are similar in Figure 13. However, they exhibit some differences after a

R1
R5 R10 R15 R20 S1 S5 S9 S13 S17 S21 S25 S29  detailed analysis. e results of the ITV regularization are the worst. e low-velocity distribution areas are smaller than the fault extension length in Figure 13(a), and there are no low-velocity distribution areas at F5. e ATV regularization results show that the low-velocity distribution areas are slightly smaller than the fault extension length. As shown in Figure 13(b), we can determine the location of the faults in the results of the ATV regularization. e results of the MATV regularization show that the low-velocity distribution areas have the best matching relationship with the fault extension length and distribution pattern (Figure 13(c)).
e YC area is the predicted geological anomaly area, and we predict the geological anomalies in the undiscovered area according to our inversion results and the actual exposed data. We think that there may be geological anomalies in the area and remind the coal mine staff to take safety measures to ensure production safety. To illustrate the performance of the MATV regularization method in the data domain, we use the MATV regularization results (Figure 13(c)) to perform forward modelling and obtain the MATV synthetic data. We compare the synthetic data with the observations, which are the first-arrival traveltime of the transmitted waves. ere are certain differences when comparing the MATV synthetic data and the observations (Figure 14(a)). To quantify the difference between the MATV synthetic data and the observation data, the residuals are obtained by subtracting the MATV synthetic data from the observation data. Figure 14(b) shows the histogram of the residuals, which are mainly distributed in the range of [− 0.01 s, +0.01 s]. Most values of the residuals are near zero, which shows the better effect of the MATV regularization method in the data space.
From the histograms of the residuals of the simulated data example and the real data example (Figures 10(b) and 14(b)), it can be seen that the distributions of the residuals may be different for different imaging models. In Figures 10(b) and 14(b), the residuals between the MATV synthetic data and the observed data are relatively small, which shows that the MATV have better performance in the data domain.

Conclusions
is paper proposes a first-arrival traveltime tomography for the surrounding rock of a coal mining face based on the MATV regularization. We first propose two-dimensional and three-dimensional MATV by using the multidirectional gradient information of the pixels. en, we establish the objective function based on the MATV regularization. e objective function is decoupled into two subproblems that are more easily solved than the original one through the split Bregman method. Next, we use the conjugate gradient method and the shrinkage operator to efficiently solve these two subproblems. Finally, we apply the MATV regularization method to simulated data and real data and compare its performance with the ITV regularization and ATV regularization. e simulated data example shows that the MATV regularization method preserves the boundaries of geological anomalies more effectively than the other methods and reduces the artifacts. In addition, the RMSE of the MATV regularization method is smaller than the RMSE of the other methods, and its UQI is larger than the UQI of the other methods, which indicates that the imaging quality is higher for the MATV regularization method than for the other methods.
e real data example shows that the MATV regularization method can better match the actual geological anomalies that are revealed and can better describe the distribution range of the geological anomalies. Moreover, the MATV proposed in this paper can be easily extended to irregular grid, while multidirectional TV using four-directional gradient information may be difficult to irregular grid.
Nevertheless, the MATV regularization method proposed in this paper still has some shortcomings. At present, we only study the two-dimensional model; further studies are required to explore the three-dimensional model. As many improved TV methods, the MATV method increases the computational cost more or less. However, the added computational cost is acceptable [11,35]. In addition, the adaptive selection principle of the regularization factor must be explored further.

Data Availability
All the data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.