Direct Method for Identification of Two Coefficients of Acoustic Equation

Abstract

We consider the coefficient inverse problem for the 2D acoustic equation. The problem is recovering the speed of sound in the medium (which depends only on the depth) and the density (function of both variables). We describe the method, based on the Gelfand–Levitan–Krein approach, which allows us to obtain both functions by solving two sets of integral equations. The main advantage of the proposed approach is that the method does not use the multiple solution of direct problems, and thus has quite low CPU time requirements. We also consider the variation of the method for the 1D case, where the variation of the wave equation is considered. We illustrate the results with numerical experiments in the 1D and 2D case and study the efficiency and stability of the approach.

1. Introduction

In this paper we consider the coefficient inverse problems for the second-order hyperbolic acoustic equation. Such problems (that can be interpreted as the recovering of the structure of the subsurface by using measurements) obtained on the daylight surface are one of the basic problems in seismology. The applied nature of the coefficient inverse problems for hyperbolic equations explains the wide range of methods, developed for its solution, from kinematic methods to reverse-time migration. While we did not plan to consider the detailed review of existing methods (one can find such survey in [1], for example), we should mention that approaches that utilize the dynamic characteristics of the wave field have tended to become more popular recently. For example, one of the rapidly developing methods is full waveform inversion [2], which is based on the reduction of the problem for the optimization of some misfit functional.

However, methods based on the optimization approach require the solution of the corresponding direct problem on each iteration of the scheme. Despite the recent progress of computational algorithms and hardware, the amount of computations required is still concerning, especially in multi-dimensional cases. The other feature of optimization approaches is their reliance on the basic structure of the medium, which, on the mathematical level, leads to the usage of prior information. In case of poor knowledge of prior information, the efficiency of the methods decreases due to the fact that the corresponding misfit functional usually have several local minima.

In this paper we use the alternate approach, developed originally in works of I.M. Gelfand, B.M. Levitan, M.G. Krein, and V.A. Marchenko [3,4,5]. During the second half of the XX century, the ideas of I.M. Gelfand, B.M. Levitan, and M.G. Krein were developed by A.S. Alekseev [6], A.S. Blagoveschenckiy [7], W. Symes [8], R. Burridge [9], F. Santosa [10], and many other researchers, and these applied for several problems of acoustics, seismics, and geoelectrics [11,12,13,14]. The multi-dimensional variation of the Gelfand–Levitan–Krein (G–L–K) approach was considered in [14,15]. Different numerical algorithms for solving G–L–K equations were considered in [16,17,18]. The G–L–K method was applied for solving coefficient inverse problems of acoustics [19], elasticity [20], and seismics [21].

G–L–K methods allow us to reduce the nonlinear inverse problem to a family of linear integral equations. Such methods that provide the inversion of the forward problem are called direct ones. Another direct method, suitable for solving such problems, is the boundary control method [22,23]. The other important feature of this approach is that one does not have to use the prior information when solving G–L–K equations. There are a few other approaches that do not rely on the usage of a priori information. Aside from boundary control, one can also mention the global convergence method [24,25,26,27,28].

In this work, we use the G–L–K method to recover both the density and the speed of sound from the acoustic equation, where density depends on two variables, while the speed of sound is the function of the depth only. The recovering of the speed of sound as the 2D function by the G–L–K approach was considered in [29], where it was reduced to quite complicated system of equations. In this work, we consider a more simple case of the structure of the speed of sound in order to recover two coefficients of the equation. One should also mention the paper [30], where the wave speed in the acoustic wave equation was estimated from boundary measurements by constructing a reduced-order model matching discrete time–domain data. The success of the proposed algorithm hinges on the data-driven Gram–Schmidt orthogonalization of the snapshots that suppresses multiple reflections and can be viewed as a discrete form of the Marchenko–Gelfand–Levitan–Krein algorithm.

There are some other papers that consider the identification of several parameters, but most of them are dealing with Maxwell’s equations and based on the Carleman estimates [31,32,33,34] and/or the optimization scheme [35]. Similar approaches were used for the system of elasticity in [36] for the problem of recovering two coefficients from the acoustic equation using interior data in [37], and for the system of the first-order acoustic equations in [38].

A method of finding the complex permittivity and permeability of a sample of isotropic material partially filling the rectangular waveguide cross section was presented in [39].

It was proved that the electromagnetic material parameters are uniquely determined by boundary measurements for the time-harmonic Maxwell equations in certain anisotropic settings [40].

In [41], the reconstruction of a complex-valued anisotropic tensor from the knowledge of several internal magnetic fields, which satisfies the anisotropic Maxwell system on a bounded domain with prescribed boundary conditions, was considered.

The recovering of conductivity, permittivity, and the electrokinetic mobility parameter in Maxwell’s equations were analyzed with internal measurements, while allowing the magnetic permeability to be a variable function [42]. It was shown that the knowledge of two internal data sets associated with well-chosen boundary electrical sources uniquely determines these parameters.

In [43], an inverse boundary value problem for the time-harmonic Maxwell equations was considered. The authors showed that the electromagnetic material parameters were determined by boundary measurements where part of the boundary data is measured on a possibly very small set.

It was justified [44] that some elements of the scheme related to the construction of the infinite system of integral equations in the case when the potential is analytic in x. In particular, the convergence of the series in these equations was proved and the conditions for the N-approximation of the system was found.

This paper has the following structure. In the introduction we have considered a brief review of different approaches related to the solution of inverse problems for hyperbolic equations. In Section 2 we formulate the inverse problem for the acoustic equation in 2D and use the G–L–K approach to recover the acoustic impedance. In Section 3 we obtain the additional family of integral equations to recover the auxiliary function, that, with the acoustic impedance, gives us the necessary information to obtain both the density of the medium (which is considered to be a 2D function) and the speed of sound (that depends only on the depth). In Section 4 we consider the 1D problem for the acoustic-type equation, which is connected to the seismic inverse problem. In order to find both the speed of sound and the velocity, we use the property of the solution of the obtained integral equation to recover the unknown functions. In Section 5 we present numerical examples of the solution of the inverse problem. In the discussion we summarize the results of numerical experiments and consider ways for further develop the approach.

2. Two-Dimensional Acoustic Inverse Problem

We consider the following equation that describes the propagation of acoustic waves in the medium:

$$\frac{1}{c^2(z,y)}{u}_{tt}=\Delta u-\nabla ln\rho (z,y)\nabla u.$$

Here $c(z,y)>0$ is the speed of the propagation of the acoustic waves, $\rho (z,y)>0$ is the density of the medium, and $u(z,y,t)$ is the exceeded pressure. As we mentioned in the introduction, we consider that the speed of sound is dependent only on the depth, and the density is dependent on both variables (we should also mention that the method proposed can be naturally extended to the 3D case with respect to the density).

Let us consider the sequence of the direct problems ($y^{\left(0\right)}\in \mathbb{R}$—the set of real numbers):

$$\begin{array}{c}\frac{1}{c^2\left(z\right)}{u}_{tt}(z,y,t)={\Delta }_{z,y}u-{\nabla }_{z,y}ln\rho (z,y){\nabla }_{z,y}u,z>0,t>0,y\in \mathbb{R},\end{array}$$

(1)

$$\begin{array}{c}{u|}_{t<0}\equiv 0,z>0,y\in \mathbb{R},\end{array}$$

(2)

$$\begin{array}{c}{u}_z{|}_{z=0}=\delta \left(t\right)\delta (y-y^{\left(0\right)}),t>0,y\in \mathbb{R}.\end{array}$$

(3)

The boundary condition (3) describes the source of the acoustic waves, which is located on the surface $z=0$ in point $y_0$ and has the form of the Dirac delta function in the time domain.

The inverse problem consists in finding the functions $c\left(z\right)$, $\rho (z,y)$ using the additional information on the surface $z=0$:

$$u(+0,y,t)=f(y,t;y^{\left(0\right)}),t>0,y\in \mathbb{R}.$$

(4)

Now, in order to rewrite Equation (1), we consider the time travel coordinates:

$$x={\int }_0^z\frac{d\xi }{c\left(\xi \right)}$$

Then, the inverse problem (1)–(4) can be rewritten as follows:

$$\begin{array}{c}{u}_{tt}=\sigma (x,y)\left[\frac{\partial }{\partial x}\left(\frac{u_x}{\sigma }\right)+c^2\left(x\right)\frac{\partial }{\partial y}\left(\frac{u_y}{\sigma }\right)\right];\end{array}$$

(5)

$$\begin{array}{c}{u|}_{t<0}\equiv 0,x>0,y\in \mathbb{R}\end{array}$$

(6)

$$\begin{array}{c}{u}_x{|}_{x=0}=c_0\delta \left(t\right)\delta (y-y^{\left(0\right)}),t>0,y\in \mathbb{R}.\end{array}$$

(7)

$$\begin{array}{c}u(+0,y,t)=f(y,t;y^{\left(0\right)}),t>0,y\in \mathbb{R}\end{array}$$

(8)

Here function $\sigma (x,y)=c\left(x\right)\rho (x,y)$ describes the acoustic impedance of the medium. We will also assume that the value $c_0=c(+0)$ is known, and the function $\rho (0,y)$ is also known and sufficiently smooth.

First, we tend to recover the function $\sigma (x,y)$. In order to do that, we apply the Fourier transform with respect to $y^{\left(0\right)}$ to the problem (5)–(8). For simplicity, we consider the case when all the functions considered are $2\pi$ periodical with respect to y. Then, we obtain the following sequence of ratios (here $k\in \mathbb{Z}$):

$$\begin{array}{c}{u}_{tt}^{\left(k\right)}=\sigma (x,y)\left[\frac{\partial }{\partial x}\left(\frac{u_x^{\left(k\right)}}{\sigma }\right)+c^2\left(x\right)\frac{\partial }{\partial y}\left(\frac{u_y^{\left(k\right)}}{\sigma }\right)\right];\end{array}$$

(9)

$$\begin{array}{c}{u}^{\left(k\right)}{|}_{t<0}\equiv 0,x>0,y\in (-\pi ,\pi );\end{array}$$

(10)

$$\begin{array}{c}{u}^{\left(k\right)}{{|}_{y=-\pi }=u^{\left(k\right)}|}_{y=\pi };\end{array}$$

(11)

$$\begin{array}{c}{u}_x^{\left(k\right)}{|}_{x=0}=c_0{e}^{iky}\delta \left(t\right),t>0,y\in \mathbb{R}.\end{array}$$

(12)

$$\begin{array}{c}{u}^{\left(k\right)}(+0,y,t)=f^{\left(k\right)}(y,t),t>0,y\in \mathbb{R}\end{array}$$

(13)

The problem is to recover $\sigma (x,y)$ by the given $f^{\left(k\right)}(y,t),k\in \mathbb{Z}$ set of integers. The problem in (9)–(13) was considered in [14,16], in case of $c\left(x\right)=1$. Using the same approach, based on a combination of the projection method and a 2D analogue of the M.G. Krein approach, we reduce the problem in (9)–(13) to the following set of systems of linear integral equations:

$${\Phi }^{\left(k\right)}(x,t)-\frac{1}{2c_0}\sum _{m\in \mathbb{Z}}\underset{-x}{\overset{x}{\int }}{f_m^{\left(k\right)}}^{\prime }(t-s){\Phi }^{\left(m\right)}(x,s)ds=\frac{1}{2c_0}\underset{-\pi }{\overset{\pi }{\int }}\frac{e^{iky}}{\rho (0,y)}dy.$$

(14)

Here $x<0$, $t\in (-x,x)$, $k\in \mathbb{Z}$, and function $f_m^{\left(k\right)}$ is the Fourier coefficient of the inverse problem’s data with respect to y:

$$\begin{array}{c}\hfill f^{\left(k\right)}(y,t)=\sum _{m\in \mathbb{Z}}{f}_m^{\left(k\right)}\left(t\right)e^{imy},t>0;f^{\left(k\right)}(y,t)=-f^{\left(k\right)}(y,-t),t<0.\end{array}$$

For every given x, the Equation (14) is a system of linear integral equations of the second kind with respect to the unknown functions ${\Phi }^{\left(k\right)}(x,t)$. These functions are connected with $\sigma (x,y)$ by the values on characteristic lines $t=x$:

$${\Phi }^{\left(m\right)}(x,x-0)=\underset{-\pi }{\overset{\pi }{\int }}\frac{e^{imy}}{2\sqrt{\sigma (x,y)\rho (0,y)}}dy.$$

Therefore, when Equation (14) is solved, the acoustic impedance $\sigma (x,y)$ can be calculated as follows:

$$\begin{array}{c}\sigma (x,y)=\frac{{\pi }^2}{\rho (0,y)}{\left[\sum _{m\in \mathbb{Z}}{\Phi }^{\left(m\right)}(x,x-0)e^{-imy}\right]}^{-2}\end{array}$$

(15)

3. Obtaining the Density and the Speed of Sound

Now we suppose that the function $\sigma (x,y)$ is known. However, the given acoustic impedance is not enough to calculate the density or the speed of sound instantly. Therefore, we consider another substitution. Let us consider the function $v(x,y,t)$, which is connected with $u(x,y,t)$ as follows:

$$u(x,y,t)=c_0\sqrt{\frac{\sigma (x,y)}{\sigma (0,y)}}v(x,y,t).$$

Let us also suppose, that the function $\sigma (x,y)$ satisfies the condition ${\sigma }_x{|}_{x=0}=0$. Then, the problem (9)–(12) can be rewritten as follows:

$$\begin{array}{c}{v}_{tt}={\Delta }_{x,y}v(x,y,t)+q(x,y)v;\\ x>0,y\in (-\pi ,\pi ),t>0,\end{array}$$

(16)

$$\begin{array}{c}{v}^{\left(k\right)}{|}_{t<0}\equiv 0,\end{array}$$

(17)

$$\begin{array}{c}{v}_x^{\left(k\right)}{|}_{x=0}=\delta \left(t\right)e^{iky}.\end{array}$$

(18)

The function $q(x,y)$ has the following structure:

$$q(x,y)=\left[\frac{1}{2}{\left(\frac{{\sigma }_x}{\sigma }\right)}_x-\frac{1}{4}{\left(\frac{{\sigma }_x}{\sigma }\right)}^2\right]+c^2\left(x\right)\left[\frac{1}{2}{\left(\frac{{\sigma }_y}{\sigma }\right)}_y-\frac{1}{4}{\left(\frac{{\sigma }_y}{\sigma }\right)}^2\right].$$

(19)

The condition (13) provides the data of the inverse problem:

$$v^{\left(k\right)}(+0,y,t)=\frac{1}{c_0}{f}^{\left(k\right)}(y,t).$$

(20)

The inverse problem (16)–(20) was extensively studied by S.I. Kabanikhin. We reduce it to the 2D analogue of the I.M. Gelfand–B.M. Levitan equation [15,17,18]:

$$\begin{array}{c}{c}_0{w}^{\left(k\right)}(x,y,t)+\underset{-x}{\overset{x}{\int }}\sum _{m\in \mathbb{Z}}{f_m^{\left(k\right)}}^{\prime }(t-s)w^{\left(m\right)}(x,y,s)ds=\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-\frac{1}{2}\left[{f^{\left(k\right)}}^{\prime }(y,t-x)+{f^{\left(k\right)}}^{\prime }(y,t+x)\right],k\in \mathbb{Z}.\end{array}$$

(21)

The connection between the function $q(x,y)$ and the solution of Equation (21) has the following form:

$$q(x,y)=-4\frac{d}{dx}{w}^{\left(0\right)}(x,y,x-0).$$

(22)

Thus, the problem of recovering functions $c\left(x\right)$ and $\rho (x,y)$ for every $x>0$ is reduced to the following procedure:

1.

Solve Equation (14) and use (15) to obtain $\sigma (x,y)$;

2.

Solve Equation (21) and use (22) to obtain $q(x,y)$;

3.

Use the representation (19) to calculate $c\left(x\right)$, when $q(x,y),\sigma (x,y)$ are known;

4.

Recover the density $\rho (x,y)=\frac{\sigma (x,y)}{c\left(x\right)}$.

When $c\left(x\right)$ is calculated, one can reverse the travel-time transform and return to the original coordinates $z={\int }_0^{x}c\left(\xi \right)d\xi$.

However, we should mention that the proposed scheme requires that function

$$q_2(x,y)={\left(\frac{{\sigma }_y}{\sigma }\right)}_y-\frac{1}{2}{\left(\frac{{\sigma }_y}{\sigma }\right)}^2$$

takes a non-zero value in at least one point y. Thus, the horizontal inhomogeneity of the density is necessary. The case when both the coefficients depend only on the depth is considered in the following section, when we have to consider a different form of equation and propose a different approach of the parameter’s reconstruction.

4. One-Dimensional Acoustic Inverse Problem

The 1D formulation of the Equation (1) has the form:

$$\frac{1}{c^2\left(z\right)}\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial z^2}+\frac{\partial ln\left(\rho \right)}{\partial z}\frac{\partial U}{\partial z}.$$

When using the travel-time transform one can obtain:

$$\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial x^2}-\frac{{\sigma }^{\prime }\left(x\right)}{\sigma \left(x\right)}\frac{\partial U}{\partial x}$$

Where, as in the previous section, $\sigma \left(x\right)=c\left(x\right)\rho \left(x\right)$. Thus, as we mentioned previously, in this case one can only recover the acoustic impedance and not the speed and the density independently. Therefore, we alter the formulation of the inverse problem and consider it as follows:

$$\begin{array}{c}\frac{1}{c^2\left(z\right)}\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial z^2}-\frac{\partial ln\rho }{\partial z}\frac{\partial U}{\partial z}-k^{2}U;\end{array}$$

(23)

$$\begin{array}{c}{U(z,t;k)|}_{t<0}\equiv 0;\end{array}$$

(24)

$$\begin{array}{c}\frac{\partial U}{\partial z}{|}_{z=0}=\delta \left(t\right);\end{array}$$

(25)

$$\begin{array}{c}{U(z,t;k)|}_{z=0}=f_k\left(t\right).\end{array}$$

(26)

Here $k\ne 0$ is some given parameter. The inverse problem (23)–(26) is to recover functions $\rho \left(z\right)$, $c\left(z\right)$ by the given function $f\left(t\right)$. Such formulation is connected with the seismic inverse problem in the case of the horizontally layered structure and was considered by A.S. Alekseev [6] in the spectral domain. The connection between the two mentioned problems is based on the special form of the sounding wave and the Hankel transform. Once again, in the travel-time coordinates, we consider the following formulation (for simplicity we also assume that $c\left(0\right)=1$):

$$\begin{array}{c}\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial x^2}-\frac{{\sigma }^{\prime }\left(x\right)}{\sigma \left(x\right)}\frac{\partial U}{\partial x}-k^2{c}^{2}U;\end{array}$$

(27)

$$\begin{array}{c}{U(x,t)|}_{t<0}\equiv 0;\end{array}$$

(28)

$$\begin{array}{c}\frac{\partial U}{\partial x}{|}_{x=0}=\delta \left(t\right);\end{array}$$

(29)

$$\begin{array}{c}{U(x,t)|}_{x=0}=f\left(t\right).\end{array}$$

(30)

The inverse problem (27)–(30) can be considered in the same manner as in the previous section. Yet, we use the 1D case to illustrate the different scheme of recovering two parameters. While the first one is based on the solution of both M.G. Krein and I.M. Gelfand–B.M. Levitan Equations (14) and (21), in this section, we show the possibility of substituting the solution of (21) with some finite-difference ratios.

Once again, we start with the recovering of the acoustic impedance. The 1D analogue of (14) gives us the M.G. Krein equation:

$$-2f_0(+0)V(x,t)-{\int }_{-x}^{x}V(x,s)f_0^{\prime }(t-s)ds=1,t\in (-x,x)$$

(31)

The solution $V(x,t)$ of Equation (31) is connected with the function $\sigma \left(x\right)$ as follows:

$$\sigma \left(x\right)=\frac{V(0,0)}{2V^2(x,x)}.$$

(32)

Thus, in order to recover $\sigma \left(x\right)$, one has to use only the one component of the solution $V(x,t)$ of Equation (31). While this property can be used to increase the efficiency of numerical algorithms [17], now we will assume that the solution of Equation (31) gives us the full set of values of $V(x,t),x>0,t\in (-x,x)$.

Using the framework of [7,14], one can obtain that the function $V(x,t)$ solves the following system:

$$\begin{array}{c}\frac{{\partial }^{2}V}{\partial t^2}=\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\sigma }^{\prime }\left(z\right)}{\sigma \left(z\right)}\frac{\partial V}{\partial x}-k^2{\int }_0^x{c}^2\left(\xi \right)\sigma \left(\xi \right)V_x(\xi ,t)d\xi ,x>0,t\in \mathbb{R};\end{array}$$

(33)

$$\begin{array}{c}{V|}_{x=0}=0;\end{array}$$

(34)

$$\begin{array}{c}{V}_x{|}_{x=0}=\delta \left(t\right).\end{array}$$

(35)

Using the initial conditions (34) and (35), one can obtain that $V(x,t)=0$ for $\left|t\right|>x$, $V(x,-t)=V(x,t),t>0$, and $V(x,x-0)=\frac{1}{2\sqrt{\sigma \left(x\right)\sigma \left(0\right)}}$. Therefore, we replace (33)–(35) with the problem:

$$\begin{array}{c}\frac{{\partial }^{2}V}{\partial t^2}=\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\sigma }^{\prime }\left(z\right)}{\sigma \left(z\right)}\frac{\partial V}{\partial x}-k^2{\int }_0^x{c}^2\left(\xi \right)\sigma \left(\xi \right)V_x(\xi ,t)d\xi ,x>0,t\in (-x,x);\\ {V|}_{t=x}{=V|}_{t=-x}=\frac{1}{2\sqrt{\sigma \left(x\right)\sigma \left(0\right)}}.\end{array}$$

Let us introduce the uniform grid $t_k=kh,x_i=ih$ and consider the finite-difference approximation of (33):

$$\frac{V_i^{k-1}-2V_i^k+V_i^{k+1}}{h^2}=\frac{V_{i-1}^k-2V_i^k+V_{i+1}^k}{h^2}+\frac{V_{i+1}^k-V_{i-1}^k}{2h}\frac{{\sigma }_{i+1}-{\sigma }_i}{h{\sigma }_i}+h{k}^2\sum _{s=1}^i{c}_s^2{P}_s^k;$$

(36)

Here ${\sigma }_i=\sigma \left(x_i\right)$, $c_s^2=c^2\left(x_s\right)$, and

$$P_s^k=\sigma \left(x_s\right)V_x(x_s,t_k)\approx {\sigma }_s\frac{V_s^k-V_{s-1}^k}{h}$$

Using (36) for $t=0$, one can obtain:

$$h{k}^2{c}_i^2{P}_i^0=\frac{1}{h^2}\left(V_i^{-1}+V_i^1-V_{i-1}^0(1-\frac{{\sigma }_{i+1}-{\sigma }_i}{2{\sigma }_i})-V_{i-1}^0(1+\frac{{\sigma }_{i+1}-{\sigma }_i}{2{\sigma }_i})\right)-h{k}^2\sum _{s=1}^{i-1}{c}_s^2{P}_s^0.$$

(37)

The ratio (37) provides the recurrent procedure, which allows us to reconstruct the value of $c^2\left(x\right)$ for each depth by using values ${\sigma }_i$ and $V_i^k$, obtained from (31) and (32).

5. Numerical Results

In this section, we consider the results of the numerical experiments of the reconstruction of velocity and density of the medium in the 1D and 2D cases. We consider the synthetic data obtained by solving the direct problem. For the 1D case, we consider the layered structure of the medium corresponded to the Yurubcheno-Tokhomskoe gas field [45]. The values of the parameters are presented in

Table 1. The 1D model—parameters of the layers.

Depth, km	0.17	0.47	0.87	1.070	1.320	1.600	2.100	2.200	2.300
$v_s\left(z\right)$, km/s	0.9	1.7	3.1	3.5	2.7	3.2	2.85	3.4	2.8
$\rho \left(z\right)$, ${10}^3$ kg/${\mathrm{m}}^3$	2.1	2.4	2.65	2.75	2.5	2.7	2.6	2.75	2.6

.

In order to solve the inverse problem, first we consider the solution of the 1D Krein Equation (31) and calculate the function $V(x,t)$. We solve the integral Equation (31) by discretizing it and reducing it to the SLAE. If $x_M=Mh$, then the system has the following form:

$$(-2f_{0}I-h{\tilde{A}}_{\left[M\right]}){\overline{V}}_{\left[M\right]}={\overline{1}}_{\left[M\right]}.$$

where I is an identity matrix, and

$${\tilde{A}}_{\left[M\right]}=\left[\begin{array}{cccc}{f}_0& f_1& \cdots & f_{2M-1}\\ f_{-1}& f_0& \cdots & f_{2M-2}\\ \ddots & \ddots & \ddots & \ddots \\ f_{-2M+2}& \cdots & f_0& f_1\\ f_{-2M+1}& \cdots & f_{-1}& f_0\end{array}\right].$$

Thus, the structure of Equation (31) allows to reduce it to the system with the Toeplitz matrix. We solve the system by using an adaptation of the Levinson–Durbin algorithm. A full description can be found in [16]. When the Krein equation is solved, we use (37) to recover the speed of sound.

The results of the reconstruction are presented in Figure 1. We used the uniform mesh with $N_t=N_x=200$ during the reconstruction. As one can see, the layered structure is fully recovered and the accuracy is acceptable.

We also reverse the travel-time transform and return to the original depth. The results of the reconstruction of the speed of sound is presented in Figure 2.

We are now moving to the 2D case. We used synthetic data during the numerical experiments, obtained from the solution of the direct problem (9)–(13) by using a combination of the projection method and the finite-difference scheme. We should mention that we used different grid parameters during the solution of the direct and inverse problem (the grid parameters were chosen as $N_{dir}=250$ and $N_{inv}=100$ for the direct and inverse problems, correspondingly). We also introduced random errors into the data, as discussed further.

First we consider the case of smooth parameters. The speed of sound was chosen as a function, closed to linear (the red line in Figure 3), while the density varies with respect to both variables and is presented in Figure 4.

The important parameter of the 2D reconstruction is the number of sources and receivers that describe the quantity of data available during the inverse problems solution. The Figure 3 and Figure 5 illustrate the dependence of the results on the amount of given sources/receivers.

We also present the result of computations for non-smooth parameters, illustrated by Figure 6. However, due to the complex connection between impedance, speed of sound, and the solution of Equation (21) provided by the formula (19), one should consider additional regularizing procedures during the simultaneous solution of Equations (14) and (21). We plan to study the restoration of parameters in the non-smooth case in future work. For now we have focused on the problem of recovering the density of the medium, while we suppose that the speed of sound is known.

During the experiment we considered the model with two inclusions into the homogeneous media (the smaller one is located beneath the larger one). As illustrated by Figure 6, it is possible to locate both inclusions when using more data obtained.

The next experiment was to study the effects of the noise in the data on the solution of the inverse problem. The noise was added to the data according to the formula:

$$f_{err}^{\left(k\right)}(y,t)=f^{\left(k\right)}(y,t)+\epsilon \alpha (y,t)f^{\left(k\right)}(y,t),$$

where $f^{\left(k\right)}(y,t)$ is the value of computed data, $\alpha$ is the random variable, which is uniformly distributed on $(-1,1)$, and $\epsilon$ is the level of noise. The results of computations with noised data are presented in Figure 7 and Figure 8.

The structure of the density of the model, used during experiments with noised data, is presented in Figure 7a. It consists of three non-homogeneous layers and two elliptic inclusions located between the layers. The next picture illustrates the result of the reconstruction for the noiseless data. We should mention that due to the relative complexity of the model, one has to use more data to obtain a relatively accurate solution. Depicted in Figure 7c is the result of adding the noise to the data. One can see the distortion caused by the noise, which becomes more observable with depth. One of the possible ways to decrease the influence of noise is, as presented in the last part of Figure 7, to reduce the number of sources and receivers used during the solution of the equations. In this case, the number of sources can be considered as a regularization parameter.

In the last series of tests we tried to compare the proposed method, based on the G–L–K approach, with the optimization approach. The latter method was broadly described in [46,47] and based on the optimization of the misfit functional by gradient-based methods. The direct and adjoint problem is solved on each iteration by using the finite volume scheme. Since the mentioned method was considered for dealing with problems of acoustic tomography, we chose another model for the last test. The model consists of a round object with several inclusions. Figure 9 provides the comparative analysis of the density reconstruction, while

Table 2. G–L–K method and optimization scheme—time–cost comparison.

Time (s)	G–L–K Method				Optimization
	2 sources	5 sources	8 sources	11 sources	8 sources
$N_x=100$	1.37	5.22	11.87	23.016	300
$N_x=200$	5.17	20.81	52.17	89.25	1700

describes the elapsed CPU time for both methods (the computations were carried out on a laptop). We should also mention that a comparative analysis of the G–L–K method and the optimal control method was carried out in [19], and the comparison with boundary control was considere in [48,49].

6. Discussion

In this paper, we considered the problem of recovering two parameters from the acoustic equation when the speed of sound in the medium depends only on the depth but the density is the function of two variables. As for the 1D case, we considered the variation of the wave equation, which allows for the reconstruction of two parameters. We used the direct method, based on the approach of G–L–K, to solve the mentioned inverse problems and provided several numerical experiments to illustrate the results.

The proposed algorithm demonstrated acceptable accuracy during the synthetic tests. The time cost of the method is also low, because the convolution type of the kernels of the obtained equations allows to use a specific numerical algorithm, based on the Toeplitz matrix inversion. The stability of the method is yet to be improved. One should mention that the issue of stability becomes more complicated due to the fact that we use the derivative of the data for the computation, and because the effects, provided by the noise, tend to accumulate with increasing depth (which is caused by the structure of Equations (14) and (21)). On the other hand, due to the over-determination of the considered inverse problem, one can decrease the impact of the noise by using the quantity of data as the regularization parameter.

We also have to mention that the comparison of the proposed approach and the approach based on the optimization scheme was carried out only on a basic level. We plan to study it in detail in future work. The efficiency of both methods depends on several factors, such as the method of optimization of the functional, the method used for solving the integral Equations (14) and (21), the complexity of the governing equation, etc. However, the main difference of the methods is based on their core structure—the direct nature of the G–L–K approach requires less computation time, while the optimization methods are more versatile and allow us to control the residual on each iteration.

We should mention, that the proposed ideas could be used for solving the inverse problems, in cases when the speed of sound can be considered as a sum of the main trend function, which depends only on the depth, and a small enough horizontal non-homogeneous part. Indeed, let us consider the inverse problem:

$$\begin{array}{c}\frac{1}{c^2(z,y)}{u}_{tt}(z,y,t)=\rho (z,y)\mathrm{div}\left(\frac{\mathrm{grad}u}{\rho }\right),z>0,t>0,y\in \mathbb{R}\end{array}$$

(38)

$$\begin{array}{c}{u|}_{t<0}\equiv 0,x>0,y\in \mathbb{R}\end{array}$$

(39)

$$\begin{array}{c}{u}_z{|}_{z=0}=\delta \left(t\right)\delta (y-y^{\left(0\right)}),t>0,y\in \mathbb{R}.\end{array}$$

(40)

$$\begin{array}{c}u(+0,y,t)=f(y,t;y^{\left(0\right)}),t>0,y\in \mathbb{R}\end{array}$$

(41)

Let us assume that

$$c(z,y)=c_0\left(z\right)+c_1(z,y).$$

In this case, one could linearize the direct problem’s operator and solve (38)–(41) in two steps:

1.

Use the G–L–K approach to recover $c_0\left(z\right)$ by solving (1)–(4);

2.

Recover $c_1(z,y)$ by solving the linearized problem

$$\begin{array}{c}\frac{1}{c_0^2\left(z\right)}{u}_{tt}^{\left(1\right)}(z,y,t)=\rho (z,y)\mathrm{div}\left(\frac{\mathrm{grad},u^{\left(1\right)}}{\rho }\right)+2c_1(z,y)Q(x,y,t),z>0,t>0,y\in \mathbb{R},\end{array}$$

(42)

$$\begin{array}{c}{u|}_{t<0}\equiv 0,x>0,y\in \mathbb{R}\end{array}$$

(43)

$$\begin{array}{c}{u}_z{|}_{z=0}=0,t>0,y\in \mathbb{R}.\end{array}$$

(44)

$$\begin{array}{c}u(+0,y,t)=\left[f-f^{\left(0\right)}\right](y,t;y^{\left(0\right)}),t>0,y\in \mathbb{R}\end{array}$$

(45)

Here the functions $Q(x,y,t)=\frac{1}{c\left(z\right)}\mathrm{div}\left(\frac{\mathrm{grad}{u}^{\left(0\right)}}{\rho }\right)$ and $f^{\left(0\right)}(y,t;y^{\left(0\right)})=u^{\left(0\right)}(+0,y,t)$ depend only on the solution $u^{\left(0\right)}(x,y,t)$ of the direct problem (1)–(3) and can be calculated, since the coefficients of the direct problem were found in the previous step.

Since the last formulated problem is linear, it can be solved by several methods, direct or iterational. In that case that the approach of G–L–K is used in the first step and provides the initial approximation of the parameters, this is improved in the next step of data processing.