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

Authors: Bo Chen; Guowei Zhu; Zhenqiang Yang
Journal: Shock and Vibration
 (2021)

---

## Abstract

The 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. The surrounding rock of coal mining faces uses principles of tomography similar to those of medical X-rays. The 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 the MATV regularization and then used the split Bregman method to solve the objective function. The 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.


---

## Body

## 1. Introduction

X-ray computed tomography (CT) has been widely used in medical diagnosis [1–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. The 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]. There 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–11], image reconstruction [12, 13], and geophysical inversion [14–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. These improvements can be divided into the following categories: (1) generalize the difference of operation in TV to a higher-order difference [18–20] or fractional-order difference [21] such as high-order TV [18–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–32]; (4) use Tikhonov regularization and TV regularization to solve the objective function in an alternative way, such as modified TV [33–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–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. The refracted waves that propagate along the surrounding rock will be generated when the seismic waves are excited in the coal mining face. The velocity of the refracted waves that reach the receivers first is the fastest. Therefore, 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. Through 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. The results of the two examples show the effectiveness of the proposed method.


## 2. Multidirectional Anisotropic Total Variation

As shown in Figure1(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]:(1)PATV=∑m,jPm,j−Pm,j+1+Pm,j−Pm+1,j=∇xP1+∇yP1,where P is the pixel matrix, Pm,j is the element of the matrix P, and ∇xP and ∇yP represent the first-order difference of the pixels in the x direction and y direction, respectively. The 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. ∇x=−10⋯0︷r−110⋯⋯0−10⋯010⋯⋮⋮⋱⋱⋯⋱⋱⋯00⋯−10⋯01r×s−1×r×s, and ∇y=−110⋯⋯0−110⋯⋮⋮⋱⋱⋯00⋯−11r×s−1×r×s.Figure 1
Schematic model. (a) Two-dimensional pixel and its 2-direction gradient map. (b) Two-dimensional pixel and its 8-direction gradient map. (c) Three-dimensional voxel and its 26-direction gradient map.
(a)(b)(c)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(2)PMATV_2d=∑m,jPm,j−Pm−1,j−1+Pm,j−Pm,j−1+Pm,j−Pm+1,j−1+Pm,j−Pm−1,j+Pm,j−Pm+1,j+Pm,j−Pm−1,j+1+Pm,j−Pm,j+1+Pm,j−Pm+1,j+1=∇1P1+∇2P1+∇3P1+∇4P1+∇5P1+∇6P1+∇7P1+∇8P1,where ∇1P, ∇2P, ∇3P, ∇4P, ∇5P, ∇6P, ∇7P, and ∇8P represent the first-order difference of the pixel in 8 adjacent directions.Generalized to three-dimensional images, as shown in Figure1(c), three-dimensional MATV uses the gradient information of a voxel and its 26 adjacent voxels, and its specific expression is defined as(3)PMATV_3d=∑m,j,nPm,j,n−Pm+1,j−1,n+1+Pm,j,n−Pm+1,j−1,n+Pm,j,n−Pm+1,j−1,n−1+Pm,j,n−Pm+1,j,n+1+Pm,j,n−Pm+1,j,n+Pm,j,n−Pm+1,j,n−1+Pm,j,n−Pm+1,j+1,n+1+Pm,j,n−Pm+1,j+1,n+Pm,j,n−Pm+1,j+1,n−1+Pm,j,n−Pm,j−1,n+1+Pm,j,n−Pm,j−1,n+Pm,j,n−Pm,j−1,n−1+Pm,j,n−Pm,j,n+1+Pm,j,n−Pm,j,n−1+Pm,j,n−Pm,j+1,n+1+Pi,j,k−Pm,j+1,n+Pm,j,n−Pm,j+1,n−1+Pm,j,n−Pm−1,j−1,n+1+Pm,j,n−Pm−1,j−1,n+Pm,j,n−Pm−1,j−1,n−1+Pm,j,n−Pm−1,j,n+1+Pm,j,n−Pm−1,j,n+Pm,j,n−Pm−1,j,n−1+Pm,j,n−Pm−1,j+1,n+1+Pm,j,n−Pm−1,j+1,n+Pm,j,n−Pm−1,j+1,n−1=∇1P1+∇2P1+∇3P1+∇4P1+∇5P1+∇6P1+∇7S1+∇8P1+∇9P1+∇10P1+∇11P1+∇12P1+∇13P1+∇14P1+∇15P1+∇16P1+∇17P1+∇18P1+∇19P1+∇20P1+∇21P1+∇22P1+∇23P1+∇24P1+∇25P1+∇26P1,where ∇1P, ∇2P, ∇3P, ∇4P, ∇5P, ∇6P, ∇7P, ∇8P, ∇9P, ∇10P, ∇11P, ∇12P, ∇13P, ∇14P, ∇15P, ∇16P, ∇17P, ∇18P, ∇19P, ∇20P, ∇21P, ∇22P, ∇23P, ∇24P, ∇25P, and ∇26P represent the first-order difference of the voxel in 26 adjacent directions.


## 3. MATV Regularization Tomography for the Surrounding Rock of a Coal Mining Face

### 3.1. Transmitted Waves of a Coal Mining Face

A typical three-layer symmetric model of rock-coal-rock is shown in Figure2(a). The velocity of the coal seam is generally lower than the velocity of the surrounding rock. The 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. Through transmitted wave analysis, we can obtain the geological anomalies of the coal mining face.Figure 2
Transmitted waves in the coal mining face. (a) The three-layer symmetric coal seam model. (b) Observation system of seismic transmission exploration in a coal mining face.S is a shot point in Roadway 1. R1–Rn are receivers in Roadway 2. Point O is the projection of S in Roadway 2. L is the distance between Roadway 1 and Roadway 2. X is the distance from point O to the receivers. (c) Ray path of different seismic waves for a symmetrical rock-coal-rock three-layer model. θ is the critical angle of the coal seam. Vrock is the velocity of the surrounding rock, and Vcoal is the velocity of the coal seam. The thickness of the coal seam is d. Both shot points and receivers are placed at d/2. r is the distance between the incident point of the refracted waves and the shot.
(a)(b)(c)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 Figure2(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. Therefore, 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. Therefore, 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. The 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(c) is from 1.15 m to 4.24 m. The traveltime from the shot to the incident point can be calculated if the velocity of the coal mine is known. Then, we can correct the first-arrival traveltime of the transmitted waves.


### 3.2. Establishment of the Objective Function

In the two-dimensional plane ofx and y, if the path of the transmitted wave is L and the velocity of the medium is Vx,y, then the traveltime T of the transmitted waves can be expressed as(4)T=∫LdLVx,y.When the medium velocity changes, the first-arrival traveltime is a nonlinear function of the medium velocity, which can be expressed as follows:(5)T=LS,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 ofLS is written as(6)LS≈LSk+L′SkS−Sk,where Sk is the slowness of the model at the k-th iteration. Let L′Sk = Ak, where Ak is the length matrix of the ray on the model grid at the k-th iteration. The backward process of linear traveltime interpolation (LTI) [42] is used to solve Ak in this paper. LSk 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 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:(7)argminSd−LS22+λSMATV_2d,where d is the first-arrival traveltime, LS is the theoretical traveltime calculated by the model’s slowness S, λ is the regularization factor, ‖ ‖2 represents the errors between d and LS, and argminS indicates the value of S when the objective function is minimized.Substituting equation (6) into the objective function (7) yields(8)argminSd−LSk+AkS−Sk22+λSMATV_2d.


### 3.3. Solution Method

The split Bregman algorithm [44] is used to solve equation (8). First, auxiliary variables are introduced, namely, ζ1=∇1S, ζ2=∇2S, ζ3=∇3S, ζ4=∇4S, ζ5=∇5S, ζ6=∇6S, ζ7=∇7S, and ζ8=∇8S. Equation (8) is transformed into(9)argminS,ζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8d−LSk+AkS−Sk22+λSMATV_2d+αζ1−∇1S22+αζ2−∇2S22+αζ3−∇3S22+αζ4−∇4S22+αζ5−∇5S22+αζ6−∇6S22+αζ7−∇7S22+αζ8−∇8S22.According to Goldstein and Osher [44], the regularization factor α is set to α=2λ.By introducing the intermediate variables ofq1, q2, q3, q4, q5, q6, q7, and q8, equation (9) is transformed into(10)argminS,ζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8d−LSk+AkS−Sk22+λSMATV_2d+αζ1−∇1S−q1k22+αζ2−∇2S−q2k22+αζ3−∇3S−q3k22+αζ4−∇4S−q4k22+αζ5−∇5S−q5k22+αζ6−∇6S−q6k22+αζ7−∇7S−q7k22++αζ8−∇8S−q8k22,where qik+1=qik+∇iSk+1−ζik+1 and qi0=0 for i=1,2,3,…,8, and k is the number of iterations.Equation (10) can be decomposed into the following two subproblems:(11)argminSd−LSk+AkS−Sk22+αζ1k−∇1S−q1k22+αζ2k−∇2S−q2k22+αζ3k−∇3S−q3k22+αζ4k−∇4S−q4k22+αζ5k−∇5S−q5k22+αζ6k−∇6S−q6k22+αζ7k−∇7S−q7k22+αζ8k−∇8S−q8k22,(12)argminζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8λSMATV_2d+αζ1−∇1S−q1k22+αζ2−∇2S−q2k22+αζ3−∇3S−q3k22+αζ4−∇4S−q4k22+αζ5−∇5S−q5k22+αζ6−∇6S−q6k22+αζ7−∇7S−q7k22+αζ8−∇8S−q8k22.We set the first derivative of the first subproblem (11) to zero [45] and then obtain(13)AΔSk=B,where A=AkTAk+α∇1T∇1+α∇2T∇2+α∇3T∇3+α∇4T∇4+α∇5T∇5+α∇6T∇6+α∇7T∇7+α∇8T∇8, ΔSk=Sk+1−Sk, and B=AkTd−LSk+α∇1Tζ1k−q1k+α∇2Tζ2k−q2k+α∇3Tζ3k−q3k+α∇4Tζ4k−q4k+α∇5Tζ5k−q5k+α∇6Tζ6k−q6k+α∇7Tζ7k−q7k+α∇8Tζ8k−q8k−α∇1T∇1+α∇2T∇2+α∇3T∇3+α∇4T∇4+α∇5T∇5+α∇6T∇6+α∇7T∇7+α∇8T∇8Sk. Equation (13) can be efficiently solved by using the conjugate gradient method.The second subproblem can be solved with a shrinkage operator [46]:(14)ζik+1=shrink∇iSk+1+qik,λ2α,i=1,2,3,…8,where shrink is defined as shrinkθ,α=signθmaxθ−α,0.The pseudocode of using the split Bregman algorithm to solve equation (8) is shown in Algorithm 1.Algorithm 1:First-arrival traveltime tomography with the MATV regularization.
Input:λ, ε, N, d, S0Output:S.(1)
Initializek=0, S0=S0, ζi=qi=0, i = 1, 2, 3, …, 8;(2)
whilek<Ndo(3)
Sk+1=argminSd−LSk+AkS−Sk22+αζ1k−∇1S−q1k22+αζ2k−∇2S−q2k22+αζ3k−∇3S−q3k22+αζ4k−∇4S−q4k22+αζ5k−∇5S−q5k22+αζ6k−∇6S−q6k22+αζ7k−∇7S−q7k22+αζ8k−∇8S−q8k22;(4)
ζik+1=shrink∇iSk+1+qik,λ/2α, i = 1, 2, 3, …, 8;(5)
qik+1=qik+∇iSk+1−ζik+1, i = 1, 2, 3, …, 8;(6)
k=k+1;(7)
ifSk+1−Sk2/Sk2>εthen(8)
Sk = Sk+1;(9)
else(10)
returnSk+1;(11)
end(12)
end


## 3.1. Transmitted Waves of a Coal Mining Face

A typical three-layer symmetric model of rock-coal-rock is shown in Figure2(a). The velocity of the coal seam is generally lower than the velocity of the surrounding rock. The 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. Through transmitted wave analysis, we can obtain the geological anomalies of the coal mining face.Figure 2
Transmitted waves in the coal mining face. (a) The three-layer symmetric coal seam model. (b) Observation system of seismic transmission exploration in a coal mining face.S is a shot point in Roadway 1. R1–Rn are receivers in Roadway 2. Point O is the projection of S in Roadway 2. L is the distance between Roadway 1 and Roadway 2. X is the distance from point O to the receivers. (c) Ray path of different seismic waves for a symmetrical rock-coal-rock three-layer model. θ is the critical angle of the coal seam. Vrock is the velocity of the surrounding rock, and Vcoal is the velocity of the coal seam. The thickness of the coal seam is d. Both shot points and receivers are placed at d/2. r is the distance between the incident point of the refracted waves and the shot.
(a)(b)(c)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 Figure2(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. Therefore, 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. Therefore, 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. The 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(c) is from 1.15 m to 4.24 m. The traveltime from the shot to the incident point can be calculated if the velocity of the coal mine is known. Then, we can correct the first-arrival traveltime of the transmitted waves.


## 3.2. Establishment of the Objective Function

In the two-dimensional plane ofx and y, if the path of the transmitted wave is L and the velocity of the medium is Vx,y, then the traveltime T of the transmitted waves can be expressed as(4)T=∫LdLVx,y.When the medium velocity changes, the first-arrival traveltime is a nonlinear function of the medium velocity, which can be expressed as follows:(5)T=LS,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 ofLS is written as(6)LS≈LSk+L′SkS−Sk,where Sk is the slowness of the model at the k-th iteration. Let L′Sk = Ak, where Ak is the length matrix of the ray on the model grid at the k-th iteration. The backward process of linear traveltime interpolation (LTI) [42] is used to solve Ak in this paper. LSk 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 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:(7)argminSd−LS22+λSMATV_2d,where d is the first-arrival traveltime, LS is the theoretical traveltime calculated by the model’s slowness S, λ is the regularization factor, ‖ ‖2 represents the errors between d and LS, and argminS indicates the value of S when the objective function is minimized.Substituting equation (6) into the objective function (7) yields(8)argminSd−LSk+AkS−Sk22+λSMATV_2d.


## 3.3. Solution Method

The split Bregman algorithm [44] is used to solve equation (8). First, auxiliary variables are introduced, namely, ζ1=∇1S, ζ2=∇2S, ζ3=∇3S, ζ4=∇4S, ζ5=∇5S, ζ6=∇6S, ζ7=∇7S, and ζ8=∇8S. Equation (8) is transformed into(9)argminS,ζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8d−LSk+AkS−Sk22+λSMATV_2d+αζ1−∇1S22+αζ2−∇2S22+αζ3−∇3S22+αζ4−∇4S22+αζ5−∇5S22+αζ6−∇6S22+αζ7−∇7S22+αζ8−∇8S22.According to Goldstein and Osher [44], the regularization factor α is set to α=2λ.By introducing the intermediate variables ofq1, q2, q3, q4, q5, q6, q7, and q8, equation (9) is transformed into(10)argminS,ζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8d−LSk+AkS−Sk22+λSMATV_2d+αζ1−∇1S−q1k22+αζ2−∇2S−q2k22+αζ3−∇3S−q3k22+αζ4−∇4S−q4k22+αζ5−∇5S−q5k22+αζ6−∇6S−q6k22+αζ7−∇7S−q7k22++αζ8−∇8S−q8k22,where qik+1=qik+∇iSk+1−ζik+1 and qi0=0 for i=1,2,3,…,8, and k is the number of iterations.Equation (10) can be decomposed into the following two subproblems:(11)argminSd−LSk+AkS−Sk22+αζ1k−∇1S−q1k22+αζ2k−∇2S−q2k22+αζ3k−∇3S−q3k22+αζ4k−∇4S−q4k22+αζ5k−∇5S−q5k22+αζ6k−∇6S−q6k22+αζ7k−∇7S−q7k22+αζ8k−∇8S−q8k22,(12)argminζ1,ζ2,ζ3,ζ4,ζ5,ζ6,ζ7,ζ8λSMATV_2d+αζ1−∇1S−q1k22+αζ2−∇2S−q2k22+αζ3−∇3S−q3k22+αζ4−∇4S−q4k22+αζ5−∇5S−q5k22+αζ6−∇6S−q6k22+αζ7−∇7S−q7k22+αζ8−∇8S−q8k22.We set the first derivative of the first subproblem (11) to zero [45] and then obtain(13)AΔSk=B,where A=AkTAk+α∇1T∇1+α∇2T∇2+α∇3T∇3+α∇4T∇4+α∇5T∇5+α∇6T∇6+α∇7T∇7+α∇8T∇8, ΔSk=Sk+1−Sk, and B=AkTd−LSk+α∇1Tζ1k−q1k+α∇2Tζ2k−q2k+α∇3Tζ3k−q3k+α∇4Tζ4k−q4k+α∇5Tζ5k−q5k+α∇6Tζ6k−q6k+α∇7Tζ7k−q7k+α∇8Tζ8k−q8k−α∇1T∇1+α∇2T∇2+α∇3T∇3+α∇4T∇4+α∇5T∇5+α∇6T∇6+α∇7T∇7+α∇8T∇8Sk. Equation (13) can be efficiently solved by using the conjugate gradient method.The second subproblem can be solved with a shrinkage operator [46]:(14)ζik+1=shrink∇iSk+1+qik,λ2α,i=1,2,3,…8,where shrink is defined as shrinkθ,α=signθmaxθ−α,0.The pseudocode of using the split Bregman algorithm to solve equation (8) is shown in Algorithm 1.Algorithm 1:First-arrival traveltime tomography with the MATV regularization.
Input:λ, ε, N, d, S0Output:S.(1)
Initializek=0, S0=S0, ζi=qi=0, i = 1, 2, 3, …, 8;(2)
whilek<Ndo(3)
Sk+1=argminSd−LSk+AkS−Sk22+αζ1k−∇1S−q1k22+αζ2k−∇2S−q2k22+αζ3k−∇3S−q3k22+αζ4k−∇4S−q4k22+αζ5k−∇5S−q5k22+αζ6k−∇6S−q6k22+αζ7k−∇7S−q7k22+αζ8k−∇8S−q8k22;(4)
ζik+1=shrink∇iSk+1+qik,λ/2α, i = 1, 2, 3, …, 8;(5)
qik+1=qik+∇iSk+1−ζik+1, i = 1, 2, 3, …, 8;(6)
k=k+1;(7)
ifSk+1−Sk2/Sk2>εthen(8)
Sk = Sk+1;(9)
else(10)
returnSk+1;(11)
end(12)
end


## 4. 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 low-velocity 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.


### 4.1. Curved Ray-Tracing

As shown in Figure3(a), the size of the model is 100 m × 100 m, and it contains two local low-velocity blocks (geological anomalies). The 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.Figure 3
Model and observation system. (a) Model with two local low-velocity blocks. (b) Model grid division and observation system with seismic transmission exploration (the red dotted lines are the positions of the low-velocity blocks).
(a)(b)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 ray-tracing, we design the model with the size of 500 m × 500 m and a velocity of 2000 m/s.Figure4(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 4
Comparison of the traveltime. (a) Contour map of the traveltime with the FSM and theoretical. (b) Error map of the traveltime with the FSM and theoretical.
(a)(b)Figure5 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.Figure 5
Comparison of the ray paths with the backward process of LTI and theoretical.Figure6(a) is the traveltime contour map of the 25th shot calculated by the FSM for Figure 3(a). The 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. The ray paths bypass the low-velocity blocks to find the minimum traveltime path and present a curved shape.Figure 6
Curved ray-tracing of the 25th shot. (a) The traveltime contours of Figure3(a) are calculated by the FSM. (b) The traveltime diagram of receivers R1–R51. (c) The ray paths traced by the backward process of LTI (the red dotted lines are the positions of the low-velocity blocks).
(a)(b)(c)


### 4.2. Comparison of Several Tomography Methods

The first-arrival traveltime data are calculated by the FSM with the model in Figure3(a). The initial value is the result of the back projection technique (BPT) (Figure 7(a)). Then, 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)).Figure 7
Traveltime tomography results with (a) BPT, (b) ITV regularization, (c) ATV regularization, and (d) MATV regularization (the white dotted lines are the positions of the low-velocity blocks).
(a)(b)(c)(d)In the TV-based method, it is important to select the appropriate regularization factors. The 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 Figure7(b)–7(d), all methods can predict the geological anomalies to a certain extent. The 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. The 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. Thus, compared with the ITV regularization and ATV regularization, the proposed MATV regularization method has high imaging accuracy and can improve image quality.To intuitively compare the several traveltime tomography results in Figure7, we extract the results at x = 28 m and x = 68 m for analysis, as shown in Figure 8. The 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.Figure 8
Comparison of the traveltime tomography results. (a) Velocity profile of the traveltime tomography results atx = 28 m. (b) Velocity profile of the traveltime tomography results at x = 68 m.
(a)(b)The tomography errors are obtained by subtracting the model shown in Figure3(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. That is to say, compared with the other two methods, the tomography results of the MATV regularization method have more points equal to the model. This indicates that the errors of the MATV regularization method are smaller than the errors of the other two methods at most points.Figure 9
Histogram curves of the tomography errors.


### 4.3. Quantitative Evaluation

To evaluate the traveltime tomography results quantitatively, the universal quality index (UQI) [47] and root mean square error (RMSE) are introduced. The 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. The RMSE is defined as(15)RMSE=∑j=1r∑m=1sVj,m−Xj,m2r×s,where V is the model velocity, X is the velocity of the traveltime tomography, and r×s is the model size.The value of UQI is in the range of [−1, 1], and when the value is closer to 1, the imaging quality will be better. The UQI is defined as(16)UQI=4σVXV¯X¯σV2+σX2V¯2+X¯2,where V¯=1/r×s∑j=1r∑m=1sVjm, X¯=1/r×s∑j=1r∑m=1sXjm, σV2=1/r×s−1∑j=1r∑m=1sVjm−V¯2, σX2=1/r×s−1∑j=1r∑m=1sXjm−X¯2, σVX=1/r×s−1∑j=1r∑m=1sVjm−V¯Xjm−X¯, 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.The RMSE and UQI of the traveltime tomography results with the ITV regularization, ATV regularization, and MATV regularization are shown in Table1. Compared with the other methods, the RMSE of the MATV regularization method is smaller, which indicates a better imaging quality than the other methods. The 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.Table 1
The RMSE and UQI of the ITV, ATV, and MATV.
MethodsITVATVMATVRMSE239.50232.34212.87UQI0.490.490.64Two 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. Therefore, 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 Figure7(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. The histogram of the residuals (Figure 10(b)) shows that the residuals are distributed in the range of [−6.0×10−4s, +6.0×10−4s], which suggests that the residuals of the MATV synthetic data and observations are very small.Figure 10
Data comparison in the data field. (a) Comparison of MATV synthetic data and observations. (b) Histogram of the residuals of MATV synthetic data and observations.
(a)(b)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.


## 4.1. Curved Ray-Tracing

As shown in Figure3(a), the size of the model is 100 m × 100 m, and it contains two local low-velocity blocks (geological anomalies). The 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.Figure 3
Model and observation system. (a) Model with two local low-velocity blocks. (b) Model grid division and observation system with seismic transmission exploration (the red dotted lines are the positions of the low-velocity blocks).
(a)(b)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 ray-tracing, we design the model with the size of 500 m × 500 m and a velocity of 2000 m/s.Figure4(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 4
Comparison of the traveltime. (a) Contour map of the traveltime with the FSM and theoretical. (b) Error map of the traveltime with the FSM and theoretical.
(a)(b)Figure5 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.Figure 5
Comparison of the ray paths with the backward process of LTI and theoretical.Figure6(a) is the traveltime contour map of the 25th shot calculated by the FSM for Figure 3(a). The 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. The ray paths bypass the low-velocity blocks to find the minimum traveltime path and present a curved shape.Figure 6
Curved ray-tracing of the 25th shot. (a) The traveltime contours of Figure3(a) are calculated by the FSM. (b) The traveltime diagram of receivers R1–R51. (c) The ray paths traced by the backward process of LTI (the red dotted lines are the positions of the low-velocity blocks).
(a)(b)(c)


## 4.2. Comparison of Several Tomography Methods

The first-arrival traveltime data are calculated by the FSM with the model in Figure3(a). The initial value is the result of the back projection technique (BPT) (Figure 7(a)). Then, 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)).Figure 7
Traveltime tomography results with (a) BPT, (b) ITV regularization, (c) ATV regularization, and (d) MATV regularization (the white dotted lines are the positions of the low-velocity blocks).
(a)(b)(c)(d)In the TV-based method, it is important to select the appropriate regularization factors. The 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 Figure7(b)–7(d), all methods can predict the geological anomalies to a certain extent. The 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. The 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. Thus, compared with the ITV regularization and ATV regularization, the proposed MATV regularization method has high imaging accuracy and can improve image quality.To intuitively compare the several traveltime tomography results in Figure7, we extract the results at x = 28 m and x = 68 m for analysis, as shown in Figure 8. The 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.Figure 8
Comparison of the traveltime tomography results. (a) Velocity profile of the traveltime tomography results atx = 28 m. (b) Velocity profile of the traveltime tomography results at x = 68 m.
(a)(b)The tomography errors are obtained by subtracting the model shown in Figure3(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. That is to say, compared with the other two methods, the tomography results of the MATV regularization method have more points equal to the model. This indicates that the errors of the MATV regularization method are smaller than the errors of the other two methods at most points.Figure 9
Histogram curves of the tomography errors.


## 4.3. Quantitative Evaluation

To evaluate the traveltime tomography results quantitatively, the universal quality index (UQI) [47] and root mean square error (RMSE) are introduced. The 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. The RMSE is defined as(15)RMSE=∑j=1r∑m=1sVj,m−Xj,m2r×s,where V is the model velocity, X is the velocity of the traveltime tomography, and r×s is the model size.The value of UQI is in the range of [−1, 1], and when the value is closer to 1, the imaging quality will be better. The UQI is defined as(16)UQI=4σVXV¯X¯σV2+σX2V¯2+X¯2,where V¯=1/r×s∑j=1r∑m=1sVjm, X¯=1/r×s∑j=1r∑m=1sXjm, σV2=1/r×s−1∑j=1r∑m=1sVjm−V¯2, σX2=1/r×s−1∑j=1r∑m=1sXjm−X¯2, σVX=1/r×s−1∑j=1r∑m=1sVjm−V¯Xjm−X¯, 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.The RMSE and UQI of the traveltime tomography results with the ITV regularization, ATV regularization, and MATV regularization are shown in Table1. Compared with the other methods, the RMSE of the MATV regularization method is smaller, which indicates a better imaging quality than the other methods. The 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.Table 1
The RMSE and UQI of the ITV, ATV, and MATV.
MethodsITVATVMATVRMSE239.50232.34212.87UQI0.490.490.64Two 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. Therefore, 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 Figure7(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. The histogram of the residuals (Figure 10(b)) shows that the residuals are distributed in the range of [−6.0×10−4s, +6.0×10−4s], which suggests that the residuals of the MATV synthetic data and observations are very small.Figure 10
Data comparison in the data field. (a) Comparison of MATV synthetic data and observations. (b) Histogram of the residuals of MATV synthetic data and observations.
(a)(b)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.


## 5. 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. Figure11 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. The shots and receivers are located in the middle of the coal seam. The 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. The first-arrival traveltime of the 12th shot is shown in Figure 12.Figure 11
Observation system of seismic transmission exploration in a coal mining face.Figure 12
(a) Transmitted waves’ record and first-arrival traveltime of the 12th shot in a coal mining face. The red dots are the first-arrival traveltime picked manually. (b) Picked the first-arrival traveltime.
(a)(b)Figure13(a)–13(c) shows the results of the first-arrival traveltime tomography with the ITV regularization, ATV regularization, and MATV regularization. The 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.Figure 13
The results of the first-arrival traveltime tomography with (a) ITV regularization, (b) ATV regularization, and (c) MATV regularization.
(a)(b)(c)The results of the several methods are similar in Figure13. However, they exhibit some differences after a detailed analysis. The results of the ITV regularization are the worst. The 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. The 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. The 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)).The 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 (Figure13(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. There 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.Figure 14
Data comparison in the data field. (a) Comparison of MATV synthetic data and observations. (b) Histogram of the residuals of MATV synthetic data and observations.
(a)(b)From the histograms of the residuals of the simulated data example and the real data example (Figures10(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.


## 6. Conclusions

This 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. Then, we establish the objective function based on the MATV regularization. The 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. The 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. The 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.



---

*Source: 1138620-2021-09-07.xml*