A SOURCE TIME REVERSAL METHOD FOR SEISMICITY INDUCED BY MINING

. In this work, we present a modiﬁed Time-Reversal Mirror (TRM) Method, called Source Time Reversal (STR), to ﬁnd the spatial distribution of a seismic source induced by mining activity. This methodology is based on a known full description of the temporal dependence of the source, the Duhamel’s principle, and the time-reverse property of the wave equation. We also provide an error estimate of the reconstruction when the measurements are acquired over the entire boundary, and we show experimentally the inﬂuence of measuring on a subdomain of the boundary. Numerical results indicate that the methodology is able to recover continuous and discontinuous sources, and it remains stable for partial boundary measurements.


1.
Introduction. As mining is a very important activity around the world and also one of the principal economic activities in several countries such as Chile, Peru, South Africa, and China, it is expected to continue for a long time. In this dangerous labor, many precautions are considered in order to prevent mining accidents. In several cases, the lack of safety ends on accidents ranging from minor injuries to death of miners. For example, India had in the 90's around 1.100 serious accidents and 230 fatal accidents per year [23]; in Spain, 70.000 fatal and non-fatal accidents occurred between 2003 and 2012 [27]; and 28.868 accidents with 47.875 deaths took place on coal mines in China between 2001 and 2011 [7]. In addition, numerous studies have been performed to measure and try to prevent accidents, e.g. [18,7,8,32]. Rock failures (or rockbursts) are the principal cause of those accidents in mines [16].
Mining is a labor that for its nature induces seismic activity. Due to the work inside mines, the rocks are stressed, often producing tremors or microseisms (in reference to its magnitude), a process that is known as block caving. This type of seismic activity is named induced seismicity. Knowing the sources of induced microseisms and studying the hazard in those zones is useful to improve the miners safety inside mines. In this work, we focus on the characterization of sources in induced seismicity.
There are no systematic differences between seismicity and induced seismicity [16]. Thus, it is possible to use the same mathematical model for describing a seism and an induced seism. As models for seismicity and induced seismicity are the same, it is natural to apply seismicity source location techniques also to the case of induced seismicity.
A widely employed technique to determine source locations in traditional seismic methods is a time-reversal method introduced by M. Fink [13]. This method takes advantage of the time symmetries of waves, and it was introduced in 2006 for source location at the global scale with the work of Larmat et. al. [21]. In this work, they reversed seismogram measurements of the Sumatra Earthquake on time, refocusing the seismic wave energy on the earthquake source location.
On the other hand, the problem of locating the tremor source in induced seismicity is often studied in a simplified form as: to find simultaneously the time when the tremor began (origin time t 0 ) and the position where it started (hypocenter x 0 , where it is also often assumed that f (x) = δ(x − x 0 )). To do that, we look for an estimate of the wave travel time at each measurement station (travel time). More precisely, we only know the time at which the wave arrives at a measurement station (arrival time). Then, we have the following identity for the i-th station: arrival time i = origin time + travel time i .
In 1912, L. Geiger [15] developed one of the most classical methods for locating sources in mining, based on an accurate propagation velocity map of the mine structure and a least-square technique. This method is able to estimate the hypocenter location using the arrival time to measurement stations. Depending on the stations geometrical distribution, this method could present problems of accuracy or unicity. Later, in the 80s, W. Spence [30] modified Geiger's least-square method to consider the P-wave arrival time difference data, which allowed to eliminate anomalies produced by velocity heterogeneities. In 1984, M. Matsu'ura [24] proposed a method for locating the spatial-dependence of the source with a Bayesian approach using prior spatial information. In addition, the origin time is eliminated of the problem via integration on a density function over the whole origin time range. In another work, F. Du et. al. [9] considered the wave propagation velocity as a random variable into the minimization method for undersea mining microseismicity. Additionally, these authors divided the mine space in three areas based on the tremor occurrence probability to predict possible future hazard. As we can see, the methods for finding the source in seismicity induced by mining are more focused on finding the origin time by estimating the travel times, and then obtaining the spatial position.
The main difference between seismicity and induced seismicity for mining is the scale. In the former, scales are in the order of tens or hundreds of kilometers, while in the latter, scales are of one hundred meters or below. Related to the scale and the tremor duration, in seismicity, the source could be considered as a Dirac delta source, while in induced seismicity it is convenient to include the source duration effect in the mathematical model. More precisely, in mathematical terms, in seismicity is often assumed that the source s(x, t) is of the form s(x, t) = f (x)δ(t − t 0 ), where t 0 is the time when the seismic movement starts. In induced seismicity, the form of the time-dependence of the source may significantly affect to the models. For this reason, in here we only assume that source s(x, t) can be described using separation of variables as s(x, t) = f (x)g(t), but we do not make any further assumption on g(t). We note that this assumption of considering the source as f (x)g(t) has been used in hyperbolic problems (e.g. [34]) and parabolic problems (e.g. [14]) to obtain stability conditions and to reconstruct the spatial dependence of the source. This approach has also been used for medical imaging, see, e.g., chapter 12 of [3], where authors summarize time-reversal methods applied to tomography techniques by considering Dirac delta functions as source.
In contrast to previous approaches, we are going to assume that the timedependence form of the source g(t) is known, and it is not necessarily given by a Dirac delta function of the form δ(t − t 0 ). Works as [26] attempt to find an accurate representation of the amplitude of a seismic wave using Ricker wavelet, and if we include any method to estimate the origin time t 0 as [15,30], it is possible to obtain a complete characterization of the temporal source term g(t). Based on a full description of g(t), in this work we look for a more accurate representation of f (x) than that obtained when assuming that g(t) = δ(t − t 0 ), and this allows us to produce a more realistic mathematical model and a better reconstruction of f (x). In this work, for simplicity of the analysis of our proposed method, we restrict to acoustic equations instead of the more complex elastic equations. The former equations are known to provide adequate estimates of compressional velocities, and they are widely used as first approximations of the true compressional and shear velocities dictated by the elastic equations (see, e.g., [1]).
The main contribution of this work is to enhance the accuracy for the spatial source term reconstruction in induced seismicity. To do that, we adapt a timereversal method in order to find an approximation of the spatial term using the temporal source information and the geophones measurements. We denote this novel method as Source Time Reversal. In addition, we study numerically the properties of such source reconstruction technique. In particular, we show numerically that continuous sources are easier to reconstruct than discontinuous ones. While classical time-reversal methods allow to accurately recover initial conditions or sources acting for an instant of time (i.e. f (x)δ(t − t 0 )), the method developed in this work is able to recover the spatial characteristics for a more general class of sources f (x)g(t) where g(t) is assumed to be known.
This document is divided into five sections. In Section 2, we state the problem and introduce a modified Time-Reversal Mirror from M. Fink point of view; we also show an error estimate for our problem following the ideas of [17], and we provide a three dimensional example of how the time-reversal mirror method works in an exact case. In Section 3, we present the Source Time Reversal (STR) methodology. Section 4 introduces four different numerical experiments to study the characteristics of the method, including a synthetic microseismic experiment, where we compare the proposed method with the standard time-reversal mirror method. Finally in Section 5 we present our conclusions and some possible extensions.

Model assumptions and the time-reversal method.
In what follows, we assume our seismic events are governed by the wave equation in an isotropic media. The source can be expressed using separation of variables as f (x)g(t), where f is compactly supported in a bounded set Ω ⊂ R n . Thus, our original problem is given by: where c ∈ C ∞ (R n ), c(x) > c 1 > 0 for all x ∈ R n is a known propagation velocity map with (1 − c(x)) having compact support, and c(x) verifies the non-trapping condition (see Definition 2.1). Also, the temporal source distribution g(t) is assumed to be known, with g ∈ W 1,∞ (0, T ).
Definition 2.1 (non-trapping condition [10]). Let c : R n → R be a C ∞ function. We define the Hamiltonian H(x, ξ) = 1 2 c 2 (x)|ξ| 2 and the following system: The solutions of (1) are called bicharacteristics, and the projection of the xcomponents into R n of a bicharacteristic is called ray. We define the non-trapping condition as: all the rays go to infinity as t → ∞, that is lim t→∞ x(t) = +∞.
Let us consider a set of measurements obtained by using geophones, given by ∂ t u(x, t) recorded over the domain boundary ∂Ω for all time t ∈ [0, T ]. The objective is to reconstruct the spatial source contribution f (x). To mathematically express these measurements, we introduce the following operator: We denote the measurement function over the boundary as m u (x, t) = Λf (x, t).

Time-reversal mirror.
To reconstruct the spatial contribution of the source, we use a time-reversal method based on M. Fink ideas, which reverses and propagates the boundary measurements to recover an initial condition. Time-reversal mirror is a method that allows to reverse waves toward its origin (see, e.g. [13]). This method has been applied in several areas of physics and engineering, including optics [6,36,35], underwater acoustic [29,19], ultrasound [12,33], and wireless communications [25,20]. Time-reversal mirror is based on the invariance on time of the wave equation. This method allows to recover the origin of the wave or the initial wave 1 through measurements performed by transducers located at some fixed points inside the domain where the waves propagate. Each transducer records a signal and acts as a special mirror (reflecting first the last measurement and last the first measurement) by sending back the recorded signal. The reflected signal moves back towards the initial wave due to the time invariance of the wave equation. For a more detailed explanation on time-reversal mirror, see [13].
To explain the time-reversal mirror process used in this work, let us consider the following problem governed by a homogeneous wave equation in R n . Let v be the solution of the problem: for large values of |x|, and c(x) verifies the non-trapping condition. It is well known that the previous problem has a unique strong solution. Moreover, if we consider the above conditions, [5]).
In the more common case, the goal of the time-reversal mirror method is to recover the initial condition ϕ(x) or ψ(x) by measuring and recording v(x, t) on the boundary. In our case, the goal is to recover the initial condition ψ(x) by measuring and recording the displacement velocity ∂ t v(x, t) with geophones located on all ∂Ω. To recover the initial condition ψ(x), we define the following operator to model the measurement: This operator records the displacement velocity over the boundary. The measurements are given by m(x, t) = Λ(ϕ, ψ)(x, t) for all x ∈ ∂Ω and for all t ∈ (0, T ). Then, we define the following boundary value problem The solution of equation (3) is w(x, t) = ∂ t v(x, t), and at time t = 0, we obtain w(x, 0) = ψ(x). Then, in order to recover exactly the initial condition ψ(x), we need to know the values of ∂ t v(x, T ) and ∂ 2 t v(x, T ) inside Ω. In practice, it is not possible to know ∂ t v(x, T ) and ∂ 2 t v(x, T ) inside Ω. However, for a sufficiently large T , due to the local energy decay (see Theorem 3.1), we can approximate ∂ t v(x, T ) and ∂ 2 t v(x, T ) by zero. It is feasible to solve problem (3) with final data equal to zero and the measurements m(x, t) as boundary condition and recover exactly ψ(x) with w(x, 0) when T → ∞. In the case when T < ∞, we can recover exactly ψ(x) only for odd dimensions n and constant propagation velocity (Huygens' principle, see e.g. [22]). For even dimensions or variable propagation velocity, it is only possible to approximate ψ(x) by assuming the initial data equal to zero in problem (3).
We now define the problem with final conditions equal to zero, where φ ε is a smooth function such that φ ε (t) = 1 for t ∈ (0, T − ε) and φ ε (t) = 0 for t ≥ T , for a given small ε > 0. This function φ ε facilitates the matching of the boundary condition with the conditions at time t = T by introducing a smooth transition. Y. Hristova, in his work [17], estimates the error in approximating the initial data ϕ(x) in problem (2) at finite time when ψ(x) = 0, and measuring v(x, t) on the boundary of Ω. We follow the ideas of Hristova to estimate the error committed in our case, when the first initial condition ϕ(x) is zero. Approximating the second initial condition ψ(x) by measuring ∂ t v(x, t) on the boundary, considering a finite time T > 0, a non-trapping c(x), and a bounded domain Ω, we obtain a similar estimate (for a proof of this result, see Subsection 3.2).
. Also, let v, w, and w be solutions of the problems (2), (3), and (4), respectively, with ϕ(x) ≡ 0. We assume the support of ψ is contained in Ω, with Ω ⊂ R n a bounded set. We also assume c(x) > c 1 > 0, (1 − c(x)) has compact support in Ω, and c(x) satisfies the nontrapping condition. Then, we have the following estimates: where η k (t) = t 1−n−k for even n (k = 3 for t ∈ (0, 1), and k = 1 for t > 1), η k (t) = e −δt for odd n, and the constant C depends on Ω and c(x).

2.2.
Example of an exact reconstruction. We now consider the following example in R 3 given by: We want to know the initial condition ψ(x) of the next problem with constant propagation velocity c(x) ≡ 1. Then, we measure m(x, t) = ∂ t u(x, t)| ∂Ω×[0,T ] over a bounded domain Ω at a time T long enough such that we register all the signal on ∂Ω. In this case, the use of φ ε is unnecessary, since all waves leave Ω after certain time T 0 < T . Hence, in this case, we send back the measurement in the homogeneous media by considering (4) as follows Thus, if we solve (6) at time T , we obtain w(x, T ) = ψ(x) for all x ∈ Ω. For this case, we consider the initial conditions (at time t = 0) instead of the conditions at final time t = T . To make this change consistent, we reverse the measurements. For a simple proof based on Huygens' principle of how time-reversal mirror works in this exact case example, see [4]. For further information about another approaches on time-reversal methods, we refer to [4,31,2]. 3. Source time reversal (STR) method. In this section, we extend the timereversal idea to a non homogeneous wave equation, as that considered in our model problem. In this extension, we introduce a new approach of a time-reversal method called Source Time Reversal (STR). This method employs the information of the temporal source dependence to improve the reconstruction of the spatial source term. We also provide an error estimate.
3.1. State of the problem. Based on the Duhamel's principle idea [11], we consider the auxiliary problem (FORWARD PR.) inside Ω, where the unknown term of the source f (x) of (ORIGINAL PR.) now appears as part of the initial conditions in Ω.
The solution v(x, t) of problem (FORWARD PR.) is related with the solution u(x, t) of the problem (ORIGINAL PR.) inside Ω by the convolution here v x (·) represents v(x, ·). This convolution allows to relate a non homogeneous problem with a homogeneous initial-value problem. Also, it gives us a relation between the boundary information in problem (FORWARD PR.) and the measurements obtained in (ORIGINAL PR.). Thus we have the following result. The proof of this Proposition is direct from the definition of u and v.

Proposition 1.
Let v be solution of problem (FORWARD PR.) and u be as in (6). Then, u is solution of problem (ORIGINAL PR.).
To perform the time-reversal mirror method on (FORWARD PR.), it is necessary to know the measurement function m v (x, t) on the boundary, this means, to known the values of ∂ t v(x, t) on ∂Ω × (0, T ). Hence, from (6), the problem is reduced to the following Volterra equation of first kind Thus, for a given x ∈ ∂Ω, we want to find ∂ t v(x, t) by knowing ∂ t u(x, t) and g(t) for all t ∈ (0, T ). From Theorem 2.2 and assuming that the v-measurements are exact, the error estimate in the reconstruction function f (x) is bounded by where η(t) = t 1−n−k for even n (k = 3 for t ∈ (0, 1), and k = 1 for t > 1) and η(t) = e −δt for odd n.
On the other hand, to solve the above Volterra equation numerically, we use the Fourier transform to obtain Then, we can write ∂ t v(x, t) in terms of the following inverse Fourier transform: for all ω such that F(g)(ω) = 0, where F −1 represents the inverse Fourier transform. Since all functions in the convolution are defined for t ≥ 0, the known property F(∂ t u x * g) = F(∂ t v x )F(g), is still valid for this type of convolution (the proof is straightforward from the definition of the convolution).
Hence, we define the operator A : L 2 (0, T ; H 1/2 (∂Ω)) × W 1,∞ (0, T ) → L 2 (0, T ; H 1/2 (∂Ω)) as follows (7) A In the above formula, we have introduced a small regularization constant c 0 > 0 to avoid a possible division by zero. We estimate the velocity measurements on the boundary in (FORWARD PR.) with the velocities measured by the geophones, by using the operator Once we obtain the measurements m v (x, t), we focus on problem (BACKWARD PR.), where the boundary condition of the equation is provided by m v (x, t).
Finally, solving the problem (BACKWARD PR.) at time t = 0, we obtain f (x) := w(x, 0) for all x ∈ Ω, where f (x) is an approximation of f (x). We summarize the above methodology in Figure 1.
Proof sketch of Theorem 2.2. We follow the proof performed in [17]. To do this, we define the error function e T (x, t) = w(x, t)− w(x, t), where w(x, t) and w(x, t) are solutions of (3) and (4), respectively. We also select φ ε (t) such that max 0≤t≤T |φ ε (t)| ≤ C 1 /ε and max 0≤t≤T |φ ε (t)| ≤ C 2 /ε 2 , where C 1 and C 2 are constants. Then, the error function is solution of the following problem: Considering the operator of harmonic extension E : H 1/2 (∂Ω) → H 1 (Ω), such that for γ ∈ H 1/2 (∂Ω), we have a unique θ := Eγ solution of the problem Then, it is possible to obtain the following inequality: We can estimate the terms on the right hand side in the above identity by using the extension operator, the trace operator (see e.g. [11]), and the bounds for the derivatives of the function φ ε (t): . Thus, we have the following estimate, As we can see in the right hand side, the maximum is considered over t moving between T − ε and T . Then, we can use Theorem 3.1 to estimate the above norms in terms of the functions η k ϕ, and ψ Let us notice that η k are decreasing functions and η k (t) ≤ η k+1 (t) for t ∈ [0, 1] and η k+1 (t) ≤ η k (t) for t ∈ [1, ∞). Then, where k = 3 for T ∈ [0, 1] and k = 1 for T > 1. C is a constant depending on Ω and c(x).

Numerical results.
In this section, we present four two-dimensional numerical experiments in different scenarios to study the STR method. The first experiment analyzes the influence of the smoothness of f (x) and g(t) on the reconstruction f (x). The second experiment studies the reconstruction behavior of f (x) when the temporal source term g(t) is not accurately estimated, which could be produced due to an incorrect estimation of the origin time t 0 or an improper calculation of the shape of g(t). The third experiment shows the influence of constant c 0 in the reconstruction of f (x). The last experiment considers a synthetic seismic event and shows the differences between the reconstruction of f (x) with the classical timereversal mirror method and with our STR method.

Numerical implementation.
For these experiments, we use an explicit finite difference scheme of centered differences. We code the experiments in MAT-LAB considering the Courtant-Friedrichs-Lewy (CFL) condition [28] to obtain the convergence of the finite difference scheme. Also, it is important to consider an appropriate regularization constant for the STR method stability. In this case, and for all the experiments, we consider the regularization constant c 0 = 0.01. Finally, for implementing the Fourier transform and its inverse, we use the standard functions of MATLAB.
In the first, second, and third experiments, we consider the following parameters: a domain Ω = (−3 m, 3 m) × (−3 m, 3 m), a total time T = 23 s, and a propagation velocity c(x) ≡ 1 m/s. For the numerical discretization, we consider the step size ∆x = 0.1 m, ∆y = 0.1 m, and ∆t = 0.025 s. In addition, for these experiments we assume the geophones are located over the entire boundary of Ω. For the fourth experiment, we utilize a domain Ω = (−300 m, 300 m) × (0 m, 600 m), where the first component represents the distance along the surface and the second one represents the depth in the Earth. In this experiment, we consider the parameters: T = 0.5 s, and c(x) ≡ 2500 m/s. For the numerical discretization, we consider ∆x = 5 m, ∆z = 5 m, ∆t = 1.4 × 10 −3 s. The geophones are located on the surface and we add noise to the measurements. For the noise, we compute the standard deviation of all information recorded on the boundary, and we add the standard deviation to the clean measurements weighted by a factor as follow m k i noise = m k i clean + f actor * std(m) * rand k i , where rand k i is a uniformly distributed random number in the interval (−1, 1) for each gephone i and for all k in the time discretization.

4.2.
Influence of the smoothness of f (x) and g(t) on the reconstruction. This experiment consists on studying the influence of the smoothness of functions f and g on the reconstruction process. For this purpose, we consider functions f (x) and g(t) with different degrees of smoothness, namely, C ∞ , C 0 , and discontinuous, as we can see in Figures 2 and 3. Figure 2. Functions selected as temporal source terms g(t). Figure 3. Functions selected as spatial source terms f (x).
In this experiment, we simulate "tremors" with synthetic sources f i (x)g j (t), where i, j ∈ {1, 2, 3}, and we use our STR methodology measuring the displacement velocity signal over the entire boundary to reconstruct the spatial source term f i (x).
In Figure 4, we show the results of these reconstructions. The first row (Figures 4a, 4b, and 4c) shows the reconstruction results of the three different spatial source terms f 1 (x), f 2 (x), and f 3 (x) under the assumption of tremors generated by sources f i (x)g 1 (t), i ∈ {1, 2, 3}, respectively. The second row (Figures 4d, 4e, and 4f) describes the reconstruction results of the three different spatial source terms under the assumption of tremors generated by sources f i (x)g 2 (t), i ∈ {1, 2, 3}. Similarly, the last row (Figures 4g, 4h, and 4i) shows the result of the spatial sources reconstruction under tremors generated by sources of type f i (x)g 3 (t), i ∈ {1, 2, 3}.
In addition, we present the relative error in L 2 -norm of these experiments in Table  1. From the results of this experiment, we can see that the best results in terms of the smallest relative error percent are obtained for a discontinuous (piecewise constant) g function. Also, we see that a change in the spatial component produce more variations in the error than a change in the temporal component. This could be due to the stability constant needed to estimate the boundary measurements m v (x, t) by the operator in (7), since without this constant the method became unstable with respect to small variations of g(t). Furthermore, we see from this experiment that the discontinuous spatial source f 3 (x) is the one that presents worst reconstruction with errors ranging from 4.1% to 8.7%.

4.3.
Sensitivity of the reconstruction with respect to g(t). In this subsection, we study the influence of employing a g(t) different than the one used for generating Figure 4. Spatial source term reconstruction for the different sources f i (x)g j (t) i, j ∈ {1, 2, 3}. Table 1. Summary of the relative error in experiment smoothness of f (x) and g(t).
the tremor when reconstructing the source term f (x). For this purpose, we simulate synthetic "tremors" with the same three spatial source terms f 1 (x), f 2 (x), and f 3 (x) used in the previous experiment (see Figure 3), and two temporal source terms g a (t) and g b (t) (see Figure 5). Then, we simulate six synthetic tremors using the source combinations f 1 (x)g a (t), f 2 (x)g a (t), f 3 (x)g a (t), f 1 (x)g b (t), f 2 (x)g b (t), and In order to test the robustness of our STR method, we select temporal source terms g(t) different to the one used for generating the tremors, g a (t) and g b (t). We in experiment sensitivity with respect to g(t).
select the following function g γ to reconstruct f i (x), for i ∈ {1, 2, 3} on each case: Figure 5. Functions selected as temporal source terms g(t) to generate tremors.
For this experiment, we reconstruct the spatial source term using the following set of values for γ: {0.6, 0.7, 0.8, 0.9, 1.0}. Figures 6 and 7 show the results of the reconstruction in this experiment, and Table 2 summarizes the relative errors obtained when reconstructing the different tremors using g γ (t) for the different values of γ. In this table, the first row indicates the source employed to generate the tremor, and the first column corresponds to the γ used in g γ (t) to reconstruct f i (x).
As we can see in Table 2, the best results are obtained for γ = 0.8. This could occur because function g γ=0.8 is the closest one to g a and g b in terms of the area contained under the curve.

4.4.
Influence of the reconstruction with respect to c 0 . We analyze the influence of constant c 0 in the reconstruction of the spatial source term. To do this, we select the three spatial source terms f 1 (x), f 2 (x), and f 3 (x) used in the previous experiments, and we add the MATLAB's phantom Modified Shepp-Logan as a fourth spatial source f 4 (x) (see Figure 9a). We also consider the three temporal source terms g 1 (t), g 2 (t), and g 3 (t) introduced in Subsection 4.2.
In the following, we study the relative error in the reconstruction of f i (x) for waves generated by the sources f i (x)g j (t) for i ∈ {1, 2, 3, 4} and j ∈ {1, 2, 3}, for different values of c 0 .
(a) c0 analyis for f1(x) and f2(x) (b) c0 analyis for f3(x) and f4(x) Figure 8. Relative error variation of the reconstruction with respect to the constant c 0 . Figure 8 shows the relative errors (in the L 2 -norm) as a function of constant c 0 for various source combinations. For each source, we distinguish two regimes as we increase c 0 . In the first part (small values of c 0 ), the reconstruction error diminishes as we increase c 0 . Then, after attaining certain optimal value of c 0 , the opposite behavior occurs. As we observe on the figures, the optimal value of c 0 is source dependent. Finding its optimal value theoretically for any source type is beyond the scope of the paper and will be further investigated in future studies. Figure 9 shows reconstructions of the MATLAB's phantom for different values of c 0 and different temporal sources. Table 3 displays the relative errors corresponding to that figure. The boldface numbers in Table 3 indicate that the best errors are achieved for a "moderate" value of c 0 that is larger than 0 but small enough, as shown in Figure 8 for different sources. For waves generates with the source f 4 (x)g 3 (t), the numerical simulation diverge if we try to reconstruct the spatial source term with c 0 = 0. Therefore, for this case we have selected a small constant c 0 = 10 −5 (Figure 9h), and we have obtained a small relative error (in comparison with the sources f 4 (x)g 1 (t) or f 4 (x)g 2 (t)). Let us recall that g 3 is the discontinuous time-function, and this result is consistent with the observations presented in Subsection 4.2.  Fig. 9h; c 0 = 10 −5 ) ( Fig. 9i; c 0 = 0.01) ( Fig. 9j; c 0 = 0.1) Table 3. Relative errors when reconstructing Phantom's source. proposed STR method. For this purpose, we place a synthetic source 300 m under the surface ground, as illustrated in Figure 10a, and we select the Ricker wavelet as the temporal source term (see Figure 10b). In this experiment, we consider two geophones distributions on the surface to compare both methods. In the first case, the geophones are located every 10 meters, and in the second case, they are located every 50 meters.
To add noise to the geophones measurement, we consider a factor of 0.5 in both cases.    Figure 11. Spatial source term reconstruction in seismic experiments.
Figures 11a and 11b (TRM and STR method respectively) describe the results of the source reconstruction for the first experiment, where the geophones are located every 10 meters; and in Figures 11c and 11d we can see the results for the experiment with the geophones distributed 50 meters apart. To study these results, we consider two aspects: the ε-support and the shape of the reconstruction. We call ε-support to {x ∈ Ω, such that | w(x, 0)| > ε}, and we call shape to { w(x, 0), such that x ∈ Ω}. As the synthetic spatial source and the reconstructions are in different scale, it is necessary to normalize them for comparison purposes. Nonetheless, notice that the figure reconstruction scales are below the original, this is because the geophones are located only on the ground surface, and we lose information in the other directions; however, STR method is able to recover more signal information than the standard TRM method, since we consider the temporal source information on the reconstruction.
From the point of view of ε-support, it is necessary to consider a smaller threshold value to obtain a better reconstruction of the source support. As an example, with a threshold ε = 0.1, we obtain a relative error of 7.5% with STR method and a relative error of 13.5% with the TRM method. This means that we significantly reduce the artifacts appearing in the reconstruction. From the point if view of the shape, we also improve the reconstruction by modifying the measurements with the Duhamenl's principle to obtain more signal and increase the accuracy of the reconstruction. If we normalize the reconstructions, the relative error obtained with the STR method is 2.9% and with the classical TRM method is 4.2% for the first geophones' distribution. Then, it is clear that the STR method reduces the artifacts on the reconstruction image, and also improves the shape accuracy of the source reconstruction.
5. Conclusions. We described a novel source reconstruction method based on a time-reversal method which is able to properly recover the spatial source dependence in microseismic events. The methodology improve the reconstruction of the spatial source term, regardless of the source time duration. This characteristic is the principal difference with classical time-reversal methods, since traditional methods lose precision when the temporal source term is different from a Dirac delta function.
We obtain an error estimate in the wave equation for the time-reversal mirror method in the case when the boundary measurements are the wave displacement velocities over the entire boundary. This estimate depends on the domain, the propagation velocity map, and the initial displacement velocity. In terms of our original problem, this estimate depends on the domain where the measurements are acquired, the propagation velocity map, and the source. The error decays with T exponentially for odd dimensions and polynomially for even dimensions.
Our proposed method properly reconstructs continuous (exhibiting a relative error below 3%) and discontinuous (relative error below 9%) spatial sources when the velocity information is available over the entire boundary. In the case with partial boundary information, the reconstruction also exhibits good results, although the reconstructed scale is significantly smaller than in the original source. Nonetheless, the relative error remains below 3% when the reconstruction is normalized. In addition, we compare the STR method with the classical approach of time-reversal mirror method, and we show quantitatively the advantages of our methodology, reducing the error by approximately a factor of two in the considered examples.
As future work, we plan to study the stability of the Fourier transforms, which in our case is obtained via the use of a regularization constant. Another way to guarantee stability could be studying the Volterra equation and searching for a numerical way to solve it. In another line of research, we shall implement the STR method for linear elasticity in order to obtain a more realistic seismic model.