Recovering Density and Speed of Sound Coefﬁcients in the 2D Hyperbolic System of Acoustic Equations of the First Order by a Finite Number of Observations

: We consider the coefﬁcient inverse problem for the ﬁrst-order hyperbolic system, which describes the propagation of the 2D acoustic waves in a heterogeneous medium. We recover both the denstity of the medium and the speed of sound by using a ﬁnite number of data measurements. We use the second-order MUSCL-Hancock scheme to solve the direct and adjoint problems, and apply optimization scheme to the coefﬁcient inverse problem. The obtained functional is minimized by using the gradient-based approach. We consider different variations of the method in order to obtain the better accuracy and stability of the appoach and present the results of numerical experiments.


Introduction
In this article, we deal with the numerical solution of the coefficient inverse problem, which corresponds to the problems of ultrasound tomography. The problems of developing methods and algorithms to use the ultrasound to recognize the genesis of breast cancer in its early stages have been studied extensively lately [1][2][3][4]. On the mathematical level, such problems are usually considered as inverse problems, when one has to recover the parameters of the model (that in our case describes the propagation of the ultrasound through the object of investigation) by using some measurement data [5][6][7][8]. The inverse problems are known for their ill-posedness and a requirement for a large number of computational resources for the numerical solution. Hence, the goal is to use a numerical algorithm that utilizes the provided data in an effective manner.
The mathematical models for ultrasound acoustics usually have the form of either the second-order equation or the first-order system of PDE equations. The models, based on the second-order wave equation are usually easier to study, and therefore, there are more ways to efficiently solve the direct problem. The first-order system of acoustics, that we consider in this paper, requires more computational resources for solving. Its advantage relies on its close connection to the physics of wave propagation, since the equations can be derived straight from the conservation laws. We mention several papers [9][10][11][12][13][14], where authors investigated inverse problems for a system of hyperbolic partial differential equations. One can also mention the following papers [15,16], where authors considered the direct problems for linear hyperbolic systems and studied their well-posedness.
Since we use the model, that based on the two-dimensional system of acoustic equations, its parameters, that we aim to recover, are the density of the medium and the speed of waves propagation-both are the functions of the two spatial variables. The reconstruction of two coefficients in the partial differential equations with a finite number of observations is very challenging.
The inverse problem of determining several coefficients in scalar hyperbolic equation by data given on the part of the boundary was investigated in [17]. The Carleman estimate was obtained to prove the uniqueness and a Lipschitz stability estimate for the coefficient inverse problem.
In [18] authors considered the reconstruction of three coefficients depended on space variables in the dynamical isotropic system of elasticity from two boundary measurements and proved the uniqueness and a Hölder stability using Carleman estimates in Sobolev spaces of negative order.
The inverse problem of recovering several coefficients in Maxwell's equations was investigated in [19][20][21] by a finite number of measurements. For the coefficient inverse problem, the Lipschitz stability estimate was proven by using the Carleman estimate.
The coefficient inverse problem for the dynamic Maxwell equations was considered [22]. Hölder stability and global Carleman estimate for the inverse problem were proved.
The coefficient inverse problem of the recovering of the magnetic permeability and dielectric permittivity of Maxwell's system in three dimensions by data of the electric field given on the part of the boundary was considered in [23]. The authors applied the Carleman estimates to get the theoretical stability. The inverse problem was reformulated as an optimization problem. The hybrid finite element and difference method was implemented for solving the direct and adjoint problem.
An inverse problem of finding two coefficients in a hyperbolic acoustic equation of the second-order by interior data was considered in [24]. The authors applied a Carleman technique estimates to obtain the Lipschitz stability estimates and therefore unique reconstruction of both coefficients was guaranteed. Numerical experiments of recovering both coefficients by data with noise were presented.
The continuation and coefficient inverse problem of recovering dielectric permeability and conductivity in application to ground penetrating radar was investigated in [25]. The inverse problems were reformulated as optimization problems. To minimize the cost functional gradient method was applied.
The numerical algorithms, based on the S. K. Godunov scheme [26,27] is applied for solving direct problems. Such kinds of methods allow us to construct the effective numerical realization of the physical process and benefits from the usage of the piecewisesmooth structure of the state variables on each time step. If one uses the methods, based on finite approximation, for solving the forward problem in the case of piecewise-smooth medium, then one has to add special conditions on the interface. Such conditions are very hard to add in the case of complex media interfaces and impossible to realize when solving the inverse problem, when one does not have the knowledge of media interfaces.
For the solution of the coefficient inverse problems gradient methods [25,[28][29][30][31] and global-convergence [6,[32][33][34][35] are applied. We should also mention the family of Newtontype methods. However, their drawback is the solution of the additional linear inverse problem, which has to be solved on each iteration. When considering the multidimensional problems, this necessity to deal with that additional linear problem tends to become too complicated.
An optimization approach to a three-dimensional acoustic inverse problem was considered in the time-domain [28]. The velocity and the density were reconstructed by minimizing an objective functional. The gradient of the cost functional was found with an explicit expression via the solution of adjoint problem. The parameters were reconstructed by the conjugate gradient method. The uniqueness of the solution was proved.
The problem of modeling the acoustic radiation pattern of source was considered in [36]. The problem is formulated in the form of control problem for the 2D first-order system of hyperbolic equations. The modelling of the acoustic radiation patterns of sources could allow to improve the resolution of acoustic tomography.

Problem Formulation
From the conservation laws of impulse in direction x and y, and the conservation law of mass [12] let us consider the direct problem of acoustic wave propagation through the 2D medium in the domain Ω = (x, y) ∈ [0, L] × [0, L]: ∂p ∂t + ρc 2 ∂u ∂x u, v, p| t=0 = 0.
Here u = u(x, y, t) is the velocity vector with respect to x, v = v(x, y, t) is the velocity vector with respect to y, p = p(x, y, t) is the exceeded pressure, ρ = ρ(x, y) is the density of the medium, c = c(x, y) is the wave speed. Such system is often used to describe the propagation of the ultrasound through the fluid medium, and the acoustic parameters of the models, that were considered during the numerical experiments are close to fluid. θ Ω (x, y) is the characteristic function of the source location, I(t) has the following form: Here, ν 0 is a frequency.

Methods for Solving the Direct Problem
In this subsection we consider the brief description of the two methods, that we used for solving the direct problem during the experiments. Both of them originates from the upwind scheme, developed by S.K. Godunov [26,27]. His approach was based on using the integral form of the problem, piecewise-constant approximation of the state variables inside the numerical cell and the solution of the Riemann problem-the initial problem with conditions represented by two constant states separated by a discontinuity. The MUSCL-Hancock scheme, first published by van Leer [37], extends the ideas behind the Godunov scheme in order to obtain higher accuracy of the method and nowadays is one of the common approaches for dealing with the hyperbolic systems. One can find more detailed reviews of papers, related to both methods, in [37][38][39].
In order to describe methods, used for solving the direct problem (1)-(4), we use a generalized formulation of the Equations (1) and (2) of the following form:

∂U ∂t
Here U = (u, v, p) is the vector of state variables, and F(U), G(U) are the vectors of fluxes correspondingly. After introducing the discretizing the computational domain into the finite number of cells, one can obtain: The Equation (7) corresponds to the numerical cell (i − 1/2, j − 1/2), where the sub-indexes indicate the state values U on the current time step, and sup-indexes-on the next time step, h x , h y , τ are the grid steps with respect to the spatial coordinates and time correspondingly. The values F, G on the each boundary of the cell considered, are the solutions of the Riemann (or discontinuity decay) problem [27]. For example, the approximation (7) of the first of the Equation (1) has the following form: The value P i,j−1/2 is the solution of the discontinuity decay problem, that arises on the boundaries of the cell considered. The formula has the following form: The two other equations of the system can be considered in the same manner. The right-hand side of the equation can also be easily taken into account. We skip the rest of the formulas, yet one can find them, for example, in [27], as well as the study of the stability of the scheme.
Another approach, that we used for solving the direct problem is MUSCL-Hancock scheme. First, we summarize the basic Godunov approach as the combination of the two following steps:

1.
Obtaining the flux values by solving the Riemann problem; 2.
Updating the state variables on the next time step, using the solution of the Riemann problem.
One should mention, that the Godunov approach based on the piecewise-constant approximation of the parameters of each of the cells. The MUSCL-Hancock scheme uses the piecewise-linear approximation of the parameters, based on the slope limited approximation. Another feature of the scheme is the solution of the Riemann problem on the half-step, which corresponded to 0.5τ. Such update allows the described scheme to second-order accuracy [39]. The workflow of the scheme can now be summarized as follows: 1.
Reconstruction of the state variables on the current time step, using piecewise-linear extrapolation; 2.
Evolution of the reconstructed state variables by conservation laws with a time 0.5τ; 3.
Obtaining the flux values for the "midpoint" time step by solving the Riemann problem for evolved state variables; 4.
Updating the state values from the current time step on the next time step by conservation laws with a time τ.
We skip the formulas for the sake of brevity. However, the precise description of the MUSCL-Hancock scheme can be found in [37][38][39].

Inverse Problems
While solving the direct problem for the system (1)-(4) we suppose, that parameters of the system are known. However, when considering the possible applications, one usually has to solve the problem of calculating these parameters using the additional data-the inverse problem.
In this paper we suppose, that we obtain the data of inverse problem by measuring the pressure inside the receivers: Here we consider the system of N receivers, each located in the corresponded domain Ω k . The inverse problem, therefore, is to recover functions c(x, y), ρ(x, y) in (1)-(4) using the additional information (10).
Inverse problems of acoustics are usually considered in the case of data, given on the part of the boundary. The theoretical study of different methods in the case of spatially distributed receivers is very challenging (even in the case of mathematical models, based on the second-order equation). Thus, we consider the formulated inverse problem from a numerical point of view.
We reformulate inverse problem in the (1)-(4), (10) in operator form Let us reduce the inverse problem (1)-(4), (10) to minimization problem of the following cost functional: We use the gradient-based approach to minimize the functional (12) by considering the following iteration scheme: Here α ∈ (0, ||A|| −2 ) is descent parameter, J (q (n) ) is the gradient of the functional. Let us note [40,41] that Here [A (q)] * is adjoint of Frechet derivative of the operator A.
One can choose the parameter of the descent differently, for example, one can fix α throughout the iterations (Landweber iterations), or search for the best parameter on each step of the optimization process (steepest descent). The initial approximation for the gradient descent is often connected with the prior information about the solution. We mention the structure of initial approximation in the next chapter. Using a priori information about the solution of the coefficient inverse problem in algorithm allowed us to decrease the number of iteration extremely [42]. The theorem of strong convergence of the Landweber iteration method for the coefficient inverse problems for hyperbolic equations was proved in [29].
The gradient of the functional can be computed as follows [11]. Let us introduce the adjoint problem [25,42]: Then the gradient J (q) = (J q 1 (q), J q 2 (q)) has the following form: Thus, in order to make one step of the gradient descent, one has to solve the direct problem (1)-(4), using the current approximation ρ n , c n of the parameters, then solve the adjoint problem (14)- (18), and after that, using the solution of both problems, compute the gradient of the cost functional, using (19) and (20). Since the adjoint problem (14)-(18) also has the form of the hyperbolic system, one can use Godunov-type methods, considered earlier, to compute the solution of the adjoint problem as well.
As we mentioned in the introduction, it is possible to consider the Newton method for solving the inverse problem, which tends to have second-order convergence, as opposed to the first order of gradient-based approaches. For operator formulation of inverse problem (11) the Newton method can be described as follows The drawback of the Newton method is the necessity to solve the additional problem on each iteration, that corresponds to the inversion of the Frechet derivative of the operator, in contrast to gradient methods (13), where one only has to solve the adjoint problems. In our case, the linear two-dimensional inverse problem has the following form: This formulation corresponds to the problem of the reconstruction of the density (yet it changes slightly when considering the reconstruction of the speed of sound). The problem is to find the functionρ(x, y) using the following additional information: Here functions u, v and p are solutions of the direct problem (1)-(4). Thus, we have to deal with the additional inverse problem on every iteration (or, in the case of quasi-Newton methods, every fixed number of iterations) of the scheme. Such an increase in the complexity of each iteration sometimes outweighs the decrease in the total number of iterations in the case of multidimensional problems in Hilbert spaces. Thus, despite the fact the question of efficiency of the Newton-type methods for the problem considered remains interesting, we leave it to future work and focus on the gradient descent in the present paper.

Numerical Results
In this section, we present the results of numerical computations. In this paper, we use synthetic data to study the proposed algorithms. Throughout the computations, we use the following setup, which consists of the test object and the water zone outside the object, where the transducers are located. The radius of the object is 0.07 m, the transducers are uniformly placed in a circle (of the radius 0.115 m) around the object. The acoustic parameters (density, speed of waves propagation) inside the object are chosen as equal to the parameters of the normal human tissue (ρ = 0.9 kg/m 3 , c = 1.2 km/s). We also consider different inclusions inside the object (their quantity and form varies from one model to another), which has different parameters (density varies from ρ = 1.1 kg/m 3 to ρ = 1.3 kg/m 3 in different inclusions , speed of sound varies from c = 1.45 km/s to c = 1.6 km/s). The initial approximation was considered as an object without any inclusions, and we seek to identify said inclusions by solving the inverse problem. Physical We start with the model, presented in Figure 1. It consists of an object and two inclusions, located one inside the other.

Model 1-density distribution
Model 1-speed of sound distribution During the first experiments, we considered the problems of recovering the density and the velocity independently. While one of the medium parameters was unknown, the other one was given exactly (in this case we have the exact values of the second parameter in the whole domain) or inexactly (in this case we do not have the information about the inclusions in the given parameter).
The computational results for this test are presented in Figures 2 and 3. We considered the system of 16 transducers during the computations, and we considered 1000 iterations of the gradient descent. We should mention, that during iteration each position of the source runs through the transducers. When the location of the source is fixed, other transducers work as the receivers. Since that, we work with the array of gradients, each corresponded to the fixed source location, on each iteration. The problem of usage of such data for the better performance of the gradient descent was considered in [14]. The descent parameters were chosen by the trial and error method. One can see, that while the result of recovering one parameter of the system, while the other is given exactly, is accurate enough (graphs a) and (c) on the (Figure 2), the situation changes drastically, when the other parameter is given inexactly (graphs b) and (d) on the ( Figure 2). The errors, introduced in the model by inexact fixed parameter, leads to the fact, that solution, obtained by the gradient descent, tends to drift from the true solution, starting from a certain number of iteration, while the residual functional still decreases. The more precise behavior of the relative errors and residual is presented in Figure 3, where red lines correspond to the exact setup of the fixed parameter, and green lines correspond to the inexact setup of the fixed parameter. Such a situation (when the number of iterations should be chosen according to the errors in the data of inverse problem or in the model itself) arises often, when solving inverse problems. However, such behavior provides additional difficulties, when recovering both parameters simultaneously, as in the following tests.
For the next test, we considered a different model, presented in the figure. In this case, we have one relatively small inclusion in the left part of the object. The model is presented in Figure 4. During the following series of experiments, we considered both the density and speed of sound to be determined. However, because we have only an initial approximation for both functions, the inaccuracy of the model increases, compared to the first series of tests. During the computations, it leads to the accumulation of inaccuracy. In order to overcome such behavior, we proposed the following scheme. We spend some fixed number of iteration, updating only one of the desired parameters (for example, the density). The other parameter (speed of sound) is fixed. After the fixed number of iterations, we switch the fixed-function and the changeable function (start to update the values of the speed of sound, while the density is fixed). Such a scheme allows us to avoid the situation, when the convergence of the descent process, based on minimization of the cost functional, is overwhelmed by the negative effect of the gradient descent with an inaccurate model (demonstrated in the first series of experiments). We considered the several variations of such approach-switch (between updating the density and updating the speed of sound) every 100, 50 and 1 iteration of the gradient descent.  The computational results are presented in Figures 5 and 6. One can see that the simultaneous update of both parameters lead to increasing oscillations in the relative errors ( Figure 5, red curves). The proposed switch scheme allows us to improve the performance of the descent process (some traces of the oscillation behaviour could be seen by considering the green curves). The regular switch between the updates of the density and the speed of sound leads to a more stable and smooth decrease of both the relative errors and the residual. When considering the reconstruction results ( Figure 6), this allows us to obtain the solution with less number of artifacts, while the inclusion is clearly visible. The results of the third series of tests are presented in Figure 7 and the Tables 1 and 2. The exact solution is presented on the upper graphs of Figure 7. It consists of several inclusions inside the objects with different forms and sizes. The diameter of the smallest inclusion is 1cm, which corresponds to the earliest stages of breast cancer genesis. According to the comparative results for the second model, we used the regular switch of the gradient descent between the updates of the density and the speed of sound. The results of the reconstruction for the exact data are presented on the lower graphs of Figure 7. One can see that all three inclusions are visible (despite the fact that the complex structure of the inclusion in the left part of the object is not clearly distinguishable).  Tables 1 and 2 contain the results of the parameter reconstruction for different mesh sizes and the noise in the data. The behaviour of the gradient descent using the noised data was considered, for example, in [14].

Discussion
In this work we considered the coefficient inverse problem of recovering two coefficients (density and the speed of sound) of the first order hyperbolic system, that describes the propagation of the 2D acoustic waves in a heterogeneous medium. We used the secondorder MUSCL-Hancock scheme to solve the direct and adjoint problems, and applied optimization scheme to the coefficient inverse problem.
The presented numerical results illustrate the acceptable accuracy and stability of the proposed algorithm. The ways to improve the method, that we plan to consider in the future work, related to the study of better methods of choosing the descent parameters, consideration of the different ways to measure the data and other related problems. We are also planning to extend the formulation of the problem by considering the attenuation of the waves in the medium.