remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \headersConvexification for the 3D Problem of Travel Time Tomography M. V. Klibanov, J. Li, V. G. Romanov, and Z. Yang
Convexification for the 3D Problem of Travel Time Tomography
Abstract
The travel time tomography problem is a coefficient inverse problem for the eikonal equation. This problem has well known applications in seismic. The eikonal equation is considered here in the circular cylinder, where point sources run along its axis and measurements of travel times are conductes on the whole surface of this cylinder. A new version of the globally convergent convexification numerical method for this problem is developed. Results of numerical studies are presented.
Key Words: eikonal equation, geodesic lines, coefficient inverse problem, travel time tomography in 3d, globally convergent numerical method, Carleman estimate, numerical solution
2020 MSC codes: 35R30.
1 Introduction
The travel time tomography problem (TTTP) has a long history since it has broad applications in seismic studies of Earth, see, e.g. [21, 24]. First, we refer here to the pioneering works of Herglotz in 1905 [7] and Wiechert and Zoeppritz in 1907 [25], who have solved this problem in the 1d case. There are some important uniqueness and stability results for the TTTP. Since we are focused here on the numerical aspect, we cite only a few of those [17, 18, 19, 21]. We also refer to [22, 26] for numerical studies of the TTTP by the methods, which are significantly different from the convexification method of this paper.
Our goal here is to present the theory and computational experiments for a new version of the globally convergent convexification numerical method for the TTTP in the 3d case. The theory of the first version of the convexification for this problem was published in [12], [14, Chapter 11] and numerical studies were published in [15]. However, it was assumed in these references that the domain of interest is a vertically elongated rectangular prism, and point sources run along an interval of a straight line, which is located in a horizontal plane “below” As a result, only horizontally oriented abnormalities were imaged in [15].
Since we do not want to be restricted to only horizontally oriented abnormalities, then we consider in this paper the case when is a vertically oriented circular cylinder with point sources running along its axis and measurements of travel times being conducted on the lateral surface, top and bottom boundaries of this cylinder, see Figure 1. This source/detector configuration has applications in seismic imaging of wells. Naturally, our source/detector configuration requires some changes in the theory of [12, 15].
Below we first formulate and prove a new theorem about the monotonicity of the solution of the eikonal equation with respect to the radius of the cylindrical coordinate system. Next, using this theorem, we present a new version of the convexification method, which is adapted to the sources/detectors configuration of Figure 1. Next, we formulate some theorems of the convergence analysis. Since their close analogs were proven in [15], then we do not prove them here. Finally, we describe our numerical studies and present their results.
Definition. We call a numerical method for a Coefficient Inverse Problem (CIP) globally convergent, if there exists a theorem claiming that this method delivers points in a sufficiently small neighborhood of the true solution of that problem without any advanced knowledge of this neighborhood. The size of this neighborhood depends only on the noise level in the input data for this CIP. In other words, a good first guess about the true solution is not required.
The first theoretical results about the convexification method for CIPs were obtained in [9, 10]. More recent works on this method for various CIPs contain both theory and numerical results, see, e.g. [1, 2, 13, 14, 15, 16] and references cited therein. The development of the convexification is caused by the desire to avoid the phenomenon of local minima and ravines of conventional least squares cost functionals for CIPs, see, e.g. [3, 4, 5]. In the convexification, a globally strongly cost functional is constructed. The key element of this functional is the presence of the Carleman Weight Function (CWF) in it. This is the function, which is used as the weight in the Carleman estimate for an associated operator. Usually this is a partial differential operator. In the case of the TTTP, however, this is a Volterra-like integral operator. Finally, it is proven that the gradient descent method of the minimization of this functional converges globally to the true solution of the corresponding CIP as long as the level of noise in the data tends to zero.
All functions considered below are real valued ones. In section 2 we pose our inverse problem. In section 3 we prove the above mentioned important theorem about the monotonicity of the travel time function with respect to the radius of the cylindrical coordinate system. In section 4 we derive an integral differential equation, which does not contain the unknown coefficient. In section 5 we present an important element of our method: the orthonormal basis of [11], [14, section 6.2.3]. In section 6 we derive a boundary value problem and also prove a Carleman estimate for the Volterra integral operator. In section 7 we introduce partial finite differences for our problem. In section 8 we present the convexification functional and formulate theorems of its global convergence analysis. Section 9 is devoted to numerical studies.
2 Statement of the Problem
Below denotes points in Let be the speed of waves propagation, , where is the refractive index. The function generates the Riemannian metric [21, Chapter 3]
(1) |
Let be a source of waves and be an observation point. The time, which the wave needs to propagate from a source to a point is called the “first arrival time” or “travel time”. In the case of a heterogeneous medium with the first arriving signal, which arrives at the arrival time to the point , travels along the so-called “geodesic line” generated by metric (1). The first arrival time is
where is the element of the euclidean length. The function satisfies the following eikonal equation [21, Chapter 3]:
(2) |
Consider two numbers and let be a small number. We define the domains as:
(3) | |||
(4) | |||
(5) |
Hence, and are vertically oriented cylinders with the radii and respectively, and is the vertical axis of both cylinders, along which sources run.
Consider the cylindrical coordinates Therefore by (3) and (4) in cylindrical coordinates
(6) | |||
(7) | |||
(8) |
We consider below the eikonal equation (2) only in the domain We assume that the refractive index satisfies the following conditions:
(9) | ||||
(10) | ||||
(11) |
In cylindrical coordinates eikonal equation (2) has the form
(12) |
By (5) denote
(13) |
Travel Time Tomography Problem (TTTP). Let conditions (3)-(13) hold. Find the refractive index for assuming that the following functions , are known:
(14) | ||||
(15) | ||||
(16) |
Everywhere below we assume that the following condition is valid:
Assumption of the Regularity of Geodesic Lines. For any pair of points there exists unique geodesic line connecting these two points.
3 The Monotonicity Property of the Function
Theorem 3.1. Assume that conditions (3)-(11) are in place. Then the following inequality is valid:
(17) |
Proof. Denote , , . Then the eikonal equation (12) becomes:
(18) |
Differentiating (18) with respect to and using the equalities
we obtain that along the geodesic line with
(19) |
(20) |
where is the element of the Riemannian length of .
Equations (19), (20) form a coupled system of ordinary differential equations. We need to figure out initial conditions for this system. Consider the geodesic line that connects the point with the point The tangent vector to this line at the point is the vector
, , . Hence, the initial conditions for this geodesic line are:
(21) | ||||
(22) |
Thus, solving Cauchy problem (19)-(22), we obtain equations of in the parametric form
Due to the regularity of geodesic lines, there exists a one-to-one correspondence between points and In particular, this means that for any fixed there exists a one-to-one correspondence between and parameters . We denote this correspondence as:
Moreover, we have along the geodesic line
(23) |
For any point the geodesic line is:
where
(24) |
By (9) is a part of a straight line. Let be the intersection of with the lateral boundary of the cylinder Then for Hence,
(25) |
Hence, by (11), (20) and (25) for
4 Integral Differential Equation
By the first line of (28)
(29) | ||||
(30) |
Substituting (28)-(30) in the eikonal equation (12), we obtain the following equation in the domain defined in the last line of (28):
(31) |
Differentiating (31) with respect to and using we obtain
(32) |
In addition, (15), (16) and (26) imply:
(33) |
(34) |
Thus, we came up with the following boundary value problem (BVP).
Boundary Value Problem 1 (BVP1). Find the function satisfying conditions (32), (33) in the domain (34).
Suppose that we have solved this problem. Then we substitute its solution in the left hand side of equation (31) and find the target function Clearly BVP1 is a very complicated one. Therefore, we construct below an approximate mathematical model for its solution, and then confirm the validity of this model computationally.
5 A Special Orthonormal Basis in
The basis described in this section was first introduced in [11]. Then it was used for the convexification method in some other works, see, e.g. [12, 14, 15]. Consider the set of functions Obviously, this set is complete in and these functions are linearly independent. Apply the Gram-Schmidt orthonormalization procedure to this set. Then we obtain the orthonormal basis in . The function has the form , where is a polynomial of the degree . Denote
It was proven in [11], [14, Theorem 6.2.1] that
(35) |
Let be an integer. Consider the matrix By (35) Hence, the matrix is invertible. We observe that neither classical orthonormal polynomials nor the basis of trigonometric functions do not provide a corresponding invertible matrix . This is because the first function in these two cases is an identical constant, implying that the first column of that analog of is formed only by zeros.
6 Boundary Value Problem 2
The first step of our approximate mathematical model mentioned in section 4 is the assumption that the function can be represented as truncated Fourier-like series,
(36) |
We do not prove convergence of our method as We point out that the absence of such proofs is a rather common place in the theory of Inverse Problems, basically due to their ill-posedness, see, e.g. [6, 8].
Consider the D vector function of coefficients of the series (36)
(37) |
We want to find the vector function in (37). Naturally, we assume that for functions in (14)-(16)
(38) |
(39) |
Denote
(40) |
Substituting (36)-(40) in (32), multiplying the obtained equality sequentially by the functions and integrating with respect to we obtain
(41) |
where is a vector function of its arguments. The fourth line in (41) means the periodicity condition with respect to the angle Thus, (41) is the boundary value problem for the nonlinear system of integral differential equations with respect to the vector function Here, the th component of the D vector function has the form
Here has the following form after the integration by parts with respect to :
(42) | ||||
And has the form after the integration by parts with respect to :
(43) | ||||
In (42), (43) the function has the form (36), and the vector function has the form (37). Thus, we have obtained Boundary Value Problem 2 (BVP2).
Boundary Value Problem 2 (BVP2). Let be the domain defined in (8). Find the D vector function satisfying conditions (41), where the function is connected with via (36), (37).
We now formulate a Carleman estimate for the Volterra integral operator occurring in (41)-(43). First, we introduce the Carleman Weight Function for this Volterra operator, where is a parameter. This function is:
(44) |
Theorem 6.1 (Carleman estimate for the Volterra integral operator). The following Carleman estimate is valid:
(45) |
7 Partial Finite Differences
To solve boundary value problem (41), we use partial finite differences with respect to the variables and without imposing an additional assumption. This is the second and the last step of our approximate mathematical model. More precisely, we assume that problem (41) in finite differences with respect to and where the grid step sizes and with respect to and satisfy the following inequality:
(47) |
where the number is fixed. Consider finite difference grids with respect to and
(48) | ||||
(49) |
Denote The vector function is turned now in the semi-discrete function, which is defined on the grid (48), (49) for all i.e. remains a continuous variable. We denote this function . The semi-discrete form of (37) is:
(50) |
And the semi-discrete form of (36) is:
(51) |
7.1 Some specifics of partial finite differences
We use the second order approximation accuracy for derivatives
(52) | ||||
(53) |
A separate question is about ensuring boundary conditions in the second and third lines of (41) as well as the periodicity condition in the fourth line of (41). To satisfy boundary conditions in the second and third lines of (41), we use in (52):
(54) | ||||
(55) |
In addition, to ensure the periodicity condition in the fourth line of (41), we set
(56) |
Let be a Hilbert space. Consider the direct product of such spaces times. If is the norm of in , then the norm of the vector function is set below as:
We now introduce some Hilbert spaces of semi-discrete functions, similarly with [15]. Denote Along with the above notation , , we also use the equivalent one, which is For the domain defined in (8) we denote
Next, using (48) and (49), we define those Hilbert spaces as:
(57) |
(58) |
In (57) derivatives are understood as in (52)-(55) being supplied by the periodicity condition (56).
7.2 Completion of the approximate mathematical model
It was stated in the end of section 4 that we construct an approximate mathematical model for our travel time tomography problem, and then we verify this model computationally. The first step of this model is the truncated Fourier-like series (36). We remind that truncated series are used quite often in developments of numerical methods for various inverse problems. On the other hand, convergencies of such methods as the number of terms tends to infinity are usually not proven, see, e.g. [6, 8]. This is caused by the ill-posed nature of inverse problems. The second and the final step of that approximate mathematical model is the assumption of partial finite differences with condition (47). This step is reflected in (48)-(56).
8 Convexification
8.1 The cost functional of the convexification
Let be an arbitrary number. Let be the number in (17). Define the set of vector functions of (50) as
(59) |
where the space is defined in (57).
Thus, we have replaced BVP2 with BVP3:
Boundary Value Problem 3 (BVP3). Find the vector function
To solve BVP3, we consider the following minimization problem:
Minimization Problem. Minimize the following functional on the set
(60) |
8.2 The global convergence
We now formulate theorems of the global convergence analysis for the Minimization Problem. Both their formulations and proofs are similar with those of [15]. Therefore, we omit proofs. We now briefly explain the role of the Carleman Weight Function in the proof of our central result, which is Theorem 8.1. Given the set in (59), this theorem claims the strong convexity of the functional on for sufficiently large values of the parameter First, consider for a moment only the quadratic functional Since the matrix is invertible, then this functional is strongly convex on the entire space for any value of However, the presence of the nonlinear term in (60) might destroy the strong convexity property. To dominate this term, we use in the proof of Theorem 8.1 the Carleman estimate of Theorem 6.1, which, roughly speaking, states that the nonlinear term is dominated by the strongly convex quadratic term
Another question is on what exactly “sufficiently large is. Philosophically this is very similar with any asymptotic theory. Indeed, such a theory usually claims that as soon as a certain parameter is sufficiently large, a certain formula is valid with a good accuracy. In a computational practice, however, one always works with specific ranges of many parameters. Therefore, only results of numerical experiments can establish which exactly values of are sufficient to ensure a good accuracy of . In particular, we demonstrate below that is an optimal value of the parameter for our computations. We also note that in all previous publications on the convexification of this research group optimal values of were , see, e.g. [13]-[16].
Theorem 8.1 (strong convexity). Assume that conditions (6)-(11) hold, and let be the number in (17) and (59). Then:
1. The functional has the Fréchet derivative at every point and for all Furthermore, the Fréchet derivative satisfies Lipschitz continuity condition on i.e. there exists a number
depending only on listed parameters such that the following estimate holds:
2. There exist a sufficiently large number and a number
(61) |
both numbers depending only on listed parameters, such that for every the functional is strongly convex on the set The strong convexity means that the following estimate holds
3. Furthermore, for every there exists unique minimizer of the functional on the set and the following inequality holds:
Below denotes different positive numbers depending only on parameters listed in (61). In the regularization theory [23], the minimizer of functional (60) is called “regularized solution”. It is important to estimate the accuracy of the regularized solution depending on the level of the noise in the data. To do this, we recall first that, following the regularization theory, we need to assume the existence of the “ideal” solution of BVP3, i.e. solution with the noiseless data. The ideal solution is also called “exact” solution. We denote this solution We denote the noiseless data in the fourth line of (59) as We assume that the exact solution where is the following analog of the set in (59)
(62) |
Let be the level of the noise in the boundary data . To simplify the presentation, we assume that the noise is not introduced in the function in (14), i.e. we assume that although this case can also be included. Suppose that there exist such extensions and of the pairs of boundary data and that for and
(63) |
Theorem 8.2 (an estimate of the accuracy of the regularized solution). Assume that conditions of Theorem 8.1 hold. Let the exact solution where the set is defined in (62). Assume that conditions (63) hold. Furthermore, assume that where the number is so small that
(64) |
Let be the number of Theorem 8.1 in the case when the number in (61) is replaced with
(65) |
For any let be the minimizer on of the functional on the set which is found in Theorem 8.1. Then the following accuracy estimate holds:
For we now construct the gradient descent method of the minimization of the functional Let be a number and let be an arbitrary point. The gradient descent method is constructed as the following sequence:
(66) |
Since by Theorem 8.1 , then all vector functions satisfy the same boundary condition as the ones in the third line of (59), see (58).
Theorem 8.3. Let in (64) and let the number Suppose that the exact solution . Let where is defined in (65). Then there exists a sufficiently small number such that for any the sequence ( 66) In addition, there exists a number such that the following convergence estimate holds
Furthermore, let the unction be the exact semi-discrete target function and let be the semi-discrete target function, which is found via the substitution of , sequentially, first in the semi-discrete analog of the first line of (37), then in (36) and finally in the semi-discrete analog of the left hand side of (31). Then the following convergence rate holds:
9 Numerical Studies
In this section we describe our numerical studies. We specify parameters of the domains and in (3) and (4) as:
The refractive index in eikonal equation (2) is taken as
(67) |
To ensure that as in (10), we smooth out near the boundaries of our tested inclusions. Then we set:
(68) | |||
(69) |
In the numerical tests below, we take which means 1.5:1, 3:1 and 5:1 inclusion background contrasts respectively, see (68). To demonstrate that our numerical method can work with inclusions, which have sophisticated non-convex shapes with voids in them, we take shapes inclusions like letters ‘’ and ‘’, This is similar with our previous works [13]-[16] on the convexification.
After choosing the function , we use the fast marching toolbox ”Toolbox Fast Marching” [20] in MATLAB to solve eikonal equation (2). Then we convert data from Cartesian coordinates to polar coordinates to gain the observation data in (14)-(16). Then we solve the Minimization Problem (60) to gain the computed solution . Finally, we convert the reconstructed solution from polar coordinates to Cartesian coordinates to exhibit results.
To solve eikonal equation (2) for generating the observation data in (14)-(16), we choose in (48) and (49), as well as we choose to generate discrete points along the direction. Because rather than in (8), then the first interval along direction is . To solve the Minimization Problem (60), we choose .
To guarantee that the solution of the Minimization Problem (60) satisfies the boundary conditions in (54)-(56), we adopt the Matlab’s built-in optimization toolbox fmincon. The iterations of fmincon stop when we get
(70) |
The starting point of iterations of fmincon is chosen as . Although the starting point does not satisfy the boundary conditions in (54)-(56), fmincon makes sure that the boundary conditions (54)-(56) are satisfied on all other iterations of fmincon.
We consider the random noise in observation data in (14)-(16) as follows:
(71) | ||||
(72) | ||||
(73) |
In (71)-(73), is the Gaussian random variable depending on variables and are the uniformly distributed random variables in the interval depending on variables . Also, and , which correspond to the and random noise levels respectively. Recall that, to simplify the presentation, we did not consider the noise in the function in our theoretical derivations, see (63). Nevertheless, we still introduce random noise in this function in our numerical studies, see (71).
To calculate the derivative and the derivative of the noisy function , , in (29)-(30), as well as the derivative of noisy functions , in (33), we firstly use the natural cubic splines to approximate the noisy data in (71)-(73). Then we use the derivatives of those splines to approximate the derivatives of corresponding noisy data.
The parameters and are the two key parameters in our numerical method. We firstly choose the optimal , when is large enough with . Then, keeping that optimal we find the optimal value of . This is done in Test 1. An important point to make here that, once chosen in Test 1, that optimal pair is kept for all other tests.
Test 1. We test the case when the inclusion in (67) has the shape of the horizontally oriented letter ‘’ with in it. The goal of this test is to find the optimal values of the parameters and . To make sure that the chosen parameter does not impact our results with different values of , we set to be large enough, i.e. . Computational results are displayed in Figure 2. We observe that the images have a low quality for . Then the quality is improved with , and the reconstruction quality significantly deteriorates at . Hence, we select as the optimal one.
Next, we fix the optimal value of and test the influence of the parameter . The results with are shown in Figure 3. We can find that the reconstruction for is of a low quality, and the reconstructions for are almost same. Hence, to reduce the computational cost we select as the optimal value of , instead of . Furthermore, we also consider the truncated expansion of the function in (26) with , which is denoted as . We have obtained that
(74) |
It is clear from (74) that is a quite informative case.
Conclusion: We choose
(75) |
as the optimal values of these parameters, and we use these values in Tests 2-5.
Test 2. We test the case when the inclusion in (67) has the shape of the vertically oriented letter ‘’ with in it. Results are presented on Figure 4. An accurate reconstruction can be observed.
Test 3. We test the cases when the inclusion in (67) has the shape of the vertically oriented letter ‘’ with and in it. Results are displayed on Figures 5 and 6 respectively. The inclusion/background contrasts in (68) are respectively and . We see that shape of the inclusion is imaged accurately in both cases. In addition, the computed inclusion/background contrasts in (69) are accurate.
Test 4. We test the case when the inclusion in (67) has the shape of the letter ‘’ elongated along the axis and with in it. Results are presented on Figure 7. We again observe an accurate reconstruction of both the shape of the inclusion and the inclusion/background contrast.
Test 5. We consider the case when the random noise is present in the data in (71)-(73) with and , i.e. with 1% and 3% noise level respectively. We test the reconstruction for the case when the inclusion in (67) has the shape of the vertically oriented letter ‘’ with in it. Results are displayed on Figure 8. The reconstructions of the shape of the inclusion as well as computed inclusion/background contrasts in (69) are accurate.
Remarks 9.1:
-
1.
Recall that by Theorem 3.1 condition (11) is a sufficient condition for our method to work. Recall also that in the data generation process we smooth out tested inclusions in small neighborhoods of their boundaries. Given these, a careful analysis of correct images of Figures 2-8 indicates that condition (11) is satisfied at least in a major part of the domain in each of the above Tests 1-5.
- 2.
Acknowledgments. The work of Li was partially supported by Guangdong Basic and Applied Basic Research Foundation 2023B1515250005. The work of Romanov was partially supported by the Mathematical Center in Akademgorodok under Agreement 075152022281 with the Ministry of Science and Higher Education of the Russian Federation. The work of Yang was partially supported by Supercomputing Center of Lanzhou University.
References
- [1] L. Baudouin, M. de Buhan, and S. Ervedoza, Convergent algorithm based on Carleman estimates for the recovery of a potential in the wave equation, SIAM J. Numer. Anal., 55 (2017), pp. 1578–1613.
- [2] L. Baudouin, M. de Buhan, S. Ervedoza, and A. Osses, Carleman-based reconstruction algorithm for the waves, SIAM J. Numer. Anal., 59 (2021), pp. 998–1039.
- [3] L. Beilina and E. Lindström, An adaptive finite element/finite difference domain decomposition method for applications in microwave imaging, Electronics, 11 (2022), p. 1359.
- [4] A. V. Goncharsky, S. Y. Romanov, and S. Y. Seryozhnikov, On mathematical problems of two-coefficient inverse problems of ultrasonic tomography, Inverse Probl., 40 (2024), p. 045026.
- [5] M. J. Grote and U. Nahum, Adaptive eigenspace for multi-parameter inverse scattering problems, Computers & Mathematics with Applications, 77 (2019), pp. 3264–3280.
- [6] J. P. Guillement and R. G. Novikov, Inversion of weighted Radon transforms via finite Fourier series weight approximation, Inverse Probl. Sci. En., 22 (2013), pp. 787–802.
- [7] G. Herglotz, Űber die Elastizitaet der Erde bei Beruecksichtigung ihrer variablen Dichte, Zeitschr. fur Math. Phys., 52 (1905), pp. 275–299.
- [8] S. I. Kabanikhin, K. K. Sabelfeld, N. S. Novikov, and M. A. Shishlenin, Numerical solution of the multidimensional Gelfand-Levitan equation, J. Inverse Ill-Posed Probl., 23 (2015), pp. 439–450.
- [9] M. V. Klibanov and O. V. Ioussoupova, Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem, SIAM J. Math. Anal, 26 (1995), pp. 147–179.
- [10] M. V. Klibanov, Global convexity in a three-dimensional inverse acoustic problem, SIAM J. Math. Anal., 28 (1997), pp. 1371–1388.
- [11] M. V. Klibanov, Convexification of restricted Dirichlet to Neumann map, J. Inverse Ill-Posed Probl., 25 (2017), pp. 669–685.
- [12] M. V. Klibanov, Travel time tomography with formally determined incomplete data in 3D, Inverse Problems and Imaging, 13 (2019), pp. 1367–1393.
- [13] M. V. Klibanov, J. Li, and W. Zhang, Convexification for an inverse parabolic problem, Inverse Probl., 36 (2020), p. 085008.
- [14] M. V. Klibanov and J. Li, Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data, De Gruyter, Berlin, 2021.
- [15] M. V. Klibanov, J. Li, and W. Zhang, Numerical solution of the 3-D travel time tomography problem, J. Computational Physics, 476 (2023), p. 111910.
- [16] M. V. Klibanov, J. Li, and Z. Yang, Convexification numerical method for the retrospective problem of mean field games, Applied Mathematics and Optimization, (2024), 90:6.
- [17] R. G. Mukhometov, The reconstruction problem of a two-dimensional Riemannian metric and integral geometry, Soviet Math. Dokl., 18 (1977), pp. 32–35.
- [18] R. G. Mukhometov and V. G. Romanov, On the problem of determining an isotropic Riemannian metric in the -dimensional space, Dokl. Acad. Sci. USSR, 19 (1978), pp. 1330–1333.
- [19] L. Pestov and G. Uhlmann, Two dimensional simple Riemannian manifolds are boundary distance rigid, Annals of Mathematics, 161 (2005), pp. 1093–1110.
- [20] G. Peyre, Toolbox fast marching, MATLAB Central File Exchange, (2023).
- [21] V. G. Romanov, Inverse Problems of Mathematical Physics, VNU Press, Utrecht, The Netherlands, 1986.
- [22] U. Schrőder and T. Schuster, An iterative method to reconstruct the refractive index of a medium from time-off-light measurements, Inverse Problems, 32 (2016), p. 085009.
- [23] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola, Numerical methods for the solution of Ill-posed problems, Kluwer, London, 1995.
- [24] L. Volgyesi and M. Moser, The inner structure of the Earth, Periodica Polytechnica Chemical Engineering, 26 (1982), pp. 155–204.
- [25] E. Wiechert and K. Zoeppritz, Uber Erdbebenwellen, Nachr. Koenigl. Geselschaft Wiss. Gottingen, (1907), pp. 415–549.
- [26] H. Zhao and Y. Zhong, A hybrid adaptive phase space method for reflection traveltime tomography, SIAM J. Imaging Sciences, 12 (2019), pp. 28–53.