On the modeling of ultrasound wave propagation in the frame of inverse problem solution

In this paper we consider the inverse problem of detecting the inclusions inside the human tissue by using the acoustic sounding wave. The problem is considered in the form of coefficient inverse problem for first-order system of PDE and we use the gradient descent approach to recover the coefficients of that system. The important part of the sceme is the solution of the direct and adjoint problem on each iteration of the descent. We consider two finite-volume methods of solving the direct problem and study their Influence on the performance of recovering the coefficients.

The system of hyperbolic equations of the first order is considered as the mathematical model of the acoustic tomography, because these equations are obtained directly from the conservation laws of continuum mechanics. This allows us to control the preservation of the basic invariants when solving direct and inverse problems. This is important for solving unstable problems, since the conservation laws of basic invariants are the only criterion for the well-posedness of the solution. Some considerations on the choice of such a mathematical model based on a first-order system and a method for solving such a system can be found in in [32,35,37,38].
We apply the numerical algorithm for solving direct problems, based on the S. K. Godunov schemes [1]. Actually, one of the most popular second-order scheme has been introduced by van Leer [4]; namely the MUSCL scheme. It is a finite volume method where the flux approximation is second-order accurate. This scheme is considered and extended in a lot of applications [6,8,9,13].
MUSCL-Hancock method was introduced in [5] as a variant of the MUSCL scheme. This scheme is full time and space second-order accurate. It differs from the other MUSCL variants [4,10,11,15] in its second-order time of accuracy. Indeed, the time accuracy rises when considering a step based on a central-like scheme [7,16]. The stability of the MUSCL-Hancock scheme was investigated in [20].
Numerical methods of an increased order of accuracy based on the Godunov scheme are widely used and there are a huge number of their variants and implementations [30].
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. The data of inverse problem is the pressure, obtained in the receivers, located in Ω k , k = 1, . . . , N : p(x, y, t) = f k (x, y, t), (x, y) ∈ Ω k , k = 1, . . . , N.
We reformulate inverse problem in the (1)-(4), (6) in operator form Let us reduce the inverse problem (1)-(4), (6) to minimization problem of the following cost functional: We use the gradient-based approach to minimize the functional (8) by considering the following iteration scheme: The gradient of the functional can be obtained by solving the adjoint problem. Explicit formulas are presented in [27,36,37]. At each iteration of gradient descent, the solution of both direct and conjugate problems is required. In [40] it was presented an approach to save twice memory on the stage of adjoint problem and gradient calculation and compare it with usual approach in memory and CPU time cost. Convergence of gradient methods for hyperbolic equations was investigated in [14,18,21]. Work [29] presented a criterion for the termination of the iterative process in the gradient method, consistent with the accumulation of machine round-off errors. What is of interest in this study is how the choice of a numerical method for solving a direct problem affects the solution of the inverse problem.
We consider two suitable methods, that are based on the approach of S.K. Godunov [1]. His approach utilises 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 [4], 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 review of papers related to both methods in [4,20,22].
The short description of the methods is presented below. Considering the direct problem (1)-(4), we use a generalized formulation of the equations (1)-(2) of the following form: 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 (10) 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 [2]. For example, the approximation (10) of the first of the equations (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 others 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 [2], as well as the study of the stability of the scheme. Another approach that we applied for solving the direct problem is MUSCL-Hancock scheme [19,23]. First, we summarize the basic Godunov approach as the combination of two following steps: (i) Obtaining the flux values by solving the Riemann problem; (ii) 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 the each 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, corresponded to 0.5τ . Such update allows the described scheme to second order accuracy [22]. The workflow of the scheme can now be summarized as follows: (i) Reconstruction of the state variables on the current time step, using piecewise-linear extrapolation ; (ii) Evolution of the reconstructed state variables by conservation laws with a time 0.5τ ; (iii) Obtaining the flux values for the "midpoint" time step by solving the Riemann problem for evolved state variables ; (iv) Updating the state values from the current time step on the next time step by conservation laws with a time τ .

The influence of the method of the direct problem solution
In order to compare the usage of two methods, we consider the test model, illustrated by the figure 1. The large circle in the center corresponds to a human tissue, small objects inside are inclusions (values of the density and speed of sound distribution inside the inclusions are higher). As one can see, the model contains several inclusions of different shape and size. The acoustic parameters of human body are taken from [3,17,26].
Density distribution Speed of sound distribution Figure 1. Test -Exact model.
We solve the inverse problems using the following settings. We use the uniform grid with 200 nodes for each space variable. We consider the system of 16 transducers, and use Pusyrev wavelet with frequency ν 0 = 100 KHz. We use the simple gradient descent method with consecutive updates of the density and velocity. During first test we considered the same parameters of descent (α ρ = 250, α c = 500) for both variations of the method. The results are presented on the figure 2.
According to the results, the the error of the solution, when one uses the second order solver of direct problem, tends to decrease faster, but shortly after start the convergence of the iterations fails. Thus, we have to use different descent parameters when changing between the first order and the second order solvers of direct problem. We find the appropriate values of the descent parameters by trial and error approach. During the next test we considered N Iter = 1000 iterations with descent parameters equal to α ρ = 100, α c = 200 when using the classic Godunov's scheme, and α ρ = 25, α c = 50, when using the MUSCL scheme. The comparative results for the residual and relative errors using the mentioned settings are presented on the figure 3. One can notice the the usage of the second order solver provides the more accurate solution. The next figure 4 provides the result of inverse problems solution. The better performance of the descent that uses the second order scheme for the direct problem solution is noticeable. The smallest inclusion was recovered better, as was the complex structure on the leftmost inclusion.    inverse problem for the system of acoustic equations. The better performance of the MUSCL scheme is related to the fact, that the second order method provides more information to the framework of the gradient descent of the functional, when dealing with inclusions of relatively small size and complex structure. In the future work we will study its influence on the quality of the inverse problems solution when using Newton-type methods for solving the inverse problem.